Skip to content

GitLab

"backup_core_local" did not exist on "e29cc4dfc909c281a77ae01116a9624efb066ea7"
Commit a5d491a1 authored by Richard Cornes's avatar Richard Cornes
Browse files

Moved Cool skin / warm layer functions to separate file

parent b523c016
Hide whitespace changes
Inline Side-by-side
Showing with 299 additions and 298 deletions
Code/AirSeaFluxCode.py
View file @ a5d491a1
......@@ -5,6 +5,7 @@ import logging
from hum_subs import (get_hum, gamma)
from util_subs import *
from flux_subs import *
from cs_wl_subs import *
class S88:
......
Code/cs_wl_subs.py 0 → 100644
View file @ a5d491a1
import numpy as np
from util_subs import (CtoK, kappa, gc)
def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
"""
Computes cool skin
following 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]
lat : float
latitude [degE]
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
g = gc(lat, None)
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*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(rho, 2))
wetc = 0.622*lv*qsea/(287.1*np.power(sst+273.16, 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))
qcol = Qnsol+Rns*fs
alq = aw*qcol+0.026*np.minimum(lhf, 0)*cpw/lv
xlamx = 6*np.ones(sst.shape)
xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
delta = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta)
dter = qcol*delta/tcw
dqer = wetc*dter
return dter, dqer, delta
# ---------------------------------------------------------------------
def delta(aw, Q, usr, lat):
"""
Computes 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]
lat : float
latitude [degE]
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*gc(lat)*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, lat):
"""
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]
lat : float
latitude [degE]
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, lat)
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, lat)
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, lat):
"""
warm layer correction following IFS Documentation cy46r1
and aerobulk (Brodeau et al., 2016)
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
bulk sst [K]
skt : float
skin sst from previous step [K]
dtc : float
cool skin correction [K]
lat : float
latitude [degE]
Returns
-------
dtwl : float
warm layer correction [K]
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
rhow, cpw, visw, rd0 = 1025, 4190, 1e-6, 3
Rns = 0.945*Rs
## Previous value of dT / warm-layer, adapted to depth:
aw = 2.1e-5*np.power(np.maximum(sst-CtoK+3.2, 0), 0.79) # thermal expansion coefficient of sea-water (SST accurate enough#)
# *** Rd = Fraction of solar radiation absorbed in warm layer (-)
a1, a2, a3 = 0.28, 0.27, 0.45
b1, b2, b3 = -71.5, -2.8, -0.06 # [m-1]
Rd = 1-(a1*np.exp(b1*rd0)+a2*np.exp(b2*rd0)+a3*np.exp(b3*rd0))
shf = rho*cp*usr*tsr
lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl
usrw = np.maximum(usr, 1e-4 )*np.sqrt(1.2/rhow) # u* in the water
zc3 = rd0*kappa*gc(lat)/np.power(1.2/rhow, 3/2)
zc4 = (0.3+1)*kappa/rd0
zc5 = (0.3+1)/(0.3*rd0)
for jwl in range(10): # itteration to solve implicitely equation for warm layer
dsst = skt-sst-dtc
## Buoyancy flux and stability parameter (zdl = -z/L) in water
ZSRD = (Qnsol+Rns*Rd)/(rhow*cpw)
ztmp = np.maximum(dsst, 0)
zdl = np.where(ZSRD > 0, 2*(np.power(usrw, 2) *
np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))+ZSRD,
np.power(usrw, 2)*np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))
usr = np.maximum(usr, 1e-4 )
zdL = zc3*aw*zdl/np.power(usr, 3)
## Stability function Phi_t(-z/L) (zdL is -z/L) :
zphi = np.where(zdL > 0, (1+(5*zdL+4*np.power(zdL, 2)) /
(1+3*zdL+0.25*np.power(zdL, 2)) +
2/np.sqrt(1-16*(-np.abs(zdL)))),
1/np.sqrt(1-16*(-np.abs(zdL))))
zz = zc4*(usrw)/zphi
zz = np.maximum(np.abs(zz) , 1e-4)*np.sign(zz)
dtwl = np.maximum(0, (zc5*ZSRD)/zz)
return dtwl
#----------------------------------------------------------------------
def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, lat, 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]
lat : float
latitude [degE]
Qs : float
radiation balance
Returns
-------
Qs : float
radiation balance
dtc : float
cool skin temperature correction [K]
"""
g = gc(lat, None)
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*g*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
#----------------------------------------------------------------------
Code/flux_subs.py
View file @ a5d491a1
import numpy as np
from util_subs import (CtoK, kappa, gc, visc_air)
from util_subs import (kappa, gc, visc_air)
# ---------------------------------------------------------------------
......@@ -258,7 +258,6 @@ def get_stabco(meth):
alpha, beta, gamma = 20, 0.25, 5
else:
raise ValueError("unknown method stabco: "+meth)
coeffs = np.zeros(3)
coeffs[0] = alpha
coeffs[1] = beta
......@@ -552,302 +551,6 @@ def psim_stab(zol, meth):
# ---------------------------------------------------------------------
def cs_C35(sst, qsea, rho, Rs, Rnl, cp, lv, delta, usr, tsr, qsr, lat):
"""
Computes cool skin
following 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]
lat : float
latitude [degE]
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
g = gc(lat, None)
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*g*cpw*np.power(rhow*visw, 3)/(np.power(tcw, 2)*np.power(rho, 2))
wetc = 0.622*lv*qsea/(287.1*np.power(sst+273.16, 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))
qcol = Qnsol+Rns*fs
alq = aw*qcol+0.026*np.minimum(lhf, 0)*cpw/lv
xlamx = 6*np.ones(sst.shape)
xlamx = np.where(alq > 0, 6/(1+(bigc*alq/usr**4)**0.75)**0.333, 6)
delta = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
delta = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), delta)
dter = qcol*delta/tcw
dqer = wetc*dter
return dter, dqer, delta
# ---------------------------------------------------------------------
def delta(aw, Q, usr, lat):
"""
Computes 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]
lat : float
latitude [degE]
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*gc(lat)*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, lat):
"""
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]
lat : float
latitude [degE]
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, lat)
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, lat)
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, lat):
"""
warm layer correction following IFS Documentation cy46r1
and aerobulk (Brodeau et al., 2016)
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
bulk sst [K]
skt : float
skin sst from previous step [K]
dtc : float
cool skin correction [K]
lat : float
latitude [degE]
Returns
-------
dtwl : float
warm layer correction [K]
"""
if (np.nanmin(sst) < 200): # if sst in Celsius convert to Kelvin
sst = sst+CtoK
rhow, cpw, visw, rd0 = 1025, 4190, 1e-6, 3
Rns = 0.945*Rs
## Previous value of dT / warm-layer, adapted to depth:
aw = 2.1e-5*np.power(np.maximum(sst-CtoK+3.2, 0), 0.79) # thermal expansion coefficient of sea-water (SST accurate enough#)
# *** Rd = Fraction of solar radiation absorbed in warm layer (-)
a1, a2, a3 = 0.28, 0.27, 0.45
b1, b2, b3 = -71.5, -2.8, -0.06 # [m-1]
Rd = 1-(a1*np.exp(b1*rd0)+a2*np.exp(b2*rd0)+a3*np.exp(b3*rd0))
shf = rho*cp*usr*tsr
lhf = rho*lv*usr*qsr
Qnsol = shf+lhf+Rnl
usrw = np.maximum(usr, 1e-4 )*np.sqrt(1.2/rhow) # u* in the water
zc3 = rd0*kappa*gc(lat)/np.power(1.2/rhow, 3/2)
zc4 = (0.3+1)*kappa/rd0
zc5 = (0.3+1)/(0.3*rd0)
for jwl in range(10): # itteration to solve implicitely equation for warm layer
dsst = skt-sst-dtc
## Buoyancy flux and stability parameter (zdl = -z/L) in water
ZSRD = (Qnsol+Rns*Rd)/(rhow*cpw)
ztmp = np.maximum(dsst, 0)
zdl = np.where(ZSRD > 0, 2*(np.power(usrw, 2) *
np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))+ZSRD,
np.power(usrw, 2)*np.sqrt(ztmp/(5*rd0*gc(lat)*aw/visw)))
usr = np.maximum(usr, 1e-4 )
zdL = zc3*aw*zdl/np.power(usr, 3)
## Stability function Phi_t(-z/L) (zdL is -z/L) :
zphi = np.where(zdL > 0, (1+(5*zdL+4*np.power(zdL, 2)) /
(1+3*zdL+0.25*np.power(zdL, 2)) +
2/np.sqrt(1-16*(-np.abs(zdL)))),
1/np.sqrt(1-16*(-np.abs(zdL))))
zz = zc4*(usrw)/zphi
zz = np.maximum(np.abs(zz) , 1e-4)*np.sign(zz)
dtwl = np.maximum(0, (zc5*ZSRD)/zz)
return dtwl
#----------------------------------------------------------------------
def cs_Beljaars(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, lat, 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]
lat : float
latitude [degE]
Qs : float
radiation balance
Returns
-------
Qs : float
radiation balance
dtc : float
cool skin temperature correction [K]
"""
g = gc(lat, None)
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*g*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_gust(beta, Ta, usr, tsrv, zi, lat):
"""
Computes gustiness
......
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