1 | MODULE trcco2_medusa |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE trcco2_medusa *** |
---|
4 | !! TOP : MEDUSA |
---|
5 | !!====================================================================== |
---|
6 | !! History : |
---|
7 | !! - ! 2010-12 (A. Yool) added for ROAM project |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | #if defined key_medusa && defined key_roam |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! MEDUSA carbonate chemistry |
---|
12 | !!---------------------------------------------------------------------- |
---|
13 | !! trc_co2_medusa : |
---|
14 | !!---------------------------------------------------------------------- |
---|
15 | USE oce_trc |
---|
16 | USE trc |
---|
17 | USE sms_medusa |
---|
18 | USE lbclnk |
---|
19 | USE prtctl_trc ! Print control for debugging |
---|
20 | USE in_out_manager ! I/O manager |
---|
21 | |
---|
22 | IMPLICIT NONE |
---|
23 | PRIVATE |
---|
24 | |
---|
25 | PUBLIC trc_co2_medusa ! called in trc_bio_medusa |
---|
26 | |
---|
27 | !!* Substitution |
---|
28 | # include "domzgr_substitute.h90" |
---|
29 | !!---------------------------------------------------------------------- |
---|
30 | !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) |
---|
31 | !! $Id$ |
---|
32 | !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | |
---|
35 | ! The following is a map of the subroutines contained within this module |
---|
36 | ! - trc_co2_medusa |
---|
37 | ! - CALLS CO2_dynamics |
---|
38 | ! - CALLS CO2DYN |
---|
39 | ! - CALLS POLYCO |
---|
40 | ! - CALLS CO2SET |
---|
41 | ! - CALLS CO2CLC |
---|
42 | ! - CALLS CaCO3_Saturation |
---|
43 | ! - CALLS Air_sea_exchange |
---|
44 | |
---|
45 | CONTAINS |
---|
46 | |
---|
47 | !======================================================================= |
---|
48 | ! |
---|
49 | SUBROUTINE trc_co2_medusa( Temp, Sal, DIC, ALK, Depth, xkw, pCO2a, & |
---|
50 | pH, pCO2w, h2co3, hco3, co3, om_cal, om_arg, co2flux, TDIC, TALK, & |
---|
51 | dcf, henry, iters ) |
---|
52 | ! |
---|
53 | !======================================================================= |
---|
54 | ! |
---|
55 | ! This file contains a set of FORTRAN subroutines that calculate the carbonate system |
---|
56 | ! at any given point in marine space time, given values for |
---|
57 | ! temperature, salinity, DIC, depth (pressure). |
---|
58 | ! This is essentially an implimentation of the Haltafall speciation code |
---|
59 | ! (Ingri et al 1967, Talanta 14, 1261 - if it ain't broke don't fix it) |
---|
60 | ! Another routine calulates the air sea exchange of CO2 given wind speed and atmospheric pCO2. |
---|
61 | ! Code developed by Jerry blackford and others at PML, based on pre-existing code. |
---|
62 | ! We accept no liability for errors or inaccuracies. |
---|
63 | ! See Zeebe & Wolf-Gladrow, 2001. CO2 in seawater: equilibrium, kinetics and isotopes. |
---|
64 | ! Elsevier Oceanography Series 65, 346. for a reasonable overview. |
---|
65 | ! Many other packages exist, replicating the same functionality in different languages. |
---|
66 | ! See http://cdiac.ornl.gov/oceans/co2rprt.html (CO2sys) |
---|
67 | ! or http://neon.otago.ac.nz/research/mfc/people/keith_hunter/software/swco2/ |
---|
68 | ! reference for prior usage of this code: Blackford & Gilbert, 2007. J Mar Sys 64, 229-241. |
---|
69 | ! |
---|
70 | ! Modifications |
---|
71 | ! 17/02/2010. Added conversion factor from per m3 to per kg (line 108-133) |
---|
72 | ! 17/02/2010. Update calculation of K1, K2, Kb to make consistant with the OCMIP protocols. |
---|
73 | ! 29/07/2011. Merged into MEDUSA with a raft of changes to this subroutine; less elsewhere |
---|
74 | ! 23/06/2015. Modified to take gas transfer velocity as an input (rather than wind speed); |
---|
75 | ! alter CO2 flux to /s rather than /d for consistency with other schemes |
---|
76 | ! |
---|
77 | ! Changes for MEDUSA include: |
---|
78 | ! - making the program a module |
---|
79 | ! - passing in input variables (obvious given the last point) |
---|
80 | ! - making alkalinity a state variable (rather than a function of salinity) |
---|
81 | ! |
---|
82 | IMPLICIT NONE |
---|
83 | |
---|
84 | REAL(wp), INTENT( in ) :: Temp ! degrees C |
---|
85 | REAL(wp), INTENT( in ) :: Sal ! PSU |
---|
86 | REAL(wp), INTENT( in ) :: DIC ! mmol / m3 |
---|
87 | REAL(wp), INTENT( in ) :: ALK ! meq / m3 |
---|
88 | REAL(wp), INTENT( in ) :: Depth ! m |
---|
89 | ! REAL(wp), INTENT( in ) :: Wnd ! m / s |
---|
90 | REAL(wp), INTENT( in ) :: xkw ! m / s |
---|
91 | REAL(wp), INTENT( in ) :: pCO2a ! uatm |
---|
92 | !---------------------------------------------------------------------- |
---|
93 | REAL(wp), INTENT( inout ) :: pH ! "normal" pH units |
---|
94 | REAL(wp), INTENT( inout ) :: pCO2w ! uatm |
---|
95 | REAL(wp), INTENT( inout ) :: h2co3 ! mmol / m3 |
---|
96 | REAL(wp), INTENT( inout ) :: hco3 ! mmol / m3 |
---|
97 | REAL(wp), INTENT( inout ) :: co3 ! mmol / m3 |
---|
98 | REAL(wp), INTENT( inout ) :: om_cal ! normalised |
---|
99 | REAL(wp), INTENT( inout ) :: om_arg ! normalised |
---|
100 | REAL(wp), INTENT( inout ) :: co2flux ! mmol / m2 / s |
---|
101 | REAL(wp), INTENT( inout ) :: TDIC ! umol / kg |
---|
102 | REAL(wp), INTENT( inout ) :: TALK ! ueq / kg |
---|
103 | REAL(wp), INTENT( inout ) :: dcf ! m3 / kg |
---|
104 | REAL(wp), INTENT( inout ) :: henry ! ? |
---|
105 | INTEGER, INTENT( inout ) :: iters ! # iterations to convergence |
---|
106 | !---------------------------------------------------------------------- |
---|
107 | ! REAL(wp) :: cco2,bicarb,carb,henry |
---|
108 | REAL(wp) :: cco2,bicarb,carb |
---|
109 | ! ====================================================================== |
---|
110 | |
---|
111 | ! write inputs to screen |
---|
112 | ! WRITE(*,*) " " |
---|
113 | ! WRITE(*,'(A28)') " .........Inputs........." |
---|
114 | ! WRITE(*,'(A24,F10.3)') " Temperature (C) = ", Temp |
---|
115 | ! WRITE(*,'(A24,F10.3)') " Salinity (psu) = ", Sal |
---|
116 | ! WRITE(*,'(A24,F10.3)') " Depth (m) = ", Depth |
---|
117 | ! WRITE(*,'(A24,F10.3)') " DIC (mmol/m3) = ", DIC |
---|
118 | ! WRITE(*,'(A24,F10.3)') " ALK (ueq/m3) = ", ALK |
---|
119 | ! WRITE(*,'(A24,F10.3)') " Wind speed (m/s) = ", Wnd |
---|
120 | ! WRITE(*,'(A24,F10.3)') " pCO2 atmos (uatm) = ", pCO2a |
---|
121 | |
---|
122 | ! AXY (07/05/13) ================================================== |
---|
123 | ! set total number of iterations to zero |
---|
124 | iters = 0 |
---|
125 | ! AXY (07/05/13) ================================================== |
---|
126 | |
---|
127 | call CO2_dynamics( Temp, Sal, Depth, DIC, ALK, pCO2a, & ! inputs |
---|
128 | pco2w, ph, h2co3, hco3, co3, henry, om_cal, om_arg, TDIC, TALK, & ! outputs |
---|
129 | dcf, iters) ! outputs |
---|
130 | |
---|
131 | ! AXY (18/08/11): only do air-sea exchange calculation if depth = 0 |
---|
132 | ! (i.e. surface calculations being performed) |
---|
133 | if (Depth .eq. 0.0) then |
---|
134 | call Air_sea_exchange( Temp, xkw, pCO2w, pCO2a, henry, dcf, & ! inputs |
---|
135 | co2flux ) ! output |
---|
136 | else |
---|
137 | co2flux = 0.0 |
---|
138 | endif |
---|
139 | |
---|
140 | ! AXY (09/01/14) ================================================== |
---|
141 | ! check for non-convergence and NaNs |
---|
142 | !! |
---|
143 | if (iters .eq. 25) then |
---|
144 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: ITERS WARNING, ', & |
---|
145 | iters, ' AT depth =', Depth |
---|
146 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: ztmp =', Temp |
---|
147 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: zsal =', Sal |
---|
148 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: zdic =', DIC |
---|
149 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: zalk =', ALK |
---|
150 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_kw660 =', xkw |
---|
151 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_ph =', ph |
---|
152 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_pco2w =', pCO2w |
---|
153 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_h2co3 =', h2co3 |
---|
154 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_hco3 =', hco3 |
---|
155 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_co3 =', co3 |
---|
156 | endif |
---|
157 | !! |
---|
158 | !! AXY (29/10/13): NaN checks |
---|
159 | if (co2flux /= co2flux) then |
---|
160 | !! if (lwp) write (numout,'(a,1pe10.2,4i6)') 'CO2FLUX-NAN', & |
---|
161 | !! & tmask(ji,jj,jk), kt, ji, jj, jk |
---|
162 | if (lwp) write (numout,'(a,a,f10.3,a,f10.3)') 'CO2FLUX-NAN', & |
---|
163 | & ' TMP', Temp, ' SAL', Sal |
---|
164 | if (lwp) write (numout,'(a,a,f10.3,a,f10.3)') 'CO2FLUX-NAN', & |
---|
165 | & ' DIC', DIC, ' ALK', ALK |
---|
166 | if (lwp) write (numout,'(a,a,f10.3,a,f10.3)') 'CO2FLUX-NAN', & |
---|
167 | & ' XKW', xkw, ' PH ', ph |
---|
168 | if (lwp) write (numout,'(a,a,i6)') 'CO2FLUX-NAN', & |
---|
169 | & ' ITERS', iters |
---|
170 | endif |
---|
171 | !! |
---|
172 | !! AXY (09/01/14): NaN fudges |
---|
173 | if (co2flux /= co2flux) then |
---|
174 | ph = 8.1 |
---|
175 | pCO2w = pCO2a |
---|
176 | h2co3 = 0.005 * DIC |
---|
177 | hco3 = 0.865 * DIC |
---|
178 | co3 = 0.130 * DIC |
---|
179 | om_cal = 4. |
---|
180 | om_arg = 2. |
---|
181 | co2flux = 0. |
---|
182 | TDIC = DIC |
---|
183 | TALK = ALK |
---|
184 | dcf = 1. |
---|
185 | henry = 1. |
---|
186 | endif |
---|
187 | ! AXY (09/01/14) ================================================== |
---|
188 | |
---|
189 | ! write outputs to screen |
---|
190 | ! WRITE(*,*) " " |
---|
191 | ! WRITE(*,'(A32,F10.3)') " ..........Outputs..........." |
---|
192 | ! WRITE(*,'(A27,F10.3)') " pH (pH) = ", pH |
---|
193 | ! WRITE(*,'(A27,F10.3)') " DIC (umol/kg) = ", TDIC |
---|
194 | ! WRITE(*,'(A27,F10.3)') " TALK (ueq/kg) = ", TALK |
---|
195 | ! WRITE(*,'(A27,F10.3)') " pco2w (uatm) = ", pco2w |
---|
196 | ! WRITE(*,'(A27,F10.3)') " carbonic acid (mmol/m3) = ", h2co3 |
---|
197 | ! WRITE(*,'(A27,F10.3)') " bicarbonate (mmol/m3) = ", bicarb |
---|
198 | ! WRITE(*,'(A27,F10.3)') " carbonate (mmol/m3) = ", carb |
---|
199 | ! WRITE(*,'(A27,F10.3)') " Omega calcite (~) = ", om_cal |
---|
200 | ! WRITE(*,'(A27,F10.3)') " Omega aragonite (~) = ", om_arg |
---|
201 | ! WRITE(*,'(A27,F10.3)') " air sea flux(mmol/m2/s) = ", flux |
---|
202 | ! WRITE(*,*) " " |
---|
203 | |
---|
204 | END SUBROUTINE trc_co2_medusa |
---|
205 | !----------------------------------------------------------------------- |
---|
206 | |
---|
207 | !======================================================================= |
---|
208 | !======================================================================= |
---|
209 | !======================================================================= |
---|
210 | |
---|
211 | !======================================================================= |
---|
212 | ! |
---|
213 | subroutine CO2_dynamics( T, S, Z, DIC, TALK, pco2a, & |
---|
214 | pco2w, ph, h2co3, bicarb, carb, henry, om_cal, om_arg, TCO2, TA, & |
---|
215 | dcf, iters ) |
---|
216 | ! |
---|
217 | !======================================================================= |
---|
218 | ! |
---|
219 | IMPLICIT NONE |
---|
220 | |
---|
221 | REAL(wp), INTENT( in ) :: T ! temperature (C) |
---|
222 | REAL(wp), INTENT( in ) :: S ! salinity (psu) |
---|
223 | REAL(wp), INTENT( in ) :: Z ! depth (metres) |
---|
224 | REAL(wp), INTENT( in ) :: DIC ! total dissolved inorganic carbon (mmol.m-3) |
---|
225 | REAL(wp), INTENT( in ) :: TALK ! total alkalinity (meq.m-3) |
---|
226 | REAL(wp), INTENT( in ) :: pco2a ! atmospheric pCO2 |
---|
227 | !---------------------------------------------------------------------- |
---|
228 | REAL(wp), INTENT( inout ) :: pco2w ! seawater pCO2 |
---|
229 | REAL(wp), INTENT( inout ) :: ph ! seawater pH |
---|
230 | REAL(wp), INTENT( inout ) :: h2co3 ! seawater carbonic acid (H2CO3) |
---|
231 | REAL(wp), INTENT( inout ) :: bicarb ! seawater bicarbonate ion (HCO3) |
---|
232 | REAL(wp), INTENT( inout ) :: carb ! seawater carbonate ion (CO3) |
---|
233 | REAL(wp), INTENT( inout ) :: henry ! Henry constant |
---|
234 | REAL(wp), INTENT( inout ) :: om_cal ! Omega calcite |
---|
235 | REAL(wp), INTENT( inout ) :: om_arg ! Omega aragonite |
---|
236 | REAL(wp), INTENT( inout ) :: TCO2 ! total dissolved inorganic carbon (mol.kg-1) |
---|
237 | REAL(wp), INTENT( inout ) :: TA ! total alkalinity (eq.kg-1) |
---|
238 | REAL(wp), INTENT( inout ) :: dcf ! density conversion factor |
---|
239 | INTEGER, INTENT( inout ) :: iters ! # iterations to convergence |
---|
240 | !---------------------------------------------------------------------- |
---|
241 | REAL(wp) :: a, b, c |
---|
242 | REAL(wp) :: ca, bc, cb |
---|
243 | REAL(wp) :: pco2water, fairco2 |
---|
244 | |
---|
245 | ! ====================================================================== |
---|
246 | |
---|
247 | ! Adjust to correct units. |
---|
248 | ! Haltafall uses mol/kg rather than umol/kg necessitating a scaling factor of /1.0D6 |
---|
249 | ! DIC (mmol/m3) needs to be converted to umol/kg via the calculation of water density at prevailing T&S |
---|
250 | ! sea-water density (Millero & Poisson, Deep-Sea Research, 1981, also known as UNESCO, 1981) |
---|
251 | ! with T: Temperature in degree Celsius; S: Salinity in practical units, density in kg/m3 |
---|
252 | ! valid for 0<T<40 and 0.5<S<43. Density Conversion factor (dcf) = * density * 1.0D3 |
---|
253 | |
---|
254 | a = 8.24493d-1 - 4.0899d-3*T + 7.6438d-5*T**2 - 8.2467d-7*T**3 + 5.3875d-9*T**4 |
---|
255 | b = -5.72466d-3 + 1.0227d-4*T - 1.6546d-6*T**2 |
---|
256 | c = 4.8314d-4 |
---|
257 | dcf = (999.842594 + 6.793952d-2*T- 9.095290d-3*T**2 + 1.001685d-4*T**3 & |
---|
258 | - 1.120083d-6*T**4 + 6.536332d-9*T**5+a*S+b*S**1.5+c*S**2)/1.0D3 |
---|
259 | |
---|
260 | TA = TALK / (1.0D6*dcf) |
---|
261 | TCO2 = DIC / (1.0D6*dcf) |
---|
262 | |
---|
263 | ! Call the parent routine for the carbonate system |
---|
264 | CALL CO2DYN ( TCO2, TA, T, S, pco2a, & ! inputs |
---|
265 | pco2water, pH, HENRY, ca, bc, cb, iters ) ! outputs |
---|
266 | |
---|
267 | ! Adjust outputs back to units used in the parent model code (e.g. mmol/m3) if appropriate |
---|
268 | |
---|
269 | pco2w = pco2water * (1.0D6) ! partial pressure of co2 in water |
---|
270 | TA = TA * (1.0D6) ! total alkalinity (umol/kg) |
---|
271 | h2co3 = ca * (1.0D6*dcf) ! carbonic acid concentration (mmol/m3) |
---|
272 | bicarb = bc * (1.0D6*dcf) ! bicarbonate ion concentration (mmol/m3) |
---|
273 | carb = cb * (1.0D6*dcf) ! carbonate ion concentration (mmol/m3) |
---|
274 | TCO2 = TCO2 * (1.0D6) ! total C or DIC in units of umol/kg |
---|
275 | |
---|
276 | ! Call carbonate saturation state subroutine to calculate calcite and aragonite calcification states |
---|
277 | |
---|
278 | CALL CaCO3_Saturation ( T, S, Z, cb, & ! inputs |
---|
279 | om_cal, om_arg ) ! outputs |
---|
280 | |
---|
281 | RETURN |
---|
282 | |
---|
283 | END SUBROUTINE CO2_dynamics |
---|
284 | !----------------------------------------------------------------------- |
---|
285 | |
---|
286 | !======================================================================= |
---|
287 | !======================================================================= |
---|
288 | !======================================================================= |
---|
289 | |
---|
290 | !======================================================================= |
---|
291 | ! |
---|
292 | SUBROUTINE Air_sea_exchange( T, xkw, pco2w, pco2a, henry, dcf, & |
---|
293 | flux ) |
---|
294 | ! |
---|
295 | !======================================================================= |
---|
296 | ! |
---|
297 | ! This routine should be called for the surface box only (like, duh) |
---|
298 | ! Uses the Nightingale and Liss parameterisation (GBC 14, 373-388 (2000). |
---|
299 | ! SC is the Schmidt number from Wanninkhof (1992, JGR 97,7373-7382), the Nightingale et al (2000) |
---|
300 | ! transfer velocity being for Sc=600 and the dependence of tranfer velocity on Sc coming |
---|
301 | ! from Jahne et al (1987, J. Geophys Res 92, 1937-1949). |
---|
302 | |
---|
303 | ! Inputs |
---|
304 | ! pCO2w partial pressure of CO2 in the water (from carbonate system subroutine call) |
---|
305 | ! pCO2a partial pressure of CO2 in the atmosphere (usually external forcing). |
---|
306 | ! T temperature (C) |
---|
307 | ! Wnd wind speed, metres (DELETED) |
---|
308 | ! xkw gas transfer velocity |
---|
309 | ! Henry henry's constant |
---|
310 | ! density the density of water for conversion between mmol/m3 and umol/kg |
---|
311 | |
---|
312 | ! Outputs are |
---|
313 | ! flux flux of CO2 in mmol C /m2/d |
---|
314 | ! +ve is in-gassing (air to sea), -ve is outgassing (sea to air). |
---|
315 | |
---|
316 | IMPLICIT NONE |
---|
317 | |
---|
318 | REAL(wp), INTENT( in ) :: T, xkw, pco2w, pco2a, henry, dcf ! INPUT PARAMETERS: |
---|
319 | !----------------------------------------------------------------------- |
---|
320 | REAL(wp), INTENT( inout ) :: flux ! OUTPUT Variables |
---|
321 | !----------------------------------------------------------------------- |
---|
322 | REAL(wp) :: sc, fwind ! LOCAL VARIABLES: |
---|
323 | |
---|
324 | ! calculate the Schmidt number and unit conversions |
---|
325 | sc = 2073.1-125.62*T+3.6276*T**2.0-0.0432190*T**3.0 |
---|
326 | ! fwind = (0.222d0 * wnd**2d0 + 0.333d0 * wnd)*(sc/660.d0)**(-0.5) |
---|
327 | fwind = xkw * (sc/660.d0)**(-0.5) |
---|
328 | fwind = fwind*24.d0/100.d0 ! convert to m/day |
---|
329 | |
---|
330 | ! flux depends on the difference in partial pressures, wind and henry |
---|
331 | ! here it is rescaled to mmol/m2/d |
---|
332 | flux = fwind * henry * ( pco2a - pco2w ) * dcf |
---|
333 | |
---|
334 | ! AXY (23/06/15): let's get it from /d to /s |
---|
335 | flux = flux / ( 86400. ) |
---|
336 | |
---|
337 | RETURN |
---|
338 | |
---|
339 | END SUBROUTINE Air_sea_exchange |
---|
340 | !----------------------------------------------------------------------- |
---|
341 | |
---|
342 | !======================================================================= |
---|
343 | !======================================================================= |
---|
344 | !======================================================================= |
---|
345 | |
---|
346 | !======================================================================= |
---|
347 | ! |
---|
348 | SUBROUTINE CO2dyn ( TCO2, TA, T, S, PCO2A, PCO2, PH, HENRY, ca, bc, cb, iters ) |
---|
349 | ! |
---|
350 | !======================================================================= |
---|
351 | ! |
---|
352 | ! This subroutine acts as an interface to the Haltafall iteration, setting options etc. |
---|
353 | |
---|
354 | REAL(wp), INTENT( in ) :: TCO2, TA, T, S, PCO2A |
---|
355 | !----------------------------------------------------------------------- |
---|
356 | REAL(wp), INTENT( inout ) :: PCO2, PH, HENRY, ca, bc, cb |
---|
357 | INTEGER, INTENT( inout ) :: iters |
---|
358 | !----------------------------------------------------------------------- |
---|
359 | REAL(wp) :: PRSS |
---|
360 | INTEGER :: MCONC, MKVAL, ICONST, ICALC |
---|
361 | |
---|
362 | PARAMETER ( MCONC = 9,MKVAL = 4 ) |
---|
363 | |
---|
364 | REAL(wp), DIMENSION(MKVAL) :: AKVAL |
---|
365 | REAL(wp), DIMENSION(MCONC) :: CONCS |
---|
366 | |
---|
367 | ICONST = 6 |
---|
368 | PRSS = 1.0D0 |
---|
369 | CONCS(1) = TCO2 |
---|
370 | CONCS(2) = TA |
---|
371 | ICALC = 1 |
---|
372 | |
---|
373 | CALL POLYCO(PRSS,T,S,CONCS,MCONC,AKVAL,MKVAL,ICALC,ICONST,iters) |
---|
374 | |
---|
375 | if(iters.eq.25) then |
---|
376 | ! AXY (13/05/13): this location has failed to converge; use |
---|
377 | ! some standard values to plug the hole |
---|
378 | CONCS(3) = PCO2A * 1E-06 ! atmospheric value |
---|
379 | CONCS(4) = 8.1 ! global average pH |
---|
380 | CONCS(5) = TCO2 * 0.005 ! some "standard" values plucked |
---|
381 | CONCS(6) = TCO2 * 0.865 ! from pages 5-6 of Zeebe & Wolf- |
---|
382 | CONCS(7) = TCO2 * 0.130 ! Gladow (2001) |
---|
383 | AKVAL(1) = 0.5E-01 ! trial-and-error value |
---|
384 | IF(lwp) THEN |
---|
385 | WRITE(numout,*) ' CO2DYN : reset failed outputs' |
---|
386 | ENDIF |
---|
387 | endif |
---|
388 | |
---|
389 | PCO2 = CONCS(3) |
---|
390 | PH = CONCS(4) |
---|
391 | ca = CONCS(5) |
---|
392 | bc = CONCS(6) |
---|
393 | cb = CONCS(7) |
---|
394 | HENRY = AKVAL(1) |
---|
395 | |
---|
396 | RETURN |
---|
397 | |
---|
398 | END SUBROUTINE |
---|
399 | !----------------------------------------------------------------------- |
---|
400 | |
---|
401 | !======================================================================= |
---|
402 | !======================================================================= |
---|
403 | !======================================================================= |
---|
404 | |
---|
405 | !======================================================================= |
---|
406 | ! |
---|
407 | SUBROUTINE POLYCO(PD,TD,SD,CONCS,NCONC,AKVAL,NKVAL,ICALC,ICONST,iters) |
---|
408 | ! |
---|
409 | !======================================================================= |
---|
410 | ! |
---|
411 | ! MASTER SUBROUTINE FOR CALCULATION OF THE CO2 SYSTEM THERMODYNAMICS |
---|
412 | ! |
---|
413 | ! EXPLANATION OF POLYCO PARAMETERS |
---|
414 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
415 | ! P - PRESSURE IN ATMOSPHERES (P<>1 CODED only for ICONST=3 and ICONST=6) |
---|
416 | ! T - TEMPERATURE IN DEG.C |
---|
417 | ! S - SALINITY IN PPT |
---|
418 | ! CONCS(1) - TOTAL C (MOL/KG) |
---|
419 | ! CONCS(2) - TOTAL ALKALINITY (MOL/KG) |
---|
420 | ! CONCS(3) - PCO2 (ATM) |
---|
421 | ! CONCS(4) - PH |
---|
422 | ! CONCS(5) - {H2CO3} (MOL/KG) |
---|
423 | ! CONCS(6) - {HCO3} (MOL/KG) |
---|
424 | ! CONCS(7) - {CO3} (MOL/KG) |
---|
425 | ! CONCS(8) - CARBONATE ALKALINITY ) FOR ICONST = 4,5,6 |
---|
426 | ! CONCS(9) - BORATE ALKALINITY ) ONLY |
---|
427 | ! NCONC - SIZE OF CONCS ARRAY (7 FOR ICONST=1,2,3; 9 FOR ICONST |
---|
428 | ! AKVAL(1) - KP (HENRY'S LAW CONSTANT) (MOL/KG/ATM) |
---|
429 | ! AKVAL(2) - K1C (H2CO3 DISSOCIATION) (MOL/KG) |
---|
430 | ! AKVAL(3) - K2C (HCO3 DISSOCIATION) (MOL/KG) |
---|
431 | ! AKVAL(4) - KB (B(OH)3 DISSOCIATION) (MOL/KG) FOR ICONST=4,5,6 |
---|
432 | ! NKVAL - SIZE OF AKVAL ARRAY (3 FOR ICONST=1,2,3; 4 FOR ICONST= |
---|
433 | ! ICALC - SELECTION OF THE TWO INPUT PARAMETERS: |
---|
434 | ! ICALC = 1 TOTAL C AND ALKALINITY |
---|
435 | ! ICALC = 2 TOTAL C AND PCO2 |
---|
436 | ! ICALC = 3 TOTAL C AND PH |
---|
437 | ! ICALC = 4 ALKALINITY AND PCO2 |
---|
438 | ! ICALC = 5 ALKALINITY AND PH |
---|
439 | ! ICALC = 6 PCO2 AND PH |
---|
440 | ! ICALC = 7 CALCULATE CONSTANTS AKVAL ONLY |
---|
441 | ! ICONST - SELECTION OF PH SCALE AND COMPONENTS: |
---|
442 | ! ICONST = 1 NBS PH SCALE |
---|
443 | ! ICONST = 2 HANSSON'S SCALE (SWS WITHOUT FLUORIDE) |
---|
444 | ! ICONST = 3 SWS PH SCALE |
---|
445 | ! ICONST = 4 AS 1 BUT INCLUDING BORATE IN THE CALCULATION |
---|
446 | ! ICONST = 5 AS 2 BUT INCLUDING BORATE IN THE CALCULATION |
---|
447 | ! ICONST = 6 AS 3 BUT INCLUDING BORATE IN THE CALCULATION |
---|
448 | |
---|
449 | ! NOTE: FOR ICONST=1,2,3 CONCS(2) REPRESENTS CARBONATE ALKALINITY SINC |
---|
450 | ! BORATE IS NOT INCLUDED IN THE CALCULATION. FOR ICONST=4,5,6 CO |
---|
451 | ! REPRESENTS TOTAL ALKALINITY (CARBONATE + BORATE), THE COMPONEN |
---|
452 | ! WHICH ARE GIVEN IN CONCS(8) AND CONCS(9) |
---|
453 | |
---|
454 | REAL(wp) :: PMIN, PMAX, SMIN, SMAX, TMIN, TMAX, & |
---|
455 | & PD, TD, SD, P, T, S, BTOT |
---|
456 | INTEGER :: MINJC, MAXJC, MINJK, MAXJK, MINCAL, MAXCAL, MINCON, & |
---|
457 | & MAXCON, NCONC, NKVAL, ICALC, ICONST, IC, iters |
---|
458 | LOGICAL :: BORON |
---|
459 | |
---|
460 | PARAMETER(MINJC=7,MAXJC=9,MINJK=3,MAXJK=4) |
---|
461 | PARAMETER(MINCAL=1,MAXCAL=7,MINCON=1,MAXCON=6) |
---|
462 | !! AXY (11/08/11): TMIN changed to -2 to stop error messages in polar regions |
---|
463 | !! AXY (09/01/14): TMIN changed to -3 to stop error messages in polar regions |
---|
464 | !! AXY (03/03/14): SMAX changed to 42 to stop error messages in salty regions (the Great Answer ...) |
---|
465 | !! AXY (05/03/14): SMAX changed to 45 to stop error messages in salty regions |
---|
466 | PARAMETER(PMIN=0.99999D0,PMAX=1.00001D0,SMIN=0.0D0, & |
---|
467 | & SMAX=45.0D0,TMIN=-3.0D0,TMAX=40.0D0) |
---|
468 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
469 | REAL(wp), DIMENSION(NCONC) :: CONCS |
---|
470 | |
---|
471 | P = PD |
---|
472 | S = SD |
---|
473 | T = TD |
---|
474 | |
---|
475 | !! AXY (09/08/11): STOP statements commented out and WARNING messages added |
---|
476 | IF(lwp) THEN |
---|
477 | IF(T.LT.TMIN.OR.T.GT.TMAX) WRITE(numout,*) ' trc_co2_medusa: T WARNING, ', T, TMIN, TMAX, S, P |
---|
478 | IF(S.LT.SMIN.OR.S.GT.SMAX) WRITE(numout,*) ' trc_co2_medusa: S WARNING, ', S, SMIN, SMAX, T, P |
---|
479 | IF(P.LT.PMIN.OR.P.GT.PMAX) WRITE(numout,*) ' trc_co2_medusa: P WARNING, ', P, PMIN, PMAX, T, S |
---|
480 | ENDIF |
---|
481 | ! IF(T.LT.TMIN.OR.T.GT.TMAX)WRITE (*,*) P, S, T, TMIN, TMAX |
---|
482 | ! IF(P.LT.PMIN.OR.P.GT.PMAX) STOP'POLYCO - PRESSURE OUT OF RANGE' |
---|
483 | ! IF(S.LT.SMIN.OR.S.GT.SMAX) STOP'POLYCO - SALINITY OUT OF RANGE' |
---|
484 | ! IF(T.LT.TMIN.OR.T.GT.TMAX) STOP'POLYCO - TEMP. OUT OF RANGE' |
---|
485 | |
---|
486 | !! AXY (17/04/13): iMARNET climate change simulation appears to be compromised |
---|
487 | !! by excessive Caspian Sea salinity (> 45 PSU); this may be |
---|
488 | !! a runoff glitch that doesn't manifest on local instance of |
---|
489 | !! NEMO; fudging this here by replacing excessively high |
---|
490 | !! salinity values with maximum salinity value for the purpose |
---|
491 | !! of carbonate chemistry calculations; if this works, a more |
---|
492 | !! sensible solution might be to kill Caspian, etc. |
---|
493 | IF(lwp) THEN |
---|
494 | IF(S.LT.SMIN) THEN |
---|
495 | WRITE(numout,*) ' trc_co2_medusa: S RESET, ', S, '->', SMIN |
---|
496 | S = SMIN |
---|
497 | ENDIF |
---|
498 | IF(S.GT.SMAX) THEN |
---|
499 | WRITE(numout,*) ' trc_co2_medusa: S RESET, ', S, '->', SMAX |
---|
500 | S = SMAX |
---|
501 | ENDIF |
---|
502 | ENDIF |
---|
503 | |
---|
504 | !! AXY (28/02/13): as above, but for T; this has been done for 1/4-degree |
---|
505 | !! simulations to provide a temporary fix to a strange |
---|
506 | !! glitch in ocean temperature; basically, T has been found |
---|
507 | !! to spike upwards arbitrarily (e.g. 24C -> 53C); this |
---|
508 | !! causes the carbonate chemistry calculations to fail; |
---|
509 | !! this fix aims to stop the problem so that we can estimate |
---|
510 | !! how frequent a problem this is; it's suggestive of a |
---|
511 | !! memory leak |
---|
512 | IF(lwp) THEN |
---|
513 | IF(T.LT.TMIN) THEN |
---|
514 | WRITE(numout,*) ' trc_co2_medusa: T RESET, ', T, '->', TMIN |
---|
515 | T = TMIN |
---|
516 | ENDIF |
---|
517 | IF(T.GT.TMAX) THEN |
---|
518 | WRITE(numout,*) ' trc_co2_medusa: T RESET, ', T, '->', TMAX |
---|
519 | T = TMAX |
---|
520 | ENDIF |
---|
521 | ENDIF |
---|
522 | |
---|
523 | IF(ICALC.LT.MINCAL.OR.ICALC.GT.MAXCAL) STOP 'POLYCO - ICALC OUT OR RANGE' |
---|
524 | IF(ICONST.LT.MINCON.OR.ICONST.GT.MAXCON) STOP 'POLYCO - ICONST OUT OF RANGE' |
---|
525 | BORON=(ICONST.GT.3) |
---|
526 | IF(BORON) THEN |
---|
527 | IC=ICONST-3 |
---|
528 | BTOT=0.0004128D0*S/35.0D0 |
---|
529 | IF(NCONC.NE.MAXJC) STOP 'POLYCO - WRONG NCONC VALUE' |
---|
530 | IF(NKVAL.NE.MAXJK) STOP 'POLYCO - WRONG NKVAL VALUE' |
---|
531 | ELSE |
---|
532 | IC=ICONST |
---|
533 | IF(NCONC.NE.MINJC) STOP 'POLYCO - WRONG NCONC VALUE' |
---|
534 | IF(NKVAL.NE.MINJK) STOP 'POLYCO - WRONG NKVAL VALUE' |
---|
535 | ENDIF |
---|
536 | |
---|
537 | CALL CO2SET(P,T,S,AKVAL,NKVAL,IC) |
---|
538 | |
---|
539 | IF(ICALC.LT.MAXCAL) & |
---|
540 | & CALL CO2CLC(CONCS,NCONC,AKVAL,NKVAL,ICALC,BORON,BTOT,iters) |
---|
541 | |
---|
542 | RETURN |
---|
543 | |
---|
544 | END SUBROUTINE |
---|
545 | !----------------------------------------------------------------------- |
---|
546 | |
---|
547 | !======================================================================= |
---|
548 | !======================================================================= |
---|
549 | !======================================================================= |
---|
550 | |
---|
551 | !======================================================================= |
---|
552 | ! |
---|
553 | SUBROUTINE CO2SET(P,T,S,AKVAL,NKVAL,IC) |
---|
554 | ! |
---|
555 | !======================================================================= |
---|
556 | ! |
---|
557 | ! Routine to calculate CO2 system constants under the conditions set by |
---|
558 | ! P,S,T (NOTE: PRESSURE <> 1ATM IS NOT YET CODED) |
---|
559 | |
---|
560 | ! I. Calculate constants at P=1 and S=0 using |
---|
561 | |
---|
562 | ! ln K0 = A + B/TK + C ln TK |
---|
563 | ! (where TK is in Kelvin) |
---|
564 | |
---|
565 | ! II. Calculate constants at P=1 and salinity S using |
---|
566 | |
---|
567 | ! ln K = ln K0 + (a0 + a1/TK + a2 ln TK) S**1/2 |
---|
568 | ! + (b0 + b1TK + b2TK**2) S |
---|
569 | |
---|
570 | ! The sources of the coefficients are as follows: |
---|
571 | |
---|
572 | ! IC= 1 2 3 |
---|
573 | ! (NBS pH scale) (SWS pH scale (SWS pH scale |
---|
574 | ! with no F) |
---|
575 | |
---|
576 | ! KP WEISS (1974) WEISS(1974) WEISS(1974 |
---|
577 | |
---|
578 | ! K1C ) MEHRBACH ACC. TO HANSSON ACC. TO HANSSON AND MEH |
---|
579 | ! K2C ) MILLERO (1979) MILLERO (1979) ACC. TO DICKS |
---|
580 | ! KB ) AND MILLERO (1 |
---|
581 | ! (K1C AND K2C |
---|
582 | ! HANSSON ACC |
---|
583 | ! MILLERO (1 |
---|
584 | ! (KB ONLY |
---|
585 | |
---|
586 | ! *** |
---|
587 | ! IMPLICIT real*8 (A-H,O-Z) |
---|
588 | |
---|
589 | ! Modified by jcb 17/02/10 to use OCMIP calculations of K1, K2, Kb. |
---|
590 | ! Differences are subtle rather than significant |
---|
591 | |
---|
592 | INTEGER :: MAXK, MAXCON, NKVAL, ICON, IC, IK |
---|
593 | ! *** |
---|
594 | PARAMETER(MAXK=4,MAXCON=3) |
---|
595 | REAL(wp), DIMENSION(MAXK) :: A, B, C |
---|
596 | REAL(wp), DIMENSION(MAXK,MAXCON) :: A0, A1, A2 |
---|
597 | REAL(wp), DIMENSION(MAXK,MAXCON) :: B0, B1, B2 |
---|
598 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
599 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
600 | ! REAL(wp) :: P,T,S,VAL,TK |
---|
601 | REAL(wp) :: P,T,S,VAL,TK, delta, kappa, Rgas |
---|
602 | REAL(wp) :: dlogTK, S2, sqrtS, S15, k1, k2, kb |
---|
603 | ! *** |
---|
604 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
605 | DATA Rgas/83.131/ |
---|
606 | |
---|
607 | DATA A/-167.8108D0, 290.9097D0, 207.6548D0, 148.0248D0/ |
---|
608 | DATA B/9345.17D0, -14554.21D0, -11843.79D0, -8966.9D0/ |
---|
609 | DATA C/23.3585D0, -45.0575D0, -33.6485D0, -24.4344D0/ |
---|
610 | DATA (A0(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
611 | DATA (A0(2,ICON),ICON=1,MAXCON) /0.0221D0, 0.5709D0, -45.8076D0/ |
---|
612 | DATA (A0(3,ICON),ICON=1,MAXCON) /0.9805D0, 1.4853D0, -39.5492D0/ |
---|
613 | DATA (A0(4,ICON),ICON=1,MAXCON) /0.0473D0, 0.5998D0, 0.5998D0/ |
---|
614 | DATA (A1(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
615 | DATA (A1(2,ICON),ICON=1,MAXCON) /34.02D0, -84.25D0, 1935.07D0/ |
---|
616 | DATA (A1(3,ICON),ICON=1,MAXCON) /-92.65D0, -192.69D0, 1590.14D0/ |
---|
617 | DATA (A1(4,ICON),ICON=1,MAXCON) /49.10D0, -75.25D0, -75.25D0/ |
---|
618 | DATA (A2(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
619 | DATA (A2(2,ICON),ICON=1,MAXCON) /2*0.0D0,6.9513D0/ |
---|
620 | DATA (A2(3,ICON),ICON=1,MAXCON) /2*0.0D0,6.1523D0/ |
---|
621 | DATA (A2(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
622 | DATA (B0(1,ICON),ICON=1,MAXCON) /3*0.023517D0/ |
---|
623 | DATA (B0(2,ICON),ICON=1,MAXCON) /0.0D0,-0.01632D0,-0.01566D0/ |
---|
624 | DATA (B0(3,ICON),ICON=1,MAXCON) /-0.03294D0,-0.05058D0,-0.04997D0/ |
---|
625 | DATA (B0(4,ICON),ICON=1,MAXCON) /0.0D0, -0.01767D0, -0.01767D0/ |
---|
626 | DATA (B1(1,ICON),ICON=1,MAXCON) /3*-2.3656D-4/ |
---|
627 | DATA (B1(2,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
628 | DATA (B1(3,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
629 | DATA (B1(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
630 | DATA (B2(1,ICON),ICON=1,MAXCON) /3*4.7036D-7/ |
---|
631 | DATA (B2(2,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
632 | DATA (B2(3,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
633 | DATA (B2(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
634 | |
---|
635 | TK=T+273.15D0 |
---|
636 | DO 100 IK=1,NKVAL |
---|
637 | VAL=A(IK) + B(IK)/TK + C(IK)*LOG(TK) |
---|
638 | VAL=VAL + (A0(IK,IC) + A1(IK,IC)/TK + A2(IK,IC)*LOG(TK))*SQRT(S) |
---|
639 | VAL=VAL + (B0(IK,IC) + B1(IK,IC)*TK + B2(IK,IC)*TK*TK)*S |
---|
640 | AKVAL(IK)=EXP(VAL) |
---|
641 | 100 CONTINUE |
---|
642 | |
---|
643 | IF (IC .EQ. 3) THEN |
---|
644 | ! Calculation of constants as used in the OCMIP process for ICONST = 3 or 6 |
---|
645 | ! see http://www.ipsl.jussieu.fr/OCMIP/ |
---|
646 | ! added jcb 17/02/10 |
---|
647 | |
---|
648 | ! Derive simple terms used more than once |
---|
649 | dlogTK = log(TK) |
---|
650 | S2 = S*S |
---|
651 | sqrtS = sqrt(S) |
---|
652 | S15 = S**1.5 |
---|
653 | ! k1 = [H][HCO3]/[H2CO3] |
---|
654 | ! k2 = [H][CO3]/[HCO3] |
---|
655 | ! Millero p.664 (1995) using Mehrbach et al. data on seawater scale |
---|
656 | k1=10**(-1*(3670.7/TK - 62.008 + 9.7944*dlogTK - & |
---|
657 | & 0.0118 * S + 0.000116*S2)) |
---|
658 | k2=10**(-1*(1394.7/TK + 4.777 - & |
---|
659 | & 0.0184*S + 0.000118*S2)) |
---|
660 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
661 | ! Correction for high pressure (from Millero 1995) |
---|
662 | ! added YA 04/10/2010 |
---|
663 | delta=-25.5+0.1271*T |
---|
664 | kappa=(-3.08+0.0877*T)/1000.0 |
---|
665 | k1=k1*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
666 | delta=-15.82-0.0219*T |
---|
667 | kappa=(1.13-0.1475*T)/1000.0 |
---|
668 | k2=k2*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
669 | |
---|
670 | ! kb = [H][BO2]/[HBO2] |
---|
671 | ! Millero p.669 (1995) using data from Dickson (1990) |
---|
672 | kb=exp((-8966.90 - 2890.53*sqrtS - 77.942*S + & |
---|
673 | & 1.728*S15 - 0.0996*S2)/TK + & |
---|
674 | & (148.0248 + 137.1942*sqrtS + 1.62142*S) + & |
---|
675 | & (-24.4344 - 25.085*sqrtS - 0.2474*S) * & |
---|
676 | & dlogTK + 0.053105*sqrtS*TK) |
---|
677 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
678 | ! Correction for high pressure (from Millero, 1995) |
---|
679 | ! added YA 04/10/2010 |
---|
680 | delta=-29.48+0.1622*T-0.002608*T**2.0 |
---|
681 | kappa=-2.84/1000.0 |
---|
682 | kb=kb*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
683 | ! Correction for high pressure of Kw (from Millero 1995) |
---|
684 | ! added YA 04/10/2010 |
---|
685 | delta=-25.60+0.2324*T-0.0036246*T**2 |
---|
686 | kappa=(-5.13+0.0794*T)/1000.0 |
---|
687 | AKVAL(1)=AKVAL(1)*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
688 | |
---|
689 | ! replace haltafall calculations with OCMIP calculations |
---|
690 | AKVAL(2) = k1 |
---|
691 | AKVAL(3) = k2 |
---|
692 | AKVAL(4) = kb |
---|
693 | END IF ! section implimenting OCMIP coefficients |
---|
694 | |
---|
695 | RETURN |
---|
696 | END SUBROUTINE |
---|
697 | !----------------------------------------------------------------------- |
---|
698 | |
---|
699 | !======================================================================= |
---|
700 | !======================================================================= |
---|
701 | !======================================================================= |
---|
702 | |
---|
703 | !======================================================================= |
---|
704 | ! |
---|
705 | SUBROUTINE CO2CLC(CONCS,NCONC,AKVAL,NKVAL,ICALC,BORON,BTOT,iters) |
---|
706 | ! |
---|
707 | !======================================================================= |
---|
708 | ! |
---|
709 | ! ROUTINE TO CARRY OUT CO2 CALCULATIONS WITH 2 FIXED PARAMETERS ACCORDI |
---|
710 | ! THE EQUATIONS GIVEN BY PARKS(1969) AND SKIRROW (1975) |
---|
711 | ! WITH ADDITIONS FOR INCLUDING BORON IF BORON=.TRUE. |
---|
712 | |
---|
713 | |
---|
714 | ! IMPLICIT real*8 (A-H,O-Z) |
---|
715 | INTEGER :: NCONC, NKVAL, ICALC, II, KARL, LQ, iters |
---|
716 | ! AXY (07/05/13) ================================================== |
---|
717 | ! put counter in to check duration in convergence loop |
---|
718 | INTEGER :: COUNTER,C_CHECK,C_SW,III |
---|
719 | ! AXY (07/05/13) ================================================== |
---|
720 | REAL(wp) :: CTOT,ALK,PCO2,PH,H2CO3,HCO3,CO3,ALKC |
---|
721 | REAL(wp) :: ALKB,AKP,AK1C,AK2C,AKB,BTOT |
---|
722 | REAL(wp) :: AKR,AHPLUS |
---|
723 | REAL(wp) :: PROD,tol1,tol2,tol3,tol4,steg,fak |
---|
724 | REAL(wp) :: STEGBY,Y,X,W,X1,Y1,X2,Y2,FACTOR,TERM,Z |
---|
725 | REAL(wp), DIMENSION(NCONC) :: CONCS |
---|
726 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
727 | REAL(wp), DIMENSION(9) :: CONCS2 |
---|
728 | REAL(wp), DIMENSION(4) :: AKVAL2 |
---|
729 | |
---|
730 | ! AXY (07/05/13) ================================================== |
---|
731 | ! setup array to store old values of concs |
---|
732 | real(wp),DIMENSION(:) :: old_CONCS(NCONC) |
---|
733 | ! AXY (07/05/13) ================================================== |
---|
734 | |
---|
735 | EQUIVALENCE (CTOT , CONCS2(1)), (ALK , CONCS2(2)), & |
---|
736 | & (PCO2 , CONCS2(3)), (PH , CONCS2(4)), & |
---|
737 | & (H2CO3 , CONCS2(5)), (HCO3 , CONCS2(6)), & |
---|
738 | & (CO3 , CONCS2(7)), (ALKC , CONCS2(8)), & |
---|
739 | & (ALKB , CONCS2(9)), & |
---|
740 | & (AKP , AKVAL2(1)), (AK1C , AKVAL2(2)), & |
---|
741 | & (AK2C , AKVAL2(3)), (AKB , AKVAL2(4)) |
---|
742 | LOGICAL :: BORON,DONE |
---|
743 | |
---|
744 | ! AXY (07/05/13) ================================================== |
---|
745 | ! DERIVING PH REQUIRES FOLLOWING LOOP TO CONVERGE. |
---|
746 | ! THIS SUBROUTINE RELIES ON CONVERGENCE. IF THE ENVIRONMENTAL |
---|
747 | ! CONDITIONS DO NOT ALLOW FOR CONVERGENCE (IN 3D MODEL THIS IS |
---|
748 | ! LIKELY TO OCCUR NEAR LOW SALINITY REGIONS) THE MODEL WILL |
---|
749 | ! BE STUCK IN THE LOOP. TO AVOID THIS A CONVERGENCE CONDITION |
---|
750 | ! IS PUT IN PLACE TO SET A FLAGG OF -99 IN THE PH VAR FOR NON CONVEGENCE. |
---|
751 | ! THE MODEL IS THEN ALLOWED TO CONTINUE. 'COUNTER, C_SW,C_CHECK' ARE |
---|
752 | ! THE LOCAL VARS USED. |
---|
753 | ! C_SW = condition of convergence 0=yes, 1= no |
---|
754 | ! COUNTER = number of iterations |
---|
755 | ! C_CHECK = maximum number of iterations |
---|
756 | |
---|
757 | ! SET COUNTER AND SWITCH TO ZERO AND OFF |
---|
758 | COUNTER=0 |
---|
759 | C_SW=0 |
---|
760 | ! FROM EXPERIENCE IF THE ITERATIONS IN THE FOLLOWING DO LOOP |
---|
761 | ! EXCEEDS 15 CONVERGENCE WILL NOT OCCUR. THE OVERHEAD OF 25 ITERATIONS |
---|
762 | ! IS OK FOR SMALL DOMAINS WITH 1/10 AND 1/15 DEG RESOLUTION. |
---|
763 | ! I RECOMMEND A LOWER VALUE OF 15 FOR HIGHER RESOLUTION OR LARGER DOMAINS. |
---|
764 | C_CHECK=25 |
---|
765 | ! AXY (07/05/13) ================================================== |
---|
766 | |
---|
767 | DO 100 II=1,NCONC |
---|
768 | CONCS2(II)=CONCS(II) |
---|
769 | |
---|
770 | ! AXY (07/05/13) ================================================== |
---|
771 | ! IF CONVERGENCE IS NOT ACHIEVED THE LOCAL ARRAY CONCS MUST BE STORED TO |
---|
772 | ! ALLOW THE MODEL TO CONTINUE. THEREFORE .... |
---|
773 | ! UPDATE OLD_CONCS |
---|
774 | old_CONCS(II)=CONCS(II) |
---|
775 | ! AXY (07/05/13) ================================================== |
---|
776 | |
---|
777 | 100 CONTINUE |
---|
778 | DO 110 II=1,NKVAL |
---|
779 | AKVAL2(II)=AKVAL(II) |
---|
780 | 110 CONTINUE |
---|
781 | AKR = AK1C/AK2C |
---|
782 | AHPLUS=10.0D0**(-PH) |
---|
783 | PROD=AKR*AKP*PCO2 |
---|
784 | |
---|
785 | IF(BORON) THEN |
---|
786 | |
---|
787 | IF(ICALC.EQ.1.OR.ICALC.EQ.4) THEN |
---|
788 | ! *** ALK, BTOT AND CTOT OR PCO2 FIXED *** |
---|
789 | ! *** ITERATIVE CALCULATION NECESSARY HERE |
---|
790 | |
---|
791 | ! SET INITIAL GUESSES AND TOLERANCE |
---|
792 | H2CO3=PCO2*AKP |
---|
793 | CO3=ALK/10.0D0 |
---|
794 | AHPLUS=1.0D-8 |
---|
795 | ALKB=BTOT |
---|
796 | TOL1=ALK/1.0D5 |
---|
797 | TOL2=H2CO3/1.0D5 |
---|
798 | TOL3=CTOT/1.0D5 |
---|
799 | TOL4=BTOT/1.0D5 |
---|
800 | |
---|
801 | ! HALTAFALL iteration to determine CO3, ALKB, AHPLUS |
---|
802 | KARL=1 |
---|
803 | STEG=2.0D0 |
---|
804 | FAK=1.0D0 |
---|
805 | STEGBY=0.4D0 |
---|
806 | 10 DONE=.TRUE. |
---|
807 | |
---|
808 | ! AXY (07/05/13) ================================================== |
---|
809 | ! SET COUNTER UPDATE. FLAG 10 IS THE POINT OF RETURN FOR |
---|
810 | ! THE CONVERGENCE CONDITION |
---|
811 | COUNTER=COUNTER+1 |
---|
812 | iters = COUNTER |
---|
813 | ! CHECK IF CONVERGENCE HAS OCCURED IN THE NUMBER OF |
---|
814 | ! ACCEPTABLE ITTERATIONS. SET C_SW TO LET MODEL KNOW |
---|
815 | ! WHAT TO DO AT THE END OF THE SUBROUTINE |
---|
816 | if(counter.ge.c_check)then |
---|
817 | IF(lwp) THEN |
---|
818 | WRITE(numout,*) ' CO2CLC : ITERS WARNING, ', iters |
---|
819 | ENDIF |
---|
820 | c_sw=1 |
---|
821 | GOTO 123 |
---|
822 | endif |
---|
823 | ! AXY (07/05/13) ================================================== |
---|
824 | |
---|
825 | IF(ICALC.EQ.4) THEN |
---|
826 | ! *** PCO2 IS FIXED *** |
---|
827 | Y=AHPLUS*AHPLUS*CO3/(AK1C*AK2C) |
---|
828 | IF(ABS(Y-H2CO3).GT.TOL2) THEN |
---|
829 | CO3=CO3*H2CO3/Y |
---|
830 | DONE=.FALSE. |
---|
831 | ENDIF |
---|
832 | ELSEIF(ICALC.EQ.1) THEN |
---|
833 | ! *** CTOT IS FIXED *** |
---|
834 | Y=CO3*(1.0D0+AHPLUS/AK2C+AHPLUS*AHPLUS/(AK1C*AK2C)) |
---|
835 | IF(ABS(Y-CTOT).GT.TOL3) THEN |
---|
836 | CO3=CO3*CTOT/Y |
---|
837 | DONE=.FALSE. |
---|
838 | ENDIF |
---|
839 | ENDIF |
---|
840 | Y=ALKB*(1.0D0+AHPLUS/AKB) |
---|
841 | IF(ABS(Y-BTOT).GT.TOL4) THEN |
---|
842 | ALKB=ALKB*BTOT/Y |
---|
843 | DONE=.FALSE. |
---|
844 | ENDIF |
---|
845 | |
---|
846 | ! Alkalinity is equivalent to -(total H+), so the sign of W is opposite |
---|
847 | ! to that normally used |
---|
848 | |
---|
849 | Y=CO3*(2.0D0+AHPLUS/AK2C)+ALKB |
---|
850 | IF(ABS(Y-ALK).GT.TOL1) THEN |
---|
851 | DONE=.FALSE. |
---|
852 | X=LOG(AHPLUS) |
---|
853 | ! W=SIGN(1.0D0,Y-ALK) |
---|
854 | IF ( Y-ALK .GE. 0.0D0 ) THEN |
---|
855 | W=1.0D0 |
---|
856 | ELSE |
---|
857 | W=-1.0D0 |
---|
858 | ENDIF |
---|
859 | IF(W.GE.0.0D0) THEN |
---|
860 | X1=X |
---|
861 | Y1=Y |
---|
862 | ELSE |
---|
863 | X2=X |
---|
864 | Y2=Y |
---|
865 | ENDIF |
---|
866 | LQ=KARL |
---|
867 | IF(LQ.EQ.1) THEN |
---|
868 | KARL=2*NINT(W) |
---|
869 | ELSEIF(IABS(LQ).EQ.2.AND.(LQ*W).LT.0.) THEN |
---|
870 | FAK=0.5D0 |
---|
871 | KARL=3 |
---|
872 | ENDIF |
---|
873 | IF(KARL.EQ.3.AND.STEG.LT.STEGBY) THEN |
---|
874 | W=(X2-X1)/(Y2-Y1) |
---|
875 | X=X1+W*(ALK-Y1) |
---|
876 | ELSE |
---|
877 | STEG=STEG*FAK |
---|
878 | X=X+STEG*W |
---|
879 | ENDIF |
---|
880 | AHPLUS=EXP(X) |
---|
881 | ENDIF |
---|
882 | IF(.NOT.DONE) GOTO 10 |
---|
883 | |
---|
884 | HCO3=CO3*AHPLUS/AK2C |
---|
885 | IF(ICALC.EQ.4) THEN |
---|
886 | CTOT=H2CO3+HCO3+CO3 |
---|
887 | ELSEIF(ICALC.EQ.1) THEN |
---|
888 | H2CO3=HCO3*AHPLUS/AK1C |
---|
889 | PCO2=H2CO3/AKP |
---|
890 | ENDIF |
---|
891 | PH=-LOG10(AHPLUS) |
---|
892 | ALKC=ALK-ALKB |
---|
893 | ELSEIF(ICALC.EQ.2) THEN |
---|
894 | ! *** CTOT, PCO2, AND BTOT FIXED *** |
---|
895 | Y=SQRT(PROD*(PROD-4.0D0*AKP*PCO2+4.0D0*CTOT)) |
---|
896 | H2CO3=PCO2*AKP |
---|
897 | HCO3=(Y-PROD)/2.0D0 |
---|
898 | CO3=CTOT-H2CO3-HCO3 |
---|
899 | ALKC=HCO3+2.0D0*CO3 |
---|
900 | AHPLUS=AK1C*H2CO3/HCO3 |
---|
901 | PH=-LOG10(AHPLUS) |
---|
902 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
903 | ALK=ALKC+ALKB |
---|
904 | ELSEIF(ICALC.EQ.3) THEN |
---|
905 | ! *** CTOT, PH AND BTOT FIXED *** |
---|
906 | FACTOR=CTOT/(AHPLUS*AHPLUS+AK1C*AHPLUS+AK1C*AK2C) |
---|
907 | CO3=FACTOR*AK1C*AK2C |
---|
908 | HCO3=FACTOR*AK1C*AHPLUS |
---|
909 | H2CO3=FACTOR*AHPLUS*AHPLUS |
---|
910 | PCO2=H2CO3/AKP |
---|
911 | ALKC=HCO3+2.0D0*CO3 |
---|
912 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
913 | ALK=ALKC+ALKB |
---|
914 | ELSEIF(ICALC.EQ.5) THEN |
---|
915 | ! *** ALK, PH AND BTOT FIXED *** |
---|
916 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
917 | ALKC=ALK-ALKB |
---|
918 | HCO3=ALKC/(1.0D0+2.0D0*AK2C/AHPLUS) |
---|
919 | CO3=HCO3*AK2C/AHPLUS |
---|
920 | H2CO3=HCO3*AHPLUS/AK1C |
---|
921 | PCO2=H2CO3/AKP |
---|
922 | CTOT=H2CO3+HCO3+CO3 |
---|
923 | ELSEIF(ICALC.EQ.6) THEN |
---|
924 | ! *** PCO2, PH AND BTOT FIXED *** |
---|
925 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
926 | H2CO3=PCO2*AKP |
---|
927 | HCO3=H2CO3*AK1C/AHPLUS |
---|
928 | CO3=HCO3*AK2C/AHPLUS |
---|
929 | CTOT=H2CO3+HCO3+CO3 |
---|
930 | ALKC=HCO3+2.0D0*CO3 |
---|
931 | ALK=ALKC+ALKB |
---|
932 | ENDIF |
---|
933 | ELSE |
---|
934 | IF(ICALC.EQ.1) THEN |
---|
935 | ! *** CTOT AND ALK FIXED *** |
---|
936 | TERM=4.0D0*ALK+CTOT*AKR-ALK*AKR |
---|
937 | Z=SQRT(TERM*TERM+4.0D0*(AKR-4.0D0)*ALK*ALK) |
---|
938 | CO3=(ALK*AKR-CTOT*AKR-4.0D0*ALK+Z)/(2.0D0*(AKR-4.0D0)) |
---|
939 | HCO3=(CTOT*AKR-Z)/(AKR-4.0D0) |
---|
940 | H2CO3=CTOT-ALK+CO3 |
---|
941 | PCO2=H2CO3/AKP |
---|
942 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
943 | ELSEIF(ICALC.EQ.2) THEN |
---|
944 | ! *** CTOT AND PCO2 FIXED *** |
---|
945 | Y=SQRT(PROD*(PROD-4.0D0*AKP*PCO2+4.0D0*CTOT)) |
---|
946 | H2CO3=PCO2*AKP |
---|
947 | HCO3=(Y-PROD)/2.0D0 |
---|
948 | CO3=CTOT-H2CO3-HCO3 |
---|
949 | ALK=HCO3+2.0D0*CO3 |
---|
950 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
951 | ELSEIF(ICALC.EQ.3) THEN |
---|
952 | ! *** CTOT AND PH FIXED *** |
---|
953 | FACTOR=CTOT/(AHPLUS*AHPLUS+AK1C*AHPLUS+AK1C*AK2C) |
---|
954 | CO3=FACTOR*AK1C*AK2C |
---|
955 | HCO3=FACTOR*AK1C*AHPLUS |
---|
956 | H2CO3=FACTOR*AHPLUS*AHPLUS |
---|
957 | PCO2=H2CO3/AKP |
---|
958 | ALK=HCO3+2.0D0*CO3 |
---|
959 | ELSEIF(ICALC.EQ.4) THEN |
---|
960 | ! *** ALK AND PCO2 FIXED *** |
---|
961 | TERM=SQRT((8.0D0*ALK+PROD)*PROD) |
---|
962 | CO3=ALK/2.0D0+PROD/8.0D0-TERM/8.0D0 |
---|
963 | HCO3=-PROD/4.0D0+TERM/4.0D0 |
---|
964 | H2CO3=PCO2*AKP |
---|
965 | CTOT=CO3+HCO3+H2CO3 |
---|
966 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
967 | ELSEIF(ICALC.EQ.5) THEN |
---|
968 | ! *** ALK AND PH FIXED *** |
---|
969 | HCO3=ALK/(1.0D0+2.0D0*AK2C/AHPLUS) |
---|
970 | CO3=HCO3*AK2C/AHPLUS |
---|
971 | H2CO3=HCO3*AHPLUS/AK1C |
---|
972 | PCO2=H2CO3/AKP |
---|
973 | CTOT=H2CO3+HCO3+CO3 |
---|
974 | ELSEIF(ICALC.EQ.6) THEN |
---|
975 | ! *** PCO2 AND PH FIXED *** |
---|
976 | H2CO3=PCO2*AKP |
---|
977 | HCO3=H2CO3*AK1C/AHPLUS |
---|
978 | CO3=HCO3*AK2C/AHPLUS |
---|
979 | CTOT=H2CO3+HCO3+CO3 |
---|
980 | ALK=HCO3+2.0D0*CO3 |
---|
981 | ENDIF |
---|
982 | ENDIF |
---|
983 | |
---|
984 | DO 120 II=1,NCONC |
---|
985 | CONCS(II)=CONCS2(II) |
---|
986 | 120 CONTINUE |
---|
987 | |
---|
988 | ! AXY (07/05/13) ================================================== |
---|
989 | ! C_SW IS SET AT 0 TO START |
---|
990 | ! THEN IF NON CONVERGENCE C_SW SET TO 1 |
---|
991 | 123 IF(C_SW.EQ.1)THEN |
---|
992 | ! IF NON CONVERGENCE, THE MODEL REQUIRES CONCS TO CONTAIN USABLE VALUES. |
---|
993 | ! BEST OFFER BEING THE OLD CONCS VALUES WHEN CONVERGENCE HAS BEEN |
---|
994 | ! ACHIEVED |
---|
995 | DO II=1,NCONC |
---|
996 | CONCS(II)=OLD_CONCS(II) |
---|
997 | END DO |
---|
998 | |
---|
999 | ! SPECIFIC CARBONATE VALUES TO PUSH CODE ON THROUGH THE |
---|
1000 | ! NON CONVERGENCE CONDITIONS |
---|
1001 | ! PCO2W = 0 SO CO2 UPTAKE WILL BE ENCOURAGED |
---|
1002 | ! CONCS(3)=O3C(III)*0.005_fp8/1.e6_fp8 |
---|
1003 | ! -99 IS A FLAG TO SHOW WHEN AND WHERE NON CONVERGENCE OCCURED |
---|
1004 | ! CONCS(4)=PH, OUTPUT OF -99 MEANS NON-CONVERGENCE |
---|
1005 | ! CONCS(4)=-99._fp8 |
---|
1006 | ! |
---|
1007 | ! AXY (10/05/13): remove this -99 flag; this is handled instead up in |
---|
1008 | ! trcbio_medusa.F90 where the iters variable is both |
---|
1009 | ! output and may be used to trigger action |
---|
1010 | ! |
---|
1011 | ! RESET SWITCH FOR NEXT CALL TO THIS SUBROUTINE |
---|
1012 | C_SW=0 |
---|
1013 | ENDIF |
---|
1014 | ! AXY (07/05/13) ================================================== |
---|
1015 | |
---|
1016 | |
---|
1017 | RETURN |
---|
1018 | |
---|
1019 | END SUBROUTINE |
---|
1020 | !----------------------------------------------------------------------- |
---|
1021 | |
---|
1022 | !======================================================================= |
---|
1023 | !======================================================================= |
---|
1024 | !======================================================================= |
---|
1025 | |
---|
1026 | !======================================================================= |
---|
1027 | ! |
---|
1028 | SUBROUTINE CaCO3_Saturation (Tc, S, D, CO3, & |
---|
1029 | Om_cal, Om_arg) |
---|
1030 | ! |
---|
1031 | !======================================================================= |
---|
1032 | ! |
---|
1033 | ! Routine to calculate the saturation state of calcite and aragonite |
---|
1034 | ! Inputs: |
---|
1035 | ! Tc Temperature (C) |
---|
1036 | ! S Salinity |
---|
1037 | ! D Depth (m) |
---|
1038 | ! CO3 Carbonate ion concentration (mol.kg-1 ie /1D6) |
---|
1039 | ! |
---|
1040 | ! Outputs |
---|
1041 | ! Om_cal Calite saturation |
---|
1042 | ! Om_arg Aragonite saturation |
---|
1043 | ! |
---|
1044 | ! Intermediates |
---|
1045 | ! K_cal Stoichiometric solubility product for calcite |
---|
1046 | ! K_arg Stoichiometric solubility product for aragonite |
---|
1047 | ! Ca Calcium 2+ concentration (mol.kg-1) |
---|
1048 | ! P Pressure (bars) |
---|
1049 | ! |
---|
1050 | ! Source |
---|
1051 | ! Zeebe & Wolf-Gladrow 2001 following Mucci (1983) |
---|
1052 | ! with pressure corrections from Millero (1995) |
---|
1053 | ! Code tested against reference values given in Z & W-G |
---|
1054 | ! Built Jerry Blackford, 2008 |
---|
1055 | !======================================================================= |
---|
1056 | |
---|
1057 | IMPLICIT None |
---|
1058 | REAL(wp), INTENT( in ) :: Tc, S, D, CO3 |
---|
1059 | REAL(wp), INTENT( inout ) :: Om_cal, Om_arg |
---|
1060 | REAL(wp) :: Tk, Kelvin, Ca |
---|
1061 | REAL(wp) :: logKspc, Kspc |
---|
1062 | REAL(wp) :: logKspa, Kspa |
---|
1063 | REAL(wp) :: tmp1, tmp2, tmp3 |
---|
1064 | REAL(wp) :: dV, dK, P, R |
---|
1065 | |
---|
1066 | ! setup |
---|
1067 | Kelvin = 273.15 |
---|
1068 | Tk = Tc + Kelvin |
---|
1069 | Ca = 0.01028 ! Currently oceanic mean value at S=25, needs refining) |
---|
1070 | Ca = 0.010279 * (S / 35.0) ! Ca varies with salinity (cf. Feeley et al., 2004; Yool et al., 2010) |
---|
1071 | R = 83.131 !(cm3.bar.mol-1.K-1) |
---|
1072 | P = D / 10.0 !pressure in bars |
---|
1073 | |
---|
1074 | ! calculate K for calcite |
---|
1075 | tmp1 = -171.9065 - (0.077993*Tk) + (2839.319/Tk) + 71.595*log10(Tk) |
---|
1076 | tmp2 = + (-0.77712 + (0.0028426*Tk) + (178.34/Tk))*SQRT(S) |
---|
1077 | tmp3 = - (0.07711*S) + (0.0041249*(S**1.5)) |
---|
1078 | logKspc = tmp1 + tmp2 + tmp3 |
---|
1079 | Kspc = 10.0**logKspc |
---|
1080 | |
---|
1081 | ! correction for pressure for calcite |
---|
1082 | IF ( D .GT. 0) THEN |
---|
1083 | dV = -48.76 + 0.5304*Tc |
---|
1084 | dK = -11.76/1.0D3 + (0.3692/1.0D3) * Tc |
---|
1085 | tmp1 = -(dV/(R*Tk))*P + (0.5*dK/(R*Tk))*P*P |
---|
1086 | Kspc = Kspc*exp(tmp1) |
---|
1087 | logKspc = log10(Kspc) |
---|
1088 | END IF |
---|
1089 | |
---|
1090 | ! calculate K for aragonite |
---|
1091 | tmp1 = -171.945 - 0.077993*Tk + 2903.293 / Tk + 71.595* log10(Tk) |
---|
1092 | tmp2 = + (-0.068393 + 0.0017276*Tk + 88.135/Tk)*SQRT(S) |
---|
1093 | tmp3 = - 0.10018*S + 0.0059415*S**1.5 |
---|
1094 | logKspa = tmp1 + tmp2 + tmp3 |
---|
1095 | Kspa = 10.0**logKspa |
---|
1096 | |
---|
1097 | ! correction for pressure for aragonite |
---|
1098 | IF ( D .GT. 0) THEN |
---|
1099 | dV = -46.00 + 0.5304*Tc |
---|
1100 | dK = -11.76/1.0D3 + (0.3692/1.0D3) * Tc |
---|
1101 | tmp1 = -(dV/(R*Tk))*P + (0.5*dK/(R*Tk))*P*P |
---|
1102 | Kspa = Kspa*exp(tmp1) |
---|
1103 | logKspa = log10(Kspa) |
---|
1104 | END IF |
---|
1105 | |
---|
1106 | ! calculate saturation states |
---|
1107 | Om_cal = (CO3 * Ca) / Kspc |
---|
1108 | Om_arg = (CO3 * Ca) / Kspa |
---|
1109 | |
---|
1110 | RETURN |
---|
1111 | |
---|
1112 | END SUBROUTINE |
---|
1113 | !----------------------------------------------------------------------- |
---|
1114 | |
---|
1115 | !======================================================================= |
---|
1116 | !======================================================================= |
---|
1117 | !======================================================================= |
---|
1118 | |
---|
1119 | # else |
---|
1120 | !!====================================================================== |
---|
1121 | !! Dummy module : No MEDUSA carbonate chemistry |
---|
1122 | !!====================================================================== |
---|
1123 | |
---|
1124 | CONTAINS |
---|
1125 | |
---|
1126 | SUBROUTINE trc_co2_medusa( kt ) ! Empty routine |
---|
1127 | |
---|
1128 | INTEGER, INTENT( in ) :: kt |
---|
1129 | |
---|
1130 | WRITE(*,*) 'trc_co2_medusa: You should not have seen this print! error?', kt |
---|
1131 | |
---|
1132 | END SUBROUTINE trc_co2_medusa |
---|
1133 | #endif |
---|
1134 | |
---|
1135 | !!====================================================================== |
---|
1136 | END MODULE trcco2_medusa |
---|