+Changed units of GV%Rho0 to [R]

  Changed the units of GV%Rho0 from [kg m-3] to [R] for dimensional consistency
testing.  This required the addition of unit_scale_type arguments to several
interfaces.  All answers are bitwise identical, but new arguments have been
added to several public interfaces and the units of an element in a public
type have changed.

config_src/mct_driver/mom_ocean_model_mct.F90

@@ -582,7 +582,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 #endif
   endif
 
-  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0)
+  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%US%R_to_kg_m3*OS%GV%Rho0)
   call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid)
 
   if (OS%use_waves) then

config_src/nuopc_driver/mom_ocean_model_nuopc.F90

@@ -570,7 +570,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
     call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average
 #endif
   endif
-  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0)
+  call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%US%R_to_kg_m3*OS%GV%Rho0)
   call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid)
 
   if (OS%use_waves) then

src/core/MOM.F90

@@ -2678,7 +2678,7 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
         call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, &
                                Rho_conv, tv%eqn_of_state)
       else
-        Rho_conv=GV%Rho0
+        Rho_conv = US%R_to_kg_m3*GV%Rho0
       endif
       IgR0 = 1.0 / (Rho_conv * GV%mks_g_Earth)
       ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0
@@ -2914,7 +2914,7 @@ subroutine extract_surface_state(CS, sfc_state)
 
        if (G%mask2dT(i,j)>0.) then
          ! instantaneous melt_potential [J m-2]
-         sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%GV%Rho0 * delT(i)
+         sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%US%R_to_kg_m3*GV%Rho0 * delT(i)
        endif
       enddo
     enddo ! end of j loop

src/core/MOM_PressureForce_Montgomery.F90

@@ -435,7 +435,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
 
   h_neglect = GV%H_subroundoff * GV%H_to_Z
   I_Rho0 = 1.0/CS%Rho0
-  G_Rho0 = GV%g_Earth/GV%Rho0
+  G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
 
   if (CS%tides) then
     !   Determine the surface height anomaly for calculating self attraction
@@ -640,7 +640,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
 
   Rho0xG = Rho0*US%L_T_to_m_s**2 * GV%g_Earth
-  G_Rho0 = GV%g_Earth / GV%Rho0
+  G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   use_EOS = associated(tv%eqn_of_state)
   z_neglect = GV%H_subroundoff*GV%H_to_Z
 

src/core/MOM_PressureForce_analytic_FV.F90

@@ -531,9 +531,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
 
   h_neglect = GV%H_subroundoff
   dz_neglect = GV%H_subroundoff * GV%H_to_Z
-  I_Rho0 = US%m_s_to_L_T**2 / GV%Rho0
+  I_Rho0 = US%m_s_to_L_T**2 / (US%R_to_kg_m3*GV%Rho0)
   g_Earth_z = US%L_T_to_m_s**2 * GV%g_Earth
-  G_Rho0 = GV%g_Earth/GV%Rho0
+  G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   rho_ref = CS%Rho0
 
   if (CS%tides) then

src/core/MOM_PressureForce_blocked_AFV.F90

@@ -515,9 +515,9 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
 
   h_neglect = GV%H_subroundoff
   dz_neglect = GV%H_subroundoff * GV%H_to_Z
-  I_Rho0 = US%m_s_to_L_T**2 / GV%Rho0
+  I_Rho0 = US%m_s_to_L_T**2 / (US%R_to_kg_m3*GV%Rho0)
   g_Earth_z = US%L_T_to_m_s**2 * GV%g_Earth
-  G_Rho0 = GV%g_Earth / GV%Rho0
+  G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   rho_ref = CS%Rho0
 
   if (CS%tides) then

src/core/MOM_barotropic.F90

@@ -724,8 +724,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
   dtbt = dt_in_T * Instep
   bebt = CS%bebt
   be_proj = CS%bebt
-  mass_accel_to_Z = US%m_to_L*US%T_to_s**2 * US%m_to_Z / GV%Rho0
-  mass_to_Z = US%m_to_Z / GV%Rho0
+  mass_accel_to_Z = US%m_to_L*US%T_to_s**2 * US%m_to_Z / (US%R_to_kg_m3*GV%Rho0)
+  mass_to_Z = US%m_to_Z / (US%R_to_kg_m3*GV%Rho0)
 
   !--- setup the weight when computing vbt_trans and ubt_trans
   if (project_velocity) then
@@ -4345,10 +4345,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
     enddo ; enddo
 ! else
 !   do j=js,je ; do I=is-1,ie
-!     CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (GV%Rho0*(G%bathyT(i+1,j) + G%bathyT(i,j)))
+!     CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (US%R_to_kg_m3*GV%Rho0*(G%bathyT(i+1,j) + G%bathyT(i,j)))
 !   enddo ; enddo
 !   do J=js-1,je ; do i=is,ie
-!     CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (GV%Rho0*(G%bathyT(i,j+1) + G%bathyT(i,j)))
+!     CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (US%R_to_kg_m3*GV%Rho0*(G%bathyT(i,j+1) + G%bathyT(i,j)))
 !   enddo ; enddo
 ! endif
 

src/core/MOM_forcing_type.F90

@@ -336,14 +336,15 @@ module MOM_forcing_type
 !! for optimization purposes. The 2d (i,j) wrapper is the next subroutine below.
 !! This routine multiplies fluxes by dt, so that the result is an accumulation of fluxes
 !! over a time step.
-subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,                           &
+subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,               &
                   FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent, &
                   h, T, netMassInOut, netMassOut, net_heat, net_salt, pen_SW_bnd, tv,   &
                   aggregate_FW, nonpenSW, netmassInOut_rate,net_Heat_Rate,      &
                   net_salt_rate, pen_sw_bnd_Rate, skip_diags)
 
   type(ocean_grid_type),    intent(in)    :: G              !< ocean grid structure
   type(verticalGrid_type),  intent(in)    :: GV             !< ocean vertical grid structure
+  type(unit_scale_type),    intent(in)    :: US             !< A dimensional unit scaling type
   type(forcing),            intent(inout) :: fluxes         !< structure containing pointers to possible
                                                             !! forcing fields. NULL unused fields.
   type(optics_type),        pointer       :: optics         !< pointer to optics
@@ -433,7 +434,7 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
   !}BGR
 
   Ih_limit  = 1.0 / FluxRescaleDepth
-  Irho0     = 1.0 / GV%Rho0
+  Irho0     = 1.0 / (US%R_to_kg_m3*GV%Rho0)
   I_Cp      = 1.0 / fluxes%C_p
   J_m2_to_H = 1.0 / (GV%H_to_kg_m2 * fluxes%C_p)
 
@@ -804,13 +805,14 @@ end subroutine extractFluxes1d
 !> 2d wrapper for 1d extract fluxes from surface fluxes type.
 !! This subroutine extracts fluxes from the surface fluxes type. It multiplies the
 !! fluxes by dt, so that the result is an accumulation of the fluxes over a time step.
-subroutine extractFluxes2d(G, GV, fluxes, optics, nsw, dt, FluxRescaleDepth, &
+subroutine extractFluxes2d(G, GV, US, fluxes, optics, nsw, dt, FluxRescaleDepth, &
                            useRiverHeatContent, useCalvingHeatContent, h, T, &
                            netMassInOut, netMassOut, net_heat, Net_salt, Pen_SW_bnd, tv, &
                            aggregate_FW)
 
   type(ocean_grid_type),            intent(in)    :: G              !< ocean grid structure
   type(verticalGrid_type),          intent(in)    :: GV             !< ocean vertical grid structure
+  type(unit_scale_type),            intent(in)    :: US             !< A dimensional unit scaling type
   type(forcing),                    intent(inout) :: fluxes         !< structure containing pointers to forcing.
   type(optics_type),                pointer       :: optics         !< pointer to optics
   integer,                          intent(in)    :: nsw            !< number of bands of penetrating SW
@@ -854,7 +856,7 @@ subroutine extractFluxes2d(G, GV, fluxes, optics, nsw, dt, FluxRescaleDepth, &
 !$OMP                                  h,T,netMassInOut,netMassOut,Net_heat,Net_salt,Pen_SW_bnd,tv, &
 !$OMP                                  aggregate_FW)
   do j=G%jsc, G%jec
-    call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,                      &
+    call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,                      &
             FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent,&
             h(:,j,:), T(:,j,:), netMassInOut(:,j), netMassOut(:,j),              &
             net_heat(:,j), net_salt(:,j), pen_SW_bnd(:,:,j), tv, aggregate_FW)
@@ -916,7 +918,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
 
   depthBeforeScalingFluxes = max( GV%Angstrom_H, 1.e-30*GV%m_to_H )
   pressure(:) = 0. ! Ignore atmospheric pressure
-  GoRho       = (GV%g_Earth*US%m_to_Z * GV%H_to_m*US%T_to_s) / GV%Rho0
+  GoRho       = (GV%g_Earth*US%m_to_Z * GV%H_to_m*US%T_to_s) / (US%R_to_kg_m3*GV%Rho0)
   start       = 1 + G%isc - G%isd
   npts        = 1 + G%iec - G%isc
 
@@ -929,7 +931,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
   ! netSalt    = salt via surface fluxes [ppt H s-1 ~> ppt m s-1 or gSalt m-2 s-1]
   ! Note that unlike other calls to extractFLuxes1d() that return the time-integrated flux
   ! this call returns the rate because dt=1
-  call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,                                 &
+  call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt,                         &
                 depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, &
                 h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW,                 &
                 netSalt, penSWbnd, tv, .false., skip_diags=skip_diags)

src/core/MOM_isopycnal_slopes.F90

@@ -121,7 +121,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
 
   present_N2_u = PRESENT(N2_u)
   present_N2_v = PRESENT(N2_v)
-  G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / GV%Rho0
+  G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0)
   if (present_N2_u) then
     do j=js,je ; do I=is-1,ie
       N2_u(I,j,1) = 0.

src/core/MOM_verticalGrid.F90

@@ -101,7 +101,7 @@ subroutine verticalGridInit( param_file, GV, US )
                  "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=US%kg_m3_to_R)
   call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, &
                  "If true, make the Boussinesq approximation.", default=.true.)
   call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_m, &
@@ -143,15 +143,15 @@ subroutine verticalGridInit( param_file, GV, US )
   GV%ke = nk
 
   if (GV%Boussinesq) then
-    GV%H_to_kg_m2 = GV%Rho0 * GV%H_to_m
+    GV%H_to_kg_m2 = US%R_to_kg_m3*GV%Rho0 * GV%H_to_m
     GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2
     GV%m_to_H = 1.0 / GV%H_to_m
     GV%Angstrom_H = GV%m_to_H * GV%Angstrom_m
     GV%H_to_MKS = GV%H_to_m
   else
     GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2
-    GV%m_to_H = GV%Rho0 * GV%kg_m2_to_H
-    GV%H_to_m = GV%H_to_kg_m2 / GV%Rho0
+    GV%m_to_H = US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H
+    GV%H_to_m = GV%H_to_kg_m2 / (US%R_to_kg_m3*GV%Rho0)
     GV%Angstrom_H = GV%Angstrom_m*1000.0*GV%kg_m2_to_H
     GV%H_to_MKS = GV%H_to_kg_m2
   endif

src/diagnostics/MOM_diagnostics.F90

@@ -842,7 +842,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS)
             z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k)
           enddo ; enddo
           call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), &
-                              z_top, z_bot, 0.0, GV%Rho0, GV%mks_g_Earth*US%Z_to_m, &
+                              z_top, z_bot, 0.0, US%R_to_kg_m3*GV%Rho0, GV%mks_g_Earth*US%Z_to_m, &
                               G%HI, G%HI, tv%eqn_of_state, dpress)
           do j=js,je ; do i=is,ie
             mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth
@@ -2006,7 +2006,7 @@ subroutine write_static_fields(G, GV, US, tv, diag)
 
   id = register_static_field('ocean_model','Rho_0', diag%axesNull, &
        'mean ocean density used with the Boussinesq approximation', &
-       'kg m-3', cmor_field_name='rhozero', &
+       'kg m-3', cmor_field_name='rhozero', conversion=US%R_to_kg_m3, &
        cmor_standard_name='reference_sea_water_density_for_boussinesq_approximation', &
        cmor_long_name='reference sea water density for boussinesq approximation')
   if (id > 0) call post_data(id, GV%Rho0, diag, .true.)

src/diagnostics/MOM_sum_output.F90

@@ -542,7 +542,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
         tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j)
       enddo ; enddo ; enddo
       mass_tot = reproducing_sum(tmp1, sums=mass_lay, EFP_sum=mass_EFP)
-      do k=1,nz ; vol_lay(k) = US%m_to_Z * (mass_lay(k) / GV%Rho0) ; enddo
+      do k=1,nz ; vol_lay(k) = US%m_to_Z * (mass_lay(k) / (US%R_to_kg_m3*GV%Rho0)) ; enddo
     endif
   endif ! Boussinesq
 
@@ -666,7 +666,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
           hint = Z_0APE(K) + (hbelow - G%bathyT(i,j))
           hbot = Z_0APE(K) - G%bathyT(i,j)
           hbot = (hbot + ABS(hbot)) * 0.5
-          PE_pt(i,j,K) = 0.5 * areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(GV%Rho0*GV%g_prime(K)) * &
+          PE_pt(i,j,K) = 0.5 * areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(US%R_to_kg_m3*GV%Rho0*GV%g_prime(K)) * &
                   (hint * hint - hbot * hbot)
         enddo
       enddo ; enddo
@@ -675,7 +675,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
         do k=nz,1,-1
           hint = Z_0APE(K) + eta(i,j,K)  ! eta and H_0 have opposite signs.
           hbot = max(Z_0APE(K) - G%bathyT(i,j), 0.0)
-          PE_pt(i,j,K) = 0.5 * (areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(GV%Rho0*GV%g_prime(K))) * &
+          PE_pt(i,j,K) = 0.5 * (areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(US%R_to_kg_m3*GV%Rho0*GV%g_prime(K))) * &
                   (hint * hint - hbot * hbot)
         enddo
       enddo ; enddo
@@ -750,7 +750,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
     CS%salt_prev_EFP = salt_EFP ; CS%net_salt_in_EFP = real_to_EFP(0.0)
     CS%heat_prev_EFP = heat_EFP ; CS%net_heat_in_EFP = real_to_EFP(0.0)
   endif
-  Irho0 = 1.0/GV%Rho0
+  Irho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0)
 
   if (CS%use_temperature) then
     Salt_chg_EFP = Salt_EFP - CS%salt_prev_EFP

src/diagnostics/MOM_wave_speed.F90

@@ -132,7 +132,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, &
   endif
 
   S => tv%S ; T => tv%T
-  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / GV%Rho0
+  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   Z_to_Pa = GV%Z_to_H * GV%H_to_Pa
   use_EOS = associated(tv%eqn_of_state)
 
@@ -600,7 +600,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
   endif ; endif
 
   S => tv%S ; T => tv%T
-  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / GV%Rho0
+  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   use_EOS = associated(tv%eqn_of_state)
   Z_to_Pa = GV%Z_to_H * GV%H_to_Pa
 

src/diagnostics/MOM_wave_structure.F90

@@ -178,7 +178,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   Pi = (4.0*atan(1.0))
 
   S => tv%S ; T => tv%T
-  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth /GV%Rho0
+  g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0)
   cg_subRO = 1e-100*US%m_s_to_L_T  ! The hard-coded value here might need to increase.
   use_EOS = associated(tv%eqn_of_state)
 
@@ -479,8 +479,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
             ! Back-calculate amplitude from energy equation
             if (Kmag2 > 0.0) then
               !### This should be simpified to use a single division.
-              KE_term = 0.25*GV%Rho0*( ((1.0 + f2/freq**2) / Kmag2)*int_dwdz2 + int_w2 )
-              PE_term = 0.25*GV%Rho0*( int_N2w2/(US%s_to_T*freq)**2 )
+              KE_term = 0.25*US%R_to_kg_m3*GV%Rho0*( ((1.0 + f2/freq**2) / Kmag2)*int_dwdz2 + int_w2 )
+              PE_term = 0.25*US%R_to_kg_m3*GV%Rho0*( int_N2w2/(US%s_to_T*freq)**2 )
               if (En(i,j) >= 0.0) then
                 W0 = sqrt( En(i,j)/(KE_term + PE_term) )
               else