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

implicit none

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

integer,intent(in) :: ngrid
character (len=*) :: nom,titre,unite
integer,intent(in) :: dim
integer,parameter :: iip1=iim+1
integer,parameter :: jjp1=jjm+1
real, dimension(ngrid,llm) :: px
real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
character (len=50) :: namebis
character (len=50), save :: firstvar
integer :: ierr,varid,nbdim,nid
integer :: meanid,sdid
integer, dimension(4)  :: id,start,size
logical, save :: firstcall=.TRUE.
integer :: l,i,j,ig0
integer,save :: index

integer, save :: step=0


if (firstcall) then
   firstcall=.false.
   firstvar=trim((nom))
   call inistats(ierr)
endif

if (firstvar==nom) then ! If we're back to the first variable
      step = step + 1
endif

if (mod(step,istats).ne.0) then
   RETURN
endif

ierr = NF_OPEN("stats.nc",NF_WRITE,nid)

namebis=trim(nom)
ierr= NF_INQ_VARID(nid,namebis,meanid)

if (ierr.ne.NF_NOERR) then

   if (firstvar==nom) then 
      index=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)
   call inivar(nid,meanid,ngrid,dim,index,px,ierr)
   namebis=trim(nom)//"_sd"
   call def_var(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr)
   call inivar(nid,sdid,ngrid,dim,index,px,ierr)

   ierr= NF_CLOSE(nid)
   return

else
   namebis=trim(nom)//"_sd"
   ierr= NF_INQ_VARID(nid,namebis,sdid)

endif

if (firstvar==nom) then 
   count(index)=count(int(index))+1
   index=index+1
   if (index>istime) then
      index=1
   endif
endif

if (count(index)==0) then
   if (dim.eq.3) then
      start=(/1,1,1,index/)
      size=(/iip1,jjp1,llm,1/)
      mean3d=0
      sd3d=0
   else if (dim.eq.2) then
      start=(/1,1,index,0/)
      size=(/iip1,jjp1,1,0/)
      mean2d=0
      sd2d=0
   endif
else
   if (dim.eq.3) then
      start=(/1,1,1,index/)
      size=(/iip1,jjp1,llm,1/)
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d)
      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d)
#else
      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d)
      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d)
#endif
      if (ierr.ne.NF_NOERR) then
         write (*,*) NF_STRERROR(ierr)
         stop ""
      endif

   else if (dim.eq.2) then
      start=(/1,1,index,0/)
      size=(/iip1,jjp1,1,0/)
#ifdef NC_DOUBLE
      ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d)
      ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d)
#else
      ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d)
      ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d)
#endif
      if (ierr.ne.NF_NOERR) then
         write (*,*) NF_STRERROR(ierr)
         stop ""
      endif
   endif
endif

if (dim.eq.3) then

!  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

   mean3d= mean3d+dx3
   sd3d= sd3d+dx3**2

#ifdef NC_DOUBLE
   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d)
   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d)
#else
   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d)
   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d)
#endif

else if (dim.eq.2) then

!    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

   mean2d= mean2d+dx2
   sd2d= sd2d+dx2**2

#ifdef NC_DOUBLE
   ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d)
   ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d)
#else
   ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d)
   ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d)
#endif

endif

ierr= NF_CLOSE(nid)

end

!======================================================
subroutine inivar(nid,varid,ngrid,dim,index,px,ierr)

implicit none

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

integer, intent(in) :: nid,varid,dim,index,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,size
real, dimension(iip1,jjp1,llm) :: dx3
real, dimension(iip1,jjp1) :: dx2

if (dim.eq.3) then

   start=(/1,1,1,index/)
   size=(/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,size,dx3)
#else
   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3)
#endif

else if (dim.eq.2) then

      start=(/1,1,index,0/)
      size=(/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,size,dx2)
#else
   ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2)
#endif

endif

end