Skip to content

GitLab

Commit 7aebd6ef authored by sbiri's avatar sbiri
Browse files

initial commit

parents
Hide whitespace changes
Inline Side-by-side
Showing with 807 additions and 0 deletions
flux_calc_v3.py 0 → 100755
View file @ 7aebd6ef
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
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
"""
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
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
# 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.isnan(wcp[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
if waveage and seastate:
print('Using seastate dependent parameterization')
logging.info('Using seastate dependent parameterization')
if waveage and ~seastate:
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))
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))
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))
# first guesses
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)
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)
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
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
monob = hh_in[0]/zetu
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))
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
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
monob[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
it = 0
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])
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)
if (np.all(np.isnan(cdn))):
break #sys.exit("cdn cannot be nan")
logging.debug('%s at iteration %s cdn negative',method,it)
zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cdn[ind]))
psim[ind] = psim_calc(hh_in[0]/monob[ind],method)
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)
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)
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
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
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)
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)
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)
cq = -usr*qsr/(dq-dqer*jcool)/wind
cqn = 1000*kappa**2/np.log(10/zo)/np.log(10/zoq)
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
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)
else:
rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht)
sensible = -1*tsr*usr*cp*rho
latent = -1*qsr*usr*lv*rho
tau = 1*rho*usr**2
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))
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
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 0 → 100755
View file @ 7aebd6ef
import numpy as np
import math
""" Conversion factor for [:math:`^\\circ` C] to [:math:`^\\circ` K] """
CtoK = 273.16 # 273.15
""" von Karman's constant """
kappa = 0.4 # NOTE: 0.41
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] = charnS
ac[:,2] = charnW
return ac
def cd_C35(u10n,wind,usr,charn,monob,Ta,hh_in,lat):
g=gc(lat,None)
zo=charn*usr**2/g+0.11*visc_air(Ta)/usr # surface roughness
rr=zo*usr/visc_air(Ta)
zoq=np.where(5.8e-5/rr**0.72>1.6e-4,1.6e-4,5.8e-5/rr**0.72) # These thermal roughness lengths give Stanton and
zot=zoq # Dalton numbers that closely approximate COARE 3.0
cdhf=kappa/(np.log(hh_in[0]/zo)-psiu_26(hh_in[0]/monob))
cthf=kappa/(np.log(hh_in[1]/zot)-psit_26(hh_in[1]/monob))
cqhf=kappa/(np.log(hh_in[2]/zoq)-psit_26(hh_in[2]/monob))
return zo,cdhf, cthf, cqhf
def cdn_calc(u10n,Ta,Tp,method="Smith80"):
if (method == "Smith80"):
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 method == "COARE4.0"):
cdn = cdn_from_roughness(u10n,Ta,None, method)
elif (method == "HEXOS"):
cdn = (0.5+0.091*u10n)*0.001 #Smith et al. 1991 #(0.27 + 0.116*u10n)*0.001 Smith et al. 1992
elif (method == "HEXOSwave"):
cdn=cdn_from_roughness(u10n,Ta,Tp, method)
elif (method == "YT96"):
cdn = np.where((u10n<6) & (u10n>=3),(0.29 + 3.1/u10n + 7.7/u10n**2)*0.001,np.where((u10n<=26) & (u10n>=6),(0.60 + 0.070*u10n)*0.001,(0.61+0.567/u10n)*0.001)) # for u<3 same as Smith80
elif (method == "LY04"):
cdn = np.where(u10n>=0.5,(0.142 + (2.7/u10n) + (u10n/13.09))*0.001,np.nan)
else:
print("unknown method cdn: "+method)
return cdn
#---------------------------------------------------------------------
def cdn_from_roughness(u10n,Ta,Tp,method="Smith88"):
g,tol = 9.812, 0.000001
cdn,ustar = np.zeros(np.asarray(u10n).shape),np.zeros(np.asarray(u10n).shape)
cdnn = (0.61 + 0.063 * u10n) * 0.001
zo,zc,zs=np.zeros(np.asarray(u10n).shape),np.zeros(np.asarray(u10n).shape),np.zeros(np.asarray(u10n).shape)
for it in range(5):
cdn = np.copy(cdnn)
ustar = np.sqrt(cdn*u10n**2)
if (method == "Smith88"):
zc = 0.011*ustar**2/g #.....Charnock roughness length (equn 4 in Smith 88)
zs = 0.11*visc_air(Ta)/ustar #.....smooth surface roughness length (equn 6 in Smith 88)
zo = zc + zs #.....equns 7 & 8 in Smith 88 to calculate new CDN
elif (method == "COARE3.0"):
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)/ ustar
zo = zc*ustar*ustar/g + zs
elif (method == "HEXOSwave"):
if (np.all(Tp==None) or np.nansum(Tp)==0):
Tp = 0.729*u10n # Taylor and Yelland 2001
cp_wave = g*Tp/2/np.pi # use input wave period
zo = 0.48*ustar**3/g/cp_wave # Smith et al. 1992
else:
print("unknown method for cdn_from_roughness "+method)
cdnn = (kappa/np.log(10/zo))**2
cdn=np.where(np.abs(cdnn-cdn) < tol, cdnn,np.nan)
return cdnn
#---------------------------------------------------------------------
def ctcqn_calc(zol, cdn, u10n, zo,Ta, method="Smith80"):
l = np.shape(u10n)
if (method == "Smith80" or method == "Smith88" or method == "YT96"):
cqn = np.ones(l)*1.20*0.001 # from Smith88, no value given by S80 or YT80
ctn = np.ones(l)*1.00*0.001
elif (method == "LP82"):
cqn = np.where((zol <= 0) & (u10n>4) & (u10n<14),1.15*0.001,np.nan)
ctn = np.where((zol <= 0) & (u10n>4) & (u10n<25), 1.13*0.001, 0.66*0.001)
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"):
ustar = (cdn * u10n**2)**0.5
rr = zo*ustar/visc_air(Ta)
zoq = 5.5e-5/rr**0.6
zoq[zoq>1.15e-4] = 1.15e-4
zot = zoq
cqn = kappa**2/np.log(10/zo)/np.log(10/zoq)
ctn = kappa**2/np.log(10/zo)/np.log(10/zot)
elif (method == "LY04"):
cqn = 34.6*0.001*cdn**0.5
ctn = np.where(zol <= 0, 32.7*0.001*cdn**0.5, 18*0.001*cdn**0.5)
else:
print("unknown method ctcqn: "+method)
return ctn, cqn
#---------------------------------------------------------------------
def cd_calc(cdn, height, ref_ht, psim):
cd = cdn*np.power(1+np.sqrt(cdn)*np.power(kappa,-1)*(np.log(height/ref_ht)-psim),-2)
return cd
#---------------------------------------------------------------------
def ctcq_calc(cdn, cd, ctn, cqn, h_t, h_q, ref_ht, psit, psiq):
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"):
coeffs = get_stabco(method)
alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (method == "COARE3.0" or method == "COARE4.0"):
psim = np.where(zol<0,psim_conv_coare3(zol,alpha,beta,gamma),psim_stab_coare3(zol,alpha,beta,gamma))
else:
psim = np.where(zol<0,psim_conv(zol,alpha,beta,gamma),psim_stab(zol,alpha,beta,gamma))
return psim
#---------------------------------------------------------------------
def psit_calc(zol, method="Smith80"):
coeffs = get_stabco(method)
alpha, beta, gamma = coeffs[0], coeffs[1], coeffs[2]
if (method == "COARE3.0" or method == "COARE4.0"):
psit = np.where(zol<0,psi_conv_coare3(zol,alpha,beta,gamma),psi_stab_coare3(zol,alpha,beta,gamma))
else:
psit = np.where(zol<0,psi_conv(zol,alpha,beta,gamma),psi_stab(zol,alpha,beta,gamma))
return psit
#---------------------------------------------------------------------
def get_stabco(method="Smith80"):
if (method == "Smith80" or method == "Smith88" or method == "LY04"):
# Smith 1980, from Dyer (1974)
alpha, beta, gamma = 16, 0.25, 5
elif (method == "LP82"):
alpha, beta, gamma = 16, 0.25, 7
elif (method == "HEXOS" or method == "HEXOSwave"):
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"):
# use separate subroutine
alpha, beta, gamma = 15, 1/3, 5 # not sure about gamma=34.15
#alpha <- NA #beta <- NA
else:
print("unknown method stabco: "+method)
coeffs = np.zeros(3)
coeffs[0] = alpha
coeffs[1] = beta
coeffs[2] = gamma
return coeffs
#---------------------------------------------------------------------
#======================================================================
def psi_conv_coare3(zol,alpha,beta,gamma):
x = (1-alpha*zol)**0.5 # Kansas unstable
psik = 2*np.log((1+x)/2.)
y = (1-34.15*zol)**beta
psic = 1.5*np.log((1+y+y*y)/3.)-(3)**0.5*np.arctan((1+2*y)/
(3)**0.5)+4*np.arctan(1)/(3)**0.5
f = zol*zol/(1.+zol*zol)
psit = (1-f)*psik+f*psic
return psit
#======================================================================
def psi_stab_coare3(zol,alpha,beta,gamma): #Stable
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
#======================================================================
def psi_conv(zol,alpha,beta,gamma):
xtmp = (1 - alpha*zol)**beta
psit = 2*np.log((1+xtmp**2)*0.5)
return psit
#======================================================================
def psi_stab(zol,alpha,beta,gamma):
psit = -gamma*zol
return psit
#======================================================================
def psim_conv_coare3(zol,alpha,beta,gamma):
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
psic = 1.5*np.log((1+y+y*y)/3.)-np.sqrt(3)*np.arctan((1+2*y)/np.sqrt(3))+4.*np.arctan(1)/np.sqrt(3)
f = zol*zol/(1+zol*zol)
psim = (1-f)*psik+f*psic
return psim
#======================================================================
def psim_stab_coare3(zol,alpha,beta,gamma):
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
#======================================================================
def psim_conv(zol,alpha,beta,gamma):
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
return psim
#======================================================================
def psim_stab(zol,alpha,beta,gamma):
psim = -gamma*zol
return psim
#======================================================================
#---------------------------------------------------------------------
def get_skin(sst,qsea,rho,jcool,Rl,Rs,Rnl,cp,lv,usr,tsr,qsr,lat):
# coded following Saunders (1967) with lambda = 6
g=gc(lat,None)
if ( np.nanmin(sst) > 200 ): # if Ta in Kelvin convert to Celsius
sst = sst-273.16
#************ cool skin constants *******
rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6 # density of water, specific heat capacity of water, water viscosity, thermal conductivity of water
Al = 2.1e-5*(sst+3.2)**0.79
be = 0.026
bigc = 16*g*cpw*(rhow*visw)**3/(tcw*tcw*rho*rho)
wetc = 0.622*lv*qsea/(287.1*(sst+273.16)**2)
Rns = 0.945*Rs # albedo correction
hsb=-rho*cp*usr*tsr
hlb=-rho*lv*usr*qsr
qout=Rnl+hsb+hlb
tkt = 0.001*np.ones(np.shape(sst))
dels=Rns*(0.065+11*tkt-6.6e-5/tkt*(1-np.exp(-tkt/8.0e-4)))
qcol=qout-dels
alq=Al*qcol+be*hlb*cpw/lv
xlamx=np.where(alq>0,6/(1+(bigc*alq/usr**4)**0.75)**0.333,6)
tkt=xlamx*visw/(np.sqrt(rho/rhow)*usr) #np.nanmin(0.01, xlamx*visw/(np.sqrt(rhoa/rhow)*usr))
tkt=np.where(alq>0,np.where(tkt > 0.01, 0.01,tkt),tkt)
dter=qcol*tkt/tcw
dqer=wetc*dter
return dter, dqer
#---------------------------------------------------------------------
def get_gust(Ta,usr,tsrv,zi,lat):
if (np.max(Ta)<200): # convert to K if in Celsius
Ta=Ta+273.16
if np.isnan(zi):
zi = 600
g=gc(lat,None)
Bf=-g/Ta*usr*tsrv
ug=np.ones(np.shape(Ta))*0.2
ug=np.where(Bf>0,1.2*(Bf*zi)**0.333,0.2)
return ug
#---------------------------------------------------------------------
def get_heights(h):
hh=np.zeros(3)
if type(h) == float or type(h) == int:
hh[0], hh[1], hh[2] = h, h, h
elif len(h)==2:
hh[0], hh[1], hh[2] = h[0], h[1], h[1]
else:
hh[0], hh[1], hh[2] = h[0], h[1], h[2]
return hh
#---------------------------------------------------------------------
def svp_calc(T):
# t is in Kelvin
# svp in mb, pure water
if ( np.nanmin(T) < 200 ): # if T in Celsius convert to Kelvin
T = T+273.16
svp = np.where(np.isnan(T), np.nan,2.1718e08*np.exp(-4157/(T-33.91-0.16)))
return svp
#---------------------------------------------------------------------
def qsea_calc(sst,pres):
# sst in Kelvin
# pres in mb
# qsea in kg/kg
if ( np.nanmin(sst) < 200 ): # if sst in Celsius convert to Kelvin
sst = sst+273.16
ed = svp_calc(sst)
e = 0.98*ed
qsea = (0.622*e)/(pres-0.378*e)
qsea = np.where(~np.isnan(sst+pres),qsea,np.nan)
return qsea
#---------------------------------------------------------------------
def rh_calc(Ta,qair,pres):
if ( np.nanmin(Ta) < 200 ): # if sst in Celsius convert to Kelvin
Ta = Ta+273.16
e = np.where(np.isnan(Ta+qair+pres),np.nan,(qair*pres)/(0.62197+qair*0.378))
ed = np.where(np.isnan(e),np.nan,svp_calc(Ta))
rh = np.where(np.isnan(ed),np.nan,e/ed*100)
return rh
#---------------------------------------------------------------------
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
if ( np.nanmin(Ta) < 200 ): # if sst in Celsius convert to Kelvin
Ta = Ta+273.15
e = np.where(np.isnan(Ta+rh+pres),np.nan,svp_calc(Ta)*rh*0.01)
qair = np.where(np.isnan(e),np.nan,((0.62197*e)/(pres-0.378*e))) # Haltiner and Martin p.24
return qair
#---------------------------------------------------------------------
def gc(lat,lon=None):
"""
computes gravity relative to latitude
inputs:
lat : latitudes in deg
lon : longitudes (optional)
output:
gc: gravity constant
"""
gamma = 9.7803267715
c1 = 0.0052790414
c2 = 0.0000232718
c3 = 0.0000001262
c4 = 0.0000000007
if lon is not None:
lon_m,lat_m = np.meshgrid(lon,lat)
else:
lat_m = lat
phi = lat_m*np.pi/180.
xx = np.sin(phi)
gc = gamma*(1+c1*np.power(xx,2)+c2*np.power(xx,4)+c3*np.power(xx,6)+c4*np.power(xx,8))
return gc
#---------------------------------------------------------------------
def PtoDepth(p,Lat):
"""
computes depth in m from pressure in db following
Saunders and Fofonoff (1976). Deep sea res.
"""
x=math.sin(math.radians(Lat))
x=x*x
# gravity variation with latitude (Anon, 1970)
gr = 9.780318 * (1 + (5.2788E-3 + 2.36E-5 * x) * x) + 1.092E-6 * p
d = (((-1.82E-15 * p + 2.279E-10) * p - 2.2512E-5) * p + 9.72659) * p;
depth = d/gr
return depth
#---------------------------------------------------------------------
def visc_air(Ta):
"""
Computes the kinematic viscosity of dry air as a function of air temperature
following Andreas (1989), CRREL Report 89-11.
input:
Ta : air temperature [Celsius]
output
visa : kinematic viscosity [m^2/s]
"""
Ta = np.asarray(Ta)
if ( np.nanmin(Ta) > 200 ): # if Ta in Kelvin convert to Celsius
Ta = Ta-273.16
visa = 1.326e-5 * (1 + 6.542e-3*Ta + 8.301e-6*Ta**2 - 4.84e-9*Ta**3)
return visa
######
# functions from coare35vn.mat
######
#------------------------------------------------------------------------------
def psit_26(zet):
"""
computes temperature structure function
"""
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
psik=2*np.log((1+x)/2)
x=(1-34.15*zet[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)
psi[k]=(1-f)*psik+f*psic
return psi
#------------------------------------------------------------------------------
def psiu_26(zet):
"""
computes velocity structure function
"""
dzet=np.where(0.35*zet > 50, 50, 0.35*zet) # 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
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
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)
psi[k]=(1-f)*psik+f*psic
return psi
#------------------------------------------------------------------------------
def psiu_40(zet):
"""
computes velocity structure function
"""
dzet=np.where(0.35*zet > 50, 50, 0.35*zet) # 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
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
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)
psi[k]=(1-f)*psik+f*psic
return psi
#------------------------------------------------------------------------------
def bucksat(T,P):
"""
computes saturation vapor pressure [mb]
given T [degC] and P [mb]
"""
T = np.asarray(T)
if ( np.nanmin(T) > 200 ): # if Ta in Kelvin convert to Celsius
T = T-CtoK
exx=6.1121*np.exp(17.502*T/(T+240.97))*(1.0007+3.46e-6*P)
return exx
#------------------------------------------------------------------------------
def qsat26sea(T,P):
"""
computes surface saturation specific humidity [g/kg]
given T [degC] and P [mb]
"""
T = np.asarray(T)
if ( np.nanmin(T) > 200 ): # if Ta in Kelvin convert to Celsius
T = T-CtoK
ex=bucksat(T,P)
es=0.98*ex # reduction at sea surface
qs=622*es/(P-0.378*es)
return qs
#------------------------------------------------------------------------------
def qsat26air(T,P,rh):
"""
computes saturation specific humidity [g/kg]
given T [degC] and P [mb]
"""
T = np.asarray(T)
if ( np.nanmin(T) > 200 ): # if Ta in Kelvin convert to Celsius
T = T-CtoK
es=bucksat(T,P)
em=0.01*rh*es
q=622*em/(P-0.378*em)
return q,em
#------------------------------------------------------------------------------
def RHcalc(T,P,Q):
"""
computes relative humidity given T,P, & Q
"""
es=6.1121*np.exp(17.502*T/(T+240.97))*(1.0007+3.46e-6*P)
em=Q*P/(0.378*Q+0.622)
RHrf=100*em/es
return RHrf
\ No newline at end of file
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