Changeset 2958 for trunk/LMDZ.GENERIC

Timestamp:

May 5, 2023, 2:23:32 PM (20 months ago)

Author:

emillour

Message:

Generic PCM:
Upgrade wstats following the Mars PCM one. It is now a module and there no
longer needs to have if (callstats) around a call to wstats as it managed
internally in the wstats routine.
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).
EM

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/aeronostd/calchim_asis.F90 (modified) (3 diffs)
• libf/phystd/callkeys_mod.F90 (modified) (1 diff)
• libf/phystd/inifis_mod.F90 (modified) (1 diff)
• libf/phystd/inistats.F (deleted)
• libf/phystd/mkstat.F90 (deleted)
• libf/phystd/physiq_mod.F90 (modified) (4 diffs)
• libf/phystd/statto_mod.F90 (deleted)
• libf/phystd/wstats_mod.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/phystd/wstats.F90) (17 diffs)

trunk/LMDZ.GENERIC/README

@@ -1801 +1801 @@
 Also added some "intent()" in optci arguments and increased length of string
 to store varaible name in writediagfi.
+
+== 05/05/2023 == EM
+Upgrade wstats following the Mars PCM one. It is now a module and there no
+longer needs to have if (callstats) around a call to wstats as it managed
+internally in the wstats routine.
+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).

trunk/LMDZ.GENERIC/libf/aeronostd/calchim_asis.F90

@@ -14 +14 @@
       use chimiedata_h
       use radcommon_h,      only: glat
+      use wstats_mod, only: wstats
 
       implicit none
@@ -339 +340 @@
       if (photochem .and. output) then
 
-         if (callstats) then
             ! photoloysis
             do i=1,nb_phot_hv_max
@@ -371 +371 @@
                end do
             endif ! end depos
-         endif ! end callstats
 
          ! photoloysis

trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90

@@ -8 +8 @@
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
-      logical,save :: callstats ! .true. to output stats.nc files
-!$OMP THREADPRIVATE(callstats)
       logical,save :: callgasvis,continuum,graybody
 !$OMP THREADPRIVATE(callgasvis,continuum,graybody)

trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

@@ -23 +23 @@
   use planetwide_mod, only: planetwide_sumval
   use callkeys_mod
+  use wstats_mod, only: callstats
   use ioipsl_getin_p_mod, only : getin_p
   use mod_phys_lmdz_para, only : is_parallel, is_master, bcast

trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

@@ -35 +35 @@
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyetat0_mod, only: phyetat0
+      use wstats_mod, only: callstats, wstats, mkstats
       use phyredem, only: physdem0, physdem1
       use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
@@ -50 +51 @@
       use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
                               calllott_nonoro, callrad, callsoil, nosurf, &
-                              callstats, &
                               calltherm, CLFvarying, co2cond, corrk, diagdtau, &
                               diurnal, enertest, fat1au, flatten, j2, &
@@ -2259 +2259 @@
 
 
-!-----------------------------------
-!        Saving statistics :
-!-----------------------------------
-
-!    Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
-!           which can later be used to make the statistic files of the run:
-!           "stats")          only possible in 3D runs !!!
-
-         
-      if (callstats) then
+! -----------------------------------------------------------------
+! WSTATS: Saving statistics
+! -----------------------------------------------------------------
+! ("stats" stores and accumulates key variables in file "stats.nc"
+!  which can later be used to make the statistic files of the run:
+!  if flag "callstats" from callphys.def is .true.)
+         
 
          call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
@@ -2343 +2340 @@
          endif
 
-         if(lastcall) then
+         if(lastcall.and.callstats) then
             write (*,*) "Writing stats..."
             call mkstats(ierr)
          endif
          
-      endif ! end of 'callstats'
 
 #ifndef MESOSCALE

trunk/LMDZ.GENERIC/libf/phystd/wstats_mod.F90

@@ +1 @@
+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(count)
+
+contains
+
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
 
-use statto_mod, only: istats,istime,count
+! 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
+                              nbp_lon, nbp_lat, nbp_lev, &
+                              grid_type, unstructured
 implicit none
 
@@ -17 +42 @@
 character (len=50) :: namebis
 character (len=50), save :: firstvar
-!$OMP THREADPRIVATE(firstvar)
+!$OMP THREADPRIVATE(mean3d,sd3d,dx3,mean2d,sd2d,dx2,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)
+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
+
+!$OMP THREADPRIVATE(stats_def,n_name_def,name_def)
 
 ! Added to work in parallel mode
@@ -44 +78 @@
 #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)
@@ -64 +137 @@
      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
@@ -74 +147 @@
    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
@@ -171 +261 @@
 
    write (*,*) "====================="
-   write (*,*) "STATS: creation de ",nom
+   write (*,*) "STATS: creating ",nom
    namebis=trim(nom)
    call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr)
@@ -245 +335 @@
          write (*,*) "wstats error reading :",trim(nom)
          write (*,*) NF_STRERROR(ierr)
-         stop ""
+         call abort_physic("wstats","Failed reading "//trim(nom),1)
       endif
 
@@ -265 +355 @@
          write (*,*) "wstats error reading :",trim(nom)
          write (*,*) NF_STRERROR(ierr)
-         stop ""
+         call abort_physic("wstats","Failed reading "//trim(nom),1)
       endif
    endif
@@ -303 +393 @@
      write (*,*) "wstats error writing :",trim(nom)
      write (*,*) NF_STRERROR(ierr)
-     stop ""
+     call abort_physic("wstats","Failed writing "//trim(nom),1)
   endif
 
@@ -325 +415 @@
      write(*,*) "sd2d:",sd2d
      write (*,*) NF_STRERROR(ierr)
-     stop ""
+     call abort_physic("wstats","Failed writing "//trim(nom),1)
   endif
 
@@ -335 +425 @@
 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
 
@@ -399 +642 @@
      write (*,*) "inivar error writing variable"
      write (*,*) NF_STRERROR(ierr)
-     stop ""
+     call abort_physic("inivar","error writing variable",1)
   endif
 
@@ -445 +688 @@
      write (*,*) "inivar error writing variable"
      write (*,*) NF_STRERROR(ierr)
-     stop ""
+     call abort_physic("inivar","error writing variable",1)
   endif
 
@@ -491 +734 @@
    write(*,*) "def_var_stats: Failed defining variable "//trim(name)
    write(*,*) NF_STRERROR(ierr)
-   stop ""
+   call abort_physic("def_var_stats","Failed defining "//trim(name),1)
 endif
 
@@ -500 +743 @@
    write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
    write(*,*) NF_STRERROR(ierr)
-   stop ""
+   call abort_physic("def_var_stats","Failed writing title for "//trim(name),1)
 endif
 
@@ -508 +751 @@
    write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
    write(*,*) NF_STRERROR(ierr)
-   stop ""
+   call abort_physic("def_var_stats","Failed writing units for "//trim(name),1)
 endif
 
@@ -516 +759 @@
 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
+   nvarid=(i-1)*2+5
+
+   ! What's the variable's name?
+   ierr=NF_INQ_VARNAME(nid,nvarid,name)
+   write(*,*) "OK variable ",name
+   ! Its units?
+   units=" "
+   ierr=NF_GET_ATT_TEXT(nid,nvarid,"units",units)
+   ! Its title?
+   title=" "
+   ierr=NF_GET_ATT_TEXT(nid,nvarid,"title",title)
+   ! Its number of dimensions?   
+   ierr=NF_INQ_VARNDIMS(nid,nvarid,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,nvarid,start,size,sum3d)
+         ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square3d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum3d)
+         ierr = NF_GET_VARA_REAL(nid,nvarid+1,start,size,square3d)
+#endif
+         ! Calculation of these nvariables
+         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,nvarid,start,size,mean3d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd3d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean3d)
+         ierr = NF_PUT_VARA_REAL(nid,nvarid+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,nvarid,start,size,sum2d)
+         ierr = NF_GET_VARA_DOUBLE(nid,nvarid+1,start,size,square2d)
+#else
+         ierr = NF_GET_VARA_REAL(nid,nvarid,start,size,sum2d)
+         ierr = NF_GET_VARA_REAL(nid,nvarid+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,nvarid,start,size,mean2d)
+         ierr = NF_PUT_VARA_DOUBLE(nid,nvarid+1,start,size,sd2d)
+#else
+         ierr = NF_PUT_VARA_REAL(nid,nvarid,start,size,mean2d)
+         ierr = NF_PUT_VARA_REAL(nid,nvarid+1,start,size,sd2d)
+#endif
+      enddo
+
+    endif 
+enddo
+
+ierr= NF_CLOSE(nid)
+
+endif ! of if (is_master)
+
+end subroutine mkstats
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+end module wstats_mod