+Remove rescaling factors from restart files

  Remove the code to account for unit rescaling within the restart files.  This
rescaling within the restart files has not been used in the code since March,
2022, and the model will work with older restart files provided that they did
not use dimensional rescaling, and even if they did they can be converted not to
use rescaling with a short run with the older code that created them.  Also
removed the publicly visible routines fix_restart_scaling and eliminated the
m_to_H_restart element of the verticalGrid_type; in any cases of non-standard
code using this element, it should be replaced with 1.0.  The various
US%..._restart elements and fix_restart_unit_scaling are being retained for now
because they are still being used in the SIS2 code.

  These changes significantly simplify the code, and they lead to a handful
of constants that are always 1 not being included in the MOM6 restart files.
All answers are bitwise identical, but a publicly visible interface has been
eliminated, as has been an element (GV%m_to_H_restart) of a transparent type.

src/core/MOM.F90

@@ -135,14 +135,12 @@ module MOM
 use MOM_tracer_flow_control,   only : tracer_flow_control_init, call_tracer_surface_state
 use MOM_tracer_flow_control,   only : tracer_flow_control_end, call_tracer_register_obc_segments
 use MOM_transcribe_grid,       only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid
-use MOM_unit_scaling,          only : unit_scale_type, unit_scaling_init
-use MOM_unit_scaling,          only : unit_scaling_end, fix_restart_unit_scaling
+use MOM_unit_scaling,          only : unit_scale_type, unit_scaling_init, unit_scaling_end
 use MOM_variables,             only : surface, allocate_surface_state, deallocate_surface_state
 use MOM_variables,             only : thermo_var_ptrs, vertvisc_type, porous_barrier_type
 use MOM_variables,             only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state
 use MOM_variables,             only : rotate_surface_state
 use MOM_verticalGrid,          only : verticalGrid_type, verticalGridInit, verticalGridEnd
-use MOM_verticalGrid,          only : fix_restart_scaling
 use MOM_verticalGrid,          only : get_thickness_units, get_flux_units, get_tr_flux_units
 use MOM_wave_interface,        only : wave_parameters_CS, waves_end, waves_register_restarts
 use MOM_wave_interface,        only : Update_Stokes_Drift
@@ -1978,7 +1976,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   real :: conv2watt            ! A conversion factor from temperature fluxes to heat
                                ! fluxes [J m-2 H-1 C-1 ~> J m-3 degC-1 or J kg-1 degC-1]
   real :: conv2salt            ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
-  real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors
   character(len=48) :: S_flux_units
 
   type(vardesc) :: vd_T, vd_S  ! Structures describing temperature and salinity variables.
@@ -3117,16 +3114,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   call register_obsolete_diagnostics(param_file, CS%diag)
 
   if (use_frazil) then
-    if (query_initialized(CS%tv%frazil, "frazil", restart_CSp)) then
-      ! Test whether the dimensional rescaling has changed for heat content.
-      if ((US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart /= 0.0) .and. &
-          (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart /= 1.0) ) then
-        QRZ_rescale = 1.0 / (US%J_kg_to_Q_restart*US%kg_m3_to_R_restart*US%m_to_Z_restart)
-        do j=js,je ; do i=is,ie
-          CS%tv%frazil(i,j) = QRZ_rescale * CS%tv%frazil(i,j)
-        enddo ; enddo
-      endif
-    else
+    if (.not.query_initialized(CS%tv%frazil, "frazil", restart_CSp)) then
       CS%tv%frazil(:,:) = 0.0
       call set_initialized(CS%tv%frazil, "frazil", restart_CSp)
     endif
@@ -3136,39 +3124,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     CS%p_surf_prev_set = query_initialized(CS%p_surf_prev, "p_surf_prev", restart_CSp)
 
     if (CS%p_surf_prev_set) then
-      ! Test whether the dimensional rescaling has changed for pressure.
-      if ((US%kg_m3_to_R_restart*US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. &
-          (US%s_to_T_restart**2 /= US%kg_m3_to_R_restart * US%m_to_L_restart**2) ) then
-        RL2_T2_rescale = US%s_to_T_restart**2 / (US%kg_m3_to_R_restart*US%m_to_L_restart**2)
-        do j=js,je ; do i=is,ie
-          CS%p_surf_prev(i,j) = RL2_T2_rescale * CS%p_surf_prev(i,j)
-        enddo ; enddo
-      endif
-
       call pass_var(CS%p_surf_prev, G%domain)
     endif
   endif
 
-  if (use_ice_shelf .and. associated(CS%Hml)) then
-    if (query_initialized(CS%Hml, "hML", restart_CSp)) then
-      ! Test whether the dimensional rescaling has changed for depths.
-      if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0) ) then
-        Z_rescale = 1.0 / US%m_to_Z_restart
-        do j=js,je ; do i=is,ie
-          CS%Hml(i,j) = Z_rescale * CS%Hml(i,j)
-        enddo ; enddo
-      endif
-    endif
-  endif
-
-  if (query_initialized(CS%ave_ssh_ibc, "ave_ssh", restart_CSp)) then
-    if ((US%m_to_Z_restart /= 0.0) .and. (US%m_to_Z_restart /= 1.0) ) then
-      Z_rescale = 1.0 / US%m_to_Z_restart
-      do j=js,je ; do i=is,ie
-        CS%ave_ssh_ibc(i,j) = Z_rescale * CS%ave_ssh_ibc(i,j)
-      enddo ; enddo
-    endif
-  else
+  if (.not.query_initialized(CS%ave_ssh_ibc, "ave_ssh", restart_CSp)) then
     if (CS%split) then
       call find_eta(CS%h, CS%tv, G, GV, US, CS%ave_ssh_ibc, eta, dZref=G%Z_ref)
     else
@@ -3195,10 +3155,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   ! initialize stochastic physics
   call stochastics_init(CS%dt_therm, CS%G, CS%GV, CS%stoch_CS, param_file, diag, Time)
 
-  !### This could perhaps go here instead of in finish_MOM_initialization?
-  ! call fix_restart_scaling(GV)
-  ! call fix_restart_unit_scaling(US)
-
   call callTree_leave("initialize_MOM()")
   call cpu_clock_end(id_clock_init)
 
@@ -3226,11 +3182,6 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
   ! Pointers for convenience
   G => CS%G ; GV => CS%GV ; US => CS%US
 
-  !### Move to initialize_MOM?
-  call fix_restart_scaling(GV, unscaled=.true.)
-  call fix_restart_unit_scaling(US, unscaled=.true.)
-
-
   if (CS%use_particles) then
     call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v)
   endif
@@ -3382,18 +3333,6 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp)
   endif
 
   ! Register scalar unit conversion factors.
-  call register_restart_field(US%m_to_Z_restart, "m_to_Z", .false., restart_CSp, &
-                              "Height unit conversion factor", "Z meter-1")
-  call register_restart_field(GV%m_to_H_restart, "m_to_H", .false., restart_CSp, &
-                              "Thickness unit conversion factor", "H meter-1")
-  call register_restart_field(US%m_to_L_restart, "m_to_L", .false., restart_CSp, &
-                              "Length unit conversion factor", "L meter-1")
-  call register_restart_field(US%s_to_T_restart, "s_to_T", .false., restart_CSp, &
-                              "Time unit conversion factor", "T second-1")
-  call register_restart_field(US%kg_m3_to_R_restart, "kg_m3_to_R", .false., restart_CSp, &
-                              "Density unit conversion factor", "R m3 kg-1")
-  call register_restart_field(US%J_kg_to_Q_restart, "J_kg_to_Q", .false., restart_CSp, &
-                              "Heat content unit conversion factor.", units="Q kg J-1")
   call register_restart_field(CS%first_dir_restart, "First_direction", .false., restart_CSp, &
                               "Indicator of the first direction in split calculations.", "nondim")
 

src/core/MOM_barotropic.F90

@@ -4318,8 +4318,6 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
                                        ! drag piston velocity.
   character(len=80)  :: wave_drag_var  ! The wave drag piston velocity variable
                                        ! name in wave_drag_file.
-  real :: vel_rescale ! A rescaling factor for horizontal velocity from the representation in
-                      ! a restart file to the internal representation in this run.
   real :: mean_SL     ! The mean sea level that is used along with the bathymetry to estimate the
                       ! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m].
   real :: det_de      ! The partial derivative due to self-attraction and loading of the reference
@@ -4788,8 +4786,6 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
   dtbt_tmp = -1.0
   if (query_initialized(CS%dtbt, "DTBT", restart_CS)) then
     dtbt_tmp = CS%dtbt
-    if ((US%s_to_T_restart /= 0.0) .and. (US%s_to_T_restart /= 1.0)) &
-      dtbt_tmp = (1.0 / US%s_to_T_restart) * CS%dtbt
   endif
 
   ! Estimate the maximum stable barotropic time step.
@@ -4948,23 +4944,13 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
     do k=1,nz ; do J=js-1,je ; do i=is,ie
       CS%vbtav(i,J) = CS%vbtav(i,J) + CS%frhatv(i,J,k) * v(i,J,k)
     enddo ; enddo ; enddo
-  elseif ((US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. &
-          (US%s_to_T_restart /= US%m_to_L_restart)) then
-    vel_rescale = US%s_to_T_restart / US%m_to_L_restart
-    do j=js,je ; do I=is-1,ie ; CS%ubtav(I,j) = vel_rescale * CS%ubtav(I,j) ; enddo ; enddo
-    do J=js-1,je ; do i=is,ie ; CS%vbtav(i,J) = vel_rescale * CS%vbtav(i,J) ; enddo ; enddo
   endif
 
   if (CS%gradual_BT_ICs) then
     if (.NOT.query_initialized(CS%ubt_IC,"ubt_IC",restart_CS) .or. &
         .NOT.query_initialized(CS%vbt_IC,"vbt_IC",restart_CS)) then
       do j=js,je ; do I=is-1,ie ; CS%ubt_IC(I,j) = CS%ubtav(I,j) ; enddo ; enddo
       do J=js-1,je ; do i=is,ie ; CS%vbt_IC(i,J) = CS%vbtav(i,J) ; enddo ; enddo
-    elseif ((US%s_to_T_restart*US%m_to_L_restart /= 0.0) .and. &
-            (US%s_to_T_restart /= US%m_to_L_restart)) then
-      vel_rescale = US%s_to_T_restart / US%m_to_L_restart
-      do j=js,je ; do I=is-1,ie ; CS%ubt_IC(I,j) = vel_rescale * CS%ubt_IC(I,j) ; enddo ; enddo
-      do J=js-1,je ; do i=is,ie ; CS%vbt_IC(i,J) = vel_rescale * CS%vbt_IC(i,J) ; enddo ; enddo
     endif
   endif
 !   Calculate other constants which are used for btstep.

src/core/MOM_dynamics_split_RK2.F90

@@ -1246,14 +1246,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=48) :: thickness_units, flux_units, eta_rest_name
-  real :: H_rescale  ! A rescaling factor for thicknesses from the representation in a
-                     ! restart file to the internal representation in this run  [various units ~> 1]
-  real :: vel_rescale ! A rescaling factor for velocities from the representation in a
-                     ! restart file to the internal representation in this run  [various units ~> 1]
-  real :: uH_rescale ! A rescaling factor for thickness transports from the representation in a
-                     ! restart file to the internal representation in this run  [various units ~> 1]
-  real :: accel_rescale ! A rescaling factor for accelerations from the representation in a
-                     ! restart file to the internal representation in this run  [various units ~> 1]
   type(group_pass_type) :: pass_av_h_uvh
   logical :: debug_truncations
   logical :: read_uv, read_h2
@@ -1410,9 +1402,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       CS%eta(i,j) = CS%eta(i,j) + h(i,j,k)
     enddo ; enddo ; enddo
     call set_initialized(CS%eta, trim(eta_rest_name), restart_CS)
-  elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then
-    H_rescale = 1.0 / GV%m_to_H_restart
-    do j=js,je ; do i=is,ie ; CS%eta(i,j) = H_rescale * CS%eta(i,j) ; enddo ; enddo
   endif
   ! Copy eta into an output array.
   do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo
@@ -1427,17 +1416,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
                               OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp)
     call set_initialized(CS%diffu, "diffu", restart_CS)
     call set_initialized(CS%diffv, "diffv", restart_CS)
-  else
-    if ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
-         (US%s_to_T_restart**2 /= US%m_to_L_restart) ) then
-      accel_rescale = US%s_to_T_restart**2 / US%m_to_L_restart
-      do k=1,nz ; do j=js,je ; do I=G%IscB,G%IecB
-        CS%diffu(I,j,k) = accel_rescale * CS%diffu(I,j,k)
-      enddo ; enddo ; enddo
-      do k=1,nz ; do J=G%JscB,G%JecB ; do i=is,ie
-        CS%diffv(i,J,k) = accel_rescale * CS%diffv(i,J,k)
-      enddo ; enddo ; enddo
-    endif
   endif
 
   if (.not. query_initialized(CS%u_av, "u2", restart_CS) .or. &
@@ -1446,11 +1424,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
     do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; CS%v_av(i,J,k) = v(i,J,k) ; enddo ; enddo ; enddo
     call set_initialized(CS%u_av, "u2", restart_CS)
     call set_initialized(CS%v_av, "v2", restart_CS)
-  elseif ( (US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
-           (US%s_to_T_restart /= US%m_to_L_restart) ) then
-    vel_rescale = US%s_to_T_restart / US%m_to_L_restart
-    do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB ; CS%u_av(I,j,k) = vel_rescale * CS%u_av(I,j,k) ; enddo ; enddo ; enddo
-    do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; CS%v_av(i,J,k) = vel_rescale * CS%v_av(i,J,k) ; enddo ; enddo ; enddo
   endif
 
   if (CS%store_CAu) then
@@ -1504,15 +1477,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       if (.not. query_initialized(CS%h_av, "h2", restart_CS)) then
         CS%h_av(:,:,:) = h(:,:,:)
         call set_initialized(CS%h_av, "h2", restart_CS)
-      elseif ((GV%m_to_H_restart /= 0.0) .and. (GV%m_to_H_restart /= 1.0)) then
-        H_rescale = 1.0 / GV%m_to_H_restart
-        do k=1,nz ; do j=js,je ; do i=is,ie ; CS%h_av(i,j,k) = H_rescale * CS%h_av(i,j,k) ; enddo ; enddo ; enddo
-      endif
-      if ( (GV%m_to_H_restart * US%s_to_T_restart * US%m_to_L_restart /= 0.0) .and. &
-           (US%s_to_T_restart /= (GV%m_to_H_restart * US%m_to_L_restart**2)) ) then
-        uH_rescale = US%s_to_T_restart / (GV%m_to_H_restart * US%m_to_L_restart**2)
-        do k=1,nz ; do j=js,je ; do I=G%IscB,G%IecB ; uh(I,j,k) = uH_rescale * uh(I,j,k) ; enddo ; enddo ; enddo
-        do k=1,nz ; do J=G%JscB,G%JecB ; do i=is,ie ; vh(i,J,k) = uH_rescale * vh(i,J,k) ; enddo ; enddo ; enddo
       endif
     endif
   endif

src/core/MOM_verticalGrid.F90

@@ -12,7 +12,7 @@ module MOM_verticalGrid
 #include <MOM_memory.h>
 
 public verticalGridInit, verticalGridEnd
-public setVerticalGridAxes, fix_restart_scaling
+public setVerticalGridAxes
 public get_flux_units, get_thickness_units, get_tr_flux_units
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
@@ -75,7 +75,7 @@ module MOM_verticalGrid
   real :: H_to_MKS      !< A constant that translates thickness units to its MKS unit
                         !! (m or kg m-2) based on GV%Boussinesq [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
 
-  real :: m_to_H_restart = 0.0 !< A copy of the m_to_H that is used in restart files.
+  real :: m_to_H_restart = 1.0 !< A copy of the m_to_H that is used in restart files.
 end type verticalGrid_type
 
 contains
@@ -187,20 +187,6 @@ subroutine verticalGridInit( param_file, GV, US )
 
 end subroutine verticalGridInit
 
-!> Set the scaling factors for restart files to the scaling factors for this run.
-subroutine fix_restart_scaling(GV, unscaled)
-  type(verticalGrid_type), intent(inout) :: GV   !< The ocean's vertical grid structure
-  logical,       optional, intent(in)    :: unscaled !< If true, set the restart factors as though the
-                                             !! model would be unscaled, which is appropriate if the
-                                             !! scaling is undone when writing a restart file.
-
-  GV%m_to_H_restart = GV%m_to_H
-  if (present(unscaled)) then ; if (unscaled) then
-    GV%m_to_H_restart = 1.0
-  endif ; endif
-
-end subroutine fix_restart_scaling
-
 !> Returns the model's thickness units, usually m or kg/m^2.
 function get_thickness_units(GV)
   character(len=48)                   :: get_thickness_units !< The vertical thickness units

src/framework/MOM_unit_scaling.F90

@@ -55,12 +55,12 @@ module MOM_unit_scaling
   real :: Pa_to_RL2_T2    !< Convert pressures from Pa to R L2 T-2                                [R L2 T-2 Pa-1 ~> 1]
   real :: Pa_to_RLZ_T2    !< Convert wind stresses from Pa to R L Z T-2                          [R L Z T-2 Pa-1 ~> 1]
 
-  ! These are used for changing scaling across restarts.
-  real :: m_to_Z_restart = 0.0 !< A copy of the m_to_Z that is used in restart files.
-  real :: m_to_L_restart = 0.0 !< A copy of the m_to_L that is used in restart files.
-  real :: s_to_T_restart = 0.0 !< A copy of the s_to_T that is used in restart files.
-  real :: kg_m3_to_R_restart = 0.0 !< A copy of the kg_m3_to_R that is used in restart files.
-  real :: J_kg_to_Q_restart = 0.0 !< A copy of the J_kg_to_Q that is used in restart files.
+  ! These are no longer used for changing scaling across restarts.
+  real :: m_to_Z_restart = 1.0 !< A copy of the m_to_Z that is used in restart files.
+  real :: m_to_L_restart = 1.0 !< A copy of the m_to_L that is used in restart files.
+  real :: s_to_T_restart = 1.0 !< A copy of the s_to_T that is used in restart files.
+  real :: kg_m3_to_R_restart = 1.0 !< A copy of the kg_m3_to_R that is used in restart files.
+  real :: J_kg_to_Q_restart = 1.0 !< A copy of the J_kg_to_Q that is used in restart files.
 end type unit_scale_type
 
 contains
@@ -233,11 +233,11 @@ subroutine fix_restart_unit_scaling(US, unscaled)
                                              !! model would be unscaled, which is appropriate if the
                                              !! scaling is undone when writing a restart file.
 
-  US%m_to_Z_restart = US%m_to_Z
-  US%m_to_L_restart = US%m_to_L
-  US%s_to_T_restart = US%s_to_T
-  US%kg_m3_to_R_restart = US%kg_m3_to_R
-  US%J_kg_to_Q_restart = US%J_kg_to_Q
+  US%m_to_Z_restart = 1.0 ! US%m_to_Z
+  US%m_to_L_restart = 1.0 ! US%m_to_L
+  US%s_to_T_restart = 1.0 ! US%s_to_T
+  US%kg_m3_to_R_restart = 1.0 ! US%kg_m3_to_R
+  US%J_kg_to_Q_restart = 1.0 ! US%J_kg_to_Q
 
   if (present(unscaled)) then ; if (unscaled) then
     US%m_to_Z_restart = 1.0