Index: trunk/LMDZ.TITAN/libf/phytitan/calchim.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calchim.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/calchim.F	(revision 1126)
@@ -1,4 +1,4 @@
       SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
-     .                   ctemp,cplay,cplev,
+     .                   ctemp,cplay,cplev,czlay,czlev,
      .                   dqyc)
       
@@ -10,12 +10,13 @@
 c            adaptation pour Titan 3D: 02/2009
 c            adaptation pour // : 04/2013
-c
+c            extension chimie jusqu a 1300km : 10/2013
+c 
 c-------------------------------------------------
 c
       use dimphy
       use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
-     .                     NLEV,NC,ND,NR
+     .                     NLEV,NLRT,NC,ND,NR
       USE comgeomphy,  only: rlatd
-      use moyzon_mod,  only: klat
+      USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
       implicit none
 #include "dimensions.h"
@@ -35,5 +36,7 @@
       REAL         ctemp(nlon,klev)      ! Temperature
       REAL         cplay(nlon,klev)      ! pression (Pa)
-      REAL         cplev(nlon,klev)      ! pression intercouches (Pa)
+      REAL         cplev(nlon,klev+1)    ! pression intercouches (Pa)
+      REAL         czlay(nlon,klev)      ! altitude (m)
+      REAL         czlev(nlon,klev+1)    ! altitude intercouches (m)
       
       REAL         dqyc(nlon,klev,NC)    ! Tendances especes chimiques 
@@ -46,17 +49,19 @@
 c variables envoyees dans la chimie: double precision
 
-      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  temp_c(NLEV)
+      REAL  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
+      REAL  cqy(NLEV,NC)    ! frac mol qui seront modifiees
+      REAL  cqy0(NLEV,NC)    ! frac mol avant modif
+      REAL  surfhaze(NLEV)
+      REAL  cprodaer(NLEV,4),cmaer(NLEV,4)
+      REAL  ccsn(NLEV,4),ccsh(NLEV,4)
+! rmil: milieu de couche, grille pour K,D,p,ct,T,y
+! rinter: intercouche (grille RA dans la chimie)
+      REAL  rmil(NLEV),rinter(NLEV),nb(NLEV)
+      REAL,save :: kedd(NLEV)
+
       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)
+      REAL  latit
       character*20 formt1,formt2,formt3
       
@@ -67,14 +72,20 @@
       SAVE    firstcal
       
-      REAL  mass(NC),duree
-      REAL  tablefluxtop(NC,jjp1,5)
-      REAL  botCH4
+      integer dec,declinint,ialt
+      REAL  declin_c                       ! declinaison en degres
+      real  factalt,factdec,krpddec,krpddecp1,krpddecm1
+      real  duree
+      REAL,save :: mass(NC)
+      REAL,save,allocatable :: md(:,:)
+      REAL,save :: botCH4
       DATA  botCH4/0.05/
+      real,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
       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
-      
+      integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
+      integer,save :: prod(200,NC),perte(2,200,NC)
+      character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
+      real  ficalt,oldq,newq,zalt
+      logical okfic
+
 c-----------------------------------------------------------------------
 c***********************************************************************
@@ -94,16 +105,7 @@
 c ************************************
 
-        allocate(krpd(15,ND+1,klev,jjp1),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
+        allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
+
+c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
 c ----------------------------------
 
@@ -114,114 +116,72 @@
 	endif
 
-c calcul de temp_c, densites et press_c au milieu de l'ensemble des points:
-c ----------------------------------------------------------------------
-
-        print*,'pression, densites et temp (chimie):'
+c calcul de temp_c, densites et press_c en moyenne planetaire :
+c --------------------------------------------------------------
+
+        print*,'pression, densites et temp (init chimie):'
         print*,'level, press_c, nb, temp_c'
         DO l=1,klev
+         rinter(l)  = (zlevmoy(l)+RA)/1000.
+         rmil(l)    = (zlaymoy(l)+RA)/1000.
 c     temp_c (K):
-         temp_c(l)  = ctemp(nlon/2+1,l)
+         temp_c(l)  = tmoy(l)
 c     press_c (mbar):
-         press_c(l) = cplay(nlon/2+1,l)/100.
+         press_c(l) = playmoy(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)
+         print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
         ENDDO
+
+c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+1,NLEV
+         rinter(l) = rinter(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+         rmil(l)   = rmil(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+       enddo        
+
+c lecture de tcp.ver, une seule fois
+c remplissage r1d,t1d,ct1d,p1d
+        open(11,file='../../INPUT/tcp.ver',status='old')
+        read(11,*) dummy
+        do i=1,131
+          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
+c         print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
+        enddo 
+        close(11)
+
+c extension pour klev+1 a NLEV avec tcp.ver
+c la jonction klev/klev+1 est brutale... Tant pis ?
+        ialt=1
+        do l=klev+1,NLEV
+           zalt=rmil(l)-RA/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
         
 c caracteristiques des composes:        
 c -----------------------------
-        mass = 0.0
-        call comp(nomqy_c,mass)
+        mass(:) = 0.0
+        call comp(nomqy_c,nb,temp_c,mass,md)
         print*,'           Mass'
         do ic=1,NC
           print*,nomqy_c(ic),mass(ic)
+c         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,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  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
@@ -244,9 +204,4 @@
 	  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
@@ -268,10 +223,4 @@
 	  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
@@ -298,28 +247,4 @@
 	  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
 
@@ -369,4 +294,29 @@
 c       enddo
 
+
+c   init kedd
+c   ---------
+c   kedd stays constant with time and space 
+c   => linked to pressure rather than altitude...
+ 
+      kedd(:) = 5e2 ! valeur mise par defaut  
+               ! UNITE ? doit etre ok pour gptitan
+      do l=1,NLEV
+       zalt=rmil(l)-RA/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
 
@@ -380,16 +330,4 @@
 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***********************************************************************
@@ -410,4 +348,23 @@
 c***********************************************************************
 
+c-----------------------------------------------------------------------
+c
+c   Distance radiale (intercouches, en km)
+c   ----------------------------------------
+
+       do l=1,klev
+         rinter(l) = (RA+czlev(j,l))/1000.
+         rmil(l)   = (RA+czlay(j,l))/1000.
+c        print*,rinter(l)
+       enddo
+
+c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
+       do l=klev+1,NLEV
+         rinter(l) = rinter(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+         rmil(l)   = rmil(klev)
+     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
+       enddo
+        
 c-----------------------------------------------------------------------
 c
@@ -423,51 +380,75 @@
                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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+c extension pour klev+1 a NLEV avec tcp.ver
+       ialt=1
+       do l=klev+1,NLEV
+           zalt=rmil(l)-RA/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
+               
+c-----------------------------------------------------------------------
+c
+c   passage krpd => krate 
+c   ---------------------
+c initialisation krate pour dissociations 
+
+      if ((declin_c*10+267).lt.14.) then
+          declinint = 0
+          dec       = 0
+      else 
+       if ((declin_c*10+267).gt.520.) then
+          declinint = 14
+          dec       = 534
+       else 
+          declinint = 1
+          dec       = 27
+          do while( (declin_c*10+267).ge.real(dec+20) )
+            dec       = dec+40
+            declinint = declinint+1
+          enddo
+       endif
+      endif
+      if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
+          factdec = ( declin_c - (dec-267)/10. ) / 4.
+      else
+          factdec = ( declin_c - (dec-267)/10. ) / 2.7
+      endif
+
+      do l=1,NLEV
+
+c INTERPOL EN ALT POUR k (krpd tous les deux km)
+        ialt    = int((rmil(l)-RA/1000.)/2.)+1
+        factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
+
+        do i=1,ND+1
+          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
+          krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
+          krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt)
+     &              + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
+
+                    ! ND+1 correspond a la dissociation de N2 par les GCR
+          if ( factdec.lt.0. ) then 
+        krate(l,i) = krpddecm1 * abs(factdec)
+     &             + krpddec   * ( 1 + factdec)
+          endif
+          if ( factdec.gt.0. ) then
+        krate(l,i) = krpddecp1 * factdec
+     &             + krpddec   * ( 1 - factdec)
+          endif
+          if ( factdec.eq.0. ) krate(l,i) = krpddec
+        enddo        
+      enddo        
 
 c-----------------------------------------------------------------------
@@ -475,11 +456,52 @@
 c   composition 
 c   ------------
-       
+
        do ic=1,NC
         do l=1,klev
-          cqy(l,ic) = qy_c(j,l,ic) 
-        enddo
-       enddo
-              
+          cqy(l,ic)  = qy_c(j,l,ic) 
+          cqy0(l,ic) = cqy(l,ic)
+        enddo
+       enddo
+
+c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
+
+      WRITE(str2,'(i2.2)') klat(j)
+      fich = "comp_"//str2//".dat"
+      inquire (file=fich,exist=okfic)
+      if (okfic) then
+       open(11,file=fich,status='old')
+c premiere ligne=declin
+       read(11,'(A15)') ficdec
+       write(curdec,'(E15.5)') declin_c
+c si la declin est la meme: ce fichier a deja ete reecrit 
+c on lit la colonne t-1 au lieu de la colonne t
+c (cas d une bande de latitude partagee par 2 procs)
+       do ic=1,NC
+        read(11,'(A10)') name
+        if (name.ne.nomqy_c(ic)) then
+          print*,"probleme lecture ",fich
+          print*,name,nomqy_c(ic)
+          stop
+        endif
+        if (ficdec.eq.curdec) then
+          do l=klev+1,NLEV
+            read(11,*) ficalt,cqy(l,ic),newq
+          enddo
+        else
+          do l=klev+1,NLEV
+            read(11,*) ficalt,oldq,cqy(l,ic)
+          enddo
+        endif
+       enddo
+       close(11)
+      else       ! le fichier n'est pas la
+       do ic=1,NC
+        do l=klev+1,NLEV
+          cqy(l,ic)=cqy(klev,ic)
+        enddo
+       enddo
+      endif
+      cqy0 = cqy
+
 c-----------------------------------------------------------------------
 c
@@ -502,7 +524,7 @@
        
        call gptitan(rinter,temp_c,nb,
-     $              nomqy_c,cqy,fluxtop,
-     $              declin_c,duree,(klat(j)-1),mass,
-     $              botCH4,krpd,krate,reactif,
+     $              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,
@@ -514,5 +536,5 @@
        do ic=1,NC
         do l=1,klev
-          dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim  ! en /s
+          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
         enddo
        enddo
@@ -535,4 +557,23 @@
 
        endif
+
+c-----------------------------------------------------------------------
+c
+c   sauvegarde compo sur NLEV 
+c   -------------------------
+
+c dans fichier compo_klat(j) (01 à 49)
+       
+      open(11,file=fich,status='replace')
+c premiere ligne=declin
+      write(11,'(E15.5)') declin_c
+      do ic=1,NC
+       write(11,'(A10)') nomqy_c(ic)
+       do l=klev+1,NLEV
+        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
+     .                                cqy0(l,ic),cqy(l,ic)
+       enddo
+      enddo
+      close(11)
 
 c***********************************************************************
Index: trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90	(revision 1126)
+++ trunk/LMDZ.TITAN/libf/phytitan/cirs_haze.F90	(revision 1126)
@@ -0,0 +1,38 @@
+subroutine cirs_haze(press,wno,taeros,taeroscat)
+IMPLICIT NONE
+
+real,intent(in)   :: press,wno
+real,intent(inout):: taeros,taeroscat
+
+!---------------------------
+! Attention: taeros est en m-1
+! mais la valeur semble devoir etre multipliee par lambda (en m) pour
+! obtenir une valeur comparable a celle fournies par Sandrine...
+! a tirer au clair...
+!---------------------------
+
+real         :: taerosold
+logical,save :: firstcall=.true.
+
+if (firstcall) then
+   print*,"CIRS HAZE"
+   firstcall=.false.
+endif
+
+if (wno.eq.600.) then
+ print*,press,wno,taeros,taeroscat
+endif
+
+taerosold = taeros
+
+! modif de taeros
+
+! taeroscat est modifie proportionnellement a taeros
+
+if (taerosold.ne.0.) then
+   taeroscat = taeroscat/taerosold*taeros
+endif
+
+! Je maintiens le cbar du calcul microphysique
+
+end subroutine cirs_haze
Index: trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90	(revision 1126)
@@ -34,7 +34,9 @@
 
 ! ancien titan_for.h
-      INTEGER, PARAMETER :: NLEV=llm,NC=44,ND=54,NR=377
+      INTEGER, PARAMETER :: NLEV=llm+70,NC=44,ND=54,NR=377
 !$OMP THREADPRIVATE(NLEV,NC,ND,NR)
 !!!  doivent etre en accord avec titan.h
+! pour l'UV (650 niveaux de 2 km)
+      INTEGER, PARAMETER :: NLRT=650
 
 ! ancien diagmuphy.h
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F	(revision 1126)
@@ -162,4 +162,7 @@
           XICLDI2(K)=TNI
  420  CONTINUE
+
+c DEBUG
+c       print*,"wnoi=",WNOI
 
 C
@@ -284,4 +287,20 @@
         TauGID(ig,:,:) = TAUGID_1pt(:,:) 
 
+c DEBUG
+c     if(ig.eq.(ngrid/2+16)) then
+c         print*,ig,'/',KLON,':'
+c         print*,'TauHID_1',TAUHID(ig,1,:)
+c         print*,'TauGID_1',TAUGID(ig,1,:)
+c         print*,'TauHID_50',TAUHID(ig,50,:)
+c         print*,'TauGID_50',TAUGID(ig,50,:)
+c         print*,"DTAUI_1=",DTAUI(ig,1,:)
+c         print*,"DTAUI_50=",DTAUI(ig,50,:)
+c         print*,'cosBI_1',COSBI(ig,1,:)
+c         print*,'cosBI_50',COSBI(ig,50,:)
+c         print*,'WBARI_1',WBARI(ig,1,:)
+c         print*,'WBARI_50',WBARI(ig,50,:)
+c         stop
+c     endif
+
 c************************************************************************
 c************************************************************************
Index: trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F	(revision 1126)
@@ -282,4 +282,6 @@
 c     CBAR=xv3(j,k)
 
+c     print*, 'HERE, CIRS AEROSOLS'
+c     call cirs_haze(PRESS(J),WNOI(K),TAEROS,TAEROSCAT,CBAR)
 
 c*********** EN TRAVAUX ***************************
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F	(revision 1126)
@@ -148,4 +148,7 @@
 c       close(1)
 
+c DEBUG
+c       print*,"wnov=",WNOV
+
       endif    ! fin initialisations premier appel
 
@@ -167,5 +170,5 @@
         else
          if (ig.eq.1)  then
-c initialisation zqaer_1pt a partir d'une look-up table (uniforme en ig)
+c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
 c boucle sur nrad=10
            open(10,file="qaer_eq_1d.dat")
@@ -220,4 +223,14 @@
         TauGVD(ig,:,:) = TAUGVD_1pt(:,:) 
 
+c DEBUG
+c     if(ig.eq.(ngrid/2+16)) then
+c         print*,ig,'/',KLON,':'
+c         print*,'TauHVD_1',TAUHVD(ig,1,:)
+c         print*,'TauGVD_1',TAUGVD(ig,1,:)
+c         print*,'TauHVD_50',TAUHVD(ig,50,:)
+c         print*,'TauGVD_50',TAUGVD(ig,50,:)
+c     stop
+c     endif
+
  101  CONTINUE
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F	(revision 1126)
@@ -231,4 +231,7 @@
 
        IF (TAEROSCAT.le.0.) CBAR=0. 
+
+c     print*, 'HERE, CIRS AEROSOLS'
+c     call cirs_haze(PRESS(J),WNOV(K),TAEROS,TAEROSCAT,CBAR)
 
 1699  FORMAT(a3,2I3,3(ES15.7,1X))
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq.F	(revision 1126)
@@ -419,5 +419,5 @@
          itap    = 0
          itaprad = 0
-         itapchim    = 0
+         itapchim    = 1
 
 c init rnuabar
@@ -1126,5 +1126,9 @@
          tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
      .                        + d_tr_mph(:,:,1:nmicro)*dtime
-
+c        call WriteField_phy('physiq_d_tr_mph01',
+c    .                          d_tr_mph(:,:,1),klev)
+c        call WriteField_phy('physiq_d_tr_mph10',
+c    .                          d_tr_mph(:,:,10),klev)
+	endif
 c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
 c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
@@ -1141,6 +1145,4 @@
         ENDDO
 
-	endif
-
 c condensation:
 c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
@@ -1149,10 +1151,11 @@
      .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
         endif
+
+c chimie:
         if ((chimi).and.(nqmax.gt.nmicro)) then
-c chimie:
          tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
 	endif
 
-      endif
+      endif  !iflag_trac=1
 
 c       ch4=0.
Index: trunk/LMDZ.TITAN/libf/phytitan/phytrac.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phytrac.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/phytrac.F	(revision 1126)
@@ -459,5 +459,5 @@
 c ------------
        CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
-     .              ztemp,zplay,zplev,
+     .              ztemp,zplay,zplev,zzlay,zzlev,
      .              pdyfi)   
 c ychim ne doit pas etre modifie, pdyfi en /s
Index: trunk/LMDZ.TITAN/libf/phytitan/profile.F
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 1121)
+++ trunk/LMDZ.TITAN/libf/phytitan/profile.F	(revision 1126)
@@ -94,7 +94,8 @@
        b2 =       3183. 
        c2 =       4737. 
+c pres est en Pa => conversion car expression veut p en hPa
        DO il=1,nlev
-         temp(il)=a1*exp(-((pres(il)-b1)/c1)**2.) 
-     .          + a2*exp(-((pres(il)-b2)/c2)**2.)
+         temp(il)=a1*exp(-((pres(il)/100.-b1)/c1)**2.) 
+     .          + a2*exp(-((pres(il)/100.-b2)/c2)**2.)
        ENDDO
        zkm(1)  = 0.0
