main time advancement update

obl_main.f90

@@ -243,6 +243,8 @@ program obl_main
     call Theta_init (Theta, Theta_dev, grid%cz, grid%dz)
     !call Theta_init (Theta, nz, dz)
     call Salin_init(Salin, Salin_dev, grid%cz, grid%dz)
+    call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
+
     call U_init(U, grid%cz)
     call V_init(V, grid%cz)
 
@@ -365,33 +367,25 @@ program obl_main
             bc%u_momentum_flux0, bc%v_momentum_flux0, bc%heat_flux0, bc%salin_flux0) 
         ! ----------------------------------------------------------------------------
 
-        call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
-#ifndef OBL_USE_GRID_DATATYPE
-        call get_N2(N2, Rho, bc%Rho_dynH, bc%rho_dyn0, grid%dz, grid%cz) 
-#else
+        !< advance turbulence closure
+        ! ----------------------------------------------------------------------------
+        !< N^2, S^2 & Ri(g)
         call get_N2_on_grid(N2, Rho, bc%Rho_dynH, bc%Rho_dyn0, grid)
-#endif
-#ifndef OBL_USE_GRID_DATATYPE
-        call get_S2(S2, U, V, &
-            bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%v_momentum_flux0, grid%dz, grid%cz)
-#else
         call get_S2_on_grid(S2, U, V, &
             bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%u_momentum_flux0, bc%u_momentum_flux0, grid)
-#endif
         call get_Rig(Ri_grad, N2, S2, grid%cz)
 
+        !< TKE production & buoyancy 
+        call get_TKE_generation(TKE_production, Km, S2, grid%cz)
+        call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
+        
         if (closure_mode.eq.1) then
             call get_eddy_viscosity(Km, Ri_grad, param_pacanowski, grid)
             call get_eddy_diffusivity(Kh, Ri_grad, param_pacanowski, grid)
-
-            call get_TKE_generation(TKE_production, Km, S2, grid%cz)
-            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
         else if (closure_mode.eq.2) then
             call get_Km_plus(Km, Ri_grad, &
                 bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, grid%cz, grid%dz, grid%height)
             call get_Kh_plus(Kh, Km, grid%cz)
-            call get_TKE_generation(TKE_production, Km, S2, grid%cz)
-            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
         else if (closure_mode.eq.4) then
             call TKE_bc(TKE, &
                 bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz)
@@ -399,8 +393,6 @@ program obl_main
                 bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz, grid%dz)
             call get_Km_k_eps (Km, TKE, EPS, grid%cz)
             call get_Kh_k_eps (Kh, Km, grid%cz)
-            call get_TKE_generation(TKE_production, Km, S2, grid%cz)
-            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid%cz)
             call TKE_bc(TKE, &
                 bc%u_momentum_fluxH, bc%v_momentum_fluxH, grid%cz)
             call eps_bc(EPS, &
@@ -413,15 +405,20 @@ program obl_main
             call solve_eps_eq_semiimplicit (EPS, Km_eps, grid%cz, grid%dz, dt, TKE_production, TKE_buoyancy, TKE)
             call limit_min_array(EPS, eps_min, grid%cz) 
         endif
+        ! ----------------------------------------------------------------------------
         
+        !< advance dynamics
+        ! ----------------------------------------------------------------------------
         call solve_scalar_eq (Theta_dev, Kh, grid%cz, grid%dz, dt, &
             bc%heat_fluxH, bc%heat_flux0, f_heat_right) 
         call solve_scalar_eq (Salin_dev, Kh, grid%cz, grid%dz, dt, &
             bc%salin_fluxH, bc%salin_flux0, f_sal_right)
+        call solve_state_eq(Rho, Theta_dev, Salin_dev, grid%cz)
         call solve_vector_eq (U, Km, grid%cz, grid%dz, dt, &
             bc%u_momentum_fluxH, bc%u_momentum_flux0, f_u_right)
         call solve_vector_eq (V, Km, grid%cz, grid%dz, dt, &
             bc%v_momentum_fluxH, bc%v_momentum_flux0, f_v_right)
+        ! ----------------------------------------------------------------------------
 
         !> advance time
         i = i + 1

obl_turb_closure.f90

@@ -17,6 +17,7 @@ module obl_turb_closure
     ! public interface
     ! --------------------------------------------------------------------------------
     public :: get_N2, get_S2, get_Rig
+    public :: get_N2_on_grid, get_S2_on_grid
     public :: get_TKE_generation, get_TKE_buoyancy
     ! --------------------------------------------------------------------------------
     
@@ -64,7 +65,6 @@ module obl_turb_closure
         
     end subroutine
 
-#ifdef OBL_USE_GRID_DATATYPE
     subroutine get_N2_on_grid(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, grid)
         !< @brief calculate square of buoyancy frequency
 
@@ -99,7 +99,6 @@ module obl_turb_closure
         N2(grid%cz) = - g/Rho_ref * dRho_tmp_surf
         
     end subroutine
-#endif
 
     subroutine get_S2(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, dz, nz)
         !< @brief calculate square of shear frequency
@@ -140,7 +139,6 @@ module obl_turb_closure
 
     end subroutine
 
-#ifdef OBL_USE_GRID_DATATYPE
     subroutine get_S2_on_grid(S2, U, V, flux_u_surf, flux_u_bot, flux_v_surf, flux_v_bot, grid)
         !< @brief calculate square of shear frequency
 
@@ -179,7 +177,6 @@ module obl_turb_closure
         call limit_min_array(S2, S2_min, grid%cz) 
 
     end subroutine
-#endif
 
     subroutine get_Rig(Ri_grad, N2, S2, nz)
         !< @brief calculate Richardson gradient number