MODULE drag3
USE param
implicit none
type, public :: meteoDataType
real :: h ! constant flux layer height [m]
real :: U ! abs(wind speed) at 'h' [m/s]
real :: dT ! difference between potential temperature at 'h' and at surface [K]
real :: Tsemi ! semi-sum of potential temperature at 'h' and at surface [K]
real :: dQ ! difference between humidity at 'h' and at surface [g/g]
real :: z0_m ! surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
type, public :: meteoDataVecType
real, allocatable :: h(:) ! constant flux layer height [m]
real, allocatable :: U(:) ! abs(wind speed) at 'h' [m/s]
real, allocatable :: dT(:) ! difference between potential temperature at 'h' and at surface [K]
real, allocatable :: Tsemi(:) ! semi-sum of potential temperature at 'h' and at surface [K]
real, allocatable :: dQ(:) ! difference between humidity at 'h' and at surface [g/g]
real, allocatable :: z0_m(:) ! surface aerodynamic roughness (should be < 0 for water bodies surface)
end type
type, public :: sfxDataType
real :: zeta ! = z/L [n/d]
real :: Rib ! bulk Richardson number [n/d]
real :: Re ! Reynolds number [n/d]
real :: lnzuzt
real :: z0_m ! aerodynamic roughness [m]
real :: z0_t ! thermal roughness [m]
real :: Rib_conv_lim ! Ri-bulk convection critical value [n/d]
real :: cm
real :: ch
real :: ct
real :: ckt
end type
type, public :: numericsType
integer :: maxiters_convection = 10 ! maximum (actual) number of iterations in convection
integer :: maxiters_charnock = 10 ! maximum (actual) number of iterations in charnock roughness
end type
type, public:: data_lutyp
integer, public :: lu_indx=1
end type
contains
subroutine surf_flux_vec(meteo, &
masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
numerics, lu1,numst)
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
real, dimension (numst) :: masout_zl
real, dimension (numst) :: masout_ri
real, dimension (numst) :: masout_re
real, dimension (numst) :: masout_lnzuzt
real, dimension (numst) :: masout_zu
real, dimension (numst) :: masout_ztout
real, dimension (numst) :: masout_rith
real, dimension (numst) :: masout_cm
real, dimension (numst) :: masout_ch
real, dimension (numst) :: masout_ct
real, dimension (numst) :: masout_ckt
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
type (data_lutyp) lu1
integer i,numst
do i = 1, numst
meteo_cell = meteoDataType(&
h = meteo%h(i), &
U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
z0_m = meteo%z0_m(i))
call surf_flux(meteo_cell, sfx_cell, numerics, lu1)
masout_zl(i)=sfx_cell%zeta
masout_ri(i)=sfx_cell%Rib
masout_re(i)=sfx_cell%Re
masout_lnzuzt(i)=sfx_cell%lnzuzt
masout_zu(i)=sfx_cell%z0_m
masout_ztout(i)=sfx_cell%z0_t
masout_rith(i)=sfx_cell%Rib_conv_lim
masout_cm(i)=sfx_cell%cm
masout_ch(i)=sfx_cell%ch
masout_ct(i)=sfx_cell%ct
masout_ckt(i)=sfx_cell%ckt
enddo
END SUBROUTINE surf_flux_vec
function f_m_conv(zeta)
real f_m_conv
real zeta
f_m_conv = (1.0 - alpha_m * zeta)**0.25
end function f_m_conv
function f_h_conv(zeta)
real f_h_conv
real zeta
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
! stability functions (neutral)
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
! ----------------------------------------------------------------------------
zeta = 0.0
psi_m = log(h0_m)
psi_h = log(h0_t) / Pr_t_0_inv
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
! --------------------------------------------------------------------------------
! stability functions (stable)
! --------------------------------------------------------------------------------
subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: Rib_coeff
real :: psi0_m, psi0_h
real :: phi
real :: c
! ----------------------------------------------------------------------------
psi0_m = alog(h0_m)
psi0_h = B / psi0_m
Rib_coeff = beta_m * Rib
c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
(2.0 * beta_m * (1.0 - Rib_coeff))
phi = beta_m * zeta
psi_m = psi0_m + phi
psi_h = (psi0_m + B) / Pr_t_0_inv + phi
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection-semistrong)
! --------------------------------------------------------------------------------
subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, x1, y0, y1
integer :: i
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t)
if (abs(B) < 1.0e-10) psi_h = psi_m
zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
y1 = (1.0 - alpha_m * zeta)**0.25e0
x1 = sqrt(1.0 - alpha_h_fix * zeta)
y0 = (1.0 - alpha_m * zeta0_m)**0.25e0
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = max(y0, 1.000001e0)
x0 = max(x0, 1.000001e0)
psi_m = log((y1 - 1.0e0)*(y0 + 1.0e0)/((y1 + 1.0e0)*(y0 - 1.0e0))) + 2.0e0*(atan(y1) - atan(y0))
psi_h = log((x1 - 1.0e0)*(x0 + 1.0e0)/((x1 + 1.0e0)*(x0 - 1.0e0))) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! stability functions (convection)
! --------------------------------------------------------------------------------
subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, x10, y10, &
h0_m, h0_t, B, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h, zeta
real, intent(in) :: Rib
real, intent(in) :: zeta_conv_lim
real, intent(in) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
integer, intent(in) :: maxiters
real :: zeta0_m, zeta0_h
real :: x0, y0
real :: p1, p0
real :: a1, b1, c, f
integer :: i
! ----------------------------------------------------------------------------
p1 = 2.0 * atan(y10) + log((y10 - 1.0) / (y10 + 1.0))
p0 = log((x10 - 1.0) / (x10 + 1.0))
zeta = zeta_conv_lim
do i = 1, maxiters + 1
zeta0_m = zeta / h0_m
zeta0_h = zeta / h0_t
if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
a1 = (zeta_conv_lim / zeta)**(1.0 / 3.0)
x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
y0 = (1.0 - alpha_m * zeta0_m)**0.25
c = log((x0 + 1.0)/(x0 - 1.0))
b1 = -2.0*atan(y0) + log((y0 + 1.0)/(y0 - 1.0))
f = 3.0 * (1.0 - a1)
psi_m = f / y10 + p1 + b1
psi_h = (f / x10 + p0 + c) / Pr_t_0_inv
if (i == maxiters + 1) exit
zeta = Rib * psi_m**2 / psi_h
end do
end subroutine
! --------------------------------------------------------------------------------
! define convection limit
! --------------------------------------------------------------------------------
subroutine get_convection_lim(zeta_lim, Rib_lim, x10, y10, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim, Rib_lim
real, intent(out) :: x10, y10
real, intent(in) :: h0_m, h0_t
real, intent(in) :: B
real :: an4, an5
real :: psi_m, psi_h
! ----------------------------------------------------------------------------
an5=(Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*(alpha_h-alpha_m))) &
/ (2.0E0*alpha_h**2)
y10 = f_m_conv(zeta_lim)
x10 = f_h_conv(zeta_lim)
! ......definition of r-prim......
an4 = zeta_lim / h0_m
an5 = zeta_lim / h0_t
if (abs(B) < 1.0e-10) an5=an4
an5 = sqrt(1.0 - alpha_h_fix * an5)
an4 = (1.0 - alpha_m * an4)**0.25
psi_m = 2.0 * (atan(y10) - atan(an4)) + alog((y10 - 1.0) * (an4 + 1.0)/((y10 + 1.0) * (an4 - 1.0)))
psi_h = alog((x10 - 1.0) * (an5 + 1.0)/((x10 + 1.0) * (an5 - 1.0))) / Pr_t_0_inv
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! --------------------------------------------------------------------------------
! charnock roughness
! --------------------------------------------------------------------------------
subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_m, u_dyn0
real, intent(in) :: U, h
integer, intent(in) :: maxiters
real :: U10m
real :: a1, c1, y1, c1min, f
integer :: i, j
! ----------------------------------------------------------------------------
U10m = U
a1 = 0.0e0
y1 = 25.0e0
c1min = log(h_charnock) / kappa
do i = 1, maxiters
f = c1_charnock - 2.0 * log(U10m)
do j = 1, maxiters
c1 = (f + 2.0 * log(y1)) / kappa
! looks like the check should use U10 instead of U
! but note that a1 should be set = 0 in cycle beforehand
if(U <= 8.0e0) a1 = log(1.0 + c2_charnock * ((y1 / U10m)**3)) / kappa
c1 = max(c1 - a1, c1min)
y1 = c1
end do
z0_m = h_charnock * exp(-c1 * kappa)
z0_m = max(z0_m,0.000015e0)
U10m = U*alog(h_charnock/z0_m)/(alog(h/z0_m))
end do
! --- define dynamic velocity in neutral conditions
u_dyn0 = U10m / c1
end subroutine
! --------------------------------------------------------------------------------
SUBROUTINE surf_flux(meteo, sfx, numerics, lu)
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
type (data_lutyp) lu
real h ! surface flux layer height [m]
real z0_m, z0_t ! aerodynamic & thermal roughness [m]
real B ! = ln(z0_m / z0_t) [n/d]
real h0_m, h0_t ! = h / z0_m, h / z0_h [n/d]
real Re ! roughness Reynolds number = u_dyn0 * z0 / nu [n/d]
real u_dyn0 ! dynamic velocity in neutral conditions [m/s]
real zeta ! = z/L [n/d]
real zeta_conv_lim ! z/L critical value for matching free convection limit [n/d]
real Rib ! bulk Richardson number
real Rib_conv_lim ! Ri-bulk critical value for matching free convection limit [n/d]
real psi_m, psi_h ! universal functions [momentum] & [thermal]
real U ! wind velocity at h [m/s]
real Tsemi !
integer lu_indx
! output ... !
real an4, o, am
real c4, c0, an0
! ...
! local unnamed !
real x10, y10
real al
real q4, t4
! .....
! just local variables
real f1, a1
! ....
integer is_ocean
integer i, j
Tsemi = meteo%Tsemi
U = meteo%U
t4 = meteo%dT
q4 = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
if(z0_m < 0.0) then
!> @brief Test - definition z0 of sea surface
!> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
is_ocean = 1
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
else
! ......parameters from viscosity sublayer......
is_ocean = 0
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
Re = u_dyn0 * z0_m / nu_air
if(Re <= Re_rough_min) then
B = B1_rough * alog(B3_rough * Re) + B2_rough
else
! B4 takes into account Re value at z' ~ O(10) z0
B = B4_rough * (Re**B2_rough)
end if
! --- apply max restriction based on surface type
if (is_ocean == 1) then
B = min(B, B_max_ocean)
else
B = min(B, B_max_land)
end if
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
! --- define relative height [thermal]
h0_t = h / z0_t
! ......humidity stratification and ri-number......
al = g / Tsemi
Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
h0_m, h0_t, B)
if (Rib > 0.0e0) then
! ......stable stratification......
!write (*,*) 'stable'
Rib = min(Rib, Rib_max)
call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
f1 = beta_m * zeta
am = 1.0 + f1
o = 1.0/Pr_t_0_inv + f1
else if (Rib < Rib_conv_lim) then
! ......strong instability.....
!write (*,*) 'instability'
call get_psi_convection(psi_m, psi_h, zeta, Rib, &
zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)
a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
am = a1 / y10
o = a1 / (Pr_t_0_inv * x10)
else if (Rib > -0.001e0) then
! ......nearly neutral......
write (*,*) 'neutral'
call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
am = 1.0e0
o = 1.0e0 / Pr_t_0_inv
else
! ......week and semistrong instability......
!write (*,*) 'semistrong'
call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
am = (1.0 - alpha_m * zeta)**(-0.25)
o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
end if
! ......computation of cu,co,k(h),alft
c4=kappa/psi_m
c0=kappa/psi_h
an4=kappa*c4*U*h/am
an0=am/o
! ......exit......
sfx%zeta = zeta
sfx%Rib = Rib
sfx%Re = Re
sfx%lnzuzt=B
sfx%z0_m = z0_m
sfx%z0_t = z0_t
sfx%Rib_conv_lim = Rib_conv_lim
sfx%cm=c4
sfx%ch=c0
sfx%ct=an4
sfx%ckt=an0
return
END SUBROUTINE surf_flux
END MODULE drag3