Skip to content

GitLab

Commit 22be86dd authored by sbiri's avatar sbiri
Browse files

UA and ERA5 algorithms added

parent 9b617551
Hide whitespace changes
Inline Side-by-side
Showing with 985 additions and 350 deletions
flux_calc_v3.py
View file @ 22be86dd
import numpy as np
import sys
import logging
from flux_subs import get_heights,cdn_calc, cd_calc, get_skin, psim_calc, \
psit_calc, ctcq_calc,ctcqn_calc, get_gust, gc,q_calc,qsea_calc,qsat26sea, qsat26air, \
visc_air, psit_26,psiu_26, psiu_40,cd_C35
from flux_subs import (kappa, CtoK, get_heights, cdn_calc, cd_calc, get_skin,
psim_calc, psit_calc, ctcq_calc, ctcqn_calc, get_gust,
gc, q_calc, qsea_calc, qsat26sea, qsat26air,
visc_air, psit_26, psiu_26, psiu_40, cd_C35,
charnock_C35)
def flux_calc_v3(spd, T, SST, lat, RH, P, hin, hout, wcp, sigH, zi = 600, \
Rl=None, Rs=None, jcool=1, method="Smith80",n=20):
"""
inputs:
flux method ("Smith80","Smith88","LP82","HEXOS","HEXOSwave","YT96","COARE3.0","COARE3.5","LY04")
spd : relative wind speed in m/s (is assumed as magnitude difference between wind
and surface current vectors)
T : air temperature in K (will convert if < 200)
SST : sea surface temperature in K (will convert if < 200)
lat : latitude
rh : relative humidity in %, default 80% for C35
P : air pressure, default 1013 for C35
hin : sensor heights in m (array of 1->3 values: 1 -> u=t=q; 2 -> u,t=q; 3 -> u,t,q ) default 10m
hout: default heights are 10m
wcp : phase speed of dominant waves (m/s) called in COARE3.5
sigH: significant wave height (m) called in COARE3.5
zi : PBL height (m) called in COARE3.5
Rl : downward longwave radiation (W/m^2)
Rs : downward shortwave radiation (W/m^2)
jcool: 0 if sst is true ocean skin temperature called in COARE3.5
n : number of iterations, for COARE3.5 set to 10
outputs:
res : which contains 1. tau 2. sensible heat 3. latent heat 4. Monin-Obhukov length (monob)
5. drag coefficient (cd) 6. neutral drag coefficient (cdn) 7. ct 8. ctn 9. cq 10. cqn
11.tsrv 12. tsr 13. qsr 14. usr 15. psim 16. psit 17. u10n 18. t10n 19. tv10n
20. q10n 21. zo 22. zot 23. zoq 24. urefs 25.trefs 26. qrefs 27. iterations
ind : the indices in the matrix for the points that did not converge after the maximum number of iterations
based on bform.f and flux_calc.R modified and ported to python by S. Biri
def flux_calc_v3_1(spd, T, SST, lat, RH, P, hin, hout, wcp, sigH, zi=600,
Rl=None, Rs=None, jcool=1, meth="S88", n=10):
""" Calculates momentum and heat fluxes using different parameterizations
Parameters
----------
meth : str
"S80","S88","LP82","YT96","UA","HEXOS","HEXOSwave",
"LY04","C30","C35","ERA5"
spd : float
relative wind speed in m/s (is assumed as magnitude difference
between wind and surface current vectors)
T : float
air temperature in K (will convert if < 200)
SST : float
sea surface temperature in K (will convert if < 200)
lat : float
latitude (deg)
RH : float
relative humidity in %
P : float
air pressure
hin : float
sensor heights in m (array of 1->3 values: 1 -> u=t=q; /
2 -> u,t=q; 3 -> u,t,q ) default 10m
hout : float
output height, default is 10m
wcp : float
phase speed of dominant waves (m/s) called in C35
sigH : float
significant wave height (m) called in C35
zi : int
PBL height (m) called in C35
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
jcool : bool
0 if sst is true ocean skin temperature called in C35
n : int
number of iterations, for C35 set to 10
out : str
format in which to save output ('output.nc' or None)
Returns
-------
res : array that contains
1. momentum flux (W/m^2)
2. sensible heat (W/m^2)
3. latent heat (W/m^2)
4. Monin-Obhukov length (mb)
5. drag coefficient (cd)
6. neutral drag coefficient (cdn)
7. heat exhange coefficient (ct)
8. neutral heat exhange coefficient (ctn)
9. moisture exhange coefficient (cq)
10. neutral moisture exhange coefficient (cqn)
11. star virtual temperature (tsrv)
12. star temperature (tsr)
13. star humidity (qsr)
14. star velocity (usr)
15. momentum stability function (psim)
16. heat stability funciton (psit)
17. 10m neutral velocity (u10n)
18. 10m neutral temperature (t10n)
19. 10m neutral virtual temperature (tv10n)
20. 10m neutral specific humidity (q10n)
21. surface roughness length (zo)
22. heat roughness length (zot)
23. moisture roughness length (zoq)
24. velocity at reference height (urefs)
25. temperature at reference height (trefs)
26. specific humidity at reference height (qrefs)
27. number of iterations until convergence
ind : int
the indices in the matrix for the points that did not converge
after the maximum number of iterations
The code is based on bform.f and flux_calc.R modified by S. Biri
"""
logging.basicConfig(filename='flux_calc.log', level=logging.INFO)
kappa,CtoK = 0.4, 273.16 # von Karman's constant; conversion factor from degC to K
logging.basicConfig(filename='flux_calc.log',
format='%(asctime)s %(message)s',level=logging.INFO)
ref_ht, tlapse = 10, 0.0098 # reference height, lapse rate
hh_in = get_heights(hin) # heights of input measurements/fields
hh_out = get_heights(hout) # desired height of output variables,default is 10
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
hh_out = get_heights(hout) # desired height of output variables
if np.all(np.isnan(lat)): # set latitude to 45deg if empty
lat=45*np.ones(np.shape(spd))
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)
# 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")
logging.debug('all input is nan')
if (np.all(np.isnan(RH)) or np.all(np.isnan(P))):
if (method=="COARE3.5"):
RH=np.ones(np.shape(spd))*80
P=np.ones(np.shape(spd))*1013
else:
sys.exit("input RH or P is empty")
if (np.all(spd[~np.isnan(spd)]== 0) and method != "COARE3.0"):
sys.exit("wind cannot be zero when method other than COARE3.0")
logging.debug('all velocity input is zero')
if (np.all(np.isnan(Rl)) and method=="COARE3.5"):
Rl=np.ones(np.shape(spd))*370 # set to default
if (np.all(np.isnan(Rs)) and method=="COARE3.5"):
Rs=np.ones(np.shape(spd))*150 # set to default
if (np.all(np.isnan(zi)) and method=="COARE3.5"):
zi=600 # set to default
### additional parameters for COARE3.5
waveage,seastate=1,1
if (np.all(np.isnan(RH)) and meth == "C35"):
RH = np.ones(np.shape(spd))*80 # if empty set to default for COARE3.5
else:
sys.exit("input RH is empty")
logging.debug('input RH is empty')
if (np.all(np.isnan(P)) and meth == "C35"):
P = np.ones(np.shape(spd))*1013 # if empty set to default for COARE3.5
else:
sys.exit("input P is empty")
logging.debug('input P is empty')
if (np.all(np.isnan(Rl)) and meth == "C35"):
Rl = np.ones(np.shape(spd))*370 # set to default for COARE3.5
if (np.all(np.isnan(Rs)) and meth == "C35"):
Rs = np.ones(np.shape(spd))*150 # set to default for COARE3.5
if (np.all(np.isnan(zi)) and meth == "C35"):
zi = 600 # set to default for COARE3.5
elif ((np.all(np.isnan(zi)) and meth == "ERA5") or
(np.all(np.isnan(zi)) and meth == "UA")):
zi = 1000
# additional parameters for C35
waveage, seastate = 1, 1
if np.isnan(wcp[0]):
wcp=np.ones(np.shape(spd))*np.nan
waveage=0
wcp = np.ones(np.shape(spd))*np.nan
waveage = 0
if np.isnan(sigH[0]):
sigH=np.ones(np.shape(spd))*np.nan
seastate=0
sigH = np.ones(np.shape(spd))*np.nan
seastate = 0
if waveage and seastate:
print('Using seastate dependent parameterization')
logging.info('Using seastate dependent parameterization')
......@@ -75,190 +135,298 @@ def flux_calc_v3(spd, T, SST, lat, RH, P, hin, hout, wcp, sigH, zi = 600, \
print('Using waveage dependent parameterization')
logging.info('Using waveage dependent parameterization')
####
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(np.nanmax(SST)<200,np.copy(SST)+CtoK,np.copy(SST))
if (method == "COARE3.5"):
qsea = qsat26sea(sst,P)/1000 # surface water specific humidity (g/kg)
Q,_ = qsat26air(T,P,RH) # specific humidity of air (g/kg)
qair=Q/1000; del Q
logging.info('method %s | qsea:%s, qair:%s',method,np.ma.median(qsea),np.ma.median(qair))
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(np.nanmax(SST) < 200, np.copy(SST)+CtoK, np.copy(SST))
if (meth == "C35"):
qsea = qsat26sea(sst, P)/1000 # surface water specific humidity (g/kg)
Q, _ = qsat26air(T, P, RH) # specific humidity of air (g/kg)
qair = Q/1000
del Q
logging.info('method %s | qsea:%s, qair:%s', meth,
np.ma.median(qsea[~np.isnan(qsea)]),
np.ma.median(qair[~np.isnan(qair)]))
else:
qsea = qsea_calc(sst,P)
qair = q_calc(Ta,RH,P)
logging.info('method %s | qsea:%s, qair:%s',method,np.ma.median(qsea),np.ma.median(qair))
qsea = qsea_calc(sst, P)
qair = q_calc(Ta, RH, P)
logging.info('method %s | qsea:%s, qair:%s', meth,
np.ma.median(qsea[~np.isnan(qsea)]),
np.ma.median(qair[~np.isnan(qair)]))
if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
print("qsea and qair cannot be nan")
logging.debug('method %s qsea and qair cannot be nan | sst:%s, Ta:%s, P:%s, RH:%s',method,np.ma.median(sst),np.ma.median(Ta),np.ma.median(P),np.ma.median(RH))
print("qsea and qair cannot be nan")
logging.debug('method %s qsea and qair cannot be nan | sst:%s, Ta:%s,'
'P:%s, RH:%s', meth, np.ma.median(sst[~np.isnan(sst)]),
np.ma.median(Ta[~np.isnan(Ta)]),
np.ma.median(P[~np.isnan(P)]),
np.ma.median(RH[~np.isnan(RH)]))
# first guesses
inan=np.where(np.isnan(spd+T+SST+lat+RH+P))
inan = np.where(np.isnan(spd+T+SST+lat+RH+P))
t10n, q10n = np.copy(Ta), np.copy(qair)
tv10n = t10n*(1 + 0.61*q10n)
rho = (0.34838*P)/(tv10n)
rhoa = P*100/(287.1*(T+CtoK)*(1+0.61*qair)) # difference with rho is that it uses T instead of Ta (which includes lapse rate)
lv = (2.501-0.00237*(SST))*1e6
dt=sst-Ta
dq=qsea-qair
if (method == "COARE3.5"):
cp=1004.67 #cpa = 1004.67, cpv = cpa*(1+0.84*Q)
rhoa = P*100/(287.1*(T+CtoK)*(1+0.61*qair))
lv = (2.501-0.00237*SST)*1e6
dt = sst - Ta
dq = qsea - qair
if (meth == "C35"):
cp = 1004.67 # cpa = 1004.67, cpv = cpa*(1+0.84*Q)
wetc = 0.622*lv*qsea/(287.1*sst**2)
ug=0.5*np.ones(np.shape(spd))
wind=np.sqrt(np.copy(spd)**2+ug**2)
dter=0.3*np.ones(np.shape(spd)) #cool skin temperature depression
dqer=wetc*dter
Rnl=0.97*(5.67e-8*(SST-dter*jcool+CtoK)**4-Rl)
u10=wind*np.log(10/1e-4)/np.log(hh_in[0]/1e-4)
ug = 0.5*np.ones(np.shape(spd))
wind = np.sqrt(np.copy(spd)**2+ug**2)
dter = 0.3*np.ones(np.shape(spd)) # cool skin temperature depression
dqer = wetc*dter
Rnl = 0.97*(5.67e-8*(SST-dter*jcool+CtoK)**4-Rl)
u10 = wind*np.log(10/1e-4)/np.log(hh_in[0]/1e-4)
usr = 0.035*u10
u10n=np.copy(u10)
zo = 0.011*usr**2/g + 0.11*visc_air(T)/usr
cdn = (kappa/np.log(10/zo))**2
cqn = 0.00115*np.ones(np.shape(spd))
ctn = cqn/np.sqrt(cdn)
cd = (kappa/np.log(hh_in[0]/zo))**2
ct = kappa/np.log(hh_in[1]/10/np.exp(kappa/ctn))
CC = kappa*ct/cd
u10n = np.copy(u10)
zo = 0.011*usr**2/g + 0.11*visc_air(T)/usr
cdn = (kappa/np.log(10/zo))**2
cqn = 0.00115*np.ones(np.shape(spd))
ctn = cqn/np.sqrt(cdn)
cd = (kappa/np.log(hh_in[0]/zo))**2
ct = kappa/np.log(hh_in[1]/10/np.exp(kappa/ctn))
CC = kappa*ct/cd
Ribcu = -hh_in[0]/zi/0.004/1.2**3
Ribu = -g*hh_in[0]/(T+CtoK)*((dt-dter*jcool)+0.61*(T+CtoK)*(dq))/wind**2
zetu=np.where(Ribu<0,CC*Ribu/(1+Ribu/Ribcu),CC*Ribu*(1+27/9*Ribu/CC))
k50=np.where(zetu>50) # stable with very thin M-O length relative to zu
Ribu = (-g*hh_in[0]/(T+CtoK)*((dt-dter*jcool)+0.61*(T+CtoK)*(dq)) /
wind**2)
zetu = np.where(Ribu < 0, CC*Ribu/(1+Ribu/Ribcu),
CC*Ribu*(1+27/9*Ribu/CC))
k50 = np.where(zetu > 50) # stable, very thin M-O relative to zu
monob = hh_in[0]/zetu
gf=wind/spd
gf = wind/spd
usr = wind*kappa/(np.log(hh_in[0]/zo)-psiu_40(hh_in[0]/monob))
tsr = -(dt-dter*jcool)*kappa/(np.log(hh_in[1]/10/np.exp(kappa/ctn))-psit_26(hh_in[1]/monob))
qsr = -(dq-wetc*dter*jcool)*kappa/(np.log(hh_in[2]/10/np.exp(kappa/ctn))-psit_26(hh_in[2]/monob))
tsrv=tsr+0.61*(T+CtoK)*qsr
ac = np.zeros((len(spd),3))
tsr = (-(dt-dter*jcool)*kappa/(np.log(hh_in[1]/10/np.exp(kappa/ctn)) -
psit_26(hh_in[1]/monob)))
qsr = (-(dq-wetc*dter*jcool)*kappa/(np.log(hh_in[2]/10 /
np.exp(kappa/ctn))-psit_26(hh_in[2]/monob)))
tsrv = tsr+0.61*(T+CtoK)*qsr
ac = np.zeros((len(spd), 3))
charn = 0.011*np.ones(np.shape(spd))
charn =np.where(wind>10,0.011+(wind-10)/(18-10)*(0.018-0.011),np.where(wind>18,0.018,0.011))
ac[:,0]=charn
charn = np.where(wind > 10, 0.011+(wind-10)/(18-10)*(0.018-0.011),
np.where(wind > 18, 0.018, 0.011))
ac[:, 0] = charn
else:
cp = 1004.67 * (1 + 0.00084*qsea)
wind,u10n=np.copy(spd),np.copy(spd)
cdn = cdn_calc(u10n,Ta,None,method)
psim, psit, psiq = np.zeros(u10n.shape), np.zeros(u10n.shape), np.zeros(u10n.shape)
psim[inan], psit[inan], psiq[inan] = np.nan,np.nan,np.nan
cd = cd_calc(cdn,hh_in[0],ref_ht,psim)
zo = 0.0001*np.ones(u10n.shape)
zo[inan]=np.nan
monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
cp = 1004.67*(1 + 0.00084*qsea)
u10n = np.copy(spd)
monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
monob[inan] = np.nan
cdn = cdn_calc(u10n, Ta, None, meth)
psim, psit, psiq = (np.zeros(u10n.shape), np.zeros(u10n.shape),
np.zeros(u10n.shape))
psim[inan], psit[inan], psiq[inan] = np.nan, np.nan, np.nan
cd = cd_calc(cdn, hh_in[0], ref_ht, psim)
tsr, tsrv = np.zeros(u10n.shape), np.zeros(u10n.shape)
qsr = np.zeros(u10n.shape)
tsr[inan], tsrv[inan], qsr[inan] = np.nan, np.nan, np.nan
if (meth == "UA"):
wind = np.sqrt(np.power(np.copy(spd),2)+np.power(0.5,2))
zo, zot = 0.0001*np.ones(u10n.shape), 0.0001*np.ones(u10n.shape)
zol = hh_in[0]/monob
elif (meth == "ERA5"):
wind = np.sqrt(np.power(np.copy(spd),2)+np.power(0.5,2))
usr = np.sqrt(cd*wind**2)
Ribu = ((g*hh_in[0]*((2*(Ta-sst))/(Ta+sst-g*hh_in[0]) +
0.61*(qair-qsea)))/np.power(wind, 2))
zo = 0.11*visc_air(Ta)/usr+0.018*np.power(usr, 2)/g
zot = 0.40*visc_air(Ta)/usr
zol = (Ribu*((np.log((hh_in[0]+zo)/zo)-psim_calc((hh_in[0]+zo) /
monob, meth)+psim_calc(zo/monob, meth)) /
(np.log((hh_in[0]+zo)/zot) -
psit_calc((hh_in[0]+zo)/monob, meth) +
psit_calc(zot/monob, meth))))
else:
wind = np.copy(spd)
zo = 0.0001*np.ones(u10n.shape)
zo[inan] = np.nan
usr = np.sqrt(cd * wind**2)
tsr, tsrv, qsr = np.zeros(u10n.shape), np.zeros(u10n.shape), np.zeros(u10n.shape)
tsr[inan], tsrv[inan], qsr[inan] = np.nan, np.nan, np.nan
tol =np.array([0.01,0.01,5e-05,0.005,0.001,5e-07]) # tolerance u,t,q,usr,tsr,qsr
# tolerance for u,t,q,usr,tsr,qsr
tol = np.array([0.01, 0.01, 5e-05, 0.005, 0.001, 5e-07])
it = 0
ind=np.where(spd>0)
ii =True
itera=np.zeros(spd.shape)*np.nan
ind = np.where(spd > 0)
ii = True
itera = np.zeros(spd.shape)*np.nan
while np.any(ii):
it += 1
# print('iter: '+str(it)+', p: '+str(np.shape(ind)[1]))
if it>n:
break
if (method == "COARE3.5"):
### This allows all iterations regardless of convergence as in original
zet=kappa*g*hh_in[0]/(T+CtoK)*(tsr +0.61*(T+CtoK)*qsr)/(usr**2)
monob=hh_in[0]/zet
zo,cd,ct,cq=cd_C35(u10n,wind,usr,ac[:,0],monob,(T+CtoK),hh_in,lat)
usr=wind*cd
qsr=-(dq-wetc*dter*jcool)*cq
tsr=-(dt-dter*jcool)*ct
tsrv=tsr+0.61*(T+CtoK)*qsr
ug = get_gust((T+CtoK),usr,tsrv,zi,lat) # gustiness
wind = np.sqrt(np.copy(spd)**2 + ug**2)
gf=wind/spd
dter, dqer=get_skin(sst,qsea,rhoa,jcool,Rl,Rs,Rnl,cp,lv,usr,tsr,qsr,lat)
Rnl=0.97*(5.67e-8*(SST-dter*jcool+CtoK)**4-Rl)
if it==1: # save first iteration solution for case of zetu>50
usr50,tsr50,qsr50,L50=usr[k50],tsr[k50],qsr[k50],monob[k50]
dter50, dqer50=dter[k50], dqer[k50]
u10n = usr/kappa/gf*np.log(10/zo)
ac=charnock_C35(wind,u10n,usr,seastate,waveage,wcp,sigH,lat)
new=np.array([np.copy(u10n),np.copy(usr),np.copy(tsr),np.copy(qsr)])
### otherwise iterations stop when u,usr,tsr,qsr converge but gives different results than above option
# old=np.array([np.copy(u10n[ind]),np.copy(usr[ind]),np.copy(tsr[ind]),np.copy(qsr[ind])])
# zet=kappa*g[ind]*hh_in[0]/(T[ind]+CtoK)*(tsr[ind] +0.61*(T[ind]+CtoK)*qsr[ind])/(usr[ind]**2)
# monob[ind]=hh_in[0]/zet
# zo[ind],cd[ind],ct[ind],cq[ind]=cd_C35(u10n[ind],wind[ind],usr[ind],ac[ind,0],monob[ind],(T[ind]+CtoK),hh_in,lat[ind])
# usr[ind]=wind[ind]*cd[ind]
# qsr[ind]=-(dq[ind]-wetc[ind]*dter[ind]*jcool)*cq[ind]
# tsr[ind]=-(dt[ind]-dter[ind]*jcool)*ct[ind]
# tsrv[ind]=tsr[ind]+0.61*(T[ind]+CtoK)*qsr[ind]
# ug[ind] = get_gust((T[ind]+CtoK),usr[ind],tsrv[ind],zi,lat[ind]) # gustiness
# wind[ind] = np.sqrt(np.copy(spd[ind])**2 + ug[ind]**2)
# gf=wind/spd
# dter[ind], dqer[ind]=get_skin(sst[ind],qsea[ind],rhoa[ind],jcool,Rl[ind],Rs[ind],Rnl[ind],cp,lv[ind],usr[ind],tsr[ind],qsr[ind],lat[ind])
# Rnl=0.97*(5.67e-8*(SST-dter*jcool+CtoK)**4-Rl)
# if it==1: # save first iteration solution for case of zetu>50
# usr50,tsr50,qsr50,L50=usr[k50],tsr[k50],qsr[k50],monob[k50]
# dter50, dqer50=dter[k50], dqer[k50]
# u10n[ind] = usr[ind]/kappa/gf[ind]*np.log(10/zo[ind])
# ac[ind][:]=charnock_C35(wind[ind],u10n[ind],usr[ind],seastate,waveage,wcp[ind],sigH[ind],lat[ind])
# new=np.array([np.copy(u10n[ind]),np.copy(usr[ind]),np.copy(tsr[ind]),np.copy(qsr[ind])])
# print(str(it)+'L= '+str(monob),file=open('/noc/mpoc/surface_data/sbiri/SAMOS/C35_int.txt','a'))
# print(str(it)+'dter= '+str(dter),file=open('/noc/mpoc/surface_data/sbiri/SAMOS/C35_int.txt','a'))
# print(str(it)+'ac= '+str(ac[:,0]),file=open('/noc/mpoc/surface_data/sbiri/SAMOS/C35_int.txt','a'))
# d=np.abs(new-old)
# ind=np.where((d[0,:]>tol[0])+(d[1,:]>tol[3])+(d[2,:]>tol[4])+(d[3,:]>tol[5]))
# itera[ind]=np.ones(1)*it
# if np.shape(ind)[0]==0:
# break
# else:
# ii = (d[0,:]>tol[0])+(d[1,:]>tol[3])+(d[2,:]>tol[4])+(d[3,:]>tol[5])
if it > n:
break
if (meth == "C35"):
old = np.array([np.copy(u10n[ind]), np.copy(usr[ind]),
np.copy(tsr[ind]), np.copy(qsr[ind])])
zet = (kappa*g[ind]*hh_in[0]/(T[ind]+CtoK)*(tsr[ind]+0.61 *
(T[ind]+CtoK)*qsr[ind])/(usr[ind]**2))
monob[ind] = hh_in[0]/zet
zo[ind], cd[ind], ct[ind], cq[ind] = cd_C35(u10n[ind], wind[ind],
usr[ind], ac[ind,0],
monob[ind],
(T[ind]+CtoK), hh_in,
lat[ind])
usr[ind] = wind[ind]*cd[ind]
qsr[ind] = -(dq[ind]-wetc[ind]*dter[ind]*jcool)*cq[ind]
tsr[ind] = -(dt[ind]-dter[ind]*jcool)*ct[ind]
tsrv[ind] = tsr[ind]+0.61*(T[ind]+CtoK)*qsr[ind]
ug[ind] = get_gust(1.2, (T[ind]+CtoK), usr[ind], tsrv[ind],
zi, lat[ind]) # gustiness
wind[ind] = np.sqrt(np.copy(spd[ind])**2 + ug[ind]**2)
gf = wind/spd
dter[ind], dqer[ind] = get_skin(sst[ind], qsea[ind], rhoa[ind],
Rl[ind], Rs[ind], Rnl[ind], cp,
lv[ind], usr[ind], tsr[ind],
qsr[ind], lat[ind])
Rnl = 0.97*(5.67e-8*(SST-dter*jcool+CtoK)**4-Rl)
if it == 1: # save first iteration solution for case of zetu>50
usr50, tsr50, qsr50, L50 = (usr[k50], tsr[k50], qsr[k50],
monob[k50])
dter50, dqer50 = dter[k50], dqer[k50]
u10n[ind] = usr[ind]/kappa/gf[ind]*np.log(10/zo[ind])
ac[ind][:] = charnock_C35(wind[ind], u10n[ind], usr[ind], seastate,
waveage, wcp[ind], sigH[ind], lat[ind])
new = np.array([np.copy(u10n[ind]), np.copy(usr[ind]),
np.copy(tsr[ind]), np.copy(qsr[ind])])
d = np.abs(new-old)
ind = np.where((d[0,:] > tol[0])+(d[1,:] > tol[3]) +
(d[2,:] > tol[4])+(d[3,:] > tol[5]))
itera[ind] = np.ones(1)*it
if np.shape(ind)[0] == 0:
break
else:
ii = ((d[0,:] > tol[0])+(d[1,:] > tol[3])+(d[2,:] > tol[4]) +
(d[3,:] > tol[5]))
else:
old=np.array([np.copy(u10n[ind]),np.copy(t10n[ind]),np.copy(q10n[ind]),np.copy(usr[ind]),np.copy(tsr[ind]),np.copy(qsr[ind])])
cdn[ind] = cdn_calc(u10n[ind],Ta[ind],None,method)
old = np.array([np.copy(u10n[ind]), np.copy(t10n[ind]),
np.copy(q10n[ind]), np.copy(usr[ind]),
np.copy(tsr[ind]), np.copy(qsr[ind])])
cdn[ind] = cdn_calc(u10n[ind], Ta[ind], None, meth)
if (np.all(np.isnan(cdn))):
break #sys.exit("cdn cannot be nan")
logging.debug('%s at iteration %s cdn negative',method,it)
break # sys.exit("cdn cannot be nan")
logging.debug('break %s at iteration %s cdn<0', meth, it)
zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cdn[ind]))
psim[ind] = psim_calc(hh_in[0]/monob[ind],method)
psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
cd[ind] = cd_calc(cdn[ind], hh_in[0], ref_ht, psim[ind])
ctn[ind], cqn[ind] = ctcqn_calc(hh_in[1]/monob[ind],cdn[ind],u10n[ind],zo[ind],Ta[ind],method)
psit[ind] = psit_calc(hh_in[1]/monob[ind],method)
psiq[ind] = psit_calc(hh_in[2]/monob[ind],method)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind], hh_in[1], hh_in[2], ref_ht, psit[ind], psiq[ind])
usr[ind] = np.sqrt(cd[ind]*wind[ind]**2)
ctn[ind], cqn[ind] = ctcqn_calc(hh_in[1]/monob[ind], cdn[ind],
u10n[ind], zo[ind], Ta[ind], meth)
psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
ct[ind], cq[ind] = ctcq_calc(cdn[ind], cd[ind], ctn[ind], cqn[ind],
hh_in[1], hh_in[2], ref_ht,
psit[ind], psiq[ind])
usr[ind] = np.sqrt(cd[ind]*wind[ind]**2)
tsr[ind] = ct[ind]*wind[ind]*(-dt[ind])/usr[ind]
qsr[ind] = cq[ind]*wind[ind]*(-dq[ind])/usr[ind]
fact = (np.log(hh_in[1]/ref_ht)-psit[ind])/kappa
t10n[ind] = Ta[ind] - (tsr[ind]*fact)
fact = (np.log(hh_in[2]/ref_ht)-psiq[ind])/kappa
q10n[ind] = qair[ind] - (qsr[ind]*fact)
if (meth == "UA"):
zot[ind] = 10/(np.exp(kappa**2/(ctn[ind]*np.log(10/zo[ind]))))
zol[ind] = hh_in[0]/monob[ind]
t10n[ind] = np.where(zol[ind]<-0.465,Ta[ind]-(tsr[ind]/kappa) *
(np.log((-0.465*monob[ind])/zot[ind]) -
psit_calc(-0.465,meth)+0.8 *
(np.power(0.465,-1/3) -
np.power(-zol[ind],-1/3))), \
np.where((zol[ind]>-0.465) & (zol[ind]<0),
Ta[ind]-(tsr[ind]/kappa) *
(np.log(hh_in[1]/zot[ind]) -
psit_calc(zol[ind],meth)), \
np.where((zol[ind]>0) & (zol[ind]<1),
Ta[ind]-(tsr[ind]/kappa) *
(np.log(hh_in[1]/zot[ind])+5*monob[ind]),
Ta[ind]-(tsr[ind]/kappa) *
(np.log(monob[ind]/zot[ind])+5 +
5*np.log(monob[ind])+monob[ind]-1))))
# Zeng et al. 1998 (11-14)
q10n[ind] = np.where(zol[ind]<-0.465,qair[ind] -
(qsr[ind]/kappa) *
(np.log((-0.465*monob[ind])/zot[ind]) -
psit_calc(-0.465,meth)+0.8 *
(np.power(0.465,-1/3) -
np.power(-zol[ind],-1/3))), \
np.where((zol[ind]>-0.465) & (zol[ind]<0),
qair[ind]-(qsr[ind]/kappa) *
(np.log(hh_in[1]/zot[ind]) -
psit_calc(zol[ind],meth)), \
np.where((zol[ind]>0) & (zol[ind]<1), \
qair[ind]-(qsr[ind]/kappa) *
(np.log(hh_in[1]/zot[ind])+5*monob[ind]),
qair[ind]-(qsr[ind]/kappa) *
(np.log(monob[ind]/zot[ind])+5 +
5*np.log(monob[ind])+monob[ind]-1))))
else:
fact = (np.log(hh_in[1]/ref_ht)-psit[ind])/kappa
t10n[ind] = Ta[ind] - (tsr[ind]*fact)
fact = (np.log(hh_in[2]/ref_ht)-psiq[ind])/kappa
q10n[ind] = qair[ind] - (qsr[ind]*fact)
tv10n[ind] = t10n[ind]*(1+0.61*q10n[ind])
tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
monob[ind] = (tv10n[ind]*usr[ind]**2)/(g[ind]*kappa*tsrv[ind])
psim[ind] = psim_calc(hh_in[0]/monob[ind],method)
psit[ind] = psit_calc(hh_in[1]/monob[ind],method)
psiq[ind] = psit_calc(hh_in[2]/monob[ind],method)
u10n[ind] = wind[ind] -(usr[ind]/kappa)*(np.log(hh_in[0]/ref_ht)-psim[ind])
u10n[u10n<0]=np.nan # 0.5*old[0,np.where(u10n<0)]
new=np.array([np.copy(u10n[ind]),np.copy(t10n[ind]),np.copy(q10n[ind]),np.copy(usr[ind]),np.copy(tsr[ind]),np.copy(qsr[ind])])
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]))
itera[ind]=np.ones(1)*it
if np.shape(ind)[0]==0:
break
if (meth == "ERA5"):
zol[ind] = (kappa*g[ind]*hh_in[0]/Ta[ind]*(tsr[ind] +
0.61*Ta[ind]*qsr[ind])/(usr[ind]**2))
monob[ind] = hh_in[0]/zol[ind]
else:
monob[ind] = (tv10n[ind]*usr[ind]**2)/(g[ind]*kappa*tsrv[ind])
psim[ind] = psim_calc(hh_in[0]/monob[ind], meth)
psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
if (meth == "UA"):
u10n[ind] = np.where(zol[ind]<-1.574,(usr[ind]/kappa) *
((np.log(-1.574*monob[ind]/zo[ind]) -
psim_calc(-1.574,meth))+1.14 *
(np.power(-monob[ind],1/3) -
np.power(1.574,1/3))), \
np.where((zol[ind]>-1.574) & (zol[ind]<0),
(usr[ind]/kappa)*(np.log(ref_ht/zo[ind]) -
psim_calc(zol[ind])), \
np.where((zol[ind]>0) & (zol[ind]<1),
(usr[ind]/kappa)*(np.log(ref_ht/zo[ind]) +
5*zol[ind]),(usr[ind]/kappa) *
(np.log(monob[ind]/zo[ind])+5 +
5*np.log(zol[ind])+zol[ind]-1))))
# Zeng et al. 1998 (7-10)
u10n[u10n<0] = np.nan
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]),2) +
np.power(get_gust(1,Ta[ind], usr[ind],
tsrv[ind],zi,lat[ind]),2))
# Zeng et al. 1998 (20)
elif (meth == "ERA5"):
wind[ind] = np.sqrt(np.copy(spd[ind])**2 +
(get_gust(1, Ta[ind], usr[ind], tsrv[ind],
zi, lat[ind]))**2)
u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0] /
ref_ht)-psim[ind]))
u10n[u10n < 0] = np.nan
else:
ii = (d[0,ind]>tol[0])+(d[1,ind]>tol[1])+(d[2,ind]>tol[2])+(d[3,ind]>tol[3])+(d[4,ind]>tol[4])+(d[5,ind]>tol[5])
u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0] /
ref_ht)-psim[ind]))
u10n[u10n < 0] = np.nan # 0.5*old[0,np.where(u10n<0)]
new = np.array([np.copy(u10n[ind]), np.copy(t10n[ind]),
np.copy(q10n[ind]), np.copy(usr[ind]),
np.copy(tsr[ind]), np.copy(qsr[ind])])
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]))
itera[ind] = np.ones(1)*it
if np.shape(ind)[0] == 0:
break
else:
ii = ((d[0, ind] > tol[0])+(d[1, ind] > tol[1]) +
(d[2, ind] > tol[2])+(d[3, ind] > tol[3]) +
(d[4, ind] > tol[4])+(d[5, ind] > tol[5]))
# calculate output parameters
if (method == "COARE3.5"):
usr[k50],tsr[k50],qsr[k50],monob[k50]=usr50,tsr50,qsr50,L50
dter[k50], dqer[k50]=dter50, dqer50
if (meth == "C35"):
usr[k50], tsr[k50], qsr[k50], monob[k50] = usr50, tsr50, qsr50, L50
dter[k50], dqer[k50] = dter50, dqer50
rhoa = P*100/(287.1*(T+CtoK)*(1+0.61*qair))
rr=zo*usr/visc_air(T+CtoK)
zoq=np.where(5.8e-5/rr**0.72>1.6e-4,1.6e-4,5.8e-5/rr**0.72)
zot=zoq
psiT=psit_26(hh_in[1]/monob)
psi10T=psit_26(10/monob)
psi=psiu_26(hh_in[0]/monob)
psirf=psiu_26(hh_out[0]/monob)
rr = zo*usr/visc_air(T+CtoK)
zoq = np.where(5.8e-5/rr**0.72 > 1.6e-4, 1.6e-4, 5.8e-5/rr**0.72)
zot = zoq
psiT = psit_26(hh_in[1]/monob)
psi10T = psit_26(10/monob)
psi = psiu_26(hh_in[0]/monob)
psirf = psiu_26(hh_out[0]/monob)
q10 = qair + qsr/kappa*(np.log(10/hh_in[2])-psi10T+psiT)
tau = rhoa*usr*usr/gf
sensible=-rhoa*cp*usr*tsr
latent=-rhoa*lv*usr*qsr
cd = tau/rhoa/wind/np.where(spd<0.1,0.1,spd)
sensible = -rhoa*cp*usr*tsr
latent = -rhoa*lv*usr*qsr
cd = tau/rhoa/wind/np.where(spd < 0.1, 0.1, spd)
cdn = 1000*kappa**2/np.log(10/zo)**2
ct = -usr*tsr/wind/(dt-dter*jcool)
ctn = 1000*kappa**2/np.log(10/zo)/np.log(10/zot)
......@@ -267,12 +435,15 @@ def flux_calc_v3(spd, T, SST, lat, RH, P, hin, hout, wcp, sigH, zi = 600, \
psim = psiu_26(hh_in[0]/monob)
psit = psit_26(hh_in[1]/monob)
u10n = u10+psiu_26(10/monob)*usr/kappa/gf
t10n = T + tsr/kappa*(np.log(10/hh_in[1])-psi10T+psiT) + tlapse*(hh_in[1]-10)+ psi10T*tsr/kappa
t10n = (T+tsr/kappa*(np.log(10/hh_in[1])-psi10T+psiT) +
tlapse*(hh_in[1]-10)+psi10T*tsr/kappa)
q10n = q10 + psi10T*qsr/kappa
tv10n = t10n*(1+0.61*q10n)
urefs = spd + usr/kappa/gf*(np.log(hh_out[0]/hh_in[0])-psirf+psi)
trefs =T + tsr/kappa*(np.log(hh_out[1]/hh_in[1])-psit_26(hh_out[1]/monob)+psiT) + tlapse*(hh_in[1]-hh_out[1])
qrefs= qair + qsr/kappa*(np.log(hh_out[2]/hh_in[2])-psit_26(hh_out[2]/monob)+psiT)
trefs = (T+tsr/kappa*(np.log(hh_out[1]/hh_in[1]) -
psit_26(hh_out[1]/monob)+psiT)+tlapse*(hh_in[1]-hh_out[1]))
qrefs = (qair+qsr/kappa*(np.log(hh_out[2]/hh_in[2]) -
psit_26(hh_out[2]/monob)+psiT))
else:
rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht)
......@@ -282,59 +453,39 @@ def flux_calc_v3(spd, T, SST, lat, RH, P, hin, hout, wcp, sigH, zi = 600, \
zo = ref_ht/np.exp(kappa/cdn**0.5)
zot = ref_ht/(np.exp(kappa**2/(ctn*np.log(ref_ht/zo))))
zoq = ref_ht/(np.exp(kappa**2/(cqn*np.log(ref_ht/zo))))
urefs=spd-(usr/kappa)*(np.log(hh_in[0]/hh_out[0])-psim+psim_calc(hh_out[0]/monob,method))
trefs=Ta-(tsr/kappa)*(np.log(hh_in[1]/hh_out[1])-psit+psit_calc(hh_out[0]/monob,method))
urefs = (spd-(usr/kappa)*(np.log(hh_in[0]/hh_out[0])-psim +
psim_calc(hh_out[0]/monob, meth)))
trefs = (Ta-(tsr/kappa)*(np.log(hh_in[1]/hh_out[1])-psit +
psit_calc(hh_out[0]/monob, meth)))
trefs = trefs-(273.16+tlapse*hh_out[1])
qrefs=qair-(qsr/kappa)*(np.log(hh_in[2]/hh_out[2])-psit+psit_calc(hh_out[2]/monob,method))
res=np.zeros((27,len(spd)))
res[0][:]=tau
res[1][:]=sensible
res[2][:]=latent
res[3][:]=monob
res[4][:]=cd
res[5][:]=cdn
res[6][:]=ct
res[7][:]=ctn
res[8][:]=cq
res[9][:]=cqn
res[10][:]=tsrv
res[11][:]=tsr
res[12][:]=qsr
res[13][:]=usr
res[14][:]=psim
res[15][:]=psit
res[16][:]=u10n
res[17][:]=t10n
res[18][:]=tv10n
res[19][:]=q10n
res[20][:]=zo
res[21][:]=zot
res[22][:]=zoq
res[23][:]=urefs
res[24][:]=trefs
res[25][:]=qrefs
res[26][:]=itera
qrefs = (qair-(qsr/kappa)*(np.log(hh_in[2]/hh_out[2]) -
psit+psit_calc(hh_out[2]/monob, meth)))
res = np.zeros((27, len(spd)))
res[0][:] = tau
res[1][:] = sensible
res[2][:] = latent
res[3][:] = monob
res[4][:] = cd
res[5][:] = cdn
res[6][:] = ct
res[7][:] = ctn
res[8][:] = cq
res[9][:] = cqn
res[10][:] = tsrv
res[11][:] = tsr
res[12][:] = qsr
res[13][:] = usr
res[14][:] = psim
res[15][:] = psit
res[16][:] = u10n
res[17][:] = t10n
res[18][:] = tv10n
res[19][:] = q10n
res[20][:] = zo
res[21][:] = zot
res[22][:] = zoq
res[23][:] = urefs
res[24][:] = trefs
res[25][:] = qrefs
res[26][:] = itera
return res, ind
def charnock_C35(wind,u10n,usr,seastate,waveage,wcp,sigH,lat):
g=gc(lat,None)
a1, a2=0.0017, -0.0050
charnC=np.where(u10n>19,a1*19+a2,a1*u10n+a2)
A, B=0.114, 0.622 #wave-age dependent coefficients
Ad, Bd=0.091, 2.0 #Sea-state/wave-age dependent coefficients
charnW=A*(usr/wcp)**B
zoS=sigH*Ad*(usr/wcp)**Bd
charnS=(zoS*g)/usr**2
charn=np.where(wind>10,0.011+(wind-10)/(18-10)*(0.018-0.011),np.where(wind>18,0.018,0.011*np.ones(np.shape(wind))))
if waveage:
if seastate:
charn=charnS
else:
charn=charnW
else:
charn=charnC
ac = np.zeros((len(wind),3))
ac[:,0] = charn
ac[:,1] = charnC
ac[:,2] = charnW
return ac
\ No newline at end of file
flux_subs.py
View file @ 22be86dd
import numpy as np
""" Conversion factor for [:math:`^\\circ` C] to [:math:`^\\circ` K] """
CtoK = 273.16 # 273.15
""" von Karman's constant """
""" Conversion factor for (^\circ\,C) to (^\\circ\\,K) """
kappa = 0.4 # NOTE: 0.41
""" von Karman's constant """
# ---------------------------------------------------------------------
def charnock_C35(wind, u10n, usr, seastate, waveage, wcp, sigH, lat):
""" Calculates Charnock number following Edson et al. 2013 based on
C35 matlab code (coare35vn.m)
Parameters
----------
wind : float
wind speed (m/s)
u10n : float
neutral 10m wind speed (m/s)
usr : float
friction velocity (m/s)
seastate : bool
0 or 1
waveage : bool
0 or 1
wcp : float
phase speed of dominant waves (m/s)
sigH : float
significant wave height (m)
lat : float
latitude (deg)
Returns
-------
ac : float
Charnock number
"""
g = gc(lat, None)
a1, a2 = 0.0017, -0.0050
charnC = np.where(u10n > 19, a1*19+a2, a1*u10n+a2)
......@@ -34,6 +62,37 @@ def charnock_C35(wind, u10n, usr, seastate, waveage, wcp, sigH, lat):
def cd_C35(u10n, wind, usr, charn, monob, Ta, hh_in, lat):
""" Calculates exchange coefficients following Edson et al. 2013 based on
C35 matlab code (coare35vn.m)
Parameters
----------
u10n : float
neutral 10m wind speed (m/s)
wind : float
wind speed (m/s)
charn : float
Charnock number
monob : float
Monin-Obukhov stability length
Ta : float
air temperature (K)
hh_in : float
input sensor's height (m)
lat : float
latitude (deg)
Returns
-------
zo : float
surface roughness (m)
cdhf : float
drag coefficient
cthf : float
heat exchange coefficient
cqhf : float
moisture exchange coefficient
"""
g = gc(lat, None)
zo = charn*usr**2/g+0.11*visc_air(Ta)/usr # surface roughness
rr = zo*usr/visc_air(Ta)
......@@ -47,15 +106,31 @@ def cd_C35(u10n, wind, usr, charn, monob, Ta, hh_in, lat):
# ---------------------------------------------------------------------
def cdn_calc(u10n, Ta, Tp, method="Smith80"):
if (method == "Smith80"):
def cdn_calc(u10n, Ta, Tp, method="S80"):
""" Calculates neutral drag coefficient
Parameters
----------
u10n : float
neutral 10m wind speed (m/s)
Ta : float
air temperature (K)
Tp : float
wave period
method : str
Returns
-------
cdn : float
"""
if (method == "S80"):
cdn = np.where(u10n <= 3, (0.61+0.567/u10n)*0.001,
(0.61+0.063*u10n)*0.001)
elif (method == "LP82"):
cdn = np.where((u10n < 11) & (u10n >= 4), 1.2*0.001,
np.where((u10n <= 25) & (u10n >= 11),
(0.49+0.065*u10n)*0.001, 1.14*0.001))
elif (method == "Smith88" or method == "COARE3.0" or
elif (method == "S88" or method == "C30" or
method == "COARE4.0" or method == "UA" or method == "ERA5"):
cdn = cdn_from_roughness(u10n, Ta, None, method)
elif (method == "HEXOS"):
......@@ -64,7 +139,7 @@ def cdn_calc(u10n, Ta, Tp, method="Smith80"):
elif (method == "HEXOSwave"):
cdn = cdn_from_roughness(u10n, Ta, Tp, method)
elif (method == "YT96"):
# for u<3 same as Smith80
# for u<3 same as S80
cdn = np.where((u10n < 6) & (u10n >= 3),
(0.29+3.1/u10n+7.7/u10n**2)*0.001,
np.where((u10n <= 26) & (u10n >= 6),
......@@ -78,7 +153,23 @@ def cdn_calc(u10n, Ta, Tp, method="Smith80"):
# ---------------------------------------------------------------------
def cdn_from_roughness(u10n, Ta, Tp, method="Smith88"):
def cdn_from_roughness(u10n, Ta, Tp, method="S88"):
""" Calculates neutral drag coefficient from roughness length
Parameters
----------
u10n : float
neutral 10m wind speed (m/s)
Ta : float
air temperature (K)
Tp : float
wave period
method : str
Returns
-------
cdn : float
"""
g, tol = 9.812, 0.000001
cdn, usr = np.zeros(Ta.shape), np.zeros(Ta.shape)
cdnn = (0.61+0.063*u10n)*0.001
......@@ -86,13 +177,13 @@ def cdn_from_roughness(u10n, Ta, Tp, method="Smith88"):
for it in range(5):
cdn = np.copy(cdnn)
usr = np.sqrt(cdn*u10n**2)
if (method == "Smith88"):
if (method == "S88"):
# .....Charnock roughness length (equn 4 in Smith 88)
zc = 0.011*np.power(usr, 2)/g
# .....smooth surface roughness length (equn 6 in Smith 88)
zs = 0.11*visc_air(Ta)/usr
zo = zc + zs # .....equns 7 & 8 in Smith 88 to calculate new CDN
elif (method == "COARE3.0"):
elif (method == "C30"):
zc = 0.011 + (u10n-10)/(18-10)*(0.018-0.011)
zc = np.where(u10n < 10, 0.011, np.where(u10n > 18, 0.018, zc))
zs = 0.11*visc_air(Ta)/usr
......@@ -116,9 +207,32 @@ def cdn_from_roughness(u10n, Ta, Tp, method="Smith88"):
# ---------------------------------------------------------------------
def ctcqn_calc(zol, cdn, u10n, zo, Ta, method="Smith80"):
if (method == "Smith80" or method == "Smith88" or method == "YT96"):
cqn = np.ones(Ta.shape)*1.20*0.001 # from Smith88
def ctcqn_calc(zol, cdn, u10n, zo, Ta, method="S80"):
""" Calculates neutral heat and moisture exchange coefficients
Parameters
----------
zol : float
height over MO length
cdn : float
neatral drag coefficient
u10n : float
neutral 10m wind speed (m/s)
zo : float
surface roughness (m)
Ta : float
air temperature (K)
method : str
Returns
-------
ctn : float
neutral heat exchange coefficient
cqn : float
neutral moisture exchange coefficient
"""
if (method == "S80" or method == "S88" or method == "YT96"):
cqn = np.ones(Ta.shape)*1.20*0.001 # from S88
ctn = np.ones(Ta.shape)*1.00*0.001
elif (method == "LP82"):
cqn = np.where((zol <= 0) & (u10n > 4) & (u10n < 14), 1.15*0.001,
......@@ -128,7 +242,7 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, method="Smith80"):
elif (method == "HEXOS" or method == "HEXOSwave"):
cqn = np.where((u10n <= 23) & (u10n >= 3), 1.1*0.001, np.nan)
ctn = np.where((u10n <= 18) & (u10n >= 3), 1.1*0.001, np.nan)
elif (method == "COARE3.0" or method == "COARE4.0"):
elif (method == "C30" or method == "COARE4.0"):
usr = (cdn*u10n**2)**0.5
rr = zo*usr/visc_air(Ta)
zoq = 5.5e-5/rr**0.6
......@@ -162,6 +276,23 @@ def ctcqn_calc(zol, cdn, u10n, zo, Ta, method="Smith80"):
def cd_calc(cdn, height, ref_ht, psim):
""" Calculates drag coefficient at reference height
Parameters
----------
cdn : float
neutral drag coefficient
height : float
original sensor height (m)
ref_ht : float
reference height (m)
psim : float
momentum stability function
Returns
-------
cd : float
"""
cd = (cdn*np.power(1+np.sqrt(cdn)*np.power(kappa, -1) *
(np.log(height/ref_ht)-psim), -2))
return cd
......@@ -169,16 +300,56 @@ def cd_calc(cdn, height, ref_ht, psim):
def ctcq_calc(cdn, cd, ctn, cqn, h_t, h_q, ref_ht, psit, psiq):
""" Calculates heat and moisture exchange coefficients at reference height
Parameters
----------
cdn : float
neutral drag coefficient
cd : float
drag coefficient at reference height
ctn : float
neutral heat exchange coefficient
cqn : float
neutral moisture exchange coefficient
h_t : float
original temperature sensor height (m)
h_q : float
original moisture sensor height (m)
ref_ht : float
reference height (m)
psit : float
heat stability function
psiq : float
moisture stability function
Returns
-------
ct : float
cq : float
"""
ct = ctn*(cd/cdn)**0.5/(1+ctn*((np.log(h_t/ref_ht)-psit)/(kappa*cdn**0.5)))
cq = cqn*(cd/cdn)**0.5/(1+cqn*((np.log(h_q/ref_ht)-psiq)/(kappa*cdn**0.5)))
return ct, cq
# ---------------------------------------------------------------------
def psim_calc(zol, method="Smith80"):
def psim_calc(zol, method="S80"):
""" Calculates momentum stability function
Parameters
----------
zol : float
height over MO length
method : str
Returns
-------
psim : float
"""
coeffs = get_stabco(method)
alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (method == "COARE3.0" or method == "COARE4.0"):
if (method == "C30" or method == "COARE4.0"):
psim = np.where(zol < 0, psim_conv_coare3(zol, alpha, beta, gamma),
psim_stab_coare3(zol, alpha, beta, gamma))
elif (method == "ERA5"):
......@@ -191,10 +362,22 @@ def psim_calc(zol, method="Smith80"):
# ---------------------------------------------------------------------
def psit_calc(zol, method="Smith80"):
def psit_calc(zol, method="S80"):
""" Calculates heat stability function
Parameters
----------
zol : float
height over MO length
method : str
Returns
-------
psit : float
"""
coeffs = get_stabco(method)
alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (method == "COARE3.0" or method == "COARE4.0"):
if (method == "C30" or method == "COARE4.0"):
psit = np.where(zol < 0, psi_conv_coare3(zol, alpha, beta, gamma),
psi_stab_coare3(zol, alpha, beta, gamma))
elif (method == "ERA5"):
......@@ -207,8 +390,18 @@ def psit_calc(zol, method="Smith80"):
# ---------------------------------------------------------------------
def get_stabco(method="Smith80"):
if (method == "Smith80" or method == "Smith88" or method == "LY04" or
def get_stabco(method="S80"):
""" Gives the coefficients \\alpha, \\beta, \\gamma for stability functions
Parameters
----------
method : str
Returns
-------
coeffs : float
"""
if (method == "S80" or method == "S88" or method == "LY04" or
method == "UA" or method == "ERA5"):
alpha, beta, gamma = 16, 0.25, 5 # Smith 1980, from Dyer (1974)
elif (method == "LP82"):
......@@ -217,7 +410,7 @@ def get_stabco(method="Smith80"):
alpha, beta, gamma = 16, 0.25, 8
elif (method == "YT96"):
alpha, beta, gamma = 20, 0.25, 5
elif (method == "COARE3.0" or method == "COARE4.0"):
elif (method == "C30" or method == "COARE4.0"):
# use separate subroutine
alpha, beta, gamma = 15, 1/3, 5 # not sure about gamma=34.15
else:
......@@ -231,6 +424,22 @@ def get_stabco(method="Smith80"):
def psi_conv_coare3(zol, alpha, beta, gamma):
""" Calculates heat stability function for unstable conditions
for method C30
Parameters
----------
zol : float
height over MO length
alpha : float
beta : float
gamma : float
constants given by get_stabco
Returns
-------
psit : float
"""
x = (1-alpha*zol)**0.5 # Kansas unstable
psik = 2*np.log((1+x)/2.)
y = (1-34.15*zol)**beta
......@@ -242,7 +451,21 @@ def psi_conv_coare3(zol, alpha, beta, gamma):
# ---------------------------------------------------------------------
def psi_stab_coare3(zol, alpha, beta, gamma): # Stable
def psi_stab_coare3(zol, alpha, beta, gamma):
""" Calculates heat stability function for stable conditions
for method C30
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psi : float
"""
c = np.where(0.35*zol > 50, 50, 0.35*zol) # Stable
psit = -((1+2*zol/3)**1.5+0.6667*(zol-14.28)/np.exp(c)+8.525)
return psit
......@@ -250,6 +473,20 @@ def psi_stab_coare3(zol, alpha, beta, gamma): # Stable
def psi_stab_era5(zol, alpha, beta, gamma):
""" Calculates heat stability function for stable conditions
for method ERA5
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psit : float
"""
# eq (3.22) p. 39 IFS Documentation cy46r1
a, b, c, d = 1, 2/3, 5, 0.35
psit = -b*(zol-c/d)*np.exp(-d*zol)-np.power(1+(2/3)*a*zol, 1.5)-(b*c)/d+1
......@@ -257,6 +494,19 @@ def psi_stab_era5(zol, alpha, beta, gamma):
# ---------------------------------------------------------------------
def psi_conv(zol, alpha, beta, gamma):
""" Calculates heat stability function for unstable conditions
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psit : float
"""
xtmp = (1-alpha*zol)**beta
psit = 2*np.log((1+xtmp**2)*0.5)
return psit
......@@ -264,30 +514,65 @@ def psi_conv(zol, alpha, beta, gamma):
def psi_stab(zol, alpha, beta, gamma):
""" Calculates heat stability function for stable conditions
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psit : float
"""
psit = -gamma*zol
return psit
# ---------------------------------------------------------------------
def psit_26(zet):
"""
computes temperature structure function as in COARE3.5
def psit_26(zol):
""" Computes temperature structure function as in C35
Parameters
----------
zol : float
height over MO length
Returns
-------
psi : float
"""
dzet = np.where(0.35*zet > 50, 50, 0.35*zet) # stable
psi = -((1+0.6667*zet)**1.5+0.6667*(zet-14.28)*np.exp(-dzet)+8.525)
k = np.where(zet < 0) # unstable
x = (1-15*zet[k])**0.5
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
psi = -((1+0.6667*zol)**1.5+0.6667*(zol-14.28)*np.exp(-dzol)+8.525)
k = np.where(zol < 0) # unstable
x = (1-15*zol[k])**0.5
psik = 2*np.log((1+x)/2)
x = (1-34.15*zet[k])**0.3333
x = (1-34.15*zol[k])**0.3333
psic = (1.5*np.log((1+x+x**2)/3)-np.sqrt(3)*np.arctan((1+2*x) /
np.sqrt(3))+4*np.arctan(1)/np.sqrt(3))
f = zet[k]**2/(1+zet[k]**2)
f = zol[k]**2/(1+zol[k]**2)
psi[k] = (1-f)*psik+f*psic
return psi
# ---------------------------------------------------------------------
def psim_conv_coare3(zol, alpha, beta, gamma):
""" Calculates momentum stability function for unstable conditions
for method C30
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psim : float
"""
x = (1-15*zol)**0.25 # Kansas unstable
psik = 2*np.log((1+x)/2)+np.log((1+x*x)/2)-2*np.arctan(x)+2*np.arctan(1)
y = (1-10.15*zol)**0.3333 # Convective
......@@ -300,6 +585,20 @@ def psim_conv_coare3(zol, alpha, beta, gamma):
def psim_stab_coare3(zol, alpha, beta, gamma):
""" Calculates momentum stability function for stable conditions
for method C30
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psim : float
"""
c = np.where(0.35*zol > 50, 50, 0.35*zol) # Stable
psim = -((1+1*zol)**1.0+0.6667*(zol-14.28)/np.exp(-c)+8.525)
return psim
......@@ -307,6 +606,20 @@ def psim_stab_coare3(zol, alpha, beta, gamma):
def psim_stab_era5(zol, alpha, beta, gamma):
""" Calculates momentum stability function for stable conditions
for method ERA5
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psim : float
"""
# eq (3.22) p. 39 IFS Documentation cy46r1
a, b, c, d = 1, 2/3, 5, 0.35
psim = -b*(zol-c/d)*np.exp(-d*zol)-a*zol-(b*c)/d
......@@ -315,6 +628,19 @@ def psim_stab_era5(zol, alpha, beta, gamma):
def psim_conv(zol, alpha, beta, gamma):
""" Calculates momentum stability function for unstable conditions
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psim : float
"""
xtmp = (1-alpha*zol)**beta
psim = (2*np.log((1+xtmp)*0.5)+np.log((1+xtmp**2)*0.5) -
2*np.arctan(xtmp)+np.pi/2)
......@@ -323,50 +649,112 @@ def psim_conv(zol, alpha, beta, gamma):
def psim_stab(zol, alpha, beta, gamma):
""" Calculates momentum stability function for stable conditions
Parameters
----------
zol : float
height over MO length
alpha, beta, gamma : float
constants given by get_stabco
Returns
-------
psim : float
"""
psim = -gamma*zol
return psim
# ---------------------------------------------------------------------
def psiu_26(zet):
def psiu_26(zol):
""" Computes velocity structure function C35
Parameters
----------
zol : float
height over MO length
Returns
-------
psi : float
"""
computes velocity structure function COARE3.5
"""
dzet = np.where(0.35*zet > 50, 50, 0.35*zet) # stable
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
a, b, c, d = 0.7, 3/4, 5, 0.35
psi = -(a*zet+b*(zet-c/d)*np.exp(-dzet)+b*c/d)
k = np.where(zet < 0) # unstable
x = (1-15*zet[k])**0.25
psi = -(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d)
k = np.where(zol < 0) # unstable
x = (1-15*zol[k])**0.25
psik = 2*np.log((1+x)/2)+np.log((1+x**2)/2)-2*np.arctan(x)+2*np.arctan(1)
x = (1-10.15*zet[k])**0.3333
x = (1-10.15*zol[k])**0.3333
psic = (1.5*np.log((1+x+x**2)/3)-np.sqrt(3)*np.arctan((1+2*x)/np.sqrt(3)) +
4*np.arctan(1)/np.sqrt(3))
f = zet[k]**2/(1+zet[k]**2)
f = zol[k]**2/(1+zol[k]**2)
psi[k] = (1-f)*psik+f*psic
return psi
# ------------------------------------------------------------------------------
def psiu_40(zet):
"""
computes velocity structure function COARE3.5
def psiu_40(zol):
""" Computes velocity structure function C35
Parameters
----------
zol : float
height over MO length
Returns
-------
psi : float
"""
dzet = np.where(0.35*zet > 50, 50, 0.35*zet) # stable
dzol = np.where(0.35*zol > 50, 50, 0.35*zol) # stable
a, b, c, d = 1, 3/4, 5, 0.35
psi = -(a*zet+b*(zet-c/d)*np.exp(-dzet)+b*c/d)
k = np.where(zet < 0) # unstable
x = (1-18*zet[k])**0.25
psi = -(a*zol+b*(zol-c/d)*np.exp(-dzol)+b*c/d)
k = np.where(zol < 0) # unstable
x = (1-18*zol[k])**0.25
psik = 2*np.log((1+x)/2)+np.log((1+x**2)/2)-2*np.arctan(x)+2*np.arctan(1)
x = (1-10*zet[k])**0.3333
x = (1-10*zol[k])**0.3333
psic = (1.5*np.log((1+x+x**2)/3)-np.sqrt(3)*np.arctan((1+2*x)/np.sqrt(3)) +
4*np.arctan(1)/np.sqrt(3))
f = zet[k]**2/(1+zet[k]**2)
f = zol[k]**2/(1+zol[k]**2)
psi[k] = (1-f)*psik+f*psic
return psi
# ---------------------------------------------------------------------
def get_skin(sst, qsea, rho, jcool, Rl, Rs, Rnl, cp, lv, usr, tsr, qsr, lat):
def get_skin(sst, qsea, rho, Rl, Rs, Rnl, cp, lv, usr, tsr, qsr, lat):
""" Computes cool skin
Parameters
----------
sst : float
sea surface temperature ($^\circ$\,C)
qsea : float
specific humidity over sea (g/kg)
rho : float
density of air (kg/m^3)
Rl : float
downward longwave radiation (W/m^2)
Rs : float
downward shortwave radiation (W/m^2)
cp : float
specific heat of air at constant pressure
lv : float
latent heat of vaporization
usr : float
friction velocity
tsr : float
star temperature
qsr : float
star humidity
lat : float
latitude
Returns
-------
dter : float
dqer : float
"""
# coded following Saunders (1967) with lambda = 6
g = gc(lat, None)
if (np.nanmin(sst) > 200): # if Ta in Kelvin convert to Celsius
......@@ -397,6 +785,27 @@ def get_skin(sst, qsea, rho, jcool, Rl, Rs, Rnl, cp, lv, usr, tsr, qsr, lat):
def get_gust(beta, Ta, usr, tsrv, zi, lat):
""" Computes gustiness
Parameters
----------
beta : float
constant
Ta : float
air temperature (K)
usr : float
friction velocity (m/s)
tsrv : float
star virtual temperature of air (K)
zi : int
scale height of the boundary layer depth (m)
lat : float
latitude
Returns
-------
ug : float
"""
if (np.max(Ta) < 200): # convert to K if in Celsius
Ta = Ta+273.16
if np.isnan(zi):
......@@ -410,6 +819,17 @@ def get_gust(beta, Ta, usr, tsrv, zi, lat):
def get_heights(h):
""" Reads input heights for velocity, temperature and humidity
Parameters
----------
h : float
input heights (m)
Returns
-------
hh : array
"""
hh = np.zeros(3)
if (type(h) == float or type(h) == int):
hh[0], hh[1], hh[2] = h, h, h
......@@ -422,10 +842,17 @@ def get_heights(h):
def svp_calc(T):
"""
calculates saturation vapour pressure
T is in Kelvin
svp in mb, pure water
""" Calculates saturation vapour pressure
Parameters
----------
T : float
temperature (K)
Returns
-------
svp : float
in mb, pure water
"""
if (np.nanmin(T) < 200): # if T in Celsius convert to Kelvin
T = T+273.16
......@@ -435,10 +862,19 @@ def svp_calc(T):
def qsea_calc(sst, pres):
"""
sst in Kelvin
pres in mb
qsea in kg/kg
""" Computes specific humidity of the sea surface air
Parameters
----------
sst : float
sea surface temperature (K)
pres : float
pressure (mb)
Returns
-------
qsea : float
(kg/kg)
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+273.16
......@@ -451,11 +887,20 @@ def qsea_calc(sst, pres):
def q_calc(Ta, rh, pres):
"""
rh in %
air in K, if not it will be converted to K
pres in mb
qair in kg/kg, as in Haltiner and Martin p.24
""" Computes specific humidity following Haltiner and Martin p.24
Parameters
----------
Ta : float
air temperature (K)
rh : float
relative humidity (%)
pres : float
air pressure (mb)
Returns
-------
qair : float, (kg/kg)
"""
if (np.nanmin(Ta) < 200): # if sst in Celsius convert to Kelvin
Ta = Ta+273.15
......@@ -466,9 +911,18 @@ def q_calc(Ta, rh, pres):
def bucksat(T, P):
"""
computes saturation vapor pressure [mb] as in COARE3.5
given T [degC] and P [mb]
""" Computes saturation vapor pressure (mb) as in C35
Parameters
----------
T : float
temperature ($^\\circ$\\,C)
P : float
pressure (mb)
Returns
-------
exx : float
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
......@@ -479,9 +933,18 @@ def bucksat(T, P):
def qsat26sea(T, P):
"""
computes surface saturation specific humidity [g/kg] as in COARE3.5
given T [degC] and P [mb]
""" Computes surface saturation specific humidity (g/kg) as in C35
Parameters
----------
T : float
temperature ($^\\circ$\\,C)
P : float
pressure (mb)
Returns
-------
qs : float
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
......@@ -494,9 +957,19 @@ def qsat26sea(T, P):
def qsat26air(T, P, rh):
"""
computes saturation specific humidity [g/kg] as in COARE3.5
given T [degC] and P [mb]
""" Computes saturation specific humidity (g/kg) as in C35
Parameters
----------
T : float
temperature ($^\circ$\,C)
P : float
pressure (mb)
Returns
-------
q : float
em : float
"""
T = np.asarray(T)
if (np.nanmin(T) > 200): # if Ta in Kelvin convert to Celsius
......@@ -509,13 +982,19 @@ def qsat26air(T, P, rh):
def gc(lat, lon=None):
"""
computes gravity relative to latitude
inputs:
lat : latitudes in deg
lon : longitudes (optional)
output:
gc: gravity constant
""" Computes gravity relative to latitude
Parameters
----------
lat : float
latitude ($^\circ$)
lon : float
longitude ($^\circ$, optional)
Returns
-------
gc : float
gravity constant (m/s^2)
"""
gamma = 9.7803267715
c1 = 0.0052790414
......@@ -535,13 +1014,18 @@ def gc(lat, lon=None):
def visc_air(Ta):
"""
Computes the kinematic viscosity of dry air as a function of air temp.
""" Computes the kinematic viscosity of dry air as a function of air temp.
following Andreas (1989), CRREL Report 89-11.
input:
Ta : air temperature [Celsius]
output
visa : kinematic viscosity [m^2/s]
Parameters
----------
Ta : float
air temperature ($^\circ$\,C)
Returns
-------
visa : float
kinematic viscosity (m^2/s)
"""
Ta = np.asarray(Ta)
if (np.nanmin(Ta) > 200): # if Ta in Kelvin convert to Celsius
......
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