Index: trunk/LMDZ.COMMON/libf/dyn3d/addfi.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/addfi.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/addfi.F	(revision 1189)
@@ -55,13 +55,19 @@
 c    -----------
 c
-      REAL pdt
+      REAL,INTENT(IN) :: pdt ! time step for the integration (s)
 c
-      REAL pvcov(ip1jm,llm),pucov(ip1jmp1,llm)
-      REAL pteta(ip1jmp1,llm),pq(ip1jmp1,llm,nqtot),pps(ip1jmp1)
+      REAL,INTENT(INOUT) :: pvcov(ip1jm,llm) ! covariant meridional wind
+      REAL,INTENT(INOUT) :: pucov(ip1jmp1,llm) ! covariant zonal wind
+      REAL,INTENT(INOUT) :: pteta(ip1jmp1,llm) ! potential temperature
+      REAL,INTENT(INOUT) :: pq(ip1jmp1,llm,nqtot) ! tracers
+      REAL,INTENT(INOUT) :: pps(ip1jmp1) ! surface pressure (Pa)
+c respective tendencies (.../s) to add
+      REAL,INTENT(IN) :: pdvfi(ip1jm,llm)
+      REAL,INTENT(IN) :: pdufi(ip1jmp1,llm)
+      REAL,INTENT(IN) :: pdqfi(ip1jmp1,llm,nqtot)
+      REAL,INTENT(IN) :: pdhfi(ip1jmp1,llm)
+      REAL,INTENT(IN) :: pdpfi(ip1jmp1)
 c
-      REAL pdvfi(ip1jm,llm),pdufi(ip1jmp1,llm)
-      REAL pdqfi(ip1jmp1,llm,nqtot),pdhfi(ip1jmp1,llm),pdpfi(ip1jmp1)
-c
-      LOGICAL leapf,forward
+      LOGICAL,INTENT(IN) :: leapf,forward ! not used
 c
 c
@@ -71,7 +77,6 @@
       REAL xpn(iim),xps(iim),tpn,tps
       INTEGER j,k,iq,ij
-      REAL qtestw, qtestt
-      PARAMETER ( qtestw = 1.0e-15 )
-      PARAMETER ( qtestt = 1.0e-40 )
+      REAL,PARAMETER :: qtestw = 1.0e-15
+      REAL,PARAMETER :: qtestt = 1.0e-40
 
       REAL SSUM
Index: trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90	(revision 1189)
@@ -9,6 +9,6 @@
   !            M.A Filiberti (04/2002)
   !
-  USE infotrac
-  USE control_mod
+  USE infotrac, ONLY: nqtot, iadv
+  USE control_mod, ONLY: iapp_tracvl, day_step
 
 
@@ -30,15 +30,17 @@
   !     Arguments
   !-------------------------------------------------------------------
+  INTEGER,INTENT(OUT) :: iapptrac
+  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
+  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
+  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
+  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
+  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
+  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
+  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
+  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
+  !-------------------------------------------------------------------
   !     Ajout PPM
   !--------------------------------------------------------
   REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
-  !--------------------------------------------------------
-  INTEGER iapptrac
-  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-  REAL q(ip1jmp1,llm,nqtot),masse(ip1jmp1,llm)
-  REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
-  REAL pk(ip1jmp1,llm)
-  REAL flxw(ip1jmp1,llm)
-
   !-------------------------------------------------------------
   !     Variables locales
Index: trunk/LMDZ.COMMON/libf/dyn3d/caldyn.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/caldyn.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/caldyn.F	(revision 1189)
@@ -1,7 +1,5 @@
 !
-! $Header$
+! $Id: $
 !
-c
-c
       SUBROUTINE caldyn
      $ (itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
@@ -10,19 +8,19 @@
       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=======================================================================
+!=======================================================================
+!
+!  Auteur :  P. Le Van
+!
+!   Objet:
+!   ------
+!
+!   Calcul des tendances dynamiques.
+!
+! Modif 04/93 F.Forget
+!=======================================================================
 
-c-----------------------------------------------------------------------
-c   0. Declarations:
-c   ----------------
+!-----------------------------------------------------------------------
+!   0. Declarations:
+!   ----------------
 
 #include "dimensions.h"
@@ -32,28 +30,36 @@
 #include "comgeom.h"
 
-c   Arguments:
-c   ----------
+!   Arguments:
+!   ----------
 
-      LOGICAL conser
+      LOGICAL,INTENT(IN) :: conser ! triggers printing some diagnostics
+      INTEGER,INTENT(IN) :: itau ! time step index
+      REAL,INTENT(IN) :: vcov(ip1jm,llm) ! covariant meridional wind
+      REAL,INTENT(IN) :: ucov(ip1jmp1,llm) ! covariant zonal wind
+      REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potential temperature
+      REAL,INTENT(IN) :: ps(ip1jmp1) ! surface pressure
+      REAL,INTENT(IN) :: phis(ip1jmp1) ! geopotential at the surface
+      REAL,INTENT(IN) :: pk(ip1jmp1,llm) ! Exner at mid-layer
+      REAL,INTENT(IN) :: pkf(ip1jmp1,llm) ! filtered Exner
+      REAL,INTENT(IN) :: tsurpk(ip1jmp1,llm) ! cpp * temperature / pk
+      REAL,INTENT(IN) :: phi(ip1jmp1,llm) ! geopotential
+      REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass
+      REAL,INTENT(OUT) :: dv(ip1jm,llm) ! tendency on vcov
+      REAL,INTENT(OUT) :: du(ip1jmp1,llm) ! tendency on ucov
+      REAL,INTENT(OUT) :: dteta(ip1jmp1,llm) ! tenddency on teta
+      REAL,INTENT(OUT) :: dp(ip1jmp1) ! tendency on ps
+      REAL,INTENT(OUT) :: w(ip1jmp1,llm) ! vertical velocity
+      REAL,INTENT(OUT) :: pbaru(ip1jmp1,llm) ! mass flux in the zonal direction
+      REAL,INTENT(OUT) :: pbarv(ip1jm,llm) ! mass flux in the meridional direction
+      REAL,INTENT(IN) :: time ! current time
 
-      INTEGER itau
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
-      REAL ps(ip1jmp1),phis(ip1jmp1)
-      REAL pk(ip1jmp1,llm),pkf(ip1jmp1,llm)
-      REAL tsurpk(ip1jmp1,llm)
+!   Local:
+!   ------
+
       REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
-      REAL phi(ip1jmp1,llm),masse(ip1jmp1,llm)
-      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
-      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 ecin(ip1jmp1,llm),convm(ip1jmp1,llm)
       REAL bern(ip1jmp1,llm)
       REAL massebxy(ip1jm,llm)
@@ -62,27 +68,43 @@
       INTEGER   ij,l
 
-c-----------------------------------------------------------------------
-c   Calcul des tendances dynamiques:
-c   --------------------------------
+!-----------------------------------------------------------------------
+!   Compute dynamical tendencies:
+!--------------------------------
 
+      ! compute contravariant winds ucont() and vcont
       CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
+      ! compute pressure p()
       CALL pression ( ip1jmp1, ap      , bp   ,  ps  , p            )
+      ! compute psexbarxy() XY-area weighted-averaged surface pressure (what for?)
       CALL psextbar (   ps   , psexbarxy                            )
+      ! compute mass in each atmospheric mesh: masse()
       CALL massdair (    p   , masse                                )
+      ! compute X and Y-averages of mass, massebx() and masseby()
       CALL massbar  (   masse, massebx , masseby                    )
+      ! compute XY-average of mass, massebxy()
       call massbarxy(   masse, massebxy                             )
+      ! compute mass fluxes pbaru() and pbarv()
       CALL flumass  ( massebx, masseby , vcont, ucont ,pbaru, pbarv )
+      ! compute dteta() , horizontal converging flux of theta
       CALL dteta1   (   teta , pbaru   , pbarv, dteta               )
+      ! compute convm(), horizontal converging flux of mass
       CALL convmas  (   pbaru, pbarv   , convm                      )
 
+      ! compute pressure variation due to mass convergence
       DO ij =1, ip1jmp1
          dp( ij ) = convm( ij,1 ) / airesurg( ij )
       ENDDO
 
+      ! compute vertical velocity w()
       CALL vitvert ( convm  , w                                  )
+      ! compute potential vorticity vorpot()
       CALL tourpot ( vcov   , ucov  , massebxy  , vorpot         )
+      ! compute rotation induced du() and dv()
       CALL dudv1   ( vorpot , pbaru , pbarv     , du     , dv    )
+      ! compute kinetic energy ecin()
       CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
+      ! compute Bernouilli function bern()
       CALL bernoui ( ip1jmp1, llm   , phi       , ecin   , bern  )
+      ! compute and add du() and dv() contributions from Bernouilli and pressure
       CALL dudv2   ( tsurpk , pkf   , bern      , du     , dv    )
 
@@ -94,23 +116,24 @@
       ENDDO
 
-
+      ! compute vertical advection contributions to du(), dv() and dteta()
       CALL advect( ang, vcov, teta, w, massebx, masseby, du, dv,dteta ) 
 
-C  WARNING probleme de peridocite de dv sur les PC/linux. Pb d'arrondi 
-C          probablement. Observe sur le code compile avec pgf90 3.0-1 
+!  WARNING probleme de peridocite de dv sur les PC/linux. Pb d'arrondi 
+!          probablement. Observe sur le code compile avec pgf90 3.0-1 
 
       DO l = 1, llm
          DO ij = 1, ip1jm, iip1
            IF( dv(ij,l).NE.dv(ij+iim,l) )  THEN
-c         PRINT *,'!!!ATTENTION!!! probleme de periodicite sur vcov',  
-c    ,   ' dans caldyn'
-c         PRINT *,' l,  ij = ', l, ij, ij+iim,dv(ij+iim,l),dv(ij,l)
+!         PRINT *,'!!!ATTENTION!!! probleme de periodicite sur vcov',  
+!    ,   ' dans caldyn'
+!         PRINT *,' l,  ij = ', l, ij, ij+iim,dv(ij+iim,l),dv(ij,l)
           dv(ij+iim,l) = dv(ij,l)
           endif
          enddo
       enddo
-c-----------------------------------------------------------------------
-c   Sorties eventuelles des variables de controle:
-c   ----------------------------------------------
+
+!-----------------------------------------------------------------------
+!   Output some control variables:
+!---------------------------------
 
       IF( conser )  THEN
@@ -120,4 +143,3 @@
       ENDIF
 
-      RETURN
       END
Index: trunk/LMDZ.COMMON/libf/dyn3d/calfis.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/calfis.F	(revision 1189)
@@ -29,6 +29,6 @@
 c    Auteur :  P. Le Van, F. Hourdin 
 c   .........
-      USE infotrac
-      USE control_mod
+      USE infotrac, ONLY: nqtot, niadv, tname
+      USE control_mod, ONLY: planet_type, nsplit_phys
       USE write_field
       USE cpdet_mod, only: t2tpot,tpot2t
@@ -106,33 +106,34 @@
 c    Arguments :
 c    -----------
-      LOGICAL  lafin
-
-
-      REAL pvcov(iip1,jjm,llm)
-      REAL pucov(iip1,jjp1,llm)
-      REAL pteta(iip1,jjp1,llm)
-      REAL pmasse(iip1,jjp1,llm)
-      REAL pq(iip1,jjp1,llm,nqtot)
-      REAL pphis(iip1,jjp1)
-      REAL pphi(iip1,jjp1,llm)
-
-      REAL pdvcov(iip1,jjm,llm)
-      REAL pducov(iip1,jjp1,llm)
-      REAL pdteta(iip1,jjp1,llm)
+      LOGICAL,INTENT(IN) ::  lafin ! .true. for the very last call to physics
+      REAL,INTENT(IN) :: jD_cur, jH_cur
+      REAL,INTENT(IN) :: pvcov(iip1,jjm,llm) ! covariant meridional velocity
+      REAL,INTENT(IN) :: pucov(iip1,jjp1,llm) ! covariant zonal velocity
+      REAL,INTENT(IN) :: pteta(iip1,jjp1,llm) ! potential temperature
+      REAL,INTENT(IN) :: pmasse(iip1,jjp1,llm) ! mass in each cell ! not used
+      REAL,INTENT(IN) :: pq(iip1,jjp1,llm,nqtot) ! tracers
+      REAL,INTENT(IN) :: pphis(iip1,jjp1) ! surface geopotential
+      REAL,INTENT(IN) :: pphi(iip1,jjp1,llm) ! geopotential
+
+      REAL,INTENT(IN) :: pdvcov(iip1,jjm,llm) ! dynamical tendency on vcov
+      REAL,INTENT(IN) :: pducov(iip1,jjp1,llm) ! dynamical tendency on ucov
+      REAL,INTENT(IN) :: pdteta(iip1,jjp1,llm) ! dynamical tendency on teta
 ! commentaire SL: pdq ne sert que pour le calcul de pcvgq,
 ! qui lui meme ne sert a rien dans la routine telle qu'elle est
 ! ecrite, et que j'ai donc commente....
-      REAL pdq(iip1,jjp1,llm,nqtot)
-
-      REAL pps(iip1,jjp1)
-      REAL pp(iip1,jjp1,llmp1)
-      REAL ppk(iip1,jjp1,llm)
-
-c TENDENCIES in */s 
-      REAL pdvfi(iip1,jjm,llm)
-      REAL pdufi(iip1,jjp1,llm)
-      REAL pdhfi(iip1,jjp1,llm)
-      REAL pdqfi(iip1,jjp1,llm,nqtot)
-      REAL pdpsfi(iip1,jjp1)
+      REAL,INTENT(IN) :: pdq(iip1,jjp1,llm,nqtot) ! dynamical tendency on tracers
+      ! NB: pdq is only used to compute pcvgq which is in fact not used...
+
+      REAL,INTENT(IN) :: pps(iip1,jjp1) ! surface pressure (Pa)
+      REAL,INTENT(IN) :: pp(iip1,jjp1,llmp1) ! pressure at mesh interfaces (Pa)
+      REAL,INTENT(IN) :: ppk(iip1,jjp1,llm) ! Exner at mid-layer
+      REAL,INTENT(IN) :: flxw(iip1,jjp1,llm)  ! Vertical mass flux on dynamics grid
+
+      ! tendencies (in */s) from the physics
+      REAL,INTENT(OUT) :: pdvfi(iip1,jjm,llm) ! tendency on covariant meridional wind
+      REAL,INTENT(OUT) :: pdufi(iip1,jjp1,llm) ! tendency on covariant zonal wind
+      REAL,INTENT(OUT) :: pdhfi(iip1,jjp1,llm) ! tendency on potential temperature (K/s)
+      REAL,INTENT(OUT) :: pdqfi(iip1,jjp1,llm,nqtot) ! tendency on tracers
+      REAL,INTENT(OUT) :: pdpsfi(iip1,jjp1) ! tendency on surface pressure (Pa/s)
 
 
@@ -173,20 +174,14 @@
 
 cIM diagnostique PVteta, Amip2
-      INTEGER ntetaSTD
-      PARAMETER(ntetaSTD=3)
-      REAL rtetaSTD(ntetaSTD)
-      DATA rtetaSTD/350., 380., 405./ ! Earth-specific values, beware !!
+      INTEGER,PARAMETER :: ntetaSTD=3
+      REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
       REAL PVteta(ngridmx,ntetaSTD)
 
-      REAL flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
       REAL flxwfi(ngridmx,llm)  ! Flux de masse verticale sur la grille physiq
       
       REAL SSUM
 
-      LOGICAL firstcal, debut
-      DATA firstcal/.true./
-      SAVE firstcal,debut
+      LOGICAL,SAVE :: firstcal=.true., debut=.true.
 !      REAL rdayvrai
-      REAL, intent(in):: jD_cur, jH_cur
 
       LOGICAL tracerdyn ! for generic/mars physics call ; possibly to get rid of
Index: trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F	(revision 1189)
@@ -174,5 +174,10 @@
       starttime = 0
       CALL getin('starttime',starttime)
-
+      
+      ! Mars: time of start for run in "start.nc" (when there are multiple time
+      !       steps stored in the file)
+      timestart=-9999 ! default value; if <0, use last stored time
+      call getin("timestart",timestart)
+      
 !Config  Key  = less1day
 !Config  Desc = Possibilite d'integrer moins d'un jour
Index: trunk/LMDZ.COMMON/libf/dyn3d/dissip.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/dissip.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/dissip.F	(revision 1189)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id: $
 !
       SUBROUTINE dissip( vcov,ucov,teta,p, dv,du,dh )
@@ -35,7 +35,12 @@
 c   ----------
 
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
-      REAL  p( ip1jmp1,llmp1 )
-      REAL dv(ip1jm,llm),du(ip1jmp1,llm),dh(ip1jmp1,llm)
+      REAL,INTENT(IN) :: vcov(ip1jm,llm) ! covariant meridional wind
+      REAL,INTENT(IN) :: ucov(ip1jmp1,llm) ! covariant zonal wind
+      REAL,INTENT(IN) :: teta(ip1jmp1,llm) ! potential temperature
+      REAL,INTENT(IN) :: p(ip1jmp1,llmp1) ! pressure
+      ! tendencies (.../s) on covariant winds and potential temperature
+      REAL,INTENT(OUT) :: dv(ip1jm,llm)
+      REAL,INTENT(OUT) :: du(ip1jmp1,llm)
+      REAL,INTENT(OUT) :: dh(ip1jmp1,llm)
 
 c   Local:
Index: trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F	(revision 1189)
@@ -12,8 +12,8 @@
       use IOIPSL
 #endif
-      USE infotrac
+      USE infotrac, ONLY: nqtot
       USE guide_mod, ONLY : guide_main
-      USE write_field
-      USE control_mod, only: planet_type,nday,day_step,iperiod,iphysiq,
+      USE write_field, ONLY: writefield
+      USE control_mod, ONLY: planet_type,nday,day_step,iperiod,iphysiq,
      &                       less1day,fractday,ndynstep,iconser,
      &                       dissip_period,offline,ip_ebil_dyn,
@@ -75,22 +75,26 @@
 ! #include "clesphys.h"
 
-      real zqmin,zqmax
-
-c   variables dynamiques
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
-      REAL ps(ip1jmp1)                       ! pression  au sol
-      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL masse(ip1jmp1,llm)                ! masse d'air
-      REAL phis(ip1jmp1)                     ! geopotentiel au sol
-      REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL w(ip1jmp1,llm)                    ! vitesse verticale
+      REAL,INTENT(IN) :: time_0 ! not used
+
+c   dynamical variables:
+      REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm)    ! zonal covariant wind
+      REAL,INTENT(INOUT) :: vcov(ip1jm,llm)      ! meridional covariant wind
+      REAL,INTENT(INOUT) :: teta(ip1jmp1,llm)    ! potential temperature
+      REAL,INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
+      REAL,INTENT(INOUT) :: masse(ip1jmp1,llm)   ! air mass
+      REAL,INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
+      REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
+
+      REAL p (ip1jmp1,llmp1  )               ! interlayer pressure
+      REAL pks(ip1jmp1)                      ! exner at the surface
+      REAL pk(ip1jmp1,llm)                   ! exner at mid-layer
+      REAL pkf(ip1jmp1,llm)                  ! filtered exner at mid-layer
+      REAL phi(ip1jmp1,llm)                  ! geopotential
+      REAL w(ip1jmp1,llm)                    ! vertical velocity
 ! ADAPTATION GCM POUR CP(T)
       REAL temp(ip1jmp1,llm)                 ! temperature  
       REAL tsurpk(ip1jmp1,llm)               ! cpp*T/pk  
+
+      real zqmin,zqmax
 
 c variables dynamiques intermediaire pour le transport
@@ -133,5 +137,4 @@
 
       REAL  SSUM
-      REAL time_0
 !     REAL finvmaold(ip1jmp1,llm)
 
@@ -410,9 +413,6 @@
 c   -------------------------------------------------------------
 
-!      IF( forward. OR . leapf )  THEN
-! Ehouarn: NB: at this point p with ps are not synchronized
-!              (whereas mass and ps are...)
-      IF((.not.forward).OR. leapf )  THEN
-        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
+      IF( forward. OR . leapf )  THEN
+! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
          CALL caladvtrac(q,pbaru,pbarv,
      *        p, masse, dq,  teta,
@@ -547,5 +547,13 @@
      $                  ucov, vcov, teta , q   ,ps ,
      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
-
+          ! since addfi updates ps(), also update p(), masse() and pk()
+          CALL pression (ip1jmp1,ap,bp,ps,p)
+          CALL massdair(p,masse)
+          if (pressure_exner) then
+            CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
+          else
+            CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
+          endif
+          
 c      Couche superieure :
 c      -------------------
@@ -600,5 +608,5 @@
           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
-
+        CALL massdair(p,masse)
 
 c-----------------------------------------------------------------------
