Index: trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1302)
@@ -225,10 +225,4 @@
       save unskap
 
-cIM diagnostique PVteta, Amip2
-      INTEGER,PARAMETER :: ntetaSTD=3
-      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
-      REAL PVteta(klon,ntetaSTD)
-      
-      
       REAL SSUM
 
@@ -267,6 +261,4 @@
       klon=klon_mpi
       
-      PVteta(:,:)=0.
-            
       IF ( firstcal )  THEN
         debut = .TRUE.
@@ -684,16 +676,4 @@
       endif
 
-
-      IF (is_sequential.and.(planet_type=="earth")) THEN
-#ifdef CPP_EARTH
-! PVtheta calls tetalevel, which is in the physics
-cIM calcul PV a teta=350, 380, 405K
-        CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
-     $           ztfi,zplay,zplev,
-     $           ntetaSTD,rtetaSTD,PVteta)
-c
-#endif
-      ENDIF
-
 c On change de grille, dynamique vers physiq, pour le flux de masse verticale
       CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
@@ -923,7 +903,5 @@
      .             zdqfi_omp,
      .             zdpsrf_omp,
-cIM diagnostique PVteta, Amip2          
-     .             pducov,
-     .             PVteta)
+     .             pducov)
 
       else if ( planet_type=="generic" ) then
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90	(revision 1302)
@@ -1,4 +1,4 @@
 !
-! $Id: ce0l.F90 1615 2012-02-10 15:42:26Z emillour $
+! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
 !
 !-------------------------------------------------------------------------------
@@ -115,12 +115,10 @@
   END IF
 
-  IF (grilles_gcm_netcdf) THEN
-     WRITE(lunout,'(//)')
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
-     WRITE(lunout,*) '  ***************************  '
-     WRITE(lunout,'(//)')
-     CALL grilles_gcm_netcdf_sub(masque,phis)
-  END IF
+  WRITE(lunout,'(//)')
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
+  WRITE(lunout,*) '  ***************************  '
+  WRITE(lunout,'(//)')
+  CALL grilles_gcm_netcdf_sub(masque,phis)
   
 #ifdef CPP_MPI
@@ -137,2 +135,3 @@
 !
 !-------------------------------------------------------------------------------
+
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1302)
@@ -2,8 +2,8 @@
 ! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
 !
-c
-c
+!
+!
       SUBROUTINE conf_gcm( tapedef, etatinit )
-c
+!
 #ifdef CPP_IOIPSL
       use IOIPSL
@@ -20,18 +20,18 @@
       use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
       IMPLICIT NONE
-c-----------------------------------------------------------------------
-c     Auteurs :   L. Fairhead , P. Le Van  .
-c
-c     Arguments :
-c
-c     tapedef   :
-c     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
-c     -metres  du zoom  avec  celles lues sur le fichier start .
-c
+!-----------------------------------------------------------------------
+!     Auteurs :   L. Fairhead , P. Le Van  .
+!
+!     Arguments :
+!
+!     tapedef   :
+!     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para- 
+!     -metres  du zoom  avec  celles lues sur le fichier start .
+!
        LOGICAL etatinit
        INTEGER tapedef
 
-c   Declarations :
-c   --------------
+!   Declarations :
+!   --------------
 #include "dimensions.h"
 #include "paramet.h"
@@ -45,8 +45,8 @@
 ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
 ! #include "clesphys.h"
-c
-c
-c   local:
-c   ------
+!
+!
+!   local:
+!   ------
 
       CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
@@ -60,32 +60,32 @@
       integer,external :: OMP_GET_NUM_THREADS
 #endif      
-c
-c  -------------------------------------------------------------------
-c
-c       .........     Version  du 29/04/97       ..........
-c
-c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
-c      tetatemp   ajoutes  pour la dissipation   .
-c
-c   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
-c
-c  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
-c    Sinon , choix de fxynew  , a derivee sinusoidale  ..
-c
-c   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
-c         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
-c                de limit.dat ( dic)                        ...........
-c           Sinon  etatinit = . FALSE .
-c
-c   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
-c    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
-c   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
-c    lectba .  
-c   Ces parmetres definissant entre autres la grille et doivent etre
-c   pareils et coherents , sinon il y aura  divergence du gcm .
-c
-c-----------------------------------------------------------------------
-c   initialisations:
-c   ----------------
+!
+!  -------------------------------------------------------------------
+!
+!       .........     Version  du 29/04/97       ..........
+!
+!   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
+!      tetatemp   ajoutes  pour la dissipation   .
+!
+!   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb ** 
+!
+!  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
+!    Sinon , choix de fxynew  , a derivee sinusoidale  ..
+!
+!   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
+!         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
+!                de limit.dat ( dic)                        ...........
+!           Sinon  etatinit = . FALSE .
+!
+!   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
+!    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
+!   celles passees  par run.def ,  au debut du gcm, apres l'appel a 
+!    lectba .  
+!   Ces parmetres definissant entre autres la grille et doivent etre
+!   pareils et coherents , sinon il y aura  divergence du gcm .
+!
+!-----------------------------------------------------------------------
+!   initialisations:
+!   ----------------
 
 !Config  Key  = lunout
@@ -290,6 +290,6 @@
        CALL getin('dissip_period',dissip_period)
 
-ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
-ccc
+!cc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
+!cc
 
 !Config  Key  = lstardis
@@ -456,5 +456,5 @@
        CALL getin('ok_guide',ok_guide)
 
-c    ...............................................................
+!     ...............................................................
 
 !Config  Key  =  read_start
@@ -632,12 +632,6 @@
       CALL getin('ok_etat0',ok_etat0)
 
-!Config  Key  = grilles_gcm_netcdf
-!Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
-!Config  Def  = n
-      grilles_gcm_netcdf = .FALSE.
-      CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
-
-c----------------------------------------
-c Parameters for zonal averages in the case of Titan
+!----------------------------------------
+! Parameters for zonal averages in the case of Titan
       moyzon_mu = .false.
       moyzon_ch = .false.
@@ -646,13 +640,13 @@
        CALL getin('moyzon_ch', moyzon_ch)
       endif
-c----------------------------------------
-
-c----------------------------------------
-ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
-c     .........   (  modif  le 17/04/96 )   .........
-c
-C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
-c
-c----------------------------------------
+!----------------------------------------
+
+!----------------------------------------
+!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
+!     .........   (  modif  le 17/04/96 )   .........
+!
+! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
+!
+!----------------------------------------
       IF( etatinit ) then
 
@@ -707,8 +701,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -758,8 +752,8 @@
        ysinus = .TRUE.
        CALL getin('ysinus',ysinus)
-c
-c----------------------------------------
+!
+!----------------------------------------
        else ! etatinit=false
-c----------------------------------------
+!----------------------------------------
 
 !Config  Key  = clon
@@ -779,6 +773,6 @@
        CALL getin('clat',clatt)
 
-c
-c
+!
+!
       IF( ABS(clat - clatt).GE. 0.001 )  THEN
         write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
@@ -834,8 +828,8 @@
 
       write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
-c
-c    alphax et alphay sont les anciennes formulat. des grossissements
-c
-c
+!
+!    alphax et alphay sont les anciennes formulat. des grossissements
+!
+!
 
 !Config  Key  = fxyhypb
@@ -862,5 +856,5 @@
          ENDIF
       ENDIF
-c
+!
 !Config  Key  = dzoomx
 !Config  Desc = extension en longitude
@@ -925,5 +919,5 @@
       ENDIF
 
-cc
+!c
       IF( .NOT.fxyhypb  )  THEN
 
@@ -955,5 +949,5 @@
 
       endif ! etatinit
-c----------------------------------------
+!----------------------------------------
 
 
@@ -1009,5 +1003,4 @@
       write(lunout,*)' ok_limit = ', ok_limit
       write(lunout,*)' ok_etat0 = ', ok_etat0
-      write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
       if (planet_type=="titan") then
        write(lunout,*)' moyzon_mu = ', moyzon_mu
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p.F	(revision 1300)
+++ 	(revision )
@@ -1,231 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
-c
-c     Auteurs :  P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c                                 -------- z                                   
-c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
-c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
-c    ( voir note de Fr.Hourdin )  ,
-c
-c    on determine successivement , du haut vers le bas des couches, les 
-c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
-c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
-c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
-c
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL unpl2k,dellta
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_hyb_p"
-c
-
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-
-c$OMP BARRIER
-
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-      
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
-
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-      unpl2k    = 1.+ 2.* kappa
-c
-      ijb=ij_begin
-      ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
-
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO     ij      = ijb,ije
-       alpha(ij,llm) = 0.
-       beta (ij,llm) = 1./ unpl2k
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
-c
-      DO l = llm -1 , 2 , -1
-c
-c$OMP DO SCHEDULE(STATIC)
-        DO ij = ijb, ije
-        dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
-        alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
-        beta (ij,l)  =   p(ij,l  ) / dellta   
-        ENDDO
-c$OMP ENDDO NOWAIT
-c
-      ENDDO
-
-c
-c  ***********************************************************************
-c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
-c
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
-     *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
-      ENDDO
-c$OMP ENDDO NOWAIT
-c
-c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
-c
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-c
-c
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-
-      RETURN
-      END
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 1302)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90	(revision 1302)
@@ -0,0 +1,187 @@
+module exner_hyb_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
+
+    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !                                 -------- z
+    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
+    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    !    on determine successivement , du haut vers le bas des couches, les 
+    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2), 
+    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches,  
+    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
+    !
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+    REAL alpha(ngrid,llm),beta(ngrid,llm)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij
+    REAL unpl2k,dellta
+
+    INTEGER ije,ijb,jje,jjb
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_hyb_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    !$OMP BARRIER
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
+
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    unpl2k    = 1.+ 2.* kappa
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    ijb=ij_begin
+    ije=ij_end
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
+
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO     ij      = ijb,ije
+       alpha(ij,llm) = 0.
+       beta (ij,llm) = 1./ unpl2k
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
+    !
+    DO l = llm -1 , 2 , -1
+       !
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij = ijb, ije
+          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
+          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
+          beta (ij,l)  =   p(ij,l  ) / dellta   
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !  ***********************************************************************
+    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
+    !
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
+            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
+    ENDDO
+    !$OMP ENDDO NOWAIT
+    !
+    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
+    !
+    DO l = 2, llm
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+       ENDDO
+       !$OMP ENDDO NOWAIT        
+    ENDDO
+
+    if (present(pkf)) then
+       !    calcul de pkf
+
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT             
+       ENDDO
+
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_hyb_p
+
+end module exner_hyb_p_m
+
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p.F	(revision 1300)
+++ 	(revision )
@@ -1,218 +1,0 @@
-!
-! $Id $
-!
-      SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
-c
-c     Auteurs :  F. Forget , Y. Wanherdrick
-c P.Le Van  , Fr. Hourdin  .
-c    ..........
-c
-c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
-c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
-c
-c   ************************************************************************
-c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
-c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
-c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
-c   ************************************************************************
-c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
-c    la pression et la fonction d'Exner  au  sol  .
-c
-c     WARNING : CECI est une version speciale de exner_hyb originale
-c               Utilise dans la version martienne pour pouvoir 
-c               tourner avec des coordonnees verticales complexe
-c              => Il ne verifie PAS la condition la proportionalite en 
-c              energie totale/ interne / potentielle (F.Forget 2001)
-c    ( voir note de Fr.Hourdin )  ,
-c
-      USE parallel_lmdz
-      IMPLICIT NONE
-c
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comgeom.h"
-#include "comvert.h"
-#include "serre.h"
-
-      INTEGER  ngrid
-      REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
-      REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
-
-c    .... variables locales   ...
-
-      INTEGER l, ij
-      REAL dum1
-
-      REAL ppn(iim),pps(iim)
-      REAL xpn, xps
-      REAL SSUM
-      EXTERNAL SSUM
-      INTEGER ije,ijb,jje,jjb
-      logical,save :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall) 
-      character(len=*),parameter :: modname="exner_milieu_p"
-
-      ! Sanity check
-      if (firstcall) then
-        ! sanity checks for Shallow Water case (1 vertical layer)
-        if (llm.eq.1) then
-          if (kappa.ne.1) then
-            call abort_gcm(modname,
-     &      "kappa!=1 , but running in Shallow Water mode!!",42)
-          endif
-          if (cpp.ne.r) then
-            call abort_gcm(modname,
-     &      "cpp!=r , but running in Shallow Water mode!!",42)
-          endif
-        endif ! of if (llm.eq.1)
-
-        firstcall=.false.
-      endif ! of if (firstcall)
-      
-c$OMP BARRIER
-
-! Specific behaviour for Shallow Water (1 vertical layer) case
-      if (llm.eq.1) then
-              
-        ! Compute pks(:),pk(:),pkf(:)
-        ijb=ij_begin
-        ije=ij_end
-!$OMP DO SCHEDULE(STATIC)
-        DO ij=ijb, ije
-          pks(ij)=(cpp/preff)*ps(ij)
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-!$OMP ENDDO
-
-!$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-          pk(ij,1) = .5*pks(ij)
-          pkf(ij,1)=pk(ij,1)
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
-          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
-        ENDDO
-      endif
-!$OMP END MASTER
-!$OMP BARRIER
-        jjb=jj_begin
-        jje=jj_end
-        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-
-        ! our work is done, exit routine
-        return
-      endif ! of if (llm.eq.1)
-
-!!!! General case:
-
-c     -------------
-c     Calcul de pks
-c     -------------
-   
-      ijb=ij_begin
-      ije=ij_end
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij  = ijb, ije
-        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-      ENDDO
-c$OMP ENDDO
-c Synchro OPENMP ici
-
-c$OMP MASTER
-      if (pole_nord) then
-        DO  ij   = 1, iim
-          ppn(ij) = aire(   ij   ) * pks(  ij     )
-        ENDDO
-        xpn      = SSUM(iim,ppn,1) /apoln
-  
-        DO ij   = 1, iip1
-          pks(   ij     )  =  xpn
-        ENDDO
-      endif
-      
-      if (pole_sud) then
-        DO  ij   = 1, iim
-          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
-        ENDDO
-        xps      = SSUM(iim,pps,1) /apols 
-  
-        DO ij   = 1, iip1
-          pks( ij+ip1jm )  =  xps
-        ENDDO
-      endif
-c$OMP END MASTER
-c$OMP BARRIER
-c
-c
-c    .... Calcul de pk  pour la couche l 
-c    --------------------------------------------
-c
-      dum1 = cpp * (2*preff)**(-kappa) 
-      DO l = 1, llm-1
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
-        ENDDO
-c$OMP ENDDO NOWAIT        
-      ENDDO
-
-c    .... Calcul de pk  pour la couche l = llm ..
-c    (on met la meme distance (en log pression)  entre Pk(llm)
-c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
-
-c$OMP DO SCHEDULE(STATIC)
-      DO   ij   = ijb, ije
-         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
-      ENDDO
-c$OMP ENDDO NOWAIT        
-
-
-c    calcul de pkf
-c    -------------
-c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
-      DO l = 1, llm
-c$OMP DO SCHEDULE(STATIC)
-         DO   ij   = ijb, ije
-           pkf(ij,l)=pk(ij,l)
-         ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-c$OMP BARRIER
-      
-      jjb=jj_begin
-      jje=jj_end
-      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
-      
-c    EST-CE UTILE ?? : calcul de beta
-c    --------------------------------
-      DO l = 2, llm
-c$OMP DO SCHEDULE(STATIC)
-        DO   ij   = ijb, ije
-          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
-        ENDDO
-c$OMP ENDDO NOWAIT             
-      ENDDO
-
-      RETURN
-      END
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 1302)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90	(revision 1302)
@@ -0,0 +1,160 @@
+module exner_milieu_p_m
+
+  IMPLICIT NONE
+
+contains
+
+  SUBROUTINE  exner_milieu_p ( ngrid, ps, p, pks, pk, pkf )
+    !
+    !     Auteurs :  F. Forget , Y. Wanherdrick
+    ! P.Le Van  , Fr. Hourdin  .
+    !    ..........
+    !
+    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
+    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
+    !
+    !   ************************************************************************
+    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des 
+    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
+    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
+    !   ************************************************************************
+    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
+    !    la pression et la fonction d'Exner  au  sol  .
+    !
+    !     WARNING : CECI est une version speciale de exner_hyb originale
+    !               Utilise dans la version martienne pour pouvoir 
+    !               tourner avec des coordonnees verticales complexe
+    !              => Il ne verifie PAS la condition la proportionalite en 
+    !              energie totale/ interne / potentielle (F.Forget 2001)
+    !    ( voir note de Fr.Hourdin )  ,
+    !
+    USE parallel_lmdz
+    !
+    include "dimensions.h"
+    include "paramet.h"
+    include "comconst.h"
+    include "comgeom.h"
+    include "comvert.h"
+    include "serre.h"
+
+    INTEGER  ngrid
+    REAL p(ngrid,llmp1),pk(ngrid,llm)
+    REAL, optional:: pkf(ngrid,llm)
+    REAL ps(ngrid),pks(ngrid)
+
+    !    .... variables locales   ...
+
+    INTEGER l, ij,ijb,ije,jjb,jje
+    REAL dum1
+
+    logical,save :: firstcall=.true.
+    !$OMP THREADPRIVATE(firstcall) 
+    character(len=*),parameter :: modname="exner_milieu_p"
+
+    ! Sanity check
+    if (firstcall) then
+       ! sanity checks for Shallow Water case (1 vertical layer)
+       if (llm.eq.1) then
+          if (kappa.ne.1) then
+             call abort_gcm(modname, &
+                  "kappa!=1 , but running in Shallow Water mode!!",42)
+          endif
+          if (cpp.ne.r) then
+             call abort_gcm(modname, &
+                  "cpp!=r , but running in Shallow Water mode!!",42)
+          endif
+       endif ! of if (llm.eq.1)
+
+       firstcall=.false.
+    endif ! of if (firstcall)
+
+    !$OMP BARRIER
+
+    ! Specific behaviour for Shallow Water (1 vertical layer) case:
+    if (llm.eq.1) then
+
+       ! Compute pks(:),pk(:),pkf(:)
+       ijb=ij_begin
+       ije=ij_end
+       !$OMP DO SCHEDULE(STATIC)
+       DO ij=ijb, ije
+          pks(ij) = (cpp/preff) * ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          if (present(pkf)) pkf(ij,1)=pk(ij,1)
+       ENDDO
+       !$OMP ENDDO
+
+       !$OMP BARRIER
+       if (present(pkf)) then
+          jjb=jj_begin
+          jje=jj_end
+          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+       end if
+
+       ! our work is done, exit routine
+       return
+    endif ! of if (llm.eq.1)
+
+    ! General case:
+
+    !     -------------
+    !     Calcul de pks
+    !     -------------
+
+    ijb=ij_begin
+    ije=ij_end
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij  = ijb, ije
+       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+    ENDDO
+    !$OMP ENDDO
+    ! Synchro OPENMP ici
+
+    !$OMP BARRIER
+    !
+    !
+    !    .... Calcul de pk  pour la couche l 
+    !    --------------------------------------------
+    !
+    dum1 = cpp * (2*preff)**(-kappa) 
+    DO l = 1, llm-1
+       !$OMP DO SCHEDULE(STATIC)
+       DO   ij   = ijb, ije
+          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
+       ENDDO
+       !$OMP ENDDO NOWAIT
+    ENDDO
+
+    !    .... Calcul de pk  pour la couche l = llm ..
+    !    (on met la meme distance (en log pression)  entre Pk(llm)
+    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
+
+    !$OMP DO SCHEDULE(STATIC)
+    DO   ij   = ijb, ije
+       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
+    ENDDO
+    !$OMP ENDDO NOWAIT        
+
+    if (present(pkf)) then
+       !    calcul de pkf
+
+       DO l = 1, llm
+          !$OMP DO SCHEDULE(STATIC)
+          DO   ij   = ijb, ije
+             pkf(ij,l)=pk(ij,l)
+          ENDDO
+          !$OMP ENDDO NOWAIT
+       ENDDO
+
+       !$OMP BARRIER
+
+       jjb=jj_begin
+       jje=jj_end
+       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+    end if
+
+  END SUBROUTINE exner_milieu_p
+
+end module exner_milieu_p_m
+
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F	(revision 1302)
@@ -516,5 +516,9 @@
 
 
-      day_end = day_ini + nday
+      if (nday>=0) then ! standard case
+        day_end=day_ini+nday
+      else ! special case when nday <0, run -nday dynamical steps
+        day_end=day_ini-nday/day_step
+      endif
       if (less1day) then
         day_end=day_ini+floor(time_0+fractday)
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90	(revision 1302)
@@ -328,4 +328,6 @@
 !=======================================================================
   SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
     USE parallel_lmdz
     USE control_mod
@@ -349,6 +351,5 @@
     REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
     ! Variables pour fonction Exner (P milieu couche)
-    REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-    REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+    REAL, DIMENSION (iip1,jjp1,llm)    :: pk
     REAL, DIMENSION (iip1,jjp1)        :: pks    
     REAL                               :: unskap
@@ -491,19 +492,19 @@
     IF (f_out) THEN
 !       Calcul niveaux pression milieu de couches 
-	CALL pression_p( ip1jmp1, ap, bp, ps, p )
-	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
-	else
-          CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+        CALL pression_p( ip1jmp1, ap, bp, ps, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
+        else
+          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
         endif
         unskap=1./kappa
-	DO l = 1, llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
-		    p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-		ENDDO
-	    ENDDO
-	ENDDO
-        CALL guide_out("P",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
+        DO l = 1, llm
+           DO j=jjb_u,jje_u
+              DO i =1, iip1
+                 p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+              ENDDO
+           ENDDO
+        ENDDO
+        CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
     ENDIF
     
@@ -517,5 +518,7 @@
         if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
-        IF (f_out) CALL guide_out("U",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
+        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(:,:)+tau*ugui2(:,:),factt)
+        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
         ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -529,5 +532,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
-        IF (f_out) CALL guide_out("T",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
         teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -541,5 +544,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
         CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
-        IF (f_out) CALL guide_out("SP",jjp1,1,f_add(1:ip1jmp1,1),factt)
+        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
         ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
         CALL pression_p(ip1jmp1,ap,bp,ps,p)
@@ -555,5 +558,5 @@
         if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
         CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
-        IF (f_out) CALL guide_out("Q",jjp1,llm,f_add(:,:),factt)
+        IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
         q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
     endif
@@ -568,5 +571,7 @@
         if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
         CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
-        IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:),factt)
+        IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
+        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1(:,:)+tau*vgui2(:,:),factt)
+        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
         vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
     endif
@@ -689,4 +694,6 @@
 !=======================================================================
   SUBROUTINE guide_interp(psi,teta)
+    use exner_hyb_p_m, only: exner_hyb_p
+    use exner_milieu_p_m, only: exner_milieu_p
   USE parallel_lmdz
   USE mod_hallo
@@ -713,6 +720,5 @@
   REAL, DIMENSION (iip1,jjm,llm)     :: pbary 
   ! Variables pour fonction Exner (P milieu couche)
-  REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
-  REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
+  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
   REAL, DIMENSION (iip1,jjp1)        :: pks    
   REAL                               :: unskap
@@ -784,25 +790,25 @@
     IF (guide_plevs.EQ.1) THEN
         DO l=1,llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
                     pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
-		ENDDO
-	    ENDDO
+                ENDDO
+            ENDDO
         ENDDO
     ELSE
-	CALL pression_p( ip1jmp1, ap, bp, psi, p )
-	if (pressure_exner) then
-          CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
+        CALL pression_p( ip1jmp1, ap, bp, psi, p )
+        if (pressure_exner) then
+          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
         else
-          CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
+          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
         endif
-	unskap=1./kappa
-	DO l = 1, llm
-	    DO j=jjb_u,jje_u
-		DO i =1, iip1
-		    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
-		ENDDO
-	    ENDDO
-	ENDDO
+        unskap=1./kappa
+        DO l = 1, llm
+            DO j=jjb_u,jje_u
+                DO i =1, iip1
+                    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
+                ENDDO
+            ENDDO
+        ENDDO
     ENDIF
 
@@ -1024,10 +1030,10 @@
         ! Calcul des nouvelles valeurs des niveaux de pression du guidage
         IF (guide_plevs.EQ.1) THEN
-	 CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
-	 CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
             DO l=1,nlevnc
                 DO j=jjb_v,jje_v
@@ -1041,10 +1047,10 @@
             ENDDO
         ELSE IF (guide_plevs.EQ.2) THEN
-	 CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
-	 CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
-	 CALL SendRequest(Req)
-	 CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
+         CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
+         CALL SendRequest(Req)
+         CALL WaitRequest(Req)
             DO l=1,nlevnc
                 DO j=jjb_v,jje_v
@@ -1807,5 +1813,5 @@
     
     ! Variables entree
-    CHARACTER, INTENT(IN)                          :: varname
+    CHARACTER*(*), INTENT(IN)                          :: varname
     INTEGER,   INTENT (IN)                         :: hsize,vsize
     REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
@@ -1817,7 +1823,10 @@
     INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
     INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
+    INTEGER       :: vid_au,vid_av
+    INTEGER       :: l
     INTEGER, DIMENSION (3) :: dim3
     INTEGER, DIMENSION (4) :: dim4,count,start
     INTEGER                :: ierr, varid
+    REAL, DIMENSION (iip1,hsize,vsize) :: field2
     
     CALL gather_field(field,iip1*hsize,vsize,0)
@@ -1834,4 +1843,5 @@
 ! Definition des dimensions
         ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu) 
+        print*,'id_lonu 1 ',id_lonu
         ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv) 
         ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu) 
@@ -1842,4 +1852,5 @@
 ! Creation des variables dimensions
         ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
+        print*,'id_lonu 2 ',id_lonu
         ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
         ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
@@ -1848,4 +1859,6 @@
         ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
         ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
+        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
+        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
         
         ierr=NF_ENDDEF(nid)
@@ -1853,4 +1866,5 @@
 ! Enregistrement des variables dimensions
 #ifdef NC_DOUBLE
+        print*,'id_lonu DOUBLE ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
@@ -1860,5 +1874,8 @@
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
 #else
+        print*,'id_lonu 3 ',id_lonu,rlonu*180./pi
         ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
         ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
@@ -1868,4 +1885,6 @@
         ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
         ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
+        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
+        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
 #endif
 ! --------------------------------------------------------------------
@@ -1875,5 +1894,5 @@
 ! Pressure (GCM)
         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = NF_DEF_VAR(nid,"P",NF_FLOAT,4,dim4,varid)
+        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
@@ -1883,5 +1902,8 @@
 ! Zonal wind
         IF (guide_u) THEN
+        print*,'id_lonu 4 ',id_lonu,varname
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1889,4 +1911,6 @@
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
+            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
+            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
             ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
         ENDIF
@@ -1912,61 +1936,41 @@
     ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
 
+    IF (varname=="SP") timestep=timestep+1
+
+    IF (varname=="SP") THEN
+      print*,'varname=SP=',varname
+    ELSE
+      print*,'varname diff SP=',varname
+    ENDIF
+
+
+    ierr = NF_INQ_VARID(nid,varname,varid)
     SELECT CASE (varname)
-    CASE ("P")
-        timestep=timestep+1
-        ierr = NF_INQ_VARID(nid,"P",varid)
+    CASE ("SP","ps")
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
-#endif
-    CASE ("SP")
-        ierr = NF_INQ_VARID(nid,"ps",varid)
-        start=(/1,1,timestep,0/)
-        count=(/iip1,jjp1,1,0/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("U")
-        ierr = NF_INQ_VARID(nid,"ucov",varid)
+    CASE ("v","va","vcov")
+        start=(/1,1,1,timestep/)
+        count=(/iip1,jjm,llm,1/)
+    CASE DEFAULT
         start=(/1,1,1,timestep/)
         count=(/iip1,jjp1,llm,1/)
+    END SELECT
+
+    SELECT CASE (varname)
+    CASE("u","ua")
+        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
+        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
+    CASE("v","va")
+        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
+    CASE DEFAULT
+        field2=field
+    END SELECT
+
 #ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
 #else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
+    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
 #endif
-    CASE ("V")
-        ierr = NF_INQ_VARID(nid,"vcov",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjm,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("T")
-        ierr = NF_INQ_VARID(nid,"teta",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    CASE ("Q")
-        ierr = NF_INQ_VARID(nid,"q",varid)
-        start=(/1,1,1,timestep/)
-        count=(/iip1,jjp1,llm,1/)
-#ifdef NC_DOUBLE
-        ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
-#else
-        ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
-#endif
-    END SELECT
  
     ierr = NF_CLOSE(nid)
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1302)
@@ -8,4 +8,8 @@
      &                    time_0)
 
+      use exner_hyb_m, only: exner_hyb
+      use exner_milieu_m, only: exner_milieu
+      use exner_hyb_p_m, only: exner_hyb_p
+      use exner_milieu_p_m, only: exner_milieu_p
        USE misc_mod
        USE parallel_lmdz
@@ -17,5 +21,5 @@
        USE vampir
        USE timer_filtre, ONLY : print_filtre_timer
-       USE infotrac
+       USE infotrac, ONLY: nqtot
        USE guide_p_mod, ONLY : guide_main
        USE getparam
@@ -165,5 +169,4 @@
       character*10 string10
 
-      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
       REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
 
@@ -236,5 +239,10 @@
       lafin=.false.
       
-      itaufin   = nday*day_step
+      if (nday>=0) then
+         itaufin   = nday*day_step
+      else
+         ! to run a given (-nday) number of dynamical steps
+         itaufin   = -nday
+      endif
       if (less1day) then
 c MODIF VENUS: to run less than one day:
@@ -292,7 +300,7 @@
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
       if (pressure_exner) then
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
       else
-        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
+        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
       endif
 c$OMP END MASTER
@@ -830,8 +838,10 @@
 c$OMP BARRIER
          if (pressure_exner) then
-           CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
+           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
          else
-           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
          endif
+! Compute geopotential (physics might need it)
+         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
 c$OMP BARRIER
            jD_cur = jD_ref + day_ini - day_ref
@@ -1040,7 +1050,7 @@
 c$OMP BARRIER
           if (pressure_exner) then
-            CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
           else
-            CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
+            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
           endif
 c$OMP BARRIER
@@ -1202,7 +1212,7 @@
 c$OMP BARRIER
         if (pressure_exner) then
-          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
         else
-          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
+          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
         endif
 c$OMP BARRIER
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h	(revision 1302)
@@ -11,5 +11,5 @@
      &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
+     &  ,ok_limit,ok_etat0,hybrid                                       &
      &  ,moyzon_mu,moyzon_ch
 
@@ -19,5 +19,5 @@
      & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
      &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
-     &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
+     &  ,ok_limit,ok_etat0
       logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
                      ! (only used if disvert_type==2)
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90	(revision 1302)
@@ -1,4 +1,4 @@
 ! 
-! $Id: mod_const_mpi.F90 1700 2012-12-20 14:43:19Z lguez $
+! $Id: mod_const_mpi.F90 2055 2014-06-04 12:33:27Z acaubel $
 !
 MODULE mod_const_mpi
@@ -17,5 +17,11 @@
     USE ioipsl_getincom, only: getin
 #endif
-
+! Use of Oasis-MCT coupler 
+#ifdef CPP_OMCT
+    USE mod_prism
+#endif
+#ifdef CPP_XIOS
+    USE wxios, only: wxios_init
+#endif
     IMPLICIT NONE
 #ifdef CPP_MPI
@@ -38,6 +44,10 @@
 #ifdef CPP_COUPLE
 !$OMP MASTER
-       CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
+#ifdef CPP_XIOS
+      CALL wxios_init("LMDZ", outcom=COMM_LMDZ, type_ocean=type_ocean)
+#else
+       CALL prism_init_comp_proto (comp_id, 'LMDZ', ierr)
        CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
+#endif
 !$OMP END MASTER
 #endif
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90	(revision 1300)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90	(revision 1302)
@@ -225,7 +225,10 @@
 #endif
 #ifdef CPP_COUPLE
+! Use of Oasis-MCT coupler 
+#if defined CPP_OMCT
+    use mod_prism
+#else
     use mod_prism_proto
 #endif
-#ifdef CPP_EARTH
 ! Ehouarn: surface_data module is in 'phylmd' ...
       use surface_data, only : type_ocean
@@ -252,4 +255,8 @@
 
       if (type_ocean == 'couple') then
+#ifdef CPP_XIOS
+    !Fermeture propre de XIOS
+      CALL wxios_close()
+#else
 #ifdef CPP_COUPLE
          call prism_terminate_proto(ierr)
@@ -258,4 +265,5 @@
          endif
 #endif 
+#endif
       else
 #ifdef CPP_XIOS
