Skip to content

WOMBATlite: Treat key biological loss terms semi-implicitly #97

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

Merged
dougiesquire merged 3 commits into from
Mar 25, 2026
Merged

WOMBATlite: Treat key biological loss terms semi-implicitly #97

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a874830
WOMBATlite: Treat key biological loss terms semi-implicitly
dougiesquire Mar 25, 2026
a443314
WOMBATlite: Add guards against some potential overflows
dougiesquire Mar 25, 2026
54b5107
Update comment as suggested in code review
dougiesquire Mar 25, 2026
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
155 changes: 97 additions & 58 deletions generic_tracers/generic_WOMBATlite.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1757,7 +1757,7 @@ subroutine generic_WOMBATlite_update_from_bottom(tracer_list, dt, tau, model_tim
do i = isc, iec
do j = jsc, jec
orgflux = wombat%det_btm(i,j) / dt * 86400 * 1e3 ! mmol C m-2 day-1
wombat%fbury(i,j) = 0.013 + 0.53 * orgflux**2.0 / (7.0 + orgflux)**2.0 ! Eq. 3 Dunne et al. 2007
wombat%fbury(i,j) = 0.013 + 0.53 * (orgflux / (7.0 + orgflux))**2.0 ! Eq. 3 Dunne et al. 2007
enddo
enddo
endif
Expand Down Expand Up @@ -1891,6 +1891,7 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
real :: biophy, biozoo, biodet, biono3, biofer, biocaco3
real :: biophyfe, biophy1, zooprefphy, zooprefdet, zooprey
real :: fbc, zval
real :: P_expl, k_loss, k_loss_zoodiss
real, parameter :: epsi = 1.0e-30
integer :: ichl
real :: par_phy_mldsum, par_z_mldsum
Expand Down Expand Up @@ -2598,28 +2599,21 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!

! NOTE: wombat%phymorq, wombat%zoomorq and wombat%detremi are now handled semi-implicitly
! and are calculated in Step 12

if (biophy>1e-3) then
wombat%phymorl(i,j,k) = wombat%phylmor * fbc * wombat%f_phy(i,j,k) ! [molC/kg/s]
wombat%phymorq(i,j,k) = wombat%phyqmor / mmol_m3_to_mol_kg * wombat%f_phy(i,j,k) * wombat%f_phy(i,j,k) ! [molC/kg/s]
else
wombat%phymorl(i,j,k) = 0.0
wombat%phymorq(i,j,k) = 0.0
endif

if (biozoo>1e-3) then
! reduce linear mortality (respiration losses) of zooplankton when there is low biomass
zoo_slmor = biozoo / (biozoo + wombat%zookz)
wombat%zoomorl(i,j,k) = wombat%zoolmor * fbc * wombat%f_zoo(i,j,k) * zoo_slmor ! [molC/kg/s]
wombat%zoomorq(i,j,k) = wombat%zooqmor / mmol_m3_to_mol_kg * wombat%f_zoo(i,j,k) * wombat%f_zoo(i,j,k) ! [molC/kg/s]
else
wombat%zoomorl(i,j,k) = 0.0
wombat%zoomorq(i,j,k) = 0.0
endif

if (wombat%f_det(i,j,k) > epsi) then
wombat%detremi(i,j,k) = wombat%reminr(i,j,k) / mmol_m3_to_mol_kg * wombat%f_det(i,j,k)**2.0 ! [molC/kg/s]
else
wombat%detremi(i,j,k) = 0.0
endif


Expand Down Expand Up @@ -2676,14 +2670,19 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!

! NOTE: wombat%zoodiss, wombat%caldiss, wombat%aradiss and wombat%pocdiss are now handled
! implicitly and are calculated in Step 12

if (do_caco3_dynamics) then
! PIC:POC ratio is a function of the substrate:inhibitor ratio, which is the
! HCO3- to free H+ ions ratio (mol/umol), following Lehmann & Bach (2024).
! We also add a T-dependent function to scale down CaCO3 production in waters colder
! than 3 degrees C based off the observation of no E hux growth beneath this (Fielding 2013; L&O)
hco3 = wombat%f_dic(i,j,k) - wombat%co3(i,j,k) - wombat%co2_star(i,j,k)
wombat%pic2poc(i,j,k) = min(0.3, (wombat%f_inorg + 10.0**(-3.0 + 4.31e-6 * &
hco3 / wombat%htotal(i,j,k))) * &
! Note the min(2.0, ...) in the exponent is just to prevent potential overflow. It does
! not affect answers as pic2poc is capped at 0.3
wombat%pic2poc(i,j,k) = min(0.3, (wombat%f_inorg + 10.0**(min(2.0, -3.0 + 4.31e-6 * &
hco3 / wombat%htotal(i,j,k)))) * &
(0.55 + 0.45 * tanh(Temp(i,j,k) - 4.0)) )
! The dissolution rate is a function of omegas for calcite and aragonite, as well the
! concentration of POC, following Kwon et al., 2024, Science Advances; Table S1, and
Expand All @@ -2698,18 +2697,6 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
wombat%dissratpoc(i,j,k) = 0.0
endif

if (wombat%f_caco3(i,j,k) > epsi) then
wombat%zoodiss(i,j,k) = wombat%zoograzdet(i,j,k) * wombat%fgutdiss * biocaco3/biodet
wombat%caldiss(i,j,k) = wombat%dissratcal(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
wombat%aradiss(i,j,k) = wombat%dissratara(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
wombat%pocdiss(i,j,k) = wombat%dissratpoc(i,j,k) * wombat%f_caco3(i,j,k) ! [mol/kg/s]
else
wombat%zoodiss(i,j,k) = 0.0
wombat%caldiss(i,j,k) = 0.0
wombat%aradiss(i,j,k) = 0.0
wombat%pocdiss(i,j,k) = 0.0
endif


!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!
Expand All @@ -2719,23 +2706,24 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &
!-----------------------------------------------------------------------!
!-----------------------------------------------------------------------!

! Nitrate equation ! [molN/kg]
!----------------------------------------------------------------------
wombat%f_no3(i,j,k) = wombat%f_no3(i,j,k) + dtsb * 16./122. * ( &
wombat%detremi(i,j,k) + &
wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) )

! Phytoplankton equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_phy(i,j,k) = wombat%f_phy(i,j,k) + dtsb * ( &
wombat%phygrow(i,j,k) - &
wombat%phymorl(i,j,k) - &
wombat%phymorq(i,j,k) - &
wombat%zoograzphy(i,j,k) )
! Here we treat phymorq semi-implicitly to ensure numerical stability
! See https://github.com/ACCESS-NRI/GFDL-generic-tracers/issues/96
P_expl = wombat%phygrow(i,j,k) - &
wombat%phymorl(i,j,k) - &
wombat%zoograzphy(i,j,k)
if (biophy>1e-3) then
! Treat phymorq semi-implicitly
k_loss = wombat%phyqmor / mmol_m3_to_mol_kg * wombat%f_phy(i,j,k)
wombat%f_phy(i,j,k) = (wombat%f_phy(i,j,k) + dtsb * P_expl) / (1.0 + dtsb * k_loss)

! Back-calculate phymorq for use below and for diagnostic output
wombat%phymorq(i,j,k) = k_loss * wombat%f_phy(i,j,k) ! [molC/kg/s]
else
wombat%f_phy(i,j,k) = wombat%f_phy(i,j,k) + dtsb * P_expl
wombat%phymorq(i,j,k) = 0.0
endif

! Phytoplankton chlorophyll equation ! [molChl/kg]
!-----------------------------------------------------------------------
Expand All @@ -2758,11 +2746,22 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &

! Zooplankton equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_zoo(i,j,k) = wombat%f_zoo(i,j,k) + dtsb * ( &
wombat%zooassiphy(i,j,k) + &
wombat%zooassidet(i,j,k) - &
wombat%zoomorl(i,j,k) - &
wombat%zoomorq(i,j,k) )
! Here we treat zoomorq semi-implicitly to ensure numerical stability
! See https://github.com/ACCESS-NRI/GFDL-generic-tracers/issues/96
P_expl = wombat%zooassiphy(i,j,k) + &
wombat%zooassidet(i,j,k) - &
wombat%zoomorl(i,j,k)
if (biozoo>1e-3) then
! Treat zoomorq semi-implicitly
k_loss = wombat%zooqmor / mmol_m3_to_mol_kg * wombat%f_zoo(i,j,k)
wombat%f_zoo(i,j,k) = (wombat%f_zoo(i,j,k) + dtsb * P_expl) / (1.0 + dtsb * k_loss)

! Back-calculate zoomorq for use below and diagnostic output
wombat%zoomorq(i,j,k) = k_loss * wombat%f_zoo(i,j,k) ! [molC/kg/s]
else
wombat%f_zoo(i,j,k) = wombat%f_zoo(i,j,k) + dtsb * P_expl
wombat%zoomorq(i,j,k) = 0.0
endif

! Zooplankton iron equation ! [molFe/kg]
!-----------------------------------------------------------------------
Expand All @@ -2779,13 +2778,34 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &

! Detritus equation ! [molC/kg]
!-----------------------------------------------------------------------
wombat%f_det(i,j,k) = wombat%f_det(i,j,k) + dtsb * ( &
wombat%zooegesphy(i,j,k) + &
wombat%zooegesdet(i,j,k) + &
wombat%phymorq(i,j,k) + &
wombat%zoomorq(i,j,k) - &
wombat%zoograzdet(i,j,k) - &
wombat%detremi(i,j,k) )
! Here we treat detremi semi-implicitly to ensure numerical stability
! See https://github.com/ACCESS-NRI/GFDL-generic-tracers/issues/96
P_expl = wombat%zooegesphy(i,j,k) + &
wombat%zooegesdet(i,j,k) + &
wombat%phymorq(i,j,k) + &
wombat%zoomorq(i,j,k) - &
wombat%zoograzdet(i,j,k)
if (wombat%f_det(i,j,k) > epsi) then
! Treat detremi semi-implicitly
k_loss = wombat%reminr(i,j,k) / mmol_m3_to_mol_kg * wombat%f_det(i,j,k)
wombat%f_det(i,j,k) = (wombat%f_det(i,j,k) + dtsb * P_expl) / (1.0 + dtsb * k_loss)

! Back-calculate detremi for use below and diagnostic output
wombat%detremi(i,j,k) = k_loss * wombat%f_det(i,j,k) ! [molC/kg/s]
else
wombat%f_det(i,j,k) = wombat%f_det(i,j,k) + dtsb * P_expl
wombat%detremi(i,j,k) = 0.0
endif

! Nitrate equation ! [molN/kg]
!----------------------------------------------------------------------
wombat%f_no3(i,j,k) = wombat%f_no3(i,j,k) + dtsb * 16./122. * ( &
wombat%detremi(i,j,k) + &
wombat%zoomorl(i,j,k) + &
wombat%zooexcrphy(i,j,k) + &
wombat%zooexcrdet(i,j,k) + &
wombat%phymorl(i,j,k) - &
wombat%phygrow(i,j,k) )

! Detrital iron equation ! [molFe/kg]
!-----------------------------------------------------------------------
Expand All @@ -2812,12 +2832,31 @@ subroutine generic_WOMBATlite_update_from_source(tracer_list, Temp, Salt, &

! Equation for CaCO3 ! [molCaCO3/kg]
!-----------------------------------------------------------------------
wombat%f_caco3(i,j,k) = wombat%f_caco3(i,j,k) + dtsb * ( &
wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) + &
wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) + &
wombat%zoograzphy(i,j,k) * (1. - wombat%fgutdiss) * wombat%pic2poc(i,j,k) - &
wombat%zoodiss(i,j,k) - wombat%caldiss(i,j,k) - &
wombat%aradiss(i,j,k) - wombat%pocdiss(i,j,k) )
! Here we treat the loss terms implicitly to ensure numerical stability
! See https://github.com/ACCESS-NRI/GFDL-generic-tracers/issues/96
! Treat the production terms explicitly
P_expl = wombat%phymorq(i,j,k) * wombat%pic2poc(i,j,k) + &
wombat%zoomorq(i,j,k) * wombat%pic2poc(i,j,k) + &
wombat%zoograzphy(i,j,k) * (1. - wombat%fgutdiss) * wombat%pic2poc(i,j,k)
if (wombat%f_caco3(i,j,k) > epsi) then
! Treat loss terms that are proportional to CaCO3 concentration implicitly
k_loss_zoodiss = wombat%zoograzdet(i,j,k) * wombat%fgutdiss / (biodet * mmol_m3_to_mol_kg)
k_loss = wombat%dissratcal(i,j,k) + wombat%dissratara(i,j,k) + wombat%dissratpoc(i,j,k) &
+ k_loss_zoodiss
wombat%f_caco3(i,j,k) = (wombat%f_caco3(i,j,k) + dtsb * P_expl) / (1.0 + dtsb * k_loss)

! Back-calculate the dissolution terms for use below and diagnostic output
wombat%zoodiss(i,j,k) = k_loss_zoodiss * wombat%f_caco3(i,j,k)
wombat%caldiss(i,j,k) = wombat%dissratcal(i,j,k) * wombat%f_caco3(i,j,k)
wombat%aradiss(i,j,k) = wombat%dissratara(i,j,k) * wombat%f_caco3(i,j,k)
wombat%pocdiss(i,j,k) = wombat%dissratpoc(i,j,k) * wombat%f_caco3(i,j,k)
else
wombat%f_caco3(i,j,k) = wombat%f_caco3(i,j,k) + dtsb * P_expl
wombat%zoodiss(i,j,k) = 0.0
wombat%caldiss(i,j,k) = 0.0
wombat%aradiss(i,j,k) = 0.0
wombat%pocdiss(i,j,k) = 0.0
end if

! Equation for DIC ! [molC/kg]
!-----------------------------------------------------------------------
Expand Down
Loading