diff --git a/src/ZSOLVER/ini_elecn.f90 b/src/ZSOLVER/ini_elecn.f90
deleted file mode 100644
index 2134cd074124f5d446c12e09d4aa547f88ab2236..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/ini_elecn.f90
+++ /dev/null
@@ -1,327 +0,0 @@
-!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(INOUT) :: PDZX    ! metric coefficient dzx
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: 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(INOUT) :: PDZX    ! metric coefficient dzx
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: 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
deleted file mode 100644
index 93725d718d648b744a1ec45aba8a91b4b0b51189..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/ini_field_elec.f90
+++ /dev/null
@@ -1,276 +0,0 @@
-!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(INOUT) ::  PDZX     ! Metric coefficients
-REAL, DIMENSION(:,:,:),  INTENT(INOUT) ::  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(INOUT) ::  PDZX  ! Metric coefficients
-REAL, DIMENSION(:,:,:),  INTENT(INOUT) ::  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
deleted file mode 100644
index cc3fa88163ca93c2245db926181e7baab89193c4..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/lap_m.f90
+++ /dev/null
@@ -1,228 +0,0 @@
-!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(INOUT)  :: PDZX      ! d*zx 
-REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: 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(INOUT)  :: PDZX      ! d*zx 
-REAL, DIMENSION(:,:,:), INTENT(INOUT)  :: 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/mass_leak.f90 b/src/ZSOLVER/mass_leak.f90
deleted file mode 100644
index 13b95970f974a9d73a67058d8310e9a84e614786..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/mass_leak.f90
+++ /dev/null
@@ -1,339 +0,0 @@
-!MNH_LIC Copyright 1995-2021 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_MASS_LEAK
-!####################
-!
-INTERFACE
-!
-          SUBROUTINE MASS_LEAK (PDXX,PDYY,HLBCX,HLBCY,PLINMASS,PRHODJ,PRUS,PRVS)
-!
-REAL, DIMENSION(:,:,:),         INTENT(IN) :: PDXX    ! metric coefficient dxx
-REAL, DIMENSION(:,:,:),         INTENT(IN) :: PDYY    ! metric coefficient dyy
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX   ! type of lateral boundary
-!                                                     ! condition (i=IB, i=IE+1)
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY   ! type of lateral boundary
-!                                                     ! condition (j=JB, j=JE+1)
-REAL,                          INTENT(IN) :: PLINMASS ! lineic mass through open
-!                                                     ! boundaries
-REAL, DIMENSION(:,:,:),        INTENT(IN)  :: PRHODJ  ! rhodref*J
-REAL, DIMENSION(:,:,:),        INTENT(INOUT) :: PRUS  ! Horizontal
-REAL, DIMENSION(:,:,:),        INTENT(INOUT) :: PRVS  ! momentum  tendencies
-!
-END SUBROUTINE MASS_LEAK
-!
-END INTERFACE
-!
-END MODULE MODI_MASS_LEAK
-!
-!     #####################################################################
-      SUBROUTINE MASS_LEAK(PDXX,PDYY,HLBCX,HLBCY,PLINMASS,PRHODJ,PRUS,PRVS)
-!     #####################################################################
-!
-!!***   *MASS_LEAK* - assures global non-divergence condition
-!!
-!!    PURPOSE
-!!    -------
-!!
-!!      This routine changes the horizontal reference dry mass fluxes through
-!!    the open boundary condition faces to set the global divergence in the
-!!    model domain to zero.
-!!
-!!**  METHOD
-!!    ------
-!!
-!!  1) The leak term is computed as:
-!!
-!!          --                     -- IE+1   --                     -- JE+1
-!!          | JE KE                 |        | IE KE                 |
-!!          | _  _                  |        | _  _                  |
-!!          | \  \    1        _ _  |        | \  \    1        _ _  |
-!!   ZLEAK= | /  /   --- PRUS dydz  |   +    | /  /   --- PRVS dxdz  |
-!!          | -  -   dxx            |        | -  -   dyy            |
-!!          | JB KB                 |        | IB KB                 |
-!!          --                     -- i=IB   --                     -- j=JB
-!!
-!!  2) Then the correction wind value ZRUSTOP is found as
-!!               ZLEAK
-!!    ZRUSTOP= ----------
-!!              PLINMASS
-!!
-!!     where PLINMASS is the lineic mass through the open lateral boundaries.
-!!
-!!  3) This horizontal wind correction is applied on the normal wind of
-!!     open lateral boundaries only.
-!!
-!!
-!!    EXTERNAL
-!!    --------
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!      Module MODD_PARAMETERS : contains declaration of parameter variables
-!!
-!!         JPHEXT   : Horizontal external points number
-!!         JPVEXT   : Vertical external points number
-!!
-!!    REFERENCE
-!!    ---------
-!!
-!!      book2
-!!
-!!    AUTHOR
-!!    ------
-!!
-!!      V. Masson   Meteo-France
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original     09/02/95
-!!      Modification 20/10/97  (J.P.Lafore) introduction of 'DAVI' type of lbc
-!!                  15/06/98  (D.Lugato, R.Guivarch) Parallelisation
-!!                   05/06    Suppression of Davies type of lbc
-!!
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-USE MODD_PARAMETERS
-!
-USE MODE_ll
-USE MODE_MPPDB
-!JUAN
-USE MODE_REPRO_SUM
-!JUAN
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK,   ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-!
-IMPLICIT NONE
-!
-!*       0.1   declarations of dummy arguments
-!
-REAL, DIMENSION(:,:,:),        INTENT(IN) :: PDXX    ! metric coefficient dxx
-REAL, DIMENSION(:,:,:),        INTENT(IN) :: PDYY    ! metric coefficient dyy
-CHARACTER(LEN=4), DIMENSION(2),INTENT(IN) :: HLBCX   ! type of lateral boundary
-!                                                    ! condition (i=IB, i=IE+1)
-CHARACTER(LEN=4), DIMENSION(2),INTENT(IN) :: HLBCY   ! type of lateral boundary
-!                                                    ! condition (j=JB, j=JE+1)
-REAL,                          INTENT(IN) :: PLINMASS! lineic mass through open
-!                                                    ! boundaries
-REAL, DIMENSION(:,:,:),        INTENT(IN) :: PRHODJ  ! rhodref*J
-REAL, DIMENSION(:,:,:),        INTENT(INOUT) :: PRUS ! Horizontal
-REAL, DIMENSION(:,:,:),        INTENT(INOUT) :: PRVS ! momentum  tendencies
-!
-!*       0.2   declarations of local variables
-!
-!JUAN16
-REAL                               :: ZLEAK     ! total leak of mass
-REAL, POINTER, CONTIGUOUS, DIMENSION (:,:) :: ZLEAK_W_2D , ZLEAK_E_2D , ZLEAK_S_2D , ZLEAK_N_2D
-!JUAN16
-
-REAL                :: ZUSTOP     ! wind correction!
-INTEGER             :: IIB        ! indice I Beginning in x direction
-INTEGER             :: IJB        ! indice J Beginning in y direction
-INTEGER             :: IKB        ! indice K Beginning in z direction
-INTEGER             :: IIE        ! indice I End       in x direction
-INTEGER             :: IJE        ! indice J End       in y direction
-INTEGER             :: IKE        ! indice K End       in z direction
-INTEGER             :: JI         ! Loop index in x direction
-INTEGER             :: JJ         ! Loop index in y direction
-INTEGER             :: JK         ! Loop index in z direction
-!
-INTEGER             :: IINFO_ll   ! return code of parallel routine
-!
-LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH
-REAL    :: ZLEAK_W,ZLEAK_E,ZLEAK_S,ZLEAK_N
-
-!$acc data present( PDXX, PDYY, PRHODJ, PRUS, PRVS )
-
-IF (MPPDB_INITIALIZED) THEN
-  !Check all IN arrays
-  CALL MPPDB_CHECK(PDXX,  "MASS_LEAK beg:PDXX")
-  CALL MPPDB_CHECK(PDYY,  "MASS_LEAK beg:PDYY")
-  CALL MPPDB_CHECK(PRHODJ,"MASS_LEAK beg:PRHODJ")
-  !Check all INOUT arrays
-  CALL MPPDB_CHECK(PRUS,"MASS_LEAK beg:PRUS")
-  CALL MPPDB_CHECK(PRVS,"MASS_LEAK beg:PRVS")
-END IF
-!-------------------------------------------------------------------------------
-!
-!*       1.    COMPUTE DIMENSIONS OF ARRAYS AND OTHER INDICES:
-!              ----------------------------------------------
-!
-CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
-!
-IKB = 1 + JPVEXT
-IKE = SIZE(PRUS,3) - JPVEXT
-!
-GWEST  = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() )
-GEAST  = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() )
-GSOUTH = ( HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() )
-GNORTH = ( HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() )
-!
-!-------------------------------------------------------------------------------
-!
-!*       2.    COMPUTATION OF LEAK
-!              -------------------
-!
-ZLEAK=0.
-ZLEAK_E=0.
-ZLEAK_W=0.
-ZLEAK_S=0.
-ZLEAK_N=0.
-!
-#ifndef MNH_OPENACC
-IF( HLBCX(1) /= 'CYCL' ) THEN
-   ALLOCATE( ZLEAK_W_2D(IIB:IIB,IJB:IJE))
-   ALLOCATE( ZLEAK_E_2D(IIE+1:IIE+1,IJB:IJE))
-END IF
-IF( HLBCY(1) /= 'CYCL' ) THEN
-   ALLOCATE( ZLEAK_S_2D(IIB:IIE,IJB:IJB))
-   ALLOCATE( ZLEAK_N_2D(IIB:IIE,IJE+1:IJE+1))
-END IF
-#else
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-IF( HLBCX(1) /= 'CYCL' ) THEN
-   CALL MNH_MEM_GET(ZLEAK_W_2D , IIB  ,IIB   , IJB,IJE )
-   CALL MNH_MEM_GET(ZLEAK_E_2D , IIE+1,IIE+1 , IJB,IJE )
-END IF
-IF( HLBCY(1) /= 'CYCL' ) THEN
-   CALL MNH_MEM_GET(ZLEAK_S_2D , IIB,IIE , IJB  ,IJB   )
-   CALL MNH_MEM_GET(ZLEAK_N_2D , IIB,IIE , IJE+1,IJE+1 )
-END IF
-#endif
-!
-IF( HLBCX(1) /= 'CYCL' ) THEN
-   !$acc kernels present(ZLEAK_W_2D) async
-   ZLEAK_W_2D(:,:) = 0.0   
-   IF( GWEST ) THEN
-      !$acc loop seq
-      DO JK=IKB,IKE
-         DO JJ=IJB,IJE
-            ZLEAK_W_2D(IIB,JJ) = ZLEAK_W_2D(IIB,JJ) - 1./PDXX(IIB,JJ,JK) *PRUS(IIB,JJ,JK)
-         END DO
-      END DO
-   END IF
-   !$acc end kernels
-   !
-   !$acc kernels present(ZLEAK_E_2D) async   
-   ZLEAK_E_2D(:,:) = 0.0
-   IF( GEAST ) THEN
-      !$acc loop seq
-      DO JK=IKB,IKE
-         DO JJ=IJB,IJE
-            ZLEAK_E_2D(IIE+1,JJ) = ZLEAK_E_2D(IIE+1,JJ) + 1./PDXX(IIE+1,JJ,JK)*PRUS(IIE+1,JJ,JK)
-         END DO
-      END DO       
-   END IF
-   !$acc end kernels
-   !
-   !$acc wait
-   !
-   !$acc update host(ZLEAK_W_2D,ZLEAK_E_2D)
-   ZLEAK_W = SUM_DD_R2_ll(ZLEAK_W_2D)
-   ZLEAK_E = SUM_DD_R2_ll(ZLEAK_E_2D)
-END IF
-!
-IF( HLBCY(1) /= 'CYCL' ) THEN
-   !
-   !$acc kernels present(ZLEAK_S_2D) async
-   ZLEAK_S_2D(:,:) = 0.0 
-   IF( GSOUTH ) THEN
-      !$acc loop seq
-      DO JK=IKB,IKE
-         DO JI=IIB,IIE  
-            ZLEAK_S_2D(JI,IJB) = ZLEAK_S_2D(JI,IJB) - 1./PDYY(JI,IJB,JK) *PRVS(JI,IJB,JK)
-         END DO
-      END DO
-   END IF
-   !$acc end kernels
-   !
-   !$acc kernels present(ZLEAK_N_2D) async
-   ZLEAK_N_2D(:,:) = 0.0
-   IF ( GNORTH ) THEN
-      !$acc loop seq
-      DO JK=IKB,IKE
-         DO JI=IIB,IIE   
-            ZLEAK_N_2D(JI,IJE+1) = ZLEAK_N_2D(JI,IJE+1)  + 1./PDYY(JI,IJE+1,JK)*PRVS(JI,IJE+1,JK)
-         END DO
-      END DO
-   END IF
-   !$acc end kernels
-   !
-   !$acc wait
-   !
-   !$acc update host(ZLEAK_S_2D,ZLEAK_N_2D)
-   ZLEAK_S = SUM_DD_R2_ll(ZLEAK_S_2D)
-   ZLEAK_N = SUM_DD_R2_ll(ZLEAK_N_2D)   
-!
-END IF
-!
-ZLEAK = ZLEAK_E + ZLEAK_W + ZLEAK_S + ZLEAK_N
-!
-!CALL REDUCESUM_ll(ZLEAK,IINFO_ll)	! we do the reducesum_ll in SUM_DD_R2_ll so we do not do it here
-!
-!-------------------------------------------------------------------------------
-!
-!*       3.    CORRECTION OF WIND ON OPEN BOUNDARIES
-!              -------------------------------------
-!
-ZUSTOP=ZLEAK
-ZUSTOP=ZUSTOP/PLINMASS
-!
-IF (HLBCX(1)=='OPEN' .AND. LWEST_ll() ) THEN
-   !$acc kernels async
-   PRUS(IIB,:,:)=PRUS(IIB,:,:)+ZUSTOP*0.5*(PRHODJ(IIB,:,:)+PRHODJ(IIB-1,:,:))
-   !$acc end kernels
-END IF
-!
-IF (HLBCX(2)=='OPEN' .AND. LEAST_ll() ) THEN
-   !$acc kernels async
-   PRUS(IIE+1,:,:)=PRUS(IIE+1,:,:)-ZUSTOP*0.5*(PRHODJ(IIE+1,:,:)+PRHODJ(IIE,:,:))
-   !$acc end kernels
-END IF
-!
-IF (HLBCY(1)=='OPEN' .AND. LSOUTH_ll() ) THEN
-   !$acc kernels async
-   PRVS(:,IJB,:)=PRVS(:,IJB,:)+ZUSTOP*0.5*(PRHODJ(:,IJB,:)+PRHODJ(:,IJB-1,:))
-   !$acc end kernels
-END IF
-!
-IF (HLBCY(2)=='OPEN' .AND. LNORTH_ll() ) THEN
-   !$acc kernels async
-   PRVS(:,IJE+1,:)=PRVS(:,IJE+1,:)-ZUSTOP*0.5*(PRHODJ(:,IJE+1,:)+PRHODJ(:,IJE,:))
-   !$acc end kernels
-END IF
-!
-!$acc wait
-!
-!
-IF (MPPDB_INITIALIZED) THEN
-  !Check all INOUT arrays
-  CALL MPPDB_CHECK(PRUS,"MASS_LEAK end:PRUS")
-  CALL MPPDB_CHECK(PRVS,"MASS_LEAK end:PRVS")
-END IF
-
-#ifndef MNH_OPENACC
-IF( HLBCX(1) /= 'CYCL' ) THEN
-   DEALLOCATE(ZLEAK_W_2D,ZLEAK_E_2D)
-END IF
-IF( HLBCY(1) /= 'CYCL' ) THEN
-   DEALLOCATE( ZLEAK_S_2D,ZLEAK_N_2D)
-END IF
-#else
-!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
-CALL MNH_MEM_RELEASE()
-#endif
-
-!$acc end data
-
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE MASS_LEAK
diff --git a/src/ZSOLVER/pressure_in_prep.f90 b/src/ZSOLVER/pressure_in_prep.f90
deleted file mode 100644
index 19687655725eb8e16437ab3455c886d92e98fe18..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/pressure_in_prep.f90
+++ /dev/null
@@ -1,308 +0,0 @@
-!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(INOUT) :: PDZX     ! metric coefficient dzx 
-REAL,DIMENSION(:,:,:), INTENT(INOUT) :: 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(INOUT) :: PDZX     ! metric coefficient dzx 
-REAL,DIMENSION(:,:,:), INTENT(INOUT) :: 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/viscosity.f90 b/src/ZSOLVER/viscosity.f90
deleted file mode 100644
index e202cf20e7c246f0a794bf94cfb0771392eb358a..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/viscosity.f90
+++ /dev/null
@@ -1,364 +0,0 @@
-!MNH_LIC Copyright 1994-2021 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(INOUT) :: PDZX
-     REAL, DIMENSION(:,:,:), INTENT(INOUT) :: 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
-!  P. Wautelet    02/2020: use the new data structures and subroutines for budgets
-!  T. Nagel       02/2021: add adhesion condition in case of an IBM-obstacle at the domain top boundary
-!-------------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-!
-  USE MODD_ARGSLIST_ll, ONLY: LIST_ll
-  use modd_budget,      only: lbudget_u,  lbudget_v,  lbudget_w,  lbudget_th, lbudget_rv,  lbudget_rc,  &
-                              lbudget_rr, lbudget_ri, lbudget_rs, lbudget_rg, lbudget_rh,  lbudget_sv,  &
-                              NBUDGET_U,  NBUDGET_V,  NBUDGET_W,  NBUDGET_TH, NBUDGET_RV,  NBUDGET_RC,  &
-                              NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH,  NBUDGET_SV1, &
-                              tbudgets
-  USE MODD_CONF
-  USE MODD_DRAG_n
-  USE MODD_PARAMETERS
-  USE MODD_VISCOSITY
-
-  use mode_budget,      only: Budget_store_init, Budget_store_end
-  USE MODE_ll
-
-  USE MODI_SHUMAN
-  USE MODI_LAP_M
-!
-!-------------------------------------------------------------------------------
-!
-  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(INOUT) :: PDZX
-     REAL, DIMENSION(:,:,:), INTENT(INOUT) :: 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)
-
-if ( lbudget_u  .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_U ), 'VISC', prus (:, :, :)    )
-if ( lbudget_v  .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_V ), 'VISC', prvs (:, :, :)    )
-if ( lbudget_w  .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_W ), 'VISC', prws (:, :, :)    )
-if ( lbudget_th .and. ovisc_th  ) call Budget_store_init( tbudgets(NBUDGET_TH), 'VISC', prths(:, :, :)    )
-if ( lbudget_rv .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RV), 'VISC', prrs (:, :, :, 1) )
-if ( lbudget_rc .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RC), 'VISC', prrs (:, :, :, 2) )
-if ( lbudget_rr .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RR), 'VISC', prrs (:, :, :, 3) )
-if ( lbudget_ri .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RI), 'VISC', prrs (:, :, :, 4) )
-if ( lbudget_rs .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RS), 'VISC', prrs (:, :, :, 5) )
-if ( lbudget_rg .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RG), 'VISC', prrs (:, :, :, 6) )
-if ( lbudget_rh .and. ovisc_r   ) call Budget_store_init( tbudgets(NBUDGET_RH), 'VISC', prrs (:, :, :, 7) )
-if ( lbudget_sv .and. ovisc_sv  ) then
-  do ik = 1, ksv
-    call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ik), 'VISC', prsvs(:, :, :, ik) )
-  end do
-end if
-
-!*       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
-!
-!-------------------------------------------------------------------------------
-!
-!*       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
-!
-!-------------------------------------------------------------------------------
-!
-!*       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
-!
-!-------------------------------------------------------------------------------
-!
-!*       4.    Viscous forcing for momentum
-!	       ----------------------------
-!
-IF (OVISC_UVW) THEN
-!
-!*       4.1   U - component
-!
-!
-      ZY1 = MXF(PUT)
-      IF (ODRAG) THEN
-         ZY1(:,:,1) = PDRAG * ZY1(:,:,2)
-!!Add adhesion condition in case of an IBM-obstacle at the domain top boundary
-!         ZY1(:,:,IKU) = PDRAG * ZY1(:,:,IKE)
-      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  .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_U ), 'VISC', prus (:, :, :)    )
-if ( lbudget_v  .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_V ), 'VISC', prvs (:, :, :)    )
-if ( lbudget_w  .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_W ), 'VISC', prws (:, :, :)    )
-if ( lbudget_th .and. ovisc_th  ) call Budget_store_end( tbudgets(NBUDGET_TH), 'VISC', prths(:, :, :)    )
-if ( lbudget_rv .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RV), 'VISC', prrs (:, :, :, 1) )
-if ( lbudget_rc .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RC), 'VISC', prrs (:, :, :, 2) )
-if ( lbudget_rr .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RR), 'VISC', prrs (:, :, :, 3) )
-if ( lbudget_ri .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RI), 'VISC', prrs (:, :, :, 4) )
-if ( lbudget_rs .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RS), 'VISC', prrs (:, :, :, 5) )
-if ( lbudget_rg .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RG), 'VISC', prrs (:, :, :, 6) )
-if ( lbudget_rh .and. ovisc_r   ) call Budget_store_end( tbudgets(NBUDGET_RH), 'VISC', prrs (:, :, :, 7) )
-if ( lbudget_sv .and. ovisc_sv  ) then
-  do ik = 1, ksv
-    call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ik), 'VISC', prsvs(:, :, :, ik) )
-  end do
-end if
-
-END SUBROUTINE VISCOSITY