C++ and Fortran Version with Legion and Kokkos

src/Makefile

@@ -5,13 +5,13 @@ OBJECT        = .
 default: all
 
 all:
-	$(MAKE) -C fortran all
+	$(MAKE) -C cc all
 
 clean:
-	$(MAKE) -C fortran clean
+	$(MAKE) -C cc clean
 
 allclean:
-	$(MAKE) -C fortran clean
+	$(MAKE) -C cc clean
 
 
 

src/cc/CanopyHydrology_FracH2OSfc.cc

+#include <algorithm>
+#include <stdio.h>     
+#include <math.h>   
+#include <iostream>
+#include <string>
+#include "landunit_varcon.h"
+#include "CanopyHydrology_cc.hh"
+
+using namespace std;
+
+namespace ELM {
+
+void CanopyHydrology_FracH2OSfc(const double& dtime,
+        const double& min_h2osfc,
+        const int& ltype,
+        const double& micro_sigma,
+        const double& h2osno,
+        double& h2osfc,
+        double& h2osoi_liq,
+        double& frac_sno,
+        double& frac_sno_eff,
+        double& qflx_h2osfc2topsoi,
+        double& frac_h2osfc)
+  {
+    bool no_update = false;
+    double shr_const_pi=4.0e0*atan(1.0e0) ;
+    bool no_update_l ;
+
+    
+    double d,fd,dfdd,sigma   ;
+
+    if (!no_update) { 
+      no_update_l = false; }
+    else { no_update_l = no_update; }
+    
+    qflx_h2osfc2topsoi = 0.0 ;
+    
+    if ( ltype  == istsoil || ltype == istcrop) {
+
+       
+
+       if (h2osfc > min_h2osfc) {
+          
+          d=0.0 ;
+
+          sigma=1.0e3 * micro_sigma ;
+          for(int l = 0 ; l < 10; l++) {
+             fd = 0.5*d*(1.00+erf(d/(sigma*sqrt(2.0)))) + sigma/sqrt(2.0*shr_const_pi)*std::exp(-std::pow(d,2)/(2.0*std::pow(sigma,2))) -h2osfc;
+             dfdd = 0.5*(1.00+erf(d/(sigma*sqrt(2.0))));
+
+             d = d - fd/dfdd;
+          }
+          
+          frac_h2osfc = 0.5*(1.00+erf(d/(sigma*sqrt(2.0)))) ;  }
+
+       else {
+          frac_h2osfc = 0.0 ;
+          h2osoi_liq = h2osoi_liq + h2osfc ;
+          qflx_h2osfc2topsoi = h2osfc/dtime ;
+          h2osfc=0.0 ;
+        }
+
+       if (!no_update_l) {
+
+          
+          if (frac_sno > (1.0 - frac_h2osfc) && h2osno > 0) {
+
+             if (frac_h2osfc > 0.010) {
+                frac_h2osfc = std::max(1.00 - frac_sno,0.010) ;
+                frac_sno = 1.00 - frac_h2osfc; }
+             else {
+                frac_sno = 1.00 - frac_h2osfc;
+              }
+             
+             frac_sno_eff=frac_sno;
+
+          }
+
+       } 
+    }  
+    else {
+
+       frac_h2osfc = 0.0;
+
+    }
+
+  }
+
+}
\ No newline at end of file

src/cc/CanopyHydrology_FracWet.cc

+#include <algorithm>
+#include <iostream>
+#include <cmath>
+#include <string>
+#include "CanopyHydrology_cc.hh"
+
+using namespace std;
+
+namespace ELM {
+
+void CanopyHydrology_FracWet(const int& frac_veg_nosno,
+        const double& h2ocan,
+        const double& elai, 
+        const double& esai,
+        const double& dewmx,
+        double& fwet,
+        double& fdry)
+{  
+
+  double vegt, dewmxi ;    
+
+  if (frac_veg_nosno == 1) {
+    if (h2ocan > 0.0) {
+        vegt    = frac_veg_nosno*(elai + esai);
+        dewmxi  = 1.00/dewmx;
+        fwet = std::pow(((dewmxi/vegt)*h2ocan), 2.0/3);
+        fwet = std::min(fwet,1.00);   
+        fdry = (1.0-fwet)*elai/(elai+esai);
+      }
+    else{
+        fwet = 0.0;
+        fdry = 0.0 ;
+      } }
+  else{
+     fwet = 0.0;
+     fdry = 0.0;
+  }
+
+}
+
+}
\ No newline at end of file

src/cc/CanopyHydrology_Interception.cc

+#include <algorithm>
+#include <stdio.h>     
+#include <math.h>   
+#include <iostream>
+#include <string>
+#include "landunit_varcon.h"
+#include "column_varcon.h"     
+#include "CanopyHydrology_cc.hh"
+
+using namespace std;
+
+namespace ELM {
+
+void CanopyHydrology_Interception(double dtime,
+        const double& forc_rain,
+        const double& forc_snow,
+        const double& irrig_rate,
+        const int& ltype, const int& ctype,
+        const bool& urbpoi, const bool& do_capsnow,
+        const double& elai, const double& esai,
+        const double& dewmx, const int& frac_veg_nosno,
+        double& h2ocan,
+        int& n_irrig_steps_left,
+        double& qflx_prec_intr,
+        double& qflx_irrig,
+        double& qflx_prec_grnd,
+        double& qflx_snwcp_liq,
+        double& qflx_snwcp_ice,
+        double& qflx_snow_grnd_patch,
+        double& qflx_rain_grnd)
+
+ {  
+
+  
+      double  fpi, xrun, h2ocanmx   ;
+      double  qflx_candrip, qflx_through_snow, qflx_through_rain ;
+      double  qflx_prec_grnd_snow;
+      double  qflx_prec_grnd_rain ;
+      double  fracsnow ;
+      double  fracrain , forc_irrig;
+
+
+      if (ltype==istsoil || ltype==istwet || urbpoi || ltype==istcrop) {
+
+         qflx_candrip = 0.0      ;
+         qflx_through_snow = 0.0 ;
+         qflx_through_rain = 0.0 ;
+         qflx_prec_intr = 0.0    ;
+         fracsnow = 0.0          ;
+         fracrain = 0.0          ;
+         forc_irrig = 0.0;
+
+
+         if (ctype != icol_sunwall && ctype != icol_shadewall) {
+            if (frac_veg_nosno == 1 && (forc_rain + forc_snow) > 0.0) {
+
+              
+               fracsnow = forc_snow/(forc_snow + forc_rain);
+               fracrain = forc_rain/(forc_snow + forc_rain);
+
+               
+               h2ocanmx = dewmx * (elai + esai);
+
+               
+               fpi = 0.250*(1.0 - exp(-0.50*(elai + esai)));
+
+              
+               qflx_through_snow = forc_snow * (1.0-fpi);
+               qflx_through_rain = forc_rain * (1.0-fpi);
+
+               
+               qflx_prec_intr = (forc_snow + forc_rain) * fpi;
+               
+
+
+               
+               h2ocan = std::max(0.0, h2ocan + dtime*qflx_prec_intr);
+
+               
+               qflx_candrip = 0.0;
+
+               
+               xrun = (h2ocan - h2ocanmx)/dtime;
+
+               
+               if (xrun > 0.0) {
+                  qflx_candrip = xrun;
+                  h2ocan = h2ocanmx;
+               }
+
+            }
+         }
+
+      else if (ltype==istice || ltype==istice_mec) {
+         
+         h2ocan            = 0.0;
+         qflx_candrip      = 0.0;
+         qflx_through_snow = 0.0;
+         qflx_through_rain = 0.0;
+         qflx_prec_intr    = 0.0;
+         fracsnow          = 0.0;
+         fracrain          = 0.0;
+
+      }
+
+      
+
+      if (ctype != icol_sunwall && ctype != icol_shadewall) {
+         if (frac_veg_nosno == 0) {
+            qflx_prec_grnd_snow = forc_snow;
+            qflx_prec_grnd_rain = forc_rain;  }
+         else{
+            qflx_prec_grnd_snow = qflx_through_snow + (qflx_candrip * fracsnow);
+            qflx_prec_grnd_rain = qflx_through_rain + (qflx_candrip * fracrain);
+          }
+      }   
+      else{
+         qflx_prec_grnd_snow = 0.;
+         qflx_prec_grnd_rain = 0.;
+        }
+
+      
+      if (n_irrig_steps_left > 0) {
+         qflx_irrig         = forc_irrig;
+         n_irrig_steps_left = n_irrig_steps_left - 1; }
+      else{
+         qflx_irrig = 0.0;
+        }
+
+      
+      qflx_prec_grnd_rain = qflx_prec_grnd_rain + qflx_irrig;
+
+      
+
+      qflx_prec_grnd = qflx_prec_grnd_snow + qflx_prec_grnd_rain;
+
+      if (do_capsnow) {
+         qflx_snwcp_liq = qflx_prec_grnd_rain;
+         qflx_snwcp_ice = qflx_prec_grnd_snow;
+
+         qflx_snow_grnd_patch = 0.0;
+         qflx_rain_grnd = 0.0;  }
+      else{
+
+         qflx_snwcp_liq = 0.0;
+         qflx_snwcp_ice = 0.0;
+         qflx_snow_grnd_patch = qflx_prec_grnd_snow   ;      //ice onto ground (mm/s)
+         qflx_rain_grnd     = qflx_prec_grnd_rain      ;   //liquid water onto ground (mm/s)
+        }
+
+    }
+  }
+
+}
\ No newline at end of file

src/cc/CanopyHydrology_SnowWater_impl.hh

+#include <stdio.h>  
+#include <math.h>     
+#include "landunit_varcon.h"    
+#include "column_varcon.h" 
+#include "clm_varpar.h"         
+#include "clm_varctl.h"     
+#include "CanopyHydrology_cc.hh"    
+
+using namespace std;
+
+namespace ELM {
+
+template<typename Array_d>
+void CanopyHydrology_SnowWater(const double& dtime,
+        const double& qflx_floodg,
+        const int& ltype,
+        const int& ctype,
+        const bool& urbpoi,
+        const bool& do_capsnow,                            
+        const int& oldfflag,
+        const double& forc_t,
+        const double& t_grnd,
+        const double& qflx_snow_grnd_col,
+        const double& qflx_snow_melt,
+        const double& n_melt,
+        const double& frac_h2osfc,
+        double& snow_depth,
+        double& h2osno,
+        double& int_snow,
+        Array_d swe_old,
+        Array_d h2osoi_liq,
+        Array_d h2osoi_ice,
+        Array_d t_soisno,
+        Array_d frac_iceold,
+        int& snow_level,
+        Array_d dz,
+        Array_d z,
+        Array_d zi,
+        int& newnode,
+        double& qflx_floodc,
+        double& qflx_snow_h2osfc,
+        double& frac_sno_eff,
+        double& frac_sno)
+{       
+  
+  
+//parameters
+  double rpi=4.0e0*atan(1.0e0)  ;
+  double tfrz=273.15;
+  double zlnd = 0.010;
+
+  
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: swe_old 
+  // real(r8), intent(inout), dimension(-nlevsno+1:0) :: h2osoi_liq, h2osoi_ice
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: t_soisno, frac_iceold
+  // real(r8), intent(inout), dimension(-nlevsno+1:0)  :: dz, z, zi
+  
+//local variables 
+  double  temp_intsnow, temp_snow_depth, z_avg, fmelt, dz_snowf, snowmelt ;
+  double  newsnow, bifall, accum_factor, fsno_new, smr ;
+  int j ;
+
+//apply gridcell flood water flux to non-lake columns
+  if (ctype != icol_sunwall && ctype != icol_shadewall) {      
+     qflx_floodc = qflx_floodg; }
+  else{
+     qflx_floodc = 0.0;
+  }
+
+//Determine snow height and snow water
+
+//Use Alta relationship, Anderson(1976); LaChapelle(1961),
+//U.S.Department of Agriculture Forest Service, Project F,
+//Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.
+
+  qflx_snow_h2osfc = 0.0;
+//set temporary variables prior to updating
+  temp_snow_depth=snow_depth;
+//save initial snow content
+  for(j = -nlevsno+1; j < snow_level; j++) {
+     swe_old[j] = 0.00;
+  }
+  for(j = snow_level+1; j < 0; j++) {
+     swe_old[j]=h2osoi_liq[j]+h2osoi_ice[j];
+  }
+
+  if (do_capsnow) {
+     dz_snowf = 0.;
+     newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
+     frac_sno=1.;
+     int_snow = 5.e2; 
+  } else {
+
+     if (forc_t > tfrz + 2.) {
+        bifall=50. + 1.7*std::pow((17.0),1.5);
+     } else if (forc_t > tfrz - 15.) {
+        bifall=50. + 1.7*std::pow((forc_t - tfrz + 15.),1.5);
+     } else {
+        bifall=50.;
+     }
+
+     // newsnow is all snow that doesn't fall on h2osfc
+     newsnow = (1. - frac_h2osfc) * qflx_snow_grnd_col * dtime;
+
+     // update int_snow
+     int_snow = std::max(int_snow,h2osno) ; //h2osno could be larger due to frost
+
+     // snowmelt from previous time step * dtime
+     snowmelt = qflx_snow_melt * dtime;
+
+     // set shape factor for accumulation of snow
+     accum_factor=0.1;
+
+     if (h2osno > 0.0) {
+
+        //======================  FSCA PARAMETERIZATIONS  ======================
+        // fsca parameterization based on *changes* in swe
+        // first compute change from melt during previous time step
+        if(snowmelt > 0.) {
+
+           smr=std::min(1.,(h2osno)/(int_snow));
+
+           frac_sno = 1. - std::pow((acos(min(1.,(2.*smr - 1.)))/rpi),(n_melt)) ;
+
+        }
+
+        // update fsca by new snow event, add to previous fsca
+        if (newsnow > 0.) {
+           fsno_new = 1. - (1. - tanh(accum_factor*newsnow))*(1. - frac_sno);
+           frac_sno = fsno_new;
+
+           // reset int_snow after accumulation events
+           temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*std::pow((1.0-std::max(frac_sno,1e-6)),(1.0/n_melt))+1.0))) ;
+           int_snow = std::min(1.e8,temp_intsnow) ;
+        }
+
+        //====================================================================
+
+        // for subgrid fluxes
+        if (subgridflag ==1 && ! urbpoi) {
+           if (frac_sno > 0.){
+              snow_depth=snow_depth + newsnow/(bifall * frac_sno);
+           } else {
+              snow_depth=0.;
+           }
+        } else {
+           // for uniform snow cover
+           snow_depth=snow_depth+newsnow/bifall;
+        }
+
+        // use original fsca formulation (n&y 07)
+        if (oldfflag == 1) { 
+           // snow cover fraction in Niu et al. 2007
+           if(snow_depth > 0.0)  {
+              frac_sno = tanh(snow_depth/(2.5*zlnd*std::pow((std::min(800.0,(h2osno+ newsnow)/snow_depth)/100.0),1.0)) ) ;
+           }
+           if(h2osno < 1.0)  {
+              frac_sno=std::min(frac_sno,h2osno);
+           }
+        }
+
+     } else { //h2osno == 0
+        // initialize frac_sno and snow_depth when no snow present initially
+        if (newsnow > 0.) { 
+           z_avg = newsnow/bifall;
+           fmelt=newsnow;
+           frac_sno = tanh(accum_factor*newsnow);
+
+           // make int_snow consistent w/ new fsno, h2osno
+           int_snow = 0. ;//reset prior to adding newsnow below
+           temp_intsnow= (h2osno + newsnow) / (0.5*(cos(rpi*std::pow((1.0-std::max(frac_sno,1e-6)),(1.0/n_melt)))+1.0));
+           int_snow = std::min(1.e8,temp_intsnow);
+
+           // update snow_depth and h2osno to be consistent with frac_sno, z_avg
+           if (subgridflag ==1 && !urbpoi) {
+              snow_depth=z_avg/frac_sno;
+           } else {
+              snow_depth=newsnow/bifall;
+           }
+           // use n&y07 formulation
+           if (oldfflag == 1) { 
+              // snow cover fraction in Niu et al. 2007
+              if(snow_depth > 0.0)  {
+                 frac_sno = tanh(snow_depth/(2.5*zlnd*std::pow((std::min(800.0,newsnow/snow_depth)/100.0),1.0)) );
+              }
+           }
+        } else {
+           z_avg = 0.;
+           snow_depth = 0.;
+           frac_sno = 0.;
+        }
+     } // end of h2osno > 0
+
+     // snow directly falling on surface water melts, increases h2osfc
+     qflx_snow_h2osfc = frac_h2osfc*qflx_snow_grnd_col;
+
+     // update h2osno for new snow
+     h2osno = h2osno + newsnow ;
+     int_snow = int_snow + newsnow;
+
+     // update change in snow depth
+     dz_snowf = (snow_depth - temp_snow_depth) / dtime;
+
+  } //end of do_capsnow construct
+
+  // set frac_sno_eff variable
+  if (ltype == istsoil || ltype == istcrop) {
+     if (subgridflag ==1) { 
+        frac_sno_eff = frac_sno;
+     } else {
+        frac_sno_eff = 1.;
+     }
+  } else {
+     frac_sno_eff = 1.;
+  }
+
+  if (ltype==istwet && t_grnd>tfrz) {
+     h2osno=0.;
+     snow_depth=0.;
+  }
+
+//When the snow accumulation exceeds 10 mm, initialize snow layer
+//Currently, the water temperature for the precipitation is simply set
+//as the surface air temperature
+  newnode = 0 ; //flag for when snow node will be initialized
+        if (snow_level == 0 && qflx_snow_grnd_col > 0.00 && frac_sno*snow_depth >= 0.010) {
+           newnode = 1;
+           snow_level = -1;
+           dz[0] = snow_depth ;                    //meter
+           z[0] = -0.50*dz[0];
+           zi[-1] = -dz[0];
+           t_soisno[0] = min(tfrz, forc_t) ;   //K
+           h2osoi_ice[0] = h2osno ;            //kg/m2
+           h2osoi_liq[0] = 0.0  ;               //kg/m2
+           frac_iceold[0] = 1.0;
+        }
+
+//The change of ice partial density of surface node due to precipitation.
+//Only ice part of snowfall is added here, the liquid part will be added
+//later.
+        if (snow_level < 0 && newnode == 0) {
+        h2osoi_ice[snow_level+1] = h2osoi_ice[snow_level+1]+newsnow;
+        dz[snow_level+1] = dz[snow_level+1]+dz_snowf*dtime;
+        }
+  }
+ }
+

src/cc/CanopyHydrology_cc.hh

+#ifndef ELM_CANOPY_HYDROLOGY_CC_HH_
+#define ELM_CANOPY_HYDROLOGY_CC_HH_
+
+namespace ELM {
+
+  void CanopyHydrology_Interception(double dtime,
+        const double& forc_rain,
+        const double& forc_snow,
+        const double& irrig_rate,
+        const int& ltype, const int& ctype,
+        const bool& urbpoi, const bool& do_capsnow,
+        const double& elai, const double& esai,
+        const double& dewmx, const int& frac_veg_nosno,
+        double& h2ocan,
+        int& n_irrig_steps_left,
+        double& qflx_prec_intr,
+        double& qflx_irrig,
+        double& qflx_prec_grnd,
+        double& qflx_snwcp_liq,
+        double& qflx_snwcp_ice,
+        double& qflx_snow_grnd_patch,
+        double& qflx_rain_grnd) ;
+
+
+
+ void CanopyHydrology_FracWet(const int& frac_veg_nosno,
+        const double& h2ocan,
+        const double& elai, 
+        const double& esai,
+        const double& dewmx,
+        double& fwet,
+        double& fdry) ;
+
+
+template<typename Array_d>
+void CanopyHydrology_SnowWater(const double& dtime,
+        const double& qflx_floodg,
+        const int& ltype,
+        const int& ctype,
+        const bool& urbpoi,
+        const bool& do_capsnow,                            
+        const int& oldfflag,
+        const double& forc_t,
+        const double& t_grnd,
+        const double& qflx_snow_grnd_col,
+        const double& qflx_snow_melt,
+        const double& n_melt,
+        const double& frac_h2osfc,
+        double& snow_depth,
+        double& h2osno,
+        double& int_snow,