diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py index 814c82d0b0a979f60c0cd84b1aa53e430c94df04..d0df1d52231e9fca260c4bc42fd38a19c1efe875 100644 --- a/AirSeaFluxCode.py +++ b/AirSeaFluxCode.py @@ -2,7 +2,7 @@ import numpy as np import pandas as pd import logging from get_init import get_init -from hum_subs import (get_hum, gamma_moist) +from hum_subs import (get_hum, gamma) from util_subs import (kappa, CtoK, get_heights) from flux_subs import (cs_C35, cs_Beljaars, cs_ecmwf, wl_ecmwf, get_gust, get_L, get_strs, psim_calc, @@ -162,8 +162,10 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, sst = np.where(SST < 200, np.copy(SST)+CtoK, np.copy(SST)) qair, qsea = get_hum(hum, T, sst, P, qmeth) Rb = np.empty(sst.shape) + # specific heat + cp = 1004.67*(1+0.00084*qsea) #1005*np.ones(spd.shape)# #lapse rate - tlapse = gamma_moist(SST, T, qair/1000) + 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, @@ -201,7 +203,6 @@ def AirSeaFluxCode(spd, T, SST, lat=None, hum=None, P=None, hin=18, hout=10, # ------------ rho = P*100/(287.1*tv10n) lv = (2.501-0.00237*(sst-CtoK))*1e6 - cp = 1005*np.ones(spd.shape)#1004.67*(1 + 0.00084*qsea) u10n = np.copy(spd) usr = 0.035*u10n cd10n = cdn_calc(u10n, usr, Ta, lat, meth) diff --git a/hum_subs.py b/hum_subs.py index cba06afcde4ff2524cb2ccde2cf5971541395e8b..89e3e1384a8cb3d395f9ed826da8b08ee9662ea5 100644 --- a/hum_subs.py +++ b/hum_subs.py @@ -400,12 +400,16 @@ def get_hum(hum, T, sst, P, qmeth): #------------------------------------------------------------------------------ -def gamma_moist(sst, t, q): +def gamma(opt, sst, t, q, cp): """ - Computes the moist adiabatic lapse-rate + Computes the adiabatic lapse-rate Parameters ---------- + opt : str + type of adiabatic lapse rate dry or "moist" + dry has options to be constant "dry_c", for dry air "dry", or + for unsaturated air with water vapor "dry_v" sst : float sea surface temperature [K] t : float @@ -423,12 +427,22 @@ def gamma_moist(sst, t, q): sst = sst+CtoK if (np.nanmin(t) < 200): # if sst in Celsius convert to Kelvin t = t+CtoK - t = np.maximum(t, 180) - q = np.maximum(q, 1e-6) - w = q/(1-q) # mixing ratio w = q/(1-q) - iRT = 1/(287.05*t) - # latent heat of vaporization of water as a function of temperature - lv = (2.501-0.00237*(sst-CtoK))*1e6 - gamma = 9.8*(1+lv*w*iRT)/(1005+np.power(lv, 2)*w*(287.05/461.495)*iRT/t) + if (opt == "moist"): + t = np.maximum(t, 180) + q = np.maximum(q, 1e-6) + w = q/(1-q) # mixing ratio w = q/(1-q) + iRT = 1/(287.05*t) + # latent heat of vaporization of water as a function of temperature + lv = (2.501-0.00237*(sst-CtoK))*1e6 + gamma = 9.8*(1+lv*w*iRT)/(1005+np.power(lv, 2)*w*(287.05/461.495) * + iRT/t) + elif (opt == "dry_c"): + gamma = 0.0098*np.ones(t.shape) + elif (opt == "dry"): + gamma = 9.81/cp + elif (opt == "dry_v"): + w = q/(1-q) # mixing ratio + f_v = 1-0.85*w #(1+w)/(1+w*) + gamma = f_v*9.81/cp return gamma #------------------------------------------------------------------------------