Skip to content
Snippets Groups Projects
obl_pph.f90 10.6 KiB
Newer Older
  • Learn to ignore specific revisions
  • #include "obl_def.fi"
    
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        !< @brief standard pacanowski-philander scheme
    
        ! --------------------------------------------------------------------------------
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! modules used
    
        ! --------------------------------------------------------------------------------
    
    #ifdef USE_CONFIG_PARSER
        use iso_c_binding, only: C_NULL_CHAR
        use config_parser
    #endif
    
    
        use obl_grid
    
        use obl_turb_common
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        ! directives list
    
        ! --------------------------------------------------------------------------------
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        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 parameters
    
        type, public :: pphParamType
    
            real :: Km_0 = 2.0 * 0.01   !< neutral eddy viscosity, [m**2 / s] *: INMCM: 7.0 * 0.01
            real :: Kh_0 = 2.0 * 0.01   !< neutral eddy diffusivity, [m**2 / s], *: INMCM: 5.0 * 0.01
            real :: alpha = 5.0         !< constant
            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 :: kappa = 0.4         !< von Karman constant, [n/d]
            real :: PrT0 = 1.0          !< neutral Prandtl number, [n/d]
    
            real :: c_S2_min = 1e-4     !< min shear frequency, [(1/s)**2]
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
        contains
    
        ! --------------------------------------------------------------------------------
    
        subroutine define_stability_functions(param, bc, grid)
    
            !< @brief advance stability functions: N**2, S**2, Ri-gradient
            ! ----------------------------------------------------------------------------
            use obl_state
    
            type(pphParamType), 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)
    
        ! --------------------------------------------------------------------------------
        subroutine init_turbulence_closure(param, bc, grid)
            !< @brief advance turbulence closure
            ! ----------------------------------------------------------------------------
            use obl_state
    
            type(pphParamType), intent(in) :: param
            type(oblBcType), intent(in) :: bc
            type (gridDataType), intent(in) :: grid
            ! --------------------------------------------------------------------------------
    
            call get_eddy_viscosity(Km, Ri_grad, param, grid)
            call get_eddy_diffusivity(Kh, Ri_grad, param, grid)
        end subroutine
    
    
        ! --------------------------------------------------------------------------------
    
        subroutine advance_turbulence_closure(param, bc, grid, dt)
    
            !< @brief advance turbulence closure
            ! ----------------------------------------------------------------------------
            use obl_state
    
            type(pphParamType), intent(in) :: param
            type(oblBcType), intent(in) :: bc
            type (gridDataType), intent(in) :: grid
            real, intent(in) :: dt
            ! --------------------------------------------------------------------------------
    
    
            call get_TKE_production(TKE_production, Km, S2, grid)
            call get_TKE_buoyancy(TKE_buoyancy, Kh, N2, grid)
    
            call get_eddy_viscosity(Km, Ri_grad, param, grid)
            call get_eddy_diffusivity(Kh, Ri_grad, param, grid)
        end subroutine
    
        ! --------------------------------------------------------------------------------
        subroutine get_eddy_viscosity(Km, Ri_grad, param, grid)
            !< @brief calculate eddy viscosity on grid
            ! ----------------------------------------------------------------------------
            type(pphParamType), intent(in) :: param
            type (gridDataType), intent(in) :: grid
    
            real, dimension(grid%cz), intent(in) :: Ri_grad      !< Richardson gradient number 
            real, dimension(grid%cz), intent(out) :: Km          !< eddy viscosity, [m**2 / s]
            ! --------------------------------------------------------------------------------
    
            call get_eddy_viscosity_vec(Km, Ri_grad, param, grid%cz)
        end subroutine
    
    
        subroutine get_eddy_viscosity_vec(Km, Ri_grad, param, n)
    
            !< @brief calculate eddy viscosity
    
            ! ----------------------------------------------------------------------------
    
            type(pphParamType), intent(in) :: param
    
            integer, intent(in) :: n                       !< vector size
            real, dimension(n), intent(in) :: Ri_grad      !< Richardson gradient number 
            real, dimension(n), intent(out) :: Km          !< eddy viscosity, [m**2 / s]
    
            ! --------------------------------------------------------------------------------
    
            do k = 1, n
    
                Km(k) = param%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, param, grid)
            !< @brief calculate eddy diffusivity on grid
    
            ! ----------------------------------------------------------------------------
    
            type(pphParamType), intent(in) :: param
    
            type (gridDataType), intent(in) :: grid
    
            real, dimension(grid%cz), intent(in) :: Ri_grad      !< Richardson gradient number 
    
            real, dimension(grid%cz), intent(out) :: Kh          !< eddy diffusivity, [m**2 / s]
            ! --------------------------------------------------------------------------------
            
            call get_eddy_diffusivity_vec(Kh, Ri_grad, param, grid%cz)
    
        end subroutine
    
        subroutine get_eddy_diffusivity_vec(Kh, Ri_grad, param, n)
            !< @brief calculate eddy diffusivity
            ! ----------------------------------------------------------------------------
    
            type(pphParamType), intent(in) :: param
    
            integer, intent(in) :: n                       !< vector size
            real, dimension(n), intent(in) :: Ri_grad      !< Richardson gradient number 
            real, dimension(n), intent(out) :: Kh          !< eddy diffusivity, [m**2 / s]
    
            ! --------------------------------------------------------------------------------
    
            do k = 1, n
    
                Kh(k) = param%Kh_0 / (1.0 + param%alpha * Ri_grad(k))**(param%n + 1.0) + param%Kh_b
    
    Daria Gladskikh's avatar
    Daria Gladskikh committed
            end do
        end subroutine
    
    
        ! --------------------------------------------------------------------------------
        subroutine set_config_param(param, tag, ierr)
            !< @brief set turbulence closure parameters
            ! ----------------------------------------------------------------------------
            type(pphParamType), 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)//".Km_0"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_float(trim(tag)//".Km_0"//C_NULL_CHAR, param%Km_0, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            endif
            call c_config_is_varname(trim(tag)//".Kh_0"//C_NULL_CHAR, status)
            if (status /= 0) then
                call c_config_get_float(trim(tag)//".Kh_0"//C_NULL_CHAR, param%Kh_0, status)
                if (status == 0) then
                    ierr = 1        ! signal ERROR
                    return
                end if
            endif
            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)//".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