Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WOMBATlite-sinking: Spatially variable sinking rates within the generic tracers version of WOMBAT-lite #15

Open
wants to merge 23 commits into
base: main
Choose a base branch
from
Open

WOMBATlite-sinking: Spatially variable sinking rates within the generic tracers version of WOMBAT-lite #15

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
eb016cc
Added:
Nov 12, 2024
6460fc6
Added:
Nov 12, 2024
58c7bbf
Added:
Nov 14, 2024
3258e89
Trying to determine why Fe-limitation spreads too much across ocean. …
Nov 15, 2024
cdfe919
Figured out why dFe was becoming so limiting:
Nov 26, 2024
e80a2f7
* Syntax corrections and removal of unnecessary code snippets
Nov 27, 2024
70788dc
* minimum values of DIC and ALK
Dec 3, 2024
ab7863b
reverting syntax in generic_tracers_utils.F90 on request of Dougie Sq…
Dec 3, 2024
9b79d47
Problems with DIC air-sea flux solved: fix was to ensure upper and lo…
Dec 9, 2024
5017183
Updated some Fe cycling parameters to increase the loss rates of Fe i…
Dec 11, 2024
f3062bb
First pass of adding burial. Seems to work well so far.
Dec 11, 2024
56b8369
Computing omega_ara and omega_cal and added diagnostics to the output…
Dec 11, 2024
09eaaf0
Considering CaCO3 dynamics:
Dec 13, 2024
b995b37
Added dicp and dicr as new tracers (DIC preformed and DIC remineralised)
Dec 13, 2024
991e993
Added sedimentary dissolution of CaCO3 dependent on Omega.
Dec 15, 2024
bab1035
Included a hard floor for DICp tracer == to dic_min parameter
Dec 19, 2024
ba65270
Included parameter control via the field_table on the dissolution of …
Jan 9, 2025
8bb59eb
Major bug fix : implicit sinking scheme was incorrectly treating sink…
Jan 16, 2025
ffdf20b
Added ballasting of CaCO3 to sinking rate
Jan 16, 2025
4779c63
* Increased the CFL criterion for vertically variable sinking, and ma…
Jan 16, 2025
0723305
"conservetracers" logical option now active. Setting this to true mea…
Jan 21, 2025
2c21bd8
Updated the redistribution of tracer to be spread over the tile in qu…
Jan 21, 2025
9bf0d6d
Added if statement to ensure no division by zero during conservetrace…
Jan 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 110 additions & 4 deletions generic_tracers/FMS_ocmip2_co2calc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
! local variables
!
integer :: isc, iec, jsc, jec
integer :: i,j
integer :: i,j,ipc
real :: alpha_internal
real :: bt
real :: co2star_internal
Expand All @@ -229,17 +229,32 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
real :: ks
real :: ksi
real :: kw
real :: Kspa
real :: Kspc
real :: calcium
real :: log100
real :: s2
real :: scl
real :: sqrtis
real :: sqrts
real :: s15
real :: st
real :: tk
real :: tk100
real :: tk1002
real :: logf_of_s
real :: prb
real :: salinity
real :: R
real :: total2free_0p, total2free, ks_0p, kf_0p, free2SWS, total2SWS, SWS2total, free2SWS_0p, total2SWS_0p
real, dimension(12) :: deltav, deltak, lnkpok0
real, dimension(12) :: a0, a1, a2, b0, b1, b2
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/
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/
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/
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/
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/
DATA b2 /12*0.0/

character(len=10) :: co2_calc_method

Expand All @@ -254,8 +269,7 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
co2_calc_method = 'ocmip2'
end if



R = 83.14472

! Set the loop indices.
isc = dope_vec%isc ; iec = dope_vec%iec
Expand Down Expand Up @@ -364,8 +378,10 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
sqrtis = sqrt(max(0.0,is))
s2 = s_in(i,j) * s_in(i,j)
sqrts = sqrt(max(0.0,s_in(i,j)))
s15 = sqrts**3
scl = s_in(i,j) / 1.80655
logf_of_s = log(1.0 - 0.001005 * s_in(i,j))
prb = zt(i,j) / 10.0
!
! k0 from Weiss 1974
!
Expand Down Expand Up @@ -478,6 +494,85 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
!
ft = 0.000067 / 18.9984 * scl
!
! Omega saturation states
!
if (present(omega_arag)) then
Kspa = 10.0**(-171.945 - 0.077993*tk + 2903.293/tk + 71.595*LOG10(tk) &
+ (-0.068393 + 0.0017276*tk + 88.135/tk)*sqrts &
- 0.10018*s_in(i,j) + 0.0059415*s15 )
endif
if (present(omega_calc)) then
Kspc = 10.0**(-171.9065 - 0.077993*tk + 2839.319/tk + 71.595*LOG10(tk) &
+ (-0.77712 + 0.0028426*tk + 178.34/tk)*sqrts &
- 0.07711*s_in(i,j) + 0.0041249*s15 )
endif

!!! ------------------------------------------- !!!
!!! Pressure corrections to the above constants !!! (Pearse J. Buchanan, July 2024)
!!! ------------------------------------------- !!! (copied from mocsy)
! This must be done to ensure that the co2 sys !
! calculations are accurate beneath the surface !

do ipc=1,12
deltav(ipc) = a0(ipc) + a1(ipc)*t_in(i,j) + a2(ipc)*t_in(i,j)*t_in(i,j)
deltak(ipc) = b0(ipc) + b1(ipc)*t_in(i,j) + b2(ipc)*t_in(i,j)*t_in(i,j)
lnkpok0(ipc) = ( -deltav(ipc) + (0.5*deltak(ipc) * prb) )*prb/(R*tk)
enddo

! Conversion factor total -> free scale at pressure zero
total2free_0p = 1.0/(1.0 + st/ks) ! Kfree = Ktotal*total2free
ks_0p = ks*1
ks = ks * EXP(lnkpok0(5))
! Conversion factor total -> free scale
total2free = 1.0/(1.0 + st/ks) ! Kfree = Ktotal*total2free

kf_0p = kf * total2free_0p
kf = kf_0p * EXP(lnkpok0(6))
kf = kf / total2free

! Convert between seawater and total hydrogen (pH) scales
free2SWS = 1.0 + st/ks + ft/(kf*total2free) ! using Kf on free scale
total2SWS = total2free * free2SWS ! KSWS = Ktotal*total2SWS
SWS2total = 1.0 / total2SWS
! Conversion at pressure zero
free2SWS_0p = 1.0 + st/ks_0p + ft/(kf_0p) ! using Kf on free scale
total2SWS_0p = total2free_0p * free2SWS_0p ! KSWS = Ktotal*total2SWS

! Convert from Total to Seawater scale before pressure correction
! Must change to SEAWATER scale: Kb
kb = kb*total2SWS_0p

! Already on SEAWATER scale: K1p, K2p, K3p, Kb, Ksi, Kw

! Other contants (keep on another scale):
! - K0 (independent of pH scale, already pressure corrected)
! - Ks (already on Free scale; already pressure corrected)
! - Kf (already on Total scale; already pressure corrected)
! - Kspc, Kspa (independent of pH scale; pressure-corrected below)

! Perform actual pressure correction (on seawater scale)
k1 = k1*EXP(lnkpok0(1))
k2 = k2*EXP(lnkpok0(2))
kb = kb*EXP(lnkpok0(3))
kw = kw*EXP(lnkpok0(4))
Kspc = Kspc*EXP(lnkpok0(7))
Kspa = Kspa*EXP(lnkpok0(8))
k1p = k1p*EXP(lnkpok0(9))
k2p = k2p*EXP(lnkpok0(10))
k3p = k3p*EXP(lnkpok0(11))
ksi = ksi*EXP(lnkpok0(12))

! Convert back to original total scale:
k1 = k1 *SWS2total
k2 = k2 *SWS2total
k1p = k1p*SWS2total
k2p = k2p*SWS2total
k3p = k3p*SWS2total
kb = kb *SWS2total
ksi = ksi*SWS2total
kw = kw *SWS2total

!
!***********************************************************************
!
! Calculate [H+] total when DIC and TA are known at T, S and 1 atm.
Expand Down Expand Up @@ -515,6 +610,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
k1 * htotal(i,j) + k1 * k2)
if (present(co2star)) co2star(i,j) = co2star_internal
if (present(co3_ion)) co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
if (present(omega_arag) .or. present(omega_calc)) then
co3_ion(i,j) = co2star_internal * k1 * k2 / htotal2
calcium = (0.02128/40.078)*s_in(i,j)/1.80655
if (present(omega_arag)) omega_arag(i,j) = (calcium * co3_ion(i,j)) / Kspa
if (present(omega_calc)) omega_calc(i,j) = (calcium * co3_ion(i,j)) / Kspc
endif
!
! Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6
! values)
Expand Down Expand Up @@ -548,7 +649,12 @@ subroutine FMS_ocmip2_co2calc(dope_vec, mask, &
if (present(pco2surf)) then !{
pCO2surf(i,j) = 0.0
endif !}

if (present(omega_arag)) then
omega_arag(i,j) = 0.0
endif
if (present(omega_calc)) then
omega_calc(i,j) = 0.0
endif

endif !}mask

Expand Down
Loading