diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
deleted file mode 100644
index 8d6b89347c4233f84fad9be47285f550dc422c82..0000000000000000000000000000000000000000
--- a/AirSeaFluxCode.py
+++ /dev/null
@@ -1,503 +0,0 @@
-import numpy as np
-import sys
-import logging
-#import metpy.constants as mpcon
-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)
-
-
-def AirSeaFluxCode(spd, T, SST, lat, RH, P, hin, hout, 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","LY04","C30","C35","C40","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
- 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 COARE
- n : int
- number of iterations
-
- 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',
- 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
- if np.all(np.isnan(lat)): # set latitude to 45deg if empty
- lat=45*np.ones(spd.shape)
- 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)) and meth == "C35"):
- RH = np.ones(spd.shape)*80 # if empty set to default for COARE3.5
- elif (np.all(np.isnan(RH))):
- sys.exit("input RH is empty")
- logging.debug('input RH is empty')
- if (np.all(np.isnan(P)) and (meth == "C30" or meth == "C40")):
- P = np.ones(spd.shape)*1015 # if empty set to default for COARE3.0
- elif ((np.all(np.isnan(P))) and meth == "C35"):
- P = np.ones(spd.shape)*1013 # if empty set to default for COARE3.5
- elif (np.all(np.isnan(P))):
- sys.exit("input P is empty")
- logging.debug('input P is empty')
- if (np.all(np.isnan(Rl)) and meth == "C30"):
- Rl = np.ones(spd.shape)*150 # set to default for COARE3.0
- elif ((np.all(np.isnan(Rl)) and meth == "C35") or
- (np.all(np.isnan(Rl)) and meth == "C40")):
- Rl = np.ones(spd.shape)*370 # set to default for COARE3.5
- if (np.all(np.isnan(Rs)) and meth == "C30"):
- Rs = np.ones(spd.shape)*370 # set to default for COARE3.0
- elif ((np.all(np.isnan(Rs))) and (meth == "C35" or meth == "C40")):
- Rs = np.ones(spd.shape)*150 # set to default for COARE3.5
- if ((np.all(np.isnan(zi))) and (meth == "C30" or meth == "C35" or
- meth == "C40")):
- zi = 600 # set to default for COARE3.5
- elif ((np.all(np.isnan(zi))) and (meth == "ERA5" or meth == "UA")):
- zi = 1000
- if (np.all(np.isnan(lat)) and (meth == "C30" or meth == "C35" or
- meth == "C40")):
- lat=45*np.ones(np.shape(spd))
- ####
- th = np.where(np.nanmax(T) < 200, (np.copy(T)+CtoK) *
- np.power(1000/P,287.1/1004.67),
- np.copy(T)*np.power(1000/P,287.1/1004.67)) # potential T
- 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 == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
- meth == "ERA5"):
- 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', 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.info('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
- dt = Ta - sst
- dq = qair - qsea
- t10n, q10n = np.copy(Ta), np.copy(qair)
- tv10n = t10n*(1 + 0.61*q10n)
- # Zeng et al. 1998
- tv=th*(1.+0.61*qair) # virtual potential T
- dtv=dt*(1.+0.61*qair)+0.61*th*dq
- # ------------
-# rho = P*100/(287.1*tv10n)
- rho = P*100/(287.1*(T+CtoK)*(1+0.61*qair))
- # rho=P*100/(287.1*sst*(1+0.61*qsea)) # Zeng et al. 1998
- lv = (2.501-0.00237*SST)*1e6
- cp = 1004.67*(1 + 0.00084*qsea)
-# cp = 1004.67 # Zeng et al. 1998, C3.0, C3.5
- u10n = np.copy(spd)
- monob = -100*np.ones(u10n.shape) # Monin-Obukhov length
- cdn = cdn_calc(u10n, Ta, None, lat, meth)
- psim, psit, psiq = (np.zeros(u10n.shape), np.zeros(u10n.shape),
- np.zeros(u10n.shape))
- 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)
- if (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)))
- usr = 0.06
- for i in range(5):
- zo = 0.013*np.power(usr,2)/g+0.11*visc_air(Ta)/usr
- usr=kappa*wind/np.log(hh_in[0]/zo)
- Rb = g*hh_in[0]*dtv/(tv*wind**2)
- zol = np.where(Rb >= 0, Rb*np.log(hh_in[0]/zo) /
- (1-5*np.where(Rb < 0.19, Rb, 0.19)),
- Rb*np.log(hh_in[0]/zo))
- monob = hh_in[0]/zol
- zo = 0.013*np.power(usr, 2)/g + 0.11*visc_air(Ta)/usr
- zot = zo/np.exp(2.67*np.power(usr*zo/visc_air(Ta), 0.25)-2.57)
- zoq = zot
- logging.info('method %s | wind:%s, usr:%s,'
- 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
- np.ma.median(wind[~np.isnan(wind)]),
- np.ma.median(usr[~np.isnan(usr)]),
- np.ma.median(zo[~np.isnan(zo)]),
- np.ma.median(zot[~np.isnan(zot)]),
- np.ma.median(zoq[~np.isnan(zoq)]),
- np.ma.median(Rb[~np.isnan(Rb)]),
- np.ma.median(monob[~np.isnan(monob)]))
- elif (meth == "ERA5"):
- wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
- usr = np.sqrt(cd*wind**2)
- Rb = ((g*hh_in[0]*((2*dt)/(Ta+sst-g*hh_in[0]) +
- 0.61*dq))/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 = (Rb*((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))))
- monob = hh_in[0]/zol
- logging.info('method %s | wind:%s, usr:%s,'
- 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
- np.ma.median(wind[~np.isnan(wind)]),
- np.ma.median(usr[~np.isnan(usr)]),
- np.ma.median(zo[~np.isnan(zo)]),
- np.ma.median(zot[~np.isnan(zot)]),
- np.ma.median(Rb[~np.isnan(Rb)]),
- np.ma.median(monob[~np.isnan(monob)]))
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
- wind = np.sqrt(np.power(np.copy(spd), 2)+np.power(0.5, 2))
- usr = 0.035*wind*np.log(10/1e-4)/np.log(hh_in[0]/1e-4)
- a = 0.011*np.ones(T.shape)
- a = np.where(wind > 10, 0.011+(wind-10)/(18-10)*(0.018-0.011),
- np.where(wind > 18, 0.018, a))
- zo = a*np.power(usr, 2)/g+0.11*visc_air(T)/usr
- rr = zo*usr/visc_air(T)
- zoq = np.minimum(5e-5/np.power(rr, 0.6), 1.15e-4)
- zot=zoq
- Rb = g*hh_in[0]*dtv/((T+CtoK)*np.power(wind, 2))
- zol = (Rb*((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))))
- monob = hh_in[0]/zol
-# wetc = 0.622*lv*qsea/(287.1*np.power(sst, 2))
- tkt = 0.001*np.ones(T.shape)
- Rnl = 0.97*(5.67e-8*np.power(sst-0.3*jcool+CtoK, 4)-Rl)
- dter, dqer = np.ones(T.shape)*0.3, np.zeros(T.shape)*np.nan
- logging.info('method %s | wind:%s, usr:%s,'
- 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
- np.ma.median(wind[~np.isnan(wind)]),
- np.ma.median(usr[~np.isnan(usr)]),
- np.ma.median(zo[~np.isnan(zo)]),
- np.ma.median(zot[~np.isnan(zot)]),
- np.ma.median(Rb[~np.isnan(Rb)]),
- np.ma.median(monob[~np.isnan(monob)]))
- else:
- wind = np.copy(spd)
- zo, zot = 0.0001*np.ones(u10n.shape), 0.0001*np.ones(u10n.shape)
- usr = np.sqrt(cd * wind**2)
- logging.info('method %s | wind:%s, usr:%s,'
- 'zo:%s, zot:%s, zoq:%s, Rb:%s, monob:%s', meth,
- np.ma.median(wind[~np.isnan(wind)]),
- np.ma.median(usr[~np.isnan(usr)]),
- np.ma.median(zo[~np.isnan(zo)]),
- np.ma.median(zot[~np.isnan(zot)]),
- np.ma.median(monob[~np.isnan(monob)]))
- # tolerance for u,t,q,usr,tsr,qsr
- tol = np.array([1e-06, 0.01, 5e-05, 1e-06, 0.001, 5e-07])
- it, ind, ii, itera = 0, np.where(spd > 0), True, np.zeros(spd.shape)*np.nan
- while np.any(ii):
- it += 1
- if it > n:
- break
- 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, lat[ind], meth)
- if (np.all(np.isnan(cdn))):
- break # sys.exit("cdn cannot be nan")
- logging.info('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], 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], 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])
- if (meth == "UA"):
- usr[ind] = np.where(hh_in[0]/monob[ind] < -1.574, kappa*wind[ind] /
- (np.log(-1.574*monob[ind]/zo[ind]) -
- psim_calc(-1.574, meth) +
- psim_calc(zo[ind]/monob[ind], meth) +
- 1.14*(np.power(-hh_in[0]/monob[ind],1/3) -
- np.power(1.574,1/3))),
- np.where((hh_in[0]/monob[ind] > -1.574) &
- (hh_in[0]/monob[ind] < 0),
- kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) -
- psim_calc(hh_in[0]/monob[ind], meth) +
- psim_calc(zo[ind]/monob[ind], meth)),
- np.where((hh_in[0]/monob[ind] > 0) &
- (hh_in[0]/monob[ind] < 1),
- kappa*wind[ind]/(np.log(hh_in[0]/zo[ind]) +
- 5*hh_in[0]/monob[ind]-5*zo[ind]/monob[ind]),
- kappa*wind[ind]/(np.log(monob[ind]/zo[ind]) +
- 5-5*zo[ind]/monob[ind] +
- 5*np.log(hh_in[0]/monob[ind]) +
- hh_in[0]/monob[ind]-1))))
- # Zeng et al. 1998 (7-10)
- tsr[ind] = np.where(hh_in[1]/monob[ind] < -0.465, kappa*dt[ind] /
- (np.log((-0.465*monob[ind])/zot[ind]) -
- psit_calc(-0.465, meth)+0.8 *
- (np.power(0.465,-1/3) -
- np.power(-hh_in[1]/monob[ind],-1/3))),
- np.where((hh_in[1]/monob[ind]>-0.465) &
- (hh_in[1]/monob[ind]<0),
- kappa*dt[ind]/(np.log(hh_in[1]/zot[ind]) -
- psit_calc(hh_in[1]/monob[ind], meth) +
- psit_calc(zot[ind]/monob[ind], meth)),
- np.where((hh_in[1]/monob[ind]>0) &
- (hh_in[1]/monob[ind]<1),
- kappa*dt[ind]/(np.log(hh_in[1]/zot[ind]) +
- 5*hh_in[1]/monob[ind]-5*zot[ind]/monob[ind]),
- kappa*dt[ind]/(np.log(monob[ind]/zot[ind])+5 -
- 5**zot[ind]/monob[ind] +
- 5*np.log(hh_in[1]/monob[ind]) +
- hh_in[1]/monob[ind]-1))))
- # Zeng et al. 1998 (11-14)
- qsr[ind] = np.where(hh_in[2]/monob[ind] < -0.465, kappa*dq[ind] /
- (np.log((-0.465*monob[ind])/zoq[ind]) -
- psit_calc(-0.465, meth) +
- psit_calc(zoq[ind]/monob[ind], meth) +
- 0.8*(np.power(0.465,-1/3) -
- np.power(-hh_in[2]/monob[ind],-1/3))),
- np.where((hh_in[2]/monob[ind]>-0.465) &
- (hh_in[2]/monob[ind]<0),
- kappa*dq[ind]/(np.log(hh_in[1]/zot[ind]) -
- psit_calc(hh_in[2]/monob[ind], meth) +
- psit_calc(zoq[ind]/monob[ind], meth)),
- np.where((hh_in[2]/monob[ind]>0) &
- (hh_in[2]/monob[ind]<1), kappa*dq[ind] /
- (np.log(hh_in[1]/zoq[ind]) +
- 5*hh_in[2]/monob[ind]-5*zoq[ind]/monob[ind]),
- kappa*dq[ind]/(np.log(monob[ind]/zoq[ind])+5 -
- 5*zoq[ind]/monob[ind] +
- 5*np.log(hh_in[2]/monob[ind]) +
- hh_in[2]/monob[ind]-1))))
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
- usr[ind] = (wind[ind]*kappa/(np.log(hh_in[0]/zo[ind]) -
- psiu_26(hh_in[0]/monob[ind], meth)))
- dter[ind], dqer[ind], tkt[ind] = get_skin(sst[ind], qsea[ind],
- rho[ind], Rl[ind],
- Rs[ind], Rnl[ind],
- cp[ind], lv[ind],
- np.copy(tkt[ind]),
- usr[ind], tsr[ind],
- qsr[ind], lat[ind])
- logging.info('method %s | dter = %s | Rnl = %s '
- '| usr = %s | tsr = %s | qsr = %s', meth,
- np.ma.median(dter[~np.isnan(dter)]),
- np.ma.median(Rnl[~np.isnan(Rnl)]),
- np.ma.median(usr[~np.isnan(usr)]),
- np.ma.median(tsr[~np.isnan(tsr)]),
- np.ma.median(qsr[~np.isnan(qsr)]))
- qsr[ind] = ((dq[ind]+dqer[ind]*jcool)*(kappa /
- (np.log(hin[2]/zoq[ind])-psit_26(hin[2]/monob[ind]))))
- tsr[ind] = ((dt[ind]+dter[ind]*jcool)*(kappa /
- (np.log(hin[1]/zot[ind])-psit_26(hin[1]/monob[ind]))))
- Rnl[ind] = 0.97*(5.67e-8*np.power(SST[ind] -
- dter[ind]*jcool+CtoK, 4)-Rl[ind])
- else:
- 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])
- if (meth == "UA"):
- tsrv[ind] = tsr[ind]*(1.+0.61*qair[ind])+0.61*th[ind]*qsr[ind]
- monob[ind] = (tv[ind]*np.power(usr[ind], 2))/(kappa*g[ind]*tsrv[ind])
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
- tsrv[ind] = tsr[ind]+0.61*(T[ind]+CtoK)*qsr[ind]
- zol[ind] = (kappa*g[ind]*hh_in[0]/(T[ind]+CtoK)*(tsr[ind] +
- 0.61*(T[ind]+CtoK)*qsr[ind])/np.power(usr[ind], 2))
- monob[ind] = hh_in[0]/zol[ind]
- elif (meth == "ERA5"):
- tsrv[ind] = tsr[ind]+0.61*t10n[ind]*qsr[ind]
- Rb[ind] = ((g[ind]*hh_in[0]*((2*dt[ind])/(Ta[ind]+sst[ind]-g[ind] *
- hh_in[0])+0.61*dq[ind]))/np.power(wind[ind], 2))
- zo[ind] = (0.11*visc_air(Ta[ind])/usr[ind]+0.018 *
- np.power(usr[ind], 2)/g[ind])
- zot[ind] = 0.40*visc_air(Ta[ind])/usr[ind]
- zol[ind] = (Rb[ind]*((np.log((hh_in[0]+zo[ind])/zo[ind]) -
- psim_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
- psim_calc(zo[ind]/monob[ind], meth)) /
- (np.log((hh_in[0]+zo[ind])/zot[ind]) -
- psit_calc((hh_in[0]+zo[ind])/monob[ind], meth) +
- psit_calc(zot[ind]/monob[ind], meth))))
- monob[ind] = hh_in[0]/zol[ind]
- else:
- 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], meth)
- psit[ind] = psit_calc(hh_in[1]/monob[ind], meth)
- psiq[ind] = psit_calc(hh_in[2]/monob[ind], meth)
- if (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(1,tv[ind], usr[ind],
- tsrv[ind], zi, lat[ind]), 2)))
- # Zeng et al. 1998 (20)
- u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
- psim[ind]))
- u10n[u10n < 0] = np.nan
- elif (meth == "C30" or meth == "C35" or meth == "C40"):
- wind[ind] = np.sqrt(np.power(np.copy(spd[ind]), 2) +
- np.power(get_gust(1.2, Ta[ind], usr[ind],
- tsrv[ind], zi, lat[ind]), 2))
- u10n[ind] = ((wind[ind] + usr[ind]/kappa*(np.log(10/hh_in[0])-
- psiu_26(10/monob[ind], meth) +
- psiu_26(hh_in[0]/monob[ind], meth)) +
- psiu_26(10/monob[ind], meth)*usr[ind]/kappa /
- (wind[ind]/spd[ind])))
- u10n[u10n < 0] = np.nan
- elif (meth == "ERA5"):
- 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))
- u10n[ind] = (spd[ind]+(usr[ind]/kappa)*(np.log(hh_in[0] /
- ref_ht)-psim[ind]))
- u10n[u10n < 0] = np.nan
- else:
- u10n[ind] = (wind[ind]-(usr[ind]/kappa)*(np.log(hh_in[0]/10) -
- psim[ind]))
- u10n[u10n < 0] = np.nan
- 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]))
- logging.info('method %s | # of iterations:%s', meth, np.ma.median(it))
- logging.info('method %s | # of points that did not converge :%s', meth,
- np.shape(ind))
- # calculate output parameters
-# rho = (0.34838*P)/(tv10n)
- rho = P*100./(287.1*(T+CtoK)*(1+0.61*qair)) # C35
- t10n = t10n-(273.16+tlapse*ref_ht)
- sensible = -1*tsr*usr*cp*rho
- latent = -1*qsr*usr*lv*rho
- if (meth == "C30" or meth == "C35" or meth == "C40" or meth == "UA" or
- meth == "ERA5"):
- tau = rho*np.power(usr, 2)*(spd/wind)
- else:
- tau = rho*np.power(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, 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, 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