obl_main.f90

#include "obl_def.fi"

program obl_main
    !< @brief main program for calculations for ocean boundary layer

#ifdef USE_CONFIG_PARSER
        use iso_c_binding, only: C_NULL_CHAR
        use config_parser
#endif

    ! modules used
    use obl_grid
#ifdef OBL_USE_TSLICE_OUTPUT
    use obl_tslice
#endif
    use obl_time_and_space
    use obl_initial_conditions
    use obl_forcing_and_boundary
    use obl_rhs
    use obl_equation
    use obl_turb_closure
    use obl_diag
    use obl_pacanowski_philander
    use obl_pacanowski_philander_plus
    use obl_k_epsilon
    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
    !use vertical_mixing, default = off
    !use vermix

    ! directives list
    implicit none

    type(gridDataType) :: grid

    ! turbulence closure parameters
    type(pacanowskiParamType) :: param_pacanowski

#ifdef OBL_USE_TSLICE_OUTPUT
    ! time slice
    type (oblTimeSlice) :: output_Theta, output_Salin
#endif

    real :: t !< time, [s]
    integer :: nt !< number of time steps 
    integer :: nf !< number of time steps for forcing
    real :: z !< depth, [m]
    integer :: nz !< number of steps in z (depth), [m]
    real :: dt !< time step [s] 
    real :: df !< time step for forcing [s] 
    real, allocatable :: dz(:) !< depth(z) step [m] 
    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] ???
    real, allocatable :: Times_flux_Sal_surf(:), Times_flux_Sal_bot(:) 
    real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:) !< [N/m**2]   
    real, allocatable :: flux_u_surf(:), flux_u_bot(:) !< [(m/s)**2]    
    real, allocatable :: Times_flux_u_surf(:), Times_flux_u_bot(:)       
    real, allocatable :: flux_v_surf(:), flux_v_bot(:) !< [(m/s)**2]       
    real, allocatable :: Times_flux_v_surf(:), Times_flux_v_bot(:)      
    real, allocatable :: f_Heat_right(:), f_Sal_right(:), f_u_right(:), f_v_right(:) !< RHS    

    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]
    real, allocatable :: eps(:) !< TKE dissipation rate, [W/kg]   
    real, allocatable :: P_TKE(:)  !< shear production of TKE, [m**2 / s**3]
    real, allocatable :: B_TKE(:)  !< buoyancy production of TKE, [m**2 / s**3]
    real, allocatable :: Km_TKE(:)  !eddy viscosity for momentum adjusted for Schmidt TKE number, [m**2 / s]
    real, allocatable :: Km_eps(:)  !eddy viscosity for momentum adjusted for Schmidt eps number, [m**2 / s]

    real :: u_dyn0, u_dynH
    real :: flux_heat_surf_res, flux_heat_bot_res, flux_sal_surf_res, flux_sal_bot_res
    real :: flux_u_surf_res, flux_u_bot_res, flux_v_surf_res, flux_v_bot_res 
 
    real :: N2_0
    
    ! command line arguments
    ! --------------------------------------------------------------------------------
    integer :: num_args
    character(len = 128) :: arg

    character(len = 128), parameter :: arg_key_help = '--help'
    character(len = 128), parameter :: arg_key_config = "--config"

    integer :: ierr
    ! --------------------------------------------------------------------------------

    ! screen output parameters
    integer, parameter :: nscreen = 1000

    ! file output parameters
    integer, parameter :: noutput = 60

    ! arrays for output
    real, allocatable           :: time_hrs(:)
    real, allocatable           :: time_hrs_tslice(:)

    real, allocatable           :: Theta_write(:,:)
    real, allocatable           :: Salin_write(:,:)
    real, allocatable           :: U_write(:,:)
    real, allocatable           :: V_write(:,:)
    real, allocatable           :: Rho_write(:,:)
    real, allocatable           :: N2_write(:,:)
    real, allocatable           :: S2_write(:,:)
    real, allocatable           :: Ri_grad_write(:,:)
    real, allocatable           :: Kh_write(:,:)
    real, allocatable           :: Km_write(:,:)
    real, allocatable           :: mld_write(:)
    real, allocatable           :: Theta_C_write(:,:)
    real, allocatable           :: lab_mld_write(:)

    closure_mode = 4 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom


    ! --- command line arguments processing
    num_args = command_argument_count()
    do i = 1, num_args

        call get_command_argument(i, arg)
        if (trim(arg) == trim(arg_key_help)) then
            write(*, *) ' obl model, usage:'
            write(*, *) ' --help'
            write(*, *) '    print usage options'
            write(*, *) ' --config [filename]'
            write(*, *) '    use configuration file'
            return
        end if

        if (trim(arg) == trim(arg_key_config)) then
            if (i == num_args) then
                write(*, *) ' FAILURE! > missing configuration file [key] argument'
                ierr = 1        ! signal ERROR
                return
            end if
        
            call get_command_argument(i + 1, arg)
#ifdef USE_CONFIG_PARSER
            call c_config_run(trim(arg)//C_NULL_CHAR, status)
            if (status == 0) then
                write(*, *) ' FAILURE! > unable to parse configuration file: ', trim(arg)
                ierr = 1        ! signal ERROR
                return
            end if

            call c_config_is_varname("domain.depth"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_float("domain.depth"//C_NULL_CHAR, z, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            end if

            call c_config_is_varname("grid.cz"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_int("grid.cz"//C_NULL_CHAR, nz, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            end if

            call c_config_is_varname("time.end"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_float("time.end"//C_NULL_CHAR, t, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            end if

            call c_config_is_varname("time.n"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_int("time.n"//C_NULL_CHAR, nt, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            end if
#endif
        endif
    enddo
    
#ifndef USE_CONFIG_PARSER
    ! set time, depth, steps and number of space
    call set_Time_Space(t, nt, z, nz)
#endif
    allocate(dz(1:nz))
    call set_dz_dt(dz, dt, nz, nt, t, z)
    call set_Time_Space_Forcing(nf, df) 
    write (*,*) t, dt, nt, z, nz, nf, df

    ! setting grid -- just an example
    !   > [zpos, height, cz, gcz]
    call set_uniform_grid(grid, 0.0, z, nz)

    ! debug grid print
    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))
    allocate(flux_Heat_surf(1:nf))
    allocate(Times_flux_Heat_surf(1:nf))
    allocate(flux_Heat_bot(1:nf))
    allocate(Times_flux_Heat_bot(1:nf))
    allocate(flux_Sal_surf(1:nf))
    allocate(Times_flux_Sal_surf(1:nf))
    allocate(flux_Sal_bot(1:nf))
    allocate(Times_flux_Sal_bot(1:nf))
    allocate(Tau_x_surf(1:nf))
    allocate(Tau_y_surf(1:nf))
    allocate(Tau_x_bot(1:nf))
    allocate(Tau_y_bot(1:nf))
    allocate(flux_u_surf(1:nf))
    allocate(Times_flux_u_surf(1:nf))
    allocate(flux_u_bot(1:nf))
    allocate(Times_flux_u_bot(1:nf))
    allocate(flux_v_surf(1:nf))
    allocate(Times_flux_v_surf(1:nf))
    allocate(flux_v_bot(1:nf))
    allocate(Times_flux_v_bot(1:nf))
    allocate(f_Heat_right(1:nz))
    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)) 
    allocate(B_TKE(1:nz)) 
    allocate(Km_TKE(1:nz))  
    allocate(Km_eps(1:nz))
    
    allocate(time_hrs(1:nt))
    allocate(time_hrs_tslice(1:nt))
    allocate(Theta_write(nz, nt))
    allocate(Salin_write(nz, nt))
    allocate(U_write(nz, nt))
    allocate(V_write(nz, nt))
    allocate(Rho_write(nz, nt))
    allocate(N2_write(nz, nt))
    allocate(S2_write(nz, nt))
    allocate(Ri_grad_write(nz, nt))
    allocate(Kh_write(nz, nt))
    allocate(Km_write(nz, nt))
    allocate(mld_write(nt))
    allocate(Theta_C_write(nz, nt))
    allocate(lab_mld_write(nt))

    ! initialization of main fields
    call Theta_init (Theta, Theta_dev, nz, dz)
    !call Theta_init (Theta, nz, dz)
    call Salin_init(Salin, Salin_dev, nz, dz)
    call U_init(U, nz)
    call V_init(V, nz)

    !set boundary conditions
    call set_Flux_heat_surf (Flux_heat_surf, Times_flux_heat_surf, nf, df)
    call set_Flux_heat_bot (Flux_heat_bot, Times_flux_heat_bot, nf, df) 
    call set_Flux_sal_surf (Flux_sal_surf, Times_flux_sal_surf, nf, df)
    call set_Flux_sal_bot (Flux_sal_bot, Times_flux_sal_bot, nf, df)
    call set_Tau_x_surf(Tau_x_surf, flux_u_surf, Times_flux_u_surf, nf, df)
    call set_Tau_y_surf(Tau_y_surf, flux_v_surf, Times_flux_v_surf, nf, df)
    call set_Tau_x_bot(Tau_x_bot, flux_u_bot, Times_flux_u_bot, nf, df)
    call set_Tau_y_bot(Tau_y_bot, flux_v_bot, Times_flux_v_bot, nf, df)
    call set_f_Heat_right (f_Heat_right, nz)  
    call set_f_Sal_right (f_Sal_right, nz) 
    call set_f_u_right (f_v_right, V, nz)
    call set_f_v_right (f_v_right, U, nz)

    
    
    !open (199, file= 'output_Daria/mld.txt', status ='replace')
    !open (20, file= 'output_Daria/surf_temp.txt', status ='replace')

    status = 0
    num = 0

    time_begin = 0.0
    time_end = t
    time_current = time_begin

    N2_0 = 0.0004

    !initialization of TKE & eps in case of k-epsilon closure
    if (closure_mode.eq.3 .or. closure_mode.eq.4) then
        call TKE_init(TKE, nz)
        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
            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) 
            call get_value_interpolate(flux_sal_bot_res, time_current, flux_sal_bot, Times_flux_sal_bot, nf)
            call get_value_interpolate(flux_u_surf_res, time_current, flux_u_surf, Times_flux_u_surf, nf) 
            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_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) 
#else
            call get_N2_on_grid(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, grid)
#endif
#ifndef OBL_USE_GRID_DATATYPE
            call get_S2(S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, dz, nz)
#else
            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)

                call get_TKE_generation(P_TKE, Km, S2, nz)
                call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
            else if (closure_mode.eq.2) then
                call get_Km_plus(Km, Ri_grad, flux_u_surf_res, flux_v_surf_res, mld, nz, dz, z)
                call get_Kh_plus(Kh, Km, nz)
                call get_TKE_generation(P_TKE, Km, S2, nz)
                call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
            else if (closure_mode.eq.4) then
                call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
                call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, nz, dz)
                call get_Km_k_eps (Km, TKE, eps, nz)
                call get_Kh_k_eps (Kh, Km, nz)
                call get_TKE_generation(P_TKE, Km, S2, nz)
                call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
                call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
                call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, nz, dz)
                call get_Km_TKE(Km_TKE, Km, nz)
                call get_Km_eps(Km_eps, Km, nz)
                !call solve_TKE_eq_semiimplicit (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps)
                call solve_TKE_eq (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps)
                call limit_min_array(TKE, TKE_min, nz)
                call solve_eps_eq_semiimplicit (eps, Km_eps, nz, dz, dt, P_TKE, B_TKE, TKE)
                call limit_min_array(eps, eps_min, nz) 
            endif
            
            call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf_res, flux_heat_bot_res, f_heat_right) 
            call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf_res, flux_sal_bot_res, f_sal_right)
            call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf_res, flux_u_bot_res, f_u_right)
            call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf_res, flux_v_bot_res, f_v_right)
        !else if (closure_mode.eq.5) then
        !    call solve_state_eq  (Rho, Theta_dev, Salin_dev, nz)
        !    call vermix_init(rho, u, v, Tau_x_surf(i), Tau_y_surf(i), Ri_grad, Kh, Km)
        !    call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right) 
        !    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)
        endif

        !do k = 1, nz
        !    Theta (k) = Theta(k) + Theta_0
        !end do
        
        Theta_write(:, i) = Theta_dev
        Salin_write(:, i) = Salin_dev
        U_write(:, i) = U
        V_write(:, i) = V
        Rho_write(:, i) = Rho !Rho_dev
        N2_write(:, i) = N2
        S2_write(:, i) = S2
        Ri_grad_write(:, i) = Ri_grad
        Kh_write(:, i) = Kh
        Km_write(:, i) = Km
        mld_write(i) = mld
       !lab_mld_write(i) = lab_mld

        !write(199,*) time_current, mld, lab_mld
        !write(20,*) time_current, Theta_dev(nz) + Theta_ref

        !> advancing time
        i = i + 1
        time_current = time_current + dt

        time_hrs(i) = time_current / 3600.0

        if (mod(i, nscreen) == 0) then
            write(*, '(a,g0)') ' mld = ', mld
            write(*, '(a,g0)') ' Theta(surface) = ', Theta_dev(grid%cz) + Theta_ref
            write(*, '(a,g0,a,g0,a)') ' current time = ', time_current / 3600.0, ' HRS [ ', & 
                (time_current / time_end) * 100.0, '% ]' 
            write(*, '(a)') '-------------------------------------------------'
        endif

        ! time slice output
#ifdef OBL_USE_TSLICE_OUTPUT
        if (mod(i, noutput) == 0) then
            call push_profile_to_tslice(output_Theta, Theta_dev, nz)
            call push_profile_to_tslice(output_Salin, Salin_dev, nz)

            time_hrs_tslice(output_Theta%num) = time_current / 3600.0
        endif
#endif
    enddo

    Theta_C_write = Theta_write - 273.15
    Theta_write = Theta_write + Theta_ref
    Salin_write = Salin_write + Salin_ref
    !Rho_write = Rho_write + Rho_ref

#ifndef OBL_USE_TSLICE_OUTPUT
    ! Writing 2D arrays (variables with depth and time)
    call write_2d_real_nc(Theta_write, 'output.nc', meta_theta)
    call write_2d_real_nc(Salin_write, 'output.nc', meta_salin)
    call write_2d_real_nc(U_write, 'output.nc', meta_u)
    call write_2d_real_nc(V_write, 'output.nc', meta_v)
    call write_2d_real_nc(Rho_write, 'output.nc', meta_rho)
    call write_2d_real_nc(N2_write, 'output.nc', meta_n2)
    call write_2d_real_nc(S2_write, 'output.nc', meta_s2)
    call write_2d_real_nc(Ri_grad_write, 'output.nc', meta_ri_grad)
    call write_2d_real_nc(Kh_write, 'output.nc', meta_kh)
    call write_2d_real_nc(Km_write, 'output.nc', meta_km)
    call write_2d_real_nc(Theta_C_write, 'output.nc', meta_theta_C)
#endif

    ! Writing 1D arrays (scalar variables over time)
    call write_1d_real_nc(mld_write, 'output.nc', meta_mld)
    !call write_1d_real_nc(Tau_x_surf, 'output.nc', meta_tau_u)
    !call write_1d_real_nc(Tau_y_surf, 'output.nc', meta_tau_v)
    call write_1d_real_nc(lab_mld_write, 'output.nc', meta_lab_mld)

#ifndef OBL_USE_TSLICE_OUTPUT
    ! Writing 2D arrays (variables with depth and time) to ASCII files
    call write_2d_real_ascii(Theta_write, 'output_Theta.txt', meta_theta)
    call write_2d_real_ascii(Salin_write, 'output_Salin.txt', meta_salin)
    call write_2d_real_ascii(U_write, 'output_U.txt', meta_u)
    call write_2d_real_ascii(V_write, 'output_V.txt', meta_v)
    call write_2d_real_ascii(Rho_write, 'output_Rho.txt', meta_rho)
    call write_2d_real_ascii(N2_write, 'output_N2.txt', meta_n2)
    call write_2d_real_ascii(S2_write, 'output_S2.txt', meta_s2)
    call write_2d_real_ascii(Ri_grad_write, 'output_Ri_grad.txt', meta_ri_grad)
    call write_2d_real_ascii(Kh_write, 'output_Kh.txt', meta_kh)
    call write_2d_real_ascii(Km_write, 'output_Km.txt', meta_km)
#endif
    
    ! Writing 1D arrays (scalar variables over time) to ASCII files
    !call write_1d_real_ascii(mld_write, 'output_mld.txt', meta_mld)
    !call write_1d_real_ascii(lab_mld_write, 'output_lab_mld.txt', meta_lab_mld)
    !call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
    !call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)


    ! Writing 1D arrays (scalar variables over time) to Tecplot files
    call write_1d_real_plt(mld_write, time_hrs, 'output_mld.plt', meta_mld)
    call write_1d_real_plt(lab_mld_write, time_hrs, 'output_lab_mld.plt', meta_lab_mld)
    !call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
    !call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)

#ifdef OBL_USE_TSLICE_OUTPUT
    ! time slice output
    call write_2d_real_ascii(output_Theta%data(:,1:output_Theta%num), 'Theta_tslice.txt', meta_theta)
    call write_2d_real_ascii(output_Salin%data(:,1:output_Salin%num), 'Salin_tslice.txt', meta_salin)

    call write_2d_real_plt(output_Theta%data(:,1:output_Theta%num), grid%z, time_hrs_tslice(1:output_Theta%num), 'Theta_tslice.plt', meta_theta)
#endif

    !close(199)
    !close(20)

    ! deallocation
    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)
    deallocate(Times_flux_Heat_bot)
    deallocate(flux_Sal_surf)
    deallocate(Times_flux_Sal_surf)
    deallocate(flux_Sal_bot)
    deallocate(Times_flux_Sal_bot)
    deallocate(Tau_x_surf)
    deallocate(Tau_y_surf)
    deallocate(Tau_x_bot)
    deallocate(Tau_y_bot)
    deallocate(flux_u_surf)
    deallocate(Times_flux_u_surf)
    deallocate(flux_u_bot)
    deallocate(Times_flux_u_bot)
    deallocate(flux_v_surf)
    deallocate(Times_flux_v_surf)
    deallocate(flux_v_bot)
    deallocate(Times_flux_v_bot)
    deallocate(f_Heat_right)
    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) 
    deallocate(B_TKE) 
    deallocate(Km_TKE)  
    deallocate(Km_eps)
    
    deallocate(time_hrs)
    deallocate(time_hrs_tslice)
    deallocate(Theta_write)
    deallocate(Salin_write)
    deallocate(U_write)
    deallocate(V_write)
    deallocate(Rho_write)
    deallocate(N2_write)
    deallocate(S2_write)
    deallocate(Ri_grad_write)
    deallocate(Kh_write)
    deallocate(Km_write)
    deallocate(mld_write)
    deallocate(Theta_C_write)
    deallocate(lab_mld_write)

    ! > removing grid data
    call deallocate_grid(grid)


end program