Merge pull request #21 from NOAA-GFDL/dev/gfdl

merge in latest dev/gfdl updates

src/diagnostics/MOM_diagnostics.F90

@@ -1128,8 +1128,17 @@ subroutine calculate_derivs(dt, G, CS)
   if (dt > 0.0) then ; Idt = 1.0/dt
   else ; return ; endif
 
+  ! Because the field is unknown, its grid index bounds are also unknown.
+  ! Additionally, two of the fields (dudt, dvdt) require calculation of spatial
+  ! derivatives when computing d(KE)/dt.  This raises issues in non-symmetric
+  ! mode, where the symmetric boundaries (west, south) may not be updated.
+
+  ! For this reason, we explicitly loop from isc-1:iec and jsc-1:jec, in order
+  ! to force boundary value updates, even though it may not be strictly valid
+  ! for all fields.  Note this assumes a halo, and that it has been updated.
+
   do m=1,CS%num_time_deriv
-    do k=1,CS%nlay(m) ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
+    do k=1,CS%nlay(m) ; do j=G%jsc-1,G%jec ; do i=G%isc-1,G%iec
       CS%deriv(m)%p(i,j,k) = (CS%var_ptr(m)%p(i,j,k) - CS%prev_val(m)%p(i,j,k)) * Idt
       CS%prev_val(m)%p(i,j,k) = CS%var_ptr(m)%p(i,j,k)
     enddo ; enddo ; enddo

src/framework/MOM_coms.F90

@@ -535,7 +535,7 @@ subroutine carry_overflow(int_sum, prec_error)
   ! This subroutine handles carrying of the overflow.
   integer :: i, num_carry
 
-  do i=ni,2,-1 ; if (abs(int_sum(i)) > prec) then
+  do i=ni,2,-1 ; if (abs(int_sum(i)) >= prec) then
     num_carry = int(int_sum(i) * I_prec)
     int_sum(i) = int_sum(i) - num_carry*prec
     int_sum(i-1) = int_sum(i-1) + num_carry
@@ -559,7 +559,7 @@ subroutine regularize_ints(int_sum)
   logical :: positive
   integer :: i, num_carry
 
-  do i=ni,2,-1 ; if (abs(int_sum(i)) > prec) then
+  do i=ni,2,-1 ; if (abs(int_sum(i)) >= prec) then
     num_carry = int(int_sum(i) * I_prec)
     int_sum(i) = int_sum(i) - num_carry*prec
     int_sum(i-1) = int_sum(i-1) + num_carry

src/framework/MOM_diag_mediator.F90

@@ -750,7 +750,7 @@ subroutine set_masks_for_axes(G, diag_cs)
       call assert(axes%nz == nk, 'set_masks_for_axes: vertical size mismatch at h-interfaces')
       call assert(.not. associated(axes%mask3d), 'set_masks_for_axes: already associated')
       allocate( axes%mask3d(G%isd:G%ied,G%jsd:G%jed,nk+1) ) ; axes%mask3d(:,:,:) = 0.
-      do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
+      do J=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
         if (h_axes%mask3d(i,j,1) > 0.) axes%mask3d(i,J,1) = 1.
         do K = 2, nk
           if (h_axes%mask3d(i,j,k-1) + h_axes%mask3d(i,j,k) > 0.) axes%mask3d(i,J,k) = 1.
@@ -787,7 +787,6 @@ subroutine set_masks_for_axes(G, diag_cs)
         if (h_axes%mask3d(i,j,k) + h_axes%mask3d(i+1,j+1,k) + &
             h_axes%mask3d(i+1,j,k) + h_axes%mask3d(i,j+1,k) > 0.) axes%mask3d(I,J,k) = 1.
       enddo ; enddo ; enddo
-
     endif
   enddo
 
@@ -1706,6 +1705,11 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask)
       endif
     endif
   endif
+
+  if (diag%fms_xyave_diag_id>0) then
+    call post_xy_average(diag_cs, diag, locfield)
+  endif
+
   if ((diag%conversion_factor /= 0.) .and. (diag%conversion_factor /= 1.) .and. dl<2) &
     deallocate( locfield )
 
@@ -1718,6 +1722,7 @@ subroutine post_xy_average(diag_cs, diag, field)
   type(diag_ctrl),   intent(in) :: diag_cs !< Diagnostics mediator control structure
   ! Local variable
   real, dimension(size(field,3)) :: averaged_field
+  logical, dimension(size(field,3)) :: averaged_mask
   logical :: staggered_in_x, staggered_in_y, used
   integer :: nz, remap_nz, coord
 
@@ -1732,7 +1737,8 @@ subroutine post_xy_average(diag_cs, diag, field)
     call horizontally_average_diag_field(diag_cs%G, diag_cs%h, &
                                          staggered_in_x, staggered_in_y, &
                                          diag%axes%is_layer, diag%v_extensive, &
-                                         diag_cs%missing_value, field, averaged_field)
+                                         diag_cs%missing_value, field, &
+                                         averaged_field, averaged_mask)
   else
     nz = size(field, 3)
     coord = diag%axes%vertical_coordinate_number
@@ -1749,11 +1755,12 @@ subroutine post_xy_average(diag_cs, diag, field)
     call horizontally_average_diag_field(diag_cs%G, diag_cs%diag_remap_cs(coord)%h, &
                                          staggered_in_x, staggered_in_y, &
                                          diag%axes%is_layer, diag%v_extensive, &
-                                         diag_cs%missing_value, field, averaged_field)
+                                         diag_cs%missing_value, field, &
+                                         averaged_field, averaged_mask)
   endif
 
   used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, &
-                   weight=diag_cs%time_int)
+                   weight=diag_cs%time_int, mask=averaged_mask)
 end subroutine post_xy_average
 
 !> This subroutine enables the accumulation of time averages over the specified time interval.