ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_TerrainPoisson_3D_K.H File Reference
#include <AMReX_FArrayBox.H>
Include dependency graph for ERF_TerrainPoisson_3D_K.H:

Go to the source code of this file.

Functions

template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_x (int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dxinv) noexcept
 
template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_y (int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dyinv) noexcept
 
template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_z (int i, int j, int k, amrex::Array4< T const > const &sol, amrex::Array4< T const > const &zp, T dxinv, T dyinv) noexcept
 
template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void terrpoisson_adotx (int i, int j, int k, amrex::Array4< T > const &y, amrex::Array4< T const > const &x, amrex::Array4< T const > const &zp, T dxinv, T dyinv) noexcept
 

Function Documentation

◆ terrpoisson_adotx()

template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void terrpoisson_adotx ( int  i,
int  j,
int  k,
amrex::Array4< T > const &  y,
amrex::Array4< T const > const &  x,
amrex::Array4< T const > const &  zp,
T  dxinv,
T  dyinv 
)
noexcept
161 {
162  using amrex::Real;
163  Real h_xi, h_eta, h_zeta;
164 
165  // *********************************************************
166  // Hi x-face
167  // *********************************************************
168  // On x-face
169  Real px_hi = (x(i+1,j,k) - x(i,j,k)) * dxinv;
170 
171  // On y-edges
172  Real pz_hi_md_hi = Real(0.5) * ( x(i+1,j ,k+1) + x(i ,j ,k+1)
173  -x(i+1,j ,k ) - x(i ,j ,k ) );
174  h_xi = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i-1,j ,k+1)
175  +zp(i+1,j+1,k+1) - zp(i-1,j+1,k+1) ) * dxinv;
176  h_zeta = Real(0.25) * ( zp(i+1,j ,k+2) - zp(i+1,j ,k )
177  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
178  pz_hi_md_hi *= h_xi / h_zeta;
179 
180  // On y-edges
181  Real pz_hi_md_lo = Real(0.5) * ( x(i+1,j ,k ) + x(i ,j ,k )
182  -x(i+1,j ,k-1) - x(i ,j ,k-1) );
183  h_xi = Real(0.25) * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
184  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
185  h_zeta = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
186  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
187  pz_hi_md_lo *= h_xi / h_zeta;
188 
189  // On x-face
190  px_hi -= Real(0.5) * ( pz_hi_md_hi + pz_hi_md_lo );
191 
192  // *********************************************************
193  // Lo x-face
194  // ********************************************************* // On x-face
195  Real px_lo = (x(i,j,k) - x(i-1,j,k)) * dxinv;
196 
197  // On y-edges
198  Real pz_lo_md_hi = Real(0.5) * ( x(i,j,k+1) + x(i-1,j,k+1)
199  -x(i,j,k ) - x(i-1,j,k ) );
200  h_xi = Real(0.25) * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
201  +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
202  h_zeta = Real(0.25) * ( zp(i,j ,k+2) - zp(i ,j ,k )
203  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
204  pz_lo_md_hi *= h_xi / h_zeta;
205 
206  // On y-edges
207  Real pz_lo_md_lo = Real(0.5) * ( x(i,j,k ) + x(i-1,j,k )
208  -x(i,j,k-1) - x(i-1,j,k-1) );
209  h_xi = Real(0.25) * ( zp(i,j ,k ) - zp(i-2,j ,k )
210  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
211  h_zeta = Real(0.25) * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
212  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
213  pz_lo_md_lo *= h_xi / h_zeta;
214 
215  // On x-face
216  px_lo -= Real(0.5) * ( pz_lo_md_hi + pz_lo_md_lo );
217 
218  // *********************************************************
219  // Hi y-face
220  // *********************************************************
221  // On y-face
222  Real py_hi = (x(i,j+1,k) - x(i,j,k)) * dyinv;
223 
224  // On x-edges
225  Real pz_md_hi_hi = Real(0.5) * ( x(i,j+1,k+1) + x(i,j,k+1)
226  -x(i,j+1,k ) - x(i,j,k ) );
227  h_eta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j-1,k+1)
228  +zp(i+1,j+1,k+1) - zp(i+1,j-1,k+1) ) * dyinv;
229  h_zeta = Real(0.25) * ( zp(i ,j+1,k+2) - zp(i ,j+1,k )
230  +zp(i+1,j+1,k+2) - zp(i+1,j+1,k ) );
231  pz_md_hi_hi *= h_eta/ h_zeta;
232 
233  // On x-edges
234  Real pz_md_hi_lo = Real(0.5) * ( x(i,j+1,k ) + x(i,j,k )
235  -x(i,j+1,k-1) - x(i,j,k-1) );
236  h_eta = Real(0.25) * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
237  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
238  h_zeta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
239  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
240  pz_md_hi_lo *= h_eta/ h_zeta;
241 
242  // On y-face
243  py_hi -= Real(0.5) * ( pz_md_hi_hi + pz_md_hi_lo );
244 
245  // *********************************************************
246  // Lo y-face
247  // *********************************************************
248  // On y-face
249  Real py_lo = (x(i,j,k) - x(i,j-1,k)) * dyinv;
250 
251  // On x-edges
252  Real pz_md_lo_hi = Real(0.5) * ( x(i ,j,k+1) + x(i ,j-1,k+1)
253  -x(i ,j,k ) - x(i ,j-1,k ) );
254  h_eta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
255  +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
256  h_zeta = Real(0.25) * ( zp(i ,j,k+2) - zp(i ,j ,k )
257  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
258  pz_md_lo_hi *= h_eta/ h_zeta;
259 
260  // On x-edges
261  Real pz_md_lo_lo = Real(0.5) * ( x(i ,j,k ) + x(i ,j-1,k )
262  -x(i ,j,k-1) - x(i ,j-1,k-1) );
263  h_eta = Real(0.25) * ( zp(i ,j,k ) - zp(i ,j-2,k )
264  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
265  h_zeta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
266  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
267  pz_md_lo_lo *= h_eta/ h_zeta;
268 
269  // On y-face
270  py_lo -= Real(0.5) * ( pz_md_lo_hi + pz_md_lo_lo );
271 
272  // *********************************************************
273  // Hi z-face
274  // *********************************************************
275  // On z-face
276  Real pz_hi = x(i,j,k+1) - x(i,j,k );
277  Real hzeta_inv_on_zhi = 8.0 / ( (zp(i,j,k+2) + zp(i+1,j,k+2) + zp(i,j+1,k+2) + zp(i+1,j+1,k+2))
278  -(zp(i,j,k ) + zp(i+1,j,k ) + zp(i,j+1,k ) + zp(i+1,j+1,k )) );
279  pz_hi *= hzeta_inv_on_zhi;
280 
281  // On corners
282  Real px_hi_md_hi = Real(0.5) * ( x(i+1,j,k+1) - x(i ,j ,k+1)
283  +x(i+1,j,k ) - x(i ,j ,k )) * dxinv;
284  Real px_lo_md_hi = Real(0.5) * ( x(i ,j,k+1) - x(i-1,j ,k+1)
285  +x(i ,j,k ) - x(i-1,j ,k )) * dxinv;
286  Real py_md_hi_hi = Real(0.5) * ( x(i,j+1,k+1) - x(i ,j ,k+1)
287  +x(i,j+1,k ) - x(i ,j ,k )) * dyinv;
288  Real py_md_lo_hi = Real(0.5) * ( x(i,j ,k+1) - x(i ,j-1,k+1)
289  +x(i,j ,k ) - x(i ,j-1,k )) * dyinv;
290 
291  // On z-face
292  Real h_xi_on_zhi = 0.5 * ( zp(i+1,j+1,k+1) + zp(i+1,j,k+1) - zp(i,j+1,k+1) - zp(i,j,k+1) ) * dxinv;
293  Real h_eta_on_zhi = 0.5 * ( zp(i+1,j+1,k+1) + zp(i,j+1,k+1) - zp(i+1,j,k+1) - zp(i,j,k+1) ) * dyinv;
294  //
295  // Note we do not need to recalculate pz_...hi here
296  //
297  pz_hi -= 0.5 * h_xi_on_zhi * ( (px_hi_md_hi + px_lo_md_hi) - (pz_hi_md_hi + pz_lo_md_hi) );
298  pz_hi -= 0.5 * h_eta_on_zhi * ( (py_md_hi_hi + py_md_lo_hi) - (pz_md_hi_hi + pz_md_lo_hi) );
299 
300  // *********************************************************
301  // Lo z-face
302  // *********************************************************
303  // On z-face
304  Real pz_lo = x(i,j,k ) - x(i,j,k-1);
305  Real hzeta_inv_on_zlo = 8.0 / ( (zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1))
306  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
307  pz_lo *= hzeta_inv_on_zlo;
308 
309  // On corners
310  Real px_hi_md_lo = Real(0.5) * ( x(i+1,j ,k ) - x(i ,j ,k )
311  +x(i+1,j ,k-1) - x(i ,j ,k-1)) * dxinv;
312  Real px_lo_md_lo = Real(0.5) * ( x(i ,j ,k ) - x(i-1,j ,k )
313  +x(i ,j ,k-1) - x(i-1,j ,k-1)) * dxinv;
314  Real py_md_hi_lo = Real(0.5) * ( x(i ,j+1,k ) - x(i ,j ,k )
315  +x(i ,j+1,k-1) - x(i ,j ,k-1)) * dyinv;
316  Real py_md_lo_lo = Real(0.5) * ( x(i ,j ,k ) - x(i ,j-1,k )
317  +x(i ,j ,k-1) - x(i ,j-1,k-1)) * dyinv;
318 
319  // On z-face
320  Real h_xi_on_zlo = 0.5 * (zp(i+1,j+1,k) + zp(i+1,j,k) - zp(i,j+1,k) - zp(i,j,k)) * dxinv;
321  Real h_eta_on_zlo = 0.5 * (zp(i+1,j+1,k) + zp(i,j+1,k) - zp(i+1,j,k) - zp(i,j,k)) * dyinv;
322  //
323  // Note we do not need to recalculate pz_...lo here
324  //
325  pz_lo -= 0.5 * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
326  pz_lo -= 0.5 * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
327 
328  // *********************************************************
329  // Version which calls flux routines
330  // *********************************************************
331  //
332  // This option uses calls to the flux routines so there is
333  // some duplicated computation
334  // This option should give the same answer as above
335  //
336  // Real px_lo = -terrpoisson_flux_x(i ,j,k,x,zp,dxinv);
337  // Real px_hi = -terrpoisson_flux_x(i+1,j,k,x,zp,dxinv);
338  // Real py_lo = -terrpoisson_flux_y(i,j ,k,x,zp,dxinv);
339  // Real py_hi = -terrpoisson_flux_y(i,j+1,k,x,zp,dyinv);
340  // Real pz_lo = -terrpoisson_flux_z(i,j,k ,x,zp,dxinv,dyinv);
341  // Real pz_hi = -terrpoisson_flux_z(i,j,k+1,x,zp,dxinv,dyinv);
342  //
343  // *********************************************************
344  // Adotx
345  // *********************************************************
346  Real invdJ = 4.0 / ( zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1)
347  -zp(i,j,k ) - zp(i+1,j,k ) - zp(i,j+1,k ) - zp(i+1,j+1,k ) );
348 
349  Real ax_lo = .5 * (zp(i ,j,k+1) + zp(i ,j+1,k+1) - zp(i ,j,k) - zp(i ,j+1,k));
350  Real ax_hi = .5 * (zp(i+1,j,k+1) + zp(i+1,j+1,k+1) - zp(i+1,j,k) - zp(i+1,j+1,k));
351  Real ay_lo = .5 * (zp(i,j ,k+1) + zp(i+1,j ,k+1) - zp(i,j ,k) - zp(i+1,j ,k));
352  Real ay_hi = .5 * (zp(i,j+1,k+1) + zp(i+1,j+1,k+1) - zp(i,j+1,k) - zp(i+1,j+1,k));
353 
354  y(i,j,k) = ( (ax_hi*px_hi - ax_lo*px_lo) * dxinv
355  +(ay_hi*py_hi - ay_lo*py_lo) * dyinv
356  +( pz_hi - pz_lo) ) * invdJ;
357 }
Coord::y
@ y
Coord::x
@ x

◆ terrpoisson_flux_x()

template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_x ( int  i,
int  j,
int  k,
amrex::Array4< T const > const &  sol,
amrex::Array4< T const > const &  zp,
T  dxinv 
)
noexcept
12 {
13  using amrex::Real;
14 
15  Real h_xi, h_zeta;
16 
17  // On x-face
18  Real px_lo = (sol(i,j,k) - sol(i-1,j,k)) * dxinv;
19 
20  // On y-edges
21  Real pz_lo_md_hi = Real(0.5) * ( sol(i,j,k+1) + sol(i-1,j,k+1)
22  -sol(i,j,k ) - sol(i-1,j,k ) );
23  h_xi = Real(0.25) * ( zp(i,j ,k+1) - zp(i-2,j ,k+1)
24  +zp(i,j+1,k+1) - zp(i-2,j+1,k+1) ) * dxinv;
25  h_zeta = Real(0.25) * ( zp(i,j ,k+2) - zp(i ,j ,k )
26  +zp(i,j+1,k+2) - zp(i ,j+1,k ) );
27  pz_lo_md_hi *= h_xi / h_zeta;
28 
29  // On y-edges
30  Real pz_lo_md_lo = Real(0.5) * ( sol(i,j,k ) + sol(i-1,j,k )
31  -sol(i,j,k-1) - sol(i-1,j,k-1) );
32 
33  h_xi = Real(0.25) * ( zp(i,j ,k ) - zp(i-2,j ,k )
34  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
35  h_zeta = Real(0.25) * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
36  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
37  pz_lo_md_lo *= h_xi / h_zeta;
38 
39  // On x-face
40  px_lo -= Real(0.5) * ( pz_lo_md_hi + pz_lo_md_lo );
41 
42  return -px_lo;
43 }

◆ terrpoisson_flux_y()

template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_y ( int  i,
int  j,
int  k,
amrex::Array4< T const > const &  sol,
amrex::Array4< T const > const &  zp,
T  dyinv 
)
noexcept
50 {
51  using amrex::Real;
52 
53  Real h_eta, h_zeta;
54 
55  Real py_lo = (sol(i,j,k) - sol(i,j-1,k)) * dyinv;
56 
57  // On x-edges
58  Real pz_md_lo_hi = Real(0.5) * ( sol(i,j,k+1) + sol(i,j-1,k+1)
59  -sol(i,j,k ) - sol(i,j-1,k ) );
60  h_eta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j-2,k+1)
61  +zp(i+1,j,k+1) - zp(i+1,j-2,k+1) ) * dyinv;
62  h_zeta = Real(0.25) * ( zp(i ,j,k+2) - zp(i ,j ,k )
63  +zp(i+1,j,k+2) - zp(i+1,j ,k ) );
64  pz_md_lo_hi *= h_eta/ h_zeta;
65 
66  // On x-edges
67  Real pz_md_lo_lo = Real(0.5) * ( sol(i,j,k ) + sol(i,j-1,k )
68  -sol(i,j,k-1) - sol(i,j-1,k-1) );
69 
70  h_eta = Real(0.25) * ( zp(i ,j,k ) - zp(i ,j-2,k )
71  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
72  h_zeta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
73  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
74  pz_md_lo_lo *= h_eta/ h_zeta;
75 
76  // On y-face
77  py_lo -= Real(0.5) * ( pz_md_lo_hi + pz_md_lo_lo );
78 
79  return -py_lo;
80 }

◆ terrpoisson_flux_z()

template<typename T >
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T terrpoisson_flux_z ( int  i,
int  j,
int  k,
amrex::Array4< T const > const &  sol,
amrex::Array4< T const > const &  zp,
T  dxinv,
T  dyinv 
)
noexcept
88 {
89  using amrex::Real;
90 
91  Real h_xi, h_eta, h_zeta;
92 
93  // On z-face
94  Real pz_lo = (sol(i,j,k ) - sol(i,j,k-1));
95  Real hzeta_inv_on_zlo = 8.0 / ( (zp(i,j,k+1) + zp(i+1,j,k+1) + zp(i,j+1,k+1) + zp(i+1,j+1,k+1))
96  -(zp(i,j,k-1) + zp(i+1,j,k-1) + zp(i,j+1,k-1) + zp(i+1,j+1,k-1)) );
97  pz_lo *= hzeta_inv_on_zlo;
98 
99  // On corners
100  Real px_hi_md_lo = Real(0.5) * ( sol(i+1,j ,k ) - sol(i ,j ,k )
101  +sol(i+1,j ,k-1) - sol(i ,j ,k-1)) * dxinv;
102  Real px_lo_md_lo = Real(0.5) * ( sol(i ,j ,k ) - sol(i-1,j ,k )
103  +sol(i ,j ,k-1) - sol(i-1,j ,k-1)) * dxinv;
104  Real py_md_hi_lo = Real(0.5) * ( sol(i ,j+1,k ) - sol(i ,j ,k )
105  +sol(i ,j+1,k-1) - sol(i ,j ,k-1)) * dyinv;
106  Real py_md_lo_lo = Real(0.5) * ( sol(i ,j ,k ) - sol(i ,j-1,k )
107  +sol(i ,j ,k-1) - sol(i ,j-1,k-1)) * dyinv;
108 
109  // On y-edges
110  Real pz_hi_md_lo = Real(0.5) * ( sol(i+1,j,k ) + sol(i,j,k )
111  -sol(i+1,j,k-1) - sol(i,j,k-1) );
112  h_xi = Real(0.25) * ( zp(i+1,j ,k ) - zp(i-1,j ,k )
113  +zp(i+1,j+1,k ) - zp(i-1,j+1,k ) ) * dxinv;
114  h_zeta = Real(0.25) * ( zp(i+1,j ,k+1) - zp(i+1,j ,k-1)
115  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
116  pz_hi_md_lo *= h_xi / h_zeta;
117 
118  // On y-edges
119  Real pz_lo_md_lo = Real(0.5) * ( sol(i,j,k ) + sol(i-1,j,k )
120  -sol(i,j,k-1) - sol(i-1,j,k-1) );
121  h_xi = Real(0.25) * ( zp(i,j ,k ) - zp(i-2,j ,k )
122  +zp(i,j+1,k ) - zp(i-2,j+1,k ) ) * dxinv;
123  h_zeta = Real(0.25) * ( zp(i,j ,k+1) - zp(i ,j ,k-1)
124  +zp(i,j+1,k+1) - zp(i ,j+1,k-1) );
125  pz_lo_md_lo *= h_xi / h_zeta;
126 
127  // On x-edges
128  Real pz_md_hi_lo = Real(0.5) * ( sol(i,j+1,k ) + sol(i,j,k )
129  -sol(i,j+1,k-1) - sol(i,j,k-1) );
130  h_eta = Real(0.25) * ( zp(i ,j+1,k ) - zp(i ,j-1,k)
131  +zp(i+1,j+1,k ) - zp(i+1,j-1,k) ) * dyinv;
132  h_zeta = Real(0.25) * ( zp(i ,j+1,k+1) - zp(i ,j+1,k-1)
133  +zp(i+1,j+1,k+1) - zp(i+1,j+1,k-1) );
134  pz_md_hi_lo *= h_eta/ h_zeta;
135 
136  // On x-edges
137  Real pz_md_lo_lo = Real(0.5) * ( sol(i,j,k ) + sol(i,j-1,k )
138  -sol(i,j,k-1) - sol(i,j-1,k-1) );
139  h_eta = Real(0.25) * ( zp(i ,j,k ) - zp(i ,j-2,k )
140  +zp(i+1,j,k ) - zp(i+1,j-2,k ) ) * dyinv;
141  h_zeta = Real(0.25) * ( zp(i ,j,k+1) - zp(i ,j ,k-1)
142  +zp(i+1,j,k+1) - zp(i+1,j ,k-1) );
143  pz_md_lo_lo *= h_eta/ h_zeta;
144 
145  // On z-face
146  Real h_xi_on_zlo = 0.5 * (zp(i+1,j+1,k) + zp(i+1,j,k) - zp(i,j+1,k) - zp(i,j,k)) * dxinv;
147  Real h_eta_on_zlo = 0.5 * (zp(i+1,j+1,k) + zp(i,j+1,k) - zp(i+1,j,k) - zp(i,j,k)) * dyinv;
148 
149  pz_lo -= 0.5 * h_xi_on_zlo * ( (px_hi_md_lo + px_lo_md_lo) - (pz_hi_md_lo + pz_lo_md_lo) );
150  pz_lo -= 0.5 * h_eta_on_zlo * ( (py_md_hi_lo + py_md_lo_lo) - (pz_md_hi_lo + pz_md_lo_lo) );
151 
152  return -pz_lo;
153 }