scale_atmos_dyn_tstep_tracer_fvm_heve.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_atmos_dyn_tstep_tracer_fvm_heve
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_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_dyn_tstep_tracer_fvm_heve_setup
  public :: atmos_dyn_tstep_tracer_fvm_heve

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains

  !-----------------------------------------------------------------------------
  subroutine atmos_dyn_tstep_tracer_fvm_heve_setup( type )
    use scale_process, only: &
       prc_mpistop
    implicit none
    character(len=*), intent(in) :: type

    if ( type /= 'FVM-HEVE' ) then
       write(*,*) 'xxx Tstep_tracer_type is not "FVM-HEVE"!'
       call prc_mpistop
    end if

    return
  end subroutine atmos_dyn_tstep_tracer_fvm_heve_setup

  subroutine atmos_dyn_tstep_tracer_fvm_heve( &
       QTRCo, & ! (out)
       qtrc, qtrc0, rhoq_t, &! (in)
       dens0, dens, & ! (in)
       mflx_hi, num_diff, & ! (in)
       gsqrt, mapf, & ! (in)
       cdz, rcdz, rcdx, rcdy, & ! (in)
       dtl, & ! (in)
       flag_fct_tracer, & ! (in)
       flag_fct_along_stream ) ! (in)
    use scale_gridtrans, only: &
       i_xyz, &
       i_xyw, &
       i_uyz, &
       i_xvz, &
       i_xy,  &
       i_uy,  &
       i_xv
    use scale_atmos_dyn_common, only: &
       atmos_dyn_fct
    use scale_atmos_dyn_fvm_flux, only: &
       atmos_dyn_fvm_fluxz_xyz_tracer, &
       atmos_dyn_fvm_fluxx_xyz_tracer, &
       atmos_dyn_fvm_fluxy_xyz_tracer
    use scale_atmos_dyn_fvm_flux_ud1, only: &
       atmos_dyn_fvm_fluxz_xyz_ud1, &
       atmos_dyn_fvm_fluxx_xyz_ud1, &
       atmos_dyn_fvm_fluxy_xyz_ud1
    implicit none
    real(RP), intent(inout) :: qtrco   (ka,ia,ja) ! could be identical to QTRC0
    real(RP), intent(in)  :: qtrc    (ka,ia,ja)
    real(RP), intent(in)  :: qtrc0   (ka,ia,ja)
    real(RP), intent(in)  :: rhoq_t  (ka,ia,ja)
    real(RP), intent(in)  :: dens0   (ka,ia,ja)
    real(RP), intent(in)  :: dens    (ka,ia,ja)
    real(RP), intent(in)  :: mflx_hi (ka,ia,ja,3)
    real(RP), intent(in)  :: num_diff(ka,ia,ja,3)
    real(RP), intent(in)  :: gsqrt   (ka,ia,ja,7)
    real(RP), intent(in)  :: mapf    (ia,ja,2)
    real(RP), intent(in)  :: cdz(ka)
    real(RP), intent(in)  :: rcdz(ka)
    real(RP), intent(in)  :: rcdx(ia)
    real(RP), intent(in)  :: rcdy(ja)
    real(RP), intent(in)  :: dtl
    logical,  intent(in)  :: flag_fct_tracer
    logical,  intent(in)  :: flag_fct_along_stream


    ! For tracer advection
    real(RP) :: qflx_hi  (ka,ia,ja,3)  ! rho * vel(x,y,z) * phi @ (u,v,w)-face high order
    real(RP) :: qflx_lo  (ka,ia,ja,3)  ! rho * vel(x,y,z) * phi,  monotone flux
    real(RP) :: qflx_anti(ka,ia,ja,3)  ! anti-diffusive flux

    integer  :: iis, iie
    integer  :: jjs, jje
    integer  :: i, j, k
    !---------------------------------------------------------------------------

#ifdef DEBUG
    qflx_hi(:,:,:,:) = undef
    qflx_lo(:,:,:,:) = undef
#endif

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

       ! at (x, y, w)
       call atmos_dyn_fvm_fluxz_xyz_tracer( qflx_hi(:,:,:,zdir), & ! (out)
            mflx_hi(:,:,:,zdir), qtrc, gsqrt(:,:,:,i_xyw), & ! (in)
            num_diff(:,:,:,zdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (u, y, z)
       call atmos_dyn_fvm_fluxx_xyz_tracer( qflx_hi(:,:,:,xdir), & ! (out)
            mflx_hi(:,:,:,xdir), qtrc, gsqrt(:,:,:,i_uyz), & ! (in)
            num_diff(:,:,:,xdir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)

       ! at (x, v, z)
       call atmos_dyn_fvm_fluxy_xyz_tracer( qflx_hi(:,:,:,ydir), & ! (out)
            mflx_hi(:,:,:,ydir), qtrc, gsqrt(:,:,:,i_xvz), & ! (in)
            num_diff(:,:,:,ydir), & ! (in)
            cdz, & ! (in)
            iis, iie, jjs, jje ) ! (in)
       
       if ( flag_fct_tracer ) then

          call atmos_dyn_fvm_fluxz_xyz_ud1( qflx_lo(:,:,:,zdir), & ! (out)
               mflx_hi(:,:,:,zdir), qtrc0, gsqrt(:,:,:,i_xyw), & ! (in)
               num_diff(:,:,:,zdir), & ! (in)
               cdz, & ! (in)
               iis-1, iie+1, jjs-1, jje+1 ) ! (in)

          call atmos_dyn_fvm_fluxx_xyz_ud1( qflx_lo(:,:,:,xdir), & ! (out)
               mflx_hi(:,:,:,xdir), qtrc0, gsqrt(:,:,:,i_uyz), & ! (in)
               num_diff(:,:,:,xdir), & ! (in)
               cdz, & ! (in)
               iis-1, iie+1, jjs-1, jje+1 ) ! (in)

          call atmos_dyn_fvm_fluxy_xyz_ud1( qflx_lo(:,:,:,ydir), & ! (out)
               mflx_hi(:,:,:,ydir), qtrc0, gsqrt(:,:,:,i_xvz), & ! (in)
               num_diff(:,:,:,ydir), & ! (in)
               cdz, & ! (in)
               iis-1, iie+1, jjs-1, jje+1 ) ! (in)
       end if

    enddo
    enddo

    if ( flag_fct_tracer ) then

       call atmos_dyn_fct( qflx_anti,            & ! (out)
                           qtrc0, dens0, dens,   & ! (in)
                           qflx_hi, qflx_lo,     & ! (in)
                           mflx_hi,              & ! (in)
                           rcdz, rcdx, rcdy,     & ! (in)
                           gsqrt(:,:,:,i_xyz),   & ! (in)
                           mapf, dtl,            & ! (in)
                           flag_fct_along_stream ) ! (in)

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

          !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
          do j = jjs, jje
          do i = iis, iie
          do k = ks, ke
             qtrco(k,i,j) = ( qtrc0(k,i,j) * dens0(k,i,j) &
                            + dtl * ( - ( ( qflx_hi(k  ,i  ,j  ,zdir) - qflx_anti(k  ,i  ,j  ,zdir) &
                                          - qflx_hi(k-1,i  ,j  ,zdir) + qflx_anti(k-1,i  ,j  ,zdir) ) * rcdz(k) &
                                        + ( qflx_hi(k  ,i  ,j  ,xdir) - qflx_anti(k  ,i  ,j  ,xdir) &
                                          - qflx_hi(k  ,i-1,j  ,xdir) + qflx_anti(k  ,i-1,j  ,xdir) ) * rcdx(i) &
                                        + ( qflx_hi(k  ,i  ,j  ,ydir) - qflx_anti(k  ,i  ,j  ,ydir) &
                                          - qflx_hi(k  ,i  ,j-1,ydir) + qflx_anti(k  ,i  ,j-1,ydir) ) * rcdy(j) &
                                      ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j,i_xyz) &
                               + rhoq_t(k,i,j) ) ) / dens(k,i,j)
          enddo
          enddo
          enddo

       enddo
       enddo

    else ! skip FCT

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

          !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
          !$omp shared(JJS,JJE,IIS,IIE,KS,KE,QTRCo,QTRC0,DENS0,dtl,qflx_hi,RCDZ,RCDX,RCDY,MAPF) &
          !$omp shared(GSQRT,RHOQ_t,DENS,I_XYZ) 
          do j = jjs, jje
          do i = iis, iie
          do k = ks, ke
             qtrco(k,i,j) = ( qtrc0(k,i,j) * dens0(k,i,j) &
                            + dtl * ( - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i  ,j  ,zdir)  ) * rcdz(k) &
                                        + ( qflx_hi(k,i,j,xdir) - qflx_hi(k  ,i-1,j  ,xdir)  ) * rcdx(i) &
                                        + ( qflx_hi(k,i,j,ydir) - qflx_hi(k  ,i  ,j-1,ydir)  ) * rcdy(j) &
                                        ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j,i_xyz) &
                            + rhoq_t(k,i,j) ) ) / dens(k,i,j)
          enddo
          enddo
          enddo

       enddo
       enddo

    end if

    return
  end subroutine atmos_dyn_tstep_tracer_fvm_heve

end module scale_atmos_dyn_tstep_tracer_fvm_heve