- Timestamp:
- 2018-11-16T15:59:30+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO4_beta_mirror/src/TOP/PISCES/SED/seddsr.F90
r9950 r10321 1 1 MODULE seddsr 2 #if defined key_sed3 2 !!====================================================================== 4 3 !! *** MODULE seddsr *** 5 !! Sediment : dissolution and reaction in pore water 4 !! Sediment : dissolution and reaction in pore water related 5 !! related to organic matter 6 6 !!===================================================================== 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sedmat ! linear system of equations 10 USE sedco3 ! carbonate ion and proton concentration 9 USE sed_oce 10 USE sedini 11 USE lib_mpp ! distribued memory computing library 12 USE lib_fortran 13 14 IMPLICIT NONE 15 PRIVATE 11 16 12 17 PUBLIC sed_dsr … … 14 19 !! * Module variables 15 20 16 REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_o2 17 REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_no3 18 REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_no3 19 REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_c13 20 REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: dens_mol_wgt ! molecular density 21 REAL(wp) :: zadsnh4 22 REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density 23 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables 24 21 25 22 26 !! $Id$ 23 27 CONTAINS 24 28 25 SUBROUTINE sed_dsr( kt )29 SUBROUTINE sed_dsr( kt, knt ) 26 30 !!---------------------------------------------------------------------- 27 31 !! *** ROUTINE sed_dsr *** … … 29 33 !! ** Purpose : computes pore water dissolution and reaction 30 34 !! 31 !! ** Methode : implicit simultaneous computation of undersaturation 32 !! resulting from diffusive pore water transport and chemical 33 !! pore water reactions. Solid material is consumed according 34 !! to redissolution and remineralisation 35 !! 36 !! ** Remarks : 37 !! - undersaturation : deviation from saturation concentration 38 !! - reaction rate : sink of undersaturation from dissolution 39 !! of solid material 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). 40 40 !! 41 41 !! History : … … 43 43 !! ! 04-10 (N. Emprin, M. Gehlen ) f90 44 44 !! ! 06-04 (C. Ethe) Re-organization 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. 45 49 !!---------------------------------------------------------------------- 46 50 !! Arguments 47 INTEGER, INTENT(in) :: kt ! number of iteration51 INTEGER, INTENT(in) :: kt, knt ! number of iteration 48 52 ! --- local variables 49 INTEGER :: ji, jk, js, jw ! dummy looop indices 50 INTEGER :: nv ! number of variables in linear tridiagonal eq 51 52 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zrearat ! reaction rate in pore water 53 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zundsat ! undersaturation ; indice jpwatp1 is for calcite 54 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zmo2_0, zmo2_1 ! temp. array for mass balance calculation 55 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zmno3_0, zmno3_1, zmno3_2 56 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zmc13_0, zmc13_1, zmc13_2, zmc13_3 57 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables 53 INTEGER :: ji, jk, js, jw, jn ! dummy looop indices 54 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 58 59 REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat 59 60 REAL(wp) :: zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 61 REAL(wp) :: zratio, zgamma, zbeta, zlimtmp, zundsat2 60 62 !! 61 63 !!---------------------------------------------------------------------- 62 64 63 IF( kt == nitsed000 ) THEN 64 WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' 65 WRITE(numsed,*) ' ' 66 ! 67 ALLOCATE( dens_mol_wgt((jpoce) ) 68 dens_mol_wgt(1:jpsol) = dens / mol_wgt(1:jpsol) 69 ! 70 ALLOCATE( cons_o2 (jpoce) ) ; ALLOCATE( cons_no3(jpoce) ) 71 ALLOCATE( sour_no3(jpoce) ) ; ALLOCATE( sour_c13(jpoce) ) 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 72 72 ENDIF 73 73 74 ! Initialization of data for mass balance calculation 75 !--------------------------------------------------- 76 77 tokbot(:,:) = 0. 78 cons_o2 (:) = 0. 79 cons_no3(:) = 0. 80 sour_no3(:) = 0. 81 sour_c13(:) = 0. 82 83 ! Initializations 84 !---------------------- 85 ALLOCATE( zmo2_0 (jpoce) ) ; ALLOCATE( zmo2_1 (jpoce) ) 86 ALLOCATE( zmno3_0(jpoce) ) ; ALLOCATE( zmno3_1(jpoce) ) ; ALLOCATE( zmno3_2(jpoce) ) 87 ALLOCATE( zmc13_0(jpoce) ) ; ALLOCATE( zmc13_1(jpoce) ) ; ALLOCATE( zmc13_2(jpoce) ) ; ALLOCATE( zmc13_3(jpoce) ) 88 89 zmo2_0 (:) = 0. ; zmo2_1 (:) = 0. 90 zmno3_0(:) = 0. ; zmno3_1(:) = 0. ; zmno3_2(:) = 0. 91 zmc13_0(:) = 0. ; zmc13_1(:) = 0. ; zmc13_2(:) = 0. ; zmc13_3(:) = 0. 74 ! Initializations 75 !---------------------- 92 76 93 ALLOCATE( zrearat(jpoce,jpksed,3) ) ; ALLOCATE( zundsat(jpoce,jpksed,3) ) 94 zrearat(:,:,:) = 0. ; zundsat(:,:,:) = 0. 95 96 97 ALLOCATE( zvolc(jpoce,jpksed,jpsol) ) 98 zvolc(:,:,:) = 0. 99 100 !-------------------------------------------------------------------- 101 ! Temporary accomodation to take account of particule rain deposition 102 !--------------------------------------------------------------------- 103 104 105 ! 1. Change of geometry 106 ! Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep 107 ! Warning : no change for dz(2) 108 !--------------------------------------------------------- 109 dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce) 110 111 112 ! New values for volw3d(:,2) and vols3d(:,2) 113 ! Warning : no change neither for volw(2) nor vols(2) 114 !------------------------------------------------------ 115 volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2) 116 vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2) 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 ) 86 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 117 96 118 97 ! Conversion of volume units … … 120 99 DO js = 1, jpsol 121 100 DO jk = 1, jpksed 122 DO ji = 1, jpoce 101 DO ji = 1, jpoce 123 102 zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / & 124 & ( volw3d(ji,jk) * 1.e-3 ) 103 & ( volw3d(ji,jk) * 1.e-3 ) 125 104 ENDDO 126 105 ENDDO 127 106 ENDDO 128 107 129 ! 2. Change of previous solid fractions (due to volum changes) for k=2130 !---------------------------------------------------------------------131 132 DO js = 1, jpsol133 DO ji = 1, jpoce134 solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2)135 ENDDO136 END DO137 138 ! 3. New solid fractions (including solid rain fractions) for k=2139 !------------------------------------------------------------------140 DO js = 1, jpsol141 DO ji = 1, jpoce142 solcp(ji,2,js) = solcp(ji,2,js) + &143 & ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) )144 ! rainrm are temporary cancel145 rainrm(ji,js) = 0.146 END DO147 ENDDO148 149 ! 4. Adjustment of bottom water concen.(pwcp(1)):150 ! We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume151 ! that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase.152 ! To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment153 ! of bottom water concentration.154 ! This adjustment is compensate at the end of routine155 !-------------------------------------------------------------156 DO jw = 1, jpwat157 DO ji = 1, jpoce158 pwcp(ji,1,jw) = pwcp(ji,1,jw) - &159 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)160 END DO161 ENDDO162 163 164 108 !---------------------------------------------------------- 165 ! 5. Beginning of Pore Water diffusion andsolid reaction109 ! 5. Beginning of solid reaction 166 110 !--------------------------------------------------------- 167 168 !-----------------------------------------------------------------------------169 ! For jk=2,jpksed, and for couple170 ! 1 : jwsil/jsopal ( SI/Opal )171 ! 2 : jsclay/jsclay ( clay/clay )172 ! 3 : jwoxy/jspoc ( O2/POC )173 ! reaction rate is a function of solid=concentration in solid reactif in [mol/l]174 ! and undersaturation in [mol/l].175 ! Solid weight fractions should be in ie [mol/l])176 ! second member and solution are in zundsat variable177 !-------------------------------------------------------------------------178 179 !number of variables180 nv = 3181 182 DO jk = 1, jpksed183 DO ji = 1, jpoce184 ! For Silicic Acid and clay185 zundsat(ji,jk,1) = sat_sil - pwcp(ji,jk,jwsil)186 zundsat(ji,jk,2) = sat_clay187 ! For O2188 zundsat(ji,jk,3) = pwcp(ji,jk,jwoxy) / so2ut189 ENDDO190 ENDDO191 192 111 193 112 ! Definition of reaction rates [rearat]=sans dim 194 113 ! For jk=1 no reaction (pure water without solid) for each solid compo 195 DO ji = 1, jpoce 196 zrearat(ji,1,:) = 0. 197 ENDDO 198 114 zrearat1(:,:) = 0. 115 zrearat2(:,:) = 0. 116 zrearat3(:,:) = 0. 117 118 zundsat(:,:) = pwcp(:,:,jwoxy) 119 120 DO jk = 2, jpksed 121 DO ji = 1, jpoce 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 ) 135 ENDDO 136 ENDDO 199 137 200 138 ! left hand side of coefficient matrix 201 DO jk = 2, jpksed 202 DO ji = 1, jpoce 203 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 204 zsolid2 = zvolc(ji,jk,jsclay) * solcp(ji,jk,jsclay) 205 zsolid3 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 206 207 zrearat(ji,jk,1) = ( reac_sil * dtsed * zsolid1 ) / & 208 & ( 1. + reac_sil * dtsed * zundsat(ji,jk,1 ) ) 209 zrearat(ji,jk,2) = ( reac_clay * dtsed * zsolid2 ) / & 210 & ( 1. + reac_clay * dtsed * zundsat(ji,jk,2 ) ) 211 zrearat(ji,jk,3) = ( reac_poc * dtsed * zsolid3 ) / & 212 & ( 1. + reac_poc * dtsed * zundsat(ji,jk,3 ) ) 213 ENDDO 214 ENDDO 215 216 217 CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 218 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 219 167 220 168 ! New solid concentration values (jk=2 to jksed) for each couple 221 DO js = 1, nv 222 DO jk = 2, jpksed 223 DO ji = 1, jpoce 224 zreasat = zrearat(ji,jk,js) * zundsat(ji,jk,js) / zvolc(ji,jk,js) 225 solcp(ji,jk,js) = solcp(ji,jk,js) - zreasat 226 ENDDO 227 ENDDO 228 ENDDO 229 ! mass of O2/NO3 before POC remin. for mass balance check 230 ! det. of o2 consomation/NO3 production Mc13 231 DO jk = 1, jpksed 232 DO ji = 1, jpoce 233 zvolw = volw3d(ji,jk) * 1.e-3 234 zmo2_0 (ji) = zmo2_0 (ji) + pwcp(ji,jk,jwoxy) * zvolw 235 zmno3_0(ji) = zmno3_0(ji) + pwcp(ji,jk,jwno3) * zvolw 236 zmc13_0(ji) = zmc13_0(ji) + pwcp(ji,jk,jwc13) * zvolw 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 237 177 ENDDO 238 178 ENDDO 239 179 240 180 ! New pore water concentrations 241 DO jk = 1, jpksed181 DO jk = 2, jpksed 242 182 DO ji = 1, jpoce 243 183 ! Acid Silicic 244 pwcp(ji,jk,jwsil) = sat_sil - zundsat(ji,jk,1) 245 ! For O2 (in mol/l) 246 pwcp(ji,jk,jwoxy) = zundsat(ji,jk,3) * so2ut 247 zreasat = zrearat(ji,jk,3) * zundsat(ji,jk,3) ! oxygen 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 248 186 ! For DIC 249 187 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 250 ! For nitrates 251 pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + zreasat * srno3 188 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 252 189 ! For Phosphate (in mol/l) 253 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 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 254 193 ! For alkalinity 255 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) - zreasat * ( srno3 + 2.* spo4r ) 256 ! For DIC13 257 pwcp(ji,jk,jwc13) = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb 258 ENDDO 259 ENDDO 260 261 262 ! Mass of O2 for mass balance check and det. of o2 consomation 263 DO jk = 1, jpksed 264 DO ji = 1, jpoce 265 zvolw = volw3d(ji,jk) * 1.e-3 266 zmo2_1 (ji) = zmo2_1 (ji) + pwcp(ji,jk,jwoxy) * zvolw 267 zmno3_1(ji) = zmno3_1(ji) + pwcp(ji,jk,jwno3) * zvolw 268 zmc13_1(ji) = zmc13_1(ji) + pwcp(ji,jk,jwc13) * zvolw 269 ENDDO 270 ENDDO 271 272 DO ji = 1, jpoce 273 cons_o2 (ji) = zmo2_0 (ji) - zmo2_1 (ji) 274 sour_no3(ji) = zmno3_1(ji) - zmno3_0(ji) 275 sour_c13(ji) = zmc13_1(ji) - zmc13_0(ji) 276 ENDDO 277 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) 199 ENDDO 200 ENDDO 278 201 279 202 !-------------------------------------------------------------------- … … 282 205 !-------------------------------------------------------------------- 283 206 284 nv = 1 285 DO jk = 1, jpksed 286 DO ji = 1, jpoce 287 zundsat(ji,jk,1) = pwcp(ji,jk,jwno3) / srDnit 288 ENDDO 289 ENDDO 290 DO jk = 2, jpksed 291 DO ji = 1, jpoce 292 IF( pwcp(ji,jk,jwoxy) < sthrO2 ) THEN 207 zrearat1(:,:) = 0. 208 zrearat2(:,:) = 0. 209 zrearat3(:,:) = 0. 210 211 zundsat(:,:) = pwcp(:,:,jwno3) 212 213 DO jk = 2, jpksed 214 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) ) 293 236 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 294 zrearat(ji,jk,1) = ( reac_no3 * dtsed * zsolid1 ) / & 295 & ( 1. + reac_no3 * dtsed * zundsat(ji,jk,1 ) ) 296 ELSE 297 zrearat(ji,jk,1) = 0. 298 ENDIF 299 END DO 300 END DO 301 302 303 ! solves tridiagonal system 304 CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 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 305 260 306 261 … … 308 263 DO jk = 2, jpksed 309 264 DO ji = 1, jpoce 310 zreasat = zrearat (ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jspoc)265 zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 311 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 312 271 ENDDO 313 272 ENDDO 314 273 315 274 ! New dissolved concentrations 316 DO jk = 1, jpksed 317 DO ji = 1, jpoce 318 zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) 275 DO jk = 2, jpksed 276 DO ji = 1, jpoce 319 277 ! For nitrates 320 pwcp(ji,jk,jwno3) = zundsat(ji,jk,1) * srDnit 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) 321 280 ! For DIC 322 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 323 283 ! For Phosphate (in mol/l) 324 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 325 289 ! For alkalinity 326 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srDnit - 2.* spo4r ) 327 ! For DIC13 328 pwcp(ji,jk,jwc13) = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb 329 ENDDO 330 ENDDO 331 332 333 ! Mass of O2 for mass balance check and det. of o2 consomation 334 DO jk = 1, jpksed 335 DO ji = 1, jpoce 336 zvolw = volw3d(ji,jk) * 1.e-3 337 zmno3_2(ji) = zmno3_2(ji) + pwcp(ji,jk ,jwno3) * zvolw 338 zmc13_2(ji) = zmc13_2(ji) + pwcp(ji,jk ,jwc13) * zvolw 339 ENDDO 340 ENDDO 341 342 DO ji = 1, jpoce 343 cons_no3(ji) = zmno3_1(ji) - zmno3_2(ji) 344 sour_c13(ji) = sour_c13(ji) + zmc13_2(ji) - zmc13_1(ji) 345 ENDDO 346 347 348 !--------------------------- 349 ! Solves PO4 diffusion 350 !---------------------------- 351 352 nv = 1 353 DO jk = 1, jpksed 354 DO ji = 1, jpoce 355 zundsat(ji,jk,1) = pwcp(ji,jk,jwpo4) 356 zrearat(ji,jk,1) = 0. 357 ENDDO 358 ENDDO 359 360 361 ! solves tridiagonal system 362 CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 363 364 365 ! New undsaturation values and dissolved concentrations 366 DO jk = 1, jpksed 367 DO ji = 1, jpoce 368 pwcp(ji,jk,jwpo4) = zundsat(ji,jk,1) 369 ENDDO 370 ENDDO 371 372 373 !--------------------------------------------------------------- 374 ! Performs CaCO3 particle deposition and redissolution (indice 9) 375 !-------------------------------------------------------------- 376 377 ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 378 379 CALL sed_co3( kt ) 380 381 382 nv = 1 383 ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 384 DO jk = 1, jpksed 385 DO ji = 1, jpoce 386 zundsat(ji,jk,1) = aksps(ji) * densSW(ji) * densSW(ji) / calcon2(ji) & 387 & - co3por(ji,jk) 388 ! positive values of undersaturation 389 zundsat(ji,jk,1) = MAX( 0., zundsat(ji,jk,1) ) 390 ENDDO 391 ENDDO 392 393 DO jk = 2, jpksed 394 DO ji = 1, jpoce 395 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 396 zrearat(ji,jk,1) = ( reac_cal * dtsed * zsolid1 ) / & 397 & ( 1. + reac_cal * dtsed * zundsat(ji,jk,1) ) 398 END DO 399 END DO 400 401 402 ! solves tridiagonal system 403 CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 404 405 406 ! New solid concentration values (jk=2 to jksed) for cacO3 407 DO jk = 2, jpksed 408 DO ji = 1, jpoce 409 zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jscal) 410 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 411 ENDDO 412 ENDDO 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 293 ENDDO 294 ENDDO 295 296 !-------------------------------------------------------------------- 297 ! Begining POC iron reduction 298 ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) 299 !-------------------------------------------------------------------- 300 301 zrearat1(:,:) = 0. 302 zrearat2(:,:) = 0. 303 zrearat3(:,:) = 0. 304 305 zundsat(:,:) = solcp(:,:,jsfeo) 306 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 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 413 370 414 371 ! New dissolved concentrations 415 DO jk = 1, jpksed 416 DO ji = 1, jpoce 417 zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) 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) 375 ! For FEOH 376 solcp(ji,jk,jsfeo) = zundsat(ji,jk) 418 377 ! For DIC 419 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 420 386 ! For alkalinity 421 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.* zreasat 422 ! For DIC13 423 pwcp(ji,jk,jwc13) = pwcp(ji,jk,jwc13) + zreasat * rc13Ca * pdb 424 ENDDO 425 ENDDO 426 427 DO jk = 1, jpksed 428 DO ji = 1, jpoce 429 zmc13_3(ji) = zmc13_3(ji) + pwcp(ji,jk,jwc13) * volw3d(ji,jk) * 1.e-3 430 ENDDO 431 ENDDO 432 433 DO ji = 1, jpoce 434 sour_c13(ji) = sour_c13(ji) + zmc13_3(ji) - zmc13_2(ji) 435 ENDDO 436 437 !------------------------------------------------- 438 ! Beginning DIC, Alkalinity and DIC13 diffusion 439 !------------------------------------------------- 440 441 nv = 3 442 DO jk = 1, jpksed 443 DO ji = 1, jpoce 444 zundsat(ji,jk,1) = pwcp(ji,jk,jwdic) 445 zundsat(ji,jk,2) = pwcp(ji,jk,jwalk) 446 zundsat(ji,jk,3) = pwcp(ji,jk,jwc13) 447 448 zrearat(ji,jk,1) = 0. 449 zrearat(ji,jk,2) = 0. 450 zrearat(ji,jk,3) = 0. 451 452 ENDDO 453 ENDDO 454 455 456 ! solves tridiagonal system 457 CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 458 459 460 ! New dissolved concentrations 461 DO jk = 1, jpksed 462 DO ji = 1, jpoce 463 pwcp(ji,jk,jwdic) = zundsat(ji,jk,1) 464 pwcp(ji,jk,jwalk) = zundsat(ji,jk,2) 465 pwcp(ji,jk,jwc13) = zundsat(ji,jk,3) 466 ENDDO 467 ENDDO 468 469 !---------------------------------- 470 ! Back to initial geometry 471 !----------------------------- 472 473 !--------------------------------------------------------------------- 474 ! 1/ Compensation for ajustement of the bottom water concentrations 475 ! (see note n° 1 about *por(2)) 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 391 ENDDO 392 ENDDO 393 476 394 !-------------------------------------------------------------------- 477 DO jw = 1, jpwat 478 DO ji = 1, jpoce 479 pwcp(ji,1,jw) = pwcp(ji,1,jw) + & 480 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji) 481 END DO 482 ENDDO 483 484 !----------------------------------------------------------------------- 485 ! 2/ Det of new rainrg taking account of the new weight fraction obtained 486 ! in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!) 487 ! This new rain (rgntg rm) will be used in advection/burial routine 488 !------------------------------------------------------------------------ 489 DO js = 1, jpsol 490 DO ji = 1, jpoce 491 rainrg(ji,js) = raintg(ji) * solcp(ji,2,js) 492 rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js) 493 END DO 494 ENDDO 495 496 ! New raintg 497 raintg(:) = 0. 498 DO js = 1, jpsol 499 DO ji = 1, jpoce 500 raintg(ji) = raintg(ji) + rainrg(ji,js) 501 END DO 502 ENDDO 503 504 !-------------------------------- 505 ! 3/ back to initial geometry 506 !-------------------------------- 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) 404 405 DO jk = 2, jpksed 406 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 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) 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 487 ! 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 491 ENDDO 492 ENDDO 493 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 ! ------------------------------------------------------------------------------ 498 499 call sed_dsr_redoxb 500 501 ! -------------------------------------------------------------- 502 ! 4/ Computation of the bioturbation coefficient 503 ! This parameterization is taken from Archer et al. (2002) 504 ! -------------------------------------------------------------- 505 507 506 DO ji = 1, jpoce 508 dz3d (ji,2) = dz(2) 509 volw3d(ji,2) = dz3d(ji,2) * por(2) 510 vols3d(ji,2) = dz3d(ji,2) * por1(2) 511 ENDDO 512 513 !---------------------------------------------------------------------- 514 ! 4/ Saving new amount of material in dzkbot for mass balance check 515 ! tokbot in [mol] (implicit *1cm*1cm for spacial dim) 516 !---------------------------------------------------------------------- 517 DO jw = 1, jpwat 518 DO ji = 1, jpoce 519 tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 520 END DO 521 ENDDO 522 523 DEALLOCATE( zmo2_0 ) ; DEALLOCATE( zmno3_1 ) ; DEALLOCATE( zmno3_2 ) 524 DEALLOCATE( zmc13_0 ) ; DEALLOCATE( zmc13_1 ) ; DEALLOCATE( zmc13_2 ) ; DEALLOCATE( zmc13_3 ) 525 526 DEALLOCATE( zrearat ) ; DEALLOCATE( zundsat ) ; DEALLOCATE( zvolc ) 527 507 db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) 508 END DO 509 510 ! ------------------------------------------------------ 511 ! Vertical variations of the bioturbation coefficient 512 ! ------------------------------------------------------ 513 IF (ln_btbz) THEN 514 DO ji = 1, jpoce 515 db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) 516 END DO 517 ELSE 518 DO jk = 1, jpksed 519 IF (profsedw(jk) > dbtbzsc) THEN 520 db(:,jk) = 0.0 521 ENDIF 522 END DO 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 532 END DO 533 ELSE 534 irrig(:,:) = 0.0 535 ENDIF 536 537 DEALLOCATE( zvolc ) 538 539 IF( ln_timing ) CALL timing_stop('sed_dsr') 540 ! 528 541 END SUBROUTINE sed_dsr 529 #else 530 !!====================================================================== 531 !! MODULE seddsr : Dummy module 532 !!====================================================================== 533 !! $Id$ 534 CONTAINS 535 SUBROUTINE sed_dsr ( kt ) 536 INTEGER, INTENT(in) :: kt 537 WRITE(*,*) 'sed_dsr: You should not have seen this print! error?', kt 538 END SUBROUTINE sed_dsr 539 #endif 540 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 575 DO ji = 1, jpoce 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) 698 END DO 699 END DO 700 701 IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') 702 703 END SUBROUTINE sed_dsr_redoxb 704 541 705 END MODULE seddsr
Note: See TracChangeset
for help on using the changeset viewer.