+Add rescaling for temperature and salinity (1)

  Added dimensional rescaling for temperature and salinity, as determined by the
new runtime parameters C_RESCALE_POWER and S_RESCALE_POWER.  With this change
there are 4 new elements in the transparent unit_scale_type, and these are
widely used in the code.  In addition, ### other files were added that had
checksum calls or diagnostics rescaled by these new factors, and where comments
were changed, but were otherwise unaltered as a result of the new dimensional
rescaling.  There will be another commit very shortly that completes the changes
and leads to fully functional dimensional rescaling for temperatures and
salinities, but these will involve more extensive code or interface changes, but
this commit will be useful for any possible git-bisection of any potential
changes that do not involve dimensional rescaling.  All solutions in existing
test cases are bitwise identical.

src/ALE/MOM_ALE.F90

@@ -318,9 +318,9 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
       'Layer Thickness before remapping', get_thickness_units(GV), conversion=GV%H_to_MKS, &
       v_extensive=.true.)
   CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, &
-      'Temperature before remapping', 'degC')
+      'Temperature before remapping', 'degC', conversion=US%C_to_degC)
   CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, &
-      'Salinity before remapping', 'PSU')
+      'Salinity before remapping', 'PSU', conversion=US%S_to_ppt)
   CS%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, Time, &
       'Interface Heights before remapping', 'm', conversion=US%Z_to_m)
 
@@ -451,8 +451,8 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
 
   if (CS%debug) then
     call hchksum(h,    "Post-ALE_main h", G%HI, haloshift=0, scale=GV%H_to_m)
-    call hchksum(tv%T, "Post-ALE_main T", G%HI, haloshift=0)
-    call hchksum(tv%S, "Post-ALE_main S", G%HI, haloshift=0)
+    call hchksum(tv%T, "Post-ALE_main T", G%HI, haloshift=0, scale=US%C_to_degC)
+    call hchksum(tv%S, "Post-ALE_main S", G%HI, haloshift=0, scale=US%S_to_ppt)
     call uvchksum("Post-ALE_main [uv]", u, v, G%HI, haloshift=0, scale=US%L_T_to_m_s)
   endif
 
@@ -1160,13 +1160,13 @@ subroutine TS_PLM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap
   type(verticalGrid_type), intent(in)    :: GV   !< Ocean vertical grid structure
   type(ALE_CS),            intent(inout) :: CS   !< module control structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                           intent(inout) :: S_t  !< Salinity at the top edge of each layer
+                           intent(inout) :: S_t  !< Salinity at the top edge of each layer [S ~> ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                           intent(inout) :: S_b  !< Salinity at the bottom edge of each layer
+                           intent(inout) :: S_b  !< Salinity at the bottom edge of each layer [S ~> ppt]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                           intent(inout) :: T_t  !< Temperature at the top edge of each layer
+                           intent(inout) :: T_t  !< Temperature at the top edge of each layer [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
-                           intent(inout) :: T_b  !< Temperature at the bottom edge of each layer
+                           intent(inout) :: T_b  !< Temperature at the bottom edge of each layer [C ~> degC]
   type(thermo_var_ptrs),   intent(in)    :: tv   !< thermodynamics structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                            intent(in)    :: h    !< layer thickness [H ~> m or kg m-2]

src/core/MOM_PressureForce_FV.F90

@@ -98,14 +98,14 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
     T_tmp, &    ! Temporary array of temperatures where layers that are lighter
-                ! than the mixed layer have the mixed layer's properties [degC].
+                ! than the mixed layer have the mixed layer's properties [C ~> degC].
     S_tmp       ! Temporary array of salinities where layers that are lighter
-                ! than the mixed layer have the mixed layer's properties [ppt].
+                ! than the mixed layer have the mixed layer's properties [S ~> ppt].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     S_t, &      ! Top and bottom edge values for linear reconstructions
-    S_b, &      ! of salinity within each layer [ppt].
+    S_b, &      ! of salinity within each layer [S ~> ppt].
     T_t, &      ! Top and bottom edge values for linear reconstructions
-    T_b         ! of temperature within each layer [degC].
+    T_b         ! of temperature within each layer [C ~> degC].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
     dza, &      ! The change in geopotential anomaly between the top and bottom
                 ! of a layer [L2 T-2 ~> m2 s-2].
@@ -467,12 +467,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
 
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
     T_tmp, &    ! Temporary array of temperatures where layers that are lighter
-                ! than the mixed layer have the mixed layer's properties [degC].
+                ! than the mixed layer have the mixed layer's properties [C ~> degC].
     S_tmp       ! Temporary array of salinities where layers that are lighter
-                ! than the mixed layer have the mixed layer's properties [ppt].
+                ! than the mixed layer have the mixed layer's properties [S ~> ppt].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
-    S_t, S_b, T_t, T_b ! Top and bottom edge values for linear reconstructions
-                       ! of salinity and temperature within each layer.
+    S_t, &      ! Top and bottom edge values for linear reconstructions
+    S_b, &      ! of salinity within each layer [S ~> ppt].
+    T_t, &      ! Top and bottom edge values for linear reconstructions
+    T_b         ! of temperature within each layer [C ~> degC].
   real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
   real :: p_ref(SZI_(G))     !   The pressure used to calculate the coordinate
                              ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
@@ -487,9 +489,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   logical :: use_ALE         ! If true, use an ALE pressure reconstruction.
   logical :: use_EOS         ! If true, density is calculated from T & S using an equation of state.
   type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
-  real :: Tl(5)              ! copy and T in local stencil [degC]
-  real :: mn_T               ! mean of T in local stencil [degC]
-  real :: mn_T2              ! mean of T**2 in local stencil [degC2]
+  real :: Tl(5)              ! copy and T in local stencil [C ~> degC]
+  real :: mn_T               ! mean of T in local stencil [C ~> degC]
+  real :: mn_T2              ! mean of T**2 in local stencil [C2 ~> degC2]
   real :: hl(5)              ! Copy of local stencil of H [H ~> m]
   real :: r_sm_H             ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
   real, parameter :: C1_6 = 1.0/6.0
@@ -699,7 +701,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
       ! where the layers are located.
       if ( use_ALE .and. CS%Recon_Scheme > 0 ) then
         if ( CS%Recon_Scheme == 1 ) then
-          call int_density_dz_generic_plm(k, tv,  T_t, T_b, S_t, S_b, e, &
+          call int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, &
                     rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
                     G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
                     useMassWghtInterp=CS%useMassWghtInterp, &
@@ -866,7 +868,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
                  "Negative values disable the scheme.", units="nondim", default=-1.0)
   if (CS%Stanley_T2_det_coeff>=0.) then
     CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs_pgf', diag%axesTL, &
-        Time, 'SGS temperature variance used in PGF', 'degC2')
+        Time, 'SGS temperature variance used in PGF', 'degC2', conversion=US%C_to_degC**2)
   endif
   if (CS%tides) then
     CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &

src/core/MOM_checksum_packages.F90

@@ -124,12 +124,12 @@ subroutine MOM_thermo_chksum(mesg, tv, G, US, haloshift)
   integer :: hs
   hs=1 ; if (present(haloshift)) hs=haloshift
 
-  if (associated(tv%T)) call hchksum(tv%T, mesg//" T", G%HI, haloshift=hs)
-  if (associated(tv%S)) call hchksum(tv%S, mesg//" S", G%HI, haloshift=hs)
+  if (associated(tv%T)) call hchksum(tv%T, mesg//" T", G%HI, haloshift=hs, scale=US%C_to_degC)
+  if (associated(tv%S)) call hchksum(tv%S, mesg//" S", G%HI, haloshift=hs, scale=US%S_to_ppt)
   if (associated(tv%frazil)) call hchksum(tv%frazil, mesg//" frazil", G%HI, haloshift=hs, &
                                           scale=US%Q_to_J_kg*US%R_to_kg_m3*US%Z_to_m)
-  if (associated(tv%salt_deficit)) &
-    call hchksum(tv%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs, scale=US%RZ_to_kg_m2)
+  if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, mesg//" salt deficit", G%HI, haloshift=hs, &
+                                                scale=US%S_to_ppt*US%RZ_to_kg_m2)
 
 end subroutine MOM_thermo_chksum
 
@@ -240,9 +240,9 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),  &
                            intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2].
   real, pointer, dimension(:,:,:),           &
-                           intent(in) :: Temp !< Temperature [degC].
+                           intent(in) :: Temp !< Temperature [C ~> degC].
   real, pointer, dimension(:,:,:),           &
-                           intent(in) :: Salt !< Salinity [ppt].
+                           intent(in) :: Salt !< Salinity [S ~> ppt].
   type(unit_scale_type),   intent(in) :: US    !< A dimensional unit scaling type
   logical,       optional, intent(in) :: allowChange !< do not flag an error
                                                      !! if the statistics change.
@@ -258,8 +258,11 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
   real :: Vol, dV    ! The total ocean volume and its change [m3] (unscaled to permit reproducing sum).
   real :: Area       ! The total ocean surface area [m2] (unscaled to permit reproducing sum).
   real :: h_minimum  ! The minimum layer thicknesses [H ~> m or kg m-2]
+  real :: T_scale    ! The scaling conversion factor for temperatures [degC C-1 ~> 1]
+  real :: S_scale    ! The scaling conversion factor for salinities [ppt S-1 ~> 1]
   logical :: do_TS   ! If true, evaluate statistics for temperature and salinity
-  type(stats) :: T, S, delT, delS
+  type(stats) :: T, delT ! Temperature statistics in unscaled units [degC]
+  type(stats) :: S, delS ! Salinity statistics in unscaled units [ppt]
 
   ! NOTE: save data is not normally allowed but we use it for debugging purposes here on the
   !       assumption we will not turn this on with threads
@@ -278,6 +281,8 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
   tmp_T(:,:) = 0.0
   tmp_S(:,:) = 0.0
 
+  T_scale = US%C_to_degC ; S_scale = US%S_to_ppt
+
   ! First collect local stats
   do j=js,je ; do i=is,ie
     tmp_A(i,j) = tmp_A(i,j) + US%L_to_m**2*G%areaT(i,j)
@@ -290,12 +295,12 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
       dV = US%L_to_m**2*G%areaT(i,j)*GV%H_to_m*h(i,j,k)
       tmp_V(i,j) = tmp_V(i,j) + dV
       if (do_TS .and. h(i,j,k)>0.) then
-        T%minimum = min( T%minimum, Temp(i,j,k) ) ; T%maximum = max( T%maximum, Temp(i,j,k) )
-        T%average = T%average + dV*Temp(i,j,k)
-        S%minimum = min( S%minimum, Salt(i,j,k) ) ; S%maximum = max( S%maximum, Salt(i,j,k) )
-        S%average = S%average + dV*Salt(i,j,k)
-        tmp_T(i,j) = tmp_T(i,j) + dV*Temp(i,j,k)
-        tmp_S(i,j) = tmp_S(i,j) + dV*Salt(i,j,k)
+        T%minimum = min( T%minimum, T_scale*Temp(i,j,k) ) ; T%maximum = max( T%maximum, T_scale*Temp(i,j,k) )
+        T%average = T%average + dV*T_scale*Temp(i,j,k)
+        S%minimum = min( S%minimum, S_scale*Salt(i,j,k) ) ; S%maximum = max( S%maximum, S_scale*Salt(i,j,k) )
+        S%average = S%average + dV*S_scale*Salt(i,j,k)
+        tmp_T(i,j) = tmp_T(i,j) + dV*T_scale*Temp(i,j,k)
+        tmp_S(i,j) = tmp_S(i,j) + dV*S_scale*Salt(i,j,k)
       endif
       if (h_minimum > h(i,j,k)) h_minimum = h(i,j,k)
     endif
@@ -343,11 +348,11 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
 
   if (do_TS .and. T%minimum<-5.0) then
     do j=js,je ; do i=is,ie
-      if (minval(Temp(i,j,:)) == T%minimum) then
+      if (minval(T_scale*Temp(i,j,:)) == T%minimum) then
         write(0,'(a,2f12.5)') 'x,y=', G%geoLonT(i,j), G%geoLatT(i,j)
         write(0,'(a3,3a12)') 'k','h','Temp','Salt'
         do k = 1, nz
-          write(0,'(i3,3es12.4)') k, h(i,j,k), Temp(i,j,k), Salt(i,j,k)
+          write(0,'(i3,3es12.4)') k, h(i,j,k), T_scale*Temp(i,j,k), S_scale*Salt(i,j,k)
         enddo
         stop 'Extremum detected'
       endif
@@ -360,7 +365,7 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
         write(0,'(a,2f12.5)') 'x,y=',G%geoLonT(i,j),G%geoLatT(i,j)
         write(0,'(a3,3a12)') 'k','h','Temp','Salt'
         do k = 1, nz
-          write(0,'(i3,3es12.4)') k, h(i,j,k), Temp(i,j,k), Salt(i,j,k)
+          write(0,'(i3,3es12.4)') k, h(i,j,k), T_scale*Temp(i,j,k), S_scale*Salt(i,j,k)
         enddo
         stop 'Negative thickness detected'
       endif

src/diagnostics/MOM_PointAccel.F90

@@ -50,8 +50,8 @@ module MOM_PointAccel
     v_av => NULL(), &       !< Time average velocity [L T-1 ~> m s-1]
     u_prev => NULL(), &     !< Previous u-velocity [L T-1 ~> m s-1]
     v_prev => NULL(), &     !< Previous v-velocity [L T-1 ~> m s-1]
-    T => NULL(), &          !< Temperature [degC]
-    S => NULL(), &          !< Salinity [ppt]
+    T => NULL(), &          !< Temperature [C ~> degC]
+    S => NULL(), &          !< Salinity [S ~> ppt]
     u_accel_bt => NULL(), & !< Barotropic u-accelerations [L T-2 ~> m s-2]
     v_accel_bt => NULL()    !< Barotropic v-accelerations [L T-2 ~> m s-2]
 end type PointAccel_CS
@@ -94,6 +94,8 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
   real    :: h_scale          ! A scaling factor for thicknesses [m H-1 ~> 1 or m3 kg-1]
   real    :: vel_scale        ! A scaling factor for velocities [m T s-1 L-1 ~> 1]
   real    :: uh_scale         ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1 or m3 kg-1]
+  real    :: temp_scale       ! A scaling factor for temperatures [degC C-1 ~> 1]
+  real    :: saln_scale       ! A scaling factor for salinities [ppt S-1 ~> 1]
   integer :: yr, mo, day, hr, minute, sec, yearday
   integer :: k, ks, ke
   integer :: nz
@@ -103,6 +105,7 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
 
   Angstrom = GV%Angstrom_H + GV%H_subroundoff
   h_scale = GV%H_to_m ; vel_scale = US%L_T_to_m_s ; uh_scale = GV%H_to_m*US%L_T_to_m_s
+  temp_scale = US%C_to_degC ; saln_scale = US%S_to_ppt
 
 !  if (.not.associated(CS)) return
   nz = GV%ke
@@ -257,15 +260,15 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
     do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo
     if (associated(CS%T)) then
       write(file,'(/,"T-:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') temp_scale*CS%T(i,j,k) ; enddo
       write(file,'(/,"T+:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i+1,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') temp_scale*CS%T(i+1,j,k) ; enddo
     endif
     if (associated(CS%S)) then
       write(file,'(/,"S-:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') saln_scale*CS%S(i,j,k) ; enddo
       write(file,'(/,"S+:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i+1,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') saln_scale*CS%S(i+1,j,k) ; enddo
     endif
 
     if (prev_avail) then
@@ -426,6 +429,8 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
   real    :: h_scale          ! A scaling factor for thicknesses [m H-1 ~> 1 or m3 kg-1]
   real    :: vel_scale        ! A scaling factor for velocities [m T s-1 L-1 ~> 1]
   real    :: uh_scale         ! A scaling factor for transport per unit length [m2 T s-1 L-1 H-1 ~> 1 or m3 kg-1]
+  real    :: temp_scale       ! A scaling factor for temperatures [degC C-1 ~> 1]
+  real    :: saln_scale       ! A scaling factor for salinities [ppt S-1 ~> 1]
   integer :: yr, mo, day, hr, minute, sec, yearday
   integer :: k, ks, ke
   integer :: nz
@@ -435,6 +440,7 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
 
   Angstrom = GV%Angstrom_H + GV%H_subroundoff
   h_scale = GV%H_to_m ; vel_scale = US%L_T_to_m_s ; uh_scale = GV%H_to_m*US%L_T_to_m_s
+  temp_scale = US%C_to_degC ; saln_scale = US%S_to_ppt
 
 !  if (.not.associated(CS)) return
   nz = GV%ke
@@ -592,15 +598,15 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
     do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo
     if (associated(CS%T)) then
       write(file,'(/,"T-:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') temp_scale*CS%T(i,j,k) ; enddo
       write(file,'(/,"T+:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%T(i,j+1,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') temp_scale*CS%T(i,j+1,k) ; enddo
     endif
     if (associated(CS%S)) then
       write(file,'(/,"S-:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') saln_scale*CS%S(i,j,k) ; enddo
       write(file,'(/,"S+:    ")', advance='no')
-      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') CS%S(i,j+1,k) ; enddo
+      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') saln_scale*CS%S(i,j+1,k) ; enddo
     endif
 
     if (prev_avail) then