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 (double 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 std::unique_ptr< MultiFab > &z_phys_cc, const MoistureComponentIndices &moisture_indices)
 

Function Documentation

◆ ComputeDiffusivityMYJ()

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

Referenced by ComputeTurbulentViscosity().

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