Skip to content

Add online mixed-layer tracer budget diagnostics #79

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

Draft
rmholmes wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Add online mixed-layer tracer budget diagnostics #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b120e46
Add online mixed-layer tracer budget diagnostics
rmholmes Jun 21, 2026
ae88ad2
Remove erroneous merge lines.
rmholmes Jun 22, 2026
86c359c
Add additional mixed-layer averaged tracer diagnostics for sigma-diff…
rmholmes Jun 30, 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
257 changes: 236 additions & 21 deletions src/mom5/ocean_core/ocean_barotropic.F90
View file
Open in desktop

Large diffs are not rendered by default.

6 changes: 3 additions & 3 deletions src/mom5/ocean_core/ocean_model.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1604,7 +1604,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, &

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

call mpp_clock_end(id_flux_adjust)
Expand Down Expand Up @@ -1718,7 +1718,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, &
! smoother is not needed for C-grid.
if(vert_coordinate_class==DEPTH_BASED) then
call mpp_clock_begin(id_surface_smooth)
call ocean_eta_smooth(Time, Thickness, Ext_mode, T_prog(1:num_prog_tracers))
call ocean_eta_smooth(Time, Thickness, Dens, Ext_mode, T_prog(1:num_prog_tracers))
call mpp_clock_end(id_surface_smooth)
endif

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

! diagnose the geodepth of new blobs and the depth of old blobs
Expand Down
223 changes: 200 additions & 23 deletions src/mom5/ocean_core/ocean_sbc.F90
View file
Open in desktop

Large diffs are not rendered by default.

214 changes: 210 additions & 4 deletions src/mom5/ocean_diag/ocean_tracer_diag.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -351,6 +351,8 @@ module ocean_tracer_diag_mod
public calc_potrho_mixed_layer
public send_tracer_variance
public diagnose_eta_tend_3dflux
public compute_budget_mld
public compute_tracer_at_mlb

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


enddo

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

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

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


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

!#######################################################################
! <SUBROUTINE NAME="compute_budget_mld">
!
! <DESCRIPTION>
!
! Compute the MLD-average of a 3D tracer tendency field.
!
! A per-layer weight wrk1(i,j,k) in [0,1] is built so that the weighted
! sum over k equals the depth-integrated tendency within the mixed layer.
! Three cases handle where the MLD falls relative to the layer interfaces
! (depth_zwt):
! 1. MLD shallower than bottom of layer 1: partial weight for k=1 only.
! 2. MLD within layer k (k>1): full weight for layers above, partial
! weight for the straddling layer (fraction of dzt spanned by MLD).
! 3. MLD deeper than all model levels: full weight everywhere.
! The depth-weighted column sum is then divided by the MLD to give a
! depth-averaged tendency.
!
! Ryan.Holmes, September 2024
! </DESCRIPTION>
!
subroutine compute_budget_mld(Time, Thickness, Dens, T_prog, tendency, tendency_2d)

type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
real, dimension(isd:,jsd:,:), intent(in) :: tendency ! 3D tendency field (any units/m)
real, dimension(isd:,jsd:), intent(inout) :: tendency_2d ! MLD-averaged tendency (same units)

integer :: i,j,k,kp1,n
integer :: tau
real, dimension(isd:ied,jsd:jed) :: mld
real, parameter :: epsln=1.0e-20 ! guard against division by zero at land points

if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error from ocean_tracer_diag_mod (compute_budget_mld): module needs initialization')
endif

tau = Time%tau

call calc_mixed_layer_depth(Thickness, &
T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), &
T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), &
Dens%rho(isd:ied,jsd:jed,:,tau), &
Dens%pressure_at_depth(isd:ied,jsd:jed,:), &
mld(:,:), smooth_mld_input=.false.)

! Build the per-layer fractional weight wrk1(i,j,k).
! Case 1: MLD is shallower than the bottom of the first layer.
wrk1(:,:,:) = 0.0
k=1
do j=jsc,jec
do i=isc,iec
if(Grd%tmask(i,j,k)==1.0) then
if(Thickness%depth_zwt(i,j,k) >= mld(i,j)) then
wrk1(i,j,1) = mld(i,j)/Thickness%depth_zwt(i,j,k)
wrk1(i,j,2:nk) = 0.0
endif
endif
enddo
enddo

! Case 2: MLD straddles a layer interface at depth k (k>1).
! Layers above are fully included; the straddling layer gets a partial weight.
do j=jsc,jec
do i=isc,iec
kloopA: do k=2,nk
if(Grd%tmask(i,j,k)==1.0) then
if(Thickness%depth_zwt(i,j,k) >= mld(i,j) .and. &
Thickness%depth_zwt(i,j,k-1) < mld(i,j)) then
kp1 = min(k+1,nk)
wrk1(i,j,1:k-1) = 1.0
wrk1(i,j,k) = (mld(i,j)-Thickness%depth_zwt(i,j,k-1))/Thickness%dzt(i,j,k)
wrk1(i,j,kp1:nk) = 0.0
exit kloopA
endif
endif
enddo kloopA
enddo
enddo

! Case 3: MLD is deeper than all model levels — include the full column.
k=nk
do j=jsc,jec
do i=isc,iec
if(Grd%tmask(i,j,k)==1.0) then
if(Thickness%depth_zwt(i,j,k) <= mld(i,j)) then
wrk1(i,j,:) = 1.0
endif
endif
enddo
enddo

! Depth-weighted column sum, then divide by MLD to get depth-average.
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + wrk1(i,j,k)*tendency(i,j,k)
enddo
enddo
enddo
tendency_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
if (Grd%tmask(i,j,1)==1.0) then
tendency_2d(i,j) = wrk1_2d(i,j)/(mld(i,j) + epsln)
endif
enddo
enddo

end subroutine compute_budget_mld
! </SUBROUTINE> NAME="compute_budget_mld"

!#######################################################################
! <SUBROUTINE NAME="compute_tracer_at_mlb">
!
! <DESCRIPTION>
!
! Return the tracer value at the base of the mixed layer by linear
! interpolation between the depth cell centres (depth_zt) that bracket
! the MLD. Three edge cases are handled:
! 1. MLD shallower than the first cell centre: use the surface value.
! 2. MLD between cell centres k and k+1: distance-weighted interpolation
! (W1 = distance above MLD from k, W2 = distance below MLD from k+1;
! value = (tracer_k * W2 + tracer_k+1 * W1) / (W1+W2)).
! 3. MLD deeper than the deepest wet cell: use the bottom value.
! Also returns the MLD itself via mld_out for downstream diagnostics.
!
! Ryan.Holmes, November 2025
! </DESCRIPTION>
!
subroutine compute_tracer_at_mlb(Time, Thickness, Dens, T_prog, tracer, tracer_at_mlb, mld_out)

type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
type(ocean_prog_tracer_type), intent(in) :: T_prog(:)
real, dimension(isd:,jsd:,:), intent(in) :: tracer ! 3D tracer field
real, dimension(isd:,jsd:), intent(inout) :: tracer_at_mlb ! tracer interpolated to MLD base
real, dimension(isd:,jsd:), intent(inout) :: mld_out ! MLD (m), passed out for callers

integer :: i,j,k,kmt
integer :: tau
real :: W1, W2, denominator_r
real, dimension(isd:ied,jsd:jed) :: mld
real, parameter :: epsln=1.0e-20 ! guard against zero-thickness layers

if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error from ocean_tracer_diag_mod (compute_tracer_at_mlb): module needs initialization')
endif

tau = Time%tau

call calc_mixed_layer_depth(Thickness, &
T_prog(index_salt)%field(isd:ied,jsd:jed,:,tau), &
T_prog(index_temp)%field(isd:ied,jsd:jed,:,tau), &
Dens%rho(isd:ied,jsd:jed,:,tau), &
Dens%pressure_at_depth(isd:ied,jsd:jed,:), &
mld(:,:), smooth_mld_input=.false.)

tracer_at_mlb(:,:) = 0.0
mld_out(:,:) = mld(:,:)

! Case 1: MLD shallower than first cell centre — use surface value.
do j=jsc,jec
do i=isc,iec
if(mld(i,j) <= Thickness%depth_zt(i,j,1)) then
tracer_at_mlb(i,j) = tracer(i,j,1)
endif
enddo
enddo

! Case 2: MLD between cell centres k and k+1 — linear interpolation.
! W1 = distance from cell k to MLD; W2 = distance from MLD to cell k+1.
do j=jsc,jec
do i=isc,iec
kloop: do k=1,nk-1
if(Thickness%depth_zt(i,j,k) < mld(i,j) .and. mld(i,j) <= Thickness%depth_zt(i,j,k+1)) then
if(Grd%tmask(i,j,k+1) > 0) then
W1 = mld(i,j) - Thickness%depth_zt(i,j,k)
W2 = Thickness%depth_zt(i,j,k+1) - mld(i,j)
denominator_r = 1.0/(W1 + W2 + epsln)
tracer_at_mlb(i,j) = (tracer(i,j,k)*W2 + tracer(i,j,k+1)*W1)*denominator_r
exit kloop
endif
endif
enddo kloop
enddo
enddo

! Case 3: MLD deeper than the deepest wet cell — use the bottom value.
do j=jsc,jec
do i=isc,iec
kmt = Grd%kmt(i,j)
if(kmt > 0) then
if(mld(i,j) > Thickness%depth_zt(i,j,kmt)) then
tracer_at_mlb(i,j) = tracer(i,j,kmt)
endif
endif
enddo
enddo

end subroutine compute_tracer_at_mlb
! </SUBROUTINE> NAME="compute_tracer_at_mlb"

!#######################################################################
! <SUBROUTINE NAME="tracer_change">
Expand Down
27 changes: 23 additions & 4 deletions src/mom5/ocean_param/lateral/ocean_mixdownslope.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,7 @@ module ocean_mixdownslope_mod
use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_sum
use ocean_util_mod, only: write_chksum_2d, write_chksum_3d, write_chksum_2d_int
use ocean_tracer_util_mod, only: diagnose_3d_rho
use ocean_tracer_diag_mod, only: compute_budget_mld
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk5, wrk1_v

implicit none
Expand Down Expand Up @@ -187,8 +188,9 @@ module ocean_mixdownslope_mod
integer :: id_tform_salt_mixdown_on_nrho =-1


integer, dimension(:), allocatable :: id_mixdownslope
integer, dimension(:), allocatable :: id_mixdownslope_on_nrho
integer, dimension(:), allocatable :: id_mixdownslope
integer, dimension(:), allocatable :: id_mixdownslope_on_nrho
integer, dimension(:), allocatable :: id_mixdownslope_in_mld ! mixdownslope tendency averaged in mixed layer


character(len=128) :: version=&
Expand Down Expand Up @@ -536,8 +538,10 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio

allocate (id_mixdownslope(num_prog_tracers))
allocate (id_mixdownslope_on_nrho(num_prog_tracers))
allocate (id_mixdownslope_in_mld(num_prog_tracers))
id_mixdownslope_on_nrho = -1
id_mixdownslope = -1
id_mixdownslope_in_mld = -1

do n=1,num_prog_tracers
if(T_prog(n)%name == 'temp') then
Expand All @@ -551,6 +555,11 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio
Dens%neutralrho_axes(1:3), Time%model_time, &
'cp*mixdownslope*rho*dzt*temp binned to neutral density', &
'Watt/m^2', missing_value=missing_value, range=(/-1.e20,1.e20/))
id_mixdownslope_in_mld(n) = register_diag_field ('ocean_model', &
'mixdownslope_'//trim(T_prog(n)%name)//'_in_mld', &
Grd%tracer_axes(1:2), Time%model_time, &
'cp*mixdownslope*rho*dzt*temp averaged in mixed layer', &
'Watt/m^3', missing_value=missing_value, range=(/-1.e9,1.e9/))
else
id_mixdownslope(n) = register_diag_field ('ocean_model', 'mixdownslope_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), Time%model_time, &
Expand All @@ -561,6 +570,11 @@ subroutine ocean_mixdownslope_init(Grid, Domain, Time, Dens, T_prog, Ocean_optio
Dens%neutralrho_axes(1:3), Time%model_time, &
'mixdownslope*rho*dzt*tracer for '//trim(T_prog(n)%name)//' binned to neutral density', &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e20,1.e20/))
id_mixdownslope_in_mld(n) = register_diag_field ('ocean_model', &
'mixdownslope_'//trim(T_prog(n)%name)//'_in_mld', &
Grd%tracer_axes(1:2), Time%model_time, &
'mixdownslope*rho*dzt*tracer averaged in MLD for '//trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units)//'/m', missing_value=missing_value, range=(/-1.e9,1.e9/))
endif
enddo

Expand Down Expand Up @@ -650,6 +664,7 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt)
real :: mass_sum, tmix
real :: mixdownslope_total, mixdownslope_total_r
real :: tendency
real, dimension(isd:ied,jsd:jed) :: tendency_in_mld

integer :: stdoutunit
stdoutunit=stdout()
Expand Down Expand Up @@ -907,10 +922,10 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt)
enddo
enddo

if(id_mixdownslope(nt) > 0 .or. id_mixdownslope_on_nrho(nt) > 0) then
if(id_mixdownslope(nt) > 0 .or. id_mixdownslope_on_nrho(nt) > 0 .or. id_mixdownslope_in_mld(nt) > 0) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = tend_mix(i,j,k)*mixdownslope_mask(i,j)*T_prog(nt)%conversion
enddo
Expand All @@ -922,6 +937,10 @@ subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt)
if (id_mixdownslope_on_nrho(nt) > 0) then
call diagnose_3d_rho(Time, Dens, id_mixdownslope_on_nrho(nt), wrk1)
endif
if (id_mixdownslope_in_mld(nt) > 0) then
call compute_budget_mld(Time, Thickness, Dens, T_prog, wrk1(:,:,:), tendency_in_mld(:,:))
call diagnose_2d(Time, Grd, id_mixdownslope_in_mld(nt), tendency_in_mld(:,:))
endif
endif

if(nt==index_temp) then
Expand Down
Loading
Loading