- Timestamp:
- 2021-02-03T16:03:34+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsr.F90
r10362 r14385 55 55 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water 56 56 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite 57 REAL(wp), DIMENSION(jpoce,jpksed) :: z kpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite57 REAL(wp), DIMENSION(jpoce,jpksed) :: zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 58 58 REAL(wp), DIMENSION(jpoce) :: zsumtot 59 59 REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat 60 REAL(wp) :: zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 61 REAL(wp) :: zratio, zgamma, zbeta, zlimtmp, zundsat2 60 REAL(wp) :: zgamma, zbeta, zlimtmp 62 61 !! 63 62 !!---------------------------------------------------------------------- … … 75 74 !---------------------- 76 75 77 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. ; zkpoc(:,:) = 0. 78 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 79 zlimso4(:,:) = 0. ; zkpor(:,:) = 0. ; zrearat3(:,:) = 0. 80 zkpos (:,:) = 0. 76 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. 77 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 78 zlimso4(:,:) = 0. ; zrearat3(:,:) = 0. 81 79 zsumtot(:) = rtrn 82 80 … … 84 82 zvolc(:,:,:) = 0. 85 83 zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 86 87 ! Inhibition terms for the different redox equations88 ! --------------------------------------------------89 DO jk = 1, jpksed90 DO ji = 1, jpoce91 zkpoc(ji,jk) = reac_pocl92 zkpos(ji,jk) = reac_pocs93 zkpor(ji,jk) = reac_pocr94 END DO95 END DO96 84 97 85 ! Conversion of volume units … … 116 104 zrearat3(:,:) = 0. 117 105 118 zundsat(:,:) = pwcp(:,:,jwoxy)106 zundsat(:,:) = MAX( pwcp(:,:,jwoxy) - rtrn, 0. ) 119 107 120 108 DO jk = 2, jpksed 121 109 DO ji = 1, jpoce 122 zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 )123 110 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 124 111 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 125 112 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 ) 135 ENDDO 136 ENDDO 137 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 167 168 ! New solid concentration values (jk=2 to jksed) for each couple 169 DO jk = 2, jpksed 170 DO ji = 1, jpoce 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 177 ENDDO 178 ENDDO 179 180 ! New pore water concentrations 181 DO jk = 2, jpksed 182 DO ji = 1, jpoce 183 ! Acid Silicic 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 113 zrearat1(ji,jk) = ( reac_pocl * dtsed2 * zsolid1 ) / ( 1. + reac_pocl * dtsed2 ) 114 zrearat2(ji,jk) = ( reac_pocs * dtsed2 * zsolid2 ) / ( 1. + reac_pocs * dtsed2 ) 115 zrearat3(ji,jk) = ( reac_pocr * dtsed2 * zsolid3 ) / ( 1. + reac_pocr * dtsed2 ) 116 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zrearat1(ji,jk) / zvolc(ji,jk,jspoc) 117 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zrearat2(ji,jk) / zvolc(ji,jk,jspos) 118 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zrearat3(ji,jk) / zvolc(ji,jk,jspor) 119 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 120 186 121 ! For DIC 187 122 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 188 123 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 124 125 ! For alkalinity 126 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 127 189 128 ! For Phosphate (in mol/l) 190 129 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 130 191 131 ! For iron (in mol/l) 192 132 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 193 ! For alkalinity 194 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) 133 195 134 ! Ammonium 196 135 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 136 197 137 ! Ligands 198 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) 138 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 139 140 ENDDO 141 ENDDO 142 143 ! left hand side of coefficient matrix 144 DO jk = 2, jpksed 145 DO ji = 1, jpoce 146 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 147 zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * zreasat 148 zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 149 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 150 zlimo2(ji,jk) = zundsat(ji,jk) / ( zundsat(ji,jk) + xksedo2 ) 151 152 ! Oxygen 153 pwcp(ji,jk,jwoxy) = zundsat(ji,jk) + MIN( pwcp(ji,jk,jwoxy), rtrn ) 154 zreasat = zreasat * zlimo2(ji,jk) 155 156 ! Ligands 157 sedligand(ji,jk) = sedligand(ji,jk) - reac_ligc * sedligand(ji,jk) 199 158 ENDDO 200 159 ENDDO … … 205 164 !-------------------------------------------------------------------- 206 165 207 zrearat1(:,:) = 0. 208 zrearat2(:,:) = 0. 209 zrearat3(:,:) = 0. 210 211 zundsat(:,:) = pwcp(:,:,jwno3) 166 zundsat(:,:) = MAX(0., pwcp(:,:,jwno3) - rtrn) 212 167 213 168 DO jk = 2, jpksed 214 169 DO ji = 1, jpoce 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 232 DO jk = 2, jpksed 233 DO ji = 1, jpoce 234 jflag2: DO jn = 1, 10 235 zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) 236 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 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 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 265 zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 266 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 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 271 ENDDO 272 ENDDO 273 274 ! New dissolved concentrations 275 DO jk = 2, jpksed 276 DO ji = 1, jpoce 170 zlimtmp = ( 1.0 - zlimo2(ji,jk) ) 171 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 172 zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * zreasat * zlimtmp 173 zgamma = - xksedno3 * pwcp(ji,jk,jwno3) 174 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 175 zlimno3(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedno3 ) 176 277 177 ! For nitrates 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) 280 ! For DIC 281 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 282 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 283 ! For Phosphate (in mol/l) 284 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 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 178 pwcp(ji,jk,jwno3) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwno3), rtrn) 179 zreasat = zreasat * zlimno3(ji,jk) 180 289 181 ! For alkalinity 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 182 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * srDnit 293 183 ENDDO 294 184 ENDDO … … 299 189 !-------------------------------------------------------------------- 300 190 301 zrearat1(:,:) = 0. 302 zrearat2(:,:) = 0. 303 zrearat3(:,:) = 0. 304 305 zundsat(:,:) = solcp(:,:,jsfeo) 191 zundsat(:,:) = MAX(0., solcp(:,:,jsfeo) - rtrn) 306 192 307 193 DO jk = 2, jpksed 308 194 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 325 326 ! DO jn = 1, 5 327 DO jk = 2, jpksed 328 DO ji = 1, jpoce 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 357 358 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 370 371 ! New dissolved concentrations 372 DO jk = 2, jpksed 373 DO ji = 1, jpoce 374 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk) 195 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) 196 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 197 zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 198 zgamma = -xksedfeo * solcp(ji,jk,jsfeo) 199 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 200 zlimfeo(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedfeo ) 201 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) 202 375 203 ! 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 204 solcp(ji,jk,jsfeo) = zundsat(ji,jk) + MIN(solcp(ji,jk,jsfeo), rtrn) 205 380 206 ! 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 207 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * 4.0 * redfep 208 386 209 ! For alkalinity 387 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) +8.0 * zreasat388 ! Ammonium 389 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4210 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 8.0 * zreasat 211 212 ! Iron 390 213 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 391 214 ENDDO 392 215 ENDDO 393 216 394 !-------------------------------------------------------------------- 395 ! Begining POC denitrification and NO3- diffusion 396 ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 397 !-------------------------------------------------------------------- 398 399 zrearat1(:,:) = 0. 400 zrearat2(:,:) = 0. 401 zrearat3(:,:) = 0. 402 403 zundsat(:,:) = pwcp(:,:,jwso4) 217 ! !-------------------------------------------------------------------- 218 ! ! Begining POC denitrification and NO3- diffusion 219 ! ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 220 ! !-------------------------------------------------------------------- 221 ! 222 zundsat(:,:) = MAX(0., pwcp(:,:,jwso4) - rtrn) 404 223 405 224 DO jk = 2, jpksed 406 225 DO ji = 1, jpoce 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 426 DO jk = 2, jpksed 427 DO ji = 1, jpoce 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 456 END DO 457 END DO 458 459 ! New solid concentration values (jk=2 to jksed) for each couple 460 DO jk = 2, jpksed 461 DO ji = 1, jpoce 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 468 ENDDO 469 ENDDO 470 ! 471 ! New dissolved concentrations 472 DO jk = 2, jpksed 473 DO ji = 1, jpoce 226 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) - zlimfeo(ji,jk) 227 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 228 zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 229 zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 230 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 231 zlimso4(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedso4 ) 232 474 233 ! 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) 478 ! For DIC 479 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 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 234 pwcp(ji,jk,jwso4) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwso4), rtrn) 235 zreasat = zreasat * zlimso4(ji,jk) 236 pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) + 0.5 * zreasat 237 487 238 ! For alkalinity 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 239 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat 491 240 ENDDO 492 241 ENDDO … … 496 245 ! Patankar scheme 497 246 ! ------------------------------------------------------------------------------ 498 499 247 call sed_dsr_redoxb 500 248 … … 564 312 !! Arguments 565 313 ! --- local variables 566 INTEGER :: ji, jk, jn ! dummy looop indices314 INTEGER :: ji, jk, jn, jw ! dummy looop indices 567 315 568 316 REAL, DIMENSION(6) :: zsedtrn, zsedtra … … 582 330 zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 583 331 zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 584 zsedtra(:) = zsedtrn(:)332 zsedtra(:) = MAX(0., zsedtrn(:) - rtrn ) 585 333 586 334 ! First pass of the scheme. At the end, it is 1st order … … 592 340 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 593 341 IF ( zalpha == 0. ) THEN 594 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 595 ELSE 596 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 342 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 343 ELSE 344 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 345 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 597 346 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 598 347 ENDIF 599 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 348 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 600 349 zsedtra(5) = zbeta - zsedtra(4) 601 350 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) … … 606 355 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 607 356 IF ( zalpha == 0. ) THEN 608 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )357 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 609 358 ELSE 610 359 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 615 364 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 616 365 ! NH4 + O2 617 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 618 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 619 IF ( zalpha == 0. ) THEN 620 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 *zadsnh4 * dtsed2 / 2.0 )621 ELSE 622 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) &623 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )624 ENDIF 625 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 626 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 366 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 367 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) /zadsnh4 368 IF ( zalpha == 0. ) THEN 369 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 370 ELSE 371 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 372 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 373 ENDIF 374 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 375 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 627 376 ! FeS - O2 628 377 zalpha = zsedtra(1) - 2.0 * zsedtra(6) … … 630 379 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 631 380 IF ( zalpha == 0. ) THEN 632 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 )381 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 633 382 ELSE 634 383 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 658 407 zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 659 408 IF ( zalpha == 0. ) THEN 660 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 )409 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_feh2s * dtsed2 ) 661 410 ELSE 662 411 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & … … 668 417 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 669 418 pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 670 419 ! Fe - H2S 671 420 zalpha = zsedtra(2) - zsedtra(4) 672 421 zbeta = zsedtra(4) + zsedtra(6) … … 686 435 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 687 436 IF (zalpha == 0.) THEN 688 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. )437 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2. ) 689 438 ELSE 690 439 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 694 443 zsedtra(4) = zbeta - zsedtra(6) 695 444 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 445 696 446 ! NH4 + O2 697 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 698 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 699 IF ( zalpha == 0.) THEN700 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0)701 ELSE 702 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) &703 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )704 ENDIF 705 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 706 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 447 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 448 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) / zadsnh4 449 IF ( zalpha == 0. ) THEN 450 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 451 ELSE 452 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 453 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 454 ENDIF 455 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 456 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 707 457 ! H2S + O2 708 458 zalpha = zsedtra(1) - 2.0 * zsedtra(2) … … 710 460 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 711 461 IF ( zalpha == 0. ) THEN 712 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )462 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 713 463 ELSE 714 464 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 718 468 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 719 469 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 470 720 471 ! Fe + O2 721 472 zalpha = zsedtra(1) - 0.25 * zsedtra(4) … … 724 475 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 725 476 IF ( zalpha == 0. ) THEN 726 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 727 ELSE 728 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 477 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 478 ELSE 479 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 480 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 729 481 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 730 482 ENDIF … … 733 485 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 734 486 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 487 488 ! Update the concentrations after the secondary reactions 489 zsedtra(:) = zsedtra(:) + MIN( zsedtrn(:), rtrn ) 735 490 pwcp(ji,jk,jwoxy) = zsedtra(1) 736 491 pwcp(ji,jk,jwh2s) = zsedtra(2)
Note: See TracChangeset
for help on using the changeset viewer.
