Index: trunk/LMDZ.TITAN/libf/phytitan/calchim.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1795)
@@ -164,4 +164,5 @@
        enddo
        rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
+! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire
 
 ! lecture de tcp.ver, une seule fois
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1795)
@@ -4,6 +4,6 @@
       logical,save :: callrad,corrk,calldifv,UseTurbDiff
 !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
-      logical,save :: calladj,callsoil,callchim
-!$OMP THREADPRIVATE(calladj,callsoil,callchim)
+      logical,save :: calladj,callsoil
+!$OMP THREADPRIVATE(calladj,callsoil)
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
@@ -14,4 +14,7 @@
       logical,save :: strictboundcorrk                                     
 !$OMP THREADPRIVATE(strictboundcorrk)
+
+      logical,save :: callchim, callmufi, callclouds
+!$OMP THREADPRIVATE(callchim,callmufi,callclouds)
 
       logical,save :: enertest
@@ -36,8 +39,15 @@
       
       integer,save :: ichim
+!$OMP THREADPRIVATE(ichim)
+
       integer,save :: iddist
       integer,save :: iradia
       integer,save :: startype
-!$OMP THREADPRIVATE(ichim,iddist,iradia,startype)
+!$OMP THREADPRIVATE(iddist,iradia,startype)
+      
+      real,save :: df_mufi, rm_mufi, rho_aer_mufi
+      real,save :: p_prod, tx_prod, rc_prod
+      real,save :: air_rad
+!$OMP THREADPRIVATE(df_mufi, rm_mufi, rho_aer_mufi,p_prod,tx_prod,rc_prod,air_rad)
 
       real,save :: Fat1AU
Index: trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1795)
+++ trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1795)
@@ -0,0 +1,155 @@
+SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
+  !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
+  !!
+  !! The subroutine computes the microphysics processes for a single vertical column.
+  !!
+  !! - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere.
+  !! - All output vectors are defined from GROUND to TOP of the atmosphere.
+  !! - Only tendencies are returned.
+  !!
+  !! @important
+  !! The method assumes global initialization of YAMMS model (and extras) has been already
+  !! done elsewhere.
+  !! 
+  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017
+  !!
+  USE MMP_GCM
+  USE tracer_h
+  USE comcstfi_mod, only : g
+  USE callkeys_mod, only : callclouds
+  IMPLICIT NONE
+
+  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
+  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev  !! Altitude levels (m).
+  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
+  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
+  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K).
+
+  REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
+  REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.kg^{-1}}\)).
+  
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as !! 0th order moment of the spherical mode (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as !! 3rd order moment of the spherical mode (\(m^{3}.m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0af !! 0th order moment of the fractal mode (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3af !! 3rd order moment of the fractal mode (\(m^{3}.m^{-2}\)).
+
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: m0n  !! 0th order moment of the CCN distribution (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: m3n  !! 3rd order moment of the CCN distribution (\(m^{3}.m^{-2}\)).
+  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3i  !! 3rd order moments of the ice components (\(m^{3}.m^{-2}\)).
+  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: gazs !! Condensible species gazs molar fraction (\(mol.mol^{-1}\)).
+
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0as !! Tendency of the 0th order moment of the spherical mode distribution (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3as !! Tendency of the 3rd order moment of the spherical mode distribution (\(m^{3}.m^{-2}\)). 
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0af !! Tendency of the 0th order moment of the fractal mode distribution (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3af !! Tendency of the 3rd order moment of the fractal mode distribution (\(m^{3}.m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0n  !! Tendency of the 0th order moment of the _CCN_ distribution (\(m^{-2}\)).
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3n  !! Tendency of the 3rd order moment of the _CCN_ distribution (\(m^{3}.m^{-2}\)).
+  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dm3i  !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-2}\)). 
+  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dgazs !! Tendencies of each condensible gaz species !(\(mol.mol^{-1}\)). 
+
+  REAL(kind=8), DIMENSION(:), ALLOCATABLE ::  int2ext
+  TYPE(error) :: err
+
+  INTEGER :: ilon, i,nices
+
+  INTEGER :: nlon,nlay
+
+  ! Read size of arrays
+  nlon  = size(play,DIM=1)
+  nlay  = size(play,DIM=2)
+  nices = size(ices_indx)
+  ! Conversion intensive to extensive
+  ALLOCATE( int2ext(nlay) )  
+
+  ! Loop on horizontal grid points
+
+  ! Allocate arrays
+  ALLOCATE( m0as(nlay) )
+  ALLOCATE( m3as(nlay) )
+  ALLOCATE( m0af(nlay) )
+  ALLOCATE( m3af(nlay) )
+  ALLOCATE( m0n(nlay) )
+  ALLOCATE( m3n(nlay) )
+  ALLOCATE( m3i(nlay,nices) )
+  ALLOCATE( gazs(nlay,nices) )
+
+  ALLOCATE( dm0as(nlay) )
+  ALLOCATE( dm3as(nlay) )
+  ALLOCATE( dm0af(nlay) )
+  ALLOCATE( dm3af(nlay) )
+  ALLOCATE( dm0n(nlay) )
+  ALLOCATE( dm3n(nlay) )
+  ALLOCATE( dm3i(nlay,nices) )
+  ALLOCATE( dgazs(nlay,nices) )
+
+  ! Initialization of zdq here since intent=out and no action performed on every tracers
+  zdq(:,:,:) = 0.0
+
+  DO ilon = 1, nlon
+    ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio )
+    ! We suppose a given order of tracers !
+    int2ext(:) = ( plev(ilon,2:nlay+1) - plev(ilon,1:nlay) ) / g
+
+    m0as(:) = pq(ilon,:,1) * int2ext(:)
+    m3as(:) = pq(ilon,:,2) * int2ext(:)
+    m0af(:) = pq(ilon,:,3) * int2ext(:)
+    m3af(:) = pq(ilon,:,4) * int2ext(:)
+    
+    if (callclouds) then ! if call clouds
+      dm0n(:) = pq(ilon,:,5) * int2ext(:)
+      dm3n(:) = pq(ilon,:,6) * int2ext(:)
+      do i=1,nices
+        dm3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
+        dgazs(:,i)    = pq(ilon,:,ices_indx(i)) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
+        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
+      enddo
+    endif
+
+
+    ! Initialize YAMMS atmospheric column
+    err = mm_column_init(plev(ilon,:),zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err)
+    ! Initialize YAMMS aerosols moments column
+    err = mm_aerosols_init(m0as,m3as,m0af,m3af) ; IF (err /= 0) call abort_program(err)
+    IF (callclouds) THEN ! call clouds
+      err = mm_clouds_init(m0n,m3n,m3i,gazs) ; IF (err /= 0) call abort_program(err)
+    ENDIF
+
+    ! Check on size (???)
+
+    ! initializes tendencies:
+    !dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
+    !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp
+
+    dm0as(:) = 0.0 ; dm3as(:) = 0.0 ; dm0af(:) = 0.0 ; dm3af(:) = 0.0
+    dm0n(:) = 0.0 ; dm3n(:) = 0.0 ; dm3i(:,:) = 0.0 ; dgazs(:,:) = 0.0
+    ! call microphysics
+
+    IF (callclouds) THEN ! call clouds
+      IF(.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af,dm0n,dm3n,dm3i,dgazs)) &
+        call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
+    ELSE
+      IF (.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af)) & 
+        call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
+    ENDIF
+
+    ! Convert tracers back to intensives ( except for gazs where we work with molar mass ratio )
+    ! We suppose a given order of tracers !
+
+    zdq(ilon,:,1) = dm0as(:) / int2ext(:)
+    zdq(ilon,:,2) = dm3as(:) / int2ext(:)
+    zdq(ilon,:,3) = dm0af(:) / int2ext(:)
+    zdq(ilon,:,4) = dm3af(:) / int2ext(:)
+    
+    if (callclouds) then ! if call clouds
+      zdq(ilon,:,5) = dm0n(:) / int2ext(:)
+      zdq(ilon,:,6) = dm3n(:) / int2ext(:)
+      do i=1,nices
+        zdq(ilon,:,6+i) = dm3i(:,nices) / int2ext(:)
+        zdq(ilon,:,ices_indx(i)) = dgazs(:,i)  / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
+        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
+      enddo
+    endif
+
+  END DO ! loop on ilon
+
+END SUBROUTINE calmufi
Index: trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 1795)
@@ -5,10 +5,6 @@
 
       ! Main directory: 'datadir':
-      ! Default for Berserker @ UChicago:
-!      character(len=300) :: datadir='/home/rwordsworth/datagcm'
-      ! Default for Gnome Idataplex:
-!      character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
       ! Default for LMD machines:
-      character(len=300),save :: datadir='/u/lmdz/WWW/planets/LMDZ.GENERIC/datagcm'
+      character(len=300),save :: datadir='datagcm'
 !$OMP THREADPRIVATE(datadir)
       
@@ -18,4 +14,7 @@
       character(len=12),parameter :: surfdir="surface_data"
       
+      ! Default directory for microphysics
+      character(LEN=300),save :: config_mufi ='datagcm/microphysics/config.cfg'
+
       end module datafile_mod
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1795)
@@ -10,5 +10,5 @@
 
   use radinc_h, only: ini_radinc_h
-  use datafile_mod, only: datadir
+  use datafile_mod, only: datadir,config_mufi
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
   use comgeomfi_h, only: totarea, totarea_planet
@@ -319,4 +319,62 @@
      call getin_p("ichim",ichim)
      write(*,*)" ichim = ",ichim
+
+! Microphysics
+
+     write(*,*) "Run with or without microphysics?"
+     callmufi=.false. ! default value
+     call getin_p("callmufi",callmufi)
+     write(*,*)" callmufi = ",callmufi
+
+     ! sanity check
+     if (callmufi.and.(.not.tracer)) then
+       print*,"You are running microphysics without tracer"
+       print*,"Please start again with tracer =.true."
+       stop
+     endif
+
+     write(*,*) "Compute clouds?"
+     callclouds=.false. ! default value
+     call getin_p("callclouds",callclouds)
+     write(*,*)" callclouds = ",callclouds
+
+     ! sanity check
+     if (callclouds.and.(.not.callmufi)) then
+       print*,"You are trying to make clouds without microphysics !"
+       print*,"Please start again with callmufi =.true."
+       stop
+     endif
+
+     write(*,*) "Fractal dimension ?"
+     df_mufi=2.0 ! default value
+     call getin_p("df_mufi",df_mufi)
+
+     write(*,*) "Monomer radius (m) ?"
+     rm_mufi=6.66e-08 ! default value
+     call getin_p("rm_mufi",rm_mufi)
+
+     write(*,*) "Aerosol density (kg.m-3)?"
+     rho_aer_mufi=1.e3 ! default value
+     call getin_p("rho_aer_mufi",rho_aer_mufi)
+
+     write(*,*) "Pressure level of aer. production (Pa) ?"
+     p_prod=1.0 ! default value
+     call getin_p("p_prod",p_prod)
+     
+     write(*,*) "Aerosol production rate (kg.m-2.s-1) ?"
+     tx_prod=3.5e-13 ! default value
+     call getin_p("tx_prod",tx_prod)
+
+     write(*,*) "Equivalent radius production (m) ?"
+     rc_prod=2.0e-8 ! default value
+     call getin_p("rc_prod",rc_prod)
+
+     write(*,*) "Radius of air (nitrogen) molecule (m) ?"
+     air_rad=1.75e-10 ! default value
+     call getin_p("air_rad",air_rad)
+
+     write(*,*) "Path to microphys. config file ?"
+     config_mufi='datagcm/microphysics/config.cfg' ! default value
+     call getin_p("config_mufi",config_mufi)
 
 ! Soil model
Index: trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 1795)
+++ trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 1795)
@@ -0,0 +1,88 @@
+subroutine inimufi(nq,ptimestep)
+
+  use mmp_gcm
+  use callkeys_mod, only : callclouds, df_mufi, &
+        rm_mufi, rho_aer_mufi, p_prod, tx_prod, rc_prod, air_rad
+  use tracer_h
+  use comcstfi_mod, only : g, rad, mugaz
+  use datafile_mod
+
+  implicit none
+
+
+  !============================================================================
+  !
+  !     Purpose
+  !     -------
+  !     This routine call mm_initialize which perform the global initialization
+  !     for the microphysical YAMMS model.
+  !     It also performs some sanity checks on microphysical tracer names. 
+  !
+  !     Authors
+  !     -------
+  !     J. Vatant d'Ollone, J.Burgalat, S.Lebonnois (09/2017)
+  !
+  !============================================================================
+
+
+  !--------------------------
+  ! 0. Declarations
+  ! -------------------------
+
+  real, intent(in)                   :: ptimestep    ! Timestep (s)
+
+  integer, intent(in)                :: nq           ! Total number of tracers 
+
+  integer :: i,idx
+  character(len=20), dimension(4), parameter :: aernames = &
+     (/"mu_m0as              ", "mu_m3as              ",   &
+       "mu_m0af              ", "mu_m3af              "/)
+  CHARACTER(len=20), DIMENSION(2), PARAMETER :: ccnnames = &
+     (/"mu_m0n               ", "mu_m3n               "/)
+  logical :: err
+
+
+  !----------------------------------------------------
+  ! 1. Call microphysical model global initialization
+  ! ---------------------------------------------------
+
+   call mmp_initialize(ptimestep,df_mufi,rm_mufi,rho_aer_mufi,p_prod,tx_prod,rc_prod, &
+        rad,g,air_rad,mugaz,callclouds,config_mufi)
+
+   ! -------------------------
+   ! 2. Check names of tracers
+   ! -------------------------
+   err = .false.
+   DO i=1,size(aernames)
+     idx = indexOfTracer(TRIM(aernames(i)),.false.)
+     IF (idx <= 0) THEN
+       WRITE(*,*) "inimufi: '"//TRIM(aernames(i))//"' not found in tracers table."
+       err = .true.
+     ENDIF
+   ENDDO
+   IF (callclouds) THEN
+     DO i=1,size(ccnnames)
+       idx = indexOfTracer(TRIM(ccnnames(i)),.false.)
+       IF (idx <= 0) THEN
+         WRITE(*,*) "inimufi: '"//TRIM(ccnnames(i))//"' not found in tracers table."
+         err = .true.
+       ENDIF
+     ENDDO
+     ALLOCATE(ices_indx(mm_nesp))
+     ices_indx(:) = -1
+     DO i=1,mm_nesp 
+       idx = indexOfTracer("mu_m3"//TRIM(mm_spcname(i)),.false.)
+       IF (idx <= 0) THEN
+         WRITE(*,*) "inimufi: 'mu_m3"//TRIM(mm_spcname(i))//"' not found in tracers table."
+         err = .true.
+       ELSE
+         ices_indx(i) = idx
+       ENDIF
+     ENDDO
+     IF (err) THEN
+       WRITE(*,*) "You loose.... Try again"
+       STOP
+     ENDIF
+   ENDIF
+
+end subroutine inimufi
Index: trunk/LMDZ.TITAN/libf/phytitan/initracer.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/initracer.F	(revision 1792)
+++ 	(revision )
@@ -1,357 +1,0 @@
-      SUBROUTINE initracer(ngrid,nq,nametrac)
-
-      use surfdat_h
-      USE tracer_h
-      IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Initialization related to tracer 
-c   (chemical species, ice...)
-c
-c   Name of the tracer
-c
-c   Test of dimension :
-c   Initialize COMMON tracer in tracer.h, using tracer names provided
-c   by the argument nametrac
-c
-c   author: F.Forget
-c   ------
-c            Ehouarn Millour (oct. 2008) identify tracers by their names
-c            Jan Vatant d'Ollone (fev. 2017) chemical species for Titan
-c=======================================================================
-
-      integer :: ngrid,nq
-
-!      real qsurf(ngrid,nq)       ! tracer on surface (e.g.  kg.m-2)
-      character(len=20) :: txt ! to store some text
-      integer iq,ig,count
-!      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
-
-      character*20 nametrac(nq)   ! name of the tracer from dynamics
-
-
-c------------------------------------------------
-c  rho_q(nq)       ! tracer densities (kg.m-3)
-c------------------------------------------------
-       nqtot=nq
-       !! we allocate once for all arrays in common in tracer_h.F90
-       !! (supposedly those are not used before call to initracer)
-       IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
-       ALLOCATE(mmol(nq))
-       ALLOCATE(rho_q(nq))
-
-! Initialization: get tracer names from the dynamics and check if we are
-!                 using 'old' tracer convention ('q01',q02',...)
-!                 or new convention (full tracer names)
-      ! check if tracers have 'old' names
-
-! copy tracer names from dynamics
-        do iq=1,nq
-          noms(iq)=nametrac(iq)
-        enddo
-
-! Identify tracers by their names: (and set corresponding values of mmol)
-      ! 0. initialize tracer indexes to zero:
-      ! NB: igcm_* indexes are commons in 'tracer.h'
-      
-      igcm_h       =0
-      igcm_h2      =0
-      igcm_ch      =0
-      igcm_ch2s    =0
-      igcm_ch2     =0
-      igcm_ch3     =0
-      igcm_ch4     =0
-      igcm_c2      =0
-      igcm_c2h     =0
-      igcm_c2h2    =0
-      igcm_c2h3    =0
-      igcm_c2h4    =0
-      igcm_c2h5    =0
-      igcm_c2h6    =0
-      igcm_c3h3    =0
-      igcm_c3h5    =0
-      igcm_c3h6    =0
-      igcm_c3h7    =0
-      igcm_c4h     =0
-      igcm_c4h3    =0
-      igcm_c4h4    =0
-      igcm_c4h2s   =0
-      igcm_ch2cch2 =0
-      igcm_ch3cch  =0
-      igcm_c3h8    =0
-      igcm_c4h2    =0
-      igcm_c4h6    =0
-      igcm_c4h10   =0
-      igcm_ac6h6   =0
-      igcm_c3h2    =0
-      igcm_c4h5    =0
-      igcm_ac6h5   =0
-      igcm_n2      =0
-      igcm_n4s     =0
-      igcm_cn      =0
-      igcm_hcn     =0
-      igcm_h2cn    =0
-      igcm_chcn    =0
-      igcm_ch2cn   =0
-      igcm_ch3cn   =0
-      igcm_c3n     =0
-      igcm_hc3n    =0
-      igcm_nccn    =0
-      igcm_c4n2    =0      
-      
-      write(*,*) 'initracer: noms() ', noms
-
-
-      ! 1. find chemistry tracers
-      count = 0.
-      
-      do iq=1,nq
-      
-        if (noms(iq).eq."h") then
-          igcm_h=iq
-          mmol(igcm_h)=1.01
-          count=count+1
-        endif
-        if (noms(iq).eq."h2") then
-          igcm_h2=iq
-          mmol(igcm_h2)=2.0158
-          count=count+1
-        endif
-        if (noms(iq).eq."ch") then
-          igcm_ch=iq
-          mmol(igcm_ch)=13.02
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2s") then
-          igcm_ch2s=iq
-          mmol(igcm_ch2s)=14.03
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2") then
-          igcm_ch2=iq
-          mmol(igcm_ch2)=14.03
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3") then
-          igcm_ch3=iq
-          mmol(igcm_ch3)=15.03
-          count=count+1
-        endif
-        if (noms(iq).eq."ch4") then
-          igcm_ch4=iq
-          mmol(igcm_ch4)=16.04
-          count=count+1
-        endif
-        if (noms(iq).eq."c2") then
-          igcm_c2=iq
-          mmol(igcm_c2)=24.02
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h") then
-          igcm_c2h=iq
-          mmol(igcm_c2h)=25.03
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h2") then
-          igcm_c2h2=iq
-          mmol(igcm_c2h2)=26.04
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h3") then
-          igcm_c2h3=iq
-          mmol(igcm_c2h3)=27.05
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h4") then
-          igcm_c2h4=iq
-          mmol(igcm_c2h4)=28.05
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h5") then
-          igcm_c2h5=iq
-          mmol(igcm_c2h5)=29.06
-          count=count+1
-        endif
-        if (noms(iq).eq."c2h6") then
-          igcm_c2h6=iq
-          mmol(igcm_c2h6)=30.07
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h3") then
-          igcm_c3h3=iq
-          mmol(igcm_c3h3)=39.06
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h5") then
-          igcm_c3h5=iq
-          mmol(igcm_c3h5)=41.07
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h6") then
-          igcm_c3h6=iq
-          mmol(igcm_c3h6)=42.08
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h7") then
-          igcm_c3h7=iq
-          mmol(igcm_c3h7)=43.09
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h") then
-          igcm_c4h=iq
-          mmol(igcm_c4h)=49.05
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h3") then
-          igcm_c4h3=iq
-          mmol(igcm_c4h3)=51.07
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h4") then
-          igcm_c4h4=iq
-          mmol(igcm_c4h4)=52.08
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h2s") then
-          igcm_c4h2s=iq
-          mmol(igcm_c4h2s)=50.06
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2cch2") then
-          igcm_ch2cch2=iq
-          mmol(igcm_ch2cch2)=40.07
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3cch") then
-          igcm_ch3cch=iq
-          mmol(igcm_ch3cch)=40.07
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h8") then
-          igcm_c3h8=iq
-          mmol(igcm_c3h8)=44.11
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h2") then
-          igcm_c4h2=iq
-          mmol(igcm_c4h2)=50.06
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h6") then
-          igcm_c4h6=iq
-          mmol(igcm_c4h6)=54.09
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h10") then
-          igcm_c4h10=iq
-          mmol(igcm_c4h10)=58.13
-          count=count+1
-        endif
-        if (noms(iq).eq."ac6h6") then
-          igcm_ac6h6=iq
-          mmol(igcm_ac6h6)=78.1136
-          count=count+1
-        endif
-        if (noms(iq).eq."c3h2") then
-          igcm_c3h2=iq
-          mmol(igcm_c3h2)=38.05
-          count=count+1
-        endif
-        if (noms(iq).eq."c4h5") then
-          igcm_c4h5=iq
-          mmol(igcm_c4h5)=53.07
-          count=count+1
-        endif
-        if (noms(iq).eq."ac6h5") then
-          igcm_ac6h5=iq
-          mmol(igcm_ac6h5)=77.1136
-          count=count+1
-        endif
-        if (noms(iq).eq."n2") then
-          igcm_n2=iq
-          mmol(igcm_n2)=28.0134
-          count=count+1
-        endif
-        if (noms(iq).eq."n4s") then
-          igcm_n4s=iq
-          mmol(igcm_n4s)=14.01
-          count=count+1
-        endif
-        if (noms(iq).eq."cn") then
-          igcm_cn=iq
-          mmol(igcm_cn)=26.02
-          count=count+1
-        endif
-        if (noms(iq).eq."hcn") then
-          igcm_hcn=iq
-          mmol(igcm_hcn)=27.04
-          count=count+1
-        endif
-        if (noms(iq).eq."h2cn") then
-          igcm_h2cn=iq
-          mmol(igcm_h2cn)=28.05
-          count=count+1
-        endif
-        if (noms(iq).eq."chcn") then
-          igcm_chcn=iq
-          mmol(igcm_chcn)=39.05
-          count=count+1
-        endif
-        if (noms(iq).eq."ch2cn") then
-          igcm_ch2cn=iq
-          mmol(igcm_cn)=40.04
-          count=count+1
-        endif
-        if (noms(iq).eq."ch3cn") then
-          igcm_ch3cn=iq
-          mmol(igcm_ch3cn)=41.05
-          count=count+1
-        endif
-        if (noms(iq).eq."c3n") then
-          igcm_c3n=iq
-          mmol(igcm_c3n)=50.04
-          count=count+1
-        endif
-        if (noms(iq).eq."hc3n") then
-          igcm_hc3n=iq
-          mmol(igcm_hc3n)=51.05
-          count=count+1
-        endif
-        if (noms(iq).eq."nccn") then
-          igcm_nccn=iq
-          mmol(igcm_nccn)=52.04
-          count=count+1
-        endif
-        if (noms(iq).eq."c4n2") then
-          igcm_c4n2=iq
-          mmol(igcm_c4n2)=76.1
-          count=count+1
-        endif
-        
-
-      enddo ! of do iq=1,nq
- 
-      ! check that we identified all tracers:
-      if (count.ne.nq) then
-        write(*,*) "initracer: found only ",count," tracers"
-        write(*,*) "               expected ",nq
-        do iq=1,count
-          write(*,*)'      ',iq,' ',trim(noms(iq))
-        enddo
-!        stop
-      else
-        write(*,*) "initracer: found all expected tracers, namely:"
-        do iq=1,nq
-          write(*,*)'      ',iq,' ',trim(noms(iq))
-        enddo
-      endif
-
-
-c------------------------------------------------------------
-c     Initialisation tracers ....
-c------------------------------------------------------------
-      call zerophys(nq,rho_q)
-
-
-      end
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1795)
@@ -22,5 +22,5 @@
       use geometry_mod, only: latitude, longitude, cell_area
       USE comgeomfi_h, only: totarea, totarea_planet
-      USE tracer_h, only: noms, mmol
+      USE tracer_h
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0
@@ -32,5 +32,5 @@
       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
       use time_phylmdz_mod, only: daysec
-      use logic_mod, only: moyzon_ch
+      use logic_mod, only: moyzon_ch, moyzon_mu
       use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
                             zplaybar, zzlevbar, zzlaybar, ztfibar
@@ -73,4 +73,5 @@
 !      V. Tracers
 !         V.1. Chemistry
+!         V.2. Microphysics
 !         V.3. Updates (pressure variations, surface budget).
 !         V.4. Surface Tracer Update.
@@ -159,4 +160,5 @@
 !   -------
 
+
       integer,intent(in) :: ngrid             ! Number of atmospheric columns.
       integer,intent(in) :: nlayer            ! Number of atmospheric layers.
@@ -277,5 +279,7 @@
       
       real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
-      real zdqmph(ngrid,nlayer,nq)    ! Microphysical tendency ( condensation routine only for now, no microphysical routines ).
+      real zdqcond(ngrid,nlayer,nq)   ! Condensation tendency ( chemistry routine ).
+      
+      real zdqmufi(ngrid,nlayer,nq)   ! Microphysical tendency.
                         
       ! For Winds : (m/s/s)
@@ -349,5 +353,7 @@
 ! ----------------------------------------------------------------------------
      
-      integer,parameter  :: nmicro=0 ! Temporary ! To be put in start/def
+      character*10,dimension(:),allocatable,save :: nomqy
+
+!$OMP THREADPRIVATE(nomqy)
 
       real ctimestep ! Chemistry timestep (s)
@@ -356,22 +362,15 @@
       real temp_eq(nlayer), press_eq(nlayer)
       real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
-!      real zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
       real ztemp(ngrid,nlayer)
       
-      real ychim(ngrid,nlayer,nq-nmicro)   
-
-      real rat_mmol(nq) ! Molar fraction ratio
-
-      ! 2D vmr tendencies ( chemistry and condensation )
-      real,dimension(:,:,:),allocatable,save :: dycchi
-      ! Must be saved since chemistry is not called every step
-      
-      real dycmph(ngrid,nlayer,nq-nmicro)       
-      ! ----------------------------------------------------------
-         
-      real,dimension(:,:),allocatable,save       :: qysat
-      
-      character*10,dimension(:),allocatable,save :: nomqy
-!$OMP THREADPRIVATE(dycchi,qysat,nomqy)
+      real , allocatable, dimension(:,:,:),save :: ychim
+
+      ! 2D vmr tendencies ( chemistry and condensation ) for chem. tracers
+      real,dimension(:,:,:),allocatable,save :: dycchi ! Saved since chemistry is not called every step
+      real dyccond(ngrid,nlayer,nq)       
+         
+      real,dimension(:,:),allocatable,save   :: qysat
+      
+!$OMP THREADPRIVATE(ychim,dycchi,qysat)
 
       real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank
@@ -379,4 +378,19 @@
 
 !-----------------------------------------------------------------------------
+    ! Interface to calmufi 
+    !   --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module 
+    !       (to have an explicit generated by the compiler).
+    !   Or one can put calmufi in MMP_GCM module (in muphytitan).
+    INTERFACE
+      SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev  !! Altitude levels (m).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K).  
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.kg^{-1}}\)).
+      END SUBROUTINE calmufi
+    END INTERFACE
       
 !==================================================================================================
@@ -390,4 +404,9 @@
 ! --------------------------------
       if (firstcall) then
+        ! initracer:
+        !   initialize nmicro, 
+        call initracer2(nq,nametrac)
+
+        ! ----------------------------------------------------------------------------
 
         ! Allocate saved arrays.
@@ -423,5 +442,5 @@
         ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
         ALLOCATE(qysat(nlayer,nq-nmicro))
-        ALLOCATE(nomqy(nq-nmicro+1))
+        ALLOCATE(nomqy(nq-nmicro+1)) ! +1 because of hv
         
         ! This is defined in comsaison_h
@@ -439,13 +458,37 @@
          zdtlw(:,:) = 0.0
 
-      
-!        Initialize tracer names, indexes and properties.
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
-         if (tracer) then
-            call initracer(ngrid,nq,nametrac)
-         endif
-
-         rat_mmol(:) = mmol(:) / mugaz
+!        Initialize names, timestep and saturation profiles for chemistry
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if ( callchim .and. (nq.gt.nmicro) ) then
+
+            ctimestep = ptimestep*REAL(ichim)
+
+            do iq=nmicro+1,nq
+               nomqy(iq-nmicro) = nametrac(iq)
+            enddo
+
+            nomqy(nq-nmicro+1) = "HV"
+            
+            ! qysat is taken at the equator ( small variations of t,p)
+            temp_eq(:)  = tmoy(:)
+            press_eq(:) = playmoy(:)/100. ! en mbar
+            
+            call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
+          
+            zdqchi(:,:,:)  = 0.0
+            zdqcond(:,:,:) = 0.0
+
+         endif
+
+!        Initialize microphysics.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~
+
+         IF ( callmufi ) THEN
+
+           call inimufi(nq,ptimestep)
+
+         ENDIF
+
 
 !        Read 'startfi.nc' file.
@@ -485,5 +528,4 @@
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 
-
          if(tlocked)then
             print*,'Planet is tidally locked at resonance n=',nres
@@ -533,30 +575,18 @@
          endif
 
+         ! Sanity check for microphysics
+         if ( ((.not.moyzon_mu).and.(callmufi)) ) then
+            print *, "moyzon_mu=",moyzon_mu," and callmufi=",callmufi
+            print *, "Please activate zonal mean to run microphysics (for now) !"
+            stop
+         endif
+
          ! Sanity check for chemistry
-         if ( ((moyzon_ch).and.(.not.callchim)) .or. ((.not.moyzon_ch).and.(callchim)) ) then
+         if ( (.not.moyzon_ch) .and. (callchim) ) then
             print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
-            print *, "This is not compatible..."
+            print *, "Please activate zonal mean to run chemistry !"
             stop
          endif
 
-         ! Initialize names, timestep and saturation profiles for chemistry
-
-         if ( callchim .and. (nq.gt.nmicro) ) then
-
-            ctimestep = ptimestep*REAL(ichim)
-
-            do iq=nmicro+1,nq
-               nomqy(iq-nmicro) = nametrac(iq)
-            enddo
-
-            nomqy(nq-nmicro+1) = "HV"
-            
-            ! qysat is taken at the equator ( small variations of t,p)
-            temp_eq(:)  = tmoy(:)
-            press_eq(:) = playmoy(:)/100. ! en mbar
-            
-            call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
-          
-         endif
          
          ! XIOS outputs
@@ -647,5 +677,5 @@
       ! -------------------------------Taken from old Titan --------------------------
       ! zonal averages needed
-      if (moyzon_ch) then
+      if (moyzon_ch .or. moyzon_mu) then
          
          zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
@@ -681,6 +711,6 @@
 
       endif  ! moyzon
+
       ! -------------------------------------------------------------------------------------
-      
       ! Compute potential temperature
       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
@@ -981,7 +1011,8 @@
   !   V.1. Chemistry
   ! -------------------------
+
          if (callchim) then
 
-            ! Utilisation de la moyenne zonale dans calchim
+            ! Using zonal mean for calchim
             zplev(:,:) = zplevbar(:,:)
             zplay(:,:) = zplaybar(:,:)
@@ -990,4 +1021,6 @@
             ztemp(:,:) = ztfibar(:,:)
 
+            allocate(ychim(ngrid,nlayer,nq-nmicro))
+
             if (nq.gt.nmicro) then
                do iq = nmicro+1,nq
@@ -997,9 +1030,9 @@
 
             ! Condensation tendency after the transport
-            do iq=1,nq-nmicro
+            do iq=nmicro+1,nq
                do l=1,nlayer
                   do ig=1,ngrid
-                     if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
-                        dycmph(ig,l,nmicro+iq)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
+                     if ( ychim(ig,l,iq-nmicro).gt.qysat(l,iq-nmicro) ) then
+                        dyccond(ig,l,iq)= ( -ychim(ig,l,iq-nmicro)+qysat(l,iq-nmicro) ) / ptimestep
                      endif
                   enddo
@@ -1007,8 +1040,9 @@
             enddo
 
-            if( mod(icount-1,ichim).eq.0.) then
+            if ( callclouds ) dyccond(:,:,ices_indx) = 0.0 ! Condensation will be calculated in the cloud microphysics
+
+            if( mod(icount-1,ichim).eq.0. ) then
                
-               print *, "On passe dans la chimie..."
-               
+               print *, "We enter in the chemistry ..."
                call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
                            ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
@@ -1018,17 +1052,43 @@
             endif
             
+            ! TEMPORARY COMMENT
+            ! These conversions as well as 2D/1D should be put in phytrac 
+            ! ( GCM-friendly tracers and tendencies -> format for photochem and microphys ) 
+
             if (nq.gt.nmicro) then
-               ! We convert tendencies back to mass mixing ratio
+               ! We convert tendencies to mass mixing ratio
                do iq=nmicro+1,nq
-                  zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro) / rat_mmol(iq)
-                  zdqmph(:,:,iq) = dycmph(:,:,iq-nmicro) / rat_mmol(iq)
+                  zdqchi(:,:,iq)  = dycchi(:,:,iq-nmicro)  / rat_mmol(iq)
+                  zdqcond(:,:,iq) = dyccond(:,:,iq) / rat_mmol(iq)
                enddo
 
                pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
-                    zdqchi(1:ngrid,1:nlayer,1:nq) + zdqmph(1:ngrid,1:nlayer,1:nq)
+                    zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
                
             endif
             
-         endif
+         endif ! end of 'callchim'
+
+  ! -------------------
+  !   V.2 Microphysics
+  ! -------------------
+ 
+         if (callmufi) then
+
+            ! Using zonal mean for microphysics
+            zplev(:,:) = zplevbar(:,:)
+            zplay(:,:) = zplaybar(:,:)
+            zzlev(:,:) = zzlevbar(:,:)
+            zzlay(:,:) = zzlaybar(:,:)
+            ztemp(:,:) = ztfibar(:,:)
+
+            ! Inside this routine we will split 2D->1D, intensive->extensive and separate different types of tracers
+            ! Should be put in phytrac
+
+            call calmufi(zplev,zzlev,zplay,zzlay,ztemp,pq,zdqmufi)
+
+            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmufi(1:ngrid,1:nlayer,1:nq)   
+
+         endif ! end of 'callmufi'
       
   ! ---------------
@@ -1346,9 +1406,7 @@
       
          !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
-         !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
          call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
          call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
          call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-         call writediagfi(ngrid,"shad","rings"," ", 2, fract)
          
 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
@@ -1403,22 +1461,12 @@
       if (tracer) then
       
-         do iq=1,nq
-            call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-                             'kg m^-2',2,qsurf_hist(1,iq) )
-                          
-!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-!                          'kg m^-2',2,qsurf(1,iq) )                          
-
-         enddo ! end of 'nq' loop
-         
+         if (callmufi) then ! For now we assume an given order for tracers !
+            call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'kg/kg',3,zq(:,:,1))
+            call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'kg/kg',3,zq(:,:,2))
+            call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'kg/kg',3,zq(:,:,3))
+            call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'kg/kg',3,zq(:,:,4))
+         endif ! end of 'callmufi'
+
        endif ! end of 'tracer'
-
-
-      ! Output spectrum.
-      if(specOLR.and.corrk)then      
-         call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
-         call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
-      endif
 
 ! XIOS outputs
Index: trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1792)
+++ trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1795)
@@ -1,76 +1,311 @@
 
-       module tracer_h
-
-       implicit none
-
-! nqtot : total number of tracers
-       INTEGER, SAVE :: nqtot
-!$OMP THREADPRIVATE(nqtot)
-
-       character*20, DIMENSION(:), ALLOCATABLE :: noms   ! name of the tracer
-       real, DIMENSION(:), ALLOCATABLE :: mmol     ! mole mass of tracer (g/mol-1) 
-       real, DIMENSION(:), ALLOCATABLE :: rho_q    ! tracer densities (kg.m-3)
-
-!$OMP THREADPRIVATE(noms,mmol,rho_q)
-
-! tracer indexes: these are initialized in initracer and should be 0 if the
-!                 corresponding tracer does not exist
-
-      ! chemistry:
-      
-      integer :: igcm_h
-      integer :: igcm_h2
-      integer :: igcm_ch
-      integer :: igcm_ch2s
-      integer :: igcm_ch2
-      integer :: igcm_ch3
-      integer :: igcm_ch4
-      integer :: igcm_c2
-      integer :: igcm_c2h
-      integer :: igcm_c2h2
-      integer :: igcm_c2h3
-      integer :: igcm_c2h4
-      integer :: igcm_c2h5
-      integer :: igcm_c2h6
-      integer :: igcm_c3h3
-      integer :: igcm_c3h5
-      integer :: igcm_c3h6
-      integer :: igcm_c3h7
-      integer :: igcm_c4h
-      integer :: igcm_c4h3
-      integer :: igcm_c4h4
-      integer :: igcm_c4h2s
-      integer :: igcm_ch2cch2
-      integer :: igcm_ch3cch
-      integer :: igcm_c3h8
-      integer :: igcm_c4h2
-      integer :: igcm_c4h6
-      integer :: igcm_c4h10
-      integer :: igcm_ac6h6
-      integer :: igcm_c3h2
-      integer :: igcm_c4h5
-      integer :: igcm_ac6h5
-      integer :: igcm_n2
-      integer :: igcm_n4s
-      integer :: igcm_cn
-      integer :: igcm_hcn
-      integer :: igcm_h2cn
-      integer :: igcm_chcn
-      integer :: igcm_ch2cn
-      integer :: igcm_ch3cn
-      integer :: igcm_c3n
-      integer :: igcm_hc3n
-      integer :: igcm_nccn
-      integer :: igcm_c4n2 
-      
-      
-!$OMP THREADPRIVATE(igcm_h,igcm_h2,igcm_ch,igcm_ch2s,igcm_ch2,igcm_ch3,igcm_ch4,   &
-      !$OMP igcm_c2,igcm_c2h,igcm_c2h2,igcm_c2h3,igcm_c2h4,igcm_c2h5,igcm_c2h6,    &
-      !$OMP igcm_c3h3,igcm_c3h5,igcm_c3h6,igcm_c3h7,igcm_c4h,igcm_c4h3,igcm_c4h4,  &
-      !$OMP igcm_c4h2s,igcm_ch2cch2,igcm_ch3cch,igcm_c3h8,igcm_c4h2,igcm_c4h6,     &
-      !$OMP igcm_c4h10,igcm_ac6h6,igcm_c3h2,igcm_c4h5,igcm_ac6h5,igcm_n2,igcm_n4s, &
-      !$OMP igcm_cn,igcm_hcn,igcm_h2cn,igcm_chcn,igcm_ch2cn,igcm_ch3cn,igcm_c3n,   &
-      !$OMP igcm_hc3n,igcm_nccn,igcm_c4n2)
-
-       end module tracer_h
-
+MODULE tracer_h
+  !! Stores data related to physics tracers.
+  !!
+  !! The module stores public global variables related to the number of tracers
+  !! available in the physics and their kind:
+  !!
+  !! Currently, tracers can be used either for chemistry process (nchimi) or
+  !! microphysics (nmicro).
+  !!
+  !! The subroutine "initracer2" initializes and performs sanity check of
+  !! the tracer definitions given in traceur.def and the required tracers in physics
+  !! (based on the run parameters).
+  !!
+  !! The module provides additional methods: 
+  !! 
+  !!   - indexOfTracer : search for the index of a tracer in the global table (tracers_h:noms) by name.
+  !!   - nameOfTracer  : get the name of tracer from a given index (of the global table).
+  !!   - dumpTracers   : print the names of all tracers indexes given in argument.
+  !!
+  IMPLICIT NONE
+
+  INTEGER, SAVE :: nqtot  = 0 !! Total number of tracers
+  INTEGER, SAVE :: nmicro = 0 !! Number of microphysics tracers.
+  INTEGER, SAVE :: nice   = 0 !! Number of microphysics ice tracers (subset of nmicro).
+  INTEGER, SAVE :: nchimi = 0 !! Number of chemical (gaz species) tracers.
+  !$OMP THREADPRIVATE(nqtot,nmicro,nice,nchimi)
+
+  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: chimi_indx  !! Indexes of all chemical species tracers 
+  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: micro_indx  !! Indexes of all microphysical tracers
+  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ices_indx   !! Indexes of all ice microphysical tracers
+
+  CHARACTER(len=20), DIMENSION(:), ALLOCATABLE, SAVE :: noms      !! name of the tracer
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: mmol      !! mole mass of tracer(g/mol-1)
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: rat_mmol  !! molar mass ratio 
+  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: rho_q     !! tracer densities (kg.m-3)
+  !$OMP THREADPRIVATE(noms,mmol,rat_mmol,rho_q)
+
+
+
+  ! tracer indexes: these are initialized in initracer and should be 0 if the
+  !                 corresponding tracer does not exist
+
+
+  CONTAINS
+
+  SUBROUTINE initracer2(nq,nametrac,talk_to_me)
+    !! Initialize tracer names and attributes.
+    !!
+    !! The method initializes the list of tracer names used in the physics from their
+    !! dynamics counterpart.
+    !!
+    !! In addition, it initializes arrays of indexes for the different sub-processs of the physics:
+    !!
+    !!   - tracers_h:micro_indxs, the array of tracers indexes used for the microphysics.
+    !!   - tracers_h:chimi_indxs, the array of tracers indexes used for the chemistry.
+    !!
+    !! The method also initializes the molar mass array (tracers_h:mmol) for the chemistry and the
+    !! molar mass ratio (tracers_h:rat_mmol).
+    !!
+    !! @note
+    !! Strict checking of chemical species name is performed here if the chemistry is activated
+    !! (see callchim variable). All the values of 'cnames' must be found in the tracers names
+    !! related to chemistry.
+    !! @note
+    !! Tests are more permissive for the microphysics and is only based on the mimimum number of
+    !! tracers expected. Strict name checking is performed in inimufi.
+    USE callkeys_mod
+    USE comcstfi_mod, only: mugaz
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                          :: nq         !! Total number of tracers (fixed at compile time)
+    character(len=20), DIMENSION(nq), INTENT(in) :: nametrac   !! name of the tracer from dynamics (from 'traceurs.def')
+    LOGICAL, INTENT(in), OPTIONAL                :: talk_to_me !! Enable verbose mode.
+
+    LOGICAL                                      :: verb,found
+    CHARACTER(len=20)                            :: str
+    !! Hard-coded chemical species for Titan chemistry
+    CHARACTER(len=10), DIMENSION(44), PARAMETER  :: cnames = &
+      (/"H         ", "H2        ", "CH        ", "CH2s      ", "CH2       ", "CH3       ", &
+        "CH4       ", "C2        ", "C2H       ", "C2H2      ", "C2H3      ", "C2H4      ", &
+        "C2H5      ", "C2H6      ", "C3H3      ", "C3H5      ", "C3H6      ", "C3H7      ", & 
+        "C4H       ", "C4H3      ", "C4H4      ", "C4H2s     ", "CH2CCH2   ", "CH3CCH    ", &
+        "C3H8      ", "C4H2      ", "C4H6      ", "C4H10     ", "AC6H6     ", "C3H2      ", &
+        "C4H5      ", "AC6H5     ", "N2        ", "N4S       ", "CN        ", "HCN       ", &
+        "H2CN      ", "CHCN      ", "CH2CN     ", "CH3CN     ", "C3N       ", "HC3N      ", &
+        "NCCN      ", "C4N2      "/)
+    !! Hard-coded chemical species molar mass (g.mol-1), shares the same indexing than cnames.
+    REAL, DIMENSION(44), PARAMETER               :: cmmol = (/ &
+        1.01   , 2.0158, 13.02, 14.03, 14.03, 15.03, 16.04  , 24.02, 25.03, 26.04  , 27.05  , &
+        28.05  , 29.06 , 30.07, 39.06, 41.07, 42.08, 43.09  , 49.05, 51.07, 52.08  , 50.06  , &
+        40.07  , 40.07 , 44.11, 50.06, 54.09, 58.13, 78.1136, 38.05, 53.07, 77.1136, 28.0134, &
+        14.01  , 26.02 , 27.04, 28.05, 39.05, 40.04, 41.05  , 50.04, 51.05, 52.04  , 76.1   /)
+
+    INTEGER :: i,j,n
+
+    verb = .true. ; IF (PRESENT(talk_to_me)) verb = talk_to_me
+
+    ! nqtot could be used everywhere in the physic :)
+    nqtot=nq
+
+    IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
+    noms(:)=nametrac(:)
+
+    ALLOCATE(rho_q(nq)) ! Defined for all tracers, currently initialized to 0.0
+    rho_q(:) = 0.0
+
+    ALLOCATE(mmol(nq),rat_mmol(nq))  ! Defined for all tracers, (actually) initialized only for chemical tracers
+    mmol(:)  = 0.0
+    rat_mmol(:) = 1.0
+
+    ! Compute number of microphysics tracers:
+    ! By convention they all have the prefix "mu_" (case sensitive !)
+    nmicro = 0
+    IF (callmufi) THEN
+      DO i=1,nq
+        str = noms(i)
+        IF (str(1:3) == "mu_") nmicro = nmicro+1
+      ENDDO
+      ! Checking the expected number of tracers:
+      !   no cloud:  4 ; w cloud :  4 + 2 + (1+)
+      ! Note that we do not make assumptions on the number of chemical species for clouds, this
+      ! will be checked in inimufi.
+      IF (callclouds) THEN
+        IF (nmicro < 7) THEN
+          WRITE(*,'((a),I3,(a))') "initracer2:error: Inconsistent number of microphysical tracers &
+            &(expected at least 7 tracers,",nmicro," given)"
+          CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
+          STOP
+        ENDIF
+      ELSE IF (nmicro < 4) THEN
+          WRITE(*,'((a),I3,(a))') "initracer2:error: Inconsistent number of microphysical tracers &
+            &(expected at least 4 tracers,",nmicro," given)"
+          CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
+      ELSE IF  (nmicro > 4) THEN
+        WRITE(*,'(a)') "initracer2:info: I was expecting only four tracers, you gave me &
+          &more. I'll just pretend nothing happen !"
+      ENDIF
+      ! microphysics indexes share the same values than original tracname.
+      ALLOCATE(micro_indx(nmicro))
+      j = 1
+      DO i=1,nq
+        str = noms(i)
+        IF (str(1:3) == "mu_") THEN
+          micro_indx(j) = i 
+          j=j+1
+        ENDIF
+      ENDDO
+    ELSE
+      ALLOCATE(micro_indx(nmicro))
+    ENDIF
+
+    ! Compute number of chemical species:
+    !   simply assume that all other tracers ARE chemical species
+    nchimi = nqtot - nmicro
+
+    ! Titan chemistry requires exactly 44 tracers:
+    ! Test should be in callchim condition
+
+    IF (callchim) THEN
+      IF (nchimi < SIZE(cnames)) THEN
+        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
+        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
+      ENDIF
+      ALLOCATE(chimi_indx(nchimi))
+      n = 0 ! counter on chimi_indx 
+      DO j=1,SIZE(cnames)
+        found = .false.
+        DO i=1,nq
+          IF (TRIM(cnames(j)) == TRIM(noms(i))) THEN
+            n = n + 1
+            chimi_indx(n) = i
+            mmol(i) = cmmol(j)
+            rat_mmol(i) = cmmol(j)/mugaz
+            found = .true.
+            EXIT
+          ENDIF
+        ENDDO
+        IF (.NOT.found) THEN
+          WRITE(*,*) "initracer2:error: "//TRIM(cnames(j))//" is missing from tracers definition file."
+        ENDIF
+      ENDDO
+      IF (n < SIZE(cnames)) THEN
+        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
+        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
+      ENDIF
+    ELSE
+      ALLOCATE(chimi_indx(0))
+    ENDIF
+    IF (verb) THEN
+      IF (callmufi.OR.callchim) WRITE(*,*) "===== INITRACER2 SPEAKING ====="
+      IF (callmufi) THEN
+         WRITE(*,*) "Found ",nmicro, "microphysical tracers"
+        call dumpTracers(micro_indx)
+        WRITE(*,*) "-------------------------------"
+      ENDIF
+      IF (callchim) THEN
+        WRITE(*,*) "Found ",nchimi, "chemical tracers"
+        call dumpTracers(chimi_indx)
+        WRITE(*,*) "-------------------------------"
+      ENDIF
+    ENDIF
+
+  END SUBROUTINE initracer2
+
+  FUNCTION indexOfTracer(name, sensitivity) RESULT(idx)
+    !! Get the index of a tracer by name.
+    !!
+    !! The function searches in the global tracer table (tracer_h:noms)
+    !! for the given name and returns the first index matching "name".
+    !! 
+    !! If no name in the table matches the given one, -1 is returned !
+    !!
+    !! @warning
+    !! initracers must be called before any use of this function.
+    IMPLICIT NONE
+    CHARACTER(len=*), INTENT(in)  :: name         !! Name of the tracer to search.
+    LOGICAL, OPTIONAL, INTENT(in) :: sensitivity  !! Case sensitivity (true by default).
+    INTEGER :: idx                                !! Index of the first tracer matching name or -1 if not found.
+    LOGICAL                  :: zsens
+    INTEGER                  :: j
+    CHARACTER(len=LEN(name)) :: zname
+    zsens = .true. ; IF(PRESENT(sensitivity)) zsens = sensitivity
+    idx = -1
+    IF (.NOT.ALLOCATED(noms)) RETURN
+    IF (zsens) THEN
+      DO j=1,SIZE(noms) 
+        IF (TRIM(noms(j)) == TRIM(name)) THEN
+          idx = j ; RETURN 
+        ENDIF 
+      ENDDO
+    ELSE
+      zname = to_lower(name)
+      DO j=1,SIZE(noms) 
+        IF (TRIM(to_lower(noms(j))) == TRIM(zname)) THEN
+          idx = j ; RETURN 
+        ENDIF 
+      ENDDO
+    ENDIF
+
+    CONTAINS
+
+    FUNCTION to_lower(istr) RESULT(ostr)
+      !! Lower case conversion function.
+      IMPLICIT NONE
+      CHARACTER(len=*), INTENT(in) :: istr
+      CHARACTER(len=LEN(istr)) :: ostr
+      INTEGER :: i, ic
+      ostr = istr
+      DO i = 1, LEN_TRIM(istr)
+        ic = ICHAR(istr(i:i))
+        IF (ic >= 65 .AND. ic < 90) ostr(i:i) = char(ic + 32)
+      ENDDO 
+    END FUNCTION to_lower
+  END FUNCTION indexOfTracer
+
+  FUNCTION nameOfTracer(indx) RESULT(name)
+    !! Get the name of a tracer by index.
+    !!
+    !! The function searches in the global tracer table (tracer_h:noms)
+    !! and returns the name of the tracer at given index.
+    !!
+    !! If the index is out of range an empty string is returned.
+    !!
+    !! @warning
+    !! initracers must be called before any use of this function.
+    IMPLICIT NONE
+    INTEGER, INTENT(in) :: indx   !! Index of the tracer name to retrieve.
+    CHARACTER(len=20)   :: name   !! Name of the tracer at given index.
+    name = '' 
+    IF (.NOT.ALLOCATED(noms)) RETURN
+    IF (indx <= 0 .OR. indx > SIZE(noms)) RETURN
+    name = noms(indx)
+  END FUNCTION nameOfTracer
+
+  SUBROUTINE dumpTracers(indexes)
+    !! Print the names of the given list of tracers indexes.
+    INTEGER, DIMENSION(:), INTENT(in) :: indexes
+    INTEGER :: i,idx,nt
+    CHARACTER(len=:), ALLOCATABLE :: suffix
+    
+    IF (.NOT.ALLOCATED(noms)) THEN
+      WRITE(*,'(a)') "[tracers_h:dump_tracers] warning: 'noms' is not allocated, tracers_h:initracer2 has not be called yet"
+      RETURN
+    ENDIF
+    nt = size(noms)
+    WRITE(*,"(a)") "local -> global : name"
+    DO i=1,size(indexes)
+      idx = indexes(i)
+      IF (idx < 1 .OR. idx >= nt) THEN
+        ! WRITE(*,'((a),I3,(a),I3,(a))') "index out of range (",idx,"/",nt,")"
+        CYCLE
+      ENDIF
+      IF (ANY(chimi_indx == idx)) THEN
+        suffix = ' (chimi)'
+      ELSE IF (ANY(micro_indx == idx)) THEN
+        suffix = ' (micro'
+        IF (ALLOCATED(ices_indx)) THEN
+          IF (ANY(ices_indx == idx)) suffix=suffix//", ice"
+        ENDIF
+        suffix=suffix//")"
+      ELSE 
+        suffix=" ()"
+      ENDIF
+      WRITE(*,'(I5,(a),I6,(a))') i," -> ",idx ," : "//TRIM(noms(i))//suffix
+    ENDDO
+  END SUBROUTINE dumpTracers
+
+
+END MODULE tracer_h
+
