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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trcco2_medusa.F90 in branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA – NEMO

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcco2_medusa.F90

Last change on this file was 5841, checked in by jpalmier, 8 years ago

JPALM --30-10-2015-- Add MOCSY and DMS to MEDUSA-NEMO3.6

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