Implemented xy-averaged diagnostics for interface data

- @nicjhan correctly disabled the xy-averaging of interface data in
  58b26a0 because it wasn't implemented.
  - Now it is.
- Interface data does not use thickness weights.
- Layer data that is extensive also does not use thickness weights.

src/framework/MOM_diag_mediator.F90

@@ -359,7 +359,8 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical)
       call define_axes_group(diag_cs, (/ id_xh, id_yh, id_zi /), diag_cs%remap_axesTi(i), &
            nz=nz, vertical_coordinate_number=i, &
            x_cell_method='mean', y_cell_method='mean', v_cell_method='point', &
-           is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true.)
+           is_h_point=.true., is_interface=.true., is_native=.false., needs_interpolating=.true., &
+           xyave_axes=diag_cs%remap_axesZi(i))
 
       !! \note Remapping for B points is not yet implemented so needs_remapping is not provided for remap_axesBi
       call define_axes_group(diag_cs, (/ id_xq, id_yq, id_zi /), diag_cs%remap_axesBi(i), &
@@ -939,6 +940,7 @@ subroutine post_xy_average(diag_cs, diag, field)
   if (diag_cs%ave_enabled) then
     call horizontally_average_diag_field(diag_cs%diag_remap_cs(diag%axes%vertical_coordinate_number), &
                                          diag_cs%G, staggered_in_x, staggered_in_y, &
+                                         diag%axes%is_layer, diag%v_extensive, &
                                          diag_cs%missing_value, field, averaged_field)
     used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, weight=diag_cs%time_int)
   endif
@@ -1192,9 +1194,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name,
              interp_method=interp_method, tile_count=tile_count)
     call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &
                              cell_methods, v_cell_method, v_extensive=v_extensive)
-  endif
-  if (is_root_pe() .and. diag_CS%doc_unit > 0) then
-    if (associated(axes%xyave_axes)) then
+    if (is_root_pe() .and. diag_CS%doc_unit > 0) then
       msg = ''
       if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"'
       call log_available_diag(fms_xyave_id>0, module_name, trim(field_name)//'_xyave', cm_string, &
@@ -1254,9 +1254,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name,
                err_msg=err_msg, interp_method=interp_method, tile_count=tile_count)
       call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, &
                                cell_methods, v_cell_method, v_extensive=v_extensive)
-    endif
-    if (is_root_pe() .and. diag_CS%doc_unit > 0) then
-      if (associated(axes%xyave_axes)) then
+      if (is_root_pe() .and. diag_CS%doc_unit > 0) then
         msg = 'native name is "'//trim(field_name)//'_xyave"'
         call log_available_diag(fms_xyave_id>0, module_name, trim(cmor_field_name)//'_xyave', cm_string, &
                                 msg, diag_CS, posted_cmor_long_name, posted_cmor_units, &

src/framework/MOM_diag_remap.F90

@@ -512,67 +512,143 @@ end subroutine vertically_interpolate_diag_field
 
 !> Horizontally average field
 subroutine horizontally_average_diag_field(remap_cs, G, staggered_in_x, staggered_in_y, &
+                                           is_layer, is_extensive, &
                                            missing_value, field, averaged_field)
   type(diag_remap_ctrl),  intent(in) :: remap_cs !< Diagnostic coodinate control structure
   type(ocean_grid_type),  intent(in) :: G !< Ocean grid structure
-  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
+  logical,                intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
+  logical,                intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
+  logical,                intent(in) :: is_layer !< True if the z-axis location is at h points
+  logical,                intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
   real,                   intent(in) :: missing_value !< A missing_value to assign land/vanished points
   real, dimension(:,:,:), intent(in) :: field !<  The diagnostic field to be remapped
   real, dimension(:),  intent(inout) :: averaged_field !< Field argument horizontally averaged
   ! Local variables
-  real, dimension(remap_cs%nz) :: vol_sum, stuff_sum
-  real :: v1, v2, total_volume, total_stuff
+  real, dimension(remap_cs%nz+1) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages
+  real :: v1, v2, total_volume, total_stuff, val
   integer :: i, j, k
 
   call assert(remap_cs%initialized, 'horizontally_average_diag_field: remap_cs not initialized.')
-  call assert(size(field, 3) == remap_cs%nz, &
+  if (is_layer) then
+    call assert(size(field, 3) == remap_cs%nz, &
+              'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.')
+    call assert(size(averaged_field, 1) == remap_cs%nz, &
+              'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.')
+  else
+    call assert(size(field, 3) == remap_cs%nz+1, &
               'horizontally_average_diag_field: Field vertical dimension does not match diag remap type.')
-  call assert(size(averaged_field, 1) == remap_cs%nz, &
+    call assert(size(averaged_field, 1) == remap_cs%nz+1, &
               'horizontally_average_diag_field: Averaged field dimension does not match diag remap type.')
+  endif
 
   if (staggered_in_x .and. .not. staggered_in_y) then
-    ! U-points
-    do k=1,remap_cs%nz
-      vol_sum(k) = 0.
-      stuff_sum(k) = 0.
-      do j=G%jsc, G%jec ; do i=G%isc, G%iec
-        v1 = G%areaCu(I,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i+1,j,k) )
-        v2 = G%areaCu(I-1,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i-1,j,k) )
-        vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
-        stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j)
-      enddo ; enddo
-    enddo
+    if (is_layer) then
+      ! U-points
+      do k=1,remap_cs%nz
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        if (is_extensive) then
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            v1 = G%areaCu(I,j)
+            v2 = G%areaCu(I-1,j)
+            vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+            stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j)
+          enddo ; enddo
+        else ! Intensive
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            v1 = G%areaCu(I,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i+1,j,k) )
+            v2 = G%areaCu(I-1,j) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i-1,j,k) )
+            vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+            stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j)
+          enddo ; enddo
+        endif
+      enddo
+    else ! Interface
+      do k=1,remap_cs%nz+1
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        do j=G%jsc, G%jec ; do i=G%isc, G%iec
+          v1 = G%areaCu(I,j)
+          v2 = G%areaCu(I-1,j)
+          vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+          stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(I,j,k) + v2 * field(I-1,j,k) ) * G%mask2dT(i,j)
+        enddo ; enddo
+      enddo
+    endif
   elseif (staggered_in_y .and. .not. staggered_in_x) then
-    ! V-points
-    do k=1,remap_cs%nz
-      vol_sum(k) = 0.
-      stuff_sum(k) = 0.
-      do j=G%jsc, G%jec ; do i=G%isc, G%iec
-        v1 = G%areaCv(i,J) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j+1,k) )
-        v2 = G%areaCv(i,J-1) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j-1,k) )
-        vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
-        stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j)
-      enddo ; enddo
-    enddo
+    if (is_layer) then
+      ! V-points
+      do k=1,remap_cs%nz
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        if (is_extensive) then
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            v1 = G%areaCv(i,J)
+            v2 = G%areaCv(i,J-1)
+            vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+            stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j)
+          enddo ; enddo
+        else ! Intensive
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            v1 = G%areaCv(i,J) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j+1,k) )
+            v2 = G%areaCv(i,J-1) * 0.5 * ( remap_cs%h(i,j,k) + remap_cs%h(i,j-1,k) )
+            vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+            stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j)
+          enddo ; enddo
+        endif
+      enddo
+    else ! Interface
+      do k=1,remap_cs%nz+1
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        do j=G%jsc, G%jec ; do i=G%isc, G%iec
+          v1 = G%areaCv(i,J)
+          v2 = G%areaCv(i,J-1)
+          vol_sum(k) = vol_sum(k) + 0.5 * ( v1 + v2 ) * G%mask2dT(i,j)
+          stuff_sum(k) = stuff_sum(k) + 0.5 * ( v1 * field(i,J,k) + v2 * field(i,J-1,k) ) * G%mask2dT(i,j)
+        enddo ; enddo
+      enddo
+    endif
   elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
-    ! H-points
-    do k=1,remap_cs%nz
-      vol_sum(k) = 0.
-      stuff_sum(k) = 0.
-      do j=G%jsc, G%jec ; do i=G%isc, G%iec
-        if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then
-          v1 = G%areaT(i,j) * remap_cs%h(i,j,k)
-          vol_sum(k) = vol_sum(k) + v1
-          stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k)
- if (abs(field(i,j,k))>1.e20) then
-  stop 'xxx'
-endif
+    if (is_layer) then
+      ! H-points
+      do k=1,remap_cs%nz
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        if (is_extensive) then
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then
+              v1 = G%areaT(i,j)
+              vol_sum(k) = vol_sum(k) + v1
+              stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k)
+            endif
+          enddo ; enddo
+        else ! Intensive
+          do j=G%jsc, G%jec ; do i=G%isc, G%iec
+            if (G%mask2dT(i,j)>0. .and. remap_cs%h(i,j,k)>0.) then
+              v1 = G%areaT(i,j) * remap_cs%h(i,j,k)
+              vol_sum(k) = vol_sum(k) + v1
+              stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k)
+            endif
+          enddo ; enddo
         endif
-      enddo ; enddo
-    enddo
+      enddo
+    else ! Interface
+      do k=1,remap_cs%nz+1
+        vol_sum(k) = 0.
+        stuff_sum(k) = 0.
+        do j=G%jsc, G%jec ; do i=G%isc, G%iec
+          val = field(i,j,k)
+          if (G%mask2dT(i,j)>0. .and. val/=missing_value) then
+            v1 = G%areaT(i,j)
+            vol_sum(k) = vol_sum(k) + v1
+            stuff_sum(k) = stuff_sum(k) + v1 * field(i,j,k)
+          endif
+        enddo ; enddo
+      enddo
+    endif
   else
-    call assert(.false., 'vertically_interpolate_diag_field: Q point averaging is not coded yet.')
+    call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.')
   endif
 
   call sum_across_PEs(vol_sum, remap_cs%nz)