/[MITgcm]/MITgcm/pkg/dic/carbon_chem.F
ViewVC logotype

Contents of /MITgcm/pkg/dic/carbon_chem.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.3 - (show annotations) (download)
Thu Oct 9 04:19:19 2003 UTC (20 years, 8 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52e_pre, checkpoint51o_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint51l_post, checkpoint51q_post, hrcube_1, branch-netcdf, checkpoint52d_pre, checkpoint51r_post, checkpoint52b_pre, checkpoint51o_post, checkpoint51p_post, checkpoint52a_pre, checkpoint51i_post, checkpoint52, checkpoint52d_post, checkpoint52a_post, checkpoint52b_post, checkpoint52f_post, checkpoint52c_post, checkpoint51l_pre, ecco_c52_e35, checkpoint52i_post, checkpoint51t_post, checkpoint51n_post, checkpoint52i_pre, checkpoint51u_post, checkpoint52h_pre, checkpoint52f_pre, hrcube_2, hrcube_3, checkpoint51m_post, checkpoint51s_post
Branch point for: branch-nonh, tg2-branch, checkpoint51n_branch, netcdf-sm0
Changes since 1.2: +4 -3 lines 
 o first check-in for the "branch-genmake2" merge
 o verification suite as run on shelley (gcc 3.2.2):

Wed Oct  8 23:42:29 EDT 2003
                T           S           U           V
G D M    c        m  s        m  s        m  s        m  s
E p a R  g  m  m  e  .  m  m  e  .  m  m  e  .  m  m  e  .
N n k u  2  i  a  a  d  i  a  a  d  i  a  a  d  i  a  a  d
2 d e n  d  n  x  n  .  n  x  n  .  n  x  n  .  n  x  n  .

OPTFILE=NONE

Y Y Y Y 13 16 16 16  0 16 16 16 16 16 16 16 16 13 12  0  0 pass  adjustment.128x64x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16  0  0 16 16  0  0 pass  adjustment.cs-32x32x1
Y Y Y Y 16 16 16 16  0 16 16 16 16 16 16 22  0 16 16 22  0 pass  adjust_nlfs.cs-32x32x1
Y Y Y Y -- 13 13 16 16 13 13 13 13 16 16 16 16 16 16 16 16 N/O   advect_cs
Y Y Y Y -- 22 16 16 16 16 16 16 13 16 16 16 16 16 16 16 16 N/O   advect_xy
Y Y Y Y -- 13 16 13 16 16 16 16 16 16 16 22 16 16 16 16 16 N/O   advect_xz
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  aim.5l_cs
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 16 16 16 16 13 16 pass  aim.5l_Equatorial_Channel
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 13 16 16 13 13 16 pass  aim.5l_LatLon
Y Y Y Y 13 16 16 16 16 16 16 16 16 16 13 12 13 13 16 13 16 pass  exp0
Y Y Y Y 14 16 16 16 16 16 16 16 22 16 16 16 13 16 16 22 16 pass  exp1
Y Y Y Y 13 13 16 13 16 16 16 16 16 13 13 16 16 13 13 13 13 pass  exp2
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  exp4
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 22 16 16 16 22 16 pass  exp5
Y Y Y Y 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 pass  front_relax
Y Y Y Y 14 16 16 13 13 16 16 13 13 16 13 13 16 12 13 13 16 pass  global_ocean.90x40x15
Y Y Y Y 10 16 16 13 13 16 13 16 16 13 13 13 13 16 16 13 16 FAIL  global_ocean.cs32x15
Y Y Y Y  6 11 12 13 13 12 13 16 13  9  9  9  9 10  9  9 11 FAIL  global_ocean_pressure
Y Y Y Y 14 16 16 13 16 16 16 13 13 13 13 13 16 12 16 13 16 pass  global_with_exf
Y Y Y Y 14 16 16 16 16 16 16 16 16 11 13 22 13 16 16  9 16 pass  hs94.128x64x5
Y Y Y Y 13 16 16 16 16 16 16 16 16 11 16 16 16 13 16 22 13 pass  hs94.1x64x5
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 16 13 13 16 16 22 13 pass  hs94.cs-32x32x5
Y Y Y Y 10 10 16 13 13 16 16 16 22 16 13 13 13 13 13 22 13 FAIL  ideal_2D_oce
Y Y Y Y  8 16 16 16 16 16 16 16 16 13 13  8 16 16 16 16 16 FAIL  internal_wave
Y Y Y Y 14 16 16 16 16 16 16 16 16 13 13 22 13 13 13 22 16 pass  inverted_barometer
Y Y Y Y 12 16 16 16 16 16 16 16 16 16 13 12 13 13 13 13 13 FAIL  lab_sea
Y Y Y Y 11 16 16 16 16 16 16 16 13 13 13 12 13 16 13 12 13 FAIL  natl_box
Y Y Y Y 16 16 16 16 16 16 16 16 22 16 16 16 16 16 16 16 16 pass  plume_on_slope
Y Y Y Y 13 16 16 16 16 13 16 16 16 16 16 16 16 13 16 16 16 pass  solid-body.cs-32x32x1
1 C $Header: $
2 C $Name: $
3
4 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5 CC Modified from calc_gs.F, 27/8/98, Mick CC
6 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
7 #include "DIC_OPTIONS.h"
8 #include "GCHEM_OPTIONS.h"
9
10 CStartOfInterFace
11 SUBROUTINE CALC_PCO2(
12 I donewt,inewtonmax,ibrackmax,
13 I t,s,diclocal,pt,sit,ta,
14 I k1local,k2local,
15 I k1plocal,k2plocal,k3plocal,
16 I kslocal,kblocal,kwlocal,
17 I ksilocal,kflocal,
18 I fflocal,btlocal,stlocal,ftlocal,
19 U pHlocal,pCO2surfloc)
20 C /==========================================================\
21 C | SUBROUTINE CALC_PCO2 |
22 C |==========================================================|
23 C | surface ocean inorganic carbon chemistry to OCMIP2 |
24 C | regulations modified from OCMIP2 code; |
25 C | Mick Follows, MIT, Oct 1999. |
26 C \==========================================================/
27 IMPLICIT NONE
28
29 C == GLobal variables ==
30 #include "SIZE.h"
31 #include "DYNVARS.h"
32 #include "EEPARAMS.h"
33 #include "PARAMS.h"
34 #include "GRID.h"
35 #include "FFIELDS.h"
36 #include "DIC_ABIOTIC.h"
37
38 C == Routine arguments ==
39 C diclocal = total inorganic carbon (mol/m^3)
40 C where 1 T = 1 metric ton = 1000 kg
41 C ta = total alkalinity (eq/m^3)
42 C pt = inorganic phosphate (mol/^3)
43 C sit = inorganic silicate (mol/^3)
44 C t = temperature (degrees C)
45 C s = salinity (PSU)
46 INTEGER donewt
47 INTEGER inewtonmax
48 INTEGER ibrackmax
49 _RL t, s, pt, sit, ta
50 _RL pCO2surfloc, diclocal, pHlocal
51 _RL fflocal, btlocal, stlocal, ftlocal
52 _RL k1local, k2local
53 _RL k1plocal, k2plocal, k3plocal
54 _RL kslocal, kblocal, kwlocal, ksilocal, kflocal
55 CEndOfInterface
56
57 C == Local variables ==
58 C INPUT
59 C phlo= lower limit of pH range
60 C phhi= upper limit of pH range
61 C atmpres = atmospheric pressure in atmospheres (1 atm==1013.25mbar)
62 C OUTPUT
63 C co2star = CO2*water (mol/m^3)
64 C pco2surf = oceanic pCO2 (ppmv)
65 C ---------------------------------------------------------------------
66 C OCMIP NOTE: Some words about units - (JCO, 4/4/1999)
67 C - Models carry tracers in mol/m^3 (on a per volume basis)
68 C - Conversely, this routine, which was written by
69 C observationalists (C. Sabine and R. Key), passes input
70 C arguments in umol/kg (i.e., on a per mass basis)
71 C - I have changed things slightly so that input arguments are in
72 C mol/m^3,
73 C - Thus, all input concentrations (diclocal, ta, pt, and st) should be
74 C given in mol/m^3; output arguments "co2star" and "dco2star"
75 C are likewise be in mol/m^3.
76 C ---------------------------------------------------------------------
77 _RL phhi
78 _RL phlo
79 _RL tk
80 _RL tk100
81 _RL tk1002
82 _RL dlogtk
83 _RL sqrtis
84 _RL sqrts
85 _RL s15
86 _RL scl
87 _RL c
88 _RL a
89 _RL a2
90 _RL da
91 _RL b
92 _RL b2
93 _RL db
94 _RL fn
95 _RL df
96 _RL deltax
97 _RL x
98 _RL x1
99 _RL x2
100 _RL x3
101 _RL xmid
102 _RL ftest
103 _RL htotal
104 _RL htotal2
105 _RL s2
106 _RL xacc
107 _RL co2star
108 _RL co2starair
109 _RL dco2star
110 _RL dpCO2
111 _RL phguess
112 _RL atmpres
113 INTEGER inewton
114 INTEGER ibrack
115 INTEGER hstep
116 _RL fni(3)
117 _RL xlo
118 _RL xhi
119 _RL xguess
120 _RL invtk
121 _RL is
122 _RL is2
123 _RL k123p
124 _RL k12p
125 _RL k12
126 c ---------------------------------------------------------------------
127 c import donewt flag
128 c set donewt = 1 for newton-raphson iteration
129 c set donewt = 0 for bracket and bisection
130 c ---------------------------------------------------------------------
131 C Change units from the input of mol/m^3 -> mol/kg:
132 c (1 mol/m^3) x (1 m^3/1024.5 kg)
133 c where the ocean's mean surface density is 1024.5 kg/m^3
134 c Note: mol/kg are actually what the body of this routine uses
135 c for calculations. Units are reconverted back to mol/m^3 at the
136 c end of this routine.
137 c ---------------------------------------------------------------------
138 c To convert input in mol/m^3 -> mol/kg
139 pt=pt*permil
140 sit=sit*permil
141 ta=ta*permil
142 diclocal=diclocal*permil
143 c ---------------------------------------------------------------------
144 c set first guess and brackets for [H+] solvers
145 c first guess (for newton-raphson)
146 phguess = phlocal
147
148
149 c bracketing values (for bracket/bisection)
150 phhi = 10.0
151 phlo = 5.0
152 c convert to [H+]...
153 xguess = 10.0**(-phguess)
154 xlo = 10.0**(-phhi)
155 xhi = 10.0**(-phlo)
156 xmid = (xlo + xhi)*0.5
157
158
159 c----------------------------------------------------------------
160 c iteratively solve for [H+]
161 c (i) Newton-Raphson method with fixed number of iterations,
162 c use previous [H+] as first guess
163
164 c select newton-raphson, inewt=1
165 c else select bracket and bisection
166
167 cQQQQQ
168 if( donewt .eq. 1)then
169 c.........................................................
170 c NEWTON-RAPHSON METHOD
171 c.........................................................
172 x = xguess
173 cdiags
174 c WRITE(0,*)'xguess ',xguess
175 cdiags
176 do 100 inewton = 1, inewtonmax
177 c set some common combinations of parameters used in
178 c the iterative [H+] solvers
179 x2=x*x
180 x3=x2*x
181 k12 = k1local*k2local
182 k12p = k1plocal*k2plocal
183 k123p = k12p*k3plocal
184 c = 1.0 + stlocal/kslocal
185 a = x3 + k1plocal*x2 + k12p*x + k123p
186 a2=a*a
187 da = 3.0*x2 + 2.0*k1plocal*x + k12p
188 b = x2 + k1local*x + k12
189 b2=b*b
190 db = 2.0*x + k1local
191
192 c Evaluate f([H+]) and f'([H+])
193 c fn = hco3+co3+borate+oh+hpo4+2*po4+silicate+hfree
194 c +hso4+hf+h3po4-ta
195 fn = k1local*x*diclocal/b +
196 & 2.0*diclocal*k12/b +
197 & btlocal/(1.0 + x/kblocal) +
198 & kwlocal/x +
199 & pt*k12p*x/a +
200 & 2.0*pt*k123p/a +
201 & sit/(1.0 + x/ksilocal) -
202 & x/c -
203 & stlocal/(1.0 + kslocal/x/c) -
204 & ftlocal/(1.0 + kflocal/x) -
205 & pt*x3/a -
206 & ta
207
208 c df = dfn/dx
209 cdiags
210 c WRITE(0,*)'values',b2,kblocal,x2,a2,c,x
211 cdiags
212 df = ((k1local*diclocal*b) - k1local*x*diclocal*db)/b2 -
213 & 2.0*diclocal*k12*db/b2 -
214 & btlocal/kblocal/(1.0+x/kblocal)**2. -
215 & kwlocal/x2 +
216 & (pt*k12p*(a - x*da))/a2 -
217 & 2.0*pt*k123p*da/a2 -
218 & sit/ksilocal/(1.0+x/ksilocal)**2. +
219 & 1.0/c +
220 & stlocal*(1.0 + kslocal/x/c)**(-2.0)*(kslocal/c/x2) +
221 & ftlocal*(1.0 + kflocal/x)**(-2.)*kflocal/x2 -
222 & pt*x2*(3.0*a-x*da)/a2
223 c evaluate increment in [H+]
224 deltax = - fn/df
225 c update estimate of [H+]
226 x = x + deltax
227 cdiags
228 c write value of x to check convergence....
229 c write(0,*)'inewton, x, deltax ',inewton, x, deltax
230 c write(6,*)
231 cdiags
232
233 100 end do
234 c end of newton-raphson method
235 c....................................................
236 else
237 c....................................................
238 C BRACKET AND BISECTION METHOD
239 c....................................................
240 c (ii) If first step use Bracket and Bisection method
241 c with fixed, large number of iterations
242 do 200 ibrack = 1, ibrackmax
243 do hstep = 1,3
244 if(hstep .eq. 1)x = xhi
245 if(hstep .eq. 2)x = xlo
246 if(hstep .eq. 3)x = xmid
247 c set some common combinations of parameters used in
248 c the iterative [H+] solvers
249
250
251 x2=x*x
252 x3=x2*x
253 k12 = k1local*k2local
254 k12p = k1plocal*k2plocal
255 k123p = k12p*k3plocal
256 c = 1.0 + stlocal/kslocal
257 a = x3 + k1plocal*x2 + k12p*x + k123p
258 a2=a*a
259 da = 3.0*x2 + 2.0*k1plocal*x + k12p
260 b = x2 + k1local*x + k12
261 b2=b*b
262 db = 2.0*x + k1local
263 c evaluate f([H+]) for bracketing and mid-value cases
264 fn = k1local*x*diclocal/b +
265 & 2.0*diclocal*k12/b +
266 & btlocal/(1.0 + x/kblocal) +
267 & kwlocal/x +
268 & pt*k12p*x/a +
269 & 2.0*pt*k123p/a +
270 & sit/(1.0 + x/ksilocal) -
271 & x/c -
272 & stlocal/(1.0 + kslocal/x/c) -
273 & ftlocal/(1.0 + kflocal/x) -
274 & pt*x3/a -
275 & ta
276 fni(hstep) = fn
277 end do
278 c now bracket solution within two of three
279 ftest = fni(1)/fni(3)
280 if(ftest .gt. 0.0)then
281 xhi = xmid
282 else
283 xlo = xmid
284 end if
285 xmid = (xlo + xhi)*0.5
286
287 cdiags
288 c write value of x to check convergence....
289 c WRITE(0,*)'bracket-bisection iteration ',ibrack, xmid
290 cdiags
291 200 end do
292 c last iteration gives value
293 x = xmid
294 c end of bracket and bisection method
295 c....................................
296 end if
297 c iterative [H+] solver finished
298 c----------------------------------------------------------------
299
300 c now determine pCO2 etc...
301 c htotal = [H+], hydrogen ion conc
302 htotal = x
303 C Calculate [CO2*] as defined in DOE Methods Handbook 1994 Ver.2,
304 C ORNL/CDIAC-74, dickson and Goyet, eds. (Ch 2 p 10, Eq A.49)
305 htotal2=htotal*htotal
306 co2star=diclocal*htotal2/(htotal2 + k1local*htotal
307 & + k1local*k2local)
308 phlocal=-log10(htotal)
309
310 c ---------------------------------------------------------------
311 c Add two output arguments for storing pCO2surf
312 c Should we be using K0 or ff for the solubility here?
313 c ---------------------------------------------------------------
314 pCO2surfloc = co2star / fflocal
315
316 C ----------------------------------------------------------------
317 C Reconvert units back to original values for input arguments
318 C no longer necessary????
319 C ----------------------------------------------------------------
320 c Reconvert from mol/kg -> mol/m^3
321 pt=pt/permil
322 sit=sit/permil
323 ta=ta/permil
324 diclocal=diclocal/permil
325
326 return
327 end
328
329 c=================================================================
330 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
331 CC New efficient pCO2 solver, Mick Follows CC
332 CC Taka Ito CC
333 CC Stephanie Dutkiewicz CC
334 CC 20 April 2003 CC
335 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
336 #include "DIC_OPTIONS.h"
337 CStartOfInterFace
338 SUBROUTINE CALC_PCO2_APPROX(
339 I t,s,diclocal,pt,sit,ta,
340 I k1local,k2local,
341 I k1plocal,k2plocal,k3plocal,
342 I kslocal,kblocal,kwlocal,
343 I ksilocal,kflocal,
344 I fflocal,btlocal,stlocal,ftlocal,
345 U pHlocal,pCO2surfloc)
346 C /==========================================================\
347 C | SUBROUTINE CALC_PCO2_APPROX |
348 C |==========================================================|
349 C | surface ocean inorganic carbon chemistry to OCMIP2 |
350 C | regulations modified from OCMIP2 code; |
351 C | Mick Follows, MIT, Oct 1999. |
352 C \==========================================================/
353 IMPLICIT NONE
354
355 C == GLobal variables ==
356 #include "SIZE.h"
357 #include "DYNVARS.h"
358 #include "EEPARAMS.h"
359 #include "PARAMS.h"
360 #include "GRID.h"
361 #include "FFIELDS.h"
362 #include "DIC_ABIOTIC.h"
363
364 C == Routine arguments ==
365 C diclocal = total inorganic carbon (mol/m^3)
366 C where 1 T = 1 metric ton = 1000 kg
367 C ta = total alkalinity (eq/m^3)
368 C pt = inorganic phosphate (mol/^3)
369 C sit = inorganic silicate (mol/^3)
370 C t = temperature (degrees C)
371 C s = salinity (PSU)
372 _RL t, s, pt, sit, ta
373 _RL pCO2surfloc, diclocal, pHlocal
374 _RL fflocal, btlocal, stlocal, ftlocal
375 _RL k1local, k2local
376 _RL k1plocal, k2plocal, k3plocal
377 _RL kslocal, kblocal, kwlocal, ksilocal, kflocal
378 CEndOfInterface
379
380 C == Local variables ==
381 _RL phguess
382 _RL cag
383 _RL bohg
384 _RL hguess
385 _RL stuff
386 _RL gamm
387 _RL hnew
388 _RL co2s
389 _RL h3po4g, h2po4g, hpo4g, po4g
390 _RL siooh3g
391
392
393 c ---------------------------------------------------------------------
394 C Change units from the input of mol/m^3 -> mol/kg:
395 c (1 mol/m^3) x (1 m^3/1024.5 kg)
396 c where the ocean's mean surface density is 1024.5 kg/m^3
397 c Note: mol/kg are actually what the body of this routine uses
398 c for calculations. Units are reconverted back to mol/m^3 at the
399 c end of this routine.
400 c To convert input in mol/m^3 -> mol/kg
401 pt=pt*permil
402 sit=sit*permil
403 ta=ta*permil
404 diclocal=diclocal*permil
405 c ---------------------------------------------------------------------
406 c set first guess and brackets for [H+] solvers
407 c first guess (for newton-raphson)
408 phguess = phlocal
409 cmick - new approx method
410 cmick - make estimate of htotal (hydrogen ion conc) using
411 cmick appromate estimate of CA, carbonate alkalinity
412 hguess = 10.0**(-phguess)
413 cmick - first estimate borate contribution using guess for [H+]
414 bohg = btlocal*kblocal/(hguess+kblocal)
415
416 cmick - first estimate of contribution from phosphate
417 cmick based on Dickson and Goyet
418 stuff = hguess*hguess*hguess
419 & + (k1plocal*hguess*hguess)
420 & + (k1plocal*k2plocal*hguess)
421 & + (k1plocal*k2plocal*k3plocal)
422 h3po4g = (pt*hguess*hguess*hguess) / stuff
423 h2po4g = (pt*k1plocal*hguess*hguess) / stuff
424 hpo4g = (pt*k1plocal*k2plocal*hguess) / stuff
425 po4g = (pt*k1plocal*k2plocal*k3plocal) / stuff
426
427 cmick - estimate contribution from silicate
428 cmick based on Dickson and Goyet
429 siooh3g = sit*ksilocal / (ksilocal + hguess)
430
431 cmick - now estimate carbonate alkalinity
432 cag = ta - bohg - (kwlocal/hguess) + hguess
433 & - hpo4g - 2.0*po4g + h3po4g
434 & - siooh3g
435
436 coldcmick - now estimate carbonate alkalinity
437 cold cag = ta - bohg - (kwlocal/hguess) + hguess
438 coldcmick - NB could add further corrections for other contributions
439 coldcmick e.g. Si, PO4, NO3 ...
440
441 cmick - now evaluate better guess of hydrogen ion conc
442 cmick htotal = [H+], hydrogen ion conc
443 gamm = diclocal/cag
444 stuff = (1.0-gamm)*(1.0-gamm)*k1local*k1local
445 & - 4.0*k1local*k2local*(1.0-2.0*gamm)
446 hnew = 0.5*( (gamm-1.0)*k1local + sqrt(stuff) )
447 cmick - now determine [CO2*]
448 co2s = diclocal/
449 & (1.0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))
450 cmick - return update pH to main routine
451 phlocal = -log10(hnew)
452
453 c ---------------------------------------------------------------
454 c surface pCO2 (following Dickson and Goyet, DOE...)
455 pCO2surfloc = co2s/fflocal
456
457 C ----------------------------------------------------------------
458 c Reconvert from mol/kg -> mol/m^3
459 pt=pt/permil
460 sit=sit/permil
461 ta=ta/permil
462 diclocal=diclocal/permil
463 return
464 end
465
466 c=================================================================
467 c *******************************************************************
468 c=================================================================
469 CStartOfInterFace
470 SUBROUTINE CARBON_COEFFS(
471 I ttemp,stemp,
472 I bi,bj,iMin,iMax,jMin,jMax)
473 C
474 C /==========================================================\
475 C | SUBROUTINE CARBON_COEFFS |
476 C | determine coefficients for surface carbon chemistry |
477 C | adapted from OCMIP2: SUBROUTINE CO2CALC |
478 C | mick follows, oct 1999 |
479 c | minor changes to tidy, swd aug 2002 |
480 C \==========================================================/
481 C INPUT
482 C diclocal = total inorganic carbon (mol/m^3)
483 C where 1 T = 1 metric ton = 1000 kg
484 C ta = total alkalinity (eq/m^3)
485 C pt = inorganic phosphate (mol/^3)
486 C sit = inorganic silicate (mol/^3)
487 C t = temperature (degrees C)
488 C s = salinity (PSU)
489 C OUTPUT
490 C IMPORTANT: Some words about units - (JCO, 4/4/1999)
491 c - Models carry tracers in mol/m^3 (on a per volume basis)
492 c - Conversely, this routine, which was written by observationalists
493 c (C. Sabine and R. Key), passes input arguments in umol/kg
494 c (i.e., on a per mass basis)
495 c - I have changed things slightly so that input arguments are in mol/m^3,
496 c - Thus, all input concentrations (diclocal, ta, pt, and st) should be
497 c given in mol/m^3; output arguments "co2star" and "dco2star"
498 c are likewise be in mol/m^3.
499 C--------------------------------------------------------------------------
500 IMPLICIT NONE
501 C == GLobal variables ==
502 #include "SIZE.h"
503 #include "DYNVARS.h"
504 #include "EEPARAMS.h"
505 #include "PARAMS.h"
506 #include "GRID.h"
507 #include "FFIELDS.h"
508 #include "DIC_ABIOTIC.h"
509 C == Routine arguments ==
510 C ttemp and stemp are local theta and salt arrays
511 C dont really need to pass T and S in, could use theta, salt in
512 C common block in DYNVARS.h, but this way keeps subroutine more
513 C general
514 _RL ttemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
515 _RL stemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
516 INTEGER bi,bj,iMin,iMax,jMin,jMax
517 CEndOfInterface
518
519 C LOCAL VARIABLES
520 _RL t
521 _RL s
522 _RL ta
523 _RL pt
524 _RL sit
525 _RL tk
526 _RL tk100
527 _RL tk1002
528 _RL dlogtk
529 _RL sqrtis
530 _RL sqrts
531 _RL s15
532 _RL scl
533 _RL x1
534 _RL x2
535 _RL s2
536 _RL xacc
537 _RL invtk
538 _RL is
539 _RL is2
540 INTEGER i
541 INTEGER j
542
543 C.....................................................................
544 C OCMIP note:
545 C Calculate all constants needed to convert between various measured
546 C carbon species. References for each equation are noted in the code.
547 C Once calculated, the constants are
548 C stored and passed in the common block "const". The original version
549 C of this code was based on the code by dickson in Version 2 of
550 C "Handbook of Methods C for the Analysis of the Various Parameters of
551 C the Carbon Dioxide System in Seawater", DOE, 1994 (SOP No. 3, p25-26).
552 C....................................................................
553
554 do i=imin,imax
555 do j=jmin,jmax
556 if (hFacC(i,j,1,bi,bj).gt.0.d0) then
557 t = ttemp(i,j,1,bi,bj)
558 s = stemp(i,j,1,bi,bj)
559 C terms used more than once
560 tk = 273.15 + t
561 tk100 = tk/100.0
562 tk1002=tk100*tk100
563 invtk=1.0/tk
564 dlogtk=log(tk)
565 is=19.924*s/(1000.-1.005*s)
566 is2=is*is
567 sqrtis=sqrt(is)
568 s2=s*s
569 sqrts=sqrt(s)
570 s15=s**1.5
571 scl=s/1.80655
572 C------------------------------------------------------------------------
573 C f = k0(1-pH2O)*correction term for non-ideality
574 C Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
575 ff(i,j,bi,bj) = exp(-162.8301 + 218.2968/tk100 +
576 & 90.9241*log(tk100) - 1.47696*tk1002 +
577 & s * (.025695 - .025225*tk100 +
578 & 0.0049867*tk1002))
579 C------------------------------------------------------------------------
580 C K0 from Weiss 1974
581 ak0(i,j,bi,bj) = exp(93.4517/tk100 - 60.2409 +
582 & 23.3585 * log(tk100) +
583 & s * (0.023517 - 0.023656*tk100 +
584 & 0.0047036*tk1002))
585 C------------------------------------------------------------------------
586 C k1 = [H][HCO3]/[H2CO3]
587 C k2 = [H][CO3]/[HCO3]
588 C Millero p.664 (1995) using Mehrbach et al. data on seawater scale
589 ak1(i,j,bi,bj)=10**(-1*(3670.7*invtk -
590 & 62.008 + 9.7944*dlogtk -
591 & 0.0118 * s + 0.000116*s2))
592 ak2(i,j,bi,bj)=10**(-1*(1394.7*invtk + 4.777 -
593 & 0.0184*s + 0.000118*s2))
594 C------------------------------------------------------------------------
595 C kb = [H][BO2]/[HBO2]
596 C Millero p.669 (1995) using data from dickson (1990)
597 akb(i,j,bi,bj)=exp((-8966.90 - 2890.53*sqrts - 77.942*s +
598 & 1.728*s15 - 0.0996*s2)*invtk +
599 & (148.0248 + 137.1942*sqrts + 1.62142*s) +
600 & (-24.4344 - 25.085*sqrts - 0.2474*s) *
601 & dlogtk + 0.053105*sqrts*tk)
602 C------------------------------------------------------------------------
603 C k1p = [H][H2PO4]/[H3PO4]
604 C DOE(1994) eq 7.2.20 with footnote using data from Millero (1974)
605 ak1p(i,j,bi,bj) = exp(-4576.752*invtk + 115.525 -
606 & 18.453*dlogtk +
607 & (-106.736*invtk + 0.69171)*sqrts +
608 & (-0.65643*invtk - 0.01844)*s)
609 C------------------------------------------------------------------------
610 C k2p = [H][HPO4]/[H2PO4]
611 C DOE(1994) eq 7.2.23 with footnote using data from Millero (1974))
612 ak2p(i,j,bi,bj) = exp(-8814.715*invtk + 172.0883 -
613 & 27.927*dlogtk +
614 & (-160.340*invtk + 1.3566) * sqrts +
615 & (0.37335*invtk - 0.05778) * s)
616 C------------------------------------------------------------------------
617 C k3p = [H][PO4]/[HPO4]
618 C DOE(1994) eq 7.2.26 with footnote using data from Millero (1974)
619 ak3p(i,j,bi,bj) = exp(-3070.75*invtk - 18.141 +
620 & (17.27039*invtk + 2.81197) *
621 & sqrts + (-44.99486*invtk - 0.09984) * s)
622 C------------------------------------------------------------------------
623 C ksi = [H][SiO(OH)3]/[Si(OH)4]
624 C Millero p.671 (1995) using data from Yao and Millero (1995)
625 aksi(i,j,bi,bj) = exp(-8904.2*invtk + 117.385 -
626 & 19.334*dlogtk +
627 & (-458.79*invtk + 3.5913) * sqrtis +
628 & (188.74*invtk - 1.5998) * is +
629 & (-12.1652*invtk + 0.07871) * is2 +
630 & log(1.0-0.001005*s))
631 C------------------------------------------------------------------------
632 C kw = [H][OH]
633 C Millero p.670 (1995) using composite data
634 akw(i,j,bi,bj) = exp(-13847.26*invtk + 148.9652 -
635 & 23.6521*dlogtk +
636 & (118.67*invtk - 5.977 + 1.0495 * dlogtk) *
637 & sqrts - 0.01615 * s)
638 C------------------------------------------------------------------------
639 C ks = [H][SO4]/[HSO4]
640 C dickson (1990, J. chem. Thermodynamics 22, 113)
641 aks(i,j,bi,bj)=exp(-4276.1*invtk + 141.328 -
642 & 23.093*dlogtk +
643 & (-13856*invtk + 324.57 - 47.986*dlogtk)*sqrtis +
644 & (35474*invtk - 771.54 + 114.723*dlogtk)*is -
645 & 2698*invtk*is**1.5 + 1776*invtk*is2 +
646 & log(1.0 - 0.001005*s))
647 C------------------------------------------------------------------------
648 C kf = [H][F]/[HF]
649 C dickson and Riley (1979) -- change pH scale to total
650 akf(i,j,bi,bj)=exp(1590.2*invtk - 12.641 + 1.525*sqrtis +
651 & log(1.0 - 0.001005*s) +
652 & log(1.0 + (0.1400/96.062)*(scl)/aks(i,j,bi,bj)))
653 C------------------------------------------------------------------------
654 C Calculate concentrations for borate, sulfate, and fluoride
655 C Uppstrom (1974)
656 bt(i,j,bi,bj) = 0.000232 * scl/10.811
657 C Morris & Riley (1966)
658 st(i,j,bi,bj) = 0.14 * scl/96.062
659 C Riley (1965)
660 ft(i,j,bi,bj) = 0.000067 * scl/18.9984
661 C------------------------------------------------------------------------
662 else
663 ff(i,j,bi,bj)=0.d0
664 ak0(i,j,bi,bj)= 0.d0
665 ak1(i,j,bi,bj)= 0.d0
666 ak2(i,j,bi,bj)= 0.d0
667 akb(i,j,bi,bj)= 0.d0
668 ak1p(i,j,bi,bj) = 0.d0
669 ak2p(i,j,bi,bj) = 0.d0
670 ak3p(i,j,bi,bj) = 0.d0
671 aksi(i,j,bi,bj) = 0.d0
672 akw(i,j,bi,bj) = 0.d0
673 aks(i,j,bi,bj)= 0.d0
674 akf(i,j,bi,bj)= 0.d0
675 bt(i,j,bi,bj) = 0.d0
676 st(i,j,bi,bj) = 0.d0
677 ft(i,j,bi,bj) = 0.d0
678 endif
679 end do
680 end do
681
682 return
683 end
684
  ViewVC Help
Powered by ViewVC 1.1.22