source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmpco2bal.F90 @ 8456

Last change on this file since 8456 was 8456, checked in by dford, 3 years ago

Add pCO2/fCO2 assimilation.

File size: 12.5 KB
Line 
1MODULE asmpco2bal
2   !Contains subroutine to convert a global pco2 increment to independent DIC and alkalinity increments
3
4   USE par_kind, ONLY: wp
5   USE par_oce, ONLY: jpi, jpj, jpk
6   USE dom_oce, ONLY: tmask
7#if defined key_hadocc
8   USE had_bgc_const, ONLY: RHO_WATER_SI
9#else
10   USE phycst, ONLY: RHO_WATER_SI => rau0
11#endif
12   
13   IMPLICIT NONE
14   PRIVATE
15   
16   PUBLIC asm_pco2_bal
17   CONTAINS
18
19   SUBROUTINE  asm_pco2_bal(pco2inc, tco2_tracer_bkg, alk_tracer_bkg, &
20      &                     tem_tracer_bkg, sal_tracer_bkg,           &
21      &                     tco2_tracer_inc, alk_tracer_inc)
22      !Converts a global pco2 increment to independent DIC and alkalinity increments
23      !Calculates the solution to the equation [DDIC, DAlk]=(DpCO2 / |Grad (pCO2)|)* Grad(pCO2); where Grad(pco2)=[dpCO2/dDIC, dpCO2/dalk]
24      !Gradient code pco2_grad is a modification of assimcd1.mf77 by Ian Totterdell.
25
26      !James While 06/09
27
28      IMPLICIT NONE
29
30      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) ::  pco2inc            !pCO2 increment
31     
32      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
33                                                  alk_tracer_bkg      !In alkalinity background
34     
35      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
36                                                  sal_tracer_bkg      !In salinity background
37     
38      REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: tco2_tracer_inc,      &    !Out: tCO2 (DIC) increment
39                                                   alk_tracer_inc      !Out alkalinity
40                     
41      !local variables
42      REAL(wp), DIMENSION(jpi,jpj) ::  Dpco2_tco2, Dpco2_alk
43      REAL(wp) :: wgt_fct
44      INTEGER :: ji,jj
45     
46     
47      !calculate the gradients with respect to DIC and Alkalinity
48      CALL pco2_grad(tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, Dpco2_tco2, Dpco2_alk)
49
50     
51      !loop over the ocean grid calculating the increments to tCO2 (DIC) and alkalinity
52      DO jj=1, jpj
53        DO ji=1, jpi
54
55          wgt_fct=pco2inc(ji,jj)/( (Dpco2_tco2(ji,jj)*Dpco2_tco2(ji,jj)) + (Dpco2_alk(ji,jj)*Dpco2_alk(ji,jj)) )
56          tco2_tracer_inc(ji,jj)=wgt_fct * Dpco2_tco2(ji,jj) * tmask(ji,jj,1)
57          alk_tracer_inc(ji,jj)=wgt_fct * Dpco2_alk(ji,jj) * tmask(ji,jj,1)
58 
59        ENDDO
60      ENDDO
61                         
62                                   
63
64   END SUBROUTINE  asm_pco2_bal
65   
66 
67   SUBROUTINE pco2_grad( tco2_tracer_bkg, alk_tracer_bkg, tem_tracer_bkg, sal_tracer_bkg, DPDC, DPDA )
68!
69!-------M---------M---------M---------M---------M---------M---------M--
70!
71!  Ian J Totterdell
72!  2007 August 08
73!
74!  Code for calculating derivatives of pCO2 wrt DIC and Alk
75
76!  J While
77!  June 2009
78!  Modiefied for use in NEMO
79
80!-------M---------M---------M---------M---------M---------M---------M--
81      IMPLICIT NONE
82!
83!
84      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tco2_tracer_bkg,      &    !In: tCO2 (DIC) background
85                                                  alk_tracer_bkg      !In alkalinity background
86     
87      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: tem_tracer_bkg,      &    !In: temperature background
88                                                  sal_tracer_bkg      !In salinity background
89!
90      REAL(wp), INTENT(OUT), DIMENSION(jpi,jpj) :: &
91        DPDC,   &
92        DPDA
93!
94!
95!  Local variables.
96!
97      INTEGER :: &
98        mxiter,  &
99        ji,jj,   &
100        iter
101!
102      REAL(wp) :: &
103        a1_co2, a2_co2, a3_co2,               &
104        b1_co2, b2_co2, b3_co2,               &
105        a1_k1, a2_k1, a3_k1,                  &
106        b1_k1, b2_k1, b3_k1, b4_k1, b5_k1,    &
107        a1_k2, a2_k2, a3_k2,                  &
108        b1_k2, b2_k2, b3_k2, b4_k2, b5_k2,    &
109        a1_kb, a2_kb, a3_kb, a4_kb, a5_kb,    &
110        b1_kb, b2_kb, b3_kb,                  &
111        c1_kb, c2_kb, c3_kb, d_kb,            &
112        a1_kw, a2_kw, a3_kw,                  &
113        b1_kw, b2_kw, b3_kw, c_kw
114!
115      REAL(wp) :: &
116        tk, logtk, ttk,             &
117        sal, sqsal,                 &
118        kone, ktwo, kboric, kwater, &
119        alkprime, qa,qb,qc,         &
120        fxa, xtrma, xtrmb,          &
121        dxda, dxdc,                 &
122        dgda, dgdc
123!
124      REAL(wp), DIMENSION(jpi,jpj) :: &
125      co2_sol,     &
126      rek,         &
127      rg1,         &
128      rg2,         &
129      rg3,         &
130      boron,       &
131      tco2,        &
132      alk_t,       &
133      xs,          &
134      x1,          &
135      x2,          &
136      x2p,         &
137      x3,          &
138      xx,          &
139      xf
140!
141!  Set parameter and constant values.
142!
143      PARAMETER ( mxiter = 10 )
144!
145!  CO2 solubility, constants for the calculation of.
146!   From Weiss (1974)
147!
148      DATA &
149      a1_co2 / -60.2409 /,  a2_co2 / 93.4517 /,   a3_co2 / 23.3585 /, &
150      b1_co2 / 0.023517 /,  b2_co2 / -0.023656 /, b3_co2 / 0.0047036 /
151!
152!  kone, constants for the calculation of.
153!   From Roy et al (1993)
154!
155      DATA a1_k1 /-2307.1266/,  a2_k1 / 2.83655 /,   a3_k1 /-1.5529413/ &
156        , b1_k1 / -4.0484 /,   b2_k1 /-0.20760841/, b3_k1 /0.08468345/ &
157        , b4_k1 /-0.00654208/, b5_k1 / -0.001005 /
158!
159!  ktwo, constants for the calculation of.
160!   From Roy et al (1993)
161!
162      DATA a1_k2 /-3351.6106/,  a2_k2 / -9.226508 /, a3_k2 /-0.2005743/ &
163        , b1_k2 / -23.9722 /,  b2_k2 /-0.106901773/,b3_k2 /0.1130822/ &
164        , b4_k2 /-0.00846934/, b5_k2 / -0.001005 /
165!
166!  kboric, constants for the calculation of.
167!   Using data from Dickson (1990);
168!
169      DATA a1_kb / -8966.90 /,  a2_kb / -2890.53 /,  a3_kb / -77.942 / &
170        , a4_kb / 1.728 /,     a5_kb /-0.0996/ &
171        , b1_kb / 148.0248 /,  b2_kb / 137.1942 /,  b3_kb / 1.62142 / &
172        , c1_kb / -24.4344 /,  c2_kb / -25.085 /,   c3_kb / -0.2474 / &
173        , d_kb / 0.053105 /
174!
175!  kwater, constants for the calculation of.
176!   From Millero (1995, p670), using composite data.
177!
178      DATA a1_kw / -13847.26 /, a2_kw / 148.9652 /,  a3_kw / -23.6521 / &
179        , b1_kw / 118.67 /,    b2_kw / -5.977 /,    b3_kw / 1.0495 / &
180        , c_kw / -0.01615 / 
181!
182!  Executable code begins here.
183!
184!  Calculate the carbon system constants.
185!
186
187      DO ji=1,jpi
188        DO jj=1,jpj
189!
190          tk = 273.15 + MAX( -2.0, MIN( 40.0, tem_tracer_bkg(ji,jj) ) )
191          logtk = LOG(tk)
192          ttk = MAX( 2.71, 0.01 * ( 273.15 + tem_tracer_bkg(ji,jj) ) )
193          sal = MAX( 0.0, MIN( 45.0, sal_tracer_bkg(ji,jj) ) )
194          sqsal = SQRT(sal)
195!
196!  Calculate the CO2 solubility from Weiss (1974)
197!   Units are mol/(kg.soln)/atm .
198!   Check value (Zeebe and Wolf-G; tc=25,sal=35): ln(co2_sol) = -3.5617
199!
200          co2_sol(ji,jj) = EXP( &
201           a1_co2 + a2_co2 / ttk + a3_co2 * LOG(ttk) + &
202           sal * ( b1_co2 + ttk * ( b2_co2 + b3_co2 * ttk ) ) )
203!
204!  Calculate the 1st constant of carbonic acid from Roy et al. (1993)
205!   Units are mol/(kg.soln) on the total hydrogen ion scale.
206!   K1 = (HCO3-).(H+)/(CO2(g))
207!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K1) = -13.4847
208!
209          kone &
210           = 1.0e-7 *( 1.0 + b5_k1*sal ) * EXP( a1_k1/tk + a2_k1 + a3_k1*logtk &
211           + ( b1_k1/tk + b2_k1 )*sqsal + b3_k1*sal + b4_k1*sal*sqsal )
212!
213!  Calculate the 2nd constant of carbonic acid from Roy et al. (1993)
214!   Units are mol/(kg.soln) on the total hydrogen ion scale.
215!   K2 = (CO3--).(H+)/(HCO3-)
216!   Check value (Dickson and Goyet; tc=25,sal=35): ln(K2) = -20.5504
217!
218          ktwo &
219           = 1.0e-7 *( 1.0 + b5_k2*sal ) * EXP( a1_k2/tk + a2_k2 + a3_k2*logtk &
220           + ( b1_k2/tk + b2_k2 )*sqsal + b3_k2*sal + b4_k2*sal*sqsal )
221!
222!  Calculate the equilibrium constant of boric acid from Dickson (1990).
223!   Units are mol/(kg.soln) on the total hydrogen ion scale.
224!   Kb = (B(OH)4-).(H+)/(B(OH)3) (with (Bt) = (B(OH)3) + (B(OH)4-) ).
225!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kb) = -19.7964
226!
227          kboric &
228           = 1.0e-7 * EXP( ( a1_kb + a2_kb*sqsal + a3_kb*sal + a4_kb*sal*sqsal &
229           + a5_kb*sal*sal )/tk + ( b1_kb + b2_kb*sqsal + b3_kb*sal ) &
230           + ( c1_kb + c2_kb*sqsal + c3_kb*sal )*logtk + d_kb*sqsal*tk )
231!
232!  Calculate the equilibrium constant of water from Millero (1994).
233!   Units are ( mol/(kg.soln) )^2 on the total hydrogen ion scale.
234!   Kw = (OH-).(H+)
235!   Check value (Dickson and Goyet; tc=25,sal=35): ln(Kw) = -30.434
236!
237          kwater &
238           = 1.0e-14  *EXP( ( a1_kw/tk + a2_kw + a3_kw*logtk ) &
239           + ( b1_kw/tk + b2_kw + b3_kw*logtk )*sqsal + c_kw*sal )
240!
241!  Calculate the combinations (as in PPCO2).
242!
243          rek(ji,jj) = SQRT( kone / ktwo )
244          rg3(ji,jj) = SQRT( kone * ktwo )
245          rg1(ji,jj) = rg3(ji,jj) / kboric
246          rg2(ji,jj) = kwater / rg3(ji,jj)
247!
248!  Set Boron concentration according to the surface salinity
249!   (Boron accounts for about 10% of the total Alkalinity).
250!   Use relationship used by Peng et al, Tellus 39b 1987,
251!   originally from Culkin 1965:
252!
253!         boron(ji,jj) = 4.16e-4 * sal / 35.0   ! Correct value but not in UM
254          boron(ji,jj) = 4.106e-4 * sal / 35.0
255!
256!  tco2 and alk_t are got from the modelled tracers at the surface
257!   (converted to units of mol/kg).
258!
259          tco2(ji,jj)  = MAX( 0.0, tco2_tracer_bkg(ji,jj) * 0.001 / RHO_WATER_SI )
260          alk_t(ji,jj) = MAX( 0.0, alk_tracer_bkg(ji,jj)  * 0.001 / RHO_WATER_SI )
261!
262!  Set initial values for iterative method.
263!
264          xs(ji,jj) = 1.0
265          x2p(ji,jj) = 0.0
266          x3(ji,jj) = 0.0
267          xx(ji,jj) = 0.0
268          xf(ji,jj) = xs(ji,jj)
269!
270        ENDDO
271      ENDDO
272
273!
274!  Use the secant method of similar triangles to calculate the H+ ion
275!   concentration, [H] (though strictly xf refers to  rg3/[H] ).
276!   To aid vectorisation this has been written to iterate a set
277!   number of times (mxiter), large enough to bring every example to
278!   within the 1e-5 fractional change criterion.
279!
280!  The method follows that of Bacastow (1990) in separating the
281!   alkalinity into two parts: the major part is the carbonate alk
282!   which uses the [H] that we are trying to find, while the minor
283!   part (xw) is sum of the borate and water dissociation terms which
284!   use the previous approximation to [H].
285!   This leaves an easy-to-solve quadratic.
286!
287      DO iter=1,mxiter
288!
289
290        DO ji=1,jpi
291          DO jj=1,jpj
292!
293            IF((abs(xf(ji,jj)-xx(ji,jj))/xf(ji,jj))<1.e-5) CYCLE   !prevent iteration collapsing to 0.0/0.0
294
295
296            x1(ji,jj) = x2p(ji,jj)
297            x2(ji,jj) = x3(ji,jj)
298            x2p(ji,jj) = xx(ji,jj)
299            x3(ji,jj) = xf(ji,jj)
300
301            IF (iter<=2) THEN
302              xx(ji,jj) = xf(ji,jj)
303            ELSE
304              xx(ji,jj) = x2p(ji,jj) + ( x3(ji,jj) - x2p(ji,jj) )*( x1(ji,jj) - x2p(ji,jj) ) &
305              /( x1(ji,jj) - x2(ji,jj) - x2p(ji,jj) + x3(ji,jj) )
306            ENDIF
307
308
309!
310!  All the above has been setting up the similar triangles for the
311!   secant method and using it to find an improved interpolated
312!   estimate; it has been completely independent of the equations
313!   we are trying to solve.
314!
315!  The next bit, however, IS dependent on the equations.
316!   alkprime is the carbonate alkalinity, which is calculated
317!   from the total alkalinity by subtracting the other, smaller, terms,
318!   which use the interpolated previous approximation to xf (xx).
319!
320
321            alkprime = alk_t(ji,jj) - rg2(ji,jj)*xx(ji,jj) + rg3(ji,jj)/xx(ji,jj) &
322            - boron(ji,jj)/( 1.0 + rg1(ji,jj)/xx(ji,jj) )       
323       
324!
325!  Set up the quadratic terms.
326!
327
328            qa = 2.0*tco2(ji,jj) - alkprime
329            qb = -rek(ji,jj)*( alkprime - tco2(ji,jj) )
330            qc = -alkprime
331
332            xf(ji,jj) = ( -qb + SQRT( qb*qb - 4.0*qa*qc ) )*0.5/qa
333
334
335          ENDDO
336        ENDDO
337!     
338      ENDDO
339
340!
341      DO ji=1,jpi
342        DO jj=1,jpj
343!
344          xtrma = 1.0/( 1.0 + rek(ji,jj)*xf(ji,jj) + xf(ji,jj)*xf(ji,jj) )
345          xtrmb = rg2(ji,jj) + rg3(ji,jj)/(xf(ji,jj)*xf(ji,jj)) &
346          + rg1(ji,jj)*boron(ji,jj)/(( rg1(ji,jj) + xf(ji,jj) )*( rg1(ji,jj) + xf(ji,jj) ))
347!
348          dxda = 1.0/( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) &
349          + rek(ji,jj)*xf(ji,jj)*xf(ji,jj) )*xtrma*xtrma + xtrmb )
350          dxdc = -dxda*xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma
351!
352          dgdc = xtrma - dxdc*tco2(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) )*xtrma*xtrma
353          dgda = ( 1.0 - dxda*( tco2(ji,jj)*( rek(ji,jj) + 4.0*xf(ji,jj) )*xtrma &
354          + xtrmb ) )/(xf(ji,jj)*( rek(ji,jj) + 2.0*xf(ji,jj) ))
355!
356          fxa = 1.0e-3*RHO_WATER_SI*co2_sol(ji,jj)
357          DPDC(ji,jj) = dgdc/fxa
358          DPDA(ji,jj) = dgda/fxa
359!
360        ENDDO
361      ENDDO
362         
363      RETURN
364      END SUBROUTINE pco2_grad
365
366END MODULE asmpco2bal
Note: See TracBrowser for help on using the repository browser.