#include "options.inc"



      subroutine diagnose
c=======================================================================
c      Diagnostic routines
c=======================================================================
      use cpflame_module
      implicit none
      integer :: i,k,j,js,je
      real :: totke,tbar
      real :: area,fx,cflcrt,cflu,cflv,cflw,cl
      character (len=32) :: stamp
c--------------------------------------------------------------
c   calculate basin integrated kinetic energy and thickness
c--------------------------------------------------------------
      js=max(2,js_pe); je = min(je_pe,jmt-1)

      if (snapshot_time_step.or.initial_time==current_time) then

      fx = dx**2*dz
      totke=0.0;area=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          area=area + fx*maskT(i,j,k)
          totke = totke + u(i,j,k,1,tau)**2/2.0*fx*masku(i,j,k)
          totke = totke + u(i,j,k,2,tau)**2/2.0*fx*maskv(i,j,k)
        enddo
       enddo
      enddo

      if (.not. enable_hydrostatic) then
       do k=2,km-1
        do j=js,je
         do i=2,imt-1
          totke = totke + u(i,j,k,3,tau)**2/2.0*fx*maskW(i,j,k)
         enddo
        enddo
       enddo
      endif

      call global_sum(area)
      call global_sum(totke)
      totke = totke/area

      tbar=0.0
      do k=2,km-1
       do j=js,je
        do i=2,imt-1
          tbar=tbar+b(i,j,k,tau)*fx*maskT(i,j,k)
        enddo
       enddo
      enddo
      call global_sum(tbar)
      tbar = tbar / area

      call set_stamp(stamp,get_current_time())
      if (my_pe==0) print'(a,i5,a,a,a,e12.7,a,f10.8,a,a,i5,a,i5)', 
     &      ' itt=',itt,'  ',stamp,
     &  ' ke=',totke,' m^2/m^2  tbar=',tbar,' kg/m^3',
     &    ' sor2D itts = ',sor2D_itts,' sor3D_itts = ',sor3D_itts
      call sub_flush(6)

      if (totke > 1e20 ) then
        if (my_pe==0) print*,''
        if (my_pe==0) print*,' ERROR: Kinetic energy diverges '
        if (my_pe==0) print*,''
        call halt_stop(' in diagnostics')
      endif

      if (totke +1.0 .eq. totke ) then
        if (my_pe==0) print*,''
        if (my_pe==0) print*,' ERROR: detected Nans in kinetic energy '
        if (my_pe==0) print*,''
        call halt_stop(' in diagnostics')
      endif

      call diag_snap
      endif
      end subroutine diagnose


      subroutine init_snap_cdf
c-----------------------------------------------------------------------
c     initialize NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,i,j,k,n
      integer :: lon_tdim,lon_udim,z_tdim,z_udim,itimedim
      integer :: lat_tdim,lat_udim,id
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16

      call def_grid_cdf('cpflame.cdf')
      iret=nf_open('cpflame.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call ncredf(ncid, iret)
      iret=nf_inq_dimid(ncid,'xt',lon_tdim)
      iret=nf_inq_dimid(ncid,'xu',lon_udim)
      iret=nf_inq_dimid(ncid,'yt',lat_tdim)
      iret=nf_inq_dimid(ncid,'yu',lat_udim)
      iret=nf_inq_dimid(ncid,'zt',z_tdim)
      iret=nf_inq_dimid(ncid,'zu',z_udim)
      iret=nf_inq_dimid(ncid,'Time',itimedim)

c     2 dim variables on t grid
      dims = (/Lon_udim,lat_tdim, z_tdim, iTimedim/)
      id  = ncvdef (ncid,'u', NCFLOAT,4,dims,iret)
      name = 'Zonal velocity         '; unit = 'm/s'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_udim, z_tdim, iTimedim/)
      id  = ncvdef (ncid,'v', NCFLOAT,4,dims,iret)
      name = 'Meridional velocity          '; unit = 'm/s'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, z_udim, iTimedim/)
      id  = ncvdef (ncid,'w', NCFLOAT,4,dims,iret)
      name = 'Vertical velocity         '; unit = 'm/s'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, z_tdim, iTimedim/)
      id  = ncvdef (ncid,'b', NCFLOAT,4,dims,iret)
      name = 'Buoyancy               '; unit = 'm/s^2'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      id  = ncvdef (ncid,'p', NCFLOAT,4,dims,iret)
      name = 'Pressure      '; unit = 'm^2/s^2'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, iTimedim,1/)
      id = ncvdef (ncid,'ps', NCFLOAT,3,dims,iret)
      if (enable_free_surface.or.enable_expl_free_surf) then
       name = 'Surface height'; unit = 'm'
      else
       name = 'Surface pressure'; unit = 'm^2/s^2'
      endif
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, z_udim, iTimedim/)
      if (enable_hydrostatic) then
       id  = ncvdef (ncid,'K_b', NCFLOAT,4,dims,iret)
       name = 'Vertical diffusivity'; unit = 'm^2/s'
       call dvcdf(ncid,id,name,24,unit,16,spval)
      else
       id  = ncvdef (ncid,'psi', NCFLOAT,4,dims,iret)
       name = 'Non-hydrostatic pressure'; unit = 'm^2/s^2'
       call dvcdf(ncid,id,name,24,unit,16,spval)
      endif

      dims = (/Lon_tdim,lat_tdim, iTimedim,1/)
      id = ncvdef (ncid,'surf_flux', NCFLOAT,3,dims,iret)
      name = 'Surface buoyancy flux'; unit = 'm^2/s^3'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, iTimedim,1/)
      id = ncvdef (ncid,'taux', NCFLOAT,3,dims,iret)
      name = 'Zonal wind stress'; unit = 'm^2/s^2'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      dims = (/Lon_tdim,lat_tdim, iTimedim,1/)
      id = ncvdef (ncid,'tauy', NCFLOAT,3,dims,iret)
      name = 'Meridional wind stress'; unit = 'm^2/s^2'
      call dvcdf(ncid,id,name,24,unit,16,spval)

      call ncclos (ncid, iret)

      end subroutine init_snap_cdf




      subroutine diag_snap
c-----------------------------------------------------------------------
c     write to NetCDF snapshot file
c-----------------------------------------------------------------------
      use cpflame_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,n,npe, corner(4), edges(4)
      real :: a(imt,js_pe:je_pe,km)
      integer :: itdimid,ilen,pid,rid,uid,wid,itimeid
      integer :: i,j,js,je,vid,sflid,kbid,id
      real :: fxa,ut,vt
      type(time_type) :: time

      js=max(2,js_pe); je = min(je_pe,jmt-1)

      do npe=0,n_pes-1
       call barrier
       if (my_pe==npe) then

        iret=nf_open('cpflame.cdf',NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)

        iret=nf_inq_varid(ncid,'b',rid)
        iret=nf_inq_varid(ncid,'p',pid)
        iret=nf_inq_varid(ncid,'u',uid)
        iret=nf_inq_varid(ncid,'v',vid)
        iret=nf_inq_varid(ncid,'w',wid)
        iret=nf_inq_varid(ncid,'ps',sflid)
        if (enable_hydrostatic) then
         iret=nf_inq_varid(ncid,'K_b',kbid)
        else
         iret=nf_inq_varid(ncid,'psi',kbid)
        endif
        iret=nf_inq_dimid(ncid,'Time',itdimid)
        iret=nf_inq_dimlen(ncid, itdimid,ilen)
        iret=nf_inq_varid(ncid,'Time',itimeid)
        if (my_pe==0) then
         ilen=ilen+1
         time = current_time-initial_time
         fxa = time%days + time%seconds/86400.
         print*,' writing a snapshot at stamp=',current_stamp,
     &          ' (days since origin : ',fxa,')',
     &          ' (time steps in file : ',ilen,')'
         iret= nf_put_vara_double(ncid,itimeid,ilen,1,fxa)
        endif

        Corner = (/1,js_pe,1,ilen/); !edges  = (/imt,j_blk,km,1/)
        edges  = (/imt,je_pe-js_pe+1,km,1/)

        a=b(:,js_pe:je_pe,:,tau)
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,rid,corner,edges,a)

        a=p_full(:,js_pe:je_pe,:,tau)
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,pid,corner,edges,a)

        a=u(:,js_pe:je_pe,:,1,tau)
        where( maskU(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,uid,corner,edges,a)

        a=u(:,js_pe:je_pe,:,2,tau)
        where( maskV(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,vid,corner,edges,a)

        a=u(:,js_pe:je_pe,:,3,tau)
        where( maskW(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,wid,corner,edges,a)

        if (enable_hydrostatic) then
         a=K_b(:,js_pe:je_pe,:)+K_v
         where( maskW(:,js_pe:je_pe,:) == 0.) a = spval
        else
         a=psi(:,js_pe:je_pe,:)
         where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        endif
        iret= nf_put_vara_double(ncid,kbid,corner,edges,a)

        Corner = (/1,js_pe,ilen,1/); edges  = (/imt,j_blk,1,1/)
        if (enable_free_surface.or.enable_expl_free_surf) then
         a(:,:,km-1)=eta(:,js_pe:je_pe,tau)
        else
         a(:,:,km-1)=p_surf(:,js_pe:je_pe)
        endif
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,sflid,corner,edges,a(:,:,km-1))

        Corner = (/1,js_pe,ilen,1/); edges  = (/imt,j_blk,1,1/)
        iret=nf_inq_varid(ncid,'surf_flux',id)
        a(:,:,km-1)=surf_flux(:,js_pe:je_pe)
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a(:,:,km-1))

        iret=nf_inq_varid(ncid,'taux',id)
        a(:,:,km-1)=surf_tau(:,js_pe:je_pe,1)
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a(:,:,km-1))

        iret=nf_inq_varid(ncid,'tauy',id)
        a(:,:,km-1)=surf_tau(:,js_pe:je_pe,2)
        where( maskT(:,js_pe:je_pe,:) == 0.) a = spval
        iret= nf_put_vara_double(ncid,id,corner,edges,a(:,:,km-1))

        call ncclos (ncid, iret)

       endif
       call barrier
      enddo
      end subroutine diag_snap




      subroutine def_grid_cdf(filename)
c-----------------------------------------------------------------------
c      Define standard grid in netcdf file
c-----------------------------------------------------------------------
      use cpflame_module
      implicit none
#include "netcdf.inc"
      character*(*) filename
      integer :: ncid,iret
      integer :: lon_tdim,lon_udim,z_tdim,z_udim,itimedim
      integer :: lon_tid,lon_uid,z_tid,z_uid,itimeid
      integer :: lat_tdim,lat_udim,lat_uid,lat_tid
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16
      character (len=200) :: text

      ncid = nccre (filename, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      lon_tdim  = ncddef(ncid, 'xt', imt, iret)
      Lon_udim  = ncddef(ncid, 'xu', imt, iret)
      lat_tdim  = ncddef(ncid, 'yt', jmt, iret)
      Lat_udim  = ncddef(ncid, 'yu', jmt, iret)
      z_tdim    = ncddef(ncid, 'zt',  km, iret)
      z_udim    = ncddef(ncid, 'zu',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      Lon_tid  = ncvdef (ncid,'xt',NCFLOAT,1,lon_tdim,iret)
      Lon_uid  = ncvdef (ncid,'xu',NCFLOAT,1,lon_udim,iret)
      Lat_tid  = ncvdef (ncid,'yt',NCFLOAT,1,lat_tdim,iret)
      Lat_uid  = ncvdef (ncid,'yu',NCFLOAT,1,lat_udim,iret)
      z_tid  = ncvdef (ncid,'zt', NCFLOAT,1,z_tdim,iret)
      z_uid  = ncvdef (ncid,'zu', NCFLOAT,1,z_udim,iret)
      itimeid  = ncvdef (ncid,'Time', NCFLOAT,1,itimedim,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'm'
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'm'
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid     '; unit = 'm'
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid     '; unit = 'm'
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
c      name = 'Depth on T grid      '; unit = 'm'
      name = 'Height on T grid      '; unit = 'm'
      call ncaptc(ncid, z_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, z_tid, 'units',     NCCHAR, 16, unit, iret) 
c      name = 'Depth on U grid      '; unit = 'm'
      name = 'Height on U grid      '; unit = 'm'
      call ncaptc(ncid, z_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, z_uid, '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,
     &  '01-JAN-1900 00:00:00', iret)


      write(text,'("CPFLAME Version ",f5.3)') version
#ifdef no_mpp
      text=text(1:len_trim(text))//' without MPP'
#else
      text=text(1:len_trim(text))//' with MPI'
#endif
      iret= nf_put_att_text(ncid,nf_global,'history',
     &                      len_trim(text),text)

#ifndef no_mpp
      iret= nf_put_att_int(ncid,nf_global,'n_pes_j',nf_int,1,n_pes_j)
      iret= nf_put_att_int(ncid,nf_global,'j_blk',nf_int,1,j_blk)
#endif


      call ncendf(ncid, iret)
      iret= nf_put_vara_double(ncid,lon_Tid,1,imt,xt)
      iret= nf_put_vara_double(ncid,lon_uid,1,imt,xu)
      iret= nf_put_vara_double(ncid,lat_Tid,1,jmt,yt)
      iret= nf_put_vara_double(ncid,lat_uid,1,jmt,yu)
      iret= nf_put_vara_double(ncid,z_tid,1,km,zt)
      iret= nf_put_vara_double(ncid,z_uid,1,km,zw)
      call ncclos (ncid, iret)
      end subroutine def_grid_cdf





      subroutine dvcdf(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c     define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character (len=*) :: name, unit
      real :: spval
      real :: vv
#include "netcdf.inc"
      vv=spval
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
c      call ncapt (ncid,ivarid, 'missing_value',NCFLOAT,1,vv,iret)
      call ncapt (ncid,ivarid, 'missing_value',NCDOUBLE,1,vv,iret)
        if (iret.ne.0) print*,nf_strerror(iret)
c      call ncapt (ncid,ivarid, '_FillValue', NCFLOAT, 1,vv, iret)
      call ncapt (ncid,ivarid, '_FillValue', NCDOUBLE, 1,vv, iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf



