Skip to content

GitLab

Commit 46e2a4c8 authored by sbiri's avatar sbiri
Browse files

temporary changes for the unification of cool skin

parent c0ee1bf6
Hide whitespace changes
Inline Side-by-side
Showing with 73 additions and 238 deletions
Code/cs_wl_subs.py
View file @ 46e2a4c8
import numpy as np
import numpy as np
from util_subs import (CtoK, kappa)
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))
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
return dter, delta
# ---------------------------------------------------------------------
def delta(aw, Q, usr, grav):
def delta(aw, Q, usr, grav, rho, opt):
"""
Compute the thickness (m) of the viscous skin layer.
......@@ -95,133 +26,99 @@ def delta(aw, Q, usr, grav):
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)
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
# ---------------------------------------------------------------------
# version that doesn't output dqer or import/output delta
# should be very similar to cs_ecmwf
# def cs_C35(sst, rho, Rs, Rnl, cp, lv, 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]
# 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 [ms^-1]
# tsr : float
# star temperature [K]
# qsr : float
# star humidity [g/kg]
# grav : float
# gravity [ms^-2]
# Returns
# -------
# dter : float
# cool skin correction [K]
# 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
# 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
# d = delta(aw, Qnsol, usr, grav)
# for i in range(4):
# fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/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)
# d = np.minimum(xlamx*visw/(np.sqrt(rho/rhow)*usr), 0.01)
# d = np.where(alq > 0, xlamx*visw/(np.sqrt(rho/rhow)*usr), d)
# dter = qcol*d/tcw
# return dter
# ---------------------------------------------------------------------
def cs_ecmwf(rho, Rs, Rnl, cp, lv, usr, tsr, qsr, sst, grav):
def cs(sst, d, rho, Rs, Rnl, cp, lv, usr, tsr, qsr, grav, opt):
"""
cool skin adjustment based on IFS Documentation cy46r1
Compute cool skin.
Based on COARE3.5 (Fairall et al. 1996, Edson et al. 2013)
Parameters
----------
sst : float
sea surface temperature [K]
d : float
cool skin thickness [m]
rho : float
density of air [kg/m^3]
Rs : float
downward solar radiation [Wm-2]
downward shortwave radiation [Wm-2]
Rnl : float
net thermal radiation [Wm-2]
net upwelling IR radiation [Wm-2]
cp : float
specific heat of air at constant pressure [J/K/kg]
specific heat of air at constant pressure [J/K/kg]
lv : float
latent heat of vaporization [J/kg]
latent heat of vaporization [J/kg]
usr : float
friction velocity [m/s]
friction velocity [ms^-1]
tsr : float
star temperature [K]
star temperature [K]
qsr : float
star humidity [g/kg]
sst : float
sea surface temperature [K]
star humidity [g/kg]
grav : float
gravity [ms^-2]
gravity [ms^-2]
opt : str
method to follow
Returns
-------
dtc : float
cool skin temperature correction [K]
dter : float
cool skin correction [K]
delta : float
cool skin thickness [m]
"""
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)
# 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
tcw = 0.6
Rns = 0.945*Rs # 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
Qnsol = shf+lhf+Rnl
if opt == "C35":
cpw = 4000
aw = 2.1e-5*np.power(sst+3.2, 0.79)
# d = delta(aw, Qnsol, usr, grav, rho, opt)
fs = 0.065+11*d-6.6e-5/d*(1-np.exp(-d/8.0e-4)) # eq. 17 F96
# in F96 first term in eq. 17 is 0.137 insted of 0.065
Q = Qnsol+Rns*fs
Qb = aw*Q+0.026*np.minimum(lhf, 0)*cpw/lv # eq. 8 F96
d = delta(aw, Qb, usr, grav, rho, opt)
elif opt == "ecmwf":
aw = np.maximum(1e-5, 1e-5*(sst-CtoK))
# d = delta(aw, Qnsol, usr, grav, rho, opt)
for jc in range(4):
# 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)) # eq. 8.153 Cy46r1
Q = Qnsol+Rns*fs
d = delta(aw, Q, usr, grav, rho, opt)
dter = Q*d/tcw # eq. 4 F96
return dter, d
# ---------------------------------------------------------------------
......@@ -303,68 +200,6 @@ 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.
......@@ -391,4 +226,4 @@ def get_dqer(dter, sst, qsea, lv):
wetc = 0.622*lv*qsea/(287.1*np.power(sst, 2))
dqer = wetc*dter
return dqer
# ----------------------------------------------------------------------
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