Switched versions of calculate_density in 22 calls

  Changed to the new interfaces for calculate_density and related calls in 22
places.  All answers are bitwise identical.

src/ALE/MOM_regridding.F90

@@ -1890,7 +1890,7 @@ subroutine convective_adjustment(G, GV, h, tv)
   do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
 
     ! Compute densities within current water column
-    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state )
+    call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state, US=G%US)
 
     ! Repeat restratification until complete
     do
@@ -1909,9 +1909,9 @@ subroutine convective_adjustment(G, GV, h, tv)
           tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0
           h(i,j,k)    = h1 ; h(i,j,k+1)    = h0
           ! Recompute densities at levels k and k+1
-          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state)
+          call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state, US=G%US)
           call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
-                                  densities(k+1), tv%eqn_of_state )
+                                  densities(k+1), tv%eqn_of_state, US=G%US )
           stratified = .false.
         endif
       enddo  ! k

src/ALE/coord_rho.F90

@@ -250,7 +250,7 @@ subroutine build_rho_column_iteratively(CS, US, remapCS, nz, depth, h, T, S, eqn
     enddo
 
     ! Compute densities within current water column
-    call calculate_density( T_tmp, S_tmp, pres, densities, 1, nz, eqn_of_state, US=US)
+    call calculate_density( T_tmp, S_tmp, pres, densities, eqn_of_state, US=US)
 
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))

src/initialization/MOM_state_initialization.F90

@@ -1601,7 +1601,7 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
   enddo
 
   call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, US=US)
-  call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, US=US)
+  call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US=US, dom=(/1,1/))
 
   if (fit_salin) then
     ! A first guess of the layers' temperatures.

src/tracer/MOM_tracer_Z_init.F90

@@ -801,9 +801,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,
     adjust_salt = .true.
     iter_loop: do itt = 1,niter
       do k=1, nz
-        call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US)
-        call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, &
-                                      eos, US=US)
+        call calculate_density(T(:,k), S(:,k), press, rho(:,k), eos, US=US, dom=(/1,nx/))
+        call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), &
+                                      eos, US=US, dom=(/1,nx/))
       enddo
       do k=k_start,nz ; do i=1,nx
 
@@ -831,9 +831,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h,
 
     if (adjust_salt .and. old_fit) then ; do itt = 1,niter
       do k=1, nz
-        call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US)
-        call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, &
-                                      eos, US=US)
+        call calculate_density(T(:,k), S(:,k), press, rho(:,k), eos, US=US, dom=(/1,nx/))
+        call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), &
+                                      eos, US=US, dom=(/1,nx/))
       enddo
       do k=k_start,nz ; do i=1,nx
 !       if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then

src/user/ISOMIP_initialization.F90

@@ -184,10 +184,10 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read
     if (just_read) return ! All run-time parameters have been read, so return.
 
     ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
-    call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=US%kg_m3_to_R)
+    call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, US)
     ! write(mesg,*) 'Surface density is:', rho_sur
     ! call MOM_mesg(mesg,5)
-    call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=US%kg_m3_to_R)
+    call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, US)
     ! write(mesg,*) 'Bottom density is:', rho_bot
     ! call MOM_mesg(mesg,5)
     rho_range = rho_bot - rho_sur
@@ -301,10 +301,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
   call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, &
                  "Salinity at the bottom (interface)", units="ppt", default=34.55, do_not_log=just_read)
 
-  call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state, scale=US%kg_m3_to_R)
+  call calculate_density(t_sur, s_sur, 0.0, rho_sur, eqn_of_state, US)
   ! write(mesg,*) 'Density in the surface layer:', rho_sur
   ! call MOM_mesg(mesg,5)
-  call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state, scale=US%kg_m3_to_R)
+  call calculate_density(t_bot, s_bot, 0.0, rho_bot, eqn_of_state, US)
   ! write(mesg,*) 'Density in the bottom layer::', rho_bot
   ! call MOM_mesg(mesg,5)
 
@@ -362,7 +362,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
           ! call MOM_mesg(mesg,5)
         enddo
 
-        call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, US)
+        call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US, dom=(/1,1/))
         ! write(mesg,*) 'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1)
         ! call MOM_mesg(mesg,5)
         call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, US=US)
@@ -374,8 +374,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
           enddo
           ! Refine the guesses for each layer.
           do itt=1,6
-            call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
-            call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US)
+            call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US)
+            call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US)
             do k=1,nz
               S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS1)
             enddo
@@ -388,8 +388,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi
           enddo
 
           do itt=1,6
-            call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
-            call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US)
+            call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US)
+            call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US)
             do k=1,nz
               T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
             enddo
@@ -521,10 +521,10 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
   enddo ; enddo
 
   ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
-  call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=US%kg_m3_to_R)
+  call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, US)
   !write (mesg,*) 'Surface density in sponge:', rho_sur
   ! call MOM_mesg(mesg,5)
-  call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=US%kg_m3_to_R)
+  call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, US)
   !write (mesg,*) 'Bottom density in sponge:', rho_bot
   ! call MOM_mesg(mesg,5)
   rho_range = rho_bot - rho_sur

src/user/benchmark_initialization.F90

@@ -162,8 +162,8 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state
 
 ! Refine the guesses for each layer.
   do itt=1,6
-    call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US)
-    call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US)
+    call calculate_density(T0, S0, pres, rho_guess, eqn_of_state, US=US)
+    call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US=US)
     do k=1,nz
       T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
     enddo
@@ -258,7 +258,7 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
 
   T0(k1) = 29.0
   call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, US=US)
-  call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, US=US)
+  call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, US=US, dom=(/k1,k1/))
 
 ! A first guess of the layers' temperatures.                         !
   do k=1,nz