diff --git a/MNH/advec_4th_order_aux.f90 b/MNH/advec_4th_order_aux.f90
index 559928dec68883f10ee9e2d4908b74d950ac7d25..9041e00ebefef30b68ad7af318522fe9a458cbfc 100644
--- a/MNH/advec_4th_order_aux.f90
+++ b/MNH/advec_4th_order_aux.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/advec_4th_order_aux.f90,v $ $Revision: 1.1.2.1.2.2 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/advec_4th_order_aux.f90,v $ $Revision: 1.1.2.1.2.2.18.3 $
 ! MASDEV4_7 adiab 2006/05/18 14:35:31
 !-----------------------------------------------------------------
 !     ###############################
@@ -99,6 +103,8 @@ END MODULE MODI_ADVEC_4TH_ORDER_AUX
 !!    MODIFICATIONS
 !!    -------------
 !!      Original   25/10/05
+!!      Modif
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !!
 !-------------------------------------------------------------------------------
 !
@@ -185,21 +191,21 @@ ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
 !$acc update device (ZHALO2_WEST,ZHALO2_EAST)
 !
 !$acc kernels
-  IF(NHALO == 1) THEN
+!!$  IF(NHALO == 1) THEN
     IW=IIB+1
     IE=IIE
 !    IE=IIE-1
-  ELSE
-    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
-    WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
-    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
-    !callabortstop
-    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
-    CALL ABORT
-    STOP
-!    IW=IIB
-!    IE=IIE
-  END IF  
+!!$  ELSE
+!!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
+!!$    WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
+!!$    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
+!!$    !callabortstop
+!!$    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+!!$    CALL ABORT
+!!$    STOP
+!!$!    IW=IIB
+!!$!    IE=IIE
+!!$  END IF  
 !
   IF(KGRID == 2) THEN
     IWF=IW-1
@@ -237,13 +243,14 @@ ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
       IW=IIB+1
     END IF
   ELSE
-    IF(NHALO == 1) THEN
+!!$    IF(NHALO == 1) THEN
       IW=IIB+1
-    ELSE
-      IW=IIB
-    ENDIF
+!!$    ELSE
+!!$      IW=IIB
+!!$    ENDIF
   ENDIF
-  IF (GEAST .OR. NHALO == 1) THEN
+!!$  IF (GEAST .OR. NHALO == 1) THEN
+  IF (GEAST) THEN
 ! T. Maric
 !    IE=IIE-1 ! original
     IE=IIE
@@ -271,14 +278,16 @@ ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
     ! PMEANX(1,:,:) = PMEANX(IWF-1,:,:)
     ! extrapolate
     !PMEANX(1,:,:) = 0.5*(3.0*PFIELDT(1,:,:) - PFIELDT(2,:,:))
-  ELSEIF (NHALO == 1) THEN
+!!$  ELSE IF (NHALO == 1) THEN
+  ELSE
     PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) -  &
                              ( PFIELDT(IW,:,:)+ZHALO2_WEST(:,:) ) )/12.0
   ENDIF
 !
   IF (GEAST) THEN
     PMEANX(IEF+1,:,:) = 0.5*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) )
-  ELSEIF (NHALO == 1) THEN
+!!$  ELSEIF (NHALO == 1) THEN
+  ELSE
     PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) -  &
                          ( ZHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
   ENDIF
@@ -309,21 +318,21 @@ ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
 !$acc kernels
 !
 !
-    IF(NHALO == 1) THEN
+!!$    IF(NHALO == 1) THEN
       IS=IJB+1
       IN=IJE
 !      IN=IJE-1
-    ELSE
-      CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
-      WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
-      WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
-!callabortstop
-      CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
-      CALL ABORT
-      STOP
-!      IS=IJB
-!      IN=IJE
-    END IF
+!!$    ELSE
+!!$      CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
+!!$      WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
+!!$      WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
+!!$!callabortstop
+!!$      CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+!!$      CALL ABORT
+!!$      STOP
+!!$!      IS=IJB
+!!$!      IN=IJE
+!!$    END IF
 !
     IF(KGRID == 3) THEN
       ISF=IS-1
@@ -361,13 +370,14 @@ ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
         IS=IJB+1
       END IF
     ELSE
-      IF(NHALO == 1) THEN
+!!$      IF(NHALO == 1) THEN
         IS=IJB+1
-      ELSE
-        IS=IJB
-      ENDIF
+!!$      ELSE
+!!$        IS=IJB
+!!$      ENDIF
     ENDIF
-    IF (GNORTH .OR. NHALO == 1) THEN
+!!$    IF (GNORTH .OR. NHALO == 1) THEN
+    IF (GNORTH) THEN
 ! T. Maric
 !      IN=IJE-1  ! original
       IN=IJE
@@ -391,14 +401,20 @@ ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
       ! PMEANY(:,1,:) = PMEANY(:,ISF-1,:)
       ! extrapolate
       !PMEANY(:,1,:) = 0.5*(3.0*PFIELDT(:,1,:) - PFIELDT(:,2,:))
-    ELSEIF (NHALO == 1) THEN
+!!$    ELSEIF (NHALO == 1) THEN
+    ELSE
+!!$      PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:)) -  &
+!!$                          ( PFIELDT(:,IS+1,:)+TPHALO2%SOUTH(:,:) ))/12.0
        PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:)) -  &
                            ( PFIELDT(:,IS,:)+ZHALO2_SOUTH(:,:) ))/12.0
     ENDIF
 !
     IF (GNORTH) THEN
       PMEANY(:,INF+1,:) = 0.5*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:) )
-    ELSEIF (NHALO == 1) THEN
+!!$    ELSEIF (NHALO == 1) THEN
+    ELSE
+!!$      PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN,:)+PFIELDT(:,IN-1,:)) -  &
+!!$                          ( TPHALO2%NORTH(:,:)+PFIELDT(:,IN-2,:) ))/12.0
        PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:)) -  &
                            ( ZHALO2_NORTH(:,:)+PFIELDT(:,IN-1,:) ))/12.0
     ENDIF
diff --git a/MNH/advec_ppm_algo.f90 b/MNH/advec_ppm_algo.f90
index 834360306225394b154ed49dd291ebe24330016b..3b95d189003e2cb934b06080bea096c7b73edfd0 100644
--- a/MNH/advec_ppm_algo.f90
+++ b/MNH/advec_ppm_algo.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/advec_ppm_algo.f90,v $ $Revision: 1.1.2.1.2.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/advec_ppm_algo.f90,v $ $Revision: 1.1.2.1.2.1.18.2 $
 ! MASDEV4_7 adiab 2007/03/27 10:07:52
 !-----------------------------------------------------------------
 !     ##########################
diff --git a/MNH/advection.f90 b/MNH/advection.f90
index 3bab126ca2137050303af285aaccd054a1bca38f..363bd919cf0e3357588eecae9c4120610ff8b6e8 100644
--- a/MNH/advection.f90
+++ b/MNH/advection.f90
@@ -1,3 +1,7 @@
+!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.
 !-----------------------------------------------------------------
 !     #####################
       MODULE MODI_ADVECTION
@@ -10,10 +14,8 @@ INTERFACE
                            PUM, PVM, PWM, PTHM, PRM, PTKEM, PSVM,              &
                            PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT,              &
                            PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY,               &
-                           PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS,        &
-                           TPHALO2MLIST, TPHALO2LIST, TPHALO2SLIST )
+                           PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS         )
 !
-USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 !
 CHARACTER(LEN=6),         INTENT(IN)    :: HMET_ADV_SCHEME, & ! Control of the 
                                            HSV_ADV_SCHEME,  & ! scheme applied 
@@ -50,10 +52,6 @@ REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS, PRTKES
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS
                                                   ! Sources terms 
 !
-! halo lists for 4th order advection
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2MLIST ! momentum variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST  ! meteorological scalar variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2SLIST ! tracer scalar variables
 !
 END SUBROUTINE ADVECTION
 !
@@ -67,93 +65,7 @@ END MODULE MODI_ADVECTION
                            PUM, PVM, PWM, PTHM, PRM, PTKEM, PSVM,              &
                            PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT,              &
                            PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY,               &
-                           PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS,        &
-                           TPHALO2MLIST, TPHALO2LIST, TPHALO2SLIST )
-
-  USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
-  USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
-
-IMPLICIT NONE
-!
-!*       0.1   Declarations of dummy arguments :
-!
-CHARACTER(LEN=6),         INTENT(IN)    :: HMET_ADV_SCHEME, & ! Control of the 
-                                           HSV_ADV_SCHEME,  & ! scheme applied 
-                                           HUVW_ADV_SCHEME     ! to the selected
-                                                              ! variables 
-!
-INTEGER,                  INTENT(IN)    :: KLITER        ! Iteration number for
-                                                         ! the MPDATA scheme
-!
-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
-!
-INTEGER,                  INTENT(IN)    :: KTCOUNT  ! iteration count
-REAL,                     INTENT(IN)    :: PTSTEP_MET !  Effective time step for
-                                                ! meteorological scalar variables 
-                                                ! (depending on advection scheme)
-REAL,                     INTENT(IN)    :: PTSTEP_SV !  Effective time step for
-                                                ! tracer scalar variables 
-                                                ! (depending on advection scheme)
-!
-REAL, DIMENSION(:,:,:),   INTENT(IN) :: PTHM, PTKEM
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM , PSVM
-                                                  ! Variables at t-dt
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUM, PVM, PWM
-                                                  ! Variables at t-dt
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT , PVT  , PWT
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT, PTKET, PRHODJ
-REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT , PSVT
-                                                  ! Variables at t
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PDXX,PDYY,PDZZ,PDZX,PDZY
-                                                  !  metric coefficients
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRUS , PRVS, PRWS
-REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS, PRTKES
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS
-                                                  ! Sources terms 
-! halo lists for 4th order advection
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2MLIST ! momentum variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST  ! meteorological scalar variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2SLIST ! tracer scalar variables
-
-INTEGER :: IZRUT,IZRVT,IZRWT,IZRUCT,IZRVCT,IZRWCT
-
-         CALL  MNH_GET_ZT3D(IZRUT,IZRVT,IZRWT,IZRUCT,IZRVCT,IZRWCT)
-
-         CALL ADVECTION_D (HUVW_ADV_SCHEME,HMET_ADV_SCHEME,HSV_ADV_SCHEME,     &
-                        &  KLITER, HLBCX, HLBCY,KRR, KSV, KTCOUNT,             &
-                        &  PTSTEP_MET, PTSTEP_SV,                              & 
-                        &  PUM, PVM, PWM, PTHM, PRM, PTKEM, PSVM,              &
-                        &  PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT,              &
-                        &  PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY,               &
-                        &  PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS,        &
-                        &  TPHALO2MLIST, TPHALO2LIST, TPHALO2SLIST             &
-#define USET3D
-#ifdef USET3D
-                        &  , ZT3D(:,:,:,IZRUT),ZT3D(:,:,:,IZRVT),ZT3D(:,:,:,IZRWT) &
-                        &  ,ZT3D(:,:,:,IZRUCT),ZT3D(:,:,:,IZRVCT),ZT3D(:,:,:,IZRWCT)  &
-#endif
-                        &  )
-
-         CALL  MNH_REL_ZT3D(IZRUT,IZRVT,IZRWT,IZRUCT,IZRVCT,IZRWCT)
-CONTAINS   
-
-      SUBROUTINE ADVECTION_D (HUVW_ADV_SCHEME,HMET_ADV_SCHEME,HSV_ADV_SCHEME,  &
-                           KLITER, HLBCX, HLBCY,KRR, KSV, KTCOUNT,             &
-                           PTSTEP_MET, PTSTEP_SV,                              & 
-                           PUM, PVM, PWM, PTHM, PRM, PTKEM, PSVM,              &
-                           PUT, PVT, PWT, PTHT, PRT, PTKET, PSVT,              &
-                           PRHODJ, PDXX, PDYY, PDZZ, PDZX, PDZY,               &
-                           PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS,        &
-                           TPHALO2MLIST, TPHALO2LIST, TPHALO2SLIST             &
-#ifdef USET3D
-                           ,ZRUT,ZRVT,ZRWT                                     &
-                            ,ZRUCT,ZRVCT,ZRWCT                   &
-#endif
-                           )
-!
+                           PRUS,PRVS, PRWS, PRTHS, PRRS, PRTKES, PRSVS         )
 !     ##########################################################################
 !
 !!****  *ADVECTION * - routine to call the specialized advection routines
@@ -217,29 +129,7 @@ CONTAINS
 !*       0.    DECLARATIONS
 !              ------------
 !
-USE MODE_ll
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
-USE MODD_CONF
-USE MODD_BLANK
-!
-USE MODI_SHUMAN
-USE MODI_CONTRAV
-USE MODI_ADVECUVW
-USE MODI_ADVECUVW_4TH
-USE MODI_ADVECMET      
-USE MODI_ADVECMET_4TH
-USE MODI_FCT_MET   
-USE MODI_MPDATA
-USE MODI_ADVECSCALAR
-USE MODI_ADVECSCALAR_4TH
-USE MODI_FCT_SCALAR 
-USE MODI_MPDATA_SCALAR
-USE MODI_PPM_MET
-USE MODI_PPM_SCALAR
-!
-USE MODI_GET_HALO
 !
-USE MODD_PARAMETERS , ONLY : JPHEXT
 !
 !-------------------------------------------------------------------------------
 !
@@ -272,14 +162,9 @@ REAL, DIMENSION(:,:,:),   INTENT(IN) :: PTHM, PTKEM
 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM , PSVM
                                                   ! Variables at t-dt
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUM, PVM, PWM
-!$acc declare pcopyin(PUM, PVM, PWM)
                                                   ! Variables at t-dt
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT , PVT  , PWT
-!$acc declare pcopyin(PUT, PVT, PWT)
-
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT, PTKET, PRHODJ
-!$acc declare present(PRHODJ)
-
 REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT , PSVT
                                                   ! Variables at t
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PDXX,PDYY,PDZZ,PDZX,PDZY
@@ -288,230 +173,9 @@ REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRUS , PRVS, PRWS
 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS, PRTKES
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS , PRSVS
                                                   ! Sources terms 
-! halo lists for 4th order advection
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2MLIST ! momentum variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST  ! meteorological scalar variables
-TYPE(HALO2LIST_ll), POINTER :: TPHALO2SLIST ! tracer scalar variables
-!
-INTEGER                                              :: IIU,IJU,IKU
-!
-!*       0.2   declarations of local variables
-!
-!
-!  
-#ifdef USET3D
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRUT 
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRVT 
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRWT
-!$acc declare present (ZRUT,ZRVT,ZRWT)
-                                                  ! cartesian 
-                                                  ! components of
-                                                  ! momentum
-!
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRUCT 
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRVCT
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRWCT
-!$acc declare present (ZRUCT,ZRVCT,ZRWCT)
-#endif 
-                                                  ! contravariant
-                                                  ! components
-                                                  ! of momentum
-!
-INTEGER                     :: IINFO_ll    ! return code of parallel routine
-TYPE(LIST_ll), POINTER      :: TZFIELDS_ll ! list of fields to exchange
-!
-#define mxm(PMXM,PA) PMXM(2:IIU,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:) ! MXM(PMXM,PA)
-#define mym(PMYM,PA) PMYM(:,2:IJU,:) = 0.5*( PA(:,2:IJU,:)+PA(:,1:IJU-1,:) ) ; PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:) ! MYM(PMYM,PA)
-#define mzm(PMZM,PA) PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) ) ; PMZM(:,:,1)    = -999. ! MZM(PMZM,PA)
-!
-!-------------------------------------------------------------------------------
-!
-!*       0.     COMPUTES THE DOMAIN DIMENSIONS
-!               ------------------------------
-CALL GET_DIM_EXT_ll('B',IIU,IJU)     
-IKU=SIZE(PUT,3)
-!
-!*       1.     COMPUTES THE CONTRAVARIANT COMPONENTS
-!	        -------------------------------------
-!
-!!$ZRUT = PUT(:,:,:) * MXM(PRHODJ)
-!!$ZRVT = PVT(:,:,:) * MYM(PRHODJ)
-!!$ZRWT = PWT(:,:,:) * MZM(PRHODJ)
-!$acc update device(PRHODJ)
-!$acc kernels 
-mxm(ZRUT,PRHODJ)
-ZRUT = PUT(:,:,:) * ZRUT
-mym(ZRVT,PRHODJ)
-ZRVT = PVT(:,:,:) * ZRVT
-mzm(ZRWT,PRHODJ)
-ZRWT = PWT(:,:,:) * ZRWT
-!$acc end kernels
-!
-CALL CONTRAV (1,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT)
-!
-NULLIFY(TZFIELDS_ll)
-IF(NHALO == 1) THEN
-  CALL GET_HALO_D(ZRUCT)
-  CALL GET_HALO_D(ZRVCT)
-  CALL GET_HALO_D(ZRWCT) 
-!!$  CALL ADD3DFIELD_ll(TZFIELDS_ll, ZRUCT)
-!!$  CALL ADD3DFIELD_ll(TZFIELDS_ll, ZRVCT)
-!!$  CALL ADD3DFIELD_ll(TZFIELDS_ll, ZRWCT)
-!!$  CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
-!!$  CALL CLEANLIST_ll(TZFIELDS_ll)
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!*       2.     CALLS THE ADVECTION ROUTINES FOR THE MOMENTUM 
-!	        ---------------------------------------------
-!
-! choose between 2nd and 4th order momentum advection.
-IF (HUVW_ADV_SCHEME=='CEN2ND' ) THEN                                      
-!
-   CALL ADVECUVW (PUT,PVT,PWT,ZRUCT,ZRVCT,ZRWCT,PRUS,PRVS,PRWS)
-!
-ELSEIF (HUVW_ADV_SCHEME=='CEN4TH') THEN
-! 
-   CALL ADVECUVW_4TH ( HLBCX, HLBCY, ZRUCT, ZRVCT, ZRWCT,            &
-                       PUT, PVT, PWT, PRUS, PRVS, PRWS, TPHALO2MLIST )                 
-!
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!*       3.     CALLS THE ADVECTION ROUTINES FOR THE METEOROLOGICAL SCALARS 
-!	        -----------------------------------------------------------
-!
-!            3.1. 2nd order scheme
-!
-IF (HMET_ADV_SCHEME=='CEN2ND') THEN
-!
-   CALL ADVECMET (KRR, PTHT,PRT,PTKET,     &
-                  ZRUCT,ZRVCT,ZRWCT,       &
-                  PRTHS,PRRS,PRTKES        )
-!
-!             3.2. 4th order scheme
-!
-ELSEIF (HMET_ADV_SCHEME =='CEN4TH' ) THEN
-!
-   CALL ADVECMET_4TH (HLBCX,HLBCY, KRR,                &
-                      ZRUCT, ZRVCT, ZRWCT,             &
-                      PTHT, PTKET, PRT,                &
-                      PRTHS, PRTKES, PRRS, TPHALO2LIST )
-!
-!             3.3. Flux-Corrected Transport scheme
-!
-ELSEIF ( HMET_ADV_SCHEME=='FCT2ND') THEN
-!
-   CALL FCT_MET  (HLBCX, HLBCY,KRR,                        &
-                  PTSTEP_MET, PRHODJ, PTHM, PRM, PTKEM,    &
-                  PTHT, PRT, PTKET,                        &
-                  ZRUCT, ZRVCT, ZRWCT,                     &
-                  PRTHS, PRRS, PRTKES                      )
-!
-!             3.4. MPDATA scheme
-!
-ELSEIF (HMET_ADV_SCHEME=='MPDATA') THEN
-!
-   CALL MPDATA (KLITER, HLBCX, HLBCY, KRR,                 &
-                PTSTEP_MET, PRHODJ, PTHM, PRM, PTKEM,      &
-                PTHT, PRT, PTKET,                          &
-                ZRUCT, ZRVCT, ZRWCT,                       &
-                PRTHS, PRRS, PRTKES                        )
-!
-!             3.5. PPM schemes
-!   
-ELSEIF (HMET_ADV_SCHEME(1:3)=='PPM') THEN
-!
-! extrapolate velocity field to t+dt/2 to use in forward in time PPM
-! advection scheme
-!
-!$acc kernels 
-   ZRUT = (1.5*PUT(:,:,:) - 0.5*PUM(:,:,:))
-   ZRVT = (1.5*PVT(:,:,:) - 0.5*PVM(:,:,:))
-   ZRWT = (1.5*PWT(:,:,:) - 0.5*PWM(:,:,:))
-!$acc end kernels 
-! calculate Courant numbers
-   CALL CONTRAV (1,ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT)
-!
-!$acc kernels 
-   ZRUCT = ZRUCT*PTSTEP_MET
-   ZRVCT = ZRVCT*PTSTEP_MET
-   ZRWCT = ZRWCT*PTSTEP_MET
-!$acc end kernels 
-
-   !$acc update device(PRT,PTHT,PRTHS,PRRS,PRHODJ)
-   CALL PPM_MET   (HLBCX,HLBCY, KRR, KTCOUNT,                  &
-                   ZRUCT, ZRVCT, ZRWCT, PTSTEP_MET, PRHODJ,    &
-                   PTHT, PTKET, PRT, PRTHS, PRTKES, PRRS,      &
-                   HMET_ADV_SCHEME                             )
-   !$acc update host(PRTHS,PRRS)
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!*       4.     CALLS THE ADVECTION ROUTINES FOR TRACERS
-!	        ----------------------------------------
-!
-!            4.1. 2nd order scheme
-!
-IF (KSV .NE. 0) THEN
-!
-IF (HSV_ADV_SCHEME=='CEN2ND') THEN
-!
-   CALL ADVECSCALAR  (KSV, PSVT, ZRUCT,ZRVCT,ZRWCT,PRSVS )             
-!
-!             4.2. 4th order scheme
-!
-ELSEIF (HSV_ADV_SCHEME =='CEN4TH' ) THEN
-!
-   CALL ADVECSCALAR_4TH (HLBCX,HLBCY, KSV,          &
-                         ZRUCT, ZRVCT, ZRWCT,       &
-                         PSVT, PRSVS, TPHALO2SLIST   )           
-!
-!             4.3. Flux-Corrected Transport scheme
-!
-ELSEIF ( HSV_ADV_SCHEME=='FCT2ND') THEN
-!
-   CALL FCT_SCALAR  (HLBCX, HLBCY, KSV,             &
-                     PTSTEP_SV, PRHODJ, PSVM,PSVT,  &
-                     ZRUCT, ZRVCT, ZRWCT, PRSVS     ) 
-!
-!             4.4. MPDATA scheme
-!
-ELSEIF (HSV_ADV_SCHEME=='MPDATA') THEN
-!
-   CALL MPDATA_SCALAR ( KLITER, HLBCX, HLBCY, KSV,           &
-                        PTSTEP_SV, PRHODJ, PSVM, PSVT,       &
-                        ZRUCT, ZRVCT, ZRWCT,  PRSVS          )           
-!
-!             4.5. PPM schemes
-!   
-ELSEIF (HSV_ADV_SCHEME(1:3)=='PPM') THEN
-!
-! extrapolate velocity field to t+dt/2 to use in forward in time PPM
-! advection scheme
-!
-   ZRUT = (1.5*PUT(:,:,:) - 0.5*PUM(:,:,:))
-   ZRVT = (1.5*PVT(:,:,:) - 0.5*PVM(:,:,:))
-   ZRWT = (1.5*PWT(:,:,:) - 0.5*PWM(:,:,:))
-! calculate Courant numbers
-   CALL CONTRAV (ZRUT,ZRVT,ZRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,ZRUCT,ZRVCT,ZRWCT)
-!
-   ZRUCT = ZRUCT*PTSTEP_SV
-   ZRVCT = ZRVCT*PTSTEP_SV
-   ZRWCT = ZRWCT*PTSTEP_SV
-
-   CALL PPM_SCALAR(HLBCX,HLBCY, KSV, KTCOUNT,               &
-                   ZRUCT, ZRVCT, ZRWCT, PTSTEP_SV, PRHODJ,  &
-                   PSVT, PRSVS, HSV_ADV_SCHEME      )                     
 !
-END IF
-END IF
 !
+! ROUTINE TO REMOVE
 !-------------------------------------------------------------------------------
 !
-END SUBROUTINE ADVECTION_D
-
 END SUBROUTINE ADVECTION
diff --git a/MNH/advecuvw_4th.f90 b/MNH/advecuvw_4th.f90
index fecb024145dc63cf61404bca7017b468514262d1..2bcf3dabb8da5559e51523f569040c1bbbd9c903 100644
--- a/MNH/advecuvw_4th.f90
+++ b/MNH/advecuvw_4th.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/advecuvw_4th.f90,v $ $Revision: 1.1.4.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/advecuvw_4th.f90,v $ $Revision: 1.1.4.1.16.1.2.4 $
 ! MASDEV4_7 adiab 2006/05/18 13:07:25
 !-----------------------------------------------------------------
 !     ###########################
@@ -140,48 +144,6 @@ CONTAINS
 !!         NBUPROCCTR   : process counter used for each budget variable
 !!         Switches for budgets activations:
 !!
-!!         LBU_RU       : logical for budget of RU (wind component along x)
-!!
-!!         LBU_RU       : logical for budget of RU (wind component along x)
-!!                        .TRUE. = budget of RU
-!!                        .FALSE. = no budget of RU
-!!         LBU_RV       : logical for budget of RV (wind component along y)
-!!                        .TRUE. = budget of RV
-!!                        .FALSE. = no budget of RV
-!!         LBU_RW        : logical for budget of RW (wind component along z)
-!!                        .TRUE. = budget of RW
-!!                        .FALSE. = no budget of RW
-!!         LBU_RTH      : logical for budget of RTH (potential temperature)
-!!                        .TRUE. = budget of RTH
-!!                        .FALSE. = no budget of RTH
-!!         LBU_RTKE     : logical for budget of RTKE (turbulent kinetic energy)
-!!                        .TRUE. = budget of RTKE
-!!                        .FALSE. = no budget of RTKE
-!!         LBU_RRV      : logical for budget of RRV (water vapor)
-!!                        .TRUE. = budget of RRV
-!!                        .FALSE. = no budget of RRV
-!!         LBU_RRC      : logical for budget of RRC (cloud water)
-!!                        .TRUE. = budget of RRC
-!!                        .FALSE. = no budget of RRC
-!!         LBU_RRR      : logical for budget of RRR (rain water)
-!!                        .TRUE. = budget of RRR
-!!                        .FALSE. = no budget of RRR
-!!         LBU_RRI      : logical for budget of RRI (ice)
-!!                        .TRUE. = budget of RRI
-!!                        .FALSE. = no budget of RRI
-!!         LBU_RRS      : logical for budget of RRS (snow)
-!!                        .TRUE. = budget of RRS
-!!                        .FALSE. = no budget of RRS
-!!         LBU_RRG      : logical for budget of RRG (graupel)
-!!                        .TRUE. = budget of RRG
-!!                        .FALSE. = no budget of RRG
-!!         LBU_RRH      : logical for budget of RRH (hail)
-!!                        .TRUE. = budget of RRH
-!!                        .FALSE. = no budget of RRH
-!!         LBU_RSV      : logical for budget of RSVx (scalar variable)
-!!                        .TRUE. = budget of RSVx
-!!                        .FALSE. = no budget of RSVx
-!!
 !!    MODULE MODD_ARGSLIST
 !!         HALO2LIST_ll : type for a list of "HALO2_lls"
 !!
@@ -196,6 +158,8 @@ CONTAINS
 !!    MODIFICATIONS
 !!    -------------
 !!      Original   25/10/05
+!!      Modif
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !!
 !-------------------------------------------------------------------------------
 !
@@ -206,11 +170,10 @@ USE MODE_ll
 !
 USE MODD_PARAMETERS
 USE MODD_CONF
-USE MODD_BUDGET
+USE MODD_GRID_n
 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 !
 USE MODI_SHUMAN
-USE MODI_BUDGET
 !
 USE MODI_ADVEC_4TH_ORDER_AUX
 !
@@ -241,6 +204,7 @@ TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST ! list for diffusion
 !
 INTEGER:: IIB,IJB        ! Begining useful area  in x,y,z directions
 INTEGER:: IIE,IJE        ! End useful area in x,y,z directions
+INTEGER :: IKU
 !
 TYPE(HALO2LIST_ll), POINTER :: TZHALO2LIST
 !
@@ -283,20 +247,22 @@ CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
 !!$CALL GET_DIM_EXT_ll('B',IIU,IJU)     
 !!$IKU=SIZE(PUT,3)
 !
+IKU=SIZE(XZHAT)
 !-------------------------------------------------------------------------------
 !
 !*       2.     CALL THE ADVEC_4TH_ORDER_ALGO ROUTINE FOR MOMENTUM
 !               --------------------------------------------------
 !
 IGRID = 2
-IF(NHALO == 1) THEN
+!!$IF(NHALO == 1) THEN
   TZHALO2LIST => TPHALO2LIST
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY, &
                             TZHALO2LIST%HALO2 )
-ELSE
-  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY)
-ENDIF
+!!$ELSE
+!!$  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
 !
+
 ! pcopy(prus) pcopyin(pruct,ZMEANX) create(ZTEMP1,ZTEMP2,ZTEMP3)
 !!$PRUS(:,:,:) = PRUS(:,:,:)                          &
 !!$             -DXM( MXF(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
@@ -308,9 +274,8 @@ dxm(ZTEMP3,ZTEMP2)
 PRUS(:,:,:) = PRUS(:,:,:) -  ZTEMP3
 !$acc end kernels   
          
-IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVX_BU_RU')
-
 !
+
 !!$PRUS(:,:,:) = PRUS(:,:,:)                          &
 !!$             -DYF( MXM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
 
@@ -321,8 +286,8 @@ dyf(ZTEMP3,ZTEMP2)
 PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP3
 !$acc end kernels            
 
-IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVY_BU_RU')
 !
+
 !!$PRUS(:,:,:) = PRUS(:,:,:)                             &
 !!$             -DZF( MXM(PRWCT(:,:,:))*MZM4(PUT(:,:,:)) )
 
@@ -334,18 +299,18 @@ dzf(ZTEMP4,ZTEMP3)
 PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP4
 !$acc end kernels            
 
-IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVZ_BU_RU')
 !
 !
 IGRID = 3
-IF(NHALO == 1) THEN
+!!$IF(NHALO == 1) THEN
   TZHALO2LIST => TZHALO2LIST%NEXT
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY, &
                             TZHALO2LIST%HALO2 )
-ELSE
-  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY)
-ENDIF
+!!$ELSE
+!!$  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
 !
+
 !!$PRVS(:,:,:) = PRVS(:,:,:)                          &
 !!$             -DXF( MYM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
 
@@ -356,8 +321,9 @@ dxf(ZTEMP3,ZTEMP2)
 PRVS(:,:,:) = PRVS(:,:,:) -  ZTEMP3   
 !$acc end kernels   
                      
-IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVX_BU_RV')
+
 !
+
 !!$PRVS(:,:,:) = PRVS(:,:,:)                          &
 !!$             -DYM( MYF(PRVCT(:,:,:))*ZMEANY(:,:,:) )
 
@@ -368,8 +334,8 @@ dym(ZTEMP3,ZTEMP2)
 PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP3
 !$acc end kernels   
 
-IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVY_BU_RV')
 !
+
 !!$PRVS(:,:,:) = PRVS(:,:,:)                             &
 !!$             -DZF( MYM(PRWCT(:,:,:))*MZM4(PVT(:,:,:)) )
 
@@ -381,19 +347,19 @@ dzf(ZTEMP4,ZTEMP3)
 PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP4
 !$acc end kernels
 
-IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVZ_BU_RV')
 !
 !
 IGRID = 4
 !
-IF(NHALO == 1) THEN
+!!$IF(NHALO == 1) THEN
   TZHALO2LIST => TZHALO2LIST%NEXT
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY, &
                             TZHALO2LIST%HALO2 )
-ELSE
-  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY)
-ENDIF
+!!$ELSE
+!!$  CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY)
+!!$ENDIF
 !
+
 !!$PRWS(:,:,:) = PRWS(:,:,:)                          &
 !!$             -DXF( MZM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
 
@@ -404,8 +370,8 @@ dxf(ZTEMP3,ZTEMP2)
 PRWS(:,:,:) = PRWS(:,:,:)  - ZTEMP3
 !$acc end kernels
 
-IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVX_BU_RW')
 !
+
 !!$PRWS(:,:,:) = PRWS(:,:,:)                          &
 !!$             -DYF( MZM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
 
@@ -416,8 +382,8 @@ dyf(ZTEMP3,ZTEMP2)
 PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP3
 !$acc end kernels
 
-IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVY_BU_RW')
 !
+
 !!$PRWS(:,:,:) = PRWS(:,:,:)                             &
 !!$             -DZM( MZF(PRWCT(:,:,:))*MZF4(PWT(:,:,:)) )
 
@@ -428,7 +394,7 @@ ZTEMP1 = ZTEMP1 * ZTEMP2
 dzm(ZTEMP4,ZTEMP1)
 PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP4
 !$acc end kernels
-IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVZ_BU_RW')
+
 !
 !acc end data 
 !-------------------------------------------------------------------------------
diff --git a/MNH/contrav.f90 b/MNH/contrav.f90
index ff2eefea39758c36b5f4139cc44ca58ea5f5568c..b785b04225e812c64e1523de305c96701ba572e8 100644
--- a/MNH/contrav.f90
+++ b/MNH/contrav.f90
@@ -1,48 +1,41 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/contrav.f90,v $ $Revision: 1.1.10.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/contrav.f90,v $ $Revision: 1.1.10.1.16.2.2.3 $
 ! MASDEV4_7 operators 2006/05/18 13:07:25
 !-----------------------------------------------------------------
 !     ####################
       MODULE MODI_CONTRAV
 !     ####################
 !
-INTERFACE CONTRAV
+INTERFACE
 !
-      SUBROUTINE CONTRAV_CPU (PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
-                         PRUCT,PRVCT,PRWCT                       )
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRUT     ! Cartesian comp along x
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRVT     ! Cartesian comp along y
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRWT     ! Cartesian comp along z
-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(OUT)   ::  PRUCT    ! Contrav comp along x-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRVCT    ! Contrav comp along y-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRWCT    ! Contrav comp along z-bar
-!
-END SUBROUTINE CONTRAV_CPU
-!
-      SUBROUTINE CONTRAV_ACC (IACC,PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
-                         PRUCT,PRVCT,PRWCT                       )
-INTEGER                                ::  IACC
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRUT     ! Cartesian comp along x
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRVT     ! Cartesian comp along y
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRWT     ! Cartesian comp along z
-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(OUT)   ::  PRUCT    ! Contrav comp along x-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRVCT    ! Contrav comp along y-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRWCT    ! Contrav comp along z-bar
-!
-END SUBROUTINE CONTRAV_ACC
+      SUBROUTINE CONTRAV(HLBCX,HLBCY,PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
+                         PRUCT,PRVCT,PRWCT,KADV_ORDER                       )
 
+
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
+CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type                         
+REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRUT       ! Cartesian comp along x
+REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRVT       ! Cartesian comp along y
+REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRWT       ! Cartesian comp along z
+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(OUT)   ::  PRUCT      ! Contrav comp along x-bar
+REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRVCT      ! Contrav comp along y-bar
+REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRWCT      ! Contrav comp along z-bar
+INTEGER,                 INTENT(IN)    ::  KADV_ORDER ! Order of the advection
+                                                      ! scheme
+!
+END SUBROUTINE CONTRAV
+!
 END INTERFACE
 !
 END MODULE MODI_CONTRAV 
@@ -50,8 +43,8 @@ END MODULE MODI_CONTRAV
 !
 !
 !     ##############################################################
-      SUBROUTINE CONTRAV_CPU(PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
-                         PRUCT,PRVCT,PRWCT                       )
+      SUBROUTINE CONTRAV(HLBCX,HLBCY,PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
+                         PRUCT,PRVCT,PRWCT,KADV_ORDER                        )
 !     ##############################################################
 !
 !!****  *CONTRAV * - computes the contravariant components from the
@@ -109,20 +102,32 @@ END MODULE MODI_CONTRAV
 !!      Original   27/07/94
 !!      Corrections 3/08/94 (by J.P. Lafore)
 !!      Corrections 17/10/94 (by J.P. Lafore) WC modified for w-advection
+!!      Corrections 19/01/11 (by J.P. Pinty) WC 4th order
+!!      Corrections 28/03/11 (by V.Masson) // of WC 4th order
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
 !              ------------
 USE MODD_CONF
 USE MODD_PARAMETERS
+USE MODD_GRID_n, ONLY: XZZ
 !
+USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
+USE MODE_ll
 !
 USE MODI_SHUMAN
+USE MODI_GET_HALO
+!
+USE MODE_MPPDB
 !
 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
+!
 REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRUT     ! Cartesian comp along x
 REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRVT     ! Cartesian comp along y
 REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRWT     ! Cartesian comp along z
@@ -137,211 +142,278 @@ REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PDZY     ! Metric coefficients
 REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRUCT    ! Contrav comp along x-bar
 REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRVCT    ! Contrav comp along y-bar
 REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRWCT    ! Contrav comp along z-bar
+INTEGER,                 INTENT(IN)    ::  KADV_ORDER ! Order of the advection
+                                                      ! scheme
 !
 !
 !*       0.2   declarations of local variables
 !              
 REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: Z1,Z2
 INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE
-!
+INTEGER :: IIU, IJU, IKU
+INTEGER:: IW,IE,IS,IN   ! Coordinate of forth order diffusion area
+!
+TYPE(LIST_ll),      POINTER :: TZFIELD_U, TZFIELD_V, TZFIELD_DZX, TZFIELD_DZY
+TYPE(HALO2LIST_ll), POINTER :: TZHALO2_U, TZHALO2_V, TZHALO2_DZX, TZHALO2_DZY
+INTEGER                     :: IINFO_ll
+!JUAN
+REAL          :: XPRECISION
 !-----------------------------------------------------------------------
 !
-!*       1.    Compute the contravariant components
-!              ------------------------------------
+!*       1.    Compute the horizontal contravariant components
+!              -----------------------------------------------
+!
+CALL MPPDB_CHECK3DM("contrav big ::PRU/V/WT",PRECISION,PRUT,PRVT,PRWT)                    
+!
+IIU= SIZE(PDXX,1)
+IJU= SIZE(PDXX,2)
+IKU= SIZE(PDXX,3)
 !
-IIB=2
-IJB=2
-IIE=SIZE(PDXX,1)-1
-IJE=SIZE(PDXX,2)-1
+CALL GET_INDICE_ll( IIB,IJB,IIE,IJE)
 !
 IKB=1+JPVEXT
-IKE=SIZE(PDXX,3) - JPVEXT
+IKE=IKU - JPVEXT
 !
 PRUCT(:,:,:) = PRUT(:,:,:) / PDXX(:,:,:)
 PRVCT(:,:,:) = PRVT(:,:,:) / PDYY(:,:,:)
-
-IF (LFLAT) THEN
-
-  PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)  
-  
-ELSE 
-
-  Z1(IIB:IIE,:,IKB:IKE+1)=                                             &
-      (PRUCT(IIB:IIE,:,IKB:IKE+1)+PRUCT(IIB:IIE,:,IKB-1:IKE) )         &
-       *PDZX(IIB:IIE,:,IKB:IKE+1) *0.25                                &
-     +(PRUCT(IIB+1:IIE+1,:,IKB:IKE+1)+PRUCT(IIB+1:IIE+1,:,IKB-1:IKE) ) &
-       *PDZX(IIB+1:IIE+1,:,IKB:IKE+1) *0.25   
-                      
-  Z2(:,IJB:IJE,IKB:IKE+1)=                                             &
-      (PRVCT(:,IJB:IJE,IKB:IKE+1)+PRVCT(:,IJB:IJE,IKB-1:IKE) )         &
-       *PDZY(:,IJB:IJE,IKB:IKE+1) *0.25                                &
-     +(PRVCT(:,IJB+1:IJE+1,IKB:IKE+1)+PRVCT(:,IJB+1:IJE+1,IKB-1:IKE) ) &
-       *PDZY(:,IJB+1:IJE+1,IKB:IKE+1) *0.25   
-
-  PRWCT=0.             
-  PRWCT(IIB:IIE,IJB:IJE,IKB:IKE+1) =                 &
-      (   PRWT(IIB:IIE,IJB:IJE,IKB:IKE+1)            &
-        - Z1(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
-        - Z2(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
-      ) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE+1)  
- 
-END IF
 !
-PRWCT(:,:,1) = - PRWCT(:,:,3)     ! Mirror hypothesis
-!
-!-----------------------------------------------------------------------
-!
-END SUBROUTINE CONTRAV_CPU
-!
-!
-!     ##############################################################
-      SUBROUTINE CONTRAV_ACC(IACC,PRUT,PRVT,PRWT,PDXX,PDYY,PDZZ,PDZX,PDZY,  &
-                         PRUCT,PRVCT,PRWCT                       )
-!     ##############################################################
-!
-!!****  *CONTRAV * - computes the contravariant components from the
-!!       cartesian components
-!!
-!!    PURPOSE
-!!    -------
-!       This routine computes the contravariant components of vector
-!     defined by its cartesian components (U,V,W) , using the following
-!     formulae:
-!     UC = U / DXX
-!     VC = V / DYY
-!               (     ----------x    ----------y )  
-!               (           ---z           ---z  )
-!           1   (            U              V    )
-!     WC = ---  ( W - DZX * ---    - DZY * ---   )
-!          DZZ  (           DXX            DYY   )
-!
-!  
-!       In the no-topography case, WC = W / DZZ
-!
-!
-!!**  METHOD
-!!    ------
-!!      We employ the Shuman operators to compute the averages. The metric
-!!    coefficients PDXX, PDYY, PDZX, PDZY, PDZZ are dummy arguments
-!!
-!!
-!!    EXTERNAL 
-!!    --------
-!!      MXF, MYF, MZM         : Shuman functions (mean operators)
-!!       
-!!      Module MODI_SHUMAN    : Interface for Shuman functions   
-!!
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-!!
-!!      Module MODD_CONF   : contains configuration variable 
-!!           LFLAT : Logical for topography
-!!                  = .TRUE.  if Zs = 0 (Flat terrain)
-!!
-!!    REFERENCE
-!!    ---------
-!!      Book2 of documentation (subroutine CONTRAV)
-!!
-!!
-!!    AUTHOR
-!!    ------
-!!      J.L. Redelsperger     * CNRM *
-!!	J.-P. Pinty      * Laboratoire d'Aerologie*
-!!
-!!    MODIFICATIONS
-!!    -------------
-!!      Original   27/07/94
-!!      Corrections 3/08/94 (by J.P. Lafore)
-!!      Corrections 17/10/94 (by J.P. Lafore) WC modified for w-advection
-!----------------------------------------------------------------------------
-!
-!*       0.    DECLARATIONS
-!              ------------
-USE MODD_CONF
-USE MODD_PARAMETERS
-!
-!
-USE MODI_SHUMAN
-!
-IMPLICIT NONE
-!
-!*       0.1   declarations of arguments    
-!
-INTEGER                                ::  IACC
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRUT     ! Cartesian comp along x
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRVT     ! Cartesian comp along y
-REAL, DIMENSION(:,:,:),  INTENT(IN)    ::  PRWT     ! Cartesian comp along z
-!
-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(OUT)   ::  PRUCT    ! Contrav comp along x-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRVCT    ! Contrav comp along y-bar
-REAL, DIMENSION(:,:,:),  INTENT(OUT)   ::  PRWCT    ! Contrav comp along z-bar
-!$acc declare pcopyin (PRUT,PRVT,PRWT,pdxx,pdyy,pdzz,pdzx,pdzy)
-!$acc declare pcopyout(PRUCT,PRVCT,PRWCT)
+IF (KADV_ORDER == 4 ) THEN
+ IF( .NOT. LFLAT) THEN 
+  NULLIFY(TZFIELD_U)
+  NULLIFY(TZFIELD_V)
+  CALL ADD3DFIELD_ll(TZFIELD_U, PRUCT)
+  CALL ADD3DFIELD_ll(TZFIELD_V, PRVCT)
+  CALL UPDATE_HALO_ll(TZFIELD_U,IINFO_ll)
+  CALL UPDATE_HALO_ll(TZFIELD_V,IINFO_ll)
+!!$ IF( NHALO==1 ) THEN 
+  NULLIFY(TZFIELD_DZX)
+  NULLIFY(TZFIELD_DZY)
+  CALL ADD3DFIELD_ll(TZFIELD_DZX, PDZX)
+  CALL ADD3DFIELD_ll(TZFIELD_DZY, PDZY)
+  NULLIFY(TZHALO2_U)
+  NULLIFY(TZHALO2_V)
+  NULLIFY(TZHALO2_DZX)
+  NULLIFY(TZHALO2_DZY)
+  CALL INIT_HALO2_ll(TZHALO2_U,1,IIU,IJU,IKU)
+  CALL INIT_HALO2_ll(TZHALO2_V,1,IIU,IJU,IKU)
+  CALL INIT_HALO2_ll(TZHALO2_DZX,1,IIU,IJU,IKU)
+  CALL INIT_HALO2_ll(TZHALO2_DZY,1,IIU,IJU,IKU)
+  CALL UPDATE_HALO2_ll(TZFIELD_U, TZHALO2_U, IINFO_ll)
+  CALL UPDATE_HALO2_ll(TZFIELD_V, TZHALO2_V, IINFO_ll)
+  CALL UPDATE_HALO2_ll(TZFIELD_DZX, TZHALO2_DZX, IINFO_ll)
+  CALL UPDATE_HALO2_ll(TZFIELD_DZY, TZHALO2_DZY, IINFO_ll)
+!!$ END IF
+ END IF
+END IF
 !
 !
-!*       0.2   declarations of local variables
-!              
-REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,2),SIZE(PDXX,3)):: Z1,Z2
-INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE
+!*       2.    Compute the vertical contravariant components (flat case)
+!              ------------------------------------
 !
-!-----------------------------------------------------------------------
+IF (LFLAT) THEN
+  PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)
+  RETURN
+END IF
 !
-!*       1.    Compute the contravariant components
+!*       3.    Compute the vertical contravariant components (general case)
 !              ------------------------------------
 !
-IIB=2
-IJB=2
-IIE=SIZE(PDXX,1)-1
-IJE=SIZE(PDXX,2)-1
-!
-IKB=1+JPVEXT
-IKE=SIZE(PDXX,3) - JPVEXT
+Z1 = 0.
+Z2 = 0.
 !
-!$acc data create (z1,z2 )
-!$acc kernels 
+IF (KADV_ORDER == 2 ) THEN
 !
-PRUCT(:,:,:) = PRUT(:,:,:) / PDXX(:,:,:)
-PRVCT(:,:,:) = PRVT(:,:,:) / PDYY(:,:,:)
-
-IF (LFLAT) THEN
-
-  PRWCT(:,:,:) = PRWT(:,:,:) / PDZZ(:,:,:)  
-  
-ELSE 
-
   Z1(IIB:IIE,:,IKB:IKE+1)=                                             &
       (PRUCT(IIB:IIE,:,IKB:IKE+1)+PRUCT(IIB:IIE,:,IKB-1:IKE) )         &
        *PDZX(IIB:IIE,:,IKB:IKE+1) *0.25                                &
      +(PRUCT(IIB+1:IIE+1,:,IKB:IKE+1)+PRUCT(IIB+1:IIE+1,:,IKB-1:IKE) ) &
        *PDZX(IIB+1:IIE+1,:,IKB:IKE+1) *0.25   
-                      
+                        
   Z2(:,IJB:IJE,IKB:IKE+1)=                                             &
       (PRVCT(:,IJB:IJE,IKB:IKE+1)+PRVCT(:,IJB:IJE,IKB-1:IKE) )         &
        *PDZY(:,IJB:IJE,IKB:IKE+1) *0.25                                &
      +(PRVCT(:,IJB+1:IJE+1,IKB:IKE+1)+PRVCT(:,IJB+1:IJE+1,IKB-1:IKE) ) &
        *PDZY(:,IJB+1:IJE+1,IKB:IKE+1) *0.25   
-
   PRWCT=0.             
   PRWCT(IIB:IIE,IJB:IJE,IKB:IKE+1) =                 &
       (   PRWT(IIB:IIE,IJB:IJE,IKB:IKE+1)            &
         - Z1(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
         - Z2(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
       ) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE+1)  
- 
+!
+ELSE IF (KADV_ORDER == 4 ) THEN
+!
+!!$   IF (NHALO == 1) THEN
+      IF ( LWEST_ll() .AND. HLBCX(1)/='CYCL' ) THEN
+         IW=IIB+2 -1
+      ELSE
+         IW=IIB+1 -1
+      END IF
+      IE=IIE-1
+!!$   ELSE
+!!$      IF (LWEST_ll()) THEN
+!!$         IW=IIB+1
+!!$      ELSE
+!!$         IW=IIB
+!!$      END IF
+!!$      IF (LEAST_ll() .AND. HLBCX(2)/='CYCL' ) THEN
+!!$         IE=IIE-1
+!!$      ELSE
+!!$         IE=IIE
+!!$      END IF
+!!$   END IF
+   !
+!!$   IF(NHALO == 1) THEN
+      IF ( LSOUTH_ll() .AND. HLBCY(1)/='CYCL' ) THEN
+         IS=IJB+2 -1
+      ELSE
+         IS=IJB+1 -1
+      END IF
+      IN=IJE-1
+!!$   ELSE
+!!$      IF (LSOUTH_ll()) THEN
+!!$         IS=IJB+1
+!!$      ELSE
+!!$         IS=IJB
+!!$      END IF
+!!$      IF (LNORTH_ll() .AND. HLBCY(2)/='CYCL' ) THEN
+!!$         IN=IJE-1
+!!$      ELSE
+!!$         IN=IJE
+!!$      END IF
+!!$   END IF
+   !
+   !
+   !*       3.1    interior of the processor subdomain
+!
+!
+   Z1(IW:IE,:,IKB:IKE+1)=                                               &
+       7.0*( (PRUCT(IW:IE,:,IKB:IKE+1)+PRUCT(IW:IE,:,IKB-1:IKE)) &
+            *( 9.0*PDZX(IW:IE,:,IKB:IKE+1)-(PDZX(IW+1:IE+1,:,IKB:IKE+1) &
+                  +PDZX(IW:IE,:,IKB:IKE+1)+PDZX(IW-1:IE-1,:,IKB:IKE+1))/3.0)/8.0 * 0.5 &
+           +(PRUCT(IW+1:IE+1,:,IKB:IKE+1)+PRUCT(IW+1:IE+1,:,IKB-1:IKE)) &
+            *( 9.0*PDZX(IW+1:IE+1,:,IKB:IKE+1)-(PDZX(IW+2:IE+2,:,IKB:IKE+1) &
+                  +PDZX(IW+1:IE+1,:,IKB:IKE+1)+PDZX(IW:IE,:,IKB:IKE+1))/3.0)/8.0 * 0.5 )/12.0 &
+          -( (PRUCT(IW-1:IE-1,:,IKB:IKE+1)+PRUCT(IW-1:IE-1,:,IKB-1:IKE)) &
+            *PDZX(IW-1:IE-1,:,IKB:IKE+1) *0.5 &
+            +(PRUCT(IW+2:IE+2,:,IKB:IKE+1)+PRUCT(IW+2:IE+2,:,IKB-1:IKE)) &
+            *PDZX(IW+2:IE+2,:,IKB:IKE+1) *0.5)/12.0
+
+!
+   Z2(:,IS:IN,IKB:IKE+1)=                                               &
+       7.0*( (PRVCT(:,IS:IN,IKB:IKE+1)+PRVCT(:,IS:IN,IKB-1:IKE)) &
+            *( 9.0*PDZY(:,IS:IN,IKB:IKE+1)-(PDZY(:,IS+1:IN+1,IKB:IKE+1) &
+                  +PDZY(:,IS:IN,IKB:IKE+1)+PDZY(:,IS-1:IN-1,IKB:IKE+1))/3.0)/8.0 * 0.5 &
+           +(PRVCT(:,IS+1:IN+1,IKB:IKE+1)+PRVCT(:,IS+1:IN+1,IKB-1:IKE)) &
+            *( 9.0*PDZY(:,IS+1:IN+1,IKB:IKE+1)-(PDZY(:,IS+2:IN+2,IKB:IKE+1) &
+                  +PDZY(:,IS+1:IN+1,IKB:IKE+1)+PDZY(:,IS:IN,IKB:IKE+1))/3.0)/8.0 * 0.5 )/12.0 &
+          -( (PRVCT(:,IS-1:IN-1,IKB:IKE+1)+PRVCT(:,IS-1:IN-1,IKB-1:IKE)) &
+            *PDZY(:,IS-1:IN-1,IKB:IKE+1) *0.5 &
+            +(PRVCT(:,IS+2:IN+2,IKB:IKE+1)+PRVCT(:,IS+2:IN+2,IKB-1:IKE)) &
+            *PDZY(:,IS+2:IN+2,IKB:IKE+1) *0.5)/12.0
+!
+!*       3.2    limits of the processor subdomain (inside the whole domain or in cyclic conditions)
+!
+!!$  IF (NHALO==1) THEN
+
+   Z1(IIE,:,IKB:IKE+1)=                                               &
+       7.0*( (PRUCT(IIE,:,IKB:IKE+1)+PRUCT(IIE,:,IKB-1:IKE)) &
+            *( 9.0*PDZX(IIE,:,IKB:IKE+1)-(PDZX(IIE+1,:,IKB:IKE+1) &
+                  +PDZX(IIE,:,IKB:IKE+1)+PDZX(IIE-1,:,IKB:IKE+1))/3.0)/8.0 * 0.5 &
+           +(PRUCT(IIE+1,:,IKB:IKE+1)+PRUCT(IIE+1,:,IKB-1:IKE)) &
+            *( 9.0*PDZX(IIE+1,:,IKB:IKE+1)-(TZHALO2_DZX%HALO2%EAST(:,IKB:IKE+1) &
+                  +PDZX(IIE+1,:,IKB:IKE+1)+PDZX(IIE,:,IKB:IKE+1))/3.0)/8.0 * 0.5 )/12.0 &
+          -( (PRUCT(IIE-1,:,IKB:IKE+1)+PRUCT(IIE-1,:,IKB-1:IKE)) &
+            *PDZX(IIE-1,:,IKB:IKE+1) *0.5 &
+            +(TZHALO2_U%HALO2%EAST(:,IKB:IKE+1)+TZHALO2_U%HALO2%EAST(:,IKB-1:IKE)) &
+            *TZHALO2_DZX%HALO2%EAST(:,IKB:IKE+1) *0.5)/12.0
+!
+   Z2(:,IJE,IKB:IKE+1)=                                               &
+       7.0*( (PRVCT(:,IJE,IKB:IKE+1)+PRVCT(:,IJE,IKB-1:IKE)) &
+            *( 9.0*PDZY(:,IJE,IKB:IKE+1)-(PDZY(:,IJE+1,IKB:IKE+1) &
+                  +PDZY(:,IJE,IKB:IKE+1)+PDZY(:,IJE-1,IKB:IKE+1))/3.0)/8.0 * 0.5 &
+           +(PRVCT(:,IJE+1,IKB:IKE+1)+PRVCT(:,IJE+1,IKB-1:IKE)) &
+            *( 9.0*PDZY(:,IJE+1,IKB:IKE+1)-(TZHALO2_DZY%HALO2%NORTH(:,IKB:IKE+1) &
+                  +PDZY(:,IJE+1,IKB:IKE+1)+PDZY(:,IJE,IKB:IKE+1))/3.0)/8.0 * 0.5 )/12.0 &
+          -( (PRVCT(:,IJE-1,IKB:IKE+1)+PRVCT(:,IJE-1,IKB-1:IKE)) &
+            *PDZY(:,IJE-1,IKB:IKE+1) *0.5 &
+            +(TZHALO2_V%HALO2%NORTH(:,IKB:IKE+1)+TZHALO2_V%HALO2%NORTH(:,IKB-1:IKE)) &
+            *TZHALO2_DZY%HALO2%NORTH(:,IKB:IKE+1) *0.5)/12.0
+!!$  END IF
+!
+!*       3.3    non-CYCLIC CASE IN THE X DIRECTION: 2nd order case
+!
+  IF (HLBCX(1)/='CYCL' .AND. LWEST_ll()) THEN
+!
+   Z1(IIB,:,IKB:IKE+1)=                                     &
+       (PRUCT(IIB,:,IKB:IKE+1)+PRUCT(IIB,:,IKB-1:IKE) )     &
+        *PDZX(IIB,:,IKB:IKE+1) *0.25                        &
+      +(PRUCT(IIB+1,:,IKB:IKE+1)+PRUCT(IIB+1,:,IKB-1:IKE) ) &
+        *PDZX(IIB+1,:,IKB:IKE+1) *0.25   
+  END IF
+!
+  IF (HLBCX(2)/='CYCL' .AND. LEAST_ll()) THEN
+!
+   Z1(IIE,:,IKB:IKE+1)=                                     &
+       (PRUCT(IIE,:,IKB:IKE+1)+PRUCT(IIE,:,IKB-1:IKE) )     &
+        *PDZX(IIE,:,IKB:IKE+1) *0.25                        &
+      +(PRUCT(IIE+1,:,IKB:IKE+1)+PRUCT(IIE+1,:,IKB-1:IKE) ) &
+        *PDZX(IIE+1,:,IKB:IKE+1) *0.25   
+  END IF
+!
+!*       3.4    non-CYCLIC CASE IN THE Y DIRECTION: 2nd order case
+!
+  IF (HLBCY(1)/='CYCL' .AND. LSOUTH_ll()) THEN
+!
+   Z2(:,IJB,IKB:IKE+1)=                                     &
+       (PRVCT(:,IJB,IKB:IKE+1)+PRVCT(:,IJB,IKB-1:IKE) )     &
+        *PDZY(:,IJB,IKB:IKE+1) *0.25                        &
+      +(PRVCT(:,IJB+1,IKB:IKE+1)+PRVCT(:,IJB+1,IKB-1:IKE) ) &
+        *PDZY(:,IJB+1,IKB:IKE+1) *0.25   
+!
+  END IF
+!
+  IF (HLBCY(2)/='CYCL' .AND. LNORTH_ll()) THEN
+!
+   Z2(:,IJE,IKB:IKE+1)=                                     &
+       (PRVCT(:,IJE,IKB:IKE+1)+PRVCT(:,IJE,IKB-1:IKE) )     &
+        *PDZY(:,IJE,IKB:IKE+1) *0.25                        &
+      +(PRVCT(:,IJE+1,IKB:IKE+1)+PRVCT(:,IJE+1,IKB-1:IKE) ) &
+        *PDZY(:,IJE+1,IKB:IKE+1) *0.25   
+!
+  END IF
+!
+!*       3.5    Vertical contyravariant wind
+!
+!
+!!$  CALL GET_HALO(Z1)
+!!$  CALL GET_HALO(Z2)
+!!$
+!!$  CALL MPPDB_CHECK3DM("contrav ::Z1/Z2/ PDZZ",PRECISION,Z1,Z2,PDZZ)                    
+  PRWCT=0.             
+  PRWCT(IIB:IIE,IJB:IJE,IKB:IKE+1) =                 &
+     (   PRWT(IIB:IIE,IJB:IJE,IKB:IKE+1)            &
+       - Z1(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
+       - Z2(IIB:IIE,IJB:IJE,IKB:IKE+1)              &
+     ) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE+1)  
+!
 END IF
 !
 PRWCT(:,:,1) = - PRWCT(:,:,3)     ! Mirror hypothesis
-!$acc end kernels
-!$acc end data
 !
+IF (KADV_ORDER == 4 ) THEN
+  CALL CLEANLIST_ll(TZFIELD_U)
+  CALL CLEANLIST_ll(TZFIELD_V)
+!!$  IF (NHALO==1) THEN
+    CALL CLEANLIST_ll(TZFIELD_DZX)
+    CALL CLEANLIST_ll(TZFIELD_DZY)
+    CALL DEL_HALO2_ll(TZHALO2_U)
+    CALL DEL_HALO2_ll(TZHALO2_V)
+    CALL DEL_HALO2_ll(TZHALO2_DZX)
+    CALL DEL_HALO2_ll(TZHALO2_DZY)
+!!$  END IF
+END IF
 !-----------------------------------------------------------------------
+CALL MPPDB_CHECK3DM("contrav end ::PRU/V/WCT",PRECISION,PRUCT,PRVCT,PRWCT)                    
 !
-END SUBROUTINE CONTRAV_ACC
+END SUBROUTINE CONTRAV
diff --git a/MNH/deallocate_model1.f90 b/MNH/deallocate_model1.f90
index 9872cb98d57513d6a3b722bf170b864c7c29deb7..cc4c251761f1900ce72d4def15ca6d6811c6a98b 100644
--- a/MNH/deallocate_model1.f90
+++ b/MNH/deallocate_model1.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/deallocate_model1.f90,v $ $Revision: 1.2.2.2.2.2.2.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/deallocate_model1.f90,v $ $Revision: 1.2.2.2.2.2.2.1.10.2.2.4 $
 ! MASDEV4_7 spawn 2007/01/12 14:12:39
 !-----------------------------------------------------------------
 !############################
@@ -64,6 +68,7 @@ END MODULE MODI_DEALLOCATE_MODEL1
 !!                   15/03/99 new PGD fields
 !!                   08/03/01 D.Gazen add chemical emission field
 !!                   01/2004 V. Masson surface externalization
+!!                   06/2012 M.Tomasini add 2D nesting ADVFRC
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -88,6 +93,16 @@ USE MODD_RAIN_ICE_PARAM
 USE MODD_RAIN_ICE_DESCR
 USE MODE_MODELN_HANDLER
 !
+! Modif 2D
+USE MODD_LATZ_EDFLX                      ! For ADVFRC and EDDY FLUXES
+USE MODD_DEF_EDDY_FLUX_n           ! For EDDY FLUXES
+USE MODD_DEF_EDDYUV_FLUX_n         ! For EDDY FLUXES
+!
+USE MODD_2D_FRC
+USE MODD_ADVFRC_n                  ! For ADVFRC and EDDY FLUXES
+USE MODD_RELFRC_n
+USE MODD_ADV_n
+USE MODD_PAST_FIELD_n
 IMPLICIT NONE
 !
 !*       0.1   declarations of arguments
@@ -105,25 +120,47 @@ CALL GOTO_MODEL(1)
 !*       1.    Module MODD_FIELD$n
 !
 IF ( KCALL==3 ) THEN
-  DEALLOCATE(XUM)
-  DEALLOCATE(XVM)
-  DEALLOCATE(XWM)
-  DEALLOCATE(XTHM)
-  IF (ASSOCIATED(XUT)) DEALLOCATE(XUT)
-  IF (ASSOCIATED(XVT)) DEALLOCATE(XVT)
-  IF (ASSOCIATED(XWT)) DEALLOCATE(XWT)
-  !!$IF (ASSOCIATED(XTHT)) DEALLOCATE(XTHT)
-  IF (ALLOCATED(XTHT)) DEALLOCATE(XTHT)
+  IF (CUVW_ADV_SCHEME(1:3)=='CEN') THEN
+    DEALLOCATE(XUM)
+    DEALLOCATE(XVM)
+    DEALLOCATE(XWM)
+    DEALLOCATE(XDUM)
+    DEALLOCATE(XDVM)
+    DEALLOCATE(XDWM)
+  END IF
+  DEALLOCATE(XUT)
+  DEALLOCATE(XVT)
+  DEALLOCATE(XWT)
+  DEALLOCATE(XTHT)
+  IF (L2D_ADV_FRC) THEN
+    IF (ASSOCIATED(XDTHFRC)) DEALLOCATE(XDTHFRC)
+    IF (ASSOCIATED(XDRVFRC)) DEALLOCATE(XDRVFRC)
+    IF (ASSOCIATED(TDTADVFRC)) DEALLOCATE(TDTADVFRC)
+  END IF
+    IF (L2D_REL_FRC) THEN
+    IF (ASSOCIATED(XTHREL)) DEALLOCATE(XTHREL)
+    IF (ASSOCIATED(XRVREL)) DEALLOCATE(XRVREL)
+    IF (ASSOCIATED(TDTRELFRC)) DEALLOCATE(TDTRELFRC)
+  END IF
+  ! DEALLOCATE EDDY FLUXES 
+  IF (LTH_FLX) THEN
+    DEALLOCATE(XVTH_FLUX_M)
+    DEALLOCATE(XWTH_FLUX_M)
+  END IF
+  IF (LUV_FLX) THEN
+    DEALLOCATE(XVU_FLUX_M)
+  END IF
 END IF
 IF ( KCALL==1 ) THEN
   DEALLOCATE(XRUS)
   DEALLOCATE(XRVS)
   DEALLOCATE(XRWS)
   DEALLOCATE(XRTHS)
+  DEALLOCATE(XRUS_PRES, XRVS_PRES, XRWS_PRES )                
+  DEALLOCATE(XRTHS_CLD )
 END IF
 !
 IF ( KCALL==3 ) THEN
-  IF (ASSOCIATED(XTKEM)) DEALLOCATE(XTKEM)
   IF (ASSOCIATED(XTKET)) DEALLOCATE(XTKET)
 END IF
 IF ( ASSOCIATED(XRTKES) .AND. KCALL==1 ) THEN
@@ -131,21 +168,18 @@ IF ( ASSOCIATED(XRTKES) .AND. KCALL==1 ) THEN
 END IF
 !
 IF ( KCALL==3 ) THEN
-  DEALLOCATE(XPABSM)
-  IF (ASSOCIATED(XPABST)) DEALLOCATE(XPABST)
+  DEALLOCATE(XPABST)
 !
-  DEALLOCATE(XRM)
-!!$  IF (ASSOCIATED(XRT)) DEALLOCATE(XRT)
-  IF (ALLOCATED(XRT)) DEALLOCATE(XRT)
+  DEALLOCATE(XRT)
 END IF
 !
 IF ( KCALL==1 ) THEN
   DEALLOCATE(XRRS)
+  DEALLOCATE(XRRS_CLD)
 END IF
 !
-IF ( ASSOCIATED(XSRCM) .AND. KCALL==3 ) THEN
-  DEALLOCATE(XSRCM)
-  IF (ASSOCIATED(XSRCT)) DEALLOCATE(XSRCT)
+IF ( ASSOCIATED(XSRCT) .AND. KCALL==3 ) THEN
+  DEALLOCATE(XSRCT)
   DEALLOCATE(XSIGS)
 END IF   
 !
@@ -154,11 +188,11 @@ IF ( ASSOCIATED(XCLDFR) .AND. KCALL==2 ) THEN
 END IF   
 !
 IF ( KCALL == 3 ) THEN
-  DEALLOCATE(XSVM)
-  IF (ASSOCIATED(XSVT)) DEALLOCATE(XSVT)
+  DEALLOCATE(XSVT)
 END IF
 IF ( KCALL == 1 ) THEN
   DEALLOCATE(XRSVS)
+  DEALLOCATE(XRSVS_CLD)
 END IF
 !
 !
@@ -242,6 +276,8 @@ IF ( KCALL == 1 ) THEN
   DEALLOCATE(XRHOM)
   DEALLOCATE(XALK)
   DEALLOCATE(XALKW)
+  DEALLOCATE(XALKBAS)
+  DEALLOCATE(XALKWBAS)
   IF ( ASSOCIATED(XKURELAX) ) THEN
     DEALLOCATE(XKURELAX)
     DEALLOCATE(XKVRELAX)
@@ -439,6 +475,11 @@ IF ( ASSOCIATED(XINPRS) .AND. KCALL == 3 ) THEN
   DEALLOCATE(XACPRG)
 END IF
 !
+IF ( ASSOCIATED(XINPRH) .AND. KCALL == 3 ) THEN
+  DEALLOCATE(XINPRH)
+  DEALLOCATE(XACPRH)
+END IF
+!
 !*     13b.     Module MODD_ELEC$n
 !
 IF ( ASSOCIATED(XNI_SDRYG) .AND. KCALL == 3 ) THEN
diff --git a/MNH/get_halo.f90 b/MNH/get_halo.f90
index c3ff9459884d8ff1197a7ca0aef1d20e8e4219c3..427622c6e151a0c0fadb5eddd69363b90d1d17f2 100644
--- a/MNH/get_halo.f90
+++ b/MNH/get_halo.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/get_halo.f90,v $ $Revision: 1.1.2.1.2.2 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/get_halo.f90,v $ $Revision: 1.1.2.1.2.2.18.2 $
 ! MASDEV4_7 newsrc 2007/03/01 13:18:33
 !-----------------------------------------------------------------
 !     ####################
diff --git a/MNH/gradient_m.f90 b/MNH/gradient_m.f90
index 59c0b435b1f12b49783eecdda2e206a27266cec9..81eb7e850786f528d35c896f94941fcdb5f1c0af 100644
--- a/MNH/gradient_m.f90
+++ b/MNH/gradient_m.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/gradient_m.f90,v $ $Revision: 1.1.10.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/gradient_m.f90,v $ $Revision: 1.1.10.1.16.1.2.2 $
 ! MASDEV4_7 operators 2006/05/18 13:07:25
 !-----------------------------------------------------------------
 !     ######################
@@ -11,7 +15,9 @@
 INTERFACE
 !
 !
-FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX)      RESULT(PGX_M_M)
+FUNCTION GX_M_M(KKA,KKU,KL,PA,PDXX,PDZZ,PDZX)      RESULT(PGX_M_M)
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDXX    ! metric coefficient dxx
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
@@ -22,8 +28,9 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
 END FUNCTION GX_M_M
 !
 !
-FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY)      RESULT(PGY_M_M)
-!
+FUNCTION GY_M_M(KKA,KKU,KL,PA,PDYY,PDZZ,PDZY)      RESULT(PGY_M_M)
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDYY    ! metric coefficient dyy
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
@@ -34,8 +41,10 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
 END FUNCTION GY_M_M
 !
 !
-FUNCTION GZ_M_M(PA,PDZZ)      RESULT(PGZ_M_M)
+FUNCTION GZ_M_M(KKA,KKU,KL,PA,PDZZ)      RESULT(PGZ_M_M)
 !
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
 !
@@ -43,11 +52,12 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
 !
 END FUNCTION GZ_M_M
 !
-      FUNCTION GX_M_U(PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+      FUNCTION GX_M_U(KKA,KKU,KL,PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
 !  
 IMPLICIT NONE
 !
-                                                          ! Metric coefficients:
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDXX                   ! d*xx
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZX                   ! d*zx 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   ! d*zz
@@ -59,11 +69,12 @@ REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGX_M_U  ! result at flux
 END FUNCTION GX_M_U
 !
 !
-      FUNCTION GY_M_V(PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+      FUNCTION GY_M_V(KKA,KKU,KL,PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
 !
 IMPLICIT NONE
 !
-                                                         ! Metric coefficients:
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDYY                   !d*yy
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZY                   !d*zy 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   !d*zz
@@ -74,11 +85,13 @@ REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGY_M_V  ! result at flux
                                                               ! side
 END FUNCTION GY_M_V
 !
-      FUNCTION GZ_M_W(PY,PDZZ) RESULT(PGZ_M_W)
+      FUNCTION GZ_M_W(KKA,KKU,KL,PY,PDZZ) RESULT(PGZ_M_W)
 !  
 IMPLICIT NONE
 !
                                                           ! Metric coefficient:
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   !d*zz
 !
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PY       ! variable at mass
@@ -95,7 +108,7 @@ END MODULE MODI_GRADIENT_M
 !
 !
 !     #######################################################
-      FUNCTION GX_M_M(PA,PDXX,PDZZ,PDZX)      RESULT(PGX_M_M)
+      FUNCTION GX_M_M(KKA,KKU,KL,PA,PDXX,PDZZ,PDZX)      RESULT(PGX_M_M)
 !     #######################################################
 !
 !!****  *GX_M_M* - Cartesian Gradient operator: 
@@ -160,6 +173,8 @@ IMPLICIT NONE
 !
 !*       0.1   declarations of arguments and result
 !
+INTEGER,                 INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,                 INTENT(IN)  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDXX    ! metric coefficient dxx
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
@@ -179,8 +194,8 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGX_M_M ! result mass point
 !
 IF (.NOT. LFLAT) THEN 
   PGX_M_M(:,:,:)= (DXF(MXM(PA(:,:,:)))   -                     &
-                   MZF(MXF(PDZX)*DZM(PA(:,:,:))/PDZZ(:,:,:)) ) &
-                  /MXF(PDXX(:,:,:))
+                   MZF(KKA,KKU,KL,MXF(PDZX)*DZM(KKA,KKU,KL,PA(:,:,:)) &
+                  /PDZZ(:,:,:)) )        /MXF(PDXX(:,:,:))
 ELSE
   PGX_M_M(:,:,:)=DXF(MXM(PA(:,:,:))) / MXF(PDXX(:,:,:)) 
 END IF
@@ -191,7 +206,7 @@ END FUNCTION GX_M_M
 !
 !
 !     #######################################################
-      FUNCTION GY_M_M(PA,PDYY,PDZZ,PDZY)      RESULT(PGY_M_M)
+      FUNCTION GY_M_M(KKA,KKU,KL,PA,PDYY,PDZZ,PDZY)      RESULT(PGY_M_M)
 !     #######################################################
 !
 !!****  *GY_M_M* - Cartesian Gradient operator: 
@@ -255,6 +270,8 @@ IMPLICIT NONE
 !
 !*       0.1   declarations of arguments and result
 !
+INTEGER,                 INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,                 INTENT(IN)  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDYY    ! metric coefficient dyy
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
@@ -274,8 +291,8 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGY_M_M ! result mass point
 !
 !
 IF (.NOT. LFLAT) THEN 
-  PGY_M_M(:,:,:)= (DYF(MYM(PA)) - MZF(MYF(PDZY)*DZM(PA)/PDZZ) ) &
-                /MYF(PDYY)
+  PGY_M_M(:,:,:)= (DYF(MYM(PA))-MZF(KKA,KKU,KL,MYF(PDZY)*DZM(KKA,KKU,KL,PA)&
+                /PDZZ) ) /MYF(PDYY)
 ELSE
   PGY_M_M(:,:,:)= DYF(MYM(PA))/MYF(PDYY)
 ENDIF  
@@ -288,7 +305,7 @@ END FUNCTION GY_M_M
 !
 !
 !     #######################################################
-      FUNCTION GZ_M_M(PA,PDZZ)      RESULT(PGZ_M_M)
+      FUNCTION GZ_M_M(KKA,KKU,KL,PA,PDZZ)      RESULT(PGZ_M_M)
 !     #######################################################
 !
 !!****  *GZ_M_M* - Cartesian Gradient operator: 
@@ -347,6 +364,8 @@ IMPLICIT NONE
 !
 !*       0.1   declarations of arguments and result
 !
+INTEGER,              INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PA      ! variable at the mass point
 REAL, DIMENSION(:,:,:),  INTENT(IN)  :: PDZZ    ! metric coefficient dzz
 !
@@ -362,7 +381,7 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PGZ_M_M ! result mass point
 !*       1.    DEFINITION of GZ_M_M
 !              --------------------
 !
-PGZ_M_M(:,:,:)= MZF( DZM(PA(:,:,:))/PDZZ(:,:,:) )
+PGZ_M_M(:,:,:)= MZF(KKA,KKU,KL, DZM(KKA,KKU,KL,PA(:,:,:))/PDZZ(:,:,:) )
 !
 !----------------------------------------------------------------------------
 !
@@ -370,7 +389,7 @@ END FUNCTION GZ_M_M
 !
 !
 !     ##################################################
-      FUNCTION GX_M_U(PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
+      FUNCTION GX_M_U(KKA,KKU,KL,PY,PDXX,PDZZ,PDZX) RESULT(PGX_M_U)
 !     ##################################################
 !
 !!****  *GX_M_U * - Compute the gradient along x for a variable localized at 
@@ -440,7 +459,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of arguments and result
 !              ------------------------------------
 !
-                                                          ! Metric coefficients:
+INTEGER,                INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,                INTENT(IN)  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDXX                   ! d*xx
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZX                   ! d*zx 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   ! d*zz
@@ -482,6 +502,7 @@ IF (.NOT. LFLAT) THEN
               ) * PDZX(JI,:,JK+1)* 0.25                                   &
             )  / PDXX(JI,:,JK)
     END DO
+
   END DO
 !$acc end region
 
@@ -511,8 +532,8 @@ IF (.NOT. LFLAT) THEN
 !
 !$acc region
   DO JI=1+JPHEXT,IIU
-    PGX_M_U(JI,:,IKU)=  ( PY(JI,:,IKU)-PY(JI-1,:,IKU)  )  / PDXX(JI,:,IKU) 
-    PGX_M_U(JI,:,1)=  -999.
+    PGX_M_U(JI,:,KKU)=  ( PY(JI,:,KKU)-PY(JI-1,:,KKU)  )  / PDXX(JI,:,KKU) 
+    PGX_M_U(JI,:,KKA)=  -999.
   END DO
 !
   PGX_M_U(1,:,:)=PGX_M_U(IIU-2*JPHEXT+1,:,:)
@@ -533,7 +554,7 @@ END FUNCTION GX_M_U
 !
 !
 !     ##################################################
-      FUNCTION GY_M_V(PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
+      FUNCTION GY_M_V(KKA,KKU,KL,PY,PDYY,PDZZ,PDZY) RESULT(PGY_M_V)
 !     ##################################################
 !
 !!****  *GY_M_V * - Compute the gradient along y for a variable localized at 
@@ -603,7 +624,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of arguments and results
 !              -------------------------------------
 !
-                                                         ! Metric coefficients:
+INTEGER,                INTENT(IN)  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,                INTENT(IN)  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDYY                   !d*yy
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZY                   !d*zy 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   !d*zz
@@ -623,24 +645,25 @@ IJU=SIZE(PY,2)
 IKU=SIZE(PY,3)
 IF (.NOT. LFLAT) THEN 
 !  PGY_M_V = (   DYM(PY)  -  MZF (   MYM(  DZM(PY) /PDZZ  ) * PDZY   )   )/PDYY
+
 !$acc region
-  DO JK=1+JPVEXT,IKU-JPVEXT
+  DO JK=1+JPVEXT_TURB,IKU-JPVEXT_TURB
     DO JJ=1+JPHEXT,IJU
         PGY_M_V(:,JJ,JK)=                                                 &
            (  PY(:,JJ,JK)-PY(:,JJ-1,JK)                                   &
-             -(  (PY(:,JJ,JK)-PY(:,JJ,JK-1))     / PDZZ(:,JJ,JK)          &
-                +(PY(:,JJ-1,JK)-PY(:,JJ-1,JK-1)) / PDZZ(:,JJ-1,JK)        &
+             -(  (PY(:,JJ,JK)-PY(:,JJ,JK-KL))     / PDZZ(:,JJ,JK)          &
+                +(PY(:,JJ-1,JK)-PY(:,JJ-KL,JK-KL)) / PDZZ(:,JJ-1,JK)        &
               ) * PDZY(:,JJ,JK)* 0.25                                     &
-             -(  (PY(:,JJ,JK+1)-PY(:,JJ,JK))     / PDZZ(:,JJ,JK+1)        &
-                +(PY(:,JJ-1,JK+1)-PY(:,JJ-1,JK)) / PDZZ(:,JJ-1,JK+1)      &
-              ) * PDZY(:,JJ,JK+1)* 0.25                                   &
+             -(  (PY(:,JJ,JK+KL)-PY(:,JJ,JK))     / PDZZ(:,JJ,JK+KL)        &
+                +(PY(:,JJ-1,JK+KL)-PY(:,JJ-1,JK)) / PDZZ(:,JJ-1,JK+KL)      &
+              ) * PDZY(:,JJ,JK+KL)* 0.25                                   &
             )  / PDYY(:,JJ,JK)
     END DO
   END DO
 !
   DO JJ=1+JPHEXT,IJU
-    PGY_M_V(:,JJ,IKU)=  ( PY(:,JJ,IKU)-PY(:,JJ-1,IKU)  )  / PDYY(:,JJ,IKU) 
-    PGY_M_V(:,JJ,1)=  -999.
+    PGY_M_V(:,JJ,KKU)=  ( PY(:,JJ,KKU)-PY(:,JJ-1,KKU)  )  / PDYY(:,JJ,KKU) 
+    PGY_M_V(:,JJ,KKA)=  -999.
   END DO
 !
   PGY_M_V(:,1,:)=PGY_M_V(:,IJU-2*JPHEXT+1,:)
@@ -659,7 +682,7 @@ END FUNCTION GY_M_V
 !
 !
 !     #########################################
-      FUNCTION GZ_M_W(PY,PDZZ) RESULT(PGZ_M_W)
+      FUNCTION GZ_M_W(KKA,KKU,KL,PY,PDZZ) RESULT(PGZ_M_W)
 !     #########################################
 !
 !!****  *GZ_M_W * - Compute the gradient along z direction for a 
@@ -718,6 +741,9 @@ IMPLICIT NONE
 !*       0.1   Declarations of arguments and results
 !              -------------------------------------
 !
+INTEGER,           INTENT(IN)     :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,           INTENT(IN)     :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+
                                                           ! Metric coefficient:
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDZZ                   !d*zz
 !
@@ -726,16 +752,21 @@ REAL, DIMENSION(:,:,:), INTENT(IN)                :: PY       ! variable at mass
 REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PGZ_M_W  ! result at flux
                                                               ! side
 !
-INTEGER  IKU
+INTEGER :: IKT,IKTB,IKTE
 !-------------------------------------------------------------------------------
 !
 !*       1.    COMPUTE THE GRADIENT ALONG Z
 !              -----------------------------
 !
-IKU=SIZE(PY,3)
-PGZ_M_W(:,:,1+JPVEXT:IKU) =  (PY(:,:,1+JPVEXT:IKU)-PY(:,:,JPVEXT:IKU-1))  &
-                           / PDZZ(:,:,1+JPVEXT:IKU)
-PGZ_M_W(:,:,1)=-999.
+IKT=SIZE(PY,3)
+IKTB=1+JPVEXT_TURB              
+IKTE=IKT-JPVEXT_TURB
+
+PGZ_M_W(:,:,IKTB:IKTE) =  (PY(:,:,IKTB:IKTE)-PY(:,:,IKTB-KL:IKTE-KL))  &
+                           / PDZZ(:,:,IKTB:IKTE)
+PGZ_M_W(:,:,KKU)=  (PY(:,:,KKU)-PY(:,:,KKU-KL))  &
+                           / PDZZ(:,:,KKU)
+PGZ_M_W(:,:,KKA)=-999.
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/MNH/init_mnh.f90 b/MNH/init_mnh.f90
index 80766f248c379944f03398a75b4221f3c9fe5351..ac3ee5281c0611aecda5a523bc189113027b3048 100644
--- a/MNH/init_mnh.f90
+++ b/MNH/init_mnh.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/init_mnh.f90,v $ $Revision: 1.1.2.1.2.1.2.4 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/init_mnh.f90,v $ $Revision: 1.1.2.1.2.1.2.4.2.3.2.3 $
 ! masdev4_7 BUG1 2007/06/15 17:47:18
 !-----------------------------------------------------------------
 !     ###############
@@ -85,19 +89,18 @@ USE MODE_FM
 !
 USE MODI_INI_CST
 USE MODI_INI_CTURB
-USE MODI_INI_CMFSHALL
 USE MODI_INI_SEG_n
 USE MODI_INI_MODEL_n
 USE MODI_INI_SIZE_n
 USE MODI_INI_SIZE_SPAWN
 USE MODI_RESET_EXSEG
 USE MODE_MODELN_HANDLER
+USE MODI_READ_ALL_NAMELISTS
+USE MODI_ALLOC_SURFEX
 !
 !JUAN
 USE MODE_SPLITTINGZ_ll
 !JUAN
-USE MODI_ALLOC_SURFEX
-USE MODI_READ_ALL_NAMELISTS
 !
 USE MODE_MNH_ZWORK
 !
@@ -111,11 +114,11 @@ CHARACTER (LEN=16), DIMENSION(JPMODELMAX) :: YLUOUT   ! Name for output-listing
                                                       ! of nested models
 CHARACTER (LEN=28), DIMENSION(JPMODELMAX) :: YINIFILE ! names of
                                                       ! the initial files
+CHARACTER (LEN=28), DIMENSION(JPMODELMAX) :: YINIFILEPGD                                                     
 INTEGER  :: ILUOUT0,IRESP                             ! Logical unit number for
                                                       ! output-listing common
                                                       ! to all models and return
                                                       ! code of file management
-REAL, DIMENSION(JPMODELMAX)            :: ZTSTEP_OLD  ! OLD Time STEP (DESFM)
 REAL, DIMENSION(JPMODELMAX)            :: ZTSTEP_ALL  ! Time STEP of ALL models
 INTEGER                                :: IINFO_ll    ! return code of // routines
 !
@@ -155,9 +158,9 @@ CALL INI_CST
 CALL INI_CTURB                      
 !
 !
-!*       1.4    initialize constants for the mass flux scheme                            
+!*       1.4    initialize constants for nebulosity computation                            
 !
-CALL INI_CMFSHALL
+CALL INI_NEB
 !
 !-------------------------------------------------------------------------------
 !
@@ -166,7 +169,7 @@ CALL INI_CMFSHALL
 !
 DO JMI=1,JPMODELMAX
   CALL GOTO_MODEL(JMI)
-  CALL INI_SEG_n(JMI,YLUOUT(JMI),YINIFILE(JMI),ZTSTEP_OLD(JMI),ZTSTEP_ALL)
+  CALL INI_SEG_n(JMI,YLUOUT(JMI),YINIFILE(JMI),YINIFILEPGD(JMI),ZTSTEP_ALL)
   IF (JMI.EQ.NMODEL) EXIT
 END DO
 !
@@ -180,7 +183,7 @@ IF (CPROGRAM=='DIAG') CALL RESET_EXSEG(YLUOUT(1))
 !
 DO JMI=1,NMODEL
   CALL GOTO_MODEL(JMI)
-  CALL INI_SIZE_n(JMI,YLUOUT(JMI),YINIFILE(JMI))
+  CALL INI_SIZE_n(JMI,YLUOUT(JMI),YINIFILE(JMI),YINIFILEPGD(JMI))
 END DO
 !
 IF (CPROGRAM=='SPAWN ') THEN 
@@ -215,7 +218,7 @@ CALL MNH_ALLOC_ZWORK(NMODEL)
 DO JMI=1,NMODEL
   CALL GO_TOMODEL_ll(JMI,IINFO_ll)
   CALL GOTO_MODEL(JMI)
-  CALL INI_MODEL_n(JMI,ZTSTEP_OLD(JMI),YLUOUT(JMI),YINIFILE(JMI))
+  CALL INI_MODEL_n(JMI,YLUOUT(JMI),YINIFILE(JMI),YINIFILEPGD(JMI))
 END DO
 !
 !-------------------------------------------------------------------------------
diff --git a/MNH/mesonh.f90 b/MNH/mesonh.f90
index f5b6c3ecc074893c46131595a0144a38a2f75645..a332d0771dd9333590ba0a4bded7ca7fb0307223 100644
--- a/MNH/mesonh.f90
+++ b/MNH/mesonh.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/mesonh.f90,v $ $Revision: 1.2.2.3.2.2.2.3 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/mesonh.f90,v $ $Revision: 1.2.2.3.2.2.2.3.4.4 $
 ! masdev4_7 BUG1 2007/06/22 10:53:31
 !-----------------------------------------------------------------
 !     ##############
@@ -111,6 +115,9 @@ INTEGER       :: IINFO_ll                     ! return code of // routines
 CALL MPPDB_INIT()
 !
 ! Switch to model 1 variables
+IF (LCHECK) THEN
+ CALL MPPDB_INIT()
+END IF
 CALL GOTO_MODEL(1)
 !
 CALL INITIO_ll()
@@ -131,7 +138,7 @@ GEXIT=.FALSE.
 DO JMODEL=1,NMODEL
   CALL GO_TOMODEL_ll(JMODEL,IINFO_ll)
   CALL GOTO_MODEL(JMODEL)
-  CSTORAGE_TYPE='MT'
+  CSTORAGE_TYPE='TT'
   CALL MODEL_n(1,GEXIT)
 END DO
 !
@@ -160,7 +167,12 @@ END DO
 !*       3.    FINALIZE THE PARALLEL SESSION
 !              -----------------------------
 !
-CALL END_PARA_ll(IINFO_ll)
+IF (LCHECK) THEN
+  CALL MPPDB_BARRIER()
+  CALL MPPDB_BARRIER()
+ELSE
+  CALL END_PARA_ll(IINFO_ll)
+END IF
 !
 !
 CALL DEALLOC_SURFEX
diff --git a/MNH/metrics.f90 b/MNH/metrics.f90
index b7685197bacefb8b73db3fb62be07126c9197b45..15ede79f530a63d4636e26aec00206f97b9f02e3 100644
--- a/MNH/metrics.f90
+++ b/MNH/metrics.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/metrics.f90,v $ $Revision: 1.1.8.1.2.1 $ $Date: 2009/04/21 07:42:51 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/metrics.f90,v $ $Revision: 1.1.8.1.2.1.16.1.2.2 $ $Date: 2014/01/09 15:01:56 $
 !-----------------------------------------------------------------
 !-----------------------------------------------------------------
 !-----------------------------------------------------------------
@@ -134,7 +138,7 @@ ELSE
   ZD1=1.
 END IF
 IF (.NOT.LCARTESIAN) THEN
-  ZDZZ(:,:,:) = MZF( 1.+ ZD1*PZZ(:,:,:)/XRADIUS)
+  ZDZZ(:,:,:) = MZF(1,IKU,1, 1.+ ZD1*PZZ(:,:,:)/XRADIUS)
   DO JK=1,IKU ; DO JJ=1,IJU ; DO JI=1,IIU
     PDXX(JI,JJ,JK) = ZDZZ(JI,JJ,JK) * PDXHAT(JI) /PMAP(JI,JJ)
     PDYY(JI,JJ,JK) = ZDZZ(JI,JJ,JK) * PDYHAT(JJ) /PMAP(JI,JJ)
@@ -166,7 +170,7 @@ PDZY(:,:,:) = DYM(PZZ(:,:,:))
 !*       4.  COMPUTE PDZZ  :
 !            -------------
 !
-PDZZ(:,:,:) = DZM(MZF(PZZ(:,:,:)))
+PDZZ(:,:,:) = DZM(1,IKU,1,MZF(1,IKU,1,PZZ(:,:,:)))
 PDZZ(:,:,IKU) = PZZ(:,:,IKU) - PZZ(:,:,IKU-1)  ! same delta z in IKU and IKU -1
 PDZZ(:,:,1)   = PDZZ(:,:,2)                    ! same delta z in 1   and 2
 !-----------------------------------------------------------------------------
diff --git a/MNH/modd_fieldn.f90 b/MNH/modd_fieldn.f90
index cddbd9efbaa541a92557d674d1d3889131c3631d..9de20f084de0d7c7234a9e4a6fbbb4f8fd3bb760 100644
--- a/MNH/modd_fieldn.f90
+++ b/MNH/modd_fieldn.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_fieldn.f90,v $ $Revision: 1.2.4.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_fieldn.f90,v $ $Revision: 1.2.4.1.18.3 $
 ! MASDEV4_7 modd 2006/06/27 14:17:24
 !-----------------------------------------------------------------
 !     ###################
@@ -43,6 +47,11 @@
 !!                     11/04/96  (J.-P. Pinty) add the ice concentration
 !!                     25/07/97  (J. Stein) Change the variable pressure
 !!                     20/05/06  Remove EPS
+!!                        11/11  (C.Lac) FIT version : Remove t-Dt fields except for 
+!!                                       radiative cooling (microphysics) +
+!!                               add pressure contribution to the tendencies for
+!!                               momentum (noted _PRES) + microphysics contrib
+!!                               for Theta and r (noted _CLD)
 !!
 !-------------------------------------------------------------------------------
 !
@@ -53,86 +62,80 @@ USE MODD_PARAMETERS, ONLY: JPMODELMAX
 IMPLICIT NONE
 
 TYPE FIELD_t
-  REAL, DIMENSION(:,:,:), POINTER :: XUM=>NULL(),XVM=>NULL(),XWM=>NULL()
-                                      ! U,V,W  at time t-dt
   REAL, DIMENSION(:,:,:), POINTER :: XUT=>NULL(),XVT=>NULL(),XWT=>NULL()
                                       ! U,V,W  at time t
   REAL, DIMENSION(:,:,:), POINTER :: XRUS=>NULL(),XRVS=>NULL(),XRWS=>NULL()
                                       ! Source of (rho U), (rho V), (rho w) 
-  REAL, DIMENSION(:,:,:), POINTER :: XTHM=>NULL()     ! (rho theta) at time t-dt
+  REAL, DIMENSION(:,:,:), POINTER :: XRUS_PRES=>NULL(),XRVS_PRES=>NULL(),XRWS_PRES=>NULL()
   REAL, DIMENSION(:,:,:), POINTER :: XTHT=>NULL()     ! (rho theta) at time t
   REAL, DIMENSION(:,:,:), POINTER :: XRTHS=>NULL()    ! Source of (rho theta)
-  REAL, DIMENSION(:,:,:), POINTER :: XTKEM=>NULL()    ! Kinetic energy 
-                                                     ! at time t-dt
+  REAL, DIMENSION(:,:,:), POINTER :: XRTHS_CLD=>NULL()    ! Source of (rho theta) from resolved_cloud
   REAL, DIMENSION(:,:,:), POINTER :: XTKET=>NULL()    ! Kinetic energy
                                                      ! at time t
   REAL, DIMENSION(:,:,:), POINTER :: XRTKES=>NULL()   ! Source of kinetic energy
                                                      ! (rho e)
-  REAL, DIMENSION(:,:,:), POINTER :: XPABSM=>NULL()   ! absolute pressure at
-                                                     ! time t-dt
   REAL, DIMENSION(:,:,:), POINTER :: XPABST=>NULL()   ! absolute pressure at
                                                      ! time t
-  REAL, DIMENSION(:,:,:,:), POINTER :: XRM=>NULL()    ! Moist variables 
-                                                     ! at time t-dt
   REAL, DIMENSION(:,:,:,:), POINTER :: XRT=>NULL()    ! Moist variables (rho Rn) 
                                                      ! at time t
   REAL, DIMENSION(:,:,:,:), POINTER :: XRRS=>NULL()   ! Source of Moist variables
                                                      ! (rho Rn) 
-  REAL, DIMENSION(:,:,:,:), POINTER :: XSVM=>NULL()   ! Additionnal scalar
-                                                     ! variables at time t-dt
+  REAL, DIMENSION(:,:,:,:), POINTER :: XRRS_CLD=>NULL()   ! Source of Moist variables
   REAL, DIMENSION(:,:,:,:), POINTER :: XSVT=>NULL()   ! Additionnal scalar
                                                      ! variables at time t  
   REAL, DIMENSION(:,:,:,:), POINTER :: XRSVS=>NULL()  ! Source of addi. scalar
                                                      !  variables (rho Sn.) 
+  REAL, DIMENSION(:,:,:,:), POINTER :: XRSVS_CLD=>NULL() ! Source of (rho Sn) from resolved_cloud
   REAL                          ::   XDRYMASST    ! Mass of dry air Md
   REAL                          ::   XDRYMASSS    ! LS sources of Md
   REAL, DIMENSION(:,:,:), POINTER :: XSRC=>NULL()     ! turbulent flux <s'Rc'>
   REAL, DIMENSION(:,:,:), POINTER :: XSIGS=>NULL()    ! =sqrt(<s's'>) for the
                                                      ! Subgrid Condensation
   REAL, DIMENSION(:,:,:), POINTER :: XCLDFR=>NULL()   ! cloud fraction
-  REAL, DIMENSION(:,:,:), POINTER :: XSRCM=>NULL()    ! turbulent flux <s'Rc'>
-                                                     ! at t- delta t
   REAL, DIMENSION(:,:,:), POINTER :: XSRCT=>NULL()    ! turbulent flux <s'Rc'>
                                                      ! at t
   REAL, DIMENSION(:,:,:), POINTER :: XCIT=>NULL()     ! Pristine ice concentration
+  REAL, DIMENSION(:,:,:), POINTER :: XTHM=>NULL()      ! Theta at Previous time step 
+  REAL, DIMENSION(:,:,:), POINTER :: XRCM=>NULL()      ! Cloud mixing ratio at Previous time step
+  REAL, DIMENSION(:,:,:), POINTER :: XPABSM=>NULL()      ! Theta at Previous time step 
 !
 END TYPE FIELD_t
 
 TYPE(FIELD_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: FIELD_MODEL
 
-REAL, DIMENSION(:,:,:), POINTER :: XUM=>NULL(),XVM=>NULL(),XWM=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XUT=>NULL(),XVT=>NULL(),XWT=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XRUS=>NULL(),XRVS=>NULL(),XRWS=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XTHM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRUS_PRES=>NULL(),XRVS_PRES=>NULL(),XRWS_PRES=>NULL()
 !!$REAL, DIMENSION(:,:,:), POINTER :: XTHT=>NULL()
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: XTHT
+REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET  :: XTHT
 !$acc declare mirror (XTHT)
 !!$REAL, DIMENSION(:,:,:), POINTER :: XRTHS=>NULL()
 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: XRTHS
 !$acc declare mirror (XRTHS)
-REAL, DIMENSION(:,:,:), POINTER :: XTKEM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRTHS_CLD=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XTKET=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XRTKES=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XPABSM=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XPABST=>NULL()
-REAL, DIMENSION(:,:,:,:), POINTER :: XRM=>NULL()
 !!$REAL, DIMENSION(:,:,:,:), POINTER :: XRT=>NULL()
-REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XRT
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: XRT
 !$acc declare mirror (XRT)
 !!$REAL, DIMENSION(:,:,:,:), POINTER :: XRRS=>NULL()
 REAL, DIMENSION(:,:,:,:), ALLOCATABLE, TARGET :: XRRS
 !$acc declare mirror (XRRS)
-REAL, DIMENSION(:,:,:,:), POINTER :: XSVM=>NULL()
+REAL, DIMENSION(:,:,:,:), POINTER :: XRRS_CLD=>NULL()
 REAL, DIMENSION(:,:,:,:), POINTER :: XSVT=>NULL()
 REAL, DIMENSION(:,:,:,:), POINTER :: XRSVS=>NULL()
+REAL, DIMENSION(:,:,:,:), POINTER :: XRSVS_CLD=>NULL()
 REAL, POINTER :: XDRYMASST=>NULL()
 REAL, POINTER :: XDRYMASSS=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XSRC=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XSIGS=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XCLDFR=>NULL()
-REAL, DIMENSION(:,:,:), POINTER :: XSRCM=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XSRCT=>NULL()
 REAL, DIMENSION(:,:,:), POINTER :: XCIT=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XTHM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XPABSM=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRCM=>NULL()
 
 CONTAINS
 
@@ -140,68 +143,68 @@ SUBROUTINE FIELD_GOTO_MODEL(KFROM, KTO)
 INTEGER, INTENT(IN) :: KFROM, KTO
 !
 ! Save current state for allocated arrays
-FIELD_MODEL(KFROM)%XUM=>XUM
-FIELD_MODEL(KFROM)%XVM=>XVM
-FIELD_MODEL(KFROM)%XWM=>XWM
 FIELD_MODEL(KFROM)%XUT=>XUT
 FIELD_MODEL(KFROM)%XVT=>XVT
 FIELD_MODEL(KFROM)%XWT=>XWT
 FIELD_MODEL(KFROM)%XRUS=>XRUS
 FIELD_MODEL(KFROM)%XRVS=>XRVS
 FIELD_MODEL(KFROM)%XRWS=>XRWS
-FIELD_MODEL(KFROM)%XTHM=>XTHM
+FIELD_MODEL(KFROM)%XRUS_PRES=>XRUS_PRES
+FIELD_MODEL(KFROM)%XRVS_PRES=>XRVS_PRES
+FIELD_MODEL(KFROM)%XRWS_PRES=>XRWS_PRES
 !!$FIELD_MODEL(KFROM)%XTHT=>XTHT
 !!$FIELD_MODEL(KFROM)%XRTHS=>XRTHS
-FIELD_MODEL(KFROM)%XTKEM=>XTKEM
+FIELD_MODEL(KFROM)%XRTHS_CLD=>XRTHS_CLD
 FIELD_MODEL(KFROM)%XTKET=>XTKET
 FIELD_MODEL(KFROM)%XRTKES=>XRTKES
-FIELD_MODEL(KFROM)%XPABSM=>XPABSM
 FIELD_MODEL(KFROM)%XPABST=>XPABST
-FIELD_MODEL(KFROM)%XRM=>XRM
 !!$FIELD_MODEL(KFROM)%XRT=>XRT
 FIELD_MODEL(KFROM)%XRRS=>XRRS
-FIELD_MODEL(KFROM)%XSVM=>XSVM
+FIELD_MODEL(KFROM)%XRRS_CLD=>XRRS_CLD
 FIELD_MODEL(KFROM)%XSVT=>XSVT
 FIELD_MODEL(KFROM)%XRSVS=>XRSVS
+FIELD_MODEL(KFROM)%XRSVS_CLD=>XRSVS_CLD
 FIELD_MODEL(KFROM)%XSRC=>XSRC
 FIELD_MODEL(KFROM)%XSIGS=>XSIGS
 FIELD_MODEL(KFROM)%XCLDFR=>XCLDFR
-FIELD_MODEL(KFROM)%XSRCM=>XSRCM
 FIELD_MODEL(KFROM)%XSRCT=>XSRCT
 FIELD_MODEL(KFROM)%XCIT=>XCIT
+FIELD_MODEL(KFROM)%XTHM=>XTHM
+FIELD_MODEL(KFROM)%XPABSM=>XPABSM
+FIELD_MODEL(KFROM)%XRCM=>XRCM
 !
 ! Current model is set to model KTO
-XUM=>FIELD_MODEL(KTO)%XUM
-XVM=>FIELD_MODEL(KTO)%XVM
-XWM=>FIELD_MODEL(KTO)%XWM
 XUT=>FIELD_MODEL(KTO)%XUT
 XVT=>FIELD_MODEL(KTO)%XVT
 XWT=>FIELD_MODEL(KTO)%XWT
 XRUS=>FIELD_MODEL(KTO)%XRUS
 XRVS=>FIELD_MODEL(KTO)%XRVS
 XRWS=>FIELD_MODEL(KTO)%XRWS
-XTHM=>FIELD_MODEL(KTO)%XTHM
+XRUS_PRES=>FIELD_MODEL(KTO)%XRUS_PRES
+XRVS_PRES=>FIELD_MODEL(KTO)%XRVS_PRES
+XRWS_PRES=>FIELD_MODEL(KTO)%XRWS_PRES
 !!$XTHT=>FIELD_MODEL(KTO)%XTHT
 !!$XRTHS=>FIELD_MODEL(KTO)%XRTHS
-XTKEM=>FIELD_MODEL(KTO)%XTKEM
+XRTHS_CLD=>FIELD_MODEL(KTO)%XRTHS_CLD
 XTKET=>FIELD_MODEL(KTO)%XTKET
 XRTKES=>FIELD_MODEL(KTO)%XRTKES
-XPABSM=>FIELD_MODEL(KTO)%XPABSM
 XPABST=>FIELD_MODEL(KTO)%XPABST
-XRM=>FIELD_MODEL(KTO)%XRM
 !!$XRT=>FIELD_MODEL(KTO)%XRT
 !!$XRRS=>FIELD_MODEL(KTO)%XRRS
-XSVM=>FIELD_MODEL(KTO)%XSVM
+XRRS_CLD=>FIELD_MODEL(KTO)%XRRS_CLD
 XSVT=>FIELD_MODEL(KTO)%XSVT
 XRSVS=>FIELD_MODEL(KTO)%XRSVS
+XRSVS_CLD=>FIELD_MODEL(KTO)%XRSVS_CLD
 XDRYMASST=>FIELD_MODEL(KTO)%XDRYMASST
 XDRYMASSS=>FIELD_MODEL(KTO)%XDRYMASSS
 XSRC=>FIELD_MODEL(KTO)%XSRC
 XSIGS=>FIELD_MODEL(KTO)%XSIGS
 XCLDFR=>FIELD_MODEL(KTO)%XCLDFR
-XSRCM=>FIELD_MODEL(KTO)%XSRCM
 XSRCT=>FIELD_MODEL(KTO)%XSRCT
 XCIT=>FIELD_MODEL(KTO)%XCIT
+XTHM=>FIELD_MODEL(KTO)%XTHM
+XPABSM=>FIELD_MODEL(KTO)%XPABSM
+XRCM=>FIELD_MODEL(KTO)%XRCM
 
 END SUBROUTINE FIELD_GOTO_MODEL
 
diff --git a/MNH/modd_metricsn.f90 b/MNH/modd_metricsn.f90
index d89188ded953f4642d7d1611a0575794de4e1f7d..293f8502afb2dcdfaa05743713a8cc7a2ed4da1e 100644
--- a/MNH/modd_metricsn.f90
+++ b/MNH/modd_metricsn.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_metricsn.f90,v $ $Revision: 1.2.4.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_metricsn.f90,v $ $Revision: 1.2.4.1.18.2 $
 ! MASDEV4_7 modd 2006/06/27 14:20:29
 !-----------------------------------------------------------------
 !     #####################
diff --git a/MNH/modd_refn.f90 b/MNH/modd_refn.f90
index 5527f5b4a4f1171fc1f009df003708e818a17a1a..f7a90afb5dc92bada59f21fc9209306449185670 100644
--- a/MNH/modd_refn.f90
+++ b/MNH/modd_refn.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_refn.f90,v $ $Revision: 1.2.4.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_refn.f90,v $ $Revision: 1.2.4.1.18.2 $
 ! MASDEV4_7 modd 2006/06/27 12:51:57
 !-----------------------------------------------------------------
 !     #################
diff --git a/MNH/ppm.f90 b/MNH/ppm.f90
index 2f678f7269f98e2aa554b0019bd54ee3c771cc6b..b1900ba100ccf66a403f14168260a2153f04e257 100644
--- a/MNH/ppm.f90
+++ b/MNH/ppm.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/ppm.f90,v $ $Revision: 1.1.2.3.2.2 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/Attic/ppm.f90,v $ $Revision: 1.1.2.3.2.3.16.1.2.3 $
 ! masdev4_7 BUG1 2007/06/15 17:47:18
 !-----------------------------------------------------------------
 !     ###############
@@ -236,6 +240,7 @@ CONTAINS
 !!    -------------
 !!    
 !!    11.5.2006.  T. Maric - original version
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !!
 !-------------------------------------------------------------------------------
 !
@@ -321,15 +326,16 @@ GEAST = LEAST_ll()
 !
 !*              initialise & update halo for PSRC
 !
-IF(NHALO /= 1) THEN
-    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
-    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
-    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
-   !callabortstop
-    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
-    CALL ABORT
-    STOP
-ENDIF
+
+!!$IF(NHALO /= 1) THEN
+!!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
+!!$    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
+!!$    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
+!!$   !callabortstop
+!!$    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+!!$    CALL ABORT
+!!$    STOP
+!!$ENDIF
 !
 CALL GET_HALO_D(PSRC,HDIR="01_X")
 !
@@ -806,6 +812,7 @@ CONTAINS
 !!    -------------
 !!    
 !!    11.5.2006.  T. Maric - original version
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !!
 !-------------------------------------------------------------------------------
 !
@@ -881,15 +888,15 @@ IIA=IIE
 GSOUTH=LSOUTH_ll()
 GNORTH=LNORTH_ll()
 
-IF(NHALO /= 1) THEN
-    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
-    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
-    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
-   !callabortstop
-    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
-    CALL ABORT
-    STOP
-ENDIF
+!!$IF(NHALO /= 1) THEN
+!!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
+!!$    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
+!!$    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
+!!$   !callabortstop
+!!$    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+!!$    CALL ABORT
+!!$    STOP
+!!$ENDIF
 !
 CALL GET_HALO_D(PSRC,HDIR="01_Y")
 !
@@ -1421,6 +1428,7 @@ REAL, DIMENSION(IIU,IJU,IKU) :: &
 !
 INTEGER:: IKB    ! Begining useful area in x,y,z directions
 INTEGER:: IKE    ! End useful area in x,y,z directions
+INTEGER:: IKU
 !
 !JUAN ACC
 INTEGER                          :: IIU,IJU,IKU
@@ -1442,6 +1450,7 @@ DQ(:,:,IKE+1) = -DQ(:,:,IKE) ! DIF2Z(DQ,PQ)
 !
 IKB = 1 + JPVEXT
 IKE = SIZE(PSRC,3) - JPVEXT
+
 !!$
 !!$IIU=size(psrc,1)
 !!$IJU=size(psrc,2)
@@ -1453,6 +1462,7 @@ IKE = SIZE(PSRC,3) - JPVEXT
 ! create/mirror (ZDMQ,ZQL0,ZQR0,ZDQ,ZQ60,ZQL,ZQR,ZQ6,ZFPOS,ZFNEG) 
 !$acc kernels  
 #endif
+
 !
 !-------------------------------------------------------------------------------
 !
@@ -1531,6 +1541,7 @@ ZDQ = ZQR - ZQL
 ZFPOS(:,:,IKB+1:IKE+1) = ZQR(:,:,IKB:IKE) - 0.5*PCR(:,:,IKB+1:IKE+1) * &            
      (ZDQ(:,:,IKB:IKE) - (1.0 - 2.0*PCR(:,:,IKB+1:IKE+1)/3.0)        &
      * ZQ6(:,:,IKB:IKE))
+
 !
 ! advection flux at open boundary when u(IKB) > 0
 ZFPOS(:,:,IKB) = (PSRC(:,:,IKB-1) - ZQR(:,:,IKB-1))*PCR(:,:,IKB) + &
@@ -1688,6 +1699,7 @@ CONTAINS
 !!    -------------
 !!    
 !!    20.6.2006.  T. Maric - original version
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !!
 !-------------------------------------------------------------------------------
 !
@@ -1780,15 +1792,15 @@ REAL, DIMENSION(SIZE(PCR,2),SIZE(PCR,3))             :: ZPSRC_HALO2_WEST
 !
 !*              initialise & update halo & halo2 for PSRC
 !
-IF(NHALO /= 1) THEN
-    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
-    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
-    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
-   !callabortstop
-    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
-    CALL ABORT
-    STOP
-ENDIF
+!!$IF(NHALO /= 1) THEN
+!!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
+!!$    WRITE(ILUOUT,*) 'ERROR : PPM ppm_met.f90 --> Juan '
+!!$    WRITE(ILUOUT,*) 'PPM not yet implemented/tested with NHALO /= 1'
+!!$   !callabortstop
+!!$    CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+!!$    CALL ABORT
+!!$    STOP
+!!$ENDIF
 !
 IF ( .NOT. GSMONOPROC ) THEN
 CALL GET_HALO2(PSRC,TZ_PSRC_HALO2_ll)
@@ -2427,6 +2439,7 @@ REAL, DIMENSION(:,:,:),INTENT(INOUT):: PR &
 !
 !*       0.2   Declarations of local variables :
 !
+
 ! advection fluxes
   &                                                   , ZFPOS, ZFNEG &
 !
@@ -3101,6 +3114,7 @@ REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: PR
 !
 INTEGER:: IIB,IJB,IKB   ! Begining useful area in x,y,z directions
 INTEGER:: IIE,IJE,IKE   ! End useful area in x,y,z directions
+INTEGER:: IKU
 !
 ! variable at cell edges
 REAL, DIMENSION(SIZE(PCR,1),SIZE(PCR,2),SIZE(PCR,3)) :: ZPHAT, ZRVT
@@ -3128,10 +3142,11 @@ INTEGER :: II, IJ, IK
 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
 IKB = 1 + JPVEXT
 IKE = SIZE(PSRC,3) - JPVEXT
+IKU =  SIZE(PSRC,3)
 !
 !-------------------------------------------------------------------------------
 !
-ZRVT = PCR/PTSTEP * MZM(PRHO)
+ZRVT = PCR/PTSTEP * MZM(1,IKU,1,PRHO)
 !
 ! calculate 4th order fluxes at cell edges in the inner domain !
 ZPHAT(:,:,IKB+1:IKE) = (7.0 * &
@@ -3238,7 +3253,7 @@ END WHERE
 !
 ! 1. calculate upwind tendency of the source
 !
-PR = PSRC*PRHO - PTSTEP*DZF(ZFUP)
+PR = PSRC*PRHO - PTSTEP*DZF(1,IKU,1,ZFUP)
 !
 !-------------------------------------------------------------------------------
 ! compute and apply the limiters
@@ -3366,7 +3381,7 @@ ZFCOR(:,:,IKB-1) = MIN( &
 !-------------------------------------------------------------------------------
 ! 6. apply the limited flux correction to scalar field
 !
-PR = PR - PTSTEP*DZF(ZFCOR)
+PR = PR - PTSTEP*DZF(1,IKU,1,ZFCOR)
 !
 !
 END FUNCTION PPM_S1_Z
diff --git a/MNH/ppm_met.f90 b/MNH/ppm_met.f90
index 78667fdb0bfdf30f0bef4be574d64d272ea0e748..282c3103d9928a907c80fb9568dfffc4012a832e 100644
--- a/MNH/ppm_met.f90
+++ b/MNH/ppm_met.f90
@@ -1,15 +1,22 @@
-MODULE MODI_PPM_MET
+!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.
+!
+!     #####################
+      MODULE MODI_PPM_MET  
+!     #####################
 !
 CONTAINS
 !
 !     ######################################################################
       SUBROUTINE PPM_MET (HLBCX,HLBCY, KRR, KTCOUNT,              &
                           PCRU, PCRV, PCRW, PTSTEP, PRHODJ,       &
-                          PTHT, PTKET, PRT,                       &
+                          PRHOX1, PRHOX2, PRHOY1, PRHOY2,         &
+                          PRHOZ1, PRHOZ2, PTHT, PTKET, PRT,       &
                           PRTHS, PRTKES, PRRS, HMET_ADV_SCHEME    )
         
         USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
-        USE MODE_MNH_ZWORK, ONLY : ZUNIT3D
 
         IMPLICIT NONE
 !
@@ -22,45 +29,44 @@ CHARACTER (LEN=6),               INTENT(IN) :: HMET_ADV_SCHEME
 INTEGER,                  INTENT(IN)    :: KRR    ! Number of moist variables
 INTEGER,                  INTENT(IN)    :: KTCOUNT! iteration count
 !
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! contravariant
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  !  components
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! of momentum
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! Courant
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  ! numbers
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! 
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ ! density
+! Temporary advected rhodj
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOX1,PRHOX2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOY1,PRHOY2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOZ1,PRHOZ2
 !
-REAL,                     INTENT(IN)    :: PTSTEP ! Time step 
+REAL,                     INTENT(IN)    :: PTSTEP ! Single Time step 
 !
-REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT, PTKET ! Vars at t
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT, PTKET        ! Vars at t
 REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT 
 !
 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS, PRTKES! Source terms
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS 
 
+INTEGER :: IZSRC
 
-INTEGER :: IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2
-
-      CALL  MNH_GET_ZT3D(IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2)
+      CALL  MNH_GET_ZT3D(IZSRC)
 
       CALL PPM_MET_D(HLBCX,HLBCY, KRR, KTCOUNT,              &
                    & PCRU, PCRV, PCRW, PTSTEP, PRHODJ,       &
-                   & PTHT, PTKET, PRT,                       &
+                   & PRHOX1, PRHOX2, PRHOY1, PRHOY2,         &
+                   & PRHOZ1, PRHOZ2, PTHT, PTKET, PRT,       &
                    & PRTHS, PRTKES, PRRS, HMET_ADV_SCHEME ,  &
-                   & ZT3D(:,:,:,IZSRC), &
-                   & ZT3D(:,:,:,IZRHOX1),ZT3D(:,:,:,IZRHOX2),ZT3D(:,:,:,IZRHOY1), &
-                   & ZT3D(:,:,:,IZRHOY2),ZT3D(:,:,:,IZRHOZ1),ZT3D(:,:,:,IZRHOZ2), &
-                   & ZUNIT3D   )
+                   & ZT3D(:,:,:,IZSRC)  )
 
-      CALL  MNH_REL_ZT3D(IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2)
+      CALL  MNH_REL_ZT3D(IZSRC)
 
 CONTAINS
 
       SUBROUTINE PPM_MET_D(HLBCX,HLBCY, KRR, KTCOUNT,                 &
                          & PCRU, PCRV, PCRW, PTSTEP, PRHODJ,          &
-                         & PTHT, PTKET, PRT,                          &
+                         & PRHOX1, PRHOX2, PRHOY1, PRHOY2,            &
+                         & PRHOZ1, PRHOZ2, PTHT, PTKET, PRT,          &
                          & PRTHS, PRTKES, PRRS, HMET_ADV_SCHEME,      &
-                         & ZSRC,                                      &
-                         & ZRHOX1,ZRHOX2,ZRHOY1,ZRHOY2,ZRHOZ1,ZRHOZ2, &
-                         & ZUNIT   )
-
+                         & ZSRC  )
 
 !     ######################################################################
 !
@@ -90,6 +96,7 @@ CONTAINS
 !!    MODIFICATIONS
 !!    -------------
 !!      Original 11.05.2006. T.Maric
+!!      Modification : 11.2011 C.Lac, V.Masson : Advection of (theta_l,r_t) 
 !!
 !-------------------------------------------------------------------------------
 !
@@ -100,17 +107,15 @@ CONTAINS
 !
 USE MODD_PARAMETERS
 USE MODD_CONF
-USE MODD_BUDGET
-USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 !
 USE MODI_SHUMAN
-USE MODI_BUDGET
 USE MODI_PPM
 USE MODI_ADVEC_PPM_ALGO
 !
 ! incorporate ADVEC_4TH_ORDER_ALG, MZF4 and MZM4
 !USE MODI_ADVEC_4TH_ORDER_AUX
 !
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -126,6 +131,10 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! contravariant
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  !  components
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! of momentum
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ ! density
+! Temporary advected rhodj
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOX1,PRHOX2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOY1,PRHOY2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOZ1,PRHOZ2
 !
 REAL,                     INTENT(IN)    :: PTSTEP ! Time step 
 !
@@ -135,7 +144,6 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PRT
 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS, PRTKES! Source terms
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS 
 !
-!
 !*       0.2   Declarations of local variables :
 !
 INTEGER :: JRR           ! Loop index for  moist variables
@@ -150,13 +158,6 @@ INTEGER :: IGRID ! localisation on the model grid
 REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZSRC
 !$acc declare present (ZSRC)
 !
-! Temporary advected rhodj
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOX1,ZRHOX2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOY1,ZRHOY2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOZ1,ZRHOZ2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZUNIT
-!$acc declare present (ZRHOX1,ZRHOX2,ZRHOY1,ZRHOY2,ZRHOZ1,ZRHOZ2,ZUNIT)
-!
 !-------------------------------------------------------------------------------
 !
 !*       1.     COMPUTES THE DOMAIN DIMENSIONS
@@ -171,21 +172,12 @@ GTKEALLOC = SIZE(PTKET,1) /= 0
 !
 IGRID = 1
 !
-! Calculate the advection of the density RHODJ to pass to the algorithm
-!
-!$acc data pcopyin(PCRU,PCRV,PCRW,PRHODJ)
-CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHODJ, PTSTEP, ZRHOX1)
-CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, ZRHOX1, PTSTEP, ZRHOY1)
-CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, ZRHOY1, PTSTEP, ZRHOZ1)
-CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHODJ, PTSTEP, ZRHOZ2)
-CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, ZRHOZ2, PTSTEP, ZRHOY2)
-CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, ZRHOY2, PTSTEP, ZRHOX2)
 !
 ! Potential temperature
 !
 !$acc data pcopyin (PTHT) pcopy(PRTHS)
 CALL ADVEC_PPM_ALGO(HMET_ADV_SCHEME, HLBCX, HLBCY, IGRID, PTHT, PRHODJ, PTSTEP, &
-                    ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+                    PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1, PRHOZ2, &
                     ZSRC, KTCOUNT, PCRU, PCRV, PCRW)
 ! add the advection to the sources
 !$acc kernels
@@ -193,47 +185,39 @@ PRTHS = PRTHS +  ZSRC
 !$acc end kernels 
 !$acc end data  
 !
-IF (LBUDGET_TH) CALL BUDGET (PRTHS,4,'ADV_BU_RTH')
 !
 ! Turbulence variables
 !
 IF (GTKEALLOC) THEN
 !$acc data pcopyin (PTKET) pcopy(PRTKES)
    CALL ADVEC_PPM_ALGO(HMET_ADV_SCHEME, HLBCX, HLBCY, IGRID, PTKET,PRHODJ,PTSTEP, &
-                       ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+                       PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1, PRHOZ2, &
                        ZSRC, KTCOUNT, PCRU, PCRV, PCRW)
    !$acc kernels
    PRTKES = PRTKES + ZSRC
    !$acc end kernels
 !$acc end data 
 !
-  IF (LBUDGET_TKE) CALL BUDGET (PRTKES,5,'ADV_BU_RTKE')
 !
 END IF
 !
+!
+!
 ! Case with KRR moist variables
 !
 DO JRR=1,KRR
+
    !$acc data pcopy(PRRS(:,:,:,JRR)) pcopyin(PRT(:,:,:,JRR))
    CALL ADVEC_PPM_ALGO(HMET_ADV_SCHEME, HLBCX, HLBCY, IGRID, PRT(:,:,:,JRR), &
                        PRHODJ, PTSTEP, &
-                       ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+                       PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1, PRHOZ2, &
                        ZSRC, KTCOUNT, PCRU, PCRV, PCRW)
    !$acc kernels
    PRRS(:,:,:,JRR) = PRRS(:,:,:,JRR) + ZSRC(:,:,:)
    !$acc end kernels
-!$acc end data
-!
-   IF (JRR==1.AND.LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1),6,'ADV_BU_RRV') 
-   IF (JRR==2.AND.LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'ADV_BU_RRC') 
-   IF (JRR==3.AND.LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3),8,'ADV_BU_RRR') 
-   IF (JRR==4.AND.LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'ADV_BU_RRI') 
-   IF (JRR==5.AND.LBUDGET_RS) CALL BUDGET (PRRS(:,:,:,5),10,'ADV_BU_RRS') 
-   IF (JRR==6.AND.LBUDGET_RG) CALL BUDGET (PRRS(:,:,:,6),11,'ADV_BU_RRG') 
-   IF (JRR==7.AND.LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'ADV_BU_RRH') 
+   !$acc end data
 !
 END DO
-!$acc end data
 !
 END SUBROUTINE PPM_MET_D
 
diff --git a/MNH/ppm_rhodj.f90 b/MNH/ppm_rhodj.f90
index 89d59a4de0789d89594f2598d0b622f05c6d1377..90e8e51917158639e781c83b123a652ce534e094 100644
--- a/MNH/ppm_rhodj.f90
+++ b/MNH/ppm_rhodj.f90
@@ -107,12 +107,12 @@ REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZUNIT
 IGRID = 1
 !
 ZUNIT = 1.0
-PRHOX1 = PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHODJ, PTSTEP)
-PRHOY1 = PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, PRHOX1, PTSTEP)
-PRHOZ1 = PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHOY1, PTSTEP)
-PRHOZ2 = PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHODJ, PTSTEP)
-PRHOY2 = PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, PRHOZ2, PTSTEP)
-PRHOX2 = PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHOY2, PTSTEP)
+CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHODJ, PTSTEP,PRHOX1)
+CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, PRHOX1, PTSTEP,PRHOY1)
+CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHOY1, PTSTEP,PRHOZ1)
+CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHODJ, PTSTEP,PRHOZ2)
+CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, PRHOZ2, PTSTEP,PRHOY2)
+CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHOY2, PTSTEP,PRHOX2)
 !
 !
 END SUBROUTINE PPM_RHODJ
diff --git a/MNH/ppm_scalar.f90 b/MNH/ppm_scalar.f90
index e5000d294c441e0660db7f170a4382fc4d5eaabe..cd4c45b688f4c5f18dacf682945971ba7b5f35a4 100644
--- a/MNH/ppm_scalar.f90
+++ b/MNH/ppm_scalar.f90
@@ -1,3 +1,7 @@
+!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.
 !
 !
 !     #####################
@@ -8,7 +12,9 @@ INTERFACE
 !
       SUBROUTINE PPM_SCALAR (HLBCX,HLBCY, KSV, KTCOUNT,   &
                      PCRU, PCRV, PCRW, PTSTEP, PRHODJ,    &
-                     PSVT, PRSVS, HSV_ADV_SCHEME  ) 
+                     PRHOX1, PRHOX2, PRHOY1, PRHOY2,      &
+                     PRHOZ1, PRHOZ2,                      &
+                     PSVT, PRSVS, HSV_ADV_SCHEME          ) 
 !
 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 !
@@ -23,6 +29,10 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! Courant
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  ! numbers
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! 
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ ! density
+! Temporary advected rhodj
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOX1,PRHOX2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOY1,PRHOY2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOZ1,PRHOZ2
 !
 REAL,                     INTENT(IN)    :: PTSTEP ! Time step 
 !
@@ -38,12 +48,14 @@ END INTERFACE
 END MODULE MODI_PPM_SCALAR
 !
 !     ######################################################################
+
       SUBROUTINE PPM_SCALAR (HLBCX,HLBCY, KSV, KTCOUNT,  &
                      PCRU, PCRV, PCRW, PTSTEP, PRHODJ,   &
+                     PRHOX1, PRHOX2, PRHOY1, PRHOY2,     &
+                     PRHOZ1, PRHOZ2,                     &
                      PSVT, PRSVS, HSV_ADV_SCHEME  ) 
 
   USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
-  USE MODE_MNH_ZWORK, ONLY : ZUNIT3D
   IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -59,6 +71,10 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! contravariant
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  !  components
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! of momentum
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ ! density
+! Temporary advected rhodj
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOX1,PRHOX2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOY1,PRHOY2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOZ1,PRHOZ2
 !
 REAL,                     INTENT(IN)    :: PTSTEP ! Time step 
 !
@@ -66,28 +82,27 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN)    :: PSVT
 !
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS         ! Source terms
 
-INTEGER :: IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2
+INTEGER :: IZSRC
 
-         CALL  MNH_GET_ZT3D(IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2)
+         CALL  MNH_GET_ZT3D(IZSRC)
 
          CALL  PPM_SCALAR_D (HLBCX,HLBCY, KSV, KTCOUNT,  &
                            & PCRU, PCRV, PCRW, PTSTEP, PRHODJ,  &
+                           & PRHOX1, PRHOX2, PRHOY1, PRHOY2,     &
+                           & PRHOZ1, PRHOZ2,                     &
                            & PSVT, PRSVS, HSV_ADV_SCHEME, &
-                           & ZT3D(:,:,:,IZSRC), &
-                           & ZT3D(:,:,:,IZRHOX1),ZT3D(:,:,:,IZRHOX2),ZT3D(:,:,:,IZRHOY1), &
-                           & ZT3D(:,:,:,IZRHOY2),ZT3D(:,:,:,IZRHOZ1),ZT3D(:,:,:,IZRHOZ2), &
-                           & ZUNIT3D ) 
+                           & ZT3D(:,:,:,IZSRC) ) 
 
-         CALL  MNH_REL_ZT3D(IZSRC,IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2)
+         CALL  MNH_REL_ZT3D(IZSRC)
 
 CONTAINS
 
       SUBROUTINE PPM_SCALAR_D (HLBCX,HLBCY, KSV, KTCOUNT,  &
                      & PCRU, PCRV, PCRW, PTSTEP, PRHODJ,   &
+                     & PRHOX1, PRHOX2, PRHOY1, PRHOY2,     &
+                     & PRHOZ1, PRHOZ2,                     &
                      & PSVT, PRSVS, HSV_ADV_SCHEME, &
-                     & ZSRC, &
-                     & ZRHOX1,ZRHOX2,ZRHOY1,ZRHOY2,ZRHOZ1,ZRHOZ2, &
-                     & ZUNIT ) 
+                     & ZSRC ) 
 
 !     ######################################################################
 !
@@ -117,6 +132,7 @@ CONTAINS
 !!    MODIFICATIONS
 !!    -------------
 !!      Original 11.05.2006. T.Maric
+!!      Modification : 11.2011 C.Lac, V.Masson : Advection of (theta_l,r_t) 
 !!
 !-------------------------------------------------------------------------------
 !
@@ -127,14 +143,13 @@ CONTAINS
 !
 USE MODD_PARAMETERS
 USE MODD_CONF
-USE MODD_BUDGET
 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
 !
 USE MODI_SHUMAN
-USE MODI_BUDGET
 USE MODI_PPM
 USE MODI_ADVEC_PPM_ALGO
 !
+!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -150,6 +165,10 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRU  ! contravariant
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRV  !  components
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PCRW  ! of momentum
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ ! density
+! Temporary advected rhodj
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOX1,PRHOX2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOY1,PRHOY2
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHOZ1,PRHOZ2
 !
 REAL,                     INTENT(IN)    :: PTSTEP ! Time step 
 !
@@ -170,13 +189,6 @@ INTEGER :: IGRID ! localisation on the model grid
 REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZSRC
 !$acc declare present (ZSRC)
 !
-! Temporary advected rhodj
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOX1,ZRHOX2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOY1,ZRHOY2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZRHOZ1,ZRHOZ2
-REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZUNIT
-!$acc declare present (ZRHOX1,ZRHOX2,ZRHOY1,ZRHOY2,ZRHOZ1,ZRHOZ2,ZUNIT)
-!
 !-------------------------------------------------------------------------------
 !
 !*       1.     CALL THE ADVEC_PPM_ALGO ROUTINE FOR EACH FIELD
@@ -184,35 +196,22 @@ REAL, DIMENSION(SIZE(PCRU,1),SIZE(PCRU,2),SIZE(PCRU,3)) :: ZUNIT
 !
 IGRID = 1
 !
-! Calculate the advection of the density RHODJ to pass to the algorithm
-!
-!ZUNIT = 1.0
-!$acc data pcopyin (PCRU,PRHODJ) pcopyin(PCRV,PCRW) 
-CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, PRHODJ, PTSTEP, ZRHOX1)
-CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, ZRHOX1, PTSTEP, ZRHOY1)
-CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, ZRHOY1, PTSTEP, ZRHOZ1)
-CALL PPM_S0_Z(IGRID, ZUNIT, PCRW, PRHODJ, PTSTEP, ZRHOZ2)
-CALL PPM_S0_Y(HLBCY, IGRID, ZUNIT, PCRV, ZRHOZ2, PTSTEP, ZRHOY2)
-CALL PPM_S0_X(HLBCX, IGRID, ZUNIT, PCRU, ZRHOY2, PTSTEP, ZRHOX2)
-!
 ! Case with KSV tracers
 !
 DO JSV=1,KSV
 !$acc data pcopyin (PSVT(:,:,:,JSV)) 
    CALL ADVEC_PPM_ALGO(HSV_ADV_SCHEME, HLBCX, HLBCY, IGRID, PSVT(:,:,:,JSV), & 
                        PRHODJ, PTSTEP, & 
-                       ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+                       PRHOX1, PRHOX2, PRHOY1, PRHOY2, PRHOZ1, PRHOZ2, &
                        ZSRC, KTCOUNT, PCRU, PCRV, PCRW)
 ! add the advection to the sources
    PRSVS(:,:,:,JSV) =  PRSVS(:,:,:,JSV) + ZSRC(:,:,:)   
 !
-   IF (LBUDGET_SV) CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'ADV_BU_RSV')
 !
 !$acc end data 
 !
 END DO
 !
-!$acc end data 
 !
 END SUBROUTINE PPM_SCALAR_D
 
diff --git a/MNH/set_ref.f90 b/MNH/set_ref.f90
index f20557afd17ec4c0161bb775e0a6c2f54becaf2b..1584f2993d582e86563b606f01f86ee104116b9a 100644
--- a/MNH/set_ref.f90
+++ b/MNH/set_ref.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/set_ref.f90,v $ $Revision: 1.2.4.3.2.1 $ $Date: 2011/07/07 14:38:37 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/set_ref.f90,v $ $Revision: 1.2.4.3.2.1.10.1.2.2 $ $Date: 2014/01/09 15:01:57 $
 !-----------------------------------------------------------------
 !###################
 MODULE MODI_SET_REF
@@ -146,6 +150,7 @@ END MODULE MODI_SET_REF
 !!                                point under the ground
 !!      Modification    03/12/02  (P. Jabouille)  add no thinshell condition
 !!      Modification    05/06     Remove the 'DAVI' type of lbc
+!!      Modification    07/13     (J.Colin) Special case for LBOUSS=T 
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -155,6 +160,7 @@ USE MODD_PARAMETERS
 USE MODD_REF
 USE MODD_CST
 USE MODD_CONF
+USE MODD_REF, ONLY : LBOUSS
 !
 USE MODE_ll
 !JUAN REALZ
@@ -381,8 +387,12 @@ CALL MPPDB_CHECK3D(PEXNREF,"SET_REF::PEXNREF",PRECISION)
 !
 !
 ZCVD_O_RD = (XCPD / XRD) - 1.
-ZRHOREF(:,:,:) = PEXNREF(:,:,:) ** ZCVD_O_RD * XP00 / ( XRD * PTHVREF(:,:,:) )
-ZRHOREF(:,:,1)=ZRHOREF(:,:,2)  ! this avoids to obtain erroneous values for
+IF (LBOUSS) THEN
+  ZRHOREF(:,:,:) = PRHODREF(:,:,:)
+ELSE
+  ZRHOREF(:,:,:) = PEXNREF(:,:,:) ** ZCVD_O_RD * XP00 / ( XRD * PTHVREF(:,:,:) )
+  ZRHOREF(:,:,1)=ZRHOREF(:,:,2)  ! this avoids to obtain erroneous values for
+END IF
                                ! rv at this last point
 !
 IF ( CEQNSYS == 'DUR' ) THEN
diff --git a/MNH/shuman.f90 b/MNH/shuman.f90
index 562d6b07d9656b2fbfc85d6344a3b3b9f10ff9dc..052b0644c8f5988c2f4dde27e4207c409d7ce557 100644
--- a/MNH/shuman.f90
+++ b/MNH/shuman.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/shuman.f90,v $ $Revision: 1.1.10.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/shuman.f90,v $ $Revision: 1.1.10.1.16.1.2.2 $
 ! MASDEV4_7 operators 2006/05/18 13:07:25
 !-----------------------------------------------------------------
 !     ##################
@@ -38,14 +42,18 @@ REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDYM   ! result at flux
                                                             ! side
 END FUNCTION DYM
 !
-FUNCTION DZF(PA)  RESULT(PDZF)
+FUNCTION DZF(KKA,KKU,KL,PA)  RESULT(PDZF)
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at flux
                                                             !  side
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF   ! result at mass
                                                             ! localization 
 END FUNCTION DZF
 !
-FUNCTION DZM(PA)  RESULT(PDZM)
+FUNCTION DZM(KKA,KKU,KL,PA)  RESULT(PDZM)
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass
                                                             ! localization
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZM   ! result at flux
@@ -63,7 +71,7 @@ FUNCTION MXM(PA)  RESULT(PMXM)
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass localization
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMXM   ! result at flux localization 
 END FUNCTION MXM
-!
+
 FUNCTION MYF(PA)  RESULT(PMYF)
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at flux
                                                             !   side
@@ -76,14 +84,18 @@ REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass l
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMYM   ! result at flux localization 
 END  FUNCTION MYM
 !
-FUNCTION MZF(PA)  RESULT(PMZF)
+FUNCTION MZF(KKA,KKU,KL,PA)  RESULT(PMZF)
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at flux
                                                             !  side
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF   ! result at mass
                                                             ! localization 
 END FUNCTION MZF
 !
-FUNCTION MZM(PA)  RESULT(PMZM)
+FUNCTION MZM(KKA,KKU,KL,PA)  RESULT(PMZM)
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass localization
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM   ! result at flux localization 
 END FUNCTION MZM
@@ -478,7 +490,7 @@ PMYM(:,1,:)    = PMYM(:,IJU-2*JPHEXT+1,:)
 !
 END FUNCTION MYM
 !     ###############################
-      FUNCTION MZF(PA)  RESULT(PMZF)
+      FUNCTION MZF(KKA,KKU,KL,PA)  RESULT(PMZF)
 !     ###############################
 !
 !!****  *MZF* -  Shuman operator : mean operator in z direction for a 
@@ -528,6 +540,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at flux
                                                             !  side
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZF   ! result at mass
@@ -572,7 +586,7 @@ END DO
 !
 END FUNCTION MZF
 !     ###############################
-      FUNCTION MZM(PA)  RESULT(PMZM)
+      FUNCTION MZM(KKA,KKU,KL,PA)  RESULT(PMZM)
 !     ###############################
 !
 !!****  *MZM* -  Shuman operator : mean operator in z direction for a 
@@ -622,6 +636,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass localization
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM   ! result at flux localization 
 !
@@ -1060,7 +1076,7 @@ PDYM(:,1,:)    =  PDYM(:,IJU-2*JPHEXT+1,:)
 !
 END FUNCTION DYM
 !     ###############################
-      FUNCTION DZF(PA)  RESULT(PDZF)
+      FUNCTION DZF(KKA,KKU,KL,PA)  RESULT(PDZF)
 !     ###############################
 !
 !!****  *DZF* -  Shuman operator : finite difference operator in z direction
@@ -1110,6 +1126,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at flux
                                                             !  side
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZF   ! result at mass
@@ -1154,7 +1172,7 @@ END DO
 !
 END FUNCTION DZF
 !     ###############################
-      FUNCTION DZM(PA)  RESULT(PDZM)
+      FUNCTION DZM(KKA,KKU,KL,PA)  RESULT(PDZM)
 !     ###############################
 !
 !!****  *DZM* -  Shuman operator : finite difference operator in z direction
@@ -1204,6 +1222,8 @@ IMPLICIT NONE
 !*       0.1   Declarations of argument and result
 !              ------------------------------------
 !
+INTEGER,              INTENT(IN)                  :: KKA, KKU ! near ground and uppest atmosphere array indexes
+INTEGER,              INTENT(IN)                  :: KL     ! +1 if grid goes from ground to atmosphere top, -1 otherwise
 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass
                                                             ! localization
 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PDZM   ! result at flux
diff --git a/MNH/update_metrics.f90 b/MNH/update_metrics.f90
index 8b663bb145e3090454ebd859214815205065af6e..bd4ef68d2830ef4be6e060f81d8da673b0d73c1f 100644
--- a/MNH/update_metrics.f90
+++ b/MNH/update_metrics.f90
@@ -1,7 +1,11 @@
+!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: /home/cvsroot/MNH-VX-Y-Z/src/MNH/update_metrics.f90,v $ $Revision: 1.1.4.1 $
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/update_metrics.f90,v $ $Revision: 1.1.4.1.18.3 $
 ! MASDEV4_7 newsrc 2006/05/18 13:07:25
 !-----------------------------------------------------------------
 !     ###################
@@ -56,6 +60,7 @@ END MODULE MODI_UPDATE_METRICS
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    april 2006
+!!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -104,7 +109,7 @@ NULLIFY(TZMETRICS_ll)
 !            -------------
 !
 !
-IF(NHALO == 1) THEN
+!!$IF(NHALO == 1) THEN
   CALL ADD3DFIELD_ll(TZMETRICS_ll,PDXX)
   CALL ADD3DFIELD_ll(TZMETRICS_ll,PDYY)
   CALL ADD3DFIELD_ll(TZMETRICS_ll,PDZX)
@@ -112,7 +117,7 @@ IF(NHALO == 1) THEN
   CALL ADD3DFIELD_ll(TZMETRICS_ll,PDZZ)
   CALL UPDATE_HALO_ll(TZMETRICS_ll,IINFO_ll)
   CALL CLEANLIST_ll(TZMETRICS_ll)
-END IF
+!!$END IF
 !
 !-------------------------------------------------------------------------------
 !