Skip to content

Commit b6f728a

Browse files
pearsebdougiesquire
pearseb
and
dougiesquire
authored
Fms ocmip2 co2calc update (#24)
* Extending "ocmip2" co2 calculations to depth: * Adding pressure effect to CO2 system calculations within ocmip2 option * Exposed omega_arag and omega_calc as outputs * Make pressure correction to ocmip2 carbon conversion constants optional Introduces a new optional argument to FMS_ocmip2_co2calc, called pcorr_ocmip2, that defaults to .false. --------- Co-authored-by: dougiesquire <[email protected]>
1 parent b66395e commit b6f728a

File tree

1 file changed

+128
-7
lines changed

1 file changed

+128
-7
lines changed

‎generic_tracers/FMS_ocmip2_co2calc.F90‎

Lines changed: 128 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -143,14 +143,17 @@ end subroutine read_mocsy_namelist
143143
!
144144
! co2_calc = carbonate formalation -- choices are 'ocmip2'(default)
145145
! or 'mocsy'
146+
! pcorr_ocmip2 = whether on not to apply pressure corrections to carbon
147+
! conversion constants. Defaults to false to preserve answers.
148+
! Only used if co2_calc='ocmip2'.
146149
!
147150
! OUTPUT
148151
! co2star = CO2*water, or H2CO3 concentration (mol/kg)
149152
! alpha = Solubility of CO2 for air (mol/kg/atm)
150153
! pco2surf = oceanic pCO2 (ppmv)
151154
! co3_ion = Carbonate ion, or CO3-- concentration (mol/kg)
152-
! omega_arag = aragonite saturation state (mol/kg ; avail. only w/ mocsy)
153-
! omega_calc = aragonite saturation state (mol/kg ; avail. only w/ mocsy)
155+
! omega_arag = aragonite saturation state (mol/kg)
156+
! omega_calc = aragonite saturation state (mol/kg)
154157
!
155158
! FILES and PROGRAMS NEEDED: drtsafe, ta_iter_1
156159
!
@@ -162,7 +165,7 @@ end subroutine read_mocsy_namelist
162165
subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
163166
t_in, s_in, dic_in, pt_in, sit_in, ta_in, htotallo, &
164167
htotalhi, htotal, co2_calc, zt, co2star, alpha, pCO2surf, &
165-
co3_ion, omega_arag, omega_calc) !{
168+
co3_ion, omega_arag, omega_calc, pcorr_ocmip2) !{
166169

167170
implicit none
168171

@@ -205,14 +208,17 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
205208
co3_ion, &
206209
omega_arag, &
207210
omega_calc
211+
logical, intent(in), optional :: pcorr_ocmip2
212+
208213
!
209214
! local variables
210215
!
211216
integer :: isc, iec, jsc, jec
212-
integer :: i,j
217+
integer :: i,j,ipc
213218
real :: alpha_internal
214219
real :: bt
215220
real :: co2star_internal
221+
real :: co3_ion_internal
216222
real :: dlogtk
217223
real :: ft
218224
real :: htotal2
@@ -229,19 +235,35 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
229235
real :: ks
230236
real :: ksi
231237
real :: kw
238+
real :: Kspa
239+
real :: Kspc
240+
real :: calcium
232241
real :: log100
233242
real :: s2
234243
real :: scl
235244
real :: sqrtis
236245
real :: sqrts
246+
real :: s15
237247
real :: st
238248
real :: tk
239249
real :: tk100
240250
real :: tk1002
241251
real :: logf_of_s
252+
real :: prb
242253
real :: salinity
254+
real :: R
255+
real :: total2free_0p, total2free, ks_0p, kf_0p, free2SWS, total2SWS, SWS2total, free2SWS_0p, total2SWS_0p
256+
real, dimension(12) :: deltav, deltak, lnkpok0
257+
real, dimension(12) :: a0, a1, a2, b0, b1, b2
258+
DATA a0 /-25.5, -15.82, -29.48, -20.02, -18.03, -9.78, -48.76, -45.96, -14.51, -23.12, -26.57, -29.48/
259+
DATA a1 /0.1271, -0.0219, 0.1622, 0.1119, 0.0466, -0.0090, 0.5304, 0.5304, 0.1211, 0.1758, 0.2020, 0.1622/
260+
DATA a2 /0.0, 0.0, -2.608e-3, -1.409e-3, 0.316e-3, -0.942e-3, 0.0, 0.0, -0.321e-3, -2.647e-3, -3.042e-3, -2.6080e-3/
261+
DATA b0 /-3.08e-3, 1.13e-3, -2.84e-3, -5.13e-3, -4.53e-3, -3.91e-3, -11.76e-3, -11.76e-3, -2.67e-3, -5.15e-3, -4.08e-3, -2.84e-3/
262+
DATA b1 /0.0877e-3, -0.1475e-3, 0.0, 0.0794e-3, 0.09e-3, 0.054e-3, 0.3692e-3, 0.3692e-3, 0.0427e-3, 0.09e-3, 0.0714e-3, 0.0/
263+
DATA b2 /12*0.0/
243264

244265
character(len=10) :: co2_calc_method
266+
logical :: pres_corr_ocmip2
245267

246268
if (present(co2_calc)) then
247269
co2_calc_method = trim(co2_calc)
@@ -254,8 +276,11 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
254276
co2_calc_method = 'ocmip2'
255277
end if
256278

257-
258-
279+
if (present(pcorr_ocmip2)) then
280+
pres_corr_ocmip2 = pcorr_ocmip2
281+
else
282+
pres_corr_ocmip2 = .false.
283+
end if
259284

260285
! Set the loop indices.
261286
isc = dope_vec%isc ; iec = dope_vec%iec
@@ -364,8 +389,10 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
364389
sqrtis = sqrt(max(0.0,is))
365390
s2 = s_in(i,j) * s_in(i,j)
366391
sqrts = sqrt(max(0.0,s_in(i,j)))
392+
s15 = sqrts**3
367393
scl = s_in(i,j) / 1.80655
368394
logf_of_s = log(1.0 - 0.001005 * s_in(i,j))
395+
prb = zt(i,j) / 10.0
369396
!
370397
! k0 from Weiss 1974
371398
!
@@ -478,6 +505,89 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
478505
!
479506
ft = 0.000067 / 18.9984 * scl
480507
!
508+
! Omega saturation states
509+
!
510+
if (present(omega_arag)) then
511+
Kspa = 10.0**(-171.945 - 0.077993*tk + 2903.293/tk + 71.595*LOG10(tk) &
512+
+ (-0.068393 + 0.0017276*tk + 88.135/tk)*sqrts &
513+
- 0.10018*s_in(i,j) + 0.0059415*s15 )
514+
endif
515+
if (present(omega_calc)) then
516+
Kspc = 10.0**(-171.9065 - 0.077993*tk + 2839.319/tk + 71.595*LOG10(tk) &
517+
+ (-0.77712 + 0.0028426*tk + 178.34/tk)*sqrts &
518+
- 0.07711*s_in(i,j) + 0.0041249*s15 )
519+
endif
520+
521+
!!! ------------------------------------------- !!!
522+
!!! Pressure corrections to the above constants !!! (Pearse J. Buchanan, July 2024)
523+
!!! ------------------------------------------- !!! (copied from mocsy)
524+
! This must be done to ensure that the co2 sys !
525+
! calculations are accurate beneath the surface !
526+
527+
if (pres_corr_ocmip2) then
528+
R = 83.14472
529+
530+
do ipc=1,12
531+
deltav(ipc) = a0(ipc) + a1(ipc)*t_in(i,j) + a2(ipc)*t_in(i,j)*t_in(i,j)
532+
deltak(ipc) = b0(ipc) + b1(ipc)*t_in(i,j) + b2(ipc)*t_in(i,j)*t_in(i,j)
533+
lnkpok0(ipc) = ( -deltav(ipc) + (0.5*deltak(ipc) * prb) )*prb/(R*tk)
534+
enddo
535+
536+
! Conversion factor total -> free scale at pressure zero
537+
total2free_0p = 1.0/(1.0 + st/ks) ! Kfree = Ktotal*total2free
538+
ks_0p = ks*1
539+
ks = ks * EXP(lnkpok0(5))
540+
! Conversion factor total -> free scale
541+
total2free = 1.0/(1.0 + st/ks) ! Kfree = Ktotal*total2free
542+
543+
kf_0p = kf * total2free_0p
544+
kf = kf_0p * EXP(lnkpok0(6))
545+
kf = kf / total2free
546+
547+
! Convert between seawater and total hydrogen (pH) scales
548+
free2SWS = 1.0 + st/ks + ft/(kf*total2free) ! using Kf on free scale
549+
total2SWS = total2free * free2SWS ! KSWS = Ktotal*total2SWS
550+
SWS2total = 1.0 / total2SWS
551+
! Conversion at pressure zero
552+
free2SWS_0p = 1.0 + st/ks_0p + ft/(kf_0p) ! using Kf on free scale
553+
total2SWS_0p = total2free_0p * free2SWS_0p ! KSWS = Ktotal*total2SWS
554+
555+
! Convert from Total to Seawater scale before pressure correction
556+
! Must change to SEAWATER scale: Kb
557+
kb = kb*total2SWS_0p
558+
559+
! Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw
560+
561+
! Other contants (keep on another scale):
562+
! - K0 (independent of pH scale, already pressure corrected)
563+
! - Ks (already on Free scale; already pressure corrected)
564+
! - Kf (already on Total scale; already pressure corrected)
565+
! - Kspc, Kspa (independent of pH scale; pressure-corrected below)
566+
567+
! Perform actual pressure correction (on seawater scale)
568+
k1 = k1*EXP(lnkpok0(1))
569+
k2 = k2*EXP(lnkpok0(2))
570+
kb = kb*EXP(lnkpok0(3))
571+
kw = kw*EXP(lnkpok0(4))
572+
Kspc = Kspc*EXP(lnkpok0(7))
573+
Kspa = Kspa*EXP(lnkpok0(8))
574+
k1p = k1p*EXP(lnkpok0(9))
575+
k2p = k2p*EXP(lnkpok0(10))
576+
k3p = k3p*EXP(lnkpok0(11))
577+
ksi = ksi*EXP(lnkpok0(12))
578+
579+
! Convert back to original total scale:
580+
k1 = k1 *SWS2total
581+
k2 = k2 *SWS2total
582+
k1p = k1p*SWS2total
583+
k2p = k2p*SWS2total
584+
k3p = k3p*SWS2total
585+
kb = kb *SWS2total
586+
ksi = ksi*SWS2total
587+
kw = kw *SWS2total
588+
endif
589+
590+
!
481591
!***********************************************************************
482592
!
483593
! Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
@@ -515,6 +625,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
515625
k1 * htotal(i,j) + k1 * k2)
516626
if (present(co2star)) co2star(i,j) = co2star_internal
517627
if (present(co3_ion)) co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
628+
if (present(omega_arag) .or. present(omega_calc)) then
629+
co3_ion_internal = co2star_internal * k1 * k2 / htotal2
630+
calcium = (0.02128/40.078)*s_in(i,j)/1.80655
631+
if (present(omega_arag)) omega_arag(i,j) = (calcium * co3_ion_internal) / Kspa
632+
if (present(omega_calc)) omega_calc(i,j) = (calcium * co3_ion_internal) / Kspc
633+
endif
518634
!
519635
! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6
520636
! values)
@@ -548,7 +664,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
548664
if (present(pco2surf)) then !{
549665
pCO2surf(i,j) = 0.0
550666
endif !}
551-
667+
if (present(omega_arag)) then
668+
omega_arag(i,j) = 0.0
669+
endif
670+
if (present(omega_calc)) then
671+
omega_calc(i,j) = 0.0
672+
endif
552673

553674
endif !}mask
554675

0 commit comments

Comments
 (0)