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
    use obl_state

    use obl_tforcing
    use obl_tslice
    use obl_tseries
    use obl_output
    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_boundary
    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

    !< 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]

    !< 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

    real :: t !< time, [s]
    integer :: nt !< number of time steps 
    real :: z !< depth, [m]
    integer :: nz !< number of steps in z (depth), [m]
    real :: dt !< time step [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 :: 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

    integer, parameter :: output_mode = 3   ! 1 -- netcdf, 2 -- ascii, 3 -- tecplot

    integer, parameter :: obl_setup = 1     ! 1 - Kato-Phillips, 2 - Papa station
	 

    real :: mld !< mixed layer depth, [m]
    real :: lab_mld !< theoretical mixed layer depth, [m]

    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]

    !< just a temporary to hold value 
    real :: fvalue

    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



    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)
    write (*,*) t, dt, nt, z, nz

    ! 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
    call allocate_state_vec(nz)


    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_TKE(1:nz))  
    allocate(Km_eps(1:nz))
    


    ! 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)

    !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

    !< setting atmospheric forcing
    ! ---------------------------------------------------------------------------- 
    if (obl_setup == 1) then
        call set_const_tforcing(sensible_hflux_surf, 0.0)
        call set_const_tforcing(latent_hflux_surf, 0.0)

        call set_const_tforcing(salin_flux_surf, 0.0)

        call set_const_tforcing(tau_x_surf, 0.1)
        call set_const_tforcing(tau_y_surf, 0.0)

    endif 
    if (obl_setup == 2) then
        call set_external_tforcing(sensible_hflux_surf, 'papa-2017/sensible_hflux.txt')
        call set_external_tforcing(latent_hflux_surf, 'papa-2017/latent_hflux.txt')

        call set_const_tforcing(salin_flux_surf, 0.0)

        call set_external_tforcing(tau_x_surf, 'PAPA_06_2017/tau-x.txt')
        call set_external_tforcing(tau_y_surf, 'PAPA_06_2017/tau-y.txt')
    endif
    ! ----------------------------------------------------------------------------

    !< setting bottom forcing
    ! ----------------------------------------------------------------------------
    call set_const_tforcing(hflux_bot, 0.0)

    call set_const_tforcing(salin_flux_bot, 0.0)

    call set_const_tforcing(tau_x_bot, 0.0)
    call set_const_tforcing(tau_y_bot, 0.0)
    ! ----------------------------------------------------------------------------

    !set right-hand side
    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

    i=1

    do while (time_current < time_end )
        ! ----------------------------------------------------------------------------
        
        !< define fluxes & dynamic scales [surface]
        ! ----------------------------------------------------------------------------
        !< 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

        !< kinematic heat flux = F / (rho_ref * cp)
        bc%heat_fluxH = bc%heat_fluxH / (Rho_ref * cp_water) 

        !< salinity flux 
        call get_value_tforcing(bc%salin_fluxH, time_current, salin_flux_surf)
        
        !< 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;

        !< U* def.:
        call get_dyn_velocity(bc%U_dynH, bc%u_momentum_fluxH, bc%v_momentum_fluxH)

        !< rho* def.:
        call get_rho_dyn(bc%rho_dynH, &
            bc%u_momentum_fluxH, bc%v_momentum_fluxH, bc%heat_fluxH, bc%salin_fluxH) 
        ! ----------------------------------------------------------------------------

        !< 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;

        !< rho* def.:
        call get_rho_dyn(bc%rho_dyn0, &
            bc%u_momentum_flux0, bc%v_momentum_flux0, bc%heat_flux0, bc%salin_flux0) 
        ! ----------------------------------------------------------------------------
        
        call solve_state_eq(Rho, Theta_dev, Salin_dev, nz)
#ifndef OBL_USE_GRID_DATATYPE
        call get_N2(N2, Rho, bc%Rho_dynH, bc%rho_dyn0, dz, nz) 
#else
        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, dz, nz)
#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, nz)

        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, nz)
            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
        else if (closure_mode.eq.2) then
            call get_Km_plus(Km, Ri_grad, &
                bc%u_momentum_fluxH, bc%v_momentum_fluxH, mld, nz, dz, z)
            call get_Kh_plus(Kh, Km, nz)
            call get_TKE_generation(TKE_production, Km, S2, nz)
            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
        else if (closure_mode.eq.4) then
            call TKE_bc(TKE, &
                bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz)
            call eps_bc(EPS, &
                bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz, dz)
            call get_Km_k_eps (Km, TKE, EPS, nz)
            call get_Kh_k_eps (Kh, Km, nz)
            call get_TKE_generation(TKE_production, Km, S2, nz)
            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, nz)
            call TKE_bc(TKE, &
                bc%u_momentum_fluxH, bc%v_momentum_fluxH, nz)
            call eps_bc(EPS, &
                bc%u_momentum_fluxH, bc%v_momentum_fluxH, 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, TKE_production, TKE_buoyancy, eps)
            call solve_TKE_eq (TKE, Km_TKE, nz, dz, dt, TKE_production, TKE_buoyancy, EPS)
            call limit_min_array(TKE, TKE_min, nz)
            call solve_eps_eq_semiimplicit (EPS, Km_eps, nz, dz, dt, TKE_production, TKE_buoyancy, TKE)
            call limit_min_array(EPS, eps_min, nz) 
        endif

        call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, &
            bc%heat_fluxH, bc%heat_flux0, f_heat_right) 
        call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, &
            bc%salin_fluxH, bc%salin_flux0, f_sal_right)
        call solve_vector_eq (U, Km, nz, dz, dt, &
            bc%u_momentum_fluxH, bc%u_momentum_flux0, f_u_right)
        call solve_vector_eq (V, Km, nz, dz, dt, &
            bc%v_momentum_fluxH, bc%v_momentum_flux0, f_v_right)

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

        !> advance screen output
        if (mod(i, nscreen) == 0) then
            call get_mld(mld, N2, dz, nz, z)
            call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, z)

            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

        !> advance file output
        if (mod(i, noutput) == 0) then
            call push_state_vec(nz)

            call get_mld(mld, N2, dz, nz, z)
            call get_lab_mld(lab_mld, bc%U_dynH, N2_0, time_current, z)

            call push_value_to_tseries(output_mld, mld)
            call push_value_to_tseries(output_mld_ref, lab_mld)

            call push_value_to_tseries(output_tau_x, bc%u_momentum_fluxH)
            call push_value_to_tseries(output_tau_y, bc%v_momentum_fluxH)

            call push_value_to_tseries(output_Theta_surf, Theta_dev(grid%cz))

            call push_value_to_tseries(output_time, time_current / 3600.0)
        endif
    enddo

    if (output_mode > 0) then
        output_Theta%data(:,1:output_Theta%num) = output_Theta%data(:,1:output_Theta%num) + Theta_ref
        output_Salin%data(:,1:output_Salin%num) = output_Salin%data(:,1:output_Salin%num) + Salin_ref

        output_TinC%data(:,1:output_TinC%num) = output_TinC%data(:,1:output_TinC%num) + Theta_ref - 273.15

        output_tau_x%data(1:output_tau_x%num) = output_tau_x%data(1:output_tau_x%num) * Rho_ref
        output_tau_y%data(1:output_tau_y%num) = output_tau_y%data(1:output_tau_y%num) * Rho_ref

        output_Theta_surf%data(1:output_Theta_surf%num) = output_Theta_surf%data(1:output_Theta_surf%num) + Theta_ref
    endif

    if (output_mode.eq.1) then
        write(*, *) ' >> writing netcdf output ...'
        call write_netcdf
    endif

    if (output_mode.eq.2) then
        write(*, *) ' >> writing ascii output ...'
        call write_ascii
    endif


    if (output_mode.eq.3) then
        write(*, *) ' >> writing tecplot output ...'
        call write_tecplot(grid%z, output_time%data)
    endif




    !> deallocation
    call deallocate_state_vec

    deallocate(dz)

    deallocate(f_Heat_right)
    deallocate(f_Sal_right)
    deallocate(f_u_right)
    deallocate(f_v_right)
    
    deallocate(Km_TKE)  
    deallocate(Km_eps)

    !> 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)

    
    !> removing time slice data
    call output_cleanup

    ! > removing grid data
    call deallocate_grid(grid)


end program