WIP add some fix to the aerosol implementation (still have issues) (#…

…1362) * enable mam4 aerosol model * fix tabs/whitespace * add erf_rad_data_dir * fix some bugs during the test, also add a lot of print output in this version * fix some aerosol radiation bugs
erf-model · Jan 4, 2024 · 0b90d43 · 0b90d43
1 parent aa250d0
commit 0b90d43
Show file tree

Hide file tree

Showing 32 changed files with 285,299 additions and 1,618 deletions.
diff --git a/Source/Radiation/Aero_rad_props.cpp b/Source/Radiation/Aero_rad_props.cpp
diff --git a/Source/Radiation/Cloud_rad_props.H b/Source/Radiation/Cloud_rad_props.H
@@ -74,8 +74,8 @@ class CloudRadProps {
                                LinInterp::InterpType& lambda_wgts);
   private:
     std::string name{"CloudRadProps"};
-    std::string liquid_file{"/ccs/home/yuanx/erf.radiation/Source/Radiation/data/F_nwvl200_mu20_lam50_res64_t298_c080428.nc"};
-    std::string ice_file{"/ccs/home/yuanx/erf.radiation/Source/Radiation/data/iceoptics_c080917.nc"};
+    std::string liquid_file;
+    std::string ice_file;
 
     int nlwbands, nswbands, nlambda, nmu, n_g_d;
 

diff --git a/Source/Radiation/Cloud_rad_props.cpp b/Source/Radiation/Cloud_rad_props.cpp
@@ -1,6 +1,7 @@
 //
 // get the cloud radiation props data
 //
+#include <cstdlib>
 #include "YAKL_netcdf.h"
 #include "Cloud_rad_props.H"
 using yakl::fortran::parallel_for;
@@ -20,6 +21,11 @@ void CloudRadProps::initialize() {
     realHost2d asm_sw_ice_h;
     realHost2d abs_lw_ice_h;
 
+    const char* erf_home = std::getenv("ERF_HOME");
+    auto erf_rad_data_dir = (std::string)erf_home + "/Source/Radiation/data/";
+    liquid_file = erf_rad_data_dir + "F_nwvl200_mu20_lam50_res64_t298_c080428.nc";
+    ice_file = erf_rad_data_dir + "iceoptics_c080917.nc";
+
     yakl::SimpleNetCDF liquid, ice;
     liquid.open(liquid_file, yakl::NETCDF_MODE_READ);
     ice.open(ice_file, yakl::NETCDF_MODE_READ);

diff --git a/Source/Radiation/Ebert_curry.H b/Source/Radiation/Ebert_curry.H
@@ -6,11 +6,12 @@ using yakl::fortran::SimpleBounds;
 
 class EbertCurry {
   public:
-  static constexpr real scalefactor = 1.;  //500._r8/917._r8
-
-  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) {
+   static constexpr real scalefactor = 1.;  //500._r8/917._r8
 
+   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) {
     real1d wavmin("wavmin", nswbands);
     real1d wavmax("wavmax", nswbands);
 
@@ -68,7 +69,6 @@ class EbertCurry {
     RadConstants::get_sw_spectral_boundaries(wavmin,wavmax,RadConstants::micrometer);
 
     for (auto ns=0; ns<nswbands; ++ns) {
-
        if(wavmax(ns) <= 0.7) {
           indxsl = 1;
        } else if(wavmax(ns) <= 1.25) {
@@ -78,7 +78,6 @@ class EbertCurry {
        } else if(wavmax(ns) > 2.38) {
           indxsl = 4;
        }
-
        parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i) {
           real tmp1i, tmp2i, tmp3i, g;
           auto abarii = abari(indxsl);
@@ -111,8 +110,9 @@ class EbertCurry {
   }
 
 
-  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) {
+  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) {
      real2d ficemr("ficemr",ncol,nlev);
      real2d cwp("cwp",ncol,nlev);
      real2d cldtau("cldtau",ncol,nlev);
@@ -142,7 +142,7 @@ class EbertCurry {
      });
 
      parallel_for(SimpleBounds<3>(nlwbands, ncol, nlev), YAKL_LAMBDA (int lwband, int icol, int ilev) {
-             abs_od(lwband,icol,ilev)=cldtau(icol,ilev);
+        abs_od(lwband,icol,ilev)=cldtau(icol,ilev);
      });
   }
 };

diff --git a/Source/Radiation/Linear_interpolate.H b/Source/Radiation/Linear_interpolate.H
@@ -48,7 +48,6 @@ class LinInterp
       real cmin, cmax;
       real dyinwrap;
       real avgdyin;
-
       //
       // Check validity of input coordinate arrays: must be monotonically increasing,
       // and have a total of at least 2 elements
@@ -71,9 +70,9 @@ class LinInterp
       interp_wgts.wgts = real1d("wgts", nout);
       interp_wgts.wgtn = real1d("wgtn", nout);
 
-      parallel_for(SimpleBounds<1>(nin), YAKL_LAMBDA (int j) {
+      parallel_for(SimpleBounds<1>(nin-1), YAKL_LAMBDA (int j) {
          if(yin(j) > yin(j+1)) {
-           printf("inputs are not monotonic!\n");
+           printf("init: inputs are not monotonic!\n");
            return;
          }
       });