#include "options.inc"




      subroutine init_diag
c
c-----------------------------------------------------------------------
c     initialize all diagnostics here
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      if (my_pe==0) then
       print*,' Initialising diagnostics '
#ifdef netcdf_diagnostics
       print*,' writing to NetCDF files '
#else 
       print*,' writing to unformatted files '
#endif
      endif
      if (enable_snapshots)   call init_snap_cdf
      if (enable_overturning) call init_over_cdf
      if (enable_heat_tr)     call init_heat_tr_cdf
                              call init_averages
                              call init_variances
      if (enable_timeseries)  call init_diag_timeserie
      if (enable_diag_press)  call init_diag_press
      if (enable_diag_float)  call diag_float_init
      if (enable_diag_blue)   call diag_blue_init

      if (my_pe==0) print*,' done'
      end subroutine init_diag


      subroutine diag
c
c-----------------------------------------------------------------------
c     call diagnostics here
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none

      if (enable_stability_tests) call stab()
      if (enable_ts_monitor.and. .not. eulerback2) 
     &               call diag_ts_monitor
      if (enable_snapshots .and. snapshot_time_step) 
     &               call diag_snap
      if (enable_overturning .and. snapshot_time_step) 
     &               call diag_over
c      called in tracer, we need adv_fn and diff_fn
c      if (enable_heat_tr .and. snapshot_time_step) 
c     &               call diag_heat_tr
      call averages
      call variances
      if (enable_timeseries .and. snapshot_time_step) 
     &               call diag_timeserie()
      if (enable_diag_press .and. snapshot_time_step)  
     &               call diag_press()

      if (enable_diag_float .and. .not. eulerback1) call diag_float
      if (enable_diag_blue  .and. .not. eulerback1) call diag_blue

      end subroutine diag



      subroutine diag_ts_monitor()
c
c-----------------------------------------------------------------------
c     diagnose some basinwide averages and write to stdout
c-----------------------------------------------------------------------
c
      use spflame_module
      use freesurf_module
      implicit none
      integer i,j,k,n,is,ie,js,je
      real :: totke=0.,tbar(nt),tvar(nt)
      real :: darea,fx,etaq
      character (len=80) :: string

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      tbar=0.; tvar=0.
      do j=js,je
        fx = rho0*0.5*csu(j)*dyu(j)
        do k=1,km
         do i=is,ie
          totke = totke + u(i,k,j,1,tau)**2*fx*dxu(i)
#ifdef partial_cell
     &              *dhu(i,k,j)
#else
     &              *dzt(k)
#endif
          totke = totke + u(i,k,j,2,tau)**2*fx*dxu(i)
#ifdef partial_cell
     &              *dhu(i,k,j)
#else
     &              *dzt(k)
#endif
         enddo
        enddo
        do n=1,nt
         do k=1,km
          do i=is,ie
           darea = dxt(i)*cst(j)*dyt(j)*tmask(i,k,j)
#ifdef partial_cell
     &              *dht(i,k,j)
#else
     &              *dzt(k)
#endif
           tbar(n)=tbar(n)+t(i,k,j,n,tau)*darea
           tvar(n)=tvar(n)+t(i,k,j,n,tau)**2*darea
         enddo
        enddo
       enddo
      enddo
      call global_sum(totke); 
      call global_sum_vec(tbar,nt)
      call global_sum_vec(tvar,nt) 
      totke  = totke  /ucellv
      tbar  = tbar/tcellv
      tvar  = tvar/tcellv - tbar**2

      if (enable_freesurf) then

       etaq=0.
       do j=js,je
         do i=is,ie
           darea = dxt(i)*cst(j)*dyt(j)*tmask(i,1,j)
           etaq=etaq+etat(i,j,3)*darea
         enddo
       enddo
       call global_sum(etaq); 
       etaq=etaq/tcella(1)

       if (my_pe==0)  print*,
     &  ' domain ',sub_domain,
     &  ' : (itt=',itt,' ',current_stamp,')',
     &  ' etab=',etaq,
     &  ' ke=',totke,' tbar=',tbar,' tvar=',tvar

      else
       if (my_pe==0) print*,
c     &  print'(a,i1,a,i8,4a,i3,a,e12.7,a,e12.7,a,e12.7)',
     &  ' domain ',sub_domain,
     &  ' : (itt=',itt,' ',current_stamp,')',
     &  ' itts=',itts_solver,
     &  ' ke=',totke,' tbar=',tbar,' tvar=',tvar
      endif

      end subroutine diag_ts_monitor



