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