[3443] | 1 | MODULE seddsr |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE seddsr *** |
---|
[10345] | 4 | !! Sediment : dissolution and reaction in pore water related |
---|
| 5 | !! related to organic matter |
---|
[3443] | 6 | !!===================================================================== |
---|
| 7 | !! * Modules used |
---|
| 8 | USE sed ! sediment global variable |
---|
[10345] | 9 | USE sed_oce |
---|
| 10 | USE sedini |
---|
| 11 | USE lib_mpp ! distribued memory computing library |
---|
| 12 | USE lib_fortran |
---|
[3443] | 13 | |
---|
[10345] | 14 | IMPLICIT NONE |
---|
| 15 | PRIVATE |
---|
| 16 | |
---|
[3443] | 17 | PUBLIC sed_dsr |
---|
| 18 | |
---|
| 19 | !! * Module variables |
---|
| 20 | |
---|
[10345] | 21 | REAL(wp) :: zadsnh4 |
---|
| 22 | REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density |
---|
| 23 | REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables |
---|
[3443] | 24 | |
---|
[10345] | 25 | |
---|
[5215] | 26 | !! $Id$ |
---|
[3443] | 27 | CONTAINS |
---|
| 28 | |
---|
[10345] | 29 | SUBROUTINE sed_dsr( kt, knt ) |
---|
[3443] | 30 | !!---------------------------------------------------------------------- |
---|
| 31 | !! *** ROUTINE sed_dsr *** |
---|
| 32 | !! |
---|
| 33 | !! ** Purpose : computes pore water dissolution and reaction |
---|
| 34 | !! |
---|
[10345] | 35 | !! ** Methode : Computation of the redox reactions in sediment. |
---|
| 36 | !! The main redox reactions are solved in sed_dsr whereas |
---|
| 37 | !! the secondary reactions are solved in sed_dsr_redoxb. |
---|
| 38 | !! A strand spliting approach is being used here (see |
---|
| 39 | !! sed_dsr_redoxb for more information). |
---|
[3443] | 40 | !! |
---|
| 41 | !! History : |
---|
| 42 | !! ! 98-08 (E. Maier-Reimer, Christoph Heinze ) Original code |
---|
| 43 | !! ! 04-10 (N. Emprin, M. Gehlen ) f90 |
---|
| 44 | !! ! 06-04 (C. Ethe) Re-organization |
---|
[10345] | 45 | !! ! 19-08 (O. Aumont) Debugging and improvement of the model. |
---|
| 46 | !! The original method is replaced by a |
---|
| 47 | !! Strand splitting method which deals |
---|
| 48 | !! well with stiff reactions. |
---|
[3443] | 49 | !!---------------------------------------------------------------------- |
---|
| 50 | !! Arguments |
---|
[10345] | 51 | INTEGER, INTENT(in) :: kt, knt ! number of iteration |
---|
[3443] | 52 | ! --- local variables |
---|
[10345] | 53 | INTEGER :: ji, jk, js, jw, jn ! dummy looop indices |
---|
[3443] | 54 | |
---|
[10345] | 55 | REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water |
---|
| 56 | REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite |
---|
| 57 | REAL(wp), DIMENSION(jpoce,jpksed) :: zkpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite |
---|
| 58 | REAL(wp), DIMENSION(jpoce) :: zsumtot |
---|
[3443] | 59 | REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat |
---|
[10345] | 60 | REAL(wp) :: zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc |
---|
| 61 | REAL(wp) :: zratio, zgamma, zbeta, zlimtmp, zundsat2 |
---|
[3443] | 62 | !! |
---|
| 63 | !!---------------------------------------------------------------------- |
---|
| 64 | |
---|
[10345] | 65 | IF( ln_timing ) CALL timing_start('sed_dsr') |
---|
| 66 | ! |
---|
| 67 | IF( kt == nitsed000 .AND. knt == 1 ) THEN |
---|
| 68 | IF (lwp) THEN |
---|
| 69 | WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' |
---|
| 70 | WRITE(numsed,*) ' ' |
---|
| 71 | ENDIF |
---|
[3443] | 72 | ENDIF |
---|
| 73 | |
---|
[10345] | 74 | ! Initializations |
---|
| 75 | !---------------------- |
---|
[3443] | 76 | |
---|
[10345] | 77 | zrearat1(:,:) = 0. ; zundsat(:,:) = 0. ; zkpoc(:,:) = 0. |
---|
| 78 | zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. |
---|
| 79 | zlimso4(:,:) = 0. ; zkpor(:,:) = 0. ; zrearat3(:,:) = 0. |
---|
| 80 | zkpos (:,:) = 0. |
---|
| 81 | zsumtot(:) = rtrn |
---|
| 82 | |
---|
| 83 | ALLOCATE( zvolc(jpoce, jpksed, jpsol) ) |
---|
| 84 | zvolc(:,:,:) = 0. |
---|
| 85 | zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) |
---|
[3443] | 86 | |
---|
[10345] | 87 | ! Inhibition terms for the different redox equations |
---|
| 88 | ! -------------------------------------------------- |
---|
| 89 | DO jk = 1, jpksed |
---|
| 90 | DO ji = 1, jpoce |
---|
| 91 | zkpoc(ji,jk) = reac_pocl |
---|
| 92 | zkpos(ji,jk) = reac_pocs |
---|
| 93 | zkpor(ji,jk) = reac_pocr |
---|
| 94 | END DO |
---|
| 95 | END DO |
---|
[3443] | 96 | |
---|
| 97 | ! Conversion of volume units |
---|
| 98 | !---------------------------- |
---|
| 99 | DO js = 1, jpsol |
---|
| 100 | DO jk = 1, jpksed |
---|
[10345] | 101 | DO ji = 1, jpoce |
---|
[3443] | 102 | zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / & |
---|
[10345] | 103 | & ( volw3d(ji,jk) * 1.e-3 ) |
---|
[3443] | 104 | ENDDO |
---|
| 105 | ENDDO |
---|
| 106 | ENDDO |
---|
| 107 | |
---|
| 108 | !---------------------------------------------------------- |
---|
[10345] | 109 | ! 5. Beginning of solid reaction |
---|
[3443] | 110 | !--------------------------------------------------------- |
---|
| 111 | |
---|
| 112 | ! Definition of reaction rates [rearat]=sans dim |
---|
| 113 | ! For jk=1 no reaction (pure water without solid) for each solid compo |
---|
[10345] | 114 | zrearat1(:,:) = 0. |
---|
| 115 | zrearat2(:,:) = 0. |
---|
| 116 | zrearat3(:,:) = 0. |
---|
[3443] | 117 | |
---|
[10345] | 118 | zundsat(:,:) = pwcp(:,:,jwoxy) |
---|
[3443] | 119 | |
---|
| 120 | DO jk = 2, jpksed |
---|
| 121 | DO ji = 1, jpoce |
---|
[10345] | 122 | zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) |
---|
| 123 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 124 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 125 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 126 | zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) |
---|
| 127 | zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) |
---|
| 128 | zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) |
---|
| 129 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 130 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 131 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 132 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 133 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 134 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
[3443] | 135 | ENDDO |
---|
| 136 | ENDDO |
---|
| 137 | |
---|
[10345] | 138 | ! left hand side of coefficient matrix |
---|
| 139 | ! DO jn = 1, 5 |
---|
| 140 | DO jk = 2, jpksed |
---|
| 141 | DO ji = 1, jpoce |
---|
| 142 | jflag1: DO jn = 1, 10 |
---|
| 143 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 144 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 145 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 146 | zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * ( zrearat1(ji,jk) & |
---|
| 147 | & + zrearat2(ji,jk) + zrearat3(ji,jk) ) |
---|
| 148 | zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) |
---|
| 149 | zundsat2 = zundsat(ji,jk) |
---|
| 150 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
| 151 | zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) |
---|
| 152 | zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) |
---|
| 153 | zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) |
---|
| 154 | zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) |
---|
| 155 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 156 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 157 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 158 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 159 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 160 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
| 161 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
| 162 | EXIT jflag1 |
---|
| 163 | ENDIF |
---|
| 164 | END DO jflag1 |
---|
| 165 | END DO |
---|
| 166 | END DO |
---|
[3443] | 167 | |
---|
| 168 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
[10345] | 169 | DO jk = 2, jpksed |
---|
[3443] | 170 | DO ji = 1, jpoce |
---|
[10345] | 171 | zreasat = zrearat1(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
| 172 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
| 173 | zreasat = zrearat2(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
| 174 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
| 175 | zreasat = zrearat3(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
| 176 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
[3443] | 177 | ENDDO |
---|
| 178 | ENDDO |
---|
| 179 | |
---|
| 180 | ! New pore water concentrations |
---|
[10345] | 181 | DO jk = 2, jpksed |
---|
[3443] | 182 | DO ji = 1, jpoce |
---|
| 183 | ! Acid Silicic |
---|
[10345] | 184 | pwcp(ji,jk,jwoxy) = zundsat(ji,jk) |
---|
| 185 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimo2(ji,jk) * zundsat(ji,jk) ! oxygen |
---|
[3443] | 186 | ! For DIC |
---|
| 187 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
[10345] | 188 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
[3443] | 189 | ! For Phosphate (in mol/l) |
---|
[10345] | 190 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
| 191 | ! For iron (in mol/l) |
---|
| 192 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
[3443] | 193 | ! For alkalinity |
---|
[10345] | 194 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) |
---|
| 195 | ! Ammonium |
---|
| 196 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
| 197 | ! Ligands |
---|
| 198 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) |
---|
[3443] | 199 | ENDDO |
---|
| 200 | ENDDO |
---|
| 201 | |
---|
| 202 | !-------------------------------------------------------------------- |
---|
| 203 | ! Begining POC denitrification and NO3- diffusion |
---|
| 204 | ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) |
---|
| 205 | !-------------------------------------------------------------------- |
---|
| 206 | |
---|
[10345] | 207 | zrearat1(:,:) = 0. |
---|
| 208 | zrearat2(:,:) = 0. |
---|
| 209 | zrearat3(:,:) = 0. |
---|
| 210 | |
---|
| 211 | zundsat(:,:) = pwcp(:,:,jwno3) |
---|
| 212 | |
---|
| 213 | DO jk = 2, jpksed |
---|
[3443] | 214 | DO ji = 1, jpoce |
---|
[10345] | 215 | zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) |
---|
| 216 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 217 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 218 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 219 | zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) |
---|
| 220 | zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) |
---|
| 221 | zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) |
---|
| 222 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 223 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 224 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 225 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 226 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 227 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
| 228 | END DO |
---|
| 229 | END DO |
---|
| 230 | |
---|
| 231 | ! DO jn = 1, 5 |
---|
[3443] | 232 | DO jk = 2, jpksed |
---|
| 233 | DO ji = 1, jpoce |
---|
[10345] | 234 | jflag2: DO jn = 1, 10 |
---|
| 235 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) |
---|
[3443] | 236 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
[10345] | 237 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 238 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 239 | zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * ( zrearat1(ji,jk) & |
---|
| 240 | & + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimtmp |
---|
| 241 | zgamma = - xksedno3 * pwcp(ji,jk,jwno3) |
---|
| 242 | zundsat2 = zundsat(ji,jk) |
---|
| 243 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
| 244 | zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) |
---|
| 245 | zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) |
---|
| 246 | zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) |
---|
| 247 | zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) |
---|
| 248 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 249 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 250 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 251 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 252 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 253 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
| 254 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
| 255 | EXIT jflag2 |
---|
| 256 | ENDIF |
---|
| 257 | END DO jflag2 |
---|
[3443] | 258 | END DO |
---|
| 259 | END DO |
---|
| 260 | |
---|
| 261 | |
---|
| 262 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
| 263 | DO jk = 2, jpksed |
---|
| 264 | DO ji = 1, jpoce |
---|
[10345] | 265 | zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
[3443] | 266 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
[10345] | 267 | zreasat = zrearat2(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
| 268 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
| 269 | zreasat = zrearat3(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
| 270 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
[3443] | 271 | ENDDO |
---|
| 272 | ENDDO |
---|
| 273 | |
---|
| 274 | ! New dissolved concentrations |
---|
[10345] | 275 | DO jk = 2, jpksed |
---|
[3443] | 276 | DO ji = 1, jpoce |
---|
| 277 | ! For nitrates |
---|
[10345] | 278 | pwcp(ji,jk,jwno3) = zundsat(ji,jk) |
---|
| 279 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimno3(ji,jk) * zundsat(ji,jk) |
---|
[3443] | 280 | ! For DIC |
---|
| 281 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
[10345] | 282 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
[3443] | 283 | ! For Phosphate (in mol/l) |
---|
| 284 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
[10345] | 285 | ! Ligands |
---|
| 286 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
| 287 | ! For iron (in mol/l) |
---|
| 288 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
[3443] | 289 | ! For alkalinity |
---|
[10345] | 290 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srDnit + srno3 * zadsnh4 - 2.* spo4r ) |
---|
| 291 | ! Ammonium |
---|
| 292 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
[3443] | 293 | ENDDO |
---|
| 294 | ENDDO |
---|
| 295 | |
---|
[10345] | 296 | !-------------------------------------------------------------------- |
---|
| 297 | ! Begining POC iron reduction |
---|
| 298 | ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) |
---|
| 299 | !-------------------------------------------------------------------- |
---|
[3443] | 300 | |
---|
[10345] | 301 | zrearat1(:,:) = 0. |
---|
| 302 | zrearat2(:,:) = 0. |
---|
| 303 | zrearat3(:,:) = 0. |
---|
[3443] | 304 | |
---|
[10345] | 305 | zundsat(:,:) = solcp(:,:,jsfeo) |
---|
[3443] | 306 | |
---|
[10345] | 307 | DO jk = 2, jpksed |
---|
| 308 | DO ji = 1, jpoce |
---|
| 309 | zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 310 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) |
---|
| 311 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 312 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 313 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 314 | zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) |
---|
| 315 | zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) |
---|
| 316 | zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) |
---|
| 317 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 318 | & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) |
---|
| 319 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 320 | & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) |
---|
| 321 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 322 | & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) |
---|
| 323 | END DO |
---|
| 324 | END DO |
---|
[3443] | 325 | |
---|
[10345] | 326 | ! DO jn = 1, 5 |
---|
| 327 | DO jk = 2, jpksed |
---|
[3443] | 328 | DO ji = 1, jpoce |
---|
[10345] | 329 | jflag3: DO jn = 1, 10 |
---|
| 330 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 331 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) |
---|
| 332 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 333 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 334 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 335 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) |
---|
| 336 | zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp |
---|
| 337 | zgamma = -xksedfeo * solcp(ji,jk,jsfeo) |
---|
| 338 | zundsat2 = zundsat(ji,jk) |
---|
| 339 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
| 340 | zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 341 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) |
---|
| 342 | zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) |
---|
| 343 | zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) |
---|
| 344 | zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) |
---|
| 345 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 346 | & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) |
---|
| 347 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 348 | & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) |
---|
| 349 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 350 | & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) |
---|
| 351 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/( MAX(0.,zundsat2)+rtrn)) < 1E-8 ) THEN |
---|
| 352 | EXIT jflag3 |
---|
| 353 | ENDIF |
---|
| 354 | END DO jflag3 |
---|
| 355 | END DO |
---|
| 356 | END DO |
---|
[3443] | 357 | |
---|
| 358 | |
---|
[10345] | 359 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
| 360 | DO jk = 2, jpksed |
---|
| 361 | DO ji = 1, jpoce |
---|
| 362 | zreasat = zrearat1(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
| 363 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
| 364 | zreasat = zrearat2(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
| 365 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
| 366 | zreasat = zrearat3(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
| 367 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
| 368 | END DO |
---|
| 369 | END DO |
---|
[3443] | 370 | |
---|
[10345] | 371 | ! New dissolved concentrations |
---|
| 372 | DO jk = 2, jpksed |
---|
[3443] | 373 | DO ji = 1, jpoce |
---|
[10345] | 374 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk) |
---|
| 375 | ! For FEOH |
---|
| 376 | solcp(ji,jk,jsfeo) = zundsat(ji,jk) |
---|
| 377 | ! For DIC |
---|
| 378 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
| 379 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
| 380 | ! For Phosphate (in mol/l) |
---|
| 381 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * ( spo4r + 4.0 * redfep ) |
---|
| 382 | ! Ligands |
---|
| 383 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
| 384 | ! For iron (in mol/l) |
---|
| 385 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
| 386 | ! For alkalinity |
---|
| 387 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + 8.0 * zreasat |
---|
| 388 | ! Ammonium |
---|
| 389 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
| 390 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 |
---|
[3443] | 391 | ENDDO |
---|
| 392 | ENDDO |
---|
| 393 | |
---|
[10345] | 394 | !-------------------------------------------------------------------- |
---|
| 395 | ! Begining POC denitrification and NO3- diffusion |
---|
| 396 | ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) |
---|
| 397 | !-------------------------------------------------------------------- |
---|
[3443] | 398 | |
---|
[10345] | 399 | zrearat1(:,:) = 0. |
---|
| 400 | zrearat2(:,:) = 0. |
---|
| 401 | zrearat3(:,:) = 0. |
---|
[3443] | 402 | |
---|
[10345] | 403 | zundsat(:,:) = pwcp(:,:,jwso4) |
---|
[3443] | 404 | |
---|
[10345] | 405 | DO jk = 2, jpksed |
---|
[3443] | 406 | DO ji = 1, jpoce |
---|
[10345] | 407 | zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 408 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
| 409 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) |
---|
| 410 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 411 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 412 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 413 | zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) |
---|
| 414 | zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) |
---|
| 415 | zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) |
---|
| 416 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 417 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 418 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 419 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 420 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 421 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
| 422 | END DO |
---|
| 423 | END DO |
---|
| 424 | ! |
---|
| 425 | ! DO jn = 1, 5 |
---|
[3443] | 426 | DO jk = 2, jpksed |
---|
| 427 | DO ji = 1, jpoce |
---|
[10345] | 428 | jflag4: DO jn = 1, 10 |
---|
| 429 | zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 430 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
| 431 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) |
---|
| 432 | zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) |
---|
| 433 | zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) |
---|
| 434 | zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) |
---|
| 435 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) |
---|
| 436 | zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp |
---|
| 437 | zgamma = - xksedso4 * pwcp(ji,jk,jwso4) |
---|
| 438 | zundsat2 = zundsat(ji,jk) |
---|
| 439 | zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 |
---|
| 440 | zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & |
---|
| 441 | & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & |
---|
| 442 | & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) |
---|
| 443 | zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) |
---|
| 444 | zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) |
---|
| 445 | zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) |
---|
| 446 | zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & |
---|
| 447 | & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) |
---|
| 448 | zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & |
---|
| 449 | & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) |
---|
| 450 | zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & |
---|
| 451 | & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) |
---|
| 452 | IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN |
---|
| 453 | EXIT jflag4 |
---|
| 454 | ENDIF |
---|
| 455 | END DO jflag4 |
---|
[3443] | 456 | END DO |
---|
| 457 | END DO |
---|
| 458 | |
---|
[10345] | 459 | ! New solid concentration values (jk=2 to jksed) for each couple |
---|
[3443] | 460 | DO jk = 2, jpksed |
---|
| 461 | DO ji = 1, jpoce |
---|
[10345] | 462 | zreasat = zrearat1(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) |
---|
| 463 | solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat |
---|
| 464 | zreasat = zrearat2(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) |
---|
| 465 | solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat |
---|
| 466 | zreasat = zrearat3(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) |
---|
| 467 | solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat |
---|
[3443] | 468 | ENDDO |
---|
| 469 | ENDDO |
---|
[10345] | 470 | ! |
---|
[3443] | 471 | ! New dissolved concentrations |
---|
[10345] | 472 | DO jk = 2, jpksed |
---|
[3443] | 473 | DO ji = 1, jpoce |
---|
[10345] | 474 | ! For sulfur |
---|
| 475 | pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) - ( zundsat(ji,jk) - pwcp(ji,jk,jwso4) ) |
---|
| 476 | pwcp(ji,jk,jwso4) = zundsat(ji,jk) |
---|
| 477 | zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimso4(ji,jk) * zundsat(ji,jk) |
---|
[3443] | 478 | ! For DIC |
---|
| 479 | pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat |
---|
[10345] | 480 | zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 |
---|
| 481 | ! For Phosphate (in mol/l) |
---|
| 482 | pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r |
---|
| 483 | ! Ligands |
---|
| 484 | sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat |
---|
| 485 | ! For iron (in mol/l) |
---|
| 486 | pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat |
---|
[3443] | 487 | ! For alkalinity |
---|
[10345] | 488 | pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + zreasat |
---|
| 489 | ! Ammonium |
---|
| 490 | pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 |
---|
[3443] | 491 | ENDDO |
---|
| 492 | ENDDO |
---|
| 493 | |
---|
[10345] | 494 | ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for |
---|
| 495 | ! There are two options here: A simple time splitting scheme and a modified |
---|
| 496 | ! Patankar scheme |
---|
| 497 | ! ------------------------------------------------------------------------------ |
---|
[3443] | 498 | |
---|
[10345] | 499 | call sed_dsr_redoxb |
---|
[3443] | 500 | |
---|
[10345] | 501 | ! -------------------------------------------------------------- |
---|
| 502 | ! 4/ Computation of the bioturbation coefficient |
---|
| 503 | ! This parameterization is taken from Archer et al. (2002) |
---|
| 504 | ! -------------------------------------------------------------- |
---|
[3443] | 505 | |
---|
[10345] | 506 | DO ji = 1, jpoce |
---|
| 507 | db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) |
---|
| 508 | END DO |
---|
[3443] | 509 | |
---|
[10345] | 510 | ! ------------------------------------------------------ |
---|
| 511 | ! Vertical variations of the bioturbation coefficient |
---|
| 512 | ! ------------------------------------------------------ |
---|
| 513 | IF (ln_btbz) THEN |
---|
[3443] | 514 | DO ji = 1, jpoce |
---|
[10345] | 515 | db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) |
---|
[3443] | 516 | END DO |
---|
[10345] | 517 | ELSE |
---|
| 518 | DO jk = 1, jpksed |
---|
| 519 | IF (profsedw(jk) > dbtbzsc) THEN |
---|
| 520 | db(:,jk) = 0.0 |
---|
| 521 | ENDIF |
---|
[3443] | 522 | END DO |
---|
[10345] | 523 | ENDIF |
---|
| 524 | |
---|
| 525 | IF (ln_irrig) THEN |
---|
| 526 | DO jk = 1, jpksed |
---|
| 527 | DO ji = 1, jpoce |
---|
| 528 | irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy) & |
---|
| 529 | & / (pwcp(ji,1,jwoxy) + 20.E-6) |
---|
| 530 | irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) |
---|
| 531 | END DO |
---|
[3443] | 532 | END DO |
---|
[10345] | 533 | ELSE |
---|
| 534 | irrig(:,:) = 0.0 |
---|
| 535 | ENDIF |
---|
| 536 | |
---|
| 537 | DEALLOCATE( zvolc ) |
---|
| 538 | |
---|
| 539 | IF( ln_timing ) CALL timing_stop('sed_dsr') |
---|
| 540 | ! |
---|
| 541 | END SUBROUTINE sed_dsr |
---|
| 542 | |
---|
| 543 | SUBROUTINE sed_dsr_redoxb |
---|
| 544 | !!---------------------------------------------------------------------- |
---|
| 545 | !! *** ROUTINE sed_dsr_redox *** |
---|
| 546 | !! |
---|
| 547 | !! ** Purpose : computes secondary redox reactions |
---|
| 548 | !! |
---|
| 549 | !! ** Methode : It uses Strand splitter algorithm proposed by |
---|
| 550 | !! Nguyen et al. (2009) and modified by Wang et al. (2018) |
---|
| 551 | !! Basically, each equation is solved analytically when |
---|
| 552 | !! feasible, otherwise numerically at t+1/2. Then |
---|
| 553 | !! the last equation is solved at t+1. The other equations |
---|
| 554 | !! are then solved at t+1 starting in the reverse order. |
---|
| 555 | !! Ideally, it's better to start from the fastest reaction |
---|
| 556 | !! to the slowest and then reverse the order to finish up |
---|
| 557 | !! with the fastest one. But random order works well also. |
---|
| 558 | !! The scheme is second order, positive and mass |
---|
| 559 | !! conserving. It works well for stiff systems. |
---|
| 560 | !! |
---|
| 561 | !! History : |
---|
| 562 | !! ! 18-08 (O. Aumont) Original code |
---|
| 563 | !!---------------------------------------------------------------------- |
---|
| 564 | !! Arguments |
---|
| 565 | ! --- local variables |
---|
| 566 | INTEGER :: ji, jk, jn ! dummy looop indices |
---|
| 567 | |
---|
| 568 | REAL, DIMENSION(6) :: zsedtrn, zsedtra |
---|
| 569 | REAL(wp) :: zalpha, zbeta, zgamma, zdelta, zepsi, zsedfer |
---|
| 570 | !! |
---|
| 571 | !!---------------------------------------------------------------------- |
---|
| 572 | |
---|
| 573 | IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') |
---|
| 574 | |
---|
[3443] | 575 | DO ji = 1, jpoce |
---|
[10345] | 576 | DO jk = 2, jpksed |
---|
| 577 | zsedtrn(1) = pwcp(ji,jk,jwoxy) |
---|
| 578 | zsedtrn(2) = MAX(0., pwcp(ji,jk,jwh2s) ) |
---|
| 579 | zsedtrn(3) = pwcp(ji,jk,jwnh4) |
---|
| 580 | zsedtrn(4) = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) |
---|
| 581 | zsedfer = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) |
---|
| 582 | zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) |
---|
| 583 | zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) |
---|
| 584 | zsedtra(:) = zsedtrn(:) |
---|
| 585 | |
---|
| 586 | ! First pass of the scheme. At the end, it is 1st order |
---|
| 587 | ! ----------------------------------------------------- |
---|
| 588 | ! Fe + O2 |
---|
| 589 | zalpha = zsedtra(1) - 0.25 * zsedtra(4) |
---|
| 590 | zbeta = zsedtra(4) + zsedtra(5) |
---|
| 591 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
| 592 | zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) |
---|
| 593 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 594 | & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) |
---|
| 595 | zsedtra(1) = zalpha + 0.25 * zsedtra(4) |
---|
| 596 | zsedtra(5) = zbeta - zsedtra(4) |
---|
| 597 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
| 598 | pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) |
---|
| 599 | ! H2S + O2 |
---|
| 600 | zalpha = zsedtra(1) - 2.0 * zsedtra(2) |
---|
| 601 | zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
| 602 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) |
---|
| 603 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 604 | & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) |
---|
| 605 | zsedtra(1) = zalpha + 2.0 * zsedtra(2) |
---|
| 606 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) |
---|
| 607 | pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) |
---|
| 608 | ! NH4 + O2 |
---|
| 609 | zalpha = zsedtra(1) - 2.0 * zsedtra(3) |
---|
| 610 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) |
---|
| 611 | zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 612 | & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 613 | zsedtra(1) = zalpha + 2.0 * zsedtra(3) |
---|
| 614 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) |
---|
| 615 | ! FeS - O2 |
---|
| 616 | zalpha = zsedtra(1) - 2.0 * zsedtra(6) |
---|
| 617 | zbeta = zsedtra(4) + zsedtra(6) |
---|
| 618 | zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) |
---|
| 619 | zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 620 | & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 621 | zsedtra(1) = zalpha + 2.0 * zsedtra(6) |
---|
| 622 | zsedtra(4) = zbeta - zsedtra(6) |
---|
| 623 | pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) |
---|
| 624 | ! ! Fe - H2S |
---|
| 625 | zalpha = zsedtra(2) - zsedtra(4) |
---|
| 626 | zbeta = zsedtra(4) + zsedtra(6) |
---|
| 627 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
| 628 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 629 | & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 630 | zsedtra(2) = zalpha + zsedtra(4) |
---|
| 631 | zsedtra(6) = zbeta - zsedtra(4) |
---|
| 632 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
| 633 | ! FEOH + H2S |
---|
| 634 | zalpha = zsedtra(5) - 2.0 * zsedtra(2) |
---|
| 635 | zbeta = zsedtra(5) + zsedtra(4) |
---|
| 636 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
| 637 | zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
| 638 | zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) |
---|
| 639 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & |
---|
| 640 | & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn ) |
---|
| 641 | zsedtra(5) = zalpha + 2.0 * zsedtra(2) |
---|
| 642 | zsedtra(4) = zbeta - zsedtra(5) |
---|
| 643 | pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) |
---|
| 644 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
| 645 | pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) |
---|
| 646 | ! Fe - H2S |
---|
| 647 | zalpha = zsedtra(2) - zsedtra(4) |
---|
| 648 | zbeta = zsedtra(4) + zsedtra(6) |
---|
| 649 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
| 650 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 651 | & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 652 | zsedtra(2) = zalpha + zsedtra(4) |
---|
| 653 | zsedtra(6) = zbeta - zsedtra(4) |
---|
| 654 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
| 655 | ! FeS - O2 |
---|
| 656 | zalpha = zsedtra(1) - 2.0 * zsedtra(6) |
---|
| 657 | zbeta = zsedtra(4) + zsedtra(6) |
---|
| 658 | zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) |
---|
| 659 | zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 660 | & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 661 | zsedtra(1) = zalpha + 2.0 * zsedtra(6) |
---|
| 662 | zsedtra(4) = zbeta - zsedtra(6) |
---|
| 663 | pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) |
---|
| 664 | ! NH4 + O2 |
---|
| 665 | zalpha = zsedtra(1) - 2.0 * zsedtra(3) |
---|
| 666 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) |
---|
| 667 | zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 668 | & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) |
---|
| 669 | zsedtra(1) = zalpha + 2.0 * zsedtra(3) |
---|
| 670 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) |
---|
| 671 | ! H2S + O2 |
---|
| 672 | zalpha = zsedtra(1) - 2.0 * zsedtra(2) |
---|
| 673 | zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) |
---|
| 674 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) |
---|
| 675 | zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 676 | & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) |
---|
| 677 | zsedtra(1) = zalpha + 2.0 * zsedtra(2) |
---|
| 678 | pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) |
---|
| 679 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) |
---|
| 680 | ! Fe + O2 |
---|
| 681 | zalpha = zsedtra(1) - 0.25 * zsedtra(4) |
---|
| 682 | zbeta = zsedtra(4) + zsedtra(5) |
---|
| 683 | zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) |
---|
| 684 | zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) |
---|
| 685 | zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & |
---|
| 686 | & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) |
---|
| 687 | zsedtra(1) = zalpha + 0.25 * zsedtra(4) |
---|
| 688 | zsedtra(5) = zbeta - zsedtra(4) |
---|
| 689 | pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) |
---|
| 690 | pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) |
---|
| 691 | pwcp(ji,jk,jwoxy) = zsedtra(1) |
---|
| 692 | pwcp(ji,jk,jwh2s) = zsedtra(2) |
---|
| 693 | pwcp(ji,jk,jwnh4) = zsedtra(3) |
---|
| 694 | pwcp(ji,jk,jwfe2) = zsedtra(4) + sedligand(ji,jk) + zsedfer |
---|
| 695 | pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) |
---|
| 696 | solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) |
---|
| 697 | solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) |
---|
[3443] | 698 | END DO |
---|
[10345] | 699 | END DO |
---|
[3443] | 700 | |
---|
[10345] | 701 | IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') |
---|
| 702 | |
---|
| 703 | END SUBROUTINE sed_dsr_redoxb |
---|
| 704 | |
---|
[3443] | 705 | END MODULE seddsr |
---|