Index: DZ.3.3/trunk/libf/dyn3d/fluxstoke.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/fluxstoke.F	(revision 203)
+++ 	(revision )
@@ -1,149 +1,0 @@
-      SUBROUTINE fluxstoke(pbaru,pbarv,masse,teta,phi,phis)
-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"
-
-
-      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)
-
-
-      INTEGER iadvtr 
-      INTEGER ij,l,irec,i,j
- 
-      SAVE iadvtr, massem,pbaruc,pbarvc,irec
-      SAVE phic,tetac
-      logical first
-      save first
-      data first/.true./
-      DATA iadvtr/0/
-
-      if(first) then
-#ifdef CRAY
-         CALL ASSIGN("assign -N ieee -F null f:fluxmass")
-#endif
-         open(47,file='fluxmass',form='unformatted',
-     s        access='direct',recl=4*(6*ijp1llm))
-         irec=1
-         first=.false.
-
-         open(77,file='fluxmass.ctl',status='unknown',form='formatted')
-
-      endif
-
-
-      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
-
-         irec=irec+1
-         write(47,rec=1) float(irec),dtvr,float(istdyn),
-     s    float(iim),float(jjm),float(llm),rlonu,rlonv,rlatu,rlatv
-     s    ,aire,phis
-         write(47,rec=irec) massem,pbarug,pbarvst,wg,tetac,phic
-
-c   on reinitialise a zero les flux de masse cumules.
-
-      write(77,'(a4,2x,a40)')
-     &       'DSET ','^fluxmass'
-
-      write(77,'(a12)') 'UNDEF 1.0E30'
-      write(77,'(a5,1x,a40)') 'TITLE ','Titre a voir'
-      call formcoord(77,iip1,rlonv,180./pi,.false.,'XDEF')
-      call formcoord(77,jjp1,rlatu,180./pi,.true.,'YDEF')
-      call formcoord(77,llm,presnivs,1.,.false.,'ZDEF')
-      write(77,'(a4,i10,a30)')
-     &       'TDEF ',irec,' LINEAR 02JAN1987 1DY '
-      write(77,'(a4,2x,i5)') 'VARS',6
-      write(77,1000) 'masse',llm,99,'masse    '
-      write(77,1000) 'pbaru',llm,99,'pbaru    '
-      write(77,1000) 'pbarv',llm,99,'pbarv    '
-      write(77,1000) 'w    ',llm,99,'w        '
-      write(77,1000) 'teta ',llm,99,'teta     '
-      write(77,1000) 'phi  ',llm,99,'phi      '
-      write(77,'(a7)') 'ENDVARS'
-
-1000  format(a5,3x,i4,i3,1x,a39)
-
-
-
-      ENDIF ! if iadvtr.EQ.istdyn
-
-      RETURN
-      END
Index: /LMDZ.3.3/trunk/libf/dyn3d/grilles_gcm_netcdf.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/grilles_gcm_netcdf.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/dyn3d/grilles_gcm_netcdf.F	(revision 204)
@@ -1,3 +1,6 @@
+c
 c $Header
+c
+
       PROGRAM create_fausse_var
 C
@@ -9,4 +12,7 @@
 #include "comconst.h"
 #include "comgeom.h"
+#include "comvert.h"
+
+      real temp(iim+1,jjm+1)
 #include "netcdf.inc"
 
@@ -19,11 +25,9 @@
         integer*4 out_latudim,out_latvdim,out_dim(3)
 
-c   champs a faire disparaitre
-        real ucov(ijp1llm),vcov(ijmllm),teta(ijp1llm)
-        real masse(ijp1llm),ps(ip1jmp1),phis(ip1jmp1)
-        real q(ijp1llm*nqmx)
-        real time_0
-
-        real clesphy0(20),pa,preff
+      INTEGER         longcles
+      PARAMETER     ( longcles = 20 )
+      REAL  clesphy0( longcles )
+
+      integer start(4),count(4)
 
 	integer status,i,j
@@ -31,4 +35,11 @@
         real rlonudeg(iip1),rlonvdeg(iip1)
 
+      real dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
+      real acoslat,dxkm,dykm,resol(iip1,jjp1)
+
+#include "serre.h"
+#include "fxyprim.h"
+
+      print*,'OK0'
 
       rad = 6400000
@@ -42,16 +53,12 @@
       pa= 50000.
 
-c     CALL dynetat0("start.nc",nqmx,vcov,ucov,
-c    .              teta,q,masse,ps,phis, time_0)
-
-
       open(99,file='run.def',status='old',form='formatted')
-      CALL defrun_new( 99, .TRUE. ,clesphy0)
+      CALL defrun_new( 99, .TRUE.,clesphy0 )
       close(99)
       CALL iniconst
-      print*,'inigeom pas OK'
       CALL inigeom
-      print*,'inigeom OK'
-
+
+
+      print*,'OK1'
       do j=1,jjp1
          rlatudeg(j)=rlatu(j)*180./pi
@@ -66,5 +73,7 @@
       enddo
 
-      print*,'  2 ----- OUVERTURE DE LA SORTIE NETCDF'
+
+      print*,'OK2'
+c  2 ----- OUVERTURE DE LA SORTIE NETCDF
 c ---------------------------------------------------
 c CREATION OUTPUT
@@ -78,4 +87,6 @@
         status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
 
+
+      print*,'OK3'
 c   Longitudes en u
         print *,'OUTID: ',ncid_out
@@ -138,7 +149,19 @@
      %  16,'Grille aux point v')
 
-
+c   ecriture de la grille u
+        out_dim(1)=out_lonvdim
+        out_dim(2)=out_latudim
+        status=NF_DEF_VAR(ncid_out,'grille_s',NF_FLOAT,2,out_dim,
+     %  out_varid)
+        call handle_err(status)
+        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units',
+     %  6,'Kelvin')
+        status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',
+     %  16,'Grille aux point u')
+
+
+      print*,'OK4'
         status=NF_ENDDEF(ncid_out)
-        print*,' 5) ----- FERMETURE DES FICHIERS NETCDF-'
+c 5) ----- FERMETURE DES FICHIERS NETCDF------------------
 c --------------------------------------------------------
 c 3-b- Ecriture de la grille pour la sortie
@@ -150,4 +173,22 @@
         status=NF_PUT_VARA_REAL(ncid_out,out_latvid,1,jjm,rlatvdeg)
 
+        start(1)=1
+        start(2)=1
+        start(3)=1
+        start(4)=1
+
+        count(1)=iim+1
+        count(2)=jjm+1
+        count(3)=1
+        count(4)=1
+
+        do j=1,jjm+1
+           do i=1,iim+1
+              temp(i,j)=mod(i,2)+mod(j,2)
+           enddo
+        enddo
+
+        status=NF_PUT_VARA_REAL(ncid_out,out_varid,start,
+     s  count,temp)
 
 c fermeture du fichier netcdf
@@ -155,4 +196,64 @@
         write(*,*) 'Fermeture: ',fich_out
 
+
+      print*,'OK5'
+c   Ecriture grads
+      open (20,file='grille.dat',form='unformatted',access='direct'
+     s      ,recl=4*ip1jmp1)
+      write(20,rec=1) ((float(mod(i,2)+mod(j,2)),i=1,iip1),j=1,jjp1)
+      write(20,rec=2) ((float(mod(i,2)*mod(j,2)),i=1,iip1),j=1,jjp1)
+      do j=2,jjm
+         dlat1(j)=180.*(rlatv(j)-rlatv(j-1))/pi
+c        dlat2(j)=180.*fyprim(float(j))/pi
+      enddo
+      do i=2,iip1
+         dlon1(i)=180.*(rlonu(i)-rlonu(i-1))/pi
+c        dlon2(i)=180.*fxprim(float(i))/pi
+      enddo
+      do j=2,jjm
+         dykm=(rlatv(j)-rlatv(j-1))*6400.
+         acoslat=6400.*cos(rlatu(j))
+         do i=2,iip1
+            dxkm=acoslat*(rlonu(i)-rlonu(i-1))
+            resol(i,j)=sqrt(dykm*dykm+dxkm*dxkm)
+         enddo
+         resol(1,j)=resol(iip1,j)
+      enddo
+      write(20,rec=3) resol
+      dlon1(1)=dlon1(iip1)
+      dlon2(1)=dlon2(iip1)
+      write(20,rec=4) ((dlon1(i),i=1,iip1),j=1,jjp1)
+      write(20,rec=5) ((dlon2(i),i=1,iip1),j=1,jjp1)
+      write(20,rec=6) ((dlat1(j),i=1,iip1),j=1,jjp1)
+      write(20,rec=7) ((dlat2(j),i=1,iip1),j=1,jjp1)
+        
+      open (21,file='grille.ctl',form='formatted')
+
+c   WARNING! on reecrase le fichier .ctl a chaque ecriture
+      write(21,'(a5,1x,a40)')
+     &       'DSET ','^grille.dat'
+
+      write(21,'(a12)') 'UNDEF 1.0E30'
+      write(21,'(a5,1x,a40)') 'TITLE ','grille'
+      call formcoord(21,iip1,rlonv,180./pi,.false.,'XDEF')
+      call formcoord(21,jjp1,rlatu,180./pi,.true.,'YDEF')
+      call formcoord(21,1,0.,1.,.false.,'ZDEF')
+      write(21,'(a4,i10,a30)')
+     &       'TDEF ',1,' LINEAR 23OCT1994 3hr '
+      write(21,'(a4,2x,i5)') 'VARS',7
+      write(21,'(a18)') 'grille 0 99 grille'
+      write(21,'(a18)') 'gril   0 99 gril  '
+      write(21,'(a29)') 'resol   0 99 resolution (km)  '
+      write(21,'(a18)') 'dlon1  0 99 dlon1 '
+      write(21,'(a18)') 'dlon2  0 99 dlon2 '
+      write(21,'(a18)') 'dlat1  0 99 dlat1 '
+      write(21,'(a18)') 'dlat2  0 99 dlat2 '
+      write(21,'(a7)') 'ENDVARS'
+
+
+
+
+
+      print*,'OK6'
 	end
 
@@ -161,4 +262,6 @@
         subroutine handle_err(status)
 #include "netcdf.inc"
+
+
         integer status
         print *,'handle code err: ',NF_NOERR
Index: DZ.3.3/trunk/libf/dyn3d/lectflux.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/lectflux.F	(revision 203)
+++ 	(revision )
@@ -1,233 +1,0 @@
-      SUBROUTINE lectflux(irec,massem,pbarun,pbarvn,wn,tetan,phin,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf,
-     s     frac_impa,frac_nucl,phis)
-
-      IMPLICIT NONE
-
-#include "dimensions.h"
-#include "paramet.h"
-
-#include "comvert.h"
-#include "comconst.h"
-#include "comgeom2.h"
-
-#include "tracstoke.h"
-
-      integer irec,nrec,i,j
-
-      integer ngridmx,ig,l
-      parameter (ngridmx=iim*(jjm-1)+2)
-      INTEGER nbsrf
-      PARAMETER (nbsrf=4) ! nombre de sous-fractions pour une maille
-
-      real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
-      real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
-      real coefkz(ngridmx,llm)
-      real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
-      real yu1(ngridmx), yv1(ngridmx)
-      real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
-      integer imfu,imfd,ien_u,ide_u,
-     s      ien_d,ide_d,
-     s      icoefkz,izu1,izv1,
-     s      itsol,ipsf,
-     s      ilei, ilec
-      parameter(imfu=1,imfd=llm+1,ien_u=2*llm+1,ide_u=3*llm+1,
-     s      ien_d=4*llm+1,ide_d=5*llm+1,
-     s      icoefkz=6*llm+1,
-     s      ilei=7*llm+1,ilec=8*llm+1,
-     s      izu1=9*llm+1,izv1=9*llm+2,
-     s      itsol=9*llm+3,ipsf=9*llm+3+nbsrf)
-      logical avant
-
-      real massefi(ngridmx,llm)
-
-      real massem(ip1jmp1,llm),tetan(ip1jmp1,llm)
-      real pbarun(iip1,jjp1,llm),pbarvn(iip1,jjm,llm)
-      real pbarvst(iip1,jjp1,llm)
-      real wn(iip1,jjp1,llm),phin(iip1,jjp1,llm)
-      real phis(iip1,jjp1)
-
-      real airefi(ngridmx)
-
-      real xlecn(ngridmx,9*llm+2+2*nbsrf)
-
-      real zcontrole(ngridmx),zmass,tmpdyn(iip1,jjp1),zflux
-
-      real ziadvtrac,zrec,ziadvtrac2,zrec2
-      real zim,zjm,zlm,zklon,zklev
-
-      real zpi
-
-      zpi=2.*asin(1.)
-
-
-c==================================================================
-c   Si le numero du record est 0 alors: INITIALISATION
-c==================================================================
-c
-      print*,'ENTREE DANS LECTFLUX'
-        print*,'IREC=',IREC
-      if(irec.eq.0) then
-
-        print*,'IREC==',0
-
-C test         call inigeom
-
-c==================================================================
-c   ouverture des fichiers
-c==================================================================
-
-c   Fichier fluxmass
-#ifdef CRAY
-         CALL ASSIGN("assign -N ieee -F null f:fluxmass")
-#endif
-      open(47,file='fluxmass',form='unformatted',
-     s     access='direct'
-     s     ,recl=4*(6*ijp1llm))
-      read(47,rec=1) zrec,dtvr,ziadvtrac,zim,zjm,zlm,
-     s   rlonu,rlonv,rlatu,rlatv,aire
-     s    ,phis
-      print*,'zrec,dtvr,ziadvtrac,zim,zjm,zlm'
-      print*,zrec,dtvr,ziadvtrac,zim,zjm,zlm
-      print*,rlonv
-
-
-
-c  Fichier physique
-c  Fichier lessivage (supprime les donnees utiles sont dans "physique")
-#ifdef CRAY
-         CALL ASSIGN("assign -N ieee -F null f:physique")
-#endif
-      open(49,file='physique',form='unformatted',
-     s     access='direct'
-     s     ,recl=4*ngridmx*(9*llm+2+2*nbsrf))
-      read(49,rec=1) zrec2,ziadvtrac2,zklon,zklev
-      print*,'Entete du fichier physique'
-      print*,zrec2,ziadvtrac2,zklon,zklev
-
-      nrec=zrec
-      print*,'nrec=',nrec
-
-      istdyn=ziadvtrac
-      istphy=ziadvtrac2
-
-c==================================================================
-c   Fin des initialisations
-      else ! irec=0
-c==================================================================
-
-
-c-----------------------------------------------------------------------
-c   Lecture des fichiers fluxmass et  physique:
-c   -----------------------------------------------------
-
-c  Variables dynamiques
-         read(47,rec=irec) massem,pbarun,pbarvst,wn,tetan,phin
-        do l=1,llm
-           do j=1,jjm
-              do i=1,iip1
-                 pbarvn(i,j,l)=pbarvst(i,j,l)
-              enddo
-           enddo
-        enddo
-
-c  Variables physiques
-         read(49,rec=irec) ((xlecn(ig,l),ig=1,ngridmx),
-     s                                    l=1,9*llm+2+2*nbsrf)
-
-       do l=1,llm
-          do ig=1,ngridmx
-             coefkz(ig,l)=xlecn(ig,icoefkz+l-1)
-             frac_impa(ig,l)=xlecn(ig,ilei+l-1)
-             frac_nucl(ig,l)=xlecn(ig,ilec+l-1)
-          enddo
-       enddo
-       do l=1,nbsrf
-          do ig=1,ngridmx
-             ftsol(ig,l)=xlecn(ig,itsol+l-1)
-             pctsrf(ig,l)=xlecn(ig,ipsf+l-1)
-          enddo
-       enddo
-       do ig=1,ngridmx
-          yv1(ig)=xlecn(ig,izv1)
-          yu1(ig)=xlecn(ig,izu1)
-       enddo
-C
-      if(avant) then
-c   Simu directe
-       do l=1,llm
-          do ig=1,ngridmx
-             zmfu(ig,l)=xlecn(ig,imfu+l-1)
-             zmfd(ig,l)=xlecn(ig,imfd+l-1)
-             zde_u(ig,l)=xlecn(ig,ide_u+l-1)
-             zen_u(ig,l)=xlecn(ig,ien_u+l-1)
-             zde_d(ig,l)=xlecn(ig,ide_d+l-1)
-             zen_d(ig,l)=xlecn(ig,ien_d+l-1)
-          enddo
-       enddo
-      else
-c   Simu retro
-       do l=1,llm
-          do ig=1,ngridmx
-             zmfd(ig,l)=-xlecn(ig,imfu+l-1)
-             zmfu(ig,l)=-xlecn(ig,imfd+l-1)
-             zen_d(ig,l)=xlecn(ig,ide_u+l-1)
-             zde_d(ig,l)=xlecn(ig,ien_u+l-1)
-             zen_u(ig,l)=xlecn(ig,ide_d+l-1)
-             zde_u(ig,l)=xlecn(ig,ien_d+l-1)
-          enddo
-       enddo
-      endif
-
-c-----------------------------------------------------------------------
-c   PETIT CONTROLE SUR LES FLUX CONVECTIFS...
-c-----------------------------------------------------------------------
-
-      print*,'Ap redec irec'
-
-         call gr_dyn_fi(llm,iip1,jjp1,ngridmx,massem,massefi)
-
-         do ig=1,ngridmx
-            zcontrole(ig)=1.
-         enddo
-c   zmass=(max(massem(ig,l),massem(ig,l-1))/airefi(ig)
-         do l=2,llm
-            do ig=1,ngridmx
-               zmass=max(massefi(ig,l),massefi(ig,l-1))/airefi(ig)
-               zflux=max(abs(zmfu(ig,l)),abs(zmfd(ig,l)))*dtphys
-               if(zflux.gt.0.9*zmass) then
-                 zcontrole(ig)=min(zcontrole(ig),0.9*zmass/zflux)
-               endif
-            enddo
-         enddo
-
-         do ig=1,ngridmx
-            if(zcontrole(ig).lt.0.99999) then
-               print*,'ATTENTION !!! on reduit les flux de masse '
-               print*,'convectifs au point ig=',ig
-            endif
-         enddo
-
-         call gr_fi_dyn(1,ngridmx,iip1,jjp1,zcontrole,tmpdyn)
-
-         do l=1,llm
-            do ig=1,ngridmx
-               zmfu(ig,l)=zmfu(ig,l)*zcontrole(ig)
-               zmfd(ig,l)=zmfd(ig,l)*zcontrole(ig)
-               zen_u(ig,l)=zen_u(ig,l)*zcontrole(ig)
-               zde_u(ig,l)=zde_u(ig,l)*zcontrole(ig)
-               zen_d(ig,l)=zen_d(ig,l)*zcontrole(ig)
-               zde_d(ig,l)=zde_d(ig,l)*zcontrole(ig)
-            enddo
-         enddo
-
-
-      endif ! irec=0
-
-
-      RETURN
-      END
-
-
Index: /LMDZ.3.3/trunk/libf/dyn3d/nudge.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/nudge.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/dyn3d/nudge.F	(revision 204)
@@ -99,7 +99,6 @@
      s   ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
 cnec         read(*,*) online
-          online=-1
-         print*,'Entrer les constantes de temps de rappel en jours'
-         print*,'alphaT,alphau,alphav,alphaP'
+cnec      online=-1
+cnec     print*,'Entrer les constantes de temps de rappel en jours'
 cnec         read(*,*)alphaT
 cnec         read(*,*)alphau
@@ -110,4 +109,6 @@
              alphav=0.1
              alphaP=1.e10
+         print*,'alphaT,alphau,alphav,alphaP'
+     s          ,alphaT,alphau,alphav,alphaP
          if(online.eq.-1) return
          print*,'alpha rappel pour T, u, v, P ',
@@ -142,5 +143,5 @@
 c   Lecture du premier etat des reanalyses.
          call read_reanalyse(1
-     s   ,ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,1)
+     s   ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1)
 
 c-----------------------------------------------------------------------
@@ -171,6 +172,6 @@
            CALL SCOPY( ijp1llm,ucovrea2, 1, ucovrea1 , 1 )
            CALL SCOPY( ijp1llm,tetarea2,1,tetarea1   , 1 )
-           CALL SCOPY( ijp1llm,masserea2,1,masserea1   , 1 )
-           CALL SCOPY( ip1jmp1,psrea2, 1,  psrea1, 1 )
+c          CALL SCOPY( ijp1llm,masserea2,1,masserea1   , 1 )
+c          CALL SCOPY( ip1jmp1,psrea2, 1,  psrea1, 1 )
 
           print*,'LECTURE REANALYSES, pas ',step_rea
@@ -179,5 +180,5 @@
            itau_test=itau
            call read_reanalyse(step_rea
-     s     ,ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,1)
+     s     ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1)
         endif
       else
@@ -205,6 +206,6 @@
             a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
             teta(ij,l)=alphaT*teta(ij,l)+(1-alphaT)*a
-            a=(1.-tau)*masserea1(ij,l)+tau*masserea2(ij,l)
-            masse(ij,l)=alphaP*masse(ij,l)+(1-alphaP)*a
+c           a=(1.-tau)*masserea1(ij,l)+tau*masserea2(ij,l)
+c           masse(ij,l)=alphaP*masse(ij,l)+(1-alphaP)*a
          enddo
          do ij=1,ip1jm
Index: DZ.3.3/trunk/libf/dyn3d/offline.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/offline.F	(revision 203)
+++ 	(revision )
@@ -1,640 +1,0 @@
-      PROGRAM offline
-      USE ioipsl
-      IMPLICIT NONE
-
-c      ......   Version  du 17/04/96    ..........
-c=======================================================================
-c
-c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
-c   -------
-c   Modif special traceur F.Forget 05/94
-c
-c   Modifs M-A Filiberti 07/99
-C
-c   Objet:
-c   ------
-c
-c   GCM LMD nouvelle grille
-c
-c=======================================================================
-
-c   ...      modification de l'integration de q ( 26/04/94 )      ....
-
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comdissnew.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "logic.h"
-#include "temps.h"
-#include "control.h"
-#include "ener.h"
-#include "description.h"
-#include "serre.h"
-#include "tracstoke.h"
-      integer ngridmx,nbsrf
-      PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4)
-c-----------------------------------------------------------------------
-
-      INTEGER*4  iday ! jour julien
-      REAL       time ! Heure de la journee en fraction d'1 jour
-      REAL dtdyn,dtav
-      INTEGER iday_step
-
-      INTEGER         longcles
-      PARAMETER     ( longcles = 20 )
-      REAL  clesphy0( longcles )
-      logical clefinit
-
-
-c   variables dynamiques
-      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
-      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
-      REAL pks(iip1,jjp1)                      ! exner (f pour filtre)
-      REAL phis(iip1,jjp1)                     ! geopotentiel au sol
-      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
-      REAL w(iip1,jjp1,llm)                    ! vitesse verticale
-C
-C   variables physiques
-      real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
-      REAL qfi(ngridmx,llm,nqmx)               ! champs advectes
-      real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
-      real coefkz(ngridmx,llm)
-      real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
-      real pphis (ngridmx)
-      REAL airefi(ngridmx)
-      REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point
-      real yu1(ngridmx), yv1(ngridmx)
-      real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
-
-C
-      character*10 file
-
-c variables dynamiques intermediaire pour le transport
-      REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse
-
-      REAL masse(iip1,jjp1,llm)
-      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
-      REAL teta(iip1,jjp1,llm)
-
-      INTEGER itau,irec,nrec,ldec,itauav
-      integer recmin,recmax,recstkmin,recstkmax
-      character*1 injecte,stoke
-      character*10 nom
-
-      integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
-      logical debut,lectstart,rnpb
-
-      EXTERNAL inidissip,iniconst,inifilr
-      EXTERNAL inigeom
-      EXTERNAL exner_hyb,pression
-      EXTERNAL defrun_new
-
-      REAL time_0 ! temps de debut de simulation (% journee)
-
-      character*2 str2
-      INTEGER iq,j,i,l,nq,mode
-
-      real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm)
-      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm)
-
-      real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat
-      integer i_stat(nqmx),j_stat(nqmx)
-      real intensite(nqmx)
-
-      integer iqmin,iqmax
-      integer npasinject
-      data npasinject/1/
-
-      integer histid, histvid, histaveid
-      real time_step, t_wrt, t_ops
-      integer itime_step
-      character*80 dynhist_file, dynhistave_file
-
-
-      logical avant
-
-      logical debutphy,redecoupage
-      DATA debutphy/.true./
-      
-      real unskap,pksurcp,smass(iip1,jjp1)
-      integer ig0,ii,jj,iii,ij
-
-     
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
-      REAL ps(iip1,jjp1),pkf(ip1jmp1,llm) 
-
-
-c-----------------------------------------------------------------------
-c   Initialisations:
-c   ----------------
-
-      descript = 'Run GCM LMDZ'
-
-c-----------------------------------------------------------------------
-
-
-      rad = 6400000
-      omeg = 7.272205e-05
-      g = 9.8
-      kappa = 0.285716
-      daysec = 86400
-      cpp = 1004.70885
-Cmaf a verifier....
-      preff = 101325.
-      pa= 50000.
-C
-
-      print*,'Nombre effectif de traceurs'
-      read(*,*) nq
-      write(*,*) nq
-      print*,'Splitting du pas de temps dynamique'
-      read(*,*) nsplit_dyn
-      write(*,*) nsplit_dyn
-      print*,'Splitting du pas de temps physique'
-      read(*,*) nsplit_phy
-      write(*,*) nsplit_phy
-      print*,'jour0, par rapport au debut du fichier fluxmass'
-      read(*,*) jour0
-      write(*,*) jour0
-      print*,'redecoupage (T) ou standard (F)'
-      read(*,'(l1)') redecoupage
-      write(*,*) redecoupage
-      print*,'Lecture fichier start (T)?'
-      read(*,'(l1)') lectstart
-      write(*,'(l1)') lectstart
-      print*,'Traceurs 1 et 2 = Rn222 et Pb'
-      read(*,'(l1)') rnpb
-      write(*,'(l1)') rnpb
-      print*,'Avant T ou arriere F?'
-      read(*,'(l1)') avant
-      write(*,'(l1)') avant
-      print*,'Index min et max des traceurs a initaliser'
-      read(*,*) iqmin,iqmax
-      write(*,*) iqmin,iqmax
-      print*,'localisation (lon, lat, intensite)'
-      if(iqmax.le.0) then
-         iqmin=iqmax+1
-      endif
-      do iq=iqmin,iqmax
-         read(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
-         write(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
-      enddo
-      print*,'Injection sur combien de pas de temps'
-      read(*,*) npasinject
-      write(*,*) npasinject
-      print*,'Records min et max pour le stokage et couche '
-     s    ,'verticale (specifique reseaux), def 0 0 1'
-      read(*,*) recstkmin,recstkmax,ldec
-      write(*,*) recstkmin,recstkmax,ldec
-
-c   redefinition du splitting
-      if(mod(nsplit_dyn,nsplit_phy).eq.0) then
-         nsplit=nsplit_phy
-      elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then
-         nsplit=nsplit_dyn
-      else
-         stop'prendre un des deux nsplit doit diviser l autre'
-      endif
-      nsplit_dyn=nsplit_dyn/nsplit
-      nsplit_phy=nsplit_phy/nsplit
-      print*,'nsplit,nsplit_dyn,nsplit_phy'
-     s       ,nsplit,nsplit_dyn,nsplit_phy
-c   Lecture eventuelle 
-      time_0=0.
-
-      if(nq.gt.nqmx) stop'surdimensioner nqmx'
-
-      print*,'av initial0'
-      call initial0(nq*ijp1llm,q)
-
-      print*,'av lectstart'
-      if(lectstart) then
-         CALL dynetat0("start.nc",nq,vcov,ucov,
-     .              teta,q,masse,ps,phis, time_0)
-         print*,'Lecture du start'
-      else
-         day_ini=0
-         anne_ini=0
-      endif
-      day_end=day_ini+nday
-
-
-      print*,'av defrun'
-c     open (99,file='run.def',status='old',form='formatted')
-      clefinit=.true.
-      CALL defrun_new( 99, clefinit, clesphy0 )
-c     close(99)
-
-      print*,'av iniconst'
-c   lecture du jour de demarrage
-c    premiere initialisation, eventuellement bidon
-      CALL iniconst
-      CALL inigeom
-
-C
-c  Lecture des entetes des fichiers de flux
-      if(redecoupage) then
-         call redecoupe(0,masse,pbaru,pbarv,w,teta,phi,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf,
-     s     frac_impa,frac_nucl)
-      else
-         call lectflux(0,masse,pbaru,pbarv,w,teta,phi,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf,
-     s     frac_impa,frac_nucl,phis)
-         call inigeom
-      endif
-
-c   controle sur les records
-      if(avant) then
-         recmin=2+(jour0-1)*iday_step
-         recmax=1+(jour0+nday-1)*iday_step
-      else
-         recmin=(jour0-nday)*iday_step+2
-         recmax=jour0*iday_step+1
-      endif
-      if(recmin.le.1.or.recmax.ge.nrec) then
-         print*,'Probleme de lecture '
-         print*,'recmin,recmax,nrec',recmin,recmax,nrec
-         stop
-      endif
-
-C on recalcule le nombre de pas de temps dans une journee
-c   istdyn est la frequence de stokage en "pas de temps dynamique" dans
-c   le modele on-line. dtvr est le pas de temps du meme modele en s.
-c iday-step est le nombre de lectures par jour (grand pas de temps). 
-      iday_step=daysec/(istdyn*dtvr)
-C
-      CALL iniconst
-c  on calcule le pas de temps 
-      dtvr    = daysec/FLOAT(iday_step)
-      dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst!
-      print*,'Pas de temps physique ',dtphys
-
-
-C
-
-      call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
-
-      print *,'Pas de temps dynamique', dtvr
-      print*,'Pas de temps physique ',dtphys
-
-c-----------------------------------------------------------------------
-C maf est ce utile ?
-
-      call suphec
-
-      CALL inifilr
-      print*,'Pas de temps physique ',dtphys
-      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
-     *                tetagdiv, tetagrot , tetatemp              )
-      print*,'Pas de temps physique ',dtphys
-
-c-----------------------------------------------------------------------
-c   temps de depart et de fin:
-c   --------------------------
-
-      itau = 0
-      iday = dayref
-      time = time_0
-         IF(time.GT.1.) THEN
-          time = time-1.
-          iday = iday+1
-         ENDIF
-      print*,'Pas de temps physique ',dtphys
-
-c=======================================================================
-c   Initialisation des fichiers de sortie
-c=======================================================================
-
-      print*,'Initialisation de wrgras'
-      file='sol'
-      call inigrads(1,iip1
-     s  ,rlonv,180./pi,-190.,190.,jjp1,rlatu,-90.,90.,180./pi
-     s  ,llm,presnivs,1.
-     s  ,dtphys,file,'gcmq2 ')
-      file='atm'
-      print*,'Pas de temps physique ',dtphys
-      call inigrads(2,iip1
-     s  ,rlonv,180./pi,-190.,190.,jjp1,rlatu,-90.,90.,180./pi
-     s  ,llm,presnivs,1.
-     s  ,dtphys,file,'gcmq2 ')
-
-      print*,'dtphys ap inigrads ',dtphys
-
-      print*,'Ecriture du fichier de conditions aux limites'
-      open (98,file='limit',form='unformatted')
-      write(98) float(im),float(jm),float(nq),float(nday)
-      write(98) (rlonv(i)*180./pi,i=1,iip1)
-      write(98) (rlatu(j)*180./pi,j=1,jjp1)
-      write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1)
-      close(98)
-
-      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nq)
-
-c-----------------------------------------------------------------------
-c   Debut de la boucle en temps:
-c   ----------------------------
-
-      do itau=1,nday*iday_step
-
-          print*,'ITAU ITAU ITAU =',itau
-         injecte='.'
-         stoke='.'
-
-c Gestion du temps
-         if (avant) then
-            dtdyn=dtvr/float(nsplit*nsplit_dyn)
-            irec=1+(jour0-1)*iday_step+itau
-         else
-            dtdyn=-dtvr/float(nsplit*nsplit_dyn)
-            irec=1+jour0*iday_step-itau+1
-         endif
-
-         debut=itau.eq.1
-
-c-----------------------------------------------------------------------
-c   Lecture des fichiers fluxmass et physique:
-c   -----------------------------------------------------
-
-         print*,'CCCCC Appel a redecoupe au pas itau=',itau
-      if(redecoupage) then
-         call redecoupe(irec,masse,pbaru,pbarv,w,teta,phi,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf, 
-     s     frac_impa,frac_nucl)
-      else
-         call lectflux(irec,masse,pbaru,pbarv,w,teta,phi,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf, 
-     s     frac_impa,frac_nucl,phis)
-      endif
-
-c    ...  ouverture du fichier de stockage netcdf ...
-C
-      IF (debut) then
-         dynhistave_file = 'dyn_hist_ave.nc'
-         day_ini=0
-         anne_ini=0
-         t_ops = periodav * daysec
-         t_wrt = periodav * daysec
-         dtav=dtvr/float(nsplit)
-         mode=1
-         CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
-     .              t_ops, t_wrt, nq, histaveid)
-         pi=2.*asin(1.)
-         do iq=iqmin,iqmax
-            do l=1,llm
-               do j=1,jjp1
-                  do i=1,iip1
-                     q(i,j,l,iq)=0.
-                  enddo
-               enddo
-            enddo
-            zlon=lon_stat(iq)*pi/180.
-            zlat=lat_stat(iq)*pi/180.
-            CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq))
-         enddo
-      ENDIF
-
-C calcul de la pression de surface.
-      do j=1,jjp1
-         do i=1,iip1
-            smass(i,j)=0.
-         enddo
-      enddo
-      do l=1,llm
-        do j=1,jjp1
-          do i=1,iip1
-             smass(i,j)=smass(i,j)+masse(i,j,l)
-          enddo
-       enddo
-      enddo
-      do j=1,jjp1
-         do i=1,iip1
-            ps (i,j)=smass(i,j)/aire(i,j)*g
-         end do
-      end do
-C
-C calcul de la pression et de pk (fonction d'Exner)
-C
-      CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-c
-C
-      do jj=1,jjp1
-         do ii=1,iip1
-           phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)*
-     s                             (pks(ii,jj)-pk(ii,jj,1))
-         end do
-      end do
-
-c=======================================================================
-c   TERMES SOURCES OU PUITS
-c=======================================================================
-c   terme source reparti sur npasinject pas
-c   On injecte tous les traceurs de la meme facon
-      if (itau.le.npasinject) then
-        injecte='X'
-         do iq=iqmin,iqmax
-            print*,'INJECTION AU POINT ',lon_stat(iq),lat_stat(iq)
-     s      ,'     I,J=',i_stat(iq),j_stat(iq),jjp1-j_stat(iq)+1
-            q(i_stat(iq),j_stat(iq),1,iq)=
-     s      q(i_stat(iq),j_stat(iq),1,iq)
-     s      +intensite(iq)/masse(i_stat(iq),j_stat(iq),1)/npasinject
-         enddo
-         print*,'fin des injections'
-      endif
-c=======================================================================
-c   FIN DE LA PARTIE TERMES SOURCES
-c=======================================================================
-c-----------------------------------------------------------------------
-c   Advection
-c   ----------
-
-c=======================================================================
-         do isplit=1,nsplit
-c=======================================================================
-CC MAF dans le cas du offline on ne fait pas la difference masse et massem
-            do iii=1,nsplit_dyn
-               do iq=1,nq
-                  CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn)
-               enddo
-
-c   mise à jour du champd de masse tenant compte de l'advection sur
-c   le pas de temps considere.
-
-            enddo
-
-            do iq=1,nq
-               write(str2,'(i2.2)') iq
-               if(mod(itau,1).eq.0) 
-     s         call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
-            enddo
-
-
-c-----------------------------------------------------------------------
-c   Physique  (lessivage dans phytrac maintenant) MAF
-c   ---------------------------------------------------------
-
-        if(physic) then
-         if(debut) then
-         latfi(1)=rlatu(1)*180./pi
-         lonfi(1)=0.
-         DO j=2,jjm
-            DO i=1,iim
-               latfi((j-2)*iim+1+i)= rlatu(j)*180./pi
-               lonfi((j-2)*iim+1+i)= rlonv(i)*180./pi
-            ENDDO
-         ENDDO
-         latfi(ngridmx)= rlatu(jjp1)*180./pi
-         lonfi(ngridmx)= 0.
-         endif
-
-c
-C Calcul de paprs et pplay et de temp
-c   -----------------------------------------------------------------
-c     .... paprs  definis aux (llm +1) interfaces des couches  ....
-c     .... pplay  definis aux (  llm )    milieux des couches  ....
-c   -----------------------------------------------------------------
-C
-      DO l = 1, llmp1
-        paprs( 1,l ) = p(1,1,l)
-        ig0 = 2
-          DO j = 2, jjm
-             DO i =1, iim
-              paprs( ig0,l ) = p(i,j,l)
-              ig0 = ig0 +1
-             ENDDO
-          ENDDO
-        paprs( ngridmx,l ) = p(1,jjp1,l)
-      ENDDO
-c
-      unskap   = 1./ kappa
-      DO l=1,llm
-
-         pksurcp     =  pk(1,1,l) / cpp
-         pplay(1,l)  =  preff * pksurcp ** unskap
-         tempfi(1,l)   =  teta(1,1,l) *  pksurcp
-         ig0         = 2
-
-         DO j = 2, jjm
-            DO i = 1, iim
-              pksurcp        = pk(i,j,l) / cpp
-              pplay(ig0,l)   = preff * pksurcp ** unskap
-              tempfi(ig0,l)    = teta(i,j,l)  * pksurcp
-              ig0            = ig0 + 1
-            ENDDO
-         ENDDO
-
-         pksurcp       = pk(1,jjp1,l) / cpp
-         pplay(ig0,l)  = preff * pksurcp ** unskap
-         tempfi (ig0,l)  = teta(1,jjp1,l)  * pksurcp
-      ENDDO
-c On passe le traceur et le geopot de surf  sur la grille physique
-             call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi)
-             call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis)
-
-             do iii=1,nsplit_phy
-C
-      print*,'dtphys avant phytrac ',dtphys
-             call phytrac(rnpb,
-     I                   ecritphy, debutphy,
-     I                   nq,
-     I                   ngridmx,llm,dtphys,
-     I                   tempfi,paprs,pplay,
-     I                   zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
-     I                   coefkz,yu1,yv1,ftsol,pctsrf,latfi,
-     I                   frac_impa, frac_nucl,
-     I                   lonfi,presnivs,airefi,pphis,
-     O                   qfi)
-C
-             debutphy= .FALSE.
-C
-             enddo ! nsplit_phy
-C
-c   On passe le traceur sur la grille dynamique.
-             call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q)
-
-             endif
-
-             itauav=(itau-1)*nsplit+isplit
-             CALL writedynav(histaveid, nq, itauav,vcov ,
-     ,                   ucov,teta,pk,phi,q,masse,ps,phis)
-
-         enddo ! isplit=1,nsplit
-
-         call histsync
-
-c-----------------------------------------------------------------------
-c    Calcul de ucov et vcov pour les sorties
-c   -------------------------------------
-
-         call massbar(masse, massebx, masseby )
-         do l=1,llm
-            do j=1,jjp1
-               do i=1,iip1
-               ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)*cu(i,j)*cu(i,j)
-     s         /istdyn
-               enddo
-            enddo
-            do j=1,jjm
-               do i=1,iip1
-               vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)*cv(i,j)*cv(i,j)
-     s         /istdyn
-               enddo
-            enddo
-         enddo
-
-         iday= dayref+itau/iday_step
-         time=
-     s   float(itau-(iday-dayref)*iday_step)/iday_step+time_0
-         IF(time.GT.1.) THEN
-            time = time-1.
-            iday = iday+1
-         ENDIF
-
-
-
-C
-c=======================================================================
-         if(irec.ge.recstkmin.and.irec.le.recstkmax) then
-            if (mod(itau,4).eq.0) then
-            stoke='X'
-               do iq=1,nq
-                  write(str2,'(i2.2)') iq
-                  nom='q'//str2
-                  call wrgrads(1,1,q(:,:,:,iq),nom,nom)
-                  call wrgrads(2,1,q(:,:,ldec,iq),nom,nom)
-               enddo
-            endif
-         endif
-
-         print*,'Record=',irec,'   stoke=',stoke,
-     s      '      injecte=',injecte
-
-      enddo ! itau
-
-      CALL dynredem1("restart.nc",0.0,
-     .          vcov,ucov,teta,q,nq,masse,ps)
-
-c=======================================================================
-C fermeture des fichiers netcdf
-          CALL histclo
-C
- 300  FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x,
-     s 'c"est a dire dujour',i7,3x,'au jour',i7,/,/)
-
-      end
-C----------------------------------------------------------------------
Index: /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F	(revision 204)
@@ -1,2 +1,5 @@
+c
+c $Header
+c
       PROGRAM offlinenc
       USE ioipsl
@@ -66,5 +69,5 @@
 	integer iiinj
       integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
-      logical debut,lectstart,rnpb
+      logical debut,lectstart,rnpb,lafin
 
       EXTERNAL inidissip,iniconst,inifilr
@@ -190,4 +193,9 @@
      .              teta,q,masse,ps,phis, time_0)
          print*,'Lecture du start'
+c  On zappe le radon et le plomb.
+         q(:,:,:,1)=q(:,:,:,3)
+         q(:,:,:,2)=q(:,:,:,4)
+c Initialisation d'un traceur ` 1 pour tester l'impact du lessivage.
+         q(:,:,:,3)=1.
       else
          day_ini=0
@@ -203,9 +211,9 @@
 C      print*,'av iniconst'
 c   lecture du jour de demarrage
-c    premiere initialisation, eventuellement bidon
-      CALL iniconst
+c    premiere cnitialisation, eventuellement bidon
+      cALL iniconst
       CALL inigeom
 C
-      print*,'ENTREE DANS redecoupenc ou lectfluxnc0
+      print*,'ENTREE DANS redecoupenc ou lectfluxnc
      s        pour irec=0'
 
@@ -219,5 +227,5 @@
      s     frac_impa,frac_nucl,phis)
       else
-         call lectfluxnc0(0,masse,pbaru,pbarv,w,teta,phi,
+         call lectfluxnc(0,masse,pbaru,pbarv,w,teta,phi,
      s     nrec,avant,airefi,pphis,
      s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
@@ -373,5 +381,5 @@
      s     frac_impa,frac_nucl,phis)
       else
-         call lectfluxnc0(irec,masse,pbaru,pbarv,w,teta,phi,
+         call lectfluxnc(irec,masse,pbaru,pbarv,w,teta,phi,
      s     nrec,avant,airefi,pphis,
      s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
@@ -380,4 +388,9 @@
       endif
 
+
+      print*,'TESTPHYS: ON PREND LA PUIS ',1./float(nsplit*nsplit_phy)
+     s   ,'  DES FRAC A IT=',itau
+      frac_impa(:,:)=frac_impa(:,:)**(1./float(nsplit*nsplit_phy))
+      frac_nucl(:,:)=frac_nucl(:,:)**(1./float(nsplit*nsplit_phy))
 c    ...  ouverture du fichier de stockage netcdf ...
 C
@@ -391,6 +404,6 @@
          mode=1
 
-         CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
-     .              t_ops, t_wrt, nq,mode, histaveid)
+c        CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
+c    .              t_ops, t_wrt, nq,mode, histaveid)
 
          pi=2.*asin(1.)
@@ -832,7 +845,11 @@
              do iii=1,nsplit_phy
 C
+      lafin=.false. ! en attendant mieux.
       print*,'dtphys avant phytrac ',dtphys
-             call phytrac(rnpb,
-     I                   ecritphy, debutphy,
+      print*,'TESTPHYS: APPEL A PHYTRAC IT=',itau
+
+             call phytrac(rnpb,ecritphy,
+c    I                   ecritphy, debutphy,
+     I                   debutphy, lafin,
      I                   nq,
      I                   ngridmx,llm,dtphys,
@@ -855,6 +872,6 @@
              itauav=(itau-1)*nsplit+isplit
 c             itauav=itau*nsplit+isplit
-             CALL writedynav(histaveid, nq,mode, itauav,vcov ,
-     ,                   ucov,teta,pk,phi,q,masse,ps,phis)
+c            CALL writedynav(histaveid, nq,mode, itauav,vcov ,
+c    ,                   ucov,teta,pk,phi,q,masse,ps,phis)
 c            qmoy(:,:,:)=qmoy(:,:,:)+q(:,:,1,:)
          do iq=1,nq
Index: /LMDZ.3.3/trunk/libf/dyn3d/read_reanalyse.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/read_reanalyse.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/dyn3d/read_reanalyse.F	(revision 204)
@@ -1,3 +1,5 @@
+c
 c $Header
+c
       subroutine read_reanalyse(timestep,u,v,t,masse,ps,mode)
 
@@ -14,5 +16,5 @@
       integer nlevnc
 cModef 11-2-99      parameter (nlevnc=15)
-      parameter (nlevnc=15)
+      parameter (nlevnc=21)
       integer timestep,mode,l
 
@@ -58,7 +60,7 @@
             varidt=NCVID(ncidt,'AIR',rcode)
             print*,'ncidt,varidt',ncidt,varidt
-            ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
-            varidps=NCVID(ncidps,'SP',rcode)
-            print*,'ncidps,varidps',ncidps,varidps
+c           ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
+c           varidps=NCVID(ncidps,'SP',rcode)
+c           print*,'ncidps,varidps',ncidps,varidps
       endif
 
@@ -106,7 +108,7 @@
       count(4)=0
 c  ps
-      status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
+c     status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
 c     call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
-      call correctbid(iim,jjp1,psnc)
+c     call correctbid(iim,jjp1,psnc)
 
 c   Transformations
@@ -187,10 +189,10 @@
 c -----------------------------------------------------------------
 
-      do j=1,jjp1
-         do i=1,iim
-            ps(i,j)=psnc(i,jjp1+1-j)
-         enddo
-         ps(iip1,j)=ps(1,j)
-      enddo
+c     do j=1,jjp1
+c        do i=1,iim
+c           ps(i,j)=psnc(i,jjp1+1-j)
+c        enddo
+c        ps(iip1,j)=ps(1,j)
+c     enddo
 
       CALL pression( ip1jmp1, ap, bp, ps, p )
@@ -200,5 +202,5 @@
       unskap=1./kappa
       prefkap =  preff  ** kappa
-      PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
+c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
       DO l = 1, llm
        DO j=1,jjp1
Index: DZ.3.3/trunk/libf/dyn3d/redecoupe.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/redecoupe.F	(revision 203)
+++ 	(revision )
@@ -1,637 +1,0 @@
-      SUBROUTINE redecoupe(irec,massemn,pbarun,pbarvn,wn,tetan,phin,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf,
-     s     frac_impa,frac_nucl,phisn)
-
-      IMPLICIT NONE
-
-#include "dimensions.h"
-#include "paramet.h"
-
-#include "comvert.h"
-#include "comconst.h"
-#include "comgeom2.h"
-
-#include "tracstoke.h"
-
-      integer irec,nrec,i,j
-
-      integer ig,l
-      integer imo,jmo,imn,jmn,ii,jj,ig
-      parameter (imn=iim,jmn=jjm,imo=imn/2,jmo=(jmn+1)/2)
-      integer ngrido,ngridn
-      parameter(ngrido=(jmo-1)*imo+2,ngridn=(jmn-1)*imn+2)
-
-
-      INTEGER nbsrf
-      PARAMETER (nbsrf=4) ! nombre de sous-fractions pour une maille
-
-      real zmfd(ngridn,llm),zde_d(ngridn,llm),zen_d(ngridn,llm)
-      real zmfu(ngridn,llm),zde_u(ngridn,llm),zen_u(ngridn,llm)
-      real coefkz(ngridn,llm)
-      real frac_impa(ngridn,llm),frac_nucl(ngridn,llm)
-      real yu1(ngridn), yv1(ngridn)
-      real ftsol(ngridn,nbsrf),pctsrf(ngridn,nbsrf)
-      integer imfu,imfd,ien_u,ide_u,
-     s      ien_d,ide_d,
-     s      icoefkz,izu1,izv1,
-     s      itsol,ipsf,
-     s      ilei, ilec
-      parameter(imfu=1,imfd=llm+1,ien_u=2*llm+1,ide_u=3*llm+1,
-     s      ien_d=4*llm+1,ide_d=5*llm+1,
-     s      icoefkz=6*llm+1,
-     s      ilei=7*llm+1,ilec=8*llm+1,
-     s      izu1=9*llm+1,izv1=9*llm+2,
-     s      itsol=9*llm+3,ipsf=9*llm+3+nbsrf)
-      logical avant
-
-
-      real massefi(ngridn,llm)
-
-      real massemn(imn+1,jmn+1,llm),tetan(imn+1,jmn+1,llm)
-      real pbarun(imn+1,jmn+1,llm),pbarvn(imn+1,jmn,llm)
-      real wn(imn+1,jmn+1,llm),phin(imn+1,jmn+1,llm)
-      real phisn(imn+1,jmn+1)
-
-      real massemo(imo+1,jmo+1,llm),tetao(imo+1,jmo+1,llm)
-      real pbaruo(imo+1,jmo+1,llm),pbarvo(imo+1,jmo,llm)
-      real wo(imo+1,jmo+1,llm),phio(imo+1,jmo+1,llm)
-      real phiso(imo+1,jmo+1)
-
-      real pbarvst(imo+1,jmo+1,llm)
-
-      real airefi(ngridn)
-
-      real xlecn(ngridn,9*llm+2+2*nbsrf),tmpn(imn+1,jmn+1)
-      real xleco(ngrido,9*llm+2+2*nbsrf),tmpo(imo+1,jmo+1)
-
-      real zcontrole(ngridn),zmass,tmpdyn(imn+1,jmn+1),zflux
-
-      real ziadvtrac,zrec,ziadvtrac2,zrec2
-      real zim,zjm,zlm,zklon,zklev
-
-      real zpi
-c  longitudes et latitudes lues
-      real rlonul(1:imo+1),rlatvl(1:jmo)
-      real rlonvl(1:imo+1),rlatul(1:jmo+1)
-c  longitudes et latitudes anciennes
-      real rlonuo(0:imo+1),rlatvo(0:jmo+1)
-c  longitudes et latitudes nouvelles
-      real rlonun(0:imn+1),rlatvn(0:jmn+1)
-      real aireo(imo+1,jmo+1)
-
-      integer ndecx(imo+1),ndecy(jmo+1)
-      real alphax(imn+1),alphay(jmn+1)
-      real alphaxo(imo+1)
-      real alpha(imn+1,jmn+1)
-
-      real aa,uu(0:imo+1),vv(imo+1,0:jmo+1)
-
-
-      integer iest(imo+1),iouest(imo+1)
-      integer jsud(jmo+1),jnord(jmo+1)
-
-      integer in,io,jn,jo,l,iin,jjn
-      integer i,j
-      real dlatm,dlatp,dlonm,dlonp
-
-
-      zpi=2.*asin(1.)
-
-
-c==================================================================
-c   Si le numero du record est 0 alors: INITIALISATION
-c==================================================================
-c
-      print*,'ENTREE DANS LECTFLUX'
-        print*,'IREC=',IREC
-      if(irec.eq.0) then
-
-        print*,'IREC==',0
-
-C test         call inigeom
-c==================================================================
-c   Definition des surdecoupages dans les deux directions
-c==================================================================
-
-      ndecx(1)=1
-      do io=2,imo
-         ndecx(io)=2
-      enddo
-      ndecx(imo+1)=1
-
-      ndecy(1)=1
-      do jo=2,jmo
-         ndecy(jo)=2
-      enddo
-      ndecy(jmo+1)=1
-
-      ii=0
-      do io=1,imo+1
-         ii=ii+ndecx(io)
-      enddo
-      if(ii.ne.iim) then
-         print*,'ii=',ii,'   iim=',iim
-         stop
-      endif
-
-      jj=0
-      do jo=1,jmo+1
-         jj=jj+ndecy(jo)
-      enddo
-      if(jj.ne.jjp1) then
-         print*,'jj=',jj,'   jjm=',jjm
-         stop
-      endif
-
-c==================================================================
-c   Calcul des jsud,... correspondant aux intersections des
-c   grilles.
-c==================================================================
-
-      iest(1)=0
-      do io=2,imo+1
-         iest(io)=iest(io-1)+ndecx(io-1)
-         iouest(io-1)=iest(io)
-      enddo
-      iouest(imo+1)=iest(imo+1)+ndecx(imo+1)
-
-      jnord(1)=0
-      do jo=2,jmo+1
-         jnord(jo)=jnord(jo-1)+ndecy(jo-1)
-         jsud(jo-1)=jnord(jo)
-      enddo
-      jsud(jmo+1)=jnord(jmo+1)+ndecy(jmo+1)
-
-c==================================================================
-c   ouverture des fichiers, lecture des entetes
-c==================================================================
-
-c   Fichier fluxmass
-#ifdef CRAY
-         CALL ASSIGN("assign -N ieee -F null f:fluxmass")
-#endif
-      open(47,file='fluxmass',form='unformatted',
-     s     access='direct'
-     s     ,recl=4*6*(imo+1)*(jmo+1)*llm)
-      read(47,rec=1) zrec,dtvr,ziadvtrac,zim,zjm,zlm
-     s   ,rlonul,rlonvl,rlatul,rlatvl,aireo
-     s    ,phiso
-
-      print*,'zrec,dtvr,ziadvtrac,zim,zjm,zlm'
-      print*,zrec,dtvr,ziadvtrac,zim,zjm,zlm
-      if((imo-nint(zim))*(jmo-nint(zjm)).ne.0) then
-        print*,'Modifier les dimensions dans redecoupe '
-        print*,'Mettre imo=',zim,'   jmo=',zjm
-        stop
-      endif
-
-
-c  Fichier physique
-c  Fichier lessivage (supprime les donnees utiles sont dans "physique")
-#ifdef CRAY
-         CALL ASSIGN("assign -N ieee -F null f:physique")
-#endif
-      open(49,file='physique',form='unformatted',
-     s     access='direct'
-     s     ,recl=4*ngrido*(9*llm+2+2*nbsrf))
-      read(49,rec=1) zrec2,ziadvtrac2,zklon,zklev
-      print*,'Entete du fichier physique'
-      print*,zrec2,ziadvtrac2,zklon,zklev
-
-      nrec=zrec
-      istdyn=ziadvtrac
-      istphy=ziadvtrac2
-
-c==================================================================
-c   Definition des anciennes latitudes et longitudes
-c   (qui pourraient etre relues plus tard)
-c==================================================================
-
-      rlonuo(0)=-zpi
-      do io=1,imo
-c        rlonuo(io)=2.*zpi/FLOAT(imo)*(io+0.5-0.5*FLOAT(imo)-1.)
-c        print*,'LON ',io,rlonuo(io),rlonul(io)
-         rlonuo(io)=rlonul(io)
-      enddo
-      rlonuo(imo+1)=zpi
-
-      rlatvo(0)=zpi/2.
-      do jo=1,jmo
-c        rlatvo(jo)=zpi/FLOAT(jmo)*(0.5*FLOAT(jmo)+1.-jo-0.5)
-c        print*,'LAT ',jo,rlatvo(jo),rlatvl(jo)
-         rlatvo(jo)=rlatvl(jo)
-      enddo
-      rlatvo(jmo+1)=-zpi/2.
-
-c     do jo=1,jmo+1
-c        do io=1,imo+1
-c           aireo(io,jo)=rad*rad
-c    s         *(rlonuo(io)-rlonuo(io-1))
-c    s         *(sin(rlatvo(jo-1))-sin(rlatvo(jo)))
-c           aireo(io,jo)=airel(io,jo)
-c        enddo
-c        aireo(1,jo)=aireo(1,jo)+aireo(imo+1,jo)
-c        aireo(imo+1,jo)=aireo(1,jo)
-c     enddo
-
-      do io=2,imo
-         alphaxo(io)=1.
-      enddo
-      alphaxo(1)=(rlonuo(1)-rlonuo(0))
-     s        /(rlonuo(1)-rlonuo(0)+rlonuo(imo+1)-rlonuo(imo))
-      alphaxo(imo+1)=1.-alphaxo(1)
-
-c==================================================================
-c    Definition des nouvelles latitudes et longitudes
-c==================================================================
-
-      rlonun(0)=-zpi
-      do io=1,imo+1
-         do iin=1,iouest(io)-iest(io)
-            in=iin+iest(io)
-            rlonun(in)=
-     s      rlonuo(io-1)+iin*(rlonuo(io)-rlonuo(io-1))
-     s      /ndecx(io)
-            alphax(in)=alphaxo(io)/ndecx(io)
-            print787,io,rlonuo(io-1)*180./zpi,in
-     s      ,iest(io),iouest(io),rlonun(in)*180./zpi,alphax(in)
-         enddo
-      enddo
-
-      rlatvn(0)=0.5*zpi
-      do jo=1,jmo+1
-         print*,'jo=',jo
-         do jjn=1,jsud(jo)-jnord(jo)
-            jn=jnord(jo)+jjn
-            rlatvn(jn)=rlatvo(jo-1)+jjn*(rlatvo(jo)-rlatvo(jo-1))
-     s      /ndecy(jo)
-            alphay(jn)=(sin(rlatvn(jn-1))-sin(rlatvn(jn)))
-     s                /(sin(rlatvo(jo-1))-sin(rlatvo(jo)))
-            print787,jo,rlatvo(jo-1)*180./zpi,jn
-     s      ,jnord(jo),jsud(jo),rlatvn(jn)*180./zpi,alphay(jn)
-         enddo
-      enddo
-
-787   format(i5,f10.2,3(i5),2(f10.2))
-      do in=1,imn
-         rlonu(in)=rlonun(in)
-         rlonv(in)=0.5*(rlonun(in)+rlonun(in-1))
-      enddo
-      rlonv(imn+1)=rlonv(1)+2.*zpi
-      rlonu(imn+1)=rlonu(1)+2.*zpi
-
-      do jn=1,jmn
-         rlatv(jn)=rlatvn(jn)
-      enddo
-      do jn=1,jmn+1
-         rlatu(jn)=0.5*(rlatvn(jn-1)+rlatvn(jn))
-      enddo
-
-      do jn=1,jmn+1
-         do in=1,imn
-            alpha(in,jn)=alphax(in)*alphay(jn)
-         enddo
-         alpha(imn+1,jn)=0.
-      enddo
-
-c     call dump2d(iip1,jjp1,alpha,'ALPHA   ')
-
-c      .  on a :  cu(i,j) = rad * COS(y) * dx/dX         .
-c      .          cv( j ) = rad          * dy/dY         .
-c   A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont
-c   affectees 4 aires entourant P , calculees respectivement aux points
-c            ( i + 1/4, j - 1/4 )    :    aireij1 (i,j)
-c            ( i + 1/4, j + 1/4 )    :    aireij2 (i,j)
-c            ( i - 1/4, j + 1/4 )    :    aireij3 (i,j)
-c            ( i - 1/4, j - 1/4 )    :    aireij4 (i,j)
-c
-c                             . V
-c
-c                 aireij4 .        . aireij1
-c
-c                   U .       . P      . U
-c
-c                 aireij3 .        . aireij2
-c
-c                             . V
-
-
-      do j=1,jjp1
-         do i=1,iim
-            dlonp=rlonun(i)-rlonv(i)
-            dlonm=rlonv(i)-rlonun(i-1)
-            dlatp=sin(rlatvn(j-1))-sin(rlatu(j))
-            dlatm=sin(rlatu(j))-sin(rlatvn(j))
-            aireij1 ( i,j ) = rad*rad*dlatp*dlonp
-            aireij2 ( i,j ) = rad*rad*dlatm*dlonp
-            aireij3 ( i,j ) = rad*rad*dlatm*dlonm
-            aireij4 ( i,j ) = rad*rad*dlatp*dlonm
-      aire    ( i,j )  = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) +
-     *                          aireij4(i,j)
-      alpha1  ( i,j )  = aireij1(i,j) / aire(i,j)
-      alpha2  ( i,j )  = aireij2(i,j) / aire(i,j)
-      alpha3  ( i,j )  = aireij3(i,j) / aire(i,j)
-      alpha4  ( i,j )  = aireij4(i,j) / aire(i,j)
-      alpha1p2( i,j )  = alpha1 (i,j) + alpha2 (i,j)
-      alpha1p4( i,j )  = alpha1 (i,j) + alpha4 (i,j)
-      alpha2p3( i,j )  = alpha2 (i,j) + alpha3 (i,j)
-      alpha3p4( i,j )  = alpha3 (i,j) + alpha4 (i,j)
-           enddo
-           aireij1(iip1,j)=aireij1(1,j)
-           aireij2(iip1,j)=aireij2(1,j)
-           aireij3(iip1,j)=aireij3(1,j)
-           aireij4(iip1,j)=aireij4(1,j)
-           aire(iip1,j)=aire(1,j)
-           alpha1(iip1,j)=alpha1(1,j)
-           alpha2(iip1,j)=alpha2(1,j)
-           alpha3(iip1,j)=alpha3(1,j)
-           alpha4(iip1,j)=alpha4(1,j)
-           alpha1p2(iip1,j)=alpha1p2(1,j)
-           alpha1p4(iip1,j)=alpha1p4(1,j)
-           alpha2p3(iip1,j)=alpha2p3(1,j)
-           alpha3p4(iip1,j)=alpha3p4(1,j)
-       enddo
-c     call dump2d(iip1,jjp1,aire,'AIRE   ')
-
-c     do jn=1,jjp1
-c        do in=1,iim
-c           aire(in,jn)=rad*rad*(sin(rlatvn(jn-1))-sin(rlatvn(jn)))
-c    s      *(rlonun(in)-rlonun(in-1))
-c           unsaire(in,jn)=1./aire(in,jn)
-c        enddo
-c        aire(iip1,jn)=aire(1,jn)
-c        unsaire(iip1,jn)=unsaire(1,jn)
-c     enddo
-c     call dump2d(iip1,jjp1,aire,'AIRE2   ')
-      DO 42 j = 1,jjp1
-      DO 41 i = 1,iim
-      unsaire(i,j) = 1./ aire(i,j)
-      aireu  (i,j) = aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) +
-     *                           aireij3(i+1,j)
-  41  CONTINUE
-      aireu  (iip1,j) = aireu  (1,j)
-      unsaire(iip1,j) = unsaire(1,j)
-  42  CONTINUE
-      DO 48 j = 1,jjm
-        DO i=1,iim
-         airev(i,j)     = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) +
-     *                           aireij4(i,j+1)
-        ENDDO
-       airev   (iip1,j) = airev(1,j)
-  48  CONTINUE
-      apoln=0.
-      apols=0.
-      do i=1,iim
-         apoln=apoln+aire(i,1)
-         apols=apols+aire(i,jjp1)
-      enddo
-
-
-
-      do jn=1,jjp1
-         do in=1,iim
-            cu(in,jn)=rad*cos(rlatu(jn))*(rlonv(in+1)-rlonv(in))
-         enddo
-         cu(iip1,jn)=cu(1,jn)
-      enddo
-      do jn=1,jjm
-         do in=1,iim+1
-            cv(in,jn)=rad*(rlatu(jn)-rlatu(jn+1))
-         enddo
-      enddo
-
-
-c==================================================================
-c   Fin des initialisations
-      else ! irec=0
-c==================================================================
-
-
-c-----------------------------------------------------------------------
-c   Lecture des fichiers fluxmass et  physique:
-c   -----------------------------------------------------
-
-c==================================================================
-c  Variables dynamiques
-c==================================================================
-
-         read(47,rec=irec) massemo,pbaruo,pbarvst,wo,tetao,phio
-        do l=1,llm
-           do j=1,jmo
-              do i=1,imo+1
-                 pbarvo(i,j,l)=pbarvst(i,j,l)
-              enddo
-           enddo
-        enddo
-
-         do l=1,llm
-            do jo=1,jmo+1
-               do io=1,imo+1
-                  do jn=jnord(jo)+1,jsud(jo)
-                     do in=iest(io)+1,iouest(io)
-                        wn(in,jn,l)=alpha(in,jn)*wo(io,jo,l)
-                        massemn(in,jn,l)=alpha(in,jn)
-     s                     *massemo(io,jo,l)
-                        tetan(in,jn,l)=tetao(io,jo,l)
-                        phin(in,jn,l)=phio(io,jo,l)
-                     enddo
-                  enddo
-               enddo
-            enddo
-            do jn=1,jmn+1
-               wn(imn+1,jn,l)=wn(1,jn,l)
-               massemn(imn+1,jn,l)=massemn(1,jn,l)
-               tetan(imn+1,jn,l)=tetan(1,jn,l)
-               phin(imn+1,jn,l)=phin(1,jn,l)
-            enddo
-         enddo
-
-         do l=1,llm
-            do jo=1,jmo+1
-               uu(imo+1)=0.5*(pbaruo(imo,jo,l)+pbaruo(imo+1,jo,l))
-               uu(0)=uu(imo+1)
-               do io=1,imo
-                  uu(io)=pbaruo(io,jo,l)
-               enddo
-               do io=1,imo+1
-                  do jn=jnord(jo)+1,jsud(jo)
-                     aa=0.
-                     do in=iest(io)+1,iouest(io)
-                        aa=aa+alphax(in)
-                        pbarun(in,jn,l)=alphay(jn)*
-     s                    (uu(io-1)+aa*(uu(io)-uu(io-1)))
-                     enddo
-                  enddo
-               enddo
-            enddo
-            do jn=1,jmn+1
-               pbarun(imn+1,jn,l)=pbarun(1,jn,l)
-            enddo
-         enddo
-
-         do l=1,llm
-            do jo=1,jmo
-               do io=1,imo+1
-                  vv(io,jo)=pbarvo(io,jo,l)
-               enddo
-            enddo
-            do io=1,imo+1
-               vv(io,0)=0.
-               vv(io,jmo+1)=0.
-            enddo
-            do jo=1,jmo+1
-               do io=1,imo+1
-                  aa=0.
-c                 do jn=jnord(jo)+1,max(jsud(jo),jmo)
-                  do jn=jnord(jo)+1,min(jsud(jo),jmn)
-                     aa=aa+alphay(jn)
-                     do in=iest(io)+1,iouest(io)
-                        pbarvn(in,jn,l)=alphax(in)*
-     s                  (vv(io,jo-1)+aa*(vv(io,jo)-vv(io,jo-1)))
-                     enddo
-                  enddo
-               enddo
-            enddo
-            do jn=1,jmn
-               pbarvn(iip1,jn,l)=pbarvn(1,jn,l)
-            enddo
-         enddo
-
-c==================================================================
-c  Variables physiques
-c==================================================================
-         read(49,rec=irec) ((xleco(ig,l),ig=1,ngrido),
-     s                                    l=1,9*llm+2+2*nbsrf)
-
-c==================================================================
-c  Passage  a la nouvelle grille
-c==================================================================
-         do l=1,9*llm+2+2*nbsrf
-c   passage aa la grille dynamique ancienne
-            do io=1,imo+1
-               tmpo(io,1)=xleco(1,l)
-               tmpo(io,jmo+1)=xleco(ngrido,l)
-            enddo
-            do jo=2,jmo
-               do io=1,imo
-                  tmpo(io,jo)=xleco((jo-2)*imo+io+1,l)
-               enddo
-               tmpo(imo+1,jo)=tmpo(1,jo)
-            enddo
-c   passage a la grillle dynamique nouvelle
-            do jo=1,jmo+1
-               do io=1,imo+1
-                  do jn=jnord(jo)+1,jsud(jo)
-                     do in=iest(io)+1,iouest(io)
-                        tmpn(in,jn)=tmpo(io,jo)
-                     enddo
-                  enddo
-               enddo
-            enddo
-c   passage a la grille physique nouvelle
-            xlecn(1,l)=tmpn(1,1)
-            xlecn(ngridn,l)=tmpn(1,jmn+1)
-            do jn=2,jmn
-               do in=1,imn
-                  xlecn((jn-2)*imn+in+1,l)=tmpn(in,jn)
-               enddo
-            enddo
-         enddo
-
-c==================================================================
-
-       do l=1,llm
-          do ig=1,ngridn
-             coefkz(ig,l)=xlecn(ig,icoefkz+l-1)
-             frac_impa(ig,l)=xlecn(ig,ilei+l-1)
-             frac_nucl(ig,l)=xlecn(ig,ilec+l-1)
-          enddo
-       enddo
-       do l=1,nbsrf
-          do ig=1,ngridn
-             ftsol(ig,l)=xlecn(ig,itsol+l-1)
-             pctsrf(ig,l)=xlecn(ig,ipsf+l-1)
-          enddo
-       enddo
-       do ig=1,ngridn
-          yv1(ig)=xlecn(ig,izv1)
-          yu1(ig)=xlecn(ig,izu1)
-       enddo
-C
-      if(avant) then
-c   Simu directe
-       do l=1,llm
-          do ig=1,ngridn
-             zmfu(ig,l)=xlecn(ig,imfu+l-1)
-             zmfd(ig,l)=xlecn(ig,imfd+l-1)
-             zde_u(ig,l)=xlecn(ig,ide_u+l-1)
-             zen_u(ig,l)=xlecn(ig,ien_u+l-1)
-             zde_d(ig,l)=xlecn(ig,ide_d+l-1)
-             zen_d(ig,l)=xlecn(ig,ien_d+l-1)
-          enddo
-       enddo
-      else
-c   Simu retro
-       do l=1,llm
-          do ig=1,ngridn
-             zmfd(ig,l)=-xlecn(ig,imfu+l-1)
-             zmfu(ig,l)=-xlecn(ig,imfd+l-1)
-             zen_d(ig,l)=xlecn(ig,ide_u+l-1)
-             zde_d(ig,l)=xlecn(ig,ien_u+l-1)
-             zen_u(ig,l)=xlecn(ig,ide_d+l-1)
-             zde_u(ig,l)=xlecn(ig,ien_d+l-1)
-          enddo
-       enddo
-      endif
-
-c-----------------------------------------------------------------------
-c   PETIT CONTROLE SUR LES FLUX CONVECTIFS...
-c-----------------------------------------------------------------------
-
-         call gr_dyn_fi(llm,iip1,jjp1,ngridn,massemn,massefi)
-
-      print*,'Ap redec irec'
-         do ig=1,ngridn
-            zcontrole(ig)=1.
-         enddo
-c   zmass=(max(massemn(ig,l),massemn(ig,l-1))/airefi(ig)
-         do l=2,llm
-            do ig=1,ngridn
-               zmass=max(massefi(ig,l),massefi(ig,l-1))/airefi(ig)
-               zflux=max(abs(zmfu(ig,l)),abs(zmfd(ig,l)))*dtphys
-               if(zflux.gt.0.9*zmass) then
-                 zcontrole(ig)=min(zcontrole(ig),0.9*zmass/zflux)
-               endif
-            enddo
-         enddo
-
-         do ig=1,ngridn
-            if(zcontrole(ig).lt.0.99999) then
-               print*,'ATTENTION !!! on reduit les flux de masse '
-               print*,'convectifs au point ig=',ig
-            endif
-         enddo
-
-         call gr_fi_dyn(1,ngridn,iip1,jjp1,zcontrole,tmpdyn)
-
-         do l=1,llm
-            do ig=1,ngridn
-               zmfu(ig,l)=zmfu(ig,l)*zcontrole(ig)
-               zmfd(ig,l)=zmfd(ig,l)*zcontrole(ig)
-               zen_u(ig,l)=zen_u(ig,l)*zcontrole(ig)
-               zde_u(ig,l)=zde_u(ig,l)*zcontrole(ig)
-               zen_d(ig,l)=zen_d(ig,l)*zcontrole(ig)
-               zde_d(ig,l)=zde_d(ig,l)*zcontrole(ig)
-            enddo
-         enddo
-
-
-      endif ! irec=0
-
-
-      RETURN
-      END
-
-
Index: /LMDZ.3.3/trunk/libf/dyn3d/redecoupenc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/redecoupenc.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/dyn3d/redecoupenc.F	(revision 204)
@@ -1,9 +1,11 @@
+c
 c $Header
+c
       SUBROUTINE redecoupenc
      s     (irec,massemn,pbarun,pbarvn,wn,tetan,phin,
-     s     nrec,avant,airefi,
-     s     zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
-     s     yu1,yv1,ftsol,pctsrf,
-     s     frac_impa,frac_nucl,phisn)
+     s     nrec,avant,airefin,phisfin,
+     s     tn,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkzn,
+     s     yu1n,yv1n,ftsoln,pctsrfn,
+     s     frac_impan,frac_nucln,phisn)
 
       IMPLICIT NONE
@@ -11,13 +13,11 @@
 #include "dimensions.h"
 #include "paramet.h"
-
 #include "comvert.h"
 #include "comconst.h"
 #include "comgeom2.h"
-
 #include "tracstoke.h"
+#include "logic.h"
 
       integer irec,nrec,i,j
-
       integer ig,l
       integer imo,jmo,imn,jmn,ii,jj,ig
@@ -32,28 +32,6 @@
       real zmfd(ngridn,llm),zde_d(ngridn,llm),zen_d(ngridn,llm)
       real zmfu(ngridn,llm),zde_u(ngridn,llm),zen_u(ngridn,llm)
-      real mfd(ngridn,llm),de_d(ngridn,llm),en_d(ngridn,llm)
-      real mfu(ngridn,llm),de_u(ngridn,llm),en_u(ngridn,llm)
 	
-	real*4 airedy(iip1,jjp1)
-      real*4 rlonu_dy(iip1,jjp1),rlonv_dy(iip1,jjm),
-     . rlatu_dy(iip1,jjp1),rlatv_dy(iip1,jjm)
-
-      real coefkz(ngridn,llm)
-      real frac_impa(ngridn,llm),frac_nucl(ngridn,llm)
-      real yu1(ngridn), yv1(ngridn)
-      real ftsol(ngridn,nbsrf),pctsrf(ngridn,nbsrf)
-      integer imfu,imfd,ien_u,ide_u,
-     s      ien_d,ide_d,
-     s      icoefkz,izu1,izv1,
-     s      itsol,ipsf,
-     s      ilei, ilec
-      parameter(imfu=1,imfd=llm+1,ien_u=2*llm+1,ide_u=3*llm+1,
-     s      ien_d=4*llm+1,ide_d=5*llm+1,
-     s      icoefkz=6*llm+1,
-     s      ilei=7*llm+1,ilec=8*llm+1,
-     s      izu1=9*llm+1,izv1=9*llm+2,
-     s      itsol=9*llm+3,ipsf=9*llm+3+nbsrf)
-      logical avant
-
+       logical avant
 
       real massefi(ngridn,llm)
@@ -63,5 +41,5 @@
       real wn(imn+1,jmn+1,llm),phin(imn+1,jmn+1,llm)
       real phisn(imn+1,jmn+1)
-      real phisfi(imn,jmn+1)
+
       real massemo(imo+1,jmo+1,llm),tetao(imo+1,jmo+1,llm)
       real pbaruo(imo+1,jmo+1,llm),pbarvo(imo+1,jmo,llm)
@@ -71,22 +49,43 @@
       real pbarvst(imo+1,jmo+1,llm)
 
-      real airefi(ngridn)
-
-      real xlecn(ngridn,9*llm+2+2*nbsrf),tmpn(imn+1,jmn+1)
-      real xleco(ngrido,9*llm+2+2*nbsrf),tmpo(imo+1,jmo+1)
+      real tmpo2(imo+1,jmo+1,llm),tmpo1(imo,jmo+1,llm)
+      real tmpo4(imo+1,jmo+1,nbsrf),tmpo3(imo,jmo+1,nbsrf)
+      real tmpo6(imo+1,jmo+1),tmpo5(imo,jmo+1)
+      real tmpn6(imn+1,jmn+1),tmpn5(imn,jmn+1)
+      real tmpn2(imn+1,jmn+1,llm),tmpn1(imn,jmn+1,llm)
+      real tmpn4(imn+1,jmn+1,nbsrf),tmpn3(imn,jmn+1,nbsrf)
+
+	real airefio(ngrido),phisfio(ngrido),
+     .    mfuo(ngrido,llm),mfdo(ngrido,llm),en_uo(ngrido,llm),
+     .    de_uo(ngrido,llm),en_do(ngrido,llm),
+     .    de_do(ngrido,llm),coefkzo(ngrido,llm),
+     .    frac_impao(ngrido,llm),frac_nuclo(ngrido,llm),
+     .    yu1o(ngrido),yv1o(ngrido),ftsolo(ngrido,nbsrf),
+     .    pctsrfo(ngrido,nbsrf),to(ngrido,llm) 
+
+      real airefin(ngridn),phisfin(ngridn),
+     .    mfun(ngridn,llm),en_un(ngridn,llm),mfdn(ngridn,llm),
+     .    de_un(ngridn,llm),en_dn(ngridn,llm),
+     .    de_dn(ngridn,llm),coefkzn(ngridn,llm),
+     .    frac_impan(ngridn,llm),frac_nucln(ngridn,llm),
+     .    ftsoln(ngridn,nbsrf),yu1n(ngridn),yv1n(ngridn),
+     .    pctsrfn(ngridn,nbsrf),tn(ngridn,llm)
 
       real zcontrole(ngridn),zmass,tmpdyn(imn+1,jmn+1),zflux
 
-      real ziadvtrac,zrec,ziadvtrac2,zrec2
-      real zim,zjm,zlm,zklon,zklev
+      real ziadvtrac,ziadvtrac2,zrec2
+      integer zim,zjm,zlm,zklon,zklev,zrec
 
       real zpi
+
 c  longitudes et latitudes lues
-      real rlonul(1:imo+1),rlatvl(1:jmo)
-      real rlonvl(1:imo+1),rlatul(1:jmo+1)
+      real rlonul(imo+1,jmo+1),rlatvl(imo+1,jmo)
+      real rlonvl(imo+1,jmo),rlatul(imo+1,jmo+1)
 c  longitudes et latitudes anciennes
       real rlonuo(0:imo+1),rlatvo(0:jmo+1)
+      real rlonvo(0:imo+1),rlatuo(0:jmo+1)
 c  longitudes et latitudes nouvelles
       real rlonun(0:imn+1),rlatvn(0:jmn+1)
+      real rlonvn(0:imn+1),rlatun(0:jmn+1)
       real aireo(imo+1,jmo+1)
 
@@ -95,5 +94,5 @@
       real alphaxo(imo+1)
       real alpha(imn+1,jmn+1)
-
+       real alphat(imn+1,jmn+1,llm)
       real aa,uu(0:imo+1),vv(imo+1,0:jmo+1)
 
@@ -105,6 +104,9 @@
       integer i,j
       real dlatm,dlatp,dlonm,dlonp
-
-
+c abd
+       character*10 file
+       character*10 nom
+       character*2 str2
+c fin ab
       zpi=2.*asin(1.)
 
@@ -161,16 +163,25 @@
 
       iest(1)=0
+	print*,'iest(1)=0'
       do io=2,imo+1
          iest(io)=iest(io-1)+ndecx(io-1)
          iouest(io-1)=iest(io)
+	print*,'iest(',io,')=',iest(io),'iouest('
+     s              ,io-1,')=',iouest(io-1)
+	
       enddo
       iouest(imo+1)=iest(imo+1)+ndecx(imo+1)
+	print*,'iouest(',imo+1,')=',iouest(imo+1)
 
       jnord(1)=0
+	print*,'jnord(1)=0'
       do jo=2,jmo+1
          jnord(jo)=jnord(jo-1)+ndecy(jo-1)
          jsud(jo-1)=jnord(jo)
+        print*,'jnord(',jo,')=',jnord(jo),'jsud('
+     s              ,jo-1,')=',jsud(jo-1)
       enddo
       jsud(jmo+1)=jnord(jmo+1)+ndecy(jmo+1)
+	        print*,'jsud(',jmo+1,')=',jsud(jmo+1)
 
 c==================================================================
@@ -182,23 +193,9 @@
         CALL read_fstoke(0,
      .   zrec,zim,zjm,zlm,
-     .   rlonu_dy,rlonv_dy,rlatu_dy,rlatv_dy,aireo,phiso,
+     .   rlonul,rlonvl,rlatul,rlatvl,aireo,phiso,
      .   massemo,pbaruo,pbarvo,wo,tetao,phio)
-
-	print*,'zrec,zdtvr,ziadvtrac,zim,zjm,zlm'
+      
+      print*,'zrec,zdtvr,ziadvtrac,zim,zjm,zlm'
       print*,zrec,zdtvr,ziadvtrac,zim,zjm,zlm
-
-      if((imo-nint(zim))*(jmo-nint(zjm)).ne.0) then
-        print*,'Modifier les dimensions dans redecoupe '
-        print*,'Mettre imo=',zim,'   jmo=',zjm
-        stop
-      endif
-
-        CALL read_pstoke(0,
-     .   zrec,zklon,zklev,airefi,phisfi,
-     .   mfu,mfd,en_u,de_u,en_d,de_d,coefkz,
-     .   frac_impa,frac_nucl,yu1,yv1,ftsol,pctsrf)
-	
-	print*,'Entete du fichier physique'
-      print*,zrec2,ziadvtrac2,zklon,zklev
 
       nrec=zrec
@@ -207,4 +204,39 @@
       istphy=ziadvtrac2
 
+      print*,'rlonul '
+      do io=1,imo+1
+         print*,io,rlonul(io,1)
+      enddo
+	print*,'rlonvl '
+      do io=1,imo+1
+         print*,io,rlonvl(io,1)
+      enddo
+	print*,'rlatul '
+      do jo=1,jmo+1
+         print*,jo,rlatul(1,jo)
+      enddo
+	print*,'rlatvl'
+      do jo=1,jmo
+         print*,jo,rlatvl(1,jo)
+      enddo
+      
+c      if((imo-nint(zim))*(jmo-nint(zjm)).ne.0) then
+c        print*,'Modifier les dimensions dans redecoupe '
+c        print*,'Mettre imo=',zim,'   jmo=',zjm
+c  abderr      stop
+c      endif
+
+c abderrahmane
+	if(physic)then
+	CALL read_pstoke(0,
+     .   zrec,zklon,zklev,airefio,phisfio,
+     .   to,mfuo,mfdo,en_uo,de_uo,en_do,de_do,coefkzo,
+     .   frac_impao,frac_nuclo,yu1o,yv1o,ftsolo,pctsrfo)
+	
+	print*,'Entete du fichier physique'
+        print*,zrec,zklon,zklev
+	endif
+
+
 c==================================================================
 c   Definition des anciennes latitudes et longitudes
@@ -212,30 +244,50 @@
 c==================================================================
 
-      rlonuo(0)=-zpi
+
       do io=1,imo
-c        rlonuo(io)=2.*zpi/FLOAT(imo)*(io+0.5-0.5*FLOAT(imo)-1.)
-c        print*,'LON ',io,rlonuo(io),rlonul(io)
-         rlonuo(io)=rlonul(io)
-      enddo
-      rlonuo(imo+1)=zpi
+         rlonuo(io)=rlonul(io,1)*zpi/180.
+	print*,'LON ',io,rlonuo(io)*180./zpi
+      enddo
+c abderr
+      rlonuo(imo+1)=0.5*(rlonul(imo,1)+rlonul(imo+1,1))*zpi/180.
+      print*,'LON ',imo+1,rlonuo(imo+1)*180./zpi
+      rlonuo(0)=rlonuo(imo+1)-2.*zpi
+      print*,'LON ',0,rlonuo(0)*180./zpi
+
+c abder
+c ATTENTION A REVOIR
+c	goto 22
+      do io=1,imo
+         rlonvo(io)=rlonvl(io,1)*zpi/180.
+        print*,'LON ',io,rlonvo(io)*180./zpi
+      enddo
+      rlonvo(imo+1)=0.5*(rlonvl(imo,1)+rlonvl(imo+1,1))*zpi/180.
+      print*,'LON ',imo+1,rlonvo(imo+1)*180./zpi
+      rlonvo(0)=rlonvo(imo+1)-2.*zpi
+      print*,'LON ',0,rlonvo(0)*180./zpi
+22    continue
+c fin ab
 
       rlatvo(0)=zpi/2.
+      print*,'LAT ',0,rlatvo(0)*180./zpi
       do jo=1,jmo
-c        rlatvo(jo)=zpi/FLOAT(jmo)*(0.5*FLOAT(jmo)+1.-jo-0.5)
-c        print*,'LAT ',jo,rlatvo(jo),rlatvl(jo)
-         rlatvo(jo)=rlatvl(jo)
-      enddo
-      rlatvo(jmo+1)=-zpi/2.
-
-c     do jo=1,jmo+1
-c        do io=1,imo+1
-c           aireo(io,jo)=rad*rad
-c    s         *(rlonuo(io)-rlonuo(io-1))
-c    s         *(sin(rlatvo(jo-1))-sin(rlatvo(jo)))
-c           aireo(io,jo)=airel(io,jo)
-c        enddo
-c        aireo(1,jo)=aireo(1,jo)+aireo(imo+1,jo)
-c        aireo(imo+1,jo)=aireo(1,jo)
-c     enddo
+        rlatvo(jo)=rlatvl(1,jo)*zpi/180.
+	print*,'LAT ',jo,rlatvo(jo)*180./zpi
+      enddo
+	rlatvo(jmo+1)=-zpi/2.
+	print*,'LAT ',jmo+1,rlatvo(jmo+1)*180./zpi
+c abd
+c ATTENTION A REVOIR
+c	goto 33
+c      rlatuo(0)=zpi/2.
+c      print*,'LAT ',0,rlatuo(0)*180./zpi
+      do jo=1,jmo+1
+        rlatuo(jo-1)=rlatul(1,jo)*zpi/180.
+        print*,'LAT ',jo-1,rlatuo(jo-1)*180./zpi
+      enddo
+        rlatuo(jmo+1)=-zpi/2.
+        print*,'LAT ',jmo+1,rlatuo(jmo+1)*180./zpi
+33    continue
+c abd
 
       do io=2,imo
@@ -250,5 +302,6 @@
 c==================================================================
 
-      rlonun(0)=-zpi
+c Nouvelles longitudes rlonun
+      rlonun(0)=rlonuo(0)
       do io=1,imo+1
          do iin=1,iouest(io)-iest(io)
@@ -259,23 +312,53 @@
             alphax(in)=alphaxo(io)/ndecx(io)
             print787,io,rlonuo(io-1)*180./zpi,in
-     s      ,iest(io),iouest(io),rlonun(in)*180./zpi,alphax(in)
-         enddo
-      enddo
-
+     s  ,iest(io),iouest(io),rlonun(in)*180./zpi,alphax(in)
+         enddo
+      enddo
+
+c Nouvelles longitudes rlonvn
+c	goto 44
+      rlonvn(0)=rlonvo(0)
+      do io=1,imo+1
+         do iin=1,iouest(io)-iest(io)
+            in=iin+iest(io)
+            rlonvn(in)=
+     s      rlonvo(io-1)+iin*(rlonvo(io)-rlonvo(io-1))
+     s      /ndecx(io)
+            alphax(in)=alphaxo(io)/ndecx(io)
+            print787,io,rlonvo(io-1)*180./zpi,in
+     s  ,iest(io),iouest(io),rlonvn(in)*180./zpi,alphax(in)
+         enddo
+      enddo
+44    continue
+
+c Nouvelles latitudes rlatvn
       rlatvn(0)=0.5*zpi
       do jo=1,jmo+1
-         print*,'jo=',jo
          do jjn=1,jsud(jo)-jnord(jo)
             jn=jnord(jo)+jjn
-            rlatvn(jn)=rlatvo(jo-1)+jjn*(rlatvo(jo)-rlatvo(jo-1))
+            rlatvn(jn)=
+     s      rlatvo(jo-1)+jjn*(rlatvo(jo)-rlatvo(jo-1))
      s      /ndecy(jo)
             alphay(jn)=(sin(rlatvn(jn-1))-sin(rlatvn(jn)))
      s                /(sin(rlatvo(jo-1))-sin(rlatvo(jo)))
-            print787,jo,rlatvo(jo-1)*180./zpi,jn
-     s      ,jnord(jo),jsud(jo),rlatvn(jn)*180./zpi,alphay(jn)
-         enddo
-      enddo
-
-787   format(i5,f10.2,3(i5),2(f10.2))
+           print*,jn,rlatvn(jn)*180./zpi
+         enddo
+      enddo
+
+c Nouvelles latitudes rlatun
+c	goto 55
+      rlatun(0)=0.5*zpi
+      do jo=1,jmo+1
+         do jjn=1,jsud(jo)-jnord(jo)
+            jn=jnord(jo)+jjn
+            rlatun(jn)=
+     s      rlatuo(jo-1)+jjn*(rlatuo(jo)-rlatuo(jo-1))
+     s      /ndecy(jo)
+           print*,jn,rlatvn(jn)*180./zpi
+         enddo
+      enddo
+55    continue
+
+787   format(i5,f10.2,3(i5),2(f12.6))
       do in=1,imn
          rlonu(in)=rlonun(in)
@@ -295,8 +378,17 @@
          do in=1,imn
             alpha(in,jn)=alphax(in)*alphay(jn)
+	    alphat(in,jn,1)=alpha(in,jn)
          enddo
          alpha(imn+1,jn)=0.
-      enddo
-
+	 alphat(imn+1,jn,1)=0.
+      enddo
+c abderr 19 4 00
+	do l=2,llm
+	 do jn=1,jmn+1
+         do in=1,imn+1
+	  alphat(in,jn,l)=alphat(in,jn,1)
+         enddo
+	 enddo
+        enddo
 c     call dump2d(iip1,jjp1,alpha,'ALPHA   ')
 
@@ -358,14 +450,4 @@
 c     call dump2d(iip1,jjp1,aire,'AIRE   ')
 
-c     do jn=1,jjp1
-c        do in=1,iim
-c           aire(in,jn)=rad*rad*(sin(rlatvn(jn-1))-sin(rlatvn(jn)))
-c    s      *(rlonun(in)-rlonun(in-1))
-c           unsaire(in,jn)=1./aire(in,jn)
-c        enddo
-c        aire(iip1,jn)=aire(1,jn)
-c        unsaire(iip1,jn)=unsaire(1,jn)
-c     enddo
-c     call dump2d(iip1,jjp1,aire,'AIRE2   ')
       DO 42 j = 1,jjp1
       DO 41 i = 1,iim
@@ -404,6 +486,17 @@
          enddo
       enddo
-
-
+	Print*,'Fin irec=0'
+	go to 435
+      file='pbur'
+      call inigrads(11,iip1
+     s  ,rlonu,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,1800.,file,'gcmq2 ')
+      file='pbvr'
+      call inigrads(12,iip1
+     s  ,rlonv,180./pi,-180.,180.,jjm,rlatv,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,1800.,file,'gcmq2 ')
+435   continue
 c==================================================================
 c   Fin des initialisations
@@ -415,17 +508,19 @@
 c   Lecture des fichiers fluxmass et  physique:
 c   -----------------------------------------------------
-	
+	print*,'Entrer dans read_fstoke a irec=',irec
 	CALL read_fstoke(irec,
      .   zrec,zim,zjm,zlm,
-     .   rlonu_dy,rlonv_dy,rlatu_dy,rlatv_dy,aireo,phiso,
+     .   rlonul,rlonvl,rlatul,rlatvl,aireo,phiso,
      .   massemo,pbaruo,pbarvo,wo,tetao,phio)
 
-	do l=1,llm
-           do j=1,jmo
-              do i=1,imo+1
-                 pbarvo(i,j,l)=pbarvst(i,j,l)
-              enddo
-           enddo
-        enddo
+	print*,'Apres read_fstoke a irec=',irec
+
+c	do l=1,llm
+c           do j=1,jmo
+c              do i=1,imo+1
+c                 pbarvo(i,j,l)=pbarvst(i,j,l)
+c              enddo
+c           enddo
+c        enddo
 
          do l=1,llm
@@ -434,11 +529,11 @@
                   do jn=jnord(jo)+1,jsud(jo)
                      do in=iest(io)+1,iouest(io)
-                        wn(in,jn,l)=alpha(in,jn)*wo(io,jo,l)
-                        massemn(in,jn,l)=alpha(in,jn)
-     s                     *massemo(io,jo,l)
+c                        wn(in,jn,l)=alpha(in,jn)*wo(io,jo,l)
+c                        massemn(in,jn,l)=alpha(in,jn)
+		       wn(in,jn,l)=alphat(in,jn,l)*wo(io,jo,l)
+		       massemn(in,jn,l)=alphat(in,jn,l)
+     s                                *massemo(io,jo,l)
                         tetan(in,jn,l)=tetao(io,jo,l)
                         phin(in,jn,l)=phio(io,jo,l)
-c marine                
-                        phisn(i,jn) = phiso(io,jo)
                      enddo
                   enddo
@@ -450,10 +545,23 @@
                tetan(imn+1,jn,l)=tetan(1,jn,l)
                phin(imn+1,jn,l)=phin(1,jn,l)
-c marine
-                phisn(imn+1,jn)=phisn(1,jn)
-
-            enddo
-         enddo
-
+            enddo
+         enddo
+c Test massemn
+	print*,'MASSE DANS LA NOUVELLE GRILLE'
+	goto 908
+	do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+	print*,'massemn(',in,jn,1,')=',massemn(in,jn,1)
+        	     enddo
+                  enddo
+               enddo
+	enddo
+	do jn=1,jmn+1
+	print*,'massemn(',imn+1,jn,1,')=',massemn(imn+1,jn,1)
+        enddo
+908     continue
+	print*,'Fin calcul de massemn pour nouv. gril.'
          do l=1,llm
             do jo=1,jmo+1
@@ -478,4 +586,5 @@
             enddo
          enddo
+	print*,'Fin calcul de pbarun'
 	
 	do l=1,llm
@@ -507,56 +616,508 @@
          enddo
 
-	
- 	CALL read_pstoke(irec,
-     .   zrec,zklon,zklev,airefi,phisfi,
-     .   mfu,mfd,en_u,de_u,en_d,de_d,coefkz,
-     .   frac_impa,frac_nucl,yu1,yv1,ftsol,pctsrf)
-
+c abd
+	go to 456
+ 	nom='pbaru'
+        call wrgrads(11,llm,pbarun(:,:,1),nom,nom)
+        nom='pbarv'
+        call wrgrads(12,llm,pbarvn(:,:,1),nom,nom)
+        nom='masse'
+        call wrgrads(11,llm,massemn(:,:,1),nom,nom)
+        nom='w'
+        call wrgrads(11,llm,wn(:,:,1),nom,nom)
+456	continue
+c fin ab
+
+	if(physic)then
+        CALL read_pstoke(irec,
+     .   zrec,zklon,zklev,airefio,phisfio,
+     .   to,mfuo,mfdo,en_uo,de_uo,en_do,de_do,coefkzo,
+     .   frac_impao,frac_nuclo,yu1o,yv1o,ftsolo,pctsrfo)
+       	print*,'OK read_pstoke pour irec=',irec
 c==================================================================
 c  Passage  a la nouvelle grille
 c==================================================================
-         do l=1,9*llm+2+2*nbsrf
-c   passage aa la grille dynamique ancienne
-            do io=1,imo+1
-               tmpo(io,1)=xleco(1,l)
-               tmpo(io,jmo+1)=xleco(ngrido,l)
-            enddo
-            do jo=2,jmo
-               do io=1,imo
-                  tmpo(io,jo)=xleco((jo-2)*imo+io+1,l)
-               enddo
-               tmpo(imo+1,jo)=tmpo(1,jo)
-            enddo
-c   passage a la grillle dynamique nouvelle
-            do jo=1,jmo+1
-               do io=1,imo+1
-                  do jn=jnord(jo)+1,jsud(jo)
-                     do in=iest(io)+1,iouest(io)
-                        tmpn(in,jn)=tmpo(io,jo)
-                     enddo
-                  enddo
-               enddo
-            enddo
-c   passage a la grille physique nouvelle
-            xlecn(1,l)=tmpn(1,1)
-            xlecn(ngridn,l)=tmpn(1,jmn+1)
-            do jn=2,jmn
-               do in=1,imn
-                  xlecn((jn-2)*imn+in+1,l)=tmpn(in,jn)
-               enddo
-            enddo
-         enddo
-
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,to,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=to(1,l)
+                tmpo2(imo+1,jmo+1,l)=to(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,tn)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+
+	call gr_fi_ecrit(llm,ngrido,imo,jmo+1,mfuo,tmpo1)
+	do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+		tmpo2(imo+1,1,l)=mfuo(1,l)
+                tmpo2(imo+1,jmo+1,l)=mfuo(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+		tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+	enddo
+	call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,mfun)
+	call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+	call initial0(llm*imo*(jmo+1),tmpo1)
+	call initial0(llm*imn*(jmn+1),tmpn1)
+	call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,mfdo,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=mfdo(1,l)
+                tmpo2(imo+1,jmo+1,l)=mfdo(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,mfdn)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,en_uo,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=en_uo(1,l)
+                tmpo2(imo+1,jmo+1,l)=en_uo(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,en_un)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,en_do,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=en_do(1,l)
+                tmpo2(imo+1,jmo+1,l)=en_do(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,en_dn)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,de_do,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=de_do(1,l)
+                tmpo2(imo+1,jmo+1,l)=de_do(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,de_dn)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,de_uo,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=de_uo(1,l)
+                tmpo2(imo+1,jmo+1,l)=de_uo(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,de_un)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,coefkzo,tmpo1)
+	do l=1,llm
+	   do jo=1,jmo+1
+	      do io=1,imo
+		tmpo2(io,jo,l)=tmpo1(io,jo,l)
+	      enddo
+	   enddo
+
+           tmpo2(imo+1,1,l)=coefkzo(1,l)
+           tmpo2(imo+1,jmo+1,l)=coefkzo(ngrido,l)
+
+           do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo2(1,jo,l)
+           enddo
+	enddo
+	      
+c   passage a la grillle dynamique nouvelle
+	do l=1,llm	
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,coefkzn)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,frac_impao,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=frac_impao(1,l)
+                tmpo2(imo+1,jmo+1,l)=frac_impao(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,frac_impan)
+        call initial0(llm*(imo+1)*(jmo+1),tmpo2)
+        call initial0(llm*imo*(jmo+1),tmpo1)
+        call initial0(llm*imn*(jmn+1),tmpn1)
+        call initial0(llm*(imn+1)*(jmn+1),tmpn2)
+        call gr_fi_ecrit(llm,ngrido,imo,jmo+1,frac_nuclo,tmpo1)
+        do l=1,llm
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo2(io,jo,l)=tmpo1(io,jo,l)
+              enddo
+           enddo
+
+                tmpo2(imo+1,1,l)=frac_nuclo(1,l)
+                tmpo2(imo+1,jmo+1,l)=frac_nuclo(ngrido,l)
+            do jo=2,jmo
+                tmpo2(imo+1,jo,l)=tmpo1(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn2(in,jn,l)=tmpo2(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn1(in,jn,l)=tmpn2(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(llm,ngridn,imn,jmn+1,tmpn1,frac_nucln)
+
+        call gr_fi_ecrit(nbsrf,ngrido,imo,jmo+1,ftsolo,tmpo3)
+        do l=1,nbsrf
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo4(io,jo,l)=tmpo3(io,jo,l)
+              enddo
+           enddo
+
+                tmpo4(imo+1,1,l)=ftsolo(1,l)
+                tmpo4(imo+1,jmo+1,l)=ftsolo(ngrido,l)
+            do jo=2,jmo
+                tmpo4(imo+1,jo,l)=tmpo3(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn4(in,jn,l)=tmpo3(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn3(in,jn,l)=tmpn4(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(nbsrf,ngridn,imn,jmn+1,tmpn3,ftsoln)
+
+        call initial0(nbsrf*(imo+1)*(jmo+1),tmpo4)
+        call initial0(nbsrf*imo*(jmo+1),tmpo3)
+        call initial0(nbsrf*imn*(jmn+1),tmpn3)
+        call initial0(nbsrf*(imn+1)*(jmn+1),tmpn4)
+        call gr_fi_ecrit(nbsrf,ngrido,imo,jmo+1,pctsrfo,tmpo3)
+        do l=1,nbsrf
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo4(io,jo,l)=tmpo3(io,jo,l)
+              enddo
+           enddo
+
+                tmpo4(imo+1,1,l)=pctsrfo(1,l)
+                tmpo4(imo+1,jmo+1,l)=pctsrfo(ngrido,l)
+            do jo=2,jmo
+                tmpo4(imo+1,jo,l)=tmpo3(1,jo,l)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn4(in,jn,l)=tmpo3(io,jo,l)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn3(in,jn,l)=tmpn4(in,jn,l)
+               enddo
+            enddo
+        enddo
+        call gr_ecrit_fi(nbsrf,ngridn,imn,jmn+1,tmpn3,pctsrfn)
+
+        call gr_fi_ecrit(1,ngrido,imo,jmo+1,yv1o,tmpo5)
+
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo6(io,jo)=tmpo5(io,jo)
+              enddo
+           enddo
+
+                tmpo6(imo+1,1)=yv1o(1)
+                tmpo6(imo+1,jmo+1)=yv1o(ngrido)
+            do jo=2,jmo
+                tmpo6(imo+1,jo)=tmpo5(1,jo)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn6(in,jn)=tmpo5(io,jo)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn5(in,jn)=tmpn6(in,jn)
+               enddo
+            enddo
+        call gr_ecrit_fi(1,ngridn,imn,jmn+1,tmpn5,yv1n)
+
+        call initial0((imo+1)*(jmo+1),tmpo6)
+        call initial0(imo*(jmo+1),tmpo5)
+        call initial0(imn*(jmn+1),tmpn5)
+        call initial0((imn+1)*(jmn+1),tmpn6)
+        call gr_fi_ecrit(1,ngrido,imo,jmo+1,yu1o,tmpo5)
+
+           do jo=1,jmo+1
+              do io=1,imo
+                tmpo6(io,jo)=tmpo5(io,jo)
+              enddo
+           enddo
+
+                tmpo6(imo+1,1)=yu1o(1)
+                tmpo6(imo+1,jmo+1)=yu1o(ngrido)
+            do jo=2,jmo
+                tmpo6(imo+1,jo)=tmpo5(1,jo)
+            enddo
+c   passage a la grillle dynamique nouvelle
+            do jo=1,jmo+1
+               do io=1,imo+1
+                  do jn=jnord(jo)+1,jsud(jo)
+                     do in=iest(io)+1,iouest(io)
+                        tmpn6(in,jn)=tmpo5(io,jo)
+                     enddo
+                  enddo
+               enddo
+            enddo
+            do jn=1,jmn+1
+               do in=1,imn
+                tmpn5(in,jn)=tmpn6(in,jn)
+               enddo
+            enddo
+        call gr_ecrit_fi(1,ngridn,imn,jmn+1,tmpn5,yu1n)
 c==================================================================
 	if (avant) then
 c Simu directe
-       do l=1,llm
+         do l=1,llm
           do ig=1,ngridn
-             zmfu(ig,l)=mfu(ig,l)
-             zmfd(ig,l)=mfd(ig,l)
-             zde_u(ig,l)=de_u(ig,l)
-             zen_u(ig,l)=en_u(ig,l)
-             zde_d(ig,l)=de_d(ig,l)
-             zen_d(ig,l)=en_d(ig,l)
+	     zmfd(ig,l)=mfdn(ig,l)
+             zmfu(ig,l)=mfun(ig,l)
+             zen_d(ig,l)=en_dn(ig,l)
+             zde_d(ig,l)=de_dn(ig,l)
+             zen_u(ig,l)=en_un(ig,l)
+             zde_u(ig,l)=de_un(ig,l)
           enddo
        enddo
@@ -565,10 +1126,10 @@
        do l=1,llm
           do ig=1,ngridn
-             zmfd(ig,l)=-mfu(ig,l)
-             zmfu(ig,l)=-mfd(ig,l)
-             zen_d(ig,l)=de_u(ig,l)
-             zde_d(ig,l)=en_u(ig,l)
-             zen_u(ig,l)=de_d(ig,l)
-             zde_u(ig,l)=en_d(ig,l)
+             zmfd(ig,l)=-mfdn(ig,l)
+             zmfu(ig,l)=-mfun(ig,l)
+             zen_d(ig,l)=en_dn(ig,l)
+             zde_d(ig,l)=de_dn(ig,l)
+             zen_u(ig,l)=en_un(ig,l)
+             zde_u(ig,l)=de_un(ig,l)
           enddo
        enddo
@@ -585,8 +1146,8 @@
             zcontrole(ig)=1.
          enddo
-c   zmass=(max(massemn(ig,l),massemn(ig,l-1))/airefi(ig)
+c   zmass=(max(massemn(ig,l),massemn(ig,l-1))/airefin(ig)
          do l=2,llm
             do ig=1,ngridn
-               zmass=max(massefi(ig,l),massefi(ig,l-1))/airefi(ig)
+               zmass=max(massefi(ig,l),massefi(ig,l-1))/airefin(ig)
                zflux=max(abs(zmfu(ig,l)),abs(zmfd(ig,l)))*dtphys
                if(zflux.gt.0.9*zmass) then
@@ -615,5 +1176,5 @@
             enddo
          enddo
-
+	endif ! physic
 
       endif ! irec=0
Index: DZ.3.3/trunk/libf/phylmd/phystoke.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/phystoke.F	(revision 203)
+++ 	(revision )
@@ -1,214 +1,0 @@
-      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
Index: /LMDZ.3.3/trunk/libf/phylmd/phystokenc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/phystokenc.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/phylmd/phystokenc.F	(revision 204)
@@ -6,5 +6,5 @@
      I                   pt,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
      I                   pcoefh,yu1,yv1,ftsol,pctsrf,
-     I                   frac_impa,frac_nucl,
+     I                   pfrac_impa,pfrac_nucl,
      I                   pphis,paire,dtime,itap)
       USE ioipsl
@@ -63,6 +63,6 @@
 c   ----------
 c
-      REAL frac_impa(klon,klev)
-      REAL frac_nucl(klon,klev)
+      REAL pfrac_impa(klon,klev)
+      REAL pfrac_nucl(klon,klev)
 c
 c Arguments necessaires pour les sources et puits de traceur
@@ -82,4 +82,7 @@
       REAL coefh(klon,klev) ! flux detraine dans le panache descendant
 	REAL t(klon,klev)
+      REAL frac_impa(klon,klev)
+      REAL frac_nucl(klon,klev)
+      REAL rain(klon)
 
       REAL pyu1(klon),pyv1(klon)
@@ -95,4 +98,5 @@
       save t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,dtcum
       save iadvtr,irec
+      save frac_impa,frac_nucl,rain
       save pyu1,pyv1,pftsol,ppsrf
 
@@ -133,4 +137,6 @@
          do k=1,klev
             do i=1,klon
+               frac_impa(i,k)=1.
+               frac_nucl(i,k)=1.
                mfu(i,k)=0.
                mfd(i,k)=0.
@@ -144,4 +150,5 @@
          enddo
          do i=1,klon
+            rain(i)=0.
             pyv1(i)=0.
             pyu1(i)=0.
@@ -159,4 +166,6 @@
       do k=1,klev
          do i=1,klon
+            frac_impa(i,k)=frac_impa(i,k)*pfrac_impa(i,k)
+            frac_nucl(i,k)=frac_nucl(i,k)*pfrac_nucl(i,k)
             mfu(i,k)=mfu(i,k)+pmfu(i,k)*pdtphys
             mfd(i,k)=mfd(i,k)+pmfd(i,k)*pdtphys
@@ -187,4 +196,6 @@
          do k=1,klev
             do i=1,klon
+c              frac_impa=frac_impa : c'est la fraction cumulee qu'on stoke
+c              frac_nucl=frac_nucl : c'est la fraction cumulee qu'on stoke
                mfu(i,k)=mfu(i,k)/dtcum
                mfd(i,k)=mfd(i,k)/dtcum
@@ -198,4 +209,5 @@
          enddo
          do i=1,klon
+            rain(i)=rain(i)/dtcum
             pyv1(i)=pyv1(i)/dtcum
             pyu1(i)=pyu1(i)/dtcum
@@ -226,7 +238,9 @@
 
 ccccc
+      print*,'AVANT ECRITURE'
          CALL gr_fi_ecrit(klev,klon,iim,jjm+1, t, zx_tmp_3d)
          CALL histwrite(physid,"t",itap,zx_tmp_3d,
      .                                   iim*(jjm+1)*klev,ndex)
+      print*,'APRES ECRITURE'
 
          CALL gr_fi_ecrit(klev,klon,iim,jjm+1, mfu, zx_tmp_3d)
@@ -279,4 +293,8 @@
          CALL gr_fi_ecrit(1,klon,iim,jjm+1, pftsol4, zx_tmp_2d)
       CALL histwrite(physid,"ftsol4",itap,zx_tmp_2d,
+     .                                   iim*(jjm+1),ndex)
+
+        CALL gr_fi_ecrit(1,klon,iim,jjm+1, rain, zx_tmp_2d)
+      CALL histwrite(physid,"rain",itap,zx_tmp_2d,
      .                                   iim*(jjm+1),ndex)
 
Index: /LMDZ.3.3/trunk/libf/phylmd/phytrac.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/phytrac.F	(revision 203)
+++ /LMDZ.3.3/trunk/libf/phylmd/phytrac.F	(revision 204)
@@ -30,6 +30,4 @@
 #include "dimphy.h"
 #include "indicesol.h"
-#include "control.h"
-#include "temps.h"
 c======================================================================
 
@@ -51,5 +49,4 @@
       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)
@@ -95,5 +92,8 @@
       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
@@ -136,5 +136,6 @@
       INTEGER nid_tra
       SAVE nid_tra
-      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
+c     REAL x(klon,klev,nbtr+2) ! traceurs 
+      INTEGER ndex(1)
       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)
@@ -164,6 +165,6 @@
 c
 c--modif convection tiedtke
-      INTEGER i, k, it
-
+      INTEGER i, k, it,itap
+	save itap
       REAL delp(klon,klev)
 c--end modif
@@ -211,8 +212,9 @@
 c        print*,'DANS PHYTRAC debutphy=',debutphy
 
-         ecrit_tra = NINT(86400./pdtphys *ecritphy)   
-         zsto = pdtphys
-         zout = pdtphys * FLOAT(ecrit_tra)
          if (debutphy) then
+
+          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
+          ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
+c         ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
 
          if(nbtr.lt.nqmax) then
@@ -226,4 +228,5 @@
          PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
          itra=0
+	 itap=0
 C         
          CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
@@ -242,9 +245,8 @@
      .                 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)
+         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
+     .                 klev, presnivs, nvert)
+         zsto = pdtphys
+         zout = pdtphys * FLOAT(ecrit_tra)
 c
          CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
@@ -255,4 +257,66 @@
      .                iim,jjm+1,nhori, 1,1,1, -99, 32,
      .                "once",  zsto,zout)
+
+	goto 666
+         CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+
+         CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "psrf1", "nature sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "psrf2", "nature sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "psrf3", "nature sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "psrf4", "nature sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst",  zsto,zout)
+         CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
+     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
+     .                "inst(X)",  zsto,zout)
+         CALL histdef(nid_tra, "pplay", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "inst(X)", zsto,zout)
+         CALL histdef(nid_tra, "t", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "inst(X)", zsto,zout)
+         CALL histdef(nid_tra, "mfu", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "mfd", "flux u decen","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "en_u", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "en_d", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "de_u", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "de_d", "flux u mont","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+         CALL histdef(nid_tra, "coefh", "turbulent coef","-",
+     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
+     .                "ave(X)", zsto,zout)
+
+666	continue
 c
          DO it=1,nqmax
@@ -274,4 +338,13 @@
          ENDDO
          CALL histend(nid_tra)
+         ndex(1) = 0
+c
+         i = NINT(zout/zsto)
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
+         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
+C
+         i = NINT(zout/zsto)
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
+         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
 
 c======================================================================
@@ -327,11 +400,75 @@
          inirnpb=.false.
       endif
+      if(nqmax.gt.2) aerosol(3)=.true.
+
+
+c  abder
+	goto 777
+            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
+	 ndex(1)=0
+	 itap=itap+1
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
+         CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
+     s                                  iim*(jjm+1),ndex)
+         
+	 CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
+         CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
+     s                                  iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
+         CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
+     s                                       iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
+         CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
+     s                                       iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
+         CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
+     s                                      iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
+         CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
+     s                                      iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
+         CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
+     s                                     iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
+         CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d, 
+     s                                     iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
+         CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
+     s                                     iim*(jjm+1),ndex)
+
+         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
+         CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
+     s                                     iim*(jjm+1),ndex)
+777	continue
 c======================================================================
 c   Calcul de l'effet de la convection
 c======================================================================
+	print*,'Avant convection'
+      do it=1,nqmax
+         WRITE(itn,'(i1)') it
+c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
+      enddo
 
       if (convection) then
 
-c     print*,'Pas de temps dans phytrac : ',pdtphys
+      print*,'Pas de temps dans phytrac : ',pdtphys
       DO it=1, nqmax
       CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
@@ -342,17 +479,26 @@
       ENDDO
       ENDDO
-      WRITE(itn,'(i1)') it
-      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
-      ENDDO
-c     print*,'apres nflxtr'
+c      WRITE(itn,'(i1)') it
+c      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
+      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 de la couche limite
 c======================================================================
-
-c     print*,'avant couchelimite'
+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
 
@@ -413,5 +559,9 @@
       endif ! couche limite
 
-c     print*,'apres couchelimite'
+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======================================================================
@@ -442,4 +592,5 @@
 c======================================================================
 
+      print*,'LESSIVAGE =',lessivage
       IF (lessivage) THEN
 
@@ -474,14 +625,20 @@
 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) - 1. )    
+               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 
@@ -517,16 +674,5 @@
       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)
+      ndex(1) = 0
       DO it=1,nqmax
       IF (it.LE.99) THEN
@@ -535,10 +681,10 @@
        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
+     .                                   iim*(jjm+1)*klev,ndex)
+c      IF (lessivage) THEN
+c      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
+c      CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
+c    .                                   iim*(jjm+1)*klev,ndex)
+c     ENDIF
       ELSE
          PRINT*, "Trop de traceurs"
@@ -546,5 +692,42 @@
       ENDIF
       ENDDO
-      if (ok_sync) call histsync(nid_tra)
+
+	goto 888
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
+        CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
+        CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
+        CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
+        CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
+        CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
+        CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
+        CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
+        CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
+        CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
+     .                  iim*(jjm+1)*klev,ndex)
+
+888	continue
+
+c	print*,'Sortie phytrac'
+c      do it=1,nqmax
+c         WRITE(itn,'(i1)') it
+c        call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys  '//itn)
+c      enddo
 
       if (lafin) then
