Index: /trunk/LMDZ.TITAN/libf/muphytitan/fsystem.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/muphytitan/fsystem.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/muphytitan/fsystem.F90	(revision 3318)
@@ -428,4 +428,5 @@
     CHARACTER(len=:), ALLOCATABLE :: opath
       !! A Fortran allocated string with the parent directory path or an empty string if method fails
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     TYPE(C_PTR) :: zpath
     IF (LEN_TRIM(path) == 0) THEN
@@ -433,5 +434,6 @@
       RETURN
     ENDIF
-    zpath = dirname_c(cstring(ADJUSTL(path)))
+    cpath = cstring(ADJUSTL(path))
+    zpath = dirname_c(cpath)
     IF (.NOT.C_ASSOCIATED(zpath)) THEN
       opath = ""
@@ -448,4 +450,5 @@
     CHARACTER(len=:), ALLOCATABLE :: opath
       !! The basename of the path or an empty string if method fails
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     TYPE(C_PTR) :: zpath
     IF (LEN_TRIM(path) == 0) THEN
@@ -453,5 +456,6 @@
       RETURN
     ENDIF
-    zpath = basename_c(cstring(ADJUSTL(path)))
+    cpath = cstring(ADJUSTL(path))
+    zpath = basename_c(cpath)
     IF (.NOT.C_ASSOCIATED(zpath)) THEN
       opath = ""
@@ -472,6 +476,8 @@
     CHARACTER(len=:), ALLOCATABLE :: opath
       !! The absolute of the path or an empty string if method fails
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     TYPE(C_PTR) :: zpath
-    zpath = realpath_c(cstring(ADJUSTL(path)))
+    cpath = cstring(ADJUSTL(path))
+    zpath = realpath_c(cpath)
     IF (.NOT.C_ASSOCIATED(zpath)) THEN
       opath = ""
@@ -490,6 +496,9 @@
                                     reldir  !! A directory path from which output should be relative to
     CHARACTER(len=:), ALLOCATABLE :: res    !! An allocated string with the resulting path
+    CHARACTER(len=:), ALLOCATABLE :: cpath1,cpath2
     TYPE(C_PTR) :: zpath
-    zpath = relpath_c(cstring(ADJUSTL(path)),cstring(ADJUSTL(reldir)))
+    cpath1 = cstring(ADJUSTL(path))
+    cpath2 = cstring(ADJUSTL(reldir))
+    zpath = relpath_c(cpath1,cpath2)
     IF (.NOT.C_ASSOCIATED(zpath)) THEN
       res = TRIM(ADJUSTL(path))
@@ -520,8 +529,12 @@
     CHARACTER(len=*), INTENT(in)  :: output !! Output file path destination.
     LOGICAL :: ret                          !! True on success, false otherwise.
+    CHARACTER(len=:), ALLOCATABLE :: cpath1,cpath2
+
     IF (LEN_TRIM(input) == 0 .OR. LEN_TRIM(output) == 0 .OR. input == output) THEN
       ret = .false.
     ELSE
-      ret = INT(copy_c(cstring(ADJUSTL(output)),cstring(ADJUSTL(input)))) == 0
+      cpath1 = cstring(ADJUSTL(output))
+      cpath2 = cstring(ADJUSTL(input))
+      ret = INT(copy_c(cpath1,cpath2)) == 0
     ENDIF
     RETURN
@@ -532,8 +545,10 @@
     CHARACTER(len=*), INTENT(in)  :: path !! A string with the (valid) file path to delete
     LOGICAL :: ret                        !! True on success, false otherwise.
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
     ELSE
-      ret = INT(remove_c(cstring(ADJUSTL(path)))) == 0
+      cpath = cstring(ADJUSTL(path))
+      ret = INT(remove_c(cpath)) == 0
     ENDIF
     RETURN
@@ -545,8 +560,11 @@
                                     new    !! A string with the new name of the path
     LOGICAL :: ret                         !! True on success, false otherwise.
+    CHARACTER(len=:), ALLOCATABLE :: cpath1,cpath2
     IF (LEN_TRIM(old) == 0.OR.LEN_TRIM(new) == 0) THEN
       ret = .false.
     ELSE
-      ret = INT(rename_c(cstring(ADJUSTL(old)),cstring(ADJUSTL(new)))) == 0
+      cpath1 = cstring(ADJUSTL(old))
+      cpath2 = cstring(ADJUSTL(new))
+      ret = INT(rename_c(cpath1,cpath2)) == 0
     ENDIF
     RETURN
@@ -559,9 +577,11 @@
     LOGICAL  :: ret                      !! True on success, false otherwise.
     INTEGER(kind=C_INT) :: zmode
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
     ELSE
       zmode = INT(oct_2_dec(mode),kind=C_INT)
-      ret = INT(chmod_c(cstring(ADJUSTL(path)), zmode)) == 0
+      cpath = cstring(ADJUSTL(path))
+      ret = INT(chmod_c(cpath, zmode)) == 0
     ENDIF
     RETURN
@@ -572,8 +592,10 @@
     CHARACTER(len=*), INTENT(in) :: path !! Path of the new working directory
     LOGICAL :: ret                       !! True on success, false otherwise.
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
     ELSE
-      ret = INT(chdir_c(cstring(ADJUSTL(path)))) == 0
+      cpath = cstring(ADJUSTL(path))
+      ret = INT(chdir_c(cpath)) == 0
     ENDIF
     RETURN
@@ -595,4 +617,6 @@
     INTEGER :: zmode
     LOGICAL :: zperm
+    CHARACTER(len=:), ALLOCATABLE :: cpath
+
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
@@ -605,9 +629,10 @@
         zmode = oct_2_dec(mode)
       ENDIF
+      cpath = cstring(ADJUSTL(path))
       zperm = .false. ; IF (PRESENT(permissive)) zperm = permissive
       IF (zperm) THEN
-        ret = INT(mkdirp_c(cstring(ADJUSTL(path)),INT(zmode,kind=C_INT))) == 0
+        ret = INT(mkdirp_c(cpath,INT(zmode,kind=C_INT))) == 0
       ELSE
-        ret = INT(mkdir_c(cstring(ADJUSTL(path)),INT(zmode,kind=C_INT))) == 0
+        ret = INT(mkdir_c(cpath,INT(zmode,kind=C_INT))) == 0
       ENDIF
     ENDIF
@@ -627,12 +652,14 @@
       !! True on success, false otherwise.
     LOGICAL :: zforce
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
     ELSE
       zforce = .false. ; IF (PRESENT(forced)) zforce = forced
+      cpath = cstring(ADJUSTL(path))
       IF (.NOT.zforce) THEN
-        ret = INT(rmdir_c(cstring(ADJUSTL(path)))) == 0
+        ret = INT(rmdir_c(cpath)) == 0
       ELSE
-        ret = INT(rmdirf_c(cstring(ADJUSTL(path)))) == 0
+        ret = INT(rmdirf_c(cpath)) == 0
       ENDIF
     ENDIF
@@ -668,4 +695,5 @@
     INTEGER(kind=c_long)          :: f
     CHARACTER(len=20,kind=C_CHAR) :: ta,tm,tc
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.; RETURN
@@ -677,5 +705,6 @@
       ! set default values
       pe=-1 ; ty=-1 ; ud=-1 ; gd=-1 ; fs=-1 ; at="" ; mt="" ; ct=""
-      ret = INT(fstat_c(cstring(ADJUSTL(path)),p,l,t,u,g,f,ta,tm,tc)) == 0
+      cpath = cstring(ADJUSTL(path))
+      ret = INT(fstat_c(cpath,p,l,t,u,g,f,ta,tm,tc)) == 0
       IF (ret) THEN
         pe=INT(p) ; ln=INT(l) ; ty=INT(t) ; ud=INT(u) ; gd=INT(g)
@@ -752,4 +781,5 @@
     LOGICAL :: ret                              !! True on success, false otherwise.
     INTEGER(kind=C_INT) :: zp
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     IF (LEN_TRIM(path) == 0) THEN
       ret = .false.
@@ -757,5 +787,6 @@
       zp = 0 ; IF (PRESENT(permission)) zp = INT(permission,kind=C_INT)
       ! Defaults are set in the C function.
-      ret = INT(access_c(cstring(ADJUSTL(path)),zp)) == 0
+      cpath = cstring(ADJUSTL(path))
+      ret = INT(access_c(cpath,zp)) == 0
     ENDIF
     RETURN
@@ -822,4 +853,5 @@
     INTEGER                       :: zmd,zt,zp
     CHARACTER(len=:), ALLOCATABLE :: b,e
+    CHARACTER(len=:), ALLOCATABLE :: cpath
     ret = .false.
     ! Checking for existence
@@ -856,5 +888,7 @@
     ENDIF
     zp = 0 ; IF(PRESENT(permissive)) THEN ; IF(permissive) zp=1 ; ENDIF
-    ret = INT(create_c(cstring(ADJUSTL(path)),INT(zmd,kind=C_INT),INT(zt,kind=C_INT),INT(zp,kind=C_INT))) == 0
+
+    cpath = cstring(ADJUSTL(path))
+    ret = INT(create_c(cpath,INT(zmd,kind=C_INT),INT(zt,kind=C_INT),INT(zp,kind=C_INT))) == 0
     RETURN
   END FUNCTION fs_create
Index: /trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90
===================================================================
--- /trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90	(revision 3318)
@@ -112,12 +112,12 @@
       mm_ccn_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
 
-      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ccn
+      ! Computes flux [kg.m-2.s-1] and precipitation [kg.m-2.iphysiq] of ccn
       mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
-      mm_ccn_prec = SUM(zdm3n*mm_dzlev)
-
-      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ices
+      mm_ccn_prec = SUM(zdm3n*mm_dzlev*mm_rhoaer)
+
+      ! Computes flux [kg.m-2.s-1] and precipitation [kg.m-2.iphysiq] of ices
       DO i = 1, mm_nesp
         mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,(3._mm_wp*mm_m3ice(:,i))/(4._mm_wp*mm_pi))
-        mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev)
+        mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev*mm_xESPS(i)%rho)
       ENDDO
 
@@ -258,4 +258,5 @@
     ! Saturation ratio
     Xsat = zvapX / qsat
+
     
     ! Gets nucleation rate (ccn radius is the monomer !)
@@ -740,6 +741,6 @@
     Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_g(t))
     
-    ! Computes settling velocity (correction factor : x2.0)
-    w = Us * Fc * 2._mm_wp
+    ! Computes settling velocity (correction factor : x3.0)
+    w = Us * Fc * 3._mm_wp
   END FUNCTION wsettle
 
Index: /trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90
===================================================================
--- /trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90	(revision 3318)
@@ -227,7 +227,7 @@
   REAL(kind=mm_wp), PARAMETER :: mm_rgas = mm_kboltz * mm_navo
   !> Desorption energy (\(J\)) (nucleation).
-  REAL(kind=mm_wp), PARAMETER :: mm_fdes = 0.288e-19_mm_wp
+  REAL(kind=mm_wp), PARAMETER :: mm_fdes = 1.519e-20_mm_wp
   !> Surface diffusion energy (\(J\)) (nucleation).
-  REAL(kind=mm_wp), PARAMETER :: mm_fdif = 0.288e-20_mm_wp
+  REAL(kind=mm_wp), PARAMETER :: mm_fdif = 1.519e-21_mm_wp
   !> Jump frequency (\(s^{-1}\)) (nucleation).
   REAL(kind=mm_wp), PARAMETER :: mm_nus = 1.e+13_mm_wp
@@ -429,5 +429,5 @@
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_drho
 
-  !> Aerosols precipitations (m).
+  !> Aerosols precipitations (kg.m-2.s-1).
   !!
   !! Aerosols precipitations take into account both spherical and fractal modes.
@@ -435,5 +435,5 @@
   REAL(kind=mm_wp), SAVE :: mm_aer_prec = 0._mm_wp
 
-  !> CCN precipitations (m).
+  !> CCN precipitations (kg.m-2.s-1).
   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
   REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
@@ -505,5 +505,5 @@
   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux
 
-  !> Ice components precipitations (m).
+  !> Ice components precipitations (kg.m-2.s-1).
   !!
   !! It is a vector of [[mm_globals(module):mm_nesp(variable)]] values which share the same indexing
@@ -1431,6 +1431,6 @@
     ! Initialization :
     Ntot = m0ccn
-    Vtot = pifac*m3ccn + SUM(m3ice)
-    Wtot = pifac*m3ccn*mm_rhoaer + SUM(m3ice*mm_xESPS(:)%rho)
+    Vtot = pifac*m3ccn + pifac*SUM(m3ice)
+    Wtot = pifac*m3ccn*mm_rhoaer + pifac*SUM(m3ice*mm_xESPS(:)%rho)
 
     IF (Ntot <= mm_m0n_min .OR. Vtot <= mm_m3cld_min) THEN
Index: /trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90
===================================================================
--- /trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90	(revision 3318)
@@ -117,5 +117,5 @@
 
       ! Computes precipitations
-      mm_aer_prec = SUM(zdm3as*mm_dzlev) + SUM(zdm3af*mm_dzlev)
+      mm_aer_prec = SUM(zdm3as*mm_dzlev*mm_rhoaer) + SUM(zdm3af*mm_dzlev*mm_rhoaer)
 
       ! Updates tendencies
Index: /trunk/LMDZ.TITAN/libf/muphytitan/mm_microphysic.f90
===================================================================
--- /trunk/LMDZ.TITAN/libf/muphytitan/mm_microphysic.f90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/muphytitan/mm_microphysic.f90	(revision 3318)
@@ -163,5 +163,5 @@
   END FUNCTION muphys_nocld
 
-  SUBROUTINE mm_diagnostics(aer_prec,aer_s_w,aer_f_w,aer_s_flux,aer_f_flux,ccn_prec,ccn_w,ccn_flux,ice_prec,ice_fluxes,gazs_sat)
+  SUBROUTINE mm_diagnostics(dt,aer_prec,aer_s_w,aer_f_w,aer_s_flux,aer_f_flux,ccn_prec,ccn_w,ccn_flux,ice_prec,ice_fluxes,gazs_sat)
     !! Get various diagnostic fields of the microphysics.
     !!
@@ -185,6 +185,7 @@
     !! __ccnprec__, __iceprec__, __icefluxes__ and __gazsat__ are always set to 0 if clouds
     !! microphysics is disabled (see [[mm_globals(module):mm_w_clouds(variable)]] documentation).
-    REAL(kind=mm_wp), INTENT(out), OPTIONAL                 :: aer_prec   !! Aerosols precipitations (both modes) (m).
-    REAL(kind=mm_wp), INTENT(out), OPTIONAL                 :: ccn_prec   !! CCN precipitations (m).
+    REAL(kind=8), INTENT(IN)                                :: dt         !! Physics timestep (s).
+    REAL(kind=mm_wp), INTENT(out), OPTIONAL                 :: aer_prec   !! Aerosols precipitations (both modes) (kg.m-2.s-1).
+    REAL(kind=mm_wp), INTENT(out), OPTIONAL                 :: ccn_prec   !! CCN precipitations (kg.m-2.s-1).
     REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:)   :: aer_s_w    !! Spherical aerosol settling velocity (\(m.s^{-1}\)).
     REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:)   :: aer_f_w    !! Fractal aerosol settling velocity (\(m.s^{-1}\)).
@@ -195,7 +196,7 @@
     REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:,:) :: ice_fluxes !! Ice sedimentation fluxes (\(kg.m^{-2}.s^{-1}\)).
     REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:,:) :: gazs_sat   !! Condensible gaz saturation ratios (--).
-    REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:)   :: ice_prec   !! Ice precipitations (m).
-
-    IF (PRESENT(aer_prec))   aer_prec   = ABS(mm_aer_prec)
+    REAL(kind=mm_wp), INTENT(out), OPTIONAL, DIMENSION(:)   :: ice_prec   !! Ice precipitations (kg.m-2.s-1).
+
+    IF (PRESENT(aer_prec))   aer_prec   = ABS(mm_aer_prec) / dt
     IF (PRESENT(aer_s_w))    aer_s_w    = -mm_m3as_vsed(mm_nla:1:-1)
     IF (PRESENT(aer_f_w))    aer_f_w    = -mm_m3af_vsed(mm_nla:1:-1)
@@ -204,6 +205,6 @@
 
     IF (mm_w_clouds) THEN
-      IF (PRESENT(ccn_prec))   ccn_prec   = ABS(mm_ccn_prec)
-      IF (PRESENT(ice_prec))   ice_prec   = ABS(mm_ice_prec)
+      IF (PRESENT(ccn_prec))   ccn_prec   = ABS(mm_ccn_prec) / dt
+      IF (PRESENT(ice_prec))   ice_prec   = ABS(mm_ice_prec) / dt
       IF (PRESENT(ccn_w))      ccn_w      = mm_ccn_vsed(mm_nla:1:-1)
       IF (PRESENT(ccn_flux))   ccn_flux   = mm_ccn_flux(mm_nla:1:-1)
Index: /trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90	(revision 3318)
@@ -82,18 +82,22 @@
       else if(trim(cnames(ic)).eq."C2H4") then
          ! not far from Fray and Schmidt (2009)
-         where (temp(:,:).lt.89.0)
-            ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
-                  * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                      &
-                  / press(:,:) * 1.01325e0 / 760.0
-         elsewhere (temp(:,:).lt.104.0)
-            ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
-                  / press(:,:) * 1013.25e0 / 760.0
-         elsewhere (temp(:,:).lt.120.0)
-            ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
-                  * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
-         elsewhere (temp(:,:).lt.155.0)
-            ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
-                  / press(:,:) * 1013.25e0 / 760.0
-         endwhere
+         !where (temp(:,:).lt.89.0)
+         !   ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
+         !         * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                      &
+         !         / press(:,:) * 1.01325e0 / 760.0
+         !elsewhere (temp(:,:).lt.104.0)
+         !   ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
+         !         / press(:,:) * 1013.25e0 / 760.0
+         !elsewhere (temp(:,:).lt.120.0)
+         !   ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
+         !         * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
+         !elsewhere (temp(:,:).lt.155.0)
+         !   ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
+         !         / press(:,:) * 1013.25e0 / 760.0
+         !endwhere
+
+         ! Fray and Schmidt (2009)
+         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.540e1 - 2.206e3/temp(:,:)                &
+                         - 1.216e4/temp(:,:)**2 + 2.843e5/temp(:,:)**3 - 2.203e6/temp(:,:)**4)
       
       ! C2H6 :
Index: /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 3318)
@@ -62,4 +62,5 @@
 
   USE, INTRINSIC :: iso_c_binding
+  use tracer_h
   USE comchem_h
   USE dimphy
@@ -334,4 +335,7 @@
        ENDDO
      endif
+
+     ! [BdBdT : On force la chimie... kedd nul dans le PCM]
+     kedd(:klev) = 1.e3 ! Default value =/= zero
  
      firstcall = .FALSE.
@@ -521,5 +525,5 @@
         DO ic=1,nkim
            DO l=1,klev
-              dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
+            dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
            ENDDO
         ENDDO
Index: /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 3318)
@@ -6,5 +6,5 @@
           fluxabs_sw,fluxtop_dn,                               &
           OLR_nu,OSR_nu,                                       &
-          int_dtaui,int_dtauv,popthi,popthv,poptci,poptcv,     &
+          int_dtaui,int_dtauv,popthi,popthv,poptti,popttv,     &
           lastcall)
 
@@ -88,8 +88,8 @@
       REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags ().
       ! Diagnostics
-      REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g).
-      REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g).
-      REAL,INTENT(OUT) :: poptci(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
-      REAL,INTENT(OUT) :: poptcv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
+      REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,8)  ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
+      REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,8)  ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
+      REAL,INTENT(OUT) :: poptti(ngrid,nlayer,L_NSPECTI,8)  ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
+      REAL,INTENT(OUT) :: popttv(ngrid,nlayer,L_NSPECTV,8)  ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
       
       
@@ -139,18 +139,16 @@
   
       ! Optical diagnostics
-      REAL*8 popt_hazei(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_hazev(L_LEVELS,L_NSPECTV,3)
-      REAL*8 popt_cloudsi(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_cloudsv(L_LEVELS,L_NSPECTV,3)
+      ! Haze
+      REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6)
+      REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6)
+      ! Clouds
+      REAL*8 diag_optti(L_LEVELS,L_NSPECTI,6)
+      REAL*8 diag_opttv(L_LEVELS,L_NSPECTV,6)
       ! Temporary optical diagnostics (clear column)
-      REAL*8 popt_hazei_cc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_hazev_cc(L_LEVELS,L_NSPECTV,3)
-      REAL*8 popt_cloudsi_cc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_cloudsv_cc(L_LEVELS,L_NSPECTV,3)
+      REAL*8 diag_optti_cc(L_LEVELS,L_NSPECTI,6)
+      REAL*8 diag_opttv_cc(L_LEVELS,L_NSPECTV,6)
       ! Temporary optical diagnostics (dark column)
-      REAL*8 popt_hazei_dc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_hazev_dc(L_LEVELS,L_NSPECTV,3)
-      REAL*8 popt_cloudsi_dc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_cloudsv_dc(L_LEVELS,L_NSPECTV,3)
+      REAL*8 diag_optti_dc(L_LEVELS,L_NSPECTI,6)
+      REAL*8 diag_opttv_dc(L_LEVELS,L_NSPECTV,6)
 
 
@@ -508,10 +506,10 @@
          call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
               dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
-              popt_hazev_cc,popt_cloudsv_cc,cdcolumn)
+              diag_opthv,diag_opttv_cc,cdcolumn)
          ! Dark column :
          cdcolumn = 1
          call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
               dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
-              popt_hazev_dc,popt_cloudsv_dc,cdcolumn)
+              diag_opthv,diag_opttv_dc,cdcolumn)
          
          ! Mean opacity, ssa and asf :
@@ -537,40 +535,21 @@
             end do
          end do
-         wbarv = (dtauv_cc*wbarv_cc + dtauv_dc*wbarv_dc) / (dtauv_cc + dtauv_dc)
-         cosbv = (dtauv_cc*wbarv_cc*cosbv_cc + dtauv_dc*wbarv_dc*cosbv_dc) / (dtauv_cc*wbarv_cc + dtauv_dc*wbarv_dc)
+         wbarv = (1-Fcloudy) * wbarv_cc + Fcloudy * wbarv_dc
+         cosbv = (1-Fcloudy) * cosbv_cc + Fcloudy * cosbv_dc
          
-         ! Diagnostics for haze :
-         where (popt_hazev_cc(:,:,1) .lt. 1.d-30)
-            popt_hazev_cc(:,:,1) = 1.d-30
-         endwhere
-         where (popt_hazev_dc(:,:,1) .lt. 1.d-30)
-            popt_hazev_dc(:,:,1) = 1.d-30
-         endwhere
-         popt_hazev(:,:,1) = -log( (1-Fcloudy)*exp(-popt_hazev_cc(:,:,1)) + Fcloudy*exp(-popt_hazev_dc(:,:,1)) )
-         popt_hazev(:,:,2) = ((popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2)) &
-                             / (popt_hazev_cc(:,:,1) + popt_hazev_dc(:,:,1)))
-         where ((popt_hazev_cc(:,:,2).gt.0.0D0) .or. (popt_hazev_dc(:,:,2).gt.0.0D0))
-            popt_hazev(:,:,3) = (popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2)*popt_hazev_cc(:,:,3) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2)*popt_hazev_dc(:,:,3)) &
-                             / (popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2))
-         elsewhere
-            popt_hazev(:,:,3) = 0.0D0
-         endwhere
          ! Diagnostics for clouds :
          if (callclouds) then
-            where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
-               popt_cloudsv_cc(:,:,1) = 1.d-30
+            where (diag_opttv_cc(:,:,1) .lt. 1.d-30)
+               diag_opttv_cc(:,:,1) = 1.d-30
             endwhere
-            where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
-               popt_cloudsv_dc(:,:,1) = 1.d-30
+            where (diag_opttv_dc(:,:,1) .lt. 1.d-30)
+               diag_opttv_dc(:,:,1) = 1.d-30
             endwhere
-            popt_cloudsv(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsv_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsv_dc(:,:,1)) )
-            popt_cloudsv(:,:,2) = ((popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)) &
-                              / (popt_cloudsv_cc(:,:,1) + popt_cloudsv_dc(:,:,1)))
-            where ((popt_cloudsv_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsv_dc(:,:,2).gt.0.0D0))
-               popt_cloudsv(:,:,3) = (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2)*popt_cloudsv_cc(:,:,3) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)*popt_cloudsv_dc(:,:,3)) &
-                              / (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2))
-            elsewhere
-               popt_cloudsv(:,:,3) = 0.0D0
-            endwhere
+            diag_opttv(:,:,1) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,1)) + Fcloudy*exp(-diag_opttv_dc(:,:,1)) )
+            diag_opttv(:,:,2) = (1-Fcloudy) * diag_opttv_cc(:,:,2) + Fcloudy * diag_opttv_dc(:,:,2)
+            diag_opttv(:,:,3) = (1-Fcloudy) * diag_opttv_cc(:,:,3) + Fcloudy * diag_opttv_dc(:,:,3)
+            diag_opttv(:,:,4) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,4)) + Fcloudy*exp(-diag_opttv_dc(:,:,4)) )
+            diag_opttv(:,:,5) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,5)) + Fcloudy*exp(-diag_opttv_dc(:,:,5)) )
+            diag_opttv(:,:,6) = -log( (1-Fcloudy)*exp(-diag_opttv_cc(:,:,6)) + Fcloudy*exp(-diag_opttv_dc(:,:,6)) )
          endif
          
@@ -629,10 +608,10 @@
          call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
               dtaui_cc,taucumi_cc,cosbi_cc,wbari_cc,taugsurfi,seashazefact,   &
-              popt_hazei_cc,popt_cloudsi_cc,cdcolumn)
+              diag_opthi,diag_optti_cc,cdcolumn)
          ! Dark column :
          cdcolumn = 1
          call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
               dtaui_dc,taucumi_dc,cosbi_dc,wbari_dc,taugsurfi,seashazefact,   &
-              popt_hazei_dc,popt_cloudsi_dc,cdcolumn)
+              diag_opthi,diag_optti_dc,cdcolumn)
 
          ! Mean opacity, ssa and asf :
@@ -654,45 +633,21 @@
             end do
          end do
-         where ((wbari_cc.gt.0.0D0) .or. (wbari_dc.gt.0.0D0))
-            wbari = (dtaui_cc*wbari_cc + dtaui_dc*wbari_dc) / (dtaui_cc + dtaui_dc)
-            cosbi = (dtaui_cc*wbari_cc*cosbi_cc + dtaui_dc*wbari_dc*cosbi_dc) / (dtaui_cc*wbari_cc + dtaui_dc*wbari_dc)
-         elsewhere
-            wbari = 0.0D0
-            cosbi = 0.0D0
-         endwhere
-
-         ! Diagnostics for haze :
-         where (popt_hazei_cc(:,:,1) .lt. 1.d-30)
-            popt_hazei_cc(:,:,1) = 1.d-30
-         endwhere
-         where (popt_hazei_dc(:,:,1) .lt. 1.d-30)
-            popt_hazei_dc(:,:,1) = 1.d-30
-         endwhere
-         popt_hazei(:,:,1) = -log( (1-Fcloudy)*exp(-popt_hazei_cc(:,:,1)) + Fcloudy*exp(-popt_hazei_dc(:,:,1)) )
-         popt_hazei(:,:,2) = (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2)) &
-                             / (popt_hazei_cc(:,:,1) + popt_hazei_dc(:,:,1))
-         where ((popt_hazei_cc(:,:,2).gt.0.0D0) .or. (popt_hazei_dc(:,:,2).gt.0.0D0))
-            popt_hazei(:,:,3) = (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2)*popt_hazei_cc(:,:,3) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2)*popt_hazei_dc(:,:,3)) &
-                             / (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2))
-         elsewhere
-            popt_hazei(:,:,3) = 0.0D0
-         endwhere
+         wbari = (1-Fcloudy) * wbari_cc + Fcloudy * wbari_dc
+         cosbi = (1-Fcloudy) * cosbi_cc + Fcloudy * cosbi_dc
+
          ! Diagnostics for clouds :
          if (callclouds) then
-            where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
-               popt_cloudsi_cc(:,:,1) = 1.d-30
+            where (diag_optti_cc(:,:,1) .lt. 1.d-30)
+               diag_optti_cc(:,:,1) = 1.d-30
             endwhere
-            where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
-               popt_cloudsi_dc(:,:,1) = 1.d-30
+            where (diag_optti_dc(:,:,1) .lt. 1.d-30)
+               diag_optti_dc(:,:,1) = 1.d-30
             endwhere
-            popt_cloudsi(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsi_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsi_dc(:,:,1)) )
-            popt_cloudsi(:,:,2) = ((popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)) &
-                              / (popt_cloudsi_cc(:,:,1) + popt_cloudsi_dc(:,:,1)))
-            where ((popt_cloudsi_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsi_dc(:,:,2).gt.0.0D0))
-               popt_cloudsi(:,:,3) = (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2)*popt_cloudsi_cc(:,:,3) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)*popt_cloudsi_dc(:,:,3)) &
-                              / (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2))
-            elsewhere
-               popt_cloudsi(:,:,3) = 0.0D0
-            endwhere
+            diag_optti(:,:,1) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,1)) + Fcloudy*exp(-diag_optti_dc(:,:,1)) )
+            diag_optti(:,:,2) = (1-Fcloudy) * diag_optti_cc(:,:,2) + Fcloudy * diag_optti_dc(:,:,2)
+            diag_optti(:,:,3) = (1-Fcloudy) * diag_optti_cc(:,:,3) + Fcloudy * diag_optti_dc(:,:,3)
+            diag_optti(:,:,4) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,4)) + Fcloudy*exp(-diag_optti_dc(:,:,4)) )
+            diag_optti(:,:,5) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,5)) + Fcloudy*exp(-diag_optti_dc(:,:,5)) )
+            diag_optti(:,:,6) = -log( (1-Fcloudy)*exp(-diag_optti_cc(:,:,6)) + Fcloudy*exp(-diag_optti_dc(:,:,6)) )
          endif
 
@@ -779,9 +734,9 @@
          do nw = 1, L_NSPECTV
             popthv(ig,lev2lay,nw,:) = 0.0d0
-            poptcv(ig,lev2lay,nw,:) = 0.0d0
+            popttv(ig,lev2lay,nw,:) = 0.0d0
             ! Optical thickness (dtau) :
-            popthv(ig,lev2lay,nw,1) = (popt_hazev(l,nw,1) + popt_hazev(l+1,nw,1)) / 2.0
+            popthv(ig,lev2lay,nw,1) = (diag_opthv(l,nw,1) + diag_opthv(l+1,nw,1)) / 2.0
             if (callclouds) then 
-               poptcv(ig,lev2lay,nw,1) = (popt_cloudsv(l,nw,1) + popt_cloudsv(l+1,nw,1)) / 2.0
+               popttv(ig,lev2lay,nw,1) = (diag_opttv(l,nw,1) + diag_opttv(l+1,nw,1)) / 2.0
             endif
             ! Opacity (tau) :
@@ -789,5 +744,5 @@
                popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
                if (callclouds) then
-                  poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
+                  popttv(ig,lev2lay,nw,2) = popttv(ig,lev2lay,nw,2) + popttv(ig,k,nw,1)
                endif
             enddo
@@ -795,12 +750,21 @@
             popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             if (callclouds) then
-               poptcv(ig,lev2lay,nw,3) = poptcv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               popttv(ig,lev2lay,nw,3) = popttv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             endif
             ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
-            popthv(ig,lev2lay,nw,4) = (popt_hazev(l,nw,2) + popt_hazev(l+1,nw,2)) / 2.0
-            popthv(ig,lev2lay,nw,5) = (popt_hazev(l,nw,3) + popt_hazev(l+1,nw,3)) / 2.0
+            popthv(ig,lev2lay,nw,4) = (diag_opthv(l,nw,2) + diag_opthv(l+1,nw,2)) / 2.0
+            popthv(ig,lev2lay,nw,5) = (diag_opthv(l,nw,3) + diag_opthv(l+1,nw,3)) / 2.0
             if (callclouds) then
-               poptcv(ig,lev2lay,nw,4) = (popt_cloudsv(l,nw,2) + popt_cloudsv(l+1,nw,2)) / 2.0
-               poptcv(ig,lev2lay,nw,5) = (popt_cloudsv(l,nw,3) + popt_cloudsv(l+1,nw,3)) / 2.0
+               popttv(ig,lev2lay,nw,4) = (diag_opttv(l,nw,2) + diag_opttv(l+1,nw,2)) / 2.0
+               popttv(ig,lev2lay,nw,5) = (diag_opttv(l,nw,3) + diag_opttv(l+1,nw,3)) / 2.0
+            endif
+            ! DRAYAER, TAUGAS, DCONT :
+            popthv(ig,lev2lay,nw,6) = (diag_opthv(l,nw,4) + diag_opthv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            popthv(ig,lev2lay,nw,7) = (diag_opthv(l,nw,5) + diag_opthv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            popthv(ig,lev2lay,nw,8) = (diag_opthv(l,nw,6) + diag_opthv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            if (callclouds) then
+               popttv(ig,lev2lay,nw,6) = (diag_opttv(l,nw,4) + diag_opttv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               popttv(ig,lev2lay,nw,7) = (diag_opttv(l,nw,5) + diag_opttv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               popttv(ig,lev2lay,nw,8) = (diag_opttv(l,nw,6) + diag_opttv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             endif
          enddo
@@ -808,9 +772,9 @@
          do nw=1,L_NSPECTI
             popthi(ig,lev2lay,nw,:) = 0.0d0
-            poptci(ig,lev2lay,nw,:) = 0.0d0
+            poptti(ig,lev2lay,nw,:) = 0.0d0
             ! Optical thickness (dtau) :
-            popthi(ig,lev2lay,nw,1) = (popt_hazei(l,nw,1) + popt_hazei(l+1,nw,1)) / 2.0
+            popthi(ig,lev2lay,nw,1) = (diag_opthi(l,nw,1) + diag_opthi(l+1,nw,1)) / 2.0
             if (callclouds) then
-               poptci(ig,lev2lay,nw,1) = (popt_cloudsi(l,nw,1) + popt_cloudsi(l+1,nw,1)) / 2.0
+               poptti(ig,lev2lay,nw,1) = (diag_optti(l,nw,1) + diag_optti(l+1,nw,1)) / 2.0
             endif
             ! Opacity (tau) :
@@ -818,5 +782,5 @@
                popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
                if (callclouds) then
-                  poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
+                  poptti(ig,lev2lay,nw,2) = poptti(ig,lev2lay,nw,2) + poptti(ig,k,nw,1)
                endif
             enddo
@@ -824,12 +788,21 @@
             popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             if (callclouds) then
-               poptci(ig,lev2lay,nw,3) = poptci(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               poptti(ig,lev2lay,nw,3) = poptti(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             endif
             ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
-            popthi(ig,lev2lay,nw,4) = (popt_hazei(l,nw,2) + popt_hazei(l+1,nw,2)) / 2.0
-            popthi(ig,lev2lay,nw,5) = (popt_hazei(l,nw,3) + popt_hazei(l+1,nw,3)) / 2.0
+            popthi(ig,lev2lay,nw,4) = (diag_opthi(l,nw,2) + diag_opthi(l+1,nw,2)) / 2.0
+            popthi(ig,lev2lay,nw,5) = (diag_opthi(l,nw,3) + diag_opthi(l+1,nw,3)) / 2.0
             if (callclouds) then
-               poptci(ig,lev2lay,nw,4) = (popt_cloudsi(l,nw,2) + popt_cloudsi(l+1,nw,2)) / 2.0
-               poptci(ig,lev2lay,nw,5) = (popt_cloudsi(l,nw,3) + popt_cloudsi(l+1,nw,3)) / 2.0
+               poptti(ig,lev2lay,nw,4) = (diag_optti(l,nw,2) + diag_optti(l+1,nw,2)) / 2.0
+               poptti(ig,lev2lay,nw,5) = (diag_optti(l,nw,3) + diag_optti(l+1,nw,3)) / 2.0
+            endif
+            ! DRAYAER, TAUGAS, DCONT :
+            popthi(ig,lev2lay,nw,6) = (diag_opthi(l,nw,4) + diag_opthi(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            popthi(ig,lev2lay,nw,7) = (diag_opthi(l,nw,5) + diag_opthi(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            popthi(ig,lev2lay,nw,8) = (diag_opthi(l,nw,6) + diag_opthi(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+            if (callclouds) then
+               poptti(ig,lev2lay,nw,6) = (diag_optti(l,nw,4) + diag_optti(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               poptti(ig,lev2lay,nw,7) = (diag_optti(l,nw,5) + diag_optti(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
+               poptti(ig,lev2lay,nw,8) = (diag_optti(l,nw,6) + diag_optti(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
             endif
          enddo
Index: /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 3318)
@@ -54,8 +54,9 @@
 !$OMP THREADPRIVATE(opt4clouds)
       real,save    :: Fcloudy
-      real,save    :: Fssa
+      real,save    :: FCVIS
+      real,save    :: FCIR
       real,save    :: FHVIS
       real,save    :: FHIR
-!$OMP THREADPRIVATE(Fcloudy,Fssa,FHVIS,FHIR)
+!$OMP THREADPRIVATE(Fcloudy,FCVIS,FCIR,FHVIS,FHIR)
       
       logical,save :: resCH4_inf
Index: /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 3318)
@@ -118,6 +118,6 @@
       m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)                  ! then mixed, even though both are above the threshold, their ratio can be nonsense
     ELSEWHERE
-      m0as(:)=0.D0
-      m3as(:)=0.D0
+      m0as(:) = 0.D0
+      m3as(:) = 0.D0
     ENDWHERE
     WHERE (pq(ilon,:,3) > 2.D-200 .AND. pq(ilon,:,4) > 2.D-200)
@@ -125,6 +125,6 @@
       m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
     ELSEWHERE
-      m0af(:)=0.D0
-      m3af(:)=0.D0
+      m0af(:) = 0.D0
+      m3af(:) = 0.D0
     ENDWHERE
     
@@ -168,5 +168,5 @@
     ENDIF
     ! save diags (if no clouds, relevant arrays will be set to 0 !)
-    call mm_diagnostics(mmd_aer_prec(ilon),mmd_aer_s_w(ilon,:),mmd_aer_f_w(ilon,:),mmd_aer_s_flux(ilon,:),mmd_aer_f_flux(ilon,:),  &
+    call mm_diagnostics(dt,mmd_aer_prec(ilon),mmd_aer_s_w(ilon,:),mmd_aer_f_w(ilon,:),mmd_aer_s_flux(ilon,:),mmd_aer_f_flux(ilon,:),  &
                         mmd_ccn_prec(ilon),mmd_ccn_w(ilon,:),mmd_ccn_flux(ilon,:),mmd_ice_prec(ilon,:),   &
                         mmd_ice_fluxes(ilon,:,:),mmd_gazs_sat(ilon,:,:))
Index: /trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 3318)
@@ -33,4 +33,5 @@
       ! Set in physiq_mod
       character(LEN=100),save :: nudging_file ='datagcm/SuperNudging.dat'
+      ! character(LEN=100),save :: nudging_file ='datagcm/NewNudging.dat'
 !$OMP THREADPRIVATE(nudging_file)
 
Index: /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 3318)
@@ -728,5 +728,5 @@
 ! Initialize methane surface tank
 ! -------------------------------
-      tankCH4=2. ! default value for tankCH4
+      tankCH4=200. ! default value for tankCH4
       
 ! Initialize tracers
Index: /trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90	(revision 3318)
@@ -76,9 +76,10 @@
 ! Parameters :
 !-------------
-real, parameter :: karman = 0.4     ! Karman constant [-]
-real, parameter :: humCH4 = 0.5     ! Imposed surface humidity for CH4 [-]
-
-real, parameter :: Flnp = 0.05      ! Fraction occupied by lakes (North Pole)
-real, parameter :: Flsp = 0.005     ! Fraction occupied by lakes (South Pole)
+real, parameter :: karman = 0.4 ! Karman constant [-]
+real, parameter :: humCH4 = 0.4 ! Imposed surface humidity for CH4 [-]
+
+real, parameter :: Flnp = 0.08  ! Fraction occupied by lakes (North Pole)
+real, parameter :: Flsp = 0.02  ! Fraction occupied by lakes (South Pole)
+real, parameter :: Flml = 0.02  ! Fraction not infiltrated into the ground (Mid latitudes)
 
 real, parameter :: mmolair = 28.e-3 ! Molar mass of air [kg.mol-1]
@@ -94,4 +95,5 @@
 
 ! Initialisation :
+real*8  :: Tlake      ! Lake temperature [K]
 real*8  :: rhoair     ! Density of air [kg.m-3]
 real*8  :: ws         ! Horizontal wind at the surface [m.s-1]
@@ -127,6 +129,6 @@
 
   ! Saturation profile of CH4 [mol.mol-1] (Fray and Schmidt 2009) :
-  qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/tsurf(ig) - 4.341e3/tsurf(ig)**2 + 1.035e5/tsurf(ig)**3 - 7.910e5/tsurf(ig)**4)
-  qsatCH4 = humCH4 * qsatCH4
+  Tlake = tsurf(ig) - 7 ! Lakes are 2-7K less than surface.
+  qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/Tlake - 4.341e3/Tlake**2 + 1.035e5/Tlake**3 - 7.910e5/Tlake**4)
   ! CH4 : 0.80 * qsat because of dissolution in N2
   qsatCH4 = 0.80 * qsatCH4
@@ -134,14 +136,17 @@
   ! Flux at the surface [kg.m-2.s-1] :
   flux = rhoair * Cd * ws
+  
+  ! Surface humidity :
+  qsatCH4 = humCH4 * qsatCH4
 
   ! <North Pole>
-  if (REAL(latitude_deg(ig)) .ge. 70.) then
+  if (REAL(latitude_deg(ig)) .ge. 70.0) then
     flux = Flnp * flux
   ! <South Pole>
-  else if (REAL(latitude_deg(ig)) .le. -70.) then
+  else if (REAL(latitude_deg(ig)) .le. -70.0) then
     flux = Flsp * flux
   ! <Mid Latitudes>
   else
-    flux = 0.0
+    flux = Flml * flux
   endif
 
@@ -151,13 +156,12 @@
     tankCH4(ig) = 1.e-30
   endif
-  
+
+!------------------------------------
+! 2. IMPLICIT SCHEME
+!------------------------------------
+
   ! Flux of CH4 at the surface [kg.m-2.s-1.mol.mol-1] :
   fluxCH4 = flux * (qsatCH4 - pqCH4(ig,1))
 
-
-!------------------------------------
-! 2. IMPLICIT SCHEME
-!------------------------------------
-
   ! Flux at the surface [kg.m-2.s-1] --> [s-1] :
   flux = flux / rhoair / (zzlev(ig,2) - zzlev(ig,1))
@@ -170,7 +174,8 @@
 
   ! New tank depth of CH4 [m]
-  if ((tankCH4(ig) + dtankCH4) .ge. 0.) then
+  if ((tankCH4(ig) + dtankCH4) .ge. 1.e-30) then
     tankCH4(ig) = tankCH4(ig) + dtankCH4
   else
+    !write(*,*) 'Evaporation CH4 : Empty lakes...'
     newpqCH4 = tankCH4(ig) / (rhoair * (zzlev(ig,2) - zzlev(ig,1)) * mmolCH4 / mmolair / rhoiCH4) + pqCH4(ig,1)
     tankCH4(ig) = 1.e-30
@@ -185,5 +190,5 @@
 !------------------------------------
 
-  ftm = (1. - tsurf(ig) / TcCH4)
+  ftm = (1. - Tlake / TcCH4)
   if(ftm .le. 1.e-3) then
     ftm = 1.e-3
@@ -191,5 +196,5 @@
   
   ! Latente heat of CH4 vaporisation [J.kg-1.mol.mol-1]
-  LvCH4 = 8.314 * 190.4 * (7.08 * ftm**0.354 + 10.95 * 1.1e-2 * ftm**0.456) / mmolCH4
+  LvCH4 = 8.314 * TcCH4 * (7.08 * ftm**0.354 + 10.95 * 1.1e-2 * ftm**0.456) / mmolCH4
 
   ! Evaporation heating rate [K.s-1]
Index: /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 3318)
@@ -218,8 +218,13 @@
      write(*,*) "Fcloudy = ", Fcloudy
 
-     write(*,*) "Adjustment of single scattering albedo for cloudy dropplet"
-     Fssa = 1.0
-     call getin_p("Fssa", Fssa)
-     write(*,*) "Fssa = ", Fssa
+     write(*,*) "Adjustment of droplet properties in the VIS"
+     FCVIS = 1.0
+     call getin_p("FCVIS", FCVIS)
+     write(*,*) "FCVIS = ", FCVIS
+
+     write(*,*) "Adjustment of droplet properties in the IR"
+     FCIR = 1.0
+     call getin_p("FCIR", FCIR)
+     write(*,*) "FCIR = ", FCIR
      
      write(*,*) "Adjustment of aerosol properties in the VIS"
Index: /trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 3318)
@@ -1,5 +1,5 @@
 subroutine optci(PQMO,NLAY,ZLEV,PLEV,TLEV,TMID,PMID, &
      DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT,&
-     POPT_HAZE,POPT_CLOUDS,CDCOLUMN)
+     DIAG_OPTH,DIAG_OPTT,CDCOLUMN)
 
   use radinc_h
@@ -11,5 +11,5 @@
   use callkeys_mod, only: continuum,graybody,corrk_recombin,              &
                           callclouds,callmufi,seashaze,uncoupl_optic_haze,&
-                          opt4clouds,FHIR
+                          opt4clouds,FHIR,FCIR
   use tracer_h, only: nmicro,nice,ices_indx
 
@@ -61,6 +61,6 @@
   REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
   REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
-  REAL*8, INTENT(OUT) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3)
-  REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3)
+  REAL*8, INTENT(OUT) :: DIAG_OPTH(L_LEVELS,L_NSPECTI,6)
+  REAL*8, INTENT(OUT) :: DIAG_OPTT(L_LEVELS,L_NSPECTI,6)
   ! ==========================================================
   
@@ -230,4 +230,10 @@
             ENDIF
 
+            ! Tuning of optical properties for haze :
+            !dtauaer_s = dtauaer_s * (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
+            !ssa_s(nw) = ssa_s(nw) / (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
+            dtauaer_f = dtauaer_f * (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
+            ssa_f(nw) = ssa_f(nw) / (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
+
             ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
             DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
@@ -241,8 +247,9 @@
                ASF_T(k,nw)   = 1.0
             ENDIF
+            
             ! Diagnostics for the haze :
-            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
-            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
-            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
+            DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
+            DIAG_OPTH(k,nw,2) = SSA_T(k,nw)   ! wbar
+            DIAG_OPTH(k,nw,3) = ASF_T(k,nw)   ! gbar
 
             !---------------------
@@ -253,5 +260,5 @@
                m0ccn = pqmo(ilay,5) / 2.0
                m3ccn = pqmo(ilay,6) / 2.0
-               m3cld  = 0.0d0
+               m3cld  = pqmo(ilay,6) / 2.0
                
                ! Clear / Dark column method :
@@ -285,6 +292,11 @@
                
                ! For small dropplets, opacity of nucleus dominates...
+               dtau_cld = (dtau_cld*m3ccn + dtau_cld*m3cld) / (m3ccn + m3cld)
                ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
-
+               
+               ! Tuning of optical properties for clouds :
+               dtau_cld = dtau_cld * (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
+               ssa_cld(nw) = ssa_cld(nw) / (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
+               
                ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
                DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
@@ -299,18 +311,14 @@
                ENDIF
 
-               ! Tuning of optical properties (now useless...) :
-               DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
-               SSA_T(k,nw)   = SSA_T(k,nw) / (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
-
                ! Diagnostics for clouds :
-               POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau
-               POPT_CLOUDS(k,nw,2) = SSA_T(k,nw)   ! wbar
-               POPT_CLOUDS(k,nw,3) = ASF_T(k,nw)   ! gbar
+               DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
+               DIAG_OPTT(k,nw,2) = SSA_T(k,nw)   ! wbar
+               DIAG_OPTT(k,nw,3) = ASF_T(k,nw)   ! gbar
 
             ELSE
                ! Diagnostics for clouds :
-               POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
-               POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
-               POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
+               DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
+               DIAG_OPTT(k,nw,2) = 1.0  ! wbar
+               DIAG_OPTT(k,nw,3) = 1.0  ! gbar
             ENDIF
             
@@ -321,11 +329,11 @@
             if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
             ! Diagnostics for the haze :
-            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
-            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
-            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
+            DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
+            DIAG_OPTH(k,nw,2) = SSA_T(k,nw)   ! wbar
+            DIAG_OPTH(k,nw,3) = ASF_T(k,nw)   ! gbar
             ! Diagnostics for clouds :
-            POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
-            POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
-            POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
+            DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
+            DIAG_OPTT(k,nw,2) = 1.0  ! wbar
+            DIAG_OPTT(k,nw,3) = 1.0  ! gbar
          ENDIF ! ENDIF callmufi
          
@@ -436,4 +444,11 @@
                            + DCONT   & ! For parameterized continuum absorption
                         + DHAZE_T(K,NW)     ! For Titan Haze
+
+         DIAG_OPTH(K,nw,4) = 0.d0
+         DIAG_OPTH(K,nw,5) = TAUGAS
+         DIAG_OPTH(K,nw,6) = DCONT
+         DIAG_OPTT(K,nw,4) = 0.d0
+         DIAG_OPTT(K,nw,5) = TAUGAS
+         DIAG_OPTT(K,nw,6) = DCONT
 
       end do ! k = L_LEVELS
Index: /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 3318)
@@ -1,5 +1,5 @@
 SUBROUTINE OPTCV(PQMO,NLAY,ZLEV,PLEV,TMID,PMID,  &
      DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT,&
-     POPT_HAZE,POPT_CLOUDS,CDCOLUMN)
+     DIAG_OPTH,DIAG_OPTT,CDCOLUMN)
 
   use radinc_h
@@ -11,5 +11,5 @@
   use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,   &
                           callclouds,callmufi,seashaze,uncoupl_optic_haze,&
-                          opt4clouds,FHVIS, Fssa
+                          opt4clouds,FHVIS,FCVIS
   use tracer_h, only: nmicro,nice,ices_indx
 
@@ -69,6 +69,6 @@
   REAL*8, INTENT(OUT) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
   REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTV,L_NGAUSS-1)
-  REAL*8, INTENT(OUT) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3)
-  REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3)
+  REAL*8, INTENT(OUT) :: DIAG_OPTH(L_LEVELS,L_NSPECTV,6)
+  REAL*8, INTENT(OUT) :: DIAG_OPTT(L_LEVELS,L_NSPECTV,6)
   ! ==========================================================
   
@@ -124,5 +124,5 @@
   real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
   real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV)
-  real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
+  real*8,save :: ssa_ccn(L_NSPECTV),asf_ccn(L_NSPECTV)
   real*8,save :: ssa_cld(L_NSPECTV),asf_cld(L_NSPECTV)
 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
@@ -191,4 +191,7 @@
    end do
   
+   DIAG_OPTH(:,:,:) = 0.D0
+   DIAG_OPTT(:,:,:) = 0.D0
+
    do NW=1,L_NSPECTV
       ! We ignore K=1...
@@ -254,4 +257,10 @@
             ENDIF
 
+            ! Tuning of optical properties for haze :
+            !dtauaer_s = dtauaer_s * (FHVIS * (1-ssa_s(nw)) + ssa_s(nw))
+            !ssa_s(nw) = ssa_s(nw) / (FHVIS * (1-ssa_s(nw)) + ssa_s(nw))
+            dtauaer_f = dtauaer_f * (FHVIS * (1-ssa_f(nw)) + ssa_f(nw))
+            ssa_f(nw) = ssa_f(nw) / (FHVIS * (1-ssa_f(nw)) + ssa_f(nw))
+
             ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
             DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
@@ -265,8 +274,9 @@
                ASF_T(k,nw)   = 1.0
             ENDIF
+
             ! Diagnostics for the haze :
-            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
-            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
-            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
+            DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
+            DIAG_OPTH(k,nw,2) = SSA_T(k,nw)   ! wbar
+            DIAG_OPTH(k,nw,3) = ASF_T(k,nw)   ! gbar
 
             !---------------------
@@ -277,5 +287,5 @@
                m0ccn = pqmo(ilay,5) / 2.0
                m3ccn = pqmo(ilay,6) / 2.0
-               m3cld = 0.0d0
+               m3cld  = pqmo(ilay,6) / 2.0
                
                ! Clear / Dark column method :
@@ -309,7 +319,11 @@
 
                ! For small dropplets, opacity of nucleus dominates
+               dtau_cld = (dtau_cld*m3ccn + dtau_cld*m3cld) / (m3ccn + m3cld)
                ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
-               ssa_cld(nw) = ssa_cld(nw) * Fssa
-
+
+               ! Tuning of optical properties for clouds :
+               dtau_cld = dtau_cld * (FCVIS * (1-ssa_cld(nw)) + ssa_cld(nw))
+               ssa_cld(nw) = ssa_cld(nw) / (FCVIS * (1-ssa_cld(nw)) + ssa_cld(nw))
+               
                ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
                DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
@@ -324,18 +338,14 @@
                ENDIF
                
-               ! Tuning of optical properties (now useless...) :
-               DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw))
-               SSA_T(k,nw)   = SSA_T(k,nw) / (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw))
-
                ! Diagnostics for clouds :
-               POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau
-               POPT_CLOUDS(k,nw,2) = SSA_T(k,nw)   ! wbar
-               POPT_CLOUDS(k,nw,3) = ASF_T(k,nw)   ! gbar
+               DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
+               DIAG_OPTT(k,nw,2) = SSA_T(k,nw)   ! wbar
+               DIAG_OPTT(k,nw,3) = ASF_T(k,nw)   ! gbar
             
             ELSE
                ! Diagnostics for clouds :
-               POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
-               POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
-               POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
+               DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
+               DIAG_OPTT(k,nw,2) = 1.0  ! wbar
+               DIAG_OPTT(k,nw,3) = 1.0  ! gbar
             ENDIF 
             
@@ -346,11 +356,11 @@
             if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
             ! Diagnostics for the haze :
-            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
-            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
-            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
+            DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
+            DIAG_OPTH(k,nw,2) = SSA_T(k,nw)   ! wbar
+            DIAG_OPTH(k,nw,3) = ASF_T(k,nw)   ! gbar
             ! Diagnostics for clouds :
-            POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
-            POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
-            POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
+            DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
+            DIAG_OPTT(k,nw,2) = 1.0  ! wbar
+            DIAG_OPTT(k,nw,3) = 1.0  ! gbar
          ENDIF ! ENDIF callmufi
          
@@ -472,6 +482,13 @@
          DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze
 
-      end do ! k = L_LEVELS
-   end do ! nw = L_NSPECTI
+         DIAG_OPTH(K,nw,4) = DRAYAER
+         DIAG_OPTH(K,nw,5) = TAUGAS
+         DIAG_OPTH(K,nw,6) = DCONT
+         DIAG_OPTT(K,nw,4) = DRAYAER
+         DIAG_OPTT(K,nw,5) = TAUGAS
+         DIAG_OPTT(K,nw,6) = DCONT
+
+      end do ! K = L_LEVELS
+   end do ! nw = L_NSPECTV
 
   !=======================================================================
Index: /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 3318)
@@ -20,4 +20,5 @@
                      inquire_dimension, inquire_dimension_length
   use callkeys_mod, only: surfalbedo,surfemis
+  use geometry_mod, only: latitude_deg
 
   implicit none
@@ -219,5 +220,5 @@
              minval(emis), maxval(emis)
 
-! Depth of methane tank (added by JVO 2017)
+! Depth of methane tank (added by BdBdT 2023)
  if (startphy_file) then
    call get_field("tankCH4",tankCH4,found,indextime)
@@ -225,8 +226,24 @@
      write(*,*) "phyetat0: Failed loading <tankCH4>"
      !  call abort
-     tankCH4(:)=2.
+     DO ig = 1, ngrid
+      if (REAL(latitude_deg(ig)) .ge. 70.0) then
+        tankCH4(ig) = 200.0 ! [m]
+      else if (REAL(latitude_deg(ig)) .le. -70.0) then
+        tankCH4(ig) = 50.0  ! [m]
+      else
+        tankCH4(ig) = 0.0   ! [m]
+      endif
+     ENDDO
    endif
  else
-   tankCH4(:)=2.
+  DO ig = 1, ngrid
+    if (REAL(latitude_deg(ig)) .ge. 70.0) then
+      tankCH4(ig) = 200.0 ! [m]
+    else if (REAL(latitude_deg(ig)) .le. -70.0) then
+      tankCH4(ig) = 50.0  ! [m]
+    else
+      tankCH4(ig) = 0.0   ! [m]
+    endif
+   ENDDO
  endif ! of if (startphy_file)
  write(*,*) "phyetat0: Depth of methane tank <tankCH4> range:", &
Index: /trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90	(revision 3318)
@@ -72,9 +72,9 @@
 !$OMP THREADPRIVATE(int_dtaui,int_dtauv) 
 
-      real,dimension(:,:,:,:),allocatable,save :: zpopthi ! IR optical properties of haze within narrowbands for diags (dtau,tau,k,wbar,gbar).
-      real,dimension(:,:,:,:),allocatable,save :: zpopthv ! VI optical properties of haze within narrowbands for diags (dtau,tau,k,wbar,gbar).
-      real,dimension(:,:,:,:),allocatable,save :: zpoptci ! IR optical properties of clouds within narrowbands for diags (dtau,tau,k,wbar,gbar).
-      real,dimension(:,:,:,:),allocatable,save :: zpoptcv ! VI optical properties of clouds within narrowbands for diags (dtau,tau,k,wbar,gbar).
-!$OMP THREADPRIVATE(zpopthi,zpopthv,zpoptci,zpoptcv)
+      real,dimension(:,:,:,:),allocatable,save :: zpopthi ! IR optical properties [haze] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont).
+      real,dimension(:,:,:,:),allocatable,save :: zpopthv ! VI optical properties [haze] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont).
+      real,dimension(:,:,:,:),allocatable,save :: zpoptti ! IR optical properties [haze+clouds] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont).
+      real,dimension(:,:,:,:),allocatable,save :: zpopttv ! VI optical properties [haze+clouds] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont).
+!$OMP THREADPRIVATE(zpopthi,zpopthv,zpoptti,zpopttv)
 
       real,dimension(:,:),allocatable,save :: u_ref ! Reference profile for the zonal wind nudging (m.s-1).
@@ -139,8 +139,8 @@
         ALLOCATE(int_dtaui(klon,klev,L_NSPECTI))
         ALLOCATE(int_dtauv(klon,klev,L_NSPECTV))
-        ALLOCATE(zpopthi(klon,klev,L_NSPECTI,5))
-        ALLOCATE(zpopthv(klon,klev,L_NSPECTV,5))
-        ALLOCATE(zpoptci(klon,klev,L_NSPECTI,5))
-        ALLOCATE(zpoptcv(klon,klev,L_NSPECTV,5))
+        ALLOCATE(zpopthi(klon,klev,L_NSPECTI,8))
+        ALLOCATE(zpopthv(klon,klev,L_NSPECTV,8))
+        ALLOCATE(zpoptti(klon,klev,L_NSPECTI,8))
+        ALLOCATE(zpopttv(klon,klev,L_NSPECTV,8))
         ALLOCATE(u_ref(49,24))
         allocate(dycchi(klon,klev,nkim))
@@ -207,6 +207,6 @@
         DEALLOCATE(zpopthi)
         DEALLOCATE(zpopthv)
-        DEALLOCATE(zpoptci)
-        DEALLOCATE(zpoptcv)
+        DEALLOCATE(zpoptti)
+        DEALLOCATE(zpopttv)
         DEALLOCATE(u_ref)
         DEALLOCATE(dycchi)
Index: /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 3318)
@@ -15,5 +15,5 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV
-      use radcommon_h, only: sigma, gzlat, grav, BWNV, WNOI, WNOV
+      use radcommon_h, only: sigma, gzlat, grav, BWNV, WAVEI, WAVEV
       use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
@@ -161,8 +161,6 @@
 !                + clean of all too-generic (ocean, water, co2 ...) routines
 !                + Titan's chemistry
-!           Microphysical moment model - J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
-!           Modified : B. de Batz de Trenquelléon (2023)
-!                + Optical improvements (haze and cloud)
-!                + Nudging of troposphere !!
+!           Microphysical moment model: J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
+!           Optics for haze and clouds: B. de Batz de Trenquelléon (2023)
 !============================================================================================
 
@@ -279,4 +277,6 @@
       real zdpsrfmr(ngrid)        ! Pressure tendency for mass_redistribution routine (Pa/s).
 
+      ! Temporal tracers :
+      real ppq(ngrid,nlayer,nq)   ! Tracers (kg/kg_of_air).
       
     
@@ -422,6 +422,6 @@
          zpopthi(:,:,:,:) = 0.D0
          zpopthv(:,:,:,:) = 0.D0
-         zpoptci(:,:,:,:) = 0.D0
-         zpoptcv(:,:,:,:) = 0.D0
+         zpoptti(:,:,:,:) = 0.D0
+         zpopttv(:,:,:,:) = 0.D0
 
 !        Initialize setup for correlated-k radiative transfer
@@ -589,8 +589,8 @@
          endif
 
-         ! Read SuperNudging.dat and initialize the nudging for pu
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         ! Read NewNudging.dat and initialize the nudging for pu
+         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          if (nudging_u) then
-            nudging_file=trim(datadir)//'/SuperNudging.dat'
+            nudging_file=trim(datadir)//'/NewNudging.dat'
             inquire(file=trim(nudging_file),exist=file_ok)
             if(.not.file_ok) then
@@ -620,5 +620,5 @@
          write(*,*) "physiq: call initialize_xios_output"
          call initialize_xios_output(pday,ptime,ptimestep,daysec, &
-                                     presnivs,pseudoalt,wnoi,wnov)
+                                     presnivs,pseudoalt,wavei,wavev)
 #endif
          write(*,*) "physiq: end of firstcall"
@@ -840,4 +840,5 @@
 
                ! standard callcorrk
+               !call callcorrk(ngrid,nlayer,ppq,nq,qsurf,zday,                      &
                call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                      &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
@@ -846,5 +847,5 @@
                               fluxsurfabs_sw,fluxtop_lw,                          &
                               fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
-                              int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptci,zpoptcv,&
+                              int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,&
                               lastcall)
 
@@ -932,5 +933,4 @@
 !  III. Vertical diffusion (turbulent mixing) :
 !  --------------------------------------------
-
       if (calldifv) then
       
@@ -1233,5 +1233,5 @@
             ! iv. Surface pseudo-evaporation
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            ! Infinite tank of CH4 (global)
+            ! Infinite tank of CH4
             if (resCH4_inf) then
                do ig=1,ngrid
@@ -1243,7 +1243,9 @@
                enddo
 
-            ! Tank at the pole (for the moment infinit...)
             else
-               tankCH4(:) = 2.
+               ! Fill lakes with precipitation :
+               tankCH4(:) = tankCH4(:) + (mmd_ice_prec(:,1) / 425. * ptimestep) ! [m]
+               
+               ! Evaporation of lakes :
                if (moyzon_ch) then
                  tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
@@ -1256,5 +1258,6 @@
             endif
 
-            pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro)
+
+            pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro) ! convert tendencies to mass mixing ratio
             
             ! v. Updates and positivity check
@@ -1391,10 +1394,74 @@
       enddo
 
-      ! Tracers : for now because problem with re-evaporation in microphysics ...
-      if (callclouds) then
-         pdq(:,1:13,gazs_indx(2)) = (2e-8 - pq(:,1:13,gazs_indx(2))) / (ptimestep * 100)  ! C2H6
-         pdq(:,1:12,gazs_indx(3)) = (1e-9 - pq(:,1:12,gazs_indx(3))) / (ptimestep * 100)  ! C2H2
-         pdq(:,1:11,gazs_indx(4)) = (1e-20 - pq(:,1:11,gazs_indx(4))) / (ptimestep * 100) ! HCN
-      endif
+      ! [BBT TEMPO] Firstcall : Initialize Tracers
+      !if (firstcall) then
+      !   do ig = 1, ngrid
+      !      ! Aerosols :
+      !      !-----------
+      !      WHERE (pq(ig,:24,3) < 3.0e7) &
+      !               pdq(ig,:24,3) = (3.0e7 - pq(ig,:24,3)) / ptimestep   ! M0 aerf
+      !      WHERE (pq(ig,:24,4) < 1.5e-12) &
+      !               pdq(ig,:24,4) = (1.5e-12 - pq(ig,:24,4)) / ptimestep ! M3 aerf
+      !      ! Lakes :
+      !      !--------
+      !      if (REAL(latitude_deg(ig)) .ge. 70.0) then
+      !         tankCH4(ig) = 100.0 ! [m]
+      !      else if (REAL(latitude_deg(ig)) .le. -70.0) then
+      !         tankCH4(ig) = 20.0  ! [m]
+      !      else
+      !         tankCH4(ig) = 0.0   ! [m]
+      !      endif
+      !      ! Species :
+      !      !----------
+      !      do iq = 1, size(ices_indx)
+      !         ! CH4 :
+      !         !------
+      !         if(trim(nameOfTracer(gazs_indx(iq))) .eq. "CH4") then
+      !            WHERE (pq(ig,:,gazs_indx(iq))/rat_mmol(gazs_indx(iq)) > 0.05) &
+      !               pdq(ig,:,gazs_indx(iq)) = (0.05 * rat_mmol(gazs_indx(iq)) - pq(ig,:,gazs_indx(iq))) / ptimestep
+      !            WHERE (pq(ig,:,gazs_indx(iq))/rat_mmol(gazs_indx(iq)) < 0.014) &
+      !               pdq(ig,:,gazs_indx(iq)) = (0.014 * rat_mmol(gazs_indx(iq)) - pq(ig,:,gazs_indx(iq))) / ptimestep
+      !         endif
+      !         ! C2H2 :
+      !         !-------
+      !         if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H2") then
+      !            WHERE (pq(ig,26:,gazs_indx(iq))/rat_mmol(gazs_indx(iq)) < 3.0e-6) &
+      !               pdq(ig,26:,gazs_indx(iq)) = (5.0e-6 * rat_mmol(gazs_indx(iq)) - pq(ig,26:,gazs_indx(iq))) / ptimestep
+      !         endif
+      !         ! C2H6 :
+      !         !-------
+      !         if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H6") then
+      !            WHERE (pq(ig,26:,gazs_indx(iq))/rat_mmol(gazs_indx(iq)) < 2.0e-5) &
+      !               pdq(ig,26:,gazs_indx(iq)) = (2.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,26:,gazs_indx(iq))) / ptimestep
+      !         endif
+      !         ! HCN :
+      !         !------
+      !         if(trim(nameOfTracer(gazs_indx(iq))) .eq. "HCN") then
+      !            WHERE (pq(ig,26:,gazs_indx(iq))/rat_mmol(gazs_indx(iq)) < 8.0e-6) &
+      !               pdq(ig,26:,gazs_indx(iq)) = (8.0e-6 * rat_mmol(gazs_indx(iq)) - pq(ig,26:,gazs_indx(iq))) / ptimestep
+      !         endif
+      !      enddo
+      !   enddo
+      !endif
+      ! [BdBdT : Forcage de la photochimie]
+      do ig = 1, ngrid
+         do iq = 1, size(ices_indx)
+            ! C2H2 :
+            !-------
+            if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H2") then
+               pdq(ig,nlayer-3:,gazs_indx(iq)) = (4.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
+            endif
+            ! C2H6 :
+            !-------
+            if(trim(nameOfTracer(gazs_indx(iq))) .eq. "C2H6") then
+               pdq(ig,nlayer-3:,gazs_indx(iq)) = (1.0e-4 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
+            endif
+            ! HCN :
+            !------
+            if(trim(nameOfTracer(gazs_indx(iq))) .eq. "HCN") then
+               pdq(ig,nlayer-3:,gazs_indx(iq)) = (4.0e-5 * rat_mmol(gazs_indx(iq)) - pq(ig,nlayer-3:,gazs_indx(iq))) / ptimestep
+            endif
+         enddo
+      enddo
       zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
 
@@ -1684,9 +1751,4 @@
       ENDIF
 
-      ! Subsurface temperatures (3D)
-      !CALL send_xios_field("tempsoil",tsoil)
-      ! Total energy balance diagnostics (2D)
-      !CALL send_xios_field("ALB",albedo_equivalent)
-
       !--------------------------------------------------------
       ! Winds trends :
@@ -1747,6 +1809,73 @@
       ! Optical diagnostics :
       !--------------------------------------------------------
-      ! Diagnostics for haze :
-      ! IR
+      ! Haze opacity :
+      CALL send_xios_field('ttauhv_14',zpopthv(:,:,14,2)) ! 14 --> 1.076 um
+      CALL send_xios_field('ttauhv_20',zpopthv(:,:,20,2)) ! 20 --> 0.671 um
+      CALL send_xios_field('ttauhv_23',zpopthv(:,:,23,2)) ! 23 --> 0.346 um
+      CALL send_xios_field('ttauhi_02',zpopthi(:,:,2,2))  ! 02 --> 175.3 um
+      CALL send_xios_field('ttauhi_17',zpopthi(:,:,17,2)) ! 17 --> 11.00 um
+      CALL send_xios_field('ttauhi_23',zpopthi(:,:,23,2)) ! 23 --> 4.849 um
+      ! Haze extinction :
+      CALL send_xios_field('kkhv_14',zpopthv(:,:,14,3))
+      CALL send_xios_field('kkhv_20',zpopthv(:,:,20,3))
+      CALL send_xios_field('kkhv_23',zpopthv(:,:,23,3))
+      CALL send_xios_field('kkhi_02',zpopthi(:,:,2,3))
+      CALL send_xios_field('kkhi_17',zpopthi(:,:,17,3))
+      CALL send_xios_field('kkhi_23',zpopthi(:,:,23,3))
+      ! Haze single scattering albedo :
+      CALL send_xios_field('wwhv_14',zpopthv(:,:,14,4))
+      CALL send_xios_field('wwhv_20',zpopthv(:,:,20,4))
+      CALL send_xios_field('wwhv_23',zpopthv(:,:,23,4))
+      CALL send_xios_field('wwhi_02',zpopthi(:,:,2,4))
+      CALL send_xios_field('wwhi_17',zpopthi(:,:,17,4))
+      CALL send_xios_field('wwhi_23',zpopthi(:,:,23,4))
+      ! Haze asymmetry parameter :
+      CALL send_xios_field('gghv_14',zpopthv(:,:,14,5))
+      CALL send_xios_field('gghv_20',zpopthv(:,:,20,5))
+      CALL send_xios_field('gghv_23',zpopthv(:,:,23,5))
+      CALL send_xios_field('gghi_02',zpopthi(:,:,2,5))
+      CALL send_xios_field('gghi_17',zpopthi(:,:,17,5))
+      CALL send_xios_field('gghi_23',zpopthi(:,:,23,5))
+
+      ! Diagnostics for haze and clouds :
+      IF (callclouds) THEN
+         ! Opacity :
+         CALL send_xios_field('ttauv_14',zpopttv(:,:,14,2)) ! 14 --> 1.076 um
+         CALL send_xios_field('ttauv_20',zpopttv(:,:,20,2)) ! 20 --> 0.671 um
+         CALL send_xios_field('ttauv_23',zpopttv(:,:,23,2)) ! 23 --> 0.346 um
+         CALL send_xios_field('ttaui_02',zpoptti(:,:,2,2))  ! 02 --> 175.3 um
+         CALL send_xios_field('ttaui_17',zpoptti(:,:,17,2)) ! 17 --> 11.00 um
+         CALL send_xios_field('ttaui_23',zpoptti(:,:,23,2)) ! 23 --> 4.849 um
+         ! Extinction :
+         CALL send_xios_field('kkv_14',zpopttv(:,:,14,3))
+         CALL send_xios_field('kkv_20',zpopttv(:,:,20,3))
+         CALL send_xios_field('kkv_23',zpopttv(:,:,23,3))
+         CALL send_xios_field('kki_02',zpoptti(:,:,2,3))
+         CALL send_xios_field('kki_17',zpoptti(:,:,17,3))
+         CALL send_xios_field('kki_23',zpoptti(:,:,23,3))
+         ! Single scattering albedo :
+         CALL send_xios_field('wwv_14',zpopttv(:,:,14,4))
+         CALL send_xios_field('wwv_20',zpopttv(:,:,20,4))
+         CALL send_xios_field('wwv_23',zpopttv(:,:,23,4))
+         CALL send_xios_field('wwi_02',zpoptti(:,:,2,4))
+         CALL send_xios_field('wwi_17',zpoptti(:,:,17,4))
+         CALL send_xios_field('wwi_23',zpoptti(:,:,23,4))
+         ! Asymmetry parameter :
+         CALL send_xios_field('ggv_14',zpopttv(:,:,14,5))
+         CALL send_xios_field('ggv_20',zpopttv(:,:,20,5))
+         CALL send_xios_field('ggv_23',zpopttv(:,:,23,5))
+         CALL send_xios_field('ggi_02',zpoptti(:,:,2,5))
+         CALL send_xios_field('ggi_17',zpoptti(:,:,17,5))
+         CALL send_xios_field('ggi_23',zpoptti(:,:,23,5))
+         ! DRAYAER, TAUGAS, DCONT :
+         CALL send_xios_field('drayaerv_20',zpopttv(:,:,20,6)) ! 20 --> 0.671um
+         CALL send_xios_field('taugasv_20',zpopttv(:,:,20,7))
+         CALL send_xios_field('dcontv_20',zpopttv(:,:,20,8))
+         CALL send_xios_field('drayaeri_17',zpoptti(:,:,17,6)) ! 17 --> 11.00um
+         CALL send_xios_field('taugasi_17',zpoptti(:,:,17,7))
+         CALL send_xios_field('dconti_17',zpoptti(:,:,17,8))
+      ENDIF      
+
+      ! Diagnostics for haze and clouds (4D) :
       CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
       CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
@@ -1754,5 +1883,4 @@
       CALL send_xios_field('whi',zpopthi(:,:,:,4))
       CALL send_xios_field('ghi',zpopthi(:,:,:,5))
-      ! VI
       CALL send_xios_field('dtauhv',zpopthv(:,:,:,1))
       CALL send_xios_field('tauhv',zpopthv(:,:,:,2))
@@ -1760,20 +1888,16 @@
       CALL send_xios_field('whv',zpopthv(:,:,:,4))
       CALL send_xios_field('ghv',zpopthv(:,:,:,5))
-
-      ! Diagnostics for haze and clouds :
       IF (callclouds) THEN
-         ! IR
-         CALL send_xios_field('dtaui',zpoptci(:,:,:,1))
-         CALL send_xios_field('taui',zpoptci(:,:,:,2))
-         CALL send_xios_field('ki',zpoptci(:,:,:,3))
-         CALL send_xios_field('wi',zpoptci(:,:,:,4))
-         CALL send_xios_field('gi',zpoptci(:,:,:,5))
-         ! VI
-         CALL send_xios_field('dtauv',zpoptcv(:,:,:,1))
-         CALL send_xios_field('tauv',zpoptcv(:,:,:,2))
-         CALL send_xios_field('kv',zpoptcv(:,:,:,3))
-         CALL send_xios_field('wv',zpoptcv(:,:,:,4))
-         CALL send_xios_field('gv',zpoptcv(:,:,:,5))
-      ENDIF      
+         CALL send_xios_field('dtaui',zpoptti(:,:,:,1))
+         CALL send_xios_field('taui',zpoptti(:,:,:,2))
+         CALL send_xios_field('ki',zpoptti(:,:,:,3))
+         CALL send_xios_field('wi',zpoptti(:,:,:,4))
+         CALL send_xios_field('gi',zpoptti(:,:,:,5))
+         CALL send_xios_field('dtauv',zpopttv(:,:,:,1))
+         CALL send_xios_field('tauv',zpopttv(:,:,:,2))
+         CALL send_xios_field('kv',zpopttv(:,:,:,3))
+         CALL send_xios_field('wv',zpopttv(:,:,:,4))
+         CALL send_xios_field('gv',zpopttv(:,:,:,5))
+      ENDIF
 
       !--------------------------------------------------------
@@ -1826,4 +1950,8 @@
       !--------------------------------------------------------
       IF (callchim) THEN
+         ! Surface (2D) :
+         CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
+         CALL send_xios_field("tankCH4",tankCH4(:))    ! CH4 tank at the surface (m)
+
          ! Atmosphere (3D) :
          ! Chemical species :
@@ -1877,8 +2005,4 @@
             CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
          ENDDO 
-
-         ! Surface (2D) :
-         CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
-         CALL send_xios_field("tankCH4",tankCH4(:))    ! CH4 tank at the surface (m)
       ENDIF ! of 'if callchim'
 
Index: /trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90	(revision 3317)
+++ /trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90	(revision 3318)
@@ -18,5 +18,5 @@
 
   SUBROUTINE initialize_xios_output(day,timeofday,dtphys,daysec,&
-                                    presnivs,pseudoalt,wnoi,wnov)
+                                    presnivs,pseudoalt,wavei,wavev)
 !  USE mod_phys_lmdz_para, only: gather, bcast, &
 !                                jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
@@ -40,12 +40,12 @@
   IMPLICIT NONE
   
-  REAL,INTENT(IN) :: day ! Number of elapsed sols since reference Ls=0.
-  REAL,INTENT(IN) :: timeofday ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
-  REAL,INTENT(IN) :: dtphys ! physics time step (s)
-  REAL,INTENT(IN) :: daysec ! lengthof a standard day (s)
-  REAL,INTENT(IN) :: presnivs(:) ! vertical grid approximate pressure (Pa)
+  REAL,INTENT(IN) :: day          ! Number of elapsed sols since reference Ls=0.
+  REAL,INTENT(IN) :: timeofday    ! "Universal time", given as fraction of sol (e.g.: 0.5 for noon).
+  REAL,INTENT(IN) :: dtphys       ! physics time step (s)
+  REAL,INTENT(IN) :: daysec       ! lengthof a standard day (s)
+  REAL,INTENT(IN) :: presnivs(:)  ! vertical grid approximate pressure (Pa)
   REAL,INTENT(IN) :: pseudoalt(:) ! vertical grid approximate altitude (km)
-  REAL,INTENT(IN) :: wnoi(:) ! Array of wavelength channels for the infrared.
-  REAL,INTENT(IN) :: wnov(:) ! Array of wavelength channels for the visible.
+  REAL,INTENT(IN) :: wavei(:)     ! Array of the wavelenght (in microns) at the center of each IR spectral interval
+  REAL,INTENT(IN) :: wavev(:)     ! Array of the wavelenght (in microns) at the center of each VI spectral interval
   
   INTEGER :: i
@@ -74,9 +74,9 @@
                             unit="Pa",positive="down")
     
-    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for IR_Wavelength"
-    CALL xios_set_axis_attr("IR_Wavelength", n_glo=size(wnoi), value=wnoi,unit="m",positive="up")
-
-    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for VI_Wavelength"
-    CALL xios_set_axis_attr("VI_Wavelength", n_glo=size(wnov), value=wnov,unit="m",positive="up")
+    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for wavelength_ir"
+    CALL xios_set_axis_attr("wavelength_ir", n_glo=size(wavei), value=wavei,unit="micrometer",positive="down")
+
+    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for wavelength_vi"
+    CALL xios_set_axis_attr("wavelength_vi", n_glo=size(wavev), value=wavev,unit="micrometer",positive="down")
                   
     ! Calculation of pseudo-altitudes is still to be done
