+Rescale vars in MOM_temp_salt_initialize_from_Z

  Works with rescaled temperatures and salinities for the internal calculations
in MOM_temp_salt_initialize_from_Z, taking advantage of the recently corrected
rescaling capabilities in horiz_interp_and_extrap_tracer.  This commit also
includes these other closely related changes.

 - Modified convert_temp_salt_for_TEOS10 to work with rescaled temperature and
   salinity units

 - Eliminated unused land_fill argument from determine_temperature

 - Slightly refactored tracer_z_init_array to avoid needing an extra set of do
   loop through the 3-d array when a scale argument is present

All answers are bitwise identical, but there are changes in the arguments to
publicly visible interfaces.

src/equation_of_state/MOM_EOS.F90

@@ -1606,9 +1606,9 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
   integer,               intent(in)    :: kd  !< The number of layers to work on
   type(hor_index_type),  intent(in)    :: HI       !< The horizontal index structure
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
-                         intent(inout) :: T   !< Potential temperature referenced to the surface [degC]
+                         intent(inout) :: T   !< Potential temperature referenced to the surface [C ~> degC]
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
-                         intent(inout) :: S   !< Salinity [ppt]
+                         intent(inout) :: S   !< Salinity [S ~> ppt]
   real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
                          intent(in)    :: mask_z !< 3d mask regulating which points to convert.
   type(EOS_type),        intent(in)    :: EOS !< Equation of state structure
@@ -1620,12 +1620,12 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
 
   do k=1,kd ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec
     if (mask_z(i,j,k) >= 1.0) then
-      S(i,j,k) = gsw_sr_from_sp(S(i,j,k))
+      S(i,j,k) = EOS%ppt_to_S*gsw_sr_from_sp(EOS%S_to_ppt*S(i,j,k))
 !     Get absolute salinity from practical salinity, converting pressures from Pascal to dbar.
 !     If this option is activated, pressure will need to be added as an argument, and it should be
 !     moved out into module that is not shared between components, where the ocean_grid can be used.
 !     S(i,j,k) = gsw_sa_from_sp(S(i,j,k),pres(i,j,k)*1.0e-4,G%geoLonT(i,j),G%geoLatT(i,j))
-      T(i,j,k) = gsw_ct_from_pt(S(i,j,k), T(i,j,k))
+      T(i,j,k) = EOS%degC_to_C*gsw_ct_from_pt(EOS%S_to_ppt*S(i,j,k), EOS%S_to_ppt*T(i,j,k))
     endif
   enddo ; enddo ; enddo
 end subroutine convert_temp_salt_for_TEOS10

src/framework/MOM_diag_mediator.F90

@@ -306,8 +306,8 @@ module MOM_diag_mediator
 
   ! Pointer to H, G and T&S needed for remapping
   real, dimension(:,:,:), pointer :: h => null() !< The thicknesses needed for remapping [H ~> m or kg m-2]
-  real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping [degC]
-  real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [ppt]
+  real, dimension(:,:,:), pointer :: T => null() !< The temperatures needed for remapping [C ~> degC]
+  real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [S ~> ppt]
   type(EOS_type),  pointer :: eqn_of_state => null() !< The equation of state type
   type(ocean_grid_type), pointer :: G => null()  !< The ocean grid type
   type(verticalGrid_type), pointer :: GV => null()  !< The model's vertical ocean grid

src/initialization/MOM_state_initialization.F90

@@ -1783,9 +1783,9 @@ subroutine initialize_temp_salt_linear(T, S, G, GV, US, param_file, just_read)
   type(ocean_grid_type),                     intent(in)  :: G  !< The ocean's grid structure
   type(verticalGrid_type),                   intent(in)  :: GV !< The ocean's vertical grid structure
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T  !< The potential temperature that is
-                                                               !! being initialized [degC]
+                                                               !! being initialized [C ~> degC]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S  !< The salinity that is
-                                                               !! being initialized [ppt]
+                                                               !! being initialized [S ~> ppt]
   type(unit_scale_type),                     intent(in)  :: US !< A dimensional unit scaling type
   type(param_file_type),                     intent(in)  :: param_file !< A structure to parse for
                                                                !! run-time parameters
@@ -2419,23 +2419,24 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
   real    :: PI_180   ! for conversion from degrees to radians
   real    :: Hmix_default ! The default initial mixed layer depth [m].
   real    :: Hmix_depth   ! The mixed layer depth in the initial condition [Z ~> m].
-  real    :: missing_value_temp, missing_value_salt
+  real    :: missing_value_temp  ! The missing value in the input temperature field
+  real    :: missing_value_salt  ! The missing value in the input salinity field
   logical :: correct_thickness
   real    :: h_tolerance ! A parameter that controls the tolerance when adjusting the
                          ! thickness to fit the bathymetry [Z ~> m].
   character(len=40) :: potemp_var, salin_var
 
   integer, parameter :: niter=10   ! number of iterations for t/s adjustment to layer density
   logical            :: adjust_temperature = .true.  ! fit t/s to target densities
-  real, parameter    :: missing_value = -1.e20
-  real, parameter    :: temp_land_fill = 0.0, salt_land_fill = 35.0
+  real    :: temp_land_fill  ! A temperature value to use for land points [C ~> degC]
+  real    :: salt_land_fill  ! A salinity value to use for land points [C ~> degC]
   logical :: reentrant_x, tripolar_n
 
   ! data arrays
   real, dimension(:), allocatable :: z_edges_in, z_in ! Interface heights [Z ~> m]
   real, dimension(:), allocatable :: Rb  ! Interface densities [R ~> kg m-3]
-  real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [degC]
-  real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [ppt]
+  real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [C ~> degC]
+  real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt]
   real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim]
   real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi   ! Interface heights [Z ~> m].
@@ -2464,8 +2465,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
                                    ! from data when finding the initial interface locations in
                                    ! layered mode from a dataset of T and S.
   character(len=64) :: remappingScheme
-  real :: tempAvg  ! Spatially averaged temperatures on a layer [degC]
-  real :: saltAvg  ! Spatially averaged salinities on a layer [ppt]
+  real :: tempAvg  ! Spatially averaged temperatures on a layer [C ~> degC]
+  real :: saltAvg  ! Spatially averaged salinities on a layer [S ~> ppt]
   logical :: do_conv_adj, ignore
   integer :: nPoints
   integer :: id_clock_routine, id_clock_ALE
@@ -2592,6 +2593,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
     return ! All run-time parameters have been read, so return.
   endif
 
+  !### These hard-coded constants should be made into runtime parameters
+  temp_land_fill = 0.0*US%degC_to_C
+  salt_land_fill = 35.0*US%ppt_to_S
+
   eps_z = GV%Angstrom_Z
   eps_rho = 1.0e-10*US%kg_m3_to_R
 
@@ -2610,35 +2615,24 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
   ! to the North/South Pole past the limits of the input data, they are extrapolated using the average
   ! value at the northernmost/southernmost latitude.
 
-  call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, &
+  call horiz_interp_and_extrap_tracer(tfilename, potemp_var, US%degC_to_C, 1, &
        G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, &
-       tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded)
+       tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, &
+       ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%degC_to_C)
 
-  call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, &
+  call horiz_interp_and_extrap_tracer(sfilename, salin_var, US%ppt_to_S, 1, &
        G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, &
-       tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded)
+       tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, &
+       ongrid=pre_gridded, tr_iter_tol=1.0e-3*US%ppt_to_S)
 
   kd = size(z_in,1)
 
   ! Convert the sign convention of Z_edges_in.
   do k=1,size(Z_edges_in,1) ; Z_edges_in(k) = -Z_edges_in(k) ; enddo
 
-  allocate(rho_z(isd:ied,jsd:jed,kd))
-
   ! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
   call convert_temp_salt_for_TEOS10(temp_z, salt_z, G%HI, kd, mask_z, eos)
 
-  press(:) = tv%P_Ref
-  EOSdom(:) = EOS_domain(G%HI)
-  do k=1,kd ; do j=js,je
-    call calculate_density(US%degC_to_C*temp_z(:,j,k), US%ppt_to_S*salt_z(:,j,k), press, rho_z(:,j,k), eos, EOSdom)
-  enddo ; enddo
-
-  call pass_var(temp_z,G%Domain)
-  call pass_var(salt_z,G%Domain)
-  call pass_var(mask_z,G%Domain)
-  call pass_var(rho_z,G%Domain)
-
   do j=js,je ; do i=is,ie
     Z_bottom(i,j) = -depth_tot(i,j)
   enddo ; enddo
@@ -2661,15 +2655,15 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
         do k = 1, nkd
           if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then
             zBottomOfCell = max( z_edges_in(k+1), Z_bottom(i,j))
-            tmpT1dIn(i,j,k) = US%degC_to_C*temp_z(i,j,k)
-            tmpS1dIn(i,j,k) = US%ppt_to_S*salt_z(i,j,k)
+            tmpT1dIn(i,j,k) = temp_z(i,j,k)
+            tmpS1dIn(i,j,k) = salt_z(i,j,k)
           elseif (k>1) then
             zBottomOfCell = Z_bottom(i,j)
             tmpT1dIn(i,j,k) = tmpT1dIn(i,j,k-1)
             tmpS1dIn(i,j,k) = tmpS1dIn(i,j,k-1)
           else ! This next block should only ever be reached over land
-            tmpT1dIn(i,j,k) = -99.9*US%degC_to_C
-            tmpS1dIn(i,j,k) = -99.9*US%ppt_to_S
+            tmpT1dIn(i,j,k) = -99.9*US%degC_to_C ! Change to temp_land_fill
+            tmpS1dIn(i,j,k) = -99.9*US%ppt_to_S  ! Change to salt_land_fill
           endif
           h1(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell)
           zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k
@@ -2679,9 +2673,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
       endif ! mask2dT
     enddo ; enddo
     deallocate( tmp_mask_in )
-    call pass_var(h1, G%Domain)
-    call pass_var(tmpT1dIn, G%Domain)
-    call pass_var(tmpS1dIn, G%Domain)
 
     ! Build the target grid (and set the model thickness to it)
     ! This call can be more general but is hard-coded for z* coordinates...  ????
@@ -2705,7 +2696,6 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
           h(i,j,:) = 0.
         endif ! mask2dT
       enddo ; enddo
-      call pass_var(h, G%Domain)
       deallocate( hTarget )
     endif
 
@@ -2756,9 +2746,18 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
 
     nkml = 0 ; if (separate_mixed_layer) nkml = GV%nkml
 
+    press(:) = tv%P_Ref
+    EOSdom(:) = EOS_domain(G%HI)
+    allocate(rho_z(isd:ied,jsd:jed,kd))
+    do k=1,kd ; do j=js,je
+      call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), eos, EOSdom)
+    enddo ; enddo
+
     call find_interfaces(rho_z, z_in, kd, Rb, Z_bottom, zi, G, GV, US, nlevs, nkml, &
                          Hmix_depth, eps_z, eps_rho, density_extrap_bug)
 
+    deallocate(rho_z)
+
     if (correct_thickness) then
       call adjustEtaToFitBathymetry(G, GV, US, zi, h, h_tolerance, dZ_ref_eta=G%Z_ref)
     else
@@ -2784,21 +2783,21 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
       endif
     endif
 
-    call tracer_z_init_array(temp_z, z_edges_in, kd, zi, missing_value, G, nz, nlevs, eps_z, &
-                             tv%T, scale=US%degC_to_C)
-    call tracer_z_init_array(salt_z, z_edges_in, kd, zi, missing_value, G, nz, nlevs, eps_z, &
-                             tv%S, scale=US%ppt_to_S)
-
-    do k=1,nz
-      nPoints = 0 ; tempAvg = 0. ; saltAvg = 0.
-      do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
-        nPoints = nPoints + 1
-        tempAvg = tempAvg + US%C_to_degC*tv%T(i,j,k)
-        saltAvg = saltAvg + US%S_to_ppt*tv%S(i,j,k)
-      endif ; enddo ; enddo
+    call tracer_z_init_array(temp_z, z_edges_in, kd, zi, temp_land_fill, G, nz, nlevs, eps_z, &
+                             tv%T)
+    call tracer_z_init_array(salt_z, z_edges_in, kd, zi, salt_land_fill, G, nz, nlevs, eps_z, &
+                             tv%S)
 
+    if (homogenize) then
       ! Horizontally homogenize data to produce perfectly "flat" initial conditions
-      if (homogenize) then
+      do k=1,nz
+        nPoints = 0 ; tempAvg = 0. ; saltAvg = 0.
+        do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
+          nPoints = nPoints + 1
+          tempAvg = tempAvg + tv%T(i,j,k)
+          saltAvg = saltAvg + tv%S(i,j,k)
+        endif ; enddo ; enddo
+
         !### These averages will not reproduce across PE layouts or grid rotation.
         call sum_across_PEs(nPoints)
         call sum_across_PEs(tempAvg)
@@ -2807,32 +2806,21 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
           tempAvg = tempAvg / real(nPoints)
           saltAvg = saltAvg / real(nPoints)
         endif
-        tv%T(:,:,k) = US%degC_to_C*tempAvg
-        tv%S(:,:,k) = US%ppt_to_S*saltAvg
-      endif
-    enddo
-
-  endif ! useALEremapping
-
-  ! Fill land values
-  do k=1,nz ; do j=js,je ; do i=is,ie
-    if (tv%T(i,j,k) == US%degC_to_C*missing_value) then
-      tv%T(i,j,k) = US%degC_to_C*temp_land_fill
-      tv%S(i,j,k) = US%ppt_to_S*salt_land_fill
+        tv%T(:,:,k) = tempAvg
+        tv%S(:,:,k) = saltAvg
+      enddo
     endif
-  enddo ; enddo ; enddo
 
+    if (adjust_temperature) then
+      ! Finally adjust to target density
+      ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1
+      call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), tv%P_Ref, niter, &
+                                 h, ks, G, GV, US, eos)
+    endif
 
-  if (adjust_temperature .and. .not. useALEremapping) then
-    ! Finally adjust to target density
-    ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1
-    call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), tv%P_Ref, niter, &
-                               missing_value, h, ks, G, GV, US, eos)
-  endif
+  endif ! useALEremapping
 
   deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
-  deallocate(rho_z)
-
 
   call pass_var(h, G%Domain)
   call pass_var(tv%T, G%Domain)
@@ -2984,7 +2972,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv)
   ! Local variables
   integer, parameter :: nk=5
   real, dimension(nk) :: T, T_t, T_b ! Temperatures [C ~> degC]
-  real, dimension(nk) :: S, S_t, S_b ! Salinities [ppt]
+  real, dimension(nk) :: S, S_t, S_b ! Salinities [S ~> ppt]
   real, dimension(nk) :: rho ! Layer density [R ~> kg m-3]
   real, dimension(nk) :: h   ! Layer thicknesses [H ~> m or kg m-2]
   real, dimension(nk) :: z   ! Height of layer center [Z ~> m]

src/tracer/MOM_tracer_Z_init.F90

@@ -42,7 +42,7 @@ function tracer_Z_init(tr, h, filename, tr_name, G, GV, US, missing_val, land_va
   real,        optional, intent(in)    :: land_val !< A value to use to fill in land points
 
 ! This include declares and sets the variable "version".
-#include "version_variable.h"
+# include "version_variable.h"
 
   real, allocatable, dimension(:,:,:) :: &
     tr_in   ! The z-space array of tracer concentrations that is read in.
@@ -283,7 +283,7 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
   integer,                    intent(in)  :: nlay  !< The number of vertical layers in the target grid
   real, dimension(SZI_(G),SZJ_(G),nlay+1), &
                               intent(in)  :: e     !< The depths of the target layer interfaces [Z ~> m] or [m]
-  real,                       intent(in)  :: land_fill !< fill in data over land [A]
+  real,                       intent(in)  :: land_fill !< fill in data over land [B]
   integer, dimension(SZI_(G),SZJ_(G)), &
                               intent(in)  :: nlevs !< The number of input levels with valid data
   real,                       intent(in)  :: eps_z !< A negligibly thin layer thickness [Z ~> m].
@@ -293,19 +293,22 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
                                                    !! input tracers [B A-1 ~> 1]
 
   ! Local variables
-  real :: tr_1d(nk_data) ! A copy of the input tracer concentrations in a column [A]
+  real :: tr_1d(nk_data) ! A copy of the input tracer concentrations in a column [B]
   real :: e_1d(nlay+1)   ! A 1-d column of interface heights, in the same units as e [Z ~> m] or [m]
-  real :: sl_tr          ! The tracer concentration slope times the layer thickness, in tracer units [A]
+  real :: sl_tr          ! The tracer concentration slope times the layer thickness, in tracer units [B]
   real :: wt(nk_data)    ! The fractional weight for each layer in the range between z1 and z2 [nondim]
   real :: z1(nk_data)    ! z1 and z2 are the fractional depths of the top and bottom
   real :: z2(nk_data)    ! limits of the part of a z-cell that contributes to a layer, relative
                          ! to the cell center and normalized by the cell thickness [nondim].
                          ! Note that -1/2 <= z1 <= z2 <= 1/2.
+  real :: scale_fac      ! A factor by which to scale the output tracers from the input tracers [B A-1 ~> 1]
   integer :: k_top, k_bot, k_bot_prev, kstart
   integer :: i, j, k, kz, is, ie, js, je
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
+  scale_fac = 1.0 ; if (present(scale)) then ; scale_fac = scale ; endif
+
   do j=js,je
     i_loop: do i=is,ie
       if (nlevs(i,j) == 0 .or. G%mask2dT(i,j) == 0.) then
@@ -314,7 +317,7 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
       endif
 
       do k=1,nk_data
-        tr_1d(k) = tr_in(i,j,k)
+        tr_1d(k) = scale_fac*tr_in(i,j,k)
       enddo
 
       do k=1,nlay+1
@@ -372,12 +375,6 @@ subroutine tracer_z_init_array(tr_in, z_edges, nk_data, e, land_fill, G, nlay, n
     enddo i_loop
   enddo
 
-  if (present(scale)) then ; if (scale /= 1.0) then
-    do k=1,nlay ; do j=js,je ; do i=is,ie
-      tr(i,j,k) = scale*tr(i,j,k)
-    enddo ; enddo ; enddo
-  endif ; endif
-
 end subroutine tracer_z_init_array
 
 !> This subroutine reads the vertical coordinate data for a field from a NetCDF file.
@@ -559,7 +556,7 @@ end function find_limited_slope
 
 !> This subroutine determines the potential temperature and salinity that
 !! is consistent with the target density using provided initial guess
-subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, k_start, G, GV, US, &
+subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, h, k_start, G, GV, US, &
                                  EOS, h_massless)
   type(ocean_grid_type),         intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type),       intent(in)    :: GV   !< The ocean's vertical grid structure.
@@ -571,7 +568,6 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,
   real,                          intent(in)    :: p_ref !< reference pressure [R L2 T-2 ~> Pa].
   integer,                       intent(in)    :: niter !< maximum number of iterations
   integer,                       intent(in)    :: k_start !< starting index (i.e. below the buffer layer)
-  real,                          intent(in)    :: land_fill !< land fill value
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                  intent(in)    :: h   !< layer thickness, used only to avoid working on
                                                       !! massless layers [H ~> m or kg m-2]