From 8222efbcfedd945f56cab17bf5d7308e885cb0e6 Mon Sep 17 00:00:00 2001
From: sbiri <sbiri@noc.ac.uk>
Date: Wed, 1 Jun 2022 11:59:36 +0000
Subject: [PATCH] add option 5 for gustiness following UA and shift C35 code to
option 6 Add all function (unified, not unified) in cs_wl_subs.py
---
AirSeaFluxCode/src/AirSeaFluxCode.py | 17 +-
AirSeaFluxCode/src/AirSeaFluxCode_dev.py | 6 +-
AirSeaFluxCode/src/cs_wl_subs.py | 305 ++++++++++++++++++++---
AirSeaFluxCode/src/flux_subs.py | 2 +-
AirSeaFluxCode/src/flux_subs_dev.py | 2 +-
5 files changed, 277 insertions(+), 55 deletions(-)
diff --git a/AirSeaFluxCode/src/AirSeaFluxCode.py b/AirSeaFluxCode/src/AirSeaFluxCode.py
index 26034b0..0f8ef48 100644
--- a/AirSeaFluxCode/src/AirSeaFluxCode.py
+++ b/AirSeaFluxCode/src/AirSeaFluxCode.py
@@ -10,14 +10,14 @@ from cs_wl_subs import *
class S88:
def _wind_iterate(self, ind):
- if self.gust[0] in range(1, 6):
+ if self.gust[0] in range(1, 7):
self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
np.power(get_gust(
self.gust[1], self.gust[2],
self.gust[3], self.theta[ind],
self.usr[ind], self.tsrv[ind],
self.grav[ind]), 2))
- if self.gust[0] in [3, 4]:
+ if self.gust[0] in [3, 4, 5]:
# self.GustFact[ind] = 1
# option to not remove GustFact
self.u10n[ind] = self.wind[ind]-self.usr[ind]/kappa*(
@@ -92,11 +92,6 @@ class S88:
self.SST[ind]), self.rho[ind], self.Rs[ind], self.Rnl[ind],
self.cp[ind], self.lv[ind], np.copy(self.tkt[ind]),
self.usr[ind], self.tsr[ind], self.qsr[ind], self.grav[ind])
-
- # self.dter[ind] = cs_C35(np.copy(
- # self.skt[ind]), self.rho[ind], self.Rs[ind], self.Rnl[ind],
- # self.cp[ind], self.lv[ind], self.usr[ind], self.tsr[ind],
- # self.qsr[ind], self.grav[ind])
elif self.skin == "ecmwf":
self.dter[ind] = cs_ecmwf(
self.rho[ind], self.Rs[ind], self.Rnl[ind], self.cp[ind],
@@ -406,6 +401,12 @@ class S88:
self.uref = self.wind-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))
+ elif self.gust[0] == 5:
+ self.u10n = self.spd-self.usr/kappa*(
+ np.log(self.h_in[0]/self.ref10)-self.psim)
+ 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))
else:
self.u10n = self.spd-self.usr/kappa/self.GFo*(
np.log(self.h_in[0]/self.ref10)-self.psim) # C.4-7
@@ -527,7 +528,7 @@ class S88:
gust = [0, 0, 0, 0]
assert np.size(gust) == 4, "gust input must be a 4x1 array"
- assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
+ assert gust[0] in range(7), "gust at position 0 must be 0 to 6"
self.gust = gust
def _class_flag(self):
diff --git a/AirSeaFluxCode/src/AirSeaFluxCode_dev.py b/AirSeaFluxCode/src/AirSeaFluxCode_dev.py
index 81af8a7..ec02db2 100644
--- a/AirSeaFluxCode/src/AirSeaFluxCode_dev.py
+++ b/AirSeaFluxCode/src/AirSeaFluxCode_dev.py
@@ -10,7 +10,7 @@ from cs_wl_subs import *
class S88:
def _wind_iterate(self, ind):
- if self.gust[0] in range(1, 6):
+ if self.gust[0] in range(1, 7):
self.wind[ind] = np.sqrt(np.power(np.copy(self.spd[ind]), 2) +
np.power(get_gust(
self.gust[1], self.gust[2],
@@ -361,8 +361,6 @@ class S88:
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
-
-
# Set the new variables (for comparison against "old")
new = np.array([np.copy(getattr(self, i)) for i in new_vars])
@@ -582,7 +580,7 @@ class S88:
gust = [0, 0, 0, 0]
assert np.size(gust) == 4, "gust input must be a 4x1 array"
- assert gust[0] in range(6), "gust at position 0 must be 0 to 5"
+ assert gust[0] in range(7), "gust at position 0 must be 0 to 6"
self.gust = gust
def _class_flag(self):
diff --git a/AirSeaFluxCode/src/cs_wl_subs.py b/AirSeaFluxCode/src/cs_wl_subs.py
index eba881e..2bbf861 100644
--- a/AirSeaFluxCode/src/cs_wl_subs.py
+++ b/AirSeaFluxCode/src/cs_wl_subs.py
@@ -1,51 +1,56 @@
import numpy as np
from util_subs import (CtoK, kappa)
+# def delta(aw, Q, usr, grav, rho, opt):
+# """
+# Compute the thickness (m) of the viscous skin layer.
-def delta(aw, Q, usr, grav, rho, opt):
- """
- Compute the thickness (m) of the viscous skin layer.
-
- Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
- eq. 8.155 p. 164
+# Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
+# eq. 8.155 p. 164
- Parameters
- ----------
- aw : float
- thermal expansion coefficient of sea-water [1/K]
- Q : float
- part of the net heat flux actually absorbed in the warm layer [W/m^2]
- usr : float
- friction velocity in the air (u*) [m/s]
- grav : float
- gravity [ms^-2]
+# Parameters
+# ----------
+# aw : float
+# thermal expansion coefficient of sea-water [1/K]
+# Q : float
+# part of the net heat flux actually absorbed in the warm layer [W/m^2]
+# usr : float
+# friction velocity in the air (u*) [m/s]
+# grav : float
+# gravity [ms^-2]
- Returns
- -------
- delta : float
- the thickness (m) of the viscous skin layer
+# Returns
+# -------
+# delta : float
+# the thickness (m) of the viscous skin layer
- """
- if opt == "C35":
- rhow, cw, visw, tcw = 1022, 4000, 1e-6, 0.6
- # second term in eq. 14
- tmp = 16*grav*Q*aw*rhow*cw*np.power(rhow*visw, 3)/(
- np.power(usr, 4)*np.power(rho, 2)*np.power(tcw, 2))
- # second term in eq. 14
- # lm = 6*np.ones(usr.shape)
- lm = np.where(Q > 0, 6/np.power(1+np.power(tmp, 3/4), 1/3), 6) # eq. 14 F96
- delta = np.minimum(lm*visw/np.sqrt(rho/rhow)/usr, 0.01) # eq. 12 F96
- delta = np.where(Q > 0, lm*visw/(np.sqrt(rho/rhow)*usr), delta)
- elif opt == "ecmwf":
- rhow, cw, visw, tcw = 1025, 4190, 1e-6, 0.6
- # u* in the water
- usr_w = np.maximum(usr, 1e-4)*np.sqrt(rho/rhow) # rhoa=1.2
- rcst_cs = 16*grav*np.power(visw, 3)/np.power(tcw, 2)#/(rhow*cw)
- # eq. 8.155 Cy46r1
- lm = 6*(1+np.maximum(Q*aw*rcst_cs/np.power(usr_w, 4), 0)**0.75)**(-1/3)
- ztmp = visw/usr_w
- delta = np.where(Q > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
- return delta
+# """
+# if opt == "C35":
+# rhow, cw, visw, tcw = 1022, 4000, 1e-6, 0.6
+# # second term in eq. 14
+# tmp = 16*grav*Q*aw*rhow*cw*np.power(rhow*visw, 3)/(
+# np.power(usr, 4)*np.power(rho, 2)*np.power(tcw, 2))
+# # second term in eq. 14
+# # lm = 6*np.ones(usr.shape)
+# lm = np.where(Q > 0, 6/np.power(1+np.power(tmp, 3/4), 1/3), 6) # eq. 14 F96
+# delta = np.minimum(lm*visw/np.sqrt(rho/rhow)/usr, 0.01) # eq. 12 F96
+# delta = np.where(Q > 0, lm*visw/(np.sqrt(rho/rhow)*usr), delta)
+# # bigc = 16*grav*cw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(
+# # rho, 2))
+# # lm = np.where(Qb > 0, 6/(1+(bigc*Qb/usr**4)**0.75)**0.333, 6)
+# # d = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+# # d = np.where(Qb > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), d)
+# elif opt == "ecmwf":
+# rhow, cw, visw, tcw = 1025, 4190, 1e-6, 0.6
+# # u* in the water
+# usr_w = np.maximum(usr, 1e-4)*np.sqrt(rho/rhow) # rhoa=1.2
+# rcst_cs = 16*grav*np.power(visw, 3)/np.power(tcw, 2)#/(rhow*cw)
+# # eq. 8.155 Cy46r1
+# # lm = 6/np.power(1+np.power(Q*aw*rcst_cs/np.power(usr_w, 4), 3/4), 1/3)
+# lm = 6*(1+np.maximum(Q*aw*rcst_cs/np.power(usr_w, 4), 0)**0.75)**(-1/3)
+# ztmp = visw/usr_w
+# delta = np.where(Q > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
+# return delta
# ---------------------------------------------------------------------
@@ -122,6 +127,162 @@ def cs(sst, d, rho, Rs, Rnl, cp, lv, usr, tsr, qsr, grav, opt):
# ---------------------------------------------------------------------
+def cs_C35(sst, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, grav):
+ """
+ Compute cool skin.
+
+ Based on COARE3.5 (Fairall et al. 1996, Edson et al. 2013)
+
+ Parameters
+ ----------
+ sst : float
+ sea surface temperature [K]
+ qsea : float
+ specific humidity over sea [g/kg]
+ rho : float
+ density of air [kg/m^3]
+ Rs : float
+ downward shortwave radiation [Wm-2]
+ Rnl : float
+ net upwelling IR radiation [Wm-2]
+ cp : float
+ specific heat of air at constant pressure [J/K/kg]
+ lv : float
+ latent heat of vaporization [J/kg]
+ delta : float
+ cool skin thickness [m]
+ usr : float
+ friction velocity [m/s]
+ tsr : float
+ star temperature [K]
+ qsr : float
+ star humidity [g/kg]
+ grav : float
+ gravity [ms^-2]
+
+ Returns
+ -------
+ dter : float
+ cool skin correction [K]
+ dqer : float
+ humidity corrction [g/kg]
+ delta : float
+ cool skin thickness [m]
+ """
+ # coded following Saunders (1967) with lambda = 6
+ if np.nanmin(sst) > 200: # if sst in Kelvin convert to Celsius
+ sst = sst-CtoK
+ # ************ cool skin constants *******
+ # density of water, specific heat capacity of water, water viscosity,
+ # thermal conductivity of water
+ rhow, cpw, visw, tcw = 1022, 4000, 1e-6, 0.6
+ for i in range(4):
+ aw = 2.1e-5*np.power(np.maximum(sst+3.2, 0), 0.79)
+ bigc = 16*grav*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(
+ rho, 2))
+ Rns = 0.945*Rs # albedo correction
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
+ Qnsol = shf+lhf+Rnl
+ fs = 0.065+11*delta-6.6e-5/delta*(1-np.exp(-delta/8.0e-4))
+ Q = Qnsol+Rns*fs
+ Qb = aw*Q+0.026*np.minimum(lhf, 0)*cpw/lv
+ xlamx = 6*np.ones(sst.shape)
+ xlamx = np.where(Qb > 0, 6/(1+(bigc*Qb/usr**4)**0.75)**0.333, 6)
+ delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
+ delta = np.where(Qb > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta)
+ dter = Q*delta/tcw
+ return dter, delta
+# ---------------------------------------------------------------------
+
+
+def delta(aw, Q, usr, grav):
+ """
+ Compute the thickness (m) of the viscous skin layer.
+
+ Based on Fairall et al., 1996 and cited in IFS Documentation Cy46r1
+ eq. 8.155 p. 164
+
+ Parameters
+ ----------
+ aw : float
+ thermal expansion coefficient of sea-water [1/K]
+ Q : float
+ part of the net heat flux actually absorbed in the warm layer [W/m^2]
+ usr : float
+ friction velocity in the air (u*) [m/s]
+ grav : float
+ gravity [ms^-2]
+
+ Returns
+ -------
+ delta : float
+ the thickness (m) of the viscous skin layer
+
+ """
+ rhow, visw, tcw = 1025, 1e-6, 0.6
+ # u* in the water
+ usr_w = np.maximum(usr, 1e-4)*np.sqrt(1.2/rhow) # rhoa=1.2
+ rcst_cs = 16*grav*np.power(visw, 3)/np.power(tcw, 2)
+ lm = 6*(1+np.maximum(Q*aw*rcst_cs/np.power(usr_w, 4), 0)**0.75)**(-1/3)
+ ztmp = visw/usr_w
+ delta = np.where(Q > 0, np.minimum(6*ztmp, 0.007), lm*ztmp)
+ return delta
+# ---------------------------------------------------------------------
+
+
+def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, grav):
+ """
+ cool skin adjustment based on IFS Documentation cy46r1
+
+ Parameters
+ ----------
+ rho : float
+ density of air [kg/m^3]
+ Rs : float
+ downward solar radiation [Wm-2]
+ Rnl : float
+ net thermal radiation [Wm-2]
+ cp : float
+ specific heat of air at constant pressure [J/K/kg]
+ lv : float
+ latent heat of vaporization [J/kg]
+ usr : float
+ friction velocity [m/s]
+ tsr : float
+ star temperature [K]
+ qsr : float
+ star humidity [g/kg]
+ sst : float
+ sea surface temperature [K]
+ grav : float
+ gravity [ms^-2]
+
+ Returns
+ -------
+ dtc : float
+ cool skin temperature correction [K]
+
+ """
+ if np.nanmin(sst) < 200: # if sst in Celsius convert to Kelvin
+ sst = sst+CtoK
+ aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
+ Rns = 0.945*Rs # (net solar radiation (albedo correction)
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
+ Qnsol = shf+lhf+Rnl # eq. 8.152
+ d = delta(aw, Qnsol, usr, grav)
+ for jc in range(4): # because implicit in terms of delta...
+ # # fraction of the solar radiation absorbed in layer delta eq. 8.153
+ # and Eq.(5) Zeng & Beljaars, 2005
+ fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8e-4))
+ Q = Qnsol+fs*Rns
+ d = delta(aw, Q, usr, grav)
+ dtc = Q*d/0.6 # (rhow*cw*kw)eq. 8.151
+ return dtc
+# ---------------------------------------------------------------------
+
+
def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, grav):
"""
Calculate warm layer correction following IFS Documentation cy46r1.
@@ -200,6 +361,68 @@ def wl_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, skt, dtc, grav):
# ----------------------------------------------------------------------
+def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, grav, Qs):
+ """
+ cool skin adjustment based on Beljaars (1997)
+ air-sea interaction in the ECMWF model
+
+ Parameters
+ ----------
+ rho : float
+ density of air [kg/m^3]
+ Rs : float
+ downward solar radiation [Wm-2]
+ Rnl : float
+ net thermal radiaion [Wm-2]
+ cp : float
+ specific heat of air at constant pressure [J/K/kg]
+ lv : float
+ latent heat of vaporization [J/kg]
+ usr : float
+ friction velocity [m/s]
+ tsr : float
+ star temperature [K]
+ qsr : float
+ star humidity [g/kg]
+ grav : float
+ gravity [ms^-2]
+ Qs : float
+ radiation balance
+
+ Returns
+ -------
+ Qs : float
+ radiation balance
+ dtc : float
+ cool skin temperature correction [K]
+
+ """
+ tcw = 0.6 # thermal conductivity of water (at 20C) [W/m/K]
+ visw = 1e-6 # kinetic viscosity of water [m^2/s]
+ rhow = 1025 # Density of sea-water [kg/m^3]
+ cpw = 4190 # specific heat capacity of water
+ aw = 3e-4 # thermal expansion coefficient [K-1]
+ Rns = 0.945*Rs # net solar radiation (albedo correction)
+ shf = rho*cp*usr*tsr
+ lhf = rho*lv*usr*qsr
+ Q = Rnl+shf+lhf+Qs
+ xt = 16*Q*grav*aw*cpw*np.power(rhow*visw, 3)/(
+ np.power(usr, 4)*np.power(rho*tcw, 2))
+ xt1 = 1+np.power(xt, 3/4)
+ # Saunders const eq. 22
+ ls = np.where(Q > 0, 6/np.power(xt1, 1/3), 6)
+ delta = np.where(Q > 0, (ls*visw)/(np.sqrt(rho/rhow)*usr),
+ np.where((ls*visw)/(np.sqrt(rho/rhow)*usr) > 0.01, 0.01,
+ (ls*visw)/(np.sqrt(rho/rhow)*usr))) # eq. 21
+ # fraction of the solar radiation absorbed in layer delta
+ fc = 0.065+11*delta-6.6e-5*(1-np.exp(-delta/0.0008))/delta
+ Qs = fc*Rns
+ Q = Rnl+shf+lhf+Qs
+ dtc = Q*delta/tcw
+ return Qs, dtc
+# ----------------------------------------------------------------------
+
+
def get_dqer(dter, sst, qsea, lv):
"""
Calculate humidity correction.
diff --git a/AirSeaFluxCode/src/flux_subs.py b/AirSeaFluxCode/src/flux_subs.py
index a5c08dc..665c7a6 100755
--- a/AirSeaFluxCode/src/flux_subs.py
+++ b/AirSeaFluxCode/src/flux_subs.py
@@ -633,7 +633,7 @@ def apply_GF(gust, spd, wind, step):
GustFact = wind*0+1
if gust[0] in [1, 2]:
GustFact = np.sqrt(wind/spd)
- elif gust[0] == 5:
+ elif gust[0] == 6:
# as in C35 matlab code
GustFact = wind/spd
elif step == "TSF":
diff --git a/AirSeaFluxCode/src/flux_subs_dev.py b/AirSeaFluxCode/src/flux_subs_dev.py
index 5609f53..cc33b3e 100644
--- a/AirSeaFluxCode/src/flux_subs_dev.py
+++ b/AirSeaFluxCode/src/flux_subs_dev.py
@@ -704,7 +704,7 @@ def apply_GF(gust, spd, wind, step):
GustFact = wind*0+1
if gust[0] in [1, 2]:
GustFact = np.sqrt(wind/spd)
- elif gust[0] == 5:
+ elif gust[0] == 6:
# as in C35 matlab code
GustFact = wind/spd
elif step == "TSF":
--
GitLab