SUBROUTINE calchim(ngrid,qy_c,declin,dtchim,            &
     ctemp,cpphi,cpphis,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 new generic-forked physics
  !               + 01/18 - 03/18 - Major transformations :
  !                   - Upper chemistry fields are now stored in startfi 
  !                     and defined on a pressure grid from Vervack profile
  !                   - 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 )
  !
  !               + 02/19 : To always have correct photodissociations rates, altitudes sent here by physiq are always 
  !                         calculated with effective g - and with reference to the body not the local surface -
  !                          even if in physiq we keep altitudes coherent with dynamics !
  !
  ! + STILL TO DO : + Replug the interaction with haze (cf titan.old) -> to see with JB.
  !                 + Use iso_c_binding for the fortran-C exchanges.
  !---------------------------------------------------------------------------------------------------------

  ! --------------------------------------------------------------------
  ! 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, INTRINSIC :: iso_c_binding
  USE comchem_h
  USE dimphy
  USE datafile_mod, ONLY: datadir
  USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
  USE geometry_mod, ONLY: latitude
#ifndef MESOSCALE  
  USE logic_mod, ONLY: moyzon_ch
  USE moyzon_mod, ONLY: tmoy, playmoy
#endif

  IMPLICIT NONE

! ------------------------------------------
! *********** 0. Declarations *************
! ------------------------------------------

  ! 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),           INTENT(IN)   :: cpphis      ! Surface 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 effective altitude (m) : ref = geoid.
  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: czlev       ! Inter-layer effective altitude (m) ref = geoid.

  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 :: md   ! Mean molecular diffusion coefficients (cm^2.s-1)
  REAL*8, DIMENSION(:),   ALLOCATABLE, SAVE :: mass ! Molar mass of the compounds (g.mol-1)
!$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'

     if (ngrid .eq. 1) then ! if 1D no dynamic mixing, we set the kedd in all column
       call check(nlaykim_tot,0,nlrt_kim,nkim)
     else
       call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
     endif

     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(nd_kim+1,nlrt_kim,15,nlat_actfluxes))

     ALLOCATE(nom_prod(nkim))
     ALLOCATE(nom_perte(nkim))
     ALLOCATE(reactif(5,nr_kim))
     ALLOCATE(prod(200,nkim))
     ALLOCATE(perte(2,200,nkim))
     
     ! 0. Deal with characters for C-interoperability
     ! ----------------------------------------------
     ! NB ( JVO 19 ) : Using iso_c_binding would do things in an even cleaner way !
     DO ic=1,nkim
       nomqy_c(ic) = trim(cnames(ic))//char(0) ! Add the C null terminator
     ENDDO
     nomqy_c(nkim+1)="HV"//char(0) ! For photodissociations

     ! 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 are no more calculated in firstcall, as I set kedd in pressure grid

     ! a. For GCM layers we just copy-paste ( assuming that physiq always send correct altitudes ! )

     PRINT*,'Init chemistry : pressure, density, temperature ... :'
     PRINT*,'level, press_c (mbar), nb (cm-3), temp_c (K)'

#ifndef MESOSCALE
!! JVO 20 : I put this CPP key because mesoscale cannot access tmoy and playmoy,
!  you should fix this when you will propeerly use chemistry !

     IF (ngrid.NE.1) THEN
       DO l=1,klev
          temp_c(l)  = tmoy(l)                              ! K
          press_c(l) = playmoy(l)/100.                      ! mbar
          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
          PRINT*, l, press_c(l), nb(l), temp_c(l)
       ENDDO
     ELSE
#endif
       DO l=1,klev
          temp_c(l)  = ctemp(1,l)                             ! K 
          press_c(l) = cplay(1,l)/100.                        ! mbar
          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
          PRINT*, l, press_c(l), nb(l), temp_c(l)
       ENDDO
#ifndef MESOSCALE
     ENDIF
#endif

     ! b. Extension in upper atmosphere with Vervack profile
     ! NB : Maybe the transition klev/klev+1 is harsh if T profile different from Vervack ...     

     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
        PRINT*, l , press_c(l), nb(l), temp_c(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 (e.g. Lavvas et al 08, Yelle 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-4,nlaykim_tot
     !   logp=-log10(press_c(l))
     !! 2E6 at 400 km ~ 10-2 mbar
     !   IF     ( logp.ge.2.0 .and. logp.le.3.0 ) THEN 
     !         kedd(l) = 2.e6 * 5.0**(logp-2.0)
     !! 1E7 at 500 km ~ 10-3 mbar
     !   ELSE IF     ( logp.ge.3.0 .and. logp.le.4.0 ) THEN 
     !         kedd(l) = 1.e7 * 3.0**(logp-3.0)
     !! 3E7 above 700 km ~ 10-4 mbar
     !   ELSEIF ( logp.gt.4.0                   ) THEN 
     !        kedd(l) = 3.e7
     !   ENDIF
     !ENDDO

     ! Kedd from (E7) in Vuitton 2019
     if (ngrid .eq. 1) then ! if 1D no dynamic mixing, we set the kedd in all column
       DO l=1,nlaykim_tot
         kedd(l) = 300.0 * ( 1.0E2 / press_c(l) )**1.5 * 3.0E7 /  &
                 ( 300.0 * ( 1.0E2 / press_c(l) )**1.5 + 3.0E7 ) 
       ENDDO
     else
       DO l=klev-4,nlaykim_tot
         ! JVO 18 : We keep the nominal profile in the GCM 5 upper layers 
         !          to have  a correct vertical mixing in the sponge layer
         kedd(l) = 300.0 * ( 1.0E2 / press_c(l) )**1.5 * 3.0E7 /  &
               ( 300.0 * ( 1.0E2 / press_c(l) )**1.5 + 3.0E7 ) 
       ENDDO
     endif
     
     if (ngrid .gt. 1) then ! not in 1D, no dynamic mixing
       ! Then adjust 10 layers profile fading to default value depending on kedd(ptop)
       DO l=klev-15,klev-5
          temp1   = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev-4)/press_c(klev-15)) )
          kedd(l) = 10.**( 3.0 + log10(kedd(klev-4)/1.e3) * temp1 )
       ENDDO
     endif
 
     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.1 : 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" !
     ! -------------------------------------------------------------------------------------------------------

#ifndef MESOSCALE
     IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.0)) ) .OR. (.NOT. moyzon_ch) ) THEN
#endif
        ! 1. Compute altitude for the grid point with hydrostat. equilib.
        ! ---------------------------------------------------------------

        ! a. For GCM layers we just copy-paste
        ! JVO 19 : Now physiq always sent correct altitudes with effective g for chemistry ( even if it's not the case in physiq )

        DO l=1,klev
           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)=( (czlay(ig,klev) + ( czlay(ig,klev) - czlev(ig,klev))) +rad )/1000. ! You shall not use czlev(zlev+1) whose value is 1.e7 !

        ! b. Extension in upper atmosphere with Vervack profile

        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 with effective g to have correct altitudes
           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

           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) + cpphis(ig) ) ) ) / 1000.0 ! z(phi) with g varying with altitude with reference to the geoid
        ENDDO

        DO l=klev+2,nlaykim_tot
          rinter(l) = 0.5*(rmil(l-1) + rmil(l)) ! should be balanced with the intermediate pressure rather than 0.5 
        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 ((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

        ! 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 ! nd_kim+1 is dissociation of N2 by GCR

                 krpddec   =   (   krpd(i,ialt  ,idec+1,klat)   * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec+1,klat)   * factalt       ) * (1.0-factlat) &
                             + (   krpd(i,ialt  ,idec+1,klat+1) * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec+1,klat+1) * factalt       ) * factlat

              if      ( factdec.lt.0. ) then 
                 krpddecm1 =   (   krpd(i,ialt  ,idec  ,klat)   * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec  ,klat)   * factalt       ) * (1.0-factlat) &
                             + (   krpd(i,ialt  ,idec  ,klat+1) * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec  ,klat+1) * factalt       ) * factlat
                 krate(l,i) = krpddecm1 * abs(factdec) + krpddec   * ( 1.0 + factdec)
              else if ( factdec.gt.0. ) then
                 krpddecp1 =   (   krpd(i,ialt  ,idec+2,klat)   * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec+2,klat)   * factalt       ) * (1.0-factlat) &
                             + (   krpd(i,ialt  ,idec+2,klat+1) * (1.0-factalt)                   &
                                 + krpd(i,ialt+1,idec+2,klat+1) * factalt       ) * factlat
                 krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)
              else if ( factdec.eq.0. ) then 
                 krate(l,i) = krpddec
              endif

           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,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

#ifndef MESOSCALE
     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
#endif

  ENDDO

END SUBROUTINE calchim
