Skip to content
Snippets Groups Projects
obl_main.f90 17.7 KiB
Newer Older
  • Learn to ignore specific revisions
  • #include "obl_def.fi"
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !< @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
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! modules used
    
    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
    
        use obl_diag
    
        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
    
    #ifdef USE_SFX
        use sfx_data, only: meteoDataType, sfxDataType
        use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
                      numericsType_most => numericsType
    
    #endif
    #ifdef USE_CONFIG_PARSER
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
       use config_parser
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !use vertical_mixing, default = off
        !use vermix
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! directives list
    
        type(gridDataType) :: grid
    
    
        ! turbulence closure parameters
    
        type(pphParamType) :: param_pph
        type(pphDynParamType) :: param_pph_dyn
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        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) :: 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] 
    
        !< 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 :: domain_height
        integer :: grid_cz
    
    
        !< output 
        type(oblOutputStx) :: scm_output
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        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
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        integer :: closure_mode  !< closure type:
                                 !1 - pacanowski-philander, 2 - pacanowski-philander+, 
                                 !3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        integer, parameter :: output_mode = 3   ! 1 -- netcdf, 2 -- ascii, 3 -- tecplot
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
        integer, parameter :: obl_setup = 1     ! 1 - Kato-Phillips, 2 - Papa station
    
    
        real, parameter :: Cd0 = 0.001           ! bottom drag coefficient
    
    
        real, parameter :: lw_albedo = 0.03
        real, parameter :: surface_emissivity = 0.97
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        real :: mld !< mixed layer depth, [m]
        real :: lab_mld !< theoretical mixed layer depth, [m]
    
        !< just a temporary to hold value 
        real :: fvalue
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        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
        ! --------------------------------------------------------------------------------
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! screen output parameters
        integer, parameter :: nscreen = 1000
    
        ! file output parameters
        integer, parameter :: noutput = 60
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        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
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
    
    
        !< 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) 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
    
        ! ---------------------------------------------------------------------------- 
        
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
        ! setting grid -- just an example
        !   > [zpos, height, cz, gcz]
    
        call set_uniform_grid(grid, 0.0, domain_height, grid_cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
        ! debug grid print
        call print_grid(grid)
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! allocation
    
        call allocate_state_vec(grid%cz)
    
        ! initialize scm
        call init_scm_vec(grid%cz)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! initialization of main fields
    
        call Theta_init (Theta, Theta_dev, grid%cz, grid%dz)
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !call Theta_init (Theta, nz, dz)
    
        call Salin_init(Salin, Salin_dev, grid%cz, grid%dz)
    
        call get_density(Rho, Theta_dev, Salin_dev, grid%cz)
    
        call U_init(U, grid%cz)
        call V_init(V, grid%cz)
    
        !initialization of TKE & eps in case of k-epsilon closure
        if (closure_mode.eq.3 .or. closure_mode.eq.4) then
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            call TKE_init(TKE, param_k_epsilon, grid%cz)
            call eps_init(EPS, param_k_epsilon, grid%cz, grid%height)
    
        !< 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)
    
            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)
    
        ! ----------------------------------------------------------------------------
    
        !< 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)
        ! ----------------------------------------------------------------------------
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !open (199, file= 'output_Daria/mld.txt', status ='replace')
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        N2_0 = 0.0004
    
    
        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
    
            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) 
    
    
            !< 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
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            bc%u_momentum_fluxH = bc%u_momentum_fluxH / Rho_ref
            bc%v_momentum_fluxH = bc%v_momentum_fluxH / Rho_ref
    
            !< 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)
    
            call get_rho_dynamic(bc%rho_dynH, &
                bc%U_dynH, 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
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            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)
    
            call get_rho_dynamic(bc%rho_dyn0, &
                bc%U_dyn0, bc%heat_flux0, bc%salin_flux0) 
    
            ! ----------------------------------------------------------------------------
    
    
            !< 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)
    
            ! ----------------------------------------------------------------------------
    
            !< 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
            i = i + 1
            time_current = time_current + dt
    
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            !> advance screen output
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            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)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
    
    
                ! *: add finite check
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
                write(*, '(a,g0)') ' mld = ', mld
    
    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
            !> 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)
    
    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
        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
    
        !> deallocate state
    
        call deallocate_state_vec
    
    
        !> 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(sw_flux_surf)
    
        call deallocate_tforcing(lw_in_surf)
        call deallocate_tforcing(lw_net_surf)
    
    
    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)