Index: /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90	(revision 1671)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90	(revision 1672)
@@ -372,6 +372,6 @@
     write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
     write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-  else if (planet_type.eq.'titan') then
-    jD_ref=1 ! Because we use old Titan starts (J.V.O 2016)
+!  else if (planet_type.eq.'titan') then
+!    jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
   else
 ! A voir pour Titan et Venus
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1671)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1672)
@@ -396,6 +396,6 @@
       write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
       write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
-      else if (planet_type.eq.'titan') then
-        jD_ref=1 ! Because we use old Titan starts (J.V.O 2016)
+!      else if (planet_type.eq.'titan') then
+!        jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
       else
 ! A voir pour Titan et Venus
Index: /trunk/LMDZ.TITAN/libf/chimtitan/disso.c
===================================================================
--- /trunk/LMDZ.TITAN/libf/chimtitan/disso.c	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/chimtitan/disso.c	(revision 1672)
@@ -179,5 +179,5 @@
    double f,**flux;
    double flact;
-   char   name[60],dir[24];
+   char   name[60],dir[27];
    FILE   *fp,*out;
 
@@ -186,7 +186,7 @@
 /* lecture des flux actiniques:
    - suppose que l'executable est dans $LMDGCM/RUN/xxx/
-   - et les moyennes dans $LMDGCM/INPUT/PHOT(NLAT)/Moy_(lat:1 a NLAT)
+   - et les moyennes dans $LMDGCM/datagcm/PHOT(NLAT)/Moy_(lat:1 a NLAT)
 */
-   strcpy( dir, "../../INPUT/PHOT" );
+   strcpy( dir, "../../datagcm/PHOT" );
    if( (*NLAT) < 10 )
    {
Index: /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1672)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1672)
@@ -0,0 +1,505 @@
+      SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad, &
+                         dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc,NLEV)
+      
+!--------------------------------------------------------------------
+!  
+!     Introduction d une routine chimique
+!
+!     Auteurs: S. Lebonnois,  01/2000 | 09/2003 
+!            adaptation pour Titan 3D: 02/2009
+!            adaptation pour // : 04/2013
+!            extension chimie jusqu a 1300km : 10/2013
+!      
+!              J. Vatant d'Ollone, 02/2017
+!               + adaptation pour le nouveau titan issu du generic
+! 
+!---------------------------------------------------------------------
+!
+      
+      use dimphy
+      
+      use datafile_mod, only: datadir
+      
+      use comcstfi_mod, only: rad, pi, kbol
+      use moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
+      use mod_grid_phy_lmdz, only: nbp_lat
+      
+      implicit none
+
+!   Parameters ( dans common_mod in old titan)
+!              + doivent etre en accord avec titan.h
+!   -------------------------------------------------
+
+      INTEGER,PARAMETER ::   NC   = 44     ! Nb de composes dans la chimie
+      INTEGER,PARAMETER ::   ND   = 54     ! Nb de photodissociations
+      INTEGER,PARAMETER ::   NR   = 377    ! Nb de reactions dans la chimie
+      INTEGER,PARAMETER ::   NLRT = 650    ! Pour l'UV 650 niveaux de 2 km 
+      
+!     Dummy parameters without microphysics
+!      + Just here to keep the whole stuff running without modif C sources
+!     ---------------------------------------------------------------------
+
+      INTEGER  :: utilaer(16)
+      INTEGER  :: aerprod = 0
+      INTEGER  :: htoh2   = 0
+
+!    Arguments
+!    ---------
+
+      INTEGER      nlon                   ! nb of horiz points
+      INTEGER      ny                     ! nb de composes (nqmax-nmicro) 
+      REAL*8         qy_c(nlon,klev,NC)     ! Especes chimiques apres adv.+diss.
+      character*10 nomqy_c(NC+1)          ! Noms des especes chimiques
+      REAL*8         declin_rad,ls_rad      ! declinaison et long solaire en radians
+      REAL*8         dtchim                 ! pas de temps chimie
+      REAL*8         ctemp(nlon,klev)       ! Temperature
+      REAL*8         cplay(nlon,klev)       ! pression (Pa)
+      REAL*8         cplev(nlon,klev+1)     ! pression intercouches (Pa)
+      REAL*8         czlay(nlon,klev)       ! altitude (m)
+      REAL*8         czlev(nlon,klev+1)     ! altitude intercouches (m)
+      
+      REAL*8         dqyc(nlon,klev,NC)     ! Tendances especes chimiques 
+      
+      INTEGER      NLEV                   ! nbp_lev+70 - a changer si =/ 55 layers ?
+      
+!    Local variables :
+!    -----------------
+
+      integer i,j,l,ic,jm1
+
+! variables envoyees dans la chimie: double precision
+
+      REAL*8  temp_c(NLEV)
+      REAL*8  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
+      REAL*8  cqy(NLEV,NC)    ! frac mol qui seront modifiees
+      REAL*8  cqy0(NLEV,NC)    ! frac mol avant modif
+      REAL*8  surfhaze(NLEV)
+      REAL*8  cprodaer(NLEV,4),cmaer(NLEV,4)
+      REAL*8  ccsn(NLEV,4),ccsh(NLEV,4)
+! rmil: milieu de couche, grille pour K,D,p,ct,T,y
+! rinter: intercouche (grille RA dans la chimie)
+      REAL*8  rmil(NLEV),rinter(NLEV),nb(NLEV)
+      REAL*8,save,allocatable :: kedd(:)
+
+      REAL*8  a,b
+      character str1*1,str2*2
+      REAL*8  latit
+      character*20 formt1,formt2,formt3
+      
+!     variables locales initialisees au premier appel
+
+      LOGICAL firstcall
+      DATA    firstcall/.true./
+      SAVE    firstcall
+      
+      integer dec,declinint,ialt
+      REAL*8  declin_c                       ! declinaison en degres
+      real*8  factalt,factdec,krpddec,krpddecp1,krpddecm1
+      real*8  duree
+      REAL*8,save :: mass(NC)
+      REAL*8,save,allocatable :: md(:,:)
+      REAL*8,save :: botCH4
+      DATA  botCH4/0.05/
+      real*8,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
+      REAL*8,save,allocatable :: krpd(:,:,:,:),krate(:,:)
+      integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
+      integer,save :: prod(200,NC),perte(2,200,NC)
+      character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
+      real*8  ficalt,oldq,newq,zalt
+      logical okfic
+
+
+!-----------------------------------------------------------------------
+!***********************************************************************
+!
+!    Initialisations :
+!    ----------------
+
+      duree = dtchim  ! passage en real*4 pour appel a routines C
+
+      IF (firstcall) THEN
+            
+        print*,'CHIMIE, premier appel'
+        
+! ************************************
+! Au premier appel, initialisation de toutes les variables 
+! pour les routines de la chimie.
+! ************************************
+
+	allocate(kedd(NLEV))
+
+        allocate(krpd(15,ND+1,NLRT,nbp_lat),krate(NLEV,NR),md(NLEV,NC))
+
+! Verification du nombre de composes: coherence common_mod et nqmax-nmicro
+! ----------------------------------
+
+        if (ny.ne.NC) then
+	   print*,'PROBLEME de coherence dans le nombre de composes:',ny,NC
+           STOP
+	endif
+
+! calcul de temp_c, densites et press_c en moyenne planetaire :
+! --------------------------------------------------------------
+
+        print*,'pression, densites et temp (init chimie):'
+        print*,'level, press_c, nb, temp_c'
+        DO l=1,klev
+         rinter(l)  = (zlevmoy(l)+rad)/1000.
+         rmil(l)    = (zlaymoy(l)+rad)/1000.
+!     temp_c (K):
+         temp_c(l)  = tmoy(l)
+!     press_c (mbar):
+         press_c(l) = playmoy(l)/100.
+!     nb (cm-3):
+         nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l))
+         print*,l,rmil(l)-rad/1000.,press_c(l),nb(l),temp_c(l)
+        ENDDO
+        rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
+
+! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+2,NLEV
+         rinter(l) = rinter(klev+1) &
+               + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
+         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
+       enddo
+       rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
+
+! lecture de tcp.ver, une seule fois
+! remplissage r1d,t1d,ct1d,p1d
+        open(11,file=TRIM(datadir)//'/tcp.ver',status='old')
+        read(11,*) dummy
+        do i=1,131
+          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
+          !print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
+        enddo 
+        close(11)
+
+! extension pour klev+1 a NLEV avec tcp.ver
+! la jonction klev/klev+1 est brutale... Tant pis ?
+        ialt=1
+        do l=klev+1,NLEV
+           zalt=rmil(l)-rad/1000.
+           do i=ialt,130
+            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
+              ialt=i
+            endif
+           enddo
+           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
+           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)  & 
+                           + log(p1d(ialt+1)) * factalt    )
+           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) &  
+                           + log(ct1d(ialt+1)) * factalt    )
+           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
+!          print*,l,zalt,press_c(l),nb(l),temp_c(l)
+        enddo
+        
+! caracteristiques des composes:        
+! -----------------------------
+        mass(:) = 0.0
+        call comp(nomqy_c,nb,temp_c,mass,md)
+        print*,'           Mass'
+        do ic=1,NC
+          print*,nomqy_c(ic),mass(ic)
+!         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
+        enddo	
+	        
+! taux de photodissociations:
+! --------------------------
+        call disso(krpd,nbp_lat) 
+
+! reactions chimiques:
+! -------------------
+        call chimie(nomqy_c,nb,temp_c,krate,reactif, &
+                    nom_perte,nom_prod,perte,prod)
+!        print*,'nom_prod, nom_perte:'
+!        do ic=1,NC
+!          print*,nom_prod(ic),nom_perte(ic)
+!        enddo
+!        print*,'premieres prod, perte(1:reaction,2:compagnon):'
+!        do ic=1,NC
+!          print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
+!        enddo
+
+!       l = klev-3
+!       print*,'krate a p=',press_c(l),' reactifs et produits:'
+!       do ic=1,ND+1
+!         print*,ic,krpd(7,ic,l,4)*nb(l),"  ",
+!    .     nomqy_c(reactif(1,ic)+1),
+!    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
+!    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
+!       enddo
+!       do ic=ND+2,NR
+!         print*,ic,krate(l,ic),"  ",
+!    .     nomqy_c(reactif(1,ic)+1),
+!    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
+!    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
+!       enddo
+
+
+!   init kedd
+!   ---------
+!   kedd stays constant with time and space 
+!   => linked to pressure rather than altitude...
+ 
+      kedd(:) = 5e2 ! valeur mise par defaut  
+               ! UNITE ? doit etre ok pour gptitan
+      do l=1,NLEV
+       zalt=rmil(l)-rad/1000.  ! z en km
+       if     ((zalt.ge.250.).and.(zalt.le.400.)) then 
+         kedd(l) = 10.**(3.+(zalt-250.)/50.)
+! 1E3 at 250 km
+! 1E6 at 400 km
+       elseif ((zalt.gt.400.).and.(zalt.le.600.)) then 
+         kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.)
+! 2E7 at 600 km
+       elseif ((zalt.gt.600.).and.(zalt.le.900.)) then 
+         kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.)
+! 1E8 above 900 km
+       elseif ( zalt.gt.900.                    ) then 
+        kedd(l) = 1.e8
+       endif
+      enddo
+
+      ENDIF  ! premier appel
+
+!***********************************************************************
+!-----------------------------------------------------------------------
+
+!   calcul declin_c (en degres)
+!   ---------------------------
+
+      declin_c = declin_rad*180./pi
+!      print*,'declinaison en degre=',declin_c
+       
+!***********************************************************************
+!***********************************************************************
+!
+!                BOUCLE SUR LES LATITUDES
+!
+      DO j=1,nlon 
+      
+      if (j.eq.1) then 
+         jm1=1
+      else
+         jm1=j-1
+      endif
+
+      if((j.eq.1).or.(klat(j).ne.klat(jm1))) then
+
+!***********************************************************************
+!***********************************************************************
+
+!-----------------------------------------------------------------------
+!
+!   Distance radiale (intercouches, en km)
+!   ----------------------------------------
+
+       do l=1,klev
+         rinter(l) = (rad+czlev(j,l))/1000.
+         rmil(l)   = (rad+czlay(j,l))/1000.
+!        print*,rinter(l)
+       enddo
+       rinter(klev+1)=(rad+czlev(j,klev+1))/1000.
+
+! au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+2,NLEV
+         rinter(l) = rinter(klev+1)  &
+               + (l-klev-1)*(3865.-rinter(klev+1))/(NLEV-klev-1)
+         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
+       enddo
+       rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
+
+!-----------------------------------------------------------------------
+!
+!   Temperature, pression (mbar), densite (cm-3)
+!   -------------------------------------------
+
+       DO l=1,klev
+!     temp_c (K):
+               temp_c(l)  = ctemp(j,l)
+!     press_c (mbar):
+               press_c(l) = cplay(j,l)/100.
+!     nb (cm-3):
+               nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l))
+       ENDDO
+! extension pour klev+1 a NLEV avec tcp.ver
+       ialt=1
+       do l=klev+1,NLEV
+           zalt=rmil(l)-rad/1000.
+           do i=ialt,130
+            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
+              ialt=i
+            endif
+           enddo
+           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
+           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt) &
+                           + log(p1d(ialt+1)) * factalt    )
+           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) &
+                           + log(ct1d(ialt+1)) * factalt    )
+           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
+       enddo
+               
+!-----------------------------------------------------------------------
+!
+!   passage krpd => krate 
+!   ---------------------
+! initialisation krate pour dissociations 
+
+      if ((declin_c*10+267).lt.14.) then
+          declinint = 0
+          dec       = 0
+      else 
+       if ((declin_c*10+267).gt.520.) then
+          declinint = 14
+          dec       = 534
+       else 
+          declinint = 1
+          dec       = 27
+          do while( (declin_c*10+267).ge.real(dec+20) )
+            dec       = dec+40
+            declinint = declinint+1
+          enddo
+       endif
+      endif
+      if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
+          factdec = ( declin_c - (dec-267)/10. ) / 4.
+      else
+          factdec = ( declin_c - (dec-267)/10. ) / 2.7
+      endif
+
+      do l=1,NLEV
+
+! INTERPOL EN ALT POUR k (krpd tous les deux km)
+        ialt    = int((rmil(l)-rad/1000.)/2.)+1
+        factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
+
+        do i=1,ND+1
+          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt) &
+                   + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
+          krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt) &
+                   + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
+          krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt) &
+                   + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
+
+                    ! ND+1 correspond a la dissociation de N2 par les GCR
+          if ( factdec.lt.0. ) then 
+        krate(l,i) = krpddecm1 * abs(factdec)  &
+                  + krpddec   * ( 1 + factdec)
+          endif
+          if ( factdec.gt.0. ) then
+        krate(l,i) = krpddecp1 * factdec       &
+                  + krpddec   * ( 1 - factdec)
+          endif
+          if ( factdec.eq.0. ) krate(l,i) = krpddec
+        enddo        
+      enddo        
+
+!-----------------------------------------------------------------------
+!
+!   composition 
+!   ------------
+
+       do ic=1,NC
+        do l=1,klev
+          cqy(l,ic)  = qy_c(j,l,ic) 
+          cqy0(l,ic) = cqy(l,ic)
+        enddo
+       enddo
+
+! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
+
+      WRITE(str2,'(i2.2)') klat(j)
+      fich = "comp_"//str2//".dat"
+      inquire (file=fich,exist=okfic)
+      if (okfic) then
+       open(11,file=fich,status='old')
+! premiere ligne=declin
+       read(11,'(A15)') ficdec
+       write(curdec,'(E15.5)') declin_c
+! si la declin est la meme: ce fichier a deja ete reecrit 
+! on lit la colonne t-1 au lieu de la colonne t
+! (cas d une bande de latitude partagee par 2 procs)
+       do ic=1,NC
+        read(11,'(A10)') name
+        if (name.ne.nomqy_c(ic)) then
+          print*,"probleme lecture ",fich
+          print*,name,nomqy_c(ic)
+          stop
+        endif
+        if (ficdec.eq.curdec) then
+          do l=klev+1,NLEV
+            read(11,*) ficalt,cqy(l,ic),newq
+          enddo
+        else
+          do l=klev+1,NLEV
+            read(11,*) ficalt,oldq,cqy(l,ic)
+          enddo
+        endif
+       enddo
+       close(11)
+      else       ! le fichier n'est pas la
+       do ic=1,NC
+        do l=klev+1,NLEV
+          cqy(l,ic)=cqy(klev,ic)
+        enddo
+       enddo
+      endif
+      cqy0 = cqy
+
+!-----------------------------------------------------------------------
+!
+!   Appel de chimietitan
+!   --------------------
+       
+       call gptitan(rinter,temp_c,nb,                          &
+                   nomqy_c,cqy,                               &
+                   duree,(klat(j)-1),mass,md,                 &
+                   kedd,botCH4,krate,reactif,                 &
+                   nom_prod,nom_perte,prod,perte,             &
+                   aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,  &
+                   htoh2,surfhaze)
+       
+!   Tendances composition 
+!   ---------------------
+     
+       do ic=1,NC
+        do l=1,klev
+          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
+        enddo
+       enddo
+
+
+!-----------------------------------------------------------------------
+!
+!   sauvegarde compo sur NLEV 
+!   -------------------------
+
+! dans fichier compo_klat(j) (01 à 49)
+       
+      open(11,file=fich,status='replace')
+! premiere ligne=declin
+      write(11,'(E15.5)') declin_c
+      do ic=1,NC
+       write(11,'(A10)') nomqy_c(ic)
+       do l=klev+1,NLEV
+        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-rad/1000.,cqy0(l,ic),cqy(l,ic)
+       enddo
+      enddo
+      close(11)
+
+!***********************************************************************
+!***********************************************************************
+!              FIN: BOUCLE SUR LES LATITUDES 
+
+      else      ! same latitude, we don't do calculations again
+        dqyc(j,:,:) = dqyc(jm1,:,:)
+      endif
+
+      ENDDO
+     
+!***********************************************************************
+!***********************************************************************
+
+      firstcall = .false.
+      
+      end subroutine calchim
Index: /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1672)
@@ -4,6 +4,6 @@
       logical,save :: callrad,corrk,calldifv,UseTurbDiff
 !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
-      logical,save :: calladj,callsoil
-!$OMP THREADPRIVATE(calladj,callsoil)
+      logical,save :: calladj,callsoil,callchim
+!$OMP THREADPRIVATE(calladj,callsoil,callchim)
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
@@ -38,10 +38,11 @@
       logical,save :: oblate
 !$OMP THREADPRIVATE(nosurf,oblate)
-
+      
+      integer,save :: ichim
       integer,save :: iddist
       integer,save :: iaervar
       integer,save :: iradia
       integer,save :: startype
-!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype)
+!$OMP THREADPRIVATE(ichim,iddist,iaervar,iradia,startype)
 
       real,save :: Fat1AU
Index: /trunk/LMDZ.TITAN/libf/phytitan/comcstfi_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/comcstfi_mod.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/comcstfi_mod.F90	(revision 1672)
@@ -11,5 +11,6 @@
       REAL,SAVE :: omeg ! planet rotation rate (rad/s)
       REAL,SAVE :: avocado ! something like 6.022e23
-!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado)
+      REAL,SAVE :: kbol ! something like 1.380649e-23
+!$OMP THREADPRIVATE(pi,rad,g,r,cpp,rcp,mugaz,omeg,avocado,kbol)
 
 END MODULE comcstfi_mod
Index: /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90	(revision 1672)
@@ -59,9 +59,4 @@
    press_PL(:)=press_PL(:)*1E-2
    
-   ! Correction of the table according to Lavvas et al. 2010 (ext coeff = constant < 80 km / 20 mbar ) because of condensation
-!   do i=1,49
-!     ext_PL(i,:) = ext_PL(50,:)
-!   enddo
-
    firstcall=.false.
 endif
@@ -101,4 +96,5 @@
 ! Lin-Log interpolation : linear on wln, logarithmic on press
 
+if ((iw.ne.nbwl_PL).and.(ip.ne.nblev_PL)) then
 taeros = ( ext_PL(ip,iw)*(1.-factw)   + ext_PL(ip,iw+1)  *factw ) ** (1.-factp) &
         *( ext_PL(ip+1,iw)*(1.-factw) + ext_PL(ip+1,iw+1)*factw ) ** factp
@@ -109,4 +105,10 @@
 cbar   = ( asf_PL(ip,iw)*(1.-factw)   + asf_PL(ip,iw+1)  *factw ) ** (1.-factp) &
         *( asf_PL(ip+1,iw)*(1.-factw) + asf_PL(ip+1,iw+1)*factw ) ** factp
+else if ((iw.eq.nbwl_PL).and.(ip.ne.nblev_PL)) then
+taeros =  ext_PL(ip,iw) ** (1.-factp) * ext_PL(ip+1,iw) ** factp
+
+ssa    =  ssa_PL(ip,iw) ** (1.-factp) * ssa_PL(ip+1,iw) ** factp
+
+cbar   =  asf_PL(ip,iw) ** (1.-factp) * asf_PL(ip+1,iw) ** factp
 
 ! In case of vertical extension over the max of the table
@@ -114,5 +116,5 @@
 ! Arbitray threshold pressure value, just to deal with the last level press=0
 
-if(ip.eq.nblev_PL) then 
+else if(ip.eq.nblev_PL) then 
 
   tmp_p = press
@@ -122,10 +124,16 @@
   endif
 
-  fact_t = log10( ( ext_PL(ip,iw)*(1.-factw) + ext_PL(ip,iw+1)  *factw )   & 
-           / ( ext_PL(ip-5,iw)*(1.-factw) + ext_PL(ip-5,iw+1)  *factw ) )
-  fact_s = log10( ( ssa_PL(ip,iw)*(1.-factw) + ssa_PL(ip,iw+1)  *factw )   & 
-           / ( ssa_PL(ip-5,iw)*(1.-factw) + ssa_PL(ip-5,iw+1)  *factw ) )
-  fact_c = log10( ( asf_PL(ip,iw)*(1.-factw) + asf_PL(ip,iw+1)  *factw )   & 
-           / ( asf_PL(ip-5,iw)*(1.-factw) + asf_PL(ip-5,iw+1)  *factw ) )
+  if (iw.ne.nbwl_PL) then
+    fact_t = log10( ( ext_PL(ip,iw)*(1.-factw) + ext_PL(ip,iw+1)  *factw )   & 
+             / ( ext_PL(ip-5,iw)*(1.-factw) + ext_PL(ip-5,iw+1)  *factw ) )
+    fact_s = log10( ( ssa_PL(ip,iw)*(1.-factw) + ssa_PL(ip,iw+1)  *factw )   & 
+             / ( ssa_PL(ip-5,iw)*(1.-factw) + ssa_PL(ip-5,iw+1)  *factw ) )
+    fact_c = log10( ( asf_PL(ip,iw)*(1.-factw) + asf_PL(ip,iw+1)  *factw )   & 
+             / ( asf_PL(ip-5,iw)*(1.-factw) + asf_PL(ip-5,iw+1)  *factw ) )
+  else
+    fact_t = log10( ext_PL(ip,iw) / ext_PL(ip-5,iw) )
+    fact_s = log10( ssa_PL(ip,iw) / ssa_PL(ip-5,iw) )
+    fact_c = log10( asf_PL(ip,iw) / asf_PL(ip-5,iw) )
+  endif            
 
   fact_t = fact_t / log10( press_PL(ip) / press_PL(ip-5) )
Index: /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90	(revision 1672)
+++ /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90	(revision 1672)
@@ -0,0 +1,164 @@
+      SUBROUTINE inicondens(ny,press,temp,nomy,yc)
+
+!==================================================================================
+!        Purpose
+!        -------
+!        Initialisation des profils de saturation des traceurs
+!
+!        Authors
+!        -------
+!        S. Lebonnois
+!        J. Vatant d'Ollone (2017) : adaptation au nouveau titan et passage en F90
+!==================================================================================
+
+        
+!-----------------------------------------------------------------------
+!   Declarations:
+!   -------------
+
+      use dimphy
+      use mod_grid_phy_lmdz, only: nbp_lev
+      IMPLICIT NONE
+
+!    Arguments :
+!    -----------
+      INTEGER, INTENT(IN)     :: ny
+      CHARACTER*10,INTENT(IN) :: nomy(ny+1)
+      REAL, INTENT(IN)        :: press(nbp_lev),temp(nbp_lev)  ! pressure in mbar !
+      REAL, INTENT(OUT)       :: yc(nbp_lev,ny)
+      
+!    Local variables :
+!    -----------------
+      INTEGER                 ::  l,ic
+      REAL                    ::  sy,x
+       
+      do ic=1,ny
+       print*, 'traceur CH(', ic, ')=', nomy(ic),'------------'
+           do l=1,nbp_lev
+
+!  Par defaut, yc est a 1 c'est a dire qu'on ne condense pas
+              yc(l,ic)=1.
+
+              if(nomy(ic).eq."CH4") then
+                 if (temp(l).lt.90.65) then
+                   yc(l,ic)= 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / TEMP(l) - &
+                        115352.19e0 ) / TEMP(l) + 4055.6016e0 ) / TEMP(l)     &
+                        + 453.92414e0 ) / TEMP(l) ) / PRESS(l) * 1013.25e0
+                 else
+                    yc(l,ic)= 10.0**(3.901408e0 - ( ( 154567.02e0 / TEMP(l) -  &
+                         1598.8512e0 ) / TEMP(l) + 437.54809e0 ) / TEMP(l))    &
+                         / PRESS(l) * 1013.25e0;
+                 endif
+! maintient a 1.4% minimum               
+                 if (yc(l,ic).lt.0.014) yc(l,ic)=0.014
+              endif
+
+              if(nomy(ic).eq."C2H2") then
+                 yc(l,ic)= 10.0**(6.09748e0-1644.1e0/TEMP(l)+7.42346e0     &
+                      * alog10(1.0e3/TEMP(l)) ) / PRESS(l)*1013.25e0/760.0
+              endif
+
+              if(nomy(ic).eq."C2H4") then
+                 if (temp(l).lt.89.0) then
+                   yc(l,ic)= 10.0**(1.5477e0 + (1.0e0/TEMP(l) - 0.011e0)     &
+                        *(16537.0e0*(1.0e0/TEMP(l) - 0.011e0) - 1038.1e0))   &
+                        / PRESS(l) * 1.01325e0 / 760.0
+                 elseif (temp(l).lt.104.0) then
+                   yc(l,ic)= 10.0**(8.724e0 - 901.6e0/(TEMP(l) - 2.555e0) )  &
+                        / PRESS(l) * 1013.25e0 / 760.0
+                 elseif (temp(l).lt.120.0) then
+                   yc(l,ic)= 10.0**(50.79e0 - 1703.0e0/TEMP(l) - 17.141e0 *  &
+                        alog10(TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0
+                 elseif (temp(l).lt.155.0) then
+                   yc(l,ic)= 10.0**(6.74756e0 - 585.0e0/(TEMP(l) - 18.16e0) ) &
+                        / PRESS(l) * 1013.25e0 / 760.0
+                 endif
+              endif
+
+              if(nomy(ic).eq."C2H6") then
+                 if (temp(l).lt.90.) then
+                   yc(l,ic)= 10.0**(10.01e0-1085.0e0/(TEMP(l)-0.561e0) )  &
+                        / PRESS(l) * 1013.25e0 / 760.0e0
+                 else
+                   yc(l,ic)= 10.0**(5.9366e0 - 1086.17e0/TEMP(l) + 3.83464e0 *  &
+                        alog10(1.0e3/TEMP(l)) ) / PRESS(l)*1013.25e0/760.0
+                 endif
+              endif
+
+              if((nomy(ic).eq."CH3CCH").or.(nomy(ic).eq."CH2CCH2")) then
+                 yc(l,ic)= 10.0**(2.8808e0 - 4.5e0*(249.9e0 - TEMP(l))   &
+                      /(1.15e0*TEMP(l) - 37.485e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C3H6")  then
+                 yc(l,ic)= 10.0**(7.4463e0 - 1028.5654e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C3H8")  then
+                 yc(l,ic)= 10.0**(7.217e0 - 994.30251e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if((nomy(ic).eq."C4H2").or.(nomy(ic).eq."C4H2s")) then
+                 yc(l,ic)= 10.0**(96.26781e0 - 4651.872e0/TEMP(l) - 31.68595e0 &
+                      *alog10(TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C4H4")  then
+                 yc(l,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(TEMP(l)-43.15e0) ) / PRESS(l)
+              endif
+
+              if(nomy(ic).eq."C4H6")  then
+                 yc(l,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - TEMP(l)) &
+                      /(1.15e0*TEMP(l) - 39.345e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C4H10")  then
+                 yc(l,ic)= 10.0**(8.446e0 - 1461.2e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C6H2")  then
+                 yc(l,ic)= 10.0**(4.666e0 - 4956e0/TEMP(l) + 25.845e0 * &
+                      alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C8H2")  then
+                 yc(l,ic)= 10.0**(3.95e0 - 6613e0/TEMP(l) + 35.055e0 * &
+                      alog10(1.0e3/TEMP(l)) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."AC6H6")  then
+                 x = 1.0e0 - TEMP(l) / 562.2e0
+                 yc(l,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x &
+                      - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/TEMP(l) ) / PRESS(l)
+              endif
+
+              if(nomy(ic).eq."HCN")  then
+                 yc(l,ic)= 10.0**(8.6165e0 - 1516.5e0/(TEMP(l) - 26.2e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."CH3CN")  then
+                 yc(l,ic)= 10.0**(8.458e0 - 1911.7e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C2H3CN")  then
+                 yc(l,ic)= 10.0**(9.3051e0 - 2782.21/(TEMP(l) - 51.15e0) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."NCCN")  then
+                 yc(l,ic)=  10.0**(7.454e0 - 1832e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."HC3N")  then
+                 yc(l,ic)= 10.0**(7.7446e0 - 1453.5609e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+              if(nomy(ic).eq."C4N2")  then
+                 yc(l,ic)= 10.0**(8.269e0 - 2155.0e0/TEMP(l) ) / PRESS(l) * 1013.25e0 / 760.0e0
+              endif
+
+           enddo
+       enddo
+
+       print*, 'inicondens end'
+       
+      end subroutine inicondens
Index: /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1672)
@@ -18,5 +18,5 @@
                               init_time, daysec, dtphys
   use comcstfi_mod, only: rad, cpp, g, r, rcp, &
-                          mugaz, pi, avocado
+                          mugaz, pi, avocado, kbol
   use planete_mod, only: nres
   use planetwide_mod, only: planetwide_sumval
@@ -84,4 +84,5 @@
   pi=2.*asin(1.)
   avocado = 6.02214179e23   ! added by RW
+  kbol = 1.38064852e-23  ! added by JVO for Titan chem
 
   ! Initialize some "temporal and calendar" related variables
@@ -162,6 +163,5 @@
      flatten = 0.0
      call getin_p("flatten", flatten)
-     write(*,*) "flatten = ", flatten
-         
+     write(*,*) "flatten = ", flatten         
 
      write(*,*) "Needed if oblate=.true.: J2"
@@ -300,4 +300,24 @@
      call getin_p("graybody",graybody)
      write(*,*)" graybody = ",graybody
+
+! Chemistry
+
+     write(*,*) "Run with or without chemistry?"
+     callchim=.false. ! default value
+     call getin_p("callchim",callchim)
+     write(*,*) " callchim = ",callchim
+
+     ! sanity check
+     if (callchim.and.(.not.tracer)) then
+       print*,"You are running chemistry without tracer"
+       print*,"Please start again with tracer =.true."
+       stop
+     endif
+
+     write(*,*)"Chemistry is computed every ichim", &
+                   " physical timestep"
+     ichim=1 ! default value
+     call getin_p("ichim",ichim)
+     write(*,*)" ichim = ",ichim
 
 ! Soil model
@@ -497,5 +517,8 @@
   PRINT*,'        or each ',iradia*dtphys,' seconds'
   PRINT*
-
+  PRINT*,'inifis: Chemistry is computed:'
+  PRINT*,'           each ',ichim,' physical time-step'
+  PRINT*,'        or each ',ichim*dtphys,' seconds'
+  PRINT*
 
 !-----------------------------------------------------------------------
Index: /trunk/LMDZ.TITAN/libf/phytitan/initracer.F
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/initracer.F	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/initracer.F	(revision 1672)
@@ -19,4 +19,5 @@
 c   ------
 c            Ehouarn Millour (oct. 2008) identify tracers by their names
+c            Jan Vatant d'Ollone (fev. 2017) chemical species for Titan
 c=======================================================================
 
@@ -26,5 +27,4 @@
       character(len=20) :: txt ! to store some text
       integer iq,ig,count
-      real r0_lift , reff_lift
 !      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
 
@@ -65,16 +65,284 @@
         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_n2=0
-
+      
+      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.    
-
+      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:
Index: /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1672)
@@ -33,4 +33,7 @@
       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
       use time_phylmdz_mod, only: daysec
+      use logic_mod, only: moyzon_ch
+      use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
+                            zplaybar, zzlevbar, zzlaybar, ztfibar
       use callkeys_mod
       use vertical_layers_mod, only: presnivs, pseudoalt
@@ -70,7 +73,8 @@
 !
 !      V. Tracers
-!         V.1. Aerosols and particles.
-!         V.2. Updates (pressure variations, surface budget).
-!         V.3. Surface Tracer Update.
+!         V.1. Chemistry
+!         V.2. Aerosols and particles.
+!         V.3. Updates (pressure variations, surface budget).
+!         V.4. Surface Tracer Update.
 !
 !      VI. Surface and sub-surface soil temperature.
@@ -140,5 +144,7 @@
 !           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
 !           Purge of the code : M. Turbet (2015)
-!           Fork for Titan and clean of all too-generic (ocean, water, co2 ...) routines : J. Vatant d'Ollone (2017)
+!           Fork for Titan : J. Vatant d'Ollone (2017)
+!                + clean of all too-generic (ocean, water, co2 ...) routines
+!                + Titan's chemistry
 !============================================================================================
 
@@ -277,4 +283,7 @@
       real zdqsed(ngrid,nlayer,nq)    ! Callsedim routine.
       real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
+      
+      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
+      real zdqmph(ngrid,nlayer,nq)    ! Microphysical tendency ( condensation routine only for now, no microphysical routines ).
                         
       ! For Winds : (m/s/s)
@@ -357,5 +366,33 @@
       real reffcol(ngrid,naerkind)
 
-
+! Local variables for Titan chemistry and microphysics (JVO 2017)
+! ----------------------------------------------------------------------------
+     
+      integer,parameter  :: nmicro=0 ! Temporary ! To be put in start/def
+
+      real ctimestep ! Chemistry timestep (s)
+ 
+      ! Grandeurs en moyennes zonales ------------------------
+      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)   
+
+      ! 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)
+
+!-----------------------------------------------------------------------------
+      
 !==================================================================================================
 
@@ -401,5 +438,8 @@
         ALLOCATE(zdtsw(ngrid,nlayer))
         ALLOCATE(tau_col(ngrid))
-
+        ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
+        ALLOCATE(qysat(nlayer,nq-nmicro))
+        ALLOCATE(nomqy(nq-nmicro+1))
+        
         ! This is defined in comsaison_h
         ALLOCATE(mu0(ngrid))
@@ -433,7 +473,6 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call phyetat0(startphy_file,                                 &
-                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
-                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)
+         call phyetat0(startphy_file,ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
+                       day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf)                       
          if (.not.startphy_file) then
            ! additionnal "academic" initialization of physics
@@ -514,4 +553,31 @@
                           ptimestep,pday+nday,time_phys,cell_area,          &
                           albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
+         endif
+
+         ! Sanity check for chemistry
+         if ( ((moyzon_ch).and.(.not.callchim)) .or. ((.not.moyzon_ch).and.(callchim)) ) then
+            print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
+            print *, "This is not compatible..."
+            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
          
@@ -601,4 +667,42 @@
       enddo     
 
+      ! -------------------------------Taken from old Titan --------------------------
+      ! zonal averages needed
+      if (moyzon_ch) then
+         
+         zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
+         ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
+         !       zzlaybar(1,:)=RG*RA*RA/(RG*RA-(zphibar(1,:)+zphisbar(1)))-RA
+         zzlevbar(1,1)=zphisbar(1)/g
+         DO l=2,nlayer
+            z1=(zplaybar(1,l-1)+zplevbar(1,l))/(zplevbar(1,l-1)-zplevbar(1,l))
+            z2=(zplevbar(1,l)  +zplaybar(1,l))/(zplevbar(1,l)  -zplaybar(1,l))
+            zzlevbar(1,l)=(z1*zzlaybar(1,l-1)+z2*zzlaybar(1,l))/(z1+z2)
+         ENDDO
+         zzlevbar(1,nlayer+1)=zzlaybar(1,nlayer)+(zzlaybar(1,nlayer)-zzlevbar(1,nlayer))
+
+         DO ig=2,ngrid
+            if (latitude(ig).ne.latitude(ig-1)) then
+               DO l=1,nlayer
+                  zzlaybar(ig,l)=(zphibar(ig,l)+zphisbar(ig))/g
+                  ! SI ON TIENT COMPTE DE LA VARIATION DE G AVEC L'ALTITUDE:
+                  !zzlaybar(ig,l)=RG*RA*RA/(RG*RA-(zphibar(ig,l)+zphisbar(ig)))-RA
+               ENDDO
+               zzlevbar(ig,1)=zphisbar(ig)/g
+               DO l=2,nlayer
+                  z1=(zplaybar(ig,l-1)+zplevbar(ig,l))/ (zplevbar(ig,l-1)-zplevbar(ig,l))
+                  z2=(zplevbar(ig,l)  +zplaybar(ig,l))/(zplevbar(ig,l)  -zplaybar(ig,l))
+                  zzlevbar(ig,l)=(z1*zzlaybar(ig,l-1)+z2*zzlaybar(ig,l))/(z1+z2)
+               ENDDO
+               zzlevbar(ig,nlayer+1)=zzlaybar(ig,nlayer)+(zzlaybar(ig,nlayer)-zzlevbar(ig,nlayer))
+            else
+               zzlaybar(ig,:)=zzlaybar(ig-1,:)
+               zzlevbar(ig,:)=zzlevbar(ig-1,:)
+            endif
+         ENDDO
+
+      endif  ! moyzon
+      ! -------------------------------------------------------------------------------------
+      
       ! Compute potential temperature
       ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
@@ -896,7 +1000,59 @@
       if (tracer) then
       
-      
   ! -------------------------
-  !   V.1. Aerosol particles
+  !   V.1. Chemistry
+  ! -------------------------
+         if (callchim) then
+
+            zplev(:,:) = zplevbar(:,:)
+            zplay(:,:) = zplaybar(:,:)
+            zzlev(:,:) = zzlevbar(:,:)
+            zzlay(:,:) = zzlaybar(:,:)
+            ztemp(:,:) = ztfibar(:,:)
+
+            if (nq.gt.nmicro) then
+               do iq = nmicro+1,nq
+                  ychim(:,:,iq-nmicro) = pq(:,:,iq) 
+               enddo
+            endif
+
+            ! Condensation tendency after the transport
+            do iq=1,nq-nmicro
+               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
+                     endif
+                  enddo
+               enddo
+            enddo
+
+            if( mod(icount-1,ichim).eq.0.) then
+               
+               print *, "On passe dans la chimie..."
+               
+               call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
+                           ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
+
+            ! JVO 2017 : NLEV = nlayer+70, en accord avec le C. Quid si nlay=/ 55 ?
+               
+            endif
+            
+            if (nq.gt.nmicro) then
+               ! We convert tendencies back to 3D and mass mixing ratio
+               do iq=nmicro+1,nq
+                  zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro)*pq(:,:,iq) / ychim(:,:,iq-nmicro)
+                  zdqmph(:,:,iq) = dycmph(:,:,iq-nmicro)*pq(:,:,iq) / ychim(:,:,iq-nmicro)
+               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)
+               
+            endif
+            
+         endif
+      
+  ! -------------------------
+  !   V.2. Aerosol particles
   ! -------------------------
 
@@ -919,5 +1075,5 @@
 
   ! ---------------
-  !   V.2. Updates
+  !   V.3 Updates
   ! ---------------
 
@@ -951,5 +1107,5 @@
 
   ! -----------------------------
-  !   V.3. Surface Tracer Update
+  !   V.4. Surface Tracer Update
   ! -----------------------------
 
Index: /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1671)
+++ /trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90	(revision 1672)
@@ -21,7 +21,60 @@
 ! 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
-!$OMP THREADPRIVATE(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
