ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_Ebert_curry.H
Go to the documentation of this file.
1 
2 #ifndef ERF_EBERT_CURRY_H_
3 #define ERF_EBERT_CURRY_H_
4 using yakl::fortran::parallel_for;
5 using yakl::fortran::SimpleBounds;
6 
7 class EbertCurry {
8  public:
9  static constexpr real scalefactor = 1.; //500._r8/917._r8
10 
11  static void ec_ice_optics_sw (int ncol, int nlev, int nswbands,
12  const real2d& cldn, const real2d& cicewp, const real2d& rei,
13  const real3d& ice_tau, const real3d& ice_tau_w,
14  const real3d& ice_tau_w_g, const real3d& ice_tau_w_f)
15  {
16  real1d wavmin("wavmin", nswbands);
17  real1d wavmax("wavmax", nswbands);
18 
19  // ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836)
20  real1d abari("abari", 4); // a coefficient for extinction optical depth
21  real1d bbari("bbari", 4); // b coefficient for extinction optical depth
22  real1d cbari("cbari", 4); // c coefficient for single scat albedo
23  real1d dbari("dbari", 4); // d coefficient for single scat albedo
24  real1d ebari("ebari", 4); // e coefficient for asymmetry parameter
25  real1d fbari("fbari", 4); // f coefficient for asymmetry parameter
26 
27  parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int i)
28  {
29  abari(1) = 3.448e-03;
30  abari(2) = 3.448e-03;
31  abari(3) = 3.448e-03;
32  abari(4) = 3.448e-03;
33  bbari(1) = 2.431;
34  bbari(2) = 2.431;
35  bbari(3) = 2.431;
36  bbari(4) = 2.431;
37  cbari(1) = 1.00e-05;
38  cbari(2) = 1.10e-04;
39  cbari(3) = 1.861e-02;
40  cbari(4) = 0.46658;
41  dbari(1) = 0.0;
42  dbari(2) = 1.405e-05;
43  dbari(3) = 8.328e-04;
44  dbari(4) = 2.05e-05;
45  ebari(1) = 0.7661;
46  ebari(2) = 0.7730;
47  ebari(3) = 0.794;
48  ebari(4) = 0.9595;
49  fbari(1) = 5.851e-04;
50  fbari(2) = 5.665e-04;
51  fbari(3) = 7.267e-04;
52  fbari(4) = 1.076e-04;
53  });
54 
55  // Minimum cloud amount (as a fraction of the grid-box area) to
56  // distinguish from clear sky
57  constexpr real cldmin = 1.0e-80;
58 
59  // Decimal precision of cloud amount (0 -> preserve full resolution;
60  // 10^-n -> preserve n digits of cloud amount)
61  constexpr real cldeps = 0.0;
62 
63  // Optical properties for ice are valid only in the range of
64  // 13 < rei < 130 micron (Ebert and Curry 92)
65  constexpr real rei_min = 13.;
66  constexpr real rei_max = 130.;
67 
68  // integer :: ns, i, k, indxsl
69  int indxsl;
70 
72 
73  for (auto ns=0; ns<nswbands; ++ns) {
74  if(wavmax(ns) <= 0.7) {
75  indxsl = 1;
76  } else if(wavmax(ns) <= 1.25) {
77  indxsl = 2;
78  } else if(wavmax(ns) <= 2.38) {
79  indxsl = 3;
80  } else if(wavmax(ns) > 2.38) {
81  indxsl = 4;
82  }
83  parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i)
84  {
85  real tmp1i, tmp2i, tmp3i, g;
86  auto abarii = abari(indxsl);
87  auto bbarii = bbari(indxsl);
88  auto cbarii = cbari(indxsl);
89  auto dbarii = dbari(indxsl);
90  auto ebarii = ebari(indxsl);
91  auto fbarii = fbari(indxsl);
92 
93  // note that optical properties for ice valid only
94  // in range of 13 > rei > 130 micron (Ebert and Curry 92)
95  if (cldn(i,k) >= cldmin && cldn(i,k) >= cldeps) {
96  tmp1i = abarii + bbarii/std::max(rei_min,std::min(scalefactor*rei(i,k),rei_max));
97  ice_tau(ns,i,k) = 1000. * cicewp(i,k) * tmp1i;
98  } else {
99  ice_tau(ns,i,k) = 0.0;
100  }
101 
102  tmp2i = 1. - cbarii - dbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max);
103  tmp3i = fbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max);
104  // Do not let single scatter albedo be 1. Delta-eddington solution
105  // for non-conservative case has different analytic form from solution
106  // for conservative case, and raddedmx is written for non-conservative case.
107  ice_tau_w(ns,i,k) = ice_tau(ns,i,k) * std::min(tmp2i,.999999);
108  g = ebarii + tmp3i;
109  ice_tau_w_g(ns,i,k) = ice_tau_w(ns,i,k) * g;
110  ice_tau_w_f(ns,i,k) = ice_tau_w(ns,i,k) * g * g;
111  });
112  } // nswbands
113  }
114 
115 
116  static void ec_ice_optics_lw (int ncol, int nlev, int nlwbands,
117  const real2d& cldn, const real2d& iclwpth,
118  const real2d& iciwpth, const real2d& rei, const real3d& abs_od)
119  {
120  real2d ficemr("ficemr",ncol,nlev);
121  real2d cwp("cwp",ncol,nlev);
122  real2d cldtau("cldtau",ncol,nlev);
123 
124  // longwave liquid absorption coeff (m**2/g)
125  //const real kabsl = 0.090361;
126 
127  // Optical properties for ice are valid only in the range of
128  // 13 < rei < 130 micron (Ebert and Curry 92)
129  const real rei_min = 13.;
130  const real rei_max = 130.;
131 
132  parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i)
133  {
134  cwp(i,k) = 1000.0 *iciwpth(i,k) + 1000.0 *iclwpth(i,k);
135  ficemr(i,k) = 1000.0*iciwpth(i,k)/(std::max(1.e-18,cwp(i,k)));
136  });
137 
138  parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i)
139  {
140  // Note from Andrew Conley:
141  // Optics for RK no longer supported, This is constructed to get
142  // close to bit for bit. Otherwise we could simply use ice water path
143  // note that optical properties for ice valid only
144  // in range of 13 > rei > 130 micron (Ebert and Curry 92)
145  auto kabsi = 0.005 + 1./std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max);
146  auto kabs = kabsi*ficemr(i,k); // kabsl*(1.-ficemr(i,k)) + kabsi*ficemr(i,k);
147  cldtau(i,k) = kabs*cwp(i,k);
148  });
149 
150  parallel_for(SimpleBounds<3>(nlwbands, ncol, nlev), YAKL_LAMBDA (int lwband, int icol, int ilev)
151  {
152  abs_od(lwband,icol,ilev)=cldtau(icol,ilev);
153  });
154  }
155 };
156 #endif
157 
Definition: ERF_Ebert_curry.H:7
static constexpr real scalefactor
Definition: ERF_Ebert_curry.H:9
static void ec_ice_optics_lw(int ncol, int nlev, int nlwbands, const real2d &cldn, const real2d &iclwpth, const real2d &iciwpth, const real2d &rei, const real3d &abs_od)
Definition: ERF_Ebert_curry.H:116
static void ec_ice_optics_sw(int ncol, int nlev, int nswbands, const real2d &cldn, const real2d &cicewp, const real2d &rei, const real3d &ice_tau, const real3d &ice_tau_w, const real3d &ice_tau_w_g, const real3d &ice_tau_w_f)
Definition: ERF_Ebert_curry.H:11
@ micrometer
Definition: ERF_Rad_constants.H:109
static void get_sw_spectral_boundaries(real1d &low_boundaries, real1d &high_boundaries, Units units)
Definition: ERF_Rad_constants.H:205