REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_gls.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] lev level to operate on
7 * @param[inout] mf_gls turbulent generic length scale
8 * @param[inout] mf_tke turbulent kinetic energy
9 * @param[in ] mf_W vertical velocity
10 * @param[in ] mf_msku land-sea mask on u points
11 * @param[in ] mf_mskv land-sea mask on v points
12 * @param[in ] nstp index of last time step in gls and tke MultiFabs
13 * @param[in ] nnew index of time step to update in gls and tke MultiFabs
14 * @param[in ] iic which time step we're on
15 * @param[in ] ntfirst what is the first time step?
16 * @param[in ] N number of vertical levels
17 * @param[in ] dt_lev time step at this level
18 */
19void
20REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
21 MultiFab& mf_W, MultiFab* mf_msku, MultiFab* mf_mskv,
22 const int nstp, const int nnew,
23 const int iic, const int ntfirst, const int N, const Real dt_lev)
24{
25 // temps: grad, gradL, XF, FX, FXL, EF, FE, FEL
26 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
27 Array4<Real> const& gls = mf_gls->array(mfi);
28 Array4<Real> const& tke = mf_tke->array(mfi);
29 Array4<Real const> const& W = mf_W.const_array(mfi);
30
31 Array4<Real const> const& Huon = vec_Huon[lev]->const_array(mfi);
32 Array4<Real const> const& Hvom = vec_Hvom[lev]->const_array(mfi);
33 Array4<Real const> const& Hz = vec_Hz[lev]->const_array(mfi);
34 Array4<Real const> const& pm = vec_pm[lev]->const_array(mfi);
35 Array4<Real const> const& pn = vec_pn[lev]->const_array(mfi);
36 Array4<Real const> const& msku = mf_msku->const_array(mfi);
37 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
38
39 Box bx = mfi.tilebox();
40 Box xbx = surroundingNodes(bx,0);
41 Box ybx = surroundingNodes(bx,1);
42
43 Box xbx_hi = growHi(xbx,0,1);
44
45 Box ybx_hi = growHi(ybx,0,1);
46
47 const Box& domain = geom[0].Domain();
48 const auto dlo = amrex::lbound(domain);
49 const auto dhi = amrex::ubound(domain);
50
51 GeometryData const& geomdata = geom[0].data();
52 bool is_periodic_in_x = geomdata.isPeriodic(0);
53 bool is_periodic_in_y = geomdata.isPeriodic(1);
54
55 int ncomp = 1;
56 Vector<BCRec> bcrs_x(ncomp);
57 Vector<BCRec> bcrs_y(ncomp);
58 amrex::setBC(xbx,domain,BCVars::xvel_bc,0,1,domain_bcs_type,bcrs_x);
59 amrex::setBC(ybx,domain,BCVars::yvel_bc,0,1,domain_bcs_type,bcrs_y);
60
61 FArrayBox fab_XF(xbx_hi, 1, amrex::The_Async_Arena()); fab_XF.template setVal<RunOn::Device>(0.);
62 FArrayBox fab_FX(xbx_hi, 1, amrex::The_Async_Arena()); fab_FX.template setVal<RunOn::Device>(0.);
63 FArrayBox fab_FXL(xbx_hi, 1, amrex::The_Async_Arena()); fab_FXL.template setVal<RunOn::Device>(0.);
64 FArrayBox fab_EF(ybx_hi, 1, amrex::The_Async_Arena()); fab_EF.template setVal<RunOn::Device>(0.);
65 FArrayBox fab_FE(ybx_hi, 1, amrex::The_Async_Arena()); fab_FE.template setVal<RunOn::Device>(0.);
66 FArrayBox fab_FEL(ybx_hi, 1, amrex::The_Async_Arena()); fab_FEL.template setVal<RunOn::Device>(0.);
67 FArrayBox fab_Hz_half(bx, 1, amrex::The_Async_Arena()); fab_Hz_half.template setVal<RunOn::Device>(0.);
68 FArrayBox fab_CF(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_CF.template setVal<RunOn::Device>(0.);
69 FArrayBox fab_FC(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FC.template setVal<RunOn::Device>(0.);
70 FArrayBox fab_FCL(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FCL.template setVal<RunOn::Device>(0.);
71
72 auto XF = fab_XF.array();
73 auto FX = fab_FX.array();
74 auto FXL = fab_FXL.array();
75 auto EF = fab_EF.array();
76 auto FE = fab_FE.array();
77 auto FEL = fab_FEL.array();
78 auto Hz_half = fab_Hz_half.array();
79 auto CF = fab_CF.array();
80 auto FC = fab_FC.array();
81 auto FCL = fab_FCL.array();
82
83 // need XF/FX/FXL from [xlo to xhi] by [ylo to yhi ] on u points
84 ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
85 {
86 Real grad_im1 = (tke(i-1,j,k,nstp) - tke(i-2,j,k,nstp)) * msku(i-1,j,0);
87 Real grad_ip1 = (tke(i+1,j,k,nstp) - tke(i ,j,k,nstp)) * msku(i+1,j,0);
88
89 Real gradL_im1 = (gls(i-1,j,k,nstp) - gls(i-2,j,k,nstp)) * msku(i-1,j,0);
90 Real gradL_ip1 = (gls(i+1,j,k,nstp) - gls(i ,j,k,nstp)) * msku(i+1,j,0);
91
92 // Adjust boundaries
93 // TODO: Make sure indices match with what ROMS does
94 if (i == dlo.x-1 && !is_periodic_in_x) {
95 grad_im1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp);
96 gradL_im1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp);
97 }
98 else if (i == dhi.x+1 && !is_periodic_in_x) {
99 grad_ip1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp);
100 gradL_ip1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp);
101 }
102 Real cff = 1.0_rt/6.0_rt;
103 XF(i,j,k) = 0.5_rt * (Huon(i,j,k) + Huon(i,j,k-1));
104 FX(i,j,k) = XF(i,j,k) * 0.5_rt * (tke(i-1,j,k,nstp) + tke(i,j,k,nstp) -
105 cff * (grad_ip1 - grad_im1));
106 FXL(i,j,k) = XF(i,j,k) * 0.5_rt * (gls(i-1,j,k,nstp) + gls(i,j,k,nstp) -
107 cff * (gradL_ip1 - gradL_im1));
108 });
109
110 // need EF/FE/FEL from [xlo to xhi ] by [ylo to yhi+1]
111 ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
112 {
113 Real grad_jm1 = (tke(i,j-1,k,nstp) - tke(i,j-2,k,nstp)) * mskv(i,j-1,0);
114 Real grad_jp1 = (tke(i,j+1,k,nstp) - tke(i,j ,k,nstp)) * mskv(i,j+1,0);
115
116 Real gradL_jm1 = (gls(i,j-1,k,nstp) - gls(i,j-2,k,nstp)) * mskv(i,j-1,0);
117 Real gradL_jp1 = (gls(i,j+1,k,nstp) - gls(i,j ,k,nstp)) * mskv(i,j+1,0);
118
119 // Adjust boundaries
120 // TODO: Make sure indices match with what ROMS does
121 if (j == dlo.y-1 && !is_periodic_in_y) {
122 grad_jm1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp);
123 gradL_jm1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp);
124 }
125 else if (j == dhi.y+1 && !is_periodic_in_y) {
126 grad_jp1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp);
127 gradL_jp1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp);
128 }
129 Real cff = 1.0_rt/6.0_rt;
130 EF(i,j,k) = 0.5_rt * (Hvom(i,j,k) + Hvom(i,j,k-1));
131 FE(i,j,k) = EF(i,j,k) * 0.5_rt * (tke(i,j-1,k,nstp) + tke(i,j,k,nstp) -
132 cff * (grad_jp1 - grad_jm1));
133 FEL(i,j,k) = EF(i,j,k) * 0.5_rt * (gls(i,j-1,k,nstp) + gls(i,j,k,nstp) -
134 cff * (gradL_jp1 - gradL_jm1));
135 });
136
137 Real gamma = 1.0_rt / 6.0_rt;
138 Real cff1, cff2, cff3;
139 int indx;
140 // Time step horizontal advection
141 if (iic == ntfirst) {
142 cff1 = 1.0_rt;
143 cff2 = 0.0_rt;
144 cff3 = 0.5_rt * dt_lev;
145 indx = nstp;
146 } else {
147 cff1 = 0.5_rt + gamma;
148 cff2 = 0.5_rt - gamma;
149 cff3 = (1.0_rt - gamma) * dt_lev;
150 indx = 1 - nstp;
151 }
152
153 // update tke, gls from [xlo to xhi ] by [ylo to yhi ]
154 // need XF/FX/FXL from [xlo to xhi+1] by [ylo to yhi ]
155 // need EF/FE/FEL from [xlo to xhi ] by [ylo to yhi+1]
156 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
157 {
158 Real cff = 0.5_rt * (Hz(i,j,k) + Hz(i,j,k-1));
159 Real cff4 = cff3 * pm(i,j,0) * pn(i,j,0);
160 Hz_half(i,j,k) = cff - cff4 * (XF(i+1,j,k)-XF(i,j,k)+EF(i,j+1,k)-EF(i,j,k));
161 tke(i,j,k,2) = cff * (cff1*tke(i,j,k,nstp) + cff2*tke(i,j,k,indx)) -
162 cff4 * (FX(i+1,j,k)-FX(i,j,k)+FE(i,j+1,k)-FE(i,j,k));
163 gls(i,j,k,2) = cff * (cff1 * gls(i,j,k,nstp) + cff2 * gls(i,j,k,indx)) -
164 cff4 * (FXL(i+1,j,k)-FXL(i,j,k)+FEL(i,j+1,k)-FEL(i,j,k));
165 tke(i,j,k,nnew) = cff * tke(i,j,k,nstp);
166 gls(i,j,k,nnew) = cff * gls(i,j,k,nstp);
167 });
168
169 // Will do a FillPatch after this, so don't need to do any ghost zones in x,y
170 // Compute vertical advection
171 ParallelFor(convert(bx,IntVect(0,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
172 {
173 // CF and FC/FCL are on rho points
174 CF(i,j,k) = 0.5_rt * (W(i,j,k+1) + W(i,j,k));
175 if (k == 0) {
176 Real cff1_vadv = 1.0_rt / 3.0_rt;
177 Real cff2_vadv = 5.0_rt / 6.0_rt;
178 Real cff3_vadv = 1.0_rt / 6.0_rt;
179 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * tke(i,j,0,nstp) +
180 cff2_vadv * tke(i,j,1,nstp) -
181 cff3_vadv * tke(i,j,2,nstp));
182 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * gls(i,j,0,nstp) +
183 cff2_vadv * gls(i,j,1,nstp) -
184 cff3_vadv * gls(i,j,2,nstp));
185 } else if (k == N) {
186 Real cff1_vadv = 1.0_rt / 3.0_rt;
187 Real cff2_vadv = 5.0_rt / 6.0_rt;
188 Real cff3_vadv = 1.0_rt / 6.0_rt;
189 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * tke(i,j,k+1, nstp) +
190 cff2_vadv * tke(i,j,k ,nstp)-
191 cff3_vadv * tke(i,j,k-1,nstp));
192 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * gls(i,j,k+1,nstp) +
193 cff2_vadv * gls(i,j,k ,nstp)-
194 cff3_vadv * gls(i,j,k-1,nstp));
195 } else {
196 Real cff1_vadv = 7.0_rt / 12.0_rt;
197 Real cff2_vadv = 1.0_rt / 12.0_rt;
198 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * (tke(i,j,k ,nstp) + tke(i,j,k+1,nstp)) -
199 cff2_vadv * (tke(i,j,k-1,nstp) + tke(i,j,k+2,nstp)));
200 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * (gls(i,j,k ,nstp) + gls(i,j,k+1,nstp)) -
201 cff2_vadv * (gls(i,j,k-1,nstp) + gls(i,j,k+2,nstp)));
202 }
203 });
204
205 // Time-step vertical advection
206 if (iic == ntfirst) {
207 cff3 = 0.5_rt * dt_lev;
208 } else {
209 cff3 = (1.0_rt - gamma) * dt_lev;
210 }
211 // DO k=1,N-1
212 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
213 {
214 Real cff4 = cff3 * pm(i,j,0) * pn(i,j,0);
215 Hz_half(i,j,k) = Hz_half(i,j,k) - cff4 * (CF(i,j,k)-CF(i,j,k-1));
216 Real cff1_loc = 1.0_rt / Hz_half(i,j,k);
217 tke(i,j,k,2) = cff1_loc * (tke(i,j,k,2) - cff4 * (FC (i,j,k) - FC (i,j,k-1)));
218 gls(i,j,k,2) = cff1_loc * (gls(i,j,k,2) - cff4 * (FCL(i,j,k) - FCL(i,j,k-1)));
219 });
220 }
221
222 for (int icomp=0; icomp<3; icomp++) {
223 FillPatch(lev, t_old[lev], *vec_tke[lev], GetVecOfPtrs(vec_tke), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
224 FillPatch(lev, t_old[lev], *vec_gls[lev], GetVecOfPtrs(vec_gls), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
225 }
226}
227
228/**
229 * @param[in ] lev level to operate on
230 * @param[inout] mf_gls turbulent generic length scale
231 * @param[inout] mf_tke turbulent kinetic energy
232 * @param[in ] mf_W vertical velocity
233 * @param[inout] mf_Akv vertical viscosity coefficient
234 * @param[inout] mf_Akt vertical diffusivity coefficients
235 * @param[inout] mf_Akk turbulent kinetic energy vertical diffusion coefficient
236 * @param[inout] mf_Akp turbulent length scale vertical diffusion coefficient
237 * @param[in ] mf_mskr land-sea mask on rho points
238 * @param[in ] mf_msku land-sea mask on u points
239 * @param[in ] mf_mskv land-sea mask on v points
240 * @param[in ] nstp index of last time step in gls and tke MultiFabs
241 * @param[in ] nnew index of time step to update in gls and tke MultiFabs
242 * @param[in ] N number of vertical levels
243 * @param[in ] dt_lev time step at this level
244 */
245void
246REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
247 MultiFab& mf_W, MultiFab* mf_Akv, MultiFab* mf_Akt,
248 MultiFab* mf_Akk, MultiFab* mf_Akp,
249 MultiFab* mf_mskr,
250 MultiFab* mf_msku, MultiFab* mf_mskv,
251 const int nstp, const int nnew,
252 const int N, const Real dt_lev)
253{
254//-----------------------------------------------------------------------
255// Compute several constants.
256//-----------------------------------------------------------------------
257 bool Lmy25 = ((solverChoice.gls_p == 0.0) &&
258 (solverChoice.gls_n == 1.0) &&
259 (solverChoice.gls_m == 1.0)) ? true : false;
260
261 Real L_sft = vonKar;
262 Real gls_sigp_cb = solverChoice.gls_sigp;
263 Real ogls_sigp = 1.0_rt/gls_sigp_cb;
264
265 Real gls_c3m = solverChoice.gls_c3m;
266 Real gls_c3p = solverChoice.gls_c3p;
267 Real gls_cmu0 = solverChoice.gls_cmu0;
268
269 Real gls_m = solverChoice.gls_m;
270 Real gls_n = solverChoice.gls_n;
271 Real gls_p = solverChoice.gls_p;
272
273 Real gls_Gh0 = solverChoice.gls_Gh0;
274 Real gls_Ghcri = solverChoice.gls_Ghcri;
275 Real gls_Ghmin = solverChoice.gls_Ghmin;
276
277 Real Akv_bak = solverChoice.Akv_bak;
278 Real Akt_bak = solverChoice.Akt_bak;
279 Real Akp_bak = solverChoice.Akp_bak;
280 Real Akk_bak = solverChoice.Akk_bak;
281
282 Real gls_c1 = solverChoice.gls_c1;
283 Real gls_c2 = solverChoice.gls_c2;
284 Real gls_E2 = solverChoice.gls_E2;
285 Real gls_sigk = solverChoice.gls_sigk;
286 auto gls_stability_type = solverChoice.gls_stability_type;
287
288 Real sqrt2 = std::sqrt(2.0_rt);
289 Real cmu_fac1 = std::pow(solverChoice.gls_cmu0,(-solverChoice.gls_p/solverChoice.gls_n));
290 Real cmu_fac2 = std::pow(solverChoice.gls_cmu0,(3.0_rt+solverChoice.gls_p/solverChoice.gls_n));
291 Real cmu_fac3 = 1.0_rt/std::pow(solverChoice.gls_cmu0,2.0_rt);
292 //Real cmu_fac4 = std::pow(1.5_rt*solverChoice.gls_sigk,(1.0_rt/3.0_rt))/std::pow(solverChoice.gls_cmu0,4.0_rt/3.0_rt);
293
294 //Real gls_fac1 = solverChoice.gls_n*std::pow(solverChoice.gls_cmu0,solverChoice.gls_p+1.0_rt);
295 Real gls_fac2 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p)*solverChoice.gls_n*std::pow(vonKar,solverChoice.gls_n);
296 Real gls_fac3 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p)*solverChoice.gls_n;
297 Real gls_fac4 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p);
298 Real gls_fac5 = std::pow(0.56_rt,0.5_rt*solverChoice.gls_n)*std::pow(solverChoice.gls_cmu0,solverChoice.gls_p);
299 Real gls_fac6 = 8.0_rt/std::pow(solverChoice.gls_cmu0,6.0_rt);
300
301 Real gls_exp1 = 1.0_rt/solverChoice.gls_n;
302 Real tke_exp1 = solverChoice.gls_m/solverChoice.gls_n;
303 Real tke_exp2 = 0.5_rt+solverChoice.gls_m/solverChoice.gls_n;
304 //Real tke_exp3 = 0.5_rt+solverChoice.gls_m;
305 Real tke_exp4 = solverChoice.gls_m+0.5_rt*solverChoice.gls_n;
306
307 Real gls_s0, gls_s1, gls_s2, gls_s4, gls_s5, gls_s6;
308 Real gls_b0, gls_b1, gls_b2, gls_b3, gls_b4, gls_b5;
309 Real my_Sm2, my_Sh1, my_Sh2, my_Sm3, my_Sm4;
310
311 // Compute parameters for Canuto et al. (2001) stability functions.
312 // (Canuto, V.M., Cheng, H.Y., and Dubovikov, M.S., 2001: Ocean
313 // turbulence. Part I: One-point closure model - momentum and
314 // heat vertical diffusivities, JPO, 1413-1426).
315
316 if (solverChoice.gls_stability_type == GLS_StabilityType::Canuto_A ||
317 solverChoice.gls_stability_type == GLS_StabilityType::Canuto_B) {
318
319 gls_s0=3.0_rt/2.0_rt*solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L5;
320 gls_s1=-solverChoice.gls_L4*(solverChoice.gls_L6+solverChoice.gls_L7)
321 +2.0_rt*solverChoice.gls_L4*solverChoice.gls_L5*
322 (solverChoice.gls_L1-1.0_rt/3.0_rt*solverChoice.gls_L2-solverChoice.gls_L3)
323 +3.0_rt/2.0_rt*
324 solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L8;
325 gls_s2=-3.0_rt/8.0_rt*solverChoice.gls_L1
326 *(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
327 gls_s4=2.0_rt*solverChoice.gls_L5;
328 gls_s5=2.0_rt*solverChoice.gls_L4;
329 gls_s6=2.0_rt/3.0_rt*solverChoice.gls_L5
330 *(3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
331 1.0_rt/2.0_rt*solverChoice.gls_L5*solverChoice.gls_L1*(3.0_rt*solverChoice.gls_L3-solverChoice.gls_L2)+
332 3.0_rt/4.0_rt*solverChoice.gls_L1*(solverChoice.gls_L6-solverChoice.gls_L7);
333 gls_b0=3.0_rt*solverChoice.gls_L5*solverChoice.gls_L5;
334 gls_b1=solverChoice.gls_L5*(7.0_rt*solverChoice.gls_L4+3.0_rt*solverChoice.gls_L8);
335 gls_b2=solverChoice.gls_L5*solverChoice.gls_L5*(3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
336 3.0_rt/4.0_rt*(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
337 gls_b3=solverChoice.gls_L4*(4.0_rt*solverChoice.gls_L4+3.0_rt*solverChoice.gls_L8);
338 gls_b5=1.0_rt/4.0_rt*(solverChoice.gls_L2*solverChoice.gls_L2-3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3)*
339 (solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
340 gls_b4=solverChoice.gls_L4*(solverChoice.gls_L2*solverChoice.gls_L6-3.0_rt*solverChoice.gls_L3*solverChoice.gls_L7-
341 solverChoice.gls_L5*(solverChoice.gls_L2*solverChoice.gls_L2-solverChoice.gls_L3*solverChoice.gls_L3))+solverChoice.gls_L5*solverChoice.gls_L8*
342 (3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2);
343 my_Sm2 = 0.0_rt;
344 my_Sm3 = 0.0_rt;
345 my_Sm4 = 0.0_rt;
346 my_Sh1 = 0.0_rt;
347 my_Sh2 = 0.0_rt;
348 } else {
349 gls_s0 = 0.0_rt;
350 gls_s1 = 0.0_rt;
351 gls_s2 = 0.0_rt;
352 gls_s4 = 0.0_rt;
353 gls_s5 = 0.0_rt;
354 gls_s6 = 0.0_rt;
355 gls_b0 = 0.0_rt;
356 gls_b1 = 0.0_rt;
357 gls_b2 = 0.0_rt;
358 gls_b3 = 0.0_rt;
359 gls_b4 = 0.0_rt;
360 gls_b5 = 0.0_rt;
361 //my_Sm1=solverChoice.my_A1*solverChoice.my_A2*((solverChoice.my_B2-3.0_rt*solverChoice.my_A2)*
362 // (1.0_rt-6.0_rt*solverChoice.my_A1/solverChoice.my_B1)-
363 // 3.0_rt*solverChoice.my_C1*(solverChoice.my_B2+6.0_rt*solverChoice.my_A1));
364 my_Sm2=9.0_rt*solverChoice.my_A1*solverChoice.my_A2;
365 my_Sm3=solverChoice.my_A1*(1.0_rt-3.0_rt*solverChoice.my_C1-6.0_rt*solverChoice.my_A1/solverChoice.my_B1);
366 my_Sm4=18.0_rt*solverChoice.my_A1*solverChoice.my_A1+9.0_rt*solverChoice.my_A1*solverChoice.my_A2;
367 my_Sh1=solverChoice.my_A2*(1.0_rt-6.0_rt*solverChoice.my_A1/solverChoice.my_B1);
368 my_Sh2=3.0_rt*solverChoice.my_A2*(6.0_rt*solverChoice.my_A1+solverChoice.my_B2);
369 }
370
371 Real Zos_min = std::max(solverChoice.Zos, 0.0001_rt);
372 Real Zos_eff = Zos_min;
373 Real Gadv = 1.0_rt/3.0_rt;
374 Real eps = 1.0e-10_rt;
375
376 const BoxArray& ba = cons_old[lev]->boxArray();
377 const DistributionMapping& dm = cons_old[lev]->DistributionMap();
378 MultiFab mf_dU(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0));
379 MultiFab mf_dV(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0));
380
381 MultiFab mf_CF(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0));
382
383 MultiFab mf_shear2(ba,dm,1,IntVect(NGROW,NGROW,0));
384 MultiFab mf_shear2_cached(ba,dm,1,IntVect(NGROW,NGROW,0));
385
386 MultiFab mf_buoy2(ba,dm,1,IntVect(NGROW,NGROW,0));
387 MultiFab mf_tmp_buoy(ba,dm,1,IntVect(NGROW,NGROW,0));
388 MultiFab mf_tmp_shear(ba,dm,1,IntVect(NGROW,NGROW,0));
389 MultiFab mf_curvK(ba,dm,1,IntVect(NGROW,NGROW,0));
390 MultiFab mf_curvP(ba,dm,1,IntVect(NGROW,NGROW,0));
391
392 MultiFab mf_FXK(ba,dm,1,IntVect(NGROW,NGROW,0));
393 MultiFab mf_FXP(ba,dm,1,IntVect(NGROW,NGROW,0));
394 MultiFab mf_FEK(ba,dm,1,IntVect(NGROW,NGROW,0));
395 MultiFab mf_FEP(ba,dm,1,IntVect(NGROW,NGROW,0));
396 MultiFab mf_FCK(ba,dm,1,IntVect(NGROW,NGROW,0));
397 MultiFab mf_FCP(ba,dm,1,IntVect(NGROW,NGROW,0));
398 MultiFab mf_BCK(ba,dm,1,IntVect(NGROW,NGROW,0));
399 MultiFab mf_BCP(ba,dm,1,IntVect(NGROW,NGROW,0));
400
401 const Box& domain = geom[0].Domain();
402 const auto dlo = amrex::lbound(domain);
403 const auto dhi = amrex::ubound(domain);
404
405 GeometryData const& geomdata = geom[0].data();
406 bool is_periodic_in_x = geomdata.isPeriodic(0);
407 bool is_periodic_in_y = geomdata.isPeriodic(1);
408
409
410 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
411 {
412 Box bx = mfi.tilebox();
413 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
414
415 Box bxD = bx;
416 bxD.makeSlab(2,0);
417 Box gbx1D = gbx1;
418 gbx1D.makeSlab(2,0);
419
420 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
421 Array4<Real> const& u = xvel_old[lev]->array(mfi);
422 Array4<Real> const& v = yvel_old[lev]->array(mfi);
423
424 auto dU = mf_dU.array(mfi);
425 auto dV = mf_dV.array(mfi);
426 auto CF = mf_CF.array(mfi);
427 auto shear2_cached = mf_shear2_cached.array(mfi);
428
429 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
430 {
431 CF(i,j,0) = 0.0_rt;
432 dU(i,j,0) = 0.0_rt;
433 dV(i,j,0) = 0.0_rt;
434 for (int k=1; k<=N; k++) {
435 Real cff = 1.0_rt / (2.0_rt * Hz(i,j,k) + Hz(i,j,k-1)*(2.0_rt - CF(i,j,k-1)));
436 CF(i,j,k) = cff * Hz(i,j,k);
437 dU(i,j,k)=cff*(3.0_rt*(u(i ,j,k)-u(i, j,k-1)+
438 u(i+1,j,k)-u(i+1,j,k-1))-Hz(i,j,k-1)*dU(i,j,k-1));
439 dV(i,j,k)=cff*(3.0_rt*(v(i,j ,k)-v(i,j ,k-1)+
440 v(i,j+1,k)-v(i,j+1,k-1))-Hz(i,j,k-1)*dV(i,j,k-1));
441 }
442 dU(i,j,N+1) = 0.0_rt;
443 dV(i,j,N+1) = 0.0_rt;
444 for (int k=N; k>=1; k--) {
445 dU(i,j,k) = dU(i,j,k) - CF(i,j,k) * dU(i,j,k+1);
446 dV(i,j,k) = dV(i,j,k) - CF(i,j,k) * dV(i,j,k+1);
447 }
448 shear2_cached(i,j,0) = 0.0_rt;
449 for (int k=1; k<=N; k++) {
450 shear2_cached(i,j,k) = dU(i,j,k) * dU(i,j,k) + dV(i,j,k) * dV(i,j,k);
451 }
452 });
453 }
454
455 // While potentially counterintuitive, this is what ROMS does for handling shear2 at all boundaries, even
456 // periodic
457 (*physbcs[lev])(mf_shear2_cached,*mf_mskr,0,1,mf_shear2_cached.nGrowVect(),t_new[lev],BCVars::foextrap_bc);
458 mf_CF.setVal(0.0_rt);
459
460 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
461 {
462 Box bx = mfi.tilebox();
463 Box xbx = surroundingNodes(bx,0);
464 Box ybx = surroundingNodes(bx,1);
465 Box gbx1 = grow(bx,IntVect(NGROW-1,NGROW-1,0));
466
467 Box bx_rho = bx;
468 bx_rho.convert(IntVect(0,0,0));
469
470 Box xbx_int = grow(xbx,IntVect(0,0,-1)); // interior points in w of xbx
471 Box ybx_int = grow(ybx,IntVect(0,0,-1)); // interior points in w of ybx
472 Box bx_growloxy = growLo(growLo(grow(bx,IntVect(0,0,-1)),0,1),1,1);
473
474 Box bxD = bx;
475 bxD.makeSlab(2,0);
476 Box gbx1D = gbx1;
477 gbx1D.makeSlab(2,0);
478
479 int ncomp = 1;
480 Vector<BCRec> bcrs_x(ncomp);
481 Vector<BCRec> bcrs_y(ncomp);
482 amrex::setBC(xbx,domain,BCVars::xvel_bc,0,1,domain_bcs_type,bcrs_x);
483 amrex::setBC(ybx,domain,BCVars::yvel_bc,0,1,domain_bcs_type,bcrs_y);
484
485 Array4<Real const> const& W = mf_W.const_array(mfi);
486 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
487 Array4<Real> const& pm = vec_pm[lev]->array(mfi);
488 Array4<Real> const& pn = vec_pn[lev]->array(mfi);
489 Array4<Real> const& Lscale = vec_Lscale[lev]->array(mfi);
490
491 Array4<Real> const& Huon = vec_Huon[lev]->array(mfi);
492 Array4<Real> const& Hvom = vec_Hvom[lev]->array(mfi);
493 Array4<Real> const& z_w = vec_z_w[lev]->array(mfi);
494
495 Array4<Real> const& tke = mf_tke->array(mfi);
496 Array4<Real> const& gls = mf_gls->array(mfi);
497
498 Array4<Real const> const& sustr = vec_sustr[lev]->const_array(mfi);
499 Array4<Real const> const& svstr = vec_svstr[lev]->const_array(mfi);
500 Array4<Real const> const& bustr = vec_bustr[lev]->const_array(mfi);
501 Array4<Real const> const& bvstr = vec_bvstr[lev]->const_array(mfi);
502 Array4<Real const> const& msku = mf_msku->const_array(mfi);
503 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
504
505 Array4<Real> const& ZoBot = vec_ZoBot[lev]->array(mfi);
506
507 FArrayBox fab_tmp_buoy(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_buoy.template setVal<RunOn::Device>(0.);
508 FArrayBox fab_tmp_shear(bx_growloxy,1, amrex::The_Async_Arena()); fab_tmp_shear.template setVal<RunOn::Device>(0.);
509
510 FArrayBox fab_curvK(gbx1,1, amrex::The_Async_Arena()); fab_curvK.template setVal<RunOn::Device>(0.);
511 FArrayBox fab_curvP(gbx1,1, amrex::The_Async_Arena()); fab_curvP.template setVal<RunOn::Device>(0.);
512
513 FArrayBox fab_FXK(xbx_int,1,amrex::The_Async_Arena()); fab_FXK.template setVal<RunOn::Device>(0.);
514 FArrayBox fab_FXP(xbx_int,1,amrex::The_Async_Arena()); fab_FXP.template setVal<RunOn::Device>(0.);
515 FArrayBox fab_FEK(ybx_int,1,amrex::The_Async_Arena()); fab_FEK.template setVal<RunOn::Device>(0.);
516 FArrayBox fab_FEP(ybx_int,1,amrex::The_Async_Arena()); fab_FEP.template setVal<RunOn::Device>(0.);
517
518 FArrayBox fab_FCK(bx_rho,1,amrex::The_Async_Arena()); fab_FCK.template setVal<RunOn::Device>(0.);
519 FArrayBox fab_FCP(bx_rho,1,amrex::The_Async_Arena()); fab_FCP.template setVal<RunOn::Device>(0.);
520 FArrayBox fab_BCK(bx_rho,1,amrex::The_Async_Arena()); fab_BCK.template setVal<RunOn::Device>(0.);
521 FArrayBox fab_BCP(bx_rho,1,amrex::The_Async_Arena()); fab_BCP.template setVal<RunOn::Device>(0.);
522
523 auto CF = mf_CF.array(mfi);
524 auto shear2 = mf_shear2.array(mfi);
525 auto shear2_cached = mf_shear2_cached.array(mfi);
526 auto buoy2 = mf_buoy2.array(mfi);
527 Array4<Real> const& bvf = vec_bvf[lev]->array(mfi);
528
529 auto tmp_buoy = fab_tmp_buoy.array();
530 auto tmp_shear = fab_tmp_shear.array();
531 auto curvK = fab_curvK.array();
532 auto curvP = fab_curvP.array();
533 auto FXK = fab_FXK.array();
534 auto FXP = fab_FXP.array();
535 auto FEK = fab_FEK.array();
536 auto FEP = fab_FEP.array();
537 auto FCK = fab_FCK.array();
538 auto FCP = fab_FCP.array();
539 auto BCK = fab_BCK.array();
540 auto BCP = fab_BCP.array();
541
542 auto Akt = mf_Akt->array(mfi);
543 auto Akv = mf_Akv->array(mfi);
544 auto Akp = mf_Akp->array(mfi);
545 auto Akk = mf_Akk->array(mfi);
546
547 ParallelFor(bx_growloxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
548 {
549 tmp_buoy(i,j,k)=0.25_rt * (bvf(i,j,k) + bvf(i+1,j,k) + bvf(i,j+1,k)+bvf(i+1,j+1,k));
550 tmp_shear(i,j,k)=0.25_rt * (shear2_cached(i,j,k) + shear2_cached(i+1,j,k) + shear2_cached(i,j+1,k)+shear2_cached(i+1,j+1,k));
551 });
552
553 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
554 {
555 buoy2(i,j,k)=0.25_rt * (tmp_buoy(i,j,k) + tmp_buoy(i-1,j,k) + tmp_buoy(i,j-1,k)+tmp_buoy(i-1,j-1,k));
556 shear2(i,j,k)=0.25_rt * (tmp_shear(i,j,k) + tmp_shear(i-1,j,k) + tmp_shear(i,j-1,k)+tmp_shear(i-1,j-1,k));
557 });
558
559 //Time step advective terms
560 ParallelFor(growLo(grow(xbx,IntVect(0,0,-1)),0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
561 {
562 Real gradK, gradK_ip1, gradP, gradP_ip1;
563
564 if (i == dlo.x-1 && !is_periodic_in_x) {
565 gradK_ip1 = tke(i+1,j,k,2)-tke(i ,j,k,2);
566 gradK = gradK_ip1;
567 gradP_ip1 = gls(i+1,j,k,2)-gls(i ,j,k,2);
568 gradP = gradP_ip1;
569 } else if (i == dhi.x+1 && !is_periodic_in_x) {
570 gradK = tke(i ,j,k,2)-tke(i-1,j,k,2);
571 gradK_ip1 = gradK;
572 gradP = gls(i ,j,k,2)-gls(i-1,j,k,2);
573 gradP_ip1 = gradP;
574 } else {
575 gradK = (tke(i ,j,k,2)-tke(i-1,j,k,2)) * msku(i ,j,0);
576 gradK_ip1 = (tke(i+1,j,k,2)-tke(i ,j,k,2)) * msku(i+1,j,0);
577 gradP = (gls(i ,j,k,2)-gls(i-1,j,k,2)) * msku(i ,j,0);
578 gradP_ip1 = (gls(i+1,j,k,2)-gls(i ,j,k,2)) * msku(i+1,j,0);
579 }
580
581 curvK(i,j,k) = gradK_ip1 - gradK;
582 curvP(i,j,k) = gradP_ip1 - gradP;
583 });
584 ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
585 {
586 Real cff = 0.5_rt * (Huon(i,j,k) + Huon(i,j,k-1));
587 Real cff1 = (cff > 0.0) ? curvK(i-1,j,k) : curvK(i,j,k);
588 Real cff2 = (cff > 0.0) ? curvP(i-1,j,k) : curvP(i,j,k);
589
590 FXK(i,j,k) = cff * 0.5_rt * (tke(i-1,j,k,2)+tke(i,j,k,2)-Gadv*cff1);
591 FXP(i,j,k) = cff * 0.5_rt * (gls(i-1,j,k,2)+gls(i,j,k,2)-Gadv*cff2);
592 });
593
594 //Time step advective terms
595 ParallelFor(growLo(grow(ybx,IntVect(0,0,-1)),1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
596 {
597 Real gradK = (tke(i,j ,k,2)-tke(i,j-1,k,2)) * mskv(i,j ,0);
598 Real gradK_jp1 = (tke(i,j+1,k,2)-tke(i,j ,k,2)) * mskv(i,j+1,0);
599 Real gradP = (gls(i,j ,k,2)-gls(i,j-1,k,2)) * mskv(i,j ,0);
600 Real gradP_jp1 = (gls(i,j+1,k,2)-gls(i,j ,k,2)) * mskv(i,j+1,0);
601
602 if (j == dlo.y-1 && !is_periodic_in_y) {
603 gradK = gradK_jp1;
604 gradP = gradP_jp1;
605 }
606 else if (j == dhi.y+1 && !is_periodic_in_y) {
607 gradK_jp1 = gradK;
608 gradP_jp1 = gradP;
609 }
610
611 curvK(i,j,k) = gradK_jp1 - gradK;
612 curvP(i,j,k) = gradP_jp1 - gradP;
613 });
614 ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
615 {
616 Real cff = 0.5_rt * (Hvom(i,j,k) + Hvom(i,j,k-1));
617 Real cff1 = (cff > 0.0) ? curvK(i,j-1,k) : curvK(i,j,k);
618 Real cff2 = (cff > 0.0) ? curvP(i,j-1,k) : curvP(i,j,k);
619
620 FEK(i,j,k) = cff * 0.5_rt * (tke(i,j-1,k,2)+tke(i,j,k,2)-Gadv*cff1);
621 FEP(i,j,k) = cff * 0.5_rt * (gls(i,j-1,k,2)+gls(i,j,k,2)-Gadv*cff2);
622 });
623
624 Real gls_Kmin = solverChoice.gls_Kmin;
625 Real gls_Pmin = solverChoice.gls_Pmin;
626 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
627 {
628 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
629 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff * (FXK(i+1,j ,k)-FXK(i,j,k)+
630 FEK(i ,j+1,k)-FEK(i,j,k));
631 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
632
633 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff * (FXP(i+1,j ,k)-FXP(i,j,k)+
634 FEP(i ,j+1,k)-FEP(i,j,k));
635 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew), gls_Pmin);
636 });
637
638
639 // Vertical advection
640 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
641 {
642 Real cff1 = 7.0_rt / 12.0_rt;
643 Real cff2 = 1.0_rt / 12.0_rt;
644 for (int k=1; k<=N-1; k++) {
645 Real cff = 0.5_rt * (W(i,j,k+1)+W(i,j,k));
646 FCK(i,j,k) = cff * (cff1 * (tke(i,j,k ,2)+tke(i,j,k+1,2))-
647 cff2 * (tke(i,j,k-1,2)+tke(i,j,k+2,2)));
648 FCP(i,j,k) = cff * (cff1 * (gls(i,j,k ,2)+gls(i,j,k+1,2))-
649 cff2 * (gls(i,j,k-1,2)+gls(i,j,k+2,2)));
650 }
651 cff1 = 1.0_rt/3.0_rt;
652 cff2 = 5.0_rt/6.0_rt;
653 Real cff3 = 1.0_rt / 6.0_rt;
654 Real cff = 0.5_rt * (W(i,j,0)+W(i,j,1));
655 FCK(i,j,0) = cff * (cff1 * tke(i,j,0,2)+cff2 * tke(i,j,1,2)-cff3 * tke(i,j,2,2));
656 FCP(i,j,0) = cff * (cff1 * gls(i,j,0,2)+cff2 * gls(i,j,1,2)-cff3 * gls(i,j,2,2));
657 cff = 0.5_rt * (W(i,j,N+1)+W(i,j,N));
658 FCK(i,j,N) = cff * (cff1 * tke(i,j,N+1,2)+cff2*tke(i,j,N,2)-cff3*tke(i,j,N-1,2));
659 FCP(i,j,N) = cff * (cff1 * gls(i,j,N+1,2)+cff2*gls(i,j,N,2)-cff3*gls(i,j,N-1,2));
660 });
661 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
662 {
663 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
664 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff*(FCK(i,j,k )-FCK(i,j,k-1));
665 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
666 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff*(FCP(i,j,k )-FCP(i,j,k-1));
667 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
668 });
669
670 // Compute vertical mixing, turbulent production and turbulent
671 // dissipation.
672 //
673 Real cff = -0.5 * dt_lev;
674 ParallelFor(convert(bx,IntVect(0,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
675 {
676 if (k==0 or k==N) {
677 FCK(i,j,k) = 0.0_rt;
678 FCP(i,j,k) = 0.0_rt;
679 } else {
680 FCK(i,j,k) = cff * (Akk(i,j,k) + Akk(i,j,k+1)) / Hz(i,j,k);
681 FCP(i,j,k) = cff * (Akp(i,j,k) + Akp(i,j,k+1)) / Hz(i,j,k);
682 }
683 });
684 // Compute production and dissipation terms.
685 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
686 {
687 // Compute shear and buoyant production of turbulent energy (m3/s3)
688 // at W-points (ignore small negative values of buoyancy).
689 Real strat2 = buoy2(i,j,k);
690 Real gls_c3 = (strat2 > 0.0) ? gls_c3m : gls_c3p;
691 Real Kprod = shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
692 strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak);
693 Real Pprod = gls_c1 * shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
694 gls_c3 * strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak);
695
696 // If negative production terms, then add buoyancy to dissipation terms
697 // (BCK and BCP) below, using "cff1" and "cff2" as the on/off switch.
698 Real cff1 = (Kprod < 0.0_rt) ? 0.0_rt : 1.0_rt;
699 Real cff2 = (Pprod < 0.0_rt) ? 0.0_rt : 1.0_rt;
700 if (Kprod < 0.0_rt) {
701 Kprod = Kprod + strat2*(Akt(i,j,k,Temp_comp)-Akt_bak);
702 }
703 if (Pprod < 0.0_rt) {
704 Pprod = Pprod + gls_c3*strat2*(Akt(i,j,k,Temp_comp)-Akt_bak);
705 }
706 // Time-step shear and buoyancy production terms.
707 Real cff_Hz = 0.5_rt * (Hz(i,j,k) + Hz(i,j,k-1));
708 tke(i,j,k,nnew) = tke(i,j,k,nnew)+dt_lev * cff_Hz * Kprod;
709 gls(i,j,k,nnew) = gls(i,j,k,nnew)+dt_lev
710 *cff_Hz*Pprod*gls(i,j,k,nstp) / std::max(tke(i,j,k,nstp),gls_Kmin);
711
712 // Compute dissipation of turbulent energy (m3/s3).
713 Real wall_fac = 1.0_rt;
714 if (Lmy25) {
715 // Parabolic wall function, L = ds db / (ds + db).
716 wall_fac=1.0_rt+gls_E2/(vonKar*vonKar)*
717 std::pow(std::pow(gls(i,j,k,nstp),( gls_exp1))*cmu_fac1*
718 std::pow(tke(i,j,k,nstp),-tke_exp1)*
719 (1.0_rt/ (z_w(i,j,k)-z_w(i,j,0))),2)+
720 0.25_rt/(vonKar*vonKar)*
721 std::pow(std::pow(gls(i,j,k,nstp), gls_exp1)*cmu_fac1*
722 std::pow(tke(i,j,k,nstp),-tke_exp1)*
723 (1.0_rt/ (z_w(i,j,N+1)-z_w(i,j,k))),2);
724 }
725 BCK(i,j,k)=cff_Hz*(1.0_rt+dt_lev*
726 std::pow(gls(i,j,k,nstp),(-gls_exp1))*cmu_fac2*
727 std::pow(tke(i,j,k,nstp), tke_exp2)+
728 dt_lev*(1.0_rt-cff1)*strat2*
729 (Akt(i,j,k,Temp_comp)-Akt_bak)/
730 tke(i,j,k,nstp))-
731 FCK(i,j,k)-FCK(i,j,k-1);
732 BCP(i,j,k)=cff_Hz*(1.0_rt+dt_lev*gls_c2*wall_fac*
733 std::pow(gls(i,j,k,nstp),-gls_exp1)*cmu_fac2*
734 std::pow(tke(i,j,k,nstp), tke_exp2)+
735 dt_lev*(1.0_rt-cff2)*gls_c3*strat2*
736 (Akt(i,j,k,Temp_comp)-Akt_bak)/
737 tke(i,j,k,nstp))-
738 FCP(i,j,k)-FCP(i,j,k-1);
739 });
740 // Compute production and dissipation terms.
741 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
742 {
743 Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
744 //----------------------------------------------------------------------
745 // Time-step dissipation and vertical diffusion terms implicitly.
746 //----------------------------------------------------------------------
747 //
748 // Set Dirichlet surface and bottom boundary conditions. Compute
749 // surface roughness from wind stress (Charnok) and set Craig and
750 // Banner wave breaking surface flux, if appropriate.
751
752
753 tke(i,j,N+1,nnew)=std::max(cmu_fac3*0.5_rt*
754 std::sqrt((sustr(i,j,0)+sustr(i+1,j,0))*(sustr(i,j,0)+sustr(i+1,j,0))+
755 (svstr(i,j,0)+svstr(i,j+1,0))*(svstr(i,j,0)+svstr(i,j+1,0))),
756 gls_Kmin);
757 tke(i,j,0,nnew)=std::max(cmu_fac3*0.5_rt*
758 std::sqrt((bustr(i,j,0)+bustr(i+1,j,0))*(bustr(i,j,0)+bustr(i+1,j,0))+
759 (bvstr(i,j,0)+bvstr(i,j+1,0))*(bvstr(i,j,0)+bvstr(i,j+1,0))),
760 gls_Kmin);
761
762 gls(i,j,N+1,nnew)=std::max(std::pow(gls_cmu0,gls_p)*
763 std::pow(tke(i,j,N+1,nnew),gls_m)*
764 std::pow(L_sft*Zos_eff,gls_n), gls_Pmin);
765 Real cff_gls = gls_fac4*std::pow(vonKar*Zob_min,gls_n);
766 gls(i,j,0,nnew)=std::max(cff_gls*std::pow(tke(i,j,0,nnew),(gls_m)), gls_Pmin);
767
768 // Solve tri-diagonal system for turbulent kinetic energy.
769 // Might be N instead of N-1?
770 Real tke_fluxt = 0.0_rt;
771 Real tke_fluxb = 0.0_rt;
772 Real cff_BCK = 1.0_rt/BCK(i,j,N);
773 CF(i,j,N)=cff_BCK*FCK(i,j,N-1);
774 tke(i,j,N,nnew)=cff_BCK*(tke(i,j,N,nnew)+tke_fluxt);
775 for (int k=N-1;k>=1;k--) {
776 cff_BCK = 1.0_rt / (BCK(i,j,k)-CF(i,j,k+1)*FCK(i,j,k));
777 CF(i,j,k) = cff_BCK * FCK(i,j,k-1);
778 tke(i,j,k,nnew) = cff_BCK * (tke(i,j,k,nnew) - FCK(i,j,k) * tke(i,j,k+1,nnew));
779 }
780 tke(i,j,1,nnew) = tke(i,j,1,nnew) - cff_BCK * tke_fluxb;
781 tke(i,j,1,nnew) = std::max(tke(i,j,1,nnew),gls_Kmin);
782 for (int k=2;k<=N;k++) {
783 tke(i,j,k,nnew) = tke(i,j,k,nnew) - CF(i,j,k) * tke(i,j,k-1,nnew);
784 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
785 }
786
787 // Solve tri-diagonal system for generic statistical field.
788 Real cff_tke = 0.5_rt * (tke(i,j,N+1,nnew) + tke(i,j,N,nnew));
789 Real gls_fluxt = dt_lev*gls_fac3*std::pow(cff_tke,gls_m)*
790 std::pow(L_sft,(gls_n))*
791 std::pow(Zos_eff+0.5_rt*Hz(i,j,N),gls_n-1.0_rt)*
792 0.5_rt*(Akp(i,j,N+1)+Akp(i,j,N));
793 cff_tke=0.5_rt*(tke(i,j,0,nnew)+tke(i,j,1,nnew));
794 Real gls_fluxb = dt_lev*gls_fac2*std::pow(cff_tke,gls_m)*
795 std::pow(0.5_rt*Hz(i,j,0)+Zob_min,gls_n-1.0_rt)*
796 0.5_rt*(Akp(i,j,0)+Akp(i,j,1));
797 Real cff_BCP = 1.0_rt / BCP(i,j,N);
798 CF(i,j,N) = cff_BCP * FCP(i,j,N-1);
799 gls(i,j,N,nnew)=cff_BCP*(gls(i,j,N,nnew)-gls_fluxt);
800 for (int k=N-1;k>=1;k--) {
801 cff_BCP = 1.0_rt / (BCP(i,j,k)-CF(i,j,k+1)*FCP(i,j,k));
802 CF(i,j,k) = cff_BCP * FCP(i,j,k-1);
803 gls(i,j,k,nnew) = cff_BCP * (gls(i,j,k,nnew) - FCP(i,j,k)*gls(i,j,k+1,nnew));
804 }
805 gls(i,j,1,nnew) = gls(i,j,1,nnew)-cff_BCP*gls_fluxb;
806 for (int k=2; k<=N; k++) {
807 gls(i,j,k,nnew) = gls(i,j,k,nnew) - CF(i,j,k) * gls(i,j,k-1,nnew);
808 }
809 });
810
811 // Compute vertical mixing coefficients (m2/s).
812 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
813 {
814 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
815 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
816 if (gls_n >= 0.0_rt) {
817 gls(i,j,k,nnew)=std::min(gls(i,j,k,nnew),gls_fac5*
818 std::pow(tke(i,j,k,nnew),tke_exp4)*
819 std::pow(std::sqrt(std::max(0.0_rt,
820 buoy2(i,j,k)))+eps,-gls_n));
821 } else {
822 gls(i,j,k,nnew)=std::max(gls(i,j,k,nnew), gls_fac5*
823 std::pow(tke(i,j,k,nnew),(tke_exp4))*
824 std::pow((std::sqrt(std::max(0.0_rt,
825 buoy2(i,j,k)))+eps),(-gls_n)));
826 }
827 Real Ls_lmt;
828 Real Ls_unlmt=std::max(eps,
829 std::pow(gls(i,j,k,nnew),( gls_exp1))*cmu_fac1*
830 std::pow(tke(i,j,k,nnew),(-tke_exp1)));
831 // Some problems are very sensitive to this condition (ultimate cause of
832 // some discrepancies in BoundaryLayer test between CPU and GPU)
833 if (buoy2(i,j,k) > 0.0_rt) {
834 Ls_lmt=std::min(Ls_unlmt,
835 std::sqrt(0.56_rt*tke(i,j,k,nnew)/
836 (std::max(0.0_rt,buoy2(i,j,k))+eps)));
837 } else {
838 Ls_lmt = Ls_unlmt;
839 }
840 //
841 // Recompute gls based on limited length scale
842 //
843 gls(i,j,k,nnew)=std::max(std::pow(gls_cmu0,gls_p)*
844 std::pow(tke(i,j,k,nnew),gls_m)*
845 std::pow(Ls_lmt,gls_n), gls_Pmin);
846
847 // Compute nondimensional stability functions for tracers (Sh) and
848 // momentum (Sm).
849 Real Sh, Sm;
850 Real Gh=std::min(gls_Gh0,-buoy2(i,j,k)*Ls_lmt*Ls_lmt/
851 (2.0_rt*tke(i,j,k,nnew)));
852 Gh=std::min(Gh,Gh-(Gh-gls_Ghcri)*(Gh-gls_Ghcri)/
853 (Gh+gls_Gh0-2.0_rt*gls_Ghcri));
854 Gh=std::max(Gh,gls_Ghmin);
855
856 if (gls_stability_type == GLS_StabilityType::Canuto_A ||
857 gls_stability_type == GLS_StabilityType::Canuto_B) {
858 //
859 // Canuto stability: Compute shear number.
860 //
861 Real Gm=(gls_b0/gls_fac6-gls_b1*Gh+gls_b3*gls_fac6*(Gh*Gh))/
862 (gls_b2-gls_b4*gls_fac6*Gh);
863 Gm=std::min(Gm,shear2(i,j,k)*Ls_lmt*Ls_lmt/
864 (2.0_rt*tke(i,j,k,nnew)));
865 /////Gm=std::min(Gm,(gls_s1*gls_fac6*Gh-gls_s0)/(gls_s2*gls_fac6));
866 //
867 // Compute stability functions
868 //
869 Real stab_cff=gls_b0-gls_b1*gls_fac6*Gh+gls_b2*gls_fac6*Gm+
870 gls_b3*gls_fac6*gls_fac6*Gh*Gh-gls_b4*gls_fac6*gls_fac6*Gh*Gm+
871 gls_b5*gls_fac6*gls_fac6*Gm*Gm;
872 Sm=(gls_s0-gls_s1*gls_fac6*Gh+gls_s2*gls_fac6*Gm)/stab_cff;
873 Sh=(gls_s4-gls_s5*gls_fac6*Gh+gls_s6*gls_fac6*Gm)/stab_cff;
874 Sm=std::max(Sm,0.0_rt);
875 Sh=std::max(Sh,0.0_rt);
876
877 //
878 // Relate Canuto stability to ROMS notation
879 //
880 Real gls_cmu0_cube = gls_cmu0 * gls_cmu0 * gls_cmu0;
881 Sm=Sm*sqrt2/(gls_cmu0_cube);
882 Sh=Sh*sqrt2/gls_cmu0_cube;
883 } else if (gls_stability_type == GLS_StabilityType::Galperin) {
884 Real cff_galperin = 1.0_rt - my_Sh2*Gh;
885 Sh = my_Sh1 / cff_galperin;
886 Sm = (my_Sm3+Sh*Gh*my_Sm4)/(1.0_rt-my_Sm2*Gh);
887 }
888
889 // Compute vertical mixing (m2/s) coefficients of momentum and
890 // tracers. Average ql over the two timesteps rather than using
891 // the new Lscale and just averaging tke.
892
893 Real ql=sqrt2*0.5_rt*(Ls_lmt*std::sqrt(tke(i,j,k,nnew))+
894 Lscale(i,j,k)*std::sqrt(tke(i,j,k,nstp)));
895 Akv(i,j,k)=Akv_bak+Sm*ql;
896 for (int n=0; n<NCONS; n++) {
897 Akt(i,j,k,n)=Akt_bak+Sh*ql;
898 }
899
900 // Compute vertical mixing (m2/s) coefficients of turbulent kinetic
901 // energy and generic statistical field.
902
903 Akk(i,j,k)=Akk_bak+Sm*ql/gls_sigk;
904 Akp(i,j,k)=Akp_bak+Sm*ql*ogls_sigp;
905
906 // Save limited length scale.
907 Lscale(i,j,k)=Ls_lmt;
908 });
909 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
910 {
911 Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
912 Akv(i,j,N+1)=Akv_bak+L_sft*Zos_eff*gls_cmu0*
913 std::sqrt(tke(i,j,N+1,nnew));
914 Akv(i,j,0)=Akv_bak+vonKar*Zob_min*gls_cmu0*
915 std::sqrt(tke(i,j,0,nnew));
916
917 Akk(i,j,N+1)=Akk_bak+Akv(i,j,N+1)/gls_sigk;
918 Akk(i,j,0)=Akk_bak+Akv(i,j,0)/gls_sigk;
919 Akp(i,j,N+1)=Akp_bak+Akv(i,j,N+1)*ogls_sigp;
920 Akp(i,j,0)=Akp_bak+Akv(i,j,0)/gls_sigp_cb;
921
922 for (int n=0; n<NCONS; n++) {
923 Akt(i,j,N+1,n) = Akt_bak;
924 Akt(i,j,0,n) = Akt_bak;
925 }
926 });
927 }
928
929 for (int icomp=0; icomp<3; icomp++) {
930 FillPatch(lev, t_old[lev], *mf_tke, GetVecOfPtrs(vec_tke), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
931 FillPatch(lev, t_old[lev], *mf_gls, GetVecOfPtrs(vec_gls), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
932 }
933 for (int icomp=0; icomp<NCONS; icomp++) {
934 FillPatch(lev, t_old[lev], *mf_Akt, GetVecOfPtrs(vec_Akt), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
935 }
936 FillPatchNoBC(lev, t_old[lev], *mf_Akv, GetVecOfPtrs(vec_Akv), BdyVars::null);
937 FillPatchNoBC(lev, t_old[lev], *mf_Akp, GetVecOfPtrs(vec_Akp), BdyVars::null);
938 FillPatchNoBC(lev, t_old[lev], *mf_Akk, GetVecOfPtrs(vec_Akk), BdyVars::null);
939}
vonKar
constexpr amrex::Real vonKar
Definition REMORA_Constants.H:17
GLS_StabilityType::Canuto_B
@ Canuto_B
GLS_StabilityType::Canuto_A
@ Canuto_A
GLS_StabilityType::Galperin
@ Galperin
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:13
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
NCONS
#define NCONS
Definition REMORA_IndexDefines.H:11
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1120
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:355
REMORA::vec_ZoBot
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ZoBot
Bottom roughness length [m], defined at rho points.
Definition REMORA.H:317
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:403
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:405
REMORA::vec_sustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sustr
Surface stress in the u direction.
Definition REMORA.H:282
REMORA::vec_Lscale
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Lscale
Vertical mixing turbulent length scale.
Definition REMORA.H:407
REMORA::vec_Hz
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hz
Width of cells in the vertical (z-) direction (3D, Hz in ROMS)
Definition REMORA.H:232
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:252
REMORA::FillPatchNoBC
void FillPatchNoBC(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=true)
Fill a new MultiFab by copying in phi from valid region and filling ghost cells without applying boun...
Definition REMORA_FillPatch.cpp:192
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:210
REMORA::gls_prestep
void gls_prestep(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int iic, const int ntfirst, const int N, const amrex::Real dt_lev)
Prestep for GLS calculation.
Definition REMORA_gls.cpp:20
REMORA::vec_bvf
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvf
Brunt-Vaisala frequency (3D)
Definition REMORA.H:390
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1108
REMORA::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=true, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=amrex::Real(0.0), const amrex::MultiFab &mf_calc=amrex::MultiFab())
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
Definition REMORA_FillPatch.cpp:30
REMORA::vec_svstr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_svstr
Surface stress in the v direction.
Definition REMORA.H:284
REMORA::vec_Huon
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Huon
u-volume flux (3D)
Definition REMORA.H:234
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:212
REMORA::t_new
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1098
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1221
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:409
REMORA::cons_old
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:208
REMORA::vec_bustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bustr
Bottom stress in the u direction.
Definition REMORA.H:320
REMORA::vec_bvstr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvstr
Bottom stress in the v direction.
Definition REMORA.H:322
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:357
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:250
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1100
REMORA::vec_z_w
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_w
z coordinates at w points (faces between z-cells)
Definition REMORA.H:264
REMORA::vec_Hvom
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hvom
v-volume flux (3D)
Definition REMORA.H:236
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:411
REMORA::gls_corrector
void gls_corrector(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_Akv, amrex::MultiFab *mf_Akt, amrex::MultiFab *mf_Akk, amrex::MultiFab *mf_Akp, amrex::MultiFab *mf_mskr, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int N, const amrex::Real dt_lev)
Corrector step for GLS calculation.
Definition REMORA_gls.cpp:246
SolverChoice::gls_stability_type
GLS_StabilityType gls_stability_type
Definition REMORA_DataStruct.H:572
SolverChoice::Akv_bak
amrex::Real Akv_bak
Definition REMORA_DataStruct.H:652
SolverChoice::gls_sigp
amrex::Real gls_sigp
Definition REMORA_DataStruct.H:647
SolverChoice::gls_m
amrex::Real gls_m
Definition REMORA_DataStruct.H:636
SolverChoice::my_B2
amrex::Real my_B2
Definition REMORA_DataStruct.H:674
SolverChoice::my_A1
amrex::Real my_A1
Definition REMORA_DataStruct.H:671
SolverChoice::gls_sigk
amrex::Real gls_sigk
Definition REMORA_DataStruct.H:646
SolverChoice::gls_L3
amrex::Real gls_L3
Definition REMORA_DataStruct.H:663
SolverChoice::gls_cmu0
amrex::Real gls_cmu0
Definition REMORA_DataStruct.H:641
SolverChoice::gls_L6
amrex::Real gls_L6
Definition REMORA_DataStruct.H:666
SolverChoice::gls_L2
amrex::Real gls_L2
Definition REMORA_DataStruct.H:662
SolverChoice::Akk_bak
amrex::Real Akk_bak
Definition REMORA_DataStruct.H:650
SolverChoice::gls_L1
amrex::Real gls_L1
Definition REMORA_DataStruct.H:661
SolverChoice::gls_n
amrex::Real gls_n
Definition REMORA_DataStruct.H:637
SolverChoice::Akt_bak
amrex::Real Akt_bak
Definition REMORA_DataStruct.H:653
SolverChoice::gls_c3m
amrex::Real gls_c3m
Definition REMORA_DataStruct.H:644
SolverChoice::my_A2
amrex::Real my_A2
Definition REMORA_DataStruct.H:672
SolverChoice::gls_Gh0
amrex::Real gls_Gh0
Definition REMORA_DataStruct.H:656
SolverChoice::gls_Ghmin
amrex::Real gls_Ghmin
Definition REMORA_DataStruct.H:658
SolverChoice::my_C1
amrex::Real my_C1
Definition REMORA_DataStruct.H:675
SolverChoice::gls_c1
amrex::Real gls_c1
Definition REMORA_DataStruct.H:642
SolverChoice::gls_L5
amrex::Real gls_L5
Definition REMORA_DataStruct.H:665
SolverChoice::gls_L8
amrex::Real gls_L8
Definition REMORA_DataStruct.H:668
SolverChoice::gls_Kmin
amrex::Real gls_Kmin
Definition REMORA_DataStruct.H:638
SolverChoice::gls_c3p
amrex::Real gls_c3p
Definition REMORA_DataStruct.H:645
SolverChoice::gls_p
amrex::Real gls_p
Definition REMORA_DataStruct.H:635
SolverChoice::gls_L4
amrex::Real gls_L4
Definition REMORA_DataStruct.H:664
SolverChoice::my_B1
amrex::Real my_B1
Definition REMORA_DataStruct.H:673
SolverChoice::gls_L7
amrex::Real gls_L7
Definition REMORA_DataStruct.H:667
SolverChoice::gls_Ghcri
amrex::Real gls_Ghcri
Definition REMORA_DataStruct.H:657
SolverChoice::gls_Pmin
amrex::Real gls_Pmin
Definition REMORA_DataStruct.H:639
SolverChoice::gls_c2
amrex::Real gls_c2
Definition REMORA_DataStruct.H:643
SolverChoice::Akp_bak
amrex::Real Akp_bak
Definition REMORA_DataStruct.H:651
SolverChoice::gls_E2
amrex::Real gls_E2
Definition REMORA_DataStruct.H:659