Skip to content
Snippets Groups Projects
obl_main.f90 18 KiB
Newer Older
  • Learn to ignore specific revisions
  • #include "obl_def.fi"
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< @brief ocean boundary layer (OBL) main
        ! --------------------------------------------------------------------------------
    
    #ifdef USE_CONFIG_PARSER
    
        use iso_c_binding, only: C_NULL_CHAR
        use config_parser
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! modules used
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! --------------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        use obl_grid
    
        use obl_state
    
        use obl_state_eq
    
        use obl_common_phys
    
        use obl_tforcing
    
        use obl_tslice
        use obl_tseries
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        use obl_output
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        use obl_fluxes
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        use obl_config
    
        use obl_diag
    
        use obl_pph, only: pphParamType, &
    
            set_config_param_pph => set_config_param, &
    
            define_stability_functions_pph => define_stability_functions, &
            init_turbulence_closure_pph => init_turbulence_closure, & 
    
            advance_turbulence_closure_pph => advance_turbulence_closure
        use obl_pph_dyn, only: pphDynParamType, &
    
            set_config_param_pph_dyn => set_config_param, &
    
            define_stability_functions_pph_dyn => define_stability_functions, &
            init_turbulence_closure_pph_dyn => init_turbulence_closure, & 
    
            advance_turbulence_closure_pph_dyn => advance_turbulence_closure
        use obl_k_epsilon, only: kepsilonParamType, &
    
            set_config_param_k_epsilon => set_config_param, &
    
            define_stability_functions_k_epsilon => define_stability_functions, &
            init_turbulence_closure_k_epsilon => init_turbulence_closure, & 
    
            advance_turbulence_closure_k_epsilon => advance_turbulence_closure, &
            TKE_init, EPS_init, TKE_bc, EPS_bc
    
        use obl_most, only: mostParamType, &
            set_config_param_most => set_config_param, &
            define_stability_functions_most => define_stability_functions, &
            init_turbulence_closure_most => init_turbulence_closure, & 
            advance_turbulence_closure_most => advance_turbulence_closure
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! directives list
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! --------------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! model data
        ! --------------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        integer :: obl_config_id    ! --- OBL builtin configs 
                                    ! = obl_config_kato: Kato-Phillips
                                    ! = obl_config_papa_fluxes: Papa-station (fluxes)
                                    ! = obl_config_papa_meteo: Papa-station (meteo)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                                    ! = obl_config_cbl: convective boundary layer (Willis exp.)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                                    ! = obl_config_cyclone: cyclone setup
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
                                    ! = obl_config_papa_long_fluxes: Papa-station 2 years (fluxes)
    
                                    ! = obl_config_papa_long_meteo: Papa-station 2 years (meteo)
    
        integer :: obl_model_id     ! --- OBL model def.
                                    ! = obl_model_pph: pacanowski-philander
                                    ! = obl_model_pph_dyn: pacanowski-philander (dynamic) 
                                    ! = obl_model_k_epsilon: k-epsilon
    
        type(gridDataType) :: grid
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< turbulence closure parameters
    
        type(pphParamType) :: param_pph
        type(pphDynParamType) :: param_pph_dyn
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        type(kepsilonParamType) :: param_k_epsilon
    
        type(mostParamType) :: param_most
    
        !< boundary conditions data
        type(oblBcType) :: bc
    
    
        !< output 
        type(oblOutputStx) :: scm_output
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< time
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        real*8 :: time_begin_real8, time_end_real8, time_current_real8
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        real :: time_begin, time_end, time_current
        real :: dt
        integer :: time_index
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< screen output parameters
        integer, parameter :: nscreen = 1000
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< file output parameters
        integer :: output_mode      ! --- OBL output mode
                                    ! = 1 -- netcdf
                                    ! = 2 -- ascii
                                    ! = 3 -- tecplot
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        integer, parameter :: noutput = 60
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< additional variables & parameters
        real, parameter :: Cd0 = 0.001      ! bottom drag coefficient
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        real :: mld, mld_ref                    ! mixed layer depth (model & reference), [m]
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        real :: eld, eld_ref                    ! entrainment layer depth (model & reference), [m]
    
        real, parameter :: N2_ref = 0.000044    ! reference N**2 (used in mld ref. estimate), [1/s**2]
                                                ! = 0.000044 (default Kato)
                                                ! = 2.0 * 1e-4 (default CBL)
        real, parameter :: Bsurf_ref = 0.52 * 1e-7      ! reference buoyancy flux (used in eld ref. estimate)
                                                        ! = 0.52 * 1e-7 (default CBL)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! command line arguments & configuration file variables
    
        ! --------------------------------------------------------------------------------
        integer :: num_args
        character(len = 128) :: arg
    
        character(len = 128), parameter :: arg_key_help = '--help'
        character(len = 128), parameter :: arg_key_config = "--config"
    
        character(len = 128), parameter :: arg_key_model = "--model"
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        integer :: i, status
    
        integer :: ierr
        ! --------------------------------------------------------------------------------
    
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< default config = obl_config_kato (Kato-Phiilps)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< possible values:
    
        !<      = obl_config_kato
        !<      = obl_config_papa_fluxes
        !<      = obl_config_papa_meteo
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !<      = obl_config_cbl
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !<      = obl_config_cyclone
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !<      = obl_config_papa_long_fluxes
        !<      = obl_config_papa_long_meteo
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        obl_config_id = obl_config_kato
    
        !< default closure = 4, k-epsilon
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< poissible values:
    
        !<      = obl_model_pph, pacanowski-philander
        !<      = obl_model_pph_dyn, pacanowski-philander (dynamic)
        !<      = obl_model_k_epsilon, k-epsilon
    
        !<      = obl_model_most, MOST scheme
    
        obl_model_id = obl_model_k_epsilon
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
        !< default output = 1, (netcdf)
        !< possible values:
        !<      = 1, netcdf
        !<      = 2, ascii
        !<      = 3, tecplot
        output_mode = 3
    
    
        ! --- 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 [key] || [filename]'
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                write(*, *) '    builtin setup [key] = kato (default) || papa-fluxes || papa-meteo || cbl || cyclone'
    
                write(*, *) '    use configuration file [filename]'
    
                write(*, *) ' --model [key]'
    
                write(*, *) '    [key] = pph || pph-dyn || k-epsilon (default) || most'
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            if (trim(arg) == trim(arg_key_config)) then
    
                if (i == num_args) then
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                    write(*, *) ' FAILURE! > missing config [key] argument'
    
                    ierr = 1        ! signal ERROR
                    return
                end if
                
                call get_command_argument(i + 1, arg)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                obl_config_id = get_obl_config_id(arg)
                if (obl_config_id == -1) then
    
    #ifdef USE_CONFIG_PARSER
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                    !< try reading configuration file
                    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
    
                    !< forcing configuration file setup
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                    obl_config_id = 999
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    #else
                    write(*, *) ' FAILURE! > unknown config [key]: ', trim(arg)
    
                    ierr = 1        ! signal ERROR
                    return
    #endif
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                end if
    
            else if (trim(arg) == trim(arg_key_model)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing model [key] argument'
                    ierr = 1        ! signal ERROR
                    return
                end if
                
                call get_command_argument(i + 1, arg)
                obl_model_id = get_obl_model_id(arg)
                if (obl_model_id == -1) then
                    write(*, *) ' FAILURE! > unknown model [key]: ', trim(arg)
                    ierr = 1        ! signal ERROR
                    return
                end if
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
    
    
        ! ---------------------------------------------------------------------------- 
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call set_grid(grid, obl_config_id, ierr)
    
        if (ierr /= 0) then
    
            write(*, *) ' FAILURE! > unable to set grid '
    
        !< debug grid print
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call print_grid(grid)
    
        ! ---------------------------------------------------------------------------- 
    
        !< setting model time
        ! ----------------------------------------------------------------------------  
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call set_time(time_begin, time_end, dt, obl_config_id, ierr)
    
        if (ierr /= 0) then
    
            write(*, *) ' FAILURE! > unable to set time '
    
        ! ----------------------------------------------------------------------------  
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< allocating state vector
        ! ----------------------------------------------------------------------------  
    
        call allocate_state_vec(grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! ----------------------------------------------------------------------------  
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !< initialize scm
        ! ----------------------------------------------------------------------------  
    
        call init_scm_vec(grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! ----------------------------------------------------------------------------  
    
    
        !< setting phys
        ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call set_phys(obl_config_id, ierr)
    
        if (ierr /= 0) then
    
            write(*, *) ' FAILURE! > unable to set phys parameters '
    
        endif 
        ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
        !< setting forcing
        ! ---------------------------------------------------------------------------- 
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call set_forcing(obl_config_id, ierr)
    
        if (ierr /= 0) then
    
            write(*, *) ' FAILURE! > unable to set forcing '
    
            return
        endif 
        ! ----------------------------------------------------------------------------
    
        !< setting turbulence closure parameters
        ! ---------------------------------------------------------------------------- 
        call set_config_param_pph(param_pph, "turbulence_model.pph", ierr)
        if (ierr /= 0) then
            write(*, *) ' FAILURE! > unable to set ''pph'' parameters '
            return
        endif 
        call set_config_param_pph_dyn(param_pph_dyn, "turbulence_model.pph_dyn", ierr)
        if (ierr /= 0) then
            write(*, *) ' FAILURE! > unable to set ''pph-dyn'' parameters '
            return
        endif 
        call set_config_param_k_epsilon(param_k_epsilon, "turbulence_model.k_epsilon", ierr)
        if (ierr /= 0) then
            write(*, *) ' FAILURE! > unable to set ''k-epsilon'' parameters '
            return
        endif 
        ! ----------------------------------------------------------------------------
    
    
        !< initialization of main fields
        ! ---------------------------------------------------------------------------- 
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        call set_initial_conditions(grid, obl_config_id, ierr)
    
        if (ierr /= 0) then
    
            write(*, *) ' FAILURE! > unable to set initial conditions '
    
        endif
    
        Theta_dev = Theta - Theta_ref
        Salin_dev = Salin - Salin_ref
    
        call get_density(Rho, Theta_dev, Salin_dev, grid%cz)
    
        ! ----------------------------------------------------------------------------
    
        !< initialization of turbulence closure
    
        ! ---------------------------------------------------------------------------- 
        !< define fluxes & dynamic scales at startup
        call advance_surface_fluxes(bc, time_begin, grid)
        call advance_bottom_fluxes(bc, time_begin, grid)
    
    
        if (obl_model_id.eq.obl_model_pph) then
    
            call define_stability_functions_pph(param_pph, bc, grid)
            call init_turbulence_closure_pph(param_pph, bc, grid)
    
        else if (obl_model_id.eq.obl_model_pph_dyn) then
    
            call define_stability_functions_pph_dyn(param_pph_dyn, bc, grid)
            call init_turbulence_closure_pph_dyn(param_pph_dyn, bc, grid)
    
        else if (obl_model_id.eq.obl_model_k_epsilon) then
    
            call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
            call init_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid)
    
        else if (obl_model_id.eq.obl_model_most) then
            call define_stability_functions_most(param_most, bc, grid)
            call init_turbulence_closure_most(param_most, bc, grid)
    
        ! ---------------------------------------------------------------------------- 
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        time_begin_real8 = time_begin
        time_current_real8 = time_begin
        time_end_real8 = time_end
    
    
    
        time_current = time_begin
        time_index = 1
    
        do while (time_current < time_end )
            ! ----------------------------------------------------------------------------
    
            
            !< define fluxes & dynamic scales [surface]
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            call advance_surface_fluxes(bc, time_current, grid)
    
            ! ----------------------------------------------------------------------------
    
            !< define fluxes & dynamic scales [bottom]
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            call advance_bottom_fluxes(bc, time_current, grid)
    
            ! ----------------------------------------------------------------------------
    
    
            !< advance turbulence closure
            ! ----------------------------------------------------------------------------
    
            if (obl_model_id.eq.obl_model_pph) then
    
                call define_stability_functions_pph(param_pph, bc, grid)
                call advance_turbulence_closure_pph(param_pph, bc, grid, dt)
    
            else if (obl_model_id.eq.obl_model_pph_dyn) 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 (obl_model_id.eq.obl_model_k_epsilon) then
    
                call define_stability_functions_k_epsilon(param_k_epsilon, bc, grid)
                call advance_turbulence_closure_k_epsilon(param_k_epsilon, bc, grid, dt)
    
            else if (obl_model_id.eq.obl_model_most) then
                call define_stability_functions_most(param_most, bc, grid)
                call advance_turbulence_closure_most(param_most, bc, grid, dt)
    
            ! ----------------------------------------------------------------------------
    
            !< advance single-column model eq.
    
            ! ----------------------------------------------------------------------------
    
            call advance_scm_eq(bc, grid, dt)
    
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            !> advance time
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
            time_index = time_index + 1
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
    
            time_current_real8 = time_current_real8 + dt
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            time_current = time_current + dt
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            !> advance screen output
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
            if (mod(time_index, nscreen) == 0) then
    
                call get_mld(mld, N2, grid%dz, grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
                ! *: add finite check
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                write(*, '(a,g0,a,g0)') ' mld = ', mld, ' eld = ', eld
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                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)') '-------------------------------------------------'
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            !> advance file output
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
            if (mod(time_index, noutput) == 0) then
    
                call push_state_vec(scm_output, grid%cz)
    
                call get_mld(mld, N2, grid%dz, grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                call get_mld_ref(mld_ref, bc%U_dynH, N2_ref, time_current, grid%height)
    
                call get_eld(eld, TKE_buoyancy, grid%dz, grid%cz)
                call get_eld_ref(eld_ref, N2_ref, Bsurf_ref, time_current, grid%height)
    
    
                call push_value_to_tseries(scm_output%mld, mld)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
                call push_value_to_tseries(scm_output%mld_ref, mld_ref)
    
                call push_value_to_tseries(scm_output%eld, eld)
                call push_value_to_tseries(scm_output%eld_ref, eld_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)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
                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)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !> writing file output
        ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        if (output_mode.eq.1) then
            write(*, *) ' >> writing netcdf output ...'
    
            call write_netcdf(scm_output, grid%z)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        endif
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        if (output_mode.eq.2) then
            write(*, *) ' >> writing ascii output ...'
    
            call write_ascii(scm_output, grid%z)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        endif
    
    
        if (output_mode.eq.3) then
            write(*, *) ' >> writing tecplot output ...'
    
            call write_tecplot(scm_output, grid%z)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        endif
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! ----------------------------------------------------------------------------
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !> model cleanup
        ! ----------------------------------------------------------------------------
    
        !> deallocate state
    
        call deallocate_state_vec
    
    
        !> deallocate scm
        call deallocate_scm_vec
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !> deallocate time-dependent forcing
        call deallocate_fluxes_data
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        !> removing time slice data
    
        call output_cleanup(scm_output)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! > removing grid data
        call deallocate_grid(grid)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! ----------------------------------------------------------------------------