SUBROUTINE phystoke (
     I                   nlon,nlev,pdtphys,
     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
     I                   pcoefh,yu1,yv1,ftsol,pctsrf,
     I                   frac_impa,frac_nucl)
      IMPLICIT none

c======================================================================
c Auteur(s) FH
c Objet: Moniteur general des tendances traceurs
c

c======================================================================
#include "dimensions.h"
#include "dimphy.h"
#include "tracstoke.h"
#include "indicesol.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
      real pdtphys ! pas d'integration pour la physique (seconde)
c
c   convection:
c   -----------
c
      REAL pmfu(klon,klev)  ! flux de masse dans le panache montant
      REAL pmfd(klon,klev)  ! flux de masse dans le panache descendant
      REAL pen_u(klon,klev) ! flux entraine dans le panache montant
      REAL pde_u(klon,klev) ! flux detraine dans le panache montant
      REAL pen_d(klon,klev) ! flux entraine dans le panache descendant
      REAL pde_d(klon,klev) ! flux detraine dans le panache descendant
c
c   Couche limite:
c   --------------
c
      REAL pcoefh(klon,klev)    ! coeff melange CL
      REAL yv1(klon)
      REAL yu1(klon)
c
c   Lessivage:
c   ----------
c
      REAL frac_impa(klon,klev)
      REAL frac_nucl(klon,klev)
c
c Arguments necessaires pour les sources et puits de traceur
C
      real ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
      real pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)

c======================================================================
c
      INTEGER i, k
c
      REAL mfu(klon,klev)  ! flux de masse dans le panache montant
      REAL mfd(klon,klev)  ! flux de masse dans le panache descendant
      REAL en_u(klon,klev) ! flux entraine dans le panache montant
      REAL de_u(klon,klev) ! flux detraine dans le panache montant
      REAL en_d(klon,klev) ! flux entraine dans le panache descendant
      REAL de_d(klon,klev) ! flux detraine dans le panache descendant
      REAL coefh(klon,klev) ! flux detraine dans le panache descendant

      REAL pyu1(klon),pyv1(klon)
      REAL pftsol(klon,nbsrf),ppsrf(klon,nbsrf)

      REAL dtcum

      integer iadvtr,irec
      real zmin,zmax

      save mfu,mfd,en_u,de_u,en_d,de_d,coefh,dtcum
      save iadvtr,irec
      save pyu1,pyv1,pftsol,ppsrf

      data iadvtr,irec/0,1/
c
c   Couche limite:
c======================================================================

      print*,'iadvtr= ',iadvtr
      print*,'istphy= ',istphy
      print*,'istdyn= ',istdyn

      IF (iadvtr.eq.0) THEN
#ifdef CRAY
         CALL ASSIGN("assign -N ieee -F null f:physique")
#endif
         open(49,file='physique',form='unformatted',
     s        access='direct',recl=4*klon*(9*klev+2+2*nbsrf))
      ENDIF
c
      iadvtr=iadvtr+1
c
      IF(mod(iadvtr,istphy).eq.0) THEN 
c
c   normalisation par le temps cumule
         do k=1,klev
            do i=1,klon
               mfu(i,k)=mfu(i,k)/dtcum
               mfd(i,k)=mfd(i,k)/dtcum
               en_u(i,k)=en_u(i,k)/dtcum
               de_u(i,k)=de_u(i,k)/dtcum
               en_d(i,k)=en_d(i,k)/dtcum
               de_d(i,k)=de_d(i,k)/dtcum
               coefh(i,k)=coefh(i,k)/dtcum
            enddo
         enddo
         do i=1,klon
            pyv1(i)=pyv1(i)/dtcum
            pyu1(i)=pyu1(i)/dtcum
         end do
         do k=1,nbsrf
             do i=1,klon
               pftsol(i,k)=pftsol(i,k)/dtcum
               ppsrf(i,k)=ppsrf(i,k)/dtcum
            enddo
         enddo
c
c   ecriture des champs
c
         irec=irec+1
         write(49,rec=1) float(irec),float(istphy),
     s    float(klon),float(klev)
         write(49,rec=irec) mfu,mfd,en_u,de_u,en_d,de_d
     s        ,coefh
     s        ,frac_impa,frac_nucl
     s        ,pyu1,pyv1,pftsol,ppsrf
c
cAA Test sur la valeur des coefficients de lessivage 
c
         zmin=1e33
         zmax=-1e33
         do k=1,klev
            do i=1,klon
                  zmax=max(zmax,frac_nucl(i,k))
                  zmin=min(zmin,frac_nucl(i,k))
            enddo
         enddo
         Print*,'------ coefs de lessivage (min et max) --------'
         Print*,'facteur de nucleation ',zmin,zmax
         zmin=1e33
         zmax=-1e33
         do k=1,klev
            do i=1,klon
                  zmax=max(zmax,frac_impa(i,k))
                  zmin=min(zmin,frac_impa(i,k))
            enddo
         enddo
         Print*,'facteur d impaction ',zmin,zmax

      ENDIF 

c   reinitialisation des champs cumules
      if (mod(iadvtr,istphy).eq.1) then
         do k=1,klev
            do i=1,klon
               mfu(i,k)=0.
               mfd(i,k)=0.
               en_u(i,k)=0.
               de_u(i,k)=0.
               en_d(i,k)=0.
               de_d(i,k)=0.
               coefh(i,k)=0.
            enddo
         enddo
         do i=1,klon
            pyv1(i)=0.
            pyu1(i)=0.
         end do
         do k=1,nbsrf
             do i=1,klon
               pftsol(i,k)=0.
               ppsrf(i,k)=0.
            enddo
         enddo

         dtcum=0.
      endif

      do k=1,klev
         do i=1,klon
            mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
            mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
            en_u(i,k)=en_u(i,k)+pen_u(i,k)*pdtphys
            de_u(i,k)=de_u(i,k)+pde_u(i,k)*pdtphys
            en_d(i,k)=en_d(i,k)+pen_d(i,k)*pdtphys
            de_d(i,k)=de_d(i,k)+pde_d(i,k)*pdtphys
            coefh(i,k)=coefh(i,k)+pcoefh(i,k)*pdtphys
         enddo
      enddo
         do i=1,klon
            pyv1(i)=pyv1(i)+yv1(i)*pdtphys
            pyu1(i)=pyu1(i)+yu1(i)*pdtphys
         end do
         do k=1,nbsrf
             do i=1,klon
               pftsol(i,k)=pftsol(i,k)+ftsol(i,k)*pdtphys
               ppsrf(i,k)=ppsrf(i,k)+pctsrf(i,k)*pdtphys
            enddo
         enddo

      dtcum=dtcum+pdtphys

      RETURN
      END