Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 712)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 719)
@@ -322,9 +322,12 @@
                                !   Qabs instead of Qext
                                !   (direct comparison with TES)
+                               
+      REAL dqdustsurf(ngridmx) ! surface q dust flux (kg/m2/s)
+      REAL dndustsurf(ngridmx) ! surface n dust flux (number/m2/s)
 
 c Test 1d/3d scavenging
       real h2otot(ngridmx)
-      REAL satu(ngridmx,nlayermx) ! satu ratio for output
-      REAL zqsat(ngridmx,nlayermx)    ! saturation
+      REAL satu(ngridmx,nlayermx)  ! satu ratio for output
+      REAL zqsat(ngridmx,nlayermx) ! saturation
 
       REAL time_phys
@@ -1474,6 +1477,23 @@
 
          if (tracer) then
+         
+           if(doubleq) then
+              do ig=1,ngrid 
+                dqdustsurf(ig) = 
+     &                dqsurf(ig,igcm_dust_mass)*tauscaling(ig)
+                dndustsurf(ig) = 
+     &                dqsurf(ig,igcm_dust_number)*tauscaling(ig)
+              enddo
+              if (scavenging) then
+                do ig=1,ngrid 
+                  dqdustsurf(ig) = dqdustsurf(ig) + 
+     &                     dqsurf(ig,igcm_ccn_mass)*tauscaling(ig)
+                  dndustsurf(ig) = dndustsurf(ig) + 
+     &                     dqsurf(ig,igcm_ccn_number)*tauscaling(ig)
+                enddo
+              endif
+           endif
+                    
            if (water) then
-
              mtot(:)=0
              icetot(:)=0
@@ -1588,15 +1608,15 @@
              if (water) then
                vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
-     &                  *mugaz/mmol(igcm_h2o_vap)
-               call wstats(ngrid,"vmr_h2ovapor",
+     &      *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
+               call wstats(ngrid,"vmr_h2ovap",
      &                    "H2O vapor volume mixing ratio","mol/mol",
      &                    3,vmr)
                vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-     &                  *mugaz/mmol(igcm_h2o_ice)
+     &      *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_ice)
                call wstats(ngrid,"vmr_h2oice",
      &                    "H2O ice volume mixing ratio","mol/mol",
      &                    3,vmr)
                vmr=zqsat(1:ngridmx,1:nlayermx)
-     &                  *mugaz/mmol(igcm_h2o_vap)
+     &      *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
                call wstats(ngrid,"vmr_h2osat",
      &                    "saturation volume mixing ratio","mol/mol",
@@ -1605,7 +1625,7 @@
      &                    "surface h2o_ice","kg/m2",
      &                    2,qsurf(1,igcm_h2o_ice))
-               call wstats(ngrid,'albedo',
-     &                         'albedo',
-     &                         '',2,albedo(1:ngridmx,1))
+c               call wstats(ngrid,'albedo',
+c     &                         'albedo',
+c     &                         '',2,albedo(1:ngridmx,1))
                call wstats(ngrid,"mtot",
      &                    "total mass of water vapor","kg/m2",
@@ -1637,4 +1657,52 @@
 
              endif ! of if (water)
+             
+             
+           if (dustbin.ne.0) then
+          
+             call wstats(ngridmx,'tau','taudust','SI',2,tau(1,1))
+             
+             if (doubleq) then
+c            call wstats(ngridmx,'qsurf','qsurf',
+c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
+c            call wstats(ngridmx,'Nsurf','N particles',
+c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
+c            call wstats(ngridmx,'dqsdev','ddevil lift',
+c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
+c            call wstats(ngridmx,'dqssed','sedimentation',
+c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
+c            call wstats(ngridmx,'dqsdif','diffusion',
+c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
+               call wstats(ngridmx,'dqsdust',
+     &                        'deposited surface dust mass',
+     &                        'kg.m-2.s-1',2,dqdustsurf)
+               call wstats(ngridmx,'dqndust',
+     &                        'deposited surface dust number',
+     &                        'number.m-2.s-1',2,dndustsurf)
+               call wstats(ngridmx,'reffdust','reffdust',
+     &                        'm',3,rdust*ref_r0)
+               call wstats(ngridmx,'dustq','Dust mass mr',
+     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
+               call wstats(ngridmx,'dustN','Dust number',
+     &                        'part/kg',3,pq(1,1,igcm_dust_number))
+             else
+               do iq=1,dustbin
+                 write(str2(1:2),'(i2.2)') iq
+                 call wstats(ngridmx,'q'//str2,'mix. ratio',
+     &                         'kg/kg',3,zq(1,1,iq))
+                 call wstats(ngridmx,'qsurf'//str2,'qsurf',
+     &                         'kg.m-2',2,qsurf(1,iq))
+               end do
+             endif ! (doubleq)
+
+             if (scavenging) then
+               call wstats(ngridmx,'ccnq','CCN mass mr',
+     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
+               call wstats(ngridmx,'ccnN','CCN number',
+     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
+             endif ! (scavenging)
+          
+          endif ! (dustbin.ne.0)
+
 
              if (thermochem.or.photochem) then
@@ -1643,5 +1711,7 @@
      $                 noms(iq) .ne. "dust_number" .and.
      $                 noms(iq) .ne. "ccn_mass" .and.
-     $                 noms(iq) .ne. "ccn_number") then
+     $                 noms(iq) .ne. "ccn_number" .and.
+     $                 noms(iq) .ne. "h2o_vap" .and.
+     $                 noms(iq) .ne. "h2o_ice") then
                     vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq)
      &                          *mmean(1:ngrid,1:nlayer)/mmol(iq)
@@ -1710,5 +1780,5 @@
         qsurfice(1:ngrid) = qsurf(1:ngrid,igcm_h2o_ice)
         vmr=1.e6 * zq(1:ngrid,1:nlayer,igcm_h2o_ice)
-     .           * mugaz / mmol(igcm_h2o_ice)
+     .           *mmean(1:ngridmx,1:nlayermx) / mmol(igcm_h2o_ice)
       ENDIF
       !! Dust quantity integration along the vertical axe
@@ -1745,25 +1815,25 @@
      &                  tsurf)
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
-        call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
+         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
      &                  co2ice)
 
-c         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
-c     &                  "K",2,zt(1,7))
-c         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
-c     &                  fluxsurf_lw)
-c         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
-c     &                  fluxsurf_sw_tot)
-c         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
-c     &                  fluxtop_lw)
-c         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
-c     &                  fluxtop_sw_tot)
-c         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
-c        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
-c        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
-c        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
-c         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
+         call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7",
+     &                  "K",2,zt(1,7))
+         call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
+     &                  fluxsurf_lw)
+         call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2,
+     &                  fluxsurf_sw_tot)
+         call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
+     &                  fluxtop_lw)
+         call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2,
+     &                  fluxtop_sw_tot)
+         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
+         call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
+         call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
+         call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
+         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
 c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
 c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
-c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
+         call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
 c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
 c    &                  zstress)
@@ -1808,6 +1878,6 @@
 !          call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer",
 !    &                     "kg/kg",2,zq(1,1,igcm_co2))
-!          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
-!    &                     "kg/kg",3,zq(1,1,igcm_co2))
+          call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio",
+     &                     "kg/kg",3,zq(1,1,igcm_co2))
         
          ! Compute co2 column
@@ -1844,30 +1914,30 @@
 #endif
 
-             CALL WRITEDIAGFI(ngridmx,'mtot',
+            CALL WRITEDIAGFI(ngridmx,'mtot',
      &                       'total mass of water vapor',
      &                       'kg/m2',2,mtot)
-             CALL WRITEDIAGFI(ngridmx,'icetot',
+            CALL WRITEDIAGFI(ngridmx,'icetot',
      &                       'total mass of water ice',
      &                       'kg/m2',2,icetot)
-c            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
-c     &                *mugaz/mmol(igcm_h2o_ice)
-c            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
-c     &                       'mol/mol',3,vmr)
-c            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
-c     &                *mugaz/mmol(igcm_h2o_vap)
-c            call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
-c     &                       'mol/mol',3,vmr)
-             CALL WRITEDIAGFI(ngridmx,'reffice',
+            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)
+     &      *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_ice)
+            call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr',
+     &                       'mol/mol',3,vmr)
+            vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)
+     &      *mmean(1:ngridmx,1:nlayermx)/mmol(igcm_h2o_vap)
+            call WRITEDIAGFI(ngridmx,'vmr_h2ovap','h2o vap vmr',
+     &                       'mol/mol',3,vmr)
+            CALL WRITEDIAGFI(ngridmx,'reffice',
      &                       'Mean reff',
      &                       'm',2,rave)
-             CALL WRITEDIAGFI(ngrid,"Nccntot",
+            CALL WRITEDIAGFI(ngrid,"Nccntot",
      &                    "condensation nuclei","Nbr/m2",
      &                    2,Nccntot)
-c             CALL WRITEDIAGFI(ngrid,"Mccntot",
-c     &                    "mass condensation nuclei","kg/m2",
-c     &                    2,Mccntot)
-c            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
-c     &                       'm',3,rice)
-             call WRITEDIAGFI(ngridmx,'h2o_ice_s',
+            CALL WRITEDIAGFI(ngrid,"Mccntot",
+     &                    "mass condensation nuclei","kg/m2",
+     &                    2,Mccntot)
+            call WRITEDIAGFI(ngridmx,'rice','Ice particle size',
+     &                       'm',3,rice)
+            call WRITEDIAGFI(ngridmx,'h2o_ice_s',
      &                       'surface h2o_ice',
      &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
@@ -1902,5 +1972,5 @@
 
          if (tracer.and.(dustbin.ne.0)) then
-c          call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
+          call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1))
            if (doubleq) then
 c            call WRITEDIAGFI(ngridmx,'qsurf','qsurf',
@@ -1914,10 +1984,16 @@
 c            call WRITEDIAGFI(ngridmx,'dqsdif','diffusion',
 c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
-c             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
-c     &                        'm',3,rdust*ref_r0)
-c             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
-c     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
-c             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
-c     &                        'part/kg',3,pq(1,1,igcm_dust_number))
+             call WRITEDIAGFI(ngridmx,'dqsdust',
+     &                        'deposited surface dust mass',
+     &                        'kg.m-2.s-1',2,dqdustsurf)
+             call WRITEDIAGFI(ngridmx,'dqndust',
+     &                        'deposited surface dust number',
+     &                        'number.m-2.s-1',2,dndustsurf)
+             call WRITEDIAGFI(ngridmx,'reffdust','reffdust',
+     &                        'm',3,rdust*ref_r0)
+             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
+     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
+             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
+     &                        'part/kg',3,pq(1,1,igcm_dust_number))
 #ifdef MESOINI
              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
@@ -1941,8 +2017,8 @@
 
            if (scavenging) then
-c             call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
-c     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
-c             call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
-c     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
+             call WRITEDIAGFI(ngridmx,'ccnq','CCN mass mr',
+     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
+             call WRITEDIAGFI(ngridmx,'ccnN','CCN number',
+     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
            endif ! (scavenging)
 
Index: trunk/LMDZ.MARS/libf/phymars/wstats.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 712)
+++ trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 719)
@@ -18,5 +18,5 @@
 character (len=50) :: namebis
 character (len=50), save :: firstvar
-integer :: ierr,varid,nbdim,nid
+integer :: ierr,ierr2,varid,nbdim,nid
 integer :: meanid,sdid
 integer, dimension(4)  :: id,start,size
@@ -24,13 +24,63 @@
 integer :: l,i,j,ig0
 integer,save :: index
+! Added to use stats.def to select output variable
+logical,save :: stats_def
+logical :: getout
+integer,save :: n_nom_def
+integer :: n
+integer,parameter :: n_nom_def_max=99
+character(len=20),save :: nom_def(n_nom_def_max)
+logical, save :: notfoundfirst=.TRUE.
 
 integer, save :: step=0
-
+      
+      
 
 if (firstcall) then
    firstcall=.false.
-   firstvar=trim((nom))
    call inistats(ierr)
-endif
+! at very first call, check if there is a "stats.def" to use and read it
+   open(99,file="stats.def",status='old',form='formatted',iostat=ierr2)
+   if(ierr2.eq.0) then
+     stats_def=.true.
+     write(*,*) "******************"
+     write(*,*) "Reading stats.def"
+     write(*,*) "******************"
+     do n=1,n_nom_def_max
+          read(99,fmt='(a)',end=88) nom_def(n)
+          write(*,*) 'Output in stats: ', trim(nom_def(n))
+     end do 
+88   continue
+     if (n.ge.n_nom_def_max) then
+        write(*,*)"n_nom_def_max too small in wstats.F90:",n
+        stop
+     end if 
+     n_nom_def=n-1
+     close(99)
+   else
+     firstvar=trim((nom))
+     stats_def=.false.
+     notfoundfirst=.false.
+   endif
+endif
+
+! find firstvar that is in stats.def
+if (notfoundfirst) then
+   do n=1,n_nom_def_max
+     if(trim((nom_def(n))).eq.trim((nom))) then
+        firstvar=trim((nom))
+        notfoundfirst=.false.
+     endif
+   end do 
+endif
+
+! get out of wstats if there is stats.def AND variable not listed
+if (stats_def) then
+   getout=.true.
+   do n=1,n_nom_def
+       if(trim(nom_def(n)).eq.nom) getout=.false.
+   end do
+   if (getout) return
+end if
 
 if (firstvar==nom) then ! If we're back to the first variable
