Index: trunk/LMDZ.TITAN/libf/phytitan/calc_condlh.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calc_condlh.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/calc_condlh.F90	(revision 3090)
@@ -21,6 +21,6 @@
 !   Lcx = Condensation latente heat [J.kg-1]
 !
-!   Author(s)
-!   ---------
+!   Author
+!   ------
 !   B. de Batz de Trenquelléon (07/2023)
 !
@@ -28,8 +28,9 @@
 !   ---------------
 !   We suppose a given order of tracers !
+!   Tracers come from microphysics (ices_indx) :
 !   CH4  --> 1
 !   C2H2 --> 2
 !   C2H6 --> 3
-!   HCN  --> 4 (TO DO)
+!   HCN  --> 4
 !====================================================================    
 
Index: trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90	(revision 3090)
@@ -7,6 +7,6 @@
 !     Average the Rayleigh scattering in each band, weighting the 
 !     average by the blackbody function at temperature tstellar.
-!     Works for an arbitrary mix of gases (currently CO2, N2 and 
-!     H2 are possible).
+!     Works for an arbitrary mix of gases (currently N2, H2
+!     and CH4 are possible).
 !     
 !     Authors
@@ -64,7 +64,7 @@
       typeknown=.false.
       do igas=1,ngasmx
-	 if(gfrac(igas,nivref).lt.1.e-2)then
+         if(gfrac(igas,nivref).lt.1.e-2)then
             print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
-            'as its mixing ratio is less than 0.01.' 
+               'as its mixing ratio is less than 0.01.' 
             ! ignore variable gas in Rayleigh calculation
             ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
@@ -78,13 +78,13 @@
                ! uses n=1.000132 from Optics, Fourth Edition.
                ! was out by a factor 2!
-	    elseif(igas.eq.igas_CH4)then
-                ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
-                ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
-                ! Values are taken from Caldas (2019), equation 12 & appendix D
-                ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
-                ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
-                ! 1E24 = (1E6)**4 because wavelength is in micron
-                ! 16.04*1.6726*1E-27 is CH4 molecular mass
-                tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
+            elseif(igas.eq.igas_CH4)then
+               ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
+               ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
+               ! Values are taken from Caldas (2019), equation 12 & appendix D
+               ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
+               ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
+               ! 1E24 = (1E6)**4 because wavelength is in micron
+               ! 16.04*1.6726*1E-27 is CH4 molecular mass
+               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
             else
                print*,'Error in calc_rayleigh: Gas species not recognised!'
Index: trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90	(revision 3090)
@@ -12,5 +12,9 @@
   !     J. Vatant d'Ollone
   !       -> 2017 : Adapt to new physics, move to F90 and compute every grid point 
-  !       -> 2018 : Update saturation profiles from recent litterature ( Vuitton 2019 ) 
+  !       -> 2018 : Update saturation profiles from recent litterature (Vuitton 2019)
+  !     Modified
+  !     --------
+  !     B. de Batz de Trenquelléon
+  !       -> 2022 : Update saturation profiles and correct the pressure unit
   !     ==============================================================================
 
@@ -40,146 +44,167 @@
   ALLOCATE(x(nlon,nlay))
 
-  !     By default, ysat=1, meaning we do not condense
+  ! By default, ysat = 1, meaning we do not condense
   ysat(:,:,:) = 1.0
 
-  !     Main loop
-
-  do ic=1,nkim
-
-     if(trim(cnames(ic)).eq."CH4") then
-
-        !where (temp(:,:).lt.90.65)                                      
-        !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                  &
-        !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                    &
-        !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
-        !elsewhere
-        !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 )   &
-        !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
-        !endwhere
-
-        ! Fray and Schmidt (2009)
-        ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)              &
-                       - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
-
-        ! Forcing CH4 to 1.4% minimum               
-        where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
-
-     else if(trim(cnames(ic)).eq."C2H2") then
-
-!        ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                        &
-!             * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
-
-         ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:)) ! Fray and Schmidt (2009)
-
-     else if(trim(cnames(ic)).eq."C2H4") then
-
-        where (temp(:,:).lt.89.0) ! not far from Fray and Schmidt, 2009
-           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) ! not far from Fray and Schmidt, 2009
-           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
-
-     else if(trim(cnames(ic)).eq."C2H6") then
-
-!        where (temp(:,:).lt.90.)
-!           ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
-!                / press(:,:) * 1013.25e0 / 760.0e0
-!           ysat(:,:,ic) = 1.E5 * exp ( 15.11 - 2207./temp(:,:) - 24110./(temp(:,:)**2)     &
-!                        + 7.744E5/(temp(:,:)**3) - 1.161E7/(temp(:,:)**4)                  &
-!                        + 6.763E7/(temp(:,:)**5) ) / press(:,:) ! Fray and Schmidt (2009)
-!        elsewhere
-!           ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                &
-!                * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
-!        endwhere
-         ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.511e1 - 2.207e3/temp(:,:)            &
-                      - 2.411e4/temp(:,:)**2 + 7.744e5/temp(:,:)**3 - 1.161e7/temp(:,:)**4  &
-                      + 6.763e7/temp(:,:)**5)
-
-     else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
-
-        ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
+  ! Main loop
+  do ic = 1, nkim
+      
+      ! CH4 :
+      !------
+      if(trim(cnames(ic)).eq."CH4") then
+         !where (temp(:,:).lt.90.65)                                      
+         !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                 &
+         !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                   &
+         !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
+         !elsewhere
+         !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( (154567.02e0 / temp(:,:) - 1598.8512e0)    &
+         !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
+         !endwhere
+
+         ! Fray and Schmidt (2009)
+         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)                &
+                        - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
+
+         ! Forcing CH4 to 1.4% minimum               
+         where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
+
+      ! C2H2 :
+      !-------
+      else if(trim(cnames(ic)).eq."C2H2") then
+         !ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                       &
+         !     * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
+
+         ! Fray and Schmidt (2009)
+         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:))
+
+      ! C2H4 :
+      !-------
+      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
+      
+      ! C2H6 :
+      !-------
+      else if(trim(cnames(ic)).eq."C2H6") then
+         where (temp(:,:).lt.90.)
+            !ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
+            !      / press(:,:) * 1013.25e0 / 760.0e0
+            ! Fray and Schmidt (2009)
+            ysat(:,:,ic) = (1.e3 / press(:,:)) * exp ( 15.11 - 2207./temp(:,:)               &
+                           - 24110./(temp(:,:)**2) + 7.744E5/(temp(:,:)**3)                  &
+                           - 1.161E7/(temp(:,:)**4) + 6.763E7/(temp(:,:)**5)) 
+         elsewhere
+            ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                 &
+                  * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
+         endwhere
+      
+      ! CH3CCH :
+      !---------
+      else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
+         ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
              / (1.15e0*temp(:,:) - 37.485e0) ) / press(:,:) * 1013.25e0 / 760.0e0
 
-     else if(trim(cnames(ic)).eq."C3H6")  then
-
-        ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
+      ! C3H6 :
+      !-------
+      else if(trim(cnames(ic)).eq."C3H6")  then
+         ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+
+      ! C3H8 :
+      !-------
      else if(trim(cnames(ic)).eq."C3H8")  then
-
         ysat(:,:,ic) = 10.0**(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
 
-     else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
-
-        ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
+      ! C4H2 :
+      !-------
+      else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
+         ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
              * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."C4H4")  then
-
-        ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
-
-     else if(trim(cnames(ic)).eq."C4H6")  then
-
-        ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
+      
+      ! C4H4 :
+      !-------
+      else if(trim(cnames(ic)).eq."C4H4")  then
+         ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
+
+      ! C4H6 :
+      !-------
+      else if(trim(cnames(ic)).eq."C4H6")  then
+         ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
              /(1.15e0*temp(:,:) - 39.345e0) ) / press(:,:) * 1013.25e0 / 760.0e0
 
-     else if(trim(cnames(ic)).eq."C4H10")  then
-
-        ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."C6H2")  then
-
-        ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
+      ! C4H10 :
+      !--------
+      else if(trim(cnames(ic)).eq."C4H10")  then
+         ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+
+      ! C6H2 :
+      !-------
+      else if(trim(cnames(ic)).eq."C6H2")  then
+         ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
              alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
 
-     else if(trim(cnames(ic)).eq."C8H2")  then
-
-        ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
+      ! C8H2 :
+      !-------
+      else if(trim(cnames(ic)).eq."C8H2")  then
+         ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
              alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
 
-     else if(trim(cnames(ic)).eq."AC6H6")  then
-
-        !x = 1.0e0 - temp(:,:) / 562.2e0
-        !ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
-        !     - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
-
-        ysat(:,:,ic) = 1.E5 * exp (17.35-5663./temp(:,:)) / press(:,:) ! Fray and Schmidt (2009)
-        
-     else if(trim(cnames(ic)).eq."HCN")  then
-
-        !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-        ysat(:,:,ic) = (1000. / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)             &
-                     - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4) ! Fray and Schmidt (2009)
-
-     else if(trim(cnames(ic)).eq."CH3CN")  then
-
-        ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."C2H3CN")  then
-
-        ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."NCCN")  then
-
-        ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."HC3N")  then
-
-        ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     else if(trim(cnames(ic)).eq."C4N2")  then
-
-        ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
-
-     endif
+      ! C6H6 :
+      !-------
+      else if(trim(cnames(ic)).eq."AC6H6")  then
+         !x = 1.0e0 - temp(:,:) / 562.2e0
+         !ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
+         !     - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
+
+         ! Fray and Schmidt (2009)
+         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (17.35-5663./temp(:,:))
+
+      ! HCN :
+      !-------
+      else if(trim(cnames(ic)).eq."HCN")  then
+         !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
+         
+         ! Fray and Schmidt (2009)
+         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)               &
+                        - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4) 
+
+      ! CH3CN :
+      !--------
+      else if(trim(cnames(ic)).eq."CH3CN")  then
+         ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+   
+      ! C2H3CN :
+      !-------
+      else if(trim(cnames(ic)).eq."C2H3CN")  then
+         ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
+
+      ! NCCN :
+      !-------
+      else if(trim(cnames(ic)).eq."NCCN")  then
+         ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+   
+      ! HC3N :
+      !-------
+      else if(trim(cnames(ic)).eq."HC3N")  then
+         ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+   
+      ! C4N2 :
+      !-------
+      else if(trim(cnames(ic)).eq."C4N2")  then
+         ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
+
+      endif
   enddo
 
Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 3090)
@@ -35,6 +35,13 @@
 !     Emmanuel 01/2001, Forget 09/2001
 !     Robin Wordsworth (2009)
-!     Jan Vatant d'Ollone (2018) -> corrk recombining case
-!     Bruno de Batz de Trenquelléon (2023) --> Optics of Titan's haze and coulds
+!     
+!     Modified
+!     --------
+!     Jan Vatant d'Ollone (2018)
+!        --> corrk recombining case
+!     B. de Batz de Trenquelléon (2023)
+!        --> Titan's haze and could optics
+!        --> clear/dark columns method
+!        --> optical diagnostics
 !
 !==================================================================
@@ -80,8 +87,9 @@
       REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags ().
       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,5)  ! IR 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).
       
       
@@ -97,50 +105,53 @@
       ! Optical values for the optci/cv subroutines
       REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
-      
+      ! IR
       REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
+      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      ! VI
+      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
+      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      
+      ! Temporary optical values for the optci/cv subroutines (clear column)
       REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
+      REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
+      REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      ! Temporary optical values for the optci/cv subroutines (dark column)
       REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
       REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      
-      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
+      REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
+      REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
+      REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
+      REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
       REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      
-      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-      
-      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
-      REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
-      
-      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
-      REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
-      REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
-      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
-      REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
-      REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
-      
+  
+      ! 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)
+      ! 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)
+      ! 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_hazei(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_hazei_cc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_hazei_dc(L_LEVELS,L_NSPECTI,3)
-
-      REAL*8 popt_cloudsv(L_LEVELS,L_NSPECTV,3)
-      REAL*8 popt_cloudsv_cc(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 popt_cloudsi(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_cloudsi_cc(L_LEVELS,L_NSPECTI,3)
-      REAL*8 popt_cloudsi_dc(L_LEVELS,L_NSPECTI,3)
+
 
       REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
@@ -546,19 +557,21 @@
          endwhere
          ! Diagnostics for clouds :
-         where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
-            popt_cloudsv_cc(:,:,1) = 1.d-30
-         endwhere
-         where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
-            popt_cloudsv_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
+         if (callclouds) then
+            where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
+               popt_cloudsv_cc(:,:,1) = 1.d-30
+            endwhere
+            where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
+               popt_cloudsv_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
+         endif
          
          if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
@@ -666,19 +679,21 @@
          endwhere
          ! Diagnostics for clouds :
-         where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
-            popt_cloudsi_cc(:,:,1) = 1.d-30
-         endwhere
-         where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
-            popt_cloudsi_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
+         if (callclouds) then
+            where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
+               popt_cloudsi_cc(:,:,1) = 1.d-30
+            endwhere
+            where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
+               popt_cloudsi_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
+         endif
 
          call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
@@ -758,8 +773,7 @@
       
       
-      ! Optical properties of haze and clouds (dtau, tau, k, w, g) :
+      ! Diagnostics : optical properties of haze and clouds (dtau, tau, k, w, g) :
       do l = 2, L_LEVELS, 2
          lev2lay = L_NLAYRAD+1 - l/2
-         
          ! Visible :
          do nw = 1, L_NSPECTV
@@ -768,20 +782,27 @@
             ! Optical thickness (dtau) :
             popthv(ig,lev2lay,nw,1) = (popt_hazev(l,nw,1) + popt_hazev(l+1,nw,1)) / 2.0
-            poptcv(ig,lev2lay,nw,1) = (popt_cloudsv(l,nw,1) + popt_cloudsv(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
+            endif
             ! Opacity (tau) :
             do k = L_NLAYRAD, lev2lay, -1
                popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
-               poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
+               if (callclouds) then
+                  poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
+               endif
             enddo
             ! Extinction (k) :
             popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
-            poptcv(ig,lev2lay,nw,3) = poptcv(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))
+            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
-            poptcv(ig,lev2lay,nw,4) = (popt_cloudsv(l,nw,2) + popt_cloudsv(l+1,nw,2)) / 2.0
             popthv(ig,lev2lay,nw,5) = (popt_hazev(l,nw,3) + popt_hazev(l+1,nw,3)) / 2.0
-            poptcv(ig,lev2lay,nw,5) = (popt_cloudsv(l,nw,3) + popt_cloudsv(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
+            endif
          enddo
-
          ! Infra-Red
          do nw=1,L_NSPECTI
@@ -790,18 +811,26 @@
             ! Optical thickness (dtau) :
             popthi(ig,lev2lay,nw,1) = (popt_hazei(l,nw,1) + popt_hazei(l+1,nw,1)) / 2.0
-            poptci(ig,lev2lay,nw,1) = (popt_cloudsi(l,nw,1) + popt_cloudsi(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
+            endif
             ! Opacity (tau) :
             do k = L_NLAYRAD, lev2lay, -1
                popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
-               poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
+               if (callclouds) then
+                  poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
+               endif
             enddo
             ! Extinction (k) :
             popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
-            poptci(ig,lev2lay,nw,3) = poptci(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))
+            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
-            poptci(ig,lev2lay,nw,4) = (popt_cloudsi(l,nw,2) + popt_cloudsi(l+1,nw,2)) / 2.0
             popthi(ig,lev2lay,nw,5) = (popt_hazei(l,nw,3) + popt_hazei(l+1,nw,3)) / 2.0
-            poptci(ig,lev2lay,nw,5) = (popt_cloudsi(l,nw,3) + popt_cloudsi(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
+            endif
          enddo
       enddo
@@ -826,5 +855,4 @@
                       real(-nfluxtopv),real(nfluxtopi) 
          close(116)
-
 
 !          USEFUL COMMENT - Do Not Remove.
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 3090)
@@ -58,4 +58,5 @@
       real,save    :: FHIR
 !$OMP THREADPRIVATE(Fcloudy,Fssa,FHVIS,FHIR)
+      
       logical,save :: resCH4_inf
 !$OMP THREADPRIVATE(resCH4_inf)
Index: trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 3090)
@@ -18,5 +18,6 @@
   !! done elsewhere.
   !! 
-  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017, B. de Batz de Trenquelléon (2023)
+  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017
+  !! Modified : B. de Batz de Trenquelléon (2023)
   !!
   USE MMP_GCM
Index: trunk/LMDZ.TITAN/libf/phytitan/cond_muphy.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/cond_muphy.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/cond_muphy.F90	(revision 3090)
@@ -24,6 +24,6 @@
 !   dtlc = Condensation heating rate [K.s-1]
 !
-!   Author(s)
-!   ---------
+!   Author
+!   ------
 !   B. de Batz de Trenquelléon (07/2023)
 !====================================================================    
Index: trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 3090)
@@ -17,4 +17,6 @@
 !     Jan Vatant d'Ollone (2016)
 !
+!  Modified :
+!     B. de Batz de Trenquelléon (2022)
 ! ==========================================================================
 
Index: trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/evapCH4.F90	(revision 3090)
@@ -44,6 +44,6 @@
 !
 !
-!   Author(s)
-!   ---------
+!   Author
+!   ------
 !   B. de Batz de Trenquelléon (10/2022)
 !
@@ -79,5 +79,5 @@
 real, parameter :: humCH4 = 0.5     ! Imposed surface humidity for CH4 [-]
 
-real, parameter :: Flnp = 0.01      ! Fraction occupied by lakes (North Pole)
+real, parameter :: Flnp = 0.05      ! Fraction occupied by lakes (North Pole)
 real, parameter :: Flsp = 0.005     ! Fraction occupied by lakes (South Pole)
 
@@ -129,6 +129,6 @@
   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
-  ! CH4 : 0.85 * qsat because of dissolution in N2
-  qsatCH4 = 0.85 * qsatCH4
+  ! CH4 : 0.80 * qsat because of dissolution in N2
+  qsatCH4 = 0.80 * qsatCH4
 
   ! Flux at the surface [kg.m-2.s-1] :
@@ -196,32 +196,4 @@
   dtsurfevap(ig) = - (fluxCH4 * LvCH4) / Cplake
 
-  ! >>> [TEMPO : BBT]
-  !open(501,file='Evap_CH4.out')
-  !if(REAL(latitude_deg(ig)) .ge. 70.) then
-  !      write(501,*) '-----------------------------------'
-  !      write(501,*) 'Latitude =', REAL(latitude_deg(ig))
-  !      write(501,*) ''
-  !      write(501,*) 'Initialisation variables :'
-  !     write(501,*) 'Air density =', rhoair, 'kg.m-3'
-  !      write(501,*) 'Horizontal wind =', ws, 'm.s-1'
-  !      write(501,*) 'Turbulent term =', Cd
-  !      write(501,*) 'CH4 saturation at the surface =', qsatCH4, 'mol.mol-1'
-  !      write(501,*) ''
-  !      write(501,*) 'Flux :'
-  !      write(501,*) 'Flux of CH4 =', fluxCH4, 'kg.m-2.s-1.mol.mol-1'
-  !      write(501,*) ''
-  !      write(501,*) 'Variation of CH4 :'
-  !      write(501,*) 'New molar fraction of CH4 =', newpqCH4, 'mol.mol-1'
-  !      write(501,*) 'Trend of CH4 =', pdqCH4(ig), 'mol.mol-1.s-1'
-  !      write(501,*) 'Variation of CH4 tank =', dtankCH4, 'm'
-  !      write(501,*) ''
-  !      write(501,*) 'Latent Heat :'
-  !      write(501,*) 'Latente heat of CH4 vaporisation =', LvCH4, 'J.kg-1.mol.mol-1.s-1'
-  !      write(501,*) 'dtsurfevap :', dtsurfevap(ig), 'K.s-1'
-  !      write(501,*) '-----------------------------------'
-  !endif
-  !close(501)
-  ! <<< [TEMPO : BBT]
-
 enddo ! End of main loop on the horizontal grid
 
Index: trunk/LMDZ.TITAN/libf/phytitan/get_haze_and_cloud_opacity.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/get_haze_and_cloud_opacity.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/get_haze_and_cloud_opacity.F90	(revision 3090)
@@ -43,6 +43,6 @@
   !
   !
-  !     Authors
-  !     -------
+  !     Author
+  !     ------
   !     B. de Batz de Trenquelleon (08/2022).
   !
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 3090)
@@ -33,9 +33,12 @@
 !   author: Frederic Hourdin 15 / 10 /93
 !   -------
-!   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
-!             Ehouarn Millour (oct. 2008) tracers are now identified
-!              by their names and may not be contiguously
-!              stored in the q(:,:,:,:) array
-!             E.M. (june 2009) use getin routine to load parameters
+!   modified: -Sebastien Lebonnois 11/06/2003 (new callphys.def)
+!             -Ehouarn Millour (oct. 2008) tracers are now identified
+!               by their names and may not be contiguously
+!               stored in the q(:,:,:,:) array
+!             -Ehouarn Millour (june 2009) use getin routine
+!               to load parameters
+!             -Bruno de Batz de Trenquelléon (2023) minor changes and
+!               new initialisations
 !
 !
Index: trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 3090)
@@ -23,5 +23,5 @@
   !     -------
   !     J. Vatant d'Ollone, J.Burgalat, S.Lebonnois (09/2017)
-  !     B. de Batz de Trenquelléon (02/2023)
+  !     Modified : B. de Batz de Trenquelléon (02/2023)
   !
   !============================================================================
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 3090)
@@ -33,8 +33,12 @@
   !     -------
   !     Adapted from the NASA Ames code by R. Wordsworth (2009)
-  !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
-  !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
-  !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
-  !     The whole routine needs to be redone...
+  !     
+  !     Modified
+  !     --------
+  !     J. Vatant d'Ollone (2016-17)
+  !         --> Clean and adaptation to Titan
+  !     B. de Batz de Trenquelléon (2022-2023)
+  !         --> Clean and correction to Titan
+  !         --> New optics added for Titan's clouds
   !     
   !==================================================================
@@ -105,5 +109,5 @@
   ! Variables for new optics
   integer iq, iw, FTYPE, CTYPE
-  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
+  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
   real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
   real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
@@ -178,8 +182,8 @@
                m3as = pqmo(ilay,2) / 2.0
                m3af = pqmo(ilay,4) / 2.0
-               ! Cut-off
-               IF (ilay .lt. 19) THEN
-                  m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                  m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+               ! Cut-off (here for p = 2.7e3Pa / alt = 70km)
+               IF (ilay .lt. 23) THEN
+                  m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                  m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                ENDIF
 
@@ -202,9 +206,9 @@
                m0as  = pqmo(ilay,1) / 2.0
                m3as  = pqmo(ilay,2) / 2.0
-               ! If not callclouds : must have a cut-off
+               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
                IF (.NOT. callclouds) THEN
-                  IF (ilay .lt. 19) THEN
-                     m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                     m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+                  IF (ilay .lt. 23) THEN
+                     m0as = pqmo(23,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                     m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                   ENDIF
                ENDIF
@@ -216,9 +220,9 @@
                m0af  = pqmo(ilay,3) / 2.0
                m3af  = pqmo(ilay,4) / 2.0
-               ! If not callclouds : must have a cut-off
+               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
                IF (.NOT. callclouds) THEN
-                  IF (ilay .lt. 19) THEN
-                     m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                     m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+                  IF (ilay .lt. 23) THEN
+                     m0af = pqmo(23,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                     m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                   ENDIF
                ENDIF
@@ -249,6 +253,5 @@
                m0ccn = pqmo(ilay,5) / 2.0
                m3ccn = pqmo(ilay,6) / 2.0
-               m3ices = 0.0d0
-               m3cld  = m3ccn
+               m3cld  = 0.0d0
                
                ! Clear / Dark column method :
@@ -261,5 +264,4 @@
                IF (CDCOLUMN == 0) THEN
                   DO iq = 2, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -269,5 +271,4 @@
                ELSEIF (CDCOLUMN == 1) THEN
                   DO iq = 1, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -278,5 +279,4 @@
                   WRITE(*,*) 'WE USE DARK COLUMN ...'
                   DO iq = 1, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -285,5 +285,4 @@
                
                ! For small dropplets, opacity of nucleus dominates...
-               !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
                ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 3090)
@@ -25,8 +25,12 @@
   !     -------
   !     Adapted from the NASA Ames code by R. Wordsworth (2009)
-  !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
-  !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
-  !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
-  !     The whole routine needs to be redone...
+  !     
+  !     Modified
+  !     --------
+  !     J. Vatant d'Ollone (2016-17)
+  !         --> Clean and adaptation to Titan
+  !     B. de Batz de Trenquelléon (2022-2023)
+  !         --> Clean and correction to Titan
+  !         --> New optics added for Titan's clouds
   !     
   !==================================================================
@@ -116,5 +120,5 @@
   ! Variables for new optics
   integer iq, iw, FTYPE, CTYPE
-  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
+  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
   real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
   real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
@@ -202,8 +206,8 @@
                m3as = pqmo(ilay,2) / 2.0
                m3af = pqmo(ilay,4) / 2.0
-               ! Cut-off
-               IF (ilay .lt. 19) THEN
-                  m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                  m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+               ! Cut-off (here for p = 2.7e3Pa / alt = 70km)
+               IF (ilay .lt. 23) THEN
+                  m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                  m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                ENDIF
                
@@ -226,9 +230,9 @@
                m0as  = pqmo(ilay,1) / 2.0
                m3as  = pqmo(ilay,2) / 2.0
-               ! If not callclouds : must have a cut-off
+               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
                IF (.NOT. callclouds) THEN
-                  IF (ilay .lt. 19) THEN
-                     m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                     m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+                  IF (ilay .lt. 23) THEN
+                     m0as = pqmo(23,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                     m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                   ENDIF
                ENDIF
@@ -240,9 +244,9 @@
                m0af  = pqmo(ilay,3) / 2.0
                m3af  = pqmo(ilay,4) / 2.0
-               ! If not callclouds : must have a cut-off
+               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
                IF (.NOT. callclouds) THEN
-                  IF (ilay .lt. 19) THEN
-                     m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
-                     m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
+                  IF (ilay .lt. 23) THEN
+                     m0af = pqmo(23,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
+                     m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
                   ENDIF
                ENDIF
@@ -273,6 +277,5 @@
                m0ccn = pqmo(ilay,5) / 2.0
                m3ccn = pqmo(ilay,6) / 2.0
-               m3ices = 0.0d0
-               m3cld  = m3ccn
+               m3cld = 0.0d0
                
                ! Clear / Dark column method :
@@ -285,5 +288,4 @@
                IF (CDCOLUMN == 0) THEN
                   DO iq = 2, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -293,5 +295,4 @@
                ELSEIF (CDCOLUMN == 1) THEN
                   DO iq = 1, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -302,5 +303,4 @@
                   WRITE(*,*) 'WE USE DARK COLUMN ...'
                   DO iq = 1, nice
-                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
                      m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
                   ENDDO
@@ -309,5 +309,4 @@
 
                ! For small dropplets, opacity of nucleus dominates
-               !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
                ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
                ssa_cld(nw) = ssa_cld(nw) * Fssa
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 3090)
@@ -15,5 +15,5 @@
  
       use radinc_h, only : L_NSPECTI,L_NSPECTV
-      use radcommon_h, only: sigma, gzlat, grav, BWNV
+      use radcommon_h, only: sigma, gzlat, grav, BWNV, WNOI, WNOV
       use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
@@ -161,7 +161,8 @@
 !                + clean of all too-generic (ocean, water, co2 ...) routines
 !                + Titan's chemistry
-!           Microphysical moment model - J.Burgalat / J.Vatant d'Ollone (2017-2018)
-!           Improvement of microphysical moment model - B. de Batz de Trenquelléon (2022-2023)
-!           Optics ofnn Titan's coulds - B. de Batz de Trenquelléon (2022-2023)
+!           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 !!
 !============================================================================================
 
@@ -258,5 +259,5 @@
       real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
       
-      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
+      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency (chemistry routine).
       
       real zdqmufi(ngrid,nlayer,nq)   ! Microphysical tendency.
@@ -271,5 +272,5 @@
       real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
       real zdhadj(ngrid,nlayer)                      ! Convadj routine.
-      real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind (m/s/s).
+      real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind.
      
       ! For Pressure and Mass :
@@ -301,5 +302,5 @@
       real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
 
-      real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u*u+v*v))
+      real zhorizwind(ngrid,nlayer) ! Horizontal Wind (sqrt(u*u+v*v))
 
       real vmr(ngrid,nlayer)                        ! volume mixing ratio
@@ -352,5 +353,5 @@
       real, dimension(ngrid,nlayer,nkim) :: ysat
       
-      ! Fraction of CH4 [mol/mol]
+      ! Temporary fraction of CH4 [mol/mol]
       real, dimension(ngrid,nlayer) :: tpq_CH4
 
@@ -467,7 +468,7 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-         if ( callchim ) then
-
-            if ( moyzon_ch .and. ngrid.eq.1 ) then
+         if (callchim) then
+
+            if (moyzon_ch .and. ngrid.eq.1) then
                print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
                print *, "Please desactivate zonal mean for 1D !"
@@ -483,5 +484,5 @@
 !        ~~~~~~~~~~~~~~~~~~~~~~~~
 
-         IF ( callmufi ) THEN
+         IF (callmufi) THEN
            ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
            call inimufi(ptimestep)
@@ -619,5 +620,5 @@
          write(*,*) "physiq: call initialize_xios_output"
          call initialize_xios_output(pday,ptime,ptimestep,daysec, &
-                                     presnivs,pseudoalt)
+                                     presnivs,pseudoalt,wnoi,wnov)
 #endif
          write(*,*) "physiq: end of firstcall"
@@ -1066,14 +1067,10 @@
       if (tracer) then
 
-
 #ifndef MESOSCALE
-!! JVO 20 : For now, no chemistry or microphysics in MESOSCALE, but why not in the future ?
-
-  ! -------------------
-  !   V.1 Microphysics
-  ! -------------------
-
-         ! JVO 05/18 : We must call microphysics before chemistry, for condensation !
- 
+   ! -------------------
+   !   V.1 Microphysics
+   ! -------------------
+
+         ! We must call microphysics before chemistry, for condensation !
          if (callmufi) then
 
@@ -1084,11 +1081,10 @@
             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
             tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
-#else
-            
+
+#else   
             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
-
             pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
 
-            ! Sanity check ( way safer to be done here rather than within YAMMS )
+            ! Sanity check (way safer to be done here rather than within YAMMS)
             ! Important : the sanity check intentionally include the former processes tendency !
             ! NB : Despite this sanity check there might be still some unphysical values going through :
@@ -1124,21 +1120,22 @@
 
             ! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
-            if ( callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0 ) then
-              zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
-              call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
-                           gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
-              ! TODO : Add a sanity check here !
+            if (callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0) then
+               zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
+               call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
+                            gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
+               ! TODO : Add a sanity check here !
             endif
          
-            ! >>> TEMPO : BBT
-            ! Default value -> no condensation [kg/kg_air/s] :
-            dmuficond(:,:,:) = 0.D0
-            do iq = 1, size(ices_indx)
-               dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
-            enddo
-            call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
-            !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
-            ! <<< TEMPO : BBT
-         endif
+            ! Condensation heating rate :
+            if (callclouds) then
+               ! Default value -> no condensation [kg/kg_air/s] :
+               dmuficond(:,:,:) = 0.D0
+               do iq = 1, size(ices_indx)
+                  dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
+               enddo
+               call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
+               !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
+            endif
+         endif ! callmufi
       
   ! -----------------
@@ -1146,8 +1143,4 @@
   ! -----------------
   !   NB : Must be call last ( brings fields back to an equilibrium )
-
-  ! Known bug ? ( JVO 18 ) : If you try to use chimi_indx instead of iq+nmicro
-  ! it leads to weird results / crash on dev mod ( ok in debug ) ... Why ? Idk ...
-
          if (callchim) then
 
@@ -1155,5 +1148,5 @@
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             do iq = 1,nkim
-               ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep ) / rat_mmol(iq+nmicro)
+               ychim(:,:,iq) = (pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep) / rat_mmol(iq+nmicro)
             enddo
 
@@ -1162,5 +1155,5 @@
             ! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal 
             ! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
-            if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
+            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
               do iq = 1,nkim
                 ychimbar(:,:,iq) =  zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
@@ -1174,11 +1167,12 @@
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             
-            call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point
+            call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point (!!p in mbar!!)
 
             dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
 
             do iq=1,nkim
-               where ( ychim(:,:,iq).gt.ysat(:,:,iq) )   &
-                     dyccond(:,:,iq+nmicro) = ( -ychim(:,:,iq)+ysat(:,:,iq) ) / ptimestep
+               where (ychim(:,:,iq).gt.ysat(:,:,iq))
+                  dyccond(:,:,iq+nmicro) = (-ychim(:,:,iq)+ysat(:,:,iq)) / ptimestep
+               endwhere
             enddo
 
@@ -1191,5 +1185,4 @@
             do iq=1,nkim
               ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
-
               pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
             enddo
@@ -1198,14 +1191,15 @@
             ! ii. 2D zonally averaged fields need to condense and evap before photochem
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
-
-              call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
-
-              dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
-              
-              do iq = 1,nkim
-                 where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )  &
-                       dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
-              enddo
+            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
+
+               call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
+
+               dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
+               
+               do iq = 1,nkim
+                     where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )
+                        dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
+                     endwhere
+               enddo
 
                if (callclouds) then
@@ -1215,57 +1209,50 @@
                endif
 
-              do iq=1,nkim
-                ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
-              enddo
-
-              ! Pseudo-evap ( forcing constant surface humidity )
-              !do ig=1,ngrid
-              !   if ( ychimbar(ig,1,7) .lt. botCH4 ) ychimbar(ig,1,7) = botCH4
-              !enddo
-
-            endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 )
+               do iq=1,nkim
+                  ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
+               enddo
+
+            endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq.0 )
 
             ! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
             if( mod(icount-1,ichim).eq.0. ) then
-               
-               print *, "We enter in the chemistry ..."
+               print *, "We enter in the photochemistry ..."
 
                if (moyzon_ch) then ! 2D zonally averaged chemistry
-
                  ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
                  call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar,  &
                               zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
+               
                else ! 3D chemistry (or 1D run)
                  call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi,  &
                               pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
                endif ! if moyzon
-
-            endif
+            endif ! if (mod(icount-1,ichim).eq.0)
             
             ! iv. Surface pseudo-evaporation
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            ! Infinite tank of CH4 (global)
             if (resCH4_inf) then
-                do ig=1,ngrid
-                        if ( (ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4 ) then ! + dycchi because ychim not yet updated
-                                dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
-                        else
-                                dycevapCH4(ig) = 0.D0
-                        endif
-                enddo
-
-            ! >>> New evaporation flux
+               do ig=1,ngrid
+                  if ((ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4) then ! + dycchi because ychim not yet updated
+                     dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
+                  else
+                     dycevapCH4(ig) = 0.D0
+                  endif
+               enddo
+
+            ! Tank at the pole (for the moment infinit...)
             else
-                tankCH4(:) = 2. ! Pour le moment reservoir infini...
-                if (moyzon_ch) then
-                    tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
-                else
-                    tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
-                endif
-                call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
-                             pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
-                zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
+               tankCH4(:) = 2.
+               if (moyzon_ch) then
+                 tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
+               else
+                 tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
+               endif
+               call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
+                            pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
+               zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
             endif
-            ! <<< New evaporation flux
 
             pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro)
@@ -1286,4 +1273,5 @@
          endif ! end of 'callchim'
 
+! END MESOSCALE
 #endif
 
@@ -1316,4 +1304,5 @@
             pdv(:,:)    = pdv(:,:)    + zdvmr(:,:)
             pdpsrf(:)   = pdpsrf(:)   + zdpsrfmr(:)
+            zdtsurf(:)  = zdtsurf(:)  + zdtsurfmr(:)
          endif
 
@@ -1362,18 +1351,14 @@
       if (nudging_u) then
          zdundg(:,:) = 0.D0
-         
          ! boucle sur les points de grille :
          do i = 1, ngrid
             ! boucle sur les latitudes :
             do j = 1, 49
-
                ! Nudging of the first 23 layers only !!
                if (ABS(REAL(u_ref(j,1)) - REAL(latitude_deg(i))) .lt. 1e-2) then
                   zdundg(i,1:23) = (u_ref(j,2:24) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
                endif
-            
             enddo
          enddo
-
          pdu(:,:) = pdu(:,:) + zdundg(:,:)
       endif
@@ -1406,10 +1391,10 @@
       enddo
 
-      ! Tracers.
-      ! >>> TEMPO : BBT
-      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
-      ! <<< TEMPO : BBT
+      ! 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
       zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
 
@@ -1501,5 +1486,4 @@
       endif
 
-
       if (is_master) print*,'--> Ls =',zls*180./pi
       
@@ -1555,6 +1539,9 @@
       call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
       call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
-
-      ! Total energy balance diagnostics
+      
+      ! Subsurface temperatures
+      !call writediagsoil(ngrid,"tempsoil","temperature soil","K",3,tsoil)
+
+     ! Total energy balance diagnostics
       if(callrad.and.(.not.newtonian))then
          call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
@@ -1666,5 +1653,5 @@
       ! Send fields to XIOS: (NB these fields must also be defined as
       ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
-      
+
       !--------------------------------------------------------
       ! General diagnostics :
@@ -1697,4 +1684,8 @@
       ENDIF
 
+      ! Subsurface temperatures (3D)
+      !CALL send_xios_field("tempsoil",tsoil)
+      ! Total energy balance diagnostics (2D)
+      !CALL send_xios_field("ALB",albedo_equivalent)
 
       !--------------------------------------------------------
@@ -1709,5 +1700,7 @@
       CALL send_xios_field("dudif",zdudif)
       CALL send_xios_field("duadj",zduadj)
-      CALL send_xios_field("dundg",zdundg)
+      IF (nudging_u) THEN
+         CALL send_xios_field("dundg",zdundg)
+      ENDIF
       
       ! zhorizwind = sqrt(u*u + v*v)
@@ -1727,5 +1720,7 @@
       CALL send_xios_field("dtdif",zdtdif)
       CALL send_xios_field("dtadj",zdtadj(:,:))
-      CALL send_xios_field("dtlc",zdtlc)
+      IF (callclouds) THEN
+         CALL send_xios_field("dtlc",zdtlc)
+      ENDIF
 
       ! dtrad = zdtsw + zdtlw
@@ -1752,53 +1747,33 @@
       ! Optical diagnostics :
       !--------------------------------------------------------
-      ! Optics diagnostics for haze (Channel VI04 : 2.80um).
-      CALL send_xios_field('tauhv_04',zpopthv(:,:,4,2))
-      CALL send_xios_field('khv_04',zpopthv(:,:,4,3))
-      CALL send_xios_field('wbarhv_04',zpopthv(:,:,4,4))
-      CALL send_xios_field('gbarhv_04',zpopthv(:,:,4,5))
-      ! Optics diagnostics for haze (Channel VI18 : 0.85um).
-      CALL send_xios_field('tauhv_18',zpopthv(:,:,18,2))
-      CALL send_xios_field('khv_18',zpopthv(:,:,18,3))
-      CALL send_xios_field('wbarhv_18',zpopthv(:,:,18,4))
-      CALL send_xios_field('gbarhv_18',zpopthv(:,:,18,5))
-      ! Optics diagnostics for haze (Channel IR07 : 29.5um).
-      CALL send_xios_field('tauhi_07',zpopthi(:,:,7,2))
-      CALL send_xios_field('khi_07',zpopthi(:,:,7,3))
-      CALL send_xios_field('wbarhi_07',zpopthi(:,:,7,4))
-      CALL send_xios_field('gbarhi_07',zpopthi(:,:,7,5))
-      ! Optics diagnostics for haze (Channel IR22 : 5.66um).
-      CALL send_xios_field('tauhi_22',zpopthi(:,:,22,2))
-      CALL send_xios_field('khi_22',zpopthi(:,:,22,3))
-      CALL send_xios_field('wbarhi_22',zpopthi(:,:,22,4))
-      CALL send_xios_field('gbarhi_22',zpopthi(:,:,22,5))
-
-      ! Optics diagnostics for haze + clouds (Channel VI04 : 2.80um).
-      CALL send_xios_field('taucv_04',zpoptcv(:,:,4,2))
-      CALL send_xios_field('kcv_04',zpoptcv(:,:,4,3))
-      CALL send_xios_field('wbarcv_04',zpoptcv(:,:,4,4))
-      CALL send_xios_field('gbarcv_04',zpoptcv(:,:,4,5))
-      ! Optics diagnostics for haze + clouds (Channel VI18 : 0.85um).
-      CALL send_xios_field('taucv_18',zpoptcv(:,:,18,2))
-      CALL send_xios_field('kcv_18',zpoptcv(:,:,18,3))
-      CALL send_xios_field('wbarcv_18',zpoptcv(:,:,18,4))
-      CALL send_xios_field('gbarcv_18',zpoptcv(:,:,18,5))
-      ! Optics diagnostics for haze + clouds (Channel IR07 : 29.5um).
-      CALL send_xios_field('tauci_07',zpoptci(:,:,7,2))
-      CALL send_xios_field('kci_07',zpoptci(:,:,7,3))
-      CALL send_xios_field('wbarci_07',zpoptci(:,:,7,4))
-      CALL send_xios_field('gbarci_07',zpoptci(:,:,7,5))
-      ! Optics diagnostics for haze + clouds (Channel IR22 : 5.66um).
-      CALL send_xios_field('tauci_22',zpoptci(:,:,22,2))
-      CALL send_xios_field('kci_22',zpoptci(:,:,22,3))
-      CALL send_xios_field('wbarci_22',zpoptci(:,:,22,4))
-      CALL send_xios_field('gbarci_22',zpoptci(:,:,22,5))
-
-      ! >>> [TEMPO : BBT]
-      !DO nw = 1, L_NSPECTV
-      !   WRITE(str2,'(i2.2)') nw
-      !   CALL send_xios_field('khv'//TRIM(nw),zpopthv(:,:,nw,3))
-      !ENDDO
-      ! <<< [TEMPO : BBT]
-      
+      ! Diagnostics for haze :
+      ! IR
+      CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
+      CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
+      CALL send_xios_field('khi',zpopthi(:,:,:,3))
+      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))
+      CALL send_xios_field('khv',zpopthv(:,:,:,3))
+      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      
 
       !--------------------------------------------------------
@@ -1831,9 +1806,6 @@
             CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
             CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
-            CALL send_xios_field("flux_iCH4",mmd_ice_fluxes(:,:,1))
-            CALL send_xios_field("flux_iC2H2",mmd_ice_fluxes(:,:,2))
-            CALL send_xios_field("flux_iC2H6",mmd_ice_fluxes(:,:,3))
-            CALL send_xios_field("flux_iHCN",mmd_ice_fluxes(:,:,3))
             DO iq = 1, size(ices_indx)
+               CALL send_xios_field('flux_i'//TRIM(nameOfTracer(gazs_indx(iq))),mmd_ice_fluxes(:,:,iq))
                CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
             ENDDO
@@ -1843,8 +1815,9 @@
          CALL send_xios_field("aer_prec",mmd_aer_prec(:))
          CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
-         CALL send_xios_field("iCH4_prec",mmd_ice_prec(:,1))
-         CALL send_xios_field("iC2H2_prec",mmd_ice_prec(:,2))
-         CALL send_xios_field("iC2H6_prec",mmd_ice_prec(:,3))
-         CALL send_xios_field("iHCN_prec",mmd_ice_prec(:,3)) 
+         IF (callclouds) THEN
+            DO iq = 1, size(ices_indx)
+               CALL send_xios_field('i'//TRIM(nameOfTracer(gazs_indx(iq)))//'_prec',mmd_ice_prec(:,iq))
+            ENDDO
+         ENDIF
       ENDIF ! of 'if callmufi'
 
@@ -1856,11 +1829,14 @@
          ! Chemical species :
          DO iq = 1, nkim
+            ! If no cloud : gzs_indx uninitialized
             CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
          ENDDO
 
          ! Condensation tendencies from microphysics (mol/mol/s) :
-         DO iq = 1, size(ices_indx)
-            CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
-         ENDDO
+         IF (callclouds) THEN
+            DO iq = 1, size(ices_indx)
+               CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
+            ENDDO
+         ENDIF
          
          ! Condensation tendencies (mol/mol/s) :
Index: trunk/LMDZ.TITAN/libf/phytitan/tpindex.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tpindex.F	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/tpindex.F	(revision 3090)
@@ -48,5 +48,5 @@
 !     -------
 !     Adapted from the NASA Ames code by R. Wordsworth (2009)
-!     Last update : B. de Batz de Trenquelleon (2023)
+!     Last update : B. de Batz de Trenquelléon (2023)
 !     
 !==================================================================
Index: trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 3090)
@@ -1,4 +1,3 @@
-! AUTHOR : J. Burgalat (27/03/2018)
-! NOTES : Add gazs_indx by B. de Batz de Trenquelléon (31/05/2022)
+! AUTHORS : J. Burgalat (2018), B. de Batz de Trenquelléon (2022)
 !
 ! OpenMP directives in this module are inconsistent:
Index: trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90	(revision 3083)
+++ trunk/LMDZ.TITAN/libf/phytitan/xios_output_mod.F90	(revision 3090)
@@ -11,5 +11,5 @@
 
  INTERFACE send_xios_field
-    MODULE PROCEDURE histwrite0d_xios,histwrite2d_xios,histwrite3d_xios
+    MODULE PROCEDURE histwrite0d_xios,histwrite2d_xios,histwrite3d_xios,histwrite4d_xios
  END INTERFACE
  
@@ -18,5 +18,5 @@
 
   SUBROUTINE initialize_xios_output(day,timeofday,dtphys,daysec,&
-                                    presnivs,pseudoalt)
+                                    presnivs,pseudoalt,wnoi,wnov)
 !  USE mod_phys_lmdz_para, only: gather, bcast, &
 !                                jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
@@ -46,5 +46,6 @@
   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.
   
   INTEGER :: i
@@ -72,5 +73,11 @@
     CALL xios_set_axis_attr("preskim_tot", n_glo=size(preskim_tot), value=preskim_tot,&
                             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")
+                  
     ! Calculation of pseudo-altitudes is still to be done
     !IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for zlaykim"
@@ -291,4 +298,41 @@
   END SUBROUTINE histwrite3d_xios
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE histwrite4d_xios(field_name, field)
+    USE dimphy, only: klon, klev
+    USE mod_phys_lmdz_para, only: gather_omp, grid1Dto2D_mpi, &
+                                  jj_nb, klon_mpi
+    USE xios, only: xios_send_field
+    USE print_control_mod, ONLY: prt_level,lunout
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon
+  
+    IMPLICIT NONE
+  
+      CHARACTER(LEN=*), INTENT(IN) :: field_name
+      REAL, DIMENSION(:,:,:), INTENT(IN) :: field ! --> field(klon,:,:)
+  
+      REAL,DIMENSION(klon_mpi,SIZE(field,2),SIZE(field,3)) :: buffer_omp
+      REAL :: Field4d(nbp_lon,jj_nb,SIZE(field,2),SIZE(field,3))
+      INTEGER :: ip, n, nlev, llev
+  
+    IF (prt_level >= 10) write(lunout,*)'Begin histrwrite4d_xios ',trim(field_name)
+  
+      !Et on.... écrit
+      IF (SIZE(field,1)/=klon) CALL abort_physic('iophy::histwrite4d','Field first DIMENSION not equal to klon',1)
+      
+      nlev=SIZE(field,2)
+      llev=SIZE(field,3)
+  
+      CALL Gather_omp(field,buffer_omp)
+  !$OMP MASTER
+      CALL grid1Dto2D_mpi(buffer_omp,field4d)
+  
+      CALL xios_send_field(field_name, Field4d(:,:,1:nlev,1:llev))
+  !$OMP END MASTER   
+  
+      IF (prt_level >= 10) write(lunout,*)'End histrwrite4d_xios ',trim(field_name)
+  END SUBROUTINE histwrite4d_xios
+  
 #endif
 
