c
c-----------------------------------------------------------------------
c     file contains diagnostic part for BARBI
c-----------------------------------------------------------------------
c
      subroutine diagi()
c=======================================================================
c     Initialization of diagnostics for BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      call def_snap
      end subroutine diagi


      subroutine diag()
c=======================================================================
c     diagnose BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      integer i,j
      real (kind=8) :: fxa,fxb,uarea,tarea
c
      if (snapshot_time_step) then
c
c      compute basin averaged Ekin and Epot
c
       ke=0.; ek=0.
       do j=2,jmt-1
        fxa=  csu(j)*dyu(j)
        fxb=  cst(j)*dyt(j)
        do i=2,imt-1
         uarea=fxa*dxu(i)*umask(i,j)
         tarea=fxb*dxt(i)*tmask(i,j)
         ke=ke+(u(i,j,1,tau)**2+u(i,j,2,tau)**2)/2.*hu(i,j)*uarea
         ek=ek+g(i,j,1,tau)*tarea
        enddo
       enddo
c
c      print that out
c
       write(6,'(a,i5,a,e8.2,a,e8.2,a,i5,a,i5)') 
     &  'itt=',itt
     &  ,' kin.en.=',ke
     &  ,' pot.en.=',ek
     &  ,' psi_itts:',mscans
     &  ,' P_itts:',mscans_p
c
c      write NetCDF snpshot
c
       call snapit()
c
      endif
      end subroutine diag


      subroutine snapit()
c=======================================================================
c     Append a snapshot to NetCDF file
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      real (kind=4) :: var(imt,jmt),spval=-9.9e12,tt,vv
      integer ncid,tauxid,tauyid,gid,psiid,pressid,fid
      integer uid,vid,ubid,vbid,ekid,keid,timeid,timedim
      integer iret,old_mode,len,n,nn
      integer, dimension(4) :: dims, corner, edges
      include "netcdf.inc"

      iret=nf_open(snap_file,NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, old_mode)
c
c     get the id of all variables
c
      iret=nf_inq_varid(ncid,'taux',tauxid)
      iret=nf_inq_varid(ncid,'tauy',tauyid)
      iret=nf_inq_varid(ncid,'G',Gid)
      iret=nf_inq_varid(ncid,'u',uid)
      iret=nf_inq_varid(ncid,'v',vid)
      iret=nf_inq_varid(ncid,'f',fid)
      iret=nf_inq_varid(ncid,'P',pressid)
      iret=nf_inq_varid(ncid,'psi',psiid)
      iret=nf_inq_varid(ncid,'ub',ubid)
      iret=nf_inq_varid(ncid,'vb',vbid)
      iret=nf_inq_varid(ncid,'ek',ekid)
      iret=nf_inq_varid(ncid,'ke',keid)
      iret=nf_inq_dimid(ncid,'Time',timedim)
      iret=nf_inq_varid(ncid,'Time',timeid)
c
c     append time  and zero dim. variables
c
      iret=nf_inq_dimlen(ncid, timedim,len)
      if (iret.ne.0) print*,nf_strerror(iret)
      len=len+1
      tt=float(itt)*dt/60./60./24.
      iret= nf_put_vara_real (ncid,timeid,len,1,tt)
      vv=ke
      iret= nf_put_vara_real (ncid,keid,len,1,vv)
      vv=ek
      iret= nf_put_vara_real (ncid,ekid,len,1,vv)
c
c     append 2-dim variables
c
      Corner = (/1,1,len,1/)
      edges =  (/imt,jmt,1,1/)
c
      var=u(:,:,1,tau)*hu
      where( umask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,uid , corner, edges,var)
c
      var=u(:,:,2,tau)*hu
      where( umask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,vid , corner, edges,var)
c
      var=press(:,:)
      where( tmask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,pressid , corner, edges,var)
c
      var=psi(:,:,1)/1e6
      where( tmask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,psiid , corner, edges,var)
c
      var=wind(:,:,2)*rho0
      where( umask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,tauyid , corner, edges,var)
c
      var=wind(:,:,1)*rho0
      where( umask == 0. ) var = spval
      iret= nf_put_vara_real (ncid,tauxid , corner, edges,var)
c
c     append 3-dim variables
c       write either raw vertical density/vel. modes or 
c       flat bottom baroclinic modes
c
      do n=1,n_order,2

       Corner = (/1,1,n/2+1,len/)
       edges =  (/imt,jmt,1,1/)
       if (write_modes) then
        var=0.
        do nn=1,n_order,2
         var=var+gscal(n/2+1,nn/2+1)*g(:,:,nn,tau)/h**(nn-1)
        enddo
       else
        where (tmask==1) var=g(:,:,n,tau)/h**(n-1)
       endif
       where( tmask == 0. ) var = spval
       iret= nf_put_vara_real (ncid,Gid , corner, edges,var)

       if (write_modes) then
        var=0.
        do nn=1,n_order,2
c        do not know if this is correct but it works for now
         var=var+gscal(n/2+1,nn/2+1)*ub(:,:,1,tau,nn)/hu**(nn-1)
     &             /(nn+1)
        enddo
       else
        where (umask==1) var=ub(:,:,1,tau,n)/hu**(n-1)
       endif
       where( umask == 0. ) var = spval
       iret= nf_put_vara_real (ncid,ubid , corner, edges,var)
c
       if (write_modes) then
        var=0.
        do nn=1,n_order,2
         var=var+gscal(n/2+1,nn/2+1)*ub(:,:,2,tau,nn)/hu**(nn-1)
     &             /(nn+1)
        enddo
       else
        where (umask==1) var=ub(:,:,2,tau,n)/hu**(n-1)
       endif
       where( umask == 0. ) var = spval
       iret= nf_put_vara_real (ncid,vbid , corner, edges,var)
c
      enddo
      call ncclos (ncid, iret)
      end subroutine  snapit


      subroutine def_snap()
c=======================================================================
c     define a netCDF file for output
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      include "netcdf.inc"
      integer, dimension(4) :: dims, corner, edges
      character (len=80) :: name, unit, text
      real (kind=4) :: spval=-9.9e12,var(imt,jmt)
      integer ncid,tauxid,tauyid,gid,psiid,pressid,fid
      integer uid,vid,ubid,vbid,ekid,keid,timeid
      integer lon_tdim,lon_udim,lat_tdim,lat_udim,timedim
      integer lon_tid,lon_uid,lat_tid,lat_uid,topoid,topouid
      integer iret,old_mode,ilen, order_dim,order_id,k
       
      ncid = nccre (snap_file, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, old_mode)
c
c     define dimensions
c
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      order_dim = ncddef(ncid, 'Order',  n_order/2+1, iret)
      Timedim   = ncddef(ncid, 'Time', nf_unlimited, iret)
c
c     define variables
c
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,lon_tdim,iret)
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,lon_udim,iret)
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,lat_tdim,iret)
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,lat_udim,iret)
      order_id = ncvdef (ncid,'Order', NCFLOAT,1,order_dim,iret)
      Timeid   = ncvdef(ncid,'Time',   NCFLOAT,1,timedim,iret)
c
      dims    = (/Lon_tdim,lat_tdim,timedim,0/)
      psiid   = ncvdef (ncid,'psi' , NCFLOAT,3,dims,iret)
      pressid = ncvdef (ncid,'P'   , NCFLOAT,3,dims,iret)
      topoid  = ncvdef (ncid,'topo', NCFLOAT,2,dims,iret)
c
      dims    = (/Lon_udim,lat_udim,timedim,0/)
      fid     = ncvdef (ncid,'f'     , NCFLOAT,2,dims,iret)
      topouid = ncvdef (ncid,'topo_u', NCFLOAT,2,dims,iret)
      uid     = ncvdef (ncid,'u'     , NCFLOAT,3,dims,iret)
      vid     = ncvdef (ncid,'v'     , NCFLOAT,3,dims,iret)
      tauxid  = ncvdef (ncid,'taux'  , NCFLOAT,3,dims,iret)
      tauyid  = ncvdef (ncid,'tauy'  , NCFLOAT,3,dims,iret)

      dims  = (/Lon_tdim,lat_tdim,order_dim,timedim/)
      Gid   = ncvdef (ncid,'G', NCFLOAT,4,dims,iret)

      dims  = (/Lon_udim,lat_udim,order_dim,timedim/)
      ubid  = ncvdef (ncid,'ub', NCFLOAT,4,dims,iret)
      vbid  = ncvdef (ncid,'vb', NCFLOAT,4,dims,iret)

      keid  = ncvdef (ncid,'ke', NCFLOAT,1,timedim,iret)
      ekid  = ncvdef (ncid,'ek', NCFLOAT,1,timedim,iret)
c
c     assign attributes
c
      name = 'Longitude on T grid';k=len_trim(name)
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, k, name, iret) 
      unit = 'degrees_W ';k=len_trim(unit)
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, k, unit, iret) 
c      
      name = 'Longitude on U grid';k=len_trim(name)
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, k, name, iret) 
      unit = 'degrees_W';k=len_trim(unit)
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, k, unit, iret) 
c      
      name = 'Latitude on T grid';k=len_trim(name)
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, k, name, iret) 
      unit = 'degrees_N ';k=len_trim(unit)
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, k, unit, iret) 
c      
      name = 'Latitude on U grid';k=len_trim(name)
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, k, name, iret) 
      unit = 'degrees_N';k=len_trim(unit)
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, k, unit, iret) 
c      
      name = 'Order';k=len_trim(name)
      call ncaptc(ncid, order_id, 'long_name', NCCHAR, k, name, iret) 
      unit = 'index';k=len_trim(unit)
      call ncaptc(ncid, order_id, 'units',     NCCHAR, k, unit, iret) 
c      
      name = 'Time';k=len_trim(name)
      call ncaptc(ncid,Timeid, 'long_name', NCCHAR, k, name, iret) 
      unit = 'days';k=len_trim(unit)
      call ncaptc(ncid,Timeid, 'units',     NCCHAR, k, unit, iret) 
      call ncaptc(ncid,Timeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c
      name = 'Zonal Transport'
      unit = 'm^2/s'
      call dvcdf(ncid,uid,name,unit,spval)

      name = 'Meridional Transport'
      unit = 'm^2/s'
      call dvcdf(ncid,vid,name,unit,spval)

      name = 'Meridional Wind Stress'
      unit = 'N/m^2'
      call dvcdf(ncid,tauyid,name,unit,spval)
c
      name = 'Zonal Wind Stress'
      unit = 'N/m^2'
      call dvcdf(ncid,tauxid,name,unit,spval)
c
      if (write_modes) then
       name = 'Baroclinic (flat bottom) density modes'
      else
       name = 'Vertical density Moments'
      endif
      unit = 'm^3/s^2'
      call dvcdf(ncid,Gid,name,unit,spval)
c
      name = 'Bottom pressure'
      unit = 'bar'
      call dvcdf(ncid,pressid,name,unit,spval)
c
      name = 'Coriolis Parameter'
      unit = '1/s'
      call dvcdf(ncid,fid,name,unit,spval)
c
      name = 'Barotropic Streamfunction'
      unit = 'Sv'
      call dvcdf(ncid,psiid,name,unit,spval)
c
      name = 'Zonal barocl. velocity moments'
      unit = 'm^4/s'
      call dvcdf(ncid,ubid,name,unit,spval)
c      
      name = 'Meridional barocl. velocity moments'
      unit = 'm^4/s'
      call dvcdf(ncid,vbid,name,unit,spval)
c      
      name = 'Depth of T-grid points'
      unit = 'm'
      call dvcdf(ncid,topoid,name,unit,spval)
c      
      name = 'Depth of U-grid points'
      unit = 'm'
      call dvcdf(ncid,topouid,name,unit,spval)
c      
      name = 'Basin averaged barotr. kinetic energy'
      unit = 'm^2/s^2'
      call dvcdf(ncid,keid,name,unit,spval)
c      
      name = 'Basin averaged potential energy'
      unit = 'm^3/s^2'
      call dvcdf(ncid,ekid,name,unit,spval)
c      
      text='BARBI I Snapshot'
      call ncaptc(ncid, NCGLOBAL, 'title',   NCCHAR, 80, text, iret)
      text='converted to NetCDF during model integration '
      call ncaptc(ncid, NCGLOBAL, 'history', NCCHAR, 80, text, iret)
      write(text,'("Horizontal viscosity = ",e8.2," m^2/s"'//
     &     '," vertical viscosity = ",e8.2," m^2/s")') Am,Av
      call ncaptc(ncid, NCGLOBAL, 'Viscosity', NCCHAR, 80,text,iret)
      write(text,'("Horizontal diffusivity = ",e8.2," m^2/s")') Ah
      call ncaptc(ncid, NCGLOBAL, 'Diffusivity', NCCHAR,80,text,iret)
      write(text,'("Time step = ",e8.2," s, N_0 = ",e8.2," 1/s")') dt,N0
      call ncaptc(ncid, NCGLOBAL, 'Time_step', NCCHAR,80,text,iret)
      call ncendf(ncid, iret)
c
      var(1:imt,1)=xt
      call ncvpt(ncid, Lon_tid, 1, imt, var, iret)
      var(1:imt,1)=xu
      call ncvpt(ncid, Lon_uid, 1, imt, var, iret)
      var(1:jmt,1)=yt
      call ncvpt(ncid, Lat_tid, 1, jmt, var, iret)
      var(1:jmt,1)=yu
      call ncvpt(ncid, Lat_uid, 1, jmt, var, iret)
      corner(1) = 1
      edges(1) = n_order/2+1
      do k=1,n_order/2+1
       var(k,1)=k*2-1
      enddo
      call ncvpt(ncid, order_id, corner, edges, var, iret)
      Corner = 1
      edges = (/imt,jmt,1,1/)
      var=h
      call ncvpt(ncid, topoid, corner, edges, var, iret)
      var=hu
      call ncvpt(ncid, topouid, corner, edges, var, iret)
      var= cori(:,:,1)
      where( umask == 0. ) var = spval
      call ncvpt(ncid, fid, corner, edges, var, iret)
      call ncclos (ncid, iret)


      end subroutine def_snap



      subroutine dvcdf(ncid,varid,name,unit,spval)
c=======================================================================
c     define some standard attributes of variable varid 
c     in NetCDF file ncid 
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      implicit none
      integer ncid,varid,iret
      character (len=*) :: name, unit
      real (kind=4) spval
      include "netcdf.inc"
      call ncaptc(ncid,varid, 'long_name', NCCHAR,
     &            len_trim(name) , name, iret) 
      call ncaptc(ncid,varid, 'units',     NCCHAR,
     &            len_trim(unit), unit, iret) 
      call ncapt (ncid,varid, 'missing_value',NCFLOAT,1,spval,iret)
      call ncapt (ncid,varid, '_FillValue', NCFLOAT, 1,spval, iret)
      end subroutine dvcdf
