MODULE seddsr !!====================================================================== !! *** MODULE seddsr *** !! Sediment : dissolution and reaction in pore water related !! related to organic matter !!===================================================================== !! * Modules used USE sed ! sediment global variable USE sed_oce USE sedini USE lib_mpp ! distribued memory computing library USE lib_fortran IMPLICIT NONE PRIVATE PUBLIC sed_dsr !! * Module variables REAL(wp) :: zadsnh4 REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables !! $Id$ CONTAINS SUBROUTINE sed_dsr( kt, knt ) !!---------------------------------------------------------------------- !! *** ROUTINE sed_dsr *** !! !! ** Purpose : computes pore water dissolution and reaction !! !! ** Methode : Computation of the redox reactions in sediment. !! The main redox reactions are solved in sed_dsr whereas !! the secondary reactions are solved in sed_dsr_redoxb. !! A strand spliting approach is being used here (see !! sed_dsr_redoxb for more information). !! !! History : !! ! 98-08 (E. Maier-Reimer, Christoph Heinze ) Original code !! ! 04-10 (N. Emprin, M. Gehlen ) f90 !! ! 06-04 (C. Ethe) Re-organization !! ! 19-08 (O. Aumont) Debugging and improvement of the model. !! The original method is replaced by a !! Strand splitting method which deals !! well with stiff reactions. !!---------------------------------------------------------------------- !! Arguments INTEGER, INTENT(in) :: kt, knt ! number of iteration ! --- local variables INTEGER :: ji, jk, js, jw, jn ! dummy looop indices REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite REAL(wp), DIMENSION(jpoce,jpksed) :: zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite REAL(wp), DIMENSION(jpoce) :: zsumtot REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat REAL(wp) :: zgamma, zbeta, zlimtmp !! !!---------------------------------------------------------------------- IF( ln_timing ) CALL timing_start('sed_dsr') ! IF( kt == nitsed000 .AND. knt == 1 ) THEN IF (lwp) THEN WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' WRITE(numsed,*) ' ' ENDIF ENDIF ! Initializations !---------------------- zrearat1(:,:) = 0. ; zundsat(:,:) = 0. zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. zlimso4(:,:) = 0. ; zrearat3(:,:) = 0. zsumtot(:) = rtrn ALLOCATE( zvolc(jpoce, jpksed, jpsol) ) zvolc(:,:,:) = 0. zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) ! Conversion of volume units !---------------------------- DO js = 1, jpsol DO jk = 1, jpksed DO ji = 1, jpoce zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / & & ( volw3d(ji,jk) * 1.e-3 ) ENDDO ENDDO ENDDO !---------------------------------------------------------- ! 5. Beginning of solid reaction !--------------------------------------------------------- ! Definition of reaction rates [rearat]=sans dim ! For jk=1 no reaction (pure water without solid) for each solid compo zrearat1(:,:) = 0. zrearat2(:,:) = 0. zrearat3(:,:) = 0. zundsat(:,:) = MAX( pwcp(:,:,jwoxy) - rtrn, 0. ) DO jk = 2, jpksed DO ji = 1, jpoce zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) zrearat1(ji,jk) = ( reac_pocl * dtsed2 * zsolid1 ) / ( 1. + reac_pocl * dtsed2 ) zrearat2(ji,jk) = ( reac_pocs * dtsed2 * zsolid2 ) / ( 1. + reac_pocs * dtsed2 ) zrearat3(ji,jk) = ( reac_pocr * dtsed2 * zsolid3 ) / ( 1. + reac_pocr * dtsed2 ) solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zrearat1(ji,jk) / zvolc(ji,jk,jspoc) solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zrearat2(ji,jk) / zvolc(ji,jk,jspos) solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zrearat3(ji,jk) / zvolc(ji,jk,jspor) zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ! For DIC pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 ! For alkalinity pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) ! For Phosphate (in mol/l) pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r ! For iron (in mol/l) pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat ! Ammonium pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 ! Ligands sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat ENDDO ENDDO ! left hand side of coefficient matrix DO jk = 2, jpksed DO ji = 1, jpoce zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * zreasat zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 zlimo2(ji,jk) = zundsat(ji,jk) / ( zundsat(ji,jk) + xksedo2 ) ! Oxygen pwcp(ji,jk,jwoxy) = zundsat(ji,jk) + MIN( pwcp(ji,jk,jwoxy), rtrn ) zreasat = zreasat * zlimo2(ji,jk) ! Ligands sedligand(ji,jk) = sedligand(ji,jk) - reac_ligc * sedligand(ji,jk) ENDDO ENDDO !-------------------------------------------------------------------- ! Begining POC denitrification and NO3- diffusion ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) !-------------------------------------------------------------------- zundsat(:,:) = MAX(0., pwcp(:,:,jwno3) - rtrn) DO jk = 2, jpksed DO ji = 1, jpoce zlimtmp = ( 1.0 - zlimo2(ji,jk) ) zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * zreasat * zlimtmp zgamma = - xksedno3 * pwcp(ji,jk,jwno3) zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 zlimno3(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedno3 ) ! For nitrates pwcp(ji,jk,jwno3) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwno3), rtrn) zreasat = zreasat * zlimno3(ji,jk) ! For alkalinity pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * srDnit ENDDO ENDDO !-------------------------------------------------------------------- ! Begining POC iron reduction ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) !-------------------------------------------------------------------- zundsat(:,:) = MAX(0., solcp(:,:,jsfeo) - rtrn) DO jk = 2, jpksed DO ji = 1, jpoce zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp zgamma = -xksedfeo * solcp(ji,jk,jsfeo) zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 zlimfeo(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedfeo ) zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) ! For FEOH solcp(ji,jk,jsfeo) = zundsat(ji,jk) + MIN(solcp(ji,jk,jsfeo), rtrn) ! For Phosphate (in mol/l) pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * 4.0 * redfep ! For alkalinity pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 8.0 * zreasat ! Iron pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 ENDDO ENDDO ! !-------------------------------------------------------------------- ! ! Begining POC denitrification and NO3- diffusion ! ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) ! !-------------------------------------------------------------------- ! zundsat(:,:) = MAX(0., pwcp(:,:,jwso4) - rtrn) DO jk = 2, jpksed DO ji = 1, jpoce zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) - zlimfeo(ji,jk) zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp zgamma = - xksedso4 * pwcp(ji,jk,jwso4) zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 zlimso4(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedso4 ) ! For sulfur pwcp(ji,jk,jwso4) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwso4), rtrn) zreasat = zreasat * zlimso4(ji,jk) pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) + 0.5 * zreasat ! For alkalinity pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat ENDDO ENDDO ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for ! There are two options here: A simple time splitting scheme and a modified ! Patankar scheme ! ------------------------------------------------------------------------------ call sed_dsr_redoxb ! -------------------------------------------------------------- ! 4/ Computation of the bioturbation coefficient ! This parameterization is taken from Archer et al. (2002) ! -------------------------------------------------------------- DO ji = 1, jpoce db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) END DO ! ------------------------------------------------------ ! Vertical variations of the bioturbation coefficient ! ------------------------------------------------------ IF (ln_btbz) THEN DO ji = 1, jpoce db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) END DO ELSE DO jk = 1, jpksed IF (profsedw(jk) > dbtbzsc) THEN db(:,jk) = 0.0 ENDIF END DO ENDIF IF (ln_irrig) THEN DO jk = 1, jpksed DO ji = 1, jpoce irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy) & & / (pwcp(ji,1,jwoxy) + 20.E-6) irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) END DO END DO ELSE irrig(:,:) = 0.0 ENDIF DEALLOCATE( zvolc ) IF( ln_timing ) CALL timing_stop('sed_dsr') ! END SUBROUTINE sed_dsr SUBROUTINE sed_dsr_redoxb !!---------------------------------------------------------------------- !! *** ROUTINE sed_dsr_redox *** !! !! ** Purpose : computes secondary redox reactions !! !! ** Methode : It uses Strand splitter algorithm proposed by !! Nguyen et al. (2009) and modified by Wang et al. (2018) !! Basically, each equation is solved analytically when !! feasible, otherwise numerically at t+1/2. Then !! the last equation is solved at t+1. The other equations !! are then solved at t+1 starting in the reverse order. !! Ideally, it's better to start from the fastest reaction !! to the slowest and then reverse the order to finish up !! with the fastest one. But random order works well also. !! The scheme is second order, positive and mass !! conserving. It works well for stiff systems. !! !! History : !! ! 18-08 (O. Aumont) Original code !!---------------------------------------------------------------------- !! Arguments ! --- local variables INTEGER :: ji, jk, jn, jw ! dummy looop indices REAL, DIMENSION(6) :: zsedtrn, zsedtra REAL(wp) :: zalpha, zbeta, zgamma, zdelta, zepsi, zsedfer !! !!---------------------------------------------------------------------- IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') DO ji = 1, jpoce DO jk = 2, jpksed zsedtrn(1) = pwcp(ji,jk,jwoxy) zsedtrn(2) = MAX(0., pwcp(ji,jk,jwh2s) ) zsedtrn(3) = pwcp(ji,jk,jwnh4) zsedtrn(4) = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) zsedfer = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) zsedtra(:) = MAX(0., zsedtrn(:) - rtrn ) ! First pass of the scheme. At the end, it is 1st order ! ----------------------------------------------------- ! Fe + O2 zalpha = zsedtra(1) - 0.25 * zsedtra(4) zbeta = zsedtra(4) + zsedtra(5) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) IF ( zalpha == 0. ) THEN zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) ELSE zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) ENDIF zsedtra(1) = zalpha + 0.25 * zsedtra(4) zsedtra(5) = zbeta - zsedtra(4) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) ! H2S + O2 zalpha = zsedtra(1) - 2.0 * zsedtra(2) zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) IF ( zalpha == 0. ) THEN zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) ELSE zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(2) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) ! NH4 + O2 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) /zadsnh4 IF ( zalpha == 0. ) THEN zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) ELSE zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 ! FeS - O2 zalpha = zsedtra(1) - 2.0 * zsedtra(6) zbeta = zsedtra(4) + zsedtra(6) zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) IF ( zalpha == 0. ) THEN zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2.0 ) ELSE zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(6) zsedtra(4) = zbeta - zsedtra(6) pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) ! ! Fe - H2S zalpha = zsedtra(2) - zsedtra(4) zbeta = zsedtra(4) + zsedtra(6) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) IF ( zalpha == 0. ) THEN zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) ELSE zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(2) = zalpha + zsedtra(4) zsedtra(6) = zbeta - zsedtra(4) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) ! FEOH + H2S zalpha = zsedtra(5) - 2.0 * zsedtra(2) zbeta = zsedtra(5) + zsedtra(4) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) IF ( zalpha == 0. ) THEN zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_feh2s * dtsed2 ) ELSE zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) ENDIF zsedtra(5) = zalpha + 2.0 * zsedtra(2) zsedtra(4) = zbeta - zsedtra(5) pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) ! Fe - H2S zalpha = zsedtra(2) - zsedtra(4) zbeta = zsedtra(4) + zsedtra(6) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) IF ( zalpha == 0. ) THEN zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) ELSE zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(2) = zalpha + zsedtra(4) zsedtra(6) = zbeta - zsedtra(4) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) ! FeS - O2 zalpha = zsedtra(1) - 2.0 * zsedtra(6) zbeta = zsedtra(4) + zsedtra(6) zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) IF (zalpha == 0.) THEN zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2. ) ELSE zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(6) zsedtra(4) = zbeta - zsedtra(6) pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) ! NH4 + O2 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) / zadsnh4 IF ( zalpha == 0. ) THEN zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) ELSE zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 ! H2S + O2 zalpha = zsedtra(1) - 2.0 * zsedtra(2) zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) IF ( zalpha == 0. ) THEN zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) ELSE zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) ENDIF zsedtra(1) = zalpha + 2.0 * zsedtra(2) pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) ! Fe + O2 zalpha = zsedtra(1) - 0.25 * zsedtra(4) zbeta = zsedtra(4) + zsedtra(5) zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) IF ( zalpha == 0. ) THEN zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) ELSE zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) ENDIF zsedtra(1) = zalpha + 0.25 * zsedtra(4) zsedtra(5) = zbeta - zsedtra(4) pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) ! Update the concentrations after the secondary reactions zsedtra(:) = zsedtra(:) + MIN( zsedtrn(:), rtrn ) pwcp(ji,jk,jwoxy) = zsedtra(1) pwcp(ji,jk,jwh2s) = zsedtra(2) pwcp(ji,jk,jwnh4) = zsedtra(3) pwcp(ji,jk,jwfe2) = zsedtra(4) + sedligand(ji,jk) + zsedfer pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) END DO END DO IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') END SUBROUTINE sed_dsr_redoxb END MODULE seddsr