Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2558)
+++ /trunk/LMDZ.MARS/README	(revision 2559)
@@ -2466,5 +2466,5 @@
 microphysics (just like water clouds). The amplitude is spantCO2, also read from callphys.def
 - The previous logical microphysco2 has been removed.
-- Cloud opacity at 1µm is now computed in the co2cloud.F routine
+- Cloud opacity at 1Âµm is now computed in the co2cloud.F routine
 - Most of the co2 ice clouds scheme writediagfi are now in co2clouds.F
 - Some cleaning of the CO2 ice clouds routine has been done. Not perfect yet!
@@ -2482,5 +2482,5 @@
   if a gravity wave propagates to the mesosphere or not (wheter the wave saturates on its way or not). 
   The sub-grid temperature distribution is cancelled if the saturation index value of the column (between 12 and 80 km) is > 0.1, and is applied (+/- 3K) otherwise. if the keyword satindexco2=.true. in callphys.def (only applies if CLFvaryingCO2=.true. anyway). If set to .false., deactivate this filter and the sub grid T distribution applies everywhere and everytime. See comments in co2clouds.F
-- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni Määttänen have the files. 
+- All that you need to launch a GCM run with co2 clouds is described in co2clouds.F comments. Ehouarn Millour and Anni MÃ¤Ã¤ttÃ¤nen have the files. 
  
 == 15/11/2017 == EM
@@ -3443,2 +3443,9 @@
 (between levels with no missing values in both files, which vary with space and time) and outputs
 the ratio tau_GCM/tau_MCS to be used to normalize the GCM dust opacity profiles when comparing to MCS
+
+== 07/09/2021 == EM
+Some code cleanup and refactoring around wstats:
+- turn wstats.F90 in a module
+- remove no useless statto_mod.F90
+- incorporate auxilliary inistats and mkstats routines in wstats_mod.F90
+- move flag "callstats" from callkeys.h to being a module variable of wstats_mod
Index: /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2558)
+++ /trunk/LMDZ.MARS/libf/aeronomars/calchim_mod.F90	(revision 2559)
@@ -31,4 +31,5 @@
       use photolysis_mod, only: init_photolysis, nphot
       use iono_h, only: temp_elect
+      use wstats_mod, only: callstats, wstats
 
       implicit none
Index: /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 2558)
+++ /trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F	(revision 2559)
@@ -8,5 +8,6 @@
      &                      igcm_ccn_number, varian, ccn_factor
       use conc_mod, only: rnew
-      USE comcstfi_h
+      use comcstfi_h, only: pi
+      use wstats_mod, only: callstats, wstats
       implicit none
 
Index: /trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2558)
+++ /trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 2559)
@@ -7,5 +7,5 @@
 
       COMMON/callkeys_l/callrad,calldifv,calladj,callcond,callsoil      &
-     &   ,season,diurnal,lwrite,calllott,callstats,calleofdump          &
+     &   ,season,diurnal,lwrite,calllott,calleofdump                    &
      &   ,callnirco2,callnlte,callthermos,callconduct,calleuv           &
      &   ,callmolvis,callmoldiff,thermochem,thermoswater,callemis       &
@@ -28,5 +28,5 @@
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
      &   season,diurnal,lwrite,calllott,calllott_nonoro                 &
-     &   ,callstats,calleofdump                                         &
+     &   ,calleofdump                                                   &
      &   ,callnirco2,callnlte,callthermos,callconduct,                  &
      &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2558)
+++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2559)
@@ -43,4 +43,5 @@
      &                      ini_scatterers,tauvis
       use datafile_mod, only: datadir
+      use wstats_mod, only: callstats
       use calchim_mod, only: ichemistry
       use co2condens_mod, only: scavco2cond
Index: unk/LMDZ.MARS/libf/phymars/inistats.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 2558)
+++ 	(revision )
@@ -1,145 +1,0 @@
-      subroutine inistats(ierr)
-
-      use statto_mod, only: istats,istime
-      use mod_phys_lmdz_para, only : is_master
-      USE vertical_layers_mod, ONLY: ap,bp,aps,bps,preff,
-     &                               pseudoalt,presnivs
-      USE nrtype, ONLY: pi
-      USE time_phylmdz_mod, ONLY: daysec,dtphys
-      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
-      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
-      implicit none
-
-      include "netcdf.inc"
-
-      integer,intent(out) :: ierr
-      integer :: nid
-      integer :: l,nsteppd
-      real, dimension(nbp_lev) ::  sig_s
-      real,allocatable :: lon_reg_ext(:) ! extended longitudes
-      integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
-      real, dimension(istime) :: lt
-      integer :: nvarid
-
-
-      IF (nbp_lon*nbp_lat==1) THEN
-        ! 1D model
-        ALLOCATE(lon_reg_ext(1))
-      ELSE
-        ! 3D model
-        ALLOCATE(lon_reg_ext(nbp_lon+1))
-      ENDIF
-      
-      write (*,*) 
-      write (*,*) '                        || STATS ||'
-      write (*,*) 
-      write (*,*) 'daysec',daysec
-      write (*,*) 'dtphys',dtphys
-      nsteppd=nint(daysec/dtphys)
-      write (*,*) 'nsteppd=',nsteppd
-      if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) then
-        call abort_physic("inistats",
-     &                    '1 sol .ne. n physics steps',1)
-      endif
-
-      if(mod(nsteppd,istime).ne.0) then
-        call abort_physic("inistats",
-     &                    '1 sol .ne. n*istime physics steps',1)
-      endif
-
-      istats=nsteppd/istime
-      write (*,*) 'istats=',istats
-      write (*,*) 'Storing ',istime,'times per day'
-      write (*,*) 'thus every ',istats,'physical timestep '
-      write (*,*) 
-
-      do l= 1, nbp_lev
-         sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
-         pseudoalt(l)=-10.*log(presnivs(l)/preff)   
-      enddo
-      
-      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
-      IF (nbp_lon*nbp_lat/=1) THEN
-        ! In 3D, add extra redundant point (180 degrees,
-        ! since lon_reg starts at -180)
-        lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
-      ENDIF
-
-      if (is_master) then
-      ! only the master needs do this
-
-      ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
-      if (ierr.ne.NF_NOERR) then
-         write (*,*) NF_STRERROR(ierr)
-        call abort_physic("inistats",
-     &                    'failed creating stats.nc',1)
-      endif
-
-      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
-      IF (nbp_lon*nbp_lat==1) THEN
-        ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
-      ELSE
-        ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
-      ENDIF
-      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
-      ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
-      ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
-
-      ierr = NF_ENDDEF(nid)
-      call def_var_stats(nid,"Time","Time",
-     &            "days since 0000-00-0 00:00:00",1,
-     &            idim_time,nvarid,ierr)
-! Time is initialised later by mkstats subroutine
-
-      call def_var_stats(nid,"latitude","latitude",
-     &            "degrees_north",1,idim_lat,nvarid,ierr)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
-#endif
-      call def_var_stats(nid,"longitude","East longitude",
-     &            "degrees_east",1,idim_lon,nvarid,ierr)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
-#endif
-
-! Niveaux verticaux, aps et bps
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
-#endif 
-      call def_var_stats(nid,"aps","hybrid pressure at midlayers"
-     & ," ",1,idim_llm,nvarid,ierr)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
-#endif
-
-      call def_var_stats(nid,"bps","hybrid sigma at midlayers"
-     & ," ",1,idim_llm,nvarid,ierr)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
-#else 
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
-#endif
-
-      ierr=NF_CLOSE(nid)
-
-      endif ! of if (is_master)
-      end subroutine inistats
-
Index: unk/LMDZ.MARS/libf/phymars/mkstat.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/mkstat.F90	(revision 2558)
+++ 	(revision )
@@ -1,197 +1,0 @@
-subroutine mkstats(ierr)
-
- 
-!
-!  This program writes a stats.nc file from sums and sums of squares
-!  to means and standard deviations and also writes netcdf style
-!  file so that the data can be viewed easily.  The data file is
-!  overwritten in place.  
-!  SRL  21 May 1996
-!  Yann W. july 2003
-
-use statto_mod, only: istime,count
-use mod_phys_lmdz_para, only : is_master
-use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
-
-implicit none
-
-include "netcdf.inc"
-
-integer :: ierr,nid,nbvar,i,ndims,lt,nvarid
-integer, dimension(4) :: id,varid,start,size
-integer, dimension(5) :: dimids
-character (len=50) :: name,nameout,units,title
-real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
-real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
-real, dimension(istime) :: time
-real, dimension(nbp_lat) :: lat
-real,allocatable :: lon(:)
-real, dimension(nbp_lev) :: alt
-logical :: lcopy=.true.
-!integer :: latid,lonid,altid,timeid
-integer :: meanid,sdid
-!integer, dimension(4) :: dimout
-
-! Incrementation of count for the last step, which is not done in wstats
-count(istime)=count(istime)+1
-
-if (is_master) then
-! only the master needs do this
-if (klon_glo>1) then
-  allocate(lon(nbp_lon+1))
-  allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
-  allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
-  allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
-  allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
-  allocate(sum2d(nbp_lon+1,nbp_lat))
-  allocate(square2d(nbp_lon+1,nbp_lat))
-  allocate(mean2d(nbp_lon+1,nbp_lat))
-  allocate(sd2d(nbp_lon+1,nbp_lat))
-else ! 1D model case
-  allocate(lon(1))
-  allocate(sum3d(1,1,nbp_lev))
-  allocate(square3d(1,1,nbp_lev))
-  allocate(mean3d(1,1,nbp_lev))
-  allocate(sd3d(1,1,nbp_lev))
-  allocate(sum2d(1,1))
-  allocate(square2d(1,1))
-  allocate(mean2d(1,1))
-  allocate(sd2d(1,1))
-endif
-
-ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
-
-! We catch the id of dimensions of the stats file
-
-ierr= NF_INQ_DIMID(nid,"latitude",id(1))
-ierr= NF_INQ_DIMID(nid,"longitude",id(2))
-ierr= NF_INQ_DIMID(nid,"altitude",id(3))
-ierr= NF_INQ_DIMID(nid,"Time",id(4))
-
-ierr= NF_INQ_VARID(nid,"latitude",varid(1))
-ierr= NF_INQ_VARID(nid,"longitude",varid(2))
-ierr= NF_INQ_VARID(nid,"altitude",varid(3))
-ierr= NF_INQ_VARID(nid,"Time",varid(4))
-
-! Time initialisation
-
-do i=1,istime
-   time(i)=i*24./istime
-#ifdef NC_DOUBLE
-   ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),i,1,time(i))
-#else
-   ierr= NF_PUT_VARA_REAL(nid,varid(4),i,1,time(i))
-#endif
-enddo
-
-! We catche the values of the variables
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
-         ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
-         ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
-#else
-         ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
-         ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
-         ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
-#endif
-
-! We catch the number of variables in the stats file
-ierr = NF_INQ_NVARS(nid,nbvar)
-
-! to catche the "real" number of variables (without the "additionnal variables")
-nbvar=(nbvar-4)/2 
-
-do i=1,nbvar
-   varid=(i-1)*2+5
-
-   ! What's the variable's name?
-   ierr=NF_INQ_VARNAME(nid,varid,name)
-   write(*,*) "OK variable ",name
-   ! Its units?
-   units=" "
-   ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
-   ! Its title?
-   title=" "
-   ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
-   ! Its number of dimensions?   
-   ierr=NF_INQ_VARNDIMS(nid,varid,ndims)
-   ! Its values?
-
-   if(ndims==4) then ! lat, lon, alt & time
-
-!      dimout(1)=lonid
-!      dimout(2)=latid
-!      dimout(3)=altid
-!      dimout(4)=timeid
-
-      if (klon_glo>1) then ! general case
-        size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
-      else ! 1D model
-        size=(/1,1,nbp_lev,1/)
-      endif
-      do lt=1,istime
-         start=(/1,1,1,lt/)
-         ! Extraction of the "source" variables
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum3d)
-         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square3d)
-#else
-         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum3d)
-         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square3d)
-#endif
-         ! Calculation of these variables
-         mean3d=sum3d/count(lt)
-         sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
-         ! Writing of the variables
-#ifdef NC_DOUBLE
-         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean3d)
-         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd3d)
-#else
-         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean3d)
-         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd3d)
-#endif
-      enddo
-
-    else if (ndims.eq.3) then
-
-!      dimout(1)=lonid
-!      dimout(2)=latid
-!      dimout(3)=timeid
-
-      if (klon_glo>1) then ! general case
-        size=(/nbp_lon+1,nbp_lat,1,0/)
-      else
-        size=(/1,1,1,0/)
-      endif
-      do lt=1,istime
-         start=(/1,1,lt,0/)
-         ! Extraction of the "source" variables
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum2d)
-         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square2d)
-#else
-         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum2d)
-         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square2d)
-#endif
-         ! Calculation of these variables
-         mean2d=sum2d/count(lt)
-         sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
-         ! Writing of the variables
-#ifdef NC_DOUBLE
-         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean2d)
-         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd2d)
-#else
-         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean2d)
-         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd2d)
-#endif
-      enddo
-
-    endif 
-enddo
-
-ierr= NF_CLOSE(nid)
-
-endif ! of if (is_master)
-
-end
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2558)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2559)
@@ -84,4 +84,5 @@
       use phyredem, only: physdem0, physdem1
       use phyetat0_mod, only: phyetat0, tab_cntrl_mod
+      use wstats_mod, only: callstats, wstats, mkstats
       use eofdump_mod, only: eofdump
       USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
Index: unk/LMDZ.MARS/libf/phymars/statto_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/statto_mod.F90	(revision 2558)
+++ 	(revision )
@@ -1,35 +1,0 @@
-MODULE statto_mod
-IMPLICIT NONE  
-
-!  statto_mod:
-!     This include file controls the production of statistics.
-!     Some variables could be set in a namelist, but it is easier to
-!     do it here since arrays can then be dimensioned using parameters
-!     and values shouldn't have to change too often.   SRL
-
-!     Calculate stats every istats physics timesteps, starting at first
-!     call.  If istats=0 then don't do statistics at all.  Check value
-!     if number of physics timesteps changes.
-	integer istats
-
-!     Calculate itime independent sums and sums of squares,
-!     example, istat=1,istime=1 gives a single time mean
-	integer, parameter :: istime=12
-
-!     Number of 2D and 3D variables on which to do statistics.
-	integer n2dvar, n3dvar
-	parameter (n2dvar = 8, n3dvar = 5)
-
-!     Units for writing stats header and data
-	integer usdata
-
-!     count tab to know the variable record
-        integer count(istime)
-
-!     Record of the number of stores made for each time.
-	integer nstore(istime)
-
-! Size of the "controle" array
-        integer, parameter :: cntrlsize=15
-
-END MODULE statto_mod
Index: unk/LMDZ.MARS/libf/phymars/wstats.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 2558)
+++ 	(revision )
@@ -1,517 +1,0 @@
-subroutine wstats(ngrid,nom,titre,unite,dim,px)
-
-use statto_mod, only: istats,istime,count
-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
-implicit none
-
-include "netcdf.inc"
-
-integer,intent(in) :: ngrid
-character (len=*),intent(in) :: nom,titre,unite
-integer,intent(in) :: dim
-real,intent(in) :: px(ngrid,nbp_lev)
-real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
-real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
-character (len=50) :: namebis
-character (len=50), save :: firstvar
-!$OMP THREADPRIVATE(firstvar)
-integer :: ierr,varid,nbdim,nid
-integer :: meanid,sdid
-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)
-
-! Added to work in parallel mode
-#ifdef CPP_PARA
-real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
-real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
-real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
-real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
-real px2(ngrid)
-real px3(ngrid,nbp_lev)
-#else
-! When not running in parallel mode:
-real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
-real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
-real px2_glop(ngrid) ! to store a 2D data set on global physics grid
-real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
-#endif
-
-! 1. Initialization (creation of stats.nc file)
-if (firstcall) then
-   firstcall=.false.
-   firstvar=trim((nom))
-   call inistats(ierr)
-   if (klon_glo>1) then ! general case, 3D GCM
-     allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
-     allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
-     allocate(dx3(nbp_lon+1,nbp_lat,nbp_lev))
-     allocate(mean2d(nbp_lon+1,nbp_lat))
-     allocate(sd2d(nbp_lon+1,nbp_lat))
-     allocate(dx2(nbp_lon+1,nbp_lat))
-   else ! 1D model
-     allocate(mean3d(1,1,nbp_lev))
-     allocate(sd3d(1,1,nbp_lev))
-     allocate(dx3(1,1,nbp_lev))
-     allocate(mean2d(1,1))
-     allocate(sd2d(1,1))
-     allocate(dx2(1,1))
-   endif
-endif
-
-if (firstvar==nom) then ! If we're back to the first variable, increment time counter
-      step = step + 1
-endif
-
-if (mod(step,istats).ne.0) then
-  ! if its not time to write to file, exit
-   RETURN
-endif
-
-! collect fields on a global physics grid
-#ifdef CPP_PARA
- if (dim.eq.3) then
-  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
-  ! Gather fieds on a "global" (without redundant longitude) array
-  call Gather(px3,px3_glop)
-!$OMP MASTER
-  if (is_mpi_root) then
-    call Grid1Dto2D_glo(px3_glop,px3_glo)
-    ! copy dx3_glo() to dx3(:) and add redundant longitude
-    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
-    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
-  endif
-!$OMP END MASTER
-!$OMP BARRIER
- else ! dim.eq.2
-  ! Gather fieds on a "global" (without redundant longitude) array
-  px2(:)=px(:,1)
-  call Gather(px2,px2_glop)
-!$OMP MASTER
-          if (is_mpi_root) then
-            call Grid1Dto2D_glo(px2_glop,px2_glo)
-            ! copy px2_glo() to dx2(:) and add redundant longitude
-            dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
-            dx2(nbp_lon+1,:)=dx2(1,:)
-          endif
-!$OMP END MASTER
-!$OMP BARRIER
- endif
-#else
-  if (dim.eq.3) then
-    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
-!  Passage variable physique -->  variable dynamique
-    DO l=1,nbp_lev
-      DO i=1,nbp_lon
-         px3_glo(i,1,l)=px(1,l)
-         px3_glo(i,nbp_lat,l)=px(ngrid,l)
-      ENDDO
-      DO j=2,nbp_lat-1
-         ig0= 1+(j-2)*nbp_lon
-         DO i=1,nbp_lon
-            px3_glo(i,j,l)=px(ig0+i,l)
-         ENDDO
-      ENDDO
-    ENDDO
-  else ! dim.eq.2
-    px2_glop(:)=px(:,1)
-!    Passage variable physique -->  physique dynamique
-   DO i=1,nbp_lon
-     px2_glo(i,1)=px(1,1)
-     px2_glo(i,nbp_lat)=px(ngrid,1)
-   ENDDO
-   DO j=2,nbp_lat-1
-     ig0= 1+(j-2)*nbp_lon
-     DO i=1,nbp_lon
-        px2_glo(i,j)=px(ig0+i,1)
-     ENDDO
-   ENDDO
-  endif
-#endif
-
-! 2. Write field to file
-
-if (is_master) then
-! only master needs do this
-
-ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
-
-namebis=trim(nom)
-
-! test: check if that variable already exists in file
-ierr= NF_INQ_VARID(nid,namebis,meanid)
-
-if (ierr.ne.NF_NOERR) then
-  ! variable not in file, create/define it
-   if (firstvar==nom) then 
-      indx=1
-      count(:)=0
-   endif
-
-
-!declaration de la variable
-
-! choix du nom des coordonnees
-   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
-   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
-   if (dim.eq.3) then
-      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
-      ierr= NF_INQ_DIMID(nid,"Time",id(4))
-      nbdim=4
-   else if (dim==2) then
-      ierr= NF_INQ_DIMID(nid,"Time",id(3))
-      nbdim=3
-   endif
-
-   write (*,*) "====================="
-   write (*,*) "STATS: creating ",nom
-   namebis=trim(nom)
-   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
-   if (dim.eq.3) then
-     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
-   else ! dim.eq.2
-     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
-   endif
-   namebis=trim(nom)//"_sd"
-   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
-   if (dim.eq.3) then
-     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
-   else ! dim.eq.2
-     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
-   endif
-
-   ierr= NF_CLOSE(nid)
-   return
-
-else
-   ! variable found in file
-   namebis=trim(nom)//"_sd"
-   ierr= NF_INQ_VARID(nid,namebis,sdid)
-
-endif
-
-if (firstvar==nom) then 
-   count(indx)=count(int(indx))+1
-   indx=indx+1
-   if (indx>istime) then
-      indx=1
-   endif
-endif
-
-if (count(indx)==0) then
-   ! very first time we write the variable to file
-   if (dim.eq.3) then
-      start=(/1,1,1,indx/)
-      if (klon_glo>1) then !general case
-        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
-      else
-        sizes=(/1,1,nbp_lev,1/)
-      endif
-      mean3d(:,:,:)=0
-      sd3d(:,:,:)=0
-   else if (dim.eq.2) then
-      start=(/1,1,indx,0/)
-      if (klon_glo>1) then !general case
-        sizes=(/nbp_lon+1,nbp_lat,1,0/)
-      else
-        sizes=(/1,1,1,0/)
-      endif
-      mean2d(:,:)=0
-      sd2d(:,:)=0
-   endif
-else
-   ! load values from file
-   if (dim.eq.3) then
-      start=(/1,1,1,indx/)
-      if (klon_glo>1) then !general case
-        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
-      else ! 1D model case
-        sizes=(/1,1,nbp_lev,1/)
-      endif
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
-      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
-#else
-      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
-      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
-#endif
-      if (ierr.ne.NF_NOERR) then
-         write (*,*) "wstats error reading :",trim(nom)
-         write (*,*) NF_STRERROR(ierr)
-         call abort_physic("wstats","Failed reading "//trim(nom),1)
-      endif
-
-   else if (dim.eq.2) then
-      start=(/1,1,indx,0/)
-      if (klon_glo>1) then ! general case
-        sizes=(/nbp_lon+1,nbp_lat,1,0/)
-      else
-        sizes=(/1,1,1,0/)
-      endif
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
-      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
-#else
-      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
-      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
-#endif
-      if (ierr.ne.NF_NOERR) then
-         write (*,*) "wstats error reading :",trim(nom)
-         write (*,*) NF_STRERROR(ierr)
-         call abort_physic("wstats","Failed reading "//trim(nom),1)
-      endif
-   endif
-endif ! of if (count(indx)==0)
-
-
-! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
-
-if (dim.eq.3) then
-  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
-  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
-    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
-  ENDIF
-else ! dim.eq.2
-  dx2(1:nbp_lon,:)=px2_glo(:,:)
-  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
-    dx2(nbp_lon+1,:)=dx2(1,:)
-  ENDIF
-endif
-
-
-! 2.2. Add current values to previously stored sums
-
-if (dim.eq.3) then
-
-   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
-   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
-   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
-   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
-#endif
-  if (ierr.ne.NF_NOERR) then
-     write (*,*) "wstats error writing :",trim(nom)
-     write (*,*) NF_STRERROR(ierr)
-     call abort_physic("wstats","Failed writing "//trim(nom),1)
-  endif
-
-else if (dim.eq.2) then
-
-   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
-   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
-   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
-   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
-#endif
-  if (ierr.ne.NF_NOERR) then
-     write (*,*) "wstats error writing :",trim(nom)
-     write(*,*) "start:",start
-     write(*,*) "sizes:",sizes
-     write(*,*) "mean2d:",mean2d
-     write(*,*) "sd2d:",sd2d
-     write (*,*) NF_STRERROR(ierr)
-     call abort_physic("wstats","Failed writing "//trim(nom),1)
-  endif
-
-endif ! of if (dim.eq.3) elseif (dim.eq.2)
-
-  ierr= NF_CLOSE(nid)
-endif ! of if (is_master)
-
-end subroutine wstats
-
-!======================================================
-subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
-use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
-
-implicit none
-
-include "netcdf.inc"
-
-integer, intent(in) :: nid,varid,dim,indx,ngrid
-real, dimension(ngrid,nbp_lev), intent(in) :: px
-integer, intent(out) :: ierr
-
-integer :: l,i,j,ig0
-integer, dimension(4) :: start,sizes
-real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
-real, dimension(nbp_lon+1,nbp_lat) :: dx2
-real :: dx3_1d(nbp_lev) ! for 1D outputs
-real :: dx2_1d ! for 1D outputs
-
-if (dim.eq.3) then
-
-   start=(/1,1,1,indx/)
-   if (klon_glo>1) then ! general 3D case
-     sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
-   else
-     sizes=(/1,1,nbp_lev,1/)
-   endif
-
-!  Passage variable physique -->  variable dynamique
-
-   if (klon_glo>1) then ! general case
-    DO l=1,nbp_lev
-      DO i=1,nbp_lon+1
-         dx3(i,1,l)=px(1,l)
-         dx3(i,nbp_lat,l)=px(ngrid,l)
-      ENDDO
-      DO j=2,nbp_lat-1
-         ig0= 1+(j-2)*nbp_lon
-         DO i=1,nbp_lon
-            dx3(i,j,l)=px(ig0+i,l)
-         ENDDO
-         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
-      ENDDO
-    ENDDO
-   else ! 1D model case
-     dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
-   endif
-
-#ifdef NC_DOUBLE
-   if (klon_glo>1) then
-     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
-   else
-     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3_1d)
-   endif
-#else
-   if (klon_glo>1) then
-     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
-   else
-     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3_1d)
-   endif
-#endif
-  if (ierr.ne.NF_NOERR) then
-     write (*,*) "inivar error writing variable"
-     write (*,*) NF_STRERROR(ierr)
-     call abort_physic("inivar","error writing variable",1)
-  endif
-
-else if (dim.eq.2) then
-
-      start=(/1,1,indx,0/)
-      if (klon_glo>1) then ! general 3D case
-        sizes=(/nbp_lon+1,nbp_lat,1,0/)
-      else
-        sizes=(/1,1,1,0/)
-      endif
-
-!    Passage variable physique -->  physique dynamique
-
-  if (klon_glo>1) then ! general case
-  DO i=1,nbp_lon+1
-     dx2(i,1)=px(1,1)
-     dx2(i,nbp_lat)=px(ngrid,1)
-  ENDDO
-  DO j=2,nbp_lat-1
-     ig0= 1+(j-2)*nbp_lon
-     DO i=1,nbp_lon
-        dx2(i,j)=px(ig0+i,1)
-     ENDDO
-     dx2(nbp_lon+1,j)=dx2(1,j)
-  ENDDO
-  else ! 1D model case
-    dx2_1d=px(1,1)
-  endif
-  
-#ifdef NC_DOUBLE
-   if (klon_glo>1) then
-     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
-   else
-     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2_1d)
-   endif
-#else
-   if (klon_glo>1) then
-     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
-   else
-     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2_1d)
-   endif
-#endif
-  if (ierr.ne.NF_NOERR) then
-     write (*,*) "inivar error writing variable"
-     write (*,*) NF_STRERROR(ierr)
-     call abort_physic("inivar","error writing variable",1)
-  endif
-
-endif
-
-end subroutine inivar
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
-
-! This subroutine defines variable 'name' in a (pre-existing and opened)
-! NetCDF file (known from its NetCDF ID 'nid').
-! The number of dimensions 'nbdim' of the variable, as well as the IDs of
-! corresponding dimensions must be set (in array 'dimids').
-! Upon successfull definition of the variable, 'nvarid' contains the
-! NetCDF ID of the variable.
-! The variables' attributes 'title' (Note that 'long_name' would be more
-! appropriate) and 'units' are also set. 
-
-implicit none
-
-include "netcdf.inc"
-
-integer,intent(in) :: nid ! NetCDF file ID
-character(len=*),intent(in) :: name ! the variable's name
-character(len=*),intent(in) :: title ! 'title' attribute of variable
-character(len=*),intent(in) :: units ! 'units' attribute of variable
-integer,intent(in) :: nbdim ! number of dimensions of the variable
-integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
-                                              ! the variable is defined along
-integer,intent(out) :: nvarid ! NetCDF ID of the variable
-integer,intent(out) :: ierr ! returned NetCDF staus code
-
-! 1. Switch to NetCDF define mode 
-ierr=NF_REDEF(nid)
-
-! 2. Define the variable
-#ifdef NC_DOUBLE
-ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
-#else
-ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
-#endif
-if(ierr/=NF_NOERR) then
-   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
-   write(*,*) NF_STRERROR(ierr)
-   call abort_physic("def_var_stats","Failed defining "//trim(name),1)
-endif
-
-! 3. Write attributes
-ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
-                     len_trim(adjustl(title)),adjustl(title))
-if(ierr/=NF_NOERR) then
-   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
-   write(*,*) NF_STRERROR(ierr)
-   call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
-endif
-
-ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
-                     len_trim(adjustl(units)),adjustl(units))
-if(ierr/=NF_NOERR) then
-   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
-   write(*,*) NF_STRERROR(ierr)
-   call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
-endif
-
-! 4. Switch out of NetCDF define mode
-ierr = NF_ENDDEF(nid)
-
-end subroutine def_var_stats
-
Index: /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90	(revision 2559)
+++ /trunk/LMDZ.MARS/libf/phymars/wstats_mod.F90	(revision 2559)
@@ -0,0 +1,892 @@
+module wstats_mod
+
+! module containing parameters and routines to generate the "stats.nc" file
+! which will contain a mean statistical day of variables extracted at set
+! times of day every day of the simulation, and the standard deviations thereof
+
+implicit none
+
+logical,save :: callstats ! global flag to trigger generating a stats.nc
+                          ! file or not. Initialized in conf_gcm()
+!$OMP THREADPRIVATE(callstats)
+
+integer,save :: istats ! calculate stats every istats physics timestep,
+                       ! starting at first call. Initialized by inistats()
+!$OMP THREADPRIVATE(istats)
+
+integer,parameter :: istime=12 ! number of time steps per sol to store
+
+integer,save :: count(istime) ! count tab to know the variable record
+!$OMP THREADPRIVATE(istats,count)
+
+contains
+
+subroutine wstats(ngrid,nom,titre,unite,dim,px)
+
+! Main routine to call to trigger writing a given field to the "stats.nc" file 
+
+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
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(in) :: ngrid
+character (len=*),intent(in) :: nom,titre,unite
+integer,intent(in) :: dim
+real,intent(in) :: px(ngrid,nbp_lev)
+real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
+real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
+character (len=50) :: namebis
+character (len=50), save :: firstvar
+!$OMP THREADPRIVATE(firstvar)
+integer :: ierr,varid,nbdim,nid
+integer :: meanid,sdid
+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)
+
+! Added to work in parallel mode
+#ifdef CPP_PARA
+real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
+real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
+real px2(ngrid)
+real px3(ngrid,nbp_lev)
+#else
+! When not running in parallel mode:
+real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
+real px2_glop(ngrid) ! to store a 2D data set on global physics grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
+#endif
+
+! 1. Initialization (creation of stats.nc file)
+if (firstcall) then
+   firstcall=.false.
+   firstvar=trim((nom))
+   call inistats(ierr)
+   if (klon_glo>1) then ! general case, 3D GCM
+     allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(dx3(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(mean2d(nbp_lon+1,nbp_lat))
+     allocate(sd2d(nbp_lon+1,nbp_lat))
+     allocate(dx2(nbp_lon+1,nbp_lat))
+   else ! 1D model
+     allocate(mean3d(1,1,nbp_lev))
+     allocate(sd3d(1,1,nbp_lev))
+     allocate(dx3(1,1,nbp_lev))
+     allocate(mean2d(1,1))
+     allocate(sd2d(1,1))
+     allocate(dx2(1,1))
+   endif
+endif
+
+if (firstvar==nom) then ! If we're back to the first variable, increment time counter
+      step = step + 1
+endif
+
+if (mod(step,istats).ne.0) then
+  ! if its not time to write to file, exit
+   RETURN
+endif
+
+! collect fields on a global physics grid
+#ifdef CPP_PARA
+ if (dim.eq.3) then
+  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
+  ! Gather fieds on a "global" (without redundant longitude) array
+  call Gather(px3,px3_glop)
+!$OMP MASTER
+  if (is_mpi_root) then
+    call Grid1Dto2D_glo(px3_glop,px3_glo)
+    ! copy dx3_glo() to dx3(:) and add redundant longitude
+    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
+  endif
+!$OMP END MASTER
+!$OMP BARRIER
+ else ! dim.eq.2
+  ! Gather fieds on a "global" (without redundant longitude) array
+  px2(:)=px(:,1)
+  call Gather(px2,px2_glop)
+!$OMP MASTER
+          if (is_mpi_root) then
+            call Grid1Dto2D_glo(px2_glop,px2_glo)
+            ! copy px2_glo() to dx2(:) and add redundant longitude
+            dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
+          endif
+!$OMP END MASTER
+!$OMP BARRIER
+ endif
+#else
+  if (dim.eq.3) then
+    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
+!  Passage variable physique -->  variable dynamique
+    DO l=1,nbp_lev
+      DO i=1,nbp_lon
+         px3_glo(i,1,l)=px(1,l)
+         px3_glo(i,nbp_lat,l)=px(ngrid,l)
+      ENDDO
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
+            px3_glo(i,j,l)=px(ig0+i,l)
+         ENDDO
+      ENDDO
+    ENDDO
+  else ! dim.eq.2
+    px2_glop(:)=px(:,1)
+!    Passage variable physique -->  physique dynamique
+   DO i=1,nbp_lon
+     px2_glo(i,1)=px(1,1)
+     px2_glo(i,nbp_lat)=px(ngrid,1)
+   ENDDO
+   DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
+        px2_glo(i,j)=px(ig0+i,1)
+     ENDDO
+   ENDDO
+  endif
+#endif
+
+! 2. Write field to file
+
+if (is_master) then
+! only master needs do this
+
+ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
+
+namebis=trim(nom)
+
+! test: check if that variable already exists in file
+ierr= NF_INQ_VARID(nid,namebis,meanid)
+
+if (ierr.ne.NF_NOERR) then
+  ! variable not in file, create/define it
+   if (firstvar==nom) then 
+      indx=1
+      count(:)=0
+   endif
+
+
+!declaration de la variable
+
+! choix du nom des coordonnees
+   ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+   ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+   if (dim.eq.3) then
+      ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+      ierr= NF_INQ_DIMID(nid,"Time",id(4))
+      nbdim=4
+   else if (dim==2) then
+      ierr= NF_INQ_DIMID(nid,"Time",id(3))
+      nbdim=3
+   endif
+
+   write (*,*) "====================="
+   write (*,*) "STATS: creating ",nom
+   namebis=trim(nom)
+   call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
+   if (dim.eq.3) then
+     call inivar(nid,meanid,size(px3_glop,1),dim,indx,px3_glop,ierr)
+   else ! dim.eq.2
+     call inivar(nid,meanid,size(px2_glop,1),dim,indx,px2_glop,ierr)
+   endif
+   namebis=trim(nom)//"_sd"
+   call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
+   if (dim.eq.3) then
+     call inivar(nid,sdid,size(px3_glop,1),dim,indx,px3_glop,ierr)
+   else ! dim.eq.2
+     call inivar(nid,sdid,size(px2_glop,1),dim,indx,px2_glop,ierr)
+   endif
+
+   ierr= NF_CLOSE(nid)
+   return
+
+else
+   ! variable found in file
+   namebis=trim(nom)//"_sd"
+   ierr= NF_INQ_VARID(nid,namebis,sdid)
+
+endif
+
+if (firstvar==nom) then 
+   count(indx)=count(int(indx))+1
+   indx=indx+1
+   if (indx>istime) then
+      indx=1
+   endif
+endif
+
+if (count(indx)==0) then
+   ! very first time we write the variable to file
+   if (dim.eq.3) then
+      start=(/1,1,1,indx/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else
+        sizes=(/1,1,nbp_lev,1/)
+      endif
+      mean3d(:,:,:)=0
+      sd3d(:,:,:)=0
+   else if (dim.eq.2) then
+      start=(/1,1,indx,0/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
+      mean2d(:,:)=0
+      sd2d(:,:)=0
+   endif
+else
+   ! load values from file
+   if (dim.eq.3) then
+      start=(/1,1,1,indx/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else ! 1D model case
+        sizes=(/1,1,nbp_lev,1/)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean3d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd3d)
+#endif
+      if (ierr.ne.NF_NOERR) then
+         write (*,*) "wstats error reading :",trim(nom)
+         write (*,*) NF_STRERROR(ierr)
+         call abort_physic("wstats","Failed reading "//trim(nom),1)
+      endif
+
+   else if (dim.eq.2) then
+      start=(/1,1,indx,0/)
+      if (klon_glo>1) then ! general case
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
+      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
+#else
+      ierr = NF_GET_VARA_REAL(nid,meanid,start,sizes,mean2d)
+      ierr = NF_GET_VARA_REAL(nid,sdid,start,sizes,sd2d)
+#endif
+      if (ierr.ne.NF_NOERR) then
+         write (*,*) "wstats error reading :",trim(nom)
+         write (*,*) NF_STRERROR(ierr)
+         call abort_physic("wstats","Failed reading "//trim(nom),1)
+      endif
+   endif
+endif ! of if (count(indx)==0)
+
+
+! 2.1. Build dx* (data on lon-lat grid, with redundant longitude)
+
+if (dim.eq.3) then
+  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
+  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
+  ENDIF
+else ! dim.eq.2
+  dx2(1:nbp_lon,:)=px2_glo(:,:)
+  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
+    dx2(nbp_lon+1,:)=dx2(1,:)
+  ENDIF
+endif
+
+
+! 2.2. Add current values to previously stored sums
+
+if (dim.eq.3) then
+
+   mean3d(:,:,:)=mean3d(:,:,:)+dx3(:,:,:)
+   sd3d(:,:,:)=sd3d(:,:,:)+dx3(:,:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd3d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean3d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "wstats error writing :",trim(nom)
+     write (*,*) NF_STRERROR(ierr)
+     call abort_physic("wstats","Failed writing "//trim(nom),1)
+  endif
+
+else if (dim.eq.2) then
+
+   mean2d(:,:)= mean2d(:,:)+dx2(:,:)
+   sd2d(:,:)=sd2d(:,:)+dx2(:,:)**2
+
+#ifdef NC_DOUBLE
+   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
+   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,sizes,sd2d)
+#else
+   ierr = NF_PUT_VARA_REAL(nid,meanid,start,sizes,mean2d)
+   ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "wstats error writing :",trim(nom)
+     write(*,*) "start:",start
+     write(*,*) "sizes:",sizes
+     write(*,*) "mean2d:",mean2d
+     write(*,*) "sd2d:",sd2d
+     write (*,*) NF_STRERROR(ierr)
+     call abort_physic("wstats","Failed writing "//trim(nom),1)
+  endif
+
+endif ! of if (dim.eq.3) elseif (dim.eq.2)
+
+  ierr= NF_CLOSE(nid)
+endif ! of if (is_master)
+
+end subroutine wstats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine inistats(ierr)
+
+! routine to initialize the stats file (i.e. create file, write
+! time independent coordinates, etc.)
+
+use mod_phys_lmdz_para, only : is_master
+USE vertical_layers_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
+USE nrtype, ONLY: pi
+USE time_phylmdz_mod, ONLY: daysec,dtphys
+USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(out) :: ierr
+
+integer :: nid
+integer :: l,nsteppd
+real, dimension(nbp_lev) ::  sig_s
+real,allocatable :: lon_reg_ext(:) ! extended longitudes
+integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
+real, dimension(istime) :: lt
+integer :: nvarid
+
+
+IF (nbp_lon*nbp_lat==1) THEN
+  ! 1D model
+  ALLOCATE(lon_reg_ext(1))
+ELSE
+  ! 3D model
+  ALLOCATE(lon_reg_ext(nbp_lon+1))
+ENDIF
+
+write (*,*) 
+write (*,*) '                        || STATS ||'
+write (*,*) 
+write (*,*) 'daysec',daysec
+write (*,*) 'dtphys',dtphys
+nsteppd=nint(daysec/dtphys)
+write (*,*) 'nsteppd=',nsteppd
+
+if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) then
+  call abort_physic("inistats",'1 sol .ne. n physics steps',1)
+endif
+
+if(mod(nsteppd,istime).ne.0) then
+  call abort_physic("inistats",'1 sol .ne. n*istime physics steps',1)
+endif
+
+istats=nsteppd/istime
+write (*,*) 'istats=',istats
+write (*,*) 'Storing ',istime,'times per day'
+write (*,*) 'thus every ',istats,'physical timestep '
+write (*,*) 
+
+do l= 1, nbp_lev
+  sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
+  pseudoalt(l)=-10.*log(presnivs(l)/preff)   
+enddo
+
+lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+IF (nbp_lon*nbp_lat/=1) THEN
+  ! In 3D, add extra redundant point (180 degrees,
+  ! since lon_reg starts at -180)
+  lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+ENDIF
+
+if (is_master) then
+  ! only the master needs do this
+  ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
+  if (ierr.ne.NF_NOERR) then
+    write (*,*) NF_STRERROR(ierr)
+    call abort_physic("inistats",'failed creating stats.nc',1)
+  endif
+
+  ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
+  IF (nbp_lon*nbp_lat==1) THEN
+    ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
+  ELSE
+    ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+  ENDIF
+  ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
+  ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
+  ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
+
+  ierr = NF_ENDDEF(nid)
+
+  call def_var_stats(nid,"Time","Time", &
+                     "days since 0000-00-0 00:00:00",1, &
+                     [idim_time],nvarid,ierr)
+  ! Time is initialised later by mkstats subroutine
+
+  call def_var_stats(nid,"latitude","latitude", &
+                     "degrees_north",1,[idim_lat],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+
+  call def_var_stats(nid,"longitude","East longitude", &
+                     "degrees_east",1,[idim_lon],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
+#endif
+
+  ! Niveaux verticaux, aps et bps
+  ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+  ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid)
+#else
+  ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid)
+#endif
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
+  ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
+  ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
+#endif 
+  call def_var_stats(nid,"aps","hybrid pressure at midlayers", &
+                     " ",1,[idim_llm],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
+#else
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
+#endif
+
+  call def_var_stats(nid,"bps","hybrid sigma at midlayers", &
+                     " ",1,[idim_llm],nvarid,ierr)
+#ifdef NC_DOUBLE
+  ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
+#else 
+  ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
+#endif
+
+  ierr=NF_CLOSE(nid)
+
+endif ! of if (is_master)
+
+end subroutine inistats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+
+! routine to initialize writing a variable in the stats file
+
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
+
+implicit none
+
+include "netcdf.inc"
+
+integer, intent(in) :: nid,varid,dim,indx,ngrid
+real, dimension(ngrid,nbp_lev), intent(in) :: px
+integer, intent(out) :: ierr
+
+integer :: l,i,j,ig0
+integer, dimension(4) :: start,sizes
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
+real, dimension(nbp_lon+1,nbp_lat) :: dx2
+real :: dx3_1d(nbp_lev) ! for 1D outputs
+real :: dx2_1d ! for 1D outputs
+
+if (dim.eq.3) then
+
+   start=(/1,1,1,indx/)
+   if (klon_glo>1) then ! general 3D case
+     sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+   else
+     sizes=(/1,1,nbp_lev,1/)
+   endif
+
+!  Passage variable physique -->  variable dynamique
+
+   if (klon_glo>1) then ! general case
+    DO l=1,nbp_lev
+      DO i=1,nbp_lon+1
+         dx3(i,1,l)=px(1,l)
+         dx3(i,nbp_lat,l)=px(ngrid,l)
+      ENDDO
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
+            dx3(i,j,l)=px(ig0+i,l)
+         ENDDO
+         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
+      ENDDO
+    ENDDO
+   else ! 1D model case
+     dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
+   endif
+
+#ifdef NC_DOUBLE
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
+   else
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3_1d)
+   endif
+#else
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
+   else
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3_1d)
+   endif
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "inivar error writing variable"
+     write (*,*) NF_STRERROR(ierr)
+     call abort_physic("inivar","error writing variable",1)
+  endif
+
+else if (dim.eq.2) then
+
+      start=(/1,1,indx,0/)
+      if (klon_glo>1) then ! general 3D case
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
+
+!    Passage variable physique -->  physique dynamique
+
+  if (klon_glo>1) then ! general case
+  DO i=1,nbp_lon+1
+     dx2(i,1)=px(1,1)
+     dx2(i,nbp_lat)=px(ngrid,1)
+  ENDDO
+  DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
+        dx2(i,j)=px(ig0+i,1)
+     ENDDO
+     dx2(nbp_lon+1,j)=dx2(1,j)
+  ENDDO
+  else ! 1D model case
+    dx2_1d=px(1,1)
+  endif
+  
+#ifdef NC_DOUBLE
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
+   else
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2_1d)
+   endif
+#else
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
+   else
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2_1d)
+   endif
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "inivar error writing variable"
+     write (*,*) NF_STRERROR(ierr)
+     call abort_physic("inivar","error writing variable",1)
+  endif
+
+endif
+
+end subroutine inivar
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
+
+! This subroutine defines variable 'name' in a (pre-existing and opened)
+! NetCDF file (known from its NetCDF ID 'nid').
+! The number of dimensions 'nbdim' of the variable, as well as the IDs of
+! corresponding dimensions must be set (in array 'dimids').
+! Upon successfull definition of the variable, 'nvarid' contains the
+! NetCDF ID of the variable.
+! The variables' attributes 'title' (Note that 'long_name' would be more
+! appropriate) and 'units' are also set. 
+
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(in) :: nid ! NetCDF file ID
+character(len=*),intent(in) :: name ! the variable's name
+character(len=*),intent(in) :: title ! 'title' attribute of variable
+character(len=*),intent(in) :: units ! 'units' attribute of variable
+integer,intent(in) :: nbdim ! number of dimensions of the variable
+integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
+                                              ! the variable is defined along
+integer,intent(out) :: nvarid ! NetCDF ID of the variable
+integer,intent(out) :: ierr ! returned NetCDF staus code
+
+! 1. Switch to NetCDF define mode 
+ierr=NF_REDEF(nid)
+
+! 2. Define the variable
+#ifdef NC_DOUBLE
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
+#else
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
+#endif
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   call abort_physic("def_var_stats","Failed defining "//trim(name),1)
+endif
+
+! 3. Write attributes
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
+                     len_trim(adjustl(title)),adjustl(title))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
+endif
+
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
+                     len_trim(adjustl(units)),adjustl(units))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
+endif
+
+! 4. Switch out of NetCDF define mode
+ierr = NF_ENDDEF(nid)
+
+end subroutine def_var_stats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine mkstats(ierr)
+
+!  This routine finalizes tha stats.nc file from sums and sums of squares
+!  to means and standard deviations. The data file is overwritten in place.
+
+use mod_phys_lmdz_para, only : is_master
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
+
+implicit none
+
+include "netcdf.inc"
+
+integer,intent(out) :: ierr ! netCDF return error code
+integer :: nid,nbvar,i,ndims,lt,nvarid
+integer, dimension(4) :: id,varid,start,size
+integer, dimension(5) :: dimids
+character (len=50) :: name,nameout,units,title
+real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
+real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
+real, dimension(istime) :: time
+real, dimension(nbp_lat) :: lat
+real,allocatable :: lon(:)
+real, dimension(nbp_lev) :: alt
+logical :: lcopy=.true.
+!integer :: latid,lonid,altid,timeid
+integer :: meanid,sdid
+!integer, dimension(4) :: dimout
+
+if (is_master) then
+! only the master needs do all this
+
+! Incrementation of count for the last step, which is not done in wstats
+count(istime)=count(istime)+1
+
+if (klon_glo>1) then
+  allocate(lon(nbp_lon+1))
+  allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sum2d(nbp_lon+1,nbp_lat))
+  allocate(square2d(nbp_lon+1,nbp_lat))
+  allocate(mean2d(nbp_lon+1,nbp_lat))
+  allocate(sd2d(nbp_lon+1,nbp_lat))
+else ! 1D model case
+  allocate(lon(1))
+  allocate(sum3d(1,1,nbp_lev))
+  allocate(square3d(1,1,nbp_lev))
+  allocate(mean3d(1,1,nbp_lev))
+  allocate(sd3d(1,1,nbp_lev))
+  allocate(sum2d(1,1))
+  allocate(square2d(1,1))
+  allocate(mean2d(1,1))
+  allocate(sd2d(1,1))
+endif
+
+ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
+
+! We catch the id of dimensions of the stats file
+
+ierr= NF_INQ_DIMID(nid,"latitude",id(1))
+ierr= NF_INQ_DIMID(nid,"longitude",id(2))
+ierr= NF_INQ_DIMID(nid,"altitude",id(3))
+ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+ierr= NF_INQ_VARID(nid,"latitude",varid(1))
+ierr= NF_INQ_VARID(nid,"longitude",varid(2))
+ierr= NF_INQ_VARID(nid,"altitude",varid(3))
+ierr= NF_INQ_VARID(nid,"Time",varid(4))
+
+! Time initialisation
+
+do i=1,istime
+   time(i)=i*24./istime
+#ifdef NC_DOUBLE
+   ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),i,1,time(i))
+#else
+   ierr= NF_PUT_VARA_REAL(nid,varid(4),i,1,time(i))
+#endif
+enddo
+
+! We catche the values of the variables
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon)
+         ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,varid(1),lat)
+         ierr = NF_GET_VAR_REAL(nid,varid(2),lon)
+         ierr = NF_GET_VAR_REAL(nid,varid(3),alt)
+#endif
+
+! We catch the number of variables in the stats file
+ierr = NF_INQ_NVARS(nid,nbvar)
+
+! to catche the "real" number of variables (without the "additionnal variables")
+nbvar=(nbvar-4)/2 
+
+do i=1,nbvar
+   varid=(i-1)*2+5
+
+   ! What's the variable's name?
+   ierr=NF_INQ_VARNAME(nid,varid,name)
+   write(*,*) "OK variable ",name
+   ! Its units?
+   units=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"units",units)
+   ! Its title?
+   title=" "
+   ierr=NF_GET_ATT_TEXT(nid,varid,"title",title)
+   ! Its number of dimensions?   
+   ierr=NF_INQ_VARNDIMS(nid,varid,ndims)
+   ! Its values?
+
+   if(ndims==4) then ! lat, lon, alt & time
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=altid
+!      dimout(4)=timeid
+
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else ! 1D model
+        size=(/1,1,nbp_lev,1/)
+      endif
+      do lt=1,istime
+         start=(/1,1,1,lt/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square3d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum3d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square3d)
+#endif
+         ! Calculation of these variables
+         mean3d=sum3d/count(lt)
+         sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd3d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean3d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd3d)
+#endif
+      enddo
+
+    else if (ndims.eq.3) then
+
+!      dimout(1)=lonid
+!      dimout(2)=latid
+!      dimout(3)=timeid
+
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        size=(/1,1,1,0/)
+      endif
+      do lt=1,istime
+         start=(/1,1,lt,0/)
+         ! Extraction of the "source" variables
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square2d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum2d)
+         ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square2d)
+#endif
+         ! Calculation of these variables
+         mean2d=sum2d/count(lt)
+         sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2))
+         ! Writing of the variables
+#ifdef NC_DOUBLE
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd2d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean2d)
+         ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd2d)
+#endif
+      enddo
+
+    endif 
+enddo
+
+ierr= NF_CLOSE(nid)
+
+endif ! of if (is_master)
+
+end subroutine mkstats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+end module wstats_mod
