ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Mam4_aer Class Reference

#include <ERF_Mam4_aero.H>

Collaboration diagram for Mam4_aer:

Public Member Functions

void initialize (int ncols, int nlevs, int top_levs, int nsw_bands, int nlw_bands)
 
void read_water_refindex (std::string infilename)
 
void modal_size_parameters (real sigma_logr_aer, const real2d &dgnumwet, const real2d &radsurf, const real2d &logradsurf, const real3d &cheb)
 
void binterp (const real3d &table, int ncol, int km, int im, int jm, const real1d &x, const real1d &y, const real1d &xtab, const real1d &ytab, const int1d &ix, const int1d &jy, const real1d &t, const real1d &u, const real2d &out)
 
void modal_aero_calcsize_diag (int list_idx, const real3d &dgnum_m)
 
void modal_aero_sw (int list_idx, real dt, int nnite, const int1d &idxnite, const bool &is_cmip6_volc, const real2d &pdeldry, const real2d &pmid, const real2d &temperature, const real2d &qt, const real2d &ext_cmip6_sw, const int1d &trop_level, const real3d &tauxar, const real3d &wa, const real3d &ga, const real3d &fa, const real2d &clear_rh)
 
void modal_aero_lw (int list_idx, real dt, const real2d &pdeldry, const real2d &pmid, const real2d &temperature, const real2d &qt, const real3d &tauxar, const real2d &clear_rh)
 

Public Attributes

int ncol
 
int nlev
 
int nspec
 
int top_lev
 
int nswbands
 
int nlwbands
 
std::string modal_optics_file
 
std::string water_refindex_file
 
int ncoef = 5
 
int prefr = 7
 
int prefi = 10
 
int nmodes = 4
 
int n_diag = 1
 
bool clim_modal_aero = true
 
real xrmin
 
real xrmax
 
real1d crefwswr
 
real1d crefwswi
 
real1d crefwlwr
 
real1d crefwlwi
 
real3d dgnumdry_m
 
real3d dgnumwet_m
 
real3d qaerwat_m
 
real3d wetdens_m
 
MamConstituents mam_consti
 

Member Function Documentation

◆ binterp()

void Mam4_aer::binterp ( const real3d &  table,
int  ncol,
int  km,
int  im,
int  jm,
const real1d &  x,
const real1d &  y,
const real1d &  xtab,
const real1d &  ytab,
const int1d &  ix,
const int1d &  jy,
const real1d &  t,
const real1d &  u,
const real2d &  out 
)
inline
168  {
169  // local variables
170  real1d tu("tu", ncol), tuc("tuc", ncol), tcu("tcu", ncol), tcuc("tcuc", ncol);
171 
172  if(ix(1) <= 0) {
173  if(im > 1) {
174  parallel_for (SimpleBounds<1>(ncol), YAKL_LAMBDA (int ic)
175  {
176  int i;
177  for (i = 1; i <= im; ++i) {
178  if(x(ic) < xtab(i)) break;
179  }
180 
181  ix(ic) = std::max(i-1,1);
182  auto ip1 = std::min(ix(ic)+1,im);
183  auto dx = (xtab(ip1)-xtab(ix(ic)));
184 
185  if(abs(dx) > 1.0-20) {
186  t(ic)=(x(ic)-xtab(ix(ic)))/dx;
187  } else {
188  t(ic)=0.;
189  }
190  });
191  } else {
192  parallel_for (SimpleBounds<1>(ncol), YAKL_LAMBDA (int i)
193  {
194  ix(i) = 1;
195  t(i) = 0.;
196  });
197  }
198 
199  if(jm > 1) {
200  parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int ic)
201  {
202  int j;
203  for (j = 1; j <= jm; ++j) {
204  if(y(ic) < ytab(j)) break;
205  }
206 
207  jy(ic) = std::max(j-1,1);
208  auto jp1 = std::min(jy(ic)+1,jm);
209  auto dy = (ytab(jp1)-ytab(jy(ic)));
210 
211  if(std::abs(dy) > 1.e-20) {
212  u(ic) = (y(ic)-ytab(jy(ic)))/dy;
213  if(u(ic) < 0. || u(ic) > 1.)
214  amrex::Print() << "u= " << u(ic) << "; y= " << y(ic) << "; ytab= "
215  << ytab(jy(ic)) << "; dy= " << dy << std::endl;
216  } else {
217  u(ic)=0.;
218  }
219  });
220  } else {
221  parallel_for (SimpleBounds<1>(ncol), YAKL_LAMBDA (int i)
222  {
223  jy(i) = 1;
224  u(i) = 0.;
225  });
226  }
227  }
228 
229  parallel_for (SimpleBounds<2>(ncol, km), YAKL_LAMBDA (int ic, int k)
230  {
231  tu(ic) = t(ic)*u(ic);
232  tuc(ic) = t(ic)-tu(ic);
233  tcuc(ic) = 1.-tuc(ic)-u(ic);
234  tcu(ic) = u(ic)-tu(ic);
235  auto jp1 = std::min(jy(ic)+1,jm);
236  auto ip1 = std::min(ix(ic)+1,im);
237  out(ic,k) = tcuc(ic)*table(k,ix(ic),jy(ic))+tuc(ic)*table(k,ip1 ,jy(ic))
238  +tu(ic)*table(k,ip1 ,jp1 )+tcu(ic)*table(k,ix(ic),jp1 );
239  });
240  } // subroutine binterp
int ncol
Definition: ERF_Mam4_aero.H:25

Referenced by modal_aero_lw(), and modal_aero_sw().

Here is the caller graph for this function:

◆ initialize()

void Mam4_aer::initialize ( int  ncols,
int  nlevs,
int  top_levs,
int  nsw_bands,
int  nlw_bands 
)
inline
58  {
59  ncol = ncols;
60  nlev = nlevs;
61  top_lev = top_levs;
62  nswbands = nsw_bands;
63  nlwbands = nlw_bands;
64  clim_modal_aero = true;
65 
66  const real rmmin = 0.01e-6;
67  const real rmmax = 25.e-6;
68  xrmin = std::log(rmmin);
69  xrmax = std::log(rmmax);
70 
71  // Check that dimension sizes in the coefficient arrays used to
72  // parameterize aerosol radiative properties are consistent between this
73  // module and the mode physprop files.
74  for(auto ilist = 0; ilist < n_diag; ++ilist) {
76  for(auto m = 0; m < nmodes; ++m) {
78  }
79  }
80 
81  auto erf_rad_data_dir = getRadiationDataDir() + "/";
82  water_refindex_file = erf_rad_data_dir + "water_refindex_rrtmg_c080910.nc";
84 
85  // Allocate dry and wet size variables
86  dgnumdry_m = real3d("dgnumdry", ncol, nlev, nmodes);
87  dgnumwet_m = real3d("dgnumwet", ncol, nlev, nmodes);
88  qaerwat_m = real3d("qaerwat" , ncol, nlev, nmodes);
89  wetdens_m = real3d("wetdens" , ncol, nlev, nmodes);
90 
91  crefwswr = real1d("crefwswr", nswbands);
92  crefwswi = real1d("crefwswi", nswbands);
93  crefwlwr = real1d("crefwlwr", nlwbands);
94  crefwlwi = real1d("crefwlwi", nlwbands);
95  }
real3d qaerwat_m
Definition: ERF_Mam4_aero.H:49
real3d dgnumdry_m
Definition: ERF_Mam4_aero.H:47
real1d crefwswr
Definition: ERF_Mam4_aero.H:41
MamConstituents mam_consti
Definition: ERF_Mam4_aero.H:53
real1d crefwlwi
Definition: ERF_Mam4_aero.H:44
void read_water_refindex(std::string infilename)
Definition: ERF_Mam4_aero.H:101
std::string water_refindex_file
Definition: ERF_Mam4_aero.H:28
real3d dgnumwet_m
Definition: ERF_Mam4_aero.H:48
int prefr
Definition: ERF_Mam4_aero.H:33
int nswbands
Definition: ERF_Mam4_aero.H:26
bool clim_modal_aero
Definition: ERF_Mam4_aero.H:37
real xrmin
Definition: ERF_Mam4_aero.H:39
real3d wetdens_m
Definition: ERF_Mam4_aero.H:50
int nlwbands
Definition: ERF_Mam4_aero.H:26
real1d crefwswi
Definition: ERF_Mam4_aero.H:42
int nmodes
Definition: ERF_Mam4_aero.H:35
int nlev
Definition: ERF_Mam4_aero.H:25
int prefi
Definition: ERF_Mam4_aero.H:34
int top_lev
Definition: ERF_Mam4_aero.H:25
real1d crefwlwr
Definition: ERF_Mam4_aero.H:43
int ncoef
Definition: ERF_Mam4_aero.H:32
int n_diag
Definition: ERF_Mam4_aero.H:36
real xrmax
Definition: ERF_Mam4_aero.H:39
void get_nmodes(int list_idx, int &nmodes) const
Definition: ERF_Mam4_constituents.H:350
void get_mode_props(int list_idx, int mode_idx, real &sigmag, real &rhcrystal, real &rhdeliques) const
Definition: ERF_Mam4_constituents.H:808

Referenced by AerRadProps::initialize().

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

◆ modal_aero_calcsize_diag()

void Mam4_aer::modal_aero_calcsize_diag ( int  list_idx,
const real3d &  dgnum_m 
)
inline
244  {
245  real2d dgncur_a("dgncur",ncol,nlev); // (pcols,pver)
246  real2d mode_num("mode_num",ncol,nlev); // mode number mixing ratio
247  real2d specmmr("specmmr",ncol,nlev); // specie mmr
248  real2d dryvol_a("dryvol_a",ncol,nlev); // interstital aerosol dry volume (cm^3/mol_air)
249 
250  for (auto n = 1; n <= nmodes; ++n) {
251  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k)
252  {
253  dgncur_a(i,k) = dgnum_m(i,k,n);
254  });
255 
256  // get mode properties
257  real dgnum, dgnumhi, dgnumlo, sigmag;
258  mam_consti.get_mode_props(list_idx, n-1, dgnum, dgnumhi, dgnumlo, sigmag);
259 
260  // get mode number mixing ratio
261  mam_consti.rad_cnst_get_mode_num(list_idx, n-1, "a", mode_num);
262 
263  yakl::memset(dgncur_a, dgnum);
264  yakl::memset(dryvol_a, 0.0);
265 
266  // compute dry volume mixrats =
267  // sum_over_components{ component_mass mixrat / density }
268  mam_consti.get_mode_nspec(list_idx, n-1, nspec);
269  for (auto l1 = 0; l1 < nspec; ++l1) {
270  real specdens = 1.0e-20;
271  mam_consti.rad_cnst_get_mam_mmr_by_idx(list_idx, n-1, l1, "a", specmmr);
272  mam_consti.get_mam_props(list_idx, n-1, l1, specdens);
273 
274  // need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air
275  real dummwdens = 1.0 / specdens;
276  top_lev = 1;
277 
278  for (auto k=top_lev; k <= nlev; ++k) {
279  for (auto i=1; i<=ncol; ++i) {
280  dryvol_a(i,k) = dryvol_a(i,k) + std::max(0.0, specmmr(i,k))*dummwdens;
281  }
282  }
283  }
284 
285  auto alnsg = std::log(sigmag);
286  auto dumfac = std::exp(4.5*std::pow(alnsg, 2))*PI/6.0;
287  auto voltonumblo = 1./((PI/6.)*(std::pow(dgnumlo, 3))*std::exp(4.5*std::pow(alnsg, 2)));
288  auto voltonumbhi = 1./((PI/6.)*(std::pow(dgnumhi, 3))*std::exp(4.5*std::pow(alnsg, 2)));
289  auto v2nmin = voltonumbhi;
290  auto v2nmax = voltonumblo;
291  auto dgnxx = dgnumhi;
292  auto dgnyy = dgnumlo;
293 
294  top_lev = 1;
295  for (auto k = top_lev; k <= nlev; ++k) {
296  for (auto i = 1; i <= ncol; ++i) {
297  auto drv_a = dryvol_a(i,k);
298  auto num_a0 = mode_num(i,k);
299  auto num_a = std::max(0.0, num_a0);
300 
301  if (drv_a > 0.0) {
302  if (num_a <= drv_a*v2nmin)
303  dgncur_a(i,k) = dgnxx;
304  else if (num_a >= drv_a*v2nmax)
305  dgncur_a(i,k) = dgnyy;
306  else
307  dgncur_a(i,k) = std::pow((drv_a/(dumfac*num_a)), 1./3.);
308  }
309  }
310  }
311  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k)
312  {
313  dgnum_m(i,k,n) = dgncur_a(i,k);
314  });
315  }
316  }
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
int nspec
Definition: ERF_Mam4_aero.H:25
void rad_cnst_get_mode_num(int list_idx, int mode_idx, const std::string &phase, real2d &num) const
Definition: ERF_Mam4_constituents.H:711
void get_mode_nspec(int list_idx, int m_idx, int &nspec) const
Definition: ERF_Mam4_constituents.H:394
void get_mam_props(int list_idx, int mode_idx, int spec_idx, real &density_aer, std::string &spectype, real &hygro_aer, real1d &refindex_real_aer_sw, real1d &refindex_im_aer_sw) const
Definition: ERF_Mam4_constituents.H:919
void rad_cnst_get_mam_mmr_by_idx(int list_idx, int mode_idx, int spec_idx, const std::string &phase, real2d &mmr) const
Definition: ERF_Mam4_constituents.H:639

Referenced by modal_aero_lw(), and modal_aero_sw().

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

◆ modal_aero_lw()

void Mam4_aer::modal_aero_lw ( int  list_idx,
real  dt,
const real2d &  pdeldry,
const real2d &  pmid,
const real2d &  temperature,
const real2d &  qt,
const real3d &  tauxar,
const real2d &  clear_rh 
)
inline
562  {
563  real sigma_logr_aer; // geometric standard deviation of number distribution
564  real alnsg_amode;
565  real1d xrad("xrad", ncol);
566  real3d cheby("cheby", ncoef, ncol, nlev); // chebychef polynomials
567  real2d mass("mass",ncol,nlev); // layer mass
568 
569  //real2d specmmr(:,:) ! species mass mixing ratio
570  real specdens; // species density (kg/m3)
571  real1d specrefrindex("specrefrindex", ncol); // species refractive index
572  real1d specrefiindex("specrefiindex", ncol);
573 
574  real1d vol("vol",ncol); // volume concentration of aerosol specie (m3/kg)
575  real1d dryvol("dryvol",ncol); // volume concentration of aerosol mode (m3/kg)
576  real1d wetvol("wetvol",ncol); // volume concentration of wet mode (m3/kg)
577  real1d watervol("watervol",ncol); // volume concentration of water in each mode (m3/kg)
578  real1d refr("refr",ncol); // real part of refractive index
579  real1d refi("refi",ncol); // imaginary part of refractive index
580  real1d crefinr("crefinr", ncol);
581  real1d crefini("crefini", ncol);
582 
583  real2d refrtablw("refrtablw",prefr,nlwbands); // table of real refractive indices for aerosols
584  real2d refitablw("refitablw",prefi,nlwbands); // table of imag refractive indices for aerosols
585  real4d absplw("absplw",ncoef,prefr,prefi,nlwbands); // specific absorption
586 
587  int1d itab("itab",ncol), jtab("jtab",ncol);
588  real1d ttab("ttab",ncol), utab("utab",ncol);
589  real2d cabs("cabs",ncol,ncoef);
590  real1d pabs("pabs",ncol); // parameterized specific absorption (m2/kg)
591  real1d dopaer("dopaer",ncol); // aerosol optical depth in layer
592 
593  real2d specmmr;
594  real2d dgnumwet; // wet number mode diameter (m)
595  real2d qaerwat; // aerosol water (g/g)
596 
597  constexpr int nerrmax_dopaer=1000;
598  int nerr_dopaer = 0;
599  real volf; // volume fraction of insoluble aerosol
600  amrex::ignore_unused(volf);
601 
602  // initialize output variables
603  yakl::memset(tauxar, 0.);
604 
605  // dry mass in each cell
606  parallel_for (SimpleBounds<2> (ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
607  {
608  mass(icol,ilev) = pdeldry(icol,ilev)*rga;
609  });
610 
611  // Calculate aerosol size distribution parameters and aerosol water uptake
612  if (clim_modal_aero) { // For prescribed aerosol codes
613  //radiation diagnostics are not supported for prescribed aerosols cases
614  if(list_idx != 0)
615  amrex::Print() << "Radiation diagnostic calls are not supported for prescribed aerosols\n";
616  // diagnostic aerosol size calculations
618  }
619 
620  // clear_rh provides alternate estimate non-cloudy relative humidity
622  mam_consti, qt, temperature, pmid, dgnumdry_m, dgnumwet_m,
623  qaerwat_m, wetdens_m, clear_rh);
624 
625  for(auto m = 0; m < nmodes; ++m) {
626  parallel_for (SimpleBounds<2> (ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
627  {
628  dgnumwet(icol, ilev) = dgnumwet_m(icol,ilev,m);
629  qaerwat(icol, ilev) = qaerwat_m(icol,ilev,m);
630  });
631 
632  // get mode properties
633  mam_consti.get_mode_props(list_idx, m, sigma_logr_aer,
634  refrtablw, refitablw, absplw);
635  // get mode info
636  mam_consti.get_mode_nspec(list_idx, m, nspec);
637 
638  // calc size parameter for all columns
639  // this is the same calculation that's done in modal_size_parameters, but there
640  // some intermediate results are saved and the chebyshev polynomials are stored
641  // in a array with different index order. Could be unified.
642  top_lev = 1;
643  for(auto k = top_lev; k <= nlev; ++k) {
644  for(auto i = 1; i <= ncol; ++i) {
645  alnsg_amode = std::log(sigma_logr_aer);
646  // convert from number diameter to surface area
647  xrad(i) = std::log(0.5*dgnumwet(i,k)) + 2.0*alnsg_amode*alnsg_amode;
648  // normalize size parameter
649  xrad(i) = std::max(xrad(i), xrmin);
650  xrad(i) = std::min(xrad(i), xrmax);
651  xrad(i) = (2*xrad(i)-xrmax-xrmin)/(xrmax-xrmin);
652  // chebyshev polynomials
653  cheby(1,i,k) = 1.0;
654  cheby(2,i,k) = xrad(i);
655  for(auto nc = 3; nc <= ncoef; ++nc) {
656  cheby(nc,i,k) = 2.0*xrad(i)*cheby(nc-1,i,k)-cheby(nc-2,i,k);
657  }
658  }
659  }
660 
661  for(auto ilw = 1; ilw <= nlwbands; ++ilw) {
662  for(auto k = top_lev; k <= nlev; ++k) {
663  // form bulk refractive index. Use volume mixing for infrared
664  yakl::memset(crefinr, 0.0);
665  yakl::memset(crefini, 0.0);
666  yakl::memset(dryvol, 0.0);
667  // aerosol species loop
668  for(auto l = 0; l < nspec; ++l) {
669  mam_consti.rad_cnst_get_mam_mmr_by_idx(list_idx, m, l, "a", specmmr);
670  mam_consti.get_mam_props_lw(list_idx, m, l, specdens, specrefrindex, specrefiindex);
671  for(auto i = 1; i <= ncol; ++i) {
672  vol(i) = specmmr(i,k)/specdens;
673  dryvol(i) = dryvol(i) + vol(i);
674  crefinr(i) = crefinr(i) + vol(i)*specrefrindex(ilw);
675  crefini(i) = crefini(i) + vol(i)*specrefiindex(ilw);
676  }
677  }
678  for(auto i = 1; i <= ncol; ++i) {
679  watervol(i) = qaerwat(i,k)/rhoh2o;
680  wetvol(i) = watervol(i) + dryvol(i);
681  if (watervol(i) < 0.0) {
682  if (abs(watervol(i)) > 1.e-1*wetvol(i)) {
683  amrex::Print() << "watervol is too large\n";
684  }
685  watervol(i) = 0.;
686  wetvol(i) = dryvol(i);
687  }
688  crefinr(i) = crefinr(i) + watervol(i)*crefwlwr(ilw);
689  crefini(i) = crefini(i) + watervol(i)*crefwlwi(ilw);
690  if (wetvol(i) > 1.e-40) {
691  crefinr(i) = crefinr(i)/wetvol(i);
692  crefini(i) = crefini(i)/wetvol(i);
693  }
694  refr(i) = crefinr(i);
695  refi(i) = crefini(i);
696  }
697 
698  // interpolate coefficients linear in refractive index
699  // first call calcs itab,jtab,ttab,utab
700  yakl::memset(itab, 0);
701 
702  real3d absplwr("absplwr", ncoef, prefr, prefi);
703  real1d refitablwr("refitablwr", prefi);
704  real1d refrtablwr("refrtablwr", prefr);
705 
706  parallel_for(SimpleBounds<3>(ncoef,prefr,prefi) , YAKL_LAMBDA (int icoef, int irefr, int irefi)
707  {
708  absplwr(icoef,irefr,irefi) = absplw(icoef,irefr,irefi,ilw);
709  refitablwr(irefi) = refitablw(irefi,ilw);
710  refrtablwr(irefr) = refrtablw(irefr,ilw);
711  });
712 
713  binterp(absplwr, ncol, ncoef, prefr, prefi, refr, refi,
714  refrtablwr, refitablwr, itab, jtab, ttab, utab, cabs);
715 
716  // parameterized optical properties
717  for(auto i = 1; i <= ncol; ++i) {
718  pabs(i) = 0.5*cabs(i,1);
719  for(auto nc = 2; nc <= ncoef; ++nc) {
720  pabs(i) = pabs(i) + cheby(nc,i,k)*cabs(i,nc);
721  }
722  pabs(i) = pabs(i)*wetvol(i)*rhoh2o;
723  pabs(i) = std::max(0.,pabs(i));
724  dopaer(i) = pabs(i)*mass(i,k);
725  }
726 
727  for(auto i = 1; i <= ncol; ++i) {
728  if ((dopaer(i) <= -1.e-10) || (dopaer(i) >= 20.)) {
729  if (dopaer(i) <= -1.e-10)
730  amrex::Print() << "ERROR: Negative aerosol optical depth in this layer.\n";
731  else
732  amrex::Print() << "WARNING: Aerosol optical depth is unreasonably high in this layer.\n";
733 
734  for(auto l = 0; l < nspec; ++l) {
735  mam_consti.rad_cnst_get_mam_mmr_by_idx(list_idx, m, l, "a", specmmr);
736  mam_consti.get_mam_props_lw(list_idx, m, l, specdens, specrefrindex, specrefiindex);
737  volf = specmmr(i,k)/specdens;
738  }
739 
740  nerr_dopaer = nerr_dopaer + 1;
741  if (nerr_dopaer >= nerrmax_dopaer || dopaer(i) < -1.e-10) {
742  amrex::Print() << "*** halting after nerr_dopaer = " << nerr_dopaer;
743  exit(EXIT_FAILURE);
744  }
745  }
746  }
747 
748  for(auto i = 1; i <= ncol; ++i) {
749  tauxar(i,k,ilw) = tauxar(i,k,ilw) + dopaer(i);
750  }
751  } // k = top_lev, nlev
752  } // nlwbands
753  } // m = 1, nmodes
754  } // modal_aero_lw
constexpr amrex::Real rga
Definition: ERF_Constants.H:71
constexpr amrex::Real rhoh2o
Definition: ERF_Constants.H:94
void binterp(const real3d &table, int ncol, int km, int im, int jm, const real1d &x, const real1d &y, const real1d &xtab, const real1d &ytab, const int1d &ix, const int1d &jy, const real1d &t, const real1d &u, const real2d &out)
Definition: ERF_Mam4_aero.H:165
void modal_aero_calcsize_diag(int list_idx, const real3d &dgnum_m)
Definition: ERF_Mam4_aero.H:243
void get_mam_props_lw(int list_idx, int mode_idx, int spec_idx, real &density_aer, real1d &refindex_real_aer_lw, real1d &refindex_im_aer_lw) const
Definition: ERF_Mam4_constituents.H:970
static void modal_aero_wateruptake_dr(int list_idx, int ncol, int nlev, int nmodes, int top_lev, MamConstituents &consti, const real2d &h2ommr, const real2d &t, const real2d &pmid, const real3d &dgnumdry_m, const real3d &dgnumwet_m, const real3d &qaerwat_m, const real3d &wetdens_m, const real2d &clear_rh_in)
Definition: ERF_Modal_aero_wateruptake.H:26
@ qt
Definition: ERF_Kessler.H:35

Referenced by AerRadProps::aer_rad_props_lw().

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

◆ modal_aero_sw()

void Mam4_aer::modal_aero_sw ( int  list_idx,
real  dt,
int  nnite,
const int1d &  idxnite,
const bool &  is_cmip6_volc,
const real2d &  pdeldry,
const real2d &  pmid,
const real2d &  temperature,
const real2d &  qt,
const real2d &  ext_cmip6_sw,
const int1d &  trop_level,
const real3d &  tauxar,
const real3d &  wa,
const real3d &  ga,
const real3d &  fa,
const real2d &  clear_rh 
)
inline
327  {
328  // local variables
329  real2d mass("mass",ncol,nlev); // layer mass
330  real2d air_density("air_density",ncol,nlev); // (kg/m3)
331 
332  real2d specmmr("specmmr",ncol,nlev); // species mass mixing ratio
333  real1d specrefindex_real("specrefindex_real",nswbands);
334  real1d specrefindex_im("specrefindex_im",nswbands);
335  std::string spectype; // species type
336  real hygro_aer;
337  real volf;
338  amrex::ignore_unused(volf);
339 
340  real2d dgnumwet("dgnumwet", ncol,nlev); // number mode wet diameter
341  real2d qaerwat("qaerwat", ncol,nlev); // aerosol water (g/g)
342  real2d wetdens("wetdens", ncol, nlev);
343 
344  real sigma_logr_aer; // geometric standard deviation of number distribution
345  real2d radsurf("radsurf", ncol, nlev); // aerosol surface mode radius
346  real2d logradsurf("logradsurf", ncol, nlev); // log(aerosol surface mode radius)
347  real3d cheb("cheb", ncoef, ncol, nlev);
348 
349  real1d refr("refr", ncol); // real part of refractive index
350  real1d refi("refi", ncol); // imaginary part of refractive index
351  real1d crefin_real("crefin_real",ncol);
352  real1d crefin_im("crefin_im",ncol);
353 
354  real2d refitabsw("refitabsw",prefi,nswbands);
355  real2d refrtabsw("refrtabsw",prefr,nswbands);
356 
357  real4d extpsw("extpsw",ncoef,prefr,prefi,nswbands);
358  real4d abspsw("abspsw",ncoef,prefr,prefi,nswbands);
359  real4d asmpsw("asmpsw",ncoef,prefr,prefi,nswbands);
360 
361  real1d vol("vol", ncol); // volume concentration of aerosol specie (m3/kg)
362  real1d dryvol("dryvol", ncol); // volume concentration of aerosol mode (m3/kg)
363  real1d watervol("watervol", ncol); // volume concentration of water in each mode (m3/kg)
364  real1d wetvol("wetvol", ncol); // volume concentration of wet mode (m3/kg)
365 
366  int1d itab("itab", ncol), jtab("jtab", ncol);
367  real1d ttab("ttab",ncol), utab("utab", ncol);
368  real2d cext("cext", ncol, ncoef), cabs("cabs", ncol,ncoef), casm("casm",ncol,ncoef);
369  real1d pext("pext",ncol); // parameterized specific extinction (m2/kg)
370  real1d specpext("specpext",ncol); // specific extinction (m2/kg)
371  real1d dopaer("dopaer",ncol); // aerosol optical depth in layer
372  real1d pabs("pabs",ncol); // parameterized specific absorption (m2/kg)
373  real1d pasm("pasm",ncol); // parameterized asymmetry factor
374  real1d palb("palb",ncol); // parameterized single scattering albedo
375 
376  // initialize output variables
377  yakl::memset(tauxar, 0.);
378  yakl::memset(wa, 0.);
379  yakl::memset(ga, 0.);
380  yakl::memset(fa, 0.);
381 
382  // zero'th layer does not contain aerosol
383  // parallel_for(SimpleBounds<2>(ncol, nswbands), YAKL_LAMBDA (int i, int ibnd) {
384  // wa(i,1,ibnd) = 0.925;
385  // ga(i,1,ibnd) = 0.850;
386  // fa(i,1,ibnd) = 0.7225;
387  // });
388 
389  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k)
390  {
391  mass(i,k) = pdeldry(i,k)*rga;
392  air_density(i,k) = pmid(i,k)/(rair*temperature(i,k));
393  });
394 
395  // Calculate aerosol size distribution parameters and aerosol water uptake
396  if (clim_modal_aero) {
397  // radiation diagnostics are not supported for prescribed aerosols cases
398  if(list_idx != 0) {
399  amrex::Print() << "Radiation diagnostic calls are not supported for prescribed aerosols\n";
400  exit(0);
401  }
402  // diagnostic aerosol size calculations
404  }
405 
406  // clear_rh provides alternate estimate non-cloudy relative humidity
408  mam_consti, qt, temperature, pmid, dgnumdry_m, dgnumwet_m,
409  qaerwat_m, wetdens_m, clear_rh);
410 
411  // loop over all aerosol modes
412  for(auto m = 1; m <= nmodes; ++m) {
413  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev)
414  {
415  dgnumwet(icol,ilev) = dgnumwet_m(icol,ilev,m);
416  qaerwat(icol,ilev) = qaerwat_m(icol,ilev,m);
417  });
418 
419  // get mode properties
420  mam_consti.get_mode_props(list_idx, m-1, sigma_logr_aer, refrtabsw,
421  refitabsw, extpsw, abspsw, asmpsw);
422 
423  // get mode info
424  mam_consti.get_mode_nspec(list_idx, m-1, nspec);
425 
426  // calc size parameter for all columns
427  modal_size_parameters(sigma_logr_aer, dgnumwet, radsurf, logradsurf, cheb);
428 
429  for(auto isw = 1; isw <= nswbands; ++isw) {
430  for(auto k = top_lev; k <= nlev; ++k) {
431  // form bulk refractive index
432  yakl::memset(crefin_real, 0.);
433  yakl::memset(crefin_im , 0.);
434  yakl::memset(dryvol , 0.);
435 
436  // aerosol species loop
437  for(auto l = 0; l < nspec; ++l) {
438  real specdens = 1.0e-20;
439  mam_consti.rad_cnst_get_mam_mmr_by_idx(list_idx, m-1, l, "a", specmmr);
440  mam_consti.get_mam_props(list_idx, m-1, l, specdens, spectype, hygro_aer,
441  specrefindex_real, specrefindex_im);
442  for(auto i = 1; i <= ncol; ++i) {
443  vol(i) = specmmr(i,k)/specdens;
444  dryvol(i) = dryvol(i) + vol(i);
445  crefin_real(i) = crefin_real(i) + vol(i)*specrefindex_real(isw);
446  crefin_im(i) = crefin_im(i) + vol(i)*specrefindex_im(isw);
447  }
448  }
449 
450  for(auto i = 1; i <= ncol; ++i) {
451  watervol(i) = qaerwat(i,k)/rhoh2o;
452  wetvol(i) = watervol(i) + dryvol(i);
453  if (watervol(i) < 0.) {
454  if (std::abs(watervol(i)) > 1.e-1*wetvol(i)) {
455  amrex::Print() << "watervol,wetvol=" << watervol(i) << "; " << wetvol(i) << std::endl;
456  }
457  watervol(i) = 0.;
458  wetvol(i) = dryvol(i);
459  }
460  // volume mixing
461  crefin_real(i) = crefin_real(i)+watervol(i)*crefwswr(isw);
462  crefin_real(i) = crefin_real(i)/std::max(wetvol(i),1.e-60);
463  crefin_im(i) = crefin_im(i)+watervol(i)*crefwswi(isw);
464  crefin_im(i) = crefin_im(i)/std::max(wetvol(i),1.e-60);
465  refr(i) = crefin_real(i);
466  refi(i) = std::abs(crefin_im(i));
467  }
468 
469  // interpolate coefficients linear in refractive index
470  // first call calcs itab,jtab,ttab,utab
471  real3d extpswr("extpswr", ncoef, prefr, prefi);
472  real3d abspswr("abspswr", ncoef, prefr, prefi);
473  real3d asmpswr("asmpswr", ncoef, prefr, prefi);
474  real1d refitabswr("refitabswr", prefi);
475  real1d refrtabswr("refrtabswr", prefr);
476 
477  parallel_for(SimpleBounds<3>(ncoef,prefr,prefi), YAKL_LAMBDA (int icoef, int irefr, int irefi)
478  {
479  extpswr(icoef,irefr,irefi) = extpsw(icoef,irefr,irefi,isw);
480  abspswr(icoef,irefr,irefi) = abspsw(icoef,irefr,irefi,isw);
481  asmpswr(icoef,irefr,irefi) = asmpsw(icoef,irefr,irefi,isw);
482  refitabswr(irefi) = refitabsw(irefi,isw);
483  refrtabswr(irefr) = refrtabsw(irefr,isw);
484  });
485 
486  yakl::memset(itab, 0);
487  binterp(extpswr, ncol, ncoef, prefr, prefi, refr, refi,
488  refrtabswr, refitabswr, itab, jtab, ttab, utab, cext);
489  binterp(abspswr, ncol, ncoef, prefr, prefi, refr, refi,
490  refrtabswr, refitabswr, itab, jtab, ttab, utab, cabs);
491  binterp(asmpswr, ncol, ncoef, prefr, prefi, refr, refi,
492  refrtabswr, refitabswr, itab, jtab, ttab, utab, casm);
493 
494  // parameterized optical properties
495  for(auto i=1; i <= ncol; ++i) {
496  if (logradsurf(i,k) <= xrmax) {
497  pext(i) = 0.5*cext(i,1);
498  for(auto nc = 2; nc <= ncoef; ++nc) {
499  pext(i) = pext(i) + cheb(nc,i,k)*cext(i,nc);
500  }
501  pext(i) = std::exp(pext(i));
502  } else {
503  pext(i) = 1.5/(radsurf(i,k)*rhoh2o); // geometric optics
504  }
505 
506  // convert from m2/kg water to m2/kg aerosol
507  specpext(i) = pext(i);
508  pext(i) = pext(i)*wetvol(i)*rhoh2o;
509  pabs(i) = 0.5*cabs(i,1);
510  pasm(i) = 0.5*casm(i,1);
511  for(auto nc = 2; nc <= ncoef; ++nc) {
512  pabs(i) = pabs(i) + cheb(nc,i,k)*cabs(i,nc);
513  pasm(i) = pasm(i) + cheb(nc,i,k)*casm(i,nc);
514  }
515  pabs(i) = pabs(i)*wetvol(i)*rhoh2o;
516  pabs(i) = std::max(0.,pabs(i));
517  pabs(i) = std::min(pext(i),pabs(i));
518  palb(i) = 1.-pabs(i)/std::max(pext(i),1.e-40);
519  dopaer(i) = pext(i)*mass(i,k);
520  }
521 
522  for (auto i = 1; i <= ncol; ++i) {
523  if ((dopaer(i) <= -1.e-10) || (dopaer(i) >= 30.)) {
524  if (dopaer(i) <= -1.e-10)
525  amrex::Print() << "ERROR: Negative aerosol optical depth in this layer.\n";
526  else {
527  // reset to the bound value
528  dopaer(i) = 25.0;
529  amrex::Print() << "WARNING: Aerosol optical depth is unreasonably high in this layer.\n";
530  }
531 
532  for(auto l = 1; l < nspec; ++l) {
533  real specdens = 1.0e-20;
534  mam_consti.rad_cnst_get_mam_mmr_by_idx(list_idx, m-1, l, "a", specmmr);
535  mam_consti.get_mam_props_sw(list_idx, m-1, l, specdens, specrefindex_real, specrefindex_im);
536  volf = specmmr(i,k)/specdens;
537  }
538  if (dopaer(i) < -1.e-10) {
539  amrex::Print() << "*** halting with error!\n";
540  exit(0);
541  }
542  }
543  }
544 
545  for(auto i=1; i <= ncol; ++i) {
546  tauxar(i,k,isw) = tauxar(i,k,isw) + dopaer(i);
547  wa(i,k,isw) = wa(i,k,isw) + dopaer(i)*palb(i);
548  ga(i,k,isw) = ga(i,k,isw) + dopaer(i)*palb(i)*pasm(i);
549  fa(i,k,isw) = fa(i,k,isw) + dopaer(i)*palb(i)*pasm(i)*pasm(i);
550  }
551  } //end do ! nlev
552  } //end do ! swbands
553  } // nmodes
554  } // modal_aero_sw
constexpr amrex::Real rair
Definition: ERF_Constants.H:69
void modal_size_parameters(real sigma_logr_aer, const real2d &dgnumwet, const real2d &radsurf, const real2d &logradsurf, const real3d &cheb)
Definition: ERF_Mam4_aero.H:130
void get_mam_props_sw(int list_idx, int mode_idx, int spec_idx, real &density_aer, real1d &refindex_real_aer_sw, real1d &refindex_im_aer_sw) const
Definition: ERF_Mam4_constituents.H:939

Referenced by AerRadProps::aer_rad_props_sw().

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

◆ modal_size_parameters()

void Mam4_aer::modal_size_parameters ( real  sigma_logr_aer,
const real2d &  dgnumwet,
const real2d &  radsurf,
const real2d &  logradsurf,
const real3d &  cheb 
)
inline
133  {
134  real1d xrad("xrad", ncol);
135 
136  auto alnsg_amode = std::log(sigma_logr_aer);
137  auto explnsigma = std::exp(2.0*alnsg_amode*alnsg_amode);
138  top_lev = 1;
139 
140  for(auto k = top_lev; k <= nlev; ++k) {
141  for(auto i = 1; i <= ncol; ++i) {
142  // convert from number mode diameter to surface area
143  radsurf(i,k) = 0.5*dgnumwet(i,k)*explnsigma;
144  logradsurf(i,k) = std::log(radsurf(i,k));
145 
146  // normalize size parameter
147  xrad(i) = std::max(logradsurf(i,k),xrmin);
148  xrad(i) = std::min(xrad(i),xrmax);
149  xrad(i) = (2.*xrad(i)-xrmax-xrmin)/(xrmax-xrmin);
150 
151  // chebyshev polynomials
152  cheb(1,i,k) = 1.;
153  cheb(2,i,k) = xrad(i);
154  for(auto nc = 3; nc <= ncoef; ++nc) {
155  cheb(nc,i,k) = 2.*xrad(i)*cheb(nc-1,i,k)-cheb(nc-2,i,k);
156  }
157  }
158  }
159  }

Referenced by modal_aero_sw().

Here is the caller graph for this function:

◆ read_water_refindex()

void Mam4_aer::read_water_refindex ( std::string  infilename)
inline
102  {
103  yakl::SimpleNetCDF water_ref_file;
104  water_ref_file.open(infilename, yakl::NETCDF_MODE_READ);
105 
106  // read the dimensions
107  int nlw_bands = water_ref_file.getDimSize( "lw_band" );
108  int nsw_bands = water_ref_file.getDimSize( "sw_band" );
109 
110  if (nswbands != nsw_bands || nlwbands != nlw_bands)
111  amrex::Print() << "ERROR - file and bandwidth values do not match: "
112  << "\n nswbands: " << nswbands << "; nsw_bands: " << nsw_bands
113  << "\n nlwbands: " << nlwbands << "; nlw_bands: " << nlw_bands << "\n";
114 
115  // Local variables
116  real1d refrwsw("refrwsw", nsw_bands), refiwsw("refiwsw", nsw_bands); // real, imaginary ref index for water visible
117  real1d refrwlw("refrwlw", nlw_bands), refiwlw("refiwlw", nlw_bands); // real, imaginary ref index for water infrared
118 
119  // read variables
120  water_ref_file.read( refrwsw, "refindex_real_water_sw");
121  water_ref_file.read( refiwsw, "refindex_im_water_sw");
122  water_ref_file.read( refrwlw, "refindex_real_water_lw");
123  water_ref_file.read( refiwlw, "refindex_im_water_lw");
124  }

Referenced by initialize().

Here is the caller graph for this function:

Member Data Documentation

◆ clim_modal_aero

bool Mam4_aer::clim_modal_aero = true

◆ crefwlwi

real1d Mam4_aer::crefwlwi

Referenced by initialize(), and modal_aero_lw().

◆ crefwlwr

real1d Mam4_aer::crefwlwr

Referenced by initialize(), and modal_aero_lw().

◆ crefwswi

real1d Mam4_aer::crefwswi

Referenced by initialize(), and modal_aero_sw().

◆ crefwswr

real1d Mam4_aer::crefwswr

Referenced by initialize(), and modal_aero_sw().

◆ dgnumdry_m

real3d Mam4_aer::dgnumdry_m

◆ dgnumwet_m

real3d Mam4_aer::dgnumwet_m

◆ mam_consti

◆ modal_optics_file

std::string Mam4_aer::modal_optics_file

◆ n_diag

int Mam4_aer::n_diag = 1

Referenced by initialize().

◆ ncoef

int Mam4_aer::ncoef = 5

◆ ncol

◆ nlev

◆ nlwbands

int Mam4_aer::nlwbands

◆ nmodes

int Mam4_aer::nmodes = 4

◆ nspec

int Mam4_aer::nspec

◆ nswbands

int Mam4_aer::nswbands

◆ prefi

int Mam4_aer::prefi = 10

◆ prefr

int Mam4_aer::prefr = 7

◆ qaerwat_m

real3d Mam4_aer::qaerwat_m

◆ top_lev

◆ water_refindex_file

std::string Mam4_aer::water_refindex_file

Referenced by initialize().

◆ wetdens_m

real3d Mam4_aer::wetdens_m

◆ xrmax

real Mam4_aer::xrmax

◆ xrmin

real Mam4_aer::xrmin

The documentation for this class was generated from the following file: