scale_atmos_phy_tb_dns.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_phy_tb_dns
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_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_dns_config
  public :: atmos_phy_tb_dns_setup
  public :: atmos_phy_tb_dns

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  real(RP), private :: atmos_phy_tb_dns_nu = 1.512e-5_rp ! [m2/s] kinematic viscosity coefficient for air at 20degC
! real(RP), private :: mu = 1.8E-5_RP   ! [m2/s] molecular diffusive coefficient for air at 20degC
  real(RP), private :: atmos_phy_tb_dns_mu = 1.512e-5_rp ! same as NU (needed based on hyposes. see Mellado 2010)

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_dns_config( &
       TYPE_TB,  &
       I_TKE_out )
    use scale_prc, only: &
       prc_abort
    implicit none

    character(len=*), intent(in)  :: TYPE_TB
    integer,          intent(out) :: I_TKE_out
    !---------------------------------------------------------------------------

    log_newline
    log_info("ATMOS_PHY_TB_dns_config",*) 'Setup'
    log_info("ATMOS_PHY_TB_dns_config",*) 'Tracers for Deardorff (1980) 1.5th TKE Model'

    if ( type_tb /= 'DNS' ) then
       log_error("ATMOS_PHY_TB_dns_config",*) 'ATMOS_PHY_TB_TYPE is not DNS. Check!'
       call prc_abort
    endif

    i_tke_out = -1

    return
  end subroutine atmos_phy_tb_dns_config
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_dns_setup( &
       CDZ, CDX, CDY, CZ )
    use scale_prc, only: &
       prc_abort
    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)

    namelist / param_atmos_phy_tb_dns / &
       atmos_phy_tb_dns_nu, &
       atmos_phy_tb_dns_mu

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

    log_newline
    log_info("ATMOS_PHY_TB_dns_setup",*) 'Setup'
    log_info("ATMOS_PHY_TB_dns_setup",*) 'Eddy Viscocity Model for DNS'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_tb_dns,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_TB_dns_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_TB_dns_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_DNS. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_tb_dns)

    return
  end subroutine atmos_phy_tb_dns_setup

  !-----------------------------------------------------------------------------
  subroutine atmos_phy_tb_dns( &
       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_atmos_grid_cartesc_index
    use scale_tracer
    use scale_const, only: &
       grav => const_grav
    use scale_atmos_grid_cartesc, only: &
       rcdz => atmos_grid_cartesc_rcdz, &
       rcdx => atmos_grid_cartesc_rcdx, &
       rcdy => atmos_grid_cartesc_rcdy, &
       rfdz => atmos_grid_cartesc_rfdz, &
       rfdx => atmos_grid_cartesc_rfdx, &
       rfdy => atmos_grid_cartesc_rfdy
    implicit none

    ! SGS flux
    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) :: 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

    real(RP) :: POTT(ka,ia,ja)

    integer :: IIS, IIE
    integer :: JJS, JJE

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

    log_progress(*) 'atmosphere / physics / turbulence / DNS'

#ifdef DEBUG
    qflx_sgs_momz(:,:,:,:)   = undef
    qflx_sgs_momx(:,:,:,:)   = undef
    qflx_sgs_momy(:,:,:,:)   = undef
    qflx_sgs_rhot(:,:,:,:)   = undef
    qflx_sgs_rhoq(:,:,:,:,:) = undef

    pott(:,:,:) = undef
#endif

    nu(:,:,:) = 0.0_rp
    ri(:,:,:) = 0.0_rp
    pr(:,:,:) = 1.0_rp

    ! potential temperature
    do j = js-1, je+1
    do i = is-1, ie+1
    do k = ks, ke
       pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
    enddo
    enddo
    enddo

    !##### Start Upadate #####
    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

       !##### momentum equation (z) #####
       ! (cell center)
       do j = jjs, jje
       do i = iis, iie
       do k = ks+1, ke-1
          qflx_sgs_momz(k,i,j,zdir) = -atmos_phy_tb_dns_nu * ( momz(k,i,j)-momz(k-1,i,j) ) * rcdz(k)
       enddo
       enddo
       enddo

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_momz(ks,i,j,zdir) = 0.0_rp ! bottom boundary
          qflx_sgs_momz(ke,i,j,zdir) = 0.0_rp ! top boundary
       enddo
       enddo

       ! (y edge)
       do j = jjs,   jje
       do i = iis-1, iie
       do k = ks, ke-1
          qflx_sgs_momz(k,i,j,xdir) = -atmos_phy_tb_dns_nu * ( momz(k,i+1,j)-momz(k,i,j) ) * rfdx(i) * mapf(i,j,1,i_xy)
       enddo
       enddo
       enddo

       ! (x edge)
       do j = jjs-1, jje
       do i = iis,   iie
       do k = ks, ke-1
          qflx_sgs_momz(k,i,j,ydir) = -atmos_phy_tb_dns_nu * ( momz(k,i,j+1)-momz(k,i,j) ) * rfdy(j) * mapf(i,j,2,i_xy)
       enddo
       enddo
       enddo

       !##### momentum equation (x) #####
       ! (y edge)
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke-1
          qflx_sgs_momx(k,i,j,zdir) = -atmos_phy_tb_dns_nu * ( momx(k+1,i,j)-momx(k,i,j) ) * rfdz(k)
       enddo
       enddo
       enddo

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
          qflx_sgs_momx(ke  ,i,j,zdir) = 0.0_rp ! top boundary
       enddo
       enddo

       ! (cell center)
       do j = jjs, jje
       do i = iis, iie+1
       do k = ks, ke
          qflx_sgs_momx(k,i,j,xdir) = -atmos_phy_tb_dns_nu * ( momx(k,i,j)-momx(k,i-1,j) ) * rcdx(i) * mapf(i,j,1,i_uy)
       enddo
       enddo
       enddo

       ! (z edge)
       do j = jjs-1, jje
       do i = iis,   iie
       do k = ks, ke
          qflx_sgs_momx(k,i,j,ydir) = -atmos_phy_tb_dns_nu * ( momx(k,i,j+1)-momx(k,i,j) ) * rfdy(j) * mapf(i,j,2,i_uy)
       enddo
       enddo
       enddo

       !##### momentum equation (y) #####

       ! (x edge)
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke-1
          qflx_sgs_momy(k,i,j,zdir) = -atmos_phy_tb_dns_nu * ( momy(k+1,i,j)-momy(k,i,j) ) * rfdz(k)
       enddo
       enddo
       enddo

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
          qflx_sgs_momy(ke  ,i,j,zdir) = 0.0_rp ! top boundary
       enddo
       enddo

       ! (z edge)
       do j = jjs,   jje
       do i = iis-1, iie
       do k = ks, ke
          qflx_sgs_momy(k,i,j,xdir) = -atmos_phy_tb_dns_nu * ( momy(k,i+1,j)-momy(k,i,j) ) * rfdx(i) * mapf(i,j,1,i_xv)
       enddo
       enddo
       enddo

       ! (z-x plane)
       do j = jjs, jje+1
       do i = iis, iie
       do k = ks, ke
          qflx_sgs_momy(k,i,j,ydir) = -atmos_phy_tb_dns_nu * ( momy(k,i,j)-momy(k,i,j-1) ) * rcdy(j) * mapf(i,j,2,i_xv)
       enddo
       enddo
       enddo

       !##### Thermodynamic Equation #####

       ! at x, y ,w
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke-1
          qflx_sgs_rhot(k,i,j,zdir) = -0.5_rp * ( dens(k+1,i,j)+dens(k,i,j) ) &
                                    * atmos_phy_tb_dns_mu * ( pott(k+1,i,j)-pott(k,i,j) ) * rfdz(k)
       enddo
       enddo
       enddo

       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_rhot(ks-1,i,j,zdir) = 0.0_rp
          qflx_sgs_rhot(ke  ,i,j,zdir) = 0.0_rp
       enddo
       enddo

       ! at u, y, z
       do j = jjs,   jje
       do i = iis-1, iie
       do k = ks, ke
          qflx_sgs_rhot(k,i,j,xdir) = -0.5_rp * ( dens(k,i+1,j)+dens(k,i,j) ) &
                                    * atmos_phy_tb_dns_mu * ( pott(k,i+1,j)-pott(k,i,j) ) * rfdx(i) * mapf(i,j,1,i_xy)
       enddo
       enddo
       enddo

       ! at x, v, z
       do j = jjs-1, jje
       do i = iis,   iie
       do k = ks, ke
          qflx_sgs_rhot(k,i,j,ydir) = -0.5_rp * ( dens(k,i,j+1)+dens(k,i,j) ) &
                                    * atmos_phy_tb_dns_mu * ( pott(k,i,j+1)-pott(k,i,j) ) * rfdy(j) * mapf(i,j,2,i_xy)
       enddo
       enddo
       enddo

    enddo
    enddo

    !##### Tracers #####
    do iq = 1, qa

    if ( .not. tracer_advc(iq) ) cycle

    do jjs = js, je, jblock
    jje = jjs+jblock-1
    do iis = is, ie, iblock
    iie = iis+iblock-1

       ! at x, y ,w
       do j = jjs, jje
       do i = iis, iie
       do k = ks, ke-1
          qflx_sgs_rhoq(k,i,j,zdir,iq) = -0.5_rp * ( dens(k+1,i,j)+dens(k,i,j) ) &
                                       * atmos_phy_tb_dns_mu * ( qtrc(k+1,i,j,iq)-qtrc(k,i,j,iq) ) * rfdz(k)
       enddo
       enddo
       enddo
       do j = jjs, jje
       do i = iis, iie
          qflx_sgs_rhoq(ks-1,i,j,zdir,iq) = 0.0_rp
          qflx_sgs_rhoq(ke  ,i,j,zdir,iq) = 0.0_rp
       enddo
       enddo

       ! at u, y, z
       do j = jjs,   jje
       do i = iis-1, iie
       do k = ks,   ke
          qflx_sgs_rhoq(k,i,j,xdir,iq) = -0.5_rp * ( dens(k,i+1,j)+dens(k,i,j) ) &
                                       * atmos_phy_tb_dns_mu * ( qtrc(k,i+1,j,iq)-qtrc(k,i,j,iq) ) * rfdx(i) * mapf(i,j,1,i_xy)
       enddo
       enddo
       enddo

       ! at x, v, z
       do j = jjs-1, jje
       do i = iis,   iie
       do k = ks,   ke
          qflx_sgs_rhoq(k,i,j,ydir,iq) = -0.5_rp * ( dens(k,i,j+1)+dens(k,i,j) ) &
                                       * atmos_phy_tb_dns_mu * ( qtrc(k,i,j+1,iq)-qtrc(k,i,j,iq) ) * rfdy(j) * mapf(i,j,2,i_xy)
       enddo
       enddo
       enddo

    enddo
    enddo

    enddo ! scalar quantities loop

    return
  end subroutine atmos_phy_tb_dns

end module scale_atmos_phy_tb_dns