+*Non-Boussinesq revision of MOM_vertvisc.F90

  Revised the code in MOM_vertvisc.F90 to work internally with thickness-based
units for the viscous coupling between layers, which eliminates any dependence
on the value of GV%Rho_0 in non-Boussinesq mode.  Various kinematic viscosities
are replaced with dynamic viscosities in non-Boussinesq configurations,
including revising the scaled units of the viscosities to [H Z T-1 ~> m2 s-1 or
Pa s].  This commit also modifies the code to explicitly use vertical distances
rather than thicknesses when calculating the vertical viscous coupling
coefficients in vertvisc_coef and find_coupling_coef.

  This commit changes the units of numerous variables to use thickness, vertical
distance, dynamic viscosity or other related units, including:

 - 14 elements in the vertvisc_CS type
 - 6 arguments to the private routine find_coupling_coef
 - 2 arguments to the private routine find_coupling_coef_gl90
 - 1 arguments ("a") to the public routine write_u_accel
 - 1 arguments ("a") to the public routine write_v_accel
 - 1 internal variable in vertvisc
 - 1 internal variable in vertvisc_remnant
 - 23 internal variables in vertvisc_coef,
 - 7+4+4+3 internal variables in find_coupling_coef,
 - 1 internal variable in find_coupling_coef_gl90

  Local variables that are no longer needed were eliminated in vertvisc and
vertvisc_remnant, while 2 new local variables were added to find_coupling_coef
and 6 new local variables were added to vertvisc_coef. In 6 places the
Boussinesq reference density was replaced with GV%H_to_RZ, which is equivalent
to the reference density in Boussinesq mode, but scales to 1 in non-Boussinesq
mode.

  The previous dimensional rescaling factor for KD_GL90 was incorrect (and
inconsistent with the correct scaling factor used when reading in the analagous
spatially varying kappa_gl90_2d); this has been corrected in this commit.

  A total of 59 GV%H_to_Z or GV%Z_to_H unit conversion factors or references to
GV%Rho_0 were eliminated with these changes, and in non-Boussinesq mode there is
no longer any dependence on the Boussinesq reference density.

  Replaced the forces argument to find_coupling_coef with an array of the
friction velocities and use find_ustar to set them.  When in non-Boussinesq
mode, this has the effect of using forces%tau_mag and tv%SpV_avg instead of
forces%ustar and GV%Rho0 when interpolating the friction velocity and stress
magnitude in find_coupling_coef.

  Revised units used to set the GL90 viscosities and rescale to convert
diagnostics of the vertical viscosities in the MOM_vert_friction module, so that
they do not depend on RHO_0 when in non-Boussinesq mode.

  To accomodate this change in vertvisc, the units of the "a" arguments to
write_u_accel and write_v_accel were also changed to use thickness-based
arguments.

  Because GV%Z_to_H is an exact power of 2 in Boussinesq mode, all answers are
bitwise identical in that mode.  In non-Boussinesq mode, the answers are changed
by the replacement of the Boussinesq reference density by expressions using the
layer-averaged specific volumes.  This commit changes the units of 2 arguments
in public (diagnostic) subroutine interfaces.

src/diagnostics/MOM_PointAccel.F90

@@ -83,7 +83,8 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
   real,                        intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1]
   real, optional,              intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
-                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1].
+                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc
+                                                 !! [H T-1 ~> m s-1 or Pa s m-1]
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
                      optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                  !! from vertvisc [H ~> m or kg m-2].
@@ -223,8 +224,8 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
                                       (vel_scale*ADp%du_other(I,j,k)) ; enddo
     endif
     if (present(a)) then
-      write(file,'(/,"a:     ",ES10.3," ")', advance='no') US%Z_to_m*a(I,j,ks)*dt
-      do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,K)*dt) ; enddo
+      write(file,'(/,"a:     ",ES10.3," ")', advance='no') h_scale*a(I,j,ks)*dt
+      do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (h_scale*a(I,j,K)*dt) ; enddo
     endif
     if (present(hv)) then
       write(file,'(/,"hvel:  ")', advance='no')
@@ -422,7 +423,8 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
   real,                        intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1]
   real, optional,              intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
-                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1].
+                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc
+                                                 !! [H T-1 ~> m s-1 or Pa s m-1]
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
                      optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                  !! from vertvisc [H ~> m or kg m-2].
@@ -566,8 +568,8 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
                                       (vel_scale*ADp%dv_other(i,J,k)) ; enddo
     endif
     if (present(a)) then
-      write(file,'(/,"a:     ",ES10.3," ")', advance='no') US%Z_to_m*a(i,J,ks)*dt
-      do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,J,K)*dt) ; enddo
+      write(file,'(/,"a:     ",ES10.3," ")', advance='no') h_scale*a(i,J,ks)*dt
+      do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (h_scale*a(i,J,K)*dt) ; enddo
     endif
     if (present(hv)) then
       write(file,'(/,"hvel:  ")', advance='no')