ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_m2005_effradius.H
Go to the documentation of this file.
1 /*-------------------------------------------------------------------------------------------
2  This subroutine is used to calculate droplet and ice crystal effective radius, which will be
3  used in the CAM radiation code. The method to calculate effective radius is taken out of the
4  Morrision two moment scheme from M2005MICRO_GRAUPEL. It is also very similar to the subroutine
5  effradius in the module of cldwat2m in the CAM source codes.
6  Adopted by Minghuai Wang (Minghuai.Wang@pnl.gov).
7  Calculate effective radius for radiation code
8  If no cloud water, default value is:
9  10 micron for droplets,
10  25 micron for cloud ice.
11  Be careful of the unit of effective radius : [micro meter]
12 
13  NOTE: this code is modified from E3SM
14 -------------------------------------------------------------------------------------------*/
15 #include "ERF_Constants.H"
16 #include "ERF_Rad_constants.H"
17 #include <cmath>
18 
19 using namespace amrex;
20 using yakl::intrinsics::size;
21 using yakl::fortran::parallel_for;
22 using yakl::fortran::SimpleBounds;
23 
24 YAKL_INLINE
25 void m2005_effradius (const real2d& ql, const real2d& nl, const real2d& qi,
26  const real2d& ni, const real2d& qs, const real2d& ns,
27  const real2d& cld, const real2d& pres, const real2d& tk,
28  const real2d& effl, const real2d& effi,
29  const real2d& deffi, const real2d& lamcrad,
30  const real2d& pgamrad, const real2d& des)
31 {
32  // Main computation
33  const real pi = 3.1415926535897932384626434;
34  const real qsmall = 1.0e-14; // in the SAM source code (module_mp_graupel)
35  const real rhow = 997.; // in module_mp_graupel, SAM
36  const real rhoi = 500.; // in both CAM and SAM
37  const real dcs = 125.e-6; // in module_mp_graupel, SAM
38  const real ci = rhoi*pi/6.;
39  const real di = 3.;
40 
41  // for snow water
42  const real rhos = 100.;
43  const real cs = rhos*pi/6.;
44  const real ds = 3.;
45  const real mincld = 0.0001;
46 
47  int ncol = pres.extent(0);
48  int nlev = pres.extent(1);
49  // Effective diameters of snow crystals
50  parallel_for (SimpleBounds<2> (ncol, nlev), YAKL_LAMBDA (int i, int j)
51  {
52  auto rho = pres(i,j)/(287.15*tk(i,j)); // air density [kg/m3]
53  auto cldm = std::max(cld(i,j), mincld);
54  auto qlic = std::min(5.e-3, std::max(0., ql(i,j)/cldm));
55  auto qiic = std::min(5.e-3, std::max(0., qi(i,j)/cldm));
56  auto nlic = std::max(nl(i,j), 0.)/cldm;
57  auto niic = std::max(ni(i,j), 0.)/cldm;
58 
59  real res;
60  if(qs(i,j) > 1.0e-7) {
61  auto lammaxs=1./10.e-6;
62  auto lammins=1./2000.e-6;
63  auto lams = std::pow((std::tgamma(1.+ds)*cs*ns(i,j)/qs(i,j)), (1./ds));
64  lams = std::min(lammaxs, std::max(lams,lammins));
65  res = 1.5/lams*1.0e6;
66  }
67  else {
68  res = 500.;
69  }
70 
71  // from Hugh Morrision: rhos/917 accounts for assumptions about
72  // ice density in the Mitchell optics.
73  des(i,j) = res*rhos/917.*2.;
74 
75  // Effective radius of cloud ice droplet
76  if( qiic >= qsmall ) {
77  niic = std::min(niic, qiic*1.e20);
78  auto lammax = 1./1.e-6; // in module_mp_graupel, SAM
79  auto lammin = 1./(2.*dcs+100.e-6); // in module_mp_graupel, SAM
80  auto lami = std::pow((std::tgamma(1.+di)*ci*niic/qiic), (1./di));
81  lami = std::min(lammax, std::max(lami, lammin));
82  effi(i,j) = 1.5/lami*1.e6;
83  }
84  else {
85  effi(i,j) = 25.;
86  }
87 
88  // hm ice effective radius for david mitchell's optics
89  // ac morrison indicates that this is effective diameter
90  // ac morrison indicates 917 (for the density of pure ice..)
91  deffi(i,j) = effi(i,j)*rhoi/917.*2.;
92 
93  // Effective radius of cloud liquid droplet
94  if( qlic >= qsmall ) {
95  // Matin et al., 1994 (JAS) formula for pgam (the same is used in both CAM and SAM).
96  // See also Morrison and Grabowski (2007, JAS, Eq. (2))
97  nlic = std::min(nlic, qlic*1.e20);
98 
99  // set the minimum droplet number as 20/cm3.
100  // nlic = max(nlic,20.e6_r8/rho) ! sghan minimum in #/cm3
101  //auto tempnc = nlic/rho/1.0e6; // #/kg --> #/cm3
102 
103  // Should be the in-cloud dropelt number calculated as nlic*rho/1.0e6_r8 ????!!!! +++mhwang
104  auto pgam = 0.0005714*(nlic*rho/1.e6) + 0.2714;
105  pgam = 1./std::pow(pgam, 2.0)-1.;
106  pgam = std::min(10., std::max(pgam,2.));
107  auto laml = std::pow((pi/6.*rhow*nlic*std::tgamma(pgam+4.)/(qlic*std::tgamma(pgam+1.))), (1./3.));
108  auto lammin = (pgam+1.)/50.e-6; // in cldwat2m, CAM
109  auto lammax = (pgam+1.)/2.e-6; // if lammax is too large, this will lead to crash in
110  // src/physics/rrtmg/cloud_rad_props.F90 because
111  // klambda-1 can be zero in gam_liquid_lw and gam_liquid_sw
112  // and g_lambda(kmu,klambda-1) will not be defined.
113  laml = std::min(std::max(laml, lammin),lammax);
114  effl(i,j) = std::tgamma(pgam+4.)/std::tgamma(pgam+3.)/laml/2.*1.e6;
115  lamcrad(i,j) = laml;
116  pgamrad(i,j) = pgam;
117  }
118  else {
119  // we chose 10. over 25, since 10 is a more reasonable value for liquid droplet. +++mhwang
120  effl(i,j) = 10.; // in cldwat2m, CAM
121  lamcrad(i,j) = 0.0;
122  pgamrad(i,j) = 0.0;
123  }
124  });
125  } // m2005_effradius
126 
constexpr amrex::Real rhos
Definition: ERF_Constants.H:29
YAKL_INLINE void m2005_effradius(const real2d &ql, const real2d &nl, const real2d &qi, const real2d &ni, const real2d &qs, const real2d &ns, const real2d &cld, const real2d &pres, const real2d &tk, const real2d &effl, const real2d &effi, const real2d &deffi, const real2d &lamcrad, const real2d &pgamrad, const real2d &des)
Definition: ERF_m2005_effradius.H:25
@ pres
Definition: ERF_Kessler.H:33
@ rho
Definition: ERF_Kessler.H:30
Definition: ERF_console_io.cpp:12