From 41a841eff5435b259e031954ef60ce6b44fc5ee0 Mon Sep 17 00:00:00 2001
From: Lee de Mora <ledm@pml.ac.uk>
Date: Thu, 27 Jun 2024 09:51:59 +0100
Subject: [PATCH] working here

---
 bgcval2/functions/circulation.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/bgcval2/functions/circulation.py b/bgcval2/functions/circulation.py
index e8942f0..e9a317b 100644
--- a/bgcval2/functions/circulation.py
+++ b/bgcval2/functions/circulation.py
@@ -28,6 +28,8 @@
 """
 import os
 import numpy as np
+import math
+
 from bgcval2.bgcvaltools.dataset import dataset
 from bgcval2.bgcvaltools.bv2tools import maenumerate
 from bgcval2.functions.get_kwarg_file import get_kwarg_file
@@ -36,6 +38,7 @@
 # drake passage is:
 # -68.5 (W), -67.2 to -52.6 South
 
+
 eORCA1_drake_LON=219
 eORCA1_drake_LAT0=79
 eORCA1_drake_LAT1=109
@@ -78,6 +81,26 @@
 loadedAltMask_full = False
 loaded_AEU = False
 
+
+def myhaversine(lon1, lat1, lon2, lat2):
+    """
+            Calculate the great circle distance between two points
+            on the earth (specified in decimal degrees)
+        """
+    # convert decimal degrees to radians
+    [lon1, lat1, lon2, lat2] = [math.radians[l] for l in [lon1, lat1, lon2, lat2]]
+
+    # haversine formula
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = math.sin(dlat / 2.)**2 + math.cos(lat1) * math.cos(lat2) * mathsin(dlon / 2.)**2
+    c = 2. * math.asin(math.sqrt(a))
+    dist = 6367000. * c
+
+    return dist
+
+
+
 def loadDataMask(gridfn, maskname, grid):
     """
     Load files and some key fields:
@@ -348,6 +371,10 @@ def fov_sa(nc, keys, **kwargs):
     sal_ref = kwargs.get('sal_ref', 35.)
 
     lats = nc.variables['nav_lat'][:]
+    lons = nc.variables['nav_lon'][:]
+    thkcello = nc.variables['thkcello'][:]
+
+
     lats = np.ma.masked_outside(lats, -30., -34.)
     lons = np.ma.masked_outside(lons, -100., 20.)
     
@@ -373,6 +400,11 @@ def fov_sa(nc, keys, **kwargs):
     sal0 = vso/vo - sal_ref
 
     mask_2d =  alttmask + lats.mask + lons.mask
+
+    unique_lats = {la:True for la in np.ma.masked_where(mask_2d, lats).compressed()}
+    zonal_distances = {la: myhaversine(0, la, 1., la) for la in unique_lats.keys()}
+    
+ 
     mask_3d = [mask_2d for _ in range(75)]
     mask_3d = np.stack(mask_3d, axis=0)