Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1311)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 1312)
@@ -5,5 +5,5 @@
      $            ,pplev,pplay,pphi
      $            ,pu,pv,pt,pq
-     $            ,pw
+     $            ,flxw
      $            ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
 
@@ -128,5 +128,5 @@
 c    ptdyn(ngrid,nlayer)    | corresponding variables
 c    pqdyn(ngrid,nlayer,nq) |
-c    pw(ngrid,?)           vertical velocity
+c    flxw(ngrid,nlayer)      vertical mass flux (kg/s) at layer lower boundary
 c
 c   output:
@@ -159,24 +159,31 @@
 c   inputs:
 c   -------
-      INTEGER ngrid,nlayer,nq
-      REAL ptimestep
-      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL pphi(ngrid,nlayer)
-      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
-      REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
-      REAL pw(ngrid,nlayer) !Mars pvervel transmit par dyn3d
-      REAL zh(ngrid,nlayer)      ! potential temperature (K)
-      LOGICAL firstcall,lastcall
-
-      REAL pday
-      REAL ptime
-      logical tracerdyn
+      INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers
+      INTEGER,INTENT(in) :: nq ! number of tracers
+      LOGICAL,INTENT(in) :: firstcall ! signals first call to physics
+      LOGICAL,INTENT(in) :: lastcall ! signals last call to physics
+      REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0
+      REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon)
+      REAL,INTENT(in) :: ptimestep ! physics timestep (s)
+      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
+      REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2)
+      REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s)
+      REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s)
+      REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
+      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
+      REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s)
+                                            ! at lower boundary of layer
 
 c   outputs:
 c   --------
 c     physical tendencies
-      REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer)
-      REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
-      REAL pdpsrf(ngrid) ! surface pressure tendency
+      REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s)
+      REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s)
+      REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
+      REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
+      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
+      LOGICAL,INTENT(out) :: tracerdyn ! signal to the dynamics to advect tracers or not
 
 
@@ -273,4 +280,5 @@
       REAL ztime_fin
       REAL zdh(ngrid,nlayer)
+      REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards)
       INTEGER length
       PARAMETER (length=100)
@@ -569,4 +577,19 @@
       endif
 #endif
+
+     ! Compute vertical velocity (m/s) from vertical mass flux
+     ! w = F / (rho*area) and rho = r*T/P
+     ! but first linearly interpolate mass flux to mid-layers
+     do l=1,nlayer-1
+       pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
+     enddo
+     pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
+     do l=1,nlayer
+       pw(1:ngrid,l)=(pw(1:ngrid,l)*pplay(1:ngrid,l)) /  &
+                     (r*pt(1:ngrid,l)*area(1:ngrid))
+       ! NB: here we use r and nor rnew since this diagnostic comes
+       ! from the dynamics
+     enddo
+
 c-----------------------------------------------------------------------
 c    2. Compute radiative tendencies :