Merge pull request mom-ocean#931 from marshallward/diag_fix

Diagnostic fixes to speed, dudt, and dvdt

src/core/MOM.F90

@@ -2748,10 +2748,10 @@ subroutine extract_surface_state(CS, sfc_state)
       sfc_state%SST(i,j) = CS%tv%T(i,j,1)
       sfc_state%SSS(i,j) = CS%tv%S(i,j,1)
     enddo ; enddo ; endif
-    do j=js,je ; do I=IscB,IecB
+    do j=js,je ; do I=is-1,ie
       sfc_state%u(I,j) = u(I,j,1)
     enddo ; enddo
-    do J=JscB,JecB ; do i=is,ie
+    do J=js-1,je ; do i=is,ie
       sfc_state%v(i,J) = v(i,J,1)
     enddo ; enddo
 
@@ -2803,12 +2803,15 @@ subroutine extract_surface_state(CS, sfc_state)
     enddo ! end of j loop
 
 !   Determine the mean velocities in the uppermost depth_ml fluid.
+    ! NOTE: Velocity loops start on `[ij]s-1` in order to update halo values
+    !       required by the speed diagnostic on the non-symmetric grid.
+    !       This assumes that u and v halos have already been updated.
     if (CS%Hmix_UV>0.) then
       !### This calculation should work in thickness (H) units instead of Z, but that
       !### would change answers at roundoff in non-Boussinesq cases.
       depth_ml = CS%Hmix_UV
       !$OMP parallel do default(shared) private(depth,dh,hv)
-      do J=jscB,jecB
+      do J=js-1,ie
         do i=is,ie
           depth(i) = 0.0
           sfc_state%v(i,J) = 0.0
@@ -2835,11 +2838,11 @@ subroutine extract_surface_state(CS, sfc_state)
 
       !$OMP parallel do default(shared) private(depth,dh,hu)
       do j=js,je
-        do I=iscB,iecB
+        do I=is-1,ie
           depth(I) = 0.0
           sfc_state%u(I,j) = 0.0
         enddo
-        do k=1,nz ; do I=iscB,iecB
+        do k=1,nz ; do I=is-1,ie
           hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%H_to_Z
           if (depth(i) + hu < depth_ml) then
             dh = hu
@@ -2852,17 +2855,17 @@ subroutine extract_surface_state(CS, sfc_state)
           depth(I) = depth(I) + dh
         enddo ; enddo
         ! Calculate the average properties of the mixed layer depth.
-        do I=iscB,iecB
+        do I=is-1,ie
           if (depth(I) < GV%H_subroundoff*GV%H_to_Z) &
               depth(I) = GV%H_subroundoff*GV%H_to_Z
           sfc_state%u(I,j) = sfc_state%u(I,j) / depth(I)
         enddo
       enddo ! end of j loop
     else ! Hmix_UV<=0.
-      do j=js,je ; do I=IscB,IecB
+      do j=js,je ; do I=is-1,ie
         sfc_state%u(I,j) = u(I,j,1)
       enddo ; enddo
-      do J=JscB,JecB ; do i=is,ie
+      do J=js-1,je ; do i=is,ie
         sfc_state%v(i,J) = v(i,J,1)
       enddo ; enddo
     endif

src/diagnostics/MOM_diagnostics.F90

@@ -1076,16 +1076,21 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, CS)
 end subroutine calculate_energy_diagnostics
 
 !> This subroutine registers fields to calculate a diagnostic time derivative.
-subroutine register_time_deriv(f_ptr, deriv_ptr, CS)
-  real, dimension(:,:,:), target :: f_ptr     !< Field whose derivative is taken.
-  real, dimension(:,:,:), target :: deriv_ptr !< Field in which the calculated time derivatives
-                                              !! will be placed.
+subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
+  integer, intent(in), dimension(3) :: lb     !< Lower index bound of f_ptr
+  real, dimension(lb(1):,lb(2):,:), target :: f_ptr
+                                              !< Time derivative operand
+  real, dimension(lb(1):,lb(2):,:), target :: deriv_ptr
+                                              !< Time derivative of f_ptr
   type(diagnostics_CS),  pointer :: CS        !< Control structure returned by previous call to
                                               !! diagnostics_init.
 
   ! This subroutine registers fields to calculate a diagnostic time derivative.
+  ! NOTE: Lower bound is required for grid indexing in calculate_derivs().
+  !       We assume that the vertical axis is 1-indexed.
 
-  integer :: m
+  integer :: m      !< New index of deriv_ptr in CS%deriv
+  integer :: ub(3)  !< Upper index bound of f_ptr, based on shape.
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "register_time_deriv: Module must be initialized before it is used.")
@@ -1098,9 +1103,11 @@ subroutine register_time_deriv(f_ptr, deriv_ptr, CS)
 
   m = CS%num_time_deriv+1 ; CS%num_time_deriv = m
 
-  CS%nlay(m) = size(f_ptr(:,:,:),3)
+  ub(:) = lb(:) + shape(f_ptr) - 1
+
+  CS%nlay(m) = size(f_ptr, 3)
   CS%deriv(m)%p => deriv_ptr
-  allocate(CS%prev_val(m)%p(size(f_ptr(:,:,:),1), size(f_ptr(:,:,:),2), CS%nlay(m)) )
+  allocate(CS%prev_val(m)%p(lb(1):ub(1), lb(2):ub(2), CS%nlay(m)))
 
   CS%var_ptr(m)%p => f_ptr
   CS%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)
@@ -1551,21 +1558,21 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       'Zonal Acceleration', 'm s-2')
   if ((CS%id_du_dt>0) .and. .not.associated(CS%du_dt)) then
     call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
-    call register_time_deriv(MIS%u, CS%du_dt, CS)
+    call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
   endif
 
   CS%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, Time, &
       'Meridional Acceleration', 'm s-2')
   if ((CS%id_dv_dt>0) .and. .not.associated(CS%dv_dt)) then
     call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
-    call register_time_deriv(MIS%v, CS%dv_dt, CS)
+    call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
   endif
 
   CS%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, Time, &
       'Thickness tendency', trim(thickness_units)//" s-1", v_extensive = .true.)
   if ((CS%id_dh_dt>0) .and. .not.associated(CS%dh_dt)) then
     call safe_alloc_ptr(CS%dh_dt,isd,ied,jsd,jed,nz)
-    call register_time_deriv(MIS%h, CS%dh_dt, CS)
+    call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS)
   endif
 
   ! layer thickness variables
@@ -2009,15 +2016,15 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS)
   if (associated(CS%dKE_dt)) then
     if (.not.associated(CS%du_dt)) then
       call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
-      call register_time_deriv(MIS%u, CS%du_dt, CS)
+      call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
     endif
     if (.not.associated(CS%dv_dt)) then
       call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
-      call register_time_deriv(MIS%v, CS%dv_dt, CS)
+      call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
     endif
     if (.not.associated(CS%dh_dt)) then
       call safe_alloc_ptr(CS%dh_dt,isd,ied,jsd,jed,nz)
-      call register_time_deriv(MIS%h, CS%dh_dt, CS)
+      call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS)
     endif
   endif