adding turbulence closure explicit initialization

obl_k_epsilon.f90

@@ -16,6 +16,7 @@ module obl_k_epsilon
 
     ! public interface
     ! --------------------------------------------------------------------------------
+    public :: init_turbulence_closure
     public :: advance_turbulence_closure
     public :: define_stability_functions
 
@@ -63,6 +64,30 @@ module obl_k_epsilon
         call get_Ri_gradient(Ri_grad, N2, S2, param%c_S2_min, grid)
     end subroutine
 
+    ! --------------------------------------------------------------------------------
+    subroutine init_turbulence_closure(param, bc, grid)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(kepsilonParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+        ! --------------------------------------------------------------------------------
+
+        call TKE_init(TKE, param, grid%cz)
+        call EPS_init(EPS, param, grid%cz, grid%height)
+
+        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_eddy_viscosity(Km, TKE, EPS, param, grid)
+        call get_eddy_diffusivity(Kh, Km, param, grid)
+
+    end subroutine
+
     ! --------------------------------------------------------------------------------
     subroutine advance_turbulence_closure(param, bc, grid, dt)
         !< @brief advance turbulence closure

obl_main.f90

@@ -27,12 +27,15 @@ program obl_main
     use obl_bc
     use obl_pph, only: pphParamType, &
         define_stability_functions_pph => define_stability_functions, &
+        init_turbulence_closure_pph => init_turbulence_closure, & 
         advance_turbulence_closure_pph => advance_turbulence_closure
     use obl_pph_dyn, only: pphDynParamType, &
         define_stability_functions_pph_dyn => define_stability_functions, &
+        init_turbulence_closure_pph_dyn => init_turbulence_closure, & 
         advance_turbulence_closure_pph_dyn => advance_turbulence_closure
     use obl_k_epsilon, only: kepsilonParamType, &
         define_stability_functions_k_epsilon => define_stability_functions, &
+        init_turbulence_closure_k_epsilon => init_turbulence_closure, & 
         advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
         TKE_init, EPS_init, TKE_bc, EPS_bc
     use obl_io_plt
@@ -191,8 +194,6 @@ program obl_main
         write(*, *) ' FAILURE! > unable to set time '
         return
     endif
-    time_current = time_begin
-    time_index = 1
     ! ----------------------------------------------------------------------------  
 
     !< allocating state vector
@@ -238,16 +239,24 @@ program obl_main
 
     !< initialization of turbulence closure
     ! ---------------------------------------------------------------------------- 
-    if (closure_mode.eq.3 .or. closure_mode.eq.4) then
-        call TKE_init(TKE, param_k_epsilon, grid%cz)
-        call EPS_init(EPS, param_k_epsilon, grid%cz, grid%height)
+    !< define fluxes & dynamic scales at startup
+    call advance_surface_fluxes(bc, time_begin, grid)
+    call advance_bottom_fluxes(bc, time_begin, grid)
 
-        !< have to define Km, Kh at init
-        Km = 0.0
-        Kh = 0.0
+    if (closure_mode.eq.1) then
+        call define_stability_functions_pph(param_pph, bc, grid)
+        call init_turbulence_closure_pph(param_pph, bc, grid)
+    else if (closure_mode.eq.2) then
+        call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
+        call init_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid)
+    else if (closure_mode.eq.4) then
+        call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
+        call init_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid)
     endif
     ! ---------------------------------------------------------------------------- 
 
+    time_current = time_begin
+    time_index = 1
     do while (time_current < time_end )
         ! ----------------------------------------------------------------------------
         

obl_pph.f90

@@ -15,6 +15,7 @@ module obl_pph
 
     ! public interface
     ! --------------------------------------------------------------------------------
+    public :: init_turbulence_closure
     public :: advance_turbulence_closure
     public :: define_stability_functions
 
@@ -58,6 +59,21 @@ module obl_pph
         call get_Ri_gradient(Ri_grad, N2, S2, param%c_S2_min, grid)
     end subroutine
 
+    ! --------------------------------------------------------------------------------
+    subroutine init_turbulence_closure(param, bc, grid)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(pphParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+        ! --------------------------------------------------------------------------------
+
+        call get_eddy_viscosity(Km, Ri_grad, param, grid)
+        call get_eddy_diffusivity(Kh, Ri_grad, param, grid)
+    end subroutine
+
     ! --------------------------------------------------------------------------------
     subroutine advance_turbulence_closure(param, bc, grid, dt)
         !< @brief advance turbulence closure

obl_pph_dyn.f90

@@ -17,6 +17,7 @@ module obl_pph_dyn
 
     ! public interface
     ! --------------------------------------------------------------------------------
+    public :: init_turbulence_closure
     public :: advance_turbulence_closure
     public :: define_stability_functions
 
@@ -60,6 +61,27 @@ module obl_pph_dyn
         call get_Ri_gradient(Ri_grad, N2, S2, param%c_S2_min, grid)
     end subroutine
 
+    ! --------------------------------------------------------------------------------
+    subroutine init_turbulence_closure(param, bc, grid)
+        !< @brief advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        use obl_state
+
+        type(pphDynParamType), intent(in) :: param
+        type(oblBcType), intent(in) :: bc
+        type (gridDataType), intent(in) :: grid
+
+        real :: mld
+        ! --------------------------------------------------------------------------------
+
+        call get_mld(mld, N2, grid%dz, grid%cz)
+
+        call get_eddy_viscosity(Km, Ri_grad, &
+            bc%U_dynH, mld, param, grid)
+        call get_eddy_diffusivity(Kh, Ri_grad, &
+            bc%U_dynH, mld, param, grid)
+    end subroutine
+
     ! --------------------------------------------------------------------------------
     subroutine advance_turbulence_closure(param, bc, grid, dt)
         !< @brief advance turbulence closure