Skip to content
Snippets Groups Projects
pbl_fo_turb.f90 5.34 KiB
Newer Older
  • Learn to ignore specific revisions
  • ! Created by Andrey Debolskiy on 29.11.2024.
    
    module pbl_fo_turb
        implicit none
    
        !constants to for stability functions
        real, parameter:: yb = 5.e0
        real, parameter:: yc = 5.e0
        real, parameter:: yd = 5.e0
        real:: ye = yb*yc*sqrt(1.e0/3.e0)
        real:: ricr = 2.e0/(3.e0*yd)
        real, parameter:: vlinfm = 160.0e0
        real:: vlinfh = vlinfm * sqrt(3.e0*yd/2.e0)
        real, parameter::  aka = 0.4e0
        !real, parameter:: r = 287.05
        !real, parameter:: cp = 1.005
        !    security parameters
        !    du2min is a minimal value to calculate wind shear
        !   cormin is a minimal value of coriolis parameter to calculate the
        !   pbl extention
        real, parameter::  du2min = 1.e-02
        real, parameter::  cormin = 5.e-05
    
        public get_fo_coeffs
    
        public get_fo_coeffs_louis
    
        public get_turb_length
    
    !    public :: default_stab_functions
    !    public :: esau_stab_functions
    !    private:: get_ph_efb
    !    private:: q_root
    
        contains
    
    
    !    subroutine get_fo_diff(turb, bl, fluid, grid)
    !        use pbl_turb_data, only: turbBLDataType
    !        use scm_state_data, only: stateBLDataType
    !        use pbl_grid, only: pblgridDataType
    !        use phys_fluid, only: fluidParamsDataType
    !    end subroutine get_fo_diff
        subroutine get_fo_coeffs(turb, bl, fluid,  grid)
    
            use pbl_turb_data, only: turbBLDataType
            use scm_state_data, only: stateBLDataType
            use pbl_grid, only: pblgridDataType
            use phys_fluid, only: fluidParamsDataType
    
            implicit none
            type(fluidParamsDataType), intent(in):: fluid
            type(stateBLDataType), intent(inout):: bl
            type(pblgridDataType), intent(in)::grid
            type(turbBLDataType), intent(inout):: turb
    
            if (turb%sf_type == 1) then !> Louis stab functions
                call get_fo_coeffs_louis(turb, bl, fluid,  grid)
            elseif (turb%sf_type == 2) then !> Esau stab functions
                ! call get_fo_coeffs_esau(turb, bl, fluid,  grid)
            elseif  (turb%sf_type == 3)    then   !> EFB stab functions
                ! call get_fo_coeffs_efb(turb, bl, fluid,  grid)
            else
                write(*,*) 'FAILURE: stab functions type in FOM not set!'
            end if
        end subroutine get_fo_coeffs
    
        subroutine get_fo_coeffs_louis(turb, bl, fluid,  grid)
    
            use pbl_turb_data, only: turbBLDataType
            use scm_state_data, only: stateBLDataType
            use pbl_grid, only: pblgridDataType
            use phys_fluid, only: fluidParamsDataType
    
            use fo_stability_functions
            implicit none
            type(fluidParamsDataType), intent(in):: fluid
            type(stateBLDataType), intent(inout):: bl
            type(pblgridDataType), intent(in)::grid
            type(turbBLDataType), intent(inout):: turb
    
            real fnriuv, fnritq,rho
            real ocean_flag
    
            integer k, kmax
    
            ocean_flag = real(bl%land_type)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            write(*,*) 'here'
            do k = grid%kmax-1, 1, -1
    
                call louis_stab_functions(fnriuv, fnritq, turb%rig(k), turb%l_turb_m(k), &
                        & turb%l_turb_m(k), grid%z_cell(k),ocean_flag)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
    
    
                turb%km(k) = fnriuv * turb%l_turb_m(k) * turb%l_turb_m(k) * turb%s2(k)
                turb%kh(k) = fnritq * turb%l_turb_h(k) * turb%l_turb_h(k) * turb%s2(k)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                !write(*,*) 'here_coeffs', k, turb%l_turb_m(k),  fnriuv, fnritq
    
            end do
            ! for compliance with old version
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            write(*,*) 'here'
    
            bl%vdcuv(grid%kmax) = 0.0
            bl%vdcuv(grid%kmax) = 0.0
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            write(*,*) 'here3.25'
            do k = grid%kmax-1,1, -1
    
                bl%vdcuv(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%km(k)
                bl%vdctq(k) = 0.5 * (bl%rho(k) + bl%rho(k+1)) * turb%kh(k)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                write(*,*)  bl%vdctq(k), bl%vdcuv(k), turb%km(k)
    
            end do
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            write(*,*) 'here3.5'
    
    
        end subroutine get_fo_coeffs_louis
    
    
        subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
            use pbl_turb_data, only: turbBLDataType
            use scm_state_data, only: stateBLDataType
            use pbl_grid, only: pblgridDataType
            use phys_fluid, only: fluidParamsDataType
            use diag_pbl
            implicit none
            type(fluidParamsDataType), intent(in):: fluid
    
            type(stateBLDataType), intent(inout):: bl
    
            type(pblgridDataType), intent(in)::grid
    
            type(turbBLDataType), intent(inout):: turb
            integer hbl_option
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            real l_neutral, lm, lh, hpbla_diag, z_surf
            integer kpbl
    
            integer k, kmax
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            kmax = grid%kmax
            z_surf = grid%z_edge(kmax)
    
            if (hbl_option ==1) then
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                call get_hpbl(hpbla_diag, kpbl, bl%theta_v, grid%z_cell, &
                        z_surf, grid%kmax , bl%cor, turb%ustr)
                bl%hpbla_diag = hpbla_diag
                bl%kpbl = kpbl
    
            else
                write(*,*) 'wrong hbl diagnostics option'
                return
            end if
    
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            do k = kmax, 1, -1
    
                l_neutral = fluid%kappa * grid%z_cell(k)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                if ( k < bl%kpbl+1) then
    
                    lm = 0.0
                    lh = 0.0
                else
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                    if (turb%rig(k) > 0.0) then
    
                        lm = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
                        lh = (1.0 - grid%z_cell(k) / bl%hpbla_diag)
                    else
                        lm = 1.0
                        lh = 1.0
                    endif
                endif
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                turb%l_turb_m(k) = l_neutral * lm
                turb%l_turb_h (k) = l_neutral * lh
    
            end do
    
    
        end subroutine get_turb_length
    
    
    end module pbl_fo_turb