turbulence closures implementation update

965ff2fc · Evgeny Mortikov · 95c99b68 · 965ff2fc · 965ff2fc · 965ff2fc
Commit 965ff2fc authored 7 months ago by Evgeny Mortikov
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,7 @@ set(SOURCES
    obl_k_epsilon.f90
    obl_pph.f90
    obl_pph_dyn.f90
+    obl_turb.f90
    obl_scm.f90
    obl_main.f90
    io_metadata.f90

--- a/obl_k_epsilon.f90
+++ b/obl_k_epsilon.f90
@@ -2,9 +2,9 @@ module obl_k_epsilon
    !< @brief k-epsilon scheme module
    
    ! modules used
-    use obl_turb_closure
-    use obl_common_phys_parameters
    use obl_math
+    use obl_grid
+    use obl_boundary

    ! directives list
    implicit none
@@ -25,6 +25,36 @@ module obl_k_epsilon

    contains

+    ! --------------------------------------------------------------------------------
+    subroutine advance_k_epsilon(param, bc, grid, dt)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(kepsilonParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+        real, intent(in) :: dt
+        ! --------------------------------------------------------------------------------
+
+        call TKE_bc(TKE, &
+            bc%U_dyn0, bc%U_dynH, param, grid%cz)
+        call EPS_bc(EPS, &
+            bc%U_dyn0, bc%U_dynH, param, grid%dz, grid%cz)
+        call get_Km_k_eps (Km, TKE, EPS, param, grid%cz)
+        call get_Kh_k_eps (Kh, Km, param, grid%cz)
+        call TKE_bc(TKE, &
+            bc%U_dyn0, bc%U_dynH, param, grid%cz)
+        call eps_bc(EPS, &
+            bc%U_dyn0, bc%U_dynH, param, grid%dz, grid%cz)
+        !call solve_TKE_eq_semiimplicit (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps, param)
+        call solve_TKE_eq (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS, param)
+        call limit_min_array(TKE, param%TKE_min, grid%cz)
+        call solve_eps_eq_semiimplicit (EPS, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE, param)
+        call limit_min_array(EPS, param%eps_min, grid%cz) 
+
+    end subroutine
+
    subroutine TKE_init (TKE, param, nz)
        !< @brief TKE initialization
        type(kepsilonParamType), intent(in) :: param

--- a/obl_main.f90
+++ b/obl_main.f90
@@ -373,31 +373,11 @@ program obl_main
        call get_TKE_buoyancy_vec(TKE_buoyancy, Kh, N2, grid%cz)
        
        if (closure_mode.eq.1) then
-            call get_eddy_viscosity(Km, Ri_grad, param_pph, grid)
-            call get_eddy_diffusivity(Kh, Ri_grad, param_pph, grid)
-
+            call advance_pph(param_pph, bc, grid, dt)
        else if (closure_mode.eq.2) then
-            call get_Km_plus(Km, Ri_grad, &
-                bc%U_dynH, mld, param_pph_dyn, grid%dz, grid%height, grid%cz)
-            call get_Kh_plus(Kh, Km, param_pph_dyn, grid%cz)
-        
+            call advance_pph_dyn(param_pph_dyn, bc, grid, dt)
        else if (closure_mode.eq.4) then
-            call TKE_bc(TKE, &
-                bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%cz)
-            call EPS_bc(EPS, &
-                bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%dz, grid%cz)
-            call get_Km_k_eps (Km, TKE, EPS, param_k_epsilon, grid%cz)
-            call get_Kh_k_eps (Kh, Km, param_k_epsilon, grid%cz)
-            call TKE_bc(TKE, &
-                bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%cz)
-            call eps_bc(EPS, &
-                bc%U_dyn0, bc%U_dynH, param_k_epsilon, grid%dz, grid%cz)
-            !call solve_TKE_eq_semiimplicit (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, eps, param_k_epsilon)
-            call solve_TKE_eq (TKE, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, EPS, param_k_epsilon)
-            call limit_min_array(TKE, param_k_epsilon%TKE_min, grid%cz)
-            call solve_eps_eq_semiimplicit (EPS, Km, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE, param_k_epsilon)
-            call limit_min_array(EPS, param_k_epsilon%eps_min, grid%cz) 
-        
+            call advance_k_epsilon(param_k_epsilon, bc, grid, dt)
        endif
        ! ----------------------------------------------------------------------------


--- a/obl_pph.f90
+++ b/obl_pph.f90
@@ -3,6 +3,7 @@ module obl_pph

    ! modules used
    use obl_grid
+    use obl_boundary

    ! directives list
    implicit none
@@ -10,6 +11,7 @@ module obl_pph

    ! public interface
    ! --------------------------------------------------------------------------------
+    public :: advance_pph
    public :: get_eddy_viscosity, get_eddy_diffusivity
    public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
    ! --------------------------------------------------------------------------------
@@ -26,6 +28,22 @@ module obl_pph

    contains

+    ! --------------------------------------------------------------------------------
+    subroutine advance_pph(param, bc, grid, dt)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(pphParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+        real, intent(in) :: dt
+        ! --------------------------------------------------------------------------------
+
+        call get_eddy_viscosity(Km, Ri_grad, param, grid)
+        call get_eddy_diffusivity(Kh, Ri_grad, param, grid)
+    end subroutine
+    
    ! --------------------------------------------------------------------------------
    subroutine get_eddy_viscosity(Km, Ri_grad, param, grid)
        !< @brief calculate eddy viscosity on grid

--- a/obl_pph_dyn.f90
+++ b/obl_pph_dyn.f90
@@ -3,6 +3,8 @@ module obl_pph_dyn
    
    ! modules used
    use obl_grid
+    use obl_boundary
+    use obl_diag
    
    ! directives list
    implicit none
@@ -19,6 +21,28 @@ module obl_pph_dyn

    contains

+    ! --------------------------------------------------------------------------------
+    subroutine advance_pph_dyn(param, bc, grid, dt)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(pphDynParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+        real, intent(in) :: dt
+
+        real :: mld
+        ! --------------------------------------------------------------------------------
+
+        call get_mld(mld, N2, grid%dz, grid%cz, grid%height)
+
+        call get_Km_plus(Km, Ri_grad, &
+            bc%U_dynH, mld, param, grid%dz, grid%height, grid%cz)
+        call get_Kh_plus(Kh, Km, param, grid%cz)
+    end subroutine
+
+    ! --------------------------------------------------------------------------------
    subroutine get_Km_plus(Km, Ri_grad, U_dynH, mld, param, dz, z, nz)
        !< @brief calculate momentum transfer coefficient in P-Ph+ closure
        type(pphDynParamType), intent(in) :: param

--- a/obl_turb.f90
+++ b/obl_turb.f90
+module obl_pph
+    !< @brief standard pacanowski-philander scheme
+
+    ! modules used
+    use obl_grid
+
+    ! directives list
+    implicit none
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    !public :: get_eddy_viscosity, get_eddy_diffusivity
+    ! --------------------------------------------------------------------------------
+
+
+    contains
+
+
+end module