subroutine wstats(ngrid,nom,titre,unite,dim,px)

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
use statto_mod, only: istats,istime,count

implicit none

#include "dimensions.h"
!#include "dimphys.h"
#include "netcdf.inc"

integer,intent(in) :: ngrid
character (len=*),intent(in) :: nom,titre,unite
integer,intent(in) :: dim
integer,parameter :: iip1=iim+1
integer,parameter :: jjp1=jjm+1
real,intent(in) :: px(ngrid,llm)
real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
real, dimension(iip1,jjp1) :: 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,llm) ! to store a 3D data set on global physics grid
real px3_glo(iim,jjp1,llm) ! 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(iim,jjp1) ! to store a 2D data set on global lonxlat grid
real px2(ngrid)
real px3(ngrid,llm)
#else
! When not running in parallel mode:
real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
real px3_glo(iim,jjp1,llm) ! 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(iim,jjp1) ! 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)
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:llm)=px(1:ngrid,1:llm)
  ! 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:iim,:,:)=px3_glo(1:iim,:,:)
    dx3(iip1,:,:)=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:iim,:)=px2_glo(1:iim,:)
            dx2(iip1,:)=dx2(1,:)
          endif
!$OMP END MASTER
!$OMP BARRIER
 endif
#else
  if (dim.eq.3) then
    px3_glop(:,1:llm)=px(:,1:llm)
!  Passage variable physique -->  variable dynamique
    DO l=1,llm
      DO i=1,iim
         px3_glo(i,1,l)=px(1,l)
         px3_glo(i,jjp1,l)=px(ngrid,l)
      ENDDO
      DO j=2,jjm
         ig0= 1+(j-2)*iim
         DO i=1,iim
            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,iim
     px2_glo(i,1)=px(1,1)
     px2_glo(i,jjp1)=px(ngrid,1)
   ENDDO
   DO j=2,jjm
     ig0= 1+(j-2)*iim
     DO i=1,iim
        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: creation de ",nom
   namebis=trim(nom)
   call def_var(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(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/)
      sizes=(/iip1,jjp1,llm,1/)
      mean3d(:,:,:)=0
      sd3d(:,:,:)=0
   else if (dim.eq.2) then
      start=(/1,1,indx,0/)
      sizes=(/iip1,jjp1,1,0/)
      mean2d(:,:)=0
      sd2d(:,:)=0
   endif
else
   ! load values from file
   if (dim.eq.3) then
      start=(/1,1,1,indx/)
      sizes=(/iip1,jjp1,llm,1/)
#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 (*,*) NF_STRERROR(ierr)
         stop ""
      endif

   else if (dim.eq.2) then
      start=(/1,1,indx,0/)
      sizes=(/iip1,jjp1,1,0/)
#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 (*,*) NF_STRERROR(ierr)
         stop ""
      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:iim,:,:)=px3_glo(:,:,:)
  dx3(iip1,:,:)=dx3(1,:,:)
else ! dim.eq.2
  dx2(1:iim,:)=px2_glo(:,:)
  dx2(iip1,:)=dx2(1,:)
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

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

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)

implicit none

include "dimensions.h"
!include "dimphys.h"
include "netcdf.inc"

integer, intent(in) :: nid,varid,dim,indx,ngrid
real, dimension(ngrid,llm), intent(in) :: px
integer, intent(out) :: ierr

integer,parameter :: iip1=iim+1
integer,parameter :: jjp1=jjm+1

integer :: l,i,j,ig0
integer, dimension(4) :: start,sizes
real, dimension(iip1,jjp1,llm) :: dx3
real, dimension(iip1,jjp1) :: dx2

if (dim.eq.3) then

   start=(/1,1,1,indx/)
   sizes=(/iip1,jjp1,llm,1/)

!  Passage variable physique -->  variable dynamique

   DO l=1,llm
      DO i=1,iip1
         dx3(i,1,l)=px(1,l)
         dx3(i,jjp1,l)=px(ngrid,l)
      ENDDO
      DO j=2,jjm
         ig0= 1+(j-2)*iim
         DO i=1,iim
            dx3(i,j,l)=px(ig0+i,l)
         ENDDO
         dx3(iip1,j,l)=dx3(1,j,l)
      ENDDO
   ENDDO

#ifdef NC_DOUBLE
   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
#else
   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
#endif

else if (dim.eq.2) then

      start=(/1,1,indx,0/)
      sizes=(/iip1,jjp1,1,0/)

!    Passage variable physique -->  physique dynamique

  DO i=1,iip1
     dx2(i,1)=px(1,1)
     dx2(i,jjp1)=px(ngrid,1)
  ENDDO
  DO j=2,jjm
     ig0= 1+(j-2)*iim
     DO i=1,iim
        dx2(i,j)=px(ig0+i,1)
     ENDDO
     dx2(iip1,j)=dx2(1,j)
  ENDDO

#ifdef NC_DOUBLE
   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
#else
   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
#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)
   stop ""
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)
   stop ""
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)
   stop ""
endif

! 4. Switch out of NetCDF define mode
ierr = NF_ENDDEF(nid)

end subroutine def_var_stats