! source file: /home/aschmitt/UVic/2.9/chris/c_n_isotopes/asc/MOBI1.5/tidal_3DG/lgm/updates/setmom.F
      subroutine setmom (is, ie, js, je)

!=======================================================================
!     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
      include "coord.h"
      include "cpolar.h"
      include "cprnts.h"
      include "cregin.h"
      include "csbc.h"
      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"
      include "isopyc.h"
      include "levind.h"
      include "mw.h"
      include "scalar.h"
      include "stab.h"
      include "state.h"
      include "switch.h"
      include "tmngr.h"
      include "vmixc.h"
      integer kmz(imt,jmt)
      include "tidal_kv.h"
      real tmpij(imtm2,jmtm2)

      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)
!      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

!     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))

!-----------------------------------------------------------------------
!     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
      fname = new_file_name ("O_tidenrg_green.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

!-----------------------------------------------------------------------
!     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 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
!-----------------------------------------------------------------------

!     initialize two dimensional fields on disk

!     initialize stream function fields in memory

      do n=1,2
        do jrow=1,jmt
          do i=1,imt
            psi(i,jrow,n) = c0
            ptd(i,jrow)   = c0
          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

!     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)
      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

!-----------------------------------------------------------------------
!     zero integrated time average accumulators
!-----------------------------------------------------------------------

      call ta_mom_tsi (0)

!-----------------------------------------------------------------------
!     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)

!-----------------------------------------------------------------------
!     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

!-----------------------------------------------------------------------
!     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 (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
      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

!     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 (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 = ')

!-----------------------------------------------------------------------
!     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) = ')

!-----------------------------------------------------------------------
!     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
      do jrow=1,jmt
        kmz(1,jrow) = kmz(imtm1,jrow)
      enddo

!---------------------------------------------------------------------
!     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)
      call findex (kmz, jmtfil, 1, jft1, jft2, imt, iszf, iezf)
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 ==')

!---------------------------------------------------------------------
!     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
          ztd(i,jrow) = c0
          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

      amix = am
      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

!     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
!-----------------------------------------------------------------------

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

      call avgi

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

      call isopi (error, am, ah)
!-----------------------------------------------------------------------
!    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

!-----------------------------------------------------------------------
!     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
      do jrow=1,jmt
        kmu(imt,jrow) = kmu(2,jrow)
      enddo

!---------------------------------------------------------------------
!     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
      real geoheatflux

      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"
      include "diag.h"
      include "state.h"
      include "dens.h"

      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

          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

        enddo
        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
        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)
             t(i,k,j,idiazc13,taup1) = t(i,k,j,idiaz,taup1)*6.625e-3
     &            *rc13std/(1. + rc13std)
           enddo
        enddo

!       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

!-----------------------------------------------------------------------
!       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.

!-----------------------------------------------------------------------
!       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)

!       bottom heat flux
        do i=1,imt
           bhf(i,jrow) = geoheatflux(xt(i),yt(jrow))
        enddo
!-----------------------------------------------------------------------
!       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 (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 (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 (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 (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
      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

      return
      end
