diff --git a/.github/actions/ubuntu-setup/action.yml b/.github/actions/ubuntu-setup/action.yml index 83d6795954..0f53a68c70 100644 --- a/.github/actions/ubuntu-setup/action.yml +++ b/.github/actions/ubuntu-setup/action.yml @@ -23,7 +23,7 @@ runs: run: | echo "::group::config.mk" cd .testing - echo "FCFLAGS_DEBUG = -g -O0 -Wextra -Wno-compare-reals -fbacktrace -ffpe-trap=invalid,zero,overflow -fcheck=bounds" >> config.mk + echo "FCFLAGS_DEBUG = -g -O0 -std=f2018 -Wextra -Wno-compare-reals -fbacktrace -ffpe-trap=invalid,zero,overflow -fcheck=bounds" >> config.mk echo "FCFLAGS_REPRO = -g -O2 -fbacktrace" >> config.mk echo "FCFLAGS_INIT = -finit-real=snan -finit-integer=2147483647 -finit-derived" >> config.mk echo "FCFLAGS_FMS = -g -fbacktrace -O0" >> config.mk diff --git a/.testing/Makefile b/.testing/Makefile index 085fea2655..adc1a20a1e 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -83,19 +83,21 @@ export FMS_URL # TODO: This needs more automated configuration MPIRUN ?= mpirun -# Generic compiler variables are pass through to the builds +# Generic compiler variables are passed through to the builds export CC export MPICC export FC export MPIFC # Builds are distinguished by FCFLAGS -FCFLAGS_DEBUG ?= -g -O0 +FCFLAGS ?= -g -O0 + +FCFLAGS_DEBUG ?= $(FCFLAGS) FCFLAGS_REPRO ?= -g -O2 FCFLAGS_OPT ?= -g -O3 -mavx -fno-omit-frame-pointer FCFLAGS_INIT ?= FCFLAGS_COVERAGE ?= -g -O0 -fbacktrace --coverage -FCFLAGS_FMS ?= $(FCFLAGS_DEBUG) +FCFLAGS_FMS ?= $(FCFLAGS) # Additional notes: # - These default values are simple, minimalist flags, supported by nearly all # compilers, and are somewhat analogous to GFDL's DEBUG and REPRO builds. 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/config_src/drivers/FMS_cap/ocean_model_MOM.F90 b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 index 45c14e73eb..9c4359bf60 100644 --- a/config_src/drivers/FMS_cap/ocean_model_MOM.F90 +++ b/config_src/drivers/FMS_cap/ocean_model_MOM.F90 @@ -1065,10 +1065,10 @@ subroutine ocean_model_data2D_get(OS, Ocean, name, array2D, isc, jsc) type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly !! visible ocean surface fields. character(len=*) , intent(in) :: name !< The name of the field to extract - real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must - !! cover only the computational domain [various] integer , intent(in) :: isc !< The starting i-index of array2D integer , intent(in) :: jsc !< The starting j-index of array2D + real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must + !! cover only the computational domain [various] integer :: g_isc, g_iec, g_jsc, g_jec, g_isd, g_ied, g_jsd, g_jed, i, j @@ -1188,10 +1188,10 @@ subroutine ocean_model_get_UV_surf(OS, Ocean, name, array2D, isc, jsc) type(ocean_public_type), intent(in) :: Ocean !< A structure containing various publicly !! visible ocean surface fields. character(len=*) , intent(in) :: name !< The name of the current (ua or va) to extract - real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must - !! cover only the computational domain [L T-1 ~> m s-1] integer , intent(in) :: isc !< The starting i-index of array2D integer , intent(in) :: jsc !< The starting j-index of array2D + real, dimension(isc:,jsc:), intent(out):: array2D !< The values of the named field, it must + !! cover only the computational domain [L T-1 ~> m s-1] type(ocean_grid_type) , pointer :: G !< The ocean's grid structure type(surface), pointer :: sfc_state !< A structure containing fields that 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_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 46d5a99666..5fbf12a0d0 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -2778,7 +2778,7 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "VISC_REM_CONT_HVEL_FIX", CS%visc_rem_hvel_fix, & "If true, velocity cell thickness h_[uv] from the continuity solver "//& "is not multiplied by visc_rem_[uv]. Default of this flag is set by "//& - "VISC_REM_BUG.", default=.not.visc_rem_bug) + "VISC_REM_BUG.", default=.False.) !, default=.not.visc_rem_bug) CS%diag => diag id_clock_reconstruct = cpu_clock_id('(Ocean continuity reconstruction)', grain=CLOCK_ROUTINE) 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..938634c1ea 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 @@ -2424,7 +2428,7 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & tol_here = 0.5*tol*(abs(SpV_avg_a(1)) + abs(SpV_avg_q(1))) test_OK = (abs(SpV_avg_a(1) - SpV_avg_q(1)) < tol_here) if (verbose) then - write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2"), tol=",ES16.8)') & + write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2,"), tol=",ES16.8)') & SpV_avg_a(1), SpV_avg_q(1), SpV_avg_a(1) - SpV_avg_q(1), & 2.0*(SpV_avg_a(1) - SpV_avg_q(1)) / (abs(SpV_avg_a(1)) + abs(SpV_avg_q(1)) + tiny(SpV_avg_a(1))), & tol_here @@ -2504,8 +2508,7 @@ logical function check_FD(val, val_fd, tol, verbose, field_name, order) check_FD = ( abs(val_fd(1) - val) < (1.2*abs(val_fd(2) - val)/2**order + abs(tol)) ) - ! write(mesg, '(ES16.8," and ",ES16.8," differ by ",ES16.8," (",ES10.2"), tol=",ES16.8)') & - write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2"), tol=",ES16.8)') & + write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2,"), tol=",ES16.8)') & val, val_fd(1), val - val_fd(1), & 2.0*(val - val_fd(1)) / (abs(val) + abs(val_fd(1)) + tiny(val)), & (1.2*abs(val_fd(2) - val)/2**order + abs(tol)) 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) diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index d937ed7b0c..81e4425be3 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -680,10 +680,10 @@ subroutine write_auto_mask_file(mask_table, layout, npes, filename) true_num_masked_blocks = layout(1) * layout(2) - npes call open_ASCII_file(file_ascii, trim(filename), action=WRITEONLY_FILE) - write(file_ascii, '(I0)'), true_num_masked_blocks - write(file_ascii, '(I0,",",I0)'), layout(1), layout(2) + write(file_ascii, '(I0)') true_num_masked_blocks + write(file_ascii, '(I0,",",I0)') layout(1), layout(2) do p = 1, true_num_masked_blocks - write(file_ascii, '(I0,",",I0)'), mask_table(p,1), mask_table(p,2) + write(file_ascii, '(I0,",",I0)') mask_table(p,1), mask_table(p,2) enddo call close_file(file_ascii) end subroutine write_auto_mask_file diff --git a/src/framework/posix.F90 b/src/framework/posix.F90 index 1087958939..a9829c510e 100644 --- a/src/framework/posix.F90 +++ b/src/framework/posix.F90 @@ -437,6 +437,7 @@ function setjmp_missing(env) result(rc) bind(c) error stop ! NOTE: compilers may expect a return value, even if it is unreachable + read env%state rc = -1 end function setjmp_missing @@ -450,6 +451,9 @@ subroutine longjmp_missing(env, val) bind(c) print '(a)', 'ERROR: longjmp() is not implemented in this build.' print '(a)', 'Recompile with autoconf or -DLONGJMP_NAME=\"\".' error stop + + read env%state + read char(val) end subroutine longjmp_missing !> Placeholder function for a missing or unconfigured sigsetjmp @@ -466,6 +470,8 @@ function sigsetjmp_missing(env, savesigs) result(rc) bind(c) error stop ! NOTE: compilers may expect a return value, even if it is unreachable + read env%state + read char(savesigs) rc = -1 end function sigsetjmp_missing @@ -478,6 +484,8 @@ subroutine siglongjmp_missing(env, val) bind(c) print '(a)', 'ERROR: siglongjmp() is not implemented in this build.' print '(a)', 'Recompile with autoconf or -DSIGLONGJMP_NAME=\"\".' + read env%state + read char(val) error stop end subroutine siglongjmp_missing diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 4f811cac87..52b1ebdcea 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -181,6 +181,8 @@ module MOM_ice_shelf !! fluxes. It will avoid large increase in sea level. logical :: constant_sea_level_misomip !< If true, constant_sea_level fluxes are applied only over !! the surface sponge cells from the ISOMIP/MISOMIP configuration + logical :: smb_diag !< If true, calculate diagnostics related to surface mass balance + logical :: bmb_diag !< If true, calculate diagnostics related to basal mass balance real :: min_ocean_mass_float !< The minimum ocean mass per unit area before the ice !! shelf is considered to float when constant_sea_level !! is used [R Z ~> kg m-2] @@ -213,7 +215,17 @@ module MOM_ice_shelf id_bdott_melt = -1, id_bdott_accum = -1, id_bdott = -1, & id_dvafdt = -1, id_g_adot = -1, id_f_adot = -1, id_adot = -1, & id_bdot_melt = -1, id_bdot_accum = -1, id_bdot = -1, & - id_t_area = -1, id_g_area = -1, id_f_area = -1 + id_t_area = -1, id_g_area = -1, id_f_area = -1, & + id_Ant_vaf = -1, id_Ant_g_adott = -1, id_Ant_f_adott = -1, id_Ant_adott = -1, & + id_Ant_bdott_melt = -1, id_Ant_bdott_accum = -1, id_Ant_bdott = -1, & + id_Ant_dvafdt = -1, id_Ant_g_adot = -1, id_Ant_f_adot = -1, id_Ant_adot = -1, & + id_Ant_bdot_melt = -1, id_Ant_bdot_accum = -1, id_Ant_bdot = -1, & + id_Ant_t_area = -1, id_Ant_g_area = -1, id_Ant_f_area = -1, & + id_Gr_vaf = -1, id_Gr_g_adott = -1, id_Gr_f_adott = -1, id_Gr_adott = -1, & + id_Gr_bdott_melt = -1, id_Gr_bdott_accum = -1, id_Gr_bdott = -1, & + id_Gr_dvafdt = -1, id_Gr_g_adot = -1, id_Gr_f_adot = -1, id_Gr_adot = -1, & + id_Gr_bdot_melt = -1, id_Gr_bdot_accum = -1, id_Gr_bdot = -1, & + id_Gr_t_area = -1, id_Gr_g_area = -1, id_Gr_f_area = -1 !>@} type(external_field) :: mass_handle @@ -270,12 +282,10 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) p_int !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & - exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1] - exch_vel_s, & !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1] - tmp, & !< Temporary field used when calculating diagnostics [various] - dh_bdott, & !< Basal melt/accumulation over a time step, used for diagnostics [Z ~> m] - dh_adott !< Surface melt/accumulation over a time step, used for diagnostics [Z ~> m] - + exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1] + exch_vel_s, & !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1] + dh_bdott, & !< Basal melt/accumulation over a time step, used for diagnostics [Z ~> m] + dh_adott !< Surface melt/accumulation over a time step, used for diagnostics [Z ~> m] real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: & mass_flux !< Total mass flux of freshwater across the ice-ocean interface. [R Z L2 T-1 ~> kg s-1] real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: & @@ -343,9 +353,8 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) character(len=160) :: mesg ! The text of an error message integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, is, ie, js, je, ied, jed, it1, it3 - real :: vaf0, vaf ! The previous and current volume above floatation [m3] - logical :: smb_diag=.false., bmb_diag=.false. ! Flags to calculate diagnostics related to surface/basal mass balance - real :: val ! Temporary value when calculating scalar diagnostics [various] + real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation [m3] + !for all ice sheets, Antarctica only, or Greenland only [m3] if (.not. associated(CS)) call MOM_error(FATAL, "shelf_calc_flux: "// & "initialize_ice_shelf must be called before shelf_calc_flux.") @@ -356,13 +365,14 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) time_step = time_step_in Itime_step = 1./time_step - if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & - CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0 ) smb_diag=.true. - if (CS%id_bdott>0 .or. CS%id_bdott_melt>0 .or. CS%id_bdott_accum>0 .or. & - CS%id_bdot >0 .or. CS%id_bdot_melt >0 .or. CS%id_bdot_accum >0) bmb_diag=.true. + dh_adott(:,:)=0.0; dh_bdott(:,:)=0.0 - if (CS%active_shelf_dynamics .and. CS%id_dvafdt > 0) & !calculate previous volume above floatation - call volume_above_floatation(CS%dCS, G, ISS, vaf0) + if (CS%active_shelf_dynamics) then + !calculate previous volumes above floatation + if (CS%id_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0) !all ice sheet + if (CS%id_Ant_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_A, hemisphere=0) !Antarctica only + if (CS%id_Gr_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_G, hemisphere=1) !Greenland only + endif is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; ied = G%ied ; jed = G%jed if (CS%data_override_shelf_fluxes .and. CS%active_shelf_dynamics) then @@ -766,9 +776,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ! Melting has been computed, now is time to update thickness and mass if ( CS%override_shelf_movement .and. (.not.CS%mass_from_file)) then - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_melt(ISS, G, US, time_step, fluxes, CS%density_ice, CS%debug) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -782,9 +792,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) ISS%dhdt_shelf(:,:) = ISS%h_shelf(:,:) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_melt(ISS, G, US, time_step, fluxes, CS%density_ice, CS%debug) - if (bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) + if (CS%bmb_diag) dh_bdott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using melt", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -792,9 +802,9 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) unscale=US%RZ_to_kg_m2) endif - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time) - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je) if (CS%debug) then call hchksum(ISS%h_shelf, "h_shelf after change thickness using surf acc", G%HI, haloshift=0, unscale=US%Z_to_m) @@ -846,69 +856,7 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf, ISS%h_shelf, CS%diag) if (CS%id_dhdt_shelf > 0) call post_data(CS%id_dhdt_shelf, ISS%dhdt_shelf, CS%diag) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask,ISS%hmask,CS%diag) - !scalars - if (CS%active_shelf_dynamics) then - if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) - call volume_above_floatation(CS%dCS, G, ISS, vaf) - if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf - if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Itime_step,CS%diag) !d(vaf)/dt - if (CS%id_adott > 0 .or. CS%id_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val) - if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) - if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Itime_step,CS%diag) - endif - if (CS%id_g_adott > 0 .or. CS%id_g_adot > 0) then !grounded only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) - if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Itime_step,CS%diag) - endif - if (CS%id_f_adott > 0 .or. CS%id_f_adot > 0) then !floating only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott,tmp) - tmp(:,:) = dh_adott(:,:) - tmp(:,:) - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) - if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) - endif - endif - if (CS%id_bdott > 0 .or. CS%id_bdot > 0) then !bottom accumulation - bottom melt - call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val) - if (CS%id_bdott > 0) call post_scalar_data(CS%id_bdott,val ,CS%diag) - if (CS%id_bdot > 0) call post_scalar_data(CS%id_bdot ,val*Itime_step,CS%diag) - endif - if (CS%id_bdott_melt > 0 .or. CS%id_bdot_melt > 0) then !bottom melt - tmp(:,:)=0.0 - do j=js,je ; do i=is,ie - if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) - enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_bdott_melt > 0) call post_scalar_data(CS%id_bdott_melt,val ,CS%diag) - if (CS%id_bdot_melt > 0) call post_scalar_data(CS%id_bdot_melt ,val*Itime_step,CS%diag) - endif - if (CS%id_bdott_accum > 0 .or. CS%id_bdot_accum > 0) then !bottom accumulation - tmp(:,:)=0.0 - do j=js,je ; do i=is,ie - if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) - enddo; enddo - call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_bdott_accum > 0) call post_scalar_data(CS%id_bdott_accum,val ,CS%diag) - if (CS%id_bdot_accum > 0) call post_scalar_data(CS%id_bdot_accum ,val*Itime_step,CS%diag) - endif - if (CS%id_t_area > 0) then - tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) - call post_scalar_data(CS%id_t_area,val,CS%diag) - endif - if (CS%id_g_area > 0 .or. CS%id_f_area > 0) then - tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) - if (CS%id_g_area > 0) then - call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) - call post_scalar_data(CS%id_g_area,val,CS%diag) - endif - if (CS%id_f_area > 0) then - call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val) - call post_scalar_data(CS%id_f_area,val,CS%diag) - endif - endif + call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) call disable_averaging(CS%diag) call cpu_clock_end(id_clock_shelf) @@ -926,20 +874,43 @@ subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS) end subroutine shelf_calc_flux -subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out) +subroutine integrate_over_ice_sheet_area(G, ISS, var, var_scale, var_out, hemisphere) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe the ice-shelf state real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< Ice variable to integrate in arbitrary units [A ~> a] real, intent(in) :: var_scale !< Dimensional scaling for variable to integrate [a A-1 ~> 1] real, intent(out) :: var_out !< Variable integrated over the area of the ice sheet in arbitrary units [a m2] + integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets + integer :: IS_ID ! local copy of hemisphere real, dimension(SZI_(G),SZJ_(G)) :: var_cell !< Variable integrated over the ice-sheet area of each cell !! in arbitrary units [a m2] + integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: i,j + if (present(hemisphere)) then + IS_ID=hemisphere + else + IS_ID=-1 + endif + + mask(:,:)=0 + if (IS_ID==0) then !Antarctica (S. Hemisphere) only + do j = G%jsc,G%jec; do i = G%isc,G%iec + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)<=0.0) mask(i,j)=1 + enddo; enddo + elseif (IS_ID==1) then !Greenland (N. Hemisphere) only + do j = G%jsc,G%jec; do i = G%isc,G%iec + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)>0.0) mask(i,j)=1 + enddo; enddo + else !All ice sheets + mask(G%isc:G%iec,G%jsc:G%jec)=ISS%hmask(G%isc:G%iec,G%jsc:G%jec) + endif + var_cell(:,:)=0.0 do j = G%jsc,G%jec; do i = G%isc,G%iec - if (ISS%hmask(i,j)>0) var_cell(i,j) = (var(i,j) * var_scale) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) + if (mask(i,j)>0) var_cell(i,j) = (var(i,j) * var_scale) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) enddo; enddo + var_out = reproducing_sum(var_cell) end subroutine integrate_over_ice_sheet_area @@ -2031,11 +2002,12 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, 'ice shelf surface mass flux deposition from atmosphere', & 'kg m-2 s-1', conversion=US%RZ_T_to_kg_m2s) endif - !scalars (area integrated) + + !scalars (area integrated over all ice sheets) CS%id_vaf = register_scalar_field('ice_shelf_model', 'int_vaf', CS%diag%axesT1, CS%Time, & 'Area integrated ice sheet volume above floatation', 'm3') CS%id_adott = register_scalar_field('ice_shelf_model', 'int_a', CS%diag%axesT1, CS%Time, & - 'Area integrated (entire ice sheet) change in ice-sheet thickness ' //& + 'Area integrated change in ice-sheet thickness ' //& 'due to surface accum+melt during a DT_THERM time step', 'm3') CS%id_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground', CS%diag%axesT1, CS%Time, & 'Area integrated change in grounded ice-sheet thickness ' //& @@ -2051,16 +2023,16 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, CS%id_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum', CS%diag%axesT1, CS%Time, & 'Area integrated basal accumulation over ice shelves during a DT_THERM a time step', 'm3') CS%id_t_area = register_scalar_field('ice_shelf_model', 'tot_area', CS%diag%axesT1, CS%Time, & - 'Total area of entire ice-sheet', 'm2') + 'Total ice-sheet area', 'm2') CS%id_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float', CS%diag%axesT1, CS%Time, & 'Total area of floating ice shelves', 'm2') CS%id_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground', CS%diag%axesT1, CS%Time, & - 'Total area of grounded ice sheet', 'm2') - !scalars (area integrated rates) + 'Total area of grounded ice sheets', 'm2') + !scalars (area integrated rates over all ice sheets) CS%id_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot', CS%diag%axesT1, CS%Time, & 'Area integrated rate of change in ice-sheet volume above floatation', 'm3 s-1') CS%id_adot = register_scalar_field('ice_shelf_model', 'int_adot', CS%diag%axesT1, CS%Time, & - 'Area integrated (full ice sheet) rate of change in ice-sheet thickness due to surface accum+melt', 'm3 s-1') + 'Area integrated rate of change in ice-sheet thickness due to surface accum+melt', 'm3 s-1') CS%id_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground', CS%diag%axesT1, CS%Time, & 'Area integrated rate of change in grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') CS%id_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float', CS%diag%axesT1, CS%Time, & @@ -2072,6 +2044,111 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, CS%id_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum', CS%diag%axesT1, CS%Time, & 'Area integrated basal accumulation rate over ice shelves', 'm3 s-1') + !scalars (area integrated over the Antarctic ice sheet) + CS%id_Ant_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_A', CS%diag%axesT1, CS%Time, & + 'Area integrated Antarctic ice sheet volume above floatation', 'm3') + CS%id_Ant_adott = register_scalar_field('ice_shelf_model', 'int_a_A', CS%diag%axesT1, CS%Time, & + 'Area integrated (Antarctic ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_bdott = register_scalar_field('ice_shelf_model', 'int_b_A', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Antarctic floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3') + CS%id_Ant_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt over Antarctic ice shelves during a DT_THERM time step', 'm3') + CS%id_Ant_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation over Antarctic ice shelves during a DT_THERM a time step', 'm3') + CS%id_Ant_t_area = register_scalar_field('ice_shelf_model', 'tot_area_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic ice sheet', 'm2') + CS%id_Ant_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic floating ice shelves', 'm2') + CS%id_Ant_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_A', CS%diag%axesT1, CS%Time, & + 'Total area of Antarctic grounded ice sheet', 'm2') + !scalars (area integrated rates over the Antarctic ice sheet) + CS%id_Ant_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-sheet volume above floatation', 'm3 s-1') + CS%id_Ant_adot = register_scalar_field('ice_shelf_model', 'int_adot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + CS%id_Ant_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_A', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Antarctic ice-shelf thickness due to basal accum+melt', 'm3 s-1') + CS%id_Ant_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt rate over Antarctic ice shelves', 'm3 s-1') + CS%id_Ant_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_A', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation rate over Antarctic ice shelves', 'm3 s-1') + + !scalars (area integrated over the Greenland ice sheet) + CS%id_Gr_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_G', CS%diag%axesT1, CS%Time, & + 'Area integrated Greenland ice sheet volume above floatation', 'm3') + CS%id_Gr_adott = register_scalar_field('ice_shelf_model', 'int_a_G', CS%diag%axesT1, CS%Time, & + 'Area integrated (Greenland ice sheet) change in ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland grounded ice-sheet thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland floating ice-shelf thickness ' //& + 'due to surface accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_bdott = register_scalar_field('ice_shelf_model', 'int_b_G', CS%diag%axesT1, CS%Time, & + 'Area integrated change in Greenland floating ice-shelf thickness '//& + 'due to basal accum+melt during a DT_THERM time step', 'm3') + CS%id_Gr_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt over Greenland ice shelves during a DT_THERM time step', 'm3') + CS%id_Gr_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation over Greenland ice shelves during a DT_THERM a time step', 'm3') + CS%id_Gr_t_area = register_scalar_field('ice_shelf_model', 'tot_area_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland ice sheet', 'm2') + CS%id_Gr_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland floating ice shelves', 'm2') + CS%id_Gr_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_G', CS%diag%axesT1, CS%Time, & + 'Total area of Greenland grounded ice sheet', 'm2') + !scalars (area integrated rates over the Greenland ice sheet) + CS%id_Gr_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-sheet volume above floatation', 'm3 s-1') + CS%id_Gr_adot = register_scalar_field('ice_shelf_model', 'int_adot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland grounded ice-sheet thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland floating ice-shelf thickness due to surface accum+melt', 'm3 s-1') + CS%id_Gr_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_G', CS%diag%axesT1, CS%Time, & + 'Area integrated rate of change in Greenland ice-shelf thickness due to basal accum+melt', 'm3 s-1') + CS%id_Gr_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal melt rate over Greenland ice shelves', 'm3 s-1') + CS%id_Gr_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_G', CS%diag%axesT1, CS%Time, & + 'Area integrated basal accumulation rate over Greenland ice shelves', 'm3 s-1') + + !Flags to calculate diagnostics related to surface/basal mass balance + if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & + CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0 .or. & + CS%id_Ant_adott>0 .or. CS%id_Ant_g_adott>0 .or. CS%id_Ant_f_adott>0 .or. & + CS%id_Ant_adot >0 .or. CS%id_Ant_g_adot >0 .or. CS%id_Ant_f_adot >0 .or. & + CS%id_Gr_adott>0 .or. CS%id_Gr_g_adott>0 .or. CS%id_Gr_f_adott>0 .or. & + CS%id_Gr_adot >0 .or. CS%id_Gr_g_adot >0 .or. CS%id_Gr_f_adot >0) then + CS%smb_diag=.true. + else + CS%smb_diag=.false. + endif + + if (CS%id_bdott>0 .or. CS%id_bdott_melt>0 .or. CS%id_bdott_accum>0 .or. & + CS%id_bdot >0 .or. CS%id_bdot_melt >0 .or. CS%id_bdot_accum >0 .or. & + CS%id_Ant_bdott>0 .or. CS%id_Ant_bdott_melt>0 .or. CS%id_Ant_bdott_accum>0 .or. & + CS%id_Ant_bdot >0 .or. CS%id_Ant_bdot_melt >0 .or. CS%id_Ant_bdot_accum >0 .or. & + CS%id_Gr_bdott>0 .or. CS%id_Gr_bdott_melt>0 .or. CS%id_Gr_bdott_accum>0 .or. & + CS%id_Gr_bdot >0 .or. CS%id_Gr_bdot_melt >0 .or. CS%id_Gr_bdot_accum >0) then + CS%bmb_diag=.true. + else + CS%bmb_diag=.false. + endif + call MOM_IS_diag_mediator_close_registration(CS%diag) if (present(fluxes_in)) call initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in) @@ -2447,11 +2524,9 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in logical :: coupled_GL ! If true the grounding line position is determined based on ! coupled ice-ocean dynamics. integer :: is, ie, js, je, i, j - real :: vaf0, vaf ! The previous and current volume above floatation [m3] - logical :: smb_diag=.false. ! Flags to calculate diagnostics related to surface/basal mass balance - real :: val ! Temporary value when calculating scalar diagnostics [various] + real :: vaf0, vaf0_A, vaf0_G !The previous volumes above floatation + !for all ice sheets, Antarctica only, or Greenland only [m3] real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & - tmp, & ! Temporary field used when calculating diagnostics [various] dh_adott_sum, & ! Surface melt/accumulation over a full time step, used for diagnostics [Z ~> m] dh_adott ! Surface melt/accumulation over a partial time step, used for diagnostics [Z ~> m] @@ -2475,14 +2550,14 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in ISS%dhdt_shelf(:,:) = ISS%h_shelf(:,:) - if (CS%id_adott>0 .or. CS%id_g_adott>0 .or. CS%id_f_adott>0 .or. & - CS%id_adot >0 .or. CS%id_g_adot >0 .or. CS%id_f_adot >0) then - smb_diag=.true. - dh_adott(:,:) = 0.0 ; dh_adott_sum(:,:) = 0.0 ; tmp(:,:) = 0.0 - endif + dh_adott(:,:)=0.0 + + if (CS%smb_diag) dh_adott_sum(:,:) = 0.0 - if (CS%id_dvafdt > 0) & !calculate previous volume above floatation - call volume_above_floatation(CS%dCS, G, ISS, vaf0) + !calculate previous volumes above floatation + if (CS%id_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0) !all ice sheet + if (CS%id_Ant_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_A, hemisphere=0) !Antarctica only + if (CS%id_Gr_dvafdt > 0) call volume_above_floatation(CS%dCS, G, ISS, vaf0_G, hemisphere=1) !Greenland only do while (remaining_time > 0.0) nsteps = nsteps+1 @@ -2497,9 +2572,9 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in call MOM_mesg("solo_step_ice_shelf: "//mesg, 5) endif - if (smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) + if (CS%smb_diag) dh_adott(is:ie,js:je) = ISS%h_shelf(is:ie,js:je) call change_thickness_using_precip(CS, ISS, G, US, fluxes_in, time_step, Time) - if (smb_diag) dh_adott_sum(is:ie,js:je) = dh_adott_sum(is:ie,js:je) + & + if (CS%smb_diag) dh_adott_sum(is:ie,js:je) = dh_adott_sum(is:ie,js:je) + & (ISS%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je)) remaining_time = remaining_time - time_step @@ -2525,47 +2600,222 @@ subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in if (CS%id_h_shelf > 0) call post_data(CS%id_h_shelf ,ISS%h_shelf ,CS%diag) if (CS%id_dhdt_shelf > 0) call post_data(CS%id_dhdt_shelf ,ISS%dhdt_shelf ,CS%diag) if (CS%id_h_mask > 0) call post_data(CS%id_h_mask ,ISS%hmask ,CS%diag) - if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) + call process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Ifull_time_step, dh_adott, dh_adott*0.0) + call disable_averaging(CS%diag) + + call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) +end subroutine solo_step_ice_shelf + +!> Post_data calls for ice-sheet scalars +subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott) + type(ice_shelf_CS), pointer :: CS !< A pointer to the ice shelf control structure + real :: vaf0 !< The previous volumes above floatation for all ice sheets [m3] + real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [m3] + real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [m3] + real :: Itime_step !< Inverse of the time step [T-1 ~> s-1] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_adott !< Surface (plus basal if solo shelf mode) + !! melt/accumulation over a time step [Z ~> m] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_bdott !< Surface (plus basal if solo shelf mode) + !! melt/accumulation over a time step [Z ~> m] + real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: tmp ! Temporary field used when calculating diagnostics [various] + real :: vaf ! The current ice-sheet volume above floatation [m3] + real :: val ! Temporary value when calculating scalar diagnostics [various] + type(ocean_grid_type), pointer :: G => NULL() ! A pointer to the ocean's grid structure + type(unit_scale_type), pointer :: US => NULL() ! Pointer to a structure containing various unit conversion factors + type(ice_shelf_state), pointer :: ISS => NULL() ! A structure with elements that describe the ice-shelf state + integer :: is, ie, js, je, i, j + + G => CS%grid + US => CS%US + ISS => CS%ISS + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + !---ALL ICE SHEET---! + if (CS%id_vaf > 0 .or. CS%id_dvafdt > 0) & !calculate current volume above floatation (vaf) call volume_above_floatation(CS%dCS, G, ISS, vaf) - if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf - if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Ifull_time_step,CS%diag) !d(vaf)/dt + if (CS%id_vaf > 0) call post_scalar_data(CS%id_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_dvafdt > 0) call post_scalar_data(CS%id_dvafdt,(vaf-vaf0)*Itime_step,CS%diag) !d(vaf)/dt if (CS%id_adott > 0 .or. CS%id_adot > 0) then !surface accumulation - surface melt - call integrate_over_ice_sheet_area(G, ISS, dh_adott_sum, US%Z_to_m, val) - if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) - if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Ifull_time_step,CS%diag) + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val) + if (CS%id_adott > 0) call post_scalar_data(CS%id_adott,val ,CS%diag) + if (CS%id_adot > 0) call post_scalar_data(CS%id_adot ,val*Itime_step,CS%diag) endif if (CS%id_g_adott > 0 .or. CS%id_g_adot > 0) then !grounded only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott_sum,tmp) + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) - if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Ifull_time_step,CS%diag) + if (CS%id_g_adott > 0) call post_scalar_data(CS%id_g_adott,val ,CS%diag) + if (CS%id_g_adot > 0) call post_scalar_data(CS%id_g_adot ,val*Itime_step,CS%diag) endif if (CS%id_f_adott > 0 .or. CS%id_f_adot > 0) then !floating only: surface accumulation - surface melt - call masked_var_grounded(G,CS%dCS,dh_adott_sum,tmp) - tmp(:,:) = dh_adott_sum(:,:) - tmp(:,:) + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) - if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) - if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Ifull_time_step,CS%diag) + if (CS%id_f_adott > 0) call post_scalar_data(CS%id_f_adott,val ,CS%diag) + if (CS%id_f_adot > 0) call post_scalar_data(CS%id_f_adot ,val*Itime_step,CS%diag) endif - if (CS%id_t_area > 0) then + if (CS%id_bdott > 0 .or. CS%id_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val) + if (CS%id_bdott > 0) call post_scalar_data(CS%id_bdott,val ,CS%diag) + if (CS%id_bdot > 0) call post_scalar_data(CS%id_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_bdott_melt > 0 .or. CS%id_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + if (CS%id_bdott_melt > 0) call post_scalar_data(CS%id_bdott_melt,val ,CS%diag) + if (CS%id_bdot_melt > 0) call post_scalar_data(CS%id_bdot_melt ,val*Itime_step,CS%diag) + endif + if (CS%id_bdott_accum > 0 .or. CS%id_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val) + if (CS%id_bdott_accum > 0) call post_scalar_data(CS%id_bdott_accum,val ,CS%diag) + if (CS%id_bdot_accum > 0) call post_scalar_data(CS%id_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_t_area > 0) then !ice sheet area tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) call post_scalar_data(CS%id_t_area,val,CS%diag) endif if (CS%id_g_area > 0 .or. CS%id_f_area > 0) then tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) - if (CS%id_g_area > 0) then + if (CS%id_g_area > 0) then !grounded only ice sheet area call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val) call post_scalar_data(CS%id_g_area,val,CS%diag) endif - if (CS%id_f_area > 0) then + if (CS%id_f_area > 0) then !floating only ice sheet area (ice shelf area) call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val) call post_scalar_data(CS%id_f_area,val,CS%diag) endif endif - call disable_averaging(CS%diag) - call IS_dynamics_post_data(full_time_step, Time, CS%dCS, G) -end subroutine solo_step_ice_shelf + !---ANTARCTICA ONLY---! + if (CS%id_Ant_vaf > 0 .or. CS%id_Ant_dvafdt > 0) & !calculate current volume above floatation (vaf) + call volume_above_floatation(CS%dCS, G, ISS, vaf, hemisphere=0) + if (CS%id_Ant_vaf > 0) call post_scalar_data(CS%id_Ant_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_Ant_dvafdt > 0) call post_scalar_data(CS%id_Ant_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt + if (CS%id_Ant_adott > 0 .or. CS%id_Ant_adot > 0) then !surface accumulation - surface melt + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_adott > 0) call post_scalar_data(CS%id_Ant_adott,val ,CS%diag) + if (CS%id_Ant_adot > 0) call post_scalar_data(CS%id_Ant_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_g_adott > 0 .or. CS%id_Ant_g_adot > 0) then !grounded only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_g_adott > 0) call post_scalar_data(CS%id_Ant_g_adott,val ,CS%diag) + if (CS%id_Ant_g_adot > 0) call post_scalar_data(CS%id_Ant_g_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_f_adott > 0 .or. CS%id_Ant_f_adot > 0) then !floating only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_f_adott > 0) call post_scalar_data(CS%id_Ant_f_adott,val ,CS%diag) + if (CS%id_Ant_f_adot > 0) call post_scalar_data(CS%id_Ant_f_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott > 0 .or. CS%id_Ant_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott > 0) call post_scalar_data(CS%id_Ant_bdott,val ,CS%diag) + if (CS%id_Ant_bdot > 0) call post_scalar_data(CS%id_Ant_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott_melt > 0 .or. CS%id_Ant_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott_melt > 0) call post_scalar_data(CS%id_Ant_bdott_melt,val ,CS%diag) + if (CS%id_Ant_bdot_melt > 0) call post_scalar_data(CS%id_Ant_bdot_melt ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_bdott_accum > 0 .or. CS%id_Ant_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=0) + if (CS%id_Ant_bdott_accum > 0) call post_scalar_data(CS%id_Ant_bdott_accum,val ,CS%diag) + if (CS%id_Ant_bdot_accum > 0) call post_scalar_data(CS%id_Ant_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_Ant_t_area > 0) then !ice sheet area + tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_t_area,val,CS%diag) + endif + if (CS%id_Ant_g_area > 0 .or. CS%id_Ant_f_area > 0) then + tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + if (CS%id_Ant_g_area > 0) then !grounded only ice sheet area + call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_g_area,val,CS%diag) + endif + if (CS%id_Ant_f_area > 0) then !floating only ice sheet area (ice shelf area) + call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=0) + call post_scalar_data(CS%id_Ant_f_area,val,CS%diag) + endif + endif + + !---GREENLAND ONLY---! + if (CS%id_Gr_vaf > 0 .or. CS%id_Gr_dvafdt > 0) & !calculate current volume above floatation (vaf) + call volume_above_floatation(CS%dCS, G, ISS, vaf, hemisphere=1) + if (CS%id_Gr_vaf > 0) call post_scalar_data(CS%id_Gr_vaf ,vaf ,CS%diag) !current vaf + if (CS%id_Gr_dvafdt > 0) call post_scalar_data(CS%id_Gr_dvafdt,(vaf-vaf0_A)*Itime_step,CS%diag) !d(vaf)/dt + if (CS%id_Gr_adott > 0 .or. CS%id_Gr_adot > 0) then !surface accumulation - surface melt + call integrate_over_ice_sheet_area(G, ISS, dh_adott, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_adott > 0) call post_scalar_data(CS%id_Gr_adott,val ,CS%diag) + if (CS%id_Gr_adot > 0) call post_scalar_data(CS%id_Gr_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_g_adott > 0 .or. CS%id_Gr_g_adot > 0) then !grounded only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_g_adott > 0) call post_scalar_data(CS%id_Gr_g_adott,val ,CS%diag) + if (CS%id_Gr_g_adot > 0) call post_scalar_data(CS%id_Gr_g_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_f_adott > 0 .or. CS%id_Gr_f_adot > 0) then !floating only: surface accumulation - surface melt + call masked_var_grounded(G,CS%dCS,dh_adott,tmp) + tmp(:,:) = dh_adott(:,:) - tmp(:,:) + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_f_adott > 0) call post_scalar_data(CS%id_Gr_f_adott,val ,CS%diag) + if (CS%id_Gr_f_adot > 0) call post_scalar_data(CS%id_Gr_f_adot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott > 0 .or. CS%id_Gr_bdot > 0) then !bottom accumulation - bottom melt + call integrate_over_ice_sheet_area(G, ISS, dh_bdott, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott > 0) call post_scalar_data(CS%id_Gr_bdott,val ,CS%diag) + if (CS%id_Gr_bdot > 0) call post_scalar_data(CS%id_Gr_bdot ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott_melt > 0 .or. CS%id_Gr_bdot_melt > 0) then !bottom melt + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott_melt > 0) call post_scalar_data(CS%id_Gr_bdott_melt,val ,CS%diag) + if (CS%id_Gr_bdot_melt > 0) call post_scalar_data(CS%id_Gr_bdot_melt ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_bdott_accum > 0 .or. CS%id_Gr_bdot_accum > 0) then !bottom accumulation + tmp(:,:)=0.0 + do j=js,je ; do i=is,ie + if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j) + enddo; enddo + call integrate_over_ice_sheet_area(G, ISS, tmp, US%Z_to_m, val, hemisphere=1) + if (CS%id_Gr_bdott_accum > 0) call post_scalar_data(CS%id_Gr_bdott_accum,val ,CS%diag) + if (CS%id_Gr_bdot_accum > 0) call post_scalar_data(CS%id_Gr_bdot_accum ,val*Itime_step,CS%diag) + endif + if (CS%id_Gr_t_area > 0) then !ice sheet area + tmp(:,:) = 1.0; call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_t_area,val,CS%diag) + endif + if (CS%id_Gr_g_area > 0 .or. CS%id_Gr_f_area > 0) then + tmp(:,:) = 1.0; call masked_var_grounded(G,CS%dCS,tmp,tmp) + if (CS%id_Gr_g_area > 0) then !grounded only ice sheet area + call integrate_over_ice_sheet_area(G, ISS, tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_g_area,val,CS%diag) + endif + if (CS%id_Gr_f_area > 0) then !floating only ice sheet area (ice shelf area) + call integrate_over_ice_sheet_area(G, ISS, 1.0-tmp, 1.0, val, hemisphere=1) + call post_scalar_data(CS%id_Gr_f_area,val,CS%diag) + endif + endif +end subroutine process_and_post_scalar_data !> \namespace mom_ice_shelf !! diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index 7816df32de..3ed262e5f3 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -50,10 +50,14 @@ module MOM_ice_shelf_dynamics !! on q-points (B grid) [L T-1 ~> m s-1] real, pointer, dimension(:,:) :: v_shelf => NULL() !< the meridional velocity of the ice shelf/sheet !! on q-points (B grid) [L T-1 ~> m s-1] - real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the driving stress of the ice shelf/sheet + real, pointer, dimension(:,:) :: taudx_shelf => NULL() !< the zonal driving stress of the ice shelf/sheet !! on q-points (C grid) [R L2 T-2 ~> Pa] - real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional stress of the ice shelf/sheet + real, pointer, dimension(:,:) :: taudy_shelf => NULL() !< the meridional driving stress of the ice shelf/sheet !! on q-points (C grid) [R L2 T-2 ~> Pa] + real, pointer, dimension(:,:) :: sx_shelf => NULL() !< the zonal surface slope of the ice shelf/sheet + !! on q-points (B grid) [nondim] + real, pointer, dimension(:,:) :: sy_shelf => NULL() !< the meridional surface slope of the ice shelf/sheet + !! on q-points (B grid) [nondim] real, pointer, dimension(:,:) :: u_face_mask => NULL() !< mask for velocity boundary conditions on the C-grid !! u-face - this is because the FEM cares about FACES THAT GET INTEGRATED OVER, !! not vertices. Will represent boundary conditions on computational boundary @@ -107,7 +111,7 @@ module MOM_ice_shelf_dynamics !! of "linearized" basal stress (Pa) [R L3 T-1 ~> kg s-1] !! The exact form depends on basal law exponent and/or whether flow is "hybridized" a la Goldberg 2011 real, pointer, dimension(:,:) :: C_basal_friction => NULL()!< Coefficient in sliding law tau_b = C u^(n_basal_fric), - !! units= Pa (m s-1)^(n_basal_fric) + !! units= Pa (s m-1)^(n_basal_fric) real, pointer, dimension(:,:) :: OD_rt => NULL() !< A running total for calculating OD_av [Z ~> m]. real, pointer, dimension(:,:) :: ground_frac_rt => NULL() !< A running total for calculating ground_frac. real, pointer, dimension(:,:) :: OD_av => NULL() !< The time average open ocean depth [Z ~> m]. @@ -164,6 +168,11 @@ module MOM_ice_shelf_dynamics real :: CFL_factor !< A factor used to limit subcycled advective timestep in uncoupled runs !! i.e. dt <= CFL_factor * min(dx / u) [nondim] + real :: min_h_shelf !< The minimum ice thickness used during ice dynamics [L ~> m]. + real :: min_basal_traction !< The minimum basal traction for grounded ice (Pa m-1 s) [R L T-1 ~> kg m-2 s-1] + real :: max_surface_slope !< The maximum allowed ice-sheet surface slope (to ignore, set to zero) [nondim] + real :: min_ice_visc !< The minimum allowed Glen's law ice viscosity (Pa s), in [R L2 T-1 ~> kg m-1 s-1]. + real :: n_glen !< Nonlinearity exponent in Glen's Law [nondim] real :: eps_glen_min !< Min. strain rate to avoid infinite Glen's law viscosity, [T-1 ~> s-1]. real :: n_basal_fric !< Exponent in sliding law tau_b = C u^(m_slide) [nondim] @@ -221,7 +230,8 @@ module MOM_ice_shelf_dynamics integer :: id_u_shelf = -1, id_v_shelf = -1, id_t_shelf = -1, & id_taudx_shelf = -1, id_taudy_shelf = -1, id_bed_elev = -1, & id_ground_frac = -1, id_col_thick = -1, id_OD_av = -1, id_float_cond = -1, & - id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1 + id_u_mask = -1, id_v_mask = -1, id_ufb_mask =-1, id_vfb_mask = -1, id_t_mask = -1, & + id_sx_shelf = -1, id_sy_shelf = -1 !>@} ! ids for outputting intermediate thickness in advection subroutine (debugging) !>@{ Diagnostic handles for debugging @@ -343,11 +353,13 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) allocate(CS%ice_visc(isd:ied,jsd:jed,CS%visc_qps), source=0.0) allocate(CS%AGlen_visc(isd:ied,jsd:jed), source=2.261e-25) ! [Pa-3 s-1] allocate(CS%basal_traction(isd:ied,jsd:jed), source=0.0) ! [R L3 T-1 ~> kg s-1] - allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (m-1 s)^n_sliding] + allocate(CS%C_basal_friction(isd:ied,jsd:jed), source=5.0e10) ! [Pa (s m-1)^n_sliding] allocate(CS%OD_av(isd:ied,jsd:jed), source=0.0) allocate(CS%ground_frac(isd:ied,jsd:jed), source=0.0) allocate(CS%taudx_shelf(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%taudy_shelf(IsdB:IedB,JsdB:JedB), source=0.0) + allocate(CS%sx_shelf(isd:ied,jsd:jed), source=0.0) + allocate(CS%sy_shelf(isd:ied,jsd:jed), source=0.0) allocate(CS%bed_elev(isd:ied,jsd:jed), source=0.0) allocate(CS%u_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0) allocate(CS%v_bdry_val(IsdB:IedB,JsdB:JedB), source=0.0) @@ -378,7 +390,7 @@ subroutine register_ice_shelf_dyn_restarts(G, US, param_file, CS, restart_CS) call register_restart_field(CS%ground_frac, "ground_frac", .true., restart_CS, & "fractional degree of grounding", "nondim") call register_restart_field(CS%C_basal_friction, "C_basal_friction", .true., restart_CS, & - "basal sliding coefficients", "Pa (m s-1)^n_sliding") + "basal sliding coefficients", "Pa (s m-1)^n_sliding") call register_restart_field(CS%AGlen_visc, "AGlen_visc", .true., restart_CS, & "ice-stiffness parameter", "Pa-3 s-1") call register_restart_field(CS%h_bdry_val, "h_bdry_val", .false., restart_CS, & @@ -491,6 +503,19 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ "The gravitational acceleration of the Earth.", & units="m s-2", default=9.80, scale=US%m_s_to_L_T**2*US%Z_to_m) + call get_param(param_file, mdl, "MIN_H_SHELF", CS%min_h_shelf, & + "min. ice thickness used during ice dynamics", & + units="m", default=0.,scale=US%m_to_L) + call get_param(param_file, mdl, "MIN_BASAL_TRACTION", CS%min_basal_traction, & + "min. allowed basal traction. Input is in [Pa m-1 yr], but is converted when read in to [Pa m-1 s]", & + units="Pa m-1 yr", default=0., scale=365.0*86400.0*US%Pa_to_RLZ_T2*US%L_T_to_m_s) + call get_param(param_file, mdl, "MAX_SURFACE_SLOPE", CS%max_surface_slope, & + "max. allowed ice-sheet surface slope. To ignore, set to zero.", & + units="none", default=0., scale=US%m_to_Z/US%m_to_L) + call get_param(param_file, mdl, "MIN_ICE_VISC", CS%min_ice_visc, & + "min. allowed Glen's law ice viscosity", & + units="Pa s", default=0., scale=US%Pa_to_RL2_T2*US%s_to_T) + call get_param(param_file, mdl, "GLEN_EXPONENT", CS%n_glen, & "nonlinearity exponent in Glen's Law", & units="none", default=3.) @@ -784,6 +809,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ 'x-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) CS%id_taudy_shelf = register_diag_field('ice_shelf_model','taudy_shelf',CS%diag%axesB1, Time, & 'y-driving stress of ice', 'kPa', conversion=1.e-3*US%RLZ_T2_to_Pa) + CS%id_sx_shelf = register_diag_field('ice_shelf_model','sx_shelf',CS%diag%axesB1, Time, & + 'x-surface slope of ice', 'none') + CS%id_sy_shelf = register_diag_field('ice_shelf_model','sy_shelf',CS%diag%axesB1, Time, & + 'y-surface slope of ice', 'none') CS%id_u_mask = register_diag_field('ice_shelf_model','u_mask',CS%diag%axesB1, Time, & 'mask for u-nodes', 'none') CS%id_v_mask = register_diag_field('ice_shelf_model','v_mask',CS%diag%axesB1, Time, & @@ -837,7 +866,7 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, US, Time) do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * ISS%h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -937,26 +966,48 @@ subroutine update_ice_shelf(CS, ISS, G, US, time_step, Time, calve_ice_shelf_ber end subroutine update_ice_shelf -subroutine volume_above_floatation(CS, G, ISS, vaf) +subroutine volume_above_floatation(CS, G, ISS, vaf, hemisphere) type(ice_shelf_dyn_CS), intent(in) :: CS !< The ice shelf dynamics control structure type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. type(ice_shelf_state), intent(in) :: ISS !< A structure with elements that describe !! the ice-shelf state real, intent(out) :: vaf !< area integrated volume above floatation [m3] + integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets + integer :: IS_ID ! local copy of hemisphere real, dimension(SZI_(G),SZJ_(G)) :: vaf_cell !< cell-wise volume above floatation [m3] + integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated integer :: is,ie,js,je,i,j real :: rhoi_rhow, rhow_rhoi if (CS%GL_couple) & call MOM_error(FATAL, "MOM_ice_shelf_dyn, volume above floatation calculation assumes GL_couple=.FALSE..") - vaf_cell(:,:)=0.0 rhoi_rhow = CS%density_ice / CS%density_ocean_avg rhow_rhoi = CS%density_ocean_avg / CS%density_ice is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + if (present(hemisphere)) then + IS_ID=hemisphere + else + IS_ID=-1 + endif + + mask(:,:)=0 + if (IS_ID==0) then !Antarctica (S. Hemisphere) only + do j = js,je; do i = is,ie + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)<=0.0) mask(i,j)=1 + enddo; enddo + elseif (IS_ID==1) then !Greenland (N. Hemisphere) only + do j = js,je; do i = is,ie + if (ISS%hmask(i,j)>0 .and. G%geoLatT(i,j)>0.0) mask(i,j)=1 + enddo; enddo + else !All ice sheets + mask(is:ie,js:je)=ISS%hmask(is:ie,js:je) + endif + + vaf_cell(:,:)=0.0 do j = js,je; do i = is,ie - if (ISS%hmask(i,j)>0) then + if (mask(i,j)>0) then if (CS%bed_elev(i,j) <= 0) then !grounded above sea level vaf_cell(i,j)= (ISS%h_shelf(i,j) * G%US%Z_to_m) * (ISS%area_shelf_h(i,j) * G%US%L_to_m**2) @@ -1007,6 +1058,8 @@ subroutine IS_dynamics_post_data(time_step, Time, CS, G) taud_y(:,:) = CS%taudy_shelf(:,:)*G%IareaBu(:,:) call post_data(CS%id_taudy_shelf, taud_y, CS%diag) endif + if (CS%id_sx_shelf > 0) call post_data(CS%id_sx_shelf, CS%sx_shelf, CS%diag) + if (CS%id_sy_shelf > 0) call post_data(CS%id_sy_shelf, CS%sy_shelf, CS%diag) if (CS%id_ground_frac > 0) call post_data(CS%id_ground_frac, CS%ground_frac, CS%diag) if (CS%id_float_cond > 0) call post_data(CS%id_float_cond, CS%float_cond, CS%diag) if (CS%id_OD_av >0) call post_data(CS%id_OD_av, CS%OD_av,CS%diag) @@ -1122,7 +1175,7 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) else call open_ASCII_file(CS%IS_fileenergy_ascii, trim(CS%IS_energyfile), action=WRITEONLY_FILE) if (abs(CS%timeunit - 86400.0) < 1.0) then - write(CS%IS_fileenergy_ascii,'(" Step,",7x,"Day,"8x,"Energy/Mass,",13x,"Total Mass")') + write(CS%IS_fileenergy_ascii,'(" Step,",7x,"Day,",8x,"Energy/Mass,",13x,"Total Mass")') write(CS%IS_fileenergy_ascii,'(12x,"[days]",10x,"[m2 s-2]",17x,"[kg]")') else if ((CS%timeunit >= 0.99) .and. (CS%timeunit < 1.01)) then @@ -1137,7 +1190,7 @@ subroutine write_ice_shelf_energy(CS, G, US, mass, area, day, time_step) write(time_units,'(9x,"[",es8.2," s] ")') CS%timeunit endif - write(CS%IS_fileenergy_ascii,'(" Step,",7x,"Time,"7x,"Energy/Mass,",13x,"Total Mass")') + write(CS%IS_fileenergy_ascii,'(" Step,",7x,"Time,",7x,"Energy/Mass,",13x,"Total Mass")') write(CS%IS_fileenergy_ascii,'(A25,3x,"[m2 s-2]",17x,"[kg]")') time_units endif endif @@ -1328,7 +1381,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (.not. CS%GL_couple) then do j=G%jsc,G%jec ; do i=G%isc,G%iec - if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) > 0) then + if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) > 0) then CS%ground_frac(i,j) = 1.0 CS%OD_av(i,j) =0.0 endif @@ -1346,7 +1399,7 @@ subroutine ice_shelf_solve_outer(CS, ISS, G, US, u_shlf, v_shlf, taudx, taudy, i if (CS%GL_regularize) then - call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node) + call interpolate_H_to_B(G, ISS%h_shelf, ISS%hmask, H_node, CS%min_h_shelf) do j=G%jsc,G%jec ; do i=G%isc,G%iec nodefloat = 0 @@ -2263,7 +2316,7 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) real :: neumann_val ! [R Z L2 T-2 ~> kg s-2] real :: dxh, dyh,Dx,Dy ! Local grid spacing [L ~> m] real :: grav ! The gravitational acceleration [L2 Z-1 T-2 ~> m s-2] - + real :: scale ! Scaling factor used to ensure surface slope magnitude does not exceed CS%max_surface_slope integer :: i, j, iscq, iecq, jscq, jecq, isd, jsd, ied, jed, is, js, iegq, jegq integer :: giec, gjec, gisc, gjsc, cnt, isc, jsc, iec, jec integer :: i_off, j_off @@ -2289,17 +2342,17 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) if (CS%GL_couple) then do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo - S(i,j) = -CS%bed_elev(i,j) + (OD(i,j) + ISS%h_shelf(i,j)) + S(i,j) = -CS%bed_elev(i,j) + (OD(i,j) + max(ISS%h_shelf(i,j),CS%min_h_shelf)) enddo enddo else ! check whether the ice is floating or grounded do j=jsc-G%domain%njhalo,jec+G%domain%njhalo do i=isc-G%domain%nihalo,iec+G%domain%nihalo - if (rhoi_rhow * ISS%h_shelf(i,j) - CS%bed_elev(i,j) <= 0) then - S(i,j) = (1 - rhoi_rhow)*ISS%h_shelf(i,j) + if (rhoi_rhow * max(ISS%h_shelf(i,j),CS%min_h_shelf) - CS%bed_elev(i,j) <= 0) then + S(i,j) = (1 - rhoi_rhow)*max(ISS%h_shelf(i,j),CS%min_h_shelf) else - S(i,j) = ISS%h_shelf(i,j)-CS%bed_elev(i,j) + S(i,j) = max(ISS%h_shelf(i,j),CS%min_h_shelf)-CS%bed_elev(i,j) endif enddo enddo @@ -2393,14 +2446,21 @@ subroutine calc_shelf_driving_stress(CS, ISS, G, US, taudx, taudy, OD) endif endif - sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sx)) - sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (ISS%h_shelf(i,j) * sy)) + if (CS%max_surface_slope>0) then + scale = min(CS%max_surface_slope/sqrt((sx**2)+(sy**2)),1.0) + sx = scale*sx; sy = scale*sy + endif + + sx_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sx)) + sy_e(i,j) = (-.25 * G%areaT(i,j)) * ((rho * grav) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) * sy)) + + CS%sx_shelf(i,j) = sx ; CS%sy_shelf(i,j) = sy !Stress (Neumann) boundary conditions if (CS%ground_frac(i,j) == 1) then - neumann_val = ((.5 * grav) * (rho * ISS%h_shelf(i,j)**2 - rhow * CS%bed_elev(i,j)**2)) + neumann_val = ((.5 * grav) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2 - rhow * CS%bed_elev(i,j)**2)) else - neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * ISS%h_shelf(i,j)**2)) + neumann_val = (.5 * grav) * ((1-rho/rhow) * (rho * max(ISS%h_shelf(i,j),CS%min_h_shelf)**2)) endif if ((CS%u_face_mask_bdry(I-1,j) == 2) .OR. & ((ISS%hmask(i-1,j) == 0 .OR. ISS%hmask(i-1,j) == 2) .AND. (CS%reentrant_x .OR. (i+i_off /= gisc)))) then @@ -2996,10 +3056,14 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) if ((ISS%hmask(i,j) == 1) .OR. (ISS%hmask(i,j) == 3)) then if (trim(CS%ice_viscosity_compute) == "CONSTANT") then - CS%ice_visc(i,j,1) = 1e15 * (US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T) * (G%areaT(i,j) * ISS%h_shelf(i,j)) + CS%ice_visc(i,j,1) = 1e15 * (US%kg_m3_to_R*US%m_to_L*US%m_s_to_L_T) * & + (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) ! constant viscocity for debugging elseif (trim(CS%ice_viscosity_compute) == "OBS") then - if (CS%AGlen_visc(i,j) >0) CS%ice_visc(i,j,1) = CS%AGlen_visc(i,j) * (G%areaT(i,j) * ISS%h_shelf(i,j)) + if (CS%AGlen_visc(i,j) >0) then + CS%ice_visc(i,j,1) = max(CS%AGlen_visc(i,j) * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)),& + CS%min_ice_visc) + endif ! Here CS%Aglen_visc(i,j) is the ice viscosity [Pa s ~> R L2 T-1] computed from obs and read from a file elseif (model_qp1) then !calculate viscosity at 1 cell-centered quadrature point per cell @@ -3027,9 +3091,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I-1,J) * CS%PhiC(6,i,j) + & v_shlf(I,J-1) * CS%PhiC(4,i,j)) - CS%ice_visc(i,j,1) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + CS%ice_visc(i,j,1) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & (US%s_to_T**2 * ((ux**2 + vy**2) + (ux*vy + 0.25*(uy+vx)**2) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & - (US%Pa_to_RL2_T2*US%s_to_T) + (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) elseif (model_qp4) then !calculate viscosity at 4 quadrature points per cell @@ -3057,9 +3121,9 @@ subroutine calc_shelf_visc(CS, ISS, G, US, u_shlf, v_shlf) (v_shlf(I,J-1) * CS%Phi(4,2*(jq-1)+iq,i,j) + & v_shlf(I-1,J) * CS%Phi(6,2*(jq-1)+iq,i,j)) - CS%ice_visc(i,j,2*(jq-1)+iq) = 0.5 * Visc_coef * (G%areaT(i,j) * ISS%h_shelf(i,j)) * & + CS%ice_visc(i,j,2*(jq-1)+iq) = max(0.5 * Visc_coef * (G%areaT(i,j) * max(ISS%h_shelf(i,j),CS%min_h_shelf)) * & (US%s_to_T**2 * ((ux**2 + vy**2) + (ux*vy + 0.25*(uy+vx)**2) + eps_min**2))**((1.-n_g)/(2.*n_g)) * & - (US%Pa_to_RL2_T2*US%s_to_T) + (US%Pa_to_RL2_T2*US%s_to_T),CS%min_ice_visc) enddo; enddo endif endif @@ -3123,7 +3187,7 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) if (CS%CoulombFriction) then !Effective pressure Hf = max((CS%density_ocean_avg/CS%density_ice) * CS%bed_elev(i,j), 0.0) - fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (ISS%h_shelf(i,j) - Hf)),CS%CF_MinN) + fN = max(fN_scale*((CS%density_ice * CS%g_Earth) * (max(ISS%h_shelf(i,j),CS%min_h_shelf) - Hf)),CS%CF_MinN) fB = alpha * (CS%C_basal_friction(i,j) / (CS%CF_Max * fN))**(CS%CF_PostPeak/CS%n_basal_fric) CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * & @@ -3134,6 +3198,8 @@ subroutine calc_shelf_taub(CS, ISS, G, US, u_shlf, v_shlf) CS%basal_traction(i,j) = ((G%areaT(i,j) * CS%C_basal_friction(i,j)) * (unorm**(CS%n_basal_fric-1))) * & (US%Pa_to_RLZ_T2*US%L_T_to_m_s) endif + + CS%basal_traction(i,j)=max(CS%basal_traction(i,j), CS%min_basal_traction * G%areaT(i,j)) endif enddo enddo @@ -3194,7 +3260,7 @@ subroutine update_OD_ffrac_uncoupled(CS, G, h_shelf) do j=jsd,jed do i=isd,ied - OD = CS%bed_elev(i,j) - rhoi_rhow * h_shelf(i,j) + OD = CS%bed_elev(i,j) - rhoi_rhow * max(h_shelf(i,j),CS%min_h_shelf) if (OD >= 0) then ! ice thickness does not take up whole ocean column -> floating CS%OD_av(i,j) = OD @@ -3640,7 +3706,7 @@ end subroutine update_velocity_masks !> Interpolate the ice shelf thickness from tracer point to nodal points, !! subject to a mask. -subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) +subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node, min_h_shelf) type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf. real, dimension(SZDI_(G),SZDJ_(G)), & intent(in) :: h_shelf !< The ice shelf thickness at tracer points [Z ~> m]. @@ -3650,6 +3716,7 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) real, dimension(SZDIB_(G),SZDJB_(G)), & intent(inout) :: H_node !< The ice shelf thickness at nodal (corner) !! points [Z ~> m]. + real, intent(in) :: min_h_shelf !< The minimum ice thickness used during ice dynamics [L ~> m]. integer :: i, j, isc, iec, jsc, jec, num_h, k, l, ic, jc real :: h_arr(2,2) @@ -3666,7 +3733,7 @@ subroutine interpolate_H_to_B(G, h_shelf, hmask, H_node) num_h = 0 do l=1,2; jc=j-1+l; do k=1,2; ic=i-1+k if (hmask(ic,jc) == 1.0 .or. hmask(ic,jc) == 3.0) then - h_arr(k,l)=h_shelf(ic,jc) + h_arr(k,l)=max(h_shelf(ic,jc),min_h_shelf) num_h = num_h + 1 else h_arr(k,l)=0.0 diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 819e77b2c2..f75707e581 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -510,7 +510,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after forcing') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after forcing', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after forcing', CS%En_sum enddo ; enddo endif @@ -537,7 +537,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 1/2 refraction') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after 1/2 refraction', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 1/2 refraction', CS%En_sum enddo ; enddo ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -567,7 +567,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af halo R", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after correct halo rotation') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after correct halo rotation', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after correct halo rotation', CS%En_sum enddo ; enddo endif @@ -598,7 +598,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af prop", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after propagate') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after propagate', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after propagate', CS%En_sum enddo ; enddo ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -640,7 +640,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides af refr2", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after 2/2 refraction') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after 2/2 refraction', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 2/2 refraction', CS%En_sum enddo ; enddo ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -696,9 +696,9 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after leak", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after background drag') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after background drag', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after background drag', CS%En_sum call sum_En(G, GV, US, CS, CS%TKE_leak_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after background drag') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: loss after background drag', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: loss after background drag', CS%En_sum enddo ; enddo ! Check for En<0 - for debugging, delete later do m=1,CS%nMode ; do fr=1,CS%Nfreq ; do a=1,CS%nAngle @@ -867,7 +867,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after wave", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: before Froude drag') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: before Froude drag', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: before Froude drag', CS%En_sum enddo ; enddo ! save loss term for online budget, may want to add a debug flag later do m=1,CS%nMode ; do fr=1,CS%nFreq @@ -941,9 +941,9 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C call hchksum(CS%En(:,:,:,1,1), "EnergyIntTides after froude", G%HI, haloshift=0, scale=HZ2_T2_to_J_m2) do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'prop_int_tide: after Froude drag') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: after Froude drag', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after Froude drag', CS%En_sum call sum_En(G, GV, US, CS, CS%TKE_Froude_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after Froude drag') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'prop_int_tide: loss after Froude drag', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: loss after Froude drag', CS%En_sum enddo ; enddo ! save loss term for online budget, may want to add a debug flag later do m=1,CS%nMode ; do fr=1,CS%nFreq @@ -1024,7 +1024,7 @@ subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_C CS%TKE_quad_loss_glo_dt(fr,m) - CS%TKE_itidal_loss_glo_dt(fr,m) - & CS%TKE_Froude_loss_glo_dt(fr,m) - CS%TKE_residual_loss_glo_dt(fr,m) - & CS%En_end_glo(fr,m) - if (is_root_pe()) write(stdout,'(A,F18.10)'), "error in Energy budget", CS%error_mode(fr,m) + if (is_root_pe()) write(stdout,'(A,F18.10)') "error in Energy budget", CS%error_mode(fr,m) enddo ; enddo endif @@ -1612,23 +1612,23 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2 enddo if (abs(verif_N -1.0) > threshold_verif) then - write(stdout,'(I5,I5,F18.10)'), i, j, verif_N + write(stdout,'(I5,I5,F18.10)') i, j, verif_N call MOM_error(FATAL, "mismatch integral for N profile") endif if (abs(verif_N2 -1.0) > threshold_verif) then - write(stdout,'(I5,I5,F18.10)'), i, j, verif_N2 + write(stdout,'(I5,I5,F18.10)') i, j, verif_N2 call MOM_error(FATAL, "mismatch integral for N2 profile") endif if (abs(verif_bbl -1.0) > threshold_verif) then - write(stdout,'(I5,I5,F18.10)'), i, j, verif_bbl + write(stdout,'(I5,I5,F18.10)') i, j, verif_bbl call MOM_error(FATAL, "mismatch integral for bbl profile") endif if (abs(verif_stl1 -1.0) > threshold_verif) then - write(stdout,'(I5,I5,F18.10)'), i, j, verif_stl1 + write(stdout,'(I5,I5,F18.10)') i, j, verif_stl1 call MOM_error(FATAL, "mismatch integral for stl1 profile") endif if (abs(verif_stl2 -1.0) > threshold_verif) then - write(stdout,'(I5,I5,F18.10)'), i, j, verif_stl2 + write(stdout,'(I5,I5,F18.10)') i, j, verif_stl2 call MOM_error(FATAL, "mismatch integral for stl2 profile") endif @@ -2108,7 +2108,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss) if (CS%debug) then do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: top of routine') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: top of routine', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: top of routine', CS%En_sum enddo ; enddo endif @@ -2180,7 +2180,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss) if (CS%debug) then do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after propagate_x') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after propagate_x', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: after propagate_x', CS%En_sum enddo ; enddo endif @@ -2191,7 +2191,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss) if (CS%debug) then do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after halo update') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after halo update', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: after halo update', CS%En_sum enddo ; enddo endif ! Apply propagation in y-direction (reflection included) @@ -2210,7 +2210,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss) if (CS%debug) then do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: after propagate_y') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: after propagate_y', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: after propagate_y', CS%En_sum enddo ; enddo endif @@ -2219,7 +2219,7 @@ subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, residual_loss) if (CS%debug) then do m=1,CS%nMode ; do fr=1,CS%Nfreq call sum_En(G, GV, US, CS, CS%En(:,:,:,fr,m), 'propagate: bottom of routine') - if (is_root_pe()) write(stdout,'(A,E18.10)'), 'propagate: bottom of routine', CS%En_sum + if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: bottom of routine', CS%En_sum enddo ; enddo endif diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index cf96017006..0d36ebf6d9 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -1953,10 +1953,10 @@ logical function test_answer(verbose, u, u_true, label, tol) if (abs(u - u_true) > tolerance) test_answer = .true. if (test_answer .or. verbose) then if (test_answer) then - print '(3(a,1pe24.16),x,a,x,a)','computed =',u,' correct =',u_true, & + print '(3(a,1pe24.16),1x,a,1x,a)','computed =',u,' correct =',u_true, & ' err=',u-u_true,' < wrong',label else - print '(2(a,1pe24.16),x,a)','computed =',u,' correct =',u_true,label + print '(2(a,1pe24.16),1x,a)','computed =',u,' correct =',u_true,label endif endif diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index f4ef901868..47550fa93d 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -344,10 +344,10 @@ end subroutine find_N2_bottom !> Returns TKE_itidal_input subroutine get_input_TKE(G, TKE_itidal_input, nFreq, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). + integer, intent(in) :: nFreq !< number of frequencies real, dimension(SZI_(G),SZJ_(G),nFreq), & intent(out) :: TKE_itidal_input !< The energy input to the internal waves !! [H Z2 T-3 ~> m3 s-3 or W m-2]. - integer, intent(in) :: nFreq !< number of frequencies type(int_tide_input_CS), target :: CS !< A pointer that is set to point to the control !! structure for the internal tide input module. integer :: i,j,fr @@ -361,9 +361,9 @@ end subroutine get_input_TKE !> Returns barotropic tidal velocities subroutine get_barotropic_tidal_vel(G, vel_btTide, nFreq, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). + integer, intent(in) :: nFreq !< number of frequencies real, dimension(SZI_(G),SZJ_(G),nFreq), & intent(out) :: vel_btTide !< Barotropic velocity read from file [L T-1 ~> m s-1]. - integer, intent(in) :: nFreq !< number of frequencies type(int_tide_input_CS), target :: CS !< A pointer that is set to point to the control !! structure for the internal tide input module. integer :: i,j,fr diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index cd781169af..147c48eebd 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -342,8 +342,8 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (CS%use_real_BL_depth .and. .not. present(Hbl)) then - call MOM_error(FATAL,"Attempting to use real boundary layer depth for ideal age tracers, & - but no valid boundary layer scheme was found") + call MOM_error(FATAL, "Attempting to use real boundary layer depth for ideal age tracers, " & + // "but no valid boundary layer scheme was found") endif if (CS%use_real_BL_depth .and. present(Hbl)) then