Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(*)Add parentheses for FMA rotational symmetry #1634

Merged
merged 30 commits into from
Aug 24, 2024
+850 −837
Merged

(*)Add parentheses for FMA rotational symmetry #1634

Changes from 1 commit
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
6628dbf
(*)Parenthesize squares of wind stresses for FMAs
Hallberg-NOAA Mar 1, 2024
f0e61f3
(*)Parenthesize continuity_PPM curv_3 expressions
Hallberg-NOAA Feb 29, 2024
4f710ef
(*)Add parentheses for oblique OBCs with FMAs
Hallberg-NOAA Mar 1, 2024
9172cd5
(*)Add parentheses for density_integrals with FMAs
Hallberg-NOAA Mar 1, 2024
24091cc
(*)Add parentheses to 4 EOS int_density routines
Hallberg-NOAA Mar 1, 2024
8066a3d
(*)Simplify density integral parentheses
Hallberg-NOAA Mar 4, 2024
99fd957
(*)Parenthesize PressureForce_Montgomery for FMAs
Hallberg-NOAA Mar 1, 2024
307a4e2
(*)Parenthesize calc_isoneutral_slopes for FMAs
Hallberg-NOAA Mar 1, 2024
ce559ce
(*)Parenthesize MOM_calc_varT for FMAs
Hallberg-NOAA Mar 1, 2024
c344a11
(*)Parenthesize tracer_hordiff for FMAs
Hallberg-NOAA Mar 1, 2024
b2beab2
(*)Parenthesize iceberg_forces for FMAs
Hallberg-NOAA Mar 1, 2024
5398e6f
(*)Parenthesize CoriolisStokes and LA_Stk for FMAs
Hallberg-NOAA Mar 1, 2024
56d053a
+(*)Add and use G%Coriolis2Bu
Hallberg-NOAA Mar 1, 2024
49419f7
(*)Parenthesize thickness_diffuse for FMAs
Hallberg-NOAA Mar 1, 2024
03dc6f9
(*)Parenthesize Zanna_Bolton for FMAs
Hallberg-NOAA Mar 1, 2024
654cd4a
(*)Parenthesize MOM_internal_tides for FMAs
Hallberg-NOAA Mar 1, 2024
ebf02a9
(*)Parenthesize find_uv_at_h for FMAs
Hallberg-NOAA Mar 1, 2024
c0bef18
(*)Parenthesize set_viscous_ML for FMAs
Hallberg-NOAA Mar 1, 2024
0b50a15
(*)Rearrange calc_kappa_shear_vertex for FMAs
Hallberg-NOAA Mar 1, 2024
f0c52dd
(*)Parenthesize MOM_set_diffusivity for FMAs
Hallberg-NOAA Mar 1, 2024
64b851c
(*)Parenthesize CorAdCalc for FMAs
Hallberg-NOAA Mar 1, 2024
46e8b66
(*)Parenthesize MOM_barotropic for FMAs
Hallberg-NOAA Mar 1, 2024
6216fa1
(*)Parenthesize MOM_lateral_mixing_coeffs for FMAs
Hallberg-NOAA Mar 2, 2024
ffef92f
(*)Parenthesize MOM_hor_visc for FMAs
Hallberg-NOAA Apr 18, 2024
fc2af28
(*)Parenthesize initialization squares for FMAs
Hallberg-NOAA Mar 3, 2024
44f1130
(*)Parenthesize parameterization squares for FMAs
Hallberg-NOAA Mar 3, 2024
182223c
(*)Parenthesize diagnostics for FMAs
Hallberg-NOAA Apr 30, 2024
e810ac5
(*)Parenthesize tracer_advect PPM edge values
Hallberg-NOAA May 5, 2024
ffa766b
(*)More parentheses in density_integrals for FMAs
Hallberg-NOAA Jul 31, 2024
fd82861
(*)Add parentheses in end_value_h4 for FMAs
Hallberg-NOAA Aug 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
(*)Add parentheses to 4 EOS int_density routines 
  Added parentheses to 140 lines in 8 int_density_dz and int_spec_vol_dp
routines for the linear, Wright, Wright_full and Wright_red equations of state
so that they will be rotationally invariant when fused-multiply-adds are
enabled.  All answers are bitwise identical in cases without FMAs, but answers
could change with FMAs.
@Hallberg-NOAA
Hallberg-NOAA committed Jul 29, 2024

Verified

This commit was signed with the committer’s verified signature. The key has expired.
addaleax Anna Henningsen
GPG key ID: D8B9F5AEAE84E4CF
Expired
Verified
Learn about vigilant mode
commit 24091ccf02f61b92e7db1f4ac6b3d3a1d3036881
72 changes: 36 additions & 36 deletions src/equation_of_state/MOM_EOS_Wright.F90
View file
Original file line number Diff line number Diff line change
@@ -568,24 +568,24 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i+1,j)+z_b(i+1,j))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
@@ -609,24 +609,24 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i,j+1)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i,j+1)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i,j+1)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i,j+1))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i,j+1))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i,j+1)+z_b(i,j+1))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave)
@@ -808,26 +808,26 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(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))
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)))
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))))

eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
@@ -849,26 +849,26 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wt_L*al0_2d(i,j) + wt_R*al0_2d(i,j+1)
p0 = wt_L*p0_2d(i,j) + wt_R*p0_2d(i,j+1)
lambda = wt_L*lambda_2d(i,j) + wt_R*lambda_2d(i,j+1)
al0 = (wt_L*al0_2d(i,j)) + (wt_R*al0_2d(i,j+1))
p0 = (wt_L*p0_2d(i,j)) + (wt_R*p0_2d(i,j+1))
lambda = (wt_L*lambda_2d(i,j)) + (wt_R*lambda_2d(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))
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)))
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))))

eps = 0.5 * dp / (p0 + p_ave) ; eps2 = eps*eps
intp(m) = (al0 + lambda / (p0 + p_ave) - spv_ref)*dp + 2.0*eps* &
72 changes: 36 additions & 36 deletions src/equation_of_state/MOM_EOS_Wright_full.F90
View file
Original file line number Diff line number Diff line change
@@ -573,24 +573,24 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i+1,j)+z_b(i+1,j))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -614,24 +614,24 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i,j+1)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i,j+1)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i,j+1)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i,j+1))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i,j+1))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i,j+1)+z_b(i,j+1))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -815,26 +815,26 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(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))
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)))
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 * I_pterm ; eps2 = eps*eps
@@ -857,26 +857,26 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wt_L*al0_2d(i,j) + wt_R*al0_2d(i,j+1)
p0 = wt_L*p0_2d(i,j) + wt_R*p0_2d(i,j+1)
lambda = wt_L*lambda_2d(i,j) + wt_R*lambda_2d(i,j+1)
al0 = (wt_L*al0_2d(i,j)) + (wt_R*al0_2d(i,j+1))
p0 = (wt_L*p0_2d(i,j)) + (wt_R*p0_2d(i,j+1))
lambda = (wt_L*lambda_2d(i,j)) + (wt_R*lambda_2d(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))
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)))
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 * I_pterm ; eps2 = eps*eps
72 changes: 36 additions & 36 deletions src/equation_of_state/MOM_EOS_Wright_red.F90
View file
Original file line number Diff line number Diff line change
@@ -575,24 +575,24 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i+1,j))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i+1,j)+z_b(i+1,j))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i+1,j)+z_b(i+1,j)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -616,24 +616,24 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, &
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i,j+1)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i,j+1)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i,j+1)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i,j+1))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i,j+1))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(i,j+1))

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1))
p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i,j+1)+z_b(i,j+1))) - z0pres)
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
p_ave = -GxRho*(0.5*((wt_L*(z_t(i,j)+z_b(i,j))) + (wt_R*(z_t(i,j+1)+z_b(i,j+1)))) - z0pres)

I_al0 = 1.0 / al0
I_Lzz = 1.0 / ((p0 + p_ave) + lambda * I_al0)
@@ -817,26 +817,26 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wtT_L*al0_2d(i,j) + wtT_R*al0_2d(i+1,j)
p0 = wtT_L*p0_2d(i,j) + wtT_R*p0_2d(i+1,j)
lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j)
al0 = (wtT_L*al0_2d(i,j)) + (wtT_R*al0_2d(i+1,j))
p0 = (wtT_L*p0_2d(i,j)) + (wtT_R*p0_2d(i+1,j))
lambda = (wtT_L*lambda_2d(i,j)) + (wtT_R*lambda_2d(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))
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)))
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 * I_pterm ; eps2 = eps*eps
@@ -859,26 +859,26 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, &
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom
else
hWt_LL = 1.0 ; hWt_LR = 0.0 ; hWt_RR = 1.0 ; hWt_RL = 0.0
endif

intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
al0 = wt_L*al0_2d(i,j) + wt_R*al0_2d(i,j+1)
p0 = wt_L*p0_2d(i,j) + wt_R*p0_2d(i,j+1)
lambda = wt_L*lambda_2d(i,j) + wt_R*lambda_2d(i,j+1)
al0 = (wt_L*al0_2d(i,j)) + (wt_R*al0_2d(i,j+1))
p0 = (wt_L*p0_2d(i,j)) + (wt_R*p0_2d(i,j+1))
lambda = (wt_L*lambda_2d(i,j)) + (wt_R*lambda_2d(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))
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)))
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 * I_pterm ; eps2 = eps*eps
64 changes: 32 additions & 32 deletions src/equation_of_state/MOM_EOS_linear.F90
View file
Original file line number Diff line number Diff line change
@@ -356,24 +356,24 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, &
raL = (Rho_T0_S0 - rho_ref) + (dRho_dT*T(i,j) + dRho_dS*S(i,j))
raR = (Rho_T0_S0 - rho_ref) + (dRho_dT*T(i+1,j) + dRho_dS*S(i+1,j))

intx_dpa(i,j) = G_e*C1_6 * (dzL*(2.0*raL + raR) + dzR*(2.0*raR + raL))
intx_dpa(i,j) = G_e*C1_6 * ((dzL*(2.0*raL + raR)) + (dzR*(2.0*raR + raL)))
else
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom

intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j))
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i+1,j) - z_b(i+1,j)))
rho_anom = (Rho_T0_S0 - rho_ref) + &
(dRho_dT * (wtT_L*T(i,j) + wtT_R*T(i+1,j)) + &
dRho_dS * (wtT_L*S(i,j) + wtT_R*S(i+1,j)))
(dRho_dT * ((wtT_L*T(i,j)) + (wtT_R*T(i+1,j))) + &
dRho_dS * ((wtT_L*S(i,j)) + (wtT_R*S(i+1,j))))
intz(m) = G_e*rho_anom*dz
enddo
! Use Boole's rule to integrate the values.
@@ -395,24 +395,24 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, &
raL = (Rho_T0_S0 - rho_ref) + (dRho_dT*T(i,j) + dRho_dS*S(i,j))
raR = (Rho_T0_S0 - rho_ref) + (dRho_dT*T(i,j+1) + dRho_dS*S(i,j+1))

inty_dpa(i,j) = G_e*C1_6 * (dzL*(2.0*raL + raR) + dzR*(2.0*raR + raL))
inty_dpa(i,j) = G_e*C1_6 * ((dzL*(2.0*raL + raR)) + (dzR*(2.0*raR + raL)))
else
hL = (z_t(i,j) - z_b(i,j)) + dz_neglect
hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom

intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1))
dz = (wt_L*(z_t(i,j) - z_b(i,j))) + (wt_R*(z_t(i,j+1) - z_b(i,j+1)))
rho_anom = (Rho_T0_S0 - rho_ref) + &
(dRho_dT * (wtT_L*T(i,j) + wtT_R*T(i,j+1)) + &
dRho_dS * (wtT_L*S(i,j) + wtT_R*S(i,j+1)))
(dRho_dT * ((wtT_L*T(i,j)) + (wtT_R*T(i,j+1))) + &
dRho_dS * ((wtT_L*S(i,j)) + (wtT_R*S(i,j+1))))
intz(m) = G_e*rho_anom*dz
enddo
! Use Boole's rule to integrate the values.
@@ -530,26 +530,26 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, &
dRho_TS = dRho_dT*T(i+1,j) + dRho_dS*S(i+1,j)
aaR = ((1.0 - Rho_T0_S0*alpha_ref) - dRho_TS*alpha_ref) / (Rho_T0_S0 + dRho_TS)

intx_dza(i,j) = C1_6 * (2.0*(dpL*aaL + dpR*aaR) + (dpL*aaR + dpR*aaL))
intx_dza(i,j) = C1_6 * (2.0*((dpL*aaL) + (dpR*aaR)) + ((dpL*aaR) + (dpR*aaL)))
else
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom

intp(1) = dza(i,j) ; intp(5) = dza(i+1,j)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i+1,j) - p_t(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)))

dRho_TS = dRho_dT*(wtT_L*T(i,j) + wtT_R*T(i+1,j)) + &
dRho_dS*(wtT_L*S(i,j) + wtT_R*S(i+1,j))
dRho_TS = dRho_dT*((wtT_L*T(i,j)) + (wtT_R*T(i+1,j))) + &
dRho_dS*((wtT_L*S(i,j)) + (wtT_R*S(i+1,j)))
! alpha_anom = 1.0/(Rho_T0_S0 + dRho_TS)) - alpha_ref
alpha_anom = ((1.0-Rho_T0_S0*alpha_ref) - dRho_TS*alpha_ref) / (Rho_T0_S0 + dRho_TS)
intp(m) = alpha_anom*dp
@@ -575,26 +575,26 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, &
dRho_TS = dRho_dT*T(i,j+1) + dRho_dS*S(i,j+1)
aaR = ((1.0 - Rho_T0_S0*alpha_ref) - dRho_TS*alpha_ref) / (Rho_T0_S0 + dRho_TS)

inty_dza(i,j) = C1_6 * (2.0*(dpL*aaL + dpR*aaR) + (dpL*aaR + dpR*aaL))
inty_dza(i,j) = C1_6 * (2.0*((dpL*aaL) + (dpR*aaR)) + ((dpL*aaR) + (dpR*aaL)))
else
hL = (p_b(i,j) - p_t(i,j)) + dP_neglect
hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect
hWght = hWght * ( (hL-hR)/(hL+hR) )**2
iDenom = 1.0 / ( hWght*(hR + hL) + hL*hR )
hWt_LL = (hWght*hL + hR*hL) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + hR*hL) * iDenom ; hWt_RL = (hWght*hL) * iDenom
iDenom = 1.0 / ( hWght*(hR + hL) + (hL*hR) )
hWt_LL = (hWght*hL + (hR*hL)) * iDenom ; hWt_LR = (hWght*hR) * iDenom
hWt_RR = (hWght*hR + (hR*hL)) * iDenom ; hWt_RL = (hWght*hL) * iDenom

intp(1) = dza(i,j) ; intp(5) = dza(i,j+1)
do m=2,4
wt_L = 0.25*real(5-m) ; wt_R = 1.0-wt_L
wtT_L = wt_L*hWt_LL + wt_R*hWt_RL ; wtT_R = wt_L*hWt_LR + wt_R*hWt_RR
wtT_L = (wt_L*hWt_LL) + (wt_R*hWt_RL) ; wtT_R = (wt_L*hWt_LR) + (wt_R*hWt_RR)

! T, S, and p are interpolated in the horizontal. The p interpolation
! is linear, but for T and S it may be thickness weighted.
dp = wt_L*(p_b(i,j) - p_t(i,j)) + wt_R*(p_b(i,j+1) - p_t(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)))

dRho_TS = dRho_dT*(wtT_L*T(i,j) + wtT_R*T(i,j+1)) + &
dRho_dS*(wtT_L*S(i,j) + wtT_R*S(i,j+1))
dRho_TS = dRho_dT*((wtT_L*T(i,j)) + (wtT_R*T(i,j+1))) + &
dRho_dS*((wtT_L*S(i,j)) + (wtT_R*S(i,j+1)))
! alpha_anom = 1.0/(Rho_T0_S0 + dRho_TS)) - alpha_ref
alpha_anom = ((1.0-Rho_T0_S0*alpha_ref) - dRho_TS*alpha_ref) / (Rho_T0_S0 + dRho_TS)
intp(m) = alpha_anom*dp