! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/source/sed/sed.F !======================================================================= subroutine sed (is, ie, js, je) !----------------------------------------------------------------------- ! 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 ! parametrize clay rain from calcite rain_clay_p(ip) = 9.1e-6*(rain_org_p(ip)*1.e6)**1.41 sbc(i,j,ibalkfx) = 0. sbc(i,j,ibdicfx) = 0. enddo !----------------------------------------------------------------------- ! zero time averages if not in an averaging period !----------------------------------------------------------------------- if (.not. timavgperts) call ta_sed_tavg (is, ie, js, je, 0) !----------------------------------------------------------------------- ! zero time step integrals if not in an averaging period !----------------------------------------------------------------------- if (.not. tsiperts) call ta_sed_tsi (is, ie, js, je, 0) !! 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 !----------------------------------------------------------------------- sed_year = 0. 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 !----------------------------------------------------------------------- sbc(:,:,ibalkfx) = 0. sbc(:,:,ibdicfx) = 0. 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 ! 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)) sbc(i,j,ibalkfx) = 2.*tmp sbc(i,j,ibdicfx) = tmp ! zero accumulators sbc(i,j,ibtemp) = 0. sbc(i,j,ibsalt) = 0. sbc(i,j,ibalk) = 0. sbc(i,j,ibdic) = 0. sbc(i,j,ibo2) = 0. sbc(i,j,ircal) = 0. sbc(i,j,irorg) = 0. enddo atsed = 0. endif return end