Index: /LMDZ4/trunk/libf/dyn3dpar/caldyn0.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/caldyn0.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/caldyn0.F	(revision 1000)
@@ -0,0 +1,89 @@
+!
+! $Header$
+!
+      SUBROUTINE caldyn0
+     $ (itau,ucov,vcov,teta,ps,masse,pk,phis ,
+     $  phi,w,pbaru,pbarv,time )
+
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c  Auteur :  P. Le Van
+c
+c   Objet:
+c   ------
+c
+c   Calcul des tendances dynamiques.
+c
+c Modif 04/93 F.Forget
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   0. Declarations:
+c   ----------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comvert.h"
+#include "comgeom.h"
+
+c   Arguments:
+c   ----------
+
+      INTEGER itau
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
+      REAL ps(ip1jmp1),phis(ip1jmp1)
+      REAL pk(iip1,jjp1,llm)
+      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
+      REAL phi(ip1jmp1,llm),masse(ip1jmp1,llm)
+      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL time
+
+c   Local:
+c   ------
+
+      REAL ang(ip1jmp1,llm),p(ip1jmp1,llmp1)
+      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),psexbarxy(ip1jm)
+      REAL vorpot(ip1jm,llm)
+      REAL w(ip1jmp1,llm),ecin(ip1jmp1,llm),convm(ip1jmp1,llm)
+      REAL bern(ip1jmp1,llm)
+      REAL massebxy(ip1jm,llm), dp(ip1jmp1)
+    
+
+      INTEGER   ij,l
+
+c-----------------------------------------------------------------------
+c   Calcul des tendances dynamiques:
+c   --------------------------------
+
+      CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
+      CALL pression ( ip1jmp1, ap      , bp   ,  ps  , p            )
+      CALL psextbar (   ps   , psexbarxy                            )
+      CALL massdair (    p   , masse                                )
+      CALL massbar  (   masse, massebx , masseby                    )
+      CALL massbarxy(   masse, massebxy                             )
+      CALL flumass  ( massebx, masseby , vcont, ucont ,pbaru, pbarv )
+      CALL convmas  (   pbaru, pbarv   , convm                      )
+
+      DO ij =1, ip1jmp1
+         dp( ij ) = convm( ij,1 ) / airesurg( ij )
+      ENDDO
+
+      CALL vitvert ( convm  , w                                  )
+      CALL tourpot ( vcov   , ucov  , massebxy  , vorpot         )
+      CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
+      CALL bernoui ( ip1jmp1, llm   , phi       , ecin   , bern  )
+
+      DO l=1,llm
+         DO ij=1,ip1jmp1
+            ang(ij,l) = ucov(ij,l) + constang(ij)
+         ENDDO
+      ENDDO
+
+        CALL sortvarc0
+     $ ( itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time,vcov )
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/calfis_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/calfis_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/calfis_p.F	(revision 1000)
@@ -34,5 +34,5 @@
       USE dimphy
       USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root 
-      USE parallel, ONLY : omp_chunk
+      USE parallel, ONLY : omp_chunk, using_mpi
       USE mod_interface_dyn_phys
       USE Write_Field
@@ -107,6 +107,7 @@
 #include "comgeom2.h"
 #include "control.h"
+#ifdef CPP_MPI
       include 'mpif.h'
-
+#endif
 c    Arguments :
 c    -----------
@@ -212,5 +213,9 @@
       REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv
       INTEGER :: ierr
+#ifdef CPP_MPI
       INTEGER,dimension(MPI_STATUS_SIZE,4) :: Status
+#else
+      INTEGER,dimension(1,4) :: Status
+#endif
       INTEGER, dimension(4) :: Req
       REAL,ALLOCATABLE,SAVE:: zdufi2(:,:),zdvfi2(:,:)
@@ -792,5 +797,7 @@
       call stop_timer(timer_physic)
 c$OMP END MASTER
-      
+
+      IF (using_mpi) THEN
+            
       if (MPI_rank>0) then
 
@@ -803,4 +810,5 @@
 
 c$OMP BARRIER
+#ifdef CPP_MPI 
 c$OMP MASTER
 !$OMP CRITICAL (MPI)
@@ -811,4 +819,5 @@
 !$OMP END CRITICAL (MPI)
 c$OMP END MASTER
+#endif
 c$OMP BARRIER
      
@@ -817,4 +826,5 @@
       if (MPI_rank<MPI_Size-1) then
 c$OMP BARRIER
+#ifdef CPP_MPI 
 c$OMP MASTER      
 !$OMP CRITICAL (MPI)
@@ -825,8 +835,11 @@
 !$OMP END CRITICAL (MPI)
 c$OMP END MASTER
-c$OMP BARRIER     
+#endif
       endif
 
 c$OMP BARRIER
+
+
+#ifdef CPP_MPI 
 c$OMP MASTER    
 !$OMP CRITICAL (MPI)
@@ -840,6 +853,11 @@
 !$OMP END CRITICAL (MPI)
 c$OMP END MASTER
+#endif
+
 c$OMP BARRIER     
 
+      ENDIF ! using_mpi
+      
+      
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       DO l=1,llm
Index: /LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/conf_gcm.F	(revision 1000)
@@ -712,4 +712,12 @@
       omp_chunk=1
       CALL getin('omp_chunk',omp_chunk)
+
+!Config key = ok_strato
+!Config  Desc = activation de la version strato
+!Config  Def  = .FALSE.
+!Config  Help = active la version stratosphérique de LMDZ de F. Lott
+
+      ok_strato=.FALSE.
+      CALL getin('ok_strato',ok_strato)
 
       write(lunout,*)' #########################################'
@@ -748,4 +756,5 @@
       write(lunout,*)' use_mpi_alloc = ', use_mpi_alloc
       write(lunout,*)' omp_chunk = ', omp_chunk
+      write(lunout,*)' ok_strato = ', ok_strato
 c
       RETURN
Index: /LMDZ4/trunk/libf/dyn3dpar/convflu.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/convflu.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/convflu.F	(revision 1000)
@@ -0,0 +1,62 @@
+!
+! $Header$
+!
+      SUBROUTINE convflu( xflu,yflu,nbniv,convfl )
+c
+c  P. Le Van
+c
+c
+c    *******************************************************************
+c  ... calcule la (convergence horiz. * aire locale)du flux ayant pour
+c      composantes xflu et yflu ,variables extensives .  ......
+c    *******************************************************************
+c      xflu , yflu et nbniv sont des arguments d'entree pour le s-pg ..
+c      convfl                est  un argument de sortie pour le s-pg .
+c
+c     njxflu  est le nombre de lignes de latitude de xflu, 
+c     ( = jjm ou jjp1 )
+c     nbniv   est le nombre de niveaux vert. de  xflu et de yflu .
+c
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+      REAL       xflu,yflu,convfl,convpn,convps
+      INTEGER    l,ij,nbniv
+      DIMENSION  xflu( ip1jmp1,nbniv ),yflu( ip1jm,nbniv ) ,
+     *         convfl( ip1jmp1,nbniv )
+c
+      REAL       SSUM
+c
+c
+#include "comgeom.h"
+c
+      DO 5 l = 1,nbniv
+c
+      DO 2  ij = iip2, ip1jm - 1
+      convfl( ij + 1,l ) =  xflu(   ij,l ) - xflu( ij +  1,l )   +
+     *                      yflu(ij +1,l ) - yflu( ij -iim,l )
+   2  CONTINUE
+c
+c
+
+c     ....  correction pour  convfl( 1,j,l)  ......
+c     ....   convfl(1,j,l)= convfl(iip1,j,l) ...
+c
+CDIR$ IVDEP
+      DO 3 ij = iip2,ip1jm,iip1
+      convfl( ij,l ) = convfl( ij + iim,l )
+   3  CONTINUE
+c
+c     ......  calcul aux poles  .......
+c
+      convpn =   SSUM( iim, yflu(     1    ,l ),  1 )
+      convps = - SSUM( iim, yflu( ip1jm-iim,l ),  1 )
+      DO 4 ij = 1,iip1
+      convfl(     ij   ,l ) = convpn * aire(   ij     ) / apoln
+      convfl( ij+ ip1jm,l ) = convps * aire( ij+ ip1jm) / apols
+   4  CONTINUE
+c
+   5  CONTINUE
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/convmas.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/convmas.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/convmas.F	(revision 1000)
@@ -0,0 +1,63 @@
+!
+! $Header$
+!
+      SUBROUTINE convmas (pbaru, pbarv, convm )
+c
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteurs:  P. Le Van , F. Hourdin  .
+c   -------
+c
+c   Objet:
+c   ------
+c
+c   ********************************************************************
+c   .... calcul de la convergence du flux de masse aux niveaux p ...
+c   ********************************************************************
+c
+c
+c     pbaru  et  pbarv  sont des arguments d'entree pour le s-pg  ....
+c      .....  convm      est  un argument de sortie pour le s-pg  ....
+c
+c    le calcul se fait de haut en bas, 
+c    la convergence de masse au niveau p(llm+1) est egale a 0. et
+c    n'est pas stockee dans le tableau convm .
+c
+c
+c=======================================================================
+c
+c   Declarations:
+c   -------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+#include "logic.h"
+
+      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm ),convm(  ip1jmp1,llm )
+      INTEGER   l,ij
+
+
+c-----------------------------------------------------------------------
+c    ....  calcul de - (d(pbaru)/dx + d(pbarv)/dy ) ......
+
+      CALL  convflu( pbaru, pbarv, llm, convm )
+
+c-----------------------------------------------------------------------
+c   filtrage:
+c   ---------
+
+       CALL filtreg( convm, jjp1, llm, 2, 2, .true., 1 )
+
+c    integration de la convergence de masse de haut  en bas ......
+
+      DO      l      = llmm1, 1, -1
+        DO    ij     = 1, ip1jmp1
+         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+        ENDDO
+      ENDDO
+c
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/covnat_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/covnat_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/covnat_p.F	(revision 1000)
@@ -66,5 +66,5 @@
      
       DO l = 1,klevel
-         DO ij = 1,ip1jm
+         DO ij = ijb,ije
             vnat( ij,l ) = vcov( ij,l ) / cv(ij)
          ENDDO
Index: /LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/create_etat0_limit.F	(revision 1000)
@@ -0,0 +1,87 @@
+!
+! $Header$
+!
+       PROGRAM create_etat0_limit
+       USE dimphy
+       USE comgeomphy
+       USE mod_phys_lmdz_para
+       USE mod_const_mpi
+       USE phys_state_var_mod      
+       IMPLICIT NONE
+c
+c
+c     Programme d'appel a etat0, creation des etats initiaux et limit_netcdf
+c   
+c
+c     interbar = .T . si appel a  interpol. barycentrique inter_barxy
+c
+c     extrap   = .T . si on fait une extrapolation de donnees , comme pour
+c       les  SST  lorsque  le fichier ne contient pas uniquement  des points 
+c     oceaniques .
+c
+c     oldice   = .T. si l'on veut garder les anciennes glaces , obtenues
+c     par  grille_m  ( grid_atob ) .
+c
+c     on cree le masque dans etat0 que l'on passe ensuite dans limit pour 
+c     garder les coherences
+
+      LOGICAL interbar, extrap , oldice
+      PARAMETER ( interbar = .true. , extrap = .FALSE. , oldice=.false.)
+#include "dimensions.h"
+#include "paramet.h"
+#include "indicesol.h"
+#include "advtrac.h"
+#include  "control.h"
+#include "clesphys.h"
+      REAL :: masque(iip1,jjp1)
+!      REAL :: pctsrf(iim*(jjm-1)+2, nbsrf)
+
+c initialisation traceurs
+      hadv_flg(:) = 0.
+      vadv_flg(:) = 0.
+      conv_flg(:) = 0.
+      pbl_flg(:)  = 0. 
+      tracnam(:)  = '        ' 
+      nprath = 1
+      nbtrac = 0
+      mmt_adj(:,:,:,:) = 1
+
+      IF (config_inca /= 'none') THEN
+#ifdef INCA
+         call init_const_lmdz(nbtrac,anneeref,dayref,iphysiq,day_step,nday)
+#endif
+         print *, 'nbtrac =' , nbtrac 
+      END IF
+
+      CALL init_mpi
+
+
+      CALL Init_Phys_lmdz(iim,jjp1,llm,nqmx-2,1,(/(jjm-1)*iim+2/))
+      PRINT *,'---> klon=',klon
+
+      IF (mpi_size>1 .OR. omp_size>1) THEN
+        CALL abort_gcm('create_etat0_limit','In parallel mode, 
+     &                 create_etat0_limit must be called only 
+     &                 for 1 process and 1 task')
+      ENDIF
+      CALL phys_state_var_init
+      call InitComgeomphy
+      
+
+      WRITE(6,*) '  *********************  '
+      WRITE(6,*) ' interbar = ',interbar
+      CALL etat0_netcdf ( interbar, masque )
+c
+      WRITE(6,1)
+      WRITE(6,*) '  *********************  '
+      WRITE(6,*) '  ***  Limit_netcdf ***  '
+      WRITE(6,*) '  *********************  '
+      WRITE(6,1)
+      
+c     
+      CALL  limit_netcdf ( interbar, extrap , oldice, masque)
+
+1     FORMAT(//)
+
+      STOP
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/disvert.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/disvert.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/disvert.F	(revision 1000)
@@ -11,4 +11,5 @@
 #include "paramet.h"
 #include "iniprint.h"
+#include "logic.h"
 c
 c=======================================================================
@@ -99,5 +100,12 @@
       DO l = 1, llm
          x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
-         dsig(l) = 1.0 + 7.0 * SIN(x)**2
+
+         IF (ok_strato) THEN
+           dsig(l) =(1.0 + 7.0 * SIN(x)**2)
+     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
+         ELSE
+           dsig(l) = 1.0 + 7.0 * SIN(x)**2
+         ENDIF
+
          snorm = snorm + dsig(l)
       ENDDO
Index: /LMDZ4/trunk/libf/dyn3dpar/disvert0.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/disvert0.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/disvert0.F	(revision 1000)
@@ -13,4 +13,5 @@
 #include "paramet.h"
 #include "iniprint.h"
+#include "logic.h"
 c
 c=======================================================================
@@ -104,5 +105,12 @@
       DO l = 1, llm
          x = 2.*asin(1.) * (FLOAT(l)-0.5) / float(llm+1)
-         dsig(l) = 1.0 + 7.0 * SIN(x)**2
+
+         IF (ok_strato) THEN
+           dsig(l) =(1.0 + 7.0 * SIN(x)**2)
+     &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2        
+         ELSE
+           dsig(l) = 1.0 + 7.0 * SIN(x)**2
+         ENDIF
+
          snorm = snorm + dsig(l)
       ENDDO
Index: /LMDZ4/trunk/libf/dyn3dpar/dynredem.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/dynredem.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/dynredem.F	(revision 1000)
@@ -0,0 +1,734 @@
+!
+! $Header$
+!
+c
+      SUBROUTINE dynredem0(fichnom,iday_end,phis,nq)
+      USE IOIPSL
+      IMPLICIT NONE
+c=======================================================================
+c Ecriture du fichier de redemarrage sous format NetCDF (initialisation)
+c=======================================================================
+c   Declarations:
+c   -------------
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comvert.h"
+#include "comgeom.h"
+#include "temps.h"
+#include "ener.h"
+#include "logic.h"
+#include "netcdf.inc"
+#include "description.h"
+#include "serre.h"
+#include "advtrac.h"
+
+c   Arguments:
+c   ----------
+      INTEGER iday_end
+      REAL phis(ip1jmp1)
+      CHARACTER*(*) fichnom
+      INTEGER nq
+
+c   Local:
+c   ------
+      INTEGER iq,l
+      INTEGER length
+      PARAMETER (length = 100)
+      REAL tab_cntrl(length) ! tableau des parametres du run
+      INTEGER ierr
+      character*20 modname
+      character*80 abort_message
+
+c   Variables locales pour NetCDF:
+c
+      INTEGER dims2(2), dims3(3), dims4(4)
+      INTEGER idim_index
+      INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_rlatv
+      INTEGER idim_s, idim_sig
+      INTEGER idim_tim
+      INTEGER nid,nvarid
+
+      REAL zan0,zjulian,hours
+      INTEGER yyears0,jjour0, mmois0
+      character*30 unites
+
+
+c-----------------------------------------------------------------------
+      modname='dynredem'
+
+      call ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
+      call ju2ymds(zjulian, yyears0, mmois0, jjour0, hours)
+        
+
+      DO l=1,length
+       tab_cntrl(l) = 0.
+      ENDDO
+       tab_cntrl(1)  = FLOAT(iim)
+       tab_cntrl(2)  = FLOAT(jjm)
+       tab_cntrl(3)  = FLOAT(llm)
+       tab_cntrl(4)  = FLOAT(day_ref)
+       tab_cntrl(5)  = FLOAT(annee_ref)
+       tab_cntrl(6)  = rad
+       tab_cntrl(7)  = omeg
+       tab_cntrl(8)  = g
+       tab_cntrl(9)  = cpp
+       tab_cntrl(10) = kappa
+       tab_cntrl(11) = daysec
+       tab_cntrl(12) = dtvr
+       tab_cntrl(13) = etot0
+       tab_cntrl(14) = ptot0
+       tab_cntrl(15) = ztot0
+       tab_cntrl(16) = stot0
+       tab_cntrl(17) = ang0
+       tab_cntrl(18) = pa
+       tab_cntrl(19) = preff
+c
+c    .....    parametres  pour le zoom      ......   
+
+       tab_cntrl(20)  = clon
+       tab_cntrl(21)  = clat
+       tab_cntrl(22)  = grossismx
+       tab_cntrl(23)  = grossismy
+c
+      IF ( fxyhypb )   THEN
+       tab_cntrl(24) = 1.
+       tab_cntrl(25) = dzoomx
+       tab_cntrl(26) = dzoomy
+       tab_cntrl(27) = 0.
+       tab_cntrl(28) = taux
+       tab_cntrl(29) = tauy
+      ELSE
+       tab_cntrl(24) = 0.
+       tab_cntrl(25) = dzoomx
+       tab_cntrl(26) = dzoomy
+       tab_cntrl(27) = 0.
+       tab_cntrl(28) = 0.
+       tab_cntrl(29) = 0.
+       IF( ysinus )  tab_cntrl(27) = 1.
+      ENDIF
+
+       tab_cntrl(30) = FLOAT(iday_end)
+       tab_cntrl(31) = FLOAT(itau_dyn + itaufin)
+c
+c    .........................................................
+c
+c Creation du fichier:
+c
+      ierr = NF_CREATE(fichnom, NF_CLOBBER, nid)
+      IF (ierr.NE.NF_NOERR) THEN
+         WRITE(6,*)" Pb d ouverture du fichier "//fichnom
+         WRITE(6,*)' ierr = ', ierr
+         CALL ABORT
+      ENDIF
+c
+c Preciser quelques attributs globaux:
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 27,
+     .                       "Fichier demmarage dynamique")
+c
+c Definir les dimensions du fichiers:
+c
+      ierr = NF_DEF_DIM (nid, "index", length, idim_index)
+      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid, "rlatu", jjp1, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "rlonv", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+      ierr = NF_DEF_DIM (nid, "sigs", llm, idim_s)
+      ierr = NF_DEF_DIM (nid, "sig", llmp1, idim_sig)
+      ierr = NF_DEF_DIM (nid, "temps", NF_UNLIMITED, idim_tim)
+c
+      ierr = NF_ENDDEF(nid) ! sortir du mode de definition
+c
+c Definir et enregistrer certains champs invariants:
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsigs)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsigs)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,nivsig)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,nivsig)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ap)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ap)
+#endif
+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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bp)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,bp)
+#endif
+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
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,presnivs)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,presnivs)
+#endif
+c
+c Coefficients de passage cov. <-> contra. <--> naturel
+c
+      ierr = NF_REDEF (nid)
+      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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+#endif
+c
+      ierr = NF_REDEF (nid)
+      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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+#endif
+c
+c Aire de chaque maille:
+c
+      ierr = NF_REDEF (nid)
+      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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+#endif
+c
+c Geopentiel au sol:
+c
+      ierr = NF_REDEF (nid)
+      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")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phis)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,phis)
+#endif
+c
+c Definir les variables pour pouvoir les enregistrer plus tard:
+c
+      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")
+      write(unites,200)yyears0,mmois0,jjour0
+200   format('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", 30,
+     .                         unites)
+
+c
+      dims4(1) = idim_rlonu
+      dims4(2) = idim_rlatu
+      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")
+c
+      dims4(1) = idim_rlonv
+      dims4(2) = idim_rlatv
+      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")
+c
+      dims4(1) = idim_rlonv
+      dims4(2) = idim_rlatu
+      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")
+c
+      dims4(1) = idim_rlonv
+      dims4(2) = idim_rlatu
+      dims4(3) = idim_s
+      dims4(4) = idim_tim
+      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
+      ENDIF
+c
+      dims4(1) = idim_rlonv
+      dims4(2) = idim_rlatu
+      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 ?")
+c
+      dims3(1) = idim_rlonv
+      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")
+c
+      ierr = NF_ENDDEF(nid) ! sortir du mode de definition
+      ierr = NF_CLOSE(nid) ! fermer le fichier
+
+      PRINT*,'iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
+      PRINT*,'rad,omeg,g,cpp,kappa',
+     ,        rad,omeg,g,cpp,kappa
+
+      RETURN
+      END
+      SUBROUTINE dynredem1(fichnom,time,
+     .                     vcov,ucov,teta,q,nq,masse,ps)
+      IMPLICIT NONE
+c=================================================================
+c  Ecriture du fichier de redemarrage sous format NetCDF
+c=================================================================
+#include "dimensions.h"
+#include "paramet.h"
+#include "description.h"
+#include "netcdf.inc"
+#include "comvert.h"
+#include "comgeom.h"
+#include "advtrac.h"
+#include "temps.h"
+#include "control.h"
+
+      INTEGER nq, l
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) 
+      REAL teta(ip1jmp1,llm)                   
+      REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
+      REAL q(ip1jmp1,llm,nq)
+      CHARACTER*(*) fichnom
+     
+      REAL time
+      INTEGER nid, nvarid, nid_trac, nvarid_trac
+      REAL trac_tmp(ip1jmp1,llm)      
+      INTEGER ierr, ierr_file 
+      INTEGER iq
+      INTEGER length
+      PARAMETER (length = 100)
+      REAL tab_cntrl(length) ! tableau des parametres du run
+      character*20 modname
+      character*80 abort_message
+c
+      INTEGER nb
+      SAVE nb
+      DATA nb / 0 /
+
+      modname = 'dynredem1'
+      ierr = NF_OPEN(fichnom, NF_WRITE, nid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Pb. d ouverture "//fichnom
+         CALL abort
+      ENDIF
+
+c  Ecriture/extension de la coordonnee temps
+
+      nb = nb + 1
+      ierr = NF_INQ_VARID(nid, "temps", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         print *, NF_STRERROR(ierr)
+         abort_message='Variable temps n est pas definie'
+         CALL abort_gcm(modname,abort_message,ierr)
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR1_DOUBLE (nid,nvarid,nb,time)
+#else
+      ierr = NF_PUT_VAR1_REAL (nid,nvarid,nb,time)
+#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
+c
+      ierr = NF_INQ_VARID(nid, "ucov", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Variable ucov n est pas definie"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ucov)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ucov)
+#endif
+
+      ierr = NF_INQ_VARID(nid, "vcov", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Variable vcov n est pas definie"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,vcov)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,vcov)
+#endif
+
+      ierr = NF_INQ_VARID(nid, "teta", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Variable teta n est pas definie"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,teta)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,teta)
+#endif
+
+      IF (config_inca /= 'none') THEN
+! Ajout Anne pour lecture valeurs traceurs dans un fichier start_trac.nc
+         ierr_file = NF_OPEN ("start_trac.nc", NF_NOWRITE,nid_trac)
+         IF (ierr_file .NE.NF_NOERR) THEN
+            write(6,*)' Pb d''ouverture du fichier start_trac.nc'
+            write(6,*)' ierr = ', ierr_file 
+         ENDIF
+      END IF
+
+      IF(nq.GE.1) THEN
+      do iq=1,nq 
+
+         IF (config_inca == 'none') THEN
+            ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
+            IF (ierr .NE. NF_NOERR) THEN
+               PRINT*, "Variable  tname(iq) n est pas definie"
+               CALL abort
+            ENDIF
+#ifdef NC_DOUBLE
+            ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
+#else
+            ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
+#endif
+        ELSE ! config_inca = 'chem' ou 'aero'
+! lecture de la valeur du traceur dans start_trac.nc
+           IF (ierr_file .ne. 2) THEN
+             ierr = NF_INQ_VARID (nid_trac, tname(iq), nvarid_trac)
+             IF (ierr .NE. NF_NOERR) THEN
+                PRINT*, tname(iq),"est absent de start_trac.nc"
+                ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
+                IF (ierr .NE. NF_NOERR) THEN
+                   PRINT*, "Variable ", tname(iq)," n est pas definie"
+                   CALL abort
+                ENDIF
+#ifdef NC_DOUBLE
+                ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
+#else
+                ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
+#endif
+                
+             ELSE
+                PRINT*, tname(iq), "est present dans start_trac.nc"
+#ifdef NC_DOUBLE
+               ierr = NF_GET_VAR_DOUBLE(nid_trac, nvarid_trac, trac_tmp)
+#else
+               ierr = NF_GET_VAR_REAL(nid_trac, nvarid_trac, trac_tmp)
+#endif
+                IF (ierr .NE. NF_NOERR) THEN
+                   PRINT*, "Lecture echouee pour", tname(iq)
+                   CALL abort
+                ENDIF
+                ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
+                IF (ierr .NE. NF_NOERR) THEN
+                   PRINT*, "Variable ", tname(iq)," n est pas definie"
+                   CALL abort
+                ENDIF
+#ifdef NC_DOUBLE
+                ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,trac_tmp)
+#else
+                ierr = NF_PUT_VAR_REAL (nid,nvarid,trac_tmp)
+#endif
+               
+             ENDIF ! IF (ierr .NE. NF_NOERR)
+! fin lecture du traceur
+          ELSE                  ! si il n'y a pas de fichier start_trac.nc
+!             print *, 'il n y a pas de fichier start_trac'
+             ierr = NF_INQ_VARID(nid, tname(iq), nvarid)
+             IF (ierr .NE. NF_NOERR) THEN
+                PRINT*, "Variable  tname(iq) n est pas definie"
+                CALL abort
+             ENDIF
+#ifdef NC_DOUBLE
+             ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q(1,1,iq))
+#else
+             ierr = NF_PUT_VAR_REAL (nid,nvarid,q(1,1,iq))
+#endif
+          ENDIF ! (ierr_file .ne. 2)
+       END IF   ! config_inca
+      
+      ENDDO
+      ENDIF
+c
+      ierr = NF_INQ_VARID(nid, "masse", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Variable masse n est pas definie"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,masse)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,masse)
+#endif
+c
+      ierr = NF_INQ_VARID(nid, "ps", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "Variable ps n est pas definie"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,ps)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,ps)
+#endif
+
+      ierr = NF_CLOSE(nid)
+c
+      RETURN
+      END
+
Index: /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F	(revision 1000)
@@ -0,0 +1,753 @@
+!
+! $Header$
+!
+c
+c
+      SUBROUTINE etat0_netcdf (interbar, masque)
+    
+      USE startvar
+      USE ioipsl
+      USE dimphy
+      USE fonte_neige_mod
+      USE pbl_surface_mod
+      USE phys_state_var_mod
+      !
+      IMPLICIT NONE
+      !
+#include "netcdf.inc"
+#include "dimensions.h"
+#include "paramet.h"
+      !
+      !
+!      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2, 
+!     .KLON=KFDIA-KIDIA+1,KLEV=llm
+      !
+#include "comgeom2.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "indicesol.h"
+#include "dimsoil.h"
+#include "temps.h"
+      !
+      LOGICAL interbar
+      REAL :: latfi(klon), lonfi(klon)
+      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
+     . psol(iip1, jjp1), phis(iip1, jjp1)
+      REAL :: p3d(iip1, jjp1, llm+1)
+      REAL :: uvent(iip1, jjp1, llm)
+      REAL :: vvent(iip1, jjm, llm)
+      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
+      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
+      REAL :: tsol(klon), qsol(klon), sn(klon)
+      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf) 
+      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
+      REAL :: alblw(klon,nbsrf)
+      REAL :: tsoil(klon,nsoilmx,nbsrf) 
+      REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
+      REAL :: rugmer(klon)
+      REAL :: qd(iip1, jjp1, llm)
+      REAL :: run_off_lic_0(klon)
+      ! declarations pour lecture glace de mer
+      REAL :: rugv(klon)
+      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
+      INTEGER :: itaul(1), fid
+      REAL :: lev(1), date
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
+      REAL :: flic_tmp(iip1, jjp1)
+      REAL :: champint(iim, jjp1)
+      !
+
+      CHARACTER*80 :: varname
+      !
+      INTEGER :: i,j, ig, l, ji,ii1,ii2
+      INTEGER :: nq
+      REAL :: xpi
+      !
+      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
+      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
+      REAL :: workvar(iip1,jjp1,llm)
+      !
+      REAL ::  prefkap, unskap
+      !
+      real :: time_step,t_ops,t_wrt
+
+#include "comdissnew.h"
+#include "control.h"
+#include "serre.h"
+#include "clesphys.h"
+
+      INTEGER  ::        longcles
+      PARAMETER      ( longcles  = 20 )
+      REAL :: clesphy0 ( longcles       )
+      REAL :: p(iip1,jjp1,llm)
+      INTEGER :: itau, iday
+      REAL :: masse(iip1,jjp1,llm)
+      REAL :: xpn,xps,xppn(iim),xpps(iim)
+      real :: time
+      REAL :: phi(ip1jmp1,llm)
+      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
+      REAL :: w(ip1jmp1,llm)
+      REAL ::phystep
+      REAL :: rugsrel(iip1*jjp1)
+      REAL :: fder(klon)
+      real zrel(iip1*jjp1),chmin,chmax
+
+      CHARACTER*80 :: visu_file
+      INTEGER :: visuid
+
+! pour la lecture du fichier masque ocean
+      integer :: nid_o2a
+      logical :: couple = .false.
+      INTEGER :: iml_omask, jml_omask
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
+      real, dimension(klon) :: ocemask_fi
+      integer :: isst(klon-2)
+      real zx_tmp_2d(iim,jjp1)
+
+      REAL :: dummy
+
+      logical              :: ok_newmicro
+      integer              :: iflag_radia
+      logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
+      LOGICAL              :: ok_ade, ok_aie, aerosol_couple
+      REAL                 :: bl95_b0, bl95_b1
+      real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
+      integer              :: iflag_cldcon
+      integer              :: iflag_ratqs
+      integer :: iflag_coupl
+      integer :: iflag_clos
+      integer :: iflag_wake
+      integer :: iflag_thermals,nsplit_thermals
+      real    :: tau_thermals
+      REAL      :: solarlong0
+      real :: seuil_inversion
+
+      !
+      !   Constantes 
+      !
+      pi     = 4. * ATAN(1.)
+      rad    = 6371229.
+      omeg   = 4.* ASIN(1.)/(24.*3600.)
+      g      = 9.8
+      daysec = 86400.
+      kappa  = 0.2857143
+      cpp    = 1004.70885
+      !
+      preff     = 101325.
+      unskap = 1./kappa
+      !
+      jmp1    = jjm + 1
+      !
+      !    Construct a grid
+      !
+
+!      CALL defrun_new(99,.TRUE.,clesphy0)
+      CALL conf_gcm( 99, .TRUE. , clesphy0 )
+      call conf_phys( ok_veget, ok_journe, ok_mensuel, ok_instan,       &
+     &                 ok_hf,                                           &
+     &                 solarlong0,seuil_inversion,                      &
+     &                 fact_cldcon, facttemps,ok_newmicro,iflag_radia,  &
+     &                 iflag_cldcon,                                    &
+     &                 iflag_ratqs,ratqsbas,ratqshaut,                  &
+     &                 ok_ade, ok_aie, aerosol_couple,                  &
+     &                 bl95_b0, bl95_b1,                                &
+     &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
+     &                 iflag_coupl,iflag_clos,iflag_wake )
+
+      dtvr   = daysec/FLOAT(day_step)
+      print*,'dtvr',dtvr
+
+      CALL inicons0()
+      CALL inigeom()
+      !
+      CALL inifilr()
+!      CALL phys_state_var_init()
+      !
+      latfi(1) = ASIN(1.0)
+      DO j = 2, jjm
+        DO i = 1, iim
+          latfi((j-2)*iim+1+i)=  rlatu(j)
+        ENDDO
+      ENDDO
+      latfi(klon) = - ASIN(1.0)
+      !
+      lonfi(1) = 0.0
+      DO j = 2, jjm
+        DO i = 1, iim
+          lonfi((j-2)*iim+1+i) =  rlonv(i)
+        ENDDO
+      ENDDO
+      lonfi(klon) = 0.0
+      !
+      xpi = 2.0 * ASIN(1.0)
+      DO ig = 1, klon
+        latfi(ig) = latfi(ig) * 180.0 / xpi
+        lonfi(ig) = lonfi(ig) * 180.0 / xpi
+      ENDDO
+      !
+      rlat(1) = ASIN(1.0)
+      DO j = 2, jjm
+        DO i = 1, iim
+          rlat((j-2)*iim+1+i)=  rlatu(j)
+        ENDDO
+      ENDDO
+      rlat(klon) = - ASIN(1.0)
+      !
+      rlon(1) = 0.0
+      DO j = 2, jjm
+        DO i = 1, iim
+          rlon((j-2)*iim+1+i) =  rlonv(i)
+        ENDDO
+      ENDDO
+      rlon(klon) = 0.0
+      !
+      xpi = 2.0 * ASIN(1.0)
+      DO ig = 1, klon
+        rlat(ig) = rlat(ig) * 180.0 / xpi
+        rlon(ig) = rlon(ig) * 180.0 / xpi
+      ENDDO
+      !
+      
+
+
+C
+C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
+C utilise pour calculer les poids et pour assurer l'adequation entre les
+C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque 
+C a partir du fichier relief
+C
+
+      write(*,*)'Essai de lecture masque ocean'
+      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
+      if (iret .ne. 0) then
+        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
+        write(*,*)'Run force'
+        varname = 'masque'
+        masque(:,:) = 0.0
+        CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
+     ,  jjm ,rlonu,rlatv , interbar )
+        WRITE(*,*) 'MASQUE construit : Masque'
+        WRITE(*,'(97I1)') nINT(masque(:,:))
+        call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
+        WHERE (zmasq(1 : klon) .LT. EPSFRA)
+            zmasq(1 : klon) = 0.
+        END WHERE 
+        WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
+            zmasq(1 : klon) = 1.
+        END WHERE 
+      else
+        couple = .true.
+        iret = nf_close(nid_o2a)
+        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
+     $    , nid_o2a)
+        if (iml_omask /= iim .or. jml_omask /= jjp1) then
+          write(*,*)'Dimensions non compatibles pour masque ocean'
+          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
+          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+          stop
+        endif
+        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(dlon_omask(iml_omask), stat=iret)
+        ALLOCATE(dlat_omask(jml_omask), stat=iret)
+        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
+        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
+        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
+     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
+        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 
+     $      ttm_tmp, 1, 1, ocetmp)
+        CALL flinclo(fid)
+        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
+        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
+        ocemask = ocetmp
+        if (dlat_omask(1) < dlat_omask(jml_omask)) then
+          do j = 1, jml_omask
+            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
+          enddo
+        endif 
+C
+C passage masque ocean a la grille physique
+C
+        write(*,*)'ocemask '
+        write(*,'(96i1)')int(ocemask)
+        ocemask_fi(1) = ocemask(1,1)
+        do j = 2, jjm
+          do i = 1, iim
+            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
+          enddo
+        enddo
+        ocemask_fi(klon) = ocemask(1,jjp1)
+        zmasq = 1. - ocemask_fi
+      endif
+
+      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
+
+      varname = 'relief'
+      ! This line needs to be replaced by a call to restget to get the values in the restart file
+      orog(:,:) = 0.0
+       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
+     , jjm ,rlonu,rlatv , interbar, masque )
+      !
+      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
+!      WRITE(*,'(49I1)') INT(orog(:,:))
+      !
+      varname = 'rugosite'
+      ! This line needs to be replaced by a call to restget to get the values in the restart file
+      rugo(:,:) = 0.0
+       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
+     , jjm, rlonu,rlatv , interbar )
+      !
+      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite' 
+!      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
+      !
+C
+C on initialise les sous surfaces
+C
+      pctsrf=0.
+c
+      varname = 'psol'
+      psol(:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
+     , jjm ,rlonu,rlatv , interbar )
+      !
+      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM 
+      !  anyway.
+      !
+!      WRITE(*,*) 'PSOL :', psol(10,20)
+!      WRITE(*,*) ap(:), bp(:)
+      CALL pression(ip1jmp1, ap, bp, psol, p3d)
+!      WRITE(*,*) 'P3D :', p3d(10,20,:)
+      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
+!      WRITE(*,*) 'PK:', pk(10,20,:)
+      !
+      !
+      !
+      prefkap =  preff  ** kappa
+!      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
+      DO l = 1, llm
+        DO j=1,jjp1
+          DO i =1, iip1
+            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+           ENDDO
+        ENDDO
+      ENDDO
+      !
+!      WRITE(*,*) 'PLS :', pls(10,20,:)
+      !
+      varname = 'surfgeo'
+      phis(:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
+     , jjm ,rlonu,rlatv, interbar )
+      !
+      varname = 'u'
+      uvent(:,:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
+     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
+      !  
+      varname = 'v'
+      vvent(:,:,:) = 0.0
+      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
+     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
+      !
+      varname = 't'
+      t3d(:,:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
+     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
+      !
+      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
+     .                          maxval(t3d(:,:,:))
+      varname = 'tpot'
+      tpot(:,:,:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
+     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
+      !
+      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
+     .                          maxval(t3d(:,:,:))
+      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
+     .                          maxval(pls(:,:,:))
+
+c Calcul de l'humidite a saturation
+      print*,'avant q_sat'
+      call q_sat(llm*jjp1*iip1,t3d,pls,qsat)
+      print*,'apres q_sat'
+
+      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
+     .                           maxval(qsat(:,:,:))
+      !
+      WRITE(*,*) 'QSAT :', qsat(10,20,:)
+      !
+      varname = 'q'
+      qd(:,:,:) = 0.0
+      q3d(:,:,:,:) = 0.0
+      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
+     .                           maxval(qsat(:,:,:))
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
+     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
+      q3d(:,:,:,1) = qd(:,:,:)
+      !
+      varname = 'tsol'
+      ! This line needs to be replaced by a call to restget to get the values in the restart file
+      tsol(:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
+     .    jjm, rlonu, rlatv , interbar )
+      !
+      WRITE(*,*) 'TSOL construit :'
+!      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
+      !
+      varname = 'qsol'
+      qsol(:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
+     .   jjm, rlonu, rlatv , interbar )
+      !
+      varname = 'snow'
+      sn(:) = 0.0
+      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
+     .    jjm, rlonu, rlatv , interbar )
+      !
+      varname = 'rads'
+      radsol(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
+     .    jjm, rlonu, rlatv , interbar )
+      !
+      varname = 'rugmer'
+      rugmer(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      !
+!      varname = 'agesno'
+!      agesno(:) = 0.0
+!      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
+!     .     jjm, rlonu, rlatv , interbar )
+
+      varname = 'zmea'
+      zmea(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+
+      varname = 'zstd'
+      zstd(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      varname = 'zsig'
+      zsig(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      varname = 'zgam'
+      zgam(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      varname = 'zthe'
+      zthe(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      varname = 'zpic'
+      zpic(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+      varname = 'zval'
+      zval(:) = 0.0
+      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
+     .     jjm, rlonu, rlatv , interbar )
+c
+      rugsrel(:) = 0.0
+      IF(ok_orodr)  THEN
+        DO i = 1, iip1* jjp1
+         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
+        ENDDO
+      ENDIF
+
+
+C
+C lecture du fichier glace de terre pour fixer la fraction de terre 
+C et de glace de terre
+C
+      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
+     $    , fid)
+      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
+      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
+      ALLOCATE(dlon_lic(iml_lic), stat=iret)
+      ALLOCATE(dlat_lic(jml_lic), stat=iret)
+      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
+      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
+     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
+      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
+     $    , 1, 1, fraclic)
+      CALL flinclo(fid)
+C
+C interpolation sur la grille T du modele
+C
+      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 
+     $    iml_lic, jml_lic
+c
+C sil les coordonnees sont en degres, on les transforme
+C
+      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
+          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
+      ENDIF 
+      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 
+          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
+      ENDIF 
+
+      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
+      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 
+C
+      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
+     $    ,iim, jjp1,
+     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
+cx$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
+      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
+C
+C passage sur la grille physique
+C
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
+     $    pctsrf(1:klon, is_lic))
+C adequation avec le maque terre/mer
+c      zmasq(157) = 0.
+      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 
+          pctsrf(1 : klon, is_lic) = 0. 
+      END WHERE
+      WHERE (zmasq( 1 : klon) .LT. EPSFRA) 
+          pctsrf(1 : klon, is_lic) = 0.
+      END WHERE 
+      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
+      DO ji = 1, klon
+        IF (zmasq(ji) .GT. EPSFRA) THEN 
+            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
+                pctsrf(ji, is_lic) = zmasq(ji)
+                pctsrf(ji, is_ter) = 0.
+            ELSE 
+                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
+                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
+                    pctsrf(ji,is_ter) = 0.
+                    pctsrf(ji, is_lic) = zmasq(ji)
+                ENDIF 
+            ENDIF 
+        ENDIF 
+      END DO 
+C
+C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
+C
+      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
+
+
+      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
+          pctsrf(1 : klon, is_oce) = 0.
+      END WHERE 
+
+      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
+
+      isst = 0
+      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
+C
+C verif que somme des sous surface = 1
+C
+      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0) 
+     $    .GT. EPSFRA)
+      IF (ji .NE. 0) THEN
+          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
+      ENDIF 
+
+!      where (pctsrf(1:klon, is_ter) >= .5) 
+!        pctsrf(1:klon, is_ter) = 1.
+!        pctsrf(1:klon, is_oce) = 0.
+!        pctsrf(1:klon, is_sic) = 0.
+!        pctsrf(1:klon, is_lic) = 0.
+!        zmasq = 1.
+!      endwhere
+!      where (pctsrf(1:klon, is_lic) >= .5) 
+!        pctsrf(1:klon, is_ter) = 0.
+!        pctsrf(1:klon, is_oce) = 0.
+!        pctsrf(1:klon, is_sic) = 0.
+!        pctsrf(1:klon, is_lic) = 1.
+!        zmasq = 1.
+!      endwhere
+!      where (pctsrf(1:klon, is_oce) >= .5) 
+!        pctsrf(1:klon, is_ter) = 0.
+!        pctsrf(1:klon, is_oce) = 1.
+!        pctsrf(1:klon, is_sic) = 0.
+!        pctsrf(1:klon, is_lic) = 0.
+!        zmasq = 0.
+!      endwhere
+!      where (pctsrf(1:klon, is_sic) >= .5) 
+!        pctsrf(1:klon, is_ter) = 0.
+!        pctsrf(1:klon, is_oce) = 0.
+!        pctsrf(1:klon, is_sic) = 1.
+!        pctsrf(1:klon, is_lic) = 0.
+!        zmasq = 0.
+!      endwhere
+!      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
+C
+C verif que somme des sous surface = 1
+C
+!      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 ) 
+!     $    .GT. EPSFRA)
+!      IF (ji .NE. 0) THEN
+!          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
+!     ENDIF 
+
+      CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
+      write(*,*)'zmasq = '
+      write(*,'(96i1)')nint(zx_tmp_2d)
+      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
+      WRITE(*,*) 'MASQUE construit : Masque'
+      WRITE(*,'(97I1)') nINT(masque(:,:))
+
+
+
+C Calcul intermediaire
+c 
+      CALL massdair( p3d, masse  )
+c
+
+      print *,' ALPHAX ',alphax
+
+      DO  l = 1, llm
+        DO  i    = 1, iim
+          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
+          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
+        ENDDO
+          xpn      = SUM(xppn)/apoln
+          xps      = SUM(xpps)/apols
+        DO i   = 1, iip1
+          masse(   i   ,   1     ,  l )   = xpn
+          masse(   i   ,   jjp1  ,  l )   = xps
+        ENDDO
+      ENDDO
+      q3d(iip1,:,:,:) = q3d(1,:,:,:)
+      phis(iip1,:) = phis(1,:)
+
+C init pour traceurs
+      call iniadvtrac(nq)
+C Ecriture
+      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
+     *                tetagdiv, tetagrot , tetatemp              )
+      print*,'sortie inidissip'
+      itau = 0
+      itau_dyn = 0
+      itau_phy = 0
+      iday = dayref +itau/day_step
+      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
+c     
+      IF(time.GT.1)  THEN
+       time = time - 1
+       iday = iday + 1
+      ENDIF
+      day_ref = dayref
+      annee_ref = anneeref
+
+      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
+      print*,'sortie geopot'
+      
+      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
+     *                phi,w, pbaru,pbarv,time+iday-dayref   )
+       print*,'sortie caldyn0'     
+      CALL dynredem0("start.nc",dayref,phis,nqmx)
+      print*,'sortie dynredem0'
+      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
+     .                            psol)
+      print*,'sortie dynredem1' 
+C
+C Ecriture etat initial physique
+C
+      write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad
+      phystep   = dtvr * FLOAT(iphysiq)
+      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
+      write(*,*)'phystep =', phystep, radpas
+cIM : lecture de co2_ppm & solaire ds physiq.def
+c     co2_ppm   = 348.0
+c     solaire   = 1365.0
+
+c
+c Initialisation 
+c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
+c
+      ftsol(:,is_ter) = tsol
+      ftsol(:,is_lic) = tsol
+      ftsol(:,is_oce) = tsol
+      ftsol(:,is_sic) = tsol
+      snsrf(:,is_ter) = sn
+      snsrf(:,is_lic) = sn
+      snsrf(:,is_oce) = sn
+      snsrf(:,is_sic) = sn
+      falb1(:,is_ter) = 0.08
+      falb1(:,is_lic) = 0.6
+      falb1(:,is_oce) = 0.5
+      falb1(:,is_sic) = 0.6
+      falb2 = falb1
+      evap(:,:) = 0.
+      qsolsrf(:,is_ter) = 150
+      qsolsrf(:,is_lic) = 150
+      qsolsrf(:,is_oce) = 150.
+      qsolsrf(:,is_sic) = 150.
+      do i = 1, nbsrf
+        do j = 1, nsoilmx
+          tsoil(:,j,i) = tsol
+        enddo
+      enddo
+      rain_fall = 0.; snow_fall = 0.
+      solsw = 165.
+      sollw = -53.
+      t_ancien = 273.15
+      q_ancien = 0.
+      agesno = 0.
+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)
+      frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
+      frugs(1:klon,is_sic) = 0.001
+      fder = 0.0
+      clwcon = 0.0
+      rnebcon = 0.0
+      ratqs = 0.0
+      run_off_lic_0 = 0.0 
+      rugoro = 0.0
+
+c
+c Avant l'appel a phyredem, on initialize les modules de surface
+c avec les valeurs qui vont etre ecrit dans startphy.nc
+c
+      dummy = 1.0
+      pbl_tke(:,:,:) = 1.e-8 
+      zmax0(:) = 40.
+      f0(:) = 1.e-5
+      ema_work1(:,:) = 0.
+      ema_work2(:,:) = 0.
+      wake_deltat(:,:) = 0.
+      wake_deltaq(:,:) = 0.
+      wake_s(:) = 0.
+      wake_cstar(:) = 0.
+      wake_fip(:) = 0.
+
+      call fonte_neige_init(run_off_lic_0)
+      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
+     $     evap, frugs, agesno, tsoil)
+
+      call phyredem("startphy.nc")
+
+
+
+C     Sortie Visu pour les champs dynamiques
+      if (1.eq.0 ) then
+      print*,'sortie visu'
+      time_step = 1.
+      t_ops = 2.
+      t_wrt = 2.
+      itau = 2.
+      visu_file='Etat0_visu.nc'
+      CALL initdynav(visu_file,dayref,anneeref,time_step,
+     .              t_ops, t_wrt, nqmx, visuid)
+      CALL writedynav(visuid, nqmx, itau,vvent ,
+     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
+      else
+         print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
+      endif
+      print*,'entree histclo'
+      CALL histclo
+      RETURN
+      !
+      END SUBROUTINE etat0_netcdf
+
Index: /LMDZ4/trunk/libf/dyn3dpar/flumass.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/flumass.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/flumass.F	(revision 1000)
@@ -0,0 +1,109 @@
+!
+! $Header$
+!
+      SUBROUTINE flumass (massebx,masseby, vcont, ucont, pbaru, pbarv )
+
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteurs:  P. Le Van, F. Hourdin  .
+c   -------
+c
+c   Objet:
+c   ------
+c
+c *********************************************************************
+c     .... calcul du flux de masse  aux niveaux s ......
+c *********************************************************************
+c   massebx,masseby,vcont et ucont sont des argum. d'entree pour le s-pg .
+c       pbaru  et pbarv            sont des argum.de sortie pour le s-pg .
+c
+c=======================================================================
+
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+
+      REAL massebx( ip1jmp1,llm ),masseby( ip1jm,llm ) ,
+     * vcont( ip1jm,llm ),ucont( ip1jmp1,llm ),pbaru( ip1jmp1,llm ),
+     * pbarv( ip1jm,llm )
+
+      REAL apbarun( iip1 ),apbarus( iip1 )
+
+      REAL sairen,saireun,saires,saireus,ctn,cts,ctn0,cts0
+      INTEGER  l,ij,i
+
+      REAL       SSUM
+
+
+      DO  5 l = 1,llm
+
+      DO  1 ij = iip2,ip1jm
+      pbaru( ij,l ) = massebx( ij,l ) * ucont( ij,l )
+   1  CONTINUE
+
+      DO 3 ij = 1,ip1jm
+      pbarv( ij,l ) = masseby( ij,l ) * vcont( ij,l )
+   3  CONTINUE
+
+   5  CONTINUE
+
+c    ................................................................
+c     calcul de la composante du flux de masse en x aux poles .......
+c    ................................................................
+c     par la resolution d'1 systeme de 2 equations .
+
+c     la premiere equat.decrivant le calcul de la divergence en 1 point i
+c     du pole,ce calcul etant itere de i=1 a i=im .
+c                 c.a.d   ,
+c     ( ( 0.5*pbaru(i)-0.5*pbaru(i-1) - pbarv(i))/aire(i)   =
+c                                           - somme de ( pbarv(n) )/aire pole
+
+c     l'autre equat.specifiant que la moyenne du flux de masse au pole est =0.
+c     c.a.d    somme de pbaru(n)*aire locale(n) = 0.
+
+c     on en revient ainsi a determiner la constante additive commune aux pbaru
+c     qui representait pbaru(0,j,l) dans l'equat.du calcul de la diverg.au pt
+c     i=1 .
+c     i variant de 1 a im
+c     n variant de 1 a im
+
+      sairen = SSUM( iim,  aire(   1     ), 1 )
+      saireun= SSUM( iim, aireu(   1     ), 1 )
+      saires = SSUM( iim,  aire( ip1jm+1 ), 1 )
+      saireus= SSUM( iim, aireu( ip1jm+1 ), 1 )
+
+      DO 20 l = 1,llm
+
+      ctn =  SSUM( iim, pbarv(    1     ,l),  1 )/ sairen
+      cts =  SSUM( iim, pbarv(ip1jmi1+ 1,l),  1 )/ saires
+
+      pbaru(    1   ,l )=   pbarv(    1     ,l ) - ctn * aire(    1    )
+      pbaru( ip1jm+1,l )= - pbarv( ip1jmi1+1,l ) + cts * aire( ip1jm+1 )
+
+      DO 11 i = 2,iim
+      pbaru(    i    ,l ) = pbaru(   i - 1   ,l )    +
+     *                      pbarv(    i      ,l ) - ctn * aire(   i    )
+
+      pbaru( i+ ip1jm,l ) = pbaru( i+ ip1jm-1,l )    -
+     *                      pbarv( i+ ip1jmi1,l ) + cts * aire(i+ ip1jm)
+  11  CONTINUE
+      DO 12 i = 1,iim
+      apbarun(i) = aireu(    i   ) * pbaru(   i    , l)
+      apbarus(i) = aireu(i +ip1jm) * pbaru(i +ip1jm, l)
+  12  CONTINUE
+      ctn0 = -SSUM( iim,apbarun,1 )/saireun
+      cts0 = -SSUM( iim,apbarus,1 )/saireus
+      DO 14 i = 1,iim
+      pbaru(   i    , l) = 2. * ( pbaru(   i    , l) + ctn0 )
+      pbaru(i+ ip1jm, l) = 2. * ( pbaru(i +ip1jm, l) + cts0 )
+  14  CONTINUE
+
+      pbaru(   iip1 ,l ) = pbaru(    1    ,l )
+      pbaru( ip1jmp1,l ) = pbaru( ip1jm +1,l )
+  20  CONTINUE
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/fluxstokenc_p.F	(revision 1000)
@@ -27,5 +27,5 @@
       REAL phis(ip1jmp1)
 
-      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
+      REAL,SAVE :: pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
       REAL massem(ip1jmp1,llm),tetac(ip1jmp1,llm),phic(ip1jmp1,llm)
 
@@ -55,4 +55,6 @@
 cym      wg(:,:)       = 0.
 
+c$OMP MASTER
+
       if(first) then
 
@@ -98,5 +100,5 @@
 	pbaruc(ijb:ije,1:llm)=0
 	
-	if (pole_sud) ije=ij_end-iip1
+	IF (pole_sud) ije=ij_end-iip1
 	pbarvc(ijb:ije,1:llm)=0
       ENDIF
@@ -134,7 +136,9 @@
       iadvtr   = iadvtr+1
 
-
+c$OMP END MASTER
+c$OMP BARRIER
 c   Test pour savoir si on advecte a ce pas de temps
       IF ( iadvtr.EQ.istdyn ) THEN
+c$OMP MASTER
 c    normalisation
       ijb=ij_begin
@@ -162,10 +166,13 @@
 c     1. calcul de w
 c     2. groupement des mailles pres du pole.
-
+c$OMP END MASTER
+c$OMP BARRIER 
         call Register_Hallo(pbaruc,ip1jmp1,llm,1,1,1,1,Req)
 	call Register_Hallo(pbarvc,ip1jm,llm,1,1,1,1,Req)
         call SendRequest(Req)
+c$OMP BARRIER
         call WaitRequest(Req)
-
+c$OMP BARRIER
+c$OMP MASTER
         CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
         
@@ -226,5 +233,5 @@
 	
 C
-
+c$OMP END MASTER
       ENDIF ! if iadvtr.EQ.istdyn
 
Index: /LMDZ4/trunk/libf/dyn3dpar/gcm.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/gcm.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/gcm.F	(revision 1000)
@@ -49,5 +49,4 @@
 c   Declarations:
 c   -------------
-      include 'mpif.h'
 #include "dimensions.h"
 #include "paramet.h"
@@ -220,5 +219,6 @@
       CALL set_bands
       CALL Init_interface_dyn_phys
-      call MPI_BARRIER(COMM_LMDZ,ierr)
+      CALL barrier
+
       if (mpi_rank==0) call WriteBands
       call SetDistrib(jj_Nb_Caldyn)
Index: /LMDZ4/trunk/libf/dyn3dpar/gr_u_scal_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/gr_u_scal_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/gr_u_scal_p.F	(revision 1000)
@@ -50,8 +50,7 @@
       ijb=ij_begin
       ije=ij_end
-      if (pole_nord) ijb=ij_begin+iip1
       
       DO l=1,nx
-         DO ij=ijb,ije
+         DO ij=ijb+1,ije
             x_scal(ij,l)=
      s      (aireu(ij)*x_u(ij,l)+aireu(ij-1)*x_u(ij-1,l))
Index: /LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/initdynav_p.F	(revision 1000)
@@ -70,4 +70,17 @@
       integer zan, dayref
       integer :: jjb,jje,jjn
+
+! definition du domaine d'ecriture pour le rebuild
+
+      INTEGER,DIMENSION(2) :: ddid
+      INTEGER,DIMENSION(2) :: dsg
+      INTEGER,DIMENSION(2) :: dsl
+      INTEGER,DIMENSION(2) :: dpf
+      INTEGER,DIMENSION(2) :: dpl
+      INTEGER,DIMENSION(2) :: dhs
+      INTEGER,DIMENSION(2) :: dhe 
+      
+      INTEGER :: dynave_domain_id
+      
       
       if (adjust) return
@@ -95,9 +108,19 @@
       jje=jj_end
       jjn=jj_nb
+
+      ddid=(/ 1,2 /)
+      dsg=(/ iip1,jjp1 /)
+      dsl=(/ iip1,jjn /)
+      dpf=(/ 1,jjb /)
+      dpl=(/ iip1,jje /)
+      dhs=(/ 0,0 /)
+      dhe=(/ 0,0 /)
+
+      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, 
+     .                 'box',dynave_domain_id)
              
-      call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc',
-     .             iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
-     .             1, iip1, 1, jjn,
-     .             tau0, zjulian, tstep, thoriid, fileid)
+      call histbeg(trim(infile),iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
+     .             1, iip1, 1, jjn,tau0, zjulian, tstep, thoriid,
+     .             fileid,dynave_domain_id)
 
 C
Index: /LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/initfluxsto_p.F	(revision 1000)
@@ -74,4 +74,19 @@
       logical ok_sync
       integer :: jjb,jje,jjn
+
+! definition du domaine d'ecriture pour le rebuild
+
+      INTEGER,DIMENSION(2) :: ddid
+      INTEGER,DIMENSION(2) :: dsg
+      INTEGER,DIMENSION(2) :: dsl
+      INTEGER,DIMENSION(2) :: dpf
+      INTEGER,DIMENSION(2) :: dpl
+      INTEGER,DIMENSION(2) :: dhs
+      INTEGER,DIMENSION(2) :: dhe 
+      
+      INTEGER :: dynu_domain_id
+      INTEGER :: dynv_domain_id
+
+
 C
 C  Initialisations
@@ -100,9 +115,19 @@
       jje=jj_end
       jjn=jj_nb
+
+      ddid=(/ 1,2 /)
+      dsg=(/ iip1,jjp1 /)
+      dsl=(/ iip1,jjn /)
+      dpf=(/ 1,jjb /)
+      dpl=(/ iip1,jje /)
+      dhs=(/ 0,0 /)
+      dhe=(/ 0,0 /)
+
+      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, 
+     .                 'box',dynu_domain_id)
        
-      call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc',
-     .             iip1, rlong(:,1), jjp1, rlat(1,jjb:jje),
-     .             1, iip1, 1, jjn,
-     .             tau0, zjulian, tstep, uhoriid, fileid)
+      call histbeg(trim(infile),iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
+     .             1, iip1, 1, jjn, tau0, zjulian, tstep, uhoriid,
+     .             fileid,dynu_domain_id)
 C
 C  Creation du fichier histoire pour la grille en V (oblige pour l'instant,
@@ -124,8 +149,18 @@
       if (pole_sud) jjn=jj_nb-1
 
-      call histbeg('fluxstokev_'//trim(int2str(mpi_rank))//'.nc',
-     .             iip1, rlong(:,1), jjm, rlat(1,jjb:jje),
-     .             1, iip1, 1, jjn,
-     .             tau0, zjulian, tstep, vhoriid, filevid)
+      ddid=(/ 1,2 /)
+      dsg=(/ iip1,jjm /)
+      dsl=(/ iip1,jjn /)
+      dpf=(/ 1,jjb /)
+      dpl=(/ iip1,jje /)
+      dhs=(/ 0,0 /)
+      dhe=(/ 0,0 /)
+
+      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, 
+     .                 'box',dynv_domain_id)
+     
+      call histbeg('fluxstokev',iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
+     .             1, iip1, 1, jjn,tau0, zjulian, tstep, vhoriid,
+     .             filevid,dynv_domain_id)
 	
       rl(1,1) = 1.	
Index: /LMDZ4/trunk/libf/dyn3dpar/inithist_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/inithist_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/inithist_p.F	(revision 1000)
@@ -70,4 +70,17 @@
       integer zan, dayref
       integer :: jjb,jje,jjn
+
+! definition du domaine d'ecriture pour le rebuild
+
+      INTEGER,DIMENSION(2) :: ddid
+      INTEGER,DIMENSION(2) :: dsg
+      INTEGER,DIMENSION(2) :: dsl
+      INTEGER,DIMENSION(2) :: dpf
+      INTEGER,DIMENSION(2) :: dpl
+      INTEGER,DIMENSION(2) :: dhs
+      INTEGER,DIMENSION(2) :: dhe 
+      
+      INTEGER :: dynu_domain_id
+      INTEGER :: dynv_domain_id
 C
 C  Initialisations
@@ -95,9 +108,20 @@
       jje=jj_end
       jjn=jj_nb
-      
-       call histbeg(trim(infile)//'_'//trim(int2str(mpi_rank))//'.nc',
-     .             iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
-     .             1, iip1, 1, jjn,
-     .             tau0, zjulian, tstep, uhoriid, fileid)
+
+
+      ddid=(/ 1,2 /)
+      dsg=(/ iip1,jjp1 /)
+      dsl=(/ iip1,jjn /)
+      dpf=(/ 1,jjb /)
+      dpl=(/ iip1,jje /)
+      dhs=(/ 0,0 /)
+      dhe=(/ 0,0 /)
+
+      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, 
+     .                 'box',dynu_domain_id)
+      
+       call histbeg(trim(infile),iip1, rlong(:,1), jjn, 
+     .              rlat(1,jjb:jje), 1, iip1, 1, jjn, tau0,
+     .              zjulian, tstep, uhoriid, fileid,dynu_domain_id)
 C
 C  Creation du fichier histoire pour la grille en V (oblige pour l'instant,
@@ -117,9 +141,19 @@
       if (pole_sud) jje=jj_end-1
       if (pole_sud) jjn=jj_nb-1
-      
-      call histbeg('dyn_histv_'//trim(int2str(mpi_rank))//'.nc', 
-     .             iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
-     .             1, iip1, 1, jjn,
-     .             tau0, zjulian, tstep, vhoriid, filevid)
+
+      ddid=(/ 1,2 /)
+      dsg=(/ iip1,jjm /)
+      dsl=(/ iip1,jjn /)
+      dpf=(/ 1,jjb /)
+      dpl=(/ iip1,jje /)
+      dhs=(/ 0,0 /)
+      dhe=(/ 0,0 /)
+
+      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe, 
+     .                 'box',dynv_domain_id)
+      
+      call histbeg('dyn_histv', iip1, rlong(:,1), jjn, rlat(1,jjb:jje),
+     .             1, iip1, 1, jjn, tau0, zjulian, tstep, vhoriid, 
+     .             filevid,dynv_domain_id)
 C
 C  Appel a histhori pour rajouter les autres grilles horizontales
Index: /LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F	(revision 1000)
@@ -76,5 +76,4 @@
 #include "advtrac.h"
       
-      include 'mpif.h'
       integer nq
 
@@ -237,8 +236,6 @@
 c$OMP MASTER
 
-!$OMP CRITICAL (MPI)
-      call MPI_BARRIER(COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
-
+      CALL barrier
+      
 c$OMP END MASTER
 c$OMP BARRIER
@@ -736,5 +733,4 @@
         call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
      *                               jj_Nb_physic,2,2,Request_physic)
-        call SetDistrib(jj_nb_Physic)
         
         call SendRequest(Request_Physic)
@@ -853,4 +849,8 @@
 c      ajout des tendances physiques:
 c      ------------------------------
+         IF (ok_strato) THEN
+           CALL top_bound_p( vcov,ucov,teta, dufi,dvfi,dtetafi)
+         ENDIF
+       
           CALL addfi_p( nqmx, dtphys, leapf, forward   ,
      $                  ucov, vcov, teta , q   ,ps ,
@@ -1289,4 +1289,5 @@
 c$OMP BARRIER
               call WaitRequest(TestRequest)
+c$OMP BARRIER
 c$OMP MASTER
               CALL writedynav_p(histaveid, nqmx, itau,vcov ,
@@ -1340,16 +1341,5 @@
         CALL writehist_p(histid,histvid, nqmx,itau,vcov, 
      s                       ucov,teta,phi,q,masse,ps,phis)
-c#else
-c       call Gather_Field(unat,ip1jmp1,llm,0)
-c       call Gather_Field(vnat,ip1jm,llm,0)
-c       call Gather_Field(teta,ip1jmp1,llm,0)
-c       call Gather_Field(ps,ip1jmp1,1,0)
-c       do iq=1,nqmx
-c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
-c       enddo
-c      
-c       if (mpi_rank==0) then
-c#include "write_grads_dyn.h"
-c       endif
+
 #endif
 c$OMP END MASTER
@@ -1445,4 +1435,5 @@
               call WaitRequest(TestRequest)
 
+c$OMP BARRIER
 c$OMP MASTER
               CALL writedynav_p(histaveid, nqmx, itau,vcov ,
Index: /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/limit_netcdf.F	(revision 1000)
@@ -0,0 +1,1286 @@
+!
+! $Header$
+!
+C
+C
+      SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
+      USE dimphy
+      use phys_state_var_mod , ONLY : pctsrf
+      IMPLICIT none
+c
+c-------------------------------------------------------------
+C Author : L. Fairhead
+C Date   : 27/01/94
+C Objet  : Construction des fichiers de conditions aux limites
+C          pour le nouveau
+C          modele a partir de fichiers de climatologie. Les deux
+C          grilles doivent etre regulieres
+c
+c Modifie par z.x.li (le23mars1994)
+c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
+c                         pour lecture netcdf dans LMDZ.3.3
+c Modifie par P;Le Van  ,  juillet 2001
+c-------------------------------------------------------------
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "control.h"
+#include "logic.h"
+#include "netcdf.inc"
+#include "comvert.h"
+#include "comgeom2.h"
+#include "comconst.h"
+cy#include "dimphy.h"
+#include "indicesol.h"
+c
+c-----------------------------------------------------------------------
+      LOGICAL interbar, extrap, oldice
+
+      REAL phy_nat(klon,360), phy_nat0(klon)
+      REAL phy_alb(klon,360)
+      REAL phy_sst(klon,360)
+      REAL phy_bil(klon,360)
+      REAL phy_rug(klon,360)
+      REAL phy_ice(klon)
+c
+      real pctsrf_t(klon,nbsrf,360)
+
+      REAL verif
+
+      REAL masque(iip1,jjp1)
+      REAL mask(iim,jjp1)
+CPB
+C newlmt indique l'utilisation de la sous-maille fractionnelle
+C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3)
+      LOGICAL newlmt, fracterre
+      PARAMETER(newlmt=.TRUE.)
+      PARAMETER(fracterre = .TRUE.) 
+
+C Declarations pour le champ de depart
+      INTEGER imdep, jmdep,lmdep
+      INTEGER  tbid
+      PARAMETER ( tbid = 60 )        ! >52 semaines
+      REAL  timecoord(tbid)
+c
+      REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:)
+      REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:)
+      REAL , ALLOCATABLE :: dlon(:), dlat(:)
+      REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
+      REAL , ALLOCATABLE :: champ_msk(:), champ(:)
+      REAL , ALLOCATABLE :: work(:,:)
+
+      CHARACTER*25 title
+
+C Declarations pour le champ interpole 2D
+      REAL champint(iim,jjp1)
+      real chmin,chmax
+
+C Declarations pour le champ interpole 3D
+      REAL champtime(iim,jjp1,tbid)
+      REAL timeyear(tbid)
+      REAL champan(iip1,jjp1,366)
+
+C Declarations pour l'inteprolation verticale
+      REAL ax(tbid), ay(tbid)
+      REAL by
+      REAL yder(tbid)
+
+
+      INTEGER ierr
+      INTEGER dimfirst(3)
+      INTEGER dimlast(3)
+c
+      INTEGER nid, ndim, ntim
+      INTEGER dims(2), debut(2), epais(2)
+      INTEGER id_tim
+      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
+CPB
+      INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
+
+      INTEGER i, j, k, l, ji
+c declarations pour lecture glace de mer
+      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
+      INTEGER :: itaul(1), fid
+      REAL :: lev(1), date, dt
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
+      REAL :: flic_tmp(iip1, jjp1)
+
+c Diverses variables locales
+      REAL time
+! pour la lecture du fichier masque ocean
+      integer :: nid_o2a
+      logical :: couple = .false.
+      INTEGER :: iml_omask, jml_omask
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
+      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
+      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
+      real, dimension(klon) :: ocemask_fi
+
+      INTEGER          longcles
+      PARAMETER      ( longcles = 20 )
+      REAL  clesphy0 ( longcles      )
+#include "serre.h"
+      INTEGER ncid,varid,ndimid(4),dimid
+      character*30 namedim
+      CHARACTER*80 :: varname
+
+cIM28/02/2002 <== PM
+      REAL tmidmonth(12)
+      SAVE tmidmonth
+      DATA tmidmonth/15,45,75,105,135,165,195,225,255,285,315,345/
+
+c initialisations:
+      CALL conf_gcm( 99, .TRUE. , clesphy0 )
+
+
+      pi     = 4. * ATAN(1.)
+      rad    = 6 371 229.
+      omeg   = 4.* ASIN(1.)/(24.*3600.)
+      g      = 9.8
+      daysec = 86400.
+      kappa  = 0.2857143
+      cpp    = 1004.70885
+      dtvr    = daysec/FLOAT(day_step)
+      CALL inigeom
+c
+C Traitement du relief au sol
+c
+      write(*,*) 'Traitement du relief au sol pour fabriquer masque'
+      ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)
+
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ALLOCATE( lonmsk_ini(imdep) )
+      ALLOCATE(   dlon_msk(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini)
+#endif
+
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ALLOCATE( latmsk_ini(jmdep) )
+      ALLOCATE(   dlat_msk(jmdep) )
+      ALLOCATE(  champ_msk(imdep*jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini)
+#endif
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
+#endif
+c
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+c
+      title='RELIEF'
+
+      CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini,
+     . dlon_msk, dlat_msk, champ_msk, interbar  )
+
+      DO i = 1, iim
+      DO j = 1, jjp1
+         mask(i,j) = masque(i,j)
+      ENDDO
+      ENDDO
+      WRITE(*,*) 'MASK:'
+      WRITE(*,'(96i1)')INT(mask)     
+      ierr = NF_CLOSE(ncid)
+c
+c
+C Traitement de la rugosite
+c
+      PRINT*, 'Traitement de la rugosite'
+      ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+
+      ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlon_ini(imdep) )
+      ALLOCATE(     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlat_ini(jmdep) )
+      ALLOCATE(     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      print*,'variable ', namedim, 'dimension ', lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+c
+      ALLOCATE( champ(imdep*jmdep) )
+
+      DO  200 l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+         print*,dimfirst,dimlast
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF 
+   
+        title = 'Rugosite Amip '
+c
+        CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                      dlon, dlat, champ, interbar          )
+
+       IF ( interbar )   THEN
+         DO j = 1, imdep * jmdep
+           champ(j) = LOG(champ(j))
+         ENDDO
+
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour la rugosite $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+        CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
+     ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
+         DO j=1,jjp1
+          DO i=1,iim
+           champint(i,j)=EXP(champint(i,j))
+          ENDDO
+         ENDDO
+
+         DO j = 1, jjp1
+           DO i = 1, iim
+             IF(NINT(mask(i,j)).NE.1)  THEN
+               champint( i,j ) = 0.001
+             ENDIF
+           ENDDO
+         ENDDO
+      ELSE
+         CALL rugosite(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint, mask)
+      ENDIF
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+200      CONTINUE
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+
+      PRINT 222, timeyear
+222   FORMAT(2x,' Time year ',10f6.1)
+c
+        
+      PRINT*, 'Interpolation temporelle dans l annee'
+
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1,j,k) = champan(1,j,k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Rugosite au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+      DO k = 1, 360
+         CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+c
+c
+C Traitement de la glace oceanique
+c
+      PRINT*, 'Traitement de la glace oceanique'
+
+      ierr = NF_OPEN('amipbc_sic_1x1.nc', NF_NOWRITE, ncid)
+      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
+
+cIM22/02/2002
+cIM07/03/2002 AMIP.nc & amip79to95.nc
+cIM   ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
+cIM07/03/2002 amipbc_sic_1x1_clim.nc & amipbc_sic_1x1.nc
+      ierr = NF_INQ_VARID(ncid,'sicbcs',varid)
+cIM22/02/2002
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr),'sicbcs'
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlon_ini(imdep) )
+      ALLOCATE (     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlat_ini(jmdep) )
+      ALLOCATE (     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, lmdep
+cIM28/02/2002
+cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
+c               Ici on suppose qu'on a 12 mois (de 30 jours).
+      IF (lmdep.NE.12) THEN
+          print *, 'Unknown AMIP file: not 12 months ?'
+          STOP
+       ENDIF
+cIM28/02/2002
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+c
+      ALLOCATE ( champ(imdep*jmdep) )
+
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+ 
+         title = 'Sea-ice Amip '
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                        dlon, dlat, champ, interbar          )
+c
+      IF( oldice )  THEN
+                 CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint )
+      ELSEIF ( interbar )  THEN
+       IF( l.EQ.1 )  THEN
+        WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour Sea-ice Amip  $$$ '
+        WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+       ENDIF
+cIM07/03/2002 
+cIM22/02/2002 : Sea-ice Amip entre 0. et 1.
+cIM    PRINT*,'SUB. limit_netcdf.F IM : Sea-ice et SST Amip_new clim' 
+cIM   DO j = 1, imdep * jmdep
+cIM28/02/2002 <==PM         champ(j) = champ(j)/100.
+cIM14/03/2002      champ(j) = max(0.0,(min(1.0, (champ(j)/100.) )))
+cIM      champ(j) = amax1(0.0,(amin1(1.0, (champ(j)/100.) )))
+cIM   ENDDO
+cIM22/02/2002
+         CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
+     .             iim, jjp1, rlonv, rlatu, champint )
+      ENDIF
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+cIM28/02/2002 <== PM  timeyear(l) = timecoord(l)
+cIM      timeyear(l) = timecoord(l)
+cIM07/03/2002      
+         timeyear(l) = tmidmonth(l)
+      ENDDO
+      PRINT 222,  timeyear
+c
+      PRINT*, 'Interpolation temporelle'
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Sea ice au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+cIM14/03/2002 : Sea-ice Amip entre 0. et 1.
+      PRINT*,'SUB. limit_netcdf.F IM : Sea-ice Amipbc '
+      DO k = 1, 360
+      DO j = 1, jjp1
+      DO i = 1, iim
+        champan(i, j, k) = 
+     $ amax1(0.0,(amin1(1.0,(champan(i, j, k)/100.))))
+      ENDDO
+        champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+      ENDDO
+cIM14/03/2002
+
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon,
+     .                  champan(1,1,k), phy_ice)
+        IF ( newlmt) THEN
+
+CPB  en attendant de mettre fraction de terre
+c
+          WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
+          WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
+c 
+          IF (fracterre ) THEN
+c            WRITE(*,*) 'passe dans cas fracterre' 
+            pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
+            pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
+            pctsrf_t(1:klon,is_sic,k) =   phy_ice(1:klon) 
+     $            - pctsrf_t(1:klon,is_lic,k)
+c Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
+            WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
+              pctsrf_t(1:klon,is_sic,k) = 0.
+            END WHERE 
+            WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
+              pctsrf_t(1:klon,is_sic,k) = 0.
+              pctsrf_t(1:klon,is_oce,k) = 0.
+            END WHERE
+            DO i = 1, klon
+              IF ( 1. - zmasq(i) .GT. EPSFRA) THEN 
+                IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
+                  pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
+                  pctsrf_t(i,is_oce,k) = 0.
+                ELSE 
+                  pctsrf_t(i,is_oce,k) = 1 - zmasq(i) 
+     $                    - pctsrf_t(i,is_sic,k)
+                  IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
+                    pctsrf_t(i,is_oce,k) = 0.
+                    pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
+                  ENDIF 
+                ENDIF
+              ENDIF  
+              if (pctsrf_t(i,is_oce,k) .lt. 0.) then
+                WRITE(*,*) 'pb sous maille au point : i,k '
+     $              , i,k,pctsrf_t(:,is_oce,k)
+              ENDIF
+              IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) + 
+     $          pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k)  - 1.) 
+     $            .GT. EPSFRA) THEN 
+                  WRITE(*,*) 'physiq : pb sous surface au point ', i, 
+     $                pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
+              ENDIF 
+            END DO
+          ELSE 
+            DO i = 1, klon
+              pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
+              IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
+                pctsrf_t(i,is_sic,k) = 0.
+                pctsrf_t(i,is_oce,k) = 0.                  
+                IF(phy_ice(i) .GE. 1.e-5) THEN
+                  pctsrf_t(i,is_lic,k) = phy_ice(i)
+                  pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k) 
+     .                                   - pctsrf_t(i,is_lic,k)
+                ELSE
+                  pctsrf_t(i,is_lic,k) = 0.
+                ENDIF 
+              ELSE
+                pctsrf_t(i,is_lic,k) = 0. 
+                IF(phy_ice(i) .GE. 1.e-5) THEN 
+                  pctsrf_t(i,is_ter,k) = 0.
+                  pctsrf_t(i,is_sic,k) = phy_ice(i)
+                  pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
+                ELSE
+                  pctsrf_t(i,is_sic,k) = 0.
+                  pctsrf_t(i,is_oce,k) = 1.                      
+                ENDIF 
+              ENDIF
+              verif = pctsrf_t(i,is_ter,k) +
+     .                pctsrf_t(i,is_oce,k) + 
+     .                pctsrf_t(i,is_sic,k) +
+     .                pctsrf_t(i,is_lic,k)
+              IF ( verif .LT. 1. - 1.e-5 .OR. 
+     $             verif .GT. 1 + 1.e-5) THEN  
+                WRITE(*,*) 'pb sous maille au point : i,k,verif '
+     $                    , i,k,verif
+              ENDIF 
+            END DO
+          ENDIF 
+        ELSE  
+          DO i = 1, klon
+            phy_nat(i,k) = phy_nat0(i)
+            IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN
+              IF (NINT(phy_nat0(i)).EQ.0) THEN
+                phy_nat(i,k) = 3.0
+              ELSE
+                phy_nat(i,k) = 2.0
+              ENDIF
+            ENDIF
+            IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
+              IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
+            ENDIF
+          END DO
+        ENDIF
+      ENDDO
+c
+
+      ierr = NF_CLOSE(ncid)
+c
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+
+477    continue
+c
+C Traitement de la sst
+c
+      PRINT*, 'Traitement de la sst'
+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
+        ierr = NF_OPEN('amipbc_sst_1x1_clim.nc', NF_NOWRITE, ncid)
+        if (ierr.ne.0) THEN
+          print *, NF_STRERROR(ierr)
+          STOP
+        endif
+      ENDIF
+
+cIM22/02/2002
+cIM   ierr = NF_INQ_VARID(ncid,'SST',varid)
+      ierr = NF_INQ_VARID(ncid,'tosbcs',varid)
+cIM22/02/2002
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable SST ', namedim,'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+ 
+      ALLOCATE( dlon_ini(imdep) )
+      ALLOCATE(     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable SST ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE( dlat_ini(jmdep) )
+      ALLOCATE(     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', lmdep
+cIM28/02/2002
+cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
+c               Ici on suppose qu'on a 12 mois (de 30 jours).
+      IF (lmdep.NE.12) THEN
+          print *, 'Unknown AMIP file: not 12 months ?'
+          STOP
+       ENDIF
+cIM28/02/2002
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+       ALLOCATE( champ(imdep*jmdep) )
+       IF( extrap )   THEN
+         ALLOCATE ( work(imdep,jmdep) )
+       ENDIF
+c
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+
+         title='Sst Amip'
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                            dlon, dlat, champ, interbar     )
+       IF ( extrap )  THEN
+        CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work)
+       ENDIF
+c
+
+      IF ( interbar )  THEN
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour la Sst Amip $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
+     .          iim, jjp1, rlonv, rlatu, champint   )
+      ENDIF
+
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i,j,l) = champint(i,j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+cIM28/02/2002 <==PM  timeyear(l) = timecoord(l)
+         timeyear(l) = tmidmonth(l)
+      ENDDO
+      print 222,  timeyear
+c
+C interpolation temporelle
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i,j,l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1,j,k) = champan(1,j,k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' SST au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+cIM14/03/2002 : SST amipbc greater then 271.38
+      PRINT*,'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 '
+      DO k = 1, 360
+      DO j = 1, jjp1
+      DO i = 1, iim
+         champan(i, j, k) = amax1(champan(i, j, k), 271.38)
+      ENDDO
+         champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+      ENDDO
+cIM14/03/2002
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon, 
+     .                  champan(1,1,k), phy_sst(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+c
+       DEALLOCATE( dlon      )
+       DEALLOCATE( dlon_ini  )
+       DEALLOCATE( dlat      )
+       DEALLOCATE( dlat_ini  )
+       DEALLOCATE( champ     )
+c
+C Traitement de l'albedo
+c
+      PRINT*, 'Traitement de l albedo'
+      ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF
+      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', imdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlon_ini(imdep) )
+      ALLOCATE (     dlon(imdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', jmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+
+      ALLOCATE ( dlat_ini(jmdep) )
+      ALLOCATE (     dlat(jmdep) )
+
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+      print*,'variable ', namedim, 'dimension ', lmdep
+      ierr = NF_INQ_VARID(ncid,namedim,dimid)
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
+#else
+      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
+#endif
+      if (ierr.ne.0) then
+        print *, NF_STRERROR(ierr)
+        STOP
+      ENDIF 
+c
+      ALLOCATE ( champ(imdep*jmdep) )
+
+      DO l = 1, lmdep
+         dimfirst(1) = 1
+         dimfirst(2) = 1
+         dimfirst(3) = l
+c
+         dimlast(1) = imdep
+         dimlast(2) = jmdep
+         dimlast(3) = 1
+c
+         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
+#else
+         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
+#endif
+         if (ierr.ne.0) then
+           print *, NF_STRERROR(ierr)
+           STOP
+         ENDIF
+
+         title='Albedo Amip'
+c
+         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
+     .                            dlon, dlat, champ, interbar      )
+c
+c
+      IF ( interbar )  THEN
+        IF( l.EQ.1 )  THEN
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
+     , ' pour l Albedo Amip $$$ '
+         WRITE(6,*) '-------------------------------------------------',
+     ,'------------------------'
+        ENDIF
+
+       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
+     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
+      ELSE
+       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
+     .          iim, jjp1, rlonv, rlatu, champint   )
+      ENDIF
+c
+         DO j = 1,jjp1
+         DO i = 1, iim
+            champtime (i, j, l) = champint(i, j)
+         ENDDO
+         ENDDO
+      ENDDO
+c
+      DO l = 1, lmdep
+         timeyear(l) = timecoord(l)
+      ENDDO
+      print 222,  timeyear
+c
+C interpolation temporelle
+      DO j = 1, jjp1
+      DO i = 1, iim
+          DO l = 1, lmdep
+            ax(l) = timeyear(l)
+            ay(l) = champtime (i, j, l)
+          ENDDO
+          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
+          DO k = 1, 360
+            time = FLOAT(k-1)
+            CALL SPLINT(ax,ay,yder,lmdep,time,by)
+            champan(i,j,k) = by
+          ENDDO
+      ENDDO
+      ENDDO
+      DO k = 1, 360
+      DO j = 1, jjp1
+         champan(iip1, j, k) = champan(1, j, k)
+      ENDDO
+        IF ( k.EQ.10 )  THEN
+          DO j = 1, jjp1
+            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
+            PRINT *,' Albedo au temps 10 ', chmin,chmax,j
+          ENDDO
+        ENDIF
+      ENDDO
+c
+      DO k = 1, 360
+         CALL gr_dyn_fi(1, iip1, jjp1, klon,
+     .                  champan(1,1,k), phy_alb(1,k))
+      ENDDO
+c
+      ierr = NF_CLOSE(ncid)
+c
+c
+      DO k = 1, 360
+      DO i = 1, klon
+         phy_bil(i,k) = 0.0
+      ENDDO
+      ENDDO
+c
+      PRINT*, 'Ecriture du fichier limit'
+c
+      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
+c
+      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
+     .                       "Fichier conditions aux limites")
+      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
+      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
+c
+      dims(1) = ndim
+      dims(2) = ntim
+c
+      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
+      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
+     .                        "Jour dans l annee")
+      IF (newlmt) THEN
+c
+        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
+     .                      "Fraction ocean")
+c
+        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
+     .                      "Fraction glace de mer")
+c
+        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
+     .                      "Fraction terre")
+c
+        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
+     .                      "Fraction land ice")
+c
+      ELSE 
+        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
+        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
+     .                      "Nature du sol (0,1,2,3)")
+      ENDIF 
+      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
+      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
+     .                      "Temperature superficielle de la mer")
+      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
+      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
+     .                        "Reference flux de chaleur au sol")
+      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
+      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
+     .                        "Albedo a la surface")
+      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
+      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
+     .                        "Rugosite")
+c
+      ierr = NF_ENDDEF(nid)
+c
+      DO k = 1, 360
+c
+      debut(1) = 1
+      debut(2) = k
+      epais(1) = klon
+      epais(2) = 1
+c
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
+c
+      IF (newlmt ) THEN
+          ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
+     $        ,pctsrf_t(1,is_oce,k))
+          ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
+     $        ,pctsrf_t(1,is_sic,k))
+          ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
+     $        ,pctsrf_t(1,is_ter,k))
+          ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
+     $        ,pctsrf_t(1,is_lic,k))
+      ELSE 
+          ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
+     $        ,phy_nat(1,k))
+      ENDIF 
+c
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
+#else
+      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
+      IF (newlmt ) THEN
+          ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
+     $        ,pctsrf_t(1,is_oce,k))
+          ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
+     $        ,pctsrf_t(1,is_sic,k))
+          ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
+     $        ,pctsrf_t(1,is_ter,k))
+          ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
+     $        ,pctsrf_t(1,is_lic,k))
+      ELSE 
+          ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
+     $        ,phy_nat(1,k))
+      ENDIF 
+      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
+      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
+#endif
+c
+      ENDDO
+c
+      ierr = NF_CLOSE(nid)
+c
+      STOP
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/logic.h
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/logic.h	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/logic.h	(revision 1000)
@@ -2,18 +2,18 @@
 ! $Header$
 !
-c
-c
-c-----------------------------------------------------------------------
-c INCLUDE 'logic.h'
+!
+!
+!-----------------------------------------------------------------------
+! INCLUDE 'logic.h'
 
-      COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,
-     .  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus
-     .  ,read_start,ok_guide
+      COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
+     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
+     &  ,read_start,ok_guide,ok_strato
 
-      LOGICAL purmats,forward,leapf,apphys,statcl,conser,
-     . apdiss,apdelq,saison,ecripar,fxyhypb,ysinus
-     .  ,read_start,ok_guide
+      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
+     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
+     &  ,read_start,ok_guide,ok_strato
 
       INTEGER iflag_phys
-c$OMP THREADPRIVATE(/logic/)
-c-----------------------------------------------------------------------
+!$OMP THREADPRIVATE(/logic/)
+!-----------------------------------------------------------------------
Index: /LMDZ4/trunk/libf/dyn3dpar/mod_const_para.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/mod_const_para.F90	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/mod_const_para.F90	(revision 1000)
@@ -11,5 +11,7 @@
 
     IMPLICIT NONE
+#ifdef CPP_MPI
     INCLUDE 'mpif.h'
+#endif
     INTEGER             :: ierr
     INTEGER             :: comp_id
@@ -30,19 +32,39 @@
        CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
 !$OMP END MASTER
+#ifdef CPP_MPI
+      COMM_LMDZ=MPI_COMM_WORLD
+      MPI_REAL_LMDZ=MPI_REAL8
+#endif
 #endif
     ELSE
+      CALL init_mpi
+    ENDIF
+
+  END SUBROUTINE Init_const_mpi
+  
+  SUBROUTINE Init_mpi
+  IMPLICIT NONE
+#ifdef CPP_MPI
+     INCLUDE 'mpif.h'
+#endif
+    INTEGER             :: ierr
+    INTEGER             :: thread_required
+    INTEGER             :: thread_provided
+
+#ifdef CPP_MPI
 !$OMP MASTER
-       thread_required=MPI_THREAD_SERIALIZED
-       CALL MPI_INIT_THREAD(thread_required,thread_provided,ierr)
-       IF (thread_provided < thread_required) THEN
-         CALL abort_gcm('The multithreaded level of MPI librairy do not provide the requiered level', &
-                        'mod_const_mpi::Init_const_mpi',1)
-       ENDIF
-       COMM_LMDZ=MPI_COMM_WORLD
+      thread_required=MPI_THREAD_SERIALIZED
+
+      CALL MPI_INIT_THREAD(thread_required,thread_provided,ierr)
+      IF (thread_provided < thread_required) THEN
+        PRINT *,'Warning : The multithreaded level of MPI librairy do not provide the requiered level',  &
+                ' in mod_const_mpi::Init_const_mpi'
+      ENDIF
+      COMM_LMDZ=MPI_COMM_WORLD
+      MPI_REAL_LMDZ=MPI_REAL8
 !$OMP END MASTER
-    END IF
+#endif
 
-    MPI_REAL_LMDZ=MPI_REAL8
-  END SUBROUTINE Init_const_mpi
-
+   END SUBROUTINE Init_mpi
+    
 END MODULE mod_const_mpi
Index: /LMDZ4/trunk/libf/dyn3dpar/mod_hallo.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/mod_hallo.F90	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/mod_hallo.F90	(revision 1000)
@@ -50,5 +50,5 @@
     MaxBufferSize_Used=0
 
-    IF (use_mpi_alloc) THEN
+    IF (use_mpi_alloc .AND. using_mpi) THEN
       CALL create_global_mpi_buffer
     ELSE 
@@ -67,8 +67,14 @@
   SUBROUTINE create_global_mpi_buffer
   IMPLICIT NONE
-  INCLUDE 'mpif.h'  
+#ifdef CPP_MPI
+  INCLUDE 'mpif.h'
+#endif  
     POINTER (Pbuffer,MPI_Buffer(MaxBufferSize))
     REAL :: MPI_Buffer
+#ifdef CPP_MPI
     INTEGER(KIND=MPI_ADDRESS_KIND) :: BS 
+#else
+    INTEGER(KIND=8) :: BS
+#endif
     INTEGER :: i,ierr
 
@@ -76,5 +82,7 @@
       Bs=8*MaxBufferSize
 !$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
       CALL MPI_ALLOC_MEM(BS,MPI_INFO_NULL,Pbuffer,ierr)
+#endif
 !$OMP END CRITICAL (MPI)
       DO i=1,MaxBufferSize
@@ -164,5 +172,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
       INTEGER :: ij,ll,offset,size,target
@@ -186,5 +193,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
       INTEGER :: ij,ll,offset,size,target
@@ -211,5 +217,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
     INTEGER :: ij,ll
@@ -260,5 +265,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
     INTEGER :: ij,ll,Up,Down
@@ -313,6 +317,7 @@
 #include "dimensions.h"
 #include "paramet.h"    
+#ifdef CPP_MPI
     include 'mpif.h'
-    
+#endif    
       INTEGER :: ij,ll
       REAL, dimension(ij,ll) :: Field
@@ -379,5 +384,7 @@
 #include "dimensions.h"
 #include "paramet.h"
+#ifdef CPP_MPI
       include 'mpif.h'
+#endif
 
       type(request),target :: a_request
@@ -428,6 +435,13 @@
     
 !$OMP CRITICAL (MPI)
-         call MPI_ISSEND(Buffer(req%Pos),SizeBuffer,MPI_REAL8,rank,a_request%tag+1000*omp_rank,     &
+         
+#ifdef CPP_MPI
+         call MPI_ISSEND(Buffer(req%Pos),SizeBuffer,MPI_REAL_LMDZ,rank,a_request%tag+1000*omp_rank,     &
                          COMM_LMDZ,Req%MSG_Request,ierr)
+#endif
+         IF (.NOT.using_mpi) THEN
+           PRINT *,'Erreur, echange MPI en mode sequentiel !!!'
+           STOP
+         ENDIF
 !         PRINT *,"-------------------------------------------------------------------"
 !         PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->"
@@ -460,6 +474,14 @@
              call allocate_buffer(SizeBuffer,Req%Index,Req%Pos)
 !$OMP CRITICAL (MPI)
-             call MPI_IRECV(Buffer(Req%Pos),SizeBuffer,MPI_REAL8,rank,a_request%tag+1000*omp_rank,     &
+
+#ifdef CPP_MPI
+             call MPI_IRECV(Buffer(Req%Pos),SizeBuffer,MPI_REAL_LMDZ,rank,a_request%tag+1000*omp_rank,     &
                            COMM_LMDZ,Req%MSG_Request,ierr)
+#endif             
+             IF (.NOT.using_mpi) THEN
+               PRINT *,'Erreur, echange MPI en mode sequentiel !!!'
+               STOP
+             ENDIF
+
 !         PRINT *,"-------------------------------------------------------------------"
 !         PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->"
@@ -480,5 +502,7 @@
 #include "dimensions.h"
 #include "paramet.h"
+#ifdef CPP_MPI
       include 'mpif.h'   
+#endif
       
       type(request),target :: a_request
@@ -486,5 +510,9 @@
       type(Hallo),pointer :: PtrHallo
       integer, dimension(2*mpi_size) :: TabRequest
+#ifdef CPP_MPI
       integer, dimension(MPI_STATUS_SIZE,2*mpi_size) :: TabStatus
+#else
+      integer, dimension(1,2*mpi_size) :: TabStatus
+#endif
       integer :: NbRequest
       integer :: i,rank,pos,ij,l,ierr
@@ -515,5 +543,7 @@
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"en attente"
 !        PRINT *,"No des requetes :",TabRequest(1:NbRequest)
+#ifdef CPP_MPI
         call MPI_WAITALL(NbRequest,TabRequest,TabStatus,ierr)
+#endif
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"complete"
 !        PRINT *,"-------------------------------------------------------------------"
@@ -567,11 +597,16 @@
 #include "dimensions.h"
 #include "paramet.h"
+#ifdef CPP_MPI
       include 'mpif.h'   
-      
+#endif      
       type(request),target :: a_request
       type(request_SR),pointer :: Req
       type(Hallo),pointer :: PtrHallo
       integer, dimension(mpi_size) :: TabRequest
+#ifdef CPP_MPI
       integer, dimension(MPI_STATUS_SIZE,mpi_size) :: TabStatus
+#else
+      integer, dimension(1,mpi_size) :: TabStatus
+#endif
       integer :: NbRequest
       integer :: i,rank,pos,ij,l,ierr
@@ -594,5 +629,7 @@
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"en attente"
 !        PRINT *,"No des requetes :",TabRequest(1:NbRequest)
+#ifdef CPP_MPI
         call MPI_WAITALL(NbRequest,TabRequest,TabStatus,ierr)
+#endif
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"complete"
 !        PRINT *,"-------------------------------------------------------------------"
@@ -617,5 +654,7 @@
 #include "dimensions.h"
 #include "paramet.h"
+#ifdef CPP_MPI
       include 'mpif.h'   
+#endif
       
       type(request),target :: a_request
@@ -623,5 +662,9 @@
       type(Hallo),pointer :: PtrHallo
       integer, dimension(mpi_size) :: TabRequest
+#ifdef CPP_MPI
       integer, dimension(MPI_STATUS_SIZE,mpi_size) :: TabStatus
+#else
+      integer, dimension(1,mpi_size) :: TabStatus
+#endif
       integer :: NbRequest
       integer :: i,rank,pos,ij,l,ierr
@@ -645,5 +688,7 @@
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"en attente"
 !        PRINT *,"No des requetes :",TabRequest(1:NbRequest)
+#ifdef CPP_MPI
         call MPI_WAITALL(NbRequest,TabRequest,TabStatus,ierr)
+#endif
 !        PRINT *,"Process de rang",mpi_rank,"Task : ",omp_rank,"--->",NbRequest,"complete"
 !        PRINT *,"-------------------------------------------------------------------"
@@ -691,5 +736,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
     INTEGER :: ij,ll,l
@@ -731,5 +775,4 @@
 #include "dimensions.h"
 #include "paramet.h"    
-    include 'mpif.h'
     
     INTEGER :: ij,ll,Up,Down
Index: /LMDZ4/trunk/libf/dyn3dpar/parallel.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/parallel.F90	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/parallel.F90	(revision 1000)
@@ -1,5 +1,8 @@
   module parallel
   USE mod_const_mpi
-   
+    
+    LOGICAL,SAVE :: using_mpi
+    LOGICAL,SAVE :: using_omp
+    
     integer, save :: mpi_size
     integer, save :: mpi_rank
@@ -25,4 +28,9 @@
     USE vampir
     implicit none
+#ifdef CPP_MPI
+      include 'mpif.h'
+#endif
+#include "dimensions.h"
+#include "paramet.h"
     
       integer :: ierr
@@ -31,5 +39,6 @@
       integer, dimension(3) :: blocklen,type
       integer :: comp_id
-#ifdef _OPENMP    
+
+#ifdef CPP_OMP    
       INTEGER :: OMP_GET_NUM_THREADS
       EXTERNAL OMP_GET_NUM_THREADS
@@ -37,11 +46,22 @@
       EXTERNAL OMP_GET_THREAD_NUM
 #endif  
-      include 'mpif.h'
-#include "dimensions.h"
-#include "paramet.h"
-
+
+#ifdef CPP_MPI
+       using_mpi=.TRUE.
+#else
+       using_mpi=.FALSE.
+#endif
+      
       call InitVampir
-      call MPI_COMM_SIZE(COMM_LMDZ,mpi_size,ierr)
-      call MPI_COMM_RANK(COMM_LMDZ,mpi_rank,ierr)
+      
+      IF (using_mpi) THEN
+#ifdef CPP_MPI
+        call MPI_COMM_SIZE(COMM_LMDZ,mpi_size,ierr)
+        call MPI_COMM_RANK(COMM_LMDZ,mpi_rank,ierr)
+#endif
+      ELSE
+        mpi_size=1
+        mpi_rank=0
+      ENDIF
   
       
@@ -57,8 +77,9 @@
           
          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"
+         print *," ---> diminuez le nombre de CPU ou augmentez la taille en lattitude"
           
-          call MPI_ABORT(COMM_LMDZ,-1, ierr)
-          
+#ifdef CPP_MPI
+          IF (using_mpi) call MPI_ABORT(COMM_LMDZ,-1, ierr)
+#endif          
         endif
         
@@ -106,5 +127,5 @@
 !$OMP PARALLEL
 
-#ifdef _OPENMP
+#ifdef CPP_OMP
 !$OMP MASTER
         omp_size=OMP_GET_NUM_THREADS()
@@ -162,8 +183,10 @@
       include "dimensions.h"
       include "paramet.h"
+#ifdef CPP_MPI
+      include 'mpif.h'
+#endif      
+
       integer :: ierr
       integer :: i
-      include 'mpif.h'
-      
       deallocate(jj_begin_para)
       deallocate(jj_end_para)
@@ -178,5 +201,7 @@
 #endif 
       else
-         call MPI_FINALIZE(ierr)
+#ifdef CPP_MPI
+         IF (using_mpi) call MPI_FINALIZE(ierr)
+#endif
       end if
       
@@ -229,4 +254,21 @@
       
     end subroutine UnPack_data
+
+    
+    SUBROUTINE barrier
+    IMPLICIT NONE
+#ifdef CPP_MPI
+    INCLUDE 'mpif.h'
+#endif
+    INTEGER :: ierr
+    
+!$OMP CRITICAL (MPI)      
+#ifdef CPP_MPI
+      IF (using_mpi) CALL MPI_Barrier(COMM_LMDZ,ierr)
+#endif
+!$OMP END CRITICAL (MPI)
+    
+    END SUBROUTINE barrier
+       
       
     subroutine exchange_hallo(Field,ij,ll,up,down)
@@ -235,6 +277,7 @@
 #include "dimensions.h"
 #include "paramet.h"    
+#ifdef CPP_MPI
     include 'mpif.h'
-    
+#endif    
       INTEGER :: ij,ll
       REAL, dimension(ij,ll) :: Field
@@ -245,5 +288,9 @@
       LOGICAL :: RecvUp,RecvDown
       INTEGER, DIMENSION(4) :: Request
+#ifdef CPP_MPI
       INTEGER, DIMENSION(MPI_STATUS_SIZE,4) :: Status
+#else
+      INTEGER, DIMENSION(1,4) :: Status
+#endif
       INTEGER :: NbRequest
       REAL, dimension(:),allocatable :: Buffer_Send_up,Buffer_Send_down
@@ -251,96 +298,106 @@
       INTEGER :: Buffer_size      
 
-!$OMP CRITICAL (MPI)      
-      call MPI_Barrier(COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
-      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)
-!$OMP CRITICAL (MPI)
-        call MPI_ISSEND(Buffer_send_up,Buffer_Size,MPI_REAL8,MPI_Rank-1,1,     &
-                        COMM_LMDZ,Request(NbRequest),ierr)
-!$OMP END CRITICAL (MPI)
-      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)
-        
-!$OMP CRITICAL (MPI)
-        call MPI_ISSEND(Buffer_send_down,Buffer_Size,MPI_REAL8,MPI_Rank+1,1,     &
-                        COMM_LMDZ,Request(NbRequest),ierr)
-!$OMP END CRITICAL (MPI)
-      ENDIF
-    
-  
-      IF (RecvUp) THEN
-        NbRequest=NbRequest+1
-        buffer_size=up*iip1*ll
-        allocate(Buffer_recv_up(Buffer_size))
+      IF (using_mpi) THEN
+
+        CALL barrier
+      
+        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)
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
+          call MPI_ISSEND(Buffer_send_up,Buffer_Size,MPI_REAL8,MPI_Rank-1,1,     &
+                          COMM_LMDZ,Request(NbRequest),ierr)
+#endif
+!$OMP END CRITICAL (MPI)
+        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)
+        
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
+          call MPI_ISSEND(Buffer_send_down,Buffer_Size,MPI_REAL8,MPI_Rank+1,1,     &
+                          COMM_LMDZ,Request(NbRequest),ierr)
+#endif
+!$OMP END CRITICAL (MPI)
+        ENDIF
+    
+  
+        IF (RecvUp) THEN
+          NbRequest=NbRequest+1
+          buffer_size=up*iip1*ll
+          allocate(Buffer_recv_up(Buffer_size))
               
 !$OMP CRITICAL (MPI)
-        call MPI_IRECV(Buffer_recv_up,Buffer_size,MPI_REAL8,MPI_Rank-1,1,  &
-                        COMM_LMDZ,Request(NbRequest),ierr)
+#ifdef CPP_MPI
+          call MPI_IRECV(Buffer_recv_up,Buffer_size,MPI_REAL8,MPI_Rank-1,1,  &
+                          COMM_LMDZ,Request(NbRequest),ierr)
+#endif
 !$OMP END CRITICAL (MPI)
      
        
-      ENDIF
-  
-      IF (RecvDown) THEN
-        NbRequest=NbRequest+1
-        buffer_size=down*iip1*ll
-        allocate(Buffer_recv_down(Buffer_size))
-        
-!$OMP CRITICAL (MPI)
-        call MPI_IRECV(Buffer_recv_down,Buffer_size,MPI_REAL8,MPI_Rank+1,1,     &
-                        COMM_LMDZ,Request(NbRequest),ierr)
-!$OMP END CRITICAL (MPI)
-      
-        
-      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)
-!$OMP CRITICAL (MPI)
-      call MPI_Barrier(COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
-
+        ENDIF
+  
+        IF (RecvDown) THEN
+          NbRequest=NbRequest+1
+          buffer_size=down*iip1*ll
+          allocate(Buffer_recv_down(Buffer_size))
+        
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
+          call MPI_IRECV(Buffer_recv_down,Buffer_size,MPI_REAL8,MPI_Rank+1,1,     &
+                          COMM_LMDZ,Request(NbRequest),ierr)
+#endif
+!$OMP END CRITICAL (MPI)
+        
+        ENDIF
+  
+#ifdef CPP_MPI
+        if (NbRequest > 0) call MPI_WAITALL(NbRequest,Request,Status,ierr)
+#endif
+        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 barrier
+      
+      ENDIF  ! using_mpi
+      
       RETURN
       
@@ -352,6 +409,7 @@
 #include "dimensions.h"
 #include "paramet.h"    
+#ifdef CPP_MPI
     include 'mpif.h'
-    
+#endif    
       INTEGER :: ij,ll,rank
       REAL, dimension(ij,ll) :: Field
@@ -362,60 +420,64 @@
       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))
+      IF (using_mpi) THEN
+
+        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))
 
 !CDIR NOVECTOR
-        do i=0,MPI_Size-1
-           
+          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
+  
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
+        call MPI_GATHERV(Buffer_send,(min(ij_end,ij)-ij_begin+1)*ll,MPI_REAL8,   &
+                          Buffer_Recv,Recv_count,displ,MPI_REAL8,rank,COMM_LMDZ,ierr)
+#endif
+!$OMP END CRITICAL (MPI)
+      
+        if (MPI_Rank==rank) then                  
+      
           if (ij==ip1jmp1) then 
-            Recv_count(i)=(jj_end_para(i)-jj_begin_para(i)+1)*ll*iip1
+            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
-            Recv_count(i)=(min(jj_end_para(i),jjm)-jj_begin_para(i)+1)*ll*iip1
-          else
-            stop 'erreur dans Gather_Field'
+            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
-                 
-          if (i==0) then 
-            displ(i)=0 
-          else
-            displ(i)=displ(i-1)+Recv_count(i-1)
-          endif
-          
-        enddo
-        
-      endif
-  
-!$OMP CRITICAL (MPI)
-      call MPI_GATHERV(Buffer_send,(min(ij_end,ij)-ij_begin+1)*ll,MPI_REAL8,   &
-                        Buffer_Recv,Recv_count,displ,MPI_REAL8,rank,COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
-      
-      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 
+      ENDIF ! using_mpi
+      
     end subroutine Gather_Field
 
@@ -425,14 +487,19 @@
 #include "dimensions.h"
 #include "paramet.h"    
+#ifdef CPP_MPI
     include 'mpif.h'
-    
+#endif    
       INTEGER :: ij,ll
       REAL, dimension(ij,ll) :: Field
       INTEGER :: ierr
       
-      call Gather_Field(Field,ij,ll,0)
-!$OMP CRITICAL (MPI)
+      IF (using_mpi) THEN
+        call Gather_Field(Field,ij,ll,0)
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
       call MPI_BCAST(Field,ij*ll,MPI_REAL8,0,COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
+#endif
+!$OMP END CRITICAL (MPI)
+      ENDIF
       
     end subroutine AllGather_Field
@@ -442,6 +509,7 @@
 #include "dimensions.h"
 #include "paramet.h"    
+#ifdef CPP_MPI
     include 'mpif.h'
-    
+#endif    
       INTEGER :: ij,ll
       REAL, dimension(ij,ll) :: Field
@@ -449,8 +517,13 @@
       INTEGER :: ierr
       
-!$OMP CRITICAL (MPI)
+      IF (using_mpi) THEN
+      
+!$OMP CRITICAL (MPI)
+#ifdef CPP_MPI
       call MPI_BCAST(Field,ij*ll,MPI_REAL8,rank,COMM_LMDZ,ierr)
-!$OMP END CRITICAL (MPI)
-      
+#endif
+!$OMP END CRITICAL (MPI)
+      
+      ENDIF
     end subroutine Broadcast_Field
         
Index: /LMDZ4/trunk/libf/dyn3dpar/startvar.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/startvar.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/startvar.F	(revision 1000)
@@ -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')
@@ -1051,5 +1054,5 @@
       REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
       REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
-       REAL, ALLOCATABLE :: varrr(:,:,:)
+!       REAL, ALLOCATABLE :: varrr(:,:,:)
       INTEGER, ALLOCATABLE :: lind(:)
     !
@@ -1059,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))
     !
     !
@@ -1170,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/trunk/libf/dyn3dpar/times.F90
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/times.F90	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/times.F90	(revision 1000)
@@ -138,24 +138,33 @@
     use parallel
     implicit none
+#ifdef CPP_MPI    
     include 'mpif.h'
+#endif
     integer :: ierr
     integer :: data_size
     real, allocatable,dimension(:,:) :: tmp_table
-    
-    if (AllTimer_IsActive) then
-    
-    allocate(tmp_table(max_size,nb_timer))
-    
-    data_size=max_size*nb_timer
-    
-    tmp_table(:,:)=timer_table(:,:,mpi_rank)
-    call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL8,timer_table(1,1,mpi_rank),data_size,MPI_REAL8,COMM_LMDZ,ierr)
-
-    tmp_table(:,:)=timer_table_sqr(:,:,mpi_rank)
-    call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL8,timer_table_sqr(1,1,mpi_rank),data_size,MPI_REAL8,COMM_LMDZ,ierr)
-    
-    deallocate(tmp_table)
-    
-    endif
+
+    IF (using_mpi) THEN    
+   
+      if (AllTimer_IsActive) then
+    
+    
+      allocate(tmp_table(max_size,nb_timer))
+    
+      data_size=max_size*nb_timer
+    
+      tmp_table(:,:)=timer_table(:,:,mpi_rank)
+#ifdef CPP_MPI 
+      call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL_LMDZ,timer_table(1,1,0),data_size,MPI_REAL_LMDZ,COMM_LMDZ,ierr)
+#endif
+      tmp_table(:,:)=timer_table_sqr(:,:,mpi_rank)
+#ifdef CPP_MPI
+      call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL_LMDZ,timer_table_sqr(1,1,0),data_size,MPI_REAL_LMDZ,COMM_LMDZ,ierr)
+#endif    
+      deallocate(tmp_table)
+    
+      endif
+      
+    ENDIF ! using_mpi
     
   end subroutine allgather_timer
@@ -164,5 +173,7 @@
     use parallel
     implicit none
+#ifdef CPP_MPI
     include 'mpif.h'
+#endif
     integer :: ierr
     integer :: data_size
@@ -170,24 +181,31 @@
     integer, allocatable,dimension(:,:),target :: tmp_iter
     integer :: istats
-    
-    if (AllTimer_IsActive) then
-    
-    allocate(tmp_table(max_size,nb_timer))
-    allocate(tmp_iter(max_size,nb_timer))
+
+    IF (using_mpi) THEN
+        
+      if (AllTimer_IsActive) then
+    
+      allocate(tmp_table(max_size,nb_timer))
+      allocate(tmp_iter(max_size,nb_timer))
    
-    data_size=max_size*nb_timer
-
-    tmp_table(:,:)=timer_average(:,:,mpi_rank)
-    call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL8,timer_average(1,1,0),data_size,MPI_REAL8,COMM_LMDZ,ierr)
-
-    tmp_table(:,:)=timer_delta(:,:,mpi_rank)
-    call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL8,timer_delta(1,1,0),data_size,MPI_REAL8,COMM_LMDZ,ierr)
-
-    tmp_iter(:,:)=timer_iteration(:,:,mpi_rank)
-    call mpi_allgather(tmp_iter(1,1),data_size,MPI_INTEGER,timer_iteration(1,1,0),data_size,MPI_INTEGER,COMM_LMDZ,ierr)
-    
-    deallocate(tmp_table)
-    
-    endif
+      data_size=max_size*nb_timer
+
+      tmp_table(:,:)=timer_average(:,:,mpi_rank)
+#ifdef CPP_MPI
+      call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL_LMDZ,timer_average(1,1,0),data_size,MPI_REAL_LMDZ,COMM_LMDZ,ierr)
+#endif
+      tmp_table(:,:)=timer_delta(:,:,mpi_rank)
+#ifdef CPP_MPI
+      call mpi_allgather(tmp_table(1,1),data_size,MPI_REAL_LMDZ,timer_delta(1,1,0),data_size,MPI_REAL_LMDZ,COMM_LMDZ,ierr)
+#endif
+      tmp_iter(:,:)=timer_iteration(:,:,mpi_rank)
+#ifdef CPP_MPI
+      call mpi_allgather(tmp_iter(1,1),data_size,MPI_INTEGER,timer_iteration(1,1,0),data_size,MPI_INTEGER,COMM_LMDZ,ierr)
+#endif    
+      deallocate(tmp_table)
+    
+      endif
+      
+    ENDIF  ! using_mpî
   end subroutine allgather_timer_average
   
Index: /LMDZ4/trunk/libf/dyn3dpar/top_bound_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/top_bound_p.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/top_bound_p.F	(revision 1000)
@@ -0,0 +1,113 @@
+      SUBROUTINE top_bound_p( vcov,ucov,teta, du,dv,dh )
+      USE parallel
+      IMPLICIT NONE
+c
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+CC#include "comgeom2.h"
+
+
+c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
+C     F. LOTT DEC. 2006
+c                                 (  10/12/06  )
+
+c=======================================================================
+c
+c   Auteur:  F. LOTT  
+c   -------
+c
+c   Objet:
+c   ------
+c
+c   Dissipation linéaire (ex top_bound de la physique)
+c
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Declarations:
+c   -------------
+
+#include "comgeom.h"
+#include "comdissipn.h"
+
+c   Arguments:
+c   ----------
+
+      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
+      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
+
+c   Local:
+c   ------
+
+      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
+      
+      INTEGER NDAMP
+      PARAMETER (NDAMP=4)
+      integer i	
+      REAL :: rdamp(llm) = 
+     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) 
+
+      INTEGER j,l,jjb,jje
+
+
+C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
+
+      jjb=jj_begin
+      jje=jj_end
+      IF (pole_sud) jje=jj_end-1
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      do l=1,llm
+        do j=jjb,jje
+          vzon(j,l)=0.
+          do i=1,iim
+            vzon(j,l)=vzon(j,l)+vcov(i,j,l)/float(iim)
+          enddo
+        enddo
+      enddo
+c$OMP END DO NOWAIT   
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      do l=1,llm
+        do j=jjb,jje
+          do i=1,iip1
+            dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
+          enddo
+        enddo
+      enddo
+c$OMP END DO NOWAIT
+
+      jjb=jj_begin
+      jje=jj_end
+      IF (pole_nord) jjb=jj_begin+1
+      IF (pole_sud)  jje=jj_end-1
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      do l=1,llm
+        do j=jjb,jje
+          uzon(j,l)=0.
+          tzon(j,l)=0.
+          do i=1,iim
+            uzon(j,l)=uzon(j,l)+ucov(i,j,l)/float(iim)
+            tzon(j,l)=tzon(j,l)+teta(i,j,l)/float(iim)
+          enddo
+        enddo
+      enddo
+c$OMP END DO NOWAIT
+
+C   AMORTISSEMENTS LINEAIRES:
+
+c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)      
+      do l=1,llm
+        do j=jjb,jje
+          do i=1,iip1
+            du(i,j,l)=du(i,j,l)-rdamp(l)*(ucov(i,j,l)-uzon(j,l))
+            dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
+          enddo
+       enddo
+      enddo
+c$OMP END DO NOWAIT
+      
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/tourpot.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/tourpot.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/tourpot.F	(revision 1000)
@@ -0,0 +1,81 @@
+!
+! $Header$
+!
+      SUBROUTINE tourpot ( vcov, ucov, massebxy, vorpot )
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteur:  P. Le Van
+c   -------
+c
+c   Objet:
+c   ------
+c
+c    *******************************************************************
+c    .........      calcul du tourbillon potentiel             .........
+c    *******************************************************************
+c
+c     vcov,ucov,fext et pbarxyfl sont des argum. d'entree pour le s-pg .
+c             vorpot            est  un argum.de sortie pour le s-pg .
+c
+c=======================================================================
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom.h"
+#include "logic.h"
+
+      REAL  rot( ip1jm,llm )
+      REAL vcov( ip1jm,llm ),ucov( ip1jmp1,llm )
+      REAL massebxy( ip1jm,llm ),vorpot( ip1jm,llm )
+
+      INTEGER l, ij
+
+
+
+
+c  ... vorpot = ( Filtre( d(vcov)/dx - d(ucov)/dy ) + fext ) /psbarxy ..
+
+
+
+c    ........  Calcul du rotationnel du vent V  puis filtrage  ........
+
+      DO 5 l = 1,llm
+
+      DO 2 ij = 1, ip1jm - 1
+      rot( ij,l ) = vcov(ij+1,l)-vcov(ij,l)+ucov(ij+iip1,l)-ucov(ij,l)
+   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, llm, 2, 1, .FALSE., 1 )
+
+
+      DO 10 l = 1, llm
+
+      DO 6 ij = 1, ip1jm - 1
+      vorpot( ij,l ) = ( rot(ij,l) + fext(ij) ) / massebxy(ij,l)
+   6  CONTINUE
+
+c    ..... correction pour  vorpot( iip1,j,l)  .....
+c    ....   vorpot(iip1,j,l)= vorpot(1,j,l) ....
+CDIR$ IVDEP
+      DO 8 ij = iip1, ip1jm, iip1
+      vorpot( ij,l ) = vorpot( ij -iim,l )
+   8  CONTINUE
+
+  10  CONTINUE
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/vitvert.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/vitvert.F	(revision 1000)
+++ /LMDZ4/trunk/libf/dyn3dpar/vitvert.F	(revision 1000)
@@ -0,0 +1,52 @@
+!
+! $Header$
+!
+      SUBROUTINE vitvert ( convm , w )
+c
+      IMPLICIT NONE
+
+c=======================================================================
+c
+c   Auteurs:  P. Le Van , F. Hourdin .
+c   -------
+c
+c   Objet:
+c   ------
+c
+c    *******************************************************************
+c  .... calcul de la vitesse verticale aux niveaux sigma  ....
+c    *******************************************************************
+c     convm   est un argument  d'entree pour le s-pg  ......
+c       w     est un argument de sortie pour le s-pg  ......
+c
+c    la vitesse verticale est orientee de  haut en bas .
+c    au sol, au niveau sigma(1),   w(i,j,1) = 0.
+c    au sommet, au niveau sigma(llm+1) , la vit.verticale est aussi
+c    egale a 0. et n'est pas stockee dans le tableau w  .
+c
+c
+c=======================================================================
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+
+      REAL w(ip1jmp1,llm),convm(ip1jmp1,llm)
+      INTEGER   l, ij
+
+
+
+      DO 2  l = 1,llmm1
+
+      DO 1 ij = 1,ip1jmp1
+      w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
+   1  CONTINUE
+
+   2  CONTINUE
+
+      DO 5 ij  = 1,ip1jmp1
+      w(ij,1)  = 0.
+5     CONTINUE
+
+      RETURN
+      END
Index: /LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/writedynav_p.F	(revision 1000)
@@ -105,9 +105,9 @@
 C  Vents V scalaire
 C
-      if (pole_sud) ije=jj_end-iip1
+      if (pole_sud) ije=ij_end-iip1
       if (pole_sud) jjn=jj_nb-1
       
       call gr_v_scal_p(llm, vnat, vs)
-      call histwrite(histid, 'v', itau_w, vs(ijb::ije,:), 
+      call histwrite(histid, 'v', itau_w, vs(ijb:ije,:), 
      .               iip1*jjn*llm, ndex3d)
 C
@@ -118,5 +118,5 @@
       jjn=jj_nb
       
-      call histwrite(histid, 'theta', itau_w, teta(ijb::ije,:), 
+      call histwrite(histid, 'theta', itau_w, teta(ijb:ije,:), 
      .                iip1*jjn*llm, ndex3d)
 C
Index: /LMDZ4/trunk/libf/dyn3dpar/writehist_p.F
===================================================================
--- /LMDZ4/trunk/libf/dyn3dpar/writehist_p.F	(revision 999)
+++ /LMDZ4/trunk/libf/dyn3dpar/writehist_p.F	(revision 1000)
@@ -90,5 +90,5 @@
       jjn=jj_nb
           
-      call histwrite(histid, 'ucov', itau_w, ucov, 
+      call histwrite(histid, 'ucov', itau_w, ucov(ijb:ije,:), 
      .               iip1*jjn*llm, ndexu)
 
@@ -96,5 +96,5 @@
 C  Vents V
 C
-      if (pole_sud) ije=jj_end-iip1
+      if (pole_sud) ije=ij_end-iip1
       if (pole_sud) jjn=jj_nb-1
       
