Modify remapping of vertically extensive diagnostics

Vertically extensive diagnostics like advective transports, diffusive
fluxes, etc. were being remapped onto a diagnostic grid prior to the
physical process that chagne the tracer and/or layer thicknesses.
However this means that the effective operator used for each component
of a tracer budget was not the same, and so budgets could never be
closed cell-by-cell in a diagnostic coordinate.

This commit fixes one part of the inconsistency by setting the target
diagnostic grid for all vertically extensive quantities to be one
constructed at the very beginning of the timestep.

However, problems with closing the budget should still be expected
1. Conservation of column integrals cannot be expected because the
   target grid can be smaller than the source grid (layer thicknesses
   at the native grid at the current point in the algorithm). This
   leads to some fluxes being 'thrown away' due to the reintegrate
   algorithm.
2. To truly have a consistent operator for all terms of a budget,
   the source grid for all vertically extensive grids should also
   be the same

src/core/MOM.F90

@@ -622,6 +622,10 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_int_in, CS, &
 
     if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n)
 
+    ! Update the vertically extensive diagnostic grids so that they are
+    ! referenced to the beginning timestep
+    call diag_update_remap_grids(CS%diag, update_intensive = .false., update_extensive = .true. )
+
     !===========================================================================
     ! This is the first place where the diabatic processes and remapping could occur.
     if (CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0) .and. do_thermo) then ! do thermodynamics.

src/framework/MOM_diag_mediator.F90

@@ -1473,14 +1473,16 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)
   logical :: staggered_in_x, staggered_in_y
   real, dimension(:,:,:), pointer :: h_diag => NULL()
 
+  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
+
+  ! For intensive variables only, we can choose to use a different diagnostic grid
+  ! to map to
   if (present(alt_h)) then
     h_diag => alt_h
   else
     h_diag => diag_cs%h
   endif
 
-  if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator)
-
   ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical
   ! grids, and post each.
   call assert(diag_field_id < diag_cs%next_free_diag_id, &
@@ -1500,10 +1502,11 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask, alt_h)
 
       if (id_clock_diag_remap>0) call cpu_clock_begin(id_clock_diag_remap)
       allocate(remapped_field(size(field,1), size(field,2), diag%axes%nz))
-      call vertically_reintegrate_diag_field( &
-              diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &
-              diag_cs%G, h_diag, staggered_in_x, staggered_in_y, &
-              diag%axes%mask3d, diag_cs%missing_value, field, remapped_field)
+      call vertically_reintegrate_diag_field(                                    &
+        diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), diag_cs%G,  &
+        diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number)%h_extensive, &
+        staggered_in_x, staggered_in_y, diag%axes%mask3d, diag_cs%missing_value, &
+        field, remapped_field)
       if (id_clock_diag_remap>0) call cpu_clock_end(id_clock_diag_remap)
       if (associated(diag%axes%mask3d)) then
         ! Since 3d masks do not vary in the vertical, just use as much as is
@@ -3191,19 +3194,25 @@ subroutine diag_set_state_ptrs(h, T, S, eqn_of_state, diag_cs)
 !> Build/update vertical grids for diagnostic remapping.
 !! \note The target grids need to be updated whenever sea surface
 !! height changes.
-subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S)
+subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S, update_intensive, update_extensive )
   type(diag_ctrl),        intent(inout) :: diag_cs      !< Diagnostics control structure
   real, target, optional, intent(in   ) :: alt_h(:,:,:) !< Used if remapped grids should be something other than
                                                         !! the current thicknesses [H ~> m or kg m-2]
   real, target, optional, intent(in   ) :: alt_T(:,:,:) !< Used if remapped grids should be something other than
                                                         !! the current temperatures
   real, target, optional, intent(in   ) :: alt_S(:,:,:) !< Used if remapped grids should be something other than
                                                         !! the current salinity
+  logical, optional,      intent(in   ) :: update_intensive !< If true (default), update the grids used for
+                                                            !! intensive diagnostics
+  logical, optional,      intent(in   ) :: update_extensive !< If true (not default), update the grids used for
+                                                            !! intensive diagnostics
   ! Local variables
   integer :: i
   real, dimension(:,:,:), pointer :: h_diag => NULL()
   real, dimension(:,:,:), pointer :: T_diag => NULL(), S_diag => NULL()
+  logical :: update_intensive_local, update_extensive_local
 
+  ! Set values based on optional input arguments
   if (present(alt_h)) then
     h_diag => alt_h
   else
@@ -3222,18 +3231,34 @@ subroutine diag_update_remap_grids(diag_cs, alt_h, alt_T, alt_S)
     S_diag => diag_CS%S
   endif
 
+  ! Defaults here are based on wanting to update intensive quantities frequently as soon as the model state changes.
+  ! Conversely, for extensive quantities, in an effort to close budgets and to be consistent with the total time
+  ! tendency, we construct the diagnostic grid at the beginning of the baroclinic timestep and remap all extensive
+  ! quantities to the same grid
+  update_intensive_local = .true.
+  if (present(update_intensive)) update_intensive_local = update_intensive
+  update_extensive_local = .false.
+  if (present(update_extensive)) update_extensive_local = update_extensive
+
   if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates)
 
   if (diag_cs%diag_grid_overridden) then
      call MOM_error(FATAL, "diag_update_remap_grids was called, but current grids in "// &
                            "diagnostic structure have been overridden")
   endif
 
-  do i=1, diag_cs%num_diag_coords
-    call diag_remap_update(diag_cs%diag_remap_cs(i), &
-                           diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, &
-                           diag_cs%eqn_of_state)
-  enddo
+  if (update_intensive_local) then
+    do i=1, diag_cs%num_diag_coords
+      call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, &
+                             diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h)
+    enddo
+  endif
+  if (update_extensive_local) then
+    do i=1, diag_cs%num_diag_coords
+      call diag_remap_update(diag_cs%diag_remap_cs(i), diag_cs%G, diag_cs%GV, diag_cs%US, h_diag, T_diag, S_diag, &
+                             diag_cs%eqn_of_state, diag_cs%diag_remap_cs(i)%h_extensive)
+    enddo
+  endif
 
 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
   ! Keep a copy of H - used to check whether grids are up-to-date

src/framework/MOM_diag_remap.F90

@@ -112,6 +112,7 @@ module MOM_diag_remap
   type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordiantes for this axes
   integer :: nz = 0 !< Number of vertical levels used for remapping
   real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses
+  real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive variables
   real, dimension(:), allocatable :: dz !< Nominal layer thicknesses
   integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces
   integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers
@@ -271,15 +272,16 @@ function diag_remap_axes_configured(remap_cs)
 !! height or layer thicknesses changes. In the case of density-based
 !! coordinates then technically we should also regenerate the
 !! target grid whenever T/S change.
-subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
-  type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure
-  type(ocean_grid_type),    pointer    :: G  !< The ocean's grid type
-  type(verticalGrid_type),  intent(in) :: GV !< ocean vertical grid structure
-  type(unit_scale_type),    intent(in) :: US !< A dimensional unit scaling type
-  real, dimension(:, :, :), intent(in) :: h  !< New thickness
-  real, dimension(:, :, :), intent(in) :: T  !< New T
-  real, dimension(:, :, :), intent(in) :: S  !< New S
+subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target)
+  type(diag_remap_ctrl), intent(inout)    :: remap_cs !< Diagnostic coordinate control structure
+  type(ocean_grid_type),    pointer       :: G  !< The ocean's grid type
+  type(verticalGrid_type),  intent(in   ) :: GV !< ocean vertical grid structure
+  type(unit_scale_type),    intent(in   ) :: US !< A dimensional unit scaling type
+  real, dimension(:, :, :), intent(in   ) :: h  !< Thicknesses used to construct new diagnostic grid
+  real, dimension(:, :, :), intent(in   ) :: T  !< Temperatures used to construct new diagnostic grid
+  real, dimension(:, :, :), intent(in   ) :: S  !< Salinity used to construct new diagnostic grid
   type(EOS_type),           pointer    :: eqn_of_state !< A pointer to the equation of state
+  real, dimension(:, :, :), intent(inout) :: h_target  !< Where to store the new diagnostic array
 
   ! Local variables
   real, dimension(remap_cs%nz + 1) :: zInterfaces
@@ -306,6 +308,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
     call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., &
                               answers_2018=remap_cs%answers_2018)
     allocate(remap_cs%h(G%isd:G%ied,G%jsd:G%jed, nz))
+    allocate(remap_cs%h_extensive(G%isd:G%ied,G%jsd:G%jed, nz))
     remap_cs%initialized = .true.
   endif
 
@@ -314,7 +317,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
   ! assumption that h, T, S has changed.
   do j=G%jsc-1, G%jec+1 ; do i=G%isc-1, G%iec+1
     if (G%mask2dT(i,j)==0.) then
-      remap_cs%h(i,j,:) = 0.
+      h_target(i,j,:) = 0.
       cycle
     endif
 
@@ -338,7 +341,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
 !                           US%Z_to_m*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces)
       call MOM_error(FATAL,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!")
     endif
-    remap_cs%h(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1)
+    h_target(i,j,:) = zInterfaces(1:nz) - zInterfaces(2:nz+1)
   enddo ; enddo
 
 end subroutine diag_remap_update