adding state vector module

CMakeLists.txt

@@ -18,6 +18,7 @@ set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules/)
 set(SOURCES
     obl_grid.f90
     obl_tslice.f90
+    obl_state.f90
     obl_time_and_space.f90
     obl_common_phys_parameters.f90
     obl_math.f90

obl_main.f90

@@ -13,6 +13,7 @@ program obl_main
 #ifdef OBL_USE_TSLICE_OUTPUT
     use obl_tslice
 #endif
+    use obl_state
     use obl_time_and_space
     use obl_initial_conditions
     use obl_forcing_and_boundary
@@ -61,14 +62,7 @@ program obl_main
     integer :: i, k  !< counters          
     integer :: status, num !< for file input/output
     real :: time_begin, time_end, time_current  !< time for output
-    real, allocatable :: Theta(:)  !< temperature, [K]
-    real, allocatable :: Theta_dev(:)  !< temperature, [K]
-    real, allocatable :: Salin(:)  !< salinity, [PSU]
-    real, allocatable :: Salin_dev(:)  !< temperature, [K]
-    real, allocatable :: U(:) !< U-velocity [m/s]	   
-    real, allocatable :: V(:) !< V-velocity [m/s]	  
-    real, allocatable :: Km(:) !< eddy viscosity for momentum, [m**2 / s]  
-    real, allocatable :: Kh(:) !eddy diffusivity for tracers, [m**2 / s]   	
+   	
     real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:) !< heat fluxes, [K*m/s] 
     real, allocatable :: Times_flux_Heat_surf(:), Times_flux_Heat_bot(:) 
     real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ???
@@ -83,12 +77,9 @@ program obl_main
     integer :: closure_mode  !< closure type:
                              !1 - pacanowski-philander, 2 - pacanowski-philander+, 
                              !3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
-    real, allocatable :: Rho(:) !< density, [kg / m**3]	 
-    real, allocatable :: Rho_dev(:) !< density, [kg / m**3]	 
+
     real :: Rho_dyn_surf, Rho_dyn_bot !< density dynamic, [kg / m**3]	 
-    real, allocatable :: N2(:) !< square of buoyancy frequency, [s**(-2)] 
-    real, allocatable :: S2(:) !< square of shear frequency, [s**(-2)]
-    real, allocatable :: Ri_grad(:)	!< Richardson gradient number 
+
     real :: mld !< mixed layer depth, [m]
     real :: lab_mld !< theoretical mixed layer depth, [m]
     real, allocatable :: TKE(:) !< TKE, [J/kg]
@@ -228,12 +219,8 @@ program obl_main
     call print_grid(grid)
 
     ! allocation
-    allocate(Theta(1:nz))
-    allocate(Theta_dev(1:nz))
-    allocate(Salin(1:nz))
-    allocate(Salin_dev(1:nz))
-    allocate(U(1:nz))
-    allocate(V(1:nz))
+    call allocate_state_vec(nz)
+
     allocate(flux_Heat_surf(1:nf))
     allocate(Times_flux_Heat_surf(1:nf))
     allocate(flux_Heat_bot(1:nf))
@@ -258,13 +245,7 @@ program obl_main
     allocate(f_Sal_right(1:nz))
     allocate(f_u_right(1:nz))
     allocate(f_v_right(1:nz))
-    allocate(Km(1:nz))
-    allocate(Kh(1:nz))
-    allocate(Rho(1:nz))
-    allocate(Rho_dev(1:nz))
-    allocate(N2(1:nz))
-    allocate(S2(1:nz))  
-    allocate(Ri_grad(1:nz)) 
+    
     allocate(TKE(1:nz)) 
     allocate(eps(1:nz)) 
     allocate(P_TKE(1:nz)) 
@@ -329,16 +310,14 @@ program obl_main
         call eps_init(eps, nz, z)
     endif
 
-    !do k = 1, nz
-    !   Km(k) =  2.0 * 0.01 !100.0
-    !  Kh(k) =  2.0 * 0.01 !100.0
-    !enddo
-
     i=1
 
     time_hrs(1) = time_current / 3600.0
     do while (time_current < time_end )
+        ! ----------------------------------------------------------------------------
         if (closure_mode.ne.5) then
+
+            !< define fluxes & dynamic scales
             call get_value_interpolate(flux_heat_surf_res, time_current, flux_heat_surf, Times_flux_heat_surf, nf) 
             call get_value_interpolate(flux_heat_bot_res, time_current, flux_heat_bot, Times_flux_heat_bot, nf)
             call get_value_interpolate(flux_sal_surf_res, time_current, flux_sal_surf, Times_flux_sal_surf, nf) 
@@ -347,8 +326,12 @@ program obl_main
             call get_value_interpolate(flux_v_surf_res, time_current, flux_v_surf, Times_flux_v_surf, nf) 
             call get_value_interpolate(flux_u_bot_res, time_current, flux_u_bot, Times_flux_u_surf, nf) 
             call get_value_interpolate(flux_v_bot_res, time_current, flux_v_bot, Times_flux_v_bot, nf)
+
+            call get_dyn_velocity(u_dynH, flux_u_surf_res, flux_v_surf_res)
             call get_rho_dyn(Rho_dyn_surf, flux_u_surf_res, flux_v_surf_res, flux_heat_surf_res, flux_sal_surf_res) 
             call get_rho_dyn(Rho_dyn_bot, flux_u_bot_res, flux_v_bot_res, flux_heat_bot_res, flux_sal_bot_res) 
+            ! ----------------------------------------------------------------------------
+            
             call solve_state_eq  (Rho, Theta_dev, Salin_dev, nz)
 #ifndef OBL_USE_GRID_DATATYPE
             call get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz) 
@@ -361,9 +344,7 @@ program obl_main
             call get_S2_on_grid(S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, grid)
 #endif
             call get_Rig(Ri_grad, N2, S2, nz)
-            call get_mld (mld, N2, dz, nz, z)
-            call get_dyn_velocity(u_dynH, flux_u_surf_res, flux_v_surf_res)
-            call get_lab_mld(lab_mld, u_dynH, N2_0, time_current, z)
+
             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)
@@ -404,6 +385,9 @@ program obl_main
         !    call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
         !    call solve_vector_eq(U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
         !    call solve_vector_eq(V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
+
+            call get_mld(mld, N2, dz, nz, z)
+            call get_lab_mld(lab_mld, u_dynH, N2_0, time_current, z)
         endif
 
         !do k = 1, nz
@@ -516,13 +500,9 @@ program obl_main
     !close(20)
 
     ! deallocation
+    call deallocate_state_vec
+
     deallocate(dz)
-    deallocate(Theta)
-    deallocate(Theta_dev)
-    deallocate(Salin)
-    deallocate(Salin_dev)
-    deallocate(U)
-    deallocate(V)
     deallocate(flux_Heat_surf)
     deallocate(Times_flux_Heat_surf)
     deallocate(flux_Heat_bot)
@@ -547,13 +527,7 @@ program obl_main
     deallocate(f_Sal_right)
     deallocate(f_u_right)
     deallocate(f_v_right)
-    deallocate(Km)
-    deallocate(Kh)
-    deallocate(Rho)
-    deallocate(Rho_dev)
-    deallocate(N2)
-    deallocate(S2)  
-    deallocate(Ri_grad) 
+    
     deallocate(TKE) 
     deallocate(eps) 
     deallocate(P_TKE) 

obl_state.f90

+module obl_state
+    !< @brief obl state def.
+    ! --------------------------------------------------------------------------------
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+
+    ! directives list
+    implicit none
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: allocate_state_vec, deallocate_state_vec
+    ! --------------------------------------------------------------------------------
+
+    !> @brief state data
+    real, allocatable, public :: Theta(:)           !< temperature, [K]
+    real, allocatable, public :: Theta_dev(:)       !< temperature dev, [K]
+
+    real, allocatable, public :: Salin(:)           !< salinity, [PSU]
+    real, allocatable, public :: Salin_dev(:)       !< salinity dev, [K]
+
+    real, allocatable, public :: Rho(:)             !< density, [kg / m**3]	 
+    real, allocatable, public :: Rho_dev(:)         !< density dev, [kg / m**3]	 
+
+    real, allocatable, public :: U(:)               !< U- current velocity, [m/s]	   
+    real, allocatable, public :: V(:)               !< V- current velocity, [m/s]
+
+    real, allocatable, public :: Km(:)              !< eddy viscosity, [m**2 / s]  
+    real, allocatable, public :: Kh(:)              !< eddy diffusivity, [m**2 / s]
+
+    real, allocatable, public :: N2(:)              !< square of buoyancy frequency, [s**(-2)] 
+    real, allocatable, public :: S2(:)              !< square of shear frequency, [s**(-2)]
+    real, allocatable, public :: Ri_grad(:)	        !< Richardson gradient number, [n/d] 
+    
+    contains
+
+    ! --------------------------------------------------------------------------------
+    subroutine allocate_state_vec(cz)
+        !> @brief allocate state vector
+        ! ----------------------------------------------------------------------------
+        integer, intent(in) :: cz
+        ! ----------------------------------------------------------------------------
+
+        allocate(Theta(cz), Theta_dev(cz))
+        allocate(Salin(cz), Salin_dev(cz))
+        allocate(Rho(cz), Rho_dev(cz))
+
+        allocate(U(cz), V(cz))
+
+        allocate(Km(cz), Kh(cz))
+
+        allocate(N2(cz), S2(cz), Ri_grad(cz))
+
+    end subroutine allocate_state_vec
+
+    ! --------------------------------------------------------------------------------
+    subroutine deallocate_state_vec
+        !> @brief deallocate state vector
+        ! ----------------------------------------------------------------------------
+        ! ----------------------------------------------------------------------------
+
+        deallocate(Theta, Theta_dev)
+        deallocate(Salin, Salin_dev)
+        deallocate(Rho, Rho_dev)
+
+        deallocate(U, V)
+        
+        deallocate(Km, Kh)
+
+        deallocate(N2, S2, Ri_grad)
+
+    end subroutine deallocate_state_vec
+
+end module 
\ No newline at end of file