diff --git a/classROMSX.html b/classROMSX.html index ac4c76df..967a6fe3 100644 --- a/classROMSX.html +++ b/classROMSX.html @@ -801,663 +801,662 @@

22  U_old.FillBoundary(geom[lev].periodicity());
23  V_old.FillBoundary(geom[lev].periodicity());
24  W_old.FillBoundary(geom[lev].periodicity());
-
25  MultiFab::Copy(S_new,S_old,0,0,S_new.nComp(),S_new.nGrowVect());
-
26  MultiFab::Copy(U_new,U_old,0,0,U_new.nComp(),U_new.nGrowVect());
-
27  MultiFab::Copy(V_new,V_old,0,0,V_new.nComp(),V_new.nGrowVect());
-
28  MultiFab::Copy(W_new,W_old,0,0,W_new.nComp(),W_new.nGrowVect());
-
29 
-
30  ////////// //pre_step3d corrections to boundaries
-
31 
-
32  const BoxArray& ba = S_old.boxArray();
-
33  const DistributionMapping& dm = S_old.DistributionMap();
-
34 
-
35  int nvars = S_old.nComp();
-
36 
-
37  const int ncomp = 1;
-
38  const int nrhs = ncomp-1;
-
39  const int nnew = ncomp-1;
-
40  const int nstp = ncomp-1;
+
25 
+
26  ////////// //pre_step3d corrections to boundaries
+
27 
+
28  const BoxArray& ba = S_old.boxArray();
+
29  const DistributionMapping& dm = S_old.DistributionMap();
+
30 
+
31  int nvars = S_old.nComp();
+
32 
+
33  const int ncomp = 1;
+
34  const int nrhs = ncomp-1;
+
35  const int nnew = ncomp-1;
+
36  const int nstp = ncomp-1;
+
37 
+
38  // Place-holder for source array -- for now just set to 0
+
39  MultiFab source(ba,dm,nvars,1);
+
40  source.setVal(0.0);
41 
-
42  // Place-holder for source array -- for now just set to 0
-
43  MultiFab source(ba,dm,nvars,1);
-
44  source.setVal(0.0);
+
42  //-----------------------------------------------------------------------
+
43  // Time step momentum equation
+
44  //-----------------------------------------------------------------------
45 
-
46  //-----------------------------------------------------------------------
-
47  // Time step momentum equation
-
48  //-----------------------------------------------------------------------
-
49 
-
50  //Only used locally, probably should be rearranged into FArrayBox declaration
-
51  MultiFab mf_AK(ba,dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
-
52  MultiFab mf_DC(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
-
53  MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
-
54  std::unique_ptr<MultiFab>& mf_Hz = vec_Hz[lev];
-
55  std::unique_ptr<MultiFab>& mf_z_r = vec_z_r[lev];
-
56  std::unique_ptr<MultiFab>& mf_z_w = vec_z_w[lev];
-
57  std::unique_ptr<MultiFab>& mf_h = vec_hOfTheConfusingName[lev];
-
58  //Consider passing these into the advance function or renaming relevant things
-
59  /*
-
60  MultiFab mf_u(ba,dm,1,IntVect(NGROW,NGROW,0));
-
61  MultiFab mf_v(ba,dm,1,IntVect(NGROW,NGROW,0));
-
62  MultiFab mf_uold(ba,dm,1,IntVect(NGROW,NGROW,0));
-
63  MultiFab mf_vold(ba,dm,1,IntVect(NGROW,NGROW,0));*/
-
64  MultiFab mf_u(U_new, amrex::make_alias, 0, 1);
-
65  MultiFab mf_v(V_new, amrex::make_alias, 0, 1);
-
66  MultiFab mf_uold(U_old, amrex::make_alias, 0, 1);
-
67  MultiFab mf_vold(V_old, amrex::make_alias, 0, 1);
-
68  MultiFab mf_w(ba,dm,1,IntVect(NGROW,NGROW,0));
-
69  MultiFab mf_pden(ba,dm,1,IntVect(NGROW,NGROW,0));
-
70  MultiFab mf_rho(ba,dm,1,IntVect(NGROW,NGROW,0));
-
71  MultiFab mf_rhoS(ba,dm,1,IntVect(NGROW,NGROW,0));
-
72  MultiFab mf_rhoA(ba,dm,1,IntVect(NGROW,NGROW,0));
-
73  std::unique_ptr<MultiFab>& mf_ru = vec_ru[lev];
-
74  std::unique_ptr<MultiFab>& mf_rv = vec_rv[lev];
-
75  std::unique_ptr<MultiFab>& mf_rufrc = vec_rufrc[lev];
-
76  std::unique_ptr<MultiFab>& mf_rvfrc = vec_rvfrc[lev];
-
77  std::unique_ptr<MultiFab>& mf_sustr = vec_sustr[lev];
-
78  std::unique_ptr<MultiFab>& mf_svstr = vec_svstr[lev];
-
79  std::unique_ptr<MultiFab>& mf_rdrag = vec_rdrag[lev];
-
80  std::unique_ptr<MultiFab>& mf_bustr = vec_bustr[lev];
-
81  std::unique_ptr<MultiFab>& mf_bvstr = vec_bvstr[lev];
-
82  std::unique_ptr<MultiFab>& mf_ubar = vec_ubar[lev];
-
83  std::unique_ptr<MultiFab>& mf_vbar = vec_vbar[lev];
-
84  MultiFab mf_temp(S_new, amrex::make_alias, Temp_comp, 1);
-
85 #ifdef ROMSX_USE_SALINITY
-
86  MultiFab mf_salt(S_new, amrex::make_alias, Salt_comp, 1);
-
87 #else
-
88  MultiFab mf_salt(S_new, amrex::make_alias, Temp_comp, 1);
-
89 #endif
-
90  MultiFab mf_tempold(S_old, amrex::make_alias, Temp_comp, 1);
-
91 #ifdef ROMSX_USE_SALINITY
-
92  MultiFab mf_saltold(S_old, amrex::make_alias, Salt_comp, 1);
-
93 #else
-
94  MultiFab mf_saltold(S_old, amrex::make_alias, Temp_comp, 1);
-
95 #endif
-
96  MultiFab mf_rw(ba,dm,1,IntVect(NGROW,NGROW,0));
-
97  MultiFab mf_W(ba,dm,1,IntVect(NGROW+1,NGROW+1,0));
-
98 
-
99  std::unique_ptr<MultiFab>& mf_visc2_p = vec_visc2_p[lev];
-
100  std::unique_ptr<MultiFab>& mf_visc2_r = vec_visc2_r[lev];
-
101  std::unique_ptr<MultiFab>& mf_diff2_temp = vec_diff2_temp[lev];
-
102  std::unique_ptr<MultiFab>& mf_diff2_salt = vec_diff2_salt[lev];
-
103 
-
104  // We need to set these because otherwise in the first call to romsx_advance we may
-
105  // read uninitialized data on ghost values in setting the bc's on the velocities
-
106  /*
-
107  mf_u.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
108  mf_v.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
109  mf_uold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
110  mf_vold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));*/
-
111  mf_pden.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
112  mf_rho.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
113  mf_rhoS.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
114  mf_rhoA.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
115  mf_w.setVal(0);
-
116  mf_DC.setVal(0);
-
117  mf_w.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
118 
-
119  MultiFab::Copy(mf_u,U_new,0,0,U_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
120  MultiFab::Copy(mf_v,V_new,0,0,V_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
121  MultiFab::Copy(mf_uold,U_old,0,0,U_old.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
122  MultiFab::Copy(mf_vold,V_old,0,0,V_old.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
123  MultiFab::Copy(mf_w,W_new,0,0,W_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
124  MultiFab::Copy(mf_W,S_old,Omega_comp,0,mf_W.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
-
125 
-
126  mf_u.FillBoundary(geom[lev].periodicity());
-
127  mf_v.FillBoundary(geom[lev].periodicity());
-
128  mf_uold.FillBoundary(geom[lev].periodicity());
-
129  mf_vold.FillBoundary(geom[lev].periodicity());
-
130  mf_w.FillBoundary(geom[lev].periodicity());
-
131  mf_W.FillBoundary(geom[lev].periodicity());
-
132  mf_tempold.FillBoundary(geom[lev].periodicity());
-
133  mf_temp.FillBoundary(geom[lev].periodicity());
-
134  mf_saltold.FillBoundary(geom[lev].periodicity());
-
135  mf_salt.FillBoundary(geom[lev].periodicity());
-
136 
-
137  mf_rw.setVal(0.0);
-
138  mf_W.setVal(0.0);
-
139  U_old.FillBoundary(geom[lev].periodicity());
-
140  V_old.FillBoundary(geom[lev].periodicity());
-
141  mf_rufrc->setVal(0);
-
142  mf_rvfrc->setVal(0);
-
143 
-
144  int iic = istep[lev];
-
145  int ntfirst = 0;
-
146  if(iic==ntfirst&&false)
-
147  MultiFab::Copy(S_new,S_old,0,0,S_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,NGROW)));
-
148  set_smflux(lev,t_old[lev]);
-
149 
-
150  auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ
-
151 
-
152  const auto prob_lo = Geom(lev).ProbLoArray();
-
153  const auto dxi = Geom(lev).InvCellSizeArray();
-
154  const auto dx = Geom(lev).CellSizeArray();
-
155  const int Mm = Geom(lev).Domain().size()[1];
-
156  auto geomdata = Geom(lev).data();
-
157 
-
158  //MFIter::allowMultipleMFIters(true);
-
159  for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-
160  {
-
161  Array4<Real> const& DC = mf_DC.array(mfi);
-
162  Array4<Real> const& Akv = (vec_Akv[lev])->array(mfi);
-
163  Array4<Real> const& h = (vec_hOfTheConfusingName[lev])->array(mfi);
-
164  Array4<Real> const& Hz = (vec_Hz[lev])->array(mfi);
-
165  Array4<Real> const& Huon = (vec_Huon[lev])->array(mfi);
-
166  Array4<Real> const& Hvom = (vec_Hvom[lev])->array(mfi);
-
167  Array4<Real> const& z_r = (mf_z_r)->array(mfi);
-
168  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
-
169  Array4<Real> const& uold = (mf_uold).array(mfi);
-
170  Array4<Real> const& vold = (mf_vold).array(mfi);
-
171  Array4<Real> const& u = (mf_u).array(mfi);
-
172  Array4<Real> const& v = (mf_v).array(mfi);
-
173  Array4<Real> const& pden = (mf_pden).array(mfi);
-
174  Array4<Real> const& rho = (mf_rho).array(mfi);
-
175  Array4<Real> const& rhoA = (mf_rhoA).array(mfi);
-
176  Array4<Real> const& rhoS = (mf_rhoS).array(mfi);
-
177  Array4<Real> const& tempold = (mf_tempold).array(mfi);
-
178  Array4<Real> const& saltold = (mf_saltold).array(mfi);
-
179  Array4<Real> const& temp = (mf_temp).array(mfi);
-
180  Array4<Real> const& salt = (mf_salt).array(mfi);
-
181  Array4<Real> const& tempstore = (vec_t3[lev])->array(mfi);
-
182  Array4<Real> const& saltstore = (vec_s3[lev])->array(mfi);
-
183  Array4<Real> const& ru = (mf_ru)->array(mfi);
-
184  Array4<Real> const& rv = (mf_rv)->array(mfi);
-
185  Array4<Real> const& rufrc = (mf_rufrc)->array(mfi);
-
186  Array4<Real> const& rvfrc = (mf_rvfrc)->array(mfi);
-
187  Array4<Real> const& W = (mf_W).array(mfi);
-
188  Array4<Real> const& sustr = (mf_sustr)->array(mfi);
-
189  Array4<Real> const& svstr = (mf_svstr)->array(mfi);
-
190  Array4<Real> const& rdrag = (mf_rdrag)->array(mfi);
-
191  Array4<Real> const& bustr = (mf_bustr)->array(mfi);
-
192  Array4<Real> const& bvstr = (mf_bvstr)->array(mfi);
-
193  Array4<Real> const& ubar = (mf_ubar)->array(mfi);
-
194  Array4<Real> const& vbar = (mf_vbar)->array(mfi);
-
195  Array4<Real> const& visc2_p = (mf_visc2_p)->array(mfi);
-
196  Array4<Real> const& visc2_r = (mf_visc2_r)->array(mfi);
-
197  Array4<Real> const& diff2_salt = (mf_diff2_salt)->array(mfi);
-
198  Array4<Real> const& diff2_temp = (mf_diff2_temp)->array(mfi);
-
199 
-
200  Box bx = mfi.tilebox();
-
201  Box gbx = mfi.growntilebox();
-
202  Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
-
203  Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
-
204  Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
-
205  //copy the tilebox
-
206  //Box gbx1 = bx;
-
207  //Box gbx11 = bx;
-
208  //Box gbx2 = bx;
-
209  //TODO: adjust for tiling
-
210  //Box gbx3uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-3,bx.smallEnd(1)-3,bx.smallEnd(2))),
-
211  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+2,bx.bigEnd(1)+2,bx.bigEnd(2))));
-
212  //Box gbx2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-2,bx.smallEnd(1)-2,bx.smallEnd(2))),
-
213  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+1,bx.bigEnd(1)+1,bx.bigEnd(2))));
-
214  Box ubx = surroundingNodes(bx,0);
-
215  Box vbx = surroundingNodes(bx,1);
-
216  //make only gbx be grown to match multifabs
-
217  //gbx2.grow(IntVect(NGROW,NGROW,0));
-
218  //gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
-
219  //gbx11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
-
220 
-
221  Box bxD = bx;
-
222  bxD.makeSlab(2,0);
-
223  Box gbx1D = gbx1;
-
224  gbx1D.makeSlab(2,0);
-
225  Box gbx2D = gbx2;
-
226  gbx2D.makeSlab(2,0);
-
227  //gbx1D.grow(IntVect(NGROW-1,NGROW-1,0));
-
228  //gbx2D.grow(IntVect(NGROW,NGROW,0));
-
229 
-
230  FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena()); //3D
-
231  FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena()); //3D
-
232  FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena()); //3D
-
233  FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
-
234  FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
-
235  FArrayBox fab_pn(gbx2D,1,amrex::The_Async_Arena());
-
236  FArrayBox fab_pm(gbx2D,1,amrex::The_Async_Arena());
-
237  FArrayBox fab_on_u(gbx2D,1,amrex::The_Async_Arena());
-
238  FArrayBox fab_om_v(gbx2D,1,amrex::The_Async_Arena());
-
239  FArrayBox fab_om_u(gbx2D,1,amrex::The_Async_Arena());
-
240  FArrayBox fab_on_v(gbx2D,1,amrex::The_Async_Arena());
-
241  FArrayBox fab_om_r(gbx2D,1,amrex::The_Async_Arena());
-
242  FArrayBox fab_on_r(gbx2D,1,amrex::The_Async_Arena());
-
243  FArrayBox fab_om_p(gbx2D,1,amrex::The_Async_Arena());
-
244  FArrayBox fab_on_p(gbx2D,1,amrex::The_Async_Arena());
-
245  FArrayBox fab_pmon_u(gbx2D,1,amrex::The_Async_Arena());
-
246  FArrayBox fab_pnom_u(gbx2D,1,amrex::The_Async_Arena());
-
247  FArrayBox fab_pmon_v(gbx2D,1,amrex::The_Async_Arena());
-
248  FArrayBox fab_pnom_v(gbx2D,1,amrex::The_Async_Arena());
-
249  FArrayBox fab_fomn(gbx2D,1,amrex::The_Async_Arena());
-
250  //FArrayBox fab_oHz(gbx11,1,amrex::The_Async_Arena());
-
251 
-
252  auto FC=fab_FC.array();
-
253  auto FX=fab_FX.array();
-
254  auto FE=fab_FE.array();
-
255  auto pn=fab_pn.array();
-
256  auto pm=fab_pm.array();
-
257  auto on_u=fab_on_u.array();
-
258  auto om_v=fab_om_v.array();
-
259  auto om_u=fab_om_u.array();
-
260  auto on_v=fab_on_v.array();
-
261  auto om_r=fab_om_r.array();
-
262  auto on_r=fab_on_r.array();
-
263  auto om_p=fab_om_p.array();
-
264  auto on_p=fab_on_p.array();
-
265  auto pmon_u=fab_pmon_u.array();
-
266  auto pnom_u=fab_pnom_u.array();
-
267  auto pmon_v=fab_pmon_v.array();
-
268  auto pnom_v=fab_pnom_v.array();
-
269  auto fomn=fab_fomn.array();
-
270 
-
271  //From ana_grid.h and metrics.F
-
272  amrex::ParallelFor(gbx2D,
-
273  [=] AMREX_GPU_DEVICE (int i, int j, int )
-
274  {
-
275  pm(i,j,0)=dxi[0];
-
276  pn(i,j,0)=dxi[1];
-
277  //defined UPWELLING
-
278  Real f0=-8.26e-5;
-
279  Real beta=0.0;
-
280  Real Esize=1000*(Mm);
-
281  Real y = prob_lo[1] + (j + 0.5) * dx[1];
-
282  Real f=f0+beta*(y-.5*Esize);
-
283  fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0)));
-
284  });
-
285 
-
286  amrex::ParallelFor(gbx2D,
-
287  [=] AMREX_GPU_DEVICE (int i, int j, int )
-
288  {
-
289  //Note: are the comment definitons right? Don't seem to match metrics.f90
-
290  om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j))
-
291  on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j))
-
292  om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j)
-
293  on_r(i,j,0)=1.0/dxi[1]; // 1/pn(i,j)
-
294  //todo: om_p on_p
-
295  om_p(i,j,0)=1.0/dxi[0]; // 4/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
-
296  on_p(i,j,0)=1.0/dxi[1]; // 4/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
-
297  on_v(i,j,0)=1.0/dxi[1]; // 2/(pn(i-1,j)+pn(i,j))
-
298  om_u(i,j,0)=1.0/dxi[0]; // 2/(pm(i-1,j)+pm(i,j))
-
299  pmon_u(i,j,0)=1.0; // (pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
-
300  pnom_u(i,j,0)=1.0; // (pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
-
301  pmon_v(i,j,0)=1.0; // (pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j))
-
302  pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j))
-
303  });
-
304 
-
305  amrex::ParallelFor(gbx2,
-
306  [=] AMREX_GPU_DEVICE (int i, int j, int k)
-
307  {
-
308  Huon(i,j,k,0)=0.0;
-
309  Hvom(i,j,k,0)=0.0;
-
310  });
-
311 
-
312  // Set bottom stress as defined in set_vbx.F
-
313  amrex::ParallelFor(gbx1D,
-
314  [=] AMREX_GPU_DEVICE (int i, int j, int )
-
315  {
-
316  //DO j=Jstr,Jend
-
317  // DO i=IstrU,Iend
-
318  bustr(i,j,0) = 0.5 * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0,nrhs));
-
319  //DO j=JstrV,Jend
-
320  // DO i=Istr,Iend
-
321  bvstr(i,j,0) = 0.5 * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0,nrhs));
-
322  });
-
323 
-
324  // updates Huon/Hvom
-
325  set_massflux_3d(gbx2,1,0,uold,Huon,Hz,on_u,nnew);
-
326  set_massflux_3d(gbx2,0,1,vold,Hvom,Hz,om_v,nnew);
-
327 
-
328  rho_eos(gbx2,temp,salt,rho,rhoA,rhoS,pden,Hz,z_w,h,nrhs,N);
-
329  }
-
330 
-
331  if(solverChoice.use_prestep) {
-
332  prestep(lev, mf_uold, mf_vold, mf_u, mf_v, mf_ru, mf_rv, mf_tempold, mf_saltold,
-
333  mf_temp, mf_salt, mf_Hz, vec_Akv[lev], vec_Akt[lev], vec_Huon[lev], vec_Hvom[lev], mf_W,
-
334  mf_DC, vec_t3[lev], vec_s3[lev], mf_z_r, mf_z_w, mf_h, mf_sustr, mf_svstr, mf_bustr,
-
335  mf_bvstr, iic, ntfirst, nnew, nstp, nrhs, N, dt_lev);
-
336  }
+
46  //Only used locally, probably should be rearranged into FArrayBox declaration
+
47  MultiFab mf_AK(ba,dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate
+
48  MultiFab mf_DC(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
+
49  MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate
+
50  std::unique_ptr<MultiFab>& mf_Hz = vec_Hz[lev];
+
51  std::unique_ptr<MultiFab>& mf_z_r = vec_z_r[lev];
+
52  std::unique_ptr<MultiFab>& mf_z_w = vec_z_w[lev];
+
53  std::unique_ptr<MultiFab>& mf_h = vec_hOfTheConfusingName[lev];
+
54  //Consider passing these into the advance function or renaming relevant things
+
55  /*
+
56  MultiFab mf_u(ba,dm,1,IntVect(NGROW,NGROW,0));
+
57  MultiFab mf_v(ba,dm,1,IntVect(NGROW,NGROW,0));
+
58  MultiFab mf_uold(ba,dm,1,IntVect(NGROW,NGROW,0));
+
59  MultiFab mf_vold(ba,dm,1,IntVect(NGROW,NGROW,0));*/
+
60  MultiFab mf_u(U_new, amrex::make_alias, 0, 1);
+
61  MultiFab mf_v(V_new, amrex::make_alias, 0, 1);
+
62  MultiFab mf_uold(U_old, amrex::make_alias, 0, 1);
+
63  MultiFab mf_vold(V_old, amrex::make_alias, 0, 1);
+
64  MultiFab mf_w(ba,dm,1,IntVect(NGROW,NGROW,0));
+
65  MultiFab mf_pden(ba,dm,1,IntVect(NGROW,NGROW,0));
+
66  MultiFab mf_rho(ba,dm,1,IntVect(NGROW,NGROW,0));
+
67  MultiFab mf_rhoS(ba,dm,1,IntVect(NGROW,NGROW,0));
+
68  MultiFab mf_rhoA(ba,dm,1,IntVect(NGROW,NGROW,0));
+
69  std::unique_ptr<MultiFab>& mf_ru = vec_ru[lev];
+
70  std::unique_ptr<MultiFab>& mf_rv = vec_rv[lev];
+
71  std::unique_ptr<MultiFab>& mf_rufrc = vec_rufrc[lev];
+
72  std::unique_ptr<MultiFab>& mf_rvfrc = vec_rvfrc[lev];
+
73  std::unique_ptr<MultiFab>& mf_sustr = vec_sustr[lev];
+
74  std::unique_ptr<MultiFab>& mf_svstr = vec_svstr[lev];
+
75  std::unique_ptr<MultiFab>& mf_rdrag = vec_rdrag[lev];
+
76  std::unique_ptr<MultiFab>& mf_bustr = vec_bustr[lev];
+
77  std::unique_ptr<MultiFab>& mf_bvstr = vec_bvstr[lev];
+
78  std::unique_ptr<MultiFab>& mf_ubar = vec_ubar[lev];
+
79  std::unique_ptr<MultiFab>& mf_vbar = vec_vbar[lev];
+
80  MultiFab mf_temp(S_new, amrex::make_alias, Temp_comp, 1);
+
81 #ifdef ROMSX_USE_SALINITY
+
82  MultiFab mf_salt(S_new, amrex::make_alias, Salt_comp, 1);
+
83 #else
+
84  MultiFab mf_salt(S_new, amrex::make_alias, Temp_comp, 1);
+
85 #endif
+
86  MultiFab mf_tempold(S_old, amrex::make_alias, Temp_comp, 1);
+
87 #ifdef ROMSX_USE_SALINITY
+
88  MultiFab mf_saltold(S_old, amrex::make_alias, Salt_comp, 1);
+
89 #else
+
90  MultiFab mf_saltold(S_old, amrex::make_alias, Temp_comp, 1);
+
91 #endif
+
92  MultiFab mf_rw(ba,dm,1,IntVect(NGROW,NGROW,0));
+
93  MultiFab mf_W(ba,dm,1,IntVect(NGROW+1,NGROW+1,0));
+
94 
+
95  std::unique_ptr<MultiFab>& mf_visc2_p = vec_visc2_p[lev];
+
96  std::unique_ptr<MultiFab>& mf_visc2_r = vec_visc2_r[lev];
+
97  std::unique_ptr<MultiFab>& mf_diff2_temp = vec_diff2_temp[lev];
+
98  std::unique_ptr<MultiFab>& mf_diff2_salt = vec_diff2_salt[lev];
+
99 
+
100  // We need to set these because otherwise in the first call to romsx_advance we may
+
101  // read uninitialized data on ghost values in setting the bc's on the velocities
+
102  /*
+
103  mf_u.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
104  mf_v.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
105  mf_uold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
106  mf_vold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));*/
+
107  mf_pden.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
108  mf_rho.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
109  mf_rhoS.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
110  mf_rhoA.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
111  mf_w.setVal(0);
+
112  mf_DC.setVal(0);
+
113  mf_w.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
114 
+
115  MultiFab::Copy(mf_u,U_new,0,0,U_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
116  MultiFab::Copy(mf_v,V_new,0,0,V_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
117  MultiFab::Copy(mf_uold,U_old,0,0,U_old.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
118  MultiFab::Copy(mf_vold,V_old,0,0,V_old.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
119  MultiFab::Copy(mf_w,W_new,0,0,W_new.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
120  MultiFab::Copy(mf_W,S_old,Omega_comp,0,mf_W.nComp(),IntVect(AMREX_D_DECL(NGROW,NGROW,0)));
+
121 
+
122  mf_u.FillBoundary(geom[lev].periodicity());
+
123  mf_v.FillBoundary(geom[lev].periodicity());
+
124  mf_uold.FillBoundary(geom[lev].periodicity());
+
125  mf_vold.FillBoundary(geom[lev].periodicity());
+
126  mf_w.FillBoundary(geom[lev].periodicity());
+
127  mf_W.FillBoundary(geom[lev].periodicity());
+
128  mf_tempold.FillBoundary(geom[lev].periodicity());
+
129  mf_temp.FillBoundary(geom[lev].periodicity());
+
130  mf_saltold.FillBoundary(geom[lev].periodicity());
+
131  mf_salt.FillBoundary(geom[lev].periodicity());
+
132 
+
133  mf_rw.setVal(0.0);
+
134  mf_W.setVal(0.0);
+
135  U_old.FillBoundary(geom[lev].periodicity());
+
136  V_old.FillBoundary(geom[lev].periodicity());
+
137  mf_rufrc->setVal(0);
+
138  mf_rvfrc->setVal(0);
+
139 
+
140  int iic = istep[lev];
+
141  int ntfirst = 0;
+
142  if(iic==ntfirst)
+
143  MultiFab::Copy(S_new,S_old,0,0,S_new.nComp(),S_new.nGrowVect());
+
144  MultiFab::Copy(U_new,U_old,0,0,U_new.nComp(),U_new.nGrowVect());
+
145  MultiFab::Copy(V_new,V_old,0,0,V_new.nComp(),V_new.nGrowVect());
+
146  MultiFab::Copy(W_new,W_old,0,0,W_new.nComp(),W_new.nGrowVect());
+
147  set_smflux(lev,t_old[lev]);
+
148 
+
149  auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ
+
150 
+
151  const auto prob_lo = Geom(lev).ProbLoArray();
+
152  const auto dxi = Geom(lev).InvCellSizeArray();
+
153  const auto dx = Geom(lev).CellSizeArray();
+
154  const int Mm = Geom(lev).Domain().size()[1];
+
155  auto geomdata = Geom(lev).data();
+
156 
+
157  //MFIter::allowMultipleMFIters(true);
+
158  for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+
159  {
+
160  Array4<Real> const& DC = mf_DC.array(mfi);
+
161  Array4<Real> const& Akv = (vec_Akv[lev])->array(mfi);
+
162  Array4<Real> const& h = (vec_hOfTheConfusingName[lev])->array(mfi);
+
163  Array4<Real> const& Hz = (vec_Hz[lev])->array(mfi);
+
164  Array4<Real> const& Huon = (vec_Huon[lev])->array(mfi);
+
165  Array4<Real> const& Hvom = (vec_Hvom[lev])->array(mfi);
+
166  Array4<Real> const& z_r = (mf_z_r)->array(mfi);
+
167  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
+
168  Array4<Real> const& uold = (mf_uold).array(mfi);
+
169  Array4<Real> const& vold = (mf_vold).array(mfi);
+
170  Array4<Real> const& u = (mf_u).array(mfi);
+
171  Array4<Real> const& v = (mf_v).array(mfi);
+
172  Array4<Real> const& pden = (mf_pden).array(mfi);
+
173  Array4<Real> const& rho = (mf_rho).array(mfi);
+
174  Array4<Real> const& rhoA = (mf_rhoA).array(mfi);
+
175  Array4<Real> const& rhoS = (mf_rhoS).array(mfi);
+
176  Array4<Real> const& tempold = (mf_tempold).array(mfi);
+
177  Array4<Real> const& saltold = (mf_saltold).array(mfi);
+
178  Array4<Real> const& temp = (mf_temp).array(mfi);
+
179  Array4<Real> const& salt = (mf_salt).array(mfi);
+
180  Array4<Real> const& tempstore = (vec_t3[lev])->array(mfi);
+
181  Array4<Real> const& saltstore = (vec_s3[lev])->array(mfi);
+
182  Array4<Real> const& ru = (mf_ru)->array(mfi);
+
183  Array4<Real> const& rv = (mf_rv)->array(mfi);
+
184  Array4<Real> const& rufrc = (mf_rufrc)->array(mfi);
+
185  Array4<Real> const& rvfrc = (mf_rvfrc)->array(mfi);
+
186  Array4<Real> const& W = (mf_W).array(mfi);
+
187  Array4<Real> const& sustr = (mf_sustr)->array(mfi);
+
188  Array4<Real> const& svstr = (mf_svstr)->array(mfi);
+
189  Array4<Real> const& rdrag = (mf_rdrag)->array(mfi);
+
190  Array4<Real> const& bustr = (mf_bustr)->array(mfi);
+
191  Array4<Real> const& bvstr = (mf_bvstr)->array(mfi);
+
192  Array4<Real> const& ubar = (mf_ubar)->array(mfi);
+
193  Array4<Real> const& vbar = (mf_vbar)->array(mfi);
+
194  Array4<Real> const& visc2_p = (mf_visc2_p)->array(mfi);
+
195  Array4<Real> const& visc2_r = (mf_visc2_r)->array(mfi);
+
196  Array4<Real> const& diff2_salt = (mf_diff2_salt)->array(mfi);
+
197  Array4<Real> const& diff2_temp = (mf_diff2_temp)->array(mfi);
+
198 
+
199  Box bx = mfi.tilebox();
+
200  Box gbx = mfi.growntilebox();
+
201  Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
+
202  Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
+
203  Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
+
204  //copy the tilebox
+
205  //Box gbx1 = bx;
+
206  //Box gbx11 = bx;
+
207  //Box gbx2 = bx;
+
208  //TODO: adjust for tiling
+
209  //Box gbx3uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-3,bx.smallEnd(1)-3,bx.smallEnd(2))),
+
210  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+2,bx.bigEnd(1)+2,bx.bigEnd(2))));
+
211  //Box gbx2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-2,bx.smallEnd(1)-2,bx.smallEnd(2))),
+
212  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+1,bx.bigEnd(1)+1,bx.bigEnd(2))));
+
213  Box ubx = surroundingNodes(bx,0);
+
214  Box vbx = surroundingNodes(bx,1);
+
215  //make only gbx be grown to match multifabs
+
216  //gbx2.grow(IntVect(NGROW,NGROW,0));
+
217  //gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
+
218  //gbx11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
+
219 
+
220  Box bxD = bx;
+
221  bxD.makeSlab(2,0);
+
222  Box gbx1D = gbx1;
+
223  gbx1D.makeSlab(2,0);
+
224  Box gbx2D = gbx2;
+
225  gbx2D.makeSlab(2,0);
+
226  //gbx1D.grow(IntVect(NGROW-1,NGROW-1,0));
+
227  //gbx2D.grow(IntVect(NGROW,NGROW,0));
+
228 
+
229  FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena()); //3D
+
230  FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena()); //3D
+
231  FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena()); //3D
+
232  FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
+
233  FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
+
234  FArrayBox fab_pn(gbx2D,1,amrex::The_Async_Arena());
+
235  FArrayBox fab_pm(gbx2D,1,amrex::The_Async_Arena());
+
236  FArrayBox fab_on_u(gbx2D,1,amrex::The_Async_Arena());
+
237  FArrayBox fab_om_v(gbx2D,1,amrex::The_Async_Arena());
+
238  FArrayBox fab_om_u(gbx2D,1,amrex::The_Async_Arena());
+
239  FArrayBox fab_on_v(gbx2D,1,amrex::The_Async_Arena());
+
240  FArrayBox fab_om_r(gbx2D,1,amrex::The_Async_Arena());
+
241  FArrayBox fab_on_r(gbx2D,1,amrex::The_Async_Arena());
+
242  FArrayBox fab_om_p(gbx2D,1,amrex::The_Async_Arena());
+
243  FArrayBox fab_on_p(gbx2D,1,amrex::The_Async_Arena());
+
244  FArrayBox fab_pmon_u(gbx2D,1,amrex::The_Async_Arena());
+
245  FArrayBox fab_pnom_u(gbx2D,1,amrex::The_Async_Arena());
+
246  FArrayBox fab_pmon_v(gbx2D,1,amrex::The_Async_Arena());
+
247  FArrayBox fab_pnom_v(gbx2D,1,amrex::The_Async_Arena());
+
248  FArrayBox fab_fomn(gbx2D,1,amrex::The_Async_Arena());
+
249  //FArrayBox fab_oHz(gbx11,1,amrex::The_Async_Arena());
+
250 
+
251  auto FC=fab_FC.array();
+
252  auto FX=fab_FX.array();
+
253  auto FE=fab_FE.array();
+
254  auto pn=fab_pn.array();
+
255  auto pm=fab_pm.array();
+
256  auto on_u=fab_on_u.array();
+
257  auto om_v=fab_om_v.array();
+
258  auto om_u=fab_om_u.array();
+
259  auto on_v=fab_on_v.array();
+
260  auto om_r=fab_om_r.array();
+
261  auto on_r=fab_on_r.array();
+
262  auto om_p=fab_om_p.array();
+
263  auto on_p=fab_on_p.array();
+
264  auto pmon_u=fab_pmon_u.array();
+
265  auto pnom_u=fab_pnom_u.array();
+
266  auto pmon_v=fab_pmon_v.array();
+
267  auto pnom_v=fab_pnom_v.array();
+
268  auto fomn=fab_fomn.array();
+
269 
+
270  //From ana_grid.h and metrics.F
+
271  amrex::ParallelFor(gbx2D,
+
272  [=] AMREX_GPU_DEVICE (int i, int j, int )
+
273  {
+
274  pm(i,j,0)=dxi[0];
+
275  pn(i,j,0)=dxi[1];
+
276  //defined UPWELLING
+
277  Real f0=-8.26e-5;
+
278  Real beta=0.0;
+
279  Real Esize=1000*(Mm);
+
280  Real y = prob_lo[1] + (j + 0.5) * dx[1];
+
281  Real f=f0+beta*(y-.5*Esize);
+
282  fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0)));
+
283  });
+
284 
+
285  amrex::ParallelFor(gbx2D,
+
286  [=] AMREX_GPU_DEVICE (int i, int j, int )
+
287  {
+
288  //Note: are the comment definitons right? Don't seem to match metrics.f90
+
289  om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j))
+
290  on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j))
+
291  om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j)
+
292  on_r(i,j,0)=1.0/dxi[1]; // 1/pn(i,j)
+
293  //todo: om_p on_p
+
294  om_p(i,j,0)=1.0/dxi[0]; // 4/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
+
295  on_p(i,j,0)=1.0/dxi[1]; // 4/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
+
296  on_v(i,j,0)=1.0/dxi[1]; // 2/(pn(i-1,j)+pn(i,j))
+
297  om_u(i,j,0)=1.0/dxi[0]; // 2/(pm(i-1,j)+pm(i,j))
+
298  pmon_u(i,j,0)=1.0; // (pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
+
299  pnom_u(i,j,0)=1.0; // (pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
+
300  pmon_v(i,j,0)=1.0; // (pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j))
+
301  pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j))
+
302  });
+
303 
+
304  amrex::ParallelFor(gbx2,
+
305  [=] AMREX_GPU_DEVICE (int i, int j, int k)
+
306  {
+
307  Huon(i,j,k,0)=0.0;
+
308  Hvom(i,j,k,0)=0.0;
+
309  });
+
310 
+
311  // Set bottom stress as defined in set_vbx.F
+
312  amrex::ParallelFor(gbx1D,
+
313  [=] AMREX_GPU_DEVICE (int i, int j, int )
+
314  {
+
315  //DO j=Jstr,Jend
+
316  // DO i=IstrU,Iend
+
317  bustr(i,j,0) = 0.5 * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0,nrhs));
+
318  //DO j=JstrV,Jend
+
319  // DO i=Istr,Iend
+
320  bvstr(i,j,0) = 0.5 * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0,nrhs));
+
321  });
+
322 
+
323  // updates Huon/Hvom
+
324  set_massflux_3d(gbx2,1,0,uold,Huon,Hz,on_u,nnew);
+
325  set_massflux_3d(gbx2,0,1,vold,Hvom,Hz,om_v,nnew);
+
326 
+
327  rho_eos(gbx2,temp,salt,rho,rhoA,rhoS,pden,Hz,z_w,h,nrhs,N);
+
328  }
+
329 
+
330  if(solverChoice.use_prestep) {
+
331  prestep(lev, mf_uold, mf_vold, mf_u, mf_v, mf_ru, mf_rv, mf_tempold, mf_saltold,
+
332  mf_temp, mf_salt, mf_Hz, vec_Akv[lev], vec_Akt[lev], vec_Huon[lev], vec_Hvom[lev], mf_W,
+
333  mf_DC, vec_t3[lev], vec_s3[lev], mf_z_r, mf_z_w, mf_h, mf_sustr, mf_svstr, mf_bustr,
+
334  mf_bvstr, iic, ntfirst, nnew, nstp, nrhs, N, dt_lev);
+
335  }
+
336 
337 
-
338 
-
339  mf_W.FillBoundary(geom[lev].periodicity());
-
340 
-
341  for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
-
342  {
-
343  Array4<Real> const& DC = mf_DC.array(mfi);
-
344  Array4<Real> const& Akv = (vec_Akv[lev])->array(mfi);
-
345  Array4<Real> const& Hz = (vec_Hz[lev])->array(mfi);
-
346  Array4<Real> const& Huon = (vec_Huon[lev])->array(mfi);
-
347  Array4<Real> const& Hvom = (vec_Hvom[lev])->array(mfi);
-
348  Array4<Real> const& z_r = (mf_z_r)->array(mfi);
-
349  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
-
350  Array4<Real> const& uold = (mf_uold).array(mfi);
-
351  Array4<Real> const& vold = (mf_vold).array(mfi);
-
352  Array4<Real> const& u = (mf_u).array(mfi);
-
353  Array4<Real> const& v = (mf_v).array(mfi);
-
354  Array4<Real> const& pden = (mf_pden).array(mfi);
-
355  Array4<Real> const& rho = (mf_rho).array(mfi);
-
356  Array4<Real> const& rhoA = (mf_rhoA).array(mfi);
-
357  Array4<Real> const& rhoS = (mf_rhoS).array(mfi);
-
358  Array4<Real> const& tempold = (mf_tempold).array(mfi);
-
359  Array4<Real> const& saltold = (mf_saltold).array(mfi);
-
360  Array4<Real> const& temp = (mf_temp).array(mfi);
-
361  Array4<Real> const& salt = (mf_salt).array(mfi);
-
362  Array4<Real> const& tempstore = (vec_t3[lev])->array(mfi);
-
363  Array4<Real> const& saltstore = (vec_s3[lev])->array(mfi);
-
364  Array4<Real> const& ru = (mf_ru)->array(mfi);
-
365  Array4<Real> const& rv = (mf_rv)->array(mfi);
-
366  Array4<Real> const& rufrc = (mf_rufrc)->array(mfi);
-
367  Array4<Real> const& rvfrc = (mf_rvfrc)->array(mfi);
-
368  Array4<Real> const& W = (mf_W).array(mfi);
-
369  Array4<Real> const& sustr = (mf_sustr)->array(mfi);
-
370  Array4<Real> const& svstr = (mf_svstr)->array(mfi);
-
371  Array4<Real> const& rdrag = (mf_rdrag)->array(mfi);
-
372  Array4<Real> const& bustr = (mf_bustr)->array(mfi);
-
373  Array4<Real> const& bvstr = (mf_bvstr)->array(mfi);
-
374  Array4<Real> const& ubar = (mf_ubar)->array(mfi);
-
375  Array4<Real> const& vbar = (mf_vbar)->array(mfi);
-
376  Array4<Real> const& visc2_p = (mf_visc2_p)->array(mfi);
-
377  Array4<Real> const& visc2_r = (mf_visc2_r)->array(mfi);
-
378  Array4<Real> const& diff2_salt = (mf_diff2_salt)->array(mfi);
-
379  Array4<Real> const& diff2_temp = (mf_diff2_temp)->array(mfi);
-
380 
-
381  Array4<Real> const& zeta = (vec_zeta[lev])->array(mfi);
-
382  Array4<Real> const& Zt_avg1 = (vec_Zt_avg1[lev])->array(mfi);
-
383 
-
384  Box bx = mfi.tilebox();
-
385  Box tbxp1 = bx;
-
386  Box tbxp2 = bx;
-
387  Box gbx = mfi.growntilebox();
-
388  Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
-
389  Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
-
390  Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
-
391 
-
392  tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
-
393  tbxp2.grow(IntVect(NGROW,NGROW,0));
-
394  //copy the tilebox
-
395  //Box gbx1 = bx;
-
396  //Box gbx11 = bx;
-
397  //Box gbx2 = bx;
-
398  //TODO: adjust for tiling
-
399  //Box gbx3uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-3,bx.smallEnd(1)-3,bx.smallEnd(2))),
-
400  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+2,bx.bigEnd(1)+2,bx.bigEnd(2))));
-
401  //Box gbx2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-2,bx.smallEnd(1)-2,bx.smallEnd(2))),
-
402  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+1,bx.bigEnd(1)+1,bx.bigEnd(2))));
-
403  Box ubx = surroundingNodes(bx,0);
-
404  Box vbx = surroundingNodes(bx,1);
-
405  //make only gbx be grown to match multifabs
-
406  //gbx2.grow(IntVect(NGROW,NGROW,0));
-
407  //gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
-
408  //gbx11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
-
409 
-
410  Box bxD = bx;
-
411  bxD.makeSlab(2,0);
-
412  Box gbx1D = gbx1;
-
413  gbx1D.makeSlab(2,0);
-
414  Box gbx2D = gbx2;
-
415  gbx2D.makeSlab(2,0);
-
416 
-
417  Box tbxp1D = tbxp1;
-
418  tbxp1D.makeSlab(2,0);
-
419  Box tbxp2D = tbxp2;
-
420  tbxp2D.makeSlab(2,0);
-
421  //gbx1D.grow(IntVect(NGROW-1,NGROW-1,0));
-
422  //gbx2D.grow(IntVect(NGROW,NGROW,0));
-
423 
-
424  FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena()); //3D
-
425  FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena()); //3D
-
426  FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena()); //3D
-
427  FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
-
428  FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
-
429  FArrayBox fab_pn(tbxp2D,1,amrex::The_Async_Arena());
-
430  FArrayBox fab_pm(tbxp2D,1,amrex::The_Async_Arena());
-
431  FArrayBox fab_on_u(tbxp2D,1,amrex::The_Async_Arena());
-
432  FArrayBox fab_om_v(tbxp2D,1,amrex::The_Async_Arena());
-
433  FArrayBox fab_om_u(tbxp2D,1,amrex::The_Async_Arena());
-
434  FArrayBox fab_on_v(tbxp2D,1,amrex::The_Async_Arena());
-
435  FArrayBox fab_om_r(tbxp2D,1,amrex::The_Async_Arena());
-
436  FArrayBox fab_on_r(tbxp2D,1,amrex::The_Async_Arena());
-
437  FArrayBox fab_om_p(tbxp2D,1,amrex::The_Async_Arena());
-
438  FArrayBox fab_on_p(tbxp2D,1,amrex::The_Async_Arena());
-
439  FArrayBox fab_pmon_u(tbxp2D,1,amrex::The_Async_Arena());
-
440  FArrayBox fab_pnom_u(tbxp2D,1,amrex::The_Async_Arena());
-
441  FArrayBox fab_pmon_v(tbxp2D,1,amrex::The_Async_Arena());
-
442  FArrayBox fab_pnom_v(tbxp2D,1,amrex::The_Async_Arena());
-
443  FArrayBox fab_fomn(tbxp2D,1,amrex::The_Async_Arena());
-
444  //FArrayBox fab_oHz(gbx11,1,amrex::The_Async_Arena());
-
445 
-
446  auto FC=fab_FC.array();
-
447  auto FX=fab_FX.array();
-
448  auto FE=fab_FE.array();
-
449  auto pn=fab_pn.array();
-
450  auto pm=fab_pm.array();
-
451  auto on_u=fab_on_u.array();
-
452  auto om_v=fab_om_v.array();
-
453  auto om_u=fab_om_u.array();
-
454  auto on_v=fab_on_v.array();
-
455  auto om_r=fab_om_r.array();
-
456  auto on_r=fab_on_r.array();
-
457  auto om_p=fab_om_p.array();
-
458  auto on_p=fab_on_p.array();
-
459  auto pmon_u=fab_pmon_u.array();
-
460  auto pnom_u=fab_pnom_u.array();
-
461  auto pmon_v=fab_pmon_v.array();
-
462  auto pnom_v=fab_pnom_v.array();
-
463  auto fomn=fab_fomn.array();
-
464 
-
465  //From ana_grid.h and metrics.F
-
466  amrex::ParallelFor(tbxp2D,
-
467  [=] AMREX_GPU_DEVICE (int i, int j, int )
-
468  {
-
469  pm(i,j,0)=dxi[0];
-
470  pn(i,j,0)=dxi[1];
-
471  //defined UPWELLING
-
472  Real f0=-8.26e-5;
-
473  Real beta=0.0;
-
474  Real Esize=1000*(Mm);
-
475  Real y = prob_lo[1] + (j + 0.5) * dx[1];
-
476  Real f=f0+beta*(y-.5*Esize);
-
477  fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0)));
-
478  });
-
479 
-
480  amrex::ParallelFor(tbxp2D,
-
481  [=] AMREX_GPU_DEVICE (int i, int j, int )
-
482  {
-
483  //Note: are the comment definitons right? Don't seem to match metrics.f90
-
484  om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j))
-
485  on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j))
-
486  om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j)
-
487  on_r(i,j,0)=1.0/dxi[1]; // 1/pn(i,j)
-
488  //todo: om_p on_p
-
489  om_p(i,j,0)=1.0/dxi[0]; // 4/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
-
490  on_p(i,j,0)=1.0/dxi[1]; // 4/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
-
491  on_v(i,j,0)=1.0/dxi[1]; // 2/(pn(i-1,j)+pn(i,j))
-
492  om_u(i,j,0)=1.0/dxi[0]; // 2/(pm(i-1,j)+pm(i,j))
-
493  pmon_u(i,j,0)=1.0; // (pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
-
494  pnom_u(i,j,0)=1.0; // (pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
-
495  pmon_v(i,j,0)=1.0; // (pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j))
-
496  pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j))
-
497  });
-
498 
-
499  amrex::ParallelFor(gbx2,
-
500  [=] AMREX_GPU_DEVICE (int i, int j, int k)
-
501  {
-
502  FC(i,j,k)=0.0;
-
503  });
-
504 
-
505  prsgrd(gbx1,ru,rv,on_u,om_v,rho,FC,Hz,z_r,z_w,nrhs,N);
-
506  if (verbose > 2) {
-
507  amrex::PrintToFile("ru_afterprsgrd").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
-
508  amrex::PrintToFile("rv_afterprsgrd").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
-
509  }
-
510 
-
511  if (verbose > 2) {
-
512  amrex::PrintToFile("temp_pret3dmix").SetPrecision(18)<<FArrayBox(temp) << std::endl;
-
513  amrex::PrintToFile("salt_pret3dmix").SetPrecision(18)<<FArrayBox(salt) << std::endl;
-
514  }
-
515  t3dmix(bx, temp, diff2_temp, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev);
-
516  t3dmix(bx, salt, diff2_salt, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev);
-
517  if (verbose > 2) {
-
518  amrex::PrintToFile("temp_postt3dmix").SetPrecision(18)<<FArrayBox(temp) << std::endl;
-
519  amrex::PrintToFile("salt_postt3dmix").SetPrecision(18)<<FArrayBox(salt) << std::endl;
-
520  }
-
521 
-
522  if (solverChoice.use_coriolis) {
-
523  //-----------------------------------------------------------------------
-
524  // coriolis
-
525  //-----------------------------------------------------------------------
-
526  //
-
527  // ru, rv updated
-
528  // In ROMS, coriolis is the first (un-ifdefed) thing to happen in rhs3d_tile, which gets called after t3dmix
-
529  coriolis(bx, gbx, uold, vold, ru, rv, Hz, fomn, nrhs, nrhs);
-
530  }
-
531  if (verbose > 2) {
-
532  amrex::PrintToFile("ru_aftercor").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
-
533  amrex::PrintToFile("rv_aftercor").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
-
534  amrex::PrintToFile("u_aftercor").SetPrecision(18)<<FArrayBox(u)<<std::endl;
-
535  amrex::PrintToFile("v_aftercor").SetPrecision(18)<<FArrayBox(v)<<std::endl;
-
536  amrex::PrintToFile("temp_aftercor").SetPrecision(18)<<FArrayBox(temp)<<std::endl;
-
537  amrex::PrintToFile("tempstore_aftercor").SetPrecision(18)<<FArrayBox(tempstore)<<std::endl;
-
538  amrex::PrintToFile("salt_aftercor").SetPrecision(18)<<FArrayBox(salt)<<std::endl;
-
539  amrex::PrintToFile("saltstore_aftercor").SetPrecision(18)<<FArrayBox(saltstore)<<std::endl;
-
540  }
-
541 
-
542  //
-
543  //-----------------------------------------------------------------------
-
544  // rhs_3d
-
545  //-----------------------------------------------------------------------
-
546  //
-
547 
-
548  ////rufrc from 3d is set to ru, then the wind stress (and bottom stress) is added, then the mixing is added
-
549  //rufrc=ru+sustr*om_u*on_u
-
550  rhs_3d(bx, gbx, uold, vold, ru, rv, rufrc, rvfrc, sustr, svstr, bustr, bvstr, Huon, Hvom, on_u, om_v, om_u, on_v, W, FC, nrhs, N);
-
551  if (verbose > 2) {
-
552  amrex::PrintToFile("ru_afterrhs").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
-
553  amrex::PrintToFile("rv_afterrhs").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
-
554  amrex::PrintToFile("u_afterrhs").SetPrecision(18)<<FArrayBox(u)<<std::endl;
-
555  amrex::PrintToFile("v_afterrhs").SetPrecision(18)<<FArrayBox(v)<<std::endl;
-
556  }
-
557 
-
558  if(solverChoice.use_uv3dmix) {
-
559  //u=u+(contributions from S-surfaces viscosity not scaled by dt)*dt*dx*dy
-
560  //rufrc=rufrc + (contributions from S-surfaces viscosity not scaled by dt*dx*dy)
-
561  //For debugging version: uv3dmix(bx, u, v, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev,false);
-
562  uv3dmix(bx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev);
-
563  if (verbose > 2) {
-
564  amrex::PrintToFile("ru_afteruvmix").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
-
565  amrex::PrintToFile("rv_afteruvmix").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
-
566  amrex::PrintToFile("u_afteruvmix").SetPrecision(18)<<FArrayBox(u)<<std::endl;
-
567  amrex::PrintToFile("v_afteruvmix").SetPrecision(18)<<FArrayBox(v)<<std::endl;
-
568  amrex::PrintToFile("Huon").SetPrecision(18)<<FArrayBox(Huon)<<std::endl;
-
569  amrex::PrintToFile("temp_uv3dmix").SetPrecision(18)<<FArrayBox(temp)<<std::endl;
-
570  amrex::PrintToFile("tempstore_uv3dmix").SetPrecision(18)<<FArrayBox(tempstore)<<std::endl;
-
571  amrex::PrintToFile("salt_uv3dmix").SetPrecision(18)<<FArrayBox(salt)<<std::endl;
-
572  amrex::PrintToFile("saltstore_uv3dmix").SetPrecision(18)<<FArrayBox(saltstore)<<std::endl;
+
338  mf_W.FillBoundary(geom[lev].periodicity());
+
339 
+
340  for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
+
341  {
+
342  Array4<Real> const& DC = mf_DC.array(mfi);
+
343  Array4<Real> const& Akv = (vec_Akv[lev])->array(mfi);
+
344  Array4<Real> const& Hz = (vec_Hz[lev])->array(mfi);
+
345  Array4<Real> const& Huon = (vec_Huon[lev])->array(mfi);
+
346  Array4<Real> const& Hvom = (vec_Hvom[lev])->array(mfi);
+
347  Array4<Real> const& z_r = (mf_z_r)->array(mfi);
+
348  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
+
349  Array4<Real> const& uold = (mf_uold).array(mfi);
+
350  Array4<Real> const& vold = (mf_vold).array(mfi);
+
351  Array4<Real> const& u = (mf_u).array(mfi);
+
352  Array4<Real> const& v = (mf_v).array(mfi);
+
353  Array4<Real> const& pden = (mf_pden).array(mfi);
+
354  Array4<Real> const& rho = (mf_rho).array(mfi);
+
355  Array4<Real> const& rhoA = (mf_rhoA).array(mfi);
+
356  Array4<Real> const& rhoS = (mf_rhoS).array(mfi);
+
357  Array4<Real> const& tempold = (mf_tempold).array(mfi);
+
358  Array4<Real> const& saltold = (mf_saltold).array(mfi);
+
359  Array4<Real> const& temp = (mf_temp).array(mfi);
+
360  Array4<Real> const& salt = (mf_salt).array(mfi);
+
361  Array4<Real> const& tempstore = (vec_t3[lev])->array(mfi);
+
362  Array4<Real> const& saltstore = (vec_s3[lev])->array(mfi);
+
363  Array4<Real> const& ru = (mf_ru)->array(mfi);
+
364  Array4<Real> const& rv = (mf_rv)->array(mfi);
+
365  Array4<Real> const& rufrc = (mf_rufrc)->array(mfi);
+
366  Array4<Real> const& rvfrc = (mf_rvfrc)->array(mfi);
+
367  Array4<Real> const& W = (mf_W).array(mfi);
+
368  Array4<Real> const& sustr = (mf_sustr)->array(mfi);
+
369  Array4<Real> const& svstr = (mf_svstr)->array(mfi);
+
370  Array4<Real> const& rdrag = (mf_rdrag)->array(mfi);
+
371  Array4<Real> const& bustr = (mf_bustr)->array(mfi);
+
372  Array4<Real> const& bvstr = (mf_bvstr)->array(mfi);
+
373  Array4<Real> const& ubar = (mf_ubar)->array(mfi);
+
374  Array4<Real> const& vbar = (mf_vbar)->array(mfi);
+
375  Array4<Real> const& visc2_p = (mf_visc2_p)->array(mfi);
+
376  Array4<Real> const& visc2_r = (mf_visc2_r)->array(mfi);
+
377  Array4<Real> const& diff2_salt = (mf_diff2_salt)->array(mfi);
+
378  Array4<Real> const& diff2_temp = (mf_diff2_temp)->array(mfi);
+
379 
+
380  Array4<Real> const& zeta = (vec_zeta[lev])->array(mfi);
+
381  Array4<Real> const& Zt_avg1 = (vec_Zt_avg1[lev])->array(mfi);
+
382 
+
383  Box bx = mfi.tilebox();
+
384  Box tbxp1 = bx;
+
385  Box tbxp2 = bx;
+
386  Box gbx = mfi.growntilebox();
+
387  Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
+
388  Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
+
389  Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
+
390 
+
391  tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
+
392  tbxp2.grow(IntVect(NGROW,NGROW,0));
+
393  //copy the tilebox
+
394  //Box gbx1 = bx;
+
395  //Box gbx11 = bx;
+
396  //Box gbx2 = bx;
+
397  //TODO: adjust for tiling
+
398  //Box gbx3uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-3,bx.smallEnd(1)-3,bx.smallEnd(2))),
+
399  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+2,bx.bigEnd(1)+2,bx.bigEnd(2))));
+
400  //Box gbx2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-2,bx.smallEnd(1)-2,bx.smallEnd(2))),
+
401  // IntVect(AMREX_D_DECL(bx.bigEnd(0)+1,bx.bigEnd(1)+1,bx.bigEnd(2))));
+
402  Box ubx = surroundingNodes(bx,0);
+
403  Box vbx = surroundingNodes(bx,1);
+
404  //make only gbx be grown to match multifabs
+
405  //gbx2.grow(IntVect(NGROW,NGROW,0));
+
406  //gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
+
407  //gbx11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
+
408 
+
409  Box bxD = bx;
+
410  bxD.makeSlab(2,0);
+
411  Box gbx1D = gbx1;
+
412  gbx1D.makeSlab(2,0);
+
413  Box gbx2D = gbx2;
+
414  gbx2D.makeSlab(2,0);
+
415 
+
416  Box tbxp1D = tbxp1;
+
417  tbxp1D.makeSlab(2,0);
+
418  Box tbxp2D = tbxp2;
+
419  tbxp2D.makeSlab(2,0);
+
420  //gbx1D.grow(IntVect(NGROW-1,NGROW-1,0));
+
421  //gbx2D.grow(IntVect(NGROW,NGROW,0));
+
422 
+
423  FArrayBox fab_FC(gbx2,1,amrex::The_Async_Arena()); //3D
+
424  FArrayBox fab_FX(gbx2,1,amrex::The_Async_Arena()); //3D
+
425  FArrayBox fab_FE(gbx2,1,amrex::The_Async_Arena()); //3D
+
426  FArrayBox fab_BC(gbx2,1,amrex::The_Async_Arena());
+
427  FArrayBox fab_CF(gbx2,1,amrex::The_Async_Arena());
+
428  FArrayBox fab_pn(tbxp2D,1,amrex::The_Async_Arena());
+
429  FArrayBox fab_pm(tbxp2D,1,amrex::The_Async_Arena());
+
430  FArrayBox fab_on_u(tbxp2D,1,amrex::The_Async_Arena());
+
431  FArrayBox fab_om_v(tbxp2D,1,amrex::The_Async_Arena());
+
432  FArrayBox fab_om_u(tbxp2D,1,amrex::The_Async_Arena());
+
433  FArrayBox fab_on_v(tbxp2D,1,amrex::The_Async_Arena());
+
434  FArrayBox fab_om_r(tbxp2D,1,amrex::The_Async_Arena());
+
435  FArrayBox fab_on_r(tbxp2D,1,amrex::The_Async_Arena());
+
436  FArrayBox fab_om_p(tbxp2D,1,amrex::The_Async_Arena());
+
437  FArrayBox fab_on_p(tbxp2D,1,amrex::The_Async_Arena());
+
438  FArrayBox fab_pmon_u(tbxp2D,1,amrex::The_Async_Arena());
+
439  FArrayBox fab_pnom_u(tbxp2D,1,amrex::The_Async_Arena());
+
440  FArrayBox fab_pmon_v(tbxp2D,1,amrex::The_Async_Arena());
+
441  FArrayBox fab_pnom_v(tbxp2D,1,amrex::The_Async_Arena());
+
442  FArrayBox fab_fomn(tbxp2D,1,amrex::The_Async_Arena());
+
443  //FArrayBox fab_oHz(gbx11,1,amrex::The_Async_Arena());
+
444 
+
445  auto FC=fab_FC.array();
+
446  auto FX=fab_FX.array();
+
447  auto FE=fab_FE.array();
+
448  auto pn=fab_pn.array();
+
449  auto pm=fab_pm.array();
+
450  auto on_u=fab_on_u.array();
+
451  auto om_v=fab_om_v.array();
+
452  auto om_u=fab_om_u.array();
+
453  auto on_v=fab_on_v.array();
+
454  auto om_r=fab_om_r.array();
+
455  auto on_r=fab_on_r.array();
+
456  auto om_p=fab_om_p.array();
+
457  auto on_p=fab_on_p.array();
+
458  auto pmon_u=fab_pmon_u.array();
+
459  auto pnom_u=fab_pnom_u.array();
+
460  auto pmon_v=fab_pmon_v.array();
+
461  auto pnom_v=fab_pnom_v.array();
+
462  auto fomn=fab_fomn.array();
+
463 
+
464  //From ana_grid.h and metrics.F
+
465  amrex::ParallelFor(tbxp2D,
+
466  [=] AMREX_GPU_DEVICE (int i, int j, int )
+
467  {
+
468  pm(i,j,0)=dxi[0];
+
469  pn(i,j,0)=dxi[1];
+
470  //defined UPWELLING
+
471  Real f0=-8.26e-5;
+
472  Real beta=0.0;
+
473  Real Esize=1000*(Mm);
+
474  Real y = prob_lo[1] + (j + 0.5) * dx[1];
+
475  Real f=f0+beta*(y-.5*Esize);
+
476  fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0)));
+
477  });
+
478 
+
479  amrex::ParallelFor(tbxp2D,
+
480  [=] AMREX_GPU_DEVICE (int i, int j, int )
+
481  {
+
482  //Note: are the comment definitons right? Don't seem to match metrics.f90
+
483  om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j))
+
484  on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j))
+
485  om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j)
+
486  on_r(i,j,0)=1.0/dxi[1]; // 1/pn(i,j)
+
487  //todo: om_p on_p
+
488  om_p(i,j,0)=1.0/dxi[0]; // 4/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j))
+
489  on_p(i,j,0)=1.0/dxi[1]; // 4/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j))
+
490  on_v(i,j,0)=1.0/dxi[1]; // 2/(pn(i-1,j)+pn(i,j))
+
491  om_u(i,j,0)=1.0/dxi[0]; // 2/(pm(i-1,j)+pm(i,j))
+
492  pmon_u(i,j,0)=1.0; // (pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
+
493  pnom_u(i,j,0)=1.0; // (pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
+
494  pmon_v(i,j,0)=1.0; // (pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j))
+
495  pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j))
+
496  });
+
497 
+
498  amrex::ParallelFor(gbx2,
+
499  [=] AMREX_GPU_DEVICE (int i, int j, int k)
+
500  {
+
501  FC(i,j,k)=0.0;
+
502  });
+
503 
+
504  prsgrd(gbx1,ru,rv,on_u,om_v,rho,FC,Hz,z_r,z_w,nrhs,N);
+
505  if (verbose > 2) {
+
506  amrex::PrintToFile("ru_afterprsgrd").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
+
507  amrex::PrintToFile("rv_afterprsgrd").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
+
508  }
+
509 
+
510  if (verbose > 2) {
+
511  amrex::PrintToFile("temp_pret3dmix").SetPrecision(18)<<FArrayBox(temp) << std::endl;
+
512  amrex::PrintToFile("salt_pret3dmix").SetPrecision(18)<<FArrayBox(salt) << std::endl;
+
513  }
+
514  t3dmix(bx, temp, diff2_temp, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev);
+
515  t3dmix(bx, salt, diff2_salt, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev);
+
516  if (verbose > 2) {
+
517  amrex::PrintToFile("temp_postt3dmix").SetPrecision(18)<<FArrayBox(temp) << std::endl;
+
518  amrex::PrintToFile("salt_postt3dmix").SetPrecision(18)<<FArrayBox(salt) << std::endl;
+
519  }
+
520 
+
521  if (solverChoice.use_coriolis) {
+
522  //-----------------------------------------------------------------------
+
523  // coriolis
+
524  //-----------------------------------------------------------------------
+
525  //
+
526  // ru, rv updated
+
527  // In ROMS, coriolis is the first (un-ifdefed) thing to happen in rhs3d_tile, which gets called after t3dmix
+
528  coriolis(bx, gbx, uold, vold, ru, rv, Hz, fomn, nrhs, nrhs);
+
529  }
+
530  if (verbose > 2) {
+
531  amrex::PrintToFile("ru_aftercor").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
+
532  amrex::PrintToFile("rv_aftercor").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
+
533  amrex::PrintToFile("u_aftercor").SetPrecision(18)<<FArrayBox(u)<<std::endl;
+
534  amrex::PrintToFile("v_aftercor").SetPrecision(18)<<FArrayBox(v)<<std::endl;
+
535  amrex::PrintToFile("temp_aftercor").SetPrecision(18)<<FArrayBox(temp)<<std::endl;
+
536  amrex::PrintToFile("tempstore_aftercor").SetPrecision(18)<<FArrayBox(tempstore)<<std::endl;
+
537  amrex::PrintToFile("salt_aftercor").SetPrecision(18)<<FArrayBox(salt)<<std::endl;
+
538  amrex::PrintToFile("saltstore_aftercor").SetPrecision(18)<<FArrayBox(saltstore)<<std::endl;
+
539  }
+
540 
+
541  //
+
542  //-----------------------------------------------------------------------
+
543  // rhs_3d
+
544  //-----------------------------------------------------------------------
+
545  //
+
546 
+
547  ////rufrc from 3d is set to ru, then the wind stress (and bottom stress) is added, then the mixing is added
+
548  //rufrc=ru+sustr*om_u*on_u
+
549  rhs_3d(bx, gbx, uold, vold, ru, rv, rufrc, rvfrc, sustr, svstr, bustr, bvstr, Huon, Hvom, on_u, om_v, om_u, on_v, W, FC, nrhs, N);
+
550  if (verbose > 2) {
+
551  amrex::PrintToFile("ru_afterrhs").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
+
552  amrex::PrintToFile("rv_afterrhs").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
+
553  amrex::PrintToFile("u_afterrhs").SetPrecision(18)<<FArrayBox(u)<<std::endl;
+
554  amrex::PrintToFile("v_afterrhs").SetPrecision(18)<<FArrayBox(v)<<std::endl;
+
555  }
+
556 
+
557  if(solverChoice.use_uv3dmix) {
+
558  //u=u+(contributions from S-surfaces viscosity not scaled by dt)*dt*dx*dy
+
559  //rufrc=rufrc + (contributions from S-surfaces viscosity not scaled by dt*dx*dy)
+
560  //For debugging version: uv3dmix(bx, u, v, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev,false);
+
561  uv3dmix(bx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev);
+
562  if (verbose > 2) {
+
563  amrex::PrintToFile("ru_afteruvmix").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
+
564  amrex::PrintToFile("rv_afteruvmix").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
+
565  amrex::PrintToFile("u_afteruvmix").SetPrecision(18)<<FArrayBox(u)<<std::endl;
+
566  amrex::PrintToFile("v_afteruvmix").SetPrecision(18)<<FArrayBox(v)<<std::endl;
+
567  amrex::PrintToFile("Huon").SetPrecision(18)<<FArrayBox(Huon)<<std::endl;
+
568  amrex::PrintToFile("temp_uv3dmix").SetPrecision(18)<<FArrayBox(temp)<<std::endl;
+
569  amrex::PrintToFile("tempstore_uv3dmix").SetPrecision(18)<<FArrayBox(tempstore)<<std::endl;
+
570  amrex::PrintToFile("salt_uv3dmix").SetPrecision(18)<<FArrayBox(salt)<<std::endl;
+
571  amrex::PrintToFile("saltstore_uv3dmix").SetPrecision(18)<<FArrayBox(saltstore)<<std::endl;
+
572  }
573  }
-
574  }
-
575 
-
576  // Set zeta components to time-averaged value in preparation for barotropic update
-
577  amrex::ParallelFor(gbx2D,
-
578  [=] AMREX_GPU_DEVICE (int i, int j, int)
-
579  {
-
580  zeta(i,j,0,0) = Zt_avg1(i,j,0);
-
581  zeta(i,j,0,1) = Zt_avg1(i,j,0);
-
582  });
-
583  } // MFIter
-
584 
-
585  // Update Akv with new depth. NOTE: this happens before set_zeta in ROMS
-
586  set_vmix(lev);
-
587 
-
588  mf_temp.FillBoundary(geom[lev].periodicity());
-
589  mf_salt.FillBoundary(geom[lev].periodicity());
-
590  mf_tempold.FillBoundary(geom[lev].periodicity());
-
591  mf_saltold.FillBoundary(geom[lev].periodicity());
-
592  vec_t3[lev]->FillBoundary(geom[lev].periodicity());
-
593  vec_s3[lev]->FillBoundary(geom[lev].periodicity());
-
594  vec_Huon[lev]->FillBoundary(geom[lev].periodicity());
-
595  vec_Hvom[lev]->FillBoundary(geom[lev].periodicity());
-
596 
-
597  if(solverChoice.use_barotropic) {
-
598  bool predictor_2d_step=true;
-
599  bool first_2d_step=(iic==ntfirst);
-
600  int nfast_counter=predictor_2d_step ? nfast : nfast-1;
-
601  nfast_counter += 1;
-
602  //Compute fast timestep from dt_lev and ratio
-
603  Real dtfast_lev=dt_lev/Real(fixed_ndtfast_ratio);
-
604  int next_indx1 = 0;
-
605  for(int my_iif = 0; my_iif < nfast_counter; my_iif++) {
-
606  first_2d_step=(my_iif==0);
-
607  //Predictor
-
608  predictor_2d_step=true;
-
609  advance_2d(lev, mf_u, mf_v, mf_rhoS, mf_rhoA, vec_ru[lev], vec_rv[lev],
-
610  vec_rufrc[lev], vec_rvfrc[lev],
-
611  vec_Zt_avg1[lev],
-
612  vec_DU_avg1[lev], vec_DU_avg2[lev],
-
613  vec_DV_avg1[lev], vec_DV_avg2[lev],
-
614  vec_rubar[lev], vec_rvbar[lev], vec_rzeta[lev],
-
615  vec_ubar[lev], vec_vbar[lev], vec_zeta[lev],
-
616  vec_hOfTheConfusingName[lev], vec_visc2_p[lev], vec_visc2_r[lev],
-
617  ncomp, dt_lev, dtfast_lev, predictor_2d_step, first_2d_step, my_iif, nfast, next_indx1);
-
618 
-
619  //Corrector. Skip it on last fast step
-
620  predictor_2d_step=false;
-
621  if (my_iif < nfast_counter - 1) {
-
622  advance_2d(lev, mf_u, mf_v, mf_rhoS, mf_rhoA, vec_ru[lev], vec_rv[lev],
-
623  vec_rufrc[lev], vec_rvfrc[lev],
-
624  vec_Zt_avg1[lev],
-
625  vec_DU_avg1[lev], vec_DU_avg2[lev],
-
626  vec_DV_avg1[lev], vec_DV_avg2[lev],
-
627  vec_rubar[lev], vec_rvbar[lev], vec_rzeta[lev],
-
628  vec_ubar[lev], vec_vbar[lev], vec_zeta[lev],
-
629  vec_hOfTheConfusingName[lev], vec_visc2_p[lev], vec_visc2_r[lev],
-
630  ncomp, dt_lev, dtfast_lev, predictor_2d_step, first_2d_step, my_iif, nfast, next_indx1);
-
631  }
+
574 
+
575  // Set zeta components to time-averaged value in preparation for barotropic update
+
576  amrex::ParallelFor(gbx2D,
+
577  [=] AMREX_GPU_DEVICE (int i, int j, int)
+
578  {
+
579  zeta(i,j,0,0) = Zt_avg1(i,j,0);
+
580  zeta(i,j,0,1) = Zt_avg1(i,j,0);
+
581  });
+
582  } // MFIter
+
583 
+
584  // Update Akv with new depth. NOTE: this happens before set_zeta in ROMS
+
585  set_vmix(lev);
+
586 
+
587  mf_temp.FillBoundary(geom[lev].periodicity());
+
588  mf_salt.FillBoundary(geom[lev].periodicity());
+
589  mf_tempold.FillBoundary(geom[lev].periodicity());
+
590  mf_saltold.FillBoundary(geom[lev].periodicity());
+
591  vec_t3[lev]->FillBoundary(geom[lev].periodicity());
+
592  vec_s3[lev]->FillBoundary(geom[lev].periodicity());
+
593  vec_Huon[lev]->FillBoundary(geom[lev].periodicity());
+
594  vec_Hvom[lev]->FillBoundary(geom[lev].periodicity());
+
595 
+
596  if(solverChoice.use_barotropic) {
+
597  bool predictor_2d_step=true;
+
598  bool first_2d_step=(iic==ntfirst);
+
599  int nfast_counter=predictor_2d_step ? nfast : nfast-1;
+
600  nfast_counter += 1;
+
601  //Compute fast timestep from dt_lev and ratio
+
602  Real dtfast_lev=dt_lev/Real(fixed_ndtfast_ratio);
+
603  int next_indx1 = 0;
+
604  for(int my_iif = 0; my_iif < nfast_counter; my_iif++) {
+
605  first_2d_step=(my_iif==0);
+
606  //Predictor
+
607  predictor_2d_step=true;
+
608  advance_2d(lev, mf_u, mf_v, mf_rhoS, mf_rhoA, vec_ru[lev], vec_rv[lev],
+
609  vec_rufrc[lev], vec_rvfrc[lev],
+
610  vec_Zt_avg1[lev],
+
611  vec_DU_avg1[lev], vec_DU_avg2[lev],
+
612  vec_DV_avg1[lev], vec_DV_avg2[lev],
+
613  vec_rubar[lev], vec_rvbar[lev], vec_rzeta[lev],
+
614  vec_ubar[lev], vec_vbar[lev], vec_zeta[lev],
+
615  vec_hOfTheConfusingName[lev], vec_visc2_p[lev], vec_visc2_r[lev],
+
616  ncomp, dt_lev, dtfast_lev, predictor_2d_step, first_2d_step, my_iif, nfast, next_indx1);
+
617 
+
618  //Corrector. Skip it on last fast step
+
619  predictor_2d_step=false;
+
620  if (my_iif < nfast_counter - 1) {
+
621  advance_2d(lev, mf_u, mf_v, mf_rhoS, mf_rhoA, vec_ru[lev], vec_rv[lev],
+
622  vec_rufrc[lev], vec_rvfrc[lev],
+
623  vec_Zt_avg1[lev],
+
624  vec_DU_avg1[lev], vec_DU_avg2[lev],
+
625  vec_DV_avg1[lev], vec_DV_avg2[lev],
+
626  vec_rubar[lev], vec_rvbar[lev], vec_rzeta[lev],
+
627  vec_ubar[lev], vec_vbar[lev], vec_zeta[lev],
+
628  vec_hOfTheConfusingName[lev], vec_visc2_p[lev], vec_visc2_r[lev],
+
629  ncomp, dt_lev, dtfast_lev, predictor_2d_step, first_2d_step, my_iif, nfast, next_indx1);
+
630  }
+
631  }
632  }
-
633  }
-
634 
-
635  // Currently, bathymetry doesn't change, so we do not need to re-set h here. Just re-do stretch, etc
-
636  // because zeta may have changed
-
637  stretch_transform(lev);
-
638 
-
639  advance_3d(lev, mf_u, mf_v, mf_tempold, mf_saltold, mf_temp, mf_salt, vec_t3[lev], vec_s3[lev], vec_ru[lev], vec_rv[lev],
-
640  vec_DU_avg1[lev], vec_DU_avg2[lev],
-
641  vec_DV_avg1[lev], vec_DV_avg2[lev],
-
642  vec_ubar[lev], vec_vbar[lev],
-
643  mf_AK, mf_DC,
-
644  mf_Hzk, vec_Akv[lev], vec_Akt[lev], vec_Hz[lev], vec_Huon[lev], vec_Hvom[lev], vec_z_w[lev], vec_hOfTheConfusingName[lev], ncomp, N, dt_lev);
-
645 
-
646  U_new.FillBoundary(geom[lev].periodicity());
-
647  V_new.FillBoundary(geom[lev].periodicity());
-
648 
-
649  U_old.FillBoundary(geom[lev].periodicity());
-
650  V_old.FillBoundary(geom[lev].periodicity());
-
651 
-
652  mf_temp.FillBoundary(geom[lev].periodicity());
-
653  mf_salt.FillBoundary(geom[lev].periodicity());
-
654 
-
655  mf_tempold.FillBoundary(geom[lev].periodicity());
-
656  mf_saltold.FillBoundary(geom[lev].periodicity());
-
657 
-
658  vec_t3[lev]->FillBoundary(geom[lev].periodicity());
-
659  vec_s3[lev]->FillBoundary(geom[lev].periodicity());
-
660 
-
661  for (int lev = 0; lev <= finest_level; ++lev) {
-
662  FillPatch(lev, t_new[lev], vars_new[lev]);
-
663  }
-
664 
-
665 #ifdef ROMSX_USE_PARTICLES
-
666  // Update tracer particles on level 0
-
667  if (lev == 0 && use_tracer_particles) {
-
668  MultiFab* Umac = &vars_new[lev][Vars::xvel];
-
669  tracer_particles->AdvectWithUmac(Umac, lev, dt_lev, *vec_z_phys_nd[0]);
-
670  }
-
671 #endif
+
633 
+
634  // Currently, bathymetry doesn't change, so we do not need to re-set h here. Just re-do stretch, etc
+
635  // because zeta may have changed
+
636  stretch_transform(lev);
+
637 
+
638  advance_3d(lev, mf_u, mf_v, mf_tempold, mf_saltold, mf_temp, mf_salt, vec_t3[lev], vec_s3[lev], vec_ru[lev], vec_rv[lev],
+
639  vec_DU_avg1[lev], vec_DU_avg2[lev],
+
640  vec_DV_avg1[lev], vec_DV_avg2[lev],
+
641  vec_ubar[lev], vec_vbar[lev],
+
642  mf_AK, mf_DC,
+
643  mf_Hzk, vec_Akv[lev], vec_Akt[lev], vec_Hz[lev], vec_Huon[lev], vec_Hvom[lev], vec_z_w[lev], vec_hOfTheConfusingName[lev], ncomp, N, dt_lev);
+
644 
+
645  U_new.FillBoundary(geom[lev].periodicity());
+
646  V_new.FillBoundary(geom[lev].periodicity());
+
647 
+
648  U_old.FillBoundary(geom[lev].periodicity());
+
649  V_old.FillBoundary(geom[lev].periodicity());
+
650 
+
651  mf_temp.FillBoundary(geom[lev].periodicity());
+
652  mf_salt.FillBoundary(geom[lev].periodicity());
+
653 
+
654  mf_tempold.FillBoundary(geom[lev].periodicity());
+
655  mf_saltold.FillBoundary(geom[lev].periodicity());
+
656 
+
657  vec_t3[lev]->FillBoundary(geom[lev].periodicity());
+
658  vec_s3[lev]->FillBoundary(geom[lev].periodicity());
+
659 
+
660  for (int lev = 0; lev <= finest_level; ++lev) {
+
661  FillPatch(lev, t_new[lev], vars_new[lev]);
+
662  }
+
663 
+
664 #ifdef ROMSX_USE_PARTICLES
+
665  // Update tracer particles on level 0
+
666  if (lev == 0 && use_tracer_particles) {
+
667  MultiFab* Umac = &vars_new[lev][Vars::xvel];
+
668  tracer_particles->AdvectWithUmac(Umac, lev, dt_lev, *vec_z_phys_nd[0]);
+
669  }
+
670 #endif
+
671 
672 
-
673 
-
674  //We are not storing computed W aka Omega
-
675  // MultiFab::Copy(W_new,mf_w,0,0,W_new.nComp(),IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
-
676  // W_new.FillBoundary();
-
677  // MultiFab::Copy(mf_W,S_old,Omega_comp,0,mf_W.nComp(),mf_w.nGrowVect());
-
678 
-
679 // set_drag(lev);
-
680 
-
681 }
+
673  //We are not storing computed W aka Omega
+
674  // MultiFab::Copy(W_new,mf_w,0,0,W_new.nComp(),IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));
+
675  // W_new.FillBoundary();
+
676  // MultiFab::Copy(mf_W,S_old,Omega_comp,0,mf_W.nComp(),mf_w.nGrowVect());
+
677 
+
678 // set_drag(lev);
+
679 
+
680 }