#include "options.inc"
c
c-----------------------------------------------------------------------
c     diagnose meridional overturning 
c     and write to either NetCDF or binary file
c-----------------------------------------------------------------------
c
      subroutine init_over_cdf
      use spflame_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer ncid,iret,i,j,k
      real :: spval=-9.9e12
      integer lat_udim,itimedim, depth_tdim,depth_tid
      integer lat_uid,itimeid,vsfid, dims(4), corner(4), edges(4)
      character name*24, unit*16, text*80
#ifdef netcdf_real4
      real (kind=4) , allocatable :: v2(:)
#else
      real          , allocatable :: v2(:)
#endif
      if (my_pe == 0) then
      print*,' initializing NetCDF output file ',over_file
      ncid = nccre (over_file, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call store_info_cdf(ncid)
c     dimensions
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
      dims=(/lat_udim,depth_tdim,itimedim,1/)
      vsfid=ncvdef (ncid,'vsf', NCFLOAT,3,dims,iret)
c     attributes of the grid
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c     attributes of variables
      name = 'Meridional overturning'; unit = 'Sv'
      call dvcdf(ncid,vsfid,name,24,unit,16,spval)
      call ncendf(ncid, iret)
      allocate( v2(max(km,jmt)) )
      corner(1) = 1; edges(1) = jmt
      v2(1:jmt)=yu
      call ncvpt(ncid, Lat_uid, corner, edges,v2, iret)
      corner(1) = 1; edges(1) = km
      v2(1:km)=zt/100.
      call ncvpt(ncid, depth_tid, corner, edges,v2, iret)
      deallocate(v2)
      call ncclos (ncid, iret)
      print*,' done'; call sub_flush(6)
      endif ! my_pe==0
#else
c     needs no intialisaztion
#endif
      end subroutine init_over_cdf


      subroutine diag_over
      use spflame_module
      implicit none
      integer i,j,k, is,ie,js,je
      real vbarx(km),vsf(js_pe:je_pe,km)
#ifdef netcdf_diagnostics
#include "netcdf.inc"

#ifdef netcdf_real4
      real (kind=4) :: var(js_pe:je_pe,km)
#else
      real          :: var(js_pe:je_pe,km)
#endif

#else
      integer :: io
      character (len=80) :: iotext
      character (len=60) :: expnam
      real :: vsf_big(jmt,km), reltim = 0.
#endif
      integer ncid,iret,npe
      real :: spval=-9.9e12, tt
      integer vsfid,timeid,timedim,len
      integer corner(4), edges(4)

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      do j=js_pe,je_pe
       do k=1,km
        vbarx=0.
        do i=is, ie
         vbarx(k) = vbarx(k) + u(i,k,j,2,tau)*csu(j)*dxu(i)*umask(i,k,j)
#ifdef partial_cell
     &                      *dhu(i,k,j)
#else
     &                      *dzt(k)
#endif
        enddo
        if (k == 1) then
         vsf(j,k) = vbarx(k)
	else
         vsf(j,k) = vsf(j,k-1) + vbarx(k)
	endif
       enddo
       call sum_along_jrow(vsf(j,:),km)
#ifdef netcdf_diagnostics
       do k=1,km
        npe=0
        do i=1,imt
         if (kmt_big(i,j)>=k) npe=1
        enddo
        vsf(j,k)=vsf(j,k)/1.e12*npe+spval*(1-npe)
       enddo
#endif
      enddo

      if (my_pe==0)
     &  print*,' --> writing overturning to file ',
     &        over_file(1:len_trim(over_file))

#ifdef netcdf_diagnostics
c     write to netcdf format here
      do npe=0,n_pes,n_pes_i
       call barrier
       if (my_pe==npe) then
        iret=nf_open(over_file,NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_varid(ncid,'vsf',vsfid)
        iret=nf_inq_varid(ncid,'Time',timeid)
        iret=nf_inq_dimid(ncid,'Time',timedim)
        iret=nf_inq_dimlen(ncid, timedim,len)
        if (my_pe==0) then
         len=len+1
         corner(1)=len
         edges(1)=1
         call read_stamp(current_stamp,tt)
         print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',tt,')',
     &          ' (time steps in file : ',len,')'
         var(js_pe,1)=tt
         iret= nf_put_vara_real (ncid,timeid, corner,edges,var(js_pe,1))
        endif
        corner=(/js_pe,1,len,1/); edges=(/je_pe-js_pe+1,km,1,1/)
        var=vsf
        iret= nf_put_vara_real (ncid,vsfid ,corner, edges,var)
        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
#else
c     write old flame format here
      if ( my_pe == 0) then
       call getunit (io, over_file,'u s a ieee')
       iotext = 'no iotext'
       expnam = 'SPFLAME, some experiment'
       write (io) current_stamp, iotext, expnam
c           reltim = prelyr
       write (io) jmt, km, reltim
       write (io) current_stamp, iotext, expnam
       write (io) zw
       write (io) current_stamp, iotext, expnam
       write (io) yu
c          vsf_big(1:je_pe,:)=vsf(1:je_pe,:)
       write (io) current_stamp, iotext, expnam
      endif
      vsf_big(js_pe:je_pe,:)=vsf
c      call pe0_recv_vec_along_jrow(vsf_big,km)
      call pe0_recv_merid_vec(vsf_big,km,1)
      if ( my_pe == 0) then
       write(io) vsf_big
       close(io)
      endif
#endif
      end subroutine diag_over

