diff --git a/MY_RUN/KTEST/013_Iroise_ideal_case_coupling/A_RUN_MNH_TOY/run_mesonh_xyz b/MY_RUN/KTEST/013_Iroise_ideal_case_coupling/A_RUN_MNH_TOY/run_mesonh_xyz
index 61d9446dc67cdeb183afe66005dfed70c7e5fec7..0e64c7a0705e944d9fe35db4186926307274a48b 100755
--- a/MY_RUN/KTEST/013_Iroise_ideal_case_coupling/A_RUN_MNH_TOY/run_mesonh_xyz
+++ b/MY_RUN/KTEST/013_Iroise_ideal_case_coupling/A_RUN_MNH_TOY/run_mesonh_xyz
@@ -25,4 +25,5 @@ ln -sf ../2_INPUT_TOY/TOYNAMELIST.nam_${TYPE_TOY} TOYNAMELIST.nam
 #~~~~~ OASIS
 ln -fs ../3_INPUT_OASIS/namcouple_${TYPE_TOY} namcouple
 
-time Mpirun -np 1 MESONH${XYZ} : -np 1 $PATH_EXETOY/toy_model
+#time Mpirun -np 1 MESONH${XYZ} : -np 1 $PATH_EXETOY/toy_model
+time Mpirun -np 1 $PATH_EXETOY/toy_model : -np 1 totalview MESONH${XYZ}  
diff --git a/src/LIB/SURCOUCHE/src/mode_scatter.f90 b/src/LIB/SURCOUCHE/src/mode_scatter.f90
index 32eaa729052794702023468e113cecabaa913449..b81b5a29422b5bdcd6e3fdb2f40eb2bd7cde2c1b 100644
--- a/src/LIB/SURCOUCHE/src/mode_scatter.f90
+++ b/src/LIB/SURCOUCHE/src/mode_scatter.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -7,6 +7,7 @@
 !  J. Escobar  10/02/2012: bug in MPI_RECV: replace MPI_STATUSES_IGNORE with MPI_STATUS_IGNORE
 !  P. Wautelet 26/04/2019: use modd_precision parameters for datatypes of MPI communications
 !  P. Wautelet 25/06/2019: added IO_Field_read for 3D integer arrays (SCATTERXX_N3 and SCATTERXY_N3)
+!  J. Escobar  21/07/2020: for reduction of MPI_BUFFER_SIZE, in SCATTERXY_X3 replace MPI_BSEND -> MPI_ISEND
 !-----------------------------------------------------------------
 
 MODULE MODE_SCATTER_ll
@@ -467,7 +468,8 @@ END IF
 END SUBROUTINE  SCATTERXY_X2
 
 SUBROUTINE SCATTERXY_X3(PSEND,PRECV,KROOT,KCOMM)
-USE MODD_IO, ONLY: ISP, ISNPROC
+USE MODD_IO,     ONLY: ISP, ISNPROC
+USE MODD_VAR_ll, ONLY: MNH_STATUSES_IGNORE
 
 REAL,DIMENSION(:,:,:),TARGET,INTENT(IN) :: PSEND
 REAL,DIMENSION(:,:,:),       INTENT(INOUT):: PRECV
@@ -475,25 +477,45 @@ INTEGER,                     INTENT(IN) :: KROOT
 INTEGER,                     INTENT(IN) :: KCOMM
 
 INTEGER :: IERR
-INTEGER :: JI
+INTEGER :: JI,JKU
 INTEGER :: IXO,IXE,IYO,IYE
 REAL,DIMENSION(:,:,:), POINTER :: TX3DP
+
+INTEGER,ALLOCATABLE,DIMENSION(:)    :: REQ_TAB
+INTEGER                             :: NB_REQ
+TYPE TX_3DP
+   REAL,DIMENSION(:,:,:), POINTER    :: X
+END TYPE TX_3DP
+TYPE(TX_3DP),ALLOCATABLE,DIMENSION(:) :: T_TX3DP
+
+JKU = SIZE(PSEND,3)
   
 IF (ISP == KROOT) THEN
-  DO JI = 1,ISNPROC
-    CALL GET_DOMREAD_ll(JI,IXO,IXE,IYO,IYE)
-    TX3DP=>PSEND(IXO:IXE,IYO:IYE,:)
-    
-    IF (ISP /= JI) THEN 
-      CALL MPI_BSEND(TX3DP,SIZE(TX3DP),MNHREAL_MPI,JI-1,199+KROOT,KCOMM&
-           & ,IERR)
-    ELSE 
-      PRECV(:,:,:) = TX3DP(:,:,:)
-    END IF
-  END DO
+   NB_REQ=0
+   ALLOCATE(REQ_TAB(ISNPROC-1))
+   ALLOCATE(T_TX3DP(ISNPROC-1))   
+   DO JI = 1,ISNPROC
+      CALL GET_DOMREAD_ll(JI,IXO,IXE,IYO,IYE)
+      TX3DP=>PSEND(IXO:IXE,IYO:IYE,:)
+      IF (ISP /= JI) THEN
+         NB_REQ = NB_REQ + 1
+         ALLOCATE(T_TX3DP(NB_REQ)%X(IXO:IXE,IYO:IYE,JKU))
+         T_TX3DP(NB_REQ)%X=TX3DP
+         CALL MPI_ISEND(T_TX3DP(NB_REQ)%X,SIZE(TX3DP),MNHREAL_MPI,JI-1,199+KROOT,KCOMM&
+               & ,REQ_TAB(NB_REQ),IERR)
+         !CALL MPI_BSEND(TX3DP,SIZE(TX3DP),MNHREAL_MPI,JI-1,199+KROOT,KCOMM&
+         !     & ,IERR)        
+      ELSE 
+         PRECV(:,:,:) = TX3DP(:,:,:)
+      END IF
+   END DO
+   IF (NB_REQ .GT.0 ) THEN
+      CALL MPI_WAITALL(NB_REQ,REQ_TAB,MNH_STATUSES_IGNORE,IERR)
+      DO JI=1,NB_REQ ;  DEALLOCATE(T_TX3DP(JI)%X) ; ENDDO
+   END IF
 ELSE
-  CALL MPI_RECV(PRECV,SIZE(PRECV),MNHREAL_MPI,KROOT-1,199+KROOT,KCOMM&
-       & ,MPI_STATUS_IGNORE,IERR)
+   CALL MPI_RECV(PRECV,SIZE(PRECV),MNHREAL_MPI,KROOT-1,199+KROOT,KCOMM&
+        & ,MPI_STATUS_IGNORE,IERR)
 END IF
 
 END SUBROUTINE  SCATTERXY_X3
diff --git a/src/MNH/call_rttov11.f90 b/src/MNH/call_rttov11.f90
index c4eb7ea2c9a4744f079f982cc3787edbe180b4a4..68864a5cef7e02c1312f67678ce31f45449c2344 100644
--- a/src/MNH/call_rttov11.f90
+++ b/src/MNH/call_rttov11.f90
@@ -82,7 +82,8 @@ SUBROUTINE CALL_RTTOV11(KDLON, KFLEV, PEMIS, PTSRAD,     &
 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
 !!      JP Chaboureau 30/05/2017 exclude the first layer when considering clouds
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+!  P. Wautelet   10/04/2019: replace ABORT and STOP calls by Print_msg
+!  JP Chaboureau 26/10/2020: calculate all IR intruments; deallocate MW tabs
 !!----------------------------------------------------------------------------
 !!
 !!*       0.    DECLARATIONS
@@ -109,7 +110,7 @@ USE MODE_POS
 !
 #ifdef MNH_RTTOV_11
 USE rttov_const, ONLY :  &
-       & sensor_id_ir, sensor_id_hi, sensor_id_mw, &
+       & sensor_id, sensor_id_ir, sensor_id_hi, sensor_id_mw, &
        & q_mixratio_to_ppmv, tmin, tmax, qmin, qmax, pmin, pmax
 USE rttov_types
 USE parkind1, ONLY: jpim, jprb, jplm
@@ -192,8 +193,8 @@ REAL :: ZZH, zdeg_to_rad, zrad_to_deg, zbeta, zalpha
 ! REAL, DIMENSION(:,:),   ALLOCATABLE :: ZCOSZEN, ZSINZEN, ZAZIMSOL
 
 ! -----------------------------------------------------------------------------
-REAL, DIMENSION(:), ALLOCATABLE :: ZANGL   !Satellite zenith angle (deg)
-REAL, DIMENSION(:), ALLOCATABLE :: ZANGS   !Solar zenith angle (deg)
+REAL, DIMENSION(1) :: ZANGL, ZLON, ZLAT   !Satellite zenith angle, longitude, latitude (deg)
+REAL :: ZANGS   !Solar zenith angle (deg)
 ! -----------------------------------------------------------------------------
 ! INDEXES AND TEMPORAL ARRAYS FOR VECTORIZATION
 INTEGER :: JIS, IBEG, IEND, IDIM, ICPT
@@ -296,9 +297,7 @@ DO JSAT=1,IJSAT ! loop over sensors
   instrument(3)=KRTTOVINFO(3,JSAT)
 ! PRINT *,' JSAT=',JSAT, instrument
 
-!!! METEOSAT, GOES, OR MSG PLATFORM
-  IF (KRTTOVINFO(1,JSAT) == 3 .OR. KRTTOVINFO(1,JSAT) == 4 &
-       .OR. KRTTOVINFO(1,JSAT) == 12) THEN
+  IF( sensor_id( instrument(3) ) /= sensor_id_mw) THEN
     opts % rt_ir % addsolar = .FALSE.         ! Do not include solar radiation
     opts % interpolation % addinterp  = .TRUE.  ! Allow interpolation of input profile
     opts % interpolation % interp_mode = 1       ! Set interpolation method
@@ -342,7 +341,7 @@ DO JSAT=1,IJSAT ! loop over sensors
 
   ALLOCATE(ZBT(IIU,IJU,nchannels))
   ZBT(:,:,:)=999.
-!  PRINT *,'ncan=',nchan,' nchannels=',nchannels
+! PRINT *,'ncan=',nchan,' nchannels=',nchannels
 
   ALLOCATE (chanprof     (nchannels))
   ALLOCATE (frequencies  (nchannels))
@@ -354,9 +353,7 @@ DO JSAT=1,IJSAT ! loop over sensors
   calcemis  = .TRUE.
   emissivity % emis_in = 0.0_JPRB
 
-!!! METEOSAT, GOES, OR MSG PLATFORM
-  IF (KRTTOVINFO(1,JSAT) == 3 .OR. KRTTOVINFO(1,JSAT) == 4 &
-       .OR. KRTTOVINFO(1,JSAT) == 12) calcemis = .FALSE.
+  IF( coef_rttov%coef% id_sensor /= sensor_id_mw)  calcemis = .FALSE.
 
 !  IF( coef_rttov%coef% id_sensor /= sensor_id_mw) THEN
 !  ! Allocate arrays for surface reflectance
@@ -426,18 +423,41 @@ DO JSAT=1,IJSAT ! loop over sensors
 
 !!  opts%interpolation%reg_limit_extrap = .TRUE.
 !!  profiles(1)%gas_units = 1 ! kg/kg over moist air
-!PRINT *,'nlev=',nlev,' tmax=',tmax,' tmin=',tmin,' qmax=',qmax,' qmin=',qmin
-!PRINT *, coef_rttov%coef % nlevels
+! PRINT *,'nlev=',nlev,' tmax=',tmax,' tmin=',tmin,' qmax=',qmax,' qmin=',qmin
+! PRINT *, coef_rttov%coef % nlevels
   DO JI=IIB,IIE
     DO JJ=IJB,IJE      
+      ZANGL = XUNDEF
+      ZLON  = XLON(JI,JJ)
+      ZLAT  = XLAT(JI,JJ)
+      IF (KRTTOVINFO(1,JSAT) == 2) THEN ! DMSP PLATFORM
+        ZANGL=53.1 ! see Saunders, 2002, RTTOV7 - science/validation rep, page 8
+      ELSEIF (KRTTOVINFO(1,JSAT) == 3) THEN ! METEOSAT PLATFORM
+        CALL DETER_ANGLE(5, 1, ZLAT, ZLON, ZANGL)
+        WHERE (ZANGL /= XUNDEF .AND. ZANGL /=0.) ZANGL=ACOS(1./ZANGL)*180./XPI
+      ELSEIF (KRTTOVINFO(1,JSAT) == 12) THEN ! MSG PLATFORM
+        CALL DETER_ANGLE(6, 1, ZLAT, ZLON, ZANGL)
+        WHERE (ZANGL /= XUNDEF .AND. ZANGL /=0.) ZANGL=ACOS(1./ZANGL)*180./XPI
+      ELSEIF (KRTTOVINFO(1,JSAT) == 4) THEN ! GOES-E PLATFORM
+        CALL DETER_ANGLE(1, 1, ZLAT, ZLON, ZANGL)
+        WHERE (ZANGL /= XUNDEF .AND. ZANGL /=0.) ZANGL=ACOS(1./ZANGL)*180./XPI
+      ELSEIF (KRTTOVINFO(1,JSAT) == 7) THEN ! TRMM PLATFORM
+        ZANGL=52.3
+      ELSE
+        ZANGL=0.
+      ENDIF
+! Coefficients computed from transmittances for 6 viewing angles in the range 
+! 0 to 63.6 deg (Saunders, 2002, RTTOV7 - science/validation rep., page 3)
+      profiles(1) % zenangle = MIN(ZANGL(1),65.)
+
       DO JK=IKB,IKE ! nlevels
         JKRAD = nlev-JK+2 !INVERSION OF VERTICAL LEVELS!
-!PRINT *,'jk=',jk,' jkrad=',jkrad
+!       PRINT *,'jk=',jk,' jkrad=',jkrad
         profiles(1) % p(JKRAD) = PPABST(JI,JJ,JK)*0.01
         profiles(1) % t(JKRAD) = MIN(tmax,MAX(tmin,ZTEMP(JI,JJ,JK)))
-!PRINT *,'jk=',JK,' ZTEMP=',ZTEMP(JI,JJ,JK),' t=',profiles(1) % t(JKRAD)
+!       PRINT *,'jk=',JK,' ZTEMP=',ZTEMP(JI,JJ,JK),' t=',profiles(1) % t(JKRAD)
         profiles(1) % q(JKRAD) = MIN(qmax,MAX(qmin,PRT(JI,JJ,JK,1)*q_mixratio_to_ppmv))
-!        PRINT *,JK,profiles(1) % p(JKRAD) ,profiles(1) % t(JKRAD) ,profiles(1) % q(JKRAD) 
+!       PRINT *,JK,profiles(1) % p(JKRAD) ,profiles(1) % t(JKRAD) ,profiles(1) % q(JKRAD) 
       END DO
       profiles(1) % elevation = 0.5*( PZZ(JI,JJ,1)+PZZ(JI,JJ,IKB) )
       profiles(1) % skin % t = MIN(tmax,MAX(tmin,PTSRAD(JI,JJ)))
@@ -468,7 +488,7 @@ DO JSAT=1,IJSAT ! loop over sensors
       ELSE
         DO JK=IKB,IKE
           JKRAD = nlev-JK+2 !INVERSION OF VERTICAL LEVELS!
-          cld_profiles(1) % ph (JKRAD) = 0.5*( PPABST(JI,JJ,JK) + PPABST(JI,JJ,JK+1) )*0.01
+          cld_profiles(1) %ph (JKRAD) = 0.5*( PPABST(JI,JJ,JK) + PPABST(JI,JJ,JK+1) )*0.01
           cld_profiles(1) %cc(JKRAD) = PCLDFR(JI,JJ,JK)
           cld_profiles(1) %clw(JKRAD) = MIN(ZRCMAX,PRT(JI,JJ,JK,2))
           cld_profiles(1) %rain(JKRAD) = MIN(ZRRMAX,PRT(JI,JJ,JK,3))
@@ -478,14 +498,14 @@ DO JSAT=1,IJSAT ! loop over sensors
           END IF
         END DO
         cld_profiles (1) % ph (nlev+1) =   profiles (1) % s2m % p
-!          PRINT *,nlev+1,' cld_profiles(1) % ph (nlev+1) =',cld_profiles(1) % ph (nlev+1) 
+!       PRINT *,nlev+1,' cld_profiles(1) % ph (nlev+1) =',cld_profiles(1) % ph (nlev+1) 
       END IF
 
       DO JCH=1,nchannels
         IF (.NOT.calcemis(JCH)) emissivity(JCH)%emis_in = PEMIS(JI,JJ)
       END DO
 
-!write(*,*) 'Calling forward model' 
+! write(*,*) 'Calling forward model' 
 
 ! Forward model run
       IF ( coef_rttov%coef% id_sensor /= sensor_id_mw) THEN
@@ -573,12 +593,18 @@ DO JSAT=1,IJSAT ! loop over sensors
     TZFIELD%NGRID      = 1
     TZFIELD%NTYPE      = TYPEREAL
     TZFIELD%NDIMS      = 2
-    TZFIELD%LTIMEDEP   = .FALSE.
+    TZFIELD%LTIMEDEP   = .TRUE.
 !    PRINT *,'YRECFM='//TRIM(TZFIELD%CMNHNAME)
     CALL IO_Field_write(TPFILE,TZFIELD,ZBT(:,:,JCH))
   END DO
-  DEALLOCATE(chanprof,frequencies,emissivity,calcemis,profiles,cld_profiles)
+  DEALLOCATE(chanprof,frequencies,emissivity,calcemis,profiles)
   DEALLOCATE(ZBT)
+  IF( coef_rttov%coef% id_sensor == sensor_id_mw) THEN
+    CALL rttov_alloc_scatt_prof(nprof, cld_profiles, nlev, .FALSE., 0_jpim)
+    CALL rttov_dealloc_scattcoeffs(coef_scatt)
+  END IF
+  DEALLOCATE(cld_profiles)
+  CALL rttov_dealloc_coefs(errorstatus, coef_rttov)
 !  IF( coef_rttov%coef% id_sensor /= sensor_id_mw) THEN
 !    DEALLOCATE(calcrefl,reflectance)
 !  END IF
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 15a627d2e58986319092e68f4617d46c543b69fb..141b9d4277a21b2f28f315993257c29e41b8acc0 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -748,7 +748,7 @@ END IF
 !             ---------------------------------------
 !
 IF (KMI == 1) THEN
-  LRED    = .FALSE.
+  LRED    = .TRUE.
   LWARM = .TRUE.
   CPRISTINE_ICE = 'PLAT'
   LSEDIC  = .TRUE.
diff --git a/src/MNH/ini_budget.f90 b/src/MNH/ini_budget.f90
index 48276d3e1703ab792625defad4983951d40e8403..3bcb826a94fbbdf1e6173f47b410a662cc12a2de 100644
--- a/src/MNH/ini_budget.f90
+++ b/src/MNH/ini_budget.f90
@@ -192,10 +192,12 @@ end subroutine Budget_preallocate
 !  P. Wautelet 02-03/2020: use the new data structures and subroutines for budgets
 !  B. Vie      02/03/2020: LIMA negativity checks after turbulence, advection and microphysics budgets
 !  P .Wautelet 09/03/2020: add missing budgets for electricity
+!  P. Wautelet 25/03/2020: add missing ove_relax_grd
 !  P. Wautelet 23/04/2020: add nid in tbudgetdata datatype
 !  P. Wautelet + Benoit Vié 11/06/2020: improve removal of negative scalar variables + adapt the corresponding budgets
 !  P. Wautelet 30/06/2020: use NADVSV when possible
 !  P. Wautelet 30/06/2020: add NNETURSV, NNEADVSV and NNECONSV variables
+!  P. Wautelet 06/07/2020: bugfix: add condition on HTURB for NETUR sources for SV budgets
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
diff --git a/src/MNH/ini_lesn.f90 b/src/MNH/ini_lesn.f90
index 8d9b49aeff6c0e69e67bc50d4267a4b76ad878ad..fb96af501633ab025f2977c1c4cf7772cd24a439 100644
--- a/src/MNH/ini_lesn.f90
+++ b/src/MNH/ini_lesn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2000-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2000-2020 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -36,6 +36,7 @@
 !!                     02/2019 (C. Lac) Add rain fraction as a LES diagnostic
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
+!  P. Wautelet 12/08/2020: bugfix: use NUNDEF instead of XUNDEF for integer variables
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -363,7 +364,7 @@ END IF
 !            --------------------
 !
 IF (ANY(NLES_LEVELS(:)/=NUNDEF)) THEN
-  NLES_K = COUNT (NLES_LEVELS(:)/=XUNDEF)
+  NLES_K = COUNT (NLES_LEVELS(:)/=NUNDEF)
   CLES_LEVEL_TYPE='K'
 ELSE
   IF (NLES_K==0) THEN
@@ -411,8 +412,8 @@ END IF
 !*      5.2  Case of model levels (highest priority)
 !            --------------------
 !
-IF (ANY(NSPECTRA_LEVELS(:)/=XUNDEF)) THEN
-  NSPECTRA_K = COUNT (NSPECTRA_LEVELS(:)/=XUNDEF)
+IF (ANY(NSPECTRA_LEVELS(:)/=NUNDEF)) THEN
+  NSPECTRA_K = COUNT (NSPECTRA_LEVELS(:)/=NUNDEF)
   CSPECTRA_LEVEL_TYPE='K'
 END IF
 !
diff --git a/src/MNH/latlon_to_xy.f90 b/src/MNH/latlon_to_xy.f90
index 320f6b1dc92ca461ff816bac677bb2315ec47805..972999064186796b7d4b1721bdd304e0d51f15c4 100644
--- a/src/MNH/latlon_to_xy.f90
+++ b/src/MNH/latlon_to_xy.f90
@@ -58,6 +58,7 @@
 !  P. Wautelet 07/02/2019: force TYPE to a known value for IO_File_add2list
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 10/04/2020: add missing initializations (LATLON_TO_XY was not working)
+!  J. Escobar  21/07/2020: missing modi_version
 !----------------------------------------------------------------------------
 !
 !*    0.     DECLARATION
@@ -85,6 +86,7 @@ use MODE_SPLITTINGZ_ll
 !
 USE MODI_INI_CST
 USE MODI_READ_HGRID
+USE MODI_VERSION
 !
 USE MODN_CONFIO,           ONLY: NAM_CONFIO
 !
diff --git a/src/MNH/lima.f90 b/src/MNH/lima.f90
index b69451377db3e82ab08c68beb622b8b53a2df952..450913968ef7105daff3b8cbe73e4cddcba5728b 100644
--- a/src/MNH/lima.f90
+++ b/src/MNH/lima.f90
@@ -112,7 +112,6 @@ use modd_budget,          only: lbu_enable,
                                 NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RR, NBUDGET_RI,  &
                                 NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1,             &
                                 tbudgets
-USE MODD_CLOUDPAR_n,      ONLY: NSPLITR, NSPLITG
 USE MODD_CST,             ONLY: XCI, XCL, XCPD, XCPV, XLSTT, XLVTT, XTT, XRHOLW, XP00, XRD
 USE MODD_IO,              ONLY: TFILEDATA
 USE MODD_NSV,             ONLY: NSV_LIMA_BEG,                                                   &
diff --git a/src/MNH/lima_adjust.f90 b/src/MNH/lima_adjust.f90
index 78e66899649101277119c1f82635163cbc9eb8d2..3dbd22694b4e552a82d05eb52b61afaa4415c723 100644
--- a/src/MNH/lima_adjust.f90
+++ b/src/MNH/lima_adjust.f90
@@ -135,9 +135,9 @@ END MODULE MODI_LIMA_ADJUST
 !!      C. Barthe  * LACy*   jan. 2014  add budgets
 !!      JP Chaboureau *LA*   March 2014  fix the calculation of icy cloud fraction
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
 !  P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
-!  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
 !  P. Wautelet 28/05/2020: bugfix: correct array start for PSVT and PSVS
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/lima_cold_slow_processes.f90 b/src/MNH/lima_cold_slow_processes.f90
index aabb1c03deeaa7e7b15639007acdaf81bcf1e62c..9fcacdd5a39b10fc469e99ffba04621f848456a7 100644
--- a/src/MNH/lima_cold_slow_processes.f90
+++ b/src/MNH/lima_cold_slow_processes.f90
@@ -80,6 +80,7 @@ END MODULE MODI_LIMA_COLD_SLOW_PROCESSES
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
 !  P. Wautelet    03/2020: use the new data structures and subroutines for budgets
+!
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
diff --git a/src/MNH/modd_prep_real.f90 b/src/MNH/modd_prep_real.f90
index d2d06321744568f7b2f1caf138f96388ff092c75..6933ae26a8f64486b93121fa44923c6a8a67774a 100644
--- a/src/MNH/modd_prep_real.f90
+++ b/src/MNH/modd_prep_real.f90
@@ -28,6 +28,7 @@
 !!      Original   05/05
 !!                 05/06 (I.Mallet) add *_SV_* variables to allow chemical
 !!                                 initialization from HCHEMFILE
+!!                 09/20 (Q. Rodier) add geopotential height for GFS GRIB read
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -70,6 +71,7 @@ REAL                                  :: XLEN2_LS ! Decay scale for small-scale
 REAL,DIMENSION(:,:),   ALLOCATABLE :: XPS_LS   ! surface pressure
 REAL,DIMENSION(:,:),   ALLOCATABLE :: XZS_LS   ! orography
 REAL,DIMENSION(:,:),   ALLOCATABLE :: XZSMT_LS ! smooth orography
+REAL,DIMENSION(:,:,:),   ALLOCATABLE :: XGH_LS   ! geopotential height
 REAL,DIMENSION(:,:,:), ALLOCATABLE :: XZFLUX_LS! altitude of pressure points
 REAL,DIMENSION(:,:,:), ALLOCATABLE :: XZMASS_LS! altitude of mass points
 REAL,DIMENSION(:,:,:), ALLOCATABLE :: XPMHP_LS ! pressure minus hyd. pressure
diff --git a/src/MNH/mode_prandtl.f90 b/src/MNH/mode_prandtl.f90
index a04e8aa6a4975ee2141b43df4ca538b98cf9d296..04dfe6155e2efdacbcf981be63e260d8046e1187 100644
--- a/src/MNH/mode_prandtl.f90
+++ b/src/MNH/mode_prandtl.f90
@@ -9,6 +9,7 @@
 !
 !* modification 08/2010  V. Masson  smoothing of the discontinuity in functions 
 !                                   used for implicitation of exchange coefficients
+!               05/2020   V. Masson and C. Lac : bug in D_PHI3DTDZ2_O_DDTDZ
 !
 USE MODD_CTURB,      ONLY : XCTV, XCSHF, XCTD, XPHI_LIM, XCPR3, XCPR4, XCPR5
 USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
@@ -272,32 +273,36 @@ IKE = SIZE(PREDTH1,3)-JPVEXT_TURB
 !
 !
 IF (HTURBDIM=='3DIM') THEN
-        !* 3DIM case
-  IF (OUSERV) THEN
-    WHERE (PPHI3(:,:,:)<=XPHI_LIM)
-    D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:)                      &
-          * PDTDZ(:,:,:)*(2.-PREDTH1(:,:,:)*(3./2.+PREDTH1+PREDR1) &
-               /((1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1)))) &
-          + (1.+PREDR1)*(PRED2THR3+PRED2TH3)                       &
-               / (PREDTH1*(1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1))) &
-          - (1./2.*PREDTH1+PREDR1 * (1.+PREDTH1+PREDR1))           &
-               / ((1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1)))
-    ELSEWHERE
-      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:) * 2. * PDTDZ(:,:,:)
-    ENDWHERE
+   ! by derivation of (phi3 dtdz) * dtdz according to dtdz we obtain:
+   D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PDTDZ * (PPHI3 +  &
+           D_PHI3DTDZ_O_DDTDZ(PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,OUSERV) )
 
-!
-  ELSE
-    WHERE (PPHI3(:,:,:)<=XPHI_LIM)
-    D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:)                  &
-          * PDTDZ(:,:,:)*(2.-PREDTH1(:,:,:)*(3./2.+PREDTH1)   &
-               /((1.+PREDTH1)*(1.+1./2.*PREDTH1)))             &
-          + PRED2TH3 / (PREDTH1*(1.+PREDTH1)*(1.+1./2.*PREDTH1)) &
-          - 1./2.*PREDTH1 / ((1.+PREDTH1)*(1.+1./2.*PREDTH1))
-    ELSEWHERE
-      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:) * 2. * PDTDZ(:,:,:)
-    ENDWHERE
-  END IF
+!        !* 3DIM case
+!  IF (OUSERV) THEN
+!    WHERE (PPHI3(:,:,:)<=XPHI_LIM)
+!    D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:)                      &
+!          * PDTDZ(:,:,:)*(2.-PREDTH1(:,:,:)*(3./2.+PREDTH1+PREDR1) &
+!               /((1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1)))) &
+!          + (1.+PREDR1)*(PRED2THR3+PRED2TH3)                       &
+!               / (PREDTH1*(1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1))) &
+!          - (1./2.*PREDTH1+PREDR1 * (1.+PREDTH1+PREDR1))           &
+!               / ((1.+PREDTH1+PREDR1)*(1.+1./2.*(PREDTH1+PREDR1)))
+!    ELSEWHERE
+!      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:) * 2. * PDTDZ(:,:,:)
+!    ENDWHERE
+!
+!!
+!  ELSE
+!    WHERE (PPHI3(:,:,:)<=XPHI_LIM)
+!    D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:)                  &
+!          * PDTDZ(:,:,:)*(2.-PREDTH1(:,:,:)*(3./2.+PREDTH1)   &
+!               /((1.+PREDTH1)*(1.+1./2.*PREDTH1)))             &
+!          + PRED2TH3 / (PREDTH1*(1.+PREDTH1)*(1.+1./2.*PREDTH1)) &
+!          - 1./2.*PREDTH1 / ((1.+PREDTH1)*(1.+1./2.*PREDTH1))
+!    ELSEWHERE
+!      D_PHI3DTDZ2_O_DDTDZ(:,:,:) = PPHI3(:,:,:) * 2. * PDTDZ(:,:,:)
+!    ENDWHERE
+!  END IF
 ELSE
         !* 1DIM case
     WHERE (PPHI3(:,:,:)<=XPHI_LIM)
diff --git a/src/MNH/rain_ice.f90 b/src/MNH/rain_ice.f90
index 47dd05e85865feb0ac12fc8bba044ba2cc7c1d6e..ec1f392049c5b4efdf295c90d3fe66342325628d 100644
--- a/src/MNH/rain_ice.f90
+++ b/src/MNH/rain_ice.f90
@@ -758,8 +758,7 @@ IF( IMICRO >= 0 ) THEN
 !
   CALL RAIN_ICE_SLOW(GMICRO, ZINVTSTEP, ZRHODREF,                      &
                      ZRCT, ZRRT, ZRIT, ZRST, ZRGT, ZRHODJ, ZZT, ZPRES, &
-                     ZLSFACT, ZLVFACT,                                 &
-                     ZSSI, PRHODJ, PTHS, PRVS,                         &
+                     ZLSFACT, ZLVFACT, ZSSI,                           &
                      ZRVS, ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZTHS,         &
                      ZAI, ZCJ, ZKA, ZDV, ZLBDAS, ZLBDAG)
 !
@@ -788,7 +787,7 @@ IF( IMICRO >= 0 ) THEN
     CALL RAIN_ICE_WARM(GMICRO, IMICRO, I1, I2, I3,                                                           &
                        ZRHODREF, ZRVT, ZRCT, ZRRT, ZHLC_HCF, ZHLC_LCF, ZHLC_HRC, ZHLC_LRC,                   &
                        ZRHODJ, ZPRES, ZZT, ZLBDAR, ZLBDAR_RF, ZLVFACT, ZCJ, ZKA, ZDV, ZRF, ZCF, ZTHT, ZTHLT, &
-                       PRHODJ, PTHS, PRVS, ZRVS, ZRCS, ZRRS, ZTHS, ZUSW, PEVAP3D)
+                       ZRVS, ZRCS, ZRRS, ZTHS, ZUSW, PEVAP3D)
   END IF
 !
 !-------------------------------------------------------------------------------
@@ -798,7 +797,7 @@ IF( IMICRO >= 0 ) THEN
 !               ----------------------------------------------
 !
   CALL RAIN_ICE_FAST_RS(PTSTEP, GMICRO, ZRHODREF, ZRVT, ZRCT, ZRRT, ZRST, ZRHODJ, ZPRES, ZZT, &
-                        ZLBDAR, ZLBDAS, ZLSFACT, ZLVFACT, ZCJ, ZKA, ZDV, PRHODJ, PTHS, &
+                        ZLBDAR, ZLBDAS, ZLSFACT, ZLVFACT, ZCJ, ZKA, ZDV, &
                         ZRCS, ZRRS, ZRSS, ZRGS, ZTHS)
 !
 !-------------------------------------------------------------------------------
@@ -809,7 +808,8 @@ IF( IMICRO >= 0 ) THEN
 !
   CALL RAIN_ICE_FAST_RG(KRR, GMICRO, ZRHODREF, ZRVT, ZRCT, ZRRT, ZRIT, ZRST, ZRGT, ZCIT, &
                         ZRHODJ, ZPRES, ZZT, ZLBDAR, ZLBDAS, ZLBDAG, ZLSFACT, ZLVFACT, &
-                        ZCJ, ZKA, ZDV, PRHODJ, PTHS, ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS, ZTHS, &
+                        ZCJ, ZKA, ZDV, &
+                        ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS, ZTHS, &
                         ZUSW, ZRDRYG, ZRWETG)
 !
 !-------------------------------------------------------------------------------
@@ -820,7 +820,7 @@ IF( IMICRO >= 0 ) THEN
 !
  IF ( KRR == 7 ) THEN
   CALL RAIN_ICE_FAST_RH(GMICRO, ZRHODREF, ZRVT, ZRCT, ZRIT, ZRST, ZRGT, ZRHT, ZRHODJ, ZPRES, &
-                        ZZT, ZLBDAS, ZLBDAG, ZLBDAH, ZLSFACT, ZLVFACT, ZCJ, ZKA, ZDV, PRHODJ, PTHS, &
+                        ZZT, ZLBDAS, ZLBDAG, ZLBDAH, ZLSFACT, ZLVFACT, ZCJ, ZKA, ZDV, &
                         ZRCS, ZRRS, ZRIS, ZRSS, ZRGS, ZRHS, ZTHS, ZUSW)
  END IF
 !
@@ -831,7 +831,7 @@ IF( IMICRO >= 0 ) THEN
 !               -------------------------------------------------------------
 !
   CALL RAIN_ICE_FAST_RI(GMICRO, ZRHODREF, ZRIT, ZRHODJ, ZZT, ZSSI, ZLSFACT, ZLVFACT, &
-                        ZAI, ZCJ, PRHODJ, PTHS, ZCIT, ZRCS, ZRIS, ZTHS)
+                        ZAI, ZCJ, ZCIT, ZRCS, ZRIS, ZTHS)
 !
 !
 !-------------------------------------------------------------------------------
diff --git a/src/MNH/rain_ice_fast_rg.f90 b/src/MNH/rain_ice_fast_rg.f90
index 8f2bec61c49b1e1015b200982fa8fb2914294ee6..a2d9689cf78b95b701e842367d1d7a33dcbc7a1a 100644
--- a/src/MNH/rain_ice_fast_rg.f90
+++ b/src/MNH/rain_ice_fast_rg.f90
@@ -22,7 +22,8 @@ CONTAINS
 
 SUBROUTINE RAIN_ICE_FAST_RG(KRR, OMICRO, PRHODREF, PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PCIT, &
                             PRHODJ, PPRES, PZT, PLBDAR, PLBDAS, PLBDAG, PLSFACT, PLVFACT, &
-                            PCJ, PKA, PDV, PRHODJ3D, PTHS3D, PRCS, PRRS, PRIS, PRSS, PRGS, PRHS, PTHS, &
+                            PCJ, PKA, PDV, &
+                            PRCS, PRRS, PRIS, PRSS, PRGS, PRHS, PTHS, &
                             PUSW, PRDRYG, PRWETG)
 
 !
@@ -66,8 +67,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PLVFACT  ! L_v/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PCJ      ! Function to compute the ventilation coefficient
 REAL,     DIMENSION(:),     intent(in)    :: PKA      ! Thermal conductivity of the air
 REAL,     DIMENSION(:),     intent(in)    :: PDV      ! Diffusivity of water vapor in the air
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRRS     ! Rain water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRIS     ! Pristine ice m.r. source
diff --git a/src/MNH/rain_ice_fast_rh.f90 b/src/MNH/rain_ice_fast_rh.f90
index 3917217e5f169a368526fc37896e903bb3079d16..1710f8b157114ebd6cd2d256fcca23d39d7e3571 100644
--- a/src/MNH/rain_ice_fast_rh.f90
+++ b/src/MNH/rain_ice_fast_rh.f90
@@ -21,7 +21,7 @@ MODULE MODE_RAIN_ICE_FAST_RH
 CONTAINS
 
 SUBROUTINE RAIN_ICE_FAST_RH(OMICRO, PRHODREF, PRVT, PRCT, PRIT, PRST, PRGT, PRHT, PRHODJ, PPRES, &
-                            PZT, PLBDAS, PLBDAG, PLBDAH, PLSFACT, PLVFACT, PCJ, PKA, PDV, PRHODJ3D, PTHS3D, &
+                            PZT, PLBDAS, PLBDAG, PLBDAH, PLSFACT, PLVFACT, PCJ, PKA, PDV, &
                             PRCS, PRRS, PRIS, PRSS, PRGS, PRHS, PTHS, PUSW)
 !
 !*      0. DECLARATIONS
@@ -62,8 +62,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PLVFACT  ! L_v/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PCJ      ! Function to compute the ventilation coefficient
 REAL,     DIMENSION(:),     intent(in)    :: PKA      ! Thermal conductivity of the air
 REAL,     DIMENSION(:),     intent(in)    :: PDV      ! Diffusivity of water vapor in the air
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRRS     ! Rain water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRIS     ! Pristine ice m.r. source
diff --git a/src/MNH/rain_ice_fast_ri.f90 b/src/MNH/rain_ice_fast_ri.f90
index 568c39de3beec373eef77db9f46b4ca53d2bc030..edb36a38b90fded9262d13f47e042247f5b448d5 100644
--- a/src/MNH/rain_ice_fast_ri.f90
+++ b/src/MNH/rain_ice_fast_ri.f90
@@ -19,7 +19,7 @@ MODULE MODE_RAIN_ICE_FAST_RI
 CONTAINS
 
 SUBROUTINE RAIN_ICE_FAST_RI(OMICRO, PRHODREF, PRIT, PRHODJ, PZT, PSSI, PLSFACT, PLVFACT, &
-                            PAI, PCJ, PRHODJ3D, PTHS3D, PCIT, PRCS, PRIS, PTHS)
+                            PAI, PCJ, PCIT, PRCS, PRIS, PTHS)
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -47,8 +47,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PLSFACT  ! L_s/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PLVFACT  ! L_v/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PAI      ! Thermodynamical function
 REAL,     DIMENSION(:),     intent(in)    :: PCJ      ! Function to compute the ventilation coefficient
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
 REAL,     DIMENSION(:),     intent(inout) :: PCIT     ! Pristine ice conc. at t
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRIS     ! Pristine ice m.r. source
diff --git a/src/MNH/rain_ice_fast_rs.f90 b/src/MNH/rain_ice_fast_rs.f90
index d5e2a8e865db312e8cd3b4ec247b76c3e6eb8d88..d5d605ba51366b03faf3883227859e6abfedd63a 100644
--- a/src/MNH/rain_ice_fast_rs.f90
+++ b/src/MNH/rain_ice_fast_rs.f90
@@ -21,7 +21,7 @@ MODULE MODE_RAIN_ICE_FAST_RS
 CONTAINS
 
 SUBROUTINE RAIN_ICE_FAST_RS(PTSTEP, OMICRO, PRHODREF, PRVT, PRCT, PRRT, PRST, PRHODJ, PPRES, PZT, &
-                            PLBDAR, PLBDAS, PLSFACT, PLVFACT, PCJ, PKA, PDV, PRHODJ3D, PTHS3D, &
+                            PLBDAR, PLBDAS, PLSFACT, PLVFACT, PCJ, PKA, PDV, &
                             PRCS, PRRS, PRSS, PRGS, PTHS)
 !
 !*      0. DECLARATIONS
@@ -62,8 +62,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PLVFACT  ! L_v/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PCJ      ! Function to compute the ventilation coefficient
 REAL,     DIMENSION(:),     intent(in)    :: PKA      ! Thermal conductivity of the air
 REAL,     DIMENSION(:),     intent(in)    :: PDV      ! Diffusivity of water vapor in the air
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRRS     ! Rain water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRSS     ! Snow/aggregate m.r. source
diff --git a/src/MNH/rain_ice_slow.f90 b/src/MNH/rain_ice_slow.f90
index 13803b25db289e5276d088bb687d8ecd23052382..8ed0b10ac71174f783ce4c47eecf16e4587fc1ae 100644
--- a/src/MNH/rain_ice_slow.f90
+++ b/src/MNH/rain_ice_slow.f90
@@ -17,11 +17,10 @@ MODULE MODE_RAIN_ICE_SLOW
 
 CONTAINS
 
-SUBROUTINE RAIN_ICE_SLOW(OMICRO, PINVTSTEP, PRHODREF,                      &
-                         PRCT, PRRT, PRIT, PRST, PRGT, PRHODJ, PZT, PPRES, &
-                         PLSFACT, PLVFACT,                                 &
-                         PSSI, PRHODJ3D, PTHS3D, PRVS3D,                   &
-                         PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, PTHS,         &
+SUBROUTINE RAIN_ICE_SLOW(OMICRO, PINVTSTEP, PRHODREF,                                      &
+                         PRCT, PRRT, PRIT, PRST, PRGT, PRHODJ, PZT, PPRES,                 &
+                         PLSFACT, PLVFACT, PSSI,                                           &
+                         PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, PTHS,                         &
                          PAI, PCJ, PKA, PDV, PLBDAS, PLBDAG)
 !
 !*      0. DECLARATIONS
@@ -55,9 +54,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PPRES    ! Pressure
 REAL,     DIMENSION(:),     intent(in)    :: PLSFACT  ! L_s/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PLVFACT  ! L_v/(Pi_ref*C_ph)
 REAL,     DIMENSION(:),     intent(in)    :: PSSI     ! Supersaturation over ice
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRVS3D   ! Water vapor m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRVS     ! Water vapor m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRRS     ! Rain water m.r. source
diff --git a/src/MNH/rain_ice_warm.f90 b/src/MNH/rain_ice_warm.f90
index e030d5ca7dcd74b40433b5f7626355ab3bd440a0..133dc888b8a8c74c6ca2708e58baf37b7587fc1a 100644
--- a/src/MNH/rain_ice_warm.f90
+++ b/src/MNH/rain_ice_warm.f90
@@ -21,7 +21,7 @@ CONTAINS
 SUBROUTINE RAIN_ICE_WARM(OMICRO, KMICRO, K1, K2, K3,                                                           &
                          PRHODREF, PRVT, PRCT, PRRT, PHLC_HCF, PHLC_LCF, PHLC_HRC, PHLC_LRC,                   &
                          PRHODJ, PPRES, PZT, PLBDAR, PLBDAR_RF, PLVFACT, PCJ, PKA, PDV, PRF, PCF, PTHT, PTHLT, &
-                         PRHODJ3D, PTHS3D, PRVS3D, PRVS, PRCS, PRRS, PTHS, PUSW, PEVAP3D)
+                         PRVS, PRCS, PRRS, PTHS, PUSW, PEVAP3D)
 !
 !*      0. DECLARATIONS
 !          ------------
@@ -68,9 +68,6 @@ REAL,     DIMENSION(:),     intent(in)    :: PRF      ! Rain fraction
 REAL,     DIMENSION(:),     intent(in)    :: PCF      ! Cloud fraction
 REAL,     DIMENSION(:),     intent(in)    :: PTHT     ! Potential temperature
 REAL,     DIMENSION(:),     intent(in)    :: PTHLT    ! Liquid potential temperature
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRHODJ3D ! Dry density * Jacobian
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PTHS3D   ! Theta source
-REAL,     DIMENSION(:,:,:), INTENT(IN)    :: PRVS3D   ! Water vapor m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRVS     ! Water vapor m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRCS     ! Cloud water m.r. source
 REAL,     DIMENSION(:),     INTENT(INOUT) :: PRRS     ! Rain water m.r. source
diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 030716ac977d58321ffc122cb863b6baa4eea661..b38be9d78b33525a200a7a18ae524c304ecea68d 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -134,6 +134,7 @@ END MODULE MODI_READ_ALL_DATA_GRIB_CASE
 !  Q. Rodier   16/09/2019: switch of GRIB number ID for orography in ARPEGE/AROME in EPyGrAM
 !  Q. Rodier   27/01/2020: switch of GRIB number ID for orography and hydrometeors in ARPEGE/AROME in EPyGrAM v1.3.7
 !  Q. Rodier   21/04/2020: correction GFS u and v wind component written in the right vertical order
+!  Q. Rodier   02/09/2020 : Read and interpol geopotential height for interpolation on isobaric surface Grid of NCEP 
 !-------------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -206,7 +207,7 @@ INTEGER                            :: ILUOUT0       ! Unit used for output msg.
 INTEGER                            :: IRESP   ! Return code of FM-routines
 INTEGER                            :: IRET          ! Return code from subroutines
 INTEGER(KIND=kindOfInt)            :: IRET_GRIB          ! Return code from subroutines
-INTEGER, PARAMETER                 :: JP_GFS=26     ! number of pressure levels for GFS model
+INTEGER, PARAMETER                 :: JP_GFS=31     ! number of pressure levels for GFS model
 REAL                               :: ZA,ZB,ZC      ! Dummy variables
 REAL                               :: ZD,ZE,ZF      !  |
 REAL                               :: ZTEMP         !  |
@@ -314,6 +315,7 @@ REAL, DIMENSION(:,:), ALLOCATABLE   :: ZPF_G    ! Pressure (flux point)
 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZPM_G    ! Pressure (mass point)
 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZEXNF_G  ! Exner fct. (flux point)
 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZEXNM_G  ! Exner fct. (mass point)
+REAL, DIMENSION(:,:), ALLOCATABLE   :: ZGH_G     ! Geopotential Height
 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZT_G     ! Temperature
 REAL, DIMENSION(:,:), ALLOCATABLE   :: ZQ_G     ! Specific humidity
 REAL, DIMENSION(:), ALLOCATABLE     :: ZH_G     ! Relative humidity
@@ -337,7 +339,7 @@ INTEGER :: IVERSION,ILEVTYPE
 LOGICAL                                       :: GFIND  ! to test if sea wave height is found
 !---------------------------------------------------------------------------------------
 IP_GFS=(/1000,975,950,925,900,850,800,750,700,650,600,550,500,450,400,350,300,&
-           250,200,150,100,70,50,30,20,10/)!
+           250,200,150,100,70,50,30,20,10,7,5,3,2,1/)!
 !
 TZFILE => NULL()
 !
@@ -528,15 +530,10 @@ SELECT CASE (IMODEL)
          END IF
        ENDIF 
   CASE(10) ! NCEP
-      DO IVAR=0,222
-        CALL SEARCH_FIELD(IGRIB,INUM_ZS,KPARAM=IVAR)
+       CALL SEARCH_FIELD(IGRIB,INUM_ZS,KDIS=0,KCAT=3,KNUMBER=5,KTFFS=1)
         IF(INUM_ZS < 0) THEN
-          WRITE (ILUOUT0,'(A)')'Orography is missing'
-        ENDIF
-      END DO
-      INUM_ZS=218
-      WRITE (ILUOUT0,*) 'lsm  ',IGRIB(350)
-      WRITE (ILUOUT0,*) 'orog ',IGRIB(INUM_ZS)     
+          WRITE (ILUOUT0,'(A)')'Orography is missing - abort'
+        ENDIF    
 END SELECT
 ZPARAM(:)=-999.
 CALL GRIB_GET(IGRIB(INUM_ZS),'Nj',INJ,IRET_GRIB)
@@ -739,7 +736,7 @@ IF (IMODEL/=10) THEN
   CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IQ,KLEV1=ISTARTLEVEL)
   IF(INUM < 0) call Print_msg( NVERB_FATAL, 'IO', 'READ_ALL_DATA_GRIB_CASE', 'atmospheric specific humidity is missing' )
 ELSE ! NCEP
-  ISTARTLEVEL=10
+  ISTARTLEVEL=1000
   IT=130
   IQ=157
   CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ISTARTLEVEL)
@@ -786,7 +783,7 @@ ELSE ! NCEP
     END IF
     CALL GRIB_GET(IGRIB(INUM),'values',ZQ_G(:,JLOOP1),IRET_GRIB)
     WRITE (ILUOUT0,*) 'Q ',ILEV1,IRET_GRIB
-    CALL SEARCH_FIELD(IGRIB,INUM,KPARAM=IT,KLEV1=ILEV1)
+    CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=0,KNUMBER=0,KLEV1=ILEV1,KTFFS=100)
     IF (INUM< 0) THEN
       WRITE(YMSG,*) 'atmospheric temperature level ',JLOOP1,' is missing'
       CALL PRINT_MSG(NVERB_FATAL,'IO','READ_ALL_DATA_GRIB_CASE',YMSG)
@@ -1025,6 +1022,43 @@ ELSE !NCEP
 END IF
   
 DEALLOCATE (ZOUT)
+
+
+!---------------------------------------------------------------------------------------
+!* 2.5.4.2 Read and interpol geopotential height for interpolation on isobaric surface Grid of NCEP 
+!---------------------------------------------------------------------------------------
+!
+ALLOCATE (ZGH_G(ISIZE,INLEVEL))
+!
+IF(IMODEL==10) THEN !NCEP with pressure grid only
+ DO JLOOP1=1, INLEVEL
+  ILEV1 = IP_GFS(JLOOP1)
+  WRITE (ILUOUT0,'(A)') ' | Searching geopotential height'
+  CALL SEARCH_FIELD(IGRIB,INUM,KDIS=0,KCAT=3,KNUMBER=5,KLEV1=ILEV1)
+    IF (INUM< 0) THEN
+    !callabortstop
+      WRITE(YMSG,*) 'Geopoential height level ',JLOOP1,' is missing'
+      CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_ALL_DATA_GRIB_CASE',YMSG)
+    END IF
+  !
+  CALL GRIB_GET(IGRIB(INUM),'values',ZGH_G(:,JLOOP1),IRET_GRIB)
+  CALL GRIB_GET(IGRIB(INUM),'Nj',INJ,IRET_GRIB)
+  !
+  END DO
+ !
+ CALL COORDINATE_CONVERSION(IMODEL,IGRIB(INUM_ZS),IIU,IJU,ZLONOUT,ZLATOUT,&
+          ZXOUT,ZYOUT,INI,ZPARAM,IINLO)
+ !
+ ALLOCATE(ZOUT(INO))
+ ALLOCATE(XGH_LS(IIU,IJU,INLEVEL))
+ !
+ DO JLOOP1=1, INLEVEL
+  CALL HORIBL(ZPARAM(3),ZPARAM(4),ZPARAM(5),ZPARAM(6),INT(ZPARAM(2)),IINLO,INI, &
+              ZGH_G(:,JLOOP1),INO,ZXOUT,ZYOUT,ZOUT,.FALSE.,PTIME_HORI,.FALSE.)
+  CALL ARRAY_1D_TO_2D (INO,ZOUT,IIU,IJU,XGH_LS(:,:,JLOOP1))
+ END DO
+ DEALLOCATE(ZOUT)
+END IF
 !
 !*  2.5.5  Compute atmospheric pressure on MESO-NH grid
 !
@@ -1867,7 +1901,7 @@ END SUBROUTINE ARRAY_1D_TO_2D
 !---------------------------------------------------------------------------------------
 !---------------------------------------------------------------------------------------
 !#################################################################################
-SUBROUTINE SEARCH_FIELD(KGRIB,KNUM,KPARAM,KDIS,KCAT,KNUMBER,KLEV1)
+SUBROUTINE SEARCH_FIELD(KGRIB,KNUM,KPARAM,KDIS,KCAT,KNUMBER,KLEV1,KTFFS)
 !#################################################################################
 ! search the grib message corresponding to KPARAM,KLTYPE,KLEV1,KLEV2 in all 
 ! the KGIRB messages
@@ -1885,13 +1919,14 @@ INTEGER,INTENT(IN),OPTIONAL     :: KDIS ! Discipline (GRIB2)
 INTEGER,INTENT(IN),OPTIONAL     :: KCAT ! Catégorie (GRIB2)
 INTEGER,INTENT(IN),OPTIONAL     :: KNUMBER ! parameterNumber (GRIB2)
 INTEGER,INTENT(IN),OPTIONAL     :: KLEV1  ! Level 
+INTEGER,INTENT(IN),OPTIONAL     :: KTFFS  ! TypeOfFirstFixedSurface 
 !
 ! Declaration of local variables
 !
 INTEGER :: IFOUND  ! Number of correct parameters
 INTEGER :: ISEARCH  ! Number of correct parameters to find
 INTEGER :: IRET    ! error code 
-INTEGER :: IPARAM,IDIS,ICAT,INUMBER
+INTEGER :: IPARAM,IDIS,ICAT,INUMBER,ITFFS
 INTEGER :: ILEV1   ! Level parameter 1
 INTEGER :: JLOOP   ! Dummy counter
 INTEGER :: IVERSION
@@ -1909,6 +1944,7 @@ IF (PRESENT(KDIS)) ISEARCH=ISEARCH+1
 IF (PRESENT(KCAT)) ISEARCH=ISEARCH+1
 IF (PRESENT(KNUMBER)) ISEARCH=ISEARCH+1
 IF (PRESENT(KLEV1)) ISEARCH=ISEARCH+1
+IF(PRESENT(KTFFS)) ISEARCH=ISEARCH+1
 !
 DO JLOOP=1,SIZE(KGRIB)
       IFOUND = 0
@@ -1921,6 +1957,23 @@ DO JLOOP=1,SIZE(KGRIB)
         WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
         CYCLE
       ENDIF
+      !
+     IF (PRESENT(KTFFS)) THEN
+        CALL GRIB_GET(KGRIB(JLOOP),'typeOfFirstFixedSurface',ITFFS,IRET_GRIB)
+        IF (IRET_GRIB >   0) THEN
+          WRITE (ILUOUT0,'(A)')' | Error encountered in the Grib file, skipping field'
+          CYCLE
+        ELSE IF (IRET_GRIB == -6) THEN
+          WRITE (ILUOUT0,'(A)')' | ECMWF pseudo-Grib data encountered, skipping field'
+          CYCLE
+        ENDIF
+        IF (ITFFS==KTFFS) THEN
+          IFOUND = IFOUND + 1
+        ELSE
+          CYCLE
+        ENDIF
+      ENDIF
+
       IF (PRESENT(KPARAM)) THEN
         IF (IVERSION == 2) THEN
           CALL GRIB_GET(KGRIB(JLOOP),'paramId',IPARAM,IRET_GRIB)
diff --git a/src/MNH/sources_neg_correct.f90 b/src/MNH/sources_neg_correct.f90
index 4844e2e1e4b4733c20611bfcf1ac03d8bab67c7c..72c00c721e0529b0350fd2710eeca8ad65956456 100644
--- a/src/MNH/sources_neg_correct.f90
+++ b/src/MNH/sources_neg_correct.f90
@@ -6,6 +6,7 @@
 ! Author: P. Wautelet 25/06/2020 (deduplication of code from advection_metsv, resolved_cloud and turb)
 ! Modifications:
 !  P. Wautelet 30/06/2020: remove non-local corrections in resolved_cloud for NEGA => new local corrections here
+!  J. Escobar  21/07/2020: bug <-> array of size(:,:,:,0) => return if krr=0
 !-----------------------------------------------------------------
 module mode_sources_neg_correct
 
@@ -52,6 +53,8 @@ integer :: jrmax
 integer :: jsv
 real, dimension(:, :, :), allocatable :: zt, zexn, zlv, zls, zcph, zcor
 
+if ( krr == 0 ) return
+
 if ( hbudname /= 'NEADV' .and. hbudname /= 'NECON' .and. hbudname /= 'NEGA' .and. hbudname /= 'NETUR' ) &
   call Print_msg( NVERB_WARNING, 'GEN', 'Sources_neg_correct', 'budget '//hbudname//' not yet tested' )
 
diff --git a/src/MNH/ver_prep_gribex_case.f90 b/src/MNH/ver_prep_gribex_case.f90
index 1c6ea2b2af6033a4dc355ca8f44fb22be3ae1cb5..8a6cc54bb401f84f1c5e5cbeafff1bb41cb6e2e9 100644
--- a/src/MNH/ver_prep_gribex_case.f90
+++ b/src/MNH/ver_prep_gribex_case.f90
@@ -85,6 +85,8 @@ END MODULE MODI_VER_PREP_GRIBEX_CASE
 !!                  May 2006                 Remove EPS
 !!                  Apr, 09 2018 (J.-P. Chaboureau) add isobaric surface 
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!!                  Sep, 02, 2020 (Q. Rodier) use of geopotential height instead of
+!!                                height above orography for isobaric surface interpolation
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -271,8 +273,9 @@ IF (HFILE(1:3)=='ATM') THEN
                                 ZU_LS,ZV_LS,                     &
                                 ZW_LS,'FLUX'                     )
   ELSE                      ! isobaric surfaces (w at mass points)
+    !Warning, in that case (NCEP only for now) ZZFLUX_LS is not correct (but not used)
     CALL VER_INTERP_TO_MIXED_GRID('ATM ',.TRUE.,XZS_LS,XZSMT_LS,    &
-                                  ZZMASS_LS,ZSV_LS,                &
+                                  XGH_LS,ZSV_LS, &
                                   ZZFLUX_LS,XPS_LS,ZPMHP_LS,       &
                                   ZTHV_LS,ZR_LS,                   &
                                   ZHU_LS,                          &
diff --git a/src/MNH/write_lesn.f90 b/src/MNH/write_lesn.f90
index 51cb134222ec969fe0aa83e17b59d62a9e1719df..09724ca6a50988a652edf5f748576c19a04c3c99 100644
--- a/src/MNH/write_lesn.f90
+++ b/src/MNH/write_lesn.f90
@@ -66,6 +66,8 @@ subroutine  Write_les_n( tpdiafile )
 !  C. Lac         02/2019: add rain fraction as a LES diagnostic
 !  P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
 !  P. Wautelet 12/10/2020: remove HLES_AVG dummy argument and group all 4 calls
+!  P. Wautelet 13/10/2020: bugfix: correct some names for LES_DIACHRO_2PT diagnostics (Ri)
+!  P. Wautelet 26/10/2020: bugfix: correct some comments and conditions + add missing RES_RTPZ
 !  P. Wautelet 26/10/2020: restructure subroutines to use tfield_metadata_base type
 ! --------------------------------------------------------------------------
 !
diff --git a/src/SURFEX/mode_read_grib.F90 b/src/SURFEX/mode_read_grib.F90
index 344ea5967b0658f1a362b4eef77366bc302691fa..2751bff06b17107df33231e213da599818c03176 100644
--- a/src/SURFEX/mode_read_grib.F90
+++ b/src/SURFEX/mode_read_grib.F90
@@ -63,11 +63,9 @@ IF (NGRIB_VERSION==1) THEN
  CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'indicatorOfParameter',IRET)
 ELSEIF (NGRIB_VERSION==2) THEN
  IF(HINMODEL=='ARPEGE') THEN
-   print*,"CALL GRIB_INDEX_CREATE"
    CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'parameterNumber',IRET)
-print*,IRET
  ELSE
-   CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'paramId',IRET)        
+   CALL GRIB_INDEX_CREATE(NIDX,HGRIB,'paramId',IRET)
  ENDIF
 ENDIF
 IF (IRET/=0) CALL ABOR1_SFX("MODE_READ_GRIB:MAKE_GRIB_INDEX: error while creating the grib index")
@@ -136,20 +134,16 @@ IRET = 0
 KFOUND=0
 !
 DO WHILE (IRET /= GRIB_END_OF_INDEX .AND. KFOUND/=3)
-print*,"===============new message=============="
   !
   IRET = 0
   KFOUND=0
   !
   IF (KLTYPE/=-2) THEN
     CALL GRIB_GET(KGRIB,'indicatorOfTypeOfLevel',ILTYPE,IRET)
-print*,IRET
     IF(IRET/=0) THEN
       CALL GRIB_GET(KGRIB,'typeOfFirstFixedSurface',ILTYPE,IRET)
     ENDIF
-print*,IRET
     CALL TEST_IRET(KLUOUT,ILTYPE,KLTYPE,IRET)
-    print*,"ILTYPE,KLTYPE,IRET",ILTYPE,KLTYPE,IRET
   ELSE
      IF (PRESENT(HTYPELEVEL)) THEN
         CALL GRIB_GET(KGRIB,'typeOfLevel',YTYPELEVEL,IRET)
@@ -170,7 +164,6 @@ print*,IRET
 
     ENDIF
     !
-print*,KFOUND
     IF (IRET.EQ.0) THEN
       !
       KFOUND = KFOUND + 1
@@ -185,7 +178,6 @@ print*,KFOUND
       !
       IF (IRET.EQ.0) KFOUND = KFOUND + 1
       !
-print*,KFOUND
     ENDIF
     !
   ENDIF
@@ -348,11 +340,9 @@ IF (NGRIB_VERSION == 1) THEN
  CALL GRIB_INDEX_SELECT(NIDX,'indicatorOfParameter',KPARAM,KRET)
 ELSEIF (NGRIB_VERSION == 2) THEN
  IF (HINMODEL=='ARPEGE') THEN
-  print*,"GRIB_INDEX_SELECT :",KPARAM
    CALL GRIB_INDEX_SELECT(NIDX,'parameterNumber',KPARAM,KRET)
-print*,KRET
  ELSE
-   CALL GRIB_INDEX_SELECT(NIDX,'paramId',KPARAM,KRET)         
+   CALL GRIB_INDEX_SELECT(NIDX,'paramId',KPARAM,KRET)
  ENDIF
 END IF
  CALL GRIB_NEW_FROM_INDEX(NIDX,IGRIB,KRET)
@@ -363,7 +353,6 @@ IF (KRET.EQ.0) THEN
    IF (PRESENT(HTYPELEVEL)) THEN
       CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND,HTYPELEVEL,ZLEV1,ZLEV2)
    ELSE
-print*,"CALL GET_GRIB_MESSAGE"
       CALL GET_GRIB_MESSAGE(KLUOUT,ILTYPE,ILEV1,ILEV2,IGRIB,IFOUND)
    ENDIF
 ENDIF
@@ -390,8 +379,6 @@ IF (PRESENT(KPARAM2)) THEN
   ENDIF
 ENDIF
 !
-print*,"IFOUND=",IFOUND
-
 IF (IFOUND==3) THEN
   !
   IF (PRESENT(KLTYPE)) KLTYPE = ILTYPE
@@ -443,8 +430,6 @@ INTEGER :: INUM_ZS,ISIZE,ICOUNT,JLOOP,IPARAM,IGRIB,IPARAM2
 CHARACTER(LEN=24) :: YLTYPELU
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !-------------------------------------------------------------------
-print*,"HINMODEL=",HINMODEL
-print*,"NGRIB_VERSION=",NGRIB_VERSION
 IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_LAND_MASK',0,ZHOOK_HANDLE)
 WRITE (KLUOUT,'(A)') 'MODE_READ_GRIB:READ_GRIB_LAND_MASK: | Reading land mask from ',HINMODEL
 !
@@ -454,14 +439,12 @@ SELECT CASE (HINMODEL)
   CASE ('NCEP  ')
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,172,IRET,PMASK) 
   CASE ('ARPEGE','ALADIN','MOCAGE')
-    print*,"NGRIB_VERSION=",NGRIB_VERSION
     IF(HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
       ILTYPE=1
-      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,0,IRET,PMASK,KLTYPE=ILTYPE)   
+      CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,0,IRET,PMASK,KLTYPE=ILTYPE)
     ELSE
       CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,81,IRET,PMASK)          
     ENDIF
-print*,"NGRIB_VERSION=",NGRIB_VERSION
   CASE ('HIRLAM')        
     ILTYPE=105
     ILEV  =0    
@@ -514,7 +497,7 @@ SELECT CASE (HINMODEL)
   CASE ('ECMWF ')
     CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,129,IRET,PZS) 
   CASE ('NCEP  ')
-    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228002,IRET,PZS)               
+    CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,228002,IRET,PZS)    
   CASE ('ARPEGE','MOCAGE')
     IF (HINMODEL=='ARPEGE' .AND. NGRIB_VERSION==2) THEN
       CALL READ_GRIB(HGRIB,HINMODEL,KLUOUT,4,IRET,PZS)
@@ -536,7 +519,9 @@ IF (IRET /= 0) THEN
 END IF
 !
 ! Datas given in archives are multiplied by the gravity acceleration
-PZS = PZS / XG
+IF(HINMODEL /= 'NCEP') THEN
+  PZS = PZS / XG
+END IF
 !
 IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_ZS',1,ZHOOK_HANDLE)
 END SUBROUTINE READ_GRIB_ZS
@@ -756,7 +741,7 @@ SELECT CASE (HINMODEL)
      CALL ABOR1_SFX('MODE_READ_GRIB:READ_GRIB_TSWATER:OPTION NOT SUPPORTED '//HINMODEL)    
 END SELECT
 !
-IF (SIZE(PMASK)==SIZE(PTS)) WHERE ((PMASK(:)/=1.) .OR. ((PMASK(:)==1.) .AND.(PTS(:)==9999.))) PTS = XUNDEF
+IF (SIZE(PMASK)==SIZE(PTS)) WHERE ((PMASK(:)/=0.) .OR. ((PMASK(:)==0.) .AND.(PTS(:)==9999.))) PTS = XUNDEF
 
 !
 IF (LHOOK) CALL DR_HOOK('MODE_READ_GRIB:READ_GRIB_TSWATER',1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/pgd_fieldin.F90 b/src/SURFEX/pgd_fieldin.F90
index 683611a1c6db9abc6a56ac5ce4510edf25811734..990727e151a2f567c1181f2701d293f124f10186 100644
--- a/src/SURFEX/pgd_fieldin.F90
+++ b/src/SURFEX/pgd_fieldin.F90
@@ -57,7 +57,7 @@ USE MODD_PGDWORK,        ONLY : XALL, NSIZE_ALL, CATYPE, NSIZE, XSUMVAL,   &
 USE MODD_SURF_PAR,       ONLY : XUNDEF
 USE MODD_PGD_GRID,       ONLY : NL
 !
-USE MODD_DATA_COVER_PAR, ONLY : NTYPE, LVEG_PRES, NVEGTYPE
+USE MODD_DATA_COVER_PAR, ONLY : NTYPE, LVEG_PRES, NVEGTYPE, NVEGTYPE_OLD
 !
 USE MODI_GET_LUOUT
 USE MODI_TREAT_FIELD
@@ -293,8 +293,8 @@ IF (LEN_TRIM(HFILE)/=0) THEN
 !
   DO JT=1,SIZE(NSIZE,2)
 
-    IF (.NOT.U%LECOSG.AND.JT>NVEGTYPE) EXIT
-    
+    IF (.NOT.U%LECOSG.AND.JT>NVEGTYPE_OLD) EXIT
+
     !multitype input file
     IF (SIZE(ZFIELD,2)>1) THEN
 
diff --git a/src/SURFEX/prep_isba_ascllv.F90 b/src/SURFEX/prep_isba_ascllv.F90
index 3b9ba798d4992d3a46205ac5e198381b7d45d8bd..561e33e051493f3eb1f668839b92efdbeed09996 100644
--- a/src/SURFEX/prep_isba_ascllv.F90
+++ b/src/SURFEX/prep_isba_ascllv.F90
@@ -91,6 +91,11 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
 IF (LHOOK) CALL DR_HOOK('PREP_ISBA_ASCLLV',0,ZHOOK_HANDLE)
 !
+IF ((.NOT.ALLOCATED(NINDEX)).AND.(HPROGRAM=='MESONH')) THEN 
+  ALLOCATE(NINDEX(U%NDIM_FULL))
+  NINDEX(:) = 0
+ENDIF
+!
 IF (.NOT.ALLOCATED(NNUM)) THEN
   ALLOCATE(NNUM(U%NDIM_FULL))
   IF (NRANK/=NPIO) THEN
diff --git a/src/SURFEX/write_bld_descriptionn.F90 b/src/SURFEX/write_bld_descriptionn.F90
index 8f71b236d9df9051db6512036c9c86664006dae1..8ab77a32c78345b28b901477e2d50e32fe13e2a1 100644
--- a/src/SURFEX/write_bld_descriptionn.F90
+++ b/src/SURFEX/write_bld_descriptionn.F90
@@ -70,7 +70,6 @@ INTEGER                         :: IRESP
 INTEGER                         :: I1, I2
 INTEGER                         :: JL
 INTEGER                         :: ITOT
- CHARACTER(LEN=LEN_HREC)         :: YRECFM
  CHARACTER(LEN=100)              :: YCOMMENT
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
@@ -92,9 +91,8 @@ ZWORK(5) = FLOAT(BDD%NDESC_ROOF_LAYER)
 ZWORK(6) = FLOAT(BDD%NDESC_ROAD_LAYER)
 ZWORK(7) = FLOAT(BDD%NDESC_FLOOR_LAYER)
 !
-YRECFM='Bld_dimensions  '
 YCOMMENT='Configuration numbers for descriptive building data'
- CALL WRITE_SURF(HSELECT, HPROGRAM,'BLD_DESC_CNF',ZWORK,IRESP,YCOMMENT,'-',YRECFM)
+ CALL WRITE_SURF(HSELECT, HPROGRAM,'BLD_DESC_CNF',ZWORK,IRESP,YCOMMENT,'-','Bld_dimensions  ')
 DEALLOCATE(ZWORK)
 !
 !-------------------------------------------------------------------------------
@@ -183,9 +181,8 @@ END DO
  CALL UP_DESC_IND_W(BDD%NDESC_AGE) ; ZWORK(I1:I2) = FLOAT(BDD%NDESC_AGE_DATE(:))
 !
 YCOMMENT='Descriptive building data'
-YRECFM='Bld_parameters  '
  CALL WRITE_SURF(HSELECT, &
-                 HPROGRAM,'BLD_DESC_DAT',ZWORK,IRESP,YCOMMENT,'-',YRECFM)
+                 HPROGRAM,'BLD_DESC_DAT',ZWORK,IRESP,YCOMMENT,'-','Bld_parameters  ')
 DEALLOCATE(ZWORK)
 !
 IF (LHOOK) CALL DR_HOOK('WRITE_BLD_DESCRIPTION_n',1,ZHOOK_HANDLE)