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