***************************************************************************** *HSox_O2_AJS calculates concentrations and sulfur isotope ratios of disolved *sulfide and sulfate in pore water. H2S IS OXIDIZED BY OXYGEN. The rate is limited *by availability of oxygen provided by the water column above, and supplied *soley through molecular diffusion and mixing (bioturbation/bioirrigation). *Boundary conditions for O2 are the concentration in the Water Column. Oxygen *is only consumed by sulfide oxidation and not by oxic respiration. *32SO4, 34SO4, 32H2S and 34H2S are treated separately resulting *in parallel oxidation and reduction reactions, as well as diffusion of *individual isotopically-substituted species. Although diffusion coefficients *of isotopically-substituted species are the same, isotopic fractionation due *to different concentration gradients is explicitly calculated. Sulfur fluxes *& concentrations are converted from delta notation to isotopic ratios, then *to 32S and 34S. Boundary conditions are seawater sulfate concentrations and *isotopic ratios at the SWI. H2S at SWI = 0. The bottom (~1m) is a no-flow *boundary. SO4 is inititalized at seawater values. H2S = 1mmol/L *Isotopic fractionation is expressed using epsilon notation, where *epsilon = (alpha - 1) * 1000. 15 layers of variable thickness. *Steady state solutions after 10,000 year are reported *Rates of oxidation and reduction are calculated in subroutine REDOXCONVERT *Rates of reduction decline exponentially with depth below the SWI. *[sord = g2 + sord0 * EXP(-mpdepth) / sorddec)]; where sord0 is the rate *of reduction at the SWI, mpdepth is the depth of the midpoint of the layer, *g2 is rate of reduction due to a very slowly decomposing fraction of the *organic matter, and sorddec is the 1/e folding depth. Rates of oxidation *are a function of concentrations of sulfide and oxygen: [hsox = hs * o2]' *where hsox is the rate of xidation, hs and o2 are concentrations of oxygen *and sulfide, respectively. This model is companion to another that calculates *the same values, but for oxidation coupled to nitrate reduction. *Uses the GAUSSND and SLOPERND solver for a nearly diagonal matrix in an *to caluculate concentrations in an implicit scheme: Reverse Euler. See *Walker, JCG, Numerical Adv. with Geochemical Cycles. 1991. Oxford. pp.192 *Multiple sensitivity runx can be made by manipulating the 'Do While' statement. *If you have questions: nsuits@atmos.colostate.edu. *This model is used in a manuscript submitted to AJS **This program runs in FORTRAN. Neil S. Suits ************************************************************************ PROGRAM Sulfur Implicit none Integer :: jlay,jrow,count,numlay=15, & numspec=5,nrow=75 c print*, numlay, numspec, nrow c stop C The following are double precision real number variables. Double Precision :: laydepth,hs320,hs340,so320, & so340,hs0,so40,dhs0,dso0,sord0,sorddec, & g2,hsox0,hsoxdec,delsord0,delhsox0, & x,delx,o20,bioturb C The following variables are dimensioned by nrow=numlay*numspec Double Precision, Dimension(75) :: unk,y,dely,yp,incind C The following variables are dimensioned by numlay Double Precision, Dimension(15) :: hsox,sord,transso,transhs, & mpdepth,delsord,layt,laysep, & so,hs,dso,dhs,rso,rhs,rsord,rhsox,so32, & so34,hs32,hs34,sord32,sord34, & hsox32,hsox34,o2,transo2,delhsox C The following are the constants used in the calculations. Ref is the C 34S/32S ratio of CDT (Ault and Jensen, 1963); secpy is the number of C seconds in a year; difsedso is the effective molecular diffusion coeffient C of dissolved sulfate in seawater at ~10 degrees C, and corrected for C porosity (units: cm^2/sec., Rowe and Howarth, 1985. Deep Sea Res.); C difsedhs is the same, but for HS, and is expressed as a ratio of the MDC C of sulfate (Jorgensen,1979, GCA); difsedsosed, difsedhssed and difsedo2sed convert C units to cm^2/year; transconso, transconhs and transcono2 correct for C the porosity; xstart is the start time in years; xend is the ending time C in years; mstep is the maximum size of the timestep. Double Precision, Parameter ::ref = .0450045, ! R-PDB & secpy = (365.25*3600.*24.), & prsty = .85, ! sediment porosity & difsedso = .0000043, !mol. Diff. Coeff. for SO4 corrected !for tortuosity D'=Theta^2 * D-meas & difsedhs= (1.7*difsedso), & difsedo2= .000001445, ! Diffusion Coefficient of O2 ! in H2O @10deg C, Himmelblau (1964) ! corr. for tortuosity .00002 * .85^2 & difsedsosed = (difsedso*secpy), & difsedhssed = (difsedhs*secpy), & difsedo2sed = (difsedo2*secpy), & transconso = (difsedsosed/prsty), & transconhs = (difsedhssed/prsty), & transcono2 = (difsedo2sed/prsty), & xstart = 1.,xend = 10000., & mstep = 500., & maxinc = .1 C ===================================================================== C thicknessses of layers are set. DATA layt/.5,.5,.5,1.,1.,1.5,2.5,3.5,5.,7.5,10.,10.,20.,20.,20./ laysep(1) = layt(1) / 2. DO jlay = 2, numlay laysep(jlay) = (layt(jlay - 1) + layt(jlay)) / 2. END DO C Effective diffusion coefficients and layer midpoints are calculated. laydepth = 0. DO jlay = 1, numlay transso(jlay) = transconso / layt(jlay) transhs(jlay) = transconhs / layt(jlay) transo2(jlay) = transcono2 / layt(jlay) laydepth = laydepth + laysep(jlay) mpdepth(jlay) = laydepth END DO c Effects of bioturbation on effective diffusion coefficients can be installed here c Bioturbation/bioirrigation cause mixing rates to increase c bioturb = 1.0 c DO jlay = 1, numlay c IF (jlay .le. 6) THEN c transso(jlay) = transso(jlay) * bioturb c transhs(jlay) = transhs(jlay) * bioturb c transo2(jlay) = transo2(jlay) * bioturb c ELSE c transso(jlay) = transso(jlay) c transhs(jlay) = transhs(jlay) c transo2(jlay) = transo2(jlay) c END IF c END DO C Values for any particular run are set here. sord0 is the maximum rate of C sulfate reduction. sorddec is the 1/e folding depth of the depth-dependent C rate of reduction. hsox0 and hsoxdec are analogous values for sulfide C oxidation. delsord0 and delhsox0 are the isotopic fractionation factors C during reduction and oxidation expressed using epsilon notation, where C epsilon = (alpha - 1) * 1000. Units are per mil. g2 is a flux of sulfide C which is independent of depth. It represents sulfide produced by slowly- C degrading fraction of organic matter. sord0 = 0.0005 sorddec = 10.0 g2 = 0. hsox0 = 0. ! this value is not used in O2 simulations hsoxdec = 0.0 ! this value is not used in O2 simulations delsord0= -40. delhsox0 = 20. C Initial and boundary values set, SO4 & 34S-SO4 = sea water, C boundary H2S = 0 mmol/L, 34S-H2S = -25 per mil C initial H2S = .1 mmol/L, 34S-H2S = -25 per mil C these are converted to 32S and 34S concentrations C Units are mmol/L and per mil vs CDT hs0 = 0.00 so40 = 28.9 dhs0 = -25. dso0 = 20. o20 = .278 DO jlay = 1, numlay delsord(jlay) = delsord0 !fractionation factors are constant delhsox(jlay) = delhsox0 !here, but can be varied by depth (jlay) dhs(jlay) = dhs0 dso(jlay) = dso0 hs(jlay) = hs0 + .1 so(jlay) = so40 o2(jlay) = .278 !approximate saturation @ 35S, 10 degrees C (Pilson, Chem.Ocean.) rhs(jlay) = ((dhs(jlay) * ref) / 1000.) + ref !isotope ratios calculated rso(jlay) = ((dso(jlay) * ref) / 1000.) + ref hs320 = 0.0 !boundary conditions hs340 = 0.0 so320 = so40 / (1. + rso(1)) so340 = (so40 * rso(1)) / (1. + rso(1)) hs32(jlay) = hs(jlay) / (1. + rhs(jlay)) !individual concentrations calculated hs34(jlay) = (hs(jlay) * rhs(jlay)) / (1. + rhs(jlay)) so32(jlay) = so40 / (1. + rso(jlay)) so34(jlay) = (so40 * rso(jlay)) / (1. + rso(jlay)) y(numspec * jlay - 4) = hs32(jlay) y(numspec * jlay - 3) = hs34(jlay) y(numspec * jlay - 2) = so32(jlay) y(numspec * jlay - 1) = so34(jlay) y(numspec * jlay) = o2(jlay) END DO c Values of incind are used in some versions to determine how CHECKSTEP c handles variables (see Walker) DO jrow = 1, numlay incind(numspec + jrow - 4) = 1 incind(numspec * jrow - 3) = 1 incind(numspec * jrow - 2) = 1 incind(numspec * jrow - 1) = 1 incind(numspec * jrow) = 1 END DO C output files are set up open(unit=10, file='AJS_sord_o2dhsoX24dat', status = 'unknown', & form='formatted') C ===================================================================== C The flow of the calculations is handled here. Subroutines are called and C the scope of the sensitivity test is determined by the initial conditions C and the "DO WHILE" statement Rates of oxidation and reduction, as well as C concentrations and fluxes are are calculated and converted to 32S and 34S C fluxes in REDOXCONVERT. The differential equations are set up in EQUATIONS. C These are solved in SLOPGAUSSND. The timestep is checked and adjusted in C CHECKSTEP. This continues until x = xend. Results are sent to PRINTER. C The initial step is 1 year. delx = 1. DO WHILE (sord0 < 500.) sord0 = sord0 * 10. c DO WHILE (o20 > 0.027) c o20 = o20 - 0.025 DO jlay = 1, numlay delsord(jlay) = delsord0 delhsox(jlay) = delhsox0 c o2(jlay) = o20 c y(numspec * jlay) = o2(jlay) END DO count = 0 x = xstart DO WHILE (x < xend) count = count + 1 CALL REDOXCONVERT (jlay,numlay,sord,g2,sord0,mpdepth, & sorddec,hsox,hsox0,hsoxdec,hs,y,numspec,so, & rsord,dso,delsord,ref,rhsox,dhs,delhsox,sord32, & sord34,hsox32,hsox34,hs32,hs34,so32,so34,o2) CALL EQUATIONS (jlay,numlay,numspec,hs320,hs340,so320, & so340,y,yp,transso,transhs,laysep,sord32,sord34,hsox32, & hsox34,o2,transo2,o20,hsox,sord,g2,sord0,mpdepth, & sorddec,hsox0,hsoxdec,hs,so,rsord,dso,delsord, & ref,rhsox,dhs,delhsox,hs32,hs34,so32,so34) CALL SLOPGAUSSND (jrow,numlay,numspec,nrow,delx,y,yp, & unk,jlay,hs320,hs340,so320,so340,transso, & transhs,laysep,sord32,sord34,hsox32,hsox34, & o2,transo2,hsox,sord,o20,g2,sord0,mpdepth, & sorddec,hsox0,hsoxdec,hs,so,rsord,dso,delsord, & ref,rhsox,dhs,delhsox,hs32,hs34,so32,so34) DO jrow = 1, nrow dely(jrow) = unk(jrow) END DO CALL CHECKSTEP (jrow,nrow,y,dely,x,delx,incind,mstep, & numlay,count,hs,so,dhs,dso,sord,hsox,sord0,hsox0,sorddec, & hsoxdec,delsord0,delhsox0,xend,mpdepth,hs32,hs34,so32,so34, & g2,o2,transo2,bioturb,maxinc,transso,transhs,laysep) END DO END DO STOP END PROGRAM Sulfur C*********************************************************************** SUBROUTINE EQUATIONS (jlay,numlay,numspec,hs320,hs340,so320,so340, & y,yp,transso,transhs,laysep,sord32,sord34,hsox32, & hsox34,o2,transo2,o20,hsox,sord,g2,sord0,mpdepth, & sorddec,hsox0,hsoxdec,hs,so,rsord,dso,delsord, & ref,rhsox,dhs,delhsox,hs32,hs34,so32,so34) Implicit none Integer :: jlay,jsc,numlay,numspec Double Precision :: hs320,hs340,so320,so340,o20, & g2,sord0,sorddec,hsox0,hsoxdec,ref & Double Precision, Dimension(75) :: y,yp Double Precision, Dimension(15) :: transso,transhs,laysep, & sord32,sord34,hsox32,hsox34,o2,transo2, & hsox,sord,mpdepth,hs,so,rsord,dso,rhsox, & delsord,dhs,delhsox,hs32,hs34,so32,so34 C ===================================================================== CALL REDOXCONVERT (jlay,numlay,sord,g2,sord0,mpdepth, & sorddec,hsox,hsox0,hsoxdec,hs,y,numspec,so, & rsord,dso,delsord,ref,rhsox,dhs,delhsox,sord32, & sord34,hsox32,hsox34,hs32,hs34,so32,so34,o2) C ===================================================================== C The differential equations that define the simulations. C equations for dissolved 32-sulfide. yp(1) = ((0. - y(1)) / laysep(1) + (y(6) - y(1)) / laysep(2)) & * transhs(1) + sord32(1) - hsox32(1) DO jlay = 2, numlay - 1 jsc = numspec * jlay - 4 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(jlay) & + (y(jsc + numspec) - y(jsc)) / laysep(jlay + 1)) & * transhs(jlay) + sord32(jlay) - hsox32(jlay) END DO jsc = numspec * numlay - 4 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(numlay)) & * transhs(numlay) + sord32(numlay) - hsox32(numlay) C ===================================================================== C equations for dissolved 34-sulfide yp(2) = ((0. - y(2)) / laysep(1) + (y(7) - y(2)) / laysep(2)) & * transhs(1) + sord34(1) - hsox34(1) DO jlay = 2, numlay - 1 jsc = numspec * jlay - 3 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(jlay) & + (y(jsc + numspec) - y(jsc)) / laysep(jlay + 1)) & * transhs(jlay) + sord34(jlay) - hsox34(jlay) END DO jsc = numspec * numlay - 3 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(numlay)) & * transhs(numlay) + sord34(numlay) - hsox34(numlay) C ===================================================================== C equations for dissolved 32-sulfate yp(3) = ((so320 - y(3)) / laysep(1) + (y(8) - y(3)) / laysep(2)) & * transso(1) - sord32(1) + hsox32(1) DO jlay = 2, numlay - 1 jsc = numspec * jlay - 2 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(jlay) & + (y(jsc + numspec) - y(jsc)) / laysep(jlay + 1)) & * transso(jlay) - sord32(jlay) + hsox32(jlay) END DO jsc = numspec * numlay - 2 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(numlay)) & * transso(numlay) - sord32(numlay) + hsox32(numlay) C ===================================================================== C equations for dissolved 34-sulfate yp(4) = ((so340 - y(4)) / laysep(1) + (y(9) - y(4)) / laysep(2)) & * transso(1) - sord34(1) + hsox34(1) DO jlay = 2, numlay - 1 jsc = numspec * jlay - 1 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(jlay) & + (y(jsc + numspec) - y(jsc)) / laysep(jlay + 1)) & * transso(jlay) - sord34(jlay) + hsox34(jlay) END DO jsc = numspec * numlay - 1 yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(numlay)) & * transso(numlay) - sord34(numlay) + hsox34(numlay) C ===================================================================== C equations for dissolved oxygen yp(5) = ((o20 - y(5)) / laysep(1) + (y(10) - y(5)) / laysep(2)) & * transo2(1) - (2. * hsox(1)) DO jlay = 2, numlay - 1 jsc = numspec * jlay yp(jsc) = ((y(jsc - numspec) - y(jsc)) / laysep(jlay) & + (y(jsc + numspec) - y(jsc)) / laysep(jlay + 1)) & * transo2(jlay) - (2. * hsox(jlay)) END DO jsc = numspec * numlay yp(jsc) = (y(jsc - numspec) - y(jsc)) / laysep(numlay) & * transo2(numlay) - (2. * hsox(numlay)) RETURN END SUBROUTINE EQUATIONS C*********************************************************************** SUBROUTINE SLOPGAUSSND (jrow,numlay,numspec,nrow,delx,y,yp, & unk,jlay,hs320,hs340,so320,so340,transso, & transhs,laysep,sord32,sord34,hsox32,hsox34, & o2,transo2,hsox,sord,o20,g2,sord0,mpdepth, & sorddec,hsox0,hsoxdec,hs,so,rsord,dso,delsord, & ref,rhsox,dhs,delhsox,hs32,hs34,so32,so34) Implicit none Integer :: jrow,jcol,jcll,jcul,jrns, & jrpl,jr,numlay,numspec,nrow,jlay Integer :: ncol Double Precision :: yinc,delx,diag,rsum, & coeff1,hs320,hs340,so320,so340,o20 & ,g2,sord0,sorddec,hsox0,hsoxdec,ref Double Precision, Dimension(15) :: transso,transhs, & laysep,sord32,sord34,hsox32,hsox34,o2, & transo2,hsox,sord,mpdepth,hs,so,rsord, & dso,delsord,rhsox,dhs,delhsox,hs32, & hs34,so32,so34 Double Precision, parameter :: dlny=.001 Double Precision, Dimension(75,76) :: sleq Double Precision, Dimension(75) :: y,yp,unk C Subroutine to calculate the coefficient matrix, SLEQ C for a nearly diagonal matrix, and to solve the simultaneous linear C equations. SLOPGAUSSND is a combination of SLOPERND and GAUSSND C (see below). C calculate the derivatives ncol = nrow + 1 DO jrow = 1, nrow sleq(jrow, ncol) = yp(jrow) END DO DO jcol = 1, nrow IF (y(jcol) == 0) THEN yinc = dlny ELSE yinc = y(jcol) * dlny END IF y(jcol) = y(jcol) + yinc CALL EQUATIONS (jlay,numlay,numspec,hs320,hs340,so320, & so340,y,yp,transso,transhs,laysep,sord32,sord34,hsox32, & hsox34,o2,transo2,o20,hsox,sord,g2,sord0,mpdepth, & sorddec,hsox0,hsoxdec,hs,so,rsord,dso,delsord, & ref,rhsox,dhs,delhsox,hs32,hs34,so32,so34) C Limits on JROW jcll = jcol - numspec jcul = jcol + numspec IF (jcll < 1) THEN jcll = 1 END IF IF (jcul > nrow) THEN jcul = nrow END IF C differentiate with respect to Y(jcol) DO jrow = jcll, jcul sleq(jrow, jcol) = -(yp(jrow) - sleq(jrow, ncol)) / yinc END DO C restore the original value y(jcol) = y(jcol) - yinc sleq(jcol, jcol) = sleq(jcol, jcol) + 1 / delx C extra term in diagonal elements END DO C ==================================================================== C GAUSSND solves a system of simultaneous linear C algebraic equations by Gaussian elimination and back substitution. C The number of equations (equal to the number of unknowns) is NROW. C The coefficients are in array SLEQ(NROW,NROW+1), where the last column C contains the constants on the right hand sides of the equations. C The answers are returned in the array UNK(NROW). C The array is nearly diagonal, with NUMSPEC nonzero elements C off the diagonal DO jrow = 1, nrow - 1 diag = sleq(jrow, jrow) jrns = jrow + numspec IF (jrns > nrow) THEN jrns = nrow END IF jrpl = jrow + 1 DO jcol = jrpl, jrns C divide by coefficient on the diagonal sleq(jrow, jcol) = sleq(jrow, jcol) / diag END DO sleq(jrow, ncol) = sleq(jrow, ncol) / diag sleq(jrow, jrow) = 1 DO jr = jrpl, jrns coeff1 = sleq(jr, jrow) C zeroes below the diagonal DO jcol = jrpl, jrns sleq(jr, jcol) = sleq(jr, jcol) - sleq(jrow, jcol) & * coeff1 END DO sleq(jr, ncol) = sleq(jr, ncol) - sleq(jrow, ncol) & * coeff1 sleq(jr, jrow) = 0 END DO END DO sleq(nrow, ncol) = sleq(nrow, ncol) / sleq(nrow, nrow) sleq(nrow, nrow) = 1 C ===================================================================== C calculate unknowns by back substitution unk(nrow) = sleq(nrow, ncol) DO jrow = nrow - 1, 1, -1 rsum = 0. jrns = jrow + numspec IF (jrns > nrow) THEN jrns = nrow END IF DO jcol = jrow + 1, jrns rsum = rsum + unk(jcol) * sleq(jrow, jcol) END DO unk(jrow) = sleq(jrow, ncol) - rsum END DO RETURN END SUBROUTINE SLOPGAUSSND C ********************************************************************** SUBROUTINE CHECKSTEP (jrow,nrow,y,dely,x,delx,incind,mstep,numlay, & count,hs,so,dhs,dso,sord,hsox,sord0,hsox0,sorddec,hsoxdec, & delsord0,delhsox0,xend,mpdepth,hs32,hs34,so32,so34,g2, & o2,transo2,bioturb,maxinc,transso,transhs,laysep) Implicit none Integer :: jrow,nrow,bigone,numlay,count Double Precision :: x,delx,mstep,maxinc, & biginc,relinc,sord0,hsox0,delsord0, & delhsox0,sorddec,hsoxdec,xend,g2, & bioturb Double Precision, Dimension(75) :: y,dely Double Precision, Dimension(15) :: hs,so,dhs,dso,sord, & incind,mpdepth,hs32,hs34,so32,so34, & hsox,o2,transo2,transso,transhs,laysep C Subroutine checks for the largest relative increment C ABS(dely/y); compares it with a specified value MAXINC, and adjusts C the step size if adjustment is needed. C If INCIND <> 1 then absolute increment is used. C In this case the largest permitted increment is INCIND. C If INCIND = 0 the variable is not tested. biginc = 0. DO jrow = 1, nrow C Find the largest relative increment or absolute increment, C depending on INCIND (INCIND is not used in this version.) IF (y(jrow) .EQ. 0.) THEN relinc = (ABS(dely(jrow)) / .005) * maxinc ELSE relinc = ABS(dely(jrow) / y(jrow)) END IF IF (biginc .le. relinc) THEN biginc = relinc END IF IF (biginc .EQ. relinc) THEN bigone = jrow END IF END DO IF ((biginc .lt. maxinc) .and. (biginc .gt. maxinc / 2.)) THEN CALL STEPPER (jrow,numlay,count,hs,so,dhs,dso,sord,hsox,x, & sord0,hsox0,sorddec,hsoxdec,delsord0,delhsox0,nrow,y, & dely,delx,xend,mpdepth,hs32,hs34,so32,so34,g2, & o2,transo2,biginc,maxinc,bioturb,mstep,transso, & transhs,laysep) !Increments Y and X and writes files ELSE IF (biginc .ge. maxinc) THEN delx = delx / 2. !reduce DELX and return without stepping forward ELSE IF (biginc .le. maxinc / 2.) THEN delx = delx * 1.4 CALL STEPPER (jrow,numlay,count,hs,so,dhs,dso,sord,hsox,x, & sord0,hsox0,sorddec,hsoxdec,delsord0,delhsox0,nrow,y, & dely,delx,xend,mpdepth,hs32,hs34,so32,so34,g2, & o2,transo2,biginc,maxinc,bioturb,mstep,transso, & transhs,laysep) !increase DELX and step forward END IF RETURN END SUBROUTINE CHECKSTEP C*********************************************************************** SUBROUTINE STEPPER (jrow,numlay,count,hs,so,dhs,dso,sord,hsox,x, & sord0,hsox0,sorddec,hsoxdec,delsord0,delhsox0,nrow,y, & dely,delx,xend,mpdepth,hs32,hs34,so32,so34,g2, & o2,transo2,biginc,maxinc,bioturb,mstep,transso, & transhs,laysep) Implicit none Integer :: jrow,nrow,numlay,count Double Precision :: x,delx,sord0,hsox0, & hsoxdec,sorddec,delhsox0,delsord0,xend, & g2,biginc,maxinc,bioturb,mstep Double Precision, Dimension(75) :: y,dely Double Precision, Dimension(15) :: hs,so,dhs,dso,sord, & hsox,mpdepth,hs32,hs34,so32,so34,o2, & transo2, transso,transhs,laysep C Subroutine steps forward in time Do jrow = 1, nrow y(jrow) = y(jrow) + dely(jrow) END DO IF (delx > mstep) THEN delx = mstep !put an upper limit on DELX END IF x = x + delx print*, count,sord0,sorddec,g2,hsox0,hsoxdec,delsord0, & delhsox0,x,delx IF (x > xend) THEN Call PRINTER (jrow,numlay,count,hs,so,dhs,dso,sord,hsox,x, & sord0,hsox0,sorddec,hsoxdec,delsord0,delhsox0,mpdepth, & hs32,hs34,so32,so34,y,g2,o2,transo2,bioturb) END IF RETURN END SUBROUTINE STEPPER C*********************************************************************** SUBROUTINE REDOXCONVERT (jlay,numlay,sord,g2,sord0,mpdepth, & sorddec,hsox,hsox0,hsoxdec,hs,y,numspec,so, & rsord,dso,delsord,ref,rhsox,dhs,delhsox,sord32, & sord34,hsox32,hsox34,hs32,hs34,so32,so34,o2) C Determines depth-dependent rates of oxidation and reduction, C and concentrations of HS and SO4. These are converted to concentrations C of 32HS,34HS,32SO4,34SO4,reductive and oxidative fluxes of 32S and C 34S, as well as delta34S-HS and delta34S-SO4 values. Implicit none Integer :: jlay,numlay,numspec Double Precision, Dimension(15) :: hsox,sord, & mpdepth,delhsox,delsord, & so,hs,rsord,rhsox,so32,dhs,dso, & so34,hs32,hs34,sord32,sord34, & hsox32,hsox34,o2,rhs,rso Double Precision, Dimension (75) :: y Double Precision :: ref,g2,sord0,hsox0, & delsord0,delhsox0,hsoxdec,sorddec DO jlay = 1, numlay hs32(jlay) = y(jlay * numspec - 4) hs34(jlay) = y(jlay * numspec - 3) so32(jlay) = y(jlay * numspec - 2) so34(jlay) = y(jlay * numspec - 1) o2(jlay) = y(jlay * numspec) hs(jlay) = hs32(jlay) + hs34(jlay) so(jlay) = so32(jlay) + so34(jlay) rhs(jlay) = hs34(jlay) / hs32(jlay) rso(jlay) = so34(jlay) / so32(jlay) dhs(jlay) = ((rhs(jlay) - ref) / ref) * 1000. dso(jlay) = ((rso(jlay) - ref) / ref) * 1000. sord(jlay) = g2 + sord0 * EXP(-mpdepth(jlay) / sorddec) !rates of reduction IF (so(jlay) < 3.) THEN sord(jlay) = sord(jlay) * so(jlay) / 3. !rates slowed at low SO4 END IF hsox(jlay) = hs(jlay) * o2(jlay) !rates of oxidation END DO C Converts delta notation to isotope ratios and then C to 32 and 34S fluxes produced during oxidation and reduction. DO jlay = 1, numlay rsord(jlay) = rso(jlay) / (((-1. * delsord(jlay)) / 1000.) + 1.) rhsox(jlay) = rhs(jlay) / (((-1. * delhsox(jlay)) / 1000.) + 1.) sord32(jlay) = sord(jlay) / (1. + rsord(jlay)) sord34(jlay) = (sord(jlay) * rsord(jlay)) / (1.+ rsord(jlay)) hsox32(jlay) = hsox(jlay) / (1. + rhsox(jlay)) hsox34(jlay) = (hsox(jlay) * rhsox(jlay)) / (1.+ rhsox(jlay)) END DO RETURN END SUBROUTINE REDOXCONVERT C*********************************************************************** SUBROUTINE PRINTER (jrow,numlay,count,hs,so,dhs,dso,sord,hsox,x, & sord0,hsox0,sorddec,hsoxdec,delsord0,delhsox0,mpdepth, & hs32,hs34,so32,so34,y,g2,o2,transo2,bioturb,o20) Implicit none Integer :: jrow,numlay,count Double Precision, Dimension(15) :: hs,so,dhs,dso,sord,hsox, & mpdepth,hs32,hs34,so32,so34,o2,transo2 Double Precision :: sord0,hsox0,sorddec,hsoxdec, & x,delsord0,delhsox0,g2,bioturb,o20 Double Precision, Dimension(75) :: y C Subroutine writes a file write(10,100) count, sord0, sorddec, hsox0, hsoxdec, delsord0, & delhsox0, o20 100 format (I8,3x,7F6.2) DO jrow = 1, numlay C write(10,'(30e15.6)')dhs(jrow),dso(jrow) c & ,hs34(jrow),so32(jrow),so34(jrow) c write(85,*)year,sum_antemp,sum_antempC3,sum_antempC4, c &sum_antempC4/sum_antemp write(10,'(105e15.6)')mpdepth(jrow),hs(jrow),so(jrow), & dhs(jrow),dso(jrow),sord(jrow),hsox(jrow),o2(jrow) END DO RETURN END SUBROUTINE PRINTER C***********************************************************************