Index: /LMDZ.3.3/trunk/libf/dyn3d/fluxstokenc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/fluxstokenc.F	(revision 55)
+++ /LMDZ.3.3/trunk/libf/dyn3d/fluxstokenc.F	(revision 55)
@@ -0,0 +1,172 @@
+      SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
+     . time_step,itau, fluxid, fluxvid,fluxdid )
+
+       USE IOIPSL
+c
+c     Auteur :  F. Hourdin
+c
+c
+ccc   ..   Modif. P. Le Van  ( 20/12/97 )  ...
+c
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "tracstoke.h"
+#include "temps.h"
+
+      REAL time_step,t_wrt, t_ops
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL masse(ip1jmp1,llm),teta(ip1jmp1,llm),phi(ip1jmp1,llm)
+      REAL phis(ip1jmp1)
+
+      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+      REAL massem(ip1jmp1,llm),tetac(ip1jmp1,llm),phic(ip1jmp1,llm)
+
+      REAL pbarug(ip1jmp1,llm),pbarvg(iip1,jjm,llm),wg(ip1jmp1,llm)
+
+      REAL pbarvst(iip1,jjp1,llm),zistdyn
+
+
+      INTEGER iadvtr
+      integer ndex1d(1)
+      integer ndex2d(ip1jmp1)
+      integer ndex3dv(ip1jm*llm),ndex3d(ip1jmp1*llm)
+      integer nscal
+      real tst(1),ist(1),istp(1)
+      INTEGER ij,l,irec,i,j,itau
+      INTEGER fluxid, fluxvid,fluxdid
+ 
+      SAVE iadvtr, massem,pbaruc,pbarvc,irec
+      SAVE phic,tetac
+      logical first
+      save first
+      data first/.true./
+      DATA iadvtr/0/
+
+      if(first) then
+
+	CALL initfluxsto( 'fluxstoke',
+     .  time_step,istdyn* time_step,istdyn* time_step,
+     . nqmx, fluxid,fluxvid,fluxdid) 
+	first = .false.
+
+      endif
+
+
+      ndex1d = 0
+      ndex2d = 0
+      ndex3dv = 0
+      ndex3d = 0
+
+
+      IF(iadvtr.EQ.0) THEN
+         CALL initial0(ijp1llm,phic)
+         CALL initial0(ijp1llm,tetac)
+         CALL initial0(ijp1llm,pbaruc)
+         CALL initial0(ijmllm,pbarvc)
+      ENDIF
+
+c   accumulation des flux de masse horizontaux
+      DO l=1,llm
+         DO ij = 1,ip1jmp1
+            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
+            tetac(ij,l) = tetac(ij,l) + teta(ij,l)
+            phic(ij,l) = phic(ij,l) + phi(ij,l)
+         ENDDO
+         DO ij = 1,ip1jm
+            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
+         ENDDO
+      ENDDO
+
+c   selection de la masse instantannee des mailles avant le transport.
+      IF(iadvtr.EQ.0) THEN
+         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
+      ENDIF
+
+      iadvtr   = iadvtr+1
+
+
+c   Test pour savoir si on advecte a ce pas de temps
+      IF ( iadvtr.EQ.istdyn ) THEN
+
+c    normalisation
+      DO l=1,llm
+         DO ij = 1,ip1jmp1
+            pbaruc(ij,l) = pbaruc(ij,l)/float(istdyn)
+            tetac(ij,l) = tetac(ij,l)/float(istdyn)
+            phic(ij,l) = phic(ij,l)/float(istdyn)
+         ENDDO
+         DO ij = 1,ip1jm
+            pbarvc(ij,l) = pbarvc(ij,l)/float(istdyn)
+         ENDDO
+      ENDDO
+
+c   traitement des flux de masse avant advection.
+c     1. calcul de w
+c     2. groupement des mailles pres du pole.
+
+        CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
+
+        do l=1,llm
+           do j=1,jjm
+              do i=1,iip1
+                 pbarvst(i,j,l)=pbarvg(i,j,l)
+              enddo
+           enddo
+           do i=1,iip1
+              pbarvst(i,jjp1,l)=0.
+           enddo
+        enddo
+
+         iadvtr=0
+
+c     write(*,*)'histwrite phis'
+      call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex2d)
+c     write(*,*)'histwrite aire'
+      call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex2d)
+	
+      nscal = 1
+      tst(1) = time_step
+c     write(*,*)'histwrite dtvr'
+      call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex1d)
+      ist(1)=istdyn
+c     write(*,*)'histwrite istdyn'
+      call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex1d)
+      istp(1)= istphy
+c     write(*,*)'histwrite istphy'
+      call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex1d)
+
+c       write(*,*)'histwrite masse'	
+	call histwrite(fluxid, 'masse', itau, masse,
+     .               iip1*jjp1*llm, ndex3d)
+	
+c       write(*,*)'histwrite pbaru'	
+	call histwrite(fluxid, 'pbaru', itau, pbarug,
+     .               iip1*jjp1*llm, ndex3d)
+	
+c       write(*,*)'histwrite pbarv'	
+	call histwrite(fluxvid, 'pbarv', itau, pbarvst,
+     .               iim*jjp1*llm, ndex3dv)
+	
+c       write(*,*)'histwrite w'	
+        call histwrite(fluxid, 'w' ,itau, wg, 
+     .             iip1*jjp1*llm, ndex3d) 
+	
+c       write(*,*)'histwrite teta'	
+	call histwrite(fluxid, 'teta' ,itau, tetac, 
+     .             iip1*jjp1*llm, ndex3d) 
+	
+c       write(*,*)'histwrite phi'	
+	call histwrite(fluxid, 'phi' ,itau, phic, 
+     .             iip1*jjp1*llm, ndex3d) 
+	
+C
+
+      ENDIF ! if iadvtr.EQ.istdyn
+
+      RETURN
+      END