scale_atmos_phy_bl_mynn_jmapplib.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_bl_mynn_jmapplib
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
 
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_bl_mynn_jmapplib_setup
  public :: atmos_phy_bl_mynn_jmapplib_finalize
  public :: atmos_phy_bl_mynn_jmapplib_tendency
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
#ifdef JMAPPLIB
  integer,                public, parameter :: ATMOS_PHY_BL_MYNN_JMAPPLIB_NTRACER = 4
  character(len=H_SHORT), public :: ATMOS_PHY_BL_MYNN_JMAPPLIB_NAME(4) = &
       (/ 'TKE_MYNN', 'TSQ_MYNN', 'QSQ_MYNN', 'COV_MYNN' /)
  character(len=H_LONG),  public :: ATMOS_PHY_BL_MYNN_JMAPPLIB_DESC(4) = &
       (/ 'doubled turbulent kinetic energy (MYNN)                                                 ', &
          'sub-grid variance of liquid water potential temperature (MYNN)                          ', &
          'sub-grid variance of total water content (MYNN)                                         ', &
          'sub-grid covariance of liquid water potential temperature and total water content (MYNN)' /)
  character(len=H_SHORT), public :: ATMOS_PHY_BL_MYNN_JMAPPLIB_UNITS(4) = &
       (/ 'm2/s2  ', 'K2     ', 'kg2/kg2', 'K kg   ' /)
#else
  integer,                public, parameter :: atmos_phy_bl_mynn_jmapplib_ntracer = 0
  character(len=H_SHORT), public :: atmos_phy_bl_mynn_jmapplib_name(1) = (/ '' /)
  character(len=H_LONG),  public :: atmos_phy_bl_mynn_jmapplib_desc(1) = (/ '' /)
  character(len=H_SHORT), public :: atmos_phy_bl_mynn_jmapplib_units(1) = (/ '' /)
#endif
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
#ifdef JMAPPLIB
  integer,  private, parameter :: i_tke = 1
  integer,  private, parameter :: i_tsq = 2
  integer,  private, parameter :: i_qsq = 3
  integer,  private, parameter :: i_cov = 4
 
  integer,  private :: ke_pbl
 
  character(len=3), private :: atmos_phy_bl_mynn_jmapplib_level           = "3"
  logical,          private :: atmos_phy_bl_mynn_jmapplib_shcu_buoy_flag  = .true.
  real(rp),         private :: atmos_phy_bl_mynn_jmapplib_pbl_max         = 1.e+99_rp
  real(rp),         private :: atmos_phy_bl_mynn_jmapplib_shcu_max        = 1.e+99_rp
  real(rp),         private :: atmos_phy_bl_mynn_jmapplib_sgm_min_fct     = 0.09_rp
 
  namelist / param_atmos_phy_bl_mynn_jmapplib / &
       atmos_phy_bl_mynn_jmapplib_level, &
       atmos_phy_bl_mynn_jmapplib_shcu_buoy_flag, &
       atmos_phy_bl_mynn_jmapplib_pbl_max, &
       atmos_phy_bl_mynn_jmapplib_shcu_max, &
       atmos_phy_bl_mynn_jmapplib_sgm_min_fct
#endif
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_bl_mynn_jmapplib_setup( &
       KA, KS, KE, &
       CZ, &
       dt, PBL_MAX, SHCU_MAX )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       pre00 => const_pre00
#ifdef JMAPPLIB
    use pbl_const, only: &
       pbl_const_ini
    use pbl_parm, only: &
       pbl_parm_ini
    use pbl_grid, only: &
       pbl_grid_ini
    use pbl_mym_option, only: &
       pbl_mym_option_ini
    use pbl_mym_option_symbol, only: &
       mymodel3, &
       mymodel25
    use pbl_mym_parm, only: &
       pbl_mym_parm_ini
    use pbl_mym_const, only: &
       pbl_mym_const_ini
#endif
    implicit none
 
    integer,  intent(in) :: ka, ks, ke
 
    real(rp), intent(in) :: cz(ka)
 
    real(dp), intent(in), optional :: dt
    real(rp), intent(in), optional :: pbl_max
    real(rp), intent(in), optional :: shcu_max
 
    real(rp) :: e_emit
 
    integer :: ke_shcu
    integer :: nz_pbl, nz_shcu
 
    integer :: ierr
    integer :: k
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_BL_MYNN_JMAPPLIB_setup",*) 'Setup'
    log_info("ATMOS_PHY_BL_MYNN_JMAPPLIB_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme implemented in the JMA Physics Process Library'
 
#ifdef JMAPPLIB
    if ( present(pbl_max) ) atmos_phy_bl_mynn_jmapplib_pbl_max = pbl_max
    if ( present(shcu_max) ) atmos_phy_bl_mynn_jmapplib_shcu_max = shcu_max
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_bl_mynn_jmapplib,iostat=ierr)
    if( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_BL_MYNN_JMAPPLIB_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN_JMAPPLIB. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_bl_mynn_jmapplib)
 
 
    ke_pbl  = ks+1
    ke_shcu = ks+1
    do k = ks+2, ke-1
       if ( atmos_phy_bl_mynn_jmapplib_pbl_max >= cz(k) ) then
          ke_pbl = k
       end if
       if ( atmos_phy_bl_mynn_jmapplib_shcu_max >= cz(k) ) then
          ke_shcu = k
       end if
    end do
    nz_pbl  = ke_pbl - ks + 1
    nz_shcu = min(ke_shcu - ks + 1, nz_pbl)
 
    e_emit = 1.0_rp
    call pbl_const_ini(pref_in = pre00, timestep_in = real(dt,rp), e_emit_in = e_emit)
    call pbl_parm_ini
    call pbl_grid_ini(nz_pbl, shcu_levels_in = nz_shcu)
 
    select case ( atmos_phy_bl_mynn_jmapplib_level )
    case ( "3" )
       call pbl_mym_option_ini(levflag_in = mymodel3, l_shcu_buoy_in = atmos_phy_bl_mynn_jmapplib_shcu_buoy_flag)
 
    case ( "2.5" )
       call pbl_mym_option_ini(levflag_in = mymodel25, l_shcu_buoy_in = atmos_phy_bl_mynn_jmapplib_shcu_buoy_flag)
    case default
       log_error("ATMOS_PHY_BL_MYNN_JMAPPLIB_setup",*) 'only level 2.5 and 3 are supported at this moment'
       call prc_abort
    end select
 
    call pbl_mym_parm_ini(my_sgm_min_fct_in = atmos_phy_bl_mynn_jmapplib_sgm_min_fct)
    call pbl_mym_const_ini
 
#else
 
    log_error("ATMOS_PHY_BL_MYNN_JMAPPLIB_setup",*) 'To use "MYNN-JMAPPLIB", compile SCALE with "SCALE_ENABLE_JMAPPLIB=T" option.'
    call prc_abort
 
#endif
 
    return
  end subroutine atmos_phy_bl_mynn_jmapplib_setup
 
  !-----------------------------------------------------------------------------
  !! Finalize
  
  subroutine atmos_phy_bl_mynn_jmapplib_finalize
    implicit none
 
    return
  end subroutine atmos_phy_bl_mynn_jmapplib_finalize
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_bl_mynn_jmapplib_tendency( &
       KA, KS, KE, IA, IS, IE, JA, JS, JE, &
       DENS, U, V, POTT, PROG,             &
       PRES, EXNER,                        &
       QDRY, QV, QC, QI,                   &
       SFC_DENS, SFC_PRES,                 &
       SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
       us, RLmo,                           &
       CZ, FZ, F2H, dt,                    &
       RHOU_t, RHOV_t, RHOT_t, RHOQV_t,    &
       RPROG_t,                            &
       Nu, Kh,                             &
       Zi, SFLX_BUOY                       )
    use scale_const, only: &
       cpdry => const_cpdry
    use scale_atmos_hydrometeor, only: &
       cp_vapor
    use scale_file_history, only: &
       file_history_in
#ifdef JMAPPLIB
    use pbl_mym_main, only: &
       pbl_mym_main_level3, &
       pbl_mym_main_level25
    use pbl_coupler, only: &
       pbl_coupler_flx_force_tend_run
    use pbl_diag, only: &
       pbl_diag_pbl_height
#endif
    implicit none
 
    integer, intent(in) :: ka, ks, ke
    integer, intent(in) :: ia, is, ie
    integer, intent(in) :: ja, js, je
 
    real(rp), intent(in) :: dens    (ka,ia,ja)
    real(rp), intent(in) :: u       (ka,ia,ja)
    real(rp), intent(in) :: v       (ka,ia,ja)
    real(rp), intent(in) :: pott    (ka,ia,ja)
    real(rp), intent(in) :: prog    (ka,ia,ja,atmos_phy_bl_mynn_jmapplib_ntracer)
    real(rp), intent(in) :: pres    (ka,ia,ja)
    real(rp), intent(in) :: exner   (ka,ia,ja)
    real(rp), intent(in) :: qdry    (ka,ia,ja)
    real(rp), intent(in) :: qv      (ka,ia,ja)
    real(rp), intent(in) :: qc      (ka,ia,ja)
    real(rp), intent(in) :: qi      (ka,ia,ja)
    real(rp), intent(in) :: sfc_dens(   ia,ja)
    real(rp), intent(in) :: sfc_pres(   ia,ja)
    real(rp), intent(in) :: sflx_mu (   ia,ja)
    real(rp), intent(in) :: sflx_mv (   ia,ja)
    real(rp), intent(in) :: sflx_sh (   ia,ja)
    real(rp), intent(in) :: sflx_qv (   ia,ja)
    real(rp), intent(in) :: us      (   ia,ja)
    real(rp), intent(in) :: rlmo    (   ia,ja)
 
    real(rp), intent(in) :: cz(  ka,ia,ja)
    real(rp), intent(in) :: fz(0:ka,ia,ja)
    real(rp), intent(in) :: f2h(ka,2,ia,ja)
    real(dp), intent(in) :: dt
 
    real(rp), intent(out) :: rhou_t (ka,ia,ja)
    real(rp), intent(out) :: rhov_t (ka,ia,ja)
    real(rp), intent(out) :: rhot_t (ka,ia,ja)
    real(rp), intent(out) :: rhoqv_t(ka,ia,ja)
    real(rp), intent(out) :: rprog_t(ka,ia,ja,atmos_phy_bl_mynn_jmapplib_ntracer)
    real(rp), intent(out) :: nu     (ka,ia,ja)
    real(rp), intent(out) :: kh     (ka,ia,ja)
    real(rp), intent(out) :: zi        (ia,ja)
    real(rp), intent(out) :: sflx_buoy (ia,ja)
 
#ifdef JMAPPLIB
    real(rp) :: l(ka,ia,ja)
 
    real(rp) :: qke(ks:ke_pbl)
    real(rp) :: qv_lc(ks:ke_pbl)
    real(rp) :: qc_lc(ks:ke_pbl)
    real(rp) :: qi_lc(ks:ke_pbl)
    real(rp) :: dens_lc(ks:ke_pbl)
    real(rp) :: dfm(ks:ke_pbl)
    real(rp) :: dfh(ks:ke_pbl)
    real(rp) :: taux_ex(ks:ke_pbl)
    real(rp) :: tauy_ex(ks:ke_pbl)
    real(rp) :: ftl_ex(ks:ke_pbl)
    real(rp) :: fqw_ex(ks:ke_pbl)
    real(rp) :: tend_qke(ks:ke_pbl)
    real(rp) :: tend_tsq(ks:ke_pbl)
    real(rp) :: tend_qsq(ks:ke_pbl)
    real(rp) :: tend_cov(ks:ke_pbl)
    real(rp) :: tend_u(ks:ke_pbl)
    real(rp) :: tend_v(ks:ke_pbl)
    real(rp) :: tend_pt(ks:ke_pbl)
    real(rp) :: tend_qv(ks:ke_pbl)
    real(rp) :: z_f(ks:ke_pbl)
    real(rp) :: dz_f(ks:ke_pbl)
    real(rp) :: rdz_f(ks:ke_pbl)
    real(rp) :: rdz_h(ks:ke_pbl)
    real(rp) :: h2f_m(ks:ke_pbl)
    real(rp) :: h2f_p(ks:ke_pbl)
    real(rp) :: fb_surf
    real(rp) :: rho_ov_rhoa
    real(rp) :: sflx_u
    real(rp) :: sflx_v
    real(rp) :: sflx_pt
    real(rp) :: sflx_q
    real(rp) :: cptot
 
    real(rp) :: diss(ka,ia,ja)
    real(rp) :: dummy(ka,ia,ja)
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    log_progress(*) "atmosphere / physics / pbl / MYNN-JMAPPLIB"
 
    h2f_m(:) = 0.5_rp
    h2f_p(:) = 0.5_rp
 
    !$omp parallel do &
    !$omp private(qke,qv_lc,qc_lc,qi_lc,dens_lc,taux_ex,tauy_ex,ftl_ex,fqw_ex,&
    !$omp         tend_qke,tend_tsq,tend_qsq,tend_cov,tend_u,tend_v,tend_pt,tend_qv, &
    !$omp         z_f,dz_f,rdz_f,rdz_h, &
    !$omp         rho_ov_rhoa,SFLX_U,SFLX_V,SFLX_PT,SFLX_Q,CPtot)
    do j = js, je
    do i = is, ie
 
       do k = ks, ke_pbl
          qke(k) = prog(k,i,j,i_tke) * 2.0_rp
          rho_ov_rhoa = 1.0_rp / ( qdry(k,i,j) + qv(k,i,j) )
          qv_lc(k) = qv(k,i,j) * rho_ov_rhoa
          qc_lc(k) = qc(k,i,j) * rho_ov_rhoa
          qi_lc(k) = qi(k,i,j) * rho_ov_rhoa
          z_f(k) = cz(k,i,j) - fz(ks-1,i,j)
          dz_f(k) = fz(k,i,j) - fz(k-1,i,j)
          rdz_f(k) = 1.0_rp / dz_f(k)
          rdz_h(k) = 1.0_rp / ( cz(k+1,i,j) - cz(k,i,j) )
       end do
 
       cptot = cpdry + sflx_qv(i,j) * ( cp_vapor - cpdry )
 
       sflx_u  = sflx_mu(i,j) / sfc_dens(i,j)
       sflx_v  = sflx_mv(i,j) / sfc_dens(i,j)
       sflx_pt = sflx_sh(i,j) / ( cptot * exner(ks,i,j) * sfc_dens(i,j) )
       sflx_q  = sflx_qv(i,j) / sfc_dens(i,j)
 
       select case ( atmos_phy_bl_mynn_jmapplib_level )
       case ( "3" )
          call pbl_mym_main_level3( &
               sflx_u, sflx_v, sflx_pt, sflx_q, & ! (in)
               us(i,j), rlmo(i,j), sfc_pres(i,j), & ! (in)
               u(ks:ke_pbl,i,j), v(ks:ke_pbl,i,j), pott(ks:ke_pbl,i,j), & ! (in)
               qv_lc, qc_lc, qi_lc, pres(ks:ke_pbl,i,j), exner(ks:ke_pbl,i,j), & ! (in)
               qke(:), prog(ks:ke_pbl,i,j,i_tsq), prog(ks:ke_pbl,i,j,i_qsq), prog(ks:ke_pbl,i,j,i_cov), & ! (in)
               z_f, dz_f, rdz_f, rdz_h, f2h(ks:ke_pbl,2,i,j), f2h(ks:ke_pbl,1,i,j), h2f_m, h2f_p, & ! (in)
               sflx_buoy(i,j), & ! (out)
               dfm(:), dfh(:), tend_qke(:), tend_tsq(:), tend_qsq(:), tend_cov(:), & ! (out)
               taux_ex, tauy_ex, ftl_ex, fqw_ex, & ! (out)
               l_e = l(ks:ke_pbl,i,j), tke_dis = diss(ks:ke_pbl,i,j), tsq_dis=dummy(ks:ke_pbl,i,j) ) ! (out)
       case ( "2.5" )
          call pbl_mym_main_level25( &
               sflx_u, sflx_v, sflx_pt, sflx_q, & ! (in)
               us(i,j), rlmo(i,j), sfc_pres(i,j), & ! (in)
               u(ks:ke_pbl,i,j), v(ks:ke_pbl,i,j), pott(ks:ke_pbl,i,j), & ! (in)
               qv_lc, qc_lc, qi_lc, pres(ks:ke_pbl,i,j), exner(ks:ke_pbl,i,j), & ! (in)
               qke(:), prog(ks:ke_pbl,i,j,i_tsq), prog(ks:ke_pbl,i,j,i_qsq), prog(ks:ke_pbl,i,j,i_cov), & ! (in)
               z_f, dz_f, rdz_f, rdz_h, f2h(ks:ke_pbl,2,i,j), f2h(ks:ke_pbl,1,i,j), h2f_m, h2f_p, & ! (in)
               sflx_buoy(i,j), & ! (out)
               dfm(:), dfh(:), tend_qke(:), tend_tsq(:), tend_qsq(:), tend_cov(:), & ! (out)
               taux_ex, tauy_ex, ftl_ex, fqw_ex, & ! (out)
               l_e = l(ks:ke_pbl,i,j), tke_dis = diss(ks:ke_pbl,i,j), tsq_dis=dummy(ks:ke_pbl,i,j) ) ! (out)
          do k = ks, ke_pbl
             tend_tsq(k) = ( tend_tsq(k) - prog(k,i,j,i_tsq) ) / dt
             tend_qsq(k) = ( tend_qsq(k) - prog(k,i,j,i_qsq) ) / dt
             tend_cov(k) = ( tend_cov(k) - prog(k,i,j,i_cov) ) / dt
          end do
       end select
 
       do k = ks, ke_pbl
          dens_lc(k) = dens(k,i,j) * ( qdry(k,i,j) + qv(k,i,j) )
       end do
 
       call pbl_coupler_flx_force_tend_run( &
            sfc_dens(i,j), sflx_u, sflx_v, sflx_pt, sflx_q, & ! (in)
            dfm(:), dfh(:), dens_lc(:),       & ! (in)
            taux_ex, tauy_ex, ftl_ex, fqw_ex, & ! (in)
            f2h(ks:ke_pbl,2,i,j), f2h(ks:ke_pbl,1,i,j), rdz_f, rdz_h, & ! (in)
            tend_u, tend_v, tend_pt, tend_qv ) ! (out)
 
       do k = ks, ke_pbl
          rhou_t(k,i,j) = tend_u(k) * dens(k,i,j)
          rhov_t(k,i,j) = tend_v(k) * dens(k,i,j)
          rhot_t(k,i,j) = tend_pt(k) * dens(k,i,j)
          rhoqv_t(k,i,j) = tend_qv(k) * dens(k,i,j) * ( qdry(k,i,j) + qv(k,i,j) )
          rprog_t(k,i,j,i_tke) = tend_qke(k) * dens(k,i,j) * 0.5_rp
          rprog_t(k,i,j,i_tsq) = tend_tsq(k) * dens(k,i,j)
          rprog_t(k,i,j,i_qsq) = tend_qsq(k) * dens(k,i,j)
          rprog_t(k,i,j,i_cov) = tend_cov(k) * dens(k,i,j)
       end do
       rhou_t(ks,i,j) = rhou_t(ks,i,j) - sflx_mu(i,j) * rdz_f(ks)
       rhov_t(ks,i,j) = rhov_t(ks,i,j) - sflx_mv(i,j) * rdz_f(ks)
       rhot_t(ks,i,j) = rhot_t(ks,i,j) - sflx_pt * sfc_dens(i,j) * rdz_f(ks)
       rhoqv_t(ks,i,j) = rhoqv_t(ks,i,j) - sflx_qv(i,j) * rdz_f(ks)
       do k = ke_pbl+1, ke
          rhou_t(k,i,j) = 0.0_rp
          rhov_t(k,i,j) = 0.0_rp
          rhot_t(k,i,j) = 0.0_rp
          rhoqv_t(k,i,j) = 0.0_rp
          rprog_t(k,i,j,i_tke) = 0.0_rp
          rprog_t(k,i,j,i_tsq) = 0.0_rp
          rprog_t(k,i,j,i_qsq) = 0.0_rp
          rprog_t(k,i,j,i_cov) = 0.0_rp
       end do
 
       nu(ks-1,i,j) = 0.0_rp
       kh(ks-1,i,j) = 0.0_rp
       do k = ks, ke_pbl-1
          nu(k,i,j) = dfm(k+1) * f2h(k,1,i,j) + dfm(k) * f2h(k,2,i,j)
          kh(k,i,j) = dfh(k+1) * f2h(k,1,i,j) + dfh(k) * f2h(k,2,i,j)
       end do
       do k = ke_pbl, ke
          nu(k,i,j) = 0.0_rp
          kh(k,i,j) = 0.0_rp
       end do
 
       call pbl_diag_pbl_height( &
            sflx_pt, us(i,j), rlmo(i,j), & ! (in)
            u(ks:ke_pbl,i,j), v(ks:ke_pbl,i,j), pott(ks:ke_pbl,i,j), qv_lc, z_f, & ! (in)
            zi(i,j) ) ! (out)
 
    end do
    end do
 
 
    l(ke_pbl+1:ke,:,:) = undef
    call file_history_in(l(:,:,:), 'L_mix_MYNN', 'minxing length', 'm', fill_halo=.true.)
 
    diss(ks:ke_pbl,:,:) = - diss(ks:ke_pbl,:,:)
    diss(ke_pbl+1:ke,:,:) = undef
    call file_history_in(diss(:,:,:), 'TKE_diss_MYNN', 'TKE dissipation', 'm2/s3', fill_halo=.true.)
#endif
 
    return
  end subroutine atmos_phy_bl_mynn_jmapplib_tendency
 
end module scale_atmos_phy_bl_mynn_jmapplib