From 64628f71083f24f827ed663647436b590d50f1bb Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 16 Sep 2024 05:42:42 -0400 Subject: [PATCH] +(*)Eliminate CORRECTION_INTXPA_5PT Eliminated the runtime parameter CORRECTION_INTXPA_5PT by always selecting the code that would have been used if it were true. This decision was based on the relative performance of the options with this set to true or false, and on the desire to maintain consistency with the 5-point Boole's rule quadrature used elsewhere in the pressure gradient force calculations. CORRECTION_INTXPA_5PT was only added in the same PR as this commit, so it can been simply removed and there is no need to obsolete it. This will change answers if CORRECTION_INTXPA or RESET_INTXPA_INTEGRAL are true and CORRECTION_INTXPA_5PT was set to false, but otherwise answers are bitwise identical. There are fewer entries in the MOM_parameter_doc files as a result of this commit. --- src/core/MOM_PressureForce_FV.F90 | 556 ++++++++++++------------------ 1 file changed, 224 insertions(+), 332 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index dedd86554c..42e6514ab9 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -49,9 +49,14 @@ module MOM_PressureForce_FV type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. integer :: MassWghtInterp !< A flag indicating whether and how to use mass weighting in T/S interpolation - logical :: correction_intxpa !< If true, apply a correction to surface intxpa under ice. - logical :: correction_intxpa_5pt !< Use 5 point quadrature to calculate surface intxpa - logical :: reset_intxpa_integral !< In the interior, reset intxpa at a trusted cell (for ice shelf) + logical :: correction_intxpa !< If true, apply a correction to the value of intxpa at a selected + !! interface under ice, using matching at the end values along with a + !! 5-point quadrature integral of the hydrostatic pressure or height + !! changes along that interface. The selected interface is either at the + !! ocean's surface or in the interior, depending on reset_intxpa_integral. + logical :: reset_intxpa_integral !< If true and the surface displacement between adjacent cells + !! exceeds the vertical grid spacing, reset intxpa at the interface below + !! a trusted interior cell. (This often applies in ice shelf cavities.) real :: h_nonvanished !< A minimal layer thickness that indicates that a layer is thick enough !! to usefully reestimate the pressure integral across the interface !! below it [H ~> m or kg m-2] @@ -235,7 +240,6 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! real :: oneatm ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa] real, parameter :: C1_6 = 1.0/6.0 ! [nondim] - real, parameter :: C1_12 = 1.0/12.0 ! A rational constant [nondim] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state @@ -461,89 +465,52 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ endif if (CS%correction_intxpa) then - if (CS%correction_intxpa_5pt) then - ! This version makes a 5 point quadrature correction for hydrostatic variations in surface - ! pressure under ice. - !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5) - do j=js,je ; do I=Isq,Ieq - intx_za_cor(I,j) = 0.0 - dp_sfc = (p(i+1,j,1) - p(i,j,1)) - ! If the changes in pressure and height anomaly were explicable by just a hydrostatic balance, - ! the implied specific volume would be SpV_implied = alpha_ref - (dza_x / dp_x) - if (dp_sfc * (alpha_ref*dp_sfc - (za(i+1,j,1)-za(i,j,1))) > 0.0) then - T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j) - S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j) - p5(1) = p(i,j,1) ; p5(5) = p(i+1,j,1) - do m=2,4 - wt_R = 0.25*real(m-1) - T5(m) = T5(1) + (T5(5)-T5(1))*wt_R - S5(m) = S5(1) + (S5(5)-S5(1))*wt_R - p5(m) = p5(1) + (p5(5)-p5(1))*wt_R - enddo !m - call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) - ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression. - intx_za_cor(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc - ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 - endif - intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + intx_za_cor(I,j) - enddo ; enddo - !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5) - do J=Jsq,Jeq ; do i=is,ie - inty_za_cor(i,J) = 0.0 - dp_sfc = (p(i,j+1,1) - p(i,j,1)) - if (dp_sfc * (alpha_ref*dp_sfc - (za(i,j+1,1)-za(i,j,1))) > 0.0) then - ! The pressure/depth relationship has a positive implied specific volume. - T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1) - S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1) - p5(1) = p(i,j,1) ; p5(5) = p(i,j+1,1) - do m=2,4 - wt_R = 0.25*real(m-1) - T5(m) = T5(1) + (T5(5)-T5(1))*wt_R - S5(m) = S5(1) + (S5(5)-S5(1))*wt_R - p5(m) = p5(1) + (p5(5)-p5(1))*wt_R - enddo !m - call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) - ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression. - inty_za_cor(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc - endif - inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + inty_za_cor(i,J) - enddo ; enddo - else - ! This version makes a parabolic correction for hydrostatic variations in surface pressure under ice. - - ! Determine surface specific volume for use in the pressure gradient corrections - do j=Jsq,Jeq+1 - call calculate_spec_vol(T_top(:,j), S_top(:,j), p(:,j,1), SpV_top(:,j), & - tv%eqn_of_state, EOSdom, spv_ref=alpha_ref) - enddo - - !$OMP parallel do default(shared) private(dp_sfc) - do j=js,je ; do I=Isq,Ieq - intx_za_cor(I,j) = 0.0 - dp_sfc = (p(i+1,j,1) - p(i,j,1)) - ! If the changes in pressure and height anomaly were explicable by just a hydrostatic balance, - ! the implied specific volume would be SpV_implied = alpha_ref - (dza_x / dp_x) - if (dp_sfc * (alpha_ref*dp_sfc - (za(i+1,j,1)-za(i,j,1))) > 0.0) then - ! The pressure/depth relationship has a positive implied specific volume. - ! In non-Bousinesq mode, no other restrictions seem to be needed, and even the test - ! above might be unnecessary, but a test for the implied specific volume being at least - ! half the average specific volume would be: - ! if ((alpha_ref - dza / dp) > 0.25*((SpV_top(i+1,j)+SpV_top(i,j)) + 2.0*alpha_ref)) & - intx_za_cor(I,j) = C1_12 * (SpV_top(i+1,j)-SpV_top(i,j)) * dp_sfc - endif - intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + intx_za_cor(I,j) - enddo ; enddo - !$OMP parallel do default(shared) private(dp_sfc) - do J=Jsq,Jeq ; do i=is,ie - inty_za_cor(i,J) = 0.0 - dp_sfc = (p(i,j+1,1) - p(i,j,1)) - if (dp_sfc * (alpha_ref*dp_sfc - (za(i,j+1,1)-za(i,j,1))) > 0.0) then - ! The pressure/depth relationship has a positive implied specific volume. - inty_za_cor(i,J) = C1_12 * (SpV_top(i,j+1)-SpV_top(i,j)) * dp_sfc - endif - inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + inty_za_cor(i,J) - enddo ; enddo - endif + ! This version makes a 5 point quadrature correction for hydrostatic variations in surface + ! pressure under ice. + !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5) + do j=js,je ; do I=Isq,Ieq + intx_za_cor(I,j) = 0.0 + dp_sfc = (p(i+1,j,1) - p(i,j,1)) + ! If the changes in pressure and height anomaly were explicable by just a hydrostatic balance, + ! the implied specific volume would be SpV_implied = alpha_ref - (dza_x / dp_x) + if (dp_sfc * (alpha_ref*dp_sfc - (za(i+1,j,1)-za(i,j,1))) > 0.0) then + T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j) + S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j) + p5(1) = p(i,j,1) ; p5(5) = p(i+1,j,1) + do m=2,4 + wt_R = 0.25*real(m-1) + T5(m) = T5(1) + (T5(5)-T5(1))*wt_R + S5(m) = S5(1) + (S5(5)-S5(1))*wt_R + p5(m) = p5(1) + (p5(5)-p5(1))*wt_R + enddo !m + call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) + ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression. + intx_za_cor(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc + ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 + endif + intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + intx_za_cor(I,j) + enddo ; enddo + !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5) + do J=Jsq,Jeq ; do i=is,ie + inty_za_cor(i,J) = 0.0 + dp_sfc = (p(i,j+1,1) - p(i,j,1)) + if (dp_sfc * (alpha_ref*dp_sfc - (za(i,j+1,1)-za(i,j,1))) > 0.0) then + ! The pressure/depth relationship has a positive implied specific volume. + T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1) + S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1) + p5(1) = p(i,j,1) ; p5(5) = p(i,j+1,1) + do m=2,4 + wt_R = 0.25*real(m-1) + T5(m) = T5(1) + (T5(5)-T5(1))*wt_R + S5(m) = S5(1) + (S5(5)-S5(1))*wt_R + p5(m) = p5(1) + (p5(5)-p5(1))*wt_R + enddo !m + call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) + ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression. + inty_za_cor(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc + endif + inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + inty_za_cor(i,J) + enddo ; enddo else ! This order of integrating upward and then downward again is necessary with ! a nonlinear equation of state, so that the surface geopotentials will go @@ -639,42 +606,27 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ endif ; enddo ; enddo endif - if (CS%correction_intxpa_5pt) then - do j=js,je - do I=Isq,Ieq - ! This expression assumes that temperature and salinity vary linearly with pressure - ! between the corners of the reference interfaces found above to get a correction to - ! intx_pa that takes nonlinearities in the equation of state into account. - ! It is derived from a 5 point quadrature estimate of the integral with a large-scale - ! linear correction so that the pressures and heights match at the end-points. It turns - ! out that this linear correction cancels out the mid-point specific volume. - ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land. - T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j) - T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j) - T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) - S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) - p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) - call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) - - ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 - intx_za_cor_ri(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * & - dp_int_x(I,j) - intx_za_nonlin(I,j) - enddo + do j=js,je + do I=Isq,Ieq + ! This expression assumes that temperature and salinity vary linearly with pressure + ! between the corners of the reference interfaces found above to get a correction to + ! intx_pa that takes nonlinearities in the equation of state into account. + ! It is derived from a 5 point quadrature estimate of the integral with a large-scale + ! linear correction so that the pressures and heights match at the end-points. It turns + ! out that this linear correction cancels out the mid-point specific volume. + ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land. + T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j) + T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j) + T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) + S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) + p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) + call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) + + ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 + intx_za_cor_ri(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * & + dp_int_x(I,j) - intx_za_nonlin(I,j) enddo - else - do j=js,je - call calculate_spec_vol(T_int_W(:,j), S_int_W(:,j), p_int_W(:,j), SpV_x_W(:,j), & - tv%eqn_of_state, EOSdom_u, spv_ref=alpha_ref) - call calculate_spec_vol(T_int_E(:,j), S_int_E(:,j), p_int_E(:,j), SpV_x_E(:,j), & - tv%eqn_of_state, EOSdom_u, spv_ref=alpha_ref) - do I=Isq,Ieq - ! This expression assumes that specific volume varies linearly with depth between the corners of the - ! reference interfaces found above to get a vertically uniform correction to intx_za. - ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land. - intx_za_cor_ri(I,j) = C1_12 * (SpV_x_E(I,j)-SpV_x_W(I,j)) * dp_int_x(I,j) - intx_za_nonlin(I,j) - enddo - enddo - endif + enddo ! Repeat the calculations above for v-velocity points. T_int_S(:,:) = 0.0 ; S_int_S(:,:) = 0.0 ; p_int_S(:,:) = 0.0 @@ -730,42 +682,27 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ endif ; enddo ; enddo endif - if (CS%correction_intxpa_5pt) then - do J=Jsq,Jeq - do i=is,ie - ! This expression assumes that temperature and salinity vary linearly with pressure - ! between the corners of the reference interfaces found above to get a correction to - ! intx_pa that takes nonlinearities in the equation of state into account. - ! It is derived from a 5 point quadrature estimate of the integral with a large-scale - ! linear correction so that the pressures and heights match at the end-points. It turns - ! out that this linear correction cancels out the mid-point specific volume. - ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land. - T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J) - T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J) - T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) - S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) - p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) - call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) - - ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 - inty_za_cor_ri(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * & - dp_int_y(i,J) - inty_za_nonlin(i,J) - enddo + do J=Jsq,Jeq + do i=is,ie + ! This expression assumes that temperature and salinity vary linearly with pressure + ! between the corners of the reference interfaces found above to get a correction to + ! intx_pa that takes nonlinearities in the equation of state into account. + ! It is derived from a 5 point quadrature estimate of the integral with a large-scale + ! linear correction so that the pressures and heights match at the end-points. It turns + ! out that this linear correction cancels out the mid-point specific volume. + ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land. + T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J) + T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J) + T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) + S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) + p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) + call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref) + + ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 + inty_za_cor_ri(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * & + dp_int_y(i,J) - inty_za_nonlin(i,J) enddo - else - do J=Jsq,Jeq - call calculate_spec_vol(T_int_S(:,J), S_int_S(:,J), p_int_S(:,J), SpV_y_S(:,J), & - tv%eqn_of_state, EOSdom_v, spv_ref=alpha_ref) - call calculate_spec_vol(T_int_N(:,J), S_int_N(:,J), p_int_N(:,J), SpV_y_N(:,J), & - tv%eqn_of_state, EOSdom_v, spv_ref=alpha_ref) - do i=is,ie - ! This expression assumes that specific volume varies linearly with depth between the corners of the - ! reference interfaces found above to get a vertically uniform correction to inty_pa. - ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land. - inty_za_cor_ri(i,J) = C1_12 * (SpV_y_N(i,J)-SpV_y_S(i,J)) * dp_int_y(i,J) - inty_za_nonlin(i,J) - enddo - enddo - endif + enddo if (CS%debug) then call uvchksum("Pre-reset int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, & @@ -1019,6 +956,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real :: dz_neglect ! A minimal thickness [Z ~> m], like e. real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]. + real :: T5(5) ! Temperatures and salinities at five quadrature points [C ~> degC] + real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] + real :: p5(5) ! Full pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa] + real :: pa5(5) ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at five quadrature points [R L2 T-2 ~> Pa]. + real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3] + real :: wt_R ! A weighting factor [nondim] + real, parameter :: C1_6 = 1.0/6.0 ! A rational constant [nondim] + real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] logical :: use_p_atm ! If true, use the atmospheric pressure. logical :: use_ALE ! If true, use an ALE pressure reconstruction. logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. @@ -1030,15 +975,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm integer, dimension(2) :: EOSdom_v ! The i-computational domain for the equation of state at v-velocity points integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb integer :: i, j, k, m, k2 - real :: T5(5) ! Temperatures and salinities at five quadrature points [C ~> degC] - real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] - real :: p5(5) ! Full pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa] - real :: pa5(5) ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at five quadrature points [R L2 T-2 ~> Pa]. - real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3] - real :: wt_R ! A weighting factor [nondim] - real, parameter :: C1_6 = 1.0/6.0 ! A rational constant [nondim] - real, parameter :: C1_12 = 1.0/12.0 ! A rational constant [nondim] - real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke nkmb=GV%nk_rho_varies @@ -1328,31 +1264,27 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm 0.25*((rho_top(i+1,j)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then ! The pressure difference is at least half the size of the difference expected by hydrostatic ! balance. This test gets rid of pressure differences that are small, e.g. open ocean. - if (CS%correction_intxpa_5pt) then - ! Use 5 point quadrature to calculate intxpa - T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j) - S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j) - pa5(1) = pa(i,j,1) ; pa5(5) = pa(i+1,j,1) - ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines. - p5(1) = -GxRho*(e(i,j,1) - Z_0p(i,j)) - p5(5) = -GxRho*(e(i+1,j,1) - Z_0p(i,j)) - do m=2,4 - wt_R = 0.25*real(m-1) - T5(m) = T5(1) + (T5(5)-T5(1))*wt_R - S5(m) = S5(1) + (S5(5)-S5(1))*wt_R - p5(m) = p5(1) + (p5(5)-p5(1))*wt_R - enddo !m - call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) - - ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure - ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule - ! quadrature to find the integrated correction to the integral of pressure along the interface. - ! The derivation for this expression is shown below in the y-direction version. - intx_pa_cor(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc - ! Note that (4.75 + 5.5/2) / 90 = 1/12, so this is consistent with the linear result below. - else ! Do not use 5-point quadrature. - intx_pa_cor(I,j) = C1_12 * (rho_top(i+1,j)-rho_top(i,j)) * dz_geo_sfc - endif + ! Use 5 point quadrature to calculate intxpa + T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j) + S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j) + pa5(1) = pa(i,j,1) ; pa5(5) = pa(i+1,j,1) + ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines. + p5(1) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p5(5) = -GxRho*(e(i+1,j,1) - Z_0p(i,j)) + do m=2,4 + wt_R = 0.25*real(m-1) + T5(m) = T5(1) + (T5(5)-T5(1))*wt_R + S5(m) = S5(1) + (S5(5)-S5(1))*wt_R + p5(m) = p5(1) + (p5(5)-p5(1))*wt_R + enddo !m + call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) + + ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure + ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule + ! quadrature to find the integrated correction to the integral of pressure along the interface. + ! The derivation for this expression is shown below in the y-direction version. + intx_pa_cor(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc + ! Note that (4.75 + 5.5/2) / 90 = 1/12, so this is consistent with the linear result below. endif endif intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) + intx_pa_cor(I,j) @@ -1367,72 +1299,67 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm 0.25*((rho_top(i,j+1)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then ! The pressure difference is at least half the size of the difference expected by hydrostatic ! balance. This test gets rid of pressure differences that are small, e.g. open ocean. - if (CS%correction_intxpa_5pt) then - ! Use 5 point quadrature to calculate intypa - T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1) - S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1) - pa5(1) = pa(i,j,1) ; pa5(5) = pa(i,j+1,1) - ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines. - p5(1) = -GxRho*(e(i,j,1) - Z_0p(i,j)) - p5(5) = -GxRho*(e(i,j+1,1) - Z_0p(i,j)) - - do m=2,4 - wt_R = 0.25*real(m-1) - T5(m) = T5(1) + (T5(5)-T5(1))*wt_R - S5(m) = S5(1) + (S5(5)-S5(1))*wt_R - p5(m) = p5(1) + (p5(5)-p5(1))*wt_R - enddo !m - call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) - - ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure - ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule - ! quadrature to find the integrated correction to the integral of pressure along the interface. - inty_pa_cor(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc - - ! The derivation of this correction follows: - - ! Make pressure curvature a difference from the linear fit of pressure between the two points - ! (which is equivalent to taking 4 trapezoidal rule integrals of the hydrostatic equation on - ! sub-segments), with a constant slope that is chosen so that the pressure anomalies at the - ! two ends of the segment agree with their known values. - ! d_geo_8 = 0.125*dz_geo_sfc - ! dpa_subseg = 0.25*(pa5(5)-pa5(1)) + & - ! 0.25*d_geo_8 * ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) - ! do m=2,4 - ! pa5(m) = pa5(m-1) + dpa_subseg - d_geo_8*(r5(m)+r5(m-1))) - ! enddo - - ! Explicitly finding expressions for the incremental pressures from the recursion relation above: - ! pa5(2) = 0.25*(3.*pa5(1) + pa5(5)) + 0.25*d_geo_8 * ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) ) - ! ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + 0.25*d_geo_8 * & - ! ! ( (r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)) + & - ! ! (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) - 4.*(r5(3)+r5(2)) ) - ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + d_geo_8 * (0.5*(r5(5)-r5(1)) + (r5(4)-r5(2)) ) - ! ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * & - ! ! (2.0*(r5(5)-r5(1)) + 4.0*(r5(4)-r5(2)) + (r5(5)+r5(1)) + & - ! ! 2.0*(r5(4)+r5(2)) + 2.0*r5(3) - 4.*(r5(4)+r5(3))) - ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) ) - ! ! pa5(5) = pa5(5) + 0.25*d_geo_8 * & - ! ! ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) + & - ! ! ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) - 4.*(r5(5)+r5(4)) ) - ! pa5(5) = pa5(5) ! As it should. - - ! From these: - ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + 0.25*d_geo_8 * & - ! ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) + (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) - ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + d_geo_8 * ( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) - - ! Get the correction from the difference between the 5-point quadrature integral of pa5 and - ! its trapezoidal rule integral as: - ! inty_pa_cor(i,J) = C1_90*(7.0*(pa5(1)+pa5(5)) + 32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 0.5*(pa5(1)+pa5(5))) - ! inty_pa_cor(i,J) = C1_90*((32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 38.0*(pa5(1)+pa5(5))) - ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ((32.0*( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) + & - ! (6.*(r5(5)-r5(1)) + 12.0*(r5(4)-r5(2)) )) - ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ( 38.0*(r5(5)-r5(1)) + 44.0*(r5(4)-r5(2)) ) - - else ! Do not use 5-point quadrature. - inty_pa_cor(i,J) = C1_12 * (rho_top(i,j+1)-rho_top(i,j)) * dz_geo_sfc - endif + ! Use 5 point quadrature to calculate intypa + T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1) + S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1) + pa5(1) = pa(i,j,1) ; pa5(5) = pa(i,j+1,1) + ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines. + p5(1) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p5(5) = -GxRho*(e(i,j+1,1) - Z_0p(i,j)) + + do m=2,4 + wt_R = 0.25*real(m-1) + T5(m) = T5(1) + (T5(5)-T5(1))*wt_R + S5(m) = S5(1) + (S5(5)-S5(1))*wt_R + p5(m) = p5(1) + (p5(5)-p5(1))*wt_R + enddo !m + call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) + + ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure + ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule + ! quadrature to find the integrated correction to the integral of pressure along the interface. + inty_pa_cor(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc + + ! The derivation of this correction follows: + + ! Make pressure curvature a difference from the linear fit of pressure between the two points + ! (which is equivalent to taking 4 trapezoidal rule integrals of the hydrostatic equation on + ! sub-segments), with a constant slope that is chosen so that the pressure anomalies at the + ! two ends of the segment agree with their known values. + ! d_geo_8 = 0.125*dz_geo_sfc + ! dpa_subseg = 0.25*(pa5(5)-pa5(1)) + & + ! 0.25*d_geo_8 * ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) + ! do m=2,4 + ! pa5(m) = pa5(m-1) + dpa_subseg - d_geo_8*(r5(m)+r5(m-1))) + ! enddo + + ! Explicitly finding expressions for the incremental pressures from the recursion relation above: + ! pa5(2) = 0.25*(3.*pa5(1) + pa5(5)) + 0.25*d_geo_8 * ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) ) + ! ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + 0.25*d_geo_8 * & + ! ! ( (r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)) + & + ! ! (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) - 4.*(r5(3)+r5(2)) ) + ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + d_geo_8 * (0.5*(r5(5)-r5(1)) + (r5(4)-r5(2)) ) + ! ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * & + ! ! (2.0*(r5(5)-r5(1)) + 4.0*(r5(4)-r5(2)) + (r5(5)+r5(1)) + & + ! ! 2.0*(r5(4)+r5(2)) + 2.0*r5(3) - 4.*(r5(4)+r5(3))) + ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) ) + ! ! pa5(5) = pa5(5) + 0.25*d_geo_8 * & + ! ! ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) + & + ! ! ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) - 4.*(r5(5)+r5(4)) ) + ! pa5(5) = pa5(5) ! As it should. + + ! From these: + ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + 0.25*d_geo_8 * & + ! ( (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) + (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) + ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + d_geo_8 * ( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) + + ! Get the correction from the difference between the 5-point quadrature integral of pa5 and + ! its trapezoidal rule integral as: + ! inty_pa_cor(i,J) = C1_90*(7.0*(pa5(1)+pa5(5)) + 32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 0.5*(pa5(1)+pa5(5))) + ! inty_pa_cor(i,J) = C1_90*((32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 38.0*(pa5(1)+pa5(5))) + ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ((32.0*( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) + & + ! (6.*(r5(5)-r5(1)) + 12.0*(r5(4)-r5(2)) )) + ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ( 38.0*(r5(5)-r5(1)) + 44.0*(r5(4)-r5(2)) ) endif endif inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa_cor(i,J) @@ -1537,42 +1464,27 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif ; enddo ; enddo endif - if (CS%correction_intxpa_5pt) then - do j=js,je - do I=Isq,Ieq - ! This expression assumes that temperature and salinity vary linearly with hieght - ! between the corners of the reference interfaces found above to get a correction to - ! intx_pa that takes nonlinearities in the equation of state into account. - ! It is derived from a 5 point quadrature estimate of the integral with a large-scale - ! linear correction so that the pressures and heights match at the end-points. It turns - ! out that this linear correction cancels out the mid-point density anomaly. - ! This can be used without masking because dgeo_x and intx_pa_nonlin are 0 over land. - T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j) - T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j) - T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) - S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) - p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) - call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) - - ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 - intx_pa_cor_ri(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_x(I,j) - & - intx_pa_nonlin(I,j) - enddo - enddo - else - do j=js,je - call calculate_density(T_int_W(:,j), S_int_W(:,j), p_int_W(:,j), rho_x_W(:,j), & - tv%eqn_of_state, EOSdom_u, rho_ref=rho_ref) - call calculate_density(T_int_E(:,j), S_int_E(:,j), p_int_E(:,j), rho_x_E(:,j), & - tv%eqn_of_state, EOSdom_u, rho_ref=rho_ref) - do I=Isq,Ieq - ! This expression assumes that density varies linearly with depth between the corners of the - ! reference interfaces found above to get a vertically uniform correction to intx_pa. - ! This can be used without masking because dgeo_x and intx_pa_nonlin are 0 over land. - intx_pa_cor_ri(I,j) = C1_12 * (rho_x_E(I,j)-rho_x_W(I,j)) * dgeo_x(I,j) - intx_pa_nonlin(I,j) - enddo + do j=js,je + do I=Isq,Ieq + ! This expression assumes that temperature and salinity vary linearly with hieght + ! between the corners of the reference interfaces found above to get a correction to + ! intx_pa that takes nonlinearities in the equation of state into account. + ! It is derived from a 5 point quadrature estimate of the integral with a large-scale + ! linear correction so that the pressures and heights match at the end-points. It turns + ! out that this linear correction cancels out the mid-point density anomaly. + ! This can be used without masking because dgeo_x and intx_pa_nonlin are 0 over land. + T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j) + T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j) + T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) + S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) + p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) + call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) + + ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 + intx_pa_cor_ri(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_x(I,j) - & + intx_pa_nonlin(I,j) enddo - endif + enddo ! Repeat the calculations above for v-velocity points. T_int_S(:,:) = 0.0 ; S_int_S(:,:) = 0.0 ; p_int_S(:,:) = 0.0 @@ -1633,42 +1545,27 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif ; enddo ; enddo endif - if (CS%correction_intxpa_5pt) then - do J=Jsq,Jeq - do i=is,ie - ! This expression assumes that temperature and salinity vary linearly with hieght - ! between the corners of the reference interfaces found above to get a correction to - ! intx_pa that takes nonlinearities in the equation of state into account. - ! It is derived from a 5 point quadrature estimate of the integral with a large-scale - ! linear correction so that the pressures and heights match at the end-points. It turns - ! out that this linear correction cancels out the mid-point density anomaly. - ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land. - T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J) - T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J) - T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) - S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) - p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) - call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) - - ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 - inty_pa_cor_ri(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_y(i,J) - & - inty_pa_nonlin(i,J) - enddo + do J=Jsq,Jeq + do i=is,ie + ! This expression assumes that temperature and salinity vary linearly with hieght + ! between the corners of the reference interfaces found above to get a correction to + ! intx_pa that takes nonlinearities in the equation of state into account. + ! It is derived from a 5 point quadrature estimate of the integral with a large-scale + ! linear correction so that the pressures and heights match at the end-points. It turns + ! out that this linear correction cancels out the mid-point density anomaly. + ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land. + T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J) + T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J) + T5(2) = 0.25*(3.0*T5(1) + T5(5)) ; T5(4) = 0.25*(3.0*T5(5) + T5(1)) ; T5(3) = 0.5*(T5(5) + T5(1)) + S5(2) = 0.25*(3.0*S5(1) + S5(5)) ; S5(4) = 0.25*(3.0*S5(5) + S5(1)) ; S5(3) = 0.5*(S5(5) + S5(1)) + p5(2) = 0.25*(3.0*p5(1) + p5(5)) ; p5(4) = 0.25*(3.0*p5(5) + p5(1)) ; p5(3) = 0.5*(p5(5) + p5(1)) + call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref) + + ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12 + inty_pa_cor_ri(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_y(i,J) - & + inty_pa_nonlin(i,J) enddo - else - do J=Jsq,Jeq - call calculate_density(T_int_S(:,J), S_int_S(:,J), p_int_S(:,J), rho_y_S(:,J), & - tv%eqn_of_state, EOSdom_v, rho_ref=rho_ref) - call calculate_density(T_int_N(:,J), S_int_N(:,J), p_int_N(:,J), rho_y_N(:,J), & - tv%eqn_of_state, EOSdom_v, rho_ref=rho_ref) - do i=is,ie - ! This expression assumes that density varies linearly with depth between the corners of the - ! reference interfaces found above to get a vertically uniform correction to inty_pa. - ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land. - inty_pa_cor_ri(i,J) = C1_12 * (rho_y_N(i,J)-rho_y_S(i,J)) * dgeo_y(i,J) - inty_pa_nonlin(i,J) - enddo - enddo - endif + enddo ! Correct intx_pa and inty_pa at each interface using vertically constant corrections. do K=1,nz+1 ; do j=js,je ; do I=Isq,Ieq @@ -1931,13 +1828,8 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL", CS%reset_intxpa_integral, & "If true, reset INTXPA to match pressures at first nonvanished cell. "//& "Includes pressure correction.", default=.false., do_not_log=.not.use_EOS) - call get_param(param_file, mdl, "CORRECTION_INTXPA_5PT", CS%correction_intxpa_5pt, & - "If true, use 5-point quadrature to calculate the corrections to intxpa or intxza. "//& - "This option only acts if CORRECTION_INTXPA = True or RESET_INTXPA_INTEGRAL = True.", & - default=.false., do_not_log=.not.use_EOS) if (.not.use_EOS) then ! These options do nothing without an equation of state. CS%correction_intxpa = .false. - CS%correction_intxpa_5pt = .false. CS%reset_intxpa_integral = .false. endif call get_param(param_file, mdl, "RESET_INTXPA_H_NONVANISHED", CS%h_nonvanished, &