From 1f082a845c9fdf169ef7b936c5240197c109f1e0 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Mon, 27 Mar 2023 16:10:06 +0200
Subject: [PATCH] Philippe 27/03/2023: pressurez: merge ZSOLVER version into
 MNH/

---
 src/MNH/pressurez.f90     |   72 ++-
 src/ZSOLVER/pressurez.f90 | 1160 -------------------------------------
 2 files changed, 61 insertions(+), 1171 deletions(-)
 delete mode 100644 src/ZSOLVER/pressurez.f90

diff --git a/src/MNH/pressurez.f90 b/src/MNH/pressurez.f90
index 8faeb1d9d..5860f6abd 100644
--- a/src/MNH/pressurez.f90
+++ b/src/MNH/pressurez.f90
@@ -288,7 +288,11 @@ USE MODI_SHUMAN
 #ifdef MNH_OPENACC
 USE MODI_SHUMAN_DEVICE
 #endif
-!
+#ifdef MNH_MGSOLVER
+USE MODI_ZCONJGRAD
+USE MODI_ZSOLVER
+#endif
+
 IMPLICIT NONE
 !
 !*       0.1   declarations of arguments
@@ -454,21 +458,21 @@ IF (MPPDB_INITIALIZED) THEN
   CALL MPPDB_CHECK(PDZZ,"PRESSUREZ beg:PDZZ")
   CALL MPPDB_CHECK(PDZX,"PRESSUREZ beg:PDZX")
   CALL MPPDB_CHECK(PRHOT,"PRESSUREZ beg:PRHOT")
-  CALL MPPDB_CHECK(PAF,"PRESSUREZ beg:PAF")
-  CALL MPPDB_CHECK(PCF,"PRESSUREZ beg:PCF")
-  CALL MPPDB_CHECK(PBF,"PRESSUREZ beg:PBF")
-  !CALL MPPDB_CHECK(PTRIGSX,"PRESSUREZ beg:PTRIGSX")
-  !CALL MPPDB_CHECK(PTRIGSY,"PRESSUREZ beg:PTRIGSY")
-  !CALL MPPDB_CHECK(KIFAXX,"PRESSUREZ beg:KIFAXX")
-  !CALL MPPDB_CHECK(KIFAXY,"PRESSUREZ beg:KIFAXY")
+!!$  CALL MPPDB_CHECK(PAF,"PRESSUREZ beg:PAF")
+!!$  CALL MPPDB_CHECK(PCF,"PRESSUREZ beg:PCF")
+!!$  CALL MPPDB_CHECK(PBF,"PRESSUREZ beg:PBF")
+!!$  CALL MPPDB_CHECK(PTRIGSX,"PRESSUREZ beg:PTRIGSX")
+!!$  CALL MPPDB_CHECK(PTRIGSY,"PRESSUREZ beg:PTRIGSY")
+!!$  CALL MPPDB_CHECK(KIFAXX,"PRESSUREZ beg:KIFAXX")
+!!$  CALL MPPDB_CHECK(KIFAXY,"PRESSUREZ beg:KIFAXY")
   CALL MPPDB_CHECK(PTHT,"PRESSUREZ beg:PTHT")
   CALL MPPDB_CHECK(PRT,"PRESSUREZ beg:PRT")
   CALL MPPDB_CHECK(PRHODREF,"PRESSUREZ beg:PRHODREF")
   CALL MPPDB_CHECK(PTHVREF,"PRESSUREZ beg:PTHVREF")
   CALL MPPDB_CHECK(PRVREF,"PRESSUREZ beg:PRVREF")
   CALL MPPDB_CHECK(PEXNREF,"PRESSUREZ beg:PEXNREF")
-  !CALL MPPDB_CHECK(PBFB,"PRESSUREZ beg:PBFB")
-  CALL MPPDB_CHECK(PBF_SXP2_YP1_Z,"PRESSUREZ beg:PBF_SXP2_YP1_Z")
+!!$  CALL MPPDB_CHECK(PBFB,"PRESSUREZ beg:PBFB")
+  !/!\ don't work not of size NX/NY/NZ CALL MPPDB_CHECK(PBF_SXP2_YP1_Z,"PRESSUREZ beg:PBF_SXP2_YP1_Z")
   !Check all INOUT arrays
   CALL MPPDB_CHECK(PRUS,"PRESSUREZ beg:PRUS")
   CALL MPPDB_CHECK(PRVS,"PRESSUREZ beg:PRVS")
@@ -487,8 +491,13 @@ END IF
 ILUOUT = TLUOUT%NLU
 !
 CALL GET_GLOBALDIMS_ll (IIMAX_ll,IJMAX_ll)
+#ifndef MNH_MGSOLVER
 IF ( ( MIN(IIMAX_ll,IJMAX_ll) < NPROC  ) .AND. ( HPRESOPT /= 'ZRESI' ) ) THEN
-   WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TO MANY PROCESSES WITHOUT CPRESOPT="ZRESI" '
+   WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TOO MANY PROCESSES WITHOUT CPRESOPT="ZRESI" '
+#else
+IF ( ( MIN(IIMAX_ll,IJMAX_ll) < NPROC  ) .AND. ( HPRESOPT /= 'ZRESI' ) .AND. ( HPRESOPT /= 'ZSOLV' )) THEN
+   WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TOO MANY PROCESSES WITHOUT CPRESOPT="ZRESI/ZSOLV" '
+#endif
    WRITE(UNIT=ILUOUT,FMT=*) 'MIN(IIMAX_ll,IJMAX_ll)=',MIN(IIMAX_ll,IJMAX_ll),' < NPROC =', NPROC
    WRITE(UNIT=ILUOUT,FMT=*) 'YOU HAVE TO SET CPRESOPT="ZRESI => JOB ABORTED '
    CALL PRINT_MSG(NVERB_FATAL,'GEN','PRESSUREZ','')
@@ -735,12 +744,28 @@ ELSE
      CALL CONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
      PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
      KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
+#ifdef MNH_MGSOLVER
+  CASE('ZGRAD')     ! Conjugate Gradient method
+     CALL ZCONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
+     PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
+     KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
+#endif
 !
   CASE('CRESI')     ! Conjugate Residual method
     CALL CONRESOL(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
     PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
     KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
 !
+#ifdef MNH_MGSOLVER
+  CASE('ZSOLV')     ! Conjugate Residual method
+    CALL ZSOLVER(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
+    PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                      &
+    KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT,                                    &
+    PAF_ZS,PBF_ZS,PCF_ZS,                                                   &
+    PDXATH_ZS,PDYATH_ZS,PRHO_ZS,PBFB,                                       &
+    A_K,B_K,C_K,D_K)
+!
+#endif
   CASE('ZRESI')     ! Conjugate Residual method
     CALL CONRESOLZ(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,      &
     PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
@@ -757,26 +782,51 @@ END IF
 !
 IF ( GWEST ) THEN
    !$acc kernels async
+#if 0
+!PW: original version
     ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB)
     ZPHIT(IIB-1,:,IKE+1) = ZPHIT(IIB,:,IKE)
+#else
+!PW: version from ZSOLVER
+   ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
+#endif
    !$acc end kernels
 END IF
 IF ( GEAST ) THEN
    !$acc kernels async
+#if 0
+!PW: original version
     ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB)
     ZPHIT(IIE+1,:,IKE+1) = ZPHIT(IIE,:,IKE)
+#else
+!PW: version from ZSOLVER
+   ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
+#endif
    !$acc end kernels
 END IF
 !
 IF ( GSOUTH ) THEN
    !$acc kernels async
+#if 0
+!PW: original version
     ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB)
     ZPHIT(:,IJB-1,IKE+1) = ZPHIT(:,IJB,IKE)    
+#else
+!PW: version from ZSOLVER
+   ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
+#endif
    !$acc end kernels
 END IF
 IF ( GNORTH ) THEN
    !$acc kernels async
+#if 0
+!PW: original version
+    ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB)
+    ZPHIT(:,IJE+1,IKE+1) = ZPHIT(:,IJE,IKE)  
+#else
+!PW: version from ZSOLVER
    ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
+#endif
    !$acc end kernels
 END IF
 IF ( GSOUTH2D ) THEN
diff --git a/src/ZSOLVER/pressurez.f90 b/src/ZSOLVER/pressurez.f90
deleted file mode 100644
index b91bb89c2..000000000
--- a/src/ZSOLVER/pressurez.f90
+++ /dev/null
@@ -1,1160 +0,0 @@
-!MNH_LIC Copyright 1994-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_PRESSUREZ
-!###################
-!
-INTERFACE
-!
-      SUBROUTINE PRESSUREZ(                                                &
-      HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI,                &
-      PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOT,               &
-      PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,                           &
-      KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0,                       &
-      PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS,                   &
-      PRUS,PRVS,PRWS,PPABST,                                               &
-      PBFB,                                                                &
-      PBF_SXP2_YP1_Z,                                                      &
-      PAF_ZS,PBF_ZS,PCF_ZS,                                                &
-      PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                                         &
-      A_K,B_K,C_K,D_K,                                                     &
-      PRESIDUAL) !JUAN  FULL ZSOLVER                                        
-!
-IMPLICIT NONE
-!
-CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCX    ! x-direction LBC type
-CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCY    ! y-direction LBC type
-!
-CHARACTER (LEN=5), INTENT(IN) :: HPRESOPT        ! choice of the pressure solver
-!
-INTEGER, INTENT(INOUT) :: KITR                   ! number of iterations for the
-                                                 ! pressure solver
-LOGICAL, INTENT(IN) :: OITRADJ                   ! switch to adjust or not KITR
-INTEGER, INTENT(IN) :: KTCOUNT                   ! counter value of the
-                                                 ! model temporal loop
-INTEGER, INTENT(IN) :: KMI                       ! Model index
-REAL, INTENT(IN)    :: PRELAX                    ! relaxation coefficient for
-                                                 ! the Richardson's method
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ     ! density of reference state
-                                                 ! * J
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX,PDZY ! metric coefficients
-!
-REAL, INTENT(IN) :: PDXHATM                     ! mean grid increment in the x
-                                                ! direction
-REAL, INTENT(IN) :: PDYHATM                     ! mean grid increment in the y
-                                                ! direction
-!
-REAL, DIMENSION (:), INTENT(IN) :: PRHOT         !  mean of XRHODJ on the plane x y
-                                                 !  localized at a mass level
-!
-REAL, DIMENSION(:), INTENT(IN)     :: PAF,PCF    ! vectors giving the nonvanishing
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF        ! elements of the tri-diag. y-slide
-                                                 ! matrix in the pressure eq.
-!
-                                                 ! arrays of sin or cos values
-                                                 ! for the FFT :
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX        ! - along x
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY        ! - along y
-!
-                                                 ! decomposition in prime
-                                                 ! numbers for the FFT:
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX      ! - along x
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY      ! - along y
-!
-!
-INTEGER,                  INTENT(IN)    :: KRR   ! Total number of water var.
-INTEGER,                  INTENT(IN)    :: KRRL  ! Number of liquid water var.
-INTEGER,                  INTENT(IN)    :: KRRI  ! Number of ice water var.
-!
-REAL,                     INTENT(IN)    :: PDRYMASST   ! Mass of dry air and of
-REAL,                     INTENT(IN)    :: PREFMASS    ! the ref. atmosphere
-REAL,                     INTENT(IN)    :: PMASS_O_PHI0 !    Mass / Phi0
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT       ! Temperature and water
-REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT        !  variables at time t
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODREF    ! dry Density
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHVREF     ! Virtual Temperature
-                                                       ! of the reference state
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRVREF      ! mixing ratio of the
-                                                       ! reference state
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF     ! Exner function
-                                                       ! of the reference state
-REAL,                     INTENT(IN)    :: PLINMASS    ! lineic mass through
-                                                       ! open boundaries
-!
-REAL,       INTENT(INOUT) :: PRUS(:,:,:)         ! source term along x
-REAL,       INTENT(INOUT) :: PRVS(:,:,:)         ! source term along y
-REAL,       INTENT(INOUT) :: PRWS(:,:,:)         ! source term along z
-!
-REAL,       INTENT(INOUT)   :: PPABST(:,:,:)        ! pressure(t)
-!
-!JUAN Z_SPLITING
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB       ! elements of the tri-diag b-slide .
-                                                 ! matrix in the pressure eq.
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z       ! elements of the tri-diag. SXP2_YP1_Z-slide 
-                                                   ! matrix in the pressure eq.
-
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(IN) :: A_K,B_K,C_K,D_K
-
-REAL, OPTIONAL                     :: PRESIDUAL
-!JUAN Z_SPLITING
-END SUBROUTINE PRESSUREZ
-!
-END INTERFACE
-!
-END MODULE MODI_PRESSUREZ
-!     ######################################################################
-      SUBROUTINE PRESSUREZ(                                                &
-      HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI,                &
-      PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOT,               &
-      PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,                           &
-      KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0,                       &
-      PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS,                   &
-      PRUS,PRVS,PRWS,PPABST,                                               &
-      PBFB,                                                                &
-      PBF_SXP2_YP1_Z,                                                      &
-      PAF_ZS,PBF_ZS,PCF_ZS,                                                &
-      PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                                         &
-      A_K,B_K,C_K,D_K,                                                     &
-      PRESIDUAL) !JUAN FULL ZSOLVER
-!     ######################################################################
-!
-!!****  *PRESSUREZ * - solve the pressure equation and add the pressure term
-!!      to the sources
-!!
-!!    PURPOSE
-!!    -------
-!       The purpose of this routine is to solve the pressure equation:
-!     with either the conjugate gradient method or the Richardson's method.
-!     The pressure gradient is added to the sources in order
-!     to nullify the divergence of the momentum* Thetavref*(1+Rvref)
-!     at the time t+dt.
-!
-!!**  METHOD
-!!    ------
-!!     The divergence of the sources  ( RHS of the pressure equation ) is
-!!    computed. The pressure equation is then solved by either CG method,
-!!    either Richardson's method, or an exact method. Finally, the pressure
-!!    gradient is added to the sources RUS, RVS, RWS.
-!!    Finally, the absolute pressure is diagnozed from the total mass
-!!    included in the simulation domain.
-!!
-!!    EXTERNAL
-!!    --------
-!!      Subroutine MASS_LEAK : assures global non-divergence condition in the
-!!                             case of open boundaries
-!!      Subroutine FLAT_INV  : solve the pressure equation for the case
-!!                             without orography
-!!      Subroutine RICHARDSON: solve the pressure equation with the
-!!                             Richardson's method
-!!      Subroutine CONJGRAD  : solve the pressure equation with the Conjugate
-!!                             Gradient algorithm
-!!      Function   GX_M_U : compute the gradient along x
-!!      Function   GY_M_V : compute the gradient along y
-!!      Function   GZ_M_W : compute the gradient along z
-!!      Subroutine GDIV     : compute J times the divergence of 1/J times a vector
-!!      Function MXM: compute an average in the x direction for a variable
-!!      at a mass localization
-!!      Function MYM: compute an average in the y direction for a variable
-!!      at a mass localization
-!!      Function MZM: compute an average in the z direction for a variable
-!!      at a mass localization
-!!      Subroutine P_ABS   : compute the constant for PABS and therefore, the
-!!      absolute pressure function
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!      Module MODD_CONF: model configuration
-!!        LFLAT: logical switch for zero orography
-!!        L2D  : logical switch for two-dimensional configuration
-!!        LCARTESIAN : logical switch for cartesian geometry
-!!      Module MODD_PARAMETERS: declaration of parameter variables
-!!        JPHEXT, JPVEXT: define the number of marginal points out of the
-!!        physical domain along horizontal and vertical directions respectively
-!!      Module MODD_CST: physical constants
-!!        XCPD
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation (subroutine PRESSURE) + Book1 (  )
-!!
-!!    AUTHOR
-!!    ------
-!!      P. Hereil and J. Stein       * Meteo France *
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original      05/07/94
-!!      Modification  03/01/95  (Lafore)  To add the absolute pressure diagnosis
-!!      Modification  31/01/95  (Stein)   Copy of the pressure function in the
-!!                                        2D case in the two outermost planes
-!!      Modification  16/02/95  (Mallet)  Add the call to MASS_LEAK
-!!      Modification  16/03/95  (Stein)  change the argument list of the
-!!                              gradient and remove R from the historical var.
-!!      Modification  30/06/95  (Stein)  Add a test not to compute the absolute
-!!                              pressure in the Boussinesq case
-!!                    16/10/95 (J. Stein) change the budget calls
-!!                    29/01/96 (J. Stein) call iterative resolution for
-!!                              non-cartessian geometry
-!!                    19/12/96 (J.-P. Pinty) update the budget calls
-!!                    14/01/97 (Stein,Lafore) New anelastic equations
-!!                    17/12/97 ( Stein )include the case of non-vanishing
-!!                              orography at the lbc
-!!                    26/03/98 (Stein,Jabouille) fix the value of the corner point
-!!                    15/06/98  (D.Lugato, R.Guivarch) Parallelisation
-!!                    25/08/99 (J.-P. Pinty) add CRESI option to CPRESOPT
-!!                    06/11/02 (V. Masson) update the budget calls
-!!                    24/08/2005 (J. escobar) BUG : remove IIE+1, IJE+1 out of bound
-!!                                references in parallel run
-!!                    08/2010 (V.Masson, C.Lac) Add UPDATE_HALO
-!!                    11/2010 (V.Masson, C.Lac) PPABST, must not be cyclic => add temp array
-!!                                             to save it before UPDATE_HALO
-!!                    07/2011 (J.escobar ) Bypass Bug with ifort11/12 on  HLBCX,HLBCY 
-!!                    09/2001 (J.escobar ) reintroduce correctly the GMAXLOC_ll call 
-!!                    11/2010 (V.Masson, C.Lac) PPABST must not be cyclic => add temp array
-!!                                             to save it before UPDATE_HALO
-!!                    02/2013 (J.Escobar ) add a test on abs(err) > 100.O for BG without controle of NAN
-!!                    2012    (V.Masson)  Modif update_halo due to CONTRAV
-!!                    2014    (C.Lac) correction for 3D run with LBOUSS=.TRUE.
-!!   J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
-!!   J.escobar : check nb proc versus ZRESI & min(DIMX,DIMY)
-!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!!  Philippe Wautelet: 22/01/2019: use standard FLUSH statement instead of non standard intrinsics
-!  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
-!  P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
-!  P. Wautelet 28/01/2020: use the new data structures and subroutines for budgets for U
-!!     JL Redelsperger 03/2021 : Shallow convection case added in LHE case: 
-!!                               working for both atmosphere and ocean cases  
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODD_ARGSLIST_ll, ONLY: LIST_ll
-use modd_budget,      only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets
-USE MODD_CST
-USE MODD_CONF
-USE MODD_DYN_n,       ONLY: LRES, XRES,LOCEAN
-USE MODD_FIELD_n,     ONLY: XPHIT
-USE MODD_IBM_PARAM_n, ONLY: LIBM, NIBM_ITR, XIBM_EPSI, XIBM_LS, XIBM_SU
-USE MODD_LUNIT_n,     ONLY: TLUOUT
-USE MODD_MPIF
-USE MODD_PARAMETERS
-use modd_precision,   only: MNHREAL_MPI
-USE MODD_REF,         ONLY: LBOUSS
-USE MODD_VAR_ll,      ONLY: NMNH_COMM_WORLD , NPROC
-!
-use mode_budget,      only: Budget_store_end
-USE MODE_ll
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,   ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-USE MODE_MPPDB
-USE MODE_MSG
-USE MODE_SUM2_ll,     ONLY: GMAXLOC_ll
-!
-#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
-USE MODI_BITREP
-#endif
-USE MODI_CONJGRAD
-USE MODI_CONRESOL
-USE MODI_CONRESOLZ
-USE MODI_FLAT_INV
-USE MODI_FLAT_INVZ
-USE MODI_GDIV
-#ifdef MNH_OPENACC
-USE MODI_GET_HALO
-#endif
-USE MODI_GRADIENT_M
-USE MODI_IBM_BALANCE
-USE MODI_MASS_LEAK
-USE MODI_P_ABS
-USE MODI_RICHARDSON
-USE MODI_SHUMAN
-#ifdef MNH_OPENACC
-USE MODI_SHUMAN_DEVICE
-#endif
-USE MODI_ZCONJGRAD
-USE MODI_ZSOLVER
-
-IMPLICIT NONE
-!
-!*       0.1   declarations of arguments
-!
-  CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCX    ! x-direction LBC type
-  CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCY    ! y-direction LBC type
-!
-CHARACTER (LEN=5), INTENT(IN) :: HPRESOPT        ! choice of the pressure solver
-!
-INTEGER, INTENT(INOUT) :: KITR                   ! number of iterations for the
-                                                 ! pressure solver
-LOGICAL, INTENT(IN) :: OITRADJ                   ! switch to adjust or not KITR
-INTEGER, INTENT(IN) :: KTCOUNT                   ! counter value of the
-                                                 ! model temporal loop
-INTEGER, INTENT(IN) :: KMI                       ! Model index
-REAL, INTENT(IN)    :: PRELAX                    ! relaxation coefficient for
-                                                 ! the Richardson's method
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ     ! density of reference state
-                                                 ! * J
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX,PDZY ! metric coefficients
-!
-REAL, INTENT(IN) :: PDXHATM                     ! mean grid increment in the x
-                                                ! direction
-REAL, INTENT(IN) :: PDYHATM                     ! mean grid increment in the y
-                                                ! direction
-!
-REAL, DIMENSION (:), INTENT(IN) :: PRHOT         !  mean of XRHODJ on the plane x y
-                                                 !  localized at a mass level
-!
-REAL, DIMENSION(:), INTENT(IN)     :: PAF,PCF    ! vectors giving the nonvanishing
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF        ! elements of the tri-diag y-slide .
-                                                 ! matrix in the pressure eq.
-!
-                                                 ! arrays of sin or cos values
-                                                 ! for the FFT :
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX        ! - along x
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY        ! - along y
-!
-                                                 ! decomposition in prime
-                                                 ! numbers for the FFT:
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX      ! - along x
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY      ! - along y
-!
-INTEGER,                  INTENT(IN)    :: KRR   ! Total number of water var.
-INTEGER,                  INTENT(IN)    :: KRRL  ! Number of liquid water var.
-INTEGER,                  INTENT(IN)    :: KRRI  ! Number of ice water var.
-!
-REAL,                     INTENT(IN)    :: PDRYMASST   ! Mass of dry air and of
-REAL,                     INTENT(IN)    :: PREFMASS    ! the ref. atmosphere
-REAL,                     INTENT(IN)    :: PMASS_O_PHI0 !    Mass / Phi0
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT       ! Temperature and water
-REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT        !  variables at time t
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODREF    ! dry Density
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHVREF     ! Virtual Temperature
-                                                       ! of the reference state
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRVREF      ! mixing ratio of the
-                                                       ! reference state
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF     ! Exner function
-                                                       ! of the reference state
-REAL,                     INTENT(IN)    :: PLINMASS    ! lineic mass through
-                                                       ! open boundaries
-!
-REAL,       INTENT(INOUT) :: PRUS(:,:,:)         ! source term along x
-REAL,       INTENT(INOUT) :: PRVS(:,:,:)         ! source term along y
-REAL,       INTENT(INOUT) :: PRWS(:,:,:)         ! source term along z
-!
-REAL,       INTENT(INOUT)   :: PPABST(:,:,:)        ! pressure(t)
-!
-!JUAN Z_SPLITING
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB       ! elements of the tri-diag b-slide .
-                                                 ! matrix in the pressure eq.
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z       ! elements of the tri-diag. SXP2_YP1_Z-slide 
-                                                   ! matrix in the pressure eq.
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(IN) :: A_K,B_K,C_K,D_K
-
-REAL, OPTIONAL                     :: PRESIDUAL
-!JUAN Z_SPLITING
-!
-!*       0.2   declarations of local variables
-!
-!                                                           Metric coefficients:
-!
-REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZDV_SOURCE
-!                                                   ! divergence of the sources
-!
-INTEGER :: IIB          ! indice I for the first inner mass point along x
-INTEGER :: IIE          ! indice I for the last inner mass point along x
-INTEGER :: IJB          ! indice J for the first inner mass point along y
-INTEGER :: IJE          ! indice J for the last inner mass point along y
-INTEGER :: IKB          ! indice K for the first inner mass point along z
-INTEGER :: IKE          ! indice K for the last inner mass point along z
-INTEGER :: ILUOUT       ! Logical unit of output listing
-INTEGER :: IRESP        ! Return code of FM routines
-!
-REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZTHETAV, &
-                        ! virtual potential temperature
-                                                 ZPHIT
-                        ! MAE + DUR => Exner function perturbation
-                        ! LHE       => Exner function perturbation * CPD * THVREF
-!
-REAL :: ZPHI0
-REAL            :: ZRV_OV_RD !  XRV / XRD
-REAL                  :: ZMAXVAL, ZMAXRES, ZMAX,ZMAX_ll ! for print
-INTEGER, DIMENSION(3) :: IMAXLOC ! purpose
-INTEGER         :: JWATER          ! loop index on water species
-INTEGER         :: IIU,IJU,IKU     ! array sizes in I,J,K
-INTEGER         :: JK              ! loop index on the vertical levels
-INTEGER         :: JI,JJ
-!
-REAL, SAVE , POINTER, CONTIGUOUS , DIMENSION(:,:) :: ZPABS_S ! local pressure on southern side
-REAL, SAVE , POINTER, CONTIGUOUS , DIMENSION(:,:) :: ZPABS_N ! local pressure on northern side
-REAL, SAVE , POINTER, CONTIGUOUS , DIMENSION(:,:) :: ZPABS_E ! local pressure on eastern side
-REAL, SAVE , POINTER, CONTIGUOUS , DIMENSION(:,:) :: ZPABS_W ! local pressure on western side
-INTEGER :: IINFO_ll,KINFO
-TYPE(LIST_ll), POINTER :: TZFIELDS_ll, TZFIELDS2_ll  ! list of fields to exchange
-!
-INTEGER :: IIB_I,IIE_I,IJB_I,IJE_I
-INTEGER :: IIMAX_ll,IJMAX_ll
-!
-#ifdef MNH_OPENACC
-REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZPRHODJ,ZMXM_PRHODJ,ZMZM_PRHODJ,ZGZ_M_W
-REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZMYM_PRHODJ
-#endif
-!
-!
-LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH
-LOGICAL :: GSOUTH2D,GNORTH2D,GPRVREF0
-!
-LOGICAL, SAVE :: GFIRST_CALL_PRESSUREZ = .TRUE.
-!
-!------------------------------------------------------------------------------
-IF (MPPDB_INITIALIZED) THEN
-  !Check all IN arrays
-  CALL MPPDB_CHECK(PRHODJ,"PRESSUREZ beg:PRHODJ")
-  CALL MPPDB_CHECK(PDXX,"PRESSUREZ beg:PDXX")
-  CALL MPPDB_CHECK(PDYY,"PRESSUREZ beg:PDYY")
-  CALL MPPDB_CHECK(PDZZ,"PRESSUREZ beg:PDZZ")
-  CALL MPPDB_CHECK(PDZX,"PRESSUREZ beg:PDZX")
-  CALL MPPDB_CHECK(PRHOT,"PRESSUREZ beg:PRHOT")
-!!$  CALL MPPDB_CHECK(PAF,"PRESSUREZ beg:PAF")
-!!$  CALL MPPDB_CHECK(PCF,"PRESSUREZ beg:PCF")
-!!$  CALL MPPDB_CHECK(PBF,"PRESSUREZ beg:PBF")
-!!$  CALL MPPDB_CHECK(PTRIGSX,"PRESSUREZ beg:PTRIGSX")
-!!$  CALL MPPDB_CHECK(PTRIGSY,"PRESSUREZ beg:PTRIGSY")
-!!$  CALL MPPDB_CHECK(KIFAXX,"PRESSUREZ beg:KIFAXX")
-!!$  CALL MPPDB_CHECK(KIFAXY,"PRESSUREZ beg:KIFAXY")
-  CALL MPPDB_CHECK(PTHT,"PRESSUREZ beg:PTHT")
-  CALL MPPDB_CHECK(PRT,"PRESSUREZ beg:PRT")
-  CALL MPPDB_CHECK(PRHODREF,"PRESSUREZ beg:PRHODREF")
-  CALL MPPDB_CHECK(PTHVREF,"PRESSUREZ beg:PTHVREF")
-  CALL MPPDB_CHECK(PRVREF,"PRESSUREZ beg:PRVREF")
-  CALL MPPDB_CHECK(PEXNREF,"PRESSUREZ beg:PEXNREF")
-!!$  CALL MPPDB_CHECK(PBFB,"PRESSUREZ beg:PBFB")
-  !/!\ don't work not of size NX/NY/NZ CALL MPPDB_CHECK(PBF_SXP2_YP1_Z,"PRESSUREZ beg:PBF_SXP2_YP1_Z")
-  !Check all INOUT arrays
-  CALL MPPDB_CHECK(PRUS,"PRESSUREZ beg:PRUS")
-  CALL MPPDB_CHECK(PRVS,"PRESSUREZ beg:PRVS")
-  CALL MPPDB_CHECK(PRWS,"PRESSUREZ beg:PRWS")
-  CALL MPPDB_CHECK(PPABST,"PRESSUREZ beg:PPABST")
-END IF
-!-------------------------------------------------------------------------------
-!
-!*       1.      PRELIMINARIES
-!                -------------
-!
-ILUOUT = TLUOUT%NLU
-!
-CALL GET_GLOBALDIMS_ll (IIMAX_ll,IJMAX_ll)
-IF ( ( MIN(IIMAX_ll,IJMAX_ll) < NPROC  ) .AND. ( HPRESOPT /= 'ZRESI' ) .AND. ( HPRESOPT /= 'ZSOLV' )) THEN
-   WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TO MANY PROCESSES WITHOUT CPRESOPT="ZRESI/ZSOLV" '
-   WRITE(UNIT=ILUOUT,FMT=*) 'MIN(IIMAX_ll,IJMAX_ll)=',MIN(IIMAX_ll,IJMAX_ll),' < NPROC =', NPROC
-   WRITE(UNIT=ILUOUT,FMT=*) 'YOU HAVE TO SET CPRESOPT="ZRESI => JOB ABORTED '
-   CALL PRINT_MSG(NVERB_FATAL,'GEN','PRESSUREZ','')
-ENDIF
-CALL GET_PHYSICAL_ll(IIB,IJB,IIE,IJE)
-CALL GET_DIM_EXT_ll('B',IIU,IJU)
-!
-IKB= 1+JPVEXT
-IKU= SIZE(PPABST,3)
-IKE= IKU - JPVEXT
-!
-GWEST  = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() )
-GEAST  = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() )
-GSOUTH = ( .NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() )
-GNORTH = ( .NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() )
-GSOUTH2D = ( L2D .AND. LSOUTH_ll() )
-GNORTH2D = ( L2D .AND. LNORTH_ll() )
-!
-GPRVREF0 =  ( SIZE(PRVREF,1) == 0 )
-!
-#ifndef MNH_OPENACC
-ALLOCATE( ZDV_SOURCE(IIU, IJU, IKU) )
-ALLOCATE( ZTHETAV   (IIU, IJU, IKU) )
-ALLOCATE( ZPHIT     (IIU, IJU, IKU) )
-ALLOCATE( ZPABS_S(IIU,IKU),ZPABS_N(IIU,IKU),ZPABS_E(IJU,IKU),ZPABS_W(IJU,IKU))
-#else
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( ZDV_SOURCE, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZTHETAV,    IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZPHIT,      IIU, IJU, IKU )
-!
-CALL MNH_MEM_GET( ZPABS_S, IIU,IKU )
-CALL MNH_MEM_GET( ZPABS_N, IIU,IKU )
-CALL MNH_MEM_GET( ZPABS_E, IJU,IKU )
-CALL MNH_MEM_GET( ZPABS_W, IJU,IKU )
-!
-CALL MNH_MEM_GET( ZPRHODJ,     IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZMXM_PRHODJ, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZMZM_PRHODJ, IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZGZ_M_W,     IIU, IJU, IKU )
-CALL MNH_MEM_GET( ZMYM_PRHODJ, IIU, IJU, IKU )
-#endif
-!
-!$acc kernels 
-ZPABS_S(:,:) = 0.
-ZPABS_N(:,:) = 0.
-ZPABS_E(:,:) = 0.
-ZPABS_W(:,:) = 0.
-!$acc end kernels
-
-! Done in model_n before call to Rad_bound
-! if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'PRES', prus(:, :, :) )
-! if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'PRES', prvs(:, :, :) )
-! if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'PRES', prws(:, :, :) )
-
-!-------------------------------------------------------------------------------
-!
-!*       3.    COMPUTE THE LINEIC MASS
-!              -----------------------
-!
-IF ( ANY(HLBCX(:)=='OPEN') .OR. ANY(HLBCY(:)=='OPEN') ) THEN                     
-  CALL MASS_LEAK(PDXX,PDYY,HLBCX,HLBCY,PLINMASS,PRHODJ,PRUS,PRVS)
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!*       4.    COMPUTE THE FORCING TERM FOR THE PRESSURE EQUATION
-!              --------------------------------------------------
-!
-IF (LIBM) THEN
-  WHERE(XIBM_LS(:,:,:,2).GT.-XIBM_EPSI)
-         !dir$ safe_address
-          PRUS(:,:,:) = 0.
-  ENDWHERE
-  WHERE(XIBM_LS(:,:,:,3).GT.-XIBM_EPSI) 
-          !dir$ safe_address
-          PRVS(:,:,:) = 0.
-  ENDWHERE
-  WHERE(XIBM_LS(:,:,:,4).GT.-XIBM_EPSI)
-         !dir$ safe_address
-          PRWS(:,:,:) = 0.
-  ENDWHERE
-ENDIF
-!
-IF (MPPDB_INITIALIZED) THEN
-   CALL MPPDB_CHECK3D(PRUS,"pressurez 4-before update_halo_ll::PRUS",PRECISION)
-   CALL MPPDB_CHECK3D(PRVS,"pressurez 4-before update_halo_ll::PRVS",PRECISION)
-   CALL MPPDB_CHECK3D(PRWS,"pressurez 4-before update_halo_ll::PRWS",PRECISION)
-END IF
-#ifndef MNH_OPENACC
-NULLIFY(TZFIELDS_ll)
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRUS, 'PRESSUREZ::PRUS' )
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRVS, 'PRESSUREZ::PRVS' )
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRWS, 'PRESSUREZ::PRWS' )
-CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
-CALL CLEANLIST_ll(TZFIELDS_ll)
-#else
-CALL GET_HALO_D(PRUS,HNAME='PRESSUREZ::PRUS' )
-CALL GET_HALO_D(PRVS,HNAME='PRESSUREZ::PRVS' )
-CALL GET_HALO_D(PRWS,HNAME='PRESSUREZ::PRWS' )
-#endif
-IF (MPPDB_INITIALIZED) THEN
-   CALL MPPDB_CHECK3D(PRUS,"pressurez 4-after update_halo_ll::PRUS",PRECISION)
-   CALL MPPDB_CHECK3D(PRVS,"pressurez 4-after update_halo_ll::PRVS",PRECISION)
-   CALL MPPDB_CHECK3D(PRWS,"pressurez 4-after update_halo_ll::PRWS",PRECISION)
-END IF
-!
-CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
-!
-IF (LIBM) THEN
-   CALL IBM_BALANCE(XIBM_LS,XIBM_SU,PRUS,PRVS,PRWS,ZDV_SOURCE)
-ENDIF  
-!
-! The non-homogenous Neuman problem is transformed in an homogenous Neuman
-! problem in the non-periodic cases
-IF ( GWEST ) THEN
-   !$acc kernels async
-   ZDV_SOURCE(IIB-1,:,:) = 0.
-   !$acc end kernels
-END IF
-IF ( GEAST ) THEN
-   !$acc kernels async
-   ZDV_SOURCE(IIE+1,:,:) = 0.
-   !$acc end kernels
-END IF
-IF ( GSOUTH ) THEN
-   !$acc kernels async
-   ZDV_SOURCE(:,IJB-1,:) = 0.
-   !$acc end kernels
-END IF
-IF ( GNORTH ) THEN
-   !$acc kernels async
-   ZDV_SOURCE(:,IJE+1,:) = 0.
-   !$acc end kernels
-END IF
-!$acc wait
-
-IF (LIBM) THEN  
-  !
-  IF (HLBCX(1) == 'CYCL') THEN
-    IF (LWEST_ll()) ZDV_SOURCE(IIB-1,:,:) = ZDV_SOURCE(IIB-1,:,:)*XIBM_SU(IIB,:,:,1)
-    IF (LEAST_ll()) ZDV_SOURCE(IIE+1,:,:) = ZDV_SOURCE(IIE+1,:,:)*XIBM_SU(IIE,:,:,1)
-  ENDIF
-  !
-  IF (HLBCY(1) == 'CYCL') THEN
-    IF (LSOUTH_ll()) ZDV_SOURCE(:,IJB-1,:) = ZDV_SOURCE(:,IJB-1,:)*XIBM_SU(:,IJB,:,1)
-    IF (LNORTH_ll()) ZDV_SOURCE(:,IJE+1,:) = ZDV_SOURCE(:,IJE+1,:)*XIBM_SU(:,IJE,:,1)
-  ENDIF
-  !
-  ZDV_SOURCE(:,:,IKB-1) = ZDV_SOURCE(:,:,IKB-1)*XIBM_SU(:,:,IKB,1)
-  ZDV_SOURCE(:,:,IKE+1) = ZDV_SOURCE(:,:,IKE+1)*XIBM_SU(:,:,IKE,1)
-  !
-ENDIF
-!
-!-------------------------------------------------------------------------------
-!
-!*       5.    SOLVE THE PRESSURE EQUATION
-!              ---------------------------
-!
-!
-!*       5.1   Compute the virtual theta and the pressure perturbation
-!              -------------------------------------------------------
-!
-IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
-  !$acc kernels 
-  IF(KRR > 0) THEN
-  !
-  !   compute the ratio : 1 + total water mass / dry air mass
-    ZRV_OV_RD = XRV / XRD
-    ZTHETAV(:,:,:) = 1. + PRT(:,:,:,1)
-    !$acc loop seq
-    DO JWATER = 2 , 1+KRRL+KRRI
-      ZTHETAV(:,:,:) = ZTHETAV(:,:,:) + PRT(:,:,:,JWATER)
-    END DO
-  !   compute the virtual potential temperature when water is present in any
-  !   form
-    ZTHETAV(:,:,:) = PTHT(:,:,:) * (1. + PRT(:,:,:,1) * ZRV_OV_RD) / ZTHETAV(:,:,:)
-  ELSE
-  !   compute the virtual potential temperature when water is absent
-    ZTHETAV(:,:,:) = PTHT(:,:,:)
-  END IF
-  !$acc end kernels
-  !
-  IF (LIBM) THEN 
-    WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI)
-      ZTHETAV(:,:,:) = PTHVREF(:,:,:)
-    ENDWHERE 
-  ENDIF
-  !
-#ifndef MNH_OPENACC   
-#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
-  ZPHIT(:,:,:)=(PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
-#else
-  ZPHIT(:,:,:)=BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD)-PEXNREF(:,:,:)
-#endif
-#else
-  !$acc kernels
-  DO CONCURRENT ( JI=1:IIU,JJ=1:IJU,JK=1:IKU )
-#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
-    ZPHIT(JI,JJ,JK)=(PPABST(JI,JJ,JK)/XP00)**(XRD/XCPD)-PEXNREF(JI,JJ,JK)
-#else
-    ZPHIT(JI,JJ,JK)=BR_POW((PPABST(JI,JJ,JK)/XP00),(XRD/XCPD))-PEXNREF(JI,JJ,JK)
-#endif
-  END DO
-  !$acc end kernels
-#endif  
-  !
-ELSEIF(CEQNSYS=='LHE') THEN
-  IF ( .NOT. LOCEAN) THEN
-    ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:))   &
-                 * XCPD * PTHVREF(:,:,:)
-  ELSE
-    ! Field at T- DT for LHE anelastic approx
-    ! not used in ocean case (flat LHE)
-    ZPHIT(:,:,:)=0.
-  END IF
-!
-END IF
-!
-IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN .AND. .NOT. LIBM) THEN
-   ! flat cartesian LHE case -> exact solution
- IF ( HPRESOPT /= "ZRESI" ) THEN
-  CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,         &
-               PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZDV_SOURCE,ZPHIT)
- ELSE
-  CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,         &
-               PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZDV_SOURCE,ZPHIT,& 
-               PBFB,&
-               PBF_SXP2_YP1_Z)
- ENDIF
-ELSE
-  SELECT CASE(HPRESOPT)
-  CASE('RICHA')     ! Richardson's method
-!
-    CALL RICHARDSON(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,      &
-    PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                        &
-    KIFAXX,KIFAXY,KITR,KTCOUNT,PRELAX,ZDV_SOURCE,ZPHIT)
-!
-   CASE('CGRAD')     ! Conjugate Gradient method
-     CALL CONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
-     PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
-     KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
-  CASE('ZGRAD')     ! Conjugate Gradient method
-     CALL ZCONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
-     PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
-     KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
-!
-  CASE('CRESI')     ! Conjugate Residual method
-    CALL CONRESOL(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
-    PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
-    KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
-!
-  CASE('ZSOLV')     ! Conjugate Residual method
-    CALL ZSOLVER(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
-    PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                      &
-    KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT,                                    &
-    PAF_ZS,PBF_ZS,PCF_ZS,                                                   &
-    PDXATH_ZS,PDYATH_ZS,PRHO_ZS,PBFB,                                       &
-    A_K,B_K,C_K,D_K)
-!
-  CASE('ZRESI')     ! Conjugate Residual method
-!$acc update host( ZDV_SOURCE, ZPHIT, ZTHETAV )
-    CALL CONRESOLZ(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV,       &
-    PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF,PTRIGSX,PTRIGSY,                       &
-    KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT,                                     &
-    PBFB,&
-    PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
-!$acc update device( ZPHIT )
-  END SELECT
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!*       6.    ADD THE PRESSURE GRADIENT TO THE SOURCES
-!              ----------------------------------------
-!
-IF (MPPDB_INITIALIZED) THEN
-   CALL MPPDB_CHECK(ZPHIT,"PRESSUREZ after zsolv :ZPHIT")
-END IF
-!
-IF ( GWEST ) THEN
-   !$acc kernels async
-   ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
-   !$acc end kernels
-END IF
-IF ( GEAST ) THEN
-   !$acc kernels async
-   ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
-   !$acc end kernels
-END IF
-IF ( GSOUTH ) THEN
-   !$acc kernels async
-   ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
-   !$acc end kernels
-END IF
-IF ( GNORTH ) THEN
-   !$acc kernels async
-   ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
-   !$acc end kernels
-END IF
-IF ( GSOUTH2D ) THEN
-   !$acc kernels async
-   ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)
-   !$acc end kernels
-END IF
-IF ( GNORTH2D ) THEN
-   !$acc kernels async
-   ZPHIT(:,IJE+1,:) = ZPHIT(:,IJB,:)
-   !$acc end kernels
-END IF
-!$acc wait
-!
-IF (LIBM) THEN
-  !
-  IF ( HLBCX(1) == 'CYCL' ) THEN
-     IF (LWEST_ll()) THEN
-       ZPHIT(IIB-1,:,:) = ZPHIT(IIB,:,:)*(1.-XIBM_SU(IIB,:,:,1))+XIBM_SU(IIB,:,:,1)*ZPHIT(IIB-1,:,:)
-     ENDIF    
-     IF (LEAST_ll()) THEN
-       ZPHIT(IIE+1,:,:) = ZPHIT(IIE,:,:)*(1.-XIBM_SU(IIE,:,:,1))+XIBM_SU(IIE,:,:,1)*ZPHIT(IIE+1,:,:)
-     ENDIF
-  ENDIF
-  !
-  IF ( HLBCY(1) == 'CYCL' ) THEN
-    IF (LSOUTH_ll()) THEN
-      ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)*(1.-XIBM_SU(:,IJB,:,1))+XIBM_SU(:,IJB,:,1)*ZPHIT(:,IJB-1,:)    
-    ENDIF
-    IF (LNORTH_ll()) THEN
-      ZPHIT(:,IJE+1,:) = ZPHIT(:,IJE,:)*(1.-XIBM_SU(:,IJE,:,1))+XIBM_SU(:,IJE,:,1)*ZPHIT(:,IJE+1,:) 
-    ENDIF  
-  ENDIF
-  !
-  !-------------Bottom Boundary conditions
-  ZPHIT(:,:,IKB-1) = ZPHIT(:,:,IKB-1)*XIBM_SU(:,:,IKB,1)+(1.-XIBM_SU(:,:,IKB,1))*ZPHIT(:,:,IKB)
-  ZPHIT(:,:,IKE+1) = ZPHIT(:,:,IKE+1)*XIBM_SU(:,:,IKE,1)+(1.-XIBM_SU(:,:,IKE,1))*ZPHIT(:,:,IKE) 
-  !
-ENDIF
-!
-#ifndef MNH_OPENACC
-ZDV_SOURCE = GX_M_U(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX)
-#else
-CALL GX_M_U_DEVICE(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX,ZDV_SOURCE)
-#endif
-!
-IF ( GWEST ) THEN
-!!!!!!!!!!!!!!!!  FUJI  compiler directive !!!!!!!!!!
-!!!!!!!!!!!!!!!!  FUJI  compiler directive !!!!!!!!!!
-   !$acc kernels async
-#ifdef MNH_COMPILER_NVHPC
-   !$acc loop independent collapse(2)
-#else
-   !$acc loop independent
-#endif
-   DO CONCURRENT (JJ=1:IJU , JK=2:IKU-1)
-      ZDV_SOURCE(IIB,JJ,JK)=                                                    &
-       (ZPHIT(IIB,JJ,JK) - ZPHIT(IIB-1,JJ,JK) - 0.5 * (                              &
-        PDZX(IIB,JJ,JK)   * (ZPHIT(IIB,JJ,JK)-ZPHIT(IIB,JJ,JK-1)) / PDZZ(IIB,JJ,JK)      &
-       +PDZX(IIB,JJ,JK+1) * (ZPHIT(IIB,JJ,JK+1)-ZPHIT(IIB,JJ,JK)) / PDZZ(IIB,JJ,JK+1)    &
-                                              )                              &
-       ) / PDXX(IIB,JJ,JK)
-   END DO
-   !$acc end kernels
-ENDIF
-  !
-IF( GEAST ) THEN
-   !$acc kernels async
-#ifdef MNH_COMPILER_NVHPC
-   !$acc loop independent collapse(2)
-#else
-   !$acc loop independent
-#endif
-   DO CONCURRENT (JJ=1:IJU , JK=2:IKU-1)
-      ZDV_SOURCE(IIE+1,JJ,JK)=                                                   &
-        (ZPHIT(IIE+1,JJ,JK) - ZPHIT(IIE+1-1,JJ,JK) - 0.5 * (                        &
-         PDZX(IIE+1,JJ,JK)   * (ZPHIT(IIE+1-1,JJ,JK)-ZPHIT(IIE+1-1,JJ,JK-1))           &
-                          / PDZZ(IIE+1-1,JJ,JK)                                  &
-        +PDZX(IIE+1,JJ,JK+1) * (ZPHIT(IIE+1-1,JJ,JK+1)-ZPHIT(IIE+1-1,JJ,JK))           &
-                          / PDZZ(IIE+1-1,JJ,JK+1)                                &
-                                                     )                        &
-        ) / PDXX(IIE+1,JJ,JK)
-   END DO
-   !$acc end kernels
-END IF
-!$acc wait
-!
-IF (MPPDB_INITIALIZED) THEN
-   CALL MPPDB_CHECK3DM("before MXM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS)
-END IF
-IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
-#ifndef MNH_OPENACC
-  PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
-  PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
-#else
-  !$acc kernels present(ZPRHODJ)
-  ZPRHODJ(:,:,:) = PRHODJ(:,:,:) * XCPD * ZTHETAV(:,:,:)
-  !$acc end kernels 
-  CALL MXM_DEVICE(ZPRHODJ, ZMXM_PRHODJ)
-  CALL MZM_DEVICE(ZPRHODJ, ZMZM_PRHODJ)
-  CALL GZ_M_W_DEVICE(1,IKU,1,ZPHIT,PDZZ,ZGZ_M_W)
-  !$acc kernels
-  !dir$ concurrent
-  PRUS(:,:,:) = PRUS(:,:,:) - ZMXM_PRHODJ(:,:,:) * ZDV_SOURCE(:,:,:)
-  !dir$ concurrent
-  PRWS(:,:,:) = PRWS(:,:,:) - ZMZM_PRHODJ(:,:,:) * ZGZ_M_W(:,:,:)
-  !$acc end kernels
-#endif
-ELSEIF(CEQNSYS=='LHE') THEN
-  PRUS = PRUS - MXM(PRHODJ) * ZDV_SOURCE
-  PRWS = PRWS - MZM(PRHODJ) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
-END IF
-!
-IF(.NOT. L2D) THEN
-!
-#ifndef MNH_OPENACC    
-   ZDV_SOURCE = GY_M_V(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY)
-#else
-   CALL GY_M_V_DEVICE(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY,ZDV_SOURCE)
-#endif
-!
-   IF ( GSOUTH ) THEN
-!!!!!!!!!!!!!!!!  FUJI  compiler directive !!!!!!!!!!
-!!!!!!!!!!!!!!!!  FUJI  compiler directive !!!!!!!!!!
-      !$acc kernels async
-#ifdef MNH_COMPILER_NVHPC
-      !$acc loop independent collapse(2)
-#else
-      !$acc loop independent
-#endif
-      DO CONCURRENT (JI=1:IIU , JK=2:IKU-1)
-         ZDV_SOURCE(JI,IJB,JK)=                                                  &
-              (ZPHIT(JI,IJB,JK) - ZPHIT(JI,IJB-1,JK) - 0.5 * (                            &
-              PDZY(JI,IJB,JK)   * (ZPHIT(JI,IJB,JK)-ZPHIT(JI,IJB,JK-1)) / PDZZ(JI,IJB,JK)    &
-              +PDZY(JI,IJB,JK+1) * (ZPHIT(JI,IJB,JK+1)-ZPHIT(JI,IJB,JK)) / PDZZ(JI,IJB,JK+1)  &
-              )                            &
-              ) / PDYY(JI,IJB,JK)
-      END DO
-      !$acc end kernels
-   END IF
-   !
-   IF ( GNORTH ) THEN
-      !$acc kernels async
-#ifdef MNH_COMPILER_NVHPC
-      !$acc loop independent collapse(2)
-#else
-      !$acc loop independent
-#endif
-      DO CONCURRENT (JI=1:IIU , JK=2:IKU-1) 
-         ZDV_SOURCE(JI,IJE+1,JK)=                                                &
-              (ZPHIT(JI,IJE+1,JK) - ZPHIT(JI,IJE+1-1,JK) - 0.5 * (                      &
-              PDZY(JI,IJE+1,JK)   * (ZPHIT(JI,IJE+1-1,JK)-ZPHIT(JI,IJE+1-1,JK-1))         &
-              / PDZZ(JI,IJE+1-1,JK)                                &
-              +PDZY(JI,IJE+1,JK+1) * (ZPHIT(JI,IJE+1-1,JK+1)-ZPHIT(JI,IJE+1-1,JK))         &
-              / PDZZ(JI,IJE+1-1,JK+1)                              &
-              )                      &
-              ) / PDYY(JI,IJE+1,JK)
-      END DO
-      !$acc end kernels
-  END IF
-!$acc wait  
-!
-IF (MPPDB_INITIALIZED) THEN  
-   CALL MPPDB_CHECK3DM("before MYM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS)
-END IF
-  IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
-#ifndef MNH_OPENACC
-    PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
-#else
-     CALL MYM_DEVICE(ZPRHODJ,ZMYM_PRHODJ)
-     !$acc kernels
-     !dir$ concurrent
-     PRVS(:,:,:) = PRVS(:,:,:) - ZMYM_PRHODJ(:,:,:) * ZDV_SOURCE(:,:,:)
-     !$acc end kernels
-#endif
-  ELSEIF(CEQNSYS=='LHE') THEN
-    PRVS = PRVS - MYM(PRHODJ) * ZDV_SOURCE
-  END IF
-END IF
-!
-!! same boundary conditions as in gdiv ... !! (provisory coding)
-!! (necessary when NVERB=1)
-!!
-    !$acc kernels
-    PRUS(:,:,IKB-1)=PRUS(:,:,IKB)
-    PRUS(:,:,IKE+1)=PRUS(:,:,IKE)
-    PRVS(:,:,IKB-1)=PRVS(:,:,IKB)
-    PRVS(:,:,IKE+1)=PRVS(:,:,IKE)
-    !$acc end kernels
-!
-#ifndef MNH_OPENACC    
-NULLIFY(TZFIELDS2_ll)
-CALL ADD3DFIELD_ll( TZFIELDS2_ll, PRUS, 'PRESSUREZ::PRUS' )
-CALL ADD3DFIELD_ll( TZFIELDS2_ll, PRVS, 'PRESSUREZ::PRVS' )
-CALL ADD3DFIELD_ll( TZFIELDS2_ll, PRWS, 'PRESSUREZ::PRWS' )
-CALL UPDATE_HALO_ll(TZFIELDS2_ll,IINFO_ll)
-CALL CLEANLIST_ll(TZFIELDS2_ll)
-#else
-CALL GET_HALO_D(PRUS,HNAME='PRESSUREZ::PRUS' )
-CALL GET_HALO_D(PRVS,HNAME='PRESSUREZ::PRVS' )
-CALL GET_HALO_D(PRWS,HNAME='PRESSUREZ::PRWS' )
-#endif
-!
-!  compute the residual divergence
-CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
-!
-IF (LIBM) THEN
-  ZDV_SOURCE(:,:,:)=ZDV_SOURCE(:,:,:)*XIBM_SU(:,:,:,2)
-ENDIF  
-!
-IF ( CEQNSYS=='DUR' ) THEN
-  IF ( GPRVREF0 ) THEN
-    !$acc kernels  
-    ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF
-    !$acc end kernels  
-  ELSE
-    !$acc kernels
-    ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF*(1.+PRVREF)
-    !$acc end kernels
-  END IF
-ELSEIF( CEQNSYS=='MAE' .OR. CEQNSYS=='LHE' ) THEN
-  ZDV_SOURCE=ZDV_SOURCE/PRHODJ*PRHODREF
-END IF
-!
-!$acc update host(ZDV_SOURCE)
-ZMAXVAL=MAX_ll(ABS(ZDV_SOURCE),IINFO_ll)
-!JUANZ
-IF (PRESENT(PRESIDUAL)) PRESIDUAL = ZMAXVAL
-!JUANZ
-IMAXLOC=GMAXLOC_ll( ABS(ZDV_SOURCE) )
-!
-WRITE(ILUOUT,*) 'residual divergence / 2 DT', ZMAXVAL,     &
-                ' located at ',   IMAXLOC
-FLUSH(unit=ILUOUT)
-IF (ABS(ZMAXVAL) .GT. 100.0 ) THEN
-   call Print_msg( NVERB_FATAL, 'GEN', 'PRESSUREZ', 'something wrong with pressure: abs(residual) > 100.0' )
-END IF
-! number of iterations adjusted
-IF (LRES) THEN
-   ZMAXRES = XRES
-ELSEIF (LFLAT .AND. LCARTESIAN) THEN
-  ZMAXRES = XRES_FLAT_CART
-ELSE
-  ZMAXRES = XRES_OTHER
-END IF
-!
-IF (OITRADJ) THEN
-  IF (ZMAXVAL>10.*ZMAXRES) THEN
-    KITR=KITR+2
-    WRITE(ILUOUT,*) 'NITR adjusted to ', KITR
-  ELSE IF (ZMAXVAL<ZMAXRES) THEN
-    KITR=MAX(KITR-1,1)
-    WRITE(ILUOUT,*) 'NITR adjusted to ', KITR
-  ENDIF
-ENDIF
-!
-IF (LIBM) THEN
-  KITR=MIN(NIBM_ITR,KITR)
-ENDIF
-!
-!*       7.    STORAGE OF THE FIELDS IN BUDGET ARRAYS
-!              --------------------------------------
-!
-if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'PRES', prus(:, :, :) )
-if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'PRES', prvs(:, :, :) )
-if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'PRES', prws(:, :, :) )
-!
-!-------------------------------------------------------------------------------
-!
-!*       8.    ABSOLUTE PRESSURE COMPUTATION
-!              -----------------------------
-!
-ZMAX = MAXVAL(ABS ( PRHODREF(:,:,IKB)-PRHODREF(:,:,IKE)) )
-CALL MPI_ALLREDUCE(ZMAX, ZMAX_ll, 1, MNHREAL_MPI, MPI_MAX,  &
-                   NMNH_COMM_WORLD, KINFO)
-!IF (      ABS(PRHODREF(IIB,IJB,IKB)-PRHODREF(IIB,IJB,IKE)) > 1.E-12 &
-!  .AND. KTCOUNT >0 ) THEN
-IF ((ZMAX_ll > 1.E-12) .AND. KTCOUNT >0 ) THEN 
-!IF (  KTCOUNT >0 .AND. .NOT.LBOUSS ) THEN
-  CALL P_ABS   ( KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
-                 PTHT, PRT, PRHODJ, PRHODREF, ZTHETAV, PTHVREF,      &
-                 PRVREF, PEXNREF, ZPHIT, ZPHI0                       )
-!
-  IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
-     !$acc kernels
-#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
-     PPABST(:,:,:)=XP00*(ZPHIT(:,:,:)+PEXNREF(:,:,:))**(XCPD/XRD)
-#else
-     DO CONCURRENT(JI=1:IIU,JJ=1:IJU,JK=1:IKU)
-        PPABST(JI,JJ,JK)=XP00*BR_POW((ZPHIT(JI,JJ,JK)+PEXNREF(JI,JJ,JK)),(XCPD/XRD))
-     END DO
-#endif
-     !$acc end kernels
-  ELSEIF(CEQNSYS=='LHE') THEN
-    IF (.NOT. LOCEAN) THEN
-      ! Deep atmosphere case : computing of PI fluctuation ; ZPHI0 (computed in P_ABS routine) is added 
-      XPHIT(:,:,:) = (ZPHIT+ZPHI0)/(XCPD*PTHVREF)
-      ! Absolute Pressure 
-      PPABST(:,:,:)=XP00*(XPHIT(:,:,:)+PEXNREF)**(XCPD/XRD)
-      ! Computing press fluctuation
-      XPHIT(:,:,:) = PPABST(:,:,:) - XP00*PEXNREF**(XCPD/XRD)
-    ELSE
-!    Shallow atmosphere ou ocean
-     XPHIT(:,:,:) = (ZPHIT+ZPHI0)*PRHODREF
-     PPABST(:,:,:)=XPHIT(:,:,:) + XP00*PEXNREF**(XCPD/XRD)
-    END IF
-  ENDIF
-!
-  IF( GWEST ) THEN
-     !$acc kernels async
-     ZPABS_W(:,:)= PPABST(IIB,:,:)
-     !$acc end kernels
-  END IF
-!
-  IF ( GEAST ) THEN
-     !$acc kernels async
-     ZPABS_E(:,:)= PPABST(IIE+1,:,:)
-     !$acc end kernels
-  END IF
-!
-  IF( GSOUTH ) THEN
-     !$acc kernels async
-     ZPABS_S(:,:)= PPABST(:,IJB,:)
-     !$acc end kernels
-  END IF
-!
-  IF ( GNORTH ) THEN
-     !$acc kernels async
-     ZPABS_N(:,:)= PPABST(:,IJE+1,:)
-     !$acc end kernels
-  END IF
-  !
-  !$acc wait
-  !
-#ifndef MNH_OPENACC  
-  CALL ADD3DFIELD_ll( TZFIELDS_ll, PPABST, 'PRESSUREZ::PPABST' )
-  CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
-  CALL CLEANLIST_ll(TZFIELDS_ll)
-#else
-  CALL GET_HALO_D( PPABST,HNAME='PRESSUREZ::PPABST' )
-#endif
-IF (MPPDB_INITIALIZED) THEN  
-   CALL MPPDB_CHECK(PPABST,"PRESSUREZ after get halo_d :PPABST")
-END IF
-!
-  IF( GWEST ) THEN
-     !$acc kernels async
-     PPABST(IIB,:,:) = ZPABS_W(:,:)
-     !$acc end kernels
-  END IF
-!
-  IF ( GEAST ) THEN
-     !$acc kernels async
-     PPABST(IIE+1,:,:) = ZPABS_E(:,:)
-     !$acc end kernels
-  END IF
-!
-  IF( GSOUTH ) THEN
-     !$acc kernels async
-     PPABST(:,IJB,:) = ZPABS_S(:,:)
-     !$acc end kernels
-  END IF
-!
-  IF ( GNORTH ) THEN
-     !$acc kernels async
-     PPABST(:,IJE+1,:) = ZPABS_N(:,:)
-     !$acc end kernels
-  END IF
-!
-!$acc wait  
-!
-END IF
-
-IF (MPPDB_INITIALIZED) THEN
-  !Check all INOUT arrays
-  CALL MPPDB_CHECK(PRUS,"PRESSUREZ end:PRUS")
-  CALL MPPDB_CHECK(PRVS,"PRESSUREZ end:PRVS")
-  CALL MPPDB_CHECK(PRWS,"PRESSUREZ end:PRWS")
-  CALL MPPDB_CHECK(PPABST,"PRESSUREZ end:PPABST")
-END IF
-!
-#ifndef MNH_OPENACC
-deallocate( ZDV_SOURCE )
-deallocate( ZTHETAV )
-deallocate( ZPHIT )
-deallocate( ZPABS_S,ZPABS_N,ZPABS_E,ZPABS_W )
-#else
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE PRESSUREZ
-- 
GitLab