Eliminate unneeded diagnostic arrays in CorAdCalc

  Eliminated 4 unnecessary 3-d allocatable arrays and 8 2-d diagnostic arrays
in CorAdCalc, and simplified the code calculating these diagnostics by using the
post_product_[uv] and post_product_sum_[uv] routines.  Also grouped the calls
that allocate the memory that is still needed for diagnostics.

  This commit also includes a few other minor changes to clean up the
documentation of variable intents and unit documentation in a handful of other
places in the code:

- Add intent declarations to the arguments to chksum2d() and chksum3d()

- Corrected incorrect scale arguments for two (untested) checksum calls.

- Corrected the documented units in several comments.

  All answers and output are bitwise identical.

src/framework/MOM_checksums.F90

@@ -1944,8 +1944,8 @@ end subroutine chksum1d
 !> chksum2d does a checksum of all data in a 2-d array.
 subroutine chksum2d(array, mesg)
 
-  real, dimension(:,:) :: array !< The array to be checksummed
-  character(len=*) :: mesg  !< An identifying message
+  real, dimension(:,:), intent(in) :: array !< The array to be checksummed
+  character(len=*),     intent(in) :: mesg  !< An identifying message
 
   integer :: xs,xe,ys,ye,i,j,sum1,bc
   real :: sum
@@ -1972,8 +1972,8 @@ end subroutine chksum2d
 !> chksum3d does a checksum of all data in a 2-d array.
 subroutine chksum3d(array, mesg)
 
-  real, dimension(:,:,:) :: array !< The array to be checksummed
-  character(len=*) :: mesg  !< An identifying message
+  real, dimension(:,:,:), intent(in) :: array !< The array to be checksummed
+  character(len=*),       intent(in) :: mesg  !< An identifying message
 
   integer :: xs,xe,ys,ye,zs,ze,i,j,k, bc,sum1
   real :: sum

src/initialization/MOM_state_initialization.F90

@@ -442,7 +442,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
 
   if (new_sim) call pass_vector(u, v, G%Domain)
   if (debug .and. new_sim) then
-    call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1, scale=US%m_s_to_L_T)
+    call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1, scale=US%L_T_to_m_s)
   endif
 
   ! Optionally convert the thicknesses from m to kg m-2.  This is particularly

src/parameterizations/vertical/MOM_CVMix_conv.F90

@@ -275,7 +275,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
     ! if (CS%id_kd_conv > 0) &
     !   call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     ! if (CS%id_kv_conv > 0) &
-    !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T)
+    !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
     call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
   endif

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -1191,7 +1191,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
   real :: Uh2   ! The squared magnitude of the difference between the velocity
                 ! integrated through the mixed layer and the velocity of the
                 ! interior layer layer times the depth of the the mixed layer
-                ! [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
+                ! [H2 L2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
   real :: htot_vel  ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
   real :: hwtot     ! Sum of the thicknesses used to calculate
                     ! the near-bottom velocity magnitude [H ~> m or kg m-2].

src/user/ISOMIP_initialization.F90

@@ -42,9 +42,9 @@ module ISOMIP_initialization
 subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US)
   type(dyn_horgrid_type),          intent(in)  :: G !< The dynamic horizontal grid type
   real, dimension(G%isd:G%ied,G%jsd:G%jed), &
-                                   intent(out) :: D !< Ocean bottom depth [m ~> Z]
+                                   intent(out) :: D !< Ocean bottom depth [Z ~> m]
   type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
-  real,                            intent(in)  :: max_depth !< Maximum model depth [m ~> Z]
+  real,                            intent(in)  :: max_depth !< Maximum model depth [Z ~> m]
   type(unit_scale_type),           intent(in)  :: US !< A dimensional unit scaling type
 
   ! Local variables

src/user/Kelvin_initialization.F90

@@ -119,7 +119,7 @@ end subroutine Kelvin_OBC_end
 subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US)
   type(dyn_horgrid_type),          intent(in)  :: G !< The dynamic horizontal grid type
   real, dimension(G%isd:G%ied,G%jsd:G%jed), &
-                                   intent(out) :: D !< Ocean bottom depth [m ~> Z]
+                                   intent(out) :: D !< Ocean bottom depth [Z ~> m]
   type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
   real,                            intent(in)  :: max_depth !< Maximum model depth [Z ~> m]
   type(unit_scale_type),           intent(in)  :: US !< A dimensional unit scaling type

src/user/Phillips_initialization.F90

@@ -325,7 +325,7 @@ end function sech
 subroutine Phillips_initialize_topography(D, G, param_file, max_depth, US)
   type(dyn_horgrid_type),          intent(in)  :: G !< The dynamic horizontal grid type
   real, dimension(G%isd:G%ied,G%jsd:G%jed), &
-                                   intent(out) :: D !< Ocean bottom depth [m ~> Z]
+                                   intent(out) :: D !< Ocean bottom depth [Z ~> m]
   type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
   real,                            intent(in)  :: max_depth !< Maximum model depth [Z ~> m]
   type(unit_scale_type),           intent(in)  :: US !< A dimensional unit scaling type