+Rescale fluxes%net_mass_src and other diagnostics (#23)

Retaing dimensional rescaling for several diagnostics:
- Apply dimensional rescaling of fluxes%net_mass_src to and the net_mass_src
  argument to get_net_mass_forcing to [R Z T-1].
- Rescaled the patm argument of convert_state_to_ocean_type to [R L2 T-2]
  and press_to_Z to [Z T2 R-1 L-2]; these are not actually exercised yet,
  so this has a very limited scope, although three other local variables were
  also dimensionally rescaled.
- Revised the line breaks in two calls to register_restart to place the units
  and conversion factos on the same line, following a widespread code pattern.
- Used the scale argument in calls to global_area_integral or global_area_mean
  for 6 diagnostics, so that 3 other primary variables can be calculated in
  scaled units and rescaled via factors specified in the register_restart calls,
  following a widespread code pattern.
All answers and output are bitwise identical.

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -243,8 +243,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   real :: delta_sss           ! temporary storage for sss diff from restoring value [ppt]
   real :: delta_sst           ! temporary storage for sst diff from restoring value [degC]
 
-  real :: kg_m2_s_conversion        ! A combination of unit conversion factors for rescaling
-                              ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1].
+  real :: kg_m2_s_conversion  ! A combination of unit conversion factors for rescaling
+                              ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1]
   real :: rhoXcp              ! Reference density times heat capacity times unit scaling
                               ! factors [Q R degC-1 ~> J m-3 degC-1]
   real :: sign_for_net_FW_bug ! Should be +1. but an old bug can be recovered by using -1.
@@ -658,14 +658,16 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: &
     rigidity_at_h, &  ! Ice rigidity at tracer points [L4 Z-1 T-1 ~> m3 s-1]
-    net_mass_src, &   ! A temporary of net mass sources [kg m-2 s-1].
+    net_mass_src, &   ! A temporary of net mass sources [R Z T-1 ~> kg m-2 s-1].
     ustar_tmp         ! A temporary array of ustar values [Z T-1 ~> m s-1].
 
   real :: I_GEarth      ! The inverse of the gravitational acceleration [T2 Z L-2 ~> s2 m-1]
   real :: Kv_rho_ice    ! (CS%Kv_sea_ice / CS%density_sea_ice) [L4 Z-2 T-1 R-1 ~> m5 s-1 kg-1]
   real :: mass_ice      ! mass of sea ice at a face [R Z ~> kg m-2]
   real :: mass_eff      ! effective mass of sea ice for rigidity [R Z ~> kg m-2]
   real :: wt1, wt2      ! Relative weights of previous and current values of ustar [nondim].
+  real :: kg_m2_s_conversion  ! A combination of unit conversion factors for rescaling
+                              ! mass fluxes [R Z s m2 kg-1 T-1 ~> 1]
 
   integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
@@ -682,6 +684,8 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
   isr = is-isd+1 ; ier  = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1
   i0 = is - isc_bnd ; j0 = js - jsc_bnd
 
+  kg_m2_s_conversion = US%kg_m2s_to_RZ_T
+
   ! allocation and initialization if this is the first time that this
   ! mechanical forcing type has been used.
   if (.not.forces%initialized) then
@@ -774,15 +778,15 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
     i0 = is - isc_bnd ; j0 = js - jsc_bnd
     do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
       if (associated(IOB%lprec)) &
-        net_mass_src(i,j) = net_mass_src(i,j) + IOB%lprec(i-i0,j-j0)
+        net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * IOB%lprec(i-i0,j-j0)
       if (associated(IOB%fprec)) &
-        net_mass_src(i,j) = net_mass_src(i,j) + IOB%fprec(i-i0,j-j0)
+        net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * IOB%fprec(i-i0,j-j0)
       if (associated(IOB%runoff)) &
-        net_mass_src(i,j) = net_mass_src(i,j) + IOB%runoff(i-i0,j-j0)
+        net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * IOB%runoff(i-i0,j-j0)
       if (associated(IOB%calving)) &
-        net_mass_src(i,j) = net_mass_src(i,j) + IOB%calving(i-i0,j-j0)
+        net_mass_src(i,j) = net_mass_src(i,j) + kg_m2_s_conversion * IOB%calving(i-i0,j-j0)
       if (associated(IOB%q_flux)) &
-        net_mass_src(i,j) = net_mass_src(i,j) - IOB%q_flux(i-i0,j-j0)
+        net_mass_src(i,j) = net_mass_src(i,j) - kg_m2_s_conversion * IOB%q_flux(i-i0,j-j0)
     endif ; enddo ; enddo
     if (wt1 <= 0.0) then
       do j=js,je ; do i=is,ie

config_src/drivers/FMS_cap/ocean_model_MOM.F90

@@ -154,8 +154,8 @@ module ocean_model_mod
 
   logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the
                               !! ocean dynamics and forcing fluxes.
-  real :: press_to_z          !< A conversion factor between pressure and ocean
-                              !! depth in m, usually 1/(rho_0*g) [m Pa-1].
+  real :: press_to_z          !< A conversion factor between pressure and ocean depth,
+                              !! usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1].
   real :: C_p                 !< The heat capacity of seawater [J degC-1 kg-1].
   logical :: offline_tracer_mode = .false. !< If false, use the model in prognostic mode
                               !! with the barotropic and baroclinic dynamics, thermodynamics,
@@ -242,16 +242,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
   ! Local variables
-  real :: Rho0        ! The Boussinesq ocean density [kg m-3].
-  real :: G_Earth     ! The gravitational acceleration [m s-2].
-  real :: HFrz        !< If HFrz > 0 (m), melt potential will be computed.
+  real :: Rho0        ! The Boussinesq ocean density [R ~> kg m-3]
+  real :: G_Earth     ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
+  real :: HFrz        !< If HFrz > 0 [Z ~> m], melt potential will be computed.
                       !! The actual depth over which melt potential is computed will
                       !! min(HFrz, OBLD), where OBLD is the boundary layer depth.
                       !! If HFrz <= 0 (default), melt potential will not be computed.
-  logical :: use_melt_pot!< If true, allocate melt_potential array
+  logical :: use_melt_pot !< If true, allocate melt_potential array
 
-! 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 = "ocean_model_init"  ! This module's name.
   character(len=48)  :: stagger ! A string indicating the staggering locations for the
                                 ! surface velocities returned to the coupler.
@@ -331,28 +331,29 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, wind_stagger, gas
                  "calculate accelerations and the mass for conservation "//&
                  "properties, or with BOUSSINSEQ false to convert some "//&
                  "parameters from vertical units of m to kg m-2.", &
-                 units="kg m-3", default=1035.0)
+                 units="kg m-3", default=1035.0, scale=OS%US%kg_m3_to_R)
   call get_param(param_file, mdl, "G_EARTH", G_Earth, &
                  "The gravitational acceleration of the Earth.", &
-                 units="m s-2", default = 9.80)
+                 units="m s-2", default=9.80, scale=OS%US%m_s_to_L_T**2*OS%US%Z_to_m)
 
   call get_param(param_file, mdl, "ICE_SHELF",  OS%use_ice_shelf, &
                  "If true, enables the ice shelf model.", default=.false.)
 
   call get_param(param_file, mdl, "ICEBERGS_APPLY_RIGID_BOUNDARY",  OS%icebergs_alter_ocean, &
                  "If true, allows icebergs to change boundary condition felt by ocean", default=.false.)
 
-  OS%press_to_z = 1.0/(Rho0*G_Earth)
+  OS%press_to_z = 1.0 / (Rho0*G_Earth)
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
   call get_param(param_file, mdl, "HFREEZE", HFrz, &
                  "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
                  "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
                  "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
-                 "melt potential will not be computed.", units="m", default=-1.0, do_not_log=.true.)
+                 "melt potential will not be computed.", &
+                 units="m", default=-1.0, scale=OS%US%m_to_Z, do_not_log=.true.)
 
-  if (HFrz .gt. 0.0) then
+  if (HFrz > 0.0) then
     use_melt_pot=.true.
   else
     use_melt_pot=.false.
@@ -655,7 +656,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
 
 ! Translate state into Ocean.
 !  call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US, &
-!                                   Ice_ocean_boundary%p, OS%press_to_z)
+!                                   OS%fluxes%p_surf_full, OS%press_to_z)
   call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, OS%US)
   Time1 = OS%Time ; if (do_dyn) Time1 = OS%Time_dyn
   call coupler_type_send_data(Ocean_sfc%fields, Time1)
@@ -816,9 +817,9 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
                                                !! have their data set here.
   type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
   type(unit_scale_type), intent(in)    :: US   !< A dimensional unit scaling type
-  real,        optional, intent(in)    :: patm(:,:)  !< The pressure at the ocean surface [Pa].
-  real,        optional, intent(in)    :: press_to_z !< A conversion factor between pressure and
-                                               !! ocean depth in m, usually 1/(rho_0*g) [m Pa-1].
+  real,        optional, intent(in)    :: patm(:,:)  !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
+  real,        optional, intent(in)    :: press_to_z !< A conversion factor between pressure and ocean
+                                               !! depth, usually 1/(rho_0*g) [Z T2 R-1 L-2 ~> m Pa-1]
   ! Local variables
   real :: IgR0
   character(len=48)  :: val_str
@@ -860,7 +861,7 @@ subroutine convert_state_to_ocean_type(sfc_state, Ocean_sfc, G, US, patm, press_
 
   if (present(patm)) then
     do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
-      Ocean_sfc%sea_lev(i,j) = US%Z_to_m * sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z
+      Ocean_sfc%sea_lev(i,j) = US%Z_to_m * (sfc_state%sea_lev(i+i0,j+j0) + patm(i,j) * press_to_z)
       Ocean_sfc%area(i,j) = US%L_to_m**2 * G%areaT(i+i0,j+j0)
     enddo ; enddo
   else