From cf19c9c75b16706fdbe32d92e9a3ba117762613e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 5 Aug 2024 09:50:01 -0400 Subject: [PATCH] *Non-Boussinesq code for RESET_INTXPA_INTEGRAL Added non-Boussinesq versions of the code that is exercised when CORRECTION_INTXPA_5PT and RESET_INTXPA_INTEGRAL are set to true. These options use the same names as in their Boussinesq forms, even though the arrays that are being adjusted are actually intx_za and inty_za. Also changed the units of the control structure variable that determines the surface displacement at which surface mass-weighting in the pressure gradient force is considered to be active from vertical distance units into thickness units. As a part of the testing of this commit, several checksums were added to the PressureForce_FV routines; these are enabled when DEBUG = True. The changes in this commit will change the answers in non-Boussinesq cases when either of those options are in use, but are bitwise identical when they are not, and all Boussinesq answers are bitwise identical. --- src/core/MOM_PressureForce_FV.F90 | 467 +++++++++++++++++++++++------- 1 file changed, 368 insertions(+), 99 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 453317399e..e7a9ba2d91 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -3,6 +3,7 @@ module MOM_PressureForce_FV ! This file is part of MOM6. See LICENSE.md for the license. +use MOM_debugging, only : hchksum, uvchksum use MOM_diag_mediator, only : post_data, register_diag_field use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe @@ -54,8 +55,9 @@ module MOM_PressureForce_FV 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] - real :: dz_MWIPG !< A vertical distance between in an interface and the depressed sea surface height - !! in a neighboring cell beyond which the mass-weighting is enabled [Z ~> m] + real :: h_MWIPG !< A thickness between an interface and the potentially ice-depressed + !! sea surface height or pressure in a neighboring cell beyond which the mass-weighting + !! is enabled [H ~> m or kg m-2] logical :: use_inaccurate_pgf_rho_anom !< If true, uses the older and less accurate !! method to calculate density anomalies, as used prior to !! March 2018. @@ -71,6 +73,7 @@ module MOM_PressureForce_FV !! for the finite volume pressure gradient calculation. !! By the default (1) is for a piecewise linear method + logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: use_SSH_in_Z0p !< If true, adjust the height at which the pressure used in the !! equation of state is 0 to account for the displacement of the sea !! surface including adjustments for atmospheric or sea-ice pressure. @@ -161,11 +164,46 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & inty_dza ! The change in inty_za through a layer [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJ_(G)) :: & + T_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC] + S_top, & ! Salinity of top layer used with correction_intxpa [S ~> ppt] SpV_top ! Specific volume anomaly of top layer used with correction_intxpa [R-1 ~> m3 kg-1] real, dimension(SZIB_(G),SZJ_(G)) :: & intx_za_cor ! Correction for curvature in intx_za [L2 T-2 ~> m2 s-2] real, dimension(SZI_(G),SZJB_(G)) :: & inty_za_cor ! Correction for curvature in inty_za [L2 T-2 ~> m2 s-2] + + ! These variables are used with reset_intxpa_integral. The values are taken from different + ! interfaces as a function of position. + real, dimension(SZIB_(G),SZJ_(G)) :: & + T_int_W, T_int_E, & ! Temperatures on the reference interface to the east and west of a u-point [C ~> degC] + S_int_W, S_int_E, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt] + p_int_W, p_int_E, & ! Pressures on the reference interface to the east and west of a u-point [R L2 T-2 ~> Pa] + SpV_x_W, SpV_x_E, & ! Specific volume anomalies on the reference interface to the east and west + ! of a u-point [R-1 ~> m3 kg-1] + intx_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface + ! from the value that would be obtained from assuming that pressure varies + ! linearly with depth along that interface [R L2 T-2 ~> Pa]. + dp_int_x, & ! The change in x in pressure along the reference interface [R L2 T-2 ~> Pa] + intx_za_cor_ri ! The correction to intx_za based on the reference interface calculations [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G)) :: & + T_int_S, T_int_N, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC] + S_int_S, S_int_N, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt] + p_int_S, p_int_N, & ! Pressures on the reference interface to the north and south of a v-point [R L2 T-2 ~> Pa] + SpV_y_S, SpV_y_N, & ! Specific volume anomalies on the reference interface to the north and south + ! of a v-point [R L2 T-2 ~> Pa] + inty_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface + ! from the value that would be obtained from assuming that pressure varies + ! linearly with depth along that interface [L2 T-2 ~> m2 s-2]. + dp_int_y, & ! The change in y in geopotenial height along the reference interface [R L2 T-2 ~> Pa] + inty_za_cor_ri ! The correction to inty_za based on the reference interface calculations [L2 T-2 ~> m2 s-2] + logical, dimension(SZIB_(G),SZJ_(G)) :: & + seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate + ! of the curvature terms in the intx_pa. + logical, dimension(SZI_(G),SZJB_(G)) :: & + seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate + ! of the curvature terms in the inty_pa. + + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & MassWt_u ! The fractional mass weighting at a u-point [nondim]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & @@ -182,26 +220,41 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ 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. + logical :: do_more_k ! If true, there are still points where a flatter interface remains to be found. type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S. real :: alpha_ref ! A reference specific volume [R-1 ~> m3 kg-1] that is used ! to reduce the impact of truncation errors. real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3]. + real :: dp_MWIPG ! A change in pressure between an interface and the potentially ice-enhanced + ! sea surface pressure in a neighboring cell beyond which the mass-weighting + ! is enabled [R L2 T-2 ~> Pa] real :: Pa_to_H ! A factor to convert from Pa to the thickness units (H) ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1]. 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) ! Pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa] + real :: SpV5(5) ! Specific volume anomalies at five quadrature points [R-1 ~> m3 kg-1] + real :: wt_R ! A weighting factor [nondim] + ! 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 - integer :: i, j, k + integer, dimension(2) :: EOSdom_u ! The i-computational domain for the equation of state at u-velocity points + integer, dimension(2) :: EOSdom_v ! The i-computational domain for the equation of state at v-velocity points + integer :: i, j, k, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke nkmb=GV%nk_rho_varies Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) + EOSdom_u(1) = Isq - (G%IsdB-1) ; EOSdom_u(2) = Ieq - (G%IsdB-1) + EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1) if (.not.CS%initialized) call MOM_error(FATAL, & "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.") @@ -217,6 +270,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ H_to_RL2_T2 = GV%g_Earth*GV%H_to_RZ dp_neglect = GV%g_Earth*GV%H_to_RZ * GV%H_subroundoff + dp_MWIPG = H_to_RL2_T2 * CS%h_MWIPG alpha_ref = 1.0 / CS%Rho0 I_gEarth = 1.0 / GV%g_Earth @@ -275,12 +329,14 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! and temperature across each layer. The subscripts 't' and 'b' refer ! to top and bottom values within each layer (these are the only degrees ! of freedom needed to know the linear profile). - if ( use_ALE ) then - if ( CS%Recon_Scheme == 1 ) then + if ( use_ALE .and. (CS%Recon_Scheme == 1) ) then call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) - elseif ( CS%Recon_Scheme == 2) then + elseif ( use_ALE .and. (CS%Recon_Scheme == 2) ) then call TS_PPM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) - endif + elseif (CS%reset_intxpa_integral) then + do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + T_b(i,j,k) = tv%T(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k) + enddo ; enddo ; enddo endif !$OMP parallel do default(shared) private(alpha_anom,dp) @@ -391,54 +447,108 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ enddo ; enddo enddo + if (CS%debug) then + call hchksum(za, "Pre-correction za", G%HI, haloshift=1, unscale=US%L_T_to_m_s**2) + call hchksum(p, "Pre-correction p", G%HI, haloshift=1, unscale=US%RL2_T2_to_Pa) + endif + ! With an ice-shelf or icebergs, this linearity condition might need to be applied ! to a sub-surface interface. if (CS%correction_intxpa) then - ! Determine surface specific volume for use in the pressure gradient corrections + ! Determine surface temperature and salinity for use in the pressure gradient corrections if (use_ALE .and. (CS%Recon_Scheme > 0)) then - do j=Jsq,Jeq+1 - call calculate_spec_vol(tv%T(:,j,1), tv%S(:,j,1), p(:,j,1), SpV_top(:,j), & - tv%eqn_of_state, EOSdom, spv_ref=alpha_ref) - enddo - elseif (use_EOS) then - do j=Jsq,Jeq+1 - call calculate_spec_vol(tv%T(:,j,1), tv%S(:,j,1), p(:,j,1), SpV_top(:,j), & - tv%eqn_of_state, EOSdom, spv_ref=alpha_ref) - enddo + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + T_top(i,j) = T_t(i,j,1) ; S_top(i,j) = S_t(i,j,1) + enddo ; enddo else - alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - SpV_top(i,j) = alpha_anom + T_top(i,j) = tv%T(i,j,1) ; S_top(i,j) = tv%S(i,j,1) enddo ; enddo endif - ! This version attempts to correct for hydrostatic variations in surface pressure under ice. - !$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 + 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 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 @@ -468,6 +578,129 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ enddo ; enddo enddo + if (CS%debug) then + call uvchksum("Prelim int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + call uvchksum("Prelim int[xy]_dza", intx_dza, inty_dza, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + endif + + if (CS%reset_intxpa_integral) then + ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is + ! reset intx_za there, then adjust intx_za throughout the water column. + ! Note: This option currently assumes that pressure varies quadratically along the bottom of first + ! nonvanished, nonMWIPG level. Possibly 5 point quadrature should be implemented as for the surface. + + ! Zero out the 2-d arrays that will be set from various reference interfaces. + T_int_W(:,:) = 0.0 ; S_int_W(:,:) = 0.0 ; p_int_W(:,:) = 0.0 + T_int_E(:,:) = 0.0 ; S_int_E(:,:) = 0.0 ; p_int_E(:,:) = 0.0 + intx_za_nonlin(:,:) = 0.0 ; intx_za_cor_ri(:,:) = 0.0 ; dp_int_x(:,:) = 0.0 + do j=js,je ; do I=Isq,Ieq + seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.) + enddo ; enddo + do k=1,nz-1 + do_more_k = .false. + do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then + ! Find the topmost layer for which both sides are nonvanished and MWIPG is not activated + if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i+1,j,k) > CS%h_nonvanished)) .and. & + (max(0., p(i,j,1)-p(i+1,j,K+1), p(i+1,j,1)-p(i,j,K+1)) <= dp_MWIPG)) then + ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at + ! the interface below this cell (it might have quadratic pressure dependence if sloped) + T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k) + S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k) + p_int_W(I,j) = p(i,j,K+1) ; p_int_E(I,j) = p(i+1,j,K+1) + + intx_za_nonlin(I,j) = intx_za(I,j,K+1) - 0.5*(za(i,j,K+1) + za(i+1,j,K+1)) + dp_int_x(I,j) = p(i+1,j,K+1)-p(i,j,K+1) + seek_x_cor(I,j) = .false. + else + do_more_k = .true. + endif + endif ; enddo ; enddo + if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward. + enddo + + 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 + + ! Repeat the calculations above for v-velocity points. + T_int_S(:,:) = 0.0 ; S_int_S(:,:) = 0.0 ; p_int_S(:,:) = 0.0 + T_int_N(:,:) = 0.0 ; S_int_N(:,:) = 0.0 ; p_int_N(:,:) = 0.0 + inty_za_nonlin(:,:) = 0.0 ; inty_za_cor_ri(:,:) = 0.0 ; dp_int_y(:,:) = 0.0 + do J=Jsq,Jeq ; do i=is,ie + seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.) + enddo ; enddo + do k=1,nz-1 + do_more_k = .false. + do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then + ! Find the topmost layer for which both sides are nonvanished and MWIPG is not activated + if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i,j+1,k) > CS%h_nonvanished)) .and. & + (max(0., p(i,j,1)-p(i,j+1,K+1), p(i,j+1,1)-p(i,j,K+1)) <= dp_MWIPG)) then + ! Store properties at the bottom of this cell to get a "good estimate" for intypa at + ! the interface below this cell (it might have quadratic pressure dependence if sloped) + T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k) + S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k) + p_int_S(i,J) = p(i,j,K+1) ; p_int_N(i,J) = p(i,j+1,K+1) + inty_za_nonlin(i,J) = inty_za(i,J,K+1) - 0.5*(za(i,j,K+1) + za(i,j+1,K+1)) + dp_int_y(i,J) = p(i,j+1,K+1) - p(i,j,K+1) + seek_y_cor(i,J) = .false. + else + do_more_k = .true. + endif + endif ; enddo ; enddo + if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward. + enddo + + 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) + 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) + 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 + + if (CS%debug) then + call uvchksum("Pre-reset int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + call uvchksum("int[xy]_za_cor", intx_za_cor_ri, inty_za_cor_ri, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + call uvchksum("int[xy]_za_nonlin", intx_za_nonlin, inty_za_nonlin, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + call uvchksum("dp_int_[xy]", dp_int_x, dp_int_y, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, unscale=US%RL2_T2_to_Pa) + endif + + ! 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 + intx_za(I,j,K) = intx_za(I,j,K) + intx_za_cor_ri(I,j) + enddo ; enddo ; enddo + + do K=1,nz+1 ; do J=Jsq,Jeq ; do i=is,ie + inty_za(i,J,K) = inty_za(i,J,K) + inty_za_cor_ri(i,J) + enddo ; enddo ; enddo + + if (CS%debug) then + call uvchksum("Post-reset int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2) + endif + + endif ! intx_za and inty_za have now been reset to reflect the properties of an unimpeded interface. + !$OMP parallel do default(shared) private(dp) do k=1,nz do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -628,23 +861,23 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm S_int_W, S_int_E, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt] p_int_W, p_int_E, & ! Pressures on the reference interface to the east and west of a u-point [R L2 T-2 ~> Pa] rho_x_W, rho_x_E, & ! Density anomalies on the reference interface to the east and west - ! of a u-point [R L2 T-2 ~> Pa] + ! of a u-point [R ~> kg m-3] intx_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface ! from the value that would be obtained from assuming that pressure varies ! linearly with depth along that interface [R L2 T-2 ~> Pa]. dgeo_x, & ! The change in x in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2] intx_pa_cor_ri ! The correction to intx_pa based on the reference interface calculations [R L2 T-2 ~> Pa] - real, dimension(SZIB_(G),SZJ_(G)) :: & + real, dimension(SZI_(G),SZJB_(G)) :: & T_int_S, T_int_N, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC] S_int_S, S_int_N, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt] p_int_S, p_int_N, & ! Pressures on the reference interface to the north and south of a v-point [R L2 T-2 ~> Pa] rho_y_S, rho_y_N, & ! Density anomalies on the reference interface to the north and south - ! of a v-point [R L2 T-2 ~> Pa] + ! of a v-point [R ~> kg m-3] inty_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface ! from the value that would be obtained from assuming that pressure varies ! linearly with depth along that interface [R L2 T-2 ~> Pa]. dgeo_y, & ! The change in y in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2] - inty_pa_cor_ri ! The correction to inty_pa based on the reference interface calculations [R L2 T-2 ~> Pa] + inty_pa_cor_ri ! The correction to inty_pa based on the reference interface calculations [R L2 T-2 ~> Pa] logical, dimension(SZIB_(G),SZJ_(G)) :: & seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate ! of the curvature terms in the intx_pa. @@ -667,6 +900,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & MassWt_v ! The fractional mass weighting at a v-point [nondim]. real, dimension(SZI_(G),SZJ_(G)) :: & + T_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC] + S_top, & ! Salinity of top layer used with correction_intxpa [S ~> ppt] rho_top ! Density anomaly of top layer used in calculating intx_pa_cor and inty_pa_cor real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & rho_pgf, rho_stanley_pgf ! Density [R ~> kg m-3] from EOS with and without SGS T variance @@ -689,6 +924,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1] real :: rho_ref ! The reference density [R ~> kg m-3]. real :: dz_neglect ! A minimal thickness [Z ~> m], like e. + real :: dz_MWIPG ! A vertical distance between an interface and the potentially ice-depressed + ! sea surface height in a neighboring cell beyond which the mass-weighting + ! is enabled [Z ~> m] 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]. logical :: use_p_atm ! If true, use the atmospheric pressure. @@ -702,16 +940,15 @@ 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), S5(5) ! Temperatures and salinities at five quadrature points [C ~> degC] and [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 :: 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] real, parameter :: C1_12 = 1.0/12.0 ! A rational constant [nondim] - real :: wt_R ! A weighting factor [nondim] - real :: rho_tr, rho_tl ! Store right and left densities in reset intxpa calculation [R ~> kg m-3] + 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 @@ -735,6 +972,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm G_Rho0 = GV%g_Earth / GV%Rho0 GxRho = GV%g_Earth * GV%Rho0 rho_ref = CS%Rho0 + dz_MWIPG = GV%H_to_Z * CS%h_MWIPG if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0 @@ -855,12 +1093,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! and temperature across each layer. The subscripts 't' and 'b' refer ! to top and bottom values within each layer (these are the only degrees ! of freedom needed to know the linear profile). - if ( use_ALE ) then - if ( CS%Recon_Scheme == 1 ) then - call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) - elseif ( CS%Recon_Scheme == 2 ) then - call TS_PPM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) - endif + if ( use_ALE .and. (CS%Recon_Scheme == 1) ) then + call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) + elseif ( use_ALE .and. (CS%Recon_Scheme == 2) ) then + call TS_PPM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap) + elseif (CS%reset_intxpa_integral) then + do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + T_b(i,j,k) = tv%T(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k) + enddo ; enddo ; enddo endif ! Set the surface boundary conditions on pressure anomaly and its horizontal @@ -959,37 +1199,38 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo if (CS%correction_intxpa) then - - ! Determine surface density for use in the pressure gradient corrections - if (use_ALE .and. CS%Recon_Scheme > 0) then - !$OMP parallel do default(shared) private(p_surf_EOS) - do j=Jsq,Jeq+1 - ! P_surf_EOS here is consistent with the pressure that is used in the int_density_dz routines. - do i=Isq,Ieq+1 ; p_surf_EOS(i) = -GxRho*(e(i,j,1) - Z_0p(i,j)) ; enddo - call calculate_density(T_t(:,j,1), S_t(:,j,1), p_surf_EOS, rho_top(:,j), & - tv%eqn_of_state, EOSdom, rho_ref=rho_ref) - enddo - elseif (use_EOS) then - !$OMP parallel do default(shared) private(p_surf_EOS) - do j=Jsq,Jeq+1 - ! P_surf_EOS here is consistent with the pressure that is used in the int_density_dz routines. - do i=Isq,Ieq+1 ; p_surf_EOS(i) = -GxRho*(e(i,j,1) - Z_0p(i,j)) ; enddo - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_surf_EOS, rho_top(:,j), & - tv%eqn_of_state, EOSdom, rho_ref=rho_ref) - enddo - else ! T and S are not state variables. - !$OMP parallel do default(shared) + ! Determine surface temperature and salinity for use in the pressure gradient corrections + if (use_ALE .and. (CS%Recon_Scheme > 0)) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - rho_top(i,j) = GV%Rlay(1) - rho_ref + T_top(i,j) = T_t(i,j,1) ; S_top(i,j) = S_t(i,j,1) + enddo ; enddo + else + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + T_top(i,j) = tv%T(i,j,1) ; S_top(i,j) = tv%S(i,j,1) enddo ; enddo endif + ! Determine surface density for use in the pressure gradient corrections + !$OMP parallel do default(shared) private(p_surf_EOS) + do j=Jsq,Jeq+1 + ! P_surf_EOS here is consistent with the pressure that is used in the int_density_dz routines. + do i=Isq,Ieq+1 ; p_surf_EOS(i) = -GxRho*(e(i,j,1) - Z_0p(i,j)) ; enddo + call calculate_density(T_top(:,j), S_top(:,j), p_surf_EOS, rho_top(:,j), & + tv%eqn_of_state, EOSdom, rho_ref=rho_ref) + enddo + + if (CS%debug) then + call hchksum(rho_top, "intx_pa rho_top", G%HI, haloshift=1, unscale=US%R_to_kg_m3) + call hchksum(e(:,:,1), "intx_pa e(1)", G%HI, haloshift=1, unscale=US%Z_to_m) + call hchksum(pa(:,:,1), "intx_pa pa(1)", G%HI, haloshift=1, unscale=US%RL2_T2_to_Pa) + endif + ! This version attempts to correct for hydrostatic variations in surface pressure under ice. !$OMP parallel do default(shared) private(dz_geo_sfc) do j=js,je ; do I=Isq,Ieq intx_pa_cor(I,j) = 0.0 dz_geo_sfc = GV%g_Earth * (e(i+1,j,1)-e(i,j,1)) - if (dz_geo_sfc * (rho_ref - (pa(i+1,j,1)-pa(i,j,1))*dz_geo_sfc) > 0.0) then + if ((dz_geo_sfc * rho_ref - (pa(i+1,j,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then ! The pressure/depth relationship has a positive implied density given by ! rho_implied = rho_ref - (pa(i+1,j,1)-pa(i,j,1)) / dz_geo_sfc if (-dz_geo_sfc * (pa(i+1,j,1)-pa(i,j,1)) > & @@ -998,8 +1239,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! 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_t(i,j,1) ; T5(5) = T_t(i+1,j,1) - S5(1) = S_t(i,j,1) ; S5(5) = S_t(i+1,j,1) + 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)) @@ -1029,7 +1270,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm do J=Jsq,Jeq ; do i=is,ie inty_pa_cor(i,J) = 0.0 dz_geo_sfc = GV%g_Earth * (e(i,j+1,1)-e(i,j,1)) - if (dz_geo_sfc * (rho_ref - (pa(i,j+1,1)-pa(i,j,1))*dz_geo_sfc) > 0.0) then + if ((dz_geo_sfc * rho_ref - (pa(i,j+1,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then ! The pressure/depth relationship has a positive implied density if (-dz_geo_sfc * (pa(i,j+1,1)-pa(i,j,1)) > & 0.25*((rho_top(i,j+1)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then @@ -1037,8 +1278,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! 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_t(i,j,1) ; T5(5) = T_t(i,j+1,1) - S5(1) = S_t(i,j,1) ; S5(5) = S_t(i,j+1,1) + 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)) @@ -1105,6 +1346,14 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa_cor(i,J) enddo ; enddo + + if (CS%debug) then + call uvchksum("int[xy]_pa_cor", intx_pa_cor, inty_pa_cor, G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%RL2_T2_to_Pa) + call uvchksum("int[xy]_pa(1)", intx_pa(:,:,1), inty_pa(:,:,1), G%HI, haloshift=0, & + symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%RL2_T2_to_Pa) + endif + else ! Set the surface boundary conditions on the horizontally integrated pressure anomaly, ! assuming that the surface pressure anomaly varies linearly in x and y. @@ -1151,7 +1400,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then ! Find the topmost layer for which both sides are nonvanished and MWIPG is not activated if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i+1,j,k) > CS%h_nonvanished)) .and. & - (max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1)) <= CS%dz_MWIPG)) then + (max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1)) <= dz_MWIPG)) then ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at ! the interface below this cell (it might have quadratic pressure dependence if sloped) T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k) @@ -1196,7 +1445,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then ! Find the topmost layer for which both sides are nonvanished and MWIPG is not activated if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i,j+1,k) > CS%h_nonvanished)) .and. & - (max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1)) <= CS%dz_MWIPG)) then + (max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1)) <= dz_MWIPG)) then ! Store properties at the bottom of this cell to get a "good estimate" for intypa at ! the interface below this cell (it might have quadratic pressure dependence if sloped) T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k) @@ -1401,6 +1650,8 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale ! temperature variance [nondim] integer :: default_answer_date ! Global answer date + logical :: use_temperature ! If true, temperature and salinity are used as state variables. + logical :: use_EOS ! If true, density calculated from T & S using an equation of state. logical :: useMassWghtInterp ! If true, use near-bottom mass weighting for T and S logical :: MassWghtInterpTop ! If true, use near-surface mass weighting for T and S under ice shelves logical :: MassWghtInterp_NonBous_bug ! If true, use a buggy mass weighting when non-Boussinesq @@ -1418,6 +1669,9 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, mdl = "MOM_PressureForce_FV" call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "DEBUG", CS%debug, & + "If true, write out verbose debugging data.", & + default=.false., debuggingParam=.true., do_not_log=.true.) call get_param(param_file, mdl, "RHO_PGF_REF", CS%Rho0, & "The reference density that is subtracted off when calculating pressure "//& "gradient forces. Its inverse is subtracted off of specific volumes when "//& @@ -1437,6 +1691,16 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, endif call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & "If true, calculate self-attraction and loading.", default=CS%tides) + + call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, & + "If true, Temperature and salinity are used as state variables.", & + default=.true., do_not_log=.true.) + call get_param(param_file, "MOM", "USE_EOS", use_EOS, & + "If true, density is calculated from temperature and "//& + "salinity with an equation of state. If USE_EOS is "//& + "true, ENABLE_THERMODYNAMICS must be true as well.", & + default=use_temperature, do_not_log=.true.) + call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", CS%use_SSH_in_Z0p, & "If true, include contributions from the sea surface height in the height-based "//& "pressure used in the equation of state calculations for the Boussinesq pressure "//& @@ -1468,23 +1732,28 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, if ((.not.GV%Boussinesq) .and. MassWghtInterp_NonBous_bug) & CS%MassWghtInterp = ibset(CS%MassWghtInterp, 3) ! Same as CS%MassWghtInterp + 8 - call get_param(param_file, mdl, "CORRECTION_INTXPA",CS%correction_intxpa, & + call get_param(param_file, mdl, "CORRECTION_INTXPA", CS%correction_intxpa, & "If true, use a correction for surface pressure curvature in intx_pa.", & - default = .false.) + default=.false., do_not_log=.not.use_EOS) call get_param(param_file, mdl, "CORRECTION_INTXPA_5PT", CS%correction_intxpa_5pt, & "If true, use 5point quadrature to calculate intxpa. This requires "//& - "CORRECTION_INTXPA = True.",default = .false.) + "CORRECTION_INTXPA = True.", default=.false., do_not_log=.not.use_EOS) 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.) + "Includes pressure correction.", 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, & "A minimal layer thickness that indicates that a layer is thick enough to usefully "//& "reestimate the pressure integral across the interface below.", & default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=.not.CS%reset_intxpa_integral) - call get_param(param_file, mdl, "RESET_INTXPA_MASS_WEIGHT_CLEARANCE", CS%dz_MWIPG, & + call get_param(param_file, mdl, "RESET_INTXPA_MASS_WEIGHT_CLEARANCE", CS%h_MWIPG, & "A vertical distance between in an interface and the depressed sea surface height "//& "in a neighboring cell beyond which the mass-weighting is enabled.", & - default=1.0e-8, units="m", scale=US%m_to_Z, do_not_log=.not.CS%reset_intxpa_integral) + default=1.0e-8, units="m", scale=GV%m_to_H, do_not_log=.not.CS%reset_intxpa_integral) call get_param(param_file, mdl, "USE_INACCURATE_PGF_RHO_ANOM", CS%use_inaccurate_pgf_rho_anom, & "If true, use a form of the PGF that uses the reference density "//& "in an inaccurate way. This is not recommended.", default=.false.)