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,BCVars::xvel_bc,0,1,domain_bcs_type,bcrs_x);
60 amrex::setBC(ybx,domain,BCVars::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>(0.);
63 FArrayBox fab_FX(xbx_hi, 1, amrex::The_Async_Arena()); fab_FX.template setVal<RunOn::Device>(0.);
64 FArrayBox fab_FXL(xbx_hi, 1, amrex::The_Async_Arena()); fab_FXL.template setVal<RunOn::Device>(0.);
65 FArrayBox fab_EF(ybx_hi, 1, amrex::The_Async_Arena()); fab_EF.template setVal<RunOn::Device>(0.);
66 FArrayBox fab_FE(ybx_hi, 1, amrex::The_Async_Arena()); fab_FE.template setVal<RunOn::Device>(0.);
67 FArrayBox fab_FEL(ybx_hi, 1, amrex::The_Async_Arena()); fab_FEL.template setVal<RunOn::Device>(0.);
68 FArrayBox fab_Hz_half(bx, 1, amrex::The_Async_Arena()); fab_Hz_half.template setVal<RunOn::Device>(0.);
69 FArrayBox fab_CF(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_CF.template setVal<RunOn::Device>(0.);
70 FArrayBox fab_FC(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FC.template setVal<RunOn::Device>(0.);
71 FArrayBox fab_FCL(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FCL.template setVal<RunOn::Device>(0.);
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 = 1.0_rt/6.0_rt;
104 XF(i,j,k) = 0.5_rt * (Huon(i,j,k) + Huon(i,j,k-1));
105 FX(i,j,k) = XF(i,j,k) * 0.5_rt * (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) * 0.5_rt * (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 = 1.0_rt/6.0_rt;
131 EF(i,j,k) = 0.5_rt * (Hvom(i,j,k) + Hvom(i,j,k-1));
132 FE(i,j,k) = EF(i,j,k) * 0.5_rt * (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) * 0.5_rt * (gls(i,j-1,k,nstp) + gls(i,j,k,nstp) -
135 cff * (gradL_jp1 - gradL_jm1));
136 });
137
138 Real gamma = 1.0_rt / 6.0_rt;
139 Real cff1, cff2, cff3;
140 int indx;
141 // Time step horizontal advection
142 if (iic == ntfirst) {
143 cff1 = 1.0_rt;
144 cff2 = 0.0_rt;
145 cff3 = 0.5_rt * dt_lev;
146 indx = nstp;
147 } else {
148 cff1 = 0.5_rt + gamma;
149 cff2 = 0.5_rt - gamma;
150 cff3 = (1.0_rt - 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 = 0.5_rt * (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) = 0.5_rt * (W(i,j,k+1) + W(i,j,k));
176 if (k == 0) {
177 Real cff1_vadv = 1.0_rt / 3.0_rt;
178 Real cff2_vadv = 5.0_rt / 6.0_rt;
179 Real cff3_vadv = 1.0_rt / 6.0_rt;
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 = 1.0_rt / 3.0_rt;
188 Real cff2_vadv = 5.0_rt / 6.0_rt;
189 Real cff3_vadv = 1.0_rt / 6.0_rt;
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 = 7.0_rt / 12.0_rt;
198 Real cff2_vadv = 1.0_rt / 12.0_rt;
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 = 0.5_rt * dt_lev;
209 } else {
210 cff3 = (1.0_rt - 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 = 1.0_rt / 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), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
225 FillPatch(lev, t_old[lev], *vec_gls[lev], GetVecOfPtrs(vec_gls), BCVars::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 == 0.0) &&
260 (solverChoice.gls_n == 1.0) &&
261 (solverChoice.gls_m == 1.0)) ? true : false;
262
263 Real L_sft = vonKar;
264 Real gls_sigp_cb = solverChoice.gls_sigp;
265 Real ogls_sigp = 1.0_rt/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 Akt_bak = solverChoice.Akt_bak;
281 Real Akp_bak = solverChoice.Akp_bak;
282 Real Akk_bak = solverChoice.Akk_bak;
283
284 Real gls_c1 = solverChoice.gls_c1;
285 Real gls_c2 = solverChoice.gls_c2;
286 Real gls_E2 = solverChoice.gls_E2;
287 Real gls_sigk = solverChoice.gls_sigk;
288 auto gls_stability_type = solverChoice.gls_stability_type;
289
290 Real sqrt2 = std::sqrt(2.0_rt);
291 Real cmu_fac1 = std::pow(solverChoice.gls_cmu0,(-solverChoice.gls_p/solverChoice.gls_n));
292 Real cmu_fac2 = std::pow(solverChoice.gls_cmu0,(3.0_rt+solverChoice.gls_p/solverChoice.gls_n));
293 Real cmu_fac3 = 1.0_rt/std::pow(solverChoice.gls_cmu0,2.0_rt);
294
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_exp4 = solverChoice.gls_m+0.5_rt*solverChoice.gls_n;
305
306 Real cmu0_exp_p = std::pow(gls_cmu0, gls_p);
307 Real gls_cmu0_cube = gls_cmu0 * gls_cmu0 * gls_cmu0;
308
309 Real gls_s0, gls_s1, gls_s2, gls_s4, gls_s5, gls_s6;
310 Real gls_b0, gls_b1, gls_b2, gls_b3, gls_b4, gls_b5;
311 Real my_Sm2, my_Sh1, my_Sh2, my_Sm3, my_Sm4;
312
313 // Compute parameters for Canuto et al. (2001) stability functions.
314 // (Canuto, V.M., Cheng, H.Y., and Dubovikov, M.S., 2001: Ocean
315 // turbulence. Part I: One-point closure model - momentum and
316 // heat vertical diffusivities, JPO, 1413-1426).
317
318 if (solverChoice.gls_stability_type == GLS_StabilityType::Canuto_A ||
319 solverChoice.gls_stability_type == GLS_StabilityType::Canuto_B) {
320
321 gls_s0=3.0_rt/2.0_rt*solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L5;
322 gls_s1=-solverChoice.gls_L4*(solverChoice.gls_L6+solverChoice.gls_L7)
323 +2.0_rt*solverChoice.gls_L4*solverChoice.gls_L5*
324 (solverChoice.gls_L1-1.0_rt/3.0_rt*solverChoice.gls_L2-solverChoice.gls_L3)
325 +3.0_rt/2.0_rt*
326 solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L8;
327 gls_s2=-3.0_rt/8.0_rt*solverChoice.gls_L1
328 *(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
329 gls_s4=2.0_rt*solverChoice.gls_L5;
330 gls_s5=2.0_rt*solverChoice.gls_L4;
331 gls_s6=2.0_rt/3.0_rt*solverChoice.gls_L5
332 *(3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
333 1.0_rt/2.0_rt*solverChoice.gls_L5*solverChoice.gls_L1*(3.0_rt*solverChoice.gls_L3-solverChoice.gls_L2)+
334 3.0_rt/4.0_rt*solverChoice.gls_L1*(solverChoice.gls_L6-solverChoice.gls_L7);
335 gls_b0=3.0_rt*solverChoice.gls_L5*solverChoice.gls_L5;
336 gls_b1=solverChoice.gls_L5*(7.0_rt*solverChoice.gls_L4+3.0_rt*solverChoice.gls_L8);
337 gls_b2=solverChoice.gls_L5*solverChoice.gls_L5*(3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
338 3.0_rt/4.0_rt*(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
339 gls_b3=solverChoice.gls_L4*(4.0_rt*solverChoice.gls_L4+3.0_rt*solverChoice.gls_L8);
340 gls_b5=1.0_rt/4.0_rt*(solverChoice.gls_L2*solverChoice.gls_L2-3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3)*
341 (solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
342 gls_b4=solverChoice.gls_L4*(solverChoice.gls_L2*solverChoice.gls_L6-3.0_rt*solverChoice.gls_L3*solverChoice.gls_L7-
343 solverChoice.gls_L5*(solverChoice.gls_L2*solverChoice.gls_L2-solverChoice.gls_L3*solverChoice.gls_L3))+solverChoice.gls_L5*solverChoice.gls_L8*
344 (3.0_rt*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2);
345 my_Sm2 = 0.0_rt;
346 my_Sm3 = 0.0_rt;
347 my_Sm4 = 0.0_rt;
348 my_Sh1 = 0.0_rt;
349 my_Sh2 = 0.0_rt;
350 } else {
351 gls_s0 = 0.0_rt;
352 gls_s1 = 0.0_rt;
353 gls_s2 = 0.0_rt;
354 gls_s4 = 0.0_rt;
355 gls_s5 = 0.0_rt;
356 gls_s6 = 0.0_rt;
357 gls_b0 = 0.0_rt;
358 gls_b1 = 0.0_rt;
359 gls_b2 = 0.0_rt;
360 gls_b3 = 0.0_rt;
361 gls_b4 = 0.0_rt;
362 gls_b5 = 0.0_rt;
363 my_Sm2=9.0_rt*solverChoice.my_A1*solverChoice.my_A2;
364 my_Sm3=solverChoice.my_A1*(1.0_rt-3.0_rt*solverChoice.my_C1-6.0_rt*solverChoice.my_A1/solverChoice.my_B1);
365 my_Sm4=18.0_rt*solverChoice.my_A1*solverChoice.my_A1+9.0_rt*solverChoice.my_A1*solverChoice.my_A2;
366 my_Sh1=solverChoice.my_A2*(1.0_rt-6.0_rt*solverChoice.my_A1/solverChoice.my_B1);
367 my_Sh2=3.0_rt*solverChoice.my_A2*(6.0_rt*solverChoice.my_A1+solverChoice.my_B2);
368 }
369
370 Real Zos_min = std::max(solverChoice.Zos, 0.0001_rt);
371 Real Zos_eff = Zos_min;
372 Real Gadv = 1.0_rt/3.0_rt;
373 Real eps = 1.0e-10_rt;
374
375 const BoxArray& ba = cons_old[lev]->boxArray();
376 const DistributionMapping& dm = cons_old[lev]->DistributionMap();
377
378 int ncomp_w = 0;
379 int dU_comp = ncomp_w++;
380 int dV_comp = ncomp_w++;
381 int CF_comp = ncomp_w++;
382
383 int ncomp = 0;
384 int shear2_comp = ncomp++;
385 int shear2_cache_comp = ncomp++;
386 int buoy2_comp = ncomp++;
387
388 MultiFab mf_w(convert(ba, IntVect(0,0,1)),dm,ncomp_w,IntVect(NGROW,NGROW,0));
389 MultiFab mf(ba,dm,ncomp,IntVect(NGROW,NGROW,0));
390
391 const Box& domain = geom[0].Domain();
392 const auto dlo = amrex::lbound(domain);
393 const auto dhi = amrex::ubound(domain);
394
395 GeometryData const& geomdata = geom[0].data();
396 bool is_periodic_in_x = geomdata.isPeriodic(0);
397 bool is_periodic_in_y = geomdata.isPeriodic(1);
398
399 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
400 {
401 Box bx = mfi.tilebox();
402 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
403
404 Box bxD = bx;
405 bxD.makeSlab(2,0);
406 Box gbx1D = gbx1;
407 gbx1D.makeSlab(2,0);
408
409 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
410 Array4<Real> const& u = xvel_old[lev]->array(mfi);
411 Array4<Real> const& v = yvel_old[lev]->array(mfi);
412
413 auto dU = mf_w.array(mfi,dU_comp);
414 auto dV = mf_w.array(mfi,dV_comp);
415 auto CF = mf_w.array(mfi,CF_comp);
416 auto shear2_cached = mf.array(mfi,shear2_cache_comp);
417
418 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
419 {
420 CF(i,j,0) = 0.0_rt;
421 dU(i,j,0) = 0.0_rt;
422 dV(i,j,0) = 0.0_rt;
423 for (int k=1; k<=N; k++) {
424 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)));
425 CF(i,j,k) = cff * Hz(i,j,k);
426 dU(i,j,k)=cff*(3.0_rt*(u(i ,j,k)-u(i, j,k-1)+
427 u(i+1,j,k)-u(i+1,j,k-1))-Hz(i,j,k-1)*dU(i,j,k-1));
428 dV(i,j,k)=cff*(3.0_rt*(v(i,j ,k)-v(i,j ,k-1)+
429 v(i,j+1,k)-v(i,j+1,k-1))-Hz(i,j,k-1)*dV(i,j,k-1));
430 }
431 dU(i,j,N+1) = 0.0_rt;
432 dV(i,j,N+1) = 0.0_rt;
433 for (int k=N; k>=1; k--) {
434 dU(i,j,k) = dU(i,j,k) - CF(i,j,k) * dU(i,j,k+1);
435 dV(i,j,k) = dV(i,j,k) - CF(i,j,k) * dV(i,j,k+1);
436 }
437 shear2_cached(i,j,0) = 0.0_rt;
438 for (int k=1; k<=N; k++) {
439 shear2_cached(i,j,k) = dU(i,j,k) * dU(i,j,k) + dV(i,j,k) * dV(i,j,k);
440 }
441 });
442 }
443
444 // While potentially counterintuitive, this is what ROMS does for handling shear2 at all boundaries, even
445 // periodic
446 (*physbcs[lev])(mf,*mf_mskr,shear2_cache_comp,1,mf.nGrowVect(),t_new[lev],BCVars::foextrap_bc);
447 mf.setVal(0.0_rt,CF_comp,1);
448
449 int ncomp_fab = 0;
450 int tmp_buoy_comp = ncomp_fab++;
451 int tmp_shear_comp = ncomp_fab++;
452 int curvK_comp = ncomp_fab++;
453 int curvP_comp = ncomp_fab++;
454 int FXK_comp = ncomp_fab++;
455 int FXP_comp = ncomp_fab++;
456 int FEK_comp = ncomp_fab++;
457 int FEP_comp = ncomp_fab++;
458 int FCK_comp = ncomp_fab++;
459 int FCP_comp = ncomp_fab++;
460 int BCK_comp = ncomp_fab++;
461 int BCP_comp = ncomp_fab++;
462
463 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
464 {
465 Box bx = mfi.tilebox();
466 Box xbx = surroundingNodes(bx,0);
467 Box ybx = surroundingNodes(bx,1);
468 Box gbx1 = grow(bx,IntVect(NGROW-1,NGROW-1,0));
469
470 Box bx_rho = bx;
471 bx_rho.convert(IntVect(0,0,0));
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 ncompbc = 1;
480 Vector<BCRec> bcrs_x(ncompbc);
481 Vector<BCRec> bcrs_y(ncompbc);
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(gbx1,ncomp_fab, amrex::The_Async_Arena()); fab.template setVal<RunOn::Device>(0.);
508
509 auto CF = mf_w.array(mfi,CF_comp);
510 auto shear2 = mf.array(mfi,shear2_comp);
511 auto shear2_cached = mf.array(mfi,shear2_cache_comp);
512 auto buoy2 = mf.array(mfi,buoy2_comp);
513 Array4<Real> const& bvf = vec_bvf[lev]->array(mfi);
514
515 auto tmp_buoy = fab.array(tmp_buoy_comp);
516 auto tmp_shear = fab.array(tmp_shear_comp);
517 auto curvK = fab.array(curvK_comp);
518 auto curvP = fab.array(curvP_comp);
519 auto FXK = fab.array(FXK_comp);
520 auto FXP = fab.array(FXP_comp);
521 auto FEK = fab.array(FEK_comp);
522 auto FEP = fab.array(FEP_comp);
523 auto FCK = fab.array(FCK_comp);
524 auto FCP = fab.array(FCP_comp);
525 auto BCK = fab.array(BCK_comp);
526 auto BCP = fab.array(BCP_comp);
527
528 auto Akt = mf_Akt->array(mfi);
529 auto Akv = mf_Akv->array(mfi);
530 auto Akp = mf_Akp->array(mfi);
531 auto Akk = mf_Akk->array(mfi);
532
533 ParallelFor(bx_growloxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
534 {
535 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));
536 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));
537 });
538
539 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
540 {
541 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));
542 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));
543 });
544
545 //Time step advective terms
546 ParallelFor(growLo(grow(xbx,IntVect(0,0,-1)),0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
547 {
548 Real gradK, gradK_ip1, gradP, gradP_ip1;
549
550 if (i == dlo.x-1 && !is_periodic_in_x) {
551 gradK_ip1 = tke(i+1,j,k,2)-tke(i ,j,k,2);
552 gradK = gradK_ip1;
553 gradP_ip1 = gls(i+1,j,k,2)-gls(i ,j,k,2);
554 gradP = gradP_ip1;
555 } else if (i == dhi.x+1 && !is_periodic_in_x) {
556 gradK = tke(i ,j,k,2)-tke(i-1,j,k,2);
557 gradK_ip1 = gradK;
558 gradP = gls(i ,j,k,2)-gls(i-1,j,k,2);
559 gradP_ip1 = gradP;
560 } else {
561 gradK = (tke(i ,j,k,2)-tke(i-1,j,k,2)) * msku(i ,j,0);
562 gradK_ip1 = (tke(i+1,j,k,2)-tke(i ,j,k,2)) * msku(i+1,j,0);
563 gradP = (gls(i ,j,k,2)-gls(i-1,j,k,2)) * msku(i ,j,0);
564 gradP_ip1 = (gls(i+1,j,k,2)-gls(i ,j,k,2)) * msku(i+1,j,0);
565 }
566
567 curvK(i,j,k) = gradK_ip1 - gradK;
568 curvP(i,j,k) = gradP_ip1 - gradP;
569 });
570 ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
571 {
572 Real cff = 0.5_rt * (Huon(i,j,k) + Huon(i,j,k-1));
573 Real cff1 = (cff > 0.0) ? curvK(i-1,j,k) : curvK(i,j,k);
574 Real cff2 = (cff > 0.0) ? curvP(i-1,j,k) : curvP(i,j,k);
575
576 FXK(i,j,k) = cff * 0.5_rt * (tke(i-1,j,k,2)+tke(i,j,k,2)-Gadv*cff1);
577 FXP(i,j,k) = cff * 0.5_rt * (gls(i-1,j,k,2)+gls(i,j,k,2)-Gadv*cff2);
578 });
579
580 //Time step advective terms
581 ParallelFor(growLo(grow(ybx,IntVect(0,0,-1)),1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
582 {
583 Real gradK = (tke(i,j ,k,2)-tke(i,j-1,k,2)) * mskv(i,j ,0);
584 Real gradK_jp1 = (tke(i,j+1,k,2)-tke(i,j ,k,2)) * mskv(i,j+1,0);
585 Real gradP = (gls(i,j ,k,2)-gls(i,j-1,k,2)) * mskv(i,j ,0);
586 Real gradP_jp1 = (gls(i,j+1,k,2)-gls(i,j ,k,2)) * mskv(i,j+1,0);
587
588 if (j == dlo.y-1 && !is_periodic_in_y) {
589 gradK = gradK_jp1;
590 gradP = gradP_jp1;
591 }
592 else if (j == dhi.y+1 && !is_periodic_in_y) {
593 gradK_jp1 = gradK;
594 gradP_jp1 = gradP;
595 }
596
597 curvK(i,j,k) = gradK_jp1 - gradK;
598 curvP(i,j,k) = gradP_jp1 - gradP;
599 });
600 ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
601 {
602 Real cff = 0.5_rt * (Hvom(i,j,k) + Hvom(i,j,k-1));
603 Real cff1 = (cff > 0.0) ? curvK(i,j-1,k) : curvK(i,j,k);
604 Real cff2 = (cff > 0.0) ? curvP(i,j-1,k) : curvP(i,j,k);
605
606 FEK(i,j,k) = cff * 0.5_rt * (tke(i,j-1,k,2)+tke(i,j,k,2)-Gadv*cff1);
607 FEP(i,j,k) = cff * 0.5_rt * (gls(i,j-1,k,2)+gls(i,j,k,2)-Gadv*cff2);
608 });
609
610 Real gls_Kmin = solverChoice.gls_Kmin;
611 Real gls_Pmin = solverChoice.gls_Pmin;
612 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
613 {
614 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
615 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff * (FXK(i+1,j ,k)-FXK(i,j,k)+
616 FEK(i ,j+1,k)-FEK(i,j,k));
617 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
618
619 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff * (FXP(i+1,j ,k)-FXP(i,j,k)+
620 FEP(i ,j+1,k)-FEP(i,j,k));
621 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew), gls_Pmin);
622 });
623
624 // Vertical advection
625 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
626 {
627 Real cff1 = 7.0_rt / 12.0_rt;
628 Real cff2 = 1.0_rt / 12.0_rt;
629 for (int k=1; k<=N-1; k++) {
630 Real cff = 0.5_rt * (W(i,j,k+1)+W(i,j,k));
631 FCK(i,j,k) = cff * (cff1 * (tke(i,j,k ,2)+tke(i,j,k+1,2))-
632 cff2 * (tke(i,j,k-1,2)+tke(i,j,k+2,2)));
633 FCP(i,j,k) = cff * (cff1 * (gls(i,j,k ,2)+gls(i,j,k+1,2))-
634 cff2 * (gls(i,j,k-1,2)+gls(i,j,k+2,2)));
635 }
636 cff1 = 1.0_rt/3.0_rt;
637 cff2 = 5.0_rt/6.0_rt;
638 Real cff3 = 1.0_rt / 6.0_rt;
639 Real cff = 0.5_rt * (W(i,j,0)+W(i,j,1));
640 FCK(i,j,0) = cff * (cff1 * tke(i,j,0,2)+cff2 * tke(i,j,1,2)-cff3 * tke(i,j,2,2));
641 FCP(i,j,0) = cff * (cff1 * gls(i,j,0,2)+cff2 * gls(i,j,1,2)-cff3 * gls(i,j,2,2));
642 cff = 0.5_rt * (W(i,j,N+1)+W(i,j,N));
643 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));
644 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));
645 });
646 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
647 {
648 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
649 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff*(FCK(i,j,k )-FCK(i,j,k-1));
650 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
651 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff*(FCP(i,j,k )-FCP(i,j,k-1));
652 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
653 });
654
655 // Compute vertical mixing, turbulent production and turbulent
656 // dissipation.
657 //
658 Real cff = -0.5 * dt_lev;
659 ParallelFor(convert(bx,IntVect(0,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
660 {
661 if (k==0 or k==N) {
662 FCK(i,j,k) = 0.0_rt;
663 FCP(i,j,k) = 0.0_rt;
664 } else {
665 FCK(i,j,k) = cff * (Akk(i,j,k) + Akk(i,j,k+1)) / Hz(i,j,k);
666 FCP(i,j,k) = cff * (Akp(i,j,k) + Akp(i,j,k+1)) / Hz(i,j,k);
667 }
668 });
669 // Compute production and dissipation terms.
670 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
671 {
672 // Compute shear and buoyant production of turbulent energy (m3/s3)
673 // at W-points (ignore small negative values of buoyancy).
674 Real strat2 = buoy2(i,j,k);
675 Real gls_c3 = (strat2 > 0.0) ? gls_c3m : gls_c3p;
676 Real Kprod = shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
677 strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak);
678 Real Pprod = gls_c1 * shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
679 gls_c3 * strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak);
680
681 // If negative production terms, then add buoyancy to dissipation terms
682 // (BCK and BCP) below, using "cff1" and "cff2" as the on/off switch.
683 Real cff1 = (Kprod < 0.0_rt) ? 0.0_rt : 1.0_rt;
684 Real cff2 = (Pprod < 0.0_rt) ? 0.0_rt : 1.0_rt;
685 Kprod = (Kprod < 0.0_rt) ? Kprod + strat2*(Akt(i,j,k,Temp_comp)-Akt_bak) : Kprod;
686 Pprod = (Pprod < 0.0_rt) ? Pprod + gls_c3*strat2*(Akt(i,j,k,Temp_comp)-Akt_bak) : Pprod;
687 // Time-step shear and buoyancy production terms.
688 Real cff_Hz = 0.5_rt * (Hz(i,j,k) + Hz(i,j,k-1));
689 tke(i,j,k,nnew) = tke(i,j,k,nnew)+dt_lev * cff_Hz * Kprod;
690 gls(i,j,k,nnew) = gls(i,j,k,nnew)+dt_lev
691 *cff_Hz*Pprod*gls(i,j,k,nstp) / std::max(tke(i,j,k,nstp),gls_Kmin);
692
693 Real gls_exp_exp1 = std::pow(gls(i,j,k,nstp),gls_exp1);
694 Real gls_exp_mexp1 = 1.0_rt / (gls_exp_exp1);
695 Real tke_exp_mexp1 = std::pow(tke(i,j,k,nstp),-tke_exp1);
696 Real tke_exp_exp2 = std::pow(tke(i,j,k,nstp),tke_exp2);
697
698 // Compute dissipation of turbulent energy (m3/s3).
699 Real wall_fac = 1.0_rt;
700 if (Lmy25) {
701 wall_fac=1.0_rt+gls_E2/(vonKar*vonKar)*
702 std::pow(gls_exp_exp1*cmu_fac1*
703 tke_exp_mexp1*
704 (1.0_rt/ (z_w(i,j,k)-z_w(i,j,0))),2)+
705 0.25_rt/(vonKar*vonKar)*
706 std::pow(gls_exp_exp1*cmu_fac1*
707 tke_exp_mexp1*
708 (1.0_rt/ (z_w(i,j,N+1)-z_w(i,j,k))),2);
709 }
710 BCK(i,j,k)=cff_Hz*(1.0_rt+dt_lev*
711 gls_exp_mexp1*cmu_fac2*
712 tke_exp_exp2+
713 dt_lev*(1.0_rt-cff1)*strat2*
714 (Akt(i,j,k,Temp_comp)-Akt_bak)/
715 tke(i,j,k,nstp))-
716 FCK(i,j,k)-FCK(i,j,k-1);
717 BCP(i,j,k)=cff_Hz*(1.0_rt+dt_lev*gls_c2*wall_fac*
718 gls_exp_mexp1*cmu_fac2*
719 tke_exp_exp2+
720 dt_lev*(1.0_rt-cff2)*gls_c3*strat2*
721 (Akt(i,j,k,Temp_comp)-Akt_bak)/
722 tke(i,j,k,nstp))-
723 FCP(i,j,k)-FCP(i,j,k-1);
724 });
725
726 // Compute production and dissipation terms.
727 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
728 {
729 Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
730 //----------------------------------------------------------------------
731 // Time-step dissipation and vertical diffusion terms implicitly.
732 //----------------------------------------------------------------------
733 //
734 // Set Dirichlet surface and bottom boundary conditions. Compute
735 // surface roughness from wind stress (Charnok) and set Craig and
736 // Banner wave breaking surface flux, if appropriate.
737
738
739 tke(i,j,N+1,nnew)=std::max(cmu_fac3*0.5_rt*
740 std::sqrt((sustr(i,j,0)+sustr(i+1,j,0))*(sustr(i,j,0)+sustr(i+1,j,0))+
741 (svstr(i,j,0)+svstr(i,j+1,0))*(svstr(i,j,0)+svstr(i,j+1,0))),
742 gls_Kmin);
743 tke(i,j,0,nnew)=std::max(cmu_fac3*0.5_rt*
744 std::sqrt((bustr(i,j,0)+bustr(i+1,j,0))*(bustr(i,j,0)+bustr(i+1,j,0))+
745 (bvstr(i,j,0)+bvstr(i,j+1,0))*(bvstr(i,j,0)+bvstr(i,j+1,0))),
746 gls_Kmin);
747
748 gls(i,j,N+1,nnew)=std::max(cmu0_exp_p*
749 std::pow(tke(i,j,N+1,nnew),gls_m)*
750 std::pow(L_sft*Zos_eff,gls_n), gls_Pmin);
751 Real cff_gls = gls_fac4*std::pow(vonKar*Zob_min,gls_n);
752 gls(i,j,0,nnew)=std::max(cff_gls*std::pow(tke(i,j,0,nnew),(gls_m)), gls_Pmin);
753
754 // Solve tri-diagonal system for turbulent kinetic energy.
755 // Might be N instead of N-1?
756 Real tke_fluxt = 0.0_rt;
757 Real tke_fluxb = 0.0_rt;
758 Real cff_BCK = 1.0_rt/BCK(i,j,N);
759 CF(i,j,N)=cff_BCK*FCK(i,j,N-1);
760 tke(i,j,N,nnew)=cff_BCK*(tke(i,j,N,nnew)+tke_fluxt);
761 for (int k=N-1;k>=1;k--) {
762 cff_BCK = 1.0_rt / (BCK(i,j,k)-CF(i,j,k+1)*FCK(i,j,k));
763 CF(i,j,k) = cff_BCK * FCK(i,j,k-1);
764 tke(i,j,k,nnew) = cff_BCK * (tke(i,j,k,nnew) - FCK(i,j,k) * tke(i,j,k+1,nnew));
765 }
766 tke(i,j,1,nnew) = tke(i,j,1,nnew) - cff_BCK * tke_fluxb;
767 tke(i,j,1,nnew) = std::max(tke(i,j,1,nnew),gls_Kmin);
768 for (int k=2;k<=N;k++) {
769 tke(i,j,k,nnew) = tke(i,j,k,nnew) - CF(i,j,k) * tke(i,j,k-1,nnew);
770 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
771 }
772
773 // Solve tri-diagonal system for generic statistical field.
774 Real cff_tke = 0.5_rt * (tke(i,j,N+1,nnew) + tke(i,j,N,nnew));
775 Real gls_fluxt = dt_lev*gls_fac3*std::pow(cff_tke,gls_m)*
776 std::pow(L_sft,(gls_n))*
777 std::pow(Zos_eff+0.5_rt*Hz(i,j,N),gls_n-1.0_rt)*
778 0.5_rt*(Akp(i,j,N+1)+Akp(i,j,N));
779 cff_tke=0.5_rt*(tke(i,j,0,nnew)+tke(i,j,1,nnew));
780 Real gls_fluxb = dt_lev*gls_fac2*std::pow(cff_tke,gls_m)*
781 std::pow(0.5_rt*Hz(i,j,0)+Zob_min,gls_n-1.0_rt)*
782 0.5_rt*(Akp(i,j,0)+Akp(i,j,1));
783 Real cff_BCP = 1.0_rt / BCP(i,j,N);
784 CF(i,j,N) = cff_BCP * FCP(i,j,N-1);
785 gls(i,j,N,nnew)=cff_BCP*(gls(i,j,N,nnew)-gls_fluxt);
786 for (int k=N-1;k>=1;k--) {
787 cff_BCP = 1.0_rt / (BCP(i,j,k)-CF(i,j,k+1)*FCP(i,j,k));
788 CF(i,j,k) = cff_BCP * FCP(i,j,k-1);
789 gls(i,j,k,nnew) = cff_BCP * (gls(i,j,k,nnew) - FCP(i,j,k)*gls(i,j,k+1,nnew));
790 }
791 gls(i,j,1,nnew) = gls(i,j,1,nnew)-cff_BCP*gls_fluxb;
792 for (int k=2; k<=N; k++) {
793 gls(i,j,k,nnew) = gls(i,j,k,nnew) - CF(i,j,k) * gls(i,j,k-1,nnew);
794 }
795 });
796
797 // Compute vertical mixing coefficients (m2/s).
798 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
799 {
800 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
801 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
802 Real gls_comparison = gls_fac5 *
803 std::pow(tke(i,j,k,nnew),tke_exp4)*
804 std::pow(std::sqrt(std::max(0.0_rt,
805 buoy2(i,j,k)))+eps,-gls_n);
806 gls(i,j,k,nnew) = (gls_n >= 0.0_rt) ? std::min(gls(i,j,k,nnew),gls_comparison) : std::max(gls(i,j,k,nnew),gls_comparison);
807 Real Ls_lmt;
808 Real Ls_unlmt=std::max(eps,
809 std::pow(gls(i,j,k,nnew),( gls_exp1))*cmu_fac1*
810 std::pow(tke(i,j,k,nnew),(-tke_exp1)));
811 // Some problems are very sensitive to this condition (ultimate cause of
812 // some discrepancies in BoundaryLayer test between CPU and GPU)
813 Ls_lmt = (buoy2(i,j,k) > 0.0_rt) ? std::min(Ls_unlmt,
814 std::sqrt(0.56_rt*tke(i,j,k,nnew)/
815 (std::max(0.0_rt,buoy2(i,j,k))+eps))) : Ls_unlmt;
816 //
817 // Recompute gls based on limited length scale
818 //
819 gls(i,j,k,nnew)=std::max(cmu0_exp_p*
820 std::pow(tke(i,j,k,nnew),gls_m)*
821 std::pow(Ls_lmt,gls_n), gls_Pmin);
822
823 // Compute nondimensional stability functions for tracers (Sh) and
824 // momentum (Sm).
825 Real Sh, Sm;
826 Real Gh=std::min(gls_Gh0,-buoy2(i,j,k)*Ls_lmt*Ls_lmt/
827 (2.0_rt*tke(i,j,k,nnew)));
828 Gh=std::min(Gh,Gh-(Gh-gls_Ghcri)*(Gh-gls_Ghcri)/
829 (Gh+gls_Gh0-2.0_rt*gls_Ghcri));
830 Gh=std::max(Gh,gls_Ghmin);
831
832 if (gls_stability_type == GLS_StabilityType::Canuto_A ||
833 gls_stability_type == GLS_StabilityType::Canuto_B) {
834 //
835 // Canuto stability: Compute shear number.
836 //
837 Real Gm=(gls_b0/gls_fac6-gls_b1*Gh+gls_b3*gls_fac6*(Gh*Gh))/
838 (gls_b2-gls_b4*gls_fac6*Gh);
839 Gm=std::min(Gm,shear2(i,j,k)*Ls_lmt*Ls_lmt/
840 (2.0_rt*tke(i,j,k,nnew)));
841 /////Gm=std::min(Gm,(gls_s1*gls_fac6*Gh-gls_s0)/(gls_s2*gls_fac6));
842 //
843 // Compute stability functions
844 //
845 Real stab_cff=gls_b0-gls_b1*gls_fac6*Gh+gls_b2*gls_fac6*Gm+
846 gls_b3*gls_fac6*gls_fac6*Gh*Gh-gls_b4*gls_fac6*gls_fac6*Gh*Gm+
847 gls_b5*gls_fac6*gls_fac6*Gm*Gm;
848 Sm=(gls_s0-gls_s1*gls_fac6*Gh+gls_s2*gls_fac6*Gm)/stab_cff;
849 Sh=(gls_s4-gls_s5*gls_fac6*Gh+gls_s6*gls_fac6*Gm)/stab_cff;
850 Sm=std::max(Sm,0.0_rt);
851 Sh=std::max(Sh,0.0_rt);
852
853 //
854 // Relate Canuto stability to ROMS notation
855 //
856 Sm=Sm*sqrt2/(gls_cmu0_cube);
857 Sh=Sh*sqrt2/gls_cmu0_cube;
858 } else if (gls_stability_type == GLS_StabilityType::Galperin) {
859 Real cff_galperin = 1.0_rt - my_Sh2*Gh;
860 Sh = my_Sh1 / cff_galperin;
861 Sm = (my_Sm3+Sh*Gh*my_Sm4)/(1.0_rt-my_Sm2*Gh);
862 }
863
864 // Compute vertical mixing (m2/s) coefficients of momentum and
865 // tracers. Average ql over the two timesteps rather than using
866 // the new Lscale and just averaging tke.
867
868 Real ql=sqrt2*0.5_rt*(Ls_lmt*std::sqrt(tke(i,j,k,nnew))+
869 Lscale(i,j,k)*std::sqrt(tke(i,j,k,nstp)));
870 Akv(i,j,k)=Akv_bak+Sm*ql;
871 for (int n=0; n<NCONS; n++) {
872 Akt(i,j,k,n)=Akt_bak+Sh*ql;
873 }
874
875 // Compute vertical mixing (m2/s) coefficients of turbulent kinetic
876 // energy and generic statistical field.
877
878 Akk(i,j,k)=Akk_bak+Sm*ql/gls_sigk;
879 Akp(i,j,k)=Akp_bak+Sm*ql*ogls_sigp;
880
881 // Save limited length scale.
882 Lscale(i,j,k)=Ls_lmt;
883 });
884
885 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
886 {
887 Real Zob_min = std::max(ZoBot(i,j,0), 0.0001_rt);
888 Akv(i,j,N+1)=Akv_bak+L_sft*Zos_eff*gls_cmu0*
889 std::sqrt(tke(i,j,N+1,nnew));
890 Akv(i,j,0)=Akv_bak+vonKar*Zob_min*gls_cmu0*
891 std::sqrt(tke(i,j,0,nnew));
892
893 Akk(i,j,N+1)=Akk_bak+Akv(i,j,N+1)/gls_sigk;
894 Akk(i,j,0)=Akk_bak+Akv(i,j,0)/gls_sigk;
895 Akp(i,j,N+1)=Akp_bak+Akv(i,j,N+1)*ogls_sigp;
896 Akp(i,j,0)=Akp_bak+Akv(i,j,0)/gls_sigp_cb;
897
898 for (int n=0; n<NCONS; n++) {
899 Akt(i,j,N+1,n) = Akt_bak;
900 Akt(i,j,0,n) = Akt_bak;
901 }
902 });
903 }
904
905 for (int icomp=0; icomp<3; icomp++) {
906 FillPatch(lev, t_old[lev], *mf_tke, GetVecOfPtrs(vec_tke), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
907 FillPatch(lev, t_old[lev], *mf_gls, GetVecOfPtrs(vec_gls), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
908 }
909 for (int icomp=0; icomp<NCONS; icomp++) {
910 FillPatch(lev, t_old[lev], *mf_Akt, GetVecOfPtrs(vec_Akt), BCVars::zvel_bc, BdyVars::null, icomp, false, false);
911 }
912 FillPatchNoBC(lev, t_old[lev], *mf_Akv, GetVecOfPtrs(vec_Akv), BdyVars::null);
913 FillPatchNoBC(lev, t_old[lev], *mf_Akp, GetVecOfPtrs(vec_Akp), BdyVars::null);
914 FillPatchNoBC(lev, t_old[lev], *mf_Akk, GetVecOfPtrs(vec_Akk), BdyVars::null);
915}
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:1178
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:366
REMORA::vec_ZoBot
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ZoBot
Bottom roughness length [m], defined at rho points.
Definition REMORA.H:328
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:414
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:416
REMORA::vec_sustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sustr
Surface stress in the u direction.
Definition REMORA.H:293
REMORA::vec_Lscale
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Lscale
Vertical mixing turbulent length scale.
Definition REMORA.H:418
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:238
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:258
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:200
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:216
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:401
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1166
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:295
REMORA::vec_Huon
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Huon
u-volume flux (3D)
Definition REMORA.H:240
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:218
REMORA::t_new
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1156
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:420
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:214
REMORA::vec_bustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bustr
Bottom stress in the u direction.
Definition REMORA.H:331
REMORA::vec_bvstr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvstr
Bottom stress in the v direction.
Definition REMORA.H:333
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:368
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:256
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1158
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:270
REMORA::vec_Hvom
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hvom
v-volume flux (3D)
Definition REMORA.H:242
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:422
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
SolverChoice::gls_stability_type
GLS_StabilityType gls_stability_type
Definition REMORA_DataStruct.H:601
SolverChoice::Akv_bak
amrex::Real Akv_bak
Definition REMORA_DataStruct.H:681
SolverChoice::gls_sigp
amrex::Real gls_sigp
Definition REMORA_DataStruct.H:676
SolverChoice::gls_m
amrex::Real gls_m
Definition REMORA_DataStruct.H:665
SolverChoice::my_B2
amrex::Real my_B2
Definition REMORA_DataStruct.H:703
SolverChoice::my_A1
amrex::Real my_A1
Definition REMORA_DataStruct.H:700
SolverChoice::gls_sigk
amrex::Real gls_sigk
Definition REMORA_DataStruct.H:675
SolverChoice::gls_L3
amrex::Real gls_L3
Definition REMORA_DataStruct.H:692
SolverChoice::gls_cmu0
amrex::Real gls_cmu0
Definition REMORA_DataStruct.H:670
SolverChoice::gls_L6
amrex::Real gls_L6
Definition REMORA_DataStruct.H:695
SolverChoice::gls_L2
amrex::Real gls_L2
Definition REMORA_DataStruct.H:691
SolverChoice::Akk_bak
amrex::Real Akk_bak
Definition REMORA_DataStruct.H:679
SolverChoice::gls_L1
amrex::Real gls_L1
Definition REMORA_DataStruct.H:690
SolverChoice::gls_n
amrex::Real gls_n
Definition REMORA_DataStruct.H:666
SolverChoice::Akt_bak
amrex::Real Akt_bak
Definition REMORA_DataStruct.H:682
SolverChoice::gls_c3m
amrex::Real gls_c3m
Definition REMORA_DataStruct.H:673
SolverChoice::my_A2
amrex::Real my_A2
Definition REMORA_DataStruct.H:701
SolverChoice::gls_Gh0
amrex::Real gls_Gh0
Definition REMORA_DataStruct.H:685
SolverChoice::gls_Ghmin
amrex::Real gls_Ghmin
Definition REMORA_DataStruct.H:687
SolverChoice::my_C1
amrex::Real my_C1
Definition REMORA_DataStruct.H:704
SolverChoice::gls_c1
amrex::Real gls_c1
Definition REMORA_DataStruct.H:671
SolverChoice::gls_L5
amrex::Real gls_L5
Definition REMORA_DataStruct.H:694
SolverChoice::gls_L8
amrex::Real gls_L8
Definition REMORA_DataStruct.H:697
SolverChoice::gls_Kmin
amrex::Real gls_Kmin
Definition REMORA_DataStruct.H:667
SolverChoice::gls_c3p
amrex::Real gls_c3p
Definition REMORA_DataStruct.H:674
SolverChoice::gls_p
amrex::Real gls_p
Definition REMORA_DataStruct.H:664
SolverChoice::gls_L4
amrex::Real gls_L4
Definition REMORA_DataStruct.H:693
SolverChoice::my_B1
amrex::Real my_B1
Definition REMORA_DataStruct.H:702
SolverChoice::gls_L7
amrex::Real gls_L7
Definition REMORA_DataStruct.H:696
SolverChoice::gls_Ghcri
amrex::Real gls_Ghcri
Definition REMORA_DataStruct.H:686
SolverChoice::gls_Pmin
amrex::Real gls_Pmin
Definition REMORA_DataStruct.H:668
SolverChoice::gls_c2
amrex::Real gls_c2
Definition REMORA_DataStruct.H:672
SolverChoice::Akp_bak
amrex::Real Akp_bak
Definition REMORA_DataStruct.H:680
SolverChoice::gls_E2
amrex::Real gls_E2
Definition REMORA_DataStruct.H:688