      subroutine setmom (is, ie, js, je)

#if defined O_mom
!=======================================================================
!     set up everything which must be done only once per run
!=======================================================================

      implicit none

      character (120) :: fname, new_file_name

      integer ie, is, je, js, i, ioun, j, jrow, n, m, indp, ip, k, kr
      integer jq, mask, kz, nv, ke, ks, iou, isle

      integer ib(10), ic(10)

      logical error, vmixset, hmixset, exists

      real snapint, snapls, snaple, snapde, trajint, cksum, checksum
      real ocnp, boxat, boxau, dvolt, dvolu, sum, pctocn, diffa, amix
      real runstep, dtatms, ahbkg, spnep, senep, c1e3

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "accel.h"
      include "calendar.h"
      logical annlev, annlevobc, initpt
# if defined O_neptune
      include "cnep.h"
# endif
      include "coord.h"
# if defined O_fourfil || defined O_firfil
      include "cpolar.h"
# endif
      include "cprnts.h"
      include "cregin.h"
      include "csbc.h"
# if defined O_shortwave
      include "cshort.h"
# endif
      include "ctmb.h"
      include "diag.h"
      include "docnam.h"
      include "emode.h"
      include "grdvar.h"
      include "hmixc.h"
      include "index.h"
      include "iounit.h"
      include "isleperim.h"
# if defined O_isopycmix
      include "isopyc.h"
# endif
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "stab.h"
      include "state.h"
      include "switch.h"
      include "tmngr.h"
      include "vmixc.h"
# if defined O_mom_tbt
      include "tbt.h"
# endif
# if defined O_fourfil || defined O_firfil
      integer kmz(imt,jmt)
# endif
# if defined O_tidal_kv
      include "tidal_kv.h"
# endif
# if defined O_fwa
      include "fwa.h"
# endif
      real tmpij(imtm2,jmtm2)

# if defined O_tidal_kv
      gravrho0r = grav*rho0r
      zetar = 1./500.e2        ! inverse vertical decay scale, zeta = 500m
      ogamma = 0.2*rho0r*zetar ! Osborn, 1980's mixing efficiency, gamma
                               ! scaled by 1./(zeta*rho0)
# endif
!      error in tracer conservation generated by wide_open_mw = .true.
!      if (jmw .eq. jmt) then
!        wide_open_mw = .true.
!      else
        wide_open_mw = .false.
!      endif

!     stability diagnostic

      call stabi

      error = .false.
      vmixset = .false.
      hmixset = .false.

      visc_cbu_limit = 1.0e6
      diff_cbt_limit = 1.0e6
# if defined O_rigid_lid_surface_pressure

      alph = 1.0
      gam = 0.0
      theta = 1.0
# endif
# if defined O_implicit_free_surface

      alph = 0.3333333
      gam = 0.3333333
      theta = 0.5
# endif

!     user specified tracer names are placed into "trname" here.

      do m=1,nt
        trname(m) = '**unknown***'
      enddo
      trname(1) = 'potentl_temp'
      trname(2) = 'salinity    '

!-----------------------------------------------------------------------
!     get i/o units needed for prognostic variables
!-----------------------------------------------------------------------

      call getunitnumber (kflds)
      call getunitnumber (latdisk(1))
      call getunitnumber (latdisk(2))
# if defined O_tidal_kv

!-----------------------------------------------------------------------
!     read in 2d array of tidal energy dissipation rate
!     similar to that of Jayne and St. Laurent, GRL 2001
!     except here we use E0 = W/(A*Nlev)

!     no need to convert MOM native cgs
!     units lesson: W = kg m^2 /s^3, therefore, W/m^2 = kg/s^3
!     therefore multiply edr by N0*1e3 to get g/s^3, where N0 is some
!     reference value of N.
!     in tidal.nc the energy flux already includes the Nb, so it's only
!     multiplied by 1e3 to convert W/m^2 in g/s^3
!-----------------------------------------------------------------------

      edrm2(:,:,:)=0.
      edrs2(:,:,:)=0.
      edrk1(:,:,:)=0.
      edro1(:,:,:)=0.
c      edr(:,:) = 0.
      ib(:) = 1
      ic(:) = 1
      ic(1) = imt
      ic(2) = jmt
      ic(3) = 1
c      fname = new_file_name ("O_tidenrg_egbert.nc")
c      fname = new_file_name ("O_tidenrg_green.nc")
      fname = new_file_name ("O_tidenrg.nc")
      inquire (file=trim(fname), exist=exists)
      if (exists) then
        c1e3 = 1.e3
        call openfile (trim(fname), iou)
        do k=1,km
           ib(3) = k
c           print*, "read tidal_egbert, k=",k
           call getvara ('M2', iou, imt*jmt, ib, ic 
     &          ,edrm2(1:imt,k,1:jmt), c1e3, c0)
           call getvara ('S2', iou, imt*jmt, ib, ic 
     &          ,edrs2(1:imt,k,1:jmt), c1e3, c0)
           call getvara ('K1', iou, imt*jmt, ib, ic 
     &          ,edrk1(1:imt,k,1:jmt), c1e3, c0)
           call getvara ('O1', iou, imt*jmt, ib, ic 
     &          ,edro1(1:imt,k,1:jmt), c1e3, c0)
c           print*, edrm2(1:imt,k,1:jmt)
        enddo
        call closefile (iou)
c        call getvara ('tidenrg', iou, imtm2*jmtm2, ib, ic, tmpij
c     &,   c1e3, c0)
c        edr(2:imtm1,2:jmtm1) = tmpij(1:imtm2,1:jmtm2)
      else
        print*, "Warning => Can not find ",trim(fname)
      endif
# endif

!-----------------------------------------------------------------------
!     compute density coefficients based on depth of grid points
!-----------------------------------------------------------------------

      call eqstate (zt, km, ro0, to, so, c, tmink, tmaxk, smink, smaxk)

      cksum = checksum(c, km, 9)
      write (stdout
     &,'(6x,"Checksum for density coefficients =",e14.7)') cksum

# if defined O_linearized_density
!     eliminate the pressure effect when using linearized density

      do k=1,km
        ro0(k) = c0
        to(k)  = c0
        so(k)  = c0
      enddo

# endif
# if defined O_linearized_advection
!     set the initial idealized stratification as a function of
!     temperature only. Salinity is fixed at 35 ppt (tbar(k,2)=0) but
!     passive tracers (if nt>2) are 1.0 at k=1 and 0.0 for k>1

      tzero = 7.5
      tone  = 10.0
      z0    = 30.0e2
      bigl  = 80.0e2
      write (stdout,'(//a/)')
     & 'Initial Tracer profile as a function of depth'
      do k=1,km
        tbarz(k,1) = tzero*(1.0-tanh((zt(k)-bigl)/z0)) +
     &                tone*(1.0-zt(k)/zt(km))
        tbarz(k,2) = c0
        if (nt .gt. 2) then
          do n=3,nt
            if (k .eq. 1) then
              tbarz(k,n) = c1
            else
              tbarz(k,n) = c0
            endif
          enddo
        endif
        write (stdout,'(1x,"k=",i3,", zt(k)=",f8.1
     &,        "cm,  tbarz(k,n),n=1,nt =",10e14.7)')
     &      k, zt(k),(tbarz(k,n),n=1,nt)
      enddo

# endif
!-----------------------------------------------------------------------
!     if the MW is not fully opened, then set time level indicators in
!     the MW ("tau-1" "tau" "tau+1") to constant values.
!-----------------------------------------------------------------------

      if (.not. wide_open_mw) then
        taum1 = -1
        tau   =  0
        taup1 = +1
      endif

!-----------------------------------------------------------------------
!     set prognostic quantities to either initial conditions or restart
!-----------------------------------------------------------------------
# if defined O_neptune

!     calculate "neptune" velocities

      call neptune
# endif

!     initialize two dimensional fields on disk

# if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface

!     initialize surface pressure fields in memory

      do jrow=1,jmt
        do i=1,imt
          ps(i,jrow,1)     = c0
          ps(i,jrow,2)     = c0
          pguess(i,jrow)   = c0
          ubar(i,jrow,1)   = c0
          ubar(i,jrow,2)   = c0
          ubarm1(i,jrow,1) = c0
          ubarm1(i,jrow,2) = c0
        enddo
      enddo
# endif
# if defined O_stream_function

!     initialize stream function fields in memory

      do n=1,2
        do jrow=1,jmt
          do i=1,imt
#  if defined O_neptune

!           initialize to "neptune" streamfunction

            psi(i,jrow,n) = pnep(i,jrow)
            ptd(i,jrow)   = c0
#  else
            psi(i,jrow,n) = c0
            ptd(i,jrow)   = c0
#  endif
          enddo
        enddo
      enddo

!     initialize stream function guess fields on disk
!     block`s 1 & 2 are for the stream function guess field on disk

      do n=1,nkflds
        call oput (kflds, nwds, n, ptd(1,1))
      enddo
# endif

!     initialize all latitude rows

      call rowi

!     initialize time step counter = 0

      itt    = 0
      irstdy = 0
      msrsdy = 0

!     initialize all prognostic quantities from the restart

      if (.not. init) then
        fname = new_file_name ("restart_mom.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) call mom_rest_in (fname,  is, ie, js, je)
# if defined O_restart_2
        fname = new_file_name ("restart_2_mom.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) call mom_rest_in (fname, is, ie, js, je)
# endif
      endif

!     construct depth arrays associated with "u" cells
      call depth_u (kmt, imt, jmt, zw, km, kmu, h, hr)

!     compute a topography checksum

      cksum = 0.0
      ocnp   = 0
      do jrow=2,jmt-1
        do i=2,imt-1
          cksum = cksum + i*jrow*kmt(i,jrow)
          ocnp   = ocnp + float(kmt(i,jrow))
        enddo
      enddo
      write (stdout,*) ' "kmt" checksum = ', cksum

# if defined O_time_averages
!-----------------------------------------------------------------------
!     zero integrated time average accumulators
!-----------------------------------------------------------------------

      call ta_mom_tsi (0)

# endif
!-----------------------------------------------------------------------
!     initialize the time manager with specified initial conditions
!     time, user reference time, model time, and how long to integrate.
!-----------------------------------------------------------------------

      call tmngri (year0, month0, day0, hour0, min0, sec0
     &,            ryear, rmonth, rday, rhour, rmin, rsec
     &,            irstdy, msrsdy, runlen, rununits, rundays, dtts)

# if defined O_stability_tests
!-----------------------------------------------------------------------
!     convert starting and ending longitudes for the stability tests
!     to nearest model grid points.
!-----------------------------------------------------------------------

      if (stabint .ge. c0) then
        iscfl  = max(indp (cflons, xt, imt), 2)
        cflons = xt(iscfl)
        iecfl  = min(indp (cflone, xt, imt), imt-1)
        cflone = xt(iecfl)
        jscfl  = max(indp (cflats, yt, jmt), 2)
        cflats = yt(jscfl)
        jecfl  = min(indp (cflate, yt, jmt), jmt-1)
        cflate = yt(jecfl)
        kscfl  = indp (cfldps, zt, km)
        cfldps = zt(kscfl)
        kecfl  = indp (cfldpe, zt, km)
        cfldpe = zt(kecfl)
      endif
# endif
# if defined O_restorst && !defined O_replacst

!-----------------------------------------------------------------------
!      damp surface tracers back to data. schematically, the restoring
!      term will be = (dampdz/(dampts*daylen)) * (data - t)
!      where dampdz is the thickness (cm) and dampts is the time
!      scale (days)
!-----------------------------------------------------------------------

      do n=1,nt
        write (stdout,'(/,1x,a,i2,a,a,/,1x,1pe14.7,a,1pe14.7,a,/)')
     &  ' Surface tracer #',n,' will be damped back to data using a'
     &, ' Newtonian restoring time scale of '
     &,  dampts(n),' days. and a level thickness =', dampdz(n),' cm.'
      enddo
# endif
# if defined O_shortwave

!-----------------------------------------------------------------------
!     Solar Shortwave energy penetrates below the ocean surface. Clear
!     water assumes energy partitions between two exponentials as
!     follows:

!     58% of the energy decays with a 35 cm e-folding scale
!     42% of the energy decays with a 23 m e-folding scale

!     if the thickness of the first ocean level "dzt(1)" is 50 meters,
!     then shortwave penetration wouldn't matter. however, for
!     dzt(1) = 10 meters, the effect can be significant. this can be
!     particularly noticeable in the summer hemisphere.

!     Paulson and Simpson (1977 Irradiance measurements in the upper
!                               ocean JPO 7, 952-956)
!     Also see ... Jerlov (1968 Optical oceanography. Elsevier)
!                  A General Circulation Model for Upper Ocean
!                  Simulation (Rosati and Miyakoda JPO vol 18,Nov 1988)
!-----------------------------------------------------------------------

      write (stdout,'(/a/)')
     &' => Shortwave penetration is a double exponential as follows:'
      rpart  = 0.58
      efold1 = 35.0e0
      efold2 = 23.0e2
      rscl1  = 1.0/efold1
      rscl2  = 1.0/efold2

!     note that pen(0) is set 0.0 instead of 1.0 to compensate for the
!     shortwave part of the total surface flux in "stf(i,1)"

      pen(0) = c0

      do k=1,km
        swarg1 = -min(zw(k)*rscl1,70.0)
        swarg2 = -min(zw(k)*rscl2,70.0)
        pen(k) = rpart*exp(swarg1) + (c1-rpart)*exp(swarg2)
        divpen(k) = (pen(k-1) - pen(k))/dzt(k)
        write (stdout,9200) k, zw(k)*1.e-2, pen(k), zt(k)*1.e-2
     &,                     divpen(k)
      enddo
      write (stdout,*) ' '
9200  format (1x,'k=',i3,' zw(k)=',f8.2,'(m)  pen(k)=',e14.7
     &,       ' zt(k)=',f8.2,'(m)  divpen(k)=',e14.7)
# endif

!-----------------------------------------------------------------------
!     compute the surface area and volume of the ocean regions. index 0
!     refers to the sum of all regions.
!     (values used in subroutine region are done in terms of meters,
!     rather than centimetres)
!-----------------------------------------------------------------------

      areag = c0
      volgt  = c0

      do k=1,km
        volgk(k) = c0
      enddo

      do n=0,numreg
        areat(n) = c0
        areau(n) = c0
        volt(n)  = c0
        volu(n)  = c0
      enddo
      do mask=1,nhreg
        areab(mask) = c0
        volbt(mask) = c0
        do k=1,km
          volbk(mask,k) = c0
        enddo
      enddo

      do jrow=2,jmtm1
        do i=2,imtm1
          mask = mskhr(i,jrow)
          kz   = kmt(i,jrow)
          if (kz .gt. 0 .and. mask .gt. 0) then
            boxat = cst(jrow) * dxt(i) * dyt(jrow)
            if (kmu(i,jrow) .ne. 0) then
              boxau = csu(jrow) * dxu(i) * dyu(jrow)
            else
              boxau = c0
            endif
            areat(0) = areat(0) +  boxat
            areau(0) = areau(0) +  boxau
            areab(mask)  = areab(mask) + boxat * 1.e-4
            do k=1,kz
              volbk(mask,k) = volbk(mask,k) + boxat * dzt(k) * 1.e-6
              dvolt   = boxat * dzt(k)
              if (kmu(i,jrow) .ge. k) then
                dvolu   = boxau * dzt(k)
              else
                dvolu   = c0
              endif
              n = nhreg*(mskvr(k)-1) + mask
              if (n .gt. 0) then
                volt(0) = volt(0) +  dvolt
                volu(0) = volu(0) +  dvolu
                volt(n) = volt(n) +  dvolt
                volu(n) = volu(n) +  dvolu
                do nv=1,nvreg
                  ks = llvreg(nv,1)
                  if (k .eq. ks) then
                    areat(n) = areat(n) +  boxat
                    areau(n) = areau(n) +  boxau
                  endif
                enddo
              endif
            enddo
          endif
        enddo
      enddo

      do mask=0,numreg
        if (areat(mask) .eq. c0) then
          rareat(mask) = c0
        else
          rareat(mask) = c1/areat(mask)
        endif

        if (areau(mask) .eq. c0) then
          rareau(mask) = c0
        else
          rareau(mask) = c1/areau(mask)
        endif

        if (volt(mask) .eq. c0) then
          rvolt(mask) = c0
        else
          rvolt(mask) = c1/volt(mask)
        endif

        if (volu(mask) .eq. c0) then
          rvolu(mask) = c0
        else
          rvolu(mask) = c1/volu(mask)
        endif
      enddo

      do mask=1,nhreg
        do k=1,km
          volbt(mask) = volbt(mask) + volbk(mask,k)
          volgk(k) = volgk(k) + volbk(mask,k)
        enddo
        areag = areag + areab(mask)
        volgt = volgt + volbt(mask)
      enddo

# if defined O_tracer_averages
      if (iotavg .ne. stdout .or. iotavg .lt. 0) then
        call getunit (iou, 'tracer_avg.dta'
     &,               'unformatted sequential append ieee')
        call reg1st (iou, .false., .true., .true., .false., .false.)
        call relunit (iou)
      endif
# endif
      if (iotavg .eq. stdout .or. iotavg .lt. 0) then
        call reg1st (stdout, .true., .true., .true., .false., .false.)
      endif

!     compute and print statistics for regions

      sum = c0
      do n=1,numreg
        sum = sum + volt(n)
      enddo
      sum    = 100.0*sum/tcellv
      pctocn = 100.0*ocnp/float((imt-2)*(jmt-2)*km)
      diffa  = 100.0 * (c1 - (tcella(1) - 10000.0*areag)/tcella(1))

      write (stdout,9342) diffa, numreg, sum, pctocn
9342  format ('  the horizontal regional masks cover',f8.3
     &, '% of the total ocean surface area.'/
     &, '  there are ', i6, ' regions over which tracer & '
     &, 'momentum balances will be computed',/,'  accounting for '
     &, f6.2,'% of the total ocean volume.'/
     &, 1x,f6.2,'% of the grid points lie within the ocean.'/)
#
!-----------------------------------------------------------------------
!     find all land mass perimeters for Poisson solvers
!-----------------------------------------------------------------------

      auto_kmt_changes = .false.
!     set mask for land mass perimeters on which to perform calculations
!     imask(-n) = .false.  [no equations ever on dry land mass n]
!     imask(0)  = .true.   [equations at all mid ocean points]
!     imask(n)  = .true./.false [controls whether there will be
!                                equations on the ocean perimeter of
!                                land mass n]
!     note: land mass 1 is the northwest-most land mass
!           usually includes the "north pole", and at low resolutions,
!           the "main continent"

      do isle=-mnisle,mnisle
        if (isle .ge. 0 .and. isle .le. nisle) then
          imask(isle) = .true.
        else
          imask(isle) = .false.
        endif
      enddo
# if defined O_symmetry

!     do not perform island integrals for land masses whose perimeters
!     touch the equator in models symmetric about the equator

      do isle=1,nisle
        do n=1,nippts(isle)
          j = jperm(iofs(isle)+n)
          if (j .eq. jmt-1) then
            imask(isle) = .false.
          endif
        enddo
      enddo
# endif

!     user-specified changes to island mask
!       imask(1) = .true.
!       imask(2) = .true.

!     there are problems if imask is set .true. for a nonexistent
!     island.

!     print diagnostic information

      do isle=-mnisle,mnisle
        if (imask(isle)) then
          if (isle .eq. 0) then
            print '(a)','=> calculations enabled for mid ocean points'
          else
            print '(2a,i3)','=> calculations enabled for ocean ',
     &                      'perimeter of land mass',isle
          endif
        endif
      enddo
      do isle=0,nisle
        if (.not. imask(isle)) then
            print '(2a,i3)','=> calculations disabled for ocean ',
     &                      'perimeter of land mass',isle
        endif
      enddo

!     imain is the land mass on which dpsi is normalized to 0
!     if imain is 0, then dpsi is not normalized.
!     default value of imain is land mass with longest perimeter

      imain = min(2,nisle)
      do isle=1,nisle
        if (nippts(isle) .gt. nippts(imain)) then
          imain = isle
        endif
      enddo

!     if any island perimeter is not calculated, imain must be one such

      do isle=1,nisle
        if (.not.imask(isle)) then
          imain = isle
        endif
      enddo
# if defined O_symmetry

!     do not normalize dpsi when using symmetry about the equator

      imain = 0
# endif

      if (imain .gt. 0 .and. imain .le. nisle) then
        print '(a,i4)', 'dpsi normalized to zero on land mass',imain
      elseif (imain .eq. 0) then
        print *, 'no normalization on dpsi'
      else
        print *, 'ERROR: illegal value for choice of normalization ',
     &           'land mass, imain =', imain
      endif
      print *,' (user may set "imain" to any valid land mass number)'

!---------------------------------------------------------------------
!     compute checksum of density coefficients
!---------------------------------------------------------------------

      print *,' '
      call print_checksum (c(1,1), km, 9
     &,                   ' density coefficient checksum = ')

# if defined O_stream_function

!-----------------------------------------------------------------------
!     checksum the starting stream function.
!-----------------------------------------------------------------------

      call print_checksum (psi(1,1,1), imt, jmt
     &, ' checksum for psi(,,1) = ')
      call print_checksum (psi(1,1,2), imt, jmt
     &, ' checksum for psi(,,2) = ')
# endif

# if defined O_fourfil || defined O_firfil

!-----------------------------------------------------------------------
!     compute an array to indicate "interior" stream function grid cells
!-----------------------------------------------------------------------

      do jrow=1,jmt
        kmz(1,jrow)   = 0
        kmz(imt,jrow) = 0
      enddo

      do i=1,imt
        kmz(i,1)   = 0
        kmz(i,jmt) = 0
      enddo

      do jrow=2,jmtm1
        do i=2,imt
          kmz(i,jrow) = min(kmu(i-1,jrow-1), kmu(i,jrow-1)
     &,                     kmu(i-1,jrow), kmu(i,jrow))
        enddo
      enddo
#  if defined O_cyclic
      do jrow=1,jmt
        kmz(1,jrow) = kmz(imtm1,jrow)
      enddo
#  endif

!---------------------------------------------------------------------
!     find and print start & end indices for filtering
!---------------------------------------------------------------------

      write (stdout,9551)
      if (lsegf.gt.11) write (stdout,9552)
      write (stdout,9553)
      call findex (kmt, jmtfil, km, jft1, jft2, imt, istf, ietf)
      write (stdout,9554)
      call findex (kmu, jmtfil, km, jfu1, jfu2, imt, isuf, ieuf)
      write (stdout,9555)
#  if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface
      call findex (kmu, jmtfil, 1, jfu1, jfu2, imt, iszf, iezf)
#  endif
#  if defined O_stream_function
      call findex (kmz, jmtfil, 1, jft1, jft2, imt, iszf, iezf)
#  endif
9551  format (/' ==== start and end indices for Fourier filtering ====')
9552  format (' only 11 sets of indices fit across the page.',
     &       '  others will not be printed.'/)
9553  format (/,' == filtering indices for t,s ==')
9554  format (/,' == filtering indices for u,v ==')
9555  format (/,' == filtering indices for stream function ==')
# endif

!---------------------------------------------------------------------
!     print the timestep multipliers
!---------------------------------------------------------------------

      write (stdout,9601) (dtxcel(k),k=1,km)
9601  format(/,' "dtxcel(km)" tracer timestep multipliers:',/,10(f9.3))

!-----------------------------------------------------------------------
!     initialize various things
!-----------------------------------------------------------------------

      do jrow=1,jmt
        do i=1,imt
# if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface
          uhat(i,jrow,1) = c0
          uhat(i,jrow,2) = c0
          divf(i,jrow) = c0
# endif
# if defined O_stream_function
          ztd(i,jrow) = c0
# endif
          zu(i,jrow,1)  = c0
          zu(i,jrow,2)  = c0
        enddo
      enddo

!     Coriolis factors

      do jrow=1,jmt
        do i=1,imt
          cori(i,jrow,1) = c2*omega*sin(ulat(i,jrow)/radian)
          cori(i,jrow,2) = -cori(i,jrow,1)
        enddo
      enddo

!     metric diffusion factors

# if defined O_consthmix
#  if defined O_biharmonic
      amix = sqrt(abs(ambi))
#  else
      amix = am
#  endif
      do jrow=1,jmt
        am3(jrow)   = amix*(c1-tng(jrow)*tng(jrow))/(radius**2)
        am4(jrow,1) = -amix*c2*sine(jrow)/(radius*csu(jrow)
     &                                     *csu(jrow))
        am4(jrow,2) = -am4(jrow,1)
      enddo
# endif

!     metric advection factors

      do jrow=1,jmt
        advmet(jrow,1) = tng(jrow)/radius
        advmet(jrow,2) = -advmet(jrow,1)
      enddo

!     diffusive flux through bottom of cells

      do j=jsmw,jemw
        do k=0,km
          do i=1,imt
            diff_fb(i,k,j) = c0
            adv_fb(i,k,j)  = c0
          enddo
        enddo
      enddo

!-----------------------------------------------------------------------
!     initialize diagnostics
!-----------------------------------------------------------------------

# if defined O_meridional_tracer_budget

!-----------------------------------------------------------------------
!     set basins and initialize arrays
!-----------------------------------------------------------------------

      if (tmbint .ge. c0) then

!       set all points to signify basin # 1

        do jrow=1,jmt
          do i=2,imtm1
            msktmb(i,jrow) = 1
          enddo
          msktmb(1,jrow)   = 0
          msktmb(imt,jrow) = 0
        enddo

!       divide "msktmb" into other basins (2 .. ntmbb) here if desired.

!       set all land points to basin # 0

        do jrow=1,jmt
          do i=1,imt
            if (kmt(i,jrow) .eq. 0) msktmb(i,jrow) = 0
          enddo
        enddo

        if (ntmbb .gt. 1) then
          write (stdout,*)
     &  ' Basin arrangement for Meridional Tracer Budget diagnostic'
          call iplot (msktmb, imt, imt, jmt)
        endif

!       write out the meridional tracer budget basin mask

        if ((iotmb .ne. stdout .or. iotmb .lt. 0) .and. itmb) then
          call getunit (iu, 'tracer_bud.dta'
     &,                'unformatted sequential append ieee')
          iotext =
     &  ' read (iotmb) imt, jmt, ((msktmb(i,j),i=1,imt),j=1,jmt)'
          write (iu) stamp, iotext, expnam
          write (iu) imt, jmt, ((msktmb(i,j),i=1,imt),j=1,jmt)
          call relunit (iu)
        endif
      endif

      numtmb = 0
      do mask=0,ntmbb
        do jrow=1,jmt
          smdvol(jrow,mask)  = c0
        enddo
      enddo
      do mask=0,ntmbb
        do n=1,nt
          do jrow=1,jmt
            tstor(jrow,n,mask) = c0
            tdiv(jrow,n,mask)  = c0
            tflux(jrow,n,mask) = c0
            tdif(jrow,n,mask)  = c0
            tsorc(jrow,n,mask) = c0
          enddo
        enddo
      enddo
# endif

# if defined O_bryan_lewis_vertical || defined O_bryan_lewis_horizontal

!-----------------------------------------------------------------------
!     initialize Bryan_Lewis tracer diffusion coefficients
!-----------------------------------------------------------------------

      call blmixi
# endif
# if defined O_time_averages

!-----------------------------------------------------------------------
!     initialize time mean "averaging" grid data
!-----------------------------------------------------------------------

      call avgi

# endif
# if defined O_xbts

!-----------------------------------------------------------------------
!     initialize XBT locations and averaging arrays
!-----------------------------------------------------------------------

      call xbti

# endif
# if defined O_mom_tbt

!-----------------------------------------------------------------------
!     initialize tbt averaging arrays
!-----------------------------------------------------------------------

      do j=1,jmt
        do i=1,imt
          do n=1,nt
            do k=1,kmt(i,j)
              do m=1,11
                tbt(i,j,k,n,m) = 0.
              enddo
            enddo
            tbtsf(i,j,n) = 0.
          enddo
        enddo
      enddo
      ntbtts = 0

# endif
# if defined O_ppvmix

!-----------------------------------------------------------------------
!     initialize pacanowski-philander vertical mixing scheme
!-----------------------------------------------------------------------

      call ppmixi (vmixset)
# endif
# if defined O_smagnlmix && !defined O_consthmix

!-----------------------------------------------------------------------
!     initialize Smagorinsky nonlinear horizontal mixing scheme
!-----------------------------------------------------------------------

      call smagnli (hmixset)
# endif
# if defined O_isopycmix

!-----------------------------------------------------------------------
!    initialize the isopycnal mixing
!-----------------------------------------------------------------------

      call isopi (error, am, ah)
# endif
# if defined O_npzd
!-----------------------------------------------------------------------
!    initialize the model of ocean biogeochemistry and isotopes (MOBI)
!-----------------------------------------------------------------------

      if (kpzd .ne. km) then
         print*, 'stop: kpzd is not equal to km'
         print*, '      currently inconsistent with use of mobi'
      endif
      do k=1,km
         if (dtxcel(k) .ne. 1.0) then
            print*, 'stop: dtxcel is not equal to one'
            print*, '      inconsistent with use of mobi'
         endif
      enddo

      call mobii 

# endif
# if defined O_fwa

!-----------------------------------------------------------------------
!     calculate areas and masks for fresh water anomalies
!-----------------------------------------------------------------------

      fname = new_file_name ("O_fwawt.nc")
      inquire (file=trim(fname), exist=exists)
      if (.not. exists) then
        print*, "Warning => ", trim(fname), " does not exist."
        fwawt(:,:) = 0.
        if (mrfwa .gt. 0 .and. mrfwa .le. nhreg) then
          where (mskhr(:,:) .eq. mrfwa) fwawt(:,:) = 1.
          print*, "          using regional mask to define area"
        else
          print*, "          using box corners to define area"
          fwawt(isfwa1:iefwa1,jsfwa:jefwa) = 1.
          fwawt(isfwa2:iefwa2,jsfwa:jefwa) = 1.
        endif
      else
        call openfile (fname, iou)
        ib(:) = 1
        ic(:) = 1
        ic(1) = imtm2
        ic(2) = jmtm2
        call getvara ('O_fwawt', iou, imtm2*jmtm2, ib, ic, tmpij
     &,   c1, c0)
        fwawt(2:imtm1,2:jmtm1) = tmpij(1:imtm2,1:jmtm2)
        call embmbc (fwawt)
      endif

      areafwc = 0.
      areafwa = 0.
      do j=2,jmtm1
        do i=2,imtm1
          if (kmt(i,j) .gt. 0) then
            areafwa = areafwa + dxt(i)*dyt(j)*cst(j)*fwawt(i,j)
            fwcwt(i,j) = 1. - fwawt(i,j)
            areafwc = areafwc + dxt(i)*dyt(j)*cst(j)*fwcwt(i,j)
          else
            fwawt(i,j) = 0.
            fwcwt(i,j) = 0.
          endif
        enddo
      enddo

# endif
!-----------------------------------------------------------------------
!     do all consistency checks last
!-----------------------------------------------------------------------

      call checks (error, vmixset, hmixset)

      return
      end

      subroutine depth_u (kmt, imt, jmt, zw, km, kmu, h, hr)

!=======================================================================
!     calculate depth arrays associated with "u" cells.

!     input:
!       kmt = number of oecan "t" cells from surface to bottom of ocean
!       imt = longitudinal dimension of arrays
!       jmt = latitudinal dimension of arrays
!       zw  = depth to bottom of "t" cells
!       km  = max number of depths

!     output:
!       kmu = number of ocean "u" cells from surface to bottom of ocean
!       h   = depth to ocean floor over "u" cells (cm)
!       hr  = reciprocal of "h"
!=======================================================================

      implicit none

      integer i, imt, jmt, jrow, km
      integer kmt(imt,jmt), kmu(imt,jmt)

      real c0, c1
      real  h(imt,jmt), hr(imt,jmt), zw(km)

!-----------------------------------------------------------------------
!     set some constants
!-----------------------------------------------------------------------

      c0 = 0.0
      c1 = 1.0

!-----------------------------------------------------------------------
!     compute number of vertical levels on the "u" grid
!-----------------------------------------------------------------------

      do jrow=1,jmt
        kmu(imt,jrow) = 0
      enddo

      do i=1,imt
        kmu(i,jmt) = 0
      enddo

      do jrow=1,jmt-1
        do i=1,imt-1
           kmu(i,jrow) = min (kmt(i,jrow), kmt(i+1,jrow), kmt(i,jrow+1)
     &,                       kmt(i+1,jrow+1))
        enddo
      enddo
# if defined O_cyclic
      do jrow=1,jmt
        kmu(imt,jrow) = kmu(2,jrow)
      enddo
# endif
# if defined O_symmetry
      do i=1,imt
        kmu(i,jmt) = kmu(i,jmt-2)
      enddo
# endif

!---------------------------------------------------------------------
!     compute depths and reciprocal depths over "u" cells
!---------------------------------------------------------------------

      do jrow=1,jmt
        do i=1,imt
          hr(i,jrow) = c0
          h(i,jrow)  = c0
          if (kmu(i,jrow) .ne. 0) then
            hr(i,jrow) = c1/zw(kmu(i,jrow))
            h (i,jrow) = zw(kmu(i,jrow))
          endif
        enddo
      enddo

      return
      end

      subroutine rowi

!-----------------------------------------------------------------------
!     initialize prognostic quantities at "tau-1" and "tau"

!     inputs:

!     kmt  = number of vertical levels on "t" cells
!     yt   = latitudes of "t" points
!     zt   = depths of "t" points
!-----------------------------------------------------------------------

      implicit none

      character(120) :: fname, new_file_name, text

      integer i, itt, iou, ib(10), ic(10), j, jrow, n, k, kz, ln

      logical exists

      real dens, drodt, drods, p1, c100, c1e3, p001, p035, ucksum
      real vcksum, tcksum, scksum, theta0, tq, sq, s, d, tins, tinsit
      real checksum, C2K, jdi
# if defined O_gthflx
      real geoheatflux
# endif

      include "size.h"
      include "param.h"
      include "npzd.h" ! rstd*
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "csbc.h"
      include "iounit.h"
      include "levind.h"
      include "mw.h"
# if defined O_tai_slh
      include "diag.h"
      include "state.h"
      include "dens.h"
# endif

      real tmpik(imtm2,km), dmsk(imt,jmt)

!-----------------------------------------------------------------------
!     update pointers to tau-1, tau, & tau+1 data on disk.
!     for latitude rows they point to latdisk(1) or latdisk(2)
!     for 2D fields they point to records on kflds
!-----------------------------------------------------------------------

      itt   = 0
      taum1disk = mod(itt+1,2) + 1
      taudisk   = mod(itt  ,2) + 1
      taup1disk = taum1disk

      p1 = 0.1
      c100 = 100.
      c1e3 = 1000.
      p001 = 0.001
      p035 = 0.035
      C2K = 273.15

!-----------------------------------------------------------------------
!     update pointers to tau-1, tau, & tau+1 data in the MW
!-----------------------------------------------------------------------

      if (wide_open_mw) then

!       rotate time levels instead of moving data

        taum1 = mod(itt+0,3) - 1
        tau   = mod(itt+1,3) - 1
        taup1 = mod(itt+2,3) - 1
      endif

      ucksum = 0.0
      vcksum = 0.0
      tcksum = 0.0
      scksum = 0.0

!-----------------------------------------------------------------------
!     initialize every latitude jrow either in the MW (when wide opened)
!     or on disk (when jmw < jmt)
!-----------------------------------------------------------------------

      do jrow=1,jmt

        if (wide_open_mw) then
          j = jrow
        else
          j = jmw
        endif
        jdi = min(max(jrow-1,1),jmtm2)

!-----------------------------------------------------------------------
!       zero out all variables. velocities are internal modes only
!-----------------------------------------------------------------------

        do k=1,km
          do i=1,imt
            u(i,k,j,1,taup1) = c0
            u(i,k,j,2,taup1) = c0
          enddo
        enddo
        do n=1,nt
          do k=1,km
            do i=1,imt
              t(i,k,j,n,taup1) = c0
            enddo
          enddo
        enddo

!-----------------------------------------------------------------------
!       set tracers
!-----------------------------------------------------------------------

        do n=1,nt

          do i=1,imt
            do k=1,kmt(i,jrow)
              if (trim(mapt(n)) .eq. 'temp') then
                t(i,k,j,n,taup1) = theta0 (yt(jrow), zt(k))
              elseif (trim(mapt(n)) .eq. 'salt') then
                t(i,k,j,n,taup1) = 0.03472 - 0.035
              elseif (trim(mapt(n)) .eq. 'dic') then
                t(i,k,j,n,taup1) = 2.315
              elseif (trim(mapt(n)) .eq. 'alk') then
                t(i,k,j,n,taup1) = 2.429
              elseif (trim(mapt(n)) .eq. 'o2') then
                t(i,k,j,n,taup1) = 0.1692
              elseif (trim(mapt(n)) .eq. 'po4') then
                if (k .eq. 1) then
                  t(i,k,j,n,taup1) = 0.543
                else
                  t(i,k,j,n,taup1) = 2.165
                endif
              elseif (trim(mapt(n)) .eq. 'dop') then
                if (k .le. 2) then
                  t(i,k,j,n,taup1) = 0.1
                elseif (k .gt. 2 .and. k .le. 12) then
                   t(i,k,j,n,taup1) = 0.025
                else
                  t(i,k,j,n,taup1) = 0.005
                endif
              elseif (trim(mapt(n)) .eq. 'phyt') then
                t(i,k,j,n,taup1) = 0.14*exp((zt(1)-zt(k))/100.e2)
              elseif (trim(mapt(n)) .eq. 'zoop') then
                t(i,k,j,n,taup1) = 0.014*exp((zt(1)-zt(k))/100.e2)
              elseif (trim(mapt(n)) .eq. 'detr') then
                t(i,k,j,n,taup1) = 1.e-4
              elseif (trim(mapt(n)) .eq. 'no3') then
                if (k .eq. 1) then
                  t(i,k,j,n,taup1) = 5.30
                else
                  t(i,k,j,n,taup1) = 30.84
                endif
              elseif (trim(mapt(n)) .eq. 'don') then
                if (k .le. 2) then
                  t(i,k,j,n,taup1) = 3.5
                elseif (k .gt. 2 .and. k .le. 12) then
                   t(i,k,j,n,taup1) = 1.5
                else
                  t(i,k,j,n,taup1) = 0.5
                endif
              elseif (trim(mapt(n)) .eq. 'diaz') then
                t(i,k,j,n,taup1) = 0.014*exp((zt(1)-zt(k))/100.e2)
              elseif (trim(mapt(n)) .eq. 'c14') then
                t(i,k,j,n,taup1) = -150 ! permil (converted later)
              elseif (trim(mapt(n)) .eq. 'cfc11') then
                t(i,k,j,n,taup1) = 0.0
              elseif (trim(mapt(n)) .eq. 'cfc12') then
                t(i,k,j,n,taup1) = 0.0
              else
                if (k .eq. 1) then
                  t(i,k,j,n,taup1) = c1
                else
                  t(i,k,j,n,taup1) = c0
                endif
              endif
            enddo
          enddo

# if !defined O_idealized_ic
          ib(:) = 1
          ib(2) = jdi
          ib(3) = 1
          ic(:) = imtm2
          ic(2) = 1
          ic(3) = km
          if (trim(mapt(n)) .eq. 'temp') then
!           expecting units of degrees K
            fname = new_file_name ("O_temp.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
              call getvara ('O_temp', iou, imtm2*km, ib, ic, tmpik
     &,         c1, c0)
              text = "C"
              call getatttext (iou, 'O_temp', 'units', text)
!             convert to model units (C)
              if (trim(text).eq."K")
     &          where ((tmpik(1:imtm2,1:km)) .lt. 1.e30)
     &            tmpik(1:imtm2,1:km) = tmpik(1:imtm2,1:km) - C2K
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'salt') then
!           expecting units of psu
            fname = new_file_name ("O_sal.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
!             convert to model units ((psu-35)/1000)
              call getvara ('O_sal', iou, imtm2*km, ib, ic, tmpik
     &,         p001, -p035)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'dic') then
!           expecting units of mol m-3 (equals model units, umol cm-3)
            fname = new_file_name ("O_totcarb.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
              call getvara ('O_totcarb', iou, imtm2*km, ib, ic, tmpik
     &,         c1, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'alk') then
!           expecting units of mol m-3 (equals model units, umol cm-3)
            fname = new_file_name ("O_alk.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
              call getvara ('O_alk', iou, imtm2*km, ib, ic, tmpik
     &,         c1, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'o2') then
!           expecting units of mol m-3 (equals model units, umol cm-3)
            fname = new_file_name ("O_o2.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
              call getvara ('O_o2', iou, imtm2*km, ib, ic, tmpik
     &,         c1, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'po4') then
!           expecting units of mol m-3
            fname = new_file_name ("O_po4.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
!             convert to model units (nmol cm-3)
              call getvara ('O_po4', iou, imtm2*km, ib, ic, tmpik
     &,         c1e3, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'no3') then
!           expecting units of mol m-3
            fname = new_file_name ("O_no3.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
!             convert to model units (nmol cm-3)
              call getvara ('O_no3', iou, imtm2*km, ib, ic, tmpik
     &,         c1e3, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          elseif (trim(mapt(n)) .eq. 'c14') then
!           expecting units of permil
            fname = new_file_name ("O_dc14.nc")
            inquire (file=trim(fname), exist=exists)
            if (exists) then
              call openfile (trim(fname), iou)
              tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1)
!             convert later
              call getvara ('O_dc14', iou, imtm2*km, ib, ic, tmpik
     &,         c1, c0)
              t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km)
              call setbcx (t(1:imt,1:km,j,n,taup1),imt,km)
            else
              if (jrow .eq. 1)
     &          print*, "Warning => Can not find ",trim(fname)
            endif
          endif

          do k=1,km
            do i=1,imt
              if (kmt(i,jrow) .lt. k)  t(i,k,j,n,taup1) = c0
            enddo
          enddo
# endif

        enddo
# if defined O_npzd_nitrogen_15
        do k=1,km
          do i=1,imt
             t(i,k,j,idin15,taup1) = 1.005*rn15std*t(i,k,j,ino3,taup1)
     &           / (1 + 1.005*rn15std)
             t(i,k,j,idon15,taup1) = 1.0045*rn15std*t(i,k,j,idon,taup1)
     &           / (1 + 1.0045*rn15std)
              t(i,k,j,idon15,taup1) = rn15std*t(i,k,j,idon,taup1)
     &           / (1 + rn15std)
             t(i,k,j,iphytn15,taup1) = rn15std*t(i,k,j,iphyt,taup1)
     &           / (1 + rn15std)
             t(i,k,j,izoopn15,taup1) = rn15std*t(i,k,j,izoop,taup1)
     &           / (1 + rn15std)
             t(i,k,j,idetrn15,taup1) = rn15std*t(i,k,j,idetr,taup1)
     &           / (1 + rn15std)
             t(i,k,j,idiazn15,taup1) = rn15std*t(i,k,j,idiaz,taup1)
     &           / (1 + rn15std)
          enddo
       enddo
# endif
# if defined O_carbon_13
        do k=1,km
           do i=1,imt
!     initialize at delta del13C = 0 permil
             t(i,k,j,idic13,taup1) = t(i,k,j,idic,taup1)*1.0004*rc13std
     &             / (1. + 1.0004*rc13std)
             t(i,k,j,idoc13,taup1) = t(i,k,j,idon,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
             t(i,k,j,iphytc13,taup1) = t(i,k,j,iphyt,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
             t(i,k,j,izoopc13,taup1) = t(i,k,j,izoop,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
             t(i,k,j,idetrc13,taup1) = t(i,k,j,idetr,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
#  if defined O_npzd_nitrogen
             t(i,k,j,idiazc13,taup1) = t(i,k,j,idiaz,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
#  endif
           enddo
        enddo
# endif
# if defined O_carbon && defined O_carbon_14

!       convert c14 from permil to model units (umol cm-3)
        do k=1,km
          do i=1,imt
            t(i,k,j,ic14,taup1) = (t(i,k,j,ic14,taup1)*0.001 + 1)
     &                            *t(i,k,j,idic,taup1)*rc14std
          enddo
        enddo
# endif
# if defined O_linearized_advection

!       initialize density profile with tbarz

        do n=1,nt
          do i=1,imt
            do k=1,km
              t(i,k,j,n,taup1) = tbarz(k,n)
            enddo
          enddo
        enddo
        if (jrow .eq. jmt) then
          write (stdout,'(a,a)')
     & '=> Note: initialized T & S with idealized profile T(z), S=35ppt'
        endif
# endif

# if defined O_tai_slh
!-----------------------------------------------------------------------
!       set sea level reference density field
!-----------------------------------------------------------------------

!       set reference density field to initial field
        do i=1,imt
          do k=1,km
            t_slh(i,jrow,k,1) = t(i,k,j,itemp,taup1)
            t_slh(i,jrow,k,2) = t(i,k,j,isalt,taup1)
          enddo
        enddo
!       read reference density field if it exists
        fname = new_file_name ("O_slhref.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) then
          call openfile (fname, iou)
          ib(:) = 1
          ib(2) = jdi
          ic(:) = 1
          ic(1) = imtm2
          ic(3) = km
          tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,itemp,taup1)
          call getvara ('O_temp', iou, imtm2*km, ib, ic, tmpik
     &,     c1, c0)
          text = "C"
          call getatttext (iou, 'O_temp', 'units', text)
!         convert to model units (C)
          if (trim(text) .eq. "K")
     &      where ((tmpik(1:imtm2,1:km)) .lt. 1.e30)
     &        tmpik(1:imtm2,1:km) = tmpik(1:imtm2,1:km) - C2K
          do i=2,imtm1
            do k=1,km
              t_slh(i,jrow,k,1) = tmpik(i-1,k)
            enddo
          enddo
          tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,isalt,taup1)
          call getvara ('O_sal', iou, imtm2*km, ib, ic, tmpik
     &,     p001, -p035)
          do i=2,imtm1
            do k=1,km
              t_slh(i,jrow,k,2) = tmpik(i-1,k)
            enddo
          enddo
        endif

!-----------------------------------------------------------------------
!       calculate reference density field
!-----------------------------------------------------------------------

        do i=2,imtm1
          do k=1,kmt(i,jrow)
            s = (t_slh(i,jrow,k,2)*1000.) + 35.
            d = zt(k)*0.01
            tins = tinsit(t_slh(i,jrow,k,1), s, d)
            d_slh(i,jrow,k) = dens(tins-to(k),t_slh(i,jrow,k,2)-so(k),k)
          enddo
        enddo
!       zero for future accumulation
        t_slh(:,:,:,:) = 0.

# endif
# if defined O_sed

!-----------------------------------------------------------------------
!       initialize bottom boundary for missing ocean tracers from data
!-----------------------------------------------------------------------
        ib(:) = 1
        ib(2) = jdi
        ib(3) = 1
        ic(:) = imtm2
        ic(2) = 1
        ic(3) = km
#  if !defined O_carbon
        fname = new_file_name ("O_totcarb.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) then
          call openfile (trim(fname), iou)
          call getvara ('O_totcarb', iou, imtm2*km, ib, ic, tmpik
     &,     c1, c0)
          do i=2,imtm1
            k = kmt(i,jrow)
            sbc(i,jrow,ibdic) = tmpik(i-1,k)
          enddo
          call setbcx (sbc(:,jrow,ibdic), imt, 1)
        else
          print*, "Error => Can not find ",trim(fname)
          stop '=> setmom sed dic'
        endif
#  endif
#  if !defined O_npzd_alk
        fname = new_file_name ("O_alk.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) then
          call openfile (trim(fname), iou)
          call getvara ('O_alk', iou, imtm2*km, ib, ic, tmpik, c1, c0)
          do i=2,imtm1
            k = kmt(i,jrow)
            sbc(i,jrow,ibalk) = tmpik(i,k)
          enddo
        call setbcx (sbc(:,jrow,ibalk), imt, 1)
        else
          print*, "Error => Can not find ",trim(fname)
          stop '=> setmom sed alk'
        endif
#  endif
#  if !defined O_npzd_o2
        fname = new_file_name ("O_o2.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) then
          call openfile (trim(fname), iou)
          call getvara ('O_o2', iou, imtm2*km, ib, ic, tmpik, c1, c0)
          do i=2,imtm1
            k = kmt(i,jrow)
            sbc(i,jrow,ibo2) = tmpik(i,k)
          enddo
        call setbcx (sbc(:,jrow,ibo2), imt, 1)
        else
          print*, "Error => Can not find ",trim(fname)
          stop '=> setmom sed o2'
        endif
#  endif
# endif

!-----------------------------------------------------------------------
!       initialize surface boundary
!-----------------------------------------------------------------------
        if (isst .ne. 0) sbc(1:imt,jrow,isst) = t(1:imt,1,j,itemp,taup1)
        if (isss .ne. 0) sbc(1:imt,jrow,isss) = t(1:imt,1,j,isalt,taup1)
        if (issdic .ne. 0)
     &    sbc(1:imt,jrow,issdic) = t(1:imt,1,j,idic,taup1)
        if (issdic13 .ne. 0)
     &    sbc(1:imt,jrow,issdic13) = t(1:imt,1,j,idic13,taup1)
        if (isso2 .ne. 0) sbc(1:imt,jrow,isso2) = t(1:imt,1,j,io2,taup1)
        if (issalk .ne. 0)
     &    sbc(1:imt,jrow,issalk) = t(1:imt,1,j,ialk,taup1)
        if (isspo4 .ne. 0)
     &    sbc(1:imt,jrow,isspo4) = t(1:imt,1,j,ipo4,taup1)
        if (issdop .ne. 0)
     &    sbc(1:imt,jrow,issdop) = t(1:imt,1,j,idop,taup1)
        if (issphyt .ne. 0)
     &    sbc(1:imt,jrow,issphyt) = t(1:imt,1,j,iphyt,taup1)
        if (isszoop .ne. 0)
     &    sbc(1:imt,jrow,isszoop) = t(1:imt,1,j,izoop,taup1)
        if (issdetr .ne. 0)
     &    sbc(1:imt,jrow,issdetr) = t(1:imt,1,j,idetr,taup1)
        if (issno3 .ne. 0)
     &    sbc(1:imt,jrow,issno3) = t(1:imt,1,j,ino3,taup1)
        if (issdon .ne. 0)
     &    sbc(1:imt,jrow,issdon) = t(1:imt,1,j,idon,taup1)
        if (issdiaz .ne. 0)
     &    sbc(1:imt,jrow,issdiaz) = t(1:imt,1,j,idiaz,taup1)
        if (issdin15 .ne. 0)
     &    sbc(1:imt,jrow,issdin15) = t(1:imt,1,j,idin15,taup1)
        if (issdon15 .ne. 0)
     &    sbc(1:imt,jrow,issdon15) = t(1:imt,1,j,idon15,taup1)
        if (issphytn15 .ne. 0)
     &    sbc(1:imt,jrow,issphytn15) = t(1:imt,1,j,iphytn15,taup1)
        if (isszoopn15 .ne. 0)
     &    sbc(1:imt,jrow,isszoopn15) = t(1:imt,1,j,izoopn15,taup1)
        if (issdetrn15 .ne. 0)
     &    sbc(1:imt,jrow,issdetrn15) = t(1:imt,1,j,idetrn15,taup1)
        if (issdiazn15 .ne. 0)
     &    sbc(1:imt,jrow,issdiazn15) = t(1:imt,1,j,idiazn15,taup1)
        if (issdoc13 .ne. 0)
     &    sbc(1:imt,jrow,issdoc13) = t(1:imt,1,j,idoc13,taup1)
        if (issphytc13 .ne. 0)
     &    sbc(1:imt,jrow,issphytc13) = t(1:imt,1,j,iphytc13,taup1)
        if (isszoopc13 .ne. 0)
     &    sbc(1:imt,jrow,isszoopc13) = t(1:imt,1,j,izoopc13,taup1)
        if (issdetrc13 .ne. 0)
     &    sbc(1:imt,jrow,issdetrc13) = t(1:imt,1,j,idetrc13,taup1)
        if (issdiazc13 .ne. 0)
     &    sbc(1:imt,jrow,issdiazc13) = t(1:imt,1,j,idiazc13,taup1)
        if (issc14 .ne. 0)
     &    sbc(1:imt,jrow,issc14) = t(1:imt,1,j,ic14,taup1)
        if (isscfc11 .ne. 0)
     &    sbc(1:imt,jrow,isscfc11) = t(1:imt,1,j,icfc11,taup1)
        if (isscfc12 .ne. 0)
     &    sbc(1:imt,jrow,isscfc12) = t(1:imt,1,j,icfc12,taup1)

# if defined O_gthflx
!       bottom heat flux
        do i=1,imt
           bhf(i,jrow) = geoheatflux(xt(i),yt(jrow))
        enddo
# endif
!-----------------------------------------------------------------------
!       zero out tracers in land points
!-----------------------------------------------------------------------

        do i=1,imt
          kz = kmt(i,jrow)
          do k=1,km
            if (k .gt. kz) then
              do n=1,nt
                t(i,k,j,n,taup1) = c0
              enddo
            endif
          enddo
        enddo

!-----------------------------------------------------------------------
!       checksum the initial conditions
!-----------------------------------------------------------------------

        ucksum = ucksum + checksum (u(1,1,j,1,taup1), imt, km)
        vcksum = vcksum + checksum (u(1,1,j,2,taup1), imt, km)
        tcksum = tcksum + checksum (t(1,1,j,1,taup1), imt, km)
        scksum = scksum + checksum (t(1,1,j,2,taup1), imt, km)

!-----------------------------------------------------------------------
!       initialize every latitude jrow either on disk (when jmw < jmt)
!       or in the MW (when the last jrow is complete and jmw = jmt)
!-----------------------------------------------------------------------

        if (wide_open_mw) then
          if (jrow .eq. jmt) then
            call copy_all_rows (taup1, tau)
            call copy_all_rows (tau, taum1)
          endif
        else
              call putrow (latdisk(taudisk),  nslab, jrow
     &,                u(1,1,j,1,taup1), t(1,1,j,1,taup1))
              call putrow (latdisk(taup1disk), nslab, jrow
     &,                u(1,1,j,1,taup1), t(1,1,j,1,taup1))
        endif

      enddo

!-----------------------------------------------------------------------
!     find inital average surface references
!-----------------------------------------------------------------------
      print*, " "
      print*, "inital average surface references: "
      dmsk(:,:) = 1.
      where (kmt(:,:) .eq. 0) dmsk(:,:) = 0.
      gaost(:) = 0.
      if (isalt .ne. 0 .and. isss .ne. 0) then
        call areaavg (sbc(1,1,isss), dmsk, gaost(isalt))
        gaost(isalt) = gaost(isalt) + 0.035
        socn = gaost(isalt)
        print*, "global average sea surface salinity (psu) = "
     &,   gaost(isalt)*1000.
      endif
      if (idic .ne. 0 .and. issdic .ne. 0) then
        call areaavg (sbc(1,1,issdic), dmsk, gaost(idic))
        print*, "global average sea surface dic (mol m-3) = "
     &,   gaost(idic)
      endif
      if (idic13 .ne. 0 .and. issdic13 .ne. 0) then
        call areaavg (sbc(1,1,issdic13), dmsk, gaost(idic13))
        print*, "global average sea surface dic 13 (mol m-3) = "
     &,   gaost(idic13)
      endif
      if (io2 .ne. 0 .and. isso2 .ne. 0) then
        call areaavg (sbc(1,1,isso2), dmsk, gaost(io2))
        print*, "global average sea surface oxygen (mol m-3) = "
     &,   gaost(io2)
      endif
      if (ialk .ne. 0 .and. issalk .ne. 0) then
        call areaavg (sbc(1,1,issalk), dmsk, gaost(ialk))
        print*, "global average sea surface alkalinity (mol m-3) = "
     &,   gaost(ialk)
      endif
      if (ipo4 .ne. 0 .and. isspo4 .ne. 0) then
        call areaavg (sbc(1,1,isspo4), dmsk, gaost(ipo4))
        print*, "global average sea surface phosphate (mol m-3) = "
     &,   gaost(ipo4)*0.001
      endif
      if (idop .ne. 0 .and. issdop .ne. 0) then
        call areaavg (sbc(1,1,issdop), dmsk, gaost(idop))
        print*, "global average sea surface DOP (mol m-3) = "
     &,   gaost(idop)*0.001
      endif
# if !defined O_npzd_no_vflux
      if (iphyt .ne. 0 .and. issphyt .ne. 0) then
        call areaavg (sbc(1,1,issphyt), dmsk, gaost(iphyt))
        print*, "global average sea surface phytoplankton (mol m-3) = "
     &,   gaost(iphyt)*0.001
      endif
      if (izoop .ne. 0 .and. isszoop .ne. 0) then
        call areaavg (sbc(1,1,isszoop), dmsk, gaost(izoop))
        print*, "global average sea surface zooplankton (mol m-3) = "
     &,   gaost(izoop)*0.001
      endif
      if (idetr .ne. 0 .and. issdetr .ne. 0) then
        call areaavg (sbc(1,1,issdetr), dmsk, gaost(idetr))
        print*, "global average sea surface detritus (mol m-3) = "
     &,   gaost(idetr)*0.001
      endif
# endif
      if (ino3 .ne. 0 .and. issno3 .ne. 0) then
        call areaavg (sbc(1,1,issno3), dmsk, gaost(ino3))
        print*, "global average sea surface nitrate (mol m-3) = "
     &,   gaost(ino3)*0.001
      endif
      if (idon .ne. 0 .and. issdon .ne. 0) then
        call areaavg (sbc(1,1,issdon), dmsk, gaost(idon))
        print*, "global average sea surface DON (mol m-3) = "
     &,   gaost(idon)*0.001
      endif
# if !defined O_npzd_no_vflux
      if (idiaz .ne. 0 .and. issdiaz .ne. 0) then
        call areaavg (sbc(1,1,issdiaz), dmsk, gaost(idiaz))
        print*, "global average sea surface diazotrophs (mol m-3) = "
     &,   gaost(idiaz)*0.001
      endif
# endif
      if (idin15 .ne. 0 .and. issdin15 .ne. 0) then
        call areaavg (sbc(1,1,issdin15), dmsk, gaost(idin15))
        print*, "global average sea surface nitrate 15 (mol m-3) = "
     &,   gaost(idin15)*0.001
      endif
      if (idon15 .ne. 0 .and. issdon15 .ne. 0) then
        call areaavg (sbc(1,1,issdon15), dmsk, gaost(idon15))
        print*, "global average sea surface DON15 (mol m-3) = "
     &,   gaost(idon15)*0.001
      endif
# if !defined O_npzd_no_vflux
      if (iphytn15 .ne. 0 .and. issphytn15 .ne. 0) then
        call areaavg (sbc(1,1,issphytn15), dmsk, gaost(iphytn15))
        print*, "global average sea surface phytoplankton n15 " 
     &,         "(mol m-3) = ", gaost(iphytn15)*0.001
      endif
      if (izoopn15 .ne. 0 .and. isszoopn15 .ne. 0) then
        call areaavg (sbc(1,1,isszoopn15), dmsk, gaost(izoopn15))
        print*, "global average sea surface zooplankton n15 (mol m-3) 
     &           = ", gaost(izoopn15)*0.001
      endif
      if (idetrn15 .ne. 0 .and. issdetrn15 .ne. 0) then
        call areaavg (sbc(1,1,issdetrn15), dmsk, gaost(idetrn15))
        print*, "global average sea surface detritus n15 (mol m-3) 
     &           = ", gaost(idetrn15)*0.001
      endif
      if (idiazn15 .ne. 0 .and. issdiazn15 .ne. 0) then
        call areaavg (sbc(1,1,issdiazn15), dmsk, gaost(idiazn15))
        print*, "global average sea surface diazotrophs n15 (mol m-3) 
     &           = ", gaost(idiazn15)*0.001
      endif
# endif
      if (idoc13 .ne. 0 .and. issdoc13 .ne. 0) then
        call areaavg (sbc(1,1,issdoc13), dmsk, gaost(idoc13))
        print*, "global average sea surface DOC13"
     &,         " (mol m-3) = ", gaost(idoc13)*0.001
      endif
# if !defined o_npzd_no_vflux
      if (iphytc13 .ne. 0 .and. issphytc13 .ne. 0) then
        call areaavg (sbc(1,1,issphytc13), dmsk, gaost(iphytc13))
        print*, "global average sea surface phytoplankton C13"
     &,         " (mol m-3) = ", gaost(iphytc13)*0.001
      endif
      if (izoopc13 .ne. 0 .and. isszoopc13 .ne. 0) then
        call areaavg (sbc(1,1,isszoopc13), dmsk, gaost(izoopc13))
        print*, "global average sea surface zooplankton C13"
     &,         " (mol m-3) = ", gaost(izoopc13)*0.001
      endif
      if (idetrc13 .ne. 0 .and. issdetrc13 .ne. 0) then
        call areaavg (sbc(1,1,issdetrc13), dmsk, gaost(idetrc13))
        print*, "global average sea surface detritus c13"
     &,         " (mol m-3) = ", gaost(idetrc13)*0.001
      endif
      if (idiazc13 .ne. 0 .and. issdiazc13 .ne. 0) then
        call areaavg (sbc(1,1,issdiazc13), dmsk, gaost(idiazc13))
        print*, "global average sea surface diazotrophs c13"
     &,         " (mol m-3) = ", gaost(idiazc13)*0.001
      endif
# endif
      if (ic14 .ne. 0 .and. issc14 .ne. 0) then
        call areaavg (sbc(1,1,issc14), dmsk, gaost(ic14))
        print*, "global average sea surface carbon 14 (mol m-3) = "
     &,   gaost(ic14)
      endif
      if (icfc11 .ne. 0 .and. isscfc11 .ne. 0) then
        call areaavg (sbc(1,1,isscfc11), dmsk, gaost(icfc11))
        print*, "global average sea surface cfc 11 (mol m-3) = "
     &,   gaost(icfc11)
      endif
      if (icfc12 .ne. 0 .and. isscfc12 .ne. 0) then
        call areaavg (sbc(1,1,isscfc12), dmsk, gaost(icfc12))
        print*, "global average sea surface cfc 12 (mol m-3) = "
     &,   gaost(icfc12)
      endif
      print*, " "

      write (stdout,*) ' I.C. checksum for t =',tcksum
      write (stdout,*) ' I.C. checksum for s =',scksum
      write (stdout,*) ' I.C. checksum for u =',ucksum
      write (stdout,*) ' I.C. checksum for v =',vcksum

      return
      end

      function theta0 (ydeg, depth)

!=======================================================================
!     this subroutine returns estimates of global mean potential
!     temperature for model initialization as a function of depth.
!     it is used to produce a reference thermal stratification for the
!     upper 2000m of the MOM`s test case.  below 2000m, the
!     potential temperature returned is 2.0 degrees C.  surface
!     values are set slightly above 18.4 degrees C at the reference
!     latitude "reflat".
!     the estimates are produced from a 7th order polynomial fit to
!     the annual mean world ocean potential temperature observations
!     of Levitus (1982).

!     input [units]:
!       a latitude (ydeg): [degrees]
!       a zt value (depth): [centimetres]
!     output [units]:
!       potential temperature estimate (est): [degrees centigrade]

!     variables:
!       coeft     = coefficients for the polynomial fit of potential
!                   temperature vs. depth
!       reflat    = reference latitude at which observed surface
!                   temperatures approximately equal coeft(1)
!       factor    = the ratio of the cosine of the latitude requested
!                   ("ydeg") to the reference latitude ("reflat")
!                   used to scale the upper 2000 meters of the vertical
!                   temperature profile
!       tmin,tmax = the minimum and maximum potential temperatures
!                   allowed at the time of model initialization

!     reference:
!       Levitus, S., Climatological atlas of the world ocean, NOAA
!     Prof. Paper 13, US Gov`t printing Office, Washington, DC, 1982.
!=======================================================================

      implicit none

      integer ndeg, nn
      parameter (ndeg=7)

      real c0, coslat, factor, pi, theta0, tmin, tmax, refcos, reflat
      real ydeg, z, depth, est, bb
      real coeft(ndeg+1)

      save coeft, factor, pi, tmin, tmax, reflat

      data coeft / 0.184231944E+02,-0.430306621E-01, 0.607121504E-04
     &           ,-0.523806281E-07, 0.272989082E-10,-0.833224666E-14
     &           , 0.136974583E-17,-0.935923382E-22/

      data tmin, tmax, reflat /2.0, 25.0, 34.0/

      c0 = 0.0
      pi = atan(1.0) * 4.0
      refcos = abs(cos(pi*reflat/180.))

      coslat = abs(cos(pi*ydeg/180.))
      factor = coslat/refcos
      z = depth * 0.01

      if (z .gt. 2000.) then
        est = 2.0
      else
        est = c0
        bb = 1.0
        do nn=1,ndeg+1
!          if (nn.gt.1) bb = z**(nn-1)
          est = est + coeft(nn)*bb
          bb = bb*z
        enddo
        est = est * factor
      endif

      if (est .gt. tmax) est = tmax
      if (est .lt. tmin) est = tmin

      theta0 = est

#endif
      return
      end
