diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index 124859864f721567f6bbbd1ceac20c28510cd3d8..154075f983de7d10fc165846eac00c79fab8c380 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -34,13 +34,13 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
hout : float
output height, default is 10m
zi : int
- PBL height (m) called in C35
+ PBL height (m)
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
jcool : int
- 0 if sst is true ocean skin temperature called in COARE, else 1
+ 0 if sst is true ocean skin temperature, else 1
n : int
number of iterations
@@ -89,6 +89,7 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
g = gc(lat, None) # acceleration due to gravity
ctn, ct, cqn, cq = (np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan,
np.zeros(spd.shape)*np.nan, np.zeros(spd.shape)*np.nan)
+ cdn = np.zeros(spd.shape)*np.nan
# if input values are nan break
if (np.all(np.isnan(spd)) or np.all(np.isnan(T)) or np.all(np.isnan(SST))):
sys.exit("input wind, T or SST is empty")
@@ -119,16 +120,13 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
zi = 600 # set to default for COARE3.5
elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
zi = 1000
- if (np.all(np.isnan(lat)) and (meth == "C30" or meth == "C35" or
- meth == "C40")):
- lat=45*np.ones(np.shape(spd))
####
- th = np.where(T < 200, (np.copy(T)+CtoK) *
+ th = np.where(np.nanmax(T) < 200, (np.copy(T)+CtoK) *
np.power(1000/P,287.1/1004.67),
np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T
- Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*hh_in[1],
+ Ta = np.where(np.nanmax(T) < 200, np.copy(T)+CtoK+tlapse*hh_in[1],
np.copy(T)+tlapse*hh_in[1]) # convert to Kelvin if needed
- sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
+ sst = np.where(np.nanmax(SST) < 200, np.copy(SST)+CtoK, np.copy(SST))
if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
meth == "ERA5"):
qsea = qsat26sea(sst, P)/1000 # surface water specific humidity (g/kg)
@@ -161,10 +159,8 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
dtv=dt*(1.+0.61*qair)+0.61*th*dq
# ------------
rho = P*100/(287.1*tv10n)
-# rho=P*100/(287.1*sst*(1+0.61*qsea)) # Zeng et al. 1998
lv = (2.501-0.00237*SST)*1e6
cp = 1004.67*(1 + 0.00084*qsea)
-# cp = 1004.67*np.ones(Ta.shape) # Zeng et al. 1998, C3.0, C3.5
u10n = np.copy(spd)
monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
cdn = cdn_calc(u10n, Ta, None, lat, meth)
@@ -176,7 +172,7 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
# if jcool == 1:
tkt = 0.001*np.ones(T.shape)
Rnl = 0.97*(5.67e-8*np.power(sst-0.3*jcool+CtoK, 4)-Rl)
- dter = np.ones(T.shape)*0.3
+ dter = np.ones(T.shape)*0.3
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
if (meth == "UA"):
wind = np.where(dtv >= 0, np.where(spd > 0.1, spd, 0.1),
@@ -262,9 +258,9 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
tsr = (dt+dter*jcool)*kappa/(np.log(hin[1]/zot) -
psit_calc(hh_in[1]/monob, meth))
qsr = (dq+dqer*jcool)*kappa/(np.log(hin[2]/zot) -
- psit_calc(hh_in[2]/monob, meth))
+ psit_calc(hh_in[2]/monob, meth))
# tolerance for u,t,q,usr,tsr,qsr
- tol = np.array([1e-06, 0.01, 5e-07, 1e-06, 0.001, 5e-07])
+ tol = np.array([0.01, 0.01, 5e-05, 0.005, 0.001, 5e-07])
it, ind = 0, np.where(spd > 0)
ii, itera = True, np.zeros(spd.shape)*np.nan
while np.any(ii):
@@ -369,17 +365,21 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa /
(np.log(hin[1]/zot[ind])-psit_26(hin[1]/monob[ind]))))
else:
- usr[ind] = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
- psim_calc(hh_in[0]/monob[ind], meth)))#np.sqrt(cd[ind]*wind[ind]**2)
+ usr[ind] = np.sqrt(cd[ind]*np.power(wind[ind], 2))
tsr[ind] = ct[ind]*wind[ind]*(dt[ind]+dter[ind]*jcool)/usr[ind]
qsr[ind] = cq[ind]*wind[ind]*(dq[ind]+dqer[ind]*jcool)/usr[ind]
- dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
- rho[ind], Rl[ind],
- Rs[ind], Rnl[ind],
- cp[ind], lv[ind],
- np.copy(tkt[ind]),
- usr[ind], tsr[ind],
- qsr[ind], lat[ind])
+ if (jcool == 1):
+ dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
+ rho[ind], Rl[ind],
+ Rs[ind], Rnl[ind],
+ cp[ind], lv[ind],
+ np.copy(tkt[ind]),
+ usr[ind], tsr[ind],
+ qsr[ind], lat[ind])
+ else:
+ dter[ind] = np.zeros(spd[ind].shape)
+ dqer[ind] = np.zeros(spd[ind].shape)
+ tkt[ind] = np.zeros(spd[ind].shape)
Rnl[ind] = 0.97*(5.67e-8*np.power(SST[ind] -
dter[ind]*jcool+CtoK, 4)-Rl[ind])
t10n[ind] = (Ta[ind] -
@@ -451,8 +451,8 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
np.copy(usr), np.copy(tsr), np.copy(qsr)])
d = np.abs(new-old)
ind = np.where((d[0, :] > tol[0])+(d[1, :] > tol[1]) +
- (d[2, :] > tol[2])+(d[3, :] > tol[3]) +
- (d[4, :] > tol[4])+(d[5, :] > tol[5]))
+ (d[2, :] > tol[2])+(d[3, :] > tol[3]) +
+ (d[4, :] > tol[4])+(d[5, :] > tol[5]))
if (ind[0].size == 0):
ii = False
else:
@@ -462,7 +462,6 @@ def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, zi=600,
ind[0].size)
# calculate output parameters
rho = (0.34838*P)/(tv10n)
-# rho = P*100./(287.1*(T+CtoK)*(1+0.61*qair)) # C35
t10n = t10n-(273.16+tlapse*ref_ht)
sensible = -rho*cp*usr*tsr
latent = -rho*lv*usr*qsr