Rescale 6 ice shelf variables

  Changed the rescaling of 6 ice shelf variables to cancel out common
conversion factors that appear in several expressions.

  The C_basal_friction argument to initialize_ice_C_basal_friction is now in
partially rescaled units, reflecting the portion that does not cancel out
fractional-power units from a power law fit with an arbitrary power.  The
C_basal_friction array in ice_shelf_dyn_CS and the C_friction variable in
initialize_ice_C_basal_friction were similarly rescaled.  There are new
scale factors in a get_param_call and a MOM_read_data call and a conversion
factor in register_restart_field call that reflect these changes.

  The KE_tot and mass_tot variables in write_ice_shelf_energy,
are kept in scaled units until they are written.

  The internal variable fN in calc_shelf_taub is kept in scaled units, but there
is now a scaling factor of US%Z_to_L in the expression for fB.  This latter
could be folded into the CF_Max element of ice_shelf_dyn_CS, but I am unsure
whether this would be physically sensible.

  All answers should be bitwise identical and no output should change, but this
has not been extensively tested yet.

src/ice_shelf/MOM_ice_shelf_dynamics.F90

@@ -111,7 +111,7 @@ module MOM_ice_shelf_dynamics
                                                             !! of "linearized" basal stress (Pa) [R Z L2 T-1 ~> kg s-1]
                 !!  The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011
   real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric),
-                                                            !!  units= [Pa (s m-1)^(n_basal_fric)]
+                               !! units of [R L Z T-2 (s m-1)^(n_basal_fric) ~> Pa (s m-1)^(n_basal_fric)]
   real, pointer, dimension(:,:) :: OD_rt => NULL()         !< A running total for calculating OD_av [Z ~> m].
   real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac.
   real, pointer, dimension(:,:) :: OD_av => NULL()         !< The time average open ocean depth [Z ~> m].
@@ -177,7 +177,7 @@ module MOM_ice_shelf_dynamics
   real :: eps_glen_min      !< Min. strain rate to avoid infinite Glen's law viscosity, [T-1 ~> s-1].
   real :: n_basal_fric      !< Exponent in sliding law tau_b = C u^(m_slide) [nondim]
   logical :: CoulombFriction !< Use Coulomb friction law (Schoof 2005, Gagliardini et al 2007)
-  real :: CF_MinN           !< Minimum Coulomb friction effective pressure [Pa]
+  real :: CF_MinN           !< Minimum Coulomb friction effective pressure [R L2 T-2 ~> Pa]
   real :: CF_PostPeak       !< Coulomb friction post peak exponent [nondim]
   real :: CF_Max            !< Coulomb friction maximum coefficient [nondim]
   real :: density_ocean_avg !< A typical ocean density [R ~> kg m-3].  This does not affect ocean
@@ -359,7 +359,8 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
     allocate(CS%ice_visc(isd:ied,jsd:jed,CS%visc_qps), source=0.0)
     allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1]
     allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0)   ! [R Z L2 T-1 ~> kg s-1]
-    allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (s m-1)^n_sliding]
+    allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10*US%Pa_to_RLZ_T2)
+             ! Units of [R L Z T-2 (s m-1)^n_sliding ~> Pa (s m-1)^n_sliding]
     allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0)
     allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0)
     allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0)
@@ -396,7 +397,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS)
     call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, &
                                 "fractional degree of grounding", "nondim")
     call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, &
-                                "basal sliding coefficients", "Pa (s m-1)^n_sliding")
+                                "basal sliding coefficients", "Pa (s m-1)^n_sliding", conversion=US%RLZ_T2_to_Pa)
     call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, &
                                 "ice-stiffness parameter", "Pa-3 s-1")
     call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, &
@@ -536,7 +537,7 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_
                  units="none", default=.false., fail_if_missing=.false.)
     call get_param(param_file, mdl, "CF_MinN", CS%CF_MinN, &
                  "Minimum Coulomb friction effective pressure", &
-                 units="Pa", default=1.0, fail_if_missing=.false.)
+                 units="Pa", default=1.0, scale=US%Pa_to_RLZ_T2, fail_if_missing=.false.)
     call get_param(param_file, mdl, "CF_PostPeak", CS%CF_PostPeak, &
                  "Coulomb friction post peak exponent", &
                  units="none", default=1.0, fail_if_missing=.false.)
@@ -1192,8 +1193,8 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step)
   ! Local variables
   type(time_type) :: dt ! A time_type version of the timestep.
   real, dimension(SZDI_(G),SZDJ_(G)) :: tmp1 ! A temporary array used in reproducing sums [various]
-  real :: KE_tot    ! THe total kinetic energy [J]
-  real :: mass_tot  ! The total mass [kg]
+  real :: KE_tot    ! The total kinetic energy [R Z L4 T-2 ~> J]
+  real :: mass_tot  ! The total mass [R Z L2 ~> kg]
   integer :: is, ie, js, je, isr, ier, jsr, jer, i, j
   character(len=32)  :: mesg_intro, time_units, day_str, n_str, date_str
   integer :: start_of_day, num_days
@@ -1250,16 +1251,15 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step)
        (((CS%v_shelf(I-1,J-1)+CS%v_shelf(I,J))+(CS%v_shelf(I,J-1)+CS%v_shelf(I-1,J)))**2))
   enddo; enddo
 
-  KE_tot = US%RZL2_to_kg*US%L_T_to_m_s**2 * &
-           reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=(US%RZL2_to_kg*US%L_T_to_m_s**2))
+  KE_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=(US%RZL2_to_kg*US%L_T_to_m_s**2))
 
   !calculate mass
   tmp1(:,:) = 0.0
   do j=js,je ; do i=is,ie
     tmp1(i,j) = mass(i,j) * area(i,j)
   enddo; enddo
 
-  mass_tot = US%RZL2_to_kg * reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=US%RZL2_to_kg)
+  mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, unscale=US%RZL2_to_kg)
 
   if (is_root_pe()) then  ! Only the root PE actually writes anything.
     if (day > CS%Start_time) then
@@ -1305,7 +1305,7 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step)
     else                        ; write(n_str, '(I10)') CS%prev_IS_energy_calls ; endif
 
     write(CS%IS_fileenergy_ascii,'(A,",",A,", En ",ES22.16,", M ",ES11.5)') &
-      trim(n_str), trim(day_str), KE_tot/mass_tot, mass_tot
+      trim(n_str), trim(day_str), US%L_T_to_m_s**2*KE_tot/mass_tot, US%RZL2_to_kg*mass_tot
   endif
 
   CS%prev_IS_energy_calls = CS%prev_IS_energy_calls + 1
@@ -3380,11 +3380,10 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
   real :: umid, vmid ! Velocities [L T-1 ~> m s-1]
   real :: eps_min ! A minimal strain rate used in the Glens flow law expression [T-1 ~> s-1]
   real :: unorm ! The magnitude of the velocity in mks units for use with fractional powers [m s-1]
-  real :: alpha !Coulomb coefficient [nondim]
+  real :: alpha ! Coulomb coefficient [nondim]
   real :: Hf !"floatation thickness" for Coulomb friction [Z ~> m]
-  real :: fN !Effective pressure (ice pressure - ocean pressure) for Coulomb friction [Pa]
+  real :: fN ! Effective pressure (ice pressure - ocean pressure) for Coulomb friction [R L2 T-2 ~> Pa]
   real :: fB !for Coulomb Friction [(T L-1)^CS%CF_PostPeak ~> (s m-1)^CS%CF_PostPeak]
-  real :: fN_scale !To convert effective pressure to mks units during Coulomb friction [Pa T2 R-1 L-2 ~> 1]
 
   isc = G%isc ; jsc = G%jsc ; iec = G%iec ; jec = G%jec
   iscq = G%iscB ; iecq = G%iecB ; jscq = G%jscB ; jecq = G%jecB
@@ -3402,7 +3401,6 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
     else
       alpha = 1.0
     endif
-    fN_scale = US%R_to_kg_m3 * US%L_T_to_m_s**2  ! Unscaling factor for use in the fractional power law.
   endif
 
   do j=jsd+1,jed
@@ -3416,16 +3414,16 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf)
         if (CS%CoulombFriction) then
           !Effective pressure
           Hf = max((CS%density_ocean_avg/CS%density_ice) * CS%bed_elev(i,j), 0.0)
-          fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)),CS%CF_MinN)
-          fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric)
+          fN = max(((CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)), CS%CF_MinN)
+          fB = alpha * (US%Z_to_L*CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric)
 
           CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * &
               (unorm**(CS%n_basal_fric-1.0) / (1.0 + fB * unorm**CS%CF_PostPeak)**(CS%n_basal_fric))) * &
-              (US%Pa_to_RLZ_T2*US%L_T_to_m_s)   ! Restore the scaling after the fractional power law.
+              US%L_T_to_m_s   ! Restore the scaling after the fractional power law.
         else
           !linear (CS%n_basal_fric=1) or "Weertman"/power-law (CS%n_basal_fric /= 1)
           CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * &
-                                   (US%Pa_to_RLZ_T2*US%L_T_to_m_s) ! Rescale after the fractional power law.
+                                   US%L_T_to_m_s ! Rescale after the fractional power law.
         endif
 
         CS%basal_traction(i,j)=max(CS%basal_traction(i,j), CS%min_basal_traction * G%areaT(i,j))

src/ice_shelf/MOM_ice_shelf_initialize.F90

@@ -556,12 +556,14 @@ end subroutine initialize_ice_shelf_boundary_from_file
 subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF)
   type(ocean_grid_type), intent(in)    :: G    !< The ocean's grid structure
   real, dimension(SZDI_(G),SZDJ_(G)), &
-                         intent(inout) :: C_basal_friction !< Ice-stream basal friction [Pa (s m-1)^(n_basal_fric)]
+                         intent(inout) :: C_basal_friction !< Ice-stream basal friction
+                                             !! in units of [R L Z T-2 (s m-1)^n_basal_fric ~> Pa (s m-1)^n_basal_fric]
   type(unit_scale_type), intent(in)    :: US !< A structure containing unit conversion factors
   type(param_file_type), intent(in)    :: PF !< A structure to parse for run-time parameters
 
 !  integer :: i, j
-  real :: C_friction
+  real :: C_friction  ! Constant ice-stream basal friction in units of
+                      ! [R L Z T-2 (s m-1)^n_basal_fric ~> Pa (s m-1)^n_basal_fric]
   character(len=40)  :: mdl = "initialize_ice_basal_friction" ! This subroutine's name.
   character(len=200) :: config
   character(len=200) :: varname
@@ -574,7 +576,7 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF)
 
   if (trim(config)=="CONSTANT") then
     call get_param(PF, mdl, "BASAL_FRICTION_COEFF", C_friction, &
-                 "Coefficient in sliding law.", units="Pa (s m-1)^(n_basal_fric)", default=5.e10)
+                 "Coefficient in sliding law.", units="Pa (s m-1)^(n_basal_fric)", default=5.e10, scale=US%RLZ_T2_to_Pa)
 
     C_basal_friction(:,:) = C_friction
   elseif (trim(config)=="FILE") then
@@ -595,7 +597,7 @@ subroutine initialize_ice_C_basal_friction(C_basal_friction, G, US, PF)
     if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
        " initialize_ice_basal_friction_from_file: Unable to open "//trim(filename))
 
-    call MOM_read_data(filename,trim(varname),C_basal_friction,G%Domain)
+    call MOM_read_data(filename, trim(varname), C_basal_friction, G%Domain, scale=US%RLZ_T2_to_Pa)
 
   endif
 end subroutine