+Non-Boussinesq thickness diagnostics in kg m-2

  Modified the units for 10 thickness or thickness tendency diagnostics to write
them in their native units (e.g., kg m-2) when run in non-Boussinesq mode. These
changes were proposed in mom-ocean#1565 and
discussed and agreed to at the MOM6 dev call on 4/25/22.  Previously these
diagnostics had been converted to approximate thicknesses in m using a constant
reference density, but this was only accurate to about 0.1% and could be
misinterpreted.  A number of other (more commonly used) thickness-related
variables were already being written in their native units, and this makes the
model's output more self consistent.  All solutions are bitwise identical, but
some diagnostics will change when run in non-Boussinesq mode.

src/ALE/MOM_ALE.F90

@@ -306,7 +306,11 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   type(diag_ctrl), target,    intent(in)  :: diag  !< Diagnostics control structure
   type(ALE_CS), pointer                   :: CS    !< Module control structure
 
+  ! Local variables
+  character(len=48)  :: thickness_units
+
   CS%diag => diag
+  thickness_units = get_thickness_units(GV)
 
   ! These diagnostics of the state variables before ALE are useful for
   ! debugging the ALE code.
@@ -315,7 +319,7 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, &
       'Meridional velocity before remapping', 'm s-1', conversion=US%L_T_to_m_s)
   CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, &
-      'Layer Thickness before remapping', get_thickness_units(GV), conversion=GV%H_to_MKS, &
+      'Layer Thickness before remapping', thickness_units, conversion=GV%H_to_MKS, &
       v_extensive=.true.)
   CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, &
       'Temperature before remapping', 'degC', conversion=US%C_to_degC)
@@ -324,14 +328,15 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   CS%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, Time, &
       'Interface Heights before remapping', 'm', conversion=US%Z_to_m)
 
-  CS%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,Time, &
+  CS%id_dzRegrid = register_diag_field('ocean_model', 'dzRegrid', diag%axesTi, Time, &
       'Change in interface height due to ALE regridding', 'm', conversion=GV%H_to_m)
-  cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', &
-      diag%axestl, time, 'layer thicknesses after ALE regridding and remapping', &
-      'm', conversion=GV%H_to_m, v_extensive=.true.)
-  cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, &
+  cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, &
+      'layer thicknesses after ALE regridding and remapping', &
+      thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
+  cs%id_vert_remap_h_tendency = register_diag_field('ocean_model', &
+      'vert_remap_h_tendency', diag%axestl, Time, &
       'Layer thicknesses tendency due to ALE regridding and remapping', &
-      'm s-1', conversion=GV%H_to_m*US%s_to_T, v_extensive=.true.)
+      trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
 
 end subroutine ALE_register_diags
 

src/diagnostics/MOM_diagnostics.F90

@@ -31,7 +31,7 @@ module MOM_diagnostics
 use MOM_unit_scaling,      only : unit_scale_type
 use MOM_variables,         only : thermo_var_ptrs, ocean_internal_state, p3d
 use MOM_variables,         only : accel_diag_ptrs, cont_diag_ptrs, surface
-use MOM_verticalGrid,      only : verticalGrid_type, get_thickness_units
+use MOM_verticalGrid,      only : verticalGrid_type, get_thickness_units, get_flux_units
 use MOM_wave_speed,        only : wave_speed, wave_speed_CS, wave_speed_init
 
 implicit none ; private
@@ -1593,11 +1593,9 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
                  "forms of the same expressions.", default=default_2018_answers)
   call get_param(param_file, mdl, "SPLIT", split, default=.true., do_not_log=.true.)
 
-  if (GV%Boussinesq) then
-    thickness_units = "m" ; flux_units = "m3 s-1" ; convert_H = GV%H_to_m
-  else
-    thickness_units = "kg m-2" ; flux_units = "kg s-1" ; convert_H = GV%H_to_kg_m2
-  endif
+  thickness_units = get_thickness_units(GV)
+  flux_units = get_flux_units(GV)
+  convert_H = GV%H_to_MKS
 
   CS%id_masscello = register_diag_field('ocean_model', 'masscello', diag%axesTL,&
       Time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2', & !### , conversion=GV%H_to_kg_m2, &
@@ -1607,11 +1605,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
       diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass')
 
   CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, &
-      long_name = 'Cell Thickness', standard_name='cell_thickness', &
+      long_name='Cell Thickness', standard_name='cell_thickness', &
       units='m', conversion=US%Z_to_m, v_extensive=.true.)
   CS%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, Time, &
-      long_name = 'Cell thickness from the previous timestep', &
-      units='m', conversion=GV%H_to_m, v_extensive=.true.)
+      long_name='Cell thickness from the previous timestep', &
+      units=thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
 
   ! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and
   ! must be registered first. We earlier stored the handle of volcello but need it here for posting
@@ -1948,10 +1946,11 @@ subroutine register_transport_diags(Time, G, GV, US, IDs, diag)
   character(len=48) :: thickness_units, accum_flux_units
 
   thickness_units = get_thickness_units(GV)
+  H_convert = GV%H_to_MKS
   if (GV%Boussinesq) then
-    H_convert = GV%H_to_m ; accum_flux_units = "m3"
+    accum_flux_units = "m3"
   else
-    H_convert = GV%H_to_kg_m2 ; accum_flux_units = "kg"
+    accum_flux_units = "kg"
   endif
 
   ! Diagnostics related to tracer and mass transport
@@ -1979,10 +1978,10 @@ subroutine register_transport_diags(Time, G, GV, US, IDs, diag)
       standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')
   IDs%id_dynamics_h = register_diag_field('ocean_model','dynamics_h',  &
       diag%axesTl, Time, 'Layer thicknesses prior to horizontal dynamics', &
-      'm', v_extensive=.true., conversion=GV%H_to_m)
+      thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
   IDs%id_dynamics_h_tendency = register_diag_field('ocean_model','dynamics_h_tendency',  &
       diag%axesTl, Time, 'Change in layer thicknesses due to horizontal dynamics', &
-      'm s-1', v_extensive=.true., conversion=GV%H_to_m*US%s_to_T)
+      trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
 
 end subroutine register_transport_diags
 

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -3253,8 +3253,9 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
   ! available only for ALE algorithm.
   ! diagnostics for tendencies of temp and heat due to frazil
   CS%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, Time, &
-      long_name='Cell thickness used during diabatic diffusion', &
-      units='m', conversion=GV%H_to_m, v_extensive=.true.)
+      'Cell thickness used during diabatic diffusion', &
+      thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
+
   if (CS%useALEalgorithm) then
     CS%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &
         'diabatic_diff_temp_tendency', diag%axesTL, Time,              &
@@ -3326,12 +3327,12 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
     ! available only for ALE algorithm.
   ! diagnostics for tendencies of temp and heat due to frazil
     CS%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, Time, &
-        long_name='Cell thickness after applying boundary forcing', &
-        units='m', conversion=GV%H_to_m, v_extensive=.true.)
+        'Cell thickness after applying boundary forcing', &
+        thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
     CS%id_boundary_forcing_h_tendency = register_diag_field('ocean_model',   &
         'boundary_forcing_h_tendency', diag%axesTL, Time,                &
         'Cell thickness tendency due to boundary forcing', &
-        'm s-1', conversion=GV%H_to_m*US%s_to_T, v_extensive=.true.)
+        trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
     if (CS%id_boundary_forcing_h_tendency > 0) then
       CS%boundary_forcing_tendency_diag = .true.
     endif
@@ -3388,7 +3389,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
   ! diagnostics for tendencies of temp and heat due to frazil
   CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, &
       long_name='Cell Thickness', standard_name='cell_thickness', &
-      units='m', conversion=GV%H_to_m, v_extensive=.true.)
+      units=thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
 
   ! diagnostic for tendency of temp due to frazil
   CS%id_frazil_temp_tend = register_diag_field('ocean_model',&

src/tracer/MOM_offline_main.F90

@@ -35,7 +35,7 @@ module MOM_offline_main
 use MOM_tracer_registry,      only : tracer_registry_type, MOM_tracer_chksum, MOM_tracer_chkinv
 use MOM_unit_scaling,         only : unit_scale_type
 use MOM_variables,            only : thermo_var_ptrs
-use MOM_verticalGrid,         only : verticalGrid_type
+use MOM_verticalGrid,         only : verticalGrid_type, get_thickness_units
 
 implicit none ; private
 
@@ -1160,7 +1160,7 @@ subroutine register_diags_offline_transport(Time, diag, CS, GV, US)
     'at the end of the offline timestep', 'm', conversion=GV%H_to_m)
   CS%id_h_redist = register_diag_field('ocean_model','h_redist', diag%axesTL, Time, &
     'Layer thicknesses before redistribution of mass fluxes', &
-    'm', conversion=GV%H_to_m)
+    get_thickness_units(GV), conversion=GV%H_to_MKS)
 
   ! Regridded/remapped input fields
   CS%id_uhtr_regrid = register_diag_field('ocean_model', 'uhtr_regrid', diag%axesCuL, Time, &