! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/sed/setsed.F
      subroutine setsed (is, ie, js, je)

!-----------------------------------------------------------------------
!     set up everything which must be done only once per run
!-----------------------------------------------------------------------

      implicit none

      integer is, ie, js, je

      character(120) :: fname, new_file_name

      integer i, ioun, j, k

      logical exists, inqvardef

      real expb, db, difo2, difc(3)

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

      dtsed = dtsedyr*yrlen*daylen
      nsedacc = max(nsedacc, 1)

      if (nsedacc .gt. 1) then
        print*, '==> Warning: sediment model accelerated and coupled'
      endif

!-----------------------------------------------------------------------
!     initialise sediment model
!-----------------------------------------------------------------------
      weathflx = 0.  ! weathering flux in umol s-1
      sed_year = 0.
      carb(:,1,:) = 20.E-6
      carb(:,2,:) = 2000.E-6
      carb(:,3,:) = 80.E-6
      dcpls(:,:,:) = 0.
      dcmin(:,:,:) = 0.
      pore(:,:) = 0.
      form(:,:) = 0.
      o2(:,:) = 150.e-6
      orgml(:,:) = 0.
      calml(:,:) = 0.
      dopls(:,:) = 0.
      domin(:,:) = 0.
      dbpls(:,:) = 0.
      dbmin(:,:) = 0.
      co3_p(:) = 0.
      atsed = 0
      ipsed = 0
      map_sed(:,:) = 0
      calgg(:,:) = 0.1
      orggg(:,:) = 0.0002
      buried_mass(:,:) = 0.
      buried_calfrac(:,:) = 0.1
      buried_mass(1,:) = 500.
      kmin = max(1,kmin)
      ttrcal(:) = 0.
      rain_cal_p(:) = 0.
      sedsa = 0.
      carblith = 0.
      do j=2,jmtm1
        do i=2,imtm1
          if (kmt(i,j) .ge. kmin) then
            ipsed = ipsed + 1
            map_sed(i,j) = ipsed
            imap(ipsed) = i
            jmap(ipsed) = j
            water_z_p(ipsed) = zw(kmt(i,j))*0.01
	    sedsa = sedsa + dxt(i)*dyt(j)*cst(j)
          endif
        enddo
      enddo

      delz(:) = 3.
      delz(1:5) = (/0.,.5,.5,1.,2./)
      kmax = 7
      zsed(1) = 0.
      do k=2,kmax
        zsed(k) = zsed(k-1) + delz(k)
      enddo
      zrct(:) = zsed(kmax)
      dissc = 1.1574e-5
      dissn = 4.5

!     depth_age is only used with the option "sed_profile".
!     make sure you have enough levels for the profile you want to simulate.
!     sediment is stored inverse to the depth_age. the latest sediment
!     is stored incrementing toward ibmax. once ibmax is reached all new
!     sediment is stored at ibmax. for example if sed_year is 120000 and
!     ibmax is 20 with a time increment of 4000, then the sediment in
!     level 1 is older than 120000 (usually the initial condition), in
!     level 2 it is between 120000 and 116000, level 3 between 116000 and
!     112000, level ibmax-1 between 76000 and 72000 and level ibmax between
!     76000 and 0. if sed_year is 20000 then the sediment in level 1 is
!     older than 20000, in level 2 it is between 20000 and 16000, in
!     level 5 between 8000 and 4000, in level 6 between 4000 and 0 and
!     levels higher than 6 would be available for future accumulation.
      do k=1,ibmax
        depth_age(k) = 4000.*(k-1)
      enddo

      expb = 3.
      difo2 = 12.1e-6
      difc(1) = 10.5e-6
      difc(2) = 6.4e-6
      difc(3) = 5.2e-6
      db = 0.15

      call set_pore (calgg, zsed, pore, kmax, ipsed, ipmax, nzmax)
      call pore_2_form (pore, form, kmax, expb, ipsed, ipmax, nzmax)
      call sldcon (orgml, orggg, 12., pore, kmax, ipsed, ipmax, nzmax)
      call sldcon (calml,calgg, 100., pore, kmax, ipsed, ipmax, nzmax)
      call calc_do2 (difo2, form, pore, delz, kmax, dopls, domin, ipsed
     &,              ipmax, nzmax)
      call calc_dc (difc, form, pore, delz, kmax, dcpls, dcmin, ipsed
     &,             ipmax, nzmax)
      call calc_db (db, pore, zsed, delz, kmax, dbpls, dbmin, ipsed
     &,             ipmax, nzmax)
!     calculate mass of mixed layer (g/cm2)
      call get_sed_ml_mass (delz, pore, sed_ml_mass, nzmax, ipsed
     &,                     ipmax, kmax)

!     zero boundary condition accumulators
      sbc(:,:,ibtemp) = 0.
      sbc(:,:,ibsalt) = 0.
      sbc(:,:,ircal) = 0.
      sbc(:,:,irorg) = 0.

      sbc(:,:,ibdic) = 0.

      sbc(:,:,ibalk) = 0.

      sbc(:,:,ibo2) = 0.

      if (.not. init) then
        fname = new_file_name ("restart_sed.nc")
        inquire (file=trim(fname), exist=exists)
        if (exists) call sed_rest_in (fname, is, ie, js, je)

      endif

!-----------------------------------------------------------------------
!     if defined, use namelist value for constant weathering flux
!-----------------------------------------------------------------------

!     convert from kg/s to umol/s
      if (weath .lt. 1.e20) weathflx = weath*1.e9/12.

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

      call ta_sed_tavg (is, ie, js, je, 0)

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

      call ta_sed_tsi (is, ie, js, je, 0)

      return
      end