From 93ffc2f0a80519746971a18f22292a02ddc3b21a Mon Sep 17 00:00:00 2001
From: Evgeny Mortikov <evgeny.mortikov@gmail.com>
Date: Sun, 17 Dec 2023 21:16:45 +0300
Subject: [PATCH] major code update
---
COMP_DRAG.txt | 2 +-
compile2.sh | 11 +-
drag3.f90 | 505 -----------------------------
makefile | 4 +-
param.f90 | 69 ----
sfx_esm.f90 | 593 ++++++++++++++++++++++++++++++++++
sfx_esm_param.f90 | 71 ++++
main_drag.f90 => sfx_main.f90 | 100 +++---
sfx_phys_const.f90 | 12 +
9 files changed, 734 insertions(+), 633 deletions(-)
delete mode 100644 drag3.f90
delete mode 100644 param.f90
create mode 100644 sfx_esm.f90
create mode 100644 sfx_esm_param.f90
rename main_drag.f90 => sfx_main.f90 (62%)
create mode 100644 sfx_phys_const.f90
diff --git a/COMP_DRAG.txt b/COMP_DRAG.txt
index 8e9d5fc..0dd43a3 100644
--- a/COMP_DRAG.txt
+++ b/COMP_DRAG.txt
@@ -1,4 +1,4 @@
-ifort -c inputdata.f90
+fort -c inputdata.f90
ifort -c param.f90
ifort -c prmt.f90
ifort -c drag3.F
diff --git a/compile2.sh b/compile2.sh
index e874355..0e6d889 100755
--- a/compile2.sh
+++ b/compile2.sh
@@ -1,9 +1,10 @@
#!/bin/bash
rm drag_ddt.exe *.o
-gfortran -c -Wuninitialized inputdata.f90
-gfortran -c -Wuninitialized param.f90
-gfortran -c -Wuninitialized drag3.f90
-gfortran -c -Wuninitialized main_drag.f90
-gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o
+gfortran -c -cpp -Wuninitialized inputdata.f90
+gfortran -c -cpp -Wuninitialized sfx_phys_const.f90
+gfortran -c -cpp -Wuninitialized sfx_esm_param.f90
+gfortran -c -cpp -Wuninitialized sfx_esm.f90
+gfortran -c -Wuninitialized sfx_main.f90
+gfortran -o drag_ddt.exe sfx_main.o sfx_esm.o inputdata.o sfx_esm_param.o sfx_phys_const.o
diff --git a/drag3.f90 b/drag3.f90
deleted file mode 100644
index c93706f..0000000
--- a/drag3.f90
+++ /dev/null
@@ -1,505 +0,0 @@
- MODULE drag3
- USE param
-
- implicit none
-
- type, public :: meteoDataType
- real :: h ! constant flux layer height [m]
- real :: U ! abs(wind speed) at 'h' [m/s]
- real :: dT ! difference between potential temperature at 'h' and at surface [K]
- real :: Tsemi ! semi-sum of potential temperature at 'h' and at surface [K]
- real :: dQ ! difference between humidity at 'h' and at surface [g/g]
- real :: z0_m ! surface aerodynamic roughness (should be < 0 for water bodies surface)
- end type
-
- type, public :: meteoDataVecType
- real, allocatable :: h(:) ! constant flux layer height [m]
- real, allocatable :: U(:) ! abs(wind speed) at 'h' [m/s]
- real, allocatable :: dT(:) ! difference between potential temperature at 'h' and at surface [K]
- real, allocatable :: Tsemi(:) ! semi-sum of potential temperature at 'h' and at surface [K]
- real, allocatable :: dQ(:) ! difference between humidity at 'h' and at surface [g/g]
- real, allocatable :: z0_m(:) ! surface aerodynamic roughness (should be < 0 for water bodies surface)
- end type
-
- type, public :: sfxDataType
- real :: zeta ! = z/L [n/d]
- real :: Rib ! bulk Richardson number [n/d]
- real :: Re ! Reynolds number [n/d]
- real :: lnzuzt
- real :: z0_m ! aerodynamic roughness [m]
- real :: z0_t ! thermal roughness [m]
- real :: Rib_conv_lim ! Ri-bulk convection critical value [n/d]
- real :: cm
- real :: ch
- real :: ct
- real :: ckt
- end type
-
- type, public :: numericsType
- integer :: maxiters_convection = 10 ! maximum (actual) number of iterations in convection
- integer :: maxiters_charnock = 10 ! maximum (actual) number of iterations in charnock roughness
- end type
-
- type, public:: data_lutyp
- integer, public :: lu_indx=1
- end type
-
- contains
-
- subroutine surf_flux_vec(meteo, &
- masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
- masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
- numerics, lu1,numst)
-
- type (meteoDataVecType), intent(in) :: meteo
- type (numericsType), intent(in) :: numerics
-
- real, dimension (numst) :: masout_zl
- real, dimension (numst) :: masout_ri
- real, dimension (numst) :: masout_re
- real, dimension (numst) :: masout_lnzuzt
- real, dimension (numst) :: masout_zu
- real, dimension (numst) :: masout_ztout
- real, dimension (numst) :: masout_rith
- real, dimension (numst) :: masout_cm
- real, dimension (numst) :: masout_ch
- real, dimension (numst) :: masout_ct
- real, dimension (numst) :: masout_ckt
-
- type (meteoDataType) meteo_cell
- type (sfxDataType) sfx_cell
-
- type (data_lutyp) lu1
- integer i,numst
-
- do i = 1, numst
- meteo_cell = meteoDataType(&
- h = meteo%h(i), &
- U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
- z0_m = meteo%z0_m(i))
-
- call surf_flux(meteo_cell, sfx_cell, numerics, lu1)
-
- masout_zl(i)=sfx_cell%zeta
- masout_ri(i)=sfx_cell%Rib
- masout_re(i)=sfx_cell%Re
- masout_lnzuzt(i)=sfx_cell%lnzuzt
- masout_zu(i)=sfx_cell%z0_m
- masout_ztout(i)=sfx_cell%z0_t
- masout_rith(i)=sfx_cell%Rib_conv_lim
- masout_cm(i)=sfx_cell%cm
- masout_ch(i)=sfx_cell%ch
- masout_ct(i)=sfx_cell%ct
- masout_ckt(i)=sfx_cell%ckt
- enddo
- END SUBROUTINE surf_flux_vec
-
- function f_m_conv(zeta)
- real f_m_conv
- real zeta
-
- f_m_conv = (1.0 - alpha_m * zeta)**0.25
-
- end function f_m_conv
-
- function f_h_conv(zeta)
- real f_h_conv
- real zeta
-
- f_h_conv = (1.0 - alpha_h * zeta)**0.5
-
- end function f_h_conv
-
-
- ! stability functions (neutral)
- ! --------------------------------------------------------------------------------
- subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h, zeta
- real, intent(in) :: h0_m, h0_t
- real, intent(in) :: B
- ! ----------------------------------------------------------------------------
-
- zeta = 0.0
- psi_m = log(h0_m)
- psi_h = log(h0_t) / Pr_t_0_inv
- if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! stability functions (stable)
- ! --------------------------------------------------------------------------------
- subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h, zeta
- real, intent(in) :: Rib
- real, intent(in) :: h0_m, h0_t
- real, intent(in) :: B
-
- real :: Rib_coeff
- real :: psi0_m, psi0_h
- real :: phi
- real :: c
- ! ----------------------------------------------------------------------------
-
- psi0_m = alog(h0_m)
- psi0_h = B / psi0_m
-
- Rib_coeff = beta_m * Rib
- c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
- zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
- (2.0 * beta_m * (1.0 - Rib_coeff))
-
- phi = beta_m * zeta
-
- psi_m = psi0_m + phi
- psi_h = (psi0_m + B) / Pr_t_0_inv + phi
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! stability functions (convection-semistrong)
- ! --------------------------------------------------------------------------------
- subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h, zeta
- real, intent(in) :: Rib
- real, intent(in) :: h0_m, h0_t
- real, intent(in) :: B
- integer, intent(in) :: maxiters
-
- real :: zeta0_m, zeta0_h
- real :: x0, x1, y0, y1
-
- integer :: i
- ! ----------------------------------------------------------------------------
-
- psi_m = log(h0_m)
- psi_h = log(h0_t)
- if (abs(B) < 1.0e-10) psi_h = psi_m
-
- zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
-
- do i = 1, maxiters + 1
- zeta0_m = zeta / h0_m
- zeta0_h = zeta / h0_t
- if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
-
- y1 = (1.0 - alpha_m * zeta)**0.25e0
- x1 = sqrt(1.0 - alpha_h_fix * zeta)
-
- y0 = (1.0 - alpha_m * zeta0_m)**0.25e0
- x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
-
- y0 = max(y0, 1.000001e0)
- x0 = max(x0, 1.000001e0)
-
- psi_m = log((y1 - 1.0e0)*(y0 + 1.0e0)/((y1 + 1.0e0)*(y0 - 1.0e0))) + 2.0e0*(atan(y1) - atan(y0))
- psi_h = log((x1 - 1.0e0)*(x0 + 1.0e0)/((x1 + 1.0e0)*(x0 - 1.0e0))) / Pr_t_0_inv
-
- if (i == maxiters + 1) exit
-
- zeta = Rib * psi_m**2 / psi_h
- end do
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! stability functions (convection)
- ! --------------------------------------------------------------------------------
- subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
- zeta_conv_lim, x10, y10, &
- h0_m, h0_t, B, maxiters)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: psi_m, psi_h, zeta
- real, intent(in) :: Rib
- real, intent(in) :: zeta_conv_lim
- real, intent(in) :: x10, y10
- real, intent(in) :: h0_m, h0_t
- real, intent(in) :: B
- integer, intent(in) :: maxiters
-
- real :: zeta0_m, zeta0_h
- real :: x0, y0
- real :: p1, p0
- real :: a1, b1, c, f
-
- integer :: i
- ! ----------------------------------------------------------------------------
-
- p1 = 2.0 * atan(y10) + log((y10 - 1.0) / (y10 + 1.0))
- p0 = log((x10 - 1.0) / (x10 + 1.0))
-
- zeta = zeta_conv_lim
- do i = 1, maxiters + 1
- zeta0_m = zeta / h0_m
- zeta0_h = zeta / h0_t
- if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
-
- a1 = (zeta_conv_lim / zeta)**(1.0 / 3.0)
- x0 = sqrt(1.0 - alpha_h_fix * zeta0_h)
- y0 = (1.0 - alpha_m * zeta0_m)**0.25
- c = log((x0 + 1.0)/(x0 - 1.0))
- b1 = -2.0*atan(y0) + log((y0 + 1.0)/(y0 - 1.0))
- f = 3.0 * (1.0 - a1)
-
- psi_m = f / y10 + p1 + b1
- psi_h = (f / x10 + p0 + c) / Pr_t_0_inv
-
- if (i == maxiters + 1) exit
- zeta = Rib * psi_m**2 / psi_h
- end do
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! define convection limit
- ! --------------------------------------------------------------------------------
- subroutine get_convection_lim(zeta_lim, Rib_lim, x10, y10, &
- h0_m, h0_t, B)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: zeta_lim, Rib_lim
- real, intent(out) :: x10, y10
-
- real, intent(in) :: h0_m, h0_t
- real, intent(in) :: B
-
- real :: an4, an5
- real :: psi_m, psi_h
- ! ----------------------------------------------------------------------------
-
- an5=(Pr_t_inf_inv / Pr_t_0_inv)**4
- zeta_lim = (2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*(alpha_h-alpha_m))) &
- / (2.0E0*alpha_h**2)
-
- y10 = f_m_conv(zeta_lim)
- x10 = f_h_conv(zeta_lim)
-
- ! ......definition of r-prim......
- an4 = zeta_lim / h0_m
- an5 = zeta_lim / h0_t
- if (abs(B) < 1.0e-10) an5=an4
-
- an5 = sqrt(1.0 - alpha_h_fix * an5)
- an4 = (1.0 - alpha_m * an4)**0.25
-
- psi_m = 2.0 * (atan(y10) - atan(an4)) + alog((y10 - 1.0) * (an4 + 1.0)/((y10 + 1.0) * (an4 - 1.0)))
- psi_h = alog((x10 - 1.0) * (an5 + 1.0)/((x10 + 1.0) * (an5 - 1.0))) / Pr_t_0_inv
-
- Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
- ! charnock roughness
- ! --------------------------------------------------------------------------------
- subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
- ! ----------------------------------------------------------------------------
- real, intent(out) :: z0_m, u_dyn0
-
- real, intent(in) :: U, h
- integer, intent(in) :: maxiters
-
- real :: U10m
- real :: a1, c1, y1, c1min, f
-
- integer :: i, j
- ! ----------------------------------------------------------------------------
-
- U10m = U
- a1 = 0.0e0
- y1 = 25.0e0
- c1min = log(h_charnock) / kappa
- do i = 1, maxiters
- f = c1_charnock - 2.0 * log(U10m)
- do j = 1, maxiters
- c1 = (f + 2.0 * log(y1)) / kappa
- ! looks like the check should use U10 instead of U
- ! but note that a1 should be set = 0 in cycle beforehand
- if(U <= 8.0e0) a1 = log(1.0 + c2_charnock * ((y1 / U10m)**3)) / kappa
- c1 = max(c1 - a1, c1min)
- y1 = c1
- end do
- z0_m = h_charnock * exp(-c1 * kappa)
- z0_m = max(z0_m,0.000015e0)
- U10m = U*alog(h_charnock/z0_m)/(alog(h/z0_m))
- end do
-
- ! --- define dynamic velocity in neutral conditions
- u_dyn0 = U10m / c1
-
- end subroutine
- ! --------------------------------------------------------------------------------
-
-
- SUBROUTINE surf_flux(meteo, sfx, numerics, lu)
- type (sfxDataType), intent(out) :: sfx
-
- type (meteoDataType), intent(in) :: meteo
- type (numericsType), intent(in) :: numerics
-
- type (data_lutyp) lu
-
- real h ! surface flux layer height [m]
-
- real z0_m, z0_t ! aerodynamic & thermal roughness [m]
- real B ! = ln(z0_m / z0_t) [n/d]
- real h0_m, h0_t ! = h / z0_m, h / z0_h [n/d]
-
- real Re ! roughness Reynolds number = u_dyn0 * z0 / nu [n/d]
- real u_dyn0 ! dynamic velocity in neutral conditions [m/s]
-
- real zeta ! = z/L [n/d]
- real zeta_conv_lim ! z/L critical value for matching free convection limit [n/d]
-
- real Rib ! bulk Richardson number
- real Rib_conv_lim ! Ri-bulk critical value for matching free convection limit [n/d]
-
- real psi_m, psi_h ! universal functions [momentum] & [thermal]
-
- real U ! wind velocity at h [m/s]
- real Tsemi !
-
-
- integer lu_indx
-
- ! output ... !
- real an4, o, am
- real c4, c0, an0
- ! ...
-
- ! local unnamed !
- real x10, y10
- real al
-
- real q4, t4
- ! .....
-
- ! just local variables
- real f1, a1
- ! ....
-
- integer is_ocean
- integer i, j
-
-
- Tsemi = meteo%Tsemi
-
- U = meteo%U
- t4 = meteo%dT
- q4 = meteo%dQ
- h = meteo%h
- z0_m = meteo%z0_m
-
- if(z0_m < 0.0) then
- !> @brief Test - definition z0 of sea surface
- !> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
- is_ocean = 1
-
- call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
-
- ! --- define relative height
- h0_m = h / z0_m
- else
- ! ......parameters from viscosity sublayer......
- is_ocean = 0
-
- ! --- define relative height
- h0_m = h / z0_m
- ! --- define dynamic velocity in neutral conditions
- u_dyn0 = U * kappa / log(h0_m)
- end if
-
- Re = u_dyn0 * z0_m / nu_air
- if(Re <= Re_rough_min) then
- B = B1_rough * alog(B3_rough * Re) + B2_rough
- else
- ! B4 takes into account Re value at z' ~ O(10) z0
- B = B4_rough * (Re**B2_rough)
- end if
-
- ! --- apply max restriction based on surface type
- if (is_ocean == 1) then
- B = min(B, B_max_ocean)
- else
- B = min(B, B_max_land)
- end if
-
- ! --- define roughness [thermal]
- z0_t = z0_m / exp(B)
- ! --- define relative height [thermal]
- h0_t = h / z0_t
-
- ! ......humidity stratification and ri-number......
- al = g / Tsemi
- Rib = al * h * (t4 + 0.61e0 * Tsemi * q4) / U**2
-
- ! --- define free convection transition zeta = z/L value
- call get_convection_lim(zeta_conv_lim, Rib_conv_lim, x10, y10, &
- h0_m, h0_t, B)
-
- if (Rib > 0.0e0) then
- ! ......stable stratification......
- !write (*,*) 'stable'
-
- Rib = min(Rib, Rib_max)
- call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
-
- f1 = beta_m * zeta
- am = 1.0 + f1
- o = 1.0/Pr_t_0_inv + f1
-
- else if (Rib < Rib_conv_lim) then
- ! ......strong instability.....
- !write (*,*) 'instability'
-
- call get_psi_convection(psi_m, psi_h, zeta, Rib, &
- zeta_conv_lim, x10, y10, h0_m, h0_t, B, numerics%maxiters_convection)
-
- a1 = (zeta_conv_lim / zeta)**(1.0/3.0)
- am = a1 / y10
- o = a1 / (Pr_t_0_inv * x10)
-
- else if (Rib > -0.001e0) then
-
- ! ......nearly neutral......
- write (*,*) 'neutral'
-
- call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
-
- am = 1.0e0
- o = 1.0e0 / Pr_t_0_inv
- else
- ! ......week and semistrong instability......
- !write (*,*) 'semistrong'
-
- call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
-
- am = (1.0 - alpha_m * zeta)**(-0.25)
- o = 1.0/(Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
- end if
-
- ! ......computation of cu,co,k(h),alft
- c4=kappa/psi_m
- c0=kappa/psi_h
- an4=kappa*c4*U*h/am
- an0=am/o
-
- ! ......exit......
- sfx%zeta = zeta
- sfx%Rib = Rib
- sfx%Re = Re
- sfx%lnzuzt=B
- sfx%z0_m = z0_m
- sfx%z0_t = z0_t
- sfx%Rib_conv_lim = Rib_conv_lim
- sfx%cm=c4
- sfx%ch=c0
- sfx%ct=an4
- sfx%ckt=an0
- return
-
- END SUBROUTINE surf_flux
-
- END MODULE drag3
\ No newline at end of file
diff --git a/makefile b/makefile
index ef11984..bd4ab1b 100644
--- a/makefile
+++ b/makefile
@@ -1,4 +1,4 @@
-RUN = drag.exe
+RUN = sfx
COMPILER ?= gnu
FC_KEYS ?=
@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
-OBJ_F90 = inputdata.o param.o prmt.o
+OBJ_F90 = sfx_phys_const.o sfx_esm_param.o sfx_esm.o
OBJ_F = DRAG.o drag3.o
OBJ = $(OBJ_F90) $(OBJ_F)
diff --git a/param.f90 b/param.f90
deleted file mode 100644
index 4a1519e..0000000
--- a/param.f90
+++ /dev/null
@@ -1,69 +0,0 @@
-module param
- implicit none
- ! acceleration due to gravity [m/s^2]
- real, parameter :: g = 9.81
- ! molecular Prandtl number (air) [n/d]
- real, parameter :: Pr_m = 0.71
- ! kinematic viscosity of air [m^2/s]
- real, parameter :: nu_air = 0.000015e0
- ! 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
- ! inverse Prandtl (turbulent) number in free convection [n/d]
- real, parameter :: Pr_t_inf_inv = 3.5
-
- ! stability function coeff. (unstable) [= g4 & g10 in deprecated code]
- real, parameter :: alpha_m = 16.0
- real, parameter :: alpha_h = 16.0
- real, parameter :: alpha_h_fix = 1.2
-
- ! stability function coeff. (stable)
- real, parameter :: beta_m = 4.7
- real, parameter :: beta_h = beta_m
-
- ! --- max Ri-bulk value in stable case ( < 1 / beta_m )
- real, parameter :: Rib_max = 0.9 / beta_m
-
- ! 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
- real, parameter :: B_max_ocean = 8.0
-
- ! Charnock parameters
- ! z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
- real, parameter :: gamma_c = 0.0144
- real, parameter :: Re_visc_min = 0.111
-
- real, parameter :: h_charnock = 10.0
- real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
- real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
-
-! AN5=(A6/A0)**4
-! D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
-! Y10=(1.0E0-G4*D1)**.25E0
-! X10=(1.0E0-G10*D1)**.5E0
-! P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
-! P0=ALOG((X10-1.0E0)/(X10+1.0E0))
-
-end module PARAM
-
-
-
-
-
-
-
-
-
diff --git a/sfx_esm.f90 b/sfx_esm.f90
new file mode 100644
index 0000000..188ab38
--- /dev/null
+++ b/sfx_esm.f90
@@ -0,0 +1,593 @@
+module sfx_esm
+ !> @brief main Earth System Model surface flux module
+
+ ! macro defs.
+ ! --------------------------------------------------------------------------------
+#define SFX_FORCE_DEPRECATED_CODE
+ ! --------------------------------------------------------------------------------
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_esm_param
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_surface_fluxes
+ public :: get_surface_fluxes_vec
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: meteoDataType
+ !> @brief meteorological input for surface flux calculation
+
+ real :: h !> constant flux layer height [m]
+ real :: U !> abs(wind speed) at 'h' [m/s]
+ real :: dT !> difference between potential temperature at 'h' and at surface [K]
+ real :: Tsemi !> semi-sum of potential temperature at 'h' and at surface [K]
+ real :: dQ !> difference between humidity at 'h' and at surface [g/g]
+ real :: z0_m !> surface aerodynamic roughness (should be < 0 for water bodies surface)
+ end type
+
+ type, public :: meteoDataVecType
+ !> @brief meteorological input for surface flux calculation
+ !> &details using arrays as input
+
+ real, allocatable :: h(:) !> constant flux layer height [m]
+ real, allocatable :: U(:) !> abs(wind speed) at 'h' [m/s]
+ real, allocatable :: dT(:) !> difference between potential temperature at 'h' and at surface [K]
+ real, allocatable :: Tsemi(:) !> semi-sum of potential temperature at 'h' and at surface [K]
+ real, allocatable :: dQ(:) !> difference between humidity at 'h' and at surface [g/g]
+ real, allocatable :: z0_m(:) !> surface aerodynamic roughness (should be < 0 for water bodies surface)
+ end type
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: sfxDataType
+ !> @brief surface flux output data
+
+ real :: zeta !> = z/L [n/d]
+ real :: Rib !> bulk Richardson number [n/d]
+ real :: Re !> Reynolds number [n/d]
+ real :: B !> = log(z0_m / z0_h) [n/d]
+ real :: z0_m !> aerodynamic roughness [m]
+ real :: z0_t !> thermal roughness [m]
+ real :: Rib_conv_lim !> Ri-bulk convection critical value [n/d]
+ real :: Cm !> transfer coefficient for momentum [n/d]
+ real :: Ct !> transfer coefficient for heat [n/d]
+ real :: Km !> eddy viscosity coeff. at h [m^2/s]
+ real :: Pr_t_inv !> inverse turbulent Prandtl number at h [n/d]
+ end type
+
+ type, public :: sfxDataVecType
+ !> @brief surface flux output data
+ !> &details using arrays as output
+
+ real, allocatable :: zeta(:) !> = z/L [n/d]
+ real, allocatable :: Rib(:) !> bulk Richardson number [n/d]
+ real, allocatable :: Re(:) !> Reynolds number [n/d]
+ real, allocatable :: B(:) !> = log(z0_m / z0_h) [n/d]
+ real, allocatable :: z0_m(:) !> aerodynamic roughness [m]
+ real, allocatable :: z0_t(:) !> thermal roughness [m]
+ real, allocatable :: Rib_conv_lim(:) !> Ri-bulk convection critical value [n/d]
+ real, allocatable :: Cm(:) !> transfer coefficient for momentum [n/d]
+ real, allocatable :: Ct(:) !> transfer coefficient for heat [n/d]
+ real, allocatable :: Km(:) !> eddy viscosity coeff. at h [m^2/s]
+ real, allocatable :: Pr_t_inv(:) !> inverse turbulent Prandtl number at h [n/d]
+ end type
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: numericsType
+ integer :: maxiters_convection = 10 !> maximum (actual) number of iterations in convection
+ integer :: maxiters_charnock = 10 !> maximum (actual) number of iterations in charnock roughness
+ end type
+ ! --------------------------------------------------------------------------------
+
+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
+#ifdef USE_DEPRECATED_CODE
+#else
+ meteo_cell = meteoDataType(&
+ h = meteo%h(i), &
+ U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
+ z0_m = meteo%z0_m(i))
+
+ call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
+
+ call push_sfx_data(sfx, sfx_cell, i)
+#endif
+ end do
+
+ end subroutine get_surface_fluxes_vec
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ subroutine push_sfx_data(sfx, sfx_cell, idx)
+ !> @brief helper subroutine for copying data in sfxDataVecType
+ ! ----------------------------------------------------------------------------
+ type (sfxDataVecType), intent(inout) :: sfx
+ type (sfxDataType), intent(in) :: sfx_cell
+ integer, intent(in) :: idx
+ ! ----------------------------------------------------------------------------
+
+ sfx%zeta(idx) = sfx_cell%zeta
+ sfx%Rib(idx) = sfx_cell%Rib
+ sfx%Re(idx) = sfx_cell%Re
+ sfx%B(idx) = sfx_cell%B
+ sfx%z0_m(idx) = sfx_cell%z0_m
+ sfx%z0_t(idx) = sfx_cell%z0_t
+ sfx%Rib_conv_lim(idx) = sfx_cell%Rib_conv_lim
+ sfx%Cm(idx) = sfx_cell%Cm
+ sfx%Ct(idx) = sfx_cell%Ct
+ sfx%Km(idx) = sfx_cell%Km
+ sfx%Pr_t_inv(idx) = sfx_cell%Pr_t_inv
+
+ end subroutine push_sfx_data
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ subroutine get_surface_fluxes(sfx, meteo, numerics)
+ !> @brief surface flux calculation for single cell
+ !> @details contains C/C++ interface
+ ! ----------------------------------------------------------------------------
+ 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)
+ ! ----------------------------------------------------------------------------
+
+ ! --- 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 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 fval !> just a shortcut for partial calculations
+ ! ----------------------------------------------------------------------------
+
+ ! --- shadowing names for clarity
+ U = meteo%U
+ Tsemi = meteo%Tsemi
+ 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 B = log(z0_m / z0_h)
+ Re = u_dyn0 * z0_m / nu_air
+ if(Re <= Re_rough_min) then
+ B = B1_rough * alog(B3_rough * Re) + B2_rough
+ else
+ ! *: B4 takes into account Re value at z' ~ O(10) z0
+ B = B4_rough * (Re**B2_rough)
+ end if
+
+ ! --- apply max restriction based on surface type
+ if (surface_type == surface_ocean) then
+ B = min(B, B_max_ocean)
+ endif
+ if (surface_type == surface_land) then
+ B = min(B, B_max_land)
+ end if
+
+ ! --- define roughness [thermal]
+ z0_t = z0_m / exp(B)
+ ! --- define relative height [thermal]
+ h0_t = h / z0_t
+
+ ! --- 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_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
+
+ fval = beta_m * zeta
+ phi_m = 1.0 + fval
+ phi_h = 1.0/Pr_t_0_inv + fval
+
+ else if (Rib < Rib_conv_lim) then
+ ! --- strong instability block
+
+ call get_psi_convection(psi_m, psi_h, zeta, Rib, &
+ zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection)
+
+ fval = (zeta_conv_lim / zeta)**(1.0/3.0)
+ phi_m = fval / f_m_conv_lim
+ phi_h = fval / (Pr_t_0_inv * f_h_conv_lim)
+
+ else if (Rib > -0.001) then
+ ! --- nearly neutral [-0.001, 0] block
+
+ call get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+
+ phi_m = 1.0
+ phi_h = 1.0 / Pr_t_0_inv
+ else
+ ! --- weak & semistrong instability block
+
+ call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)
+
+ phi_m = (1.0 - alpha_m * zeta)**(-0.25)
+ phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
+ end if
+
+ ! --- define transfer coeff. (momentum) & (heat)
+ Cm = kappa / psi_m
+ Ct = kappa / psi_h
+
+ ! --- 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)
+
+ end subroutine get_surface_fluxes
+ ! --------------------------------------------------------------------------------
+
+ ! 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
+ ! --------------------------------------------------------------------------------
+
+ ! universal functions
+ ! --------------------------------------------------------------------------------
+ subroutine get_psi_neutral(psi_m, psi_h, zeta, h0_m, h0_t, B)
+ !> @brief universal functions (momentum) & (heat): neutral case
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: psi_m, psi_h !> universal functions
+ real, intent(out) :: zeta !> = z/L
+
+ real, intent(in) :: h0_m, h0_t !> = z/z0_m, z/z0_h
+ real, intent(in) :: B !> = log(z0_m / z0_h)
+ ! ----------------------------------------------------------------------------
+
+ zeta = 0.0
+ psi_m = log(h0_m)
+ psi_h = log(h0_t) / Pr_t_0_inv
+ if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv
+
+ end subroutine
+
+ subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
+ !> @brief universal functions (momentum) & (heat): stable case
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
+ real, intent(out) :: zeta !> = z/L [n/d]
+
+ real, intent(in) :: Rib !> bulk Richardson number [n/d]
+ 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 :: Rib_coeff
+ real :: psi0_m, psi0_h
+ real :: phi
+ real :: c
+ ! ----------------------------------------------------------------------------
+
+ psi0_m = alog(h0_m)
+ psi0_h = B / psi0_m
+
+ Rib_coeff = beta_m * Rib
+ c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
+ zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
+ (2.0 * beta_m * (1.0 - Rib_coeff))
+
+ phi = beta_m * zeta
+
+ psi_m = psi0_m + phi
+ psi_h = (psi0_m + B) / Pr_t_0_inv + phi
+
+ end subroutine
+
+ subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
+ !> @brief universal functions (momentum) & (heat): semi-strong convection case
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
+ real, intent(out) :: zeta !> = z/L [n/d]
+
+ real, intent(in) :: Rib !> bulk Richardson number [n/d]
+ 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]
+ integer, intent(in) :: maxiters !> maximum number of iterations
+
+ ! --- local variables
+ real :: zeta0_m, zeta0_h
+ real :: f0_m, f0_h
+ real :: f_m, f_h
+
+ integer :: i
+ ! ----------------------------------------------------------------------------
+
+ psi_m = log(h0_m)
+ psi_h = log(h0_t)
+ if (abs(B) < 1.0e-10) psi_h = psi_m
+
+ zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h
+
+ do i = 1, maxiters + 1
+ zeta0_m = zeta / h0_m
+ zeta0_h = zeta / h0_t
+ if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
+
+ f_m = (1.0 - alpha_m * zeta)**0.25e0
+ f_h = sqrt(1.0 - alpha_h_fix * zeta)
+
+ f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0
+ f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
+
+ f0_m = max(f0_m, 1.000001e0)
+ f0_h = max(f0_h, 1.000001e0)
+
+ psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m))
+ psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / Pr_t_0_inv
+
+ ! *: don't update zeta = z/L at last iteration
+ if (i == maxiters + 1) exit
+
+ zeta = Rib * psi_m**2 / psi_h
+ end do
+
+ end subroutine
+
+ subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
+ zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, &
+ h0_m, h0_t, B, maxiters)
+ !> @brief universal functions (momentum) & (heat): fully convective case
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: psi_m, psi_h !> universal functions [n/d]
+ real, intent(out) :: zeta !> = z/L [n/d]
+
+ real, intent(in) :: Rib !> bulk Richardson number [n/d]
+ 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]
+ integer, intent(in) :: maxiters !> maximum number of iterations
+
+ real, intent(in) :: zeta_conv_lim !> convective limit zeta
+ real, intent(in) :: f_m_conv_lim, f_h_conv_lim !> universal function shortcuts in limiting case
+ ! ----------------------------------------------------------------------------
+
+ ! --- local variables
+ real :: zeta0_m, zeta0_h
+ real :: f0_m, f0_h
+ real :: p_m, p_h
+ real :: a_m, a_h
+ real :: c_lim, f
+
+ integer :: i
+ ! ----------------------------------------------------------------------------
+
+ p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0))
+ p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0))
+
+ zeta = zeta_conv_lim
+ do i = 1, maxiters + 1
+ zeta0_m = zeta / h0_m
+ zeta0_h = zeta / h0_t
+ if (abs(B) < 1.0e-10) zeta0_h = zeta0_m
+
+ f0_m = (1.0 - alpha_m * zeta0_m)**0.25
+ f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)
+
+ a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0))
+ a_h = log((f0_h + 1.0)/(f0_h - 1.0))
+
+ c_lim = (zeta_conv_lim / zeta)**(1.0 / 3.0)
+ f = 3.0 * (1.0 - c_lim)
+
+ psi_m = f / f_m_conv_lim + p_m + a_m
+ psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv
+
+ ! *: don't update zeta = z/L at last iteration
+ if (i == maxiters + 1) exit
+
+ zeta = Rib * psi_m**2 / psi_h
+ end do
+
+ end subroutine
+ ! --------------------------------------------------------------------------------
+
+ ! convection limit definition
+ ! --------------------------------------------------------------------------------
+ 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
+ ! --------------------------------------------------------------------------------
+
+ ! charnock roughness definition
+ ! --------------------------------------------------------------------------------
+ subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: z0_m !> aerodynamic roughness [m]
+ real, intent(out) :: u_dyn0 !> dynamic velocity in neutral conditions [m/s]
+
+ real, intent(in) :: h !> constant flux layer height [m]
+ real, intent(in) :: U !> abs(wind speed) [m/s]
+ integer, intent(in) :: maxiters !> maximum number of iterations
+ ! ----------------------------------------------------------------------------
+
+ ! --- local variables
+ real :: Uc ! wind speed at h_charnock [m/s]
+
+ real :: a, b, c, c_min
+ real :: f
+
+ integer :: i, j
+ ! ----------------------------------------------------------------------------
+
+ Uc = U
+ a = 0.0
+ b = 25.0
+ c_min = log(h_charnock) / kappa
+
+ do i = 1, maxiters
+ f = c1_charnock - 2.0 * log(Uc)
+ do j = 1, maxiters
+ c = (f + 2.0 * log(b)) / kappa
+ ! looks like the check should use U10 instead of U
+ ! but note that a1 should be set = 0 in cycle beforehand
+ if (U <= 8.0e0) a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
+ c = max(c - a, c_min)
+ b = c
+ end do
+ z0_m = h_charnock * exp(-c * kappa)
+ z0_m = max(z0_m, 0.000015e0)
+ Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
+ end do
+
+ ! --- define dynamic velocity in neutral conditions
+ u_dyn0 = Uc / c
+
+ end subroutine
+ ! --------------------------------------------------------------------------------
+
+end module sfx_esm
\ No newline at end of file
diff --git a/sfx_esm_param.f90 b/sfx_esm_param.f90
new file mode 100644
index 0000000..fad0d75
--- /dev/null
+++ b/sfx_esm_param.f90
@@ -0,0 +1,71 @@
+module sfx_esm_param
+ !> @brief ESM surface flux model parameters
+ !> @details all in SI units
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+ integer, public, parameter :: surface_ocean = 0 !> ocean surface
+ integer, public, parameter :: surface_land = 1 !> land surface
+
+ !> 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
+ !> inverse Prandtl (turbulent) number in free convection [n/d]
+ real, parameter :: Pr_t_inf_inv = 3.5
+
+ !> stability function coeff. (unstable) [= g4 & g10 in deprecated code]
+ real, parameter :: alpha_m = 16.0
+ real, parameter :: alpha_h = 16.0
+ real, parameter :: alpha_h_fix = 1.2
+
+ !> stability function coeff. (stable)
+ real, parameter :: beta_m = 4.7
+ real, parameter :: beta_h = beta_m
+
+ !> --- max Ri-bulk value in stable case ( < 1 / beta_m )
+ real, parameter :: Rib_max = 0.9 / beta_m
+
+ !> 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
+ real, parameter :: B_max_ocean = 8.0
+
+ !> Charnock parameters
+ !> z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
+ real, parameter :: gamma_c = 0.0144
+ real, parameter :: Re_visc_min = 0.111
+
+ real, parameter :: h_charnock = 10.0
+ real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
+ real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
+
+end module sfx_esm_param
+
+
+
+
+
+
+
+
+
diff --git a/main_drag.f90 b/sfx_main.f90
similarity index 62%
rename from main_drag.f90
rename to sfx_main.f90
index 776bee4..40ac4b8 100644
--- a/main_drag.f90
+++ b/sfx_main.f90
@@ -1,8 +1,9 @@
PROGRAM main_ddt
-
- USE param
+
+ use sfx_phys_const
+ USE sfx_esm_param
USE inputdata
- USE drag3
+ USE sfx_esm
type (meteoDataType):: data_in1
@@ -11,7 +12,6 @@
type (numericsType) :: data_par1
- type (data_lutyp) :: data_lutyp1
integer :: numst, i
real :: cflh, z0in
character(len = 50) :: filename_in
@@ -39,21 +39,21 @@
! lu_indx - 1 for land, 2 for sea, 3 for lake
! test - file input
- type :: datatype_outMAS1
- real, allocatable :: masout_zl(:)
- real, allocatable :: masout_ri(:)
- real, allocatable :: masout_re(:)
- real, allocatable :: masout_lnzuzt(:)
- real, allocatable :: masout_zu(:)
- real, allocatable :: masout_ztout(:)
- real, allocatable :: masout_rith(:)
- real, allocatable :: masout_cm(:)
- real, allocatable :: masout_ch(:)
- real, allocatable :: masout_ct(:)
- real, allocatable :: masout_ckt(:)
- end type
-
- type(datatype_outMAS1) :: data_outMAS
+ !type :: datatype_outMAS1
+ ! real, allocatable :: masout_zl(:)
+ ! real, allocatable :: masout_ri(:)
+ ! real, allocatable :: masout_re(:)
+ ! real, allocatable :: masout_lnzuzt(:)
+ ! real, allocatable :: masout_zu(:)
+ ! real, allocatable :: masout_ztout(:)
+ ! real, allocatable :: masout_rith(:)
+ ! real, allocatable :: masout_cm(:)
+ ! real, allocatable :: masout_ch(:)
+ ! real, allocatable :: masout_ct(:)
+ ! real, allocatable :: masout_ckt(:)
+ !end type
+
+ type(sfxDataVecType) :: data_outMAS
!output
!masout_zl - non-dimensional cfl hight
@@ -101,18 +101,17 @@
allocate(meteo%z0_m(numst))
-
- allocate(data_outMAS%masout_zl(numst))
- allocate(data_outMAS%masout_ri(numst))
- allocate(data_outMAS%masout_re(numst))
- allocate(data_outMAS%masout_lnzuzt(numst))
- allocate(data_outMAS%masout_zu(numst))
- allocate(data_outMAS%masout_ztout(numst))
- allocate(data_outMAS%masout_rith(numst))
- allocate(data_outMAS%masout_cm(numst))
- allocate(data_outMAS%masout_ch(numst))
- allocate(data_outMAS%masout_ct(numst))
- allocate(data_outMAS%masout_ckt(numst))
+ allocate(data_outMAS%zeta(numst))
+ allocate(data_outMAS%Rib(numst))
+ allocate(data_outMAS%Re(numst))
+ allocate(data_outMAS%B(numst))
+ allocate(data_outMAS%z0_m(numst))
+ allocate(data_outMAS%z0_t(numst))
+ allocate(data_outMAS%Rib_conv_lim(numst))
+ allocate(data_outMAS%Cm(numst))
+ allocate(data_outMAS%Ct(numst))
+ allocate(data_outMAS%Km(numst))
+ allocate(data_outMAS%Pr_t_inv(numst))
open (11, file=filename_in2, status ='old')
open (1, file= filename_in, status ='old')
@@ -128,11 +127,11 @@
meteo%z0_m(i)=z0in
enddo
- CALL surf_flux_vec(meteo, &
- data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,&
- data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,&
- data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,&
- data_par1, data_lutyp1,numst)
+ CALL get_surface_fluxes_vec(data_outMAS, meteo, &
+ !data_outMAS%zeta, data_outMAS%Rib, data_outMAS%Re, data_outMAS%B,&
+ !data_outMAS%z0_m,data_outMAS%z0_t,data_outMAS%Rib_conv_lim,data_outMAS%Cm,&
+ !data_outMAS%Ct,data_outMAS%Km,data_outMAS%Pr_t_inv,&
+ data_par1, numst)
!CALL surf_fluxMAS(meteo%mas_w, meteo%mas_dt, meteo%mas_st, meteo%mas_dq,&
! meteo%mas_cflh, meteo%mas_z0in,&
@@ -142,9 +141,9 @@
! data_par1, data_lutyp1,numst)
do i=1,numst
- write (2,20) data_outMAS%masout_zl(i), data_outMAS%masout_ri(i), data_outMAS%masout_re(i), data_outMAS%masout_lnzuzt(i),&
- data_outMAS%masout_zu(i), data_outMAS%masout_ztout(i), data_outMAS%masout_rith(i), data_outMAS%masout_cm(i),&
- data_outMAS%masout_ch(i), data_outMAS%masout_ct(i), data_outMAS%masout_ckt(i)
+ write (2,20) data_outMAS%zeta(i), data_outMAS%Rib(i), data_outMAS%Re(i), data_outMAS%B(i),&
+ data_outMAS%z0_m(i), data_outMAS%z0_t(i), data_outMAS%Rib_conv_lim(i), data_outMAS%Cm(i),&
+ data_outMAS%Ct(i), data_outMAS%Km(i), data_outMAS%Pr_t_inv(i)
enddo
@@ -155,18 +154,17 @@
deallocate(meteo%dQ)
deallocate(meteo%z0_m)
-
- deallocate(data_outMAS%masout_zl)
- deallocate(data_outMAS%masout_ri)
- deallocate(data_outMAS%masout_re)
- deallocate(data_outMAS%masout_lnzuzt)
- deallocate(data_outMAS%masout_zu)
- deallocate(data_outMAS%masout_ztout)
- deallocate(data_outMAS%masout_rith)
- deallocate(data_outMAS%masout_cm)
- deallocate(data_outMAS%masout_ch)
- deallocate(data_outMAS%masout_ct)
- deallocate(data_outMAS%masout_ckt)
+ deallocate(data_outMAS%zeta)
+ deallocate(data_outMAS%Rib)
+ deallocate(data_outMAS%Re)
+ deallocate(data_outMAS%B)
+ deallocate(data_outMAS%z0_m)
+ deallocate(data_outMAS%z0_t)
+ deallocate(data_outMAS%Rib_conv_lim)
+ deallocate(data_outMAS%Cm)
+ deallocate(data_outMAS%Ct)
+ deallocate(data_outMAS%Km)
+ deallocate(data_outMAS%Pr_t_inv)
10 format (f8.4,2x,f8.4)
diff --git a/sfx_phys_const.f90 b/sfx_phys_const.f90
new file mode 100644
index 0000000..8c8ec1c
--- /dev/null
+++ b/sfx_phys_const.f90
@@ -0,0 +1,12 @@
+module sfx_phys_const
+ !> @brief general physics constants
+ !> @details all in SI units
+ implicit none
+ public
+
+ real, parameter :: g = 9.81 !> acceleration due to gravity [m/s^2]
+
+ real, parameter :: nu_air = 0.000015e0 !> kinematic viscosity of air [m^2/s]
+ real, parameter :: Pr_m = 0.71 !> molecular Prandtl number (air) [n/d]
+
+end module sfx_phys_const
--
GitLab