obl_pph_dyn.f90

module obl_pph_dyn
    !< @brief dynamic pacanowski-philander scheme
    ! --------------------------------------------------------------------------------

    
    ! modules used
    ! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
    use iso_c_binding, only: C_NULL_CHAR
    use config_parser
#endif

    use obl_grid
    use obl_bc
    use obl_turb_common
    use obl_diag
    
    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: init_turbulence_closure
    public :: advance_turbulence_closure
    public :: define_stability_functions

    public :: get_eddy_viscosity, get_eddy_diffusivity
    public :: get_eddy_viscosity_vec, get_eddy_diffusivity_vec

    public :: set_config_param
    ! --------------------------------------------------------------------------------

    !> @brief Pacanowski-Philander dynamic parameters
    type, public :: pphDynParamType
        real :: alpha = 4.0         !< constant, *: default = 5.0
        real :: n = 2.0             !< constant
        real :: Km_b = 0.0001       !< background eddy viscosity, [m**2 / s], *: INMCM: 5 * 10**(-6)  
        real :: Kh_b = 0.00001      !< background eddy diffusivity, [m**2 / s]

        real :: C_l = 4.0 / 9.0     !< length scale tuning constant, [n/d], *: = (4.0 / pi^2) in Obukhov length scale
        
        real :: kappa = 0.4         !< von Karman constant
        real :: PrT0 = 1.0          !< turbulent Prandtl number, *: probably redundant

        real :: c_S2_min = 1e-5     !< min shear frequency, [(1/s)**2]
    end type

    contains


    ! --------------------------------------------------------------------------------
    subroutine define_stability_functions(param, bc, grid)
        !< @brief advance stability functions: N**2, S**2, Ri-gradient
        ! ----------------------------------------------------------------------------
        use obl_state

        type(pphDynParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid
        ! --------------------------------------------------------------------------------

        call get_N2(N2, Rho, bc%Rho_dyn0, bc%Rho_dynH, &
            param%kappa, param%PrT0, grid)
        call get_S2(S2, U, V, bc%U_dyn0, bc%U_dynH, &
            param%kappa, grid)
        call get_Ri_gradient(Ri_grad, N2, S2, param%c_S2_min, grid)
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine init_turbulence_closure(param, bc, grid)
        !< @brief advance turbulence closure
        ! ----------------------------------------------------------------------------
        use obl_state

        type(pphDynParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid

        real :: mld
        ! --------------------------------------------------------------------------------

        call get_mld(mld, N2, grid%dz, grid%cz)

        call get_eddy_viscosity(Km, Ri_grad, &
            bc%U_dynH, mld, param, grid)
        call get_eddy_diffusivity(Kh, Ri_grad, &
            bc%U_dynH, mld, param, grid)
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine advance_turbulence_closure(param, bc, grid, dt)
        !< @brief advance turbulence closure
        ! ----------------------------------------------------------------------------
        use obl_state

        type(pphDynParamType), intent(in) :: param
        type(oblBcType), intent(in) :: bc
        type (gridDataType), intent(in) :: grid
        real, intent(in) :: dt

        real :: mld
        ! --------------------------------------------------------------------------------

        call get_TKE_production(TKE_production, Km, S2, grid)
        call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid)

        call get_mld(mld, N2, grid%dz, grid%cz)

        call get_eddy_viscosity(Km, Ri_grad, &
            bc%U_dynH, mld, param, grid)
        call get_eddy_diffusivity(Kh, Ri_grad, &
            bc%U_dynH, mld, param, grid)
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine get_eddy_viscosity(Km, Ri_grad, U_dynH, mld, param, grid)
        !< @brief calculate eddy viscosity
        ! ----------------------------------------------------------------------------
        type(pphDynParamType), intent(in) :: param
        type (gridDataType), intent(in) :: grid

        real, dimension(grid%cz), intent(out) :: Km     !< eddy viscosity, [m**2 / s]
        real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number, [n/d] 

        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
        real, intent(in) :: mld         !< mixed layer depth, [m]
        ! ----------------------------------------------------------------------------
        
        call get_eddy_viscosity_vec(Km, Ri_grad, U_dynH, mld, param, grid%cz)
    end subroutine

    subroutine get_eddy_viscosity_vec(Km, Ri_grad, U_dynH, mld, param, cz)
        !< @brief calculate eddy viscosity
        ! ----------------------------------------------------------------------------
        type(pphDynParamType), intent(in) :: param

        integer, intent(in) :: cz                   !< grid size 
        real, dimension(cz), intent(out) :: Km      !< eddy viscosity, [m**2 / s]
        real, dimension(cz), intent(in) :: Ri_grad  !< Richardson gradient number, [n/d] 

        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
        real, intent(in) :: mld         !< mixed layer depth, [m]

        real :: Km_0                    !< neutral eddy viscosity, [m**2 / s]
        integer :: k
        ! ----------------------------------------------------------------------------
        
        Km_0 = param%C_l * U_dynH * param%kappa * mld / 2.0

        do k = 1, cz
            Km(k) = Km_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n) + param%Km_b
        end do
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine get_eddy_diffusivity(Kh, Ri_grad, U_dynH, mld, param, grid)
        !< @brief calculate eddy diffusivity
        ! ----------------------------------------------------------------------------
        type(pphDynParamType), intent(in) :: param
        type (gridDataType), intent(in) :: grid

        real, dimension(grid%cz), intent(out) :: Kh     !< eddy diffusivity, [m**2 / s]
        real, dimension(grid%cz), intent(in) :: Ri_grad !< Richardson gradient number, [n/d]

        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
        real, intent(in) :: mld         !< mixed layer depth, [m]
        ! ----------------------------------------------------------------------------

        call get_eddy_diffusivity_vec(Kh, Ri_grad, U_dynH, mld, param, grid%cz)
    end subroutine

    subroutine get_eddy_diffusivity_vec(Kh, Ri_grad, U_dynH, mld, param, cz)
        !< @brief calculate eddy diffusivity
        ! ----------------------------------------------------------------------------
        type(pphDynParamType), intent(in) :: param

        integer, intent(in) :: cz                   !< grid size 
        real, dimension(cz), intent(out) :: Kh      !< eddy diffusivity, [m**2 / s]
        real, dimension(cz), intent(in) :: Ri_grad  !< Richardson gradient number, [n/d]

        real, intent(in) :: U_dynH      !< dynamic velocity, [m/s]
        real, intent(in) :: mld         !< mixed layer depth, [m]

        real :: Kh_0                    !< neutral eddy diffusivity, [m**2 / s]
        integer :: k
        ! ----------------------------------------------------------------------------

        Kh_0 = param%C_l * U_dynH * param%kappa * mld / 2.0

        do k = 1, cz
            Kh(k) = Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b
        end do
    end subroutine

    ! --------------------------------------------------------------------------------
    subroutine set_config_param(param, tag, ierr)
        !< @brief set turbulence closure parameters
        ! ----------------------------------------------------------------------------
        type(pphDynParamType), intent(out) :: param
        integer, intent(out) :: ierr

        character(len = *), intent(in) :: tag

        integer :: status
        ! --------------------------------------------------------------------------------

        ierr = 0        ! = OK

#ifdef USE_CONFIG_PARSER
call c_config_is_varname(trim(tag)//".alpha"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".alpha"//C_NULL_CHAR, param%alpha, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".n"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".n"//C_NULL_CHAR, param%n, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".Km_b"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".Km_b"//C_NULL_CHAR, param%Km_b, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".Kh_b"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".Kh_b"//C_NULL_CHAR, param%Kh_b, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".C_l"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".C_l"//C_NULL_CHAR, param%C_l, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".kappa"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".kappa"//C_NULL_CHAR, param%kappa, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".PrT0"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".PrT0"//C_NULL_CHAR, param%PrT0, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
        call c_config_is_varname(trim(tag)//".c_S2_min"//C_NULL_CHAR, status)
        if (status /= 0) then
            call c_config_get_float(trim(tag)//".c_S2_min"//C_NULL_CHAR, param%c_S2_min, status)
            if (status == 0) then
                ierr = 1        ! signal ERROR
                return
            end if
        endif
#else
        !> *: just skipping setup without configuration file
#endif

    end subroutine

end module