Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • inmcm-mirror/sfx
1 result
Show changes
Commits on Source (49)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 37100 additions and 34378 deletions
CMakeLists.txt
View file @ 28e6903d
......@@ -83,15 +83,20 @@ set(SOURCES_F
srcF/sfx_log_param.f90
srcF/sfx_run.f90
srcF/sfx_phys_const.f90
srcF/sfx_z0m_all_surface.f90
srcF/sfx_z0t_all_surface.f90
srcF/sfx_surface.f90
srcF/sfx_thermal_roughness.f90
srcF/sfx_most.f90
srcF/sfx_most_param.f90
srcF/sfx_most_snow.f90
srcF/sfx_most_snow_param.f90
srcF/sfx_sheba.f90
srcF/sfx_sheba_noniterative.f90
srcF/sfx_sheba_param.f90
srcF/sfx_sheba_noniterative.f90
srcF/sfx_sheba_noit_param.f90
srcF/sfx_sheba_coare.f90
srcF/sfx_sheba_coare_param.f90
srcF/sfx_fc_wrapper.F90
srcF/sfx_api_inmcm.f90
srcF/sfx_api_term.f90
......
api-examples/sfx_inmcm_ex.f90
View file @ 28e6903d
......@@ -30,7 +30,7 @@ program sfx_inmcm_ex
AR1(6) = 1.0e-3 ! surface roughness parameter
! --- converting AR input to SFX format -> meteo cell
call inmcm_to_sfx_in_cell(meteo_cell, AR1)
call inmcm_to_sfx_in_cell(meteo_cell, AR1, IVEG_sfx, depth_inm, lai_inm)
! --- calculating fluxes
call get_surface_fluxes_esm(sfx_cell, meteo_cell, numerics)
! --- converting SFX cell output to AR format
......
data/MOSAiC.txt
View file @ 28e6903d
This diff is collapsed.
diag/pbldia_new_sfx.f90
View file @ 28e6903d
module pbldia_new_sfx
implicit none
private
public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
pbldia_new_most,pbldia_new_esm
real,parameter,private::kappa=0.4,R=287.,appa=0.287
contains
! in this module interpolation is performed using integration of universal functions from z1 to z2,
! where z1 is measurement or lowest model level height and z2 is required height
! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx)
!C* DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND),
!C* OF POTENTIAL TEMPERATURE DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE),
!C* OF SPECIFIC HUMIDITY DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE)
!C* INPUT DATA USED:
!C* ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE)
!C* AR2(1) - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH)
!C* HEIGHT OF CONSTANT FLUX LAYER TOP
!C* AR2(8) - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM
!C* AR2(9) - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY
!C* ARRAY ARDIN
!C* ARDIN(1) - MODULE OF WIND VELOCITY AT THE TOP OF CONSTANT FLUX LAYER
!C* ARDIN(2) - POTENTIAL TEMPERATURE DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(3) - SPECIFIC HUMIDITY DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(4) - DIMENSIONAL HEIGHT OF CONSTANT FLUX LAYER TOP
!C* ARDIN(5) - = HWIND
!C* ARDIN(6) - = HTEMP
!C* OUTPUT DATA (ARRAY ARDOUT):
!C* ARDOUT(1) - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT
!C* ARDOUT(2) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
!C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
!C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
!C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
use sfx_sheba, only: &
get_psi_sheba => get_psi
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) zeta=zetalim
L = HIN/zeta
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_sheba(psi_m,psi_h,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_sheba(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_sheba
subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT)
use sfx_sheba_noniterative, only: &
get_psi_stable_sheba => get_psi_stable
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) zeta=zetalim
L = HIN/zeta
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_stable_sheba(psi_m,psi_h,HTEMP/L,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else ! unstable stratification (functions from ESM)
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
endif
else
! neutral stratification (z/L=0)
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_sheba_noit
subroutine pbldia_new_most(AR2,ARDIN,ARDOUT)
use sfx_most, only: &
get_psi_most => get_psi
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) zeta=zetalim
L = HIN/zeta
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_most(psi_m,psi_h,HWIND/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_most(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_most
subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT)
real,intent(in):: AR2(11),ARDIN(6)
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 2. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) zeta=zetalim
L = HIN/zeta
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_esm
subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
use sfx_esm_param
real, intent(out) :: psi_m, psi_h
real, intent(in) :: z1,z2,L
real an5,d1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2
an5 = (Pr_t_inf_inv/Pr_t_0_inv)**4
d1=(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)
zeta1 = z1 / L
zeta2 = z2 / L
if(zeta2 .ge. 0.) then
psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)
else
gmz2 = sqrt(sqrt(1.E0 - alpha_m * zeta2))
gmz1 = sqrt(sqrt(1.E0 - alpha_m * zeta1))
gmst = sqrt(sqrt(1.E0 - alpha_m * D1))
gtz2 = sqrt(1.E0 - alpha_h * zeta2)
gtz1 = sqrt(1.E0 - alpha_h * zeta1)
gtst = sqrt(1.E0 - alpha_h * D1)
if(zeta2 .ge. d1) then
psi_m = fim(gmz2) - fim(gmz1)
psi_h = fit(gtz2) - fit(gtz1)
else
psi_m = (3.E0/gmst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fim(gmst) - fim(gmz1)
psi_h = (3.E0/gtst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fit(gtst) - fit(gtz1)
endif
endif
end subroutine get_psi_esm1
! functions fim,fit were copied from the original pbldia
REAL FUNCTION FIM(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X)
RETURN
END function
REAL FUNCTION FIT(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIT = ALOG((X-1.E0)/(X+1.E0))
RETURN
END function
end module pbldia_new_sfx
\ No newline at end of file
implicit none
private
public :: pbldia_new_sheba,pbldia_new_sheba_noit,&
pbldia_new_most,pbldia_new_esm,pbldia_new_log,pbldia_new_sheba_coare
real,parameter,private::kappa=0.4,R=287.,appa=0.287
contains
! in this module interpolation is performed using integration of universal functions from z1 to z2,
! where z1 is measurement or lowest model level height and z2 is required height
! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx)
!C* DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND),
!C* OF POTENTIAL TEMPERATURE DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE),
!C* OF SPECIFIC HUMIDITY DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE)
!C* INPUT DATA USED:
!C* ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE)
!C* AR2(1) - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH)
!C* HEIGHT OF CONSTANT FLUX LAYER TOP
!C* AR2(8) - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM
!C* AR2(9) - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY
!C* ARRAY ARDIN
!C* ARDIN(1) - MODULE OF WIND VELOCITY AT THE TOP OF CONSTANT FLUX LAYER
!C* ARDIN(2) - POTENTIAL TEMPERATURE DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(3) - SPECIFIC HUMIDITY DEFICIT IN CONSTANT FLUX LAYER
!C* ARDIN(4) - DIMENSIONAL HEIGHT OF CONSTANT FLUX LAYER TOP
!C* ARDIN(5) - = HWIND
!C* ARDIN(6) - = HTEMP
!C* OUTPUT DATA (ARRAY ARDOUT):
!C* ARDOUT(1) - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT
!C* ARDOUT(2) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface
!C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface
!C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY
!C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS
subroutine pbldia_new_log(AR2,ARDIN,ARDOUT,lat)
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real hwind,htemp,ustar,tstar,qstar
real dteta,dq,wind,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
WIND = (USTAR/kappa)*ALOG(hwind/HIN)
DTETA = (TSTAR/kappa)*ALOG(htemp/HIN)
DQ = (QSTAR/kappa)*ALOG(htemp/HIN)
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_log
subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT,lat)
use sfx_sheba, only: get_psi_sheba => get_psi
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for
! stable SL for wind
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
real Rib,ksi,wg,f,L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) L_lim = HIN/zetalim
L = HIN/zeta
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L_lim)
call get_psi_sheba(psi_m,psi_h,HWIND/L_lim)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_sheba(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_sheba
subroutine pbldia_new_sheba_coare(AR2,ARDIN,ARDOUT,lat)
use sfx_sheba_coare, only: &
get_psi_coare => get_psi_a, &
get_psi_stable_sheba => get_psi_stable
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta,L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) L_lim = HIN/zetalim
L = HIN/zeta
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
else
call get_psi_coare(psi_m_hs,psi_h_hs,HIN/L,HIN/L)
call get_psi_coare(psi_m,psi_h,HWIND/L,HTEMP/L)
endif
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_sheba_coare
subroutine pbldia_new_sheba_noit(AR2,ARDIN,ARDOUT,lat)
use sfx_sheba_noniterative, only: &
get_psi_stable_sheba => get_psi_stable
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
real Rib,wg,ksi,f, L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) L_lim=HIN/zetalim
L = HIN/zeta
if(zeta.gt.0)then ! stable stratification
call get_psi_stable_sheba(psi_m_hs,psi_h_hs,HIN/L_lim,HIN/L)
call get_psi_stable_sheba(psi_m,psi_h,HWIND/L_lim,HTEMP/L)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else ! unstable stratification (functions from ESM)
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
endif
else
! neutral stratification (z/L=0)
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_sheba_noit
subroutine pbldia_new_most(AR2,ARDIN,ARDOUT,lat)
use sfx_most, only: &
get_psi_most => get_psi
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L,HIN,zeta
real Rib,ksi,wg,f,L_lim
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
Rib = AR2(2)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) L_lim=HIN/zetalim
L = HIN/zeta
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L_lim)
call get_psi_most(psi_m,psi_h,HWIND/L_lim)
UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs)
call get_psi_most(psi_m_hs,psi_h_hs,HIN/L)
call get_psi_most(psi_m,psi_h,HTEMP/L)
UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs)
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_most
subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT,lat)
real,intent(in):: AR2(11),ARDIN(6),lat
real,intent(out)::ARDOUT(3)
real,parameter::zetalim = 1. !maximum value of z/L for stable SL
real psi_m,psi_h,psi_m_hs,psi_h_hs
real hwind,htemp,ustar,tstar,qstar,&
& ufwind,uftemp,dteta,dq,wind,L_lim,L,HIN,zeta
HWIND = ARDIN(5)
HTEMP = ARDIN(6)
HIN = ardin(4)
zeta = AR2(1)
USTAR = AR2(8) * ARDIN(1)
TSTAR = AR2(9) * ARDIN(2)
QSTAR = AR2(9) * ARDIN(3)
if(zeta.ne.0)then
if(zeta.gt.zetalim) L_lim=HIN/zetalim
L = HIN/zeta
call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L_lim)
UFWIND = psi_m
call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L)
UFTEMP = psi_h
else
UFWIND = ALOG(HWIND/HIN)
UFTEMP = ALOG(HTEMP/HIN)
endif
WIND = (USTAR/kappa) * UFWIND
DTETA = (TSTAR/kappa) * UFTEMP
DQ = (QSTAR/kappa) * UFTEMP
ARDOUT(1) = ARDIN(1)+WIND
ARDOUT(2) = DTETA+ardin(2)
ARDOUT(3) = DQ+ardin(3)
return
end subroutine pbldia_new_esm
subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L)
use sfx_esm_param
real, intent(out) :: psi_m, psi_h
real, intent(in) :: z1,z2,L
real an5,d1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2
an5 = (Pr_t_inf_inv/Pr_t_0_inv)**4
d1=(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)
zeta1 = z1 / L
zeta2 = z2 / L
if(zeta2 .ge. 0.) then
psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1)
psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1)
else
gmz2 = sqrt(sqrt(1.E0 - alpha_m * zeta2))
gmz1 = sqrt(sqrt(1.E0 - alpha_m * zeta1))
gmst = sqrt(sqrt(1.E0 - alpha_m * D1))
gtz2 = sqrt(1.E0 - alpha_h * zeta2)
gtz1 = sqrt(1.E0 - alpha_h * zeta1)
gtst = sqrt(1.E0 - alpha_h * D1)
if(zeta2 .ge. d1) then
psi_m = fim(gmz2) - fim(gmz1)
psi_h = fit(gtz2) - fit(gtz1)
else
psi_m = (3.E0/gmst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fim(gmst) - fim(gmz1)
psi_h = (3.E0/gtst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fit(gtst) - fit(gtz1)
endif
endif
end subroutine get_psi_esm1
! functions fim,fit were copied from the original pbldia
REAL FUNCTION FIM(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X)
RETURN
END function
REAL FUNCTION FIT(XX)
IMPLICIT NONE
REAL X, XX
X = AMAX1(XX,1.000001E0)
FIT = ALOG((X-1.E0)/(X+1.E0))
RETURN
END function
end module pbldia_new_sfx
\ No newline at end of file
makefile
View file @ 28e6903d
......@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_surface.o sfx_log_param.o sfx_log.o sfx_most_param.o sfx_most.o sfx_sheba_param.o sfx_sheba.o sfx_esm_param.o sfx_esm.o sfx_main.o
OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_io.o sfx_data.o sfx_z0t_all_surface.o sfx_z0m_all_surface.o sfx_z0m_all.o sfx_z0t_all.o sfx_surface.o sfx_config.o sfx_log_param.o sfx_log.o sfx_most_param.o sfx_most.o sfx_most_snow_param.o sfx_most_snow.o sfx_sheba_param.o sfx_sheba.o sfx_esm_param.o sfx_esm.o sfx_sheba_noit_param.o sfx_sheba_noniterative.o sfx_sheba_coare_param.o sfx_sheba_coare.o sfx_run.o sfx_main.o
OBJ_F =
OBJ = $(OBJ_F90) $(OBJ_F)
......
srcF/module_z0t_lake.f90 deleted 100644 → 0
View file @ b5be41e8
module module_z0t_lake
!< @brief surface thermal roughness parameterizations for ocean
implicit none
public :: get_thermal_roughness_kl
public :: get_thermal_roughness_ca
public :: get_thermal_roughness_zm
public :: get_thermal_roughness_br
public :: get_thermal_roughness_re
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
real, parameter, private :: Pr_m = 0.71 !< molecular Prandtl number (air) [n/d]
!< Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3
!< roughness model coeff. [n/d]
!< --- transitional mode
!< B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m
!< --- fully rough mode (Re > Re_rough_min)
!< B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
real, parameter :: B_max_lake = 8.0
contains
! thermal roughness definition by Kazakov, Lykosov
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_kl(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
! ----------------------------------------------------------------------------
!--- define B = log(z0_m / z0_t)
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
B = min(B, B_max_lake)
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_ca(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-3.8 !4-Cahill et al.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition z0_t = C*z0_m
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_zm(z0_t, B, &
z0_m, Czm)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Czm !< proportionality coefficient
z0_t =Czm*z0_m
B=log(z0_m / z0_t)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Brutsaert W., 2003.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_br(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-2.0 !Brutsaert
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! thermal roughness definition by Repina, 2023.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_re(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=alog(-0.56*(4.0*(Re)**(0.5)-3.4)) !Repina, 2023
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
end module module_z0t_lake
srcF/module_z0t_land.f90 deleted 100644 → 0
View file @ b5be41e8
module module_z0t
!< @brief surface thermal roughness parameterizations for land
implicit none
public :: get_thermal_roughness_kl
public :: get_thermal_roughness_cz
public :: get_thermal_roughness_zi
public :: get_thermal_roughness_ca
public :: get_thermal_roughness_zm
public :: get_thermal_roughness_ot
public :: get_thermal_roughness_du
public :: get_thermal_roughness_mix
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
real, parameter, private :: Pr_m = 0.71 !< molecular Prandtl number (air) [n/d]
!< Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3
!< roughness model coeff. [n/d]
!< --- transitional mode
!< B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m
!< --- fully rough mode (Re > Re_rough_min)
!< B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
real, parameter :: B_max_land = 2.0
contains
! thermal roughness definition by Kazakov, Lykosov
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_kl(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
! ----------------------------------------------------------------------------
!--- define B = log(z0_m / z0_t)
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
B = min(B, B_max_land)
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Chen, F., Zhang, Y., 2009.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_cz(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45) !Chen and Zhang
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Zilitinkevich, S., 1995.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_zi(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=0.1*kappa*(Re**0.5) !6-Zilitinkevich
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
! It is better to use for dynamic surfaces such as sand
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_ca(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-3.8 !4-Cahill et al.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Owen P. R., Thomson W. R., 1963.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_ot(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=kappa*(Re**0.45) !Owen P. R., Thomson W. R.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Duynkerke P. G., 1992.
!It is better to use for surfaces wiht forest
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_du(z0_t, B, &
z0_m, u_dyn, LAI)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: u_dyn !< dynamic velocity [m/s]
real, intent(in) :: LAI !< leaf-area index
B=(13*u_dyn**0.4)/LAI+0.85 !Duynkerke P. G., 1992.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition z0_t = C*z0_m
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_zm(z0_t, B, &
z0_m, Czm)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Czm !< proportionality coefficient
z0_t =Czm*z0_m
B=log(z0_m / z0_t)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Chen and Zhang and Zilitinkevich
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_mix(z0_t, B, &
z0_m, u_dyn, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: u_dyn !< dynamic velocity [m/s]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
real, parameter :: u_dyn_th=0.17 !< dynamic velocity treshhold [m/s]
if (u_dyn <= u_dyn_th) then
B=0.1*kappa*(Re**0.5) !Zilitinkevich
else
B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45) !Chen and Zhang
end if
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
end module module_z0t
srcF/module_z0t_ocean.f90 deleted 100644 → 0
View file @ b5be41e8
module module_z0t_ocean
!< @brief surface thermal roughness parameterizations for ocean
implicit none
public :: get_thermal_roughness_kl
public :: get_thermal_roughness_ca
public :: get_thermal_roughness_zm
public :: get_thermal_roughness_br
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
real, parameter, private :: Pr_m = 0.71 !< molecular Prandtl number (air) [n/d]
!< Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3
!< roughness model coeff. [n/d]
!< --- transitional mode
!< B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m
!< --- fully rough mode (Re > Re_rough_min)
!< B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
real, parameter :: B_max_ocean = 8.0
contains
! thermal roughness definition by Kazakov, Lykosov
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_kl(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
! ----------------------------------------------------------------------------
!--- define B = log(z0_m / z0_t)
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
B = min(B, B_max_ocean)
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_ca(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-3.8 !4-Cahill et al.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition z0_t = C*z0_m
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_zm(z0_t, B, &
z0_m, Czm)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Czm !< proportionality coefficient
z0_t =Czm*z0_m
B=log(z0_m / z0_t)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Brutsaert W., 2003.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_br(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-2.0 !Brutsaert
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
end module module_z0t_ocean
srcF/module_z0t_snow.f90 deleted 100644 → 0
View file @ b5be41e8
module module_z0t_snow
!< @brief surface thermal roughness parameterizations for snow
implicit none
public :: get_thermal_roughness_kl
public :: get_thermal_roughness_ca
public :: get_thermal_roughness_zm
public :: get_thermal_roughness_br
! --------------------------------------------------------------------------------
real, parameter, private :: kappa = 0.40 !< von Karman constant [n/d]
real, parameter, private :: Pr_m = 0.71 !< molecular Prandtl number (air) [n/d]
!< Re fully roughness minimum value [n/d]
real, parameter :: Re_rough_min = 16.3
!< roughness model coeff. [n/d]
!< --- transitional mode
!< B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
real, parameter :: B1_rough = 5.0 / 6.0
real, parameter :: B2_rough = 0.45
real, parameter :: B3_rough = kappa * Pr_m
!< --- fully rough mode (Re > Re_rough_min)
!< B = B4 * Re^(B2)
real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
real, parameter :: B_max_snow = 8.0
contains
! thermal roughness definition by Kazakov, Lykosov
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_kl(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
! ----------------------------------------------------------------------------
!--- define B = log(z0_m / z0_t)
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
B = min(B, B_max_snow)
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_ca(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-3.8 !4-Cahill et al.
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition z0_t = C*z0_m
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_zm(z0_t, B, &
z0_m, Czm)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Czm !< proportionality coefficient
z0_t =Czm*z0_m
B=log(z0_m / z0_t)
end subroutine
! --------------------------------------------------------------------------------
! thermal roughness definition by Brutsaert W., 2003.
! --------------------------------------------------------------------------------
subroutine get_thermal_roughness_br(z0_t, B, &
z0_m, Re)
! ----------------------------------------------------------------------------
real, intent(out) :: z0_t !< thermal roughness [m]
real, intent(out) :: B !< = log(z0_m / z0_t) [n/d]
real, intent(in) :: z0_m !< aerodynamic roughness [m]
real, intent(in) :: Re !< roughness Reynolds number [n/d]
B=2.46*(Re**0.25)-2.0 !Brutsaert
! --- define roughness [thermal]
z0_t = z0_m / exp(B)
end subroutine
end module module_z0t_snow
srcF/sfx_api_inmcm.f90
View file @ 28e6903d
......@@ -20,12 +20,16 @@ module sfx_api_inmcm
contains
! --------------------------------------------------------------------------------
subroutine inmcm_to_sfx_in_cell(meteo, arg)
subroutine inmcm_to_sfx_in_cell(meteo, arg, IVEG_sfx, depth_inm, lai_inm)
!> @brief converts legacy arg [AR1 INMCM format] array to sfx meteo input
! ----------------------------------------------------------------------------
implicit none
type (meteoDataType), intent(inout) :: meteo
real, dimension(6), intent(in) :: arg
integer,intent(in) :: IVEG_sfx
real,intent(in) :: depth_inm
real,intent(in) :: lai_inm
! ----------------------------------------------------------------------------
......@@ -35,7 +39,10 @@ contains
meteo%dQ = arg(4)
meteo%h = arg(5)
meteo%z0_m = arg(6)
meteo%depth = depth_inm
meteo%lai = lai_inm
meteo%surface_type = IVEG_sfx
!write(*,*) 'surface_type, IVEG_sfx', meteo%surface_type, IVEG_sfx
end subroutine inmcm_to_sfx_in_cell
! --------------------------------------------------------------------------------
......
srcF/sfx_config.f90
View file @ 28e6903d
#include "../includeF/sfx_def.fi"
!> @brief surface flux model config subroutines
module sfx_config
module sfx_config
! modules used
! --------------------------------------------------------------------------------
......@@ -11,7 +11,7 @@ module sfx_config
! --------------------------------------------------------------------------------
implicit none
! --------------------------------------------------------------------------------
public :: set_dataset
public
!> @brief model enum def.
......@@ -20,12 +20,16 @@ module sfx_config
integer, parameter :: model_most = 2 !< MOST model
integer, parameter :: model_sheba = 3 !< SHEBA model
integer, parameter :: model_sheba_noit = 4 !< SHEBA model noniterative
integer, parameter :: model_most_snow = 5 !< MOST_SNOW model
integer, parameter :: model_sheba_coare = 6 !< MOST_SNOW model
character(len = 16), parameter :: model_esm_tag = 'esm'
character(len = 16), parameter :: model_log_tag = 'log'
character(len = 16), parameter :: model_most_tag = 'most'
character(len = 16), parameter :: model_sheba_tag = 'sheba'
character(len = 16), parameter :: model_sheba_noit_tag = 'sheba_noit'
character(len = 16), parameter :: model_most_snow_tag = 'most_snow'
character(len = 16), parameter :: model_sheba_coare_tag = 'sheba_coare'
!> @brief dataset enum def.
integer, parameter :: dataset_mosaic = 1 !< MOSAiC campaign
......@@ -49,8 +53,8 @@ module sfx_config
character(len = 256) :: filename
integer :: nmax
integer :: surface
real :: h, z0_m, z0_h
integer :: surface, surface_type
real :: h, z0_m, z0_h, lai, depth
end type
......@@ -72,7 +76,11 @@ contains
id = model_sheba
else if (trim(tag) == trim(model_sheba_noit_tag)) then
id = model_sheba_noit
end if
else if (trim(tag) == trim(model_most_snow_tag)) then
id = model_most_snow
else if (trim(tag) == trim(model_sheba_coare_tag)) then
id = model_sheba_coare
end if
end function
......@@ -92,6 +100,10 @@ contains
tag = model_sheba_tag
else if (id == model_sheba_noit) then
tag = model_sheba_noit_tag
else if (id == model_most_snow) then
tag = model_most_snow_tag
else if (id == model_sheba_coare) then
tag = model_sheba_coare_tag
end if
end function
......@@ -159,12 +171,16 @@ contains
dataset%filename = get_dataset_filename(id)
dataset%nmax = 0
dataset%surface = surface_land
dataset%surface = surface_snow
dataset%surface_type = surface_snow
dataset%z0_h = -1.0
dataset%lai = 1.0
dataset%depth = 10.0
if (id == dataset_mosaic) then
dataset%h = 3.8
dataset%z0_m = 0.01
dataset%h = 6.0
dataset%z0_m = 0.0078
dataset%surface_type = surface_snow
else if (id == dataset_irgason) then
dataset%h = 9.2
dataset%z0_m = 0.01
......@@ -178,15 +194,15 @@ contains
dataset%z0_m = -1.0
else if (id == dataset_papa) then
dataset%h = 10.0
dataset%surface = surface_ocean
dataset%surface = surface_lake
dataset%z0_m = -1.0
else if (id == dataset_toga) then
! *: check & fix values
dataset%h = 15.0
dataset%surface = surface_ocean
dataset%surface = surface_lake
dataset%z0_m = -1.0
end if
write (*,*) dataset%surface, 'dataset%surface'
end subroutine
function get_dataset_filename(id) result(filename)
......
srcF/sfx_data.f90
View file @ 28e6903d
......@@ -35,6 +35,9 @@ module sfx_data
real(C_FLOAT) :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
real(C_FLOAT) :: dQ !< difference between humidity at 'h' and at surface [g/g]
real(C_FLOAT) :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
real(C_FLOAT) :: depth
real(C_FLOAT) :: lai
integer(C_INT) :: surface_type
end type
!> @brief meteorological input for surface flux calculation
......@@ -46,6 +49,9 @@ module sfx_data
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)
real, allocatable :: depth(:)
real, allocatable :: lai(:)
integer, allocatable :: surface_type(:)
end type
#if defined(INCLUDE_CXX)
......@@ -57,6 +63,10 @@ module sfx_data
type(C_PTR) :: Tsemi !< semi-sum of potential temperature at 'h' and at surface [K]
type(C_PTR) :: dQ !< difference between humidity at 'h' and at surface [g/g]
type(C_PTR) :: z0_m !< surface aerodynamic roughness (should be < 0 for water bodies surface)
type(C_PTR) :: depth
type(C_PTR) :: lai
type(C_PTR) :: surface_type
end type
#endif
! --------------------------------------------------------------------------------
......@@ -157,6 +167,9 @@ contains
allocate(meteo%Tsemi(n))
allocate(meteo%dQ(n))
allocate(meteo%z0_m(n))
allocate(meteo%depth(n))
allocate(meteo%lai(n))
allocate(meteo%surface_type(n))
end subroutine allocate_meteo_vec
......@@ -173,6 +186,9 @@ contains
meteo_C%Tsemi = c_loc(meteo%Tsemi)
meteo_C%dQ = c_loc(meteo%dQ)
meteo_C%z0_m = c_loc(meteo%z0_m)
meteo_C%depth = c_loc(meteo%depth)
meteo_C%lai = c_loc(meteo%lai)
meteo_C%surface_type = c_loc(meteo%surface_type)
end subroutine set_meteo_vec_c
#endif
! --------------------------------------------------------------------------------
......@@ -190,6 +206,9 @@ contains
deallocate(meteo%Tsemi)
deallocate(meteo%dQ)
deallocate(meteo%z0_m)
deallocate(meteo%depth)
deallocate(meteo%lai)
deallocate(meteo%surface_type)
end subroutine deallocate_meteo_vec
! --------------------------------------------------------------------------------
......
srcF/sfx_esm.f90
View file @ 28e6903d
......@@ -27,6 +27,8 @@ module sfx_esm
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
......@@ -139,7 +141,7 @@ contains
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))
z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
......@@ -174,6 +176,9 @@ contains
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)
real :: depth
real :: lai
integer :: surface_type
! ----------------------------------------------------------------------------
! --- local variables
......@@ -202,10 +207,9 @@ contains
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
integer surface_type !< surface type = (ocean || land)
real fval !< just a shortcut for partial calculations
real z0_m1
#ifdef SFX_CHECK_NAN
real NaN
#endif
......@@ -231,32 +235,31 @@ contains
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define thermal roughness & B = log(z0_m / z0_h)
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
write (*,*) surface_type, 'esm'
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
!write (*,*) surface_type, z0_m, z0m_id, z0t_id, z0_t, 'hi3'
!write (*,*) u_dyn0, 'u_dyn0'
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
h0_m = h / z0_m
! --- define relative height [thermal]
h0_t = h / z0_t
......
srcF/sfx_log.f90
View file @ 28e6903d
......@@ -23,6 +23,8 @@ module sfx_log
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
......@@ -55,7 +57,7 @@ contains
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))
z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
......@@ -88,6 +90,8 @@ contains
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)
real :: depth
real :: lai
! ----------------------------------------------------------------------------
! --- local variables
......@@ -138,31 +142,25 @@ contains
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define thermal roughness & B = log(z0_m / z0_h)
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
h0_m = h / z0_m
! --- define relative height [thermal]
h0_t = h / z0_t
......
srcF/sfx_most.f90
View file @ 28e6903d
......@@ -24,6 +24,8 @@ module sfx_most
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
public :: get_psi
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
......@@ -56,7 +58,7 @@ contains
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))
z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
......@@ -89,6 +91,8 @@ contains
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)
real :: depth
real lai
! ----------------------------------------------------------------------------
! --- local variables
......@@ -113,7 +117,7 @@ contains
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
integer surface_type !< surface type = (ocean || land)
real z0_m1
#ifdef SFX_CHECK_NAN
real NaN
......@@ -140,33 +144,31 @@ contains
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define thermal roughness & B = log(z0_m / z0_h)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
! --- define relative height [thermal]
h0_t = h / z0_t
! --- define Ri-bulk
......
srcF/sfx_most_snow.f90 0 → 100644
View file @ 28e6903d
#include "../includeF/sfx_def.fi"
module sfx_most_snow
!< @brief MOST surface flux module
! modules used
! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
use sfx_data
use sfx_surface
use sfx_most_snow_param
use sfx_config
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
public :: get_psi
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: numericsType
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
integer :: maxiters_snow = 10 !< maximum (actual) number of iterations in snow roughness
end type
! --------------------------------------------------------------------------------
contains
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface
! ----------------------------------------------------------------------------
type (sfxDataVecType), intent(inout) :: sfx
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
integer, intent(in) :: n
! ----------------------------------------------------------------------------
! --- local variables
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
integer i
! ----------------------------------------------------------------------------
do i = 1, n
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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
end do
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes(sfx, meteo, numerics)
!< @brief surface flux calculation for single cell
!< @details contains C/C++ interface
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use ieee_arithmetic
#endif
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
! --- meteo derived datatype name shadowing
! ----------------------------------------------------------------------------
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)
real :: depth
real :: lai
! ----------------------------------------------------------------------------
! --- local variables
! ----------------------------------------------------------------------------
real z0_t !< 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 u_dyn0 !< dynamic velocity in neutral conditions [m/s]
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real zeta !< = z/L [n/d]
real Rib !< bulk Richardson number
real Udyn, Tdyn, Qdyn !< dynamic scales
real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d]
real Km !< eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !< invese Prandt number [n/d]
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
real z0_m1
! --- local additional variables
! ----------------------------------------------------------------------------
!real phi_m, phi_m
real hfx, mfx
real S_mean, Lsnow
real z0_s, h_salt
!real S_salt
integer surface_type !< surface type = (ocean || land)
#ifdef SFX_CHECK_NAN
real NaN
#endif
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
! --- checking if arguments are finite
if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
.and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then
!NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN
sfx = sfxDataType(zeta = NaN, Rib = NaN, &
Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
Rib_conv_lim = NaN, &
Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
return
end if
#endif
! --- shadowing names for clarity
U = meteo%U
Tsemi = meteo%Tsemi
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
h0_m = h / z0_m
! --- define relative height [thermal]
h0_t = h / z0_t
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
! --- get the fluxes
! ----------------------------------------------------------------------------
call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
Lsnow, S_mean, h_salt, surface_type, &
U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
! ----------------------------------------------------------------------------
call get_phi(phi_m, phi_h, zeta)
! ----------------------------------------------------------------------------
! --- define transfer coeff. (momentum) & (heat)
Cm = 0.0
if (U > 0.0) then
Cm = Udyn / U
end if
Ct = 0.0
if (abs(dT) > 0.0) then
Ct = Tdyn / dT
end if
! --- define eddy viscosity & inverse Prandtl number
Km = kappa * Cm * U * h / phi_m
Pr_t_inv = phi_m / phi_h
! --- define heat flux and momentum flux
hfx=-rho_air*U*dT*Cm*Ct
mfx=-rho_air*Cm*Cm*U*U
h_salt=h_s
! --- setting output
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
Rib_conv_lim = 0.0, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
! --- setting additional output
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
!< @brief get dynamic scales
! --------------------------------------------------------------------------------
subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
Lsnow, S_mean, h_salt, surface_type, &
U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales
real, intent(out) :: zeta !< = z/L
real, intent(out) :: Lsnow, S_mean
real, intent(out) :: h_salt
integer, intent(in) :: surface_type
real, intent(in) :: U !< abs(wind speed) at z
real, intent(in) :: Tsemi !< semi-sum of temperature at z and at surface
real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface
real, intent(in) :: z !< constant flux layer height
real, intent(in) :: z0_m, z0_t !< roughness parameters
real, intent(in) :: beta !< buoyancy parameter
integer, intent(in) :: maxiters !< maximum number of iterations
! ----------------------------------------------------------------------------
! --- local variables
real, parameter :: gamma = 0.61
real :: psi_m, psi_h
real :: psi0_m, psi0_h
real :: Linv
real :: w_snow, sigma_m, S_salt, d_s, a, b
integer :: i, n
! ----------------------------------------------------------------------------
Udyn = kappa * U / log(z / z0_m)
Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t)
Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t)
zeta = 0.0
! --- no wind
if (Udyn < 1e-5) return
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
! --- near neutral case
if (Linv < 1e-5) return
do i = 1, maxiters
call get_psi(psi_m, psi_h, zeta)
call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)
Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
if (Udyn < 1e-5) exit
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
!S_salt =0.0004
call S_mean_fun(S_salt, Udyn)
if (surface_type==3) then
if (Udyn>u_thsnow) then
call simpson_integration(d_s, 0.00000886, 0.0001, 1000);
call get_S_mean(S_mean, S_salt, h_s, z)
call get_sigma_m(sigma_m, rho_air, rho_s)
call get_w_snow(w_snow, sigma_m, g, d_s, nu_air)
Linv=Linv*((1-S_mean)/(1+sigma_m*S_mean))+(g*w_snow*sigma_m*S_mean/(Udyn**3.0))/(1+sigma_m*S_mean)
zeta = z * Linv
Lsnow=1/Linv
endif
else
Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
endif
end do
end subroutine get_dynamic_scales
! --------------------------------------------------------------------------------
! stability functions
! --------------------------------------------------------------------------------
subroutine get_phi(phi_m, phi_h, zeta)
!< @brief stability functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< stability functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
if (zeta >= 0.0) then
phi_m = 1.0 + beta_m * zeta
phi_h = 1.0 + beta_h * zeta
else
phi_m = (1.0 - alpha_m * zeta)**(-0.25)
phi_h = (1.0 - alpha_h * zeta)**(-0.5)
end if
end subroutine
! --------------------------------------------------------------------------------
subroutine get_S_mean(S_mean, S_salt, h_s, z)
!< @brief function for snow consentration
! ----------------------------------------------------------------------------
real, intent(out) :: S_mean !< snow consentration
real, intent(in) :: S_salt, h_s, z
! ----------------------------------------------------------------------------
S_mean = S_salt * h_s/z
end subroutine
! --------------
subroutine get_sigma_m(sigma_m, rho_air, rho_s)
!< @brief function for
! ----------------------------------------------------------------------------
real, intent(out) :: sigma_m !< s
real, intent(in) :: rho_air, rho_s
! ----------------------------------------------------------------------------
sigma_m = (rho_s - rho_air)/rho_air
end subroutine
subroutine get_w_snow(w_snow, sigma_m, g, d_s, nu_air)
!< @brief function for smow velosity
! ----------------------------------------------------------------------------
real, intent(out) :: w_snow !<
real, intent(in) :: sigma_m, g, d_s, nu_air
! ----------------------------------------------------------------------------
w_snow = sigma_m * g * d_s * d_s / (18.0 * nu_air);
end subroutine
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi(psi_m, psi_h, zeta)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
! ----------------------------------------------------------------------------
if (zeta >= 0.0) then
psi_m = -beta_m * zeta;
psi_h = -beta_h * zeta;
else
x_m = (1.0 - alpha_m * zeta)**(0.25)
x_h = (1.0 - alpha_h * zeta)**(0.25)
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
end if
end subroutine
subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta_m, zeta_h !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
! ----------------------------------------------------------------------------
if (zeta_m >= 0.0) then
psi_m = -beta_m * zeta_m
else
x_m = (1.0 - alpha_m * zeta_m)**(0.25)
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
end if
if (zeta_h >= 0.0) then
psi_h = -beta_h * zeta_h;
else
x_h = (1.0 - alpha_h * zeta_h)**(0.25)
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
end if
end subroutine
! --------------------------------------------------------------------------------
subroutine S_mean_fun(S_salt, Udyn)
real, intent(in) :: Udyn
real, intent(out) :: S_salt
real S_salt1
if (Udyn <= u_thsnow) then
S_salt1 = 0.0
else
S_salt1 = (Udyn*Udyn - u_thsnow*u_thsnow) / (3.25 * Udyn * 9.8 * 0.08436 * Udyn**1.27)
end if
S_salt = S_salt1 / (S_salt1 + rho_s / rho_air)
end subroutine
subroutine simpson_integration(d_s, a, b, n)
real, intent(in) :: a, b
integer, intent(in) :: n
real, intent(out) :: d_s
real h, sum
integer i
h = (b - a) / n
sum=(exp(-0.5 * a * a)/sqrt(2 * 3.14))+(exp(-0.5 * b * b)/sqrt(2 * 3.14))
do i = 1, n - 1, 2
sum = sum + 4 * (exp(-0.5 * (a + i * h) * (a + i * h))/sqrt(2 * 3.14))
end do
do i = 2, n - 2, 2
sum = sum + 2 * (exp(-0.5 * (a + i * h) * (a + i * h))/sqrt(2 * 3.14))
end do
d_s = sum * (h / 3)
end subroutine
end module sfx_most_snow
\ No newline at end of file
srcF/sfx_most_snow_param.f90 0 → 100644
View file @ 28e6903d
module sfx_most_snow_param
!< @brief MOST BD71 surface flux model parameters
!< @details all in SI units
! modules used
! --------------------------------------------------------------------------------
use sfx_phys_const
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
! --------------------------------------------------------------------------------
!< von Karman constant [n/d]
real, parameter :: kappa = 0.40
!< inverse Prandtl (turbulent) number in neutral conditions [n/d]
real, parameter :: Pr_t_0_inv = 1.15
!< stability function coeff. (unstable)
real, parameter :: alpha_m = 16.0
real, parameter :: alpha_h = 16.0
!< stability function coeff. (stable)
real, parameter :: beta_m = 4.7
real, parameter :: beta_h = beta_m
!< snow parameters
real, parameter :: rho_s =900
!real, parameter :: d_s=0.00000886
real, parameter :: h_s=0.05
real, parameter :: u_thsnow=0.25
!real, parameter :: S_salt=0.000004
real, parameter :: rho_air=1.2
end module sfx_most_snow_param
srcF/sfx_run.f90
View file @ 28e6903d
......@@ -43,6 +43,12 @@ contains
use sfx_sheba_noniterative, only: &
get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
numericsType_sheba_noit => numericsType
use sfx_most_snow, only: &
get_surface_fluxes_vec_most_snow => get_surface_fluxes_vec, &
numericsType_most_snow => numericsType
use sfx_sheba_coare, only: &
get_surface_fluxes_vec_sheba_coare => get_surface_fluxes_vec, &
numericsType_sheba_coare => numericsType
! --------------------------------------------------------------------------------
character(len=*), intent(in) :: filename_out
......@@ -63,7 +69,9 @@ contains
type(numericsType_most) :: numerics_most !< surface flux module (MOST) numerics parameters
type(numericsType_sheba) :: numerics_sheba !< surface flux module (SHEBA) numerics parameters
type(numericsType_sheba_noit) :: numerics_sheba_noit !< surface flux module (SHEBA) numerics parameters
type(numericsType_most_snow) :: numerics_most_snow !< surface flux module (MOST_SNOW) numerics parameters
type(numericsType_most_snow) :: numerics_sheba_coare !< surface flux module (MOST_SNOW) numerics parameters
integer :: num !< number of 'cells' in input
! --------------------------------------------------------------------------------
......@@ -85,6 +93,9 @@ contains
write(*, '(a,g0)') ' h = ', dataset%h
write(*, '(a,g0)') ' z0(m) = ', dataset%z0_m
write(*, '(a,g0)') ' z0(h) = ', dataset%z0_h
write(*, '(a,g0)') ' lai = ', dataset%lai
write(*, '(a,g0)') ' depth = ', dataset%depth
write(*, '(a,g0)') ' surface_type = ', dataset%surface_type
! --- define number of elements
......@@ -121,6 +132,9 @@ contains
! --- setting height & roughness
meteo_cell%h = dataset%h
meteo_cell%z0_m = dataset%z0_m
meteo_cell%depth = dataset%depth
meteo_cell%lai = dataset%lai
meteo_cell%surface_type = dataset%surface_type
! --- read input data
open(newunit = io, file = dataset%filename, iostat = status, status = 'old')
......@@ -133,11 +147,16 @@ contains
read(io, *) meteo_cell%U, meteo_cell%dT, meteo_cell%Tsemi, meteo_cell%dQ
meteo%h(i) = meteo_cell%h
meteo%depth(i)=meteo_cell%depth
meteo%lai(i)=meteo_cell%lai
meteo%surface_type(i)=meteo_cell%surface_type
meteo%U(i) = meteo_cell%U
meteo%dT(i) = meteo_cell%dT
meteo%Tsemi(i) = meteo_cell%Tsemi
meteo%dQ(i) = meteo_cell%dQ
meteo%z0_m(i) = meteo_cell%z0_m
enddo
close(io)
......@@ -153,6 +172,10 @@ contains
call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, num)
else if (model == model_sheba_noit) then
call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, num)
else if (model == model_most_snow) then
call get_surface_fluxes_vec_most_snow(sfx, meteo, numerics_most_snow, num)
else if (model == model_sheba_coare) then
call get_surface_fluxes_vec_sheba_coare(sfx, meteo, numerics_sheba_coare, num)
end if
......@@ -239,7 +262,7 @@ contains
write(*, *) ' --help'
write(*, *) ' print usage options'
write(*, *) ' --model [key]'
write(*, *) ' key = esm (default) || log || most || sheba || sheba_noit'
write(*, *) ' key = esm (default) || log || most || sheba || sheba_noit || most_snow || sheba_coare'
write(*, *) ' --dataset [key]'
write(*, *) ' key = mosaic (default) || irgason || sheba'
write(*, *) ' = lake || papa || toga || user [filename] [param]'
......@@ -413,6 +436,7 @@ contains
!< save nmax if previously set
nmax = dataset%nmax
call set_dataset(dataset, id)
dataset%nmax = nmax
call c_config_is_varname("dataset.filename"//C_NULL_CHAR, status)
......
srcF/sfx_sheba.f90
View file @ 28e6903d
......@@ -29,6 +29,8 @@ module sfx_sheba
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
public :: get_psi
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
......@@ -138,7 +140,7 @@ contains
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))
z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
......@@ -171,6 +173,8 @@ contains
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)
real :: lai
real :: depth
! ----------------------------------------------------------------------------
! --- local variables
......@@ -223,32 +227,22 @@ contains
dQ = meteo%dQ
h = meteo%h
z0_m = meteo%z0_m
! --- define surface type
if (z0_m < 0.0) then
surface_type = surface_ocean
else
surface_type = surface_land
end if
if (surface_type == surface_ocean) then
! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
! --- define relative height
h0_m = h / z0_m
endif
if (surface_type == surface_land) then
! --- define relative height
h0_m = h / z0_m
! --- define dynamic velocity in neutral conditions
u_dyn0 = U * kappa / log(h0_m)
end if
! --- define thermal roughness & B = log(z0_m / z0_h)
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
! --- define relative height [thermal]
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
h0_t = h / z0_t
! --- define Ri-bulk
......
srcF/sfx_sheba_coare.f90 0 → 100644
View file @ 28e6903d
#include "../includeF/sfx_def.fi"
module sfx_sheba_coare
! modules used
! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
use sfx_data
use sfx_surface
use sfx_sheba_coare_param
#if defined(INCLUDE_CXX)
use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
use C_FUNC
#endif
! --------------------------------------------------------------------------------
! directives list
! --------------------------------------------------------------------------------
implicit none
private
! --------------------------------------------------------------------------------
! public interface
! --------------------------------------------------------------------------------
public :: get_surface_fluxes
public :: get_surface_fluxes_vec
public :: get_psi_stable
public :: get_psi_a
public :: get_psi_convection
public :: get_psi_BD
integer z0m_id
integer z0t_id
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
type, public :: numericsType
integer :: maxiters_charnock = 10 !< maximum (actual) number of iterations in charnock roughness
integer :: maxiters_convection = 10 !< maximum (actual) number of iterations in charnock roughness
end type
! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
type, BIND(C), public :: sfx_sheba_coare_param_C
real(C_FLOAT) :: kappa
real(C_FLOAT) :: Pr_t_0_inv
real(C_FLOAT) :: Pr_t_inf_inv
real(C_FLOAT) :: alpha_m
real(C_FLOAT) :: alpha_h
real(C_FLOAT) :: gamma
real(C_FLOAT) :: zeta_a
real(C_FLOAT) :: a_m
real(C_FLOAT) :: b_m
real(C_FLOAT) :: a_h
real(C_FLOAT) :: b_h
real(C_FLOAT) :: c_h
real(C_FLOAT) :: beta_m
real(C_FLOAT) :: beta_h
end type
type, BIND(C), public :: sfx_sheba_coare_numericsType_C
integer(C_INT) :: maxiters_convection
integer(C_INT) :: maxiters_charnock
end type
INTERFACE
SUBROUTINE c_sheba_coare_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, &
name="c_sheba_coare_compute_flux")
use sfx_data
use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
Import :: sfx_sheba_coare_param_C, sfx_sheba_coare_numericsType_C
implicit none
integer(C_INT) :: grid_size
type(C_PTR), value :: sfx
type(C_PTR), value :: meteo
type(sfx_sheba_coare_param_C) :: model_param
type(sfx_surface_sheba_coare_param) :: surface_param
type(sfx_sheba_coare_numericsType_C) :: numerics
type(sfx_phys_constants) :: constants
END SUBROUTINE c_sheba_coare_compute_flux
END INTERFACE
#endif
contains
! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
subroutine set_c_struct_sfx_sheba_coare_param_values(sfx_model_param)
type (sfx_sheba_coare_param_C), intent(inout) :: sfx_model_param
sfx_model_param%kappa = kappa
sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
sfx_model_param%Pr_t_inf_inv = Pr_t_inf_inv
sfx_model_param%alpha_m = alpha_m
sfx_model_param%alpha_h = alpha_h
sfx_model_param%gamma = gamma
sfx_model_param%zeta_a = zeta_a
sfx_model_param%a_m = a_m
sfx_model_param%b_m = b_m
sfx_model_param%a_h = a_h
sfx_model_param%b_h = b_h
sfx_model_param%c_h = c_h
sfx_model_param%c3 = beta_m
sfx_model_param%c4 = beta_h
end subroutine set_c_struct_sfx_sheba_coare_param_values
#endif
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
!< @brief surface flux calculation for array data
!< @details contains C/C++ & CUDA interface
! ----------------------------------------------------------------------------
type (sfxDataVecType), intent(inout) :: sfx
type (meteoDataVecType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
integer, intent(in) :: n
! ----------------------------------------------------------------------------
! --- local variables
type (meteoDataType) meteo_cell
type (sfxDataType) sfx_cell
integer i
! ----------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
type (meteoDataVecTypeC), target :: meteo_c !< meteorological data (input)
type (sfxDataVecTypeC), target :: sfx_c !< surface flux data (output)
type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
type (sfx_sheba_coare_param_C) :: model_param
type (sfx_surface_param) :: surface_param
type (sfx_sheba_coare_numericsType_C) :: numerics_c
type (sfx_phys_constants) :: phys_constants
numerics_c%maxiters_convection = numerics%maxiters_convection
numerics_c%maxiters_charnock = numerics%maxiters_charnock
phys_constants%Pr_m = Pr_m;
phys_constants%nu_air = nu_air;
phys_constants%g = g;
call set_c_struct_sfx_sheba_coare_param_values(model_param)
call set_c_struct_sfx_surface_param_values(surface_param)
call set_meteo_vec_c(meteo, meteo_c)
call set_sfx_vec_c(sfx, sfx_c)
meteo_c_ptr = C_LOC(meteo_c)
sfx_c_ptr = C_LOC(sfx_c)
call c_sheba_coare_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
do i = 1, n
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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
call push_sfx_data(sfx, sfx_cell, i)
end do
#endif
end subroutine get_surface_fluxes_vec
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
subroutine get_surface_fluxes(sfx, meteo, numerics)
!< @brief surface flux calculation for single cell
!< @details contains C/C++ interface
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
use ieee_arithmetic
#endif
type (sfxDataType), intent(out) :: sfx
type (meteoDataType), intent(in) :: meteo
type (numericsType), intent(in) :: numerics
! ----------------------------------------------------------------------------
! --- meteo derived datatype name shadowing
! ----------------------------------------------------------------------------
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)
!real :: hpbl
real :: depth
real :: lai
integer :: surface_type
! ----------------------------------------------------------------------------
! --- local variables
! ----------------------------------------------------------------------------
real z0_t !< 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 u_dyn0 !< dynamic velocity in neutral conditions [m/s]
real Re !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
real zeta !< = z/L [n/d]
real Rib !< bulk Richardson number
real zeta_conv_lim !< z/L critical value for matching free convection limit [n/d]
real Rib_conv_lim !< Ri-bulk critical value for matching free convection limit [n/d]
real f_m_conv_lim !< stability function (momentum) value in free convection limit [n/d]
real f_h_conv_lim !< stability function (heat) value in free convection limit [n/d]
real psi_m, psi_h !< universal functions (momentum) & (heat) [n/d]
real psi0_m, psi0_h !< universal functions (momentum) & (heat) [n/d]
real z0_m1
! real psi_m_BD, psi_h_BD !< universal functions (momentum) & (heat) [n/d]
! real psi0_m_BD, psi0_h_BD !< universal functions (momentum) & (heat) [n/d]
! real psi_m_conv, psi_h_conv !< universal functions (momentum) & (heat) [n/d]
! real psi0_m_conv, psi0_h_conv !< universal functions (momentum) & (heat) [n/d]
real Udyn, Tdyn, Qdyn !< dynamic scales
real phi_m, phi_h !< stability functions (momentum) & (heat) [n/d]
real Km !< eddy viscosity coeff. at h [m^2/s]
real Pr_t_inv !< invese Prandt number [n/d]
real Cm, Ct !< transfer coeff. for (momentum) & (heat) [n/d]
!integer surface_type !< surface type = (ocean || land)
real c_wdyn
#ifdef SFX_CHECK_NAN
real NaN
#endif
! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
! --- checking if arguments are finite
if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
.and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then
NaN = ieee_value(0.0, ieee_quiet_nan) ! setting NaN
sfx = sfxDataType(zeta = NaN, Rib = NaN, &
Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
Rib_conv_lim = NaN, &
Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
!Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN, c_wdyn = NaN)
return
end if
#endif
! --- shadowing names for clarity
U = meteo%U
Tsemi = meteo%Tsemi
dT = meteo%dT
dQ = meteo%dQ
h = meteo%h
z0_m1 = meteo%z0_m
depth = meteo%depth
lai = meteo%lai
surface_type=meteo%surface_type
!hpbl = meteo%hpbl
call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
Re = u_dyn0 * z0_m / nu_air
call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
! --- define relative height
h0_m = h / z0_m
! --- define relative height [thermal]
h0_t = h / z0_t
! --- define Ri-bulk
Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
! --- define free convection transition zeta = z/L value
call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
h0_m, h0_t, B)
! --- get the fluxes
! ----------------------------------------------------------------------------
if (Rib > 0.0) then
! --- stable stratification block
! --- restrict bulk Ri value
! *: note that this value is written to output
! Rib = min(Rib, Rib_max)
call get_zeta_stable(zeta, Rib, h, z0_m, z0_t)
call get_psi_stable(psi_m, psi_h, zeta, zeta)
call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
else if (Rib <= -0.001)then
call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, psi_m, psi_h, psi0_m, psi0_h,&
U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), numerics%maxiters_convection)
call get_phi_a(phi_m,phi_h,zeta)
!! call get_phi_a2(phi_m,phi_h,zeta)
!! call get_phi_a3(phi_m,phi_h,zeta)
!! print *, zeta,psi_m,psi_h,phi_m,phi_h
psi_m = (log(h / z0_m) - (psi_m - psi0_m))
psi_h = (log(h / z0_t) - (psi_h - psi0_h))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!non-iterative version below is not debugged yet!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! call get_zeta_conv(zeta,Rib,h,z0_m,z0_t)
!
! call get_psi_a(psi_m, psi_h, zeta,zeta)
! call get_psi_a(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
!
! Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
! Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
!
! call get_phi_a(phi_m,phi_h,zeta)
!print *, zeta,psi_m,psi_h,phi_m,phi_h
!
! psi_m = (log(h / z0_m) - (psi_m - psi0_m))
! psi_h = (log(h / z0_t) - (psi_h - psi0_h))
else
! --- nearly neutral [-0.001, 0] block
call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
zeta = 0.0
phi_m = 1.0
phi_h = 1.0 / Pr_t_0_inv
Udyn = kappa * U / log(h / z0_m)
Tdyn = kappa * dT * Pr_t_0_inv / log(h / z0_t)
end if
! ----------------------------------------------------------------------------
! --- define transfer coeff. (momentum) & (heat)
if(Rib > 0)then
Cm = 0.0
if (U > 0.0) then
Cm = Udyn / U
end if
Ct = 0.0
if (abs(dT) > 0.0) then
Ct = Tdyn / dT
end if
else
Cm = kappa / psi_m
Ct = kappa / psi_h
end if
! --- define eddy viscosity & inverse Prandtl number
Km = kappa * Cm * U * h / phi_m
Pr_t_inv = phi_m / phi_h
! --- setting output
sfx = sfxDataType(zeta = zeta, Rib = Rib, &
Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
Rib_conv_lim = Rib_conv_lim, &
Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
! Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv, c_wdyn = 0)
end subroutine get_surface_fluxes
! --------------------------------------------------------------------------------
! universal functions
! --------------------------------------------------------------------------------
subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h
real, intent(in) :: B !< = log(z0_m / z0_h)
! ----------------------------------------------------------------------------
psi_m = log(h0_m)
psi_h = log(h0_t) / Pr_t_0_inv
!*: this looks redundant z0_t = z0_m in case |B| ~ 0
if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
end subroutine
subroutine get_zeta_stable(zeta, Rib, h, z0_m, z0_t)
real,intent(out) :: zeta
real,intent(in) :: Rib, h, z0_m, z0_t
real :: Ribl, C1, A1, A2, lne, lnet
real :: psi_m, psi_h, psi0_m, psi0_h
Ribl = (Rib * Pr_t_0_inv) * (1 - z0_t / h) / ((1 - z0_m / h)**2)
call get_psi_stable(psi_m, psi_h, zeta_a, zeta_a)
call get_psi_stable(psi0_m, psi0_h, zeta_a * z0_m / h, zeta_a * z0_t / h)
lne = log(h/z0_m)
lnet = log(h/z0_t)
C1 = (lne**2)/lnet
A1 = ((lne - psi_m + psi0_m)**(2*(gamma-1))) &
& / ((zeta_a**(gamma-1))*((lnet-(psi_h-psi0_h)*Pr_t_0_inv)**(gamma-1)))
A2 = ((lne - psi_m + psi0_m)**2) / (lnet-(psi_h-psi0_h)*Pr_t_0_inv) - C1
zeta = C1 * Ribl + A1 * A2 * (Ribl**gamma)
end subroutine get_zeta_stable
subroutine get_psi_stable(psi_m, psi_h, zeta_m, zeta_h)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta_m, zeta_h !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
real :: q_m, q_h
! ----------------------------------------------------------------------------
q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
x_m = (1.0 + zeta_m)**(1.0 / 3.0)
psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
2.0 * sqrt(3.0) * (&
atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
q_h = sqrt(c_h * c_h - 4.0)
x_h = zeta_h
psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
log((c_h - q_h) / (c_h + q_h)))
end subroutine get_psi_stable
subroutine get_psi_convection(psi_m, psi_h, zeta_m, zeta_h)
!< @brief Carl et al. 1973 with Grachev et al. 2000 corrections of beta_m, beta_h
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions [n/d]
real, intent(in) :: zeta_m, zeta_h !< = z/L [n/d]
real y
! ----------------------------------------------------------------------------
! beta_m = 10, beta_h = 34
y = (1.0 - beta_m * zeta_m)**(1/3.)
psi_m = 3.0 * 0.5 *log((y*y + y + 1.0)/3.) - sqrt(3.0) *atan((2.0*y + 1)/sqrt(3.0)) + pi/sqrt(3.0)
y = (1.0 - beta_h * zeta_h)**(1/3.)
psi_h = 3.0 * 0.5 *log((y*y + y + 1.0)/3.) - sqrt(3.0) *atan((2.0*y + 1)/sqrt(3.0)) + pi/sqrt(3.0)
end subroutine
subroutine get_psi_BD(psi_m, psi_h, zeta_m, zeta_h)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta_m,zeta_h !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h
! ----------------------------------------------------------------------------
x_m = (1.0 - alpha_m * zeta_m)**(0.25)
x_h = (1.0 - alpha_h * zeta_h)**(0.25)
psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
end subroutine
subroutine get_phi_a(phi_m, phi_h, zeta)
!< @brief universal functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< universal functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: x_m, x_h, y
real :: psi_m_bd,psi_h_bd,psi_m_conv,psi_h_conv
real :: dpsi_m_bd,dpsi_h_bd,dpsi_m_conv,dpsi_h_conv
! ----------------------------------------------------------------------------
call get_psi_BD(psi_m_bd,psi_h_bd,zeta,zeta)
call get_psi_convection(psi_m_conv,psi_h_conv,zeta,zeta)
x_m = (1.0 - alpha_m * zeta)**(0.25)
x_h = (1.0 - alpha_h * zeta)**(0.25)
dpsi_m_bd = -(alpha_m/(2.0*(x_m**3))) * (1/(1+x_m) + (x_m-1)/(1+x_m**2))
dpsi_h_bd = -alpha_h / ((x_h**2)*(1+x_h**2))
y = (1 - beta_m * zeta)**(1/3.)
dpsi_m_conv = -beta_m/(y*(y**2 + y + 1))
y = (1 - beta_h * zeta)**(1/3.)
dpsi_h_conv = -beta_h/(y*(y**2 + y + 1))
phi_m = 1.0 - zeta * (dpsi_m_bd/(1.0+zeta**2) - psi_m_bd*2.0*zeta/((1.0+zeta**2)**2) + &
dpsi_m_conv/(1.0+1.0/(zeta**2)) + 2.0*psi_m_conv/((zeta**3)*((1.0+1.0/(zeta**2))**2)))
phi_h = 1.0 - zeta * (dpsi_h_bd/(1.0+zeta**2) - psi_h_bd*2.0*zeta/((1.0+zeta**2)**2) + &
dpsi_h_conv/(1.0+1.0/(zeta**2)) + 2.0*psi_h_conv/((zeta**3)*((1.0+1.0/(zeta**2))**2)))
end subroutine
subroutine get_phi_a2(phi_m,phi_h,zeta)
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< universal functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: phi_m_bd,phi_h_bd,phi_m_conv,phi_h_conv
call get_phi_convection(phi_m_conv, phi_h_conv, zeta)
call get_phi_BD(phi_m_BD, phi_h_BD, zeta)
phi_m = (phi_m_BD + (zeta**2) * phi_m_conv) / (1 + zeta**2)
phi_h = (phi_h_BD + (zeta**2) * phi_h_conv) / (1 + zeta**2)
end subroutine
subroutine get_phi_a3(phi_m,phi_h,zeta)
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< universal functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: phi_m_bd,phi_h_bd,phi_m_conv,phi_h_conv
real :: psi_m_a,psi_h_a,psi_m_conv,psi_h_conv,psi_m_BD,psi_h_BD
call get_phi_convection(phi_m_conv, phi_h_conv, zeta)
call get_phi_BD(phi_m_BD, phi_h_BD, zeta)
call get_psi_convection(psi_m_conv, psi_h_conv, zeta,zeta)
call get_psi_BD(psi_m_BD, psi_h_BD, zeta,zeta)
! call get_psi_a(psi_m_a, psi_h_a, zeta,zeta)
phi_m = (1-phi_m_BD)/(zeta*(1+zeta**2)) + 2*zeta*(psi_m_conv-psi_m_BD)/((1+zeta**2)**2) &
+ zeta*(1-phi_m_conv)/((1+zeta**2)**2)
phi_h = (1-phi_h_BD)/(zeta*(1+zeta**2)) + 2*zeta*(psi_h_conv-psi_h_BD)/((1+zeta**2)**2) &
+ zeta*(1-phi_h_conv)/((1+zeta**2)**2)
end subroutine
subroutine get_phi_BD(phi_m, phi_h, zeta)
!< @brief stability functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< stability functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
phi_m = (1.0 - alpha_m * zeta)**(-0.25)
phi_h = (1.0 - alpha_h * zeta)**(-0.5)
end subroutine
subroutine get_phi_convection(phi_m, phi_h, zeta)
!< @brief stability functions (momentum) & (heat): neutral case
! ----------------------------------------------------------------------------
real, intent(out) :: phi_m, phi_h !< stability functions
real, intent(in) :: zeta !< = z/L
! ----------------------------------------------------------------------------
phi_m = (1.0 - beta_m * zeta)**(-1.0/3.0)
phi_h = (1.0 - beta_h * zeta)**(-1.0/3.0)
end subroutine
!< @brief get dynamic scales
! --------------------------------------------------------------------------------
subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, psi_m, psi_h, &
psi0_m,psi0_h, U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
! ----------------------------------------------------------------------------
real, intent(out) :: Udyn, Tdyn, Qdyn !< dynamic scales
real, intent(out) :: zeta !< = z/L
real, intent(out) :: psi_m,psi_h,psi0_m,psi0_h
real, intent(in) :: U !< abs(wind speed) at z
real, intent(in) :: Tsemi !< semi-sum of temperature at z and at surface
real, intent(in) :: dT, dQ !< temperature & humidity difference between z and at surface
real, intent(in) :: z !< constant flux layer height
real, intent(in) :: z0_m, z0_t !< roughness parameters
real, intent(in) :: beta !< buoyancy parameter
integer, intent(in) :: maxiters !< maximum number of iterations
! ----------------------------------------------------------------------------
! --- local variables
real :: Linv
integer :: i
! ----------------------------------------------------------------------------
Udyn = kappa * U / log(z / z0_m)
Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t)
Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t)
zeta = 0.0
! --- no wind
if (Udyn < 1e-5) return
Linv = kappa * beta * (Tdyn + 0.61 * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
psi_m = log(z/z0_m)
psi_h = log(z/z0_t) / Pr_t_0_inv
psi0_m = log(z/z0_m)
psi0_h = log(z/z0_t) / Pr_t_0_inv
! --- near neutral case
! if (Linv < 1e-5) return
do i = 1, maxiters
call get_psi_a(psi_m, psi_h, zeta, zeta)
call get_psi_a(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)
Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
if (Udyn < 1e-5) exit
Linv = kappa * beta * (Tdyn + 0.61 * Qdyn * Tsemi) / (Udyn * Udyn)
zeta = z * Linv
end do
end subroutine get_dynamic_scales
subroutine get_psi_a(psi_m,psi_h,zeta_m, zeta_h)
! ----------------------------------------------------------------------------
real, intent(out) :: psi_m, psi_h !< universal functions
real, intent(in) :: zeta_m,zeta_h !< = z/L
! ----------------------------------------------------------------------------
! --- local variables
real :: psi_m_bd,psi_h_bd,psi_m_conv,psi_h_conv
call get_psi_convection(psi_m_conv, psi_h_conv, zeta_m, zeta_h)
call get_psi_BD(psi_m_BD, psi_h_BD, zeta_m, zeta_h)
psi_m = (psi_m_BD + (zeta_m**2) * psi_m_conv) / (1 + zeta_m**2)
psi_h = (psi_h_BD + (zeta_h**2) * psi_h_conv) / (1 + zeta_h**2)
end subroutine
subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, &
h0_m, h0_t, B)
! ----------------------------------------------------------------------------
real, intent(out) :: zeta_lim !< limiting value of z/L
real, intent(out) :: Rib_lim !< limiting value of Ri-bulk
real, intent(out) :: f_m_lim, f_h_lim !< limiting values of universal functions shortcuts
real, intent(in) :: h0_m, h0_t !< = z/z0_m, z/z0_h [n/d]
real, intent(in) :: B !< = log(z0_m / z0_h) [n/d]
! ----------------------------------------------------------------------------
! --- local variables
real :: psi_m, psi_h
real :: f_m, f_h
real :: c
! ----------------------------------------------------------------------------
! --- define limiting value of zeta = z / L
c = (Pr_t_inf_inv / Pr_t_0_inv)**4
zeta_lim = (2.0 * alpha_h - c * alpha_m - &
sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2)
f_m_lim = f_m_conv(zeta_lim)
f_h_lim = f_h_conv(zeta_lim)
! --- universal functions
f_m = zeta_lim / h0_m
f_h = zeta_lim / h0_t
if (abs(B) < 1.0e-10) f_h = f_m
f_m = (1.0 - alpha_m * f_m)**0.25
f_h = sqrt(1.0 - alpha_h_fix * f_h)
psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)))
psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv
! --- bulk Richardson number
Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
end subroutine
! convection universal functions shortcuts
! --------------------------------------------------------------------------------
function f_m_conv(zeta)
! ----------------------------------------------------------------------------
real :: f_m_conv
real, intent(in) :: 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, intent(in) :: zeta
! ----------------------------------------------------------------------------
f_h_conv = (1.0 - alpha_h * zeta)**0.5
end function f_h_conv
subroutine get_zeta_conv(zeta,Rib,z,z0m,z0t)
!< @brief Srivastava and Sharan 2017, Abdella and Assefa 2005
! ----------------------------------------------------------------------------
real, intent(out) :: zeta !< = z/L [n/d]
real, intent(in) :: Rib !
real, intent(in) :: z,z0m,z0t !
real A,a0,a1,a2,r,q,s1,s2,theta,delta
real ksi_m,ksi_h,ksi_m_0,ksi_m_inf,ksi_h_0,ksi_h_inf
real f_m_inf,f_h_inf
real psi_m_zeta,psi_m_zeta0,psi_h_zeta,psi_h_zeta0
! ----------------------------------------------------------------------------
A = ( 1 / Pr_t_0_inv ) * ( (1 - z0m/z)**2) * log(z/z0t) / ( (1 - z0t/z) * ((log(z/z0m))**2) )
call get_psi_convection(psi_m_zeta,psi_h_zeta,Rib/A, Rib/A)
call get_psi_convection(psi_m_zeta0,psi_h_zeta0, (z0m/z) * (Rib/A),(z0t/z) * (Rib/A))
f_m_inf = 1 - (psi_m_zeta - psi_m_zeta0) / log(z/z0m)
f_h_inf = 1 - (psi_h_zeta - psi_h_zeta0) / log(z/z0t)
ksi_m_0 = ((z0m / z) - 1.0) / log(z0m/z)
ksi_h_0 = ((z0t / z) - 1.0) / log(z0t/z)
ksi_m_inf = (A / (beta_m * Rib)) * (1.0 - 1.0 / (f_m_inf**4))
ksi_h_inf = (A / (beta_h * Rib)) * (1.0 - ((1.0 / Pr_t_0_inv)**2) / (f_h_inf**2))
ksi_m = cc1 * ksi_m_inf + cc2 * ksi_m_0
ksi_h = cc3 * ksi_h_inf + cc4 * ksi_h_0
a0 = (1 / (beta_m * ksi_m)) * ((Rib / A)**2)
a1 = -1 * (beta_h * ksi_h / (beta_m * ksi_m)) * ((Rib / A)**2)
a2 = -1 / (beta_m * ksi_m)
r = (9.0 * (a1 * a2 - 3 * a0) - 2.0 * (a2**3)) / 54.0
q = (3.0 * a1 - a2*a2) / 9.0
delta = q**3 + r**2
s1 = (r + sqrt(delta))**(1./3.)
s2 = (r - sqrt(delta))**(1./3.)
theta = 1.0 / cos(r / sqrt(-1 * (q**3)))
if(delta <= 0.0)then
zeta = 2.0 * sqrt(-1.0 * q) * cos((theta + 2.0 * pi)/3.0) + 1 /(3.0 * beta_m * ksi_m)
else
s1 = cmplx(r + delta**0.5)**(1./3.)
s2 = cmplx(r - delta**0.5)**(1./3.)
zeta = -1*(real(real(s1)) + real(real(s2)) + 1.0 /(3.0 * beta_m * ksi_m))
endif
end subroutine
end module sfx_sheba_coare
\ No newline at end of file