From 8912b134006b7cc18c50dd4e10f4926e2afe1e93 Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Fri, 19 Feb 2021 22:47:31 -0500
Subject: [PATCH 01/18] Fix the bug that Rayleigh_[uv] was not initialized in
 MOM_barotropic

---
 src/core/MOM_barotropic.F90 | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 628f8e1c39..443905dccf 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -948,11 +948,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
   do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
     Cor_ref_u(I,j) = 0.0 ; BT_force_u(I,j) = 0.0 ; ubt(I,j) = 0.0
     Datu(I,j) = 0.0 ; bt_rem_u(I,j) = 0.0 ; uhbt0(I,j) = 0.0
+    Rayleigh_u(I,j) = 0.0
   enddo ; enddo
   !$OMP parallel do default(shared)
   do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
     Cor_ref_v(i,J) = 0.0 ; BT_force_v(i,J) = 0.0 ; vbt(i,J) = 0.0
     Datv(i,J) = 0.0 ; bt_rem_v(i,J) = 0.0 ; vhbt0(i,J) = 0.0
+    Rayleigh_v(i,J) = 0.0
   enddo ; enddo
 
   ! Copy input arrays into their wide-halo counterparts.

From 42eb7c7dfe9e08d5b2037f17e78474f80414fabc Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Wed, 24 Feb 2021 12:17:42 -0500
Subject: [PATCH 02/18] 1. Initialize Rayleigh_[uv] to separate loops 2. Split
 the loop that addes Rayleigh drag term to [uv]_accel_bt, so that the if-block
 is now outside the loop

---
 src/core/MOM_barotropic.F90 | 88 +++++++++++++++++++++++++++----------
 1 file changed, 64 insertions(+), 24 deletions(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 443905dccf..2c496c77bc 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -948,15 +948,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
   do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
     Cor_ref_u(I,j) = 0.0 ; BT_force_u(I,j) = 0.0 ; ubt(I,j) = 0.0
     Datu(I,j) = 0.0 ; bt_rem_u(I,j) = 0.0 ; uhbt0(I,j) = 0.0
-    Rayleigh_u(I,j) = 0.0
   enddo ; enddo
   !$OMP parallel do default(shared)
   do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
     Cor_ref_v(i,J) = 0.0 ; BT_force_v(i,J) = 0.0 ; vbt(i,J) = 0.0
     Datv(i,J) = 0.0 ; bt_rem_v(i,J) = 0.0 ; vhbt0(i,J) = 0.0
-    Rayleigh_v(i,J) = 0.0
   enddo ; enddo
 
+  if (CS%linear_wave_drag) then
+    !$OMP parallel do default(shared)
+    do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
+      Rayleigh_u(I,j) = 0.0
+    enddo ; enddo
+    !$OMP parallel do default(shared)
+    do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
+      Rayleigh_v(i,J) = 0.0
+    enddo ; enddo
+  endif
+
   ! Copy input arrays into their wide-halo counterparts.
   if (interp_eta_PF) then
     !$OMP parallel do default(shared)
@@ -1902,21 +1911,28 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
         endif ; enddo ; enddo
         !$OMP end do nowait
       endif
+
       !$OMP do schedule(static)
       do J=jsv-1,jev ; do i=isv-1,iev+1
         vel_prev = vbt(i,J)
         vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + &
                     dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
         vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev
+      enddo ; enddo
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv-1,iev+1
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * &
               ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
-        else
+        enddo ; enddo
+      else
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv-1,iev+1
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J))
-        endif
-      enddo ; enddo
-
+        enddo ; enddo
+      endif
+      
       if (integral_BT_cont) then
         !$OMP do schedule(static)
         do J=jsv-1,jev ; do i=isv-1,iev+1
@@ -1971,22 +1987,31 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
         endif ; enddo ; enddo
         !$OMP end do nowait
       endif
+
       !$OMP do schedule(static)
       do j=jsv,jev ; do I=isv-1,iev
         vel_prev = ubt(I,j)
         ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + &
              dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j)))
         if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0
-        ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev
+        ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev   
+      enddo ; enddo
+      !$OMP end do nowait
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do j=jsv,jev ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * &
-              ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
-        else
+             ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
+        enddo ; enddo
+        !$OMP end do nowait
+      else
+        !$OMP do schedule(static)
+        do j=jsv,jev ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j))
-        endif
-      enddo ; enddo
-      !$OMP end do nowait
+        enddo ; enddo
+        !$OMP end do nowait
+      endif
 
       if (integral_BT_cont) then
         !$OMP do schedule(static)
@@ -2048,14 +2073,20 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
              dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j)))
         if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0
         ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev
+      enddo ; enddo
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do j=jsv-1,jev+1 ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * &
               ((Cor_u(I,j) + PFu(I,j)) - ubt(I,j)*Rayleigh_u(I,j))
-        else
+        enddo ; enddo
+      else
+        !$OMP do schedule(static)
+        do j=jsv-1,jev+1 ; do I=isv-1,iev
           u_accel_bt(I,j) = u_accel_bt(I,j) + wt_accel(n) * (Cor_u(I,j) + PFu(I,j))
-        endif
-      enddo ; enddo
+        enddo ; enddo
+      endif
 
       if (integral_BT_cont) then
         !$OMP do schedule(static)
@@ -2130,15 +2161,24 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
              dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
         if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0
         vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev
+      enddo ; enddo
+      !$OMP end do nowait
 
-        if (CS%linear_wave_drag) then
+      if (CS%linear_wave_drag) then
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv,iev
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * &
-              ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
-        else
+             ((Cor_v(i,J) + PFv(i,J)) - vbt(i,J)*Rayleigh_v(i,J))
+        enddo ; enddo
+        !$OMP end do nowait
+      else
+        !$OMP do schedule(static)
+        do J=jsv-1,jev ; do i=isv,iev
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J))
-        endif
-      enddo ; enddo
-      !$OMP end do nowait
+        enddo ; enddo
+        !$OMP end do nowait
+      endif
+
       if (integral_BT_cont) then
         !$OMP do schedule(static)
         do J=jsv-1,jev ; do i=isv,iev

From a0b296e45b30224a05fa285442080c52814323b8 Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Wed, 24 Feb 2021 12:58:33 -0500
Subject: [PATCH 03/18] Remove the trailing space

---
 src/core/MOM_barotropic.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 2c496c77bc..962df7d0a2 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1994,7 +1994,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
         ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + &
              dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j)))
         if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0
-        ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev   
+        ubt_trans(I,j) = trans_wt1*ubt(I,j) + trans_wt2*vel_prev
       enddo ; enddo
       !$OMP end do nowait
 

From 9544de2de4c7f6d0a82a718a9cb93570a80592b5 Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Wed, 24 Feb 2021 13:05:01 -0500
Subject: [PATCH 04/18] Remove blank line spaces

---
 src/core/MOM_barotropic.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 962df7d0a2..545b7fabf7 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1932,7 +1932,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
           v_accel_bt(i,J) = v_accel_bt(i,J) + wt_accel(n) * (Cor_v(i,J) + PFv(i,J))
         enddo ; enddo
       endif
-      
+
       if (integral_BT_cont) then
         !$OMP do schedule(static)
         do J=jsv-1,jev ; do i=isv-1,iev+1

From fb882c9a278b5ff7b05e4c7d146b1c7215996e5a Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Thu, 25 Feb 2021 10:27:14 -0500
Subject: [PATCH 05/18] Fix a missing line to reset vbt if smaller than
 vel_underflow in MOM_barotropic

---
 src/core/MOM_barotropic.F90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 545b7fabf7..f30b86f57a 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -1916,7 +1916,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       do J=jsv-1,jev ; do i=isv-1,iev+1
         vel_prev = vbt(i,J)
         vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + &
-                    dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
+             dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J)))
+        if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0
         vbt_trans(i,J) = trans_wt1*vbt(i,J) + trans_wt2*vel_prev
       enddo ; enddo
 

From 0393546f781e8ee3f01d1728b0912712fa68d50f Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Thu, 25 Feb 2021 12:19:42 -0500
Subject: [PATCH 06/18] Depth-integrated momentum budget diagnostics

---
 src/core/MOM_barotropic.F90         | 11 ++++
 src/core/MOM_dynamics_split_RK2.F90 | 95 ++++++++++++++++++++++++++++-
 src/core/MOM_variables.F90          |  3 +
 3 files changed, 108 insertions(+), 1 deletion(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 545b7fabf7..cd3e6d8a5c 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -2674,6 +2674,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
     enddo ; enddo ; enddo
   endif
 
+  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hu))) then
+    do k=1,nz ; do j=js,je ; do I=is-1,ie
+      ADp%diag_hu(I,j,k) = BT_cont%h_u(I,j,k)
+    enddo ; enddo ; enddo
+  endif
+  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_h_v))) then
+    do k=1,nz ; do J=js-1,je ; do i=is,ie
+      ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k)
+    enddo ; enddo ; enddo
+  endif
+
   if (G%nonblocking_updates) then
     if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain)
     call complete_group_pass(CS%pass_ubta_uhbta, G%Domain)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index a89514cde4..37d2739945 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -161,14 +161,17 @@ module MOM_dynamics_split_RK2
   integer :: id_CAu    = -1, id_CAv    = -1
   ! integer :: id_hf_PFu    = -1, id_hf_PFv    = -1
   integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1
+  integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1
   ! integer :: id_hf_CAu    = -1, id_hf_CAv    = -1
   integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1
+  integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1
 
   ! Split scheme only.
   integer :: id_uav        = -1, id_vav        = -1
   integer :: id_u_BT_accel = -1, id_v_BT_accel = -1
   ! integer :: id_hf_u_BT_accel    = -1, id_hf_v_BT_accel    = -1
   integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1
+  integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1
   !>@}
 
   type(diag_ctrl), pointer       :: diag !< A structure that is used to regulate the
@@ -334,7 +337,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   real, allocatable, dimension(:,:) :: &
     hf_PFu_2d, hf_PFv_2d, & ! Depth integeral of hf_PFu, hf_PFv [L T-2 ~> m s-2].
     hf_CAu_2d, hf_CAv_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2].
-    hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
+    hf_u_BT_accel_2d, hf_v_BT_accel_2d, & ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
+    intz_PFu_2d, intz_PFv_2d, & ! Depth integeral of PFu, PFv [L2 T-2 ~> m2 s-2].
+    intz_CAu_2d, intz_CAv_2d, & ! Depth integeral of CAu, CAv [L2 T-2 ~> m2 s-2].
+    intz_u_BT_accel_2d, intz_v_BT_accel_2d ! Depth integeral of u_BT_accel, v_BT_accel
 
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
@@ -897,6 +903,25 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  enddo ; enddo ; enddo
   !  call post_data(CS%id_hf_PFv, hf_PFv, CS%diag)
   !endif
+  if (CS%id_intz_PFu_2d > 0) then
+    allocate(intz_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+    intz_PFu_2d(:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      intz_PFu_2d(I,j) = intz_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_PFu_2d, intz_PFu_2d, CS%diag)
+    deallocate(intz_PFu_2d)
+  endif
+  if (CS%id_intz_PFv_2d > 0) then
+    allocate(intz_PFv_2d(G%isd:G%ied,G%JsdB:G%JedB))
+    intz_PFv_2d(:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      intz_PFv_2d(i,J) = intz_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_PFv_2d, intz_PFv_2d, CS%diag)
+    deallocate(intz_PFv_2d)
+  endif
+
   if (CS%id_hf_PFu_2d > 0) then
     allocate(hf_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     hf_PFu_2d(:,:) = 0.0
@@ -930,6 +955,25 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  enddo ; enddo ; enddo
   !  call post_data(CS%id_hf_CAv, hf_CAv, CS%diag)
   !endif
+  if (CS%id_intz_CAu_2d > 0) then
+    allocate(intz_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+    intz_CAu_2d(:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      intz_CAu_2d(I,j) = intz_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_CAu_2d, intz_CAu_2d, CS%diag)
+    deallocate(intz_CAu_2d)
+  endif
+  if (CS%id_intz_CAv_2d > 0) then
+    allocate(intz_CAv_2d(G%isd:G%ied,G%JsdB:G%JedB))
+    intz_CAv_2d(:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      intz_CAv_2d(i,J) = intz_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_CAv_2d, intz_CAv_2d, CS%diag)
+    deallocate(intz_CAv_2d)
+  endif
+
   if (CS%id_hf_CAu_2d > 0) then
     allocate(hf_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     hf_CAu_2d(:,:) = 0.0
@@ -963,6 +1007,25 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  enddo ; enddo ; enddo
   !  call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag)
   !endif
+  if (CS%id_intz_u_BT_accel_2d > 0) then
+    allocate(intz_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+    intz_u_BT_accel_2d(:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      intz_u_BT_accel_2d(I,j) = intz_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_u_BT_accel_2d, intz_u_BT_accel_2d, CS%diag)
+    deallocate(intz_u_BT_accel_2d)
+  endif
+  if (CS%id_intz_v_BT_accel_2d > 0) then
+    allocate(intz_v_BT_accel_2d(G%isd:G%ied,G%JsdB:G%JedB))
+    intz_v_BT_accel_2d(:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      intz_v_BT_accel_2d(i,J) = intz_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_v_BT_accel_2d, intz_v_BT_accel_2d, CS%diag)
+    deallocate(intz_v_BT_accel_2d)
+  endif
+
   if (CS%id_hf_u_BT_accel_2d > 0) then
     allocate(hf_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     hf_u_BT_accel_2d(:,:) = 0.0
@@ -1388,6 +1451,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       conversion=US%L_T2_to_m_s2)
   if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
+  CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
+      'Depth-integral of Zonal Pressure Force Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
+      'Depth-integral of Meridional Pressure Force Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
@@ -1398,6 +1471,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       conversion=US%L_T2_to_m_s2)
   if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
+  CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
+      'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
+      'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
       'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)
   CS%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, Time, &
@@ -1428,6 +1511,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       conversion=US%L_T2_to_m_s2)
   if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
+  CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
+      'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
+
+  CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
+      'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
   id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=CLOCK_MODULE)
   id_clock_pres       = cpu_clock_id('(Ocean pressure force)',           grain=CLOCK_MODULE)
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index d81cf28e17..eb54e59d06 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -188,6 +188,9 @@ module MOM_variables
 
   real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points
   real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points
+  real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
+  real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
+  
 end type accel_diag_ptrs
 
 !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.

From 73addc4c2b05aedf26e239f716f1dacb17bd2526 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Thu, 25 Feb 2021 13:54:20 -0500
Subject: [PATCH 07/18] Depth-integrated momentum budget diagnostics

---
 src/core/MOM_CoriolisAdv.F90                  | 96 +++++++++++++++++--
 src/core/MOM_dynamics_split_RK2.F90           | 24 ++---
 src/core/MOM_variables.F90                    |  2 +-
 .../lateral/MOM_hor_visc.F90                  | 38 +++++++-
 4 files changed, 135 insertions(+), 25 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 8dab711d32..4a39d0e839 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -76,8 +76,10 @@ module MOM_CoriolisAdv
   integer :: id_rvxu = -1, id_rvxv = -1
   ! integer :: id_hf_gKEu    = -1, id_hf_gKEv    = -1
   integer :: id_hf_gKEu_2d = -1, id_hf_gKEv_2d = -1
+  integer :: id_intz_gKEu_2d = -1, id_intz_gKEv_2d = -1
   ! integer :: id_hf_rvxu    = -1, id_hf_rvxv    = -1
   integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1
+  integer :: id_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1
   !>@}
 end type CoriolisAdv_CS
 
@@ -218,7 +220,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
 ! Diagnostics for fractional thickness-weighted terms
   real, allocatable, dimension(:,:) :: &
     hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2].
-    hf_rvxu_2d, hf_rvxv_2d    ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
+    hf_rvxu_2d, hf_rvxv_2d, & ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
+    intz_gKEu_2d, intz_gKEv_2d, & ! Depth integral of gKEu, gKEv [L T-2 ~> m s-2].
+    intz_rvxu_2d, intz_rvxv_2d    ! Depth integral of rvxu, rvxv [L T-2 ~> m s-2].
   !real, allocatable, dimension(:,:,:) :: &
   !  hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness  [L T-2 ~> m s-2].
   !  hf_rvxu, hf_rvxv    ! accel. due to RV x fract. thickness [L T-2 ~> m s-2].
@@ -883,6 +887,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
       deallocate(hf_gKEv_2d)
     endif
 
+    if (CS%id_intz_gKEu_2d > 0) then
+      allocate(intz_gKEu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+      intz_gKEu_2d(:,:) = 0.0
+      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+        intz_gKEu_2d(I,j) = intz_gKEu_2d(I,j) + AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_intz_gKEu_2d, intz_gKEu_2d, CS%diag)
+      deallocate(intz_gKEu_2d)
+    endif
+
+    if (CS%id_intz_gKEv_2d > 0) then
+      allocate(intz_gKEv_2d(G%isd:G%ied,G%JsdB:G%JedB))
+      intz_gKEv_2d(:,:) = 0.0
+      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+        intz_gKEv_2d(i,J) = intz_gKEv_2d(i,J) + AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_intz_gKEv_2d, intz_gKEv_2d, CS%diag)
+      deallocate(intz_gKEv_2d)
+    endif
+
     !if (CS%id_hf_rvxv > 0) then
     !  allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
     !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
@@ -917,6 +941,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
       enddo ; enddo ; enddo
       call post_data(CS%id_hf_rvxu_2d, hf_rvxu_2d, CS%diag)
       deallocate(hf_rvxu_2d)
+    endif
+
+        if (CS%id_intz_rvxv_2d > 0) then
+      allocate(intz_rvxv_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+      intz_rvxv_2d(:,:) = 0.0
+      do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+        intz_rvxv_2d(I,j) = intz_rvxv_2d(I,j) + AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_intz_rvxv_2d, intz_rvxv_2d, CS%diag)
+      deallocate(intz_rvxv_2d)
+    endif
+
+    if (CS%id_intz_rvxu_2d > 0) then
+      allocate(intz_rvxu_2d(G%isd:G%ied,G%JsdB:G%JedB))
+      intz_rvxu_2d(:,:) = 0.0
+      do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+        intz_rvxu_2d(i,J) = intz_rvxu_2d(i,J) + AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k)
+      enddo ; enddo ; enddo
+      call post_data(CS%id_intz_rvxu_2d, intz_rvxu_2d, CS%diag)
+      deallocate(intz_rvxu_2d)
     endif
   endif
 
@@ -1163,24 +1207,24 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'Potential Vorticity', 'm-1 s-1', conversion=GV%m_to_H*US%s_to_T)
 
   CS%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, Time, &
-     'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'Zonal Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_gKEu > 0) call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, Time, &
-     'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'Meridional Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_gKEv > 0) call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
 
   CS%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, Time, &
-     'Meridional Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'Meridional Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_rvxu > 0) call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
 
   CS%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, Time, &
-     'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'Zonal Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_rvxv > 0) call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
 
   !CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, &
   !   'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
-  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
+  !   'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_gKEu > 0) then
   !  call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
   !  call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
@@ -1188,7 +1232,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   !CS%id_hf_gKEv = register_diag_field('ocean_model', 'hf_gKEv', diag%axesCvL, Time, &
   !   'Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
-  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
+  !   'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_gKEv > 0) then
   !  call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
   !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
@@ -1196,7 +1240,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, Time, &
      'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
-     'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_gKEu_2d > 0) then
     call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
@@ -1204,12 +1248,28 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_hf_gKEv_2d = register_diag_field('ocean_model', 'hf_gKEv_2d', diag%axesCv1, Time, &
      'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
-     'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
   endif
 
+  CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
+     'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
+     'm2 s-2', conversion=US%L_T_to_m_s**2)
+  if (CS%id_intz_gKEu_2d > 0) then
+    call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
+  CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
+     'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
+     'm2 s-2', conversion=US%L_T_to_m_s**2)
+  if (CS%id_intz_gKEv_2d > 0) then
+    call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+  endif
+
   !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
   !   'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
   !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
@@ -1227,7 +1287,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   !endif
 
   CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, Time, &
-     'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
+     'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
      'm-1 s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
@@ -1242,6 +1302,22 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
     call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
   endif
 
+  CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
+     'Depth-integral of Meridional Acceleration from Relative Vorticity', &
+     'm-2 s-2', conversion=US%L_T_to_m_s**2)
+  if (CS%id_intz_rvxu_2d > 0) then
+    call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+  endif
+
+  CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
+     'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
+     'm-2 s-2', conversion=US%L_T_to_m_s**2)
+  if (CS%id_intz_rvxv_2d > 0) then
+    call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
+    call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
+  endif
+
 end subroutine CoriolisAdv_init
 
 !> Destructor for coriolisadv_cs
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 37d2739945..1387e71e01 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1452,13 +1452,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
-      'Depth-integral of Zonal Pressure Force Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
-      'Depth-integral of Meridional Pressure Force Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
@@ -1472,13 +1472,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
-      'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
-      'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
@@ -1512,13 +1512,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
-      'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
-      'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm s-2', &
-      conversion=US%L_T2_to_m_s2)
+      'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
   if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90
index eb54e59d06..886ee77510 100644
--- a/src/core/MOM_variables.F90
+++ b/src/core/MOM_variables.F90
@@ -190,7 +190,7 @@ module MOM_variables
   real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points
   real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
   real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
-  
+
 end type accel_diag_ptrs
 
 !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index fc6d97040a..9f2ae0ed7a 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -188,6 +188,7 @@ module MOM_hor_visc
   integer :: id_diffu     = -1, id_diffv         = -1
   ! integer :: id_hf_diffu  = -1, id_hf_diffv      = -1
   integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
+  integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1
   integer :: id_Ah_h      = -1, id_Ah_q          = -1
   integer :: id_Kh_h      = -1, id_Kh_q          = -1
   integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
@@ -276,8 +277,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2]
     boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]
 
-  real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2]
-  real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2]
+  real, allocatable, dimension(:,:) :: hf_diffu_2d, hf_diffv_2d ! Depth sum of hf_diffu, hf_diffv [L T-2 ~> m s-2]
+  real, allocatable, dimension(:,:) :: intz_diffu_2d, intz_diffv_2d ! Integral of diffu, diffv [L2 T-2 ~> m2 s-2]
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
@@ -1662,6 +1663,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     deallocate(hf_diffv_2d)
   endif
 
+  if (present(ADp) .and. (CS%id_intz_diffu_2d > 0)) then
+    allocate(intz_diffu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+    intz_diffu_2d(:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      intz_diffu_2d(I,j) = intz_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hu(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_diffu_2d, intz_diffu_2d, CS%diag)
+    deallocate(intz_diffu_2d)
+  endif
+  if (present(ADp) .and. (CS%id_intz_diffv_2d > 0)) then
+    allocate(intz_diffv_2d(G%isd:G%ied,G%JsdB:G%JedB))
+    intz_diffv_2d(:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      intz_diffv_2d(i,J) = intz_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hv(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_intz_diffv_2d, intz_diffv_2d, CS%diag)
+    deallocate(intz_diffv_2d)
+  endif
+
 end subroutine horizontal_viscosity
 
 !> Allocates space for and calculates static variables used by horizontal_viscosity().
@@ -2378,6 +2398,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
     call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif
 
+  CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, &
+      'Depth-integral of Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
+  if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
+  endif
+
+  CS%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, &
+      'Depth-integral of Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
+      conversion=US%L_T_to_m_s**2)
+  if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
+    call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
+  endif
+
   if (CS%biharmonic) then
     CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time,    &
         'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, &

From 5f827ed1dbcda5b708f5bfb6af44336f0eff1edd Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Thu, 25 Feb 2021 14:00:31 -0500
Subject: [PATCH 08/18] typo correction

---
 src/core/MOM_barotropic.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index cd3e6d8a5c..e3741c2d12 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -2679,7 +2679,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       ADp%diag_hu(I,j,k) = BT_cont%h_u(I,j,k)
     enddo ; enddo ; enddo
   endif
-  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_h_v))) then
+  if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hv))) then
     do k=1,nz ; do J=js-1,je ; do i=is,ie
       ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k)
     enddo ; enddo ; enddo

From 77d44cbb2f62f77e0e5da76f359a4ce969205f9b Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Mon, 1 Mar 2021 11:29:31 -0500
Subject: [PATCH 09/18] Dimensional consistency test

---
 src/core/MOM_CoriolisAdv.F90 | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 4a39d0e839..58cc03ce58 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1256,7 +1256,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
      'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%L_T_to_m_s**2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEu_2d > 0) then
     call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
@@ -1264,7 +1264,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%L_T_to_m_s**2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1304,7 +1304,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Relative Vorticity', &
-     'm-2 s-2', conversion=US%L_T_to_m_s**2)
+     'm-2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1312,7 +1312,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
      'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
-     'm-2 s-2', conversion=US%L_T_to_m_s**2)
+     'm-2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxv_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)

From 2c4137b24abff6598aef3db03f6cf5d385daac76 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Mon, 1 Mar 2021 11:54:53 -0500
Subject: [PATCH 10/18] dimensions

---
 src/core/MOM_CoriolisAdv.F90 | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 58cc03ce58..f26a5c1c0f 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -221,8 +221,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
   real, allocatable, dimension(:,:) :: &
     hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2].
     hf_rvxu_2d, hf_rvxv_2d, & ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
-    intz_gKEu_2d, intz_gKEv_2d, & ! Depth integral of gKEu, gKEv [L T-2 ~> m s-2].
-    intz_rvxu_2d, intz_rvxv_2d    ! Depth integral of rvxu, rvxv [L T-2 ~> m s-2].
+    intz_gKEu_2d, intz_gKEv_2d, & ! Depth integral of gKEu, gKEv [L2 T-2 ~> m2 s-2].
+    intz_rvxu_2d, intz_rvxv_2d    ! Depth integral of rvxu, rvxv [L2 T-2 ~> m2 s-2].
   !real, allocatable, dimension(:,:,:) :: &
   !  hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness  [L T-2 ~> m s-2].
   !  hf_rvxu, hf_rvxv    ! accel. due to RV x fract. thickness [L T-2 ~> m s-2].
@@ -1256,7 +1256,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
      'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEu_2d > 0) then
     call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
@@ -1264,7 +1264,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1288,7 +1288,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, Time, &
      'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
-     'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
@@ -1296,7 +1296,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, Time, &
      'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
-     'm-1 s-2', conversion=US%L_T2_to_m_s2)
+     'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_rvxv_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
@@ -1304,7 +1304,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Relative Vorticity', &
-     'm-2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1312,7 +1312,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
      'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
-     'm-2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxv_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)

From 2165ebb6fb71d7e455e4cf254fb8d89fe746e3ee Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Mon, 1 Mar 2021 12:16:18 -0500
Subject: [PATCH 11/18] dimen.

---
 src/core/MOM_CoriolisAdv.F90                   |  8 ++++----
 src/core/MOM_dynamics_split_RK2.F90            | 12 ++++++------
 src/parameterizations/lateral/MOM_hor_visc.F90 |  4 ++--
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index f26a5c1c0f..60d4f0c00a 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1256,7 +1256,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
      'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEu_2d > 0) then
     call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
@@ -1264,7 +1264,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1304,7 +1304,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Relative Vorticity', &
-     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1312,7 +1312,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
      'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
-     'm2 s-2', conversion=US%H_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxv_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 1387e71e01..0431318384 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1453,12 +1453,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
@@ -1473,12 +1473,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
@@ -1513,12 +1513,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
       'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 9f2ae0ed7a..efe4a71bdd 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -2400,14 +2400,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
 
   CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
   endif
 
   CS%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=US%L_T_to_m_s**2)
+      conversion=US%Z_to_m*US%L_T2_to_m_s2)
   if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif

From 7ed57b019d6ecd687c7ecaf50f987c1b5fe4ed15 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Mon, 1 Mar 2021 14:38:21 -0500
Subject: [PATCH 12/18] H_to_m

---
 src/core/MOM_CoriolisAdv.F90                   |  8 ++++----
 src/core/MOM_dynamics_split_RK2.F90            | 12 ++++++------
 src/parameterizations/lateral/MOM_hor_visc.F90 |  4 ++--
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 60d4f0c00a..2be8e28b65 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1256,7 +1256,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
      'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEu_2d > 0) then
     call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
@@ -1264,7 +1264,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1304,7 +1304,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
      'Depth-integral of Meridional Acceleration from Relative Vorticity', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
     call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
@@ -1312,7 +1312,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
 
   CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
      'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
-     'm2 s-2', conversion=US%Z_to_m*US%L_T2_to_m_s2)
+     'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxv_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz)
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 0431318384..87b9d9148a 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1453,12 +1453,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
@@ -1473,12 +1473,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
@@ -1513,12 +1513,12 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
 
   CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
       'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index efe4a71bdd..8b63740687 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -2400,14 +2400,14 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
 
   CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
   endif
 
   CS%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
-      conversion=US%Z_to_m*US%L_T2_to_m_s2)
+      conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
     call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
   endif

From e52eabf0952cbd5129657c8d7055a076f4998d81 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 2 Mar 2021 11:31:41 -0500
Subject: [PATCH 13/18] Change to stack arrays for 2D diagnostics

---
 src/core/MOM_CoriolisAdv.F90 | 19 +++++++------------
 1 file changed, 7 insertions(+), 12 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 2be8e28b65..b1cd6209ab 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -220,15 +220,18 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
 ! Diagnostics for fractional thickness-weighted terms
   real, allocatable, dimension(:,:) :: &
     hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2].
-    hf_rvxu_2d, hf_rvxv_2d, & ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
-    intz_gKEu_2d, intz_gKEv_2d, & ! Depth integral of gKEu, gKEv [L2 T-2 ~> m2 s-2].
-    intz_rvxu_2d, intz_rvxv_2d    ! Depth integral of rvxu, rvxv [L2 T-2 ~> m2 s-2].
+    hf_rvxu_2d, hf_rvxv_2d    ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
+
   !real, allocatable, dimension(:,:,:) :: &
   !  hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness  [L T-2 ~> m s-2].
   !  hf_rvxu, hf_rvxv    ! accel. due to RV x fract. thickness [L T-2 ~> m s-2].
   ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
   ! The code is retained for degugging purposes in the future.
 
+! Diagnostics for depth-integrated momentum budget terms     
+  real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJB_(G)) :: intz_gKEv_2d, intz_rvxu_2d ! [L2 T-2 ~> m2 s-2].
+
 ! To work, the following fields must be set outside of the usual
 ! is to ie range before this subroutine is called:
 !   v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
@@ -888,23 +891,19 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
     endif
 
     if (CS%id_intz_gKEu_2d > 0) then
-      allocate(intz_gKEu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
       intz_gKEu_2d(:,:) = 0.0
       do k=1,nz ; do j=js,je ; do I=Isq,Ieq
         intz_gKEu_2d(I,j) = intz_gKEu_2d(I,j) + AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k)
       enddo ; enddo ; enddo
       call post_data(CS%id_intz_gKEu_2d, intz_gKEu_2d, CS%diag)
-      deallocate(intz_gKEu_2d)
     endif
 
     if (CS%id_intz_gKEv_2d > 0) then
-      allocate(intz_gKEv_2d(G%isd:G%ied,G%JsdB:G%JedB))
       intz_gKEv_2d(:,:) = 0.0
       do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
         intz_gKEv_2d(i,J) = intz_gKEv_2d(i,J) + AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k)
       enddo ; enddo ; enddo
       call post_data(CS%id_intz_gKEv_2d, intz_gKEv_2d, CS%diag)
-      deallocate(intz_gKEv_2d)
     endif
 
     !if (CS%id_hf_rvxv > 0) then
@@ -943,24 +942,20 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
       deallocate(hf_rvxu_2d)
     endif
 
-        if (CS%id_intz_rvxv_2d > 0) then
-      allocate(intz_rvxv_2d(G%IsdB:G%IedB,G%jsd:G%jed))
+    if (CS%id_intz_rvxv_2d > 0) then
       intz_rvxv_2d(:,:) = 0.0
       do k=1,nz ; do j=js,je ; do I=Isq,Ieq
         intz_rvxv_2d(I,j) = intz_rvxv_2d(I,j) + AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k)
       enddo ; enddo ; enddo
       call post_data(CS%id_intz_rvxv_2d, intz_rvxv_2d, CS%diag)
-      deallocate(intz_rvxv_2d)
     endif
 
     if (CS%id_intz_rvxu_2d > 0) then
-      allocate(intz_rvxu_2d(G%isd:G%ied,G%JsdB:G%JedB))
       intz_rvxu_2d(:,:) = 0.0
       do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
         intz_rvxu_2d(i,J) = intz_rvxu_2d(i,J) + AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k)
       enddo ; enddo ; enddo
       call post_data(CS%id_intz_rvxu_2d, intz_rvxu_2d, CS%diag)
-      deallocate(intz_rvxu_2d)
     endif
   endif
 

From b895f1e3ac914bef26f521ce5f43e35d4ad08f90 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Tue, 2 Mar 2021 12:08:55 -0500
Subject: [PATCH 14/18] Changed to stack arrays

---
 src/core/MOM_CoriolisAdv.F90                  |  4 ++--
 src/core/MOM_dynamics_split_RK2.F90           | 23 ++++++-------------
 .../lateral/MOM_hor_visc.F90                  |  7 ++----
 3 files changed, 11 insertions(+), 23 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index b1cd6209ab..e13bd492b8 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -228,8 +228,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
   ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
   ! The code is retained for degugging purposes in the future.
 
-! Diagnostics for depth-integrated momentum budget terms     
-  real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [L2 T-2 ~> m2 s-2]
+! Diagnostics for depth-integrated momentum budget terms
+  real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [L2 T-2 ~> m2 s-2].
   real, dimension(SZI_(G),SZJB_(G)) :: intz_gKEv_2d, intz_rvxu_2d ! [L2 T-2 ~> m2 s-2].
 
 ! To work, the following fields must be set outside of the usual
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 87b9d9148a..66f142eb21 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -337,10 +337,13 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   real, allocatable, dimension(:,:) :: &
     hf_PFu_2d, hf_PFv_2d, & ! Depth integeral of hf_PFu, hf_PFv [L T-2 ~> m s-2].
     hf_CAu_2d, hf_CAv_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2].
-    hf_u_BT_accel_2d, hf_v_BT_accel_2d, & ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
-    intz_PFu_2d, intz_PFv_2d, & ! Depth integeral of PFu, PFv [L2 T-2 ~> m2 s-2].
-    intz_CAu_2d, intz_CAv_2d, & ! Depth integeral of CAu, CAv [L2 T-2 ~> m2 s-2].
-    intz_u_BT_accel_2d, intz_v_BT_accel_2d ! Depth integeral of u_BT_accel, v_BT_accel
+    hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
+
+  real, dimension(SZIB_(G),SZJ_(G)) :: &
+    intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [L2 T-2 ~> m2 s-2].
+
+  real, dimension(SZI_(G),SZJB_(G)) :: &
+    intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [L2 T-2 ~> m2 s-2].
 
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
@@ -904,22 +907,18 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  call post_data(CS%id_hf_PFv, hf_PFv, CS%diag)
   !endif
   if (CS%id_intz_PFu_2d > 0) then
-    allocate(intz_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     intz_PFu_2d(:,:) = 0.0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
       intz_PFu_2d(I,j) = intz_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_PFu_2d, intz_PFu_2d, CS%diag)
-    deallocate(intz_PFu_2d)
   endif
   if (CS%id_intz_PFv_2d > 0) then
-    allocate(intz_PFv_2d(G%isd:G%ied,G%JsdB:G%JedB))
     intz_PFv_2d(:,:) = 0.0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
       intz_PFv_2d(i,J) = intz_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_PFv_2d, intz_PFv_2d, CS%diag)
-    deallocate(intz_PFv_2d)
   endif
 
   if (CS%id_hf_PFu_2d > 0) then
@@ -956,22 +955,18 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  call post_data(CS%id_hf_CAv, hf_CAv, CS%diag)
   !endif
   if (CS%id_intz_CAu_2d > 0) then
-    allocate(intz_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     intz_CAu_2d(:,:) = 0.0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
       intz_CAu_2d(I,j) = intz_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_CAu_2d, intz_CAu_2d, CS%diag)
-    deallocate(intz_CAu_2d)
   endif
   if (CS%id_intz_CAv_2d > 0) then
-    allocate(intz_CAv_2d(G%isd:G%ied,G%JsdB:G%JedB))
     intz_CAv_2d(:,:) = 0.0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
       intz_CAv_2d(i,J) = intz_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_CAv_2d, intz_CAv_2d, CS%diag)
-    deallocate(intz_CAv_2d)
   endif
 
   if (CS%id_hf_CAu_2d > 0) then
@@ -1008,22 +1003,18 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   !  call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag)
   !endif
   if (CS%id_intz_u_BT_accel_2d > 0) then
-    allocate(intz_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     intz_u_BT_accel_2d(:,:) = 0.0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
       intz_u_BT_accel_2d(I,j) = intz_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_u_BT_accel_2d, intz_u_BT_accel_2d, CS%diag)
-    deallocate(intz_u_BT_accel_2d)
   endif
   if (CS%id_intz_v_BT_accel_2d > 0) then
-    allocate(intz_v_BT_accel_2d(G%isd:G%ied,G%JsdB:G%JedB))
     intz_v_BT_accel_2d(:,:) = 0.0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
       intz_v_BT_accel_2d(i,J) = intz_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_v_BT_accel_2d, intz_v_BT_accel_2d, CS%diag)
-    deallocate(intz_v_BT_accel_2d)
   endif
 
   if (CS%id_hf_u_BT_accel_2d > 0) then
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 8b63740687..5bbc495e93 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -278,7 +278,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]
 
   real, allocatable, dimension(:,:) :: hf_diffu_2d, hf_diffv_2d ! Depth sum of hf_diffu, hf_diffv [L T-2 ~> m s-2]
-  real, allocatable, dimension(:,:) :: intz_diffu_2d, intz_diffv_2d ! Integral of diffu, diffv [L2 T-2 ~> m2 s-2]
+  real, dimension(SZIB_(G),SZJ_(G)) :: intz_diffu_2d ! Depth-integral of diffu [L2 T-2 ~> m2 s-2]
+  real, dimension(SZI_(G),SZJB_(G)) :: intz_diffv_2d ! Depth-integral of diffv [L2 T-2 ~> m2 s-2]
 
   real, dimension(SZIB_(G),SZJB_(G)) :: &
     dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
@@ -1664,22 +1665,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   endif
 
   if (present(ADp) .and. (CS%id_intz_diffu_2d > 0)) then
-    allocate(intz_diffu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
     intz_diffu_2d(:,:) = 0.0
     do k=1,nz ; do j=js,je ; do I=Isq,Ieq
       intz_diffu_2d(I,j) = intz_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hu(I,j,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_diffu_2d, intz_diffu_2d, CS%diag)
-    deallocate(intz_diffu_2d)
   endif
   if (present(ADp) .and. (CS%id_intz_diffv_2d > 0)) then
-    allocate(intz_diffv_2d(G%isd:G%ied,G%JsdB:G%JedB))
     intz_diffv_2d(:,:) = 0.0
     do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
       intz_diffv_2d(i,J) = intz_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hv(i,J,k)
     enddo ; enddo ; enddo
     call post_data(CS%id_intz_diffv_2d, intz_diffv_2d, CS%diag)
-    deallocate(intz_diffv_2d)
   endif
 
 end subroutine horizontal_viscosity

From ba7bbabf43885afd72802a670252459ade634ee6 Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 13:04:31 -0500
Subject: [PATCH 15/18] Typo corrected - Change from Jsd to JsdB in array
 allocations

---
 src/core/MOM_CoriolisAdv.F90                   | 10 +++++-----
 src/core/MOM_dynamics_split_RK2.F90            | 18 +++++++++---------
 .../vertical/MOM_vert_friction.F90             |  4 ++--
 3 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index e13bd492b8..c358373505 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1246,7 +1246,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, &
@@ -1262,7 +1262,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_gKEv_2d > 0) then
     call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
@@ -1270,7 +1270,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_rvxu > 0) then
   !  call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   !CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, &
@@ -1286,7 +1286,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm s-2', conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, Time, &
@@ -1302,7 +1302,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
      'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_rvxu_2d > 0) then
     call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz)
   endif
 
   CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, &
diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90
index 66f142eb21..ef7da5c291 100644
--- a/src/core/MOM_dynamics_split_RK2.F90
+++ b/src/core/MOM_dynamics_split_RK2.F90
@@ -1420,7 +1420,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
   !    'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
@@ -1430,7 +1430,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
@@ -1440,7 +1440,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', &
@@ -1450,7 +1450,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+  if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
@@ -1460,7 +1460,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
       'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
@@ -1470,7 +1470,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
       'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+  if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
       'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)
@@ -1490,7 +1490,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
   !    'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
   !    v_extensive=.true., conversion=US%L_T2_to_m_s2)
-  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
@@ -1500,7 +1500,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, Time, &
       'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
       conversion=US%L_T2_to_m_s2)
-  if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
   CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
       'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
@@ -1510,7 +1510,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
       'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
       conversion=GV%H_to_m*US%L_T2_to_m_s2)
-  if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)
+  if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
   id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=CLOCK_MODULE)
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index 7e2c2c6926..c7c9355f97 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -1817,7 +1817,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   !if (CS%id_hf_dv_dt_visc > 0) then
   !  call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz)
   !  call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, Time, &
@@ -1833,7 +1833,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
       conversion=US%L_T2_to_m_s2)
   if (CS%id_hf_dv_dt_visc_2d > 0) then
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
-    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &

From 075be64b521b40b34c7bc3feb57d36d3faba7acd Mon Sep 17 00:00:00 2001
From: Hemant Khatri <hemantkhatri1091@gmail.com>
Date: Fri, 5 Mar 2021 13:16:46 -0500
Subject: [PATCH 16/18] Typo correction - Change from Jsd to JsdB in array
 allocations

---
 src/core/MOM_CoriolisAdv.F90                           | 2 +-
 src/diagnostics/MOM_diagnostics.F90                    | 4 ++--
 src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +-
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index c358373505..7cbc1eb262 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -1230,7 +1230,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   !   'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
   !if (CS%id_hf_gKEv > 0) then
   !  call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
-  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, Time, &
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 47d322dfa0..bb5243892d 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -1767,7 +1767,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   !    call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
   !    call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
   !  endif
-  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   !endif
 
   CS%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, Time, &
@@ -1787,7 +1787,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
       call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
     endif
-    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
   endif
 
   ! layer thickness variables
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index 649fc725de..c96edb785c 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -2982,7 +2982,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
     CS%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, Time, &
         'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', &
         'm s-2', conversion=US%L_T2_to_m_s2)
-    if (CS%id_hf_dvdt_dia_2d > 0)  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
+    if (CS%id_hf_dvdt_dia_2d > 0)  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
 
     if ((CS%id_dudt_dia > 0) .or. (CS%id_hf_dudt_dia_2d > 0)) &
       call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz)

From 54033aab43071a667b93cf9fe1bf4089ec1b3c52 Mon Sep 17 00:00:00 2001
From: He Wang <He.Wang@noaa.gov>
Date: Mon, 8 Mar 2021 21:46:10 -0500
Subject: [PATCH 17/18] Bugfix: empty restart_input_dir for *_solo.res This
 commit fixes the issue that dirs%restart_input_dir is not initiated before
 "ocean_solo.res" and "ice_solo.res" are read.

---
 config_src/ice_solo_driver/ice_shelf_driver.F90 | 3 ++-
 config_src/solo_driver/MOM_driver.F90           | 4 ++++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/config_src/ice_solo_driver/ice_shelf_driver.F90 b/config_src/ice_solo_driver/ice_shelf_driver.F90
index bd64050a6f..959e4676d0 100644
--- a/config_src/ice_solo_driver/ice_shelf_driver.F90
+++ b/config_src/ice_solo_driver/ice_shelf_driver.F90
@@ -185,6 +185,8 @@ program Shelf_main
     endif
   endif
 
+  call Get_MOM_Input(param_file, dirs)
+
   ! Read ocean_solo restart, which can override settings from the namelist.
   if (file_exists(trim(dirs%restart_input_dir)//'ice_solo.res')) then
     call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ice_solo.res', action=READONLY_FILE)
@@ -215,7 +217,6 @@ program Shelf_main
     Start_time = real_to_time(0.0)
   endif
 
-  call Get_MOM_Input(param_file, dirs)
   ! Determining the internal unit scaling factors for this run.
   call unit_scaling_init(param_file, US)
 
diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90
index 9c222bb0bb..8edad7fa05 100644
--- a/config_src/solo_driver/MOM_driver.F90
+++ b/config_src/solo_driver/MOM_driver.F90
@@ -250,6 +250,10 @@ program MOM_main
   ! This call sets the number and affinity of threads with openMP.
   !$  call set_MOM_thread_affinity(ocean_nthreads, use_hyper_thread)
 
+  ! This call is required to initiate dirs%restart_input_dir for ocean_solo.res
+  ! The contents of dirs will be reread in initialize_MOM.
+  call get_MOM_input(dirs=dirs)
+
   ! Read ocean_solo restart, which can override settings from the namelist.
   if (file_exists(trim(dirs%restart_input_dir)//'ocean_solo.res')) then
     call open_ASCII_file(unit, trim(dirs%restart_input_dir)//'ocean_solo.res', action=READONLY_FILE)

From bb5eb2527da481826721c047182dced61fe2835d Mon Sep 17 00:00:00 2001
From: Robert Hallberg <Robert.Hallberg@noaa.gov>
Date: Tue, 9 Mar 2021 18:34:27 -0500
Subject: [PATCH 18/18] Explicitly set (1x1) io_domain as a default

  Added code to explicitly set a (1x1) io_domain when no other io_layout is
specified, complying with changing requirements for 2020 and later versions of
FMS, and following the default behavior of previous versions.  Also corrected
three module use statements in infra/FMS1 to eliminate any dependencies of code
in config_src/infra/FMS1 on code in src/framework.  These use statements
eventually point to the same place as before, but with less indirection.  This
change should facilitate later steps to compile everything in and under the
config_src/infra directories as libraries.  All answers are bitwise identical.
---
 config_src/infra/FMS1/MOM_diag_manager_infra.F90 | 4 ++--
 config_src/infra/FMS1/MOM_domain_infra.F90       | 4 +++-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/config_src/infra/FMS1/MOM_diag_manager_infra.F90 b/config_src/infra/FMS1/MOM_diag_manager_infra.F90
index 702c464814..18c80cf24c 100644
--- a/config_src/infra/FMS1/MOM_diag_manager_infra.F90
+++ b/config_src/infra/FMS1/MOM_diag_manager_infra.F90
@@ -20,9 +20,9 @@ module MOM_diag_manager_infra
 use diag_manager_mod, only : register_diag_field_fms => register_diag_field
 use diag_manager_mod, only : register_static_field_fms => register_static_field
 use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id
-use time_manager_mod, only : time_type
+use MOM_time_manager, only : time_type
 use MOM_domain_infra, only : MOM_domain_type
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
+use MOM_error_infra,  only : MOM_error => MOM_err, FATAL, WARNING
 
 implicit none ; private
 
diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90
index 86e85e60a6..fc39777a2f 100644
--- a/config_src/infra/FMS1/MOM_domain_infra.F90
+++ b/config_src/infra/FMS1/MOM_domain_infra.F90
@@ -4,7 +4,7 @@ module MOM_domain_infra
 ! This file is part of MOM6. See LICENSE.md for the license.
 
 use MOM_coms_infra,  only : PE_here, root_PE, num_PEs
-use MOM_cpu_clock,   only : cpu_clock_begin, cpu_clock_end
+use MOM_cpu_clock_infra, only : cpu_clock_begin, cpu_clock_end
 use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, WARNING, FATAL
 
 use mpp_domains_mod, only : domain2D, domain1D
@@ -1689,6 +1689,8 @@ subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, &
   if ((MD_in%io_layout(1) + MD_in%io_layout(2) > 0) .and. &
       (MD_in%layout(1)*MD_in%layout(2) > 1)) then
     call mpp_define_io_domain(mpp_domain, MD_in%io_layout)
+  else
+    call mpp_define_io_domain(mpp_domain, (/ 1, 1 /) )
   endif
 
 end subroutine clone_MD_to_d2D