Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2562)
+++ /trunk/LMDZ.MARS/README	(revision 2563)
@@ -3446,5 +3446,5 @@
 == 07/09/2021 == EM
 Some code cleanup and refactoring around wstats:
-- turn wstats.F90�in a module
+- turn wstats.F90 into a module
 - remove no useless statto_mod.F90
 - incorporate auxilliary inistats and mkstats routines in wstats_mod.F90
@@ -3465,2 +3465,14 @@
  - add tracers ccnco2_h2o_mass_ice, ccnco2_h2o_mass_ccn, ccnco2_h2o_number in place of mem_Nccn_co2, mem_Qccn_co2
  - modification of updaterice_microco2 to take into account water ice as CN and not only dust as CN
+ 
+== 20/09/2021 == EM
+Improve wstats: the callstats flag is now embedded within wstats and checked
+internaly so no need to have some "if (callstats)" conditions around calls
+to wstats anymore.
+In addition: wstats now looks for an (optional) stats.def file in the
+directory where the GCM is run to know which variable should be included
+in the stats.nc file. The stats.def ASCII file should simply contain 
+one variable name per line, in the same way as the diagfi.def file for
+diagfi outputs. If there is no stats.def file then all variables sent to
+wstats will be in the stats.nc file (which matches the behaviour prior to
+this improvement).
Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2562)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2563)
@@ -31,5 +31,5 @@
       use photolysis_mod, only: init_photolysis, nphot
       use iono_h, only: temp_elect
-      use wstats_mod, only: callstats, wstats
+      use wstats_mod, only: wstats
 
       implicit none
@@ -953,5 +953,4 @@
 !            endif
             
-            if (callstats) then
                call wstats(ngrid,'jo3','j o3->o1d',       &
                            's-1',3,jo3_3d(1,1))
@@ -962,5 +961,4 @@
                call wstats(ngrid,'mmean','mean molecular mass',       &
                            'g.mole-1',3,mmean(1,1))
-            endif
          end if ! of if (ngrid.gt.1)
       end if ! of if (output)
Index: /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 2562)
+++ /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 2563)
@@ -9,5 +9,5 @@
       use conc_mod, only: rnew
       use comcstfi_h, only: pi
-      use wstats_mod, only: callstats, wstats
+      use wstats_mod, only: wstats
       implicit none
 
@@ -114,10 +114,9 @@
 ! write diagnostics in micron2/cm3
       
-      if (callstats) then
-        call wstats(ngrid,"surfdust", "Dust surface area",
+      call wstats(ngrid,"surfdust", "Dust surface area",
      $            "micron2 cm-3",3,surfdust*1.e6)
-        call wstats(ngrid,"surfice", "Ice cloud surface area",
+      call wstats(ngrid,"surfice", "Ice cloud surface area",
      $            "micron2 cm-3",3,surfice*1.e6)
-      endif
+
       call writediagfi(ngrid,"surfdust", "Dust surface area",
      $            "micron2 cm-3",3,surfdust*1.e6)
@@ -125,4 +124,3 @@
      $            "micron2 cm-3",3,surfice*1.e6)
 
-      return
       end
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2562)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2563)
@@ -2602,10 +2602,8 @@
 c        WSTATS: Saving statistics
 c        -----------------------------------------------------------------
-c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        ("stats" stores and accumulates key variables in file "stats.nc"
 c        which can later be used to make the statistic files of the run:
-c        "stats")          only possible in 3D runs !
+c        if flag "callstats" from callphys.def is .true.)
          
-       IF (callstats) THEN
-
         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
         call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
@@ -2858,10 +2856,9 @@
            end if ! of if (tracer)
 
-           IF(lastcall) THEN
+           IF(lastcall.and.callstats) THEN
              write (*,*) "Writing stats..."
              call mkstats(ierr)
            ENDIF
 
-         ENDIF !if callstats
 
 c        (Store EOF for Mars Climate database software)
Index: /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90	(revision 2562)
+++ /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90	(revision 2563)
@@ -28,5 +28,6 @@
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
 use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
-                              nbp_lon, nbp_lat, nbp_lev
+                              nbp_lon, nbp_lat, nbp_lev, &
+                              grid_type, unstructured
 implicit none
 
@@ -46,9 +47,16 @@
 integer, dimension(4)  :: id,start,sizes
 logical, save :: firstcall=.TRUE.
-integer :: l,i,j,ig0
 integer,save :: indx
-
 integer, save :: step=0
 !$OMP THREADPRIVATE(firstcall,indx,step)
+integer :: l,i,j,ig0,n
+
+! Added to read an optional stats.def file to select outputs
+logical,save :: stats_def ! .true. if there is a stats.def file
+integer,save :: n_name_def ! number of fields to output in stats.nc
+! NB: stats_def and n_name_def do not need be threadprivate
+integer,parameter :: n_name_def_max=199 ! max number of fields to output
+character(len=120),save :: name_def(n_name_def_max)
+logical :: getout ! to trigger an early exit if variable not in output list
 
 ! Added to work in parallel mode
@@ -68,7 +76,46 @@
 #endif
 
+! 0. Preliminary stuff 
+if (callstats.eqv..false.) then
+  ! exit because creating/writing stats.nc not requested by user
+  return 
+endif
+
+if (grid_type==unstructured) then
+  ! exit because non-structured grid case is not handled
+  return
+endif
+
 ! 1. Initialization (creation of stats.nc file)
 if (firstcall) then
    firstcall=.false.
+
+!$OMP MASTER
+  ! open stats.def definition file if there is one
+  open(99,file="stats.def",status='old',form='formatted',&
+       iostat=ierr)
+  if (ierr.eq.0) then
+    stats_def=.true. ! yes there is a stats.def file
+    write(*,*) "*****************"
+    write(*,*) "Reading stats.def"
+    write(*,*) "*****************"
+    do n=1,n_name_def_max
+      read(99,fmt='(a)',end=88) name_def(n)
+      write(*,*) 'Output in stats: ', trim(name_def(n))
+    enddo
+88  continue
+    ! check there is no overflow
+    if (n.ge.n_name_def_max) then
+      write(*,*) "n_name_def_max too small in wstats:",n
+      call abort_physic("wstats","n_name_def_max too small",1)
+    endif
+    n_name_def=n-1
+    close(99)
+  else
+    stats_def=.false. ! no stats.def file; output all fields sent to wstats
+  endif ! of if (ierr.eq.0)
+!$OMP END MASTER
+!$OMP BARRIER
+
    firstvar=trim((nom))
    call inistats(ierr)
@@ -88,5 +135,5 @@
      allocate(dx2(1,1))
    endif
-endif
+endif ! of if (firstcall)
 
 if (firstvar==nom) then ! If we're back to the first variable, increment time counter
@@ -98,4 +145,21 @@
    RETURN
 endif
+
+! Exit if there is a stats.def file and the variable is not in the list
+if (stats_def) then
+  getout=.true.
+  do n=1,n_name_def
+    ! look for the variable's name in the list
+    if (trim(name_def(n)).eq.nom) then
+      getout=.false.
+      ! found it, no need to scan the rest of the list exit loop
+      exit
+    endif
+  enddo
+  if (getout) then
+    ! variable not in the list so exit routine now
+    return
+  endif
+endif ! of if (stats_def) 
 
 ! collect fields on a global physics grid
