Merge pull request #1311 from xyuan/xyuan/rad_develop

fix some issues during test

Exec/Radiation/inputs_radiation

@@ -46,7 +46,6 @@ erf.use_coriolis = false
 erf.use_rayleigh_damping = false
 
 erf.moisture_model = "SAM"
-
 erf.les_type = "Deardorff"
 #erf.les_type = "None"
 #

Source/Radiation/Albedo.cpp

@@ -7,7 +7,10 @@ void set_albedo(const real1d& coszrs, real2d& albedo_dir, real2d& albedo_dif)
   // Albedos for land type I (Briegleb)
   auto nswbands = albedo_dir.extent(0);
   auto ncol     = albedo_dif.extent(1);
-  yakl::c::parallel_for(yakl::c::Bounds<2>(nswbands,ncol), YAKL_LAMBDA (int ibnd, int icol) {
+  using yakl::fortran::parallel_for;
+  using yakl::fortran::SimpleBounds;
+
+  parallel_for(SimpleBounds<2>(nswbands,ncol), YAKL_LAMBDA (int ibnd, int icol) {
     albedo_dir(ibnd, icol) = 1.4 * 0.24 / ( 1. + 0.8 * coszrs(icol));
     albedo_dif(ibnd, icol) = 1.2 * 0.24;
   });

Source/Radiation/Cloud_rad_props.H

@@ -77,8 +77,7 @@ class 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"};
 
-    int nlw_band, nsw_band, nmu, nlambda;
-    int nlwbands, nswbands, n_g_d;
+    int nlwbands, nswbands, nlambda, nmu, n_g_d;
 
     real1d g_mu;         // mu samples on grid
     real2d g_lambda;     // lambda scale samples on grid

Source/Radiation/Cloud_rad_props.cpp

@@ -3,6 +3,8 @@
 //
 #include "YAKL_netcdf.h"
 #include "Cloud_rad_props.H"
+using yakl::fortran::parallel_for;
+using yakl::fortran::SimpleBounds;
 
 void CloudRadProps::initialize() {
     realHost1d g_mu_h;         // mu samples on grid
@@ -23,9 +25,9 @@ void CloudRadProps::initialize() {
     ice.open(ice_file, yakl::NETCDF_MODE_READ);
 
     // read the dimensions
-    nlw_band = liquid.getDimSize( "lw_band" );
-    nsw_band = liquid.getDimSize( "sw_band" );
-    nmu      = liquid.getDimSize( "nmu"     );
+    nlwbands = liquid.getDimSize( "lw_band" );
+    nswbands = liquid.getDimSize( "sw_band" );
+    nmu      = liquid.getDimSize( "mu"     );
     nlambda  = liquid.getDimSize( "lambda_scale" );
 
     liquid.read( g_mu_h, "mu");
@@ -35,6 +37,13 @@ void CloudRadProps::initialize() {
     liquid.read( asm_sw_liq_h, "asm_sw");
     liquid.read( abs_lw_liq_h, "k_abs_lw");
 
+    g_mu = real1d("gmu", nmu);
+    g_lambda = real2d("g_lambda", nmu, nlambda);
+    ext_sw_liq = real3d("ext_sw_liq", nmu, nlambda, nswbands);
+    ssa_sw_liq = real3d("ssa_sw_liq", nmu, nlambda, nswbands);
+    asm_sw_liq = real3d("asm_sw_liq", nmu, nlambda, nswbands);
+    abs_lw_liq = real3d("abs_lw_liq", nmu, nlambda, nlwbands);
+
     g_mu_h.deep_copy_to(g_mu);
     g_lambda_h.deep_copy_to(g_lambda);
     ext_sw_liq_h.deep_copy_to(ext_sw_liq);
@@ -43,11 +52,11 @@ void CloudRadProps::initialize() {
     abs_lw_liq_h.deep_copy_to(abs_lw_liq);
 
    // I forgot to convert kext from m^2/Volume to m^2/Kg
-   yakl::c::parallel_for(yakl::c::Bounds<3>(nmu,nlambda,nsw_band) , YAKL_LAMBDA (int i, int j, int k) {
+   parallel_for(SimpleBounds<3>(nmu,nlambda,nswbands) , YAKL_LAMBDA (int i, int j, int k) {
      ext_sw_liq(i,j,k) = ext_sw_liq(i,j,k) / 0.9970449e3;
    });
 
-   yakl::c::parallel_for(yakl::c::Bounds<3>(nmu,nlambda,nlw_band) , YAKL_LAMBDA (int i, int j, int k) {
+   parallel_for(SimpleBounds<3>(nmu,nlambda,nlwbands) , YAKL_LAMBDA (int i, int j, int k) {
      abs_lw_liq(i,j,k) = abs_lw_liq(i,j,k) / 0.9970449e3;
    });
 
@@ -62,6 +71,12 @@ void CloudRadProps::initialize() {
    ice.read( asm_sw_ice_h, "sw_asm");
    ice.read( abs_lw_ice_h, "lw_abs");
 
+   g_d_eff = real1d("g_d_eff",n_g_d);
+   ext_sw_ice = real2d("ext_sw_ice", n_g_d, nswbands);
+   ssa_sw_ice = real2d("ssa_sw_ice", n_g_d, nswbands);
+   asm_sw_ice = real2d("asm_sw_ice", n_g_d, nswbands);
+   abs_lw_ice = real2d("abs_lw_ice", n_g_d, nlwbands);
+
    g_d_eff_h.deep_copy_to(g_d_eff);
    ext_sw_ice_h.deep_copy_to(ext_sw_ice);
    ssa_sw_ice_h.deep_copy_to(ssa_sw_ice);
@@ -83,9 +98,9 @@ void CloudRadProps::gammadist_liq_optics_sw(const int& ncol,
   real1d tauw1d("tauw1d", nswbands);
   real1d tauwg1d("tauwg1d", nswbands);
   real1d tauwf1d("tauwf1d", nswbands);
-  yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
+  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
      if(lamc(i,k) > 0.) { // This seems to be clue from microphysics of no cloud
-        for (auto iband=0; iband<nswbands; ++iband) {
+        for (auto iband=1; iband<=nswbands; ++iband) {
           tau1d(iband)   = tau(iband,i,k);
           tauw1d(iband)  = tau_w(iband,i,k);
           tauwg1d(iband) = tau_w_g(iband,i,k);
@@ -94,7 +109,7 @@ void CloudRadProps::gammadist_liq_optics_sw(const int& ncol,
         gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), tau1d, tauw1d, tauwg1d, tauwf1d);
      }
      else {
-       for (auto iband=0; iband<nswbands; ++iband) {
+       for (auto iband=1; iband<=nswbands; ++iband) {
          tau(iband,i,k) = 0.;
          tau_w(iband,i,k) = 0.;
          tau_w_g(iband,i,k) = 0.;
@@ -111,13 +126,13 @@ void CloudRadProps::gammadist_liq_optics_lw(const int& ncol,
                                             const real2d& pgam,
                                             real3d& abs_od) {
    auto abs_od_1d = real1d("abs_od_1d", nlwbands);
-   yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
+   parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
       if(lamc(i,k) > 0.) { // This seems to be the clue for no cloud from microphysics formulation
-        for (auto ib=0; ib<nlwbands; ++ib) abs_od_1d(ib) = abs_od(ib,i,k);
+        for (auto ib=1; ib<=nlwbands; ++ib) abs_od_1d(ib) = abs_od(ib,i,k);
         gam_liquid_lw(iclwpth(i,k), lamc(i,k), pgam(i,k), abs_od_1d);
       }
       else {
-        for(auto j=0; j<nlwbands; ++j) abs_od(j,i,k) = 0.;
+        for(auto j=1; j<=nlwbands; ++j) abs_od(j,i,k) = 0.;
      }
   });
 }
@@ -131,8 +146,6 @@ void CloudRadProps::mitchell_ice_optics_sw(const int& ncol,
                                            real3d& tau_w,
                                            real3d& tau_w_g,
                                            real3d& tau_w_f) {
-  LinInterp::InterpType dei_wgts;
-
   real1d ext("ext", nswbands),
          ssa("ssa", nswbands),
          assm("assm", nswbands);
@@ -142,14 +155,14 @@ void CloudRadProps::mitchell_ice_optics_sw(const int& ncol,
          asm_sw_ice_1d("asm_sw_ice_1d", n_g_d);
 
   real1d dei_1d("dei_1d",1);
-  real1d ext_1d("ext_1d",1);
-  real1d ssa_1d("ssa_1d",1);
-  real1d assm_1d("assm_1d",1);
+  auto ext_1d  = real1d("ext_1d",1);
+  auto ssa_1d  = real1d("ssa_1d",1);
+  auto assm_1d = real1d("assm_1d",1);
 
-  yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
+  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
      if( iciwpth(i,k) < 1.e-80 || dei(i,k) == 0.) {
         // if ice water path is too small
-        for(auto ib=0; ib<nswbands; ++ib) {
+        for(auto ib=1; ib<=nswbands; ++ib) {
           tau    (ib,i,k) = 0.;
           tau_w  (ib,i,k) = 0.;
           tau_w_g(ib,i,k) = 0.;
@@ -158,23 +171,24 @@ void CloudRadProps::mitchell_ice_optics_sw(const int& ncol,
      }
      else {
        // for each cell interpolate to find weights in g_d_eff grid.
-       dei_1d(0) = dei(i,k);
+       dei_1d(1) = dei(i,k);
+       LinInterp::InterpType dei_wgts;
        LinInterp::init(g_d_eff, n_g_d, dei_1d, 1, LinInterp::extrap_method_bndry, dei_wgts);
        // interpolate into grid and extract radiative properties
-       for (auto is=0; is<nswbands; ++is) {
-          for(auto ig=0; ig<n_g_d; ++ig) {
+       for (auto is=1; is<=nswbands; ++is) {
+          for(auto ig=1; ig<=n_g_d; ++ig) {
              ext_sw_ice_1d(ig) = ext_sw_ice(ig, is);
              ssa_sw_ice_1d(ig) = ssa_sw_ice(ig, is);
              asm_sw_ice_1d(ig) = asm_sw_ice(ig, is);
           }
-          ext_1d(0) = ext(is);
-          ssa_1d(0) = ssa(is);
-          assm_1d(0) = assm(is);
-          LinInterp::interp1d(ext_sw_ice_1d, n_g_d, ext_1d, 1, dei_wgts);
-          LinInterp::interp1d(ssa_sw_ice_1d, n_g_d, ssa_1d, 1, dei_wgts);
+          ext_1d(1)  = ext(is);
+          ssa_1d(1)  = ssa(is);
+          assm_1d(1) = assm(is);
+          LinInterp::interp1d(ext_sw_ice_1d, n_g_d, ext_1d,  1, dei_wgts);
+          LinInterp::interp1d(ssa_sw_ice_1d, n_g_d, ssa_1d,  1, dei_wgts);
           LinInterp::interp1d(asm_sw_ice_1d, n_g_d, assm_1d, 1, dei_wgts);
        }
-       for (auto is=0; is<nswbands; ++is) {
+       for (auto is=1; is<=nswbands; ++is) {
          tau    (is,i,k) = iciwpth(i,k) * ext(is);
          tau_w  (is,i,k) = tau(is,i,k) * ssa(is);
          tau_w_g(is,i,k) = tau_w(is,i,k) * assm(is);
@@ -189,38 +203,35 @@ void CloudRadProps::mitchell_ice_optics_lw(const int& ncol,
                                            const real2d& iciwpth,
                                            const real2d& dei,
                                            real3d& abs_od) {
-  LinInterp::InterpType dei_wgts;
-
   real1d absor("absor", nlwbands);
   real1d dei_1d("dei_1d",1);
   real1d abs_lw_ice_1d("abs_lw_ice_1d",n_g_d);
 
   real1d absor_1d("absor_1d",1);
 
-  for(auto i=0; i<ncol; ++i) {
-    for(auto k=0; k<nlev; ++k) {
+  parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) {
      // if ice water path is too small, OD := 0
       if(iciwpth(i,k) < 1.e-80 || dei(i,k) == 0.) {
-        for (auto lb=0; lb<nlwbands; ++lb) {
+        for (auto lb=1; lb<=nlwbands; ++lb) {
           abs_od (lb,i,k) = 0.;
         }
       }
       else {
         // for each cell interpolate to find weights in g_d_eff grid.
-        dei_1d(0) = dei(i,k);
+        dei_1d(1) = dei(i,k);
+        LinInterp::InterpType dei_wgts;
         LinInterp::init(g_d_eff, n_g_d, dei_1d, 1, LinInterp::extrap_method_bndry, dei_wgts);
         // interpolate into grid and extract radiative properties
-        for(auto lwband = 0; lwband < nlwbands; ++lwband) {
-          absor_1d(0) = absor(lwband);
+        for(auto lwband = 1; lwband <= nlwbands; ++lwband) {
+          absor_1d(1) = absor(lwband);
           LinInterp::interp1d(abs_lw_ice_1d, n_g_d, absor_1d, 1, dei_wgts);
-          absor(lwband) = absor_1d(0);
+          absor(lwband) = absor_1d(1);
         }
-        for (auto lwband = 0; lwband < nlwbands; ++lwband) {
+        for (auto lwband = 1; lwband <= nlwbands; ++lwband) {
            abs_od(lwband,i,k) = iciwpth(i,k) * absor(lwband);
         }
      }
-   }
- }
+ });
 }
 
 
@@ -241,18 +252,18 @@ void CloudRadProps::gam_liquid_lw(const real& clwptn,
 
   get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts);
 
-  for(auto lwband = 1; lwband < nlwbands; ++lwband) {
-    for (auto imu=0; imu<nmu; ++imu) {
-       for (auto ilb=0; ilb<nlambda; ++ilb) {
+  parallel_for(SimpleBounds<1>(nlwbands), YAKL_LAMBDA (int lwband) {
+    for (auto imu=1; imu<=nmu; ++imu) {
+       for (auto ilb=1; ilb<=nlambda; ++ilb) {
            abs_lw_liq_2d(imu, ilb) = abs_lw_liq(imu, ilb, lwband);
        }
     }
-    abs_od_1d(0) = abs_od(lwband);
+    abs_od_1d(1) = abs_od(lwband);
     LinInterp::interp2d1d(abs_lw_liq_2d, nmu, nlambda, abs_od_1d, 1, mu_wgts, lambda_wgts);
-    abs_od(lwband) = clwptn*abs_od_1d(0);
-  }
+    abs_od(lwband) = clwptn*abs_od_1d(1);
+ });
 
-  yakl::c::parallel_for(yakl::c::Bounds<1>(nlwbands) , YAKL_LAMBDA (int iband) {
+  parallel_for(SimpleBounds<1>(nlwbands) , YAKL_LAMBDA (int iband) {
     abs_od(iband) = clwptn * abs_od(iband);
   });
 }
@@ -285,9 +296,9 @@ void CloudRadProps::gam_liquid_sw(const real& clwptn,
 
   get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts);
 
-  for(auto swband = 1; swband < nswbands; ++swband) {
-     for (auto imu=0; imu<nmu; ++imu) {
-        for (auto lb=0; lb<nlambda; ++lb) {
+  parallel_for(SimpleBounds<1>(nswbands), YAKL_LAMBDA (int swband) {
+     for (auto imu=1; imu<=nmu; ++imu) {
+        for (auto lb=1; lb<=nlambda; ++lb) {
            ext_sw_liq_2d(imu,lb) = ext_sw_liq(imu,lb,swband);
            ssa_sw_liq_2d(imu,lb) = ssa_sw_liq(imu,lb,swband);
            asm_sw_liq_2d(imu,lb) = asm_sw_liq(imu,lb,swband);
@@ -296,10 +307,10 @@ void CloudRadProps::gam_liquid_sw(const real& clwptn,
      LinInterp::interp2d1d(ext_sw_liq_2d, nmu, nlambda, extb, 1, mu_wgts, lambda_wgts);
      LinInterp::interp2d1d(ssa_sw_liq_2d, nmu, nlambda, ssab, 1, mu_wgts, lambda_wgts);
      LinInterp::interp2d1d(asm_sw_liq_2d, nmu, nlambda, asmb, 1, mu_wgts, lambda_wgts);
-  }
+  });
 
   // compute radiative properties
-  yakl::c::parallel_for(yakl::c::Bounds<1>(nswbands), YAKL_LAMBDA (int iband) {
+  parallel_for(SimpleBounds<1>(nswbands), YAKL_LAMBDA (int iband) {
     tau(iband)     = clwptn * extb(iband);
     tau_w(iband)   = tau(iband) * ssab(iband);
     tau_w_g(iband) = tau_w(iband) * asmb(iband);
@@ -327,10 +338,10 @@ void CloudRadProps::get_mu_lambda_weights(const real& lamc,
   LinInterp::init(g_mu0, nmu, pgam1d, 1, LinInterp::extrap_method_bndry, mu_wgts);
 
   // Use mu weights to interpolate to a row in the lambda table.
-  for(auto i=0; i<nlambda; ++i) {
-     for (auto im=0; im<nmu; ++im) g_lambda1d(im) = g_lambda(im, i);
+  parallel_for(SimpleBounds<1>(nlambda), YAKL_LAMBDA (int i) {
+     for (auto im=1; im<=nmu; ++im) g_lambda1d(im) = g_lambda(im, i);
      LinInterp::interp1d(g_lambda1d, nmu, g_lambda_interp, 1, mu_wgts);
-  }
+  });
 
   // Make interpolation weights for lambda.
   LinInterp::init(g_lambda_interp, nlambda, lamc1d, 1, LinInterp::extrap_method_bndry, lambda_wgts);

Source/Radiation/Ebert_curry.H

@@ -1,6 +1,8 @@
 
 #ifndef ERF_EBERT_CURRY_H_
 #define ERF_EBERT_CURRY_H_
+using yakl::fortran::parallel_for;
+using yakl::fortran::SimpleBounds;
 
 class EbertCurry {
   public:

Source/Radiation/Init_rrtmgp.cpp

@@ -10,15 +10,10 @@
 // Rrtmgp
 #include "Rrtmgp.H"
 
-void Rrtmgp::initialize(int num_gas, std::vector<std::string> active_gas_names,
+void Rrtmgp::initialize(int num_gas, const std::vector<std::string>& active_gas_names,
                         const char* rrtmgp_coefficients_file_sw,
                         const char* rrtmgp_coefficients_file_lw)
 {
-    // First, make sure yakl has been initialized
-    if (!yakl::isInitialized()) {
-        yakl::init();
-    }
-
     // Read gas optics coefficients from file
     // Need to initialize available_gases here! The only field of the
     // available_gases type that is used in the kdist initialize is
@@ -28,17 +23,14 @@ void Rrtmgp::initialize(int num_gas, std::vector<std::string> active_gas_names,
     // impossible from this initialization routine because I do not think the
     // rad_cnst objects are setup yet.
     // the other tasks!
-    ngas      = num_gas;
-    for (auto i = 0; i < ngas; ++i)
-       gas_names[i] = active_gas_names[i].c_str();
-
+    ngas  = num_gas;
     coefficients_file_sw = rrtmgp_coefficients_file_sw;
     coefficients_file_lw = rrtmgp_coefficients_file_lw;
 
-    auto active_gases = string1d("active_gases", ngas);
-    for (int igas=0; igas<ngas; igas++) {
-        active_gases(igas+1) = gas_names[igas];
-    }
+    active_gases = string1d("active_gases", ngas);
+    for (int igas=0; igas<ngas; igas++)
+        active_gases(igas+1) = active_gas_names[igas];
+
     GasConcs available_gases;
     available_gases.init(active_gases, 1, 1);
     load_and_init(k_dist_sw, coefficients_file_sw, available_gases);