Skip to content

Commit 12f398b

Browse files
rmholmesclaude
rmholmes
and
claude
committed
Add online mixed-layer tracer budget diagnostics
Integrates three-dimensional tracer budget process terms across the mixed layer at every model time step and divides by the mixed layer depth, enabling accurate online computation of the mixed-layer tracer budget (Holmes and Malan, article submitted to JAMES). Uses some existing mixed layer routines in ocean_tracer_diag_mod. Two core routines added to ocean_tracer_diag_mod: compute_budget_mld — MLD-averages a 3D tendency field compute_tracer_at_mlb — linearly interpolates a tracer to the MLD base Mixed layer budget tracer diagnostics: temp/salt_in_mld — MLD-averaged tracer concentration (times rho0) temp/salt_at_mlb — tracer at the MLD base (C_ent) MLD-averaged tracer budget process diagnostics: temp/salt_tendency_in_mld — total tracer tendency temp/salt_advection_in_mld — resolved advective convergence temp/salt_submeso_in_mld — submesoscale parameterized advection temp/salt_vdiffuse_diff_cbt_in_mld — vertical diffusion temp/salt_vdiffuse_diff_cbt_conv_in_mld — convective vertical diffusion temp/salt_vdiffuse_k33_in_mld — along-isopycnal K33 mixing temp/salt_vdiffuse_sbc_in_mld — surface boundary condition flux temp/salt_nonlocal_KPP_in_mld — nonlocal KPP temp/salt_rivermix_in_mld — river mixing neutral_physics_ndiffuse_in_mld — neutral diffusion neutral_physics_gm_in_mld — GM mesoscale eddy parameterization frazil_3d_in_mld — frazil ice formation sw_heat_in_mld — shortwave penetration below surface swflx/lw_heat/sens_heat/evap_heat_in_mld — individual surface heat flux components net_sfc_heating_in_mld — total surface heat flux sfc_hflux_pme_in_mld — P-E+R tracer flux at surface pme_river_in_mld — P-E+R mass flux within MLD temp/salt_eta_smooth_in_mld — SSH smoother tendency within MLD GVC-to-Eulerian correction terms (see Holmes and Malan appendix): eta_t_tendency_times_temp/salt_in_mld — ∂η/∂t × MLD-avg tracer pme_river_times_temp/salt_in_mld — P-E+R × MLD-avg tracer eta_smoother_times_temp/salt_in_mld — SSH smoother × MLD-avg tracer s_surf_ent_temp/salt — z* coordinate entrainment at MLD base: (1-H/(D+η)) × ∂η/∂t × C_ent Testing diagnostics: eta_t_tendency_times_temp/salt_at_mlb — ∂η/∂t × tracer at MLD base For more information including analysis code see: - https://github.com/rmholmes/access-om2-sst-budget, in particular the Theory_and_Diagnostics.ipynb notebook. - Holmes and Malan, journal article submitted to JAMES. The changes in this commit were developed on an older version of MOM5. Claude Sonnet 4.6 assisted in porting them onto the current ACCESS-NRI MOM5 codebase (upstream/upstream-master) by extracting a patch from the development branch, applying it with conflict resolution, and improving inline Fortran documentation in the new code. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
1 parent 59e0aad commit 12f398b

13 files changed

Lines changed: 976 additions & 65 deletions

‎src/mom5/ocean_core/ocean_barotropic.F90‎

Lines changed: 236 additions & 21 deletions
Large diffs are not rendered by default.

‎src/mom5/ocean_core/ocean_model.F90‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1604,7 +1604,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, &
16041604

16051605
! compute "flux adjustments" (e.g., surface tracer restoring, flux correction)
16061606
call mpp_clock_begin(id_flux_adjust)
1607-
call flux_adjust(Time, T_diag(1:num_diag_tracers), Dens, Ext_mode, &
1607+
call flux_adjust(Time, T_diag(1:num_diag_tracers), Dens, Thickness, Ext_mode, &
16081608
T_prog(1:num_prog_tracers), Velocity, river, melt, pme)
16091609

16101610
call mpp_clock_end(id_flux_adjust)
@@ -1718,7 +1718,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, &
17181718
! smoother is not needed for C-grid.
17191719
if(vert_coordinate_class==DEPTH_BASED) then
17201720
call mpp_clock_begin(id_surface_smooth)
1721-
call ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog(1:num_prog_tracers))
1721+
call ocean_eta_smooth(Time, Thickness, Dens, Ext_mode, T_prog(1:num_prog_tracers))
17221722
call mpp_clock_end(id_surface_smooth)
17231723
endif
17241724

@@ -1879,7 +1879,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, &
18791879
! diagnose time=taup1 ocean free surface height or bottom pressure.
18801880
! also diagnose geodepth_zt and geodepth_zwt.
18811881
call mpp_clock_begin(id_eta_and_pbot_diagnose)
1882-
call eta_and_pbot_diagnose(Time, Dens, Thickness, patm, pme, river, Ext_mode, Lagrangian_system, use_blobs)
1882+
call eta_and_pbot_diagnose(Time, Dens, Thickness, T_prog, patm, pme, river, Ext_mode, Lagrangian_system, use_blobs)
18831883
call mpp_clock_end(id_eta_and_pbot_diagnose)
18841884

18851885
! diagnose the geodepth of new blobs and the depth of old blobs

‎src/mom5/ocean_core/ocean_sbc.F90‎

Lines changed: 200 additions & 23 deletions
Large diffs are not rendered by default.

‎src/mom5/ocean_diag/ocean_tracer_diag.F90‎

Lines changed: 210 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -351,6 +351,8 @@ module ocean_tracer_diag_mod
351351
public calc_potrho_mixed_layer
352352
public send_tracer_variance
353353
public diagnose_eta_tend_3dflux
354+
public compute_budget_mld
355+
public compute_tracer_at_mlb
354356

355357
private compute_subduction
356358
private compute_tracer_mld
@@ -694,8 +696,6 @@ subroutine ocean_tracer_diag_init(Grid, Domain, Time, Time_steps, Thickness, T_p
694696
'Vertically integrated tracer [sum(rho_dzt*tracer)] in mixed layer for '//trim(T_prog(n)%name),&
695697
'kg/m^2 * '//trim(T_prog(n)%units), missing_value=missing_value, range=(/-1.e20,1.e20/))
696698
if(id_tracer_mld(n) > 0) compute_tracer_mld_diags = .true.
697-
698-
699699
enddo
700700

701701
id_kappa_sort = register_diag_field ('ocean_model','kappa_sort', &
@@ -1072,7 +1072,7 @@ subroutine calc_mixed_layer_depth(Thickness, salinity, theta, rho, pressure, &
10721072
real, dimension(isd:,jsd:), intent(inout) :: hmxl
10731073
logical, optional, intent(in) :: smooth_mld_input
10741074

1075-
real, parameter :: epsln=1.0e-20 ! for divisions
1075+
real, parameter :: epsln=1.0e-20 ! for divisions
10761076
integer :: i, j, k, km1, kb
10771077
logical :: smooth_mld_routine
10781078

@@ -1672,7 +1672,6 @@ subroutine compute_tracer_mld(Time, Thickness, Dens, T_prog, salinity, theta)
16721672
integer :: i,j,k,kp1,n
16731673
integer :: taup1
16741674
real, dimension(isd:ied,jsd:jed) :: mld
1675-
16761675

16771676
if (.not.module_is_initialized) then
16781677
call mpp_error(FATAL, &
@@ -1755,7 +1754,214 @@ subroutine compute_tracer_mld(Time, Thickness, Dens, T_prog, salinity, theta)
17551754
end subroutine compute_tracer_mld
17561755
! </SUBROUTINE> NAME="compute_tracer_mld"
17571756

1757+
!#######################################################################
1758+
! <SUBROUTINE NAME="compute_budget_mld">
1759+
!
1760+
! <DESCRIPTION>
1761+
!
1762+
! Compute the MLD-average of a 3D tracer tendency field.
1763+
!
1764+
! A per-layer weight wrk1(i,j,k) in [0,1] is built so that the weighted
1765+
! sum over k equals the depth-integrated tendency within the mixed layer.
1766+
! Three cases handle where the MLD falls relative to the layer interfaces
1767+
! (depth_zwt):
1768+
! 1. MLD shallower than bottom of layer 1: partial weight for k=1 only.
1769+
! 2. MLD within layer k (k>1): full weight for layers above, partial
1770+
! weight for the straddling layer (fraction of dzt spanned by MLD).
1771+
! 3. MLD deeper than all model levels: full weight everywhere.
1772+
! The depth-weighted column sum is then divided by the MLD to give a
1773+
! depth-averaged tendency.
1774+
!
1775+
! Ryan.Holmes, September 2024
1776+
! </DESCRIPTION>
1777+
!
1778+
subroutine compute_budget_mld(Time, Thickness, Dens, T_prog, tendency, tendency_2d)
1779+
1780+
type(ocean_time_type), intent(in) :: Time
1781+
type(ocean_thickness_type), intent(in) :: Thickness
1782+
type(ocean_density_type), intent(in) :: Dens
1783+
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
1784+
real, dimension(isd:,jsd:,:), intent(in) :: tendency ! 3D tendency field (any units/m)
1785+
real, dimension(isd:,jsd:), intent(inout) :: tendency_2d ! MLD-averaged tendency (same units)
1786+
1787+
integer :: i,j,k,kp1,n
1788+
integer :: tau
1789+
real, dimension(isd:ied,jsd:jed) :: mld
1790+
real, parameter :: epsln=1.0e-20 ! guard against division by zero at land points
1791+
1792+
if (.not.module_is_initialized) then
1793+
call mpp_error(FATAL, &
1794+
'==>Error from ocean_tracer_diag_mod (compute_budget_mld): module needs initialization')
1795+
endif
1796+
1797+
tau = Time%tau
1798+
1799+
call calc_mixed_layer_depth(Thickness, &
1800+
T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), &
1801+
T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), &
1802+
Dens%rho(isd:ied,jsd:jed,:,tau), &
1803+
Dens%pressure_at_depth(isd:ied,jsd:jed,:), &
1804+
mld(:,:), smooth_mld_input=.false.)
1805+
1806+
! Build the per-layer fractional weight wrk1(i,j,k).
1807+
! Case 1: MLD is shallower than the bottom of the first layer.
1808+
wrk1(:,:,:) = 0.0
1809+
k=1
1810+
do j=jsc,jec
1811+
do i=isc,iec
1812+
if(Grd%tmask(i,j,k)==1.0) then
1813+
if(Thickness%depth_zwt(i,j,k) >= mld(i,j)) then
1814+
wrk1(i,j,1) = mld(i,j)/Thickness%depth_zwt(i,j,k)
1815+
wrk1(i,j,2:nk) = 0.0
1816+
endif
1817+
endif
1818+
enddo
1819+
enddo
1820+
1821+
! Case 2: MLD straddles a layer interface at depth k (k>1).
1822+
! Layers above are fully included; the straddling layer gets a partial weight.
1823+
do j=jsc,jec
1824+
do i=isc,iec
1825+
kloopA: do k=2,nk
1826+
if(Grd%tmask(i,j,k)==1.0) then
1827+
if(Thickness%depth_zwt(i,j,k) >= mld(i,j) .and. &
1828+
Thickness%depth_zwt(i,j,k-1) < mld(i,j)) then
1829+
kp1 = min(k+1,nk)
1830+
wrk1(i,j,1:k-1) = 1.0
1831+
wrk1(i,j,k) = (mld(i,j)-Thickness%depth_zwt(i,j,k-1))/Thickness%dzt(i,j,k)
1832+
wrk1(i,j,kp1:nk) = 0.0
1833+
exit kloopA
1834+
endif
1835+
endif
1836+
enddo kloopA
1837+
enddo
1838+
enddo
1839+
1840+
! Case 3: MLD is deeper than all model levels — include the full column.
1841+
k=nk
1842+
do j=jsc,jec
1843+
do i=isc,iec
1844+
if(Grd%tmask(i,j,k)==1.0) then
1845+
if(Thickness%depth_zwt(i,j,k) <= mld(i,j)) then
1846+
wrk1(i,j,:) = 1.0
1847+
endif
1848+
endif
1849+
enddo
1850+
enddo
1851+
1852+
! Depth-weighted column sum, then divide by MLD to get depth-average.
1853+
wrk1_2d(:,:) = 0.0
1854+
do k=1,nk
1855+
do j=jsc,jec
1856+
do i=isc,iec
1857+
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)*tendency(i,j,k)
1858+
enddo
1859+
enddo
1860+
enddo
1861+
tendency_2d(:,:) = 0.0
1862+
do j=jsc,jec
1863+
do i=isc,iec
1864+
if (Grd%tmask(i,j,1)==1.0) then
1865+
tendency_2d(i,j) = wrk1_2d(i,j)/(mld(i,j) + epsln)
1866+
endif
1867+
enddo
1868+
enddo
1869+
1870+
end subroutine compute_budget_mld
1871+
! </SUBROUTINE> NAME="compute_budget_mld"
1872+
1873+
!#######################################################################
1874+
! <SUBROUTINE NAME="compute_tracer_at_mlb">
1875+
!
1876+
! <DESCRIPTION>
1877+
!
1878+
! Return the tracer value at the base of the mixed layer by linear
1879+
! interpolation between the depth cell centres (depth_zt) that bracket
1880+
! the MLD. Three edge cases are handled:
1881+
! 1. MLD shallower than the first cell centre: use the surface value.
1882+
! 2. MLD between cell centres k and k+1: distance-weighted interpolation
1883+
! (W1 = distance above MLD from k, W2 = distance below MLD from k+1;
1884+
! value = (tracer_k * W2 + tracer_k+1 * W1) / (W1+W2)).
1885+
! 3. MLD deeper than the deepest wet cell: use the bottom value.
1886+
! Also returns the MLD itself via mld_out for downstream diagnostics.
1887+
!
1888+
! Ryan.Holmes, November 2025
1889+
! </DESCRIPTION>
1890+
!
1891+
subroutine compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tracer, tracer_at_mlb, mld_out)
1892+
1893+
type(ocean_time_type), intent(in) :: Time
1894+
type(ocean_thickness_type), intent(in) :: Thickness
1895+
type(ocean_density_type), intent(in) :: Dens
1896+
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
1897+
real, dimension(isd:,jsd:,:), intent(in) :: tracer ! 3D tracer field
1898+
real, dimension(isd:,jsd:), intent(inout) :: tracer_at_mlb ! tracer interpolated to MLD base
1899+
real, dimension(isd:,jsd:), intent(inout) :: mld_out ! MLD (m), passed out for callers
1900+
1901+
integer :: i,j,k,kmt
1902+
integer :: tau
1903+
real :: W1, W2, denominator_r
1904+
real, dimension(isd:ied,jsd:jed) :: mld
1905+
real, parameter :: epsln=1.0e-20 ! guard against zero-thickness layers
1906+
1907+
if (.not.module_is_initialized) then
1908+
call mpp_error(FATAL, &
1909+
'==>Error from ocean_tracer_diag_mod (compute_tracer_at_mlb): module needs initialization')
1910+
endif
1911+
1912+
tau = Time%tau
1913+
1914+
call calc_mixed_layer_depth(Thickness, &
1915+
T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), &
1916+
T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), &
1917+
Dens%rho(isd:ied,jsd:jed,:,tau), &
1918+
Dens%pressure_at_depth(isd:ied,jsd:jed,:), &
1919+
mld(:,:), smooth_mld_input=.false.)
1920+
1921+
tracer_at_mlb(:,:) = 0.0
1922+
mld_out(:,:) = mld(:,:)
1923+
1924+
! Case 1: MLD shallower than first cell centre — use surface value.
1925+
do j=jsc,jec
1926+
do i=isc,iec
1927+
if(mld(i,j) <= Thickness%depth_zt(i,j,1)) then
1928+
tracer_at_mlb(i,j) = tracer(i,j,1)
1929+
endif
1930+
enddo
1931+
enddo
1932+
1933+
! Case 2: MLD between cell centres k and k+1 — linear interpolation.
1934+
! W1 = distance from cell k to MLD; W2 = distance from MLD to cell k+1.
1935+
do j=jsc,jec
1936+
do i=isc,iec
1937+
kloop: do k=1,nk-1
1938+
if(Thickness%depth_zt(i,j,k) < mld(i,j) .and. mld(i,j) <= Thickness%depth_zt(i,j,k+1)) then
1939+
if(Grd%tmask(i,j,k+1) > 0) then
1940+
W1 = mld(i,j) - Thickness%depth_zt(i,j,k)
1941+
W2 = Thickness%depth_zt(i,j,k+1) - mld(i,j)
1942+
denominator_r = 1.0/(W1 + W2 + epsln)
1943+
tracer_at_mlb(i,j) = (tracer(i,j,k)*W2 + tracer(i,j,k+1)*W1)*denominator_r
1944+
exit kloop
1945+
endif
1946+
endif
1947+
enddo kloop
1948+
enddo
1949+
enddo
1950+
1951+
! Case 3: MLD deeper than the deepest wet cell — use the bottom value.
1952+
do j=jsc,jec
1953+
do i=isc,iec
1954+
kmt = Grd%kmt(i,j)
1955+
if(kmt > 0) then
1956+
if(mld(i,j) > Thickness%depth_zt(i,j,kmt)) then
1957+
tracer_at_mlb(i,j) = tracer(i,j,kmt)
1958+
endif
1959+
endif
1960+
enddo
1961+
enddo
17581962

1963+
end subroutine compute_tracer_at_mlb
1964+
! </SUBROUTINE> NAME="compute_tracer_at_mlb"
17591965

17601966
!#######################################################################
17611967
! <SUBROUTINE NAME="tracer_change">

0 commit comments

Comments
 (0)