module FElib / Atmosphere / Physics / boundary layer turbulence More...

Functions/Subroutines
subroutine, public	atm_phy_bl_dgm_mynn_lv2_init (mesh)
	Initialize a module of MYNN Level 2 PBL turbulence parameterization.
subroutine, public	atm_phy_bl_dgm_mynn_lv2_final ()
	Finalize a module of MYNN Level 2 PBL turbulence parameterization.
subroutine, public	atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef (nu, kh, tke, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, rtot, pres, pt, dz, lift, lmesh, elem, is_bound)

Detailed Description

module FElib / Atmosphere / Physics / boundary layer turbulence

Description: A module to provide routines to calculate turbulent diffusion coefficients based on the Mellor-Yamada-Nakanishi-Niino (MYNN) Level 2 closure.

This implementation uses the Level 2 diagnostic closure only. Prognostic turbulent kinetic energy, nonlocal mixing, and the MY Level 2.5 equations are not included.

Reference: Mellor, G. L. and T. Yamada, 1982: Development of a turbulence closure model for geophysical fluid problems. Reviews of Geophysics and Space Physics, 20, 851-875.

Nakanishi, M. and H. Niino, 2009: Development of an improved turbulence closure model for the atmospheric boundary layer. Journal of the Meteorological Society of Japan, 87, 895-912.

Blackadar, A. K., 1962: The vertical distribution of wind and turbulent exchange in a neutral atmosphere. Journal of Geophysical Research, 67, 3095-3102.

Author: Yuta Kawai, Team SCALE

NAMELIST

PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2

name type default value comment

L_INF real(RP) 100.0_RP

L_MIN real(RP) 1.0E-6_RP

History Output: No history output

Function/Subroutine Documentation

◆ atm_phy_bl_dgm_mynn_lv2_init()

subroutine, public scale_atm_phy_bl_dgm_mynn_lv2::atm_phy_bl_dgm_mynn_lv2_init ( class(meshbase3d), intent(in) mesh )

Initialize a module of MYNN Level 2 PBL turbulence parameterization.

Definition at line 104 of file scale_atm_phy_bl_dgm_mynn_lv2.F90.

    implicit none    
    class(MeshBase3D), intent(in) :: mesh
 
    namelist / param_atmos_phy_bl_dgm_mynn_lv2 / &
      l_inf, l_min
    
    integer :: ierr    
    !--------------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Setup'
    log_info("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'MYNN Level 2 PBL turbulence parameterization'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_bl_dgm_mynn_lv2,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_bl_dgm_mynn_lv2)
 
    !-
 
    g2 = ( 2.0_rp * a1 * ( 3.0_rp - 2.0_rp * c2 )  &
         + b2 * ( 1.0_rp - c3 )                    &
         ) / b1
 
    f2 = b1 * ( g1 + g2 ) &
        - 3.0_rp * a1 * ( 1.0_rp - c2 )
 
    rf2 = b1 * g1 / f2
 
    rfc = g1 / ( g1 + g2 )
    return

◆ atm_phy_bl_dgm_mynn_lv2_final()

subroutine, public scale_atm_phy_bl_dgm_mynn_lv2::atm_phy_bl_dgm_mynn_lv2_final

Finalize a module of MYNN Level 2 PBL turbulence parameterization.

Definition at line 146 of file scale_atm_phy_bl_dgm_mynn_lv2.F90.

    implicit none    
    !--------------------------------------------------------------------------------
    return

◆ atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef()

subroutine, public scale_atm_phy_bl_dgm_mynn_lv2::atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	nu,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	kh,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	tke,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens_,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pt,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	lift,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		logical, dimension(elem%nfptot,lmesh%ne), intent(in)	is_bound )

Parameters

[out]	nu	Vertical eddy viscosity
[out]	kh	Vertical eddy diffusivity
[out]	tke	Turbulent kinetic energy
[in]	ddens_	Density perturbation
[in]	momx_	Momentum in x1 direction
[in]	momy_	Momentum in x2 direction
[in]	momz_	Momentum in x3 direction
[in]	drhot_	Density x potential temperature perturbation
[in]	dens_hyd	Reference density in hydrostatic balance
[in]	pres_hyd	Reference pressure in hydrostatic balance
[in]	rtot	Gas constant
[in]	pres	Pressure
[in]	pt	Potential temperature
[in]	is_bound	Flag whether nodes are located at domain boundaries

Definition at line 153 of file scale_atm_phy_bl_dgm_mynn_lv2.F90.

    implicit none
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    real(RP), intent(out) :: Nu(elem%Np,lmesh%NeA)       !< Vertical eddy viscosity
    real(RP), intent(out) :: Kh(elem%Np,lmesh%NeA)       !< Vertical eddy diffusivity
    real(RP), intent(out) :: TKE(elem%Np,lmesh%NeA)      !< Turbulent kinetic energy
    real(RP), intent(in)  :: DDENS_(elem%Np,lmesh%NeA)   !< Density perturbation
    real(RP), intent(in)  :: MOMX_ (elem%Np,lmesh%NeA)   !< Momentum in x1 direction
    real(RP), intent(in)  :: MOMY_ (elem%Np,lmesh%NeA)   !< Momentum in x2 direction
    real(RP), intent(in)  :: MOMZ_ (elem%Np,lmesh%NeA)   !< Momentum in x3 direction
    real(RP), intent(in)  :: DRHOT_(elem%Np,lmesh%NeA)   !< Density x potential temperature perturbation
    real(RP), intent(in)  :: DENS_hyd(elem%Np,lmesh%NeA) !< Reference density in hydrostatic balance
    real(RP), intent(in)  :: PRES_hyd(elem%Np,lmesh%NeA) !< Reference pressure in hydrostatic balance
    real(RP), intent(in)  :: Rtot(elem%Np,lmesh%NeA)     !< Gas constant
    real(RP), intent(in)  :: PRES(elem%Np,lmesh%NeA)     !< Pressure
    real(RP), intent(in)  :: PT(elem%Np,lmesh%NeA)       !< Potential temperature
    type(SparseMat), intent(in) :: Dz, Lift
    logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) !< Flag whether nodes are located at domain boundaries
 
    real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: DENS(elem%Np), RDENS(elem%Np), RHOT(elem%Np), Q(elem%Np)
    real(RP) :: DdensDz(elem%Np), DVelDz(elem%Np,3), DptDz(elem%Np), DrtotDz(elem%Np)
 
    real(RP) :: del_flux_rho (elem%NfpTot,lmesh%Ne)
    real(RP) :: del_flux_mom (elem%NfpTot,lmesh%Ne,3)
    real(RP) :: del_flux_rhot(elem%NfpTot,lmesh%Ne)
    real(RP) :: del_flux_Rtot(elem%NfpTot,lmesh%Ne)
 
    real(RP) :: S2(elem%Np) ! Square of vertical shear of horizontal wind
    real(RP) :: N2          ! Square of Brunt-Vaisala frequency
    real(RP) :: Ri          ! Richardson number
    real(RP) :: Rf(elem%Np) ! Flux Richardson number
 
    real(RP) :: A2_loc
    real(RP) :: F1
    real(RP) :: Rf1
    real(RP) :: AF12
    real(RP) :: discriminant
    real(RP) :: denom_h
    real(RP) :: denom_m
    real(RP) :: S_M(elem%Np), S_H(elem%Np)     ! Stability functions for momentum and heat
 
    real(RP) :: mixlen(elem%Np)                    ! Mixing length
    real(RP) :: zsfc(elem%Nnode_h1D**2,lmesh%Ne2D) ! Surface height
    real(RP) :: dz1(elem%Nnode_h1D**2,lmesh%Ne2D)  ! Vertical grid spacing at the first layer
    integer :: HSlice_b(elem%Nnode_h1D**2)
    integer :: HSlice_t(elem%Nnode_h1D**2)
  
    real(RP) :: kz
 
    integer :: ke, ke2D
    integer :: p, ph, pz
    class(LocalMesh2D), pointer :: lmesh2D
    class(ElementBase2D), pointer :: elem2D
 
    real(RP), parameter :: EPS_RF   = 1.0e-10_rp
    real(RP), parameter :: EPS_DISC = 1.0e-14_rp    
    !--------------------------------------------------------------------------------
 
    ! For standard case of NN2009, these parameter are independent of Ri.
    ! If we will implement the K2010 correction in the future, 
    ! these parameter should be calculated in the ke loop because they depend on Ri (i.e., A2_loc = A_2/(1+Ri)).
    a2_loc = a2
 
    f1 = b1 * ( g1 - c1 )                                    &
       + 2.0_rp * a1 * ( 3.0_rp - 2.0_rp * c2 )              &
       + 3.0_rp * a2_loc * ( 1.0_rp - c3 ) * ( 1.0_rp - c5 )
    
    rf1 = b1 * ( g1  - c1 ) / f1
 
    af12 = a1 * f1 / ( a2_loc * f2 )
 
    lmesh2d => lmesh%lcmesh2D
    elem2d  => lmesh2d%refElem2D
 
    hslice_b(:) = elem%Hslice(:,1)
    hslice_t(:) = elem%Hslice(:,elem%Nnode_v)
 
    !---------------------------------
 
    call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, del_flux_rtot, & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, rtot, dens_hyd, pres_hyd, pt,              & ! (in)
      lmesh%normal_fn(:,:,3), lmesh%Fscale, lmesh%vmapM, lmesh%vmapP,                 & ! (in)
      lmesh, elem, is_bound )                                                           ! (in)
 
    !$omp parallel &
    !$omp private( ke2D, Fz, LiftDelFlx, DENS, RDENS, RHOT, Q, DdensDz, DVelDz, DptDz, DrtotDz, &
    !$omp N2, S2, Ri, Rf, discriminant, denom_m, denom_h, S_M, S_H, mixlen, kz                  )
 
    !$omp do
    do ke2d=lmesh2d%NeS, lmesh2d%NeE
      zsfc(:,ke2d) = lmesh%zlev(hslice_b,ke2d)
      dz1(:,ke2d)  =  ( lmesh%zlev(hslice_t,ke2d) - lmesh%zlev(hslice_b,ke2d) ) / real(elem%Nnode_v,kind=rp) * 0.5_rp      
    end do
    !$omp end do
    !$omp do
    do ke=lmesh%NeS, lmesh%NeE
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
      rdens(:) = 1.0_rp / dens(:)
      rhot(:) = dens(:) * pt(:,ke)
 
      ! gradient of density
      call sparsemat_matmul( dz, dens, fz )
      call sparsemat_matmul( lift, del_flux_rho(:,ke), liftdelflx )
      ddensdz(:) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
      ! gradient of u
      q(:) = momx_(:,ke) * rdens(:) 
      call sparsemat_matmul( dz, momx_(:,ke), fz )
      call sparsemat_matmul( lift, del_flux_mom(:,ke,1), liftdelflx )
      dveldz(:,1) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdz(:) ) * rdens(:)
 
      ! gradient of v
      q(:) = momy_(:,ke) * rdens(:) 
      call sparsemat_matmul( dz, momy_(:,ke), fz )
      call sparsemat_matmul( lift, del_flux_mom(:,ke,2), liftdelflx )
      dveldz(:,2) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdz(:) ) * rdens(:)
    
      ! gradient of pt
      q(:) = rhot(:) * rdens(:)
      call sparsemat_matmul( dz, rhot(:) , fz )
      call sparsemat_matmul( lift, del_flux_rhot(:,ke), liftdelflx )
      dptdz(:) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdz(:) ) * rdens(:)
 
      ! gradient of Rtot
      call sparsemat_matmul( dz, rtot(:,ke), fz )
      call sparsemat_matmul( lift, del_flux_rtot(:,ke), liftdelflx )
      drtotdz(:) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
 
     ! Calculate flux Richardson number: Rf
      do p=1, elem%Np
        n2 = grav * ( dptdz(p) / pt(p,ke) + drtotdz(p) / rtot(p,ke) )
 
        s2(p) = dveldz(p,1)**2 + dveldz(p,2)**2
        ri = n2 / max(s2(p), eps)
 
        discriminant = ri * ri &
                     + 2.0_rp * af12 * ( rf1 - 2.0_rp * rf2 ) * ri &
                     + ( af12 * rf1 )**2
        discriminant = max( discriminant, 0.0_rp )
 
        rf(p) = 0.5_rp / af12 * ( ri + af12 * rf1 - sqrt(discriminant) )
        rf(p) = min( rf(p), rfc - eps_rf )
      end do
 
      ! Calculate stability function: S_M and S_H
      do p=1, elem%Np
        denom_h = max( 1.0_rp - rf(p), eps_rf )
        s_h(p) = 3.0_rp * a2_loc * ( g1 + g2 ) * ( rfc - rf(p) ) / denom_h
        
        denom_m = rf2 - rf(p)
        if ( abs(denom_m) < eps_rf ) then
            denom_m = sign(eps_rf, denom_m)
        end if
        s_m(p) = s_h(p) * af12 * ( rf1 - rf(p) ) / denom_m
 
        s_m(p) = max( s_m(p), 0.0_rp )
        s_h(p) = max( s_h(p), 0.0_rp )
      end do
 
      ! Calculate mixing length
 
      ke2d = lmesh%EMap3Dto2D(ke)
      do pz=1, elem%Nnode_v
      do ph=1, elem%Nnode_h1D**2
        p = ph + (pz-1)*elem%Nnode_h1D**2
        kz = karman * max( lmesh%zlev(p,ke) - zsfc(ph,ke2d), dz1(ph,ke2d) )
 
        mixlen(p) = kz * l_inf  / max( kz + l_inf, l_min )
      end do
      end do
 
      ! Calculate eddy viscosity and diffusivity
      do p=1, elem%Np
        ! q^2
        q(p) = b1 * mixlen(p)**2 &
             * s_m(p) * max( 1.0_rp - rf(p), 0.0_rp ) * s2(p)
        tke(p,ke) = 0.5_rp * max( q(p), 0.0_rp )
 
        q(p) = sqrt( 2.0_rp * tke(p,ke) )
        kh(p,ke) = mixlen(p) * q(p) * s_h(p)
        nu(p,ke) = mixlen(p) * q(p) * s_m(p)
      end do
    end do
    !$omp end do
    !$omp end parallel
 
    return

name	type	default value	comment
L_INF	real(RP)	100.0_RP
L_MIN	real(RP)	1.0E-6_RP

Functions/Subroutines

Detailed Description

Function/Subroutine Documentation

◆ atm_phy_bl_dgm_mynn_lv2_init()

◆ atm_phy_bl_dgm_mynn_lv2_final()

◆ atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef()