advec_ppm_algo.f90

!MNH_LIC Copyright 2007-2022 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
!     ##########################
      MODULE MODI_ADVEC_PPM_ALGO
!     ##########################
!
INTERFACE
!
      SUBROUTINE  ADVEC_PPM_ALGO(HMET_ADV_SCHEME, HLBCX, HLBCY, KGRID, PFIELDT,&
                                 PRHODJ, PTSTEP, PTSTEP_PPM,                   &
                                 PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1,PRHOZ2,&
                                 PSRC,  TPDTCUR, PCRU, PCRV, PCRW)
!
USE MODD_TIME, ONLY: DATE_TIME
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
CHARACTER (LEN=6),               INTENT(IN) :: HMET_ADV_SCHEME
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PFIELDT      ! variable at t
INTEGER,                INTENT(IN)  :: KGRID        ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCRU, PCRV, PCRW ! Courant numbers 
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHODJ  ! density
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOX1, PRHOX2
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOY1, PRHOY2
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOZ1, PRHOZ2
REAL,                   INTENT(IN)  :: PTSTEP ! Time step model  
REAL,                   INTENT(IN)  :: PTSTEP_PPM ! Time Step PPM
TYPE (DATE_TIME),       INTENT(IN)  :: TPDTCUR ! current date and time
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRC    ! source term after advection
!
END SUBROUTINE  ADVEC_PPM_ALGO
!
END INTERFACE
!
END MODULE MODI_ADVEC_PPM_ALGO
!
!
!     ##########################################################################
      SUBROUTINE  ADVEC_PPM_ALGO(HMET_ADV_SCHEME, HLBCX, HLBCY, KGRID, PFIELDT, &
                                 PRHODJ, PTSTEP, PTSTEP_PPM, &
                                 PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1,PRHOZ2,&
                                 PSRC, TPDTCUR, PCRU, PCRV, PCRW)
!     ##########################################################################
!!
!!****  *ADVEC_PPM_ALGO* - interface for 3D advection with PPM type scheme
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    MODULE MODD_ARGSLIST
!!         HALO2LIST_ll : type for a list of "HALO2_lls"
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    MODIFICATIONS
!!    -------------
!       10/16 (C.Lac) : Correction on the flag for Strang splitting to insure
!                       reproducibility between START and RESTA
!
!
USE MODD_TYPE_DATE
#ifdef MNH_OPENACC
USE MODE_DEVICE
USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
#endif
use mode_mppdb
#ifdef MNH_OPENACC
use mode_msg
#endif
!
USE MODI_SHUMAN
USE MODI_PPM
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
CHARACTER (LEN=6),               INTENT(IN) :: HMET_ADV_SCHEME
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PFIELDT      ! variable at t
INTEGER,                INTENT(IN)  :: KGRID        ! C grid localisation
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PCRU, PCRV, PCRW ! Courant numbers 
!
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHODJ  ! density
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOX1, PRHOX2
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOY1, PRHOY2
REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRHOZ1, PRHOZ2
REAL,                   INTENT(IN)  :: PTSTEP ! Time step model  
REAL,                   INTENT(IN)  :: PTSTEP_PPM ! Time Step PPM
TYPE (DATE_TIME),       INTENT(IN)  :: TPDTCUR ! current date and time
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRC    ! source term after advection
!
!TYPE(HALO2_ll), OPTIONAL, POINTER   :: TPHALO2 ! halo2 for the field at t
!
!*       0.2   Declarations of local variables :
!
LOGICAL  :: GFLAG   ! Logical flag
#ifdef MNH_OPENACC
REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPPM   ! temp PPM output
#endif
!
!-------------------------------------------------------------------------------
!
!$acc data present( PFIELDT, PCRU, PCRV, PCRW, PRHODJ, PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1, PRHOZ2, PSRC )

IF (MPPDB_INITIALIZED) THEN
  !Check all IN arrays
  CALL MPPDB_CHECK(PFIELDT,"ADVEC_PPM_ALGO beg:PFIELDT")
  CALL MPPDB_CHECK(PCRU,"ADVEC_PPM_ALGO beg:PCRU")
  CALL MPPDB_CHECK(PCRV,"ADVEC_PPM_ALGO beg:PCRV")
  CALL MPPDB_CHECK(PCRW,"ADVEC_PPM_ALGO beg:PCRW")
  CALL MPPDB_CHECK(PRHODJ,"ADVEC_PPM_ALGO beg:PRHODJ")
  CALL MPPDB_CHECK(PRHOX1,"ADVEC_PPM_ALGO beg:PRHOX1")
  CALL MPPDB_CHECK(PRHOX2,"ADVEC_PPM_ALGO beg:PRHOX2")
  CALL MPPDB_CHECK(PRHOY1,"ADVEC_PPM_ALGO beg:PRHOY1")
  CALL MPPDB_CHECK(PRHOY2,"ADVEC_PPM_ALGO beg:PRHOY2")
  CALL MPPDB_CHECK(PRHOZ1,"ADVEC_PPM_ALGO beg:PRHOZ1")
  CALL MPPDB_CHECK(PRHOZ2,"ADVEC_PPM_ALGO beg:PRHOZ2")
END IF
! The scalar PFIELDT is first advected by U*dt first in X, then the resulting
! field is a advected in Y and finally in Z direction. The advection steps are
! stored in PSRC which is finally passed back to the model as a source after
! advection. 
!
!*       0.     INITIAL STEP
!               ------------
!
!$acc kernels present(PSRC,PFIELDT)
PSRC = PFIELDT
!$acc end kernels
GFLAG = ABS(MOD(TPDTCUR%xtime/PTSTEP,2.)-1.) .LE. 0.5
!
SELECT CASE (HMET_ADV_SCHEME)
!
! unlimited scheme (Skamarock notation)
!
CASE('PPM_00')
#ifdef MNH_OPENACC
  call Print_msg( NVERB_WARNING, 'GEN', 'ADVEC_PPM_ALGO', 'OpenACC: PPM_00 not yet tested' )
#endif
!
!  IF (MODULO(KTCOUNT,2) .EQ. 0) THEN ! JUANTEST50
   IF (GFLAG ) THEN 
!
!*       1.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S0_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PTSTEP_PPM)
#else
      CALL PPM_S0_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PTSTEP_PPM, PSRC )
#endif
      !$acc kernels present(PSRC,PRHOX1)
      PSRC = PSRC / PRHOX1
      !$acc end kernels
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S0_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PTSTEP_PPM)
#else
      CALL PPM_S0_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOY1)
      PSRC = PSRC / PRHOY1
      !$acc end kernels
!
!*       3.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S0_Z(KGRID, PSRC, PCRW, PRHOY1, PTSTEP_PPM)
#else
      CALL PPM_S0_Z(KGRID, PSRC, PCRW, PRHOY1, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOZ1)
      PSRC = PSRC / PRHOZ1
      !$acc end kernels
!
   ELSE
!
!
!*       1.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S0_Z(KGRID, PSRC, PCRW, PRHODJ, PTSTEP_PPM)
#else
      CALL PPM_S0_Z(KGRID, PSRC, PCRW, PRHODJ, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOZ2)
      PSRC = PSRC / PRHOZ2
      !$acc end kernels
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S0_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PTSTEP_PPM)
#else
      CALL PPM_S0_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOY2)
      PSRC = PSRC / PRHOY2
      !$acc end kernels
!
!*       3.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
     PSRC = PPM_S0_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PTSTEP_PPM)
#else
      CALL PPM_S0_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOX2)
      PSRC = PSRC / PRHOX2
      !$acc end kernels
!
   END IF
!
! classic (Colella) monotonic scheme
!
CASE('PPM_01')
!
#ifdef MNH_OPENACC
  !Pin positions in the pools of MNH memory
  CALL MNH_MEM_POSITION_PIN()

  CALL MNH_MEM_GET( ZPPM, SIZE( PSRC, 1 ), SIZE( PSRC, 2 ), SIZE( PSRC, 3 ) )

  !$acc data present( ZPPM )

  CALL INIT_ON_HOST_AND_DEVICE(ZPPM,PVALUE=-1e99,HNAME='ADVEC_PPM_ALGO::ZPPM')
#endif
   IF (GFLAG ) THEN
!
!*       1.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHODJ) - &
           PPM_01_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PTSTEP_PPM)
      PSRC = PSRC / PRHOX1
#else
      CALL PPM_01_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHODJ,ZPPM,PRHOX1)
      PSRC = ( PSRC * PRHODJ ) - ZPPM
      PSRC = PSRC / PRHOX1
      !$acc end kernels
#endif
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHOX1) - &
           PPM_01_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PTSTEP_PPM)
      PSRC = PSRC / PRHOY1
#else
      CALL PPM_01_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHOX1,ZPPM,PRHOY1)
      PSRC = (PSRC * PRHOX1) - ZPPM
      PSRC = PSRC / PRHOY1
      !$acc end kernels
#endif
!
!*       3.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHOY1) - &
           PPM_01_Z(KGRID, PSRC, PCRW, PRHOY1, PTSTEP_PPM)
      PSRC = PSRC / PRHOZ1
#else
      CALL PPM_01_Z(KGRID, PSRC, PCRW, PRHOY1, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHOY1,ZPPM,PRHOZ1)
      PSRC = (PSRC * PRHOY1) - ZPPM
      PSRC = PSRC / PRHOZ1
      !$acc end kernels
#endif
!
   ELSE
!
!*       1.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHODJ) - &
           PPM_01_Z(KGRID, PSRC, PCRW, PRHODJ, PTSTEP_PPM)
      PSRC = PSRC / PRHOZ2
#else
      CALL PPM_01_Z(KGRID, PSRC, PCRW, PRHODJ, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHODJ,ZPPM,PRHOZ2)
      PSRC = (PSRC * PRHODJ) - ZPPM
      PSRC = PSRC / PRHOZ2
      !$acc end kernels
#endif
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHOZ2) - &
           PPM_01_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PTSTEP_PPM)
      PSRC = PSRC / PRHOY2
#else
       CALL PPM_01_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHOZ2,ZPPM,PRHOY2)
      PSRC = (PSRC * PRHOZ2) - ZPPM
      PSRC = PSRC / PRHOY2
      !$acc end kernels
#endif
!
!*       3.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = (PSRC * PRHOY2) - &
           PPM_01_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PTSTEP_PPM)
      PSRC = PSRC / PRHOX2
#else
      CALL PPM_01_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PTSTEP_PPM, ZPPM)
      !$acc kernels present(PSRC,PRHOY2,ZPPM,PRHOX2)
      PSRC = (PSRC * PRHOY2) - ZPPM
      PSRC = PSRC / PRHOX2
      !$acc end kernels
#endif
!
   END IF

   !$acc end data

#ifdef MNH_OPENACC
  !Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
  CALL MNH_MEM_RELEASE()
#endif
!
! monotonic scheme (Skamarock notation)
!
CASE('PPM_02')
#ifdef MNH_OPENACC
  call Print_msg( NVERB_WARNING, 'GEN', 'ADVEC_PPM_ALGO', 'OpenACC: PPM_02 not yet tested' )
#endif
!
   IF (GFLAG ) THEN 
!
!*       1.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PRHOX1, PTSTEP_PPM)
#else
      CALL PPM_S1_X(HLBCX, KGRID, PSRC, PCRU, PRHODJ, PRHOX1, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOX1)
      PSRC = PSRC / PRHOX1
      !$acc end kernels
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PRHOY1, PTSTEP_PPM)
#else
      CALL PPM_S1_Y(HLBCY, KGRID, PSRC, PCRV, PRHOX1, PRHOY1, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOY1)
      PSRC = PSRC / PRHOY1
      !$acc end kernels
!
!*       3.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_Z(KGRID, PSRC, PCRW, PRHOY1, PRHOZ1, PTSTEP_PPM)
#else
      CALL PPM_S1_Z(KGRID, PSRC, PCRW, PRHOY1, PRHOZ1, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOZ1)
      PSRC = PSRC / PRHOZ1
      !$acc end kernels
!
   ELSE
!
!*       1.     ADVECTION IN Z DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_Z(KGRID, PSRC, PCRW, PRHODJ, PRHOZ2, PTSTEP_PPM)
#else
      CALL PPM_S1_Z(KGRID, PSRC, PCRW, PRHODJ, PRHOZ2, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOZ2)
      PSRC = PSRC / PRHOZ2
      !$acc end kernels
!
!*       2.     ADVECTION IN Y DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PRHOY2, PTSTEP_PPM)
#else
      CALL PPM_S1_Y(HLBCY, KGRID, PSRC, PCRV, PRHOZ2, PRHOY2, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOY2)
      PSRC = PSRC / PRHOY2
      !$acc end kernels
!
!*       3.     ADVECTION IN X DIRECTION
!               ------------------------
!
#ifndef MNH_OPENACC
      PSRC = PPM_S1_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PRHOX2, PTSTEP_PPM)
#else
      CALL PPM_S1_X(HLBCX, KGRID, PSRC, PCRU, PRHOY2, PRHOX2, PTSTEP_PPM, PSRC)
#endif
      !$acc kernels present(PSRC,PRHOX2)
      PSRC = PSRC / PRHOX2
      !$acc end kernels
!
   END IF
!
END SELECT
!
!*       4.     CREATE THE FORCING TERM TO PASS BACK TO THE MODEL
!               -------------------------------------------------
!
! PSRC now contains the field advected in one dt. 
! To create a forcing term
! compatible to the rest of the model forcings, we need to substract the
! initial field, devide by dt and muliplty by RHODJ
!
!$acc kernels present(PSRC,PFIELDT,PRHODJ)
PSRC = (PSRC - PFIELDT)*PRHODJ/PTSTEP_PPM
!$acc end kernels
!
IF (MPPDB_INITIALIZED) THEN
  !Check all OUT arrays
  CALL MPPDB_CHECK(PSRC,"ADVEC_PPM_ALGO end:PSRC")
END IF
!
!$acc end data

END SUBROUTINE ADVEC_PPM_ALGO