diff --git a/obl_k_epsilon.f90 b/obl_k_epsilon.f90
index 78b2874c8ee923515a93e1bed625f3855bb5e5db..182dd0642a652381936bb84a36e1bfccfc924017 100644
--- a/obl_k_epsilon.f90
+++ b/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
diff --git a/obl_main.f90 b/obl_main.f90
index 3a340feaa2c0bd3ec7687ae4bb3660226ec2b18c..9a4eb5e3a5883a84d3832fd3049ac45aeebb1f19 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -26,13 +26,16 @@ program obl_main
     use obl_diag
     use obl_bc
     use obl_pph, only: pphParamType, &
-        define_stability_functions_pph => define_stability_functions, & 
+        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, & 
+        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, & 
+        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
@@ -190,9 +193,7 @@ program obl_main
     if (ierr /= 0) then
         write(*, *) ' FAILURE! > unable to set time '
         return
-    endif 
-    time_current = time_begin
-    time_index = 1
+    endif
     ! ----------------------------------------------------------------------------  
 
     !< allocating state vector
@@ -237,17 +238,25 @@ 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)
-
-        !< have to define Km, Kh at init
-        Km = 0.0
-        Kh = 0.0
+    ! ---------------------------------------------------------------------------- 
+    !< define fluxes & dynamic scales at startup
+    call advance_surface_fluxes(bc, time_begin, grid)
+    call advance_bottom_fluxes(bc, time_begin, grid)
+
+    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 )
         ! ----------------------------------------------------------------------------
         
diff --git a/obl_pph.f90 b/obl_pph.f90
index 53c33da6f0bbf9ad139209f74fa7f51392f109f4..fa63efb63f3d49b6915b4898bf758c6cb612fbb9 100644
--- a/obl_pph.f90
+++ b/obl_pph.f90
@@ -15,6 +15,7 @@ module obl_pph
 
     ! public interface
     ! --------------------------------------------------------------------------------
+    public :: init_turbulence_closure
     public :: advance_turbulence_closure
     public :: define_stability_functions
 
@@ -39,7 +40,7 @@ module obl_pph
 
     contains
 
-    
+
     ! --------------------------------------------------------------------------------
     subroutine define_stability_functions(param, bc, grid)
         !< @brief advance stability functions: N**2, S**2, Ri-gradient
@@ -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
diff --git a/obl_pph_dyn.f90 b/obl_pph_dyn.f90
index a94143fe0bfce1118c6e28013f8819ace4447807..ebb94268382557a4a4a54a84c35196b3662f1458 100644
--- a/obl_pph_dyn.f90
+++ b/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