source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcco2_medusa.F90 @ 5707

Last change on this file since 5707 was 5707, checked in by acc, 5 years ago

JPALM —25-08-2015 — add MEDUSA in the branch. MEDUSA version already up-to-date with this trunk revision

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