diff --git a/src/MNH/boundaries.f90 b/src/MNH/boundaries.f90
index e40949d5e7a522683d6df30670dd4a28c6a8703d..46fa4bc36998240959a92fdc7de2c4bf6018c38c 100644
--- a/src/MNH/boundaries.f90
+++ b/src/MNH/boundaries.f90
@@ -880,14 +880,14 @@ IF (CCLOUD == 'LIMA' .AND. IMI == 1) THEN
   DO JSV=NSV_LIMA_CCN_FREE,NSV_LIMA_CCN_FREE+NMOD_CCN-1 ! LBC for CCN from MACC
     IF (GLIMABOUNDARY(JSV-NSV_LIMA_CCN_FREE+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO 
   DO JSV=NSV_LIMA_IFN_FREE,NSV_LIMA_IFN_FREE+NMOD_IFN-1 ! LBC for IFN from MACC
     IF (GLIMABOUNDARY(JSV-NSV_LIMA_IFN_FREE+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO 
@@ -910,7 +910,7 @@ IF (LUSECHEM .AND. IMI == 1) THEN
   DO JSV=NSV_CHEMBEG,NSV_CHEMEND
     IF (GCHBOUNDARY(JSV-NSV_CHEMBEG+1))  THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -933,7 +933,7 @@ IF (LUSECHIC .AND. IMI == 1) THEN
   DO JSV=NSV_CHICBEG,NSV_CHICEND
     IF (GICBOUNDARY(JSV-NSV_CHICBEG+1))  THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -955,7 +955,7 @@ IF (LORILAM .AND. IMI == 1) THEN
   DO JSV=NSV_AERBEG,NSV_AEREND
     IF (GAERBOUNDARY(JSV-NSV_AERBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -978,7 +978,7 @@ IF (LDUST .AND. IMI == 1) THEN
   DO JSV=NSV_DSTBEG,NSV_DSTEND
     IF (GDSTBOUNDARY(JSV-NSV_DSTBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO 
@@ -1001,7 +1001,7 @@ IF (LSALT .AND. IMI == 1) THEN
   DO JSV=NSV_SLTBEG,NSV_SLTEND
     IF (GSLTBOUNDARY(JSV-NSV_SLTBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -1024,7 +1024,7 @@ IF ( LPASPOL .AND. IMI == 1) THEN
   DO JSV=NSV_PPBEG,NSV_PPEND
     IF (GPPBOUNDARY(JSV-NSV_PPBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-       CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -1047,7 +1047,7 @@ IF ( LCONDSAMP .AND. IMI == 1) THEN
   DO JSV=NSV_CSBEG,NSV_CSEND
     IF (GCSBOUNDARY(JSV-NSV_CSBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-       CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
@@ -1071,9 +1071,7 @@ IF ( LFOREFIRE .AND. IMI == 1) THEN
   DO JSV=NSV_FFBEG,NSV_FFEND
     IF (GFFBOUNDARY(JSV-NSV_FFBEG+1)) THEN
       IF (SIZE(PSVT)>0) THEN
-       CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV))
-      ELSE
-!!$        CALL CH_BOUNDARIES (HLBCX,HLBCY,PSVM(:,:,:,JSV),PUM,PVM,PSVM(:,:,:,JSV))
+        CALL CH_BOUNDARIES (HLBCX,HLBCY,PUT,PVT,PSVT(:,:,:,JSV),XSVMIN(JSV))
       ENDIF
     ENDIF
   ENDDO
diff --git a/src/MNH/ch_boundaries.f90 b/src/MNH/ch_boundaries.f90
index 0f34c9577a154d8c645b31782fa2abb49da48d1f..33fdb27af5a9fa3e23dee73a1fce12e6261b4573 100644
--- a/src/MNH/ch_boundaries.f90
+++ b/src/MNH/ch_boundaries.f90
@@ -15,11 +15,12 @@ MODULE MODI_CH_BOUNDARIES
 INTERFACE
 !
       SUBROUTINE CH_BOUNDARIES (HLBCX,HLBCY,                 &
-                                PUT,PVT,PSVBT          )
+                                PUT,PVT,PSVBT,PSVMIN         )
 !
 CHARACTER(LEN=4), DIMENSION(2), INTENT(IN)  :: HLBCX,HLBCY  
 REAL, DIMENSION(:,:,:),       INTENT(INOUT) :: PSVBT
 REAL, DIMENSION(:,:,:),         INTENT(IN) :: PUT,PVT
+REAL                                        :: PSVMIN
 !
 END SUBROUTINE CH_BOUNDARIES
 !
@@ -30,7 +31,7 @@ END MODULE MODI_CH_BOUNDARIES
 !
 !     ####################################################################
 SUBROUTINE CH_BOUNDARIES (HLBCX,HLBCY,                           &
-                          PUT,PVT,PSVBT                    )
+                          PUT,PVT,PSVBT,PSVMIN                   )
 !     ####################################################################
 !
 !!****  *CH_BOUNDARIES* - routine to prepare the Lateral Boundary Conditions for
@@ -69,6 +70,7 @@ SUBROUTINE CH_BOUNDARIES (HLBCX,HLBCY,                           &
 !!     Original        06/06/00
 !!     06/06/00 (C. Mari) embedded into mesonh routines
 !!     15/02/01 (P. Tulet) update for MOCAGE lateral boundary conditions
+!!     10/02/17 (M. Leriche) prevent negative values
 !!	
 !-------------------------------------------------------------------------------
 !
@@ -90,6 +92,7 @@ IMPLICIT NONE
 CHARACTER(LEN=4), DIMENSION(2), INTENT(IN)  :: HLBCX,HLBCY  
 REAL, DIMENSION(:,:,:),     INTENT(INOUT)   :: PSVBT
 REAL, DIMENSION(:,:,:),          INTENT(IN) :: PUT,PVT
+REAL                                        :: PSVMIN
 !
 !
 !
@@ -170,12 +173,12 @@ IF (LWEST_ll( ) .AND. HLBCX(1)=='OPEN') THEN
     DO IK= 1,IKU
       IF ( PUT(IIB,IJ,IK) <= 0. ) THEN         !  OUTFLOW condition
         PSVBT(IIB-1,IJ,IK) = &
-             2.*PSVBT (IIB,IJ,IK) -PSVBT (IIB+1,IJ,IK)
+             MAX(PSVMIN,2.*PSVBT (IIB,IJ,IK) -PSVBT (IIB+1,IJ,IK))
       ELSE                                     !  INFLOW  condition
         IF (ZZZW(1,IJ,IK) > -999.) THEN
-          PSVBT(IIB-1,IJ,IK) = PSVBT(IIB+1,IJ,IZZW(1,IJ,IK))+&
+          PSVBT(IIB-1,IJ,IK) = MAX(PSVMIN,PSVBT(IIB+1,IJ,IZZW(1,IJ,IK))+&
                (PSVBT(IIB+1,IJ,IZZW(1,IJ,IK)+1)-&
-               PSVBT(IIB+1,IJ,IZZW(1,IJ,IK))) * ZZZW(1,IJ,IK)
+               PSVBT(IIB+1,IJ,IZZW(1,IJ,IK))) * ZZZW(1,IJ,IK))
           !
         ELSE
           PSVBT(IIB-1,IJ,IK) =  PSVBT(IIB+1,IJ,IK)
@@ -218,13 +221,13 @@ IF (LEAST_ll( ) .AND. HLBCX(1)=='OPEN') THEN
     DO IK=1,IKU
       IF ( PUT(IIE+1,IJ,IK) >= 0. ) THEN         !  OUTFLOW condition
         PSVBT(IIE+1,IJ,IK) = &
-             2.*PSVBT (IIE,IJ,IK) -PSVBT (IIE-1,IJ,IK)
+             MAX(PSVMIN,2.*PSVBT (IIE,IJ,IK) -PSVBT (IIE-1,IJ,IK))
       ELSE                                       !  INFLOW  condition
         IF (ZZZE(1,IJ,IK) > -999.) THEN
           PSVBT(IIE+1,IJ,IK) = &
-               PSVBT(IIE-1,IJ,IZZE(1,IJ,IK))+&
+               MAX(PSVMIN,PSVBT(IIE-1,IJ,IZZE(1,IJ,IK))+&
                (PSVBT(IIE-1,IJ,IZZE(1,IJ,IK)+1)-&
-               PSVBT(IIE-1,IJ,IZZE(1,IJ,IK))) * ZZZE(1,IJ,IK)
+               PSVBT(IIE-1,IJ,IZZE(1,IJ,IK))) * ZZZE(1,IJ,IK))
           !
         ELSE
           PSVBT(IIE+1,IJ,IK) =  PSVBT(IIE-1,IJ,IK)
@@ -268,13 +271,13 @@ IF (LSOUTH_ll( ) .AND. HLBCY(1)=='OPEN') THEN
     DO IK=1,IKU
       IF ( PVT(II,IJB,IK) <= 0. ) THEN         !  OUTFLOW condition
         PSVBT(II,IJB-1,IK) = &
-             2.*PSVBT (II,IJB,IK) -PSVBT (II,IJB+1,IK)
+             MAX(PSVMIN,2.*PSVBT (II,IJB,IK) -PSVBT (II,IJB+1,IK))
       ELSE                                     !  INFLOW  condition
         IF (ZZZS(II,1,IK) > -999.) THEN
           PSVBT(II,IJB-1,IK) = &
-               PSVBT(II,IJB+1,IZZS(II,1,IK))+&
+               MAX(PSVMIN,PSVBT(II,IJB+1,IZZS(II,1,IK))+&
                (PSVBT(II,IJB+1,IZZS(II,1,IK)+1)-&
-               PSVBT(II,IJB+1,IZZS(II,1,IK))) * ZZZS(II,1,IK)
+               PSVBT(II,IJB+1,IZZS(II,1,IK))) * ZZZS(II,1,IK))
           !
         ELSE
           PSVBT(II,IJB-1,IK) =  PSVBT(II,IJB+1,IK)
@@ -319,13 +322,13 @@ IF (LNORTH_ll( ) .AND. HLBCY(2)=='OPEN') THEN
     DO IK=1,IKU
       IF ( PVT(II,IJE+1,IK) >= 0. ) THEN         !  OUTFLOW condition
         PSVBT(II,IJE+1,IK) = &
-             2.*PSVBT (II,IJE,IK) -PSVBT (II,IJE-1,IK)
+             MAX(PSVMIN,2.*PSVBT (II,IJE,IK) -PSVBT (II,IJE-1,IK))
       ELSE                                       !  INFLOW  condition
         IF (ZZZN(II,1,IK) > -999.) THEN
           PSVBT(II,IJE+1,IK) = &
-               PSVBT(II,IJE-1,IZZN(II,1,IK))+&
+               MAX(PSVMIN,PSVBT(II,IJE-1,IZZN(II,1,IK))+&
                (PSVBT(II,IJE-1,IZZN(II,1,IK)+1)-&
-               PSVBT(II,IJE-1,IZZN(II,1,IK))) * ZZZN(II,1,IK)
+               PSVBT(II,IJE-1,IZZN(II,1,IK))) * ZZZN(II,1,IK))
           !
         ELSE
           PSVBT(II,IJE+1,IK) =  PSVBT(II,IJE-1,IK)
diff --git a/src/MNH/ini_mean_field.f90 b/src/MNH/ini_mean_field.f90
index 2a5e5eca1b8214deca9eb357c0cb7b3f8c848c12..36eafb4586599ef980f804c73ac1674cea7448c8 100644
--- a/src/MNH/ini_mean_field.f90
+++ b/src/MNH/ini_mean_field.f90
@@ -47,9 +47,7 @@ END MODULE MODI_INI_MEAN_FIELD
 !!    MODIFICATIONS
 !!    -------------
 !!      Original        11/12/09
-!!      Modifications   10/2016 (C.Lac) Add max values
-!!                      04/2017 (P. Wautelet) Initialize MAX variables to lowest possible value
-!!
+!!                      10/2016 (C.Lac) Add max values
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -62,10 +60,6 @@ USE MODD_PARAM_n
 
 IMPLICIT NONE
 !
-REAL :: ZMIN !Largest real negative value
-!
-ZMIN = -HUGE(ZMIN)
-!
 MEAN_COUNT = 0
 
 XUM_MEAN  = 0.0
@@ -83,12 +77,12 @@ XTH2_MEAN = 0.0
 XTEMP2_MEAN = 0.0
 XPABS2_MEAN = 0.0
 
-XUM_MAX  = ZMIN
-XVM_MAX  = ZMIN
-XWM_MAX  = ZMIN
-XTHM_MAX = ZMIN
-XTEMPM_MAX = ZMIN
-IF (CTURB /= 'NONE') XTKEM_MAX = ZMIN
-XPABSM_MAX = ZMIN
+XUM_MAX  = -1.E20
+XVM_MAX  = -1.E20
+XWM_MAX  = -1.E20
+XTHM_MAX = 0.0
+XTEMPM_MAX = 0.0
+IF (CTURB /= 'NONE') XTKEM_MAX = 0.0
+XPABSM_MAX = 0.0
 
 END SUBROUTINE INI_MEAN_FIELD
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index c4bed29b516a198818ea5505a892cbbfbdddca5b..46cb888c3842ae85a202508353934ff5e893a48b 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -270,9 +270,9 @@ END MODULE MODI_INI_MODEL_n
 !!                   M.Leriche 2016 Chemistry
 !!                   10/2016 M.Mazoyer New KHKO output fields
 !!                      10/2016 (C.Lac) Add max values
-!!       F. Brosse   Oct.  2016 add prod/loss terms computation for chemistry
-!!                   01/2017 (G.Delautier) bug chemistry : modify test for prod/loss terms computation
-!!                   Apr. 2017 (P. Wautelet) allocate MAX variables if LMEAN_FIELD and call INI_MEAN_FIELD
+!!       F. Brosse   Oct.  2016 add prod/loss terms computation for chemistry       
+!!                   M.Leriche 2016 Chemistry
+!!                   M.Leriche 10/02/17 prevent negative values in LBX(Y)SVS 
 !---------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -686,31 +686,34 @@ IF (LMEAN_FIELD) THEN
 !
   MEAN_COUNT = 0
 !
-  ALLOCATE(XUM_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XVM_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XWM_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XTHM_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XTEMPM_MEAN(IIU,IJU,IKU))
-  IF (CTURB/='NONE') ALLOCATE(XTKEM_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XPABSM_MEAN(IIU,IJU,IKU))
-!
-  ALLOCATE(XU2_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XV2_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XW2_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XTH2_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XTEMP2_MEAN(IIU,IJU,IKU))
-  ALLOCATE(XPABS2_MEAN(IIU,IJU,IKU))
-!
-  ALLOCATE(XUM_MAX(IIU,IJU,IKU))
-  ALLOCATE(XVM_MAX(IIU,IJU,IKU))
-  ALLOCATE(XWM_MAX(IIU,IJU,IKU))
-  ALLOCATE(XTHM_MAX(IIU,IJU,IKU))
-  ALLOCATE(XTEMPM_MAX(IIU,IJU,IKU))
-  IF (CTURB/='NONE') ALLOCATE(XTKEM_MAX(IIU,IJU,IKU))
-  ALLOCATE(XPABSM_MAX(IIU,IJU,IKU))
-!
-  CALL INI_MEAN_FIELD()
-!
+  ALLOCATE(XUM_MEAN(IIU,IJU,IKU))      ; XUM_MEAN  = 0.0
+  ALLOCATE(XVM_MEAN(IIU,IJU,IKU))      ; XVM_MEAN  = 0.0
+  ALLOCATE(XWM_MEAN(IIU,IJU,IKU))      ; XWM_MEAN  = 0.0
+  ALLOCATE(XTHM_MEAN(IIU,IJU,IKU))     ; XTHM_MEAN = 0.0
+  ALLOCATE(XTEMPM_MEAN(IIU,IJU,IKU))   ; XTEMPM_MEAN = 0.0
+  IF (CTURB/='NONE') THEN
+     ALLOCATE(XTKEM_MEAN(IIU,IJU,IKU))    
+     XTKEM_MEAN = 0.0
+  END IF
+  ALLOCATE(XPABSM_MEAN(IIU,IJU,IKU))   ; XPABSM_MEAN = 0.0
+!
+  ALLOCATE(XU2_MEAN(IIU,IJU,IKU))      ; XU2_MEAN  = 0.0
+  ALLOCATE(XV2_MEAN(IIU,IJU,IKU))      ; XV2_MEAN  = 0.0
+  ALLOCATE(XW2_MEAN(IIU,IJU,IKU))      ; XW2_MEAN  = 0.0
+  ALLOCATE(XTH2_MEAN(IIU,IJU,IKU))     ; XTH2_MEAN = 0.0
+  ALLOCATE(XTEMP2_MEAN(IIU,IJU,IKU))   ; XTEMP2_MEAN = 0.0
+  ALLOCATE(XPABS2_MEAN(IIU,IJU,IKU))   ; XPABS2_MEAN = 0.0
+!
+  ALLOCATE(XUM_MAX(IIU,IJU,IKU))      ; XUM_MAX  = -1.E20
+  ALLOCATE(XVM_MAX(IIU,IJU,IKU))      ; XVM_MAX  = -1.E20
+  ALLOCATE(XWM_MAX(IIU,IJU,IKU))      ; XWM_MAX  = -1.E20
+  ALLOCATE(XTHM_MAX(IIU,IJU,IKU))     ; XTHM_MAX = 0.0
+  ALLOCATE(XTEMPM_MAX(IIU,IJU,IKU))   ; XTEMPM_MAX = 0.0
+  IF (CTURB/='NONE') THEN
+     ALLOCATE(XTKEM_MAX(IIU,IJU,IKU))    
+     XTKEM_MAX = 0.0
+  END IF
+  ALLOCATE(XPABSM_MAX(IIU,IJU,IKU))   ; XPABSM_MAX = 0.0
 END IF
 !
 IF ((CUVW_ADV_SCHEME(1:3)=='CEN') .AND. (CTEMP_SCHEME == 'LEFR') ) THEN
@@ -1708,6 +1711,59 @@ IF ((KMI==1).AND.(.NOT. LSTEADYLS)) THEN
                XLBXUS,XLBXVS,XLBXWS,XLBXTHS,XLBXTKES,XLBXRS,XLBXSVS,          &
                XLBYUS,XLBYVS,XLBYWS,XLBYTHS,XLBYTKES,XLBYRS,XLBYSVS           )
   CALL MPPDB_CHECK3D(XUT,"INI_MODEL_N-after ini_cpl::XUT",PRECISION)
+!
+      DO JSV=NSV_CHEMBEG,NSV_CHEMEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_LNOXBEG,NSV_LNOXEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_AERBEG,NSV_AEREND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_DSTBEG,NSV_DSTEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_DSTDEPBEG,NSV_DSTDEPEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_SLTBEG,NSV_SLTEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_SLTDEPBEG,NSV_SLTDEPEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+      DO JSV=NSV_PPBEG,NSV_PPEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+#ifdef MNH_FOREFIRE
+      DO JSV=NSV_FFBEG,NSV_FFEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+      !
+#endif
+      DO JSV=NSV_CSBEG,NSV_CSEND
+        XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+        XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+      ENDDO
+!
 END IF
 !
 IF ( KMI > 1) THEN
@@ -2183,11 +2239,12 @@ IF ( LFOREFIRE ) THEN
 		, TDTCUR%TDATE%YEAR, TDTCUR%TDATE%MONTH, TDTCUR%TDATE%DAY, TDTCUR%TIME, XTSTEP)
 END IF
 #endif
+
 !-------------------------------------------------------------------------------
 !
 !*      30.   Total production/Loss for chemical species
 !
-IF (LCHEMDIAG)  THEN
+IF (LUSECHEM.OR.LCHEMDIAG)  THEN
         CALL CH_INIT_PRODLOSSTOT_n(ILUOUT)
         IF (NEQ_PLT>0) THEN
                 ALLOCATE(XPROD(IIU,IJU,IKU,NEQ_PLT))
@@ -2207,7 +2264,7 @@ END IF
 !
 !*     31. Extended production/loss terms for chemical species
 !
-IF (LCHEMDIAG) THEN
+IF (LUSECHEM.OR.LCHEMDIAG) THEN
         CALL CH_INIT_BUDGET_n(ILUOUT)
         IF (NEQ_BUDGET>0) THEN
                 ALLOCATE(IINDEX(2,NNONZEROTERMS))
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index becf21fdf5e42b50ecb8526ab9ac2b8a368ad703..39fad63c3e76e1ef085c24f34a582113f3a6b451 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -223,6 +223,7 @@ END MODULE MODI_PHYS_PARAM_n
 !!                       2014  (M.Faivre)
 !!  06/2016     (G.Delautier) phasage surfex 8
 !!  2016 B.VIE LIMA
+!!      M. Leriche 02/2017 Avoid negative fluxes if sv=0 outside the physics domain
 !!-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -1292,7 +1293,12 @@ IF ( CTURB == 'TKEL' ) THEN
     ZSFRV(IIB-1,:)=ZSFRV(IIB,:)
     ZSFU(IIB-1,:)=ZSFU(IIB,:)
     ZSFV(IIB-1,:)=ZSFV(IIB,:)
-    IF (NSV>0)           ZSFSV(IIB-1,:,:)=ZSFSV(IIB,:,:)
+    IF (NSV>0)  THEN
+      ZSFSV(IIB-1,:,:)=ZSFSV(IIB,:,:)
+      WHERE ((ZSFSV(IIB-1,:,:).LT.0.).AND.(XSVT(IIB-1,:,IKB,:).EQ.0.))
+          ZSFSV(IIB-1,:,:) = 0.
+      END WHERE
+    ENDIF
     ZSFCO2(IIB-1,:)=ZSFCO2(IIB,:)
   END IF
   IF ( CLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
@@ -1300,7 +1306,12 @@ IF ( CTURB == 'TKEL' ) THEN
     ZSFRV(IIE+1,:)=ZSFRV(IIE,:)
     ZSFU(IIE+1,:)=ZSFU(IIE,:)
     ZSFV(IIE+1,:)=ZSFV(IIE,:)
-    IF (NSV>0)           ZSFSV(IIE+1,:,:)=ZSFSV(IIE,:,:)
+    IF (NSV>0) THEN
+      ZSFSV(IIE+1,:,:)=ZSFSV(IIE,:,:)
+      WHERE ((ZSFSV(IIE+1,:,:).LT.0.).AND.(XSVT(IIE+1,:,IKB,:).EQ.0.))
+          ZSFSV(IIE+1,:,:) = 0.
+      END WHERE
+    ENDIF
     ZSFCO2(IIE+1,:)=ZSFCO2(IIE,:)
   END IF
   IF ( CLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
@@ -1308,7 +1319,12 @@ IF ( CTURB == 'TKEL' ) THEN
     ZSFRV(:,IJB-1)=ZSFRV(:,IJB)
     ZSFU(:,IJB-1)=ZSFU(:,IJB)
     ZSFV(:,IJB-1)=ZSFV(:,IJB)
-    IF (NSV>0)           ZSFSV(:,IJB-1,:)=ZSFSV(:,IJB,:)
+    IF (NSV>0) THEN
+      ZSFSV(:,IJB-1,:)=ZSFSV(:,IJB,:)
+      WHERE ((ZSFSV(:,IJB-1,:).LT.0.).AND.(XSVT(:,IJB-1,IKB,:).EQ.0.))
+          ZSFSV(:,IJB-1,:) = 0.
+      END WHERE
+    ENDIF
     ZSFCO2(:,IJB-1)=ZSFCO2(:,IJB)
   END IF
   IF ( CLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
@@ -1316,7 +1332,12 @@ IF ( CTURB == 'TKEL' ) THEN
     ZSFRV(:,IJE+1)=ZSFRV(:,IJE)
     ZSFU(:,IJE+1)=ZSFU(:,IJE)
     ZSFV(:,IJE+1)=ZSFV(:,IJE)
-    IF (NSV>0)           ZSFSV(:,IJE+1,:)=ZSFSV(:,IJE,:)
+    IF (NSV>0) THEN
+      ZSFSV(:,IJE+1,:)=ZSFSV(:,IJE,:)
+      WHERE ((ZSFSV(:,IJE+1,:).LT.0.).AND.(XSVT(:,IJE+1,IKB,:).EQ.0.))
+          ZSFSV(:,IJE+1,:) = 0.
+      END WHERE
+    ENDIF
     ZSFCO2(:,IJE+1)=ZSFCO2(:,IJE)
   END IF
 !
@@ -1454,6 +1475,9 @@ XTIME_LES_BU_PROCESS = 0.
 IF (LUSECHEM) THEN
   CALL CH_MONITOR_n(ZWETDEPAER,KTCOUNT,XTSTEP, ILUOUT, NVERB)
 END IF
+DO JSV=1,NSV
+ XRSVS(:,:,:,JSV)   = MAX((XRSVS(:,:,:,JSV)),XSVMIN(JSV))
+END DO
 !
 ! For inert aerosol (dust and sea salt) => aer_monitor_n
 IF ((LDUST).OR.(LSALT)) THEN
diff --git a/src/MNH/prep_real_case.f90 b/src/MNH/prep_real_case.f90
index 659452539c6c43ec95affaa87ab732ce71cb1bd0..4df2194990a5e8e4f28d140829cd763f9f441368 100644
--- a/src/MNH/prep_real_case.f90
+++ b/src/MNH/prep_real_case.f90
@@ -370,6 +370,9 @@
 !!                  Jun   01, 2002 (O.Nuissier) filtering of tropical cyclone
 !!                  Aou   09, 2005 (D.Barbary) add CDADATMFILE CDADBOGFILE
 !!                   May   2006    Remove KEPS
+!!                  Feb   02, 2012 (C. Mari) interpolation from MOZART
+!!                                  add call to READ_CHEM_NETCDF_CASE &
+!!                                  VER_PREP_NETCDF_CASE 
 !!                  Mar   2012    Add NAM_NCOUT for netcdf output
 !!                  July  2013     (Bosseur & Filippi) Adds Forefire
 !!                  Mars  2014     (J.Escobar) Missing 'full' UPDATE_METRICS for arp2lfi // run
@@ -422,6 +425,8 @@ USE MODI_MNHREAD_ZS_DUMMY_n
 USE MODI_MNHWRITE_ZS_DUMMY_n
 USE MODI_COMPARE_DAD 
 USE MODI_PREP_SURF_MNH
+USE MODI_READ_CHEM_DATA_NETCDF_CASE
+USE MODI_VER_PREP_NETCDF_CASE
 !
 USE MODD_CONF            ! declaration modules
 USE MODD_CONF_n
@@ -739,6 +744,8 @@ IF(LEN_TRIM(YCHEMFILE)>0)THEN
   IF (GFOUND) READ(IPRE_REAL1,NAM_AERO_CONF)
   IF(YCHEMFILETYPE=='GRIBEX') &
   CALL READ_ALL_DATA_GRIB_CASE('CHEM',YPRE_REAL1,YCHEMFILE,YPGDFILE,ZHORI,NVERB,LDUMMY_REAL)
+  IF (YCHEMFILETYPE=='NETCDF') &
+  CALL READ_CHEM_DATA_NETCDF_CASE(YPRE_REAL1,YCHEMFILE,YPGDFILE,ZHORI,NVERB,LDUMMY_REAL)
 END IF
 !
 CALL CLOSE_ll(YPRE_REAL1, IOSTAT=IRESP)
@@ -876,6 +883,9 @@ END IF
 IF (LEN_TRIM(YCHEMFILE)>0 .AND. YCHEMFILETYPE=='GRIBEX') THEN
   CALL VER_PREP_GRIBEX_CASE('CHEM',ZDG)
 END IF
+IF (LEN_TRIM(YCHEMFILE)>0 .AND. YCHEMFILETYPE=='NETCDF') THEN
+  CALL VER_PREP_NETCDF_CASE(ZDG)
+END IF
 !
 CALL SECOND_MNH(ZTIME2)
 ZPREP = ZTIME2 - ZTIME1 - ZDG
diff --git a/src/MNH/profilern.f90 b/src/MNH/profilern.f90
index 64ed6ad1b478cb83c2d5643dccaf71e392798736..2c6bb61080c6fcf33910fab9928471b8a3e79364 100644
--- a/src/MNH/profilern.f90
+++ b/src/MNH/profilern.f90
@@ -407,13 +407,13 @@ ZTHV(:,:,:) = PTH(:,:,:) / (1.+WATER_SUM(PR(:,:,:,:)))*(1.+PR(:,:,:,1)/ZRDSRV)
 ZTEMPV(:,:,:)=ZTHV(:,:,:)*(PP(:,:,:)/ XP00) **(XRD/XCPD)
 CALL GPS_ZENITH_GRID(PR(:,:,:,1),ZTEMP,PP,ZZTD,ZZHD,ZZWD)
 ! Kunkel formulation
-IF (SIZE(PR,2) >= 2) THEN
+IF (SIZE(PR,4) >= 2) THEN
   WHERE ( PR(:,:,:,2) /=0 )
     ZVISIKUN(:,:,:) =0.027/(PR(:,:,:,2)*PRHODREF(:,:,:))**0.88
   END WHERE
 END IF
 ! Gultepe formulation
-IF ((SIZE(PR,2) >= 2) .AND. NSV_C2R2END /= 0 ) THEN 
+IF ((SIZE(PR,4) >= 2) .AND. NSV_C2R2END /= 0 ) THEN 
   WHERE ( (PR(:,:,:,2) /=0. ) .AND. (PSV(:,:,:,NSV_C2R2BEG+1) /=0. ) )
     ZVISI(:,:,:) =1.002/(PR(:,:,:,2)*PRHODREF(:,:,:)*PSV(:,:,:,NSV_C2R2BEG+1))**0.6473
   END WHERE
diff --git a/src/MNH/rain_c2r2_khko.f90 b/src/MNH/rain_c2r2_khko.f90
index 834f410cb629e2c1a832f7a3d95bc662bfa1fb39..77ed5a59085503f2edd108d9734e8d564847f868 100644
--- a/src/MNH/rain_c2r2_khko.f90
+++ b/src/MNH/rain_c2r2_khko.f90
@@ -218,6 +218,7 @@ END MODULE MODI_RAIN_C2R2_KHKO
 !!                            activation by cooling (OACTIT : mis en commentaires)
 !!      M.Mazoyer : 04/2016 : Add supersaturation diagnostics
 !!      C.Lac     : 07/2016 : Add droplet deposition
+!!      C.Lac     : 01/2017 : Correction on droplet deposition
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -1741,6 +1742,9 @@ INTEGER                           :: JL       ! and PACK intrinsics
 !  optimization by looking for locations where
 !  the precipitating fields are larger than a minimal value only !!!
 !
+IF (OSEDC) PINPRC (:,:) = 0.
+IF (LDEPOC) PINDEP (:,:) = 0.
+!
 DO JN = 1 , KSPLITR
   GSEDIM(:,:,:) = .FALSE.
   IF( OSEDC ) THEN
@@ -1933,10 +1937,12 @@ IF (LDEPOC) THEN
   GDEP(IIB:IIE,IJB:IJE) =    PRCS(IIB:IIE,IJB:IJE,2) >0 .AND. &
                                   PCCS(IIB:IIE,IJB:IJE,2) >0
   WHERE (GDEP)
-     PRCS(:,:,2) = PRCS(:,:,2) - XVDEPOC * PRCT(:,:,2) * PRHODJ(:,:,2)
-     PCCS(:,:,2) = PCCS(:,:,2) - XVDEPOC * PCCT(:,:,2) * PRHODJ(:,:,2)
-     PINPRC(:,:) = PINPRC(:,:) + XVDEPOC * PRCT(:,:,2) * PRHODJ(:,:,2) /XRHOLW 
-     PINDEP(:,:) = XVDEPOC * PRCT(:,:,2) * PRHODJ(:,:,2) /XRHOLW 
+     PRCS(:,:,2) = PRCS(:,:,2) - XVDEPOC * PRCS(:,:,2)*PTSTEP / &
+                  ( PZZ(:,:,3) - PZZ(:,:,2))
+     PCCS(:,:,2) = PCCS(:,:,2) - XVDEPOC * PCCS(:,:,2)*PTSTEP / &          
+                  ( PZZ(:,:,3) - PZZ(:,:,2))
+     PINPRC(:,:) = PINPRC(:,:) + XVDEPOC * PRCS(:,:,2) * PTSTEP * PRHODREF(:,:,2) /XRHOLW
+     PINDEP(:,:) = XVDEPOC * PRCS(:,:,2) * PTSTEP * PRHODREF(:,:,2) /XRHOLW
   END WHERE
 END IF
 !
diff --git a/src/MNH/rain_ice.f90 b/src/MNH/rain_ice.f90
index e2f9fd9196b28fe58f31914d89fe0ffd29b6926d..2647d09e05c73b2f53d5275a25e02b6187ae16a4 100644
--- a/src/MNH/rain_ice.f90
+++ b/src/MNH/rain_ice.f90
@@ -238,7 +238,8 @@ END MODULE MODI_RAIN_ICE
 !!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
 !!                                      aircraft, ballon and profiler
 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
-!!      C.LAc : 10/2016 : add droplets depposition
+!!      C.Lac : 10/2016 : add droplet deposition
+!!      C.Lac : 01/2017 : correction on droplet deposition
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -1738,9 +1739,9 @@ IF (LDEPOSC) THEN
   GDEP(:,:) = .FALSE.
   GDEP(IIB:IIE,IJB:IJE) =    PRCS(IIB:IIE,IJB:IJE,2) >0 
   WHERE (GDEP)
-     PRCS(:,:,2) = PRCS(:,:,2) - XVDEPOSC * PRCT(:,:,2) * PRHODJ(:,:,2)
-     PINPRC(:,:) = PINPRC(:,:) + XVDEPOSC * PRCT(:,:,2) * PRHODJ(:,:,2) /XRHOLW 
-     PINDEP(:,:) = XVDEPOSC * PRCT(:,:,2) * PRHODJ(:,:,2) /XRHOLW 
+     PRCS(:,:,2) = PRCS(:,:,2) - XVDEPOSC * PRCS(:,:,2)*PTSTEP / PDZZ(:,:,2)
+     PINPRC(:,:) = PINPRC(:,:) + XVDEPOSC * PRCS(:,:,2) * PTSTEP * PRHODREF(:,:,2) /XRHOLW 
+     PINDEP(:,:) = XVDEPOSC * PRCS(:,:,2) * PTSTEP * PRHODREF(:,:,2) /XRHOLW 
   END WHERE
 END IF
 !
diff --git a/src/MNH/read_chem_data_netcdf_case.f90 b/src/MNH/read_chem_data_netcdf_case.f90
new file mode 100644
index 0000000000000000000000000000000000000000..247462703e1c50e022daf10ccd3196c26c5a4eb5
--- /dev/null
+++ b/src/MNH/read_chem_data_netcdf_case.f90
@@ -0,0 +1,821 @@
+!     ################################
+      MODULE MODI_READ_CHEM_DATA_NETCDF_CASE
+!     #################################
+INTERFACE
+SUBROUTINE READ_CHEM_DATA_NETCDF_CASE(HPRE_REAL1,HFILE,HPGDFILE,      &
+                    PTIME_HORI,KVERB,ODUMMY_REAL                         ) 
+!
+CHARACTER(LEN=28),  INTENT(IN) :: HPRE_REAL1 ! name of the PRE_REAL1 file
+CHARACTER(LEN=28),  INTENT(IN) :: HFILE    ! name of the NETCDF file
+CHARACTER(LEN=28),  INTENT(IN) :: HPGDFILE ! name of the physiographic data file
+INTEGER,           INTENT(IN)  :: KVERB    ! verbosity level
+LOGICAL,           INTENT(IN)  :: ODUMMY_REAL! flag to interpolate dummy fields
+REAL,           INTENT(INOUT)  :: PTIME_HORI ! time spent in hor. interpolations
+END SUBROUTINE READ_CHEM_DATA_NETCDF_CASE
+!
+END INTERFACE
+END MODULE MODI_READ_CHEM_DATA_NETCDF_CASE
+!     ##########################################################################
+      SUBROUTINE READ_CHEM_DATA_NETCDF_CASE(HPRE_REAL1,HFILE,HPGDFILE,      &
+                       PTIME_HORI,KVERB,ODUMMY_REAL                            )
+!     ##########################################################################
+!
+!!****  *READ_CHEM_DATA_NETCDF_CASE* - reads data for the initialization of real cases.
+!!
+!!    PURPOSE
+!!    -------
+!     This routine reads the two input files :
+!       The PGD which is closed after reading
+!       The NETCDF file
+!     Projection is read in READ_LFIFM_PGD (MODD_GRID).
+!     Grid and definition of large domain are read in PGD file and 
+!           NETCDF files.
+!     The PGD files are also read in READ_LFIFM_PGD.
+!     The PGD file is closed.
+!     Vertical grid is defined in READ_VER_GRID.
+!     PGD fields are stored on MESO-NH domain (in TRUNC_PGD).
+!!
+!!**  METHOD
+!!    ------
+!!  0. Declarations
+!!    1. Declaration of arguments
+!!    2. Declaration of local variables
+!!  1. Read PGD file
+!!    1. Domain restriction
+!!    2. Coordinate conversion to lat,lon system
+!!  2. Read Netcdf fields
+!!  3. Vertical grid
+!!  4. Free all temporary allocations
+!!
+!!    EXTERNAL
+!!    --------
+!!    subroutine READ_LFIFM_PGD    : to read PGD file
+!!    subroutine READ_VER_GRID     : to read the vertical grid in namelist file.
+!!    subroutine HORIBL            : horizontal bilinear interpolation
+!!    subroutine XYTOLATLON        : projection from conformal to lat,lon
+!!
+!!    function   FMLOOK_ll         : to retrieve the logical unit associated with a file
+!!
+!!    Module     MODI_READ_VER_GRID     : interface for subroutine READ_VER_GRID
+!!    Module     MODI_HORIBL            : interface for subroutine HORIBL
+!!    Module     MODI_XYTOLATLON        : interface for subroutine XYTOLATLON
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!
+!!      Module MODD_CONF      : contains configuration variables for all models.
+!!         NVERB : verbosity level for output-listing
+!!      Module MODD_LUNIT     : contains logical unit names for all models
+!!         CLUOUT0 : name of output-listing
+!!      Module MODD_PGDDIM    : contains dimension of PGD fields
+!!         NPGDIMAX: dimension along x (no external point)
+!!         NPGDJMAX: dimension along y (no external point)
+!!      Module MODD_PARAMETERS
+!!         JPHEXT
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    23/01/12 (C. Mari) 
+!!      A. Berger   20/03/12 adapt whatever the chemical mechanism in BASIC
+!-------------------------------------------------------------------------------
+!
+!*      0. DECLARATIONS
+!------------
+USE MODE_FM
+USE MODE_IO_ll
+USE MODE_TIME
+!
+USE MODI_READ_HGRID_n
+USE MODI_READ_VER_GRID
+USE MODI_XYTOLATLON
+USE MODI_HORIBL
+USE MODI_INI_NSV
+USE MODI_CH_INIT_SCHEME_n
+USE MODI_CH_AER_INIT_SOA
+!
+USE MODD_CONF
+USE MODD_CONF_n
+USE MODD_CST
+USE MODD_LUNIT
+USE MODD_PARAMETERS
+USE MODD_GRID
+USE MODD_GRID_n
+USE MODD_DIM_n
+USE MODD_PARAM_n, ONLY : CTURB
+USE MODD_TIME
+USE MODD_TIME_n
+USE MODD_CH_MNHC_n, ONLY : LUSECHEM,LUSECHAQ,LUSECHIC,LCH_PH
+USE MODD_CH_M9_n,   ONLY : NEQ ,  CNAMES
+USE MODD_CH_AEROSOL, ONLY: CORGANIC, NCARB, NSOA, NSP, LORILAM,&
+                           JPMODE, LVARSIGI, LVARSIGJ,CAERONAMES
+USE MODD_NSV  
+USE MODD_PREP_REAL
+USE MODE_MODELN_HANDLER
+!JUAN REALZ
+USE MODE_MPPDB
+!JUAN REALZ
+USE MODI_CH_OPEN_INPUT  
+USE MODE_THERMO
+!
+USE MODD_BLANK
+!
+IMPLICIT NONE
+!
+include 'netcdf.inc'
+!
+!* 0.1. Declaration of arguments
+!       ------------------------
+!
+CHARACTER(LEN=28),      INTENT(IN) :: HPRE_REAL1 ! name of the PRE_REAL1 file
+CHARACTER(LEN=28),      INTENT(IN) :: HFILE      ! name of the GRIB file
+CHARACTER(LEN=28),      INTENT(IN) :: HPGDFILE   ! name of the physiographic data file
+INTEGER,               INTENT(IN)  :: KVERB    ! verbosity level
+LOGICAL,               INTENT(IN)  :: ODUMMY_REAL! flag to interpolate dummy fields
+REAL,                INTENT(INOUT) :: PTIME_HORI ! time spent in hor. interpolations
+!
+!* 0.2 Declaration of local variables
+!      ------------------------------
+! General purpose variables
+INTEGER                            :: ILUOUT0       ! Unit used for output msg.
+INTEGER                            :: IRESP   ! Return code of FM-routines
+INTEGER                            :: IRET          ! Return code from subroutines
+INTEGER                            :: JI,JJ,JK      ! Dummy counters
+INTEGER                            :: JLOOP1        !  |
+INTEGER                            :: JNCHEM, JNAER ! conters of chemical species in BASIC
+! Variables used by the PGD reader
+CHARACTER(LEN=28)                  :: YPGD_NAME     ! not used - dummy argument
+CHARACTER(LEN=28)                  :: YPGD_DAD_NAME ! not used - dummy argument
+CHARACTER(LEN=2)                   :: YPGD_TYPE     ! not used - dummy argument
+! PGD Grib definition variables
+INTEGER                            :: INO           ! Number of points of the grid
+INTEGER                            :: IIU           ! Number of points along X
+INTEGER                            :: IJU           ! Number of points along Y
+REAL, DIMENSION(:), ALLOCATABLE    :: ZLONOUT       ! mapping PGD -> Grib (lon.)
+REAL, DIMENSION(:), ALLOCATABLE    :: ZLATOUT       ! mapping PGD -> Grib (lat.)
+REAL, DIMENSION(:,:), ALLOCATABLE  :: ZXM           ! X of PGD mass points
+REAL, DIMENSION(:,:), ALLOCATABLE  :: ZYM           ! Y of PGD mass points
+REAL, DIMENSION(:,:), ALLOCATABLE  :: ZLATM         ! Lat of PGD mass points
+REAL, DIMENSION(:,:), ALLOCATABLE  :: ZLONM         ! Lon of PGD mass points
+! Variable involved in the task of reading the netcdf  file
+REAL,DIMENSION(:,:),ALLOCATABLE    :: ZVALUE        ! Intermediate array
+REAL,DIMENSION(:),ALLOCATABLE      :: ZVALUE1D        ! Intermediate array
+REAL,DIMENSION(:,:),ALLOCATABLE    :: ZOUT          ! Intermediate arrays
+REAL,DIMENSION(:),ALLOCATABLE      :: ZOUT1D          ! Intermediate arrays
+INTEGER                            :: ind_netcdf    ! Indice for netcdf var.
+!chemistry field infile MOZ1.nam
+INTEGER                                       :: ICHANNEL
+CHARACTER(LEN=8)                              :: YMOZ="MOZ1.nam"
+integer                                       :: IMOZ
+CHARACTER(LEN=68)                             :: YFORMAT
+CHARACTER(LEN=40), DIMENSION(:), ALLOCATABLE  :: YSPCMNH
+integer, dimension(:), ALLOCATABLE            :: ISPCMOZ
+CHARACTER(LEN=9)                              :: YA
+REAL,DIMENSION(:,:),ALLOCATABLE               :: ZCOEFMOZART
+CHARACTER(LEN=18),dimension(:,:),ALLOCATABLE  :: YCHANGE
+type TZMOZ
+real                                          :: ZCOEFMOZ
+character(16)                                 :: YSPCMOZ
+end type TZMOZ
+type(TZMOZ), DIMENSION(:,:),ALLOCATABLE       :: TZSTOC
+! model indice
+INTEGER                           :: IMI
+!
+! For netcdf 
+!
+integer                               :: status, ncid, varid
+integer :: lat_varid, lon_varid, lev_varid, time_varid 
+integer :: hyam_varid, hybm_varid, p0_varid, t_varid, q_varid, ps_varid 
+integer :: recid, latid, lonid, levid, timeid
+integer :: latlen, lonlen, levlen, nrecs,timelen
+integer :: itimeindex, KILEN, jrec
+CHARACTER(LEN=40)                     :: recname
+REAL, DIMENSION(:), ALLOCATABLE       :: lats
+REAL, DIMENSION(:), ALLOCATABLE       :: lons 
+REAL, DIMENSION(:), ALLOCATABLE       :: levs 
+INTEGER, DIMENSION(:), ALLOCATABLE    :: count3d, start3d
+INTEGER, DIMENSION(:), ALLOCATABLE    :: count2d, start2d 
+REAL, DIMENSION(:), ALLOCATABLE       :: time, hyam, hybm 
+REAL                                  :: p0 
+INTEGER, DIMENSION(:), ALLOCATABLE    :: kinlo 
+REAL, DIMENSION(:,:,:), ALLOCATABLE   :: vartemp3d,vartemp3dbis,vartemp3dter 
+REAL, DIMENSION(:,:,:), ALLOCATABLE   :: vartemp3dquater 
+REAL, DIMENSION(:,:,:), ALLOCATABLE   :: ZCHEMMOZ, TMOZ, QMOZ
+REAL, DIMENSION(:,:), ALLOCATABLE     :: PSMOZ 
+
+real ::a,b
+
+!----------------------------------------------------------------------
+!
+IMI = GET_CURRENT_MODEL_INDEX()
+!
+!* 1. READ PGD FILE
+!     -------------
+!
+CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRET)
+CALL READ_HGRID_n(HPGDFILE,YPGD_NAME,YPGD_DAD_NAME,YPGD_TYPE)
+!
+! 1.1 Domain restriction
+!
+CALL GET_DIM_EXT_ll('B',IIU,IJU)
+INO = IIU * IJU
+!
+!
+! 1.2 Coordinate conversion to lat,lon system
+!
+ALLOCATE (ZXM(IIU,IJU))
+ALLOCATE (ZYM(IIU,IJU))
+ALLOCATE (ZLONM(IIU,IJU))
+ALLOCATE (ZLATM(IIU,IJU))
+ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
+ZXM(IIU,1)     = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
+ZXM(:,2:IJU)   = SPREAD(ZXM(:,1),2,IJU-1)
+ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
+ZYM(1,IJU)     = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
+ZYM(2:IIU,:)   = SPREAD(ZYM(1,:),1,IIU-1)
+CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
+                      IIU,IJU)
+ALLOCATE (ZLONOUT(INO))
+ALLOCATE (ZLATOUT(INO))
+JLOOP1 = 0
+DO JJ = 1, IJU
+  ZLONOUT(JLOOP1+1:JLOOP1+IIU) = ZLONM(1:IIU,JJ)
+  ZLATOUT(JLOOP1+1:JLOOP1+IIU) = ZLATM(1:IIU,JJ)
+  JLOOP1 = JLOOP1 + IIU
+ENDDO
+DEALLOCATE (ZYM)
+DEALLOCATE (ZXM)
+!
+!
+!* 2. READ NETCDF FIELDS
+!     ------------------
+!
+! 2.1 Open netcdf files
+print*,'Open netcdf files:',HFILE
+!
+status = nf_open(HFILE, nf_nowrite, ncid) 
+if (status /= nf_noerr) call handle_err(status)
+!
+! 2.2 Read netcdf files
+!
+! get dimension IDs
+!
+!* get dimension ID of unlimited variable in netcdf file
+status = nf_inq_unlimdim(ncid, recid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimid(ncid, "lat", latid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimid(ncid, "lon", lonid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimid(ncid, "lev", levid)
+if (status /= nf_noerr) call handle_err(status)
+!
+! get dimensions
+!
+!* get dimension and name of unlimited variable in netcdf file
+status = nf_inq_dim(ncid, recid, recname, nrecs)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimlen(ncid, latid, latlen)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimlen(ncid, lonid, lonlen)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_dimlen(ncid, levid, levlen)
+if (status /= nf_noerr) call handle_err(status)
+print*, latlen, lonlen, levlen, nrecs
+!
+! get variable IDs
+!
+status = nf_inq_varid(ncid, "lat", lat_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "lon", lon_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "lev", lev_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "time", time_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "P0", p0_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "hyam", hyam_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "hybm", hybm_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "T", t_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "Q", q_varid)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_inq_varid(ncid, "PS", ps_varid)
+if (status /= nf_noerr) call handle_err(status)
+!
+KILEN = latlen * lonlen
+!
+! 2.3 Read data.
+!
+ALLOCATE (count3d(4))
+ALLOCATE (start3d(4))
+ALLOCATE (count2d(3))
+ALLOCATE (start2d(3))
+ALLOCATE (lats(latlen))
+ALLOCATE (lons(lonlen))
+ALLOCATE (levs(levlen))
+ALLOCATE (time(nrecs))
+ALLOCATE (kinlo(latlen))
+kinlo(:) = lonlen
+ALLOCATE (vartemp3d(lonlen,latlen,levlen))
+ALLOCATE (vartemp3dbis(lonlen,latlen,levlen))
+ALLOCATE (vartemp3dter(lonlen,latlen,levlen))
+ALLOCATE (vartemp3dquater(lonlen,latlen,levlen))
+ALLOCATE (ZCHEMMOZ(lonlen,latlen,levlen))
+ALLOCATE (TMOZ(lonlen,latlen,levlen))
+ALLOCATE (QMOZ(lonlen,latlen,levlen))
+ALLOCATE (PSMOZ(lonlen,latlen))
+ALLOCATE (XA_SV_LS(levlen))
+ALLOCATE (hyam(levlen))
+ALLOCATE (XB_SV_LS(levlen))
+ALLOCATE (hybm(levlen))
+ALLOCATE (XT_SV_LS(IIU,IJU,levlen))
+ALLOCATE (XQ_SV_LS(IIU,IJU,levlen,1))
+ALLOCATE (XPS_SV_LS(IIU,IJU))
+ALLOCATE (XZS_SV_LS(IIU,IJU))
+! take the orography from ECMWF
+XZS_SV_LS(:,:) = XZS_LS(:,:)
+!
+! get values of variables
+!
+status = nf_get_var_double(ncid, lat_varid, lats(:))
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, lon_varid, lons(:))
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, lev_varid, levs(:))
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, time_varid, time(:))
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, hyam_varid, hyam)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, hybm_varid, hybm)
+if (status /= nf_noerr) call handle_err(status)
+status = nf_get_var_double(ncid, p0_varid, p0)
+if (status /= nf_noerr) call handle_err(status)
+XP00_SV_LS = p0
+!
+! hyam and hybm coefficients for pressure calculations have to be reversed 
+! from top-bottom to bottom-up direction
+do JJ = 1, levlen
+  XA_SV_LS(JJ) = hyam(levlen+1-JJ)
+  XB_SV_LS(JJ) = hybm(levlen+1-JJ)
+end do
+!
+!
+!     Read 1 record of lon*lat*lev values, starting at the
+!     beginning of the record (the (1, 1, 1, rec) element in the netCDF
+!     file).
+ count3d(1) = lonlen
+ count3d(2) = latlen
+ count3d(3) = levlen
+ count3d(4) = 1
+ start3d(1) = 1
+ start3d(2) = 1
+ start3d(3) = 1
+! Choose time index according to the chosen time in namelist
+! 1 for 00 - 2 for 06 - 3 for 12 - 4 for 18
+IF (CDUMMY1=="00") THEN 
+        itimeindex=1
+ELSEIF (CDUMMY1=="06") THEN
+        itimeindex=2
+ELSEIF (CDUMMY1=="12") THEN
+        itimeindex=3
+ELSEIF (CDUMMY1=="18") THEN
+        itimeindex=4
+ENDIF
+ start3d(4) = itimeindex
+!
+  status = nf_get_vara_double(ncid, t_varid, start3d, count3d, &
+           vartemp3d)
+  if (status /= nf_noerr) call handle_err(status)
+!
+do JJ=1,levlen
+! lev, lat, lon
+ TMOZ(:,:,JJ) = vartemp3d(:,:,levlen+1-JJ)
+enddo
+!
+  status = nf_get_vara_double(ncid, q_varid, start3d, count3d, &
+           vartemp3d)
+  if (status /= nf_noerr) call handle_err(status)
+!
+do JJ=1,levlen
+! lev, lat, lon
+ QMOZ(:,:,JJ) = vartemp3d(:,:,levlen+1-JJ)
+enddo
+!
+ count2d(1) = lonlen
+ count2d(2) = latlen
+ count2d(3) = 1
+ start2d(1) = 1
+ start2d(2) = 1
+ start2d(3) = itimeindex
+  status = nf_get_vara_double(ncid, ps_varid, start2d, count2d, PSMOZ(:,:))
+  if (status /= nf_noerr) call handle_err(status)
+
+  
+!------------------------------------------------------------------------
+!* 3 Interpolation of MOZART variable
+!---------------------------------------------------------------------
+  ! Always initialize chemical scheme variables before INI_NSV call !
+  CALL CH_INIT_SCHEME_n(IMI,LUSECHAQ,LUSECHIC,LCH_PH,ILUOUT0,KVERB)
+  LUSECHEM = .TRUE.
+  IF (LORILAM) THEN
+    CORGANIC = "MPMPO"
+    LVARSIGI = .TRUE.
+    LVARSIGJ = .TRUE.
+    CALL CH_AER_INIT_SOA(ILUOUT0, KVERB)
+  END IF
+  ! initialise NSV_* variables
+  CALL INI_NSV(1)
+    DEALLOCATE(XSV_LS)
+    ALLOCATE (XSV_LS(IIU,IJU,levlen,NSV))
+   XSV_LS(:,:,:,:) = 0.
+!
+  WRITE (ILUOUT0,'(A,A4,A)') ' | Reading MOZART species (ppp) from ',HFILE,' file'
+
+where (ZLONOUT(:) < 0.) ZLONOUT(:) = ZLONOUT(:) + 360.
+!
+ALLOCATE(ZVALUE(levlen,KILEN))
+ALLOCATE(ZOUT(levlen,INO))
+ALLOCATE(ZVALUE1D(KILEN))
+ALLOCATE(ZOUT1D(INO))
+
+!
+!*       2.6.1  read MOZART species from the file MOZ1.nam
+!
+! open input file
+CALL CH_OPEN_INPUT(YMOZ,"MOZ2MESONH",ICHANNEL,ILUOUT0,KVERB)
+!
+!read number of mocage species to transfer into mesonh
+READ(ICHANNEL, *) IMOZ
+IF (KVERB >= 5) WRITE (ILUOUT0,*) "number of mozart species to transfer into &
+& mesonh : ", IMOZ
+!
+!read data input format
+READ(ICHANNEL,"(A)") YFORMAT
+YFORMAT=UPCASE(YFORMAT)
+IF (KVERB >= 5) WRITE (ILUOUT0,*) "input format is: ", YFORMAT
+!
+!allocate fields
+ALLOCATE(YSPCMNH(IMOZ))      !MESONH species
+ALLOCATE(TZSTOC(IMOZ,4))     !MOZART coefficient and MOZART species associated
+ALLOCATE(ISPCMOZ(IMOZ))      !MOZART species number into MESONH species
+ALLOCATE(ZCOEFMOZART(IMOZ,4))!Coef stoich of each MOZART species
+ALLOCATE(YCHANGE(IMOZ,4))    !MOZART species with _VMR_inst
+!read MESONH variable names and MOZART variable names associated 
+DO JI = 1,IMOZ               !for every MNH species existing in MOZ1.nam                              
+  READ(ICHANNEL,YFORMAT) YSPCMNH(JI), ISPCMOZ(JI), TZSTOC(JI,1)%ZCOEFMOZ,& !reading line by line
+                 TZSTOC(JI,1)%YSPCMOZ, TZSTOC(JI,2)%ZCOEFMOZ,&             !of string
+                 TZSTOC(JI,2)%YSPCMOZ, TZSTOC(JI,3)%ZCOEFMOZ,&
+                 TZSTOC(JI,3)%YSPCMOZ, TZSTOC(JI,4)%ZCOEFMOZ,&
+                 TZSTOC(JI,4)%YSPCMOZ
+  WRITE(ILUOUT0,YFORMAT) YSPCMNH(JI), ISPCMOZ(JI),&                        !writing in arrays
+                TZSTOC(JI,1)%ZCOEFMOZ, TZSTOC(JI,1)%YSPCMOZ,&
+                TZSTOC(JI,2)%ZCOEFMOZ, TZSTOC(JI,2)%YSPCMOZ,&
+                TZSTOC(JI,3)%ZCOEFMOZ, TZSTOC(JI,3)%YSPCMOZ,&
+                TZSTOC(JI,4)%ZCOEFMOZ, TZSTOC(JI,4)%YSPCMOZ
+!
+  ZCOEFMOZART(JI,1) =  (TZSTOC(JI,1)%ZCOEFMOZ) !coef stoich of each MOZART species set into an array 
+  ZCOEFMOZART(JI,2) =  (TZSTOC(JI,2)%ZCOEFMOZ) 
+  ZCOEFMOZART(JI,3) =  (TZSTOC(JI,3)%ZCOEFMOZ)
+  ZCOEFMOZART(JI,4) =  (TZSTOC(JI,4)%ZCOEFMOZ)
+! 
+  YA="_VMR_inst"
+  YCHANGE(JI,1)=trim(TZSTOC(JI,1)%YSPCMOZ)//YA !set into an array MOZART species with _VMR_inst
+  YCHANGE(JI,2)=trim(TZSTOC(JI,2)%YSPCMOZ)//YA 
+  YCHANGE(JI,3)=trim(TZSTOC(JI,3)%YSPCMOZ)//YA
+  YCHANGE(JI,4)=trim(TZSTOC(JI,4)%YSPCMOZ)//YA
+!
+!* exchange mozart values onto prognostic variables XSV_LS
+! and convert MOZART fields to 2D for use in horizontal interpolation 
+! routine HORIBL.f90
+!
+  DO JNCHEM = NSV_CHEMBEG, NSV_CHEMEND  !loop on all MNH species
+    IF (trim(CNAMES(JNCHEM-NSV_CHEMBEG+1))==trim(YSPCMNH(JI))) THEN !MNH mechanism species
+       IF (ISPCMOZ(JI)==1) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==2) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ) + &
+                               ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ) 
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==3) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,3)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dter)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)+&
+                            ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ)+&
+                            ZCOEFMOZART(JI,3)*vartemp3dter(:,:,levlen+1-JJ)
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==4) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,3)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dter)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,4)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                         vartemp3dquater)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,3)*vartemp3dter(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,4)*vartemp3dquater(:,:,levlen+1-JJ)
+         ENDDO
+       ENDIF
+       DO JK = 1, levlen
+         JLOOP1 = 0
+         DO JJ = 1, latlen
+           ZVALUE(JK,JLOOP1+1:JLOOP1+lonlen) = ZCHEMMOZ(1:lonlen,JJ,JK)
+           JLOOP1 = JLOOP1+lonlen
+         ENDDO                                                                                           
+         CALL HORIBL(lats(1),lons(1),lats(latlen),lons(lonlen), &
+                     latlen,kinlo,KILEN,                        &
+                     ZVALUE(JK,:),INO,ZLONOUT,ZLATOUT,          &
+                     ZOUT(JK,:),.FALSE.,PTIME_HORI,.TRUE.)
+         CALL ARRAY_1D_TO_2D(INO,ZOUT(JK,:),IIU,IJU, &
+                             XSV_LS(:,:,JK,JNCHEM)   )
+       ENDDO  ! levlen
+    ENDIF      
+
+  ENDDO ! JNCHEM
+  DO JNAER = NSV_AERBEG, NSV_AEREND 
+    IF (trim(CAERONAMES(JNAER-NSV_AERBEG+1))==trim(YSPCMNH(JI))) THEN !MNH mechanism species
+       IF (ISPCMOZ(JI)==1) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==2) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ) + &
+                               ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ) 
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==3) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,3)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dter)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)+&
+                            ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ)+&
+                            ZCOEFMOZART(JI,3)*vartemp3dter(:,:,levlen+1-JJ)
+         ENDDO
+       ELSE IF (ISPCMOZ(JI)==4) THEN
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,1)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                               vartemp3d)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,2)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dbis)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,3)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                            vartemp3dter)
+         if (status /= nf_noerr) call  handle_err(status)
+         status = nf_inq_varid(ncid, trim(YCHANGE(JI,4)), ind_netcdf)
+         if (status /= nf_noerr) call handle_err(status)
+         status = nf_get_vara_double(ncid,ind_netcdf, start3d, count3d, &
+                                                         vartemp3dquater)
+         if (status /= nf_noerr) call  handle_err(status)
+         DO JJ=1,levlen ! lev, lat, lon
+           ZCHEMMOZ(:,:,JJ)=ZCOEFMOZART(JI,1)*vartemp3d(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,2)*vartemp3dbis(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,3)*vartemp3dter(:,:,levlen+1-JJ)+&
+                               ZCOEFMOZART(JI,4)*vartemp3dquater(:,:,levlen+1-JJ)
+         ENDDO
+       ENDIF
+       DO JK = 1, levlen
+         JLOOP1 = 0
+         DO JJ = 1, latlen
+           ZVALUE(JK,JLOOP1+1:JLOOP1+lonlen) = ZCHEMMOZ(1:lonlen,JJ,JK)
+           JLOOP1 = JLOOP1+lonlen
+         ENDDO                                                                                           
+         CALL HORIBL(lats(1),lons(1),lats(latlen),lons(lonlen), &
+                     latlen,kinlo,KILEN,                        &
+                     ZVALUE(JK,:),INO,ZLONOUT,ZLATOUT,          &
+                     ZOUT(JK,:),.FALSE.,PTIME_HORI,.TRUE.)
+         CALL ARRAY_1D_TO_2D(INO,ZOUT(JK,:),IIU,IJU, &
+                             XSV_LS(:,:,JK,JNAER)   )
+       ENDDO  ! levlen
+    ENDIF         
+  ENDDO ! JNAER
+ENDDO  ! JIDO JNCHEM = NSV_CHEMBEG, NSV_CHEMEND  !loop on all MNH species
+DEALLOCATE(YSPCMNH) 
+DEALLOCATE(TZSTOC)
+DEALLOCATE(ISPCMOZ) 
+DEALLOCATE(ZCOEFMOZART)
+DEALLOCATE(YCHANGE)
+!
+XSV_LS(:,:,:,:) = MAX(XSV_LS(:,:,:,:),0.)
+!
+DO JK = 1, levlen
+  JLOOP1 = 0
+  DO JJ = 1, latlen
+    ZVALUE(JK,JLOOP1+1:JLOOP1+lonlen) = TMOZ(1:lonlen,JJ,JK)
+    JLOOP1 = JLOOP1 + lonlen
+  ENDDO
+  CALL HORIBL(lats(1),lons(1),lats(latlen),lons(lonlen), &
+              latlen,kinlo,KILEN,                        &
+              ZVALUE(JK,:),INO,ZLONOUT,ZLATOUT,          &
+              ZOUT(JK,:),.FALSE.,PTIME_HORI,.FALSE.)
+!
+  CALL ARRAY_1D_TO_2D(INO,ZOUT(JK,:),IIU,IJU, &
+                      XT_SV_LS(:,:,JK))
+ENDDO 
+XT_SV_LS(:,:,:) = MAX(XT_SV_LS(:,:,:),0.)
+!
+DO JK = 1, levlen
+  JLOOP1 = 0
+  DO JJ = 1, latlen
+    ZVALUE(JK,JLOOP1+1:JLOOP1+lonlen) = QMOZ(1:lonlen,JJ,JK)
+    JLOOP1 = JLOOP1 + lonlen
+  ENDDO
+  CALL HORIBL(lats(1),lons(1),lats(latlen),lons(lonlen), &
+              latlen,kinlo,KILEN,                                &
+              ZVALUE(JK,:),INO,ZLONOUT,ZLATOUT,                  &
+              ZOUT(JK,:),.FALSE.,PTIME_HORI,.FALSE.)
+!
+   CALL ARRAY_1D_TO_2D(INO,ZOUT(JK,:),IIU,IJU,                    &
+                       XQ_SV_LS(:,:,JK,1))
+ENDDO 
+XQ_SV_LS(:,:,:,1) = MAX(XQ_SV_LS(:,:,:,1),0.)
+!
+JLOOP1 = 0
+DO JJ = 1, latlen
+  ZVALUE1D(JLOOP1+1:JLOOP1+lonlen) = PSMOZ(1:lonlen,JJ)
+  JLOOP1 = JLOOP1 + lonlen
+ENDDO
+CALL HORIBL(lats(1),lons(1),lats(latlen),lons(lonlen), &
+            latlen,kinlo,KILEN,                                &
+            ZVALUE1D(:),INO,ZLONOUT,ZLATOUT,                  &
+            ZOUT1D(:),.FALSE.,PTIME_HORI,.FALSE.)
+!
+CALL ARRAY_1D_TO_2D(INO,ZOUT1D(:),IIU,IJU,                    &
+                    XPS_SV_LS(:,:))
+XPS_SV_LS(:,:) = MAX(XPS_SV_LS(:,:),0.)
+!
+!
+!
+! close the netcdf file
+!nf_close
+status = nf_close(ncid) 
+if (status /= nf_noerr) call handle_err(status)
+!
+  DEALLOCATE (ZVALUE)
+  DEALLOCATE (ZOUT)
+  DEALLOCATE (ZVALUE1D) 
+  DEALLOCATE (ZOUT1D)
+!!
+
+! close
+! file
+CALL CLOSE_ll(YMOZ,IOSTAT=IRET)
+
+
+!-------------------------------------------------------------
+!
+!* 4. VERTICAL GRID
+!
+!* 4.1 Read VERTICAL GRID
+!
+WRITE (ILUOUT0,'(A)') ' | Reading of vertical grid in progress'
+CALL READ_VER_GRID(HPRE_REAL1)
+!
+!--------------------------------------------------------------
+!
+!* 4.2 Interpolate on Meso-NH VERTICAL GRID
+!
+!* 4.3 Free all temporary allocations
+!
+DEALLOCATE (ZLATOUT)
+DEALLOCATE (ZLONOUT)
+DEALLOCATE (hyam)
+DEALLOCATE (hybm)
+DEALLOCATE (vartemp3d)
+DEALLOCATE (vartemp3dbis)
+DEALLOCATE (vartemp3dter)
+DEALLOCATE (vartemp3dquater)
+!
+WRITE (ILUOUT0,'(A,A4,A)') ' -- netcdf decoder for ',HFILE,' file ended successfully'
+!
+!
+CONTAINS
+!
+!     #############################
+      SUBROUTINE HANDLE_ERR(STATUS)
+!     #############################
+     INTEGER STATUS
+     IF (STATUS .NE. NF_NOERR) THEN
+        PRINT *, NF_STRERROR(STATUS)
+     STOP 'Stopped'
+     ENDIF
+     END SUBROUTINE HANDLE_ERR
+!
+!
+!     #############################################
+      SUBROUTINE ARRAY_1D_TO_2D (KN1,P1,KL1,KL2,P2)
+!     #############################################
+!
+!       Small routine used to store a linear array into a 2 dimension array
+!
+IMPLICIT NONE
+INTEGER,                INTENT(IN)  :: KN1
+REAL,DIMENSION(KN1),    INTENT(IN)  :: P1
+INTEGER,                INTENT(IN)  :: KL1
+INTEGER,                INTENT(IN)  :: KL2
+REAL,DIMENSION(KL1,KL2),INTENT(OUT) :: P2
+INTEGER                 :: JLOOP1_A1T2
+INTEGER                 :: JLOOP2_A1T2
+INTEGER                 :: JPOS_A1T2
+!
+IF (KN1 < KL1*KL2) THEN
+  WRITE (ILUOUT0,'(A)') ' | Error in "ARRAY_1D_TO_2D", sizes do not match - abort'
+ !callabortstop
+  CALL CLOSE_ll(CLUOUT0,IOSTAT=IRESP)
+  CALL ABORT
+  STOP
+END IF
+JPOS_A1T2 = 1
+DO JLOOP2_A1T2 = 1, KL2
+  DO JLOOP1_A1T2 = 1, KL1
+    P2(JLOOP1_A1T2,JLOOP2_A1T2) = P1(JPOS_A1T2)
+    JPOS_A1T2 = JPOS_A1T2 + 1
+  END DO
+END DO
+END SUBROUTINE ARRAY_1D_TO_2D
+!
+END SUBROUTINE READ_CHEM_DATA_NETCDF_CASE
diff --git a/src/MNH/shallow_mf_pack.f90 b/src/MNH/shallow_mf_pack.f90
index 421a54d6b3598337e2c498090c7bedbd62f2d385..a91a1d9722befa1e1af18c90a00d0707c8783ca8 100644
--- a/src/MNH/shallow_mf_pack.f90
+++ b/src/MNH/shallow_mf_pack.f90
@@ -114,6 +114,7 @@ END MODULE MODI_SHALLOW_MF_PACK
 !!      V.Masson 09/2010
 !!      Modification R. Honnert 07/2012 : introduction of vertical wind 
 !!                                        for the height of the thermal
+!!                   M. Leriche 02/2017 : avoid negative values for sv tendencies
 !! --------------------------------------------------------------------------
 !
 !*      0. DECLARATIONS
@@ -365,8 +366,8 @@ PRVS(:,:,:)   = PRVS(:,:,:)  +MYM(  &
 
 DO JSV=1,ISV 
   IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
-  PRSVS(:,:,:,JSV)   = PRSVS(:,:,:,JSV)  +    &
-                  PRHODJ(:,:,:)*ZDSVDT(:,:,:,JSV)
+  PRSVS(:,:,:,JSV)   = MAX((PRSVS(:,:,:,JSV)  +    &
+                  PRHODJ(:,:,:)*ZDSVDT(:,:,:,JSV)),XSVMIN(JSV))
 END DO     
 
 !!! 7. call to MesoNH budgets
diff --git a/src/MNH/turb_ver_sv_flux.f90 b/src/MNH/turb_ver_sv_flux.f90
index c3ec5d032a87c555859d442c25e112aae11d74f0..35aef367a7328d33656d5d733539d57819e2a331 100644
--- a/src/MNH/turb_ver_sv_flux.f90
+++ b/src/MNH/turb_ver_sv_flux.f90
@@ -267,6 +267,9 @@ END MODULE MODI_TURB_VER_SV_FLUX
 !!                                              change of YCOMMENT
 !!                     Feb 2012(Y. Seity) add possibility to run with reversed 
 !!                                              vertical levels
+!!                     Feb 2017(M. Leriche) add initialisation of ZSOURCE
+!!                                   to avoid unknwon values outside physical domain
+!!                                   and avoid negative values in sv tendencies
 !!--------------------------------------------------------------------------
 !       
 !*      0. DECLARATIONS
@@ -277,7 +280,7 @@ USE MODD_CTURB
 USE MODD_PARAMETERS
 USE MODD_LES
 USE MODD_CONF
-USE MODD_NSV, ONLY : NSV_LGBEG,NSV_LGEND
+USE MODD_NSV, ONLY : XSVMIN,NSV_LGBEG,NSV_LGEND
 !
 USE MODI_GRADIENT_U
 USE MODI_GRADIENT_V
@@ -397,6 +400,7 @@ DO JSV=1,ISV
   ZA(:,:,:)    = -PTSTEP*XCHF*PPSI_SV(:,:,:,JSV) *   &
                  ZKEFF * MZM(KKA,KKU,KKL,PRHODJ) /   &
                  PDZZ**2
+  ZSOURCE(:,:,:) = 0.
 !
 ! Compute the sources for the JSVth scalar variable
 
@@ -421,8 +425,8 @@ DO JSV=1,ISV
   CALL TRIDIAG(KKA,KKU,KKL,PSVM(:,:,:,JSV),ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,ZSOURCE,ZRES)
 !
 !  Compute the equivalent tendency for the JSV scalar variable
-  PRSVS(:,:,:,JSV)= PRSVS(:,:,:,JSV)+    &
-                    PRHODJ(:,:,:)*(ZRES(:,:,:)-PSVM(:,:,:,JSV))/PTSTEP
+  PRSVS(:,:,:,JSV)= MAX((PRSVS(:,:,:,JSV)+    &
+                    PRHODJ(:,:,:)*(ZRES(:,:,:)-PSVM(:,:,:,JSV))/PTSTEP),XSVMIN(JSV))
 !
   IF ( (OTURB_FLX .AND. OCLOSE_OUT) .OR. LLES_CALL ) THEN
     ! Diagnostic of the cartesian vertical flux
diff --git a/src/MNH/ver_prep_netcdf_case.f90 b/src/MNH/ver_prep_netcdf_case.f90
new file mode 100644
index 0000000000000000000000000000000000000000..097f009c5498ca5cebe9e64755f2f8f92382ad07
--- /dev/null
+++ b/src/MNH/ver_prep_netcdf_case.f90
@@ -0,0 +1,211 @@
+!     ################################
+      MODULE MODI_VER_PREP_NETCDF_CASE
+!     ################################
+INTERFACE
+      SUBROUTINE VER_PREP_NETCDF_CASE(PDIAG)
+!
+REAL, INTENT(OUT)                 :: PDIAG    ! diagnostics computing time
+!
+END SUBROUTINE VER_PREP_NETCDF_CASE
+END INTERFACE
+END MODULE MODI_VER_PREP_NETCDF_CASE
+!     ####################################################################
+      SUBROUTINE VER_PREP_NETCDF_CASE(PDIAG)
+!     ####################################################################
+!
+!!****  *VER_PREP_NETCDF_CASE* - monitors the preparation to orographic change
+!!
+!!    PURPOSE
+!!    -------
+!!    This routine monitors the preparation of variables to future change
+!!    of orography, according to the type of input file.
+!!
+!!**  METHOD
+!!    ------
+!!
+!!
+!!    EXTERNAL
+!!    --------
+!!
+!!    function MZF
+!!    function FMLOOK  :to retrieve a logical unit number associated with a file
+!!    routine VER_INTERP_TO_MIXED_GRID
+!!    routine CHANGE_GRIBEX_VAR
+!!
+!!    module MODI_SHUMAN
+!!    module MODI_VER_INTERP_TO_MIXED_GRID
+!!    module MODI_CHANGE_GRIBEX_VAR
+!!
+!!    IMPLICIT ARGUMENTS
+!!    ------------------
+!!
+!!      Module MODD_CONF1     : contains configuration variables for all models.
+!!         NVERB      : verbosity level for output-listing
+!!      Module MODD_LUNIT     :  contains logical unit names for all models
+!!         CLUOUT0 : name of output-listing
+!!      Module MODD_CST       : contains physical constants
+!!         XRD : gas constant for dry air
+!!         XRV : gas constant for vapor
+!!         XP00: reference pressure
+!!         XCPD: specific heat for dry air
+!!         XG  : gravity constant
+!!         XRADIUS : earth radius
+!!
+!!    REFERENCE
+!!    ---------
+!!
+!!      Book 2
+!!
+!!    AUTHOR
+!!    ------
+!!
+!!      V.Masson  Meteo-France
+!!
+!!    MODIFICATIONS
+!!    -------------
+!!      Original    14/12/94
+!!                  Jan, 31 1996 (V. Masson) duplication of the routine
+!!                               to accept different input fields
+!!                  May, 25 1996 (V. Masson) take into account the upper level
+!!                  Aug, 20 1996 (V. Masson) correction on theta
+!!                  Oct, 20 1996 (V. Masson) add deallocations
+!!                  Dec, 06 1996 (V. Masson) add air temperature at ground
+!!                  Dec, 12 1996 (V. Masson) add vertical wind velocity
+!!                  May, 07 1997 (V. Masson) add null tke
+!!                  Jun, 10 1997 (V. Masson) add null difference between
+!!                                           pressure and hydrostatic pressure
+!!                  Jul, 11 1997 (V. Masson) add null scalar variables
+!!                  Nov, 22 2000 (I. Mallet) add scalar variables
+!!                  Nov, 22 2000 (P. Jabouille) change routine name
+!!                  May 2006                 Remove EPS
+!-------------------------------------------------------------------------------
+!
+!*       0.    DECLARATIONS
+!              ------------
+!
+USE MODE_THERMO
+USE MODE_FM
+!
+USE MODI_SHUMAN         ! interface modules
+USE MODI_CHANGE_GRIBEX_VAR
+USE MODI_VER_INTERP_TO_MIXED_GRID
+USE MODI_RMS_AT_Z
+USE MODI_COMPUTE_EXNER_FROM_TOP
+USE MODI_WATER_SUM
+!
+USE MODD_CONF           ! declaration modules
+USE MODD_CONF_n
+USE MODD_LUNIT
+USE MODD_CST
+USE MODD_PREP_REAL
+USE MODD_PARAMETERS, ONLY : JPVEXT, XUNDEF
+!
+IMPLICIT NONE
+!
+!*       0.1   Declaration of arguments
+!              ------------------------
+!
+REAL, INTENT(OUT)                 :: PDIAG    ! diagnostics computing time
+!
+!*       0.2   Declaration of local variables
+!              ------------------------------
+INTEGER                            :: IRESP, ILUOUT0
+INTEGER                            :: IIU,IJU,ILU
+REAL                               :: ZTIME1, ZTIME2
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZTH_LS    ! potential temperature
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZTH_MX    ! potential temperature
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZPMASS_MX    ! pressure
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZHEXNFLUX_MX ! pressure function
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZHEXNMASS_MX ! pressure function
+!
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZZFLUX_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZZMASS_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZPMHP_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZTHV_LS
+REAL,DIMENSION(:,:,:,:),ALLOCATABLE:: ZR_LS
+REAL,DIMENSION(:,:,:,:),ALLOCATABLE:: ZSV_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZHU_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZU_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZV_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZW_LS
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZTKE_LS
+INTEGER                            :: JRR     ! loop counter
+INTEGER                            :: JSV     ! loop counter
+INTEGER                            :: JK      ! loop counter
+!-------------------------------------------------------------------------------
+!
+CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
+!
+!*       1.    CHANGING OF VARIABLES
+!              ---------------------
+!
+  IIU=SIZE(XT_SV_LS,1)
+  IJU=SIZE(XT_SV_LS,2)
+  ILU=SIZE(XT_SV_LS,3)
+!
+!
+  ALLOCATE(XPMASS_SV_LS(IIU,IJU,ILU))
+  ALLOCATE(XZMASS_SV_LS(IIU,IJU,ILU),XZFLUX_SV_LS(IIU,IJU,ILU))
+  ALLOCATE(XTHV_SV_LS(IIU,IJU,ILU),XR_SV_LS(IIU,IJU,ILU,NRR),XHU_SV_LS(IIU,IJU,ILU))
+  CALL CHANGE_GRIBEX_VAR(XA_SV_LS,XB_SV_LS,XP00_SV_LS,XPS_SV_LS,XZS_SV_LS, &
+                         XT_SV_LS,XQ_SV_LS,XPMASS_SV_LS,XZFLUX_SV_LS,XZMASS_SV_LS, &
+                         XTHV_SV_LS,XR_SV_LS,XHU_SV_LS                      )
+!
+!-------------------------------------------------------------------------------
+!
+!*       2.    INTERPOLATION TO MIXED GRID AND DIAGNOSTIC VARIABLES
+!              ----------------------------------------------------
+!* Add extra points below and above grids, in order to use MESONH linear
+! vertical interpolation programs with all ILU physical points
+!
+ALLOCATE(ZZMASS_LS(IIU,IJU,ILU+2*JPVEXT))
+ALLOCATE(ZSV_LS(IIU,IJU,ILU+2*JPVEXT,SIZE(XSV_LS,4)))
+!
+ZZMASS_LS (:,:,JPVEXT+1:JPVEXT+ILU) = XZMASS_SV_LS(:,:,:)
+DO JK=1,JPVEXT
+   ZZMASS_LS(:,:,           JK) = XZMASS_SV_LS(:,:,1)   - (XZMASS_SV_LS(:,:,2)  -XZMASS_SV_LS(:,:,1)    )*(JPVEXT+1-JK)
+   ZZMASS_LS(:,:,ILU+JPVEXT+JK) = XZMASS_SV_LS(:,:,ILU) + (XZMASS_SV_LS(:,:,ILU)-XZMASS_SV_LS(:,:,ILU-1))*          JK
+END DO
+!
+!ZSV_LS    = XUNDEF
+ZSV_LS    = -999.
+!
+DO JSV=1,SIZE(XSV_LS,4)
+  ZSV_LS    (:,:,JPVEXT+1:JPVEXT+ILU,JSV) = XSV_LS (:,:,:,JSV)
+END DO
+!
+  CALL VER_INTERP_TO_MIXED_GRID('CHEM',.TRUE.,XZS_SV_LS,XZS_SV_LS,&
+                                ZZMASS_LS,ZSV_LS                 )
+!
+DEALLOCATE(ZZMASS_LS)
+DEALLOCATE(ZSV_LS)
+!-------------------------------------------------------------------------------
+!
+!*       3.    ERROR CONTROL
+!              -------------
+!
+CALL SECOND_MNH(ZTIME1)
+PDIAG = ZTIME2 - ZTIME1
+!
+!-------------------------------------------------------------------------------
+!
+!*       4.    DEALLOCATIONS
+!              -------------
+!
+  DEALLOCATE(XA_SV_LS)
+  DEALLOCATE(XB_SV_LS)
+  DEALLOCATE(XT_SV_LS)
+  DEALLOCATE(XQ_SV_LS)
+  DEALLOCATE(XZMASS_SV_LS)
+  DEALLOCATE(XZFLUX_SV_LS)
+  DEALLOCATE(XTHV_SV_LS)
+  DEALLOCATE(XR_SV_LS)
+  DEALLOCATE(XHU_SV_LS)
+  DEALLOCATE(XSV_LS)
+!
+!
+!-------------------------------------------------------------------------------
+!
+WRITE(ILUOUT0,*) 'Routine VER_PREP_NETCDF_CASE completed'
+!
+END SUBROUTINE VER_PREP_NETCDF_CASE