!
! $Header$
!
c
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                    fm_therm,
     I                    entr_therm,
     I                    yu1,
     I                    yv1,
     I                    ftsol,
     I                    pctsrf,
     I                    xlat,
     I                    frac_impa,
     I                    frac_nucl,
     I                    xlon,
     I                    presnivs,
     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                    da,
     I                    phi,
     I                    mp,
     I                    upwd,
     I                    dnwd,
#ifdef INCA
     I                    flxmass_w,
#if defined(INCA_AER) && defined(CPP_COUPLE)
     I                   tau_inca, 
     I                   piz_inca, 
     I                   cg_inca,
     I                   ccm,
     I                   rfname,
#endif
#endif
     O                    tr_seri)

      USE ioipsl
      USE dimphy
      USE mod_grid_phy_lmdz
      USE mod_phys_lmdz_para
      USE comgeomphy
      USE iophy
      USE vampir

      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"
cym#include "dimphy.h"
#include "indicesol.h"
#include "clesphys.h"
#include "temps.h"
#include "paramet.h"
#include "control.h"
cym#include "comgeomphy.h"
#include "advtrac.h"
#include "thermcell.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 tr_seri(nlon,nlev,nbtr) ! traceur  
      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 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 pphis(klon)
      REAL presnivs(klev)
      logical debutphy       ! le flag de l'initialisation de la physique
      logical lafin          ! le flag de la fin de la physique
c Olivia      
      integer nsplit
      REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)   !--lessivage convection
      REAL prfl(klon,klev+1),   psfl(klon,klev+1)     !--lessivage large-scale

#ifdef INCA
      REAL flxmass_w(klon,klev)
      CHARACTER(len=8) :: solsym(nqmax)
#if defined(INCA_AER) && defined(CPP_COUPLE)
      integer la
      REAL              ::    tau_inca(klon,klev,9,2)
      REAL              ::    piz_inca(klon,klev,9,2)
      REAL              ::    cg_inca(klon,klev,9,2)
      character*4       ::    rfname(9) 
      REAL              ::    ccm(klon,klev,2) 
#endif
#endif
c      integer iflag_con

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

c
c   thermiques:
c   -----------
c
      real fm_therm(klon,klev+1),entr_therm(klon,klev)
	real fm_therm1(klon,klev)
c
      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 KE
      real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
      REAL upwd(nlon,nlev)      ! saturated updraft mass flux
      REAL dnwd(nlon,nlev)      ! saturated downdraft mass flux

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,nqmax)       ! a voir lorsque le flux est prescrit 
cAA 
cAA Pour la source de radon et son reservoir de sol
cAA ................................................
 
      REAL,save,allocatable :: trs(:,:)    ! Conc. radon ds le sol
c$OMP THREADPRIVATE(trs)
cym      SAVE trs
      REAL :: trs_tmp(klon_glo)
      
      REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur
c                            Masque de l'echange avec la surface
c                           (1 = reservoir) ou (possible => 1 )
c$OMP THREADPRIVATE(masktr)
cym      SAVE masktr
      REAL,save,allocatable :: fshtr(:,:)  ! Flux surfacique dans le reservoir de sol
c$OMP THREADPRIVATE(fshtr)
cym      SAVE fshtr
      REAL,save,allocatable :: hsoltr(:)      ! Epaisseur equivalente du reservoir de sol
c$OMP THREADPRIVATE(hsoltr)
cym      SAVE hsoltr
      REAL,save,allocatable :: tautr(:)       ! Constante de decroissance radioactive
c$OMP THREADPRIVATE(tautr)
cym      SAVE tautr
      REAL,save,allocatable :: vdeptr(:)      ! Vitesse de depot sec dans la couche Brownienne
c$OMP THREADPRIVATE(vdeptr)
cym      SAVE vdeptr
      REAL,save,allocatable :: scavtr(:)      ! Coefficient de lessivage
c$OMP THREADPRIVATE(scavtr)
cym      SAVE scavtr
cAA
      CHARACTER*2 itn
C maf ioipsl
      CHARACTER*2 str2
      INTEGER nhori, nvert
      REAL zsto, zout, zjulian
      INTEGER nid_tra
      SAVE nid_tra
c$OMP THREADPRIVATE(nid_tra)
#ifdef INCA_AER
      INTEGER nid_tra2,nid_tra3
      SAVE nid_tra2,nid_tra3
c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
#endif
c     REAL x(klon,klev,nbtr+2) ! traceurs 
      INTEGER ndex(1)
      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
      REAL zx_tmp_fi2d(klon)      ! variable temporaire grille physique
      REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D 
      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)
c
      integer itau_w   ! pas de temps ecriture = nstep + itau_phy
c

C
C Variables liees a l'ecriture de la bande histoire : phytrac.nc
c
c      INTEGER ecrit_tra
c      SAVE ecrit_tra 

      logical ok_sync
      parameter (ok_sync = .true.)
C
C nature du traceur
c
      logical,save,allocatable :: aerosol(:)  ! Nature du traceur
c                            ! aerosol(it) = true  => aerosol 
c                            ! aerosol(it) = false => gaz 
c                            ! nat_trac(it) = 1. aerosol
      logical,save,allocatable :: clsol(:)    ! clsol(it) = true => CL sol calculee
      logical,save,allocatable :: radio(:)    ! radio(it)=true => decroisssance radioactive
c$OMP THREADPRIVATE(aerosol,clsol,radio)  
cym      save aerosol,clsol,radio
C
c======================================================================
c
c Declaration des procedures appelees
c
c--modif convection tiedtke
      INTEGER i, k, it
      INTEGER iq, iiq
      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,nbtr) ! tendance de traceurs  couche limite
      REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs  conv pour chq traceur
      REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
      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 
      real zmasse(klon,klev)
      real ztra_th(klon,klev)

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,
     s        sorties,inirnpb
c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
c$OMP+              sorties,inirnpb)
c      data first,couchelimite,convection,lessivage,sorties
c     s     /.true.,.true.,.false.,.true.,.true./
c Olivia
       data first,couchelimite,convection,lessivage,
     s      sorties
     s     /.true.,.true.,.true.,.true.,.true./


#ifdef INCA
      INTEGER           :: lastgas
      INTEGER           :: ncsec
      INTEGER           :: prt_flag_ts(nbtr)  

      REAL, PARAMETER   :: dry_mass = 28.966
      REAL, POINTER     :: hbuf(:)
      REAL, ALLOCATABLE :: obuf(:)
      REAL              :: calday
      REAL              :: pdel(klon,klev)
      REAL              :: dummy(klon,klev)
#endif
c
c======================================================================

#ifdef INCA
      prt_flag_ts(:)=(/
#ifdef INCA_CH4
     .             1,1,1,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
#ifdef INCA_AER 
     .             ,1,1,1,1,0,1,1,1,1,0,
     .             1,1,1,1,1,1,0,1,0,1,
     .             1,1,1,1,0,1,0,1,1,1
#endif
#endif
#ifdef INCA_NMHC
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1
#ifdef INCA_AER 
     .             ,1,1,1,1,0,1,1,1,1,0,
     .             1,1,1,1,1,1,0,1,0,1,
     .             1,1,1,1,0,1,0,1,1,1
#endif
#endif
#if defined(INCA_AER) && !defined(INCA_CH4) && !defined(INCA_NMHC)
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1,1,
     .             1,1,1,1,1,1,1,1,1
#endif
#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
     .             1,1,1,1,1,1,1,1,1,1,
     .             1                                         
#endif

     .             /)
      dummy(:,:) = 0.

#endif
         modname='phytrac'

         ps(:)=paprs(:,1)

         if (debutphy) then
           allocate( trs(klon,nbtr) )
	   allocate( masktr(klon,nbtr)) 
           allocate( fshtr(klon,nbtr) )
           allocate( hsoltr(nbtr))
           allocate( tautr(nbtr))
           allocate( vdeptr(nbtr))
           allocate( scavtr(nbtr))
           allocate( aerosol(nbtr))
           allocate( clsol(nbtr))
           allocate( radio(nbtr))


c jg: c'est ca qu'on veut?????           
          ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy)) 
          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra

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

         inirnpb=rnpb
         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
C         
c=============================================================
c   Initialisation des sorties
c=============================================================

#ifdef CPP_IOIPSL
#include "ini_histrac.h"
#endif

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
         trs(:,:) = 0.
c$OMP MASTER         
       if (is_mpi_root) then
         trs_tmp(:)=0.
	 open (99,file='starttrac',status='old',
     .         err=999,form='formatted')
         read(99,*) (trs_tmp(i),i=1,klon_glo)
999      close(99)
       endif
c$OMP END MASTER
       call Scatter(trs_tmp,trs(:,1))

c         print*, 'apres starttrac'

c Initialisation de la fraction d'aerosols lessivee
c
         d_tr_lessi_impa(:,:,:) = 0.
         d_tr_lessi_nucl(:,:,:) = 0. 
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
      call VTe(VTphysiq)
      call VTb(VTinca)
!======================================================================
!     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,
     $                 pdel,
     $                 prfl,
     $                 pmflxr,
     $                 psfl,
     $                 pmflxs,
     $                 pmfu,
     $                 itop_con,
     $                 ibas_con,
     $                 pphi,
     $                 airephy, ! paire,
     $                 nstep,
     $                 rneb,         ! for chimiaq
     $                 t_seri,       ! for chimiaq
     $                 rh,
#ifdef CPP_COUPLE
     $                 tau_inca,
     $                 piz_inca,
     $                 cg_inca,
     $                 rfname,
     $                 ccm,
#endif
     $                 lafin)
#endif

        CALL chemmain (tr_seri,    !mmr
     $                 nstep,      !nstep
     $                 calday,     !calday
     $                 julien,     !ncdate
     $                 ncsec,      !ncsec
     $                 1,          !lat
     $                 pdtphys,    !delt
     $                 paprs(1,1), !ps
     $                 pplay,      !pmid
     $                 pdel,       !pdel
     $                 airephy,
     $                 pctsrf(1,1),!oro
     $                 ftsol,      !tsurf
     $                 albsol,     !albs
     $                 pphi,       !zma
     $                 pphis,      !phis
     $                 cldfra,     !cldfr
     $                 rneb,       !cldfr_st
     $                 diafra,     !cldfr_cv
     $                 itop_con,   !cldtop
     $                 ibas_con,   !cldbot
     $                 cldliq,     !cwat
     $                 prfl,       !flxrst
     $                 pmflxr,     !flxrcv
     $                 psfl,       !flxsst
     $                 pmflxs,     !flxscv
     $                 pmfu,       !flxupd
     $                 flxmass_w,  !flxmass_w
     $                 t_seri,     !tfld
     $                 sh,         !sh
     $                 rh,         !rh
     $                 .false.,    !wrhstts
     $                 hbuf,       !hbuf
     $                 obuf,       !obuf
     $                 iip1,       !nx
     $                 jjp1,       !ny
     $                 source,
     $                 solsym)
#ifdef INCAINFO
#ifdef INCA_AER

c Budget calculation for aerosol species
CALL tbudget(airephy,pdtphys,nstep,tr_seri,.false.)

c-- summary info----------------------------------------------------------------

if (MOD(nstep,nint(86400./pdtphys)) .eq. 0) then
print *, "global aerosol optical thickness "

write (form,'(A,I2,A)') "(A,",trnx-trmx+1,"A10)"
print form,"lamba [nm] ", (solsym(it),it=trmx,trnx)

write (form,'(A,I2,A)') "(I11,",trnx-trmx+1,"F10.4)"
do i=1,las
print form,int(lambda(i)),(sum(tausum(:,i,it)*airephy)/sum(airephy),it=trmx,trnx)
enddo

print *,"global mean angstroem component ", sum(angst*airephy)/sum(airephy)
endif
#endif
#endif

#ifdef INCAINFO
      PRINT *, 'OK.'
      DO it = 1, nbtr
      PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
     $                     MAXVAL(tr_seri(:,:,it))
      END DO
#endif
      call VTe(VTinca)
      call VTb(VTphysiq)
#else

c Abder
ctestmaf      if(nqmax.gt.2) aerosol(3)=.true.

       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
c Abder
#endif
c======================================================================
c   Calcul de l'effet de la convection
c======================================================================
c     print*,'Avant convection'
c      do it=1,nqmax
c         WRITE(itn,'(i2)') it
c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
c      enddo

      if (convection) then

c      print*,'Pas de temps dans phytrac : ',pdtphys
      DO it=1, nqmax
#ifdef INCA
      IF ( conv_flg(it) == 0 ) CYCLE
#endif
      if (iflag_con.eq.2) then
c tiedke
      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
      else if (iflag_con.eq.3) then
c KE
      call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
     .           upwd,dnwd,d_tr_cv(1,1,it))
      endif

      DO k = 1, nlev
      DO i = 1, klon
         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
      ENDDO
      ENDDO
#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 des thermiques
c======================================================================

      do k=1,klev
         do i=1,klon
            zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
         enddo
      enddo

c      print*,'masse dans ph ',zmasse
      do it=1,nqmax
         do k=1,klev
            do i=1,klon
               d_tr_th(i,k,it)=0.
               tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
               tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
            enddo
         enddo
      enddo

      if (iflag_thermals.gt.0) then
c        print*,'calcul de leffet des thermiques'
        nsplit=10
        DO it=1, nqmax
c        WRITE(itn,'(i1)') it
c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
c	     print*,'avant dqthermiquesretro'
c             call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI      ')

         do isplit=1,nsplit
c  Abderr 25 11 02
C Thermiques 
c	print*,'Avant dans phytrac'
            call dqthermcell(klon,klev,pdtphys/nsplit
     .       ,fm_therm,entr_therm,zmasse
     .       ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)

            do k=1,klev
               do i=1,klon
                  d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
                  d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
                  tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
               enddo
            enddo
          enddo ! nsplit
c          print*,'apres thermiques'
c             call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th      ')
c            do k=1,klev
c	print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
c	   enddo

c        WRITE(itn,'(i1)') it
c        CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
       ENDDO ! it
       endif ! Thermiques
c       print*,'ATTENTION: sdans thermniques'
      
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(1,1,it),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,it)
            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
#ifndef INCA

Cmaf provisoire source / traceur a creer
        DO i=1, klon
          source(i,it) = 0.0 ! pas de source, pour l'instant
        ENDDO
C
#endif
          CALL cltrac(pdtphys, coefh,t_seri,
     s               tr_seri(1,1,it), source(:,it),
     e               paprs, pplay, delp,
     s               d_tr_cl(1,1,it))
            DO k = 1, nlev
               DO i = 1, klon
                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
               ENDDO
            ENDDO
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
c       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======================================================================
c   Calcul de l'effet de la precipitation
c======================================================================

c      print*,'LESSIVAGE =',lessivage
      IF (lessivage) THEN

c     print*,'avant lessivage'

          d_tr_lessi_nucl(:,:,:) = 0. 
          d_tr_lessi_impa(:,:,:) = 0. 
          flestottr(:,:,:) = 0. 
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        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

c=============================================================
c   Ecriture des sorties
c=============================================================

#ifdef CPP_IOIPSL
#include "write_histrac.h"
#endif

c=============================================================

      if (lafin) then
         print*, 'c est la fin de la physique'
	 call Gather(trs(:,1),trs_tmp)
c$OMP MASTER	 
	 if (is_mpi_root) then
         
	   open (99,file='restarttrac',  form='formatted')
           do i=1,klon_glo
               write(99,*) trs_tmp(i)
           enddo
           PRINT*, 'Ecriture du fichier restarttrac'
           close(99)
	 endif
c$OMP END MASTER
      else
c         print*, 'physique pas fini'
      endif


      RETURN
      END
