From bb9632b573744f2ee3ec7105052fafaadc07c16c Mon Sep 17 00:00:00 2001
From: sbiri <sbiri@noc.ac.uk>
Date: Mon, 14 Jun 2021 14:30:59 +0100
Subject: [PATCH] changed gamma_moist in hum_subs to gamma for different
 options of adiabatic lapse rate (dry-constant, definition, or for unsaturated
 air with water vapor-, moist) adjusted it in ASFC

---
 AirSeaFluxCode.py |  7 ++++---
 hum_subs.py       | 32 +++++++++++++++++++++++---------
 2 files changed, 27 insertions(+), 12 deletions(-)

diff --git a/AirSeaFluxCode.py b/AirSeaFluxCode.py
index 814c82d..d0df1d5 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 cba06af..89e3e13 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
 #------------------------------------------------------------------------------
-- 
GitLab