- Timestamp:
- 2021-09-28T11:20:56+02: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
r15127 r15297 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce 9 10 USE sedini 10 11 USE lib_mpp ! distribued memory computing library … … 23 24 CONTAINS 24 25 25 SUBROUTINE sed_dsr( ji )26 SUBROUTINE sed_dsr( accmask ) 26 27 !!---------------------------------------------------------------------- 27 28 !! *** ROUTINE sed_dsr *** … … 45 46 !!---------------------------------------------------------------------- 46 47 !! Arguments 47 INTEGER, INTENT(in) :: ji ! number of iteration48 48 ! --- local variables 49 INTEGER :: jk, js, jw, jn ! dummy looop indices 50 51 REAL(wp), DIMENSION(jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 49 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 50 INTEGER :: ji, jk, js, jw, jn ! dummy looop indices 51 52 REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 52 53 REAL(wp) :: zreasat 53 54 !! … … 59 60 !---------------------- 60 61 61 zlimo2 (: ) = 0. ; zlimno3(:) = 0.62 zlimso4(: ) = 0.62 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. 63 zlimso4(:,:) = 0. 63 64 64 65 !---------------------------------------------------------- … … 73 74 ! -------------------------------------------------------------- 74 75 DO jk = 2, jpksed 75 zreasat = rearatpom(ji,jk) 76 ! For alkalinity 77 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 78 ! For Phosphate (in mol/l) 79 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 80 81 ! For iron (in mol/l) 82 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radsfe2(jk) 76 DO ji = 1, jpoce 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 83 89 ENDDO 84 90 85 91 ! Computes SMS of oxygen 86 92 DO jk = 2, jpksed 87 zlimo2(jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 88 zreasat = rearatpom(ji,jk) * zlimo2(jk) 89 ! Acid Silicic 90 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 93 DO ji = 1, jpoce 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 101 END DO 91 102 ENDDO 92 103 … … 95 106 !-------------------------------------------------------------------- 96 107 DO jk = 2, jpksed 97 zlimno3(jk) = ( 1.0 - zlimo2(jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 98 zreasat = rearatpom(ji,jk) * zlimno3(jk) 99 ! For nitrates 100 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat * srDnit 101 ! For alkalinity 102 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * srDnit 108 DO ji = 1, jpoce 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 118 END DO 103 119 ENDDO 104 120 … … 109 125 110 126 DO jk = 2, jpksed 111 zlimfeo(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 112 zreasat = rearatpom(ji,jk) * zlimfeo(jk) 113 ! For FEOH 114 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 4.0 * zreasat / volc(ji,jk,jsfeo) 115 ! For PO4 116 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * 4.0 * redfep 117 ! For alkalinity 118 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 8.0 * zreasat 119 ! Iron 120 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * 4.0 * radsfe2(jk) 127 DO ji = 1, jpoce 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 141 END DO 121 142 ENDDO 122 143 … … 126 147 127 148 DO jk = 2, jpksed 128 zlimso4(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) - zlimfeo(jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 129 zreasat = rearatpom(ji,jk) * zlimso4(jk) 130 ! For sulfur 131 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - 0.5 * zreasat 132 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + 0.5 * zreasat 133 ! For alkalinity 134 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 149 DO ji = 1, jpoce 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 160 END DO 135 161 ENDDO 136 162 … … 138 164 ! ------------------------- 139 165 140 call sed_dsr_redoxb( ji)166 call sed_dsr_redoxb( accmask ) 141 167 142 168 IF( ln_timing ) CALL timing_stop('sed_dsr') … … 144 170 END SUBROUTINE sed_dsr 145 171 146 SUBROUTINE sed_dsr_redoxb( ji)172 SUBROUTINE sed_dsr_redoxb( accmask ) 147 173 !!---------------------------------------------------------------------- 148 174 !! *** ROUTINE sed_dsr_redox *** … … 155 181 !! Arguments 156 182 ! --- local variables 157 INTEGER, INTENT(IN) :: ji158 INTEGER :: j ni, jnj, jk ! dummy looop indices183 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 184 INTEGER :: ji, jni, jnj, jk ! dummy looop indices 159 185 160 186 REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zreasat … … 165 191 166 192 DO jk = 2, jpksed 167 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 168 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 169 170 ! First pass of the scheme. At the end, it is 1st order 171 ! ----------------------------------------------------- 172 ! Fe (both adsorbed and solute) + O2 173 zreasat = reac_fe2 * dtsed * zsedfer / radsfe2(jk) * pwcp(ji,jk,jwoxy) 174 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 175 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 176 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 177 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 178 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 179 180 ! H2S + O2 181 zreasat = reac_h2s * dtsed * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 182 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 183 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 184 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 185 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 186 187 ! NH4 + O2 188 zreasat = reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) / radsnh4(jk) * pwcp(ji,jk,jwoxy) 189 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radsnh4(jk) 190 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 191 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 192 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 193 194 ! FeS - O2 195 zreasat = reac_feso * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 196 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 197 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 198 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radsfe2(jk) 199 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 200 201 ! FEOH + H2S 202 zreasat = reac_feh2s * dtsed * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 203 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 204 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 205 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radsfe2(jk) 206 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 207 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 208 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 209 210 ! Fe + H2S 211 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 212 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 213 IF ( zexcess >= 0.0 ) THEN 214 zreasat = reac_fesp * zexcess * dtsed 215 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 216 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 217 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 218 ELSE 219 zreasat = reac_fesd * zexcess * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 220 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 221 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 222 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 223 ENDIF 224 ! For alkalinity 225 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 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 197 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 254 ENDIF 255 END DO 226 256 END DO 227 257
Note: See TracChangeset
for help on using the changeset viewer.
