From ed4623b43fbd2d7d65cd2de80e8be0902e68a425 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 12 Mar 2023 12:27:55 -0400 Subject: [PATCH] *Refactor calculate_specific_vol_wright_full Refactored the specific volume calculations for the WRIGHT_FULL and WRIGHT_RED equations of states for simplicity or to reduce the impacts of roundoff when removing a reference value. Also added code to multiply by the reciprocal of the denominator rather than dividing in several places in the int_spec_vol_dp routines for these same two equations of state, both for efficiency and greater consistency across optimization levels. These changes are mathematically equivalent but will change answers at roundoff with these two equations of state, but they are so new that they can not have been used yet. --- src/equation_of_state/MOM_EOS_Wright_full.F90 | 50 ++++++++++++------- src/equation_of_state/MOM_EOS_Wright_red.F90 | 50 ++++++++++++------- 2 files changed, 64 insertions(+), 36 deletions(-) diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 8b7fe6751d..3f00a92cef 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -179,20 +179,30 @@ subroutine calculate_spec_vol_array_wright(T, S, pressure, specvol, start, npts, ! Local variables real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] real :: p0 ! The pressure offset in the Wright EOS [Pa] - real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2], perhaps with + ! an offset to account for spv_ref + real :: al_TS ! The contributions of temperature and salinity to al0 [m3 kg-1] + real :: p_TSp ! A combination of the pressure and the temperature and salinity contributions to p0 [Pa] + real :: lam_000 ! A corrected offset to lambda, including contributions from spv_ref [m2 s-2] integer :: j - do j=start,start+npts-1 - al0 = a0 + (a1*T(j) + a2*S(j)) - p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) - lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) - - if (present(spv_ref)) then - specvol(j) = (lambda + (al0 - spv_ref)*(pressure(j) + p0)) / (pressure(j) + p0) - else - specvol(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) - endif - enddo + if (present(spv_ref)) then + lam_000 = c0 + (a0 - spv_ref)*b0 + do j=start,start+npts-1 + al_TS = a1*T(j) + a2*S(j) + p_TSp = pressure(j) + (b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j)))) + lambda = lam_000 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + ! This is equivalent to the expression below minus spv_ref, but less sensitive to roundoff. + specvol(j) = al_TS + (lambda + (a0 - spv_ref)*p_TSp) / (b0 + p_TSp) + enddo + else + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + specvol(j) = al0 + lambda / (pressure(j) + p0) + enddo + endif end subroutine calculate_spec_vol_array_wright !> Return the thermal/haline expansion coefficients for 1-d array inputs and outputs @@ -793,6 +803,7 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & 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 :: I_pterm ! The inverse of p0 plus p_ave [T2 R-1 L-2 ~> Pa-1]. 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]. @@ -866,9 +877,10 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dp = p_b(i,j) - p_t(i,j) p_ave = 0.5*(p_t(i,j)+p_b(i,j)) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - alpha_anom = (al0 - spv_ref) + lambda / (p0 + p_ave) + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + alpha_anom = (al0 - spv_ref) + lambda * I_pterm rem = (lambda * eps2) * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) dza(i,j) = alpha_anom*dp + 2.0*eps*rem if (present(intp_dza)) & @@ -906,9 +918,10 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) p_ave = 0.5*(wt_L*(p_t(i,j)+p_b(i,j)) + wt_R*(p_t(i+1,j)+p_b(i+1,j))) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - intp(m) = ((al0 - spv_ref) + lambda / (p0 + p_ave))*dp + 2.0*eps* & + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + intp(m) = ((al0 - spv_ref) + lambda * I_pterm)*dp + 2.0*eps* & lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) enddo ! Use Boole's rule to integrate the values. @@ -947,9 +960,10 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) p_ave = 0.5*(wt_L*(p_t(i,j)+p_b(i,j)) + wt_R*(p_t(i,j+1)+p_b(i,j+1))) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - intp(m) = ((al0 - spv_ref) + lambda / (p0 + p_ave))*dp + 2.0*eps* & + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + intp(m) = ((al0 - spv_ref) + lambda * I_pterm)*dp + 2.0*eps* & lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) enddo ! Use Boole's rule to integrate the values. diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index 4d5de35a1f..cf78ce2211 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -179,20 +179,30 @@ subroutine calculate_spec_vol_array_wright(T, S, pressure, specvol, start, npts, ! Local variables real :: al0 ! The specific volume at 0 lambda in the Wright EOS [m3 kg-1] real :: p0 ! The pressure offset in the Wright EOS [Pa] - real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2] + real :: lambda ! The sound speed squared at 0 alpha in the Wright EOS [m2 s-2], perhaps with + ! an offset to account for spv_ref + real :: al_TS ! The contributions of temperature and salinity to al0 [m3 kg-1] + real :: p_TSp ! A combination of the pressure and the temperature and salinity contributions to p0 [Pa] + real :: lam_000 ! A corrected offset to lambda, including contributions from spv_ref [m2 s-2] integer :: j - do j=start,start+npts-1 - al0 = a0 + (a1*T(j) + a2*S(j)) - p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) - lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) - - if (present(spv_ref)) then - specvol(j) = (lambda + (al0 - spv_ref)*(pressure(j) + p0)) / (pressure(j) + p0) - else - specvol(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) - endif - enddo + if (present(spv_ref)) then + lam_000 = c0 + (a0 - spv_ref)*b0 + do j=start,start+npts-1 + al_TS = a1*T(j) + a2*S(j) + p_TSp = pressure(j) + (b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j)))) + lambda = lam_000 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + ! This is equivalent to the expression below minus spv_ref, but less sensitive to roundoff. + specvol(j) = al_TS + (lambda + (a0 - spv_ref)*p_TSp) / (b0 + p_TSp) + enddo + else + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + specvol(j) = al0 + lambda / (pressure(j) + p0) + enddo + endif end subroutine calculate_spec_vol_array_wright !> Return the thermal/haline expansion coefficients for 1-d array inputs and outputs @@ -793,6 +803,7 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & 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 :: I_pterm ! The inverse of p0 plus p_ave [T2 R-1 L-2 ~> Pa-1]. 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]. @@ -866,9 +877,10 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dp = p_b(i,j) - p_t(i,j) p_ave = 0.5*(p_t(i,j)+p_b(i,j)) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - alpha_anom = (al0 - spv_ref) + lambda / (p0 + p_ave) + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + alpha_anom = (al0 - spv_ref) + lambda * I_pterm rem = (lambda * eps2) * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) dza(i,j) = alpha_anom*dp + 2.0*eps*rem if (present(intp_dza)) & @@ -906,9 +918,10 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(i+1,j)) p_ave = 0.5*(wt_L*(p_t(i,j)+p_b(i,j)) + wt_R*(p_t(i+1,j)+p_b(i+1,j))) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - intp(m) = ((al0 - spv_ref) + lambda / (p0 + p_ave))*dp + 2.0*eps* & + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + intp(m) = ((al0 - spv_ref) + lambda * I_pterm)*dp + 2.0*eps* & lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) enddo ! Use Boole's rule to integrate the values. @@ -947,9 +960,10 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(i,j+1)) p_ave = 0.5*(wt_L*(p_t(i,j)+p_b(i,j)) + wt_R*(p_t(i,j+1)+p_b(i,j+1))) + I_pterm = 1.0 / (p0 + p_ave) - eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps - intp(m) = ((al0 - spv_ref) + lambda / (p0 + p_ave))*dp + 2.0*eps* & + eps = 0.5 * dp * I_pterm ; eps2 = eps*eps + intp(m) = ((al0 - spv_ref) + lambda * I_pterm)*dp + 2.0*eps* & lambda * eps2 * (C1_3 + eps2*(0.2 + eps2*(C1_7 + C1_9*eps2))) enddo ! Use Boole's rule to integrate the values.