Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/calfis.F	(revision 232)
@@ -319,13 +319,13 @@
 c
       DO l=1,llm
-        pvervel(1,l)=pw(1,1,l)/apoln
+        pvervel(1,l)=pw(1,1,l) * g /apoln
         ig0=2
        DO j=2,jjm
            DO i = 1, iim
-              pvervel(ig0,l) = pw(i,j,l) * unsaire(i,j)
+              pvervel(ig0,l) = pw(i,j,l) * g * unsaire(i,j)
               ig0 = ig0 + 1
            ENDDO
        ENDDO
-        pvervel(ig0,l)=pw(1,jjp1,l)/apols
+        pvervel(ig0,l)=pw(1,jjp1,l) * g /apols
       ENDDO
 
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/comvert.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/comvert.h	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/comvert.h	(revision 232)
@@ -3,5 +3,5 @@
 
       COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm) ,
-     ,               pa,preff,nivsigs(llm),nivsig(llmp1)
+     ,               pa,preff,nivsigs(llm),nivsig(llm+1)
 
       REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/create_limit.F	(revision 232)
@@ -305,5 +305,4 @@
         ENDIF 
       END DO 
-
 c
 c
@@ -338,5 +337,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -354,5 +357,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -370,5 +377,9 @@
         STOP
       ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -387,5 +398,9 @@
          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
          print*,dimfirst,dimlast
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
          if (ierr.ne.0) then
            print *, NF_STRERROR(ierr)
@@ -471,5 +486,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -487,5 +506,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -503,5 +526,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -519,5 +546,9 @@
 c
          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
          if (ierr.ne.0) then
            print *, NF_STRERROR(ierr)
@@ -698,5 +729,10 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
+#endif
+
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -714,5 +750,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -730,5 +770,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -746,5 +790,9 @@
 c
          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
          if (ierr.ne.0) then
            print *, NF_STRERROR(ierr)
@@ -825,5 +873,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -841,5 +893,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -857,5 +913,9 @@
         STOP
       ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
       ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
       if (ierr.ne.0) then
         print *, NF_STRERROR(ierr)
@@ -873,5 +933,9 @@
 c
          PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
          ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
          if (ierr.ne.0) then
            print *, NF_STRERROR(ierr)
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/defrun_new.F	(revision 232)
@@ -1,2 +1,5 @@
+c
+c $Header$
+c
       SUBROUTINE defrun_new( tapedef, etatinit, clesphy0 )
 c
@@ -36,5 +39,5 @@
       INTEGER   tapeout
       REAL clonn,clatt,grossismxx,grossismyy
-      REAL dzoomxx,dzoomyy
+      REAL dzoomxx,dzoomyy,tauxx,tauyy
       LOGICAL  fxyhypbb, ysinuss
       INTEGER i
@@ -258,9 +261,8 @@
       WRITE(tapeout,*)    clonn
       IF( ABS(clon - clonn).GE. 0.001 )  THEN
-        PRINT *,' La valeur de clon passee par run.def est differente de 
-     *  celle lue sur le fichier  start '
+       WRITE(tapeout,*) ' La valeur de clon passee par run.def est diffe
+     *rente de  celle lue sur le fichier  start '
         STOP
       ENDIF
-c
 c
       READ (tapedef,9001) ch1,ch4
@@ -270,6 +272,6 @@
 
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
-        PRINT *,' La valeur de clat passee par run.def est differente de 
-     *  celle lue sur le fichier  start '
+       WRITE(tapeout,*) ' La valeur de clat passee par run.def est diffe
+     *rente de  celle lue sur le fichier  start '
         STOP
       ENDIF
@@ -281,6 +283,6 @@
 
       IF( ABS(grossismx - grossismxx).GE. 0.001 )  THEN
-        PRINT *,' La valeur de grossismx passee par run.def est differente  
-     *  de celle lue sur le fichier  start '
+       WRITE(tapeout,*) ' La valeur de grossismx passee par run.def est
+     , differente de celle lue sur le fichier  start '
         STOP
       ENDIF
@@ -292,11 +294,11 @@
 
       IF( ABS(grossismy - grossismyy).GE. 0.001 )  THEN
-        PRINT *,' La valeur de grossismy passee par run.def est differen
-     * te de celle lue sur le fichier  start '
+       WRITE(tapeout,*) ' La valeur de grossismy passee par run.def est
+     , differente de celle lue sur le fichier  start '
         STOP
       ENDIF
       
       IF( grossismx.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+        WRITE(tapeout,*) ' ***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -306,5 +308,5 @@
 
       IF( grossismy.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+        WRITE(tapeout,*) ' ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -312,5 +314,4 @@
       ENDIF
 
-      PRINT *,' alphax alphay defrun ',alphax,alphay
 c
 c    alphax et alphay sont les anciennes formulat. des grossissements
@@ -324,14 +325,14 @@
       IF( .NOT.fxyhypb )  THEN
            IF( fxyhypbb )     THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' *** fxyhypb lu sur le fichier start est F ',
-     *       'alors  qu il est  T  sur  run.def  ***'
+            WRITE(tapeout,*) ' ********  PBS DANS  DEFRUN  ******** '
+            WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est F'
+     *,      '                   alors  qu il est  T  sur  run.def  ***'
               STOP
            ENDIF
       ELSE
            IF( .NOT.fxyhypbb )   THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' ***  fxyhypb lu sur le fichier start est T ',
-     *        'alors  qu il est  F  sur  run.def  ****  '
+            WRITE(tapeout,*) ' ********  PBS DANS  DEFRUN  ******** '
+            WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est t'
+     *,      '                   alors  qu il est  F  sur  run.def  ***'
               STOP
            ENDIF
@@ -343,12 +344,4 @@
       WRITE(tapeout,*)    dzoomxx
 
-      IF( fxyhypb )  THEN
-       IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
-        PRINT *,' La valeur de dzoomx passee par run.def est differente 
-     *  de celle lue sur le fichier  start '
-        STOP
-       ENDIF
-      ENDIF
-
       READ (tapedef,9001) ch1,ch4
       READ (tapedef,*)    dzoomyy
@@ -356,10 +349,40 @@
       WRITE(tapeout,*)    dzoomyy
 
+      READ (tapedef,9001) ch1,ch4
+      READ (tapedef,*)    tauxx
+      WRITE(tapeout,9001) ch1,'taux'
+      WRITE(tapeout,*)    tauxx
+
+      READ (tapedef,9001) ch1,ch4
+      READ (tapedef,*)    tauyy
+      WRITE(tapeout,9001) ch1,'tauy'
+      WRITE(tapeout,*)    tauyy
+
       IF( fxyhypb )  THEN
+
+       IF( ABS(dzoomx - dzoomxx).GE. 0.001 )  THEN
+        WRITE(tapeout,*)' La valeur de dzoomx passee par run.def est dif
+     *ferente de celle lue sur le fichier  start '
+        CALL ABORT
+       ENDIF
+
        IF( ABS(dzoomy - dzoomyy).GE. 0.001 )  THEN
-        PRINT *,' La valeur de dzoomy passee par run.def est differente 
-     * de celle lue sur le fichier  start '
-        STOP
+        WRITE(tapeout,*)' La valeur de dzoomy passee par run.def est dif
+     *ferente de celle lue sur le fichier  start '
+        CALL ABORT
        ENDIF
+
+       IF( ABS(taux - tauxx).GE. 0.001 )  THEN
+        WRITE(6,*)' La valeur de taux passee par run.def est differente
+     *  de celle lue sur le fichier  start '
+        CALL ABORT
+       ENDIF
+
+       IF( ABS(tauy - tauyy).GE. 0.001 )  THEN
+        WRITE(6,*)' La valeur de tauy passee par run.def est differente
+     *  de celle lue sur le fichier  start '
+        CALL ABORT
+       ENDIF
+
       ENDIF
       
@@ -374,14 +397,14 @@
         IF( .NOT.ysinus )  THEN
            IF( ysinuss )     THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' *** ysinus lu sur le fichier start est F ',
-     *       'alors  qu il est  T  sur  run.def  ***'
+              WRITE(6,*) ' ********  PBS DANS  DEFRUN  ******** '
+              WRITE(tapeout,*)'** ysinus lu sur le fichier start est F',
+     *       ' alors  qu il est  T  sur  run.def  ***'
               STOP
            ENDIF
         ELSE
            IF( .NOT.ysinuss )   THEN
-              PRINT *,' ********  PBS DANS  DEFRUN  ******** '
-              PRINT *,' ***  ysinus lu sur le fichier start est T ',
-     *        'alors  qu il est  F  sur  run.def  ****  '
+              WRITE(6,*) ' ********  PBS DANS  DEFRUN  ******** '
+              WRITE(tapeout,*)'** ysinus lu sur le fichier start est T',
+     *       ' alors  qu il est  F  sur  run.def  ***'
               STOP
            ENDIF
@@ -389,5 +412,8 @@
       ENDIF
 c
-      close(tapedef)
+      WRITE(6,*) ' alphax alphay defrun ',alphax,alphay
+
+      CLOSE(tapedef)
+
       RETURN
 c   ...............................................
@@ -416,5 +442,5 @@
 
       IF( grossismx.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismx < 1 .   *** '
+        WRITE(tapeout,*) '***  ATTENTION !! grossismx < 1 .   *** '
          STOP
       ELSE
@@ -423,5 +449,5 @@
 
       IF( grossismy.LT.1. )  THEN
-        PRINT *,' ***  ATTENTION !! grossismy < 1 .   *** '
+        WRITE(tapeout,*) ' ***  ATTENTION !! grossismy < 1 .   *** '
          STOP
       ELSE
@@ -429,6 +455,4 @@
       ENDIF
 
-      PRINT *,' alphax alphay defrun ',alphax,alphay
-
 c
       READ (tapedef,9001) ch1,ch4
@@ -446,5 +470,15 @@
       WRITE(tapeout,9001) ch1,'dzoomy'
       WRITE(tapeout,*)    dzoomy
-c
+
+      READ (tapedef,9001) ch1,ch4
+      READ (tapedef,*)    taux
+      WRITE(tapeout,9001) ch1,'taux'
+      WRITE(tapeout,*)    taux
+c
+      READ (tapedef,9001) ch1,ch4
+      READ (tapedef,*)    tauy
+      WRITE(tapeout,9001) ch1,'tauy'
+      WRITE(tapeout,*)    tauy
+
       READ (tapedef,9001) ch1,ch4
       READ (tapedef,*)    ysinus
@@ -452,9 +486,11 @@
       WRITE(tapeout,*)    ysinus
        
+      WRITE(tapeout,*) ' alphax alphay defrun ',alphax,alphay
 c
 9000  FORMAT(3(/,a72))
 9001  FORMAT(/,a72,/,a12)
 cc
-      close(tapedef)
+      CLOSE(tapedef)
+
       RETURN
       END
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynetat0.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynetat0.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynetat0.F	(revision 232)
@@ -1,2 +1,5 @@
+c
+c $Header$
+c
       SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
      .                    teta,q,masse,ps,phis,time)
@@ -105,4 +108,6 @@
         dzoomx   = tab_cntrl(25)
         dzoomy   = tab_cntrl(26)
+        taux     = tab_cntrl(28)
+        tauy     = tab_cntrl(29)
       ELSE
         fxyhypb = . FALSE .
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynredem.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynredem.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/dynredem.F	(revision 232)
@@ -94,4 +94,7 @@
        tab_cntrl(25) = dzoomx
        tab_cntrl(26) = dzoomy
+       tab_cntrl(27) = 0.
+       tab_cntrl(28) = taux
+       tab_cntrl(29) = tauy
       ELSE
        tab_cntrl(24) = 0.
@@ -99,5 +102,7 @@
        tab_cntrl(26) = dzoomy
        tab_cntrl(27) = 0.
-       IF( ysinus )  tab_cntrl(27) = 1. 
+       tab_cntrl(28) = 0.
+       tab_cntrl(29) = 0.
+       IF( ysinus )  tab_cntrl(27) = 1.
       ENDIF
 c
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/etat0_netcdf.F	(revision 232)
@@ -274,4 +274,5 @@
       q3d(:,:,:,:) = 0.0
       qd(:,:,:) = 0.0
+      q3d(:,:,:,:) = 0.0
       WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
      .                           maxval(qsat(:,:,:))
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/fluxstokenc.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/fluxstokenc.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/fluxstokenc.F	(revision 232)
@@ -1,4 +1,4 @@
       SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
-     . time_step,itau, fluxid, fluxvid,fluxdid )
+     . time_step,itau )
 
        USE IOIPSL
@@ -30,10 +30,7 @@
 
       REAL pbarvst(iip1,jjp1,llm),zistdyn
+	real dtcum
 
-
-      INTEGER iadvtr
-      integer ndex1d(1)
-      integer ndex2d(ip1jmp1)
-      integer ndex3dv(ip1jm*llm),ndex3d(ip1jmp1*llm)
+      INTEGER iadvtr,ndex(1) 
       integer nscal
       real tst(1),ist(1),istp(1)
@@ -53,13 +50,21 @@
      .  time_step,istdyn* time_step,istdyn* time_step,
      . nqmx, fluxid,fluxvid,fluxdid) 
+	
+	ndex(1) = 0
+        call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex)
+        call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex)
+	
+	ndex(1) = 0
+        nscal = 1
+        tst(1) = time_step
+        call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex)
+        ist(1)=istdyn
+        call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex)
+        istp(1)= istphy
+        call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex)
+	
 	first = .false.
 
       endif
-
-
-      ndex1d = 0
-      ndex2d = 0
-      ndex3dv = 0
-      ndex3d = 0
 
 
@@ -93,5 +98,4 @@
 c   Test pour savoir si on advecte a ce pas de temps
       IF ( iadvtr.EQ.istdyn ) THEN
-
 c    normalisation
       DO l=1,llm
@@ -124,44 +128,23 @@
 
          iadvtr=0
-
-c     write(*,*)'histwrite phis'
-      call histwrite(fluxid, 'phis', 1, phis, iip1*jjp1, ndex2d)
-c     write(*,*)'histwrite aire'
-      call histwrite(fluxid, 'aire', 1, aire, iip1*jjp1, ndex2d)
+	Print*,'ITAU auqel on stoke les fluxmasses',itau
 	
-      nscal = 1
-      tst(1) = time_step
-c     write(*,*)'histwrite dtvr'
-      call histwrite(fluxdid, 'dtvr', 1, tst, nscal, ndex1d)
-      ist(1)=istdyn
-c     write(*,*)'histwrite istdyn'
-      call histwrite(fluxdid, 'istdyn', 1, ist, nscal, ndex1d)
-      istp(1)= istphy
-c     write(*,*)'histwrite istphy'
-      call histwrite(fluxdid, 'istphy', 1, istp, nscal, ndex1d)
-
-c       write(*,*)'histwrite masse'	
-	call histwrite(fluxid, 'masse', itau, masse,
-     .               iip1*jjp1*llm, ndex3d)
+	call histwrite(fluxid, 'masse', itau, massem,
+     .               iip1*jjp1*llm, ndex)
 	
-c       write(*,*)'histwrite pbaru'	
 	call histwrite(fluxid, 'pbaru', itau, pbarug,
-     .               iip1*jjp1*llm, ndex3d)
+     .               iip1*jjp1*llm, ndex)
 	
-c       write(*,*)'histwrite pbarv'	
-	call histwrite(fluxvid, 'pbarv', itau, pbarvst,
-     .               iim*jjp1*llm, ndex3dv)
+	call histwrite(fluxvid, 'pbarv', itau, pbarvg,
+     .               iip1*jjm*llm, ndex)
 	
-c       write(*,*)'histwrite w'	
         call histwrite(fluxid, 'w' ,itau, wg, 
-     .             iip1*jjp1*llm, ndex3d) 
+     .             iip1*jjp1*llm, ndex) 
 	
-c       write(*,*)'histwrite teta'	
 	call histwrite(fluxid, 'teta' ,itau, tetac, 
-     .             iip1*jjp1*llm, ndex3d) 
+     .             iip1*jjp1*llm, ndex) 
 	
-c       write(*,*)'histwrite phi'	
 	call histwrite(fluxid, 'phi' ,itau, phic, 
-     .             iip1*jjp1*llm, ndex3d) 
+     .             iip1*jjp1*llm, ndex) 
 	
 C
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F	(revision 232)
@@ -1,4 +1,5 @@
        SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
-     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025)
+     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
+     , champmin,champmax                                               )
 
 c      Auteur :  P. Le Van 
@@ -11,9 +12,12 @@
 c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
 c     dzoom  etant  la distance totale de la zone du zoom
-c     tau  la transition ,   normalement  = 1  .
-c
+c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
+c
+c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
+c   ********************************************************************
+
 
        INTEGER nmax, nmax2
-       PARAMETER (  nmax = 50000, nmax2 = 2*nmax )
+       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
 
 #include "dimensions.h"
@@ -23,63 +27,112 @@
 c
        REAL xzoomdeg,dzoom,tau,grossism
+
+c    ......   arguments  de  sortie  ......
+
        REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
      ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
 
-c    ......   arguments  de  sortie  ......
-c
-       REAL xlon(iip1),xprimm(iip1),xuv
-       REAL xtild(0:nmax2)
-       REAL fhyp(0:nmax),ffdx(0:nmax),beta,Xprimt(0:nmax2)
-       REAL Xf(0:nmax2),xxpr(0:nmax2)
-       REAL xvrai(iip1),xxprim(iip1) 
-       REAL pi,depi,epsilon,xzoom
+c     .... variables locales  ....
+c
+       REAL*8 xlon(iip1),xprimm(iip1),xuv
+       REAL*8 xtild(0:nmax2)
+       REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
+       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
+       REAL*8 xvrai(iip1),xxprim(iip1) 
+       REAL*8 pi,depi,epsilon,xzoom,fa,fb
+       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
        INTEGER i,it,ik,iter,ii,idif
-       REAL xi,xo1,xint,xmoy,xlon2,fxm,Xprimin
-       REAL champmin,champmax
+       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
+       REAL*8 champmin,champmax
        INTEGER is2
        SAVE is2
-       REAL dlon1(iip1),dlon2(iip1),dlon3(iip1)
+
+       REAL*8 heavyside
+       EXTERNAL coefpoly,heavyside
 
        pi       = 2. * ASIN(1.)
        depi     = 2. * pi
-       epsilon  = 1.e-6
+       epsilon  = 1.e-3
        xzoom    = xzoomdeg * pi/180. 
 
-
+       IF( dzoom.LT.1.)  THEN
+         dzoom = dzoom * depi
+       ELSEIF( dzoom.LT. 25. ) THEN
+         WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
+     ,menter et relancer ! '
+         STOP 1
+       ELSE
+         dzoom = dzoom * pi/180.
+       ENDIF
+
+       WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
+       WRITE(6,24) xzoom,grossism,tau,dzoom
 
        DO i = 0, nmax2 
-         xtild(i) = FLOAT(i) /nmax2
-        IF( xtild(i).EQ. 0.5 )  xtild(i) = xtild(i) + 1.e-6
-       ENDDO
-
-       DO i = 1, nmax
-        fhyp(i) = TANH ( ( xtild(i) - 0.5*(1.- dzoom) )          /
-     ,                 ( tau * xtild(i) * ( 0.5 -xtild(i))) )
-       ENDDO
-
-        fhyp(   0  ) = - 1.
-        fhyp( nmax ) =   1.
+        xtild(i) = - pi + FLOAT(i) * depi /nmax2
+       ENDDO
+
+       DO i = 0, nmax2
+
+       fa  = tau*  ( dzoom/2.  - xtild(i) )
+       fb  = xtild(i) *  ( pi - xtild(i) )
+
+
+       IF( 200.* fb .LT. - fa )   THEN
+         fhyp ( i) = - 1.
+       ELSEIF( 200. * fb .LT. fa ) THEN
+         fhyp ( i) =   1.
+       ELSE
+         fhyp(i) =  TANH ( fa/fb )
+       ENDIF
+
+       IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
+       IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
+
+       ENDDO
 
 cc  ....  Calcul  de  beta  ....
-c
-       ffdx( 0 ) = 0.
-
-       DO i = 1, nmax
-        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
-        fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
-     ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
-        ffdx(i) = ffdx(i-1) + fxm * ( xtild(i) - xtild(i-1) )
-       ENDDO
-
-        beta  = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 )
+c   ............................
+
+       ffdx = 0.
+
+       DO i = nmax +1,nmax2
+
+       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
+       fa  = tau*  ( dzoom/2.  - xmoy )
+       fb  = xmoy *  ( pi - xmoy )
+
+       IF( 200.* fb .LT. - fa )   THEN
+         fxm = - 1.
+       ELSEIF( 200. * fb .LT. fa ) THEN
+         fxm =   1.
+       ELSE
+         fxm =  TANH ( fa/fb )
+       ENDIF
+
+       IF ( xmoy.EQ. 0. )  fhyp(i) =  1.
+       IF ( xmoy.EQ. pi )  fhyp(i) = -1.
+       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
+
+       ENDDO
+
+        beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
+
+       IF( 2.*beta - grossism.LE. 0.)  THEN
+        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
+     ,tine fxhyp est mauvaise ! '
+        WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
+     , ' et relancer ! ***  '
+        CALL ABORT
+       ENDIF
 c
 c   .....  calcul  de  Xprimt   .....
 c
        
-       DO i = 0, nmax
+       DO i = nmax, nmax2
         Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
        ENDDO
 c   
-       DO i = 0, nmax
+       DO i =  nmax+1, nmax2
         Xprimt( nmax2 - i ) = Xprimt( i )
        ENDDO
@@ -88,14 +141,28 @@
 c   .....  Calcul  de  Xf     ........
 
-        Xf(0) = 0.
-       DO i = 1, nmax
-        xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
-        fxm     = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) )    /
-     ,                 ( tau * xmoy * ( 0.5 -xmoy)) )
-        xxpr(i)    = beta + ( grossism - beta ) * fxm
-       ENDDO
-
-       DO i = 1,nmax
-         xxpr(nmax2-i+1) = xxpr(i)
+       Xf(0) = - pi
+
+       DO i =  nmax +1, nmax2
+
+       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
+       fa  = tau*  ( dzoom/2.  - xmoy )
+       fb  = xmoy *  ( pi - xmoy )
+
+       IF( 200.* fb .LT. - fa )   THEN
+         fxm = - 1.
+       ELSEIF( 200. * fb .LT. fa ) THEN
+         fxm =   1.
+       ELSE
+         fxm =  TANH ( fa/fb )
+       ENDIF
+
+       IF ( xmoy.EQ. 0. )  fxm =  1.
+       IF ( xmoy.EQ. pi )  fxm = -1.
+       xxpr(i)    = beta + ( grossism - beta ) * fxm
+
+       ENDDO
+
+       DO i = nmax+1, nmax2
+        xxpr(nmax2-i+1) = xxpr(i)
        ENDDO
 
@@ -103,10 +170,4 @@
          Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
         ENDDO
-        do i=1,nmax2
-        xf(i)=xf(i)/xf(nmax2)
-        enddo
-
-
-       PRINT *,' XF ',xf(0),xf(nmax),xf(nmax2)
 
 
@@ -117,4 +178,5 @@
 c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
 c
+      WRITE(6,18)
 c
       DO 5000  ik = 1, 4
@@ -130,66 +192,68 @@
        ENDIF
 
+      xo1   = 0.
 
       DO 1500 i = 1, iim
 
-      xlon2 = (  FLOAT(i) + xuv - 0.75) / FLOAT(iim) 
-
-      xo1   = 0.
-      xi    = xlon2
-c
-      DO 500 iter = 1,300
-
+      xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim) 
+
+      Xfi    = xlon2
+c
       DO 250 it =  nmax2,0,-1
-      IF( xi.GE.xtild(it))  GO TO 350
+      IF( Xfi.GE.Xf(it))  GO TO 350
 250   CONTINUE
 
       it = 0
-      xi = xtild(it)
 
 350   CONTINUE
-       IF(it.EQ.nmax2)  THEN
-        it       = nmax2 -1
-        xf(it+1) = 1.
-       ENDIF
-c  .................................................................
-c  ....  Interpolation entre  xi(it) et xi(it+1)   pour avoir X(xi)  
-c      .....           et   X'(xi)                             .....
-c  .................................................................
-
-       xint   = ( Xf(it+1)-Xf(it) ) / ( xtild(it+1)-xtild(it) )      * 
-     +                    ( xi-xtild(it) )  +  Xf(it)
-      Xprimin = ( Xprimt(it+1)-Xprimt(it) )/ ( xtild(it+1)-xtild(it) ) *
-     +                    ( xi-xtild(it) )  +  Xprimt(it)
-
-       xi = xi - (xint-xlon2)/Xprimin
-
-      IF( ABS(xi-xo1).LE.epsilon) GO TO 550
-      xo1 = xi
-c
+
+c    ......  Calcul de   Xf(xi)    ...... 
+c
+      xi  = xtild(it)
+
+      IF(it.EQ.nmax2)  THEN
+       it       = nmax2 -1
+       Xf(it+1) = pi
+      ENDIF
+c  .....................................................................
+c
+c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
+c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
+c          et (Xf(it+1),xtild(it+1) )
+
+       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
+     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
+
+       Xf1     = Xf(it)
+       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
+
+       DO 500 iter = 1,300
+        xi = xi - ( Xf1 - Xfi )/ Xprimin
+
+        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
+         xo1      = xi
+         xi2      = xi * xi
+         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
+         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
 500   CONTINUE
-      PRINT *,' ***   PAS DE SOLUTION  ****  ',i,xlon2
-        STOP 4
+        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
+          STOP 6
 550   CONTINUE
 
-       xxprim(i)   = depi/( FLOAT(iim) * Xprimin)
-       xvrai(i)    = depi *  (xi - 0.5) + xzoom
-
+       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
+       xvrai(i)  =  xi + xzoom
 
 1500   CONTINUE
 
        DO i = 1 , iim
-        xlon  (i)   = xvrai(i)
+        xlon(i)     = xvrai(i)
         xprimm(i)   = xxprim(i)
-cc        xxlon(i) = xlon(i)*180./pi
        ENDDO
         
-cc      PRINT *,' XLON avant '
-cc      PRINT 68,(xxlon(i),i=1,iim)
-       
        DO i = 1, iim -1
         IF( xvrai(i+1). LT. xvrai(i) )  THEN
-         PRINT *,' PBS.  avec  rlonu(',i+1,' plus petit que rlonu(',i,
-     ,   ')'
-         STOP
+         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
+     ,  ')'
+        STOP 7
         ENDIF
        ENDDO
@@ -205,12 +269,9 @@
        ENDDO
 
-       PRINT *,' LONGITUDES  min max ',champmin,champmax
-
       IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
                 GO TO 1600
       ELSE
-        PRINT 18
-        PRINT *,'Reorganisation des longitudes pour avoir entre - pi ',
-     ,  ' et  pi '
+       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
+     ,  ' et pi '
 c
         IF( xzoom.LE.0.)  THEN
@@ -219,6 +280,6 @@
            IF( xvrai(i).GE. - pi )  GO TO 80
           ENDDO
-            PRINT *, ' PBS. 1 '
-            STOP
+            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
+            STOP 8
  80       CONTINUE
           is2 = i
@@ -240,10 +301,9 @@
             IF( xvrai(i).LE. pi )  GO TO 90
            ENDDO
-             PRINT *,' PBS.  2 '
-              STOP
+             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
+              STOP 9
  90        CONTINUE
             is2 = i
           ENDIF
-cc           PRINT *,' IS2 ',is2
            idif = iim -is2
            DO ii = 1, is2
@@ -270,84 +330,76 @@
          ENDDO
 
-
          IF( ik.EQ.1 )  THEN
-         PRINT *, ' XLON aux pts. V-0.25   apres ( en  deg. ) '
-         PRINT 18
-         PRINT 68,xvrai
-         PRINT *,' XPRIM '
-         PRINT 68, xprimm
+c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
+c          WRITE(6,18) 
+c          WRITE(6,68) xvrai
+ccc          WRITE(6,*) ' XPRIM k ',ik
+ccc          WRITE(6,566)  xprimm
+
            DO i = 1,iim + 1
              rlonm025(i) = xlon( i )
             xprimm025(i) = xprimm(i)
            ENDDO
+
          ELSE IF( ik.EQ.2 )  THEN
-         PRINT 18
-         PRINT *, ' XLON aux pts. V   apres ( en  deg. ) '
-         PRINT 68,xvrai
-         PRINT *,' XPRIM '
-         PRINT 68, xprimm
+c          WRITE(6,18) 
+c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
+c          WRITE(6,68) xvrai
+ccc          WRITE(6,*) ' XPRIM k ',ik
+ccc          WRITE(6,566)  xprimm
+
            DO i = 1,iim + 1
              rlonv(i) = xlon( i )
             xprimv(i) = xprimm(i)
            ENDDO
-         ELSE IF( ik.EQ.3 )  THEN
-         PRINT 18
-         PRINT *, ' XLON aux pts. U   apres ( en  deg. ) '
-         PRINT 68,xvrai
-         PRINT *,' XPRIM '
-         PRINT 68, xprimm
+
+         ELSE IF( ik.EQ.3)   THEN
+c          WRITE(6,18) 
+c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
+c          WRITE(6,68) xvrai
+ccc          WRITE(6,*) ' XPRIM ik ',ik
+ccc          WRITE(6,566)  xprimm
+
            DO i = 1,iim + 1
              rlonu(i) = xlon( i )
             xprimu(i) = xprimm(i)
            ENDDO
+
          ELSE IF( ik.EQ.4 )  THEN
-         PRINT 18
-         PRINT *, ' XLON aux pts. V+0.25   apres ( en  deg. ) '
-         PRINT 68,xvrai
-         PRINT *,' XPRIM '
-         PRINT 68, xprimm
+c          WRITE(6,18) 
+c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
+c          WRITE(6,68) xvrai
+ccc          WRITE(6,*) ' XPRIM ik ',ik
+ccc          WRITE(6,566)  xprimm
+
            DO i = 1,iim + 1
              rlonp025(i) = xlon( i )
             xprimp025(i) = xprimm(i)
            ENDDO
+
          ENDIF
 
 5000    CONTINUE
 c
-c       ...........  fin  de la boucle  do 5000      ............
-
-c
-        DO i = 1, iim + 1
-         dlon1(i) = rlonm025(i) - rlonv(i)
-         dlon2(i) = rlonm025(i) - rlonp025(i)
-         dlon3(i) = rlonm025(i) - rlonu(i)
+       WRITE(6,18)
+c
+c    ...........  fin  de la boucle  do 5000      ............
+
+        DO i = 1, iim
+         xlon(i) = rlonv(i+1) - rlonv(i)
         ENDDO
-
-        DO i = 1, iim + 1
-         rlonm025(i) = rlonm025(i) + dlon1(i)
+        champmin =  1.e12
+        champmax = -1.e12
+        DO i = 1, iim
+         champmin = MIN( champmin, xlon(i) )
+         champmax = MAX( champmax, xlon(i) )
         ENDDO
-
-        DO i = 1, iim + 1
-         rlonv(i)    = rlonm025(i) - dlon1(i)
-         rlonp025(i) = rlonm025(i) - dlon2(i)
-         rlonu(i)    = rlonm025(i) - dlon3(i)
-        ENDDO
-
-        DO i = 1, iim
-         xprimu   (i) = rlonu(i+1) - rlonu(i)
-         xprimv   (i) = rlonv(i+1) - rlonv(i)
-         xprimm025(i) = rlonm025(i+1) - rlonm025(i)
-         xprimp025(i) = rlonp025(i+1) - rlonp025(i)
-        ENDDO
-         xprimu   (iip1) = xprimu   (1)
-         xprimv   (iip1) = xprimv   (1)
-         xprimm025(iip1) = xprimm025(1)
-         xprimp025(iip1) = xprimp025(1)
-
-
-18       FORMAT(/)
-68       FORMAT(1x,7f9.2)
-
-
-         RETURN
-         END
+         champmin = champmin * 180./pi
+         champmax = champmax * 180./pi
+
+18     FORMAT(/)
+24     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
+68     FORMAT(1x,7f9.2)
+
+       RETURN
+       END
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxyhyper.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxyhyper.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxyhyper.F	(revision 232)
@@ -1,4 +1,7 @@
-       SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy,deltay  , 
-     ,                      xzoom, grossx, dzoomx,taux,
+c
+c $Header$
+c
+       SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy  ,   
+     ,                       xzoom, grossx, dzoomx,taux  ,
      , rlatu,yprimu,rlatv,yprimv,rlatu1,  yprimu1,  rlatu2,  yprimu2  , 
      , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025)
@@ -19,13 +22,9 @@
 c     a) le grossissement du zoom  :  grossy  ( en y ) et grossx ( en x )
 c     b) l' extension     du zoom  :  dzoomy  ( en y ) et dzoomx ( en x )
-c     c) la raideur       du zoom  :   tau   , ici  =  1.
+c     c) la raideur de la transition du zoom  :   taux et tauy   
 c
-c     N.B :   le produit ( grossissement * extension )  doit etre <  1.
-c    *******
-c     En plus , il y a un autre parametre , moins important mais quand
-c     meme utile , c'est  deltay , deplacement de la zone du zoom en
-c     latitude   :  Si  on deplace  de 10. degres vers le nord  ( deltay
-c      = 10.  dans inigeom ) .
-c
+c  N.B : Il vaut mieux avoir   :   grossx * dzoomx <  pi    ( radians )
+c ******
+c                  et              grossy * dzoomy <  pi/2  ( radians )
 c
 #include "dimensions.h"
@@ -40,5 +39,5 @@
        REAL rlonu(iip1),xprimu(iip1),rlonv(iip1),xprimv(iip1),
      , rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),xprimp025(iip1)
-       REAL deltay
+       REAL*8  dxmin, dxmax , dymin, dymax
 
 c   ....   var. locales   .....
@@ -47,30 +46,30 @@
 c
 
-       CALL fyhyp ( yzoom, grossy, dzoomy,tauy, deltay , 
-     ,  rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 )
+       CALL fyhyp ( yzoom, grossy, dzoomy,tauy  , 
+     ,  rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
+     ,  dymin,dymax                                               )
+
+       CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
+     , xprimv,rlonu,xprimu,rlonp025,xprimp025 , dxmin,dxmax         )
 
 
-       CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv,
-     , xprimv,rlonu,xprimu,rlonp025,xprimp025 )
-
-
-        DO i = 1, iim
+        DO i = 1, iip1
           IF(rlonp025(i).LT.rlonv(i))  THEN
-           PRINT *,' Attention !  rlonp025 < rlonv',i
+           WRITE(6,*) ' Attention !  rlonp025 < rlonv',i
             STOP
           ENDIF
 
           IF(rlonv(i).LT.rlonm025(i))  THEN 
-           PRINT *,' Attention !  rlonm025 > rlonv',i
+           WRITE(6,*) ' Attention !  rlonm025 > rlonv',i
             STOP
           ENDIF
 
           IF(rlonp025(i).GT.rlonu(i))  THEN
-           PRINT *,' Attention !  rlonp025 > rlonu',i
+           WRITE(6,*) ' Attention !  rlonp025 > rlonu',i
             STOP
           ENDIF
         ENDDO
 
-        PRINT *,'  *** TEST DE COHERENCE  OK    POUR   FX **** '
+        WRITE(6,*) '  *** TEST DE COHERENCE  OK    POUR   FX **** '
 
 c
@@ -78,30 +77,30 @@
 c
        IF(rlatu1(j).LE.rlatu2(j))   THEN
-         PRINT *,' Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
+         WRITE(6,*)'Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j
          STOP 13
        ENDIF
 c
        IF(rlatu2(j).LE.rlatu(j+1))  THEN
-        PRINT *,' Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
+        WRITE(6,*)'Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j
         STOP 14
        ENDIF
 c
        IF(rlatu(j).LE.rlatu1(j))    THEN
-        PRINT *,' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
+        WRITE(6,*)' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j
         STOP 15
        ENDIF
 c
        IF(rlatv(j).LE.rlatu2(j))    THEN
-        PRINT *,' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j
+        WRITE(6,*)' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j
         STOP 16
        ENDIF
 c
        IF(rlatv(j).ge.rlatu1(j))    THEN
-        PRINT *,' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j
+        WRITE(6,*)' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j
         STOP 17
        ENDIF
 c
        IF(rlatv(j).ge.rlatu(j))     THEN
-        PRINT *,'Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j
+        WRITE(6,*) ' Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j
         STOP 18
        ENDIF
@@ -109,7 +108,28 @@
        ENDDO
 c
-       PRINT *,'  *** TEST DE COHERENCE  OK    POUR   FY **** '
-
+       WRITE(6,*) '  *** TEST DE COHERENCE  OK    POUR   FY **** '
 c
+        WRITE(6,18)
+        WRITE(6,*) '  Latitudes  '
+        WRITE(6,*) ' *********** '
+        WRITE(6,18)
+        WRITE(6,3)  dymin, dymax
+        WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
+     ,ametres  grossism , tau , dzoom pour Y et repasser ! '
+c
+        WRITE(6,18)
+        WRITE(6,*) '  Longitudes  '
+        WRITE(6,*) ' ************ '
+        WRITE(6,18)
+        WRITE(6,3)  dxmin, dxmax
+        WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
+     ,ametres  grossism , tau , dzoom pour Y et repasser ! '
+        WRITE(6,18)
+c
+3      Format(1x, ' Au centre du zoom , la longueur de la maille est',
+     ,  ' d environ ',f8.2 ,' degres  ',
+     , ' alors que la maille en dehors de la zone du zoom est d environ
+     , ', f8.2,' degres ' )
+18      FORMAT(/)
 
        RETURN
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/fyhyp.F	(revision 232)
@@ -1,5 +1,10 @@
-       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau,deltay ,  
-     ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ) 
-
+c
+c $Header$
+c
+       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau  ,  
+     ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
+     ,  champmin,champmax                                            ) 
+
+cc    ...  Version du 01/04/2001 ....
 
        IMPLICIT NONE
@@ -13,14 +18,10 @@
 c
 c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc)
-c     dzoom  etant  la distance totale de la zone du zoom
-c     tau  la transition ,   normalement  = 1  .
-
-c  N.B :  on doit avoir :  grossism * dzoom  <   1 .
-c         **************
-c
-c    Pour Indoex , on a pris :
-c         *******
-c    grossism = 2.5 , dzoom = 7/24  en x et  y  , pour iim = 128 et jjm=64
-c    yzoomdeg = 0.  , tau = 1.  et delaty = 10.
+c     dzoom  etant  la distance totale de la zone du zoom ( en radians )
+c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom   
+c
+c
+c N.B : Il vaut mieux avoir : grossism * dzoom  <  pi/2  (radians) ,en lati.
+c      ********************************************************************
 c
 c
@@ -29,10 +30,11 @@
 
        INTEGER      nmax , nmax2
-       PARAMETER (  nmax = 50000, nmax2 = 2*nmax )
+       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
 c
 c
 c     .......  arguments  d'entree    .......
 c
-       REAL yzoomdeg, grossism,dzoom,tau , deltay
+       REAL yzoomdeg, grossism,dzoom,tau 
+c         ( rentres  par  run.def )
 
 c     .......  arguments  de sortie   .......
@@ -42,17 +44,18 @@
 
 c
-c    ..... Champs  locaux    .....
+c     .....     champs  locaux    .....
 c
      
-       REAl ylat(jjp1), yprim(jjp1)
-       REAL yuv
-       REAL ytild(0:nmax2)
-       REAL fhyp(0:nmax),ffdx(0:nmax),beta,Ytprim(0:nmax2)
-       SAVE Ytprim, ytild,Yf
-       REAL Yf(0:nmax2),yypr(0:nmax2)
-       REAL yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
-       REAL pi,depi,pis2,epsilon,yzoom
-       REAL yo1,yi,ylon2,fxm,ymoy,yint,Yprimin
-       REAL ypn,deply,y00
+       REAL*8 ylat(jjp1), yprim(jjp1)
+       REAL*8 yuv
+       REAL*8 yt(0:nmax2)
+       REAL*8 fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2)
+       SAVE Ytprim, yt,Yf
+       REAL*8 Yf(0:nmax2),yypr(0:nmax2)
+       REAL*8 yvrai(jjp1), yprimm(jjp1),ylatt(jjp1)
+       REAL*8 pi,depi,pis2,epsilon,y0,pisjm
+       REAL*8 yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax
+       REAL*8 yfi,Yf1,ffdy
+       REAL*8 ypn,deply,y00
        SAVE y00, deply
 
@@ -60,71 +63,119 @@
        INTEGER jpn,jjpn
        SAVE jpn
-
+       REAL*8 a0,a1,a2,a3,yi2,heavyy0,heavyy0m
+       REAL*8 fa(0:nmax2),fb(0:nmax2)
+       REAL y0min,y0max
+
+       REAL*8     heavyside
+       EXTERNAL   heavyside
 
        pi       = 2. * ASIN(1.)
        depi     = 2. * pi
        pis2     = pi/2.
-       epsilon  = 1.e-6
-       yzoom    = yzoomdeg * pi/180. 
-
-
+       pisjm    = pi/ FLOAT(jjm)
+       epsilon  = 1.e-3
+       y0       =  yzoomdeg * pi/180. 
+
+       IF( dzoom.LT.1.)  THEN
+         dzoom = dzoom * pi
+       ELSEIF( dzoom.LT. 12. ) THEN
+         WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug
+     ,menter et relancer ! '
+         STOP 1
+       ELSE
+         dzoom = dzoom * pi/180.
+       ENDIF
+
+       WRITE(6,18)
+       WRITE(6,*) ' yzoom( rad.),grossism,tau,dzoom (radians)'
+       WRITE(6,24) y0,grossism,tau,dzoom
 
        DO i = 0, nmax2 
-        ytild(i) = FLOAT(i) /nmax2
-       IF( ytild(i).EQ.0.5 )  ytild(i) = ytild(i) + 1.e-6
-       ENDDO
-
-       DO i = 1, nmax
-        fhyp(i) = TANH ( ( ytild(i) - 0.5*(1.- dzoom) )          /
-     ,                 ( tau * ytild(i) * ( 0.5 -ytild(i))) )
-       ENDDO
-
-       fhyp(   0  ) = - 1.
-       fhyp( nmax ) =   1.
+        yt(i) = - pis2  + FLOAT(i)* pi /nmax2
+       ENDDO
+
+       heavyy0m = heavyside( -y0 )
+       heavyy0  = heavyside(  y0 )
+       y0min    = 2.*y0*heavyy0m - pis2
+       y0max    = 2.*y0*heavyy0  + pis2
+
+       DO i = 0, nmax2 
+        IF( yt(i).LT.y0 )  THEN
+         fa (i) = tau*  (yt(i)-y0+dzoom/2. )
+         fb(i) =   (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) )
+        ELSEIF ( yt(i).GT.y0 )  THEN
+         fa(i) =   tau *(y0-yt(i)+dzoom/2. )
+         fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 
+       ENDIF
+        
+       IF( 200.* fb(i) .LT. - fa(i) )   THEN
+         fhyp ( i) = - 1.
+       ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
+         fhyp ( i) =   1.
+       ELSE  
+         fhyp(i) =  TANH ( fa(i)/fb(i) )
+       ENDIF
+
+       IF( yt(i).EQ.y0 )  fhyp(i) = 1.
+       IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1.
+
+       ENDDO
 
 cc  ....  Calcul  de  beta  ....
 c
-       ffdx( 0 ) = 0.
-
-       DO i = 1, nmax
-        ymoy    = 0.5 * ( ytild(i-1) + ytild( i ) )
-        fxm     = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) )    /
-     ,                 ( tau * ymoy * ( 0.5 -ymoy)) )
-        ffdx(i) = ffdx(i-1) + fxm * ( ytild(i) - ytild(i-1) )
-       ENDDO
-
-        beta  = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 )
+       ffdy   = 0.
+
+       DO i = 1, nmax2
+        ymoy    = 0.5 * ( yt(i-1) + yt( i ) )
+        IF( ymoy.LT.y0 )  THEN
+         fa(i)= tau * ( ymoy-y0+dzoom/2.) 
+         fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy )
+        ELSEIF ( ymoy.GT.y0 )  THEN
+         fa(i)= tau * ( y0-ymoy+dzoom/2. ) 
+         fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 )
+        ENDIF
+
+        IF( 200.* fb(i) .LT. - fa(i) )    THEN
+         fxm ( i) = - 1.
+        ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN
+         fxm ( i) =   1.
+        ELSE
+         fxm(i) =  TANH ( fa(i)/fb(i) )
+        ENDIF
+         IF( ymoy.EQ.y0 )  fxm(i) = 1.
+         IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1.
+         ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) )
+
+        ENDDO
+
+        beta  = ( grossism * ffdy - pi ) / ( ffdy - pi )
+
+       IF( 2.*beta - grossism.LE. 0.)  THEN
+
+        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
+     ,tine fyhyp est mauvaise ! '
+        WRITE(6,*)'Modifier les valeurs de  grossismy ,tauy ou dzoomy',
+     , ' et relancer ! ***  '
+        CALL ABORT
+
+       ENDIF
 c
 c   .....  calcul  de  Ytprim   .....
 c
        
-       DO i = 0, nmax
+       DO i = 0, nmax2
         Ytprim(i) = beta  + ( grossism - beta ) * fhyp(i)
        ENDDO
-c   
-       DO i = 0, nmax
-        Ytprim( nmax2 - i ) = Ytprim( i )
-       ENDDO
-c
 
 c   .....  Calcul  de  Yf     ........
 
-        Yf(0) = 0.
-       DO i = 1, nmax
-        ymoy    = 0.5 * ( ytild(i-1) + ytild( i ) )
-        fxm     = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) )    /
-     ,                 ( tau * ymoy * ( 0.5 -ymoy)) )
-        yypr(i)    = beta + ( grossism - beta ) * fxm
-       ENDDO
-
-       DO i = 1,nmax
-         yypr(nmax2-i+1) = yypr(i)
-       ENDDO
-
-        DO i=1,nmax2
-         Yf(i)   = Yf(i-1) + yypr(i) * ( ytild(i) - ytild(i-1) )
-        ENDDO
-c
-c
+       Yf(0) = - pis2
+       DO i = 1, nmax2
+        yypr(i)    = beta + ( grossism - beta ) * fxm(i)
+       ENDDO
+
+       DO i=1,nmax2
+        Yf(i)   = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) )
+       ENDDO
 
 c    ****************************************************************
@@ -133,4 +184,5 @@
 c   .....   yuv  = 0.5  si calcul des latitudes  aux pts.  V  .....
 c
+      WRITE(6,18)
 c
       DO 5000  ik = 1,4
@@ -149,28 +201,21 @@
          jlat = jjm 
        ENDIF
-        
-c
+c
+       yo1   = 0.
        DO 1500 j =  1,jlat
-
-        ylon2 =  ( FLOAT(j)  + yuv  -1.) / FLOAT(jjm) 
-
         yo1   = 0.
-        yi    = ylon2
-
-c
-       DO 500 iter = 1,300
-
+        ylon2 =  - pis2 + pisjm * ( FLOAT(j)  + yuv  -1.)  
+        yfi    = ylon2
+c
        DO 250 it =  nmax2,0,-1
-        IF( yi.GE.ytild(it))  GO TO 350
+        IF( yfi.GE.Yf(it))  GO TO 350
 250    CONTINUE
-
        it = 0
-       yi = ytild(it)
-
 350    CONTINUE
 
+       yi = yt(it)
        IF(it.EQ.nmax2)  THEN
         it       = nmax2 -1
-        Yf(it+1) = 1.
+        Yf(it+1) = pis2
        ENDIF
 c  .................................................................
@@ -179,41 +224,42 @@
 c  .................................................................
 
-       yint   = ( Yf(it+1)-Yf(it) ) / ( ytild(it+1)-ytild(it) )      * 
-     +                    ( yi-ytild(it) )  +  Yf(it)
-      Yprimin = ( Ytprim(it+1)-Ytprim(it) )/ ( ytild(it+1)-ytild(it) ) *
-     +                    ( yi-ytild(it) )  +  Ytprim(it)
-       yi = yi - (yint-ylon2)/Yprimin
-
-      IF( ABS(yi-yo1).LE.epsilon) GO TO 550
-      yo1 = yi
-c
+       CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1),   
+     ,                  yt(it),yt(it+1) ,   a0,a1,a2,a3   )      
+
+       Yf1     = Yf(it)
+       Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi
+
+       DO 500 iter = 1,300
+         yi = yi - ( Yf1 - yfi )/ Yprimin
+
+        IF( ABS(yi-yo1).LE.epsilon)  GO TO 550
+         yo1      = yi
+         yi2      = yi * yi
+         Yf1      = a0 +  a1 * yi +     a2 * yi2  +     a3 * yi2 * yi
+         Yprimin  =       a1      + 2.* a2 *  yi  + 3.* a3 * yi2
 500   CONTINUE
-      PRINT *,' ***   PAS DE SOLUTION  ****  ',j,ylon2,iter
-        STOP 4
+        WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter
+         STOP 2
 550   CONTINUE
-
-       yprim(j)    = pi /( FLOAT(jjm) *  Yprimin)
-       yvrai(j)    = pi *  (yi - 0.5) + yzoom
+c
+       Yprimin   = a1  + 2.* a2 *  yi   + 3.* a3 * yi* yi
+       yprim(j)  = pi / ( jjm * Yprimin )
+       yvrai(j)  = yi 
 
 1500    CONTINUE
-
-cc          print *,' LAT avant reorgan '
-cc         print 68,(yyvrai(j),j=1,jlat)
 
        DO j = 1, jlat -1
         IF( yvrai(j+1). LT. yvrai(j) )  THEN
-         PRINT *,' PBS.  avec  rlat(',j+1,' plus petit que rlat(',j,
-     ,   ')'
-         STOP
-        ENDIF
-       ENDDO
-
-        PRINT 18
-        PRINT *,'Reorganisation des latitudes pour avoir entre - pi/2 ',
-     ,  ' et  pi/2 '
-c
-           
+         WRITE(6,*) ' PBS. avec  rlat(',j+1,') plus petit que rlat(',j,
+     ,  ')'
+         STOP 3
+        ENDIF
+       ENDDO
+
+       WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2'
+     , ,' et  pi/2 '
+c
         IF( ik.EQ.1 )   THEN
-           ypn = pis2 - deltay * pi/180.
+           ypn = pis2 
           DO j = jlat,1,-1
            IF( yvrai(j).LE. ypn ) GO TO 1502
@@ -239,11 +285,7 @@
          ENDDO
 
-
-
 c      ***********   Fin de la reorganisation     *************
 c
  1600   CONTINUE
-
-
 
        DO j = 1, jlat
@@ -256,33 +298,35 @@
         ENDDO
 
-
         IF( ik.EQ.1 )  THEN
-        PRINT 18
-        PRINT *, ' YLAT  en U   apres ( en  deg. ) '
-        PRINT 68,(yvrai(j),j=1,jlat)
-        PRINT *,' YPRIM '
-        PRINT 68,( yprim(j),j=1,jlat)
+c         WRITE(6,18) 
+c         WRITE(6,*)  ' YLAT  en U   apres ( en  deg. ) '
+c         WRITE(6,68) (yvrai(j),j=1,jlat)
+cc         WRITE(6,*) ' YPRIM '
+cc         WRITE(6,445) ( yprim(j),j=1,jlat)
+
           DO j = 1, jlat
             rrlatu(j) =  ylat( j )
            yyprimu(j) = yprim( j )
           ENDDO
-c
+
         ELSE IF ( ik.EQ. 2 )  THEN
-        PRINT 18
-        PRINT *, ' YLAT   en V  apres ( en  deg. ) '
-        PRINT 68,(yvrai(j),j=1,jlat)
-        PRINT *,' YPRIM '
-        PRINT 68,( yprim(j),j=1,jlat)
+c         WRITE(6,18) 
+c         WRITE(6,*) ' YLAT   en V  apres ( en  deg. ) '
+c         WRITE(6,68) (yvrai(j),j=1,jlat)
+cc         WRITE(6,*)' YPRIM '
+cc         WRITE(6,445) ( yprim(j),j=1,jlat)
+
           DO j = 1, jlat
             rrlatv(j) =  ylat( j )
            yyprimv(j) = yprim( j )
           ENDDO
-c
+
         ELSE IF ( ik.EQ. 3 )  THEN
-        PRINT 18
-        PRINT *, ' YLAT  en U + 0.75  apres ( en  deg. ) '
-        PRINT 68,(yvrai(j),j=1,jlat)
-        PRINT *,' YPRIM '
-        PRINT 68,( yprim(j),j=1,jlat)
+c         WRITE(6,18) 
+c         WRITE(6,*)  ' YLAT  en U + 0.75  apres ( en  deg. ) '
+c         WRITE(6,68) (yvrai(j),j=1,jlat)
+cc         WRITE(6,*) ' YPRIM '
+cc         WRITE(6,445) ( yprim(j),j=1,jlat)
+
           DO j = 1, jlat
             rlatu2(j) =  ylat( j )
@@ -291,19 +335,36 @@
 
         ELSE IF ( ik.EQ. 4 )  THEN
-        PRINT 18
-        PRINT *, ' YLAT en U + 0.25  apres ( en  deg. ) '
-        PRINT 68,(yvrai(j),j=1,jlat)
-        PRINT *,' YPRIM '
-        PRINT 68,( yprim(j),j=1,jlat)
+c         WRITE(6,18) 
+c         WRITE(6,*)  ' YLAT en U + 0.25  apres ( en  deg. ) '
+c         WRITE(6,68)(yvrai(j),j=1,jlat)
+cc         WRITE(6,*) ' YPRIM '
+cc         WRITE(6,68) ( yprim(j),j=1,jlat)
+
           DO j = 1, jlat
             rlatu1(j) =  ylat( j )
            yprimu1(j) = yprim( j )
           ENDDO
+
         ENDIF
 
 5000   CONTINUE
 c
+        WRITE(6,18)
+c
 c  .....     fin de la boucle  do 5000 .....
 
+        DO j = 1, jjm
+         ylat(j) = rrlatu(j) - rrlatu(j+1)
+        ENDDO
+        champmin =  1.e12
+        champmax = -1.e12
+        DO j = 1, jjm
+         champmin = MIN( champmin, ylat(j) )
+         champmax = MAX( champmax, ylat(j) )
+        ENDDO
+         champmin = champmin * 180./pi
+         champmax = champmax * 180./pi
+
+24     FORMAT(2x,'Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
 18      FORMAT(/)
 68      FORMAT(1x,7f9.2)
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F	(revision 232)
@@ -22,11 +22,13 @@
 c
 c=======================================================================
-
-c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv 
-c        et possibilite d'appeler une fonction f(y)  a derivee tangente 
-c        hyperbolique a la  place de la fonction a derivee sinusoidale.         
-
-c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
-c         q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
+c
+c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
+c      et possibilite d'appeler une fonction f(y)  a derivee tangente
+c      hyperbolique a la  place de la fonction a derivee sinusoidale.
+
+c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
+c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
+c
+c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
 c
 c-----------------------------------------------------------------------
@@ -135,4 +137,7 @@
       LOGICAL true_calendar
       PARAMETER (true_calendar = .false.)
+C Run nudge
+      LOGICAL ok_nudge
+      PARAMETER (ok_nudge = .false.)
 
 c-----------------------------------------------------------------------
@@ -160,4 +165,6 @@
 c     iadv = 2    shema   amont
 c     iadv = 3    shema  Van-leer
+c     iadv = 4    schema  Van-leer + humidite specifique
+c                        Modif F.Codron
 c 
 c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
@@ -171,5 +178,20 @@
       ENDDO
 c
-      DO  iq = 1, nqmx
+       IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
+     * ,' pour la vapeur d''eau'
+       IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
+     * ,' vapeur d''eau '
+       IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
+     * ,'la vapeur d''eau'
+       IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + humidite'
+     * ,' specifique pour la vapeur d''eau'
+c
+       IF( iadv(1).LE.0.OR.iadv(1).GT.4 )   THEN
+        PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
+     * ,  car  iadv(1) = ', iadv(1)
+         CALL ABORT
+       ENDIF
+
+      DO  iq = 2, nqmx
        IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
      * ,' pour le traceur no ', iq
@@ -178,8 +200,16 @@
        IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
      * ,'le traceur no ', iq
-       IF( iadv(iq).LE.0.OR.iadv(iq).GT.3 )   THEN
+
+       IF( iadv(iq).EQ.4 )  THEN
+         PRINT *,' Le shema  Van-Leer + humidite specifique ',
+     * ' est  uniquement pour la vapeur d eau .'
+         PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
+         CALL ABORT
+       ENDIF
+
+       IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
        PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
      * . '
-         STOP 
+         CALL ABORT 
        ENDIF
       ENDDO
@@ -188,5 +218,5 @@
       numvanle = nqmx + 1
       DO  iq = 1, nqmx
-        IF( iadv(iq).EQ.3.AND.first ) THEN
+        IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
           numvanle = iq
           first    = .FALSE.
@@ -195,10 +225,10 @@
 c
       DO  iq = 1, nqmx
-        IF( iadv(iq).NE.3.AND.iq.GT.numvanle )   THEN
+      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
           PRINT *,' Il y a discontinuite dans le choix du shema de ',
      *    'Van-leer pour les traceurs . Corriger et repasser . '
-           STOP
-        ENDIF
-        IF( iadv(iq).LT.1.OR.iadv(iq).GT.3 )    THEN
+           CALL ABORT
+      ENDIF
+        IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 )    THEN
           PRINT *,' Le choix de  iadv  est errone pour le traceur ',
      *    iq
@@ -297,5 +327,8 @@
 
    1  CONTINUE
-c     call nudge(itau,ucov,vcov,teta,masse,ps)
+
+      if (ok_nudge) then
+        call nudge(itau,ucov,vcov,teta,masse,ps)
+      endif
 c
 c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
@@ -379,5 +412,7 @@
 c
             CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
-     *                             p,  masse, dq                  )
+     *                      p, masse, dq,  iadv(1), teta, pk     )
+c
+c                   ...  Modif  F.Codron  ....
 c
          ENDIF
@@ -389,5 +424,6 @@
 
             CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
-     .   time_step, itau,fluxid, fluxvid,fluxdid )
+     .   time_step, itau)
+
 
          ENDIF
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/inigeom.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/inigeom.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/inigeom.F	(revision 232)
@@ -1,22 +1,18 @@
+c
+c $Header$
+c
       SUBROUTINE inigeom
 c
 c     Auteur :  P. Le Van
-c    .....................
-c
-c   ............      Version  du 20/12/98     ........................
+c
+c   ............      Version  du 01/04/2001     ........................
 c
 c  Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4  aux memes en-
 c     endroits que les aires aireij1,..aireij4 .
+
 c  Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol.
 c
 c
       IMPLICIT NONE
-c
-      REAL        deltay, tauy, taux
-      PARAMETER ( deltay =  0. , tauy = 1. , taux = 1. )
-c
-c     deltay  est ( en degres ) le deplacement eventuel en Y du zoom
-cc    taux  et  tauy  sont  les raideurs   du  zoom
-c
 c
 #include "dimensions.h"
@@ -50,10 +46,6 @@
       SAVE rlonm025,xprimm025,rlonp025,xprimp025
 
-
-
-c----------------------------------------------------------------------
       REAL      SSUM
       EXTERNAL  SSUM
-c
 c
 c
@@ -62,5 +54,4 @@
 c   -    calcul des coeff. ( cu, cv , 1./cu**2,  1./cv**2  )         -
 c   -                                                                -
-c   ------------------------------------------------------------------
 c   ------------------------------------------------------------------
 c
@@ -168,5 +159,5 @@
 c
 c
-      PRINT 3
+      WRITE(6,3) 
  3    FORMAT( // 10x,' ....  INIGEOM  date du 01/06/98   ..... ',
      * //5x,'   Calcul des elongations cu et cv  comme sommes des 4 ' /
@@ -191,10 +182,10 @@
       ENDIF
 
-      PRINT *,' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
+      WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,
      *  nitergdiv,nitergrot,niterh
 c
       pi    = 2.* ASIN(1.)
 c
-      PRINT 990
+      WRITE(6,990) 
 
 c     ----------------------------------------------------------------
@@ -205,5 +196,5 @@
        IF( ysinus )  THEN
 c
-        PRINT *,' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
+        WRITE(6,*) ' ***  Inigeom ,  Y = Sinus ( Latitude ) *** '
 c
 c   .... utilisation de f(x,y )  avec  y  =  sinus de la latitude  .....
@@ -215,5 +206,5 @@
        ELSE
 c
-        PRINT *,' *** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
+        WRITE(6,*) '*** Inigeom ,  Y = Latitude  , der. sinusoid . ***'
 
 c  .... utilisation  de f(x,y) a tangente sinusoidale , y etant la latit. ...
@@ -267,11 +258,11 @@
       ELSE
 c
-c   ....  utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
+c   ....  Utilisation  de fxyhyper , f(x,y) a derivee tangente hyperbol.
 c   .....................................................................
 
-       PRINT *,'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
-
-       CALL fxyhyper( clat, grossismy, dzoomy, tauy, deltay  ,
-     ,               clon, grossismx, dzoomx, taux          ,
+      WRITE(6,*)'*** Inigeom , Y = Latitude  , der.tg. hyperbolique ***'
+ 
+       CALL fxyhyper( clat, grossismy, dzoomy, tauy    , 
+     ,                clon, grossismx, dzoomx, taux    ,
      , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2  ,
      , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 )
@@ -387,9 +378,8 @@
 c
       IF ( j. eq. jjp1 )  THEN
-
        yprp               = yprimu2(j-1)
        rlatp              = rlatu2 (j-1)
-cc       yprp             = fyprim( FLOAT(j) - 0.25 )
-cc       rlatp            = fy    ( FLOAT(j) - 0.25 )
+ccc       yprp             = fyprim( FLOAT(j) - 0.25 )
+ccc       rlatp            = fy    ( FLOAT(j) - 0.25 )
 c
       coslatp             = COS( rlatp )
@@ -425,8 +415,8 @@
         rlatm    = rlatu1 (  j  )
         yprm     = yprimu1(  j  )
-c         rlatp    = fy    ( FLOAT(j) - 0.25 )
-c         yprp     = fyprim( FLOAT(j) - 0.25 )
-c         rlatm    = fy    ( FLOAT(j) + 0.25 )
-c         yprm     = fyprim( FLOAT(j) + 0.25 )
+cc         rlatp    = fy    ( FLOAT(j) - 0.25 )
+cc         yprp     = fyprim( FLOAT(j) - 0.25 )
+cc         rlatm    = fy    ( FLOAT(j) + 0.25 )
+cc         yprm     = fyprim( FLOAT(j) + 0.25 )
 
          coslatm  = COS( rlatm )
@@ -490,5 +480,4 @@
   36  CONTINUE
 c
-c   ....  Modif  P. Le Van  ( 4/07/96 )  .....
 c
       aire    (iip1,j) = aire    (1,j)
@@ -620,5 +609,5 @@
          aiuscv2gam(iip1,j)  = aiuscv2gam(1,j)
       ENDDO
-c
+
 c
 c   calcul des aires aux  poles :
@@ -666,30 +655,43 @@
 c-----------------------------------------------------------------------
 c
-       PRINT *,' INIGEOM  RLONV '
+       WRITE(6,*) '   ***  Coordonnees de la grille  *** '
+       WRITE(6,995)
+c
+       WRITE(6,*) '   LONGITUDES  aux pts.   V  ( degres )  '
+       WRITE(6,995)
         DO i=1,iip1
          rlonvv(i) = rlonv(i)*180./pi
         ENDDO
-       PRINT 400,rlonvv
-c
-       PRINT *,' RLATV '
+       WRITE(6,400) rlonvv
+c
+       WRITE(6,995)
+       WRITE(6,*) '   LATITUDES   aux pts.   V  ( degres )  '
+       WRITE(6,995)
         DO i=1,jjm
          rlatuu(i)=rlatv(i)*180./pi
         ENDDO
-       PRINT 400,(rlatuu(i),i=1,jjm)
+       WRITE(6,400) (rlatuu(i),i=1,jjm)
 c
         DO i=1,iip1
           rlonvv(i)=rlonu(i)*180./pi
         ENDDO
-       PRINT *,' RLONU '
-       PRINT 400,rlonvv
-c
-       PRINT *,' RLATU '
+       WRITE(6,995)
+       WRITE(6,*) '   LONGITUDES  aux pts.   U  ( degres )  '
+       WRITE(6,995)
+       WRITE(6,400) rlonvv
+       WRITE(6,995)
+
+       WRITE(6,*) '   LATITUDES   aux pts.   U  ( degres )  '
+       WRITE(6,995)
         DO i=1,jjp1
          rlatuu(i)=rlatu(i)*180./pi
         ENDDO
-       PRINT 400,(rlatuu(i),i=1,jjp1)
-c
+       WRITE(6,400) (rlatuu(i),i=1,jjp1)
+       WRITE(6,995)
+c
+444    format(f10.3,f6.0)
 400    FORMAT(1x,8f8.2)
 990    FORMAT(//)
+995    FORMAT(/)
 c
       RETURN
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/paramet.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/paramet.h	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/paramet.h	(revision 232)
@@ -7,5 +7,6 @@
       INTEGER jcfil,jcfllm
 
-      PARAMETER( iip1= iim+1, iip2= iim+2, iip3= iim+3, jjp1= jjm+1 )
+      PARAMETER( iip1= iim+1-1/iim,iip2=iim+2,iip3=iim+3
+     s    ,jjp1=jjm+1-1/jjm)
       PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
       PARAMETER( kftd  = iim/2 -ndm )
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/physdem.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/physdem.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/physdem.F	(revision 232)
@@ -527,5 +527,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -540,5 +539,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -553,5 +551,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -566,5 +563,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -579,5 +575,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zpic)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -592,5 +587,4 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,zval)
 #endif
-c
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/serre.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/serre.h	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/serre.h	(revision 232)
@@ -1,6 +1,9 @@
+c
+c $Header$
+c
 c..include serre.h
 c
        REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo,
-     ,  grossismx, grossismy, dzoomx, dzoomy
+     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
        COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo , 
-     ,  grossismx, grossismy, dzoomx, dzoomy
+     ,  grossismx, grossismy, dzoomx, dzoomy,taux,tauy
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/tracvl.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/tracvl.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/tracvl.F	(revision 232)
@@ -1,4 +1,4 @@
       SUBROUTINE tracvl(numvanle,iapp_tracvl,nq,pbaru,pbarv ,
-     *                            p,  masse , q, iapptrac    )
+     *                    p, masse , q, iapptrac, iadv1, teta, pk  )
 c
 c     Auteur :  F. Hourdin
@@ -6,4 +6,6 @@
 c
 ccc   ..   Modif. P. Le Van  ( 20/12/97 )  ...
+c                 F. Codron     (10/99)
+
 c
       IMPLICIT NONE
@@ -15,10 +17,15 @@
 #include "comgeom.h"
 
-      INTEGER nq,iapp_tracvl
+c     .... Arguments  ....
+c
+      INTEGER numvanle, nq, iapp_tracvl, iapptrac, iadv1
 
       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
       REAL q(ip1jmp1,llm,nq),masse(ip1jmp1,llm)
-      REAL p( ip1jmp1,llmp1 )
+      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
+      REAL pk(ip1jmp1,llm)
 
+c     ....  var. locales  .....
+c
       REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
       REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
@@ -26,13 +33,13 @@
       REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
 
-      INTEGER iapptrac
 
-      INTEGER iadvtr, numvanle
+      INTEGER iadvtr, numvan
       INTEGER ij,l,iq
       REAL zdpmin, zdpmax
       EXTERNAL  minmax
- 
-      SAVE iadvtr, massem,pbaruc,pbarvc
+      SAVE iadvtr, massem, pbaruc, pbarvc, numvan
       DATA iadvtr/0/
+
+      numvan = numvanle
 
       IF(iadvtr.EQ.0) THEN
@@ -99,5 +106,14 @@
 
 c   Advection proprement dite.
-         DO iq = numvanle, nq
+c
+c   test sur iadv1 pour le schema de vapeur d'eau
+c
+         IF (numvanle.EQ.1.AND.iadv1.EQ.4) THEN
+           CALL vlspltqs( q(1,1,1), 2., massem, wg ,
+     *                 pbarug,pbarvg,dtvr,p,pk,teta )
+           numvan = 2
+         ENDIF
+
+         DO iq = numvan, nq
           CALL vlsplt( q(1,1,iq), 2. ,massem,wg,pbarug,pbarvg,dtvr )
          ENDDO
Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/vanleer.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/vanleer.F	(revision 231)
+++ LMDZ.3.3/branches/rel-LF/libf/dyn3d/vanleer.F	(revision 232)
@@ -1,8 +1,10 @@
       SUBROUTINE vanleer(numvanle,iapp_tracvl,nq,q,pbaru,pbarv ,
-     *                                    p ,masse,  dq          )
+     *                     p ,masse, dq ,  iadv1, teta, pk      )
 c
       IMPLICIT NONE
 c
-c     Auteurs:   F.Hourdin , P.Le Van, F.Forget 
+c     Auteurs:   F.Hourdin , P.Le Van, F.Forget, F.Codron  
+c
+c     F.Codron (10/99) : ajout humidite specifique pour eau vapeur
 c=======================================================================
 c
@@ -18,8 +20,8 @@
 c   Arguments:
 c   ----------
-      INTEGER nq,numvanle,iapp_tracvl 
+      INTEGER nq, numvanle, iapp_tracvl, iadv1 
       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
       REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nq),dq( ip1jmp1,llm,nq )
-
+      REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
 c  ..................................................................
 c
@@ -28,4 +30,6 @@
 c     numvanle est le numero du premier traceur qui appelle le
 c       shema de Van-Leer  (  1<=  numvanle <=  nqmx )  .
+c     si numvanle=1, iadv1=3 si Van-Leer seul, iadv1=4 si 
+c     Van-Leer+humidite specifique.
 c  ..................................................................
 c
@@ -61,5 +65,5 @@
 
       CALL tracvl( numvanle,iapp_tracvl,nq,pbaru,pbarv,p , masse,q  ,
-     *                                iapptrac                       )
+     *                      iapptrac, iadv1, teta ,pk              )
 
       IF( numvanle.EQ.1 ) THEN
