From 1de9b15ad6d12d22814d9ef4c8c18eea76bf03df Mon Sep 17 00:00:00 2001
From: Evgeny Mortikov <evgeny.mortikov@gmail.com>
Date: Mon, 18 Dec 2023 04:39:00 +0300
Subject: [PATCH] separating 'surface' module for aerodynamics and thermal
roughness calculation
---
compile.sh | 4 +-
makefile | 2 +-
srcF/sfx_data.f90 | 4 --
srcF/sfx_def.fi | 2 +-
srcF/sfx_esm.f90 | 23 ++-----
srcF/sfx_esm_param.f90 | 19 ------
srcF/sfx_log.f90 | 21 +-----
srcF/sfx_log_param.f90 | 19 ------
srcF/{sfx_roughness.f90 => sfx_surface.f90} | 73 +++++++++++++++++++--
9 files changed, 79 insertions(+), 88 deletions(-)
rename srcF/{sfx_roughness.f90 => sfx_surface.f90} (56%)
diff --git a/compile.sh b/compile.sh
index c04fab7..146d20b 100755
--- a/compile.sh
+++ b/compile.sh
@@ -4,7 +4,7 @@ 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_surface.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_log_param.f90
gfortran -c -cpp -Wuninitialized srcF/sfx_log.f90
@@ -13,5 +13,5 @@ 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_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
+gfortran -o sfx.exe sfx_main.o sfx_log.o sfx_log_param.o sfx_esm.o sfx_esm_param.o sfx_surface.o sfx_data.o sfx_common.o sfx_phys_const.o
rm *.o *.mod
diff --git a/makefile b/makefile
index 913ada0..f940f58 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_data.o sfx_roughness.o sfx_log_param.o sfx_log.o sfx_esm_param.o sfx_esm.o sfx_main.o
+OBJ_F90 = sfx_phys_const.o sfx_common.o sfx_data.o sfx_surface.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
index 5e252d7..c13ddc2 100644
--- a/srcF/sfx_data.f90
+++ b/srcF/sfx_data.f90
@@ -1,10 +1,6 @@
module sfx_data
!> @brief surface flux module data
- ! macro defs.
- ! --------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------
-
! modules used
! --------------------------------------------------------------------------------
! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_def.fi b/srcF/sfx_def.fi
index 6d90e99..2fc04b6 100644
--- a/srcF/sfx_def.fi
+++ b/srcF/sfx_def.fi
@@ -1,4 +1,4 @@
! sfx model macro definitions
-!#define SFX_FORCE_DEPRECATED_CODE
+!#define SFX_FORCE_DEPRECATED_ESM_CODE
!#define SFX_CHECK_NAN
\ No newline at end of file
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index dd38ce4..d1fef61 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -9,7 +9,7 @@ module sfx_esm
use sfx_common
#endif
use sfx_data
- use sfx_roughness
+ use sfx_surface
use sfx_esm_param
! --------------------------------------------------------------------------------
@@ -53,7 +53,7 @@ contains
! ----------------------------------------------------------------------------
do i = 1, n
-#ifdef SFX_FORCE_DEPRECATED_CODE
+#ifdef SFX_FORCE_DEPRECATED_ESM_CODE
#else
meteo_cell = meteoDataType(&
h = meteo%h(i), &
@@ -171,25 +171,10 @@ contains
u_dyn0 = U * kappa / log(h0_m)
end if
- ! --- define B = log(z0_m / z0_h)
+ ! --- define thermal roughness & 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
+ call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
- ! --- define roughness [thermal]
- z0_t = z0_m / exp(B)
! --- define relative height [thermal]
h0_t = h / z0_t
diff --git a/srcF/sfx_esm_param.f90 b/srcF/sfx_esm_param.f90
index 0f9e253..76003ed 100644
--- a/srcF/sfx_esm_param.f90
+++ b/srcF/sfx_esm_param.f90
@@ -12,9 +12,6 @@ module sfx_esm_param
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
@@ -35,20 +32,4 @@ module sfx_esm_param
!> --- 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
-
end module sfx_esm_param
diff --git a/srcF/sfx_log.f90 b/srcF/sfx_log.f90
index 225eb37..02e93b0 100644
--- a/srcF/sfx_log.f90
+++ b/srcF/sfx_log.f90
@@ -9,7 +9,7 @@ module sfx_log
use sfx_common
#endif
use sfx_data
- use sfx_roughness
+ use sfx_surface
use sfx_log_param
! --------------------------------------------------------------------------------
@@ -160,25 +160,10 @@ contains
u_dyn0 = U * kappa / log(h0_m)
end if
- ! --- define B = log(z0_m / z0_h)
+ ! --- define thermal roughness & 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
+ call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
- ! --- define roughness [thermal]
- z0_t = z0_m / exp(B)
! --- define relative height [thermal]
h0_t = h / z0_t
diff --git a/srcF/sfx_log_param.f90 b/srcF/sfx_log_param.f90
index 3918880..16d6e8d 100644
--- a/srcF/sfx_log_param.f90
+++ b/srcF/sfx_log_param.f90
@@ -12,29 +12,10 @@ module sfx_log_param
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_roughness.f90 b/srcF/sfx_surface.f90
similarity index 56%
rename from srcF/sfx_roughness.f90
rename to srcF/sfx_surface.f90
index e5e5463..628626f 100644
--- a/srcF/sfx_roughness.f90
+++ b/srcF/sfx_surface.f90
@@ -1,6 +1,6 @@
#include "sfx_def.fi"
-module sfx_roughness
+module sfx_surface
!> @brief surface roughness parameterizations
! modules used
@@ -14,13 +14,20 @@ module sfx_roughness
private
! --------------------------------------------------------------------------------
+ ! public interface
! --------------------------------------------------------------------------------
- real, parameter :: kappa = 0.40 !> von Karman constant [n/d]
+ public :: get_charnock_roughness
+ public :: get_thermal_roughness
! --------------------------------------------------------------------------------
- ! public interface
+
+ ! *: add surface type as input value
+ integer, public, parameter :: surface_ocean = 0 !> ocean surface
+ integer, public, parameter :: surface_land = 1 !> land surface
+ integer, public, parameter :: surface_lake = 2 !> lake surface
+
! --------------------------------------------------------------------------------
- public :: get_charnock_roughness
+ real, parameter :: kappa = 0.40 !> von Karman constant [n/d]
! --------------------------------------------------------------------------------
!> Charnock parameters
@@ -34,6 +41,23 @@ module sfx_roughness
real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
! --------------------------------------------------------------------------------
+ !> 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
+ real, parameter :: B_max_lake = 8.0
+
contains
! charnock roughness definition
@@ -83,4 +107,43 @@ contains
end subroutine
! --------------------------------------------------------------------------------
-end module sfx_roughness
+ ! thermal roughness definition
+ ! --------------------------------------------------------------------------------
+ subroutine get_thermal_roughness(z0_t, B, &
+ z0_m, Re, surface_type)
+ ! ----------------------------------------------------------------------------
+ real, intent(out) :: z0_t !> thermal roughness [m]
+ real, intent(out) :: B !> = log(z0_m / z0_t) [n/d]
+
+ real, intent(in) :: z0_m !> aerodynamic roughness [m]
+ real, intent(in) :: Re !> roughness Reynolds number [n/d]
+ integer, intent(in) :: surface_type !> = [ocean] || [land] || [lake]
+ ! ----------------------------------------------------------------------------
+
+ ! --- local variables
+ ! ----------------------------------------------------------------------------
+
+ ! --- define B = log(z0_m / z0_t)
+ if (Re <= Re_rough_min) then
+ B = B1_rough * alog(B3_rough * Re) + B2_rough
+ else
+ ! *: B4 takes into account Re value at z' ~ O(10) z0
+ B = B4_rough * (Re**B2_rough)
+ end if
+
+ ! --- apply max restriction based on surface type
+ if (surface_type == surface_ocean) then
+ B = min(B, B_max_ocean)
+ else if (surface_type == surface_lake) then
+ B = min(B, B_max_lake)
+ else if (surface_type == surface_land) then
+ B = min(B, B_max_land)
+ end if
+
+ ! --- define roughness [thermal]
+ z0_t = z0_m / exp(B)
+
+ end subroutine
+ ! --------------------------------------------------------------------------------
+
+end module sfx_surface
--
GitLab