- Timestamp:
- 2020-08-03T11:46:38+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zche.F90
r13295 r13373 188 188 ! CHEMICAL CONSTANTS - SURFACE LAYER 189 189 ! ---------------------------------- 190 !CDIR NOVERRCHK 191 DO jj = 1, jpj 192 !CDIR NOVERRCHK 193 DO ji = 1, jpi 194 ! ! SET ABSOLUTE TEMPERATURE 195 ztkel = tempis(ji,jj,1) + 273.15 196 zt = ztkel * 0.01 197 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 198 ! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 199 ! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 200 zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel & 201 & + 0.0047036e-4*ztkel**2) 202 chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 203 chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 204 chemc(ji,jj,3) = 57.7 - 0.118*ztkel 205 ! 206 END DO 207 END DO 190 DO_2D( 1, 1, 1, 1 ) 191 ! ! SET ABSOLUTE TEMPERATURE 192 ztkel = tempis(ji,jj,1) + 273.15 193 zt = ztkel * 0.01 194 zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35. 195 ! ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980) 196 ! ! AND FOR THE ATMOSPHERE FOR NON IDEAL GAS 197 zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel & 198 & + 0.0047036e-4*ztkel**2) 199 chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm) 200 chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3 201 chemc(ji,jj,3) = 57.7 - 0.118*ztkel 202 END_2D 208 203 209 204 ! OXYGEN SOLUBILITY - DEEP OCEAN 210 205 ! ------------------------------- 211 !CDIR NOVERRCHK 212 DO jk = 1, jpk 213 !CDIR NOVERRCHK 214 DO jj = 1, jpj 215 !CDIR NOVERRCHK 216 DO ji = 1, jpi 217 ztkel = tempis(ji,jj,jk) + 273.15 218 zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35. 219 zsal2 = zsal * zsal 220 ztgg = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature 221 ztgg2 = ztgg * ztgg 222 ztgg3 = ztgg2 * ztgg 223 ztgg4 = ztgg3 * ztgg 224 ztgg5 = ztgg4 * ztgg 225 226 zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 & 227 & + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 & 228 & - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) & 229 & - 3.11680e-7 * zsal2 230 chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm) 231 END DO 232 END DO 233 END DO 206 DO_3D( 1, 1, 1, 1, 1, jpk ) 207 ztkel = tempis(ji,jj,jk) + 273.15 208 zsal = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35. 209 zsal2 = zsal * zsal 210 ztgg = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel ) ! Set the GORDON & GARCIA scaled temperature 211 ztgg2 = ztgg * ztgg 212 ztgg3 = ztgg2 * ztgg 213 ztgg4 = ztgg3 * ztgg 214 ztgg5 = ztgg4 * ztgg 215 216 zoxy = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3 & 217 & + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3 & 218 & - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 ) & 219 & - 3.11680e-7 * zsal2 220 chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox ! mol/(L atm) 221 END_3D 234 222 235 223 ! CHEMICAL CONSTANTS - DEEP OCEAN 236 224 ! ------------------------------- 237 !CDIR NOVERRCHK 238 DO jk = 1, jpk 239 !CDIR NOVERRCHK 240 DO jj = 1, jpj 241 !CDIR NOVERRCHK 242 DO ji = 1, jpi 243 244 ! SET PRESSION ACCORDING TO SAUNDER (1980) 245 zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 246 zc1 = 5.92E-3 + zplat**2 * 5.25E-3 247 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6 248 zpres = zpres / 10.0 249 250 ! SET ABSOLUTE TEMPERATURE 251 ztkel = tempis(ji,jj,jk) + 273.15 252 zsal = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 253 zsqrt = SQRT( zsal ) 254 zsal15 = zsqrt * zsal 255 zlogt = LOG( ztkel ) 256 ztr = 1. / ztkel 257 zis = 19.924 * zsal / ( 1000.- 1.005 * zsal ) 258 zis2 = zis * zis 259 zisqrt = SQRT( zis ) 260 ztc = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 261 262 ! CHLORINITY (WOOSTER ET AL., 1969) 263 zcl = zsal / 1.80655 264 265 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 266 zst = 0.14 * zcl /96.062 267 268 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 269 zft = 0.000067 * zcl /18.9984 270 271 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 272 zcks = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt & 273 & + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 274 & + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis & 275 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 276 & + LOG(1.0 - 0.001005 * zsal)) 277 278 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 279 zckf = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt & 280 & + LOG(1.0d0 - 0.001005d0*zsal) & 281 & + LOG(1.0d0 + zst/zcks)) 282 283 ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 284 zckb= (-8966.90 - 2890.53*zsqrt - 77.942*zsal & 285 & + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr & 286 & + (148.0248 + 137.1942*zsqrt + 1.62142*zsal) & 287 & + (-24.4344 - 25.085*zsqrt - 0.2474*zsal) & 288 & * zlogt + 0.053105*zsqrt*ztkel 289 290 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 291 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 292 zck1 = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt & 293 - 0.011555*zsal + 0.0001152*zsal*zsal) 294 zck2 = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt & 295 - 0.01781*zsal + 0.0001122*zsal*zsal) 296 297 ! PKW (H2O) (MILLERO, 1995) from composite data 298 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 299 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 300 301 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 302 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 303 & + (-106.736*ztr + 0.69171) * zsqrt & 304 & + (-0.65643*ztr - 0.01844) * zsal 305 306 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 307 & + (-160.340*ztr + 1.3566)*zsqrt & 308 & + (0.37335*ztr - 0.05778)*zsal 309 310 zck3p = -3070.75*ztr - 18.126 & 311 & + (17.27039*ztr + 2.81197) * zsqrt & 312 & + (-44.99486*ztr - 0.09984) * zsal 313 314 ! CONSTANT FOR SILICATE, MILLERO (1995) 315 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 316 & + (-458.79*ztr + 3.5913) * zisqrt & 317 & + (188.74*ztr - 1.5998) * zis & 318 & + (-12.1652*ztr + 0.07871) * zis2 & 319 & + LOG(1.0 - 0.001005*zsal) 320 321 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 322 ! (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 323 zaksp0 = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel ) & 324 & + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt & 325 & - 0.07711*zsal + 0.0041249*zsal15 326 327 ! CONVERT FROM DIFFERENT PH SCALES 328 total2free = 1.0/(1.0 + zst/zcks) 329 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 330 total2SWS = total2free * free2SWS 331 SWS2total = 1.0 / total2SWS 332 333 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 334 zak1 = 10**(zck1) * total2SWS 335 zak2 = 10**(zck2) * total2SWS 336 zakb = EXP( zckb ) * total2SWS 337 zakw = EXP( zckw ) 338 zaksp1 = 10**(zaksp0) 339 zak1p = exp( zck1p ) 340 zak2p = exp( zck2p ) 341 zak3p = exp( zck3p ) 342 zaksi = exp( zcksi ) 343 zckf = zckf * total2SWS 344 345 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 346 ! (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE 347 ! IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 348 ! TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN 349 ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 350 ! MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 351 ! WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 352 ! & GIESKES (1970), P. 1285-1286 (THE SMALL 353 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 354 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 355 zcpexp = zpres / (rgas*ztkel) 356 zcpexp2 = zpres * zcpexp 357 358 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 359 ! CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 360 ! (CF. BROECKER ET AL., 1982) 361 362 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 363 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 364 ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 365 366 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 367 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 368 ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 369 370 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 371 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 372 akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 373 374 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 375 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 376 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 377 378 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 379 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 380 aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 381 382 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 383 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 384 akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 385 386 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 387 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 388 ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 389 390 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 391 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 392 ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 393 394 zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 395 zbuf2 = 0.5 * ( devk49 + devk59 * ztc ) 396 ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 397 398 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 399 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 400 aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 401 402 ! CONVERT FROM DIFFERENT PH SCALES 403 total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 404 free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 405 total2SWS = total2free * free2SWS 406 SWS2total = 1.0 / total2SWS 407 408 ! Convert to total scale 409 ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total 410 ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total 411 akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total 412 akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total 413 ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 414 ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 415 ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 416 aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 417 akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free 418 419 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 420 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 421 ! (P. 1285) AND BERNER (1976) 422 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 423 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 424 aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 425 426 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 427 borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 428 sulfat(ji,jj,jk) = zst 429 fluorid(ji,jj,jk) = zft 430 431 ! Iron and SIO3 saturation concentration from ... 432 sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 433 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 434 435 ! Liu and Millero (1999) only valid 5 - 50 degC 436 ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16 437 fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1) 438 fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 ) 439 fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 ) 440 fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 ) 441 fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 ) 442 END DO 443 END DO 444 END DO 225 DO_3D( 1, 1, 1, 1, 1, jpk ) 226 ! SET PRESSION ACCORDING TO SAUNDER (1980) 227 zplat = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) ) 228 zc1 = 5.92E-3 + zplat**2 * 5.25E-3 229 zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6 230 zpres = zpres / 10.0 231 232 ! SET ABSOLUTE TEMPERATURE 233 ztkel = tempis(ji,jj,jk) + 273.15 234 zsal = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35. 235 zsqrt = SQRT( zsal ) 236 zsal15 = zsqrt * zsal 237 zlogt = LOG( ztkel ) 238 ztr = 1. / ztkel 239 zis = 19.924 * zsal / ( 1000.- 1.005 * zsal ) 240 zis2 = zis * zis 241 zisqrt = SQRT( zis ) 242 ztc = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20. 243 244 ! CHLORINITY (WOOSTER ET AL., 1969) 245 zcl = zsal / 1.80655 246 247 ! TOTAL SULFATE CONCENTR. [MOLES/kg soln] 248 zst = 0.14 * zcl /96.062 249 250 ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln] 251 zft = 0.000067 * zcl /18.9984 252 253 ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990) 254 zcks = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt & 255 & + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt & 256 & + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis & 257 & - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2 & 258 & + LOG(1.0 - 0.001005 * zsal)) 259 260 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 261 zckf = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt & 262 & + LOG(1.0d0 - 0.001005d0*zsal) & 263 & + LOG(1.0d0 + zst/zcks)) 264 265 ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE 266 zckb= (-8966.90 - 2890.53*zsqrt - 77.942*zsal & 267 & + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr & 268 & + (148.0248 + 137.1942*zsqrt + 1.62142*zsal) & 269 & + (-24.4344 - 25.085*zsqrt - 0.2474*zsal) & 270 & * zlogt + 0.053105*zsqrt*ztkel 271 272 ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 273 ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale 274 zck1 = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt & 275 - 0.011555*zsal + 0.0001152*zsal*zsal) 276 zck2 = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt & 277 - 0.01781*zsal + 0.0001122*zsal*zsal) 278 279 ! PKW (H2O) (MILLERO, 1995) from composite data 280 zckw = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr & 281 - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal 282 283 ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995) 284 zck1p = -4576.752*ztr + 115.540 - 18.453*zlogt & 285 & + (-106.736*ztr + 0.69171) * zsqrt & 286 & + (-0.65643*ztr - 0.01844) * zsal 287 288 zck2p = -8814.715*ztr + 172.1033 - 27.927*zlogt & 289 & + (-160.340*ztr + 1.3566)*zsqrt & 290 & + (0.37335*ztr - 0.05778)*zsal 291 292 zck3p = -3070.75*ztr - 18.126 & 293 & + (17.27039*ztr + 2.81197) * zsqrt & 294 & + (-44.99486*ztr - 0.09984) * zsal 295 296 ! CONSTANT FOR SILICATE, MILLERO (1995) 297 zcksi = -8904.2*ztr + 117.400 - 19.334*zlogt & 298 & + (-458.79*ztr + 3.5913) * zisqrt & 299 & + (188.74*ztr - 1.5998) * zis & 300 & + (-12.1652*ztr + 0.07871) * zis2 & 301 & + LOG(1.0 - 0.001005*zsal) 302 303 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER 304 ! (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983) 305 zaksp0 = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel ) & 306 & + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt & 307 & - 0.07711*zsal + 0.0041249*zsal15 308 309 ! CONVERT FROM DIFFERENT PH SCALES 310 total2free = 1.0/(1.0 + zst/zcks) 311 free2SWS = 1. + zst/zcks + zft/(zckf*total2free) 312 total2SWS = total2free * free2SWS 313 SWS2total = 1.0 / total2SWS 314 315 ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?) 316 zak1 = 10**(zck1) * total2SWS 317 zak2 = 10**(zck2) * total2SWS 318 zakb = EXP( zckb ) * total2SWS 319 zakw = EXP( zckw ) 320 zaksp1 = 10**(zaksp0) 321 zak1p = exp( zck1p ) 322 zak2p = exp( zck2p ) 323 zak3p = exp( zck3p ) 324 zaksi = exp( zcksi ) 325 zckf = zckf * total2SWS 326 327 ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970) 328 ! (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE 329 ! IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS 330 ! TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres IN 331 ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS 332 ! MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION 333 ! WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND 334 ! & GIESKES (1970), P. 1285-1286 (THE SMALL 335 ! FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE 336 ! SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)) 337 zcpexp = zpres / (rgas*ztkel) 338 zcpexp2 = zpres * zcpexp 339 340 ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE 341 ! CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968) 342 ! (CF. BROECKER ET AL., 1982) 343 344 zbuf1 = - ( devk10 + devk20 * ztc + devk30 * ztc * ztc ) 345 zbuf2 = 0.5 * ( devk40 + devk50 * ztc ) 346 ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 347 348 zbuf1 = - ( devk11 + devk21 * ztc + devk31 * ztc * ztc ) 349 zbuf2 = 0.5 * ( devk41 + devk51 * ztc ) 350 ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 351 352 zbuf1 = - ( devk12 + devk22 * ztc + devk32 * ztc * ztc ) 353 zbuf2 = 0.5 * ( devk42 + devk52 * ztc ) 354 akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 355 356 zbuf1 = - ( devk13 + devk23 * ztc + devk33 * ztc * ztc ) 357 zbuf2 = 0.5 * ( devk43 + devk53 * ztc ) 358 akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 359 360 zbuf1 = - ( devk14 + devk24 * ztc + devk34 * ztc * ztc ) 361 zbuf2 = 0.5 * ( devk44 + devk54 * ztc ) 362 aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 363 364 zbuf1 = - ( devk15 + devk25 * ztc + devk35 * ztc * ztc ) 365 zbuf2 = 0.5 * ( devk45 + devk55 * ztc ) 366 akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 367 368 zbuf1 = - ( devk17 + devk27 * ztc + devk37 * ztc * ztc ) 369 zbuf2 = 0.5 * ( devk47 + devk57 * ztc ) 370 ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 371 372 zbuf1 = - ( devk18 + devk28 * ztc + devk38 * ztc * ztc ) 373 zbuf2 = 0.5 * ( devk48 + devk58 * ztc ) 374 ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 375 376 zbuf1 = - ( devk19 + devk29 * ztc + devk39 * ztc * ztc ) 377 zbuf2 = 0.5 * ( devk49 + devk59 * ztc ) 378 ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 379 380 zbuf1 = - ( devk110 + devk210 * ztc + devk310 * ztc * ztc ) 381 zbuf2 = 0.5 * ( devk410 + devk510 * ztc ) 382 aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 383 384 ! CONVERT FROM DIFFERENT PH SCALES 385 total2free = 1.0/(1.0 + zst/aks3(ji,jj,jk)) 386 free2SWS = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk) 387 total2SWS = total2free * free2SWS 388 SWS2total = 1.0 / total2SWS 389 390 ! Convert to total scale 391 ak13(ji,jj,jk) = ak13(ji,jj,jk) * SWS2total 392 ak23(ji,jj,jk) = ak23(ji,jj,jk) * SWS2total 393 akb3(ji,jj,jk) = akb3(ji,jj,jk) * SWS2total 394 akw3(ji,jj,jk) = akw3(ji,jj,jk) * SWS2total 395 ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total 396 ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total 397 ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total 398 aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total 399 akf3(ji,jj,jk) = akf3(ji,jj,jk) / total2free 400 401 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 402 ! AS FUNCTION OF PRESSURE FOLLOWING MILLERO 403 ! (P. 1285) AND BERNER (1976) 404 zbuf1 = - ( devk16 + devk26 * ztc + devk36 * ztc * ztc ) 405 zbuf2 = 0.5 * ( devk46 + devk56 * ztc ) 406 aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 407 408 ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L] 409 borat(ji,jj,jk) = 0.0002414 * zcl / 10.811 410 sulfat(ji,jj,jk) = zst 411 fluorid(ji,jj,jk) = zft 412 413 ! Iron and SIO3 saturation concentration from ... 414 sio3eq(ji,jj,jk) = EXP( LOG( 10.) * ( 6.44 - 968. / ztkel ) ) * 1.e-6 415 fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 416 ! Liu and Millero (1999) only valid 5 - 50 degC 417 ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16 418 fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1) 419 fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 ) 420 fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 ) 421 fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 ) 422 fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 ) 423 END_3D 445 424 ! 446 425 IF( ln_timing ) CALL timing_stop('p4z_che')
Note: See TracChangeset
for help on using the changeset viewer.
