scale_atm_phy_bl_dgm_mynn_lv2.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics / boundary layer turbulence
!!
!! @par 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.
!!
!! @par 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
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_bl_dgm_mynn_lv2
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    eps  => const_eps,     &
    grav => const_grav,    &
    rdry => const_rdry,    &
    cpdry => const_cpdry,  &
    cvdry => const_cvdry,  &
    pres00 => const_pre00, &
    karman  => const_karman, &
    rplanet => const_radius    
  use scale_sparsemat  
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_phy_bl_dgm_mynn_lv2_init
  public :: atm_phy_bl_dgm_mynn_lv2_final
  public :: atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
  !- MYNN closure constants -----------
 
  real(RP), parameter :: A1 = 1.18_rp
  real(RP), parameter :: A2 = 0.665_rp
 
  real(RP), parameter :: B1 = 24.0_rp
  real(RP), parameter :: B2 = 15.0_rp
 
  real(RP), parameter :: C1 = 0.137_rp
  real(RP), parameter :: C2 = 0.75_rp
  real(RP), parameter :: C3 = 0.352_rp
  real(RP), parameter :: C5 = 0.2_rp
 
  real(RP), parameter :: G1 = 0.235_rp
 
  ! Derived constants
  real(RP) :: G2
  real(RP) :: F2
  real(RP) :: Rf2
  real(RP) :: RFc
 
  ! Limiter with mixing length
  real(RP) :: L_INF = 100.0_rp
  real(RP) :: L_MIN = 1.0e-6_rp
 
contains
  !> Initialize a module of MYNN Level 2 PBL turbulence parameterization
!OCL SERIAL
  subroutine atm_phy_bl_dgm_mynn_lv2_init( mesh )
    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
  end subroutine atm_phy_bl_dgm_mynn_lv2_init
 
  !> Finalize a module of MYNN Level 2 PBL turbulence parameterization
!OCL SERIAL
  subroutine atm_phy_bl_dgm_mynn_lv2_final()
    implicit none    
    !--------------------------------------------------------------------------------
    return
  end subroutine atm_phy_bl_dgm_mynn_lv2_final
 
!OCL SERIAL
  subroutine atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef( &
    Nu, Kh, TKE,                                             & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
    rtot, pres, pt,                                          & ! (in)
    dz, lift, lmesh, elem, is_bound  )                         ! (in)
    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
  end subroutine atm_phy_bl_dgm_mynn_lv2_cal_vviscdiffcoef
 
!-- private --------------------------------------------------------
 
!OCL SERIAL  
  subroutine 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)
    nz, fscale, vmapm, vmapp, lmesh, elem, is_bound                                       ) ! (in)
 
    implicit none
 
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem  
    real(rp), intent(out) ::  del_flux_rho(elem%nfptot*lmesh%ne)
    real(rp), intent(out) ::  del_flux_mom(elem%nfptot*lmesh%ne,3)
    real(rp), intent(out) ::  del_flux_rhot(elem%nfptot*lmesh%ne)
    real(rp), intent(out) ::  del_flux_rtot(elem%nfptot*lmesh%ne)
    real(rp), intent(in) ::  ddens_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  momx_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  momy_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  momz_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  drhot_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  rtot(elem%np*lmesh%nea)
    real(rp), intent(in) ::  dens_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pres_hyd(elem%np*lmesh%nea)
    real(rp), intent(in) ::  pt_(elem%np*lmesh%nea)
    real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
    real(rp), intent(in) :: fscale(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
    logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
    
    integer :: i, ip, im
    real(rp) :: densm, densp
    real(rp) :: facz
    !------------------------------------------------------------------------
    
    !$omp parallel do private ( iM, iP, densM, densP, facz )
    do i=1, elem%NfpTot * lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
 
      densm = ddens_(im) + dens_hyd(im)
      densp = ddens_(ip) + dens_hyd(ip)
 
      if ( is_bound(i) ) then
        facz = 0.5_rp * fscale(i) * nz(i)
      else
        ! facz = ( 1.0_RP - sign(1.0_RP,nz(i)) ) * Fscale(i) * nz(i)
        facz = 0.5_rp * fscale(i) * nz(i)
      end if
 
      del_flux_rho(i) = facz * ( densp - densm )
 
      del_flux_mom(i,1) = facz * ( momx_(ip) - momx_(im) )
      del_flux_mom(i,2) = facz * ( momy_(ip) - momy_(im) )
      del_flux_mom(i,3) = facz * ( momz_(ip) - momz_(im) )
 
      del_flux_rhot(i) = facz * ( densp * pt_(ip) - densm * pt_(im) )
      del_flux_rtot(i) = facz * ( rtot(ip) - rtot(im) )
    end do
 
    return
  end subroutine cal_del_flux_grad
end module scale_atm_phy_bl_dgm_mynn_lv2