diff --git a/MY_RUN/KTEST/009_ICARTT/003_mesonh/EXSEG1.nam.src b/MY_RUN/KTEST/009_ICARTT/003_mesonh/EXSEG1.nam.src
index 501672e664980d2b5e16751d7f21166033fadcf0..60fd1369bab362a593fddbc76eb19fec0a1145db 100644
--- a/MY_RUN/KTEST/009_ICARTT/003_mesonh/EXSEG1.nam.src
+++ b/MY_RUN/KTEST/009_ICARTT/003_mesonh/EXSEG1.nam.src
@@ -22,7 +22,6 @@
 &NAM_CH_MNHCn       LUSECHEM = T,
                     LCH_CONV_LINOX = T,
                     LCH_INIT_FIELD = F,
-                    LCH_SURFACE_FLUX = T,
                     LCH_CONV_SCAV = T,
                     CCHEM_INPUT_FILE  = "ReLACS_poet.nam",
                     CCH_TDISCRETIZATION = "SPLIT"
diff --git a/MY_RUN/KTEST/011_KW78CHEM/002_mesonh/EXSEG1.nam b/MY_RUN/KTEST/011_KW78CHEM/002_mesonh/EXSEG1.nam
index 7aef871115a3f6ad509b555952fde79b79ef25ea..ae748b4f8898f5dbd7da3066fbfe1943478edc4e 100644
--- a/MY_RUN/KTEST/011_KW78CHEM/002_mesonh/EXSEG1.nam
+++ b/MY_RUN/KTEST/011_KW78CHEM/002_mesonh/EXSEG1.nam
@@ -14,7 +14,6 @@
                     LUSECHIC = T,
                     LCH_CONV_LINOX = T,
                     LCH_INIT_FIELD = T,
-                    LCH_SURFACE_FLUX = F,
                     LCH_CONV_SCAV = T,
                     CCHEM_INPUT_FILE  = "MNHC.input",
                     LCH_RET_ICE = F,
diff --git a/src/MNH/ch_aer_driver.f90 b/src/MNH/ch_aer_driver.f90
index dc81a8e5928d7d52649dce6054a89b73faa668e1..819eb0ae1dc7a596adf413576e4949e517908041 100644
--- a/src/MNH/ch_aer_driver.f90
+++ b/src/MNH/ch_aer_driver.f90
@@ -56,6 +56,7 @@ SUBROUTINE CH_AER_DRIVER(PM, PSIG0, PRG0, PN0, PCTOTG, PCTOTA,&
 !!    -------------
 !!    Original
 !!       M.Leriche 2015 Calcul de la fraction massique entre les modes
+!!       M.Leriche 08/16 suppress moments index declaration already in modd_aerosol
 !!
 !!    EXTERNAL
 !!    --------
@@ -103,13 +104,6 @@ REAL, DIMENSION(SIZE(PM,1))   :: ZPKM, ZPKH2O, ZSUM
 !-----------------------------------------------------------------------------
 
 ZDT=PDTACT
-! Moments index
-NM0(1) = 1
-NM3(1) = 2
-NM6(1) = 3
-NM0(2) = 4
-NM3(2) = 5
-NM6(2) = 6
 
 !*************************************************************
 ! Calcul de la fraction massique entre les modes
diff --git a/src/MNH/ch_aer_solv.f90 b/src/MNH/ch_aer_solv.f90
index 6bfd75ce3fef93ae5c9bb8151cf522cbc86660fe..22d84df38ff5c22925e18d6d71de7c0c3dc6f97f 100644
--- a/src/MNH/ch_aer_solv.f90
+++ b/src/MNH/ch_aer_solv.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$ $Date$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/ch_aer_solv.f90,v $ $Revision: 1.1.2.1.2.1.16.2.2.1.2.1 $ $Date: 2015/12/01 15:26:23 $
 !-----------------------------------------------------------------
 !-----------------------------------------------------------------
 !!   #######################
@@ -62,10 +62,14 @@ END MODULE MODI_CH_AER_SOLV
 !!    P. Tulet organic condensation
 !!    P. Tulet thermodynamic equilibrium for each mode
 !!    P. Tulet add third mode
+!!    M. Leriche 2015 correction bug
+!!    M. Leriche 08/16 suppress moments index declaration already in modd_aerosol
+!!    M. Leriche 08/16 add an other particular case for the M0 resolution to
+!!               avoid a division by zero (when ZK = 1)
 !!
 !!    EXTERNAL
 !!    --------
-!!    M.Leriche 2015 correction bug
+!!
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -131,16 +135,6 @@ LOGICAL, SAVE               :: GPHYSLIM = .TRUE. ! flag
   ZPMIN(5) = ZPMIN(4) * (ZINIRADIUSJ**3)*EXP(4.5 * LOG(XSIGJMIN)**2) 
   ZPMIN(6) = ZPMIN(4) * (ZINIRADIUSJ**6)*EXP(18. * LOG(XSIGJMIN)**2)
   !
-! Moments index
-NM0(1) = 1
-NM3(1) = 2
-NM6(1) = 3
-NM0(2) = 4
-NM3(2) = 5
-NM6(2) = 6
-!
-
-
 !write(*,*)
 !write(*,*) '******************************************'
 !write(*,*) '         Debut Solveur Aerosol            '
@@ -170,20 +164,22 @@ ZA(:)=PDMINTRA(:,NM0(JI))
 ZB(:)=PDMINTER(:,NM0(JI))
 ZC(:)=PDMCOND(:,NM0(JI))
 
+
 DO JK=1,SIZE(PM,1)
- IF  (ZB(JK) == 0. .AND. ZC(JK)/PM(JK,NM0(JI)) <= 1.e-10)  THEN
+ IF ((ZB(JK) == 0. .AND. ZC(JK)/PM(JK,NM0(JI)) <= 1.e-10).OR. &
+     (ZC(JK) <= 1.e-10 .AND. ZB(JK)/ZA(JK) <= 1.e-3))  THEN
+! type dY/dt=-AY^2
    Z0(JK)=PM(JK,NM0(JI)) 
    PM(JK,NM0(JI))=Z0(JK)/(1.+ZA(JK)*Z0(JK)*PDT)
  ELSE
-   ZD(JK)=SQRT(ZB(JK)**2+4.*ZA(JK)*ZC(JK))
-
    ZCONST1(JK)=ZB(JK)/(2.*ZA(JK))
-   ZCONST2(JK)=ZD(JK)/(2.*ABS(ZA(JK)))
    Z0(JK)=PM(JK,NM0(JI))+ZCONST1(JK)
-  
    IF (((ZB(JK)**2+4.*ZA(JK)*ZC(JK))) < 0.) THEN
+     ZD(JK)=SQRT(ABS(ZB(JK)**2+4.*ZA(JK)*ZC(JK)))
      PM(JK,NM0(JI))=-ZCONST1(JK)+ZD(JK)*TAN(ATAN(Z0(JK)/ZD(JK))-ZA(JK)*ZD(JK)*PDT)
    ELSE
+     ZD(JK)=SQRT(ZB(JK)**2+4.*ZA(JK)*ZC(JK))
+     ZCONST2(JK)=ZD(JK)/(2.*ABS(ZA(JK)))
      ZKEXP(JK)=EXP(-2.*ZA(JK)*ZCONST2(JK)*PDT)
      ZK(JK)=(Z0(JK)-ZCONST2(JK))/(Z0(JK)+ZCONST2(JK))*ZKEXP(JK)
      PM(JK,NM0(JI))=-ZCONST1(JK)+ZCONST2(JK)*(1.+ZK(JK))/(1.-ZK(JK))
diff --git a/src/MNH/ch_aqueous_sedim1mom.f90 b/src/MNH/ch_aqueous_sedim1mom.f90
index c9683e2636286727a2c2eb6695723c93c34a38ac..9bf9109b022306bf0ef99d6fd47b8d62230da921 100644
--- a/src/MNH/ch_aqueous_sedim1mom.f90
+++ b/src/MNH/ch_aqueous_sedim1mom.f90
@@ -7,12 +7,12 @@
 !      ################################
 !
 INTERFACE
-      SUBROUTINE CH_AQUEOUS_SEDIM1MOM (PTIME, HCLOUD, OUSECHIC, PTSTEP, &
-                                       PZZ, PRHODREF, PRHODJ, PRRS,     &
-                                       PRSS, PRGS, PRRSVS, PSGRSVS      )
+      SUBROUTINE CH_AQUEOUS_SEDIM1MOM (KSPLITR, HCLOUD, OUSECHIC, PTSTEP,  &
+                                       PZZ, PRHODREF, PRHODJ, PRRS,        &
+                                       PRSS, PRGS, PRRSVS, PSGRSVS, PINPRR )
 !
 CHARACTER (LEN=4),        INTENT(IN)    :: HCLOUD  ! Cloud parameterization
-REAL,                     INTENT(IN)    :: PTIME  ! Current time
+INTEGER,                  INTENT(IN)    :: KSPLITR ! Current time
 REAL,                     INTENT(IN)    :: PTSTEP  ! Time step          
 LOGICAL,                  INTENT(IN)    :: OUSECHIC ! flag for ice chem.
 !
@@ -24,16 +24,17 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRSS    ! Snow m.r. source
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRGS    ! Graupel m.r. source
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRSVS  ! Rainwater aq. species source
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSGRSVS ! Precip. ice species source
+REAL, DIMENSION(:,:),     INTENT(OUT)   :: PINPRR  ! instantaneaous precip.
 !
 END SUBROUTINE CH_AQUEOUS_SEDIM1MOM
 END INTERFACE
 END MODULE MODI_CH_AQUEOUS_SEDIM1MOM
 !
-!     ###################################################################
-      SUBROUTINE CH_AQUEOUS_SEDIM1MOM (PTIME, HCLOUD, OUSECHIC, PTSTEP, &
-                                       PZZ, PRHODREF, PRHODJ, PRRS,     &
-                                       PRSS, PRGS, PRRSVS, PSGRSVS      )
-!     ###################################################################
+!     ######################################################################
+      SUBROUTINE CH_AQUEOUS_SEDIM1MOM (KSPLITR, HCLOUD, OUSECHIC, PTSTEP,  &
+                                       PZZ, PRHODREF, PRHODJ, PRRS,        &
+                                       PRSS, PRGS, PRRSVS, PSGRSVS, PINPRR )
+!     ######################################################################
 !
 !!****  * -  compute the explicit microphysical sources 
 !!
@@ -77,6 +78,7 @@ END MODULE MODI_CH_AQUEOUS_SEDIM1MOM
 !!    04/11/08 (M Leriche) add ICE3    
 !!    17/09/10 (M Leriche) add LUSECHIC flag
 !!    J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!!    16/12/15 (M Leriche) compute instantaneous rain at the surface
 !!
 !-------------------------------------------------------------------------------
 !
@@ -85,6 +87,7 @@ END MODULE MODI_CH_AQUEOUS_SEDIM1MOM
 !
 USE MODD_PARAMETERS,      ONLY : JPHEXT, JPVEXT
 USE MODD_CONF
+USE MODD_CST,             ONLY : XRHOLW
 USE MODD_CLOUDPAR,        ONLY : VCEXVT=>XCEXVT, XCRS, XCEXRS
 USE MODD_RAIN_ICE_DESCR,  ONLY : WCEXVT=>XCEXVT, WRTMIN=>XRTMIN
 USE MODD_RAIN_ICE_PARAM,  ONLY : XFSEDR, XEXSEDR, &
@@ -97,7 +100,7 @@ IMPLICIT NONE
 !
 !
 CHARACTER (LEN=4),        INTENT(IN)    :: HCLOUD  ! Cloud parameterization
-REAL,                     INTENT(IN)    :: PTIME  ! Current time
+INTEGER,                  INTENT(IN)    :: KSPLITR ! Current time
 REAL,                     INTENT(IN)    :: PTSTEP  ! Time step          
 LOGICAL,                  INTENT(IN)    :: OUSECHIC ! flag for ice chem.
 !
@@ -109,6 +112,7 @@ REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRSS    ! Snow m.r. source
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRGS    ! Graupel m.r. source
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRSVS  ! Rainwater aq. species source
 REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSGRSVS ! Precip. ice species source
+REAL, DIMENSION(:,:),     INTENT(OUT)   :: PINPRR  ! instantaneaous precip.
 !
 !*       0.2   Declarations of local variables :
 !
@@ -162,9 +166,8 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZRHODREF, & ! RHO Dry REFerence
                                    ZZW         ! Work array
 REAL, DIMENSION(7), SAVE        :: Z_XRTMIN
 !
-REAL                            :: ZVTRMAX, ZDZMIN, ZT
+REAL                            :: ZVTRMAX, ZT
 LOGICAL, SAVE                   :: GSFIRSTCALL = .TRUE.
-INTEGER, SAVE                   :: ISPLITR
 REAL,    SAVE                   :: ZFSEDR, ZEXSEDR, ZCEXVT
 !
 INTEGER , DIMENSION(SIZE(GSEDIMR)) :: IR1,IR2,IR3 ! Used to replace the COUNT
@@ -180,6 +183,7 @@ INTEGER                           :: JL       ! and PACK intrinsics
 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
 IKB=1+JPVEXT
 IKE=SIZE(PZZ,3) - JPVEXT
+PINPRR(:,:) = 0. ! initialize instantaneous precip.
 !
 !-------------------------------------------------------------------------------
 !
@@ -197,7 +201,7 @@ ENDIF
 !*       3.     COMPUTE THE SEDIMENTATION (RS) SOURCE
 !	        -------------------------------------
 !
-!*       3.1    splitting factor for high Courant number C=v_fall*(del_Z/del_T)
+!*       3.1    Initialize some constants
 !  
 firstcall : IF (GSFIRSTCALL) THEN
   GSFIRSTCALL = .FALSE.
@@ -209,13 +213,6 @@ firstcall : IF (GSFIRSTCALL) THEN
     CASE('ICE4')
       ZVTRMAX = 40.
   END SELECT
-  ZDZMIN = MINVAL(PZZ(IIB:IIE,IJB:IJE,IKB+1:IKE+1)-PZZ(IIB:IIE,IJB:IJE,IKB:IKE))
-  ISPLITR = 1
-  SPLIT : DO
-    ZT = PTSTEP / FLOAT(ISPLITR)
-    IF ( ZT * ZVTRMAX / ZDZMIN .LT. 1.) EXIT SPLIT
-    ISPLITR = ISPLITR + 1
-  END DO SPLIT
 !
   SELECT CASE ( HCLOUD )  ! constants for rain sedimentation
     CASE('KESS')
@@ -233,7 +230,7 @@ END IF firstcall
 !
 !*       3.2    time splitting loop initialization
 !
-ZTSPLITR = PTSTEP / FLOAT(ISPLITR)       ! Small time step
+ZTSPLITR = PTSTEP / FLOAT(KSPLITR)       ! Small time step
 !
 !*       3.3    compute the fluxes
 !
@@ -245,7 +242,7 @@ IF (HCLOUD(1:3) == 'ICE') THEN
   ZSV_SEDIM_FACTS(:,:,:) = 1.0
   ZSV_SEDIM_FACTG(:,:,:) = 1.0
 ENDIF
-DO JN = 1 , ISPLITR
+DO JN = 1 , KSPLITR
   IF( JN==1 ) THEN
     ZW(:,:,:) = 0.0
     DO JK = IKB , IKE-1
@@ -275,6 +272,7 @@ DO JN = 1 , ISPLITR
       ZRR_SEDIM(:,:,JK) = ZW(:,:,JK)*(ZWSED(:,:,JK+1)-ZWSED(:,:,JK))
     END DO
     ZZRRS(:,:,:) = ZZRRS(:,:,:) + ZRR_SEDIM(:,:,:)
+    PINPRR(:,:) = PINPRR(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
 !
     ZZW(:) = ZFSEDR * ZZZRRS(:)**(ZEXSEDR-1.0) * ZRHODREF(:)**(ZEXSEDR-ZCEXVT)
     ZWSED(:,:,:) = UNPACK( ZZW(:),MASK=GSEDIMR(:,:,:),FIELD=0.0 )
diff --git a/src/MNH/ch_aqueous_tmicice.f90 b/src/MNH/ch_aqueous_tmicice.f90
index 042fefda1364c47dda471c50daceae2a442677aa..969909f6082f338238d6c926bb15eb3b7becdc6b 100644
--- a/src/MNH/ch_aqueous_tmicice.f90
+++ b/src/MNH/ch_aqueous_tmicice.f90
@@ -128,6 +128,7 @@ USE MODD_RAIN_ICE_PARAM,  ONLY : XTIMAUTC, XCRIAUTC, XFCACCR, XEXCACCR, &
                                  XKER_RDRYG, XLBRDRYG1, XLBRDRYG2, XLBRDRYG3,   &
                                  XCOLIG, XCOLEXIG, XCOLSG, XCOLEXSG
 USE MODD_CH_ICE                              ! value of retention coefficient
+USE MODD_CH_ICE_n                            ! index for ice phase chemistry with IC3/4
 !
 #ifdef MNH_PGI
 USE MODE_PACK_PGI
@@ -185,9 +186,6 @@ INTEGER :: IKB
 INTEGER :: IKE           
 !
 INTEGER :: IMICRO        ! case number of r_x>0 locations
-INTEGER, DIMENSION(SIZE(PSGRSVS,4))  ::  INDEXGI  ! index array for ice phase chemistry
-INTEGER, DIMENSION(SIZE(PSGRSVS,4))  ::  INDEXWI  ! index array for ice phase chemistry
-INTEGER, DIMENSION(SIZE(PRRSVS,4))   ::  INDEXWG  ! index array for degassing when freezing
 LOGICAL, DIMENSION(SIZE(PRCT,1),SIZE(PRCT,2),SIZE(PRCT,3))   &
                                 :: GMICRO   ! where to compute mic. processes
 REAL,    DIMENSION(SIZE(PRCT,1),SIZE(PRCT,2),SIZE(PRCT,3))   &
@@ -455,50 +453,10 @@ IF( IMICRO >= 1 ) THEN
 !
 !-------------------------------------------------------------------------------
 !
-!*       5.     PREPARE INDEX ARRAY FOR ICE PHASE CHEMISTRY
-!               -------------------------------------------
-!
-IF (OUSECHIC) THEN
-  DO JLI = 1, SIZE(PSGRSVS,4)
-    DO JLG = 1, SIZE(PGRSVS,4)
-      IF ( TRIM(HICNAMES(JLI)(4:32)) == TRIM(HNAMES(JLG)) ) THEN
-         INDEXGI(JLI) = JLG
-         EXIT
-      ELSE
-         INDEXGI(JLI) = 0
-      ENDIF
-    ENDDO
-    DO JLW = KEQ-KEQAQ+1, KEQ-KEQAQ/2  ! loop over cloud chem. species
-      IF ( TRIM(HICNAMES(JLI)(4:32)) == TRIM(HNAMES(JLW)(4:32))) THEN
-        INDEXWI(JLI) = JLW - (KEQ-KEQAQ)
-        EXIT
-      ELSE
-        INDEXWI(JLI) = 0
-      ENDIF
-    ENDDO
-  ENDDO
-ELSE
-  IF (.NOT.(OCH_RET_ICE)) THEN
-    DO JLW = KEQ-KEQAQ+1, KEQ-KEQAQ/2  ! loop over cloud chem. species
-      DO JLG = 1, SIZE(PGRSVS,4)
-        IF ( TRIM(HNAMES(JLW)(4:32)) == TRIM(HNAMES(JLG)) ) THEN
-          INDEXWG(JLW-(KEQ-KEQAQ)) = JLG
-          EXIT
-        ELSE
-          INDEXWG(JLW-(KEQ-KEQAQ)) = 0
-        ENDIF
-      ENDDO
-    ENDDO
-  ENDIF
-ENDIF
-!
-!
-!-------------------------------------------------------------------------------
-!
-!*       6.     COMPUTES THE SLOW COLD PROCESS SOURCES
+!*       5.     COMPUTES THE SLOW COLD PROCESS SOURCES
 !               --------------------------------------
 !
-!*       6.1    compute the spontaneous freezing source: RRHONG
+!*       5.1    compute the spontaneous freezing source: RRHONG
 !
  ZZW(:) = 0.0
  ZZW2(:,:) = 0.0
@@ -512,33 +470,35 @@ ENDIF
      IF (OUSECHIC) THEN
      DO JLI = 1, SIZE(PSGRSVS,4)
        IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-          .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))           
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+          .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+          .OR. NINDEXGI(JLI).EQ.0) THEN
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))           
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))  
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))  
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
        ELSE
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
        ENDIF
      ENDDO
      ELSE
        IF (.NOT.(OCH_RET_ICE)) THEN
          DO JLW = 1, SIZE(PRRSVS,4)
-            IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-              ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+            IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+              ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
             ENDIF
          ENDDO
        ENDIF
@@ -549,10 +509,10 @@ ENDIF
 !
 !-------------------------------------------------------------------------------
 !
-!*       7.     COMPUTES THE FAST COLD PROCESS SOURCES
+!*       6.     COMPUTES THE FAST COLD PROCESS SOURCES
 !               --------------------------------------
 !
-!*       7.1    compute the slope parameter Lbda_s and Lbda_g
+!*       6.1    compute the slope parameter Lbda_s and Lbda_g
 !
  WHERE ( ZRST(:)>0.0 )
    ZLBDAS(:)  = MIN( XLBDAS_MAX,                                           &
@@ -563,7 +523,7 @@ ENDIF
    ZLBDAG(:)  = XLBG*( ZRHODREF(:)*MAX( ZRGT(:),PRTMIN_AQ*1.e3/ZRHODREF(:)))**XLBEXG
  END WHERE
 !
-!*       7.2    cloud droplet riming of the aggregates
+!*       6.2    cloud droplet riming of the aggregates
 !
  ZZW1(:,:) = 0.0
  ZZW(:) = 0.0
@@ -576,18 +536,18 @@ ENDIF
 !
  IF( IGRIM>0 ) THEN
 !
-!        7.2.0  allocations
+!        6.2.0  allocations
 !
    ALLOCATE(ZVEC1(IGRIM))
    ALLOCATE(ZVEC2(IGRIM))
    ALLOCATE(IVEC1(IGRIM))
    ALLOCATE(IVEC2(IGRIM))
 !
-!        7.2.1  select the ZLBDAS
+!        6.2.1  select the ZLBDAS
 !
    ZVEC1(:) = PACK( ZLBDAS(:),MASK=GRIM(:) )
 !
-!        7.2.2  find the next lower indice for the ZLBDAS in the geometrical
+!        6.2.2  find the next lower indice for the ZLBDAS in the geometrical
 !               set of Lbda_s used to tabulate some moments of the incomplete
 !               gamma function
 !
@@ -596,14 +556,14 @@ ENDIF
    IVEC2(1:IGRIM) = INT( ZVEC2(1:IGRIM) )
    ZVEC2(1:IGRIM) = ZVEC2(1:IGRIM) - FLOAT( IVEC2(1:IGRIM) )
 !
-!        7.2.3  perform the linear interpolation of the normalized
+!        6.2.3  perform the linear interpolation of the normalized
 !               "2+XDS"-moment of the incomplete gamma function
 !
    ZVEC1(1:IGRIM) =   XGAMINC_RIM1( IVEC2(1:IGRIM)+1 )* ZVEC2(1:IGRIM)      &
                     - XGAMINC_RIM1( IVEC2(1:IGRIM)   )*(ZVEC2(1:IGRIM) - 1.0)
    ZZW(:) = UNPACK( VECTOR=ZVEC1(:),MASK=GRIM,FIELD=0.0 )
 !
-!        7.2.4  riming of the small sized aggregates
+!        6.2.4  riming of the small sized aggregates
 !
    ZZW2(:,:) = 0.0
    DO JL = 1,IMICRO
@@ -616,33 +576,35 @@ ENDIF
        IF (OUSECHIC) THEN
        DO JLI = 1, SIZE(PSGRSVS,4)
          IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-            .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+            .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+            .OR. NINDEXGI(JLI).EQ.0) THEN
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                  TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
          ELSE
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
          ENDIF
        ENDDO
        ELSE
          IF (.NOT.(OCH_RET_ICE)) THEN
            DO JLW = 1, SIZE(PCRSVS,4)
-             IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-               ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+             IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+               ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
              ENDIF
            ENDDO
          ENDIF
@@ -650,7 +612,7 @@ ENDIF
      ENDIF
    ENDDO
 !
-!        7.2.5  riming-conversion of the large sized aggregates into graupel
+!        6.2.5  riming-conversion of the large sized aggregates into graupel
 !
    ZZW2(:,:) = 0.0
    DO JL = 1,IMICRO
@@ -663,33 +625,35 @@ ENDIF
        IF (OUSECHIC) THEN
        DO JLI = 1, SIZE(PSGRSVS,4)
          IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-            .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+            .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+            .OR. NINDEXGI(JLI).EQ.0) THEN
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                  TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
          ELSE
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
          ENDIF
        ENDDO
        ELSE
          IF (.NOT.(OCH_RET_ICE)) THEN
            DO JLW = 1, SIZE(PCRSVS,4)
-             IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-               ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+             IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+               ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
              ENDIF
            ENDDO
          ENDIF
@@ -704,7 +668,7 @@ ENDIF
  END IF
  DEALLOCATE(GRIM)
 !
-!*       7.3    rain accretion onto the aggregates
+!*       6.3    rain accretion onto the aggregates
 !
  ZZW(:) = 0.0
  ZZW1(:,2:3) = 0.0
@@ -716,7 +680,7 @@ ENDIF
 !
  IF( IGACC>0 ) THEN
 !
-!        7.3.0  allocations
+!        6.3.0  allocations
 !
    ALLOCATE(ZVEC1(IGACC))
    ALLOCATE(ZVEC2(IGACC))
@@ -724,12 +688,12 @@ ENDIF
    ALLOCATE(IVEC1(IGACC))
    ALLOCATE(IVEC2(IGACC))
 !
-!        7.3.1  select the (ZLBDAS,ZLBDAR) couplet
+!        6.3.1  select the (ZLBDAS,ZLBDAR) couplet
 !
    ZVEC1(:) = PACK( ZLBDAS(:),MASK=GACC(:) )
    ZVEC2(:) = PACK( ZLBDAR(:),MASK=GACC(:) )
 !
-!        7.3.2  find the next lower indice for the ZLBDAS and for the ZLBDAR
+!        6.3.2  find the next lower indice for the ZLBDAS and for the ZLBDAR
 !               in the geometrical set of (Lbda_s,Lbda_r) couplet use to
 !               tabulate the RACCSS-kernel
 !
@@ -743,7 +707,7 @@ ENDIF
    IVEC2(1:IGACC) = INT( ZVEC2(1:IGACC) )
    ZVEC2(1:IGACC) = ZVEC2(1:IGACC) - FLOAT( IVEC2(1:IGACC) )
 !
-!        7.3.3  perform the bilinear interpolation of the normalized
+!        6.3.3  perform the bilinear interpolation of the normalized
 !               RACCSS-kernel
 !
    DO JJ = 1,IGACC
@@ -756,7 +720,7 @@ ENDIF
    END DO
    ZZW(:) = UNPACK( VECTOR=ZVEC3(:),MASK=GACC,FIELD=0.0 )
 !
-!        7.3.4  raindrop accretion on the small sized aggregates
+!        6.3.4  raindrop accretion on the small sized aggregates
 !
    ZZW2(:,:) = 0.0
    DO JL = 1,IMICRO
@@ -773,33 +737,35 @@ ENDIF
        IF (OUSECHIC) THEN
        DO JLI = 1, SIZE(PSGRSVS,4)
          IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-            .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+            .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+            .OR. NINDEXGI(JLI).EQ.0) THEN
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                  TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
          ELSE
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
          ENDIF
        ENDDO
        ELSE
          IF (.NOT.(OCH_RET_ICE)) THEN
            DO JLW = 1, SIZE(PRRSVS,4)
-             IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-               ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+             IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+               ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
              ENDIF
            ENDDO
          ENDIF
@@ -807,7 +773,7 @@ ENDIF
      ENDIF
    ENDDO
 !
-!        7.3.4b perform the bilinear interpolation of the normalized
+!        6.3.4b perform the bilinear interpolation of the normalized
 !               RACCS-kernel
 !
    DO JJ = 1,IGACC
@@ -820,7 +786,7 @@ ENDIF
    END DO
    ZZW1(:,2) = ZZW1(:,2)*UNPACK( VECTOR=ZVEC3(:),MASK=GACC(:),FIELD=0.0 )
 !
-!        7.3.5  raindrop accretion-conversion of the large sized aggregates
+!        6.3.5  raindrop accretion-conversion of the large sized aggregates
 !               into graupeln
 !
    ZZW2(:,:) = 0.0
@@ -835,33 +801,35 @@ ENDIF
        IF (OUSECHIC) THEN
        DO JLI = 1, SIZE(PSGRSVS,4)
          IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-            .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+            .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+            .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+            .OR. NINDEXGI(JLI).EQ.0) THEN
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
             .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
          ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
             .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                  TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
          ELSE
-           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-           ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                     (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+           ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+           ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                     (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
          ENDIF
        ENDDO
        ELSE
          IF (.NOT.(OCH_RET_ICE)) THEN
            DO JLW = 1, SIZE(PRRSVS,4)
-             IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-               ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+             IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+               ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
              ENDIF
            ENDDO
          ENDIF
@@ -877,7 +845,7 @@ ENDIF
  END IF
  DEALLOCATE(GACC)
 !
-!*       7.4    rain contact freezing
+!*       6.4    rain contact freezing
 !
  ZZW1(:,4) = 0.0
  ZZW2(:,:) = 0.0
@@ -894,33 +862,35 @@ ENDIF
      IF (OUSECHIC) THEN
      DO JLI = 1, SIZE(PSGRSVS,4)
        IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-          .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+          .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+          .OR. NINDEXGI(JLI).EQ.0) THEN
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
        ELSE
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
        ENDIF
      ENDDO
      ELSE
        IF (.NOT.(OCH_RET_ICE)) THEN
          DO JLW = 1, SIZE(PRRSVS,4)
-           IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-             ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+           IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+             ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
            ENDIF
          ENDDO
        ENDIF
@@ -928,7 +898,7 @@ ENDIF
    ENDIF
  ENDDO
 !
-!*       7.5    compute the Dry growth case of graupel
+!*       6.5    compute the Dry growth case of graupel
 !
  ZZW(:) = 0.0
  ZZW1(:,:) = 0.0
@@ -944,7 +914,7 @@ ENDIF
                                     * ZRIT(:) * ZZW(:) )             ! RIDRYG
  END WHERE
 !
-!        7.5.1  accretion of aggregates on the graupeln
+!        6.5.1  accretion of aggregates on the graupeln
 !
  ALLOCATE(GDRY(IMICRO))
  GDRY(:) = (ZRST(:)>PRTMIN_AQ*1.e3/ZRHODREF(:)) .AND. &
@@ -953,7 +923,7 @@ ENDIF
 !
  IF( IGDRY>0 ) THEN
 !
-!        7.5.2  allocations
+!        6.5.2  allocations
 !
    ALLOCATE(ZVEC1(IGDRY))
    ALLOCATE(ZVEC2(IGDRY))
@@ -961,12 +931,12 @@ ENDIF
    ALLOCATE(IVEC1(IGDRY))
    ALLOCATE(IVEC2(IGDRY))
 !
-!        7.5.3  select the (ZLBDAG,ZLBDAS) couplet
+!        6.5.3  select the (ZLBDAG,ZLBDAS) couplet
 !
    ZVEC1(:) = PACK( ZLBDAG(:),MASK=GDRY(:) )
    ZVEC2(:) = PACK( ZLBDAS(:),MASK=GDRY(:) )
 !
-!        7.5.4  find the next lower indice for the ZLBDAG and for the ZLBDAS
+!        6.5.4  find the next lower indice for the ZLBDAG and for the ZLBDAS
 !               in the geometrical set of (Lbda_g,Lbda_s) couplet use to
 !               tabulate the SDRYG-kernel
 !
@@ -980,7 +950,7 @@ ENDIF
    IVEC2(1:IGDRY) = INT( ZVEC2(1:IGDRY) )
    ZVEC2(1:IGDRY) = ZVEC2(1:IGDRY) - FLOAT( IVEC2(1:IGDRY) )
 !
-!        7.5.5  perform the bilinear interpolation of the normalized
+!        6.5.5  perform the bilinear interpolation of the normalized
 !               SDRYG-kernel
 !
    DO JJ = 1,IGDRY
@@ -1009,7 +979,7 @@ ENDIF
    DEALLOCATE(ZVEC1)
  END IF
 !
-!        7.5.6  accretion of raindrops on the graupeln
+!        6.5.6  accretion of raindrops on the graupeln
 !
  GDRY(:) = (ZRRT(:)>PRTMIN_AQ*1.e3/ZRHODREF(:)) .AND. &
            (ZRGT(:)>PRTMIN_AQ*1.e3/ZRHODREF(:)) .AND. (ZZRRS(:)>0.0)
@@ -1017,7 +987,7 @@ ENDIF
 !
  IF( IGDRY>0 ) THEN
 !
-!        7.5.7  allocations
+!        6.5.7  allocations
 !
    ALLOCATE(ZVEC1(IGDRY))
    ALLOCATE(ZVEC2(IGDRY))
@@ -1025,12 +995,12 @@ ENDIF
    ALLOCATE(IVEC1(IGDRY))
    ALLOCATE(IVEC2(IGDRY))
 !
-!        7.5.8  select the (ZLBDAG,ZLBDAR) couplet
+!        6.5.8  select the (ZLBDAG,ZLBDAR) couplet
 !
    ZVEC1(:) = PACK( ZLBDAG(:),MASK=GDRY(:) )
    ZVEC2(:) = PACK( ZLBDAR(:),MASK=GDRY(:) )
 !
-!        7.5.9  find the next lower indice for the ZLBDAG and for the ZLBDAR
+!        6.5.9  find the next lower indice for the ZLBDAG and for the ZLBDAR
 !               in the geometrical set of (Lbda_g,Lbda_r) couplet use to
 !               tabulate the RDRYG-kernel
 !
@@ -1044,7 +1014,7 @@ ENDIF
    IVEC2(1:IGDRY) = INT( ZVEC2(1:IGDRY) )
    ZVEC2(1:IGDRY) = ZVEC2(1:IGDRY) - FLOAT( IVEC2(1:IGDRY) )
 !
-!        7.5.10 perform the bilinear interpolation of the normalized
+!        6.5.10 perform the bilinear interpolation of the normalized
 !               RDRYG-kernel
 !
    DO JJ = 1,IGDRY
@@ -1075,7 +1045,7 @@ ENDIF
  ZRDRYG(:) = ZZW1(:,1) + ZZW1(:,2) + ZZW1(:,3) + ZZW1(:,4)
  DEALLOCATE(GDRY)
 !
-!*       7.6    compute the Wet growth case of the graupel
+!*       6.6    compute the Wet growth case of the graupel
 !
  ZZW(:) = 0.0
  ZRWETG(:) = 0.0
@@ -1104,7 +1074,7 @@ ENDIF
                                 ( ZRHODREF(:)*(XLMTT-XCL*(XTT-ZZT(:))) )   )
  END WHERE
 !
-!*       7.7    Select Wet or Dry case for the growth of the graupel
+!*       6.7    Select Wet or Dry case for the growth of the graupel
 !
  ZZW(:) = 0.0
  ZZW2(:,:) = 0.0
@@ -1121,26 +1091,28 @@ ENDIF
      ZZW3(:,:) = 0.0
      DO JLI = 1, SIZE(PSGRSVS,4)
        IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-          .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,INDEXWI(JLI))
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+          .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+          .OR. NINDEXGI(JLI).EQ.0) THEN
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * ZZW2(JL,NINDEXWI(JLI))
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETHP) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETHP) * ZZW2(JL,NINDEXWI(JLI))
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETSU) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETSU) * ZZW2(JL,NINDEXWI(JLI))
        ELSE
-         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,INDEXWI(JLI))
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) +                        &
-                                   (1. - XRETDF) * ZZW2(JL,INDEXWI(JLI))
+         ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF * ZZW2(JL,NINDEXWI(JLI))
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) +                        &
+                                   (1. - XRETDF) * ZZW2(JL,NINDEXWI(JLI))
        ENDIF
      ENDDO
      IF (ZRST(JL)>0.0) THEN
@@ -1148,14 +1120,14 @@ ENDIF
        ZZW3(JL,:) = MAX(MIN(ZZW3(JL,:),(ZSGSVT(JL,:)/PTSTEP)),0.0)
        ZSGRSVS(JL,:) = ZSGRSVS(JL,:) - ZZW3(JL,:) !snow->rain
        DO JLI = 1, SIZE(PSGRSVS,4)
-         ZRRSVS(JL,INDEXWI(JLI)) = ZRRSVS(JL,INDEXWI(JLI)) + ZZW3(JL,JLI)
+         ZRRSVS(JL,NINDEXWI(JLI)) = ZRRSVS(JL,NINDEXWI(JLI)) + ZZW3(JL,JLI)
        ENDDO
      ENDIF
      ELSE
        IF (.NOT.(OCH_RET_ICE)) THEN
          DO JLW = 1, SIZE(PRRSVS,4)
-           IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-             ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW)
+           IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+             ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW)
            ENDIF
          ENDDO
        ENDIF
@@ -1176,37 +1148,39 @@ ENDIF
      IF (OUSECHIC) THEN
      DO JLI = 1, SIZE(PSGRSVS,4)
        IF (TRIM(HICNAMES(JLI)) == 'IC_HNO3' .OR. TRIM(HICNAMES(JLI)) == 'IC_SULF' &
-          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. HICNAMES(JLI)(1:4) == 'IC_A' &
-          .OR. HICNAMES(JLI)(1:4) == 'IC_B' ) THEN
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_H2SO4' &
+          .OR. TRIM(HICNAMES(JLI)) == 'IC_NH3' .OR. TRIM(HICNAMES(JLI)) == 'IC_HCL' &
+          .OR. HICNAMES(JLI)(1:4) == 'IC_A' .OR. HICNAMES(JLI)(1:4) == 'IC_B' &
+          .OR. NINDEXGI(JLI).EQ.0) THEN
          ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETNA * (        &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_H2O2' .OR. TRIM(HICNAMES(JLI)) == 'IC_HO2' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HONO' .OR. TRIM(HICNAMES(JLI)) == 'IC_HNO4'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_HCHO' .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA1'&
           .OR. TRIM(HICNAMES(JLI)) == 'IC_ORA2') THEN
          ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETHP *  (       &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) + (1. - XRETHP) * (  &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) + (1. - XRETHP) * (  &
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
        ELSE IF (TRIM(HICNAMES(JLI)) == 'IC_SO2' .OR. TRIM(HICNAMES(JLI)) == 'IC_OH' &
           .OR. TRIM(HICNAMES(JLI)) == 'IC_MO2' .OR. &
                TRIM(HICNAMES(JLI)) == 'IC_OP1') THEN
          ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETSU *  (       &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) + (1. - XRETSU) * (  &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) + (1. - XRETSU) * (  &
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
        ELSE
          ZSGRSVS(JL,JLI) = ZSGRSVS(JL,JLI) + XRETDF *  (       &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
-         ZGRSVS(JL,INDEXGI(JLI)) = ZGRSVS(JL,INDEXGI(JLI)) + (1. - XRETDF) * (  &
-                           ZZW2(JL,INDEXWI(JLI)) +  ZZW4(JL,INDEXWI(JLI)) )
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
+         ZGRSVS(JL,NINDEXGI(JLI)) = ZGRSVS(JL,NINDEXGI(JLI)) + (1. - XRETDF) * (  &
+                           ZZW2(JL,NINDEXWI(JLI)) +  ZZW4(JL,NINDEXWI(JLI)) )
        ENDIF
      ENDDO
      ELSE
        IF (.NOT.(OCH_RET_ICE)) THEN
          DO JLW = 1, SIZE(PRRSVS,4)
-           IF (.NOT.(INDEXWG(JLW).EQ.0)) THEN
-             ZGRSVS(JL,INDEXWG(JLW)) = ZGRSVS(JL,INDEXWG(JLW)) + ZZW2(JL,JLW) &
+           IF (.NOT.(NINDEXWG(JLW).EQ.0)) THEN
+             ZGRSVS(JL,NINDEXWG(JLW)) = ZGRSVS(JL,NINDEXWG(JLW)) + ZZW2(JL,JLW) &
                                                               + ZZW4(JL,JLW)
            ENDIF
          ENDDO
@@ -1215,7 +1189,7 @@ ENDIF
    ENDIF
  ENDDO
 !
-!*       7.8    Melting of the graupel
+!*       6.8    Melting of the graupel
 !
  IF (OUSECHIC) THEN
  ZZW(:) = 0.0
@@ -1238,7 +1212,7 @@ ENDIF
      ZZW3(JL,:) = MAX(MIN(ZZW3(JL,:),(ZSGSVT(JL,:)/PTSTEP)),0.0)
      ZSGRSVS(JL,:) = ZSGRSVS(JL,:) - ZZW3(JL,:) !graupel->rain
      DO JLI = 1, SIZE(PSGRSVS,4)
-       ZRRSVS(JL,INDEXWI(JLI)) = ZRRSVS(JL,INDEXWI(JLI)) + ZZW3(JL,JLI)
+       ZRRSVS(JL,NINDEXWI(JLI)) = ZRRSVS(JL,NINDEXWI(JLI)) + ZZW3(JL,JLI)
      ENDDO
    ENDIF
  ENDDO
@@ -1247,7 +1221,7 @@ ENDIF
 !
 !-------------------------------------------------------------------------------
 !
-!*       8.     UNPACK RESULTS AND DEALLOCATE ARRAYS
+!*       7.     UNPACK RESULTS AND DEALLOCATE ARRAYS
 !               ------------------------------------
 
 
diff --git a/src/MNH/ch_f77.fx90 b/src/MNH/ch_f77.fx90
index e6094aa4f9adb0dc6dc5954ecc0c693091e53ca9..824465342ac8cd71327380963ff933e7924e45aa 100644
--- a/src/MNH/ch_f77.fx90
+++ b/src/MNH/ch_f77.fx90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$ $Date$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/ch_f77.fx90,v $ $Revision: 1.2.2.1.2.2.2.1.8.2.2.3 $ $Date: 2014/06/19 15:18:13 $
 !-----------------------------------------------------------------
 C**FILE:     svode.f
 C**AUTHOR:   Karsten Suhre
@@ -4788,6 +4788,7 @@ C
 *   alsurf = surface albedo, wavelength independent
 *   psurf = surface pressure, mbar.  Set to negative value to use
 *           US Standard Atmosphere, 1976 (USSA76)
+      psurf = -1.
 * Column amounts of absorbers (in Dobson Units, from surface to space):
 *          Vertical profile for O3 from USSA76.  For SO2 and NO2, vertical
 *          concentration profile is 2.69e10 molec cm-3 between 0 and 
@@ -8399,7 +8400,7 @@ C     srayl(iw) = 3.90e-28/(wmicrn)**xx
 
 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 c RCS version control information:
-c $Header$
+c $Header: /home/cvsroot/MNH-VX-Y-Z/src/MNH/ch_f77.fx90,v 1.2.2.1.2.2.2.1.8.2.2.3 2014/06/19 15:18:13 escj Exp $
 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
       SUBROUTINE SUNAE( YEAR, DAY, HOUR, LAT, LONG, lrefr,
diff --git a/src/MNH/ch_monitorn.f90 b/src/MNH/ch_monitorn.f90
index 65972ab681fc9192d835f5f42cb7287a059afb8b..5aa06bc42adbb419dab9cb020bc829343dd8cb6f 100644
--- a/src/MNH/ch_monitorn.f90
+++ b/src/MNH/ch_monitorn.f90
@@ -107,6 +107,9 @@ END MODULE MODI_CH_MONITOR_n
 !!                         imply transfer H2SO4 AP in aqueous phase if aq.chem.
 !!    04/2014 (C.Lac) Remove GCENTER with FIT temporal scheme
 !!    06/11/14 (M Leriche) Bug in pH computing
+!!    11/12/15 (M. Leriche & P. Tulet) add ch_init_ice initialise index for ice chem.
+!!    18/01/16 (M Leriche) for sedimentation fusion C2R2 and khko
+!!    15/02/16 (M Leriche) call ch_init_rosenbrock only one time
 !!
 !!    EXTERNAL
 !!    --------
@@ -117,13 +120,13 @@ USE MODI_CH_SET_PHOTO_RATES
 USE MODI_CH_SOLVER_n
 USE MODI_CH_UPDATE_JVALUES
 USE MODI_BUDGET
+USE MODI_CH_INIT_ICE
 USE MODI_CH_AQUEOUS_TMICICE
 USE MODI_CH_AQUEOUS_TMICKESS
 USE MODI_CH_AQUEOUS_TMICC2R2
 USE MODI_CH_AQUEOUS_TMICKHKO
 USE MODI_CH_AQUEOUS_SEDIM1MOM
-USE MODI_CH_AQUEOUS_SEDIMC2R2
-USE MODI_CH_AQUEOUS_SEDIMKHKO
+USE MODI_CH_AQUEOUS_SEDIM2MOM
 USE MODI_CH_AQUEOUS_CHECK
 USE MODI_FM_ll
 USE MODI_SUM_ll
@@ -147,7 +150,7 @@ USE MODD_CST, ONLY : XMNH_TINY
 USE MODD_BUDGET
 USE MODD_LUNIT_n
 USE MODD_NSV, ONLY : NSV_CHEMBEG,NSV_CHEMEND,NSV_CHEM,& ! index for chemical SV
-                     NSV_CHACBEG,NSV_CHACEND,         & ! index for aqueous SV
+                     NSV_CHACBEG,NSV_CHACEND,NSV_CHAC,& ! index for aqueous SV
                      NSV_CHGSBEG,NSV_CHGSEND,         & ! index for gas phase SV
                      NSV_CHICBEG,NSV_CHICEND,         & ! index for ice phase SV
                      NSV_C2R2BEG,                     & ! index for number concentration
@@ -195,7 +198,7 @@ USE MODD_FIELD_n,   ONLY: XSVT,      &! scalar variable at t
 USE MODD_REF_n,     ONLY: XRHODREF,  &! dry density for ref. state
                           XRHODJ      ! ( rhod J ) = dry density
 !
-USE MODD_TIME 
+USE MODD_TIME,      ONLY: TDTEXP 
 !
 USE MODD_TIME_n,    ONLY: TDTCUR      ! Current Time and Date
 !
@@ -228,6 +231,10 @@ USE MODD_DYN_n,     ONLY: XTSTEP      ! time step of MesoNH
 USE MODD_PRECIP_n, ONLY: XEVAP3D
 USE MODD_CLOUDPAR_n, ONLY: NSPLITR  ! Nb of required small time step integration
 !
+!variables used by microphysical mass transfer - sedimentation
+!
+USE MODD_CLOUDPAR_n, ONLY: NSPLITR
+!
 !variables used by rosenbrock solver
 !
 USE MODD_CH_ROSENBROCK_n, ONLY: NSPARSEDIM,   & ! Dim of NSPARSE_xxx vectors
@@ -329,7 +336,7 @@ INTEGER                :: IMI            ! model index
 !-------------------------------------------------------------------------------
 !   variables for the aerosol module
 !
-REAL                   :: ZTIME                ! time beginning at TDTEXP%TIME
+REAL                   :: ZTIME                ! current time 
 REAL, ALLOCATABLE, DIMENSION(:,:)   :: ZM, ZSIG0, ZN0, ZRG0, &   ! work array
                                        ZCTOTG, ZSEDA, ZFRAC, ZMI ! for aerosols
 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZCTOTA, ZCCTOT
@@ -348,6 +355,7 @@ REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NSV_AER) :: ZCWETAERO
 INTEGER                :: JRR          ! Loop index for the moist variables
 REAL,DIMENSION(SIZE(XRT,1),SIZE(XRT,2),SIZE(XRT,3),SIZE(XRT,4))     :: ZRT_VOL
                                        ! liquid content in vol/vol
+REAL, DIMENSION(SIZE(XRT,1), SIZE(XRT,2))     :: ZINPRR! Rain instant precip
 !
 !-------------------------------------------------------------------------------
 !
@@ -529,6 +537,11 @@ IF (KTCOUNT == 1) THEN
   ALLOCATE(LU_DIM_SPECIES(ISVECNPT))
   LU_DIM_SPECIES(:) = NEQ
 !
+!        1.1.3 determine index for ice phase chemistry or degassing with ICE3/4
+  IF ((LUSECHAQ).AND.((CCLOUD=='ICE3' .OR. CCLOUD=='ICE4'))) THEN 
+     CALL CH_INIT_ICE(LUSECHIC,LCH_RET_ICE,CNAMES,CICNAMES,NEQ,NEQAQ) 
+  ENDIF
+!
 ENDIF  ! first time step
 !
 !*       1.2   calculate timestep variables
@@ -614,13 +627,10 @@ END IF
 !*       2.    UPDATE PHOTOLYSIS RATES
 !              -----------------------
 !
-ZTIME  = TDTCUR%TIME
-!
 IF (KTCOUNT==1 .OR. &
     (MOD(ISTCOUNT, MAX(1, INT(XCH_TUV_TUPDATE/XTSTEP)) ) .EQ. 0)) THEN
 !
-  WRITE(KLUOUT,*)"TIME: ", (TDTCUR%TIME - TDTEXP%TIME)                  &
-                     + 86400.*(TDTCUR%TDATE%DAY - TDTEXP%TDATE%DAY)
+  WRITE(KLUOUT,*)"TIME call update jvalue: ",TDTCUR%TIME
 !
   IF (.NOT.ASSOCIATED(XJVALUES)) &
              ALLOCATE(XJVALUES(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),JPJVMAX))
@@ -644,6 +654,7 @@ ISTCOUNT = ISTCOUNT + 1
 !*       3.1 sedimentation term and wet deposition for aerosols tendency (XSEDA)
 !
 IF (LORILAM) THEN
+  ZTIME  = TDTCUR%TIME ! need for ch_orilam
   XSEDA(:,:,:,:) = 0.
   ZSEDA(:,:) = 0.
 ! dry sedimentation
@@ -748,29 +759,23 @@ IF (LUSECHAQ.AND.(NRRL>=2) ) THEN
   IF (MAXVAL(ZRT_VOL(:,:,:,3))>XRTMIN_AQ) THEN
     SELECT CASE ( CCLOUD )
       CASE ('KESS','ICE3','ICE4')
-        CALL CH_AQUEOUS_SEDIM1MOM(TDTCUR%TIME, CCLOUD, LUSECHIC,                &
+        CALL CH_AQUEOUS_SEDIM1MOM(NSPLITR, CCLOUD, LUSECHIC,                    &
                                   PTSTEP , XZZ, XRHODREF,                       &
                                   XRHODJ, XRRS(:,:,:,3), XRRS(:,:,:,5),         &
                                   XRRS(:,:,:,6),                                &
                                   XRSVS(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND), &
-                                  XRSVS(:,:,:,NSV_CHICBEG:NSV_CHICEND)          )
+                                  XRSVS(:,:,:,NSV_CHICBEG:NSV_CHICEND),         &
+                                  ZINPRR(:,:)                                   )
 
-      CASE ('C2R2','C3R5')
-        CALL CH_AQUEOUS_SEDIMC2R2(TDTCUR%TIME, PTSTEP, XRTMIN_AQ,              &
+      CASE ('C2R2','C3R5','KHKO')
+        CALL CH_AQUEOUS_SEDIM2MOM(NSPLITR, CCLOUD, PTSTEP, XRTMIN_AQ,         &
                                   XZZ, XRHODREF, XRHODJ,                       &
                                   XRT(:,:,:,3),XRRS(:,:,:,3),                  &
                                   XSVT(:,:,:,NSV_C2R2BEG+2),                   &
                                   XRSVS(:,:,:,NSV_C2R2BEG+2),                  &
                                   XSVT(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND), &
-                                  XRSVS(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND) )
-
-      CASE ('KHKO')
-        CALL CH_AQUEOUS_SEDIMKHKO(PTSTEP , XZZ, XRHODREF, XRHODJ,                &
-                                  XRT(:,:,:,3), XRRS(:,:,:,3),                   &
-                                  XSVT(:,:,:,NSV_C2R2BEG+2),                     &
-                                  XRSVS(:,:,:,NSV_C2R2BEG+2),                    &
-                                  XSVT(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND),   &
-                                  XRSVS(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND)   )
+                                  XRSVS(:,:,:,NSV_CHACBEG+NEQAQ/2:NSV_CHACEND),&
+                                  ZINPRR(:,:)                                  )
     END SELECT
   END IF
 ELSE IF (LUSECHAQ.AND.(NRRL==1) ) THEN
@@ -1104,6 +1109,27 @@ DO JL=1,ISVECNMASK
   ENDIF
 END DO
 !
+!*        4.9  compute accumalated concentrations in rain at the surface
+!
+IF (CCLOUD /= 'REVE' ) THEN
+  IF (LUSECHAQ) THEN
+    DO JSV=1,NSV_CHAC/2
+      WHERE((XRRS(:,:,IKB,3) .GT. 0.).AND.(XRSVS(:,:,IKB,JSV+NSV_CHACBEG+NSV_CHAC/2-1).GT.0.))
+          XACPRAQ(:,:,JSV) = XACPRAQ(:,:,JSV) + &
+              (XRSVS(:,:,IKB,JSV+NSV_CHACBEG+NSV_CHAC/2-1))/ (XMD*XRRS(:,:,IKB,3))*& ! moles i  / kg eau
+               1E3*ZINPRR(:,:) * XTSTEP ! moles i / m2
+      END WHERE
+    ENDDO
+    IF (LCH_PH) THEN
+      WHERE ((ZINPRR(:,:)>0.).AND.(XPHR(:,:,IKB)>0.))
+      ! moles of H+ / m2
+        XACPHR(:,:) =  XACPHR(:,:) + 1E3*ZINPRR(:,:) * XTSTEP * &
+                     10**(-XPHR(:,:,IKB))
+      END WHERE
+    END IF
+  END IF
+END IF
+
 !
 IF (LBUDGET_SV) THEN
   DO JSV=NSV_CHEMBEG,NSV_CHEMEND
@@ -1235,17 +1261,7 @@ END DO
 ! system dimensions.
 !
     IF (KTCOUNT == 1) THEN
-      IF( JL>1 ) THEN
-        DEALLOCATE(NSPARSE_IROW)
-        DEALLOCATE(NSPARSE_ICOL)
-        DEALLOCATE(NSPARSE_CROW)
-        DEALLOCATE(NSPARSE_DIAG)
-        DEALLOCATE(NSPARSE_IROW_NAQ)
-        DEALLOCATE(NSPARSE_ICOL_NAQ)
-        DEALLOCATE(NSPARSE_CROW_NAQ)
-        DEALLOCATE(NSPARSE_DIAG_NAQ)
-      END IF
-      CALL CH_INIT_ROSENBROCK(IMI,KLUOUT)
+      IF (JL==1) CALL CH_INIT_ROSENBROCK(IMI,KLUOUT)
       IF( ASSOCIATED(LU_DIM_SPECIES) ) THEN
         DEALLOCATE(LU_DIM_SPECIES)
       END IF
diff --git a/src/MNH/ch_orilam.f90 b/src/MNH/ch_orilam.f90
index 723cf9aef74721bf515a6920664f2bfa147f95bb..3e4e1219db5865ec837ded090bd7eb733a73ceaa 100644
--- a/src/MNH/ch_orilam.f90
+++ b/src/MNH/ch_orilam.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/ch_orilam.f90,v $ $Revision: 1.1.2.1.2.1.18.1 $
 ! MASDEV4_7 chimie 2007/03/02 13:59:37
 !-----------------------------------------------------------------
 !!   #########################
@@ -80,6 +80,7 @@ SUBROUTINE CH_ORILAM(PAERO, PCHEM, PM, PSIG0, PRG0, PN0, PCTOTG, PCTOTA,&
 !!    MODIFICATIONS
 !!    -------------
 !!    Original
+!!    M. Leriche (08/16) add initialization of ZMASK
 !!
 !!    EXTERNAL
 !!    --------
@@ -117,6 +118,9 @@ CHARACTER(LEN=10),                      INTENT(IN)    :: GSCHEME
 REAL, DIMENSION(SIZE(PAERO,1),JPMODE)                 :: ZMASK
 !
 !-------------------------------------------------------------------------------
+!initialize ZMASK
+ZMASK(:,:) = 1.
+!
 ! transfer gas phase variables into aerosol variables
 CALL CH_AER_TRANS(0, PM, PSIG0, PRG0, PN0, PRHOP0,PAERO, PCHEM, PCTOTG, PCTOTA, PCCTOT,&
                          PFRAC, PMI, ZMASK,GSCHEME)
diff --git a/src/MNH/ch_solve_ph.f90 b/src/MNH/ch_solve_ph.f90
index 0c0c430f7d6555f06215d73d8bf9ea3397db9af4..f9b994cdf4c09e497aaa3dc9d9e887cb9bbd09c0 100644
--- a/src/MNH/ch_solve_ph.f90
+++ b/src/MNH/ch_solve_ph.f90
@@ -131,7 +131,7 @@ C0 = 0.                             !NH3
 C1 = 0.                             !CO2
 C2 = 0.                             !SO2
 C3 = 0.                             !HCOOH = ORA1
-C4 = 0.                             !HNO3 + 2 x H2SO4 = strong acid
+C4 = 0.                             !HNO3 + 2 x H2SO4 + HCL = strong acid
 SOM = 0.
 IORDER = 8 !polynomial order
 ALLOCATE(ZCOEFS(KLW,IORDER+1))
@@ -153,6 +153,7 @@ SELECT CASE (KRR)
       IF (TRIM(CNAMES(JJ))=='WC_HNO3') C4(:)= C4(:)+PCONC(:,JI)/(ZFACT(:))
       IF ((TRIM(CNAMES(JJ))=='WC_SULF') .OR. (TRIM(CNAMES(JJ))=='WC_H2SO4')) &
           C4(:)= C4(:)+2.*PCONC(:,JI)/(ZFACT(:))
+      IF (TRIM(CNAMES(JJ))=='WC_HCL') C4(:)= C4(:)+PCONC(:,JI)/(ZFACT(:))
       IF (CNAMES(JJ)(1:4)=='WC_A') SOM(:) = SOM(:) + PCONC(:,JI)/(ZFACT(:))
       IF (CNAMES(JJ)(1:4)=='WC_B') SOM(:) = SOM(:) + 2.*PCONC(:,JI)/(ZFACT(:))
     END DO
@@ -166,6 +167,7 @@ SELECT CASE (KRR)
       IF (TRIM(CNAMES(JJ))=='WR_HNO3') C4(:)= C4(:)+PCONC(:,JI)/(ZFACT(:))
       IF ((TRIM(CNAMES(JJ))=='WR_SULF') .OR. (TRIM(CNAMES(JJ))=='WR_H2SO4')) &
           C4(:)= C4(:)+2.*PCONC(:,JI)/(ZFACT(:))
+      IF (TRIM(CNAMES(JJ))=='WR_HCL') C4(:)= C4(:)+PCONC(:,JI)/(ZFACT(:))
       IF (CNAMES(JJ)(1:4)=='WR_A') SOM(:) = SOM(:) + PCONC(:,JI)/(ZFACT(:))
       IF (CNAMES(JJ)(1:4)=='WR_B') SOM(:) = SOM(:) + 2.*PCONC(:,JI)/(ZFACT(:))
     END DO
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index c99ba20e770585e62f48dfc7aeb59936e515c970..f1640aa5e49ebe2be38be93491f742f2563a2c82 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -208,6 +208,9 @@ END MODULE MODI_DEFAULT_DESFM_n
 !!                   07/2013  (C.Lac) add WENO, LCHECK              
 !!                   07/2013  (Bosseur & Filippi) adds Forefire
 !!                   08/2015  (Redelsperger & Pianezze) add XPOND coefficient for LBC
+!!      Modification 24/03/16 (Leriche) remove LCH_SURFACE_FLUX 
+!!                                      put NCH_VEC_LENGTH = 50 instead of 1000
+!!
 !!                   04/2016 (C.LAC) negative contribution to the budget splitted between advection, turbulence and microphysics for KHKO/C2R2
 !-------------------------------------------------------------------------------
 !
@@ -1059,7 +1062,6 @@ LUSECHEM            = .FALSE.
 LUSECHAQ            = .FALSE.
 LUSECHIC            = .FALSE.
 LCH_INIT_FIELD      = .FALSE.
-LCH_SURFACE_FLUX    = .FALSE.
 LCH_CONV_SCAV       = .FALSE.
 LCH_CONV_LINOX      = .FALSE.
 LCH_PH              = .FALSE.
@@ -1076,7 +1078,7 @@ XCH_TUV_ALBNEW      = -1.
 XCH_TUV_DOBNEW      = -1.
 XCH_TUV_TUPDATE     = 600.
 CCH_VEC_METHOD      = 'MAX'
-NCH_VEC_LENGTH      = 1000
+NCH_VEC_LENGTH      = 50
 XCH_TS1D_TSTEP      = 600.
 CCH_TS1D_COMMENT    = 'no comment'
 CCH_TS1D_FILENAME   = 'IO1D'
diff --git a/src/MNH/exchange.f90 b/src/MNH/exchange.f90
index ac14a5b9b9edfe1c603f73528a7987ff161b2728..b96fb5154aad1a534d056a61752d31b3cfee1e8a 100644
--- a/src/MNH/exchange.f90
+++ b/src/MNH/exchange.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$ $Date$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/exchange.f90,v $ $Revision: 1.2.2.2.2.2.16.1.2.5.2.1 $ $Date: 2015/12/01 15:26:23 $
 !-----------------------------------------------------------------
 !-----------------------------------------------------------------
 !-----------------------------------------------------------------
@@ -82,13 +82,18 @@ END MODULE MODI_EXCHANGE
 !!    -------------
 !!
 !!    original     18/09/98
-!!                 05/2006   Remove KEPS
-!!                 10/2009 (C.Lac) FIT for variables advected by PPM
-!!                 05/2014 (C.Lac) Correction of negative values of chemical
+!!      05/2006   Remove KEPS
+!!      10/2009 (C.Lac) FIT for variables advected by PPM
+!!      05/2014 (C.Lac) Correction of negative values of chemical
 !!                   tracers moved from ch_monitor to the end of the time step
-!!                 11/2014 (G.Delautier) Call correction of negative values only
-!!                         if LUSECHEM 
-!------------------------------------------------------------------------------
+!!      11/2014 (G.Delautier) Call correction of negative values only if LUSECHEM 
+!!      16/02/16 (M. Leriche) conserve total mass for gas phase chem.
+!!                   species only, remove negative values without mass 
+!!                   corrections for aq. phase and ice phase (lost mass neglig.)
+!!      25/08/16 (M.Leriche) remove negative values for aerosols and conserve
+!!                   total mass for chemical species in aerosols
+!!            
+!-----------------------------------------------------------------------------------------
 !
 !*      0.   DECLARATIONS
 !            ------------
@@ -104,7 +109,8 @@ USE MODD_LUNIT_n,     ONLY : CLUOUT
 USE MODI_SHUMAN
 USE MODI_SUM_ll
 USE MODI_BUDGET
-USE MODD_CH_MNHC_n, ONLY : LUSECHEM
+USE MODD_CH_MNHC_n,   ONLY : LUSECHEM, LUSECHAQ, LUSECHIC
+USE MODD_CH_AEROSOL,  ONLY : LORILAM, NM6_AER
 !
 IMPLICIT NONE
 !
@@ -158,10 +164,10 @@ IF (SIZE(PRTKES,1) /= 0) PRTKES(:,:,:) = PRTKES(:,:,:)*PTSTEP/PRHODJ
 !      REMOVE NEGATIVE VALUES OF CHEM SCALAR
 !
 IF (LUSECHEM) THEN
-  DO JSV=NSV_CHEMBEG,NSV_CHEMEND
+  DO JSV=NSV_CHGSBEG,NSV_CHGSEND
     IF ( MIN_ll( PRSVS(:,:,:,JSV), IINFO_ll) < 0.0 ) THEN
 !
-! compute the total water mass computation
+! compute the total mass 
 !
       ZMASSTOT = MAX( 0. , SUM3D_ll( PRSVS(:,:,:,JSV), IINFO_ll ) )
 !
@@ -178,16 +184,79 @@ IF (LUSECHEM) THEN
       ZRATIO = ZMASSTOT / ZMASSPOS
       PRSVS(:,:,:,JSV) = PRSVS(:,:,:,JSV) * ZRATIO
 !
-      WRITE(ILUOUT,*)'DUE TO CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
-      WRITE(ILUOUT,*)'SOURCES IS CORRECTED BY RATIO',ZRATIO
+      WRITE(ILUOUT,*)'CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
+      WRITE(ILUOUT,*)'GAS SOURCES IS CORRECTED BY RATIO',ZRATIO
     END IF
   END DO
+  IF (LUSECHAQ) THEN
+    DO JSV =  NSV_CHACBEG, NSV_CHACEND
+      IF ( MIN_ll( PRSVS(:,:,:,JSV), IINFO_ll) < 0.0 ) THEN
+! remove the negative values
+        PRSVS(:,:,:,JSV) = MAX(0., PRSVS(:,:,:,JSV) )
+        WRITE(ILUOUT,*)'CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
+        WRITE(ILUOUT,*)'CLOUD OR RAIN SOURCE IS SET TO ZERO'
+      END IF
+    END DO
+  ENDIF
+!
+  IF (LUSECHIC) THEN
+    DO JSV =  NSV_CHICBEG, NSV_CHICEND
+      IF ( MIN_ll( PRSVS(:,:,:,JSV), IINFO_ll) < 0.0 ) THEN
+! remove the negative values
+        PRSVS(:,:,:,JSV) = MAX(0., PRSVS(:,:,:,JSV) )
+        WRITE(ILUOUT,*)'CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
+        WRITE(ILUOUT,*)'ICE PHASE SOURCE IS SET TO ZERO'
+      END IF
+    END DO
+  ENDIF
 !
   IF (LBUDGET_SV) THEN
     DO JSV=NSV_CHEMBEG,NSV_CHEMEND
       CALL BUDGET(PRSVS(:,:,:,JSV),JSV+12,'NEGA_BU_RSV')
     ENDDO
   ENDIF
+!
+! aerosol sv
+  IF (LORILAM) THEN
+    DO JSV=NSV_AERBEG,NSV_AEREND-2-NM6_AER ! keep chem. species only
+      IF ( MIN_ll( PRSVS(:,:,:,JSV), IINFO_ll) < 0.0 ) THEN
+!
+! compute the total mass
+!
+        ZMASSTOT = MAX( 0. , SUM3D_ll( PRSVS(:,:,:,JSV), IINFO_ll ) )
+!
+! remove the negative values
+!
+        PRSVS(:,:,:,JSV) = MAX(0., PRSVS(:,:,:,JSV) )
+!
+! compute the new total mass
+!
+        ZMASSPOS = MAX(XMNH_TINY,SUM3D_ll( PRSVS(:,:,:,JSV), IINFO_ll ) )
+!
+! correct again in such a way to conserve the total mass 
+!
+        ZRATIO = ZMASSTOT / ZMASSPOS
+        PRSVS(:,:,:,JSV) = PRSVS(:,:,:,JSV) * ZRATIO
+!
+        WRITE(ILUOUT,*)'CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
+        WRITE(ILUOUT,*)'AP SOURCES IS CORRECTED BY RATIO',ZRATIO
+      END IF
+    END DO
+!
+    DO JSV=NSV_AEREND-2-NM6_AER+1,NSV_AEREND
+      IF ( MIN_ll( PRSVS(:,:,:,JSV), IINFO_ll) < 0.0 ) THEN
+! remove the negative values for M0 and M6
+         PRSVS(:,:,:,JSV) = MAX(0., PRSVS(:,:,:,JSV) )
+         WRITE(ILUOUT,*)'CHEMISTRY',JSV,'HAS NEGATIVE VALUES'
+         WRITE(ILUOUT,*)'AP MOMENT SOURCES IS SET TO ZERO'
+      END IF
+    END DO
+    IF (LBUDGET_SV) THEN
+      DO JSV=NSV_AERBEG,NSV_AEREND
+        CALL BUDGET(PRSVS(:,:,:,JSV),JSV+12,'NEGA_BU_RSV')
+      ENDDO
+    ENDIF
+  ENDIF
 ENDIF
 !
 DO JSV=1,KSV
diff --git a/src/MNH/goto_model_wrapper.f90 b/src/MNH/goto_model_wrapper.f90
index 24127f9a2d4e5713dadf2280a481e659ce050d4a..da89c864a1e2ef7e7ea0c1632bdbe52a0653e2d0 100644
--- a/src/MNH/goto_model_wrapper.f90
+++ b/src/MNH/goto_model_wrapper.f90
@@ -13,6 +13,8 @@
 !!      06/12 (Tomasini) Grid-nesting of ADVFRC and EDDY_FLUX
 !!      07/13 (Bosseur & Filippi) adds Forefire
 !!      2014 (Faivre)
+!!      2016  (Leriche) Add MODD_CH_ICE Suppress MODD_CH_DEP_n
+!!      Modification    01/2016  (JP Pinty) Add LIMA
 !-----------------------------------------------------------------
 MODULE MODI_GOTO_MODEL_WRAPPER
 
@@ -29,7 +31,6 @@ SUBROUTINE GOTO_MODEL_WRAPPER(KFROM, KTO)
 USE MODD_ADV_n
 USE MODD_BIKHARDT_n
 USE MODD_CH_AERO_n
-USE MODD_CH_DEP_n
 USE MODD_CH_JVALUES_n
 USE MODD_CH_MNHC_n
 USE MODD_CH_SOLVER_n
@@ -95,6 +96,7 @@ USE MODD_TIMEZ
 USE MODD_SUB_PASPOL_n
 USE MODD_SUB_ELEC_n
 USE MODD_CH_PH_n
+USE MODD_CH_ICE_n
 USE MODD_CH_M9_n
 USE MODD_CH_ROSENBROCK_n
 USE MODD_RBK90_Global_n
@@ -113,7 +115,6 @@ INTEGER,INTENT(IN) :: KFROM, KTO
 CALL ADV_GOTO_MODEL(KFROM, KTO)
 CALL BIKHARDT_GOTO_MODEL(KFROM, KTO)
 CALL CH_AERO_GOTO_MODEL(KFROM,KTO)
-CALL CH_DEP_GOTO_MODEL(KFROM, KTO)
 CALL CH_JVALUES_GOTO_MODEL(KFROM, KTO)
 CALL CH_MNHC_GOTO_MODEL(KFROM, KTO)
 CALL CH_SOLVER_GOTO_MODEL(KFROM, KTO)
@@ -178,6 +179,7 @@ CALL TIME_GOTO_MODEL(KFROM, KTO)
 CALL TURB_GOTO_MODEL(KFROM, KTO)
 CALL TIMEZ_GOTO_MODEL(KFROM, KTO)
 CALL CH_PH_GOTO_MODEL(KFROM, KTO)
+CALL CH_ICE_GOTO_MODEL(KFROM, KTO)
 CALL CH_M9_GOTO_MODEL(KFROM, KTO)
 CALL CH_ROSENBROCK_GOTO_MODEL(KFROM, KTO)
 CALL RBK90_Global_GOTO_MODEL(KFROM, KTO)
diff --git a/src/MNH/ground_paramn.f90 b/src/MNH/ground_paramn.f90
index ef0cccfbc9a72388aa40c92c607c942069469af6..08f4679261df37ab1b9218d03f8ee33c5c9dd1a0 100644
--- a/src/MNH/ground_paramn.f90
+++ b/src/MNH/ground_paramn.f90
@@ -104,6 +104,7 @@ END MODULE MODI_GROUND_PARAM_n
 !!     (D.Gazen)              01/12/03  change emissions handling for surf. externalization
 !!     (J.escobar)            18/10/2012 missing USE MODI_COUPLING_SURF_ATM_n & MODI_DIAG_SURF_ATM_n
 !      (J.escobar)            2/2014 add Forefire coupling
+!!      (M.Leriche)            24/03/16 remove flag for chemical surface fluxes
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -112,7 +113,7 @@ END MODULE MODI_GROUND_PARAM_n
 USE MODD_CST,        ONLY : XP00, XCPD, XRD, XRV,XRHOLW, XDAY, XPI, XLVTT, XMD, XAVOGADRO
 USE MODD_PARAMETERS, ONLY : JPVEXT, XUNDEF
 USE MODD_DYN_n,      ONLY : XTSTEP
-USE MODD_CH_MNHC_n,  ONLY : LCH_SURFACE_FLUX
+USE MODD_CH_MNHC_n,  ONLY : LUSECHEM
 USE MODD_FIELD_n,    ONLY : XUT, XVT, XWT, XTHT, XRT, XPABST, XSVT, XTKET
 USE MODD_METRICS_n,  ONLY : XDXX, XDYY, XDZZ
 USE MODD_DIM_n,      ONLY : NKMAX
@@ -133,11 +134,11 @@ USE MODD_GRID,       ONLY : XLON0, XRPK, XBETA
 USE MODD_PARAM_ICE,  ONLY : LSEDIC
 USE MODD_PARAM_C2R2, ONLY : LSEDC
 USE MODD_DIAG_IN_RUN
-USE MODD_DUST,       ONLY : LDUST, CDUSTNAMES
-USE MODD_SALT,       ONLY : LSALT, CSALTNAMES
-USE MODD_CH_AEROSOL
-USE MODD_CSTS_DUST
-USE MODD_CSTS_SALT
+USE MODD_DUST,       ONLY : LDUST 
+USE MODD_SALT,       ONLY : LSALT
+USE MODD_CH_AEROSOL, ONLY : LORILAM
+USE MODD_CSTS_DUST,  ONLY : XMOLARWEIGHT_DUST
+USE MODD_CSTS_SALT,  ONLY : XMOLARWEIGHT_SALT
 !
 USE MODI_NORMAL_INTERPOL
 USE MODI_ROTATE_WIND
@@ -581,7 +582,7 @@ END IF
 !
 !* conversion from chemistry flux (molec/m2/s) to (ppp.m.s-1)
 !
-IF (LCH_SURFACE_FLUX) THEN
+IF (LUSECHEM) THEN
   DO JSV=NSV_CHEMBEG,NSV_CHEMEND
     PSFSV(:,:,JSV) = ZSFTS(:,:,JSV) * XMD / ( XAVOGADRO * XRHODREF(:,:,IKB)) 
   END DO
diff --git a/src/MNH/ini_lb.f90 b/src/MNH/ini_lb.f90
index 86a1a71b81d4942d7e6ca75308ff6efe428a75fe..2022a8a1dd3d81d90555d73690a6203666fd940e 100644
--- a/src/MNH/ini_lb.f90
+++ b/src/MNH/ini_lb.f90
@@ -126,6 +126,7 @@ SUBROUTINE INI_LB(HINIFILE,HLUOUT,OLSOURCE,KSV,                    &
 !!      M.Leriche       16/07/10    Add ice phase chemical species
 !!      Pialat/tulet    15/02/12    Add ForeFire scalars 
 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
+!!      M.Leriche       09/02/16    Treat gas and aq. chemicals separately
 !!      J.Escobar : 27/04/2016 : bug , test only on ANY(HGETSVM({{1:KSV}})=='READ'
 !-------------------------------------------------------------------------------
 !
@@ -709,12 +710,12 @@ DO JSV = NSV_ELECBEG, NSV_ELECEND
     IF ( SIZE(PLBYSVM,1) /= 0 ) PLBYSVM(:,:,:,JSV) = 0.
   END SELECT
 END DO
-! Chemical scalar variables
-DO JSV = NSV_CHEMBEG, NSV_CHEMEND
+! Chemical gas phase scalar variables
+DO JSV = NSV_CHGSBEG, NSV_CHGSEND
   SELECT CASE(HGETSVM(JSV))
   CASE ('READ')
     IF ( KSIZELBXSV_ll /= 0 ) THEN
-      YRECFM = 'LBX_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHEMBEG+1)))
+      YRECFM = 'LBX_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHGSBEG+1)))
       YDIRLB='LBX'
       CALL FMREAD_LB(HINIFILE,YRECFM,HLUOUT,YDIRLB,PLBXSVM(:,:,:,JSV),IRIMX,IL3DX,&
            & IGRID,ILENCH,YCOMMENT,IRESP)
@@ -724,7 +725,7 @@ DO JSV = NSV_CHEMBEG, NSV_CHEMEND
             PLBXSVM(:,:,:,JSV)=PLBXSVMM(:,:,:,JSV)
             WRITE(ILUOUT,*) 'Chemical PLBXSVM   will be initialized to 0'
           ELSE
-            WRITE(ILUOUT,*) 'Pb to initialize Chemical PLBXSVM '
+            WRITE(ILUOUT,*) 'Pb to initialize gas phase Chemical PLBXSVM '
 !callabortstop
             CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
             CALL ABORT
@@ -735,7 +736,7 @@ DO JSV = NSV_CHEMBEG, NSV_CHEMEND
     END IF
 !
     IF (KSIZELBYSV_ll  /= 0 ) THEN
-      YRECFM = 'LBY_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHEMBEG+1)))
+      YRECFM = 'LBY_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHGSBEG+1)))
       YDIRLB='LBY'
       CALL FMREAD_LB(HINIFILE,YRECFM,HLUOUT,YDIRLB,PLBYSVM(:,:,:,JSV),IRIMY,IL3DY,&
            & IGRID,ILENCH,YCOMMENT,IRESP)
@@ -745,7 +746,57 @@ DO JSV = NSV_CHEMBEG, NSV_CHEMEND
             PLBYSVM(:,:,:,JSV)=PLBYSVMM(:,:,:,JSV)
             WRITE(ILUOUT,*) 'Chemical PLBYSVM   will be initialized to 0'
           ELSE
-            WRITE(ILUOUT,*) 'Pb to initialize Chemical PLBYSVM '
+            WRITE(ILUOUT,*) 'Pb to initialize gas phase Chemical PLBYSVM '
+!callabortstop
+            CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
+            CALL ABORT
+            STOP
+          ENDIF
+        END IF
+      END IF
+    END IF
+  CASE('INIT')
+    IF ( SIZE(PLBXSVM,1) /= 0 ) PLBXSVM(:,:,:,JSV) = 0.
+    IF ( SIZE(PLBYSVM,1) /= 0 ) PLBYSVM(:,:,:,JSV) = 0.
+  END SELECT
+END DO
+! Chemical aqueous phase scalar variables
+DO JSV = NSV_CHACBEG, NSV_CHACEND
+  SELECT CASE(HGETSVM(JSV))
+  CASE ('READ')
+    IF ( KSIZELBXSV_ll /= 0 ) THEN
+      YRECFM = 'LBX_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHACBEG+NSV_CHGS+1)))
+      YDIRLB='LBX'
+      CALL FMREAD_LB(HINIFILE,YRECFM,HLUOUT,YDIRLB,PLBXSVM(:,:,:,JSV),IRIMX,IL3DX,&
+           & IGRID,ILENCH,YCOMMENT,IRESP)
+      IF ( SIZE(PLBXSVM,1) /= 0 ) THEN
+        IF (IRESP/=0) THEN
+          IF (PRESENT(PLBXSVMM)) THEN
+            PLBXSVM(:,:,:,JSV)=PLBXSVMM(:,:,:,JSV)
+            WRITE(ILUOUT,*) 'Chemical PLBXSVM   will be initialized to 0'
+          ELSE
+            WRITE(ILUOUT,*) 'Pb to initialize aqueous phase chemical PLBXSVM '
+!callabortstop
+            CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
+            CALL ABORT
+            STOP
+          ENDIF
+        END IF
+      END IF
+    END IF
+!
+    IF (KSIZELBYSV_ll  /= 0 ) THEN
+      YRECFM = 'LBY_'//TRIM(UPCASE(CNAMES(JSV-NSV_CHACBEG+NSV_CHGS+1)))
+      YDIRLB='LBY'
+      CALL FMREAD_LB(HINIFILE,YRECFM,HLUOUT,YDIRLB,PLBYSVM(:,:,:,JSV),IRIMY,IL3DY,&
+           & IGRID,ILENCH,YCOMMENT,IRESP)
+      IF ( SIZE(PLBYSVM,1) /= 0 ) THEN
+        IF (IRESP/=0) THEN
+          IF (PRESENT(PLBYSVMM)) THEN
+            PLBYSVM(:,:,:,JSV)=PLBYSVMM(:,:,:,JSV)
+            WRITE(ILUOUT,*) 'Chemical PLBYSVM   will be initialized to 0'
+          ELSE
+            WRITE(ILUOUT,*) 'Pb to initialize aqueous phase chemical PLBYSVM '
 !callabortstop
             CALL CLOSE_ll(HLUOUT,IOSTAT=IRESP)
             CALL ABORT
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 7b4983029728b9cacb107e4731a83a1895c1a8c9..ee0eedb62c97c0321140fde8c67fde61d1798a5d 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1,4 +1,3 @@
-
 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
@@ -266,6 +265,7 @@ END MODULE MODI_INI_MODEL_n
 !!       V. Masson     Feb 2015 replaces, for aerosols, cover fractions by sea, town, bare soil fractions
 !!                   J.Escobar : 19/04/2016 : Pb IOZ/NETCDF , missing OPARALLELIO=.FALSE. for PGD files
 !!                   J.Escobar : 01/06/2016 : correct check limit of NRIM versus local subdomain size IDIM
+!!                   M.Leriche 2016 Chemistry
 !---------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -1449,10 +1449,14 @@ IF (LUSECHAQ.AND.(CPROGRAM == 'DIAG  '.OR.CPROGRAM == 'MESONH')) THEN
     ALLOCATE(XPHC(IIU,IJU,IKU))
     IF (NRRL==2) THEN
       ALLOCATE(XPHR(IIU,IJU,IKU))
+      ALLOCATE(XACPHR(IIU,IJU))
+      XACPHR(:,:) =  0.
     ENDIF
   ENDIF
-  ALLOCATE(XACPRAQ(IIU,IJU,NSV_CHAC/2))
-  XACPRAQ(:,:,:) = 0.
+  IF (NRRL==2) THEN
+    ALLOCATE(XACPRAQ(IIU,IJU,NSV_CHAC/2))
+    XACPRAQ(:,:,:) = 0.
+  ENDIF
 ENDIF
 !
 !-------------------------------------------------------------------------------
diff --git a/src/MNH/modd_ch_aerosol.f90 b/src/MNH/modd_ch_aerosol.f90
index c6ff57f2edeff8de23f595ffae0f3c987e2724a1..e0019a581fbb165a4ee631c4bca6c7918c18cd96 100644
--- a/src/MNH/modd_ch_aerosol.f90
+++ b/src/MNH/modd_ch_aerosol.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_ch_aerosol.f90,v $ $Revision: 1.1.2.2.2.1.2.1.2.1.2.2 $
 ! MASDEV4_7 modd 2007/03/02 13:59:38
 !-----------------------------------------------------------------
 !!     ######################
@@ -34,6 +34,7 @@
 !!     MODIFICATIONS
 !!     -------------
 !!     (30-01-01) P.Tulet (LA) * modifications for secondary biogenics aerosols
+!!     (25-08-16) M.Leriche (LA) * NM6_AER is now in SAVE and assign in ini_nsv
 !!
 !!--------------------------------------------------------------------
 !!     DECLARATIONS
@@ -80,7 +81,7 @@ INTEGER, PARAMETER :: JP_AER_DST = 7
 
 INTEGER            :: NSOA = 10    ! number of condensable species that may form
                                    ! secondary aerosols
-INTEGER            :: NM6_AER = 2  ! number of condensable species that may form
+INTEGER, SAVE      :: NM6_AER ! number of mode for which M6 is computed define in ini_sv
                                    ! secondary aerosols
 INTEGER            :: JP_AER_SOA1 = 8 
 INTEGER            :: JP_AER_SOA2 = 9
diff --git a/src/MNH/modd_ch_mnhcn.f90 b/src/MNH/modd_ch_mnhcn.f90
index a6d7afc0f5f858836f0097b149073f9f507fbb28..ddf17da788dd3e17866d10f9508354609bb770f3 100644
--- a/src/MNH/modd_ch_mnhcn.f90
+++ b/src/MNH/modd_ch_mnhcn.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$ $Date$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modd_ch_mnhcn.f90,v $ $Revision: 1.2.4.1.2.1.12.2 $ $Date: 2014/01/09 15:01:56 $
 !-----------------------------------------------------------------
 !!    #####################
       MODULE MODD_CH_MNHC_n
@@ -35,6 +35,7 @@
 !!    25/04/08 (M. Leriche) add threshold for aqueous phase chemistry
 !!    16/09/10 (M. Leriche) add flag for ice phase chemistry
 !!    13/01/11 (M. Leriche) add flag for retention in ice 
+!!    24/03/16 (M. Leriche) remove surface option -> manage them in SURFEX
 !!
 !!    IMPLICIT ARGUMENTS
 !!    ------------------
@@ -64,21 +65,15 @@ TYPE CH_MNHC_t
 !* Initialization
 !
   LOGICAL :: LCH_INIT_FIELD ! flag indicating whether initialization
-		 ! of chemical fields shall be done during MesoNH run using
-		 ! CH_INIT_FIELD (overwrites initial values from FM-files)
-		 ! or not
-!
-!* Surface options
-!
-  LOGICAL :: LCH_SURFACE_FLUX  ! flag indicating whether surface flux
-		 ! for chemical species shall be calculated using 
-		 ! CH_SURFACE_FLUX or not (dry deposition and emission)
+                 ! of chemical fields shall be done during MesoNH run using
+                 ! CH_INIT_FIELD (overwrites initial values from FM-files)
+                 ! or not
 !
 !* Scavenging in parameterized convective clouds
 !
   LOGICAL :: LCH_CONV_SCAV 
                  ! flag for calculation of scavenging 
-		 ! by convective precipitations (active only if LCHTRANS=.TRUE.)
+                 ! by convective precipitations (active only if LCHTRANS=.TRUE.)
 !
 !* pH calculation
 !
@@ -104,23 +99,23 @@ TYPE CH_MNHC_t
 !
   CHARACTER(LEN=80) :: CCHEM_INPUT_FILE 
                  ! name of general 
-		 ! purpose ASCII input file (handeled by CH_OPEN_INPUT)
+                 ! purpose ASCII input file (handeled by CH_OPEN_INPUT)
 !
   CHARACTER(LEN=10) :: CCH_TDISCRETIZATION 
-		 ! temporal discretization:
+                 ! temporal discretization:
                  ! "SPLIT"  : use time-splitting, input fields for solver are
-		 !            scalar variables at t+dt (derived from XRSVS)
-		 ! "CENTER" : input fields for solver are 
-		 !            scalar variables at t (XSVT)
-		 ! "LAGGED" : input fields for solver are 
-		 !            scalar variables at t-dt (XSVM)
+                 !            scalar variables at t+dt (derived from XRSVS)
+                 ! "CENTER" : input fields for solver are 
+                 !            scalar variables at t (XSVT)
+                 ! "LAGGED" : input fields for solver are 
+                 !            scalar variables at t-dt (XSVM)
 !
   INTEGER           :: NCH_SUBSTEPS
                  ! number of chemical timesteps to be taken during one 
-		 ! double timestep of MesoNH (MesoNH integrates with timesteps
-		 ! of lenght 2*XTSTEP using leapfrog), the timestep of the 
-		 ! solver will be calculated as 
-		 ! ZDTSOLVER = 2*XTSTEP / NCH_SUBSTEPS
+                 ! double timestep of MesoNH (MesoNH integrates with timesteps
+                 ! of lenght 2*XTSTEP using leapfrog), the timestep of the 
+                 ! solver will be calculated as 
+                 ! ZDTSOLVER = 2*XTSTEP / NCH_SUBSTEPS
 !* LiNOx
 !
   LOGICAL :: LCH_CONV_LINOX
@@ -129,8 +124,8 @@ TYPE CH_MNHC_t
 !* photolysis rates (TUV)
 !
   LOGICAL      :: LCH_TUV_ONLINE  ! switch online/lookup table
-  CHARACTER*80 :: CCH_TUV_LOOKUP  ! name of lookup table file
-  CHARACTER*4  :: CCH_TUV_CLOUDS  ! method for calculating the
+  CHARACTER(LEN=80) :: CCH_TUV_LOOKUP  ! name of lookup table file
+  CHARACTER(LEN=4)  :: CCH_TUV_CLOUDS  ! method for calculating the
                                 ! impact of clouds on radiation
                                 ! "FOUQ" (model clouds, only 1-D)
   REAL :: XCH_TUV_ALBNEW  ! surface albedo (if negative the albedo
@@ -144,18 +139,18 @@ TYPE CH_MNHC_t
 !
 !* vectorization
 !
-  CHARACTER*3 :: CCH_VEC_METHOD          ! type of vectorization mask
+  CHARACTER(LEN=3) :: CCH_VEC_METHOD          ! type of vectorization mask
                                        ! 'MAX' take NCH_VEC_LENGTH points
                                        ! 'TOT' take all grid points
                                        ! 'HOR' take horizontal layers
                                        ! 'VER' take vertical columns
-  INTEGER     :: NCH_VEC_LENGTH          ! number of points for 'MAX' option
+  INTEGER          :: NCH_VEC_LENGTH          ! number of points for 'MAX' option
 !
 !* 1-D time series
 !
-  REAL         :: XCH_TS1D_TSTEP         ! time between two call to write_ts1d
-  CHARACTER*80 :: CCH_TS1D_COMMENT       ! comment for write_ts1d
-  CHARACTER*80 :: CCH_TS1D_FILENAME      ! filename for write_ts1d files
+  REAL              :: XCH_TS1D_TSTEP         ! time between two call to write_ts1d
+  CHARACTER(LEN=80) :: CCH_TS1D_COMMENT       ! comment for write_ts1d
+  CHARACTER(LEN=80) :: CCH_TS1D_FILENAME      ! filename for write_ts1d files
 !
 END TYPE CH_MNHC_t
 
@@ -165,7 +160,6 @@ LOGICAL, POINTER :: LUSECHEM=>NULL()
 LOGICAL, POINTER :: LUSECHAQ=>NULL()
 LOGICAL, POINTER :: LUSECHIC=>NULL()
 LOGICAL, POINTER :: LCH_INIT_FIELD=>NULL()
-LOGICAL, POINTER :: LCH_SURFACE_FLUX=>NULL()
 LOGICAL, POINTER :: LCH_CONV_SCAV=>NULL()
 LOGICAL, POINTER :: LCH_PH=>NULL()
 LOGICAL, POINTER :: LCH_RET_ICE=>NULL()
@@ -177,16 +171,16 @@ CHARACTER(LEN=10), POINTER :: CCH_TDISCRETIZATION=>NULL()
 INTEGER, POINTER :: NCH_SUBSTEPS=>NULL()
 LOGICAL, POINTER :: LCH_CONV_LINOX=>NULL()
 LOGICAL, POINTER :: LCH_TUV_ONLINE=>NULL()
-CHARACTER*80, POINTER :: CCH_TUV_LOOKUP=>NULL()
-CHARACTER*4, POINTER :: CCH_TUV_CLOUDS=>NULL()
+CHARACTER(LEN=80), POINTER :: CCH_TUV_LOOKUP=>NULL()
+CHARACTER(LEN=4), POINTER :: CCH_TUV_CLOUDS=>NULL()
 REAL, POINTER :: XCH_TUV_ALBNEW=>NULL()
 REAL, POINTER :: XCH_TUV_DOBNEW=>NULL()
 REAL, POINTER :: XCH_TUV_TUPDATE=>NULL()
-CHARACTER*3, POINTER :: CCH_VEC_METHOD=>NULL()
+CHARACTER(LEN=3), POINTER :: CCH_VEC_METHOD=>NULL()
 INTEGER, POINTER :: NCH_VEC_LENGTH=>NULL()
 REAL, POINTER :: XCH_TS1D_TSTEP=>NULL()
-CHARACTER*80, POINTER :: CCH_TS1D_COMMENT=>NULL()
-CHARACTER*80, POINTER :: CCH_TS1D_FILENAME=>NULL()
+CHARACTER(LEN=80), POINTER :: CCH_TS1D_COMMENT=>NULL()
+CHARACTER(LEN=80), POINTER :: CCH_TS1D_FILENAME=>NULL()
 
 CONTAINS
 
@@ -200,7 +194,6 @@ LUSECHEM=>CH_MNHC_MODEL(KTO)%LUSECHEM
 LUSECHAQ=>CH_MNHC_MODEL(KTO)%LUSECHAQ
 LUSECHIC=>CH_MNHC_MODEL(KTO)%LUSECHIC
 LCH_INIT_FIELD=>CH_MNHC_MODEL(KTO)%LCH_INIT_FIELD
-LCH_SURFACE_FLUX=>CH_MNHC_MODEL(KTO)%LCH_SURFACE_FLUX
 LCH_CONV_SCAV=>CH_MNHC_MODEL(KTO)%LCH_CONV_SCAV
 LCH_PH=>CH_MNHC_MODEL(KTO)%LCH_PH
 LCH_RET_ICE=>CH_MNHC_MODEL(KTO)%LCH_RET_ICE
diff --git a/src/MNH/modd_ch_phn.f90 b/src/MNH/modd_ch_phn.f90
index 8e5f6230b69b119c95092942926aa1003a0f68ce..f028017567cd79e4cca06538764411afc176f037 100644
--- a/src/MNH/modd_ch_phn.f90
+++ b/src/MNH/modd_ch_phn.f90
@@ -28,6 +28,7 @@
 !!    -------------
 !!    Original 01/06/07
 !!       P. Tulet      Nov 2014 accumulated moles of aqueous species that fall at the surface   
+!!       P. Tulet & M. Leriche Nov 2015 add pH in rain at the surface
 !!
 !!    IMPLICIT ARGUMENTS
 !!    ------------------
@@ -42,10 +43,12 @@ IMPLICIT NONE
 
 TYPE CH_PH_t
 !
+
   REAL, POINTER, DIMENSION(:,:,:) :: XPHC ! cloud
   REAL, POINTER, DIMENSION(:,:,:) :: XPHR ! rain
   REAL, POINTER, DIMENSION(:,:,:) :: XACPRAQ ! sum of aqueous chemical species fall at the surface by rain
                                              ! in moles i / m2 (ratio with XACPRR for concentration
+  REAL, POINTER, DIMENSION(:,:) :: XACPHR !  mean PH in accumulated surface rain
 !
 !-----------------------------------------------------------------------------
 END TYPE CH_PH_t
@@ -54,6 +57,7 @@ TYPE(CH_PH_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: CH_PH_MODEL
 
 REAL, POINTER, DIMENSION(:,:,:) :: XPHC=>NULL()
 REAL, POINTER, DIMENSION(:,:,:) :: XPHR=>NULL()
+REAL, POINTER, DIMENSION(:,:) :: XACPHR=>NULL()
 REAL, POINTER, DIMENSION(:,:,:) :: XACPRAQ=>NULL()
 
 CONTAINS
@@ -64,11 +68,13 @@ INTEGER, INTENT(IN) :: KFROM, KTO
 ! Save current state for allocated arrays
 CH_PH_MODEL(KFROM)%XPHC=>XPHC
 CH_PH_MODEL(KFROM)%XPHR=>XPHR
+CH_PH_MODEL(KFROM)%XACPHR=>XACPHR
 CH_PH_MODEL(KFROM)%XACPRAQ=>XACPRAQ
 !
 ! Current model is set to model KTO
 XPHC=>CH_PH_MODEL(KTO)%XPHC
 XPHR=>CH_PH_MODEL(KTO)%XPHR
+XACPHR=>CH_PH_MODEL(KTO)%XACPHR
 XACPRAQ=>CH_PH_MODEL(KTO)%XACPRAQ
 
 END SUBROUTINE CH_PH_GOTO_MODEL
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index 3e2fc2332f8351e84f2d932db2cf7381551e43d2..08911c4e1a22dbb68a132cead29c91d0d486f736 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -237,6 +237,8 @@ END MODULE MODI_MODEL_n
 !!                              of write_phys_param
 !!      J.Escobar : 19/04/2016 : Pb IOZ/NETCDF , missing OPARALLELIO=.FALSE. for PGD files
 !!      M.Mazoyer : 04/2016      DTHRAD used for radiative cooling when LACTIT
+!!      M.Leriche : 03/2016 Move computation of accumulated chem. in rain to ch_monitor
+!!                  09/2016 Add filter on negative values on AERDEP SV before relaxation
 !!-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -303,7 +305,6 @@ USE MODD_SERIES_n, ONLY: NFREQSERIES
 USE MODD_CH_AERO_n,    ONLY: XSOLORG, XMI
 USE MODD_CH_MNHC_n,    ONLY: LUSECHEM,LCH_CONV_LINOX,LUSECHAQ,LUSECHIC, &
                              LCH_INIT_FIELD
-USE MODD_CH_PH_n
 USE MODD_CST, ONLY: XMD
 USE MODD_NUDGING_n
 USE MODD_PARAM_MFSHALL_n
@@ -1218,6 +1219,9 @@ END DO
 DO JSV = NSV_SLTDEPBEG,NSV_SLTDEPEND
   XRSVS(:,:,:,JSV) = MAX(XRSVS(:,:,:,JSV),0.)
 END DO
+DO JSV = NSV_AERDEPBEG,NSV_AERDEPEND
+  XRSVS(:,:,:,JSV) = MAX(XRSVS(:,:,:,JSV),0.)
+END DO
 
 IF (CELEC .NE. 'NONE') THEN
   XRSVS(:,:,:,NSV_ELECBEG) = MAX(XRSVS(:,:,:,NSV_ELECBEG),0.)
@@ -1732,15 +1736,6 @@ IF (CCLOUD /= 'NONE' .AND. CELEC == 'NONE') THEN
 !
   IF (CCLOUD /= 'REVE' ) THEN
     XACPRR = XACPRR + XINPRR * XTSTEP
-      IF (LUSECHAQ) THEN
-      DO JSV=1,NSV_CHAC/2
-      WHERE(XRT(:,:,IKB,3) .GT. 0.)
-      XACPRAQ(:,:,JSV) = XACPRAQ(:,:,JSV) + &
-              (XSVT(:,:,IKB,JSV+NSV_CHACBEG+NSV_CHAC/2-1))/ (XMD*XRT(:,:,IKB,3))*& ! moles i  / kg eau
-               XINPRR(:,:) * XTSTEP ! moles i / m2
-      END WHERE
-      END DO
-      END IF
     IF ((CCLOUD(1:3) == 'ICE' .AND. LSEDIC ) .OR.                       &
         ((CCLOUD == 'C2R2' .OR. CCLOUD == 'C3R5' .OR. CCLOUD == 'KHKO') &
                               .AND. LSEDC  )      )   THEN                  
@@ -1785,8 +1780,8 @@ IF (CELEC /= 'NONE' .AND. (CCLOUD(1:3) == 'ICE')) THEN
     CALL MNHGET_SURF_PARAM_n (PSEA=ZSEA(:,:),PTOWN=ZTOWN(:,:))
     CALL RESOLVED_ELEC_n (CCLOUD, CSCONV, CMF_CLOUD,                     &
                           NRR, NSPLITR, IMI, KTCOUNT, OEXIT,             &
-                          CLBCX, CLBCY, CRAD, CTURBDIM,                  &
-                          LSUBG_COND, LSIGMAS,VSIGQSAT,CSUBG_AUCV,       &
+                          CLBCX, CLBCY, YFMFILE, CLUOUT, CRAD, CTURBDIM, &
+                          GCLOSE_OUT, LSUBG_COND, LSIGMAS,VSIGQSAT,CSUBG_AUCV,   &
                           XTSTEP, XZZ, XRHODJ, XRHODREF, XEXNREF,        &
                           XPABST, XTHT, XRTHS, XWT,  XRT, XRRS,          &
                           XSVT, XRSVS, XCIT,                             &
@@ -1799,8 +1794,8 @@ IF (CELEC /= 'NONE' .AND. (CCLOUD(1:3) == 'ICE')) THEN
   ELSE
     CALL RESOLVED_ELEC_n (CCLOUD, CSCONV, CMF_CLOUD,                     &
                           NRR, NSPLITR, IMI, KTCOUNT, OEXIT,             &
-                          CLBCX, CLBCY, CRAD, CTURBDIM,                  &
-                          LSUBG_COND, LSIGMAS,VSIGQSAT, CSUBG_AUCV,      &
+                          CLBCX, CLBCY, YFMFILE, CLUOUT, CRAD, CTURBDIM, &
+                          GCLOSE_OUT, LSUBG_COND, LSIGMAS,VSIGQSAT, CSUBG_AUCV,   &
                           XTSTEP, XZZ, XRHODJ, XRHODREF, XEXNREF,        &
                           XPABST, XTHT, XRTHS, XWT,                      &
                           XRT, XRRS, XSVT, XRSVS, XCIT,                  &
diff --git a/src/MNH/modn_ch_mnhcn.f90 b/src/MNH/modn_ch_mnhcn.f90
index 892876874d078fe63d2cca65b4046b12c2bd622e..4450b9a8b5544b85c0c535a5dc618c8a685ad084 100644
--- a/src/MNH/modn_ch_mnhcn.f90
+++ b/src/MNH/modn_ch_mnhcn.f90
@@ -5,7 +5,7 @@
 !-----------------------------------------------------------------
 !--------------- special set of characters for RCS information
 !-----------------------------------------------------------------
-! $Source$ $Revision$ $Date$
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/modn_ch_mnhcn.f90,v $ $Revision: 1.2.4.1.2.1.12.2 $ $Date: 2014/01/09 15:01:56 $
 !-----------------------------------------------------------------
 !!    #####################
       MODULE MODN_CH_MNHC_n
@@ -32,6 +32,7 @@
 !!    25/04/08 (M. Leriche) add threshold for aqueous phase chemistry
 !!    16/09/10 (M. Leriche) add logical for ice phase chemistry
 !!    13/01/11 (M. Leriche) add logical for retention in ice
+!!    24/03/16 (M. Leriche) remove surface option -> manage them in SURFEX
 !!
 !!    IMPLICIT ARGUMENTS
 !!    ------------------
@@ -40,7 +41,6 @@ USE MODD_CH_MNHC_n, ONLY: &
          LUSECHAQ_n => LUSECHAQ, &
          LUSECHIC_n => LUSECHIC, &
          LCH_INIT_FIELD_n => LCH_INIT_FIELD, &
-         LCH_SURFACE_FLUX_n => LCH_SURFACE_FLUX, &
          LCH_CONV_SCAV_n => LCH_CONV_SCAV, &
          LCH_CONV_LINOX_n => LCH_CONV_LINOX, &
          LCH_PH_n => LCH_PH, &
@@ -68,7 +68,6 @@ LOGICAL  :: LUSECHEM
 LOGICAL  :: LUSECHAQ
 LOGICAL  :: LUSECHIC
 LOGICAL  :: LCH_INIT_FIELD
-LOGICAL  :: LCH_SURFACE_FLUX
 LOGICAL  :: LCH_CONV_SCAV
 LOGICAL  :: LCH_CONV_LINOX
 LOGICAL  :: LCH_PH
@@ -90,13 +89,13 @@ REAL  :: XCH_TS1D_TSTEP
 CHARACTER*80  :: CCH_TS1D_COMMENT
 CHARACTER*80  :: CCH_TS1D_FILENAME
 !
-NAMELIST/NAM_CH_MNHCn/LUSECHEM,LUSECHAQ,LUSECHIC,LCH_INIT_FIELD,LCH_SURFACE_FLUX,&
-                      LCH_CONV_SCAV,LCH_CONV_LINOX,LCH_PH,LCH_RET_ICE,XCH_PHINIT,&
-                      XRTMIN_AQ,CCHEM_INPUT_FILE,CCH_TDISCRETIZATION,NCH_SUBSTEPS,&
-                      LCH_TUV_ONLINE,CCH_TUV_LOOKUP,CCH_TUV_CLOUDS,XCH_TUV_ALBNEW,&
-                      XCH_TUV_DOBNEW,XCH_TUV_TUPDATE,CCH_VEC_METHOD,&
-                      NCH_VEC_LENGTH,XCH_TS1D_TSTEP,CCH_TS1D_COMMENT,&
-                      CCH_TS1D_FILENAME
+NAMELIST/NAM_CH_MNHCn/LUSECHEM,LUSECHAQ,LUSECHIC,LCH_INIT_FIELD,LCH_CONV_SCAV,&
+                      LCH_CONV_LINOX,LCH_PH,LCH_RET_ICE,XCH_PHINIT,XRTMIN_AQ, &
+                      CCHEM_INPUT_FILE,CCH_TDISCRETIZATION,NCH_SUBSTEPS,      &
+                      LCH_TUV_ONLINE,CCH_TUV_LOOKUP,CCH_TUV_CLOUDS,           &
+                      XCH_TUV_ALBNEW,XCH_TUV_DOBNEW,XCH_TUV_TUPDATE,          &
+                      CCH_VEC_METHOD,NCH_VEC_LENGTH,XCH_TS1D_TSTEP,           &
+                      CCH_TS1D_COMMENT,CCH_TS1D_FILENAME
 !
 CONTAINS
 !
@@ -105,7 +104,6 @@ SUBROUTINE INIT_NAM_CH_MNHCn
   LUSECHAQ = LUSECHAQ_n
   LUSECHIC = LUSECHIC_n
   LCH_INIT_FIELD = LCH_INIT_FIELD_n
-  LCH_SURFACE_FLUX = LCH_SURFACE_FLUX_n
   LCH_CONV_SCAV = LCH_CONV_SCAV_n
   LCH_CONV_LINOX = LCH_CONV_LINOX_n
   LCH_PH = LCH_PH_n
@@ -133,7 +131,6 @@ SUBROUTINE UPDATE_NAM_CH_MNHCn
   LUSECHAQ_n = LUSECHAQ
   LUSECHIC_n = LUSECHIC
   LCH_INIT_FIELD_n = LCH_INIT_FIELD
-  LCH_SURFACE_FLUX_n = LCH_SURFACE_FLUX
   LCH_PH_n = LCH_PH
   LCH_RET_ICE_n = LCH_RET_ICE
   XCH_PHINIT_n = XCH_PHINIT
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index 049dab242b00264804ff5d995e56eb29059fbcb6..881d74b3b9e610593814513ab68b7a23e9a2e231 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -280,6 +280,7 @@ END MODULE MODI_READ_EXSEG_n
 !!      Modification   01/2015   (C. Barthe) add explicit LNOx
 !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
 !!      M.Leriche 18/12/2015 : bug chimie glace dans prep_real_case
+!!      Modification   02/2016   (M.Leriche) treat gas and aq. chemicals separately
 !!------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -927,8 +928,8 @@ LUSETKE(KMI) = (CTURB /= 'NONE')
 !*       2.3     Chemical and NSV_* variables initializations
 !
 CALL UPDATE_NAM_PARAMN
-CALL UPDATE_NAM_CH_MNHCN
 CALL UPDATE_NAM_DYNN
+CALL UPDATE_NAM_CONFN
 !
 IF (LORILAM .AND. .NOT. LUSECHEM) THEN
   WRITE(UNIT=ILUOUT,FMT=9002) KMI
@@ -1012,7 +1013,7 @@ IF (LUSECHEM) THEN
 END IF
 !
 
-CALL UPDATE_NAM_CONFN
+CALL UPDATE_NAM_CH_MNHCN
 CALL INI_NSV(KMI)
 !
 ! From this point, all NSV* variables contain valid values for model KMI
@@ -1376,18 +1377,31 @@ IF (CELEC /= 'NONE' .AND. LLNOX_EXPLICIT) THEN
   END IF
 END IF
 !
-! Chemical SV case (including aqueous chemical species)
+! Chemical SV case (excluding aqueous chemical species)
 !
 IF (LUSECHEM) THEN
   IF (OUSECHEM) THEN
-    CGETSVT(NSV_CHEMBEG:NSV_CHEMEND)='READ'
-    IF(CCONF=='START' .AND. LCH_INIT_FIELD ) CGETSVT(NSV_CHEMBEG:NSV_CHEMEND)='INIT'
+    CGETSVT(NSV_CHGSBEG:NSV_CHGSEND)='READ'
+    IF(CCONF=='START' .AND. LCH_INIT_FIELD ) CGETSVT(NSV_CHGSBEG:NSV_CHGSEND)='INIT'
   ELSE
     WRITE(UNIT=ILUOUT,FMT=9001) KMI
     WRITE(UNIT=ILUOUT,FMT='("THERE IS NO SCALAR VARIABLES FOR CHEMICAL &
          &SCHEME IN INITIAL FMFILE",/,&
          & "THE CHEMICAL VARIABLES HAVE BEEN INITIALIZED TO ZERO ")') 
-    CGETSVT(NSV_CHEMBEG:NSV_CHEMEND)='INIT'
+    CGETSVT(NSV_CHGSBEG:NSV_CHGSEND)='INIT'
+  END IF
+END IF
+! add aqueous chemical species
+IF (LUSECHAQ) THEN
+  IF (OUSECHAQ) THEN
+    CGETSVT(NSV_CHACBEG:NSV_CHACEND)='READ'
+!    IF(CCONF=='START') CGETSVT(NSV_CHACBEG:NSV_CHACEND)='INIT'
+  ELSE
+    WRITE(UNIT=ILUOUT,FMT=9001) KMI
+    WRITE(UNIT=ILUOUT,FMT='("THERE IS NO SCALAR VARIABLES FOR CHEMICAL &
+         &SCHEME IN AQUEOUS PHASE IN INITIAL FMFILE",/,&
+         & "THE AQUEOUS PHASE CHEMICAL VARIABLES HAVE BEEN INITIALIZED TO ZERO ")') 
+    CGETSVT(NSV_CHACBEG:NSV_CHACEND)='INIT'
   END IF
 END IF
 ! add ice phase chemical species
diff --git a/src/MNH/read_field.f90 b/src/MNH/read_field.f90
index e95c80676abe4e5e872cbaf9521d88933c2f9bfd..07a16e6f83bc6989d2786205a545585f155d842a 100644
--- a/src/MNH/read_field.f90
+++ b/src/MNH/read_field.f90
@@ -230,6 +230,7 @@ END MODULE MODI_READ_FIELD
 !!          Bosseur & Filippi 07/13 Adds Forefire
 !!          M. Leriche  11/14     correct bug in pH initialization
 !!          C.Lac       12/14     correction for reproducibility START/RESTA
+!!          M. Leriche  02/16     treat gas and aq. chemicals separately
 !!-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -666,11 +667,24 @@ DO JSV = NSV_ELECBEG,NSV_ELECEND
   END SELECT
 END DO
 !
-DO JSV = NSV_CHEMBEG,NSV_CHEMEND
+DO JSV = NSV_CHGSBEG,NSV_CHGSEND
   SELECT CASE(HGETSVT(JSV))
   CASE ('READ')
-    CNAMES(JSV-NSV_CHEMBEG+1) = UPCASE(CNAMES(JSV-NSV_CHEMBEG+1))
-    YRECFM=TRIM(CNAMES(JSV-NSV_CHEMBEG+1))//'T'
+    CNAMES(JSV-NSV_CHGSBEG+1) = UPCASE(CNAMES(JSV-NSV_CHGSBEG+1))
+    YRECFM=TRIM(CNAMES(JSV-NSV_CHGSBEG+1))//'T'
+    CALL FMREAD(HINIFILE,YRECFM,HLUOUT,YDIR,Z3D,IGRID,ILENCH,  &
+         YCOMMENT,IRESP)
+    PSVT(:,:,:,JSV) = Z3D(:,:,:)
+  CASE ('INIT')
+    PSVT(:,:,:,JSV) = 0.
+  END SELECT    
+END DO
+!
+DO JSV = NSV_CHACBEG,NSV_CHACEND
+  SELECT CASE(HGETSVT(JSV))
+  CASE ('READ')
+    CNAMES(JSV-NSV_CHACBEG+NSV_CHGS+1) = UPCASE(CNAMES(JSV-NSV_CHACBEG+NSV_CHGS+1))
+    YRECFM=TRIM(CNAMES(JSV-NSV_CHACBEG+NSV_CHGS+1))//'T'
     CALL FMREAD(HINIFILE,YRECFM,HLUOUT,YDIR,Z3D,IGRID,ILENCH,  &
          YCOMMENT,IRESP)
     PSVT(:,:,:,JSV) = Z3D(:,:,:)
diff --git a/src/MNH/write_aircraft_balloon.f90 b/src/MNH/write_aircraft_balloon.f90
index 3060afd81857f75e954c360420119884f28e8f1a..8771a99fd74bc2a9c7731d2bd60ae78897603fb6 100644
--- a/src/MNH/write_aircraft_balloon.f90
+++ b/src/MNH/write_aircraft_balloon.f90
@@ -61,8 +61,9 @@ END MODULE MODI_WRITE_AIRCRAFT_BALLOON
 !!     Original 15/05/2000
 !!     10/01/2011 adding IMI, the model number
 !!     March, 2013 :  C.Lac : add vertical profiles
-!!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
+!!     July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
 !!                                      aircraft, ballon and profiler
+!!     August 2016 (M.Leriche) Add mass concentration of aerosol species
 !!
 !! --------------------------------------------------------------------------
 !       
@@ -75,7 +76,13 @@ USE MODD_PARAMETERS
 !
 USE MODD_AIRCRAFT_BALLOON
 USE MODD_CH_M9_n,         ONLY: CNAMES
-USE MODD_CH_AEROSOL,      ONLY: CAERONAMES, LORILAM, JPMODE
+USE MODD_CH_AEROSOL,      ONLY: CAERONAMES, LORILAM, NSP, NCARB, NSOA,    & 
+                                JPMODE, JP_AER_BC, JP_AER_OC, JP_AER_DST, &
+                                JP_AER_H2O, JP_AER_SO4, JP_AER_NO3,       &
+                                JP_AER_NH3, JP_AER_SOA1, JP_AER_SOA2,     &
+                                JP_AER_SOA3, JP_AER_SOA4, JP_AER_SOA5,    &
+                                JP_AER_SOA6, JP_AER_SOA7, JP_AER_SOA8,    &
+                                JP_AER_SOA9, JP_AER_SOA10
 USE MODD_RAIN_C2R2_DESCR, ONLY: C2R2NAMES
 USE MODD_ICE_C1R3_DESCR,  ONLY: C1R3NAMES
 USE MODD_ELEC_DESCR,      ONLY: CELECNAMES
@@ -172,6 +179,7 @@ REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ZW6    ! contains temporal serie to
 REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ZWORKZ6! contains temporal serie
 REAL, DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: ZWZ6   ! contains temporal serie
 REAL, DIMENSION(:,:,:,:),     ALLOCATABLE :: ZSV, ZN0, ZSIG, ZRG
+REAL, DIMENSION(:,:,:,:,:),   ALLOCATABLE :: ZPTOTA
 REAL, DIMENSION(:,:,:),       ALLOCATABLE :: ZRHO
 !
 INTEGER, DIMENSION(:),            ALLOCATABLE :: IGRID    ! grid indicator
@@ -460,6 +468,7 @@ IF (SIZE(TPFLYER%SV,2)>=1) THEN
     ALLOCATE (ZN0(1,1,SIZE(TPFLYER%TIME),JPMODE)) 
     ALLOCATE (ZRG(1,1,SIZE(TPFLYER%TIME),JPMODE)) 
     ALLOCATE (ZSIG(1,1,SIZE(TPFLYER%TIME),JPMODE)) 
+    ALLOCATE (ZPTOTA(1,1,SIZE(TPFLYER%TIME),NSP+NCARB+NSOA,JPMODE))    
     ZSV(1,1,:,1:NSV_AER) = TPFLYER%SV(:,NSV_AERBEG:NSV_AEREND)
     IF (SIZE(TPFLYER%R,2) >0) THEN
       ZRHO(1,1,:) = 0.
@@ -473,7 +482,15 @@ IF (SIZE(TPFLYER%SV,2)>=1) THEN
     ENDIF
     ZRHO(1,1,:) =  TPFLYER%P(:) / &
                   (XRD *ZRHO(1,1,:) *((TPFLYER%P(:)/XP00)**(XRD/XCPD))  )
-    CALL PPP2AERO(ZSV,ZRHO, PSIG3D=ZSIG, PRG3D=ZRG, PN3D=ZN0)
+    ZSIG = 0.
+    ZRG = 0.
+    ZN0 = 0.
+    ZPTOTA = 0.
+    DO JPT=1,SIZE(TPFLYER%TIME) ! prevent division by zero if ZSV = 0.
+      IF (ALL(ZSV(1,1,JPT,:)/=0.)) THEN
+        CALL PPP2AERO(ZSV,ZRHO, PSIG3D=ZSIG, PRG3D=ZRG, PN3D=ZN0, PCTOTA=ZPTOTA)
+      ENDIF
+    ENDDO
     DO JSV=1,JPMODE
       ! mean radius
       JPROC = JPROC+1
@@ -493,9 +510,114 @@ IF (SIZE(TPFLYER%SV,2)>=1) THEN
       YUNIT    (JPROC) = '  '
       WRITE(YCOMMENT(JPROC),'(A13,I1,A6)')'N0 AERO MODE ',JSV,' (1/m3)'
       ZWORK6 (1,1,1,:,1,JPROC) = ZN0(1,1,:,JSV)
+      ! mass concentration in microg/m3
+      ! sulfate
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MSO4',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SO4 AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SO4,JSV)
+      ! nitrate
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MNO3',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS NO3 AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_NO3,JSV)
+      ! amoniac
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MNH3',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS NH3 AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_NH3,JSV)
+      ! water
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MH2O',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS H2O AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_H2O,JSV)
+      IF (NSOA .EQ. 10) THEN
+        ! SOA1
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA1',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA1 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA1,JSV)
+        ! SOA2
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA2',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA2 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA2,JSV)
+        ! SOA3
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA3',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA3 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA3,JSV)
+        ! SOA4
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA4',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA4 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA4,JSV)
+        ! SOA5
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA5',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA5 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA5,JSV)
+        ! SOA6
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA6',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA6 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA6,JSV)
+        ! SOA7
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA7',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA7 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA7,JSV)
+        ! SOA8
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA8',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA8 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA8,JSV)
+        ! SOA9
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA9',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA9 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA9,JSV)
+        ! SOA10
+        JPROC = JPROC + 1
+        WRITE(YTITLE(JPROC),'(A4,I1)')'MSOA10',JSV
+        YUNIT    (JPROC) = 'ug/m3'
+        WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS SOA10 AEROSOL MODE ',JSV,'(ug/m3)'
+        ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_SOA10,JSV)
+      ENDIF
+      ! OC
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MOC',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS OC AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_OC,JSV)
+      ! BC
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MBC',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS BC AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_BC,JSV)
+      ! dust
+      JPROC = JPROC + 1
+      WRITE(YTITLE(JPROC),'(A4,I1)')'MDUST',JSV
+      YUNIT    (JPROC) = 'ug/m3'
+      WRITE(YCOMMENT(JPROC),'(A22,I1,A5)')'MASS DUST AEROSOL MODE ',JSV,'(ug/m3)'
+      ZWORK6(1,1,1,:,1,JPROC) = ZPTOTA(1,1,:,JP_AER_DST,JSV)      
     ENDDO
     DEALLOCATE (ZSV,ZRHO) 
-    DEALLOCATE (ZN0,ZRG,ZSIG) 
+    DEALLOCATE (ZN0,ZRG,ZSIG,ZPTOTA) 
   END IF
 ! dust scalar variables
   DO JSV = NSV_DSTBEG,NSV_DSTEND
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index bf3cba65a067be919f74a840f314193da417dd9e..38918f6659015d3d30e6458f1cf25c1bfecedfba 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -160,6 +160,7 @@ END MODULE MODI_WRITE_LFIFM_n
 !!       C.Lac         Dec.2014 writing past wind fields for centred advection
 !!       J.-P. Pinty   Jan 2015 add LNOx and flash map diagnostics
 !!       J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!!       P. Tulet & M. Leriche    Nov 2015 add mean pH value in the rain at the surface
 !!       J.escobar     04/08/2015 suit Pb with writ_lfin JSA increment , modif in ini_nsv to have good order initialization
 !!                   
 !-------------------------------------------------------------------------------
@@ -198,7 +199,7 @@ USE MODD_NESTING
 USE MODD_PARAMETERS
 USE MODD_GR_FIELD_n
 USE MODD_CH_MNHC_n,       ONLY: LUSECHEM,LCH_CONV_LINOX, &
-                                LUSECHAQ,LUSECHIC,LCH_PH
+                                LUSECHAQ,LUSECHIC,LCH_PH, XCH_PHINIT
 USE MODD_CH_PH_n
 USE MODD_CH_M9_n
 USE MODD_RAIN_C2R2_DESCR, ONLY: C2R2NAMES
@@ -1171,6 +1172,20 @@ IF (NSV >=1) THEN
         YCOMMENT='X_Y_Z_PHR'
         ILENCH=LEN(YCOMMENT)
         CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,XPHR,IGRID,ILENCH,YCOMMENT,IRESP)
+        ! compute mean pH in accumulated surface water
+        !ZWORK2D(:,:) = 10**(-XCH_PHINIT)
+        WHERE (XACPRR > 0.)
+        ZWORK2D(:,:) =  XACPHR(:,:) *1E3 / XACPRR(:,:) ! moles of H+ / l of water 
+        ELSE WHERE
+        ZWORK2D(:,:) = XUNDEF
+        END WHERE
+        WHERE ((ZWORK2D(:,:) < 1E-1).AND.(ZWORK2D(:,:) > 1E-14))
+        ZWORK2D(:,:) = -ALOG10(ZWORK2D(:,:))           ! mean pH of surface water
+        END WHERE
+        YRECFM = 'MEANPHR'
+        YCOMMENT='X_Y_MEAN_PH'
+        ILENCH=LEN(YCOMMENT)
+        CALL FMWRIT(HFMFILE,YRECFM,CLUOUT,YDIR,ZWORK2D,IGRID,ILENCH,YCOMMENT,IRESP)        
       ENDIF
     ENDIF
   ELSE IF (LCH_CONV_LINOX) THEN
diff --git a/src/SURFEX/ch_init_depconst.F90 b/src/SURFEX/ch_init_depconst.F90
index ca3513297ee871cfef745f7d1ef3b515b992ed08..e2cbd8c49959ba1f15f8aa9c904ba053388c16d4 100644
--- a/src/SURFEX/ch_init_depconst.F90
+++ b/src/SURFEX/ch_init_depconst.F90
@@ -3,7 +3,7 @@
 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !SURFEX_LIC for details. version 1.
 !     #########
-       SUBROUTINE CH_INIT_DEPCONST(KCH,KLUOUT,HSV)
+       SUBROUTINE CH_INIT_DEPCONST(HPROGRAM,HCHEM_SURF_FILE,KLUOUT,HSV)
 !!    ##################################################
 !!
 !!*** *CH_INIT_DEPCONST*
@@ -40,6 +40,9 @@
 !!    --------
 !!
 ! open the general purpose ASCII input file
+USE MODI_OPEN_NAMELIST
+USE MODI_CLOSE_NAMELIST
+
 USE MODI_CH_OPEN_INPUTB
 USE MODD_CH_SURF
 !!
@@ -57,7 +60,8 @@ IMPLICIT NONE
 !
 !*      0.1    declarations of arguments
 !
-INTEGER,                  INTENT(IN)  :: KCH      ! chemistry input namelist logical unit
+CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! Program name
+CHARACTER(LEN=28), INTENT(IN)  :: HCHEM_SURF_FILE ! ascii file for chemistry aggregation
 INTEGER,                  INTENT(IN)  :: KLUOUT   ! output listing channel
  CHARACTER(LEN=*), DIMENSION(:),  INTENT(IN)  :: HSV      ! name of chemical species
 !
@@ -83,6 +87,8 @@ INTEGER :: IHENRY         ! number of chemical Henry constant to be read
 REAL             , DIMENSION(:,:), ALLOCATABLE :: ZHENRYVAL
                           !chemical Henry constant value
 !
+INTEGER           :: ICH      ! unit of input chemical file
+
 INTEGER :: JI, JN, JNREAL ! loop control variables
 INTEGER :: INACT          ! array pointer
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -97,20 +103,21 @@ IF(.NOT. ALLOCATED(XSREALMASSMOLVAL)) ALLOCATE( XSREALMASSMOLVAL(SIZE(HSV,1)) )
 IF(.NOT. ALLOCATED(XSREALREACTVAL)  ) ALLOCATE( XSREALREACTVAL(SIZE(HSV,1)) )
 IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
 !
+CALL OPEN_NAMELIST(HPROGRAM,ICH,HFILE=HCHEM_SURF_FILE)
 !
 !*       2.  read chemical molecular diffusivity MASS_MOL
 !
 ! open input file
   WRITE(KLUOUT,*) &
        "CH_INIT_CONST: reading  molar mass"  
-  CALL CH_OPEN_INPUTB("MASS_MOL", KCH, KLUOUT)
+  CALL CH_OPEN_INPUTB("MASS_MOL", ICH, KLUOUT)
 !
 ! read number of molecular diffusivity IMASS
-  READ(KCH, *) IMASS
+  READ(ICH, *) IMASS
   WRITE(KLUOUT,*) "number of molecular diffusivity: ", IMASS
 !
 ! read data input format
-  READ(KCH,"(A)") YFORMAT
+  READ(ICH,"(A)") YFORMAT
   WRITE(KLUOUT,*) "input format is: ", YFORMAT
 !
 ! allocate fields
@@ -119,7 +126,7 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
 !
 ! read molecular diffusivity
   DO JI = 1, IMASS
-    READ(KCH,YFORMAT) YMASSMOLNAME(JI), ZMASSMOLVAL(JI)
+    READ(ICH,YFORMAT) YMASSMOLNAME(JI), ZMASSMOLVAL(JI)
     WRITE(KLUOUT,YFORMAT) YMASSMOLNAME(JI), ZMASSMOLVAL(JI)
   END DO
 !
@@ -150,14 +157,14 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
 ! open input file
    WRITE(KLUOUT,*) &
        "CH_INIT_CONST: reading  reactivity factor "  
-  CALL CH_OPEN_INPUTB("REA_FACT", KCH, KLUOUT)
+  CALL CH_OPEN_INPUTB("REA_FACT", ICH, KLUOUT)
 !
 ! read number of molecular diffusivity IREACT
-  READ(KCH, *) IREACT
+  READ(ICH, *) IREACT
   WRITE(KLUOUT,*) "number of reactivity factor : ", IREACT
 !
 ! read data input format
-  READ(KCH,"(A)") YFORMAT
+  READ(ICH,"(A)") YFORMAT
   WRITE(KLUOUT,*) "input format is: ", YFORMAT
 !
 ! allocate fields
@@ -165,7 +172,7 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
   ALLOCATE(ZREACTVAL(IREACT))
 ! read reactivity factor 
   DO JI = 1, IREACT
-    READ(KCH,YFORMAT) YREACTNAME(JI), ZREACTVAL(JI)
+    READ(ICH,YFORMAT) YREACTNAME(JI), ZREACTVAL(JI)
     WRITE(KLUOUT,YFORMAT) YREACTNAME(JI), ZREACTVAL(JI)
   END DO
 !
@@ -197,14 +204,14 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
   WRITE(KLUOUT,*) &
        "CH_INIT_CONST: reading effective Henry constant", &
        " and its temperature correction "  
-  CALL CH_OPEN_INPUTB("HENRY_SP", KCH, KLUOUT)
+  CALL CH_OPEN_INPUTB("HENRY_SP", ICH, KLUOUT)
 !
 ! read number of molecular diffusivity IHENRY
-  READ(KCH, *) IHENRY
+  READ(ICH, *) IHENRY
   WRITE(KLUOUT,*) "number of reactivity factor : ", IHENRY
 !
 ! read data input format
-  READ(KCH,"(A)") YFORMAT
+  READ(ICH,"(A)") YFORMAT
   WRITE(KLUOUT,*) "input format is: ", YFORMAT
 !
 ! allocate fields
@@ -213,7 +220,7 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
 !
 ! read reactivity factor 
   DO JNREAL = 1, IHENRY
-    READ(KCH,YFORMAT) YHENRYNAME(JNREAL), ZHENRYVAL(JNREAL,1),&
+    READ(ICH,YFORMAT) YHENRYNAME(JNREAL), ZHENRYVAL(JNREAL,1),&
                              ZHENRYVAL(JNREAL,2)  
     WRITE(KLUOUT,YFORMAT) YHENRYNAME(JNREAL), ZHENRYVAL(JNREAL,1),&
                              ZHENRYVAL(JNREAL,2)  
@@ -241,6 +248,9 @@ IF(.NOT. ALLOCATED(XSREALHENRYVAL)  ) ALLOCATE( XSREALHENRYVAL(SIZE(HSV,1),2) )
                       XSREALHENRYVAL(JNREAL,1),&
                       XSREALHENRYVAL(JNREAL,2)  
   END DO
+
+CALL CLOSE_NAMELIST(HPROGRAM,ICH)
+
 IF (LHOOK) CALL DR_HOOK('CH_INIT_DEPCONST',1,ZHOOK_HANDLE)
 !
 END SUBROUTINE CH_INIT_DEPCONST
diff --git a/src/SURFEX/ch_init_emissionn.F90 b/src/SURFEX/ch_init_emissionn.F90
index 556fa52a97d8c56a92506cf6d8f31820ff71109c..6d4ee9f3f775f9be0a66d036b89d809e4c736312 100644
--- a/src/SURFEX/ch_init_emissionn.F90
+++ b/src/SURFEX/ch_init_emissionn.F90
@@ -3,7 +3,7 @@
 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !SURFEX_LIC for details. version 1.
 !     #########
-      SUBROUTINE CH_INIT_EMISSION_n(HPROGRAM,KLU,KCH,PRHOA)
+      SUBROUTINE CH_INIT_EMISSION_n(HPROGRAM,KLU,HINIT,PRHOA,HCHEM_SURF_FILE)
 !     #######################################
 !
 !!****  *CH_INIT_EMIISION_n* - routine to initialize chemical emissions data structure
@@ -26,6 +26,7 @@
 !!      Original        08/03/2001
 !!      D.Gazen  01/12/03  change emissions handling for surf. externalization
 !!      P.Tulet  01/01/04  introduction of rhodref for externalization
+!!      M.Leriche & V. Masson 05/16 bug in write emis fields for nest
 !-----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -36,6 +37,9 @@ USE MODI_GET_LUOUT
 USE MODI_BUILD_EMISSTAB_n
 USE MODI_BUILD_PRONOSLIST_n
 USE MODI_READ_SURF
+USE MODI_OPEN_NAMELIST
+USE MODI_CLOSE_NAMELIST
+USE MODI_READ_SURF_FIELD2D
 !
 !
 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
@@ -49,8 +53,13 @@ IMPLICIT NONE
 !
  CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! Program name
 INTEGER,           INTENT(IN)  :: KLU      ! number of points
-INTEGER,           INTENT(IN)  :: KCH      ! logical unit of input chemistry file
+CHARACTER(LEN=3),  INTENT(IN)  :: HINIT    ! Flag to know if one initializes:
+!                                          ! 'ALL' : all variables for a run
+!                                          ! 'PRE' : only variables to build 
+!                                          !         an initial file
+
 REAL, DIMENSION(:),INTENT(IN)  :: PRHOA    ! air density
+CHARACTER(LEN=28), INTENT(IN)  :: HCHEM_SURF_FILE ! ascii file for chemistry aggregation
 !
 !*       0.2   declarations of local variables
 !
@@ -69,7 +78,9 @@ INTEGER,DIMENSION(:),ALLOCATABLE  :: IOFFNDX ! index array of offline emission s
 INTEGER                           :: INBTS   ! number of emission times for a species
 INTEGER                           :: INBOFF  ! Number of offline emissions
 INTEGER                           :: IVERB   ! verbose level
- CHARACTER(LEN=3)                  :: YSURF   ! surface type
+INTEGER                           :: ICH      ! logical unit of input chemistry file
+CHARACTER(LEN=3)                  :: YSURF   ! surface type
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK2D ! work array to read emission fields
 !
 INTEGER           :: IVERSION       ! version of surfex file being read
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -82,7 +93,7 @@ WRITE(ILUOUT,*) '------ Beginning of CH_INIT_EMISSION ------'
 YRECFM='VERSION'
  CALL READ_SURF(HPROGRAM,YRECFM,IVERSION,IRESP)
 !
-!*      2.     Chemical Emission fields
+!*      1.     Chemical Emission fields
 !              ------------------------
 !
 ! Read the total number of emission files 
@@ -114,6 +125,11 @@ END IF
 
 IF (.NOT. ASSOCIATED(CEMIS_AREA))   ALLOCATE(CEMIS_AREA(NEMISPEC_NBR))
 IF (.NOT. ASSOCIATED(NEMIS_TIME))   ALLOCATE(NEMIS_TIME(NEMIS_NBR))
+
+IF (HINIT/='ALL') THEN
+  ALLOCATE(XEMIS_FIELDS(KLU,NEMIS_NBR))
+  ALLOCATE(CEMIS_COMMENT(NEMIS_NBR))
+END IF
 !
 ALLOCATE(ITIMES(NEMIS_NBR))
 ALLOCATE(INBTIMES(NEMISPEC_NBR))
@@ -168,6 +184,16 @@ DO JSPEC = 1,NEMISPEC_NBR ! Loop on the number of species
   CEMIS_NAME(JSPEC) = YSPEC_NAME
   CEMIS_AREA(JSPEC) = YSURF
 ! 
+!*      2.     Simple reading of emission fields
+
+  IF (HINIT /= "ALL") THEN
+    YRECFM='E_'//TRIM(ADJUSTL(YSPEC_NAME))
+    ALLOCATE(ZWORK2D(KLU,INBTS))
+    CALL READ_SURF_FIELD2D(HPROGRAM,ZWORK2D(:,:),YRECFM,YCOMMENT)
+    XEMIS_FIELDS(:,IIND1:IIND2) = ZWORK2D(:,:)
+    CEMIS_COMMENT(IIND1:IIND2) = YCOMMENT
+    DEALLOCATE(ZWORK2D)
+  END IF
 END DO
 !
 WRITE(ILUOUT,*) '---- Nunmer of OFFLINE species = ',INBOFF
@@ -176,22 +202,27 @@ WRITE(ILUOUT,*) 'IOFFNDX=',IOFFNDX
 
 IVERB=6
 
-IF (INBOFF > 0) THEN
-  ALLOCATE(TSEMISS(INBOFF))
-  ALLOCATE(YEMIS_NAME(INBOFF))
+!*      3.     Conversion and aggregation
 
-  CALL BUILD_EMISSTAB_n(HPROGRAM,KCH,CEMIS_NAME,INBTIMES,NEMIS_TIME,&
+IF (HINIT == "ALL") THEN
+  IF (INBOFF > 0) THEN
+    CALL OPEN_NAMELIST(HPROGRAM,ICH,HFILE=HCHEM_SURF_FILE)
+    ALLOCATE(TSEMISS(INBOFF))
+    ALLOCATE(YEMIS_NAME(INBOFF))
+
+    CALL BUILD_EMISSTAB_n(HPROGRAM,ICH,CEMIS_NAME,INBTIMES,NEMIS_TIME,&
          IOFFNDX,TSEMISS,KLU,ILUOUT,IVERB,PRHOA)  
-  DO JSPEC = 1,INBOFF ! Loop on the number of species
-    YEMIS_NAME(JSPEC) = TSEMISS(JSPEC)%CNAME(1:12)
-  END DO
-  CALL BUILD_PRONOSLIST_n(SIZE(TSEMISS),YEMIS_NAME,TSPRONOSLIST,KCH,ILUOUT,IVERB)
-  DEALLOCATE(YEMIS_NAME)
-ELSE
-  ALLOCATE(TSEMISS(0))
-  NULLIFY(TSPRONOSLIST)
+    DO JSPEC = 1,INBOFF ! Loop on the number of species
+      YEMIS_NAME(JSPEC) = TSEMISS(JSPEC)%CNAME(1:12)
+    END DO
+    CALL BUILD_PRONOSLIST_n(SIZE(TSEMISS),YEMIS_NAME,TSPRONOSLIST,ICH,ILUOUT,IVERB)
+    DEALLOCATE(YEMIS_NAME)
+    CALL CLOSE_NAMELIST(HPROGRAM,ICH)
+  ELSE
+    ALLOCATE(TSEMISS(0))
+    NULLIFY(TSPRONOSLIST)
+  END IF
 END IF
-
 DEALLOCATE(ITIMES,INBTIMES,IOFFNDX)
 WRITE(ILUOUT,*) '------ Leaving CH_INIT_EMISSION ------'
 IF (LHOOK) CALL DR_HOOK('CH_INIT_EMISSION_N',1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/ch_init_snapn.F90 b/src/SURFEX/ch_init_snapn.F90
index ca708798c3a5c0d70a31df1e9a7893f112231a24..7cc563d9c3ffcf87fe40992051ec4e137268d5a5 100644
--- a/src/SURFEX/ch_init_snapn.F90
+++ b/src/SURFEX/ch_init_snapn.F90
@@ -3,7 +3,7 @@
 !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !SURFEX_LIC for details. version 1.
 !     #########
-      SUBROUTINE CH_INIT_SNAP_n(HPROGRAM,KLU,HINIT,KCH,PRHOA)
+      SUBROUTINE CH_INIT_SNAP_n(HPROGRAM,KLU,HINIT,PRHOA,HCHEM_SURF_FILE)
 !     #######################################
 !
 !!****  *CH_INIT_EMIISION_TEMP_n* - routine to initialize chemical emissions data structure
@@ -25,6 +25,7 @@
 !!    -------------
 !!      Original        11/2011
 !!        M.Moge    01/2016  using READ_SURF_FIELD2D for 2D surfex fields reads
+!!      M.Leriche & V. Masson 05/16 move open namelist for reading ascii chemi.file
 !!-----------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -38,6 +39,8 @@ USE MODI_ABOR1_SFX
 USE MODI_CH_CONVERSION_FACTOR
 USE MODI_BUILD_PRONOSLIST_n
 USE MODI_CH_OPEN_INPUTB
+USE MODI_OPEN_NAMELIST
+USE MODI_CLOSE_NAMELIST
 !
 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
 USE PARKIND1  ,ONLY : JPRB
@@ -53,8 +56,8 @@ INTEGER,           INTENT(IN)  :: KLU      ! number of points
 !                                          ! 'ALL' : all variables for a run
 !                                          ! 'PRE' : only variables to build 
 !                                          !         an initial file
-INTEGER,           INTENT(IN)  :: KCH      ! logical unit of input chemistry file
 REAL, DIMENSION(:),INTENT(IN)  :: PRHOA    ! air density
+CHARACTER(LEN=28), INTENT(IN)  :: HCHEM_SURF_FILE ! ascii file for chemistry aggregation
 !
 !*       0.2   declarations of local variables
 !
@@ -69,6 +72,7 @@ INTEGER             :: JSNAP                 ! Loop index for SNAP categories
 !
 INTEGER             :: IVERSION       ! version of surfex file being read
 INTEGER             :: IBUG           ! version of SURFEX bugfix
+INTEGER           :: ICH      ! unit of input chemical file
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !-------------------------------------------------------------------------------
 IF (LHOOK) CALL DR_HOOK('CH_INIT_SNAP_N',0,ZHOOK_HANDLE)
@@ -133,7 +137,7 @@ DO JSPEC = 1,NEMIS_NBR ! Loop on the number of species
 ! 
 ! Read  the potential emission of species for each snap
   DO JSNAP=1,NEMIS_SNAP
-    WRITE(YRECFM,'("SNAP",I2.2,"_",A3)') JSNAP,CEMIS_NAME(JSPEC)
+    WRITE(YRECFM,'("SN",I2.2,"_",A7)') JSNAP,CEMIS_NAME(JSPEC)
     CALL READ_SURF(HPROGRAM,YRECFM,XEMIS_FIELDS_SNAP(:,JSNAP,JSPEC),IRESP,YCOMMENT)
   END DO
 !
@@ -150,10 +154,11 @@ END DO
 !              -----------------
 !
 IF (HINIT=='ALL') THEN
-  CALL CH_OPEN_INPUTB("EMISUNIT", KCH, ILUOUT)
+  CALL OPEN_NAMELIST(HPROGRAM,ICH,HFILE=HCHEM_SURF_FILE)
+  CALL CH_OPEN_INPUTB("EMISUNIT", ICH, ILUOUT)
 !
 ! read unit identifier
-  READ(KCH,'(A3)') CCONVERSION
+  READ(ICH,'(A3)') CCONVERSION
 !
   ALLOCATE (XCONVERSION(KLU))
 ! determine the conversion factor
@@ -162,7 +167,8 @@ IF (HINIT=='ALL') THEN
 !*      4.     List of emissions to be aggregated into atm. chemical species
 !              -------------------------------------------------------------
 !
-  CALL BUILD_PRONOSLIST_n(NEMIS_NBR,CEMIS_NAME,TSPRONOSLIST,KCH,ILUOUT,6)
+  CALL BUILD_PRONOSLIST_n(NEMIS_NBR,CEMIS_NAME,TSPRONOSLIST,ICH,ILUOUT,6)
+  CALL CLOSE_NAMELIST(HPROGRAM,ICH)
 !
 !-------------------------------------------------------------------------------
 END IF
diff --git a/src/SURFEX/init_surf_atmn.F90 b/src/SURFEX/init_surf_atmn.F90
index 71a3d5bac34b8d33bc288e99a1273f57bbcc2d9f..a999d362a20235bdd7b9e5856e98b2d391a52a55 100644
--- a/src/SURFEX/init_surf_atmn.F90
+++ b/src/SURFEX/init_surf_atmn.F90
@@ -50,6 +50,7 @@ SUBROUTINE INIT_SURF_ATM_n(HPROGRAM,HINIT, OLAND_USE,                   &
 !!     (B. Decharme)   2013   Read grid only once in AROME case
 !!     (G. Tanguy)     2013   Add IF(ALLOCATED(NMASK_FULL))  before deallocate
 !!     (J.Durand)      2014   add activation of chemical deposition if LCH_EMIS=F
+!!      M.Leriche & V. Masson 05/16 bug in write emis fields for nest
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -70,7 +71,7 @@ USE MODD_SURF_ATM_SSO_n, ONLY : CROUGH, XAOSIP, XAOSIM, XAOSJP, XAOSJM, &
                                 XHO2IP, XHO2IM, XHO2JP, XHO2JM,         &
                                 XZ0EFFIP, XZ0EFFIM, XZ0EFFJP, XZ0EFFJM, &
                                 XZ0REL, XZ0EFFJPDIR, XFRACZ0, XCOEFBE
-USE MODD_CH_SURF_n,      ONLY : CCH_NAMES, LCH_EMIS, LRW_CH_EMIS, &
+USE MODD_CH_SURF_n,      ONLY : CCH_NAMES, LCH_EMIS,              &
                                 LCH_SURF_EMIS, CCHEM_SURF_FILE, CAER_NAMES,&
                                 CCH_EMIS
 USE MODD_SV_n,           ONLY : NBEQ, CSV, NSV_CHSBEG, NSV_CHSEND, &
@@ -128,8 +129,6 @@ USE MODI_INIT_CHEMICAL_n
 USE MODI_CH_INIT_DEPCONST
 USE MODI_CH_INIT_EMISSION_n
 USE MODI_CH_INIT_SNAP_n
-USE MODI_OPEN_NAMELIST
-USE MODI_CLOSE_NAMELIST
 USE MODI_READ_PRECIP_n
 USE MODI_ABOR1_SFX
 USE MODI_ALLOC_DIAG_SURF_ATM_n
@@ -196,7 +195,6 @@ INTEGER           :: ISWB     ! number of shortwave bands
 INTEGER           :: JTILE    ! loop counter on tiles
 INTEGER           :: IRESP    ! error return code
 INTEGER           :: ILUOUT   ! unit of output listing file
-INTEGER           :: ICH      ! unit of input chemical file
 INTEGER           :: IVERSION, IBUGFIX       ! surface version
 !
 REAL, DIMENSION(:,:), ALLOCATABLE                       :: ZFRAC_TILE     ! fraction of each surface type
@@ -405,21 +403,11 @@ IF (LCH_EMIS) THEN
     CALL READ_SURF(HPROGRAM,'CH_EMIS_OPT',CCH_EMIS,IRESP)
   END IF
   !
-  IF (CCH_EMIS=='AGGR') LRW_CH_EMIS = .TRUE.
-  !
-  IF (NBEQ > 0) THEN
-    !
-    CALL OPEN_NAMELIST(HPROGRAM,ICH,HFILE=CCHEM_SURF_FILE)
-    !
-    IF (LCH_SURF_EMIS) THEN
-      IF (CCH_EMIS=='AGGR') THEN
-        CALL CH_INIT_EMISSION_n(HPROGRAM,NSIZE_FULL,ICH,PRHOA) 
-      ELSE
-        CALL CH_INIT_SNAP_n(HPROGRAM,NSIZE_FULL,HINIT,ICH,PRHOA)
-      END IF
-    ENDIF
-    CALL CLOSE_NAMELIST(HPROGRAM,ICH)
-  ENDIF
+  IF (CCH_EMIS=='AGGR') THEN
+     CALL CH_INIT_EMISSION_n(HPROGRAM,NSIZE_FULL,HINIT,PRHOA,CCHEM_SURF_FILE) 
+  ELSE
+     CALL CH_INIT_SNAP_n(HPROGRAM,NSIZE_FULL,HINIT,PRHOA,CCHEM_SURF_FILE)
+  END IF
   !
 END IF
     !
@@ -427,11 +415,8 @@ END IF
     !    
 !
 IF (NBEQ .GT. 0) THEN
- CALL OPEN_NAMELIST(HPROGRAM,ICH,HFILE=CCHEM_SURF_FILE)
 
- IF (HINIT=='ALL') CALL CH_INIT_DEPCONST(ICH,ILUOUT,CSV(NSV_CHSBEG:NSV_CHSEND))
-!
- CALL CLOSE_NAMELIST(HPROGRAM,ICH)
+ IF (HINIT=='ALL') CALL CH_INIT_DEPCONST(HPROGRAM, CCHEM_SURF_FILE,ILUOUT,CSV(NSV_CHSBEG:NSV_CHSEND))
 !
 END IF
 !
diff --git a/src/SURFEX/writesurf_snapn.F90 b/src/SURFEX/writesurf_snapn.F90
index 252a9dc4dcc4e3b31a39048d13439cfae47947a5..1cc48df2b17feb78f9019438f5be0d7cdcbc87b9 100644
--- a/src/SURFEX/writesurf_snapn.F90
+++ b/src/SURFEX/writesurf_snapn.F90
@@ -78,7 +78,7 @@ DO JSPEC=1,NEMIS_NBR
   CALL WRITE_SURF_FIELD2D(HPROGRAM,XSNAP_HOURLY(:,:,JSPEC),YRECFM,YCOMMENT,YCOMMENTUNIT,HDIR='-')
 ! Writes the potential emission of species for each snap
   DO JSNAP=1,NEMIS_SNAP
-    WRITE(YRECFM,'("SNAP",I2.2,"_",A3)') JSNAP,CEMIS_NAME(JSPEC)
+    WRITE(YRECFM,'("SN",I2.2,"_",A7)') JSNAP,CEMIS_NAME(JSPEC)
     CALL WRITE_SURF(HPROGRAM,YRECFM,XEMIS_FIELDS_SNAP(:,JSNAP,JSPEC),IRESP,YCOMMENT)
   END DO
 !