moving input/output data, roughness parameterizations in separate modules;...

moving input/output data, roughness parameterizations in separate modules; adding options to choose surface flux model

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

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)
 

srcF/sfx_data.f90

+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

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

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
-
-
-
-
-
-
-
-
-

srcF/sfx_log.f90

+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

srcF/sfx_log_param.f90

+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

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

srcF/sfx_roughness.f90

+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