!=======================================================================
      subroutine sed (is, ie, js, je)

#if defined O_sed
!-----------------------------------------------------------------------
!     sediment model
!-----------------------------------------------------------------------

      implicit none

      integer ie, is, je, js, i, ip, i_year, j, k, n_control, n_debug
      integer n_year, nsed

      real rtsed, r13or, r14or, r13ca, r14ca, rainr, tmp

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "tmngr.h"
      include "switch.h"
      include "sed.h"
      include "csbc.h"
      include "grdvar.h"
      include "levind.h"

      integer n_buried_depths(ipmax)

!     results from sediment
      real ttrtc(ipmax), ttral(ipmax), difal(ipmax), diftc(ipmax)

!     internal variables to sediment
      real temp_p(ipmax), sal_p(ipmax), alk_p(ipmax), tco2_p(ipmax)
      real o2_bw_p(ipmax), rain_clay_p(ipmax), co2_p(ipmax)
      real hco3_p(ipmax), resp_c(ipmax,nzmax,3), cal_c(ipmax,nzmax)

      real dmsk(imt,jmt), tmpij(imt,jmt), tmpip(ipmax)

!-----------------------------------------------------------------------
!     n_control = 0  initialize without setting calcite concentration
!     n_control = 1  find steady state calcite concentration
!     n_control = 2  step ahead time-dependent calcite concentration
!     n_control = 3  maintain constant calcite concentration
!     n_control = 4  return all calcite rain to the overlying water
!     n_control = 5  use old reaction rates to save time but update bury
!-----------------------------------------------------------------------
      n_control = 2

!-----------------------------------------------------------------------
!     you must uncomment lines in sed.f and sediment.f when debugging.
!     Look for: "!!     uncomment for debugging"
!     n_debug = 0  if things are working well
!     n_debug = 1  watch the convergence in stdout
!     n_debug = 2  when convergence is slow, see one a site converging
!     n_debug = 3  when code is blowing up, see everything for all sites
!     n_debug > 3  specific site index
!-----------------------------------------------------------------------
      n_debug = 0

!-----------------------------------------------------------------------
!     transfer 2d ocean bottom variables to 1d arrays
!     convert concentrations from mol m-3 (or umol cm-3) to mol l-1
!-----------------------------------------------------------------------
      rtsed = 1.
!     set the ratio between calcite and organic rain rates
      if (atsed .gt. 0) rtsed = 1./atsed
      do ip=1,ipsed
        i = imap(ip)
        j = jmap(ip)
        temp_p(ip) = sbc(i,j,ibtemp)*rtsed
!       convert from mom units of salinity to psu
        sal_p(ip) = sbc(i,j,ibsalt)*rtsed*1000. + 35.
        alk_p(ip) = sbc(i,j,ibalk)*rtsed*0.001
        tco2_p(ip) = sbc(i,j,ibdic)*rtsed*0.001
        o2_bw_p(ip) = max(sbc(i,j,ibo2)*rtsed*0.001, 1.e-6)
!       convert from umol/cm2/s to mol/cm2/dtsed
        rain_cal_p(ip) = sbc(i,j,ircal)*rtsed*1.e-6*dtsed
        rain_org_p(ip) = sbc(i,j,irorg)*rtsed*1.e-6*dtsed
# if defined O_sed_constrain_rainr
!       set limits on rain ratio by limiting effective organic rain
        rain_org_p(ip) = min(rain_org_p(ip),rain_cal_p(ip)*10.0)
        rain_org_p(ip) = max(rain_org_p(ip),rain_cal_p(ip)*0.1)
!       bring the rain ratio closer to 1 by scaling effective organic rain
        if (rain_org_p(ip) .lt. rain_cal_p(ip)) then
          rain_org_p(ip) = (rain_org_p(ip) - rain_cal_p(ip))*.5
     &                   + rain_cal_p(ip)
        else
          rain_org_p(ip) = 1./((1./rain_org_p(ip)-1./rain_cal_p(ip))*.5
     &                   + 1./rain_cal_p(ip))
        endif
# endif
!       parametrize clay rain from calcite
        rain_clay_p(ip) = 9.1e-6*(rain_org_p(ip)*1.e6)**1.41
# if defined O_npzd_alk
        sbc(i,j,ibalkfx) = 0.
# endif
# if defined O_carbon
        sbc(i,j,ibdicfx) = 0.
# endif
      enddo
# if defined O_time_averages

!-----------------------------------------------------------------------
!     zero time averages if not in an averaging period
!-----------------------------------------------------------------------
      if (.not. timavgperts) call ta_sed_tavg (is, ie, js, je, 0)
# endif
# if defined O_time_step_monitor

!-----------------------------------------------------------------------
!     zero time step integrals if not in an averaging period
!-----------------------------------------------------------------------
      if (.not. tsiperts) call ta_sed_tsi (is, ie, js, je, 0)
# endif

!!     uncomment for debugging
!      if (n_debug .gt. 3) then
!        print*, "input at: ", n_debug
!        print*, "temp_p: ", temp_p(n_debug)
!        print*, "sal_p: ", sal_p(n_debug)
!        print*, "alk_p: ", alk_p(n_debug)
!        print*, "tco2_p: ", tco2_p(n_debug)
!        print*, "o2_bw_p: ", o2_bw_p(n_debug)
!        print*, "rain_cal_p: ", rain_cal_p(n_debug)
!        print*, "rain_org_p: ", rain_org_p(n_debug)
!        print*, "rain_clay_p: ", rain_clay_p(n_debug)
!      endif

!-----------------------------------------------------------------------
!     update sediment year for sediment profiles
!-----------------------------------------------------------------------
# if defined O_sed_profile
      sed_year = sed_year + dtsed
# else
      sed_year = 0.
# endif
      n_year = nint(sed_year)

!-----------------------------------------------------------------------
!     loop through sediment model nsedacc times (>1 during spinup)
!-----------------------------------------------------------------------
      do nsed=1,nsedacc

        call calc_k (temp_p, sal_p, water_z_p, k1, k2, k3, csat, ipsed)

        call calc_buff (alk_p, tco2_p, sal_p, k1, k2, k3, co2_p, hco3_p
     &,                 co3_p, ipsed)

        call setup_pw (co2_p, hco3_p, co3_p, o2_bw_p, carb, o2, csat
     &,                kmax, ipsed, ipmax, nzmax)

        call estimate_rc (rain_org_p, rc, ipsed)

        if (n_control .le. 0) then

!-----------------------------------------------------------------------
!         then initialize without setting calcite concentration
!-----------------------------------------------------------------------
          do ip=1,ipsed
            buried_mass(1,ip) = 500.
            buried_calfrac(1,ip) = calgg(kmax,ip)
            do k=2,ibmax
               buried_mass(k,ip) = 0.
               buried_calfrac(k,ip) = calgg(kmax,ip)
            enddo
          enddo
          call set_pore (calgg, zsed, pore, kmax, ipsed, ipmax, nzmax)
          call get_sed_ml_mass (delz, pore, sed_ml_mass, nzmax, ipsed
     &,                         ipmax, kmax)

        elseif (n_control .eq. 1) then

!-----------------------------------------------------------------------
!         find steady state calcite concentration
!-----------------------------------------------------------------------
!!         uncomment for debugging
!          if (n_debug .ge. 1) write(6,*) 'Starting calss'

          call sed_ss (zsed, delz, form, pore, kmax, o2, zrct, carb
     &,                orgml, orggg, calml, calgg, rain_cal_p
     &,                rain_org_p, rain_clay_p, r13or, r14or, r13ca
     &,                r14ca, rc, dissc, dissn, csat, k1, k2, dopls
     &,                domin, dcpls, dcmin, dbpls, dbmin, resp_c, cal_c
     &,                ttrorg, ttrcal, ttrtc, ttral, diftc, difal
     &,                c_advect, buried_mass, buried_calfrac
     &,                n_buried_depths, ipsed, ipmax, nzmax, ibmax
     &,                n_control, n_debug)

!         reinitialize without setting calcite concentration
          do ip=1,ipsed
            buried_mass(1,ip) = 500.
            buried_calfrac(1,ip) = calgg(kmax,ip)
            do k=2,ibmax
               buried_mass(k,ip) = 0.
               buried_calfrac(k,ip) = calgg(kmax,ip)
            enddo
          enddo
          call set_pore (calgg, zsed, pore, kmax, ipsed, ipmax, nzmax)
          call get_sed_ml_mass (delz, pore, sed_ml_mass, nzmax, ipsed
     &,                         ipmax, kmax)

        elseif (n_control .eq. 2) then

!-----------------------------------------------------------------------
!         step ahead time-dependent calcite concentration
!-----------------------------------------------------------------------
          call sed_const_cal (zsed, delz, form, pore, kmax, o2, zrct
     &,                       carb, orgml, orggg, calml, calgg
     &,                       rain_cal_p, rain_org_p, rain_clay_p, r13or
     &,                       r14or, r13ca, r14ca, rc, dissc, dissn
     &,                       csat, k1, k2, dopls, domin, dcpls, dcmin
     &,                       dbpls, dbmin, resp_c, cal_c, ttrorg
     &,                       ttrcal, ttrtc, ttral, diftc, difal
     &,                       c_advect, ipsed, ipmax, nzmax, n_control
     &,                       n_debug)

          call bury (n_control, zsed, delz, pore, kmax, calgg, orggg
     &,              sed_ml_mass, rain_cal_p, ttrcal, c_advect
     &,              rain_clay_p, n_year, buried_mass, buried_calfrac
     &,              depth_age, ipsed, ipmax, nzmax, ibmax)

        elseif (n_control .eq. 3) then

!-----------------------------------------------------------------------
!        maintain constant calcite concentration (update diss but don't bury)
!-----------------------------------------------------------------------
         call sed_const_cal (zsed, delz, form, pore, kmax, o2, zrct
     &,                      carb, orgml, orggg, calml, calgg
     &,                      rain_cal_p, rain_org_p, rain_clay_p, r13or
     &,                      r14or, r13ca, r14ca, rc, dissc, dissn
     &,                      csat, k1, k2, dopls, domin, dcpls, dcmin
     &,                      dbpls, dbmin, resp_c, cal_c, ttrorg
     &,                      ttrcal, ttrtc, ttral, diftc, difal
     &,                      c_advect, ipsed, ipmax, nzmax, n_control
     &,                      n_debug)

        elseif (n_control .eq. 4) then

!-----------------------------------------------------------------------
!         return all calcite rain to the overlying water
!-----------------------------------------------------------------------
          do ip=1,ipsed
            ttrcal(ip) = rain_cal_p(ip)
          enddo

        elseif (n_control .eq. 5) then

!-----------------------------------------------------------------------
!         use old reaction rates but update bury
!-----------------------------------------------------------------------
          call bury (n_control, zsed, delz, pore, kmax, calgg, orggg
     &,              sed_ml_mass, rain_cal_p, ttrcal, c_advect
     &,              rain_clay_p, n_year, buried_mass, buried_calfrac
     &,              depth_age, ipsed, ipmax, nzmax, ibmax)

        endif

      enddo

!-----------------------------------------------------------------------
!     set variables for coupling back to the ocean model
!-----------------------------------------------------------------------
# if defined O_sed_weath_diag
      weathflx = 0.
# endif
# if defined O_npzd_alk
      sbc(:,:,ibalkfx) = 0.
# endif
# if defined O_carbon
      sbc(:,:,ibdicfx) = 0.
# endif
      if (atsed .gt. 0) then
        rtsed = 1.e6/dtsed
        do ip=1,ipsed
          i = imap(ip)
          j = jmap(ip)
          tmp =  (ttrcal(ip) - rain_cal_p(ip))*rtsed
# if defined O_sed_weath_diag
!         diagnose weathering flux as he total burial flux (umol/s)
          weathflx = weathflx - tmp*dxt(i)*dyt(j)*cst(j)
# endif
!         set source terms, convert to umol cm-3 s-1
!         the source terms are really a correction to assumed zero burial
!         during the previous coupling time and is applied over the next
          tmp = tmp*dztr(kmt(i,j))
# if defined O_npzd_alk
          sbc(i,j,ibalkfx) = 2.*tmp
# endif
# if defined O_carbon
          sbc(i,j,ibdicfx) = tmp
# endif
!         zero accumulators
          sbc(i,j,ibtemp) = 0.
          sbc(i,j,ibsalt) = 0.
# if defined O_npzd_alk
          sbc(i,j,ibalk) = 0.
# endif
# if defined O_carbon
          sbc(i,j,ibdic) = 0.
# endif
# if defined O_npzd_o2
          sbc(i,j,ibo2) = 0.
# endif
          sbc(i,j,ircal) = 0.
          sbc(i,j,irorg) = 0.
        enddo
        atsed = 0.
      endif
#endif

      return
      end
