Skip to content

GitLab

Commit e3c3222c authored by Richard Cornes's avatar Richard Cornes
Browse files

Changed to class structure

parent 0c4ac7f6
Hide whitespace changes
Inline Side-by-side
Showing with 543 additions and 471 deletions
Code/AirSeaFluxCode.py
View file @ e3c3222c
import numpy as np
import pandas as pd
import logging
from get_init import get_init
from hum_subs import (get_hum, gamma)
from util_subs import (kappa, CtoK, get_heights, gc)
from flux_subs import (cs_C35, cs_Beljaars, cs_ecmwf, wl_ecmwf,
get_gust, get_L, get_strs, rough, psim_calc,
psit_calc, cdn_calc, cd_calc, ctcq_calc, ctcqn_calc)
from util_subs import *
from flux_subs import *
class S80:
def get_heights(self, hin, hout=10):
self.hout = hout
self.hin = hin
self.h_in = get_heights(hin, len(self.spd))
self.h_out = get_heights(self.hout,1)
def get_specHumidity(self,qmeth="Buck2"):
self.qair, self.qsea = get_hum(self.hum, self.T, self.SST, self.P, qmeth)
if (np.all(np.isnan(self.qsea)) or np.all(np.isnan(self.qair))):
raise ValueError("qsea and qair cannot be nan")
self.dq = self.qair - self.qsea
# Set lapse rate and Potential Temperature (now we have humdity)
self._get_potentialT()
self._get_lapse()
def _get_potentialT(self):
self.cp = 1004.67*(1+0.00084*self.qsea)
self.th = np.where(self.T < 200, (np.copy(self.T)+CtoK) *
np.power(1000/self.P, 287.1/self.cp),
np.copy(self.T)*np.power(1000/self.P, 287.1/self.cp)) # potential T
def _get_lapse(self):
self.tlapse = gamma("dry", self.SST, self.T, self.qair/1000, self.cp)
self.Ta = np.where(self.T < 200, np.copy(self.T)+CtoK+self.tlapse*self.h_in[1],
np.copy(self.T)+self.tlapse*self.h_in[1]) # convert to Kelvin if needed
self.dt = self.Ta - self.SST
def _fix_coolskin_warmlayer(self, wl, cskin, skin, Rl, Rs):
assert wl in [0,1], "wl not valid"
assert cskin in [0,1], "cskin not valid"
assert skin in ["C35", "ecmwf" or "Beljaars"], "Skin value not valid"
if ((cskin == 1 or wl == 1) and (np.all(Rl == None) or np.all(np.isnan(Rl)))
and ((np.all(Rs == None) or np.all(np.isnan(Rs))))):
print("Cool skin/warm layer is switched ON; Radiation input should not be empty")
raise
self.wl = wl
self.cskin = cskin
self.Rs = np.ones(self.spd.shape)*np.nan if Rs is None else Rs
self.Rl = np.ones(self.spd.shape)*np.nan if Rl is None else Rl
def set_coolskin_warmlayer(self, wl=0, cskin=0, skin="C35", Rl=None, Rs=None):
wl = 0 if wl is None else wl
self._fix_coolskin_warmlayer(wl, cskin, skin, Rl, Rs)
def _first_guess(self):
# reference height1
self.ref_ht = 10
# first guesses
self.t10n, self.q10n = np.copy(self.Ta), np.copy(self.qair)
self.tv10n = self.t10n*(1+0.6077*self.q10n)
# Zeng et al. 1998
tv = self.th*(1+0.6077*self.qair) # virtual potential T
self.dtv = self.dt*(1+0.6077*self.qair)+0.6077*self.th*self.dq
# Rb eq. 11 Grachev & Fairall 1997
Rb = self.g*10*(self.dtv)/(np.where(self.T < 200, np.copy(self.T)+CtoK, np.copy(self.T)) * np.power(self.wind, 2))
#Rb = self.g*10*(self.dtv)/(np.copy(self.T) * np.power(self.wind, 2))
self.monob = 1/Rb # eq. 12 Grachev & Fairall 1997
# ------------
self.rho = self.P*100/(287.1*self.tv10n)
self.lv = (2.501-0.00237*(self.SST-CtoK))*1e6 # J/kg
self.dter = np.full(self.T.shape, -0.3)*self.msk
self.tkt = np.full(self.T.shape, 0.001)*self.msk
self.dqer = self.dter*0.622*self.lv*self.qsea/(287.1*np.power(self.SST, 2))
self.Rnl = 0.97*(self.Rl-5.67e-8*np.power(self.SST-0.3*self.cskin, 4))
self.Qs = 0.945*self.Rs
self.dtwl = np.full(self.T.shape,0.3)*self.msk
self.skt = np.copy(self.SST)
# Apply the gustiness adjustment if defined for this class
try:
self._adjust_gust()
except AttributeError:
pass
self.u10n = self.wind*np.log(10/1e-4)/np.log(self.hin[0]/1e-4)
self.usr = 0.035*self.u10n
self.cd10n = cdn_calc(self.u10n, self.usr, self.Ta, self.lat, self.meth)
self.tsr = np.zeros(self.arr_shp)*self.msk
self.tsrv = np.copy(self.tsr)
self.qsr = np.copy(self.tsr)
self.psim = np.copy(self.tsr)
self.psit = np.copy(self.tsr)
self.psiq = np.copy(self.tsr)
self.cd = cd_calc(self.cd10n, self.h_in[0], self.ref_ht, self.psim)
self.usr = np.sqrt(self.cd*np.power(self.wind, 2))
self.zo = np.full(self.arr_shp,1e-4)*self.msk
self.zot = np.full(self.arr_shp,1e-4)*self.msk
self.zoq = np.full(self.arr_shp,1e-4)*self.msk
self.ct10n = np.power(kappa, 2)/(np.log(self.h_in[0]/self.zo)*np.log(self.h_in[1]/self.zot))
self.cq10n = np.power(kappa, 2)/(np.log(self.h_in[0]/self.zo)*np.log(self.h_in[2]/self.zoq))
self.ct = np.power(kappa, 2)/((np.log(self.h_in[0]/self.zo)-self.psim) * (np.log(self.h_in[1]/self.zot)-self.psit))
self.cq = np.power(kappa, 2)/((np.log(self.h_in[0]/self.zo)-self.psim) * (np.log(self.h_in[2]/self.zoq)-self.psiq))
self.tsr = (self.dt-self.dter*self.cskin-self.dtwl*self.wl)*kappa/(np.log(self.h_in[1]/self.zot) - psit_calc(self.h_in[1]/self.monob, self.meth))
self.qsr = (self.dq-self.dqer*self.cskin)*kappa/(np.log(self.h_in[2]/self.zoq) - psit_calc(self.h_in[2]/self.monob, self.meth))
self.tau = np.full(self.arr_shp,0.05)*self.msk
self.sensible = np.full(self.arr_shp,-5)*self.msk
self.latent = np.full(self.arr_shp,-65)*self.msk
def _zo_calc(self, ref_ht, cd10n):
zo = ref_ht/np.exp(kappa/np.sqrt(cd10n))
return zo
def iterate(self,n=10, tol=None):
n = 5 if n < 5 else n
tol = ['all', 0.01, 0.01, 1e-05, 1e-3, 0.1, 0.1] if tol is None else tol
assert tol[0] in ['flux', 'ref', 'all'], "unknown tolerance input"
self.itera = np.full(self.arr_shp,-1)*self.msk
ind = np.where(self.spd > 0)
it = 0
# Generate the first guess values
self._first_guess()
# iteration loop
ii = True
while ii:
it += 1
if it > n: break
if (tol[0] == 'flux'):
old = np.array([np.copy(self.tau), np.copy(self.sensible), np.copy(self.latent)])
elif (tol[0] == 'ref'):
old = np.array([np.copy(self.u10n), np.copy(self.t10n), np.copy(self.q10n)])
elif (tol[0] == 'all'):
old = np.array([np.copy(self.u10n), np.copy(self.t10n), np.copy(self.q10n),
np.copy(self.tau), np.copy(self.sensible), np.copy(self.latent)])
self.cd10n[ind] = cdn_calc(self.u10n[ind], self.usr[ind], self.Ta[ind], self.lat[ind], self.meth)
if (np.all(np.isnan(self.cd10n))):
break
logging.info('break %s at iteration %s cd10n<0', meth, it)
self.zo[ind] = self.ref_ht/np.exp(kappa/np.sqrt(self.cd10n[ind]))
#self.zo[ind] = self._zo_calc(self.ref_ht, self.cd10n[ind])
self.psim[ind] = psim_calc(self.h_in[0, ind]/self.monob[ind], self.meth)
self.cd[ind] = cd_calc(self.cd10n[ind], self.h_in[0, ind], self.ref_ht, self.psim[ind])
self.ct10n[ind], self.cq10n[ind] = ctcqn_calc(self.h_in[1, ind]/self.monob[ind], self.cd10n[ind],
self.usr[ind], self.zo[ind], self.Ta[ind], self.meth)
self.zot[ind] = self.ref_ht/(np.exp(np.power(kappa, 2) / (self.ct10n[ind]*np.log(self.ref_ht/self.zo[ind]))))
self.zoq[ind] = self.ref_ht/(np.exp(np.power(kappa, 2) / (self.cq10n[ind]*np.log(self.ref_ht/self.zo[ind]))))
self.psit[ind] = psit_calc(self.h_in[1, ind]/self.monob[ind], self.meth)
self.psiq[ind] = psit_calc(self.h_in[2, ind]/self.monob[ind], self.meth)
self.ct[ind], self.cq[ind] = ctcq_calc(self.cd10n[ind], self.cd[ind], self.ct10n[ind],self.cq10n[ind], self.h_in[:, ind],
[self.ref_ht, self.ref_ht, self.ref_ht], self.psit[ind], self.psiq[ind])
if (self.meth == "NCAR"):
self.cd = np.maximum(np.copy(self.cd), 1e-4)
self.ct = np.maximum(np.copy(self.ct), 1e-4)
self.cq = np.maximum(np.copy(self.cq), 1e-4)
self.zo = np.minimum(np.copy(self.zo), 0.0025)
self.usr[ind],self.tsr[ind], self.qsr[ind] = get_strs(self.h_in[:, ind], self.monob[ind],
self.wind[ind], self.zo[ind], self.zot[ind],
self.zoq[ind], self.dt[ind], self.dq[ind],
self.dter[ind], self.dqer[ind],
self.dtwl[ind], self.ct[ind], self.cq[ind],
self.cskin, self.wl, self.meth)
self.dter[ind] = np.zeros(self.SST[ind].shape)
self.dqer[ind] = np.zeros(self.SST[ind].shape)
self.tkt[ind] = np.full(self.T[ind].shape,0.001)
# Logging output
log_vars = {"dter":2,"dqer":7,"tkt":2,"Rnl":2,"usr":3,"tsr":4,"qsr":7}
log_vars = [np.round(np.nanmedian(getattr(self,V)),R) for V,R in log_vars.items()]
log_vars.insert(0,self.meth)
logging.info('method {} | dter = {} | dqer = {} | tkt = {} | Rnl = {} | usr = {} | tsr = {} | qsr = {}'.format(*log_vars))
self.Rnl[ind] = 0.97*(self.Rl[ind]-5.67e-8*np.power(self.SST[ind] + self.dter[ind]*self.cskin, 4))
self.t10n[ind] = (self.Ta[ind] - self.tsr[ind]/kappa*(np.log(self.h_in[1, ind]/self.ref_ht)-self.psit[ind]))
self.q10n[ind] = (self.qair[ind] - self.qsr[ind]/kappa*(np.log(self.h_in[2, ind]/self.ref_ht)-self.psiq[ind]))
self.tv10n[ind] = self.t10n[ind]*(1+0.6077*self.q10n[ind])
self.tsrv[ind], self.monob[ind], self.Rb[ind] = get_L(self.L, self.lat[ind], self.usr[ind], self.tsr[ind], self.qsr[ind], self.h_in[:, ind], self.Ta[ind],
(self.SST[ind]+self.dter[ind]*self.cskin + self.dtwl[ind]*self.wl), self.qair[ind], self.qsea[ind], self.wind[ind],
np.copy(self.monob[ind]), self.zo[ind], self.zot[ind], self.psim[ind], self.meth)
self.psim[ind] = psim_calc(self.h_in[0, ind]/self.monob[ind], self.meth)
self.psit[ind] = psit_calc(self.h_in[1, ind]/self.monob[ind], self.meth)
self.psiq[ind] = psit_calc(self.h_in[2, ind]/self.monob[ind], self.meth)
if (self.gust[0] == 1):
self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) + np.power(get_gust(self.gust[1], self.Ta[ind], self.usr[ind],
self.tsrv[ind], self.gust[2], self.lat[ind]), 2))
elif (self.gust[0] == 0):
self.wind[ind] = np.copy(self.spd[ind])
self.u10n[ind] = self.wind[ind]-self.usr[ind]/kappa*(np.log(self.h_in[0, ind]/10) - self.psim[ind])
# make sure you allow small negative values convergence
if (it < 4): self.u10n = np.where(self.u10n < 0, 0.5, self.u10n)
self.utmp = np.copy(self.u10n)
self.utmp = np.where(self.utmp < 0, np.nan, self.utmp)
self.itera[ind] = np.ones(1)*it
self.tau = self.rho*np.power(self.usr, 2)*(self.spd/self.wind)
self.sensible = self.rho*self.cp*self.usr*self.tsr
self.latent = self.rho*self.lv*self.usr*self.qsr
if (tol[0] == 'flux'):
new = np.array([np.copy(self.tau), np.copy(self.sensible), np.copy(self.latent)])
elif (tol[0] == 'ref'):
new = np.array([np.copy(self.utmp), np.copy(self.t10n), np.copy(self.q10n)])
elif (tol[0] == 'all'):
new = np.array([np.copy(self.utmp), np.copy(self.t10n), np.copy(self.q10n),
np.copy(self.tau), np.copy(self.sensible), np.copy(self.latent)])
if (it > 2): # force at least two iterations
d = np.abs(new-old)
if (tol[0] == 'flux'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + (d[2, :] > tol[3]))
elif (tol[0] == 'ref'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + (d[2, :] > tol[3]))
elif (tol[0] == 'all'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) + (d[2, :] > tol[3])+(d[3, :] > tol[4]) +
(d[4, :] > tol[5])+(d[5, :] > tol[6]))
ii = False if (ind[0].size == 0) else True
# End of iteration loop
self.itera[ind] = -1
self.itera = np.where(self.itera > n, -1, self.itera)
self.ind = ind
logging.info('method %s | # of iterations:%s', self.meth, it)
logging.info('method %s | # of points that did not converge :%s \n', self.meth, self.ind[0].size)
def _get_humidity(self):
"RH only used for flagging purposes"
if ((self.hum[0] == 'rh') or (self.hum[0] == 'no')):
self.rh = self.hum[1]
elif (self.hum[0] == 'Td'):
Td = self.hum[1] # dew point temperature (K)
Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
T = np.where(T < 200, np.copy(self.T)+CtoK, np.copy(self.T))
#T = np.copy(self.T)
esd = 611.21*np.exp(17.502*((Td-CtoK)/(Td-32.19)))
es = 611.21*np.exp(17.502*((T-CtoK)/(T-32.19)))
self.rh = 100*esd/es
def _flag(self,out=0):
"Set the general flags"
flag = np.full(self.arr_shp, "n",dtype="object")
if (self.hum[0] == 'no'):
self.flag = np.where(np.isnan(self.spd+self.T+self.SST+self.P), "m", flag)
else:
flag = np.where(np.isnan(self.spd+self.T+self.SST+self.hum[1]+self.P), "m", flag)
flag = np.where(self.rh > 100, "r", flag)
flag = np.where((self.u10n < 0) & (flag == "n"), "u",
np.where((self.u10n < 0) &
(np.char.find(flag.astype(str), 'u') == -1),
flag+[","]+["u"], flag))
flag = np.where((self.q10n < 0) & (flag == "n"), "q",
np.where((self.q10n < 0) & (flag != "n"), flag+[","]+["q"],
flag))
flag = np.where(((self.Rb < -0.5) | (self.Rb > 0.2) | ((self.hin[0]/self.monob) > 1000)) &
(flag == "n"), "l",
np.where(((self.Rb < -0.5) | (self.Rb > 0.2) |
((self.hin[0]/self.monob) > 1000)) &
(flag != "n"), flag+[","]+["l"], flag))
if (out == 1):
flag = np.where((self.itera == -1) & (flag == "n"), "i",
np.where((self.itera == -1) &
((flag != "n") &
(np.char.find(flag.astype(str), 'm') == -1)),
flag+[","]+["i"], flag))
else:
flag = np.where((self.itera == -1) & (flag == "n"), "i",
np.where((self.itera == -1) &
((flag != "n") &
(np.char.find(flag.astype(str), 'm') == -1) &
(np.char.find(flag.astype(str), 'u') == -1)),
flag+[","]+["i"], flag))
self.flag = flag
def _class_flag(self):
"A flag specific to this class"
self.flag = np.where(((self.utmp < 6) | (self.utmp > 22)) & (self.flag == "n"), "o",
np.where(((self.utmp < 6) | (self.utmp > 22)) &
((self.flag != "n") &
(np.char.find(self.flag.astype(str), 'u') == -1) &
(np.char.find(self.flag.astype(str), 'q') == -1)),
self.flag+[","]+["o"], self.flag))
def get_output(self,out=0):
assert out in [0,1], "out must be either 0 or 1"
# calculate output parameters
rho = (0.34838*self.P)/(self.tv10n)
self.t10n = self.t10n-(273.16+self.tlapse*self.ref_ht)
# solve for zo from cd10n
zo = self.ref_ht/np.exp(kappa/np.sqrt(self.cd10n))
# adjust neutral cdn at any output height
self.cdn = np.power(kappa/np.log(self.hout/zo), 2)
self.cd = cd_calc(self.cdn, self.h_out[0], self.h_out[0], self.psim)
# solve for zot, zoq from ct10n, cq10n
zot = self.ref_ht/(np.exp(kappa**2/(self.ct10n*np.log(self.ref_ht/zo))))
zoq = self.ref_ht/(np.exp(kappa**2/(self.cq10n*np.log(self.ref_ht/zo))))
# adjust neutral ctn, cqn at any output height
self.ctn = np.power(kappa, 2)/(np.log(self.h_out[0]/zo)*np.log(self.h_out[1]/zot))
self.cqn = np.power(kappa, 2)/(np.log(self.h_out[0]/zo)*np.log(self.h_out[2]/zoq))
self.ct, self.cq = ctcq_calc(self.cdn, self.cd, self.ctn, self.cqn, self.h_out, self.h_out, self.psit, self.psiq)
self.uref = (self.spd-self.usr/kappa*(np.log(self.h_in[0]/self.h_out[0])-self.psim + psim_calc(self.h_out[0]/self.monob, self.meth)))
tref = (self.Ta-self.tsr/kappa*(np.log(self.h_in[1]/self.h_out[1])-self.psit + psit_calc(self.h_out[0]/self.monob, self.meth)))
# Changed tref to use input temperature (T)
#self.tref = (self.T-self.tsr/kappa*(np.log(self.h_in[1]/self.h_out[1])-self.psit + psit_calc(self.h_out[0]/self.monob, self.meth)))
self.tref = tref-(CtoK+self.tlapse*self.h_out[1])
self.qref = (self.qair-self.qsr/kappa*(np.log(self.h_in[2]/self.h_out[2]) - self.psit+psit_calc(self.h_out[2]/self.monob, self.meth)))
if (self.wl == 0): self.dtwl = np.zeros(self.T.shape)*self.msk # reset to zero if not used
# Get the Relative humidity
self._get_humidity()
# Do not calculate lhf if a measure of humidity is not input
if (self.hum[0] == 'no'):
self.latent = np.ones(self.SST.shape)*np.nan
self.qsr = np.ones(self.SST.shape)*np.nan
self.q10n = np.ones(self.SST.shape)*np.nan
self.qref = np.ones(self.SST.shape)*np.nan
self.qair = np.ones(self.SST.shape)*np.nan
self.rh = np.ones(self.SST.shape)*np.nan
self.wind_spd = np.sqrt(np.power(self.wind, 2)-np.power(self.spd, 2))
# Combine all output variables into a pandas array
res_vars = ("tau","sensible","latent","monob","cd","cdn","ct","ctn","cq","cqn","tsrv","tsr","qsr","usr","psim","psit",
"psiq","u10n","t10n","tv10n","q10n","zo","zot","zoq","uref","tref","qref","itera","dter","dqer","dtwl",
"qair","qsea","Rl","Rs","Rnl","wind_spd","Rb","rh","tkt","lv")
res = np.zeros((len(res_vars), len(self.spd)))
for i, value in enumerate(res_vars): res[i][:] = getattr(self, value)
if (out == 0):
res[:, self.ind] = np.nan
# set missing values where data have non acceptable values
if (self.hum[0] != 'no'):
res = np.asarray([np.where(self.q10n < 0, np.nan, res[i][:]) for i in range(41)]) # FIXME: why 41?
res = np.asarray([np.where(self.u10n < 0, np.nan, res[i][:]) for i in range(41)])
else:
warnings.warn("Warning: the output will contain values for points that have not converged and negative values (if any) for u10n/q10n")
resAll = pd.DataFrame(data=res.T, index=range(self.nlen), columns=res_vars)
# Get flags
self._flag(out=out)
# Get class specific flags
try:
self._class_flag()
except AttributeError:
pass
resAll["flag"] = self.flag
return resAll
def add_variables(self, spd, T, SST, lat=None, hum=None, P=None, L=None):
# Add the mandatory variables
assert type(spd)==type(T)==type(SST)==np.ndarray, "input type of spd, T and SST should be numpy.ndarray"
self.L="tsrv" if L is None else L
self.arr_shp = spd.shape
self.nlen = len(spd)
self.spd = spd
#self.T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T))
self.T = T
self.hum = ['no', np.full(SST.shape,80)] if hum is None else hum
self.SST = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
self.lat = np.full(self.arr_shp,45) if lat is None else lat
self.g = gc(self.lat)
self.P = np.full(n, 1013) if P is None else P
# mask to preserve missing values when initialising variables
self.msk=np.empty(SST.shape)
self.msk = np.where(np.isnan(spd+T+SST), np.nan, 1)
self.Rb = np.empty(SST.shape)*self.msk
self.dtwl = np.full(T.shape,0.3)*self.msk
# Set the wind array
if (self.gust[0] == 1):
self.wind = np.sqrt(np.power(np.copy(self.spd), 2)+np.power(0.5, 2))
elif (self.gust[0] == 0):
self.wind = np.copy(spd)
def add_gust(self,gust=None):
if np.all(gust == None):
try:
gust = self.default_gust
except AttributeError:
gust = [1,1.2,800]
elif ((np.size(gust) < 3) and (gust == 0)):
gust = [0, 0, 0]
assert np.size(gust) == 3, "gust input must be a 3x1 array"
self.gust = gust
def __init__(self):
self.meth = "S80"
class S88(S80):
def __init__(self):
self.meth = "S88"
class YT96(S80):
def _class_flag(self):
self.flag = np.where(((self.utmp < 0) | (self.utmp > 26)) & (self.flag == "n"), "o",
np.where(((self.utmp < 3) | (self.utmp > 26)) &
((self.flag != "n") &
(np.char.find(self.flag.astype(str), 'u') == -1) &
(np.char.find(self.flag.astype(str), 'q') == -1)),
self.flag+[","]+["o"], self.flag))
def __init__(self):
self.meth = "YT96"
class LP82(S80):
def _class_flag(self):
self.flag = np.where(((self.utmp < 3) | (self.utmp > 25)) & (self.flag == "n"), "o",
np.where(((self.utmp < 3) | (self.utmp > 25)) &
((self.flag != "n") &
(np.char.find(self.flag.astype(str), 'u') == -1) &
(np.char.find(self.flag.astype(str), 'q') == -1)),
self.flag+[","]+["o"], self.flag))
def __init__(self):
self.meth = "LP82"
class NCAR(S80):
def _class_flag(self):
self.flag = np.where((self.utmp < 0.5) & (self.flag == "n"), "o",
np.where((self.utmp < 0.5) &
((self.flag != "n") &
(np.char.find(self.flag.astype(str), 'u') == -1) &
(np.char.find(self.flag.astype(str), 'q') == -1)),
self.flag+[","]+["o"], self.flag))
def _zo_calc(self, ref_ht, cd10n):
"Special z0 calculation for NCAR"
zo = ref_ht/np.exp(kappa/np.sqrt(cd10n))
zo = np.minimum(np.copy(zo), 0.0025)
return zo
def __init__(self):
self.meth = "NCAR"
class UA(S80):
def _class_flag(self):
self.flag = np.where((self.utmp > 18) & (self.flag == "n"), "o",
np.where((self.utmp > 18) &
((self.flag != "n") &
(np.char.find(self.flag.astype(str), 'u') == -1) &
(np.char.find(self.flag.astype(str), 'q') == -1)),
self.flag+[","]+["o"], self.flag))
def _adjust_gust(self):
# gustiness adjustment
if (self.gust[0] == 1):
self.wind = np.where(self.dtv >= 0, np.where(self.spd > 0.1, self.spd, 0.1),
np.sqrt(np.power(np.copy(self.spd), 2)+np.power(0.5, 2)))
def __init__(self):
self.meth = "UA"
self.default_gust = [1,1,1000]
class C30(S80):
def set_coolskin_warmlayer(self, wl=0, cskin=1, skin="C35", Rl=None, Rs=None):
self._fix_coolskin_warmlayer(wl, cskin, skin, Rl, Rs)
def __init__(self):
self.meth = "C30"
self.default_gust = [1,1.2,600]
class C35(C30):
def __init__(self):
self.meth = "C35"
self.default_gust = [1,1.2,600]
class ecmwf(C30):
def set_coolskin_warmlayer(self, wl=0, cskin=1, skin="ecmwf", Rl=None, Rs=None):
self._fix_coolskin_warmlayer(wl, cskin, skin, Rl, Rs)
def __init__(self):
self.meth = "ecmwf"
self.default_gust = [1,1,1000]
class Beljaars(C30):
def set_coolskin_warmlayer(self, wl=0, cskin=1, skin="Beljaars", Rl=None, Rs=None):
self._fix_coolskin_warmlayer(wl, cskin, skin, Rl, Rs)
def __init__(self):
self.meth = "Beljaars"
self.default_gust = [1,1,1000]
def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
Rl=None, Rs=None, cskin=None, skin="C35", wl=0, gust=None,
......@@ -141,470 +669,14 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10,
logging.basicConfig(filename='flux_calc.log', filemode="w",
format='%(asctime)s %(message)s', level=logging.INFO)
logging.captureWarnings(True)
# check input values and set defaults where appropriate
lat, hum, P, Rl, Rs, cskin, skin, wl, gust, tol, L, n = get_init(spd, T,
SST, lat,
hum, P,
Rl, Rs,
cskin,
skin,
wl, gust,
L, tol, n,
meth,
qmeth)
ref_ht = 10 # reference height
h_in = get_heights(hin, len(spd)) # heights of input measurements/fields
h_out = get_heights(hout, 1) # desired height of output variables
logging.info('method %s, inputs: lat: %s | P: %s | Rl: %s |'
' Rs: %s | gust: %s | cskin: %s | L : %s', meth,
np.nanmedian(lat), np.round(np.nanmedian(P), 2),
np.round(np.nanmedian(Rl), 2), np.round(np.nanmedian(Rs), 2),
gust, cskin, L)
# set up/calculate temperatures and specific humidities
sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST))
qair, qsea = get_hum(hum, T, sst, P, qmeth)
# mask to preserve missing values when initialising variables
msk = np.empty(sst.shape)
msk = np.where(np.isnan(spd+T+SST), np.nan, 1)
Rb = np.empty(sst.shape)*msk
g = gc(lat)
# specific heat
cp = 1004.67*(1+0.00084*qsea)
th = np.where(T < 200, (np.copy(T)+CtoK) *
np.power(1000/P, 287.1/cp),
np.copy(T)*np.power(1000/P, 287.1/cp)) # potential T
# lapse rate
tlapse = gamma("dry", SST, T, qair/1000, cp)
Ta = np.where(T < 200, np.copy(T)+CtoK+tlapse*h_in[1],
np.copy(T)+tlapse*h_in[1]) # convert to Kelvin if needed
logging.info('method %s and q method %s | qsea:%s, qair:%s', meth, qmeth,
np.round(np.nanmedian(qsea), 7),
np.round(np.nanmedian(qair), 7))
if (np.all(np.isnan(qsea)) or np.all(np.isnan(qair))):
print("qsea and qair cannot be nan")
if ((hum[0] == 'rh') or (hum[0] == 'no')):
rh = hum[1]
elif (hum[0] == 'Td'):
Td = hum[1] # dew point temperature (K)
Td = np.where(Td < 200, np.copy(Td)+CtoK, np.copy(Td))
T = np.where(T < 200, np.copy(T)+CtoK, np.copy(T))
esd = 611.21*np.exp(17.502*((Td-CtoK)/(Td-32.19)))
es = 611.21*np.exp(17.502*((T-CtoK)/(T-32.19)))
rh = 100*esd/es
flag = np.empty(spd.shape, dtype="object")
flag[:] = "n"
if (hum[0] == 'no'):
if (cskin == 1):
flag = np.where(np.isnan(spd+T+SST+P+Rs+Rl), "m", flag)
else:
flag = np.where(np.isnan(spd+T+SST+P), "m", flag)
else:
if (cskin == 1):
flag = np.where(np.isnan(spd+T+SST+hum[1]+P+Rs+Rl), "m", flag)
else:
flag = np.where(np.isnan(spd+T+SST+hum[1]+P), "m", flag)
flag = np.where(rh > 100, "r", flag)
dt = Ta - sst
dq = qair - qsea
# first guesses
t10n, q10n = np.copy(Ta), np.copy(qair)
tv10n = t10n*(1+0.6077*q10n)
# Zeng et al. 1998
tv = th*(1+0.6077*qair) # virtual potential T
dtv = dt*(1+0.6077*qair)+0.6077*th*dq
# ------------
rho = P*100/(287.1*tv10n)
lv = (2.501-0.00237*(sst-CtoK))*1e6 # J/kg
# cskin parameters
tkt = 0.001*np.ones(T.shape)*msk
dter = -0.3*np.ones(T.shape)*msk
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
Rnl = 0.97*(Rl-5.67e-8*np.power(sst-0.3*cskin, 4))
Qs = 0.945*Rs
dtwl = np.ones(T.shape)*0.3*msk
skt = np.copy(sst)
# gustiness adjustment
if (gust[0] == 1 and meth == "UA"):
wind = np.where(dtv >= 0, np.where(spd > 0.1, spd, 0.1),
np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2)))
elif (gust[0] == 1):
wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
elif (gust[0] == 0):
wind = np.copy(spd)
u10n = wind*np.log(10/1e-4)/np.log(hin[0]/1e-4)
usr = 0.035*u10n
cd10n = cdn_calc(u10n, usr, Ta, lat, meth)
psim, psit, psiq = (np.zeros(spd.shape)*msk, np.zeros(spd.shape)*msk,
np.zeros(spd.shape)*msk)
cd = cd_calc(cd10n, h_in[0], ref_ht, psim)
tsr, tsrv, qsr = (np.zeros(spd.shape)*msk, np.zeros(spd.shape)*msk,
np.zeros(spd.shape)*msk)
usr = np.sqrt(cd*np.power(wind, 2))
zo, zot, zoq = (1e-4*np.ones(spd.shape)*msk, 1e-4*np.ones(spd.shape)*msk,
1e-4*np.ones(spd.shape)*msk)
ct10n = np.power(kappa, 2)/(np.log(h_in[0]/zo)*np.log(h_in[1]/zot))
cq10n = np.power(kappa, 2)/(np.log(h_in[0]/zo)*np.log(h_in[2]/zoq))
ct = np.power(kappa, 2)/((np.log(h_in[0]/zo)-psim) *
(np.log(h_in[1]/zot)-psit))
cq = np.power(kappa, 2)/((np.log(h_in[0]/zo)-psim) *
(np.log(h_in[2]/zoq)-psiq))
Rb = g*10*(dtv)/(np.where(T < 200, np.copy(T)+CtoK, np.copy(T)) *
np.power(wind, 2)) # Rb eq. 11 Grachev & Fairall 1997
monob = 1/Rb # eq. 12 Grachev & Fairall 1997
tsr = (dt-dter*cskin-dtwl*wl)*kappa/(np.log(h_in[1]/zot) -
psit_calc(h_in[1]/monob, meth))
qsr = (dq-dqer*cskin)*kappa/(np.log(h_in[2]/zoq) -
psit_calc(h_in[2]/monob, meth))
it, ind = 0, np.where(spd > 0)
ii, itera = True, -1*np.ones(spd.shape)*msk
tau = 0.05*np.ones(spd.shape)*msk
sensible = -5*np.ones(spd.shape)*msk
latent = -65*np.ones(spd.shape)*msk
# iteration loop
while np.any(ii):
it += 1
if it > n:
break
if (tol[0] == 'flux'):
old = np.array([np.copy(tau), np.copy(sensible), np.copy(latent)])
elif (tol[0] == 'ref'):
old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n)])
elif (tol[0] == 'all'):
old = np.array([np.copy(u10n), np.copy(t10n), np.copy(q10n),
np.copy(tau), np.copy(sensible), np.copy(latent)])
cd10n[ind] = cdn_calc(u10n[ind], usr[ind], Ta[ind], lat[ind], meth)
if (np.all(np.isnan(cd10n))):
break
logging.info('break %s at iteration %s cd10n<0', meth, it)
zo[ind] = ref_ht/np.exp(kappa/np.sqrt(cd10n[ind]))
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
cd[ind] = cd_calc(cd10n[ind], h_in[0, ind], ref_ht, psim[ind])
ct10n[ind], cq10n[ind] = ctcqn_calc(h_in[1, ind]/monob[ind],
cd10n[ind], usr[ind], zo[ind],
Ta[ind], meth)
zot[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
(ct10n[ind]*np.log(ref_ht/zo[ind]))))
zoq[ind] = ref_ht/(np.exp(np.power(kappa, 2) /
(cq10n[ind]*np.log(ref_ht/zo[ind]))))
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
ct[ind], cq[ind] = ctcq_calc(cd10n[ind], cd[ind], ct10n[ind],
cq10n[ind], h_in[:, ind],
[ref_ht, ref_ht, ref_ht], psit[ind],
psiq[ind])
if (meth == "NCAR"):
cd = np.maximum(np.copy(cd), 1e-4)
ct = np.maximum(np.copy(ct), 1e-4)
cq = np.maximum(np.copy(cq), 1e-4)
zo = np.minimum(np.copy(zo), 0.0025)
usr[ind], tsr[ind], qsr[ind] = get_strs(h_in[:, ind], monob[ind],
wind[ind], zo[ind], zot[ind],
zoq[ind], dt[ind], dq[ind],
dter[ind], dqer[ind],
dtwl[ind], ct[ind], cq[ind],
cskin, wl, meth)
if ((cskin == 1) and (wl == 0)):
if (skin == "C35"):
dter[ind], dqer[ind], tkt[ind] = cs_C35(np.copy(sst[ind]),
qsea[ind],
rho[ind], Rs[ind],
Rnl[ind],
cp[ind], lv[ind],
np.copy(tkt[ind]),
usr[ind], tsr[ind],
qsr[ind], lat[ind])
elif (skin == "ecmwf"):
dter[ind] = cs_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), lat[ind])
dqer[ind] = (dter[ind]*0.622*lv[ind]*qsea[ind] /
(287.1*np.power(sst[ind], 2)))
elif (skin == "Beljaars"):
Qs[ind], dter[ind] = cs_Beljaars(rho[ind], Rs[ind], Rnl[ind],
cp[ind], lv[ind], usr[ind],
tsr[ind], qsr[ind], lat[ind],
np.copy(Qs[ind]))
dqer = dter*0.622*lv*qsea/(287.1*np.power(sst, 2))
skt = np.copy(sst)+dter
elif ((cskin == 1) and (wl == 1)):
if (skin == "C35"):
dter[ind], dqer[ind], tkt[ind] = cs_C35(sst[ind], qsea[ind],
rho[ind], Rs[ind],
Rnl[ind],
cp[ind], lv[ind],
np.copy(tkt[ind]),
usr[ind], tsr[ind],
qsr[ind], lat[ind])
dtwl[ind] = wl_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
skt = np.copy(sst)+dter+dtwl
elif (skin == "ecmwf"):
dter[ind] = cs_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
sst[ind], lat[ind])
dtwl[ind] = wl_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
skt = np.copy(sst)+dter+dtwl
dqer[ind] = (dter[ind]*0.622*lv[ind]*qsea[ind] /
(287.1*np.power(skt[ind], 2)))
elif (skin == "Beljaars"):
Qs[ind], dter[ind] = cs_Beljaars(rho[ind], Rs[ind], Rnl[ind],
cp[ind], lv[ind], usr[ind],
tsr[ind], qsr[ind], lat[ind],
np.copy(Qs[ind]))
dtwl[ind] = wl_ecmwf(rho[ind], Rs[ind], Rnl[ind], cp[ind],
lv[ind], usr[ind], tsr[ind], qsr[ind],
np.copy(sst[ind]), np.copy(skt[ind]),
np.copy(dter[ind]), lat[ind])
skt = np.copy(sst)+dter+dtwl
dqer = dter*0.622*lv*qsea/(287.1*np.power(skt, 2))
else:
dter[ind] = np.zeros(sst[ind].shape)
dqer[ind] = np.zeros(sst[ind].shape)
tkt[ind] = 0.001*np.ones(T[ind].shape)
logging.info('method %s | dter = %s | dqer = %s | tkt = %s | Rnl = %s '
'| usr = %s | tsr = %s | qsr = %s', meth,
np.round(np.nanmedian(dter), 2),
np.round(np.nanmedian(dqer), 7),
np.round(np.nanmedian(tkt), 2),
np.round(np.nanmedian(Rnl), 2),
np.round(np.nanmedian(usr), 3),
np.round(np.nanmedian(tsr), 4),
np.round(np.nanmedian(qsr), 7))
Rnl[ind] = 0.97*(Rl[ind]-5.67e-8*np.power(sst[ind] +
dter[ind]*cskin, 4))
t10n[ind] = (Ta[ind] -
tsr[ind]/kappa*(np.log(h_in[1, ind]/ref_ht)-psit[ind]))
q10n[ind] = (qair[ind] -
qsr[ind]/kappa*(np.log(h_in[2, ind]/ref_ht)-psiq[ind]))
tv10n[ind] = t10n[ind]*(1+0.6077*q10n[ind])
tsrv[ind], monob[ind], Rb[ind] = get_L(L, lat[ind], usr[ind], tsr[ind],
qsr[ind], h_in[:, ind], Ta[ind],
(sst[ind]+dter[ind]*cskin +
dtwl[ind]*wl), qair[ind],
qsea[ind], wind[ind],
np.copy(monob[ind]), zo[ind],
zot[ind], psim[ind], meth)
psim[ind] = psim_calc(h_in[0, ind]/monob[ind], meth)
psit[ind] = psit_calc(h_in[1, ind]/monob[ind], meth)
psiq[ind] = psit_calc(h_in[2, ind]/monob[ind], meth)
if (gust[0] == 1 and meth == "UA"):
wind[ind] = np.where(dtv[ind] >= 0, np.where(spd[ind] > 0.1,
spd[ind], 0.1),
np.sqrt(np.power(np.copy(spd[ind]), 2) +
np.power(get_gust(gust[1], tv[ind],
usr[ind], tsrv[ind],
gust[2], lat[ind]),
2))) # Zeng et al. 1998 (20)
elif (gust[0] == 1):
wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
np.power(get_gust(gust[1], Ta[ind], usr[ind],
tsrv[ind], gust[2],
lat[ind]), 2))
elif (gust[0] == 0):
wind[ind] = np.copy(spd[ind])
u10n[ind] = wind[ind]-usr[ind]/kappa*(np.log(h_in[0, ind]/10) -
psim[ind])
if (it < 4): # make sure you allow small negative values convergence
u10n = np.where(u10n < 0, 0.5, u10n)
utmp = np.copy(u10n)
utmp = np.where(utmp < 0, np.nan, utmp)
itera[ind] = np.ones(1)*it
tau = rho*np.power(usr, 2)*(spd/wind)
sensible = rho*cp*usr*tsr
latent = rho*lv*usr*qsr
if (tol[0] == 'flux'):
new = np.array([np.copy(tau), np.copy(sensible), np.copy(latent)])
elif (tol[0] == 'ref'):
new = np.array([np.copy(utmp), np.copy(t10n), np.copy(q10n)])
elif (tol[0] == 'all'):
new = np.array([np.copy(utmp), np.copy(t10n), np.copy(q10n),
np.copy(tau), np.copy(sensible), np.copy(latent)])
if (it > 2): # force at least two iterations
d = np.abs(new-old)
if (tol[0] == 'flux'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3]))
elif (tol[0] == 'ref'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3]))
elif (tol[0] == 'all'):
ind = np.where((d[0, :] > tol[1])+(d[1, :] > tol[2]) +
(d[2, :] > tol[3])+(d[3, :] > tol[4]) +
(d[4, :] > tol[5])+(d[5, :] > tol[6]))
if (ind[0].size == 0):
ii = False
else:
ii = True
itera[ind] = -1
itera = np.where(itera > n, -1, itera)
logging.info('method %s | # of iterations:%s', meth, it)
logging.info('method %s | # of points that did not converge :%s \n', meth,
ind[0].size)
# calculate output parameters
rho = (0.34838*P)/(tv10n)
t10n = t10n-(273.16+tlapse*ref_ht)
# solve for zo from cd10n
zo = ref_ht/np.exp(kappa/np.sqrt(cd10n))
# adjust neutral cdn at any output height
cdn = np.power(kappa/np.log(hout/zo), 2)
cd = cd_calc(cdn, h_out[0], h_out[0], psim)
# solve for zot, zoq from ct10n, cq10n
zot = ref_ht/(np.exp(kappa**2/(ct10n*np.log(ref_ht/zo))))
zoq = ref_ht/(np.exp(kappa**2/(cq10n*np.log(ref_ht/zo))))
# adjust neutral ctn, cqn at any output height
ctn = np.power(kappa, 2)/(np.log(h_out[0]/zo)*np.log(h_out[1]/zot))
cqn = np.power(kappa, 2)/(np.log(h_out[0]/zo)*np.log(h_out[2]/zoq))
ct, cq = ctcq_calc(cdn, cd, ctn, cqn, h_out, h_out, psit, psiq)
uref = (spd-usr/kappa*(np.log(h_in[0]/h_out[0])-psim +
psim_calc(h_out[0]/monob, meth)))
tref = (Ta-tsr/kappa*(np.log(h_in[1]/h_out[1])-psit +
psit_calc(h_out[0]/monob, meth)))
tref = tref-(CtoK+tlapse*h_out[1])
qref = (qair-qsr/kappa*(np.log(h_in[2]/h_out[2]) -
psit+psit_calc(h_out[2]/monob, meth)))
if (wl == 0):
dtwl = np.zeros(T.shape)*msk # reset to zero if not used
flag = np.where((u10n < 0) & (flag == "n"), "u",
np.where((u10n < 0) &
(np.char.find(flag.astype(str), 'u') == -1),
flag+[","]+["u"], flag))
flag = np.where((q10n < 0) & (flag == "n"), "q",
np.where((q10n < 0) & (flag != "n"), flag+[","]+["q"],
flag))
flag = np.where(((Rb < -0.5) | (Rb > 0.2) | ((hin[0]/monob) > 1000)) &
(flag == "n"), "l",
np.where(((Rb < -0.5) | (Rb > 0.2) |
((hin[0]/monob) > 1000)) &
(flag != "n"), flag+[","]+["l"], flag))
if (out == 1):
flag = np.where((itera == -1) & (flag == "n"), "i",
np.where((itera == -1) &
((flag != "n") &
(np.char.find(flag.astype(str), 'm') == -1)),
flag+[","]+["i"], flag))
else:
flag = np.where((itera == -1) & (flag == "n"), "i",
np.where((itera == -1) &
((flag != "n") &
(np.char.find(flag.astype(str), 'm') == -1) &
(np.char.find(flag.astype(str), 'u') == -1)),
flag+[","]+["i"], flag))
if (meth == "S80"):
flag = np.where(((utmp < 6) | (utmp > 22)) & (flag == "n"), "o",
np.where(((utmp < 6) | (utmp > 22)) &
((flag != "n") &
(np.char.find(flag.astype(str), 'u') == -1) &
(np.char.find(flag.astype(str), 'q') == -1)),
flag+[","]+["o"], flag))
elif (meth == "LP82"):
flag = np.where(((utmp < 3) | (utmp > 25)) & (flag == "n"), "o",
np.where(((utmp < 3) | (utmp > 25)) &
((flag != "n") &
(np.char.find(flag.astype(str), 'u') == -1) &
(np.char.find(flag.astype(str), 'q') == -1)),
flag+[","]+["o"], flag))
elif (meth == "YT96"):
flag = np.where(((utmp < 0) | (utmp > 26)) & (flag == "n"), "o",
np.where(((utmp < 3) | (utmp > 26)) &
((flag != "n") &
(np.char.find(flag.astype(str), 'u') == -1) &
(np.char.find(flag.astype(str), 'q') == -1)),
flag+[","]+["o"], flag))
elif (meth == "UA"):
flag = np.where((utmp > 18) & (flag == "n"), "o",
np.where((utmp > 18) &
((flag != "n") &
(np.char.find(flag.astype(str), 'u') == -1) &
(np.char.find(flag.astype(str), 'q') == -1)),
flag+[","]+["o"], flag))
elif (meth == "NCAR"):
flag = np.where((utmp < 0.5) & (flag == "n"), "o",
np.where((utmp < 0.5) &
((flag != "n") &
(np.char.find(flag.astype(str), 'u') == -1) &
(np.char.find(flag.astype(str), 'q') == -1)),
flag+[","]+["o"], flag))
# Do not calculate lhf if a measure of humidity is not input
if (hum[0] == 'no'):
latent = np.ones(sst.shape)*np.nan
qsr = np.ones(sst.shape)*np.nan
q10n = np.ones(sst.shape)*np.nan
qref = np.ones(sst.shape)*np.nan
qair = np.ones(sst.shape)*np.nan
rh = np.ones(sst.shape)*np.nan
res = np.zeros((41, 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][:] = psiq
res[17][:] = u10n
res[18][:] = t10n
res[19][:] = tv10n
res[20][:] = q10n
res[21][:] = zo
res[22][:] = zot
res[23][:] = zoq
res[24][:] = uref
res[25][:] = tref
res[26][:] = qref
res[27][:] = itera
res[28][:] = dter
res[29][:] = dqer
res[30][:] = dtwl
res[31][:] = qair
res[32][:] = qsea
res[33][:] = Rl
res[34][:] = Rs
res[35][:] = Rnl
res[36][:] = np.sqrt(np.power(wind, 2)-np.power(spd, 2))
res[37][:] = Rb
res[38][:] = rh
res[39][:] = tkt
res[40][:] = lv
if (out == 0):
res[:, ind] = np.nan
# set missing values where data have non acceptable values
if (hum[0] != 'no'):
res = np.asarray([np.where(q10n < 0, np.nan,
res[i][:]) for i in range(41)])
res = np.asarray([np.where(u10n < 0, np.nan,
res[i][:]) for i in range(41)])
elif (out == 1):
print("Warning: the output will contain values for points that have"
" not converged and negative values (if any) for u10n/q10n")
# output with pandas
resAll = pd.DataFrame(data=res.T, index=range(len(spd)),
columns=["tau", "shf", "lhf", "L", "cd", "cdn", "ct",
"ctn", "cq", "cqn", "tsrv", "tsr", "qsr",
"usr", "psim", "psit", "psiq", "u10n",
"t10n", "tv10n", "q10n", "zo", "zot", "zoq",
"uref", "tref", "qref", "iteration", "dter",
"dqer", "dtwl", "qair", "qsea", "Rl", "Rs",
"Rnl", "ug", "Rib", "rh", "delta", "lv"])
resAll["flag"] = flag
iclass = globals()[meth]()
iclass.add_gust(gust)
iclass.add_variables(spd, T, SST, lat=lat, hum=hum, P=P, L=L)
iclass.get_heights(hin, hout)
iclass.get_specHumidity(qmeth=qmeth)
iclass.set_coolskin_warmlayer(wl=wl, cskin=cskin,skin=skin,Rl=Rl,Rs=Rs)
iclass.iterate(n=n,tol=tol)
resAll = iclass.get_output(out=out)
return resAll
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