minor k-epsilon update

obl_k_epsilon.f90

 module obl_k_epsilon
     !< @brief k-epsilon scheme module
+    ! --------------------------------------------------------------------------------
     
     ! modules used
+    ! --------------------------------------------------------------------------------
     use obl_math
     use obl_grid
     use obl_bc
     use obl_turb_common
 
     ! directives list
+    ! --------------------------------------------------------------------------------
     implicit none
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: advance_turbulence_closure
+    public :: define_stability_functions
+
+    public :: get_eddy_viscosity, get_eddy_diffusivity
+    public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec
+
+    public :: TKE_init, EPS_init
+    public :: TKE_bc, EPS_bc
+    ! --------------------------------------------------------------------------------
 
     !> @brief k-epsilon parameters
     type, public :: kepsilonParamType
@@ -28,6 +44,7 @@ module obl_k_epsilon
 
     contains
 
+
     ! --------------------------------------------------------------------------------
     subroutine define_stability_functions(param, bc, grid)
         !< @brief advance stability functions: N**2, S**2, Ri-gradient
@@ -65,8 +82,8 @@ module obl_k_epsilon
             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 get_eddy_viscosity(Km, TKE, EPS, param, grid)
+        call get_eddy_diffusivity(Kh, Km, param, grid)
         call TKE_bc(TKE, &
             bc%U_dyn0, bc%U_dynH, param, grid%cz)
         call eps_bc(EPS, &
@@ -91,7 +108,7 @@ module obl_k_epsilon
         end do     
     end subroutine
 
-    subroutine eps_init (eps, param, nz, z)
+    subroutine EPS_init(eps, param, nz, z)
         !< @brief TKE dissipation rate initialization
         type(kepsilonParamType), intent(in) :: param
         integer, intent(in) :: nz  !< number of z-steps
@@ -137,36 +154,74 @@ module obl_k_epsilon
 
     end subroutine EPS_bc
 
-    subroutine get_Km_k_eps (Km, TKE, eps, param, nz)
-        !< @brief calculate momentum transfer coefficient in k-epsilon closure
+    ! --------------------------------------------------------------------------------
+    subroutine get_eddy_viscosity(Km, TKE, EPS, param, grid)
+        !< @brief calculate eddy viscosity
+        ! ----------------------------------------------------------------------------
+        type(kepsilonParamType), intent(in) :: param
+        type (gridDataType), intent(in) :: grid
+        
+        real, dimension(grid%cz), intent(out) :: Km          !< eddy viscosity, [m**2 / s] 
+        real, dimension(grid%cz), intent(in) :: TKE, EPS     !<  
+        
+        ! ----------------------------------------------------------------------------
+
+        call get_eddy_viscosity_vec(Km, TKE, EPS, param, grid%cz)
+    end subroutine get_eddy_viscosity 
+
+    subroutine get_eddy_viscosity_vec(Km, TKE, EPS, param, cz)
+        !< @brief calculate eddy viscosity
+        ! ----------------------------------------------------------------------------
         type(kepsilonParamType), intent(in) :: param
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(inout) :: TKE(nz), eps(nz) 
-        real, intent(out) :: Km(nz) !eddy viscosity for momentum, [m**2 / s]
-        integer :: k  !< counter
         
-        do k = 1, nz
+        integer, intent(in) :: cz                       !< grid size
+        real, dimension(cz), intent(out) :: Km          !< eddy viscosity, [m**2 / s] 
+        real, dimension(cz), intent(in) :: TKE, EPS     !<  
+        
+        integer :: k
+        ! ----------------------------------------------------------------------------
+
+        do k = 1, cz
             if (eps(k) < param%eps_min) then
                 Km(k) = param%C_mu * TKE(k) * TKE(k) / param%eps_min
             else
                 Km(k) = param%C_mu * TKE(k) * TKE(k) / eps(k)
             endif
         end do
-    end subroutine get_Km_k_eps 
+    end subroutine get_eddy_viscosity_vec 
 
-    subroutine get_Kh_k_eps (Kh, Km, param, nz)
-        !< @brief calculate scalar transfer coefficient in k-epsilon closure
+    ! --------------------------------------------------------------------------------
+    subroutine get_eddy_diffusivity(Kh, Km, param, grid)
+        !< @brief calculate eddy diffusivity
+        ! ----------------------------------------------------------------------------
         type(kepsilonParamType), intent(in) :: param
-        integer, intent(in) :: nz !< number of z-steps
-        real, intent(in) :: Km(nz) !eddy viscosity for momentum, [m**2 / s]
-        real, intent(out) :: Kh(nz) !eddy diffusivity for tracers, [m**2 / s]
-        integer :: k  !< counter
+        type (gridDataType), intent(in) :: grid
          
-        do k = 1, nz
+        real, dimension(grid%cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
+        real, dimension(grid%cz), intent(in) :: Km       !< eddy viscosity, [m**2 / s]
+        ! ----------------------------------------------------------------------------
+
+        call get_eddy_diffusivity_vec(Kh, Km, param, grid%cz)
+    end subroutine get_eddy_diffusivity
+
+    subroutine get_eddy_diffusivity_vec(Kh, Km, param, cz)
+        !< @brief calculate eddy diffusivity
+        ! ----------------------------------------------------------------------------
+        type(kepsilonParamType), intent(in) :: param
+        
+        integer, intent(in) :: cz                   !< grid size 
+        real, dimension(cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
+        real, dimension(cz), intent(in) :: Km       !< eddy viscosity, [m**2 / s]
+        
+        integer :: k
+        ! ----------------------------------------------------------------------------
+
+        do k = 1, cz
             Kh(k) = Km(k) / param%PrT0
         end do
-    end subroutine get_Kh_k_eps 
+    end subroutine get_eddy_diffusivity_vec
 
+    ! --------------------------------------------------------------------------------
     subroutine solve_TKE_eq (Field, Km, nz, dz, dt, P_TKE, B_TKE, eps, param)
         !< @brief solver for TKE equation (explicit)
 

obl_pph_dyn.f90

@@ -92,8 +92,8 @@ module obl_pph_dyn
         type(pphDynParamType), intent(in) :: param
         type (gridDataType), intent(in) :: grid
 
-        real, intent(out) :: Km(grid%cz)     !< eddy viscosity, [m**2 / s]
-        real, intent(in) :: Ri_grad(grid%cz) !< Richardson gradient number, [n/d] 
+        real, dimension(grid%cz), intent(out) :: Km     !< eddy viscosity, [m**2 / s]
+        real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number, [n/d] 
 
         real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
         real, intent(in) :: mld         !< mixed layer depth, [m]
@@ -108,8 +108,8 @@ module obl_pph_dyn
         type(pphDynParamType), intent(in) :: param
 
         integer, intent(in) :: cz                   !< grid size 
-        real, intent(out) :: Km(cz)     !< eddy viscosity, [m**2 / s]
-        real, intent(in) :: Ri_grad(cz) !< Richardson gradient number, [n/d] 
+        real, dimension(cz), intent(out) :: Km      !< eddy viscosity, [m**2 / s]
+        real, dimension(cz), intent(in) :: Ri_grad  !< Richardson gradient number, [n/d] 
 
         real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
         real, intent(in) :: mld         !< mixed layer depth, [m]
@@ -132,8 +132,8 @@ module obl_pph_dyn
         type(pphDynParamType), intent(in) :: param
         type (gridDataType), intent(in) :: grid
 
-        real, intent(out) :: Kh(grid%cz)     !< eddy diffusivity, [m**2 / s]
-        real, intent(in) :: Ri_grad(grid%cz) !< Richardson gradient number, [n/d]
+        real, dimension(grid%cz), intent(out) :: Kh     !< eddy diffusivity, [m**2 / s]
+        real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number, [n/d]
 
         real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
         real, intent(in) :: mld         !< mixed layer depth, [m]
@@ -148,8 +148,8 @@ module obl_pph_dyn
         type(pphDynParamType), intent(in) :: param
 
         integer, intent(in) :: cz                   !< grid size 
-        real, intent(out) :: Kh(cz)     !< eddy diffusivity, [m**2 / s]
-        real, intent(in) :: Ri_grad(cz) !< Richardson gradient number, [n/d]
+        real, dimension(cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
+        real, dimension(cz), intent(in) :: Ri_grad  !< Richardson gradient number, [n/d]
 
         real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
         real, intent(in) :: mld         !< mixed layer depth, [m]