From 415a6bce2e95fd55cd84e21987a52f53a34a677c Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Apr 2020 15:48:22 +0000 Subject: [PATCH 1/4] Adds unit conversions to EOS type - As part of the process of disconnected EOS from the unit_scaling_type this adds the necessary unit conversions to the EOS_type. - Initialization is currently donne by passing US to MOM_init() but ultimately it seems passing p_scaling, etc., to MOM_init() would remove all dependency on US. - No APIs other than EOS_init() have been changed yet. --- src/core/MOM.F90 | 2 +- src/equation_of_state/MOM_EOS.F90 | 130 +++++++++++++++++------------- 2 files changed, 73 insertions(+), 59 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index efd4a80a52..f07fd6a1c4 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2242,7 +2242,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Use the Wright equation of state by default, unless otherwise specified ! Note: this line and the following block ought to be in a separate ! initialization routine for tv. - if (use_EOS) call EOS_init(param_file, CS%tv%eqn_of_state) + if (use_EOS) call EOS_init(param_file, CS%tv%eqn_of_state, US) if (use_temperature) then allocate(CS%tv%TempxPmE(isd:ied,jsd:jed)) ; CS%tv%TempxPmE(:,:) = 0.0 if (use_geothermal) then diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 9788c84338..dc74e1dcf7 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -112,6 +112,13 @@ module MOM_EOS real :: dTFr_dS !< The derivative of freezing point with salinity [degC ppt-1] real :: dTFr_dp !< The derivative of freezing point with pressure [degC Pa-1] +! Unit conversion factors (normally used for dimensional testing but could also allow for +! change of units of arguments to functions) + real :: m_to_Z !< A constant that translates distances in meters to the units of depth. + real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density. + real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed. + real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa. + ! logical :: test_EOS = .true. ! If true, test the equation of state end type EOS_type @@ -161,7 +168,7 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, US, scale if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_scalar called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -180,7 +187,7 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, US, scale end select if (present(US) .or. present(scale)) then - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale rho = rho_scale * rho endif @@ -210,7 +217,7 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_array called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -247,7 +254,7 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then do j=start,start+npts-1 ; rho(j) = rho_scale * rho(j) ; enddo @@ -281,15 +288,15 @@ subroutine calculate_density_HI_1d(T, S, pressure, rho, HI, EOS, US, halo) 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 + if (EOS%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 + do i=is,ie ; pres(i) = EOS%RL2_T2_to_Pa * pressure(i) ; enddo call calculate_density_array(T, S, pres, rho, start, npts, EOS) endif - if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie - rho(i) = US%kg_m3_to_R * rho(i) + if (EOS%kg_m3_to_R /= 1.0) then ; do i=is,ie + rho(i) = EOS%kg_m3_to_R * rho(i) enddo ; endif end subroutine calculate_density_HI_1d @@ -361,18 +368,18 @@ subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, US, scale if (.not.associated(EOS)) call MOM_error(FATAL, & "calc_spec_vol_scalar called with an unassociated EOS_type EOS.") - pres(1) = pressure ; if (present(US)) pres(1) = US%RL2_T2_to_Pa*pressure + pres(1) = pressure ; if (present(US)) pres(1) = EOS%RL2_T2_to_Pa*pressure Ta(1) = T ; Sa(1) = S if (present(spv_ref)) then - spv_reference = spv_ref ; if (present(US)) spv_reference = US%kg_m3_to_R*spv_ref + spv_reference = spv_ref ; if (present(US)) spv_reference = EOS%kg_m3_to_R*spv_ref call calculate_spec_vol_array(Ta, Sa, pres, spv, 1, 1, EOS, spv_reference) else call calculate_spec_vol_array(Ta, Sa, pres, spv, 1, 1, EOS) endif specvol = spv(1) - spv_scale = 1.0 ; if (present(US)) spv_scale = US%R_to_kg_m3 + spv_scale = 1.0 ; if (present(US)) spv_scale = EOS%R_to_kg_m3 if (present(scale)) spv_scale = spv_scale * scale if (spv_scale /= 1.0) then specvol = spv_scale * specvol @@ -454,20 +461,20 @@ subroutine calc_spec_vol_HI_1d(T, S, pressure, specvol, HI, EOS, US, halo, spv_r 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) .and. (US%R_to_kg_m3 == 1.0)) then + if ((EOS%RL2_T2_to_Pa == 1.0) .and. (EOS%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 + do i=is,ie ; pres(i) = EOS%RL2_T2_to_Pa * pressure(i) ; enddo + spv_reference = EOS%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 + do i=is,ie ; pres(i) = EOS%RL2_T2_to_Pa * pressure(i) ; enddo call calculate_spec_vol_array(T, S, pres, specvol, start, 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) + if (EOS%R_to_kg_m3 /= 1.0) then ; do i=is,ie + specvol(i) = EOS%R_to_kg_m3 * specvol(i) enddo ; endif end subroutine calc_spec_vol_HI_1d @@ -578,7 +585,7 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -615,7 +622,7 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then ; do j=start,start+npts-1 drho_dT(j) = rho_scale * drho_dT(j) @@ -652,16 +659,16 @@ subroutine calculate_density_derivs_HI_1d(T, S, pressure, drho_dT, drho_dS, HI, 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 + if (EOS%RL2_T2_to_Pa == 1.0) then call calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS) else - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + do i=is,ie ; pres(i) = EOS%RL2_T2_to_Pa * pressure(i) ; enddo call calculate_density_derivs_array(T, S, pres, drho_dT, drho_dS, start, 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) + if (EOS%kg_m3_to_R /= 1.0) then ; do i=is,ie + drho_dT(i) = EOS%kg_m3_to_R * drho_dT(i) + drho_dS(i) = EOS%kg_m3_to_R * drho_dS(i) enddo ; endif end subroutine calculate_density_derivs_HI_1d @@ -690,7 +697,7 @@ subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -704,7 +711,7 @@ subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS call MOM_error(FATAL, "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.") end select - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then drho_dT = rho_scale * drho_dT @@ -746,7 +753,7 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -779,7 +786,7 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then ; do j=start,start+npts-1 drho_dS_dS(j) = rho_scale * drho_dS_dS(j) @@ -827,7 +834,7 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -843,7 +850,7 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr call MOM_error(FATAL, "calculate_density_derivs: EOS%form_of_EOS is not valid.") end select - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then drho_dS_dS = rho_scale * drho_dS_dS @@ -983,16 +990,16 @@ subroutine calc_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, U 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 + if (EOS%RL2_T2_to_Pa == 1.0) then call calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS) else - do i=is,ie ; press(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + do i=is,ie ; press(i) = EOS%RL2_T2_to_Pa * pressure(i) ; enddo call calculate_spec_vol_derivs_array(T, S, press, dSV_dT, dSV_dS, start, 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) + if (EOS%R_to_kg_m3 /= 1.0) then ; do i=is,ie + dSV_dT(i) = EOS%R_to_kg_m3 * dSV_dT(i) + dSV_dS(i) = EOS%R_to_kg_m3 * dSV_dS(i) enddo ; endif end subroutine calc_spec_vol_derivs_HI_1d @@ -1141,8 +1148,8 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) if (present(US)) then - call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, US%kg_m3_to_R*EOS%Rho_T0_S0, & - US%kg_m3_to_R*EOS%dRho_dT, US%kg_m3_to_R*EOS%dRho_dS, dza, & + call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, & + EOS%kg_m3_to_R*EOS%dRho_dT, EOS%kg_m3_to_R*EOS%dRho_dS, dza, & intp_dza, intx_dza, inty_dza, halo_size, & bathyP, dP_tiny, useMassWghtInterp) else @@ -1155,7 +1162,7 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & if (present(US)) then call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa) else call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp) @@ -1227,7 +1234,7 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R if (rho_scale /= 1.0) then call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & rho_scale*EOS%Rho_T0_S0, rho_scale*EOS%dRho_dT, rho_scale*EOS%dRho_dS, & @@ -1238,8 +1245,8 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) endif case (EOS_WRIGHT) - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R - pres_scale = 1.0 ; if (present(US)) pres_scale = US%RL2_T2_to_Pa + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + pres_scale = 1.0 ; if (present(US)) pres_scale = EOS%RL2_T2_to_Pa if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & @@ -1267,9 +1274,11 @@ logical function query_compressible(EOS) end function query_compressible !> Initializes EOS_type by allocating and reading parameters -subroutine EOS_init(param_file, EOS) +subroutine EOS_init(param_file, EOS, US) type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + optional :: US ! Local variables #include "version_variable.h" character(len=40) :: mdl = "MOM_EOS" ! This module's name. @@ -1363,6 +1372,11 @@ subroutine EOS_init(param_file, EOS) "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .") endif + ! Unit conversions + EOS%m_to_Z = 1. ; if (present(US)) EOS%m_to_Z = US%m_to_Z + EOS%kg_m3_to_R = 1. ; if (present(US)) EOS%kg_m3_to_R = US%kg_m3_to_R + EOS%R_to_kg_m3 = 1. ; if (present(US)) EOS%R_to_kg_m3 = US%R_to_kg_m3 + EOS%RL2_T2_to_Pa = 1. ; if (present(US)) EOS%RL2_T2_to_Pa = US%RL2_T2_to_Pa end subroutine EOS_init @@ -1521,9 +1535,9 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 do_massWeight = .false. @@ -1748,9 +1762,9 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & Isq = HIO%IscB ; Ieq = HIO%IecB ; Jsq = HIO%JscB ; Jeq = HIO%JecB - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 massWeightToggle = 0. if (present(useMassWghtInterp)) then @@ -2020,7 +2034,7 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t Pa_left = P_t - P_tgt ! Pa_left < 0 F_r = 1. Pa_right = P_b - P_tgt ! Pa_right > 0 - Pa_tol = GxRho * 1.0e-5*US%m_to_Z + Pa_tol = GxRho * 1.0e-5*EOS%m_to_Z if (present(z_tol)) Pa_tol = GxRho * z_tol F_guess = F_l - Pa_left / (Pa_right - Pa_left) * (F_r - F_l) @@ -2197,9 +2211,9 @@ subroutine int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, & is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 + rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 ! ============================= @@ -2639,9 +2653,9 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - SV_scale = 1.0 ; if (present(US)) SV_scale = US%R_to_kg_m3 - RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = US%RL2_T2_to_Pa - alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * US%kg_m3_to_R + SV_scale = 1.0 ; if (present(US)) SV_scale = EOS%R_to_kg_m3 + RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = EOS%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -2863,9 +2877,9 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, do_massWeight = .false. if (present(useMassWghtInterp)) do_massWeight = useMassWghtInterp - SV_scale = 1.0 ; if (present(US)) SV_scale = US%R_to_kg_m3 - RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = US%RL2_T2_to_Pa - alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * US%kg_m3_to_R + SV_scale = 1.0 ; if (present(US)) SV_scale = EOS%R_to_kg_m3 + RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = EOS%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) From 18d520e3c0da50955c9f0276c708c67bdb97952e Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Apr 2020 17:16:49 +0000 Subject: [PATCH 2/4] Remove optional US from most MOM_EOS functions - Removes US as an argument wherever it was optional since the unit conversion factors are not stored in the EOS type. --- src/ALE/MOM_regridding.F90 | 6 +- src/ALE/coord_hycom.F90 | 6 +- src/ALE/coord_rho.F90 | 11 +- src/ALE/coord_slight.F90 | 12 +- src/core/MOM.F90 | 2 +- src/core/MOM_PressureForce_Montgomery.F90 | 20 +- src/core/MOM_PressureForce_analytic_FV.F90 | 22 +-- src/core/MOM_PressureForce_blocked_AFV.F90 | 18 +- src/core/MOM_forcing_type.F90 | 2 +- src/core/MOM_interface_heights.F90 | 4 +- src/core/MOM_isopycnal_slopes.F90 | 4 +- src/diagnostics/MOM_diagnostics.F90 | 12 +- src/diagnostics/MOM_wave_speed.F90 | 4 +- src/diagnostics/MOM_wave_structure.F90 | 2 +- src/equation_of_state/MOM_EOS.F90 | 178 +++++++----------- src/framework/MOM_diag_remap.F90 | 2 +- src/ice_shelf/MOM_ice_shelf.F90 | 4 +- .../MOM_coord_initialization.F90 | 6 +- .../MOM_state_initialization.F90 | 24 +-- .../lateral/MOM_mixed_layer_restrat.F90 | 8 +- .../lateral/MOM_thickness_diffuse.F90 | 8 +- .../vertical/MOM_bulk_mixed_layer.F90 | 8 +- .../vertical/MOM_diabatic_aux.F90 | 16 +- .../vertical/MOM_diabatic_driver.F90 | 2 +- .../vertical/MOM_diapyc_energy_req.F90 | 6 +- .../vertical/MOM_entrain_diffusive.F90 | 8 +- .../vertical/MOM_full_convection.F90 | 6 +- .../vertical/MOM_geothermal.F90 | 6 +- .../vertical/MOM_internal_tide_input.F90 | 2 +- .../vertical/MOM_kappa_shear.F90 | 2 +- .../vertical/MOM_regularize_layers.F90 | 4 +- .../vertical/MOM_set_diffusivity.F90 | 10 +- .../vertical/MOM_set_viscosity.F90 | 12 +- src/tracer/MOM_neutral_diffusion.F90 | 20 +- src/tracer/MOM_tracer_Z_init.F90 | 8 +- src/tracer/MOM_tracer_hor_diff.F90 | 2 +- src/user/DOME_initialization.F90 | 8 +- src/user/ISOMIP_initialization.F90 | 12 +- src/user/RGC_initialization.F90 | 2 +- src/user/benchmark_initialization.F90 | 16 +- src/user/user_change_diffusivity.F90 | 4 +- 41 files changed, 231 insertions(+), 278 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 1000ba0d32..5ef65342e5 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1386,7 +1386,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) ! Local depth (G%bathyT is positive) nominalDepth = G%bathyT(i,j)*GV%Z_to_H - call build_rho_column(CS%rho_CS, US, nz, nominalDepth, h(i, j, :), & + call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i, j, :), & tv%T(i, j, :), tv%S(i, j, :), tv%eqn_of_state, zNew, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) @@ -1501,7 +1501,7 @@ subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS ) ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) enddo - call build_hycom1_column(CS%hycom_CS, US, tv%eqn_of_state, GV%ke, depth, & + call build_hycom1_column(CS%hycom_CS, tv%eqn_of_state, GV%ke, depth, & h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, & z_col, z_col_new, zScale=GV%Z_to_H, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) @@ -1635,7 +1635,7 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) ( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref ) enddo - call build_slight_column(CS%slight_CS, US, tv%eqn_of_state, GV%H_to_RZ*GV%g_Earth, & + call build_slight_column(CS%slight_CS, tv%eqn_of_state, GV%H_to_RZ*GV%g_Earth, & GV%H_subroundoff, nz, depth, h(i, j, :), & tv%T(i, j, :), tv%S(i, j, :), p_col, z_col, z_col_new, & h_neglect=h_neglect, h_neglect_edge=h_neglect_edge) diff --git a/src/ALE/coord_hycom.F90 b/src/ALE/coord_hycom.F90 index bfcff9005c..064860301d 100644 --- a/src/ALE/coord_hycom.F90 +++ b/src/ALE/coord_hycom.F90 @@ -4,7 +4,6 @@ module coord_hycom ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, FATAL -use MOM_unit_scaling, only : unit_scale_type use MOM_EOS, only : EOS_type, calculate_density use regrid_interp, only : interp_CS_type, build_and_interpolate_grid @@ -96,10 +95,9 @@ subroutine set_hycom_params(CS, max_interface_depths, max_layer_thickness, inter end subroutine set_hycom_params !> Build a HyCOM coordinate column -subroutine build_hycom1_column(CS, US, eqn_of_state, nz, depth, h, T, S, p_col, & +subroutine build_hycom1_column(CS, eqn_of_state, nz, depth, h, T, S, p_col, & z_col, z_col_new, zScale, h_neglect, h_neglect_edge) type(hycom_CS), intent(in) :: CS !< Coordinate control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(EOS_type), pointer :: eqn_of_state !< Equation of state structure integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom (positive [H ~> m or kg m-2]) @@ -133,7 +131,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, 1, nz, eqn_of_state) ! 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..0da2a33554 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -5,7 +5,6 @@ module coord_rho use MOM_error_handler, only : MOM_error, FATAL use MOM_remapping, only : remapping_CS, remapping_core_h -use MOM_unit_scaling, only : unit_scale_type use MOM_EOS, only : EOS_type, calculate_density use regrid_interp, only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX @@ -88,10 +87,9 @@ end subroutine set_rho_params !! !! 1. Density profiles are calculated on the source grid. !! 2. Positions of target densities (for interfaces) are found by interpolation. -subroutine build_rho_column(CS, US, nz, depth, h, T, S, eqn_of_state, z_interface, & +subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & h_neglect, h_neglect_edge) type(rho_CS), intent(in) :: CS !< coord_rho control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) real, intent(in) :: depth !< Depth of ocean bottom (positive in m) real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] @@ -126,7 +124,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, 1, nz, eqn_of_state) do k = 1,count_nonzero_layers densities(k) = densities(mapping(k)) enddo @@ -185,10 +183,9 @@ end subroutine build_rho_column !! 4. T & S are remapped onto the new grid. !! 5. Return to step 1 until convergence or until the maximum number of !! iterations is reached, whichever comes first. -subroutine build_rho_column_iteratively(CS, US, remapCS, nz, depth, h, T, S, eqn_of_state, & +subroutine build_rho_column_iteratively(CS, remapCS, nz, depth, h, T, S, eqn_of_state, & zInterface, h_neglect, h_neglect_edge, dev_tol) type(rho_CS), intent(in) :: CS !< Regridding control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options integer, intent(in) :: nz !< Number of levels real, intent(in) :: depth !< Depth of ocean bottom [Z ~> m] @@ -250,7 +247,7 @@ subroutine build_rho_column_iteratively(CS, US, remapCS, nz, depth, h, T, S, eqn enddo ! Compute densities within current water column - call calculate_density( T_tmp, S_tmp, pres, densities, 1, nz, eqn_of_state, US=US) + call calculate_density( T_tmp, S_tmp, pres, densities, 1, nz, eqn_of_state) do k = 1,count_nonzero_layers densities(k) = densities(mapping(k)) diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90 index 000315bae8..409b78c37c 100644 --- a/src/ALE/coord_slight.F90 +++ b/src/ALE/coord_slight.F90 @@ -4,7 +4,6 @@ module coord_slight ! This file is part of MOM6. See LICENSE.md for the license. use MOM_error_handler, only : MOM_error, FATAL -use MOM_unit_scaling, only : unit_scale_type use MOM_EOS, only : EOS_type, calculate_compress use MOM_EOS, only : calculate_density, calculate_density_derivs use regrid_interp, only : interp_CS_type, regridding_set_ppolys @@ -178,11 +177,10 @@ subroutine set_slight_params(CS, max_interface_depths, max_layer_thickness, & end subroutine set_slight_params !> Build a SLight coordinate column -subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, & +subroutine build_slight_column(CS, eqn_of_state, H_to_pres, H_subroundoff, & nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, & h_neglect, h_neglect_edge) type(slight_CS), intent(in) :: CS !< Coordinate control structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(EOS_type), pointer :: eqn_of_state !< Equation of state structure real, intent(in) :: H_to_pres !< A conversion factor from thicknesses to !! scaled pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] @@ -253,7 +251,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, 1, nz, eqn_of_state) ! Find the locations of the target potential densities, flagging ! locations in apparently unstable regions as not reliable. @@ -375,11 +373,11 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, & 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) + eqn_of_state) call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, & - eqn_of_state, US) + eqn_of_state) 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) + call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state) else do K=2,nz ; drho_dp(K) = 0.0 ; enddo endif diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f07fd6a1c4..f25b8792f9 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2905,7 +2905,7 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) 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) + tv%eqn_of_state) 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_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 0d8cf27dad..c8e94ca7d8 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -187,7 +187,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb !$OMP parallel do default(shared) do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), & - 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, US=US) + 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1) enddo !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 @@ -229,7 +229,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -246,7 +246,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb !$OMP parallel do default(shared) private(rho_in_situ) do k=1,nz ; do j=Jsq,Jeq+1 call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo enddo ; enddo endif ! use_EOS @@ -485,7 +485,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -506,7 +506,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, !$OMP parallel do default(shared) do k=1,nz+1 ; do j=Jsq,Jeq+1 call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo enddo ; enddo endif ! use_EOS @@ -663,7 +663,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) press(i) = -Rho0xG*e(i,j,1) enddo call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z enddo @@ -674,7 +674,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k)) enddo call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * & ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * & @@ -760,7 +760,7 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ) do j=Jsq,Jeq+1 call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) pbce(i,j,nz) = dP_dH / (rho_in_situ(i)) @@ -770,9 +770,9 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1)) S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1)) enddo - call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, start, npts, tv%eqn_of_state, US=US) + call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, start, npts, tv%eqn_of_state) call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * & ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + & diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index f0a4485399..bb5af350cb 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -229,7 +229,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -266,21 +266,21 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, & tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), & - useMassWghtInterp=CS%useMassWghtInterp, US=US) + useMassWghtInterp=CS%useMassWghtInterp) elseif ( CS%Recon_Scheme == 2 ) then call MOM_error(FATAL, "PressureForce_AFV_nonBouss: "//& "int_spec_vol_dp_generic_ppm does not exist yet.") ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), & ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), & - ! intx_dza(:,:,k), inty_dza(:,:,k), US=US) + ! intx_dza(:,:,k), inty_dza(:,:,k)) endif else call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), & p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & - useMassWghtInterp=CS%useMassWghtInterp, US=US) + useMassWghtInterp=CS%useMassWghtInterp) endif else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref @@ -335,7 +335,7 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p !$OMP parallel do default(shared) private(rho_in_situ) do j=Jsq,Jeq+1 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j)) @@ -579,7 +579,7 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -602,10 +602,10 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at do j=Jsq,Jeq+1 if (use_p_atm) then call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) else call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) endif do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) @@ -670,17 +670,17 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k),& e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, & - useMassWghtInterp=CS%useMassWghtInterp, US=US) + useMassWghtInterp=CS%useMassWghtInterp) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, & - intz_dpa, intx_dpa, inty_dpa, US=US) + intz_dpa, intx_dpa, inty_dpa) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, & - intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US) + intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp) endif !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90 index 5ac7831479..9c04ad9684 100644 --- a/src/core/MOM_PressureForce_blocked_AFV.F90 +++ b/src/core/MOM_PressureForce_blocked_AFV.F90 @@ -225,7 +225,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb) @@ -249,7 +249,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & - useMassWghtInterp=CS%useMassWghtInterp, US=US) + useMassWghtInterp=CS%useMassWghtInterp) else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -303,7 +303,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, !$OMP parallel do default(shared) private(rho_in_situ) do j=Jsq,Jeq+1 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j)) @@ -566,7 +566,7 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k) enddo ; enddo call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) do k=nkmb+1,nz ; do i=Isq,Ieq+1 if (GV%Rlay(k) < Rho_cv_BL(i)) then @@ -589,10 +589,10 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, do j=Jsq,Jeq+1 if (use_p_atm) then call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) else call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) endif do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) @@ -670,18 +670,18 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, & G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - useMassWghtInterp=CS%useMassWghtInterp, US=US) + useMassWghtInterp=CS%useMassWghtInterp) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, & GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & - intx_dpa_bk, inty_dpa_bk, US=US) + intx_dpa_bk, inty_dpa_bk) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, & dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US) + G%bathyT, dz_neglect, CS%useMassWghtInterp) endif do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*GV%Z_to_H diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 12f165372b..d5a1fe3c79 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -954,7 +954,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt ! Density derivatives call calculate_density_derivs(Temp(:,j,1), Salt(:,j,1), pressure, dRhodT, dRhodS, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) ! 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/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index bfb9ad2703..ea529d42c5 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -106,7 +106,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !$OMP do do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & - 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US) + 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo) enddo !$OMP do do j=jsv,jev @@ -208,7 +208,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !$OMP do do k = 1, nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, & - G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US) + G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo) enddo !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 78fdc51077..6cb7e049a6 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -176,7 +176,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1))) enddo call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, (is-IsdB+1)-1, ie-is+2, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) endif do I=is-1,ie @@ -262,7 +262,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & 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, is, ie-is+1, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) endif do i=is,ie if (use_EOS) then diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 94f6acc9c3..d653ddec6c 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -360,7 +360,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & 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) + tv%eqn_of_state) 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 @@ -466,7 +466,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & !$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) + tv%eqn_of_state, 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,7 @@ 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), G%HI, tv%eqn_of_state) enddo ; enddo if (CS%id_rhopot0 > 0) call post_data(CS%id_rhopot0, Rcv, CS%diag) endif @@ -596,7 +596,7 @@ 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), G%HI, tv%eqn_of_state) enddo ; enddo if (CS%id_rhopot2 > 0) call post_data(CS%id_rhopot2, Rcv, CS%diag) endif @@ -606,7 +606,7 @@ 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), G%HI, tv%eqn_of_state) pressure_1d(:) = pressure_1d(:) + 0.5 * h(:,j,k) * (GV%H_to_RZ*GV%g_Earth) ! Pressure at bottom of layer k enddo enddo @@ -839,7 +839,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k) enddo ; enddo call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, & - G%HI, G%HI, tv%eqn_of_state, dpress, US=US) + G%HI, G%HI, tv%eqn_of_state, dpress) do j=js,je ; do i=is,ie mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth enddo ; enddo diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 85dbcdc13b..ce36835c1a 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -243,7 +243,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, & S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + kf(i)-1, tv%eqn_of_state) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. @@ -738,7 +738,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + kf(i)-1, tv%eqn_of_state) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index ceb6fd6c4f..69c5bcb44f 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -278,7 +278,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo S_int(k) = 0.5*(Sf(k,i)+Sf(k-1,i)) enddo call calculate_density_derivs(T_int, S_int, pres, drho_dT, drho_dS, 2, & - kf(i)-1, tv%eqn_of_state, US) + kf(i)-1, tv%eqn_of_state) ! Sum the reduced gravities to find out how small a density difference ! is negligibly small. diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index dc74e1dcf7..11763e066a 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -114,10 +114,11 @@ module MOM_EOS ! Unit conversion factors (normally used for dimensional testing but could also allow for ! change of units of arguments to functions) - real :: m_to_Z !< A constant that translates distances in meters to the units of depth. - real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density. - real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed. - real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa. + real :: m_to_Z !< A constant that translates distances in meters to the units of depth. + real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density. + real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed. + real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa. + real :: L_T_to_m_s !< Convert lateral velocities from L T-1 to m s-1. ! logical :: test_EOS = .true. ! If true, test the equation of state end type EOS_type @@ -151,14 +152,13 @@ module MOM_EOS !! If rho_ref is present, the anomaly with respect to rho_ref is returned. The pressure and !! density can be rescaled with the US. If both the US and scale arguments are present the density !! scaling uses the product of the two scaling factors. -subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, US, scale) +subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, scale) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] real, intent(out) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] type(EOS_type), pointer :: EOS !< Equation of state structure real, optional, intent(in) :: rho_ref !< A reference density [kg m-3] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density in !! combination with scaling given by US [various] @@ -168,7 +168,7 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, US, scale if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_scalar called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -186,17 +186,15 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS, rho_ref, US, scale call MOM_error(FATAL, "calculate_density_scalar: EOS is not valid.") end select - if (present(US) .or. present(scale)) then - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R - if (present(scale)) rho_scale = rho_scale * scale - rho = rho_scale * rho - endif + rho_scale = EOS%kg_m3_to_R + if (present(scale)) rho_scale = rho_scale * scale + rho = rho_scale * rho end subroutine calculate_density_scalar !> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs. !! If rho_ref is present, the anomaly with respect to rho_ref is returned. -subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref, US, scale) +subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, 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 [Pa] or [R L2 T-2 ~> Pa] @@ -205,7 +203,6 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re integer, intent(in) :: npts !< Number of point to compute type(EOS_type), pointer :: EOS !< Equation of state structure real, optional, intent(in) :: rho_ref !< A reference density [kg m-3] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density !! in combination with scaling given by US [various] @@ -217,7 +214,7 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_array called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -254,7 +251,7 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then do j=start,start+npts-1 ; rho(j) = rho_scale * rho(j) ; enddo @@ -265,14 +262,13 @@ 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. !! 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) +subroutine calculate_density_HI_1d(T, S, pressure, rho, HI, EOS, 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. ! Local variables @@ -350,14 +346,13 @@ end subroutine calculate_spec_vol_array !> Calls the appropriate subroutine to calculate specific volume of sea water !! for scalar inputs. -subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, US, scale) +subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] real, intent(out) :: specvol !< In situ? specific volume [m3 kg-1] or [R-1 ~> m3 kg-1] type(EOS_type), pointer :: EOS !< Equation of state structure real, optional, intent(in) :: spv_ref !< A reference specific volume [m3 kg-1] or [R-1 m3 kg-1] - type(unit_scale_type), optional, 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] @@ -368,18 +363,18 @@ subroutine calc_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, US, scale if (.not.associated(EOS)) call MOM_error(FATAL, & "calc_spec_vol_scalar called with an unassociated EOS_type EOS.") - pres(1) = pressure ; if (present(US)) pres(1) = EOS%RL2_T2_to_Pa*pressure + pres(1) = EOS%RL2_T2_to_Pa*pressure Ta(1) = T ; Sa(1) = S if (present(spv_ref)) then - spv_reference = spv_ref ; if (present(US)) spv_reference = EOS%kg_m3_to_R*spv_ref + spv_reference = EOS%kg_m3_to_R*spv_ref call calculate_spec_vol_array(Ta, Sa, pres, spv, 1, 1, EOS, spv_reference) else call calculate_spec_vol_array(Ta, Sa, pres, spv, 1, 1, EOS) endif specvol = spv(1) - spv_scale = 1.0 ; if (present(US)) spv_scale = EOS%R_to_kg_m3 + spv_scale = EOS%R_to_kg_m3 if (present(scale)) spv_scale = spv_scale * scale if (spv_scale /= 1.0) then specvol = spv_scale * specvol @@ -436,14 +431,13 @@ 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) +subroutine calc_spec_vol_HI_1d(T, S, pressure, specvol, HI, EOS, 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] @@ -561,7 +555,7 @@ subroutine calculate_TFreeze_array(S, pressure, T_fr, start, npts, EOS, pres_sca end subroutine calculate_TFreeze_array !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. -subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, US, scale) +subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS, 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 [Pa] or [R L2 T-2 ~> Pa] @@ -572,7 +566,6 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star 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), optional, intent(in) :: US !< A dimensional unit scaling type real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density !! in combination with scaling given by US [various] @@ -585,7 +578,7 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -622,7 +615,7 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then ; do j=start,start+npts-1 drho_dT(j) = rho_scale * drho_dT(j) @@ -633,7 +626,7 @@ 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) +subroutine calculate_density_derivs_HI_1d(T, S, pressure, drho_dT, drho_dS, HI, EOS, 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] @@ -643,7 +636,6 @@ subroutine calculate_density_derivs_HI_1d(T, S, pressure, drho_dT, drho_dS, HI, 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. ! Local variables @@ -676,7 +668,7 @@ end subroutine calculate_density_derivs_HI_1d !> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar !! to a one-element array -subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, US, scale) +subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, scale) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] @@ -685,10 +677,8 @@ subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS real, intent(out) :: drho_dS !< The partial derivative of density with salinity, !! in [kg m-3 ppt-1] or [R ppt-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 real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density !! in combination with scaling given by US [various] - ! Local variables 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] @@ -697,7 +687,7 @@ subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -711,7 +701,7 @@ subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS call MOM_error(FATAL, "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.") end select - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then drho_dT = rho_scale * drho_dT @@ -722,7 +712,7 @@ end subroutine calculate_density_derivs_scalar !> Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs. subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, & - drho_dS_dP, drho_dT_dP, start, npts, EOS, US, scale) + drho_dS_dP, drho_dT_dP, start, npts, EOS, 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 [Pa] or [R L2 T-2 ~> Pa] @@ -739,10 +729,8 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh 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), optional, intent(in) :: US !< A dimensional unit scaling type 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(size(pressure)) :: 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] @@ -753,7 +741,7 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then select case (EOS%form_of_EOS) @@ -786,7 +774,7 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh end select endif - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then ; do j=start,start+npts-1 drho_dS_dS(j) = rho_scale * drho_dS_dS(j) @@ -808,7 +796,7 @@ end subroutine calculate_density_second_derivs_array !> Calls the appropriate subroutine to calculate density second derivatives for scalar nputs. subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, & - drho_dS_dP, drho_dT_dP, EOS, US, scale) + drho_dS_dP, drho_dT_dP, EOS, scale) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] @@ -823,7 +811,6 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure !! [kg m-3 degC-1 Pa-1] or [R degC-1 Pa-1 ~> kg m-3 degC-1 Pa-1] type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type real, optional, intent(in) :: scale !< A multiplicative factor by which to scale density !! in combination with scaling given by US [various] ! Local variables @@ -834,7 +821,7 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_density_derivs called with an unassociated EOS_type EOS.") - p_scale = 1.0 ; if (present(US)) p_scale = EOS%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -850,7 +837,7 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr call MOM_error(FATAL, "calculate_density_derivs: EOS%form_of_EOS is not valid.") end select - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (present(scale)) rho_scale = rho_scale * scale if (rho_scale /= 1.0) then drho_dS_dS = rho_scale * drho_dS_dS @@ -964,7 +951,7 @@ 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) +subroutine calc_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, 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] @@ -974,7 +961,6 @@ subroutine calc_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, U 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 @@ -1007,7 +993,7 @@ end subroutine calc_spec_vol_derivs_HI_1d !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array !! inputs. If US is present, the units of the inputs and outputs are rescaled. -subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS, US) +subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS) real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [PSU] real, dimension(:), intent(in) :: press !< Pressure [Pa] or [R L2 T-2 ~> Pa] @@ -1018,7 +1004,6 @@ subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS, 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), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables real, dimension(size(press)) :: pressure ! Pressure converted to [Pa] @@ -1028,11 +1013,7 @@ subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS, "calculate_compress called with an unassociated EOS_type EOS.") is = start ; ie = is + npts - 1 - if (present(US)) then - do i=is,ie ; pressure(i) = US%RL2_T2_to_Pa * press(i) ; enddo - else - do i=is,ie ; pressure(i) = press(i) ; enddo - endif + do i=is,ie ; pressure(i) = EOS%RL2_T2_to_Pa * press(i) ; enddo select case (EOS%form_of_EOS) case (EOS_LINEAR) @@ -1050,21 +1031,19 @@ subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS, call MOM_error(FATAL, "calculate_compress: EOS%form_of_EOS is not valid.") end select - if (present(US)) then - if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie - rho(i) = US%kg_m3_to_R * rho(i) - enddo ; endif - if (US%L_T_to_m_s /= 1.0) then ; do i=is,ie - drho_dp(i) = US%L_T_to_m_s**2 * drho_dp(i) - enddo ; endif - endif + if (EOS%kg_m3_to_R /= 1.0) then ; do i=is,ie + rho(i) = EOS%kg_m3_to_R * rho(i) + enddo ; endif + if (EOS%L_T_to_m_s /= 1.0) then ; do i=is,ie + drho_dp(i) = EOS%L_T_to_m_s**2 * drho_dp(i) + enddo ; endif end subroutine calculate_compress_array !> Calculate density and compressibility for a scalar. This just promotes the scalar to an array !! with a singleton dimension and calls calculate_compress_array. If US is present, the units of !! the inputs and outputs are rescaled. -subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS, US) +subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] @@ -1072,7 +1051,6 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS, US) real, intent(out) :: drho_dp !< The partial derivative of density with pressure (also the !! inverse of the square of sound speed) [s2 m-2] or [T2 L-2] type(EOS_type), pointer :: EOS !< Equation of state structure - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa @@ -1080,7 +1058,7 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS, US) "calculate_compress called with an unassociated EOS_type EOS.") Ta(1) = T ; Sa(1) = S; pa(1) = pressure - call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS, US) + call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS) rho = rhoa(1) ; drho_dp = drho_dpa(1) end subroutine calculate_compress_scalar @@ -1094,7 +1072,7 @@ end subroutine calculate_compress_scalar !! series for log(1-eps/1+eps) that assumes that |eps| < . 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) + bathyP, dP_tiny, useMassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -1131,8 +1109,6 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type - ! Local variables real :: pres_scale ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] real :: SV_scale ! A multiplicative factor by which to scale specific @@ -1144,33 +1120,21 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & if (EOS%EOS_quadrature) then call int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, US) + bathyP, dP_tiny, useMassWghtInterp) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - if (present(US)) then - call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, & - EOS%kg_m3_to_R*EOS%dRho_dT, EOS%kg_m3_to_R*EOS%dRho_dS, dza, & - intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) - else - call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%Rho_T0_S0, & - EOS%dRho_dT, EOS%dRho_dS, dza, intp_dza, & - intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp) - endif + call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, & + EOS%kg_m3_to_R*EOS%dRho_dT, EOS%kg_m3_to_R*EOS%dRho_dS, dza, & + intp_dza, intx_dza, inty_dza, halo_size, & + bathyP, dP_tiny, useMassWghtInterp) case (EOS_WRIGHT) - if (present(US)) then - call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & - SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa) - else - call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp) - endif + call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & + inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & + SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa) case default call int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, US) + bathyP, dP_tiny, useMassWghtInterp) end select ; endif end subroutine int_specific_vol_dp @@ -1179,7 +1143,7 @@ end subroutine int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1219,8 +1183,6 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type - ! Local variables real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the ! desired units [R m3 kg-1 ~> 1] @@ -1231,10 +1193,10 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp if (EOS%EOS_quadrature) then call int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R + rho_scale = EOS%kg_m3_to_R if (rho_scale /= 1.0) then call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & rho_scale*EOS%Rho_T0_S0, rho_scale*EOS%dRho_dT, rho_scale*EOS%dRho_dS, & @@ -1245,8 +1207,8 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) endif case (EOS_WRIGHT) - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R - pres_scale = 1.0 ; if (present(US)) pres_scale = EOS%RL2_T2_to_Pa + rho_scale = EOS%kg_m3_to_R + pres_scale = EOS%RL2_T2_to_Pa if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & @@ -1258,7 +1220,7 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp endif case default call int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) end select ; endif end subroutine int_density_dz @@ -1377,6 +1339,7 @@ subroutine EOS_init(param_file, EOS, US) EOS%kg_m3_to_R = 1. ; if (present(US)) EOS%kg_m3_to_R = US%kg_m3_to_R EOS%R_to_kg_m3 = 1. ; if (present(US)) EOS%R_to_kg_m3 = US%R_to_kg_m3 EOS%RL2_T2_to_Pa = 1. ; if (present(US)) EOS%RL2_T2_to_Pa = US%RL2_T2_to_Pa + EOS%L_T_to_m_s = 1. ; if (present(US)) EOS%L_T_to_m_s = US%L_T_to_m_s end subroutine EOS_init @@ -1459,7 +1422,7 @@ end subroutine EOS_use_linear !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & EOS, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, US) + bathyT, dz_neglect, useMassWghtInterp) type(hor_index_type), intent(in) :: HII !< Horizontal index type for input variables. type(hor_index_type), intent(in) :: HIO !< Horizontal index type for output variables. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1499,8 +1462,6 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type - ! Local variables real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa] @@ -1535,9 +1496,9 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 + rho_scale = EOS%kg_m3_to_R + GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 do_massWeight = .false. @@ -1668,7 +1629,7 @@ end subroutine int_density_dz_generic !! T and S are linear profiles. subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, US) + intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1708,7 +1669,6 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & !! divided by the y grid spacing [R L2 T-2 ~> Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -1762,9 +1722,9 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & Isq = HIO%IscB ; Ieq = HIO%IecB ; Jsq = HIO%JscB ; Jeq = HIO%JecB - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 + rho_scale = EOS%kg_m3_to_R + GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 massWeightToggle = 0. if (present(useMassWghtInterp)) then @@ -2102,7 +2062,7 @@ real function frac_dp_at_pos(T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, G_e, pos, EO T5(n) = top_weight * T_t + bottom_weight * T_b p5(n) = ( top_weight * z_t + bottom_weight * z_b ) * ( G_e * rho_ref ) enddo - call calculate_density_array(T5, S5, p5, rho5, 1, 5, EOS, US=US) + call calculate_density_array(T5, S5, p5, rho5, 1, 5, EOS) rho5(:) = rho5(:) !- rho_ref ! Work with anomalies relative to rho_ref ! Use Bode's rule to estimate the average density diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 83a7ce207c..cadd74950a 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -327,7 +327,7 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state) GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then !### I think that the conversion factor in the 2nd line should be GV%Z_to_H - call build_rho_column(get_rho_CS(remap_cs%regrid_cs), US, G%ke, & + call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, & US%Z_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 858af4e1ea..ad133cc4ab 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -376,9 +376,9 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! 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) + CS%eqn_of_state) call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, dR0_dT, dR0_dS, G%HI, & - CS%eqn_of_state, US) + CS%eqn_of_state) 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_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index b0155ae603..a70e761fa8 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -240,7 +240,7 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state ! The uppermost layer's density is set here. Subsequent layers' ! ! densities are determined from this value and the g values. ! ! T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref - call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, US=US) + call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state) ! These statements set the layer densities. ! do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo @@ -291,7 +291,7 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_s ! These statements set the interface reduced gravities. ! g_prime(1) = g_fs do k=1,nz ; Pref(k) = P_Ref ; enddo - call calculate_density(T0, S0, Pref, Rlay, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, Pref, Rlay, 1, nz, eqn_of_state) do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo call callTree_leave(trim(mdl)//'()') @@ -371,7 +371,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_sta g_prime(1) = g_fs do k=1,nz ; Pref(k) = P_Ref ; enddo - call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state, US=US) + call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state) ! Extrapolate target densities for the variable density mixed and buffer layers. do k=k_light-1,1,-1 Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 2efceb5991..4ff22d202b 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -960,7 +960,7 @@ subroutine convert_thickness(h, G, GV, US, tv) 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) + tv%eqn_of_state) do i=is,ie p_bot(i,j) = p_top(i,j) + HR_to_pres * (h(i,j,k) * rho(i)) enddo @@ -968,10 +968,10 @@ subroutine convert_thickness(h, G, GV, US, tv) do itt=1,max_itt 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) + tv%eqn_of_state, dz_geo) 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) + tv%eqn_of_state) ! 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 @@ -1600,8 +1600,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P T0(k) = T_Ref enddo - call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, US=US) + call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state) if (fit_salin) then ! A first guess of the layers' temperatures. @@ -1610,8 +1610,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P enddo ! Refine the guesses for each layer. do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS(k)) enddo @@ -1622,8 +1622,8 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1) enddo do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo @@ -1869,7 +1869,7 @@ 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), G%HI, tv%eqn_of_state) enddo call set_up_sponge_ML_density(tmp_2d, G, CSp) @@ -2188,7 +2188,7 @@ 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), G%HI, eos) enddo ; enddo call pass_var(temp_z,G%Domain) @@ -2449,7 +2449,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) S(k) = 35. + (0. * I_z_scale)*z(k) S_b(k) = 35. - (0. * I_z_scale)*e(k+1) call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*US%m_to_Z*z(k), & - rho(k), tv%eqn_of_state, US=US) + rho(k), tv%eqn_of_state) P_tot = P_tot + GV%g_Earth * rho(k) * GV%H_to_Z*h(k) enddo diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 744a801391..37f9cf2684 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -208,7 +208,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var 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) + tv%eqn_of_state, halo=1) deltaRhoAtK(:) = 0. MLD_fast(:,j) = 0. do k = 2, nz @@ -217,7 +217,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var ! 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) + tv%eqn_of_state, 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,7 @@ 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(:), G%HI, tv%eqn_of_state, halo=1) line_is_empty = .true. do i=is-1,ie+1 if (htot_fast(i,j) < MLD_fast(i,j)) then @@ -646,7 +646,7 @@ 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(:), G%HI, tv%eqn_of_state, 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..39884968f3 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -779,7 +779,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1))) enddo call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, & - drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, US=US) + drho_dS_u, (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state) endif do I=is-1,ie @@ -1030,7 +1030,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV 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) + tv%eqn_of_state) endif do i=is,ie if (calc_derivatives) then @@ -1262,7 +1262,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV S_u(I) = 0.5*(S(i,j,1) + S(i+1,j,1)) enddo call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u, drho_dS_u, & - (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state, US=US) + (is-IsdB+1)-1, ie-is+2, tv%eqn_of_state) endif do I=is-1,ie uhD(I,j,1) = -uhtot(I,j) @@ -1292,7 +1292,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV 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) + tv%eqn_of_state) endif do i=is,ie vhD(i,J,1) = -vhtot(i,J) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index cc3a6e3f69..ff915222bf 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -466,12 +466,12 @@ 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, dR0_dT, dR0_dS, G%HI, tv%eqn_of_state) call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) 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), G%HI, tv%eqn_of_state) + call calculate_density(T(:,k), S(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state) 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..8eb702d45d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -446,7 +446,7 @@ 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), G%HI, tv%eqn_of_state) enddo ! First, try to find an interior layer where inserting all the salt @@ -457,7 +457,7 @@ subroutine insert_brine(h, tv, G, GV, US, fluxes, nkmb, CS, dt, id_brine_lay) if ((G%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then s_new = S(i,k) + salt(i) / (GV%H_to_RZ * h_2d(i,k)) t0 = T(i,k) - call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state, US=US) + call calculate_density(t0, s_new, tv%P_Ref, R_new, tv%eqn_of_state) if (R_new < 0.5*(Rcv(i,k)+Rcv(i,k+1)) .and. s_new0.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 +1006,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), G%HI, tv%eqn_of_state) do i=is,ie ; dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; enddo enddo pen_TKE_2d(:,:) = 0.0 @@ -1348,7 +1348,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! Density derivatives call calculate_density_derivs(T2d(:,1), tv%S(:,j,1), SurfPressure, dRhodT, dRhodS, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) ! 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..5301fb5603 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2683,7 +2683,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e !$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) + tv%eqn_of_state) 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..21eb272d70 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, 1, nz, tv%eqn_of_state) do k=1,nz dMass = GV%H_to_RZ * h_tr(k) @@ -939,7 +939,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & N2(1) = 0.0 ; N2(nz+1) = 0.0 do K=2,nz call calculate_density(0.5*(T0(k-1) + T0(k)), 0.5*(S0(k-1) + S0(k)), & - pres(K), rho_here, tv%eqn_of_state, US=US) + pres(K), rho_here, tv%eqn_of_state) N2(K) = ((US%L_to_Z**2*GV%g_Earth) * rho_here / (0.5*GV%H_to_Z*(h_tr(k-1) + h_tr(k)))) * & ( 0.5*(dSV_dT(k-1) + dSV_dT(k)) * (T0(k-1) - T0(k)) + & 0.5*(dSV_dS(k-1) + dSV_dS(k)) * (S0(k-1) - S0(k)) ) @@ -950,7 +950,7 @@ subroutine diapyc_energy_req_calc(h_in, T_in, S_in, Kd, energy_Kd, dt, tv, & N2(1) = 0.0 ; N2(nz+1) = 0.0 do K=2,nz call calculate_density(0.5*(Tf(k-1) + Tf(k)), 0.5*(Sf(k-1) + Sf(k)), & - pres(K), rho_here, tv%eqn_of_state, US=US) + pres(K), rho_here, tv%eqn_of_state) N2(K) = ((US%L_to_Z**2*GV%g_Earth) * rho_here / (0.5*GV%H_to_Z*(h_tr(k-1) + h_tr(k)))) * & ( 0.5*(dSV_dT(k-1) + dSV_dT(k)) * (Tf(k-1) - Tf(k)) + & 0.5*(dSV_dS(k-1) + dSV_dS(k)) * (Sf(k-1) - Sf(k)) ) diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 100e79aba2..a44dc5c744 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -700,7 +700,7 @@ 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, G%HI, tv%eqn_of_state) 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 +784,7 @@ 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, G%HI, tv%eqn_of_state) 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 @@ -850,7 +850,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & endif enddo call calculate_density_derivs(T_eos, S_eos, pressure, dRho_dT, dRho_dS, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) do i=is,ie if ((k>kmb) .and. (k 0) then call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_Ref(:), Rcv_BL(:), isj, iej-isj+1, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) else Rcv_BL(:) = -1.0 endif @@ -245,11 +245,11 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, US, CS, halo) Rcv = 0.0 ; dRcv_dT = 0.0 ! Is this OK? else call calculate_density(tv%T(i,j,k), tv%S(i,j,k), tv%P_Ref, & - Rcv, tv%eqn_of_state, US=US) + Rcv, tv%eqn_of_state) T2(1) = tv%T(i,j,k) ; S2(1) = tv%S(i,j,k) T2(2) = tv%T(i,j,k_tgt) ; S2(2) = tv%S(i,j,k_tgt) call calculate_density_derivs(T2(:), S2(:), p_Ref(:), dRcv_dT_, dRcv_dS_, 1, 2, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) dRcv_dT = 0.5*(dRcv_dT_(1) + dRcv_dT_(2)) endif diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 3ba5520117..d366cb93d8 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -210,7 +210,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) enddo call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:), dRho_dS(:), G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) do i=is,ie dRho_int(i,K) = max(dRho_dT(i)*(T_f(i,j,k) - T_f(i,j,k-1)) + & dRho_dS(i)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 5e11ecee60..73a170eae2 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -911,7 +911,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & Sal_int(K) = 0.5*(Sal(k-1) + Sal(k)) enddo call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, 2, nzc-1, & - tv%eqn_of_state, US=US, scale=-g_R0) + tv%eqn_of_state, scale=-g_R0) else do K=1,nzc+1 ; dbuoy_dT(K) = -g_R0 ; dbuoy_dS(K) = 0.0 ; enddo endif diff --git a/src/parameterizations/vertical/MOM_regularize_layers.F90 b/src/parameterizations/vertical/MOM_regularize_layers.F90 index 43ba5211f7..e037697353 100644 --- a/src/parameterizations/vertical/MOM_regularize_layers.F90 +++ b/src/parameterizations/vertical/MOM_regularize_layers.F90 @@ -312,7 +312,7 @@ subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS) do j=js,je ; if (do_j(j)) then ! call cpu_clock_begin(id_clock_EOS) -! 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_cv, dRcv_dT, dRcv_dS, G%HI, tv%eqn_of_state) ! call cpu_clock_end(id_clock_EOS) do k=1,nz ; do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ; enddo ; enddo @@ -444,7 +444,7 @@ subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS) if (det_any) then call cpu_clock_begin(id_clock_EOS) do k=1,nkmb - call calculate_density(T_2d(:,k), S_2d(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state, US) + call calculate_density(T_2d(:,k), S_2d(:,k), p_ref_cv, Rcv(:,k), G%HI, tv%eqn_of_state) enddo call cpu_clock_end(id_clock_EOS) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index a65be62b3a..b333f04a62 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -715,10 +715,10 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ; enddo do k=1,nz call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_0, rho_0(:,k), G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) enddo call calculate_density(tv%T(:,j,kmb), tv%S(:,j,kmb), p_ref, Rcv_kmb, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) kb_min = kmb+1 do i=is,ie @@ -907,7 +907,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, & Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) enddo call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:,K), dRho_dS(:,K), G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) do i=is,ie dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) @@ -1066,7 +1066,7 @@ subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, US, CS, Kd_T_dd, Kd_S_dd) Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) enddo call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT, dRho_dS, G%HI, & - tv%eqn_of_state, US) + tv%eqn_of_state) do i=is,ie alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) @@ -1819,7 +1819,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) eps = 0.1 do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo do k=1,kmb - 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), G%HI, tv%eqn_of_state) enddo do i=is,ie if (kb(i) <= nz-1) then diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 4cd6a64684..6a2ce60f96 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -317,7 +317,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) !$OMP parallel do default(shared) do k=1,nkmb ; do j=Jsq,Jeq+1 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rml(:,j,k), start, npts, & - tv%eqn_of_state, US=US) + tv%eqn_of_state) enddo ; enddo endif @@ -574,7 +574,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) press(i) = press(i) + (GV%H_to_RZ*GV%g_Earth) * h_vel(i,k) enddo ; enddo call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, & - is-G%IsdB+1, ie-is+1, tv%eqn_of_state, US=US) + is-G%IsdB+1, ie-is+1, tv%eqn_of_state) endif do i=is,ie ; if (do_i(i)) then @@ -1278,7 +1278,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri S_EOS(I) = (h(i,j,k2)*tv%S(i,j,k2) + h(i+1,j,k2)*tv%S(i+1,j,k2)) * I_2hlay enddo call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, & - Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state, US=US) + Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state) endif do I=Isq,Ieq ; if (do_i(I)) then @@ -1399,7 +1399,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri if (use_EOS) then call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, & - Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state, US=US) + Isq-G%IsdB+1, Ieq-Isq+1, tv%eqn_of_state) endif do I=Isq,Ieq ; if (do_i(I)) then @@ -1515,7 +1515,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri S_EOS(i) = (h(i,j,k2)*tv%S(i,j,k2) + h(i,j+1,k2)*tv%S(i,j+1,k2)) * I_2hlay enddo call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, & - is-G%IsdB+1, ie-is+1, tv%eqn_of_state, US=US) + is-G%IsdB+1, ie-is+1, tv%eqn_of_state) endif do i=is,ie ; if (do_i(i)) then @@ -1636,7 +1636,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetri if (use_EOS) then call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, & - is-G%IsdB+1, ie-is+1, tv%eqn_of_state, US=US) + is-G%IsdB+1, ie-is+1, tv%eqn_of_state) endif do i=is,ie ; if (do_i(i)) then diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 6a109b7cba..e193ac9023 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -8,7 +8,7 @@ module MOM_neutral_diffusion use MOM_domains, only : pass_var use MOM_diag_mediator, only : diag_ctrl, time_type use MOM_diag_mediator, only : post_data, register_diag_field -use MOM_EOS, only : EOS_type, EOS_manual_init, calculate_compress, calculate_density_derivs +use MOM_EOS, only : EOS_type, EOS_manual_init, calculate_density_derivs use MOM_EOS, only : calculate_density, calculate_density_second_derivs use MOM_EOS, only : extract_member_EOS, EOS_LINEAR, EOS_TEOS10, EOS_WRIGHT use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe @@ -391,18 +391,18 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS) do k = 1, G%ke+1 if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k) call calculate_density_derivs(CS%Tint(:,j,k), CS%Sint(:,j,k), ref_pres, & - CS%dRdT(:,j,k), CS%dRdS(:,j,k), G%isc-1, G%iec-G%isc+3, CS%EOS, US) + CS%dRdT(:,j,k), CS%dRdS(:,j,k), G%isc-1, G%iec-G%isc+3, CS%EOS) enddo else ! Discontinuous reconstruction do k = 1, G%ke if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k) ! Calculate derivatives for the top interface call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, & - CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS, US) + CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, CS%EOS) if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k+1) ! Calcualte derivatives at the bottom interface call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, & - CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS, US) + CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, CS%EOS) enddo endif enddo @@ -1767,19 +1767,19 @@ subroutine calc_delta_rho_and_derivs(CS, US, T1, S1, p1_in, T2, S2, p2_in, drho, ! Use the full linear equation of state to calculate the difference in density (expensive!) if (TRIM(CS%delta_rho_form) == 'full') then pmid = 0.5 * (p1 + p2) - call calculate_density( T1, S1, pmid, rho1, CS%EOS, US=US ) - call calculate_density( T2, S2, pmid, rho2, CS%EOS, US=US ) + call calculate_density( T1, S1, pmid, rho1, CS%EOS) + call calculate_density( T2, S2, pmid, rho2, CS%EOS) drho = rho1 - rho2 ! Use the density derivatives at the average of pressures and the differentces int temperature elseif (TRIM(CS%delta_rho_form) == 'mid_pressure') then pmid = 0.5 * (p1 + p2) if (CS%ref_pres>=0) pmid = CS%ref_pres - call calculate_density_derivs(T1, S1, pmid, drdt1, drds1, CS%EOS, US) - call calculate_density_derivs(T2, S2, pmid, drdt2, drds2, CS%EOS, US) + call calculate_density_derivs(T1, S1, pmid, drdt1, drds1, CS%EOS) + call calculate_density_derivs(T2, S2, pmid, drdt2, drds2, CS%EOS) drho = delta_rho_from_derivs( T1, S1, p1, drdt1, drds1, T2, S2, p2, drdt2, drds2) elseif (TRIM(CS%delta_rho_form) == 'local_pressure') then - call calculate_density_derivs(T1, S1, p1, drdt1, drds1, CS%EOS, US) - call calculate_density_derivs(T2, S2, p2, drdt2, drds2, CS%EOS, US) + call calculate_density_derivs(T1, S1, p1, drdt1, drds1, CS%EOS) + call calculate_density_derivs(T2, S2, p2, drdt2, drds2, CS%EOS) drho = delta_rho_from_derivs( T1, S1, p1, drdt1, drds1, T2, S2, p2, drdt2, drds2) else call MOM_error(FATAL, "delta_rho_form is not recognized") diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 76ca2dac4a..c2748544c8 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -801,9 +801,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, adjust_salt = .true. iter_loop: do itt = 1,niter do k=1, nz - call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos) call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & - eos, US=US) + eos) enddo do k=k_start,nz ; do i=1,nx @@ -831,9 +831,9 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, if (adjust_salt .and. old_fit) then ; do itt = 1,niter do k=1, nz - call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, US=US) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos) call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & - eos, US=US) + eos) enddo do k=k_start,nz ; do i=1,nx ! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 6064589019..30a71951ba 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -701,7 +701,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, & !$OMP parallel do default(shared) do k=1,nkmb ; do j=js-2,je+2 call calculate_density(tv%T(:,j,k),tv%S(:,j,k), p_ref_cv, rho_coord(:,j,k), G%HI, & - tv%eqn_of_state, US, halo=2) + tv%eqn_of_state, halo=2) enddo ; enddo do j=js-2,je+2 ; do i=is-2,ie+2 diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 315e56051c..8330078555 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -359,13 +359,13 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, param_file, tr_Reg) ! target density and a salinity of 35 psu. This code is taken from ! USER_initialize_temp_sal. pres(:) = tv%P_Ref ; S0(:) = 35.0 ; T0(1) = 25.0 - call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), tv%eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, tv%eqn_of_state, US=US) + call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), tv%eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, tv%eqn_of_state) do k=1,nz ; T0(k) = T0(1) + (GV%Rlay(k)-rho_guess(1)) / drho_dT(1) ; enddo do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, tv%eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, tv%eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, tv%eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, tv%eqn_of_state) do k=1,nz ; T0(k) = T0(k) + (GV%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo enddo diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index c189cf0490..9f677f4e98 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -362,10 +362,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi ! call MOM_mesg(mesg,5) enddo - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, US) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state) ! write(mesg,*) 'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1) ! call MOM_mesg(mesg,5) - call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, US=US) + call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state) if (fit_salin) then ! A first guess of the layers' salinity. @@ -374,8 +374,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi enddo ! Refine the guesses for each layer. do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS1) enddo @@ -388,8 +388,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi enddo do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90 index 6dbee6cea7..61ccbf51ff 100644 --- a/src/user/RGC_initialization.F90 +++ b/src/user/RGC_initialization.F90 @@ -213,7 +213,7 @@ 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), G%HI, tv%eqn_of_state) enddo call set_up_sponge_ML_density(tmp, G, CSp) diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index 766474b364..492e51374c 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -152,8 +152,8 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state pres(k) = P_Ref ; S0(k) = 35.0 enddo T0(k1) = 29.0 - call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, US=US) + call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state) ! A first guess of the layers' temperatures. do k=1,nz @@ -162,8 +162,8 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state ! Refine the guesses for each layer. do itt=1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo @@ -257,8 +257,8 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & enddo T0(k1) = 29.0 - call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, US=US) + call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state) ! A first guess of the layers' temperatures. ! do k=1,nz @@ -267,8 +267,8 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & ! Refine the guesses for each layer. ! do itt = 1,6 - call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, US=US) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, US=US) + call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state) + call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state) do k=1,nz T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index f33d772352..834a624b7e 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -107,11 +107,11 @@ 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), G%HI, tv%eqn_of_state) 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), G%HI, tv%eqn_of_state) enddo endif From 0caf9ccabee4f0c9a351ddc6868c263eedf99c45 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Apr 2020 18:56:35 +0000 Subject: [PATCH 3/4] Missed a conflict resolution --- src/equation_of_state/MOM_EOS.F90 | 41 ++++++++++++++----------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index e74ace5de8..0623f955cf 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -607,7 +607,7 @@ end subroutine calculate_density_derivs_array !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. -subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, US, dom, scale) +subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, 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] @@ -616,7 +616,6 @@ subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, US 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 @@ -636,7 +635,7 @@ subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, US is = 1 ; ie = size(drho_dT) ; npts = 1 + ie - is endif - p_scale = 1.0 ; if (present(US)) p_scale = US%RL2_T2_to_Pa + p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then call calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, is, npts, EOS) @@ -645,7 +644,7 @@ subroutine calculate_density_derivs_1d(T, S, pressure, drho_dT, drho_dS, EOS, US call calculate_density_derivs_array(T, S, pres, drho_dT, drho_dS, is, npts, EOS) endif - rho_scale = US%kg_m3_to_R + rho_scale = EOS%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) @@ -897,7 +896,7 @@ end subroutine calculate_spec_vol_derivs_array !> 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) +subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, 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] @@ -906,7 +905,6 @@ subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, US, dom, 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 @@ -926,7 +924,7 @@ subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, US, dom, else 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 + p_scale = EOS%RL2_T2_to_Pa if (p_scale == 1.0) then call calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, is, npts, EOS) @@ -935,7 +933,7 @@ subroutine calc_spec_vol_derivs_1d(T, S, pressure, dSV_dT, dSV_dS, EOS, US, dom, call calculate_spec_vol_derivs_array(T, S, press, dSV_dT, dSV_dS, is, npts, EOS) endif - spv_scale = US%R_to_kg_m3 + spv_scale = EOS%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) @@ -2051,7 +2049,7 @@ end function frac_dp_at_pos !! are parabolic profiles subroutine int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, & z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & - EOS, dpa, intz_dpa, intx_dpa, inty_dpa, US) + EOS, dpa, intz_dpa, intx_dpa, inty_dpa) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays @@ -2091,7 +2089,6 @@ subroutine int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, & optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [R L2 T-2 ~> Pa] - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -2144,9 +2141,9 @@ subroutine int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, & is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - rho_scale = 1.0 ; if (present(US)) rho_scale = EOS%kg_m3_to_R - GxRho = G_e * rho_0 ; if (present(US)) GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * EOS%R_to_kg_m3 + rho_scale = EOS%kg_m3_to_R + GxRho = EOS%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref * EOS%R_to_kg_m3 I_Rho = 1.0 / rho_0 ! ============================= @@ -2509,7 +2506,7 @@ end subroutine evaluate_shape_quadratic !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, useMassWghtInterp, US) + bathyP, dP_neglect, useMassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature of the layer [degC] @@ -2547,7 +2544,6 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, & !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2586,9 +2582,9 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - SV_scale = 1.0 ; if (present(US)) SV_scale = EOS%R_to_kg_m3 - RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = EOS%RL2_T2_to_Pa - alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R + SV_scale = EOS%R_to_kg_m3 + RL2_T2_to_Pa = EOS%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -2724,7 +2720,7 @@ end subroutine int_spec_vol_dp_generic !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, & dP_neglect, bathyP, HI, EOS, dza, & - intp_dza, intx_dza, inty_dza, useMassWghtInterp, US) + intp_dza, intx_dza, inty_dza, useMassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T_t !< Potential temperature at the top of the layer [degC] @@ -2765,7 +2761,6 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2810,9 +2805,9 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, do_massWeight = .false. if (present(useMassWghtInterp)) do_massWeight = useMassWghtInterp - SV_scale = 1.0 ; if (present(US)) SV_scale = EOS%R_to_kg_m3 - RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = EOS%RL2_T2_to_Pa - alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R + SV_scale = EOS%R_to_kg_m3 + RL2_T2_to_Pa = EOS%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref * EOS%kg_m3_to_R do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) From 371ed85ca4abc1a6c865573cd11bf361393538cd Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 17 Apr 2020 19:29:12 +0000 Subject: [PATCH 4/4] Initialize scaling params without EOS_init() --- src/equation_of_state/MOM_EOS.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 0623f955cf..034912b9ff 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -115,11 +115,11 @@ module MOM_EOS ! Unit conversion factors (normally used for dimensional testing but could also allow for ! change of units of arguments to functions) - real :: m_to_Z !< A constant that translates distances in meters to the units of depth. - real :: kg_m3_to_R !< A constant that translates kilograms per meter cubed to the units of density. - real :: R_to_kg_m3 !< A constant that translates the units of density to kilograms per meter cubed. - real :: RL2_T2_to_Pa !< Convert pressures from R L2 T-2 to Pa. - real :: L_T_to_m_s !< Convert lateral velocities from L T-1 to m s-1. + real :: m_to_Z = 1. !< A constant that translates distances in meters to the units of depth. + real :: kg_m3_to_R = 1. !< A constant that translates kilograms per meter cubed to the units of density. + real :: R_to_kg_m3 = 1. !< A constant that translates the units of density to kilograms per meter cubed. + real :: RL2_T2_to_Pa = 1.!< Convert pressures from R L2 T-2 to Pa. + real :: L_T_to_m_s = 1. !< Convert lateral velocities from L T-1 to m s-1. ! logical :: test_EOS = .true. ! If true, test the equation of state end type EOS_type