diff --git a/src/ZSOLVER/advection_metsv.f90 b/src/ZSOLVER/advection_metsv.f90
new file mode 100644
index 0000000000000000000000000000000000000000..4a69e687f32445631e359ff30b6794ddf9cda828
--- /dev/null
+++ b/src/ZSOLVER/advection_metsv.f90
@@ -0,0 +1,1126 @@
+!MNH_LIC Copyright 1994-2019 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_ADVECTION_METSV
+! ###########################
+!
+INTERFACE
+ SUBROUTINE ADVECTION_METSV (TPFILE, OCLOSE_OUT,HUVW_ADV_SCHEME, &
+ HMET_ADV_SCHEME,HSV_ADV_SCHEME, HCLOUD, KSPLIT, &
+ OSPLIT_CFL, PSPLIT_CFL, OCFL_WRIT, &
+ HLBCX, HLBCY, KRR, KSV, TPDTCUR, PTSTEP, &
+ PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT, PPABST, &
+ PTHVREF, PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRTHS, PRRS, PRTKES, PRSVS, &
+ PRTHS_CLD, PRRS_CLD, PRSVS_CLD, PRTKES_ADV )
+!
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_TYPE_DATE, ONLY: DATE_TIME
+!
+LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous
+ ! file opening
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+CHARACTER(LEN=6), INTENT(IN) :: HMET_ADV_SCHEME, & ! Control of the
+ HSV_ADV_SCHEME, & ! scheme applied
+ HUVW_ADV_SCHEME
+CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of cloud parameterization
+!
+INTEGER, INTENT(INOUT):: KSPLIT ! Number of time splitting
+ ! for PPM advection
+LOGICAL, INTENT(IN) :: OSPLIT_CFL ! flag to automatically chose number of iterations
+REAL, INTENT(IN) :: PSPLIT_CFL ! maximum CFL to automatically chose number of iterations
+LOGICAL, INTENT(IN) :: OCFL_WRIT ! flag to write CFL fields in output files
+!
+CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
+!
+INTEGER, INTENT(IN) :: KRR ! Number of moist variables
+INTEGER, INTENT(IN) :: KSV ! Number of Scalar Variables
+!
+TYPE (DATE_TIME), INTENT(IN) :: TPDTCUR ! current date and time
+REAL, INTENT(IN) :: PTSTEP
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT , PVT , PWT
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET, PRHODJ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT , PSVT
+ ! Variables at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
+ ! of the reference state
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
+ ! metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS, PRTKES
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS
+ ! Sources terms
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRTHS_CLD
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRRS_CLD,PRSVS_CLD
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRTKES_ADV ! Advection TKE source term
+!
+END SUBROUTINE ADVECTION_METSV
+!
+END INTERFACE
+!
+END MODULE MODI_ADVECTION_METSV
+! ##########################################################################
+ SUBROUTINE ADVECTION_METSV (TPFILE, OCLOSE_OUT,HUVW_ADV_SCHEME, &
+ HMET_ADV_SCHEME,HSV_ADV_SCHEME, HCLOUD, KSPLIT, &
+ OSPLIT_CFL, PSPLIT_CFL, OCFL_WRIT, &
+ HLBCX, HLBCY, KRR, KSV, TPDTCUR, PTSTEP, &
+ PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT, PPABST, &
+ PTHVREF, PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRTHS, PRRS, PRTKES, PRSVS, &
+ PRTHS_CLD, PRRS_CLD, PRSVS_CLD, PRTKES_ADV )
+! ##########################################################################
+!
+!!**** *ADVECTION_METSV * - routine to call the specialized advection routines
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to control the advection routines.
+!! For that, it is first necessary to compute the metric coefficients
+!! and the contravariant components of the momentum.
+!!
+!!** METHOD
+!! ------
+!! Once the scheme is selected, it is applied to the following group of
+!! variables: METeorologicals (temperature, water substances, TKE,
+!! dissipation TKE) and Scalar Variables. It is possible to select different
+!! advection schemes for each group of variables.
+!!
+!! EXTERNAL
+!! --------
+!! CONTRAV : computes the contravariant components.
+!! ADVECUVW : computes the advection terms for momentum.
+!! ADVECSCALAR : computes the advection terms for scalar fields.
+!! ADD3DFIELD_ll : add a field to 3D-list
+!! ADVEC_4TH_ORDER : 4th order advection scheme
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book1 and book2 ( routine ADVECTION )
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!! J.-P. Lafore * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 06/07/94
+!! 01/04/95 (Ph. Hereil J. Nicolau) add the model number
+!! 23/10/95 (J. Vila and JP Lafore) advection schemes scalar
+!! 16/01/97 (JP Pinty) change presentation
+!! 30/04/98 (J. Stein P Jabouille) extrapolation for the cyclic
+!! case and parallelisation
+!! 24/06/99 (P Jabouille) case of NHALO>1
+!! 25/10/05 (JP Pinty) 4th order scheme
+!! 24/04/06 (C.Lac) Split scalar and passive
+!! tracer routines
+!! 08/06 (T.Maric) PPM scheme
+!! 04/2011 (V.Masson & C. Lac) splits the routine and add time splitting
+!! 04/2014 (C.Lac) adaptation of time
+!! splitting for L1D and L2D
+!! 09/2014 (G.Delautier) close OUTPUT_LISTING before STOP
+!! 04/2015 (J.Escobar) remove/commente some NHALO=1 test
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! J.Escobar : 01/10/2015 : add computation of CFL for L1D case
+!! 04/2016 (C.Lac) : correction of negativity for KHKO
+!! 10/2016 (C.Lac) Correction on the flag for Strang splitting
+!! to insure reproducibility between START and RESTA
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 07/2017 (V. Vionnet) : add advection of 2D variables at
+!! the surface for the blowing snow scheme
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+#ifdef MNH_OPENACC
+USE MODD_MPIF
+use modd_precision, only: MNHREAL_MPI
+USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
+#endif
+USE MODD_BUDGET
+USE MODD_CST
+USE MODD_CTURB, ONLY: XTKEMIN
+USE MODD_CONF, ONLY: LNEUTRAL,NHALO,L1D, L2D
+USE MODD_IO, ONLY: TFILEDATA
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODD_PARAM_n
+USE MODD_TYPE_DATE, ONLY: DATE_TIME
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
+USE MODD_PARAMETERS
+!
+use mode_argslist_ll, only: ADD3DFIELD_ll, ADD4DFIELD_ll, CLEANLIST_ll, LIST_ll
+use mode_exchange_ll, only: UPDATE_HALO_ll
+USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
+USE MODE_IO_FIELD_WRITE, only: IO_Field_write
+use mode_mppdb
+USE MODE_MSG
+use mode_sum_ll, only: MAX_ll
+use mode_tools_ll, only: GET_INDICE_ll, lnorth_ll, lsouth_ll, least_ll, lwest_ll
+!
+USE MODI_ADV_BOUNDARIES
+USE MODI_BUDGET
+USE MODI_CONTRAV
+USE MODI_GET_HALO
+USE MODI_PPM_RHODJ
+USE MODI_PPM_MET
+USE MODI_PPM_SCALAR
+!
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+USE MODE_MNH_ZWORK, ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D, MNH_ALLOCATE_ZT3DP , MNH_ALLOCATE_ZT2D, &
+ MNH_ALLOCATE_ZT4D , MNH_REL_ZT4D, &
+ MNH_CHECK_IN_ZT3D,MNH_CHECK_OUT_ZT3D
+USE MODI_GET_HALO
+#endif
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
+!
+!-------------------------------------------------------------------------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for synchronous
+ ! file opening
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
+CHARACTER(LEN=6), INTENT(IN) :: HMET_ADV_SCHEME, & ! Control of the
+ HSV_ADV_SCHEME, & ! scheme applied
+ HUVW_ADV_SCHEME
+CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of cloud parameterization
+!
+INTEGER, INTENT(INOUT):: KSPLIT ! Number of time splitting
+ ! for PPM advection
+LOGICAL, INTENT(IN) :: OSPLIT_CFL ! flag to automatically chose number of iterations
+REAL, INTENT(IN) :: PSPLIT_CFL ! maximum CFL to automatically chose number of iterations
+LOGICAL, INTENT(IN) :: OCFL_WRIT ! flag to write CFL fields in output files
+!
+CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
+!
+INTEGER, INTENT(IN) :: KRR ! Number of moist variables
+INTEGER, INTENT(IN) :: KSV ! Number of Scalar Variables
+!
+TYPE (DATE_TIME), INTENT(IN) :: TPDTCUR ! current date and time
+REAL, INTENT(IN) :: PTSTEP
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT , PVT , PWT
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET, PRHODJ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT , PSVT
+ ! Variables at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
+ ! of the reference state
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
+ ! metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS, PRTKES
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS
+ ! Sources terms
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRTHS_CLD
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRRS_CLD, PRSVS_CLD
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRTKES_ADV ! Advection TKE source term
+!
+!
+!* 0.2 declarations of local variables
+!
+!
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRUCPPM
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRVCPPM
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRWCPPM
+ ! contravariant
+ ! components
+ ! of momentum
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFLU
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFLV
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFLW
+! ! CFL numbers on each direction
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFL
+! ! CFL number
+INTEGER :: IZRUCPPM,IZRVCPPM,IZRWCPPM,IZCFLU,IZCFLV,IZCFLW,IZCFL
+!
+REAL :: ZCFLU_MAX, ZCFLV_MAX, ZCFLW_MAX, ZCFL_MAX ! maximum CFL numbers
+!
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZTH
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZTKE
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTHS_OTHER
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTKES_OTHER
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTHS_PPM
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTKES_PPM
+INTEGER :: IZTH,IZTKE,IZRTHS_OTHER,IZRTKES_OTHER,IZRTHS_PPM,IZRTKES_PPM
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZR
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSV
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSNWC
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSNWC_INIT
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS
+INTEGER :: IZR,IZSV,IZSNWC,IZSNWC_INIT,IZRSNWCS
+! Guess at the sub time step
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRRS_OTHER
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSVS_OTHER
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS_OTHER
+INTEGER :: IZRRS_OTHER,IZRSVS_OTHER,IZRSNWCS_OTHER
+! Tendencies since the beginning of the time step
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRRS_PPM
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSVS_PPM
+REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS_PPM
+INTEGER :: IZRRS_PPM,IZRSVS_PPM,IZRSNWCS_PPM
+! Guess at the end of the sub time step
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOX1,ZRHOX2
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOY1,ZRHOY2
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOZ1,ZRHOZ2
+REAL, DIMENSION(:,:,:),pointer , contiguous :: ZT,ZEXN,ZLV,ZLS,ZCPH
+INTEGER :: IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2 &
+ ,IZT,IZEXN,IZLV,IZLS,IZCPH
+
+! Temporary advected rhodj for PPM routines
+!
+INTEGER :: JS,JR,JSV,JSPL, JI, JJ ! Loop index
+REAL :: ZTSTEP_PPM ! Sub Time step
+LOGICAL :: GTKE
+!
+INTEGER :: IINFO_ll ! return code of parallel routine
+TYPE(LIST_ll), POINTER :: TZFIELDS0_ll ! list of fields to exchange
+TYPE(LIST_ll), POINTER :: TZFIELDS1_ll ! list of fields to exchange
+!
+!
+INTEGER :: IRESP ! Return code of FM routines
+INTEGER :: ILUOUT ! logical unit
+INTEGER :: ISPLIT_PPM ! temporal time splitting
+INTEGER :: IIB, IIE, IJB, IJE,IKB,IKE
+#ifdef MNH_OPENACC
+INTEGER :: IZ1, IZ2
+#endif
+TYPE(TFIELDDATA) :: TZFIELD
+!
+INTEGER :: JIU,JJU,JKU
+INTEGER :: JK
+!-------------------------------------------------------------------------------
+!$acc data present( PUT, PVT, PWT, PTHT, PTKET, PRHODJ, PPABST, PRT, PSVT, PTHVREF, &
+!$acc & PDXX, PDYY, PDZZ, PDZX, PDZY, PRTHS, PRTKES, PRRS, PRSVS, PRTHS_CLD, PRRS_CLD, PRSVS_CLD, PRTKES_ADV )
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all IN arrays
+ CALL MPPDB_CHECK(PUT,"ADVECTION_METSV beg:PUT")
+ CALL MPPDB_CHECK(PVT,"ADVECTION_METSV beg:PVT")
+ CALL MPPDB_CHECK(PWT,"ADVECTION_METSV beg:PWT")
+ CALL MPPDB_CHECK(PTHT,"ADVECTION_METSV beg:PTHT")
+ CALL MPPDB_CHECK(PTKET,"ADVECTION_METSV beg:PTKET")
+ CALL MPPDB_CHECK(PRHODJ,"ADVECTION_METSV beg:PRHODJ")
+ CALL MPPDB_CHECK(PPABST,"ADVECTION_METSV beg:PPABST")
+ CALL MPPDB_CHECK(PRT,"ADVECTION_METSV beg:PRT")
+ CALL MPPDB_CHECK(PSVT,"ADVECTION_METSV beg:PSVT")
+ CALL MPPDB_CHECK(PTHVREF,"ADVECTION_METSV beg:PTHVREF")
+ CALL MPPDB_CHECK(PDXX,"ADVECTION_METSV beg:PDXX")
+ CALL MPPDB_CHECK(PDYY,"ADVECTION_METSV beg:PDYY")
+ CALL MPPDB_CHECK(PDZZ,"ADVECTION_METSV beg:PDZZ")
+ CALL MPPDB_CHECK(PDZX,"ADVECTION_METSV beg:PDZX")
+ CALL MPPDB_CHECK(PDZY,"ADVECTION_METSV beg:PDZY")
+ CALL MPPDB_CHECK(PRTHS_CLD,"ADVECTION_METSV beg:PRTHS_CLD")
+ CALL MPPDB_CHECK(PRRS_CLD,"ADVECTION_METSV beg:PRRS_CLD")
+ CALL MPPDB_CHECK(PRSVS_CLD,"ADVECTION_METSV beg:PRSVS_CLD")
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRTHS,"ADVECTION_METSV beg:PRTHS")
+ CALL MPPDB_CHECK(PRTKES,"ADVECTION_METSV beg:PRTKES")
+ CALL MPPDB_CHECK(PRRS,"ADVECTION_METSV beg:PRRS")
+ CALL MPPDB_CHECK(PRSVS,"ADVECTION_METSV beg:PRSVS")
+END IF
+
+JIU = size(PUT, 1 )
+JJU = size(PUT, 2 )
+JKU = size(PUT, 3 )
+
+
+#ifndef MNH_OPENACC
+allocate( ZRUCPPM ( JIU,JJU,JKU ) )
+allocate( ZRVCPPM ( JIU,JJU,JKU ) )
+allocate( ZRWCPPM ( JIU,JJU,JKU ) )
+allocate( ZCFLU ( JIU,JJU,JKU ) )
+allocate( ZCFLV ( JIU,JJU,JKU ) )
+allocate( ZCFLW ( JIU,JJU,JKU ) )
+allocate( ZCFL ( JIU,JJU,JKU ) )
+allocate( ZTH ( JIU,JJU,JKU ) )
+allocate( ZTKE ( JIU,JJU,SIZE(PTKET,3)) )
+allocate( ZRTHS_OTHER ( JIU,JJU,JKU ) )
+allocate( ZRTKES_OTHER ( JIU,JJU,SIZE(PTKET,3)) )
+allocate( ZRTHS_PPM ( JIU,JJU,JKU ) )
+allocate( ZRTKES_PPM ( JIU,JJU,SIZE(PTKET,3)) )
+allocate( ZR ( JIU,JJU,JKU, SIZE(PRT, 4) ) )
+allocate( ZSV ( JIU,JJU,JKU, SIZE(PSVT,4) ) )
+allocate( ZSNWC ( JIU,JJU,JKU, NBLOWSNOW_2D ) )
+allocate( ZSNWC_INIT ( JIU,JJU,JKU, NBLOWSNOW_2D ) )
+allocate( ZRSNWCS ( JIU,JJU,JKU, NBLOWSNOW_2D ) )
+allocate( ZRRS_OTHER ( JIU,JJU,JKU, SIZE(PRT, 4) ) )
+allocate( ZRSVS_OTHER ( JIU,JJU,JKU, SIZE(PSVT,4) ) )
+allocate( ZRSNWCS_OTHER( JIU,JJU,JKU, NBLOWSNOW_2D ) )
+allocate( ZRRS_PPM ( JIU,JJU,JKU, SIZE(PRT, 4) ) )
+allocate( ZRSVS_PPM ( JIU,JJU,JKU, SIZE(PSVT,4) ) )
+allocate( ZRSNWCS_PPM ( JIU,JJU,JKU, NBLOWSNOW_2D ) )
+allocate( ZRHOX1 ( JIU,JJU,JKU ) )
+allocate( ZRHOX2 ( JIU,JJU,JKU ) )
+allocate( ZRHOY1 ( JIU,JJU,JKU ) )
+allocate( ZRHOY2 ( JIU,JJU,JKU ) )
+allocate( ZRHOZ1 ( JIU,JJU,JKU ) )
+allocate( ZRHOZ2 ( JIU,JJU,JKU ) )
+allocate( ZT ( JIU,JJU,JKU ) )
+allocate( ZEXN ( JIU,JJU,JKU ) )
+allocate( ZLV ( JIU,JJU,JKU ) )
+allocate( ZLS ( JIU,JJU,JKU ) )
+allocate( ZCPH ( JIU,JJU,JKU ) )
+#else
+CALL MNH_CHECK_IN_ZT3D("ADVECTION_METSV")
+IZRUCPPM = MNH_ALLOCATE_ZT3D( ZRUCPPM , JIU,JJU,JKU )
+IZRVCPPM = MNH_ALLOCATE_ZT3D( ZRVCPPM , JIU,JJU,JKU )
+IZRWCPPM = MNH_ALLOCATE_ZT3D( ZRWCPPM , JIU,JJU,JKU )
+IZCFLU = MNH_ALLOCATE_ZT3D( ZCFLU , JIU,JJU,JKU )
+IZCFLV = MNH_ALLOCATE_ZT3D( ZCFLV , JIU,JJU,JKU )
+IZCFLW = MNH_ALLOCATE_ZT3D( ZCFLW , JIU,JJU,JKU )
+IZCFL = MNH_ALLOCATE_ZT3D( ZCFL , JIU,JJU,JKU )
+IZTH = MNH_ALLOCATE_ZT3D( ZTH , JIU,JJU,JKU )
+IZTKE = MNH_ALLOCATE_ZT3D( ZTKE , JIU,JJU,SIZE(PTKET,3) )
+IZRTHS_OTHER = MNH_ALLOCATE_ZT3D( ZRTHS_OTHER , JIU,JJU,JKU )
+IZRTKES_OTHER = MNH_ALLOCATE_ZT3D( ZRTKES_OTHER , JIU,JJU,SIZE(PTKET,3) )
+IZRTHS_PPM = MNH_ALLOCATE_ZT3D( ZRTHS_PPM , JIU,JJU,JKU )
+IZRTKES_PPM = MNH_ALLOCATE_ZT3D( ZRTKES_PPM , JIU,JJU,SIZE(PTKET,3) )
+IZR = MNH_ALLOCATE_ZT4D( ZR , JIU,JJU,JKU, SIZE(PRT, 4) )
+IZSV = MNH_ALLOCATE_ZT4D( ZSV , JIU,JJU,JKU, SIZE(PSVT,4) )
+IZSNWC = MNH_ALLOCATE_ZT4D( ZSNWC , JIU,JJU,JKU, NBLOWSNOW_2D )
+IZSNWC_INIT = MNH_ALLOCATE_ZT4D( ZSNWC_INIT , JIU,JJU,JKU, NBLOWSNOW_2D )
+IZRSNWCS = MNH_ALLOCATE_ZT4D( ZRSNWCS , JIU,JJU,JKU, NBLOWSNOW_2D )
+IZRRS_OTHER = MNH_ALLOCATE_ZT4D( ZRRS_OTHER , JIU,JJU,JKU, SIZE(PRT, 4) )
+IZRSVS_OTHER = MNH_ALLOCATE_ZT4D( ZRSVS_OTHER , JIU,JJU,JKU, SIZE(PSVT,4) )
+IZRSNWCS_OTHER = MNH_ALLOCATE_ZT4D( ZRSNWCS_OTHER, JIU,JJU,JKU, NBLOWSNOW_2D )
+IZRRS_PPM = MNH_ALLOCATE_ZT4D( ZRRS_PPM , JIU,JJU,JKU, SIZE(PRT, 4) )
+IZRSVS_PPM = MNH_ALLOCATE_ZT4D( ZRSVS_PPM , JIU,JJU,JKU, SIZE(PSVT,4) )
+IZRSNWCS_PPM = MNH_ALLOCATE_ZT4D( ZRSNWCS_PPM , JIU,JJU,JKU, NBLOWSNOW_2D )
+IZRHOX1 = MNH_ALLOCATE_ZT3D( ZRHOX1 , JIU,JJU,JKU )
+IZRHOX2 = MNH_ALLOCATE_ZT3D( ZRHOX2 , JIU,JJU,JKU )
+IZRHOY1 = MNH_ALLOCATE_ZT3D( ZRHOY1 , JIU,JJU,JKU )
+IZRHOY2 = MNH_ALLOCATE_ZT3D( ZRHOY2 , JIU,JJU,JKU )
+IZRHOZ1 = MNH_ALLOCATE_ZT3D( ZRHOZ1 , JIU,JJU,JKU )
+IZRHOZ2 = MNH_ALLOCATE_ZT3D( ZRHOZ2 , JIU,JJU,JKU )
+IZT = MNH_ALLOCATE_ZT3D( ZT , JIU,JJU,JKU )
+IZEXN = MNH_ALLOCATE_ZT3D( ZEXN , JIU,JJU,JKU )
+IZLV = MNH_ALLOCATE_ZT3D( ZLV , JIU,JJU,JKU )
+IZLS = MNH_ALLOCATE_ZT3D( ZLS , JIU,JJU,JKU )
+IZCPH = MNH_ALLOCATE_ZT3D( ZCPH , JIU,JJU,JKU )
+#endif
+
+!$acc data present( ZRUCPPM, ZRVCPPM, ZRWCPPM, ZCFLU, ZCFLV, ZCFLW, ZCFL, ZTH, &
+!$acc & ZTKE, ZRTHS_OTHER, ZRTKES_OTHER, ZRTHS_PPM, ZRTKES_PPM, &
+!$acc & ZR, ZSV, ZSNWC, ZSNWC_INIT, ZRSNWCS, ZRRS_OTHER, ZRSVS_OTHER, ZRSNWCS_OTHER, &
+!$acc & ZRRS_PPM, ZRSVS_PPM, ZRSNWCS_PPM, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+!$acc & ZT, ZEXN, ZLV, ZLS, ZCPH )
+
+!
+!* 0. INITIALIZATION
+! --------------
+!
+ILUOUT = TLUOUT%NLU
+!
+CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+IKB=1+JPVEXT
+IKE=SIZE(PSVT,3) - JPVEXT
+
+!
+GTKE=(SIZE(PTKET)/=0)
+!
+#ifdef MNH_OPENACC
+CALL INIT_ON_HOST_AND_DEVICE(PRTKES_ADV,PVALUE=-1e99,HNAME='ADVECTION_METSV::PRTKES_ADV')
+CALL INIT_ON_HOST_AND_DEVICE(ZRUCPPM,PVALUE=-1e90,HNAME='ADVECTION_METSV::ZRUCPPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZRVCPPM,PVALUE=-1e91,HNAME='ADVECTION_METSV::ZRVCPPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZRWCPPM,PVALUE=-1e92,HNAME='ADVECTION_METSV::ZRWCPPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZCFLU,PVALUE=-1e99,HNAME='ADVECTION_METSV::ZCFLU')
+CALL INIT_ON_HOST_AND_DEVICE(ZCFLV,PVALUE=-1e99,HNAME='ADVECTION_METSV::ZCFLV')
+CALL INIT_ON_HOST_AND_DEVICE(ZCFLW,PVALUE=-1e99,HNAME='ADVECTION_METSV::ZCFLW')
+CALL INIT_ON_HOST_AND_DEVICE(ZCFL, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZCFL')
+CALL INIT_ON_HOST_AND_DEVICE(ZTH, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZTH')
+CALL INIT_ON_HOST_AND_DEVICE(ZTKE, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZTKE')
+CALL INIT_ON_HOST_AND_DEVICE(ZRTHS_OTHER, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRTHS_OTHER')
+CALL INIT_ON_HOST_AND_DEVICE(ZRTKES_OTHER,PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRTKES_OTHER')
+CALL INIT_ON_HOST_AND_DEVICE(ZRTHS_PPM, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRTHS_PPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZRTKES_PPM, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRTKES_PPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZR, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZR')
+CALL INIT_ON_HOST_AND_DEVICE(ZSV, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZSV')
+CALL INIT_ON_HOST_AND_DEVICE(ZRRS_OTHER, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRRS_OTHER')
+CALL INIT_ON_HOST_AND_DEVICE(ZRSVS_OTHER, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRSVS_OTHER')
+CALL INIT_ON_HOST_AND_DEVICE(ZRRS_PPM, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRRS_PPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZRSVS_PPM, PVALUE=-1e99,HNAME='ADVECTION_METSV::ZRSVS_PPM')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOX1,PVALUE=-1e93,HNAME='ADVECTION_METSV::ZRHOX1')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOX2,PVALUE=-1e94,HNAME='ADVECTION_METSV::ZRHOX2')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOY1,PVALUE=-1e95,HNAME='ADVECTION_METSV::ZRHOY1')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOY2,PVALUE=-1e96,HNAME='ADVECTION_METSV::ZRHOY2')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOZ1,PVALUE=-1e97,HNAME='ADVECTION_METSV::ZRHOZ1')
+CALL INIT_ON_HOST_AND_DEVICE(ZRHOZ2,PVALUE=-1e98,HNAME='ADVECTION_METSV::ZRHOZ2')
+!
+CALL MNH_GET_ZT3D(IZ1, IZ2)
+#endif
+!
+IF(LBLOWSNOW) THEN ! Put 2D Canopy blowing snow variables into a 3D array for advection
+#ifdef MNH_OPENACC
+ call Print_msg( NVERB_ERROR, 'GEN', 'ADVECTION_METSV', 'OpenACC: LBLOWSNOW not yet implemented' )
+#endif
+ ZSNWC_INIT(:,:,:,:) = 0.
+ ZRSNWCS(:,:,:,:) = 0.
+
+ DO JSV=1,(NBLOWSNOW_2D)
+ ZSNWC_INIT(:,:,IKB,JSV) = XSNWCANO(:,:,JSV)
+ ZRSNWCS(:,:,IKB,JSV) = XRSNWCANOS(:,:,JSV)
+ END DO
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTES THE CONTRAVARIANT COMPONENTS (FOR PPM ONLY)
+! --------------------------------------
+!
+!* 2.1 computes contravariant components
+!
+!Update on host of ZRUCPPM,ZRVCPPM,ZRWCPPM is done in CONTRAV_DEVICE
+#ifndef MNH_OPENACC
+IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN
+ CALL CONTRAV (HLBCX,HLBCY,PUT,PVT,PWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCPPM,ZRVCPPM,ZRWCPPM,2)
+ELSE
+ CALL CONTRAV (HLBCX,HLBCY,PUT,PVT,PWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCPPM,ZRVCPPM,ZRWCPPM,4)
+END IF
+#else
+IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN
+ CALL CONTRAV_DEVICE (HLBCX,HLBCY,PUT,PVT,PWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCPPM,ZRVCPPM,ZRWCPPM,2, &
+ ZT3D(:,:,:,IZ1),ZT3D(:,:,:,IZ2),ODATA_ON_DEVICE=.TRUE.)
+ELSE
+ CALL CONTRAV_DEVICE (HLBCX,HLBCY,PUT,PVT,PWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCPPM,ZRVCPPM,ZRWCPPM,4, &
+ ZT3D(:,:,:,IZ1),ZT3D(:,:,:,IZ2),ODATA_ON_DEVICE=.TRUE.)
+END IF
+#endif
+!
+!
+!* 2.2 computes CFL numbers
+!
+!PW: not necessary: data already on device due to contrav_device !$acc update device(ZRUCPPM,ZRVCPPM,ZRWCPPM)
+! acc kernels
+IF (.NOT. L1D) THEN
+ !$acc kernels
+ ZCFLU(:,:,:) = 0.0 ; ZCFLV(:,:,:) = 0.0 ; ZCFLW(:,:,:) = 0.0
+ ZCFLU(IIB:IIE,IJB:IJE,:) = ABS(ZRUCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
+ ZCFLV(IIB:IIE,IJB:IJE,:) = ABS(ZRVCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
+ ZCFLW(IIB:IIE,IJB:IJE,:) = ABS(ZRWCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
+ !$acc end kernels
+#ifndef MNH_BITREP
+ IF (.NOT. L2D) THEN
+ !$acc kernels
+ ZCFL(:,:,:) = SQRT(ZCFLU(:,:,:)**2+ZCFLV(:,:,:)**2+ZCFLW(:,:,:)**2)
+ !$acc end kernels
+ ELSE
+ !$acc kernels
+ ZCFL(:,:,:) = SQRT(ZCFLU(:,:,:)**2+ZCFLW(:,:,:)**2)
+ !$acc end kernels
+ END IF
+#else
+ IF (.NOT. L2D) THEN
+ !$acc kernels
+ !$acc loop independent collapse(3)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU )
+ ZCFL(JI,JJ,JK) = SQRT(BR_P2(ZCFLU(JI,JJ,JK))+BR_P2(ZCFLV(JI,JJ,JK))+BR_P2(ZCFLW(JI,JJ,JK)))
+ END DO
+ !$acc end kernels
+ ELSE
+ !$acc kernels
+ !$acc loop independent collapse(3)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU )
+ ZCFL(JI,JJ,JK) = SQRT(BR_P2(ZCFLU(JI,JJ,JK))+BR_P2(ZCFLW(JI,JJ,JK)))
+ END DO
+ !$acc end kernels
+ END IF
+#endif
+ELSE
+ !$acc kernels
+ ZCFLU(:,:,:) = 0.0 ; ZCFLV(:,:,:) = 0.0 ; ZCFLW(:,:,:) = 0.0
+ ZCFLW(IIB:IIE,IJB:IJE,:) = ABS(ZRWCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
+#ifndef MNH_BITREP
+ ZCFL(:,:,:) = SQRT(ZCFLW(:,:,:)**2)
+#else
+ !$acc loop independent collapse(3)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU )
+ ZCFL(JI,JJ,JK) = SQRT(BR_P2(ZCFLW(JI,JJ,JK)))
+ END DO
+#endif
+ !$acc end kernels
+END IF
+! acc end kernels
+!
+!* prints in the file the 3D Courant numbers (one should flag this)
+!
+IF (OCLOSE_OUT .AND. OCFL_WRIT .AND. (.NOT. L1D)) THEN
+!$acc update host(ZCFLU,ZCFLV,ZCFLW,ZCFL)
+ TZFIELD%CMNHNAME = 'CFLU'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'CFLU'
+ TZFIELD%CUNITS = '1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_CFLU'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZCFLU)
+!
+ IF (.NOT. L2D) THEN
+ TZFIELD%CMNHNAME = 'CFLV'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'CFLV'
+ TZFIELD%CUNITS = '1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_CFLV'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZCFLV)
+ END IF
+!
+ TZFIELD%CMNHNAME = 'CFLW'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'CFLW'
+ TZFIELD%CUNITS = '1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_CFLW'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZCFLW)
+!
+ TZFIELD%CMNHNAME = 'CFL'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'CFL'
+ TZFIELD%CUNITS = '1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_CFL'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_Field_write(TPFILE,TZFIELD,ZCFL)
+END IF
+!
+!* prints in the output file the maximum CFL
+!
+#ifndef MNH_OPENACC
+ZCFLU_MAX = MAX_ll(ZCFLU,IINFO_ll)
+ZCFLV_MAX = MAX_ll(ZCFLV,IINFO_ll)
+ZCFLW_MAX = MAX_ll(ZCFLW,IINFO_ll)
+ZCFL_MAX = MAX_ll(ZCFL,IINFO_ll)
+#else
+CALL GET_INDICE_ll( IIB,IJB,IIE,IJE)
+!
+IKB=1+JPVEXT
+IKE=SIZE(ZCFLU,3)-JPVEXT
+!
+!$acc kernels
+ZCFLU_MAX = MAXVAL(ZCFLU(IIB:IIE,IJB:IJE,IKB:IKE))
+ZCFLV_MAX = MAXVAL(ZCFLV(IIB:IIE,IJB:IJE,IKB:IKE))
+ZCFLW_MAX = MAXVAL(ZCFLW(IIB:IIE,IJB:IJE,IKB:IKE))
+ZCFL_MAX = MAXVAL(ZCFL (IIB:IIE,IJB:IJE,IKB:IKE))
+!$acc end kernels
+!
+CALL MPI_ALLREDUCE(MPI_IN_PLACE,ZCFLU_MAX,1,MNHREAL_MPI,MPI_MAX,NMNH_COMM_WORLD,IINFO_ll)
+CALL MPI_ALLREDUCE(MPI_IN_PLACE,ZCFLV_MAX,1,MNHREAL_MPI,MPI_MAX,NMNH_COMM_WORLD,IINFO_ll)
+CALL MPI_ALLREDUCE(MPI_IN_PLACE,ZCFLW_MAX,1,MNHREAL_MPI,MPI_MAX,NMNH_COMM_WORLD,IINFO_ll)
+CALL MPI_ALLREDUCE(MPI_IN_PLACE,ZCFL_MAX,1,MNHREAL_MPI,MPI_MAX,NMNH_COMM_WORLD,IINFO_ll)
+#endif
+!
+WRITE(ILUOUT,FMT='(A24,F10.2,A5,F10.2,A5,F10.2,A9,F10.2)') &
+ 'Max. CFL number for U : ',ZCFLU_MAX, &
+ ' V : ',ZCFLV_MAX,' W : ', ZCFLW_MAX,&
+ 'global : ',ZCFL_MAX
+!
+!
+!* 2.3 updates time step splitting loop
+!
+IF (OSPLIT_CFL .AND. (.NOT.L1D) ) THEN
+!
+ ISPLIT_PPM = INT(ZCFL_MAX/PSPLIT_CFL)+1
+ IF ( KSPLIT /= ISPLIT_PPM ) &
+ WRITE(ILUOUT,FMT='(A37,I2,A4,I2,A11)') &
+ 'PPM time spliting loop changed from ', &
+ KSPLIT,' to ',ISPLIT_PPM, ' iterations'
+!
+ KSPLIT = ISPLIT_PPM
+!
+END IF
+! ---------------------------------------------------------------
+IF (( (ZCFLU_MAX>=3.) .AND. (.NOT.L1D) ) .OR. &
+ ( (ZCFLV_MAX>=3.) .AND. (.NOT.L1D) .AND. (.NOT.L2D) ) .OR. &
+ ( (ZCFLW_MAX>=8.) .AND. (.NOT.L1D) ) ) THEN
+ WRITE(ILUOUT,*) ' '
+ WRITE(ILUOUT,*) ' +---------------------------------------------------+'
+ WRITE(ILUOUT,*) ' | MODEL ERROR |'
+ WRITE(ILUOUT,*) ' +---------------------------------------------------+'
+ WRITE(ILUOUT,*) ' | |'
+ WRITE(ILUOUT,*) ' | The model wind speed becomes too high |'
+ WRITE(ILUOUT,*) ' | |'
+ IF ( ZCFLU_MAX>=3. .OR. ZCFLV_MAX>=3. ) &
+ WRITE(ILUOUT,*) ' | The horizontal CFL value reaches 3. or more |'
+ IF ( ZCFLW_MAX>=8. ) &
+ WRITE(ILUOUT,*) ' | The vertical CFL value reaches 8. or more |'
+ WRITE(ILUOUT,*) ' | |'
+ WRITE(ILUOUT,*) ' | This can be due either to : |'
+ WRITE(ILUOUT,*) ' | - a numerical explosion of the model |'
+ WRITE(ILUOUT,*) ' | - or a too high wind speed for an |'
+ WRITE(ILUOUT,*) ' | acceptable accuracy of the advection |'
+ WRITE(ILUOUT,*) ' | |'
+ WRITE(ILUOUT,*) ' | Please decrease your time-step |'
+ WRITE(ILUOUT,*) ' | |'
+ WRITE(ILUOUT,*) ' +---------------------------------------------------+'
+ WRITE(ILUOUT,*) ' '
+ WRITE(ILUOUT,*) ' +---------------------------------------------------+'
+ WRITE(ILUOUT,*) ' | MODEL STOPS |'
+ WRITE(ILUOUT,*) ' +---------------------------------------------------+'
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','ADVECTION_METSV','')
+END IF
+!
+!
+ZTSTEP_PPM = PTSTEP / REAL(KSPLIT)
+!
+!
+!* 2.4 normalized contravariant components for splitted PPM time-step
+!
+!$acc kernels
+ZRUCPPM(:,:,:) = ZRUCPPM(:,:,:)*ZTSTEP_PPM
+ZRVCPPM(:,:,:) = ZRVCPPM(:,:,:)*ZTSTEP_PPM
+ZRWCPPM(:,:,:) = ZRWCPPM(:,:,:)*ZTSTEP_PPM
+!
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 3. COMPUTES THE TENDENCIES SINCE THE BEGINNING OF THE TIME STEP
+! ------------------------------------------------------------
+!
+!* This represent the effects of all OTHER processes
+! Clouds related processes from previous time-step are taken into account in PRTHS_CLD
+! Advection related processes from previous time-step will be taken into account in ZRTHS_PPM
+!
+ZRTHS_OTHER(:,:,:) = PRTHS(:,:,:) - PTHT(:,:,:) * PRHODJ(:,:,:) / PTSTEP
+IF (GTKE) ZRTKES_OTHER(:,:,:) = PRTKES(:,:,:) - PTKET(:,:,:) * PRHODJ(:,:,:) / PTSTEP
+DO JR = 1, KRR
+ ZRRS_OTHER(:,:,:,JR) = PRRS(:,:,:,JR) - PRT(:,:,:,JR) * PRHODJ(:,:,:) / PTSTEP
+END DO
+DO JSV = 1, KSV
+ ZRSVS_OTHER(:,:,:,JSV) = PRSVS(:,:,:,JSV) - PSVT(:,:,:,JSV) * PRHODJ / PTSTEP
+END DO
+!$acc end kernels
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ ZRSNWCS_OTHER(:,:,:,JSV) = ZRSNWCS(:,:,:,JSV) - ZSNWC_INIT(:,:,:,JSV) * PRHODJ / PTSTEP
+ END DO
+ENDIF
+!
+! Top and bottom Boundaries
+!
+#ifdef MNH_OPENACC
+CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZRTHS_OTHER)
+IF (GTKE) CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZRTKES_OTHER)
+DO JR = 1, KRR
+ CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZRRS_OTHER(:,:,:,JR))
+END DO
+DO JSV = 1, KSV
+ CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZRSVS_OTHER(:,:,:,JSV))
+END DO
+!Already done in ADV_BOUNDARIES_DEVICE !$acc update self(ZRTHS_OTHER,ZRTKES_OTHER,ZRRS_OTHER(:,:,:,1:KRR),ZRSVS_OTHER(:,:,:,1:KSV))
+#else
+CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRTHS_OTHER)
+IF (GTKE) CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRTKES_OTHER)
+DO JR = 1, KRR
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRRS_OTHER(:,:,:,JR))
+END DO
+DO JSV = 1, KSV
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRSVS_OTHER(:,:,:,JSV))
+END DO
+#endif
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRSNWCS_OTHER(:,:,:,JSV))
+ END DO
+END IF
+!
+! Exchanges on processors
+!
+#ifndef MNH_OPENACC
+NULLIFY(TZFIELDS0_ll)
+!!$IF(NHALO == 1) THEN
+ CALL ADD3DFIELD_ll( TZFIELDS0_ll, ZRTHS_OTHER, 'ADVECTION_METSV::ZRTHS_OTHER' )
+ IF (GTKE) CALL ADD3DFIELD_ll( TZFIELDS0_ll, ZRTKES_OTHER, 'ADVECTION_METSV::ZRTKES_OTHER' )
+ IF ( KRR>0 ) CALL ADD4DFIELD_ll( TZFIELDS0_ll, ZRRS_OTHER(:,:,:,1:KRR), 'ADVECTION_METSV::ZRRS_OTHER' )
+ IF ( KSV>0 ) CALL ADD4DFIELD_ll( TZFIELDS0_ll, ZRSVS_OTHER(:,:,:,1:KSV), 'ADVECTION_METSV::ZRSVS_OTHER' )
+ IF(LBLOWSNOW) CALL ADD4DFIELD_ll( TZFIELDS0_ll, ZRSNWCS_OTHER(:,:,:,1:NBLOWSNOW_2D), 'ADVECTION_METSV::ZRSNWCS_OTHER' )
+ CALL UPDATE_HALO_ll(TZFIELDS0_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS0_ll)
+!!$END IF
+#else
+ CALL GET_HALO_D(ZRTHS_OTHER, HNAME='ADVECTION_METSV::ZRTHS_OTHER')
+ IF (GTKE) CALL GET_HALO_D(ZRTKES_OTHER,HNAME='ADVECTION_METSV::ZRTKES_OTHER')
+ DO JR=1,KRR
+ CALL GET_HALO_D(ZRRS_OTHER(:,:,:,JR),HNAME='ADVECTION_METSV::ZRRS_OTHER')
+ END DO
+ DO JSV = 1, KSV
+ CALL GET_HALO_D(ZRSVS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSVS_OTHER')
+ END DO
+ DO JSV = 1,NBLOWSNOW_2D
+ CALL GET_HALO_D(ZRSNWCS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSNWCS_OTHER')
+ END DO
+ !PW: TODO: update only what is needed...
+ ! acc update device(ZRTHS_OTHER,ZRTKES_OTHER,ZRRS_OTHER,ZRSVS_OTHER)
+#endif
+
+
+!
+!
+
+!-------------------------------------------------------------------------------
+!
+!* 4. CALLS THE PPM ADVECTION INSIDE A TIME SPLITTING
+! --------------------------------------
+!
+CALL PPM_RHODJ(HLBCX,HLBCY, ZRUCPPM, ZRVCPPM, ZRWCPPM, &
+ ZTSTEP_PPM, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, &
+ ZRHOZ1, ZRHOZ2 )
+!
+!* values of the fields at the beginning of the time splitting loop
+!$acc kernels
+ZTH(:,:,:) = PTHT(:,:,:)
+IF (KRR /=0 ) ZR(:,:,:,:) = PRT(:,:,:,:)
+IF (KSV /=0 ) ZSV(:,:,:,:) = PSVT(:,:,:,:)
+!
+IF (GTKE) THEN
+ PRTKES_ADV(:,:,:) = 0.
+ ZTKE(:,:,:) = PTKET(:,:,:)
+END IF
+!$acc end kernels
+!
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ ZSNWC(:,:,:,JSV) = ZRSNWCS(:,:,:,JSV)* PTSTEP/ PRHODJ
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSNWC(:,:,:,JSV))
+ END DO
+ ZSNWC_INIT(:,:,:,:)=ZSNWC(:,:,:,:)
+ENDIF
+!
+!* time splitting loop
+DO JSPL=1,KSPLIT
+!
+ !ZRTHS_PPM(:,:,:) = 0.
+ !ZRTKES_PPM(:,:,:) = 0.
+ !IF (KRR /=0) ZRRS_PPM(:,:,:,:) = 0.
+ !IF (KSV /=0) ZRSVS_PPM(:,:,:,:) = 0.
+!
+ IF (LNEUTRAL) THEN
+ !Must be done in a kernels region
+!$acc kernels
+ ZTH(:,:,:)=ZTH(:,:,:)-PTHVREF(:,:,:) !* To be removed with the new PPM scheme ?
+!$acc end kernels
+ END IF
+ CALL PPM_MET (HLBCX,HLBCY, KRR, TPDTCUR, ZRUCPPM, ZRVCPPM, ZRWCPPM, PTSTEP,ZTSTEP_PPM, &
+ PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+ ZTH, ZTKE, ZR, ZRTHS_PPM, ZRTKES_PPM, ZRRS_PPM, HMET_ADV_SCHEME)
+ IF (LNEUTRAL) THEN
+ !Must be done in a kernels region
+!$acc kernels
+ ZTH(:,:,:) = ZTH(:,:,:) + PTHVREF(:,:,:) !* To be removed with the new PPM scheme ?
+!$acc end kernels
+ END IF
+!
+ CALL PPM_SCALAR (HLBCX,HLBCY, KSV, TPDTCUR, ZRUCPPM, ZRVCPPM, ZRWCPPM, PTSTEP, &
+ ZTSTEP_PPM, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+ ZSV, ZRSVS_PPM, HSV_ADV_SCHEME )
+!
+! Tendencies of PPM
+!
+! acc kernels
+ !$acc kernels
+ PRTHS(:,:,:) = PRTHS (:,:,:) + ZRTHS_PPM (:,:,:) / KSPLIT
+ IF (GTKE) PRTKES_ADV(:,:,:) = PRTKES_ADV(:,:,:) + ZRTKES_PPM(:,:,:) / KSPLIT
+ IF (KRR /=0) PRRS (:,:,:,:) = PRRS (:,:,:,:) + ZRRS_PPM (:,:,:,:) / KSPLIT
+ IF (KSV /=0 ) PRSVS (:,:,:,:) = PRSVS (:,:,:,:) + ZRSVS_PPM (:,:,:,:) / KSPLIT
+ !$acc end kernels
+!
+ IF (JSPL<KSPLIT) THEN
+!
+! Guesses of the field inside the time splitting loop
+!
+ !$acc kernels
+ ZTH(:,:,:) = ZTH(:,:,:) + ( ZRTHS_PPM(:,:,:) + ZRTHS_OTHER(:,:,:) + PRTHS_CLD(:,:,:)) * &
+ ZTSTEP_PPM / PRHODJ(:,:,:)
+ !$acc end kernels
+ IF (GTKE) THEN
+ !$acc kernels
+ ZTKE(:,:,:) = ZTKE(:,:,:) + ( ZRTKES_PPM(:,:,:) + ZRTKES_OTHER(:,:,:) ) * ZTSTEP_PPM / PRHODJ(:,:,:)
+ !$acc end kernels
+ END IF
+ !$acc kernels
+ !$acc loop independent collapse(4)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU, JR=1:KRR )
+ ZR(JI,JJ,JK,JR) = ZR(JI,JJ,JK,JR) + ( ZRRS_PPM(JI,JJ,JK,JR) + ZRRS_OTHER(JI,JJ,JK,JR) + PRRS_CLD(JI,JJ,JK,JR) ) &
+ * ZTSTEP_PPM / PRHODJ(JI,JJ,JK)
+ END DO !CONCURRENT
+ !$acc loop seq
+ DO JSV = 1, KSV
+ !$acc loop independent collapse(3)
+ DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
+ ZSV(JI,JJ,JK,JSV) = ZSV(JI,JJ,JK,JSV) + ( ZRSVS_PPM(JI,JJ,JK,JSV) + ZRSVS_OTHER(JI,JJ,JK,JSV) + &
+ PRSVS_CLD(JI,JJ,JK,JSV) ) * ZTSTEP_PPM / PRHODJ(JI,JJ,JK)
+ END DO !CONCURRENT
+ END DO
+ !$acc end kernels
+ END IF
+! acc end kernels
+
+!PW: bug PGI 18.10: not necessary for PRRS,PRSVS but error with decriptor not present
+!$acc update self(PRRS,PRSVS)
+!
+! Top and bottom Boundaries and LBC for the guesses
+!
+ IF (JSPL<KSPLIT) THEN
+#ifndef MNH_OPENACC
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZTH, PTHT )
+ IF (GTKE) CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZTKE, PTKET)
+ DO JR = 1, KRR
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZR(:,:,:,JR), PRT(:,:,:,JR))
+ END DO
+ DO JSV = 1, KSV
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSV(:,:,:,JSV), PSVT(:,:,:,JSV))
+ END DO
+#else
+ CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZTH, PTHT )
+ IF (GTKE) CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZTKE, PTKET)
+ DO JR = 1, KRR
+ CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZR(:,:,:,JR), PRT(:,:,:,JR))
+ END DO
+ DO JSV = 1, KSV
+ CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZSV(:,:,:,JSV), PSVT(:,:,:,JSV))
+ END DO
+!Already done in ADV_BOUNDARIES_DEVICE !$acc update self(ZTH,ZTKE,ZR,ZSV)
+#endif
+
+ IF(LBLOWSNOW) THEN ! Advection of Canopy mass at the 1st atmospheric level
+ ZRSNWCS_PPM(:,:,:,:) = 0.
+ !
+
+ CALL PPM_SCALAR (HLBCX,HLBCY, NBLOWSNOW_2D, TPDTCUR, ZRUCPPM, ZRVCPPM, ZRWCPPM,PTSTEP, &
+ ZTSTEP_PPM, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+ ZSNWC, ZRSNWCS_PPM, HSV_ADV_SCHEME)
+
+
+! Tendencies of PPM
+ ZRSNWCS(:,:,:,:) = ZRSNWCS(:,:,:,:) + ZRSNWCS_PPM (:,:,:,:) / KSPLIT
+! Guesses of the field inside the time splitting loop
+ DO JSV = 1, ( NBLOWSNOW_2D)
+ ZSNWC(:,:,:,JSV) = ZSNWC(:,:,:,JSV) + ZRSNWCS_PPM(:,:,:,JSV)*ZTSTEP_PPM/ PRHODJ(:,:,:)
+ END DO
+
+! Top and bottom Boundaries and LBC for the guesses
+ DO JSV = 1, (NBLOWSNOW_2D)
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSNWC(:,:,:,JSV), ZSNWC_INIT(:,:,:,JSV))
+ END DO
+ END IF
+!
+! Exchanges fields between processors
+!
+#ifndef MNH_OPENACC
+ NULLIFY(TZFIELDS1_ll)
+!!$ IF(NHALO == 1) THEN
+ CALL ADD3DFIELD_ll( TZFIELDS1_ll, ZTH, 'ZTH' )
+ IF (GTKE) CALL ADD3DFIELD_ll( TZFIELDS1_ll, ZTKE, 'ADVECTION_METSV::ZTKE' )
+ IF ( KRR>0 ) CALL ADD4DFIELD_ll( TZFIELDS1_ll, ZR (:,:,:,1:KRR), 'ADVECTION_METSV::ZR' )
+ IF ( KSV>0 ) CALL ADD4DFIELD_ll( TZFIELDS1_ll, ZSV(:,:,:,1:KSV), 'ADVECTION_METSV::ZSV' )
+ IF ( LBLOWSNOW ) CALL ADD4DFIELD_ll( TZFIELDS1_ll, ZSNWC(:,:,:,1:NBLOWSNOW_2D), 'ADVECTION_METSV::ZSNWC' )
+ CALL UPDATE_HALO_ll(TZFIELDS1_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS1_ll)
+!!$ END IF
+#else
+ CALL GET_HALO_D(ZTH,HNAME='ZTH')
+ IF (GTKE) CALL GET_HALO_D(ZTKE,HNAME='ADVECTION_METSV::ZTKE')
+ DO JR=1,KRR
+ CALL GET_HALO_D(ZR(:,:,:,JR),HNAME='ADVECTION_METSV::ZR')
+ END DO
+ DO JSV = 1, KSV
+ CALL GET_HALO_D(ZSV(:,:,:,JSV),HNAME='ADVECTION_METSV::ZSV')
+ END DO
+ DO JSV = 1,NBLOWSNOW_2D
+ CALL GET_HALO_D(ZSNWC(:,:,:,JSV),HNAME='ADVECTION_METSV::ZSNWC')
+ END DO
+#endif
+
+ END IF
+!
+END DO
+!
+!-------------------------------------------------------------------------------
+!
+! TKE special case: advection is the last process for TKE
+!
+! TKE must be greater than its minimum value
+! (previously done in tke_eps_sources)
+!
+IF (GTKE) THEN
+!$acc kernels
+ PRTKES(:,:,:) = PRTKES(:,:,:) + PRTKES_ADV(:,:,:)
+ PRTKES(:,:,:) = MAX (PRTKES(:,:,:) , XTKEMIN * PRHODJ(:,:,:) / PTSTEP )
+!$acc end kernels
+END IF
+!
+!
+!-------------------------------------------------------------------------------
+! Update tendency for cano variables : from 3D to 2D
+!
+IF(LBLOWSNOW) THEN
+
+ DO JSV=1,(NBLOWSNOW_2D)
+ DO JI=1,SIZE(PSVT,1)
+ DO JJ=1,SIZE(PSVT,2)
+ XRSNWCANOS(JI,JJ,JSV) = SUM(ZRSNWCS(JI,JJ,IKB:IKE,JSV))
+ END DO
+ END DO
+ END DO
+IF(LWEST_ll()) XRSNWCANOS(IIB,:,:) = ZRSNWCS(IIB,:,IKB,:)
+IF(LEAST_ll()) XRSNWCANOS(IIE,:,:) = ZRSNWCS(IIE,:,IKB,:)
+IF(LSOUTH_ll()) XRSNWCANOS(:,IJB,:) = ZRSNWCS(:,IJB,IKB,:)
+IF(LNORTH_ll()) XRSNWCANOS(:,IJE,:) = ZRSNWCS(:,IJE,IKB,:)
+
+END IF
+!-------------------------------------------------------------------------------
+!
+!* 5. BUDGETS
+! -------
+!
+IF (LBUDGET_TH) THEN
+!$acc update self(PRTHS)
+ CALL BUDGET (PRTHS,4,'ADV_BU_RTH')
+END IF
+IF (LBUDGET_TKE) THEN
+!$acc update self(PRTKES)
+ CALL BUDGET (PRTKES,5,'ADV_BU_RTKE')
+END IF
+IF (KRR>=1.AND.LBUDGET_RV) THEN
+!$acc update self(PRRS(:,:,:,1))
+ CALL BUDGET (PRRS(:,:,:,1),6,'ADV_BU_RRV')
+END IF
+IF (KRR>=2.AND.LBUDGET_RC) THEN
+!$acc update self(PRRS(:,:,:,2))
+ CALL BUDGET (PRRS(:,:,:,2),7,'ADV_BU_RRC')
+END IF
+IF (KRR>=3.AND.LBUDGET_RR) THEN
+!$acc update self(PRRS(:,:,:,3))
+ CALL BUDGET (PRRS(:,:,:,3),8,'ADV_BU_RRR')
+END IF
+IF (KRR>=4.AND.LBUDGET_RI) THEN
+!$acc update self(PRRS(:,:,:,4))
+ CALL BUDGET (PRRS(:,:,:,4),9,'ADV_BU_RRI')
+END IF
+IF (KRR>=5.AND.LBUDGET_RS) THEN
+!$acc update self(PRRS(:,:,:,5))
+ CALL BUDGET (PRRS(:,:,:,5),10,'ADV_BU_RRS')
+END IF
+IF (KRR>=6.AND.LBUDGET_RG) THEN
+!$acc update self(PRRS(:,:,:,6))
+ CALL BUDGET (PRRS(:,:,:,6),11,'ADV_BU_RRG')
+END IF
+IF (KRR>=7.AND.LBUDGET_RH) THEN
+!$acc update self(PRRS(:,:,:,7))
+ CALL BUDGET (PRRS(:,:,:,7),12,'ADV_BU_RRH')
+END IF
+IF (LBUDGET_SV) THEN
+!$acc update self(PRSVS)
+ DO JSV=1,KSV
+ CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'ADV_BU_RSV')
+ END DO
+END IF
+!
+IF ((HCLOUD == 'KHKO') .OR. (HCLOUD == 'C2R2')) THEN
+#ifdef MNH_OPENACC
+ call Print_msg( NVERB_ERROR, 'GEN', 'ADVECTION_METSV', 'OpenACC: HCLOUD=''KHKO'' OR ''C2R2'' not yet implemented' )
+#endif
+ ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD)
+ ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:)
+ ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT)
+ ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT)
+ ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1)
+! CALL GET_HALO(PRRS(:,:,:,2))
+! CALL GET_HALO(PRSVS(:,:,:,2))
+! CALL GET_HALO(PRSVS(:,:,:,3))
+ WHERE (PRRS(:,:,:,2) < 0. .OR. PRSVS(:,:,:,2) < 0.)
+ PRSVS(:,:,:,1) = 0.0
+ END WHERE
+ DO JSV = 2, 3
+ WHERE (PRRS(:,:,:,JSV) < 0. .OR. PRSVS(:,:,:,JSV) < 0.)
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,JSV)
+ PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,JSV) * ZLV(:,:,:) / &
+ ZCPH(:,:,:) / ZEXN(:,:,:)
+ PRRS(:,:,:,JSV) = 0.0
+ PRSVS(:,:,:,JSV) = 0.0
+ END WHERE
+ END DO
+!
+ IF (LBUDGET_TH) CALL BUDGET (PRTHS(:,:,:) , 4,'NEADV_BU_RTH')
+ IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NEADV_BU_RRV')
+ IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NEADV_BU_RRC')
+
+END IF
+!
+#ifdef MNH_OPENACC
+CALL MNH_REL_ZT3D(IZ1, IZ2)
+#endif
+!-------------------------------------------------------------------------------
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRTHS,"ADVECTION_METSV end:PRTHS")
+ CALL MPPDB_CHECK(PRTKES,"ADVECTION_METSV end:PRTKES")
+ CALL MPPDB_CHECK(PRRS,"ADVECTION_METSV end:PRRS")
+ CALL MPPDB_CHECK(PRSVS,"ADVECTION_METSV end:PRSVS")
+ !Check all OUT arrays
+ CALL MPPDB_CHECK(PRTKES_ADV,"ADVECTION_METSV end:PRTKES_ADV")
+END IF
+
+!$acc end data
+
+#ifndef MNH_OPENACC
+deallocate ( ZRUCPPM, ZRVCPPM, ZRWCPPM, ZCFLU, ZCFLV, ZCFLW, ZCFL, ZTH, &
+ ZTKE, ZRTHS_OTHER, ZRTKES_OTHER, ZRTHS_PPM, ZRTKES_PPM, &
+ ZR, ZSV, ZSNWC, ZSNWC_INIT, ZRSNWCS, ZRRS_OTHER, ZRSVS_OTHER, ZRSNWCS_OTHER, &
+ ZRRS_PPM, ZRSVS_PPM, ZRSNWCS_PPM, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+ ZT, ZEXN, ZLV, ZLS, ZCPH )
+#else
+CALL MNH_REL_ZT3D ( IZRHOX1, IZRHOX2, IZRHOY1, IZRHOY2, IZRHOZ1, IZRHOZ2, &
+ IZT, IZEXN, IZLV, IZLS, IZCPH )
+
+CALL MNH_REL_ZT4D ( NBLOWSNOW_2D , IZRSNWCS_PPM )
+CALL MNH_REL_ZT4D ( SIZE(PSVT,4) , IZRSVS_PPM )
+CALL MNH_REL_ZT4D ( SIZE(PRT, 4) , IZRRS_PPM )
+CALL MNH_REL_ZT4D ( NBLOWSNOW_2D , IZRSNWCS_OTHER )
+CALL MNH_REL_ZT4D ( SIZE(PSVT,4) , IZRSVS_OTHER )
+CALL MNH_REL_ZT4D ( SIZE(PRT, 4) , IZRRS_OTHER )
+CALL MNH_REL_ZT4D ( NBLOWSNOW_2D , IZRSNWCS )
+CALL MNH_REL_ZT4D ( NBLOWSNOW_2D , IZSNWC_INIT )
+CALL MNH_REL_ZT4D ( NBLOWSNOW_2D , IZSNWC )
+CALL MNH_REL_ZT4D ( SIZE(PSVT,4) , IZSV )
+CALL MNH_REL_ZT4D ( SIZE(PRT, 4) , IZR )
+
+CALL MNH_REL_ZT3D ( IZRUCPPM, IZRVCPPM, IZRWCPPM, IZCFLU, IZCFLV, IZCFLW, IZCFL, IZTH, &
+ IZTKE, IZRTHS_OTHER, IZRTKES_OTHER, IZRTHS_PPM, IZRTKES_PPM )
+CALL MNH_CHECK_OUT_ZT3D("ADVECTION_METSV")
+#endif
+
+!$acc end data
+
+END SUBROUTINE ADVECTION_METSV
diff --git a/src/ZSOLVER/advection_uvw_cen.f90 b/src/ZSOLVER/advection_uvw_cen.f90
new file mode 100644
index 0000000000000000000000000000000000000000..31cdeeae716c3ab246577ab203467280886f714e
--- /dev/null
+++ b/src/ZSOLVER/advection_uvw_cen.f90
@@ -0,0 +1,368 @@
+!MNH_LIC Copyright 2013-2020 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_ADVECTION_UVW_CEN
+! #####################
+!
+INTERFACE
+ SUBROUTINE ADVECTION_UVW_CEN(HUVW_ADV_SCHEME, &
+ HLBCX, HLBCY, &
+ PTSTEP, KTCOUNT, &
+ PUM, PVM, PWM, &
+ PDUM, PDVM, PDWM, &
+ PUT, PVT, PWT, &
+ PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRUS,PRVS, PRWS, &
+ TPHALO2MLIST )
+!
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+!
+CHARACTER(LEN=6), INTENT(IN) :: HUVW_ADV_SCHEME
+!
+CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
+!
+REAL, INTENT(IN) :: PTSTEP! time step
+INTEGER, INTENT(IN) :: KTCOUNT
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM, PVM, PWM
+ ! Variables at t-dt
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDUM, PDVM, PDWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT , PVT , PWT, PRHODJ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
+ ! metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS , PRVS , PRWS
+ ! Sources terms
+!
+! halo lists for 4th order advection
+TYPE(HALO2LIST_ll), POINTER :: TPHALO2MLIST ! momentum variables
+!
+END SUBROUTINE ADVECTION_UVW_CEN
+!
+END INTERFACE
+!
+END MODULE MODI_ADVECTION_UVW_CEN
+! ##########################################################################
+ SUBROUTINE ADVECTION_UVW_CEN(HUVW_ADV_SCHEME, &
+ HLBCX, HLBCY, &
+ PTSTEP, KTCOUNT, &
+ PUM, PVM, PWM, &
+ PDUM, PDVM, PDWM, &
+ PUT, PVT, PWT, &
+ PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRUS,PRVS, PRWS, &
+ TPHALO2MLIST )
+! ##########################################################################
+!
+!!**** *ADVECTION * - routine to call the specialized advection routines
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to control the advection routines.
+!! For that, it is first necessary to compute the metric coefficients
+!! and the contravariant components of the momentum.
+!!
+!!** METHOD
+!! ------
+!! The advection of momenta is calculated using a centred (second order)
+!! scheme.
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! NONE
+!!
+!! REFERENCE
+!! ---------
+!! Book1 and book2 ( routine ADVECTION )
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 01/2013 (from ADVECTION routine)
+!! Modif
+!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODE_ll
+use mode_mppdb
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
+USE MODD_CONF
+USE MODD_PARAMETERS
+USE MODD_GRID_n
+!
+#ifndef MNH_OPENACC
+USE MODI_SHUMAN
+#else
+USE MODI_SHUMAN_DEVICE
+#endif
+USE MODI_CONTRAV
+USE MODI_ADVECUVW_2ND
+USE MODI_ADVECUVW_4TH
+!
+USE MODD_BUDGET
+USE MODI_BUDGET
+!
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
+use mode_msg
+#endif
+!
+!-------------------------------------------------------------------------------
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+CHARACTER(LEN=6), INTENT(IN) :: HUVW_ADV_SCHEME
+!
+CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
+!
+REAL, INTENT(IN) :: PTSTEP! time step
+INTEGER, INTENT(IN) :: KTCOUNT
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUM, PVM, PWM
+ ! Variables at t-dt
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDUM, PDVM, PDWM
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT , PVT , PWT, PRHODJ
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
+ ! metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS , PRVS , PRWS
+ ! Sources terms
+!
+! halo lists for 4th order advection
+TYPE(HALO2LIST_ll), POINTER :: TPHALO2MLIST ! momentum variables
+!
+!
+!* 0.2 declarations of local variables
+!
+!
+REAL, DIMENSION(:,:,:), allocatable :: ZUS
+REAL, DIMENSION(:,:,:), allocatable :: ZVS
+REAL, DIMENSION(:,:,:), allocatable :: ZWS
+ ! guess of cartesian components of
+ ! momentum at future (+PTSTEP) timestep
+REAL, DIMENSION(:,:,:), allocatable :: ZRUS
+REAL, DIMENSION(:,:,:), allocatable :: ZRVS
+REAL, DIMENSION(:,:,:), allocatable :: ZRWS
+ ! cartesian components of
+ ! rhodJ times the tendency of
+ ! momentum from previous (-PTSTEP)
+ ! to future (+PTSTEP) timestep
+!
+REAL, DIMENSION(:,:,:), allocatable :: ZRUT
+REAL, DIMENSION(:,:,:), allocatable :: ZRVT
+REAL, DIMENSION(:,:,:), allocatable :: ZRWT
+ ! cartesian
+ ! components of
+ ! momentum
+!
+REAL, DIMENSION(:,:,:), allocatable :: ZRUCT
+REAL, DIMENSION(:,:,:), allocatable :: ZRVCT
+REAL, DIMENSION(:,:,:), allocatable :: ZRWCT
+ ! contravariant
+ ! components
+ ! of momentum
+REAL, DIMENSION(:,:,:), allocatable :: ZMXM_RHODJ
+REAL, DIMENSION(:,:,:), allocatable :: ZMYM_RHODJ
+REAL, DIMENSION(:,:,:), allocatable :: ZMZM_RHODJ
+!
+INTEGER :: IINFO_ll ! return code of parallel routine
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+#ifdef MNH_OPENACC
+INTEGER :: IZ1, IZ2
+#endif
+!
+!-------------------------------------------------------------------------------
+!$acc data present( PUM, PVM, PWM, PDUM, PDVM, PDWM, PUT, PVT, PWT, PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, PRUS, PRVS, PRWS )
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all IN arrays
+ CALL MPPDB_CHECK(PUM,"ADVECTION_UVW_CEN beg:PUM")
+ CALL MPPDB_CHECK(PVM,"ADVECTION_UVW_CEN beg:PVM")
+ CALL MPPDB_CHECK(PWM,"ADVECTION_UVW_CEN beg:PWM")
+ CALL MPPDB_CHECK(PDUM,"ADVECTION_UVW_CEN beg:PDUM")
+ CALL MPPDB_CHECK(PDVM,"ADVECTION_UVW_CEN beg:PDVM")
+ CALL MPPDB_CHECK(PDWM,"ADVECTION_UVW_CEN beg:PDWM")
+ CALL MPPDB_CHECK(PUT,"ADVECTION_UVW_CEN beg:PUT")
+ CALL MPPDB_CHECK(PVT,"ADVECTION_UVW_CEN beg:PVT")
+ CALL MPPDB_CHECK(PWT,"ADVECTION_UVW_CEN beg:PWT")
+ CALL MPPDB_CHECK(PRHODJ,"ADVECTION_UVW_CEN beg:PRHODJ")
+ CALL MPPDB_CHECK(PDXX,"ADVECTION_UVW_CEN beg:PDXX")
+ CALL MPPDB_CHECK(PDYY,"ADVECTION_UVW_CEN beg:PDYY")
+ CALL MPPDB_CHECK(PDZZ,"ADVECTION_UVW_CEN beg:PDZZ")
+ CALL MPPDB_CHECK(PDZX,"ADVECTION_UVW_CEN beg:PDZX")
+ CALL MPPDB_CHECK(PDZY,"ADVECTION_UVW_CEN beg:PDZY")
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRUS,"ADVECTION_UVW_CEN beg:PRUS")
+ CALL MPPDB_CHECK(PRVS,"ADVECTION_UVW_CEN beg:PRVS")
+ CALL MPPDB_CHECK(PRWS,"ADVECTION_UVW_CEN beg:PRWS")
+END IF
+
+allocate( zus ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zvs ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zws ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrus ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrvs ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrws ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrut ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrvt ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrwt ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zruct ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrvct ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zrwct ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zmxm_rhodj ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zmym_rhodj ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+allocate( zmzm_rhodj ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
+
+!$acc data create( zus, zvs, zws, zrus, zrvs, zrws, zrut, zrvt, zrwt, &
+!$acc & zruct, zrvct, zrwct, zmxm_rhodj, zmym_rhodj, zmzm_rhodj )
+
+#ifdef MNH_OPENACC
+CALL INIT_ON_HOST_AND_DEVICE(ZUS,-1e99,'ADVECTION_UVW_CEN::ZUS')
+CALL INIT_ON_HOST_AND_DEVICE(ZVS,-2e99,'ADVECTION_UVW_CEN::ZVS')
+CALL INIT_ON_HOST_AND_DEVICE(ZWS,-3e99,'ADVECTION_UVW_CEN::ZWS')
+CALL INIT_ON_HOST_AND_DEVICE(ZRUS,-1e99,'ADVECTION_UVW_CEN::ZRUS')
+CALL INIT_ON_HOST_AND_DEVICE(ZRVS,-2e99,'ADVECTION_UVW_CEN::ZRVS')
+CALL INIT_ON_HOST_AND_DEVICE(ZRWS,-3e99,'ADVECTION_UVW_CEN::ZRWS')
+CALL INIT_ON_HOST_AND_DEVICE(ZRUT,-1e99,'ADVECTION_UVW_CEN::ZRUT')
+CALL INIT_ON_HOST_AND_DEVICE(ZRVT,-2e99,'ADVECTION_UVW_CEN::ZRVT')
+CALL INIT_ON_HOST_AND_DEVICE(ZRWT,-3e99,'ADVECTION_UVW_CEN::ZRWT')
+CALL INIT_ON_HOST_AND_DEVICE(ZRUCT,-1e98,'ADVECTION_UVW_CEN::ZRUCT')
+CALL INIT_ON_HOST_AND_DEVICE(ZRVCT,-2e98,'ADVECTION_UVW_CEN::ZRVCT')
+CALL INIT_ON_HOST_AND_DEVICE(ZRWCT,-3e98,'ADVECTION_UVW_CEN::ZRWCT')
+CALL INIT_ON_HOST_AND_DEVICE(ZMXM_RHODJ,-1e97,'ADVECTION_UVW_CEN::ZMXM_RHODJ')
+CALL INIT_ON_HOST_AND_DEVICE(ZMYM_RHODJ,-2e97,'ADVECTION_UVW_CEN::ZMYM_RHODJ')
+CALL INIT_ON_HOST_AND_DEVICE(ZMZM_RHODJ,-3e97,'ADVECTION_UVW_CEN::ZMZM_RHODJ')
+!
+CALL MNH_GET_ZT3D(IZ1, IZ2)
+#endif
+!
+#ifndef MNH_OPENACC
+ZMXM_RHODJ = MXM(PRHODJ)
+ZMYM_RHODJ = MYM(PRHODJ)
+ZMZM_RHODJ = MZM(PRHODJ)
+#else
+CALL MXM_DEVICE(PRHODJ,ZMXM_RHODJ)
+CALL MYM_DEVICE(PRHODJ,ZMYM_RHODJ)
+CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ)
+#endif
+!
+!* 1. COMPUTES THE CONTRAVARIANT COMPONENTS
+! -------------------------------------
+!
+!$acc kernels present(ZRUT,ZRVT,ZRWT,PUT,PVT,PWT,ZMXM_RHODJ,ZMYM_RHODJ,ZMZM_RHODJ)
+ZRUT(:,:,:) = PUT(:,:,:) * ZMXM_RHODJ(:,:,:)
+ZRVT(:,:,:) = PVT(:,:,:) * ZMYM_RHODJ(:,:,:)
+ZRWT(:,:,:) = PWT(:,:,:) * ZMZM_RHODJ(:,:,:)
+!$acc end kernels
+!
+#ifndef MNH_OPENACC
+IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN
+ CALL CONTRAV (HLBCX,HLBCY,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT,2)
+ELSEIF (HUVW_ADV_SCHEME=='CEN4TH') THEN
+ CALL CONTRAV (HLBCX,HLBCY,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT,4)
+END IF
+#else
+IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN
+ CALL CONTRAV_DEVICE (HLBCX,HLBCY,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT,2, &
+ ZT3D(:,:,:,IZ1),ZT3D(:,:,:,IZ2),ODATA_ON_DEVICE=.TRUE.)
+ELSEIF (HUVW_ADV_SCHEME=='CEN4TH') THEN
+ CALL CONTRAV_DEVICE (HLBCX,HLBCY,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT,4, &
+ ZT3D(:,:,:,IZ1),ZT3D(:,:,:,IZ2),ODATA_ON_DEVICE=.TRUE.)
+END IF
+!Not necessary: already done in contrav_device !$acc update self(ZRUCT,ZRVCT,ZRWCT)
+#endif
+!
+NULLIFY(TZFIELDS_ll)
+!!$IF(NHALO == 1) THEN
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRWCT, 'ADVECTION_UVW_CEN::ZRWCT' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRUCT, 'ADVECTION_UVW_CEN::ZRUCT' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRVCT, 'ADVECTION_UVW_CEN::ZRVCT' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ !$acc update device(ZRUCT, ZRVCT, ZRWCT)
+!!$END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. TERM FROM PREVIOUS TIME-STEP (from initial_guess)
+! ----------------------------
+!
+!$acc kernels present(ZRUS,ZRVS,ZRWS,PUM,PVM,PWM,ZMXM_RHODJ,ZMYM_RHODJ,ZMZM_RHODJ)
+ZRUS(:,:,:) = PUM(:,:,:) * ZMXM_RHODJ(:,:,:)/(2.*PTSTEP)
+ZRVS(:,:,:) = PVM(:,:,:) * ZMYM_RHODJ(:,:,:)/(2.*PTSTEP)
+ZRWS(:,:,:) = PWM(:,:,:) * ZMZM_RHODJ(:,:,:)/(2.*PTSTEP)
+!$acc end kernels
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. CALLS THE ADVECTION ROUTINES FOR THE MOMENTUM
+! ---------------------------------------------
+!
+! choose between 2nd and 4th order momentum advection.
+IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN
+!
+#ifdef MNH_OPENACC
+ call Print_msg( NVERB_ERROR, 'GEN', 'ADVECTION_UVW_CEN', 'OpenACC: CEN2ND not yet implemented' )
+#endif
+ CALL ADVECUVW_2ND (PUT,PVT,PWT,ZRUCT,ZRVCT,ZRWCT,ZRUS,ZRVS,ZRWS)
+!
+ELSEIF (HUVW_ADV_SCHEME=='CEN4TH') THEN
+!
+ CALL ADVECUVW_4TH ( HLBCX, HLBCY, ZRUCT, ZRVCT, ZRWCT, &
+ PUT, PVT, PWT, ZRUS, ZRVS, ZRWS, TPHALO2MLIST )
+!
+END IF
+!
+!$acc kernels
+ZUS(:,:,:) = ZRUS(:,:,:)/ZMXM_RHODJ(:,:,:)*2.*PTSTEP
+ZVS(:,:,:) = ZRVS(:,:,:)/ZMYM_RHODJ(:,:,:)*2.*PTSTEP
+ZWS(:,:,:) = ZRWS(:,:,:)/ZMZM_RHODJ(:,:,:)*2.*PTSTEP
+!-------------------------------------------------------------------------------
+!
+!* 5. Extracts the variation between current and future time step
+! -----------------------------------------------------------
+!
+PRUS(:,:,:) = PRUS(:,:,:) + ( ZUS(:,:,:) - PUM(:,:,:) - 0.5* PDUM(:,:,:)) * ZMXM_RHODJ(:,:,:)/(PTSTEP)
+PRVS(:,:,:) = PRVS(:,:,:) + ( ZVS(:,:,:) - PVM(:,:,:) - 0.5* PDVM(:,:,:)) * ZMYM_RHODJ(:,:,:)/(PTSTEP)
+PRWS(:,:,:) = PRWS(:,:,:) + ( ZWS(:,:,:) - PWM(:,:,:) - 0.5* PDWM(:,:,:)) * ZMZM_RHODJ(:,:,:)/(PTSTEP)
+!
+PDUM(:,:,:) = ZUS(:,:,:) - PUM(:,:,:)
+PDVM(:,:,:) = ZVS(:,:,:) - PVM(:,:,:)
+PDWM(:,:,:) = ZWS(:,:,:) - PWM(:,:,:)
+!$acc end kernels
+!
+IF (LBUDGET_U) CALL BUDGET (PRUS,1,'ADV_BU_RU')
+IF (LBUDGET_V) CALL BUDGET (PRVS,2,'ADV_BU_RV')
+IF (LBUDGET_W) CALL BUDGET (PRWS,3,'ADV_BU_RW')
+!
+#ifdef MNH_OPENACC
+CALL MNH_REL_ZT3D(IZ1, IZ2)
+#endif
+
+IF (MPPDB_INITIALIZED) THEN
+ !Check all INOUT arrays
+ CALL MPPDB_CHECK(PRUS,"ADVECTION_UVW_CEN end:PRUS")
+ CALL MPPDB_CHECK(PRVS,"ADVECTION_UVW_CEN end:PRVS")
+ CALL MPPDB_CHECK(PRWS,"ADVECTION_UVW_CEN end:PRWS")
+END IF
+
+!$acc end data
+
+!$acc end data
+
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE ADVECTION_UVW_CEN
diff --git a/src/ZSOLVER/conjgrad.f90 b/src/ZSOLVER/conjgrad.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9adf9ed36ed2cb0a64ca545ba5c9c378d35b5224
--- /dev/null
+++ b/src/ZSOLVER/conjgrad.f90
@@ -0,0 +1,291 @@
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 solver 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ####################
+ MODULE MODI_CONJGRAD
+! ####################
+!
+INTERFACE
+!
+ SUBROUTINE CONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI)
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t
+!
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+END SUBROUTINE CONJGRAD
+!
+END INTERFACE
+!
+END MODULE MODI_CONJGRAD
+!
+!
+!
+! #########################################################################
+ SUBROUTINE CONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI)
+! #########################################################################
+!
+!!**** *CONJGRAD * - solve an elliptic equation by the conjugate gradient
+!! method
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to solve an elliptic equation using
+! the preditioned conjugate gradient (CG) method. This is a version of the
+! CG called ORTHOMIN (Young and Jea 1980).
+!
+!!** METHOD
+!! ------
+!! The equation to be solved reads:
+!!
+!! Q (PHI) = Y
+!!
+!! where Q is the quasi-Laplacian ( subroutine QLAP ) and PHI the pressure
+!! function.
+!! We precondition the problem by the operator F :
+!! -1 -1
+!! F * Q (PHI) = F (Y)
+!! F represents the flat Laplacian ie. without orography. Its inversion is
+!! realized in the routine FLAT_INV. This equation is solved with a Conjugate
+!! Gradient method.
+!! The initial guess is given by the pressure at the previous time step.
+!! The resolution stops after ITR iterations of the solver.
+!!
+!! EXTERNAL
+!! --------
+!! Subroutine GDIV: compute J times the divergence of 1/J times a vector
+!! Function QLAP: compute the complete quasi-Laplacian Q
+!! Subroutine FLAT_INV : invert the flat quasi-laplacien F
+!! Function DOTPROD: compute the dot product of 2 vectors
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODI_GDIV: interface for the subroutine GDIV
+!! Module MODI_QLAP: interface for the function QLAP
+!! Module MODI_FLAT_INV: interface for the subroutine FLAT_INV
+!! Module MODI_DOTPROD: interface for the function DOTPROD
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine CONJGRAD)
+!! Kapitza and Eppel (1992) Beit. Physik ...
+!! Young and Jea (1980) ....
+!!
+!! AUTHOR
+!! ------
+!! P. HÅreil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/07/94
+!!
+!! 14/01/97 Durran anelastic equation (Stein,Lafore)
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_GDIV
+USE MODI_QLAP
+USE MODI_FLAT_INV
+USE MODI_DOTPROD
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of 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
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t
+!
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: JM ! loop index
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZDELTA
+ ! array containing the auxilary field DELTA of the CG method
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZP
+ ! array containing the auxilary field P of the CG method
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZWORK ! work
+ ! array containing the source term to be multiplied by the F inverse
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZWORKD ! work
+ ! array containing the result of the F inversion * Q (DELTA)
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZWORKP ! work
+ ! array containing the result of the F inversion * Q (P)
+!
+REAL :: ZALPHA, ZLAMBDA ! amplitude of the descent in the Conjugate
+ ! directions
+REAL :: ZDOTPP ! dot product of ZWORKP by itself
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+ZLAMBDA = 0.
+ZP = 0.
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. ITERATIVE LOOP
+! --------------
+!
+DO JM = 1,KITR
+!
+!* 2.1 compute the new pressure function
+!
+ PPHI = PPHI + ZLAMBDA * ZP ! the case JM =0 is special because
+ ! PPHI is not changed
+!
+!* 2.2 compute the auxiliary field DELTA
+!
+! -1
+! compute the vector DELTA = F * ( Y - Q ( PHI ) )
+!
+ ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI)
+ ! Q (PHI)
+!
+ ZWORK = PY - ZWORK ! Y - Q (PHI)
+!
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZDELTA) ! F (Y - Q (PHI)))
+!
+!* 2.3 compute the auxiliary field P
+!
+! -1
+! compute the vector P = DELTA + alpha F * Q ( DELTA )
+!
+ IF (JM == 1) THEN
+ ZP = ZDELTA ! P = DELTA at the first solver iteration
+ ELSE
+ ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY, &
+ PDZZ,PRHODJ,PTHETAV,ZDELTA) ! Q ( DELTA )
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKD) ! F * Q ( DELTA )
+!
+ ZALPHA = - DOTPROD(ZWORKD,ZWORKP,HLBCX,HLBCY)/ZDOTPP ! ZWORKP,ZDOTPP come
+ ! from the previous solver iteration (section 2.4)
+ ZP = ZDELTA + ZALPHA * ZP ! new vector P
+!
+ END IF
+!
+!* 2.4 compute LAMBDA
+!
+!
+ ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,&
+ PDZZ,PRHODJ,PTHETAV,ZP) ! Q ( P )
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,& ! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKP) ! F * Q ( P )
+!
+! store the scalar product to compute lambda and next P
+ ZDOTPP = DOTPROD(ZWORKP,ZWORKP,HLBCX,HLBCY)
+!
+ ZLAMBDA = DOTPROD(ZDELTA,ZWORKP,HLBCX,HLBCY) / ZDOTPP ! lambda
+!
+!
+END DO ! end of the loop for the iterative solver
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. COMPUTE THE FINAL PRESSURE FUNCTION
+! -----------------------------------
+!
+PPHI = PPHI + ZLAMBDA * ZP
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE CONJGRAD
diff --git a/src/ZSOLVER/conresol.f90 b/src/ZSOLVER/conresol.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c0e103e9a6b06282ccc14e37f90535404b147340
--- /dev/null
+++ b/src/ZSOLVER/conresol.f90
@@ -0,0 +1,272 @@
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 solver 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ####################
+ MODULE MODI_CONRESOL
+! ####################
+!
+INTERFACE
+!
+ SUBROUTINE CONRESOL(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI)
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual pot. temp. at time t
+!
+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) :: PRHOM ! XRHODJ mean on the X Y plane
+ ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+END SUBROUTINE CONRESOL
+!
+END INTERFACE
+!
+END MODULE MODI_CONRESOL
+!
+!
+!
+! #########################################################################
+ SUBROUTINE CONRESOL(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI)
+! #########################################################################
+!
+!!**** *CONRESOL * - solve an elliptic equation by the conjugate residual
+!! method
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to solve an elliptic equation using
+! the preconditioned conjugate residual (CR) method. This is a version
+! of the scheme proposed by Skamarock, Smolarkiewicz and Klemp (MWR, 1997).
+!
+!!** METHOD
+!! ------
+!! The equation to be solved reads:
+!!
+!! Q (PHI) = Y
+!!
+!! where Q is the quasi-Laplacian ( subroutine QLAP ) and PHI the pressure
+!! function.
+!! We precondition the problem by the operator F :
+!! -1 -1
+!! F * Q (PHI) = F (Y)
+!! F represents the flat Laplacian ie. without orography. Its inversion is
+!! realized in the routine FLAT_INV. This equation is solved with a Conjugate
+!! Residual method.
+!! The initial guess is given by the pressure at the previous time step.
+!! The resolution stops after ITR iterations of the solver.
+!!
+!! EXTERNAL
+!! --------
+!! Subroutine GDIV: compute J times the divergence of 1/J times a vector
+!! Function QLAP: compute the complete quasi-Laplacian Q
+!! Subroutine FLAT_INV : invert the flat quasi-laplacien F
+!! Function DOTPROD: compute the dot product of 2 vectors
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODI_GDIV: interface for the subroutine GDIV
+!! Module MODI_QLAP: interface for the function QLAP
+!! Module MODI_FLAT_INV: interface for the subroutine FLAT_INV
+!! Module MODI_DOTPROD: interface for the function DOTPROD
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine CONRESOL)
+!! Skamarock, Smolarkiewicz and Klemp (1997) MWR
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty *Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/08/99
+!! J.-P. Pinty & P. Jabouille
+!! 11/07/00 bug in ZALPHA
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_GDIV
+USE MODI_QLAP
+USE MODI_FLAT_INV
+USE MODI_DOTPROD
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of 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
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual pot. temp. at time t
+!
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: JM ! loop index
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZDELTA, ZKSI
+ ! array containing the auxilary fields DELTA and KSI of the CR method
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZP, ZQ
+ ! array containing the auxilary fields P and Q of the CR method
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZRESIDUE
+ ! array containing the error field at each iteration Q(PHI) - Y
+!
+REAL :: ZALPHA, ZLAMBDA ! amplitude of the descent in the Conjugate
+ ! directions
+REAL :: ZDOT_DELTA ! dot product of ZDELTA by itself
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!
+!* 1.1 compute the vector: r^(0) = Q(PHI) - Y
+!
+ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY
+!
+!* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y )
+!
+CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZP)
+!
+!* 1.3 compute the vector: delta^(0) = Q ( p^(0) )
+!
+ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. ITERATIVE LOOP
+! --------------
+!
+DO JM = 1,KITR
+!
+!* 2.1 compute the step LAMBDA
+!
+ ZDOT_DELTA = DOTPROD(ZDELTA, ZDELTA,HLBCX,HLBCY) ! norm of DELTA
+ ZLAMBDA = - DOTPROD(ZRESIDUE,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA
+!
+!* 2.2 update the pressure function PHI
+!
+ PPHI = PPHI + ZLAMBDA * ZP
+!
+!
+ IF( JM == KITR ) EXIT
+!
+!
+!* 2.3 update the residual error: r
+!
+ ZRESIDUE = ZRESIDUE + ZLAMBDA * ZDELTA
+!
+!* 2.4 compute the vector: q = F^(-1)*( Q(PHI) - Y )
+!
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZQ)
+!
+!* 2.5 compute the auxiliary field: ksi = Q ( q )
+!
+ ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ)
+! -1
+!* 2.6 compute the step ALPHA
+!
+ ZALPHA = - DOTPROD(ZKSI,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA ! lambda
+!
+!* 2.7 update p and DELTA
+!
+ ZP = ZQ + ZALPHA * ZP
+ ZDELTA = ZKSI + ZALPHA * ZDELTA
+!
+END DO ! end of the loop for the iterative solver
+!
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE CONRESOL
diff --git a/src/ZSOLVER/conresolz.f90 b/src/ZSOLVER/conresolz.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3624bb2eb58b095c96245e1894b2ad8f4f9c71e8
--- /dev/null
+++ b/src/ZSOLVER/conresolz.f90
@@ -0,0 +1,295 @@
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$ $Date$
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
+! ####################
+ MODULE MODI_CONRESOLZ
+! ####################
+!
+INTERFACE
+!
+ SUBROUTINE CONRESOLZ(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI, &
+ PBFB,&
+ PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual pot. temp. at time t
+!
+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) :: PRHOM ! XRHODJ mean on the X Y plane
+ ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+!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.
+!JUAN Z_SPLITING
+END SUBROUTINE CONRESOLZ
+!
+END INTERFACE
+!
+END MODULE MODI_CONRESOLZ
+!
+!
+!
+! #########################################################################
+ SUBROUTINE CONRESOLZ(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI, &
+ PBFB,&
+ PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+! #########################################################################
+!
+!!**** *CONRESOLZ * - solve an elliptic equation by the conjugate residual
+!! method
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to solve an elliptic equation using
+! the preconditioned conjugate residual (CR) method. This is a version
+! of the scheme proposed by Skamarock, Smolarkiewicz and Klemp (MWR, 1997).
+!
+!!** METHOD
+!! ------
+!! The equation to be solved reads:
+!!
+!! Q (PHI) = Y
+!!
+!! where Q is the quasi-Laplacian ( subroutine QLAP ) and PHI the pressure
+!! function.
+!! We precondition the problem by the operator F :
+!! -1 -1
+!! F * Q (PHI) = F (Y)
+!! F represents the flat Laplacian ie. without orography. Its inversion is
+!! realized in the routine FLAT_INVZ. This equation is solved with a Conjugate
+!! Residual method.
+!! The initial guess is given by the pressure at the previous time step.
+!! The resolution stops after ITR iterations of the solver.
+!!
+!! EXTERNAL
+!! --------
+!! Subroutine GDIV: compute J times the divergence of 1/J times a vector
+!! Function QLAP: compute the complete quasi-Laplacian Q
+!! Subroutine FLAT_INVZ : invert the flat quasi-laplacien F
+!! Function DOTPROD: compute the dot product of 2 vectors
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODI_GDIV: interface for the subroutine GDIV
+!! Module MODI_QLAP: interface for the function QLAP
+!! Module MODI_FLAT_INVZ: interface for the subroutine FLAT_INVZ
+!! Module MODI_DOTPROD: interface for the function DOTPROD
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine CONRESOL)
+!! Skamarock, Smolarkiewicz and Klemp (1997) MWR
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty *Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/08/99
+!! J.-P. Pinty & P. Jabouille
+!! 11/07/00 bug in ZALPHA
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODI_GDIV
+USE MODI_QLAP
+USE MODI_FLAT_INVZ
+USE MODI_DOTPROD
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of 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
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual pot. temp. at time t
+!
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!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.
+!JUAN Z_SPLITING
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: JM ! loop index
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZDELTA, ZKSI
+ ! array containing the auxilary fields DELTA and KSI of the CR method
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZP, ZQ
+ ! array containing the auxilary fields P and Q of the CR method
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZRESIDUE
+ ! array containing the error field at each iteration Q(PHI) - Y
+!
+REAL :: ZALPHA, ZLAMBDA ! amplitude of the descent in the Conjugate
+ ! directions
+REAL :: ZDOT_DELTA ! dot product of ZDELTA by itself
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!
+!* 1.1 compute the vector: r^(0) = Q(PHI) - Y
+!
+ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY
+!
+!* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y )
+!
+CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZP,&
+ PBFB,&
+ PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+!JUAN print*, "size ZP=",SIZE(ZP)
+!JUAN print*, "size ZRESIDUE=",SIZE(ZRESIDUE)
+!
+!* 1.3 compute the vector: delta^(0) = Q ( p^(0) )
+!
+ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. ITERATIVE LOOP
+! --------------
+!
+DO JM = 1,KITR
+!
+!* 2.1 compute the step LAMBDA
+!
+ ZDOT_DELTA = DOTPROD(ZDELTA, ZDELTA,HLBCX,HLBCY) ! norm of DELTA
+ ZLAMBDA = - DOTPROD(ZRESIDUE,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA
+!
+!* 2.2 update the pressure function PHI
+!
+ PPHI = PPHI + ZLAMBDA * ZP
+!
+!
+ IF( JM == KITR ) EXIT
+!
+!
+!* 2.3 update the residual error: r
+!
+ ZRESIDUE = ZRESIDUE + ZLAMBDA * ZDELTA
+!
+!* 2.4 compute the vector: q = F^(-1)*( Q(PHI) - Y )
+!
+CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZQ, &
+ PBFB,&
+ PBF_SXP2_YP1_Z) !JUAN Z_SPLITTING
+!
+!* 2.5 compute the auxiliary field: ksi = Q ( q )
+!
+ ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ)
+! -1
+!* 2.6 compute the step ALPHA
+!
+ ZALPHA = - DOTPROD(ZKSI,ZDELTA,HLBCX,HLBCY) / ZDOT_DELTA ! lambda
+!
+!* 2.7 update p and DELTA
+!
+ ZP = ZQ + ZALPHA * ZP
+ ZDELTA = ZKSI + ZALPHA * ZDELTA
+!
+END DO ! end of the loop for the iterative solver
+!
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE CONRESOLZ
diff --git a/src/ZSOLVER/ini_elecn.f90 b/src/ZSOLVER/ini_elecn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..27ed168cdda2b4f23db4406621869b9161434049
--- /dev/null
+++ b/src/ZSOLVER/ini_elecn.f90
@@ -0,0 +1,327 @@
+!MNH_LIC Copyright 2009-2019 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_INI_ELEC_n
+! ######################
+!
+INTERFACE
+ SUBROUTINE INI_ELEC_n (KLUOUT, HELEC, HCLOUD, TPINIFILE, &
+ PTSTEP, PZZ, &
+ PDXX, PDYY, PDZZ, PDZX, PDZY )
+!
+USE MODD_IO, ONLY : TFILEDATA
+!
+INTEGER, INTENT(IN) :: KLUOUT ! Logical unit number for prints
+CHARACTER (LEN=4), INTENT(IN) :: HELEC ! atmospheric electricity scheme
+CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! microphysics scheme
+TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE! Initial file
+REAL, INTENT(IN) :: PTSTEP ! Time STEP
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+END SUBROUTINE INI_ELEC_n
+END INTERFACE
+END MODULE MODI_INI_ELEC_n
+!
+! #########################################################
+ SUBROUTINE INI_ELEC_n(KLUOUT, HELEC, HCLOUD, TPINIFILE, &
+ PTSTEP, PZZ, &
+ PDXX, PDYY, PDZZ, PDZX, PDZY )
+! #########################################################
+!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to initialize the variables
+! of the atmospheric electricity scheme
+!
+!! METHOD
+!! ------
+!! The initialization of the scheme is performed as follows :
+!!
+!! EXTERNAL
+!! --------
+!! CLEANLIST_ll : deaalocate a list
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! C. Barthe * Laboratoire de l'Atmosphère et des Cyclones *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 09/11/09
+!! M. Chong 13/05/11 Add computation of specific parameters for solving
+!! the electric field equation (elements of tri-diag
+!! matrix)
+!! J.-P. Pinty 13/04/12 Add elec_trid to initialise the tridiagonal syst.
+!! J.-P. Pinty 01/07/12 Add a non-homogeneous Neuman fair-weather
+!! boundary condition at the top
+!! J.-P. Pinty 15/11/13 Initialize the flash maps
+!! 10/2016 (C.Lac) Add droplet deposition
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 14/02/2019: remove CLUOUT/CLUOUT0 and associated variables
+! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CLOUDPAR_n, ONLY : NSPLITR
+USE MODD_CONF, ONLY : CEQNSYS,CCONF,CPROGRAM
+USE MODD_CONF_n, ONLY : NRR
+USE MODD_CST
+USE MODD_DIM_n, ONLY : NIMAX_ll, NJMAX_ll
+USE MODD_DYN
+USE MODD_DYN_n, ONLY : XRHOM, XTRIGSX, XTRIGSY, XAF, XCF, XBFY, XBFB, XDXHATM, &
+ XDYHATM, NIFAXX, NIFAXY, XBF_SXP2_YP1_Z
+USE MODD_ELEC_DESCR
+USE MODD_ELEC_FLASH
+USE MODD_ELEC_n, ONLY : XRHOM_E, XAF_E, XCF_E, XBFY_E, XBFB_E, XBF_SXP2_YP1_Z_E
+USE MODD_GET_n, ONLY : CGETINPRC, CGETINPRR, CGETINPRS, CGETINPRG, CGETINPRH, &
+ CGETCLOUD, CGETSVT
+USE MODD_GRID_n, ONLY : XMAP, XDXHAT, XDYHAT
+USE MODD_IO, ONLY : TFILEDATA
+USE MODD_LBC_n, ONLY : CLBCX, CLBCY
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODD_PARAM_C2R2, ONLY : LDEPOC
+USE MODD_PARAMETERS, ONLY : JPVEXT, JPHEXT
+USE MODD_PARAM_ICE, ONLY : LDEPOSC
+USE MODD_PRECIP_n, ONLY : XINPRR, XACPRR, XINPRS, XACPRS, XINPRG, XACPRG, &
+ XINPRH, XACPRH, XINPRC, XACPRC, XINPRR3D, XEVAP3D,&
+ XINDEP,XACDEP
+USE MODD_REF
+USE MODD_REF_n, ONLY : XRHODJ, XTHVREF
+USE MODD_TIME
+!
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODE_ll
+use mode_msg
+!
+USE MODI_ELEC_TRIDZ
+USE MODI_INI_CLOUD
+USE MODI_INI_FIELD_ELEC
+USE MODI_INI_FLASH_GEOM_ELEC
+USE MODI_INI_PARAM_ELEC
+USE MODI_INI_RAIN_ICE_ELEC
+USE MODI_READ_PRECIP_FIELD
+!
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of dummy arguments
+!
+INTEGER, INTENT(IN) :: KLUOUT ! Logical unit number for prints
+CHARACTER (LEN=4), INTENT(IN) :: HELEC ! atmospheric electricity scheme
+CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! microphysics scheme
+TYPE(TFILEDATA), INTENT(IN) :: TPINIFILE! Initial file
+REAL, INTENT(IN) :: PTSTEP ! Time STEP
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height z
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: ILUOUT ! Logical unit number of output-listing
+!
+INTEGER :: IIU ! Upper dimension in x direction (local)
+INTEGER :: IJU ! Upper dimension in y direction (local)
+INTEGER :: IKU ! Upper dimension in z direction
+INTEGER :: IKB, IKE
+INTEGER :: JK ! Loop vertical index
+INTEGER :: IINFO_ll ! Return code of // routines
+INTEGER :: IINTVL ! Number of intervals to integrate the kernels
+REAL :: ZFDINFTY ! Factor used to define the "infinite" diameter
+!
+REAL :: ZRHO00 ! Surface reference air density
+REAL :: ZDZMIN
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDZ ! mesh size
+CHARACTER (LEN=3) :: YEQNSYS
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. PROLOGUE
+! --------
+!
+ILUOUT = TLUOUT%NLU
+!
+CALL GET_DIM_EXT_ll('B',IIU,IJU)
+IKU = SIZE(PZZ,3)
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. ALLOCATE Module MODD_PRECIP_n
+! -----------------------------
+!
+IF (HCLOUD(1:3) == 'ICE') THEN
+ ALLOCATE( XINPRR(IIU,IJU) )
+ ALLOCATE( XINPRR3D(IIU,IJU,IKU) )
+ ALLOCATE( XEVAP3D(IIU,IJU,IKU) )
+ ALLOCATE( XACPRR(IIU,IJU) )
+ XINPRR(:,:) = 0.0
+ XACPRR(:,:) = 0.0
+ XINPRR3D(:,:,:) = 0.0
+ XEVAP3D(:,:,:) = 0.0
+ ALLOCATE( XINPRC(IIU,IJU) )
+ ALLOCATE( XACPRC(IIU,IJU) )
+ XINPRC(:,:) = 0.0
+ XACPRC(:,:) = 0.0
+ ALLOCATE( XINPRS(IIU,IJU) )
+ ALLOCATE( XACPRS(IIU,IJU) )
+ XINPRS(:,:) = 0.0
+ XACPRS(:,:) = 0.0
+ ALLOCATE( XINPRG(IIU,IJU) )
+ ALLOCATE( XACPRG(IIU,IJU) )
+ XINPRG(:,:) = 0.0
+ XACPRG(:,:) = 0.0
+END IF
+!
+IF (HCLOUD == 'ICE4') THEN
+ ALLOCATE( XINPRH(IIU,IJU) )
+ ALLOCATE( XACPRH(IIU,IJU) )
+ XINPRH(:,:) = 0.0
+ XACPRH(:,:) = 0.0
+ELSE
+ ALLOCATE( XINPRH(0,0) )
+ ALLOCATE( XACPRH(0,0) )
+END IF
+!
+IF ( LDEPOSC) THEN
+ ALLOCATE(XINDEP(IIU,IJU))
+ ALLOCATE(XACDEP(IIU,IJU))
+ XINDEP(:,:)=0.0
+ XACDEP(:,:)=0.0
+ELSE
+ ALLOCATE(XINDEP(0,0))
+ ALLOCATE(XACDEP(0,0))
+END IF
+!
+IF(SIZE(XINPRR) == 0) RETURN
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. Initialize MODD_PRECIP_n variables
+! -----------------------------------
+!
+CALL READ_PRECIP_FIELD (TPINIFILE, CPROGRAM, CCONF, &
+ CGETINPRC,CGETINPRR,CGETINPRS,CGETINPRG,CGETINPRH, &
+ XINPRC,XACPRC,XINDEP,XACDEP,XINPRR,XINPRR3D,XEVAP3D, &
+ XACPRR, XINPRS, XACPRS, XINPRG, XACPRG, XINPRH, XACPRH)
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. INITIALIZE THE PARAMETERS
+!* FOR THE MICROPHYSICS AND THE ELECTRICITY
+! ----------------------------------------
+!
+!* 3.1 Compute the minimun vertical mesh size
+!
+ALLOCATE( ZDZ(IIU,IJU,IKU) )
+ZDZ(:,:,:) = 0.
+!
+IKB = 1 + JPVEXT
+IKE = SIZE(PZZ,3) - JPVEXT
+!
+DO JK = IKB, IKE
+ ZDZ(:,:,JK) = PZZ(:,:,JK+1) - PZZ(:,:,JK)
+END DO
+ZDZMIN = MIN_ll (ZDZ,IINFO_ll,1,1,IKB,NIMAX_ll+2*JPHEXT,NJMAX_ll+2*JPHEXT,IKE )
+!
+DEALLOCATE(ZDZ)
+!
+!
+IF (HELEC(1:3) == 'ELE') THEN
+!
+!
+!* 3.2 initialize the parameters for the mixed-phase microphysics
+!* and the electrification
+!
+ CALL INI_RAIN_ICE_ELEC (KLUOUT, PTSTEP, ZDZMIN, NSPLITR, HCLOUD, &
+ IINTVL, ZFDINFTY)
+!
+!
+!* 3.3 initialize the electrical parameters
+!
+ ZRHO00 = XP00 / (XRD * XTHVREFZ(IKB))
+!
+ CALL INI_PARAM_ELEC (TPINIFILE, CGETSVT, ZRHO00, NRR, IINTVL, &
+ ZFDINFTY, IIU, IJU, IKU)
+!
+!
+!* 3.4 initialize the parameters for the electric field
+!
+ IF (LINDUCTIVE .OR. ((.NOT. LOCG) .AND. LELEC_FIELD)) THEN
+ CALL INI_FIELD_ELEC (PDXX, PDYY, PDZZ, PDZX, PDZY, PZZ)
+ END IF
+!
+!
+!* 3.5 initialize the parameters for the lightning flashes
+!
+ IF (.NOT. LOCG) THEN
+ IF (LFLASH_GEOM) THEN
+ CALL INI_FLASH_GEOM_ELEC
+ ELSE
+ call Print_msg( NVERB_FATAL, 'GEN', 'INI_ELEC_n', 'INI_LIGHTNING_ELEC not yet developed' )
+ END IF
+ END IF
+!
+ELSE IF (HELEC /= 'NONE') THEN
+ call Print_msg( NVERB_FATAL, 'GEN', 'INI_ELEC_n', 'not yet developed for CELEC='//trim(HELEC) )
+END IF
+!
+!* 3.6 initialize the parameters for the resolution of the electric field
+!
+YEQNSYS = CEQNSYS
+CEQNSYS = 'LHE'
+! Force any CEQNSYS (DUR, MAE, LHE) to LHE to obtain a unique set of coefficients
+! for the flat laplacian operator and Return to the original CEQNSYS
+
+ALLOCATE (XRHOM_E(SIZE(XRHOM)))
+ALLOCATE (XAF_E(SIZE(XAF)))
+ALLOCATE (XCF_E(SIZE(XCF)))
+ALLOCATE (XBFY_E(SIZE(XBFY,1),SIZE(XBFY,2),SIZE(XBFY,3)))
+ALLOCATE (XBFB_E(SIZE(XBFB,1),SIZE(XBFB,2),SIZE(XBFB,3)))
+ALLOCATE (XBF_SXP2_YP1_Z_E(SIZE(XBF_SXP2_YP1_Z,1),SIZE(XBF_SXP2_YP1_Z,2),&
+ SIZE(XBF_SXP2_YP1_Z,3)))
+!
+CALL ELEC_TRIDZ (CLBCX,CLBCY, &
+ XMAP,XDXHAT,XDYHAT,XDXHATM,XDYHATM,XRHOM_E,XAF_E, &
+ XCF_E,XTRIGSX,XTRIGSY,NIFAXX,NIFAXY, &
+ XRHODJ,XTHVREF,PZZ,XBFY_E,XEPOTFW_TOP, &
+ XBFB_E,XBF_SXP2_YP1_Z_E)
+!
+CEQNSYS=YEQNSYS
+!
+!* 3.7 initialize the flash maps
+!
+ALLOCATE( NMAP_TRIG_IC(IIU,IJU) ); NMAP_TRIG_IC(:,:) = 0
+ALLOCATE( NMAP_IMPACT_CG(IIU,IJU) ); NMAP_IMPACT_CG(:,:) = 0
+ALLOCATE( NMAP_2DAREA_IC(IIU,IJU) ); NMAP_2DAREA_IC(:,:) = 0
+ALLOCATE( NMAP_2DAREA_CG(IIU,IJU) ); NMAP_2DAREA_CG(:,:) = 0
+ALLOCATE( NMAP_3DIC(IIU,IJU,IKU) ); NMAP_3DIC(:,:,:) = 0
+ALLOCATE( NMAP_3DCG(IIU,IJU,IKU) ); NMAP_3DCG(:,:,:) = 0
+!
+!-------------------------------------------------------------------------------
+!
+!
+END SUBROUTINE INI_ELEC_n
diff --git a/src/ZSOLVER/ini_field_elec.f90 b/src/ZSOLVER/ini_field_elec.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c5dcbb79a30d93761faecc41871205996584c396
--- /dev/null
+++ b/src/ZSOLVER/ini_field_elec.f90
@@ -0,0 +1,276 @@
+!MNH_LIC Copyright 2002-2020 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_INI_FIELD_ELEC
+! ##########################
+!
+INTERFACE
+!
+ SUBROUTINE INI_FIELD_ELEC (PDXX, PDYY, PDZZ, PDZX, PDZY, PZZ)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! vertical grid
+!
+END SUBROUTINE INI_FIELD_ELEC
+END INTERFACE
+END MODULE MODI_INI_FIELD_ELEC
+!
+! ############################################################
+ SUBROUTINE INI_FIELD_ELEC(PDXX, PDYY, PDZZ, PDZX, PDZY, PZZ)
+! ############################################################
+!
+!
+!!**** *INI_FIELD_ELEC* - routine to initialize the electric field
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to initialize the variables
+! of the electric field computation
+!
+!!** METHOD
+!! ------
+!! The initialization of the scheme is performed as follows :
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aérologie *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 29/11/02
+!! C. Barthe 10/11/09 phasage en version 4.8.1
+!! M. Chong 26/01/10 Small ions parameters
+!! + Fair weather field from Helsdon-Farley
+!! (JGR, 1987, 5661-5675)
+!! J.-P. Pinty 01/07/12 Add a non-homogeneous Neuman fair-weather
+!! boundary condition at the top
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!!
+!!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+USE MODD_CONF
+USE MODD_CST
+USE MODD_REF_n, ONLY: XRHODREF
+USE MODD_ELEC_DESCR
+USE MODD_ELEC_n
+!
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+!
+USE MODI_GDIV
+USE MODI_SHUMAN
+!
+USE MODE_ll
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of dummy arguments
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! Metric coefficients
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! vertical grid
+!
+!* 0.2 Declaration of local variables
+!
+!
+CHARACTER(LEN=4), DIMENSION(2) :: ZLBCX ! x-direction LBC type
+CHARACTER(LEN=4), DIMENSION(2) :: ZLBCY ! y-direction LBC type
+!
+INTEGER :: JK ! loop over the vertical levels
+INTEGER :: IINFO_ll !
+!
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZZMASS, ZWORK, ZWORK1, ZWORK2
+!
+TYPE(LIST_ll),POINTER :: TZFIELDS_ll ! list of fields to exchange
+!
+!
+!------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+ZLBCX = 'OPEN' ! forced LBC
+ZLBCY = 'OPEN' ! forced LBC
+!
+NULLIFY(TZFIELDS_ll)
+!
+! Allocations
+!
+ALLOCATE( XEFIELDU(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XEFIELDV(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XEFIELDW(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XESOURCEFW(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( ZZMASS(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) ) !Alt at mass point
+ALLOCATE( ZWORK(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( ZWORK1(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+IF( .NOT. LCOSMIC_APPROX ) THEN
+ ALLOCATE( ZWORK2(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+END IF
+ALLOCATE( XIONSOURCEFW(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XCION_POS_FW(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XCION_NEG_FW(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XMOBIL_POS(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+ALLOCATE( XMOBIL_NEG(SIZE(PDZZ,1),SIZE(PDZZ,2),SIZE(PDZZ,3)) )
+!
+!++ jpp
+ALLOCATE( XEPOTFW_TOP(SIZE(PDZZ,1),SIZE(PDZZ,2)) )
+!-- jpp
+!
+!------------------------------------------------------------------------------
+!
+!* 2. FAIR WEATHER ELECTRIC FIELD
+! ---------------------------
+!
+! The vertical component of the electric field is given by : E=E_0 * exp(k_e*z)
+! where E_0 = -100 V m^-1 and k_e = -292e-6 m^-1
+! We define the electric field as : E = - rhodJ * Nabla V
+!
+! Helsdon-Farley: E=E_0 (b1 exp(-a1 z) + b2 exp(-a2 z) + b3 exp(-a3 z)
+XEFIELDU(:,:,:) = 0.
+XEFIELDV(:,:,:) = 0.
+!
+! Initialization of Fair Weather Electric Field at W-point
+IF( .NOT. LFW_HELFA ) THEN
+ XEFIELDW(:,:,:) = XE_0 * EXP(XKEF * PZZ(:,:,:))
+ELSE
+ XEFIELDW(:,:,:) = XE0_HF * (XB1_HF*EXP(-XA1_HF*PZZ(:,:,:)) &
+ + XB2_HF*EXP(-XA2_HF*PZZ(:,:,:)) &
+ + XB3_HF*EXP(-XA3_HF*PZZ(:,:,:)))
+END IF
+!++ jpp
+XEPOTFW_TOP(:,:) = XEFIELDW(:,:,SIZE(PDZZ,3)) ! used in the top boundary
+ ! condition when inverting the
+ ! Gauss equation in V (here EPOT)
+!-- jpp
+XEFIELDW(:,:,SIZE(PDZZ,3)) = 2. * XEFIELDW(:,:,SIZE(PDZZ,3)-1) - &
+ XEFIELDW(:,:,SIZE(PDZZ,3)-2)
+
+! Computing the mobility of small positive (negative) ions at Mass-point
+ZZMASS = MZF( PZZ ) ! altitude at mass point
+
+DO JK = 2,SIZE(PZZ,3)-1
+ XMOBIL_POS(:,:,JK) = XF_POS * EXP( XEXPMOB* ZZMASS(:,:,JK) )
+ XMOBIL_NEG(:,:,JK) = XF_NEG * EXP( XEXPMOB* ZZMASS(:,:,JK) )
+END DO
+
+XMOBIL_POS(:,:,1) = 2.0 * XMOBIL_POS(:,:,2) - XMOBIL_POS(:,:,3)
+XMOBIL_POS(:,:,SIZE(PDZZ,3)) = 2. * XMOBIL_POS(:,:,SIZE(PDZZ,3)-1) - &
+ XMOBIL_POS(:,:,SIZE(PDZZ,3)-2)
+XMOBIL_NEG(:,:,1) = 2.0*XMOBIL_NEG(:,:,2) - XMOBIL_NEG(:,:,3)
+XMOBIL_NEG(:,:,SIZE(PDZZ,3)) = 2. * XMOBIL_NEG(:,:,SIZE(PDZZ,3)-1) - &
+ XMOBIL_NEG(:,:,SIZE(PDZZ,3)-2)
+!
+! Initial number concentrations of small positive (negative) free ions
+IF( .NOT. LFW_HELFA ) THEN
+ ZWORK(:,:,:) = XE_0 * EXP(XKEF * ZZMASS(:,:,:))
+ ZWORK1(:,:,:) = XE_0 * XKEF * EXP(XKEF * ZZMASS(:,:,:))
+ IF(.NOT. LCOSMIC_APPROX) THEN
+ ZWORK2(:,:,:) = XE_0 * XKEF * XKEF * EXP(XKEF * ZZMASS(:,:,:))
+ END IF
+ELSE
+ ZWORK(:,:,:)= XE0_HF * (XB1_HF*EXP(-XA1_HF*ZZMASS(:,:,:)) &
+ + XB2_HF*EXP(-XA2_HF*ZZMASS(:,:,:)) &
+ + XB3_HF*EXP(-XA3_HF*ZZMASS(:,:,:)))
+ ZWORK1(:,:,:)= XE0_HF * (-XB1_HF*XA1_HF*EXP(-XA1_HF*ZZMASS(:,:,:)) &
+ -XB2_HF*XA2_HF*EXP(-XA2_HF*ZZMASS(:,:,:)) &
+ -XB3_HF*XA3_HF*EXP(-XA3_HF*ZZMASS(:,:,:)))
+ IF(.NOT. LCOSMIC_APPROX) THEN
+ ZWORK2(:,:,:)= XE0_HF * (XB1_HF*XA1_HF*XA1_HF*EXP(-XA1_HF*ZZMASS(:,:,:)) &
+ +XB2_HF*XA2_HF*XA2_HF*EXP(-XA2_HF*ZZMASS(:,:,:)) &
+ +XB3_HF*XA3_HF*XA3_HF*EXP(-XA3_HF*ZZMASS(:,:,:)))
+ END IF
+END IF
+!
+XCION_POS_FW(:,:,:) = (XMOBIL_NEG(:,:,:) * XEPSILON * ZWORK1(:,:,:) + &
+ XJCURR_FW / ZWORK(:,:,:)) / &
+ (XECHARGE * (XMOBIL_POS(:,:,:) + XMOBIL_NEG(:,:,:)))
+XCION_NEG_FW(:,:,:) = XCION_POS_FW - XEPSILON * ZWORK1(:,:,:) / XECHARGE
+XCION_POS_FW(:,:,SIZE(PDZZ,3)) = 2. * XCION_POS_FW(:,:,SIZE(PDZZ,3)-1) - &
+ XCION_POS_FW(:,:,SIZE(PDZZ,3)-2)
+XCION_NEG_FW(:,:,SIZE(PDZZ,3)) = 2. * XCION_NEG_FW(:,:,SIZE(PDZZ,3)-1) - &
+ XCION_NEG_FW(:,:,SIZE(PDZZ,3)-2)
+!
+WHERE(XCION_NEG_FW < 0.) XCION_NEG_FW = 0.
+!
+! Computing the ion source from cosmic rays
+XIONSOURCEFW(:,:,:) = XIONCOMB * XCION_POS_FW(:,:,:) * XCION_NEG_FW(:,:,:)
+!
+IF ( .NOT. LCOSMIC_APPROX ) THEN
+ XIONSOURCEFW(:,:,:) = XIONSOURCEFW(:,:,:) + &
+ XMOBIL_POS(:,:,:) * XMOBIL_NEG(:,:,:) * XEPSILON * &
+ (XEXPMOB * ZWORK(:,:,:) * ZWORK1(:,:,:) + &
+ ZWORK1(:,:,:) * ZWORK1(:,:,:) + &
+ ZWORK(:,:,:) * ZWORK2(:,:,:)) / &
+ (XECHARGE * (XMOBIL_POS(:,:,:) + XMOBIL_NEG(:,:,:)))
+
+ XIONSOURCEFW(:,:,1) = 0.
+ XIONSOURCEFW(:,:,SIZE(PDZZ,3)) = 2. * XIONSOURCEFW(:,:,SIZE(PDZZ,3)-1) - &
+ XIONSOURCEFW(:,:,SIZE(PDZZ,3)-2)
+END IF
+!
+! Transform ion concentration into ion mixing ratio (Number/kg of air)
+
+XCION_POS_FW(:,:,:) = XCION_POS_FW(:,:,:) / XRHODREF(:,:,:)
+XCION_NEG_FW(:,:,:) = XCION_NEG_FW(:,:,:) / XRHODREF(:,:,:)
+XCION_POS_FW(:,:,SIZE(PDZZ,3)) = 2. * XCION_POS_FW(:,:,SIZE(PDZZ,3)-1) - &
+ XCION_POS_FW(:,:,SIZE(PDZZ,3)-2)
+XCION_NEG_FW(:,:,SIZE(PDZZ,3)) = 2. * XCION_NEG_FW(:,:,SIZE(PDZZ,3)-1) - &
+ XCION_NEG_FW(:,:,SIZE(PDZZ,3)-2)
+!
+XEFIELDW(:,:,1) = 0. ! Electric field null in a conductor
+XEFIELDW(:,:,SIZE(PDZZ,3)) = 0. ! either the ground or the ionosphere!
+!
+!
+!------------------------------------------------------------------------------
+!
+!* 3. FAIR WEATHER SPACE CHARGE
+! -------------------------
+!
+CALL ADD3DFIELD_ll( TZFIELDS_ll, XEFIELDU, 'INI_FIELD_ELEC::XEFIELDU' )
+CALL ADD3DFIELD_ll( TZFIELDS_ll, XEFIELDV, 'INI_FIELD_ELEC::XEFIELDV' )
+CALL ADD3DFIELD_ll( TZFIELDS_ll, XEFIELDW, 'INI_FIELD_ELEC::XEFIELDW' )
+CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+CALL GDIV (ZLBCX, ZLBCY, &
+ PDXX, PDYY, PDZX, PDZY, PDZZ, &
+ XEFIELDU,XEFIELDV,XEFIELDW, &
+ XESOURCEFW )
+!
+XESOURCEFW(:,:,:) = XESOURCEFW(:,:,:) * XEPSILON ! Nabla E * epsilon = + rho
+ ! C / m^3
+XESOURCEFW(:,:,SIZE(PZZ,3)) = XESOURCEFW(:,:,SIZE(PZZ,3)-1)
+!
+DEALLOCATE(ZZMASS)
+DEALLOCATE(ZWORK)
+DEALLOCATE(ZWORK1)
+IF( .NOT. LCOSMIC_APPROX) THEN
+ DEALLOCATE(ZWORK2)
+END IF
+!
+!
+!------------------------------------------------------------------------------
+!
+END SUBROUTINE INI_FIELD_ELEC
diff --git a/src/ZSOLVER/lap_m.f90 b/src/ZSOLVER/lap_m.f90
new file mode 100644
index 0000000000000000000000000000000000000000..f1936c828237a3c61b7f54fcb0ab6ee86376396a
--- /dev/null
+++ b/src/ZSOLVER/lap_m.f90
@@ -0,0 +1,228 @@
+!MNH_LIC Copyright 2007-2020 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_LAP_M
+! #################
+!
+INTERFACE
+!
+ FUNCTION LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PY) &
+ RESULT(PLAP_M)
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density_reference * J
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components
+!
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PLAP_M ! final divergence
+!
+END FUNCTION LAP_M
+!
+END INTERFACE
+!
+END MODULE MODI_LAP_M
+!
+!
+!
+! #########################################################################
+ FUNCTION LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PY) &
+ RESULT(PLAP_M)
+! #########################################################################
+!
+!!**** *LAP_M * - compute the Laplacian of a field PY
+!!
+!! PURPOSE
+!! -------
+! This function computes laplacian of a scalar field PY
+! localized at mass points, with bottom topography.
+! The result is localized at a mass point and defined by:
+! ( ( GX_M_U (PY) ) )
+! PLAP_M = GDIV ( rho * J ( GX_M_V (PY) ) )
+! ( ( GX_M_W (PY) ) )
+!
+! Where GX_M_.. are the cartesian components of the gradient of PY and
+! GDIV is the operator acting on a vector AA:
+!
+! GDIV ( AA ) = J * divergence (1/J AA )
+!
+!!** METHOD
+!! ------
+!! First, we compute the gradients along x, y , z of the P field. The
+!! result is multiplied by rhod.
+!! Then, the pseudo-divergence ( J * DIV (1/J o ) ) is computed by the
+!! subroutine GDIV. The result is localized at a mass point.
+!!
+!! EXTERNAL
+!! --------
+!! 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
+!! 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 GDIV : compute J times the divergence of 1/J times a vector
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_PARAMETERS: JPHEXT, JPVEXT
+!! Module MODD_CONF: L2D
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! P. Hereil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 01/02/07 simplified from QLAP function, T.Maric
+!! Modification
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODE_ll
+!
+USE MODD_PARAMETERS
+USE MODD_CONF
+!USE MODD_CST
+USE MODI_GDIV
+!USE MODI_GDIV_M
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of 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
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components
+!
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PLAP_M ! final divergence
+!
+!* 0.2 declarations of local variables
+!
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZU ! rho*J*gradient along x
+!
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZV ! rho*J*gradient along y
+!
+REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZW ! rho*J*gradient along z
+!
+INTEGER :: IIU,IJU,IKU ! I,J,K array sizes
+INTEGER :: JK,JJ,JI ! vertical loop index
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. COMPUTE THE GRADIENT COMPONENTS
+! -------------------------------
+!
+!
+CALL GET_DIM_EXT_ll('B',IIU,IJU)
+IKU=SIZE(PY,3)
+!
+ZU = GX_M_U(1,IKU,1,PY,PDXX,PDZZ,PDZX)
+!
+IF ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() ) THEN
+ DO JK=2,IKU-1
+ DO JJ=1,IJU
+ ZU(2,JJ,JK)= (PY(2,JJ,JK) - PY(1,JJ,JK) - 0.5 * ( &
+ PDZX(2,JJ,JK) * (PY(2,JJ,JK)-PY(2,JJ,JK-1)) / PDZZ(2,JJ,JK) &
+ +PDZX(2,JJ,JK+1) * (PY(2,JJ,JK+1)-PY(2,JJ,JK)) / PDZZ(2,JJ,JK+1) &
+ ) ) / PDXX(2,JJ,JK)
+ END DO
+ END DO
+END IF
+!
+IF ( HLBCX(1) /= 'CYCL' .AND. LEAST_ll() ) THEN
+ DO JK=2,IKU-1
+ DO JJ=1,IJU
+ ZU(IIU,JJ,JK)= (PY(IIU,JJ,JK) - PY(IIU-1,JJ,JK) - 0.5 * ( &
+ PDZX(IIU,JJ,JK) * (PY(IIU-1,JJ,JK)-PY(IIU-1,JJ,JK-1)) / PDZZ(IIU-1,JJ,JK) &
+ +PDZX(IIU,JJ,JK+1) * (PY(IIU-1,JJ,JK+1)-PY(IIU-1,JJ,JK)) / PDZZ(IIU-1,JJ,JK+1)&
+ ) ) / PDXX(IIU,JJ,JK)
+ END DO
+ END DO
+END IF
+!
+IF(.NOT. L2D) THEN
+!
+ ZV = GY_M_V(1,IKU,1,PY,PDYY,PDZZ,PDZY)
+!
+ IF ( HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() ) THEN
+ DO JK=2,IKU-1
+ DO JI=1,IIU
+ ZV(JI,2,JK)= (PY(JI,2,JK) - PY(JI,1,JK) - 0.5 * ( &
+ PDZY(JI,2,JK) * (PY(JI,2,JK)-PY(JI,2,JK-1)) / PDZZ(JI,2,JK) &
+ +PDZY(JI,2,JK+1) * (PY(JI,2,JK+1)-PY(JI,2,JK)) / PDZZ(JI,2,JK+1) &
+ ) ) / PDYY(JI,2,JK)
+ END DO
+ END DO
+ END IF
+ IF ( HLBCY(1) /= 'CYCL' .AND. LNORTH_ll() ) THEN
+!
+ DO JK=2,IKU-1
+ DO JI=1,IIU
+ ZV(JI,IJU,JK)= (PY(JI,IJU,JK) - PY(JI,IJU-1,JK) - 0.5 * ( &
+ PDZY(JI,IJU,JK) * (PY(JI,IJU-1,JK)-PY(JI,IJU-1,JK-1)) / PDZZ(JI,IJU-1,JK) &
+ +PDZY(JI,IJU,JK+1) * (PY(JI,IJU-1,JK+1)-PY(JI,IJU-1,JK)) / PDZZ(JI,IJU-1,JK+1)&
+ ) ) / PDYY(JI,IJU,JK)
+ END DO
+ END DO
+ END IF
+!
+ELSE
+ ZV=0.
+ENDIF
+!
+ZU = MXM(PRHODJ) * ZU
+!
+IF(.NOT. L2D) THEN
+ ZV = MYM(PRHODJ) * ZV
+ENDIF
+!
+ZW = MZM(PRHODJ) * GZ_M_W(1,IKU,1,PY,PDZZ)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTE THE DIVERGENCE
+! ----------------------
+!
+CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PLAP_M)
+!
+PLAP_M(:,:,1) = 0.
+!
+!-------------------------------------------------------------------------------
+!
+END FUNCTION LAP_M
diff --git a/src/ZSOLVER/pressure.f90 b/src/ZSOLVER/pressure.f90
new file mode 100644
index 0000000000000000000000000000000000000000..f06d79f0a5b295008267b2929473caa4db73f913
--- /dev/null
+++ b/src/ZSOLVER/pressure.f90
@@ -0,0 +1,687 @@
+!MNH_LIC Copyright 1994-2020 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_PRESSURE
+!###################
+!
+INTERFACE
+!
+ SUBROUTINE PRESSURE( &
+ HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI, &
+ PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOM, &
+ PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PPABSM, &
+ KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0, &
+ PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS, &
+ PRUS,PRVS,PRWS,PPABST,PRESIDUAL)
+!
+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,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) :: PRHOM ! 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.
+ ! 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
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABSM ! pressure (t-dt)
+!
+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
+REAL, OPTIONAL :: PRESIDUAL
+!JUAN
+!
+END SUBROUTINE PRESSURE
+!
+END INTERFACE
+!
+END MODULE MODI_PRESSURE
+! ######################################################################
+ SUBROUTINE PRESSURE( &
+ HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI, &
+ PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOM, &
+ PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PPABSM, &
+ KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0, &
+ PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS, &
+ PRUS,PRVS,PRWS,PPABST,PRESIDUAL)
+! ######################################################################
+!
+!!**** *PRESSURE * - 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
+!! 06/2011 (J.escobar ) Bypass Bug with ifort11/12 on HLBCX,HLBCY
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAMETERS
+USE MODD_BUDGET
+USE MODD_CONF
+USE MODD_CST
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODI_MASS_LEAK
+USE MODI_GDIV
+USE MODI_FLAT_INV
+USE MODI_RICHARDSON
+USE MODI_CONJGRAD
+USE MODI_CONRESOL
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN
+USE MODI_P_ABS
+USE MODI_BUDGET
+!
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODE_ll
+!
+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,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) :: PRHOM ! 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.
+ ! 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
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABSM ! pressure (t-dt)
+!
+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
+REAL, OPTIONAL :: PRESIDUAL
+!JUAN
+!
+!
+!* 0.2 declarations of local variables
+!
+! Metric coefficients:
+!
+REAL, DIMENSION(SIZE(PPABSM,1),SIZE(PPABSM,2),SIZE(PPABSM,3)) :: 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
+!
+REAL, DIMENSION(SIZE(PPABSM,1),SIZE(PPABSM,2),SIZE(PPABSM,3)) :: ZTHETAV, &
+ ! virtual potential temperature
+ ZPHIT
+ ! MAE + DUR => Exner function perturbation
+ ! LHE => Exner function perturbation * CPD * THVREF
+!
+REAL :: ZRV_OV_RD ! XRV / XRD
+REAL :: ZMAXVAL, ZMAXRES ! 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, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,3)) :: ZPABS_S ! local pressure on southern side
+REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,3)) :: ZPABS_N ! local pressure on northern side
+REAL, DIMENSION(SIZE(PDYY,2),SIZE(PDXX,3)) :: ZPABS_E ! local pressure on eastern side
+REAL, DIMENSION(SIZE(PDYY,2),SIZE(PDXX,3)) :: ZPABS_W ! local pressure on western side
+INTEGER :: IINFO_ll
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+TYPE(LIST_ll), POINTER :: TZFIELDS_2_ll ! list of fields to exchange
+!
+!
+!------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
+NULLIFY(TZFIELDS_ll)
+NULLIFY(TZFIELDS_2_ll)
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+ILUOUT = TLUOUT%NLU
+!
+CALL GET_PHYSICAL_ll(IIB,IJB,IIE,IJE)
+CALL GET_DIM_EXT_ll('B',IIU,IJU)
+!
+IKB= 1+JPVEXT
+IKU= SIZE(PPABSM,3)
+IKE= IKU - JPVEXT
+!
+ZPABS_S(:,:) = 0.
+ZPABS_N(:,:) = 0.
+ZPABS_E(:,:) = 0.
+ZPABS_W(:,:) = 0.
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 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
+! --------------------------------------------------
+!
+!
+CALL ADD3DFIELD_ll( TZFIELDS_ll, PRUS, 'PRESSURE::PRUS' )
+CALL ADD3DFIELD_ll( TZFIELDS_ll, PRVS, 'PRESSURE::PRVS' )
+CALL ADD3DFIELD_ll( TZFIELDS_ll, PRWS, 'PRESSURE::PRWS' )
+CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
+!
+! The non-homogenous Neuman problem is transformed in an homogenous Neuman
+! problem in the non-periodic cases
+IF (HLBCX(1) /= 'CYCL') THEN
+ IF (LWEST_ll()) ZDV_SOURCE(IIB-1,:,:) = 0.
+ IF (LEAST_ll()) ZDV_SOURCE(IIE+1,:,:) = 0.
+ENDIF
+!
+IF (.NOT. L2D .AND. HLBCY(1) /= 'CYCL') THEN
+ IF (LSOUTH_ll()) ZDV_SOURCE(:,IJB-1,:) = 0.
+ IF (LNORTH_ll()) ZDV_SOURCE(:,IJE+1,:) = 0.
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+!* 5. SOLVE THE PRESSURE EQUATION
+! ---------------------------
+!
+!
+!* 5.1 Compute the virtual theta and the pressure perturbation
+! -------------------------------------------------------
+!
+IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+ IF(KRR > 0) THEN
+ !
+ ! compute the ratio : 1 + total water mass / dry air mass
+ ZRV_OV_RD = XRV / XRD
+ ZTHETAV(:,:,:) = 1. + PRT(:,:,:,1)
+ 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
+ !
+ ZPHIT(:,:,:)=(PPABSM(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
+ !
+ELSEIF(CEQNSYS=='LHE') THEN
+ ZPHIT(:,:,:)= ((PPABSM(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) &
+ * XCPD * PTHVREF(:,:,:)
+ !
+END IF
+!
+IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN) THEN
+ ! flat cartesian LHE case -> exact solution
+ !
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZDV_SOURCE,ZPHIT)
+ELSE
+ SELECT CASE(HPRESOPT)
+ CASE('RICHA') ! Richardson's method
+!
+ CALL RICHARDSON(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
+ PDXHATM,PDYHATM,PRHOM,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,PRHOM,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,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
+ KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
+ END SELECT
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 6. ADD THE PRESSURE GRADIENT TO THE SOURCES
+! ----------------------------------------
+!
+IF ( HLBCX(1) /= 'CYCL' ) THEN
+ IF(LWEST_ll()) ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
+ IF(LEAST_ll()) ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
+ENDIF
+IF ( HLBCY(1) /= 'CYCL' ) THEN
+ IF (LSOUTH_ll()) ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
+ IF (LNORTH_ll()) ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
+ENDIF
+!
+IF ( L2D ) THEN
+ IF (LSOUTH_ll()) ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)
+ IF (LNORTH_ll()) ZPHIT(:,IJE+1,:) = ZPHIT(:,IJB,:)
+END IF
+!
+ZDV_SOURCE = GX_M_U(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX)
+!
+IF ( HLBCX(1) /= 'CYCL' ) THEN
+ IF(LWEST_ll()) THEN
+!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+ DO JK=2,IKU-1
+ ZDV_SOURCE(2,:,JK)= &
+ (ZPHIT(2,:,JK) - ZPHIT(1,:,JK) - 0.5 * ( &
+ PDZX(2,:,JK) * (ZPHIT(2,:,JK)-ZPHIT(2,:,JK-1)) / PDZZ(2,:,JK) &
+ +PDZX(2,:,JK+1) * (ZPHIT(2,:,JK+1)-ZPHIT(2,:,JK)) / PDZZ(2,:,JK+1) &
+ ) &
+ ) / PDXX(2,:,JK)
+ END DO
+ ENDIF
+ !
+ IF(LEAST_ll()) THEN
+ DO JK=2,IKU-1
+ ZDV_SOURCE(IIU,:,JK)= &
+ (ZPHIT(IIU,:,JK) - ZPHIT(IIU-1,:,JK) - 0.5 * ( &
+ PDZX(IIU,:,JK) * (ZPHIT(IIU-1,:,JK)-ZPHIT(IIU-1,:,JK-1)) &
+ / PDZZ(IIU-1,:,JK) &
+ +PDZX(IIU,:,JK+1) * (ZPHIT(IIU-1,:,JK+1)-ZPHIT(IIU-1,:,JK)) &
+ / PDZZ(IIU-1,:,JK+1) &
+ ) &
+ ) / PDXX(IIU,:,JK)
+ END DO
+ END IF
+END IF
+!
+IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+ PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
+ PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
+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
+!
+ ZDV_SOURCE = GY_M_V(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY)
+!
+ IF ( HLBCY(1) /= 'CYCL' ) THEN
+ IF (LSOUTH_ll()) THEN
+!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
+ DO JK=2,IKU-1
+ ZDV_SOURCE(:,2,JK)= &
+ (ZPHIT(:,2,JK) - ZPHIT(:,1,JK) - 0.5 * ( &
+ PDZY(:,2,JK) * (ZPHIT(:,2,JK)-ZPHIT(:,2,JK-1)) / PDZZ(:,2,JK) &
+ +PDZY(:,2,JK+1) * (ZPHIT(:,2,JK+1)-ZPHIT(:,2,JK)) / PDZZ(:,2,JK+1) &
+ ) &
+ ) / PDYY(:,2,JK)
+ END DO
+ END IF
+ !
+ IF (LNORTH_ll()) THEN
+ DO JK=2,IKU-1
+ ZDV_SOURCE(:,IJU,JK)= &
+ (ZPHIT(:,IJU,JK) - ZPHIT(:,IJU-1,JK) - 0.5 * ( &
+ PDZY(:,IJU,JK) * (ZPHIT(:,IJU-1,JK)-ZPHIT(:,IJU-1,JK-1)) &
+ / PDZZ(:,IJU-1,JK) &
+ +PDZY(:,IJU,JK+1) * (ZPHIT(:,IJU-1,JK+1)-ZPHIT(:,IJU-1,JK)) &
+ / PDZZ(:,IJU-1,JK+1) &
+ ) &
+ ) / PDYY(:,IJU,JK)
+ END DO
+ END IF
+ END IF
+!
+ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+ PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
+ 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)
+DO JJ = IJB,IJE ! copy the horizontal components under
+ DO JI = IIB,IIE
+ PRUS(JI,JJ,IKB-1)=PRUS(JI,JJ,IKB)
+ PRUS(JI,JJ,IKE+1)=PRUS(JI,JJ,IKE)
+ END DO
+END DO
+!
+DO JJ = IJB,IJE
+ DO JI = IIB,IIE ! the ground and above the top
+ PRVS(JI,JJ,IKB-1)=PRVS(JI,JJ,IKB)
+ PRVS(JI,JJ,IKE+1)=PRVS(JI,JJ,IKE)
+ END DO
+END DO
+!!
+!
+! compute the residual divergence
+CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
+!
+IF ( CEQNSYS=='DUR' ) THEN
+ IF ( SIZE(PRVREF,1) == 0 ) THEN
+ ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF
+ ELSE
+ ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF*(1.+PRVREF)
+ END IF
+ELSEIF( CEQNSYS=='MAE' .OR. CEQNSYS=='LHE' ) THEN
+ ZDV_SOURCE=ZDV_SOURCE/PRHODJ*PRHODREF
+END IF
+!
+ZMAXVAL=MAX_ll(ABS(ZDV_SOURCE),IINFO_ll)
+IF (PRESENT(PRESIDUAL)) PRESIDUAL = ZMAXVAL
+IMAXLOC=MAXLOC( ABS(ZDV_SOURCE(IIB:IIE,IJB:IJE,IKB:IKE))) !provisory coding one one processor only
+!
+WRITE(ILUOUT,*) 'residual divergence / 2 DT', ZMAXVAL, &
+ ' located at ', IMAXLOC
+! number of iterations adjusted
+IF (LFLAT .AND. LCARTESIAN) THEN
+ ZMAXRES = 1.E-12
+ELSE
+ ZMAXRES = 1.E-9
+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
+!
+!* 7. STORAGE OF THE FIELDS IN BUDGET ARRAYS
+! --------------------------------------
+!
+IF (LBUDGET_U) CALL BUDGET (PRUS,1,'PRES_BU_RU')
+IF (LBUDGET_V) CALL BUDGET (PRVS,2,'PRES_BU_RV')
+IF (LBUDGET_W) CALL BUDGET (PRWS,3,'PRES_BU_RW')
+!
+!-------------------------------------------------------------------------------
+!
+!* 8. ABSOLUTE PRESSURE COMPUTATION
+! -----------------------------
+!
+!IF ( ABS(PRHODREF(IIB,IJB,IKB)-PRHODREF(IIB,IJB,IKE)) > 1.E-16 &
+IF ( ABS(PRHODREF(IIB,IJB,IKB)-PRHODREF(IIB,IJB,IKE)) > 1.E-12 &
+ .AND. KTCOUNT >0 ) THEN
+ CALL P_ABS ( KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
+ PTHT, PRT, PRHODJ, PRHODREF, ZTHETAV, PTHVREF, &
+ PRVREF, PEXNREF, ZPHIT )
+!
+ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
+ PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
+ ELSEIF(CEQNSYS=='LHE') THEN
+ PPABST(:,:,:)=XP00*(ZPHIT/(XCPD*PTHVREF)+PEXNREF)**(XCPD/XRD)
+ ENDIF
+!
+ IF( HLBCX(1) == 'CYCL' ) THEN
+ IF (LWEST_ll()) THEN
+ ZPABS_W(:,:)= PPABST(IIB,:,:)
+ END IF
+!
+ IF (LEAST_ll()) THEN
+ ZPABS_E(:,:)= PPABST(IIE+1,:,:)
+ END IF
+!
+ END IF
+!
+ IF( HLBCY(1) == 'CYCL' ) THEN
+ IF (LSOUTH_ll()) THEN
+ ZPABS_S(:,:)= PPABST(:,IJB,:)
+ END IF
+!
+ IF (LNORTH_ll()) THEN
+ ZPABS_N(:,:)= PPABST(:,IJE+1,:)
+ END IF
+!
+ END IF
+!
+ CALL ADD3DFIELD_ll( TZFIELDS_2_ll, PPABST, 'PRESSURE::PPABST' )
+ CALL UPDATE_HALO_ll(TZFIELDS_2_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_2_ll)
+!
+ IF( HLBCX(1) == 'CYCL' ) THEN
+ IF (LWEST_ll()) THEN
+ PPABST(IIB,:,:) = ZPABS_W(:,:)
+ END IF
+!
+ IF (LEAST_ll()) THEN
+ PPABST(IIE+1,:,:) = ZPABS_E(:,:)
+ END IF
+!
+ END IF
+!
+ IF( HLBCY(1) == 'CYCL' ) THEN
+ IF (LSOUTH_ll()) THEN
+ PPABST(:,IJB,:) = ZPABS_S(:,:)
+ END IF
+!
+ IF (LNORTH_ll()) THEN
+ PPABST(:,IJE+1,:) = ZPABS_N(:,:)
+ END IF
+!
+ END IF
+!
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE PRESSURE
diff --git a/src/ZSOLVER/pressure_in_prep.f90 b/src/ZSOLVER/pressure_in_prep.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6219e352f27f06e8a418d8073fdd8aed4031802d
--- /dev/null
+++ b/src/ZSOLVER/pressure_in_prep.f90
@@ -0,0 +1,308 @@
+!MNH_LIC Copyright 1998-2020 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_PRESSURE_IN_PREP
+! ######################
+!
+INTERFACE
+!
+ SUBROUTINE PRESSURE_IN_PREP(PDXX,PDYY,PDZX,PDZY,PDZZ)
+!
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+END SUBROUTINE PRESSURE_IN_PREP
+!
+END INTERFACE
+!
+END MODULE MODI_PRESSURE_IN_PREP
+!
+! #####################################################
+ SUBROUTINE PRESSURE_IN_PREP(PDXX,PDYY,PDZX,PDZY,PDZZ)
+! #####################################################
+!
+!!**** *PRESSURE_IN_PREP* - calls the pressure solver in prep_real_case and
+!! checks the result
+!!
+!! PURPOSE
+!! -------
+!!
+!!
+!!
+!!** METHOD
+!! ------
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!!
+!! Book 2
+!!
+!! AUTHOR
+!! ------
+!!
+!! V.Masson Meteo-France
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 22/12/98
+!! parallelization 18/06/00 (Jabouille)
+!! 2014 M.Faivre
+!! 08/2015 M.Moge removing UPDATE_HALO_ll on XUT, XVT, XRHODJ in part 4
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODE_ll
+USE MODE_MSG
+!
+USE MODI_ANEL_BALANCE_n
+USE MODI_GDIV
+USE MODI_SHUMAN
+!
+USE MODD_CONF ! declaration modules
+USE MODD_CONF_n
+USE MODD_LUNIT
+USE MODD_DIM_n
+USE MODD_GRID_n
+USE MODD_LBC_n
+USE MODD_PARAMETERS
+USE MODD_FIELD_n, ONLY: XUT,XVT,XWT
+USE MODD_DYN_n
+USE MODD_REF_n
+USE MODD_CST
+USE MODE_MPPDB
+USE MODE_EXTRAPOL
+!
+IMPLICIT NONE
+!
+!* 0.1 Declaration of dummy arguments
+! ------------------------------
+!
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDXX ! metric coefficient dxx
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDYY ! metric coefficient dyy
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZX ! metric coefficient dzx
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZY ! metric coefficient dzy
+REAL,DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! metric coefficient dzz
+!
+!* 0.2 Declaration of local variables
+! ------------------------------
+!
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZU ! U
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZV ! V
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZW ! W
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZRU ! U * rho * J
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZRV ! V * rho * J
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZRW ! W * rho * J
+REAL,DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: ZDIV ! residual divergence
+!
+!* file management variables and counters
+!
+INTEGER :: ILUOUT0 ! logical unit for listing file
+INTEGER :: IINFO_ll
+REAL :: ZMAXRES
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+REAL :: ZMAXVAL,ZRESIDUAL
+INTEGER, DIMENSION(3) :: IMAXLOC
+INTEGER :: I,J,K
+!-------------------------------------------------------------------------------
+!
+!* 1. Initialisations
+! ---------------
+!
+ILUOUT0 = TLUOUT0%NLU
+!
+ZU(:,:,:) = XUT(:,:,:)
+ZV(:,:,:) = XVT(:,:,:)
+ZW(:,:,:) = XWT(:,:,:)
+!
+NULLIFY(TZFIELDS_ll)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. Loop
+! ----
+!
+DO
+ XUT(:,:,:) = ZU(:,:,:)
+ XVT(:,:,:) = ZV(:,:,:)
+ XWT(:,:,:) = ZW(:,:,:)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, XUT, 'PRESSURE_IN_PREP::XUT' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, XVT, 'PRESSURE_IN_PREP::XVT' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, XWT, 'PRESSURE_IN_PREP::XWT' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+
+ CALL MPPDB_CHECK3D(XUT,"PREP::XUM",PRECISION)
+ CALL MPPDB_CHECK3D(XVT,"PREP::XVM",PRECISION)
+ CALL MPPDB_CHECK3D(XWT,"PREP::XWM",PRECISION)
+ CALL MPPDB_CHECK3D(XRHODJ,"PREP::XRHODJ",PRECISION)
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. ANELASTIC CORRECTION
+! --------------------
+!
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) 'CPRESOPT = ',CPRESOPT
+ WRITE(ILUOUT0,*) 'NITR = ',NITR
+ IF (CPRESOPT=='RICHA') &
+ WRITE(ILUOUT0,*) 'XRELAX = ',XRELAX
+!
+ CALL ANEL_BALANCE_n(ZRESIDUAL)
+!
+!-------------------------------------------------------------------------------
+!
+!* 4. compute the residual divergence
+! -------------------------------
+!20140225 forgot this update_halo
+!20131112 check 1st XUT
+CALL MPPDB_CHECK3D(XUT,"PressInP4-beforeupdhalo::XUT",PRECISION)
+CALL MPPDB_CHECK3D(XVT,"PressInP4-beforeupdhalo::XVT",PRECISION)
+!CALL ADD3DFIELD_ll( TZFIELDS_ll, XUT, 'PRESSURE_IN_PREP::XUT' )
+!CALL ADD3DFIELD_ll( TZFIELDS_ll, XVT, 'PRESSURE_IN_PREP::XVT' )
+!CALL ADD3DFIELD_ll( TZFIELDS_ll, XRHODJ, 'PRESSURE_IN_PREP::XRHODJ' )
+! CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+! CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+ ZRU(:,:,:) = XUT(:,:,:) * MXM(XRHODJ)
+ ZRV(:,:,:) = XVT(:,:,:) * MYM(XRHODJ)
+ ZRW(:,:,:) = XWT(:,:,:) * MZM(XRHODJ)
+!
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRU, 'PRESSURE_IN_PREP::ZRU' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRV, 'PRESSURE_IN_PREP::ZRV' )
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRW, 'PRESSURE_IN_PREP::ZRW' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+ CALL GDIV(CLBCX,CLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZRU,ZRV,ZRW,ZDIV)
+CALL MPPDB_CHECK3D(XUT,"PressInP4-afterupdhalo::XUT",PRECISION)
+CALL MPPDB_CHECK3D(XVT,"PressInP4-afterupdhalo::XVT",PRECISION)
+!
+!20131125 add extrapol on ZRU
+CALL EXTRAPOL('W',ZRU)
+CALL MPPDB_CHECK3D(ZRU,"PressInP4-afterextrapol W::ZRU",PRECISION)
+!
+!20131126 add extrapol on ZRV
+CALL EXTRAPOL('S',ZRV)
+!
+ IF ( CEQNSYS=='DUR' ) THEN
+ IF ( SIZE(XRVREF,1) == 0 ) THEN
+ ZDIV=ZDIV/XRHODJ/XTH00*XRHODREF*XTHVREF
+ ELSE
+ ZDIV=ZDIV/XRHODJ/XTH00*XRHODREF*XTHVREF*(1.+XRVREF)
+ END IF
+ ELSEIF( CEQNSYS=='MAE' .OR. CEQNSYS=='LHE' ) THEN
+ ZDIV=ZDIV/XRHODJ*XRHODREF
+ END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 5. modifies the solver parameters if necessary
+! -------------------------------------------
+!
+ IF (LRES) THEN
+ ZMAXRES = XRES
+ ELSE
+ ZMAXRES = XRES_PREP
+ END IF
+!
+ ZMAXVAL = ZRESIDUAL
+ IF ( ZMAXVAL > ZMAXRES) THEN
+!!$ IF (LRES) THEN
+!!$ ZMAXRES = XRES
+!!$ ELSE
+!!$ ZMAXRES = 1.E-08
+!!$ END IF
+!!$!
+!!$ IF ( MAX_ll( ABS(ZDIV),IINFO_ll) > ZMAXRES ) THEN
+ IF (CPRESOPT=='RICHA') THEN
+ IF (XRELAX>0.75) THEN
+ XRELAX = XRELAX - 0.1
+ ELSE
+ XRELAX = 1.
+ NITR = NITR + 4
+ END IF
+ ELSE
+ NITR = NITR + 4
+ END IF
+!
+ IF (NITR>40) THEN
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) '* *'
+ IF (CPROGRAM=='REAL ') WRITE(ILUOUT0,*) &
+ '* WARNING in PREP_REAL_CASE *'
+ IF (CPROGRAM=='IDEAL ') WRITE(ILUOUT0,*) &
+ '* WARNING in PREP_IDEAL_CASE *'
+ WRITE(ILUOUT0,*) '* *'
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) 'The pressure solver was unable to converge for your case.'
+ WRITE(ILUOUT0,*) 'This can be due to : '
+ WRITE(ILUOUT0,*) ' a locally very steep orography'
+ WRITE(ILUOUT0,*) ' a too high vertical stretching'
+ WRITE(ILUOUT0,*) ' a too high horizontal stretching on the sphere (large domains)'
+ WRITE(ILUOUT0,*) ' a too high horizontal stretching in conformal plane'
+ WRITE(ILUOUT0,*) ' an error in your input velocity and thermodynamical fields'
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) 'The model will probably not be able to run with this initial or coupling file.'
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) '******************************************************************************'
+ WRITE(ILUOUT0,*) ' '
+ !callabortstop
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','PRESSURE_IN_PREP','')
+ END IF
+ ELSE
+!* 7. Happy conclusion
+! ----------------
+!
+ WRITE(ILUOUT0,*) ' '
+ IF (.NOT. LCARTESIAN) THEN
+ WRITE(ILUOUT0,*) 'Horizontal stretching in conformal plane: '
+ WRITE(ILUOUT0,*) ' map factor varies between ', MINVAL(XMAP),' and ', MAXVAL(XMAP)
+ WRITE(ILUOUT0,*) ' '
+ ENDIF
+ WRITE(ILUOUT0,*) '***************************************************'
+ WRITE(ILUOUT0,*) 'The variables CPRESOPT = ',CPRESOPT
+ WRITE(ILUOUT0,*) ' NITR = ',NITR
+ IF (CPRESOPT=='RICHA') &
+ WRITE(ILUOUT0,'(A27,F3.1)') ' XRELAX = ',XRELAX
+ WRITE(ILUOUT0,*) ' '
+ WRITE(ILUOUT0,*) 'LOOK correct for the pressure problem in your case.'
+ WRITE(ILUOUT0,*) 'They are stored in the coupling file, and will be'
+ WRITE(ILUOUT0,*) 'the default for the model run.'
+ WRITE(ILUOUT0,*) '***************************************************'
+ WRITE(ILUOUT0,*) ' '
+ EXIT
+ END IF
+!
+!
+!* 6. Next try
+! --------
+!
+END DO
+!
+!
+END SUBROUTINE PRESSURE_IN_PREP
diff --git a/src/ZSOLVER/richardson.f90 b/src/ZSOLVER/richardson.f90
new file mode 100644
index 0000000000000000000000000000000000000000..492454276f2ffcfaaa52e9c4e0b917b6ea70b963
--- /dev/null
+++ b/src/ZSOLVER/richardson.f90
@@ -0,0 +1,243 @@
+!MNH_LIC Copyright 1994-2014 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.
+!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source$ $Revision$
+! MASDEV4_7 solver 2006/05/18 13:07:25
+!-----------------------------------------------------------------
+! ######spl
+ MODULE MODI_RICHARDSON
+! ######################
+!
+INTERFACE
+!
+ SUBROUTINE RICHARDSON(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ, &
+ PRHODJ,PTHETAV,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
+ KIFAXX,KIFAXY,KITR,KTCOUNT,PRELAX,PY,PPHI)
+!
+IMPLICIT NONE
+!
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at
+ ! time t
+!
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+INTEGER, INTENT(IN) :: KTCOUNT ! present iteration of the
+ ! model temporal counter
+!
+REAL, INTENT(IN) :: PRELAX ! relaxation coefficient
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+END SUBROUTINE RICHARDSON
+!
+END INTERFACE
+!
+END MODULE MODI_RICHARDSON
+! ######spl
+ SUBROUTINE RICHARDSON(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ, &
+ PRHODJ,PTHETAV,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
+ KIFAXX,KIFAXY,KITR,KTCOUNT,PRELAX,PY,PPHI)
+! ########################################################################
+!
+!!**** *RICHARDSON * - solve the elliptic pressure equation by the Richardson's
+!! method
+!!
+!! PURPOSE
+!! -------
+! The purpose of this routine is to solve the elliptic pressure equation
+! by the Richardson's method preconditioned by the flat Laplacian.
+!
+!!** METHOD
+!! ------
+!!
+!! The equation to be solved reads:
+!!
+!! Q (PHI) = Y
+!!
+!! where Q is the quasi-Laplacian ( subroutine QLAP ) and PHI the pressure
+!! function.
+!! We precondition the problem by the operator F :
+!! -1 -1
+!! F * Q (PHI) = F (Y)
+!! F represents the flat quasi-laplacian ie. without orography. Its
+!! inversion is realized in the routine FLAT_INV.
+!! This equation is solved with the Richardson's method:
+!!
+!! (n+1) (n) ( -1 -1 (n) )
+!! PHI = PHI + RELAX ( F ( Y ) - F * Q ( PHI ) )
+!! ( )
+!!
+!! EXTERNAL
+!! --------
+!! Subroutine GDIV: compute J times the divergence of 1/J times a vector
+!! Function QLAP: compute the complete quasi-Laplacian Q
+!! Subroutine FLAT_INV : invert the flat quasi-laplacien F
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CONF: model configuration
+!! CCONF: option for the first time step
+!! of the segment ( start or restart)
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine RICHARDSON)
+!!
+!! AUTHOR
+!! ------
+!! P. HÅreil and J. Stein * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 25/07/94
+!! 14/01/97 Durran anelastic equation (Stein,Lafore)
+!! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_CONF
+USE MODI_GDIV
+USE MODI_QLAP
+USE MODI_FLAT_INV
+!
+IMPLICIT NONE
+!
+!* 0.1 declarations of 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
+!
+ ! Metric coefficients:
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX ! d*zx
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY ! d*zy
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at
+ ! time t
+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) :: PRHOM ! 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.
+ ! 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) :: KITR ! number of iterations for the
+ ! pressure solver
+INTEGER, INTENT(IN) :: KTCOUNT ! present iteration of the
+ ! model temporal counter
+!
+REAL, INTENT(IN) :: PRELAX ! relaxation coefficient
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: JM ! loop index
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZF_1_Y
+ ! RHS of the preconditioned problem
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZCORREC
+ ! iterative correction to the solution
+!
+REAL, DIMENSION(SIZE(PPHI,1),SIZE(PPHI,2),SIZE(PPHI,3)) :: ZWORK
+ ! quasi-laplacien Q of the iterative solution
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE THE RHS OF THE PRECONDITIONED PROBLEM
+! ---------------------------------------------
+!
+CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,ZF_1_Y) ! F ( Y )
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. ITERATIVE LOOP
+! --------------
+!
+IF ( KTCOUNT < 1 .OR. ( KTCOUNT == 1 .AND. CCONF == 'START') ) THEN
+ PPHI = ZF_1_Y ! when no first guess is available, we take the solution
+ ! for the flat problem
+END IF
+!
+DO JM = 1,KITR
+!
+ ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI)
+ ! Q (PHI)
+!
+ ZCORREC = 0.
+!
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZCORREC ) ! F * Q (PHI)
+!
+! -1 -1
+! update the iterative solution PHI = PHI + relax* (F (Y) - F * Q (PHI))
+ PPHI = PPHI + PRELAX * (ZF_1_Y - ZCORREC)
+!
+END DO
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE RICHARDSON
diff --git a/src/ZSOLVER/viscosity.f90 b/src/ZSOLVER/viscosity.f90
new file mode 100644
index 0000000000000000000000000000000000000000..d02759b83919be70fb614b6fe82bae858675325e
--- /dev/null
+++ b/src/ZSOLVER/viscosity.f90
@@ -0,0 +1,341 @@
+!MNH_LIC Copyright 1994-2020 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_VISCOSITY
+! #####################
+!
+INTERFACE
+!
+ SUBROUTINE VISCOSITY(HLBCX, HLBCY, KRR, KSV, PNU, PPRANDTL, &
+ OVISC_UVW, OVISC_TH, OVISC_SV, OVISC_R, &
+ ODRAG, &
+ PUT, PVT, PWT, PTHT, PRT, PSVT, &
+ PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRUS, PRVS, PRWS, PRTHS, PRRS, PRSVS,PDRAG )
+!
+ IMPLICIT NONE
+!
+!* 0.1 Declarations of arguments:
+!
+! X and Y lateral boundary conditions
+ CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY
+!
+ INTEGER, INTENT(IN) :: KRR ! number of moist variables
+ INTEGER, INTENT(IN) :: KSV ! number of scalar variables
+!
+ REAL, INTENT(IN) :: PNU ! viscosity coefficient
+ REAL, INTENT(IN) :: PPRANDTL ! Parandtl number
+!
+!
+! logical switches
+ LOGICAL, INTENT(IN) :: OVISC_UVW ! momentum
+ LOGICAL, INTENT(IN) :: OVISC_TH ! theta
+ LOGICAL, INTENT(IN) :: OVISC_SV ! scalar tracer
+ LOGICAL, INTENT(IN) :: OVISC_R ! moisture
+ LOGICAL, INTENT(IN) :: ODRAG ! noslip/freeslip
+!
+!
+! input variables at time t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PVT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PWT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT
+!
+!
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! rho_ref * Jacobian
+!
+ REAL, DIMENSION(:,:), INTENT(IN) :: PDRAG ! Array -1/1 defining where the no-slipcondition is applied
+! metric coefficients
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY
+!
+! output source terms
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS, PRSVS
+!
+ END SUBROUTINE VISCOSITY
+!
+END INTERFACE
+!
+END MODULE MODI_VISCOSITY
+!
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE VISCOSITY(HLBCX, HLBCY, KRR, KSV, PNU, PPRANDTL, &
+ OVISC_UVW, OVISC_TH, OVISC_SV, OVISC_R, &
+ ODRAG, &
+ PUT, PVT, PWT, PTHT, PRT, PSVT, &
+ PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY, &
+ PRUS, PRVS, PRWS, PRTHS, PRRS, PRSVS,PDRAG )
+!
+! IMPLICIT ARGUMENTS
+! ------------------
+! Module MODD_PARAMETERS: JPHEXT, JPVEXT
+! Module MODD_CONF: L2D
+!
+!! AUTHOR
+!! ------
+!! Jeanne Colin * CNRM *
+!!
+!! MODIFICATIONS
+!! -------------
+!! 01/18 (C.Lac) Add budgets
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! P. Wautelet 08/11/2019: corrected wrong budget name VISC_BU_RU -> VISC_BU_RTH
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+ USE MODI_LAP_M
+ USE MODI_SHUMAN
+ USE MODD_PARAMETERS
+ USE MODD_CONF
+ USE MODD_VISCOSITY
+ USE MODD_DRAG_n
+ USE MODD_BUDGET
+ USE MODE_ll
+ USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+ USE MODI_BUDGET
+!
+!-------------------------------------------------------------------------------
+!
+ IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+!
+! X and Y lateral boundary conditions
+ CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY
+!
+ INTEGER, INTENT(IN) :: KRR ! number of moist variables
+ INTEGER, INTENT(IN) :: KSV ! number of scalar variables
+!
+ REAL, INTENT(IN) :: PPRANDTL ! Parandtl number
+ REAL, INTENT(IN) :: PNU ! viscous diffusion rate
+!
+! logical switches
+ LOGICAL, INTENT(IN) :: OVISC_UVW ! momentum
+ LOGICAL, INTENT(IN) :: OVISC_TH ! theta
+ LOGICAL, INTENT(IN) :: OVISC_SV ! scalar tracer
+ LOGICAL, INTENT(IN) :: OVISC_R ! moisture
+ LOGICAL, INTENT(IN) :: ODRAG ! noslip/freeslip
+!
+! input variables at time t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PVT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PWT
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT
+!
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! rho_ref * Jacobian
+!
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PDRAG ! Array -1/1 defining where the no-slip condition is applied
+
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZX
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZY
+!
+! output source terms
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS, PRSVS
+!
+!
+!* 0.2 Declarations of local variables
+!
+ INTEGER :: IK ! counter
+ INTEGER :: IKB, IKE
+!
+ REAL, DIMENSION(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)) :: ZTMP ! temp storage
+ REAL, DIMENSION(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)) :: ZLAPu ! temp storage
+ REAL, DIMENSION(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)) :: ZLAPv ! temp storage
+ REAL, DIMENSION(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)) :: ZY1 ! temp storage
+ REAL, DIMENSION(SIZE(PWT,1),SIZE(PWT,2),SIZE(PWT,3)) :: ZY2 ! temp storage
+!
+!
+INTEGER :: IIU,IJU,IKU ! I,J,K array sizes
+!
+INTEGER :: JI,JJ,JK ! I loop index
+INTEGER :: IINFO_ll
+TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+!
+!
+!-------------------------------------------------------------------------------
+!
+IIU=SIZE(PWT,1)
+IJU=SIZE(PWT,2)
+IKU=SIZE(PWT,3)
+
+!* 1. Viscous forcing for potential temperature
+! -----------------------------------------
+!
+!
+IF (OVISC_TH) THEN
+!
+!
+ PRTHS = PRTHS + PNU/PPRANDTL * &
+ LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHT)
+!
+!
+END IF
+!
+IF (LBUDGET_TH) CALL BUDGET (PRTHS,4,'VISC_BU_RTH')
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. Viscous forcing for moisture
+! ----------------------------
+!
+IF (OVISC_R .AND. (SIZE(PRT,1) > 0)) THEN
+!
+!
+ DO IK = 1, KRR
+ PRRS(:,:,:,IK) = PRRS(:,:,:,IK) + PNU/PPRANDTL * &
+ LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PRT(:,:,:,IK))
+ END DO
+!
+!
+END IF
+!
+IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1),6,'VISC_BU_RRV')
+IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'VISC_BU_RRC')
+IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3),8,'VISC_BU_RRR')
+IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'VISC_BU_RRI')
+IF (LBUDGET_RS) CALL BUDGET (PRRS(:,:,:,5),10,'VISC_BU_RRS')
+IF (LBUDGET_RG) CALL BUDGET (PRRS(:,:,:,6),11,'VISC_BU_RRG')
+IF (LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'VISC_BU_RRH')
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. Viscous forcing for passive scalars
+! -----------------------------------
+!
+IF (OVISC_SV .AND. (SIZE(PSVT,1) > 0)) THEN
+!
+!
+ DO IK = 1, KSV
+ PRSVS(:,:,:,IK) = PRSVS(:,:,:,IK) + PNU/PPRANDTL * &
+ LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PSVT(:,:,:,IK))
+ END DO
+!
+END IF
+!
+IF (LBUDGET_SV) THEN
+ DO IK = 1, KSV
+ CALL BUDGET (PRSVS(:,:,:,IK), 12+IK, 'VISC_BU_RSV')
+ END DO
+END IF
+!
+
+!-------------------------------------------------------------------------------
+!
+!* 4. Viscous forcing for momentum
+! ----------------------------
+!
+IF (OVISC_UVW) THEN
+!
+!* 4.1 U - component
+!
+!
+ ZY1 = MXF(PUT)
+ IF (ODRAG) THEN
+ ZY1(:,:,1) = PDRAG * ZY1(:,:,2)
+ ENDIF
+!
+!
+ ZLAPu = LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX, &
+ PDZY,PDZZ,PRHODJ,ZY1)
+!! Update halo to compute the source term
+ NULLIFY(TZFIELDS_ll)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZLAPu, 'VISCOSITY::ZLAPu' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+ PRUS = PRUS + MXM(PNU*ZLAPu)
+!
+!* 4.2 V - component
+! -------------
+!
+ IF (.NOT. L2D) THEN
+
+ ZY2 = MYF(PVT)
+ IF (ODRAG) THEN
+ ZY2(:,:,1) = PDRAG * ZY2(:,:,2)
+ ENDIF
+!
+ ZLAPv = LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX, &
+ PDZY,PDZZ,PRHODJ,ZY2)
+!! Update halo to compute the source term
+!
+ NULLIFY(TZFIELDS_ll)
+ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZLAPv, 'VISCOSITY::ZLAPv' )
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+!
+ PRVS = PRVS + MYM(PNU*ZLAPv)
+
+ENDIF
+
+!
+!* 4.3 W - component
+! -------------
+!
+ IKB = JPVEXT + 1
+ IKE = SIZE(PWT,3) - JPVEXT
+
+ ZTMP = MZF(PWT)
+!
+ IF (ODRAG) THEN
+ WHERE (PDRAG==-1)
+ ZTMP(:,:,IKB) = 0.
+ ENDWHERE
+ ENDIF
+!
+ DO IK = 1,JPVEXT
+ ZTMP(:,:,IK) = ZTMP(:,:,IKB)
+ ZTMP(:,:,IKE+IK) = ZTMP(:,:,IKE)
+ END DO
+!
+ ZTMP = MZM( PNU * &
+ LAP_M(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTMP) )
+!
+ DO IK = 1,JPVEXT
+ ZTMP(:,:,IK) = ZTMP(:,:,IKB)
+ ZTMP(:,:,IKE+IK) = ZTMP(:,:,IKE)
+ END DO
+ PRWS = PRWS + ZTMP
+!
+!!! Debug provisoire dans le cas ou le noslip est applique jusqu'au bord de
+!sortie de flux en OPEN
+ IF ( LWEST_ll().AND.(ODRAG)) THEN
+ IF ( MINVAL(PDRAG(IIU,:))== -1) THEN
+ DO JK=1,IKU
+ WHERE(PDRAG(IIU,:)== -1)
+ PRUS(IIU,:,JK) = PRUS(IIU-1,:,JK)
+ PRVS(IIU,:,JK) = PRVS(IIU-1,:,JK)
+ PRWS(IIU,:,JK) = PRWS(IIU-1,:,JK)
+ ENDWHERE
+ END DO
+ ENDIF
+ ENDIF
+END IF
+!
+IF (LBUDGET_U) CALL BUDGET (PRUS,1,'VISC_BU_RU')
+IF (LBUDGET_V) CALL BUDGET (PRVS,2,'VISC_BU_RV')
+IF (LBUDGET_W) CALL BUDGET (PRWS,2,'VISC_BU_RW')
+!
+END SUBROUTINE VISCOSITY