SUBROUTINE phytrac (rnpb,
     I                   debutphy,
     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)
      USE ioipsl
      USE histcom

      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 "dimphy.h"
#include "indicesol.h"
#include "control.h"
#include "temps.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
      real pdtphys  ! pas d'integration pour la physique (seconde)
      real t_seri(nlon,nlev) ! temperature
      real tr_seri(nlon,nlev,nbtr) ! traceur  
      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
      real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
      real znivsig(klev) ! niveaux sigma
      real paire(klon)
      real pphis(klon)
      logical debutphy       ! le flag de l'initialisation 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)

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*1 itn
C maf ioipsl
      CHARACTER*2 str2
      INTEGER nhori, nvert
      REAL zsto, zout, zjulian
      INTEGER nid_tra
      SAVE nid_tra
      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
      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

      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./
c
c======================================================================
c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
         modname='phytrac'

c        print*,'DANS PHYTRAC debutphy=',debutphy

         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
         zsto = pdtphys
         zout = pdtphys * FLOAT(ecrit_tra)
         if (debutphy) then

         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
         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
         itra=0
C         
         CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
         zjulian = zjulian + day_ini
c
         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, 0, zjulian, pdtphys,
     .                 nhori, nid_tra)
         call histvert(nid_tra, 'sig_s', 'Niveaux sigma','-',
     .              klev, znivsig, nvert)
C
C         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
C     .                 klev, presnivs, nvert)
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)
c
         DO it=1,nqmax
         IF (it.LE.99) THEN
         WRITE(str2,'(i2.2)') it
C champ 3D
         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
         ELSE
         PRINT*, "Trop de traceurs"
         CALL abort
         ENDIF
         ENDDO
         CALL histend(nid_tra)

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
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
c======================================================================
c   Calcul de l'effet de la convection
c======================================================================

      if (convection) then

c     print*,'Pas de temps dans phytrac : ',pdtphys
      DO it=1, nqmax
      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
      WRITE(itn,'(i1)') it
      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
      ENDDO
c     print*,'apres nflxtr'


      endif ! convection

c======================================================================
c   Calcul de l'effet de la couche limite
c======================================================================

c     print*,'avant couchelimite'
      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
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
Cmaf provisoire source / traceur a creer
        DO i=1, klon
          source(i) = 0.0 ! pas de source, pour l'instant
        ENDDO
C
          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======================================================================
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
        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======================================================================

      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
       DO it = 1, nqmax
         IF (aerosol(it)) THEN
           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) - 1. )    
              ENDDO
           ENDDO
         ENDIF
       ENDDO
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

C
C Sorties IOIPSL
      ndex2d = 0
      ndex3d = 0
c
c     write(*,*)'sorties ioipsl phytrac',zsto,zout
      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
      CALL histwrite(nid_tra,"phis",itra,zx_tmp_2d,iim*(jjm+1),ndex2d)
C
      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),ndex2d)
      DO it=1,nqmax
      IF (it.LE.99) THEN
       WRITE(str2,'(i2.2)') it
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,ndex3d)
       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,ndex3d)
      ENDIF
      ELSE
         PRINT*, "Trop de traceurs"
         CALL abort
      ENDIF
      ENDDO
      if (ok_sync) call histsync(nid_tra)

      RETURN
      END