Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

trcco2_medusa.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 50.0 KB

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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcco2_medusa.F90 @ 11738

Download in other formats: