diff --git a/src/bitsea/validation/deliverables/metrics.py b/src/bitsea/validation/deliverables/metrics.py index 7b7ed0c4..d80d419e 100644 --- a/src/bitsea/validation/deliverables/metrics.py +++ b/src/bitsea/validation/deliverables/metrics.py @@ -30,32 +30,30 @@ def find_DCM(Chl_profile, zlev): return CM, DCM -def MLD(Temperature, Salinity, Pres): - """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) +def MLD(Temperature, Salinity, Pres, insitu_T=True): + """ + Calculation of Mixed Layer Depth based on temperature difference of 0.2 + inputs are Temperature ([C], in-situ, optionally conservative temperature), Salinity ([psu]), Pressure ([kg/m^3]) + 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 and POTENTIAL DENSITY at mld """ th = 10 # threshold of depth minimum i_10 = np.abs(Pres - th).argmin() - T = Temperature - S = Salinity - D1000 = Pres[ - (Pres < 1000) & (Pres > th) - ] # CONSIDER VALUES the reference level at 10m (de Boyer-Montegut 2004) - T1000 = T[(Pres < 1000) & (Pres > th)] - S1000 = S[(Pres < 1000) & (Pres > th)] - SA1000 = gsw.SA_from_SP(S1000, D1000, 0, 0) - CT1000 = gsw.CT_from_t(SA1000, T1000, D1000) - Dens1000 = gsw.rho(SA1000, CT1000, D1000) - 1000 # DENSITY # SIGMA - PDens1000 = ( - gsw.sigma0(SA1000, CT1000) - ) # POTENTIAL DENSITY # SIGMA THETA - for ip, p in enumerate(T1000): - # abs_diff=abs(p-T1000[0]) - # if abs_diff > 0.2: - # abs_diff=abs(p-np.mean(T[Pres th)] + CT1000 = CT[(Pres < 1000) & (Pres > th)] + SA1000 = SA[(Pres < 1000) & (Pres > th)] + # in situ density anomaly [kg/m^3] + Dens1000 = gsw.rho(SA1000, CT1000, D1000) - 1000 # in situ density anomaly [kg/m^3] DENSITY # SIGMA + # potential density anomaly with reference pressure of 0 dbar + PDens1000 = (gsw.sigma0(SA1000, CT1000)) + for ip, p in enumerate(CT1000): + abs_diff = abs(p - CT[i_10]) if abs_diff > 0.2: break MixedLayerDepth = D1000[ip]