diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index f56fba1802c75250268300cd213bbcb0eaa2662c..d322648bde4331ed8275ce2fc7d83d4ff8d130c0 100644
--- a/AirSeaFluxCode.py
+++ b/AirSeaFluxCode.py
@@ -215,9 +215,9 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
np.zeros(spd.shape)*msk)
# cskin parameters
tkt = 0.001*np.ones(T.shape)*msk
- dter = np.ones(T.shape)*0.3*msk
+ dter = -0.3*np.ones(T.shape)*msk
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
- Rnl = 0.97*(5.67e-8*np.power(sst-0.3*cskin, 4)-Rl)
+ Rnl = 0.97*(Rl-5.67e-8*np.power(sst-0.3*cskin, 4))
Qs = 0.945*Rs
dtwl = np.ones(T.shape)*0.3*msk
skt = np.copy(sst)
@@ -240,16 +240,16 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
cq = np.power(kappa, 2)/((np.log(h_in[0]/zo)-psim) *
(np.log(h_in[2]/zoq)-psiq))
monob = -100*np.ones(spd.shape)*msk # Monin-Obukhov length
- tsr = (dt+dter*cskin-dtwl*wl)*kappa/(np.log(h_in[1]/zot) -
+ tsr = (dt-dter*cskin-dtwl*wl)*kappa/(np.log(h_in[1]/zot) -
psit_calc(h_in[1]/monob, meth))
- qsr = (dq+dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
+ qsr = (dq-dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
psit_calc(h_in[2]/monob, meth))
# set-up to feed into iteration loop
it, ind = 0, np.where(spd > 0)
ii, itera = True, -1*np.ones(spd.shape)*msk
tau = 0.05*np.ones(spd.shape)*msk
- sensible = 5*np.ones(spd.shape)*msk
- latent = 65*np.ones(spd.shape)*msk
+ sensible = -5*np.ones(spd.shape)*msk
+ latent = -65*np.ones(spd.shape)*msk
# iteration loop
while np.any(ii):
it += 1
@@ -326,7 +326,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
- skt = np.copy(sst)-dter+dtwl
+ skt = np.copy(sst)+dter+dtwl
elif (skin == "ecmwf"):
dter[ind] = cs_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
@@ -335,7 +335,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
- skt = np.copy(sst)-dter+dtwl
+ skt = np.copy(sst)+dter+dtwl
dqer[ind] = (dter[ind]*0.622*lv[ind]*qsea[ind] /
(287.1*np.power(skt[ind], 2)))
elif (skin == "Beljaars"):
@@ -347,7 +347,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
- skt = np.copy(sst)-dter+dtwl
+ skt = np.copy(sst)+dter+dtwl
dqer = dter*0.622*lv*qsea/(287.1*np.power(skt, 2))
else:
dter[ind] = np.zeros(sst[ind].shape)
@@ -362,8 +362,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
np.round(np.nanmedian(usr), 3),
np.round(np.nanmedian(tsr), 4),
np.round(np.nanmedian(qsr), 7))
- Rnl[ind] = 0.97*(5.67e-8*np.power(sst[ind] -
- dter[ind]*cskin, 4)-Rl[ind])
+ Rnl[ind] = 0.97*(Rl[ind]-5.67e-8*np.power(sst[ind] +
+ dter[ind]*cskin, 4))
t10n[ind] = (Ta[ind] -
tsr[ind]/kappa*(np.log(h_in[1, ind]/ref_ht)-psit[ind]))
q10n[ind] = (qair[ind] -
@@ -371,11 +371,10 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
tv10n[ind] = t10n[ind]*(1+0.6077*q10n[ind])
tsrv[ind], monob[ind], Rb[ind] = get_L(L, lat[ind], usr[ind], tsr[ind],
qsr[ind], h_in[:, ind], Ta[ind],
- sst[ind]-dter[ind]*cskin+dtwl[ind]*wl,
+ sst[ind]+dter[ind]*cskin+dtwl[ind]*wl,
qair[ind], qsea[ind], wind[ind],
np.copy(monob[ind]), zo[ind],
zot[ind], psim[ind], meth)
- # sst[ind]-dter[ind]*cskin+dtwl[ind]*wl
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
@@ -403,8 +402,8 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
utmp = np.copy(u10n)
utmp = np.where(utmp < 0, np.nan, utmp)
itera[ind] = np.ones(1)*it
- sensible = -rho*cp*usr*tsr
- latent = -rho*lv*usr*qsr
+ sensible = rho*cp*usr*tsr
+ latent = rho*lv*usr*qsr
if (gust[0] == 1):
tau = rho*np.power(usr, 2)*(spd/wind)
elif (gust[0] == 0):
@@ -521,7 +520,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
flag+[","]+["o"], flag))
- res = np.zeros((39, len(spd)))
+ res = np.zeros((40, len(spd)))
res[0][:] = tau
res[1][:] = sensible
res[2][:] = latent
@@ -561,14 +560,15 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
res[36][:] = np.sqrt(np.power(wind, 2)-np.power(spd, 2))
res[37][:] = Rb
res[38][:] = rh
+ res[39][:] = tkt
if (out == 0):
res[:, ind] = np.nan
# set missing values where data have non acceptable values
res = np.asarray([np.where(q10n < 0, np.nan,
- res[i][:]) for i in range(39)])
+ res[i][:]) for i in range(40)])
res = np.asarray([np.where(u10n < 0, np.nan,
- res[i][:]) for i in range(39)])
+ res[i][:]) for i in range(40)])
# output with pandas
resAll = pd.DataFrame(data=res.T, index=range(len(spd)),
columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct",
@@ -577,7 +577,7 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
"t10n", "tv10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "iteration", "dter",
"dqer", "dtwl", "qair", "qsea", "Rl", "Rs",
- "Rnl", "ug", "Rib", "rh"])
+ "Rnl", "ug", "Rib", "rh", "delta"])
resAll["flag"] = flag
return resAll
diff --git a/flux_subs.py b/flux_subs.py
index 93af1e5b7344da1cf4673306d890828585a47b76..14182fdd965f8810d1fab702ee4aa43f8f997726 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -563,7 +563,7 @@ def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
Rs : float
downward shortwave radiation [Wm-2]
Rnl : float
- upwelling IR radiation [Wm-2]
+ net upwelling IR radiation [Wm-2]
cp : float
specific heat of air at constant pressure [J/K/kg]
lv : float
@@ -601,12 +601,12 @@ def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
bigc = 16*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(rho, 2))
wetc = 0.622*lv*qsea/(287.1*np.power(sst+273.16, 2))
Rns = 0.945*Rs # albedo correction
- shf = -rho*cp*usr*tsr
- lhf = -rho*lv*usr*qsr
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl
fs = 0.065+11*delta-6.6e-5/delta*(1-np.exp(-delta/8.0e-4))
- qcol = Qnsol-Rns*fs
- alq = aw*qcol+0.026*lhf*cpw/lv
+ qcol = Qnsol+Rns*fs
+ alq = aw*qcol+0.026*np.minimum(lhf, 0)*cpw/lv
xlamx = 6*np.ones(sst.shape)
xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
@@ -617,7 +617,7 @@ def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
# ---------------------------------------------------------------------
-def delta(aw, Qnsol, usr, lat):
+def delta(aw, Q, usr, lat):
"""
Computes the thickness (m) of the viscous skin layer.
Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
@@ -627,7 +627,7 @@ def delta(aw, Qnsol, usr, lat):
----------
aw : float
thermal expansion coefficient of sea-water [1/K]
- Qnsol : float
+ Q : float
part of the net heat flux actually absorbed in the warm layer [W/m^2]
usr : float
friction velocity in the air (u*) [m/s]
@@ -644,11 +644,9 @@ def delta(aw, Qnsol, usr, lat):
# u* in the water
usr_w = np.maximum(usr, 1e-4)*np.sqrt(1.2/rhow) # rhoa=1.2
rcst_cs = 16*gc(lat)*np.power(visw, 3)/np.power(tcw, 2)
- lm = 6/np.power(1+np.power(np.maximum(Qnsol*aw*rcst_cs /
- np.power(usr_w, 4), 0), 3/4),
- 1/3)
+ lm = 6*(1+np.maximum(Q*aw*rcst_cs/np.power(usr_w, 4), 0)**0.75)**(-1/3)
ztmp = visw/usr_w
- delta = np.where(Qnsol > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
+ delta = np.where(Q > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
return delta
# ---------------------------------------------------------------------
@@ -689,16 +687,16 @@ def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, lat):
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
- Rns = 0.945*Rs # (1-0.066)*Rs net solar radiation (albedo correction)
- shf = -rho*cp*usr*tsr
- lhf = -rho*lv*usr*qsr
+ Rns = 0.945*Rs # net solar radiation (albedo correction)
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl # eq. 8.152
d = delta(aw, Qnsol, usr, lat)
for jc in range(4): # because implicit in terms of delta...
# # fraction of the solar radiation absorbed in layer delta eq. 8.153
# and Eq.(5) Zeng & Beljaars, 2005
fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8e-4))
- Q = Qnsol-fs*Rns
+ Q = Qnsol+fs*Rns
d = delta(aw, Q, usr, lat)
dtc = Q*d/0.6 # (rhow*cw*kw)eq. 8.151
return dtc
@@ -753,17 +751,17 @@ def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, lat):
a1, a2, a3 = 0.28, 0.27, 0.45
b1, b2, b3 = -71.5, -2.8, -0.06 # [m-1]
Rd = 1-(a1*np.exp(b1*rd0)+a2*np.exp(b2*rd0)+a3*np.exp(b3*rd0))
- shf = -rho*cp*usr*tsr
- lhf = -rho*lv*usr*qsr
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl
usrw = np.maximum(usr, 1e-4 )*np.sqrt(1.2/rhow) # u* in the water
zc3 = rd0*kappa*gc(lat)/np.power(1.2/rhow, 3/2)
zc4 = (0.3+1)*kappa/rd0
zc5 = (0.3+1)/(0.3*rd0)
for jwl in range(10): # itteration to solve implicitely equation for warm layer
- dsst = skt-sst+dtc
+ dsst = skt-sst-dtc
## Buoyancy flux and stability parameter (zdl = -z/L) in water
- ZSRD = (Qnsol-Rns*Rd)/(rhow*cpw)
+ ZSRD = (Qnsol+Rns*Rd)/(rhow*cpw)
ztmp = np.maximum(dsst, 0)
zdl = np.where(ZSRD > 0, 2*(np.power(usrw, 2) *
np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))+ZSRD,
@@ -825,8 +823,8 @@ def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, lat, Qs):
cpw = 4190 # specific heat capacity of water
aw = 3e-4 # thermal expansion coefficient [K-1]
Rns = 0.945*Rs # net solar radiation (albedo correction)
- shf = -rho*cp*usr*tsr
- lhf = -rho*lv*usr*qsr
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
Q = Rnl+shf+lhf+Qs
xt = (16*Q*g*aw*cpw*np.power(rhow*visw, 3))/(np.power(usr, 4) *
np.power(rho*tcw, 2))