Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/gr_kim_vervack.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/gr_kim_vervack.F90	(revision 1889)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/gr_kim_vervack.F90	(revision 1889)
@@ -0,0 +1,151 @@
+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/dynphy_lonlat/phytitan/lect_start_archive.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 1887)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/lect_start_archive.F	(revision 1889)
@@ -5,4 +5,5 @@
 
 !      USE surfdat_h
+      USE comchem_h, ONLY: nlaykimold, preskimold
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
       USE infotrac, ONLY: tname, nqtot
@@ -46,5 +47,4 @@
 c------------------------------------
       INTEGER   imold,jmold,lmold,nsoilold,nqold
-
 
 c Variables pour les lectures des fichiers "ini" 
@@ -147,5 +147,5 @@
       real, dimension(:,:), allocatable :: emisold
       real, dimension(:,:,:,:), allocatable :: qold
-
+      
       real tab_cntrl(100)
 
@@ -158,4 +158,5 @@
       logical :: therminertia_3D=.true. ! flag
 ! therminertia_3D=.true. if thermal inertia is 3D and read from datafile
+
 c Variable intermediaires iutilise pour l'extrapolation verticale 
 c----------------------------------------------------------------
@@ -247,8 +248,15 @@
       endif
 
-
       if (nsoilold.ne.nsoilmx) then ! interpolation will be required
         depthinterpol=.true.
       endif
+
+! 1.2.2 find out the # of upper chemistry layers
+
+      ierr= NF_INQ_DIMID(nid,"upper_chemistry_layers",dimid)
+      ierr= NF_INQ_DIMLEN(nid,dimid,nlaykimold)
+
+      ! NB : The regriding, if needed cannot be done here since the new
+      ! pressure grid is only computed at the end of newstart
 
 ! 1.3 Report dimensions
@@ -266,4 +274,5 @@
       write(*,*) '  Otherwise, set nsoilmx -in dimphys.h- to: ',nsoilold
       endif
+      write(*,*) "upper_chemistry_layers: ",nlaykimold
       write(*,*) "time lenght: ",timelen
       write(*,*) 
@@ -295,4 +304,6 @@
       allocate(mlayerold(nsoilold))
       allocate(qsurfold(imold+1,jmold+1,nqtot))
+      
+      allocate(preskimold(nlaykimold))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -512,7 +523,26 @@
        endif
       endif
-
-c-----------------------------------------------------------------------
-c 3.6 Lecture geopotentiel au sol
+      
+c-----------------------------------------------------------------------
+c 3.6 Read upper chemistry mid-layer pressure
+c-----------------------------------------------------------------------  
+    
+      ierr=NF_INQ_VARID(nid,"preskim",nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <preskim> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, preskimold)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, preskimold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <preskim>"
+         CALL abort
+      ENDIF
+
+c-----------------------------------------------------------------------
+c 3.7 Lecture geopotentiel au sol
 c-----------------------------------------------------------------------
 c
Index: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F
===================================================================
--- trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1887)
+++ trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/newstart.F	(revision 1889)
@@ -17,4 +17,5 @@
      &                              is_master
       use infotrac, only: infotrac_init, nqtot, tname
+      USE comchem_h, ONLY: nlaykim_up, nlaykimold
       USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
       USE surfdat_h, ONLY: phisfi, albedodat,
@@ -130,5 +131,5 @@
       real tab_cntrl(100)
       real tab_cntrl_bis(100)
-
+      
 c variables diverses
 c-------------------
@@ -258,5 +259,4 @@
 
       endif
-
 
 c=======================================================================
@@ -536,6 +536,6 @@
         CALL lect_start_archive(ngridmx,llm,
      &   date,tsurf,tsoil,emis,q2,
-     &   t,ucov,vcov,ps,teta,phisold_newgrid,q,qsurf,
-     &   surfith,nid)
+     &   t,ucov,vcov,ps,teta,phisold_newgrid,
+     &   q,qsurf,surfith,nid)
         write(*,*) "OK, read start_archive file"
         ! copy soil thermal inertia
@@ -572,10 +572,4 @@
       write(*,*) 'q=x : give a specific uniform value to one tracer'
       write(*,*) 'q=profile    : specify a profile for a tracer'
-!      write(*,*) 'ini_q : tracers initialisation for chemistry, water an
-!     $d ice   '
-!      write(*,*) 'ini_q-H2O : tracers initialisation for chemistry and 
-!     $ice '
-!      write(*,*) 'ini_q-iceH2O : tracers initialisation for chemistry on
-!     $ly '
       write(*,*) 'isotherm  : Isothermal Temperatures, wind set to zero'
       write(*,*) 'radequi   : Earth-like radiative equilibrium temperature
@@ -1071,4 +1065,37 @@
       CALL inifilr 
       CALL pression(ip1jmp1, ap, bp, ps, p3d)
+      
+
+c=========================================================================
+c  Calcul de la dimension verticale pour la chimie  - JVO 2017
+c  start_archive seulement, la grille verticale pouvant avoir ete modifiee
+c==========================================================================
+
+      IF (choix_1.eq.0) THEN
+
+        ! Calculate the # of upper chemistry layers with the "new" pressure grid
+        ! For this we use Vervack profile for upper atmosphere with dz=10km
+
+        CALL gr_kim_vervack
+
+        WRITE(*,*)
+        WRITE(*,*) " With the compiled vertical grid we found :"
+        WRITE(*,*) " Number of upper chemistry layers =", nlaykim_up
+        
+        ! Regriding is then done, if needed
+        
+        IF (nlaykimold.ne.nlaykim_up) THEN
+
+          WRITE(*,*) " Warning, nlaykimold=", nlaykimold
+          WRITE(*,*) ' which implies that a regriding on upper chemistry
+     & will be performed.'
+          WRITE(*,*)
+          
+!          CALL  regrid_kim
+          
+        ENDIF
+ 
+      endif ! of if (choix_1.eq.0)      
+ 
 
 c-----------------------------------------------------------------------
@@ -1125,4 +1152,5 @@
          day_ini=int(date)
       endif
+
 c
       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
