From ab9d8e95ffc695f2d5c9feed388ec708c4563bf7 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 17 Apr 2020 11:05:11 -0400 Subject: [PATCH] +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. --- src/ALE/coord_hycom.F90 | 2 +- src/ALE/coord_rho.F90 | 2 +- src/ALE/coord_slight.F90 | 10 +- src/core/MOM.F90 | 6 +- src/core/MOM_forcing_type.F90 | 6 +- src/diagnostics/MOM_diagnostics.F90 | 19 +- src/equation_of_state/MOM_EOS.F90 | 399 +++++++++--------- src/ice_shelf/MOM_ice_shelf.F90 | 10 +- .../MOM_state_initialization.F90 | 16 +- .../lateral/MOM_mixed_layer_restrat.F90 | 16 +- .../lateral/MOM_thickness_diffuse.F90 | 10 +- .../vertical/MOM_CVMix_ddiff.F90 | 3 +- .../vertical/MOM_bulk_mixed_layer.F90 | 15 +- .../vertical/MOM_diabatic_aux.F90 | 23 +- .../vertical/MOM_diabatic_driver.F90 | 6 +- .../vertical/MOM_diapyc_energy_req.F90 | 2 +- .../vertical/MOM_energetic_PBL.F90 | 2 - .../vertical/MOM_entrain_diffusive.F90 | 15 +- .../vertical/MOM_full_convection.F90 | 15 +- .../vertical/MOM_internal_tide_input.F90 | 6 +- .../vertical/MOM_kappa_shear.F90 | 6 +- .../vertical/MOM_regularize_layers.F90 | 8 +- .../vertical/MOM_set_diffusivity.F90 | 21 +- .../vertical/MOM_set_viscosity.F90 | 18 +- src/tracer/MOM_neutral_diffusion.F90 | 16 +- src/tracer/MOM_tracer_hor_diff.F90 | 6 +- src/user/RGC_initialization.F90 | 5 +- src/user/user_change_diffusivity.F90 | 8 +- 28 files changed, 337 insertions(+), 334 deletions(-) diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index bfcff9005c..bc25377dbe 100644 --- 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 index 565656ecb0..7c6a00e714 100644 --- 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 index 000315bae8..de21e7027e 100644 --- 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 index efd4a80a52..4ca4682ef4 100644 --- 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 index 12f165372b..c54b26a6ab 100644 --- 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 index 94f6acc9c3..8d4c5dfa9b 100644 --- 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,7 +588,8 @@ 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 @@ -596,7 +597,8 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & 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 diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 9788c84338..b57677d8fa 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -40,7 +40,7 @@ module MOM_EOS public calculate_compress, calculate_density, query_compressible public calculate_density_derivs, calculate_specific_vol_derivs public calculate_density_second_derivs -public EOS_init, EOS_manual_init, EOS_end, EOS_allocate +public EOS_init, EOS_manual_init, EOS_end, EOS_allocate, EOS_domain public EOS_use_linear, calculate_spec_vol public int_density_dz, int_specific_vol_dp public int_density_dz_generic_plm, int_density_dz_generic_ppm @@ -59,23 +59,24 @@ module MOM_EOS !> Calculates density of sea water from T, S and P interface calculate_density - module procedure calculate_density_scalar, calculate_density_array, calculate_density_HI_1d + module procedure calculate_density_scalar, calculate_density_array, calculate_density_1d end interface calculate_density !> Calculates specific volume of sea water from T, S and P interface calculate_spec_vol - module procedure calc_spec_vol_scalar, calculate_spec_vol_array, calc_spec_vol_HI_1d, calc_spec_vol_US + module procedure calc_spec_vol_scalar, calculate_spec_vol_array, & + calc_spec_vol_1d end interface calculate_spec_vol !> Calculate the derivatives of density with temperature and salinity from T, S, and P interface calculate_density_derivs module procedure calculate_density_derivs_scalar, calculate_density_derivs_array, & - calculate_density_derivs_HI_1d + calculate_density_derivs_1d end interface calculate_density_derivs !> Calculate the derivatives of specific volume with temperature and salinity from T, S, and P interface calculate_specific_vol_derivs - module procedure calculate_spec_vol_derivs_array, calc_spec_vol_derivs_US, calc_spec_vol_derivs_HI_1d + module procedure calculate_spec_vol_derivs_array, calc_spec_vol_derivs_1d end interface calculate_specific_vol_derivs !> Calculates the second derivatives of density with various combinations of temperature, @@ -255,44 +256,64 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re end subroutine calculate_density_array -!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs -!! using array extents determined from a hor_index_type. +!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs, +!! potentially limiting the domain of indices that are worked on. !! If rho_ref is present, the anomaly with respect to rho_ref is returned. -subroutine calculate_density_HI_1d(T, S, pressure, rho, HI, EOS, US, halo) - type(hor_index_type), intent(in) :: HI !< The horizontal index structure - real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] - real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] - real, dimension(HI%isd:HI%ied), intent(inout) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3] - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. +subroutine calculate_density_1d(T, S, pressure, rho, EOS, US, dom, rho_ref, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3] + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density + !! in combination with scaling given by US [various] ! Local variables - real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] - integer :: i, is, ie, start, npts, halo_sz + real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1] + real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1] + real :: rho_unscale ! A factor to convert density from R to kg m-3 [kg m-3 R-1 ~> 1] + real :: rho_reference ! rho_ref converted to [kg m-3] + real, dimension(size(rho)) :: pres ! Pressure converted to [Pa] + integer :: i, is, ie, npts if (.not.associated(EOS)) call MOM_error(FATAL, & - "calculate_density_HI_1d called with an unassociated EOS_type EOS.") + "calculate_density_1d called with an unassociated EOS_type EOS.") - halo_sz = 0 ; if (present(halo)) halo_sz = halo + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is + else + is = 1 ; ie = size(rho) ; npts = 1 + ie - is + endif - start = HI%isc - (HI%isd-1) - halo_sz - npts = HI%iec - HI%isc + 1 + 2*halo_sz - is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + rho_unscale = 1.0 ; if (present(US)) rho_unscale = US%R_to_kg_m3 + + if ((p_scale == 1.0) .and. (rho_unscale == 1.0)) then + call calculate_density_array(T, S, pressure, rho, is, npts, EOS, rho_ref=rho_ref) + elseif (present(rho_ref)) then ! This is the same as above, but with some extra work to rescale variables. + do i=is,ie ; pres(i) = p_scale * pressure(i) ; enddo + rho_reference = rho_unscale*rho_ref + call calculate_density_array(T, S, pres, rho, is, npts, EOS, rho_ref=rho_reference) + else ! There is rescaling of variables, but rho_ref is not present. Passing a 0 value of rho_ref + ! changes answers at roundoff for some equations of state, like Wright and UNESCO. + do i=is,ie ; pres(i) = p_scale * pressure(i) ; enddo + call calculate_density_array(T, S, pres, rho, is, npts, EOS) + endif - if (US%RL2_T2_to_Pa == 1.0) then - call calculate_density_array(T, S, pressure, rho, start, npts, EOS) - else ! There is rescaling of variables, including pressure. - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_density_array(T, S, pres, rho, start, npts, EOS) + if (present(US) .or. present(scale)) then + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + if (present(scale)) rho_scale = rho_scale * scale + if (rho_scale /= 1.0) then ; do i=is,ie + rho(i) = rho_scale * rho(i) + enddo ; endif endif - if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie - rho(i) = US%kg_m3_to_R * rho(i) - enddo ; endif +end subroutine calculate_density_1d -end subroutine calculate_density_HI_1d !> Calls the appropriate subroutine to calculate the specific volume of sea water !! for 1-D array inputs. @@ -380,97 +401,64 @@ subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, US, scale end subroutine calc_spec_vol_scalar -!> Calls the appropriate subroutine to calculate the specific volume of sea water -!! for 1-D array inputs with dimensional rescaling. -subroutine calc_spec_vol_US(T, S, pressure, specvol, start, npts, EOS, US, spv_ref, scale) - real, dimension(:), intent(in) :: T !< potential temperature relative to the surface [degC] - real, dimension(:), intent(in) :: S !< salinity [ppt] - real, dimension(:), intent(in) :: pressure !< pressure [Pa] or [R L2 T-2 ~> Pa] - real, dimension(:), intent(inout) :: specvol !< in situ specific volume [kg m-3] or [R-1 ~> m3 kg-1] - integer, intent(in) :: start !< the starting point in the arrays. - integer, intent(in) :: npts !< the number of values to calculate. - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1] - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific - !! volume in combination with scaling given by US [various] - - ! Local variables - real, dimension(size(pressure)) :: pres ! Pressure converted to [Pa] - real :: spv_reference ! spv_ref converted to [m3 kg-1] - real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg R-1 m-3 ~> 1] - integer :: i, is, ie - - if (.not.associated(EOS)) call MOM_error(FATAL, & - "calculate_spec_vol_array called with an unassociated EOS_type EOS.") - - is = start ; ie = is + npts - 1 - - if ((US%RL2_T2_to_Pa == 1.0) .and. (US%R_to_kg_m3 == 1.0)) then - call calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref) - elseif (present(spv_ref)) then ! This is the same as above, but with some extra work to rescale variables. - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - spv_reference = US%kg_m3_to_R*spv_ref - call calculate_spec_vol_array(T, S, pres, specvol, start, npts, EOS, spv_reference) - else ! There is rescaling of variables, but spv_ref is not present. Passing a 0 value of spv_ref - ! changes answers at roundoff for some equations of state, like Wright and UNESCO. - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_spec_vol_array(T, S, pres, specvol, start, npts, EOS) - endif - - spv_scale = US%R_to_kg_m3 - if (present(scale)) spv_scale = spv_scale * scale - if (spv_scale /= 1.0) then ; do i=is,ie - specvol(i) = spv_scale * specvol(i) - enddo ; endif - -end subroutine calc_spec_vol_US - - !> Calls the appropriate subroutine to calculate the specific volume of sea water for 1-D array -!! inputs using array extents determined from a hor_index_type. -subroutine calc_spec_vol_HI_1d(T, S, pressure, specvol, HI, EOS, US, halo, spv_ref) - type(hor_index_type), intent(in) :: HI !< The horizontal index structure - real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] - real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] - real, dimension(HI%isd:HI%ied), intent(inout) :: specvol !< In situ specific volume [R-1 ~> m3 kg-1] - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. - real, optional, intent(in) :: spv_ref !< A reference specific volume [R-1 ~> m3 kg-1] +!! inputs, potentially limiting the domain of indices that are worked on. +subroutine calc_spec_vol_1d(T, S, pressure, specvol, EOS, US, dom, spv_ref, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: specvol !< In situ specific volume [R-1 ~> m3 kg-1] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: spv_ref !< A reference specific volume [R-1 ~> m3 kg-1] + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale + !! output specific volume in combination with + !! scaling given by US [various] ! Local variables - real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] + real, dimension(size(specvol)) :: pres ! Pressure converted to [Pa] + real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1] + real :: spv_unscale ! A factor to convert specific volume from R-1 to m3 kg-1 [m3 kg-1 R ~> 1] + real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] real :: spv_reference ! spv_ref converted to [m3 kg-1] - integer :: i, is, ie, start, npts, halo_sz + integer :: i, is, ie, npts if (.not.associated(EOS)) call MOM_error(FATAL, & - "calc_spec_vol_HI_1d called with an unassociated EOS_type EOS.") + "calc_spec_vol_1d called with an unassociated EOS_type EOS.") - halo_sz = 0 ; if (present(halo)) halo_sz = halo + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is + else + is = 1 ; ie = size(specvol) ; npts = 1 + ie - is + endif - start = HI%isc - (HI%isd-1) - halo_sz - npts = HI%iec - HI%isc + 1 + 2*halo_sz - is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + spv_unscale = 1.0 ; if (present(US)) spv_unscale = US%kg_m3_to_R - if ((US%RL2_T2_to_Pa == 1.0) .and. (US%R_to_kg_m3 == 1.0)) then - call calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref) + if ((p_scale == 1.0) .and. (spv_unscale == 1.0)) then + call calculate_spec_vol_array(T, S, pressure, specvol, is, npts, EOS, spv_ref) elseif (present(spv_ref)) then ! This is the same as above, but with some extra work to rescale variables. - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - spv_reference = US%kg_m3_to_R*spv_ref - call calculate_spec_vol_array(T, S, pres, specvol, start, npts, EOS, spv_reference) + do i=is,ie ; pres(i) = p_scale * pressure(i) ; enddo + spv_reference = spv_unscale*spv_ref + call calculate_spec_vol_array(T, S, pres, specvol, is, npts, EOS, spv_reference) else ! There is rescaling of variables, but spv_ref is not present. Passing a 0 value of spv_ref ! changes answers at roundoff for some equations of state, like Wright and UNESCO. - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_spec_vol_array(T, S, pres, specvol, start, npts, EOS) + do i=is,ie ; pres(i) = p_scale * pressure(i) ; enddo + call calculate_spec_vol_array(T, S, pres, specvol, is, npts, EOS) endif - if (US%R_to_kg_m3 /= 1.0) then ; do i=is,ie - specvol(i) = US%R_to_kg_m3 * specvol(i) - enddo ; endif + if (present(US) .or. present(scale)) then + spv_scale = 1.0 ; if (present(US)) spv_scale = US%R_to_kg_m3 + if (present(scale)) spv_scale = spv_scale * scale + if (spv_scale /= 1.0) then ; do i=is,ie + specvol(i) = spv_scale * specvol(i) + enddo ; endif + endif + +end subroutine calc_spec_vol_1d -end subroutine calc_spec_vol_HI_1d !> Calls the appropriate subroutine to calculate the freezing point for scalar inputs. subroutine calculate_TFreeze_scalar(S, pressure, T_fr, EOS, pres_scale) @@ -626,45 +614,55 @@ end subroutine calculate_density_derivs_array !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. -subroutine calculate_density_derivs_HI_1d(T, S, pressure, drho_dT, drho_dS, HI, EOS, US, halo) - type(hor_index_type), intent(in) :: HI !< The horizontal index structure - real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] - real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] - real, dimension(HI%isd:HI%ied), intent(inout) :: drho_dT !< The partial derivative of density with potential - !! temperature [R degC-1 ~> kg m-3 degC-1] - real, dimension(HI%isd:HI%ied), intent(inout) :: drho_dS !< The partial derivative of density with salinity - !! [R degC-1 ~> kg m-3 ppt-1] - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. +subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, US, dom, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: drho_dT !< The partial derivative of density with potential + !! temperature [R degC-1 ~> kg m-3 degC-1] + real, dimension(:), intent(inout) :: drho_dS !< The partial derivative of density with salinity + !! [R degC-1 ~> kg m-3 ppt-1] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density + !! in combination with scaling given by US [various] ! Local variables - real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] - integer :: i, is, ie, start, npts, halo_sz + real, dimension(size(drho_dT)) :: pres ! Pressure converted to [Pa] + real :: rho_scale ! A factor to convert density from kg m-3 to the desired units [R m3 kg-1 ~> 1] + real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1] + integer :: i, is, ie, npts if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - halo_sz = 0 ; if (present(halo)) halo_sz = halo + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is + else + is = 1 ; ie = size(drho_dT) ; npts = 1 + ie - is + endif - start = HI%isc - (HI%isd-1) - halo_sz - npts = HI%iec - HI%isc + 1 + 2*halo_sz - is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa - if (US%RL2_T2_to_Pa == 1.0) then - call calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS) + if (p_scale == 1.0) then + call calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, is, npts, EOS) else - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_density_derivs_array(T, S, pres, drho_dT, drho_dS, start, npts, EOS) + do i=is,ie ; pres(i) = p_scale * pressure(i) ; enddo + call calculate_density_derivs_array(T, S, pres, drho_dT, drho_dS, is, npts, EOS) endif - if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie - drho_dT(i) = US%kg_m3_to_R * drho_dT(i) - drho_dS(i) = US%kg_m3_to_R * drho_dS(i) - enddo ; endif + if (present(US) .or. present(scale)) then + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + if (present(scale)) rho_scale = rho_scale * scale + if (rho_scale /= 1.0) then ; do i=is,ie + drho_dT(i) = rho_scale * drho_dT(i) + drho_dS(i) = rho_scale * drho_dS(i) + enddo ; endif + endif -end subroutine calculate_density_derivs_HI_1d +end subroutine calculate_density_derivs_1d !> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar @@ -912,90 +910,56 @@ subroutine calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start end subroutine calculate_spec_vol_derivs_array - -!> Calls the appropriate subroutine to calculate specific volume derivatives for an array with unit scaling. -subroutine calc_spec_vol_derivs_US(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, US, scale) - real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(:), intent(in) :: S !< Salinity [ppt] - real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] - real, dimension(:), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential - !! temperature [R-1 degC-1 ~> m3 kg-1 degC-1] - real, dimension(:), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity - !! [R-1 ppt-1 ~> m3 kg-1 ppt-1] - integer, intent(in) :: start !< Starting index within the array - integer, intent(in) :: npts !< The number of values to calculate - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific - !! volume in combination with scaling given by US [various] +!> Calls the appropriate subroutine to calculate specific volume derivatives for 1-d array inputs, +!! potentially limiting the domain of indices that are worked on. +subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, US, dom, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(:), intent(in) :: S !< Salinity [ppt] + real, dimension(:), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential + !! temperature [R-1 degC-1 ~> m3 kg-1 degC-1] + real, dimension(:), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity + !! [R-1 ppt-1 ~> m3 kg-1 ppt-1] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific + !! volume in combination with scaling given by US [various] ! Local variables - real, dimension(size(T)) :: press ! Pressure converted to [Pa] + real, dimension(size(dSV_dT)) :: press ! Pressure converted to [Pa] real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg R-1 m-3 ~> 1] - integer :: i, is, ie + real :: p_scale ! A factor to convert pressure to units of Pa [Pa T2 R-1 L-2 ~> 1] + integer :: i, is, ie, npts if (.not.associated(EOS)) call MOM_error(FATAL, & - "calculate_spec_vol_derivs_array called with an unassociated EOS_type EOS.") - - is = start ; ie = is + npts - 1 + "calculate_spec_vol_derivs_1d called with an unassociated EOS_type EOS.") - if (US%RL2_T2_to_Pa == 1.0) then - call calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS) + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is else - do i=is,ie ; press(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_spec_vol_derivs_array(T, S, press, dSV_dT, dSV_dS, start, npts, EOS) + is = 1 ; ie = size(dSV_dT) ; npts = 1 + ie - is endif + p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa - spv_scale = US%R_to_kg_m3 - if (present(scale)) spv_scale = spv_scale * scale - if (spv_scale /= 1.0) then ; do i=is,ie - dSV_dT(i) = spv_scale * dSV_dT(i) - dSV_dS(i) = spv_scale * dSV_dS(i) - enddo ; endif - -end subroutine calc_spec_vol_derivs_US - -!> Calls the appropriate subroutine to calculate specific volume derivatives for array inputs -!! using array extents determined from a hor_index_type.. -subroutine calc_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, US, halo) - type(hor_index_type), intent(in) :: HI !< The horizontal index structure - real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] - real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] - real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] - real, dimension(HI%isd:HI%ied), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential - !! temperature [R-1 degC-1 ~> m3 kg-1 degC-1] - real, dimension(HI%isd:HI%ied), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity - !! [R-1 ppt-1 ~> m3 kg-1 ppt-1] - type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. - - ! Local variables - real, dimension(HI%isd:HI%ied) :: press ! Pressure converted to [Pa] - integer :: i, is, ie, start, npts, halo_sz - - if (.not.associated(EOS)) call MOM_error(FATAL, & - "calculate_spec_vol_derivs_HI_1d called with an unassociated EOS_type EOS.") - - halo_sz = 0 ; if (present(halo)) halo_sz = halo - - start = HI%isc - (HI%isd-1) - halo_sz - npts = HI%iec - HI%isc + 1 + 2*halo_sz - is = HI%isc - halo_sz ; ie = HI%iec + halo_sz - - if (US%RL2_T2_to_Pa == 1.0) then - call calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS) + if (p_scale == 1.0) then + call calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, is, npts, EOS) else - do i=is,ie ; press(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo - call calculate_spec_vol_derivs_array(T, S, press, dSV_dT, dSV_dS, start, npts, EOS) + do i=is,ie ; press(i) = p_scale * pressure(i) ; enddo + call calculate_spec_vol_derivs_array(T, S, press, dSV_dT, dSV_dS, is, npts, EOS) endif - if (US%R_to_kg_m3 /= 1.0) then ; do i=is,ie - dSV_dT(i) = US%R_to_kg_m3 * dSV_dT(i) - dSV_dS(i) = US%R_to_kg_m3 * dSV_dS(i) - enddo ; endif + if (present(US) .or. present(scale)) then + spv_scale = 1.0 ; if (present(US)) spv_scale = US%R_to_kg_m3 + if (present(scale)) spv_scale = spv_scale * scale + if (spv_scale /= 1.0) then ; do i=is,ie + dSV_dT(i) = spv_scale * dSV_dT(i) + dSV_dS(i) = spv_scale * dSV_dS(i) + enddo ; endif + endif -end subroutine calc_spec_vol_derivs_HI_1d +end subroutine calc_spec_vol_derivs_1d !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array @@ -1079,12 +1043,31 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS, US) end subroutine calculate_compress_scalar +!> This subroutine returns a two point integer array indicating the domain of i-indices +!! to work on in EOS calls based on information from a hor_index type +function EOS_domain(HI, halo) result(EOSdom) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + integer, dimension(2) :: EOSdom !< The index domain that the EOS will work on, taking into account + !! that the arrays inside the EOS routines will start at 1. + + ! Local variables + integer :: halo_sz + + halo_sz = 0 ; if (present(halo)) halo_sz = halo + + EOSdom(1) = HI%isc - (HI%isd-1) - halo_sz + EOSdom(2) = HI%iec - (HI%isd-1) + halo_sz + +end function EOS_domain + + !> Calls the appropriate subroutine to calculate analytical and nearly-analytical !! integrals in pressure across layers of geopotential anomalies, which are !! required for calculating the finite-volume form pressure accelerations in a !! non-Boussinesq model. There are essentially no free assumptions, apart from the !! use of Bode's rule to do the horizontal integrals, and from a truncation in the -!! series for log(1-eps/1+eps) that assumes that |eps| < . +!! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & bathyP, dP_tiny, useMassWghtInterp, US) @@ -3039,7 +3022,7 @@ subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS) intent(in) :: mask_z !< 3d mask regulating which points to convert. type(EOS_type), pointer :: EOS !< Equation of state structure - integer :: i,j,k + integer :: i, j, k real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp real :: p diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 858af4e1ea..8f1587ab8b 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -35,7 +35,7 @@ module MOM_ice_shelf use MOM_forcing_type, only : mech_forcing, allocate_mech_forcing, MOM_mech_forcing_chksum use MOM_forcing_type, only : copy_common_forcing_fields use MOM_get_input, only : directories, Get_MOM_input -use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_TFreeze +use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_TFreeze, EOS_domain use MOM_EOS, only : EOS_type, EOS_init use MOM_ice_shelf_dynamics, only : ice_shelf_dyn_CS, update_ice_shelf use MOM_ice_shelf_dynamics, only : register_ice_shelf_dyn_restarts, initialize_ice_shelf_dyn @@ -375,10 +375,10 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) 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(state%sst(:,j), state%sss(:,j), p_int, Rhoml(:), G%HI, & - CS%eqn_of_state, US) - call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, dR0_dT, dR0_dS, G%HI, & - CS%eqn_of_state, US) + call calculate_density(state%sst(:,j), state%sss(:,j), p_int, Rhoml(:), & + CS%eqn_of_state, US, dom=EOS_domain(G%HI)) + call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, dR0_dT, dR0_dS, & + CS%eqn_of_state, US, dom=EOS_domain(G%HI)) do i=is,ie if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2efceb5991..2ac8ac47bf 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -41,7 +41,7 @@ module MOM_state_initialization use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type use MOM_ALE, only : pressure_gradient_plm -use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type, EOS_domain use MOM_EOS, only : int_specific_vol_dp, convert_temp_salt_for_TEOS10 use user_initialization, only : user_initialize_thickness, user_initialize_velocity use user_initialization, only : user_init_temperature_salinity @@ -959,8 +959,8 @@ subroutine convert_thickness(h, G, GV, US, tv) do k=1,nz do j=js,je do i=is,ie ; p_top(i,j) = p_bot(i,j) ; enddo - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, G%HI, & - tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) do i=is,ie p_bot(i,j) = p_top(i,j) + HR_to_pres * (h(i,j,k) * rho(i)) enddo @@ -970,8 +970,8 @@ subroutine convert_thickness(h, G, GV, US, tv) call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, 0.0, G%HI, & tv%eqn_of_state, dz_geo, US=US) if (itt < max_itt) then ; do j=js,je - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, G%HI, & - tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_bot(:,j), rho, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) ! Use Newton's method to correct the bottom value. ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed. do i=is,ie @@ -1869,7 +1869,8 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, C call MOM_read_data(filename, salin_var, tmp2(:,:,:), G%Domain) do j=js,je - call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), G%HI, tv%eqn_of_state, US) + call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) enddo call set_up_sponge_ML_density(tmp_2d, G, CSp) @@ -2188,7 +2189,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param press(:) = tv%P_Ref do k=1,kd ; do j=js,je - call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), G%HI, eos, US) + call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), eos, US, & + dom=EOS_domain(G%HI)) enddo ; enddo call pass_var(temp_z,G%Domain) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 744a801391..fa1504b431 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -19,7 +19,7 @@ module MOM_mixed_layer_restrat use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density +use MOM_EOS, only : calculate_density, EOS_domain implicit none ; private @@ -207,8 +207,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var pRef_MLD(:) = 0. do j = js-1, je+1 dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, G%HI, & - tv%eqn_of_state, US, halo=1) + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, & + US, dom=EOS_domain(G%HI, halo=1)) deltaRhoAtK(:) = 0. MLD_fast(:,j) = 0. do k = 2, nz @@ -216,8 +216,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K ! Mixed-layer depth, using sigma-0 (surface reference pressure) deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, G%HI, & - tv%eqn_of_state, US, halo=1) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, & + US, dom=EOS_domain(G%HI, halo=1)) do i = is-1,ie+1 deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface enddo @@ -322,7 +322,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var h_avail(i,j,k) = max(I4dt*G%areaT(i,j)*(h(i,j,k)-GV%Angstrom_H),0.0) enddo if (keep_going) then - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), G%HI, tv%eqn_of_state, US, halo=1) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, & + US, dom=EOS_domain(G%HI, halo=1)) line_is_empty = .true. do i=is-1,ie+1 if (htot_fast(i,j) < MLD_fast(i,j)) then @@ -646,7 +647,8 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) htot(i,j) = 0.0 ; Rml_av(i,j) = 0.0 enddo do k=1,nkml - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, Rho0(:), G%HI, tv%eqn_of_state, US, halo=1) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, Rho0(:), tv%eqn_of_state, & + US, dom=EOS_domain(G%HI, halo=1)) do i=is-1,ie+1 Rml_av(i,j) = Rml_av(i,j) + h(i,j,k)*Rho0(i) htot(i,j) = htot(i,j) + h(i,j,k) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 4d42f05629..62f1dad445 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -9,7 +9,7 @@ module MOM_thickness_diffuse use MOM_diag_mediator, only : diag_update_remap_grids use MOM_domains, only : pass_var, CORNER, pass_vector use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe -use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_interface_heights, only : find_eta @@ -1029,8 +1029,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1))) S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1))) enddo - call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, G%HI, & - tv%eqn_of_state, US) + call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) endif do i=is,ie if (calc_derivatives) then @@ -1291,8 +1291,8 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV T_v(i) = 0.5*(T(i,j,1) + T(i,j+1,1)) S_v(i) = 0.5*(S(i,j,1) + S(i,j+1,1)) enddo - call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, G%HI, & - tv%eqn_of_state, US) + call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v, drho_dS_v, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) endif do i=is,ie vhD(i,J,1) = -vhtot(i,J) diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index 1bcb6a1266..f169147d03 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -236,8 +236,7 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS) pRef = pRef + GV%H_to_Pa * h(i,j,k-1) enddo ! k-loop finishes - call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, & - G%ke, TV%EQN_OF_STATE) + call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dT, drho_dS, tv%eqn_of_state) ! The "-1.0" below is needed so that the following criteria is satisfied: ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger" diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index cc3a6e3f69..2a6dd66c20 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -15,7 +15,7 @@ module MOM_bulk_mixed_layer use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain implicit none ; private @@ -466,12 +466,15 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C do k=1,CS%nkml ; do i=is,ie p_ref(i) = p_ref(i) + 0.5*(GV%H_to_RZ*GV%g_Earth)*h(i,k) enddo ; enddo - call calculate_density_derivs(T(:,1), S(:,1), p_ref, dR0_dT, dR0_dS, G%HI, tv%eqn_of_state, US) - call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, G%HI, & - tv%eqn_of_state, US) + call calculate_density_derivs(T(:,1), S(:,1), p_ref, dR0_dT, dR0_dS, tv%eqn_of_state, & + US, dom=EOS_domain(G%HI)) + call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, tv%eqn_of_state, & + US, dom=EOS_domain(G%HI)) do k=1,nz - call calculate_density(T(:,k), S(:,k), p_ref, R0(:,k), G%HI, tv%eqn_of_state, US) - call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state, US) + call calculate_density(T(:,k), S(:,k), p_ref, R0(:,k), tv%eqn_of_state, & + US, dom=EOS_domain(G%HI)) + call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), tv%eqn_of_state, & + US, dom=EOS_domain(G%HI)) enddo if (id_clock_EOS>0) call cpu_clock_end(id_clock_EOS) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 3cde9ce91e..ae7adc05eb 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -9,7 +9,7 @@ module MOM_diabatic_aux use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl, time_type -use MOM_EOS, only : calculate_density, calculate_TFreeze +use MOM_EOS, only : calculate_density, calculate_TFreeze, EOS_domain use MOM_EOS, only : calculate_specific_vol_derivs, calculate_density_derivs use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint @@ -446,7 +446,8 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay) h_2d(i,k) = MAX(h(i,j,k), GV%Angstrom_H) enddo - call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state, US) + call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) enddo ! First, try to find an interior layer where inserting all the salt @@ -766,7 +767,8 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, pRef_MLD(:) = 0.0 do j=js,je do i=is,ie ; dK(i) = 0.5 * h(i,j,1) * GV%H_to_Z ; enddo ! Depth of center of surface layer - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, G%HI, tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) do i=is,ie deltaRhoAtK(i) = 0. MLD(i,j) = 0. @@ -807,7 +809,8 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, ! Mixed-layer depth, using sigma-0 (surface reference pressure) do i=is,ie ; deltaRhoAtKm1(i) = deltaRhoAtK(i) ; enddo ! Store value from previous iteration of K - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, G%HI, tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) do i = is, ie deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i) @@ -830,8 +833,10 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, ! T_deeper(i) = tv%T(i,j,nz) ; S_deeper(i) = tv%S(i,j,nz) ! N2_region_set(i) = .true. ! endif - call calculate_density(T_subML, S_subML, pRef_N2, rho_subML, G%HI, tv%eqn_of_state, US) - call calculate_density(T_deeper, S_deeper, pRef_N2, rho_deeper, G%HI, tv%eqn_of_state, US) + call calculate_density(T_subML, S_subML, pRef_N2, rho_subML, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) + call calculate_density(T_deeper, S_deeper, pRef_N2, rho_deeper, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) do i=is,ie ; if ((G%mask2dT(i,j)>0.5) .and. N2_region_set(i)) then subMLN2(i,j) = gE_rho0 * (rho_deeper(i) - rho_subML(i)) / (GV%H_to_z * dH_N2(i)) endif ; enddo @@ -1006,7 +1011,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t pres(i) = pres(i) + d_pres(i) enddo call calculate_specific_vol_derivs(T2d(:,k), tv%S(:,j,k), p_lay(:), & - dSV_dT(:,j,k), dSV_dS(:,j,k), G%HI, tv%eqn_of_state, US=US) + dSV_dT(:,j,k), dSV_dS(:,j,k), tv%eqn_of_state, US=US, dom=EOS_domain(G%HI)) do i=is,ie ; dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; enddo enddo pen_TKE_2d(:,:) = 0.0 @@ -1347,8 +1352,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t do i=is,ie ; do nb=1,nsw ; netPen_rate(i) = netPen_rate(i) + pen_SW_bnd_rate(nb,i) ; enddo ; enddo ! Density derivatives - call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dRhodT, dRhodS, G%HI, & - tv%eqn_of_state, US) + call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dRhodT, dRhodS, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) ! 1. Adjust netSalt to reflect dilution effect of FW flux ! 2. Add in the SW heating for purposes of calculating the net ! surface buoyancy flux affecting the top layer. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index bdfb6b7a9e..607d02722f 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -34,7 +34,7 @@ module MOM_diabatic_driver use MOM_energetic_PBL, only : energetic_PBL_get_MLD use MOM_entrain_diffusive, only : entrainment_diffusive, entrain_diffusive_init use MOM_entrain_diffusive, only : entrain_diffusive_end, entrain_diffusive_CS -use MOM_EOS, only : calculate_density, calculate_TFreeze +use MOM_EOS, only : calculate_density, calculate_TFreeze, EOS_domain use MOM_error_handler, only : MOM_error, FATAL, WARNING, callTree_showQuery,MOM_mesg use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_version, param_file_type, read_param @@ -2682,8 +2682,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo !$OMP parallel do default(shared) do j=js,je - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), G%HI, & - tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) enddo call apply_sponge(h, dt, G, GV, US, ea, eb, CS%sponge_CSp, Rcv_ml) else diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index cde4b9e484..edad667592 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -298,7 +298,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & ! Solve the tridiagonal equations for new temperatures. - call calculate_specific_vol_derivs(T0, S0, p_lay, dSV_dT, dSV_dS, 1, nz, tv%eqn_of_state, US=US) + call calculate_specific_vol_derivs(T0, S0, p_lay, dSV_dT, dSV_dS, tv%eqn_of_state, US=US) do k=1,nz dMass = GV%H_to_RZ * h_tr(k) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 8b5da8565b..48b265a0e2 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -17,8 +17,6 @@ module MOM_energetic_PBL use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only: wave_parameters_CS, Get_Langmuir_Number -! use MOM_EOS, only : calculate_density, calculate_density_derivs - implicit none ; private #include diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 100e79aba2..7e8d306ff0 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -12,7 +12,7 @@ module MOM_entrain_diffusive use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density, calculate_density_derivs +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain implicit none ; private @@ -700,7 +700,8 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & call determine_dSkb(h_bl, Sref, Ent_bl, eakb, is, ie, kmb, G, GV, & .true., dS_kb, dS_anom_lim=dS_anom_lim) do k=nz-1,kb_min,-1 - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, G%HI, tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) do i=is,ie if ((k>kb(i)) .and. (F(i,k) > 0.0)) then ! Within a time step, a layer may entrain no more than its @@ -784,7 +785,8 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & else ! not bulkmixedlayer do k=K2,nz-1; - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, G%HI, tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, Rcv, tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) do i=is,ie ; if (F(i,k) > 0.0) then ! Within a time step, a layer may entrain no more than ! its thickness for correction. This limitation should @@ -849,8 +851,8 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & S_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) endif enddo - call calculate_density_derivs(T_eos, S_eos, pressure, dRho_dT, dRho_dS, G%HI, & - tv%eqn_of_state, US) + call calculate_density_derivs(T_eos, S_eos, pressure, dRho_dT, dRho_dS, & + tv%eqn_of_state, US, dom=EOS_domain(G%HI)) do i=is,ie if ((k>kmb) .and. (k @@ -213,7 +213,8 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) do i=is-1,ie ; pres(i) = tv%P_Ref ; enddo do j=js,je - call calculate_density(T(:,j,1), S(:,j,1), pres, tmp(:,j), G%HI, tv%eqn_of_state, US) + call calculate_density(T(:,j,1), S(:,j,1), pres, tmp(:,j), tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) enddo call set_up_sponge_ML_density(tmp, G, CSp) diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index f33d772352..92b17e13fb 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -10,7 +10,7 @@ module user_change_diffusivity use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density +use MOM_EOS, only : calculate_density, EOS_domain implicit none ; private @@ -107,11 +107,13 @@ subroutine user_change_diff(h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_i do j=js,je if (present(T_f) .and. present(S_f)) then do k=1,nz - call calculate_density(T_f(:,j,k), S_f(:,j,k), p_ref, Rcv(:,k), G%HI, tv%eqn_of_state, US) + call calculate_density(T_f(:,j,k), S_f(:,j,k), p_ref, Rcv(:,k), tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) enddo else do k=1,nz - call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rcv(:,k), G%HI, tv%eqn_of_state, US) + call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rcv(:,k), tv%eqn_of_state, US, & + dom=EOS_domain(G%HI)) enddo endif