Index: LMDZ4/branches/pre_V3/libf/dyn3d/PVtheta.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/PVtheta.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/PVtheta.F	(revision 701)
@@ -0,0 +1,196 @@
+      SUBROUTINE PVtheta(ilon,ilev,pucov,pvcov,pteta,
+     $           ztfi,zplay,zplev,
+     $           nbteta,theta,PVteta)
+      IMPLICIT none
+
+c=======================================================================
+c
+c   Auteur:  I. Musat
+c   -------
+c
+c   Objet:
+c   ------
+c
+c    *******************************************************************
+c    Calcul de la vorticite potentielle PVteta sur des iso-theta selon
+c    la methodologie du NCEP/NCAR :
+c    1) on calcule la stabilite statique N**2=g/T*(dT/dz+g/cp) sur les
+c       niveaux du modele => N2
+c    2) on interpole les vents, la temperature et le N**2 sur des isentropes
+c       (en fait sur des iso-theta) lineairement en log(theta) =>
+c       ucovteta, vcovteta, N2teta
+c    3) on calcule la vorticite absolue sur des iso-theta => vorateta
+c    4) on calcule la densite rho sur des iso-theta => rhoteta 
+c
+c       rhoteta = (T/theta)**(cp/R)*p0/(R*T)
+c
+c    5) on calcule la vorticite potentielle sur des iso-theta => PVteta
+c
+c       PVteta = (vorateta * N2 * theta)/(g * rhoteta) ! en PVU
+c
+c       NB: 1PVU=10**(-6) K*m**2/(s * kg)
+c
+c       PVteta =  vorateta * N2/(g**2 * rhoteta) ! en 1/(Pa*s)
+c
+c
+c    *******************************************************************
+c
+c
+c     Variables d'entree : ilon,ilev,pucov,pvcov,pteta,ztfi,zplay,zplev,nbteta,theta
+c                       -> sur la grille dynamique
+c     Variable de sortie : PVteta
+c                       -> sur la grille physique 
+c=======================================================================
+
+#include "dimensions.h"
+#include "paramet.h"
+c
+c variables Input
+c
+      INTEGER ilon, ilev
+      REAL pvcov(iip1,jjm,ilev)
+      REAL pucov(iip1,jjp1,ilev)
+      REAL pteta(iip1,jjp1,ilev)
+      REAL ztfi(ilon,ilev)
+      REAL zplay(ilon,ilev), zplev(ilon,ilev+1)
+      INTEGER nbteta
+      REAL theta(nbteta)
+c
+c variable Output
+c
+      REAL PVteta(ilon,nbteta)
+c
+c variables locales
+c
+      INTEGER i, j, l, ig0
+      REAL SSUM
+      REAL teta(ilon, ilev)
+      REAL ptetau(ip1jmp1, ilev), ptetav(ip1jm, ilev)
+      REAL ucovteta(ip1jmp1,ilev), vcovteta(ip1jm,ilev)
+      REAL N2(ilon,ilev-1), N2teta(ilon,nbteta)
+      REAL ztfiteta(ilon,nbteta)
+      REAL rhoteta(ilon,nbteta)
+      REAL vorateta(iip1,jjm,nbteta)
+      REAL voratetafi(ilon,nbteta), vorpol(iim)
+c
+#include "comgeom2.h"
+#include "comconst.h"
+#include "comvert.h"
+c
+c projection teta sur la grille physique
+c
+      DO l=1,llm
+       teta(1,l)   =  pteta(1,1,l)
+       ig0         = 2
+       DO j = 2, jjm
+        DO i = 1, iim
+         teta(ig0,l)    = pteta(i,j,l)
+         ig0            = ig0 + 1
+        ENDDO
+       ENDDO
+       teta(ig0,l)    = pteta(1,jjp1,l)
+      ENDDO
+c
+c calcul pteta sur les grilles U et V
+c
+      DO l=1, llm
+       DO j=1, jjp1
+        DO i=1, iip1
+         ig0=i+(j-1)*iip1
+         ptetau(ig0,l)=pteta(i,j,l)
+        ENDDO !i
+       ENDDO !j
+       DO j=1, jjm
+        DO i=1, iip1
+         ig0=i+(j-1)*iip1
+         ptetav(ig0,l)=0.5*(pteta(i,j,l)+pteta(i,j+1,l))
+        ENDDO !i
+       ENDDO !j
+      ENDDO !l
+c
+c projection pucov, pvcov sur une surface de theta constante
+c
+      DO l=1, nbteta
+cIM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l),
+       CALL tetalevel(ip1jmp1,llm,.true.,ptetau,theta(l),
+     .                pucov,ucovteta(:,l))
+cIM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l),
+       CALL tetalevel(ip1jm,llm,.true.,ptetav,theta(l),
+     .                pvcov,vcovteta(:,l))
+      ENDDO !l
+c
+c calcul vorticite absolue sur une iso-theta : vorateta
+c
+      CALL tourabs(nbteta,vcovteta,ucovteta,vorateta)
+c
+c projection vorateta sur la grille physique => voratetafi
+c
+      DO l=1,nbteta
+       DO j=2,jjm
+        ig0=1+(j-2)*iim
+        DO i=1,iim
+         voratetafi(ig0+i+1,l) = vorateta( i ,j-1,l) * alpha4(i+1,j) +
+     $                           vorateta(i+1,j-1,l) * alpha1(i+1,j) +
+     $                           vorateta(i  ,j  ,l) * alpha3(i+1,j) +
+     $                           vorateta(i+1,j  ,l) * alpha2(i+1,j)
+        ENDDO
+        voratetafi(ig0 +1,l) = voratetafi(ig0 +1+ iim,l)
+       ENDDO
+      ENDDO
+c
+      DO l=1,nbteta
+       DO i=1,iim
+        vorpol(i)  = vorateta(i,1,l)*aire(i,1)
+       ENDDO
+       voratetafi(1,l)= SSUM(iim,vorpol,1)/apoln
+      ENDDO
+c
+      DO l=1,nbteta
+       DO i=1,iim
+        vorpol(i)  = vorateta(i,jjm,l)* aire(i,jjm +1)
+       ENDDO
+       voratetafi(ilon,l)= SSUM(iim,vorpol,1)/apols
+      ENDDO
+c 
+c calcul N**2 sur la grille physique => N2
+c
+      DO l=1, llm-1 
+       DO i=1, ilon
+        N2(i,l) = (g**2 * zplay(i,l) * 
+     $            (ztfi(i,l+1)-ztfi(i,l)) )/
+     $            (R*ztfi(i,l)*ztfi(i,l)*
+     $            (zplev(i,l)-zplev(i,l+1)) )+
+     $            (g**2)/(ztfi(i,l)*CPP)
+       ENDDO !i
+      ENDDO !l
+c
+c calcul N2 sur une iso-theta => N2teta 
+c
+      DO l=1, nbteta
+       CALL tetalevel(ilon,llm-1,.true.,teta,theta(l),
+     $                N2,N2teta(:,l))
+       CALL tetalevel(ilon,llm,.true.,teta,theta(l),
+     $                ztfi,ztfiteta(:,l))
+      ENDDO !l=1, nbteta
+c
+c calcul rho et PV sur une iso-theta : rhoteta, PVteta
+c
+      DO l=1, nbteta
+       DO i=1, ilon
+        rhoteta(i,l)=(ztfiteta(i,l)/theta(l))**(CPP/R)*
+     $  (preff/(R*ztfiteta(i,l)))
+c
+c PVteta en PVU
+c
+        PVteta(i,l)=(theta(l)*g*voratetafi(i,l)*N2teta(i,l))/
+     $              (g**2*rhoteta(i,l))
+c
+c PVteta en 1/(Pa*s)
+c
+        PVteta(i,l)=(voratetafi(i,l)*N2teta(i,l))/
+     $              (g**2*rhoteta(i,l))
+       ENDDO !i
+      ENDDO !l
+c
+      RETURN
+      END 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/advtrac.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/advtrac.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/advtrac.F	(revision 701)
@@ -4,5 +4,5 @@
 c
 c
-#ifdef INCA_CH4
+#ifdef INCA
       SUBROUTINE advtrac(pbaru,pbarv ,
      *                   p,  masse,q,iapptrac,teta,
@@ -16,5 +16,4 @@
      *                  pk)
 #endif
-
 c     Auteur :  F. Hourdin
 c
@@ -51,7 +50,7 @@
       REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
       REAL pk(ip1jmp1,llm)
-#ifdef INCA_CH4
-      INTEGER            :: hadv_flg(nq)
-      REAL               :: mmt_adj(ip1jmp1,llm)
+#ifdef INCA
+      INTEGER            :: hadv_flg(nqmx)
+      REAL               :: mmt_adj(ip1jmp1,llm,1)
       REAL               :: flxw(ip1jmp1,llm)
 #endif
@@ -129,5 +128,5 @@
         CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
 
-#ifdef INCA_CH4
+#ifdef INCA
       ! ... Flux de masse diaganostiques traceurs
       flxw = wg / FLOAT(iapp_tracvl)
@@ -214,5 +213,5 @@
         else if (iadv(iq).eq.20) then
             call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
-#ifdef INCA_CH4
+#ifdef INCA
        do iiq = iq+1, iq+3
          q(:,:,iiq)=q(:,:,iiq)*mmt_adj(:,:,1)
@@ -232,5 +231,5 @@
            call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
      s                     n,dtbon)
-#ifdef INCA_CH4
+#ifdef INCA
        do iiq = iq+1, iq+9
          q(:,:,iiq)=q(:,:,iiq)*mmt_adj(:,:,1)
@@ -338,6 +337,6 @@
 
 c                  call clock(t_final)
-                  tps_cpu=t_final-t_initial
-                  cpuadv(iq)=cpuadv(iq)+tps_cpu
+cym                  tps_cpu=t_final-t_initial
+cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
 
        end DO
Index: LMDZ4/branches/pre_V3/libf/dyn3d/advx.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/advx.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/advx.F	(revision 701)
@@ -106,5 +106,6 @@
          DO j = 1,jjp1
             DO i = 1,iim
-               sqi = sqi + S0(i,j,l,9)
+cIM 240305            sqi = sqi + S0(i,j,l,9)
+               sqi = sqi + S0(i,j,l,ntra)
             ENDDO
          ENDDO
@@ -483,5 +484,6 @@
         DO j = 1, jjp1
           DO i = 1, iim
-             sqf = sqf + S0(i,j,l,9)
+cIM 240405          sqf = sqf + S0(i,j,l,9)
+             sqf = sqf + S0(i,j,l,ntra)
           END DO  
         END DO
Index: LMDZ4/branches/pre_V3/libf/dyn3d/advz.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/advz.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/advz.F	(revision 701)
@@ -105,5 +105,6 @@
          DO j = 1,jjp1
             DO i = 1,iim
-               sqi = sqi + S0(i,j,l,9)
+cIM 240305            sqi = sqi + S0(i,j,l,9)
+               sqi = sqi + S0(i,j,l,ntra)
             ENDDO
          ENDDO
@@ -307,5 +308,6 @@
          DO j = 1,jjp1
             DO i = 1,iim
-               sqf = sqf + S0(i,j,l,9)
+cIM 240305            sqf = sqf + S0(i,j,l,9)
+               sqf = sqf + S0(i,j,l,ntra)
             ENDDO
          ENDDO
Index: LMDZ4/branches/pre_V3/libf/dyn3d/bilan_dyn.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/bilan_dyn.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/bilan_dyn.F	(revision 701)
@@ -2,7 +2,4 @@
 ! $Header$
 !
-c
-c $Header$
-c
       SUBROUTINE bilan_dyn (ntrac,dt_app,dt_cum,
      s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
@@ -63,6 +60,9 @@
 
 
-      character*6 nom(nQ)
-      character*6 unites(nQ)
+cym      character*6 nom(nQ)
+cym      character*6 unites(nQ)
+      character*6,save :: nom(nQ)
+      character*6,save :: unites(nQ)
+
       character*10 file
       integer ifile
@@ -116,7 +116,11 @@
       parameter (ntr=5)
 
-      character*10 znom(ntr,nQ)
-      character*20 znoml(ntr,nQ)
-      character*10 zunites(ntr,nQ)
+cym      character*10 znom(ntr,nQ)
+cym      character*20 znoml(ntr,nQ)
+cym      character*10 zunites(ntr,nQ)
+      character*10,save :: znom(ntr,nQ)
+      character*20,save :: znoml(ntr,nQ)
+      character*10,save :: zunites(ntr,nQ)
+
       integer iave,itot,immc,itrs,istn
       data iave,itot,immc,itrs,istn/1,2,3,4,5/
@@ -163,4 +167,6 @@
       time=time+dt_app
       itau=itau+1
+cIM
+      ndex3d=0
 
       if (first) then
@@ -301,4 +307,5 @@
 
 c   énergie cinétique
+      ucont(:,:,:)=0
       CALL covcont(llm,ucov,vcov,ucont,vcont)
       CALL enercin(vcov,ucov,vcont,ucont,ecin)
@@ -315,5 +322,5 @@
       Q(:,:,:,iang)=ang(:,:,:)
       Q(:,:,:,iu)=unat(:,:,:)
-      Q(:,:,:,iovap)=q(:,:,:,1)
+      Q(:,:,:,iovap)=trac(:,:,:,1)
       Q(:,:,:,iun)=1.
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/caladvtrac.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/caladvtrac.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/caladvtrac.F	(revision 701)
@@ -4,5 +4,5 @@
 c
 c
-#ifdef INCA_CH4
+#ifdef INCA
             SUBROUTINE caladvtrac(q,pbaru,pbarv ,
      *                   p ,masse, dq ,  teta,
@@ -41,7 +41,7 @@
       REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqmx),dq( ip1jmp1,llm,2 )
       REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
-#ifdef INCA_CH4
-      INTEGER            :: hadv_flg(nq)
-      REAL               :: mmt_adj(iip1,jjp1,llm)
+#ifdef INCA
+      INTEGER            :: hadv_flg(nqmx)
+      REAL               :: mmt_adj(iip1,jjp1,llm,1)
       REAL               :: flxw(ip1jmp1,llm)
 #endif
@@ -73,5 +73,5 @@
 c   advection
 
-#ifdef INCA_CH4
+#ifdef INCA
       CALL advtrac( pbaru,pbarv, 
      *             p,  masse,q,iapptrac, teta,
Index: LMDZ4/branches/pre_V3/libf/dyn3d/calfis.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/calfis.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/calfis.F	(revision 701)
@@ -23,5 +23,5 @@
      $                  pdq,
      $                  pw,
-#ifdef INCA_CH4
+#ifdef INCA
      $                  flxw,
 #endif
@@ -162,6 +162,13 @@
       REAL zsinbis(iim),zcosbis(iim),z1bis(iim)
       REAL unskap, pksurcp
-
-#ifdef INCA_CH4
+c
+cIM diagnostique PVteta, Amip2
+      INTEGER ntetaSTD
+      PARAMETER(ntetaSTD=3)
+      REAL rtetaSTD(ntetaSTD)
+      DATA rtetaSTD/350., 380., 405./
+      REAL PVteta(ngridmx,ntetaSTD)
+c
+#ifdef INCA
       REAL flxw(iip1,jjp1,llm)
       REAL flxwfi(ngridmx,llm)
@@ -430,7 +437,11 @@
 
       ENDDO
-
-
-#ifdef INCA_CH4
+c
+cIM calcul PV a teta=350, 380, 405K
+      CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
+     $           ztfi,zplay,zplev,
+     $           ntetaSTD,rtetaSTD,PVteta)
+c
+#ifdef INCA
       CALL gr_dyn_fi(llm,iip1,jjp1,ngridmx,flxw,flxwfi)
 #endif
@@ -461,5 +472,5 @@
      .             zqfi,
      .             pvervel,
-#ifdef INCA_CH4
+#ifdef INCA
      .             flxwfi,
 #endif
@@ -468,5 +479,8 @@
      .             zdtfi,
      .             zdqfi,
-     .             zdpsrf)
+     .             zdpsrf,
+cIM diagnostique PVteta, Amip2          
+     .             pducov,
+     .             PVteta)
 
 500   CONTINUE
Index: LMDZ4/branches/pre_V3/libf/dyn3d/conf_dat2d.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/conf_dat2d.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/conf_dat2d.F	(revision 701)
@@ -215,4 +215,7 @@
        ENDDO
 
+      deallocate(xtemp)
+      deallocate(ytemp)
+
       RETURN
       END
Index: LMDZ4/branches/pre_V3/libf/dyn3d/conf_gcm.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/conf_gcm.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/conf_gcm.F	(revision 701)
@@ -149,7 +149,7 @@
 !Config  Key  = iapp_tracvl
 !Config  Desc = frequence du groupement des flux 
-!Config  Def  = 5
+!Config  Def  = iperiod
 !Config  Help = frequence du groupement des flux (en pas de temps) 
-       iapp_tracvl = 5
+       iapp_tracvl = iperiod
        CALL getin('iapp_tracvl',iapp_tracvl)
 
@@ -365,4 +365,16 @@
        iflag_con = 2
        CALL getin('iflag_con',iflag_con)
+!
+!Config  Key  = ip_ebil_dyn
+!Config  Desc = PRINT level for energy conserv. diag.
+!Config  Def  = 0
+!Config  Help = PRINT level for energy conservation diag. ;
+!               les options suivantes existent :
+!Config         0 pas de print
+!Config         1 pas de print
+!Config         2 print,
+       ip_ebil_dyn = 0
+       CALL getin('ip_ebil_dyn',ip_ebil_dyn)
+!
 
       DO i = 1, longcles
@@ -574,4 +586,12 @@
       ENDIF
 c
+!Config  Key  = offline
+!Config  Desc = Nouvelle eau liquide
+!Config  Def  = n
+!Config  Help = Permet de mettre en route la
+!Config         nouvelle parametrisation de l'eau liquide !
+       offline = .FALSE.
+       CALL getin('offline',offline)
+
 
       write(lunout,*)' #########################################'
@@ -615,4 +635,5 @@
       write(lunout,*)' tauxx = ', tauxx
       write(lunout,*)' tauyy = ', tauyy
+      write(lunout,*)' offline = ', offline
 
       RETURN
@@ -720,4 +741,13 @@
        CALL getin('ysinus',ysinus)
 c
+!Config  Key  = offline
+!Config  Desc = Nouvelle eau liquide
+!Config  Def  = n
+!Config  Help = Permet de mettre en route la
+!Config         nouvelle parametrisation de l'eau liquide !
+       offline = .FALSE.
+       CALL getin('offline',offline)
+      write(lunout,*)' offline = ', offline
+
 
       write(lunout,*)' #########################################'
@@ -761,4 +791,5 @@
       write(lunout,*)' taux = ', taux
       write(lunout,*)' tauy = ', tauy
+      write(lunout,*)' offline = ', offline
 c
       RETURN
Index: LMDZ4/branches/pre_V3/libf/dyn3d/control.h
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/control.h	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/control.h	(revision 701)
@@ -6,11 +6,13 @@
 
       COMMON/control/nday,day_step,
-     $              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq ,
-     $              periodav,ecritphy,iecrimoy,dayref,anneeref,
-     $              raz_date
+     .              iperiod,iapp_tracvl,iconser,iecri,idissip,iphysiq ,
+     .              periodav,ecritphy,iecrimoy,dayref,anneeref,
+     .              raz_date,offline,ip_ebil_dyn
 
       INTEGER   nday,day_step,iperiod,iapp_tracvl,iconser,iecri,
-     $          idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date
+     .          idissip,iphysiq,iecrimoy,dayref,anneeref, raz_date
+     .          ,ip_ebil_dyn
       REAL periodav, ecritphy
+      logical offline
 
 c-----------------------------------------------------------------------
Index: LMDZ4/branches/pre_V3/libf/dyn3d/dynetat0.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/dynetat0.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/dynetat0.F	(revision 701)
@@ -106,8 +106,8 @@
       IF ( tab_cntrl(24).EQ.1. )  THEN
         fxyhypb  = . TRUE .
-        dzoomx   = tab_cntrl(25)
-        dzoomy   = tab_cntrl(26)
-        taux     = tab_cntrl(28)
-        tauy     = tab_cntrl(29)
+c        dzoomx   = tab_cntrl(25)
+c        dzoomy   = tab_cntrl(26)
+c        taux     = tab_cntrl(28)
+c        tauy     = tab_cntrl(29)
       ELSE
         fxyhypb = . FALSE .
Index: LMDZ4/branches/pre_V3/libf/dyn3d/dynredem.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/dynredem.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/dynredem.F	(revision 701)
@@ -144,5 +144,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"controle",NF_DOUBLE,1,idim_index,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"controle",NF_FLOAT,1,idim_index,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
      .                       "Parametres de controle")
@@ -155,5 +161,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlonu",NF_DOUBLE,1,idim_rlonu,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"rlonu",NF_FLOAT,1,idim_rlonu,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
      .                       "Longitudes des points U")
@@ -166,5 +178,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlatu",NF_DOUBLE,1,idim_rlatu,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"rlatu",NF_FLOAT,1,idim_rlatu,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
      .                       "Latitudes des points U")
@@ -177,5 +195,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlonv",NF_DOUBLE,1,idim_rlonv,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"rlonv",NF_FLOAT,1,idim_rlonv,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
      .                       "Longitudes des points V")
@@ -188,5 +212,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"rlatv",NF_DOUBLE,1,idim_rlatv,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"rlatv",NF_FLOAT,1,idim_rlatv,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
      .                       "Latitudes des points V")
@@ -199,5 +229,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"nivsigs",NF_DOUBLE,1,idim_s,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"nivsigs",NF_FLOAT,1,idim_s,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 28,
      .                       "Numero naturel des couches s")
@@ -210,5 +246,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"nivsig",NF_DOUBLE,1,idim_sig,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"nivsig",NF_FLOAT,1,idim_sig,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 32,
      .                       "Numero naturel des couches sigma")
@@ -221,5 +263,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"ap",NF_DOUBLE,1,idim_sig,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"ap",NF_FLOAT,1,idim_sig,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
      .                       "Coefficient A pour hybride")
@@ -232,5 +280,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"bp",NF_DOUBLE,1,idim_sig,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"bp",NF_FLOAT,1,idim_sig,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
      .                       "Coefficient B pour hybride")
@@ -243,5 +297,11 @@
 c
       ierr = NF_REDEF (nid)
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"presnivs",NF_DOUBLE,1,idim_s,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"presnivs",NF_FLOAT,1,idim_s,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
@@ -256,5 +316,11 @@
       dims2(1) = idim_rlonu
       dims2(2) = idim_rlatu
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"cu",NF_DOUBLE,2,dims2,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"cu",NF_FLOAT,2,dims2,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29,
      .                       "Coefficient de passage pour U")
@@ -269,5 +335,11 @@
       dims2(1) = idim_rlonv
       dims2(2) = idim_rlatv
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"cv",NF_DOUBLE,2,dims2,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"cv",NF_FLOAT,2,dims2,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 29,
      .                       "Coefficient de passage pour V")
@@ -284,5 +356,11 @@
       dims2(1) = idim_rlonv
       dims2(2) = idim_rlatu
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"aire",NF_DOUBLE,2,dims2,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"aire",NF_FLOAT,2,dims2,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 22,
      .                       "Aires de chaque maille")
@@ -299,5 +377,11 @@
       dims2(1) = idim_rlonv
       dims2(2) = idim_rlatu
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"phisinit",NF_DOUBLE,2,dims2,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"phisinit",NF_FLOAT,2,dims2,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
      .                       "Geopotentiel au sol")
@@ -313,5 +397,11 @@
       ierr = NF_REDEF (nid) ! entrer dans le mode de definition
 c
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"temps",NF_DOUBLE,1,idim_tim,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"temps",NF_FLOAT,1,idim_tim,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
      .                       "Temps de simulation")
@@ -326,5 +416,11 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"ucov",NF_DOUBLE,4,dims4,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"ucov",NF_FLOAT,4,dims4,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
      .                       "Vitesse U")
@@ -334,5 +430,11 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"vcov",NF_DOUBLE,4,dims4,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"vcov",NF_FLOAT,4,dims4,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
      .                       "Vitesse V")
@@ -342,5 +444,11 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"teta",NF_DOUBLE,4,dims4,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"teta",NF_FLOAT,4,dims4,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 11,
      .                       "Temperature")
@@ -352,5 +460,11 @@
       IF(nq.GE.1) THEN
       DO iq=1,nq
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,tname(iq),NF_DOUBLE,4,dims4,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,tname(iq),NF_FLOAT,4,dims4,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,ttext(iq))
       ENDDO
@@ -361,5 +475,11 @@
       dims4(3) = idim_s
       dims4(4) = idim_tim
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"masse",NF_DOUBLE,4,dims4,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"masse",NF_FLOAT,4,dims4,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,
      .                       "C est quoi ?")
@@ -368,5 +488,11 @@
       dims3(2) = idim_rlatu
       dims3(3) = idim_tim
+cIM 220306 BEG
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid,"ps",NF_DOUBLE,3,dims3,nvarid)
+#else
       ierr = NF_DEF_VAR (nid,"ps",NF_FLOAT,3,dims3,nvarid)
+#endif
+cIM 220306 END
       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 15,
      .                       "Pression au sol")
@@ -394,4 +520,5 @@
 #include "comgeom.h"
 #include "advtrac.h"
+#include "temps.h"
 
       INTEGER nq, l
@@ -406,4 +533,7 @@
       INTEGER ierr
       INTEGER iq
+      INTEGER length
+      PARAMETER (length = 100)
+      REAL tab_cntrl(length) ! tableau des parametres du run
       character*20 modname
       character*80 abort_message
@@ -435,4 +565,25 @@
 #endif
       PRINT*, "Enregistrement pour ", nb, time
+
+c
+c  Re-ecriture du tableau de controle, itaufin n'est plus defini quand
+c  on passe dans dynredem0
+      ierr = NF_INQ_VARID (nid, "controle", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         abort_message="dynredem1: Le champ <controle> est absent"
+         ierr = 1
+         CALL abort_gcm(modname,abort_message,ierr)
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
+#else
+      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
+#endif
+       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
 
 c  Ecriture des champs
Index: LMDZ4/branches/pre_V3/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/etat0_netcdf.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/etat0_netcdf.F	(revision 701)
@@ -43,5 +43,7 @@
       REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
       REAL :: solsw(klon), sollw(klon), fder(klon)
-      REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon,nbsrf)
+cIM "slab" ocean
+      REAL :: tslab(klon), seaice(klon)
+      REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
       REAL :: rugmer(klon)
       REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
@@ -51,4 +53,5 @@
       REAL :: pctsrf(klon, nbsrf)
       REAL :: t_ancien(klon,klev), q_ancien(klon,klev)      !
+      REAL :: run_off_lic_0(klon)
       real :: clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
       ! declarations pour lecture glace de mer
@@ -364,7 +367,13 @@
      .    jjm, rlonu, rlatv , interbar )
       !
-      varname = 'deltat'
-      deltat(:) = 0.0
-      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,deltat,0.0,
+cIM "slab" ocean
+      varname = 'tslab'
+      tslab(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,tslab,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+c
+      varname = 'seaice'
+      seaice(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,seaice,0.0,
      .     jjm, rlonu, rlatv , interbar )
       !
@@ -650,5 +659,8 @@
       q_ancien = 0.
       agesno = 0.
-      deltat = 0.
+cIM "slab" ocean
+      tslab(1:klon) = tsolsrf(1:klon,is_oce)
+      seaice = 0. !kg m-2
+c
       frugs(1:klon,is_oce) = rugmer(1:klon)
       frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
@@ -659,13 +671,16 @@
       rnebcon = 0.0
       ratqs = 0.0
+      run_off_lic_0 = 0.0
 
 cIM   call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
       call phyredem("startphy.nc",phystep,radpas,
-     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, qsol, 
+     $    latfi, lonfi, pctsrf, tsolsrf, tsoil, tslab, seaice, 
+     $    qsolsrf, qsol,
      $    snsrf, 
      $    albe, alblw, evap, rain_fall, snow_fall, solsw, sollw, fder,
      $    radsol, frugs,  agesno, 
      $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, 
-     $    t_ancien, q_ancien, rnebcon, ratqs, clwcon)
+     $    t_ancien, q_ancien, rnebcon, ratqs, clwcon,
+     $    run_off_lic_0)
       print*,'sortie phyredem'
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/fluxstokenc.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/fluxstokenc.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/fluxstokenc.F	(revision 701)
@@ -0,0 +1,162 @@
+      SUBROUTINE fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
+     . time_step,itau )
+
+       USE IOIPSL
+c
+c     Auteur :  F. Hourdin
+c
+c
+ccc   ..   Modif. P. Le Van  ( 20/12/97 )  ...
+c
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "tracstoke.h"
+#include "temps.h"
+
+      REAL time_step,t_wrt, t_ops
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL masse(ip1jmp1,llm),teta(ip1jmp1,llm),phi(ip1jmp1,llm)
+      REAL phis(ip1jmp1)
+
+      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+      REAL massem(ip1jmp1,llm),tetac(ip1jmp1,llm),phic(ip1jmp1,llm)
+
+      REAL pbarug(ip1jmp1,llm),pbarvg(iip1,jjm,llm),wg(ip1jmp1,llm)
+
+      REAL pbarvst(iip1,jjp1,llm),zistdyn
+	real dtcum
+
+      INTEGER iadvtr,ndex(1) 
+      integer nscal
+      real tst(1),ist(1),istp(1)
+      INTEGER ij,l,irec,i,j,itau
+      INTEGER, SAVE :: fluxid, fluxvid,fluxdid
+ 
+      SAVE iadvtr, massem,pbaruc,pbarvc,irec
+      SAVE phic,tetac
+      logical first
+      save first
+      data first/.true./
+      DATA iadvtr/0/
+
+
+c AC initialisations
+      pbarug(:,:)   = 0.
+      pbarvg(:,:,:) = 0.
+      wg(:,:)       = 0.
+      
+
+      if(first) then
+
+	CALL initfluxsto( 'fluxstoke',
+     .  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
+
+
+      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
+	Print*,'ITAU auqel on stoke les fluxmasses',itau
+	
+	call histwrite(fluxid, 'masse', itau, massem,
+     .               iip1*jjp1*llm, ndex)
+	
+	call histwrite(fluxid, 'pbaru', itau, pbarug,
+     .               iip1*jjp1*llm, ndex)
+	
+	call histwrite(fluxvid, 'pbarv', itau, pbarvg,
+     .               iip1*jjm*llm, ndex)
+	
+        call histwrite(fluxid, 'w' ,itau, wg, 
+     .             iip1*jjp1*llm, ndex) 
+	
+	call histwrite(fluxid, 'teta' ,itau, tetac, 
+     .             iip1*jjp1*llm, ndex) 
+	
+	call histwrite(fluxid, 'phi' ,itau, phic, 
+     .             iip1*jjp1*llm, ndex) 
+	
+C
+
+      ENDIF ! if iadvtr.EQ.istdyn
+
+      RETURN
+      END
Index: LMDZ4/branches/pre_V3/libf/dyn3d/fxhyp.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/fxhyp.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/fxhyp.F	(revision 701)
@@ -4,5 +4,5 @@
 c
 c
-       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
+       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau ,
      , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
      , champmin,champmax                                               )
@@ -38,5 +38,5 @@
 c     ......  arguments  d'entree   .......
 c
-       REAL xzoomdeg,dzoom,tau,grossism
+       REAL xzoomdeg,dzooma,tau,grossism
 
 c    ......   arguments  de  sortie  ......
@@ -47,4 +47,5 @@
 c     .... variables locales  ....
 c
+       REAL   dzoom
        REAL*8 xlon(iip1),xprimm(iip1),xuv
        REAL*8 xtild(0:nmax2)
@@ -74,12 +75,12 @@
        WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx
 c
-       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
+       IF( dzooma.LT.1.)  THEN
+         dzoom = dzooma * depi
+       ELSEIF( dzooma.LT. 25. ) THEN
+         WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
      ,menter et relancer ! '
          STOP 1
        ELSE
-         dzoom = dzoom * pi/180.
+         dzoom = dzooma * pi/180.
        ENDIF
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/fyhyp.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/fyhyp.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/fyhyp.F	(revision 701)
@@ -4,5 +4,5 @@
 c
 c
-       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau  ,  
+       SUBROUTINE fyhyp ( yzoomdeg, grossism, dzooma,tau  ,  
      ,  rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ,
      ,  champmin,champmax                                            ) 
@@ -37,5 +37,5 @@
 c     .......  arguments  d'entree    .......
 c
-       REAL yzoomdeg, grossism,dzoom,tau 
+       REAL yzoomdeg, grossism,dzooma,tau 
 c         ( rentres  par  run.def )
 
@@ -49,4 +49,5 @@
 c
      
+       REAL   dzoom
        REAL*8 ylat(jjp1), yprim(jjp1)
        REAL*8 yuv
@@ -78,12 +79,12 @@
        y0       =  yzoomdeg * pi/180. 
 
-       IF( dzoom.LT.1.)  THEN
-         dzoom = dzoom * pi
-       ELSEIF( dzoom.LT. 12. ) THEN
+       IF( dzooma.LT.1.)  THEN
+         dzoom = dzooma * pi
+       ELSEIF( dzooma.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.
+         dzoom = dzooma * pi/180.
        ENDIF
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/gcm.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/gcm.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/gcm.F	(revision 701)
@@ -56,6 +56,5 @@
 #include "com_io_dyn.h"
 #include "iniprint.h"
-
-c#include "tracstoke.h"
+#include "tracstoke.h"
 
 
@@ -113,12 +112,6 @@
       REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
       CHARACTER*15 ztit
-      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
-      SAVE      ip_ebil_dyn
-      DATA      ip_ebil_dyn/0/
 c-jld 
 
-
-      LOGICAL offline  ! Controle du stockage ds "fluxmass"
-      PARAMETER (offline=.false.)
 
       character*80 dynhist_file, dynhistave_file
@@ -368,4 +361,11 @@
 #endif
 
+c  Choix des frequences de stokage pour le offline
+c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
+c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
+      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
+      istphy=istdyn/iphysiq     
+
+
 c
 c-----------------------------------------------------------------------
@@ -380,5 +380,6 @@
 
 
-      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0)
+      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
+     .              time_0)
 
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/gradsdef.h
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/gradsdef.h	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/gradsdef.h	(revision 701)
@@ -3,5 +3,5 @@
 !
       integer nfmx,imx,jmx,lmx,nvarmx
-      parameter(nfmx=10,imx=200,jmx=150,lmx=20,nvarmx=1000)
+      parameter(nfmx=10,imx=200,jmx=150,lmx=200,nvarmx=1000)
 
       real xd(imx,nfmx),yd(jmx,nfmx),zd(lmx,nfmx),dtime(nfmx)
Index: LMDZ4/branches/pre_V3/libf/dyn3d/grid_noro.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/grid_noro.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/grid_noro.F	(revision 701)
@@ -53,4 +53,5 @@
 c             zval:  Minimum altitude
 C=======================================================================
+
       IMPLICIT INTEGER (I,J)
       IMPLICIT REAL(X,Z) 
@@ -79,4 +80,5 @@
       REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
       REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
+      REAL zmea0(imar+1,jmar) ! GK211005 (CG)
       REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
       REAL zpic(imar+1,jmar),zval(imar+1,jmar)
@@ -309,4 +311,5 @@
 C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
 
+       zmea0(:,:) = zmea(:,:) ! GK211005 (CG) on sauvegarde la topo non lissee
        CALL MVA9(zmea,iim+1,jjm+1)
        CALL MVA9(zstd,iim+1,jjm+1)
@@ -353,5 +356,8 @@
 c angle theta:
            zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask_tmp(ii,jj)
-           zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
+           ! GK211005 (CG) ne pas forcement lisser la topo
+           ! zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
+           zphi(ii,jj)=zmea0(ii,jj)*mask_tmp(ii,jj)
+           !
            zmea(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
            zpic(ii,jj)=zpic(ii,jj)*mask_tmp(ii,jj)
Index: LMDZ4/branches/pre_V3/libf/dyn3d/grilles_gcm_netcdf.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/grilles_gcm_netcdf.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/grilles_gcm_netcdf.F	(revision 701)
@@ -73,6 +73,6 @@
 
       do i=1,iip1
-         rlonudeg(i)=rlonu(i)*180./pi
-         rlonvdeg(i)=rlonv(i)*180./pi
+         rlonudeg(i)=rlonu(i)*180./pi + 360.
+         rlonvdeg(i)=rlonv(i)*180./pi + 360.
       enddo
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/guide.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/guide.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/guide.F	(revision 701)
@@ -81,10 +81,10 @@
       INTEGER step_rea,count_no_rea
 
-      real aire_min,aire_max
+cIM 180305   real aire_min,aire_max
       integer ilon,ilat
       real factt,ztau(ip1jmp1)
 
       INTEGER itau,ij,l,i,j
-      integer ncidt,varidpl,nlev,status
+      integer ncidpl,varidpl,nlev,status
       integer rcod,rid 
       real ditau,tau,a
@@ -208,16 +208,32 @@
          step_rea=1
          count_no_rea=0
+         ncidpl=-99
 
 c    itau_test    montre si l'importation a deja ete faite au rang itau
 c lecture d'un fichier netcdf pour determiner le nombre de niveaux
-         ncidt=NCOPN('T.nc',NCNOWRIT,rcod)
+         if (guide_u) then
+           if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
+         endif
+c
+         if (guide_v) then
+           if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
+         endif
+c
+         if (guide_T) then
+           if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
+         endif
+c
+         if (guide_Q) then
+           if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
+         endif
+c
          if (ncep) then
-          status=NF_INQ_DIMID(ncidt,'LEVEL',rid)
+          status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
          else
-          status=NF_INQ_DIMID(ncidt,'PRESSURE',rid)
-         endif
-          status=NF_INQ_DIMLEN(ncidt,rid,nlev)
+          status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
+         endif
+          status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
          print *,'nlev', nlev 
-          call ncclos(ncidt,rcod)
+          call ncclos(ncidpl,rcod)
 c   Lecture du premier etat des reanalyses.
          call read_reanalyse(1,ps
@@ -296,6 +312,4 @@
       tau=tau-aint(tau)
 
-      print*,'ATTENTION !!!! ON NE GUIDE QUE JUSQU A 15N'
-
 c  ucov
       if (guide_u) then
@@ -485,6 +499,9 @@
                zlat=rlatv(j)*180./pi
             endif
+         if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
+c  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
+             alpha(i,j)=alphamin
+         else
             xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
-c  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
             xi=min(xi,1.)
             if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
@@ -493,4 +510,5 @@
                alpha(i,j)=0.
             endif
+         endif
          enddo
       enddo
Index: LMDZ4/branches/pre_V3/libf/dyn3d/ini_paramLMDZ_dyn.h
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/ini_paramLMDZ_dyn.h	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/ini_paramLMDZ_dyn.h	(revision 701)
@@ -0,0 +1,259 @@
+c
+      dt_cum = dtvr*day_step
+
+      zan = annee_ref
+      dayref = day_ref
+      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
+      tau0 = itau_dyn
+c
+       pi = 4.0 * ATAN(1.0)
+       degres = 180./pi
+       rlong = rlonu * degres
+       rlatg = rlatu * degres
+c
+      CALL histbeg("paramLMDZ_dyn.nc", 
+     .                 iip1,rlong, jjp1,rlatg,
+     .                 1,1,1,1,
+     .                 tau0, zjulian, dt_cum,
+     .                 thoriid, nid_ctesGCM)
+c
+         CALL histdef(nid_ctesGCM, "prt_level", 
+     .        "Niveau impression debuggage dynamique",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "dayref", 
+     .        "Jour de l etat initial ( = 350  si 20 Decembre par ex.)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "anneeref", 
+     .        "Annee de l etat initial",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "anneelim", 
+     .        "Annee du fichier limitxxxx.nc  si  ok_limitvrai =y",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "raz_date", 
+     .   "Remise a zero (raz) date init.: 0 pas de raz;1=date gcm.def",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "nday", 
+     .   "Nombre de jours d integration",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "day_step", 
+     .   "nombre de pas par jour pour dt = 1 min",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iperiod", 
+     .   "periode pour le pas Matsuno (en pas de temps)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iapp_tracvl", 
+     .   "frequence du groupement des flux (en pas de temps)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iconser", 
+     .  "periode de sortie des variables de controle (en pas de temps)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iecri", 
+     .  "periode d ecriture du fichier histoire (en jour)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "periodav", 
+     .  "periode de stockage fichier histmoy (en jour)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "idissip", 
+     .  "periode de la dissipation (en pas) ... a completer",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "lstardis", 
+     .  "choix de l operateur de dissipation: 1= star,0=non-star ??",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "nitergdiv", 
+     .  "nombre d iterations de l operateur de dissipation gradiv",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "nitergrot", 
+     .  "nombre d iterations de l operateur de dissipation nxgradrot",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "niterh", 
+     .  "nombre d iterations de l operateur de dissipation divgrad",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "tetagdiv", 
+     ."temps dissipation des + petites long. d ondes pour u,v (gradiv)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "tetagrot", 
+     ."temps diss. des + petites long. d ondes pour u,v (nxgradrot)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "tetatemp", 
+     ."temps diss. des + petites long. d ondes pour h (divgrad)",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "coefdis", 
+     ."coefficient pour gamdissip",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "purmats", 
+     ."Choix schema integration temporel: 1=Matsuno,0=Matsuno-leapfrog",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ok_guide", 
+     ."Guidage: 1=true ,0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "true_calendar", 
+     ."Choix du calendrier: 1=gregorien ,0=calen. a 360 j",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "guide_calend", 
+     ."Guidage calendrier gregorien: 1=oui ,0=non",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iflag_phys", 
+     ."Permet de faire tourner le modele sans physique: 1=avec ,0=sans",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iphysiq", 
+     ."Periode de la physique en pas de temps de la dynamique",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ecritphy", 
+     ."frequence  de l ecriture du fichier histphy en jours",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "cycle_diurne", 
+     ."Avoir un cycle diurne: 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "soil_model", 
+     ."Choix du modele de sol (Thermique ?): 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "new_oliq", 
+     ."Nouvelle eau liquide: 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ok_orodr", 
+     ."Orography drag : 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ok_orolf", 
+     ."Orography lift : 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ok_limitvrai", 
+     ."Force la lecture de la bonne annee: 1=true, 0=false",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "nbapp_rad", 
+     ."Nombre  d appels des routines de rayonnements par jour",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "iflag_con", 
+     ."Flag  pour la convection: 1=LMD, 2=Tiedtke, 3=KE, 4=KE vect.",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "clon", 
+     ."longitude en degres du centre du zoom",
+     .                "deg",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "clat", 
+     ."latitude en degres du centre du zoom",
+     .                "deg",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "grossismx", 
+     ."facteur de grossissement du zoom, selon la longitude",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "grossismy", 
+     ."facteur de grossissement du zoom, selon la latitude",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "fxyhypb", 
+     ."Fonction f(y) hyperbolique  si true=1, sinusoidale si false=0",
+     .                "-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "dzoomx", 
+     ."extension en longitude de la zone du zoom (fraction zone totale)"
+     .                ,"-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "dzoomy", 
+     ."extension en latitude de la zone du zoom (fraction zone totale)"
+     .                ,"-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "taux", 
+     ."raideur du zoom en  X"
+     .                ,"-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "tauy", 
+     ."raideur du zoom en  Y"
+     .                ,"-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ysinus", 
+     ."ysinus=1: Ftion f(y) avec y=Sin(latit.)/ ysinus=0: y = latit"
+     .                ,"-",iip1,jjp1,thoriid, 1,1,1, -99, 32, 
+     .                "once", dt_cum,dt_cum)
+c
+         CALL histdef(nid_ctesGCM, "ip_ebil_dyn", 
+     ."PRINTlevel for energy conservation diag.; 0/1= pas de print,
+     . 2= print","-",iip1,jjp1,thoriid, 1,1,1, -99, 32,
+     .                "once", dt_cum,dt_cum)
+c
+c=================================================================
+c
+         CALL histend(nid_ctesGCM)
+c
+c=================================================================
Index: LMDZ4/branches/pre_V3/libf/dyn3d/iniadvtrac.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/iniadvtrac.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/iniadvtrac.F	(revision 701)
@@ -8,4 +8,5 @@
 #ifdef INCA
       USE transport_controls, only : hadv_flg, vadv_flg
+      USE species_names
       USE chemshut
 #endif
@@ -125,25 +126,43 @@
      s     iostat=ierr)
       if(ierr.eq.0) then
-       print*,'ouverture de traceur.def ok'
-       read(90,*) nq
-       print*,'nombre de traceurs ',nq
-       if (nq.gt.nqmx) then
-       print*,'nombre de traceurs trop important'
-       print*,'verifier traceur.def'
-       stop
-       endif
+        print*,'ouverture de traceur.def ok'
+        read(90,*) nq
+        print*,'nombre de traceurs ',nq
+        if (nq.gt.nqmx) then
+          print*,'nombre de traceurs trop important'
+          print*,'verifier traceur.def'
+          stop
+        endif
 C
-       do iq=1,nq
-         read(90,999) hadv(iq),vadv(iq),tnom(iq)
-       end do
-       close(90)  
-       PRINT*,'lecture de traceur.def :'   
-       do iq=1,nq
-         write(*,*) hadv(iq),vadv(iq),tnom(iq)
-       end do       
+        do iq=1,nq
+          read(90,999) hadv(iq),vadv(iq),tnom(iq)
+        end do
+        close(90)  
+        PRINT*,'lecture de traceur.def :'   
+        do iq=1,nq
+          write(*,*) hadv(iq),vadv(iq),tnom(iq)
+        end do       
       else
-       print*,'pb ouverture traceur.def'
-       print*,'ATTENTION on prend des valeurs par defaut'
-      endif
+        print*,'pb ouverture traceur.def'
+        print*,'ATTENTION on prend des valeurs par defaut'
+        nq = 4
+        hadv(1) = 14
+        vadv(1) = 14
+        tnom(1) = 'H2Ov'
+        hadv(2) = 10
+        vadv(2) = 10
+        tnom(2) = 'H2Ol'
+        hadv(3) = 10
+        vadv(3) = 10
+        tnom(3) = 'RN'
+        hadv(4) = 10
+        vadv(4) = 10
+        tnom(4) = 'PB'
+      ENDIF
+      PRINT*,'Valeur de traceur.def :'
+      do iq=1,nq
+        write(*,*) hadv(iq),vadv(iq),tnom(iq)
+      end do       
+
 #endif
 c a partir du nom court du traceur et du schema d'advection au detemine le nom long.
@@ -174,5 +193,9 @@
          str1=tnom(iq)
          tname(iiq)=tnom(iq)
-         ttext(iiq)=str1(1:lnblnk(str1))//descrq(iadv(iiq))
+         IF (iadv(iiq).eq.0) THEN
+           ttext(iiq)=str1(1:lnblnk(str1))
+         ELSE
+           ttext(iiq)=str1(1:lnblnk(str1))//descrq(iadv(iiq))
+         endif
          str2=ttext(iiq)
 c   schemas tenant compte des moments d'ordre superieur.
Index: LMDZ4/branches/pre_V3/libf/dyn3d/iniinterp_horiz.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/iniinterp_horiz.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/iniinterp_horiz.F	(revision 701)
@@ -0,0 +1,179 @@
+C 
+C $Header$
+C
+      subroutine iniinterp_horiz (imo,jmo,imn,jmn ,kllm,
+     &       rlonuo,rlatvo,rlonun,rlatvn,
+     &       ktotal,iik,jjk,jk,ik,intersec,airen)
+   
+      implicit none
+
+
+
+c ---------------------------------------------------------
+c Prepare l' interpolation des variables d'une grille LMDZ
+c  dans une autre grille LMDZ en conservant la quantite
+c  totale pour les variables intensives (/m2) : ex : Pression au sol
+c
+c   (Pour chaque case autour d'un point scalaire de la nouvelle
+c    grille, on calcule la surface (en m2)en intersection avec chaque
+c    case de l'ancienne grille , pour la future interpolation)
+c
+c on calcule aussi l' aire dans la nouvelle grille 
+c
+c
+c   Auteur:  F.Forget 01/1995
+c   -------
+c
+c ---------------------------------------------------------
+c   Declarations:
+c ==============
+c
+c  ARGUMENTS
+c  """""""""
+c INPUT
+       integer imo, jmo ! dimensions ancienne grille
+       integer imn,jmn  ! dimensions nouvelle grille
+       integer kllm ! taille du tableau des intersections
+       real rlonuo(imo+1)     !  Latitude et
+       real rlatvo(jmo)       !  longitude des
+       real rlonun(imn+1)     !  bord des
+       real rlatvn(jmn)     !  cases "scalaires" (input)
+
+c OUTPUT
+       integer ktotal ! nombre totale d'intersections reperees
+       integer iik(kllm), jjk(kllm),jk(kllm),ik(kllm)
+       real intersec(kllm)  ! surface des intersections (m2)
+       real airen (imn+1,jmn+1) ! aire dans la nouvelle grille
+
+
+       
+ 
+c Autres variables
+c """"""""""""""""
+       integer i,j,ii,jj,k
+       real a(imo+1),b(imo+1),c(jmo+1),d(jmo+1)
+       real an(imn+1),bn(imn+1),cn(jmn+1),dn(jmn+1)
+       real aa, bb,cc,dd
+       real pi
+
+       pi      = 2.*ASIN( 1. )
+
+
+
+c On repere les frontieres des cases :
+c =================================== 
+c
+c Attention, on ruse avec des latitudes = 90 deg au pole.
+
+
+c  ANcienne grile
+c  """"""""""""""
+      a(1) =   - rlonuo(imo+1)
+      b(1) = rlonuo(1)
+      do i=2,imo+1
+         a(i) = rlonuo(i-1)
+         b(i) =  rlonuo(i)
+      end do
+
+      d(1) = pi/2 
+      do j=1,jmo
+         c(j) = rlatvo(j) 
+         d(j+1) = rlatvo(j)
+      end do
+      c(jmo+1) = -pi/2 
+      
+
+c  Nouvelle grille
+c  """""""""""""""
+      an(1) =  - rlonun(imn+1)
+      bn(1) = rlonun(1)
+      do i=2,imn+1
+         an(i) = rlonun(i-1)
+         bn(i) =  rlonun(i)
+      end do
+
+      dn(1) = pi/2 
+      do j=1,jmn
+         cn(j) = rlatvn(j)
+         dn(j+1) = rlatvn(j)
+      end do
+      cn(jmn+1) = -pi/2 
+
+c Calcul de la surface des cases scalaires de la nouvelle grille
+c ==============================================================
+      do ii=1,imn + 1
+        do jj = 1,jmn+1
+               airen(ii,jj) = (bn(ii)-an(ii))*(sin(dn(jj))-sin(cn(jj)))
+        end do
+      end do
+
+c Calcul de la surface des intersections
+c ======================================
+
+c     boucle sur la nouvelle grille
+c     """"""""""""""""""""""""""""
+      ktotal = 0
+      do jj = 1,jmn+1
+       do j=1, jmo+1
+          if((cn(jj).lt.d(j)).and.(dn(jj).gt.c(j)))then
+              do ii=1,imn + 1
+                do i=1, imo +1
+                    if (  ((an(ii).lt.b(i)).and.(bn(ii).gt.a(i)))
+     &        .or. ((an(ii).lt.b(i)-2*pi).and.(bn(ii).gt.a(i)-2*pi)
+     &             .and.(b(i)-2*pi.lt.-pi) )
+     &        .or. ((an(ii).lt.b(i)+2*pi).and.(bn(ii).gt.a(i)+2*pi)
+     &             .and.(a(i)+2*pi.gt.pi) )
+     &                     )then
+                      ktotal = ktotal +1
+                      iik(ktotal) =ii
+                      jjk(ktotal) =jj
+                      ik(ktotal) =i
+                      jk(ktotal) =j
+
+                      dd = min(d(j), dn(jj))
+                      cc = cn(jj)
+                      if (cc.lt. c(j))cc=c(j)
+                      if((an(ii).lt.b(i)-2*pi).and.
+     &                  (bn(ii).gt.a(i)-2*pi)) then 
+                          bb = min(b(i)-2*pi,bn(ii))
+                          aa = an(ii)
+                          if (aa.lt.a(i)-2*pi) aa=a(i)-2*pi
+                      else if((an(ii).lt.b(i)+2*pi).and.
+     &                       (bn(ii).gt.a(i)+2*pi)) then
+                          bb = min(b(i)+2*pi,bn(ii))
+                          aa = an(ii)
+                          if (aa.lt.a(i)+2*pi) aa=a(i)+2*pi
+                      else 
+                          bb = min(b(i),bn(ii))
+                          aa = an(ii)
+                          if (aa.lt.a(i)) aa=a(i)
+                      end if
+                      intersec(ktotal)=(bb-aa)*(sin(dd)-sin(cc))
+                     end if
+                end do
+               end do
+             end if
+         end do
+       end do       
+
+
+
+c     TEST  INFO
+c     DO k=1,ktotal 
+c      ii = iik(k) 
+c      jj = jjk(k)
+c      i = ik(k)
+c      j = jk(k)
+c      if ((ii.eq.10).and.(jj.eq.10).and.(i.eq.10).and.(j.eq.10))then
+c      if (jj.eq.2.and.(ii.eq.1))then
+c          write(*,*) '**************** jj=',jj,'ii=',ii
+c          write(*,*) 'i,j =',i,j
+c          write(*,*) 'an,bn,cn,dn', an(ii), bn(ii), cn(jj),dn(jj)
+c          write(*,*) 'a,b,c,d', a(i), b(i), c(j),d(j)
+c          write(*,*) 'intersec(k)',intersec(k)
+c          write(*,*) 'airen(ii,jj)=',airen(ii,jj)
+c      end if
+c     END DO 
+
+      return
+      end
Index: LMDZ4/branches/pre_V3/libf/dyn3d/interp_horiz.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/interp_horiz.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/interp_horiz.F	(revision 701)
@@ -0,0 +1,154 @@
+c
+c $Header$
+c
+      subroutine interp_horiz (varo,varn,imo,jmo,imn,jmn,lm,
+     &  rlonuo,rlatvo,rlonun,rlatvn)  
+
+c===========================================================
+c  Interpolation Horizontales des variables d'une grille LMDZ
+c (des points SCALAIRES au point SCALAIRES)
+c  dans une autre grille LMDZ en conservant la quantite
+c  totale pour les variables intensives (/m2) : ex : Pression au sol
+c
+c Francois Forget (01/1995)
+c===========================================================
+
+      IMPLICIT NONE 
+
+c   Declarations:
+c ==============
+c
+c  ARGUMENTS
+c  """""""""
+        
+       integer imo, jmo ! dimensions ancienne grille (input)
+       integer imn,jmn  ! dimensions nouvelle grille (input)
+
+       real rlonuo(imo+1)     !  Latitude et
+       real rlatvo(jmo)       !  longitude des
+       real rlonun(imn+1)     !  bord des 
+       real rlatvn(jmn)     !  cases "scalaires" (input)
+
+       integer lm ! dimension verticale (input)
+       real varo (imo+1, jmo+1,lm) ! var dans l'ancienne grille (input)
+       real varn (imn+1,jmn+1,lm) ! var dans la nouvelle grille (output)
+
+c Autres variables
+c """"""""""""""""
+       real airetest(imn+1,jmn+1)
+       integer ii,jj,l
+
+       real airen (imn+1,jmn+1) ! aire dans la nouvelle grille
+c    Info sur les ktotal intersection entre les cases new/old grille
+       integer kllm, k, ktotal
+       parameter (kllm = 400*200*10)
+       integer iik(kllm), jjk(kllm),jk(kllm),ik(kllm)
+       real intersec(kllm)
+       real R
+       real totn, tots
+
+       logical firstcall, firsttest, aire_ok
+       save firsttest
+       data firsttest /.true./
+       data aire_ok /.true./
+
+       
+
+
+
+c initialisation
+c --------------
+c Si c'est le premier appel, on prepare l'interpolation
+c en calculant pour chaque case autour d'un point scalaire de la
+c nouvelle grille, la surface  de intersection avec chaque
+c    case de l'ancienne grille.
+
+
+        call iniinterp_horiz (imo,jmo,imn,jmn ,kllm,
+     &       rlonuo,rlatvo,rlonun,rlatvn,
+     &          ktotal,iik,jjk,jk,ik,intersec,airen)
+
+      do l=1,lm
+       do jj =1 , jmn+1
+        do ii=1, imn+1
+          varn(ii,jj,l) =0.
+        end do
+       end do
+      end do 
+       
+c Interpolation horizontale
+c -------------------------
+c boucle sur toute les ktotal intersections entre les cases
+c de l'ancienne et la  nouvelle grille
+c
+      PRINT *, 'ktotal 1 = ', ktotal
+     
+      do k=1,ktotal
+        do l=1,lm
+         varn(iik(k),jjk(k),l) = varn(iik(k),jjk(k),l) 
+     &        + varo(ik(k), jk(k),l)*intersec(k)/airen(iik(k),jjk(k))
+        end do
+      end do
+
+c Une seule valeur au pole pour les variables ! :
+c -----------------------------------------------
+       do l=1, lm
+         totn =0.
+         tots =0.
+           do ii =1, imn+1
+             totn = totn + varn(ii,1,l)
+             tots = tots + varn (ii,jmn+1,l)
+           end do 
+           do ii =1, imn+1
+             varn(ii,1,l) = totn/float(imn+1)
+             varn(ii,jmn+1,l) = tots/float(imn+1)
+           end do 
+       end do
+           
+
+c---------------------------------------------------------------
+c  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST  TEST 
+!!       if (.not.(firsttest)) goto 99
+!!       firsttest = .false.
+!! !     write (*,*) 'INTERP. HORIZ. : TEST SUR LES AIRES:'
+!!       do jj =1 , jmn+1
+!!         do ii=1, imn+1
+!!           airetest(ii,jj) =0.
+!!         end do
+!!       end do 
+!!       PRINT *, 'ktotal = ', ktotal
+!!       PRINT *, 'jmn+1 =', jmn+1, 'imn+1', imn+1
+!! 
+!!       do k=1,ktotal
+!!          airetest(iik(k),jjk(k))= airetest(iik(k),jjk(k)) +intersec(k) 
+!!       end DO
+!! 
+!! 
+!!       PRINT *, 'fin boucle'
+!!       do jj =1 , jmn+1
+!!        do ii=1, imn+1
+!!          r = airen(ii,jj)/airetest(ii,jj)
+!!          if ((r.gt.1.001).or.(r.lt.0.999)) then
+!! !             write (*,*) '********** PROBLEME D'' AIRES !!!',
+!! !     &                   ' DANS L''INTERPOLATION HORIZONTALE'
+!! !             write(*,*)'ii,jj,airen,airetest',
+!! !     &          ii,jj,airen(ii,jj),airetest(ii,jj)
+!!              aire_ok = .false.
+!!          end if
+!!        end do
+!!       end do
+!! !      if (aire_ok) write(*,*) 'INTERP. HORIZ. : AIRES OK'
+!!  99   continue
+
+c FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST  FIN TEST
+c---------------------------------------------------------------
+
+
+
+
+
+
+
+
+        return
+        end
Index: LMDZ4/branches/pre_V3/libf/dyn3d/leapfrog.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/leapfrog.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/leapfrog.F	(revision 701)
@@ -4,9 +4,13 @@
 c
 c
-      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0)
+      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
+     &                    time_0)
 
 #ifdef INCA
       USE transport_controls, ONLY : hadv_flg, mmt_adj
 #endif
+
+cIM : pour sortir les param. du modele dans un fis. netcdf 110106
+      USE IOIPSL 
 
       IMPLICIT NONE
@@ -62,4 +66,6 @@
 #include "academic.h"
 
+#include "clesphys.h"
+
       integer nq
 
@@ -117,5 +123,6 @@
       REAL time_0 , finvmaold(ip1jmp1,llm)
 
-      LOGICAL lafin
+cym      LOGICAL  lafin
+      LOGICAL :: lafin=.false.
       INTEGER ij,iq,l
       INTEGER ik
@@ -125,4 +132,18 @@
       REAL rdayvrai,rdaym_ini
       LOGICAL first,callinigrads
+cIM : pour sortir les param. du modele dans un fis. netcdf 110106
+      save first
+      data first/.true./
+      real dt_cum, zjulian
+      character*10 infile
+      integer zan, tau0, thoriid
+      integer nid_ctesGCM
+      save nid_ctesGCM
+      real degres
+      real rlong(iip1), rlatg(jjp1)
+      real zx_tmp_2d(iip1,jjp1)
+      integer ndex2d(iip1*jjp1)
+      logical ok_sync
+      parameter (ok_sync = .true.) 
 
       data callinigrads/.true./
@@ -130,5 +151,5 @@
 
       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
-#ifdef INCA_CH4
+#ifdef INCA
       REAL :: flxw(ip1jmp1,llm)
 #endif
@@ -144,11 +165,8 @@
       REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
       CHARACTER*15 ztit
-      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
-      SAVE      ip_ebil_dyn
-      DATA      ip_ebil_dyn/0/
+!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
+!IM   SAVE      ip_ebil_dyn
+!IM   DATA      ip_ebil_dyn/0/
 c-jld 
-
-      LOGICAL offline  ! Controle du stockage ds "fluxmass"
-      PARAMETER (offline=.false.)
 
       character*80 dynhist_file, dynhistave_file
@@ -199,6 +217,9 @@
 
 #ifdef CPP_IOIPSL
-      if (ok_guide) then
+      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
         call guide(itau,ucov,vcov,teta,q,masse,ps)
+      else
+        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
+     .    '6 dernieres heures'
       endif
 #endif
@@ -272,5 +293,5 @@
 
 c
-#ifdef INCA_CH4
+#ifdef INCA
              CALL caladvtrac(q,pbaru,pbarv,
      *                      p, masse, dq,  teta,
@@ -289,6 +310,6 @@
 
 #ifdef CPP_IOIPSL
-c           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
-c    .   time_step, itau)
+           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
+     .   dtvr, itau)
 #endif
 
@@ -352,9 +373,17 @@
       ENDIF 
 c-jld
+cIM : pour sortir les param. du modele dans un fis. netcdf 110106
+      IF (first) THEN
+       first=.false.
+#include "ini_paramLMDZ_dyn.h"
+      ENDIF
+c
+#include "write_paramLMDZ_dyn.h"
+c
 
         CALL calfis( nq, lafin ,rdayvrai,time  ,
      $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
      $               du,dv,dteta,dq,w,
-#ifdef INCA_CH4
+#ifdef INCA
      $               flxw,
 #endif
@@ -486,11 +515,11 @@
 
             IF( itau. EQ. itaufinp1 ) then  
-c      write(79,*) 'ucov',ucov
-c      write(80,*) 'vcov',vcov
-c      write(81,*) 'teta',teta
-c      write(82,*) 'ps',ps
-c      write(83,*) 'q',q
-c      WRITE(85,*) 'q1 = ',q(:,:,1)
-c      WRITE(86,*) 'q3 = ',q(:,:,3)
+c$$$       write(79,*) 'ucov',ucov
+c$$$       write(80,*) 'vcov',vcov
+c$$$       write(81,*) 'teta',teta
+c$$$       write(82,*) 'ps',ps
+c$$$       write(83,*) 'q',q
+c$$$       WRITE(85,*) 'q1 = ',q(:,:,1)
+c$$$       WRITE(86,*) 'q3 = ',q(:,:,3)
 
               abort_message = 'Simulation finished'
Index: LMDZ4/branches/pre_V3/libf/dyn3d/limit_netcdf.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/limit_netcdf.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/limit_netcdf.F	(revision 701)
@@ -452,7 +452,10 @@
 
       ierr = NF_OPEN('amipbc_sic_1x1.nc', NF_NOWRITE, ncid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
+      if (ierr.ne.0) THEN
+        ierr = NF_OPEN('amipbc_sic_1x1_clim.nc', NF_NOWRITE, ncid)
+        if (ierr.ne.0) THEN
+          print *, NF_STRERROR(ierr)
+          STOP
+        endif
       ENDIF
 
@@ -775,7 +778,10 @@
 c     ierr = NF_OPEN('AMIP_SST.nc', NF_NOWRITE, ncid)
       ierr = NF_OPEN('amipbc_sst_1x1.nc', NF_NOWRITE, ncid)
-      if (ierr.ne.0) then
-        print *, NF_STRERROR(ierr)
-        STOP
+      if (ierr.ne.0) THEN
+        ierr = NF_OPEN('amipbc_sst_1x1_clim.nc', NF_NOWRITE, ncid)
+        if (ierr.ne.0) THEN
+          print *, NF_STRERROR(ierr)
+          STOP
+        endif
       ENDIF
 
Index: LMDZ4/branches/pre_V3/libf/dyn3d/parallel.F90
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/parallel.F90	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/parallel.F90	(revision 701)
@@ -0,0 +1,463 @@
+  module parallel
+  
+    integer, save :: mpi_size
+    integer, save :: mpi_rank
+    integer, save :: jj_begin
+    integer, save :: jj_end
+    integer, save :: jj_nb
+    integer, save :: ij_begin
+    integer, save :: ij_end
+    logical, save :: pole_nord
+    logical, save :: pole_sud
+    
+    integer, allocatable, save, dimension(:) :: jj_begin_para
+    integer, allocatable, save, dimension(:) :: jj_end_para
+    integer, allocatable, save, dimension(:) :: jj_nb_para
+    
+ contains
+
+    subroutine init_parallel
+
+#ifdef CPP_PARALLEL
+
+    USE vampir
+    implicit none
+    
+      integer :: ierr
+      integer :: i,j
+      integer :: type_size
+      integer, dimension(3) :: blocklen,type
+      
+      
+      include 'mpif.h'
+#include "dimensions90.h"
+#include "paramet90.h"
+      
+      call MPI_INIT(ierr)
+      call InitVampir
+      call MPI_COMM_SIZE(MPI_COMM_WORLD,mpi_size,ierr)
+      call MPI_COMM_RANK(MPI_COMM_WORLD,mpi_rank,ierr)
+  
+      
+      allocate(jj_begin_para(0:mpi_size-1))
+      allocate(jj_end_para(0:mpi_size-1))
+      allocate(jj_nb_para(0:mpi_size-1))
+      
+      do i=0,mpi_size-1
+        jj_nb_para(i)=(jjm+1)/mpi_size
+        if ( i < MOD((jjm+1),mpi_size) ) jj_nb_para(i)=jj_nb_para(i)+1
+        
+        if (jj_nb_para(i) <= 2 ) then
+          
+         print *,"Arret : le nombre de bande de lattitude par process est trop faible (<2)."
+          print *," ---> diminuez le nombre de CPU ou augmentez la taille en lattitude"
+          
+          call MPI_ABORT(MPI_COMM_WORLD,-1, ierr)
+          
+        endif
+        
+      enddo
+      
+!      jj_nb_para(0)=11
+!      jj_nb_para(1)=25
+!      jj_nb_para(2)=25
+!      jj_nb_para(3)=12      
+
+      j=1
+      
+      do i=0,mpi_size-1 
+        
+        jj_begin_para(i)=j
+        jj_end_para(i)=j+jj_Nb_para(i)-1
+        j=j+jj_Nb_para(i)
+      
+      enddo
+      
+      jj_begin = jj_begin_para(mpi_rank)
+      jj_end   = jj_end_para(mpi_rank)
+      jj_nb    = jj_nb_para(mpi_rank)
+      
+      ij_begin=(jj_begin-1)*iip1+1
+      ij_end=jj_end*iip1
+      
+      if (mpi_rank.eq.0) then
+        pole_nord=.TRUE.
+      else 
+        pole_nord=.FALSE.
+      endif
+      
+      if (mpi_rank.eq.mpi_size-1) then
+        pole_sud=.TRUE.
+      else 
+        pole_sud=.FALSE.
+      endif
+        
+      print *,"jj_begin",jj_begin
+      print *,"jj_end",jj_end
+      print *,"ij_begin",ij_begin
+      print *,"ij_end",ij_end
+    
+#endif
+    
+    end subroutine init_parallel
+
+    
+    subroutine SetDistrib(jj_Nb_New)
+
+#ifdef CPP_PARALLEL
+
+    implicit none
+
+#include "dimensions90.h"
+#include "paramet90.h"
+
+      INTEGER,dimension(0:MPI_Size-1) :: jj_Nb_New
+      INTEGER :: i  
+  
+      jj_Nb_Para=jj_Nb_New
+      
+      jj_begin_para(0)=1
+      jj_end_para(0)=jj_Nb_Para(0)
+      
+      do i=1,mpi_size-1 
+        
+        jj_begin_para(i)=jj_end_para(i-1)+1
+        jj_end_para(i)=jj_begin_para(i)+jj_Nb_para(i)-1
+      
+      enddo
+      
+      jj_begin = jj_begin_para(mpi_rank)
+      jj_end   = jj_end_para(mpi_rank)
+      jj_nb    = jj_nb_para(mpi_rank)
+      
+      ij_begin=(jj_begin-1)*iip1+1
+      ij_end=jj_end*iip1
+
+#endif
+
+    end subroutine SetDistrib
+    
+    subroutine Finalize_parallel
+
+#ifdef CPP_PARALLEL
+
+    implicit none
+
+#include "dimensions90.h"
+#include "paramet90.h"
+      integer :: ierr
+      integer :: i
+      include 'mpif.h'
+      
+      deallocate(jj_begin_para)
+      deallocate(jj_end_para)
+      deallocate(jj_nb_para)
+      
+      call MPI_FINALIZE(ierr)
+ 
+#endif
+     
+    end subroutine Finalize_parallel
+    
+    subroutine Pack_Data(Field,ij,ll,row,Buffer)
+    implicit none
+
+#include "dimensions90.h"
+#include "paramet90.h"
+
+      integer, intent(in) :: ij,ll,row
+      real,dimension(ij,ll),intent(in) ::Field
+      real,dimension(ll*iip1*row), intent(out) :: Buffer 
+            
+      integer :: Pos
+      integer :: i,l
+      
+      Pos=0
+      do l=1,ll
+        do i=1,row*iip1
+          Pos=Pos+1
+          Buffer(Pos)=Field(i,l)
+        enddo
+      enddo
+      
+    end subroutine Pack_data 
+     
+    subroutine Unpack_Data(Field,ij,ll,row,Buffer)
+    implicit none
+
+#include "dimensions90.h"
+#include "paramet90.h"
+
+      integer, intent(in) :: ij,ll,row
+      real,dimension(ij,ll),intent(out) ::Field
+      real,dimension(ll*iip1*row), intent(in) :: Buffer 
+            
+      integer :: Pos
+      integer :: i,l
+      
+      Pos=0
+      
+      do l=1,ll
+        do i=1,row*iip1
+          Pos=Pos+1
+          Field(i,l)=Buffer(Pos)
+        enddo
+      enddo
+      
+    end subroutine UnPack_data
+      
+    subroutine exchange_hallo(Field,ij,ll,up,down)
+
+#ifdef CPP_PARALLEL
+
+    USE Vampir
+    implicit none
+#include "dimensions90.h"
+#include "paramet90.h"    
+    include 'mpif.h'
+    
+      INTEGER :: ij,ll
+      REAL, dimension(ij,ll) :: Field
+      INTEGER :: up,down
+      
+      INTEGER :: ierr
+      LOGICAL :: SendUp,SendDown
+      LOGICAL :: RecvUp,RecvDown
+      INTEGER, DIMENSION(4) :: Request
+      INTEGER, DIMENSION(MPI_STATUS_SIZE,4) :: Status
+      INTEGER :: NbRequest
+      REAL, dimension(:),allocatable :: Buffer_Send_up,Buffer_Send_down
+      REAL, dimension(:),allocatable :: Buffer_Recv_up,Buffer_Recv_down
+      INTEGER :: Buffer_size      
+      
+      call MPI_Barrier(MPI_COMM_WORLD,ierr)
+      call VTb(VThallo)
+      
+      SendUp=.TRUE.
+      SendDown=.TRUE.
+      RecvUp=.TRUE.
+      RecvDown=.TRUE.
+        
+      IF (pole_nord) THEN
+        SendUp=.FALSE.
+        RecvUp=.FALSE.
+      ENDIF
+  
+      IF (pole_sud) THEN
+        SendDown=.FALSE.
+        RecvDown=.FALSE.
+      ENDIF
+      
+      if (up.eq.0) then
+        SendDown=.FALSE.
+        RecvUp=.FALSE.
+      endif
+      
+      if (down.eq.0) then
+        SendUp=.FALSE.
+        RecvDown=.FALSE.
+      endif
+      
+      NbRequest=0
+  
+      IF (SendUp) THEN
+        NbRequest=NbRequest+1
+        buffer_size=down*iip1*ll
+        allocate(Buffer_Send_up(Buffer_size))
+        call PACK_Data(Field(ij_begin,1),ij,ll,down,Buffer_Send_up)
+        call MPI_ISSEND(Buffer_send_up,Buffer_Size,MPI_REAL8,MPI_Rank-1,1,     &
+                        MPI_COMM_WORLD,Request(NbRequest),ierr)
+      ENDIF
+  
+      IF (SendDown) THEN
+        NbRequest=NbRequest+1
+        
+        buffer_size=up*iip1*ll
+        allocate(Buffer_Send_down(Buffer_size))
+        call PACK_Data(Field(ij_end+1-up*iip1,1),ij,ll,up,Buffer_send_down)
+        
+        call MPI_ISSEND(Buffer_send_down,Buffer_Size,MPI_REAL8,MPI_Rank+1,1,     &
+                        MPI_COMM_WORLD,Request(NbRequest),ierr)
+      ENDIF
+    
+  
+      IF (RecvUp) THEN
+        NbRequest=NbRequest+1
+        buffer_size=up*iip1*ll
+        allocate(Buffer_recv_up(Buffer_size))
+              
+        call MPI_IRECV(Buffer_recv_up,Buffer_size,MPI_REAL8,MPI_Rank-1,1,  &
+                        MPI_COMM_WORLD,Request(NbRequest),ierr)
+     
+       
+      ENDIF
+  
+      IF (RecvDown) THEN
+        NbRequest=NbRequest+1
+        buffer_size=down*iip1*ll
+        allocate(Buffer_recv_down(Buffer_size))
+        
+        call MPI_IRECV(Buffer_recv_down,Buffer_size,MPI_REAL8,MPI_Rank+1,1,     &
+                        MPI_COMM_WORLD,Request(NbRequest),ierr)
+      
+        
+      ENDIF
+  
+      if (NbRequest > 0) call MPI_WAITALL(NbRequest,Request,Status,ierr)
+      IF (RecvUp)  call Unpack_Data(Field(ij_begin-up*iip1,1),ij,ll,up,Buffer_Recv_up)
+      IF (RecvDown) call Unpack_Data(Field(ij_end+1,1),ij,ll,down,Buffer_Recv_down)  
+
+      call VTe(VThallo)
+      call MPI_Barrier(MPI_COMM_WORLD,ierr)
+      RETURN
+
+#endif
+      
+    end subroutine exchange_Hallo
+    
+    
+    subroutine Gather_Field(Field,ij,ll,rank)
+
+
+
+    implicit none
+#include "dimensions90.h"
+#include "paramet90.h"    
+
+#ifdef CPP_PARALLEL
+    include 'mpif.h'
+#endif
+    
+      INTEGER :: ij,ll,rank
+      REAL, dimension(ij,ll) :: Field
+#ifdef CPP_PARALLEL
+      REAL, dimension(:),allocatable :: Buffer_send   
+      REAL, dimension(:),allocatable :: Buffer_Recv
+      INTEGER, dimension(0:MPI_Size-1) :: Recv_count, displ
+      INTEGER :: ierr
+      INTEGER ::i
+      
+      if (ij==ip1jmp1) then 
+         allocate(Buffer_send(iip1*ll*(jj_end-jj_begin+1)))
+         call Pack_Data(Field(ij_begin,1),ij,ll,jj_end-jj_begin+1,Buffer_send)
+      else if (ij==ip1jm) then
+         allocate(Buffer_send(iip1*ll*(min(jj_end,jjm)-jj_begin+1)))
+         call Pack_Data(Field(ij_begin,1),ij,ll,min(jj_end,jjm)-jj_begin+1,Buffer_send)
+      else
+         print *,ij
+	 stop 'erreur dans Gather_Field'
+      endif
+      
+      if (MPI_Rank==rank) then
+        allocate(Buffer_Recv(ij*ll))
+        do i=0,MPI_Size-1
+          
+          if (ij==ip1jmp1) then 
+            Recv_count(i)=(jj_end_para(i)-jj_begin_para(i)+1)*ll*iip1
+          else if (ij==ip1jm) then
+            Recv_count(i)=(min(jj_end_para(i),jjm)-jj_begin_para(i)+1)*ll*iip1
+          else
+            stop 'erreur dans Gather_Field'
+          endif
+          
+          if (i==0) then 
+            displ(i)=0 
+          else
+            displ(i)=displ(i-1)+Recv_count(i-1)
+          endif
+          
+        enddo
+        
+      endif
+      
+      call MPI_GATHERV(Buffer_send,(min(ij_end,ij)-ij_begin+1)*ll,MPI_REAL8,   &
+                        Buffer_Recv,Recv_count,displ,MPI_REAL8,rank,MPI_COMM_WORLD,ierr)
+      
+      if (MPI_Rank==rank) then                  
+      
+        if (ij==ip1jmp1) then 
+          do i=0,MPI_Size-1
+            call Unpack_Data(Field((jj_begin_para(i)-1)*iip1+1,1),ij,ll,                 &
+                             jj_end_para(i)-jj_begin_para(i)+1,Buffer_Recv(displ(i)+1))
+          enddo
+        else if (ij==ip1jm) then
+          do i=0,MPI_Size-1
+             call Unpack_Data(Field((jj_begin_para(i)-1)*iip1+1,1),ij,ll,                       &
+                             min(jj_end_para(i),jjm)-jj_begin_para(i)+1,Buffer_Recv(displ(i)+1))
+          enddo
+        endif
+      
+      endif
+
+#endif
+      
+    end subroutine Gather_Field
+    
+    subroutine AllGather_Field(Field,ij,ll)
+
+#ifdef CPP_PARALLEL
+
+    implicit none
+#include "dimensions90.h"
+#include "paramet90.h"    
+    include 'mpif.h'
+    
+      INTEGER :: ij,ll
+      REAL, dimension(ij,ll) :: Field
+      INTEGER :: ierr
+      
+      call Gather_Field(Field,ij,ll,0)
+      call MPI_BCAST(Field,ij*ll,MPI_REAL8,0,MPI_COMM_WORLD)
+
+#endif
+      
+    end subroutine AllGather_Field
+    
+   subroutine Broadcast_Field(Field,ij,ll,rank)
+
+#ifdef CPP_PARALLEL
+
+    implicit none
+#include "dimensions90.h"
+#include "paramet90.h"    
+    include 'mpif.h'
+    
+      INTEGER :: ij,ll
+      REAL, dimension(ij,ll) :: Field
+      INTEGER :: rank
+      INTEGER :: ierr
+      
+      call MPI_BCAST(Field,ij*ll,MPI_REAL8,rank,MPI_COMM_WORLD)
+
+#endif
+      
+    end subroutine Broadcast_Field
+        
+   
+    /*  
+  Subroutine verif_hallo(Field,ij,ll,up,down)
+
+#ifdef CPP_PARALLEL
+
+    implicit none
+#include "dimensions90.h"
+#include "paramet90.h"    
+    include 'mpif.h'
+    
+      INTEGER :: ij,ll
+      REAL, dimension(ij,ll) :: Field
+      INTEGER :: up,down 
+      
+      REAL,dimension(ij,ll): NewField
+      
+      NewField=0
+      
+      ijb=ij_begin
+      ije=ij_end
+      if (pole_nord) 
+      NewField(ij_be       
+
+#endif
+
+*/
+  end module parallel
Index: LMDZ4/branches/pre_V3/libf/dyn3d/ppm3d.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/ppm3d.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/ppm3d.F	(revision 701)
@@ -737,5 +737,5 @@
 c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
       Q(I,  2,k,IC) = Q(I,  1,k,IC)
-      Q(I,JMR,k,IC) = Q(I,JMP,k,IC)
+      Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
 400   CONTINUE
       endif
Index: LMDZ4/branches/pre_V3/libf/dyn3d/read_reanalyse.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/read_reanalyse.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/read_reanalyse.F	(revision 701)
@@ -43,4 +43,5 @@
       integer varidpl,ncidQ,varidQ
       save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
+      save ncidpl
       save varidpl,ncidQ,varidQ
 
@@ -48,5 +49,5 @@
       real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
       real*4 Qnc(iip1,jjp1,nlevnc)
-      real*4 pl(nlevnc),presnc(iip1,jjp1,nlevnc)
+      real*4 pl(nlevnc)
 
       integer start(4),count(4),status
@@ -80,5 +81,5 @@
             varidv=NCVID(ncidv,'VWND',rcode)
             print*,'ncidv,varidv',ncidv,varidv
-            if (ncidpl.eq.-99) ncidpl=ncidu
+            if (ncidpl.eq.-99) ncidpl=ncidv
             endif
 
@@ -88,5 +89,5 @@
             varidt=NCVID(ncidt,'AIR',rcode)
             print*,'ncidt,varidt',ncidt,varidt
-            if (ncidpl.eq.-99) ncidpl=ncidu
+            if (ncidpl.eq.-99) ncidpl=ncidt
             endif
 
@@ -96,5 +97,5 @@
             varidQ=NCVID(ncidQ,'RH',rcode)
             print*,'ncidQ,varidQ',ncidQ,varidQ
-            if (ncidpl.eq.-99) ncidpl=ncidu
+            if (ncidpl.eq.-99) ncidpl=ncidQ
             endif
 
@@ -106,5 +107,5 @@
             endif
 
-c Coordonnee vertcale
+c Coordonnee verticale
             if (ncep) then
                print*,'Vous etes entrain de lire des donnees NCEP'
@@ -115,16 +116,8 @@
             endif
             print*,'ncidu,varidpl',ncidu,varidpl
-
-      endif
-
+      endif
       print*,'ok1'
 
-c -----------------------------------------------------------------
-c   lecture des champs u, v, T, ps
-c -----------------------------------------------------------------
-
-c  niveaux de pression
-c  -------------------
-
+c Niveaux de pression
       print*,'WARNING!!! Il n y a pas de test de coherence'
       print*,'sur le nombre de niveaux verticaux dans le fichier nc'
@@ -136,14 +129,13 @@
 c  passage en pascal
       pl(:)=100.*pl(:)
-c  passage des presions sur une grille 3D
-      do l=1,nlevnc
-         presnc(:,:,l)=pl(l)
-      enddo
-      if(first) then
-         do l=1,nlevnc
-            print*,'PL(',l,')=',pl(l)
-         enddo
-      endif
-
+      if (first) then
+       do l=1,nlevnc
+          print*,'PL(',l,')=',pl(l)
+       enddo
+      endif
+
+c -----------------------------------------------------------------
+c   lecture des champs u, v, T, ps
+c -----------------------------------------------------------------
 
 c  dimensions pour les champs scalaires et le vent zonal
@@ -159,4 +151,11 @@
       count(3)=nlevnc
       count(4)=1
+
+c mise a zero des tableaux 
+c ------------------------
+       unc(:,:,:)=0.
+       vnc(:,:,:)=0.
+       tnc(:,:,:)=0.
+       Qnc(:,:,:)=0.
 
 c  Vent zonal
Index: LMDZ4/branches/pre_V3/libf/dyn3d/startvar.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/startvar.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/startvar.F	(revision 701)
@@ -521,5 +521,8 @@
             CASE ('snow')
                   champ(:) = 0.0
-            CASE ('deltat')
+cIM "slab" ocean
+            CASE ('tslab')
+                   champ(:) = 0.0
+            CASE ('seaice')
                   champ(:) = 0.0
             CASE ('rugmer')
@@ -558,9 +561,11 @@
     !  LOCAL
     !
-      REAL :: lev(1), date, dt
+!ac     REAL :: lev(1), date, dt
+      REAL :: date, dt
+      REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini
+!ac
       INTEGER :: itau(1)
       INTEGER ::  llm_tmp, ttm_tmp
       INTEGER :: i, j
-      INTEGER :: iret
     !
       CHARACTER*25 title
@@ -579,12 +584,20 @@
      . ttm_tmp, fid_phys)
     !
-    !
       ALLOCATE (lat_phys(iml_phys,jml_phys))
       ALLOCATE (lon_phys(iml_phys,jml_phys))
+!ac
+      ALLOCATE (levphys_ini(llm_tmp))
+    !
+!      CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 
+!     . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp, 
+!     . itau, date, dt, fid_phys)
     !
       CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 
-     . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp, 
+     . llm_tmp, lon_phys, lat_phys, levphys_ini, ttm_tmp, 
      . itau, date, dt, fid_phys)
     !
+      DEALLOCATE (levphys_ini)
+!ac
+    !
     ! Allocate the space we will need to get the data out of this file
     !
@@ -593,7 +606,4 @@
     !   In case we have a file which is in degrees we do the transformation
     !
-      if (allocated(lon_rad)) DEALLOCATE(lon_rad, stat=iret)
-      if (allocated(lon_ini)) DEALLOCATE(lon_ini, stat=iret)
-
       ALLOCATE(lon_rad(iml_phys))
       ALLOCATE(lon_ini(iml_phys))
@@ -1044,5 +1054,5 @@
       REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
       REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
-       REAL, ALLOCATABLE :: varrr(:,:,:)
+!       REAL, ALLOCATABLE :: varrr(:,:,:)
       INTEGER, ALLOCATABLE :: lind(:)
     !
@@ -1052,5 +1062,5 @@
           ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
       ENDIF
-          ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn))
+!          ALLOCATE(varrr(iml_dyn, jml_dyn, llm_dyn))
     !
     !
@@ -1163,5 +1173,8 @@
 
       DEALLOCATE(lon_rad)
+      DEALLOCATE(lon_ini)
       DEALLOCATE(lat_rad)
+      DEALLOCATE(lat_ini)
+      DEALLOCATE(lev_dyn)
       DEALLOCATE(var_tmp2d)
       DEALLOCATE(var_tmp3d)
Index: LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j.F	(revision 701)
@@ -0,0 +1,137 @@
+c================================================================
+c================================================================
+      SUBROUTINE tetaleveli1j(ilon,ilev,lnew,pgcm,pres,Qgcm,Qpres)
+c================================================================
+c================================================================
+
+      IMPLICIT none
+
+#include "dimensions.h"
+#include "dimphy.h"
+
+c================================================================
+c
+c Interpoler des champs 3-D u, v et g du modele a un niveau de
+c pression donnee (pres)
+c
+c INPUT:  ilon ----- nombre de points
+c         ilev ----- nombre de couches
+c         lnew ----- true si on doit reinitialiser les poids
+c         pgcm ----- pressions modeles
+c         pres ----- pression vers laquelle on interpolle
+c         Qgcm ----- champ GCM
+c         Qpres ---- champ interpolle au niveau pres
+c
+c================================================================
+c
+c   arguments :
+c   -----------
+
+      INTEGER ilon, ilev
+      logical lnew
+
+      REAL pgcm(ilon,ilev)
+      REAL Qgcm(ilon,ilev)
+      real pres
+      REAL Qpres(ilon)
+
+c   local :
+c   -------
+
+cIM 211004
+c     INTEGER lt(klon), lb(klon)
+c     REAL ptop, pbot, aist(klon), aisb(klon)
+c
+#include "paramet.h"
+c
+      INTEGER lt(ip1jm), lb(ip1jm)
+      REAL ptop, pbot, aist(ip1jm), aisb(ip1jm)
+cMI 211004
+      save lt,lb,ptop,pbot,aist,aisb
+
+      INTEGER i, k
+c
+c     PRINT*,'tetalevel pres=',pres
+c=====================================================================
+      if (lnew) then
+c   on réinitialise les réindicages et les poids
+c=====================================================================
+
+
+c Chercher les 2 couches les plus proches du niveau a obtenir
+c
+c Eventuellement, faire l'extrapolation a partir des deux couches
+c les plus basses ou les deux couches les plus hautes:
+      DO 130 i = 1, ilon
+cIM      IF ( ABS(pres-pgcm(i,ilev) ) .LT.
+         IF ( ABS(pres-pgcm(i,ilev) ) .GT.
+     .        ABS(pres-pgcm(i,1)) ) THEN
+            lt(i) = ilev     ! 2
+            lb(i) = ilev-1   ! 1
+         ELSE
+            lt(i) = 2
+            lb(i) = 1
+         ENDIF
+cIM   PRINT*,'i, ABS(pres-pgcm),ABS(pres-pgcm)',
+cIM  .i, ABS(pres-pgcm(i,ilev)),ABS(pres-pgcm(i,1))
+  130 CONTINUE
+      DO 150 k = 1, ilev-1
+         DO 140 i = 1, ilon
+            pbot = pgcm(i,k)
+            ptop = pgcm(i,k+1)
+cIM         IF (ptop.LE.pres .AND. pbot.GE.pres) THEN
+            IF (ptop.GE.pres .AND. pbot.LE.pres) THEN
+               lt(i) = k+1
+               lb(i) = k
+            ENDIF
+  140    CONTINUE
+  150 CONTINUE
+c
+c Interpolation lineaire:
+c
+      DO i = 1, ilon
+c interpolation en logarithme de pression:
+c
+c ...   Modif . P. Le Van    ( 20/01/98) ....
+c       Modif Frédéric Hourdin (3/01/02)
+
+        IF(pgcm(i,lb(i)).EQ.0.OR.
+     $     pgcm(i,lt(i)).EQ.0.) THEN
+c
+        PRINT*,'i,lb,lt,2pgcm,pres',i,lb(i),
+     .  lt(i),pgcm(i,lb(i)),pgcm(i,lt(i)),pres
+c
+        ENDIF 
+c
+        aist(i) = LOG( pgcm(i,lb(i))/ pres )
+     .       / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i)) )
+        aisb(i) = LOG( pres / pgcm(i,lt(i)) )
+     .       / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i)))
+      enddo
+
+
+      endif ! lnew
+
+c======================================================================
+c    inteprollation
+c======================================================================
+
+      do i=1,ilon
+         Qpres(i)= Qgcm(i,lb(i))*aisb(i)+Qgcm(i,lt(i))*aist(i)
+cIM      PRINT*,'i,Qgcm,Qpres',i,Qgcm(i,lb(i)),aisb(i),
+cIM  $   Qgcm(i,lt(i)),aist(i),Qpres(i)
+      enddo
+c
+c Je mets les vents a zero quand je rencontre une montagne
+      do i = 1, ilon
+cIM      if (pgcm(i,1).LT.pres) THEN
+         if (pgcm(i,1).GT.pres) THEN
+c           Qpres(i)=1e33
+            Qpres(i)=1e+20
+cIM         PRINT*,'i,pgcm(i,1),pres =',i,pgcm(i,1),pres
+         endif
+      enddo
+
+c
+      RETURN
+      END
Index: LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j1.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j1.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/tetaleveli1j1.F	(revision 701)
@@ -0,0 +1,137 @@
+c================================================================
+c================================================================
+      SUBROUTINE tetaleveli1j1(ilon,ilev,lnew,pgcm,pres,Qgcm,Qpres)
+c================================================================
+c================================================================
+
+      IMPLICIT none
+
+#include "dimensions.h"
+#include "dimphy.h"
+
+c================================================================
+c
+c Interpoler des champs 3-D u, v et g du modele a un niveau de
+c pression donnee (pres)
+c
+c INPUT:  ilon ----- nombre de points
+c         ilev ----- nombre de couches
+c         lnew ----- true si on doit reinitialiser les poids
+c         pgcm ----- pressions modeles
+c         pres ----- pression vers laquelle on interpolle
+c         Qgcm ----- champ GCM
+c         Qpres ---- champ interpolle au niveau pres
+c
+c================================================================
+c
+c   arguments :
+c   -----------
+
+      INTEGER ilon, ilev
+      logical lnew
+
+      REAL pgcm(ilon,ilev)
+      REAL Qgcm(ilon,ilev)
+      real pres
+      REAL Qpres(ilon)
+
+c   local :
+c   -------
+
+cIM 211004
+c     INTEGER lt(klon), lb(klon)
+c     REAL ptop, pbot, aist(klon), aisb(klon)
+c
+#include "paramet.h"
+c
+      INTEGER lt(ip1jmp1), lb(ip1jmp1)
+      REAL ptop, pbot, aist(ip1jmp1), aisb(ip1jmp1)
+cMI 211004
+      save lt,lb,ptop,pbot,aist,aisb
+
+      INTEGER i, k
+c
+c     PRINT*,'tetalevel pres=',pres
+c=====================================================================
+      if (lnew) then
+c   on réinitialise les réindicages et les poids
+c=====================================================================
+
+
+c Chercher les 2 couches les plus proches du niveau a obtenir
+c
+c Eventuellement, faire l'extrapolation a partir des deux couches
+c les plus basses ou les deux couches les plus hautes:
+      DO 130 i = 1, ilon
+cIM      IF ( ABS(pres-pgcm(i,ilev) ) .LT.
+         IF ( ABS(pres-pgcm(i,ilev) ) .GT.
+     .        ABS(pres-pgcm(i,1)) ) THEN
+            lt(i) = ilev     ! 2
+            lb(i) = ilev-1   ! 1
+         ELSE
+            lt(i) = 2
+            lb(i) = 1
+         ENDIF
+cIM   PRINT*,'i, ABS(pres-pgcm),ABS(pres-pgcm)',
+cIM  .i, ABS(pres-pgcm(i,ilev)),ABS(pres-pgcm(i,1))
+  130 CONTINUE
+      DO 150 k = 1, ilev-1
+         DO 140 i = 1, ilon
+            pbot = pgcm(i,k)
+            ptop = pgcm(i,k+1)
+cIM         IF (ptop.LE.pres .AND. pbot.GE.pres) THEN
+            IF (ptop.GE.pres .AND. pbot.LE.pres) THEN
+               lt(i) = k+1
+               lb(i) = k
+            ENDIF
+  140    CONTINUE
+  150 CONTINUE
+c
+c Interpolation lineaire:
+c
+      DO i = 1, ilon
+c interpolation en logarithme de pression:
+c
+c ...   Modif . P. Le Van    ( 20/01/98) ....
+c       Modif Frédéric Hourdin (3/01/02)
+
+        IF(pgcm(i,lb(i)).EQ.0.OR.
+     $     pgcm(i,lt(i)).EQ.0.) THEN
+c
+        PRINT*,'i,lb,lt,2pgcm,pres',i,lb(i),
+     .  lt(i),pgcm(i,lb(i)),pgcm(i,lt(i)),pres
+c
+        ENDIF 
+c
+        aist(i) = LOG( pgcm(i,lb(i))/ pres )
+     .       / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i)) )
+        aisb(i) = LOG( pres / pgcm(i,lt(i)) )
+     .       / LOG( pgcm(i,lb(i))/ pgcm(i,lt(i)))
+      enddo
+
+
+      endif ! lnew
+
+c======================================================================
+c    inteprollation
+c======================================================================
+
+      do i=1,ilon
+         Qpres(i)= Qgcm(i,lb(i))*aisb(i)+Qgcm(i,lt(i))*aist(i)
+cIM      PRINT*,'i,Qgcm,Qpres',i,Qgcm(i,lb(i)),aisb(i),
+cIM  $   Qgcm(i,lt(i)),aist(i),Qpres(i)
+      enddo
+c
+c Je mets les vents a zero quand je rencontre une montagne
+      do i = 1, ilon
+cIM      if (pgcm(i,1).LT.pres) THEN
+         if (pgcm(i,1).GT.pres) THEN
+c           Qpres(i)=1e33
+            Qpres(i)=1e+20
+cIM         PRINT*,'i,pgcm(i,1),pres =',i,pgcm(i,1),pres
+         endif
+      enddo
+
+c
+      RETURN
+      END
Index: LMDZ4/branches/pre_V3/libf/dyn3d/tourabs.F
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/tourabs.F	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/tourabs.F	(revision 701)
@@ -0,0 +1,98 @@
+      SUBROUTINE tourabs ( ntetaSTD,vcov, ucov, vorabs )
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Modif:  I. Musat (28/10/04)
+c   -------
+c   adaptation du code tourpot.F pour le calcul de la vorticite absolue
+c   cf. P. Le Van
+c
+c   Objet: 
+c   ------
+c
+c    *******************************************************************
+c    .............  calcul de la vorticite absolue     .................
+c    *******************************************************************
+c
+c     ntetaSTD, vcov,ucov      sont des argum. d'entree pour le s-pg .
+c             vorabs            est  un argum.de sortie pour le s-pg .
+c
+c=======================================================================
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+#include "logic.h"
+#include "comconst.h"
+c
+      INTEGER ntetaSTD
+      REAL vcov( ip1jm,ntetaSTD ), ucov( ip1jmp1,ntetaSTD )
+      REAL vorabs( ip1jm,ntetaSTD )
+c
+c variables locales
+      INTEGER l, ij, i, j
+      REAL  rot( ip1jm,ntetaSTD )
+
+
+
+c  ... vorabs = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) ..
+
+
+
+c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
+
+      DO 5 l = 1,ntetaSTD
+
+      DO 2 i = 1, iip1
+      DO 2 j = 1, jjm
+c
+       ij=i+(j-1)*iip1
+       IF(ij.LE.ip1jm - 1) THEN
+c
+        IF(cv(ij).EQ.0..OR.cv(ij+1).EQ.0..OR.
+     $     cu(ij).EQ.0..OR.cu(ij+iip1).EQ.0.) THEN
+         rot( ij,l ) = 0.
+         continue
+        ELSE
+         rot( ij,l ) = (vcov(ij+1,l)/cv(ij+1)-vcov(ij,l)/cv(ij))/
+     $                 (2.*pi*RAD*cos(rlatv(j)))*float(iim)
+     $                +(ucov(ij+iip1,l)/cu(ij+iip1)-ucov(ij,l)/cu(ij))/
+     $                 (pi*RAD)*(float(jjm)-1.)
+c
+        ENDIF
+       ENDIF !(ij.LE.ip1jm - 1) THEN
+   2  CONTINUE
+
+c    ....  correction pour  rot( iip1,j,l )  .....
+c    ....     rot(iip1,j,l) = rot(1,j,l)    .....
+
+CDIR$ IVDEP
+
+      DO 3 ij = iip1, ip1jm, iip1
+      rot( ij,l ) = rot( ij -iim, l )
+   3  CONTINUE
+
+   5  CONTINUE
+
+
+      CALL  filtreg( rot, jjm, ntetaSTD, 2, 1, .FALSE., 1 )
+
+
+      DO 10 l = 1, ntetaSTD
+
+      DO 6 ij = 1, ip1jm - 1
+      vorabs( ij,l ) = ( rot(ij,l) + fext(ij)*unsairez(ij) )
+   6  CONTINUE
+
+c    ..... correction pour  vorabs( iip1,j,l)  .....
+c    ....   vorabs(iip1,j,l)= vorabs(1,j,l) ....
+CDIR$ IVDEP
+      DO 8 ij = iip1, ip1jm, iip1
+      vorabs( ij,l ) = vorabs( ij -iim,l )
+   8  CONTINUE
+
+  10  CONTINUE
+
+      RETURN
+      END
Index: LMDZ4/branches/pre_V3/libf/dyn3d/vlsplt.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/vlsplt.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/vlsplt.F	(revision 701)
@@ -352,5 +352,5 @@
       ENDDO
 
-      IF(n0.gt.1) THEN
+      IF(n0.gt.0) THEN
       PRINT*,'Nombre de points pour lesquels on advect plus que le'
      &       ,'contenu de la maille : ',n0
Index: LMDZ4/branches/pre_V3/libf/dyn3d/vlspltqs.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/vlspltqs.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/vlspltqs.F	(revision 701)
@@ -389,5 +389,5 @@
       ENDDO
 
-      IF(n0.gt.1) THEN
+      IF(n0.gt.0) THEN
 ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
 ccc     &       ,'contenu de la maille : ',n0
Index: LMDZ4/branches/pre_V3/libf/dyn3d/wrgrads.F
===================================================================
--- LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/wrgrads.F	(revision 525)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/wrgrads.F	(revision 701)
@@ -29,5 +29,5 @@
       writectl=.false.
 
-c     print*,if,iid(if),jid(if),ifd(if),jfd(if)
+      print*,if,iid(if),jid(if),ifd(if),jfd(if)
       iii=iid(if)
       iji=jid(if)
@@ -38,6 +38,6 @@
       lm=lmd(if)
 
-c      print*,'im,jm,lm,name,firsttime(if)'
-c      print*,im,jm,lm,name,firsttime(if)
+      print*,'im,jm,lm,name,firsttime(if)'
+      print*,im,jm,lm,name,firsttime(if)
 
       if(firsttime(if)) then
@@ -81,6 +81,6 @@
       endif
 
-c     print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
-c     print*,ivar(if),nvar(if),var(ivar(if),if),writectl
+      print*,'ivar(if),nvar(if),var(ivar(if),if),writectl'
+      print*,ivar(if),nvar(if),var(ivar(if),if),writectl
       do l=1,nl
          irec(if)=irec(if)+1
Index: LMDZ4/branches/pre_V3/libf/dyn3d/write_paramLMDZ_dyn.h
===================================================================
--- LMDZ4/branches/pre_V3/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 701)
+++ LMDZ4/branches/pre_V3/libf/dyn3d/write_paramLMDZ_dyn.h	(revision 701)
@@ -0,0 +1,240 @@
+c
+      ndex2d = 0
+      itau_dyn=itau_dyn+itau
+      print*,'write_paramLMDZ_dyn :itau_dyn=',itau_dyn
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(prt_level) 
+      CALL histwrite(nid_ctesGCM, "prt_level", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(dayref)
+      CALL histwrite(nid_ctesGCM, "dayref", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(anneeref)
+      CALL histwrite(nid_ctesGCM, "anneeref", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(raz_date)
+      CALL histwrite(nid_ctesGCM, "raz_date", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nday)
+      CALL histwrite(nid_ctesGCM, "nday", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(day_step)
+      CALL histwrite(nid_ctesGCM, "day_step", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iperiod)
+      CALL histwrite(nid_ctesGCM, "iperiod", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iapp_tracvl)
+      CALL histwrite(nid_ctesGCM, "iapp_tracvl", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iconser)
+      CALL histwrite(nid_ctesGCM, "iconser", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iecri)
+      CALL histwrite(nid_ctesGCM, "iecri", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=periodav
+      CALL histwrite(nid_ctesGCM, "periodav", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(idissip)
+      CALL histwrite(nid_ctesGCM, "idissip", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(lstardis) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "lstardis", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nitergdiv)
+      CALL histwrite(nid_ctesGCM, "nitergdiv", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nitergrot)
+      CALL histwrite(nid_ctesGCM, "nitergrot", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(niterh) 
+      CALL histwrite(nid_ctesGCM, "niterh", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=tetagdiv
+      CALL histwrite(nid_ctesGCM, "tetagdiv", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=tetagrot
+      CALL histwrite(nid_ctesGCM, "tetagrot", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=tetatemp
+      CALL histwrite(nid_ctesGCM, "tetatemp", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=coefdis
+      CALL histwrite(nid_ctesGCM, "coefdis", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(purmats) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "purmats", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(ok_guide) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "ok_guide", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(true_calendar) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "true_calendar", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iflag_phys)
+      CALL histwrite(nid_ctesGCM, "iflag_phys", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iphysiq)
+      CALL histwrite(nid_ctesGCM, "iphysiq", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(ecritphy)
+      CALL histwrite(nid_ctesGCM, "ecritphy", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(cycle_diurne) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "cycle_diurne", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(soil_model) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "soil_model", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(new_oliq) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "new_oliq", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(ok_orodr) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "ok_orodr", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(ok_orolf) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "ok_orolf", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(ok_limitvrai) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "ok_limitvrai", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(nbapp_rad)
+      CALL histwrite(nid_ctesGCM, "nbapp_rad", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(iflag_con)
+      CALL histwrite(nid_ctesGCM, "iflag_con", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=clon
+      CALL histwrite(nid_ctesGCM, "clon", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=clat
+      CALL histwrite(nid_ctesGCM, "clat", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=grossismx
+      CALL histwrite(nid_ctesGCM, "grossismx", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=grossismy
+      CALL histwrite(nid_ctesGCM, "grossismy", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(fxyhypb) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "fxyhypb", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=dzoomx
+      CALL histwrite(nid_ctesGCM, "dzoomx", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=dzoomy
+      CALL histwrite(nid_ctesGCM, "dzoomy", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=taux
+      CALL histwrite(nid_ctesGCM, "taux", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=tauy
+      CALL histwrite(nid_ctesGCM, "tauy", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      IF(ysinus) THEN
+       zx_tmp_2d(1:iip1,1:jjp1)=1.
+      ELSE
+       zx_tmp_2d(1:iip1,1:jjp1)=0.
+      ENDIF
+      CALL histwrite(nid_ctesGCM, "ysinus", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+      zx_tmp_2d(1:iip1,1:jjp1)=FLOAT(ip_ebil_dyn)
+      CALL histwrite(nid_ctesGCM, "ip_ebil_dyn", itau_dyn,
+     .               zx_tmp_2d,iip1*jjp1,ndex2d)
+c
+c=================================================================
+c
+      if (ok_sync) then
+        call histsync(nid_ctesGCM)
+      endif
+c
+c=================================================================
