.. _program_listing_file_src_modules_Richard_albedo.cpp:

Program Listing for File Richard_albedo.cpp
===========================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_modules_Richard_albedo.cpp>` (``src/modules/Richard_albedo.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   //
   // Canadian Hydrological Model - The Canadian Hydrological Model (CHM) is a novel
   // modular unstructured mesh based approach for hydrological modelling
   // Copyright (C) 2018 Christopher Marsh
   //
   // This file is part of Canadian Hydrological Model.
   //
   // Canadian Hydrological Model is free software: you can redistribute it and/or
   // modify
   // it under the terms of the GNU General Public License as published by
   // the Free Software Foundation, either version 3 of the License, or
   // (at your option) any later version.
   //
   // Canadian Hydrological Model is distributed in the hope that it will be useful,
   // but WITHOUT ANY WARRANTY; without even the implied warranty of
   // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   // GNU General Public License for more details.
   //
   // You should have received a copy of the GNU General Public License
   // along with Canadian Hydrological Model.  If not, see
   // <http://www.gnu.org/licenses/>.
   //
   
   #include "Richard_albedo.hpp"
   REGISTER_MODULE_CPP(Richard_albedo);
   
   Richard_albedo::Richard_albedo(config_file cfg)
   : module_base("Richard_albedo", parallel::data, cfg)
   {
   
       depends("swe");
       depends("T_s_0"); // snow temp
       depends("p_snow");
   
       provides("snow_albedo");
       provides("melting_albedo");
   
   }
   Richard_albedo::~Richard_albedo()
   {
   
   }
   
   void Richard_albedo::checkpoint(mesh& domain,  netcdf& chkpt)
   {
   
       chkpt.create_variable1D("Richard_albedo:albedo", domain->size_local_faces());
   
       for (size_t i = 0; i < domain->size_local_faces(); i++)
       {
           auto face = domain->face(i);
           chkpt.put_var1D("Richard_albedo:albedo",i,face->get_module_data<Richard_albedo::data>(ID).albedo);
       }
   }
   
   void Richard_albedo::load_checkpoint(mesh& domain,  netcdf& chkpt)
   {
       for (size_t i = 0; i < domain->size_local_faces(); i++)
       {
           auto face = domain->face(i);
           face->get_module_data<Richard_albedo::data>(ID).albedo = chkpt.get_var1D("Richard_albedo:albedo",i);
       }
   }
   
   void Richard_albedo::run(mesh_elem &face)
   {
       if(is_water(face))
       {
           set_all_nan_on_skip(face);
           return;
       }
       double albedo = face->get_module_data<Richard_albedo::data>(ID).albedo;
   
       double swe = (*face)["swe"_s];
   
       double dt = global_param->dt();
       if(global_param->first_time_step)
       {
           if(swe > 1.0)
           {
               albedo = albedo_snow;
           }
           else
           {
               albedo = albedo_bare;
           }
       }
   
       if ( swe > 0.)
       {
           if ((*face)["T_s_0"_s] >= 273.)  //melting snow, T_s_0???
           {
               albedo = (albedo - amin) * exp(-dt/a2) + amin;
               (*face)["melting_albedo"_s]=1;
           }
           else
           {
               albedo = albedo - dt/a1; // cold snow decay
               (*face)["melting_albedo"_s]=0;
           }
   
           double psnow = (*face)["p_snow"_s];
           albedo = albedo + (amax - albedo) * (psnow )/min_swe_refresh; //* dt
   
           albedo = std::max(albedo,amin);
           albedo = std::min(albedo,amax);
       }
       else
       {
           (*face)["melting_albedo"_s]=0;
           albedo = albedo_bare;
       }
   
   
       (*face)["snow_albedo"_s]=albedo;
       face->get_module_data<Richard_albedo::data>(ID).albedo = albedo;
   }
   
   void Richard_albedo::init(mesh& domain)
   {
   
       //these
       amin = cfg.get("albedo_min",0.5);
       amax = cfg.get("albedo_max",0.84);
       a1 = cfg.get("a1",1.08e7); //cold snow decay const
       a2 = cfg.get("a2",7.2e5);  //melting snow decay const
       min_swe_refresh = cfg.get("min_swe_refresh",1.); //min swe to refresh albedo
       albedo_snow = cfg.get("init_albedo_snow",0.85); // intial snow albedo
       albedo_bare = cfg.get("init_albedo_bare",0.17); //initial bare ground albedo
   
   
   #pragma omp parallel for
       for (size_t i=0; i < domain->size_local_faces(); ++i)
       {
   
           auto face = domain->face(i);
           auto& d = face->make_module_data<data>(ID);
           d.albedo = albedo_bare;
   
       }
   
   }