From 7632abe2d78f74040dd7b75965ab4b4f8de92a50 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 5 Apr 2020 09:06:13 -0400 Subject: [PATCH] +Add optional SV_scale arg to int_specific_vol_dp Optionally rescale the units of the specific volume integrals. Added new optional SV_scale and pres_scale arguments to int_specific_vol_dp, int_spec_vol_dp_generic, and int_spec_vol_dp_Wright. All answers are bitwise identical, but there are new optional arguments to public interfaces. --- src/equation_of_state/MOM_EOS.F90 | 219 ++++++++++++++--------- src/equation_of_state/MOM_EOS_Wright.F90 | 70 +++++--- src/equation_of_state/MOM_EOS_linear.F90 | 50 +++--- 3 files changed, 202 insertions(+), 137 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index d1fc2a917b..5603246ace 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -800,43 +800,50 @@ 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) + bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) 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] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: S !< Salinity [ppt] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_t !< Pressure at the top of the layer [Pa]. + intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_b !< Pressure at the bottom of the layer [Pa]. + intent(in) :: p_b !< Pressure at the bottom of the layer [R L2 T-2 ~> Pa] or [Pa]. real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out - !! to reduce the magnitude of each of the integrals, m3 kg-1. The - !! calculation is mathematically identical with different values of + !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. + !! The calculation is mathematically identical with different values of !! alpha_ref, but this reduces the effects of roundoff. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(inout) :: dza !< The change in the geopotential anomaly across - !! the layer [m2 s-2]. + !! the layer [T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the !! geopotential anomaly relative to the anomaly at the bottom of the - !! layer [Pa m2 s-2]. + !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & optional, intent(inout) :: intx_dza !< The integral in x of the difference between the !! geopotential anomaly at the top and bottom of the layer divided by - !! the x grid spacing [m2 s-2]. + !! the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & optional, intent(inout) :: inty_dza !< The integral in y of the difference between the !! geopotential anomaly at the top and bottom of the layer divided by - !! the y grid spacing [m2 s-2]. + !! the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(in) :: bathyP !< The pressure at the bathymetry [Pa] + optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa] real, optional, intent(in) :: dP_tiny !< A miniscule pressure change with - !! the same units as p_t (Pa?) + !! 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. + real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific + !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure + !! into Pa [Pa T2 R-1 L-2 ~> 1]. + + ! Local variables + real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] if (.not.associated(EOS)) call MOM_error(FATAL, & "int_specific_vol_dp called with an unassociated EOS_type EOS.") @@ -844,21 +851,29 @@ 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) + bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%Rho_T0_S0, & + if (present(SV_scale)) then + rho_scale = 1.0 / SV_scale + call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, rho_scale*EOS%Rho_T0_S0, & + rho_scale*EOS%dRho_dT, rho_scale*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 case (EOS_WRIGHT) 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) + bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) 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) + bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) end select ; endif end subroutine int_specific_vol_dp @@ -1176,7 +1191,7 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure !! into Pa [Pa T2 R-1 L-2 ~> 1]. ! Local variables - real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] + 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] real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3] or [kg m-3] real :: rho_anom ! The depth averaged density anomaly [R ~> kg m-3] or [kg m-3]. @@ -2198,44 +2213,48 @@ 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) + bathyP, dP_neglect, useMassWghtInterp, SV_scale, pres_scale) 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]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: S !< Salinity of the layer [ppt]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_t !< Pressure atop the layer [Pa]. + intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_b !< Pressure below the layer [Pa]. - real, intent(in) :: alpha_ref !< A mean specific volume that is - !! subtracted out to reduce the magnitude of each of the - !! integrals [m3 kg-1]. The calculation is mathematically - !! identical with different values of alpha_ref, but alpha_ref - !! alters the effects of roundoff, and answers do change. + intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]. + real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out + !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. + !! The calculation is mathematically identical with different values of + !! alpha_ref, but alpha_ref alters the effects of roundoff, and + !! answers do change. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(inout) :: dza !< The change in the geopotential anomaly - !! across the layer [m2 s-2]. + !! across the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(inout) :: intp_dza !< The integral in pressure through the - !! layer of the geopotential anomaly relative to the anomaly - !! at the bottom of the layer [Pa m2 s-2]. + optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of + !! the geopotential anomaly relative to the anomaly at the bottom of the + !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(inout) :: intx_dza !< The integral in x of the difference - !! between the geopotential anomaly at the top and bottom of - !! the layer divided by the x grid spacing [m2 s-2]. + optional, intent(inout) :: intx_dza !< The integral in x of the difference between + !! the geopotential anomaly at the top and bottom of the layer divided + !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(inout) :: inty_dza !< The integral in y of the difference - !! between the geopotential anomaly at the top and bottom of - !! the layer divided by the y grid spacing [m2 s-2]. + optional, intent(inout) :: inty_dza !< The integral in y of the difference between + !! the geopotential anomaly at the top and bottom of the layer divided + !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(in) :: bathyP !< The pressure at the bathymetry [Pa] + optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with - !! the same units as p_t (Pa?) + !! 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. + real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific + !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure + !! into Pa [Pa T2 R-1 L-2 ~> 1]. ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2244,19 +2263,24 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & ! Bode's rule to do the horizontal integrals, and from a truncation in the ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. - real :: T5(5), S5(5), p5(5), a5(5) - real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1]. - real :: dp ! The pressure change through a layer [Pa]. -! real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa]. - real :: hWght ! A pressure-thickness below topography [Pa]. - real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Pa]. - real :: iDenom ! The inverse of the denominator in the weights [Pa-2]. + ! Local variables + real :: T5(5) ! Temperatures at five quadrature points [degC] + real :: S5(5) ! Salinities at five quadrature points [ppt] + real :: p5(5) ! Pressures at five quadrature points, scaled back to Pa if necessary [Pa] + real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: alpha_anom ! The depth averaged specific density anomaly [R-1 ~> m3 kg-1]. + real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]. + real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]. + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]. + real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1] + real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]. real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]. real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. real :: intp(5) ! The integrals of specific volume with pressure at the - ! 5 sub-column locations [m2 s-2]. + ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. + real :: RL2_T2_to_Pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_90 = 1.0/90.0 ! A rational constant. integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, halo @@ -2267,6 +2291,9 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & 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 + RL2_T2_to_Pa = 1.0 ; if (present(pres_scale)) RL2_T2_to_Pa = pres_scale + alpha_ref_mks = alpha_ref ; if (present(SV_scale)) alpha_ref_mks = alpha_ref / SV_scale + do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then do_massWeight = .true. @@ -2280,9 +2307,9 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dp = p_b(i,j) - p_t(i,j) do n=1,5 T5(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = p_b(i,j) - 0.25*real(n-1)*dp + p5(n) = RL2_T2_to_Pa * (p_b(i,j) - 0.25*real(n-1)*dp) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref) + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) ! Use Bode's rule to estimate the interface height anomaly change. alpha_anom = C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) @@ -2318,15 +2345,15 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & ! T, S, and p are interpolated in the horizontal. The p interpolation ! is linear, but for T and S it may be thickness wekghted. - p5(1) = wt_L*p_b(i,j) + wt_R*p_b(i+1,j) + p5(1) = RL2_T2_to_Pa * (wt_L*p_b(i,j) + wt_R*p_b(i+1,j)) dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - 0.25*dp + T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - RL2_T2_to_Pa * 0.25*dp enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref) + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) ! Use Bode's rule to estimate the interface height anomaly change. intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & @@ -2362,14 +2389,14 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & ! T, S, and p are interpolated in the horizontal. The p interpolation ! is linear, but for T and S it may be thickness wekghted. - p5(1) = wt_L*p_b(i,j) + wt_R*p_b(i,j+1) + p5(1) = RL2_T2_to_Pa * (wt_L*p_b(i,j) + wt_R*p_b(i,j+1)) dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) do n=2,5 - T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - 0.25*dp + T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = RL2_T2_to_Pa * (p5(n-1) - 0.25*dp) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref) + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) ! Use Bode's rule to estimate the interface height anomaly change. intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & @@ -2388,7 +2415,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) + intp_dza, intx_dza, inty_dza, useMassWghtInterp, SV_scale, pres_scale) 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]. @@ -2399,36 +2426,40 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: S_b !< Salinity at the bottom the layer [ppt]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_t !< Pressure atop the layer [Pa]. + intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_b !< Pressure below the layer [Pa]. - real, intent(in) :: alpha_ref !< A mean specific volume that is - !! subtracted out to reduce the magnitude of each of the - !! integrals [m3 kg-1]. The calculation is mathematically - !! identical with different values of alpha_ref, but alpha_ref - !! alters the effects of roundoff, and answers do change. - real, intent(in) :: dP_neglect !< A miniscule pressure change with - !! the same units as p_t (Pa?) + intent(in) :: p_b !< Pressure below the layer [R L2 T-2 ~> Pa] or [Pa]. + real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out + !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. + !! The calculation is mathematically identical with different values of + !! alpha_ref, but alpha_ref alters the effects of roundoff, and + !! answers do change. + real, intent(in) :: dP_neglect ! Pa] or [Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: bathyP !< The pressure at the bathymetry [Pa] + intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa] type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(inout) :: dza !< The change in the geopotential anomaly - !! across the layer [m2 s-2]. + !! across the layer [L2 T-2 ~> m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(inout) :: intp_dza !< The integral in pressure through the - !! layer of the geopotential anomaly relative to the anomaly - !! at the bottom of the layer [Pa m2 s-2]. + optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of + !! the geopotential anomaly relative to the anomaly at the bottom of the + !! layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & - optional, intent(inout) :: intx_dza !< The integral in x of the difference - !! between the geopotential anomaly at the top and bottom of - !! the layer divided by the x grid spacing [m2 s-2]. + optional, intent(inout) :: intx_dza !< The integral in x of the difference between + !! the geopotential anomaly at the top and bottom of the layer divided + !! by the x grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & - optional, intent(inout) :: inty_dza !< The integral in y of the difference - !! between the geopotential anomaly at the top and bottom of - !! the layer divided by the y grid spacing [m2 s-2]. + optional, intent(inout) :: inty_dza !< The integral in y of the difference between + !! the geopotential anomaly at the top and bottom of the layer divided + !! 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. + real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific + !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure + !! into Pa [Pa T2 R-1 L-2 ~> 1]. ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2437,23 +2468,31 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! Bode's rule to do the horizontal integrals, and from a truncation in the ! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. - real, dimension(5) :: T5, S5, p5, a5 - real, dimension(15) :: T15, S15, p15, a15 - real :: wt_t(5), wt_b(5) + real :: T5(5) ! Temperatures at five quadrature points [degC] + real :: S5(5) ! Salinities at five quadrature points [ppt] + real :: p5(5) ! Pressures at five quadrature points, scaled back to Pa as necessary [Pa] + real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: T15(15) ! Temperatures at fifteen interior quadrature points [degC] + real :: S15(15) ! Salinities at fifteen interior quadrature points [ppt] + real :: p15(15) ! Pressures at fifteen quadrature points, scaled back to Pa as necessary [Pa] + real :: a15(15) ! Specific volumes at fifteen quadrature points [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: wt_t(5), wt_b(5) ! Weights of top and bottom values at quadrature points [nondim] real :: T_top, T_bot, S_top, S_bot, P_top, P_bot real :: alpha_anom ! The depth averaged specific density anomaly [m3 kg-1]. - real :: dp ! The pressure change through a layer [Pa]. - real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [Pa]. - real :: hWght ! A pressure-thickness below topography [Pa]. - real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Pa]. - real :: iDenom ! The inverse of the denominator in the weights [Pa-2]. + real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]. + real :: dp_90(2:4) ! The pressure change through a layer divided by 90 [R L2 T-2 ~> Pa]. + real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]. + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]. + real :: alpha_ref_mks ! The reference specific volume in MKS units, never rescaled from m3 kg-1 [m3 kg-1] + real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]. real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]. real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. real :: intp(5) ! The integrals of specific volume with pressure at the - ! 5 sub-column locations [m2 s-2]. + ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. + real :: RL2_T2_to_Pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant. logical :: do_massWeight ! Indicates whether to do mass weighting. integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos @@ -2463,6 +2502,10 @@ 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 + RL2_T2_to_Pa = 1.0 ; if (present(pres_scale)) RL2_T2_to_Pa = pres_scale + alpha_ref_mks = alpha_ref ; if (present(SV_scale)) alpha_ref_mks = alpha_ref / SV_scale + + do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) wt_b(n) = 1.0 - wt_t(n) @@ -2474,11 +2517,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, do j=Jsq,Jeq+1; do i=Isq,Ieq+1 dp = p_b(i,j) - p_t(i,j) do n=1,5 ! T, S and p are linearly interpolated in the vertical. - p5(n) = wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j) + p5(n) = RL2_T2_to_Pa * (wt_t(n) * p_t(i,j) + wt_b(n) * p_b(i,j)) S5(n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) T5(n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref) + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) ! Use Bode's rule to estimate the interface height anomaly change. alpha_anom = C1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3)) @@ -2529,13 +2572,13 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! Salinity, temperature and pressure with linear interpolation in the vertical. pos = (m-2)*5 do n=1,5 - p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot + p15(pos+n) = RL2_T2_to_Pa * (wt_t(n) * P_top + wt_b(n) * P_bot) S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref) + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) do m=2,4 @@ -2588,13 +2631,13 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! Salinity, temperature and pressure with linear interpolation in the vertical. pos = (m-2)*5 do n=1,5 - p15(pos+n) = wt_t(n) * P_top + wt_b(n) * P_bot + p15(pos+n) = RL2_T2_to_Pa * (wt_t(n) * P_top + wt_b(n) * P_bot) S15(pos+n) = wt_t(n) * S_top + wt_b(n) * S_bot T15(pos+n) = wt_t(n) * T_top + wt_b(n) * T_bot enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref) + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) do m=2,4 diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 39d1dd26d4..cd590aa611 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -631,7 +631,7 @@ end subroutine int_density_dz_wright !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, useMassWghtInterp) + bathyP, dP_neglect, useMassWghtInterp, SV_scale, pres_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -639,53 +639,66 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: S !< Salinity [PSU]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_t !< Pressure at the top of the layer [Pa]. + intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_b !< Pressure at the top of the layer [Pa]. + intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]. real, intent(in) :: spv_ref !< A mean specific volume that is subtracted out - !! to reduce the magnitude of each of the integrals [m3 kg-1]. The calculation is - !! mathematically identical with different values of spv_ref, but this reduces the - !! effects of roundoff. + !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. + !! The calculation is mathematically identical with different values of + !! spv_ref, but this reduces the effects of roundoff. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(inout) :: dza !< The change in the geopotential anomaly across - !! the layer [m2 s-2]. + !! the layer [T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of !! the geopotential anomaly relative to the anomaly - !! at the bottom of the layer [Pa m2 s-2]. + !! at the bottom of the layer [R L4 T-4 ~> Pa m2 s-2] + !! or [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & optional, intent(inout) :: intx_dza !< The integral in x of the difference between the !! geopotential anomaly at the top and bottom of - !! the layer divided by the x grid spacing [m2 s-2]. + !! the layer divided by the x grid spacing + !! [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & optional, intent(inout) :: inty_dza !< The integral in y of the difference between the !! geopotential anomaly at the top and bottom of - !! the layer divided by the y grid spacing [m2 s-2]. + !! the layer divided by the y grid spacing + !! [L2 T-2 ~> m2 s-2] or [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate !! dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(in) :: bathyP !< The pressure at the bathymetry [Pa] + optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with - !! the same units as p_t [Pa] + !! 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. + real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific + !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure + !! into Pa [Pa T2 R-1 L-2 ~> 1]. ! Local variables real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d - real :: al0, p0, lambda - real :: p_ave - real :: rem, eps, eps2 - real :: alpha_anom ! The depth averaged specific volume anomaly [m3 kg-1]. - real :: dp ! The pressure change through a layer [Pa]. - real :: hWght ! A pressure-thickness below topography [Pa]. - real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Pa]. - real :: iDenom ! The inverse of the denominator in the weights [Pa-2]. + real :: al0 ! A term in the Wright EOS [R-1 ~> m3 kg-1] + real :: p0 ! A term in the Wright EOS [R L2 T-2 ~> Pa] + real :: lambda ! A term in the Wright EOS [L2 T-2 ~> m2 s-2] + real :: al0_scale ! Scaling factor to convert al0 from MKS units [R-1 kg m-3 ~> 1] + real :: p0_scale ! Scaling factor to convert p0 from MKS units [R L2 T-2 Pa-1 ~> 1] + real :: lam_scale ! Scaling factor to convert lambda from MKS units [L2 s2 T-2 m-2 ~> 1] + real :: p_ave ! The layer average pressure [R L2 T-2 ~> Pa] + real :: rem ! [L2 T-2 ~> m2 s-2] + real :: eps, eps2 ! A nondimensional ratio and its square [nondim] + real :: alpha_anom ! The depth averaged specific volume anomaly [R-1 ~> m3 kg-1]. + real :: dp ! The pressure change through a layer [R L2 T-2 ~> Pa]. + real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa]. + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa]. + real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2]. real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]. real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. real :: intp(5) ! The integrals of specific volume with pressure at the - ! 5 sub-column locations [m2 s-2]. + ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. @@ -697,6 +710,14 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, 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 + + al0_scale = 1.0 ; if (present(SV_scale)) al0_scale = SV_scale + p0_scale = 1.0 + if (present(pres_scale)) then ; if (pres_scale /= 1.0) then + p0_scale = 1.0 / pres_scale + endif ; endif + lam_scale = al0_scale * p0_scale + do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then do_massWeight = .true. @@ -706,10 +727,11 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & ! "dP_neglect must be present if useMassWghtInterp is present and true.") endif ; endif + ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) do j=jsh,jeh ; do i=ish,ieh - al0_2d(i,j) = (a0 + a1*T(i,j)) + a2*S(i,j) - p0_2d(i,j) = (b0 + b4*S(i,j)) + T(i,j) * (b1 + T(i,j)*((b2 + b3*T(i,j))) + b5*S(i,j)) - lambda_2d(i,j) = (c0 +c4*S(i,j)) + T(i,j) * (c1 + T(i,j)*((c2 + c3*T(i,j))) + c5*S(i,j)) + al0_2d(i,j) = al0_scale * ( (a0 + a1*T(i,j)) + a2*S(i,j) ) + p0_2d(i,j) = p0_scale * ( (b0 + b4*S(i,j)) + T(i,j) * (b1 + T(i,j)*((b2 + b3*T(i,j))) + b5*S(i,j)) ) + lambda_2d(i,j) = lam_scale * ( (c0 + c4*S(i,j)) + T(i,j) * (c1 + T(i,j)*((c2 + c3*T(i,j))) + c5*S(i,j)) ) al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dp = p_b(i,j) - p_t(i,j) diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 2c19b617c6..623db27ad3 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -510,56 +510,56 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: S !< Salinity [PSU]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_t !< Pressure at the top of the layer [Pa]. + intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - intent(in) :: p_b !< Pressure at the top of the layer [Pa]. - real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out - !! to reduce the magnitude of each of the integrals, m3 kg-1. The calculation is - !! mathematically identical with different values of alpha_ref, but this reduces the - !! effects of roundoff. - real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3]. + intent(in) :: p_b !< Pressure at the top of the layer [R L2 T-2 ~> Pa] or [Pa]. + real, intent(in) :: alpha_ref !< A mean specific volume that is subtracted out + !! to reduce the magnitude of each of the integrals [R-1 ~> m3 kg-1]. + !! The calculation is mathematically identical with different values of + !! alpha_ref, but this reduces the effects of roundoff. + real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [R ~> kg m-3] or [kg m-3]. real, intent(in) :: dRho_dT !< The derivative of density with temperature - !! [kg m-3 degC-1]. + !! [R degC-1 ~> kg m-3 degC-1] or [kg m-3 degC-1]. real, intent(in) :: dRho_dS !< The derivative of density with salinity, - !! in [kg m-3 ppt-1]. + !! in [R ppt-1 ~> kg m-3 ppt-1] or [kg m-3 ppt-1]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(out) :: dza !< The change in the geopotential anomaly across - !! the layer [m2 s-2]. + !! the layer [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(out) :: intp_dza !< The integral in pressure through the layer of - !! the geopotential anomaly relative to the anomaly - !! at the bottom of the layer [Pa m2 s-2]. + optional, intent(out) :: intp_dza !< The integral in pressure through the layer of the + !! geopotential anomaly relative to the anomaly at the + !! bottom of the layer [R L4 T-4 ~> Pa m2 s-2] or [Pa m2 s-2]. real, dimension(HI%IsdB:HI%IedB,HI%jsd:HI%jed), & optional, intent(out) :: intx_dza !< The integral in x of the difference between the !! geopotential anomaly at the top and bottom of !! the layer divided by the x grid spacing - !! [m2 s-2]. + !! [L2 T-2 ~> m2 s-2] or [m2 s-2]. real, dimension(HI%isd:HI%ied,HI%JsdB:HI%JedB), & optional, intent(out) :: inty_dza !< The integral in y of the difference between the !! geopotential anomaly at the top and bottom of !! the layer divided by the y grid spacing - !! [m2 s-2]. + !! [L2 T-2 ~> m2 s-2] or [m2 s-2]. integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & - optional, intent(in) :: bathyP !< The pressure at the bathymetry [Pa] + optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] or [Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with - !! the same units as p_t [Pa] + !! 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. ! Local variables - real :: dRho_TS ! The density anomaly due to T and S [kg m-3]. - real :: alpha_anom ! The specific volume anomaly from 1/rho_ref [m3 kg-1]. - real :: aaL, aaR ! rho_anom to the left and right [kg m-3]. - real :: dp, dpL, dpR ! Layer pressure thicknesses [Pa]. - real :: hWght ! A pressure-thickness below topography [Pa]. - real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Pa]. - real :: iDenom ! The inverse of the denominator in the weights [Pa-2]. + real :: dRho_TS ! The density anomaly due to T and S [R ~> kg m-3] or [kg m-3]. + real :: alpha_anom ! The specific volume anomaly from 1/rho_ref [R-1 ~> m3 kg-1] or [m3 kg-1]. + real :: aaL, aaR ! The specific volume anomaly to the left and right [R-1 ~> m3 kg-1] or [m3 kg-1]. + real :: dp, dpL, dpR ! Layer pressure thicknesses [R L2 T-2 ~> Pa] or [Pa]. + real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] or [Pa]. + real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa] or [Pa]. + real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-2 ~> Pa-2] or [Pa-2]. real :: hWt_LL, hWt_LR ! hWt_LA is the weighted influence of A on the left column [nondim]. real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. real :: intp(5) ! The integrals of specific volume with pressure at the - ! 5 sub-column locations [m2 s-2]. + ! 5 sub-column locations [L2 T-2 ~> m2 s-2] or [m2 s-2]. logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_6 = 1.0/6.0, C1_90 = 1.0/90.0 ! Rational constants. integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo