ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_buoyancy_utils.H
Go to the documentation of this file.
1 #ifndef ERF_BUOYANCY_UTILS_H_
2 #define ERF_BUOYANCY_UTILS_H_
3 
4 #include <ERF_EOS.H>
5 #include <ERF_Constants.H>
6 
7 AMREX_GPU_DEVICE
8 AMREX_FORCE_INLINE
9 amrex::Real
11  int& j,
12  int& k,
13  amrex::Real const& grav_gpu,
14  amrex::Real const& rd_over_cp,
15  const amrex::Array4<const amrex::Real>& p0_arr,
16  const amrex::Array4<const amrex::Real>& r0_arr,
17  const amrex::Array4<const amrex::Real>& cell_data)
18 {
19  amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
20  amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
21  amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
22  amrex::Real qplus = (t_hi-t0_hi)/t0_hi;
23 
24  amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
25  amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
26  amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
27  amrex::Real qminus = (t_lo-t0_lo)/t0_lo;
28 
29  amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
30  return (-r0_q_avg * grav_gpu);
31 }
32 
33 
34 AMREX_GPU_DEVICE
35 AMREX_FORCE_INLINE
36 amrex::Real
38  int& j,
39  int& k,
40  amrex::Real const& grav_gpu,
41  amrex::Real const& rv_over_rd,
42  const amrex::Array4<const amrex::Real>& p0_arr,
43  const amrex::Array4<const amrex::Real>& r0_arr,
44  const amrex::Array4<const amrex::Real>& cell_data)
45 {
46  amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
47  amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp)/r0_arr(i,j,k-1);
48  amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp)/r0_arr(i,j,k-1);
49  amrex::Real qt_lo = qv_lo + qc_lo;
50  amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);
51 
52  amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k);
53  amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp)/r0_arr(i,j,k);
54  amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp)/r0_arr(i,j,k);
55  amrex::Real qt_hi = qv_hi + qc_hi;
56  amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);
57 
58  amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi);
59 
60  amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1);
61  amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k )) / r0_arr(i,j,k );
62  amrex::Real theta_v_0_lo = theta_d_0_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);
63  amrex::Real theta_v_0_hi = theta_d_0_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);
64 
65  amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi);
66 
67  return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface);
68 }
69 
70 AMREX_GPU_DEVICE
71 AMREX_FORCE_INLINE
72 amrex::Real
74  int& j,
75  int& k,
76  amrex::Real const& grav_gpu,
77  amrex::Real const& rd_over_cp,
78  const amrex::Array4<const amrex::Real>& p0_arr,
79  const amrex::Array4<const amrex::Real>& r0_arr,
80  const amrex::Array4<const amrex::Real>& cell_data)
81 {
82  amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
83  amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
84  amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp));
85  amrex::Real qplus = (t_hi-t0_hi)/t0_hi;
86 
87  amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
88  amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
89  amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
90  amrex::Real qminus = (t_lo-t0_lo)/t0_lo;
91 
92  amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
93  return (-r0_q_avg * grav_gpu);
94 }
95 
96 AMREX_GPU_DEVICE
97 AMREX_FORCE_INLINE
98 amrex::Real
100  int& j,
101  int& k,
102  const int& n_qstate,
103  amrex::Real const& grav_gpu,
104  const amrex::Array4<const amrex::Real>& r0_arr,
105  const amrex::Array4<const amrex::Real>& cell_data)
106 {
107  amrex::Real rhop_hi = cell_data(i,j,k ,Rho_comp) - r0_arr(i,j,k );
108  amrex::Real rhop_lo = cell_data(i,j,k-1,Rho_comp) - r0_arr(i,j,k-1);
109  for (int q_offset(0); q_offset<n_qstate; ++q_offset) {
110  rhop_hi += cell_data(i,j,k ,RhoQ1_comp+q_offset);
111  rhop_lo += cell_data(i,j,k-1,RhoQ1_comp+q_offset);
112  }
113  return( grav_gpu * amrex::Real(0.5) * ( rhop_hi + rhop_lo ) );
114 }
115 
116 
117 AMREX_GPU_DEVICE
118 AMREX_FORCE_INLINE
119 amrex::Real
121  int& j,
122  int& k,
123  const int& n_qstate,
124  amrex::Real const& grav_gpu,
125  amrex::Real* rho_d_ptr,
126  amrex::Real* theta_d_ptr,
127  amrex::Real* qv_d_ptr,
128  amrex::Real* qc_d_ptr,
129  amrex::Real* qp_d_ptr,
130  const amrex::Array4<const amrex::Real>& cell_prim,
131  const amrex::Array4<const amrex::Real>& cell_data)
132 {
133  amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]);
134  amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);
135 
136  amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp),
137  cell_data(i,j,k ,RhoTheta_comp),
138  cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp));
139  amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
140  cell_data(i,j,k-1,RhoTheta_comp),
141  cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));
142 
143  amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
144  amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;
145 
146  amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
147  amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;
148 
149  amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
150  amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;
151 
152  amrex::Real qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) -
153  ( qc_plus - qc_d_ptr[k] +
154  qp_plus - qp_d_ptr[k] )
155  + (tempp3d-tempp1d)/tempp1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]);
156 
157  amrex::Real qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) -
158  ( qc_minus - qc_d_ptr[k-1] +
159  qp_minus - qp_d_ptr[k-1] )
160  + (tempm3d-tempm1d)/tempm1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]);
161 
162  amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
163  amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);
164 
165  return (-qavg * r0avg * grav_gpu);
166 }
167 
168 
169 AMREX_GPU_DEVICE
170 AMREX_FORCE_INLINE
171 amrex::Real
173  int& j,
174  int& k,
175  const int& n_qstate,
176  amrex::Real const& grav_gpu,
177  amrex::Real* rho_d_ptr,
178  amrex::Real* theta_d_ptr,
179  amrex::Real* qv_d_ptr,
180  const amrex::Array4<const amrex::Real>& cell_prim,
181  const amrex::Array4<const amrex::Real>& cell_data)
182 {
183  amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]);
184  amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);
185 
186  amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp),
187  cell_data(i,j,k ,RhoTheta_comp),
188  cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp));
189  amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
190  cell_data(i,j,k-1,RhoTheta_comp),
191  cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));
192 
193  amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
194  amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;
195 
196  amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
197  amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;
198 
199  amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
200  amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;
201 
202  amrex::Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d;
203 
204  amrex::Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d;
205 
206  amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
207  amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);
208 
209  return ( -qavg * r0avg * grav_gpu );
210 }
211 
212 
213 AMREX_GPU_DEVICE
214 AMREX_FORCE_INLINE
215 amrex::Real
217  int& j,
218  int& k,
219  const int& n_qstate,
220  amrex::Real const& grav_gpu,
221  amrex::Real* rho_d_ptr,
222  amrex::Real* theta_d_ptr,
223  amrex::Real* qv_d_ptr,
224  amrex::Real* qc_d_ptr,
225  amrex::Real* qp_d_ptr,
226  const amrex::Array4<const amrex::Real>& cell_prim,
227  const amrex::Array4<const amrex::Real>& cell_data)
228 {
229  amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
230  amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;
231 
232  amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
233  amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;
234 
235  amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
236  amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;
237 
238  amrex::Real qplus = amrex::Real(0.61) * ( qv_plus - qv_d_ptr[k] ) -
239  ( qc_plus - qc_d_ptr[k] +
240  qp_plus - qp_d_ptr[k] )
241  + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ];
242 
243  amrex::Real qminus = amrex::Real(0.61) * ( qv_minus - qv_d_ptr[k-1] ) -
244  ( qc_minus - qc_d_ptr[k-1] +
245  qp_minus - qp_d_ptr[k-1] )
246  + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1];
247 
248  amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
249  amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);
250 
251  return (-qavg * r0avg * grav_gpu);
252 }
253 #endif
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:144
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0)
Definition: ERF_EOS.H:42
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenPandTh(const amrex::Real P, const amrex::Real Th, const amrex::Real rdOcp)
Definition: ERF_EOS.H:30
#define PrimQ1_comp
Definition: ERF_IndexDefines.H:53
#define PrimQ2_comp
Definition: ERF_IndexDefines.H:54
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
#define RhoQ2_comp
Definition: ERF_IndexDefines.H:43
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42
#define PrimQ3_comp
Definition: ERF_IndexDefines.H:55
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_default(int &i, int &j, int &k, amrex::Real const &grav_gpu, amrex::Real const &rd_over_cp, const amrex::Array4< const amrex::Real > &p0_arr, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:73
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_type2(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, amrex::Real *rho_d_ptr, amrex::Real *theta_d_ptr, amrex::Real *qv_d_ptr, amrex::Real *qc_d_ptr, amrex::Real *qp_d_ptr, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:120
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_moist_anelastic(int &i, int &j, int &k, amrex::Real const &grav_gpu, amrex::Real const &rv_over_rd, const amrex::Array4< const amrex::Real > &p0_arr, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:37
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_type3(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, amrex::Real *rho_d_ptr, amrex::Real *theta_d_ptr, amrex::Real *qv_d_ptr, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:172
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_type4(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, amrex::Real *rho_d_ptr, amrex::Real *theta_d_ptr, amrex::Real *qv_d_ptr, amrex::Real *qc_d_ptr, amrex::Real *qp_d_ptr, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:216
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_type1(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:99
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_anelastic(int &i, int &j, int &k, amrex::Real const &grav_gpu, amrex::Real const &rd_over_cp, const amrex::Array4< const amrex::Real > &p0_arr, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_buoyancy_utils.H:10