From 9bda0f47e6146995bf293f031e48e334fdfe8358 Mon Sep 17 00:00:00 2001 From: sbiri <sbiri@noc.ac.uk> Date: Mon, 29 Mar 2021 13:06:25 +0100 Subject: [PATCH] - removed restriction for rh > 100 to lead to NaN - set default method S88 --- AirSeaFluxCode.py | 7 +++---- hum_subs.py | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py index e5d77a2..297f444 100644 --- a/AirSeaFluxCode.py +++ b/AirSeaFluxCode.py @@ -12,7 +12,7 @@ from flux_subs import (cs_C35, cs_Beljaars, cs_ecmwf, wl_ecmwf, def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, Rl=None, Rs=None, cskin=None, skin="C35", wl=0, gust=None, - meth="S80", qmeth="Buck2", tol=None, n=10, out=0, L=None): + meth="S88", qmeth="Buck2", tol=None, n=10, out=0, L=None): """ Calculates turbulent surface fluxes using different parameterizations Calculates height adjusted values for spd, T, q @@ -379,7 +379,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, wind[ind] = np.copy(spd[ind]) u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) - psim[ind]) - # usr[ind]/np.sqrt(cd10n[ind]) if (it < 4): # make sure you allow small negative values convergence u10n = np.where(u10n < 0, 0.5, u10n) flag = np.where((u10n < 0) & (flag == "n"), "u", @@ -490,7 +489,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, rh = np.ones(sst.shape)*80 elif (hum[0] == 'rh'): rh = hum[1] - rh = np.where(rh > 100, np.nan, rh) + # rh = np.where(rh > 100, np.nan, rh) elif (hum[0] == 'Td'): Td = hum[1] # dew point temperature (K) Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td)) @@ -498,7 +497,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, esd = 611.21*np.exp(17.502*((Td-CtoK)/(Td-32.19))) es = 611.21*np.exp(17.502*((T-CtoK)/(T-32.19))) rh = 100*esd/es - rh = np.where(rh > 100, np.nan, rh) + # rh = np.where(rh > 100, np.nan, rh) res = np.zeros((39, len(spd))) res[0][:] = tau diff --git a/hum_subs.py b/hum_subs.py index 98591c1..cba06af 100644 --- a/hum_subs.py +++ b/hum_subs.py @@ -380,7 +380,7 @@ def get_hum(hum, T, sst, P, qmeth): if (np.all(RH < 1)): sys.exit("input relative humidity units should be \%") qair, qsea = np.nan, np.nan - RH = np.where(RH > 100, np.nan, RH) # ensure RH <=100 + # RH = np.where(RH > 100, np.nan, RH) # ensure RH <=100 qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (kg/kg) qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (kg/kg) elif (hum[0] == 'q'): @@ -393,7 +393,7 @@ def get_hum(hum, T, sst, P, qmeth): esd = 611.21*np.exp(17.502*((Td-273.16)/(Td-32.19))) es = 611.21*np.exp(17.502*((T-273.16)/(T-32.19))) RH = 100*esd/es - RH = np.where(RH > 100, np.nan, RH) # ensure RH <=100 + # RH = np.where(RH > 100, np.nan, RH) # ensure RH <=100 qair = qsat_air(T, P, RH, qmeth)/1000 # q of air (kg/kg) qsea = qsat_sea(sst, P, qmeth)/1000 # surface water q (kg/kg) return qair, qsea -- GitLab