Skip to content

GitLab

Commit 0deb28e4 authored by sbiri's avatar sbiri
Browse files

changed sign convention for heat fluxes 

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