diff --git a/src/MNH/anel_balancen.f90 b/src/MNH/anel_balancen.f90
index a0ed07162cfa75de94cc3ad5755d6ac8fe53a9df..d9388e7dce829900bc880cc61fce1c9e86ff292d 100644
--- a/src/MNH/anel_balancen.f90
+++ b/src/MNH/anel_balancen.f90
@@ -174,12 +174,6 @@ INTEGER                       ::  IMI          ! model index
 INTEGER                               ::  IIU_B,IJU_B,IKU
 INTEGER                               ::  IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll
 REAL, DIMENSION(:,:,:), ALLOCATABLE   ::  ZBFB,ZBF_SXP2_YP1_Z
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), ALLOCATABLE   :: ZAF_ZS,ZBF_ZS,ZCF_ZS
-REAL, DIMENSION(:,:)  , ALLOCATABLE   :: ZDXATH_ZS,ZDYATH_ZS
-REAL, DIMENSION(:,:,:), ALLOCATABLE   :: ZRHO_ZS
-REAL, DIMENSION(:),     ALLOCATABLE   :: ZA_K,ZB_K,ZC_K,ZD_K
-#endif
 !
 INTEGER :: IINFO_ll
 TYPE(LIST_ll), POINTER :: TZFIELDS_ll=>NULL()   ! list of fields to exchange
@@ -202,20 +196,6 @@ IKU=SIZE(XRHODJ,3)
 CALL GET_DIM_EXT_ll('B',IIU_B,IJU_B)
 ALLOCATE(ZBFB(IIU_B,IJU_B,IKU))
 !
-#ifdef MNH_MGSOLVER
-IF ( CPRESOPT == 'ZSOLV' ) THEN
-  ALLOCATE(ZAF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(ZBF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(ZCF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(ZDXATH_ZS(IIU_B,IJU_B))
-  ALLOCATE(ZDYATH_ZS(IIU_B,IJU_B))
-  ALLOCATE(ZRHO_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(ZA_K(IKU))
-  ALLOCATE(ZB_K(IKU))
-  ALLOCATE(ZC_K(IKU))
-  ALLOCATE(ZD_K(IKU))
-END IF
-#endif
 !
 CALL GET_DIM_EXTZ_ll('SXP2_YP1_Z',IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll)
 ALLOCATE(ZBF_SXP2_YP1_Z(IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll))
@@ -232,14 +212,7 @@ CALL MPPDB_CHECK3D(XUT,"anel_balancen1-::XUT",PRECISION)
 CALL TRIDZ(CLBCX,CLBCY,XMAP,XDXHAT,XDYHAT,CPRESOPT,                     &
           ZDXHATM,ZDYHATM,ZRHOM,                                        &
           ZAF,ZCF,ZTRIGSX,ZTRIGSY,IIFAXX,IIFAXY,XRHODJ,XTHVREF,XZZ,ZBFY,&
-#ifndef MNH_MGSOLVER
           ZBFB,ZBF_SXP2_YP1_Z) 
-#else
-          ZBFB,ZBF_SXP2_YP1_Z,                                          &
-          ZAF_ZS,ZBF_ZS,ZCF_ZS,                                         &
-          ZDXATH_ZS,ZDYATH_ZS,ZRHO_ZS,                      &
-          ZA_K,ZB_K,ZC_K,ZD_K) !JUAN  FULL ZSOLVER
-#endif
 CALL MPPDB_CHECK3D(XRHODJ,"anel_balancen1-after TRIDZ::XRHODJ",PRECISION)
 !
 !-------------------------------------------------------------------------------
diff --git a/src/MNH/ini_dynamics.f90 b/src/MNH/ini_dynamics.f90
index eed1f93895eb28c5fbcb6cdcb9b0ed3eb1849462..788b346f454621bdc441029e56b102613b613e04 100644
--- a/src/MNH/ini_dynamics.f90
+++ b/src/MNH/ini_dynamics.f90
@@ -27,14 +27,7 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ,                   &
                OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX,                     &
                PDK2U,PDK4U,PDK2TH,PDK4TH,PDK2SV,PDK4SV,OZDIFFU,PZDIFFU_HALO2,& 
                PBFB,                                                         &
-#ifndef MNH_MGSOLVER
                PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
-#else
-               PBF_SXP2_YP1_Z,                                               &
-               PAF_ZS,PBF_ZS,PCF_ZS,                                         &
-               PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                                  &
-               A_K,B_K,C_K,D_K) !JUAN FULL ZSOLVER
-#endif
 !  intent in arguments
 !
 USE MODE_TYPE_ZDIFFU
@@ -180,12 +173,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements of the tri-diag matrix
                                             ! on an b-slice of global physical domain
 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide 
                                                    ! matrix in the pressure eq.
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(OUT) :: A_K,B_K,C_K,D_K
-#endif
 
 END SUBROUTINE INI_DYNAMICS
 !
@@ -213,14 +200,7 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ,                   &
                OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX,                     &
                PDK2U,PDK4U,PDK2TH,PDK4TH,PDK2SV,PDK4SV,OZDIFFU,PZDIFFU_HALO2,& 
                PBFB,                                                         &
-#ifndef MNH_MGSOLVER
                PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
-#else
-               PBF_SXP2_YP1_Z,                                               &
-               PAF_ZS,PBF_ZS,PCF_ZS,                                         &
-               PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                                  &
-               A_K,B_K,C_K,D_K) !JUAN FULL ZSOLVER
-#endif
 !     ######################################################################
 !
 !!****  *INI_DYNAMICS* - routine to initialize the parameters for the dynamics
@@ -470,12 +450,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements of the tri-diag matrix
                                             ! on an b-slice of global physical domain
 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide 
                                                    ! matrix in the pressure eq.
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(OUT) :: A_K,B_K,C_K,D_K
-#endif
 !
 !*       0.2   declarations of local variables
 !
@@ -537,14 +511,7 @@ IF (.NOT.L1D) THEN
             PDXHATM,PDYHATM,PRHOM,PAF,                                    &
             PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,                            &
             PRHODJ,PTHVREF,PZZ,PBFY,PBFB,                                 &
-#ifndef MNH_MGSOLVER
             PBF_SXP2_YP1_Z)
-#else
-            PBF_SXP2_YP1_Z,                                               &
-            PAF_ZS,PBF_ZS,PCF_ZS,                                         &
-            PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                                  &
-            A_K,B_K,C_K,D_K) !JUAN  FULL ZSOLVER
-#endif
 END IF
 !
 !
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 956abcf442326b3bd84f202c1cb539661cc3c977..746403eaf1d53aaa543ee141a006fdb240b9fa94 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1093,31 +1093,6 @@ IF ( CPRESOPT == 'ZRESI' .OR. CPRESOPT == 'ZSOLV' ) THEN
   ALLOCATE(XBFB(IIU_B,IJU_B,IKU))
 END IF
 
-#ifdef MNH_MGSOLVER
-IF ( CPRESOPT == 'ZSOLV' ) THEN
-  ALLOCATE(XAF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(XBF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(XCF_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(XDXATH_ZS(IIU_B,IJU_B))
-  ALLOCATE(XDYATH_ZS(IIU_B,IJU_B))
-  ALLOCATE(XRHO_ZS(IIU_B,IJU_B,IKU))
-  ALLOCATE(XA_K(IKU))
-  ALLOCATE(XB_K(IKU))
-  ALLOCATE(XC_K(IKU))
-  ALLOCATE(XD_K(IKU))
-ELSE
-  ALLOCATE(XAF_ZS(0,0,0))
-  ALLOCATE(XBF_ZS(0,0,0))
-  ALLOCATE(XCF_ZS(0,0,0))
-  ALLOCATE(XDXATH_ZS(0,0))
-  ALLOCATE(XDYATH_ZS(0,0))
-  ALLOCATE(XRHO_ZS(0,0,0))
-  ALLOCATE(XA_K(0))
-  ALLOCATE(XB_K(0))
-  ALLOCATE(XC_K(0))
-  ALLOCATE(XD_K(0))
-END IF
-#endif
 IF ( CPRESOPT == 'ZRESI' ) THEN
   CALL GET_DIM_EXTZ_ll('SXP2_YP1_Z',IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll)
   ALLOCATE(XBF_SXP2_YP1_Z(IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll))
@@ -2274,14 +2249,7 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT,            &
              LMASK_RELAX,XKURELAX,XKVRELAX,XKWRELAX,                          &
              XDK2U,XDK4U,XDK2TH,XDK4TH,XDK2SV,XDK4SV,                         &
              LZDIFFU,XZDIFFU_HALO2,                                           &
-#ifndef MNH_MGSOLVER
              XBFB,XBF_SXP2_YP1_Z                                              )
-#else
-             XBFB,XBF_SXP2_YP1_Z,                                             &
-             XAF_ZS,XBF_ZS,XCF_ZS,                                            &
-             XDXATH_ZS,XDYATH_ZS,XRHO_ZS,                                     &
-             XA_K,XB_K,XC_K,XD_K                                              )
-#endif
 
 !$acc update device( XRHOM, XAF, XBFY, XCF, XTRIGSX, XTRIGSY, NIFAXX, NIFAXY, XBFB, XBF_SXP2_YP1_Z )
 
diff --git a/src/MNH/ini_spectren.f90 b/src/MNH/ini_spectren.f90
index 2428e79798f69c778e2430f2e43152536dcdf5f6..924ce7aff37cf43daf974f984a88b83854fe0e9e 100644
--- a/src/MNH/ini_spectren.f90
+++ b/src/MNH/ini_spectren.f90
@@ -913,14 +913,7 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT,            &
              LMASK_RELAX,XKURELAX,XKVRELAX,XKWRELAX,                          &
              XDK2U,XDK4U,XDK2TH,XDK4TH,XDK2SV,XDK4SV,                         &
              LZDIFFU,XZDIFFU_HALO2,                                           &
-#ifndef MNH_MGSOLVER
              XBFB,XBF_SXP2_YP1_Z                                              )
-#else
-             XBFB,XBF_SXP2_YP1_Z,                                             &
-             XAF_ZS,XBF_ZS,XCF_ZS,                                            &
-             XDXATH_ZS,XDYATH_ZS,XRHO_ZS,                                     &
-             XA_K,XB_K,XC_K,XD_K) 
-#endif
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/modd_dynn.f90 b/src/MNH/modd_dynn.f90
index 71f373d8b0734494c6a94869258f782682e3b182..37772d57529f2ecd4bdae097c22ee0355a6b2a12 100644
--- a/src/MNH/modd_dynn.f90
+++ b/src/MNH/modd_dynn.f90
@@ -66,12 +66,6 @@ TYPE DYN_t
   REAL, DIMENSION(:,:,:), POINTER :: XBFB=>NULL()     ! Vectors giving the non
   REAL, DIMENSION(:,:,:), POINTER :: XBF_SXP2_YP1_Z=>NULL()     ! Vectors giving the non
   REAL, DIMENSION(:),     POINTER :: XAF=>NULL(),XCF=>NULL() ! vanishing elements of the tri-diag matrix in the pressure equation
-#ifdef MNH_MGSOLVER
-  REAL, DIMENSION(:,:,:), POINTER :: XAF_ZS=>NULL(), XBF_ZS=>NULL(), XCF_ZS=>NULL() ! coef for Zsolver
-  REAL, DIMENSION(:,:)  , POINTER :: XDXATH_ZS=>NULL(), XDYATH_ZS=>NULL()
-  REAL, DIMENSION(:,:,:), POINTER :: XRHO_ZS=>NULL()
-  REAL, DIMENSION(:),     POINTER :: XA_K=>NULL(), XB_K=>NULL(), XC_K=>NULL(), XD_K=>NULL()
-#endif
 !
                                                    ! Arrays of sinus or cosinus
                                                    ! values for the FFT 
@@ -192,12 +186,6 @@ REAL, DIMENSION(:,:,:), POINTER :: XBFY=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XBFB=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XBF_SXP2_YP1_Z=>NULL()
 REAL, DIMENSION(:),     POINTER :: XAF=>NULL(),XCF=>NULL()
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), POINTER :: XAF_ZS=>NULL(), XBF_ZS=>NULL(), XCF_ZS=>NULL()
-REAL, DIMENSION(:,:)  , POINTER :: XDXATH_ZS=>NULL(), XDYATH_ZS=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XRHO_ZS=>NULL()
-REAL, DIMENSION(:),     POINTER :: XA_K=>NULL(), XB_K=>NULL(), XC_K=>NULL(), XD_K=>NULL()
-#endif
 REAL, DIMENSION(:),   POINTER :: XTRIGSX=>NULL()
 REAL, DIMENSION(:),   POINTER :: XTRIGSY=>NULL()
 INTEGER, DIMENSION(:), POINTER :: NIFAXX=>NULL()
@@ -281,20 +269,6 @@ DYN_MODEL(KFROM)%XBF_SXP2_YP1_Z=>XBF_SXP2_YP1_Z
 DYN_MODEL(KFROM)%XAF=>XAF
 DYN_MODEL(KFROM)%XCF=>XCF
 
-#ifdef MNH_MGSOLVER
-DYN_MODEL(KFROM)%XAF_ZS=>XAF_ZS
-DYN_MODEL(KFROM)%XBF_ZS=>XBF_ZS
-DYN_MODEL(KFROM)%XCF_ZS=>XCF_ZS
-
-DYN_MODEL(KFROM)%XDXATH_ZS=>XDXATH_ZS
-DYN_MODEL(KFROM)%XDYATH_ZS=>XDYATH_ZS
-DYN_MODEL(KFROM)%XRHO_ZS=>XRHO_ZS
-DYN_MODEL(KFROM)%XA_K=>XA_K
-DYN_MODEL(KFROM)%XB_K=>XB_K
-DYN_MODEL(KFROM)%XC_K=>XC_K
-DYN_MODEL(KFROM)%XD_K=>XD_K
-#endif
-
 DYN_MODEL(KFROM)%XTRIGSX=>XTRIGSX
 DYN_MODEL(KFROM)%XTRIGSY=>XTRIGSY
 DYN_MODEL(KFROM)%XRHOM=>XRHOM
@@ -316,20 +290,6 @@ XBF_SXP2_YP1_Z=>DYN_MODEL(KTO)%XBF_SXP2_YP1_Z
 XAF=>DYN_MODEL(KTO)%XAF
 XCF=>DYN_MODEL(KTO)%XCF
 
-#ifdef MNH_MGSOLVER
-XAF_ZS=>DYN_MODEL(KTO)%XAF_ZS
-XBF_ZS=>DYN_MODEL(KTO)%XBF_ZS
-XCF_ZS=>DYN_MODEL(KTO)%XCF_ZS
-
-XDXATH_ZS=>DYN_MODEL(KFROM)%XDXATH_ZS
-XDYATH_ZS=>DYN_MODEL(KFROM)%XDYATH_ZS
-XRHO_ZS=>DYN_MODEL(KFROM)%XRHO_ZS
-XA_K=>DYN_MODEL(KFROM)%XA_K
-XB_K=>DYN_MODEL(KFROM)%XB_K
-XC_K=>DYN_MODEL(KFROM)%XC_K
-XD_K=>DYN_MODEL(KFROM)%XD_K
-#endif
-
 XTRIGSX=>DYN_MODEL(KTO)%XTRIGSX
 XTRIGSY=>DYN_MODEL(KTO)%XTRIGSY
 NIFAXX=>DYN_MODEL(KTO)%NIFAXX
diff --git a/src/MNH/tridz.f90 b/src/MNH/tridz.f90
index 93bf853351331290e8e11759d82db5d21abb16ba..badb66af95bb4fcec5d7f9f3635a3be06d5d5805 100644
--- a/src/MNH/tridz.f90
+++ b/src/MNH/tridz.f90
@@ -14,14 +14,7 @@ INTERFACE
                       PDXHATM,PDYHATM,PRHOM,                            &
                       PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,            &
                       PRHODJ,PTHVREF,PZZ,PBFY,PBFB,                     &
-#ifndef MNH_MGSOLVER
                       PBF_SXP2_YP1_Z)!JUAN Z_SPLITING
-#else
-                      PBF_SXP2_YP1_Z,                                   &
-                      PAF_ZS,PBF_ZS,PCF_ZS,                             &
-                      PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                      &
-                      A_K,B_K,C_K,D_K) !JUAN  FULL ZSOLVER
-#endif
 !
 IMPLICIT NONE
 !
@@ -57,12 +50,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB     ! elements (bsplit slide) of the
                                                 ! matrix in the pressure eq. 
 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide 
                                                    ! matrix in the pressure eq
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(OUT) :: A_K,B_K,C_K,D_K
-#endif
 !JUAN
 !
                                                  ! arrays of sin or cos values
@@ -88,14 +75,7 @@ END MODULE MODI_TRIDZ
                       PDXHATM,PDYHATM,PRHOM,                            &
                       PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,            &
                       PRHODJ,PTHVREF,PZZ,PBFY,PBFB,                     &
-#ifndef MNH_MGSOLVER
                       PBF_SXP2_YP1_Z)!JUAN Z_SPLITING
-#else
-                      PBF_SXP2_YP1_Z,                                   &
-                      PAF_ZS,PBF_ZS,PCF_ZS,                             &
-                      PDXATH_ZS,PDYATH_ZS,PRHO_ZS,                      &
-                      A_K,B_K,C_K,D_K) !JUAN  FULL ZSOLVER
-#endif
 !    ####################################################################
 !
 !!****  *TRIDZ * - Compute coefficients for the flat operator
@@ -256,12 +236,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB     ! elements (bsplit slide) of the
                                                 ! matrix in the pressure eq. 
 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide 
                                                    ! matrix in the pressure eq.
-#ifdef MNH_MGSOLVER
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
-REAL, DIMENSION(:,:)  , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
-REAL, DIMENSION(:)    , INTENT(OUT) :: A_K,B_K,C_K,D_K
-#endif
 !JUAN
 !                                 
                                                  ! arrays of sin or cos values
@@ -326,7 +300,8 @@ INTEGER :: IXMODE_SXP2_YP1_Z_ll,IYMODE_SXP2_YP1_Z_ll ! number of modes in the x
 !TYPE(DOUBLE_DOUBLE)  , DIMENSION (SIZE(PZZ,3)) :: ZRHOM_ll   , ZDZM_ll
 REAL, ALLOCATABLE, DIMENSION(:,:)              :: ZRHOM_2D   , ZDZM_2D
 #ifdef MNH_MGSOLVER
-REAL, ALLOCATABLE, DIMENSION(:,:,:)            :: ZDZM_ZS
+REAL, DIMENSION(:),     ALLOCATABLE :: ZA_K, ZB_K, ZC_K, ZD_K
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDZM_ZS
 #endif
 !JUAN16
 !
@@ -369,6 +344,15 @@ IIE_ll = IIMAX_ll  + JPHEXT
 IJB_ll = 1 + JPHEXT
 IJE_ll = IJMAX_ll + JPHEXT
 !
+#ifdef MNH_MGSOLVER
+IF ( HPRESOPT == 'ZSOLV' ) THEN
+  ALLOCATE( ZA_K(IKU) )
+  ALLOCATE( ZB_K(IKU) )
+  ALLOCATE( ZC_K(IKU) )
+  ALLOCATE( ZD_K(IKU) )
+END IF
+#endif
+!
 !     the use of local array ZEIGENX and ZEIGEN would require some technical modifications
 !
 ALLOCATE (ZEIGENX_ll(IIMAX_ll+2*JPHEXT))
@@ -399,9 +383,6 @@ ZGWNY = XPI/REAL(IJMAX_ll)
 ZINVMEAN = 1./REAL(IIMAX_ll*IJMAX_ll)
 !JUAN16
 ALLOCATE(ZRHOM_2D(IIB:IIE, IJB:IJE))
-#ifdef MNH_MGSOLVER
-IF ( HPRESOPT == 'ZSOLV' ) PRHO_ZS(:,:,:) = 1.0
-#endif
 !
 DO JK = 1,SIZE(PZZ,3)
    IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN
@@ -417,9 +398,6 @@ DO JK = 1,SIZE(PZZ,3)
          END DO
       END DO
    END IF
-#ifdef MNH_MGSOLVER
-   IF ( HPRESOPT == 'ZSOLV' ) PRHO_ZS(IIB:IIE,IJB:IJE,JK) = ZRHOM_2D(:,:)
-#endif
    !  global sum
    PRHOM(JK) =  SUM_DD_R2_ll(ZRHOM_2D) * ZINVMEAN
 END DO
@@ -480,13 +458,6 @@ PDXHATM =0.
 !     . local sum
 IF (LCARTESIAN) THEN
   PDXHATM = SUM_1DFIELD_ll ( PDXHAT,'X',IIB_ll,IIE_ll,IINFO_ll)
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    DO JJ=1,SIZE(PDXATH_ZS,2)
-      PDXATH_ZS(:,JJ) = PDXHAT(:)
-    END DO
-  END IF
-#endif
 ELSE
   ! Extraction of x-slice PMAP at j=(IJB_ll+IJE_ll)/2
   CALL GET_SLICE_ll (PMAP,'X',(IJB_ll+IJE_ll)/2,ZXMAP(IIB:IIE) &
@@ -494,13 +465,6 @@ ELSE
   ! initialize the work array = PDXHAT/ZXMAP
   ZWORKX(IIB:IIE) = PDXHAT(IIB:IIE)/ ZXMAP (IIB:IIE)
   PDXHATM = SUM_1DFIELD_ll ( ZWORKX,'X',IIB_ll,IIE_ll,IINFO_ll)
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    DO JJ=1,SIZE(PDXATH_ZS,2)
-      PDXATH_ZS(:,JJ) = PDXHAT(:) / PMAP(:,JJ)
-    END DO
-  END IF
-#endif
 END IF
 !    . division to complete sum
 PDXHATM = PDXHATM /  REAL(IIMAX_ll)
@@ -513,13 +477,6 @@ PDXHATM = PDXHATM /  REAL(IIMAX_ll)
 PDYHATM = 0.
 IF (LCARTESIAN) THEN
   PDYHATM = SUM_1DFIELD_ll ( PDYHAT,'Y',IJB_ll,IJE_ll,IINFO_ll)
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    DO JI=1,SIZE(PDYATH_ZS,1)
-      PDYATH_ZS(JI,:) = PDYHAT(:)
-    END DO
-  END IF
-#endif
 ELSE
   ! Extraction of y-slice PMAP at i=IIB_ll+IIE_ll/2 
   CALL GET_SLICE_ll (PMAP,'Y',(IIB_ll+IIE_ll)/2,ZYMAP(IJB:IJE) &
@@ -527,13 +484,6 @@ ELSE
   ! initialize the work array = PDYHAT / ZYMAP
   ZWORKY(IJB:IJE) = PDYHAT(IJB:IJE) / ZYMAP (IJB:IJE)
   PDYHATM = SUM_1DFIELD_ll ( ZWORKY,'Y',IJB_ll,IJE_ll,IINFO_ll)
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    DO JI=1,SIZE(PDYATH_ZS,1)
-      PDYATH_ZS(JI,:) = PDYHAT(:) / PMAP(JI,:)
-    END DO
-  END IF
-#endif
 END IF
 !    . division to complete sum
 PDYHATM= PDYHATM / REAL(IJMAX_ll)
@@ -545,11 +495,9 @@ PDYHATM= PDYHATM / REAL(IJMAX_ll)
 !
 #ifdef MNH_MGSOLVER
 IF ( HPRESOPT == 'ZSOLV' ) THEN
-  PAF_ZS = 1.0
-  PCF_ZS = 1.0
-  A_K = 0.0
-  B_K = 0.0
-  C_K = 0.0
+  ZA_K = 0.0
+  ZB_K = 0.0
+  ZC_K = 0.0
 END IF
 #endif
 IF ( HPRESOPT /= 'ZSOLV' ) THEN
@@ -562,14 +510,9 @@ END IF
 #ifdef MNH_MGSOLVER
 IF ( HPRESOPT == 'ZSOLV' ) THEN
   DO JK = IKB,IKE
-    PAF_ZS(IIB:IIE,IJB:IJE,JK) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,JK-1) + PRHO_ZS(IIB:IIE,IJB:IJE,JK)   ) &
-                                       / ZDZM_ZS(IIB:IIE,IJB:IJE,JK)   **2 
-    PCF_ZS(IIB:IIE,IJB:IJE,JK) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,JK) + PRHO_ZS(IIB:IIE,IJB:IJE,JK+1)   ) &
-                                       / ZDZM_ZS(IIB:IIE,IJB:IJE,JK+1) **2 
-
-    D_K(JK) = PRHOM(JK) ! / ZDZM(JK)
-    B_K(JK) = ( 0.5 * ( PRHOM(JK)   + PRHOM(JK+1) ) / ZDZM(JK+1) **2 ) / D_K(JK)
-    C_K(JK) = ( 0.5 * ( PRHOM(JK-1) + PRHOM(JK)   ) / ZDZM(JK)   **2 ) / D_K(JK)
+    ZD_K(JK) = PRHOM(JK) ! / ZDZM(JK)
+    ZB_K(JK) = ( 0.5 * ( PRHOM(JK)   + PRHOM(JK+1) ) / ZDZM(JK+1) **2 ) / ZD_K(JK)
+    ZC_K(JK) = ( 0.5 * ( PRHOM(JK-1) + PRHOM(JK)   ) / ZDZM(JK)   **2 ) / ZD_K(JK)
   END DO
 END IF
 #endif
@@ -582,8 +525,8 @@ IF ( HPRESOPT /= 'ZSOLV' ) THEN
 END IF
 #ifdef MNH_MGSOLVER
 IF ( HPRESOPT == 'ZSOLV' ) THEN
-  D_K(IKE+1) = PRHOM(IKE+1) ! / ZDZM(IKE+1)
-  C_K(IKE+1) = ( -0.5 * ( PRHOM(IKE)   + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2 ) / D_K(IKE+1)
+  ZD_K(IKE+1) = PRHOM(IKE+1) ! / ZDZM(IKE+1)
+  ZC_K(IKE+1) = ( -0.5 * ( PRHOM(IKE)   + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2 ) / ZD_K(IKE+1)
 END IF
 #endif
 
@@ -592,37 +535,29 @@ IF ( HPRESOPT /= 'ZSOLV' ) THEN
 END IF
 #ifdef MNH_MGSOLVER
 IF ( HPRESOPT == 'ZSOLV' ) THEN
-  D_K(IKB-1) = PRHOM(IKB-1) ! / ZDZM(IKB-1)
-  B_K(IKB-1) = ( 0.5 * ( PRHOM(IKB-1) + PRHOM(IKB)   ) / ZDZM(IKB)**2 ) / D_K(IKB-1)
+  ZD_K(IKB-1) = PRHOM(IKB-1) ! / ZDZM(IKB-1)
+  ZB_K(IKB-1) = ( 0.5 * ( PRHOM(IKB-1) + PRHOM(IKB)   ) / ZDZM(IKB)**2 ) / ZD_K(IKB-1)
 END IF
 #endif
 !
 IF ( HPRESOPT /= 'ZSOLV' ) PAF(IKB-1) =  0.0 ! 0.5 * ( PRHOM(IKB-1) + PRHOM(IKB)   ) / ZDZM(IKB)   **2
 #ifdef MNH_MGSOLVER
-IF ( HPRESOPT == 'ZSOLV' ) C_K(IKB-1) =  0.0
+IF ( HPRESOPT == 'ZSOLV' ) ZC_K(IKB-1) =  0.0
 
 #endif
 IF ( HPRESOPT /= 'ZSOLV' ) PCF(IKE+1) =  0.0 ! 0.5 * ( PRHOM(IKE)   + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2
 #ifdef MNH_MGSOLVER
 IF ( HPRESOPT == 'ZSOLV' ) THEN
-  B_K(IKE+1) =  0.0
+  ZB_K(IKE+1) =  0.0
 
-  PAF_ZS(IIB:IIE,IJB:IJE,IKE+1) = -0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKE) + PRHO_ZS(IIB:IIE,IJB:IJE,IKE+1)   ) &
-                                       / ZDZM_ZS(IIB:IIE,IJB:IJE,IKE+1)   **2 
-  PCF_ZS(IIB:IIE,IJB:IJE,IKB-1) =  0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKB-1) + PRHO_ZS(IIB:IIE,IJB:IJE,IKB)   ) &
-                                       / ZDZM_ZS(IIB:IIE,IJB:IJE,IKB) **2 
-
-  PAF_ZS(IIB:IIE,IJB:IJE,IKB-1) = 0.0
-  PCF_ZS(IIB:IIE,IJB:IJE,IKE+1) = 0.0
-
-   if ( mppdb_initialized ) then  
-      call Mppdb_check( A_K,       "Tridz zsolv beg:A_K"       )
-      call Mppdb_check( B_K,       "Tridz zsolv beg:B_K"       )
-      call Mppdb_check( C_K,       "Tridz zsolv beg:C_K"       )
-      call Mppdb_check( D_K,       "Tridz zsolv beg:D_K"       )
+   if ( mppdb_initialized ) then
+      call Mppdb_check( ZA_K,       "Tridz zsolv beg:ZA_K"       )
+      call Mppdb_check( ZB_K,       "Tridz zsolv beg:ZB_K"       )
+      call Mppdb_check( ZC_K,       "Tridz zsolv beg:ZC_K"       )
+      call Mppdb_check( ZD_K,       "Tridz zsolv beg:ZD_K"       )
    end if
    call mg_main_mnh_init(IIMAX_ll,IKU,PDXHATM*IIMAX_ll,ZDZM(IKB)*IKU,&
-        A_K,B_K,C_K,D_K)
+        ZA_K,ZB_K,ZC_K,ZD_K)
 END IF
 #endif
 !------------------------------------------------------------------------------
@@ -774,9 +709,6 @@ ELSE
 !JUAN Z_SPLITTING
 !JUAN for Z splitting we need to do the vertical substitution
 !JUAN in Bsplitting slides so need PBFB 
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) PBF_ZS(:,:,:) = 1.0
-#endif
   IF ( HPRESOPT == 'ZRESI' .OR. HPRESOPT == 'ZSOLV' ) THEN
   DO JK= IKB,IKE
     DO JJ= IJB,IJE
@@ -788,19 +720,6 @@ ELSE
     END DO
   END DO
   END IF
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    DO JK= IKB,IKE
-      DO JJ= IJB,IJE
-        DO JI= IIB, IIE
-          PBF_ZS(JI,JJ,JK) = PRHO_ZS(JI,JJ,JK)* ( -2.0 / PDXATH_ZS(JI,JJ)**2  -2.0 /PDYATH_ZS(JI,JJ)**2  )   - 0.5 * &
-                             ( ( PRHO_ZS(JI,JJ,JK-1) + PRHO_ZS(JI,JJ,JK)   ) / ZDZM_ZS(JI,JJ,JK)  **2                &
-                              +( PRHO_ZS(JI,JJ,JK)   + PRHO_ZS(JI,JJ,JK+1) ) / ZDZM_ZS(JI,JJ,JK+1)**2 )
-        END DO
-      END DO
-    END DO
-  END IF
-#endif
 ! at the upper and lower levels PBFB is computed using the Neumann
 ! condition
 !
@@ -809,15 +728,6 @@ ELSE
   !
   PBFB(IIB:IIE,IJB:IJE,IKE+1) = + 0.5 * ( PRHOM(IKE)   + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2
   END IF
-#ifdef MNH_MGSOLVER
-  IF ( HPRESOPT == 'ZSOLV' ) THEN
-    PBF_ZS(IIB:IIE,IJB:IJE,IKB-1) = - 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKB-1) + PRHO_ZS(IIB:IIE,IJB:IJE,IKB)   ) &
-                                    / ZDZM_ZS(IIB:IIE,IJB:IJE,IKB)**2
-
-    PBF_ZS(IIB:IIE,IJB:IJE,IKE+1) =   0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKE)   + PRHO_ZS(IIB:IIE,IJB:IJE,IKE+1) ) &
-                                    / ZDZM_ZS(IIB:IIE,IJB:IJE,IKE+1)**2
-  END IF
-#endif
 !
 IF (HLBCX(1) == 'CYCL' .AND. .NOT.(L2D) ) THEN
    !JUAN