Changeset 14323 for NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90
- Timestamp:
- 2021-01-21T00:33:21+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90
r12360 r14323 23 23 CONTAINS 24 24 25 SUBROUTINE sed_inorg( kt )25 SUBROUTINE sed_inorg( kt, knt ) 26 26 !!---------------------------------------------------------------------- 27 27 !! *** ROUTINE sed_inorg *** … … 46 46 !!---------------------------------------------------------------------- 47 47 !! Arguments 48 INTEGER, INTENT(in) :: kt ! number of iteration48 INTEGER, INTENT(in) :: kt, knt ! number of iteration 49 49 ! --- local variables 50 50 INTEGER :: ji, jk, js, jw ! dummy looop indices … … 54 54 REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc ! temp. variables 55 55 REAL(wp), DIMENSION(jpoce) :: zsieq 56 REAL(wp) :: zsolid1, zvolw, zreasat 56 REAL(wp) :: zsolid1, zvolw, zreasat, zbeta, zgamma 57 57 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 58 58 !! … … 61 61 IF( ln_timing ) CALL timing_start('sed_inorg') 62 62 ! 63 IF( kt == nitsed000 ) THEN63 IF( kt == nitsed000 .AND. knt == 1 ) THEN 64 64 IF (lwp) THEN 65 65 WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' … … 73 73 74 74 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. 75 zrearat2(:,:) = 0. ; zrearat2(:,:) = 0.75 zrearat2(:,:) = 0. 76 76 zco3eq(:) = rtrn 77 77 zvolc(:,:,:) = 0. … … 86 86 zsolcpsi = 0.0 87 87 DO jk = 1, jpksed 88 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk)89 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)88 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 89 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 90 90 END DO 91 91 zsolcpsi = MAX( zsolcpsi, rtrn ) … … 135 135 DO jk = 2, jpksed 136 136 DO ji = 1, jpoce 137 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 138 zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) ) 139 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 140 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji) 141 zrearat1(ji,jk) = ( reac_sil * znusil * dtsed * zsolid1 ) / & 142 & ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) ) 143 ENDDO 144 ENDDO 145 146 CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 147 148 ! New solid concentration values (jk=2 to jksed) for each couple 149 DO jk = 2, jpksed 150 DO ji = 1, jpoce 151 zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) ) 152 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat 153 ENDDO 154 ENDDO 155 156 ! New pore water concentrations 157 DO jk = 1, jpksed 158 DO ji = 1, jpoce 159 ! Acid Silicic 160 pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 137 IF ( zundsat(ji,jk) > 0. ) THEN 138 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 139 zsatur = zundsat(ji,jk) / zsieq(ji) 140 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 141 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji) 142 zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1 143 zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk) 144 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 145 solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal) 146 pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 147 ENDIF 161 148 ENDDO 162 149 ENDDO … … 175 162 zco3eq(ji) = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 176 163 zco3eq(ji) = MAX( rtrn, zco3eq(ji) ) 177 zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk))164 zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji) 178 165 ENDDO 179 166 ENDDO … … 181 168 DO jk = 2, jpksed 182 169 DO ji = 1, jpoce 183 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 184 zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / & 185 & ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) ) 170 IF ( zundsat(ji,jk) > 0. ) THEN 171 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 172 zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1 173 zgamma = -reac_cal * dtsed2 * zundsat(ji,jk) 174 zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 175 saturco3(ji,jk) = zundsat(ji,jk) 176 zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) ) 177 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal) 178 ! For DIC 179 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 180 ! For alkalinity 181 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 182 ENDIF 186 183 END DO 187 184 END DO 188 189 ! solves tridiagonal system190 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )191 192 ! New solid concentration values (jk=2 to jksed) for cacO3193 DO jk = 2, jpksed194 DO ji = 1, jpoce195 zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)196 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat197 ENDDO198 ENDDO199 200 ! New dissolved concentrations201 DO jk = 1, jpksed202 DO ji = 1, jpoce203 zreasat = zrearat1(ji,jk) * zundsat(ji,jk)204 ! For DIC205 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat206 ! For alkalinity207 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat208 ENDDO209 ENDDO210 211 !-------------------------------------------------212 ! Beginning DIC, Alkalinity213 !-------------------------------------------------214 215 DO jk = 1, jpksed216 DO ji = 1, jpoce217 zundsat(ji,jk) = pwcp(ji,jk,jwdic)218 zrearat1(ji,jk) = 0.219 ENDDO220 ENDDO221 222 ! solves tridiagonal system223 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )224 225 ! New dissolved concentrations226 DO jk = 1, jpksed227 DO ji = 1, jpoce228 pwcp(ji,jk,jwdic) = zundsat(ji,jk)229 ENDDO230 ENDDO231 232 !-------------------------------------------------233 ! Beginning DIC, Alkalinity234 !-------------------------------------------------235 236 DO jk = 1, jpksed237 DO ji = 1, jpoce238 zundsat(ji,jk) = pwcp(ji,jk,jwalk)239 zrearat1(ji,jk) = 0.240 ENDDO241 ENDDO242 !243 ! ! solves tridiagonal system244 CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )245 !246 ! ! New dissolved concentrations247 DO jk = 1, jpksed248 DO ji = 1, jpoce249 pwcp(ji,jk,jwalk) = zundsat(ji,jk)250 ENDDO251 ENDDO252 253 !----------------------------------254 ! Back to initial geometry255 !-----------------------------256 257 !---------------------------------------------------------------------258 ! 1/ Compensation for ajustement of the bottom water concentrations259 ! (see note n° 1 about *por(2))260 !--------------------------------------------------------------------261 DO jw = 1, jpwat262 DO ji = 1, jpoce263 pwcp(ji,1,jw) = pwcp(ji,1,jw) + &264 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)265 END DO266 ENDDO267 268 !-----------------------------------------------------------------------269 ! 2/ Det of new rainrg taking account of the new weight fraction obtained270 ! in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)271 ! This new rain (rgntg rm) will be used in advection/burial routine272 !------------------------------------------------------------------------273 DO js = 1, jpsol274 DO ji = 1, jpoce275 rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)276 rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)277 END DO278 ENDDO279 280 ! New raintg281 raintg(:) = 0.282 DO js = 1, jpsol283 DO ji = 1, jpoce284 raintg(ji) = raintg(ji) + rainrg(ji,js)285 END DO286 ENDDO287 288 !--------------------------------289 ! 3/ back to initial geometry290 !--------------------------------291 DO ji = 1, jpoce292 dz3d (ji,2) = dz(2)293 volw3d(ji,2) = dz3d(ji,2) * por(2)294 vols3d(ji,2) = dz3d(ji,2) * por1(2)295 ENDDO296 185 297 186 IF( ln_timing ) CALL timing_stop('sed_inorg')
Note: See TracChangeset
for help on using the changeset viewer.
