Merge pull request mom-ocean#1070 from ashao/remap_v_extensive

Diagnostic remapping of vertically extensive fields

src/core/MOM.F90

@@ -660,6 +660,10 @@ subroutine step_MOM(forces_in, fluxes_in, 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

@@ -333,6 +333,9 @@ module MOM_diag_mediator
   !> Number of checksum-only diagnostics
   integer :: num_chksum_diags
 
+  real, dimension(:,:,:), allocatable :: h_begin !< Layer thicknesses at the beginning of the timestep used
+                                                 !! for remapping of extensive variables
+
 end type diag_ctrl
 
 !>@{ CPU clocks
@@ -456,6 +459,10 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical)
     ! For each possible diagnostic coordinate
     call diag_remap_configure_axes(diag_cs%diag_remap_cs(i), GV, US, param_file)
 
+    ! Allocate these arrays since the size of the diagnostic array is now known
+    allocate(diag_cs%diag_remap_cs(i)%h(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz))
+    allocate(diag_cs%diag_remap_cs(i)%h_extensive(G%isd:G%ied,G%jsd:G%jed, diag_cs%diag_remap_cs(i)%nz))
+
     ! This vertical coordinate has been configured so can be used.
     if (diag_remap_axes_configured(diag_cs%diag_remap_cs(i))) then
 
@@ -1475,14 +1482,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, &
@@ -1502,10 +1511,12 @@ 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%h_begin,                                                         &
+        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
@@ -3064,6 +3075,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
   diag_cs%S => null()
   diag_cs%eqn_of_state => null()
 
+  allocate(diag_cs%h_begin(G%isd:G%ied,G%jsd:G%jed,nz))
 #if defined(DEBUG) || defined(__DO_SAFETY_CHECKS__)
   allocate(diag_cs%h_old(G%isd:G%ied,G%jsd:G%jed,nz))
   diag_cs%h_old(:,:,:) = 0.0
@@ -3192,19 +3204,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() ! The layer thickneses for diagnostics [H ~> m or kg m-2]
   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
@@ -3223,17 +3241,35 @@ 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
+    diag_cs%h_begin(:,:,:) = diag_cs%h(:,:,:)
+    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,7 +272,7 @@ 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)
+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
@@ -280,6 +281,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
   real, dimension(:, :, :), intent(in) :: T  !< New temperatures [degC]
   real, dimension(:, :, :), intent(in) :: S  !< New salinities [ppt]
   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 ! Interface positions [H ~> m or kg m-2]
@@ -305,7 +307,6 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
     ! Initialize remapping and regridding on the first call
     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))
     remap_cs%initialized = .true.
   endif
 
@@ -314,7 +315,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
 
@@ -339,7 +340,9 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state)
 !                           GV%Z_to_H*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)
+    do k = 1,nz
+      h_target(i,j,k) = zInterfaces(k) - zInterfaces(k+1)
+    enddo
   enddo ; enddo
 
 end subroutine diag_remap_update
@@ -485,11 +488,12 @@ subroutine diag_remap_calc_hmask(remap_cs, G, mask)
 end subroutine diag_remap_calc_hmask
 
 !> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid.
-subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, &
+subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered_in_x, staggered_in_y, &
                                              mask, missing_value, field, reintegrated_field)
   type(diag_remap_ctrl),  intent(in) :: remap_cs !< Diagnostic coodinate control structure
-  type(ocean_grid_type),  intent(in) :: G !< Ocean grid structure
-  real, dimension(:,:,:), intent(in) :: h   !< The current thicknesses
+  type(ocean_grid_type),  intent(in) :: G        !< Ocean grid structure
+  real, dimension(:,:,:), intent(in) :: h        !< The thicknesses of the source grid
+  real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid
   logical,                intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
   logical,                intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
   real, dimension(:,:,:), pointer    :: mask !< A mask for the field
@@ -526,7 +530,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta
           if (mask(I,j,1) == 0.) cycle
         endif
         h_src(:) = 0.5 * (h(i_lo,j,:) + h(i_hi,j,:))
-        h_dest(:) = 0.5 * (remap_cs%h(i_lo,j,:) + remap_cs%h(i_hi,j,:))
+        h_dest(:) = 0.5 * (h_target(i_lo,j,:) + h_target(i_hi,j,:))
         call reintegrate_column(nz_src, h_src, field(I1,j,:), &
                                 nz_dest, h_dest, 0., reintegrated_field(I1,j,:))
       enddo
@@ -541,7 +545,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta
           if (mask(i,J,1) == 0.) cycle
         endif
         h_src(:) = 0.5 * (h(i,j_lo,:) + h(i,j_hi,:))
-        h_dest(:) = 0.5 * (remap_cs%h(i,j_lo,:) + remap_cs%h(i,j_hi,:))
+        h_dest(:) = 0.5 * (h_target(i,j_lo,:) + h_target(i,j_hi,:))
         call reintegrate_column(nz_src, h_src, field(i,J1,:), &
                                 nz_dest, h_dest, 0., reintegrated_field(i,J1,:))
       enddo
@@ -554,7 +558,7 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, staggered_in_x, sta
           if (mask(i,j,1) == 0.) cycle
         endif
         h_src(:) = h(i,j,:)
-        h_dest(:) = remap_cs%h(i,j,:)
+        h_dest(:) = h_target(i,j,:)
         call reintegrate_column(nz_src, h_src, field(i,j,:), &
                                 nz_dest, h_dest, 0., reintegrated_field(i,j,:))
       enddo