79 new varmod poc

src/bitsea/instruments/superfloat.py

@@ -7,6 +7,7 @@
 from bitsea.instruments.instrument import Instrument, Profile
 from scipy.optimize import curve_fit
 from bitsea.instruments.var_conversions import FLOATVARS as conversion
+from bitsea.validation.deliverables.metrics import MLD
 
 mydtype= np.dtype([
           ('file_name','S200'),
@@ -138,15 +139,39 @@ def read(self,var,var_mod=None):
                 ii=Profile<=0
                 Profile[ii] = 0.0
 
+        #if (var_mod == "POC"):
+        #    POC = Profile *  52779.37 - 3.57 # Bellacicco 2019
+        #    if ii.sum() > 0 :
+        #        shift=POC[ii].mean()
+        #        print( "POC: adding a shift of " + str(shift))
+        #        Profile = POC - shift
+        #        ii=Profile<=0
+        #        Profile[ii] = 0.0
         if (var_mod == "POC"):
-            POC = Profile *  52779.37 - 3.57 # Bellacicco 2019
-            if ii.sum() > 0 :
-                shift=POC[ii].mean()
-                print( "POC: adding a shift of " + str(shift))
-                Profile = POC - shift
-                ii=Profile<=0
-                Profile[ii] = 0.0
-
+            # POC/bbp700 ratio after Gali' et al. 2022
+            Pres_PSAL, PSAL, Qc_PSAL = self.read_raw('PSAL')
+            Pres_TEMP, TEMP, Qc_TEMP = self.read_raw('TEMP')
+            z_mld = MLD(TEMP, PSAL, Pres_TEMP, insitu_T=True)[0]
+            # Interpolate TEMP @ Z = Pres, because Pres and Pres_TEMP are different
+            TEMP = np.interp(Pres, Pres_TEMP, TEMP) #piecewise linear interpolation
+            PSAL = np.interp(Pres, Pres_PSAL, PSAL) #piecewise linear interpolation
+            # Gali et al. 2022 parameters
+            #c = 1000.0      # [mmol_C m^-2] best value globally
+            c = 500.0      # [mmol_C m^-2] best fit in subtropical biomes
+            c = c * 12.011  # [mg_C m^-2]
+            b = 6.57        # [m^-1]
+            a = 41305.0     # [mg_C m^-2]
+            z_surf_Med = 14.0 # [m]
+            #
+            Z_term = np.maximum(0.0, Pres - z_surf_Med)
+            PB_ratio = c + a * np.exp(-0.001 * b * Z_term)
+            #
+            z_ref = np.max([z_surf_Med, z_mld])
+            z_term_i = np.max([0.0, z_ref - z_surf_Med])
+            PB_zref = c + a * np.exp(-0.001 * b * z_term_i)
+            PB_ratio[Pres<=z_ref] = PB_zref
+            # POC / bbp700 ratio
+            Profile = Profile * PB_ratio
 
         return Pres, Profile, Qc
 

src/bitsea/validation/deliverables/metrics.py

@@ -30,22 +30,19 @@ def find_DCM(Chl_profile, zlev):
     return CM, DCM
 
 
-def MLD(Temperature, Salinity, Pres):
+def MLD(Temperature, Salinity, Pres, insitu_T=False):
     """Calculation of Mixed Layer Depth based on temperature difference of 0.2
     mld is defined as the level where the difference of temperature with respect the reference level of 10m
     is of 0.2C
     It resurns also DENSITY (SIGMA) and POTENTIAL DENSITY (SIGMA THETA)
     """
     th = 10  # threshold of depth minimum
-    i_less_than_10 = Pres < th
-    i_10 = (
-        i_less_than_10[-1] + 1
-    )  # (this is the index corresponding approx ~10m of the zlevel array)
-    i_10 = (
-        i_less_than_10[-1] + 1
-    )  # (this is the index corresponding approx ~10m of the zlevel array)
+    i_10 = np.abs(Pres - th).argmin()
     T = Temperature
     S = Salinity
+    # if T is in-situ (e.g. float data) convert to potential temperature 
+    if insitu_T: 
+        T = sw.ptmp(S, T, Pres)
     D1000 = Pres[
         (Pres < 1000) & (Pres > th)
     ]  # CONSIDER VALUES the reference level at 10m (de Boyer-Montegut 2004)