+Rescaled optics%opacity_band

  Rescaled the units of optics%opacity_band from [m-1] to [Z-1 ~> m-1], with the
opacity values returned from extract_optics_slice also rescaled by the same
factor, which can be offset by compensating changes to the opacity_scale
argument.  Also rescaled 4 other internal variables and documented the units on
3 more.  One uncommon parameter (SW_1ST_EXP_RATIO) listed the wrong units in its
get_param call, and this was corrected, but turned out not to have been logged
for any of the MOM6-examples test cases.  Some compensating changes were also
made in the MOM_generic_tracer module, which directly accesses the contents of
the optics_type (thereby preventing it from being opaque).  All answers and
output in the MOM6-examples test suite are bitwise identical.

src/parameterizations/vertical/MOM_bulk_mixed_layer.F90

@@ -424,7 +424,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
       eps(i,k) = 0.0 ; if (k > nkmb) eps(i,k) = GV%Angstrom_H
       T(i,k) = tv%T(i,j,k) ; S(i,k) = tv%S(i,j,k)
     enddo ; enddo
-    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacity_band, opacity_scale=GV%H_to_m)
+    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacity_band, opacity_scale=GV%H_to_Z)
 
     do k=1,nz ; do i=is,ie
       d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -1105,7 +1105,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
   ! To accommodate vanishing upper layers, we need to allow for an instantaneous
   ! distribution of forcing over some finite vertical extent. The bulk mixed layer
   ! code handles this issue properly.
-  H_limit_fluxes = max(GV%Angstrom_H, 1.E-30*GV%m_to_H)
+  H_limit_fluxes = max(GV%Angstrom_H, 1.0e-30*GV%m_to_H)
 
   ! diagnostic to see if need to create mass to avoid grounding
   if (CS%id_createdH>0) CS%createdH(:,:) = 0.
@@ -1160,7 +1160,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
     ! Nothing more is done on this j-slice if there is no buoyancy forcing.
     if (.not.associated(fluxes%sw)) cycle
 
-    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%m_to_H))
+    if (nsw>0) call extract_optics_slice(optics, j, G, GV, opacity=opacityBand, opacity_scale=(1.0/GV%Z_to_H))
 
     ! The surface forcing is contained in the fluxes type.
     ! We aggregate the thermodynamic forcing for a time step into the following:

src/parameterizations/vertical/MOM_opacity.F90

@@ -25,7 +25,7 @@ module MOM_opacity
 type, public :: optics_type
   integer :: nbands     !< The number of penetrating bands of SW radiation
 
-  real, allocatable :: opacity_band(:,:,:,:) !< SW optical depth per unit thickness [m-1]
+  real, allocatable :: opacity_band(:,:,:,:) !< SW optical depth per unit thickness [Z-1 ~> m-1]
                         !! The number of radiation bands is most rapidly varying (first) index.
 
   real, allocatable :: sw_pen_band(:,:,:) !< shortwave radiation [Q R Z T-1 ~> W m-2]
@@ -55,15 +55,15 @@ module MOM_opacity
                              !! water properties into the opacity (i.e., the e-folding depth) and
                              !! (perhaps) the number of bands of penetrating shortwave radiation to use.
   real :: pen_sw_scale       !<   The vertical absorption e-folding depth of the
-                             !! penetrating shortwave radiation [m].
+                             !! penetrating shortwave radiation [Z ~> m].
   real :: pen_sw_scale_2nd   !<   The vertical absorption e-folding depth of the
-                             !! (2nd) penetrating shortwave radiation [m].
-  real :: SW_1ST_EXP_RATIO   !< Ratio for 1st exp decay in Two Exp decay opacity
+                             !! (2nd) penetrating shortwave radiation [Z ~> m].
+  real :: SW_1ST_EXP_RATIO   !< Ratio for 1st exp decay in Two Exp decay opacity [nondim]
   real :: pen_sw_frac        !<   The fraction of shortwave radiation that is
-                             !! penetrating with a constant e-folding approach.
+                             !! penetrating with a constant e-folding approach [nondim]
   real :: blue_frac          !<   The fraction of the penetrating shortwave
                              !! radiation that is in the blue band [nondim].
-  real :: opacity_land_value !< The value to use for opacity over land [m-1].
+  real :: opacity_land_value !< The value to use for opacity over land [Z-1 ~> m-1].
                              !! The default is 10 m-1 - a value for muddy water.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                              !! regulate the timing of diagnostic output.
@@ -107,15 +107,15 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_
 
   ! Local variables
   integer :: i, j, k, n, is, ie, js, je, nz
-  real :: inv_sw_pen_scale  ! The inverse of the e-folding scale [m-1].
+  real :: inv_sw_pen_scale  ! The inverse of the e-folding scale [Z-1 ~> m-1].
   real :: Inv_nbands        ! The inverse of the number of bands of penetrating
                             ! shortwave radiation [nondim]
   logical :: call_for_surface  ! if horizontal slice is the surface layer
-  real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array.
+  real :: tmp(SZI_(G),SZJ_(G),SZK_(GV)) ! A 3-d temporary array for diagnosing opacity [Z-1 ~> m-1]
   real :: chl(SZI_(G),SZJ_(G),SZK_(GV)) ! The concentration of chlorophyll-A [mg m-3].
   real :: Pen_SW_tot(SZI_(G),SZJ_(G))   ! The penetrating shortwave radiation
                                         ! summed across all bands [Q R Z T-1 ~> W m-2].
-  real :: op_diag_len       ! A tiny lengthscale [m] used to remap opacity
+  real :: op_diag_len       ! A tiny lengthscale [Z ~> m] used to remap diagnostics of opacity
                             ! from op to 1/op_diag_len * tanh(op * op_diag_len)
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
@@ -128,14 +128,14 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_
     else ; Inv_nbands = 1.0 / real(optics%nbands) ; endif
 
     ! Make sure there is no division by 0.
-    inv_sw_pen_scale = 1.0 / max(CS%pen_sw_scale, 0.1*GV%Angstrom_m, &
-                                 GV%H_to_m*GV%H_subroundoff)
+    inv_sw_pen_scale = 1.0 / max(CS%pen_sw_scale, 0.1*GV%Angstrom_Z, &
+                                 GV%H_to_Z*GV%H_subroundoff)
     if ( CS%Opacity_scheme == DOUBLE_EXP ) then
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js,je ; do i=is,ie
         optics%opacity_band(1,i,j,k) = inv_sw_pen_scale
         optics%opacity_band(2,i,j,k) = 1.0 / max(CS%pen_sw_scale_2nd, &
-             0.1*GV%Angstrom_m,GV%H_to_m*GV%H_subroundoff)
+             0.1*GV%Angstrom_Z, GV%H_to_Z*GV%H_subroundoff)
       enddo ; enddo ; enddo
       if (.not.associated(sw_total) .or. (CS%pen_SW_scale <= 0.0)) then
         !$OMP parallel do default(shared)
@@ -199,7 +199,7 @@ subroutine set_opacity(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir, sw_
       call post_data(CS%id_sw_vis_pen, Pen_SW_tot, CS%diag)
     endif
     do n=1,optics%nbands ; if (CS%id_opacity(n) > 0) then
-      op_diag_len = 1e-10 ! A dimensional depth to constrain the range of opacity [m]
+      op_diag_len = 1.0e-10*US%m_to_Z ! A minimal extinction depth to constrain the range of opacity [Z ~> m]
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js,je ; do i=is,ie
         ! Remap opacity (op) to 1/L * tanh(op * L) where L is one Angstrom.
@@ -375,18 +375,18 @@ subroutine opacity_from_chl(optics, sw_total, sw_vis_dir, sw_vis_dif, sw_nir_dir
             enddo
           else
             ! Band 1 is Manizza blue.
-            optics%opacity_band(1,i,j,k) = 0.0232 + 0.074*chl_data(i,j)**0.674
+            optics%opacity_band(1,i,j,k) = (0.0232 + 0.074*chl_data(i,j)**0.674) * US%Z_to_m
             if (nbands >= 2) &  !  Band 2 is Manizza red.
-              optics%opacity_band(2,i,j,k) = 0.225 + 0.037*chl_data(i,j)**0.629
+              optics%opacity_band(2,i,j,k) = (0.225 + 0.037*chl_data(i,j)**0.629) * US%Z_to_m
             ! All remaining bands are NIR, for lack of something better to do.
-            do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86 ; enddo
+            do n=3,nbands ; optics%opacity_band(n,i,j,k) = 2.86*US%Z_to_m ; enddo
           endif
         enddo ; enddo
       case (MOREL_88)
         do j=js,je ; do i=is,ie
           optics%opacity_band(1,i,j,k) = CS%opacity_land_value
           if (G%mask2dT(i,j) > 0.5) &
-            optics%opacity_band(1,i,j,k) = opacity_morel(chl_data(i,j))
+            optics%opacity_band(1,i,j,k) = US%Z_to_m * opacity_morel(chl_data(i,j))
 
           do n=2,optics%nbands
             optics%opacity_band(n,i,j,k) = optics%opacity_band(1,i,j,k)
@@ -460,7 +460,8 @@ subroutine extract_optics_slice(optics, j, G, GV, opacity, opacity_scale, penSW_
   type(ocean_grid_type),   intent(in)  :: G      !< The ocean's grid structure.
   type(verticalGrid_type), intent(in)  :: GV     !< The ocean's vertical grid structure.
   real, dimension(max(optics%nbands,1),SZI_(G),SZK_(GV)), &
-                 optional, intent(out) :: opacity   !< The opacity in each band, i-point, and layer
+                 optional, intent(out) :: opacity   !< The opacity in each band, i-point, and layer [Z-1 ~> m-1],
+                                                    !! but with units that can be altered by opacity_scale.
   real,          optional, intent(in)  :: opacity_scale !< A factor by which to rescale the opacity.
   real, dimension(max(optics%nbands,1),SZI_(G)), &
                  optional, intent(out) :: penSW_top !< The shortwave radiation [Q R Z T-1 ~> W m-2]
@@ -866,7 +867,7 @@ subroutine sumSWoverBands(G, GV, US, h, nsw, optics, j, dt, &
       if (h(i,k) > 0.0) then
         do n=1,nsw ; if (Pen_SW_bnd(n,i) > 0.0) then
           ! SW_trans is the SW that is transmitted THROUGH the layer
-          opt_depth = h(i,k)*GV%H_to_m * optics%opacity_band(n,i,j,k)
+          opt_depth = h(i,k)*GV%H_to_Z * optics%opacity_band(n,i,j,k)
           exp_OD = exp(-opt_depth)
           SW_trans = exp_OD
 
@@ -1015,19 +1016,18 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
       call MOM_mesg('opacity_init: opacity scheme set to "'//trim(tmpstr)//'".', 5)
     endif
     call get_param(param_file, mdl, "PEN_SW_SCALE", CS%pen_sw_scale, &
-                 "The vertical absorption e-folding depth of the "//&
-                 "penetrating shortwave radiation.", units="m", default=0.0)
+                 "The vertical absorption e-folding depth of the penetrating shortwave radiation.", &
+                 units="m", default=0.0, scale=US%m_to_Z)
     !BGR/ Added for opacity_scheme==double_exp read in 2nd exp-decay and fraction
     if (CS%Opacity_scheme == DOUBLE_EXP ) then
       call get_param(param_file, mdl, "PEN_SW_SCALE_2ND", CS%pen_sw_scale_2nd, &
                  "The (2nd) vertical absorption e-folding depth of the "//&
-                 "penetrating shortwave radiation "//&
-                 "(use if SW_EXP_MODE==double.)",&
-                 units="m", default=0.0)
+                 "penetrating shortwave radiation (use if SW_EXP_MODE==double.)", &
+                 units="m", default=0.0, scale=US%m_to_Z)
       call get_param(param_file, mdl, "SW_1ST_EXP_RATIO", CS%sw_1st_exp_ratio, &
                  "The fraction of 1st vertical absorption e-folding depth "//&
                  "penetrating shortwave radiation if SW_EXP_MODE==double.",&
-                  units="m", default=0.0)
+                  units="nondim", default=0.0)
     elseif (CS%OPACITY_SCHEME == Single_Exp) then
       !/Else disable 2nd_exp scheme
       CS%pen_sw_scale_2nd = 0.0
@@ -1094,12 +1094,12 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
 
   call get_param(param_file, mdl, "OPACITY_LAND_VALUE", CS%opacity_land_value, &
                  "The value to use for opacity over land. The default is "//&
-                 "10 m-1 - a value for muddy water.", units="m-1", default=10.0)
+                 "10 m-1 - a value for muddy water.", units="m-1", default=10.0, scale=US%Z_to_m)
 
   CS%warning_issued = .false.
 
   if (.not.allocated(optics%opacity_band)) &
-    allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz))
+    allocate(optics%opacity_band(optics%nbands,isd:ied,jsd:jed,nz), source=0.0)
   if (.not.allocated(optics%sw_pen_band)) &
     allocate(optics%sw_pen_band(optics%nbands,isd:ied,jsd:jed))
   allocate(CS%id_opacity(optics%nbands), source=-1)
@@ -1114,7 +1114,7 @@ subroutine opacity_init(Time, G, GV, US, param_file, diag, CS, optics)
     longname = 'Opacity for shortwave radiation in band '//trim(adjustl(bandnum)) &
       // ', saved as L^-1 tanh(Opacity * L) for L = 10^-10 m'
     CS%id_opacity(n) = register_diag_field('ocean_model', shortname, diag%axesTL, Time, &
-      longname, 'm-1')
+      longname, 'm-1', conversion=US%m_to_Z)
   enddo
 
 end subroutine opacity_init

src/tracer/MOM_generic_tracer.F90

@@ -54,7 +54,7 @@ module MOM_generic_tracer
 
   implicit none ; private
 
-  !> An state hidden in module data that is very much not allowed in MOM6
+  !> A state hidden in module data that is very much not allowed in MOM6
   ! ### This needs to be fixed
   logical :: g_registered = .false.
 
@@ -83,13 +83,8 @@ module MOM_generic_tracer
     !> Pointer to the first element of the linked list of generic tracers.
     type(g_tracer_type), pointer :: g_tracer_list => NULL()
 
-    integer :: H_to_m !< Auxiliary to access GV%H_to_m in routines that do not have access to GV
-
   end type MOM_generic_tracer_CS
 
-! This include declares and sets the variable "version".
-#include "version_variable.h"
-
 contains
 
   !> Initializes the generic tracer packages and adds their tracers to the list
@@ -104,9 +99,12 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
                                                            !! advection and diffusion module.
     type(MOM_restart_CS), target, intent(inout)  :: restart_CS !< MOM restart control struct
 
-! Local variables
+    ! Local variables
     logical :: register_MOM_generic_tracer
 
+    ! This include declares and sets the variable "version".
+#   include "version_variable.h"
+
     character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer'
     character(len=200) :: inputdir ! The directory where NetCDF input files are.
     ! These can be overridden later in via the field manager?
@@ -381,8 +379,6 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file,
     call g_tracer_set_csdiag(CS%diag)
 #endif
 
-    CS%H_to_m = GV%H_to_m
-
   end subroutine initialize_MOM_generic_tracer
 
   !>  Column physics for generic tracers.
@@ -503,7 +499,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     !
     !Calculate tendencies (i.e., field changes at dt) from the sources / sinks
     !
-    if ((G%US%L_to_m == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. (G%US%s_to_T == 1.0)) then
+    if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. &
+        (G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0)) then
       ! Avoid unnecessary copies when no unit conversion is needed.
       call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
                G%areaT, get_diag_time_end(CS%diag), &
@@ -512,7 +509,9 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml,
     else
       call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
                G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
-               optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, &
+               optics%nbands, optics%max_wavelength_band, &
+               sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), &
+               opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), &
                internal_heat=G%US%RZ_to_kg_m2*tv%internal_heat(:,:), &
                frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
     endif
@@ -864,7 +863,7 @@ subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS)
     !nnz: fake rho0
     rho0=1.0
 
-    dzt(:,:,:) = CS%H_to_m * h(:,:,:)
+    dzt(:,:,:) = GV%H_to_m * h(:,:,:)
 
     sosga = global_area_mean(sfc_state%SSS, G)