+Add rescaling for temperature and salinity (2)

  This commit completes the dimensional rescaling for temperature and salinity,
and it has been confirmed that the solutions for the existing test cases pass
these tests.  There are new unit_scale_type arguments to several publicly
visible interfaces, mostly related to temperature initialization.  There is also
a new optional argument, conc_scale to the register_tracer calls to specify the
conversion that should be done to tracer concentrations during output.
Additionally, there are new entries in the incorrect units on some runtime
parameters for user code were corrected in the MOM_parameter_doc files for some
test cases.  All answers are bitwise identical in the MOM6 regression suite,
including when the temperature and salinity rescaling are enabled.

src/core/MOM_variables.F90

@@ -37,7 +37,7 @@ module MOM_variables
 end type p2d
 
 !> Pointers to various fields which may be used describe the surface state of MOM, and which
-!! will be returned to a the calling program
+!! will be returned to the calling program
 type, public :: surface
   real, allocatable, dimension(:,:) :: &
     SST, &         !< The sea surface temperature [degC].
@@ -81,24 +81,23 @@ module MOM_variables
 !! potential temperature, salinity, heat capacity, and the equation of state control structure.
 type, public :: thermo_var_ptrs
   ! If allocated, the following variables have nz layers.
-  real, pointer :: T(:,:,:) => NULL() !< Potential temperature [degC].
-  real, pointer :: S(:,:,:) => NULL() !< Salinity [PSU] or [gSalt/kg], generically [ppt].
+  real, pointer :: T(:,:,:) => NULL() !< Potential temperature [C ~> degC].
+  real, pointer :: S(:,:,:) => NULL() !< Salinity [PSU] or [gSalt/kg], generically [S ~> ppt].
   real, pointer :: p_surf(:,:) => NULL() !< Ocean surface pressure used in equation of state
                          !! calculations [R L2 T-2 ~> Pa]
   type(EOS_type), pointer :: eqn_of_state => NULL() !< Type that indicates the
                          !! equation of state to use.
   real :: P_Ref          !<   The coordinate-density reference pressure [R L2 T-2 ~> Pa].
                          !! This is the pressure used to calculate Rml from
                          !! T and S when eqn_of_state is associated.
-  real :: C_p            !<   The heat capacity of seawater [Q degC-1 ~> J degC-1 kg-1].
+  real :: C_p            !<   The heat capacity of seawater [Q C-1 ~> J degC-1 kg-1].
                          !! When conservative temperature is used, this is
                          !! constant and exactly 3991.86795711963 J degC-1 kg-1.
   logical :: T_is_conT = .false. !< If true, the temperature variable tv%T is
                          !! actually the conservative temperature [degC].
   logical :: S_is_absS = .false. !< If true, the salinity variable tv%S is
                          !! actually the absolute salinity in units of [gSalt kg-1].
-  real :: min_salinity = 0.01 !< The minimum value of salinity when BOUND_SALINITY=True [ppt].
-                         !! The default is 0.01 for backward compatibility but should be 0.
+  real :: min_salinity   !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt].
   ! These arrays are accumulated fluxes for communication with other components.
   real, dimension(:,:), pointer :: frazil => NULL()
                          !< The energy needed to heat the ocean column to the
@@ -111,19 +110,19 @@ module MOM_variables
   real, dimension(:,:), pointer :: TempxPmE => NULL()
                          !<   The net inflow of water into the ocean times the
                          !! temperature at which this inflow occurs since the
-                         !! last call to calculate_surface_state [degC R Z ~> degC kg m-2].
+                         !! last call to calculate_surface_state [C R Z ~> degC kg m-2].
                          !! This should be prescribed in the forcing fields, but
                          !! as it often is not, this is a useful heat budget diagnostic.
   real, dimension(:,:), pointer :: internal_heat => NULL()
                          !< Any internal or geothermal heat sources that
                          !! have been applied to the ocean since the last call to
-                         !! calculate_surface_state [degC R Z ~> degC kg m-2].
+                         !! calculate_surface_state [C R Z ~> degC kg m-2].
   ! The following variables are most normally not used but when they are they
   ! will be either set by parameterizations or prognostic.
-  real, pointer :: varT(:,:,:) => NULL() !< SGS variance of potential temperature [degC2].
-  real, pointer :: varS(:,:,:) => NULL() !< SGS variance of salinity [ppt2].
+  real, pointer :: varT(:,:,:) => NULL() !< SGS variance of potential temperature [C2 ~> degC2].
+  real, pointer :: varS(:,:,:) => NULL() !< SGS variance of salinity [S2 ~> ppt2].
   real, pointer :: covarTS(:,:,:) => NULL() !< SGS covariance of salinity and potential
-                                  !! temperature [degC ppt].
+                                  !! temperature [C S ~> degC ppt].
 end type thermo_var_ptrs
 
 !> Pointers to all of the prognostic variables allocated in MOM_variables.F90 and MOM.F90.
@@ -133,8 +132,8 @@ module MOM_variables
 !! they refer to in MOM.F90.
 type, public :: ocean_internal_state
   real, pointer, dimension(:,:,:) :: &
-    T => NULL(), & !< Pointer to the temperature state variable [degC]
-    S => NULL(), & !< Pointer to the salinity state variable [ppt ~> PSU or g/kg]
+    T => NULL(), & !< Pointer to the temperature state variable [C ~> degC]
+    S => NULL(), & !< Pointer to the salinity state variable [S ~> ppt] (i.e., PSU or g/kg)
     u => NULL(), & !< Pointer to the zonal velocity [L T-1 ~> m s-1]
     v => NULL(), & !< Pointer to the meridional velocity [L T-1 ~> m s-1]
     h => NULL()    !< Pointer to the layer thicknesses [H ~> m or kg m-2]
@@ -580,15 +579,15 @@ subroutine MOM_thermovar_chksum(mesg, tv, G, US)
   ! counts, there must be no redundant points, so all variables use is..ie
   ! and js...je as their extent.
   if (associated(tv%T)) &
-    call hchksum(tv%T, mesg//" tv%T", G%HI)
+    call hchksum(tv%T, mesg//" tv%T", G%HI, scale=US%C_to_degC)
   if (associated(tv%S)) &
-    call hchksum(tv%S, mesg//" tv%S", G%HI)
+    call hchksum(tv%S, mesg//" tv%S", G%HI, scale=US%S_to_ppt)
   if (associated(tv%frazil)) &
     call hchksum(tv%frazil, mesg//" tv%frazil", G%HI, scale=US%Q_to_J_kg*US%RZ_to_kg_m2)
   if (associated(tv%salt_deficit)) &
-    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, scale=US%RZ_to_kg_m2)
+    call hchksum(tv%salt_deficit, mesg//" tv%salt_deficit", G%HI, scale=US%RZ_to_kg_m2*US%S_to_ppt)
   if (associated(tv%TempxPmE)) &
-    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=US%RZ_to_kg_m2)
+    call hchksum(tv%TempxPmE, mesg//" tv%TempxPmE", G%HI, scale=US%RZ_to_kg_m2*US%C_to_degC)
 end subroutine MOM_thermovar_chksum
 
 end module MOM_variables

src/diagnostics/MOM_sum_output.F90

@@ -683,7 +683,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
   if (CS%use_temperature) then
     Temp_int(:,:) = 0.0 ; Salt_int(:,:) = 0.0
     do k=1,nz ; do j=js,je ; do i=is,ie
-      Salt_int(i,j) = Salt_int(i,j) + tv%S(i,j,k) * &
+      Salt_int(i,j) = Salt_int(i,j) + US%S_to_ppt*tv%S(i,j,k) * &
                       (h(i,j,k)*(HL2_to_kg * areaTm(i,j)))
       Temp_int(i,j) = Temp_int(i,j) + (US%Q_to_J_kg*tv%C_p * tv%T(i,j,k)) * &
                       (h(i,j,k)*(HL2_to_kg * areaTm(i,j)))
@@ -756,9 +756,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
   mass_chg = EFP_to_real(mass_chg_EFP)
 
   if (CS%use_temperature) then
-    salin = Salt / mass_tot ; salin_anom = Salt_anom / mass_tot
+    salin = Salt / mass_tot
+    salin_anom = Salt_anom / mass_tot
    ! salin_chg = Salt_chg / mass_tot
-    temp = heat / (mass_tot*US%Q_to_J_kg*tv%C_p) ; temp_anom = Heat_anom / (mass_tot*US%Q_to_J_kg*tv%C_p)
+    temp = heat / (mass_tot*US%Q_to_J_kg*US%degC_to_C*tv%C_p)
+    temp_anom = Heat_anom / (mass_tot*US%Q_to_J_kg*US%degC_to_C*tv%C_p)
   endif
   En_mass = toten / mass_tot
 
@@ -995,18 +997,18 @@ subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
     ! smg: old code
     if (associated(tv%TempxPmE)) then
       do j=js,je ; do i=is,ie
-        heat_in(i,j) = heat_in(i,j) + (fluxes%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%TempxPmE(i,j)
+        heat_in(i,j) = heat_in(i,j) + (tv%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%TempxPmE(i,j)
       enddo ; enddo
     elseif (associated(fluxes%evap)) then
       do j=js,je ; do i=is,ie
-        heat_in(i,j) = heat_in(i,j) + (US%Q_to_J_kg*fluxes%C_p * sfc_state%SST(i,j)) * FW_in(i,j)
+        heat_in(i,j) = heat_in(i,j) + (US%Q_to_J_kg*tv%C_p * US%degC_to_C*sfc_state%SST(i,j)) * FW_in(i,j)
       enddo ; enddo
     endif
 
     ! The following heat sources may or may not be used.
     if (associated(tv%internal_heat)) then
       do j=js,je ; do i=is,ie
-        heat_in(i,j) = heat_in(i,j) + (fluxes%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%internal_heat(i,j)
+        heat_in(i,j) = heat_in(i,j) + (tv%C_p * QRZL2_to_J*G%areaT(i,j)) * tv%internal_heat(i,j)
       enddo ; enddo
     endif
     if (associated(tv%frazil)) then ; do j=js,je ; do i=is,ie

src/equation_of_state/MOM_EOS.F90

@@ -1536,10 +1536,10 @@ subroutine EOS_init(param_file, EOS, US)
   EOS%R_to_kg_m3 = 1. ; if (present(US)) EOS%R_to_kg_m3 = US%R_to_kg_m3
   EOS%RL2_T2_to_Pa = 1. ; if (present(US)) EOS%RL2_T2_to_Pa = US%RL2_T2_to_Pa
   EOS%L_T_to_m_s = 1. ; if (present(US)) EOS%L_T_to_m_s = US%L_T_to_m_s
-  EOS%degC_to_C = 1. !### ; if (present(US)) EOS%degC_to_C = US%degC_to_C
-  EOS%C_to_degC = 1. !### ; if (present(US)) EOS%C_to_degC = US%C_to_degC
-  EOS%ppt_to_S = 1. !### ; if (present(US)) EOS%ppt_to_S = US%ppt_to_S
-  EOS%S_to_ppt = 1. !### ; if (present(US)) EOS%S_to_ppt = US%S_to_ppt
+  EOS%degC_to_C = 1. ; if (present(US)) EOS%degC_to_C = US%degC_to_C
+  EOS%C_to_degC = 1. ; if (present(US)) EOS%C_to_degC = US%C_to_degC
+  EOS%ppt_to_S = 1. ; if (present(US)) EOS%ppt_to_S = US%ppt_to_S
+  EOS%S_to_ppt = 1. ; if (present(US)) EOS%S_to_ppt = US%S_to_ppt
 
 end subroutine EOS_init
 

src/ice_shelf/MOM_ice_shelf.F90

@@ -244,9 +244,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
   real, dimension(SZI_(CS%grid)) :: &
     Rhoml, &   !< Ocean mixed layer density [R ~> kg m-3].
     dR0_dT, &  !< Partial derivative of the mixed layer density
-               !< with temperature [R degC-1 ~> kg m-3 degC-1].
+               !< with temperature [R C-1 ~> kg m-3 degC-1].
     dR0_dS, &  !< Partial derivative of the mixed layer density
-               !< with salinity [R ppt-1 ~> kg m-3 ppt-1].
+               !< with salinity [R S-1 ~> kg m-3 ppt-1].
     p_int      !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa].
 
   real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
@@ -426,10 +426,10 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
     do i=is,ie ; p_int(i) = CS%g_Earth * ISS%mass_shelf(i,j) ; enddo
 
     ! Calculate insitu densities and expansion coefficients
-    call calculate_density(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, Rhoml(:), &
-                                 CS%eqn_of_state, EOSdom)
-    call calculate_density_derivs(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, dR0_dT, dR0_dS, &
+    call calculate_density(US%degC_to_C*sfc_state%sst(:,j), US%ppt_to_S*sfc_state%sss(:,j), p_int, Rhoml(:), &
                                  CS%eqn_of_state, EOSdom)
+    call calculate_density_derivs(US%degC_to_C*sfc_state%sst(:,j), US%ppt_to_S*sfc_state%sss(:,j), p_int, &
+                                  dR0_dT, dR0_dS, CS%eqn_of_state, EOSdom)
 
     do i=is,ie
       if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. &
@@ -451,8 +451,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step, CS)
           hBL_neut_h_molec = ZETA_N * ((hBL_neut * ustar_h) / (5.0 * CS%kv_molec))
 
           ! Determine the mixed layer buoyancy flux, wB_flux.
-          dB_dS = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dS(i)
-          dB_dT = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * dR0_dT(i)
+          dB_dS = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * US%ppt_to_S*dR0_dS(i)
+          dB_dT = (US%L_to_Z**2*CS%g_Earth / Rhoml(i)) * US%degC_to_C*dR0_dT(i)
           ln_neut = 0.0 ; if (hBL_neut_h_molec > 1.0) ln_neut = log(hBL_neut_h_molec)
 
           if (CS%find_salt_root) then

src/initialization/MOM_coord_initialization.F90

@@ -204,8 +204,8 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
                                                       !! [R L2 T-2 ~> Pa].
 
   ! Local variables
-  real :: T_ref   ! Reference temperature
-  real :: S_ref   ! Reference salinity
+  real :: T_ref   ! Reference temperature [C ~> degC]
+  real :: S_ref   ! Reference salinity [S ~> ppt]
   real :: g_int   ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2].
   real :: g_fs    ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
   character(len=40)  :: mdl = "set_coord_from_TS_ref" ! This subroutine's name.
@@ -214,11 +214,11 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
 
   call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90")
 
-  call get_param(param_file, mdl, "T_REF", T_Ref, &
-                 "The initial temperature of the lightest layer.", units="degC", &
-                 fail_if_missing=.true.)
-  call get_param(param_file, mdl, "S_REF", S_Ref, &
-                 "The initial salinities.", units="PSU", default=35.0)
+  call get_param(param_file, mdl, "T_REF", T_ref, &
+                 "The initial temperature of the lightest layer.", &
+                 units="degC", scale=US%degC_to_C, fail_if_missing=.true.)
+  call get_param(param_file, mdl, "S_REF", S_ref, &
+                 "The initial salinities.", units="PSU", default=35.0, scale=US%ppt_to_S)
   call get_param(param_file, mdl, "GFS", g_fs, &
                  "The reduced gravity at the free surface.", units="m s-2", &
                  default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m)
@@ -254,7 +254,9 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
                                                    !! [R L2 T-2 ~> Pa].
 
   ! Local variables
-  real, dimension(GV%ke) :: T0, S0,  Pref
+  real, dimension(GV%ke) :: T0   ! A profile of temperatures [C ~> degC]
+  real, dimension(GV%ke) :: S0   ! A profile of salinities [S ~> ppt]
+  real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa]
   real :: g_fs    ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
   integer :: k, nz
   character(len=40)  :: mdl = "set_coord_from_TS_profile" ! This subroutine's name.
@@ -274,8 +276,8 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s
   filename = trim(slasher(inputdir))//trim(coord_file)
   call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename)
 
-  call MOM_read_data(filename, "PTEMP", T0(:))
-  call MOM_read_data(filename, "SALT", S0(:))
+  call MOM_read_data(filename, "PTEMP", T0(:), scale=US%degC_to_C)
+  call MOM_read_data(filename, "SALT", S0(:), scale=US%ppt_to_S)
 
   if (.not.file_exists(filename)) call MOM_error(FATAL, &
       " set_coord_from_TS_profile: Unable to open " //trim(filename))
@@ -301,9 +303,13 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
                                                    !! [R L2 T-2 ~> Pa].
 
   ! Local variables
-  real, dimension(GV%ke) :: T0, S0,  Pref
-  real :: S_Ref, S_Light, S_Dense ! Salinity range parameters [ppt].
-  real :: T_Ref, T_Light, T_Dense ! Temperature range parameters [degC].
+  real, dimension(GV%ke) :: T0   ! A profile of temperatures [C ~> degC]
+  real, dimension(GV%ke) :: S0   ! A profile of salinities [S ~> ppt]
+  real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa]
+  real :: S_Ref   ! Default salinity range parameters [ppt].
+  real :: T_Ref   ! Default temperature range parameters [degC].
+  real :: S_Light, S_Dense ! Salinity range parameters [S ~> ppt].
+  real :: T_Light, T_Dense ! Temperature range parameters [C ~> degC].
   real :: res_rat ! The ratio of density space resolution in the denser part
                   ! of the range to that in the lighter part of the range.
                   ! Setting this greater than 1 increases the resolution for
@@ -321,19 +327,19 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
                  "The default initial temperatures.", units="degC", default=10.0)
   call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", T_Light, &
                  "The initial temperature of the lightest layer when "//&
-                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref)
+                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref, scale=US%degC_to_C)
   call get_param(param_file, mdl, "TS_RANGE_T_DENSE", T_Dense, &
                  "The initial temperature of the densest layer when "//&
-                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref)
+                 "COORD_CONFIG is set to ts_range.", units="degC", default=T_Ref, scale=US%degC_to_C)
 
   call get_param(param_file, mdl, "S_REF", S_Ref, &
                  "The default initial salinities.", units="PSU", default=35.0)
   call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_Light, &
                  "The initial lightest salinities when COORD_CONFIG "//&
-                 "is set to ts_range.", default = S_Ref, units="PSU")
+                 "is set to ts_range.", default = S_Ref, units="PSU", scale=US%ppt_to_S)
   call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_Dense, &
                  "The initial densest salinities when COORD_CONFIG "//&
-                 "is set to ts_range.", default = S_Ref, units="PSU")
+                 "is set to ts_range.", default = S_Ref, units="PSU", scale=US%ppt_to_S)
 
   call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, &
                  "The ratio of density space resolution in the densest "//&
@@ -352,7 +358,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta
   k_light = GV%nk_rho_varies + 1
 
   ! Set T0(k) to range from T_LIGHT to T_DENSE, and simliarly for S0(k).
-  T0(k_light) = T_light ; S0(k_light) = S_light
+  T0(k_light) = T_Light ; S0(k_light) = S_Light
   a1 = 2.0 * res_rat / (1.0 + res_rat)
   do k=k_light+1,nz
     k_frac = real(k-k_light)/real(nz-k_light)