diff --git a/CMakeLists.txt b/CMakeLists.txt
index e37ac3d631992abdbcc38c940c3d5f401b229ff3..b09f16278f3aa758a0a8ec869e0be36d57307126 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -47,6 +47,7 @@ set(SOURCES
     obl_k_epsilon.f90
     obl_pph.f90
     obl_pph_dyn.f90
+    obl_fluxes.f90
     obl_scm.f90
     obl_main.f90
     io_metadata.f90
diff --git a/obl_common_phys.f90 b/obl_common_phys.f90
index 3abf9434687444e2cef36c797ff7afeee977b81b..a3bb256fd151379801fa3f5a0ed240a443b1ea13 100644
--- a/obl_common_phys.f90
+++ b/obl_common_phys.f90
@@ -22,6 +22,10 @@ module obl_common_phys
 
     real, public, parameter :: R_d = 287.0              !< dry air gas constant [J / (K * kg)]
     real, public, parameter :: R_v = 461.0              !< water vapor gas constant [J / (K * kg)]
+
+    real, public, parameter :: rho_air = 1.22           !< air reference density [kg / m**3]
+
+    real, public, parameter :: L_v = 2.5 * 1e6          !< latenth heat of vaporization [J/kg]
     ! --------------------------------------------------------------------------------
 
 
diff --git a/obl_fluxes.f90 b/obl_fluxes.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b53dcd03cc6f0ad0cc8e9d34ffc8d377c5bdf78a
--- /dev/null
+++ b/obl_fluxes.f90
@@ -0,0 +1,278 @@
+module obl_fluxes
+    !< @brief obl fluxes calculation & setup
+    ! --------------------------------------------------------------------------------
+
+    ! TO DO:
+    !   -- ***
+
+    ! modules used
+    ! --------------------------------------------------------------------------------
+    use obl_tforcing
+
+#ifdef USE_SFX
+    use sfx_data, only: meteoDataType, sfxDataType
+    use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
+                  numericsType_most => numericsType
+#endif
+
+    ! directives list
+    ! --------------------------------------------------------------------------------
+    implicit none
+    private
+
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: advance_surface_fluxes
+    public :: advance_bottom_fluxes
+    public :: deallocate_fluxes_data
+    ! --------------------------------------------------------------------------------
+
+    !> @brief setup mode [bool]
+    integer, public :: is_meteo_setup
+
+    !> @brief surface forcing 
+    type(timeForcingDataType), public :: sensible_hflux_surf, latent_hflux_surf     !< heat fluxes, [W/m^2] 
+    type(timeForcingDataType), public :: salin_flux_surf                            !< salinity flux, [PSU*m/s] 
+    type(timeForcingDataType), public :: tau_x_surf, tau_y_surf                     !< momentum flux, [N/m**2]
+
+    type(timeForcingDataType), public :: Ua, Va                                  !< air wind, [m/s]
+    type(timeForcingDataType), public :: Ta                                      !< air temperature, [K]
+    type(timeForcingDataType), public :: Pa                                      !< air pressure, [Pa]
+    type(timeForcingDataType), public :: Qa                                      !< air humidity, [kg/kg]
+    type(timeForcingDataType), public :: RHa                                     !< air relative humidity, [%]
+
+    type(timeForcingDataType), public :: sw_flux_surf                    !< shortwave radiation flux IN, [W/m^2]
+    type(timeForcingDataType), public :: lw_in_surf, lw_net_surf         !< longwave radiation flux IN/NET, [W/m^2] 
+
+    !> @brief bottom forcing
+    type(timeForcingDataType), public :: hflux_bot                  !< heat flux, [W/m^2]
+    type(timeForcingDataType), public :: salin_flux_bot             !< salinity flux, [PSU*m/s] 
+    type(timeForcingDataType), public :: tau_x_bot, tau_y_bot       !< momentum flux, [N/m**2]
+
+    real, parameter :: lw_albedo = 0.03
+    real, parameter :: surface_emissivity = 0.97
+
+#ifdef USE_SFX
+    type(meteoDataType) :: meteo
+    type(sfxDataType) :: sfx
+    type(numericsType_most) :: sfx_numerics
+#endif
+    
+    contains
+
+
+    ! --------------------------------------------------------------------------------
+    subroutine advance_surface_fluxes(bc, time_current, grid)
+        !< @brief set fluxes & dynamic scales b.c. (surface)
+        ! ----------------------------------------------------------------------------
+        use obl_grid
+        use obl_state
+        use obl_bc
+        use obl_state_eq
+        use obl_common_phys
+
+        type(oblBcType), intent(out) :: bc
+        real, intent(in) :: time_current
+        type (gridDataType), intent(in) :: grid
+
+        real :: fvalue
+
+        real :: lw_net_t
+        ! --------------------------------------------------------------------------------
+
+        !< define LW net radiation 
+        lw_net_t = 0.0
+        if (lw_net_surf%num > 0) then
+            call get_value_tforcing(lw_net_t, time_current, lw_net_surf)
+        else if (lw_in_surf%num > 0) then
+            call get_value_tforcing(fvalue, time_current, lw_in_surf)
+            lw_net_t = (1.0 - lw_albedo) * fvalue
+
+            !< adding LW outgoing flux to balance
+            block
+                real :: Theta_surf 
+                real :: lw_out_surf
+            
+                Theta_surf = Theta_dev(grid%cz) + Theta_ref
+                lw_out_surf = surface_emissivity * stefan_boltzmann_const * & 
+                    Theta_surf * Theta_surf * Theta_surf * Theta_surf
+
+                lw_net_t = lw_net_t - lw_out_surf
+            end block  
+        endif
+
+        if (is_meteo_setup == 0) then
+            !< using 'flux' type forcing
+            ! --------------------------------------------------------------------------------
+
+            !< heat flux  
+            bc%heat_fluxH = 0.0
+            call get_value_tforcing(fvalue, time_current, sensible_hflux_surf)
+            bc%heat_fluxH = bc%heat_fluxH + fvalue
+            call get_value_tforcing(fvalue, time_current, latent_hflux_surf)
+            bc%heat_fluxH = bc%heat_fluxH + fvalue
+
+            bc%heat_fluxH = bc%heat_fluxH + lw_net_t
+
+            !< momentum flux 
+            call get_value_tforcing(bc%u_momentum_fluxH, time_current, tau_x_surf)
+            call get_value_tforcing(bc%v_momentum_fluxH, time_current, tau_y_surf)
+
+        else if (is_meteo_setup == 1) then
+#ifdef USE_SFX
+            !< using 'meteo' type forcing
+            ! --------------------------------------------------------------------------------
+
+            block
+                real :: Ua_t, Va_t, Ta_t, Pa_t, Qa_t, RHa_t, e_sat_a
+                real :: Theta_surf, Q_surf, e_sat
+
+                call get_value_tforcing(Ua_t, time_current, Ua)
+                call get_value_tforcing(Va_t, time_current, Va)
+                call get_value_tforcing(Ta_t, time_current, Ta)
+                call get_value_tforcing(Pa_t, time_current, Pa)
+
+                !< define air specific humidity
+                if (Qa%num > 0) then 
+                    call get_value_tforcing(Qa_t, time_current, Qa)
+                else if (RHa%num > 0) then 
+                    call get_value_tforcing(RHa_t, time_current, RHa)
+
+                    call get_water_saturation_vapor_pressure_in_kPa(e_sat_a, Ta_t)
+                    call get_specific_humidity(Qa_t, e_sat_a, Pa_t / 1000.0, R_d, R_v)
+
+                    Qa_t = Qa_t * (RHa_t / 100.0)
+                endif
+
+                Theta_surf = Theta_dev(grid%cz) + Theta_ref
+
+                !< define saturation specific humidity
+                call get_water_saturation_vapor_pressure_in_kPa(e_sat, Theta_surf)
+                call get_specific_humidity(Q_surf, e_sat, Pa_t / 1000.0, R_d, R_v)
+
+                meteo%U = sqrt(Ua_t * Ua_t + Va_t * Va_t)
+                meteo%dT = Ta_t - Theta_surf
+                meteo%Tsemi = 0.5 * (Ta_t + Theta_surf)
+                meteo%dQ = Qa_t - Q_surf
+                meteo%h = 2.0
+                meteo%z0_m = -1.0
+
+                call get_surface_fluxes_most(sfx, meteo, sfx_numerics)
+
+                !< heat flux  
+                bc%heat_fluxH = &
+                    rho_air * cp_air * sfx%Cm * sfx%Ct * meteo%U * meteo%dT + &
+                    rho_air * L_v * sfx%Cm * sfx%Ct * meteo%U * meteo%dQ 
+
+                bc%heat_fluxH = bc%heat_fluxH + lw_net_t
+
+                !< momentum flux 
+                bc%u_momentum_fluxH = rho_air * sfx%Cm * sfx%Cm * meteo%U * Ua_t
+                bc%v_momentum_fluxH = rho_air * sfx%Cm * sfx%Cm * meteo%U * Va_t
+
+            end block 
+#else
+            write(*, *) ' >> FAILURE! Meteo setup requires SFX'
+            stop
+#endif
+        endif
+
+        !< kinematic heat flux = F / (rho_ref * cp)
+        bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water) 
+
+        !< kinematic momentum flux = tau / rho_ref
+        bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
+        bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
+
+        !< salinity flux 
+        call get_value_tforcing(bc%salin_fluxH, time_current, salin_flux_surf)
+
+        !< shortwave radiation flux [W/m^2]
+        call get_value_tforcing(bc%sw_fluxH, time_current, sw_flux_surf)
+
+        !< U* def.:
+        call get_u_dynamic(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)
+
+        !< rho* def.:
+        call get_rho_dynamic(bc%rho_dynH, &
+            bc%U_dynH, bc%heat_fluxH, bc%salin_fluxH) 
+
+    end subroutine
+
+    ! --------------------------------------------------------------------------------
+    subroutine advance_bottom_fluxes(bc, time_current, grid)
+        !< @brief set fluxes & dynamic scales b.c. (bottom)
+        ! ----------------------------------------------------------------------------
+        use obl_grid
+        use obl_state
+        use obl_bc
+        use obl_state_eq
+        use obl_common_phys
+
+        type(oblBcType), intent(out) :: bc
+        real, intent(in) :: time_current
+        type (gridDataType), intent(in) :: grid
+        ! --------------------------------------------------------------------------------
+
+        !< heat flux  
+        call get_value_tforcing(bc%heat_flux0, time_current, hflux_bot)
+
+        !< kinematic heat flux = F / (rho_ref * cp)
+        bc%heat_flux0 = bc%heat_flux0 / (Rho_ref * cp_water)
+
+        !< salinity flux 
+        call get_value_tforcing(bc%salin_flux0, time_current, salin_flux_bot)
+
+        !< momentum flux 
+        call get_value_tforcing(bc%u_momentum_flux0, time_current, tau_x_bot)
+        call get_value_tforcing(bc%v_momentum_flux0, time_current, tau_y_bot)
+
+        !< kinematic momentum flux = tau / rho_ref
+        bc%u_momentum_flux0 = bc%u_momentum_flux0 / Rho_ref
+        bc%v_momentum_flux0 = bc%v_momentum_flux0 / Rho_ref
+
+        !< U* def.:
+        call get_u_dynamic(bc%U_dyn0, bc%u_momentum_flux0, bc%v_momentum_flux0)
+
+        !< rho* def.:
+        call get_rho_dynamic(bc%rho_dyn0, &
+            bc%U_dyn0, bc%heat_flux0, bc%salin_flux0) 
+
+    end subroutine
+
+    ! --------------------------------------------------------------------------------
+    subroutine deallocate_fluxes_data
+        !> @brief deallocate fluxes data
+        ! ----------------------------------------------------------------------------
+        ! ----------------------------------------------------------------------------
+
+        call deallocate_tforcing(sensible_hflux_surf)
+        call deallocate_tforcing(latent_hflux_surf)
+    
+        call deallocate_tforcing(salin_flux_surf)
+    
+        call deallocate_tforcing(tau_x_surf)
+        call deallocate_tforcing(tau_y_surf)
+    
+        call deallocate_tforcing(hflux_bot)
+        call deallocate_tforcing(salin_flux_bot)
+    
+        call deallocate_tforcing(tau_x_bot)
+        call deallocate_tforcing(tau_y_bot)
+    
+        call deallocate_tforcing(Ua)
+        call deallocate_tforcing(Va)
+    
+        call deallocate_tforcing(Ta)
+        call deallocate_tforcing(Pa)
+        call deallocate_tforcing(Qa)
+        call deallocate_tforcing(RHa)
+    
+        call deallocate_tforcing(sw_flux_surf)
+    
+        call deallocate_tforcing(lw_in_surf)
+        call deallocate_tforcing(lw_net_surf)
+
+    end subroutine deallocate_fluxes_data
+
+end module 
diff --git a/obl_main.f90 b/obl_main.f90
index 8c5f535a562d1b94bf4934050dbdb612977aa5a3..04209fc13a117da2b6eb5e8d8be9710f381a1b8f 100644
--- a/obl_main.f90
+++ b/obl_main.f90
@@ -18,6 +18,7 @@ program obl_main
     use obl_tseries
     use obl_output
     use obl_init
+    use obl_fluxes
     use obl_scm
     use obl_diag
     use obl_bc
@@ -34,11 +35,6 @@ program obl_main
     use obl_io_plt
     use io
     use io_metadata
-#ifdef USE_SFX
-    use sfx_data, only: meteoDataType, sfxDataType
-    use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
-                  numericsType_most => numericsType
-#endif
 #ifdef USE_CONFIG_PARSER
    use config_parser
 #endif
@@ -48,12 +44,6 @@ program obl_main
     ! directives list
     implicit none
 
-#ifdef USE_SFX
-    type(meteoDataType) :: meteo
-    type(sfxDataType) :: sfx
-    type(numericsType_most) :: sfx_numerics
-#endif
-
     type(gridDataType) :: grid
 
     ! turbulence closure parameters
@@ -61,27 +51,6 @@ program obl_main
     type(pphDynParamType) :: param_pph_dyn
     type(kepsilonParamType) :: param_k_epsilon
 
-    !< surface forcing 
-    type(timeForcingDataType) :: sensible_hflux_surf, latent_hflux_surf     !< heat fluxes, [W/m^2] 
-    type(timeForcingDataType) :: salin_flux_surf                            !< salinity flux, [PSU*m/s] 
-    type(timeForcingDataType) :: tau_x_surf, tau_y_surf                     !< momentum flux, [N/m**2]
-
-    type(timeForcingDataType) :: Ua, Va                                  !< air wind, [m/s]
-    type(timeForcingDataType) :: Ta                                      !< air temperature, [K]
-    type(timeForcingDataType) :: Pa                                      !< air pressure, [Pa]
-    type(timeForcingDataType) :: Qa                                      !< air humidity, [kg/kg]
-    type(timeForcingDataType) :: RHa                                     !< air relative humidity, [%]
-
-    type(timeForcingDataType) :: sw_flux_surf                    !< shortwave radiation flux IN, [W/m^2]
-    type(timeForcingDataType) :: lw_in_surf, lw_net_surf         !< longwave radiation flux IN/NET, [W/m^2] 
-
-    integer :: is_meteo_setup
-
-    !< bottom forcing
-    type(timeForcingDataType) :: hflux_bot                  !< heat flux, [W/m^2]
-    type(timeForcingDataType) :: salin_flux_bot             !< salinity flux, [PSU*m/s] 
-    type(timeForcingDataType) :: tau_x_bot, tau_y_bot       !< momentum flux, [N/m**2]
-
     !< boundary conditions data
     type(oblBcType) :: bc
 
@@ -108,16 +77,11 @@ program obl_main
 
     real, parameter :: Cd0 = 0.001           ! bottom drag coefficient
 
-    real, parameter :: lw_albedo = 0.03
-    real, parameter :: surface_emissivity = 0.97
 	 
 
     real :: mld !< mixed layer depth, [m]
     real :: lab_mld !< theoretical mixed layer depth, [m]
 
-    !< just a temporary to hold value 
-    real :: fvalue
-
     real :: N2_0
     
     ! command line arguments
@@ -392,165 +356,12 @@ program obl_main
         
         !< define fluxes & dynamic scales [surface]
         ! ----------------------------------------------------------------------------
-        if (is_meteo_setup == 0) then
-
-            !< heat flux  
-            bc%heat_fluxH = 0.0
-            call get_value_tforcing(fvalue, time_current, sensible_hflux_surf)
-            bc%heat_fluxH = bc%heat_fluxH + fvalue
-            call get_value_tforcing(fvalue, time_current, latent_hflux_surf)
-            bc%heat_fluxH = bc%heat_fluxH + fvalue
-
-            if (lw_net_surf%num > 0) then
-                call get_value_tforcing(fvalue, time_current, lw_net_surf)
-                bc%heat_fluxH = bc%heat_fluxH + fvalue 
-            else if (lw_in_surf%num > 0) then
-                call get_value_tforcing(fvalue, time_current, lw_in_surf)
-                bc%heat_fluxH = bc%heat_fluxH + (1.0 - lw_albedo) * fvalue
-
-                !< adding LW outgoing flux to balance
-                block
-                    real :: Theta_surf 
-                    real :: lw_out_surf
-                
-                    Theta_surf = Theta_dev(grid%cz) + Theta_ref
-                    lw_out_surf = surface_emissivity * stefan_boltzmann_const * & 
-                        Theta_surf * Theta_surf * Theta_surf * Theta_surf
-
-                    bc%heat_fluxH = bc%heat_fluxH - lw_out_surf
-                end block  
-            endif
-
-            !< kinematic heat flux = F / (rho_ref * cp)
-            bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water) 
-        
-            !< momentum flux 
-            call get_value_tforcing(bc%u_momentum_fluxH, time_current, tau_x_surf)
-            call get_value_tforcing(bc%v_momentum_fluxH, time_current, tau_y_surf)
-
-            !< kinematic momentum flux = tau / rho_ref
-            bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
-            bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
-
-        else if (is_meteo_setup == 1) then
-#ifdef USE_SFX
-            block
-                real :: Ua_t, Va_t, Ta_t, Pa_t, Qa_t, RHa_t, e_sat_a
-                real :: Theta_surf, Q_surf, e_sat
-
-                real :: rho_a 
-                real :: Lv
-
-                rho_a = 1.22
-                Lv = 2.5 * 1e6
-
-                call get_value_tforcing(Ua_t, time_current, Ua)
-                call get_value_tforcing(Va_t, time_current, Va)
-                call get_value_tforcing(Ta_t, time_current, Ta)
-                call get_value_tforcing(Pa_t, time_current, Pa)
-                if (Qa%num > 0) then 
-                    call get_value_tforcing(Qa_t, time_current, Qa)
-                else if (RHa%num > 0) then 
-                    call get_value_tforcing(RHa_t, time_current, RHa)
-
-                    call get_water_saturation_vapor_pressure_in_kPa(e_sat_a, Ta_t)
-                    call get_specific_humidity(Qa_t, e_sat_a, Pa_t / 1000.0, R_d, R_v)
-
-                    Qa_t = Qa_t * (RHa_t / 100.0)
-                endif
-
-                Theta_surf = Theta_dev(grid%cz) + Theta_ref
-                call get_water_saturation_vapor_pressure_in_kPa(e_sat, Theta_surf)
-                call get_specific_humidity(Q_surf, e_sat, Pa_t / 1000.0, R_d, R_v)
-
-                meteo%U = sqrt(Ua_t * Ua_t + Va_t * Va_t)
-                meteo%dT = Ta_t - Theta_surf
-                meteo%Tsemi = 0.5 * (Ta_t + Theta_surf)
-                meteo%dQ = Qa_t - Q_surf
-                meteo%h = 2.0
-                meteo%z0_m = -1.0
-
-                call get_surface_fluxes_most(sfx, meteo, sfx_numerics)
-
-                !< heat flux  
-                bc%heat_fluxH = &
-                    rho_a * cp_air * sfx%Cm * sfx%Ct * meteo%U * meteo%dT + &
-                    rho_a * Lv * sfx%Cm * sfx%Ct * meteo%U * meteo%dQ 
-
-                if (lw_net_surf%num > 0) then
-                    call get_value_tforcing(fvalue, time_current, lw_net_surf)
-                    bc%heat_fluxH = bc%heat_fluxH + fvalue 
-                else if (lw_in_surf%num > 0) then
-                    call get_value_tforcing(fvalue, time_current, lw_in_surf)
-                    bc%heat_fluxH = bc%heat_fluxH + (1.0 - lw_albedo) * fvalue
-        
-                    !< adding LW outgoing flux to balance
-                    block
-                        real :: Theta_surf 
-                        real :: lw_out_surf
-                        
-                        Theta_surf = Theta_dev(grid%cz) + Theta_ref
-                        lw_out_surf = surface_emissivity * stefan_boltzmann_const * & 
-                            Theta_surf * Theta_surf * Theta_surf * Theta_surf
-        
-                        bc%heat_fluxH = bc%heat_fluxH - lw_out_surf
-                    end block  
-                endif
-        
-                !< kinematic heat flux = F / (rho_ref * cp)
-                bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water) 
-
-                !< momentum flux 
-                bc%u_momentum_fluxH = rho_a * sfx%Cm * sfx%Cm * meteo%U * Ua_t
-                bc%v_momentum_fluxH = rho_a * sfx%Cm * sfx%Cm * meteo%U * Va_t
-
-                !< kinematic momentum flux = tau / rho_ref
-                bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
-                bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
-
-            end block 
-#endif
-        endif
-
-        !< salinity flux 
-        call get_value_tforcing(bc%salin_fluxH, time_current, salin_flux_surf)
-
-        !< shortwave radiation flux [W/m^2]
-        call get_value_tforcing(bc%sw_fluxH, time_current, sw_flux_surf)
-
-        !< U* def.:
-        call get_u_dynamic(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)
-
-        !< rho* def.:
-        call get_rho_dynamic(bc%rho_dynH, &
-            bc%U_dynH, bc%heat_fluxH, bc%salin_fluxH) 
+        call advance_surface_fluxes(bc, time_current, grid)
         ! ----------------------------------------------------------------------------
 
         !< define fluxes & dynamic scales [bottom]
         ! ----------------------------------------------------------------------------
-        !< heat flux  
-        call get_value_tforcing(bc%heat_flux0, time_current, hflux_bot)
-
-        !< kinematic heat flux = F / (rho_ref * cp)
-        bc%heat_flux0 = bc%heat_flux0 / (Rho_ref * cp_water)
-
-        !< salinity flux 
-        call get_value_tforcing(bc%salin_flux0, time_current, salin_flux_bot)
-
-        !< momentum flux 
-        call get_value_tforcing(bc%u_momentum_flux0, time_current, tau_x_bot)
-        call get_value_tforcing(bc%v_momentum_flux0, time_current, tau_y_bot)
-
-        !< kinematic momentum flux = tau / rho_ref
-        bc%u_momentum_flux0 = bc%u_momentum_flux0 / Rho_ref
-        bc%v_momentum_flux0 = bc%v_momentum_flux0 / Rho_ref
-
-        !< U* def.:
-        call get_u_dynamic(bc%U_dyn0, bc%u_momentum_flux0, bc%v_momentum_flux0)
-
-        !< rho* def.:
-        call get_rho_dynamic(bc%rho_dyn0, &
-            bc%U_dyn0, bc%heat_flux0, bc%salin_flux0) 
+        call advance_bottom_fluxes(bc, time_current, grid)
         ! ----------------------------------------------------------------------------
 
         !< advance turbulence closure
@@ -634,34 +445,8 @@ program obl_main
     !> deallocate scm
     call deallocate_scm_vec
 
-    !> removing time-dependent forcing data
-    call deallocate_tforcing(sensible_hflux_surf)
-    call deallocate_tforcing(latent_hflux_surf)
-
-    call deallocate_tforcing(salin_flux_surf)
-
-    call deallocate_tforcing(tau_x_surf)
-    call deallocate_tforcing(tau_y_surf)
-
-    call deallocate_tforcing(hflux_bot)
-    call deallocate_tforcing(salin_flux_bot)
-
-    call deallocate_tforcing(tau_x_bot)
-    call deallocate_tforcing(tau_y_bot)
-
-    call deallocate_tforcing(Ua)
-    call deallocate_tforcing(Va)
-
-    call deallocate_tforcing(Ta)
-    call deallocate_tforcing(Pa)
-    call deallocate_tforcing(Qa)
-    call deallocate_tforcing(RHa)
-
-    call deallocate_tforcing(sw_flux_surf)
-
-    call deallocate_tforcing(lw_in_surf)
-    call deallocate_tforcing(lw_net_surf)
-
+    !> deallocate time-dependent forcing
+    call deallocate_fluxes_data
     
     !> removing time slice data
     call output_cleanup(scm_output)