scale_atmos_dyn_tstep_tracer_fvm_heve.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_dyn_tstep_tracer_fvm_heve
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_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_prc, only: &
       prc_abort
    implicit none
    character(len=*), intent(in) :: type
 
    if ( type /= 'FVM-HEVE' ) then
       log_error("ATMOS_DYN_Tstep_tracer_fvm_heve_setup",*) 'Tstep_tracer_type is not "FVM-HEVE"!'
       call prc_abort
    end if
 
    return
  end subroutine atmos_dyn_tstep_tracer_fvm_heve_setup
 
  subroutine atmos_dyn_tstep_tracer_fvm_heve( &
       QTRCo, & ! (out)
       qflx_hi,  & ! (out)
       qtrc, qtrc0, rhoq_t, &! (in)
       dens0, dens, & ! (in)
       mflx_hi, num_diff, & ! (in)
       gsqrt, mapf, & ! (in)
       cdz, rcdz, rcdx, rcdy, & ! (in)
       twod, & ! (in)
       dtl, & ! (in)
       flag_fct_tracer, & ! (in)
       flag_fct_along_stream ) ! (in)
    use scale_atmos_dyn_fvm_fct, only: &
       atmos_dyn_fvm_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(out) :: qflx_hi (ka,ia,ja,3) ! rho * vel(x,y,z) * phi @ (u,v,w)-face high order
    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)
    logical,  intent(in)  :: twod
    real(rp), intent(in)  :: dtl
    logical,  intent(in)  :: flag_fct_tracer
    logical,  intent(in)  :: flag_fct_along_stream
 
 
    ! For tracer advection
    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  :: iis0, jjs0
    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)
       if ( .not. twod ) &
       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)
 
          if ( .not. twod ) &
          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_fvm_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, twod, 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-1, ke
             qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) - qflx_anti(k,i,j,zdir)
          end do
          end do
          end do
 
          if ( .not. twod ) then
             if ( iis == is ) then
                iis0 = iis-1
             else
                iis0 = iis
             end if
             !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
             do j = jjs, jje
             do i = iis0, iie
             do k = ks, ke
                qflx_hi(k,i,j,xdir) = qflx_hi(k,i,j,xdir) - qflx_anti(k,i,j,xdir)
             end do
             end do
             end do
          end if
 
          if ( jjs == js ) then
             jjs0 = jjs-1
          else
             jjs0 = jjs
          end if
          !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
          do j = jjs0, jje
          do i = iis, iie
          do k = ks, ke
             qflx_hi(k,i,j,ydir) = qflx_hi(k,i,j,ydir) - qflx_anti(k,i,j,ydir)
          end do
          end do
          end do
 
 
          if ( twod ) then
             !$omp parallel do private(i,j,k) OMP_SCHEDULE_
             do j = jjs, jje
             do k = ks, ke
                qtrco(k,is,j) = ( qtrc0(k,is,j) * dens0(k,is,j) &
                                + dtl * ( - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j  ,zdir)  ) * rcdz(k) &
                                            + ( qflx_hi(k,is,j,ydir) - qflx_hi(k  ,is,j-1,ydir)  ) * rcdy(j) &
                                          ) * mapf(is,j,2) / gsqrt(k,is,j,i_xyz) &
                                + rhoq_t(k,is,j) ) ) / dens(k,is,j)
             enddo
             enddo
          else
             !$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_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
          end if
 
       enddo
       enddo
 
    else ! skip FCT
 
       do jjs = js, je, jblock
       jje = jjs+jblock-1
       do iis = is, ie, iblock
       iie = iis+iblock-1
 
          if ( twod ) then
             !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ &
             !$omp shared(JJS,JJE,IS,KS,KE,QTRCo,QTRC0,DENS0,dtl,qflx_hi,RCDZ,RCDY,MAPF) &
             !$omp shared(GSQRT,RHOQ_t,DENS,I_XYZ) 
             do j = jjs, jje
             do k = ks, ke
                qtrco(k,is,j) = ( qtrc0(k,is,j) * dens0(k,is,j) &
                                + dtl * ( - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j  ,zdir)  ) * rcdz(k) &
                                            + ( qflx_hi(k,is,j,ydir) - qflx_hi(k  ,is,j-1,ydir)  ) * rcdy(j) &
                                           ) * mapf(is,j,2) / gsqrt(k,is,j,i_xyz) &
                               + rhoq_t(k,is,j) ) ) / dens(k,is,j)
             enddo
             enddo
          else
             !$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
          end if
 
       enddo
       enddo
 
    end if
 
    return
  end subroutine atmos_dyn_tstep_tracer_fvm_heve
 
end module scale_atmos_dyn_tstep_tracer_fvm_heve