c     Dump NetCDF RI statistics for SOHC data
      implicit none
      integer k,kk,kkk,nx,ny,i,j,nn,obs,obs3
      parameter (nx= 641,ny=nyyyy,nn=3)
      integer ncid,rcode,count(nn),start(nn),id
      real ohc(nx,ny,1),ssha(nx,ny,1),omld(nx,ny,1)
      real sst(nx,ny,1),missing,min,max
      double precision mean,std

      include '/data/starfs1/libs/netcdf-3.6.3/include/netcdf.inc'
c    ***************************************************************
c      open(8, file='data/sohc_all_dump.dat',form='formatted')
      open(9, file='data/sohc_all_statistics.dat',form='formatted')
c    ***************************************************************
      missing=-999
c    ***************************************************************

      start(1)=1
      start(2)=1
      start(3)=1
      count(1)=nx
      count(2)=ny
      count(3)=1

      ncid=ncopn('sohc_all.nc',ncnowrit,rcode)

      id=ncvid(ncid,'ohc',rcode)
      call ncvgt(ncid,id,start,count,ohc,rcode)
      write(91,'(10f10.2)') ohc

      id=ncvid(ncid,'sst',rcode)
      call ncvgt(ncid,id,start,count,sst,rcode)
      write(92,'(10f10.2)') sst

      id=ncvid(ncid,'ssha',rcode)
      call ncvgt(ncid,id,start,count,ssha,rcode)
      write(93,'(10f10.2)') ssha

      id=ncvid(ncid,'omld',rcode)
      call ncvgt(ncid,id,start,count,omld,rcode)
      write(94,'(10f10.2)') omld


c    ***************************************************************
      call sub_stat(nx,ny,ohc,obs,obs3,std,mean,min,max,missing)
      write(9,101) float(obs),float(obs3),std,mean,min,max

      call sub_stat(nx,ny,sst,obs,obs3,std,mean,min,max,missing)
      write(9,101) float(obs),float(obs3),std,mean,min,max

      call sub_stat(nx,ny,ssha,obs,obs3,std,mean,min,max,missing)
      write(9,101) float(obs),float(obs3),std,mean,min,max

      call sub_stat(nx,ny,omld,obs,obs3,std,mean,min,max,missing)
      write(9,101) float(obs),float(obs3),std,mean,min,max
c    ***************************************************************

 101  format(6f14.5)

      stop 
      end