! source file: /Users/oschlies/UVIC/master/source/mom/diag.F
      subroutine diag (joff, js, je, is, ie)

!=======================================================================
!     calculate diagnostics

!     input:

!      joff   = offset between row j in the MW and latitude jrow on disk
!      js     = starting row for calculations
!      je     = ending row for calculations
!      is     = starting longitude index for calculations
!      ie     = ending longitude index for calculations
!=======================================================================

      include "param.h"
      include "coord.h"
      include "cregin.h"
      include "diag.h"
      include "diaga.h"
      include "docnam.h"
      include "grdvar.h"
      include "iounit.h"
      include "isopyc.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"
      include "vmixc.h"
      include "levind.h"
      include "emode.h"
      include "npzd.h"

      character(120) :: fname
      integer ntrec

      real time, tmp_t(imt,km,nt), tmp_stf(imt,nt)

      dimension vbarx(km)
      dimension aibuf(imt,km)

!-----------------------------------------------------------------------
!     bail out if starting row exceeds ending row
!-----------------------------------------------------------------------

      if (js .gt. je) return

!-----------------------------------------------------------------------
!     limit longitudes
!-----------------------------------------------------------------------

      istrt  = max(2,is)
      iend   = min(imt-1,ie)

      if (tsiperts .and. .not. euler2 .and. joff .eq. 0)
     &  nv_otsf = nv_otsf + 1
      if (tsiperts .and. .not. euler2 .and. joff .eq. 0)
     &  nt_slh = nt_slh + 1
      do j=js,je
        jrow = joff + j
!-----------------------------------------------------------------------
!       diagnostic: accumulate "tau" data for time step integrals
!-----------------------------------------------------------------------

        if (tsiperts .and. .not. euler2) then

          if (jrow .ge. jsot .and. jrow .le. jeot) then
            if (mrot .gt. 0 .and. mrot .le. nhreg) then
              do i=2,imtm1
                if (mskhr(i,jrow) .eq. mrot) then
                  do k=1,kmu(i,jrow)
                    v_otsf(jrow,k) = v_otsf(jrow,k) + u(i,k,j,2,tau)*
     &                               dxu(i)
                  enddo
                endif
              enddo
            else
              do i=isot1,ieot1
                do k=1,kmu(i,jrow)
                  v_otsf(jrow,k) = v_otsf(jrow,k) + u(i,k,j,2,tau)*
     &                             dxu(i)
                enddo
              enddo
              do i=isot2,ieot2
                do k=1,kmu(i,jrow)
                  v_otsf(jrow,k) = v_otsf(jrow,k) + u(i,k,j,2,tau)*
     &                             dxu(i)
                enddo
              enddo
            endif
          endif
          do i=2,imtm1
            do k=1,kmt(i,jrow)
              t_slh(i,jrow,k,1) = t_slh(i,jrow,k,1) + t(i,k,j,1,tau)
              t_slh(i,jrow,k,2) = t_slh(i,jrow,k,2) + t(i,k,j,2,tau)
            enddo
          enddo

        endif

!-----------------------------------------------------------------------
!       diagnostic: accumulate "tau" data for time means
!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

        if (timavgperts .and. .not. euler2) then
          if (istrt .ne. 2 .and. iend .ne. imt-1) then
            write (stdout,*) '=>Error: istrt = ',istrt,' and iend ='
     &,     iend,' are not allowed when calling "avgvar"'
            stop '=>diag'
          else
            call avgvar (j, jrow, adv_vbt(1,1,j), u(1,1,1,1,tau)
     &,                  t(1,1,1,1,tau), stf, smf, mapt)
          endif
        endif

!-----------------------------------------------------------------------
!       diagnostic: write instantaneous sample of MOM data
!       based on code by: R. C. Pacanowski and A. Rosati
!-----------------------------------------------------------------------

        if (snapts .and. .not. euler2) then
          avgper = 0.
          is = 1
          ie = imt
          tmp_t(1:imt,1:km,1:nt) = t(1:imt,1:km,j,1:nt,tau)
          tmp_stf(1:imt,1:nt) = stf(1:imt,j,1:nt)
          time = relyr + year0
          if (jrow .eq. 2) call def_snap
          call def_snap_mom (fname)
          call mom_snap_out (fname, is, ie, jrow, jrow, imt, jmt, km, nt
     &,                      xt, yt, zt, xu, yu, zw, dxt, dyt, dzt, dxu
     &,                      dyu, dzw, avgper, time, stamp, mapt, tmp_t
     &,                      u(is,1,j,1,tau), u(is,1,j,2,tau)
     &,                      adv_vbt(is,1,j), tmp_stf, smf(is,j,1)
     &,                      smf(is,j,2)
     &,                      adv_vetiso(is,1,j), adv_vntiso(is,1,j)
     &,                      adv_vbtiso(is,1,j)
     &,                      dc14(is,1,j)
     &,                      diff_cbt(is,1,j)
     &,                      kpzd, rnpp(is,1,j), rgraz(is,1,j)
     &,                      rmorp(is,1,j), rmorpt(is,1,j)
     &,                      rmorz(is,1,j), rremi(is,1,j)
     &,                      rexcr(is,1,j), rexpo(is,1,j)
     &,                      rnpp_D(is,1,j), rgraz_D(is,1,j)
     &,                      rmorp_D(is,1,j), rnfix(is,1,j)
     &,                      rdeni(is,1,j)
     &,                      totalk(is,j), vdepth(is,j), pe(is,j)
     &,                      psi(is,jrow,1)
     &,                      kmt(is,jrow), mskhr(is,jrow)
     &,                      tmask(is,1,j), umask(is,1,j)
     &,                      tlat(is,jrow), tlon(is,jrow)
     &,                      ulat(is,jrow), ulon(is,jrow), ntrec)

          if (jrow .eq. 2)
     &      write (*,'(a,i4,a,a,a,i10,a,a)') '=> Ocn snapshot #', ntrec
     &,        ' written to ',trim(fname),' on ts = ',itt, ', ', stamp
        endif

!-----------------------------------------------------------------------
!       diagnostic: compute stability diagnostics
!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

        if (stabts .and. eots) then
          if (istrt .ne. 2 .and. iend .ne. imt-1) then
            write (stdout,*) '=>Error: istrt = ',istrt,' and iend ='
     &,     iend,' are not allowed when calling "stab"'
            stop '=>diag'
          else
            call stab (j, jrow)
          endif
        endif

!-----------------------------------------------------------------------
!       construct meridional overturning of mass
!     based on code by: R. C. Pacanowski
!-----------------------------------------------------------------------

        if (jrow .lt. jmtm1 .and. vmsfts .and. eots) then
          do k=1,km
            vbarx(k) = c0
          enddo

          do k=1,km
            do i=istrt,iend
              vbarx(k) = vbarx(k) + u(i,k,j,2,tau)*csu(jrow)*dxu(i)
            enddo
            if (k .eq. 1) then
              vmsf(jrow,k) = vbarx(k)*dzt(k)
            else
              vmsf(jrow,k) = vmsf(jrow,k-1) + vbarx(k)*dzt(k)
            endif
          enddo
        endif

      enddo

      return
      end