Index: trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/aerosol_mod.F90	(revision 3193)
@@ -38,9 +38,7 @@
   use radinc_h, only: naerkind
   use tracer_h, only: n_rgcs, nqtot, is_rgcs
-  use callkeys_mod, only: aeron2, dusttau, aeroh2so4, &
+  use callkeys_mod, only: aeron2, dusttau, &
                           aeroback2lay, aeronh3, nlayaero, aeronlay, &
-                          aeroaurora, aerogeneric, &
-                          aerovenus1, aerovenus2, aerovenus2p, &
-                          aerovenus3, aerovenusUV
+                          aeroaurora, aerogeneric
 
   IMPLICIT NONE
@@ -89,10 +87,5 @@
   if (is_master) write(*,*) '--- Dust aerosol = ', iaero_dust
 
-  if (aeroh2so4) then
-     ia=ia+1
-     iaero_h2so4=ia
-  endif
-  if (is_master) write(*,*) '--- H2SO4 aerosol = ', iaero_h2so4
-      
+     
   if (aeroback2lay) then
      ia=ia+1
Index: trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90	(revision 3193)
@@ -37,5 +37,5 @@
                               diagdtau,kastprof,strictboundcorrk,specOLR, &
                               tplanckmin,tplanckmax,global1d, &
-                              generic_condensation, aerovenus
+                              generic_condensation
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -808,17 +808,9 @@
       plevrad(1) = 0.
 !      plevrad(2) = 0.   !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. 
-      if (aerovenus) then
-!!  GG19 modified below after SL routines
-        plevrad(2) = 0.
-      endif
-
+      
       tlevrad(1) = tlevrad(2)
       tlevrad(2*nlayer+1)=tsurf(ig)
       
       pmid(1) = pplay(ig,nlayer)/scalep   
-      if (aerovenus) then
-!! GG19 modified below after SL routines
-        pmid(1) = max(pgasmin,0.0001*plevrad(3))
-      endif
       pmid(2) =  pmid(1)
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90	(revision 3193)
@@ -6,6 +6,6 @@
       logical,save :: calladj,calltherm,n2cond,callsoil
 !$OMP THREADPRIVATE(calladj,calltherm,n2cond,callsoil)
-      logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
-!$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
+      logical,save :: season,diurnal,lwrite
+!$OMP THREADPRIVATE(season,diurnal,lwrite)
       logical,save :: callgasvis,continuum,graybody
 !$OMP THREADPRIVATE(callgasvis,continuum,graybody)
@@ -37,18 +37,9 @@
       logical,save :: sedimentation
       logical,save :: generic_condensation
-      logical,save :: generic_rain
 !$OMP THREADPRIVATE(varactive,varfixed,sedimentation,generic_condensation,generic_rain)
-      logical,save :: water ,watercond, waterrain, moistadjustment
-!$OMP THREADPRIVATE(water, watercond, waterrain, moistadjustment)
       logical,save :: aeron2, aeroh2o, aeroh2so4, aeroback2lay
 !$OMP THREADPRIVATE(aeron2, aeroh2o, aeroh2so4, aeroback2lay)
       logical,save :: aeronh3, aeronlay, aeroaurora
 !$OMP THREADPRIVATE(aeronh3,aeronlay,aeroaurora)
-
-      logical,save :: aerovenus ! master flag for "Venus-like" aerosol additions
-!$OMP THREADPRIVATE(aerovenus)
-      ! detailed sub-options when with "Venus-like" aerosol additions
-      logical,save :: aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV
-!$OMP THREADPRIVATE(aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV)
 
       logical,save :: aerofixn2
@@ -63,8 +54,43 @@
       logical,save :: jonline
       logical,save :: depos
-      logical,save :: haze
 !$OMP THREADPRIVATE(photoheat,jonline,depos)
-      logical,save :: calllott_nonoro
-!$OMP THREADPRIVATE(calllott_nonoro)
+
+!! Pluto-specific variables
+      logical,save :: haze_radproffix,haze_proffix 
+!$OMP THREADPRIVATE(haze_radproffix,haze_proffix)      
+      logical,save :: haze,fasthaze,changeti,changetid,aerohaze         
+!$OMP THREADPRIVATE(haze,fasthaze,changeti,changetid,aerohaze)      
+      logical,save :: fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel
+!$OMP THREADPRIVATE(fast,metcloud,monoxcloud,glaflow,triton,paleo,nlte,strobel      
+      logical,save :: kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice 
+!$OMP THREADPRIVATE(kbo,nbsub,cooling,mode_n2,thres_n2ice,thres_coice)      
+      logical,save :: deltab,nb_monomer                                 
+!$OMP THREADPRIVATE(deltab,nb_monomer)      
+      logical,save :: metlateq,tholateq,metls1,metls2                   
+!$OMP THREADPRIVATE(metlateq,tholateq,metls1,metls2)      
+      logical,save :: mode_ch4,mode_tholins                             
+!$OMP THREADPRIVATE(mode_ch4,mode_tholins)      
+      logical,save :: source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4      
+!$OMP THREADPRIVATE(source_haze,ch4lag,latlag,vmrlag,kfix,hazeconservch4)      
+      logical,save :: fracsource,latsource,lonsource                    
+!$OMP THREADPRIVATE(fracsource,latsource,lonsource)      
+      logical,save :: spelon1,spelon2,spelat1,spelat2,specalb           
+!$OMP THREADPRIVATE(spelon1,spelon2,spelat1,spelat2,specalb)      
+      logical,save :: assymflux,mode_hs,tsurfmax,albmin_ch4
+!$OMP THREADPRIVATE(assymflux,mode_hs,tsurfmax,albmin_ch4)      
+      logical,save :: feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats   
+!$OMP THREADPRIVATE(feedback_met,thres_ch4ice,fdch4_latn,fdch4_lats)      
+      logical,save :: globmean1d,kmix_proffix,rad_haze          
+!$OMP THREADPRIVATE(globmean1d,kmix_proffix,rad_haze)      
+      logical,save :: fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb
+!$OMP THREADPRIVATE(fdch4_lone,fdch4_lonw,fdch4_maxalb,fdch4_depalb,fdch4_finalb      
+      logical,save :: tholatn,tholats,tholone,tholonw                   
+!$OMP THREADPRIVATE(tholatn,tholats,tholone,tholonw)      
+      logical,save :: fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff   
+!$OMP THREADPRIVATE(fdch4_ampl,fdch4_maxice,condmetsurf,condcosurf,vertdiff)      
+      logical,save :: convergeps,conservn2,patchflux,condensn2,no_n2frost
+!$OMP THREADPRIVATE(convergeps,conservn2,patchflux,condensn2,no_n2frost      
+
+
       logical,save :: global1d
       real,save    :: szangle
Index: trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/callsedim.F	(revision 3193)
@@ -8,5 +8,4 @@
       USE tracer_h, only : igcm_n2_ice,radius,rho_q
       use comcstfi_mod, only: g
-      use callkeys_mod, only : water
 
       IMPLICIT NONE
@@ -75,9 +74,4 @@
         firstcall=.false.
         ! add some tests on presence of required tracers/aerosols:
-    !     if (water.and.(igcm_h2o_ice.eq.0)) then
-    !         write(*,*) "callsedim error: water=.true.",
-    !  &                 " but igcm_h2o_ice=0"
-    !       stop
-    !     endif
         ! allocate "naerkind" size local arrays (which are also
         ! "saved" so that this is done only once in for all even if
@@ -120,13 +114,5 @@
 ! Sedimentation 
 !======================================================================
-! Water          
-!        if (water.and.(iaero_h2o.ne.0).and.(iq.eq.igcm_h2o_ice)) then
-!             ! compute radii for h2o_ice 
-!              call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
-!      &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
-!             ! call sedimentation for h2o_ice
-!              call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
-!      &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
-!      &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq,iq)
+! Water      !AF24: deleted
 
 ! ! General Case
Index: trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90	(revision 3193)
@@ -1,526 +1,875 @@
-      subroutine condense_n2(ngrid,nlayer,nq,ptimestep, &
-          pcapcal,pplay,pplev,ptsrf,pt,                  &
-          pdt,pdtsrf,pq,pdq,                             &
-          pqsurf,pdqsurfc,albedo,pemisurf,               &
-          albedo_bareground,albedo_n2_ice_SPECTV,       &
-          pdtc,pdtsrfc,pdpsrfc,pdqc)
-
-      use radinc_h, only : L_NSPECTV
-      use gases_h, only: gfrac, igas_n2
-      use radii_mod, only : n2_reffrad
-      use aerosol_mod, only : iaero_n2
-      USE surfdat_h, only: emisice, emissiv
-      USE geometry_mod, only: latitude ! in radians
-      USE tracer_h, only: noms, rho_n2
-      use comcstfi_mod, only: g, r, cpp
-
-      implicit none
+subroutine condens_n2(klon,nlayer,nq,ptimestep, &
+      pcapcal,pplay,pplev,ptsrf,pt, &
+      pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
+      picen2,psolaralb,pemisurf, &
+      pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
+      pdqc,pdicen2)
+
+  use radinc_h, only : naerkind
+  use comgeomfi_h
+  use comcstfi_mod, only: g, r, cpp
+  USE surfdat_h, only: phisfi
+  USE tracer_h, only: noms, igcm_n2
+  USE callkeys_mod, only: fast,ch4lag,latlag,lw_n2,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
+  USE dimphy, only : klon,klev
+
+
+  implicit none
 
 !==================================================================
 !     Purpose
 !     -------
-!     Condense and/or sublime N2 ice on the ground and in the atmosphere, and sediment the ice.
+!     Condense and/or sublime N2 ice on the ground and in the
+!     atmosphere, and sediment the ice.
 !  
 !     Inputs
-!     ------
-!     ngrid                 Number of vertical columns.
-!     nlayer                Number of vertical layers.
-!     nq                    Number of tracers.
-!     ptimestep             Duration of the physical timestep (s).
-!     pplay(ngrid,nlayer)   Pressure layers (Pa).
-!     pplev(ngrid,nlayer+1) Pressure levels (Pa).
-!     pt(ngrid,nlayer)      Atmospheric Temperatures (K).
-!     ptsrf(ngrid)          Surface temperatures (K).
-!     pq(ngrid,nlayer,nq)   Atmospheric tracers mixing ratios (kg/kg of air).
-!     pqsurf(ngrid,nq)      Surface tracers (kg/m2).
+!     ------ 
+!     klon                 Number of vertical columns
+!     nlayer                Number of layers
+!     pplay(klon,nlayer)   Pressure layers
+!     pplev(klon,nlayer+1) Pressure levels 
+!     pt(klon,nlayer)      Temperature (in K)
+!     ptsrf(klon)          Surface temperature
 !     
-!     pdt(ngrid,nlayer)     Time derivative before condensation/sublimation of pt.
-!     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf.
-!     pdq(ngrid,nlayer,nq)  Time derivative before condensation/sublimation of
-!
-!     albedo_bareground(ngrid)           Albedo of the bare ground.
-!     albedo_n2_ice_SPECTV(L_NSPECTV)   Spectral albedo of N2 ice.
+!     pdt(klon,nlayermx)   Time derivative before condensation/sublimation of pt
+!     pdtsrf(klon)         Time derivative before condensation/sublimation of ptsrf
+!     picen2(klon)         n2 ice at the surface (kg/m2)
 !     
 !     Outputs
 !     -------
-!     pdpsrfc(ngrid)          \       Contribution of condensation/sublimation 
-!     pdtc(ngrid,nlayer)       \            to the time derivatives of 
-!     pdtsrfc(ngrid)           /     Surface Pressure, Atmospheric Temperatures,
-!     pdqsurfc(ngrid)         /         Surface Temperatures, Surface Tracers,
-!     pdqc(ngrid,nlayer,nq)  /                and Atmospheric Tracers.*
-!
-!     pemisurf(ngrid)              Emissivity of the surface.
+!     pdpsrf(klon)         \  Contribution of condensation/sublimation
+!     pdtc(klon,nlayermx)  /  to the time derivatives of Ps, pt, and ptsrf
+!     pdtsrfc(klon)       /
+!     pdicen2(klon)         Tendancy n2 ice at the surface (kg/m2)
 !     
 !     Both
 !     ----
-!     albedo(ngrid,L_NSPECTV)      Spectral albedo of the surface.
+!     psolaralb(klon)      Albedo at the surface
+!     pemisurf(klon)       Emissivity of the surface
 !
 !     Authors
 !     ------- 
-!     Francois Forget (1996)
+!     Francois Forget (1996,2013)
 !     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
-!     Includes simplifed nucleation by J. Leconte (2011)
+!     
 !     
 !==================================================================
 
-!--------------------------
-!        Arguments
-!--------------------------
-
-
-      INTEGER,INTENT(IN) :: ngrid
-      INTEGER,INTENT(IN) :: nlayer
-      INTEGER,INTENT(IN) :: nq
-      REAL,INTENT(IN) :: ptimestep 
-      REAL,INTENT(IN) :: pcapcal(ngrid)
-      REAL,INTENT(IN) :: pplay(ngrid,nlayer)
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
-      REAL,INTENT(IN) :: ptsrf(ngrid)
-      REAL,INTENT(IN) :: pt(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdt(ngrid,nlayer)
-      REAL,INTENT(IN) :: pdtsrf(ngrid)
-      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
-      REAL,INTENT(IN) :: pqsurf(ngrid,nq)
-      REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
-      REAL,INTENT(IN) :: albedo_bareground(ngrid)
-      REAL,INTENT(IN) :: albedo_n2_ice_SPECTV(L_NSPECTV)
-      REAL,INTENT(INOUT) :: albedo(ngrid,L_NSPECTV)
-      REAL,INTENT(OUT) :: pemisurf(ngrid)
-      REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
-      REAL,INTENT(OUT) :: pdtsrfc(ngrid)
-      REAL,INTENT(OUT) :: pdpsrfc(ngrid)
-      REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
-      REAL,INTENT(OUT) :: pdqsurfc(ngrid)
-
-!------------------------------
-!       Local variables
-!------------------------------
-
-      INTEGER l,ig,icap,ilay,iq,nw,igas,it
-
-      REAL reffrad(ngrid,nlayer) ! Radius (m) of the N2 ice particles.
-      REAL*8 zt(ngrid,nlayer)    ! Updated Atmospheric Temperatures (K).
-      REAL ztsrf(ngrid)          ! Updated Surface Temperatures (K).
-      REAL zq(ngrid,nlayer,nq)   ! Updated Atmospheric tracers mixing ratios (kg/kg of air).
-      REAL picen2(ngrid)        ! Updated Surface Tracer (kg/m2).
-      REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of N2.
-      REAL ztnuc (ngrid,nlayer)  ! Atmospheric Nucleation Temperatures.
-      REAL ztcondsol(ngrid)      ! Temperatures of condensation of N2 at the surface.
-      REAL zcondices(ngrid)      ! Condensation rate on the ground (kg/m2/s).
-      REAL zfallice(ngrid)       ! Flux of ice falling on the surface (kg/m2/s).
-      REAL Mfallice(ngrid)       ! Total amount of ice fallen to the ground during the timestep (kg/m2).
-      REAL wq(ngrid,nlayer+1)    ! Total amount of ice fallen to the ground during the timestep (kg/m2).
-      REAL subptimestep          ! Duration of the subtimestep (s) for the sedimentation.
-      Integer Ntime              ! Number of subtimesteps.
-      REAL masse (ngrid,nlayer)  ! Mass of atmospheric layers (kg/m2)
-      REAL w(ngrid,nlayer,nq)    !
-      REAL vstokes,reff          !
-      REAL ppn2                 !
-
-
-!------------------------------------------
-!         Saved local variables
-!------------------------------------------
-
-
-      REAL,SAVE :: latcond=2.5e5
-      REAL,SAVE :: ccond
-      REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
-!$OMP THREADPRIVATE(latcond,ccond,emisref)
-      LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-      INTEGER,SAVE :: i_n2ice=0      ! n2 ice
-!$OMP THREADPRIVATE(i_n2ice)
-      CHARACTER(LEN=20) :: tracername ! to temporarily store text
-
-
-!------------------------------------------------
-!       Initialization at the first call
-!------------------------------------------------
-
-
-      IF (firstcall) THEN
-
-!         ALLOCATE(emisref(ngrid))
-!         ! Find N2 ice tracer.
-!         do iq=1,nq
-!            tracername=noms(iq)
-!            if (tracername.eq."n2_ice") then
-!               i_n2ice=iq
-!            endif
-!         enddo
-         
-!         write(*,*) "condense_n2: i_n2ice=",i_n2ice       
-
-!         if((i_n2ice.lt.1))then
-!            print*,'In condens_cloud but no N2 ice tracer, exiting.'
-!            print*,'Still need generalisation to arbitrary species!'
-!            stop
-!         endif
-
-         ccond=cpp/(g*latcond)
-         print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
-
-         ! Prepare special treatment if gas is not pure N2
-         ! if (addn2) then
-         !    m_n2   = 44.01E-3 ! N2 molecular mass (kg/mol)   
-         !    m_non2 = 28.02E-3 ! N2 molecular mass (kg/mol)  
-         ! Compute A and B coefficient use to compute
-         ! mean molecular mass Mair defined by
-         ! 1/Mair = q(in2)/m_n2 + (1-q(in2))/m_non2
-         ! 1/Mair = A*q(in2) + B
-         ! A = (1/m_n2 - 1/m_non2)
-         ! B = 1/m_non2
-         ! endif
-
-         ! Minimum N2 mixing ratio below which mixing occurs with layer above : qn2min =0.75  
-
-           firstcall=.false.
-      ENDIF
-
-
-!------------------------------------------------
-!        Tendencies initially set to 0
-!------------------------------------------------
-
-
-      pdqc(1:ngrid,1:nlayer,1:nq)     = 0.
-      pdtc(1:ngrid,1:nlayer)          = 0.
-      zq(1:ngrid,1:nlayer,1:nq)       = 0.
-      zt(1:ngrid,1:nlayer)            = 0.
-      Mfallice(1:ngrid)               = 0.      
-      zfallice(1:ngrid)               = 0.
-      zcondices(1:ngrid)              = 0.
-      pdtsrfc(1:ngrid)                = 0.
-      pdpsrfc(1:ngrid)                = 0.
-      pdqsurfc(1:ngrid)               = 0.
-
-
-!----------------------------------
+!-----------------------------------------------------------------------
+!     Arguments
+
+  INTEGER klon, nlayer, nq
+
+  REAL ptimestep 
+  REAL pplay(klon,nlayer),pplev(klon,nlayer+1)
+  REAL pcapcal(klon)
+  REAL pt(klon,nlayer)
+  REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon)
+  REAL pphi(klon,nlayer)
+  REAL pdt(klon,nlayer),pdtsrf(klon),pdtc(klon,nlayer)
+  REAL pdtsrfc(klon),pdpsrf(klon)
+  REAL picen2(klon),psolaralb(klon),pemisurf(klon)
+ 
+
+  REAL pu(klon,nlayer) ,  pv(klon,nlayer)
+  REAL pdu(klon,nlayer) , pdv(klon,nlayer)
+  REAL pduc(klon,nlayer) , pdvc(klon,nlayer)
+  REAL pq(ngridmx,nlayer,nq),pdq(klon,nlayer,nq)
+  REAL pdqc(klon,nlayer,nq)
+
+!-----------------------------------------------------------------------
+!     Local variables
+
+  INTEGER l,ig,ilay,it,iq,ich4_gas
+
+  REAL*8 zt(ngridmx,nlayermx)
+  REAL tcond_n2
+  REAL pcond_n2
+  REAL glob_average2d         ! function
+  REAL zqn2(ngridmx,nlayermx) ! N2 MMR used to compute Tcond/zqn2
+  REAL ztcond (ngridmx,nlayermx)
+  REAL ztcondsol(ngridmx),zfallheat
+  REAL pdicen2(ngridmx)
+  REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
+  REAL zfallice(ngridmx,nlayermx+1)
+  REAL zmflux(nlayermx+1)
+  REAL zu(nlayermx),zv(nlayermx)
+  REAL zq(nlayermx,nqmx),zq1(nlayermx)
+  REAL ztsrf(ngridmx)
+  REAL ztc(nlayermx), ztm(nlayermx+1) 
+  REAL zum(nlayermx+1) , zvm(nlayermx+1)
+  REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
+  LOGICAL condsub(ngridmx)
+  REAL subptimestep
+  Integer Ntime
+  real masse (nlayermx),w(nlayermx+1)
+  real wq(ngridmx,nlayermx+1)
+  real vstokes,reff
+  real dWtotsn2
+  real condnconsn2(ngridmx)
+  real nconsMAXn2
+!     Special diagnostic variables
+  real tconda1(ngridmx,nlayermx)
+  real tconda2(ngridmx,nlayermx)
+  real zdtsig (ngridmx,nlayermx)
+  real zdtlatent (ngridmx,nlayermx)
+  real zdt (ngridmx,nlayermx)
+  REAL albediceF(ngridmx)
+  SAVE albediceF
+  INTEGER nsubtimestep,itsub    !number of subtimestep when calling vl1d
+  REAL subtimestep        !ptimestep/nsubtimestep
+  REAL dtmax              
+
+  REAL zplevhist(ngridmx)
+  REAL zplevnew(ngridmx)
+  REAL zplev(ngridmx)
+  REAL zpicen2(ngridmx)
+  REAL ztsrfhist(ngridmx)
+  REAL zdtsrf(ngridmx)
+  REAL globzplevnew
+
+  REAL vmrn2(ngridmx)
+  SAVE vmrn2
+  REAL stephan   
+  DATA stephan/5.67e-08/  ! Stephan Boltzman constant
+  SAVE stephan
+!-----------------------------------------------------------------------
+!     Saved local variables
+
+!      REAL latcond
+  REAL ccond
+  REAL cpice  ! for atm condensation
+  SAVE cpice
+!      SAVE latcond,ccond
+  SAVE ccond
+
+  LOGICAL firstcall
+  SAVE firstcall
+  REAL SSUM
+  EXTERNAL SSUM
+
+!      DATA latcond /2.5e5/
+!      DATA latcond /1.98e5/
+  DATA cpice /1300./
+  DATA firstcall/.true./
+
+  INTEGER,SAVE :: i_n2ice=0       ! n2 ice
+  CHARACTER(LEN=20) :: tracername ! to temporarily store text
+  logical olkin  ! option to prevent N2 ice effect in the south
+  DATA olkin/.false./
+  save olkin
+
+  CHARACTER(len=10) :: tname
+
+!-----------------------------------------------------------------------
+
+!     Initialisation
+  IF (firstcall) THEN
+     ccond=cpp/(g*lw_n2)
+     print*,'In condens_n2cloud: ccond=',ccond,' latcond=',lw_n2
+
+     ! calculate global mean surface pressure for the fast mode
+     IF (fast) THEN
+        DO ig=1,ngridmx
+           kp(ig) = exp(-phisfi(ig)/(r*38.))
+        ENDDO
+        p00=glob_average2d(kp) ! mean pres at ref level
+     ENDIF
+
+     vmrn2(:) = 1.
+     IF (ch4lag) then
+        DO ig=1,ngridmx
+           if (lati(ig)*180./pi.ge.latlag) then
+              vmrn2(ig) = vmrlag
+           endif
+        ENDDO
+     ENDIF   
+     IF (no_n2frost) then
+        DO ig=1,ngridmx
+           if (picen2(ig).eq.0.) then
+              vmrn2(ig) = 1.e-15
+           endif
+        ENDDO
+     ENDIF   
+     firstcall=.false.
+  ENDIF
+
+!-----------------------------------------------------------------------
+!     Calculation of n2 condensation / sublimation 
+
+!     Variables used:
+!     picen2(klon)            : amount of n2 ice on the ground     (kg/m2)
+!     zcondicea(klon,nlayermx): condensation rate in layer l       (kg/m2/s)
+!     zcondices(klon)         : condensation rate on the ground    (kg/m2/s)
+!     zfallice(klon,nlayermx) : amount of ice falling from layer l (kg/m2/s)
+!     zdtlatent(klon,nlayermx): dT/dt due to phase changes         (K/s)
+
+!     Tendencies initially set to 0
+  zcondices(1:klon) = 0.
+  pdtsrfc(1:klon) = 0.
+  pdpsrf(1:klon) = 0.
+  ztsrfhist(1:klon) = 0.
+  condsub(1:klon) = .false.
+  pdicen2(1:klon) = 0.
+  zfallheat=0
+  pdqc(1:klon,1:nlayer,1:nq)=0
+  pdtc(1:klon,1:nlayer)=0
+  pduc(1:klon,1:nlayer)=0
+  pdvc(1:klon,1:nlayer)=0
+  zfallice(1:klon,1:nlayer+1)=0
+  zcondicea(1:klon,1:nlayer)=0
+  zdtlatent(1:klon,1:nlayer)=0
+  zt(1:klon,1:nlayer)=0.
+
+!-----------------------------------------------------------------------
 !     Atmospheric condensation
-!----------------------------------
-
-
-!     Compute N2 Volume mixing ratio
-!     -------------------------------
-!      if (addn2) then
-!         DO l=1,nlayer
-!            DO ig=1,ngrid
-!              qn2=pq(ig,l,in2)+pdq(ig,l,in2)*ptimestep
-!              Mean air molecular mass = 1/(q(in2)/m_n2 + (1-q(in2))/m_non2)
-!              mmean=1/(A*qn2 +B)
-!              vmr_n2(ig,l) = qn2*mmean/m_n2 
-!            ENDDO
-!         ENDDO
-!      else
-!         DO l=1,nlayer
-!            DO ig=1,ngrid
-!              vmr_n2(ig,l)=0.5
-!            ENDDO
-!         ENDDO
-!      end if
-
-
-      ! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'.
-!      DO l=1,nlayer
-!         DO ig=1,ngrid
-!            ppn2=gfrac(igas_N2)*pplay(ig,l)
-!            call get_tcond_n2(ppn2,ztcond(ig,l))
-!            call get_tnuc_n2(ppn2,ztnuc(ig,l))
-!         ENDDO
-!      ENDDO
+
+!     Condensation / sublimation in the atmosphere
+!     --------------------------------------------
+!      (calcul of zcondicea , zfallice and pdtc)
+
+  zt(1:klon,1:nlayer)=pt(1:klon,1:nlayer)+ pdt(1:klon,1:nlayer)*ptimestep
+  if (igcm_n2.ne.0) then
+     zqn2(1:klon,1:nlayer) = 1. ! & temporaire
+!        zqn2(1:klon,1:nlayer)= pq(1:klon,1:nlayer,igcm_n2) + pdq(1:klon,1:nlayer,igcm_n2)*ptimestep
+  else
+     zqn2(1:klon,1:nlayer) = 1.
+  end if
+  
+  if (.not.fast) then
+!     Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
+    DO l=1,nlayer
+      DO ig=1,klon
+          ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l)) 
+      ENDDO
+    ENDDO
+
+    DO l=nlayer,1,-1
+      DO ig=1,klon
+       pdtc(ig,l)=0.  ! final tendancy temperature set to 0
+
+       IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
+           condsub(ig)=.true.  !condensation at level l
+           IF (zfallice(ig,l+1).gt.0) then 
+             zfallheat=zfallice(ig,l+1)*&  
+            (pphi(ig,l+1)-pphi(ig,l) +&
+           cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
+           ELSE
+               zfallheat=0.
+           ENDIF
+           zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
+           zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
+                         *ccond*zdtlatent(ig,l)- zfallheat
+!              Case when the ice from above sublimes entirely
+!              """""""""""""""""""""""""""""""""""""""""""""""
+           IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
+                .AND. (zfallice(ig,l+1).gt.0)) THEN
+
+              zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&     
+               (ccond*(pplev(ig,l)-pplev(ig,l+1)))
+              zcondicea(ig,l)= -zfallice(ig,l+1)
+           END IF
+
+           zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)
       
-      ! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf.
-!      DO ig=1,ngrid
-!         picen2(ig)=pqsurf(ig,i_n2ice)
-!         DO l=1,nlayer
-!            zt(ig,l)=pt(ig,l)
-!            zq(ig,l,i_n2ice)=pq(ig,l,i_n2ice)
-!            IF( zq(ig,l,i_n2ice).lt.-1.e-6 ) THEN
-!               print*,'Uh-oh, zq = ',zq(ig,l,i_n2ice),'at ig,l=',ig,l
-!               if(l.eq.1)then
-!                  print*,'Perhaps the atmosphere is collapsing on surface...?'
-!               endif
-!            END IF
-!         ENDDO
-!      ENDDO
-
-      ! Calculate the mass of each atmospheric layer (kg.m-2)
-!      do  ilay=1,nlayer
-!         DO ig=1,ngrid
-!            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
-!         end do
-!      end do
+        END IF
+     
+      ENDDO
+    ENDDO
+  endif  ! not fast
+
+!-----------------------------------------------------------------------
+!     Condensation/sublimation on the ground
+!     (calculation of zcondices and pdtsrfc)
+
+!     Adding subtimesteps : in fast version, pressures too low lead to
+!     instabilities.
+  IF (fast) THEN 
+     IF (pplev(1,1).gt.0.3) THEN
+         nsubtimestep= 1  
+     ELSE
+         nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1) 
+     ENDIF
+  ELSE
+     nsubtimestep= 1 ! if more then kp and p00 have to be calculated
+                     ! for nofast mode
+  ENDIF
+  subtimestep=ptimestep/float(nsubtimestep)
+
+  do itsub=1,nsubtimestep
+     ! first loop : getting zplev, ztsurf
+     IF (itsub.eq.1) then
+       DO ig=1,klon
+          zplev(ig)=pplev(ig,1)
+          ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
+          ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep    !!
+          zpicen2(ig)=picen2(ig)
+       ENDDO
+     ELSE
+     ! additional loop : 
+     ! 1)  pressure updated
+     ! 2)  direct redistribution of pressure over the globe
+     ! 3)  modification pressure for unstable cases
+     ! 4)  pressure update to conserve mass
+     ! 5)  temperature updated with radiative tendancies
+       DO ig=1,klon
+          zplevhist(ig)=zplev(ig)
+          zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep   ! 1)
+          !IF (zplevnew(ig).le.0.0001) then
+          !   zplevnew(ig)=0.0001*kp(ig)/p00
+          !ENDIF
+       ENDDO
+       ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
+       globzplevnew=glob_average2d(zplevnew)
+       DO ig=1,klon
+         zplev(ig)=kp(ig)*globzplevnew/p00  ! 2)
+       ENDDO
+       DO ig=1,klon     ! 3) unstable case condensation and sublimation: zplev=zplevhist
+          IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or.  &
+          ((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
+             zplev(ig)=zplevhist(ig)
+          ENDIF
+          zplevhist(ig)=zplev(ig)
+       ENDDO
+       zplev=zplev*globzplevnew/glob_average2d(zplevhist)   ! 4)
+       DO ig=1,klon    ! 5)
+         zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4)
+         ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)*subtimestep+zdtsrf(ig)*subtimestep
+         zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
+       ENDDO
+     ENDIF  ! (itsub=1)
+   
+   DO ig=1,klon
+     ! forecast of frost temperature ztcondsol
+     !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
+     ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
+
+!     Loop over where we have condensation / sublimation
+     IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &    ! ground cond
+                ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &    ! ground sublim
+               (zpicen2(ig) .GT. 0.))) THEN
+        condsub(ig) = .true.    ! condensation or sublimation
+
+!     Condensation or partial sublimation of N2 ice
+        if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
+!             Include a correction to account for the cooling of air near 
+!             the surface before condensing:
+         zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
+         /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep)
+        else                                    ! sublimation
+          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
+          /(lw_n2*subtimestep) 
+        end if
+
+        pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep
+
+!     partial sublimation of N2 ice
+!     If the entire N_2 ice layer sublimes
+!     (including what has just condensed in the atmosphere)
+        IF((zpicen2(ig)/subtimestep).LE. &
+            -zcondices(ig))THEN
+           zcondices(ig) = -zpicen2(ig)/subtimestep 
+           pdtsrfc(ig)=(lw_n2/pcapcal(ig))*       &
+               (zcondices(ig))
+        END IF
+
+!     Changing N2 ice amount and pressure
+
+        pdicen2(ig) = zcondices(ig) 
+        pdpsrf(ig)   = -pdicen2(ig)*g
+    !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
+        IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
+            pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
+            pdicen2(ig)=-pdpsrf(ig)/g
+        ENDIF
+
+     ELSE  ! no condsub
+        pdpsrf(ig)=0.
+        pdicen2(ig)=0.
+        pdtsrfc(ig)=0.
+     ENDIF
+   ENDDO ! ig
+  enddo ! subtimestep
+
+!     Updating pressure, temperature and ice reservoir
+  DO ig=1,klon
+     pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
+     ! Two options here : 1 ok, 2 is wrong
+     pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
+     !pdicen2(ig)=-pdpsrf(ig)/g
+
+     pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep
+
+! security
+     if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
+       write(*,*) 'WARNING in condense_n2:'
+       write(*,*) picen2(ig),pdicen2(ig)*ptimestep
+       pdicen2(ig)= -picen2(ig)/ptimestep
+       pdpsrf(ig)=-pdicen2(ig)*g
+     endif
+
+     if(.not.picen2(ig).ge.0.) THEN
+!           if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
+          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep  
+!              pdicen2(ig)= -picen2(ig)/ptimestep
+!           else
+          picen2(ig)=0.0
+!           endif
+     endif
+  ENDDO
+
+! ***************************************************************
+!  Correction to account for redistribution between sigma or hybrid 
+!  layers when changing surface pressure (and warming/cooling
+!  of the n2 currently changing phase).
+! *************************************************************
+  if (.not.fast) then
+   DO ig=1,klon
+    if (condsub(ig)) then
+
+!  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
+!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+        zmflux(1) = -zcondices(ig)
+        DO l=1,nlayer
+         zmflux(l+1) = zmflux(l) -zcondicea(ig,l)  &
+         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))     
+! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld  
+      if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0.
+        END DO
+
+! Mass of each layer
+! ------------------ 
+       DO l=1,nlayer
+         masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
+       END DO
+
+
+!  Corresponding fluxes for T,U,V,Q
+!  """"""""""""""""""""""""""""""""
+!           averaging operator for TRANSPORT  
+!           """"""""""""""""""""""""""""""""
+
+!           Subtimestep loop to perform the redistribution gently and simultaneously with
+!           the other tendencies  
+!           Estimation of subtimestep (using only  the first layer, the most critical)
+        dtmax=ptimestep
+        if (zmflux(1).gt.1.e-20) then 
+            dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
+        endif
+        nsubtimestep= max(nint(ptimestep/dtmax),nint(2.)) 
+        subtimestep=ptimestep/float(nsubtimestep)
+
+!          New flux for each subtimestep
+       do l=1,nlayer+1
+            w(l)=-zmflux(l)*subtimestep
+       enddo
+!          initializing variables that will vary during subtimestep:
+       do l=1,nlayer
+           ztc(l)  =pt(ig,l)   
+           zu(l)   =pu(ig,l) 
+           zv(l)   =pv(ig,l)  
+           do iq=1,nqmx
+              zq(l,iq) = pq(ig,l,iq) 
+           enddo
+       end do
+
+!          loop over nsubtimestep
+!          """"""""""""""""""""""
+       do itsub=1,nsubtimestep
+!             Progressively adding tendancies from other processes. 
+          do l=1,nlayer
+             ztc(l)  =ztc(l)   +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
+             zu(l)   =zu(l)   +pdu( ig,l) * subtimestep
+             zv(l)   =zv(l)   +pdv( ig,l) * subtimestep
+             do iq=1,nqmx
+                zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
+             enddo
+           end do
+
+!             Change of mass in each layer
+          do l=1,nlayer
+             masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))&
+                                             /(g*pplev(ig,1))
+          end do
+
+           ztm(2:nlayermx+1)=0.
+           zum(2:nlayermx+1)=0.
+           zvm(2:nlayermx+1)=0.
+           zqm1(1:nlayermx+1)=0.
+
+!             Van Leer scheme:
+          call vl1d(nlayer,ztc,2.,masse,w,ztm) 
+          call vl1d(nlayer,zu ,2.,masse,w,zum) 
+          call vl1d(nlayer,zv ,2.,masse,w,zvm) 
+         do iq=1,nqmx
+           do l=1,nlayer
+            zq1(l)=zq(l,iq)
+           enddo
+           zqm1(1)=zqm(1,iq)
+           call vl1d(nlayer,zq1,2.,masse,w,zqm1)
+           do l=2,nlayer
+            zqm(l,iq)=zqm1(l)
+           enddo
+          enddo
+
+!             Correction to prevent negative value for qn2
+          if (igcm_n2.ne.0) then
+            zqm(1,igcm_n2)=1.
+            do l=1,nlayer-1
+              if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then 
+                 zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2),    &
+                 (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) )
+              else 
+                exit
+              endif
+            end do
+           end if
+
+!             Value transfert at the surface interface when condensation sublimation:
+
+          if (zmflux(1).lt.0) then 
+!               Surface condensation
+            zum(1)= zu(1)
+            zvm(1)= zv(1) 
+            ztm(1) = ztc(1)
+          else 
+!               Surface sublimation:
+            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
+            zum(1) = 0  
+            zvm(1) = 0  
+          end if
+          do iq=1,nqmx
+             zqm(1,iq)=0. ! most tracer do not condense !
+          enddo
+!             Special case if the tracer is n2 gas
+          if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.
+
+          !!! Source haze: 0.02 pourcent when n2 sublimes
+          IF (source_haze) THEN
+           IF (pdicen2(ig).lt.0) THEN
+            DO iq=1,nq
+             tname=noms(iq)
+             if (tname(1:4).eq."haze") then 
+                   !zqm(1,iq)=0.02
+                   !zqm(1,iq)=-pdicen2(ig)*0.02
+                   zqm(1,iq)=-pdicen2(ig)*ptimestep*0.02
+                   !zqm(10,iq)=-pdicen2(ig)*ptimestep*100.
+                   !zqm(1,iq)=-pdicen2(ig)*1000000.
+
+             endif
+            ENDDO
+           ENDIF      
+          ENDIF      
+          ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... 
+          zum(nlayer+1)= zu(nlayer)  ! should not be used, but...
+          zvm(nlayer+1)= zv(nlayer)  ! should not be used, but...
+          do iq=1,nqmx
+           zqm(nlayer+1,iq)= zq(nlayer,iq)
+          enddo
+
+!             Tendencies on T, U, V, Q 
+!             """""""""""""""""""""""
+          DO l=1,nlayer
+
+!               Tendencies on T
+            zdtsig(ig,l) = (1/masse(l)) *   &
+           ( zmflux(l)*(ztm(l) - ztc(l))    &
+           - zmflux(l+1)*(ztm(l+1) - ztc(l)) &
+           + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
+
+!               Tendencies on U
+              pduc(ig,l)   = (1/masse(l)) *   &
+           ( zmflux(l)*(zum(l) - zu(l))&
+           - zmflux(l+1)*(zum(l+1) - zu(l)) )
+
+!               Tendencies on V
+              pdvc(ig,l)   = (1/masse(l)) *   &
+           ( zmflux(l)*(zvm(l) - zv(l))   &
+           - zmflux(l+1)*(zvm(l+1) - zv(l)) )
+
+          END DO
+
+!             Tendencies on Q
+          do iq=1,nqmx
+            if (iq.eq.igcm_n2) then
+!                 SPECIAL Case when the tracer IS N2 :
+              DO l=1,nlayer
+              pdqc(ig,l,iq)= (1/masse(l)) *        &
+               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))    &
+               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
+               + zcondicea(ig,l)*(zq(l,iq)-1.) )
+              END DO
+            else
+              DO l=1,nlayer
+                pdqc(ig,l,iq)= (1/masse(l)) *         &
+               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))     &
+               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
+               + zcondicea(ig,l)*zq(l,iq) )
+              END DO
+            end if
+          enddo
+!             Update variables at the end of each subtimestep.
+          do l=1,nlayer
+            ztc(l)  =ztc(l)  + zdtsig(ig,l)  *subtimestep
+            zu(l)   =zu(l)   + pduc(ig,l)  *subtimestep
+            zv(l)   =zv(l)   + pdvc(ig,l)  *subtimestep
+            do iq=1,nqmx
+              zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
+            enddo
+          end do
+       enddo   ! loop on nsubtimestep
+!          Recomputing Total tendencies
+       do l=1,nlayer
+         pdtc(ig,l)   = (ztc(l) - zt(ig,l) )/ptimestep
+         pduc(ig,l)   = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
+         pdvc(ig,l)   = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
+         do iq=1,nqmx
+           pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep
+
+
+!           Correction temporaire
+            if (iq.eq.igcm_n2) then
+             if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
+                     .lt.0.01) then ! if n2 < 1 %  !
+               pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)  
+             end if
+            end if
+
+         enddo
+       end do
+! *******************************TEMPORAIRE ******************
+       if (klon.eq.1) then 
+              write(*,*) 'nsubtimestep=' ,nsubtimestep
+           write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
+              write(*,*) 'masse apres' , masse(1) 
+              write(*,*) 'zmflux*DT, l=1' ,  zmflux(1)*ptimestep
+              write(*,*) 'zmflux*DT, l=2' ,  zmflux(2)*ptimestep
+          write(*,*) 'pq, l=1,2,3' ,  pq(1,1,1), pq(1,2,1),pq(1,3,1)
+          write(*,*) 'zq, l=1,2,3' ,  zq(1,1), zq(2,1),zq(3,1)
+              write(*,*) 'dq*Dt l=1' ,  pdq(1,1,1)*ptimestep
+              write(*,*) 'dqcond*Dt l=1' ,  pdqc(1,1,1)*ptimestep
+       end if
+
+! ***********************************************************
+    end if ! if (condsub)
+   END DO  ! loop on ig 
+  endif ! not fast
+
+! ************ Option Olkin to prevent N2 effect in the south********
+112   continue
+  if (olkin) then
+  DO ig=1,klon
+     if (lati(ig).lt.0) then
+       pdtsrfc(ig) = max(0.,pdtsrfc(ig))
+       pdpsrf(ig) = 0.
+       pdicen2(ig)  = 0.
+       do l=1,nlayer
+         pdtc(ig,l)  =  max(0.,zdtlatent(ig,l))
+         pduc(ig,l)  = 0.
+         pdvc(ig,l)  = 0.
+         do iq=1,nqmx
+           pdqc(ig,l,iq) = 0
+         enddo
+       end do
+     end if
+  END DO
+  end if
+! *******************************************************************
+
+! ***************************************************************
+! Ecriture des diagnostiques
+! ***************************************************************
+
+!     DO l=1,nlayer
+!        DO ig=1,klon
+!         Taux de cond en kg.m-2.pa-1.s-1
+!          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
+!          Taux de cond en kg.m-3.s-1
+!          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
+!        END DO
+!     END DO
+!     call WRITEDIAGFI(ngridmx,'tconda1',              &
+!     'Taux de condensation N2 atmospherique /Pa',    &
+!      'kg.m-2.Pa-1.s-1',3,tconda1)
+!     call WRITEDIAGFI(ngridmx,'tconda2',              &
+!     'Taux de condensation N2 atmospherique /m',     &
+!      'kg.m-3.s-1',3,tconda2)
+
+
+  return
+  end subroutine condens_n2
+
+!-------------------------------------------------------------------------
+
+  real function tcond_n2(p,vmr)
+!   Calculates the condensation temperature for N2 at pressure P and vmr
+  implicit none
+  real, intent(in):: p,vmr
+
+!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
+  IF (p.ge.0.529995) then
+!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
+  ! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e5*0.98)*log(1./(0.125570*1.e5)*p*vmr))
+    tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr))
+  ELSE
+!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
+  ! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e5*1.09)*log(1./(0.508059)*p*vmr))
+    tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr))
+  ENDIF
+  return
+  end function tcond_n2
+
+!-------------------------------------------------------------------------
+
+  real function pcond_n2(t,vmr)
+!   Calculates the condensation pressure for N2 at temperature T and vmr
+  implicit none
+  real, intent(in):: t,vmr
+
+!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
+  IF (t.ge.35.6) then
+!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB 
+  ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t))
+    pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t))
+  ELSE
+!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
+  ! pcond_n2 = 0.508059/vmr*exp((2.5e5*1.09)/296.925*(1./35.6-1./t))
+    pcond_n2 = 0.508059/vmr*exp(917.7401701*(1./35.6-1./t))
+  ENDIF
+  return
+  end function pcond_n2
+
+!-------------------------------------------------------------------------
+
+  real function glob_average2d(var)
+!   Calculates the global average of variable var     
+  use comgeomfi_h
+  implicit none
+
+!      INTEGER klon
+  REAL var(ngridmx)
+  INTEGER ig
+
+  glob_average2d = 0.
+  DO ig=2,ngridmx-1
+       glob_average2d = glob_average2d + var(ig)*area(ig)
+  END DO
+  glob_average2d = glob_average2d + var(1)*area(1)*iim
+  glob_average2d = glob_average2d + var(ngridmx)*area(ngridmx)*iim
+  glob_average2d = glob_average2d/(totarea+(area(1)+area(ngridmx))*(iim-1))
+
+  end function glob_average2d
+
+! *****************************************************************
+
+  subroutine vl1d(nlayer,q,pente_max,masse,w,qm)
+!    
+!     Operateur de moyenne inter-couche pour calcul de transport type
+!     Van-Leer " pseudo amont " dans la verticale
+!    q,w sont des arguments d'entree  pour le s-pg ....
+!    masse : masse de la couche Dp/g
+!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
+!    pente_max = 2 conseillee
+!   --------------------------------------------------------------------
+  IMPLICIT NONE
+
+!   Arguments:
+!   ----------
+  integer nlayer
+  real masse(nlayer),pente_max
+  REAL q(nlayer),qm(nlayer+1)
+  REAL w(nlayer+1)
 !
-
-!-----------------------------------------------------------
-!     START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
-!-----------------------------------------------------------
-
-
-!      Ntime =  20  ! number of sub-timestep 
-!      subptimestep = ptimestep/float(Ntime)           
-
-      ! Add the tendencies from other physical processes at each subtimstep.
-!      DO it=1,Ntime
-!         DO l=1,nlayer
-!            DO ig=1,ngrid
-!               zt(ig,l)   = zt(ig,l)   + pdt(ig,l)   * subptimestep
-!               zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdq(ig,l,i_n2ice) * subptimestep
-!            END DO
-!         END DO
-
-         ! Gravitational sedimentation starts.
-            
-         ! Sedimentation computed from radius computed from q in module radii_mod.
-!	 call n2_reffrad(ngrid,nlayer,nq,zq,reffrad)
-!	 
-!         DO  ilay=1,nlayer
-!            DO ig=1,ngrid
-
-!               reff = reffrad(ig,ilay)
-
-!              call stokes                      &
-!                   (pplev(ig,ilay),pt(ig,ilay), &
-!                    reff,vstokes,rho_n2)
-
-!               !w(ig,ilay,i_n2ice) = 0.0
-!               w(ig,ilay,i_n2ice) = vstokes *  subptimestep * &
-!                   pplev(ig,ilay)/(r*pt(ig,ilay))
-
-!            END DO
-!         END DO
-
-         ! Computing q after sedimentation
-!         call vlz_fi(ngrid,nlayer,zq(1,1,i_n2ice),2.,masse,w(1,1,i_n2ice),wq)
-
-
-         ! Progressively accumulating the flux to the ground.
-         ! Mfallice is the total amount of ice fallen to the ground.
-!         DO ig=1,ngrid
-!            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_n2ice)
-!         END DO
-
-!----------------------------------------------------------
-!       Condensation / sublimation in the atmosphere
-!----------------------------------------------------------
-!     (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation of N2 into tracer i_n2ice)
-
-
-!         DO l=nlayer , 1, -1
-!            DO ig=1,ngrid
-!               pdtc(ig,l)=0.
-
-               ! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
-!               IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_n2ice).gt.1.E-10)) THEN 
-!                  pdtc(ig,l)   = (ztcond(ig,l) - zt(ig,l))/subptimestep
-!                  pdqc(ig,l,i_n2ice) = pdtc(ig,l)*ccond*g
-
-                  ! Case when the ice from above sublimes entirely
-!                  IF ((zq(ig,l,i_n2ice).lt.-pdqc(ig,l,i_n2ice)*subptimestep) &
-!                      .AND. (zq(ig,l,i_n2ice).gt.0)) THEN
-
- !                    pdqc(ig,l,i_n2ice) = -zq(ig,l,i_n2ice)/subptimestep
- !                    pdtc(ig,l)   =-zq(ig,l,i_n2ice)/(ccond*g*subptimestep)
-
-!                  END IF
-
-                  ! Temperature and q after condensation
-!                  zt(ig,l)   = zt(ig,l)   + pdtc(ig,l)   * subptimestep
-!                  zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdqc(ig,l,i_n2ice) * subptimestep
-!               END IF
-
-!            ENDDO
-!         ENDDO
-         
-!      ENDDO! end of subtimestep loop.
-
-      ! Computing global tendencies after the subtimestep.
- !     DO l=1,nlayer
- !        DO ig=1,ngrid
- !           pdtc(ig,l) = &
- !             (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
- !           pdqc(ig,l,i_n2ice) = &
- !             (zq(ig,l,i_n2ice)-(pq(ig,l,i_n2ice)+pdq(ig,l,i_n2ice)*ptimestep))/ptimestep
- !        END DO
- !     END DO
- !     DO ig=1,ngrid
-!         zfallice(ig) = Mfallice(ig)/ptimestep
-!      END DO
-
-
-!-----------------------------------------------------------------------
-!              Condensation/sublimation on the ground
-!-----------------------------------------------------------------------
-
-
-      ! Forecast of ground temperature ztsrf and frost temperature ztcondsol.
-      DO ig=1,ngrid
-         picen2(ig)=pqsurf(ig,i_n2ice)
-         ppn2=gfrac(igas_N2)*pplay(ig,1)
-         call get_tcond_n2(ppn2,ztcondsol(ig))
-         
-         ztsrf(ig) = ptsrf(ig)
-
-         if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
-            print*,'N2 is condensing on the surface in 1D. This atmosphere is doomed.'
-            print*,'T_surf = ',ztsrf,'K'
-            print*,'T_cond = ',ztcondsol,'K'
-            open(116,file='surf_vals.out')
-            write(116,*) 0.0, pplev(1,1), 0.0, 0.0
-            close(116)
-            call abort
+!      Local 
+!   ---------
+!
+  INTEGER l
+!
+  real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax
+  real sigw, Mtot, MQtot
+  integer m 
+
+
+!    On oriente tout dans le sens de la pression 
+!     W > 0 WHEN DOWN !!!!!!!!!!!!!
+
+  do l=2,nlayer
+        dzqw(l)=q(l-1)-q(l)
+        adzqw(l)=abs(dzqw(l))
+  enddo
+
+  do l=2,nlayer-1
+        if(dzqw(l)*dzqw(l+1).gt.0.) then
+            dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
+        else
+            dzq(l)=0.
+        endif
+        dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
+        dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
+  enddo
+
+     dzq(1)=0.
+     dzq(nlayer)=0.
+
+   do l = 1,nlayer-1
+
+!         Regular scheme (transfered mass < layer mass)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
+         sigw=w(l+1)/masse(l+1)
+         qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
+      else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
+         sigw=w(l+1)/masse(l)
+         qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))
+
+!         Extended scheme (transfered mass > layer mass)
+!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      else if(w(l+1).gt.0.) then
+         m=l+1
+         Mtot = masse(m)
+         MQtot = masse(m)*q(m)
+         !if (m.lt.nlayer) then ! because some compilers will have problems
+         !                      ! evaluating masse(nlayer+1)
+         do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1))))
+            m=m+1
+            Mtot = Mtot + masse(m)
+            MQtot = MQtot + masse(m)*q(m)
+         !   if (m.eq.nlayer) exit
+         end do
+         !endif
+         if (m.lt.nlayer) then
+            sigw=(w(l+1)-Mtot)/masse(m+1)
+            qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* &
+           (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
+         else
+            w(l+1) = Mtot
+            qm(l+1) = Mqtot / Mtot
+            write(*,*) 'top layer is disapearing !'
+            stop
+         end if
+      else      ! if(w(l+1).lt.0) 
+         m = l-1 
+         Mtot = masse(m+1)
+         MQtot = masse(m+1)*q(m+1)
+         if (m.gt.0) then ! because some compilers will have problems
+                          ! evaluating masse(0)
+          do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
+            m=m-1
+            Mtot = Mtot + masse(m+1)
+            MQtot = MQtot + masse(m+1)*q(m+1)
+            if (m.eq.0) exit
+          end do
          endif
-
-         ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
-
-      ENDDO
-     
-      DO ig=1,ngrid
-      
-         IF(ig.GT.ngrid/2+1) THEN
-            icap=2
-         ELSE
-            icap=1
-         ENDIF
-
-         ! Loop over where we have condensation / sublimation
-         IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &        ! ground condensation
-                    (zfallice(ig).NE.0.) .OR.              &        ! falling snow
-                    ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &        ! ground sublimation
-                    ((picen2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
-
-
-            ! Condensation or partial sublimation of N2 ice
-            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
-                /(latcond*ptimestep)         
-            pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-
-            ! If the entire CO_2 ice layer sublimes
-            ! (including what has just condensed in the atmosphere)
-            IF((picen2(ig)/ptimestep+zfallice(ig)).LE. &
-                -zcondices(ig))THEN
-               zcondices(ig) = -picen2(ig)/ptimestep - zfallice(ig)
-               pdtsrfc(ig)=(latcond/pcapcal(ig))*       &
-                   (zcondices(ig))
-            END IF
-
-            ! Changing N2 ice amount and pressure
-            picen2(ig)  =  picen2(ig)   + pdqsurfc(ig)*ptimestep
-            pdqsurfc(ig) =  zcondices(ig) + zfallice(ig)
-            pdpsrfc(ig)  = -pdqsurfc(ig)*g
-
-            IF(ABS(pdpsrfc(ig)*ptimestep).GT.pplev(ig,1)) THEN
-               PRINT*,'STOP in condens in condense_n2'
-               PRINT*,'condensing more than total mass'
-               PRINT*,'Grid point ',ig
-               PRINT*,'Ps = ',pplev(ig,1)
-               PRINT*,'d Ps = ',pdpsrfc(ig)
-               STOP
-            ENDIF
-         END IF
-         
-      ENDDO ! end of ngrid loop.
-
-
-!---------------------------------------------------------------------------------------------
-!     Surface albedo and emissivity of the ground below the snow (emisref)
-!---------------------------------------------------------------------------------------------
-
-
-      DO ig=1,ngrid
-      
-!         IF(latitude(ig).LT.0.) THEN
-!            icap=2 ! Southern Hemisphere
-!         ELSE
-!            icap=1 ! Nortnern hemisphere
-!         ENDIF
-
-         if(.not.picen2(ig).ge.0.) THEN
-            if(picen2(ig).le.-1.e-8) print*,   &
-              'WARNING : in condense_n2cloud: picen2(',ig,')=', picen2(ig)
-            picen2(ig)=0.
-         endif
-!         if (picen2(ig) .gt. 1.) then  ! N2 Albedo condition changed to ~1 mm coverage. Change by MT2015.
-!	    DO nw=1,L_NSPECTV
-!               albedo(ig,nw) = albedo_n2_ice_SPECTV(nw)
-!	    ENDDO
-!            emisref(ig)   = emisice(icap)
-!         else
-!            DO nw=1,L_NSPECTV
-!               albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.
-!	    ENDDO
-!            emisref(ig)   = emissiv
-!            pemisurf(ig)  = emissiv
-!         end if
-         
-      END DO
-
-      return
-      
-      end subroutine condense_n2
-
-
-
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-
-
-
-      subroutine get_tcond_n2(p,tcond)  ! Calculates the condensation temperature for N2
-     
-
-      implicit none
-
-      real p, tcond
-
-      IF (p.ge.0.529995) then
-     ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
-         tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr))
-      ELSE
-     ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
-         tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr))
-      ENDIF
-      return
-
-      end subroutine get_tcond_n2
-
-
-
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-!-------------------------------------------------------------------------
-
-
-
-      subroutine get_tnuc_n2(p,tnuc)
-      ! Calculates the nucleation temperature for N2, based on a simple super saturation criterion. JL 2011.
-
-      use callkeys_mod, only: n2supsat
-
-      implicit none
-
-      real p, peff, tnuc
-      real, parameter :: ptriple=518000.0
-
-      peff=p/n2supsat
-
-      !! AF24: this code below is for CO2, needs to be udpated for N2
-
-      ! if(peff.lt.ptriple) then
-      !    tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
-      ! else
-      !    tnuc = 684.2-92.3*log(peff)+4.32*log(peff)**2 
-      !    ! liquid-vapour transition (based on CRC handbook 2003 data)
-      ! endif
-      
-      return
-
-      end subroutine get_tnuc_n2
+         if (m.gt.0) then
+            sigw=(w(l+1)+Mtot)/masse(m)
+            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* & 
+           (q(m)-0.5*(1.+sigw)*dzq(m))  )
+         else
+            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
+         end if   
+      end if
+   enddo
+
+   return
+   end  subroutine vl1d
+
Index: trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90	(revision 3193)
@@ -21,4 +21,7 @@
       character(LEN=18),parameter :: aerdir="aerosol_properties" 
 
+      ! Data haze properties
+      character(len=300),save :: hazeprop 
+      
       end module datafile_mod
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90	(revision 3193)
@@ -13,5 +13,5 @@
   use radcommon_h, only: ini_radcommon_h
   use radii_mod, only: radfixed, Nmix_n2
-  use datafile_mod, only: datadir
+  use datafile_mod, only: datadir, hazeprop
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
   use comgeomfi_h, only: totarea, totarea_planet
@@ -136,4 +136,10 @@
      if (is_master) write(*,*) trim(rname)//": datadir = ",trim(datadir)
 
+     if (is_master) write(*,*) "Haze optical properties datafile"
+     hazeprop="../../../datadir"  ! default path
+     call getin_p("hazeprop",hazeprop) 
+     if (is_master) write(*,*) trim(rname)//" hazeprop = ",trim(hazeprop)
+
+
      if (is_master) write(*,*) trim(rname)//&
        ": Run with or without tracer transport ?"
@@ -788,25 +794,4 @@
      if (is_master) write(*,*)trim(rname)//": generic_condensation = ",generic_condensation
      
-     if (is_master) write(*,*)trim(rname)//": Generic Condensation of tracers ?"
-     generic_rain=.false. !default value 
-     call getin_p("generic_rain",generic_rain)
-     if (is_master) write(*,*)trim(rname)//": generic_rain = ",generic_rain
-
-     if (is_master) write(*,*)trim(rname)//": Compute water cycle ?"
-     water=.false. ! default value
-     call getin_p("water",water)
-     if (is_master) write(*,*)trim(rname)//": water = ",water
-         
-! Test of incompatibility:
-! if water is true, there should be at least a tracer
-     if (water.and.(.not.tracer)) then
-        call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
-     endif
-
-     if (is_master) write(*,*)trim(rname)//": Include water condensation ?"
-     watercond=.false. ! default value
-     call getin_p("watercond",watercond)
-     if (is_master) write(*,*)trim(rname)//": watercond = ",watercond
-
 ! Test of incompatibility:
      if (is_master) write(*,*)trim(rname)//": Spectral Dependant albedo ?"
Index: trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/initracer.F90	(revision 3193)
@@ -520,4 +520,38 @@
       rho_n2=1000.   ! N2 ice density (kg.m-3)
 
+      lw_co=274000.
+      lw_ch4=586700.
+      lw_n2=2.5e5
+      
+      if (haze) then
+        ! the sedimentation radius remains radius(igcm_haze)
+        if (fractal) then
+           nmono=nb_monomer
+        else
+           nmono=1
+        endif  
+
+        ia=0
+        if (aerohaze) then
+           ia=ia+1
+           iaero_haze=ia
+           write(*,*) '--- number of haze aerosol = ', iaero_haze
+
+           block=0  ! Only one type of haze is active : the first one set in traceur.def
+           do iq=1,nqmx
+             tracername=noms(iq)
+             write(*,*) "--> tracername ",iq,'/',nqmx,' = ',tracername
+             if (tracername(1:4).eq."haze".and.block.eq.0) then
+               i_haze=iq
+               block=1
+               write(*,*) "i_haze=",i_haze
+               write(*,*) "Careful: if you set many haze traceurs in 
+    & traceur.def,only ",tracername," will be radiatively active 
+    & (first one in traceur.def)"
+             endif
+           enddo
+        endif
+     endif
+
 !     Initialization for water vapor !AF24: removed
 
Index: trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90
===================================================================
--- trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90	(revision 3187)
+++ trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90	(revision 3193)
@@ -10,5 +10,5 @@
       use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
       use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
-      use datafile_mod, only: datadir, aerdir
+      use datafile_mod, only: datadir, aerdir, hazeprop
 
       ! outputs
@@ -136,139 +136,30 @@
         print*, 'naerkind= 0'
       endif
+      
+      
       do iaer=1,naerkind
-       if (iaer.eq.iaero_n2) then
-        forgetit=.true.
-          if (.not.noaero) then
-              print*, 'naerkind= n2', iaer
-          end if
-!     visible
-        if(forgetit)then
-           file_id(iaer,1) = 'optprop_n2_vis_ff.dat' ! Francois' values
-        else
-           file_id(iaer,1) = 'optprop_n2ice_vis_n50.dat'
-        endif
-        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
-
-!     infrared
-        if(forgetit)then
-           file_id(iaer,2) = 'optprop_n2_ir_ff.dat' ! Francois' values
-        else
-           file_id(iaer,2) = 'optprop_n2ice_ir_n50.dat'
-        endif
-        lamrefir(iaer)=12.1E-6   ! Dummy in generic phys. (JVO 20)
-       endif ! N2 aerosols  
-! NOTE: these lamref's are currently for dust, not N2 ice.
-! JB thinks this shouldn't matter too much, but it needs to be 
-! fixed / decided for the final version.
-
-       if (iaer.eq.iaero_dust) then
-        print*, 'naerkind= dust', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
-         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
-         lamrefvis(iaer)=0.67e-6
-!     infrared
-         file_id(iaer,2) = 'optprop_dustir_n50.dat'
-         lamrefir(iaer)=9.3E-6     ! Dummy in generic phys. (JVO 20)
-       endif 
-
-       if (iaer.eq.iaero_h2so4) then
-         print*, 'naerkind= h2so4', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
-         lamrefvis(iaer)=1.5E-6   ! no idea, must find
-!     infrared
-         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
-         lamrefir(iaer)=9.3E-6 ! ! Dummy in generic phys. (JVO 20)
-! added by LK
-       endif
-
-       if (iaer.eq.iaero_back2lay) then
-         print*, 'naerkind= back2lay', iaer
+         !     naerkind=1: haze
+         if (iaer.eq.1) then
+
+         ! Only one table of optical properties : 
+         write(*,*)'Suaer haze optical properties, using: ', &
+                           TRIM(hazeprop)
+
+         ! visible
+         file_id(iaer,1)=TRIM(hazeprop)
+         ! infrared
+         file_id(iaer,2)=file_id(iaer,1)
          
-!     visible
-         file_id(iaer,1) = TRIM(optprop_back2lay_vis)
-         lamrefvis(iaer)=0.8E-6  ! This is the important one.
-!     infrared
-         file_id(iaer,2) = TRIM(optprop_back2lay_ir)
-         lamrefir(iaer)=6.E-6    ! This is dummy.
-! added by SG
-       endif
-      
-       if (iaer.eq.iaero_nh3) then
-         print*, 'naerkind= nh3', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_nh3ice_vis.dat'
-         lamrefvis(iaer)=0.756E-6  ! 
-!     infrared
-         file_id(iaer,2) = 'optprop_nh3ice_ir.dat'
-         lamrefir(iaer)=6.E-6  ! dummy 
-! added by SG
-       endif
-
-       do ia=1,nlayaero
-         if (iaer.eq.iaero_nlay(ia)) then
-           print*, 'naerkind= nlay', iaer
-           
-!       visible
-           file_id(iaer,1) = TRIM(optprop_aeronlay_vis(ia))
-           lamrefvis(iaer)=aeronlay_lamref(ia)
-!       infrared
-           file_id(iaer,2) = TRIM(optprop_aeronlay_ir(ia))
-           lamrefir(iaer)=6.E-6 ! Dummy value
+         lamrefvis(iaer)=0.607E-6   ! reference wavelength for opacity vis pivot wavelength Cheng et al 2017
+         lamrefir(iaer)=2.E-6   !  reference wavelength for opacity IR (in the LEISA range)
+
          endif
-       enddo
-! added by JVO
-      
-       if (iaer.eq.iaero_aurora) then
-         print*, 'naerkind= aurora', iaer
-
-!     visible
-         file_id(iaer,1) = 'optprop_aurora_vis.dat'
-         lamrefvis(iaer)=0.3E-6  ! 
-!     infrared
-         file_id(iaer,2) = 'optprop_aurora_ir.dat'
-         lamrefir(iaer)=6.E-6  ! dummy 
-! added by SG
-       endif
-
-       
-! the following was added by LT
-       do ia=1,aerogeneric ! Read Radiative Generic Condensable Species data
-         if (iaer .eq. iaero_generic(ia)) then 
-            if (index(noms(i_rgcs_ice(ia)),'Fe') .ne. 0) then 
-               print*,"Reading Fe file"
-               file_id(iaer,1)='Fe.ocst.txt'
-               file_id(iaer,2)='Fe.ocst.txt'
-               lamrefvis(iaer) = 1.0E-6 !random pick
-               lamrefir(iaer) = 1.0E-6 !dummy but random pick
-            else if (index(noms(i_rgcs_ice(ia)),'Mn') .ne. 0) then 
-               print*,"Reading MnS file"
-               file_id(iaer,1)='MnS.ocst.txt'
-               file_id(iaer,2)='MnS.ocst.txt'
-               lamrefvis(iaer) = 1.0E-6 !random pick
-               lamrefir(iaer) = 1.0E-6 !dummy but random pick   
-            else if (index(noms(i_rgcs_ice(ia)),'Mg') .ne. 0) then  
-               print*,"Reading Mg2SiO4 file" 
-               file_id(iaer,1)='Mg2SiO4.ocst.txt'
-               file_id(iaer,2)='Mg2SiO4.ocst.txt'
-               lamrefvis(iaer) = 1.0E-6 !random pick
-               lamrefir(iaer) = 1.0E-6 !dummy but random pick  
-            else if (index(noms(i_rgcs_ice(ia)),'Cr') .ne. 0) then
-               print*,"Reading Cr file"
-               file_id(iaer,1)='Cr.ocst.txt'
-               file_id(iaer,2)='Cr.ocst.txt'
-               lamrefvis(iaer) = 1.0E-6 !random pick
-               lamrefir(iaer) = 1.0E-6 !dummy but random pick
-            else
-! If you want to add another specie, copy,paste & adapt the elseif block up here to your new specie (LT 2022)
-               call abort_physic("suaer_corrk", "Unknown specie in radiative generic condensable species",1)
-            endif
-         endif
-       enddo ! ia=1,aerogeneric
-      enddo ! of do iaer=1,naerkind
+      enddo   
+
+      IF (naerkind .gt. 1) THEN
+         print*,'naerkind = ',naerkind
+         print*,'but we only have data for 1 type, exiting.'
+         call abort
+      ENDIF
       
 !------------------------------------------------------------------
