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