Index: LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/phystokenc.F90
===================================================================
--- LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/phystokenc.F	(revision 1436)
+++ LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/phystokenc.F90	(revision 1447)
@@ -1,435 +1,362 @@
-!
-c
-c
-      SUBROUTINE phystokenc (
-     I                   nlon,nlev,pdtphys,rlon,rlat,
-     I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
-     I                   pfm_therm,pentr_therm,
-     I                   cdragh, pcoefh,yu1,yv1,ftsol,pctsrf,
-     I                   frac_impa,frac_nucl,
-     I                   pphis,paire,dtime,itap)
-      USE ioipsl
-      USE dimphy
-      USE infotrac, ONLY : nqtot
-      USE iophy
-      USE control_mod
-
-      IMPLICIT none
-
-c======================================================================
-c Auteur(s) FH
-c Objet: Moniteur general des tendances traceurs
-c
-
-c======================================================================
-#include "dimensions.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
-      integer physid, itap
-      save physid
-c$OMP THREADPRIVATE(physid)
-      integer ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
-
-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
-      real pt(klon,klev)
-      REAL,allocatable,save :: t(:,:)
-c$OMP THREADPRIVATE(t)
-c
-      REAL rlon(klon), rlat(klon), dtime
-      REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
-
-c   Couche limite:
-c   --------------
-c
-      REAL cdragh(klon)          ! cdrag
-      REAL pcoefh(klon,klev)     ! coeff melange CL
-      REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
-      REAL yv1(klon)
-      REAL yu1(klon),pphis(klon),paire(klon)
-
-c   Les Thermiques : (Abderr 25 11 02)
-c   ---------------
-      REAL pfm_therm(klon,klev+1)
-      real fm_therm1(klon,klev)
-      REAL pentr_therm(klon,klev)
-    
-      REAL,allocatable,save :: entr_therm(:,:)
-      REAL,allocatable,save :: fm_therm(:,:)
-c$OMP THREADPRIVATE(entr_therm)
-c$OMP THREADPRIVATE(fm_therm)
-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,allocatable,save :: mfu(:,:)  ! flux de masse dans le panache montant
-      REAL,allocatable,save :: mfd(:,:)  ! flux de masse dans le panache descendant
-      REAL,allocatable,save :: en_u(:,:) ! flux entraine dans le panache montant
-      REAL,allocatable,save :: de_u(:,:) ! flux detraine dans le panache montant
-      REAL,allocatable,save :: en_d(:,:) ! flux entraine dans le panache descendant
-      REAL,allocatable,save :: de_d(:,:) ! flux detraine dans le panache descendant
-      REAL,allocatable,save :: coefh(:,:) ! flux detraine dans le panache descendant
-
-      REAL,allocatable,save :: pyu1(:)
-      REAL,allocatable,save :: pyv1(:)
-      REAL,allocatable,save :: pftsol(:,:)
-      REAL,allocatable,save :: ppsrf(:,:)
-c$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
-c$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
-      real pftsol1(klon),pftsol2(klon),pftsol3(klon),pftsol4(klon)
-      real ppsrf1(klon),ppsrf2(klon),ppsrf3(klon),ppsrf4(klon)
-
-      REAL dtcum
-
-      integer iadvtr,irec
-      real zmin,zmax
-      logical ok_sync
- 
-      save dtcum
-      save iadvtr,irec
-c$OMP THREADPRIVATE(dtcum,iadvtr,irec)
-      data iadvtr,irec/0,1/
-      logical,save :: first=.true.
-c$OMP THREADPRIVATE(first)
-c
-c   Couche limite:
-c======================================================================
-
-c Dans le meme vecteur on recombine le drag et les coeff d'echange
-      pcoefh_buf(:,1)      = cdragh(:)
-      pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
-
-      ok_sync = .true.
-	print*,'Dans phystokenc.F'
-      print*,'iadvtr= ',iadvtr
-      print*,'istphy= ',istphy
-      print*,'istdyn= ',istdyn
-
-      if (first) then
-      
-        allocate( t(klon,klev))
-        allocate( mfu(klon,klev))  
-        allocate( mfd(klon,klev))  
-        allocate( en_u(klon,klev)) 
-        allocate( de_u(klon,klev)) 
-        allocate( en_d(klon,klev)) 
-        allocate( de_d(klon,klev)) 
-        allocate( coefh(klon,klev)) 
-        allocate( entr_therm(klon,klev))
-        allocate( fm_therm(klon,klev))
-        allocate( pyu1(klon))
-        allocate( pyv1(klon))
-        allocate( pftsol(klon,nbsrf))
-        allocate( ppsrf(klon,nbsrf))
-  
-        first=.false.
-      endif
-      
-      IF (iadvtr.eq.0) THEN
-	
-	CALL initphysto('phystoke',
-     . rlon,rlat,dtime, dtime*istphy,dtime*istphy,nqtot,physid)
-  	
-	write(*,*) 'apres initphysto ds phystokenc' 
-
-	
-      ENDIF
-c
-      ndex2d = 0
-      ndex3d = 0
-      i=itap 
-cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
-      CALL histwrite_phy(physid,"phis",i,pphis)
-c
-      i=itap
-cym      CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
-      CALL histwrite_phy(physid,"aire",i,paire)
-
-      iadvtr=iadvtr+1
-c
-      if (mod(iadvtr,istphy).eq.1.or.istphy.eq.1) then
-	print*,'reinitialisation des champs cumules 
-     s          a iadvtr=',iadvtr
-         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.
-                t(i,k)=0.
-		fm_therm(i,k)=0.
-               entr_therm(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_buf(i,k)*pdtphys
-                t(i,k)=t(i,k)+pt(i,k)*pdtphys
-       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
-       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(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
-
-      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
-c Unitel a enlever
-	      t(i,k)=t(i,k)/dtcum	
-               fm_therm(i,k)=fm_therm(i,k)/dtcum
-	       entr_therm(i,k)=entr_therm(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
-               pftsol1(i) = pftsol(i,1)
-               pftsol2(i) = pftsol(i,2)
-               pftsol3(i) = pftsol(i,3)
-               pftsol4(i) = pftsol(i,4)
-
-               ppsrf(i,k)=ppsrf(i,k)/dtcum
-               ppsrf1(i) = ppsrf(i,1)
-               ppsrf2(i) = ppsrf(i,2)
-               ppsrf3(i) = ppsrf(i,3)
-               ppsrf4(i) = ppsrf(i,4)
-
-            enddo
-         enddo
-c
-c   ecriture des champs
-c
-         irec=irec+1
-
-ccccc
-cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
-         CALL histwrite_phy(physid,"t",itap,t)
-
-cym         CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
-      CALL histwrite_phy(physid,"mfu",itap,mfu)
-cym	CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfd, zx_tmp_3d)
-      CALL histwrite_phy(physid,"mfd",itap,mfd)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_u, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_u",itap,en_u)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_u, zx_tmp_3d)
-      CALL histwrite_phy(physid,"de_u",itap,de_u)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, en_d, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_d",itap,en_d)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, de_d, zx_tmp_3d)       
-      CALL histwrite_phy(physid,"de_d",itap,de_d)
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1, coefh, zx_tmp_3d)         
-      CALL histwrite_phy(physid,"coefh",itap,coefh)	
-
-c ajou...
-	do k=1,klev
-           do i=1,klon
-	 fm_therm1(i,k)=fm_therm(i,k)	
-	   enddo
-	enddo
-
-cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, fm_therm1, zx_tmp_3d)
-      CALL histwrite_phy(physid,"fm_th",itap,fm_therm1)
-c
-cym      CALL gr_fi_ecrit(klev,klon,iim,jjm+1, entr_therm, zx_tmp_3d)
-      CALL histwrite_phy(physid,"en_th",itap,entr_therm)
-cccc
-cym       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_impa,zx_tmp_3d)
-        CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
-
-cym        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,frac_nucl,zx_tmp_3d)
-        CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
- 
-cym        CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyu1,zx_tmp_2d)
-      CALL histwrite_phy(physid,"pyu1",itap,pyu1)
-	
-cym	CALL gr_fi_ecrit(1, klon,iim,jjm+1, pyv1,zx_tmp_2d)
-      CALL histwrite_phy(physid,"pyv1",itap,pyv1)
-	
-cym	CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol1, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol1",itap,pftsol1)
-cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol2, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol2",itap,pftsol2)
-cym          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol3, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol3",itap,pftsol3)
-cym         CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
-      CALL histwrite_phy(physid,"ftsol4",itap,pftsol4)
-
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf1, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf1",itap,ppsrf1)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf2, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf2",itap,ppsrf2)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf3, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf3",itap,ppsrf3)
-cym        CALL gr_fi_ecrit(1,klon,iim,jjm+1, ppsrf4, zx_tmp_2d)
-      CALL histwrite_phy(physid,"psrf4",itap,ppsrf4)
-
-c$OMP MASTER
-      if (ok_sync) call histsync(physid)
-c$OMP END MASTER
-c     if (ok_sync) call histsync
-	
-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
-	go to 768
-      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.
-	       t(i,k)=0.
-               fm_therm(i,k)=0.
-	       entr_therm(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_buf(i,k)*pdtphys
-		t(i,k)=t(i,k)+pt(i,k)*pdtphys
-       fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
-       entr_therm(i,k)=entr_therm(i,k)+pentr_therm(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
-768   continue
-
-      RETURN
-      END
+SUBROUTINE phystokenc (nlon,nlev,pdtphys,rlon,rlat, &
+     pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
+     pfm_therm,pentr_therm, &
+     cdragh, pcoefh,yu1,yv1,ftsol,pctsrf, &
+     frac_impa,frac_nucl, &
+     pphis,paire,dtime,itap, &
+     psh, pda, pphi, pmp, pupwd, pdnwd)
+  
+  USE ioipsl
+  USE dimphy
+  USE infotrac, ONLY : nqtot
+  USE iophy
+  USE control_mod
+  
+  IMPLICIT NONE
+  
+!======================================================================
+! Auteur(s) FH
+! Objet: Ecriture des variables pour transport offline
+!
+!======================================================================
+  INCLUDE "dimensions.h"
+  INCLUDE "tracstoke.h"
+  INCLUDE "indicesol.h"
+  INCLUDE "iniprint.h"
+!======================================================================
+
+! Arguments:
+!
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: psh   ! humidite specifique
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pda
+  REAL,DIMENSION(klon,klev,klev), INTENT(IN):: pphi
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pmp
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pupwd ! saturated updraft mass flux
+  REAL,DIMENSION(klon,klev), INTENT(IN)     :: pdnwd ! saturated downdraft mass flux
+
+!   EN ENTREE:
+!   ==========
+!
+!   divers:
+!   -------
+!
+  INTEGER nlon ! nombre de points horizontaux
+  INTEGER nlev ! nombre de couches verticales
+  REAL pdtphys ! pas d'integration pour la physique (seconde)
+  INTEGER itap
+  INTEGER, SAVE :: physid
+!$OMP THREADPRIVATE(physid)
+
+!   convection:
+!   -----------
+!
+  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
+  REAL pt(klon,klev)
+  REAL,ALLOCATABLE,SAVE :: t(:,:)
+!$OMP THREADPRIVATE(t)
+!
+  REAL rlon(klon), rlat(klon), dtime
+  REAL zx_tmp_3d(iim,jjm+1,klev),zx_tmp_2d(iim,jjm+1)
+
+!   Couche limite:
+!   --------------
+!
+  REAL cdragh(klon)          ! cdrag
+  REAL pcoefh(klon,klev)     ! coeff melange CL
+  REAL pcoefh_buf(klon,klev) ! coeff melange CL + cdrag
+  REAL yv1(klon)
+  REAL yu1(klon),pphis(klon),paire(klon)
+
+!   Les Thermiques : (Abderr 25 11 02)
+!   ---------------
+  REAL, INTENT(IN) ::  pfm_therm(klon,klev+1)
+  REAL pentr_therm(klon,klev)
+  
+  REAL,ALLOCATABLE,SAVE :: entr_therm(:,:)
+  REAL,ALLOCATABLE,SAVE :: fm_therm(:,:)
+!$OMP THREADPRIVATE(entr_therm)
+!$OMP THREADPRIVATE(fm_therm)
+!
+!   Lessivage:
+!   ----------
+!
+  REAL frac_impa(klon,klev)
+  REAL frac_nucl(klon,klev)
+!
+! Arguments necessaires pour les sources et puits de traceur
+!
+  REAL ftsol(klon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
+  REAL pctsrf(klon,nbsrf) ! Pourcentage de sol f(nature du sol)
+!======================================================================
+!
+  INTEGER i, k, kk
+  REAL,ALLOCATABLE,SAVE :: mfu(:,:)  ! flux de masse dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: mfd(:,:)  ! flux de masse dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: en_u(:,:) ! flux entraine dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: de_u(:,:) ! flux detraine dans le panache montant
+  REAL,ALLOCATABLE,SAVE :: en_d(:,:) ! flux entraine dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: de_d(:,:) ! flux detraine dans le panache descendant
+  REAL,ALLOCATABLE,SAVE :: coefh(:,:) ! flux detraine dans le panache descendant
+  
+  REAL,ALLOCATABLE,SAVE :: pyu1(:)
+  REAL,ALLOCATABLE,SAVE :: pyv1(:)
+  REAL,ALLOCATABLE,SAVE :: pftsol(:,:)
+  REAL,ALLOCATABLE,SAVE :: ppsrf(:,:)
+!$OMP THREADPRIVATE(mfu,mfd,en_u,de_u,en_d,de_d,coefh)
+!$OMP THREADPRIVATE(pyu1,pyv1,pftsol,ppsrf)
+
+
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: sh  
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: da
+  REAL,DIMENSION(:,:,:), ALLOCATABLE,SAVE   :: phi
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: mp
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: upwd
+  REAL,DIMENSION(:,:), ALLOCATABLE,SAVE     :: dnwd
+  
+  REAL, SAVE :: dtcum
+  INTEGER, SAVE:: iadvtr=0
+!$OMP THREADPRIVATE(dtcum,iadvtr)
+  REAL zmin,zmax
+  LOGICAL ok_sync
+  CHARACTER(len=12) :: nvar
+!
+!======================================================================
+
+  iadvtr=iadvtr+1
+
+! Dans le meme vecteur on recombine le drag et les coeff d'echange
+  pcoefh_buf(:,1)      = cdragh(:)
+  pcoefh_buf(:,2:klev) = pcoefh(:,2:klev)
+  
+  ok_sync = .TRUE.
+
+! Initialization done only once
+!======================================================================
+  IF (iadvtr==1) THEN
+     ALLOCATE( t(klon,klev))
+     ALLOCATE( mfu(klon,klev))  
+     ALLOCATE( mfd(klon,klev))  
+     ALLOCATE( en_u(klon,klev)) 
+     ALLOCATE( de_u(klon,klev)) 
+     ALLOCATE( en_d(klon,klev)) 
+     ALLOCATE( de_d(klon,klev)) 
+     ALLOCATE( coefh(klon,klev)) 
+     ALLOCATE( entr_therm(klon,klev))
+     ALLOCATE( fm_therm(klon,klev))
+     ALLOCATE( pyu1(klon))
+     ALLOCATE( pyv1(klon))
+     ALLOCATE( pftsol(klon,nbsrf))
+     ALLOCATE( ppsrf(klon,nbsrf))
+     
+     ALLOCATE(sh(klon,klev)) 
+     ALLOCATE(da(klon,klev)) 
+     ALLOCATE(phi(klon,klev,klev)) 
+     ALLOCATE(mp(klon,klev)) 
+     ALLOCATE(upwd(klon,klev)) 
+     ALLOCATE(dnwd(klon,klev)) 
+
+     CALL initphysto('phystoke', dtime, dtime*istphy,dtime*istphy,physid)
+     
+     ! Write field phis and aire only once
+     CALL histwrite_phy(physid,"phis",itap,pphis)
+     CALL histwrite_phy(physid,"aire",itap,paire)
+     CALL histwrite_phy(physid,"longitudes",itap,rlon)
+     CALL histwrite_phy(physid,"latitudes",itap,rlat)
+
+  END IF
+  
+  
+! Set to zero cumulating fields
+!======================================================================
+  IF (MOD(iadvtr,istphy)==1.OR.istphy==1) THEN
+     WRITE(lunout,*)'reinitialisation des champs cumules a iadvtr=',iadvtr
+     mfu(:,:)=0.
+     mfd(:,:)=0.
+     en_u(:,:)=0.
+     de_u(:,:)=0.
+     en_d(:,:)=0.
+     de_d(:,:)=0.
+     coefh(:,:)=0.
+     t(:,:)=0.
+     fm_therm(:,:)=0.
+     entr_therm(:,:)=0.
+     pyv1(:)=0.
+     pyu1(:)=0.
+     pftsol(:,:)=0.
+     ppsrf(:,:)=0.
+     sh(:,:)=0.
+     da(:,:)=0.
+     phi(:,:,:)=0.
+     mp(:,:)=0.
+     upwd(:,:)=0.
+     dnwd(:,:)=0.
+     dtcum=0.
+  ENDIF
+  
+
+! Cumulate fields at each time step
+!======================================================================
+  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_buf(i,k)*pdtphys
+        t(i,k)=t(i,k)+pt(i,k)*pdtphys
+        fm_therm(i,k)=fm_therm(i,k)+pfm_therm(i,k)*pdtphys
+        entr_therm(i,k)=entr_therm(i,k)+pentr_therm(i,k)*pdtphys
+        sh(i,k) = sh(i,k) + psh(i,k)*pdtphys
+        da(i,k) = da(i,k) + pda(i,k)*pdtphys
+        mp(i,k) = mp(i,k) + pmp(i,k)*pdtphys
+        upwd(i,k) = upwd(i,k) + pupwd(i,k)*pdtphys
+        dnwd(i,k) = dnwd(i,k) + pdnwd(i,k)*pdtphys
+     ENDDO
+  ENDDO
+
+  DO kk=1,klev
+     DO k=1,klev
+        DO i=1,klon
+           phi(i,k,kk) = phi(i,k,kk) + pphi(i,k,kk)*pdtphys
+        END DO
+     END DO
+  END DO
+
+  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
+  
+! Add time step to cumulated time
+  dtcum=dtcum+pdtphys
+  
+
+! Write fields to file, if it is time to do so
+!======================================================================
+  IF(MOD(iadvtr,istphy)==0) THEN 
+
+     ! normalize with time period
+     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
+           t(i,k)=t(i,k)/dtcum	
+           fm_therm(i,k)=fm_therm(i,k)/dtcum
+           entr_therm(i,k)=entr_therm(i,k)/dtcum
+           sh(i,k)=sh(i,k)/dtcum
+           da(i,k)=da(i,k)/dtcum
+           mp(i,k)=mp(i,k)/dtcum
+           upwd(i,k)=upwd(i,k)/dtcum
+           dnwd(i,k)=dnwd(i,k)/dtcum
+        ENDDO
+     ENDDO
+     DO kk=1,klev
+        DO k=1,klev
+           DO i=1,klon
+              phi(i,k,kk) = phi(i,k,kk)/dtcum
+           END DO
+        END DO
+     END DO
+     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
+
+     ! write fields
+     CALL histwrite_phy(physid,"t",itap,t)
+     CALL histwrite_phy(physid,"mfu",itap,mfu)
+     CALL histwrite_phy(physid,"mfd",itap,mfd)
+     CALL histwrite_phy(physid,"en_u",itap,en_u)
+     CALL histwrite_phy(physid,"de_u",itap,de_u)
+     CALL histwrite_phy(physid,"en_d",itap,en_d)
+     CALL histwrite_phy(physid,"de_d",itap,de_d)
+     CALL histwrite_phy(physid,"coefh",itap,coefh)	
+     CALL histwrite_phy(physid,"fm_th",itap,fm_therm)
+     CALL histwrite_phy(physid,"en_th",itap,entr_therm)
+     CALL histwrite_phy(physid,"frac_impa",itap,frac_impa)
+     CALL histwrite_phy(physid,"frac_nucl",itap,frac_nucl)
+     CALL histwrite_phy(physid,"pyu1",itap,pyu1)
+     CALL histwrite_phy(physid,"pyv1",itap,pyv1)
+     CALL histwrite_phy(physid,"ftsol1",itap,pftsol(:,1))
+     CALL histwrite_phy(physid,"ftsol2",itap,pftsol(:,2))
+     CALL histwrite_phy(physid,"ftsol3",itap,pftsol(:,3))
+     CALL histwrite_phy(physid,"ftsol4",itap,pftsol(:,4))
+     CALL histwrite_phy(physid,"psrf1",itap,ppsrf(:,1))
+     CALL histwrite_phy(physid,"psrf2",itap,ppsrf(:,2))
+     CALL histwrite_phy(physid,"psrf3",itap,ppsrf(:,3))
+     CALL histwrite_phy(physid,"psrf4",itap,ppsrf(:,4))
+     CALL histwrite_phy(physid,"sh",itap,sh)
+     CALL histwrite_phy(physid,"da",itap,da)
+     CALL histwrite_phy(physid,"mp",itap,mp)
+     CALL histwrite_phy(physid,"upwd",itap,upwd)
+     CALL histwrite_phy(physid,"dnwd",itap,dnwd)
+
+
+! phi
+     DO k=1,klev
+        IF (k<10) THEN
+           WRITE(nvar,'(i1)') k
+        ELSE IF (k<100) THEN
+           WRITE(nvar,'(i2)') k
+        ELSE
+           WRITE(nvar,'(i3)') k
+        END IF
+        nvar='phi_lev'//trim(nvar)
+        
+        CALL histwrite_phy(physid,nvar,itap,phi(:,:,k))
+     END DO
+     
+     ! Syncronize file
+!$OMP MASTER
+     IF (ok_sync) CALL histsync(physid)
+!$OMP END MASTER
+     
+     
+     ! Calculate min and max values for some fields (coefficients de lessivage)
+     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
+     WRITE(lunout,*)'------ coefs de lessivage (min et max) --------'
+     WRITE(lunout,*)'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
+     WRITE(lunout,*)'facteur d impaction ',zmin,zmax
+     
+  ENDIF ! IF(MOD(iadvtr,istphy)==0)
+
+END SUBROUTINE phystokenc
