Index: /trunk/LMDZ.TITAN/README
===================================================================
--- /trunk/LMDZ.TITAN/README	(revision 1907)
+++ /trunk/LMDZ.TITAN/README	(revision 1908)
@@ -1399,4 +1399,5 @@
 + Upper chemistry fields ( and chemistry tracers if needed ) in newstart are
  regrided in vert_regrid_kim.F90
++ Interface with chemistry routines itself will be done in a later commit (r1908)
 
 == 13/01/18 == JVO
@@ -1423,4 +1424,16 @@
 == 25/02/18 == JVO
 Fix two bugs, one related to zonal means and one to rat_mmol
-+ Anyway zonal means apparently lead to crashes in microphysics and chemistry.
-+ Chemistry is still on the bench, but for microphysics, please use full 3D !!
++ Zonal means may lead to crashes in microphysics, please use full 3D, which is very quick !!
+
+== 08/03/18 == JVO
+Making chemistry handling more flexible - Major and Final Step (hopefully) !
+After preliminary commits r1871-86-87-91-94-98, here comes the major update of the interface
+with photochemical module + fix how tendencies for chem and mufi tracers are managed in physiq_mod !
++ Major modifs in : 
+++ calchim.F90 to comply with flexible resolution, parallelism, upper pressure grid ...
+++ physiq_mod.F90 where there was a bug on the update of the tracers and their tendencies for calchim
+   and calmufi ( we actually were sending non-updated fields to these processes )
+   We also now put the same tendency on all longitudes within a lat band and not
+   relative tendency if 2D chemistry ( and we set to zero if ever negs are created )
++ Also modifs to have chemistry in 1D in rcm1d ( and moved gr_kim_vervack in phytitan to be accessible for 1d )
++ In chemistry added a check.c to verify coherence of sizes between C and Fortran
Index: /trunk/LMDZ.TITAN/libf/chimtitan/check.c
===================================================================
--- /trunk/LMDZ.TITAN/libf/chimtitan/check.c	(revision 1908)
+++ /trunk/LMDZ.TITAN/libf/chimtitan/check.c	(revision 1908)
@@ -0,0 +1,51 @@
+/* Sanity check routine about coherence between C and Fortran parameters for chemistry */
+/* Author : Jan Vatant d'Ollone - 2018 */
+
+#include "titan.h"
+
+void check_(
+     int *NLEV_fort, int *NLD_fort,
+     int *NLRT_fort, int *NC_fort   )
+{
+
+   if( (*NLEV_fort) != NLEV )
+   {
+    printf("\n");
+    printf("!! Coherence problem between C and Fortran !!\n");
+    printf("In C, NLEV=%d whereas in Fortran nlaykim_tot=%d\n",NLEV,*(NLEV_fort));
+    printf("You surely didn't modified titan.h according to your startfile ! \n");
+    printf("I abort ...\n");
+    exit(0);
+   }
+
+   if( (*NLD_fort) != NLD )
+   {
+    printf("\n");
+    printf("!! Coherence problem between C and Fortran !!\n");
+    printf("The number of levels with diffusion must be Numbers-of-levels-of-GCM minus 15 !\n");
+    printf("In C, NLD=%d whereas in Fortran klev-15=%d\n",NLD,*(NLD_fort));
+    printf("You surely didn't modified titan.h according to your startfile ! \n");
+    printf("I abort ...\n");
+    exit(0);
+   }
+
+   if( (*NLRT_fort) != NLRT )
+   {
+    printf("\n");
+    printf("!! Coherence problem between C and Fortran for the number of UV levels !!\n");
+    printf("In C in titan.h, NLRT=%d whereas in Fortran comchem_h, nlrt_kim=%d\n",NLRT,*(NLRT_fort));
+    printf("I abort ...\n");
+    exit(0);
+   }
+
+   if( (*NC_fort) != NC )
+   {
+    printf("\n");
+    printf("!! Coherence problem between C and Fortran for number of compounds !!\n");
+    printf("In C, NC=%d whereas in Fortran nkim=%d\n",NC,*(NC_fort));
+    printf("You have maybe changed the chemistry but running with old startfiles ! \n");
+    printf("I abort ...\n");
+    exit(0);
+   }
+
+}
Index: /trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c
===================================================================
--- /trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c	(revision 1908)
@@ -1,7 +1,4 @@
 /* gptitan: photochimie */
 /* GCCM */
-
-/* tout est passe en simple precision */
-/* sauf pour l'inversion de la matrice */
 
 /* nitriles et hydrocarbures separes pour l'inversion */
@@ -14,5 +11,5 @@
      double *RA, double *TEMP, double *NB, 
      char CORPS[][10], double Y[][NLEV],
-     double *FIN, int *LAT, double *MASS, double MD[][NLEV],
+     double *FIN, double *LAT, double *MASS, double MD[][NLEV],
      double *KEDD, double *botCH4, double KRATE[][NLEV],
      int reactif[][5], int *nom_prod, int *nom_perte, 
@@ -53,5 +50,5 @@
 
 /* DEBUG */
-      printf("CHIMIE: lat=%d\n",(*LAT)+1);
+      printf("CHIMIE: lat=%g\n",(*LAT));
 /**/
 
@@ -65,5 +62,5 @@
    strcat( outlog, ".log" );
    out = fopen( outlog, "w" );
-   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
+   fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
    fclose( out );
 
@@ -79,5 +76,5 @@
 /*  DEBUG 
             out = fopen( outlog, "a" );
-   fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
+   fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
             fclose( out );
 */
@@ -236,5 +233,5 @@
              }
 /* DEBUG
-printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
+printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
 if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
       printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
@@ -286,5 +283,5 @@
 
 /* DEBUG 
-printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
+printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
 if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
 printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
@@ -511,5 +508,5 @@
                   {
                      out = fopen( outlog, "a" );
-                     fprintf( out, "Lat no %d;", (*LAT)+1);
+                     fprintf( out, "Latitude %g;", (*LAT));
                      fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
                      fclose( out );
@@ -666,5 +663,5 @@
              }
 /* DEBUG
-printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
+printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
 if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
       printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
@@ -716,5 +713,5 @@
 
 /* DEBUG 
-printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
+printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
 if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
 printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
@@ -810,5 +807,5 @@
                   {
                      out = fopen( outlog, "a" );
-                     fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
+                     fprintf( out, "Latitude: %g; declin:%e;", (*LAT), (*DECLIN) );
                      fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
                      fclose( out );
Index: /trunk/LMDZ.TITAN/libf/chimtitan/titan.h
===================================================================
--- /trunk/LMDZ.TITAN/libf/chimtitan/titan.h	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/chimtitan/titan.h	(revision 1908)
@@ -7,6 +7,6 @@
 
 #define R0    (double)(2575.0) /* Titan's radius */
-#define NLEV  (int)(125)  /* Nbre de niv verticaux - =llm+70 dans common_mod */
-#define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff */
+#define NLEV  (int)(133)  /* Nbre de niv verticaux - =llm+70 dans common_mod -> Need to be coherent with the vertical grid used !! */
+#define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff -> Need to be coherent with the vertical grid used !! */
 #define NLRT  (int)(650)  /* Nbre de niv verticaux dans table fmoy - aussi dans common_mod */
 
Index: unk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/gr_kim_vervack.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/gr_kim_vervack.F90	(revision 1907)
+++ 	(revision )
@@ -1,151 +1,0 @@
-SUBROUTINE gr_kim_vervack
-
-  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  !
-  ! Purpose : * Calculates the number of layers needed for upper chemistry
-  ! -------   based on the GCM vertical grid size, depending on Ptop=ap.
-  !           * Generates also the pressure grid at mid-layer for upper levels.
-  !           * We use an upper atmosphere profile based on Voyager 1 data 
-  !           (Vervack et al, 2004) to have dz=10km between Ptop and 1300km.
-  !
-  ! Author : Jan Vatant d'Ollone (2017)
-  ! ------
-  !
-  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  USE comchem_h, ONLY: nlaykim_up, preskim
-  USE comvert_mod, ONLY: ap, aps
-  USE datafile_mod, ONLY: datadir
-
-  IMPLICIT NONE
-
-  INCLUDE "dimensions.h"
-
-  REAL*8  :: ptop, zptop_ver, zkim
-  REAL*8  :: factz
-  REAL*8  :: grkim_dz=10 ! Discretisation of chemistry grid (km)
-
-
-  REAL*8  :: rver(131),tver(131),ctver(131),pver(131) 
-
-  INTEGER :: i, j, itop, jlay, jold
-  
-  LOGICAL :: foundver
-
-  LOGICAL :: crop
-
-  ! ------------------------
-  ! 1. Read Vervack profile
-  ! ------------------------
-
-  ! Check the file exists
-  INQUIRE(file=TRIM(datadir)//'/tcp.ver',exist=foundver)
-  IF(.NOT.foundver) THEN
-    WRITE(*,*)'The file ',TRIM(datadir)//'/tcp.ver'
-    WRITE(*,*)'was not found by gr_kim_vervack.F90, exiting'
-    WRITE(*,*)'Check that your path to datagcm:',trim(datadir)
-    WRITE(*,*)'is correct. You can change it in callphys.def with:'
-    WRITE(*,*)'datadir = /absolute/path/to/datagcm'
-    WRITE(*,*)'Also check that file tcp.ver exists there.'
-  ENDIF
-
-  ! Read file
-  OPEN(11,file=TRIM(datadir)//'/tcp.ver',status='old')
-  READ(11,*)
-  DO i=1,131
-     READ(11,*) rver(i),tver(i),ctver(i),pver(i)
-     pver(i) = pver(i)*100.0 ! mbar->Pa
-  ENDDO
-  CLOSE(11)
-  
-  ! --------------------------------------------------------------
-  ! 2. Define ptop as the value aps should have if it wasn't zero
-  ! assuming ap(llm)-aps(llm) half pressure thickness of top-layer
-  ! --------------------------------------------------------------
-  
-  ! NB : At the top of the model we are in pure pressure coord. -> ap
-  
-  ptop = 2.0*aps(llm) - ap(llm)
-  
-  ! --------------------------------------------
-  ! 3. Interpolate Ptop and equivalent altitude 
-  ! --------------------------------------------
-  
-  itop=1
-  
-  DO i=2,131
-    IF ( ptop .LT. pver(i) ) THEN
-      itop=i
-    ENDIF
-  ENDDO
-  
-  ! Crazy case in a far far away future where GCM top will reach 1300km
-  IF ( itop .eq. 131 ) THEN
-    WRITE(*,*) " You've reach the bounds of Vervack profile ... "
-    WRITE(*,*) " Congrats but it wasn't predicted in 2017 !"
-    WRITE(*,*) " I'll stop here ... "
-    CALL abort
-  ENDIF
-    
-  factz = ( ptop - pver(itop) ) / ( pver(itop+1) - pver(itop) )
-  
-  zptop_ver = rver(itop)*(1.-factz) + rver(itop+1)*factz
-  
-  ! ---------------------------------------------------------
-  ! 4. Find zkim max assuming dz=10km and hence nlaykim_up
-  ! ---------------------------------------------------------
-  
-  zkim = zptop_ver
-  i=0
-  
-  DO WHILE ( zkim.lt.rver(131) )
-
-    zkim = zkim + grkim_dz
-    i=i+1
-  
-  ENDDO
-  
-  ! We want the ceiling at 1300km sharp, so we will either crop
-  ! the last layer or remove it and enlarge the "n-1"th
-  IF ( zkim .GT. rver(131)+0.5*grkim_dz ) THEN 
-    nlaykim_up = i-1
-    crop = .FALSE.
-  ELSE
-    nlaykim_up = i
-    crop = .TRUE.
-  ENDIF
-  
-  ! -----------------------------------------------------------------
-  ! 5. Calculates preskim grid interpolating back on Vervack profile
-  ! -----------------------------------------------------------------
-
-  ALLOCATE(preskim(nlaykim_up))
-  
-  jold=2
-
-  zkim = zptop_ver + 0.5*grkim_dz ! We want values at mid-layer here !
-  
-  DO i=1,nlaykim_up
-  
-    DO j=jold,131
-      IF ( zkim .GT. rver(j) ) THEN
-        jlay = j
-      ENDIF
-    ENDDO
-
-    jold = jlay ! keep in memory where we were in the above loop
-
-    ! We want the ceiling at 1300km sharp, we readjust the size of last layer
-    IF (i.eq.nlaykim_up) THEN
-      zkim = 0.5* ( zkim - 0.5*grkim_dz +  rver(131) )
-    ENDIF
-    
-    factz = ( zkim - rver(jlay) ) / ( rver(jlay+1) - rver(jlay) )
-    preskim(i) = pver(jlay)**(1.-factz) * pver(jlay+1)**factz 
-    
-    zkim = zkim + grkim_dz 
-
-  ENDDO
-
-  
-END SUBROUTINE gr_kim_vervack
Index: /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calchim.F90	(revision 1908)
@@ -1,488 +1,542 @@
-      SUBROUTINE calchim(ngrid,qy_c,declin_rad,ls_rad, &
-                         dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc)
-      
-!--------------------------------------------------------------------
-!  
-!     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 comchem_h
-      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
+SUBROUTINE calchim(ngrid,qy_c,declin,dtchim,            &
+     ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
+
+  !---------------------------------------------------------------------------------------------------------
+  !  
+  ! Purpose : Interface subroutine to photochemical C model for Titan GCM.
+  ! -------
+  !           The subroutine computes the chemical processes for a single vertical column.
+  !
+  ! - Only tendencies are returned.
+  ! - With moyzon_ch=.true. and input vectors zonally averaged
+  !   the calculation is done only once per lat. band
+  !
+  ! Authors: + S. Lebonnois : 01/2000 | 09/2003 
+  ! -------                  adaptation for Titan 3D : 02/2009
+  !                          adaptation for // : 04/2013
+  !                          extension chemistry up to 1300km : 10/2013
+  !      
+  !          + J. Vatant d'Ollone
+  !               + 02/17 - adaptation for the issu du generic
+  !               + 01/18 - 03/18 - Major transformations :
+  !                   - Upper chemistry fields are now stored in startfi 
+  !                     and defined on a pressure grid from Vervack profile
+  !                   - Altitudes above GCM top are computed using correct g(z) 
+  !                     to be coherent with obs and with dz=10km for UV processes
+  !                     but this should be even better if all physiq "under" was doing this !
+  !                   - These modifs enables to run chemistry with others resolution than 32x48x55 !
+  !                   - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
+  !                   - Chemistry can still be done in 2D 
+  !                      -> Calcul. once per band lat and put same tendency in all longi. 
+  !                         Check for negs in physiq_mod.
+  !                      -> If procs sharing a lat band, no problem, the calcul will just be done twice.
+  !                      -> Will not work with Dynamico, where the chemistry will have to be done in 3D.
+  !                          ( and there'll be work to do to get rid of averaged fields )
+  !
+  ! + STILL TO DO : Replug the interaction with haze (cf titan.old) -> to see with JB.
+  !---------------------------------------------------------------------------------------------------------
+
+  ! --------------------------------------------------------------------
+  ! Structure :
+  ! -----------
+  !   0.  Declarations
+  !   I.  Init and firstcall
+  !         1. Read and store Vervack profile
+  !         2. Compute planetar averaged atm. properties
+  !         3. Init compound caracteristics
+  !         4. Init photodissociations rates from actinic fluxes
+  !         5. Init chemical reactions
+  !         6. Init eddy diffusion coeff
+  !   II. Loop on latitudes/grid-points
+  !         0. Check on 2D chemistry
+  !           1. Compute atm. properties at grid points
+  !           2. Interpolate photodissociation rates at lat,alt,dec
+  !           3. Read composition
+  !           4. Call main solver gptitan C routine
+  !           5. Calculate output tendencies on advected tracers
+  !           6. Update upper chemistry fields
+  !         0bis. If 2D chemsitry, don't recalculate if needed
+  ! -----------------------------------------------------------------
+
+  USE comchem_h
+  USE dimphy
+  USE datafile_mod, ONLY: datadir
+  USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
+  USE geometry_mod, ONLY: latitude
+  USE logic_mod, ONLY: moyzon_ch
+  USE moyzon_mod, ONLY: tmoy, playmoy, phimoy, zlaymoy, zlevmoy
+
+  IMPLICIT NONE
 
 ! ------------------------------------------
 ! *********** 0. Declarations *************
 ! ------------------------------------------
-      
-!    Arguments
-!    ---------
-
-      INTEGER      ngrid                   ! nb of horiz points
-      REAL*8         qy_c(ngrid,klev,nkim)     ! Especes chimiques apres adv.+diss.
-      REAL*8         declin_rad,ls_rad      ! declinaison et long solaire en radians
-      REAL*8         dtchim                 ! pas de temps chimie
-      REAL*8         ctemp(ngrid,klev)       ! Temperature
-      REAL*8         cplay(ngrid,klev)       ! pression (Pa)
-      REAL*8         cplev(ngrid,klev+1)     ! pression intercouches (Pa)
-      REAL*8         czlay(ngrid,klev)       ! altitude (m)
-      REAL*8         czlev(ngrid,klev+1)     ! altitude intercouches (m)
-      
-      REAL*8         dqyc(ngrid,klev,nkim)     ! Tendances especes chimiques
-      
-      
-!    Local variables :
-!    -----------------
-
-      integer i,j,l,ic,jm1
-
-! variables envoyees dans la chimie: double precision
-
-      REAL*8  temp_c(nlaykim_tot)
-      REAL*8  press_c(nlaykim_tot)   ! T,p(mbar) a 1 lat donnee
-      REAL*8  cqy(nlaykim_tot,nkim)    ! frac mol qui seront modifiees
-      REAL*8  cqy0(nlaykim_tot,nkim)    ! frac mol avant modif
-      REAL*8  surfhaze(nlaykim_tot)
-      REAL*8  cprodaer(nlaykim_tot,4),cmaer(nlaykim_tot,4)
-      REAL*8  ccsn(nlaykim_tot,4),ccsh(nlaykim_tot,4)
-! rmil: milieu de couche, grille pour K,D,p,ct,T,y
-! rinter: intercouche (grille RA dans la chimie)
-      REAL*8  rmil(nlaykim_tot),rinter(nlaykim_tot),nb(nlaykim_tot)
-      REAL*8,save,allocatable :: kedd(:)
-
-      character str1*1,str2*2
-      
-!     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(nkim)
-      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_kim),nom_prod(nkim),nom_perte(nkim)
-      integer,save :: prod(200,nkim),perte(2,200,nkim)
-      character  fich*7,ficdec*15,curdec*15,name*10
-      real*8  ficalt,oldq,newq,zalt
-      logical okfic
-
-!     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
-
-!-----------------------------------------------------------------------
-!***********************************************************************
-!
-!    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(nlaykim_tot))
-
-        allocate(krpd(15,nd_kim+1,nlrt_kim,nbp_lat),krate(nlaykim_tot,nr_kim),md(nlaykim_tot,nkim))
-
-
-! 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'
+
+  ! Arguments
+  ! ---------
+
+  INTEGER, INTENT(IN)                              :: ngrid       ! Number of atmospheric columns.
+  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(IN)   :: qy_c        ! Chemical species on GCM layers after adv.+diss. (mol/mol).
+  REAL*8, INTENT(IN)                               :: declin      ! Solar declination (rad).
+  REAL*8, INTENT(IN)                               :: dtchim      ! Chemistry timsetep (s).
+  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: ctemp       ! Mid-layer temperature (K).
+  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cpphi       ! Mid-layer geopotential (m2.s-2).
+  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cplay       ! Mid-layer pressure (Pa).
+  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: cplev       ! Inter-layer pressure (Pa).
+  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: czlay       ! Mid-layer altitude (m).
+  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: czlev       ! Inter-layer altitude (m).
+
+  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT)  :: dqyc        ! Chemical species tendencies on GCM layers (mol/mol/s).
+
+  ! Local variables :
+  ! -----------------
+
+  INTEGER :: i , l, ic, ig, igm1
+
+  INTEGER :: dec, idec, ipres, ialt, klat
+
+  REAL*8  :: declin_c  ! Declination (deg).
+  REAL*8  :: factp, factalt, factdec, factlat, krpddec, krpddecp1, krpddecm1
+  REAL*8  :: temp1, logp
+
+  ! Variables sent into chemistry module (must be in double precision)
+  ! ------------------------------------------------------------------
+
+  REAL*8, DIMENSION(nlaykim_tot) :: temp_c  ! Temperature (K).
+  REAL*8, DIMENSION(nlaykim_tot) :: press_c ! Pressure (Pa).
+  REAL*8, DIMENSION(nlaykim_tot) :: phi_c   ! Geopotential (m2.s-2) - actually not sent in chem. module but used to compute alts.
+  REAL*8, DIMENSION(nlaykim_tot) :: nb      ! Density (cm-3).
+
+  REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy   ! Chemical species in whole column (mol/mol) sent to chem. module.
+  REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy0  !     "      "     "   "      "        "     before modifs.
+
+  REAL*8  :: surfhaze(nlaykim_tot)
+  REAL*8  :: cprodaer(nlaykim_tot,4), cmaer(nlaykim_tot,4)
+  REAL*8  :: ccsn(nlaykim_tot,4), ccsh(nlaykim_tot,4)
+
+  REAL*8, DIMENSION(nlaykim_tot) :: rmil   ! Mid-layer distance (km) to planetographic center.
+  REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
+  ! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
+
+  ! Saved variables initialized at firstcall
+  ! ----------------------------------------
+
+  LOGICAL, SAVE :: firstcall = .TRUE.
+!$OMP THREADPRIVATE(firstcall)
+
+  REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: kedd ! Eddy mixing coefficient for upper chemistry (cm^2.s-1)
+!$OMP THREADPRIVATE(kedd)
+
+  REAL*8, DIMENSION(:),   ALLOCATABLE, SAVE :: mass ! Molar mass of the compunds (g.mol-1)
+  REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: md
+!$OMP THREADPRIVATE(mass,md)
+
+  REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: r1d, ct1d, p1d, t1d ! Vervack profile
+  ! JVO 18 : No threadprivate for those as they'll be read in tcp.ver by master 
+
+  REAL*8, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: krpd  ! Photodissociations rate table
+  REAL*8, DIMENSION(:,:)    , ALLOCATABLE, SAVE :: krate ! Reactions rate ( photo + chem )
+!$OMP THREADPRIVATE(krpd,krate)
+
+  INTEGER, DIMENSION(:),     ALLOCATABLE, SAVE :: nom_prod, nom_perte
+  INTEGER, DIMENSION(:,:),   ALLOCATABLE, SAVE :: reactif
+  INTEGER, DIMENSION(:,:),   ALLOCATABLE, SAVE :: prod
+  INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: perte
+!$OMP THREADPRIVATE(nom_prod,nom_perte,reactif,prod,perte)
+
+  ! TEMPORARY : Dummy parameters without microphysics
+  ! Here to keep the whole stuff running without modif chem. module
+  ! ---------------------------------------------------------------
+
+  INTEGER  :: utilaer(16)
+  INTEGER  :: aerprod = 0
+  INTEGER  :: htoh2   = 0
+
+  ! -----------------------------------------------------------------------
+  ! ***************** I. Initialisations and Firstcall ********************
+  ! -----------------------------------------------------------------------
+
+  IF (firstcall) THEN
+
+     PRINT*, 'CHIMIE, premier appel'
+
+     call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
+
+     ALLOCATE(r1d(131))
+     ALLOCATE(ct1d(131))
+     ALLOCATE(p1d(131))
+     ALLOCATE(t1d(131))
+
+     ALLOCATE(md(nlaykim_tot,nkim))
+     ALLOCATE(mass(nkim))
+     
+     ALLOCATE(kedd(nlaykim_tot))
+     
+     ALLOCATE(krate(nlaykim_tot,nr_kim))
+     ALLOCATE(krpd(15,nd_kim+1,nlrt_kim,nlat_actfluxes))
+
+     ALLOCATE(nom_prod(nkim))
+     ALLOCATE(nom_perte(nkim))
+     ALLOCATE(reactif(5,nr_kim))
+     ALLOCATE(prod(200,nkim))
+     ALLOCATE(perte(2,200,nkim))
+     
+     ! 1. Read Vervack profile "tcp.ver", once for all
+     ! -----------------------------------------------
+
+!$OMP MASTER
+     OPEN(11,FILE=TRIM(datadir)//'/tcp.ver',STATUS='old')
+     READ(11,*)
+     DO i=1,131
+        READ(11,*) r1d(i), t1d(i), ct1d(i), p1d(i)
+        ! For debug :
+        ! -----------
+        ! PRINT*, "tcp.ver", r1d(i), t1d(i), ct1d(i), p1d(i)
+     ENDDO
+     CLOSE(11)
+!$OMP END MASTER
+!$OMP BARRIER
+
+     ! 2. Calculation of temp_c, densities and altitudes in planetary average
+     ! ----------------------------------------------------------------------
+     
+     ! JVO18 : altitudes calculated in firstcall are just for diag, as I set kedd in pressure grid
+
+     ! a. For GCM layers we just copy-paste
+
+!$OMP MASTER
+     PRINT*,'Init chemistry : pressure, density, temperature ... :'
+     PRINT*,'level, rmil (km), press_c (mbar), nb (cm-3), temp_c(K)'
+!$OMP END MASTER
+!$OMP BARRIER
+     
+     IF (ngrid.NE.1) THEN
+       DO l=1,klev
+          rinter(l)  = (zlevmoy(l)+rad)/1000.0              ! km
+          rmil(l)    = (zlaymoy(l)+rad)/1000.0              ! km
+          temp_c(l)  = tmoy(l)                              ! K
+          phi_c(l)   = phimoy(l)                            ! m2.s-2
+          press_c(l) = playmoy(l)/100.                      ! mbar
+          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
+!$OMP MASTER
+          PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
+!$OMP END MASTER
+!$OMP BARRIER
+       ENDDO
+       rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
+     ELSE
+       DO l=1,klev
+          rinter(l)  = (czlev(1,l)+rad)/1000.0                ! km
+          rmil(l)    = (czlay(1,l)+rad)/1000.0                ! km
+          temp_c(l)  = ctemp(1,l)                             ! K 
+          phi_c(l)   = cpphi(1,l)                             ! m2.s-2
+          press_c(l) = cplay(1,l)/100.                        ! mbar
+          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
+          PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
+       ENDDO
+       rinter(klev+1)=(czlev(1,klev+1)+rad)/1000.
+     ENDIF ! if ngrid.ne.1
+
+     ! b. Extension in upper atmosphere with Vervack profile
+     !    NB : Maybe the transition klev/klev+1 is harsh ...     
+
+     ipres=1
+     DO l=klev+1,nlaykim_tot
+        press_c(l) = preskim(l-klev) / 100.0
+        DO i=ipres,130
+           IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
+              ipres=i
+           ENDIF
+        ENDDO
+        factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
+
+        nb(l)      = exp( log(ct1d(ipres))*(1-factp) + log(ct1d(ipres+1))* factp )
+        temp_c(l)  = t1d(ipres)*(1-factp) + t1d(ipres+1)*factp
+     ENDDO
+
+     ! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
+     ! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
+
+     DO l=klev+1,nlaykim_tot
+
+        ! Compute geopotential on the upper grid, to have correct altitude
+        ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
+        ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
+
+        temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
+        phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
+     
+        rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
+
+!$OMP MASTER
+        PRINT*, l , rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
+!$OMP END MASTER
+!$OMP BARRIER
+     ENDDO
+
+     DO l=klev+2,nlaykim_tot
+        rinter(l) = 0.5*(rmil(l-1) + rmil(l))
+     ENDDO
+
+     ! 3. Compounds caracteristics
+     ! ---------------------------
+     mass(:) = 0.0
+     call comp(nomqy_c,nb,temp_c,mass,md)
+     PRINT*,'           Mass'
+     DO ic=1,nkim
+        PRINT*, nomqy_c(ic), mass(ic)
+     ENDDO
+
+     ! 4. Photodissociation rates
+     ! --------------------------
+     call disso(krpd,nlat_actfluxes) 
+
+     ! 5. Init. chemical reactions with planetary average T prof.
+     ! ----------------------------------------------------------
+
+     !  NB : Chemical reactions rate are assumed to be constant within the T range of Titan's atm
+     !  so we fill their krate once for all but krate for photodiss will be filled at each timestep
+     
+     call chimie(nomqy_c,nb,temp_c,krate,reactif, &
+          nom_perte,nom_prod,perte,prod)
+
+     ! 6. Eddy mixing coefficients (constant with time and space)
+     ! ----------------------------------------------------------
+
+     kedd(:) = 1.e3 ! Default value =/= zero
+
+     ! NB : Eddy coeffs from 1D models (e.g. Lavvas et al 08) in altitude but they're rather linked to pressure
+     !      Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
+
+     ! First calculate kedd for upper chemistry layers
+     DO l=klev+1,nlaykim_tot
+        logp=-log10(press_c(l))
+     ! 1E6 at 300 km ~ 10-1 mbar
+        IF     ( logp.ge.1.0 .and. logp.le.4.0 ) THEN 
+              kedd(l) = 10.**(6.0+1.1*(logp-1.0))
+     ! 2E7 at 600 km ~ 10-4 mbar
+        ELSEIF ( logp.gt.4.0 .and. logp.le.6.0 ) THEN
+              kedd(l) = 10.**(7.3+0.35*(logp-4.0))
+     ! 1E8 above 900 km ~ 10-6 mbar
+        ELSEIF ( logp.gt.6.0                   ) THEN 
+             kedd(l) = 1.e8
+        ENDIF
+     ENDDO
+     
+     ! Then adjust 15 last layers profile fading to default value depending on kedd(ptop)
+     DO l=klev-15,klev
+        temp1   = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev+1)/press_c(klev-15)) )
+        kedd(l) = 10.**( 3.0 + log10(kedd(klev+1)/1.e3) * temp1 )
+     ENDDO
+     
+     firstcall = .FALSE.
+  ENDIF  ! firstcall
+  
+  declin_c = declin*180./pi
+
+  ! -----------------------------------------------------------------------
+  ! *********************** II. Loop on latitudes *************************
+  ! -----------------------------------------------------------------------
+  
+  DO ig=1,ngrid 
+
+    IF (ig.eq.1) THEN 
+        igm1=1
+     ELSE
+        igm1=ig-1
+     ENDIF
+
+     ! If 2D chemistry, trick to do the calculation only once per latitude band within the chunk
+     ! NB1 : Will be obsolete with DYNAMICO, the chemistry will necessarly be 3D
+     ! NB2 : Test of same latitude with dlat=0.5 : I think that if you run sims better than 1/10th degree then
+     ! either it's with Dynamico and doesn't apply OR it is more than enough in terms of "preco / calc time" !
+     ! -------------------------------------------------------------------------------------------------------
+
+     IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.)) ) .OR. (.NOT. moyzon_ch) ) THEN
+
+        ! 1. Compute altitude for the grid point with hydrostat. equilib.
+        ! ---------------------------------------------------------------
+
+        ! a. For GCM layers we just copy-paste
+
         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(nlaykim_tot)=1290+2575 km.
-       do l=klev+2,nlaykim_tot
-         rinter(l) = rinter(klev+1) &
-               + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
-         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
-       enddo
-       rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-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
-! remplissage r1d,t1d,ct1d,p1d
-        open(11,file=TRIM(datadir)//'/tcp.ver',status='old')
-        read(11,*)
-        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 nlaykim_tot avec tcp.ver
-! la jonction klev/klev+1 est brutale... Tant pis ?
-        ialt=1
-        do l=klev+1,nlaykim_tot
-           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,nkim
-          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,nkim
-!          print*,nom_prod(ic),nom_perte(ic)
-!        enddo
-!        print*,'premieres prod, perte(1:reaction,2:compagnon):'
-!        do ic=1,nkim
-!          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_kim+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_kim+2,nr_kim
-!         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...
+           rinter(l)  = (czlev(ig,l)+rad)/1000.0             ! km
+           rmil(l)    = (czlay(ig,l)+rad)/1000.0             ! km
+           temp_c(l)  = ctemp(ig,l)                          ! K
+           phi_c(l)   = cpphi(ig,l)                          ! m2.s-2
+           press_c(l) = cplay(ig,l)/100.                     ! mbar
+           nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
+        ENDDO
+        rinter(klev+1)=(czlev(ig,klev+1)+rad)/1000.
+
+        ! b. Extension in upper atmosphere with Vervack profile
+        !    NB : The transition klev/klev+1 is maybe harsh if we have correct g above but not under  
+
+        ipres=1
+        DO l=klev+1,nlaykim_tot
+          press_c(l) = preskim(l-klev) / 100.0
+          DO i=ipres,130
+              IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
+                ipres=i
+              ENDIF
+          ENDDO
+          factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
+
+          nb(l)      = exp( log(ct1d(ipres))*(1-factp) + log(ct1d(ipres+1))* factp )
+          temp_c(l)  = t1d(ipres)*(1-factp) + t1d(ipres+1)*factp
+        ENDDO
  
-      kedd(:) = 5e2 ! valeur mise par defaut  
-               ! UNITE ? doit etre ok pour gptitan
-      do l=1,nlaykim_tot
-       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
-!
-!                * Permet de faire le calcul une seule fois par lat
-!
-      DO j=1,ngrid 
-      
-      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(nlaykim_tot)=1290+2575 km.
-       do l=klev+2,nlaykim_tot
-         rinter(l) = rinter(klev+1)  &
-               + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
-         rmil(l-1) = (rinter(l-1)+rinter(l))/2.
-       enddo
-       rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-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 nlaykim_tot avec tcp.ver
-       ialt=1
-       do l=klev+1,nlaykim_tot
-           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,nlaykim_tot
-
-! 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_kim+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_kim+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,nkim
-        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 nlaykim_tot
-
-      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,nkim
-        read(11,'(A10)') name
-        if (name.ne.nomqy_c(ic)) then
-          print*,"probleme lecture ",fich
-          print*,name,nomqy_c(ic)
-          stop
+        ! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
+        ! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
+
+        DO l=klev+1,nlaykim_tot
+
+           ! Compute geopotential on the upper grid, to have correct altitude
+           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
+           ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
+
+           temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
+           phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
+
+           rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
+        ENDDO
+
+        DO l=klev+2,nlaykim_tot
+          rinter(l) = 0.5*(rmil(l-1) + rmil(l))
+        ENDDO
+
+        ! 2. From krpd, compute krate for dissociations (declination-latitude-altitude interpolation)
+        ! -------------------------------------------------------------------------------------------
+
+        ! a. Calculate declination dependence
+
+        if ((declin_c*10+267).lt.14.) then
+           idec = 0
+           dec  = 0
+        else 
+           if ((declin_c*10+267).gt.520.) then
+              idec = 14
+              dec  = 534
+           else 
+              idec = 1
+              dec  = 27
+              do while( (declin_c*10+267).ge.real(dec+20) )
+                 dec  = dec+40
+                 idec = idec+1
+              enddo
+           endif
         endif
-        if (ficdec.eq.curdec) then
-          do l=klev+1,nlaykim_tot
-            read(11,*) ficalt,cqy(l,ic),newq
-          enddo
+        if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
+           factdec = ( declin_c - (dec-267)/10. ) / 4.
         else
-          do l=klev+1,nlaykim_tot
-            read(11,*) ficalt,oldq,cqy(l,ic)
-          enddo
+           factdec = ( declin_c - (dec-267)/10. ) / 2.7
         endif
-       enddo
-       close(11)
-      else       ! le fichier n'est pas la
-       do ic=1,nkim
-        do l=klev+1,nlaykim_tot
-          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,nkim
-        do l=1,klev
-          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
-        enddo
-       enddo
-
-
-!-----------------------------------------------------------------------
-!
-!   sauvegarde compo sur nlaykim_tot 
-!   -------------------------
-
-! 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,nkim
-       write(11,'(A10)') nomqy_c(ic)
-       do l=klev+1,nlaykim_tot
-        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, only adjust longitudinal variations
-        dqyc(j,:,:) = dqyc(jm1,:,:)/qy_c(jm1,:,:)*qy_c(j,:,:)
-      endif
-
-      ENDDO
-     
-!***********************************************************************
-!***********************************************************************
-
-      firstcall = .false.
-      
-      end subroutine calchim
+
+        ! b. Calculate klat for interpolation on fixed latitudes of actinic fluxes input
+
+        klat=1
+        DO i=1,nlat_actfluxes
+          IF (latitude(ig).LT.lat_actfluxes(i)) klat=i
+        ENDDO
+        IF (klat==nlat_actfluxes) THEN ! avoid rounding problems
+          klat    = nlat_actfluxes-1
+          factlat = 1.0
+        ELSE
+          factlat = (latitude(ig)-lat_actfluxes(klat))/(lat_actfluxes(klat+1)-lat_actfluxes(klat))
+        ENDIF
+
+        ! c. Altitude loop
+
+        DO l=1,nlaykim_tot
+
+           ! Calculate ialt for interpolation in altitude (krpd every 2 km)
+           ialt    = int((rmil(l)-rad/1000.)/2.)+1
+           factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
+
+           ! Altitude can go above top limit of UV levels - in this case we keep the 1310km top fluxes
+           IF (ialt.GT.nlrt_kim-1) THEN
+             ialt     = nlrt_kim-1 ! avoid out-of-bound array
+             factalt  = 1.0
+           ENDIF
+
+           DO i=1,nd_kim+1
+
+              krpddecm1 =   (   krpd(idec  ,i,ialt  ,klat)   * (1.0-factalt)                   &
+                              + krpd(idec  ,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                          * (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                              + krpd(idec  ,i,ialt+1,klat+1) * factalt       ) * factlat
+              krpddec   =   (   krpd(idec+1,i,ialt  ,klat)   * (1.0-factalt)                   &
+                              + krpd(idec+1,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                          * (   krpd(idec+1,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                              + krpd(idec+1,i,ialt+1,klat+1) * factalt       ) * factlat
+              krpddecp1 =   (   krpd(idec+2,i,ialt  ,klat)   * (1.0-factalt)                   &
+                              + krpd(idec+2,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
+                          * (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
+                              + krpd(idec+2,i,ialt+1,klat+1) * factalt       ) * factlat
+
+              ! nd_kim+1 is dissociation of N2 by GCR
+              if ( factdec.lt.0. ) then 
+                 krate(l,i) = krpddecm1 * abs(factdec) + krpddec   * ( 1.0 + factdec)
+              endif
+              if ( factdec.gt.0. ) then
+                 krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)
+              endif
+              if ( factdec.eq.0. ) krate(l,i) = krpddec
+
+           ENDDO ! i=1,nd_kim+1
+        ENDDO ! l=1,nlaykim_tot
+
+        ! 3. Read composition 
+        ! -------------------
+
+        DO ic=1,nkim
+           DO l=1,klev
+              cqy(l,ic) = qy_c(ig,l,ic) ! advected tracers for the GCM part converted to molar frac.
+           ENDDO
+           
+           DO l=1,nlaykim_up
+              cqy(klev+l,ic) = ykim_up(ic,ig,l) ! ykim_up for the upper atm.
+           ENDDO
+        ENDDO
+
+        cqy0(:,:) = cqy(:,:) ! Stores compo. before modifs
+
+        ! 4. Call main Titan chemistry C routine
+        ! --------------------------------------
+
+        call gptitan(rinter,temp_c,nb,                  &
+             nomqy_c,cqy,                               &
+             dtchim,latitude(ig)*180./pi,mass,md,       &
+             kedd,botCH4,krate,reactif,                 &
+             nom_prod,nom_perte,prod,perte,             &
+             aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,  &
+             htoh2,surfhaze)
+
+        ! 5. Calculates tendencies on composition for advected tracers 
+        ! ------------------------------------------------------------
+        DO ic=1,nkim
+           DO l=1,klev
+              dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
+           ENDDO
+        ENDDO
+
+        ! 6. Update ykim_up
+        ! -----------------
+        DO ic=1,nkim
+           DO l=1,nlaykim_up
+              ykim_up(ic,ig,l) = cqy(klev+l,ic)
+           ENDDO
+        ENDDO
+        ! NB: The full vertical composition grid will be created only for the outputs
+
+
+     ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again ! 
+        dqyc(ig,:,:)    = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
+        ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations
+     ENDIF
+
+  ENDDO
+
+END SUBROUTINE calchim
Index: /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1908)
@@ -1,8 +1,8 @@
 
 
-SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdq)
+SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdqfi, zdq)
   !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
   !!
-  !! The subroutine computes the microphysics processes for a si:le vertical column.
+  !! 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.
@@ -30,6 +30,9 @@
   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(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
+  REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(kg.kg^{-1}}\)).
+  REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(kg.kg^{-1}}\)).
+  
+  REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: zq !! Local tracers updated from former processes (\(kg.kg^{-1}}\)).
   
   REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as !! 0th order moment of the spherical mode (\(m^{-2}\)).
@@ -56,7 +59,8 @@
 
   INTEGER :: ilon, i,nices
-  INTEGER :: nlon,nlay
+  INTEGER :: nq,nlon,nlay
 
   ! Read size of arrays
+  nq    = size(pq,DIM=3)
   nlon  = size(play,DIM=1)
   nlay  = size(play,DIM=2)
@@ -65,7 +69,7 @@
   ALLOCATE( int2ext(nlay) )  
 
-  ! Loop on horizontal grid points
+  ! Allocate arrays
+  ALLOCATE( zq(nlon,nlay,nq) )
 
-  ! Allocate arrays
   ALLOCATE( m0as(nlay) )
   ALLOCATE( m3as(nlay) )
@@ -89,19 +93,23 @@
   zdq(:,:,:) = 0.0
 
+  ! Initialize tracers updated with former processes from physics
+  zq(:,:,:) = pq(:,:,:) + zdqfi(:,:,:)*dt
+
+  ! Loop on horizontal grid points
   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,1:nlay)-plev(ilon,2:nlay+1) ) / g
-    m0as(:) = pq(ilon,:,1) * int2ext(:)
-    m3as(:) = pq(ilon,:,2) * int2ext(:)
-    m0af(:) = pq(ilon,:,3) * int2ext(:)
-    m3af(:) = pq(ilon,:,4) * int2ext(:)
+    m0as(:) = zq(ilon,:,1) * int2ext(:)
+    m3as(:) = zq(ilon,:,2) * int2ext(:)
+    m0af(:) = zq(ilon,:,3) * int2ext(:)
+    m3af(:) = zq(ilon,:,4) * int2ext(:)
     
     if (callclouds) then ! if call clouds
-      m0n(:) = pq(ilon,:,5) * int2ext(:)
-      m3n(:) = pq(ilon,:,6) * int2ext(:)
+      m0n(:) = zq(ilon,:,5) * int2ext(:)
+      m3n(:) = zq(ilon,:,6) * int2ext(:)
       do i=1,nices
-        m3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
-        gazs(:,i)    = pq(ilon,:,ices_indx(i)) / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
+        m3i(:,nices) = zq(ilon,:,6+i) * int2ext(:)
+        gazs(:,i)    = zq(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
Index: /trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/phytitan/comchem_h.F90	(revision 1908)
@@ -1,10 +1,12 @@
 MODULE comchem_h
 
-! ----------------------------------------------------------------------------
-! Purpose : Stores data relative to chemistry in the GCM and upper chemistry
-! ------- 
-!           Please note that upper fields ykim_up are in molar fraction !
+! -------------------------------------------------------------------------------------------------------------
+! Purpose : + Stores data relative to :  * 1. Chemistry in the GCM
+! -------                                * 2. Upper chemistry pressure grid
+!                                        * 3. Coupling with C photochem. module ( cf calchim.F90)
+!           
+!           + Also contains a routine of initialization for chemistry in the GCM.
 !
-!          NB : For newstart there is a specific comchem_newstart_h module.
+!           + NB : For newstart there is a specific comchem_newstart_h module.
 !
 ! Author : Jan Vatant d'Ollone (2017-18)
@@ -12,12 +14,19 @@
 !
 ! NB : A given order is assumed for the 44 chemistry tracers :
-! --   H, H2, CH, CH2s, CH2, CH3, CH4, C2, C2H, C2H2, C2H3, C2H4, C2H5, 
+!      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
-        
-! ----------------------------------------------------------------------------
+!
+!        
+! IMPORTANT : Advected chem. tracers are in MASS fraction but upper fields ykim_up are in MOLAR fraction !
+! 
+! -------------------------------------------------------------------------------------------------------------
 
 IMPLICIT NONE  
+
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~
+   ! 1. Chemistry in the GCM
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
    !! Hard-coded number of chemical species for Titan chemistry
@@ -50,30 +59,47 @@
        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   /)
+
+   !! Hard-coded molar fraction of surface methane
+   REAL, PARAMETER :: botCH4 = 0.048
    
-   ! !!!!!!!!!!!!!!!!!!!!
-   !  Upper chemistry 
-   ! !!!!!!!!!!!!!!!!!!!!
-   
-   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! These parameters as well as nkim above, MUST match titan.h in chimtitan
-   INTEGER, PARAMETER :: nd_kim   = 54     ! Number of photodissociations
-   INTEGER, PARAMETER :: nr_kim   = 377    ! Number of reactions in chemistry scheme
-   INTEGER, PARAMETER :: nlrt_kim = 650    ! For the UV rad. transf., 650 levels of 2 km
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   
+
+
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
+   !  2. Upper chemistry grid 
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
    INTEGER, SAVE :: nlaykim_up   ! Number of upper atm. layers for chemistry from GCM top to 4.5E-5 Pa (1300km)
    INTEGER, SAVE :: nlaykim_tot  ! Number of total layers for chemistry from surface to 4.5E-5 Pa (1300km)
 !$OMP_THREADPRIVATE(nlaykim_up,nlay_kim_tot)
 
-  ! NB : For the startfile we use nlaykim_up grid (upper atm) and for outputs we use nlaykim_tot grid (all layers)
+   ! NB : For the startfile we use nlaykim_up grid (upper atm) and for outputs we use nlaykim_tot grid (all layers)
+  
+   REAL*8, PARAMETER :: grkim_dz = 10.0 ! Vertical discretization of the upper chemsitry grid (km)
 
-   REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: preskim  ! Pressure (Pa) of upper chemistry (mid)-layers
-!$OMP_THREADPRIVATE(preskim)
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:)     :: preskim  ! Pressure (Pa) of upper chemistry (mid)-layers
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: zlay_kim ! Altitude (km) of all chemistry (mid)-layers
+   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_up  ! Upper chemistry fields (mol/mol)
+!$OMP_THREADPRIVATE(preskim,zlay_kim,ykim_up)
 
-   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: zlay_kim  ! Altitude (km) of all chemistry (mid)-layers
-!$OMP_THREADPRIVATE(zlay_kim)
 
-   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_up ! Upper chemistry fields (mol/mol)
-!$OMP_THREADPRIVATE(ykim_up)
+
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   !  3. Interface with photochemical module (cf calchim.F90)
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   
+   ! These 3 parameters as well as nkim above, MUST match titan.h in chimtitan !!
+   INTEGER, PARAMETER :: nd_kim   = 54     ! Number of photodissociations
+   INTEGER, PARAMETER :: nr_kim   = 377    ! Number of reactions in chemistry scheme
+   INTEGER, PARAMETER :: nlrt_kim = 650    ! For the UV rad. transf., 650 levels of 2 km
+
+   !! Hardcoded latitude discretisation for actinic fluxes - MUST be coherent with disso.c input files !!
+   INTEGER, PARAMETER                         :: nlat_actfluxes = 49
+   REAL, DIMENSION(nlat_actfluxes), PARAMETER :: lat_actfluxes  = (/ &
+      90.0 ,  86.25,  82.5 ,  78.75,  75.0 ,  71.25,  67.5 ,  63.75,  60.0 ,  56.25,  52.5 ,  48.75, & 
+      45.0 ,  41.25,  37.5 ,  33.75,  30.0 ,  26.25,  22.5 ,  18.75,  15.0 ,  11.25,   7.5 ,   3.75, &
+       0.0 ,  -3.75,  -7.5 , -11.25, -15.0 , -18.75, -22.5 , -26.25, -30.0 , -33.75, -37.5 , -41.25, &
+     -45.0 , -48.75, -52.5 , -56.25, -60.0 , -63.75, -67.5 , -71.25, -75.0 , -78.75, -82.5 , -86.25, &
+     -90.0 /)
+    
 
 CONTAINS
Index: /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F	(revision 1908)
@@ -12,4 +12,5 @@
      &                     dtemisice
       use comdiurn_h, only: sinlat, coslat, sinlon, coslon
+      use comchem_h, only: nkim, nlaykim_up, ykim_up
       use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
       use phyredem, only: physdem0,physdem1
@@ -23,5 +24,5 @@
      &                            nday, iphysiq
       use callkeys_mod, only: tracer,check_cpp_match,rings_shadow,
-     &	                      specOLR,pceil
+     &	                      specOLR,pceil,callchim
       USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
      &                       presnivs,pseudoalt,scaleheight
@@ -670,6 +671,4 @@
          ENDDO
          
-
-
          DO ilayer=1,nlayer
 !     zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
@@ -762,4 +761,48 @@
         
       endif ! end if tracer
+      
+! Initialize upper atmosphere & chemistry
+! ---------------------------------------
+
+      ! Calculate the # of upper chemistry layers
+      call gr_kim_vervack
+      
+      write(*,*)
+      write(*,*) " With the compiled vertical grid we found :"
+      write(*,*) " Number of upper chemistry layers =", nlaykim_up
+        
+      if (callchim) then
+            
+        if (.NOT.allocated(ykim_up)) then
+          allocate(ykim_up(nkim,1,nlaykim_up))
+        endif
+      
+        ! Inquire if we have upper chemistry input profiles
+        
+        do iq=1,nkim
+        
+          file_prof = "profile_"//trim(noms(iq))//"_up"
+
+          INQUIRE(file=trim(file_prof),exist=findprof)
+        
+          IF (findprof) THEN
+            ! Read profile for upper chemistry specie.
+            PRINT *, "Input file found for upper chemistry specie ", 
+     &                trim(noms(iq))//"_up"
+            OPEN(11,file=trim(file_prof),status='old',form='formatted')
+            PRINT *, "---> reading upper profile from input file ..."
+            DO ilayer=1,nlaykim_up
+              READ (11,*) ykim_up(iq,1,ilayer)
+            ENDDO
+          ELSE
+            PRINT *, "No input profile found for ",
+     &                trim(noms(iq))//"_up"
+            PRINT *, "---> we set it to zero - that's bold !!"
+          ENDIF
+
+        enddo ! nkim loop
+        
+      endif ! of if callchim
+      
 
 c  Write a "startfi" file
Index: /trunk/LMDZ.TITAN/libf/phytitan/gr_kim_vervack.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/gr_kim_vervack.F90	(revision 1908)
+++ /trunk/LMDZ.TITAN/libf/phytitan/gr_kim_vervack.F90	(revision 1908)
@@ -0,0 +1,154 @@
+SUBROUTINE gr_kim_vervack
+
+  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !
+  ! Purpose : * Calculates the number of layers needed for upper chemistry
+  ! -------   based on the GCM vertical grid size, depending on Ptop=ap.
+  !           * Generates also the pressure grid at mid-layer for upper levels.
+  !           * We use an upper atmosphere profile based on Voyager 1 data 
+  !           (Vervack et al, 2004) to have dz=10km between Ptop and 1300km.
+  !
+  ! Author : Jan Vatant d'Ollone (2017)
+  ! ------
+  !
+  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  USE comchem_h, ONLY: nlaykim_up, preskim, grkim_dz
+  USE comvert_mod, ONLY: ap, aps, bp, bps, preff 
+  USE datafile_mod, ONLY: datadir
+
+  IMPLICIT NONE
+
+  INCLUDE "dimensions.h"
+
+  REAL*8  :: ptop, zptop_ver, zkim
+  REAL*8  :: factz
+
+  REAL*8  :: rver(131),tver(131),ctver(131),pver(131) 
+
+  INTEGER :: i, j, itop, jlay, jold
+  
+  LOGICAL :: foundver
+
+  LOGICAL :: crop
+
+  ! ------------------------
+  ! 1. Read Vervack profile
+  ! ------------------------
+
+  ! Check the file exists
+  INQUIRE(file=TRIM(datadir)//'/tcp.ver',exist=foundver)
+  IF(.NOT.foundver) THEN
+    WRITE(*,*)'The file ',TRIM(datadir)//'/tcp.ver'
+    WRITE(*,*)'was not found by gr_kim_vervack.F90, exiting'
+    WRITE(*,*)'Check that your path to datagcm:',trim(datadir)
+    WRITE(*,*)'is correct. You can change it in callphys.def with:'
+    WRITE(*,*)'datadir = /absolute/path/to/datagcm'
+    WRITE(*,*)'Also check that file tcp.ver exists there.'
+  ENDIF
+
+  ! Read file
+  OPEN(11,file=TRIM(datadir)//'/tcp.ver',status='old')
+  READ(11,*)
+  DO i=1,131
+     READ(11,*) rver(i),tver(i),ctver(i),pver(i)
+     pver(i) = pver(i)*100.0 ! mbar->Pa
+  ENDDO
+  CLOSE(11)
+  
+  ! --------------------------------------------------------------
+  ! 2. Define ptop as the value aps should have if it wasn't zero
+  ! assuming ap(llm)-aps(llm) half pressure thickness of top-layer
+  ! --------------------------------------------------------------
+  
+  ! NB : At the top of the model we are in pure pressure coord. -> ap
+  ! ( except for 1D where we have only bp )
+  
+  IF (jjm.GT.1) THEN
+    ptop = 2.0*aps(llm) - ap(llm)
+  ELSE
+    ptop = preff*(2.0*bps(llm) - bp(llm))
+  ENDIF
+  
+  ! --------------------------------------------
+  ! 3. Interpolate Ptop and equivalent altitude 
+  ! --------------------------------------------
+  
+  itop=1
+  
+  DO i=2,131
+    IF ( ptop .LT. pver(i) ) THEN
+      itop=i
+    ENDIF
+  ENDDO
+  
+  ! Crazy case in a far far away future where GCM top will reach 1300km
+  IF ( itop .eq. 131 ) THEN
+    WRITE(*,*) " You've reach the bounds of Vervack profile ... "
+    WRITE(*,*) " Congrats but it wasn't predicted in 2017 !"
+    WRITE(*,*) " I'll stop here ... "
+    CALL abort
+  ENDIF
+    
+  factz = ( ptop - pver(itop) ) / ( pver(itop+1) - pver(itop) )
+  
+  zptop_ver = rver(itop)*(1.-factz) + rver(itop+1)*factz
+  
+  ! ---------------------------------------------------------
+  ! 4. Find zkim max assuming dz=10km and hence nlaykim_up
+  ! ---------------------------------------------------------
+  
+  zkim = zptop_ver
+  i=0
+  
+  DO WHILE ( zkim.lt.rver(131) )
+
+    zkim = zkim + grkim_dz
+    i=i+1
+  
+  ENDDO
+  
+  ! We want the ceiling at 1300km sharp, so we will either crop
+  ! the last layer or remove it and enlarge the "n-1"th
+  IF ( zkim .GT. rver(131)+0.5*grkim_dz ) THEN 
+    nlaykim_up = i-1
+    crop = .FALSE.
+  ELSE
+    nlaykim_up = i
+    crop = .TRUE.
+  ENDIF
+  
+  ! -----------------------------------------------------------------
+  ! 5. Calculates preskim grid interpolating back on Vervack profile
+  ! -----------------------------------------------------------------
+
+  ALLOCATE(preskim(nlaykim_up))
+  
+  jold=2
+
+  zkim = zptop_ver + 0.5*grkim_dz ! We want values at mid-layer here !
+  
+  DO i=1,nlaykim_up
+  
+    DO j=jold,131
+      IF ( zkim .GT. rver(j) ) THEN
+        jlay = j
+      ENDIF
+    ENDDO
+
+    jold = jlay ! keep in memory where we were in the above loop
+
+    ! We want the ceiling at 1300km sharp, we readjust the size of last layer
+    IF (i.eq.nlaykim_up) THEN
+      zkim = 0.5* ( zkim - 0.5*grkim_dz +  rver(131) )
+    ENDIF
+    
+    factz = ( zkim - rver(jlay) ) / ( rver(jlay+1) - rver(jlay) )
+    preskim(i) = pver(jlay)**(1.-factz) * pver(jlay+1)**factz 
+    
+    zkim = zkim + grkim_dz 
+
+  ENDDO
+
+  
+END SUBROUTINE gr_kim_vervack
Index: /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90	(revision 1907)
+++ /trunk/LMDZ.TITAN/libf/phytitan/inicondens.F90	(revision 1908)
@@ -159,5 +159,3 @@
        enddo
 
-       print*, 'inicondens end'
-       
       end subroutine inicondens
