[3443] | 1 | MODULE seddsr |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE seddsr *** |
---|
[10222] | 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 |
---|
[10222] | 9 | USE sed_oce |
---|
| 10 | USE sedini |
---|
| 11 | USE lib_mpp ! distribued memory computing library |
---|
| 12 | USE lib_fortran |
---|
[3443] | 13 | |
---|
[10222] | 14 | IMPLICIT NONE |
---|
| 15 | PRIVATE |
---|
| 16 | |
---|
[3443] | 17 | PUBLIC sed_dsr |
---|
| 18 | |
---|
| 19 | !! * Module variables |
---|
| 20 | |
---|
[10222] | 21 | REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density |
---|
[3443] | 22 | |
---|
[5215] | 23 | !! $Id$ |
---|
[3443] | 24 | CONTAINS |
---|
| 25 | |
---|
[15450] | 26 | SUBROUTINE sed_dsr( accmask ) |
---|
[3443] | 27 | !!---------------------------------------------------------------------- |
---|
| 28 | !! *** ROUTINE sed_dsr *** |
---|
| 29 | !! |
---|
| 30 | !! ** Purpose : computes pore water dissolution and reaction |
---|
| 31 | !! |
---|
[10222] | 32 | !! ** Methode : Computation of the redox reactions in sediment. |
---|
| 33 | !! The main redox reactions are solved in sed_dsr whereas |
---|
| 34 | !! the secondary reactions are solved in sed_dsr_redoxb. |
---|
| 35 | !! A strand spliting approach is being used here (see |
---|
| 36 | !! sed_dsr_redoxb for more information). |
---|
[3443] | 37 | !! |
---|
| 38 | !! History : |
---|
| 39 | !! ! 98-08 (E. Maier-Reimer, Christoph Heinze ) Original code |
---|
| 40 | !! ! 04-10 (N. Emprin, M. Gehlen ) f90 |
---|
| 41 | !! ! 06-04 (C. Ethe) Re-organization |
---|
[10222] | 42 | !! ! 19-08 (O. Aumont) Debugging and improvement of the model. |
---|
| 43 | !! The original method is replaced by a |
---|
| 44 | !! Strand splitting method which deals |
---|
| 45 | !! well with stiff reactions. |
---|
[3443] | 46 | !!---------------------------------------------------------------------- |
---|
| 47 | !! Arguments |
---|
| 48 | ! --- local variables |
---|
[15450] | 49 | INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask |
---|
[10222] | 50 | INTEGER :: ji, jk, js, jw, jn ! dummy looop indices |
---|
[3443] | 51 | |
---|
[15450] | 52 | REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite |
---|
| 53 | REAL(wp) :: zreasat |
---|
[3443] | 54 | !! |
---|
| 55 | !!---------------------------------------------------------------------- |
---|
| 56 | |
---|
[10222] | 57 | IF( ln_timing ) CALL timing_start('sed_dsr') |
---|
| 58 | ! |
---|
| 59 | ! Initializations |
---|
| 60 | !---------------------- |
---|
[3443] | 61 | |
---|
[15450] | 62 | zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. |
---|
| 63 | zlimso4(:,:) = 0. |
---|
[10222] | 64 | |
---|
[3443] | 65 | !---------------------------------------------------------- |
---|
[10222] | 66 | ! 5. Beginning of solid reaction |
---|
[3443] | 67 | !--------------------------------------------------------- |
---|
| 68 | |
---|
[15450] | 69 | ! Computes the SMS of the species which do not affect the remin |
---|
| 70 | ! processes (Alkalinity, PO4, NH4, and the release of Fe from |
---|
| 71 | ! biogenic Fe |
---|
| 72 | ! In the following, all loops start from jpk = 2 as reactions |
---|
| 73 | ! only occur in the sediment |
---|
| 74 | ! -------------------------------------------------------------- |
---|
[3443] | 75 | DO jk = 2, jpksed |
---|
| 76 | DO ji = 1, jpoce |
---|
[15450] | 77 | IF (accmask(ji) == 0 ) THEN |
---|
| 78 | zreasat = rearatpom(ji,jk) |
---|
| 79 | ! For alkalinity |
---|
| 80 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) |
---|
| 81 | ! For Phosphate (in mol/l) |
---|
| 82 | pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r |
---|
| 83 | ! For Ammonium |
---|
| 84 | pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) + zreasat * srno3 * radssol(jk,jwnh4) |
---|
| 85 | ! For iron (in mol/l) |
---|
| 86 | pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radssol(jk,jwfe2) |
---|
| 87 | ENDIF |
---|
| 88 | END DO |
---|
[3443] | 89 | ENDDO |
---|
| 90 | |
---|
[15450] | 91 | ! Computes SMS of oxygen |
---|
[10222] | 92 | DO jk = 2, jpksed |
---|
| 93 | DO ji = 1, jpoce |
---|
[15450] | 94 | IF (accmask(ji) == 0 ) THEN |
---|
| 95 | zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) |
---|
| 96 | zlimo2(ji,jk) = MAX( 0., zlimo2(ji,jk) ) |
---|
| 97 | zreasat = rearatpom(ji,jk) * zlimo2(ji,jk) |
---|
| 98 | ! Acid Silicic |
---|
| 99 | pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat |
---|
| 100 | ENDIF |
---|
[10222] | 101 | END DO |
---|
[3443] | 102 | ENDDO |
---|
| 103 | |
---|
| 104 | !-------------------------------------------------------------------- |
---|
[15450] | 105 | ! Denitrification |
---|
[3443] | 106 | !-------------------------------------------------------------------- |
---|
[10222] | 107 | DO jk = 2, jpksed |
---|
[3443] | 108 | DO ji = 1, jpoce |
---|
[15450] | 109 | IF (accmask(ji) == 0 ) THEN |
---|
| 110 | zlimno3(ji,jk) = ( 1.0 - zlimo2(ji,jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) |
---|
| 111 | zlimno3(ji,jk) = MAX(0., zlimno3(ji,jk) ) |
---|
| 112 | zreasat = rearatpom(ji,jk) * zlimno3(ji,jk) * srDnit |
---|
| 113 | ! For nitrates |
---|
| 114 | pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat |
---|
| 115 | ! For alkalinity |
---|
| 116 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat |
---|
| 117 | ENDIF |
---|
[3443] | 118 | END DO |
---|
| 119 | ENDDO |
---|
| 120 | |
---|
| 121 | |
---|
[10222] | 122 | !-------------------------------------------------------------------- |
---|
| 123 | ! Begining POC iron reduction |
---|
| 124 | !-------------------------------------------------------------------- |
---|
[3443] | 125 | |
---|
[10222] | 126 | DO jk = 2, jpksed |
---|
| 127 | DO ji = 1, jpoce |
---|
[15450] | 128 | IF (accmask(ji) == 0 ) THEN |
---|
| 129 | zlimfeo(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) |
---|
| 130 | zlimfeo(ji,jk) = MAX(0., zlimfeo(ji,jk) ) |
---|
| 131 | zreasat = 4.0 * rearatpom(ji,jk) * zlimfeo(ji,jk) |
---|
| 132 | ! For FEOH |
---|
| 133 | solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - zreasat / volc(ji,jk,jsfeo) |
---|
| 134 | ! For PO4 |
---|
| 135 | pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * redfep |
---|
| 136 | ! For alkalinity |
---|
| 137 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat |
---|
| 138 | ! Iron |
---|
| 139 | pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) |
---|
| 140 | ENDIF |
---|
[10222] | 141 | END DO |
---|
[3443] | 142 | ENDDO |
---|
| 143 | |
---|
[10222] | 144 | !-------------------------------------------------------------------- |
---|
| 145 | ! Begining POC denitrification and NO3- diffusion |
---|
| 146 | !-------------------------------------------------------------------- |
---|
[3443] | 147 | |
---|
[10222] | 148 | DO jk = 2, jpksed |
---|
[3443] | 149 | DO ji = 1, jpoce |
---|
[15450] | 150 | IF (accmask(ji) == 0 ) THEN |
---|
| 151 | zlimso4(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) - zlimfeo(ji,jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) |
---|
| 152 | zlimso4(ji,jk) = MAX(0., zlimso4(ji,jk) ) |
---|
| 153 | zreasat = 0.5 * rearatpom(ji,jk) * zlimso4(ji,jk) |
---|
| 154 | ! For sulfur |
---|
| 155 | pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - zreasat |
---|
| 156 | pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + zreasat |
---|
| 157 | ! For alkalinity |
---|
| 158 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat |
---|
| 159 | ENDIF |
---|
[3443] | 160 | END DO |
---|
| 161 | ENDDO |
---|
| 162 | |
---|
[15450] | 163 | ! Secondary redox reactions |
---|
| 164 | ! ------------------------- |
---|
[3443] | 165 | |
---|
[15450] | 166 | call sed_dsr_redoxb( accmask ) |
---|
[3443] | 167 | |
---|
[10222] | 168 | IF( ln_timing ) CALL timing_stop('sed_dsr') |
---|
| 169 | ! |
---|
| 170 | END SUBROUTINE sed_dsr |
---|
| 171 | |
---|
[15450] | 172 | SUBROUTINE sed_dsr_redoxb( accmask ) |
---|
[10222] | 173 | !!---------------------------------------------------------------------- |
---|
| 174 | !! *** ROUTINE sed_dsr_redox *** |
---|
| 175 | !! |
---|
| 176 | !! ** Purpose : computes secondary redox reactions |
---|
| 177 | !! |
---|
| 178 | !! History : |
---|
| 179 | !! ! 18-08 (O. Aumont) Original code |
---|
| 180 | !!---------------------------------------------------------------------- |
---|
| 181 | !! Arguments |
---|
| 182 | ! --- local variables |
---|
[15450] | 183 | INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask |
---|
| 184 | INTEGER :: ji, jni, jnj, jk ! dummy looop indices |
---|
[10222] | 185 | |
---|
[15450] | 186 | REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zreasat |
---|
[10222] | 187 | !! |
---|
| 188 | !!---------------------------------------------------------------------- |
---|
| 189 | |
---|
| 190 | IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') |
---|
| 191 | |
---|
[15450] | 192 | DO jk = 2, jpksed |
---|
| 193 | DO ji = 1, jpoce |
---|
| 194 | IF (accmask(ji) == 0 ) THEN |
---|
| 195 | zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 |
---|
| 196 | zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 |
---|
[10222] | 197 | |
---|
[15450] | 198 | ! First pass of the scheme. At the end, it is 1st order |
---|
| 199 | ! ----------------------------------------------------- |
---|
| 200 | ! Fe (both adsorbed and solute) + O2 |
---|
| 201 | zreasat = reac_fe2 * zsedfer / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) |
---|
| 202 | pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) |
---|
| 203 | pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat |
---|
| 204 | pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat |
---|
| 205 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat |
---|
| 206 | solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) |
---|
| 207 | |
---|
| 208 | ! H2S + O2 |
---|
| 209 | zreasat = reac_h2s * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) |
---|
| 210 | pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat |
---|
| 211 | pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat |
---|
| 212 | pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat |
---|
| 213 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat |
---|
| 214 | |
---|
| 215 | ! NH4 + O2 |
---|
| 216 | zreasat = reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) * pwcp(ji,jk,jwoxy) |
---|
| 217 | pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radssol(jk,jwnh4) |
---|
| 218 | pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat |
---|
| 219 | pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat |
---|
| 220 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat |
---|
| 221 | |
---|
| 222 | ! FeS - O2 |
---|
| 223 | zreasat = reac_feso * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) |
---|
| 224 | solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) |
---|
| 225 | pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat |
---|
| 226 | pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) |
---|
| 227 | pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat |
---|
| 228 | |
---|
| 229 | ! FEOH + H2S |
---|
| 230 | zreasat = reac_feh2s * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) |
---|
| 231 | solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) |
---|
| 232 | pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat |
---|
| 233 | pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radssol(jk,jwfe2) |
---|
| 234 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat |
---|
| 235 | pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat |
---|
| 236 | pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat |
---|
| 237 | |
---|
| 238 | ! Fe + H2S |
---|
| 239 | zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) |
---|
| 240 | zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 |
---|
| 241 | IF ( zexcess >= 0.0 ) THEN |
---|
| 242 | zreasat = reac_fesp * zexcess |
---|
| 243 | pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) |
---|
| 244 | solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) |
---|
| 245 | pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat |
---|
| 246 | ELSE |
---|
| 247 | zreasat = reac_fesd * zexcess * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) |
---|
| 248 | pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) |
---|
| 249 | solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) |
---|
| 250 | pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat |
---|
| 251 | ENDIF |
---|
| 252 | ! For alkalinity |
---|
| 253 | pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 |
---|
[10362] | 254 | ENDIF |
---|
[3443] | 255 | END DO |
---|
[15450] | 256 | END DO |
---|
[3443] | 257 | |
---|
[10222] | 258 | IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') |
---|
| 259 | |
---|
| 260 | END SUBROUTINE sed_dsr_redoxb |
---|
| 261 | |
---|
[3443] | 262 | END MODULE seddsr |
---|