Index: LMDZ.3.3/trunk/libf/dyn3d/calfis.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/calfis.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/dyn3d/calfis.F	(revision 344)
@@ -1,6 +1,31 @@
-      SUBROUTINE calfis(nq, lafin, rdayvrai,rday_ecri, heure,
-     $            pucov,pvcov,pteta,pq,pmasse,pps,pp,ppk,pphis,pphi,
-     $            pducov,pdvcov,pdteta,pdq,pw, clesphy0,
-     $            pdufi,pdvfi,pdhfi,pdqfi,pdpsfi                   )
+      SUBROUTINE calfis(nq, 
+     $                  lafin,
+     $                  rdayvrai,
+     $                  rday_ecri, 
+     $                  heure,
+     $                  pucov,
+     $                  pvcov,
+     $                  pteta,
+     $                  pq,
+     $                  pmasse,
+     $                  pps,
+     $                  pp,
+     $                  ppk,
+     $                  pphis,
+     $                  pphi,
+     $                  pducov,
+     $                  pdvcov,
+     $                  pdteta,
+     $                  pdq,
+     $                  pw, 
+#ifdef INCA_CH4
+     $                  flxw,   
+#endif
+     $                  clesphy0,
+     $                  pdufi,
+     $                  pdvfi,
+     $                  pdhfi,
+     $                  pdqfi,
+     $                  pdpsfi)
 c
 c    Auteur :  P. Le Van, F. Hourdin 
@@ -132,5 +157,9 @@
       REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
       REAL unskap, pksurcp
-c
+ 
+#ifdef INCA_CH4
+      REAL flxw(iip1,jjp1,llm)
+      REAL flxwfi(ngridmx,llm)
+#endif
       
       EXTERNAL gr_dyn_fi,gr_fi_dyn
@@ -421,4 +450,8 @@
       ENDDO
 
+#ifdef INCA_CH4
+      CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
+#endif
+
 c-----------------------------------------------------------------------
 c   Appel de la physique:
@@ -426,11 +459,33 @@
 
 
-      CALL physiq (ngridmx,llm,nq,debut,lafin,
-     ,   rdayvrai,rday_ecri,heure,dtphys,zplev,zplay,zphi,zphis,airefi,
-     ,             presnivs,clesphy0,  zufi, zvfi,ztfi, zqfi,  
-ccc     ,             pcvgu, pcvgv, pcvgt, pcvgq, 
-     ,             pvervel,
-C - sorties
-     s             zdufi, zdvfi, zdtfi, zdqfi,zdpsrf              )
+      CALL physiq (ngridmx,
+     .             llm,
+     .             nq,
+     .             debut,
+     .             lafin,
+     .             rdayvrai,
+     .             rday_ecri,
+     .             heure,
+     .             dtphys,
+     .             zplev,
+     .             zplay,
+     .             zphi,
+     .             zphis,
+     .             airefi,
+     .             presnivs,
+     .             clesphy0, 
+     .             zufi,
+     .             zvfi,
+     .             ztfi,
+     .             zqfi,  
+     .             pvervel,
+#ifdef INCA_CH4
+     .             flxwfi,
+#endif
+     .             zdufi, 
+     .             zdvfi, 
+     .             zdtfi,
+     .             zdqfi,
+     .             zdpsrf)
 
 500   CONTINUE
Index: LMDZ.3.3/trunk/libf/dyn3d/gcm.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/gcm.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/dyn3d/gcm.F	(revision 344)
@@ -2,4 +2,7 @@
       PROGRAM gcm
 
+#ifdef INCA
+      USE transport_controls, ONLY : adv_flg, mmt_adj
+#endif  
       USE IOIPSL
 
@@ -127,4 +130,8 @@
       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
 
+#ifdef INCA_CH4
+      REAL :: flxw(ip1jmp1,llm)
+#endif
+
       LOGICAL offline  ! Controle du stockage ds "fluxmass"
       PARAMETER (offline=.true.)
@@ -409,6 +416,23 @@
 cccc      de masse pour  Van-Leer dans la routine  tracvl  .
 c
-            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
-     *                      p, masse, dq,  iadv(1), teta, pk     )
+      CALL vanleer(numvanle,
+     .             iapp_tracvl,
+     .             nqmx,
+     .             q,
+     .             pbaru,
+     .             pbarv,
+     .             p, 
+     .             masse, 
+     .             dq,  
+     .             iadv(1), 
+     .             teta, 
+#ifdef INCA_CH4
+     .             flxw, 
+     .             pk, 
+     .             mmt_adj, 
+     .             adv_flg)
+#else
+     .             pk) 
+#endif
 c
 c                   ...  Modif  F.Codron  ....
@@ -470,7 +494,33 @@
            ENDIF
 c
-        CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
-     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
-     $     du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi  )
+        CALL calfis( nqmx, 
+     $               lafin,
+     $               rdayvrai,
+     $               rday_ecri,
+     $               time,
+     $               ucov,
+     $               vcov,
+     $               teta,
+     $               q,
+     $               masse,
+     $               ps,
+     $               p,
+     $               pk,
+     $               phis,
+     $               phi,
+     $               du,
+     $               dv,
+     $               dteta,
+     $               dq,
+     $               w,
+#ifdef INCA_CH4
+     $               flxw,
+#endif
+     $               clesphy0,
+     $               dufi,
+     $               dvfi,
+     $               dhfi,
+     $               dqfi,
+     $               dpfi)
 
 c      ajout des tendances physiques:
Index: LMDZ.3.3/trunk/libf/dyn3d/tracvl.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/tracvl.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/dyn3d/tracvl.F	(revision 344)
@@ -1,4 +1,21 @@
-      SUBROUTINE tracvl(numvanle,iapp_tracvl,nq,pbaru,pbarv ,
-     *                    p, masse , q, iapptrac, iadv1, teta, pk  )
+      SUBROUTINE tracvl(numvanle,
+     *                  iapp_tracvl,
+     *                  nq,
+     *                  pbaru,
+     *                  pbarv,
+     *                  p, 
+     *                  masse, 
+     *                  q, 
+     *                  iapptrac, 
+     *                  iadv1, 
+     *                  teta, 
+#ifdef INCA_CH4
+     *                  flxw,
+     *                  pk, 
+     *                  mmt_adj,
+     *                  adv_flg)
+#else
+     *                  pk) 
+#endif
 c
 c     Auteur :  F. Hourdin
@@ -20,10 +37,16 @@
 c
       INTEGER numvanle, nq, iapp_tracvl, iapptrac, iadv1
-
       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
       REAL q(ip1jmp1,llm,nq),masse(ip1jmp1,llm)
       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
       REAL pk(ip1jmp1,llm)
-
+#ifdef INCA_CH4
+      INTEGER, PARAMETER :: ntra   = 1
+      INTEGER, PARAMETER :: nprath = 1
+      INTEGER            :: adv_flg(nq)
+      REAL               :: mmt_adj(ip1jmp1,llm,nprath)
+      REAL, SAVE         :: qpente(ip1jmp1,llm,10,nprath)
+      REAL               :: flxw(ip1jmp1,llm)
+#endif
 c     ....  var. locales  .....
 c
@@ -82,4 +105,8 @@
         CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
 
+#ifdef INCA_CH4
+      ! ... Flux de masse diaganostiques traceurs
+      flxw = wg / FLOAT(iapp_tracvl) 
+#endif
 
 c  test sur l'eventuelle creation de valeurs negatives de la masse
@@ -115,9 +142,29 @@
          ENDIF
 
-         DO iq = numvan, nq
-          CALL vlsplt( q(1,1,iq), 2. ,massem,wg,pbarug,pbarvg,dtvr )
-         ENDDO
+#ifdef INCA_CH4
+      do iq = 2, 10
+        qpente(:,:,iq,1)=qpente(:,:,iq,1)*mmt_adj(:,:,1)
+      enddo
+#endif
 
-         iadvtr=0
+      DO iq = numvan, 2 
+#ifdef INCA
+      IF (adv_flg(iq) == 0) CYCLE
+#endif
+      CALL vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+      ENDDO
+
+#ifdef INCA_CH4
+!     CALL prather(q(1,1,3),wg,massem,pbarug,pbarvg,ntra,qpente(1,1,1,1))
+#endif
+
+      DO  iq =3,nq
+#ifdef INCA
+      IF (adv_flg(iq) == 0) CYCLE
+#endif
+      CALL vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
+      ENDDO
+
+      iadvtr=0
 
 c   on reinitialise a zero les flux de masse cumules.
Index: LMDZ.3.3/trunk/libf/dyn3d/vanleer.F
===================================================================
--- LMDZ.3.3/trunk/libf/dyn3d/vanleer.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/dyn3d/vanleer.F	(revision 344)
@@ -1,4 +1,21 @@
-      SUBROUTINE vanleer(numvanle,iapp_tracvl,nq,q,pbaru,pbarv ,
-     *                     p ,masse, dq ,  iadv1, teta, pk      )
+      SUBROUTINE vanleer(numvanle,
+     *                   iapp_tracvl,
+     *                   nq,
+     *                   q,
+     *                   pbaru,
+     *                   pbarv,
+     *                   p,
+     *                   masse, 
+     *                   dq, 
+     *                   iadv1,
+     *                   teta,
+#ifdef INCA_CH4
+     *                   flxw, 
+     *                   pk,
+     *                   mmt_adj, 
+     *                   adv_flg)
+#else
+     *                   pk)
+#endif
 c
       IMPLICIT NONE
@@ -24,4 +41,10 @@
       REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nq),dq( ip1jmp1,llm,nq )
       REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
+#ifdef INCA_CH4
+      INTEGER, PARAMETER :: nprath=1
+      INTEGER            :: adv_flg(nq)
+      REAL               :: mmt_adj(iip1,jjp1,llm,nprath)
+      REAL               :: flxw(ip1jmp1,llm)
+#endif
 c  ..................................................................
 c
@@ -64,6 +87,23 @@
 c   advection
 
-      CALL tracvl( numvanle,iapp_tracvl,nq,pbaru,pbarv,p , masse,q  ,
-     *                      iapptrac, iadv1, teta ,pk              )
+      CALL tracvl( numvanle,
+     .             iapp_tracvl,
+     .             nq,
+     .             pbaru,
+     .             pbarv,
+     .             p, 
+     .             masse,
+     .             q,
+     .             iapptrac, 
+     .             iadv1, 
+     .             teta,
+#ifdef INCA_CH4
+     .             flxw,
+     .             pk,
+     .             mmt_adj,
+     .             adv_flg)
+#else
+     .             pk)
+#endif
 
       IF( numvanle.EQ.1 ) THEN
Index: LMDZ.3.3/trunk/libf/phylmd/physiq.F
===================================================================
--- LMDZ.3.3/trunk/libf/phylmd/physiq.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/phylmd/physiq.F	(revision 344)
@@ -1,11 +1,38 @@
 c $Header$
 c
-      SUBROUTINE physiq (nlon,nlev,nqmax  ,
-     .            debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
-     .            paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
-     .            u,v,t,qx,
-     .            omega,
-     .            d_u, d_v, d_t, d_qx, d_ps)
+      SUBROUTINE physiq (nlon,
+     .                   nlev,
+     .                   nqmax,
+     .                   debut,
+     .                   lafin,
+     .                   rjourvrai,
+     .                   rjour_ecri,
+     .                   gmtime,
+     .                   pdtphys,
+     .                   paprs,
+     .                   pplay,
+     .                   pphi,
+     .                   pphis,
+     .                   paire,
+     .                   presnivs,
+     .                   clesphy0,
+     .                   u,
+     .                   v,
+     .                   t,
+     .                   qx,
+     .                   omega,
+#ifdef INCA_CH4
+     .                   flxmass_w,
+#endif
+     .                   d_u, 
+     .                   d_v, 
+     .                   d_t, 
+     .                   d_qx, 
+     .                   d_ps)
       USE ioipsl
+#ifdef INCA
+      USE chemshut 
+      USE obs_pos
+#endif
       IMPLICIT none
 c======================================================================
@@ -71,5 +98,9 @@
       INTEGER npas, nexca, itimestep
       logical rnpb
+#ifdef INCA
+      parameter(rnpb=.false.)
+#else
       parameter(rnpb=.true.)
+#endif
       PARAMETER (npas=1440)
       PARAMETER (nexca=48)
@@ -174,4 +205,6 @@
       REAL omega(klon,klev)
 
+      REAL flxmass_w(klon,klev)
+
       REAL d_u(klon,klev)
       REAL d_v(klon,klev)
@@ -326,4 +359,9 @@
       REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
       REAL frac_nucl(klon,klev) ! idem (nucleation)
+
+#ifdef INCA
+      REAL          :: calday
+#endif
+
 cAA
       REAL rain_fall(klon) ! pluie
@@ -403,4 +441,5 @@
 c Variables locales
 c
+      REAL rhcl(klon,klev)    ! humiditi relative ciel clair
       REAL dialiq(klon,klev)  ! eau liquide nuageuse
       REAL diafra(klon,klev)  ! fraction nuageuse
@@ -1401,4 +1440,30 @@
          ENDIF
 c
+
+#ifdef INCA
+           calday = FLOAT(MOD(NINT(xjour),360)) + gmtime
+           print *, 'initial time ', xjour, calday
+#ifdef INCAINFO
+           PRINT *, 'Appel CHEMINI ...'
+#endif
+           CALL chemini( rpi,
+     $                   rg,
+     $                   ra,
+     $                   paire,
+     $                   rlat,
+     $                   rlon,
+     $                   calday,
+     $                   tracnam,
+     $                   natsnam,
+     $                   mxoutflds,
+     $                   outinst,
+     $                   outtimav,
+     $                   klon,
+     $                   nqmax)
+#ifdef INCAINFO
+           PRINT *, 'OK.'
+#endif
+#endif
+
       ENDIF
 c
@@ -1848,5 +1913,5 @@
      .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
      .           frac_impa, frac_nucl,
-     .           prfl, psfl)
+     .           prfl, psfl, rhcl)
       DO k = 1, klev
       DO i = 1, klon
@@ -1929,4 +1994,47 @@
       ENDDO
       ENDDO
+
+#ifdef INCA
+           calday = FLOAT(julien) + gmtime
+
+#ifdef INCA_AER
+      call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs,
+     &   prfl,psfl,pctsrf,paire)
+#endif
+
+#ifdef INCAINFO
+           PRINT *, 'Appel CHEMHOOK_BEGIN ...'
+#endif
+           CALL chemhook_begin (calday,
+     $                          pctsrf(1,3),
+     $                          rlat,
+     $                          rlon,
+     $                          paire,
+     $                          paprs,
+     $                          pplay,
+     $                          ycoefh,
+     $                          pphi,
+     $                          t_seri,
+     $                          u,
+     $                          v,
+     $                          wo,
+     $                          q_seri,
+     $                          zxtsol,
+     $                          zxsnow,
+     $                          solsw,
+     $                          albsol,
+     $                          rain_fall,
+     $                          snow_fall,
+     $                          itop_con, 
+     $                          ibas_con, 
+     $                          cldfra,
+     $                          iim, 
+     $                          jjm, 
+     $                          tr_seri)
+#ifdef INCAINFO
+           PRINT *, 'OK.'
+#endif
+#endif
+
 c
 c Calculer les parametres optiques des nuages et quelques
@@ -2124,12 +2232,18 @@
 C
       call phytrac (rnpb,
+     I                   itap, julien, gmtime, 
      I                   debut,lafin,
      I                   nqmax-2,
      I                   nlon,nlev,dtime,
-     I                   t,paprs,pplay,
+     I                   u,v,t,paprs,pplay,
      I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      I                   ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
      I                   frac_impa, frac_nucl,
-     I                   rlon,presnivs,paire,pphis,
+     I                   rlon,paire,pphis,pphi,
+     I                   albsol, 
+     I                   qx(1,1,1), rhcl, 
+     I                   cldfra, rneb, diafra, cldliq, itop_con,
+     I                   ibas_con,
+     I                   pmflxr,pmflxs,prfl,psfl,flxmass_w,
      O                   tr_seri)
       ENDIF 
@@ -2822,4 +2936,24 @@
 cc         CALL ecrirega(84,d_q_lsc)
       ENDIF
+
+#ifdef INCA
+#ifdef INCAINFO
+           PRINT *, 'Appel CHEMHOOK_END ...'
+#endif
+           CALL chemhook_end (calday,
+     $                        dtime,
+     $                        pplay,
+     $                        t_seri,
+     $                        tr_seri,
+     $                        nbtr,
+     $                        paprs,
+     $                        anne_ini,
+     $                        day_ini,
+     $                        xjour)
+#ifdef INCAINFO
+           PRINT *, 'OK.'
+#endif
+#endif
+
 c
 c Convertir les incrementations en tendances
Index: LMDZ.3.3/trunk/libf/phylmd/phytrac.F
===================================================================
--- LMDZ.3.3/trunk/libf/phylmd/phytrac.F	(revision 343)
+++ LMDZ.3.3/trunk/libf/phylmd/phytrac.F	(revision 344)
@@ -3,14 +3,67 @@
 c
       SUBROUTINE phytrac (rnpb,
-     I                   debutphy,lafin,
-     I                   nqmax,
-     I                   nlon,nlev,pdtphys,
-     I                   t_seri,paprs,pplay,
-     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     I                   coefh,yu1,yv1,ftsol,pctsrf,xlat,
-     I                   frac_impa, frac_nucl,
-     I                   xlon,presnivs,paire,pphis,
-     O                   tr_seri)
+     I                    nstep, 
+     I                    julien, 
+     I                    gmtime,
+     I                    debutphy,
+     I                    lafin,
+     I                    nqmax,
+     I                    nlon,
+     I                    nlev,
+     I                    pdtphys,
+     I                    u,
+     I                    v,
+     I                    t_seri,
+     I                    paprs,
+     I                    pplay,
+     I                    pmfu, 
+     I                    pmfd, 
+     I                    pen_u, 
+     I                    pde_u, 
+     I                    pen_d, 
+     I                    pde_d,
+     I                    coefh,
+     I                    yu1,
+     I                    yv1,
+     I                    ftsol,
+     I                    pctsrf,
+     I                    xlat,
+     I                    frac_impa, 
+     I                    frac_nucl,
+     I                    xlon,
+     I                    paire,
+     I                    pphis,
+     I                    pphi,
+     I                    albsol,
+     I                    sh,
+     I                    rh,
+     I                    cldfra, 
+     I                    rneb, 
+     I                    diafra, 
+     I                    cldliq, 
+     I                    itop_con,
+     I                    ibas_con,
+     I                    pmflxr,
+     I                    pmflxs,
+     I                    prfl,
+     I                    psfl,
+     I                    flxmass_w,
+     O                    tr_seri)
       USE ioipsl
+
+#ifdef INCA
+      USE sflx
+      USE chem_tracnm
+      USE species_names
+      USE chem_mods
+      USE pht_tables, ONLY : jrates
+      USE transport_controls, ONLY : conv_flg, pbl_flg
+      USE airplane_src, ONLY : ptrop
+      USE lightning, ONLY : prod_light
+#ifdef INCA_AER
+      USE AEROSOL_DIAG
+#endif
+#endif
+
       IMPLICIT none
 c======================================================================
@@ -28,6 +81,8 @@
 #include "YOMCST.h"
 #include "dimensions.h"
+#include "paramet.h"
 #include "dimphy.h"
 #include "indicesol.h"
+#include "comvert.h"
 #include "temps.h"
 #include "control.h"
@@ -42,15 +97,34 @@
 c   -------
 c
-      integer nlon  ! nombre de points horizontaux
-      integer nlev  ! nombre de couches verticales
-      integer nqmax ! nombre de traceurs auxquels on applique la physique
+      integer nlon   ! nombre de points horizontaux
+      integer nlev   ! nombre de couches verticales
+      integer nqmax  ! nombre de traceurs auxquels on applique la physique
+      integer nstep  ! appel physique
+      integer julien !jour julien
+      integer itop_con(nlon) 
+      integer ibas_con(nlon) 
+      real gmtime
       real pdtphys  ! pas d'integration pour la physique (seconde)
       real t_seri(nlon,nlev) ! temperature
+      real u(klon,klev)
+      real v(klon,klev)
+      real sh(nlon,nlev)     ! humidite specifique
+      real rh(nlon,nlev)     ! humidite relative
+      real cldliq(nlon,nlev) ! eau liquide nuageuse
+      real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
+      real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
+      real rneb(nlon,nlev)   ! fraction nuageuse (grande echelle)
+      real albsol(nlon)      ! albedo surface
       real tr_seri(nlon,nlev,nbtr) ! traceur  
       real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
+      real ps(nlon)  ! pression surface
       real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
-      real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
+      real pphi(nlon,klev) ! geopotentiel
+      real znivsig(klev) ! niveaux sigma
       real paire(klon)
       real pphis(klon)
+      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
+      REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale
+      REAL flxmass_w(klon,klev)
       logical debutphy       ! le flag de l'initialisation de la physique
       logical lafin          ! le flag de la fin de la physique
@@ -131,8 +205,8 @@
       SAVE scavtr
 cAA
-      CHARACTER*1 itn
+      CHARACTER*2 itn
 C maf ioipsl
       CHARACTER*2 str2
-      INTEGER nhori, nvert
+      INTEGER nhori, nvert, nverta, nvertb, nverts
       REAL zsto, zout, zjulian
       INTEGER nid_tra
@@ -207,8 +281,34 @@
       data first,couchelimite,convection,lessivage,sorties
      s     /.true.,.true.,.true.,.true.,.true./
+
+#ifdef INCA
+      INTEGER           :: ncsec
+      INTEGER           :: prt_flag_ts(nbtr)=(/1,1,1
+#ifdef INCA_CH4
+     .                                              ,0,0,1,1,1,1,1,
+     .                                         0,1,0,0,0,0,0,1,0,0,
+     .                                         0,1,1,1,1,0,1,1,1,0,
+     .                                         1,1,1,1,1,1,1,1,1,1,
+     .                                         1,0,0
+#endif
+#ifdef INCA_AER
+     .                                              ,1,1,1,1,1,1
+#endif
+     .                                         /)
+      REAL, PARAMETER   :: dry_mass = 28.966
+      REAL, POINTER     :: hbuf(:)
+      REAL, ALLOCATABLE :: obuf(:)
+      REAL              :: calday
+      REAL              :: pdel(klon,klev)
+      REAL              :: dummy(klon,klev) = 0.
+#endif
+
 c
 c======================================================================
 c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
          modname='phytrac'
+         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
+!DH      print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
+         ps(:) = paprs(:,1) 
 
 c        print*,'DANS PHYTRAC debutphy=',debutphy
@@ -228,11 +328,14 @@
 
          inirnpb=rnpb
-         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
-         itra=0
-	 itap=0
+!DH      PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
+
 C         
          CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
-         zjulian = zjulian + day_ini
-c
+         itra = NINT(FLOAT(day_ini)*86400./pdtphys)
+         itap = NINT(FLOAT(day_ini)*86400./pdtphys)
+!        itra=0
+!        itap=0
+!        zjulian = zjulian + day_ini
+ 
          CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
          DO i = 1, iim
@@ -240,13 +343,21 @@
             zx_lon(i,jjm+1) = xlon(i+1)
          ENDDO
-c         DO ll=1,klev
-c            znivsig(ll)=float(ll)
-c         ENDDO
+         DO ll=1,klev
+            znivsig(ll)=float(ll)
+         ENDDO
          CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
          CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
-     .                 1,iim,1,jjm+1, 0, zjulian, pdtphys,
+     .                 1,iim,1,jjm+1, itra, zjulian, pdtphys,
      .                 nhori, nid_tra)
-         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
+
+         call histvert(nid_tra, "presnivs", "presnivs", "mb",
      .                 klev, presnivs, nvert)
+#ifdef INCA
+!        call histvert(nid_tra, "ap", "Hybrid A parameter", "-",
+!    .                 klev+1, ap, nverta)
+!        call histvert(nid_tra, "bp", "Hybrid B parameter", "-",
+!    .                 klev+1, bp, nvertb)
+#endif
+
          zsto = pdtphys
          zout = pdtphys * FLOAT(ecrit_tra)
@@ -260,70 +371,138 @@
      .                "once",  zsto,zout)
 
-	goto 666
-         CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-
-         CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "psrf1", "nature sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "psrf2", "nature sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "psrf3", "nature sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "psrf4", "nature sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst",  zsto,zout)
-         CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
-     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
-     .                "inst(X)",  zsto,zout)
-         CALL histdef(nid_tra, "pplay", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "inst(X)", zsto,zout)
-         CALL histdef(nid_tra, "t", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "inst(X)", zsto,zout)
-         CALL histdef(nid_tra, "mfu", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "mfd", "flux u decen","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "en_u", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "en_d", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "de_u", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "de_d", "flux u mont","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-         CALL histdef(nid_tra, "coefh", "turbulent coef","-",
-     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
-     .                "ave(X)", zsto,zout)
-
-666	continue
-c
+#ifdef INCA
+         CALL histdef(nid_tra, "ps", "Surface pressure", "Pa",
+     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "ptrop", "Tropopause pressure", "Pa",
+     .                iim,jjm+1,nhori, 1,1,1,-99, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "temp", "Air temperature", "K",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "u", "zonal wind component", "m/s",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "v", "zonal wind component", "m/s",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "h2o", "Specific Humidity", "MMR",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "pmid", "Pressure", "Pa",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+         CALL histdef(nid_tra, "pdel", "Delta Pressure", "Pa",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+#ifdef INCA_CH4
+#ifdef INCAINFO
+         DO it=1, phtcnt
+         WRITE(str2,'(i2.2)') it
+         CALL histdef(nid_tra, "j"//str2,"j"//str2, "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         ENDDO
+         DO it=1, hetcnt
+         WRITE(str2,'(i2.2)') it
+         CALL histdef(nid_tra, "w"//str2,"w"//str2, "S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         ENDDO
+
+         DO it=1, extcnt
+         WRITE(str2,'(i2.2)') it
+        CALL histdef(nid_tra, "ext"//str2,"ext"//str2, "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         ENDDO
+
+         DO it=1, nfs
+         WRITE(str2,'(i2.2)') it
+         CALL histdef(nid_tra, "INV"//str2, "INV"//str2, "CM-3",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         ENDDO
+
+#else
+         CALL histdef(nid_tra, "jO3","jO3", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "jNO2","jNO2", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "jH2O2","jH2O2", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "wHNO3","wHNO3", "S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "kN2O5", "kN2O5","CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "LghtNO","LghtNO", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+#endif
+
+         DO it=1, grpcnt
+         CALL histdef(nid_tra, grpsym(it), grpsym(it), "VMR",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         ENDDO
+#endif
+#endif
+
+#ifdef INCA_AER
+C   3d FIELDS
+
+        CALL histdef(nid_tra, "scavcoef_st","scavcoef_st", "S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+        CALL histdef(nid_tra, "scavcoef_cv","scavcoef_cv", "S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+#endif
+ 
          DO it=1,nqmax
          IF (it.LE.99) THEN
          WRITE(str2,'(i2.2)') it
-C champ 3D
+C champ 2D
+
+#ifdef INCA
+         IF ( prt_flag_ts(it) == 0 ) CYCLE
+
+         CALL histdef(nid_tra, "Emi_"//solsym(it), "Emi_"//solsym(it),
+     .           "kg/m2/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "Dep_"//solsym(it), "Dep_"//solsym(it),
+     .           "cm/s", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+#ifdef INCA_AER
+        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
+     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
+     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
+         CALL histdef(nid_tra, "Sed_"//solsym(it), "Sed_"//solsym(it),
+     .           "kg/m2", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+        endif
+        IF (solsym(it) .eq. 'CIDUSTM') THEN
+         CALL histdef(nid_tra, "OD_"//solsym(it), "OD_"//solsym(it),
+     .           "opt. depth", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+        endif
+
+#endif
+         CALL histdef(nid_tra, solsym(it), solsym(it), "VMR",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+#else
          CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
      .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
@@ -334,4 +513,5 @@
      .                "ave(X)", zsto,zout)
          ENDIF
+#endif
          ELSE
          PRINT*, "Trop de traceurs"
@@ -339,4 +519,38 @@
          ENDIF
          ENDDO
+
+#ifdef INCA_CH4
+         CALL histdef(nid_tra, "O3_column", "O3_column",
+     .           "DU", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "CO_column", "CO_column",
+     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "CH4_column", "CH4_column",
+     .           "10^18 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "NO2_column", "NO2_column",
+     .           "10^15 CM-2", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "O3_ste", "O3_ste",
+     .           "CM-2 S-1", iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .           "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "O3_prod", "O3_prod", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "O3_loss", "O3_loss", "CM-3 S-1",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+!        Special variables for daytime averaging
+!        CALL histdef(nid_tra, "day_cnt", "day_cnt", "-",
+!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+!    .                "t_sum(X)", zsto,zout)
+!        CALL histdef(nid_tra, "NO_day", "NO_day", "VMR",
+!    .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+!    .                "t_sum(X)", zsto,zout)
+
+#endif
+
          CALL histend(nid_tra)
          ndex(1) = 0
@@ -367,5 +581,5 @@
          read(99,*) (trs(i,1),i=1,klon)
 999      close(99)
-         print*, 'apres starttrac'
+!DH      print*, 'apres starttrac'
 
 c Initialisation de la fraction d'aerosols lessivee
@@ -402,67 +616,91 @@
          inirnpb=.false.
       endif
-      if(nqmax.gt.2) aerosol(3)=.true.
-
-
-c  abder
-	goto 777
-            do i=1,nlon
-               pftsol1(i) = ftsol(i,1)
-               pftsol2(i) = ftsol(i,2)
-               pftsol3(i) = ftsol(i,3)
-               pftsol4(i) = ftsol(i,4)
-
-               ppsrf1(i) = pctsrf(i,1)
-               ppsrf2(i) = pctsrf(i,2)
-               ppsrf3(i) = pctsrf(i,3)
-               ppsrf4(i) = pctsrf(i,4)
-
-            enddo
-	 ndex(1)=0
-	 itap=itap+1
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
-         CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
-     s                                  iim*(jjm+1),ndex)
-         
-	 CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
-         CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
-     s                                  iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
-         CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
-     s                                       iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
-         CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
-     s                                       iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
-         CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
-     s                                      iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
-         CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
-     s                                      iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
-         CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
-     s                                     iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
-         CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d, 
-     s                                     iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
-         CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
-     s                                     iim*(jjm+1),ndex)
-
-         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
-         CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
-     s                                     iim*(jjm+1),ndex)
-777	continue
+
+#ifdef INCA
+!======================================================================
+!     Chimie
+!======================================================================
+
+        calday = FLOAT(julien) + gmtime
+        ncsec  = NINT (86400.*gmtime)
+
+        DO k = 1, nlev
+        pdel(:,k) = paprs(:,k) - paprs (:,k+1)
+        END DO
+
+#ifdef INCAINFO
+        PRINT *, 'CHEMMAIN @ ', calday, ' ... '
+        DO it = 1, nbtr
+        PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
+     $                       MAXVAL(tr_seri(:,:,it))
+      END DO
+#endif
+
+
+#ifdef INCA_AER
+        CALL aerosolmain (tr_seri,
+     $                 pdtphys,
+     $                 pplay,
+     $                 prfl,
+     $                 pmflxr,
+     $                 psfl,
+     $                 pmflxs,
+     $                 pmfu,
+     $                 itop_con,
+     $                 ibas_con,
+     $                 pphi,
+     $                 paire,
+     $                 nstep)
+#endif 
+
+        CALL chemmain (tr_seri,
+     $                 nas,
+     $                 nstep,
+     $                 calday,
+     $                 julien,
+     $                 ncsec,
+     $                 1,
+     $                 pdtphys,
+     $                 paprs(1,1),
+     $                 pplay,
+     $                 pdel,
+     $                 pctsrf(1,3),
+     $                 ftsol,
+     $                 albsol,
+     $                 pphi,
+     $                 pphis,
+     $                 cldfra,
+     $                 rneb,
+     $                 diafra,
+     $                 itop_con,
+     $                 cldliq,
+     $                 prfl,
+     $                 pmflxr,
+     $                 psfl,
+     $                 pmflxs,
+     $                 pmfu,
+     $                 flxmass_w,
+     $                 t_seri,
+     $                 sh,
+     $                 rh,
+     $                 .false.,
+     $                 hbuf,
+     $                 obuf,
+     $                 iip1,
+     $                 jjp1)
+#ifdef INCAINFO
+      PRINT *, 'OK.'
+      DO it = 1, nbtr
+      PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
+     $                     MAXVAL(tr_seri(:,:,it))
+      END DO
+#endif
+#endif
+C
+
 c======================================================================
 c   Calcul de l'effet de la convection
 c======================================================================
-	print*,'Avant convection'
+!DH   print*,'Avant convection'
       do it=1,nqmax
          WRITE(itn,'(i1)') it
@@ -472,6 +710,9 @@
       if (convection) then
 
-      print*,'Pas de temps dans phytrac : ',pdtphys
+!DH   print*,'Pas de temps dans phytrac : ',pdtphys
       DO it=1, nqmax
+#ifdef INCA
+      IF ( conv_flg(it) == 0 ) CYCLE
+#endif
       CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      .            pplay, paprs, tr_seri(1,1,it), d_tr_cv)
@@ -481,9 +722,12 @@
       ENDDO
       ENDDO
-c      WRITE(itn,'(i1)') it
-c      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
-      ENDDO
-c      print*,'apres nflxtr'
-
+c     WRITE(itn,'(i2)') it
+#ifdef INCA
+      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//solsym(it))
+#else
+      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
+#endif
+      ENDDO
+c     print*,'apres nflxtr'
 
       endif ! convection
@@ -513,5 +757,8 @@
 C maf modif pour tenir compte du cas rnpb + traceur
       DO it=1, nqmax
-c     print *,'it',it,clsol(it)
+#ifdef INCA
+      IF ( pbl_flg(it) == 0 ) CYCLE
+#endif
+C     print *,'it',it,clsol(it)
       if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
           CALL cltracrn(it, pdtphys, yu1, yv1,
@@ -539,9 +786,15 @@
 C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
       else ! couche limite avec flux prescrit
+
+#ifdef INCA
+        DO k =  1, klon
+          source(k) = eflux(k,it)-dflux(k,it)
+        END DO
+#else
 Cmaf provisoire source / traceur a creer
         DO i=1, klon
           source(i) = 0.0 ! pas de source, pour l'instant
         ENDDO
-C
+#endif
           CALL cltrac(pdtphys, coefh,t_seri,
      s               tr_seri(1,1,it), source,
@@ -574,4 +827,5 @@
 c si radio=true mais pour l'instant radiornpb propre au cas rnpb
       if(rnpb) then
+        print *, 'decroissance radiactive activee'
         call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
 C
@@ -589,10 +843,10 @@
 c
       endif ! rnpb decroissance  radioactive
-C
+
 c======================================================================
 c   Calcul de l'effet de la precipitation
 c======================================================================
 
-      print*,'LESSIVAGE =',lessivage
+!DH   print*,'LESSIVAGE =',lessivage
       IF (lessivage) THEN
 
@@ -675,18 +929,199 @@
          ENDDO
       ENDDO
+
+
+
       itra=itra+1
       ndex(1) = 0
+
+      i = NINT(zout/zsto)
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, paire,zx_tmp_2d)
+      CALL histwrite(nid_tra,"aire",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+#ifdef INCA
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ps,zx_tmp_2d)
+      CALL histwrite(nid_tra,"ps",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ptrop,zx_tmp_2d)
+      CALL histwrite(nid_tra,"ptrop",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri, zx_tmp_3d)
+      CALL histwrite(nid_tra,"temp",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,u, zx_tmp_3d)
+      CALL histwrite(nid_tra,"u",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,v, zx_tmp_3d)
+      CALL histwrite(nid_tra,"v",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,sh, zx_tmp_3d)
+      CALL histwrite(nid_tra,"h2o",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pdel, zx_tmp_3d)
+      CALL histwrite(nid_tra,"pdel",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay, zx_tmp_3d)
+      CALL histwrite(nid_tra,"pmid",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+#ifdef INCA_CH4
+#ifdef INCAINFO
+      DO it=1, phtcnt
+      WRITE(str2,'(i2.2)') it
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,it),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"j"//str2,itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+      ENDDO
+
+      DO it=1, hetcnt
+      WRITE(str2,'(i2.2)') it
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,it),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"w"//str2,itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+      ENDDO
+
+ 
+      DO it=1, extcnt
+      WRITE(str2,'(i2.2)') it
+
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,it),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"ext"//str2,itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+      ENDDO
+
+      DO it=1, nfs
+      WRITE(str2,'(i2.2)') it
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,invariants(1,1,it),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"INV"//str2,itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+      ENDDO
+
+
+#else
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,2),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"jO3",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,4),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"jNO2",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,jrates(1,1,13),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"jH2O2",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,hrates(1,1,1),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"wHNO3",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,krates(1,1,1),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"kN2O5",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,extflx(1,1,1),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"LghtNO",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+#endif
+      DO it=1, grpcnt
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,nas(1,1,it),zx_tmp_3d)
+      zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(it)
+      CALL histwrite(nid_tra,grpsym(it),itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+      ENDDO
+#endif
+
+#endif
+#ifdef INCA_AER
+C   3d FIELDS
+
+      it = id_CIDUSTM
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_st(1,1,it),
+     .                  zx_tmp_3d)
+       CALL histwrite(nid_tra,"scavcoef_st",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,scavcoef_cv(1,1,it),
+     .                  zx_tmp_3d)
+       CALL histwrite(nid_tra,"scavcoef_cv",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+
+#endif 
+
       DO it=1,nqmax
       IF (it.LE.99) THEN
        WRITE(str2,'(i2.2)') it
+
+#ifdef INCA
+      IF ( prt_flag_ts(it) == 0 ) CYCLE
+C champs 2D
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, eflux(1,it),zx_tmp_2d)
+      CALL histwrite(nid_tra,"Emi_"//solsym(it),itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, dvel(1,it),zx_tmp_2d)
+      CALL histwrite(nid_tra,"Dep_"//solsym(it),itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+#ifdef INCA_AER
+        IF (solsym(it) .eq. 'CIDUSTM'.or.solsym(it) .eq. 'CIN'
+     .  .or.solsym(it) .eq. 'CSSSM'  .or.solsym(it) .eq. 'CSN'
+     .  .or.solsym(it) .eq. 'ASSSM'  .or.solsym(it) .eq. 'ASN' ) THEN
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1,sflux(1,it),zx_tmp_2d)
+      CALL histwrite(nid_tra,"Sed_"//solsym(it),itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+        endif
+        IF (solsym(it) .eq. 'CIDUSTM') THEN
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1,tausum(1),zx_tmp_2d)
+      CALL histwrite(nid_tra,"OD_"//solsym(it),itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+        endif
+
+#endif 
+C champs 3D
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
+
+       !Prefer vmr to mmr for transported species
+       if( adv_mass(it) /= 0. ) then
+       zx_tmp_3d = zx_tmp_3d * dry_mass / adv_mass(it)
+       else
+#ifdef INCA_CH4
+       if ( solsym(it) == 'OX' ) then
+       zx_tmp_3d = zx_tmp_3d * dry_mass / nadv_mass(id_o3)
+       end if
+#endif
+       end if
+
+       CALL histwrite(nid_tra,solsym(it),itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+#else
+
 C champs 3D
        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
        CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
      .                                   iim*(jjm+1)*klev,ndex)
-c      IF (lessivage) THEN
-c      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
-c      CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
-c    .                                   iim*(jjm+1)*klev,ndex)
-c     ENDIF
+       IF (lessivage) THEN
+       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
+       CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+       ENDIF
+#endif
       ELSE
          PRINT*, "Trop de traceurs"
@@ -695,44 +1130,51 @@
       ENDDO
 
-	goto 888
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
-        CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
-        CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
-        CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
-        CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
-        CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
-        CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
-        CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
-        CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
-        CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
-     .                  iim*(jjm+1)*klev,ndex)
-
-888	continue
-
-c	print*,'Sortie phytrac'
-c      do it=1,nqmax
-c         WRITE(itn,'(i1)') it
-c        call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys  '//itn)
-c      enddo
+#ifdef INCA_CH4
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_tr_col(1), zx_tmp_2d)
+      CALL histwrite(nid_tra,"O3_column",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, co_tr_col(1), zx_tmp_2d)
+      CALL histwrite(nid_tra,"CO_column",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, ch4_tr_col(1), zx_tmp_2d)
+      CALL histwrite(nid_tra,"CH4_column",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, no2_tr_col(1), zx_tmp_2d)
+      CALL histwrite(nid_tra,"NO2_column",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(1, klon,iim,jjm+1, o3_st_flx(1), zx_tmp_2d)
+      CALL histwrite(nid_tra,"O3_ste",itra,zx_tmp_2d,
+     .     iim*(jjm+1),ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_prod(1,1),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"O3_prod",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,o3_loss(1,1),
+     .     zx_tmp_3d)
+      CALL histwrite(nid_tra,"O3_loss",itra,zx_tmp_3d,
+     .                                   iim*(jjm+1)*klev,ndex)
+
+!     ... Special section for daytime averaging
+!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,day_cnt(1,1),
+!    .       zx_tmp_3d)
+!       CALL histwrite(nid_tra,"day_cnt",itra,zx_tmp_3d,
+!    .                                  iim*(jjm+1)*klev,ndex)
+!       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,no_daytime(1,1),
+!    .       zx_tmp_3d)
+!       CALL histwrite(nid_tra,"NO_day",itra,zx_tmp_3d,
+!    .                                  iim*(jjm+1)*klev,ndex)
+
+#endif
+
+      if (ok_sync) call histsync(nid_tra)
 
       if (lafin) then
-         print*, 'c est la fin de la physique'
+!DH      print*, 'c est la fin de la physique'
          open (99,file='restarttrac',  form='formatted')
          do i=1,klon
@@ -742,5 +1184,5 @@
          close(99)
       else
-         print*, 'physique pas fini'
+!DH      print*, 'physique pas fini'
       endif
 
