diff --git a/flux_calc_v3.py b/flux_calc_v3.py
index dcc963a6eb08aca3f58d760992ff1ac186534cbd..aff60cc450ca66131780da391160b16ab3956fa2 100755
--- a/flux_calc_v3.py
+++ b/flux_calc_v3.py
@@ -1,73 +1,133 @@
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
diff --git a/flux_subs.py b/flux_subs.py
index adc074d198d3b319ac484458aa737e628575e9e5..94dc9f543cd1f0c56d10b4d7b45b236a1898e2df 100755
--- a/flux_subs.py
+++ b/flux_subs.py
@@ -1,13 +1,41 @@
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