diff --git a/obl_main.f90 b/obl_main.f90
index 35fa113e19292db991921b6e3573a017c5ef3e1c..b88073dad87b02df460dca2680d541fc606444eb 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -1,7 +1,8 @@
 #include "obl_def.fi"
 
 program obl_main
-    !< @brief main program for calculations for ocean boundary layer
+    !< @brief ocean boundary layer (OBL) main
+    ! --------------------------------------------------------------------------------
 
 #ifdef USE_CONFIG_PARSER
     use iso_c_binding, only: C_NULL_CHAR
@@ -9,6 +10,7 @@ program obl_main
 #endif
 
     ! modules used
+    ! --------------------------------------------------------------------------------
     use obl_grid
     use obl_state
     use obl_state_eq
@@ -20,6 +22,7 @@ program obl_main
     use obl_init
     use obl_fluxes
     use obl_scm
+    use obl_config
     use obl_diag
     use obl_bc
     use obl_pph, only: pphParamType, &
@@ -36,17 +39,27 @@ program obl_main
     use io
     use io_metadata
  
-    use obl_config
-
-    !use vertical_mixing, default = off
-    !use vermix
-
     ! directives list
+    ! --------------------------------------------------------------------------------
     implicit none
 
+    ! model data
+    ! --------------------------------------------------------------------------------
+    integer :: obl_setup    ! --- OBL builtin setup 
+                            ! = setup_kato: Kato-Phillips
+                            ! = setup_papa_fluxes: Papa-station (fluxes)
+                            ! = setup_papa_meteo: Papa-station (meteo)
+
+    integer :: closure_mode     ! --- OBL closure def.
+                                ! = 1 - pacanowski-philander
+                                ! = 2 - pacanowski-philander dynamic 
+                                ! = 3 - k-epsilon (explicit)
+                                ! = 4 - k-epsilon (semiimplicit)
+                                ! = 5 - inmom
+
     type(gridDataType) :: grid
 
-    ! turbulence closure parameters
+    !< turbulence closure parameters
     type(pphParamType) :: param_pph
     type(pphDynParamType) :: param_pph_dyn
     type(kepsilonParamType) :: param_k_epsilon
@@ -54,37 +67,32 @@ program obl_main
     !< boundary conditions data
     type(oblBcType) :: bc
 
-    real :: domain_height
-    integer :: grid_cz
-
     !< output 
     type(oblOutputStx) :: scm_output
 
+    !< time
+    real :: time_begin, time_end, time_current
+    real :: dt
+    integer :: time_index
 
-    real :: dt !< time step [s] 
-    integer :: i, k  !< counters          
-    integer :: status, num !< for file input/output
-    real :: time_begin, time_end, time_current  !< time for output
-   	 
-
-    integer :: closure_mode  !< closure type:
-                             !1 - pacanowski-philander, 2 - pacanowski-philander+, 
-                             !3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
-
-    integer, parameter :: output_mode = 3   ! 1 -- netcdf, 2 -- ascii, 3 -- tecplot
-
-    integer :: obl_setup     ! 1 - Kato-Phillips, 2 - Papa station (fluxes), 3 - Papa station (meteo)
+    !< screen output parameters
+    integer, parameter :: nscreen = 1000
 
-    real, parameter :: Cd0 = 0.001           ! bottom drag coefficient
+    !< file output parameters
+    integer :: output_mode      ! --- OBL output mode
+                                ! = 1 -- netcdf
+                                ! = 2 -- ascii
+                                ! = 3 -- tecplot
 
-	 
+    integer, parameter :: noutput = 60
 
-    real :: mld !< mixed layer depth, [m]
-    real :: lab_mld !< theoretical mixed layer depth, [m]
+    !< additional variables & parameters
+    real, parameter :: Cd0 = 0.001      ! bottom drag coefficient
 
-    real :: N2_0
+    real :: mld, mld_ref                    ! mixed layer depth (model & reference), [m]
+    real, parameter :: N2_ref = 0.00044     ! reference N**2 (using in mld ref. estimate), [1/s**2]
     
-    ! command line arguments
+    ! command line arguments & configuration file variables
     ! --------------------------------------------------------------------------------
     integer :: num_args
     character(len = 128) :: arg
@@ -93,20 +101,32 @@ program obl_main
     character(len = 128), parameter :: arg_key_setup = '--setup'
     character(len = 128), parameter :: arg_key_config = "--config"
 
+    integer :: i, status
     integer :: ierr
     ! --------------------------------------------------------------------------------
 
-    ! screen output parameters
-    integer, parameter :: nscreen = 1000
-
-    ! file output parameters
-    integer, parameter :: noutput = 60
-
-
-
-    obl_setup = setup_kato !< ! 0 - Kato-Phillips (default), 1 - Papa station (fluxes), 2 - Papa station (meteo)
-    closure_mode = 4 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
 
+    !< default setup = setup_kato (Kato-Phiilps)
+    !< possible values:
+    !<      setup_kato
+    !<      setup_papa_fluxes
+    !<      setup_papa_meteo
+    obl_setup = setup_kato
+
+    !< default closure = 4, k-epsilon (implicit)
+    !< poissible values:
+    !<      = 1, Pacanowski-Philander
+    !<      = 2, Pacanowski-Philander (dynamic)
+    !<      = 3, k-epsilon (explicit) *: DEPRECATED
+    !<      = 4, k-epsilon (implicit)
+    closure_mode = 4
+
+    !< default output = 1, (netcdf)
+    !< possible values:
+    !<      = 1, netcdf
+    !<      = 2, ascii
+    !<      = 3, tecplot
+    output_mode = 3
 
     ! --- command line arguments processing
     num_args = command_argument_count()
@@ -155,7 +175,7 @@ program obl_main
             end if
 
             !< forcing configuration file setup
-            obl_setup = -1
+            obl_setup = 999
 #endif
         endif
     enddo
@@ -178,13 +198,18 @@ program obl_main
         return
     endif 
     time_current = time_begin
+    time_index = 1
     ! ----------------------------------------------------------------------------  
 
-    ! allocation
+    !< allocating state vector
+    ! ----------------------------------------------------------------------------  
     call allocate_state_vec(grid%cz)
+    ! ----------------------------------------------------------------------------  
 
-    ! initialize scm
+    !< initialize scm
+    ! ----------------------------------------------------------------------------  
     call init_scm_vec(grid%cz)
+    ! ----------------------------------------------------------------------------  
 
     !< setting phys
     ! ----------------------------------------------------------------------------
@@ -215,21 +240,13 @@ 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)
     endif
     ! ---------------------------------------------------------------------------- 
 
-    !< setting reference data
-    ! ----------------------------------------------------------------------------
-    N2_0 = 0.00044
-    ! ----------------------------------------------------------------------------
-
-    status = 0
-    num = 0
-
-    i=1
     do while (time_current < time_end )
         ! ----------------------------------------------------------------------------
         
@@ -263,13 +280,15 @@ program obl_main
         ! ----------------------------------------------------------------------------
 
         !> advance time
-        i = i + 1
+        ! ----------------------------------------------------------------------------
+        time_index = time_index + 1
         time_current = time_current + dt
+        ! ----------------------------------------------------------------------------
 
         !> advance screen output
-        if (mod(i, nscreen) == 0) then
+        ! ----------------------------------------------------------------------------
+        if (mod(time_index, nscreen) == 0) then
             call get_mld(mld, N2, grid%dz, grid%cz)
-            call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
 
             ! *: add finite check
 
@@ -279,16 +298,18 @@ program obl_main
                 (time_current / time_end) * 100.0, '% ]' 
             write(*, '(a)') '-------------------------------------------------'
         endif
+        ! ----------------------------------------------------------------------------
 
         !> advance file output
-        if (mod(i, noutput) == 0) then
+        ! ----------------------------------------------------------------------------
+        if (mod(time_index, noutput) == 0) then
             call push_state_vec(scm_output, grid%cz)
 
             call get_mld(mld, N2, grid%dz, grid%cz)
-            call get_mld_ref(lab_mld, bc%U_dynH, N2_0, time_current, grid%height)
+            call get_mld_ref(mld_ref, bc%U_dynH, N2_ref, time_current, grid%height)
 
             call push_value_to_tseries(scm_output%mld, mld)
-            call push_value_to_tseries(scm_output%mld_ref, lab_mld)
+            call push_value_to_tseries(scm_output%mld_ref, mld_ref)
 
             call push_value_to_tseries(scm_output%tau_x, bc%u_momentum_fluxH * Rho_ref)
             call push_value_to_tseries(scm_output%tau_y, bc%v_momentum_fluxH * Rho_ref)
@@ -297,8 +318,11 @@ program obl_main
 
             call push_value_to_tseries(scm_output%time, time_current / 3600.0)
         endif
+        ! ----------------------------------------------------------------------------
     enddo
 
+    !> writing file output
+    ! ----------------------------------------------------------------------------
     if (output_mode.eq.1) then
         write(*, *) ' >> writing netcdf output ...'
         call write_netcdf(scm_output, grid%z)
@@ -314,10 +338,10 @@ program obl_main
         write(*, *) ' >> writing tecplot output ...'
         call write_tecplot(scm_output, grid%z)
     endif
+    ! ----------------------------------------------------------------------------
 
-
-
-
+    !> model cleanup
+    ! ----------------------------------------------------------------------------
     !> deallocate state
     call deallocate_state_vec
 
@@ -332,6 +356,6 @@ program obl_main
 
     ! > removing grid data
     call deallocate_grid(grid)
-
+    ! ----------------------------------------------------------------------------
 
 end program