diff --git a/compile.sh b/compile.sh
index 107662793561f45fcded8574cf3e2f2f41fa07fb..c04fab77b69f92e2f4813facd9f5a81304bdf331 100755
--- a/compile.sh
+++ b/compile.sh
@@ -2,8 +2,16 @@
gfortran -c -cpp -Wuninitialized srcF/sfx_phys_const.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_common.f90
+gfortran -c -cpp -Wuninitialized srcF/sfx_data.f90
+
+gfortran -c -cpp -Wuninitialized srcF/sfx_roughness.f90
+
+gfortran -c -cpp -Wuninitialized srcF/sfx_log_param.f90
+gfortran -c -cpp -Wuninitialized srcF/sfx_log.f90
+
gfortran -c -cpp -Wuninitialized srcF/sfx_esm_param.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_esm.f90
+
gfortran -c -cpp -Wuninitialized srcF/sfx_main.f90
-gfortran -o sfx.exe sfx_main.o sfx_esm.o sfx_esm_param.o sfx_common.o sfx_phys_const.o
+gfortran -o sfx.exe sfx_main.o sfx_log.o sfx_log_param.o sfx_esm.o sfx_esm_param.o sfx_roughness.o sfx_data.o sfx_common.o sfx_phys_const.o
rm *.o *.mod
diff --git a/makefile b/makefile
index 05e2436e4687416b83d0c1a259264bc626949fed..913ada0ce2fdb39f87657db9b8ec57fcae3ca142 100644
--- a/makefile
+++ b/makefile
@@ -10,7 +10,7 @@ ifeq ($(COMPILER),gnu)
FC = gfortran
endif
-OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_esm_param.o sfx_esm.o sfx_main.o
+OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_data.o sfx_roughness.o sfx_log_param.o sfx_log.o sfx_esm_param.o sfx_esm.o sfx_main.o
OBJ_F =
OBJ = $(OBJ_F90) $(OBJ_F)
diff --git a/srcF/sfx_data.f90 b/srcF/sfx_data.f90
new file mode 100644
index 0000000000000000000000000000000000000000..7ab691c9f4fa73904f1e96290d4856cf6ff25d22
--- /dev/null
+++ b/srcF/sfx_data.f90
@@ -0,0 +1,82 @@
+module sfx_data
+ !> @brief surface flux module data
+
+ ! macro defs.
+ ! --------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ 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
+ ! --------------------------------------------------------------------------------
+
+end module sfx_data
\ No newline at end of file
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index a611bd2e4a159e88e4b914e0b6062e31b4db2833..2633a8481953c3891745f781e034b615491631a8 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -12,6 +12,8 @@ module sfx_esm
#ifdef SFX_CHECK_NAN
use sfx_common
#endif
+ use sfx_data
+ use sfx_roughness
use sfx_esm_param
! --------------------------------------------------------------------------------
@@ -27,66 +29,6 @@ module sfx_esm
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
@@ -569,51 +511,4 @@ contains
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/srcF/sfx_esm_param.f90 b/srcF/sfx_esm_param.f90
index 3aec59df790a95e63a61dec30cc353c1a5ae24cd..0f9e253919ff842374631fd9bb79243efcec1836 100644
--- a/srcF/sfx_esm_param.f90
+++ b/srcF/sfx_esm_param.f90
@@ -51,22 +51,4 @@ module sfx_esm_param
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/srcF/sfx_log.f90 b/srcF/sfx_log.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a3c34655dd813972c3d95d1b4bc92f7697c49a86
--- /dev/null
+++ b/srcF/sfx_log.f90
@@ -0,0 +1,258 @@
+module sfx_log
+ !> @brief simple log-roughness surface flux module
+
+ ! macro defs.
+ ! --------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+#ifdef SFX_CHECK_NAN
+ use sfx_common
+#endif
+ use sfx_data
+ use sfx_roughness
+ use sfx_log_param
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_surface_fluxes
+ public :: get_surface_fluxes_vec
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ type, public :: numericsType
+ 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
+ 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)
+ 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
+ ! ----------------------------------------------------------------------------
+#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)
+ ! ----------------------------------------------------------------------------
+
+ ! --- 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 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)
+
+
+#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_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
+
+ ! --- get the fluxes
+ ! ----------------------------------------------------------------------------
+ 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
+
+ ! --- 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 = 0.0, &
+ Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+
+ end subroutine get_surface_fluxes
+ ! --------------------------------------------------------------------------------
+
+ ! 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
+ ! --------------------------------------------------------------------------------
+
+end module sfx_log
\ No newline at end of file
diff --git a/srcF/sfx_log_param.f90 b/srcF/sfx_log_param.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3918880c9093be206167a72917d920491cd9d42a
--- /dev/null
+++ b/srcF/sfx_log_param.f90
@@ -0,0 +1,40 @@
+module sfx_log_param
+ !> @brief log-roughness surface flux model parameters
+ !> @details all in SI units
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ ! --------------------------------------------------------------------------------
+
+ ! *: add lake surface & add surface type as input value
+ 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
+
+ !> 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
+
+end module sfx_log_param
diff --git a/srcF/sfx_main.f90 b/srcF/sfx_main.f90
index 96a435cf954d1dd8e4535822c0d6b716c994fe0c..e65c55ee3efb7c3bb58736fa52bfa1ead3ddcdfe 100644
--- a/srcF/sfx_main.f90
+++ b/srcF/sfx_main.f90
@@ -4,9 +4,14 @@ program sfx_main
! --------------------------------------------------------------------------------
use sfx_phys_const
use sfx_common
-
- use sfx_esm_param
- use sfx_esm
+ use sfx_data
+
+ use sfx_esm, only: &
+ get_surface_fluxes_vec_esm => get_surface_fluxes_vec, &
+ numericsType_esm => numericsType
+ use sfx_log, only: &
+ get_surface_fluxes_vec_log => get_surface_fluxes_vec, &
+ numericsType_log => numericsType
! --------------------------------------------------------------------------------
! directives list
@@ -15,17 +20,19 @@ program sfx_main
! --------------------------------------------------------------------------------
!> dataset ID:
- !> = 0, USER
- !> = 1, MOSAiC dataset
- !> = 2, IRGASON dataset
- !> = 3, SHEBA dataset
integer :: dataset_id
character(len = 256) :: dataset_name
- integer, parameter :: dataset_MOSAiC = 1
- integer, parameter :: dataset_IRGASON = 2
+ integer, parameter :: dataset_MOSAiC = 1 !> MOSAiC campaign
+ integer, parameter :: dataset_IRGASON = 2 !> IRGASON data
integer, parameter :: dataset_SHEBA = 3 !> please spell 'SHIBA'
integer, parameter :: dataset_USER = 4 !> used defined dataset
+ !> sfx model ID:
+ integer :: model_id
+ character(len = 256) :: model_name
+ integer, parameter :: model_esm = 0 !> ESM model
+ integer, parameter :: model_log = 1 !> LOG simplified model
+
! input/output data
! --------------------------------------------------------------------------------
type(meteoDataVecType) :: meteo !> meteorological data (input)
@@ -33,7 +40,8 @@ program sfx_main
type(sfxDataVecType) :: sfx !> surface fluxes (output)
- type(numericsType) :: numerics !> surface flux module numerics parameters
+ type(numericsType_esm) :: numerics_esm !> surface flux module (ESM) numerics parameters
+ type(numericsType_log) :: numerics_log !> surface flux module (LOG) numerics parameters
integer :: num !> number of 'cells' in input
@@ -49,11 +57,15 @@ program sfx_main
integer :: num_args
character(len = 128) :: arg
+ character(len = 128), parameter :: arg_key_model = '--model'
character(len = 128), parameter :: arg_key_dataset = '--dataset'
character(len = 128), parameter :: arg_key_output = '--output'
character(len = 128), parameter :: arg_key_nmax = '--nmax'
character(len = 128), parameter :: arg_key_help = '--help'
+ character(len = 128), parameter :: arg_key_esm = 'esm'
+ character(len = 128), parameter :: arg_key_log = 'log'
+
character(len = 128), parameter :: arg_key_mosaic = 'mosaic'
character(len = 128), parameter :: arg_key_irgason = 'irgason'
character(len = 128), parameter :: arg_key_sheba = 'sheba'
@@ -70,7 +82,8 @@ program sfx_main
! --------------------------------------------------------------------------------
- !> @brief define dataset
+ !> @brief define model & dataset
+ model_id = model_esm !> default = ESM
dataset_id = dataset_MOSAiC !> default = MOSAiC
is_output_set = 0
@@ -82,6 +95,8 @@ program sfx_main
write(*, *) ' sfx model, usage:'
write(*, *) ' --help '
write(*, *) ' print usage options '
+ write(*, *) ' --model [key]'
+ write(*, *) ' key = esm || log'
write(*, *) ' --dataset [key]'
write(*, *) ' key = mosaic || irgason || sheba || user [files]'
write(*, *) ' files = in-common-file in-file out-file'
@@ -91,6 +106,21 @@ program sfx_main
write(*, *) ' max number of data points > 0 '
stop
end if
+ if (trim(arg) == trim(arg_key_model)) then
+ if (i == num_args) then
+ write(*, *) ' FAILURE! > missing model [key] argument'
+ stop
+ end if
+ call get_command_argument(i + 1, arg)
+ if (trim(arg) == trim(arg_key_esm)) then
+ model_id = model_esm
+ else if (trim(arg) == trim(arg_key_log)) then
+ model_id = model_log
+ else
+ write(*, *) ' FAILURE! > unknown model [key]: ', trim(arg)
+ stop
+ end if
+ end if
if (trim(arg) == trim(arg_key_dataset)) then
if (i == num_args) then
write(*, *) ' FAILURE! > missing dataset [key] argument'
@@ -144,7 +174,17 @@ program sfx_main
end do
- !> @brief set filenames & data for specific dataset
+ !> @brief set name for specific model
+ if (model_id == model_esm) then
+ model_name = "ESM"
+ else if (model_id == model_log) then
+ model_name = "LOG"
+ else
+ write(*, *) ' FAILURE! > unknown model id: ', model_id
+ stop
+ end if
+
+ !> @brief set name & filenames for specific dataset
if (dataset_id == dataset_MOSAiC) then
dataset_name = 'MOSAiC'
@@ -166,9 +206,11 @@ program sfx_main
dataset_name = 'USER'
else
write(*, *) ' FAILURE! > unknown dataset id: ', dataset_id
+ stop
end if
write(*, *) ' Running SFX model'
+ write(*, *) ' model = ', trim(model_name)
write(*, *) ' dataset = ', trim(dataset_name)
write(*, *) ' filename[IN-COMMON] = ', trim(filename_in_common)
write(*, *) ' filename[IN] = ', trim(filename_in)
@@ -239,7 +281,11 @@ program sfx_main
!> @brief calling flux module
- call get_surface_fluxes_vec(sfx, meteo, numerics, num)
+ if (model_id == model_esm) then
+ call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, num)
+ else if (model_id == model_log) then
+ call get_surface_fluxes_vec_log(sfx, meteo, numerics_log, num)
+ end if
!> @brief write output data
diff --git a/srcF/sfx_roughness.f90 b/srcF/sfx_roughness.f90
new file mode 100644
index 0000000000000000000000000000000000000000..8e01cc9438299d1371783caa4ba81b51c12de45d
--- /dev/null
+++ b/srcF/sfx_roughness.f90
@@ -0,0 +1,88 @@
+module sfx_roughness
+ !> @brief surface roughness parameterizations
+
+ ! macro defs.
+ ! --------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------
+
+ ! modules used
+ ! --------------------------------------------------------------------------------
+ use sfx_phys_const
+ ! --------------------------------------------------------------------------------
+
+ ! directives list
+ ! --------------------------------------------------------------------------------
+ implicit none
+ private
+ ! --------------------------------------------------------------------------------
+
+ ! --------------------------------------------------------------------------------
+ real, parameter :: kappa = 0.40 !> von Karman constant [n/d]
+ ! --------------------------------------------------------------------------------
+
+ ! public interface
+ ! --------------------------------------------------------------------------------
+ public :: get_charnock_roughness
+ ! --------------------------------------------------------------------------------
+
+ !> 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
+ ! --------------------------------------------------------------------------------
+
+contains
+
+ ! 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_roughness