[5726] | 1 | MODULE trcco2_medusa |
---|
| 2 | !!====================================================================== |
---|
| 3 | !! *** MODULE trcco2_medusa *** |
---|
| 4 | !! TOP : MEDUSA |
---|
| 5 | !!====================================================================== |
---|
| 6 | !! History : |
---|
| 7 | !! - ! 2010-12 (A. Yool) added for ROAM project |
---|
| 8 | !!---------------------------------------------------------------------- |
---|
| 9 | #if defined key_medusa && defined key_roam |
---|
| 10 | !!---------------------------------------------------------------------- |
---|
| 11 | !! MEDUSA carbonate chemistry |
---|
| 12 | !!---------------------------------------------------------------------- |
---|
| 13 | !! trc_co2_medusa : |
---|
| 14 | !!---------------------------------------------------------------------- |
---|
| 15 | USE oce_trc |
---|
| 16 | USE trc |
---|
| 17 | USE sms_medusa |
---|
| 18 | USE lbclnk |
---|
| 19 | USE prtctl_trc ! Print control for debugging |
---|
| 20 | USE in_out_manager ! I/O manager |
---|
| 21 | |
---|
| 22 | IMPLICIT NONE |
---|
| 23 | PRIVATE |
---|
| 24 | |
---|
| 25 | PUBLIC trc_co2_medusa ! called in trc_bio_medusa |
---|
| 26 | |
---|
| 27 | !!* Substitution |
---|
| 28 | # include "domzgr_substitute.h90" |
---|
| 29 | !!---------------------------------------------------------------------- |
---|
| 30 | !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) |
---|
| 31 | !! $Id$ |
---|
| 32 | !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) |
---|
| 33 | !!---------------------------------------------------------------------- |
---|
| 34 | |
---|
| 35 | ! The following is a map of the subroutines contained within this module |
---|
| 36 | ! - trc_co2_medusa |
---|
| 37 | ! - CALLS CO2_dynamics |
---|
| 38 | ! - CALLS CO2DYN |
---|
| 39 | ! - CALLS POLYCO |
---|
| 40 | ! - CALLS CO2SET |
---|
| 41 | ! - CALLS CO2CLC |
---|
| 42 | ! - CALLS CaCO3_Saturation |
---|
| 43 | ! - CALLS Air_sea_exchange |
---|
| 44 | |
---|
| 45 | CONTAINS |
---|
| 46 | |
---|
| 47 | !======================================================================= |
---|
| 48 | ! |
---|
[5841] | 49 | SUBROUTINE trc_co2_medusa( Temp, Sal, DIC, ALK, Depth, xkw, pCO2a, & |
---|
[5726] | 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 |
---|
[5841] | 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 |
---|
[5726] | 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 |
---|
[5841] | 89 | ! REAL(wp), INTENT( in ) :: Wnd ! m / s |
---|
| 90 | REAL(wp), INTENT( in ) :: xkw ! m / s |
---|
[5726] | 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 |
---|
[5841] | 100 | REAL(wp), INTENT( inout ) :: co2flux ! mmol / m2 / s |
---|
[5726] | 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 |
---|
[5841] | 134 | call Air_sea_exchange( Temp, xkw, pCO2w, pCO2a, henry, dcf, & ! inputs |
---|
[5726] | 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 |
---|
[5841] | 150 | IF(lwp) WRITE(numout,*) ' trc_co2_medusa: f_kw660 =', xkw |
---|
[5726] | 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', & |
---|
[5841] | 167 | & ' XKW', xkw, ' PH ', ph |
---|
[5726] | 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 |
---|
[5841] | 201 | ! WRITE(*,'(A27,F10.3)') " air sea flux(mmol/m2/s) = ", flux |
---|
[5726] | 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 | ! |
---|
[5841] | 292 | SUBROUTINE Air_sea_exchange( T, xkw, pco2w, pco2a, henry, dcf, & |
---|
[5726] | 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) |
---|
[5841] | 307 | ! Wnd wind speed, metres (DELETED) |
---|
| 308 | ! xkw gas transfer velocity |
---|
[5726] | 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 | |
---|
[5841] | 318 | REAL(wp), INTENT( in ) :: T, xkw, pco2w, pco2a, henry, dcf ! INPUT PARAMETERS: |
---|
[5726] | 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 |
---|
[5841] | 326 | ! fwind = (0.222d0 * wnd**2d0 + 0.333d0 * wnd)*(sc/660.d0)**(-0.5) |
---|
| 327 | fwind = xkw * (sc/660.d0)**(-0.5) |
---|
[5726] | 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 | |
---|
[5841] | 334 | ! AXY (23/06/15): let's get it from /d to /s |
---|
| 335 | flux = flux / ( 86400. ) |
---|
| 336 | |
---|
[5726] | 337 | RETURN |
---|
| 338 | |
---|
| 339 | END SUBROUTINE Air_sea_exchange |
---|
| 340 | !----------------------------------------------------------------------- |
---|
| 341 | |
---|
| 342 | !======================================================================= |
---|
| 343 | !======================================================================= |
---|
| 344 | !======================================================================= |
---|
| 345 | |
---|
| 346 | !======================================================================= |
---|
| 347 | ! |
---|
| 348 | SUBROUTINE CO2dyn ( TCO2, TA, T, S, PCO2A, PCO2, PH, HENRY, ca, bc, cb, iters ) |
---|
| 349 | ! |
---|
| 350 | !======================================================================= |
---|
| 351 | ! |
---|
| 352 | ! This subroutine acts as an interface to the Haltafall iteration, setting options etc. |
---|
| 353 | |
---|
| 354 | REAL(wp), INTENT( in ) :: TCO2, TA, T, S, PCO2A |
---|
| 355 | !----------------------------------------------------------------------- |
---|
| 356 | REAL(wp), INTENT( inout ) :: PCO2, PH, HENRY, ca, bc, cb |
---|
| 357 | INTEGER, INTENT( inout ) :: iters |
---|
| 358 | !----------------------------------------------------------------------- |
---|
| 359 | REAL(wp) :: PRSS |
---|
| 360 | INTEGER :: MCONC, MKVAL, ICONST, ICALC |
---|
| 361 | |
---|
| 362 | PARAMETER ( MCONC = 9,MKVAL = 4 ) |
---|
| 363 | |
---|
| 364 | REAL(wp), DIMENSION(MKVAL) :: AKVAL |
---|
| 365 | REAL(wp), DIMENSION(MCONC) :: CONCS |
---|
| 366 | |
---|
| 367 | ICONST = 6 |
---|
| 368 | PRSS = 1.0D0 |
---|
| 369 | CONCS(1) = TCO2 |
---|
| 370 | CONCS(2) = TA |
---|
| 371 | ICALC = 1 |
---|
| 372 | |
---|
| 373 | CALL POLYCO(PRSS,T,S,CONCS,MCONC,AKVAL,MKVAL,ICALC,ICONST,iters) |
---|
| 374 | |
---|
| 375 | if(iters.eq.25) then |
---|
| 376 | ! AXY (13/05/13): this location has failed to converge; use |
---|
| 377 | ! some standard values to plug the hole |
---|
| 378 | CONCS(3) = PCO2A * 1E-06 ! atmospheric value |
---|
| 379 | CONCS(4) = 8.1 ! global average pH |
---|
| 380 | CONCS(5) = TCO2 * 0.005 ! some "standard" values plucked |
---|
| 381 | CONCS(6) = TCO2 * 0.865 ! from pages 5-6 of Zeebe & Wolf- |
---|
| 382 | CONCS(7) = TCO2 * 0.130 ! Gladow (2001) |
---|
| 383 | AKVAL(1) = 0.5E-01 ! trial-and-error value |
---|
| 384 | IF(lwp) THEN |
---|
| 385 | WRITE(numout,*) ' CO2DYN : reset failed outputs' |
---|
| 386 | ENDIF |
---|
| 387 | endif |
---|
| 388 | |
---|
| 389 | PCO2 = CONCS(3) |
---|
| 390 | PH = CONCS(4) |
---|
| 391 | ca = CONCS(5) |
---|
| 392 | bc = CONCS(6) |
---|
| 393 | cb = CONCS(7) |
---|
| 394 | HENRY = AKVAL(1) |
---|
| 395 | |
---|
| 396 | RETURN |
---|
| 397 | |
---|
| 398 | END SUBROUTINE |
---|
| 399 | !----------------------------------------------------------------------- |
---|
| 400 | |
---|
| 401 | !======================================================================= |
---|
| 402 | !======================================================================= |
---|
| 403 | !======================================================================= |
---|
| 404 | |
---|
| 405 | !======================================================================= |
---|
| 406 | ! |
---|
| 407 | SUBROUTINE POLYCO(PD,TD,SD,CONCS,NCONC,AKVAL,NKVAL,ICALC,ICONST,iters) |
---|
| 408 | ! |
---|
| 409 | !======================================================================= |
---|
| 410 | ! |
---|
| 411 | ! MASTER SUBROUTINE FOR CALCULATION OF THE CO2 SYSTEM THERMODYNAMICS |
---|
| 412 | ! |
---|
| 413 | ! EXPLANATION OF POLYCO PARAMETERS |
---|
| 414 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
| 415 | ! P - PRESSURE IN ATMOSPHERES (P<>1 CODED only for ICONST=3 and ICONST=6) |
---|
| 416 | ! T - TEMPERATURE IN DEG.C |
---|
| 417 | ! S - SALINITY IN PPT |
---|
| 418 | ! CONCS(1) - TOTAL C (MOL/KG) |
---|
| 419 | ! CONCS(2) - TOTAL ALKALINITY (MOL/KG) |
---|
| 420 | ! CONCS(3) - PCO2 (ATM) |
---|
| 421 | ! CONCS(4) - PH |
---|
| 422 | ! CONCS(5) - {H2CO3} (MOL/KG) |
---|
| 423 | ! CONCS(6) - {HCO3} (MOL/KG) |
---|
| 424 | ! CONCS(7) - {CO3} (MOL/KG) |
---|
| 425 | ! CONCS(8) - CARBONATE ALKALINITY ) FOR ICONST = 4,5,6 |
---|
| 426 | ! CONCS(9) - BORATE ALKALINITY ) ONLY |
---|
| 427 | ! NCONC - SIZE OF CONCS ARRAY (7 FOR ICONST=1,2,3; 9 FOR ICONST |
---|
| 428 | ! AKVAL(1) - KP (HENRY'S LAW CONSTANT) (MOL/KG/ATM) |
---|
| 429 | ! AKVAL(2) - K1C (H2CO3 DISSOCIATION) (MOL/KG) |
---|
| 430 | ! AKVAL(3) - K2C (HCO3 DISSOCIATION) (MOL/KG) |
---|
| 431 | ! AKVAL(4) - KB (B(OH)3 DISSOCIATION) (MOL/KG) FOR ICONST=4,5,6 |
---|
| 432 | ! NKVAL - SIZE OF AKVAL ARRAY (3 FOR ICONST=1,2,3; 4 FOR ICONST= |
---|
| 433 | ! ICALC - SELECTION OF THE TWO INPUT PARAMETERS: |
---|
| 434 | ! ICALC = 1 TOTAL C AND ALKALINITY |
---|
| 435 | ! ICALC = 2 TOTAL C AND PCO2 |
---|
| 436 | ! ICALC = 3 TOTAL C AND PH |
---|
| 437 | ! ICALC = 4 ALKALINITY AND PCO2 |
---|
| 438 | ! ICALC = 5 ALKALINITY AND PH |
---|
| 439 | ! ICALC = 6 PCO2 AND PH |
---|
| 440 | ! ICALC = 7 CALCULATE CONSTANTS AKVAL ONLY |
---|
| 441 | ! ICONST - SELECTION OF PH SCALE AND COMPONENTS: |
---|
| 442 | ! ICONST = 1 NBS PH SCALE |
---|
| 443 | ! ICONST = 2 HANSSON'S SCALE (SWS WITHOUT FLUORIDE) |
---|
| 444 | ! ICONST = 3 SWS PH SCALE |
---|
| 445 | ! ICONST = 4 AS 1 BUT INCLUDING BORATE IN THE CALCULATION |
---|
| 446 | ! ICONST = 5 AS 2 BUT INCLUDING BORATE IN THE CALCULATION |
---|
| 447 | ! ICONST = 6 AS 3 BUT INCLUDING BORATE IN THE CALCULATION |
---|
| 448 | |
---|
| 449 | ! NOTE: FOR ICONST=1,2,3 CONCS(2) REPRESENTS CARBONATE ALKALINITY SINC |
---|
| 450 | ! BORATE IS NOT INCLUDED IN THE CALCULATION. FOR ICONST=4,5,6 CO |
---|
| 451 | ! REPRESENTS TOTAL ALKALINITY (CARBONATE + BORATE), THE COMPONEN |
---|
| 452 | ! WHICH ARE GIVEN IN CONCS(8) AND CONCS(9) |
---|
| 453 | |
---|
| 454 | REAL(wp) :: PMIN, PMAX, SMIN, SMAX, TMIN, TMAX, & |
---|
| 455 | & PD, TD, SD, P, T, S, BTOT |
---|
| 456 | INTEGER :: MINJC, MAXJC, MINJK, MAXJK, MINCAL, MAXCAL, MINCON, & |
---|
| 457 | & MAXCON, NCONC, NKVAL, ICALC, ICONST, IC, iters |
---|
| 458 | LOGICAL :: BORON |
---|
| 459 | |
---|
| 460 | PARAMETER(MINJC=7,MAXJC=9,MINJK=3,MAXJK=4) |
---|
| 461 | PARAMETER(MINCAL=1,MAXCAL=7,MINCON=1,MAXCON=6) |
---|
| 462 | !! AXY (11/08/11): TMIN changed to -2 to stop error messages in polar regions |
---|
| 463 | !! AXY (09/01/14): TMIN changed to -3 to stop error messages in polar regions |
---|
| 464 | !! AXY (03/03/14): SMAX changed to 42 to stop error messages in salty regions (the Great Answer ...) |
---|
| 465 | !! AXY (05/03/14): SMAX changed to 45 to stop error messages in salty regions |
---|
| 466 | PARAMETER(PMIN=0.99999D0,PMAX=1.00001D0,SMIN=0.0D0, & |
---|
| 467 | & SMAX=45.0D0,TMIN=-3.0D0,TMAX=40.0D0) |
---|
| 468 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
| 469 | REAL(wp), DIMENSION(NCONC) :: CONCS |
---|
| 470 | |
---|
| 471 | P = PD |
---|
| 472 | S = SD |
---|
| 473 | T = TD |
---|
| 474 | |
---|
| 475 | !! AXY (09/08/11): STOP statements commented out and WARNING messages added |
---|
| 476 | IF(lwp) THEN |
---|
| 477 | IF(T.LT.TMIN.OR.T.GT.TMAX) WRITE(numout,*) ' trc_co2_medusa: T WARNING, ', T, TMIN, TMAX, S, P |
---|
| 478 | IF(S.LT.SMIN.OR.S.GT.SMAX) WRITE(numout,*) ' trc_co2_medusa: S WARNING, ', S, SMIN, SMAX, T, P |
---|
| 479 | IF(P.LT.PMIN.OR.P.GT.PMAX) WRITE(numout,*) ' trc_co2_medusa: P WARNING, ', P, PMIN, PMAX, T, S |
---|
| 480 | ENDIF |
---|
| 481 | ! IF(T.LT.TMIN.OR.T.GT.TMAX)WRITE (*,*) P, S, T, TMIN, TMAX |
---|
| 482 | ! IF(P.LT.PMIN.OR.P.GT.PMAX) STOP'POLYCO - PRESSURE OUT OF RANGE' |
---|
| 483 | ! IF(S.LT.SMIN.OR.S.GT.SMAX) STOP'POLYCO - SALINITY OUT OF RANGE' |
---|
| 484 | ! IF(T.LT.TMIN.OR.T.GT.TMAX) STOP'POLYCO - TEMP. OUT OF RANGE' |
---|
| 485 | |
---|
| 486 | !! AXY (17/04/13): iMARNET climate change simulation appears to be compromised |
---|
| 487 | !! by excessive Caspian Sea salinity (> 45 PSU); this may be |
---|
| 488 | !! a runoff glitch that doesn't manifest on local instance of |
---|
| 489 | !! NEMO; fudging this here by replacing excessively high |
---|
| 490 | !! salinity values with maximum salinity value for the purpose |
---|
| 491 | !! of carbonate chemistry calculations; if this works, a more |
---|
| 492 | !! sensible solution might be to kill Caspian, etc. |
---|
| 493 | IF(lwp) THEN |
---|
| 494 | IF(S.LT.SMIN) THEN |
---|
| 495 | WRITE(numout,*) ' trc_co2_medusa: S RESET, ', S, '->', SMIN |
---|
| 496 | S = SMIN |
---|
| 497 | ENDIF |
---|
| 498 | IF(S.GT.SMAX) THEN |
---|
| 499 | WRITE(numout,*) ' trc_co2_medusa: S RESET, ', S, '->', SMAX |
---|
| 500 | S = SMAX |
---|
| 501 | ENDIF |
---|
| 502 | ENDIF |
---|
| 503 | |
---|
| 504 | !! AXY (28/02/13): as above, but for T; this has been done for 1/4-degree |
---|
| 505 | !! simulations to provide a temporary fix to a strange |
---|
| 506 | !! glitch in ocean temperature; basically, T has been found |
---|
| 507 | !! to spike upwards arbitrarily (e.g. 24C -> 53C); this |
---|
| 508 | !! causes the carbonate chemistry calculations to fail; |
---|
| 509 | !! this fix aims to stop the problem so that we can estimate |
---|
| 510 | !! how frequent a problem this is; it's suggestive of a |
---|
| 511 | !! memory leak |
---|
| 512 | IF(lwp) THEN |
---|
| 513 | IF(T.LT.TMIN) THEN |
---|
| 514 | WRITE(numout,*) ' trc_co2_medusa: T RESET, ', T, '->', TMIN |
---|
| 515 | T = TMIN |
---|
| 516 | ENDIF |
---|
| 517 | IF(T.GT.TMAX) THEN |
---|
| 518 | WRITE(numout,*) ' trc_co2_medusa: T RESET, ', T, '->', TMAX |
---|
| 519 | T = TMAX |
---|
| 520 | ENDIF |
---|
| 521 | ENDIF |
---|
| 522 | |
---|
| 523 | IF(ICALC.LT.MINCAL.OR.ICALC.GT.MAXCAL) STOP 'POLYCO - ICALC OUT OR RANGE' |
---|
| 524 | IF(ICONST.LT.MINCON.OR.ICONST.GT.MAXCON) STOP 'POLYCO - ICONST OUT OF RANGE' |
---|
| 525 | BORON=(ICONST.GT.3) |
---|
| 526 | IF(BORON) THEN |
---|
| 527 | IC=ICONST-3 |
---|
| 528 | BTOT=0.0004128D0*S/35.0D0 |
---|
| 529 | IF(NCONC.NE.MAXJC) STOP 'POLYCO - WRONG NCONC VALUE' |
---|
| 530 | IF(NKVAL.NE.MAXJK) STOP 'POLYCO - WRONG NKVAL VALUE' |
---|
| 531 | ELSE |
---|
| 532 | IC=ICONST |
---|
| 533 | IF(NCONC.NE.MINJC) STOP 'POLYCO - WRONG NCONC VALUE' |
---|
| 534 | IF(NKVAL.NE.MINJK) STOP 'POLYCO - WRONG NKVAL VALUE' |
---|
| 535 | ENDIF |
---|
| 536 | |
---|
| 537 | CALL CO2SET(P,T,S,AKVAL,NKVAL,IC) |
---|
| 538 | |
---|
| 539 | IF(ICALC.LT.MAXCAL) & |
---|
| 540 | & CALL CO2CLC(CONCS,NCONC,AKVAL,NKVAL,ICALC,BORON,BTOT,iters) |
---|
| 541 | |
---|
| 542 | RETURN |
---|
| 543 | |
---|
| 544 | END SUBROUTINE |
---|
| 545 | !----------------------------------------------------------------------- |
---|
| 546 | |
---|
| 547 | !======================================================================= |
---|
| 548 | !======================================================================= |
---|
| 549 | !======================================================================= |
---|
| 550 | |
---|
| 551 | !======================================================================= |
---|
| 552 | ! |
---|
| 553 | SUBROUTINE CO2SET(P,T,S,AKVAL,NKVAL,IC) |
---|
| 554 | ! |
---|
| 555 | !======================================================================= |
---|
| 556 | ! |
---|
| 557 | ! Routine to calculate CO2 system constants under the conditions set by |
---|
| 558 | ! P,S,T (NOTE: PRESSURE <> 1ATM IS NOT YET CODED) |
---|
| 559 | |
---|
| 560 | ! I. Calculate constants at P=1 and S=0 using |
---|
| 561 | |
---|
| 562 | ! ln K0 = A + B/TK + C ln TK |
---|
| 563 | ! (where TK is in Kelvin) |
---|
| 564 | |
---|
| 565 | ! II. Calculate constants at P=1 and salinity S using |
---|
| 566 | |
---|
| 567 | ! ln K = ln K0 + (a0 + a1/TK + a2 ln TK) S**1/2 |
---|
| 568 | ! + (b0 + b1TK + b2TK**2) S |
---|
| 569 | |
---|
| 570 | ! The sources of the coefficients are as follows: |
---|
| 571 | |
---|
| 572 | ! IC= 1 2 3 |
---|
| 573 | ! (NBS pH scale) (SWS pH scale (SWS pH scale |
---|
| 574 | ! with no F) |
---|
| 575 | |
---|
| 576 | ! KP WEISS (1974) WEISS(1974) WEISS(1974 |
---|
| 577 | |
---|
| 578 | ! K1C ) MEHRBACH ACC. TO HANSSON ACC. TO HANSSON AND MEH |
---|
| 579 | ! K2C ) MILLERO (1979) MILLERO (1979) ACC. TO DICKS |
---|
| 580 | ! KB ) AND MILLERO (1 |
---|
| 581 | ! (K1C AND K2C |
---|
| 582 | ! HANSSON ACC |
---|
| 583 | ! MILLERO (1 |
---|
| 584 | ! (KB ONLY |
---|
| 585 | |
---|
| 586 | ! *** |
---|
| 587 | ! IMPLICIT real*8 (A-H,O-Z) |
---|
| 588 | |
---|
| 589 | ! Modified by jcb 17/02/10 to use OCMIP calculations of K1, K2, Kb. |
---|
| 590 | ! Differences are subtle rather than significant |
---|
| 591 | |
---|
| 592 | INTEGER :: MAXK, MAXCON, NKVAL, ICON, IC, IK |
---|
| 593 | ! *** |
---|
| 594 | PARAMETER(MAXK=4,MAXCON=3) |
---|
| 595 | REAL(wp), DIMENSION(MAXK) :: A, B, C |
---|
| 596 | REAL(wp), DIMENSION(MAXK,MAXCON) :: A0, A1, A2 |
---|
| 597 | REAL(wp), DIMENSION(MAXK,MAXCON) :: B0, B1, B2 |
---|
| 598 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
| 599 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
| 600 | ! REAL(wp) :: P,T,S,VAL,TK |
---|
| 601 | REAL(wp) :: P,T,S,VAL,TK, delta, kappa, Rgas |
---|
| 602 | REAL(wp) :: dlogTK, S2, sqrtS, S15, k1, k2, kb |
---|
| 603 | ! *** |
---|
| 604 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
| 605 | DATA Rgas/83.131/ |
---|
| 606 | |
---|
| 607 | DATA A/-167.8108D0, 290.9097D0, 207.6548D0, 148.0248D0/ |
---|
| 608 | DATA B/9345.17D0, -14554.21D0, -11843.79D0, -8966.9D0/ |
---|
| 609 | DATA C/23.3585D0, -45.0575D0, -33.6485D0, -24.4344D0/ |
---|
| 610 | DATA (A0(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 611 | DATA (A0(2,ICON),ICON=1,MAXCON) /0.0221D0, 0.5709D0, -45.8076D0/ |
---|
| 612 | DATA (A0(3,ICON),ICON=1,MAXCON) /0.9805D0, 1.4853D0, -39.5492D0/ |
---|
| 613 | DATA (A0(4,ICON),ICON=1,MAXCON) /0.0473D0, 0.5998D0, 0.5998D0/ |
---|
| 614 | DATA (A1(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 615 | DATA (A1(2,ICON),ICON=1,MAXCON) /34.02D0, -84.25D0, 1935.07D0/ |
---|
| 616 | DATA (A1(3,ICON),ICON=1,MAXCON) /-92.65D0, -192.69D0, 1590.14D0/ |
---|
| 617 | DATA (A1(4,ICON),ICON=1,MAXCON) /49.10D0, -75.25D0, -75.25D0/ |
---|
| 618 | DATA (A2(1,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 619 | DATA (A2(2,ICON),ICON=1,MAXCON) /2*0.0D0,6.9513D0/ |
---|
| 620 | DATA (A2(3,ICON),ICON=1,MAXCON) /2*0.0D0,6.1523D0/ |
---|
| 621 | DATA (A2(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 622 | DATA (B0(1,ICON),ICON=1,MAXCON) /3*0.023517D0/ |
---|
| 623 | DATA (B0(2,ICON),ICON=1,MAXCON) /0.0D0,-0.01632D0,-0.01566D0/ |
---|
| 624 | DATA (B0(3,ICON),ICON=1,MAXCON) /-0.03294D0,-0.05058D0,-0.04997D0/ |
---|
| 625 | DATA (B0(4,ICON),ICON=1,MAXCON) /0.0D0, -0.01767D0, -0.01767D0/ |
---|
| 626 | DATA (B1(1,ICON),ICON=1,MAXCON) /3*-2.3656D-4/ |
---|
| 627 | DATA (B1(2,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 628 | DATA (B1(3,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 629 | DATA (B1(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 630 | DATA (B2(1,ICON),ICON=1,MAXCON) /3*4.7036D-7/ |
---|
| 631 | DATA (B2(2,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 632 | DATA (B2(3,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 633 | DATA (B2(4,ICON),ICON=1,MAXCON) /3*0.0D0/ |
---|
| 634 | |
---|
| 635 | TK=T+273.15D0 |
---|
| 636 | DO 100 IK=1,NKVAL |
---|
| 637 | VAL=A(IK) + B(IK)/TK + C(IK)*LOG(TK) |
---|
| 638 | VAL=VAL + (A0(IK,IC) + A1(IK,IC)/TK + A2(IK,IC)*LOG(TK))*SQRT(S) |
---|
| 639 | VAL=VAL + (B0(IK,IC) + B1(IK,IC)*TK + B2(IK,IC)*TK*TK)*S |
---|
| 640 | AKVAL(IK)=EXP(VAL) |
---|
| 641 | 100 CONTINUE |
---|
| 642 | |
---|
| 643 | IF (IC .EQ. 3) THEN |
---|
| 644 | ! Calculation of constants as used in the OCMIP process for ICONST = 3 or 6 |
---|
| 645 | ! see http://www.ipsl.jussieu.fr/OCMIP/ |
---|
| 646 | ! added jcb 17/02/10 |
---|
| 647 | |
---|
| 648 | ! Derive simple terms used more than once |
---|
| 649 | dlogTK = log(TK) |
---|
| 650 | S2 = S*S |
---|
| 651 | sqrtS = sqrt(S) |
---|
| 652 | S15 = S**1.5 |
---|
| 653 | ! k1 = [H][HCO3]/[H2CO3] |
---|
| 654 | ! k2 = [H][CO3]/[HCO3] |
---|
| 655 | ! Millero p.664 (1995) using Mehrbach et al. data on seawater scale |
---|
| 656 | k1=10**(-1*(3670.7/TK - 62.008 + 9.7944*dlogTK - & |
---|
| 657 | & 0.0118 * S + 0.000116*S2)) |
---|
| 658 | k2=10**(-1*(1394.7/TK + 4.777 - & |
---|
| 659 | & 0.0184*S + 0.000118*S2)) |
---|
| 660 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
| 661 | ! Correction for high pressure (from Millero 1995) |
---|
| 662 | ! added YA 04/10/2010 |
---|
| 663 | delta=-25.5+0.1271*T |
---|
| 664 | kappa=(-3.08+0.0877*T)/1000.0 |
---|
| 665 | k1=k1*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
| 666 | delta=-15.82-0.0219*T |
---|
| 667 | kappa=(1.13-0.1475*T)/1000.0 |
---|
| 668 | k2=k2*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
| 669 | |
---|
| 670 | ! kb = [H][BO2]/[HBO2] |
---|
| 671 | ! Millero p.669 (1995) using data from Dickson (1990) |
---|
| 672 | kb=exp((-8966.90 - 2890.53*sqrtS - 77.942*S + & |
---|
| 673 | & 1.728*S15 - 0.0996*S2)/TK + & |
---|
| 674 | & (148.0248 + 137.1942*sqrtS + 1.62142*S) + & |
---|
| 675 | & (-24.4344 - 25.085*sqrtS - 0.2474*S) * & |
---|
| 676 | & dlogTK + 0.053105*sqrtS*TK) |
---|
| 677 | ! AXY (16/08/11): Yuri change for non-surface carbonate chemistry |
---|
| 678 | ! Correction for high pressure (from Millero, 1995) |
---|
| 679 | ! added YA 04/10/2010 |
---|
| 680 | delta=-29.48+0.1622*T-0.002608*T**2.0 |
---|
| 681 | kappa=-2.84/1000.0 |
---|
| 682 | kb=kb*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
| 683 | ! Correction for high pressure of Kw (from Millero 1995) |
---|
| 684 | ! added YA 04/10/2010 |
---|
| 685 | delta=-25.60+0.2324*T-0.0036246*T**2 |
---|
| 686 | kappa=(-5.13+0.0794*T)/1000.0 |
---|
| 687 | AKVAL(1)=AKVAL(1)*exp((-delta+0.5*kappa*P)*P/(Rgas*TK)) |
---|
| 688 | |
---|
| 689 | ! replace haltafall calculations with OCMIP calculations |
---|
| 690 | AKVAL(2) = k1 |
---|
| 691 | AKVAL(3) = k2 |
---|
| 692 | AKVAL(4) = kb |
---|
| 693 | END IF ! section implimenting OCMIP coefficients |
---|
| 694 | |
---|
| 695 | RETURN |
---|
| 696 | END SUBROUTINE |
---|
| 697 | !----------------------------------------------------------------------- |
---|
| 698 | |
---|
| 699 | !======================================================================= |
---|
| 700 | !======================================================================= |
---|
| 701 | !======================================================================= |
---|
| 702 | |
---|
| 703 | !======================================================================= |
---|
| 704 | ! |
---|
| 705 | SUBROUTINE CO2CLC(CONCS,NCONC,AKVAL,NKVAL,ICALC,BORON,BTOT,iters) |
---|
| 706 | ! |
---|
| 707 | !======================================================================= |
---|
| 708 | ! |
---|
| 709 | ! ROUTINE TO CARRY OUT CO2 CALCULATIONS WITH 2 FIXED PARAMETERS ACCORDI |
---|
| 710 | ! THE EQUATIONS GIVEN BY PARKS(1969) AND SKIRROW (1975) |
---|
| 711 | ! WITH ADDITIONS FOR INCLUDING BORON IF BORON=.TRUE. |
---|
| 712 | |
---|
| 713 | |
---|
| 714 | ! IMPLICIT real*8 (A-H,O-Z) |
---|
| 715 | INTEGER :: NCONC, NKVAL, ICALC, II, KARL, LQ, iters |
---|
| 716 | ! AXY (07/05/13) ================================================== |
---|
| 717 | ! put counter in to check duration in convergence loop |
---|
| 718 | INTEGER :: COUNTER,C_CHECK,C_SW,III |
---|
| 719 | ! AXY (07/05/13) ================================================== |
---|
| 720 | REAL(wp) :: CTOT,ALK,PCO2,PH,H2CO3,HCO3,CO3,ALKC |
---|
| 721 | REAL(wp) :: ALKB,AKP,AK1C,AK2C,AKB,BTOT |
---|
| 722 | REAL(wp) :: AKR,AHPLUS |
---|
| 723 | REAL(wp) :: PROD,tol1,tol2,tol3,tol4,steg,fak |
---|
| 724 | REAL(wp) :: STEGBY,Y,X,W,X1,Y1,X2,Y2,FACTOR,TERM,Z |
---|
| 725 | REAL(wp), DIMENSION(NCONC) :: CONCS |
---|
| 726 | REAL(wp), DIMENSION(NKVAL) :: AKVAL |
---|
| 727 | REAL(wp), DIMENSION(9) :: CONCS2 |
---|
| 728 | REAL(wp), DIMENSION(4) :: AKVAL2 |
---|
| 729 | |
---|
| 730 | ! AXY (07/05/13) ================================================== |
---|
| 731 | ! setup array to store old values of concs |
---|
| 732 | real(wp),DIMENSION(:) :: old_CONCS(NCONC) |
---|
| 733 | ! AXY (07/05/13) ================================================== |
---|
| 734 | |
---|
| 735 | EQUIVALENCE (CTOT , CONCS2(1)), (ALK , CONCS2(2)), & |
---|
| 736 | & (PCO2 , CONCS2(3)), (PH , CONCS2(4)), & |
---|
| 737 | & (H2CO3 , CONCS2(5)), (HCO3 , CONCS2(6)), & |
---|
| 738 | & (CO3 , CONCS2(7)), (ALKC , CONCS2(8)), & |
---|
| 739 | & (ALKB , CONCS2(9)), & |
---|
| 740 | & (AKP , AKVAL2(1)), (AK1C , AKVAL2(2)), & |
---|
| 741 | & (AK2C , AKVAL2(3)), (AKB , AKVAL2(4)) |
---|
| 742 | LOGICAL :: BORON,DONE |
---|
| 743 | |
---|
| 744 | ! AXY (07/05/13) ================================================== |
---|
| 745 | ! DERIVING PH REQUIRES FOLLOWING LOOP TO CONVERGE. |
---|
| 746 | ! THIS SUBROUTINE RELIES ON CONVERGENCE. IF THE ENVIRONMENTAL |
---|
| 747 | ! CONDITIONS DO NOT ALLOW FOR CONVERGENCE (IN 3D MODEL THIS IS |
---|
| 748 | ! LIKELY TO OCCUR NEAR LOW SALINITY REGIONS) THE MODEL WILL |
---|
| 749 | ! BE STUCK IN THE LOOP. TO AVOID THIS A CONVERGENCE CONDITION |
---|
| 750 | ! IS PUT IN PLACE TO SET A FLAGG OF -99 IN THE PH VAR FOR NON CONVEGENCE. |
---|
| 751 | ! THE MODEL IS THEN ALLOWED TO CONTINUE. 'COUNTER, C_SW,C_CHECK' ARE |
---|
| 752 | ! THE LOCAL VARS USED. |
---|
| 753 | ! C_SW = condition of convergence 0=yes, 1= no |
---|
| 754 | ! COUNTER = number of iterations |
---|
| 755 | ! C_CHECK = maximum number of iterations |
---|
| 756 | |
---|
| 757 | ! SET COUNTER AND SWITCH TO ZERO AND OFF |
---|
| 758 | COUNTER=0 |
---|
| 759 | C_SW=0 |
---|
| 760 | ! FROM EXPERIENCE IF THE ITERATIONS IN THE FOLLOWING DO LOOP |
---|
| 761 | ! EXCEEDS 15 CONVERGENCE WILL NOT OCCUR. THE OVERHEAD OF 25 ITERATIONS |
---|
| 762 | ! IS OK FOR SMALL DOMAINS WITH 1/10 AND 1/15 DEG RESOLUTION. |
---|
| 763 | ! I RECOMMEND A LOWER VALUE OF 15 FOR HIGHER RESOLUTION OR LARGER DOMAINS. |
---|
| 764 | C_CHECK=25 |
---|
| 765 | ! AXY (07/05/13) ================================================== |
---|
| 766 | |
---|
| 767 | DO 100 II=1,NCONC |
---|
| 768 | CONCS2(II)=CONCS(II) |
---|
| 769 | |
---|
| 770 | ! AXY (07/05/13) ================================================== |
---|
| 771 | ! IF CONVERGENCE IS NOT ACHIEVED THE LOCAL ARRAY CONCS MUST BE STORED TO |
---|
| 772 | ! ALLOW THE MODEL TO CONTINUE. THEREFORE .... |
---|
| 773 | ! UPDATE OLD_CONCS |
---|
| 774 | old_CONCS(II)=CONCS(II) |
---|
| 775 | ! AXY (07/05/13) ================================================== |
---|
| 776 | |
---|
| 777 | 100 CONTINUE |
---|
| 778 | DO 110 II=1,NKVAL |
---|
| 779 | AKVAL2(II)=AKVAL(II) |
---|
| 780 | 110 CONTINUE |
---|
| 781 | AKR = AK1C/AK2C |
---|
| 782 | AHPLUS=10.0D0**(-PH) |
---|
| 783 | PROD=AKR*AKP*PCO2 |
---|
| 784 | |
---|
| 785 | IF(BORON) THEN |
---|
| 786 | |
---|
| 787 | IF(ICALC.EQ.1.OR.ICALC.EQ.4) THEN |
---|
| 788 | ! *** ALK, BTOT AND CTOT OR PCO2 FIXED *** |
---|
| 789 | ! *** ITERATIVE CALCULATION NECESSARY HERE |
---|
| 790 | |
---|
| 791 | ! SET INITIAL GUESSES AND TOLERANCE |
---|
| 792 | H2CO3=PCO2*AKP |
---|
| 793 | CO3=ALK/10.0D0 |
---|
| 794 | AHPLUS=1.0D-8 |
---|
| 795 | ALKB=BTOT |
---|
| 796 | TOL1=ALK/1.0D5 |
---|
| 797 | TOL2=H2CO3/1.0D5 |
---|
| 798 | TOL3=CTOT/1.0D5 |
---|
| 799 | TOL4=BTOT/1.0D5 |
---|
| 800 | |
---|
| 801 | ! HALTAFALL iteration to determine CO3, ALKB, AHPLUS |
---|
| 802 | KARL=1 |
---|
| 803 | STEG=2.0D0 |
---|
| 804 | FAK=1.0D0 |
---|
| 805 | STEGBY=0.4D0 |
---|
| 806 | 10 DONE=.TRUE. |
---|
| 807 | |
---|
| 808 | ! AXY (07/05/13) ================================================== |
---|
| 809 | ! SET COUNTER UPDATE. FLAG 10 IS THE POINT OF RETURN FOR |
---|
| 810 | ! THE CONVERGENCE CONDITION |
---|
| 811 | COUNTER=COUNTER+1 |
---|
| 812 | iters = COUNTER |
---|
| 813 | ! CHECK IF CONVERGENCE HAS OCCURED IN THE NUMBER OF |
---|
| 814 | ! ACCEPTABLE ITTERATIONS. SET C_SW TO LET MODEL KNOW |
---|
| 815 | ! WHAT TO DO AT THE END OF THE SUBROUTINE |
---|
| 816 | if(counter.ge.c_check)then |
---|
| 817 | IF(lwp) THEN |
---|
| 818 | WRITE(numout,*) ' CO2CLC : ITERS WARNING, ', iters |
---|
| 819 | ENDIF |
---|
| 820 | c_sw=1 |
---|
| 821 | GOTO 123 |
---|
| 822 | endif |
---|
| 823 | ! AXY (07/05/13) ================================================== |
---|
| 824 | |
---|
| 825 | IF(ICALC.EQ.4) THEN |
---|
| 826 | ! *** PCO2 IS FIXED *** |
---|
| 827 | Y=AHPLUS*AHPLUS*CO3/(AK1C*AK2C) |
---|
| 828 | IF(ABS(Y-H2CO3).GT.TOL2) THEN |
---|
| 829 | CO3=CO3*H2CO3/Y |
---|
| 830 | DONE=.FALSE. |
---|
| 831 | ENDIF |
---|
| 832 | ELSEIF(ICALC.EQ.1) THEN |
---|
| 833 | ! *** CTOT IS FIXED *** |
---|
| 834 | Y=CO3*(1.0D0+AHPLUS/AK2C+AHPLUS*AHPLUS/(AK1C*AK2C)) |
---|
| 835 | IF(ABS(Y-CTOT).GT.TOL3) THEN |
---|
| 836 | CO3=CO3*CTOT/Y |
---|
| 837 | DONE=.FALSE. |
---|
| 838 | ENDIF |
---|
| 839 | ENDIF |
---|
| 840 | Y=ALKB*(1.0D0+AHPLUS/AKB) |
---|
| 841 | IF(ABS(Y-BTOT).GT.TOL4) THEN |
---|
| 842 | ALKB=ALKB*BTOT/Y |
---|
| 843 | DONE=.FALSE. |
---|
| 844 | ENDIF |
---|
| 845 | |
---|
| 846 | ! Alkalinity is equivalent to -(total H+), so the sign of W is opposite |
---|
| 847 | ! to that normally used |
---|
| 848 | |
---|
| 849 | Y=CO3*(2.0D0+AHPLUS/AK2C)+ALKB |
---|
| 850 | IF(ABS(Y-ALK).GT.TOL1) THEN |
---|
| 851 | DONE=.FALSE. |
---|
| 852 | X=LOG(AHPLUS) |
---|
| 853 | ! W=SIGN(1.0D0,Y-ALK) |
---|
| 854 | IF ( Y-ALK .GE. 0.0D0 ) THEN |
---|
| 855 | W=1.0D0 |
---|
| 856 | ELSE |
---|
| 857 | W=-1.0D0 |
---|
| 858 | ENDIF |
---|
| 859 | IF(W.GE.0.0D0) THEN |
---|
| 860 | X1=X |
---|
| 861 | Y1=Y |
---|
| 862 | ELSE |
---|
| 863 | X2=X |
---|
| 864 | Y2=Y |
---|
| 865 | ENDIF |
---|
| 866 | LQ=KARL |
---|
| 867 | IF(LQ.EQ.1) THEN |
---|
| 868 | KARL=2*NINT(W) |
---|
| 869 | ELSEIF(IABS(LQ).EQ.2.AND.(LQ*W).LT.0.) THEN |
---|
| 870 | FAK=0.5D0 |
---|
| 871 | KARL=3 |
---|
| 872 | ENDIF |
---|
| 873 | IF(KARL.EQ.3.AND.STEG.LT.STEGBY) THEN |
---|
| 874 | W=(X2-X1)/(Y2-Y1) |
---|
| 875 | X=X1+W*(ALK-Y1) |
---|
| 876 | ELSE |
---|
| 877 | STEG=STEG*FAK |
---|
| 878 | X=X+STEG*W |
---|
| 879 | ENDIF |
---|
| 880 | AHPLUS=EXP(X) |
---|
| 881 | ENDIF |
---|
| 882 | IF(.NOT.DONE) GOTO 10 |
---|
| 883 | |
---|
| 884 | HCO3=CO3*AHPLUS/AK2C |
---|
| 885 | IF(ICALC.EQ.4) THEN |
---|
| 886 | CTOT=H2CO3+HCO3+CO3 |
---|
| 887 | ELSEIF(ICALC.EQ.1) THEN |
---|
| 888 | H2CO3=HCO3*AHPLUS/AK1C |
---|
| 889 | PCO2=H2CO3/AKP |
---|
| 890 | ENDIF |
---|
| 891 | PH=-LOG10(AHPLUS) |
---|
| 892 | ALKC=ALK-ALKB |
---|
| 893 | ELSEIF(ICALC.EQ.2) THEN |
---|
| 894 | ! *** CTOT, PCO2, AND BTOT FIXED *** |
---|
| 895 | Y=SQRT(PROD*(PROD-4.0D0*AKP*PCO2+4.0D0*CTOT)) |
---|
| 896 | H2CO3=PCO2*AKP |
---|
| 897 | HCO3=(Y-PROD)/2.0D0 |
---|
| 898 | CO3=CTOT-H2CO3-HCO3 |
---|
| 899 | ALKC=HCO3+2.0D0*CO3 |
---|
| 900 | AHPLUS=AK1C*H2CO3/HCO3 |
---|
| 901 | PH=-LOG10(AHPLUS) |
---|
| 902 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
| 903 | ALK=ALKC+ALKB |
---|
| 904 | ELSEIF(ICALC.EQ.3) THEN |
---|
| 905 | ! *** CTOT, PH AND BTOT FIXED *** |
---|
| 906 | FACTOR=CTOT/(AHPLUS*AHPLUS+AK1C*AHPLUS+AK1C*AK2C) |
---|
| 907 | CO3=FACTOR*AK1C*AK2C |
---|
| 908 | HCO3=FACTOR*AK1C*AHPLUS |
---|
| 909 | H2CO3=FACTOR*AHPLUS*AHPLUS |
---|
| 910 | PCO2=H2CO3/AKP |
---|
| 911 | ALKC=HCO3+2.0D0*CO3 |
---|
| 912 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
| 913 | ALK=ALKC+ALKB |
---|
| 914 | ELSEIF(ICALC.EQ.5) THEN |
---|
| 915 | ! *** ALK, PH AND BTOT FIXED *** |
---|
| 916 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
| 917 | ALKC=ALK-ALKB |
---|
| 918 | HCO3=ALKC/(1.0D0+2.0D0*AK2C/AHPLUS) |
---|
| 919 | CO3=HCO3*AK2C/AHPLUS |
---|
| 920 | H2CO3=HCO3*AHPLUS/AK1C |
---|
| 921 | PCO2=H2CO3/AKP |
---|
| 922 | CTOT=H2CO3+HCO3+CO3 |
---|
| 923 | ELSEIF(ICALC.EQ.6) THEN |
---|
| 924 | ! *** PCO2, PH AND BTOT FIXED *** |
---|
| 925 | ALKB=BTOT/(1.0D0+AHPLUS/AKB) |
---|
| 926 | H2CO3=PCO2*AKP |
---|
| 927 | HCO3=H2CO3*AK1C/AHPLUS |
---|
| 928 | CO3=HCO3*AK2C/AHPLUS |
---|
| 929 | CTOT=H2CO3+HCO3+CO3 |
---|
| 930 | ALKC=HCO3+2.0D0*CO3 |
---|
| 931 | ALK=ALKC+ALKB |
---|
| 932 | ENDIF |
---|
| 933 | ELSE |
---|
| 934 | IF(ICALC.EQ.1) THEN |
---|
| 935 | ! *** CTOT AND ALK FIXED *** |
---|
| 936 | TERM=4.0D0*ALK+CTOT*AKR-ALK*AKR |
---|
| 937 | Z=SQRT(TERM*TERM+4.0D0*(AKR-4.0D0)*ALK*ALK) |
---|
| 938 | CO3=(ALK*AKR-CTOT*AKR-4.0D0*ALK+Z)/(2.0D0*(AKR-4.0D0)) |
---|
| 939 | HCO3=(CTOT*AKR-Z)/(AKR-4.0D0) |
---|
| 940 | H2CO3=CTOT-ALK+CO3 |
---|
| 941 | PCO2=H2CO3/AKP |
---|
| 942 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
| 943 | ELSEIF(ICALC.EQ.2) THEN |
---|
| 944 | ! *** CTOT AND PCO2 FIXED *** |
---|
| 945 | Y=SQRT(PROD*(PROD-4.0D0*AKP*PCO2+4.0D0*CTOT)) |
---|
| 946 | H2CO3=PCO2*AKP |
---|
| 947 | HCO3=(Y-PROD)/2.0D0 |
---|
| 948 | CO3=CTOT-H2CO3-HCO3 |
---|
| 949 | ALK=HCO3+2.0D0*CO3 |
---|
| 950 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
| 951 | ELSEIF(ICALC.EQ.3) THEN |
---|
| 952 | ! *** CTOT AND PH FIXED *** |
---|
| 953 | FACTOR=CTOT/(AHPLUS*AHPLUS+AK1C*AHPLUS+AK1C*AK2C) |
---|
| 954 | CO3=FACTOR*AK1C*AK2C |
---|
| 955 | HCO3=FACTOR*AK1C*AHPLUS |
---|
| 956 | H2CO3=FACTOR*AHPLUS*AHPLUS |
---|
| 957 | PCO2=H2CO3/AKP |
---|
| 958 | ALK=HCO3+2.0D0*CO3 |
---|
| 959 | ELSEIF(ICALC.EQ.4) THEN |
---|
| 960 | ! *** ALK AND PCO2 FIXED *** |
---|
| 961 | TERM=SQRT((8.0D0*ALK+PROD)*PROD) |
---|
| 962 | CO3=ALK/2.0D0+PROD/8.0D0-TERM/8.0D0 |
---|
| 963 | HCO3=-PROD/4.0D0+TERM/4.0D0 |
---|
| 964 | H2CO3=PCO2*AKP |
---|
| 965 | CTOT=CO3+HCO3+H2CO3 |
---|
| 966 | PH=-LOG10(AK1C*H2CO3/HCO3) |
---|
| 967 | ELSEIF(ICALC.EQ.5) THEN |
---|
| 968 | ! *** ALK AND PH FIXED *** |
---|
| 969 | HCO3=ALK/(1.0D0+2.0D0*AK2C/AHPLUS) |
---|
| 970 | CO3=HCO3*AK2C/AHPLUS |
---|
| 971 | H2CO3=HCO3*AHPLUS/AK1C |
---|
| 972 | PCO2=H2CO3/AKP |
---|
| 973 | CTOT=H2CO3+HCO3+CO3 |
---|
| 974 | ELSEIF(ICALC.EQ.6) THEN |
---|
| 975 | ! *** PCO2 AND PH FIXED *** |
---|
| 976 | H2CO3=PCO2*AKP |
---|
| 977 | HCO3=H2CO3*AK1C/AHPLUS |
---|
| 978 | CO3=HCO3*AK2C/AHPLUS |
---|
| 979 | CTOT=H2CO3+HCO3+CO3 |
---|
| 980 | ALK=HCO3+2.0D0*CO3 |
---|
| 981 | ENDIF |
---|
| 982 | ENDIF |
---|
| 983 | |
---|
| 984 | DO 120 II=1,NCONC |
---|
| 985 | CONCS(II)=CONCS2(II) |
---|
| 986 | 120 CONTINUE |
---|
| 987 | |
---|
| 988 | ! AXY (07/05/13) ================================================== |
---|
| 989 | ! C_SW IS SET AT 0 TO START |
---|
| 990 | ! THEN IF NON CONVERGENCE C_SW SET TO 1 |
---|
| 991 | 123 IF(C_SW.EQ.1)THEN |
---|
| 992 | ! IF NON CONVERGENCE, THE MODEL REQUIRES CONCS TO CONTAIN USABLE VALUES. |
---|
| 993 | ! BEST OFFER BEING THE OLD CONCS VALUES WHEN CONVERGENCE HAS BEEN |
---|
| 994 | ! ACHIEVED |
---|
| 995 | DO II=1,NCONC |
---|
| 996 | CONCS(II)=OLD_CONCS(II) |
---|
| 997 | END DO |
---|
| 998 | |
---|
| 999 | ! SPECIFIC CARBONATE VALUES TO PUSH CODE ON THROUGH THE |
---|
| 1000 | ! NON CONVERGENCE CONDITIONS |
---|
| 1001 | ! PCO2W = 0 SO CO2 UPTAKE WILL BE ENCOURAGED |
---|
| 1002 | ! CONCS(3)=O3C(III)*0.005_fp8/1.e6_fp8 |
---|
| 1003 | ! -99 IS A FLAG TO SHOW WHEN AND WHERE NON CONVERGENCE OCCURED |
---|
| 1004 | ! CONCS(4)=PH, OUTPUT OF -99 MEANS NON-CONVERGENCE |
---|
| 1005 | ! CONCS(4)=-99._fp8 |
---|
| 1006 | ! |
---|
| 1007 | ! AXY (10/05/13): remove this -99 flag; this is handled instead up in |
---|
| 1008 | ! trcbio_medusa.F90 where the iters variable is both |
---|
| 1009 | ! output and may be used to trigger action |
---|
| 1010 | ! |
---|
| 1011 | ! RESET SWITCH FOR NEXT CALL TO THIS SUBROUTINE |
---|
| 1012 | C_SW=0 |
---|
| 1013 | ENDIF |
---|
| 1014 | ! AXY (07/05/13) ================================================== |
---|
| 1015 | |
---|
| 1016 | |
---|
| 1017 | RETURN |
---|
| 1018 | |
---|
| 1019 | END SUBROUTINE |
---|
| 1020 | !----------------------------------------------------------------------- |
---|
| 1021 | |
---|
| 1022 | !======================================================================= |
---|
| 1023 | !======================================================================= |
---|
| 1024 | !======================================================================= |
---|
| 1025 | |
---|
| 1026 | !======================================================================= |
---|
| 1027 | ! |
---|
| 1028 | SUBROUTINE CaCO3_Saturation (Tc, S, D, CO3, & |
---|
| 1029 | Om_cal, Om_arg) |
---|
| 1030 | ! |
---|
| 1031 | !======================================================================= |
---|
| 1032 | ! |
---|
| 1033 | ! Routine to calculate the saturation state of calcite and aragonite |
---|
| 1034 | ! Inputs: |
---|
| 1035 | ! Tc Temperature (C) |
---|
| 1036 | ! S Salinity |
---|
| 1037 | ! D Depth (m) |
---|
| 1038 | ! CO3 Carbonate ion concentration (mol.kg-1 ie /1D6) |
---|
| 1039 | ! |
---|
| 1040 | ! Outputs |
---|
| 1041 | ! Om_cal Calite saturation |
---|
| 1042 | ! Om_arg Aragonite saturation |
---|
| 1043 | ! |
---|
| 1044 | ! Intermediates |
---|
| 1045 | ! K_cal Stoichiometric solubility product for calcite |
---|
| 1046 | ! K_arg Stoichiometric solubility product for aragonite |
---|
| 1047 | ! Ca Calcium 2+ concentration (mol.kg-1) |
---|
| 1048 | ! P Pressure (bars) |
---|
| 1049 | ! |
---|
| 1050 | ! Source |
---|
| 1051 | ! Zeebe & Wolf-Gladrow 2001 following Mucci (1983) |
---|
| 1052 | ! with pressure corrections from Millero (1995) |
---|
| 1053 | ! Code tested against reference values given in Z & W-G |
---|
| 1054 | ! Built Jerry Blackford, 2008 |
---|
| 1055 | !======================================================================= |
---|
| 1056 | |
---|
| 1057 | IMPLICIT None |
---|
| 1058 | REAL(wp), INTENT( in ) :: Tc, S, D, CO3 |
---|
| 1059 | REAL(wp), INTENT( inout ) :: Om_cal, Om_arg |
---|
| 1060 | REAL(wp) :: Tk, Kelvin, Ca |
---|
| 1061 | REAL(wp) :: logKspc, Kspc |
---|
| 1062 | REAL(wp) :: logKspa, Kspa |
---|
| 1063 | REAL(wp) :: tmp1, tmp2, tmp3 |
---|
| 1064 | REAL(wp) :: dV, dK, P, R |
---|
| 1065 | |
---|
| 1066 | ! setup |
---|
| 1067 | Kelvin = 273.15 |
---|
| 1068 | Tk = Tc + Kelvin |
---|
| 1069 | Ca = 0.01028 ! Currently oceanic mean value at S=25, needs refining) |
---|
| 1070 | Ca = 0.010279 * (S / 35.0) ! Ca varies with salinity (cf. Feeley et al., 2004; Yool et al., 2010) |
---|
| 1071 | R = 83.131 !(cm3.bar.mol-1.K-1) |
---|
| 1072 | P = D / 10.0 !pressure in bars |
---|
| 1073 | |
---|
| 1074 | ! calculate K for calcite |
---|
| 1075 | tmp1 = -171.9065 - (0.077993*Tk) + (2839.319/Tk) + 71.595*log10(Tk) |
---|
| 1076 | tmp2 = + (-0.77712 + (0.0028426*Tk) + (178.34/Tk))*SQRT(S) |
---|
| 1077 | tmp3 = - (0.07711*S) + (0.0041249*(S**1.5)) |
---|
| 1078 | logKspc = tmp1 + tmp2 + tmp3 |
---|
| 1079 | Kspc = 10.0**logKspc |
---|
| 1080 | |
---|
| 1081 | ! correction for pressure for calcite |
---|
| 1082 | IF ( D .GT. 0) THEN |
---|
| 1083 | dV = -48.76 + 0.5304*Tc |
---|
| 1084 | dK = -11.76/1.0D3 + (0.3692/1.0D3) * Tc |
---|
| 1085 | tmp1 = -(dV/(R*Tk))*P + (0.5*dK/(R*Tk))*P*P |
---|
| 1086 | Kspc = Kspc*exp(tmp1) |
---|
| 1087 | logKspc = log10(Kspc) |
---|
| 1088 | END IF |
---|
| 1089 | |
---|
| 1090 | ! calculate K for aragonite |
---|
| 1091 | tmp1 = -171.945 - 0.077993*Tk + 2903.293 / Tk + 71.595* log10(Tk) |
---|
| 1092 | tmp2 = + (-0.068393 + 0.0017276*Tk + 88.135/Tk)*SQRT(S) |
---|
| 1093 | tmp3 = - 0.10018*S + 0.0059415*S**1.5 |
---|
| 1094 | logKspa = tmp1 + tmp2 + tmp3 |
---|
| 1095 | Kspa = 10.0**logKspa |
---|
| 1096 | |
---|
| 1097 | ! correction for pressure for aragonite |
---|
| 1098 | IF ( D .GT. 0) THEN |
---|
| 1099 | dV = -46.00 + 0.5304*Tc |
---|
| 1100 | dK = -11.76/1.0D3 + (0.3692/1.0D3) * Tc |
---|
| 1101 | tmp1 = -(dV/(R*Tk))*P + (0.5*dK/(R*Tk))*P*P |
---|
| 1102 | Kspa = Kspa*exp(tmp1) |
---|
| 1103 | logKspa = log10(Kspa) |
---|
| 1104 | END IF |
---|
| 1105 | |
---|
| 1106 | ! calculate saturation states |
---|
| 1107 | Om_cal = (CO3 * Ca) / Kspc |
---|
| 1108 | Om_arg = (CO3 * Ca) / Kspa |
---|
| 1109 | |
---|
| 1110 | RETURN |
---|
| 1111 | |
---|
| 1112 | END SUBROUTINE |
---|
| 1113 | !----------------------------------------------------------------------- |
---|
| 1114 | |
---|
| 1115 | !======================================================================= |
---|
| 1116 | !======================================================================= |
---|
| 1117 | !======================================================================= |
---|
| 1118 | |
---|
| 1119 | # else |
---|
| 1120 | !!====================================================================== |
---|
| 1121 | !! Dummy module : No MEDUSA carbonate chemistry |
---|
| 1122 | !!====================================================================== |
---|
| 1123 | |
---|
| 1124 | CONTAINS |
---|
| 1125 | |
---|
| 1126 | SUBROUTINE trc_co2_medusa( kt ) ! Empty routine |
---|
| 1127 | |
---|
| 1128 | INTEGER, INTENT( in ) :: kt |
---|
| 1129 | |
---|
| 1130 | WRITE(*,*) 'trc_co2_medusa: You should not have seen this print! error?', kt |
---|
| 1131 | |
---|
| 1132 | END SUBROUTINE trc_co2_medusa |
---|
| 1133 | #endif |
---|
| 1134 | |
---|
| 1135 | !!====================================================================== |
---|
| 1136 | END MODULE trcco2_medusa |
---|