ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ERF_PBLModels.H
Go to the documentation of this file.
1 #ifndef ERF_PBLMODELS_H_
2 #define ERF_PBLMODELS_H_
3 
4 #include "ERF_DataStruct.H"
5 
6 /**
7  * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
8  * Mellor-Yamada-Nakanishi-Niino Level 2.5 PBL scheme
9  *
10  * @param[in] xvel Velocity in x-dir
11  * @param[in] yvel Velocity in y-dir
12  * @param[in] cons_in Cell center conserved quantities
13  * @param[out] eddyViscosity Holds turbulent viscosity
14  * @param[in] geom Problem geometry
15  * @param[in] turbChoice Container with turbulence parameters
16  * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
17  * @param[in] use_moisture If we have microphysics enabled
18  * @param[in] level Current level
19  * @param[in] bc_ptr Pointer to array with boundary condition info
20  * @param[in] vert_only Only compute vertical eddy diffusivities
21  * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
22  */
23 void
24 ComputeDiffusivityMYNN25 (const amrex::MultiFab& xvel,
25  const amrex::MultiFab& yvel,
26  const amrex::MultiFab& cons_in,
27  amrex::MultiFab& eddyViscosity,
28  const amrex::Geometry& geom,
29  const TurbChoice& turbChoice,
30  std::unique_ptr<SurfaceLayer>& SurfLayer,
31  bool use_terrain_fitted_coords,
32  bool use_moisture,
33  int level,
34  const amrex::BCRec* bc_ptr,
35  bool /*vert_only*/,
36  const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
37  const MoistureComponentIndices& moisture_indices);
38 
39 /**
40  * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
41  * Mellor-Yamada-Nakanishi-Niino Level 2.5 PBL scheme
42  *
43  * @param[in] xvel Velocity in x-dir
44  * @param[in] yvel Velocity in y-dir
45  * @param[in] cons_in Cell center conserved quantities
46  * @param[out] eddyViscosity Holds turbulent viscosity
47  * @param[in] geom Problem geometry
48  * @param[in] turbChoice Container with turbulence parameters
49  * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
50  * @param[in] use_moisture If we have microphysics enabled
51  * @param[in] level Current level
52  * @param[in] bc_ptr Pointer to array with boundary condition info
53  * @param[in] vert_only Only compute vertical eddy diffusivities
54  * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
55  */
56 void
57 ComputeDiffusivityMYNNEDMF (const amrex::MultiFab& xvel,
58  const amrex::MultiFab& yvel,
59  const amrex::MultiFab& cons_in,
60  amrex::MultiFab& eddyViscosity,
61  const amrex::Geometry& geom,
62  const TurbChoice& turbChoice,
63  std::unique_ptr<SurfaceLayer>& SurfLayer,
64  bool use_terrain_fitted_coords,
65  bool use_moisture,
66  int level,
67  const amrex::BCRec* bc_ptr,
68  bool /*vert_only*/,
69  const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
70  const MoistureComponentIndices& moisture_indices);
71 
72 /**
73  * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
74  * Yonsei University PBL scheme
75  *
76  * @param[in] xvel Velocity in x-dir
77  * @param[in] yvel Velocity in y-dir
78  * @param[in] cons_in Cell center conserved quantities
79  * @param[out] eddyViscosity Holds turbulent viscosity
80  * @param[in] geom Problem geometry
81  * @param[in] turbChoice Container with turbulence parameters
82  * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
83  * @param[in] use_moisture If we have microphysics enabled
84  * @param[in] level Current level
85  * @param[in] bc_ptr Pointer to array with boundary condition info
86  * @param[in] vert_only Only compute vertical eddy diffusivities
87  * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
88  */
89 void
90 ComputeDiffusivityYSU (const amrex::MultiFab& xvel,
91  const amrex::MultiFab& yvel,
92  const amrex::MultiFab& cons_in,
93  amrex::MultiFab& eddyViscosity,
94  const amrex::Geometry& geom,
95  const TurbChoice& turbChoice,
96  std::unique_ptr<SurfaceLayer>& SurfLayer,
97  bool use_terrain_fitted_coords,
98  bool use_moisture,
99  int level,
100  const amrex::BCRec* bc_ptr,
101  bool /*vert_only*/,
102  const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
103  const MoistureComponentIndices& moisture_indices);
104 /**
105  * Compute eddy diffusivities of momentum (eddy viscosity) and heat using the
106  * older version of Yonsei University PBL scheme
107  *
108  * @param[in] xvel Velocity in x-dir
109  * @param[in] yvel Velocity in y-dir
110  * @param[in] cons_in Cell center conserved quantities
111  * @param[out] eddyViscosity Holds turbulent viscosity
112  * @param[in] geom Problem geometry
113  * @param[in] turbChoice Container with turbulence parameters
114  * @param[in] SurfLayer Pointer to Monin-Obukhov class if instantiated
115  * @param[in] use_moisture If we have microphysics enabled
116  * @param[in] level Current level
117  * @param[in] bc_ptr Pointer to array with boundary condition info
118  * @param[in] vert_only Only compute vertical eddy diffusivities
119  * @param[in] z_phys_nd Physical location of grid nodes, if terrain (or grid stretching) is enabled
120  */
121 void
122 ComputeDiffusivityMRF (const amrex::MultiFab& xvel,
123  const amrex::MultiFab& yvel,
124  const amrex::MultiFab& cons_in,
125  amrex::MultiFab& eddyViscosity,
126  const amrex::Geometry& geom,
127  const TurbChoice& turbChoice,
128  std::unique_ptr<SurfaceLayer>& SurfLayer,
129  bool use_terrain_fitted_coords,
130  bool use_moisture,
131  int level,
132  const amrex::BCRec* bc_ptr,
133  bool /*vert_only*/,
134  const std::unique_ptr<amrex::MultiFab>& z_phys_nd,
135  const MoistureComponentIndices& moisture_indices);
136 
137 /**
138  * Function for computing vertical derivatives for use in PBL model
139  *
140  * @param[in] u velocity in x-dir
141  * @param[in] v velocity in y-dir
142  * @param[in] cell_data conserved cell center vars
143  */
144 AMREX_GPU_DEVICE
145 AMREX_FORCE_INLINE
146 void
147 ComputeVerticalDerivativesPBL (int i, int j, int k,
148  const amrex::Array4<const amrex::Real>& uvel,
149  const amrex::Array4<const amrex::Real>& vvel,
150  const amrex::Array4<const amrex::Real>& cell_data,
151  const int izmin,
152  const int izmax,
153  const amrex::Real& dz_inv,
154  const bool c_ext_dir_on_zlo,
155  const bool c_ext_dir_on_zhi,
156  const bool u_ext_dir_on_zlo,
157  const bool u_ext_dir_on_zhi,
158  const bool v_ext_dir_on_zlo,
159  const bool v_ext_dir_on_zhi,
160  amrex::Real& dthetadz,
161  amrex::Real& dudz,
162  amrex::Real& dvdz,
163  const MoistureComponentIndices& moisture_indices)
164 {
165  if ( k==izmax && c_ext_dir_on_zhi ) {
166  dthetadz = (1.0/3.0)*( -GetThetav(i,j,k-1,cell_data,moisture_indices)
167  - 3.0 * GetThetav(i,j,k ,cell_data,moisture_indices)
168  + 4.0 * GetThetav(i,j,k+1,cell_data,moisture_indices) )*dz_inv;
169  } else if ( k==izmin && c_ext_dir_on_zlo ) {
170  dthetadz = (1.0/3.0)*( GetThetav(i,j,k+1,cell_data,moisture_indices)
171  + 3.0 * GetThetav(i,j,k ,cell_data,moisture_indices)
172  - 4.0 * GetThetav(i,j,k-1,cell_data,moisture_indices) )*dz_inv;
173  } else {
174  dthetadz = 0.5*( GetThetav(i,j,k+1,cell_data,moisture_indices)
175  - GetThetav(i,j,k-1,cell_data,moisture_indices) )*dz_inv;
176  }
177 
178  if ( k==izmax && u_ext_dir_on_zhi ) {
179  dudz = (1.0/6.0)*( (-uvel(i ,j,k-1) - 3.0 * uvel(i ,j,k ) + 4.0 * uvel(i ,j,k+1))
180  + (-uvel(i+1,j,k-1) - 3.0 * uvel(i+1,j,k ) + 4.0 * uvel(i+1,j,k+1)) )*dz_inv;
181  } else if ( k==izmin && u_ext_dir_on_zlo ) {
182  dudz = (1.0/6.0)*( (uvel(i ,j,k+1) + 3.0 * uvel(i ,j,k ) - 4.0 * uvel(i ,j,k-1))
183  + (uvel(i+1,j,k+1) + 3.0 * uvel(i+1,j,k ) - 4.0 * uvel(i+1,j,k-1)) )*dz_inv;
184  } else {
185  dudz = 0.25*( uvel(i,j,k+1) - uvel(i,j,k-1) + uvel(i+1,j,k+1) - uvel(i+1,j,k-1) )*dz_inv;
186  }
187 
188  if ( k==izmax && v_ext_dir_on_zhi ) {
189  dvdz = (1.0/6.0)*( (-vvel(i,j ,k-1) - 3.0 * vvel(i,j ,k ) + 4.0 * vvel(i,j ,k+1))
190  + (-vvel(i,j+1,k-1) - 3.0 * vvel(i,j+1,k ) + 4.0 * vvel(i,j+1,k+1)) )*dz_inv;
191  } else if ( k==izmin && v_ext_dir_on_zlo ) {
192  dvdz = (1.0/6.0)*( (vvel(i,j ,k+1) + 3.0 * vvel(i,j ,k ) - 4.0 * vvel(i,j ,k-1))
193  + (vvel(i,j+1,k+1) + 3.0 * vvel(i,j+1,k ) - 4.0 * vvel(i,j+1,k-1)) )*dz_inv;
194  } else {
195  dvdz = 0.25*( vvel(i,j,k+1) - vvel(i,j,k-1) + vvel(i,j+1,k+1) - vvel(i,j+1,k-1) )*dz_inv;
196  }
197 }
198 
199 /**
200  * Function for computing the QKE source terms (NN09, Eqn. 5).
201  *
202  * @param[in] u velocity in x-dir
203  * @param[in] v velocity in y-dir
204  * @param[in] cell_data conserved cell center vars
205  * @param[in] cell_prim primitive cell center vars
206  * @param[in] K_turb turbulent viscosity
207  * @param[in] cellSizeInv inverse cell size array
208  * @param[in] domain box of the whole domain
209  * @param[in] pbl_mynn_B1_l a parameter
210  * @param[in] theta_mean average theta
211  */
212 AMREX_GPU_DEVICE
213 AMREX_FORCE_INLINE
214 amrex::Real
215 ComputeQKESourceTerms (int i, int j, int k,
216  const amrex::Array4<const amrex::Real>& uvel,
217  const amrex::Array4<const amrex::Real>& vvel,
218  const amrex::Array4<const amrex::Real>& cell_data,
219  const amrex::Array4<const amrex::Real>& cell_prim,
220  const amrex::Array4<const amrex::Real>& K_turb,
221  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSizeInv,
222  const amrex::Box& domain,
223  amrex::Real pbl_mynn_B1_l,
224  const amrex::Real theta_mean,
225  const MoistureComponentIndices& moisture_indices,
226  bool c_ext_dir_on_zlo,
227  bool c_ext_dir_on_zhi,
228  bool u_ext_dir_on_zlo,
229  bool u_ext_dir_on_zhi,
230  bool v_ext_dir_on_zlo,
231  bool v_ext_dir_on_zhi,
232  const amrex::Real met_h_zeta=1.0)
233 {
234  // Compute some relevant derivatives
235  amrex::Real dthetadz, dudz, dvdz;
236  amrex::Real source_term = 0.0;
237 
238  amrex::Real dz_inv = cellSizeInv[2];
239  int izmin = domain.smallEnd(2);
240  int izmax = domain.bigEnd(2);
241 
243  uvel, vvel, cell_data, izmin, izmax, dz_inv/met_h_zeta,
244  c_ext_dir_on_zlo, c_ext_dir_on_zhi,
245  u_ext_dir_on_zlo, u_ext_dir_on_zhi,
246  v_ext_dir_on_zlo, v_ext_dir_on_zhi,
247  dthetadz, dudz, dvdz,
248  moisture_indices);
249 
250  // Notes:
251  // - We transport TKE = 0.5*QKE rather than QKE, so the RHS terms do not
252  // have a factor of two.
253  // - Transport terms due to turbulence and pressure are included when
254  // DiffusionSrcForState_* is called from ERF_slow_rhs_post.
255  // - Eddy diffusivities are updated at the beginning of each time step only.
256 
257  // Second-order turbulent fluxes, e.g.:
258  // -<uw> = L q SM dU/dz (NN09, Eqn. 18)
259  // = Kmv/rho dU/dz
260 
261  // Shear Production
262  source_term += K_turb(i,j,k,EddyDiff::Mom_v) * (dudz*dudz + dvdz*dvdz);
263 
264  // Buoyancy Production
265  source_term -= (CONST_GRAV/theta_mean)*K_turb(i,j,k,EddyDiff::Theta_v)*dthetadz;
266 
267  // Dissipation (NN09, Eqn. 12)
268  amrex::Real qke = 2.0 * cell_prim(i,j,k,PrimKE_comp);
269  if (std::abs(qke) > 0.0) {
270  source_term -= cell_data(i,j,k,Rho_comp) * std::pow(qke,1.5) /
271  (pbl_mynn_B1_l * K_turb(i,j,k,EddyDiff::Turb_lengthscale));
272  }
273 
274  return source_term;
275 }
276 #endif
constexpr amrex::Real CONST_GRAV
Definition: ERF_Constants.H:21
TurbChoice turbChoice
Definition: ERF_DiffSetup.H:2
#define Rho_comp
Definition: ERF_IndexDefines.H:36
#define PrimKE_comp
Definition: ERF_IndexDefines.H:51
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real GetThetav(int i, int j, int k, const amrex::Array4< amrex::Real const > &cell_data, const MoistureComponentIndices &moisture_indices)
Definition: ERF_MoistUtils.H:72
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real ComputeQKESourceTerms(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 amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &K_turb, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &cellSizeInv, const amrex::Box &domain, amrex::Real pbl_mynn_B1_l, const amrex::Real theta_mean, const MoistureComponentIndices &moisture_indices, bool c_ext_dir_on_zlo, bool c_ext_dir_on_zhi, bool u_ext_dir_on_zlo, bool u_ext_dir_on_zhi, bool v_ext_dir_on_zlo, bool v_ext_dir_on_zhi, const amrex::Real met_h_zeta=1.0)
Definition: ERF_PBLModels.H:215
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:147
void ComputeDiffusivityMYNN25(const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
void ComputeDiffusivityYSU(const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
void ComputeDiffusivityMYNNEDMF(const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
void ComputeDiffusivityMRF(const amrex::MultiFab &xvel, const amrex::MultiFab &yvel, const amrex::MultiFab &cons_in, amrex::MultiFab &eddyViscosity, const amrex::Geometry &geom, const TurbChoice &turbChoice, std::unique_ptr< SurfaceLayer > &SurfLayer, bool use_terrain_fitted_coords, bool use_moisture, int level, const amrex::BCRec *bc_ptr, bool, const std::unique_ptr< amrex::MultiFab > &z_phys_nd, const MoistureComponentIndices &moisture_indices)
@ Theta_v
Definition: ERF_IndexDefines.H:176
@ Turb_lengthscale
Definition: ERF_IndexDefines.H:180
@ Mom_v
Definition: ERF_IndexDefines.H:175
@ xvel
Definition: ERF_IndexDefines.H:141
@ yvel
Definition: ERF_IndexDefines.H:142
Definition: ERF_DataStruct.H:99
Definition: ERF_TurbStruct.H:41