c
c $Header$
c
      SUBROUTINE phytrac (rnpb,
     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======================================================================
c Auteur(s) FH
c Objet: Moniteur general des tendances traceurs
c
cAA Remarques en vrac:
cAA--------------------
cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien 
cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
cAA 2/ Le choix du radon et du pb se fait juste avec un data 
cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir 
cAA    une variable "type de traceur" 
c======================================================================
#include "YOMCST.h"
#include "dimensions.h"
#include "paramet.h"
#include "dimphy.h"
#include "indicesol.h"
#include "comvert.h"
#include "temps.h"
#include "control.h"
c======================================================================

c Arguments:
c
c   EN ENTREE:
c   ==========
c
c   divers:
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 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 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

      integer ll
c
cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
c
c   convection:
c   -----------
c
      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
      REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
      REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
      REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
c
c   Couche limite:
c   --------------
c
      REAL coefh(nlon,nlev) ! coeff melange CL
      REAL yu1(nlon)        ! vents au premier niveau
      REAL yv1(nlon)        ! vents au premier niveau
      REAL xlat(nlon)       ! latitudes pour chaque point 
      REAL xlon(nlon)       ! longitudes pour chaque point 

c
c   Lessivage:
c   ----------
c
c pour le ON-LINE
c
      REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
      REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
c
cAA
cAA Arguments necessaires pour les sources et puits de traceur:
cAA ----------------
cAA
      real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
      real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
c abder
      real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
      real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
c fin
cAA ----------------------------
cAA  VARIABLES LOCALES TRACEURS
cAA ----------------------------
cAA
cAA Sources et puits des traceurs:
cAA ------------------------------
cAA
cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.

      REAL source(klon)       ! a voir lorsque le flux est prescrit 
      REAL topuits(klev,nbtr)  ! a voir 
cAA 
cAA Pour la source de radon et son reservoir de sol
cAA ................................................
 
      REAL trs(klon,nbtr)    ! Conc. radon ds le sol
      SAVE trs

      REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
c                            Masque de l'echange avec la surface
c                           (1 = reservoir) ou (possible => 1 )
      SAVE masktr
      REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
      SAVE fshtr
      REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
      SAVE hsoltr
      REAL tautr(nbtr)       ! Constante de decroissance radioactive
      SAVE tautr
      REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
      SAVE vdeptr
      REAL scavtr(nbtr)      ! Coefficient de lessivage
      SAVE scavtr
cAA
      CHARACTER*2 itn
C maf ioipsl
      CHARACTER*2 str2
      INTEGER nhori, nvert, nverta, nvertb, nverts
      REAL zsto, zout, zjulian
      integer idayref
      INTEGER nid_tra
      SAVE nid_tra
c     REAL x(klon,klev,nbtr+2) ! traceurs 
      INTEGER ndex(1)
      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
      INTEGER itra
      SAVE itra                   ! compteur pour la physique
C
C Variables liees a l'ecriture de la bande histoire : phytrac.nc
c
      INTEGER ecrit_tra
      SAVE ecrit_tra   
      logical ok_sync
      parameter (ok_sync = .true.)
C
C nature du traceur
c
      logical aerosol(nbtr)  ! Nature du traceur
c                            ! aerosol(it) = true  => aerosol 
c                            ! aerosol(it) = false => gaz 
c                            ! nat_trac(it) = 1. aerosolc
      logical clsol(nbtr)    ! clsol(it) = true => CL sol calculee
      logical radio(nbtr)    ! radio(it)=true => decroisssance radioactive
      save aerosol,clsol,radio
C
c======================================================================
c
c Declaration des procedures appelees
c
c--modif convection tiedtke
      INTEGER i, k, it,itap
	save itap
      REAL delp(klon,klev)
c--end modif
c
c Variables liees a l'ecriture de la bande histoire physique
c
c Variables locales pour effectuer les appels en serie
c----------------------------------------------------
c
      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs 
      REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection 
      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance 
c                                   ! radioactive du rn - > pb 
      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage 
c                                          ! par impaction
      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage 
c                                          ! par nucleation 
      REAL fluxrn(klon,klev) 
      REAL fluxpb(klon,klev) 
      REAL pbimpa(klon,klev) 
      REAL pbnucl(klon,klev) 
      REAL rn(klon,klev) 
      REAL pb(klon,klev) 
      REAL flestottr(klon,klev,nbtr) ! flux de lessivage 
c                                    ! dans chaque couche 

C
      character*20 modname
      character*80 abort_message
c
c   Controles
c-------------
      logical first,couchelimite,convection,lessivage,sorties,
     s        rnpb,inirnpb
      save first,couchelimite,convection,lessivage,sorties,
     s     inirnpb
      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

         if (debutphy) then

          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
          ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
c         ecrit_tra = NINT(86400./pdtphys) ! tous les 24H

         if(nbtr.lt.nqmax) then
           print*,'NQMAX=',nqmax
           print*,'NBTR=',nbtr
           abort_message='See above'
           call abort_gcm(modname,abort_message,1)
         endif

         inirnpb=rnpb
!DH      PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra

C         
         idayref = 1
         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
         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
            zx_lon(i,1) = xlon(i+1)
            zx_lon(i,jjm+1) = xlon(i+1)
         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, itra, zjulian, pdtphys,
     .                 nhori, nid_tra)

         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)
c
         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
     .                "once",  zsto,zout)
c
         CALL histdef(nid_tra, "aire", "Grid area", "-",
     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
     .                "once",  zsto,zout)

#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 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,
     .                "ave(X)", zsto,zout)
         IF (lessivage) THEN
         CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2, 
     .                "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
     .                "ave(X)", zsto,zout)
         ENDIF
#endif
         ELSE
         PRINT*, "Trop de traceurs"
         CALL abort
         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
c
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
C
         i = NINT(zout/zsto)
         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)

c======================================================================
c   Initialisation de certaines variables pour le Rn et le Pb 
c======================================================================

c Initialisation du traceur dans le sol (couche limite radonique)
c
c        print*,'valeur de debut dans phytrac :',debutphy
         do it=1,nqmax
            do i=1,klon
               trs(i,it) = 0.
            enddo
         END DO

         open (99,file='starttrac',status='old',
     .         err=999,form='formatted')
         read(99,*) (trs(i,1),i=1,klon)
999      close(99)
!DH      print*, 'apres starttrac'

c Initialisation de la fraction d'aerosols lessivee
c
         DO it=1,nqmax
           DO k = 1, nlev
              DO i = 1, klon
                 d_tr_lessi_impa(i,k,it) =  0.0
                 d_tr_lessi_nucl(i,k,it) =  0.0 
             ENDDO
           ENDDO
         ENDDO
c
c Initialisation de la nature des traceurs
c
         DO it = 1, nqmax
            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
         ENDDO
c
      ENDIF  ! fin debutphy 
c Initialisation du traceur dans le sol (couche limite radonique)
      if(inirnpb) THEN
c
         radio(1)= .true.
         radio(2)= .true.
         clsol(1)= .true.
         clsol(2)= .true.
         aerosol(2) = .TRUE. ! le Pb est un aerosol 
c
         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
     .                   ,vdeptr,scavtr)
         inirnpb=.false.
      endif

#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======================================================================
!DH   print*,'Avant convection'
      do it=1,nqmax
         WRITE(itn,'(i1)') it
c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
      enddo

      if (convection) then

!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)
      DO k = 1, nlev
      DO i = 1, klon
         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
      ENDDO
      ENDDO
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
c        print*,'Apres convection'
c      do it=1,nqmax
c         WRITE(itn,'(i1)') it
c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
c      enddo

c======================================================================
c   Calcul de l'effet de la couche limite
c======================================================================
c	print *,'Avant couchelimite'
c      do it=1,nqmax
c         WRITE(itn,'(i1)') it
c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
c      enddo

      if (couchelimite) then

      DO k = 1, nlev
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO

C maf modif pour tenir compte du cas rnpb + traceur
      DO it=1, nqmax
#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,
     e                    coefh,t_seri,ftsol,pctsrf,
     e                    tr_seri(1,1,it),trs(1,it),
     e                    paprs, pplay, delp,
     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
     e                    tautr(it),vdeptr(it),
     e                    xlat,
     s                    d_tr_cl,d_trs)
          DO k = 1, nlev
            DO i = 1, klon
              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
            ENDDO
          ENDDO
c
c Traceur ds sol
c
          DO i = 1, klon
            trs(i,it) = trs(i,it) + d_trs(i)
          END DO
C
C maf provisoire suppression des prints
C         WRITE(itn,'(i1)') it
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
#endif
          CALL cltrac(pdtphys, coefh,t_seri,
     s               tr_seri(1,1,it), source,
     e               paprs, pplay, delp,
     s               d_tr )
            DO k = 1, nlev
               DO i = 1, klon
                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
               ENDDO
            ENDDO
Cmaf provisoire suppression des prints
Cmaf          WRITE(itn,'(i1)') it
cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
      endif
      ENDDO
c
      endif ! couche limite

c      print*,'Apres couchelimite'
c      do it=1,nqmax
c         WRITE(itn,'(i1)') it
c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
c      enddo

c======================================================================
c   Calcul de l'effet du puits radioactif
c======================================================================

C MAF il faudrait faire une modification pour passer dans radiornpb 
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
        DO it=1,nqmax
            if(radio(it)) then
            DO k = 1, nlev
               DO i = 1, klon
                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
               ENDDO
            ENDDO
            WRITE(itn,'(i1)') it
            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
            endif
        ENDDO
c
      endif ! rnpb decroissance  radioactive

c======================================================================
c   Calcul de l'effet de la precipitation
c======================================================================

!DH   print*,'LESSIVAGE =',lessivage
      IF (lessivage) THEN

c     print*,'avant lessivage'

       DO it = 1, nqmax
           DO k = 1, nlev
              DO i = 1, klon
               d_tr_lessi_nucl(i,k,it) = 0.0 
               d_tr_lessi_impa(i,k,it) = 0.0 
               flestottr(i,k,it) = 0.0 
              ENDDO
           ENDDO
       ENDDO
c
c tendance des aerosols nuclees et impactes 
c
       DO it = 1, nqmax
         IF (aerosol(it)) THEN
           DO k = 1, nlev
              DO i = 1, klon
               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
              ENDDO
           ENDDO
         ENDIF
       ENDDO
c
c Mises a jour des traceurs + calcul des flux de lessivage 
c Mise a jour due a l'impaction et a la nucleation
c
c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
       DO it = 1, nqmax
c         print*,'IT=',it,aerosol(it)
         IF (aerosol(it)) THEN
c           print*,'IT=',it,' On lessive'
           DO k = 1, nlev
              DO i = 1, klon
               tr_seri(i,k,it)=tr_seri(i,k,it)
     s         *frac_impa(i,k)*frac_nucl(i,k)
              ENDDO
           ENDDO
         ENDIF
       ENDDO
c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
c
c Flux lessivage total 
c
      DO it = 1, nqmax
           DO k = 1, nlev
            DO i = 1, klon
               flestottr(i,k,it) = flestottr(i,k,it) -
     s                   ( d_tr_lessi_nucl(i,k,it)   +
     s                     d_tr_lessi_impa(i,k,it) ) *
     s                   ( paprs(i,k)-paprs(i,k+1) ) / 
     s                   (RG * pdtphys)
            ENDDO
           ENDDO
c
Cmaf suppression provisoire
Cmaf        WRITE(itn,'(i1)') it
Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
      ENDDO
c
c     print*,'apres lessivage'
      ENDIF
Cc
      DO k = 1, nlev
         DO i = 1, klon
            fluxrn(i,k) = flestottr(i,k,1)
            fluxpb(i,k) = flestottr(i,k,2)
            rn(i,k) = tr_seri(i,k,1)
            pb(i,k) = tr_seri(i,k,2)
            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
         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)
       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"
         CALL abort
      ENDIF
      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
!DH      print*, 'c est la fin de la physique'
         open (99,file='restarttrac',  form='formatted')
         do i=1,klon
             write(99,*) trs(i,1)
         enddo
         PRINT*, 'Ecriture du fichier restarttrac'
         close(99)
      else
!DH      print*, 'physique pas fini'
      endif


      RETURN
      END