+Added dom interface to calculate_density

Removed the interfaces using G%HI to calculate_density and related routines and added a new variant with a new optional argument 'dom' specifying a two- element integer array with the start and end values of the domain to compute on for 1-d arrays starting at 1. If this array is not present in this new variant, calculations are done over the entire output array extent. All calls to the interfaces from before the rescale_pressure pull request will still work. Also added the new function EOS_domain that sets this two-element domain array from a horiz_index_type. This new interface is used throughout the code where the old, removed form was in use. All answers are bitwise identical.
NOAA-GFDL · Apr 17, 2020 · ab9d8e9 · ab9d8e9
1 parent 620a97c
commit ab9d8e9
Show file tree

Hide file tree

Showing 28 changed files with 337 additions and 334 deletions.
diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90
@@ -133,7 +133,7 @@ subroutine build_hycom1_column(CS, US, eqn_of_state, nz, depth, h, T, S, p_col,
   z_scale = 1.0 ; if (present(zScale)) z_scale = zScale
 
   ! Work bottom recording potential density
-  call calculate_density(T, S, p_col, rho_col, 1, nz, eqn_of_state, US=US)
+  call calculate_density(T, S, p_col, rho_col, eqn_of_state, US=US)
   ! This ensures the potential density profile is monotonic
   ! although not necessarily single valued.
   do k = nz-1, 1, -1

diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90
@@ -126,7 +126,7 @@ subroutine build_rho_column(CS, US, nz, depth, h, T, S, eqn_of_state, z_interfac
 
     ! Compute densities on source column
     pres(:) = CS%ref_pressure
-    call calculate_density(T, S, pres, densities, 1, nz, eqn_of_state, US=US)
+    call calculate_density(T, S, pres, densities, eqn_of_state, US=US)
     do k = 1,count_nonzero_layers
       densities(k) = densities(mapping(k))
     enddo

diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90
@@ -253,7 +253,7 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, &
     dz = (z_col(nz+1) - z_col(1)) / real(nz)
     do K=2,nz ; z_col_new(K) = z_col(1) + dz*real(K-1) ; enddo
   else
-    call calculate_density(T_col, S_col, p_col, rho_col, 1, nz, eqn_of_state, US=US)
+    call calculate_density(T_col, S_col, p_col, rho_col, eqn_of_state, US=US)
 
     ! Find the locations of the target potential densities, flagging
     ! locations in apparently unstable regions as not reliable.
@@ -374,10 +374,10 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, &
       enddo
       T_int(nz+1) = T_f(nz) ; S_int(nz+1) = S_f(nz)
       p_IS(nz+1) = z_col(nz+1) * H_to_pres
-      call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, 2, nz-1, &
-                                    eqn_of_state, US)
-      call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, &
-                                    eqn_of_state, US)
+      call calculate_density_derivs(T_int, S_int, p_IS, drhoIS_dT, drhoIS_dS, &
+                                    eqn_of_state, US, dom=(/2,nz/))
+      call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, &
+                                    eqn_of_state, US, dom=(/2,nz/))
       if (CS%compressibility_fraction > 0.0) then
         call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state, US)
       else

diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
@@ -74,7 +74,7 @@ module MOM
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
 use MOM_dyn_horgrid,           only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
-use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze
+use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze, EOS_domain
 use MOM_fixed_initialization,  only : MOM_initialize_fixed
 use MOM_forcing_type,          only : allocate_forcing_type, allocate_mech_forcing
 use MOM_forcing_type,          only : deallocate_mech_forcing, deallocate_forcing_type
@@ -2904,8 +2904,8 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
     ! Correct the output sea surface height for the contribution from the ice pressure.
     do j=js,je
       if (calc_rho) then
-        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, G%HI, &
-                               tv%eqn_of_state, US)
+        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), Rho_conv, &
+                               tv%eqn_of_state, US, dom=EOS_domain(G%HI))
         do i=is,ie
           IgR0 = US%Z_to_m / (Rho_conv(i) * GV%g_Earth)
           ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0

diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
@@ -10,7 +10,7 @@ module MOM_forcing_type
 use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_alloc, query_averaging_enabled
 use MOM_diag_mediator, only : enable_averages, enable_averaging, disable_averaging
 use MOM_error_handler, only : MOM_error, FATAL, WARNING
-use MOM_EOS,           only : calculate_density_derivs
+use MOM_EOS,           only : calculate_density_derivs, EOS_domain
 use MOM_file_parser,   only : get_param, log_param, log_version, param_file_type
 use MOM_grid,          only : ocean_grid_type
 use MOM_opacity,       only : sumSWoverBands, optics_type, extract_optics_slice, optics_nbands
@@ -953,8 +953,8 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt
                       H_limit_fluxes, .true., penSWbnd, netPen)
 
   ! Density derivatives
-  call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, G%HI, &
-                                tv%eqn_of_state, US)
+  call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, &
+                                tv%eqn_of_state, US, dom=EOS_domain(G%HI))
 
   ! Adjust netSalt to reflect dilution effect of FW flux
   netSalt(G%isc:G%iec) = netSalt(G%isc:G%iec) - Salt(G%isc:G%iec,j,1) * netH(G%isc:G%iec) ! ppt H/s

diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
@@ -15,7 +15,7 @@ module MOM_diagnostics
 use MOM_diag_mediator,     only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
 use MOM_domains,           only : create_group_pass, do_group_pass, group_pass_type
 use MOM_domains,           only : To_North, To_East
-use MOM_EOS,               only : calculate_density, int_density_dz
+use MOM_EOS,               only : calculate_density, int_density_dz, EOS_domain
 use MOM_EOS,               only : gsw_sp_from_sr, gsw_pt_from_ct
 use MOM_error_handler,     only : MOM_error, FATAL, WARNING
 use MOM_file_parser,       only : get_param, log_version, param_file_type
@@ -359,8 +359,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
             pressure_1d(i) = pressure_1d(i) + 0.5*(GV%H_to_RZ*GV%g_Earth)*h(i,j,k)
           enddo
           ! Store in-situ density [R ~> kg m-3] in work_3d
-          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rho_in_situ, G%HI, &
-                                 tv%eqn_of_state, US)
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rho_in_situ, &
+                                 tv%eqn_of_state, US, dom=EOS_domain(G%HI))
           do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d
             work_3d(i,j,k) = (GV%H_to_RZ*h(i,j,k)) / rho_in_situ(i)
           enddo
@@ -465,8 +465,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       pressure_1d(:) = tv%P_Ref
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js-1,je+1
-        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), G%HI, &
-                               tv%eqn_of_state, US, halo=1)
+        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), tv%eqn_of_state, &
+                               US, dom=EOS_domain(G%HI, halo=1))
       enddo ; enddo
     else ! Rcv should not be used much in this case, so fill in sensible values.
       do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
@@ -588,15 +588,17 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
       pressure_1d(:) = 0.
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js,je
-        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), G%HI, tv%eqn_of_state, US)
+        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
+                                tv%eqn_of_state, US, dom=EOS_domain(G%HI))
       enddo ; enddo
       if (CS%id_rhopot0 > 0) call post_data(CS%id_rhopot0, Rcv, CS%diag)
     endif
     if (CS%id_rhopot2 > 0) then
       pressure_1d(:) = 2.0e7*US%kg_m3_to_R*US%m_s_to_L_T**2 ! 2000 dbars
       !$OMP parallel do default(shared)
       do k=1,nz ; do j=js,je
-        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), G%HI, tv%eqn_of_state, US)
+        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
+                                tv%eqn_of_state, US, dom=EOS_domain(G%HI))
       enddo ; enddo
       if (CS%id_rhopot2 > 0) call post_data(CS%id_rhopot2, Rcv, CS%diag)
     endif
@@ -606,7 +608,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
         pressure_1d(:) = 0. ! Start at p=0 Pa at surface
         do k=1,nz
           pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure in middle of layer k
-          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), G%HI, tv%eqn_of_state, US)
+          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, Rcv(:,j,k), &
+                                tv%eqn_of_state, US, dom=EOS_domain(G%HI))
           pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure at bottom of layer k
         enddo
       enddo