Incorrect area weights used in denominator of downsampled MMS diagnostics

[hdrake]

Downsampled diagnostics with MMS cell methods (i.e. "xh:mean yh:mean z_l:sum) are currently bugged, because incorrect weights are used in the weighted-averaging operation.

In the following lines, the weights are the product of a 3D mask array and a 2D area array.

MOM6/src/framework/MOM_diag_mediator.F90

Lines 4444 to 4457 in 4622d6a

	elseif (method == MMP .or. method == MMS) then !e.g., T_advection_xy
	do k=ks,ke ; do j=jsv_d,jev_d ; do i=isv_d,iev_d
	i0 = isv_o+dl*(i-isv_d)
	j0 = jsv_o+dl*(j-jsv_d)
	ave = 0.0
	total_weight = 0.0
	do jj=j0,j0+dl-1 ; do ii=i0,i0+dl-1
	! do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1
	weight = mask(ii,jj,k) * diag_cs%G%areaT(ii,jj)
	total_weight = total_weight + weight
	ave = ave+field_in(ii,jj,k)*weight
	enddo ; enddo
	field_out(i,j,k) = ave / (total_weight+eps_area) !Avoid zero mask at all aggregating cells where ave=0.0
	enddo ; enddo ; enddo

However, the 3D mask used to downsample the diagnostics seems to be a 2D wet mask that is simply projected downwards for all vertical levels:

MOM6/src/framework/MOM_diag_mediator.F90

Lines 3665 to 3705 in 4622d6a

	subroutine diag_masks_set(G, nz, diag_cs)
	type(ocean_grid_type), target, intent(in) :: G !< The ocean grid type.
	integer, intent(in) :: nz !< The number of layers in the model's native grid.
	type(diag_ctrl), pointer :: diag_cs !< A pointer to a type with many variables
	!! used for diagnostics
	! Local variables
	integer :: k
	! 2d masks point to the model masks since they are identical
	diag_cs%mask2dT => G%mask2dT
	diag_cs%mask2dBu => G%mask2dBu
	diag_cs%mask2dCu => G%mask2dCu
	diag_cs%mask2dCv => G%mask2dCv
	! 3d native masks are needed by diag_manager but the native variables
	! can only be masked 2d - for ocean points, all layers exists.
	allocate(diag_cs%mask3dTL(G%isd:G%ied,G%jsd:G%jed,1:nz))
	allocate(diag_cs%mask3dBL(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz))
	allocate(diag_cs%mask3dCuL(G%IsdB:G%IedB,G%jsd:G%jed,1:nz))
	allocate(diag_cs%mask3dCvL(G%isd:G%ied,G%JsdB:G%JedB,1:nz))
	do k=1,nz
	diag_cs%mask3dTL(:,:,k) = diag_cs%mask2dT(:,:)
	diag_cs%mask3dBL(:,:,k) = diag_cs%mask2dBu(:,:)
	diag_cs%mask3dCuL(:,:,k) = diag_cs%mask2dCu(:,:)
	diag_cs%mask3dCvL(:,:,k) = diag_cs%mask2dCv(:,:)
	enddo
	allocate(diag_cs%mask3dTi(G%isd:G%ied,G%jsd:G%jed,1:nz+1))
	allocate(diag_cs%mask3dBi(G%IsdB:G%IedB,G%JsdB:G%JedB,1:nz+1))
	allocate(diag_cs%mask3dCui(G%IsdB:G%IedB,G%jsd:G%jed,1:nz+1))
	allocate(diag_cs%mask3dCvi(G%isd:G%ied,G%JsdB:G%JedB,1:nz+1))
	do k=1,nz+1
	diag_cs%mask3dTi(:,:,k) = diag_cs%mask2dT(:,:)
	diag_cs%mask3dBi(:,:,k) = diag_cs%mask2dBu(:,:)
	diag_cs%mask3dCui(:,:,k) = diag_cs%mask2dCu(:,:)
	diag_cs%mask3dCvi(:,:,k) = diag_cs%mask2dCv(:,:)
	enddo
	!Allocate and initialize the downsampled masks
	call downsample_diag_masks_set(G, nz, diag_cs)
	end subroutine diag_masks_set

In vertical coordinates that include masked land cells, such as below topography in depth coordinates, however, this 3D area weight array will incorrectly have non-zero areas for any cells that are "wet" according to the 2D mask but actually should be "dry" in 3D because they are under the sea floor.

Since the numerator in the weighted average is a product of the weight with an MMS diagnostic, which should be zero under the seafloor anyways, there is no problem there. The problem arises in the total weight used in the denominator, which will erroneously include all of these sub-seafloor areas in the total.

[hdrake]

Here is some additional context as to the impact of this bug, and how it can be fixed in post-processing.

In the CM4X (0.125º ocean) simulations, vertically-resolved heat and salt budget diagnostics were only saved at the d2 downsampled resolution in depth coordinates. Heat budgets computed from these diagnostics appear to still "close", in the sense that the LHS and RHS are equal, but only because all terms in the budget have been scaled by the same incorrectly coarsened areas. Inspecting individual terms in the budget, however, reveals the problem. For instance, vertical advective or diffusive heat flux tendencies should integrate to zero over the entire column. Since this bug produces errors for cells adjacent to the sea floor, the downsampled diagnostics no longer correctly sum to zero (see left panel of Figure below).

However, these errors can be easily corrected in offline postprocessing by multiplying any 3D MMS diagnostics by a downsampled 3D wet mask that reflect the fraction of native cells that are wet, for each of the coarsened cells (middle and right panels of the figure).

[StephenGriffies]

So to be clear, is the problem fully resolvable via your offline fix, so that there is no need to rerun the simulation with the corrected code?

[hdrake]

Correct, the offline fix works, so there is no need to re-run any of the CM4X simulations.

Separately, we should decide how to communicate to users of the data that they are bugged but fixable offline. Can this fix be incorporated into the post-processing algorithms that generate the /pp/ directories?

[StephenGriffies]

Potentially, we could rerun the relevant pp/d2 generated files, once the pp code is fixed. The obvious person to champion this effort is @raphaeldussin. I am unsure if he can add this to his list before he goes on furlough.

Plan B is to rely on good communication/documentation. For that purpose, perhaps we should produce a CM4X Github repo that has documentation and further analysis code. I created this one for the CM4X paper

https://github.com/StephenGriffies/CM4X

But I think another one that is owned by you would be much more of a living repo, particularly since I am in transition to Europe and will not be on top of things for some months.

[hdrake]

Plan B sounds good for now. I use this repo (https://github.com/hdrake/CM4Xutils) for all of my personal CM4X analysis. It includes various convenience functions for loading in diagnostics efficiently but also some post-processing steps like extracting sea ice melt from the precipitation diagnostics and downsampling or transforming to different vertical coordinates. I will implement the fix in there.

[StephenGriffies]

Great.

Can you put a pointer to my CM4X paper repo? That will allow us to build up a web of repos for CM4X analysis.