ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_ComputeDiffusivityMYJ.cpp File Reference
#include "ERF_SurfaceLayer.H"
#include "ERF_DirectionSelector.H"
#include "ERF_Diffusion.H"
#include "ERF_Constants.H"
#include "ERF_TurbStruct.H"
#include "ERF_PBLModels.H"
#include <math.h>
Include dependency graph for ERF_ComputeDiffusivityMYJ.cpp:

Functions

void ComputeDiffusivityMYJ (Real dt, const MultiFab &xvel, const MultiFab &yvel, MultiFab &cons_in, MultiFab &eddyViscosity, const Geometry &geom, const TurbChoice &, std::unique_ptr< SurfaceLayer > &, bool use_terrain_fitted_coords, bool, int, const BCRec *bc_ptr, bool, const std::unique_ptr< MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
 

Function Documentation

◆ ComputeDiffusivityMYJ()

void ComputeDiffusivityMYJ ( Real  dt,
const MultiFab &  xvel,
const MultiFab &  yvel,
MultiFab &  cons_in,
MultiFab &  eddyViscosity,
const Geometry &  geom,
const TurbChoice &  ,
std::unique_ptr< SurfaceLayer > &  ,
bool  use_terrain_fitted_coords,
bool  ,
int  ,
const BCRec *  bc_ptr,
bool  ,
const std::unique_ptr< MultiFab > &  z_phys_nd,
const MoistureComponentIndices &  moisture_indices 
)
28 {
29  // Dirichlet flags to switch derivative stencil
30  bool c_ext_dir_on_zlo = ( (bc_ptr[BCVars::cons_bc].lo(2) == ERFBCType::ext_dir) );
31  bool c_ext_dir_on_zhi = ( (bc_ptr[BCVars::cons_bc].lo(5) == ERFBCType::ext_dir) );
32  bool u_ext_dir_on_zlo = ( (bc_ptr[BCVars::xvel_bc].lo(2) == ERFBCType::ext_dir) );
33  bool u_ext_dir_on_zhi = ( (bc_ptr[BCVars::xvel_bc].lo(5) == ERFBCType::ext_dir) );
34  bool v_ext_dir_on_zlo = ( (bc_ptr[BCVars::yvel_bc].lo(2) == ERFBCType::ext_dir) );
35  bool v_ext_dir_on_zhi = ( (bc_ptr[BCVars::yvel_bc].lo(5) == ERFBCType::ext_dir) );
36 
37  // Expose constants
38  Real d_kappa = KAPPA;
39 
40  // Closure coefficients (from Janjic (2002), NCEP Office Note 437)
41  Real EPS1 = 1.0e-12;
42  Real EPS2 = 0.;
43  Real EPSRS = 1.0e-7;
44  Real EPSRU = 1.0e-7;
45  Real EPSTRB = 1.0e-24;
46  Real EPSL = 0.32;
47  Real EPSQ2 = 0.2;
48  Real EPSQ1 = std::sqrt(EPSQ2);
49  Real ESQHF = 2.5;
50  Real FH = 1.01;
51 
52  Real G = CONST_GRAV;
53  Real ALPHA = 0.3;
54  Real BETA = 1./273.;
55  Real EL0MAX = 1000.;
56  Real EL0MIN = 1.;
57  Real ELFC = 0.23*0.5;
58 
59  Real A1 = 0.659888514560862645;
60  Real A2x = 0.6574209922667784586;
61  Real B1 = 11.87799326209552761;
62  Real B2 = 7.226971804046074028;
63  Real C1 = 0.000830955950095854396;
64 
65  Real BTG = BETA*G;
66  Real RB1 = 1./B1;
67 
68  Real ADNH = 9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG;
69  Real ADNM = 18.*A1*A1*A2x*(B2-3.*A2x)*BTG;
70  Real ANMH = -9.*A1*A2x*A2x*BTG*BTG;
71  Real ANMM = -3.*A1*A2x*(3.*A2x+3.*B2*C1+18.*A1*C1-B2)*BTG;
72  Real BDNH = 3.*A2x*(7.*A1+B2)*BTG;
73  Real BDNM = 6.*A1*A1;
74  Real BEQH = A2x*B1*BTG+3.*A2x*(7.*A1+B2)*BTG;
75  Real BEQM = -A1*B1*(1.-3.*C1)+6.*A1*A1;
76  Real BNMH = -A2x*BTG;
77  Real BNMM = A1*(1.-3.*C1);
78  Real BSHH = 9.*A1*A2x*A2x*BTG;
79  Real BSHM = 18.*A1*A1*A2x*C1;
80  Real BSMH = -3.*A1*A2x*(3.*A2x+3.*B2*C1+12.*A1*C1-B2)*BTG;
81  Real CESH = A2x;
82  Real CESM = A1*(1.-3.*C1);
83 
84  Real AEQH = 9.*A1*A2x*A2x*B1*BTG*BTG
85  + 9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG;
86  Real AEQM = 3.*A1*A2x*B1*(3.*A2x+3.*B2*C1+18.*A1*C1-B2)*BTG
87  + 18.*A1*A1*A2x*(B2-3.*A2x)*BTG;
88 
89  Real REQU = -AEQH/AEQM;
90  Real EPSGH = 1.E-9;
91  Real EPSGM = REQU*EPSGH;
92 
93  Real UBRYL = (18.*REQU*A1*A1*A2x*B2*C1*BTG + 9.*A1*A2x*A2x*B2*BTG*BTG)
94  / (REQU*ADNM+ADNH);
95  Real UBRY = (1.+EPSRS)*UBRYL;
96  Real UBRY3 = 3.*UBRY;
97 
98  Real AUBH = 27.*A1*A2x*A2x*B2*BTG*BTG-ADNH*UBRY3;
99  Real AUBM = 54.*A1*A1*A2x*B2*C1*BTG-ADNM*UBRY3;
100  Real BUBH = (9.*A1*A2x+3.*A2x*B2)*BTG-BDNH*UBRY3;
101  Real BUBM = 18.*A1*A1*C1-BDNM*UBRY3;
102  Real CUBR = 1.-UBRY3;
103  Real RCUBR = 1./CUBR;
104 
105 #ifdef _OPENMP
106 #pragma omp parallel if (Gpu::notInLaunchRegion())
107 #endif
108  for (MFIter mfi(eddyViscosity,false); mfi.isValid(); ++mfi) {
109 
110  const Box& bx = mfi.validbox();
111  const Array4<Real >& cell_data = cons_in.array(mfi);
112  const Array4<Real >& K_turb = eddyViscosity.array(mfi);
113  const Array4<Real const>& uvel = xvel.array(mfi);
114  const Array4<Real const>& vvel = yvel.array(mfi);
115 
116  // Ensure the box spans the vertical domain
117  const Box& dbx = geom.Domain();
118  AMREX_ALWAYS_ASSERT(bx.smallEnd(2) == dbx.smallEnd(2) && bx.bigEnd(2) == dbx.bigEnd(2));
119 
120  // Create a plane box
121  int klo = bx.smallEnd(2);
122  int khi = bx.bigEnd(2);
123  Box planexy = makeSlab(bx,2,klo);
124 
125  // Expose for GPU capture
126  const GeometryData gdata = geom.data();
127 
128  // Allocate space for integrals
129  const Box xybx = PerpendicularBox<ZDir>(bx, IntVect{0,0,0});
130  FArrayBox qturb(bx,1);
131  FArrayBox qintegral(xybx,2);
132  IArrayBox pbl_k(xybx,1);
133  qintegral.setVal<RunOn::Device>(0.0);
134  pbl_k.setVal<RunOn::Device>(khi);
135  const Array4<Real> qint = qintegral.array();
136  const Array4<Real> qvel = qturb.array();
137  const Array4<int> k_arr = pbl_k.array();
138 
139  // Terrain and gradient calcs
140  const Array4<Real const> &z_nd_arr = z_phys_nd->array(mfi);
141  const auto& dxInv = geom.InvCellSizeArray();
142  int izmin = geom.Domain().smallEnd(2);
143  int izmax = geom.Domain().bigEnd(2);
144 
145  // Ustar for BC
146  //MultiFab* ustar = SurfLayer->get_u_star(level);
147  //const Array4<Real const>& ustar_arr = ustar->array(mfi);
148 
149  // Vertical integrals to compute l0
150  if (use_terrain_fitted_coords) {
151  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
152  {
153  // Locate PBL k index and set qvel
154  for (int k(klo); k<=khi; ++k) {
155  Real q2 = 2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
156  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(q2 > 0.0, "KE must have a positive value");
157  qvel(i,j,k) = std::sqrt(q2);
158  if (q2<=EPSQ2*FH) {
159  k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
160  }
161  }
162 
163  // Perform integral over PBL height
164  for (int k(klo); k<=k_arr(i,j,0); ++k) {
165  const Real dz = Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr);
166  const Real Zval = Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr);
167  Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k)*dz);
168  Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k)*dz);
169  }
170  });
171  } else {
172  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
173  {
174  // Locate PBL k index and set qvel
175  for (int k(klo); k<=khi; ++k) {
176  Real q2 = 2.0 * cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
177  AMREX_ALWAYS_ASSERT_WITH_MESSAGE(q2 > 0.0, "KE must have a positive value");
178  qvel(i,j,k) = std::sqrt(q2);
179  if (q2<=EPSQ2*FH) {
180  k_arr(i,j,0) = std::min(k,k_arr(i,j,0));
181  }
182  }
183 
184  // Perform integral over PBL height
185  for (int k(klo); k<=k_arr(i,j,0); ++k) {
186  // Not multiplying by dz: it's constant and would fall out when we divide qint0/qint1 anyway
187  const Real Zval = gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
188  Gpu::Atomic::Add(&qint(i,j,0,0), Zval*qvel(i,j,k));
189  Gpu::Atomic::Add(&qint(i,j,0,1), qvel(i,j,k));
190  }
191  });
192  }
193 
194  // Main work to fill diffusivities
195  ParallelFor(planexy, [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/) noexcept
196  {
197  // Get the PBL k index
198  int kpbl = k_arr(i,j,0);
199 
200  // Compute the integral length scale
201  Real l0 = std::max(std::min(ALPHA*qint(i,j,0,0)/qint(i,j,0,1),EL0MAX),EL0MIN);
202 
203  // Compute diffusivities in each column
204  for (int k(klo); k<=khi; ++k) {
205  // Gradients for shear and buoy production
206  const Real met_h_zeta = use_terrain_fitted_coords ? Compute_h_zeta_AtCellCenter(i,j,k,dxInv,z_nd_arr) : 1.0;
207  Real dthetavdz, dudz, dvdz;
208  ComputeVerticalDerivativesPBL(i, j, k,
209  uvel, vvel, cell_data, izmin, izmax, dxInv[2]/met_h_zeta,
210  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
211  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
212  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
213  dthetavdz, dudz, dvdz,
214  moisture_indices);
215 
216  // Calculate dimensional production terms
217  Real GML = std::max(dudz*dudz + dvdz*dvdz, EPSGM);
218  // NOTE: model uses BTG = beta*g in coeffs above
219  // NOTE: sign convention follows code but theory differs
220  Real GHL = dthetavdz;
221  if (std::fabs(GHL)<=EPSGH) { GHL=EPSGH; }
222 
223  // Find the maximum mixing length
224  Real ELM;
225  if (GHL >= EPSGH) {
226  if (GML/GHL <= REQU) {
227  ELM = EPSL;
228  } else {
229  Real AUBR = (AUBM*GML+AUBH*GHL)*GHL;
230  Real BUBR = BUBM*GML+BUBH*GHL;
231  Real QOL2ST = (-0.5*BUBR+std::sqrt(BUBR*BUBR*0.25-AUBR*CUBR))*RCUBR;
232  Real ELOQ2X = 1./QOL2ST;
233  ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
234  }
235  } else {
236  Real ADEN = (ADNM*GML+ADNH*GHL)*GHL;
237  Real BDEN = BDNM*GML+BDNH*GHL;
238  Real QOL2UN = -0.5*BDEN+std::sqrt(BDEN*BDEN*0.25-ADEN);
239  Real ELOQ2X = 1./(QOL2UN+EPSRU);
240  ELM = std::max(std::sqrt(ELOQ2X*qvel(i,j,k)*qvel(i,j,k)),EPSL);
241  }
242 
243  // Compute master length scale
244  Real L;
245  if (k>kpbl) {
246  L = std::min((met_h_zeta/dxInv[2])*ELFC, ELM);
247  } else {
248  const Real zval = use_terrain_fitted_coords ? Compute_Zrel_AtCellCenter(i,j,k,z_nd_arr)
249  : gdata.ProbLo(2) + (k + 0.5)*gdata.CellSize(2);
250  L = std::min(l0*d_kappa*zval / (d_kappa*zval + l0), ELM);
251  }
252 
253  // Update qvel from production and dissipation
254  Real AEQU = (AEQM*GML+AEQH*GHL)*GHL;
255  Real BEQU = BEQM*GML+BEQH*GHL;
256 
257  Real EQOL2 = -0.5*BEQU+std::sqrt(BEQU*BEQU*0.25-AEQU);
258 
259  if ( ((GML+GHL*GHL)<=EPSTRB) ||
260  ((GHL>=EPSGH) && ((GML/GHL)<=REQU)) ||
261  (EQOL2<=EPS2) ) {
262  L = EPSL;
263  qvel(i,j,k) = EPSQ1;
264  } else {
265  Real ANUM=(ANMM*GML+ANMH*GHL)*GHL;
266  Real BNUM= BNMM*GML+BNMH*GHL;
267 
268  Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
269  Real BDEN= BDNM*GML+BDNH*GHL;
270  Real CDEN= 1.;
271 
272  Real ARHS=-(ANUM*BDEN-BNUM*ADEN)*2.;
273  Real BRHS=- ANUM*4.;
274  Real CRHS=- BNUM*2.;
275 
276  Real DLOQ1=L/qvel(i,j,k);
277 
278  Real ELOQ21=1./EQOL2;
279  Real ELOQ11=std::sqrt(ELOQ21);
280  Real ELOQ31=ELOQ21*ELOQ11;
281  Real ELOQ41=ELOQ21*ELOQ21;
282  Real ELOQ51=ELOQ21*ELOQ31;
283 
284  Real RDEN1=1./(ADEN*ELOQ41+BDEN*ELOQ21+CDEN);
285 
286  Real RHSP1=(ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1;
287 
288  Real DTTURBL=dt;
289  Real ELOQ12=std::max(ELOQ11+(DLOQ1-ELOQ11)*exp(RHSP1*DTTURBL),EPS1);
290 
291  Real ELOQ22=ELOQ12*ELOQ12;
292  Real ELOQ32=ELOQ22*ELOQ12;
293  Real ELOQ42=ELOQ22*ELOQ22;
294  Real ELOQ52=ELOQ22*ELOQ32;
295 
296  Real RDEN2=1./(ADEN*ELOQ42+BDEN*ELOQ22+CDEN);
297  Real RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1;
298  Real RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2;
299  Real RHST2=RHS2/RHSP2;
300 
301  Real ELOQ13=std::max(ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*exp(RHSP2*DTTURBL),EPS1);
302 
303  Real ELOQN=ELOQ13;
304  if (ELOQN>EPS1) {
305  qvel(i,j,k) = std::max(L/ELOQN,EPSQ1);
306  if (qvel(i,j,k)==EPSQ1) {L = EPSL; }
307  } else {
308  L = EPSL;
309  qvel(i,j,k) = EPSQ1;
310  }
311  }
312  /*
313  // Boundary condition
314  if (k==klo) {
315  Real q2 = std::pow(B1,(2./3.))*ustar_arr(i,j,k)*ustar_arr(i,j,k);
316  Real q = std::max(std::sqrt(q2),EPSQ1);
317  qvel(i,j,k) = 0.5 * (q + qvel(i,j,k));
318  }
319  */
320  cell_data(i,j,k,RhoKE_comp) = 0.5*cell_data(i,j,k,Rho_comp)*qvel(i,j,k)*qvel(i,j,k);
321 
322  // L^n/Q^n
323  Real ELOQ2 = L*L/(qvel(i,j,k)*qvel(i,j,k));
324  Real ELOQ4 = ELOQ2*ELOQ2;
325 
326  // COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
327  Real ADEN=(ADNM*GML+ADNH*GHL)*GHL;
328  Real BDEN= BDNM*GML+BDNH*GHL;
329  Real CDEN= 1.;
330 
331  // COEFFICIENTS FOR THE SM DETERMINANT
332  Real BESM=BSMH*GHL;
333 
334  // COEFFICIENTS FOR THE SH DETERMINANT
335  Real BESH=BSHM*GML+BSHH*GHL;
336 
337  // 1./DENOMINATOR
338  Real RDEN=1./(ADEN*ELOQ4+BDEN*ELOQ2+CDEN);
339 
340  // SM, SH, SQ
341  Real SM=(BESM*ELOQ2+CESM)*RDEN;
342  Real SH=(BESH*ELOQ2+CESH)*RDEN;
343  Real SQ=ESQHF*SH;
344 
345  // Finally, compute the eddy viscosity/diffusivities
346  const Real rho = cell_data(i,j,k,Rho_comp);
347  K_turb(i,j,k,EddyDiff::Mom_v ) = rho * L * qvel(i,j,k) * SM;
348  K_turb(i,j,k,EddyDiff::Theta_v) = rho * L * qvel(i,j,k) * SH;
349  K_turb(i,j,k,EddyDiff::KE_v ) = rho * L * qvel(i,j,k) * SQ;
350  K_turb(i,j,k,EddyDiff::Q_v ) = rho * L * qvel(i,j,k) * SH;
351  K_turb(i,j,k,EddyDiff::Turb_lengthscale) = L;
352 
353  // NOTE: Ghost cells are handled at end of ERF_ComputeTurbulentViscosity.cpp
354 
355  } // for k
356  }); // ParFor
357  } // mfi
358 }
KAPPA
constexpr amrex::Real KAPPA
Definition: ERF_Constants.H:20
CONST_GRAV
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoKE_comp
#define RhoKE_comp
Definition: ERF_IndexDefines.H:38
ComputeVerticalDerivativesPBL
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void ComputeVerticalDerivativesPBL(int i, int j, int k, const amrex::Array4< const amrex::Real > &uvel, const amrex::Array4< const amrex::Real > &vvel, const amrex::Array4< const amrex::Real > &cell_data, const int izmin, const int izmax, const amrex::Real &dz_inv, const bool c_ext_dir_on_zlo, const bool c_ext_dir_on_zhi, const bool u_ext_dir_on_zlo, const bool u_ext_dir_on_zhi, const bool v_ext_dir_on_zlo, const bool v_ext_dir_on_zhi, amrex::Real &dthetadz, amrex::Real &dudz, amrex::Real &dvdz, const MoistureComponentIndices &moisture_indices)
Definition: ERF_PBLModels.H:181
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Compute_h_zeta_AtCellCenter
AMREX_FORCE_INLINE AMREX_GPU_DEVICE amrex::Real Compute_h_zeta_AtCellCenter(const int &i, const int &j, const int &k, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:47
Compute_Zrel_AtCellCenter
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real Compute_Zrel_AtCellCenter(const int &i, const int &j, const int &k, const amrex::Array4< const amrex::Real > &z_nd)
Definition: ERF_TerrainMetrics.H:390
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:88
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:87
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:209
EddyDiff::Theta_v
@ Theta_v
Definition: ERF_IndexDefines.H:176
EddyDiff::Turb_lengthscale
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
EddyDiff::Q_v
@ Q_v
Definition: ERF_IndexDefines.H:179
EddyDiff::Mom_v
@ Mom_v
Definition: ERF_IndexDefines.H:175
EddyDiff::KE_v
@ KE_v
Definition: ERF_IndexDefines.H:177
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142

Referenced by ComputeTurbulentViscosity().

Here is the call graph for this function:
Here is the caller graph for this function: