      subroutine tracer (joff, js, je, is, ie)

#if defined O_mom
!=======================================================================
!     compute tracers at "tau+1" for rows js through je in the MW.

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW
!=======================================================================

      implicit none

      character(120) :: fname, new_file_name

      integer istrt, iend, i, k, j, ip, kr, jq, n, jp, jrow, iou, js
      integer je, limit, joff, is, ie, kmx, m, kb, idiag, index
      integer it(10), iu(10), ib(10), ic(10), nfnpzd, mfnpzd, mxfnpzd
      integer id_time, id_xt, id_yt, id_zt,fe_jlo,fe_m,fe_k,fe_n

      parameter (fe_n = 14)

      logical inqvardef, exists

      real rctheta, declin, gl, impo, expo, npp, time
      real remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac
      real graz_Det, graz_Z, avej, avej_D, gmax, no3P, po4P, po4_D
      real phin, dz, prca, dprca, nud, bct, tap, fo2, so2, ai, hi, hs
      real npp_D, graz_D, morp_D, no3flag, deni, nfix, felimit
      real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r
      real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso
      real adv_tziso, diff_tx, diff_ty, diff_tz, zmax, cont, drho
      real drhom1, wt, ahbi_cstr, ahbi_csu_dyur, gamma, rrstd, fy, fyz
      real fe_dy, fe_conc, fe_x(fe_n), fe_y(fe_n), bctz, felimit_D

      real expofe, impofe, feorgads, remife, thetamax, deffe, fecol
      real thetachl, chl, chl_D, feprime, fesed, bfe, sgb

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "accel.h"
      include "coord.h"
      include "cregin.h"
      include "csbc.h"
      include "emode.h"
      include "grdvar.h"
      include "hmixc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "timeavgs.h"
      include "tmngr.h"
      include "vmixc.h"
# if defined O_save_convection || defined O_carbon_14
      include "diaga.h"
# endif
# if defined O_ice
#  if defined O_ice_cpts
      include "cpts.h"
#  endif
      include "ice.h"
# endif
# if defined O_npzd
#  if defined O_embm
      include "atm.h"
#  endif
      include "npzd.h"
# endif
# if defined O_carbon_fnpzd
      include "calendar.h"
# endif

      parameter (istrt=2, iend=imt-1)

      real twodt(km)
# if defined O_plume
      real  work(imt,km,jsmw:jemw)
# endif
# if defined O_npzd
      real snpzd(ntnpzd), tnpzd(ntnpzd)
# endif
# if defined O_npzd || defined O_carbon_14
      real src(imt,km,jsmw:jemw,nsrc)
# endif
# if defined O_carbon_fnpzd
      real tmpijk(imtm2,jmtm2,km), tmpi(imtm2), tmpj(jmtm2)
      real, allocatable :: fnpzd(:,:,:,:)
      save fnpzd, nfnpzd, mfnpzd, mxfnpzd
# endif
# if defined O_isopycmix
      include "isopyc.h"
# endif
      include "fdift.h"

# if defined O_carbon_fnpzd
!-----------------------------------------------------------------------
!     create file for writing npzd 3d fluxes of carbon and alkalinity
!-----------------------------------------------------------------------
      if (.not. allocated (fnpzd)) then
        allocate ( fnpzd(imt,jmt,km,2) )
        fnpzd(:,:,:,:) = 0.
!       initialize time record counter
        mfnpzd = 0
!       initialize records written counter
        nfnpzd = 0
!       set the repeating time interval to one year
        mxfnpzd = nint(daylen*yrlen/dtts)
        fname = new_file_name ("O_fnpzd.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) then
!         limit time interval to what exists
          call openfile (fname, iou)
          call getdimlen ('time', iou, nfnpzd)
          mxfnpzd = min(mxfnpzd, nfnpzd)
#  if defined O_carbon
          if (inqvardef('O_fnpzd1', iou)) nfnpzd = mxfnpzd
#  endif
#  if defined O_npzd_alk
          if (inqvardef('O_fnpzd2', iou)) nfnpzd = mxfnpzd
#  endif
        else
!         define file if not there
          call openfile (fname, iou)
          call redef (iou)
          call defdim ('time', iou, 0, id_time)
          call defdim ('longitude', iou, imtm2, id_xt)
          call defdim ('latitude', iou, jmtm2, id_yt)
          call defdim ('depth', iou, km, id_zt)
          it(1) = id_time
          call defvar ('time', iou, 1, it, c0, c0, 'T', 'D'
     &,     'time', 'time', 'years since 0-1-1')
          call putatttext (iou, 'time', 'calendar', calendar)
          it(1) = id_xt
          call defvar ('longitude', iou, 1, it, c0, c0, 'X', 'D'
     &,     'longitude', 'longitude', 'degrees')
          it(1) = id_yt
          call defvar ('latitude', iou, 1, it, c0, c0, 'Y', 'D'
     &,     'latitude', 'latitude', 'degrees')
          it(1) = id_zt
          call defvar ('depth', iou, 1, it, c0, c0, 'Z', 'D'
     &,     'depth', 'depth', 'm')
          it(1) = id_xt
          it(2) = id_yt
          it(3) = id_zt
          it(4) = id_time
#  if defined O_carbon
          call defvar ('O_fnpzd1', iou, 4, it, c0, c0, ' ', 'D'
     &,       ' ', ' ', ' ')
#  endif
#  if defined O_npzd_alk
          call defvar ('O_fnpzd2', iou, 4, it, c0, c0, ' ', 'D'
     &,       ' ', ' ', ' ')
#  endif
          call enddef (iou)
          ib(:) = 1
          ic(:) = imtm2
          tmpi(1:imtm2) = xt(2:imtm1)
          call putvara ('longitude', iou, imtm2, ib, ic, xt, c1, c0)
          ib(:) = 1
          ic(:) = jmtm2
          tmpj(1:jmtm2) = yt(2:jmtm1)
          call putvara ('latitude', iou, jmtm2, ib, ic, yt, c1, c0)
          ib(:) = 1
          ic(:) = km
          call putvara ('depth', iou, km, ib, ic, zt, c1, c0)
        endif
      endif
# endif
!-----------------------------------------------------------------------
!     bail out if starting row exceeds ending row
!-----------------------------------------------------------------------

      if (js .gt. je) return

!-----------------------------------------------------------------------
!     limit the longitude indices based on those from the argument list
!     Note: this is currently bypassed. istrt and iend are set as
!           parameters to optimize performance
!-----------------------------------------------------------------------

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

!-----------------------------------------------------------------------
!     build coefficients to minimize advection and diffusion computation
!-----------------------------------------------------------------------

#if defined O_fourth_order_tracer_advection || defined O_fct || defined O_quicker || defined O_pressure_gradient_average || defined O_biharmonic || defined O_isopycmix
      limit = min(je+1+joff,jmt) - joff
      do j=js,limit
# else
      do j=js,je
# endif
        jrow = j + joff
        do i=istrt-1,iend
          cstdxtr(i,j)    = cstr(jrow)*dxtr(i)
          cstdxt2r(i,j)   = cstr(jrow)*dxtr(i)*p5
          cstdxur(i,j)    = cstr(jrow)*dxur(i)
# if defined O_consthmix && !defined O_bryan_lewis_horizontal && !defined O_biharmonic
          ah_cstdxur(i,j) = diff_cet*cstr(jrow)*dxur(i)
# endif
        enddo
      enddo
# if defined O_plume

!-----------------------------------------------------------------------
!     find depth of penetration for the llnl simple plume model
!-----------------------------------------------------------------------

!     set parameters
      zmax = 160.0e2 ! use constant depth penetration (if cont=0)
!      cont = 0.0    ! use rho contrast for penetration (if cont>0)
      cont = 0.4e-3  ! 0.4 kg/m3 = 0.4e-3g/cm3

      do j=js,je
        jrow = j + joff
        do i=is,ie
          subz(i,jrow) = zmax
        enddo
        if (cont .gt. c0) then
          call state_ref (t(1,1,1,1,tau), t(1,1,1,2,tau)
     &,                   work(1,1,jsmw), j, j, 1, imt, 1)
          do i=is,ie
            kmx = kmt(i,jrow)
            subz(i,jrow) = zw(1)
            do k=2,kmx
              drho = work(i,k,j) - work(i,1,j)
              drhom1 = work(i,k-1,j) - work(i,1,j)
              if (drho .ge. cont .and. drhom1 .lt. cont)
     &          subz(i,jrow) = zt(k-1) + (zt(k) - zt(k-1))
     &                      *(cont - drhom1)/(drho - drhom1)
            enddo
            if (drho .le. cont) subz(i,jrow) = zw(kmx)
          enddo
        endif

!-----------------------------------------------------------------------
!     calculate "3d stf" multiplier for the llnl simple plume model
!-----------------------------------------------------------------------

        do i=is,ie
          if (kmt(i,jrow) .gt. 0) then
            subz(i,jrow) = max(min(subz(i,jrow),zw(kmt(i,jrow))),zw(1))
            work(i,1,j) = c1/subz(i,jrow)
            do k=2,km
              wt =  min(max(c0,(subz(i,jrow) - zw(k-1))),dzt(k))*dztr(k)
              work(i,k,j) = wt*work(i,1,j)
            enddo
          else
            do k=1,km
              work(i,k,j) = c0
            enddo
          endif
        enddo
      enddo
# endif
# if defined O_npzd

!-----------------------------------------------------------------------
!     calculation of biological interactions
!-----------------------------------------------------------------------

      declin = sin((mod(relyr,1.) - 0.22)*2.*pi)*0.4   ! declination

      do k=1,km
        twodt(k) = c2dtts*dtxcel(k)
        nbio(k) = twodt(k)/dtnpzd
        dtbio(k) = twodt(k)/nbio(k)
        rdtts(k) = 1./twodt(k)
        rnbio(k) = 1./nbio(k)
      enddo


# if defined O_npzd_iron
      tap = 2.*par   ! alpha is iron dependent and calculated in npzd_src
# else
      tap = 2.*alpha*par
# endif

      do j=js,je
        jrow = j + joff

        do i=is,ie

          if (kmt(i,jrow) .gt. 0) then

#  if defined O_ice
#   if defined O_ice_cpts
            ai = 0.
            hi = 0.
            hs = 0.
            do n=1,ncat
              ai =  ai + A(i,jrow,2,n)
              hi =  hi + heff(i,jrow,2,n)
              hs =  hs + hseff(i,jrow,2,n)
            enddo
#   else
            ai = aice(i,jrow,2)
            hi = hice(i,jrow,2)
            hs = hsno(i,jrow,2)
#   endif
#  else
            ai = 0.
            hi = 0.
            hs = 0.
#  endif
!           calculate day fraction and incoming solar
!           angle of incidence = lat - declin, refraction index = 1.33
            rctheta = max(-1.5, min(1.5, tlat(i,jrow)/radian - declin))
            rctheta = kw/sqrt(1. - (1. - cos(rctheta)**2.)/1.33**2.)
            dayfrac = min( 1., -tan(tlat(i,jrow)/radian)*tan(declin))
            dayfrac = max(1e-12, acos(max(-1., dayfrac))/pi)
#  if defined O_embm
            swr = dnswr(i,jrow)*1e-3*(1. + ai*(exp(-ki*(hi + hs)) - 1.))
#  else
            swr = 200.
#  endif
            expo = 0.0
            impo = 0.0
# if defined O_npzd_iron
            expofe = 0.0
            impofe = 0.0
# endif
            phin = 0.0 ! integrated phytoplankton
            prca = 0.0 ! integrated production of calcite

            kmx = min(kmt(i,jrow), kpzd)
            do k=1,kmx

!-----------------------------------------------------------------------
!             limit tracers to positive values
!-----------------------------------------------------------------------

              tnpzd(1) = max(t(i,k,j,ipo4,taum1), trcmin)
              tnpzd(2) = max(t(i,k,j,iphyt,taum1), trcmin)
              tnpzd(3) = max(t(i,k,j,izoop,taum1), trcmin)
              tnpzd(4) = max(t(i,k,j,idetr,taum1), trcmin)
#  if defined O_npzd_nitrogen
              tnpzd(5) = max(t(i,k,j,ino3,taum1), trcmin)
              tnpzd(6) = max(t(i,k,j,idiaz,taum1), trcmin)
#  endif
#  if defined O_npzd_iron
              tnpzd(7) = max(t(i,k,j,idfe,taum1), trcmin)
              tnpzd(8) = max(t(i,k,j,idetrfe,taum1), trcmin)
#  endif

              swr = swr*exp(-kc*phin)

# if defined O_npzd_nitrogen
              phin = phin + (tnpzd(6) + tnpzd(2))*dzt(k)
# else
              phin = phin + tnpzd(2)*dzt(k)
# endif
              gl = tap*swr*exp(ztt(k)*rctheta)
              impo = expo*dztr(k)
# if defined O_npzd_iron
              impofe = expofe*dztr(k)
# endif
              bct = bbio**(cbio*t(i,k,j,itemp,taum1))
# if defined O_zoop_graz_upper_temp_limit
              if (t(i,k,j,itemp,taum1).gt.20) then
                 bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1))
     &                  *bbio**(cbio*20)
              else
                 bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1))
     &                  *bct
              endif
# else
                 bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1))
     &                  *bct
# endif
#  if defined O_npzd_fe_limitation
              if (k.le.3) then
! create x (time) and y (data) arrays for fe interpolation
! first zero them out just in case
                 fe_y(:) = 0
                 fe_x(:) = 0
! write in values for days 0 and 365 for interp. bounds
                 fe_y(1) = fe_dissolved(i,j,k,1)
                 fe_x(1) = 0
                 fe_y(14) = fe_dissolved(i,j,k,12)
                 fe_x(14) = 365
! now create the rest of the array based on the BLING
! data which is from day 16 of each month
                 do m=2,13
                    fe_y(m) = fe_dissolved(i,j,k,m-1)
                    if (m.eq.2) then
                       fe_x(2) = 16
                    else if (m.eq.3) then
                       fe_x(3) = 47
                    else if (m.eq.4) then
                       fe_x(4) = 75
                    else if (m.eq.5) then
                       fe_x(5) = 106
                    else if (m.eq.6) then
                       fe_x(6) = 136
                    else if (m.eq.7)  then
                       fe_x(7) = 167
                    else if (m.eq.8) then
                       fe_x(8) = 197
                    else if (m.eq.9) then
                       fe_x(9) = 228
                    else if (m.eq.10) then
                       fe_x(10) = 259
                    else if (m.eq.11) then
                       fe_x(11) = 289
                    else if (m.eq.12) then
                       fe_x(12) = 320
                    else if (m.eq.13) then
                       fe_x(13) = 350
                    endif
                 enddo 
                    fe_jlo = 2
                    fe_m = 4
! find fe data day index
                    call hunt (fe_x,fe_n,dayoyr,fe_jlo)
! initialize the fe data array at the right day 
                    fe_k = min(max(fe_jlo-(fe_m-1)/2,1),fe_n+1-fe_m)
! interpolate the fe data, note this does not use the whole array
                    call polint (fe_x(fe_k),fe_y(fe_k)
     &,                          fe_m,dayoyr,fe_conc,fe_dy)
! calculate the fe limitation term                
                    felimit = fe_conc/(kfe + fe_conc)
                    felimit_D = fe_conc/(kfe_D + fe_conc)
                 else
                    felimit = 1
                    felimit_D = 1
                 end if
# endif
#  if defined O_npzd_o2
!             decrease remineralisation rate in oxygen minimum zone
              nud = nud0*(0.65+0.35*tanh(t(i,k,j,io2,taum1)*1000.-6.))
#  else
              nud = nud0
#  endif

!-----------------------------------------------------------------------
!             call the npzd model
!-----------------------------------------------------------------------
              call npzd_src (tnpzd, nbio(k), dtbio(k), gl, bct, impo
     &,                      dzt(k), dayfrac, wd(k), rkwz(k), nud
     &,                      snpzd, expo, graz, morp, morz, graz_Det
     &,                      graz_Z
#  if defined O_save_npzd
     &,                      npp, morpt, remi, excr
#   if defined O_npzd_nitrogen
     &,                      npp_D, graz_D, morp_D, nfix
#   endif
#   if defined O_npzd_extra_diagnostics
     &,                      avej, avej_D, gmax, no3P, po4P, po4_D
#   endif
#  endif
#  if defined O_npzd_fe_limitation
     &,                    felimit, felimit_D
#  endif
     &,                    bctz
# if defined O_npzd_iron
     &,                    expofe, impofe, remife
# if defined O_npzd_iron_diagnostics
     &,                    feorgads, thetamax, deffe, feprime
     &,                    fecol
# endif
     &,                    t(i,k,j,io2,taum1)
#if defined O_npzd_chl
     &,                    chl, thetachl
#if defined O_npzd_nitrogen
     &,                    chl_D
# endif
# endif
# endif
     &                      )
! These are source/sink terms 
              snpzd(1:4) = snpzd(1:4)*rdtts(k)
#  if defined O_npzd_nitrogen
              snpzd(5:6) = snpzd(5:6)*rdtts(k)
#  endif
#  if defined O_npzd_iron
              snpzd(7:8) = snpzd(7:8)*rdtts(k)
              expofe=expofe*rnbio(k)
#  endif
              expo = expo*rnbio(k)
#  if defined O_save_npzd
              rexpo(i,k,j) = expo
              rgraz(i,k,j) = graz*rnbio(k) 
              rgraz_Det(i,k,j) = graz_Det*rnbio(k)
              rgraz_Z(i,k,j) = graz_Z*rnbio(k)
              rmorp(i,k,j) = morp*rnbio(k)
              rmorz(i,k,j) = morz*rnbio(k)
              rnpp(i,k,j) = npp*rnbio(k)
              rmorpt(i,k,j) = morpt*rnbio(k)
              rremi(i,k,j) = remi*rnbio(k)
              rexcr(i,k,j) = excr*rnbio(k)
#   if defined O_npzd_nitrogen
              rnpp_D(i,k,j) = npp_D*rnbio(k)
              rgraz_D(i,k,j) = graz_D*rnbio(k)
              rmorp_D(i,k,j) = morp_D*rnbio(k)
              rnfix(i,k,j) = nfix*rnbio(k)
#   endif
#   if defined O_npzd_iron
              rexpofe(i,k,j) = expofe
              rremife(i,k,j) = remife*rnbio(k)
#   if defined O_npzd_iron_diagnostics
              rfeorgads(i,k,j) = feorgads*rnbio(k)
              rdeffe(i,k,j) = deffe*rnbio(k)
              rfeprime(i,k,j) = feprime*rnbio(k)
              rfecol(i,k,j) = fecol*rnbio(k)
#   endif
#  if defined O_npzd_chl
              rchl(i,k,j) = chl*rnbio(k)
#  if defined O_npzd_nitrogen
              rchl_D(i,k,j) = chl_D*rnbio(k)
#  endif
#  endif
#   endif
#   if defined O_npzd_extra_diagnostics
              ravej(i,k,j) = avej*rnbio(k)
              ravej_D(i,k,j) = avej_D*rnbio(k)
              rgmax(i,k,j) = gmax*rnbio(k)
              rno3P(i,k,j) = no3P*rnbio(k)
              rpo4P(i,k,j) = po4P*rnbio(k)
              rpo4_D(i,k,j) = po4_D*rnbio(k)
#   endif
#  endif
!-----------------------------------------------------------------------
!             calculate detritus at the bottom and remineralize
!-----------------------------------------------------------------------

#  if defined O_npzd_subgridbathy
              sgb = sg_bathy(i,j,k)
              if (sgb .gt. 0) then
#  if defined O_sed
                if (addflxo .and. eots)
     &            sbc(i,jrow,irorg) = sbc(i,jrow,irorg)
     &                              + expo*redctn*twodt(1)*dzt(k)
#  endif
#  if defined O_save_npzd
                rremi(i,k,j) = rremi(i,k,j) + expo*sgb
#  endif
                snpzd(1) = snpzd(1) + redptn*expo*sgb
#  if defined O_npzd_nitrogen
                snpzd(5) = snpzd(5) + expo*sgb
#  endif
#  if defined O_npzd_iron
                fesed=fetopsed*bct*redptn*expo*sgb
                bfe=fesed
                snpzd(7) = snpzd(7) + fesed
                if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then
                  snpzd(7) = snpzd(7) + expofe*sgb
                  bfe=bfe+expofe*sgb
                  expofe = expofe - sgb * expofe
                endif
                rremife(i,k,j) = rremife(i,k,j) + bfe
#  if defined O_npzd_iron_diagnostics
                rfesed=fesed
                rbfe=bfe
#  endif
#  endif
                expo = expo - sgb * expo
              endif

#  else

              if (k .eq. kmt(i,jrow)) then
#  if defined O_sed
                if (addflxo .and. eots)
     &            sbc(i,jrow,irorg) = sbc(i,jrow,irorg)
     &                              + expo*redctn*twodt(1)*dzt(k)
#  endif
#  if defined O_save_npzd
                rremi(i,k,j) = rremi(i,k,j) + expo
#  endif
                snpzd(1) = snpzd(1) + redptn*expo
#  if defined O_npzd_nitrogen
                snpzd(5) = snpzd(5) + expo
#  endif
#  if defined O_npzd_iron
                fesed=fetopsed*bct*redptn*expo
                bfe=fesed
                snpzd(7) = snpzd(7) + fesed
                if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then
                  snpzd(7) = snpzd(7) + expofe
                  bfe=bfe+expofe
                endif
                rremife(i,k,j) = rremife(i,k,j) + bfe
#  if defined O_npzd_iron_diagnostics
                rfesed=fesed
                rbfe=bfe
#  endif
#  endif
              endif

#  endif

!-----------------------------------------------------------------------
!             set source/sink terms
!-----------------------------------------------------------------------

              src(i,k,j,ispo4) = snpzd(1)
              src(i,k,j,isphyt) = snpzd(2)
              src(i,k,j,iszoop) = snpzd(3)
              src(i,k,j,isdetr) = snpzd(4)
#  if defined O_npzd_nitrogen
              src(i,k,j,isno3) = snpzd(5)
              src(i,k,j,isdiaz) = snpzd(6)
#  endif
#  if defined O_npzd_iron
              if ( k .eq. 1 ) then
                snpzd(7) = snpzd(7) +    ! comes in mol/m2/s gets converted to mmol/m3/s
     &                   sbc(i,j,idfeadep)*1000./(dzt(1)/100.)  !dzt is in cm
              endif

              src(i,k,j,isdfe) = snpzd(7)

              src(i,k,j,isdetrfe) = snpzd(8)
#  endif
!             production of calcite
              dprca = (morp+morz+(graz+graz_Z)*(1.-gamma1))
     &              *capr*redctn*rnbio(k)
              prca = prca + dprca*dzt(k)

!             These are sources and sinks of DIC (i.e. remin - pp)
!             all are based on po4 uptake and remineralization
!             dprca is a correction term
#  if defined O_carbon
              src(i,k,j,isdic) = (snpzd(1)*redctp - dprca)
#  endif
#  if defined O_npzd_alk
              src(i,k,j,isalk) = (-snpzd(1)*redntp*1.e-3 - 2.*dprca)
#  endif
#  if defined O_time_averages && defined O_save_npzd

!-----------------------------------------------------------------------
!             accumulate time averages
!-----------------------------------------------------------------------

              if (timavgperts .and. .not. euler2) then
                ta_rnpp(i,k,jrow) = ta_rnpp(i,k,jrow) + rnpp(i,k,j)
                ta_rgraz(i,k,jrow) = ta_rgraz(i,k,jrow) + rgraz(i,k,j)
                ta_rgraz_Z(i,k,jrow) = ta_rgraz_Z(i,k,jrow) 
     &                               + rgraz_Z(i,k,j)
                ta_rgraz_Det(i,k,jrow) = ta_rgraz_Det(i,k,jrow)
     &                                 + rgraz_Det(i,k,j)
                ta_rmorp(i,k,jrow) = ta_rmorp(i,k,jrow) + rmorp(i,k,j)
                ta_rmorpt(i,k,jrow)= ta_rmorpt(i,k,jrow) + rmorpt(i,k,j)
                ta_rmorz(i,k,jrow) = ta_rmorz(i,k,jrow) + rmorz(i,k,j)
                ta_rexcr(i,k,jrow) = ta_rexcr(i,k,jrow) + rexcr(i,k,j)
#   if defined O_npzd_nitrogen
                ta_rnpp_D(i,k,jrow) = ta_rnpp_D(i,k,jrow)
     &                              + rnpp_D(i,k,j)
                ta_rgraz_D(i,k,jrow) = ta_rgraz_D(i,k,jrow)
     &                               + rgraz_D(i,k,j)
                ta_rmorp_D(i,k,jrow) = ta_rmorp_D(i,k,jrow)
     &                               + rmorp_D(i,k,j)
                ta_rnfix(i,k,jrow) = ta_rnfix(i,k,jrow) + rnfix(i,k,j)
#   endif
#   if defined O_npzd_iron
#   if defined O_npzd_iron_diagnostics
                ta_rfeorgads(i,k,jrow) = ta_rfeorgads(i,k,jrow)
     &                              + rfeorgads(i,k,j)
                ta_rdeffe(i,k,jrow) = ta_rdeffe(i,k,jrow)
     &                              + rdeffe(i,k,j)
                ta_rfeprime(i,k,jrow) = ta_rfeprime(i,k,jrow)
     &                              + rfeprime(i,k,j)
                ta_rfesed(i,k,jrow) = ta_rfesed(i,k,jrow)
     &                              + rfesed(i,k,j)
                ta_rbfe(i,k,jrow) = ta_rbfe(i,k,jrow)
     &                              + rbfe(i,k,j)
                ta_rfecol(i,k,jrow) = ta_rfecol(i,k,jrow)
     &                              + rfecol(i,k,j)
#   endif
# if defined O_npzd_chl
                ta_rchl(i,k,jrow) = ta_rchl(i,k,jrow)
     &                              + rchl(i,k,j)
# if defined O_npzd_nitrogen
                ta_rchl_D(i,k,jrow) = ta_rchl_D(i,k,jrow)
     &                              + rchl_D(i,k,j)
#   endif
#   endif
#   endif
#   if defined O_npzd_extra_diagnostics
                ta_ravej(i,k,jrow) = ta_ravej(i,k,jrow)
     &                              + ravej(i,k,j)
                ta_ravej_D(i,k,jrow) = ta_ravej_D(i,k,jrow)
     &                              + ravej_D(i,k,j)
                ta_rgmax(i,k,jrow) = ta_rgmax(i,k,jrow)
     &                              + rgmax(i,k,j)
                ta_rno3P(i,k,jrow) = ta_rno3P(i,k,jrow)
     &                              + rno3P(i,k,j)
                ta_rpo4P(i,k,jrow) = ta_rpo4P(i,k,jrow)
     &                              + rpo4P(i,k,j)
                ta_rpo4_D(i,k,jrow) = ta_rpo4_D(i,k,jrow)
     &                              + rpo4_D(i,k,j)
#   endif
              endif
#  endif
!             calculate total export to get total import for next layer
              expo = expo*dzt(k)
#  if defined O_npzd_iron
              expofe = expofe*dzt(k)
#  endif
            enddo

#  if defined O_npzd_o2
            kmx = kmt(i,jrow)
            do k=1,kmx
!             limit oxygen consumption below concentrations of
!             5umol/kg as recommended in OCMIP
              fo2 = 0.5*tanh(t(i,k,j,io2,taum1)*1000. - 5.)
!             sink of oxygen
              so2 = src(i,k,j,ispo4)*redotp
              src(i,k,j,iso2) = -so2*(0.5 + fo2)
#   if defined O_npzd_nitrogen
!             add denitrification as source term for NO3
              no3flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin)
!             800 = 0.8*1000 = (elec/mol O2)/(elec/mol NO3)*(mmol/mol)
              deni = 800.*no3flag*so2*(0.5 - fo2)
              src(i,k,j,isno3) = src(i,k,j,isno3) - deni
#    if defined O_save_npzd
              rdeni(i,k,jrow) = deni
#    endif
#   endif
            enddo
#  endif

!-----------------------------------------------------------------------
!           remineralize calcite
!-----------------------------------------------------------------------
            kmx = kmt(i,jrow)
            do k=1,kmx-1
#  if defined O_carbon
              src(i,k,j,isdic) = src(i,k,j,isdic) + prca*rcak(k)
#  endif
#  if defined O_npzd_alk
              src(i,k,j,isalk) = src(i,k,j,isalk) + 2.*prca*rcak(k)
#  endif
            enddo
#  if defined O_sed
            if (addflxo .and. eots)
     &        sbc(i,jrow,ircal) = sbc(i,jrow,ircal) + (prca*rcab(kmx)
     &                          - prca*rcak(kmx))*twodt(1)*dzt(kmx)
#  endif
#  if defined O_carbon
            src(i,kmx,j,isdic) = src(i,kmx,j,isdic) + prca*rcab(kmx)
#  endif
#  if defined O_npzd_alk
            src(i,kmx,j,isalk) = src(i,kmx,j,isalk) + 2.*prca*rcab(kmx)
#  endif
#  if defined O_time_averages && defined O_save_npzd

!-----------------------------------------------------------------------
!           accumulate time averages for full depth variables
!-----------------------------------------------------------------------
            if (timavgperts .and. .not. euler2) then
              kmx = kmt(i,jrow)
              expo = prca
              ta_rprocal(i,jrow) =  ta_rprocal(i,jrow) + prca
              do k=1,kmx
                expo = expo*dztr(k)
                ta_rremi(i,k,jrow) = ta_rremi(i,k,jrow) + rremi(i,k,j)
                ta_rexpo(i,k,jrow) = ta_rexpo(i,k,jrow) + rexpo(i,k,j)
                expo = expo - prca*rcak(k)
                ta_rexpocal(i,k,jrow) = ta_rexpocal(i,k,jrow) + expo
                expo = expo*dzt(k)
#   if defined O_npzd_nitrogen && defined O_npzd_o2
                ta_rdeni(i,k,jrow) = ta_rdeni(i,k,jrow) + rdeni(i,k,j)
#   endif
#   if defined O_npzd_iron
                ta_rremife(i,k,jrow) = ta_rremife(i,k,jrow) 
     &               + rremife(i,k,j)
                ta_rexpofe(i,k,jrow) = ta_rexpofe(i,k,jrow)
     &               + rexpofe(i,k,j)
#   endif
              enddo
            endif
#  endif

          endif

        enddo

      enddo
# endif
# if defined O_carbon_fnpzd

!-----------------------------------------------------------------------
!    option to write and read fixed biological fluxes for carbon and
!    alkalinity. code assumes that restarts line up with the first
!    record. the repeating time period and restarts should probably
!    line up with the start of the year.
!-----------------------------------------------------------------------

      do j=js,je
        jrow = j + joff

!       adjust counters and read fluxes, if first row of slab
        if (jrow .eq. 2) then
          mfnpzd = mfnpzd + 1
          if (mfnpzd .gt. mxfnpzd) mfnpzd = 1
          if (nfnpzd .le. mxfnpzd) nfnpzd = nfnpzd + 1
!         read fluxes, if enough records written
          if (nfnpzd .gt. mxfnpzd) then
            fname = new_file_name ("O_fnpzd.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              ib(:) = 1
              ic(:) = 1
              ic(1) = imtm2
              ic(2) = jmtm2
              ic(3) = km
              ib(4) = mfnpzd
              print*, "Reading fnpzd from: ", trim(fname),
     &          " record:", mfnpzd
              call openfile (fname, iou)
#  if defined O_carbon
              call getvara ('O_fnpzd1', iou, imtm2*jmtm2*km, ib, ic
     &,         tmpijk, c1, c0)
              fnpzd(2:imtm1,2:jmtm1,1:km,1)=tmpijk(1:imtm2,1:jmtm2,1:km)
#  endif
#  if defined O_npzd_alk
              call getvara ('O_fnpzd2', iou, imtm2*jmtm2*km, ib, ic
     &,         tmpijk, c1, c0)
              fnpzd(2:imtm1,2:jmtm1,1:km,2)=tmpijk(1:imtm2,1:jmtm2,1:km)
#  endif
            endif
          endif
        endif

        if (nfnpzd .le. mxfnpzd) then

!         store fluxes, if enough records written
          do i=is,ie
            do k=1,km
#  if defined O_carbon
              fnpzd(i,jrow,k,1) = src(i,k,j,isdic)
#  endif
#  if defined O_npzd_alk
              fnpzd(i,jrow,k,2) = src(i,k,j,isalk)
#  endif
            enddo
          enddo
!         write fluxes, if last row of slab
          if (jrow .eq. jmt-1) then
            fname = new_file_name ("O_fnpzd.nc")
            call openfile (fname, iou)
            print*, "Writing fnpzd to: ", trim(fname),
     &          " record:", mfnpzd
            ib(:) = 1
            ic(:) = 1
            ic(1) = imtm2
            ic(2) = jmtm2
            ic(3) = km
            ib(4) = mfnpzd
            time = relyr
            call putvars ('time', iou, mfnpzd, time, c1, c0)
#  if defined O_carbon
            tmpijk(1:imtm2,1:jmtm2,1:km) = fnpzd(2:imtm1,2:jmtm1,1:km,1)
            call putvara ('O_fnpzd1', iou, imtm2*jmtm2*km, ib, ic
     &,       tmpijk, c1, c0)
#  endif
#  if defined O_npzd_alk
            tmpijk(1:imtm2,1:jmtm2,1:km) = fnpzd(2:imtm1,2:jmtm1,1:km,2)
            call putvara ('O_fnpzd2', iou, imtm2*jmtm2*km, ib, ic
     &,       tmpijk, c1, c0)
#  endif
          endif

        else

!         set fluxes to source terms
          do i=is,ie
            do k=1,km
#  if defined O_carbon
              src(i,k,j,isdic) = fnpzd(i,jrow,k,1)
#  endif
#  if defined O_npzd_alk
              src(i,k,j,isalk) = fnpzd(i,jrow,k,2)
#  endif
            enddo
          enddo

        endif

      enddo
# endif
# if defined O_carbon && defined O_carbon_14

!-----------------------------------------------------------------------
!     set source for c14
!-----------------------------------------------------------------------
      do j=js,je
        jrow = j + joff
        do i=is,ie
          if (kmt(i,jrow) .gt. 0) then
            do k=1,kmt(i,jrow)
#  if defined O_npzd
              src(i,k,j,isc14) = src(i,k,j,isdic)*rstd
     &                         - 3.836e-12*t(i,k,j,ic14,taum1)
#  else
              src(i,k,j,isc14) = - 3.836e-12*t(i,k,j,ic14,taum1)
#  endif
            enddo
          endif
        enddo
      enddo
# endif
# if defined O_sed && !defined O_sed_uncoupled

!-----------------------------------------------------------------------
!     set source terms from sediment model
!-----------------------------------------------------------------------
      do j=js,je
        jrow = j + joff
        do i=is,ie
          if (kmt(i,jrow) .gt. 0) then
            k = kmt(i,jrow)
#  if defined O_carbon && defined O_npzd
            src(i,k,j,isdic) = src(i,k,j,isdic) + sbc(i,j,ibdicfx)
#   if defined O_global_sums
            if (addflxo .and. eots) dtoic = dtoic - sbc(i,j,ibdicfx)
     &        *c2dtts*dtxcel(k)*dxt(i)*dyt(jrow)*cst(jrow)*dzt(k)
#   endif
#  endif
#  if defined O_npzd_alk
            src(i,k,j,isalk) = src(i,k,j,isalk) + sbc(i,j,ibalkfx)
#  endif
          endif
        enddo
      enddo
# endif

!-----------------------------------------------------------------------
!     solve for one tracer at a time
!     n = 1 => temperature
!     n = 2 => salinity
!     n > 2 => other tracers (if applicable)
!-----------------------------------------------------------------------

      do n=1,nt

!-----------------------------------------------------------------------
!       calculate advective tracer flux
!-----------------------------------------------------------------------

        call adv_flux (joff, js, je, is, ie, n)

!-----------------------------------------------------------------------
!       calculate diffusive flux across eastern and northern faces
!       of "T" cells due to various parameterizations for diffusion.
!-----------------------------------------------------------------------

# if defined O_consthmix
#  if !defined O_biharmonic || defined O_bryan_lewis_horizontal

!       diffusive flux on eastern face of "T" cells

        do j=js,je
          do k=1,km
            do i=istrt-1,iend
              diff_fe(i,k,j) =
#   if defined O_bryan_lewis_horizontal
     &                         diff_cet(k)*cstdxur(i,j)*
#   else
     &                         ah_cstdxur(i,j)*
#   endif
     &                         (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
        enddo
#   if defined O_isopycmix

!       diffusive flux on northern face of "T" cells
!       (background for isopycnal mixing)

        do j=js-1,je
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              diff_fn(i,k,j) =
#    if defined O_bryan_lewis_horizontal
     &                         diff_cnt(k)*
#    else
     &                         diff_cnt*
#    endif
     &           csu_dyur(jrow)*(t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
        enddo
#   else

!       diffusive flux on northern face of "T" cells is not calculated.
!       instead, it is incorporated into DIFF_Ty for performance issues.

#   endif
#  else

!-----------------------------------------------------------------------
!       diffusive flux across eastern face of "T" cell
!       diffusive flux across northern face of "T" cell
!-----------------------------------------------------------------------

        m = 2 + n
        do j=js,je
          jrow = j + joff
          ahbi_cstr = diff_cet*cstr(jrow)
          do k=1,km
            do i=istrt-1,iend
              diff_fe(i,k,j) = ahbi_cstr*dxur(i)*
     &                         (del2(i+1,k,j,m)-del2(i,k,j,m))
            enddo
          enddo
        enddo
        do j=js-1,je
          jrow = j + joff
          ahbi_csu_dyur = diff_cnt*csu(jrow)*dyur(jrow)
          do k=1,km
            do i=istrt,iend
              diff_fn(i,k,j) = ahbi_csu_dyur*
     &                        (del2(i,k,j+1,m) - del2(i,k,j,m))
            enddo
          enddo
        enddo
#  endif
# else
#  if defined O_smagnlmix

!       diffusive flux on eastern and northern faces of "T" cells

        do j=js,je
          do k=1,km
            do i=istrt-1,iend
              diff_fe(i,k,j) = diff_cet(i,k,j)*cstdxur(i,j)*
     &                         (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
        enddo
        do j=js-1,je
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              diff_fn(i,k,j) = diff_cnt(i,k,j)*csu_dyur(jrow)*
     &                     (t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1))
            enddo
          enddo
        enddo
#  endif
# endif
!-----------------------------------------------------------------------
!       calculate diffusive flux across bottom face of "T" cells
!-----------------------------------------------------------------------

        do j=js,je
          do k=1,km-1
            do i=istrt,iend
              diff_fb(i,k,j) = diff_cbt(i,k,j)*dzwr(k)*
     &                         (t(i,k,j,n,taum1) - t(i,k+1,j,n,taum1))
            enddo
          enddo
        enddo

# if defined O_isopycmix

!-----------------------------------------------------------------------
!       compute isopycnal diffusive flux through east, north,
!       and bottom faces of T cells.
!-----------------------------------------------------------------------

        call isoflux (joff, js, je, is, ie, n)
# endif

!-----------------------------------------------------------------------
!       set surface and bottom vert b.c. on "T" cells for diffusion
!       and advection. for isopycnal diffusion, set adiabatic boundary
!       conditions.
!       note: the b.c. at adv_fb(i,k=bottom,j) is set by the above code.
!             However, it is not set when k=km so it is set below.
!             adv_fb(i,km,j) is always zero (to within roundoff).
!-----------------------------------------------------------------------

        do j=js,je
          jrow   = j + joff
          do i=istrt,iend
            kb              = kmt(i,jrow)
# if defined O_replacst
            diff_fb(i,0,j)  = c0
# else
            diff_fb(i,0,j)  = stf(i,j,n)
# endif
            diff_fb(i,kb,j) = btf(i,j,n)
            adv_fb(i,0,j)   = adv_vbt(i,0,j)*(t(i,1,j,n,tau) +
     &                                        t(i,1,j,n,tau))
            adv_fb(i,km,j)  = adv_vbt(i,km,j)*t(i,km,j,n,tau)
          enddo
        enddo

# if defined O_source_term || defined O_npzd || defined O_carbon_14

!-----------------------------------------------------------------------
!       set source term for "T" cells
!-----------------------------------------------------------------------

        source(:,:,:) = c0

#  if defined O_npzd || defined O_carbon_14
        if (itrc(n) .ne. 0) then
          do j=js,je
            do k=1,km
              do i=istrt,iend
                source(i,k,j) = src(i,k,j,itrc(n))
              enddo
            enddo
          enddo
        endif
#  endif
#  if defined O_shortwave

!-----------------------------------------------------------------------
!       incorporate short wave penetration into source
!-----------------------------------------------------------------------

        if (n .eq. 1) then
          call swflux0 (joff, js, je, istrt, iend, source)
        endif
#  endif
# endif

!-----------------------------------------------------------------------
!       solve for "tau+1" tracer using statement functions to represent
!       each component of the calculation
!-----------------------------------------------------------------------

!       1st: solve using all components which are treated explicitly

        do j=js,je
          jrow   = j + joff
          do k=1,km
            twodt(k) = c2dtts*dtxcel(k)
            do i=istrt,iend
              t(i,k,j,n,taup1) = t(i,k,j,n,taum1) + twodt(k)*(
     &          DIFF_Tx(i,k,j) + DIFF_Ty(i,k,j,jrow,n) + DIFF_Tz(i,k,j)
     &          - ADV_Tx(i,k,j) -  ADV_Ty(i,k,j,jrow,n) -  ADV_Tz(i,k,j)
# if defined O_isopycmix && defined O_gent_mcwilliams && !defined O_fct && !defined O_quicker
     &          - ADV_Txiso(i,k,j,n) - ADV_Tyiso(i,k,j,jrow,n)
     &          - ADV_Tziso(i,k,j)
# endif
# if defined O_source_term || defined O_npzd || defined O_carbon_14
     &          + source(i,k,j)
# endif
# if defined O_plume
     &          + work(i,k,j)*subflux(i,jrow,n)
# endif
     &          )*tmask(i,k,j)
            enddo
          enddo
        enddo
# if defined O_implicitvmix || defined O_isopycmix || defined O_redi_diffusion

!       2nd: add in portion of vertical diffusion handled implicitly

        call ivdift (joff, js, je, istrt, iend, n, twodt)
# endif

        do j=js,je
          call setbcx (t(1,1,j,n,taup1), imt, km)
        enddo

!-----------------------------------------------------------------------
!       construct diagnostics associated with tracer "n"
!-----------------------------------------------------------------------

        call diagt1 (joff, js, je, istrt, iend, n, twodt)

!-----------------------------------------------------------------------
!       end of tracer component "n" loop
!-----------------------------------------------------------------------

      enddo

# if defined O_replacst

      gamma = zw(1)/twodt(1)
      do j=js,je
        jrow = j + joff
        do i=istrt,iend
          stf(i,j,1) = gamma*(sbc(i,jrow,isst) - t(i,1,j,1,taup1))
          stf(i,j,2) = gamma*(sbc(i,jrow,isss) - t(i,1,j,2,taup1))
          t(i,1,j,1,taup1) = sbc(i,jrow,isst)
          t(i,1,j,2,taup1) = sbc(i,jrow,isss)
        enddo
      enddo
# endif
# if defined O_convect_brine

!-----------------------------------------------------------------------
!     explicit convection: adjust column if gravitationally unstable
!     convect brine rejected under all ice categories
!-----------------------------------------------------------------------

      call convect_brine (joff, js, je, is, ie)
# else
#  if !defined O_implicitvmix || defined O_isopycmix

!-----------------------------------------------------------------------
!     explicit convection: adjust column if gravitationally unstable
!-----------------------------------------------------------------------

#   if defined O_fullconvect
      call convct2 (t(1,1,1,1,taup1), joff, js, je, is, ie, kmt)
      do j=js,je
        do n=1,nt
          call setbcx (t(1,1,j,n,taup1), imt, km)
        enddo
      enddo
#   else
      call convct (t(1,1,1,1,taup1), ncon, joff, js, je, istrt, iend
     &,            kmt)
#   endif
#  endif
# endif
# if defined O_save_convection
      if (timavgperts .and. eots) then
        if (joff .eq. 0) nta_conv = nta_conv + 1
        do j=js,je
          jrow = j + joff
          do i=istrt,iend
            ta_totalk(i,jrow) = ta_totalk(i,jrow) + totalk(i,j)
            ta_vdepth(i,jrow) = ta_vdepth(i,jrow) + vdepth(i,j)
            ta_pe(i,jrow) = ta_pe(i,jrow) + pe(i,j)
          enddo
        enddo
      endif
# endif

!-----------------------------------------------------------------------
!     construct diagnostics after convection
!-----------------------------------------------------------------------

      idiag = 10
      call diagt2 (joff, js, je, istrt, iend, idiag)

# if defined O_fourfil || defined O_firfil

!-----------------------------------------------------------------------
!     filter tracers at high latitudes
!-----------------------------------------------------------------------

      if (istrt .eq. 2 .and. iend .eq. imt-1) then
        call filt (joff, js, je)
      else
        write (stdout,'(a)')
     &  'Error: filtering requires is=2 and ie=imt-1 in tracer'
        stop '=>tracer'
      endif
# endif
      do n=1,nt
        do j=js,je
          call setbcx (t(1,1,j,n,taup1), imt, km)
        enddo
      enddo

!-----------------------------------------------------------------------
!     construct diagnostics after filtering (for total dT/dt)
!-----------------------------------------------------------------------

      idiag = 1
      call diagt2 (joff, js, je, istrt, iend, idiag)

# if !defined O_replacst

!-----------------------------------------------------------------------
!     if needed, construct the Atmos S.B.C.(surface boundary conditions)
!     averaged over this segment
!     eg: SST and possibly SSS
!-----------------------------------------------------------------------

      call asbct (joff, js, je, istrt, iend, isst, itemp)
      call asbct (joff, js, je, istrt, iend, isss, isalt)
# endif
# if defined O_carbon
      call asbct (joff, js, je, istrt, iend, issdic, idic)
#  if defined O_carbon_14
      call asbct (joff, js, je, istrt, iend, issc14, ic14)
#  endif
# endif
# if defined O_npzd_alk
      call asbct (joff, js, je, istrt, iend, issalk, ialk)
# endif
# if defined O_npzd_o2
      call asbct (joff, js, je, istrt, iend, isso2, io2)
# endif
# if defined O_npzd
      call asbct (joff, js, je, istrt, iend, isspo4, ipo4)
#  if !defined O_npzd_no_vflux
      call asbct (joff, js, je, istrt, iend, issphyt, iphyt)
      call asbct (joff, js, je, istrt, iend, isszoop, izoop)
      call asbct (joff, js, je, istrt, iend, issdetr, idetr)
#  endif
#  if defined O_npzd_nitrogen
      call asbct (joff, js, je, istrt, iend, issno3, ino3)
#   if !defined O_npzd_no_vflux
      call asbct (joff, js, je, istrt, iend, issdiaz, idiaz)
#   endif
#  endif
#  if defined O_npzd_iron
      call asbct (joff, js, je, istrt, iend, issdfe, idfe)
#   if !defined O_npzd_no_vflux
      call asbct (joff, js, je, istrt, iend, issdetrfe, idetrfe)
#   endif
#  endif
# endif
# if defined O_cfcs_data_transient
      call asbct (joff, js, je, istrt, iend, isscfc11, icfc11)
      call asbct (joff, js, je, istrt, iend, isscfc12, icfc12)
# endif
# if defined O_sed
!-----------------------------------------------------------------------
!     accumulate bottom tracer values for the sediment model
!-----------------------------------------------------------------------

      if (addflxo .and. eots) then
        if (joff .eq. 0) atsed = atsed + twodt(1)
        do j=js,je
          jrow = j + joff
          do i=is,ie
            if (kmt(i,jrow) .gt. 0) then
              k = kmt(i,jrow)
              sbc(i,jrow,ibtemp) = sbc(i,jrow,ibtemp)
     &                           + t(i,k,j,itemp,taup1)*twodt(1)
              sbc(i,jrow,ibsalt) = sbc(i,jrow,ibsalt)
     &                           + t(i,k,j,isalt,taup1)*twodt(1)
#  if defined O_carbon
              sbc(i,jrow,ibdic) = sbc(i,jrow,ibdic)
     &                          + t(i,k,j,idic,taup1)*twodt(1)
#  endif
#  if defined O_npzd_alk
              sbc(i,jrow,ibalk) = sbc(i,jrow,ibalk)
     &                          + t(i,k,j,ialk,taup1)*twodt(1)
#  endif
#  if defined O_npzd_o2
              sbc(i,jrow,ibo2) = sbc(i,jrow,ibo2)
     &                         + t(i,k,j,io2,taup1)*twodt(1)
#  endif
            endif
          enddo
        enddo
      endif
# endif
# if defined O_carbon && defined O_carbon_14

!-----------------------------------------------------------------------
!     calculate diagnostic delta carbon 14
!-----------------------------------------------------------------------

      if (tsiperts .or. timavgperts) then
        rrstd = 1000./rstd
        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              dc14(i,k,j) = (rrstd*t(i,k,j,ic14,taup1)
     &                      /(t(i,k,j,idic,taup1) + epsln) - 1000.)
     &                      *tmask(i,k,j)
            enddo
          enddo
        enddo
      endif
      if (tsiperts .and. eots) then
        if (js+joff .eq. 2) dc14bar = 0.
        do j=js,je
          jrow = j + joff
          fy = cst(jrow)*dyt(jrow)
          do k=1,km
            fyz = fy*dzt(k)
            do i=istrt,iend
              dc14bar = dc14bar + dc14(i,k,j)*dxt(i)*fyz*tmask(i,k,j)
            enddo
          enddo
        enddo
      endif
      if (timavgperts .and. .not. euler2) then
        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=istrt,iend
              ta_dc14(i,k,jrow) = ta_dc14(i,k,jrow) + dc14(i,k,j)
            enddo
          enddo
        enddo
      endif
# endif

      return
      end

      subroutine diagt1 (joff, js, je, is, ie, n, twodt)

!-----------------------------------------------------------------------
!     construct diagnostics associated with tracer component "n"

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       n     = (1,2) = (u,v) velocity component
!       twodt = (2*dtts,dtts) on (leapfrog,mixing) time steps
!-----------------------------------------------------------------------

      implicit none

      integer i, k, j, ip, kr, jq, n, jp, jrow, js, je, joff, is, ie
      integer mask, m

      real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r
      real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso
      real adv_tziso, diff_tx, diff_ty, diff_tz, dtdx, dtdy, dtdz
      real r2dt, cosdyt, fx, darea, boxar, rtwodt, sumdx, delx
      real sumdxr, dxdy, dxdydz

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "accel.h"
      include "coord.h"
      include "cregin.h"
      include "csbc.h"
# if defined O_tracer_averages
      include "ctavg.h"
# endif
      include "diag.h"
      include "diaga.h"
      include "emode.h"
      include "grdvar.h"
      include "hmixc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "vmixc.h"

# if defined O_meridional_tracer_budget
      include "ctmb.h"
      real stor(imt,km), div(imt,km), sorc(imt,km), dif(imt,km)
      real t1(imt), t2(imt), t3(imt), t4(imt)
# endif
# if defined O_time_step_monitor
      real temp1(imt,km), temp2(imt,km), temp3(imt,km)
# endif
      real twodt(km)

# if defined O_isopycmix
      include "isopyc.h"
# endif
      include "fdift.h"

# if defined O_save_mixing_coeff

!-----------------------------------------------------------------------
!     diagnostic: estimate mixing coefficients on east, north, and
!                 bottom face of T cells from the flux
!-----------------------------------------------------------------------

      if (cmixts .and. n .eq. 1 .and. eots) then
        do j=js,je
          jrow = j + joff
          do k=1,km
            do i=2,imt-1
              dtdx = (t(i+1,k,j,1,taum1)-t(i,k,j,1,taum1))
     &                *tmask(i+1,k,j)*tmask(i,k,j)
     &                *cstr(jrow)*dxur(i) + epsln
              ce(i,k,j,2) = diff_fe(i,k,j)/dtdx
#  if !defined O_consthmix || defined O_biharmonic || defined O_isopycmix
              dtdy = (t(i,k,j+1,1,taum1)-t(i,k,j,1,taum1))
     &                *tmask(i,k,j+1)*tmask(i,k,j)
     &                *dyur(jrow) + epsln
              cn(i,k,j,2) = csur(jrow)*diff_fn(i,k,j)/dtdy
#  else
              cn(i,k,j,2) = csur(jrow)*ah
#  endif
            enddo
          enddo
        enddo
        do j=js,je
          jrow = j + joff
          do k=1,km-1
            do i=2,imt-1
              dtdz = (t(i,k,j,1,taum1)-t(i,k+1,j,1,taum1))
     &                *tmask(i,k,j)*tmask(i,k+1,j)
     &                *dzwr(k) + epsln
              cb(i,k,j,2) = diff_fb(i,k,j)/dtdz
#  if defined O_isopycmix
     &                      + diff_fbiso(i,k,j)/dtdz
#  endif
            enddo
          enddo
          do i=2,imt-1
            cb(i,km,j,2) = 0.0
          enddo
        enddo

        do j=js,je
          call setbcx (ce(1,1,j,2), imt, km)
          call setbcx (cn(1,1,j,2), imt, km)
          call setbcx (cb(1,1,j,2), imt, km)
        enddo
      endif
# endif

# if defined O_save_convection_full

!-----------------------------------------------------------------------
!     diagnostic: initialize temperature before convection
!-----------------------------------------------------------------------

      if (exconvts .and. n .eq. 1 .and. eots) then
        do j=js,je
          do k=1,km
            do i=1,imt
              excnv0(i,k,j) = t(i,k,j,1,taup1)
            enddo
          enddo
        enddo
      endif

# endif

# if defined O_time_step_monitor

!-----------------------------------------------------------------------
!     diagnostic: integrate |d(tracer)/dt|  and tracer variance on "tau"
!                 globally
!-----------------------------------------------------------------------

      if (tsiperts .and. eots) then
        do j=js,je
          jrow = j + joff
          r2dt    = c1/c2dtts
          cosdyt  = cst(jrow)*dyt(jrow)
          do k=1,km
            fx = r2dt/dtxcel(k)
            do i=is,ie
              darea      = dzt(k)*dxt(i)*cosdyt*tmask(i,k,j)
              temp3(i,k) = t(i,k,j,n,tau)*darea
              temp1(i,k) = t(i,k,j,n,tau)**2*darea
              temp2(i,k) = abs(t(i,k,j,n,taup1)-t(i,k,j,n,taum1))*
     &                     darea*fx
            enddo
            do i=is,ie
              tbar(k,n,jrow)   = tbar(k,n,jrow) + temp3(i,k)
              travar(k,n,jrow) = travar(k,n,jrow) + temp1(i,k)
              dtabs(k,n,jrow)  = dtabs(k,n,jrow) + temp2(i,k)
            enddo
          enddo
        enddo
      endif
# endif

# if defined O_tracer_averages

!-----------------------------------------------------------------------
!     diagnostic: accumulate tracers for averages under horizontal
!                 regions (use units of meters, rather than cm)
!-----------------------------------------------------------------------

      if (tavgts .and. eots) then
        do j=js,je
          jrow = j + joff
          do i=is,ie
            mask = mskhr(i,jrow)
            if (mask .ne. 0) then
              boxar = cst(jrow)*dxt(i)*dyt(jrow)*tmask(i,1,j)*0.0001
              sumbf(mask,n) = sumbf(mask,n) + stf(i,j,n)*boxar
              do k=1,km
                sumbk(mask,k,n) = sumbk(mask,k,n) + t(i,k,j,n,tau)
     &                             *boxar*dzt(k)*tmask(i,k,j)*0.01
             enddo
            endif
          enddo
        enddo
      endif
# endif

# if defined O_tracer_yz

!----------------------------------------------------------------------
!     diagnostic: integrate tracer equations in longitude
!----------------------------------------------------------------------

      if (tyzts .and. eots) then
        do j=js,je
          jrow = j + joff
          do k=1,km
            rtwodt = c1/twodt(k)
            sumdx = c0
            do m=1,5
              tyz(jrow,k,n,m) = c0
            enddo
            do i=is,ie
              delx  = dxt(i)*tmask(i,k,j)
              sumdx = sumdx + delx
              tyz(jrow,k,n,1) = tyz(jrow,k,n,1)+ delx*t(i,k,j,n,tau)
              tyz(jrow,k,n,2) = tyz(jrow,k,n,2) + rtwodt*delx*
     &          (t(i,k,j,n,taup1) - t(i,k,j,n,taum1))
              tyz(jrow,k,n,3) =  tyz(jrow,k,n,3)
     &          - delx*(ADV_Ty(i,k,j,jrow,n) + ADV_Tz(i,k,j))
              tyz(jrow,k,n,4) =  tyz(jrow,k,n,4)
     &          + delx*(DIFF_Ty(i,k,j,jrow,n) + DIFF_Tz(i,k,j))
#  if defined O_source_term || defined O_npzd || defined O_carbon_14
              tyz(jrow,k,n,5) = tyz(jrow,k,n,5) + delx*source(i,k,j)
#  endif
            enddo
            sumdxr = c1/(sumdx + epsln)
            do m=1,5
              tyz(jrow,k,n,m) = tyz(jrow,k,n,m)*sumdxr
            enddo
          enddo
        enddo
        if (js+joff .eq. 2) then
          do m=1,5
            do k=1,km
              tyz(1,k,n,m)   = c0
              tyz(jmt,k,n,m) = c0
            enddo
          enddo
        endif
      endif
# endif

# if defined O_meridional_tracer_budget

!----------------------------------------------------------------------
!     diagnostic: integrate equations in depth, lon, and time for
!                 each basin and jrow. basin #  0 is used to catch
!                 values in land areas
!----------------------------------------------------------------------

      if (tmbperts .and. eots) then
        do j=js,je
          jrow = j + joff
          if (jrow .eq. 2 .and. n .eq. 1) numtmb = numtmb + 1
          cosdyt = cst(jrow)*dyt(jrow)
          do k=1,km
            do i=is,ie
              stor(i,k) = 0.0
              div(i,k)  = 0.0
              sorc(i,k) = 0.0
              dif(i,k)  = 0.0
            enddo
          enddo
          do k=1,km
            rtwodt = c1/twodt(k)
            do i=is,ie
              dxdy            = cosdyt*dxt(i)*tmask(i,k,j)
              dxdydz          = dxdy*dzt(k)
              stor(i,k) = stor(i,k) + rtwodt*dxdydz*
     &                            (t(i,k,j,n,taup1) - t(i,k,j,n,taum1))
              div(i,k)  = div(i,k)  - dxdydz*ADV_Ty(i,k,j,jrow,n)
#  if defined O_source_term || defined O_npzd || defined O_carbon_14
              sorc(i,k) = sorc(i,k) + dxdydz*source(i,k,j)
#  endif
              dif(i,k)  = dif(i,k)  + dxdydz*DIFF_Ty(i,k,j,jrow,n)
            enddo
          enddo

!         the accumulation is done this way for issues of speed

          do i=is,ie
            t1(i) = stor(i,1)
            t2(i) = div(i,1)
            t3(i) = sorc(i,1)
            t4(i) = dif(i,1)
          enddo

          do k=2,km
            do i=is,ie
              t1(i) = t1(i) + stor(i,k)
              t2(i) = t2(i)  + div(i,k)
              t3(i) = t3(i) + sorc(i,k)
              t4(i) = t4(i)  + dif(i,k)
            enddo
          enddo

          do i=is,ie
            m               = msktmb(i,jrow)
            tstor(jrow,n,m) = tstor(jrow,n,m) + t1(i)
            tdiv(jrow,n,m)  = tdiv(jrow,n,m)  + t2(i)
            tsorc(jrow,n,m) = tsorc(jrow,n,m) + t3(i)
            tdif(jrow,n,m)  = tdif(jrow,n,m)  + t4(i)
          enddo

          k = 1
          do i=is,ie
            m               = msktmb(i,jrow)
            dxdy            = cosdyt*dxt(i)*tmask(i,k,j)
            tflux(jrow,n,m) = tflux(jrow,n,m) + dxdy*stf(i,j,n)
          enddo
          if (n .eq. 1) then
            do k=1,km
              do i=is,ie
                m               = msktmb(i,jrow)
                dxdy            = cosdyt*dxt(i)*tmask(i,k,j)
                dxdydz          = dxdy*dzt(k)
                smdvol(jrow,m)  = smdvol(jrow,m)  + dxdydz
              enddo
            enddo
          endif
        enddo
      endif
# endif

# if defined O_gyre_components

!-----------------------------------------------------------------------
!     diagnostic: compute the northward transport components of
!                 each tracer
!-----------------------------------------------------------------------

      if (gyrets .and. eots)  call gyre (joff, js, je, is, ie, n)
# endif

# if defined O_term_balances

!-----------------------------------------------------------------------
!     diagnostic: integrate r.h.s. terms in the tracer equations
!                 over specified regional volumes.
!-----------------------------------------------------------------------

      if (trmbts .and. eots)  call ttb1 (joff, js, je, is, ie, n)
# endif

# if defined O_xbts

!-----------------------------------------------------------------------
!     diagnostic: accumulate r.h.s. terms in the tracer equations
!-----------------------------------------------------------------------

      if (xbtperts .and. eots) call txbt1 (joff, js, je, n)
# endif
# if defined O_mom_tbt

!-----------------------------------------------------------------------
!     diagnostic: accumulate r.h.s. terms in the tracer equations
!-----------------------------------------------------------------------

      if (tbtperts .and. eots) call tbt1 (joff, js, je, n)
# endif

      return
      end

      subroutine diagt2 (joff, js, je, is, ie, idiag)

!-----------------------------------------------------------------------
!     construct d(tracer)/dt diagnostics

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       idiag = 1  => total tracer change
!       idiag = 10 => change of tracer due to filtering(also convection)
!-----------------------------------------------------------------------

      implicit none

      integer idiag, j, js, je, k, i, joff, iocv, jrow, is, ie

      real rdt, reltim, period

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "diaga.h"
      include "iounit.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"
      include "tmngr.h"
      include "timeavgs.h"

# if defined O_save_convection_full

!-----------------------------------------------------------------------
!     diagnostic: save tracer time change due to explicit convection
!     idiag = 10 signifies diagnostics immediately after convection
!-----------------------------------------------------------------------

      if (exconvts .and. idiag .eq. 10 .and. eots) then
        rdt = c1/c2dtts

!       excnv1 = epsln over land cells and d(convect)/dt over ocean

        do j=js,je
          do k=1,km
            do i=1,imt
              excnv1(i,k,j) = tmask(i,k,j)*
     &         (t(i,k,j,1,taup1)-excnv0(i,k,j))*rdt
     &                       + (c1-tmask(i,k,j))*epsln
            enddo
          enddo
        enddo

        reltim = relyr
        if (joff + js .eq. 2) then
          write (stdout,*) ' =>Writing explicit convection at ts=',itt
     &  , ' ',stamp
          call getunit (iocv, 'cvct.dta'
     &,                'unformatted sequential append ieee')

          period = 0.0
          iotext = 'read(iocv) reltim, period, imt, jmt, km, flag'
          write (iocv) stamp, iotext, expnam
          write (iocv) reltim, period, imt, jmt, km, epsln

          iotext = 'read(iocv) (xt(i),i=1,imt)'
          write (iocv) stamp, iotext, expnam
          call wrufio (iocv, xt, imt)

          iotext = 'read(iocv) (yt(j),j=1,jmt)'
          write (iocv) stamp, iotext, expnam
          call wrufio (iocv, yt, jmt)

          iotext = 'read(iocv) (zt(k),k=1,km)'
          write (iocv) stamp, iotext, expnam
          call wrufio (iocv, zt, km)

          call relunit (iocv)
        endif

        call getunit (iocv, 'cvct.dta'
     &,                'unformatted sequential append ieee')

        do j=js,je
          jrow = j+joff

          write(iotext,'(a10,i4)') ' for jrow=',jrow
          iotext(15:) = ': read (iocv) (convect(i,k),i=1,imt),k=1,km)'
          write (iocv) stamp, iotext, expnam
          call wrufio (iocv, excnv1(1,1,j), imt*km)

        enddo
        call relunit(iocv)
      endif
# endif

# if defined O_term_balances

!-----------------------------------------------------------------------
!     diagnostic: integrate d/dt(tracer) over specified regional volumes
!                  after convection and filtering
!-----------------------------------------------------------------------

      if (trmbts .and. eots) call ttb2 (joff, js, je, is, ie, idiag)
# endif

# if defined O_xbts

!-----------------------------------------------------------------------
!     diagnostic: accumulate d/dt(tracer) after convection
!                 and filtering
!-----------------------------------------------------------------------

      if (xbtperts .and. eots) call txbt2 (joff, js, je, idiag)
# endif
# if defined O_mom_tbt

!-----------------------------------------------------------------------
!     diagnostic: accumulate d/dt(tracer) after convection
!                 and filtering

!-----------------------------------------------------------------------

      if (tbtperts .and. eots) call tbt2 (joff, js, je, idiag)
# endif

      return
      end

      subroutine asbct (joff, js, je, is, ie, isbc, itr)

!-----------------------------------------------------------------------
!     construct the Atmos S.B.C. (surface boundary conditions)

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       isbc  = index for sbc
!       itr   = index for tracer
!-----------------------------------------------------------------------

      implicit none

      integer isbc, itr, j, js, je, jrow, joff, i, is, ie

      real rts

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "switch.h"

!     initialize the Atmos S.B.C. at the start of each ocean segment
!     (do not alter values in land)

      if (isbc .le. 0 .or. itr .le. 0) return

      if (eots .and. osegs) then
        do j=js,je
          jrow   = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0) sbc(i,jrow,isbc) = c0
          enddo
        enddo
      endif

!     accumulate surface tracers for the Atmos S.B.C. every time step

      if (eots) then
        do j=js,je
          jrow = j + joff
          do i=is,ie
            sbc(i,jrow,isbc) = sbc(i,jrow,isbc)+t(i,1,j,itr,taup1)
          enddo
        enddo
      endif

!     average the surface tracers for the Atmos S.B.C. at the end of
!     each ocean segment. (do not alter values in land)

      if (eots .and. osege) then
        rts = c1/ntspos
        do j=js,je
          jrow   = j + joff
          do i=is,ie
            if (kmt(i,jrow) .ne. 0)
     &        sbc(i,jrow,isbc) = rts*sbc(i,jrow,isbc)
          enddo
        enddo
      endif

      return
      end

      subroutine ivdift (joff, js, je, is, ie, n, twodt)

# if defined O_implicitvmix || defined O_isopycmix || defined O_redi_diffusion
!-----------------------------------------------------------------------
!     solve vertical diffusion of tracers implicitly

!     input:
!       joff  = offset relating "j" in the MW to latitude "jrow"
!       js    = starting row in the MW
!       je    = ending row in the MW
!       is    = starting longitude index in the MW
!       ie    = ending longitude index in the MW
!       n     = tracer component
!       twodt = (2*dtts, dtts) on (leapfrog, mixing) time steps
!-----------------------------------------------------------------------

      implicit none

      integer j, js, je, k, i, is, ie, n, joff

      real rc2dt

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "levind.h"
      include "mw.h"
      include "switch.h"
      include "vmixc.h"

      real twodt(km)

!     store terms to compute implicit vertical mixing on
!     diagnostic time steps

#  if defined O_xbts || defined O_mom_tbt
      if (eots) then
        do j=js,je
          do k=1,km
            do i=is,ie
              zzi(i,k,j) = t(i,k,j,n,taup1)
            enddo
          enddo
        enddo
      endif
#  else
#   if defined O_term_balances
      if (trmbts .and. eots) then
        do j=js,je
          do k=1,km
            do i=is,ie
              zzi(i,k,j) = t(i,k,j,n,taup1)
            enddo
          enddo
        enddo
      endif
#   endif
#  endif

      call invtri (t(1,1,1,n,taup1), stf(1,1,n), btf(1,1,n)
     &, diff_cbt(1,1,jsmw), twodt, kmt, tmask(1,1,1), is, ie
     &, joff, js, je)

!     compute residual implicit vertical mixing

#  if defined O_xbts || defined O_mom_tbt
      if (eots) then
        do j=js,je
          do k=1,km
            rc2dt = c1/twodt(k)
            do i=is,ie
              zzi(i,k,j) = rc2dt*(t(i,k,j,n,taup1) - zzi(i,k,j))
            enddo
          enddo
        enddo
      endif
#  else
#   if defined O_term_balances
      if (trmbts .and. eots) then
        do j=js,je
          do k=1,km
            rc2dt = c1/twodt(k)
            do i=is,ie
              zzi(i,k,j) = rc2dt*(t(i,k,j,n,taup1) - zzi(i,k,j))
            enddo
          enddo
        enddo
      endif
#   endif
#  endif
# endif

      return
      end

      subroutine swflux0 (joff, js, je, is, ie, source)

# if defined O_mom && defined O_shortwave
!=======================================================================
!       incorporate short wave penetration via the "source"
!       term. note that the divergence of shortwave for the first
!       level "divpen(1)" is compensating for the effect of having
!       the shortwave component already included in the total
!       surface tracer flux "stf(i,j,1)"

!       incorporating the penetrative shortwave radiative flux into
!       the vertical diffusive flux term directly is not correct when
!       vertical diffusion is solved implicitly.

!     input:
!       joff = offset relating "j" in the MW to latitude "jrow"
!       js   = starting row in the MW
!       je   = ending row in the MW
!       is   = starting longitude index in the MW
!       ie   = ending longitude index in the MW

!     output:
!       source = source term penetrative short wave
!=======================================================================

      implicit none

      integer j, js, je, jrow, joff, k, i, is, ie

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "cshort.h"

      real source(imt,km,jsmw:jemw)

      if (iswr .ne. 0) then
        do j=js,je
          jrow   = j + joff
          do k=1,km
            do i=is,ie
              source(i,k,j) = source(i,k,j)
     &                       + sbc(i,jrow,ipsw)*divpen(k)
            enddo
          enddo
        enddo
      endif
# endif
#endif

      return
      end
!
# if defined O_npzd_fe_limitation
      subroutine polint (xa,ya,n,x,y,dy)
!
      implicit none
!
      integer n,i,m,ns, nmax
      parameter (nmax = 10) ! largest anticipated value of n
      real dy,x,y,xa(n),ya(n)
      real den,dif,dift,ho,hp,w,c(nmax),d(nmax)
!
! Given arrays xa and ya, each of length n, and a given value x, this routine
! returns a value y, and an error estimate dy. If P(x) is the polynomial of
! degree N-1 such that P(xai) = yai, i=1,....,n then the returned value y=P(x)
!
!
      ns=1
      dif=abs(x-xa(1))
!
      do i=1,n
         dift = abs(x-xa(i))
         IF (dift.lt.dif) THEN
            ns = i
            dif = dift
         END IF
!
       c(i) = ya(i)
       d(i) = ya(i)
!
      enddo
!
      y = ya (ns)        ! this is the initial approximation to y
      ns = ns-1
!
      do m=1,n-1
         do i=1,n-m
            ho = xa(i) - x
            hp = xa (i+m) -x
            w = c(i+1) - d(i)
            den = ho - hp
            IF (den.eq.0.) PAUSE 'failure in polint'
            den = w/den
            d(i) = hp*den
            c(i) = ho*den
         enddo
         IF (2*ns.lt.n-m) THEN
            dy = c(ns+1)
         ELSE
            dy = d(ns)
            ns = ns - 1
         END IF
         y = y + dy
      enddo
      return
      end
!
      subroutine hunt (xa,n,x,jlo)
!
      integer jlo,n,inc,jhi,jm
      real x, xa(n)
      logical ascnd
!
! Given an array xa(1:n), and given a value x, returns a value jlo
! such that x is between xa(jlo) and xa(jlo +1). xa(1:n) must be
! monotonic, either increasing or decreasing. jlo=0 or jlo=n is
! returned to indicate that x is out of range. jlo on input is taken
! as the initial guess for jlo on output
!
      ascnd = xa(n).ge.xa(1)    ! True if ascending order of table, false otherwise
      IF (jlo.le.0.or.jlo.gt.n) THEN ! Input guess not useful. Go to bisection
         jlo = 0
         jhi = n+1
         goto 3
      END IF
      inc = 1                   ! Set the hunting increment
      IF (x.ge.xa(jlo).eqv.ascnd) THEN ! Hunting up:
 1       jhi = jlo + inc
         IF (jhi.gt.n) THEN     ! Done hunting, since off end of table
            jhi = n + 1
         ELSE IF (x.ge.xa(jhi).eqv.ascnd) THEN ! Not done hunting
            jlo = jhi
            inc = inc + inc     ! so double the increment
      goto 1              ! and try again
      END IF                    ! Done hunting, value bracketed.
      ELSE                      ! Hunt down:
         jhi = jlo
 2       jlo = jhi-inc
         IF (jlo.lt.1) THEN     !Done hunting, since off end of table
            jlo = 0
         ELSE IF (x.lt.xa(jlo).eqv.ascnd) THEN ! Not done hunting
            jhi = jlo
            inc =  inc + inc    ! so double the increment
            goto 2              ! and try again
         END IF                 ! Done hunting, value bracketed
      END IF                    ! Hunt is done, so begin the final bisection phase:
!
 3    IF (jhi-jlo.eq.1) THEN
         IF (x.eq.xa(n)) jlo = n-1
         IF (x.eq.xa(1)) jlo = 1
         RETURN
      END IF
      jm = (jhi + jlo)/2
      IF (x.ge.xa(jm).eqv.ascnd) THEN
         jlo = jm
      ELSE
         jhi = jm
      END IF
!
      GOTO 3
!
      end
#endif
