<commit_before>/*
  Functions for calculating the snow / surface energy balance.
*/

#include <cmath>
#include <algorithm>
#include <boost/math/tools/roots.hpp>

#include "dbc.hh"

#include "seb_physics_funcs.hh"

namespace Amanzi {
namespace SurfaceBalance {
namespace SEBPhysics {

#define SWE_EPS 1.e-12

void UpdateIncomingRadiation(const SEB& seb, EnergyBalance& eb, bool debug) {
  // Calculate incoming short-wave radiation
  eb.fQswIn = (1 - seb.in.surf.albedo) * seb.in.met.QswIn;

  // Calculate incoming long-wave radiation
  const ThermoProperties& vp_air = seb.in.met.vp_air;
  double e_air = std::pow(10*vp_air.actual_vaporpressure, vp_air.temp / 2016.);
  e_air = 1.08 * (1 - std::exp(-e_air));
  eb.fQlwIn = e_air * seb.params.stephB * std::pow(vp_air.temp,4);

  // Calculate D_h, D_e, 
  eb.Dhe = std::pow(seb.params.VKc,2) * seb.in.met.Us
                       / std::pow(std::log(seb.params.Zr / seb.in.surf.Zo), 2);

  if (debug) {
    std::cout << "Incoming Radiation Energy Terms:" << "\n"
              << "  windspeed, Zo: " << seb.in.met.Us <<"  "<<seb.in.surf.Zo << "\n"
              << "  fQswIn   = " << eb.fQswIn << "\n"
              << "  fQlwIn   = " << eb.fQlwIn << std::endl;
  }

}

void UpdateEvapResistance(const SEB& seb, const ThermoProperties& vp_surf, EnergyBalance& eb, bool debug) {
   const ThermoProperties& vp_air = seb.in.met.vp_air;
   
   double Vaper_direction = vp_air.actual_vaporpressure - vp_surf.actual_vaporpressure;  //
   
// Equation for reduced vapor diffusivity See Sakagucki and Zeng 2009 eqaution (9) and Moldrup et al., 2004. 
   double Clab_Horn_b = 1;
   double actual_porosity = 0.9;  // Hard coded for moss Fix this to pass in form ATS ~AA
   double Surface_Vap_Diffusion = std::pow((1-(0.0556/actual_porosity)),(2+3*Clab_Horn_b));
   Surface_Vap_Diffusion = 0.000022 * (std::pow(actual_porosity,2)) * Surface_Vap_Diffusion;
// Sakagucki and Zeng 2009 eqaution (10)
   double cell_dimension = 0.01/2; // This is from cell center to the boundary.
   double VWC = seb.in.surf.saturation_liquid * actual_porosity;
   double L_Rsoil = std::exp(std::pow((1-(VWC/actual_porosity)),5));
   L_Rsoil = cell_dimension * (L_Rsoil -1) * (1/(2.718-1));
   double Rsoil = 0.0;
   
   if(Vaper_direction <= 0){
      Rsoil = L_Rsoil/Surface_Vap_Diffusion;
      }else{
      Rsoil = 0.0;
   }
   double Rair = 1/eb.Dhe;
   eb.Evap_Resistance = Rair;

   if((seb.in.snow_old.ht==0)&&(seb.in.surf.saturation_liquid < 1)) {
      eb.Evap_Resistance = Rair + Rsoil;
   }else{  
      eb.Evap_Resistance = Rair;
   }

   if (debug) {
   std::cout<<"Snow ht: "<<seb.out.snow_new.ht<<"  Saturation: "<<seb.in.surf.saturation_liquid<<std::endl;
   std::cout<<"Air_Temp: "<<vp_air.temp<<" Ground Temp: "<<vp_surf.temp<< std::endl;
   std::cout<<"Vapor Direction: "<<Vaper_direction<<" Air Vapor Pres: "<<vp_air.actual_vaporpressure<<" Surf Vapor Pres: "<<vp_surf.actual_vaporpressure << std::endl;
   std::cout<<"Surface_Vap_Diffusion: "<< Surface_Vap_Diffusion <<"  L_Rsoil:  "<<L_Rsoil<<"  Rsoil: "<<Rsoil<<"   VWC: "<<VWC<<std::endl;
   std::cout<<"Air Resistance: "<<Rair<<std::endl;
   std::cout<<"Evap Resistance: "<<eb.Evap_Resistance<<std::endl; 
   }
}


void UpdateEnergyBalance(const SEB& seb, const ThermoProperties& vp_surf, EnergyBalance& eb, bool debug) {
  const ThermoProperties& vp_air = seb.in.met.vp_air;

  // Calculate outgoing long-wave radiation
  eb.fQlwOut = -seb.in.surf.emissivity*seb.params.stephB*std::pow(vp_surf.temp,4);

  double Sqig;          // Stability function for use in e and h, precalculated for efficiency
  double air_temp = seb.in.met.vp_air.temp;
  double Ri  = seb.params.gravity * seb.params.Zr * (air_temp - vp_surf.temp)
      / (air_temp * std::pow(seb.in.met.Us,2));
  if (Ri >= 0.) {
    // stable condition or snow?
    Sqig = 1 / (1 + 10*Ri);
  } else {
    // Unstable condition
    Sqig = (1-10*Ri);
  }

  // Calculate sensible heat flux
  eb.fQh = seb.params.density_air * seb.params.Cp * eb.Dhe * Sqig * (vp_air.temp - vp_surf.temp);

  // Calculate latent heat flux
  double LatenHeatOf = 0.0; 
  if (seb.in.snow_old.ht > 0.){
    LatenHeatOf = seb.params.Ls;
   }else{
    LatenHeatOf = seb.params.Le;
  } 

  // With New Evaporation resistance term
  eb.fQe = vp_surf.porosity * seb.params.density_air * LatenHeatOf * (1/eb.Evap_Resistance) * Sqig * 0.622
      * (vp_air.actual_vaporpressure - vp_surf.actual_vaporpressure) / seb.params.Apa;
  if (debug) {
    std::cout<<"Porosity: "<<vp_surf.porosity<<"  LatentHeatOf: "<<LatenHeatOf<<"  Sqig: "<<Sqig<<std::endl;
    std::cout<<"Evap_Resistance: "<<eb.Evap_Resistance<<" Air vapor pres: "<<vp_air.actual_vaporpressure<<" Snow vapor pres: "<<vp_surf.actual_vaporpressure<<std::endl;
  }

  // Calculate heat conducted to ground, if snow
  if (seb.in.snow_old.ht > 0.) {
    double Ks = 2.9e-6 * std::pow(seb.in.snow_old.density,2);
    eb.fQc = Ks * (vp_surf.temp - seb.in.vp_ground.temp) / seb.in.snow_old.ht;
  }


  if (debug) {
    std::cout << "Energy Balance Terms (ht_snow = " << seb.in.snow_old.ht << "):" << "\n"
              << "  SnowSurfaceTemp  = " << vp_surf.temp << "\n"
              << "  fQlwOut  = " << eb.fQlwOut << "\n"
              << "  fQh      = " << eb.fQh << "\n"
              << "  fQe      = " << eb.fQe << "\n"
              << "  fQc      = " << eb.fQc << std::endl;
  }

}


void UpdateMassBalance(const SEB& seb, MassBalance& mb, EnergyBalance& eb, SnowProperties& snow_new, bool debug) {
  // this dt is the max timestep that may be taken to conserve snow mass
  mb.dt = seb.in.dt;

  // Melt rate given by available energy rate divided by heat of fusion.
  mb.Mm = eb.fQm / (seb.in.vp_ground.density_w * seb.params.Hf);

  // Condensation rate
  mb.Me = eb.fQe / (seb.in.vp_ground.density_w * seb.params.Ls);

  // Snow balance
  if (seb.in.snow_old.ht > 0.) {
    double swe_old = seb.in.snow_old.ht * seb.in.snow_old.density / seb.in.vp_ground.density_w;
    double swe_new = swe_old + (seb.in.met.Ps - mb.Mm + mb.Me)*seb.in.dt;
    std::cout<<"OLD_SWE: "<<swe_old<<std::endl;    
    // First do a pass to ensure we are not melting or sublimation ALL
    // of the available snow.  If so, adjust dt
    if (swe_new < 0.) {
      // Must adjust... we are melting or sublimating all of the available snow.
      ASSERT(seb.in.met.Ps - mb.Mm + mb.Me < 0.);
      mb.dt = -swe_old / (seb.in.met.Ps - mb.Mm + mb.Me);<commit_msg>remove print statement from seb<commit_after>/*
  Functions for calculating the snow / surface energy balance.
*/

#include <cmath>
#include <algorithm>
#include <boost/math/tools/roots.hpp>

#include "dbc.hh"

#include "seb_physics_funcs.hh"

namespace Amanzi {
namespace SurfaceBalance {
namespace SEBPhysics {

#define SWE_EPS 1.e-12

void UpdateIncomingRadiation(const SEB& seb, EnergyBalance& eb, bool debug) {
  // Calculate incoming short-wave radiation
  eb.fQswIn = (1 - seb.in.surf.albedo) * seb.in.met.QswIn;

  // Calculate incoming long-wave radiation
  const ThermoProperties& vp_air = seb.in.met.vp_air;
  double e_air = std::pow(10*vp_air.actual_vaporpressure, vp_air.temp / 2016.);
  e_air = 1.08 * (1 - std::exp(-e_air));
  eb.fQlwIn = e_air * seb.params.stephB * std::pow(vp_air.temp,4);

  // Calculate D_h, D_e, 
  eb.Dhe = std::pow(seb.params.VKc,2) * seb.in.met.Us
                       / std::pow(std::log(seb.params.Zr / seb.in.surf.Zo), 2);

  if (debug) {
    std::cout << "Incoming Radiation Energy Terms:" << "\n"
              << "  windspeed, Zo: " << seb.in.met.Us <<"  "<<seb.in.surf.Zo << "\n"
              << "  fQswIn   = " << eb.fQswIn << "\n"
              << "  fQlwIn   = " << eb.fQlwIn << std::endl;
  }

}

void UpdateEvapResistance(const SEB& seb, const ThermoProperties& vp_surf, EnergyBalance& eb, bool debug) {
   const ThermoProperties& vp_air = seb.in.met.vp_air;
   
   double Vaper_direction = vp_air.actual_vaporpressure - vp_surf.actual_vaporpressure;  //
   
// Equation for reduced vapor diffusivity See Sakagucki and Zeng 2009 eqaution (9) and Moldrup et al., 2004. 
   double Clab_Horn_b = 1;
   double actual_porosity = 0.9;  // Hard coded for moss Fix this to pass in form ATS ~AA
   double Surface_Vap_Diffusion = std::pow((1-(0.0556/actual_porosity)),(2+3*Clab_Horn_b));
   Surface_Vap_Diffusion = 0.000022 * (std::pow(actual_porosity,2)) * Surface_Vap_Diffusion;
// Sakagucki and Zeng 2009 eqaution (10)
   double cell_dimension = 0.01/2; // This is from cell center to the boundary.
   double VWC = seb.in.surf.saturation_liquid * actual_porosity;
   double L_Rsoil = std::exp(std::pow((1-(VWC/actual_porosity)),5));
   L_Rsoil = cell_dimension * (L_Rsoil -1) * (1/(2.718-1));
   double Rsoil = 0.0;
   
   if(Vaper_direction <= 0){
      Rsoil = L_Rsoil/Surface_Vap_Diffusion;
      }else{
      Rsoil = 0.0;
   }
   double Rair = 1/eb.Dhe;
   eb.Evap_Resistance = Rair;

   if((seb.in.snow_old.ht==0)&&(seb.in.surf.saturation_liquid < 1)) {
      eb.Evap_Resistance = Rair + Rsoil;
   }else{  
      eb.Evap_Resistance = Rair;
   }

   if (debug) {
   std::cout<<"Snow ht: "<<seb.out.snow_new.ht<<"  Saturation: "<<seb.in.surf.saturation_liquid<<std::endl;
   std::cout<<"Air_Temp: "<<vp_air.temp<<" Ground Temp: "<<vp_surf.temp<< std::endl;
   std::cout<<"Vapor Direction: "<<Vaper_direction<<" Air Vapor Pres: "<<vp_air.actual_vaporpressure<<" Surf Vapor Pres: "<<vp_surf.actual_vaporpressure << std::endl;
   std::cout<<"Surface_Vap_Diffusion: "<< Surface_Vap_Diffusion <<"  L_Rsoil:  "<<L_Rsoil<<"  Rsoil: "<<Rsoil<<"   VWC: "<<VWC<<std::endl;
   std::cout<<"Air Resistance: "<<Rair<<std::endl;
   std::cout<<"Evap Resistance: "<<eb.Evap_Resistance<<std::endl; 
   }
}


void UpdateEnergyBalance(const SEB& seb, const ThermoProperties& vp_surf, EnergyBalance& eb, bool debug) {
  const ThermoProperties& vp_air = seb.in.met.vp_air;

  // Calculate outgoing long-wave radiation
  eb.fQlwOut = -seb.in.surf.emissivity*seb.params.stephB*std::pow(vp_surf.temp,4);

  double Sqig;          // Stability function for use in e and h, precalculated for efficiency
  double air_temp = seb.in.met.vp_air.temp;
  double Ri  = seb.params.gravity * seb.params.Zr * (air_temp - vp_surf.temp)
      / (air_temp * std::pow(seb.in.met.Us,2));
  if (Ri >= 0.) {
    // stable condition or snow?
    Sqig = 1 / (1 + 10*Ri);
  } else {
    // Unstable condition
    Sqig = (1-10*Ri);
  }

  // Calculate sensible heat flux
  eb.fQh = seb.params.density_air * seb.params.Cp * eb.Dhe * Sqig * (vp_air.temp - vp_surf.temp);

  // Calculate latent heat flux
  double LatenHeatOf = 0.0; 
  if (seb.in.snow_old.ht > 0.){
    LatenHeatOf = seb.params.Ls;
   }else{
    LatenHeatOf = seb.params.Le;
  } 

  // With New Evaporation resistance term
  eb.fQe = vp_surf.porosity * seb.params.density_air * LatenHeatOf * (1/eb.Evap_Resistance) * Sqig * 0.622
      * (vp_air.actual_vaporpressure - vp_surf.actual_vaporpressure) / seb.params.Apa;
  if (debug) {
    std::cout<<"Porosity: "<<vp_surf.porosity<<"  LatentHeatOf: "<<LatenHeatOf<<"  Sqig: "<<Sqig<<std::endl;
    std::cout<<"Evap_Resistance: "<<eb.Evap_Resistance<<" Air vapor pres: "<<vp_air.actual_vaporpressure<<" Snow vapor pres: "<<vp_surf.actual_vaporpressure<<std::endl;
  }

  // Calculate heat conducted to ground, if snow
  if (seb.in.snow_old.ht > 0.) {
    double Ks = 2.9e-6 * std::pow(seb.in.snow_old.density,2);
    eb.fQc = Ks * (vp_surf.temp - seb.in.vp_ground.temp) / seb.in.snow_old.ht;
  }


  if (debug) {
    std::cout << "Energy Balance Terms (ht_snow = " << seb.in.snow_old.ht << "):" << "\n"
              << "  SnowSurfaceTemp  = " << vp_surf.temp << "\n"
              << "  fQlwOut  = " << eb.fQlwOut << "\n"
              << "  fQh      = " << eb.fQh << "\n"
              << "  fQe      = " << eb.fQe << "\n"
              << "  fQc      = " << eb.fQc << std::endl;
  }

}


void UpdateMassBalance(const SEB& seb, MassBalance& mb, EnergyBalance& eb, SnowProperties& snow_new, bool debug) {
  // this dt is the max timestep that may be taken to conserve snow mass
  mb.dt = seb.in.dt;

  // Melt rate given by available energy rate divided by heat of fusion.
  mb.Mm = eb.fQm / (seb.in.vp_ground.density_w * seb.params.Hf);

  // Condensation rate
  mb.Me = eb.fQe / (seb.in.vp_ground.density_w * seb.params.Ls);

  // Snow balance
  if (seb.in.snow_old.ht > 0.) {
    double swe_old = seb.in.snow_old.ht * seb.in.snow_old.density / seb.in.vp_ground.density_w;
    double swe_new = swe_old + (seb.in.met.Ps - mb.Mm + mb.Me)*seb.in.dt;
    // First do a pass to ensure we are not melting or sublimation ALL
    // of the available snow.  If so, adjust dt
    if (swe_new < 0.) {
      // Must adjust... we are melting or sublimating all of the available snow.
      ASSERT(seb.in.met.Ps - mb.Mm + mb.Me < 0.);
      mb.dt = -swe_old / (seb.in.met.Ps - mb.Mm + mb.Me);
{'fraction_non_alphanumeric': 0.15975039001560062, 'fraction_numerical': 0.015288611544461778, 'mean_word_length': 3.025431711145997, 'pattern_counts': {'":': 0, '<': 335, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2406155', 'n_tokens_mistral': 5452, 'n_tokens_neox': 5108, 'n_words': 1043}