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

            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)
            !write(*,*) 'here_coeffs', k, turb%l_turb_m(k),  fnriuv, fnritq
        end do
        ! for compliance with old version
        write(*,*) 'here'
        bl%vdcuv(grid%kmax) = 0.0
        bl%vdcuv(grid%kmax) = 0.0
        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)
            write(*,*)  bl%vdctq(k), bl%vdcuv(k), turb%km(k)
        end do
        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
        real l_neutral, lm, lh, hpbla_diag, z_surf
        integer kpbl

        integer k, kmax
        kmax = grid%kmax
        z_surf = grid%z_edge(kmax)
        if (hbl_option ==1) then
            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

        do k = kmax, 1, -1
            l_neutral = fluid%kappa * grid%z_cell(k)
            if ( k < bl%kpbl+1) then
                lm = 0.0
                lh = 0.0
            else
                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
            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