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_state_eq
    use obl_common_phys
    use obl_tforcing
    use obl_tslice
    use obl_tseries
    use obl_output
    use obl_init
    use obl_fluxes
    use obl_scm
    use obl_diag
    use obl_bc
    use obl_pph, only: pphParamType, &
        define_stability_functions_pph => define_stability_functions, & 
        advance_turbulence_closure_pph => advance_turbulence_closure
    use obl_pph_dyn, only: pphDynParamType, &
        define_stability_functions_pph_dyn => define_stability_functions, & 
        advance_turbulence_closure_pph_dyn => advance_turbulence_closure
    use obl_k_epsilon, only: kepsilonParamType, &
        define_stability_functions_k_epsilon => define_stability_functions, & 
        advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
        TKE_init, EPS_init, TKE_bc, EPS_bc
    use obl_io_plt
    use io
    use io_metadata
#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(pphParamType) :: param_pph
    type(pphDynParamType) :: param_pph_dyn
    type(kepsilonParamType) :: param_k_epsilon

    !< boundary conditions data
    type(oblBcType) :: bc

    real :: domain_height
    integer :: grid_cz

    !< output 
    type(oblOutputStx) :: scm_output


    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, parameter :: obl_setup = 1     ! 1 - Kato-Phillips, 2 - Papa station (fluxes), 3 - Papa station (meteo)

    real, parameter :: Cd0 = 0.001           ! bottom drag coefficient

	 

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

    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, domain_height, 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, grid_cz, 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, time_end, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            end if

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

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

    !< setting model time
    ! ----------------------------------------------------------------------------  
#ifndef USE_CONFIG_PARSER
    if (obl_setup == 1) then
        time_begin = 0.0
        time_end = 300.0 * 3600.0
        dt = 1.0

        domain_height = 100.0
        grid_cz = 32
    endif
    if ((obl_setup == 2).or.(obl_setup == 3)) then
        time_begin = 0.0
        time_end = 431.0 * 3600.0
        dt = 1.0

        domain_height = 128.0
        grid_cz = 32
    endif
#endif

    time_current = time_begin
    ! ---------------------------------------------------------------------------- 
    
    ! setting grid -- just an example
    !   > [zpos, height, cz, gcz]
    call set_uniform_grid(grid, 0.0, domain_height, grid_cz)

    ! debug grid print
    call print_grid(grid)

    ! allocation
    call allocate_state_vec(grid%cz)

    ! initialize scm
    call init_scm_vec(grid%cz)

    !< setting phys
    ! ----------------------------------------------------------------------------
    if (obl_setup == 1) then
        f = 0.0
    endif 
    if ((obl_setup == 2).or.(obl_setup == 3)) then
        f = 1.116 * 1e-4

        a_band_ratio = 0.67
        a_extinction_coeff = 1.0
        b_extinction_coeff = 1.0 / 17.0
    endif
    ! ----------------------------------------------------------------------------


    !< initialization of main fields
    ! ---------------------------------------------------------------------------- 
    if (obl_setup == 1) then
        call set_linear_profile(Theta, 330.0, 0.3, grid)
        call set_const_profile(Salin, 35.0, grid)
    endif 
    if ((obl_setup == 2).or.(obl_setup == 3)) then
        call set_external_profile(Theta, 'PAPA_06_2017/Theta.dat', grid)
        call set_external_profile(Salin, 'PAPA_06_2017/Salin.dat', grid)
    endif

    Theta_dev = Theta - Theta_ref
    Salin_dev = Salin - Salin_ref
    call get_density(Rho, Theta_dev, Salin_dev, grid%cz)

    call set_const_profile(U, 0.0, grid)
    call set_const_profile(V, 0.0, grid)

    !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, param_k_epsilon, grid%cz)
        call eps_init(EPS, param_k_epsilon, grid%cz, grid%height)
    endif
    ! ---------------------------------------------------------------------------- 

    !< setting atmospheric forcing
    ! ---------------------------------------------------------------------------- 
    is_meteo_setup = 0
    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)

        call set_const_tforcing(sw_flux_surf, 0.0)

        call set_const_tforcing(lw_net_surf, 0.0)

    endif 
    if (obl_setup == 2) then
        call set_external_tforcing(sensible_hflux_surf, 'PAPA_06_2017/sensible_hflux.dat')
        call set_external_tforcing(latent_hflux_surf, 'PAPA_06_2017/latent_hflux.dat')

        call set_const_tforcing(salin_flux_surf, 0.0)

        call set_external_tforcing(tau_x_surf, 'PAPA_06_2017/tau-x.dat')
        call set_external_tforcing(tau_y_surf, 'PAPA_06_2017/tau-y.dat')

        call set_external_tforcing(sw_flux_surf, 'PAPA_06_2017/SW-down.dat')

        call set_external_tforcing(lw_in_surf, 'PAPA_06_2017/LW-down.dat')

        !< normalize time in external forcing: hrs -> sec
        call normalize_time_tforcing(sensible_hflux_surf, 3600.0)
        call normalize_time_tforcing(latent_hflux_surf, 3600.0)

        call normalize_time_tforcing(tau_x_surf, 3600.0)
        call normalize_time_tforcing(tau_y_surf, 3600.0)

        call normalize_time_tforcing(sw_flux_surf, 3600.0)

        call normalize_time_tforcing(lw_in_surf, 3600.0)
    endif
    if (obl_setup == 3) then
        is_meteo_setup = 1

        call set_external_tforcing(Ua, 'PAPA_06_2017/u-wind.dat')
        call set_external_tforcing(Va, 'PAPA_06_2017/v-wind.dat')

        call set_const_tforcing(salin_flux_surf, 0.0)

        call set_external_tforcing(Ta, 'PAPA_06_2017/Tair.dat')
        call set_external_tforcing(Pa, 'PAPA_06_2017/Pair.dat')
        call set_external_tforcing(RHa, 'PAPA_06_2017/RHair.dat')

        call set_external_tforcing(sw_flux_surf, 'PAPA_06_2017/SW-down.dat')

        call set_external_tforcing(lw_in_surf, 'PAPA_06_2017/LW-down.dat')

        !< normalize time in external forcing: hrs -> sec
        call normalize_time_tforcing(Ua, 3600.0)
        call normalize_time_tforcing(Va, 3600.0)

        call normalize_time_tforcing(Ta, 3600.0)
        call normalize_time_tforcing(Pa, 3600.0)
        call normalize_time_tforcing(RHa, 3600.0)

        call normalize_time_tforcing(sw_flux_surf, 3600.0)

        call normalize_time_tforcing(lw_in_surf, 3600.0)
    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)
    ! ----------------------------------------------------------------------------

    !< setting reference data
    ! ----------------------------------------------------------------------------
    N2_0 = 0.00044
    ! ----------------------------------------------------------------------------

    status = 0
    num = 0

    i=1
    do while (time_current < time_end )
        ! ----------------------------------------------------------------------------
        
        !< define fluxes & dynamic scales [surface]
        ! ----------------------------------------------------------------------------
        call advance_surface_fluxes(bc, time_current, grid)
        ! ----------------------------------------------------------------------------

        !< define fluxes & dynamic scales [bottom]
        ! ----------------------------------------------------------------------------
        call advance_bottom_fluxes(bc, time_current, grid)
        ! ----------------------------------------------------------------------------

        !< advance turbulence closure
        ! ----------------------------------------------------------------------------
        if (closure_mode.eq.1) then
            call define_stability_functions_pph(param_pph, bc, grid)
            call advance_turbulence_closure_pph(param_pph, bc, grid, dt)
        else if (closure_mode.eq.2) then
            call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
            call advance_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid, dt)
        else if (closure_mode.eq.4) then
            call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
            call advance_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid, dt)
        endif
        ! ----------------------------------------------------------------------------

        !< advance single-column model eq.
        ! ----------------------------------------------------------------------------
        call advance_scm_eq(bc, grid, dt)
        ! ----------------------------------------------------------------------------

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

        !> advance screen output
        if (mod(i, 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

            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(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 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%tau_x, bc%u_momentum_fluxH * Rho_ref)
            call push_value_to_tseries(scm_output%tau_y, bc%v_momentum_fluxH * Rho_ref)

            call push_value_to_tseries(scm_output%Theta_surf, Theta_dev(grid%cz) + Theta_ref)

            call push_value_to_tseries(scm_output%time, time_current / 3600.0)
        endif
    enddo

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

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


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




    !> deallocate state
    call deallocate_state_vec

    !> deallocate scm
    call deallocate_scm_vec

    !> deallocate time-dependent forcing
    call deallocate_fluxes_data
    
    !> removing time slice data
    call output_cleanup(scm_output)

    ! > removing grid data
    call deallocate_grid(grid)


end program