← back to index

src/core/MOM_PressureForce_FV.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Finite volume pressure gradient (integrated by quadrature or analytically)
6module MOM_PressureForce_FV
7
8use MOM_debugging, only : hchksum, uvchksum
9use MOM_diag_mediator, only : post_data, register_diag_field
10use MOM_diag_mediator, only : safe_alloc_ptr, diag_ctrl, time_type
11use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
12use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
13use MOM_grid, only : ocean_grid_type
14use MOM_PressureForce_Mont, only : set_pbce_Bouss, set_pbce_nonBouss
15use MOM_self_attr_load, only : calc_SAL, SAL_CS
16use MOM_tidal_forcing, only : calc_tidal_forcing, tidal_forcing_CS
17use MOM_tidal_forcing, only : calc_tidal_forcing_legacy
18use MOM_unit_scaling, only : unit_scale_type
19use MOM_variables, only : thermo_var_ptrs, accel_diag_ptrs
20use MOM_verticalGrid, only : verticalGrid_type
21use MOM_EOS, only : calculate_density, calculate_spec_vol, EOS_domain
22use MOM_density_integrals, only : int_density_dz, int_specific_vol_dp
23use MOM_density_integrals, only : int_density_dz_generic_plm, int_density_dz_generic_ppm
24use MOM_density_integrals, only : int_spec_vol_dp_generic_plm
25use MOM_density_integrals, only : int_density_dz_generic_pcm, int_spec_vol_dp_generic_pcm
26use MOM_density_integrals, only : diagnose_mass_weight_Z, diagnose_mass_weight_p
27use MOM_ALE, only : TS_PLM_edge_values, TS_PPM_edge_values, TS_PLM_WLS_edge_values, ALE_CS
28
29implicit none ; private
30
31#include <MOM_memory.h>
32
33public PressureForce_FV_init
34public PressureForce_FV_Bouss, PressureForce_FV_nonBouss
35
36! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
37! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
38! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
39! vary with the Boussinesq approximation, the Boussinesq variant is given first.
40
41#ifdef __NVCOMPILER_OPENMP_GPU
42integer, parameter :: default_nkblock = 0 !< Default k block size for PLM density integrals [nondim]
43integer, parameter :: default_njblock_plm = 0 !< Default j block size for PLM density integrals [nondim]
44#else
45integer, parameter :: default_nkblock = 1 !< Default k block size for PLM density integrals [nondim]
46integer, parameter :: default_njblock_plm = 1 !< Default j block size for PLM density integrals [nondim]
47#endif
48integer, parameter :: default_niblock_plm = 0 !< Default i block size for PLM density integrals [nondim]
49
50!> Finite volume pressure gradient control structure
51type, public :: PressureForce_FV_CS ; private
52 logical :: initialized = .false. !< True if this control structure has been initialized.
53 logical :: calculate_SAL = .false. !< If true, calculate self-attraction and loading.
54 logical :: sal_use_bpa = .false. !< If true, use bottom pressure anomaly instead of SSH
55 !! to calculate SAL.
56 logical :: tides = .false. !< If true, apply tidal momentum forcing.
57 real :: rho_ref !< The reference density that is subtracted off when calculating pressure
58 !! gradient forces [R ~> kg m-3].
59 logical :: rho_ref_bug !< If true, recover a bug that mixes GV%Rho0 and CS%rho_ref in Boussinesq mode.
60 real :: GFS_scale !< A scaling of the surface pressure gradients to
61 !! allow the use of a reduced gravity model [nondim].
62 type(time_type), pointer :: Time !< A pointer to the ocean model's clock.
63 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
64 !! timing of diagnostic output.
65 integer :: MassWghtInterp !< A flag indicating whether and how to use mass weighting in T/S interpolation
66 logical :: correction_intxpa !< If true, apply a correction to the value of intxpa at a selected
67 !! interface under ice, using matching at the end values along with a
68 !! 5-point quadrature integral of the hydrostatic pressure or height
69 !! changes along that interface. The selected interface is either at the
70 !! ocean's surface or in the interior, depending on reset_intxpa_integral.
71 logical :: reset_intxpa_integral !< If true and the surface displacement between adjacent cells
72 !! exceeds the vertical grid spacing, reset intxpa at the interface below
73 !! a trusted interior cell. (This often applies in ice shelf cavities.)
74 logical :: MassWghtInterpVanOnly !< If true, don't do mass weighting of T/S interpolation unless vanished
75 logical :: reset_intxpa_flattest !< If true, use flattest interface rather than top for reset integral
76 !! in cases where no best nonvanished interface
77 real :: h_nonvanished !< A minimal layer thickness that indicates that a layer is thick enough
78 !! to usefully reestimate the pressure integral across the interface
79 !! below it [H ~> m or kg m-2]
80 logical :: use_inaccurate_pgf_rho_anom !< If true, uses the older and less accurate
81 !! method to calculate density anomalies, as used prior to
82 !! March 2018.
83 logical :: boundary_extrap !< Indicate whether high-order boundary
84 !! extrapolation should be used within boundary cells
85
86 logical :: reconstruct !< If true, polynomial profiles of T & S will be
87 !! reconstructed and used in the integrals for the
88 !! finite volume pressure gradient calculation.
89 !! The default depends on whether regridding is being used.
90
91 integer :: Recon_Scheme !< Order of the polynomial of the reconstruction of T & S
92 !! for the finite volume pressure gradient calculation.
93 !! By the default (1) is for a piecewise linear method
94
95 logical :: debug !< If true, write verbose checksums for debugging purposes.
96 logical :: use_SSH_in_Z0p !< If true, adjust the height at which the pressure used in the
97 !! equation of state is 0 to account for the displacement of the sea
98 !! surface including adjustments for atmospheric or sea-ice pressure.
99 logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
100 logical :: bq_sal_tides = .false. !< If true, use an alternative method for SAL and tides
101 !! in Boussinesq mode
102 integer :: tides_answer_date = 99991231 !< Recover old answers with tides
103 integer :: id_e_tide = -1 !< Diagnostic identifier
104 integer :: id_e_tidal_eq = -1 !< Diagnostic identifier
105 integer :: id_e_tidal_sal = -1 !< Diagnostic identifier
106 integer :: id_e_sal = -1 !< Diagnostic identifier
107 integer :: id_rho_pgf = -1 !< Diagnostic identifier
108 integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier
109 integer :: id_p_stanley = -1 !< Diagnostic identifier
110 integer :: id_MassWt_u = -1 !< Diagnostic identifier
111 integer :: id_MassWt_v = -1 !< Diagnostic identifier
112 integer :: id_sal_u = -1 !< Diagnostic identifier
113 integer :: id_sal_v = -1 !< Diagnostic identifier
114 integer :: id_tides_u = -1 !< Diagnostic identifier
115 integer :: id_tides_v = -1 !< Diagnostic identifier
116 integer :: nkblock !< Vertical block size for PLM density integrals [nondim]
117 integer :: niblock_plm !< i-tile size for PLM density integrals, 0 = full domain [nondim]
118 integer :: njblock_plm !< j-tile size for PLM density integrals, 0 = full domain [nondim]
119 type(SAL_CS), pointer :: SAL_CSp => NULL() !< SAL control structure
120 type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure
121end type PressureForce_FV_CS
122
123contains
124
125!> \brief Non-Boussinesq analytically-integrated finite volume form of pressure gradient
126!!
127!! Determines the acceleration due to hydrostatic pressure forces, using
128!! the analytic finite volume form of the Pressure gradient, and does not
129!! make the Boussinesq approximation.
130!!
131!! To work, the following fields must be set outside of the usual (is:ie,js:je)
132!! range before this subroutine is called:
133!! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
1340subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp, p_atm, pbce, eta)
135 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
136 type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
137 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
138 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> kg m-2]
139 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
140 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: PFu !< Zonal acceleration [L T-2 ~> m s-2]
141 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2]
142 type(PressureForce_FV_CS), intent(in) :: CS !< Finite volume PGF control structure
143 type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
144 type(accel_diag_ptrs), pointer :: ADp !< Acceleration diagnostic pointers
145 real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean
146 !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
147 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure
148 !! anomaly in each layer due to eta anomalies
149 !! [L2 T-2 H-1 ~> m4 s-2 kg-1].
150 real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The total column mass used to
151 !! calculate PFu and PFv [H ~> kg m-2].
152 ! Local variables
1530 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: p ! Interface pressure [R L2 T-2 ~> Pa].
154 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
1550 T_tmp, & ! Temporary array of temperatures where layers that are lighter
156 ! than the mixed layer have the mixed layer's properties [C ~> degC].
1570 S_tmp ! Temporary array of salinities where layers that are lighter
158 ! than the mixed layer have the mixed layer's properties [S ~> ppt].
159 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1600 S_t, S_b, & ! Top and bottom edge values for linear reconstructions
161 ! of salinity within each layer [S ~> ppt].
1620 T_t, T_b ! Top and bottom edge values for linear reconstructions
163 ! of temperature within each layer [C ~> degC].
164 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1650 dza, & ! The change in geopotential anomaly between the top and bottom
166 ! of a layer [L2 T-2 ~> m2 s-2].
1670 intp_dza ! The vertical integral in depth of the pressure anomaly less
168 ! the pressure anomaly at the top of the layer [R L4 T-4 ~> Pa m2 s-2].
169 real, dimension(SZI_(G),SZJ_(G)) :: &
1700 dp, & ! The (positive) change in pressure across a layer [R L2 T-2 ~> Pa].
1710 SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if
172 ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].
1730 pbot, & ! Total bottom pressure for self-attraction and loading. Used if
174 ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].
1750 e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
1760 e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources [Z ~> m].
1770 e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
178 ! specific to tides [Z ~> m].
1790 e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing, used for recovering
180 ! old answers only [Z ~> m].
1810 dM ! The barotropic adjustment to the Montgomery potential to
182 ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
183 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1840 za ! The geopotential anomaly (i.e. g*e + alpha_0*pressure) at the
185 ! interfaces [L2 T-2 ~> m2 s-2].
186
1870 real, dimension(SZI_(G)) :: Rho_cv_BL ! The coordinate potential density in the deepest variable
188 ! density near-surface layer [R ~> kg m-3].
189 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
1900 intx_za ! The zonal integral of the geopotential anomaly along the
191 ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2].
192 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
1930 intx_dza ! The change in intx_za through a layer [L2 T-2 ~> m2 s-2].
194 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
1950 inty_za ! The meridional integral of the geopotential anomaly along the
196 ! interfaces, divided by the grid spacing [L2 T-2 ~> m2 s-2].
197 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
1980 inty_dza ! The change in inty_za through a layer [L2 T-2 ~> m2 s-2].
199 real, dimension(SZI_(G),SZJ_(G)) :: &
2000 T_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC]
2010 S_top ! Salinity of top layer used with correction_intxpa [S ~> ppt]
202 real, dimension(SZIB_(G),SZJ_(G)) :: &
2030 intx_za_cor ! Correction for curvature in intx_za [L2 T-2 ~> m2 s-2]
204 real, dimension(SZI_(G),SZJB_(G)) :: &
2050 inty_za_cor ! Correction for curvature in inty_za [L2 T-2 ~> m2 s-2]
206
207 ! These variables are used with reset_intxpa_integral. The values are taken from different
208 ! interfaces as a function of position.
209 real, dimension(SZIB_(G),SZJ_(G)) :: &
2100 T_int_W, T_int_E, & ! Temperatures on the reference interface to the east and west of a u-point [C ~> degC]
2110 S_int_W, S_int_E, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt]
2120 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]
2130 intx_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface
214 ! from the value that would be obtained from assuming that pressure varies
215 ! linearly with depth along that interface [R L2 T-2 ~> Pa].
2160 dp_int_x, & ! The change in x in pressure along the reference interface [R L2 T-2 ~> Pa]
2170 intx_za_cor_ri ! The correction to intx_za based on the reference interface calculations [L2 T-2 ~> m2 s-2]
218 real, dimension(SZI_(G),SZJB_(G)) :: &
2190 T_int_S, T_int_N, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC]
2200 S_int_S, S_int_N, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt]
2210 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]
2220 inty_za_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface
223 ! from the value that would be obtained from assuming that pressure varies
224 ! linearly with depth along that interface [L2 T-2 ~> m2 s-2].
2250 dp_int_y, & ! The change in y in geopotenial height along the reference interface [R L2 T-2 ~> Pa]
2260 inty_za_cor_ri ! The correction to inty_za based on the reference interface calculations [L2 T-2 ~> m2 s-2]
227 logical, dimension(SZIB_(G),SZJ_(G)) :: &
2280 seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate
229 ! of the curvature terms in the intx_pa.
230 logical, dimension(SZI_(G),SZJB_(G)) :: &
2310 seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate
232 ! of the curvature terms in the inty_pa.
233 real, dimension(SZIB_(G),SZJ_(G)) :: &
2340 delta_p_x ! If using flattest interface for reset integral, store x interface
235 ! differences [R L2 T-2 ~> Pa]
236 real, dimension(SZI_(G),SZJB_(G)) :: &
2370 delta_p_y ! If using flattest interface for reset integral, store y interface
238 ! differences [R L2 T-2 ~> Pa]
239 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
2400 MassWt_u ! The fractional mass weighting at a u-point [nondim].
241 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
2420 MassWt_v ! The fractional mass weighting at a v-point [nondim].
2430 real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
244 ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
245 real :: dp_sfc ! The change in surface pressure between adjacent cells [R L2 T-2 ~> Pa]
246
247 real :: dp_neglect ! A thickness that is so small it is usually lost
248 ! in roundoff and can be neglected [R L2 T-2 ~> Pa].
249 real :: p_nonvanished ! nonvanshed pressure [R L2 T-2 ~> Pa]
250 real :: I_gEarth ! The inverse of GV%g_Earth [T2 Z L-2 ~> s2 m-1]
251 real :: alpha_anom ! The in-situ specific volume, averaged over a
252 ! layer, less alpha_ref [R-1 ~> m3 kg-1].
253 logical :: use_p_atm ! If true, use the atmospheric pressure.
254 logical :: use_ALE ! If true, use an ALE pressure reconstruction.
255 logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
256 logical :: do_more_k ! If true, there are still points where a flatter interface remains to be found.
2570 type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
258
259 real :: alpha_ref ! A reference specific volume [R-1 ~> m3 kg-1] that is used
260 ! to reduce the impact of truncation errors.
2610 real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
262 real :: Pa_to_H ! A factor to convert from Pa to the thickness units (H)
263 ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1].
264 real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure
265 ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
266 real :: T5(5) ! Temperatures and salinities at five quadrature points [C ~> degC]
267 real :: S5(5) ! Salinities at five quadrature points [S ~> ppt]
268 real :: p5(5) ! Pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa]
269 real :: SpV5(5) ! Specific volume anomalies at five quadrature points [R-1 ~> m3 kg-1]
270 real :: wt_R ! A weighting factor [nondim]
271
272 ! real :: oneatm ! 1 standard atmosphere of pressure in [R L2 T-2 ~> Pa]
273 real, parameter :: C1_6 = 1.0/6.0 ! [nondim]
274 real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
275 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
276 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
277 integer, dimension(2) :: EOSdom_u ! The i-computational domain for the equation of state at u-velocity points
278 integer, dimension(2) :: EOSdom_v ! The i-computational domain for the equation of state at v-velocity points
279 integer :: i, j, k, m
280
2810 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2820 nkmb=GV%nk_rho_varies
2830 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
2840 EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)
2850 EOSdom_u(1) = Isq - (G%IsdB-1) ; EOSdom_u(2) = Ieq - (G%IsdB-1)
2860 EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1)
287
2880 if (.not.CS%initialized) call MOM_error(FATAL, &
2890 "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
290
2910 if (CS%use_stanley_pgf) call MOM_error(FATAL, &
292 "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet "//&
2930 "implemented in non-Boussinesq mode.")
294
2950 use_p_atm = associated(p_atm)
2960 use_EOS = associated(tv%eqn_of_state)
2970 use_ALE = .false.
2980 if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS
299
3000 H_to_RL2_T2 = GV%g_Earth*GV%H_to_RZ
3010 dp_neglect = GV%g_Earth*GV%H_to_RZ * GV%H_subroundoff
3020 alpha_ref = 1.0 / CS%rho_ref
3030 I_gEarth = 1.0 / GV%g_Earth
3040 p_nonvanished = GV%g_Earth*GV%H_to_RZ*CS%h_nonvanished
305
3060 if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then
3070 MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0
308 endif
309
3100 if (use_p_atm) then
311 !$OMP parallel do default(shared)
3120 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
3130 p(i,j,1) = p_atm(i,j)
314 enddo ; enddo
315 else
316 ! oneatm = 101325.0 * US%Pa_to_RL2_T2 ! 1 atm scaled to [R L2 T-2 ~> Pa]
317 !$OMP parallel do default(shared)
3180 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
3190 p(i,j,1) = 0.0 ! or oneatm
320 enddo ; enddo
321 endif
322 !$OMP parallel do default(shared)
3230 do j=Jsq,Jeq+1 ; do k=2,nz+1 ; do i=Isq,Ieq+1
3240 p(i,j,K) = p(i,j,K-1) + H_to_RL2_T2 * h(i,j,k-1)
325 enddo ; enddo ; enddo
326
3270 if (use_EOS) then
328 ! With a bulk mixed layer, replace the T & S of any layers that are
329 ! lighter than the buffer layer with the properties of the buffer
330 ! layer. These layers will be massless anyway, and it avoids any
331 ! formal calculations with hydrostatically unstable profiles.
3320 if (nkmb>0) then
3330 tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp
3340 tv_tmp%eqn_of_state => tv%eqn_of_state
3350 do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
336 !$OMP parallel do default(shared) private(Rho_cv_BL)
3370 do j=Jsq,Jeq+1
3380 do k=1,nkmb ; do i=Isq,Ieq+1
3390 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
340 enddo ; enddo
341 call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
3420 tv%eqn_of_state, EOSdom)
3430 do k=nkmb+1,nz ; do i=Isq,Ieq+1
3440 if (GV%Rlay(k) < Rho_cv_BL(i)) then
3450 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
346 else
3470 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
348 endif
349 enddo ; enddo
350 enddo
351 else
3520 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
3530 tv_tmp%eqn_of_state => tv%eqn_of_state
354 endif
355 endif
356
357 ! If regridding is activated, do a linear reconstruction of salinity
358 ! and temperature across each layer. The subscripts 't' and 'b' refer
359 ! to top and bottom values within each layer (these are the only degrees
360 ! of freedom needed to know the linear profile).
3610 if ( use_ALE .and. (CS%Recon_Scheme == 1) ) then
3620 call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap)
3630 elseif ( use_ALE .and. (CS%Recon_Scheme == 2) ) then
3640 call TS_PPM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap)
3650 elseif ( use_ALE .and. (CS%Recon_Scheme == 3) ) then
3660 call TS_PLM_WLS_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
3670 elseif (CS%reset_intxpa_integral) then
3680 do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
3690 T_b(i,j,k) = tv%T(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k)
370 enddo ; enddo ; enddo
371 endif
372
373 !$OMP parallel do default(shared) private(alpha_anom,dp)
3740 do k=1,nz
375 ! Calculate 4 integrals through the layer that are required in the
376 ! subsequent calculation.
3770 if (use_EOS) then
3780 if ( use_ALE .and. CS%Recon_Scheme > 0 ) then
3790 if ( CS%Recon_Scheme == 1 .or. CS%Recon_Scheme == 3 ) then
380 call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), &
381 p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, &
382 tv%eqn_of_state, US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), &
383 P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp, &
3840 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
3850 elseif ( CS%Recon_Scheme == 2 ) then
386 call MOM_error(FATAL, "PressureForce_FV_nonBouss: "//&
3870 "int_spec_vol_dp_generic_ppm does not exist yet.")
388 ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), &
389 ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), &
390 ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), &
391 ! intx_dza(:,:,k), inty_dza(:,:,k), P_surf=p(:,:,1), MassWghtInterp=CS%MassWghtInterp)
392 endif
393 else
394 call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), &
395 p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, &
396 US, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), &
397 inty_dza(:,:,k), bathyP=p(:,:,nz+1), P_surf=p(:,:,1), dP_tiny=dp_neglect, &
398 MassWghtInterp=CS%MassWghtInterp, &
3990 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
400 endif
4010 if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) &
402 call diagnose_mass_weight_p(p(:,:,K), p(:,:,K+1), p(:,:,nz+1), p(:,:,1), dp_neglect, CS%MassWghtInterp, &
403 G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k), &
4040 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, p_nv=p_nonvanished)
405 else
4060 alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref
4070 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4080 dp(i,j) = H_to_RL2_T2 * h(i,j,k)
4090 dza(i,j,k) = alpha_anom * dp(i,j)
4100 intp_dza(i,j,k) = 0.5 * alpha_anom * dp(i,j)**2
411 enddo ; enddo
4120 do j=js,je ; do I=Isq,Ieq
4130 intx_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i+1,j))
414 enddo ; enddo
4150 do J=Jsq,Jeq ; do i=is,ie
4160 inty_dza(i,j,k) = 0.5 * alpha_anom * (dp(i,j)+dp(i,j+1))
417 enddo ; enddo
418 endif
419 enddo
420
421 ! The bottom geopotential anomaly is calculated first so that the increments
422 ! to the geopotential anomalies can be reused. Alternately, the surface
423 ! geopotential could be calculated directly with separate calls to
424 ! int_specific_vol_dp with alpha_ref=0, and the anomalies used going
425 ! downward, which would relieve the need for dza, intp_dza, intx_dza, and
426 ! inty_dza to be 3-D arrays.
427
428 ! Sum vertically to determine the surface geopotential anomaly.
429 !$OMP parallel do default(shared)
4300 do j=Jsq,Jeq+1
4310 do i=Isq,Ieq+1
4320 za(i,j,nz+1) = alpha_ref*p(i,j,nz+1) - GV%g_Earth*G%bathyT(i,j)
433 enddo
4340 do k=nz,1,-1 ; do i=Isq,Ieq+1
4350 za(i,j,K) = za(i,j,K+1) + dza(i,j,k)
436 enddo ; enddo
437 enddo
438
439 ! Calculate and add self-attraction and loading (SAL) geopotential height anomaly to interface height.
4400 if (CS%calculate_SAL) then
4410 if (CS%sal_use_bpa) then
442 !$OMP parallel do default(shared)
4430 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4440 pbot(i,j) = p(i,j,nz+1)
445 enddo ; enddo
4460 call calc_SAL(pbot, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
447 else
448 !$OMP parallel do default(shared)
4490 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4500 SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref
451 ! Remove above sea level topography at floodable cells
4520 SSH(i,j) = SSH(i,j) - max(-G%bathyT(i,j)-G%meanSL(i,j), 0.0)
453 enddo ; enddo
4540 call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
455 endif
456
457 ! This gives new answers after the change of separating SAL from tidal forcing module.
4580 if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then
459 !$OMP parallel do default(shared)
4600 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4610 za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal(i,j)
462 enddo ; enddo
463 endif
464 endif
465
466 ! Calculate and add tidal geopotential height anomaly to interface height.
4670 if (CS%tides) then
4680 if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq)) then
4690 call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
470 !$OMP parallel do default(shared)
4710 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4720 za(i,j,1) = za(i,j,1) - GV%g_Earth * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
473 enddo ; enddo
474 else ! This block recreates older answers with tides.
4750 if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0
476 call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &
4770 G, US, CS%tides_CSp)
478 !$OMP parallel do default(shared)
4790 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4800 za(i,j,1) = za(i,j,1) - GV%g_Earth * e_sal_and_tide(i,j)
481 enddo ; enddo
482 endif
483 endif
484
485 ! Find the height anomalies at the interfaces. If there are no tides and no SAL,
486 ! there is no need to correct za, but omitting this changes answers at roundoff.
4870 do k=1,nz
488 !$OMP parallel do default(shared)
4890 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
4900 za(i,j,K+1) = za(i,j,K) - dza(i,j,k)
491 enddo ; enddo
492 enddo
493
4940 if (CS%debug) then
4950 call hchksum(za, "Pre-correction za", G%HI, haloshift=1, unscale=US%L_T_to_m_s**2)
4960 call hchksum(p, "Pre-correction p", G%HI, haloshift=1, unscale=US%RL2_T2_to_Pa)
497 endif
498
499 ! With an ice-shelf or icebergs, this linearity condition might need to be applied
500 ! to a sub-surface interface.
5010 if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then
502 ! Determine surface temperature and salinity for use in the pressure gradient corrections
5030 if (use_ALE .and. (CS%Recon_Scheme > 0)) then
5040 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
5050 T_top(i,j) = T_t(i,j,1) ; S_top(i,j) = S_t(i,j,1)
506 enddo ; enddo
507 else
5080 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
5090 T_top(i,j) = tv%T(i,j,1) ; S_top(i,j) = tv%S(i,j,1)
510 enddo ; enddo
511 endif
512 endif
513
5140 if (CS%correction_intxpa) then
515 ! This version makes a 5 point quadrature correction for hydrostatic variations in surface
516 ! pressure under ice.
517 !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5)
5180 do j=js,je ; do I=Isq,Ieq
5190 intx_za_cor(I,j) = 0.0
5200 dp_sfc = (p(i+1,j,1) - p(i,j,1))
521 ! If the changes in pressure and height anomaly were explicable by just a hydrostatic balance,
522 ! the implied specific volume would be SpV_implied = alpha_ref - (dza_x / dp_x)
5230 if (dp_sfc * (alpha_ref*dp_sfc - (za(i+1,j,1)-za(i,j,1))) > 0.0) then
5240 T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j)
5250 S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j)
5260 p5(1) = p(i,j,1) ; p5(5) = p(i+1,j,1)
5270 do m=2,4
5280 wt_R = 0.25*real(m-1)
5290 T5(m) = T5(1) + (T5(5)-T5(1))*wt_R
5300 S5(m) = S5(1) + (S5(5)-S5(1))*wt_R
5310 p5(m) = p5(1) + (p5(5)-p5(1))*wt_R
532 enddo !m
5330 call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref)
534 ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression.
5350 intx_za_cor(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc
536 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12
537 endif
5380 intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1)) + intx_za_cor(I,j)
539 enddo ; enddo
540 !$OMP parallel do default(shared) private(dp_sfc,T5,S5,p5,wt_R,SpV5)
5410 do J=Jsq,Jeq ; do i=is,ie
5420 inty_za_cor(i,J) = 0.0
5430 dp_sfc = (p(i,j+1,1) - p(i,j,1))
5440 if (dp_sfc * (alpha_ref*dp_sfc - (za(i,j+1,1)-za(i,j,1))) > 0.0) then
545 ! The pressure/depth relationship has a positive implied specific volume.
5460 T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1)
5470 S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1)
5480 p5(1) = p(i,j,1) ; p5(5) = p(i,j+1,1)
5490 do m=2,4
5500 wt_R = 0.25*real(m-1)
5510 T5(m) = T5(1) + (T5(5)-T5(1))*wt_R
5520 S5(m) = S5(1) + (S5(5)-S5(1))*wt_R
5530 p5(m) = p5(1) + (p5(5)-p5(1))*wt_R
554 enddo !m
5550 call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref)
556 ! See the Boussinesq calculation of inty_pa_cor for the derivation of the following expression.
5570 inty_za_cor(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * dp_sfc
558 endif
5590 inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1)) + inty_za_cor(i,J)
560 enddo ; enddo
561 else
562 ! This order of integrating upward and then downward again is necessary with
563 ! a nonlinear equation of state, so that the surface geopotentials will go
564 ! linearly between the values at thickness points, but the bottom geopotentials
565 ! will not now be linear at the sub-grid-scale. Doing this ensures no motion
566 ! with flat isopycnals, even with a nonlinear equation of state.
567 !$OMP parallel do default(shared)
5680 do j=js,je ; do I=Isq,Ieq
5690 intx_za(I,j,1) = 0.5*(za(i,j,1) + za(i+1,j,1))
570 enddo ; enddo
571 !$OMP parallel do default(shared)
5720 do J=Jsq,Jeq ; do i=is,ie
5730 inty_za(i,J,1) = 0.5*(za(i,j,1) + za(i,j+1,1))
574 enddo ; enddo
575 endif
576
5770 do k=1,nz
578 !$OMP parallel do default(shared)
5790 do j=js,je ; do I=Isq,Ieq
5800 intx_za(I,j,K+1) = intx_za(I,j,K) - intx_dza(I,j,k)
581 enddo ; enddo
582 enddo
5830 do k=1,nz
584 !$OMP parallel do default(shared)
5850 do J=Jsq,Jeq ; do i=is,ie
5860 inty_za(i,J,K+1) = inty_za(i,J,K) - inty_dza(i,J,k)
587 enddo ; enddo
588 enddo
589
5900 if (CS%debug) then
591 call uvchksum("Prelim int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, &
5920 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
593 call uvchksum("Prelim int[xy]_dza", intx_dza, inty_dza, G%HI, haloshift=0, &
5940 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
595 endif
596
5970 if (CS%reset_intxpa_integral) then
598 ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is
599 ! reset intx_za there, then adjust intx_za throughout the water column.
600
601 ! Zero out the 2-d arrays that will be set from various reference interfaces.
6020 T_int_W(:,:) = 0.0 ; S_int_W(:,:) = 0.0 ; p_int_W(:,:) = 0.0
6030 T_int_E(:,:) = 0.0 ; S_int_E(:,:) = 0.0 ; p_int_E(:,:) = 0.0
6040 intx_za_nonlin(:,:) = 0.0 ; intx_za_cor_ri(:,:) = 0.0 ; dp_int_x(:,:) = 0.0
6050 do j=js,je ; do I=Isq,Ieq
6060 seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.)
6070 delta_p_x(I,j) = 0.0
608 enddo ; enddo
609
6100 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
6110 if ((p(i+1,j,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i+1,j,1))) then
612 ! This is the typical case in the open ocean, so use the topmost interface.
6130 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
6140 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
6150 p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
6160 intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
6170 dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
6180 seek_x_cor(I,j) = .false.
619 endif
620 endif ; enddo ; enddo
621
6220 do k=1,nz
6230 do_more_k = .false.
6240 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
625 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not
626 ! activated in the subgrid interpolation.
6270 if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i+1,j,k) > CS%h_nonvanished)) .and. &
628 (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
629 ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at
630 ! the interface below this cell (it might have quadratic pressure dependence if sloped)
6310 T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
6320 S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
6330 p_int_W(I,j) = p(i,j,K+1) ; p_int_E(I,j) = p(i+1,j,K+1)
634
6350 intx_za_nonlin(I,j) = intx_za(I,j,K+1) - 0.5*(za(i,j,K+1) + za(i+1,j,K+1))
6360 dp_int_x(I,j) = p(i+1,j,K+1)-p(i,j,K+1)
6370 seek_x_cor(I,j) = .false.
638 else
6390 do_more_k = .true.
640 endif
641 endif ; enddo ; enddo
6420 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.
643 enddo
644
6450 if (do_more_k) then
6460 if (CS%reset_intxpa_flattest) then
647 ! There are still points where a correction is needed, so use flattest interface
6480 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
649 ! choose top layer first
6500 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
6510 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
6520 p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
6530 intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
6540 dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
6550 delta_p_x(I,j) = abs(p(i+1,j,1)-p(i,j,1))
6560 do k=1,nz
6570 if (abs(p(i+1,j,k+1)-p(i,j,k+1)) < delta_p_x(I,j)) then
658 ! bottom of layer is less sloped than top. Use this layer
6590 delta_p_x(I,j) = abs(p(i+1,j,k+1)-p(i,j,k+1))
6600 T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
6610 S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
6620 p_int_W(I,j) = p(i,j,K+1) ; p_int_E(I,j) = p(i+1,j,K+1)
6630 intx_za_nonlin(I,j) = intx_za(I,j,K+1) - 0.5*(za(i,j,K+1) + za(i+1,j,K+1))
6640 dp_int_x(I,j) = p(i+1,j,K+1)-p(i,j,K+1)
665 endif
666 enddo
6670 seek_x_cor(I,j) = .false.
668 endif ; enddo ; enddo
669 else
670 ! There are still points where a correction is needed, so use the top interface.
6710 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
6720 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
6730 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
6740 p_int_W(I,j) = p(i,j,1) ; p_int_E(I,j) = p(i+1,j,1)
6750 intx_za_nonlin(I,j) = intx_za(I,j,1) - 0.5*(za(i,j,1) + za(i+1,j,1))
6760 dp_int_x(I,j) = p(i+1,j,1)-p(i,j,1)
6770 seek_x_cor(I,j) = .false.
678 endif ; enddo ; enddo
679 endif
680 endif
681
6820 do j=js,je
6830 do I=Isq,Ieq
684 ! This expression assumes that temperature and salinity vary linearly with pressure
685 ! between the corners of the reference interfaces found above to get a correction to
686 ! intx_pa that takes nonlinearities in the equation of state into account.
687 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale
688 ! linear correction so that the pressures and heights match at the end-points. It turns
689 ! out that this linear correction cancels out the mid-point specific volume.
690 ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land.
6910 T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j)
6920 T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j)
6930 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))
6940 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))
6950 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))
6960 call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref)
697
698 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12
699 intx_za_cor_ri(I,j) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * &
7000 dp_int_x(I,j) - intx_za_nonlin(I,j)
701 enddo
702 enddo
703
704 ! Repeat the calculations above for v-velocity points.
7050 T_int_S(:,:) = 0.0 ; S_int_S(:,:) = 0.0 ; p_int_S(:,:) = 0.0
7060 T_int_N(:,:) = 0.0 ; S_int_N(:,:) = 0.0 ; p_int_N(:,:) = 0.0
7070 inty_za_nonlin(:,:) = 0.0 ; inty_za_cor_ri(:,:) = 0.0 ; dp_int_y(:,:) = 0.0
7080 do J=Jsq,Jeq ; do i=is,ie
7090 seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.)
7100 delta_p_y(i,J) = 0.0
711 enddo ; enddo
712
7130 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
7140 if ((p(i,j+1,2) >= p(i,j,1)) .and. (p(i,j,2) >= p(i,j+1,1))) then
715 ! This is the typical case in the open ocean, so use the topmost interface.
7160 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
7170 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
7180 p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
7190 inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
7200 dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
7210 seek_y_cor(i,J) = .false.
722 endif
723 endif ; enddo ; enddo
724
7250 do k=1,nz
7260 do_more_k = .false.
7270 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
728 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not
729 ! activated in the subgrid interpolation.
7300 if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i,j+1,k) > CS%h_nonvanished)) .and. &
731 (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
732 ! Store properties at the bottom of this cell to get a "good estimate" for intypa at
733 ! the interface below this cell (it might have quadratic pressure dependence if sloped)
7340 T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
7350 S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
7360 p_int_S(i,J) = p(i,j,K+1) ; p_int_N(i,J) = p(i,j+1,K+1)
7370 inty_za_nonlin(i,J) = inty_za(i,J,K+1) - 0.5*(za(i,j,K+1) + za(i,j+1,K+1))
7380 dp_int_y(i,J) = p(i,j+1,K+1) - p(i,j,K+1)
7390 seek_y_cor(i,J) = .false.
740 else
7410 do_more_k = .true.
742 endif
743 endif ; enddo ; enddo
7440 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.
745 enddo
746
7470 if (do_more_k) then
7480 if (CS%reset_intxpa_flattest) then
749 ! There are still points where a correction is needed, so use flattest interface.
7500 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
751 ! choose top interface first
7520 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
7530 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
7540 p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
7550 inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
7560 dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
7570 delta_p_y(i,J) = abs(p(i,j+1,1)-p(i,j,1))
7580 do k=1,nz
7590 if (abs(p(i,j+1,k+1)-p(i,j,k+1)) < delta_p_y(i,J)) then
760 ! bottom of layer is less sloped than top. Use this layer
7610 delta_p_y(i,J) = abs(p(i,j+1,k+1)-p(i,j,k+1))
7620 T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
7630 S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
7640 p_int_S(i,J) = p(i,j,K+1) ; p_int_N(i,J) = p(i,j+1,K+1)
7650 inty_za_nonlin(i,J) = inty_za(i,J,K+1) - 0.5*(za(i,j,K+1) + za(i,j+1,K+1))
7660 dp_int_y(i,J) = p(i,j+1,K+1) - p(i,j,K+1)
767 endif
768 enddo
7690 seek_y_cor(i,J) = .false.
770 endif ; enddo ; enddo
771 else
772 ! There are still points where a correction is needed, so use the top interface.
7730 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
7740 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
7750 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
7760 p_int_S(i,J) = p(i,j,1) ; p_int_N(i,J) = p(i,j+1,1)
7770 inty_za_nonlin(i,J) = inty_za(i,J,1) - 0.5*(za(i,j,1) + za(i,j+1,1))
7780 dp_int_y(i,J) = p(i,j+1,1) - p(i,j,1)
7790 seek_y_cor(i,J) = .false.
780 endif ; enddo ; enddo
781 endif
782 endif
783
7840 do J=Jsq,Jeq
7850 do i=is,ie
786 ! This expression assumes that temperature and salinity vary linearly with pressure
787 ! between the corners of the reference interfaces found above to get a correction to
788 ! intx_pa that takes nonlinearities in the equation of state into account.
789 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale
790 ! linear correction so that the pressures and heights match at the end-points. It turns
791 ! out that this linear correction cancels out the mid-point specific volume.
792 ! This can be used without masking because dp_int_x and intx_za_nonlin are 0 over land.
7930 T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J)
7940 T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J)
7950 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))
7960 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))
7970 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))
7980 call calculate_spec_vol(T5, S5, p5, SpV5, tv%eqn_of_state, spv_ref=alpha_ref)
799
800 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12
801 inty_za_cor_ri(i,J) = C1_90 * (4.75*(SpV5(5)-SpV5(1)) + 5.5*(SpV5(4)-SpV5(2))) * &
8020 dp_int_y(i,J) - inty_za_nonlin(i,J)
803 enddo
804 enddo
805
8060 if (CS%debug) then
807 call uvchksum("Pre-reset int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, &
8080 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
809 call uvchksum("int[xy]_za_cor", intx_za_cor_ri, inty_za_cor_ri, G%HI, haloshift=0, &
8100 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
811 call uvchksum("int[xy]_za_nonlin", intx_za_nonlin, inty_za_nonlin, G%HI, haloshift=0, &
8120 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
813 call uvchksum("dp_int_[xy]", dp_int_x, dp_int_y, G%HI, haloshift=0, &
8140 symmetric=G%Domain%symmetric, unscale=US%RL2_T2_to_Pa)
815 endif
816
817 ! Correct intx_pa and inty_pa at each interface using vertically constant corrections.
8180 do K=1,nz+1 ; do j=js,je ; do I=Isq,Ieq
8190 intx_za(I,j,K) = intx_za(I,j,K) + intx_za_cor_ri(I,j)
820 enddo ; enddo ; enddo
821
8220 do K=1,nz+1 ; do J=Jsq,Jeq ; do i=is,ie
8230 inty_za(i,J,K) = inty_za(i,J,K) + inty_za_cor_ri(i,J)
824 enddo ; enddo ; enddo
825
8260 if (CS%debug) then
827 call uvchksum("Post-reset int[xy]_za", intx_za, inty_za, G%HI, haloshift=0, &
8280 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%L_T_to_m_s**2)
829 endif
830
831 endif ! intx_za and inty_za have now been reset to reflect the properties of an unimpeded interface.
832
833 !$OMP parallel do default(shared) private(dp)
8340 do k=1,nz
8350 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
8360 dp(i,j) = H_to_RL2_T2 * h(i,j,k)
837 enddo ; enddo
838
839 ! Find the horizontal pressure gradient accelerations.
840 ! These expressions for the accelerations have been carefully checked in
841 ! a set of idealized cases, and should be bug-free.
8420 do j=js,je ; do I=Isq,Ieq
843 PFu(I,j,k) = ( ((za(i,j,K+1)*dp(i,j) + intp_dza(i,j,k)) - &
844 (za(i+1,j,K+1)*dp(i+1,j) + intp_dza(i+1,j,k))) + &
845 ((dp(i+1,j) - dp(i,j)) * intx_za(I,j,K+1) - &
846 (p(i+1,j,K) - p(i,j,K)) * intx_dza(I,j,k)) ) * &
8470 (2.0*G%IdxCu(I,j) / ((dp(i,j) + dp(i+1,j)) + dp_neglect))
848 enddo ; enddo
849
8500 do J=Jsq,Jeq ; do i=is,ie
851 PFv(i,J,k) = (((za(i,j,K+1)*dp(i,j) + intp_dza(i,j,k)) - &
852 (za(i,j+1,K+1)*dp(i,j+1) + intp_dza(i,j+1,k))) + &
853 ((dp(i,j+1) - dp(i,j)) * inty_za(i,J,K+1) - &
854 (p(i,j+1,K) - p(i,j,K)) * inty_dza(i,J,k))) * &
8550 (2.0*G%IdyCv(i,J) / ((dp(i,j) + dp(i,j+1)) + dp_neglect))
856 enddo ; enddo
857 enddo
858
8590 if (CS%GFS_scale < 1.0) then
860 ! Adjust the Montgomery potential to make this a reduced gravity model.
8610 if (use_EOS) then
862 !$OMP parallel do default(shared) private(rho_in_situ)
8630 do j=Jsq,Jeq+1
864 call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, &
8650 tv%eqn_of_state, EOSdom)
866
8670 do i=Isq,Ieq+1
8680 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j,1))
869 enddo
870 enddo
871 else
872 !$OMP parallel do default(shared)
8730 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
8740 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/GV%Rlay(1) - alpha_ref) + za(i,j,1))
875 enddo ; enddo
876 endif
877
878 !$OMP parallel do default(shared)
8790 do k=1,nz
8800 do j=js,je ; do I=Isq,Ieq
8810 PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j)
882 enddo ; enddo
8830 do J=Jsq,Jeq ; do i=is,ie
8840 PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IdyCv(i,J)
885 enddo ; enddo
886 enddo
887 endif
888
8890 if (present(pbce)) then
8900 call set_pbce_nonBouss(p, tv_tmp, G, GV, US, CS%GFS_scale, pbce)
891 endif
892
8930 if (present(eta)) then
8940 Pa_to_H = 1.0 / (GV%g_Earth * GV%H_to_RZ)
8950 if (use_p_atm) then
896 !$OMP parallel do default(shared)
8970 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
8980 eta(i,j) = (p(i,j,nz+1) - p_atm(i,j))*Pa_to_H ! eta has the same units as h.
899 enddo ; enddo
900 else
901 !$OMP parallel do default(shared)
9020 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
9030 eta(i,j) = p(i,j,nz+1)*Pa_to_H ! eta has the same units as h.
904 enddo ; enddo
905 endif
906 endif
907
9080 if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag)
9090 if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag)
910
911 ! Diagnostics for tidal forcing and SAL height anomaly
9120 if (CS%id_e_tide>0) then
913 ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.
914 ! New diagnostics are given for each individual field.
9150 if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
9160 e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)
917 enddo ; enddo ; endif
9180 call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag)
919 endif
9200 if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag)
9210 if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag)
9220 if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag)
923
924 ! Diagnostics for tidal forcing and SAL horizontal gradients
9250 if (CS%calculate_SAL .and. (associated(ADp%sal_u) .or. associated(ADp%sal_v))) then
9260 if (CS%tides) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
9270 e_sal(i,j) = e_sal(i,j) + e_tidal_sal(i,j)
928 enddo ; enddo ; endif
9290 if (associated(ADp%sal_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
9300 ADp%sal_u(I,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j)
931 enddo ; enddo ; enddo ; endif
9320 if (associated(ADp%sal_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
9330 ADp%sal_v(i,J,k) = (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J)
934 enddo ; enddo ; enddo ; endif
9350 if (CS%id_sal_u>0) call post_data(CS%id_sal_u, ADp%sal_u, CS%diag)
9360 if (CS%id_sal_v>0) call post_data(CS%id_sal_v, ADp%sal_v, CS%diag)
937 endif
938
9390 if (CS%tides .and. (associated(ADp%tides_u) .or. associated(ADp%tides_v))) then
9400 if (associated(ADp%tides_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
9410 ADp%tides_u(I,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * GV%g_Earth * G%IdxCu(I,j)
942 enddo ; enddo ; enddo ; endif
9430 if (associated(ADp%tides_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
9440 ADp%tides_v(i,J,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * GV%g_Earth * G%IdyCv(i,J)
945 enddo ; enddo ; enddo ; endif
9460 if (CS%id_tides_u>0) call post_data(CS%id_tides_u, ADp%tides_u, CS%diag)
9470 if (CS%id_tides_v>0) call post_data(CS%id_tides_v, ADp%tides_v, CS%diag)
948 endif
9490end subroutine PressureForce_FV_nonBouss
950
951!> \brief Boussinesq analytically-integrated finite volume form of pressure gradient
952!!
953!! Determines the acceleration due to hydrostatic pressure forces, using
954!! the finite volume form of the terms and analytic integrals in depth.
955!!
956!! To work, the following fields must be set outside of the usual (is:ie,js:je)
957!! range before this subroutine is called:
958!! h(isB:ie+1,jsB:je+1), T(isB:ie+1,jsB:je+1), and S(isB:ie+1,jsB:je+1).
95948subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp, p_atm, pbce, eta)
960 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
961 type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
962 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
963 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m]
964 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
965 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: PFu !< Zonal acceleration [L T-2 ~> m s-2]
966 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2]
967 type(PressureForce_FV_CS), intent(in) :: CS !< Finite volume PGF control structure
968 type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
969 type(accel_diag_ptrs), pointer :: ADp !< Acceleration diagnostic pointers
970 real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean
971 !! or atmosphere-ocean interface [R L2 T-2 ~> Pa].
972 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure
973 !! anomaly in each layer due to eta anomalies
974 !! [L2 T-2 H-1 ~> m s-2].
975 real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: eta !< The sea-surface height used to
976 !! calculate PFu and PFv [H ~> m], with any
977 !! tidal contributions.
978 ! Local variables
97948 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: e ! Interface height in depth units [Z ~> m].
980 real, dimension(SZI_(G),SZJ_(G)) :: &
98148 e_sal_and_tide, & ! The summation of self-attraction and loading and tidal forcing [Z ~> m].
98248 e_sal, & ! The bottom geopotential anomaly due to self-attraction and loading [Z ~> m].
98348 e_tidal_eq, & ! The bottom geopotential anomaly due to tidal forces from astronomical sources
984 ! [Z ~> m].
98548 e_tidal_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
986 ! specific to tides [Z ~> m].
98748 Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m]
98848 SSH, & ! Sea surfae height anomaly for self-attraction and loading. Used if
989 ! CALCULATE_SAL is True and SAL_USE_BPA is False [Z ~> m].
99048 pbot, & ! Total bottom pressure for self-attraction and loading. Used if
991 ! CALCULATE_SAL is True and SAL_USE_BPA is True [R L2 T-2 ~> Pa].
99248 dM ! The barotropic adjustment to the Montgomery potential to
993 ! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
994 real, dimension(SZI_(G)) :: &
99548 Rho_cv_BL ! The coordinate potential density in the deepest variable
996 ! density near-surface layer [R ~> kg m-3].
997 real, dimension(SZI_(G),SZJ_(G)) :: &
99848 dz_geo ! The change in geopotential thickness through a layer [L2 T-2 ~> m2 s-2].
999 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
100048 pa ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at the
1001 ! the interface atop a layer [R L2 T-2 ~> Pa].
1002 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
100348 dpa, & ! The change in pressure anomaly between the top and bottom
1004 ! of a layer [R L2 T-2 ~> Pa].
100548 intz_dpa ! The vertical integral in depth of the pressure anomaly less the
1006 ! pressure anomaly at the top of the layer [H R L2 T-2 ~> m Pa].
1007 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
100848 intx_pa ! The zonal integral of the pressure anomaly along the interface
1009 ! atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
1010 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
101148 intx_dpa ! The change in intx_pa through a layer [R L2 T-2 ~> Pa].
1012 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
101348 inty_pa ! The meridional integral of the pressure anomaly along the
1014 ! interface atop a layer, divided by the grid spacing [R L2 T-2 ~> Pa].
1015 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
101648 inty_dpa ! The change in inty_pa through a layer [R L2 T-2 ~> Pa].
1017 real, dimension(SZIB_(G),SZJ_(G)) :: &
101848 intx_pa_cor ! Correction for curvature in intx_pa [R L2 T-2 ~> Pa]
1019 real, dimension(SZI_(G),SZJB_(G)) :: &
102048 inty_pa_cor ! Correction for curvature in inty_pa [R L2 T-2 ~> Pa]
1021
1022 ! These variables are used with reset_intxpa_integral. The values are taken from different
1023 ! interfaces as a function of position.
1024 real, dimension(SZIB_(G),SZJ_(G)) :: &
102548 T_int_W, T_int_E, & ! Temperatures on the reference interface to the east and west of a u-point [C ~> degC]
102648 S_int_W, S_int_E, & ! Salinities on the reference interface to the east and west of a u-point [S ~> ppt]
102748 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]
102848 intx_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface
1029 ! from the value that would be obtained from assuming that pressure varies
1030 ! linearly with depth along that interface [R L2 T-2 ~> Pa].
103148 dgeo_x, & ! The change in x in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2]
103248 intx_pa_cor_ri ! The correction to intx_pa based on the reference interface calculations [R L2 T-2 ~> Pa]
1033 real, dimension(SZI_(G),SZJB_(G)) :: &
103448 T_int_S, T_int_N, & ! Temperatures on the reference interface to the north and south of a v-point [C ~> degC]
103548 S_int_S, S_int_N, & ! Salinities on the reference interface to the north and south of a v-point [S ~> ppt]
103648 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]
103748 inty_pa_nonlin, & ! Deviations in the previous version of intx_pa for the reference interface
1038 ! from the value that would be obtained from assuming that pressure varies
1039 ! linearly with depth along that interface [R L2 T-2 ~> Pa].
104048 dgeo_y, & ! The change in y in geopotenial height along the reference interface [L2 T-2 ~> m2 s-2]
104148 inty_pa_cor_ri ! The correction to inty_pa based on the reference interface calculations [R L2 T-2 ~> Pa]
1042 logical, dimension(SZIB_(G),SZJ_(G)) :: &
104348 seek_x_cor ! If true, try to find a u-point interface that would provide a better estimate
1044 ! of the curvature terms in the intx_pa.
1045 logical, dimension(SZI_(G),SZJB_(G)) :: &
104648 seek_y_cor ! If true, try to find a v-point interface that would provide a better estimate
1047 ! of the curvature terms in the inty_pa.
1048 real, dimension(SZIB_(G),SZJ_(G)) :: &
104948 delta_z_x ! If using flattest interface for reset integral, store x interface differences [Z ~> m]
1050 real, dimension(SZI_(G),SZJB_(G)) :: &
105148 delta_z_y ! If using flattest interface for reset integral, store y interface differences [Z ~> m]
1052
1053 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: &
105448 T_tmp, & ! Temporary array of temperatures where layers that are lighter
1055 ! than the mixed layer have the mixed layer's properties [C ~> degC].
105648 S_tmp ! Temporary array of salinities where layers that are lighter
1057 ! than the mixed layer have the mixed layer's properties [S ~> ppt].
1058 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
105948 S_t, S_b, & ! Top and bottom edge values for linear reconstructions
1060 ! of salinity within each layer [S ~> ppt].
106148 T_t, T_b ! Top and bottom edge values for linear reconstructions
1062 ! of temperature within each layer [C ~> degC].
1063 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
106448 MassWt_u ! The fractional mass weighting at a u-point [nondim].
1065 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
106648 MassWt_v ! The fractional mass weighting at a v-point [nondim].
1067 real, dimension(SZI_(G),SZJ_(G)) :: &
106848 T_top, & ! Temperature of top layer used with correction_intxpa [C ~> degC]
106948 S_top, & ! Salinity of top layer used with correction_intxpa [S ~> ppt]
107048 rho_top ! Density anomaly of top layer used in calculating intx_pa_cor and inty_pa_cor
1071 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
107248 rho_pgf, rho_stanley_pgf ! Density [R ~> kg m-3] from EOS with and without SGS T variance
1073 ! in Stanley parameterization.
1074 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
107548 p_stanley ! Pressure [R L2 T-2 ~> Pa] estimated with Rho_0
107648 real :: zeros(SZI_(G)) ! An array of zero values that can be used as an argument [various]
107748 real :: rho_in_situ(SZI_(G),SZJ_(G)) ! The in situ density [R ~> kg m-3].
107848 real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
1079 ! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
108048 real :: p_surf_EOS(SZI_(G)) ! The pressure at the ocean surface determined from the surface height,
1081 ! consistent with what is used in the density integral routines [R L2 T-2 ~> Pa]
108248 real :: p0(SZI_(G), SZJ_(G)) ! An array of zeros to use for pressure [R L2 T-2 ~> Pa].
1083 real :: dz_geo_sfc ! The change in surface geopotential height between adjacent cells [L2 T-2 ~> m2 s-2]
1084 real :: GxRho0 ! The gravitational acceleration times mean ocean density [R L2 Z-1 T-2 ~> Pa m-1]
1085 real :: GxRho_ref ! The gravitational acceleration times reference density [R L2 Z-1 T-2 ~> Pa m-1]
1086 real :: rho0_int_density ! Rho0 used in int_density_dz_* subroutines [R ~> kg m-3]
1087 real :: rho0_set_pbce ! Rho0 used in set_pbce_Bouss subroutine [R ~> kg m-3]
1088 real :: h_neglect ! A thickness that is so small it is usually lost
1089 ! in roundoff and can be neglected [H ~> m].
1090 real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
1091 real :: G_Rho0 ! G_Earth / Rho_0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
1092 real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]
1093 real :: rho_ref ! The reference density [R ~> kg m-3].
1094 real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
1095 real :: dz_nonvanished ! A small thickness considered to be vanished for mass weighting [Z ~> m]
1096 real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure
1097 ! units [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1].
1098 real :: T5(5) ! Temperatures and salinities at five quadrature points [C ~> degC]
1099 real :: S5(5) ! Salinities at five quadrature points [S ~> ppt]
1100 real :: p5(5) ! Full pressures at five quadrature points for use with the equation of state [R L2 T-2 ~> Pa]
1101 real :: pa5(5) ! The pressure anomaly (i.e. pressure + g*RHO_0*e) at five quadrature points [R L2 T-2 ~> Pa].
1102 real :: r5(5) ! Densities at five quadrature points [R ~> kg m-3]
1103 real :: wt_R ! A weighting factor [nondim]
1104 real, parameter :: C1_6 = 1.0/6.0 ! A rational constant [nondim]
1105 real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]
1106 logical :: use_p_atm ! If true, use the atmospheric pressure.
1107 logical :: use_ALE ! If true, use an ALE pressure reconstruction.
1108 logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
1109 logical :: do_more_k ! If true, there are still points where a flatter interface remains to be found.
111024 type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
1111 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1112 integer, dimension(2) :: EOSdom_h ! The i-computational domain for the equation of state at tracer points
1113 integer, dimension(2) :: EOSdom_u ! The i-computational domain for the equation of state at u-velocity points
1114 integer, dimension(2) :: EOSdom_v ! The i-computational domain for the equation of state at v-velocity points
1115 integer :: EOSdom2d(2,2) ! The 2D compute domain for the equation of state
1116 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
1117 integer :: i, j, k, m, kstart, kend
1118 integer :: nkblock
1119
112024 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
112124 nkblock = CS%nkblock ; if (nkblock <= 0) nkblock = nz
112224 nkmb=GV%nk_rho_varies
112324 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
112424 EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)
112524 EOSdom_u(1) = Isq - (G%IsdB-1) ; EOSdom_u(2) = Ieq - (G%IsdB-1)
112624 EOSdom_v(1) = is - (G%isd-1) ; EOSdom_v(2) = ie - (G%isd-1)
1127
112872 EOSdom2d(1,:) = EOSdom(:)
112972 EOSdom2d(2,:) = [Jsq - (G%isd - 1), (je + 1) - (G%jsd - 1)]
1130
1131 ! TODO: This would be done outside of the function!
1132
113324 if (.not.CS%initialized) call MOM_error(FATAL, &
11340 "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")
1135
113624 use_p_atm = associated(p_atm)
113724 use_EOS = associated(tv%eqn_of_state)
113824 use_ALE = .false.
113924 if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS
1140
114124 h_neglect = GV%H_subroundoff
114224 dz_neglect = GV%dZ_subroundoff
114324 dz_nonvanished = GV%H_to_Z*CS%h_nonvanished
114424 I_Rho0 = 1.0 / GV%Rho0
114524 G_Rho0 = GV%g_Earth / GV%Rho0
114624 GxRho0 = GV%g_Earth * GV%Rho0
114724 rho_ref = CS%rho_ref
1148
114924 if (CS%rho_ref_bug) then
115024 rho0_int_density = rho_ref
115124 rho0_set_pbce = rho_ref
115224 GxRho_ref = GxRho0
115324 I_g_rho = 1.0 / (rho_ref * GV%g_Earth)
1154 else
11550 rho0_int_density = GV%Rho0
11560 rho0_set_pbce = GV%Rho0
11570 GxRho_ref = GV%g_Earth * rho_ref
11580 I_g_rho = 1.0 / (GV%rho0 * GV%g_Earth)
1159 endif
1160
116124 if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then
11620 MassWt_u(:,:,:) = 0.0 ; MassWt_v(:,:,:) = 0.0
1163 endif
1164
1165 !$omp target enter data map(alloc: e, T_t, T_b, S_t, S_b)
1166
11672952 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
1168184488 e(i,j,nz+1) = -G%bathyT(i,j)
1169 enddo
1170
1171 ! The following two if-blocks are used to recover old answers for self-attraction and loading
1172 ! (SAL) and tides only. The old algorithm moves interface heights before density calculations,
1173 ! and therefore is incorrect without SSH_IN_EOS_PRESSURE_FOR_PGF=True (added in August 2024).
1174 ! See the code right after Pa calculation loop for the new algorithm.
1175
1176 ! Calculate and add SAL geopotential anomaly to interface height (old answers)
117724 if (CS%calculate_SAL .and. CS%tides_answer_date<=20250131) then
1178 !$omp target update from(e(:,:,nz+1))
1179
1180 !$OMP parallel do default(shared)
11810 do j=Jsq,Jeq+1
11820 do i=Isq,Ieq+1
11830 SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0)
1184 enddo
11850 do k=1,nz ; do i=Isq,Ieq+1
11860 SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
1187 enddo ; enddo
1188 enddo
11890 call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
1190
11910 if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]
1192 !$OMP parallel do default(shared)
11930 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
11940 e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
1195 enddo ; enddo
1196 endif
1197
1198 !$omp target update to(e(:,:,nz+1))
1199 endif
1200
1201 ! Calculate and add tidal geopotential anomaly to interface height (old answers)
120224 if (CS%tides .and. CS%tides_answer_date<=20250131) then
1203 !$omp target update from(e(:,:,nz+1))
1204
12050 if (CS%tides_answer_date>20230630) then ! answers_date between [20230701, 20250131]
12060 call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
1207 !$OMP parallel do default(shared)
12080 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
12090 e(i,j,nz+1) = e(i,j,nz+1) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1210 enddo ; enddo
1211 else ! answers_date before 20230701
12120 if (.not.CS%calculate_SAL) e_sal(:,:) = 0.0
1213 call calc_tidal_forcing_legacy(CS%Time, e_sal, e_sal_and_tide, e_tidal_eq, e_tidal_sal, &
12140 G, US, CS%tides_CSp)
1215 !$OMP parallel do default(shared)
12160 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
12170 e(i,j,nz+1) = e(i,j,nz+1) - e_sal_and_tide(i,j)
1218 enddo ; enddo
1219 endif
1220
1221 !$omp target update to(e(:,:,nz+1))
1222 endif
1223
122424 do concurrent(j=Jsq:Jeq+1)
1225113112 do k=nz,1,-1
12261488 do concurrent (i=Isq:Ieq+1)
122713726800 e(i,j,K) = e(i,j,K+1) + h(i,j,k)*GV%H_to_Z
1228 enddo
1229 enddo
1230 enddo
1231
123224 if (use_EOS) then
123324 if (nkmb>0) then
1234 ! With a bulk mixed layer, replace the T & S of any layers that are lighter than the buffer
1235 ! layer with the properties of the buffer layer. These layers will be massless anyway, and
1236 ! it avoids any formal calculations with hydrostatically unstable profiles.
12370 tv_tmp%T => T_tmp ; tv_tmp%S => S_tmp
12380 tv_tmp%eqn_of_state => tv%eqn_of_state
1239
12400 do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
1241 !$OMP parallel do default(shared) private(Rho_cv_BL)
12420 do j=Jsq,Jeq+1
12430 do k=1,nkmb ; do i=Isq,Ieq+1
12440 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
1245 enddo ; enddo
1246 call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
12470 tv%eqn_of_state, EOSdom)
1248
12490 do k=nkmb+1,nz ; do i=Isq,Ieq+1
12500 if (GV%Rlay(k) < Rho_cv_BL(i)) then
12510 tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
1252 else
12530 tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
1254 endif
1255 enddo ; enddo
1256 enddo
1257 else
125824 tv_tmp%T => tv%T ; tv_tmp%S => tv%S
125924 tv_tmp%eqn_of_state => tv%eqn_of_state
1260 endif
1261 !$omp target enter data map(to: tv_tmp, tv_tmp%T, tv_tmp%S)
1262 endif
1263
1264 ! If regridding is activated, do a linear reconstruction of salinity
1265 ! and temperature across each layer. The subscripts 't' and 'b' refer
1266 ! to top and bottom values within each layer (these are the only degrees
1267 ! of freedom needed to know the linear profile).
126824 if ( use_ALE .and. (CS%Recon_Scheme == 1) ) then
126924 call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap)
12700 elseif ( use_ALE .and. (CS%Recon_Scheme == 2) ) then
12710 call TS_PPM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, CS%boundary_extrap)
12720 elseif ( use_ALE .and. (CS%Recon_Scheme == 3) ) then
12730 call TS_PLM_WLS_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h)
12740 elseif (CS%reset_intxpa_integral) then
12750 do k=1,nz ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
12760 T_b(i,j,k) = tv%T(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k)
1277 enddo ; enddo ; enddo
1278 endif
1279
1280 !$omp target enter data map(alloc: pa)
1281
1282 ! Set the surface boundary conditions on pressure anomaly and its horizontal
1283 ! integrals, assuming that the surface pressure anomaly varies linearly
1284 ! in x and y.
128524 if (use_p_atm) then
12862952 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
1287184488 pa(i,j,1) = GxRho_ref * (e(i,j,1) - G%Z_ref) + p_atm(i,j)
1288 enddo
1289 else
12900 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
12910 pa(i,j,1) = GxRho_ref * (e(i,j,1) - G%Z_ref)
1292 enddo
1293 endif
1294
129524 if (use_EOS) then
1296 !$omp target enter data map(alloc: Z_0p)
129724 if (CS%use_SSH_in_Z0p .and. use_p_atm) then
12980 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
12990 Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho
1300 enddo
130124 elseif (CS%use_SSH_in_Z0p) then
13020 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
13030 Z_0p(i,j) = e(i,j,1)
1304 enddo
1305 else
13062952 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
1307184488 Z_0p(i,j) = G%meanSL(i,j)
1308 enddo
1309 endif
1310 !$omp target update from(Z_0p) &
1311 !$omp if((use_ALE .and. CS%Recon_Scheme == 2) .or. CS%reset_intxpa_integral .or. CS%correction_intxpa)
1312 endif
1313
1314 ! Calculate 4 integrals through the layer that are required in the
1315 ! subsequent calculation.
1316 !$omp target enter data map(alloc: dpa, intx_dpa, inty_dpa, intz_dpa)
1317
131824 if (use_EOS) then
1319 !$omp target update from(e) if( (use_ALE .and. CS%Recon_Scheme == 2) .or. &
1320 !$omp (CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0))
1321
1322 ! The following routine computes the integrals that are needed to
1323 ! calculate the pressure gradient force. Linear profiles for T and S are
1324 ! assumed when regridding is activated. Otherwise, the previous version
1325 ! is used, whereby densities within each layer are constant no matter
1326 ! where the layers are located.
1327 !$omp target enter data map(to: tv)
13281824 do kstart=1,nz,nkblock
13291800 kend = min(kstart+nkblock-1, nz)
13301800 if ( use_ALE .and. CS%Recon_Scheme > 0 ) then
13311800 if ( CS%Recon_Scheme == 1 .or. CS%Recon_Scheme == 3 ) then
1332 call int_density_dz_generic_plm(kstart, kend, tv, T_t, T_b, S_t, S_b, e, &
1333 rho_ref, rho0_int_density, GV%g_Earth, dz_neglect, G%bathyT, &
1334 G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, &
1335 intx_dpa, inty_dpa, &
1336 MassWghtInterp=CS%MassWghtInterp, &
1337 use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=Z_0p, &
1338 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished, &
13391800 niblock=CS%niblock_plm, njblock=CS%njblock_plm)
13400 elseif ( CS%Recon_Scheme == 2 ) then
13410 do k=kstart,kend
1342 call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
1343 rho_ref, rho0_int_density, GV%g_Earth, dz_neglect, G%bathyT, &
1344 G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
1345 intx_dpa(:,:,k), inty_dpa(:,:,k), &
1346 MassWghtInterp=CS%MassWghtInterp, Z_0p=Z_0p, &
13470 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished)
1348 enddo
1349 endif
1350 else
13510 do k=kstart,kend
1352 call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
1353 rho_ref, rho0_int_density, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), &
1354 intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, e(:,:,1), dz_neglect, &
1355 CS%MassWghtInterp, Z_0p=Z_0p, &
13560 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=dz_nonvanished)
1357 enddo
1358 endif
13591800 if (GV%Z_to_H /= 1.0) then
13600 do concurrent (k=kstart:kend, j=Jsq:Jeq+1, i=Isq:Ieq+1)
13610 intz_dpa(i,j,k) = intz_dpa(i,j,k)*GV%Z_to_H
1362 enddo
1363 endif
13641824 if ((CS%id_MassWt_u > 0) .or. (CS%id_MassWt_v > 0)) then
13650 do k=kstart,kend
1366 call diagnose_mass_weight_Z(e(:,:,K), e(:,:,K+1), G%bathyT, e(:,:,1), dz_neglect, CS%MassWghtInterp, &
1367 G%HI, MassWt_u(:,:,k), MassWt_v(:,:,k), &
13680 MassWghtInterpVanOnly=CS%MassWghtInterpVanOnly, h_nv=CS%h_nonvanished)
1369 enddo
1370 endif
1371 enddo
1372 !$omp target exit data map(release: tv)
1373 else
1374 !$omp target data map(alloc: dz_geo)
13750 do k=1,nz
13760 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
13770 dz_geo(i,j) = GV%g_Earth * GV%H_to_Z*h(i,j,k)
13780 dpa(i,j,k) = (GV%Rlay(k) - rho_ref) * dz_geo(i,j)
13790 intz_dpa(i,j,k) = 0.5*(GV%Rlay(k) - rho_ref) * dz_geo(i,j)*h(i,j,k)
1380 enddo
13810 do concurrent (j=js:je, I=Isq:Ieq)
13820 intx_dpa(I,j,k) = 0.5*(GV%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i+1,j))
1383 enddo
13840 do concurrent (J=Jsq:Jeq, i=is:ie)
13850 inty_dpa(i,J,k) = 0.5*(GV%Rlay(k) - rho_ref) * (dz_geo(i,j) + dz_geo(i,j+1))
1386 enddo
1387 enddo
1388 !$omp end target data
1389 endif
1390
1391 ! Set the pressure anomalies at the interfaces.
139224 do concurrent (j=Jsq:Jeq+1)
1393113112 do k=1,nz
13941488 do concurrent (i=Isq:Ieq+1)
139513726800 pa(i,j,K+1) = pa(i,j,K) + dpa(i,j,k)
1396 enddo
1397 enddo
1398 enddo
1399
1400 ! Calculate and add SAL geopotential anomaly to interface height (new answers)
140124 if (CS%calculate_SAL .and. CS%tides_answer_date>20250131) then
14020 if (CS%sal_use_bpa) then
1403 !$OMP parallel do default(shared)
14040 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14050 pbot(i,j) = pa(i,j,nz+1) - GxRho_ref * (e(i,j,nz+1) - G%Z_ref)
1406 enddo ; enddo
14070 call calc_SAL(pbot, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
1408 else
1409 !$OMP parallel do default(shared)
14100 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14110 SSH(i,j) = e(i,j,1) - G%Z_ref
1412 ! Remove above sea level topography at floodable cells
14130 SSH(i,j) = SSH(i,j) - max(-G%bathyT(i,j)-G%meanSL(i,j), 0.0)
1414 enddo ; enddo
14150 call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
1416 endif
14170 if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
1418 !$OMP parallel do default(shared)
14190 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14200 e(i,j,K) = e(i,j,K) - e_sal(i,j)
14210 pa(i,j,K) = pa(i,j,K) - GxRho_ref * e_sal(i,j)
1422 enddo ; enddo
1423 enddo ; endif
1424 endif
1425
1426 ! Calculate and add tidal geopotential anomaly to interface height (new answers)
142724 if (CS%tides .and. CS%tides_answer_date>20250131) then
14280 call calc_tidal_forcing(CS%Time, e_tidal_eq, e_tidal_sal, G, US, CS%tides_CSp)
14290 if (.not.CS%bq_sal_tides) then ; do K=1,nz+1
1430 !$OMP parallel do default(shared)
14310 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14320 e(i,j,K) = e(i,j,K) - (e_tidal_eq(i,j) + e_tidal_sal(i,j))
14330 pa(i,j,K) = pa(i,j,K) - GxRho_ref * (e_tidal_eq(i,j) + e_tidal_sal(i,j))
1434 enddo ; enddo
1435 enddo ; endif
1436 endif
1437
143824 if (CS%correction_intxpa .or. CS%reset_intxpa_integral) then
1439 ! Determine surface temperature and salinity for use in the pressure gradient corrections
14400 if (use_ALE .and. (CS%Recon_Scheme > 0)) then
14410 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14420 T_top(i,j) = T_t(i,j,1) ; S_top(i,j) = S_t(i,j,1)
1443 enddo ; enddo
1444 else
14450 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
14460 T_top(i,j) = tv%T(i,j,1) ; S_top(i,j) = tv%S(i,j,1)
1447 enddo ; enddo
1448 endif
1449 endif
1450
1451 !$omp target enter data map(alloc: intx_pa, inty_pa)
145224 if (CS%correction_intxpa) then
1453 ! TODO needs to be moved to GPU
1454
1455 ! Determine surface density for use in the pressure gradient corrections
1456 !$OMP parallel do default(shared) private(p_surf_EOS)
14570 do j=Jsq,Jeq+1
1458 ! P_surf_EOS here is consistent with the pressure that is used in the int_density_dz routines.
14590 do i=Isq,Ieq+1 ; p_surf_EOS(i) = -GxRho0*(e(i,j,1) - Z_0p(i,j)) ; enddo
1460 call calculate_density(T_top(:,j), S_top(:,j), p_surf_EOS, rho_top(:,j), &
14610 tv%eqn_of_state, EOSdom, rho_ref=rho_ref)
1462 enddo
1463
14640 if (CS%debug) then
14650 call hchksum(rho_top, "intx_pa rho_top", G%HI, haloshift=1, unscale=US%R_to_kg_m3)
14660 call hchksum(e(:,:,1), "intx_pa e(1)", G%HI, haloshift=1, unscale=US%Z_to_m)
14670 call hchksum(pa(:,:,1), "intx_pa pa(1)", G%HI, haloshift=1, unscale=US%RL2_T2_to_Pa)
1468 endif
1469
1470 ! This version attempts to correct for hydrostatic variations in surface pressure under ice.
1471 !$OMP parallel do default(shared) private(dz_geo_sfc)
14720 do j=js,je ; do I=Isq,Ieq
14730 intx_pa_cor(I,j) = 0.0
14740 dz_geo_sfc = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
14750 if ((dz_geo_sfc * rho_ref - (pa(i+1,j,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then
1476 ! The pressure/depth relationship has a positive implied density given by
1477 ! rho_implied = rho_ref - (pa(i+1,j,1)-pa(i,j,1)) / dz_geo_sfc
14780 if (-dz_geo_sfc * (pa(i+1,j,1)-pa(i,j,1)) > &
1479 0.25*((rho_top(i+1,j)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then
1480 ! The pressure difference is at least half the size of the difference expected by hydrostatic
1481 ! balance. This test gets rid of pressure differences that are small, e.g. open ocean.
1482 ! Use 5 point quadrature to calculate intxpa
14830 T5(1) = T_top(i,j) ; T5(5) = T_top(i+1,j)
14840 S5(1) = S_top(i,j) ; S5(5) = S_top(i+1,j)
14850 pa5(1) = pa(i,j,1) ; pa5(5) = pa(i+1,j,1)
1486 ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines.
14870 p5(1) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
14880 p5(5) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
14890 do m=2,4
14900 wt_R = 0.25*real(m-1)
14910 T5(m) = T5(1) + (T5(5)-T5(1))*wt_R
14920 S5(m) = S5(1) + (S5(5)-S5(1))*wt_R
14930 p5(m) = p5(1) + (p5(5)-p5(1))*wt_R
1494 enddo !m
14950 call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
1496
1497 ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure
1498 ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule
1499 ! quadrature to find the integrated correction to the integral of pressure along the interface.
1500 ! The derivation for this expression is shown below in the y-direction version.
15010 intx_pa_cor(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc
1502 ! Note that (4.75 + 5.5/2) / 90 = 1/12, so this is consistent with the linear result below.
1503 endif
1504 endif
15050 intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1)) + intx_pa_cor(I,j)
1506 enddo ; enddo
1507 !$OMP parallel do default(shared) private(dz_geo_sfc)
15080 do J=Jsq,Jeq ; do i=is,ie
15090 inty_pa_cor(i,J) = 0.0
15100 dz_geo_sfc = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
15110 if ((dz_geo_sfc * rho_ref - (pa(i,j+1,1)-pa(i,j,1)))*dz_geo_sfc > 0.0) then
1512 ! The pressure/depth relationship has a positive implied density
15130 if (-dz_geo_sfc * (pa(i,j+1,1)-pa(i,j,1)) > &
1514 0.25*((rho_top(i,j+1)+rho_top(i,j))-2.0*rho_ref) * dz_geo_sfc**2) then
1515 ! The pressure difference is at least half the size of the difference expected by hydrostatic
1516 ! balance. This test gets rid of pressure differences that are small, e.g. open ocean.
1517 ! Use 5 point quadrature to calculate intypa
15180 T5(1) = T_top(i,j) ; T5(5) = T_top(i,j+1)
15190 S5(1) = S_top(i,j) ; S5(5) = S_top(i,j+1)
15200 pa5(1) = pa(i,j,1) ; pa5(5) = pa(i,j+1,1)
1521 ! Pressure input to density EOS is consistent with the pressure used in the int_density_dz routines.
15220 p5(1) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
15230 p5(5) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
1524
15250 do m=2,4
15260 wt_R = 0.25*real(m-1)
15270 T5(m) = T5(1) + (T5(5)-T5(1))*wt_R
15280 S5(m) = S5(1) + (S5(5)-S5(1))*wt_R
15290 p5(m) = p5(1) + (p5(5)-p5(1))*wt_R
1530 enddo !m
15310 call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
1532
1533 ! Use a trapezoidal rule integral of the hydrostatic equation to determine the pressure
1534 ! anomalies at 5 equally spaced points along the interface, and then use Boole's rule
1535 ! quadrature to find the integrated correction to the integral of pressure along the interface.
15360 inty_pa_cor(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dz_geo_sfc
1537
1538 ! The derivation of this correction follows:
1539
1540 ! Make pressure curvature a difference from the linear fit of pressure between the two points
1541 ! (which is equivalent to taking 4 trapezoidal rule integrals of the hydrostatic equation on
1542 ! sub-segments), with a constant slope that is chosen so that the pressure anomalies at the
1543 ! two ends of the segment agree with their known values.
1544 ! d_geo_8 = 0.125*dz_geo_sfc
1545 ! dpa_subseg = 0.25*(pa5(5)-pa5(1)) + &
1546 ! 0.25*d_geo_8 * ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)))
1547 ! do m=2,4
1548 ! pa5(m) = pa5(m-1) + dpa_subseg - d_geo_8*(r5(m)+r5(m-1)))
1549 ! enddo
1550
1551 ! Explicitly finding expressions for the incremental pressures from the recursion relation above:
1552 ! 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)) )
1553 ! ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + 0.25*d_geo_8 * &
1554 ! ! ( (r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3)) + &
1555 ! ! (r5(5)-3.*r5(1)) + 2.0*((r5(4)-r5(2)) + r5(3)) - 4.*(r5(3)+r5(2)) )
1556 ! pa5(3) = 0.5*(pa5(1) + pa5(5)) + d_geo_8 * (0.5*(r5(5)-r5(1)) + (r5(4)-r5(2)) )
1557 ! ! pa5(4) = 0.25*(pa5(1) + 3.0*pa5(5)) + 0.25*d_geo_8 * &
1558 ! ! (2.0*(r5(5)-r5(1)) + 4.0*(r5(4)-r5(2)) + (r5(5)+r5(1)) + &
1559 ! ! 2.0*(r5(4)+r5(2)) + 2.0*r5(3) - 4.*(r5(4)+r5(3)))
1560 ! 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)) )
1561 ! ! pa5(5) = pa5(5) + 0.25*d_geo_8 * &
1562 ! ! ( (3.*r5(5)-r5(1)) + 2.0*((r5(4)-r5(2)) - r5(3)) + &
1563 ! ! ((r5(5)+r5(1)) + 2.0*((r5(4)+r5(2)) + r5(3))) - 4.*(r5(5)+r5(4)) )
1564 ! pa5(5) = pa5(5) ! As it should.
1565
1566 ! From these:
1567 ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + 0.25*d_geo_8 * &
1568 ! ( (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))
1569 ! pa5(2) + pa5(4) = (pa5(1) + pa5(5)) + d_geo_8 * ( (r5(5)-r5(1)) + (r5(4)-r5(2)) )
1570
1571 ! Get the correction from the difference between the 5-point quadrature integral of pa5 and
1572 ! its trapezoidal rule integral as:
1573 ! 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)))
1574 ! inty_pa_cor(i,J) = C1_90*((32.0*(pa5(2)+pa5(4)) + 12.0*pa5(3)) - 38.0*(pa5(1)+pa5(5)))
1575 ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ((32.0*( (r5(5)-r5(1)) + (r5(4)-r5(2)) ) + &
1576 ! (6.*(r5(5)-r5(1)) + 12.0*(r5(4)-r5(2)) ))
1577 ! inty_pa_cor(i,J) = C1_90*d_geo_8 * ( 38.0*(r5(5)-r5(1)) + 44.0*(r5(4)-r5(2)) )
1578 endif
1579 endif
15800 inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1)) + inty_pa_cor(i,J)
1581 enddo ; enddo
1582
15830 if (CS%debug) then
1584 call uvchksum("int[xy]_pa_cor", intx_pa_cor, inty_pa_cor, G%HI, haloshift=0, &
15850 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%RL2_T2_to_Pa)
1586 call uvchksum("int[xy]_pa(1)", intx_pa(:,:,1), inty_pa(:,:,1), G%HI, haloshift=0, &
15870 symmetric=G%Domain%symmetric, scalar_pair=.true., unscale=US%RL2_T2_to_Pa)
1588 endif
1589
1590 ! TODO until moved to GPU
1591 !$omp target update to(intx_pa, inty_pa)
1592 else
1593 ! Set the surface boundary conditions on the horizontally integrated pressure anomaly,
1594 ! assuming that the surface pressure anomaly varies linearly in x and y.
1595 ! If there is an ice-shelf or icebergs, this linear variation would need to be applied
1596 ! to an interior interface.
159724 do concurrent (j=js:je, I=Isq:Ieq)
1598177168 intx_pa(I,j,1) = 0.5*(pa(i,j,1) + pa(i+1,j,1))
1599 enddo
160024 do concurrent (J=Jsq:Jeq, i=is:ie)
1601178584 inty_pa(i,J,1) = 0.5*(pa(i,j,1) + pa(i,j+1,1))
1602 enddo
1603 endif
1604
160524 do concurrent (j=js:je)
1606109464 do k=1,nz
16071440 do concurrent (I=Isq:Ieq)
160813176000 intx_pa(I,j,K+1) = intx_pa(I,j,K) + intx_dpa(I,j,k)
1609 enddo
1610 enddo
1611 enddo
161224 do concurrent (J=Jsq:Jeq)
1613111288 do k=1,nz
16141464 do concurrent (i=is:ie)
161513285800 inty_pa(i,J,K+1) = inty_pa(i,J,K) + inty_dpa(i,J,k)
1616 enddo
1617 enddo
1618 enddo
1619
162024 if (CS%reset_intxpa_integral) then
1621 ! Having stored the pressure gradient info, we can work out where the first nonvanished layers is
1622 ! reset intxpa there, then adjust intxpa throughout the water column.
1623
1624 ! TODO temporarily move back to CPU
1625 !$omp target update from(intx_pa, inty_pa)
1626
1627 ! Zero out the 2-d arrays that will be set from various reference interfaces.
16280 T_int_W(:,:) = 0.0 ; S_int_W(:,:) = 0.0 ; p_int_W(:,:) = 0.0
16290 T_int_E(:,:) = 0.0 ; S_int_E(:,:) = 0.0 ; p_int_E(:,:) = 0.0
16300 intx_pa_nonlin(:,:) = 0.0 ; dgeo_x(:,:) = 0.0 ; intx_pa_cor_ri(:,:) = 0.0
16310 do j=js,je ; do I=Isq,Ieq
16320 seek_x_cor(I,j) = (G%mask2dCu(I,j) > 0.)
16330 delta_z_x(I,j) = 0.0
1634 enddo ; enddo
1635
16360 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
16370 if ((e(i+1,j,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i+1,j,1))) then
1638 ! This is a typical case in the open ocean, so use the topmost interface.
16390 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
16400 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
16410 p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
16420 p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
16430 intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
16440 dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
16450 seek_x_cor(I,j) = .false.
1646 endif
1647 endif ; enddo ; enddo
1648
16490 do k=1,nz
16500 do_more_k = .false.
16510 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
1652 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not
1653 ! activated in the subgrid interpolation.
16540 if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i+1,j,k) > CS%h_nonvanished)) .and. &
1655 (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
1656 ! Store properties at the bottom of this cell to get a "good estimate" for intxpa at
1657 ! the interface below this cell (it might have quadratic pressure dependence if sloped)
16580 T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
16590 S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
1660 ! These pressures are only used for the equation of state, and are only a function of
1661 ! height, consistent with the expressions in the int_density_dz routines.
16620 p_int_W(I,j) = -GxRho0*(e(i,j,K+1) - Z_0p(i,j))
16630 p_int_E(I,j) = -GxRho0*(e(i+1,j,K+1) - Z_0p(i,j))
1664
16650 intx_pa_nonlin(I,j) = intx_pa(I,j,K+1) - 0.5*(pa(i,j,K+1) + pa(i+1,j,K+1))
16660 dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,K+1)-e(i,j,K+1))
16670 seek_x_cor(I,j) = .false.
1668 else
16690 do_more_k = .true.
1670 endif
1671 endif ; enddo ; enddo
16720 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.
1673 enddo
1674
16750 if (do_more_k) then
16760 if (CS%reset_intxpa_flattest) then
1677 ! There are still points where a correction is needed, so use flattest interface
16780 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
1679 ! choose top layer first
16800 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
16810 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
16820 p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
16830 p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
16840 intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
16850 dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
16860 delta_z_x(I,j) = abs(e(i+1,j,1)-e(i,j,1))
16870 do k=1,nz
16880 if (abs(e(i+1,j,k+1)-e(i,j,k+1)) < delta_z_x(I,j)) then
1689 ! bottom of layer is less sloped than top. Use this layer
16900 delta_z_x(I,j) = abs(e(i+1,j,k+1)-e(i,j,k+1))
16910 T_int_W(I,j) = T_b(i,j,k) ; T_int_E(I,j) = T_b(i+1,j,k)
16920 S_int_W(I,j) = S_b(i,j,k) ; S_int_E(I,j) = S_b(i+1,j,k)
16930 p_int_W(I,j) = -GxRho0*(e(i,j,K+1) - Z_0p(i,j))
16940 p_int_E(I,j) = -GxRho0*(e(i+1,j,K+1) - Z_0p(i,j))
16950 intx_pa_nonlin(I,j) = intx_pa(I,j,K+1) - 0.5*(pa(i,j,K+1) + pa(i+1,j,K+1))
16960 dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,K+1)-e(i,j,K+1))
1697 endif
1698 enddo
16990 seek_x_cor(I,j) = .false.
1700 endif ; enddo ; enddo
1701 else
1702 ! There are still points where a correction is needed, so use the top interface for lack of a better idea?
17030 do j=js,je ; do I=Isq,Ieq ; if (seek_x_cor(I,j)) then
17040 T_int_W(I,j) = T_top(i,j) ; T_int_E(I,j) = T_top(i+1,j)
17050 S_int_W(I,j) = S_top(i,j) ; S_int_E(I,j) = S_top(i+1,j)
17060 p_int_W(I,j) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
17070 p_int_E(I,j) = -GxRho0*(e(i+1,j,1) - Z_0p(i,j))
17080 intx_pa_nonlin(I,j) = intx_pa(I,j,1) - 0.5*(pa(i,j,1) + pa(i+1,j,1))
17090 dgeo_x(I,j) = GV%g_Earth * (e(i+1,j,1)-e(i,j,1))
17100 seek_x_cor(I,j) = .false.
1711 endif ; enddo ; enddo
1712 endif
1713 endif
1714
17150 do j=js,je
17160 do I=Isq,Ieq
1717 ! This expression assumes that temperature and salinity vary linearly with hieght
1718 ! between the corners of the reference interfaces found above to get a correction to
1719 ! intx_pa that takes nonlinearities in the equation of state into account.
1720 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale
1721 ! linear correction so that the pressures and heights match at the end-points. It turns
1722 ! out that this linear correction cancels out the mid-point density anomaly.
1723 ! This can be used without masking because dgeo_x and intx_pa_nonlin are 0 over land.
17240 T5(1) = T_Int_W(I,j) ; S5(1) = S_Int_W(I,j) ; p5(1) = p_Int_W(I,j)
17250 T5(5) = T_Int_E(I,j) ; S5(5) = S_Int_E(I,j) ; p5(5) = p_Int_E(I,j)
17260 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))
17270 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))
17280 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))
17290 call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
1730
1731 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12
1732 intx_pa_cor_ri(I,j) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_x(I,j) - &
17330 intx_pa_nonlin(I,j)
1734 enddo
1735 enddo
1736
1737 ! Repeat the calculations above for v-velocity points.
17380 T_int_S(:,:) = 0.0 ; S_int_S(:,:) = 0.0 ; p_int_S(:,:) = 0.0
17390 T_int_N(:,:) = 0.0 ; S_int_N(:,:) = 0.0 ; p_int_N(:,:) = 0.0
17400 inty_pa_nonlin(:,:) = 0.0 ; dgeo_y(:,:) = 0.0 ; inty_pa_cor_ri(:,:) = 0.0
17410 do J=Jsq,Jeq ; do i=is,ie
17420 seek_y_cor(i,J) = (G%mask2dCv(i,J) > 0.)
17430 delta_z_y(i,J) = 0.0
1744 enddo ; enddo
1745
17460 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
17470 if ((e(i,j+1,2) <= e(i,j,1)) .and. (e(i,j,2) <= e(i,j+1,1))) then
1748 ! This is a typical case in the open ocean, so use the topmost interface.
17490 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
17500 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
17510 p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
17520 p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
17530 inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
17540 dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
17550 seek_y_cor(i,J) = .false.
1756 endif
1757 endif ; enddo ; enddo
1758
17590 do k=1,nz
17600 do_more_k = .false.
17610 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
1762 ! Find the topmost layer for which both sides are nonvanished and mass-weighting is not
1763 ! activated in the subgrid interpolation.
17640 if (((h(i,j,k) > CS%h_nonvanished) .and. (h(i,j+1,k) > CS%h_nonvanished)) .and. &
1765 (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
1766 ! Store properties at the bottom of this cell to get a "good estimate" for intypa at
1767 ! the interface below this cell (it might have quadratic pressure dependence if sloped)
17680 T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
17690 S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
1770 ! These pressures are only used for the equation of state, and are only a function of
1771 ! height, consistent with the expressions in the int_density_dz routines.
17720 p_int_S(i,J) = -GxRho0*(e(i,j,K+1) - Z_0p(i,j))
17730 p_int_N(i,J) = -GxRho0*(e(i,j+1,K+1) - Z_0p(i,j))
17740 inty_pa_nonlin(i,J) = inty_pa(i,J,K+1) - 0.5*(pa(i,j,K+1) + pa(i,j+1,K+1))
17750 dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,K+1)-e(i,j,K+1))
17760 seek_y_cor(i,J) = .false.
1777 else
17780 do_more_k = .true.
1779 endif
1780 endif ; enddo ; enddo
17810 if (.not.do_more_k) exit ! All reference interfaces have been found, so stop working downward.
1782 enddo
1783
17840 if (do_more_k) then
17850 if (CS%reset_intxpa_flattest) then
1786 ! There are still points where a correction is needed, so use flattest interface.
17870 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
1788 ! choose top interface first
17890 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
17900 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
17910 p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
17920 p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
17930 inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
17940 dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
17950 delta_z_y(i,J) = abs(e(i,j+1,1)-e(i,j,1))
17960 do k=1,nz
17970 if (abs(e(i,j+1,k+1)-e(i,j,k+1)) < delta_z_y(i,J)) then
1798 ! bottom of layer is less sloped than top. Use this layer
17990 delta_z_y(i,J) = abs(e(i,j+1,k+1)-e(i,j,k+1))
18000 T_int_S(i,J) = T_b(i,j,k) ; T_int_N(i,J) = T_b(i,j+1,k)
18010 S_int_S(i,J) = S_b(i,j,k) ; S_int_N(i,J) = S_b(i,j+1,k)
18020 p_int_S(i,J) = -GxRho0*(e(i,j,k+1) - Z_0p(i,j))
18030 p_int_N(i,J) = -GxRho0*(e(i,j+1,k+1) - Z_0p(i,j))
18040 inty_pa_nonlin(i,J) = inty_pa(i,J,k+1) - 0.5*(pa(i,j,k+1) + pa(i,j+1,k+1))
18050 dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,k+1)-e(i,j,k+1))
1806 endif
1807 enddo
18080 seek_y_cor(i,J) = .false.
1809 endif ; enddo ; enddo
1810 else
1811 ! There are still points where a correction is needed, so use the top interface for lack of a better idea?
18120 do J=Jsq,Jeq ; do i=is,ie ; if (seek_y_cor(i,J)) then
18130 T_int_S(i,J) = T_top(i,j) ; T_int_N(i,J) = T_top(i,j+1)
18140 S_int_S(i,J) = S_top(i,j) ; S_int_N(i,J) = S_top(i,j+1)
18150 p_int_S(i,J) = -GxRho0*(e(i,j,1) - Z_0p(i,j))
18160 p_int_N(i,J) = -GxRho0*(e(i,j+1,1) - Z_0p(i,j))
18170 inty_pa_nonlin(i,J) = inty_pa(i,J,1) - 0.5*(pa(i,j,1) + pa(i,j+1,1))
18180 dgeo_y(i,J) = GV%g_Earth * (e(i,j+1,1)-e(i,j,1))
18190 seek_y_cor(i,J) = .false.
1820 endif ; enddo ; enddo
1821 endif
1822 endif
1823
18240 do J=Jsq,Jeq
18250 do i=is,ie
1826 ! This expression assumes that temperature and salinity vary linearly with hieght
1827 ! between the corners of the reference interfaces found above to get a correction to
1828 ! intx_pa that takes nonlinearities in the equation of state into account.
1829 ! It is derived from a 5 point quadrature estimate of the integral with a large-scale
1830 ! linear correction so that the pressures and heights match at the end-points. It turns
1831 ! out that this linear correction cancels out the mid-point density anomaly.
1832 ! This can be used without masking because dgeo_y and inty_pa_nonlin are 0 over land.
18330 T5(1) = T_Int_S(i,J) ; S5(1) = S_Int_S(i,J) ; p5(1) = p_Int_S(i,J)
18340 T5(5) = T_Int_N(i,J) ; S5(5) = S_Int_N(i,J) ; p5(5) = p_Int_N(i,J)
18350 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))
18360 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))
18370 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))
18380 call calculate_density(T5, S5, p5, r5, tv%eqn_of_state, rho_ref=rho_ref)
1839
1840 ! Note the consistency with the linear form below because (4.75 + 5.5/2) / 90 = 1/12
1841 inty_pa_cor_ri(i,J) = C1_90 * (4.75*(r5(5)-r5(1)) + 5.5*(r5(4)-r5(2))) * dgeo_y(i,J) - &
18420 inty_pa_nonlin(i,J)
1843 enddo
1844 enddo
1845
1846 ! Correct intx_pa and inty_pa at each interface using vertically constant corrections.
18470 do K=1,nz+1 ; do j=js,je ; do I=Isq,Ieq
18480 intx_pa(I,j,K) = intx_pa(I,j,K) + intx_pa_cor_ri(I,j)
1849 enddo ; enddo ; enddo
1850
18510 do K=1,nz+1 ; do J=Jsq,Jeq ; do i=is,ie
18520 inty_pa(i,J,K) = inty_pa(i,J,K) + inty_pa_cor_ri(i,J)
1853 enddo ; enddo ; enddo
1854
1855 ! TODO temporarily move back to CPU
1856 !$omp target update to(intx_pa, inty_pa)
1857 endif ! intx_pa and inty_pa have now been reset to reflect the properties of an unimpeded interface.
1858
1859 ! NOTE: None of these `enter data` statements should be here. They are only
1860 ! temporarily being used to set up the data regions below.
1861 ! Eventually, they should be set up *outside* of the function.
1862
1863 !$omp target enter data if(use_EOS) &
1864 !$omp map(to: tv, EOSdom2d)
1865
1866 ! NOTE: e_sal condition could be sharpened, but this is close enough.
1867 !$omp target enter data map(to: e_tidal_eq, e_tidal_sal, e_sal_and_tide) if (CS%tides)
1868 !$omp target enter data map(to: e_sal) if (CS%calculate_SAL)
1869
1870 !$omp target data map(to: h)
1871
187224 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
1873 PFu(I,j,k) = (((pa(i,j,K)*h(i,j,k) + intz_dpa(i,j,k)) - &
1874 (pa(i+1,j,K)*h(i+1,j,k) + intz_dpa(i+1,j,k))) + &
1875 ((h(i+1,j,k) - h(i,j,k)) * intx_pa(I,j,K) - &
1876 (e(i+1,j,K+1) - e(i,j,K+1)) * intx_dpa(I,j,k) * GV%Z_to_H)) * &
1877 ((2.0*I_Rho0*G%IdxCu(I,j)) / &
187813245168 ((h(i,j,k) + h(i+1,j,k)) + h_neglect))
1879 enddo
1880
1881 ! Compute pressure gradient in y direction
188224 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
1883 PFv(i,J,k) = (((pa(i,j,K)*h(i,j,k) + intz_dpa(i,j,k)) - &
1884 (pa(i,j+1,K)*h(i,j+1,k) + intz_dpa(i,j+1,k))) + &
1885 ((h(i,j+1,k) - h(i,j,k)) * inty_pa(i,J,K) - &
1886 (e(i,j+1,K+1) - e(i,j,K+1)) * inty_dpa(i,J,k) * GV%Z_to_H)) * &
1887 ((2.0*I_Rho0*G%IdyCv(i,J)) / &
188813354584 ((h(i,j,k) + h(i,j+1,k)) + h_neglect))
1889 enddo
1890
1891 !$omp end target data
1892
189324 if (CS%tides_answer_date>20230630 .and. CS%bq_sal_tides) then
1894 ! Calculate SAL geopotential anomaly and add its gradient to pressure
1895 ! gradient force
18960 if (CS%calculate_SAL) then
18970 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
18980 PFu(I,j,k) = PFu(I,j,k) + (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j)
1899 enddo
19000 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
19010 PFv(i,J,k) = PFv(i,J,k) + (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J)
1902 enddo
1903 endif
1904
1905 ! Calculate tidal geopotential anomaly and add its gradient to pressure
1906 ! gradient force
19070 if (CS%tides) then
19080 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
1909 PFu(I,j,k) = PFu(I,j,k) + ((e_tidal_eq(i+1,j) + e_tidal_sal(i+1,j)) &
19100 - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdxCu(I,j)
1911 enddo
19120 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
1913 PFv(i,J,k) = PFv(i,J,k) + ((e_tidal_eq(i,j+1) + e_tidal_sal(i,j+1)) &
19140 - (e_tidal_eq(i,j) + e_tidal_sal(i,j))) * GV%g_Earth * G%IdyCv(i,J)
1915 enddo
1916 endif
1917 endif
1918
191924 if (CS%GFS_scale < 1.0) then
1920 ! Adjust the Montgomery potential to make this a reduced gravity model.
1921 !$omp target data map(alloc: dM)
1922
19230 if (use_EOS) then
1924 !$omp target data map(alloc: rho_in_situ)
1925
19260 if (use_p_atm) then
1927 call calculate_density(tv_tmp%T(:,:,1), tv_tmp%S(:,:,1), p_atm, rho_in_situ, &
19280 tv%eqn_of_state, EOSdom2d)
1929 else
1930 !$omp target data map(alloc: p0)
1931
1932 !$omp target
1933 !$omp parallel loop collapse(2)
19340 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
19350 p0(i,j) = 0.
1936 enddo ; enddo
1937 !$omp end target
1938
1939 call calculate_density(tv_tmp%T(:,:,1), tv_tmp%S(:,:,1), p0, rho_in_situ, &
19400 tv%eqn_of_state, EOSdom2d)
1941 !$omp end target data
1942 endif
1943
1944 !$omp target
1945 !$omp parallel loop collapse(2)
19460 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
19470 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i,j)) * (e(i,j,1) - G%Z_ref)
1948 enddo ; enddo
1949 !$omp end target
1950
1951 !$omp end target data
1952 else
19530 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
19540 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * GV%Rlay(1)) * (e(i,j,1) - G%Z_ref)
1955 enddo
1956 endif
1957
19580 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
19590 PFu(I,j,k) = PFu(I,j,k) - (dM(i+1,j) - dM(i,j)) * G%IdxCu(I,j)
1960 enddo
19610 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
19620 PFv(i,J,k) = PFv(i,J,k) - (dM(i,j+1) - dM(i,j)) * G%IdyCv(i,J)
1963 enddo
1964
1965 !$omp end target data
1966 endif
1967
196824 if (present(pbce)) then
196924 call set_pbce_Bouss(e, tv_tmp, G, GV, US, rho0_set_pbce, CS%GFS_scale, pbce, nkblock=nkblock)
1970 endif
1971
1972 !$omp target exit data if (use_EOS) &
1973 !$omp map(delete:tv_tmp, tv_tmp%T, tv_tmp%S, tv, tv%eqn_of_state, EOSdom2d)
1974
1975 !$omp target exit data map(delete: Z_0p) if (use_EOS)
1976
1977 !$omp target exit data &
1978 !$omp map(delete: pa, dpa) &
1979 !$omp map(delete: intx_pa, inty_pa, intx_dpa, inty_dpa, intz_dpa)
1980
198124 if (present(eta)) then
1982 ! eta is the sea surface height relative to a time-invariant geoid, for comparison with
1983 ! what is used for eta in btstep. See how e was calculated about 200 lines above.
19842952 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
1985184488 eta(i,j) = e(i,j,1)*GV%Z_to_H
1986 enddo
1987
198824 if (CS%tides .and. (.not.CS%bq_sal_tides)) then
19890 if (CS%tides_answer_date>20230630) then
19900 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
19910 eta(i,j) = eta(i,j) + (e_tidal_eq(i,j)+e_tidal_sal(i,j))*GV%Z_to_H
1992 enddo
1993 else
19940 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
19950 eta(i,j) = eta(i,j) + e_sal_and_tide(i,j)*GV%Z_to_H
1996 enddo
1997 endif
1998 endif
1999
200024 if (CS%calculate_SAL .and. (CS%tides_answer_date>20230630) .and. (.not.CS%bq_sal_tides)) then
20010 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
20020 eta(i,j) = eta(i,j) + e_sal(i,j)*GV%Z_to_H
2003 enddo
2004 endif
2005 endif
2006
2007 !$omp target exit data map(delete: e, T_t, T_b, S_t, S_b)
2008 !$omp target exit data map(delete: e_tidal_eq, e_tidal_sal, e_sal_and_tide) if (CS%tides)
2009 !$omp target exit data map(delete: e_sal) if (CS%calculate_SAL)
2010
201124 if (CS%use_stanley_pgf) then
2012 ! Calculated diagnostics related to the Stanley parameterization
20130 zeros(:) = 0.0
20140 EOSdom_h(:) = EOS_domain(G%HI)
20150 if ((CS%id_p_stanley>0) .or. (CS%id_rho_pgf>0) .or. (CS%id_rho_stanley_pgf>0)) then
2016 ! Find the pressure at the mid-point of each layer.
20170 H_to_RL2_T2 = GV%g_Earth*GV%H_to_RZ
20180 if (use_p_atm) then
20190 do j=js,je ; do i=is,ie
20200 p_stanley(i,j,1) = 0.5*h(i,j,1) * H_to_RL2_T2 + p_atm(i,j)
2021 enddo ; enddo
2022 else
20230 do j=js,je ; do i=is,ie
20240 p_stanley(i,j,1) = 0.5*h(i,j,1) * H_to_RL2_T2
2025 enddo ; enddo
2026 endif
20270 do k=2,nz ; do j=js,je ; do i=is,ie
20280 p_stanley(i,j,k) = p_stanley(i,j,k-1) + 0.5*(h(i,j,k-1) + h(i,j,k)) * H_to_RL2_T2
2029 enddo ; enddo ; enddo
2030 endif
20310 if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag)
20320 if (CS%id_rho_pgf>0) then
20330 do k=1,nz ; do j=js,je
2034 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), zeros, &
20350 zeros, zeros, rho_pgf(:,j,k), tv%eqn_of_state, EOSdom_h)
2036 enddo ; enddo
20370 call post_data(CS%id_rho_pgf, rho_pgf, CS%diag)
2038 endif
20390 if (CS%id_rho_stanley_pgf>0) then
20400 do k=1,nz ; do j=js,je
2041 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_stanley(:,j,k), tv%varT(:,j,k), &
20420 zeros, zeros, rho_stanley_pgf(:,j,k), tv%eqn_of_state, EOSdom_h)
2043 enddo ; enddo
20440 call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag)
2045 endif
2046 endif
2047
204824 if (CS%id_MassWt_u>0) call post_data(CS%id_MassWt_u, MassWt_u, CS%diag)
204924 if (CS%id_MassWt_v>0) call post_data(CS%id_MassWt_v, MassWt_v, CS%diag)
2050
205124 if (CS%id_rho_pgf>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag)
205224 if (CS%id_rho_stanley_pgf>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag)
205324 if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag)
2054
2055 ! Diagnostics for tidal forcing and SAL height anomaly
205624 if (CS%id_e_tide>0) then
2057 ! To be consistent with old runs, tidal forcing diagnostic also includes total SAL.
2058 ! New diagnostics are given for each individual field.
20590 if (CS%tides_answer_date>20230630) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
20600 e_sal_and_tide(i,j) = e_sal(i,j) + e_tidal_eq(i,j) + e_tidal_sal(i,j)
2061 enddo ; enddo ; endif
20620 call post_data(CS%id_e_tide, e_sal_and_tide, CS%diag)
2063 endif
206424 if (CS%id_e_sal>0) call post_data(CS%id_e_sal, e_sal, CS%diag)
206524 if (CS%id_e_tidal_eq>0) call post_data(CS%id_e_tidal_eq, e_tidal_eq, CS%diag)
206624 if (CS%id_e_tidal_sal>0) call post_data(CS%id_e_tidal_sal, e_tidal_sal, CS%diag)
2067
2068 ! Diagnostics for tidal forcing and SAL horizontal gradients
206924 if (CS%calculate_SAL .and. ((associated(ADp%sal_u) .or. associated(ADp%sal_v)))) then
20700 if (CS%tides) then ; do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
20710 e_sal(i,j) = e_sal(i,j) + e_tidal_sal(i,j)
2072 enddo ; enddo ; endif
20730 if (CS%bq_sal_tides) then
2074 ! sal_u = ( e(i+1) - e(i) ) * g / dx
20750 if (associated(ADp%sal_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
20760 ADp%sal_u(I,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * GV%g_Earth * G%IdxCu(I,j)
2077 enddo ; enddo ; enddo ; endif
20780 if (associated(ADp%sal_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
20790 ADp%sal_v(i,J,k) = (e_sal(i,j+1) - e_sal(i,j)) * GV%g_Earth * G%IdyCv(i,J)
2080 enddo ; enddo ; enddo ; endif
2081 else
2082 ! sal_u = ( e(i+1) - e(i) ) * g / dx * (rho(k) / rho0)
20830 if (associated(ADp%sal_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
2084 ADp%sal_u(I,j,k) = (e_sal(i+1,j) - e_sal(i,j)) * G%IdxCu(I,j) * I_Rho0 * &
20850 (2.0 * intx_dpa(I,j,k) * GV%Z_to_H / ((h(i,j,k) + h(i+1,j,k)) + h_neglect) + GxRho_ref)
2086 enddo ; enddo ; enddo ; endif
20870 if (associated(ADp%sal_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
2088 ADp%sal_v(i,J,k) = (e_sal(i,j+1) - e_sal(i,j)) * G%IdyCv(i,J) * I_Rho0 * &
20890 (2.0 * inty_dpa(i,J,k) * GV%Z_to_H / ((h(i,j,k) + h(i,j+1,k)) + h_neglect) + GxRho_ref)
2090 enddo ; enddo ; enddo ; endif
2091 endif
20920 if (CS%id_sal_u>0) call post_data(CS%id_sal_u, ADp%sal_u, CS%diag)
20930 if (CS%id_sal_v>0) call post_data(CS%id_sal_v, ADp%sal_v, CS%diag)
2094 endif
2095
209624 if (CS%tides .and. ((associated(ADp%tides_u) .or. associated(ADp%tides_v)))) then
20970 if (CS%bq_sal_tides) then
2098 ! tides_u = ( e(i+1) - e(i) ) * g / dx
20990 if (associated(ADp%tides_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
21000 ADp%tides_u(I,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * GV%g_Earth * G%IdxCu(I,j)
2101 enddo ; enddo ; enddo ; endif
21020 if (associated(ADp%tides_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
21030 ADp%tides_v(i,J,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * GV%g_Earth * G%IdyCv(i,J)
2104 enddo ; enddo ; enddo ; endif
2105 else
2106 ! tides_u = ( e(i+1) - e(i) ) * g / dx * (rho(k) / rho0)
21070 if (associated(ADp%tides_u)) then ; do k=1,nz ; do j=js,je ; do I=Isq,Ieq
2108 ADp%tides_u(I,j,k) = (e_tidal_eq(i+1,j) - e_tidal_eq(i,j)) * G%IdxCu(I,j) * I_Rho0 * &
21090 (2.0 * intx_dpa(I,j,k) * GV%Z_to_H / ((h(i,j,k) + h(i+1,j,k)) + h_neglect) + GxRho_ref)
2110 enddo ; enddo ; enddo ; endif
21110 if (associated(ADp%tides_v)) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
2112 ADp%tides_v(i,J,k) = (e_tidal_eq(i,j+1) - e_tidal_eq(i,j)) * G%IdyCv(i,J) * I_Rho0 * &
21130 (2.0 * inty_dpa(i,J,k) * GV%Z_to_H / ((h(i,j,k) + h(i,j+1,k)) + h_neglect) + GxRho_ref)
2114 enddo ; enddo ; enddo ; endif
2115 endif
21160 if (CS%id_tides_u>0) call post_data(CS%id_tides_u, ADp%tides_u, CS%diag)
21170 if (CS%id_tides_v>0) call post_data(CS%id_tides_v, ADp%tides_v, CS%diag)
2118 endif
211972end subroutine PressureForce_FV_Bouss
2120
2121!> Initializes the finite volume pressure gradient control structure
21221subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL_CSp, tides_CSp)
2123 type(time_type), target, intent(in) :: Time !< Current model time
2124 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2125 type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2126 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2127 type(param_file_type), intent(in) :: param_file !< Parameter file handles
2128 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
2129 type(PressureForce_FV_CS), intent(inout) :: CS !< Finite volume PGF control structure
2130 type(accel_diag_ptrs), pointer :: ADp !< Acceleration diagnostic pointers
2131 type(SAL_CS), intent(in), target, optional :: SAL_CSp !< SAL control structure
2132 type(tidal_forcing_CS), intent(in), target, optional :: tides_CSp !< Tides control structure
2133
2134 ! Local variables
2135 real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
2136 ! temperature variance [nondim]
2137 integer :: default_answer_date ! Global answer date
2138 logical :: use_temperature ! If true, temperature and salinity are used as state variables.
2139 logical :: use_EOS ! If true, density calculated from T & S using an equation of state.
2140 logical :: useMassWghtInterp ! If true, use near-bottom mass weighting for T and S
2141 logical :: MassWghtInterpTop ! If true, use near-surface mass weighting for T and S under ice shelves
2142 logical :: MassWghtInterp_NonBous_bug ! If true, use a buggy mass weighting when non-Boussinesq
2143 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
2144 ! recreate the bugs, or if false bugs are only used if actively selected.
2145 ! This include declares and sets the variable "version".
2146# include "version_variable.h"
2147 character(len=40) :: mdl ! This module's name.
2148 logical :: use_ALE ! If true, use the Vertical Lagrangian Remap algorithm
2149 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
2150
21511 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
21521 IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
2153
21541 CS%initialized = .true.
21551 CS%diag => diag ; CS%Time => Time
21561 if (present(tides_CSp)) &
21571 CS%tides_CSp => tides_CSp
21581 if (present(SAL_CSp)) &
21591 CS%SAL_CSp => SAL_CSp
2160
21611 mdl = "MOM_PressureForce_FV"
21621 call log_version(param_file, mdl, version, "")
2163 call get_param(param_file, mdl, "DEBUG", CS%debug, &
2164 "If true, write out verbose debugging data.", &
21651 default=.false., debuggingParam=.true., do_not_log=.true.)
2166 call get_param(param_file, mdl, "RHO_PGF_REF", CS%rho_ref, &
2167 "The reference density that is subtracted off when calculating pressure "//&
2168 "gradient forces. Its inverse is subtracted off of specific volumes when "//&
2169 "in non-Boussinesq mode. The default is RHO_0.", &
21701 units="kg m-3", default=GV%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R)
2171 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
21721 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
2173 call get_param(param_file, mdl, "RHO_PGF_REF_BUG", CS%rho_ref_bug, &
2174 "If true, recover a bug that RHO_0 (the mean seawater density in Boussinesq mode) "//&
2175 "and RHO_PGF_REF (the subtracted reference density in finite volume pressure "//&
2176 "gradient forces) are incorrectly interchanged in several instances in Boussinesq mode.", &
21771 default=enable_bugs)
2178 call get_param(param_file, mdl, "TIDES", CS%tides, &
21791 "If true, apply tidal momentum forcing.", default=.false.)
21801 if (CS%tides) then
2181 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2182 "This sets the default value for the various _ANSWER_DATE parameters.", &
21830 default=99991231)
2184 call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, "The vintage of "//&
2185 "self-attraction and loading (SAL) and tidal forcing calculations. Setting "//&
2186 "dates before 20230701 recovers old answers (Boussinesq and non-Boussinesq "//&
2187 "modes) when SAL is part of the tidal forcing calculation. The answer "//&
2188 "difference is only at bit level and due to a reordered summation. Setting "//&
2189 "dates before 20250201 recovers answers (Boussinesq mode) that interface "//&
2190 "heights are modified before pressure force integrals are calculated.", &
21910 default=default_answer_date, do_not_log=(.not.CS%tides))
2192 endif
2193 call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
21941 "If true, calculate self-attraction and loading.", default=CS%tides)
21951 if (CS%calculate_SAL) &
2196 call get_param(param_file, '', "SAL_USE_BPA", CS%sal_use_bpa, default=.false., &
21970 do_not_log=.true.)
21981 if ((CS%tides .or. CS%calculate_SAL) .and. GV%Boussinesq) &
2199 call get_param(param_file, mdl, "BOUSSINESQ_SAL_TIDES", CS%bq_sal_tides, "If true, "//&
2200 "in Boussinesq mode, use an alternative method to include self-attraction "//&
2201 "and loading (SAL) and tidal forcings in pressure gradient, in which their "//&
2202 "gradients are calculated separately, instead of adding geopotential "//&
2203 "anomalies as corrections to the interface height. This alternative method "//&
2204 "elimates a baroclinic component of the SAL and tidal forcings.", &
22050 default=.false.)
2206 call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
2207 "If true, Temperature and salinity are used as state variables.", &
22081 default=.true., do_not_log=.true.)
2209 call get_param(param_file, "MOM", "USE_EOS", use_EOS, &
2210 "If true, density is calculated from temperature and "//&
2211 "salinity with an equation of state. If USE_EOS is "//&
2212 "true, ENABLE_THERMODYNAMICS must be true as well.", &
22131 default=use_temperature, do_not_log=.true.)
2214
2215 call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", CS%use_SSH_in_Z0p, &
2216 "If true, include contributions from the sea surface height in the height-based "//&
2217 "pressure used in the equation of state calculations for the Boussinesq pressure "//&
2218 "gradient forces, including adjustments for atmospheric or sea-ice pressure.", &
22191 default=.false., do_not_log=.not.GV%Boussinesq)
22201 if (CS%tides .and. CS%tides_answer_date<=20250131 .and. CS%use_SSH_in_Z0p) &
2221 call MOM_error(FATAL, trim(mdl) // ", PressureForce_FV_init: SSH_IN_EOS_PRESSURE_FOR_PGF "//&
22220 "needs to be FALSE to recover tide answers before 20250131.")
2223
2224 call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
2225 "If True, use the ALE algorithm (regridding/remapping). "//&
22261 "If False, use the layered isopycnal algorithm.", default=.false. )
2227 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT", useMassWghtInterp, &
2228 "If true, use mass weighting when interpolating T/S for integrals "//&
2229 "near the bathymetry in FV pressure gradient calculations.", &
22301 default=.false.)
2231 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP", MassWghtInterpTop, &
2232 "If true and MASS_WEIGHT_IN_PRESSURE_GRADIENT is true, use mass weighting when "//&
2233 "interpolating T/S for integrals near the top of the water column in FV "//&
2234 "pressure gradient calculations. ", &
22351 default=useMassWghtInterp)
2236 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_NONBOUS_BUG", MassWghtInterp_NonBous_bug, &
2237 "If true, use a masking bug in non-Boussinesq calculations with mass weighting "//&
2238 "when interpolating T/S for integrals near the bathymetry in FV pressure "//&
2239 "gradient calculations.", &
22401 default=.false., do_not_log=(GV%Boussinesq .or. (.not.useMassWghtInterp)))
2241 call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_VANISHED_ONLY", CS%MassWghtInterpVanOnly, &
2242 "If true, use mass weighting when interpolating T/S for integrals "//&
2243 "only if one side is vanished according to RESET_INTXPA_H_NONVANISHED. ", &
22441 default=.false.)
2245
22461 CS%MassWghtInterp = 0
22471 if (useMassWghtInterp) &
22481 CS%MassWghtInterp = ibset(CS%MassWghtInterp, 0) ! Same as CS%MassWghtInterp + 1
22491 if (MassWghtInterpTop) &
22501 CS%MassWghtInterp = ibset(CS%MassWghtInterp, 1) ! Same as CS%MassWghtInterp + 2
22511 if ((.not.GV%Boussinesq) .and. MassWghtInterp_NonBous_bug) &
22520 CS%MassWghtInterp = ibset(CS%MassWghtInterp, 3) ! Same as CS%MassWghtInterp + 8
2253
2254 call get_param(param_file, mdl, "CORRECTION_INTXPA", CS%correction_intxpa, &
2255 "If true, use a correction for surface pressure curvature in intx_pa.", &
22561 default=.false., do_not_log=.not.use_EOS)
2257 call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL", CS%reset_intxpa_integral, &
2258 "If true, reset INTXPA to match pressures at first nonvanished cell. "//&
22591 "Includes pressure correction.", default=.false., do_not_log=.not.use_EOS)
2260 call get_param(param_file, mdl, "RESET_INTXPA_INTEGRAL_FLATTEST", CS%reset_intxpa_flattest, &
2261 "If true, use flattest interface as reference interface where there is no "//&
2262 "better choice for RESET_INTXPA_INTEGRAL. Otherwise, use surface interface.", &
22631 default=.false., do_not_log=.not.use_EOS)
22641 if (.not.use_EOS) then ! These options do nothing without an equation of state.
22650 CS%correction_intxpa = .false.
22660 CS%reset_intxpa_integral = .false.
22670 CS%reset_intxpa_flattest = .false.
2268 endif
2269 call get_param(param_file, mdl, "RESET_INTXPA_H_NONVANISHED", CS%h_nonvanished, &
2270 "A minimal layer thickness that indicates that a layer is thick enough to usefully "//&
2271 "reestimate the pressure integral across the interface below.", &
22721 default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=.not.CS%reset_intxpa_integral)
2273 call get_param(param_file, mdl, "USE_INACCURATE_PGF_RHO_ANOM", CS%use_inaccurate_pgf_rho_anom, &
2274 "If true, use a form of the PGF that uses the reference density "//&
22751 "in an inaccurate way. This is not recommended.", default=.false.)
2276 call get_param(param_file, mdl, "RECONSTRUCT_FOR_PRESSURE", CS%reconstruct, &
2277 "If True, use vertical reconstruction of T & S within "//&
2278 "the integrals of the FV pressure gradient calculation. "//&
2279 "If False, use the constant-by-layer algorithm. "//&
2280 "The default is set by USE_REGRIDDING.", &
22811 default=use_ALE )
2282 call get_param(param_file, mdl, "PRESSURE_RECONSTRUCTION_SCHEME", CS%Recon_Scheme, &
2283 "Order of vertical reconstruction of T/S to use in the "//&
2284 "integrals within the FV pressure gradient calculation.\n"//&
2285 " 0: PCM or no reconstruction.\n"//&
2286 " 1: PLM reconstruction.\n"//&
2287 " 2: PPM reconstruction.\n"//&
22881 " 3: PLM with least squares slope.", default=1)
2289 call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION_PRESSURE", CS%boundary_extrap, &
2290 "If true, the reconstruction of T & S for pressure in "//&
2291 "boundary cells is extrapolated, rather than using PCM "//&
2292 "in these cells. If true, the same order polynomial is "//&
22931 "used as is used for the interior cells.", default=.true.)
2294 call get_param(param_file, mdl, "USE_STANLEY_PGF", CS%use_stanley_pgf, &
2295 "If true, turn on Stanley SGS T variance parameterization "// &
22961 "in PGF code.", default=.false.)
2297 call get_param(param_file, mdl, "PGF_PLM_NKBLOCK", CS%nkblock, &
2298 "Vertical block size for the PLM pressure gradient density integral. "//&
2299 "0 processes all layers in a single block.", &
23001 default=default_nkblock)
2301 call get_param(param_file, mdl, "PGF_PLM_NIBLOCK", CS%niblock_plm, &
2302 "i-tile size for the PLM pressure gradient density integral. "//&
2303 "0 uses the full i compute domain.", &
23041 default=default_niblock_plm)
2305 call get_param(param_file, mdl, "PGF_PLM_NJBLOCK", CS%njblock_plm, &
2306 "j-tile size for the PLM pressure gradient density integral. "//&
2307 "0 uses the full j compute domain.", &
23081 default=default_njblock_plm)
23091 if (CS%use_stanley_pgf) then
2310 call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
2311 "Coefficient correlating the temperature gradient and SGS T variance.", &
23120 units="nondim", default=-1.0, do_not_log=.true.)
23130 if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
23140 "STANLEY_COEFF must be set >= 0 if USE_STANLEY_PGF is true.")
2315
2316 CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, &
23170 Time, 'rho in PGF', 'kg m-3', conversion=US%R_to_kg_m3)
2318 CS%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, &
23190 Time, 'rho in PGF with Stanley correction', 'kg m-3', conversion=US%R_to_kg_m3)
2320 CS%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, &
23210 Time, 'p in PGF with Stanley correction', 'Pa', conversion=US%RL2_T2_to_Pa)
2322 endif
23231 if (CS%calculate_SAL) then
2324 CS%id_e_sal = register_diag_field('ocean_model', 'e_sal', diag%axesT1, Time, &
23250 'Self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m)
2326 CS%id_sal_u = register_diag_field('ocean_model', 'SAL_u', diag%axesCuL, Time, &
23270 'Zonal Acceleration due to self-attraction and loading', 'm s-2', conversion=US%L_T2_to_m_s2)
2328 CS%id_sal_v = register_diag_field('ocean_model', 'SAL_v', diag%axesCvL, Time, &
23290 'Meridional Acceleration due to self-attraction and loading', 'm s-2', conversion=US%L_T2_to_m_s2)
23300 if (CS%id_sal_u > 0) &
23310 call safe_alloc_ptr(ADp%sal_u, IsdB, IedB, jsd, jed, nz)
23320 if (CS%id_sal_v > 0) &
23330 call safe_alloc_ptr(ADp%sal_v, isd, ied, JsdB, JedB, nz)
2334 endif
23351 if (CS%tides) then
2336 CS%id_e_tide = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, Time, &
23370 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=US%Z_to_m)
2338 CS%id_e_tidal_eq = register_diag_field('ocean_model', 'e_tide_eq', diag%axesT1, Time, &
23390 'Equilibrium tides height anomaly', 'meter', conversion=US%Z_to_m)
2340 CS%id_e_tidal_sal = register_diag_field('ocean_model', 'e_tide_sal', diag%axesT1, Time, &
23410 'Read-in tidal self-attraction and loading height anomaly', 'meter', conversion=US%Z_to_m)
2342 CS%id_tides_u = register_diag_field('ocean_model', 'tides_u', diag%axesCuL, Time, &
23430 'Zonal Acceleration due to tidal forcing', 'm s-2', conversion=US%L_T2_to_m_s2)
2344 CS%id_tides_v = register_diag_field('ocean_model', 'tides_v', diag%axesCvL, Time, &
23450 'Meridional Acceleration due to tidal forcing', 'm s-2', conversion=US%L_T2_to_m_s2)
23460 if (CS%id_tides_u > 0) &
23470 call safe_alloc_ptr(ADp%tides_u, IsdB, IedB, jsd, jed, nz)
23480 if (CS%id_tides_v > 0) &
23490 call safe_alloc_ptr(ADp%tides_v, isd, ied, JsdB, JedB, nz)
2350 endif
2351
2352 CS%id_MassWt_u = register_diag_field('ocean_model', 'MassWt_u', diag%axesCuL, Time, &
23531 'The fractional mass weighting at u-point PGF calculations', 'nondim')
2354 CS%id_MassWt_v = register_diag_field('ocean_model', 'MassWt_v', diag%axesCvL, Time, &
23551 'The fractional mass weighting at v-point PGF calculations', 'nondim')
2356
23571 CS%GFS_scale = 1.0
23581 if (GV%g_prime(1) /= GV%g_Earth) CS%GFS_scale = GV%g_prime(1) / GV%g_Earth
2359
2360 !$omp target update to (CS)
2361
23621 call log_param(param_file, mdl, "GFS / G_EARTH", CS%GFS_scale, units="nondim")
2363
23641end subroutine PressureForce_FV_init
2365
2366!> \namespace mom_pressureforce_fv
2367!!
2368!! Provides the Boussinesq and non-Boussinesq forms of horizontal accelerations
2369!! due to pressure gradients using a vertically integrated finite volume form,
2370!! as described by Adcroft et al., 2008. Integration in the vertical is made
2371!! either by quadrature or analytically.
2372!!
2373!! This form eliminates the thermobaric instabilities that had been a problem with
2374!! previous forms of the pressure gradient force calculation, as described by
2375!! Hallberg, 2005.
2376!!
2377!! Adcroft, A., R. Hallberg, and M. Harrison, 2008: A finite volume discretization
2378!! of the pressure gradient force using analytic integration. Ocean Modelling, 22,
2379!! 106-113. http://doi.org/10.1016/j.ocemod.2008.02.001
2380!!
2381!! Hallberg, 2005: A thermobaric instability of Lagrangian vertical coordinate
2382!! ocean models. Ocean Modelling, 8, 279-300.
2383!! http://dx.doi.org/10.1016/j.ocemod.2004.01.001
2384
23850end module MOM_PressureForce_FV