scale_atmos_phy_tb_hybrid.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_phy_tb_hybrid
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer

#ifdef DEBUG
  use scale_debug, only: &
     check
  use scale_const, only: &
     undef  => const_undef, &
     iundef => const_undef2
#endif
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_tb_hybrid_config
  public :: atmos_phy_tb_hybrid_setup
  public :: atmos_phy_tb_hybrid

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  abstract interface
     subroutine tb( &
       qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
       qflx_sgs_rhot, qflx_sgs_rhoq,                &
       RHOQ_t, nu_C, Ri, Pr,                        &
       MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2,      &
       SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q,  &
       GSQRT, J13G, J23G, J33G, MAPF, dt            )
       use scale_precision
       use scale_grid_index
       use scale_tracer
       implicit none

       real(RP), intent(out)   :: qflx_sgs_momz(KA,IA,JA,3)
       real(RP), intent(out)   :: qflx_sgs_momx(KA,IA,JA,3)
       real(RP), intent(out)   :: qflx_sgs_momy(KA,IA,JA,3)
       real(RP), intent(out)   :: qflx_sgs_rhot(KA,IA,JA,3)
       real(RP), intent(out)   :: qflx_sgs_rhoq(KA,IA,JA,3,QA)

       real(RP), intent(inout) :: RHOQ_t       (KA,IA,JA,QA) ! tendency of rho * QTRC

       real(RP), intent(out)   :: nu_C         (KA,IA,JA)    ! eddy viscosity (center)
       real(RP), intent(out)   :: Ri           (KA,IA,JA)    ! Richardson number
       real(RP), intent(out)   :: Pr           (KA,IA,JA)    ! Prantle number

       real(RP), intent(in)    :: MOMZ         (KA,IA,JA)
       real(RP), intent(in)    :: MOMX         (KA,IA,JA)
       real(RP), intent(in)    :: MOMY         (KA,IA,JA)
       real(RP), intent(in)    :: RHOT         (KA,IA,JA)
       real(RP), intent(in)    :: DENS         (KA,IA,JA)
       real(RP), intent(in)    :: QTRC         (KA,IA,JA,QA)
       real(RP), intent(in)    :: N2           (KA,IA,JA)

       real(RP), intent(in)    :: SFLX_MW      (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_Q       (IA,JA,QA)

       real(RP), intent(in)    :: GSQRT        (KA,IA,JA,7)
       real(RP), intent(in)    :: J13G         (KA,IA,JA,7)
       real(RP), intent(in)    :: J23G         (KA,IA,JA,7)
       real(RP), intent(in)    :: J33G
       real(RP), intent(in)    :: MAPF         (IA,JA,2,4)
       real(DP), intent(in)    :: dt
     end subroutine tb

     subroutine su( &
          CDZ, CDX, CDY, CZ )
       use scale_precision
       use scale_grid_index
       use scale_tracer
       implicit none

       real(RP), intent(in) :: CDZ(KA)
       real(RP), intent(in) :: CDX(IA)
       real(RP), intent(in) :: CDY(JA)
       real(RP), intent(in) :: CZ (KA,IA,JA)
     end subroutine su
  end interface

  procedure(tb), pointer :: sgs_tb       => null()
  procedure(tb), pointer :: pbl_tb       => null()

  procedure(su), pointer :: sgs_tb_setup => null()
  procedure(su), pointer :: pbl_tb_setup => null()

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private              :: atmos_phy_tb_hybrid_sgs_dx = 100.0_rp 
  real(RP), private              :: atmos_phy_tb_hybrid_pbl_dx = 500.0_rp 

  real(RP), private, allocatable :: frac_sgs    (:,:)
  real(RP), private, allocatable :: frac_pbl    (:,:)
  real(RP), private, allocatable :: frac_sgs_tke(:,:)
  real(RP), private, allocatable :: frac_pbl_tke(:,:)

  integer,  private              :: i_tke_sgs, i_tke_pbl

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_hybrid_config( &
       TB_TYPE,  &
       I_TKE_out )
    use scale_process, only: &
       prc_mpistop
    use scale_atmos_phy_tb_smg, only: &
       atmos_phy_tb_smg_config, &
       atmos_phy_tb_smg_setup, &
       atmos_phy_tb_smg
    use scale_atmos_phy_tb_mynn, only: &
       atmos_phy_tb_mynn_config, &
       atmos_phy_tb_mynn_setup, &
       atmos_phy_tb_mynn
    implicit none

    character(len=*), intent(in)  :: tb_type
    integer,          intent(out) :: i_tke_out

    character(len=H_SHORT) :: atmos_phy_tb_hybrid_sgs_type = 'SMAGORINSKY'
    character(len=H_SHORT) :: atmos_phy_tb_hybrid_pbl_type = 'MYNN'

    namelist / param_atmos_phy_tb_hybrid / &
       atmos_phy_tb_hybrid_sgs_dx, &
       atmos_phy_tb_hybrid_pbl_dx, &
       atmos_phy_tb_hybrid_sgs_type, &
       atmos_phy_tb_hybrid_pbl_type

    integer  :: ierr
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
    if( io_l ) write(io_fid_log,*) '*** Tracers for SGS-parameterization hybrid Model'

    if ( tb_type /= 'HYBRID' ) then
       write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not HYBRID. Check!'
       call prc_mpistop
    endif

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_tb_hybrid,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_HYBRID. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_tb_hybrid)

    select case( atmos_phy_tb_hybrid_sgs_type )
    case('SMAGORINSKY')
       call atmos_phy_tb_smg_config( &
            atmos_phy_tb_hybrid_sgs_type, &
            i_tke_sgs )
       sgs_tb       => atmos_phy_tb_smg
       sgs_tb_setup => atmos_phy_tb_smg_setup
    case default
       write(*,*) 'xxx ATMOS_PHY_TB_HYBRID_SGS_TYPE is invalid'
       call prc_mpistop
    end select

    select case( atmos_phy_tb_hybrid_pbl_type )
    case('MYNN')
       call atmos_phy_tb_mynn_config( &
            atmos_phy_tb_hybrid_pbl_type, &
            i_tke_pbl )
       pbl_tb       => atmos_phy_tb_mynn
       pbl_tb_setup => atmos_phy_tb_mynn_setup
    case default
       write(*,*) 'xxx ATMOS_PHY_TB_HYBRID_PBL_TYPE is invalid'
       call prc_mpistop
    end select

    i_tke_out = i_tke_pbl

    return
  end subroutine atmos_phy_tb_hybrid_config

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_hybrid_setup( &
       CDZ, CDX, CDY, CZ )
    implicit none

    real(RP), intent(in) :: cdz(ka)
    real(RP), intent(in) :: cdx(ia)
    real(RP), intent(in) :: cdy(ja)
    real(RP), intent(in) :: cz (ka,ia,ja)

    real(RP) :: dxy

    integer  :: i, j
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
    if( io_l ) write(io_fid_log,*) '*** SGS-parameterization hybrid Model'

    call sgs_tb_setup( cdz, cdx, cdy, cz )
    call pbl_tb_setup( cdz, cdx, cdy, cz )

    allocate( frac_sgs(ia,ja) )
    allocate( frac_pbl(ia,ja) )
    allocate( frac_sgs_tke(ia,ja) )
    allocate( frac_pbl_tke(ia,ja) )

    do j = 1, ja
    do i = 1, ia
       dxy = sqrt( 0.5_rp * ( cdx(i)**2 + cdy(j)**2 ) )

       frac_pbl(i,j) = ( dxy - atmos_phy_tb_hybrid_sgs_dx ) / ( atmos_phy_tb_hybrid_pbl_dx - atmos_phy_tb_hybrid_sgs_dx )
       frac_pbl(i,j) = min( 1.0_rp, max( 0.0_rp, frac_pbl(i,j) ) )
       frac_sgs(i,j) = 1.0_rp - frac_pbl(i,j)
    end do
    end do

    return
  end subroutine atmos_phy_tb_hybrid_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_hybrid( &
       qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
       qflx_sgs_rhot, qflx_sgs_rhoq,                &
       RHOQ_t,                                      &
       Nu, Ri, Pr,                                  &
       MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2,      &
       SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q,  &
       GSQRT, J13G, J23G, J33G, MAPF, dt            )
    use scale_precision
    use scale_grid_index
    use scale_tracer
    use scale_const, only: &
       grav => const_grav
    implicit none

    real(RP), intent(out)   :: qflx_sgs_momz(ka,ia,ja,3)
    real(RP), intent(out)   :: qflx_sgs_momx(ka,ia,ja,3)
    real(RP), intent(out)   :: qflx_sgs_momy(ka,ia,ja,3)
    real(RP), intent(out)   :: qflx_sgs_rhot(ka,ia,ja,3)
    real(RP), intent(out)   :: qflx_sgs_rhoq(ka,ia,ja,3,qa)

    real(RP), intent(inout) :: rhoq_t       (ka,ia,ja,qa) ! tendency of rho * QTRC

    real(RP), intent(out)   :: nu           (ka,ia,ja) ! eddy viscosity (center)
    real(RP), intent(out)   :: pr           (ka,ia,ja) ! Plandtle number
    real(RP), intent(out)   :: ri           (ka,ia,ja) ! Richardson number

    real(RP), intent(in)    :: momz         (ka,ia,ja)
    real(RP), intent(in)    :: momx         (ka,ia,ja)
    real(RP), intent(in)    :: momy         (ka,ia,ja)
    real(RP), intent(in)    :: rhot         (ka,ia,ja)
    real(RP), intent(in)    :: dens         (ka,ia,ja)
    real(RP), intent(in)    :: qtrc         (ka,ia,ja,qa)
    real(RP), intent(in)    :: n2           (ka,ia,ja)

    real(RP), intent(in)    :: sflx_mw      (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_q       (ia,ja,qa)

    real(RP), intent(in)    :: gsqrt        (ka,ia,ja,7)
    real(RP), intent(in)    :: j13g         (ka,ia,ja,7)
    real(RP), intent(in)    :: j23g         (ka,ia,ja,7)
    real(RP), intent(in)    :: j33g
    real(RP), intent(in)    :: mapf         (ia,ja,2,4)
    real(DP), intent(in)    :: dt

    real(RP) :: w_qflx_sgs_momz(ka,ia,ja,3,2)
    real(RP) :: w_qflx_sgs_momx(ka,ia,ja,3,2)
    real(RP) :: w_qflx_sgs_momy(ka,ia,ja,3,2)
    real(RP) :: w_qflx_sgs_rhot(ka,ia,ja,3,2)
    real(RP) :: w_qflx_sgs_rhoq(ka,ia,ja,3,qa,2)

    real(RP) :: w_nu           (ka,ia,ja,2)
    real(RP) :: w_ri           (ka,ia,ja,2)
    real(RP) :: w_pr           (ka,ia,ja,2)

    integer  :: k, i, j, iq
    !---------------------------------------------------------------------------

    call sgs_tb( w_qflx_sgs_momz(:,:,:,:,1), w_qflx_sgs_momx(:,:,:,:,1), & ! [OUT]
                 w_qflx_sgs_momy(:,:,:,:,1), w_qflx_sgs_rhot(:,:,:,:,1), & ! [OUT]
                 w_qflx_sgs_rhoq(:,:,:,:,:,1),                           & ! [OUT]
                 rhoq_t,                                                 & ! [INOUT]
                 w_nu(:,:,:,1), w_ri(:,:,:,1), w_pr(:,:,:,1),            & ! [OUT]
                 momz, momx, momy, rhot, dens, qtrc, n2,                 & ! [IN]
                 sflx_mw, sflx_mu, sflx_mv, sflx_sh, sflx_q,             & ! [IN]
                 gsqrt, j13g, j23g, j33g, mapf, dt                       ) ! [IN]

    call pbl_tb( w_qflx_sgs_momz(:,:,:,:,2), w_qflx_sgs_momx(:,:,:,:,2), & ! [OUT]
                 w_qflx_sgs_momy(:,:,:,:,2), w_qflx_sgs_rhot(:,:,:,:,2), & ! [OUT]
                 w_qflx_sgs_rhoq(:,:,:,:,:,2),                           & ! [OUT]
                 rhoq_t,                                                 & ! [INOUT]
                 w_nu(:,:,:,2), w_ri(:,:,:,2), w_pr(:,:,:,2),            & ! [OUT]
                 momz, momx, momy, rhot, dens, qtrc, n2,                 & ! [IN]
                 sflx_mw, sflx_mu, sflx_mv, sflx_sh, sflx_q,             & ! [IN]
                 gsqrt, j13g, j23g, j33g, mapf, dt                       ) ! [IN]

    !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
    !$omp shared(JA,IA,KS,KE,qflx_sgs_momz,w_qflx_sgs_momz,frac_sgs,frac_pbl)
    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       qflx_sgs_momz(k,i,j,zdir) = w_qflx_sgs_momz(k,i,j,zdir,1) * frac_sgs(i,j) &
                                 + w_qflx_sgs_momz(k,i,j,zdir,2) * frac_pbl(i,j)
       qflx_sgs_momz(k,i,j,xdir) = w_qflx_sgs_momz(k,i,j,xdir,1)
       qflx_sgs_momz(k,i,j,ydir) = w_qflx_sgs_momz(k,i,j,ydir,1)
    end do
    end do
    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       qflx_sgs_momx(k,i,j,zdir) = w_qflx_sgs_momx(k,i,j,zdir,1) * frac_sgs(i,j) &
                                 + w_qflx_sgs_momx(k,i,j,zdir,2) * frac_pbl(i,j)
       qflx_sgs_momx(k,i,j,xdir) = w_qflx_sgs_momx(k,i,j,xdir,1)
       qflx_sgs_momx(k,i,j,ydir) = w_qflx_sgs_momx(k,i,j,ydir,1)
    end do
    end do
    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       qflx_sgs_momy(k,i,j,zdir) = w_qflx_sgs_momy(k,i,j,zdir,1) * frac_sgs(i,j) &
                                 + w_qflx_sgs_momy(k,i,j,zdir,2) * frac_pbl(i,j)
       qflx_sgs_momy(k,i,j,xdir) = w_qflx_sgs_momy(k,i,j,xdir,1)
       qflx_sgs_momy(k,i,j,ydir) = w_qflx_sgs_momy(k,i,j,ydir,1)
    end do
    end do
    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       qflx_sgs_rhot(k,i,j,zdir) = w_qflx_sgs_rhot(k,i,j,zdir,1) * frac_sgs(i,j) &
                                 + w_qflx_sgs_rhot(k,i,j,zdir,2) * frac_pbl(i,j)
       qflx_sgs_rhot(k,i,j,xdir) = w_qflx_sgs_rhot(k,i,j,xdir,1)
       qflx_sgs_rhot(k,i,j,ydir) = w_qflx_sgs_rhot(k,i,j,ydir,1)
    end do
    end do
    end do

    !$omp parallel do default(none)                                                        &
    !$omp shared(JA,IA,KS,KE,qflx_sgs_rhoq,w_qflx_sgs_rhoq,frac_sgs,frac_pbl,QA,I_TKE_SGS) &
    !$omp shared(I_TKE_PBL,TRACER_ADVC)                                                    &
    !$omp private(i,j,k,iq) OMP_SCHEDULE_
    do iq = 1, qa

       if ( iq == i_tke_sgs .or. iq == i_tke_pbl ) then
          qflx_sgs_rhoq(:,:,:,:,iq) = 0.0_rp
          cycle
       end if
       if ( .not. tracer_advc(iq) ) cycle

       do j = 1, ja
       do i = 1, ia
       do k = ks, ke
          qflx_sgs_rhoq(k,i,j,zdir,iq) = w_qflx_sgs_rhoq(k,i,j,zdir,iq,1) * frac_sgs(i,j) &
                                       + w_qflx_sgs_rhoq(k,i,j,zdir,iq,2) * frac_pbl(i,j)
          qflx_sgs_rhoq(k,i,j,xdir,iq) = w_qflx_sgs_rhoq(k,i,j,xdir,iq,1)
          qflx_sgs_rhoq(k,i,j,ydir,iq) = w_qflx_sgs_rhoq(k,i,j,ydir,iq,1)
       end do
       end do
       end do

    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       nu(k,i,j) = w_nu(k,i,j,1) * frac_sgs(i,j) &
                 + w_nu(k,i,j,2) * frac_pbl(i,j)
    end do
    end do
    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       ri(k,i,j) = w_ri(k,i,j,1) * frac_sgs(i,j) &
                 + w_ri(k,i,j,2) * frac_pbl(i,j)
    end do
    end do
    end do

    do j = 1, ja
    do i = 1, ia
    do k = ks, ke
       pr(k,i,j) = w_pr(k,i,j,1) * frac_sgs(i,j) &
                 + w_pr(k,i,j,2) * frac_pbl(i,j)
    end do
    end do
    end do

    return
  end subroutine atmos_phy_tb_hybrid

end module scale_atmos_phy_tb_hybrid