File indexing completed on 2024-06-04 05:10:25 UTC
view on githubraw file Latest commit 6518e4ec on 2024-06-03 19:02:13 UTC
c0d1c06c15 Matt*0001 #include "BLING_OPTIONS.h"
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 subroutine CALC_PCO2(
0012 I donewt,inewtonmax,ibrackmax,
0013 I t,s,diclocal,pt,sit,ta,
0014 I k1local,k2local,
0015 I k1plocal,k2plocal,k3plocal,
0016 I kslocal,kblocal,kwlocal,
0017 I ksilocal,kflocal,
0018 I k0local, fugflocal,
0019 I fflocal,btlocal,stlocal,ftlocal,
0020 U pHlocal,pCO2surfloc,
0021 I i,j,k,bi,bj,myIter,myThid )
0022
4ac06494d5 Matt*0023
0024
0025
0026
0027
0028
c0d1c06c15 Matt*0029
e0f9a7ba0b Matt*0030 IMPLICIT NONE
0031
c0d1c06c15 Matt*0032
0033 #include "SIZE.h"
0034 #include "DYNVARS.h"
0035 #include "EEPARAMS.h"
0036 #include "PARAMS.h"
0037 #include "GRID.h"
0038 #include "FFIELDS.h"
0039 #include "BLING_VARS.h"
0040
0041
0042
0043
0044
0045
0046
0047
ba0b047096 Mart*0048
c0d1c06c15 Matt*0049 INTEGER donewt
0050 INTEGER inewtonmax
0051 INTEGER ibrackmax
0052 _RL t, s, pt, sit, ta
0053 _RL pCO2surfloc, diclocal, pHlocal
0054 _RL fflocal, btlocal, stlocal, ftlocal
0055 _RL k1local, k2local
0056 _RL k1plocal, k2plocal, k3plocal
0057 _RL kslocal, kblocal, kwlocal, ksilocal, kflocal
0058 _RL k0local, fugflocal
0059 INTEGER i,j,k,bi,bj,myIter
0060 INTEGER myThid
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084 _RL phhi
0085 _RL phlo
0086 _RL c
0087 _RL a
0088 _RL a2
0089 _RL da
0090 _RL b
0091 _RL b2
0092 _RL db
0093 _RL fn
0094 _RL df
0095 _RL deltax
0096 _RL x
0097 _RL x2
0098 _RL x3
0099 _RL xmid
0100 _RL ftest
0101 _RL htotal
0102 _RL htotal2
0103 _RL co2star
0104 _RL phguess
0105 _RL fco2
0106 INTEGER inewton
0107 INTEGER ibrack
0108 INTEGER hstep
0109 _RL fni(3)
0110 _RL xlo
0111 _RL xhi
0112 _RL xguess
0113 _RL k123p
0114 _RL k12p
0115 _RL k12
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129 pt=pt*permil
0130 sit=sit*permil
0131 ta=ta*permil
0132 diclocal=diclocal*permil
0133
0134
0135
0136 phguess = phlocal
0137
0138
0139 phhi = 10.0
0140 phlo = 5.0
0141
0142 xguess = 10.0**(-phguess)
0143 xlo = 10.0**(-phhi)
0144 xhi = 10.0**(-phlo)
0145 xmid = (xlo + xhi)*0.5
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156 if( donewt .eq. 1)then
0157
0158
0159
0160 x = xguess
0161
0162
0163
0164 do inewton = 1, inewtonmax
0165
0166
0167 x2=x*x
0168 x3=x2*x
0169 k12 = k1local*k2local
0170 k12p = k1plocal*k2plocal
0171 k123p = k12p*k3plocal
0172 c = 1.0 + stlocal/kslocal
0173 a = x3 + k1plocal*x2 + k12p*x + k123p
0174 a2=a*a
0175 da = 3.0*x2 + 2.0*k1plocal*x + k12p
0176 b = x2 + k1local*x + k12
0177 b2=b*b
0178 db = 2.0*x + k1local
0179
0180
0181
0182
0183 fn = k1local*x*diclocal/b +
0184 & 2.0*diclocal*k12/b +
0185 & btlocal/(1.0 + x/kblocal) +
0186 & kwlocal/x +
0187 & pt*k12p*x/a +
0188 & 2.0*pt*k123p/a +
0189 & sit/(1.0 + x/ksilocal) -
0190 & x/c -
0191 & stlocal/(1.0 + kslocal/x/c) -
0192 & ftlocal/(1.0 + kflocal/x) -
0193 & pt*x3/a -
0194 & ta
0195
0196
0197
0198
0199
0200 df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -
0201 & 2.0*diclocal*k12*db/b2 -
0202 & btlocal/kblocal/(1.0+x/kblocal)**2. -
0203 & kwlocal/x2 +
0204 & (pt*k12p*(a - x*da))/a2 -
0205 & 2.0*pt*k123p*da/a2 -
0206 & sit/ksilocal/(1.0+x/ksilocal)**2. +
0207 & 1.0/c +
0208 & stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +
0209 & ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -
0210 & pt*x2*(3.0*a-x*da)/a2
0211
0212 deltax = - fn/df
0213
0214 x = x + deltax
0215
0216
0217
0218
0219
0220
0221 end do
0222
0223
0224 else
0225
0226
0227
0228
0229
0230 do ibrack = 1, ibrackmax
0231 do hstep = 1,3
0232 if(hstep .eq. 1)x = xhi
0233 if(hstep .eq. 2)x = xlo
0234 if(hstep .eq. 3)x = xmid
0235
0236
0237
0238 x2=x*x
0239 x3=x2*x
0240 k12 = k1local*k2local
0241 k12p = k1plocal*k2plocal
0242 k123p = k12p*k3plocal
0243 c = 1.0 + stlocal/kslocal
0244 a = x3 + k1plocal*x2 + k12p*x + k123p
0245 a2=a*a
0246 da = 3.0*x2 + 2.0*k1plocal*x + k12p
0247 b = x2 + k1local*x + k12
0248 b2=b*b
0249 db = 2.0*x + k1local
0250
0251 fn = k1local*x*diclocal/b +
0252 & 2.0*diclocal*k12/b +
0253 & btlocal/(1.0 + x/kblocal) +
0254 & kwlocal/x +
0255 & pt*k12p*x/a +
0256 & 2.0*pt*k123p/a +
0257 & sit/(1.0 + x/ksilocal) -
0258 & x/c -
0259 & stlocal/(1.0 + kslocal/x/c) -
0260 & ftlocal/(1.0 + kflocal/x) -
0261 & pt*x3/a -
0262 & ta
0263 fni(hstep) = fn
0264 end do
0265
0266 ftest = fni(1)/fni(3)
0267 if(ftest .gt. 0.0)then
0268 xhi = xmid
0269 else
0270 xlo = xmid
0271 end if
0272 xmid = (xlo + xhi)*0.5
0273
0274
0275
0276
0277
0278 end do
0279
0280 x = xmid
0281
0282
0283 end if
0284
0285
0286
0287
0288
0289 htotal = x
0290
0291
0292 htotal2=htotal*htotal
0293 co2star=diclocal*htotal2/(htotal2 + k1local*htotal
0294 & + k1local*k2local)
0295 phlocal=-log10(htotal)
0296
0297
0298
0299
0300
0301 fco2 = co2star / k0local
0302 pCO2surfloc = fco2/fugflocal
0303
0304
0305
0306
0307
0308
0309 pt=pt/permil
0310 sit=sit/permil
0311 ta=ta/permil
0312 diclocal=diclocal/permil
0313
0314 RETURN
0315 END
0316
0317
0318
0319
0320 subroutine CALC_PCO2_APPROX(
0321 I t,s,diclocal,pt,sit,ta,
0322 I k1local,k2local,
0323 I k1plocal,k2plocal,k3plocal,
0324 I kslocal,kblocal,kwlocal,
0325 I ksilocal,kflocal,
0326 I k0local, fugflocal,
0327 I fflocal,btlocal,stlocal,ftlocal,
0328 U pHlocal,pCO2surfloc,co3local,
0329 I i,j,k,bi,bj,myIter,myThid )
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
e0f9a7ba0b Matt*0344 IMPLICIT NONE
0345
c0d1c06c15 Matt*0346
0347 #include "SIZE.h"
0348 #include "DYNVARS.h"
0349 #include "EEPARAMS.h"
0350 #include "PARAMS.h"
0351 #include "GRID.h"
0352 #include "FFIELDS.h"
0353 #include "BLING_VARS.h"
0354
0355
0356
0357
0358
0359
0360
0361
ba0b047096 Mart*0362
c0d1c06c15 Matt*0363 _RL t, s, pt, sit, ta
0364 _RL pCO2surfloc, diclocal, pHlocal
0365 _RL fflocal, btlocal, stlocal, ftlocal
0366 _RL k1local, k2local
0367 _RL k1plocal, k2plocal, k3plocal
0368 _RL kslocal, kblocal, kwlocal, ksilocal, kflocal
0369 _RL k0local, fugflocal
0370 _RL co3local
0371 INTEGER i,j,k,bi,bj,myIter
0372 INTEGER myThid
0373
0374
0375
0376 _RL phguess
0377 _RL cag
0378 _RL bohg
0379 _RL hguess
0380 _RL stuff
0381 _RL gamm
0382 _RL hnew
0383 _RL co2s
0384 _RL h3po4g, h2po4g, hpo4g, po4g
0385 _RL siooh3g
0386 _RL fco2
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396 pt=pt*permil
0397 sit=sit*permil
0398 ta=ta*permil
0399 diclocal=diclocal*permil
0400
0401
0402
0403 phguess = phlocal
0404
0405
0406
0407 hguess = 10.0**(-phguess)
0408
0409 bohg = btlocal*kblocal/(hguess+kblocal)
0410
0411
0412
0413 stuff = hguess*hguess*hguess
0414 & + (k1plocal*hguess*hguess)
0415 & + (k1plocal*k2plocal*hguess)
0416 & + (k1plocal*k2plocal*k3plocal)
0417 h3po4g = (pt*hguess*hguess*hguess) / stuff
0418 h2po4g = (pt*k1plocal*hguess*hguess) / stuff
0419 hpo4g = (pt*k1plocal*k2plocal*hguess) / stuff
0420 po4g = (pt*k1plocal*k2plocal*k3plocal) / stuff
0421
0422
0423
0424 siooh3g = sit*ksilocal / (ksilocal + hguess)
0425
1c4959ef9f Matt*0426
0427
0428 IF (siooh3g.GT.0.2 _d 0*ta) siooh3g = 0.2 _d 0*ta
0429
c0d1c06c15 Matt*0430
0431 cag = ta - bohg - (kwlocal/hguess) + hguess
0432 & - hpo4g - 2.0 _d 0*po4g + h3po4g
0433 & - siooh3g
e0f9a7ba0b Matt*0434
0435 IF ( cag .LT. (0.1 _d 0*ta) ) THEN
0436 print *,'CMMcag1', cag,i,j,k,xC(i,j,1,1),yC(i,j,1,1)
0437
0438 print *,'CMMcag2',ta,bohg,(kwlocal/hguess),hguess,hpo4g,
0439 & 2.0 _d 0*po4g,h3po4g, siooh3g
0440 cag = 0.1 _d 0*ta
0441 ENDIF
0442
c0d1c06c15 Matt*0443
0444
0445 gamm = diclocal/cag
0446 stuff = (1.0 _d 0-gamm)*(1.0 _d 0-gamm)*k1local*k1local
0447 & - 4.0 _d 0*k1local*k2local*(1.0 _d 0-2.0 _d 0*gamm)
0448 hnew = 0.5 _d 0*( (gamm-1.0 _d 0)*k1local + sqrt(stuff) )
e0f9a7ba0b Matt*0449
0450 IF ( hnew .GT. (0.00001 _d 0) ) THEN
0451
0452 hnew = 0.00000001 _d 0
0453 ENDIF
0454 IF ( hnew .LT. (0.0000000001 _d 0) ) THEN
0455
0456 hnew = 0.00000001 _d 0
0457 ENDIF
0458
c0d1c06c15 Matt*0459
0460 co2s = diclocal/
0461 & (1.0 _d 0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))
0462
0463 phlocal = -log10(hnew)
0464
0465
0466
0467
0468 co3local = k1local*k2local*diclocal /
0469 & (hnew*hnew + k1local*hnew + k1local*k2local)
0470
0471
0472
0473 fco2 = co2s/k0local
0474 pco2surfloc = fco2/fugflocal
0475
0476
0477
0478 pt=pt/permil
0479 sit=sit/permil
0480 ta=ta/permil
0481 diclocal=diclocal/permil
0482
0483 RETURN
0484 END
0485
0486
0487
0488 subroutine CARBON_COEFFS(
0489 I ttemp,stemp,
0490 I bi,bj,iMin,iMax,jMin,jMax,myThid)
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
ba0b047096 Mart*0507
c0d1c06c15 Matt*0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
e0f9a7ba0b Matt*0522 IMPLICIT NONE
0523
c0d1c06c15 Matt*0524
0525 #include "SIZE.h"
0526 #include "DYNVARS.h"
0527 #include "EEPARAMS.h"
0528 #include "PARAMS.h"
0529 #include "GRID.h"
0530 #include "FFIELDS.h"
aba82ffbc6 Matt*0531 #include "BLING_VARS.h"
4ac06494d5 Matt*0532
c0d1c06c15 Matt*0533
0534
0535
0536
0537
0538 _RL ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0539 _RL stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0540 INTEGER bi,bj,iMin,iMax,jMin,jMax
0541 INTEGER myThid
0542
0543
0544
0545 _RL t
0546 _RL s
0547 _RL tk
0548 _RL tk100
0549 _RL tk1002
0550 _RL dlogtk
0551 _RL sqrtis
0552 _RL sqrts
0553 _RL s15
0554 _RL scl
0555 _RL s2
0556 _RL invtk
0557 _RL is
0558 _RL is2
0559
0560 _RL P1atm
0561 _RL Rgas
0562 _RL RT
0563 _RL delta
0564 _RL B1
0565 _RL B
e0f9a7ba0b Matt*0566 #ifdef CARBONCHEM_TOTALPHSCALE
0567
0568 _RL total2free
0569 _RL free2sw
0570 _RL total2sw
0571 #endif
c0d1c06c15 Matt*0572 INTEGER i
0573 INTEGER j
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586 do i=imin,imax
0587 do j=jmin,jmax
0588 if (hFacC(i,j,1,bi,bj).gt.0. _d 0) then
0589 t = ttemp(i,j)
0590 s = stemp(i,j)
e0f9a7ba0b Matt*0591
0592
c0d1c06c15 Matt*0593 tk = 273.15 _d 0 + t
0594 tk100 = tk/100. _d 0
0595 tk1002=tk100*tk100
0596 invtk=1.0 _d 0/tk
0597 dlogtk=log(tk)
e0f9a7ba0b Matt*0598
c0d1c06c15 Matt*0599 is=19.924 _d 0*s/(1000. _d 0-1.005 _d 0*s)
0600 is2=is*is
0601 sqrtis=sqrt(is)
e0f9a7ba0b Matt*0602
c0d1c06c15 Matt*0603 s2=s*s
0604 sqrts=sqrt(s)
0605 s15=s**1.5 _d 0
0606 scl=s/1.80655 _d 0
0607
0608
0609
0610
0611
0612 P1atm = 1.01325 _d 0
0613 Rgas = 83.1451 _d 0
0614 RT = Rgas*tk
0615 delta = (57.7 _d 0 - 0.118 _d 0*tk)
0616 B1 = -1636.75 _d 0 + 12.0408 _d 0*tk - 0.0327957 _d 0*tk*tk
0617 B = B1 + 3.16528 _d 0*tk*tk*tk*(0.00001 _d 0)
4ac06494d5 Matt*0618 fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) *
e0f9a7ba0b Matt*0619 & AtmosP(i,j,bi,bj) * P1atm / RT)
4ac06494d5 Matt*0620
c0d1c06c15 Matt*0621
0622
0623
0624 ff(i,j,bi,bj) = exp(-162.8301 _d 0 + 218.2968 _d 0/tk100 +
0625 & 90.9241 _d 0*log(tk100) - 1.47696 _d 0*tk1002 +
0626 & s * (.025695 _d 0 - .025225 _d 0*tk100 +
0627 & 0.0049867 _d 0*tk1002))
0628
0629
e0f9a7ba0b Matt*0630
c0d1c06c15 Matt*0631 ak0(i,j,bi,bj) = exp(93.4517 _d 0/tk100 - 60.2409 _d 0 +
0632 & 23.3585 _d 0 * log(tk100) +
0633 & s * (0.023517 _d 0 - 0.023656 _d 0*tk100 +
0634 & 0.0047036 _d 0*tk1002))
0635
0636
0637
0638
e0f9a7ba0b Matt*0639
c0d1c06c15 Matt*0640 ak1(i,j,bi,bj)=10.**(-1. _d 0*(3670.7 _d 0*invtk -
0641 & 62.008 _d 0 + 9.7944 _d 0*dlogtk -
0642 & 0.0118 _d 0 * s + 0.000116 _d 0*s2))
b4daa24319 Shre*0643 ak2(i,j,bi,bj)=10.**(-1. _d 0*(1394.7 _d 0*invtk +4.777 _d 0
0644 & -0.0184 _d 0*s + 0.000118 _d 0*s2))
c0d1c06c15 Matt*0645
0646
0647
e0f9a7ba0b Matt*0648
c0d1c06c15 Matt*0649 akb(i,j,bi,bj)=exp((-8966.90 _d 0- 2890.53 _d 0*sqrts -
0650 & 77.942 _d 0*s + 1.728 _d 0*s15 - 0.0996 _d 0*s2)*invtk +
0651 & (148.0248 _d 0 + 137.1942 _d 0*sqrts + 1.62142 _d 0*s) +
0652 & (-24.4344 _d 0 - 25.085 _d 0*sqrts - 0.2474 _d 0*s) *
0653 & dlogtk + 0.053105 _d 0*sqrts*tk)
0654
0655
0656
e0f9a7ba0b Matt*0657
0658
0659
c0d1c06c15 Matt*0660 ak1p(i,j,bi,bj) = exp(-4576.752 _d 0*invtk + 115.525 _d 0 -
0661 & 18.453 _d 0*dlogtk +
0662 & (-106.736 _d 0*invtk + 0.69171 _d 0)*sqrts +
0663 & (-0.65643 _d 0*invtk - 0.01844 _d 0)*s)
0664
0665
0666
e0f9a7ba0b Matt*0667
0668
0669
c0d1c06c15 Matt*0670 ak2p(i,j,bi,bj) = exp(-8814.715 _d 0*invtk + 172.0883 _d 0 -
0671 & 27.927 _d 0*dlogtk +
0672 & (-160.340 _d 0*invtk + 1.3566 _d 0) * sqrts +
0673 & (0.37335 _d 0*invtk - 0.05778 _d 0) * s)
0674
0675
0676
e0f9a7ba0b Matt*0677
0678
0679
c0d1c06c15 Matt*0680 ak3p(i,j,bi,bj) = exp(-3070.75 _d 0*invtk - 18.141 _d 0 +
0681 & (17.27039 _d 0*invtk + 2.81197 _d 0) *
0682 & sqrts + (-44.99486 _d 0*invtk - 0.09984 _d 0) * s)
0683
0684
0685
e0f9a7ba0b Matt*0686
0687
0688
0689
c0d1c06c15 Matt*0690 aksi(i,j,bi,bj) = exp(-8904.2 _d 0*invtk + 117.385 _d 0 -
0691 & 19.334 _d 0*dlogtk +
0692 & (-458.79 _d 0*invtk + 3.5913 _d 0) * sqrtis +
0693 & (188.74 _d 0*invtk - 1.5998 _d 0) * is +
0694 & (-12.1652 _d 0*invtk + 0.07871 _d 0) * is2 +
0695 & log(1.0 _d 0-0.001005 _d 0*s))
0696
0697
0698
e0f9a7ba0b Matt*0699
0700
0701
c0d1c06c15 Matt*0702 akw(i,j,bi,bj) = exp(-13847.26 _d 0*invtk + 148.9652 _d 0 -
0703 & 23.6521 _d 0*dlogtk +
0704 & (118.67 _d 0*invtk - 5.977 _d 0 + 1.0495 _d 0 * dlogtk)
0705 & * sqrts - 0.01615 _d 0 * s)
0706
0707
0708
e0f9a7ba0b Matt*0709
0710
c0d1c06c15 Matt*0711 aks(i,j,bi,bj)=exp(-4276.1 _d 0*invtk + 141.328 _d 0 -
0712 & 23.093 _d 0*dlogtk +
0713 & (-13856. _d 0*invtk + 324.57 _d 0 - 47.986 _d 0*dlogtk)*sqrtis+
0714 & (35474. _d 0*invtk - 771.54 _d 0 + 114.723 _d 0*dlogtk)*is -
0715 & 2698. _d 0*invtk*is**1.5 _d 0 + 1776. _d 0*invtk*is2 +
0716 & log(1.0 _d 0 - 0.001005 _d 0*s))
0717
0718
e0f9a7ba0b Matt*0719
0720
c0d1c06c15 Matt*0721 akf(i,j,bi,bj)=exp(1590.2 _d 0*invtk - 12.641 _d 0 +
0722 & 1.525 _d 0*sqrtis + log(1.0 _d 0 - 0.001005 _d 0*s) +
0723 & log(1.0 _d 0 + (0.1400 _d 0/96.062 _d 0)*(scl)/aks(i,j,bi,bj)))
0724
0725
0726
0727 bt(i,j,bi,bj) = 0.000232 _d 0 * scl/10.811 _d 0
0728
0729 st(i,j,bi,bj) = 0.14 _d 0 * scl/96.062 _d 0
0730
0731 ft(i,j,bi,bj) = 0.000067 _d 0 * scl/18.9984 _d 0
0732
e0f9a7ba0b Matt*0733 #ifdef CARBONCHEM_TOTALPHSCALE
0734
0735
0736 total2free = 1. _d 0/
0737 & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
0738 free2sw = 1. _d 0
0739 & + st(i,j,bi,bj)/ aks(i,j,bi,bj)
0740 & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free)
0741 total2sw = total2free * free2sw
0742 ak1(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw
0743 ak2(i,j,bi,bj) = ak2(i,j,bi,bj)/total2sw
0744 aks(i,j,bi,bj) = aks(i,j,bi,bj)/total2free
0745 #endif
c0d1c06c15 Matt*0746 else
0747
0748 fugf(i,j,bi,bj)=0. _d 0
0749 ff(i,j,bi,bj)=0. _d 0
0750 ak0(i,j,bi,bj)= 0. _d 0
0751 ak1(i,j,bi,bj)= 0. _d 0
0752 ak2(i,j,bi,bj)= 0. _d 0
0753 akb(i,j,bi,bj)= 0. _d 0
0754 ak1p(i,j,bi,bj) = 0. _d 0
0755 ak2p(i,j,bi,bj) = 0. _d 0
0756 ak3p(i,j,bi,bj) = 0. _d 0
0757 aksi(i,j,bi,bj) = 0. _d 0
0758 akw(i,j,bi,bj) = 0. _d 0
0759 aks(i,j,bi,bj)= 0. _d 0
0760 akf(i,j,bi,bj)= 0. _d 0
0761 bt(i,j,bi,bj) = 0. _d 0
0762 st(i,j,bi,bj) = 0. _d 0
0763 ft(i,j,bi,bj) = 0. _d 0
0764 endif
0765 end do
0766 end do
0767
0768 RETURN
0769 END
0770
0771
0772
0773
0774 subroutine CARBON_COEFFS_PRESSURE_DEP(
0775 I ttemp,stemp,
0776 I bi,bj,iMin,iMax,jMin,jMax,
0777 I Klevel,myThid)
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
ba0b047096 Mart*0796
c0d1c06c15 Matt*0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
e0f9a7ba0b Matt*0814 IMPLICIT NONE
0815
c0d1c06c15 Matt*0816
0817 #include "SIZE.h"
0818 #include "DYNVARS.h"
0819 #include "EEPARAMS.h"
0820 #include "PARAMS.h"
0821 #include "GRID.h"
0822 #include "FFIELDS.h"
aba82ffbc6 Matt*0823 #include "BLING_VARS.h"
0824
c0d1c06c15 Matt*0825
0826
0827
0828
0829
0830 _RL ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0831 _RL stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0832 INTEGER bi,bj,iMin,iMax,jMin,jMax
0833
0834 INTEGER Klevel
0835 INTEGER myThid
0836
0837
0838
0839 _RL t
0840 _RL s
0841 _RL tk
0842 _RL tk100
0843 _RL tk1002
0844 _RL dlogtk
0845 _RL sqrtis
0846 _RL sqrts
0847 _RL s15
0848 _RL scl
0849 _RL s2
0850 _RL invtk
0851 _RL is
0852 _RL is2
0853
0854 _RL P1atm
0855 _RL Rgas
0856 _RL RT
0857 _RL delta
0858 _RL B1
0859 _RL B
0860 INTEGER i
0861 INTEGER j
0862 INTEGER k
0863 _RL bdepth
0864 _RL cdepth
0865 _RL pressc
0866 _RL Ksp_T_Calc
0867 _RL Ksp_T_Arag
0868 _RL xvalue
0869 _RL zdum
0870 _RL tmpa1
0871 _RL tmpa2
0872 _RL tmpa3
0873 _RL logKspc
0874 _RL logKspa
0875 _RL dv
0876 _RL dk
0877 _RL pfactor
0878 _RL bigR
e0f9a7ba0b Matt*0879 #ifdef CARBONCHEM_TOTALPHSCALE
0880
0881 _RL total2free_surf
0882 _RL free2sw_surf
0883 _RL total2sw_surf
0884 _RL total2free
0885 _RL free2sw
0886 _RL total2sw
0887 #endif
c0d1c06c15 Matt*0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
0900
0901
0902
0903
0904
0905 bdepth = 0.0d0
0906 cdepth = 0.0d0
0907 pressc = 1.0d0
0908 do k = 1,Klevel
0909 cdepth = bdepth + 0.5d0*drF(k)
0910 bdepth = bdepth + drF(k)
0911 pressc = 1.0d0 + 0.1d0*cdepth
0912 end do
0913
6518e4ecc5 aver*0914 bigR = 83.145 _d 0
0915 P1atm = 1.01325 _d 0
0916 Rgas = 83.1451 _d 0
0917
c0d1c06c15 Matt*0918 do i=imin,imax
0919 do j=jmin,jmax
0920 if (hFacC(i,j,Klevel,bi,bj).gt.0.d0) then
0921 t = ttemp(i,j)
0922 s = stemp(i,j)
e0f9a7ba0b Matt*0923
0924
c0d1c06c15 Matt*0925 tk = 273.15 + t
0926 tk100 = tk/100.0
0927 tk1002=tk100*tk100
0928 invtk=1.0/tk
0929 dlogtk=log(tk)
e0f9a7ba0b Matt*0930
c0d1c06c15 Matt*0931 is=19.924*s/(1000.-1.005*s)
0932 is2=is*is
0933 sqrtis=sqrt(is)
e0f9a7ba0b Matt*0934
c0d1c06c15 Matt*0935 s2=s*s
0936 sqrts=sqrt(s)
0937 s15=s**1.5
0938 scl=s/1.80655
0939
0940
0941
0942
0943
0944
0945
0946 RT = Rgas*tk
0947 delta = (57.7 _d 0 - 0.118 _d 0*tk)
0948 B1 = -1636.75 _d 0 + 12.0408 _d 0*tk - 0.0327957 _d 0*tk*tk
0949 B = B1 + 3.16528 _d 0*tk*tk*tk*(0.00001 _d 0)
e0f9a7ba0b Matt*0950 fugf(i,j,bi,bj) = exp( (B+2. _d 0*delta) *
0951 & AtmosP(i,j,bi,bj) * P1atm / RT)
c0d1c06c15 Matt*0952
0953
0954
0955 ff(i,j,bi,bj) = exp(-162.8301 + 218.2968/tk100 +
0956 & 90.9241*log(tk100) - 1.47696*tk1002 +
0957 & s * (.025695 - .025225*tk100 +
0958 & 0.0049867*tk1002))
0959
0960
e0f9a7ba0b Matt*0961
c0d1c06c15 Matt*0962 ak0(i,j,bi,bj) = exp(93.4517/tk100 - 60.2409 +
0963 & 23.3585 * log(tk100) +
0964 & s * (0.023517 - 0.023656*tk100 +
0965 & 0.0047036*tk1002))
0966
0967
0968
0969
e0f9a7ba0b Matt*0970
b4daa24319 Shre*0971 ak1(i,j,bi,bj)=10.**(-1.*(3670.7*invtk -
c0d1c06c15 Matt*0972 & 62.008 + 9.7944*dlogtk -
0973 & 0.0118 * s + 0.000116*s2))
b4daa24319 Shre*0974 ak2(i,j,bi,bj)=10.**(-1.*(1394.7*invtk + 4.777 -
c0d1c06c15 Matt*0975 & 0.0184*s + 0.000118*s2))
0976
0977
0978
0979
0980 ak1(i,j,bi,bj) = ak1(i,j,bi,bj)*
0981 & exp( (24.2-0.085*t)*(pressc-1.0)/(83.143*tk) )
0982
0983
0984
0985
0986
0987 ak2(i,j,bi,bj) = ak2(i,j,bi,bj)*
0988 & exp( (16.4-0.040*t)*(pressc-1.0)/(83.143*tk) )
0989
0990
0991
e0f9a7ba0b Matt*0992
c0d1c06c15 Matt*0993 akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s +
0994 & 1.728*s15 - 0.0996*s2)*invtk +
0995 & (148.0248 + 137.1942*sqrts + 1.62142*s) +
0996 & (-24.4344 - 25.085*sqrts - 0.2474*s) *
0997 & dlogtk + 0.053105*sqrts*tk)
e0f9a7ba0b Matt*0998 #ifdef CARBONCHEM_TOTALPHSCALE
0999
1000
c0d1c06c15 Matt*1001
1002
1003
1004
e0f9a7ba0b Matt*1005
1006
1007
1008
1009
1010 #else
1011
1012
1013
1014
1015 dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
1016 dk = -2.84 _d -3
c0d1c06c15 Matt*1017 pfactor = - (dv/(bigR*tk))*pressc
e0f9a7ba0b Matt*1018 & + (0.5 _d 0*dk/(bigR*tk))*pressc*pressc
c0d1c06c15 Matt*1019 akb(i,j,bi,bj) = akb(i,j,bi,bj)*exp(pfactor)
e0f9a7ba0b Matt*1020 #endif
c0d1c06c15 Matt*1021
1022
1023
e0f9a7ba0b Matt*1024
1025
1026
c0d1c06c15 Matt*1027 ak1p(i,j,bi,bj) = exp(-4576.752*invtk + 115.525 -
1028 & 18.453*dlogtk +
1029 & (-106.736*invtk + 0.69171)*sqrts +
1030 & (-0.65643*invtk - 0.01844)*s)
1031
1032
1033
e0f9a7ba0b Matt*1034
1035
1036
c0d1c06c15 Matt*1037 ak2p(i,j,bi,bj) = exp(-8814.715*invtk + 172.0883 -
1038 & 27.927*dlogtk +
1039 & (-160.340*invtk + 1.3566) * sqrts +
1040 & (0.37335*invtk - 0.05778) * s)
1041
1042
1043
e0f9a7ba0b Matt*1044
1045
1046
c0d1c06c15 Matt*1047 ak3p(i,j,bi,bj) = exp(-3070.75*invtk - 18.141 +
1048 & (17.27039*invtk + 2.81197) *
1049 & sqrts + (-44.99486*invtk - 0.09984) * s)
1050
1051
1052
e0f9a7ba0b Matt*1053
1054
1055
1056
c0d1c06c15 Matt*1057 aksi(i,j,bi,bj) = exp(-8904.2*invtk + 117.385 -
1058 & 19.334*dlogtk +
1059 & (-458.79*invtk + 3.5913) * sqrtis +
1060 & (188.74*invtk - 1.5998) * is +
1061 & (-12.1652*invtk + 0.07871) * is2 +
1062 & log(1.0-0.001005*s))
1063
1064
1065
e0f9a7ba0b Matt*1066
1067
1068
c0d1c06c15 Matt*1069 akw(i,j,bi,bj) = exp(-13847.26*invtk + 148.9652 -
1070 & 23.6521*dlogtk +
1071 & (118.67*invtk - 5.977 + 1.0495 * dlogtk) *
1072 & sqrts - 0.01615 * s)
1073
1074
1075
e0f9a7ba0b Matt*1076
1077
c0d1c06c15 Matt*1078 aks(i,j,bi,bj)=exp(-4276.1*invtk + 141.328 -
1079 & 23.093*dlogtk +
1080 & (-13856*invtk + 324.57 - 47.986*dlogtk)*sqrtis +
1081 & (35474*invtk - 771.54 + 114.723*dlogtk)*is -
1082 & 2698*invtk*is**1.5 + 1776*invtk*is2 +
1083 & log(1.0 - 0.001005*s))
1084
1085
e0f9a7ba0b Matt*1086
1087
c0d1c06c15 Matt*1088 akf(i,j,bi,bj)=exp(1590.2*invtk - 12.641 + 1.525*sqrtis +
1089 & log(1.0 - 0.001005*s) +
1090 & log(1.0 + (0.1400/96.062)*(scl)/aks(i,j,bi,bj)))
1091
1092
1093
1094 bt(i,j,bi,bj) = 0.000232 * scl/10.811
1095
1096 st(i,j,bi,bj) = 0.14 * scl/96.062
1097
1098 ft(i,j,bi,bj) = 0.000067 * scl/18.9984
1099
e0f9a7ba0b Matt*1100
1101 #ifdef CARBONCHEM_TOTALPHSCALE
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111 total2free_surf = 1. _d 0/
1112 & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
1113
1114 free2sw_surf = 1. _d 0
1115 & + st(i,j,bi,bj)/ aks(i,j,bi,bj)
1116 & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free_surf)
1117
1118 total2sw_surf = total2free_surf * free2sw_surf
1119
1120
1121
1122 dv=-18.03 _d 0 + 0.0466 _d 0*t + 0.316 _d -3 *t*t
1123 dk=-4.53 _d -3 + 0.09 _d -3 *t
1124 pfactor = -(dv/(bigR*tk))*pressc
1125 & +(0.5*dk/(bigR*tk))*pressc*pressc
1126 aks(i,j,bi,bj) = aks(i,j,bi,bj)*exp(pfactor)
1127
1128 total2free = 1. _d 0/
1129 & (1. _d 0 + st(i,j,bi,bj)/aks(i,j,bi,bj))
1130
1131 free2sw = 1. _d 0
1132 & + st(i,j,bi,bj)/ aks(i,j,bi,bj)
1133
1134 aks(i,j,bi,bj) = aks(i,j,bi,bj)/total2free
1135
1136
1137
1138 akf(i,j,bi,bj) = akf(i,j,bi,bj) * total2free_surf
1139 dv = -9.78 _d 0 -0.0090 _d 0 *t -0.942 _d -3 *t*t
1140 dk = -3.91 _d -3 + 0.054 _d -3 *t
1141 pfactor = -(dv/(bigR*tk))*pressc
1142 & +(0.5*dk/(bigR*tk))*pressc*pressc
1143 akf(i,j,bi,bj) = akf(i,j,bi,bj) * exp(pfactor)
1144 akf(i,j,bi,bj) = akf(i,j,bi,bj)/total2free
1145
1146 free2sw = free2sw
1147 & + ft(i,j,bi,bj)/(akf(i,j,bi,bj)*total2free)
1148 total2sw = total2free * free2sw
1149
1150
1151
1152 ak1(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw
1153 ak2(i,j,bi,bj) = ak1(i,j,bi,bj)/total2sw
1154
1155
1156
1157 dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
1158 dk = -2.84 _d -3
1159 pfactor = - (dv/(bigR*tk))*pressc
1160 & + (0.5*dk/(bigR*tk))*pressc*pressc
1161 akb(i,j,bi,bj) = total2sw_surf*akb(i,j,bi,bj)*exp(pfactor)
1162 akb(i,j,bi,bj) = akb(i,j,bi,bj)/total2sw
1163
1164
1165
1166 dv = -14.51 _d 0 + 0.1211 _d 0*t -0.321 _d -3*t*t
1167 dk = -2.67 _d -3 + 0.0427 _d -3*t
1168 pfactor = - (dv/(bigR*tk))*pressc
1169 & + (0.5*dk/(bigR*tk))*pressc*pressc
1170 ak1p(i,j,bi,bj) = total2sw_surf*ak1p(i,j,bi,bj)*exp(pfactor)
1171 ak1p(i,j,bi,bj) = ak1p(i,j,bi,bj)/total2sw
1172
1173
1174
1175 dv = -23.12 _d 0 + 0.1758 _d 0*t -2.647 _d -3*t*t
1176 dk = -5.15 _d -3 + 0.09 _d -3*t
1177 pfactor = - (dv/(bigR*tk))*pressc
1178 & + (0.5*dk/(bigR*tk))*pressc*pressc
1179 ak2p(i,j,bi,bj) = total2sw_surf*ak2p(i,j,bi,bj)*exp(pfactor)
1180 ak2p(i,j,bi,bj) = ak2p(i,j,bi,bj)/total2sw
1181
1182
1183
1184 dv = -26.57 _d 0 + 0.2020 _d 0*t -3.042 _d -3*t*t
1185 dk = -4.08 _d -3 + 0.0714 _d -3*t
1186 pfactor = - (dv/(bigR*tk))*pressc
1187 & + (0.5*dk/(bigR*tk))*pressc*pressc
1188 ak3p(i,j,bi,bj) = total2sw_surf*ak3p(i,j,bi,bj)*exp(pfactor)
1189 ak3p(i,j,bi,bj) = ak3p(i,j,bi,bj)/total2sw
1190
1191
1192
1193 dv = -20.02 _d 0 + 0.1119 _d 0*t -1.409 _d -3*t*t
1194 dk = -5.13 _d -3 + 0.0794 _d -3 *t
1195 pfactor = - (dv/(bigR*tk))*pressc
1196 & + (0.5*dk/(bigR*tk))*pressc*pressc
1197 akw(i,j,bi,bj) = total2sw_surf*akw(i,j,bi,bj)*exp(pfactor)
1198 akw(i,j,bi,bj) = akw(i,j,bi,bj)/total2sw
1199
1200
1201
1202
1203 dv = -29.48 _d 0 + 0.1622 _d 0*t + 2.608 _d -3*t*t
1204 dk = -2.84 _d -3
1205 pfactor = - (dv/(bigR*tk))*pressc
1206 & + (0.5*dk/(bigR*tk))*pressc*pressc
1207 aksi(i,j,bi,bj) = total2sw_surf*aksi(i,j,bi,bj)*exp(pfactor)
1208 aksi(i,j,bi,bj) = aksi(i,j,bi,bj)/total2sw
1209 #endif
1210
1211
c0d1c06c15 Matt*1212
e0f9a7ba0b Matt*1213
c0d1c06c15 Matt*1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227 tmpa1 = - 171.9065 - (0.077993*tk) + (2839.319/tk)
1228 & + (71.595*log10(tk))
1229 tmpa2 = +(-0.77712 + (0.0028426*tk) + (178.34/tk) )*sqrts
1230 tmpa3 = -(0.07711*s) + (0.0041249*s15)
1231 logKspc = tmpa1 + tmpa2 + tmpa3
1232 Ksp_T_Calc = 10.0**logKspc
1233
1234 tmpa1 = - 171.945 - (0.077993*tk) + (2903.293/tk)
1235 & + (71.595*log10(tk))
1236 tmpa2 = +(-0.068393 + (0.0017276*tk) + (88.135/tk) )*sqrts
1237 tmpa3 = -(0.10018*s) + (0.0059415*s15)
1238 logKspa = tmpa1 + tmpa2 + tmpa3
1239 Ksp_T_Arag = 10.0**logKspa
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256 zdum = (pressc*10.0d0 - 10.0d0)/10.0d0
1257 xvalue = ( (48.8d0 - 0.53d0*t)*zdum
1258 & + (-0.00588d0 + 0.0001845d0*t)*zdum*zdum)
1259 & / (188.93d0*(t + 273.15d0))
1260
b4daa24319 Shre*1261 Ksp_TP_Calc(i,j,bi,bj) = Ksp_T_Calc*10.**(xvalue)
1262 Ksp_TP_Arag(i,j,bi,bj) = Ksp_T_Arag*10.**(xvalue)
c0d1c06c15 Matt*1263
1264
1265 else
e0f9a7ba0b Matt*1266 fugf(i,j,bi,bj) = 0.d0
1267 ff(i,j,bi,bj) = 0.d0
c0d1c06c15 Matt*1268 ak0(i,j,bi,bj)= 0.d0
1269 ak1(i,j,bi,bj)= 0.d0
1270 ak2(i,j,bi,bj)= 0.d0
1271 akb(i,j,bi,bj)= 0.d0
1272 ak1p(i,j,bi,bj) = 0.d0
1273 ak2p(i,j,bi,bj) = 0.d0
1274 ak3p(i,j,bi,bj) = 0.d0
1275 aksi(i,j,bi,bj) = 0.d0
1276 akw(i,j,bi,bj) = 0.d0
1277 aks(i,j,bi,bj)= 0.d0
1278 akf(i,j,bi,bj)= 0.d0
1279 bt(i,j,bi,bj) = 0.d0
1280 st(i,j,bi,bj) = 0.d0
1281 ft(i,j,bi,bj) = 0.d0
1282 Ksp_TP_Calc(i,j,bi,bj) = 0.d0
1283 Ksp_TP_Arag(i,j,bi,bj) = 0.d0
1284 endif
1285 end do
1286 end do
1287
1288 return
1289 end