      SUBROUTINE calchim(ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
     .                   ctemp,cplay,cplev,
     .                   dqyc)
      
c-------------------------------------------------
c  
c     Introduction d une routine chimique
c
c     Auteur: S. Lebonnois,  01/2000 | 09/2003 
c            adaptation pour Titan 3D: 02/2009
c
c-------------------------------------------------
c
      use dimphy
      implicit none
#include "dimensions.h"
#include "clesphys.h"
#include "paramet.h"
#include "YOMCST.h"

#include "titan_for.h"  
!!!  doit etre en accord avec titan.h
#include "aerprod.h"

c    Arguments
c    ---------

      INTEGER      ny                     ! nb de composes (nqmax-nmicro) 
      REAL         qy_c(jjm+1,klev,NC)    ! Especes chimiques apres adv.+diss.
      character*10 nomqy_c(NC+1)          ! Noms des especes chimiques
      REAL         declin_rad,ls_rad      ! declinaison et long solaire en radians
      REAL         dtchim                 ! pas de temps chimie
      REAL         ctemp(jjm+1,klev)      ! Temperature
      REAL         cplay(jjm+1,klev)      ! pression (Pa)
      REAL         cplev(jjm+1,klev)      ! pression intercouches (Pa)
      
      REAL         dqyc(jjm+1,klev,NC)    ! Tendances especes chimiques 
      
c    Local variables :
c    -----------------
c variables envoyees dans la chimie: double precision

      integer i,j,l,ic
      REAL  temp_c(klev),press_c(klev)     ! T,p(mbar) a 1 lat donnee
      REAL  declin_c                       ! declinaison en degres
      REAL  cqy(klev,NC)                   ! frac mol qui seront modifiees
      REAL  surfhaze(klev)
      REAL  cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4)
      REAL  rinter(klev),nb(klev)
      REAL  a,b
      character str1*1,str2*2
      integer ntotftop
      parameter (ntotftop=50)
      integer nftop(ntotftop),isaisonflux
      REAL  fluxtop(NC),latit,tabletmp(ntotftop),factflux
      character*10 nomsftop(ntotftop+1)
      character*20 formt1,formt2,formt3
      
c     variables locales initialisees au premier appel

      LOGICAL firstcal
      DATA    firstcal/.true./
      SAVE    firstcal
      
      REAL  mass(NC),duree
      REAL  tablefluxtop(NC,jjm+1,5)
      REAL  botCH4
      DATA  botCH4/0.05/
      REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
      integer reactif(5,NR),nom_prod(NC),nom_perte(NC)
      integer prod(200,NC),perte(2,200,NC)
      SAVE    mass,tablefluxtop,botCH4
      SAVE    reactif,nom_prod,nom_perte,prod,perte
      
c-----------------------------------------------------------------------
c***********************************************************************
c
c    Initialisations :
c    ----------------

      duree = dtchim  ! passage en real*4 pour appel a routines C

      IF (firstcal) THEN
            
        print*,'CHIMIE, premier appel'
        
c ************************************
c Au premier appel, initialisation de toutes les variables 
c pour les routines de la chimie.
c ************************************

        allocate(krpd(15,ND+1,klev,jjm+1),krate(klev,NR))

c Verification dimension verticale: coherence titan_for.h et klev
c --------------------------------

        if (klev.ne.NLEV) then
	   print*,'PROBLEME de coherence dans la dimension verticale:'
     .           ,klev,NLEV
           STOP
	endif

c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro
c ----------------------------------

        if (ny.ne.NC) then
	   print*,'PROBLEME de coherence dans le nombre de composes:'
     .           ,ny,NC
           STOP
	endif

c calcul de temp_c, densites et press_c a l'equateur:
c --------------------------------------------------

        print*,'pression, densites et temp a l equateur (chimie):'
        print*,'level, press_c, nb, temp_c'
        DO l=1,klev
c     temp_c (K):
         temp_c(l)  = ctemp(jjm/2+1,l)
c     press_c (mbar):
         press_c(l) = cplay(jjm/2+1,l)/100.
c     nb (cm-3):
         nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
         print*,l,press_c(l),nb(l),temp_c(l)
        ENDDO
        
c caracteristiques des composes:        
c -----------------------------
        mass = 0.0
        call comp(nomqy_c,mass)
        print*,'           Mass'
        do ic=1,NC
          print*,nomqy_c(ic),mass(ic)
        enddo
	

c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON... 

c flux dans la derniere couche:(issu du modele 1D, a 500 km) 
c -----------------------------
c       !!  en cm-2.s-1 !!
c  ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54 
c 
c Lecture des tables de flux en fonction lat et saison
c ----------------------------------------------------

	WRITE(str2(1:2),'(i2.2)') ntotftop
        formt1 = '(A10,'//str2//'(3X,A10))'
        formt2 = '(F12.6,'//str2//'(2X,F13.6))'
        formt3 = '(F13.6,'//str2//'(2X,F13.6))'
	
	do j=1,jjp1
	  do ic=1,NC
	    do l=1,5
	     tablefluxtop(ic,j,l) = 0.
	    enddo
	  enddo
	enddo
	
        do l=1,4
	  WRITE(str1(1:1),'(i1.1)') l
c hx1 -> Ls=RPI/2 / hx4 -> Ls=0
          open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1)
c         open(10,file="FLUXTOP/flux500.hs"//str1)
c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ...

          read(10,formt1) nomsftop
	  do j=1,ntotftop
           do i=1,10   !justification a gauche...
             if (nomsftop(j+1)(i:i).ne.' ') then
                nomsftop(j) = nomsftop(j+1)(i:)
                goto 100
             endif
           enddo
100        continue
c         print*,nomsftop(j)
           nftop(j) = 0
           do i=1,NC
	    if (nomqy_c(i).eq.nomsftop(j)) then
             nftop(j) = i
	    endif
	   enddo
c          if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j))
	  enddo

	 
c lecture des flux. Les formats sont un peu alambiques... 
c     a revoir eventuellement
	  do j=1,jjp1/2+1
	    read(10,formt2) latit,tabletmp
	    do i=1,ntotftop
              if (nftop(i).ne.0) then
	       tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
              endif
	    enddo
	  enddo
	  do j=jjp1/2+2,jjp1
	    read(10,formt3) latit,tabletmp
	    do i=1,ntotftop
              if (nftop(i).ne.0) then
	       tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
              endif
	    enddo
	  enddo

          close(10)
	enddo  ! l
	
	do j=1,jjp1
	  do ic=1,NC
	     tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5)
	  enddo
	enddo
	
c  Stockage des composes utilises dans la prod d'aerosols
c     (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
c     ! decalage de 1 car utilise dans le c !
c ------------------------------------------
c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN
c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10
c ------------------------------------------

        do ic=1,NC

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c         if (nomqy_c(ic).eq."CH4") then
c           do l=1,llm
c            do j=1,ip1jmp1
c              if (qy_c(j,l,ic).le.0.015) qy_c(j,l,ic) = 0.015 
c            enddo
c           enddo
c         endif
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          
	  if (nomqy_c(ic).eq."C4H2") then
            utilaer(10) = ic-1
	  endif
	  if (nomqy_c(ic).eq."HCN") then
            utilaer(6)  = ic-1
	    do j=1,jjp1
              do i=1,5
             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6.
              enddo
            enddo
	  endif
	  if (nomqy_c(ic).eq."HC3N") then
            utilaer(7)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."NCCN") then
            utilaer(14) = ic-1
	  endif
	  if (nomqy_c(ic).eq."CH3CN") then
            utilaer(15) = ic-1
            utilaer(16) = ic-1 ! si pas C2H3CN, CH3CN utilise, mais reac. nulle
	  endif
	  if (nomqy_c(ic).eq."H") then
            utilaer(1)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."H2") then
            utilaer(2)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."C2H2") then
            utilaer(3)  = ic-1
	    do j=1,jjp1
              do i=1,5
c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3
             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10.
              enddo
            enddo
	  endif
	  if (nomqy_c(ic).eq."AC6H6") then
            utilaer(13) = ic-1
	  endif
	  if (nomqy_c(ic).eq."C2H3CN") then
            utilaer(16) = ic-1
	  endif
	  if (nomqy_c(ic).eq."C2") then
            utilaer(4)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."C2H") then
            utilaer(5)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."C3N") then
            utilaer(8)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."H2CN") then
            utilaer(9)  = ic-1
	  endif
	  if (nomqy_c(ic).eq."C4H3") then
            utilaer(11) = ic-1
	  endif
	  if (nomqy_c(ic).eq."AC6H5") then
            utilaer(12) = ic-1
	  endif

	  if (nomqy_c(ic).eq."C2H6") then
	    do j=1,jjp1
              do i=1,5
c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640.
             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200.
              enddo
            enddo
	  endif
	  if (nomqy_c(ic).eq."C3H8") then
	    do j=1,jjp1
              do i=1,5
c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.)
              enddo
            enddo
	  endif
	  if (nomqy_c(ic).eq."C4H10") then
	    do j=1,jjp1
              do i=1,5
             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
              enddo
            enddo
	  endif

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c         if ((nomqy_c(ic).eq."HC3N").or.
c    $        (nomqy_c(ic).eq."C3N")) then
c     DO j=1,ip1jmp1 
c       do l=1,34  ! p>~ 1 mbar
c         qy_c(j,l,ic) = 1.e-30
c       enddo
c     ENDDO
c         endif
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        enddo
	        
c taux de photodissociations:
c --------------------------
        call disso(krpd,jjp1) 

c reactions chimiques:
c -------------------
        call chimie(nomqy_c,nb,temp_c,krate,reactif,
     .               nom_perte,nom_prod,perte,prod)
c        print*,'nom_prod, nom_perte:'
c        do ic=1,NC
c          print*,nom_prod(ic),nom_perte(ic)
c        enddo
c        print*,'premieres prod, perte(1:reaction,2:compagnon):'
c        do ic=1,NC
c          print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
c        enddo

c       l = klev-3
c       print*,'krate a p=',press_c(l),' reactifs et produits:'
c       do ic=1,ND+1
c         print*,ic,krpd(7,ic,l,4)*nb(l),"  ",
c    .     nomqy_c(reactif(1,ic)+1),
c    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
c    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
c       enddo
c       do ic=ND+2,NR
c         print*,ic,krate(l,ic),"  ",
c    .     nomqy_c(reactif(1,ic)+1),
c    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
c    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
c       enddo

      ENDIF  ! premier appel

c***********************************************************************
c-----------------------------------------------------------------------

c   calcul declin_c (en degres)
c   ---------------------------

      declin_c = declin_rad*180./RPI
c      print*,'declinaison en degre=',declin_c
       
c-----------------------------------------------------------------------
c
c   calcul du facteur d'interpolation entre deux saisons pour flux
c   --------------------------------------------------------------

       isaisonflux = int(ls_rad*2./RPI)+1
       if (isaisonflux.eq.5) then
         isaisonflux = 1
       endif
       factflux =    (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/
     .   (sin(isaisonflux*RPI/2.)-sin((isaisonflux-1)*RPI/2.))

c***********************************************************************
c***********************************************************************
c
c                BOUCLE SUR LES LATITUDES
c
      DO j=1,jjp1 
      
c***********************************************************************
c***********************************************************************

c-----------------------------------------------------------------------
c
c   Temperature, pression (mbar), densite (cm-3)
c   -------------------------------------------

       DO l=1,klev
c     temp_c (K):
               temp_c(l)  = ctemp(j,l)
c     press_c (mbar):
               press_c(l) = cplay(j,l)/100.
c     nb (cm-3):
               nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
       ENDDO
       
c-----------------------------------------------------------------------
c
c   Distances radiales (intercouches, en km)
c   ----------------------------------------

       rinter(1) = RA/1000.  
       do l=2,klev
c A REVOIR: g doit varier avec r !
         rinter(l) = rinter(l-1) + 
     .    (cplev(j,l-1)-cplev(j,l))/cplay(j,l)*(RD*temp_c(l)/RG)/1000.
c        print*,rinter(l)
       enddo

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c FLUX AU TOP: interpolation en saison et
c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C) 
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        do ic=1,NC
          fluxtop(ic) = tablefluxtop(ic,j,isaisonflux)  *(1-factflux)
     .                + tablefluxtop(ic,j,isaisonflux+1)*  factflux
          fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5)
        enddo

c TEST: sortie de l'un des flux avec saveLs et factflux
c       dans un fichier special pour tracer avec gnuplot
c       if (j.eq.2) then
c       open(unit=11,file="flux_80N.txt",status='old',position='append')
c         write(11,*) ls_rad*180./RPI, factflux, 
c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
c         write(11,*) ls_rad*180./RPI, factflux, 
c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
c       close(11)
c       endif
c       if (j.eq.17) then
c       open(unit=11,file="flux_eq.txt",status='old',position='append')
c         write(11,*) ls_rad*180./RPI, factflux, 
c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
c         write(11,*) ls_rad*180./RPI, factflux, 
c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
c       close(11)
c       endif
c FIN TEST: sortie de l'un des flux avec ls et factflux

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       if (firstcal.and.(j.eq.1)) then
         print*,'Alt, densites et temp au pole (chimie):'
         print*,'level, z_bas, nb, temp_c'
         do l=1,klev
          print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l)
         enddo
       endif

       if (firstcal.and.(j.eq.jjm/2)) then
c         print*,'g,mugaz'
c         print*,g,mugaz
         print*,'Alt, densites et temp a l equateur (chimie):'
         print*,'level, z_bas, nb, temp_c'
         do l=1,klev
          print*,l,rinter(l)-RA/1000.,nb(l),temp_c(l)
         enddo
       endif
       
c-----------------------------------------------------------------------
c
c   composition 
c   ------------
       
       do ic=1,NC
        do l=1,klev
          cqy(l,ic) = qy_c(j,l,ic) 
        enddo
       enddo
              
c-----------------------------------------------------------------------
c
c   total haze area (um2/cm3)
c   -------------------------

       if (htoh2.eq.1) then
        do l=1,klev
! ATTENTION, INVERSION PAR RAPPORT A pg2.F !!!
         surfhaze(l) = psurfhaze(j,klev+1-l) 
c        if (j.eq.25)
c    .    print*,'psurfhaze en um2/cm3:',surfhaze(l)
        enddo
       endif

c-----------------------------------------------------------------------
c
c   Appel de chimietitan
c   --------------------
       
       call gptitan(jjp1,rinter,temp_c,nb,
     $              nomqy_c,cqy,fluxtop,
     $              declin_c,duree,(j-1),mass,
     $              botCH4,krpd,krate,reactif,
     $              nom_prod,nom_perte,prod,perte,
     $              aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
     $              htoh2,surfhaze)
       
c       if ( j.eq.jjm/2 ) 
c    $        print*,cqy(1,1),cqy(klev,1),cqy(1,2),cqy(klev,2)
c       if ( j.eq.jjm/2 ) 
c    $  print*,qy_c(j,1,1),qy_c(j,klev,1),qy_c(j,1,2),qy_c(j,klev,2)

c       stop

c   Tendances composition 
c   ---------------------
     
       do ic=1,NC
        do l=1,klev
          dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim  ! en /s
        enddo
       enddo

c-----------------------------------------------------------------------
c
c   production aer 
c   --------------
       
       if (aerprod.eq.1) then

       do ic=1,4
        do l=1,klev
          prodaer(j,l,ic) = cprodaer(l,ic) 
          maer(j,l,ic)    = cmaer(l,ic) 
          csn(j,l,ic)     = ccsn(l,ic) 
          csh(j,l,ic)     = ccsh(l,ic) 
        enddo
       enddo

       endif

c***********************************************************************
c***********************************************************************
c
c              FIN: BOUCLE SUR LES LATITUDES 
c
      ENDDO
     
c***********************************************************************
c***********************************************************************


      firstcal = .false.
      RETURN
      END
