Index: trunk/LMDZ.COMMON/libf/dyn3dpar/addfi_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/addfi_p.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/addfi_p.F	(revision 1189)
@@ -55,13 +55,19 @@
 c    -----------
 c
-      REAL pdt
-c
-      REAL pvcov(ip1jm,llm),pucov(ip1jmp1,llm)
-      REAL pteta(ip1jmp1,llm),pq(ip1jmp1,llm,nqtot),pps(ip1jmp1)
-c
-      REAL pdvfi(ip1jm,llm),pdufi(ip1jmp1,llm)
-      REAL pdqfi(ip1jmp1,llm,nqtot),pdhfi(ip1jmp1,llm),pdpfi(ip1jmp1)
-c
-      LOGICAL leapf,forward
+      REAL,INTENT(IN) :: pdt ! time step for the integration (s)
+c
+      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
+      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/dyn3dpar/advtrac_p.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 1189)
@@ -16,6 +16,6 @@
   USE Vampir
   USE times
-  USE infotrac
-  USE control_mod
+  USE infotrac, ONLY: nqtot, iadv
+  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
   IMPLICIT NONE
   !
@@ -34,15 +34,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/dyn3dpar/caldyn_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/caldyn_p.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/caldyn_p.F	(revision 1189)
@@ -1,7 +1,5 @@
 !
-! $Header$
-!
-c
-c
+! $Id: $
+!
 #undef DEBUG_IO
 c#define DEBUG_IO
@@ -15,19 +13,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=======================================================================
-
-c-----------------------------------------------------------------------
-c   0. Declarations:
-c   ----------------
+!=======================================================================
+!
+!  Auteur :  P. Le Van
+!
+!   Objet:
+!   ------
+!
+!   Calcul des tendances dynamiques.
+!
+! Modif 04/93 F.Forget
+!=======================================================================
+
+!-----------------------------------------------------------------------
+!   0. Declarations:
+!   ----------------
 
 #include "dimensions.h"
@@ -37,25 +35,32 @@
 #include "comgeom.h"
 
-c   Arguments:
-c   ----------
-
-      LOGICAL conser
-
-      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)
+!   Arguments:
+!   ----------
+
+      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
+
+!   Local:
+!   ------
+
       REAL,SAVE :: 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 w(ip1jmp1,llm)
-      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
-      REAL time
-
-c   Local:
-c   ------
-
       REAL,SAVE :: ang(ip1jmp1,llm)
       REAL,SAVE :: p(ip1jmp1,llmp1)
@@ -70,16 +75,25 @@
       INTEGER   ij,l,ijb,ije,ierr
 
-c-----------------------------------------------------------------------
-c   Calcul des tendances dynamiques:
-c   --------------------------------
+!-----------------------------------------------------------------------
+!   Compute dynamical tendencies:
+!--------------------------------
+
+      ! compute contravariant winds ucont() and vcont
       CALL covcont_p  ( llm    , ucov    , vcov , ucont, vcont        )
+      ! compute pressure p()
       CALL pression_p ( ip1jmp1, ap      , bp   ,  ps  , p            )
 cym      CALL psextbar (   ps   , psexbarxy                          )
 c$OMP BARRIER
+      ! compute mass in each atmospheric mesh: masse()
       CALL massdair_p (    p   , masse                                )
+      ! compute X and Y-averages of mass, massebx() and masseby()
       CALL massbar_p  (   masse, massebx , masseby                    )
+      ! compute XY-average of mass, massebxy()
       call massbarxy_p(   masse, massebxy                             )
+      ! compute mass fluxes pbaru() and pbarv()
       CALL flumass_p  ( massebx, masseby , vcont, ucont ,pbaru, pbarv )
+      ! compute dteta() , horizontal converging flux of theta
       CALL dteta1_p   (   teta , pbaru   , pbarv, dteta               )
+      ! compute convm(), horizontal converging flux of mass
       CALL convmas1_p  (   pbaru, pbarv   , convm                      )
 c$OMP BARRIER      
@@ -108,5 +122,5 @@
       ijb=ij_begin
       ije=ij_end
-            
+      ! compute pressure variation due to mass convergence
       DO ij =ijb, ije
          dp( ij ) = convm( ij,1 ) / airesurg( ij )
@@ -115,6 +129,10 @@
 c$OMP BARRIER
 c$OMP FLUSH
+
+      ! compute vertical velocity w()
       CALL vitvert_p ( convm  , w                                  )
+      ! compute potential vorticity vorpot()
       CALL tourpot_p ( vcov   , ucov  , massebxy  , vorpot         )
+      ! compute rotation induced du() and dv()
       CALL dudv1_p   ( vorpot , pbaru , pbarv     , du     , dv    )
 
@@ -129,6 +147,10 @@
 c$OMP BARRIER
 #endif      
+      
+      ! compute kinetic energy ecin()
       CALL enercin_p ( vcov   , ucov  , vcont     , ucont  , ecin  )
+      ! compute Bernouilli function bern()
       CALL bernoui_p ( ip1jmp1, llm   , phi       , ecin   , bern  )
+      ! compute and add du() and dv() contributions from Bernouilli and pressure
       CALL dudv2_p   ( tsurpk , pkf   , bern      , du     , dv    )
 
@@ -159,4 +181,5 @@
 c$OMP END DO
 
+      ! compute vertical advection contributions to du(), dv() and dteta()
       CALL advect_new_p(ang,vcov,teta,w,massebx,masseby,du,dv,dteta) 
 
@@ -190,4 +213,3 @@
       ENDIF
 
-      RETURN
       END
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F	(revision 1189)
@@ -37,6 +37,6 @@
       Use Write_field_p
       USE Times
-      USE infotrac
-      USE control_mod
+      USE infotrac, ONLY: nqtot, niadv, tname
+      USE control_mod, ONLY: planet_type, nsplit_phys
       USE cpdet_mod, only: tpot2t_p, t2tpot_p
 
@@ -116,33 +116,34 @@
 c    Arguments :
 c    -----------
-      LOGICAL  lafin
-!      REAL heure
-      REAL, intent(in):: jD_cur, jH_cur
-      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 flxw(iip1,jjp1,llm)  ! Flux de masse verticale sur la grille dynamique
-
-      REAL pps(iip1,jjp1)
-      REAL pp(iip1,jjp1,llmp1)
-      REAL ppk(iip1,jjp1,llm)
-
-      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)
 
 #ifdef CPP_PHYS
@@ -225,8 +226,6 @@
 
 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(klon,ntetaSTD)
       
@@ -234,7 +233,5 @@
       REAL SSUM
 
-      LOGICAL firstcal, debut
-      DATA firstcal/.true./
-      SAVE firstcal,debut
+      LOGICAL,SAVE :: firstcal=.true., debut=.true.
 c$OMP THREADPRIVATE(firstcal,debut)
       
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F	(revision 1189)
@@ -202,4 +202,9 @@
       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/dyn3dpar/dissip_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/dissip_p.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/dissip_p.F	(revision 1189)
@@ -1,2 +1,5 @@
+!
+! $Id: $
+!
       SUBROUTINE dissip_p( vcov,ucov,teta,p, dv,du,dh )
 c
@@ -34,7 +37,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) ! potentail 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/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1129)
+++ trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F	(revision 1189)
@@ -83,22 +83,26 @@
 #include "academic.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,SAVE :: p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL,SAVE :: pks(ip1jmp1)                      ! exner au  sol
-      REAL,SAVE :: pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL,SAVE :: pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL :: masse(ip1jmp1,llm)                ! masse d'air
-      REAL :: phis(ip1jmp1)                     ! geopotentiel au sol
-      REAL,SAVE :: phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL,SAVE :: 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,SAVE :: p (ip1jmp1,llmp1  )       ! interlayer pressure
+      REAL,SAVE :: pks(ip1jmp1)              ! exner at the surface
+      REAL,SAVE :: pk(ip1jmp1,llm)           ! exner at mid-layer
+      REAL,SAVE :: pkf(ip1jmp1,llm)          ! filtered exner at mid-layer
+      REAL,SAVE :: phi(ip1jmp1,llm)          ! geopotential
+      REAL,SAVE :: w(ip1jmp1,llm)            ! vertical velocity
 ! ADAPTATION GCM POUR CP(T)
       REAL,SAVE :: temp(ip1jmp1,llm)                 ! temperature  
       REAL,SAVE :: tsurpk(ip1jmp1,llm)               ! cpp*T/pk  
+
+      real zqmin,zqmax
 
 c variables dynamiques intermediaire pour le transport
@@ -144,6 +148,5 @@
       REAL       time 
 
-      REAL  SSUM
-      REAL time_0 
+      REAL  SSUM 
 !      REAL,SAVE :: finvmaold(ip1jmp1,llm)
 
@@ -721,9 +724,6 @@
 c   -------------------------------------------------------------
 
-!      IF( forward. OR . leapf )  THEN
-      IF((.not.forward).OR. leapf )  THEN
-        ! Ehouarn: gather mass fluxes during backward Matsuno or LF step
-cc$OMP PARALLEL DEFAULT(SHARED) 
-c
+      IF( forward. OR . leapf )  THEN
+! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
          CALL advtrac_p( pbaru,pbarv, 
      *             p,  masse,q,iapptrac, teta,
@@ -737,5 +737,4 @@
 
       ENDIF ! of IF( forward. OR . leapf )
-cc$OMP END PARALLEL
 
 c-----------------------------------------------------------------------
@@ -1037,4 +1036,15 @@
      $                  ucov, vcov, teta , q   ,ps ,
      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
+          ! since addfi updates ps(), also update p(), masse() and pk()
+          CALL pression_p(ip1jmp1,ap,bp,ps,p)
+c$OMP BARRIER
+          CALL massdair_p(p,masse)
+c$OMP BARRIER
+          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)
+          endif
+c$OMP BARRIER
 
 c      Couche superieure :
@@ -1195,4 +1205,6 @@
           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
         endif
+c$OMP BARRIER
+        CALL massdair_p(p,masse)
 c$OMP BARRIER
 
