diff --git a/.testing/tc1/MOM_input b/.testing/tc1/MOM_input index 151c093ff9..04204bd67f 100644 --- a/.testing/tc1/MOM_input +++ b/.testing/tc1/MOM_input @@ -278,6 +278,12 @@ BOUND_CORIOLIS = True ! [Boolean] default = False ! have no effect on the SADOURNY Coriolis scheme if it ! were possible to use centered difference thickness fluxes. +! === module MOM_PressureForce_FV === +MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False + ! If true, use mass weighting when interpolating T/S for integrals + ! near the bathymetry in FV pressure gradient calculations. + + ! === module MOM_hor_visc === AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0 ! The velocity scale which is multiplied by the cube of diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index c7d2a35aa6..e142c64ff8 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -302,11 +302,16 @@ PGF_STANLEY_T2_DET_COEFF = -1.0 ! [nondim] default = -1.0 ! gradient in the deterministic part of the Stanley form of the Brankart ! correction. Negative values disable the scheme. +! === module MOM_PressureForce_FV === +MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False + ! If true, use mass weighting when interpolating T/S for integrals + ! near the bathymetry in FV pressure gradient calculations. + ! === module MOM_hor_visc === LAPLACIAN = True KH_VEL_SCALE = 0.05 SMAGORINSKY_KH = True ! [Boolean] default = False -SMAG_LAP_CONST = 0.06 ! [nondim] default = 0.0 +SMAG_LAP_CONST = 0.06 ! [nondim] default = 0.0 AH_VEL_SCALE = 0.05 ! [m s-1] default = 0.0 ! The velocity scale which is multiplied by the cube of ! the grid spacing to calculate the Laplacian viscosity. diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index 591ed4c788..c4ef8475a9 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -269,6 +269,9 @@ BOUND_CORIOLIS = True ! [Boolean] default = False ! === module MOM_PressureForce === ! === module MOM_PressureForce_FV === +MASS_WEIGHT_IN_PRESSURE_GRADIENT = True ! [Boolean] default = False + ! If true, use mass weighting when interpolating T/S for integrals + ! near the bathymetry in FV pressure gradient calculations. RECONSTRUCT_FOR_PRESSURE = False ! [Boolean] default = True ! If True, use vertical reconstruction of T & S within the integrals of the FV ! pressure gradient calculation. If False, use the constant-by-layer algorithm. diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 7b970f5686..42e6514ab9 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 @@ -15,7 +16,7 @@ module MOM_PressureForce_FV use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain +use MOM_EOS, only : calculate_density, calculate_spec_vol, EOS_domain use MOM_density_integrals, only : int_density_dz, int_specific_vol_dp use MOM_density_integrals, only : int_density_dz_generic_plm, int_density_dz_generic_ppm use MOM_density_integrals, only : int_spec_vol_dp_generic_plm @@ -48,6 +49,17 @@ 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 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] 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. @@ -63,6 +75,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. @@ -152,12 +165,54 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2]. 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)) :: & MassWt_v ! The fractional mass weighting at a v-point [nondim]. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). + real :: dp_sfc ! The change in surface pressure between adjacent cells [R L2 T-2 ~> Pa] real :: dp_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [R L2 T-2 ~> Pa]. @@ -167,6 +222,7 @@ 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 @@ -176,16 +232,27 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ ! [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_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.") @@ -259,12 +326,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) @@ -277,25 +346,25 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, & tv%eqn_of_state, US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), & - MassWghtInterp=CS%MassWghtInterp) + P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp) elseif ( CS%Recon_Scheme == 2 ) then call MOM_error(FATAL, "PressureForce_FV_nonBouss: "//& "int_spec_vol_dp_generic_ppm does not exist yet.") ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), & ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), & - ! intx_dza(:,:,k), inty_dza(:,:,k)) + ! intx_dza(:,:,k), inty_dza(:,:,k), P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp) endif else call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), & p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & - inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & + inty_dza(:,:,k), bathyP=p(:,:,nz+1), P_surf=p(:,:,1), dP_tiny=dp_neglect, & MassWghtInterp=CS%MassWghtInterp) endif if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) & - call diagnose_mass_weight_p(p(:,:,K), p(:,:,K+1), dp_neglect, p(:,:,nz+1), G%HI, & - MassWt_u(:,:,k), MassWt_v(:,:,k)) + call diagnose_mass_weight_p(p(:,:,K), p(:,:,K+1), p(:,:,nz+1), p(:,:,1), dp_neglect, CS%MassWghtInterp, & + G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k)) else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -375,28 +444,95 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ enddo ; enddo enddo - ! This order of integrating upward and then downward again is necessary with - ! a nonlinear equation of state, so that the surface geopotentials will go - ! linearly between the values at thickness points, but the bottom geopotentials - ! will not now be linear at the sub-grid-scale. Doing this ensures no motion - ! with flat isopycnals, even with a nonlinear equation of state. + 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. - !$OMP parallel do default(shared) - do j=js,je ; do I=Isq,Ieq - intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) - enddo ; enddo + if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then + ! 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 + 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 + endif + + if (CS%correction_intxpa) 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 order of integrating upward and then downward again is necessary with + ! a nonlinear equation of state, so that the surface geopotentials will go + ! linearly between the values at thickness points, but the bottom geopotentials + ! will not now be linear at the sub-grid-scale. Doing this ensures no motion + ! with flat isopycnals, even with a nonlinear equation of state. + !$OMP parallel do default(shared) + do j=js,je ; do I=Isq,Ieq + intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + enddo ; enddo + !$OMP parallel do default(shared) + do J=Jsq,Jeq ; do i=is,ie + inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + enddo ; enddo + endif + do k=1,nz !$OMP parallel do default(shared) do j=js,je ; do I=Isq,Ieq intx_za(I,j,K+1) = intx_za(I,j,K) - intx_dza(I,j,k) enddo ; enddo enddo - - !$OMP parallel do default(shared) - do J=Jsq,Jeq ; do i=is,ie - inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) - enddo ; enddo do k=1,nz !$OMP parallel do default(shared) do J=Jsq,Jeq ; do i=is,ie @@ -404,6 +540,197 @@ 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. + + ! 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 j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then + if ((p(i+1,j,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i+1,j,1))) then + ! This is the typical case in the open ocean, so use the topmost interface. + T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j) + S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j) + p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1) + intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1)) + dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1) + seek_x_cor(I,j) = .false. + endif + endif ; enddo ; enddo + + do k=1,nz + 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 mass-weighting is not + ! activated in the subgrid interpolation. + 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)) <= 0.0)) 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 + + if (do_more_k) then + ! There are still points where a correction is needed, so use the top interface. + do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then + T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j) + S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j) + p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1) + intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1)) + dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1) + seek_x_cor(I,j) = .false. + endif ; enddo ; enddo + endif + + 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 + 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 J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then + if ((p(i,j+1,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i,j+1,1))) then + ! This is the typical case in the open ocean, so use the topmost interface. + T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1) + S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1) + p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1) + inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1)) + dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1) + seek_y_cor(i,J) = .false. + endif + endif ; enddo ; enddo + + do k=1,nz + 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 mass-weighting is not + ! activated in the subgrid interpolation. + 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)) <= 0.0)) 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 + + if (do_more_k) then + ! There are still points where a correction is needed, so use the top interface. + do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then + T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1) + S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1) + p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1) + inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1)) + dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1) + seek_y_cor(i,J) = .false. + endif ; enddo ; enddo + endif + + 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 + 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 @@ -552,6 +879,41 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm ! interface atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & inty_dpa ! The change in inty_pa through a layer [R L2 T-2 ~> Pa]. + real, dimension(SZIB_(G),SZJ_(G)) :: & + intx_pa_cor ! Correction for curvature in intx_pa [R L2 T-2 ~> Pa] + real, dimension(SZI_(G),SZJB_(G)) :: & + inty_pa_cor ! Correction for curvature in inty_pa [R L2 T-2 ~> Pa] + + ! 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] + rho_x_W, rho_x_E, & ! Density anomalies on the reference interface to the east and west + ! 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(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 ~> 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] + 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(SZI_(G),SZJ_(G),SZK_(GV)), target :: & T_tmp, & ! Temporary array of temperatures where layers that are lighter @@ -567,6 +929,10 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm MassWt_u ! The fractional mass weighting at a u-point [nondim]. 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 ! in Stanley parameterization. @@ -576,7 +942,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3]. real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar). + real :: p_surf_EOS(SZI_(G)) ! The pressure at the ocean surface determined from the surface height, + ! consistent with what is used in the density integral routines [R L2 T-2 ~> Pa] real :: p0(SZI_(G)) ! An array of zeros to use for pressure [R L2 T-2 ~> Pa]. + real :: dz_geo_sfc ! The change in surface geopotential height between adjacent cells [L2 T-2 ~> m2 s-2] + real :: GxRho ! The gravitational acceleration times density [R L2 Z-1 T-2 ~> Pa m-1] real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m]. real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]. @@ -586,20 +956,32 @@ 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. + 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, parameter :: C1_6 = 1.0/6.0 ! [nondim] integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer, dimension(2) :: EOSdom_h ! The i-computational domain for the equation of state at tracer points + 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 :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb - integer :: i, j, k + integer :: i, j, k, m, k2 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_Bouss: Module must be initialized before it is used.") @@ -614,6 +996,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm dz_neglect = GV%dZ_subroundoff I_Rho0 = 1.0 / GV%Rho0 G_Rho0 = GV%g_Earth / GV%Rho0 + GxRho = GV%g_Earth * GV%Rho0 rho_ref = CS%Rho0 if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then @@ -735,12 +1118,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 @@ -749,12 +1134,12 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (use_p_atm) then !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j,1) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref) + p_atm(i,j) + pa(i,j,1) = GxRho*(e(i,j,1) - G%Z_ref) + p_atm(i,j) enddo ; enddo else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - pa(i,j,1) = (rho_ref*GV%g_Earth)*(e(i,j,1) - G%Z_ref) + pa(i,j,1) = GxRho*(e(i,j,1) - G%Z_ref) enddo ; enddo endif @@ -800,7 +1185,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), & - intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, dz_neglect, & + intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, e(:,:,1), dz_neglect, & CS%MassWghtInterp, Z_0p=Z_0p) endif if (GV%Z_to_H /= 1.0) then @@ -810,8 +1195,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo endif if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) & - call diagnose_mass_weight_Z(e(:,:,K), e(:,:,K+1), dz_neglect, G%bathyT, G%HI, & - MassWt_u(:,:,k), MassWt_v(:,:,k)) + call diagnose_mass_weight_Z(e(:,:,K), e(:,:,K+1), G%bathyT, e(:,:,1), dz_neglect, CS%MassWghtInterp, & + G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k)) else !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -838,25 +1223,176 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo enddo - ! Set the surface boundary conditions on the horizontally integrated pressure anomaly, - ! assuming that the surface pressure anomaly varies linearly in x and y. - ! If there is an ice-shelf or icebergs, this linear variation would need to be applied - ! to an interior interface. - !$OMP parallel do default(shared) - do j=js,je ; do I=Isq,Ieq - intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) - enddo ; enddo + if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then + ! 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 + 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 + endif + + if (CS%correction_intxpa) then + ! 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 + ! 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)) > & + 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. + ! 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) + enddo ; enddo + !$OMP parallel do default(shared) private(dz_geo_sfc) + 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 + ! 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 + ! 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. + ! 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) + 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. + ! If there is an ice-shelf or icebergs, this linear variation would need to be applied + ! to an interior interface. + !$OMP parallel do default(shared) + do j=js,je ; do I=Isq,Ieq + intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) + enddo ; enddo + !$OMP parallel do default(shared) + do J=Jsq,Jeq ; do i=is,ie + inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + enddo ; enddo + endif + do k=1,nz !$OMP parallel do default(shared) do j=js,je ; do I=Isq,Ieq intx_pa(I,j,K+1) = intx_pa(I,j,K) + intx_dpa(I,j,k) enddo ; enddo enddo - - !$OMP parallel do default(shared) - do J=Jsq,Jeq ; do i=is,ie - inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) - enddo ; enddo do k=1,nz !$OMP parallel do default(shared) do J=Jsq,Jeq ; do i=is,ie @@ -864,6 +1400,183 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm enddo ; enddo enddo + if (CS%reset_intxpa_integral) then + ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is + ! reset intxpa there, then adjust intxpa throughout the water column. + + ! 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_pa_nonlin(:,:) = 0.0 ; dgeo_x(:,:) = 0.0 ; intx_pa_cor_ri(:,:) = 0.0 + do j=js,je ; do I=Isq,Ieq + seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.) + enddo ; enddo + + do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then + if ((e(i+1,j,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i+1,j,1))) then + ! This is a typical case in the open ocean, so use the topmost interface. + T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j) + S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j) + p_int_W(I,j) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p_int_E(I,j) = -GxRho*(e(i+1,j,1) - Z_0p(i,j)) + intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1)) + dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1)) + seek_x_cor(I,j) = .false. + endif + endif ; enddo ; enddo + + do k=1,nz + 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 mass-weighting is not + ! activated in the subgrid interpolation. + 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)) <= 0.0)) 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) + ! These pressures are only used for the equation of state, and are only a function of + ! height, consistent with the expressions in the int_density_dz routines. + p_int_W(I,j) = -GxRho*(e(i,j,K+1) - Z_0p(i,j)) + p_int_E(I,j) = -GxRho*(e(i+1,j,K+1) - Z_0p(i,j)) + + intx_pa_nonlin(I,j) = intx_pa(I,j,K+1) - 0.5*(pa(i,j,K+1) + pa(i+1,j,K+1)) + dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,K+1)-e(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 + + if (do_more_k) then + ! There are still points where a correction is needed, so use the top interface for lack of a better idea? + do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then + T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j) + S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j) + p_int_W(I,j) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p_int_E(I,j) = -GxRho*(e(i+1,j,1) - Z_0p(i,j)) + intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1)) + dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1)) + seek_x_cor(I,j) = .false. + endif ; enddo ; enddo + endif + + 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 + + ! 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_pa_nonlin(:,:) = 0.0 ; dgeo_y(:,:) = 0.0 ; inty_pa_cor_ri(:,:) = 0.0 + do J=Jsq,Jeq ; do i=is,ie + seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.) + enddo ; enddo + + do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then + if ((e(i,j+1,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i,j+1,1))) then + ! This is a typical case in the open ocean, so use the topmost interface. + T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1) + S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1) + p_int_S(i,J) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p_int_N(i,J) = -GxRho*(e(i,j+1,1) - Z_0p(i,j)) + inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1)) + dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1)) + seek_y_cor(i,J) = .false. + endif + endif ; enddo ; enddo + + do k=1,nz + 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 mass-weighting is not + ! activated in the subgrid interpolation. + 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)) <= 0.0)) 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) + ! These pressures are only used for the equation of state, and are only a function of + ! height, consistent with the expressions in the int_density_dz routines. + p_int_S(i,J) = -GxRho*(e(i,j,K+1) - Z_0p(i,j)) + p_int_N(i,J) = -GxRho*(e(i,j+1,K+1) - Z_0p(i,j)) + inty_pa_nonlin(i,J) = inty_pa(i,J,K+1) - 0.5*(pa(i,j,K+1) + pa(i,j+1,K+1)) + dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,K+1)-e(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 + + if (do_more_k) then + ! There are still points where a correction is needed, so use the top interface for lack of a better idea? + do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then + T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1) + S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1) + p_int_S(i,J) = -GxRho*(e(i,j,1) - Z_0p(i,j)) + p_int_N(i,J) = -GxRho*(e(i,j+1,1) - Z_0p(i,j)) + inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1)) + dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1)) + seek_y_cor(i,J) = .false. + endif ; enddo ; enddo + endif + + 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 + 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 + intx_pa(I,j,K) = intx_pa(I,j,K) + intx_pa_cor_ri(I,j) + enddo ; enddo ; enddo + + do K=1,nz+1 ; do J=Jsq,Jeq ; do i=is,ie + inty_pa(i,J,K) = inty_pa(i,J,K) + inty_pa_cor_ri(i,J) + enddo ; enddo ; enddo + endif ! intx_pa and inty_pa have now been reset to reflect the properties of an unimpeded interface. + ! Compute pressure gradient in x direction !$OMP parallel do default(shared) do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -1027,7 +1740,10 @@ 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 ! This include declares and sets the variable "version". # include "version_variable.h" @@ -1043,6 +1759,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 "//& @@ -1062,6 +1781,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 "//& @@ -1072,9 +1801,14 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, "If True, use the ALE algorithm (regridding/remapping). "//& "If False, use the layered isopycnal algorithm.", default=.false. ) call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", useMassWghtInterp, & - "If true, use mass weighting when interpolating T/S for "//& - "integrals near the bathymetry in FV pressure gradient "//& - "calculations.", default=.false.) + "If true, use mass weighting when interpolating T/S for integrals "//& + "near the bathymetry in FV pressure gradient calculations.", & + default=.false.) + call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP", MassWghtInterpTop, & + "If true and MASS_WEIGHT_IN_PRESSURE_GRADIENT is true, use mass weighting when "//& + "interpolating T/S for integrals near the top of the water column in FV "//& + "pressure gradient calculations. ", & + default=.false.) !### Change Default to MASS_WEIGHT_IN_PRESSURE_GRADIENT? call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_NONBOUS_BUG", MassWghtInterp_NonBous_bug, & "If true, use a masking bug in non-Boussinesq calculations with mass weighting "//& "when interpolating T/S for integrals near the bathymetry in FV pressure "//& @@ -1083,8 +1817,25 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp, CS%MassWghtInterp = 0 if (useMassWghtInterp) & CS%MassWghtInterp = ibset(CS%MassWghtInterp, 0) ! Same as CS%MassWghtInterp + 1 + if (MassWghtInterpTop) & + CS%MassWghtInterp = ibset(CS%MassWghtInterp, 1) ! Same as CS%MassWghtInterp + 2 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, & + "If true, use a correction for surface pressure curvature in intx_pa.", & + 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., 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%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, "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.) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 8b820594aa..90994dd073 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -40,7 +40,7 @@ module MOM_density_integrals !! required for calculating the finite-volume form pressure accelerations in a !! Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -77,6 +77,8 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, !! layer divided by the y grid spacing [R L2 T-2 ~> Pa] real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] + real, dimension(SZI_(HI),SZJ_(HI)), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -85,10 +87,10 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, if (EOS_quadrature(EOS)) then call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p=Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p=Z_0p) else call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p=Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif end subroutine int_density_dz @@ -97,7 +99,7 @@ end subroutine int_density_dz !> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, use_inaccurate_form, Z_0p) type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables. real, dimension(SZI_(HI),SZJ_(HI)), & @@ -135,6 +137,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! layer divided by the y grid spacing [R L2 T-2 ~> Pa] real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] + real, dimension(SZI_(HI),SZJ_(HI)), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -169,7 +173,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim] real :: intz(5) ! The gravitational acceleration times the integrals of density ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] - logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: do_massWeight ! Indicates whether to do mass weighting near bathymetry + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation ! of density anomalies. integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state @@ -198,13 +203,16 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form endif - do_massWeight = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (do_massWeight) then - if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - "bathyT must be present if MassWghtInterp is present and true.") - if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - "dz_neglect must be present if MassWghtInterp is present and true.") + do_massWeight = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + if (do_massWeight .and. .not.present(bathyT)) call MOM_error(FATAL, & + "int_density_dz_generic: bathyT must be present if near-bottom mass weighting is in use.") + if (top_massWeight .and. .not.present(SSH)) call MOM_error(FATAL, & + "int_density_dz_generic: SSH must be present if near-surface mass weighting is in use.") + if ((do_massWeight .or. top_massWeight) .and. .not.present(dz_neglect)) call MOM_error(FATAL, & + "int_density_dz_generic: dz_neglect must be present if mass weighting is in use.") endif ! Set the loop ranges for equation of state calculations at various points. @@ -248,6 +256,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect @@ -314,6 +324,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect @@ -470,11 +482,13 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & real :: dz(HI%iscB:HI%iecB+1) ! Layer thicknesses at tracer points [Z ~> m] real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m] real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m] - real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim] + real :: massWeightToggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim] + real :: TopWeightToggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim] real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [C ~> degC] real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [S ~> ppt] real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thickness weight [Z ~> m] + real :: hWghtTop ! An ice draft limited thickness weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] logical :: use_rho_ref ! Pass rho_ref to the equation of state for more accurate calculation @@ -496,8 +510,11 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & else z0pres(:,:) = 0.0 endif - massWeightToggle = 0. - if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif + massWeightToggle = 0. ; TopWeightToggle = 0. + if (present(MassWghtInterp)) then + if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. + if (BTEST(MassWghtInterp, 1)) TopWeightToggle = 1. + endif use_rho_ref = .true. if (present(use_inaccurate_form)) use_rho_ref = .not. use_inaccurate_form @@ -592,6 +609,17 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & ! this distance by the layer thickness to replicate other models. hWght = massWeightToggle * & max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K)) + ! CY: The below code just uses top interface, which may be bad in high res open ocean + ! We want something like if (pa(i+1,k+1) 0.) then hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff @@ -688,6 +716,17 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & ! this distance by the layer thickness to replicate other models. hWght = massWeightToggle * & max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K)) + ! CY: The below code just uses top interface, which may be bad in high res open ocean + ! We want something like if (pa(j+1,k+1) 0.) then hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff @@ -873,7 +912,8 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: dz ! Layer thicknesses at tracer points [Z ~> m] real :: dz_x(5,HI%iscB:HI%iecB) ! Layer thicknesses along an x-line of subgrid locations [Z ~> m] real :: dz_y(5,HI%isc:HI%iec) ! Layer thicknesses along a y-line of subgrid locations [Z ~> m] - real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim] + real :: massWeightToggle ! A non-dimensional toggle factor for near-bottom mass weighting (0 or 1) [nondim] + real :: TopWeightToggle ! A non-dimensional toggle factor for near-surface mass weighting (0 or 1) [nondim] real :: Ttl, Tbl, Tml, Ttr, Tbr, Tmr ! Temperatures at the velocity cell corners [C ~> degC] real :: Stl, Sbl, Sml, Str, Sbr, Smr ! Salinities at the velocity cell corners [S ~> ppt] real :: s6 ! PPM curvature coefficient for S [S ~> ppt] @@ -882,6 +922,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S [S ~> ppt] real :: z0pres(HI%isd:HI%ied,HI%jsd:HI%jed) ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thickness weight [Z ~> m] + real :: hWghtTop ! A surface displacement limited thickness weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state @@ -902,8 +943,11 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & else z0pres(:,:) = 0.0 endif - massWeightToggle = 0. - if (present(MassWghtInterp)) then ; if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. ; endif + massWeightToggle = 0. ; TopWeightToggle = 0. + if (present(MassWghtInterp)) then + if (BTEST(MassWghtInterp, 0)) massWeightToggle = 1. + if (BTEST(MassWghtInterp, 1)) TopWeightToggle = 1. + endif ! In event PPM calculation is bypassed with use_PPM=False s6 = 0. @@ -990,6 +1034,9 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! this distance by the layer thickness to replicate other models. hWght = massWeightToggle * & max(0., -bathyT(i,j)-e(i+1,j,K), -bathyT(i+1,j)-e(i,j,K)) + hWghtTop = TopWeightToggle * & + max(0., e(i+1,j,K+1)-e(i,j,1), e(i,j,K+1)-e(i+1,j,1)) + hWght = max(hWght, hWghtTop) if (hWght > 0.) then hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff hR = (e(i+1,j,K) - e(i+1,j,K+1)) + dz_subroundoff @@ -1095,6 +1142,9 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! this distance by the layer thickness to replicate other models. hWght = massWeightToggle * & max(0., -bathyT(i,j)-e(i,j+1,K), -bathyT(i,j+1)-e(i,j,K)) + hWghtTop = TopWeightToggle * & + max(0., e(i,j+1,K+1)-e(i,j,1), e(i,j,K+1)-e(i,j+1,1)) + hWght = max(hWght, hWghtTop) if (hWght > 0.) then hL = (e(i,j,K) - e(i,j,K+1)) + dz_subroundoff hR = (e(i,j+1,K) - e(i,j+1,K+1)) + dz_subroundoff @@ -1196,7 +1246,7 @@ end subroutine int_density_dz_generic_ppm !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, MassWghtInterp) + bathyP, P_surf, dP_tiny, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1230,6 +1280,8 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_tiny !< A minuscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -1238,11 +1290,11 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & if (EOS_quadrature(EOS)) then call int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, MassWghtInterp) + bathyP, P_surf, dP_tiny, MassWghtInterp) else call analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, MassWghtInterp) + bathyP, P_surf, dP_tiny, MassWghtInterp) endif end subroutine int_specific_vol_dp @@ -1254,7 +1306,7 @@ end subroutine int_specific_vol_dp !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, MassWghtInterp) + bathyP, P_surf, dP_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [C ~> degC] @@ -1289,6 +1341,8 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(SZI_(HI),SZJ_(HI)), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A minuscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -1325,6 +1379,7 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state @@ -1338,14 +1393,17 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set - if (do_massWeight) then - if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& - "bathyP must be present if MassWghtInterp is present and true.") - if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& - "dP_neglect must be present if MassWghtInterp is present and true.") + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + if (do_massWeight .and. .not.present(bathyP)) call MOM_error(FATAL, & + "int_spec_vol_dp_generic_pcm: bathyP must be present if near-bottom mass weighting is in use.") + if (top_massWeight .and. .not.present(P_surf)) call MOM_error(FATAL, & + "int_spec_vol_dp_generic_pcm: P_surf must be present if near-surface mass weighting is in use.") + if ((do_massWeight .or. top_massWeight) .and. .not.present(dP_neglect)) call MOM_error(FATAL, & + "int_spec_vol_dp_generic_pcm: dP_neglect must be present if mass weighting is in use.") endif ! Set the loop ranges for equation of state calculations at various points. @@ -1390,6 +1448,8 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -1448,6 +1508,8 @@ subroutine int_spec_vol_dp_generic_pcm(T, S, p_t, p_b, alpha_ref, HI, EOS, US, d elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect @@ -1503,7 +1565,7 @@ end subroutine int_spec_vol_dp_generic_pcm !! no free assumptions, apart from the use of Boole's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, & dP_neglect, bathyP, HI, EOS, US, dza, & - intp_dza, intx_dza, inty_dza, MassWghtInterp) + intp_dza, intx_dza, inty_dza, P_surf, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T_t !< Potential temperature at the top of the layer [C ~> degC] @@ -1543,6 +1605,8 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, optional, intent(inout) :: inty_dza !< The integral in y of the difference between !! the geopotential anomaly at the top and bottom of the layer divided !! by the y grid spacing [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(HI),SZJ_(HI)), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -1581,6 +1645,7 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, ! 5 sub-column locations [L2 T-2 ~> m2 s-2] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting integer, dimension(2) :: EOSdom_h5 ! The 5-point h-point i-computational domain for the equation of state integer, dimension(2) :: EOSdom_q15 ! The 3x5-point q-point i-computational domain for the equation of state @@ -1589,9 +1654,14 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, Isq = HI%IscB ; Ieq = HI%IecB ; Jsq = HI%JscB ; Jeq = HI%JecB - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + if (top_massWeight .and. .not.present(P_surf)) call MOM_error(FATAL, & + "int_spec_vol_dp_generic_plm: P_surf must be present if near-surface mass weighting is in use.") + endif do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) @@ -1639,6 +1709,8 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -1703,6 +1775,8 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect @@ -1760,71 +1834,87 @@ end subroutine int_spec_vol_dp_generic_plm !> Diagnose the fractional mass weighting in a layer that might be used with a Boussinesq calculation. -subroutine diagnose_mass_weight_Z(z_t, z_b, dz_neglect, bathyT, HI, MassWt_u, MassWt_v) +subroutine diagnose_mass_weight_Z(z_t, z_b, bathyT, SSH, dz_neglect, MassWghtInterp, HI, & + MassWt_u, MassWt_v) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: z_t !< Height at the top of the layer in depth units [Z ~> m] real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: z_b !< Height at the bottom of the layer [Z ~> m] - real, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] - real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: SSH !< The sea surface height [Z ~> m] + real, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] + integer, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, dimension(SZIB_(HI),SZJ_(HI)), & - optional, intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] + intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] real, dimension(SZI_(HI),SZJB_(HI)), & - optional, intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] + intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] ! Local variables real :: hWght ! A pressure-thickness below topography [Z ~> m] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m] real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2] + logical :: do_massWeight ! Indicates whether to do mass weighting near bathymetry + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface integer :: Isq, Ieq, Jsq, Jeq, i, j Isq = HI%IscB ; Ieq = HI%IecB Jsq = HI%JscB ; Jeq = HI%JecB - if (present(MassWt_u)) then - do j=HI%jsc,HI%jec ; do I=Isq,Ieq - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, like thickness weighting. + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + + ! Calculate MassWt_u + do j=HI%jsc,HI%jec ; do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = 0.0 + if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) - if (hWght > 0.) then - 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 ) - MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom - else - MassWt_u(I,j) = 0.0 - endif - enddo ; enddo - endif + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) + if (hWght > 0.) then + 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 ) + MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_u(I,j) = 0.0 + endif + enddo ; enddo - if (present(MassWt_v)) then - do J=Jsq,Jeq ; do i=HI%isc,HI%iec - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, like thickness weighting. + ! Calculate MassWt_v + do J=Jsq,Jeq ; do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = 0.0 + if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) - if (hWght > 0.) then - 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 ) - MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom - else - MassWt_v(i,J) = 0.0 - endif - enddo ; enddo - endif + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) + if (hWght > 0.) then + 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 ) + MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_v(i,J) = 0.0 + endif + enddo ; enddo end subroutine diagnose_mass_weight_Z !> Diagnose the fractional mass weighting in a layer that might be used with a non-Boussinesq calculation. -subroutine diagnose_mass_weight_p(p_t, p_b, dP_neglect, bathyP, HI, MassWt_u, MassWt_v) +subroutine diagnose_mass_weight_p(p_t, p_b, bathyP, P_surf, dP_neglect, MassWghtInterp, HI, & + MassWt_u, MassWt_v) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: p_t !< Pressure atop the layer [R L2 T-2 ~> Pa] @@ -1834,55 +1924,78 @@ subroutine diagnose_mass_weight_p(p_t, p_b, dP_neglect, bathyP, HI, MassWt_u, Ma !! the same units as p_t [R L2 T-2 ~> Pa] real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] + integer, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use + !! mass weighting to interpolate T/S in integrals real, dimension(SZIB_(HI),SZJ_(HI)), & - optional, intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] + intent(inout) :: MassWt_u !< The fractional mass weighting at u-points [nondim] real, dimension(SZI_(HI),SZJB_(HI)), & - optional, intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] + intent(inout) :: MassWt_v !< The fractional mass weighting at v-points [nondim] ! Local variables real :: hWght ! A pressure-thickness below topography [R L2 T-2 ~> Pa] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [R L2 T-2 ~> Pa] real :: iDenom ! The inverse of the denominator in the weights [T4 R-2 L-4 ~> Pa-2] + logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface + logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting integer :: Isq, Ieq, Jsq, Jeq, i, j Isq = HI%IscB ; Ieq = HI%IecB Jsq = HI%JscB ; Jeq = HI%JecB - if (present(MassWt_u)) then - do j=HI%jsc,HI%jec ; do I=Isq,Ieq - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, like thickness weighting. + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + + ! Calculate MassWt_u + do j=HI%jsc,HI%jec ; do I=Isq,Ieq + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = 0.0 + if (do_massWeight .and. massWeight_bug) then + hWght = max(0., bathyP(i,j)-p_t(i+1,j), bathyP(i+1,j)-p_t(i,j)) + elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) - if (hWght > 0.) then - 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 ) - MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom - else - MassWt_u(I,j) = 0.0 - endif - enddo ; enddo - endif + endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) + if (hWght > 0.) then + 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 ) + MassWt_u(I,j) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_u(I,j) = 0.0 + endif + enddo ; enddo - if (present(MassWt_v)) then - do J=Jsq,Jeq ; do i=HI%isc,HI%iec - ! hWght is the distance measure by which the cell is violation of - ! hydrostatic consistency. For large hWght we bias the interpolation - ! of T,S along the top and bottom integrals, like thickness weighting. + ! Calculate MassWt_v + do J=Jsq,Jeq ; do i=HI%isc,HI%iec + ! hWght is the distance measure by which the cell is violation of + ! hydrostatic consistency. For large hWght we bias the interpolation + ! of T,S along the top and bottom integrals, like thickness weighting. + hWght = 0.0 + if (do_massWeight .and. massWeight_bug) then + hWght = max(0., bathyP(i,j)-p_t(i,j+1), bathyP(i,j+1)-p_t(i,j)) + elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) - if (hWght > 0.) then - 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 ) - MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom - else - MassWt_v(i,J) = 0.0 - endif - enddo ; enddo - endif + endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) + if (hWght > 0.) then + 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 ) + MassWt_v(i,J) = (hWght*hR + hWght*hL) * iDenom + else + MassWt_v(i,J) = 0.0 + endif + enddo ; enddo end subroutine diagnose_mass_weight_p diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 74f540f64f..bfab3f5719 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1226,7 +1226,7 @@ end function EOS_domain !! series for log(1-eps/1+eps) that assumes that |eps| < 0.34. subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, MassWghtInterp) + bathyP, P_surf, dP_tiny, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1259,6 +1259,8 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_tiny !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -1280,20 +1282,20 @@ subroutine analytic_int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%kg_m3_to_R*EOS%Rho_T0_S0, & dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, dza, & intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, MassWghtInterp) + bathyP, P_surf, dP_tiny, MassWghtInterp) case (EOS_WRIGHT) call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & + inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case (EOS_WRIGHT_FULL) call int_spec_vol_dp_wright_full(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & + inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case (EOS_WRIGHT_REDUCED) call int_spec_vol_dp_wright_red(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & - inty_dza, halo_size, bathyP, dP_tiny, MassWghtInterp, & + inty_dza, halo_size, bathyP, P_surf, dP_tiny, MassWghtInterp, & SV_scale=EOS%R_to_kg_m3, pres_scale=EOS%RL2_T2_to_Pa, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt) case default @@ -1306,7 +1308,7 @@ end subroutine analytic_int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp, Z_0p) + intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] @@ -1342,6 +1344,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, !! layer divided by the y grid spacing [R L2 T-2 ~> Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -1367,23 +1371,23 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (dRdT_scale /= 1.0) .or. (dRdS_scale /= 1.0)) then call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & rho_scale*EOS%Rho_T0_S0, dRdT_scale*EOS%dRho_dT, dRdS_scale*EOS%dRho_dS, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp) else call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, MassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, MassWghtInterp) endif case (EOS_WRIGHT) rho_scale = EOS%kg_m3_to_R pres_scale = EOS%RL2_T2_to_Pa if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case (EOS_WRIGHT_FULL) @@ -1391,12 +1395,12 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, pres_scale = EOS%RL2_T2_to_Pa if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case (EOS_WRIGHT_REDUCED) @@ -1404,12 +1408,12 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, pres_scale = EOS%RL2_T2_to_Pa if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0) .or. (EOS%C_to_degC /= 1.0) .or. (EOS%S_to_ppt /= 1.0)) then call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, rho_scale, pres_scale, & temp_scale=EOS%C_to_degC, saln_scale=EOS%S_to_ppt, Z_0p=Z_0p) else call int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, & dz_neglect, MassWghtInterp, Z_0p=Z_0p) endif case default diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 11fa57644d..874d3e784e 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -387,7 +387,7 @@ end subroutine EoS_fit_range_buggy_Wright !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, & MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -424,6 +424,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -481,6 +483,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m @@ -531,14 +534,11 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - ! if (do_massWeight) then - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if MassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if MassWghtInterp is present and true.") - ! endif + do_massWeight = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = (a0 + a1s*T(i,j)) + a2s*S(i,j) @@ -571,6 +571,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect @@ -613,6 +615,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect @@ -656,7 +660,7 @@ end subroutine int_density_dz_wright !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & - intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & + intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, & MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -693,6 +697,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & !! dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -743,6 +749,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] @@ -780,15 +787,12 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set -! if (do_massWeight) then -! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if MassWghtInterp is present and true.") -! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if MassWghtInterp is present and true.") -! endif + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + endif ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) do j=jsh,jeh ; do i=ish,ieh @@ -818,6 +822,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -862,6 +868,8 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 6dba8444dd..4be5f2940e 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -393,7 +393,7 @@ end subroutine EoS_fit_range_Wright_full !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, & MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -430,6 +430,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -487,6 +489,7 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m @@ -537,14 +540,11 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - ! if (do_massWeight) then - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if MassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if MassWghtInterp is present and true.") - ! endif + do_massWeight = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = a0 + (a1s*T(i,j) + a2s*S(i,j)) @@ -576,6 +576,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect @@ -618,6 +620,8 @@ subroutine int_density_dz_wright_full(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect @@ -661,7 +665,7 @@ end subroutine int_density_dz_wright_full !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & - intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & + intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, & MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -698,6 +702,8 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & !! dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -749,6 +755,7 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] @@ -786,15 +793,12 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set -! if (do_massWeight) then -! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if MassWghtInterp is present and true.") -! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if MassWghtInterp is present and true.") -! endif + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + endif ! alpha = (lambda + al0*(pressure + p0)) / (pressure + p0) do j=jsh,jeh ; do i=ish,ieh @@ -825,6 +829,8 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -870,6 +876,8 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index 87d6a16dba..1635f9e809 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -395,7 +395,7 @@ end subroutine EoS_fit_range_Wright_red !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, SSH, dz_neglect, & MassWghtInterp, rho_scale, pres_scale, temp_scale, saln_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -432,6 +432,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -489,6 +491,7 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m @@ -539,14 +542,11 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - ! if (do_massWeight) then - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if MassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if MassWghtInterp is present and true.") - ! endif + do_massWeight = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 al0_2d(i,j) = a0 + (a1s*T(i,j) + a2s*S(i,j)) @@ -578,6 +578,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i+1,j) - z_b(i+1,j)) + dz_neglect @@ -620,6 +622,8 @@ subroutine int_density_dz_wright_red(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) if (hWght > 0.) then hL = (z_t(i,j) - z_b(i,j)) + dz_neglect hR = (z_t(i,j+1) - z_b(i,j+1)) + dz_neglect @@ -663,7 +667,7 @@ end subroutine int_density_dz_wright_red !! rule to do the horizontal integrals, and from a truncation in the series for log(1-eps/1+eps) !! that assumes that |eps| < 0.34. subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & - intp_dza, intx_dza, inty_dza, halo_size, bathyP, dP_neglect, & + intp_dza, intx_dza, inty_dza, halo_size, bathyP, P_surf, dP_neglect, & MassWghtInterp, SV_scale, pres_scale, temp_scale, saln_scale) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & @@ -700,6 +704,8 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & !! dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -751,6 +757,7 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & real :: c4s ! Partly rescaled version of c4 [m2 s-2 S-1 ~> m2 s-2 PSU-1] real :: c5s ! Partly rescaled version of c5 [m2 s-2 C-1 S-1 ~> m2 s-2 degC-1 PSU-1] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants [nondim] real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants [nondim] @@ -788,15 +795,12 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & c4s = c4s * saln_scale ; c5s = c5s * saln_scale endif ; endif - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set -! if (do_massWeight) then -! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if MassWghtInterp is present and true.") -! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if MassWghtInterp is present and true.") -! endif + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + endif ! alpha(j) = (lambda + al0*(pressure(j) + p0)) / (pressure(j) + p0) do j=jsh,jeh ; do i=ish,ieh @@ -827,6 +831,8 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i+1,j) - p_t(i+1,j)) + dP_neglect @@ -872,6 +878,8 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght > 0.) then hL = (p_b(i,j) - p_t(i,j)) + dP_neglect hR = (p_b(i,j+1) - p_t(i,j+1)) + dP_neglect diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index f5673ba5f2..e443970535 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -258,7 +258,7 @@ end subroutine set_params_linear !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, MassWghtInterp) + bathyT, SSH, dz_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -299,6 +299,8 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & !! layer divided by the y grid spacing [R L2 T-2 ~> Pa] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyT !< The depth of the bathymetry [Z ~> m]. + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: SSH !< The sea surface height [Z ~> m] real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use !! mass weighting to interpolate T/S in integrals @@ -317,6 +319,7 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & real :: intz(5) ! The integrals of density with height at the ! 5 sub-column locations [R L2 T-2 ~> Pa] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface real, parameter :: C1_6 = 1.0/6.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]. integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j, m @@ -327,14 +330,11 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & is = HI%isc ; ie = HI%iec js = HI%jsc ; je = HI%jec - do_massWeight = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - ! if (do_massWeight) then - ! if (.not.present(bathyT)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "bathyT must be present if MassWghtInterp is present and true.") - ! if (.not.present(dz_neglect)) call MOM_error(FATAL, "int_density_dz_generic: "//& - ! "dz_neglect must be present if MassWghtInterp is present and true.") - ! endif + do_massWeight = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + endif do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 dz = z_t(i,j) - z_b(i,j) @@ -350,6 +350,8 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i+1,j), -bathyT(i+1,j)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i+1,j)-SSH(i,j), z_b(i,j)-SSH(i+1,j)) if (hWght <= 0.0) then dzL = z_t(i,j) - z_b(i,j) ; dzR = z_t(i+1,j) - z_b(i+1,j) @@ -389,6 +391,8 @@ subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HI, & hWght = 0.0 if (do_massWeight) & hWght = max(0., -bathyT(i,j)-z_t(i,j+1), -bathyT(i,j+1)-z_t(i,j)) + if (top_massWeight) & + hWght = max(hWght, z_b(i,j+1)-SSH(i,j), z_b(i,j)-SSH(i,j+1)) if (hWght <= 0.0) then dzL = z_t(i,j) - z_b(i,j) ; dzR = z_t(i,j+1) - z_b(i,j+1) @@ -429,7 +433,7 @@ end subroutine int_density_dz_linear !! model. Specific volume is assumed to vary linearly between adjacent points. subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & dRho_dT, dRho_dS, dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, MassWghtInterp) + bathyP, P_surf, dP_neglect, MassWghtInterp) type(hor_index_type), intent(in) :: HI !< The ocean's horizontal index type. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -469,6 +473,8 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & integer, optional, intent(in) :: halo_size !< The width of halo points on which to calculate dza. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(in) :: bathyP !< The pressure at the bathymetry [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & + optional, intent(in) :: P_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] real, optional, intent(in) :: dP_neglect !< A miniscule pressure change with !! the same units as p_t [R L2 T-2 ~> Pa] integer, optional, intent(in) :: MassWghtInterp !< A flag indicating whether and how to use @@ -488,6 +494,7 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2] logical :: do_massWeight ! Indicates whether to do mass weighting. + logical :: top_massWeight ! Indicates whether to do mass weighting the sea surface logical :: massWeight_bug ! If true, use an incorrect expression to determine where to apply mass weighting real, parameter :: C1_6 = 1.0/6.0, C1_90 = 1.0/90.0 ! Rational constants [nondim]. integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, halo @@ -498,15 +505,12 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh) ; endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh) ; endif - do_massWeight = .false. ; massWeight_bug = .false. - if (present(MassWghtInterp)) do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values - if (present(MassWghtInterp)) massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set -! if (do_massWeight) then -! if (.not.present(bathyP)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "bathyP must be present if MassWghtInterp is present and true.") -! if (.not.present(dP_neglect)) call MOM_error(FATAL, "int_spec_vol_dp_generic: "//& -! "dP_neglect must be present if MassWghtInterp is present and true.") -! endif + do_massWeight = .false. ; massWeight_bug = .false. ; top_massWeight = .false. + if (present(MassWghtInterp)) then + do_massWeight = BTEST(MassWghtInterp, 0) ! True for odd values + top_massWeight = BTEST(MassWghtInterp, 1) ! True if the 2 bit is set + massWeight_bug = BTEST(MassWghtInterp, 3) ! True if the 8 bit is set + endif do j=jsh,jeh ; do i=ish,ieh dp = p_b(i,j) - p_t(i,j) @@ -527,6 +531,8 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & elseif (do_massWeight) then hWght = max(0., p_t(i+1,j)-bathyP(i,j), p_t(i,j)-bathyP(i+1,j)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i+1,j), P_surf(i+1,j)-p_b(i,j)) if (hWght <= 0.0) then dpL = p_b(i,j) - p_t(i,j) ; dpR = p_b(i+1,j) - p_t(i+1,j) @@ -575,6 +581,8 @@ subroutine int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, Rho_T0_S0, & elseif (do_massWeight) then hWght = max(0., p_t(i,j+1)-bathyP(i,j), p_t(i,j)-bathyP(i,j+1)) endif + if (top_massWeight) & + hWght = max(hWght, P_surf(i,j)-p_b(i,j+1), P_surf(i,j+1)-p_b(i,j)) if (hWght <= 0.0) then dpL = p_b(i,j) - p_t(i,j) ; dpR = p_b(i,j+1) - p_t(i,j+1)