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