+Standardize diagnostic rescaling and get_param calls

src/diagnostics/MOM_diagnostics.F90

@@ -1293,32 +1293,27 @@ subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv
 
   real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
   real, dimension(SZI_(G),SZJ_(G)) :: &
-    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m]
+    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [Z ~> m]
   real :: I_time_int    ! The inverse of the time interval [T-1 ~> s-1].
-  real :: zos_area_mean ! Global area mean sea surface height [m]
+  real :: zos_area_mean ! Global area mean sea surface height [Z ~> m]
   real :: volo          ! Total volume of the ocean [m3]
-  real :: ssh_ga        ! Global ocean area weighted mean sea seaface height [m]
+  real :: ssh_ga        ! Global ocean area weighted mean sea seaface height [Z ~> m]
   integer :: i, j, is, ie, js, je
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
   ! area mean SSH
   if (IDs%id_ssh_ga > 0) then
-    ssh_ga = global_area_mean(ssh, G, scale=US%Z_to_m)
+    ssh_ga = global_area_mean(ssh, G, tmp_scale=US%Z_to_m)
     call post_data(IDs%id_ssh_ga, ssh_ga, diag)
   endif
 
   ! post the dynamic sea level, zos, and zossq.
-  ! zos is ave_ssh with sea ice inverse barometer removed,
-  ! and with zero global area mean.
+  ! zos is ave_ssh with sea ice inverse barometer removed, and with zero global area mean.
   if (IDs%id_zos > 0 .or. IDs%id_zossq > 0) then
-    zos(:,:) = 0.0
+    zos_area_mean = global_area_mean(ssh_ibc, G, tmp_scale=US%Z_to_m)
     do j=js,je ; do i=is,ie
-      zos(i,j) = US%Z_to_m*ssh_ibc(i,j)
-    enddo ; enddo
-    zos_area_mean = global_area_mean(zos, G)
-    do j=js,je ; do i=is,ie
-      zos(i,j) = zos(i,j) - G%mask2dT(i,j)*zos_area_mean
+      zos(i,j) = ssh_ibc(i,j) - G%mask2dT(i,j)*zos_area_mean
     enddo ; enddo
     if (IDs%id_zos > 0) call post_data(IDs%id_zos, zos, diag, mask=G%mask2dT)
     if (IDs%id_zossq > 0) then
@@ -1844,14 +1839,14 @@ subroutine register_surface_diags(Time, G, US, IDs, diag, tv)
       standard_name='sea_water_volume')
   IDs%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, Time,&
       standard_name = 'sea_surface_height_above_geoid',                   &
-      long_name= 'Sea surface height above geoid', units='m')
+      long_name= 'Sea surface height above geoid', units='m', conversion=US%Z_to_m)
   IDs%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, Time,&
       standard_name='square_of_sea_surface_height_above_geoid',             &
-      long_name='Square of sea surface height above geoid', units='m2')
+      long_name='Square of sea surface height above geoid', units='m2', conversion=US%Z_to_m**2)
   IDs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, Time, &
       'Sea Surface Height', 'm', conversion=US%Z_to_m)
   IDs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', Time, diag,&
-      long_name='Area averaged sea surface height', units='m',            &
+      long_name='Area averaged sea surface height', units='m', conversion=US%Z_to_m, &
       standard_name='area_averaged_sea_surface_height')
   IDs%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, Time, &
       'Sea Surface Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)

src/framework/MOM_diag_mediator.F90

@@ -2775,7 +2775,7 @@ end subroutine attach_cell_methods
 function register_scalar_field(module_name, field_name, init_time, diag_cs, &
             long_name, units, missing_value, range, standard_name, &
             do_not_log, err_msg, interp_method, cmor_field_name, &
-            cmor_long_name, cmor_units, cmor_standard_name)
+            cmor_long_name, cmor_units, cmor_standard_name, conversion)
   integer :: register_scalar_field !< An integer handle for a diagnostic array.
   character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model"
                                               !! or "ice_shelf_model"
@@ -2796,6 +2796,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, &
   character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field
   character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field
   character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field
+  real,             optional, intent(in) :: conversion !< A value to multiply data by before writing to file
 
   ! Local variables
   real :: MOM_missing_value
@@ -2826,6 +2827,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, &
     call assert(associated(diag), 'register_scalar_field: diag allocation failed')
     diag%fms_diag_id = fms_id
     diag%debug_str = trim(module_name)//"-"//trim(field_name)
+    if (present(conversion)) diag%conversion_factor = conversion
   endif
 
   if (present(cmor_field_name)) then
@@ -2856,6 +2858,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, &
       call alloc_diag_with_id(dm_id, diag_cs, cmor_diag)
       cmor_diag%fms_diag_id = fms_id
       cmor_diag%debug_str = trim(module_name)//"-"//trim(cmor_field_name)
+      if (present(conversion)) cmor_diag%conversion_factor = conversion
     endif
   endif
 

src/ice_shelf/MOM_ice_shelf.F90

@@ -136,7 +136,7 @@ module MOM_ice_shelf
 
   !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!!
 
-  real :: time_step    !< this is the shortest timestep that the ice shelf sees, and
+  real :: time_step    !< this is the shortest timestep that the ice shelf sees [T ~> s], and
                        !! is equal to the forcing timestep (it is passed in when the shelf
                        !! is initialized - so need to reorganize MOM driver.
                        !! it will be the prognistic timestep ... maybe.
@@ -153,8 +153,9 @@ module MOM_ice_shelf
   real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
   real :: T0                !< temperature at ocean surface in the restoring region [degC]
   real :: S0                !< Salinity at ocean surface in the restoring region [ppt].
-  real :: input_flux        !< Ice volume flux at an upstream open boundary [m3 s-1].
-  real :: input_thickness   !< Ice thickness at an upstream open boundary [m].
+  real :: input_flux        !< The vertically integrated inward ice thickness flux per
+                            !! unit face length at an upstream boundary [Z L T-1 ~> m2 s-1]
+  real :: input_thickness   !< Ice thickness at an upstream open boundary [Z ~> m].
 
   type(time_type) :: Time                !< The component's time.
   type(EOS_type) :: eqn_of_state         !< Type that indicates the
@@ -368,7 +369,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
   CS%Time = Time
 
   if (CS%override_shelf_movement) then
-    CS%time_step = time_step
+    CS%time_step = US%s_to_T*time_step
     ! update shelf mass
     if (CS%mass_from_file) call update_shelf_mass(G, US, CS, ISS, Time)
   endif
@@ -1109,7 +1110,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
 
     ! take into account changes in mass (or thickness) when imposing ice shelf mass
     if (CS%override_shelf_movement .and. CS%mass_from_file) then
-      dTime = real_to_time(CS%time_step)
+      dTime = real_to_time(US%T_to_s*CS%time_step)
 
       ! Compute changes in mass after at least one full time step
       if (CS%Time > dTime) then
@@ -1144,8 +1145,8 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
             delta_float_mass(i,j) = 0.0
           endif
         enddo ; enddo
-        delta_mass_shelf = US%kg_m2s_to_RZ_T*(global_area_integral(delta_float_mass, G, scale=US%RZ_to_kg_m2, &
-                                                                   area=ISS%area_shelf_h) / CS%time_step)
+        delta_mass_shelf = global_area_integral(delta_float_mass, G, tmp_scale=US%RZ_to_kg_m2, &
+                                                area=ISS%area_shelf_h) / CS%time_step
       else! first time step
         delta_mass_shelf = 0.0
       endif
@@ -1166,8 +1167,8 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
 
     balancing_area = global_area_integral(bal_frac, G)
     if (balancing_area > 0.0) then
-      balancing_flux = ( US%kg_m2s_to_RZ_T*global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, &
-                                                                area=ISS%area_shelf_h) + &
+      balancing_flux = ( global_area_integral(ISS%water_flux, G, tmp_scale=US%RZ_T_to_kg_m2s, &
+                                              area=ISS%area_shelf_h) + &
                          delta_mass_shelf ) / balancing_area
     else
       balancing_flux = 0.0
@@ -1184,7 +1185,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes)
     enddo ; enddo
 
     if (CS%debug) then
-      write(mesg,*) 'Balancing flux (kg/(m^2 s)), dt = ', balancing_flux*US%RZ_T_to_kg_m2s, CS%time_step
+      write(mesg,*) 'Balancing flux (kg/(m^2 s)), dt = ', balancing_flux*US%RZ_T_to_kg_m2s, US%T_to_s*CS%time_step
       call MOM_mesg(mesg)
       call MOM_forcing_chksum("After constant sea level", fluxes, G, CS%US, haloshift=0)
     endif
@@ -1487,15 +1488,15 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
   if (CS%constant_sea_level) CS%min_ocean_mass_float = dz_ocean_min_float*CS%Rho_ocn
 
   call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", CS%flux_factor, &
-                 "Non-dimensional factor applied to shelf thermodynamic "//&
-                 "fluxes.", units="none", default=1.0)
+                 "Non-dimensional factor applied to shelf thermodynamic fluxes.", &
+                 units="none", default=1.0)
 
   call get_param(param_file, mdl, "KV_ICE", CS%kv_ice, &
                  "The viscosity of the ice.", &
                  units="m2 s-1", default=1.0e10, scale=US%Z_to_L**2*US%m_to_L**2*US%T_to_s)
   call get_param(param_file, mdl, "KV_MOLECULAR", CS%kv_molec, &
-                 "The molecular kinimatic viscosity of sea water at the "//&
-                 "freezing temperature.", units="m2 s-1", default=1.95e-6, scale=US%m2_s_to_Z2_T)
+                 "The molecular kinimatic viscosity of sea water at the freezing temperature.", &
+                 units="m2 s-1", default=1.95e-6, scale=US%m2_s_to_Z2_T)
   call get_param(param_file, mdl, "ICE_SHELF_SALINITY", CS%Salin_ice, &
                  "The salinity of the ice inside the ice shelf.", units="psu", &
                  default=0.0)
@@ -1511,7 +1512,7 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
   call get_param(param_file, mdl, "DT_FORCING", CS%time_step, &
                  "The time step for changing forcing, coupling with other "//&
                  "components, or potentially writing certain diagnostics. "//&
-                 "The default value is given by DT.", units="s", default=0.0)
+                 "The default value is given by DT.", units="s", default=0.0, scale=US%s_to_T)
 
   call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", col_thick_melt_thresh, &
                  "The minimum ocean column thickness where melting is allowed.", &
@@ -1571,9 +1572,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces_in,
                  "A typical density of ice.", units="kg m-3", default=917.0, scale=US%kg_m3_to_R)
 
     call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", CS%input_flux, &
-                 "volume flux at upstream boundary", units="m2 s-1", default=0.)
+                 "volume flux at upstream boundary", units="m2 s-1", default=0., scale=US%m_to_Z*US%m_s_to_L_T)
     call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", CS%input_thickness, &
-                 "flux thickness at upstream boundary", units="m", default=1000.)
+                 "flux thickness at upstream boundary", units="m", default=1000., scale=US%m_to_Z)
   else
     ! This is here because of inconsistent defaults.  I don't know why.  RWH
     call get_param(param_file, mdl, "DENSITY_ICE", CS%density_ice, &
@@ -2005,7 +2006,7 @@ end subroutine initialize_shelf_mass
 !> This subroutine applies net accumulation/ablation at the top surface to the dynamic ice shelf.
 !>>acc_rate[m-s]=surf_mass_flux/density_ice is ablation/accumulation rate
 !>>positive for accumulation negative for ablation
-subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes,time_step, Time)
+subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time)
   type(ice_shelf_CS),    intent(in)    :: CS  !< A pointer to the ice shelf control structure
   type(ocean_grid_type), intent(inout) :: G  !< The ocean's grid structure.
   type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
@@ -2113,8 +2114,8 @@ end subroutine update_shelf_mass
 subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes)
   type(ice_shelf_CS),         pointer    :: CS !< ice shelf control structure
   type(ocean_grid_type), intent(in)      :: G  !< A pointer to an ocean grid control structure.
-  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: frac_shelf_h !< Ice shelf area fraction [nodim].
-  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: mass_shelf !<Ice shelf mass [R Z -> kg m-2]
+  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: frac_shelf_h !< Ice shelf area fraction [nondim].
+  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: mass_shelf !<Ice shelf mass [R Z ~> kg m-2]
   logical, optional                      :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using
                                                !! the data_override capability (only for MOSAIC grids)
 

src/ice_shelf/MOM_marine_ice.F90

@@ -176,8 +176,9 @@ subroutine marine_ice_init(Time, G, param_file, diag, CS)
   type(param_file_type),   intent(in)    :: param_file !< Runtime parameter handles
   type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
   type(marine_ice_CS),     pointer       :: CS   !< Pointer to the control structure for MOM_marine_ice
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   character(len=40)  :: mdl = "MOM_marine_ice"  ! This module's name.
 
   if (associated(CS)) then
@@ -197,8 +198,8 @@ subroutine marine_ice_init(Time, G, param_file, diag, CS)
                  "The latent heat of fusion.", units="J/kg", default=hlf, scale=G%US%J_kg_to_Q)
   call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", CS%berg_area_threshold, &
                  "Fraction of grid cell which iceberg must occupy, so that fluxes "//&
-                 "below berg are set to zero. Not applied for negative "//&
-                 "values.", units="non-dim", default=-1.0)
+                 "below berg are set to zero. Not applied for negative values.", &
+                 units="nondim", default=-1.0)
 
 end subroutine marine_ice_init
 

src/initialization/MOM_shared_initialization.F90

@@ -184,7 +184,7 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US)
   type(unit_scale_type),            intent(in)    :: US !< A dimensional unit scaling type
 
   ! Local variables
-  real, dimension(:), allocatable :: new_depth ! The new values of the depths [m]
+  real, dimension(:), allocatable :: new_depth ! The new values of the depths [Z ~> m]
   integer, dimension(:), allocatable :: ig, jg ! The global indicies of the points to modify
   character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
   character(len=40)  :: mdl = "apply_topography_edits_from_file" ! This subroutine's name.
@@ -247,22 +247,22 @@ subroutine apply_topography_edits_from_file(D, G, param_file, US)
   ! Read iEdit, jEdit and zEdit
   call read_variable(topo_edits_file, 'iEdit', ig, ncid_in=ncid)
   call read_variable(topo_edits_file, 'jEdit', jg, ncid_in=ncid)
-  call read_variable(topo_edits_file, 'zEdit', new_depth, ncid_in=ncid)
+  call read_variable(topo_edits_file, 'zEdit', new_depth, ncid_in=ncid, scale=US%m_to_Z)
   call close_file_to_read(ncid, topo_edits_file)
 
   do n = 1, n_edits
     i = ig(n) - G%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1
     j = jg(n) - G%jsd_global + 2
     if (i>=G%isc .and. i<=G%iec .and. j>=G%jsc .and. j<=G%jec) then
-      if (new_depth(n)*US%m_to_Z /= mask_depth) then
+      if (new_depth(n) /= mask_depth) then
         write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
-          'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)*US%Z_to_m, '->', abs(new_depth(n)), i, j
-        D(i,j) = abs(US%m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
+          'Ocean topography edit: ', n, ig(n), jg(n), D(i,j)*US%Z_to_m, '->', abs(US%Z_to_m*new_depth(n)), i, j
+        D(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
       else
         if (topo_edits_change_mask) then
           write(stdout,'(a,3i5,f8.2,a,f8.2,2i4)') &
-            'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)*US%Z_to_m,'->',abs(new_depth(n)),i,j
-            D(i,j) = abs(US%m_to_Z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
+            'Ocean topography edit: ',n,ig(n),jg(n),D(i,j)*US%Z_to_m,'->',abs(US%Z_to_m*new_depth(n)),i,j
+            D(i,j) = abs(new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
         else
           call MOM_error(FATAL, ' apply_topography_edits_from_file: '//&
             "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
@@ -454,8 +454,8 @@ subroutine set_rotation_planetary(f, G, param_file, US)
   call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
 
   call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, &
-                 "The rotation rate of the earth.", units="s-1", &
-                 default=7.2921e-5, scale=US%T_to_s)
+                 "The rotation rate of the earth.", &
+                 units="s-1", default=7.2921e-5, scale=US%T_to_s)
   PI = 4.0*atan(1.0)
 
   do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB
@@ -893,13 +893,13 @@ subroutine reset_face_lengths_list(G, param_file, US)
     allocate(v_line_used(num_lines), source=0)
     allocate(v_line_no(num_lines), source=0)
 
-    allocate(Dmin_u(num_lines))    ; Dmin_u(:) = 0.0
-    allocate(Dmax_u(num_lines))    ; Dmax_u(:) = 0.0
-    allocate(Davg_u(num_lines))    ; Davg_u(:) = 0.0
+    allocate(Dmin_u(num_lines), source=0.0)
+    allocate(Dmax_u(num_lines), source=0.0)
+    allocate(Davg_u(num_lines), source=0.0)
 
-    allocate(Dmin_v(num_lines))    ; Dmin_v(:) = 0.0
-    allocate(Dmax_v(num_lines))    ; Dmax_v(:) = 0.0
-    allocate(Davg_v(num_lines))    ; Davg_v(:) = 0.0
+    allocate(Dmin_v(num_lines), source=0.0)
+    allocate(Dmax_v(num_lines), source=0.0)
+    allocate(Davg_v(num_lines), source=0.0)
 
     ! Actually read the lines.
     if (is_root_pe()) then

src/initialization/MOM_state_initialization.F90

@@ -1179,7 +1179,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
                  fail_if_missing=.not.just_read, do_not_log=just_read)
   call get_param(PF, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
                  "The initial condition variable for the surface pressure exerted by ice.", &
-                 units="Pa", default="", do_not_log=just_read)
+                 default="", do_not_log=just_read)
   call get_param(PF, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
   filename = trim(slasher(inputdir))//trim(p_surf_file)
   if (.not.just_read) call log_param(PF,  mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)