#include "PotentialEvaporationCalc.h"

//PotentialEvaporationCalc.cpp functions compute potential evaporation, evapotranspiration, and sublimation

//References: 
//Allen, R.G. (1998). Chapter 2 - FAO Penman-Monteith equation. https://www.fao.org/3/x0490e/x0490e06.htm
//Chin, D. A. (2021). Water Resources Engineering, Fourth Edition. Hoboken, NJ: Pearson Education.
//Fassnacht, S.R. 2004 Estimating..Snowpack Simulation, Hydrologic Processes 18:3481-3492 https://onlinelibrary.wiley.com/doi/10.1002/hyp.5806
//Jensen, M. E., & Allen, R. G. (2015). Evaporation, Evapotranspiration, and Irrigation Water Requirements: Manual of Practice 70 (2nd ed.). Reston, VA: American Society of Civil Engineers.
//Light, P. (1941). Analysis of high rates of snow-melting. Eos, Transactions American Geophysical Union, 22(1), 195-205. doi:https://doi.org/10.1029/TR022i001p00195
//Lundberg, A., Calder, I., & Harding, R. (1998). Evaporation of intercepted snow: measurement and modelling. Journal of Hydrology, 206(3-4), 151-163. 
//McCutcheon, S. C., Martin, J. L., & Barnwell, T. O. J. (1993). Water Quality. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 11.11-11.73). New York: McGraw-Hill.
//Shuttleworth, J. W. (1993). Evaporation. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 4.1-4.5.3). New York: McGraw-Hill.
//Stull, R. B. (2000). Meteorology for Scientists and Engineers (2nd ed.). Pacific Grove, CA: Brooks/Cole.

//Note: The following are not directly used but may prove useful in code refactoring
//Sanow, J. E., Fassnacht, S. R., Kamin, D. J., Sexstone, G. A., Bauerle, W. L., & Oprea, I. (2018). Geometric Versus Anemometric Surface Roughness for a Shallow Accumulating Snowpack. Geosciences, 8(12), 463. 
//Note: Sanow et al. (2018) demonstrate methods to adjust snow roughness length based on snow depth; median value z0_snow = 0.024 used by UCAR Community Land Model
//Sexstone, G. A., Clow, D. W., Fassnacht, S. R., Liston, G. E., Hiemstra, C. A., Knowles, J. F., & Penn, C. A. (2018). Snow Sublimation in Mountain Environments and Its Sensitivity to Forest Disturbance and Climate Warming. Water Resources Research, 54(2), 1191-1211. doi:https://doi.org/10.1002/2017WR021172
//Note: Sexstone et al. (2018) use SnowModel and demonstrate relative rates of sublimation, melt, and wind on snow ablation for ground and canopy. Canopy loss not always > Ground loss.

//Note: Theory in email to T. Endreny on 2/4/17 from Rick Allen, an expert in Penman Monteith evaporation methods, author of FAO and ASCE manuals. See Jensen et al. (2015).
//Note:		Allen points to Chapter 3, 5, 7, and 11 in Jensen et al. (2015)
//Note: ... Allen states: Trees exercise higher surface resistance, rs, than shorter plants like shrubs due to:
//Note: ...		a) larger roughness, b) larger leaf area; and c) therefore the need to increase rs to conserve water under a wide range of water availabilities.
//Note: ... Allen states: For trees and shrubs, I would estimate rs to be higher than estimated by the rs ~ 200/LAI relation.  
//Note: ...		The rs ~ 200/LAI relation is mostly for annual agricultural crops that have been bred to have fully open stomates as frequently as possible.  
//Note: ...		Most trees and shrubs have evolved to be more in 'survival mode' than production mode, and therefore tend to reduce their stomatal openings more of the time.  
//Note: ...		In addition, the larger roughness of trees and shrubs makes them more vulnerable to loss of water and therefore, again, they tend to exercise more stomatal control.
//Note: Theory of Farquhar Photosynthesis Model is alternative for transpiration, used in LiBry (Phlipp Porada) and several Earth System Models
//Note: ...     https://biocycle.atmos.colostate.edu/shiny/photosynthesis/#:~:text=Photosynthetic%20Responses-,MODEL%20DESCRIPTION,(1993)%20and%20Bonan%202008.
//Method_PenmanMonteith function will use Penman Monteith to compute potential evaporation and transpiration
void PotentialEvaporationCalc::Method_PenmanMonteith(Inputs *input, CompactRagged* beC, WeatherProcessor* WeatherPro, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
    //If timeStep equals size of evaporation vector, PotentialEvapotranspiration_TreeCover_m.size(), then proceed
    //Note: Conditional prohibits computing and appending evaporation multiple times per time step; a better gateway conditional may be useful
    if (timeStep == input->PotentialEvapotranspiration_TreeCover_m.size()) {

        //Note: StatisticalHydro model when MultipleStation is true assigns input-> to use WeightedInputMap[MapPixel_ID][Station_Variable_name] 

        //Temperature_Air_C (C) is Tair_C (C) 
        double Temperature_Air_C = input->Tair_C[timeStep];
        //Temperature_DewPoint_C (C) is Tdew_C (C) 
        double Temperature_DewPoint_C = input->Tdew_C[timeStep];
        //Pressure_Atmosphere_kPa (kPa) is defined from input AtmPres_kPa
        double Pressure_Atmosphere_kPa = input->AtmPres_kPa[timeStep];
        //WindSpeed_Station_m_p_s (m/s) defined from input WindSpd_mps (m/s
        double WindSpeed_Station_m_p_s = input->WindSpd_mps[timeStep];
        //Height_Tree_m (m) is average height of feature for roughness in turbulence exchange calculations
        double Height_Tree_m = input->InputXml["Height_Avg_TreeCover_m"];
        //Height_SVeg_m (m) is average height of feature for roughness in turbulence exchange calculations
        double Height_SVeg_m = input->InputXml["Height_Avg_SVegCover_m"];
        //Height_WindMeasurement_m (m) is reference height for wind speed measurements
        double Height_WeatherStationWindSensor_m = input->SimulationNumericalParams["Height_Sensor_Wind_m"];
        //LAI_Tree_m2_p_m2 (m2/m2) defined from input LAI_Tree_m2_p_m2
        double LAI_Tree_m2_p_m2 = input->LAI_Tree_m2_p_m2[timeStep];

        //Radiation_Net_Wpm2 (W/m2) is energy from downwelling shortwave and longwave radiation (W/m^2)
        double Radiation_Net_Wpm2 = input->RadiationNet_Wpm2[timeStep];
        //Radiation_Shortwave_Direct_Wpm2 (W/m2) is energy from direct shortwave radiation (W/m^2)
        double Radiation_Shortwave_Direct_Wpm2 = input->Radiation_Shortwave_Direct_Wpm2[timeStep];
        //Radiation_Shortwave_Diffuse_Wpm2 (W/m2) is energy from diffuse shortwave radiation (W/m^2)
        double Radiation_Shortwave_Diffuse_Wpm2 = input->Radiation_Shortwave_Diffuse_Wpm2[timeStep];
        //Radiation_Longwave_Upwelling_Wpm2 (W/m2) is energy in upwelling longwave radiation (W/m^2)
        double Radiation_Longwave_Upwelling_Wpm2 = input->Radiation_Longwave_Upwelling_Wpm2[timeStep];


        //SoilAreaIndex_m2_p_m2 for trees, shrubs, and ground is estimated at 1
        double SoilAreaIndex_m2_p_m2 = 1;

        //Initialize local variable names for potential evaporation and sublimation (m/s) as computed, before conversion to m per time step
        double PotentialEvapotranspiration_TreeCover_m_p_s = 0;
        double PotentialEvaporation_WaterOnGround_m_p_s = 0;
        double PotentialEvaporation_WaterOnTree_m_p_s = 0;
        double PotentialSublimation_SnowOnGround_m_p_s = 0;
        double PotentialSublimation_SnowOnTree_m_p_s = 0;
        //Initialize local variable names for potential evaporation and sublimation (m/s) as m per time step
        double PotentialEvapotranspiration_TreeCover_m = 0;
        double PotentialEvaporation_WaterOnGround_m = 0;
        double PotentialEvaporation_WaterOnTree_m = 0;
        double PotentialSublimation_SnowOnGround_m = 0;
        double PotentialSublimation_SnowOnTree_m = 0;

        //Height_TreeWindSensor_m created to avoid negative logarithms below, ensuring tree height is below wind sensor height
        //Note: Error Handling: If users enter wind sensor height = 2 m, tree height = 5 m, then aerodynamic resistance equations fail
        //Note: ... Aerodynamic resistance equations take log of difference of sensor height and zero plane displacement height of object
        //Note: ... tree height * zero plane displacement fraction = 5 * 0.67 = 3.35; ln(2-3.35) = No Solution
        double Height_TreeWindSensor_m = Height_WeatherStationWindSensor_m + Height_Tree_m;
        double Height_TreeZeroPlaneDisplacement_m = Height_Tree_m * ZeroPlaneDisplacementHeight_frac;
        //Height_Ground_m is set to 0.1 m to represent surface irregularity
        //Note: If Height_Ground_m converted to roughness length with Zom scale factor = 0.123 then it has classification between open and smooth 
        double Height_Ground_m = 0.1;
        double Height_GroundWindSensor_m = Height_WeatherStationWindSensor_m + Height_Ground_m;
        double Height_GroundZeroPlaneDisplacement_m = Height_Ground_m * ZeroPlaneDisplacementHeight_frac;

        double WindSpeed_TreeCanopy_m_p_s;
        double DensityWater_kgpm3;
        double DensityAir_kg_p_m3;
        double VaporPressure_Actual_kPa;
        double VaporPressure_Saturated_kPa;
        double VaporPressure_Deficit_kPa;
        double VaporPressure_Deficit_Pa;
        double VaporPressure_Deficit_OverWater_Pa;

        double LeafAreaIndex_Active_Tree_m2_p_m2;
        double ResistanceSurface_Tree_and_Soil_s_p_m;
        double ResistanceSurface_Water_s_p_m;
        double ResistanceAerodynamic_TreeCanopy_s_p_m;
        double ResistanceAerodynamic_SurfaceWater_s_p_m;
        double ResistanceAerodynamic_SnowTreeCanopy_s_p_m;
        double ResistanceAerodynamic_SnowGround_s_p_m;
        double LAI_plus_BAI_plus_GAI_SVeg_to_Tree_ratio;
        double LAI_plus_BAI_plus_GAI_Soil_to_Tree_ratio;

        double LatentHeat_Vaporization_JpkgK;
        double LatentHeatOfVaporizaton_MJpkgK;
        double PsychrometricConstant_Pa_p_C;
        double PsychrometricConstant_kPa_p_C;
        double VaporPressureGradient_Pa;
        double VaporPressureGradient_kPa;
        double SpecificHeat_MoistAir_MJ_p_kg_p_C;
        double GroundHeatFlux_W_p_m2;
        double ResistanceSurface_Vegetation_s_p_m;
        double BulkStomatalResistance_Adjustment_frac;
        double HeatStorage_Water_J_p_m2_p_s;
        double ResistanceAerodynamic_Snow_to_Rain_ratio;


        //VaporPressure_Saturated_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
        VaporPressure_Saturated_kPa = 0.6108 * exp((17.27 * Temperature_Air_C) / (237.3 + Temperature_Air_C));
        //VaporPressure_Actual_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
        VaporPressure_Actual_kPa = 0.6108 * exp((17.27 * Temperature_DewPoint_C) / (237.3 + Temperature_DewPoint_C));
        //VaporPressure_Deficit_kPa is vapor pressure deficit (kPa), using previously computed saturated vapor pressure SatVPrsKpa (kPa)
        VaporPressure_Deficit_kPa = (VaporPressure_Saturated_kPa - VaporPressure_Actual_kPa);
        //VaporPressure_Deficit_Pa (Pa) is vapor pressure deficit in fundamental units used in Penman-Monteith Eq 13.1 of Chin (2021) 
        VaporPressure_Deficit_Pa = VaporPressure_Deficit_kPa * Ratio_Pa_to_kPa;
        //VaporPressure_Deficit_OverWater_Pa (Pa) is vapor pressure deficit over water, initially equal to at station value
        VaporPressure_Deficit_OverWater_Pa = VaporPressure_Deficit_Pa;

        //Note: 202212 This option currently remains off until further evidence that it is needed.
        //Note: Adjustment to approximate theory of Jensen and Allen (2015) Chapter 6, Evaporation from Water Surfaces:
        //Note: ... "Evaporation can be low during summer when shortwave radiative energy is absorbed by the cold water bodies and
        //Note: ... when vapor pressure gradients above the water are small due to high humidity of air caused by regional ET.
        //Note: ... Evaporation is high during winter as dry air from nontranspiring frozen regions is coupled with saturation
        //Note: ... vapor pressure at the water surface that may be much higher than the air."
        //If TempC > freezing then reduce vapor pressure deficit due to high humidity of air by evaporation; otherwise no adjustment
        /*
        if (weatherData[i].TempC > 0) {
            //VaporPressure_Deficit_Reduction_frac set to 0.75, but should be adjusted to fit evidence of altered deficit
            double VaporPressure_Deficit_Reduction_frac = 0.75;
            //VaporPressure_Deficit_Adjustment_frac (frac) formula derived here; should be adjusted to fit observations; 50 degrees used as asymptote
            double VaporPressure_Deficit_Adjustment_frac = 1 - (Math.Pow(weatherData[i].TempC / 50, 2) * VaporPressure_Deficit_Reduction_frac);
            //If VaporPressure_Deficit_Adjustment_frac goes below VaporPressure_Deficit_Reduction_frac, then set to that limit
            if (VaporPressure_Deficit_Adjustment_frac < VaporPressure_Deficit_Reduction_frac) {
                VaporPressure_Deficit_Adjustment_frac = VaporPressure_Deficit_Reduction_frac; }
            //VaporPressure_Deficit_OverWater_Pa (Pa) reduced to VaporPressure_Deficit_Adjustment_frac its magnitude
            VaporPressure_Deficit_OverWater_Pa = VaporPressure_Deficit_OverWater_Pa * VaporPressure_Deficit_Adjustment_frac;
        }
        */

        //If WindSpeed_Station_m_p_s < = zero, then set to 0.1 m/s to avoid division by zero
        if (WindSpeed_Station_m_p_s <= 0) { WindSpeed_Station_m_p_s = 0.1; }
        //WindSpeed_TreeCanopy_m_p_s defined as wind speed at tree top (m/s) from Eq 4.14b of Stull (2000)
        WindSpeed_TreeCanopy_m_p_s = WindSpeed_Station_m_p_s * (log(Height_TreeWindSensor_m / RoughnessLength_Airport_m) /
            log(Height_WeatherStationWindSensor_m / RoughnessLength_Airport_m));

        //DensityAir_kg_p_m3 (kg/m3) is density of air from Eq 4.2.4 from Shuttleworth (1993), w correction in conversion to C from K
        //Note: Eq 4.24 used incorrect denominator of 275 rather than 273.15 to convert from C to K; see Chin (2021) Eq 13.51
        //Note: Chin Eq 13.51 which unfortunately uses 3.45 in place of the correct 3.486.
        //Note: This tested well against values of air density from the EngineeringToolbox.com for air temperature from 0 to 50 C at standard atmospheric pressure
        DensityAir_kg_p_m3 = 3.486 * (Pressure_Atmosphere_kPa / (273.15 + Temperature_Air_C));

        //Calculate HeatMetrics with call to HeatMetrics_Calc functions
        //HeatMetrics_Calc HeatMetricsCalc(input) creates pointer to access HeatMetrics_Calc functions with sending input
        HeatMetrics_Calc HeatMetricsCalc(input);

        //DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"))
        DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

        //LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(Tair_K)
        LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));
        //LatentHeatOfVaporizaton_MJpkgK = LatentHeat_Vaporization_JpkgK * Ratio_MJ_to_J
        LatentHeatOfVaporizaton_MJpkgK = LatentHeat_Vaporization_JpkgK * Ratio_MJ_to_J;

        //PsychrometricConstant_kPa_p_C (kPa/C) from Eq 13.37 Chin (2021) is psychrometric constant
        //Note: Eq 13.37 should use specific heat with units of MJ/kg/C, but incorrectly states units are kJ/kg/C
        SpecificHeat_MoistAir_MJ_p_kg_p_C = SpecificHeat_MoistAir_J_p_kg_p_C * Ratio_MJ_to_J;
        PsychrometricConstant_kPa_p_C = (SpecificHeat_MoistAir_MJ_p_kg_p_C * Pressure_Atmosphere_kPa) /
            (Ratio_MolecularWeightVapor_to_MolecularWeightDryAir * LatentHeatOfVaporizaton_MJpkgK);
        //PsychrometricConstant_Pa_p_C (Pa/C) is converted from kPa to Pa
        PsychrometricConstant_Pa_p_C = PsychrometricConstant_kPa_p_C * Ratio_Pa_to_kPa;

        //VaporPressureGradient_kPa (kPa) from Eq 13.43 Chin (2021) is gradient of saturation pressure vs temperature curve
        VaporPressureGradient_kPa = 4098 * (0.6108 * exp((17.27 * Temperature_Air_C) / (Temperature_Air_C + 237.3))) /
            pow((Temperature_Air_C + 237.3), 2);
        //VaporPressureGradient_Pa (Pa) is converted from kPa to Pa
        VaporPressureGradient_Pa = VaporPressureGradient_kPa * Ratio_Pa_to_kPa;


        //GroundHeatFlux_W_p_m2 (W/m2) defined for two vegetation heights after Eq 13.32 by Chin (2021); also Jensen et. al., (2015) Chapter 5
        //Note: Jensen and Allen (2015) "Standardized Reference Evapotranspiration for Short Reference ETos: Reference ET calculated
        //Note: ... for a short crop having height of 0.12 m (similar to grass), albedo of 0.23, surface resistance of 70 sm−1 for
        //Note: ... 24-h calculation time steps, and 50 sm−1 for hourly or shorter periods during daytime and 200 sm−1 during"
        //Note: Jensen and Allen (2015) "Standardized Reference Evapotranspiration for Tall Reference ETrs: Reference ET calculated
        //Note: ... for a tall crop having height of 0.50m (similar to alfalfa), albedo of 0.23, surface resistance of 45 sm−1 for
        //Note: ... 24-h calculation time steps, and 30 sm−1 for hourly or shorter periods during daytime and 200 sm−1 during nighttime
        //Note: ... with daytime G=Rn =0.04 and nighttime G=Rn =0.2."
        //Note: Short crop coefficients 0.1 during day, 0.5 during night. Tall crop coefficients 0.04 during day, 0.2 during night. 
        //Note: Jensen et al. (2015) generally suggest using reference crop coefficients and then adjusting to actual conditions
        //If Radiation_Net_Wpm2 (W/m2) positive, presume during day then GroundHeatFlux_W_p_m2 takes one form
        if (Radiation_Net_Wpm2 > 0) {
            //GroundHeatFlux_W_p_m2 (W/m2) defined for tall crop after Jensen et. al., (2015) Chapter 5
            GroundHeatFlux_W_p_m2 = 0.04 * Radiation_Net_Wpm2;
        }
        //Else If Radiation_Net_Wpm2 (W/m2) negative, presume during night then GroundHeatFlux_W_p_m2 takes another form
        else {
            //GroundHeatFlux_W_p_m2 (W/m2) defined for tall crop after Jensen et. al., (2015) Chapter 5
            GroundHeatFlux_W_p_m2 = 0.2 * Radiation_Net_Wpm2;
        }

        //LeafAreaIndex_Active_Tree_m2_p_m2 (m2/m2) Eq 13.11 from Chin (2021)
        //Note: Chin (2021) cites Allen (1998), explaining crops typically have 50% of LAI active; trees may have less
        LeafAreaIndex_Active_Tree_m2_p_m2 = LAI_Tree_m2_p_m2 * 0.5;

        //BulkStomatalResistance_s_p_m (s/m) = 200 based on Shuttleworth (1990) Eq 4.2.22 and interpretation of Allen (1998), Chin (2021) around Eq 13.9
        //Note: ResistanceSurface_Tree_and_Soil_s_p_m derived from BulkStomatalResistance_s_p_m by division with LAI
        //Note: Minimum values of ResistanceSurface_Tree_and_Soil_s_p_m ~ 100 s/m in Table 2 of Liang et al. (1994) VIC model and ...
        //Note: ... Box 5 of Allen et al. (1998) Chapter 2 - FAO Penman-Monteith equation https://www.fao.org/3/x0490e/x0490e06.htm, ...
        //Note: ... and minimum values of ResistanceSurface_Tree_and_Soil_s_p_m ~ 40 s/m are found in Chpt 11 of Jensen and Allen (2015)
        ResistanceSurface_Vegetation_s_p_m = 200;

        //BulkStomatalResistance_Adjustment_frac based on suggestion of Jensen and Allen (2015) Chapter 4 Energy Balance:
        //Note: ... "In addition, a recommendation by Allen et al. (2006a) to use the same 50 sm−1 surface resistance for hourly
        //Note: ... or shorter periods during daytime and 200 sm−1 during nighttime with the FAO-56 Penman-Monteith method has
        //Note: ... made the FAO ETo and ETos references equivalent for hourly time steps and for 24-h periods."
        //Note: Allen (personal communication) suggests larger stomatal resistance for vegetation in a landscape,
        //Note: ... i.e., not irrigated reference crop, stating it evolved to survive by restricting water loss. 
        //Note: Equation derived here and coefficient 0.9 was adjusted to fit expected trends in evapotranspiration 
        BulkStomatalResistance_Adjustment_frac = pow((1 - VaporPressure_Deficit_kPa / VaporPressure_Saturated_kPa), 0.9);
        ResistanceSurface_Vegetation_s_p_m = ResistanceSurface_Vegetation_s_p_m / BulkStomatalResistance_Adjustment_frac;

        //If LAI_Tree_m2_p_m2 < 1 then prohibit ResistanceSurface_Tree_and_Soil_s_p_m from going toward infinity with division by LAI < 1
        if (LAI_Tree_m2_p_m2 < 1) {
            ResistanceSurface_Tree_and_Soil_s_p_m = ResistanceSurface_Vegetation_s_p_m / 1.0;
        }
        //Else If LAI_Tree_m2_p_m2 >= 1 then ResistanceSurface_Tree_and_Soil_s_p_m divided by LeafAreaIndex_Active_Tree_m2_p_m2 
        else {
            //ResistanceSurface_Tree_and_Soil_s_p_m (s/m) defined with Eq 13.9 in Chin (2021) and Eq 4.2.22 in Shuttleworth (1993) 
            ResistanceSurface_Tree_and_Soil_s_p_m = ResistanceSurface_Vegetation_s_p_m / LeafAreaIndex_Active_Tree_m2_p_m2;
        }
        //ResistanceSurface_Water_s_p_m defined as 0, which is the case for water surface, defined by Chin (2021) after Eq 13.12
        ResistanceSurface_Water_s_p_m = 0;

        //LAI_plus_BAI_plus_GAI_SVeg_to_Tree_ratio is ratio of (LAI + BAI + GAI) for tree to (LAI + BAI + GAI) for short veg
        //Note: LAI = leaf area index, BAI = bark area index, GAI = ground area index, where GAI = 1 for bare and snow covered ground
        LAI_plus_BAI_plus_GAI_SVeg_to_Tree_ratio = (input->LAI_BAI_Tree_m2_p_m2[timeStep] + SoilAreaIndex_m2_p_m2) /
            (input->LAI_BAI_SVeg_m2_p_m2[timeStep] + SoilAreaIndex_m2_p_m2);
        //LAI_plus_BAI_plus_GAI_Soil_to_Tree_ratio is ratio of (LAI + BAI + GAI) for ground to (LAI + BAI + GAI) for tree
        //Note: (LAI + BAI + GAI) for ground = (0 + 0 + GAI) for ground, despite overhanging vegetation
        LAI_plus_BAI_plus_GAI_Soil_to_Tree_ratio = (0 + 0 + SoilAreaIndex_m2_p_m2) /
            (input->LAI_BAI_Tree_m2_p_m2[timeStep] + SoilAreaIndex_m2_p_m2);

        //ResistanceAerodynamic_TreeCanopy_s_p_m (s/m) from Eq 13.2 of Chin (2021) or Eq 4.2.25 of Shuttleworth (1993) 
        //Note: Eq 4.2.25 should use the station windspeed, not a windspeed estimated at the height of any lower object
        //Note: Height_TreeWindSensor_m is used in place of Height_WeatherStationWindSensor_m to ensure real values from log
        ResistanceAerodynamic_TreeCanopy_s_p_m = (log((Height_TreeWindSensor_m - Height_TreeZeroPlaneDisplacement_m) /
            (Zom_RoughnessLengthMomentumTransferCoefficient * Height_Tree_m)) * log((Height_TreeWindSensor_m - Height_TreeZeroPlaneDisplacement_m) /
                (Zom_RoughnessLengthHeatVaporTransferCoefficient * Height_Tree_m))) / (pow(vonKarman_Constant, 2) * WindSpeed_Station_m_p_s);

        //ResistanceAerodynamic_SurfaceWater_s_p_m is Eq 13.8 Chin (2021) or Eq 4.2.29 from Shuttleworth (1993) for open water
        //Note: RoughnessLength_Water_m = 1.37 mm from Chin (2021)
        //Note: Height_WeatherStationWindSensor_m can be any height with accompanying WindSpeed_Station_m_p_s
        ResistanceAerodynamic_SurfaceWater_s_p_m = (4.72 * pow((log(Height_WeatherStationWindSensor_m / RoughnessLength_Water_m)), 2)) /
            (1 + 0.536 * WindSpeed_Station_m_p_s);

        //ResistanceAerodynamic_Snow_to_Rain_ratio set to 10 according to Lundberg et al. (1998) Eq 6 and 7
        //Note: ResistanceAerodynamic_Snow_to_Rain_ratio is ratio of aerodynamic resistance for snow to rain 
        ResistanceAerodynamic_Snow_to_Rain_ratio = 10;

        //ResistanceAerodynamic_SnowTreeCanopy_s_p_m contains resistance terms in Eq 4 of Fassnacht (2004) & Eq 3 of Light (1941)
        ResistanceAerodynamic_SnowTreeCanopy_s_p_m = (log(Height_TreeWindSensor_m / RoughnessLength_Snow_m) *
            log(Height_TreeWindSensor_m / RoughnessLength_Snow_m)) / (pow(vonKarman_Constant, 2) * WindSpeed_TreeCanopy_m_p_s);

        //ResistanceAerodynamic_SnowGround_s_p_m contains resistance terms in Eq 4 of Fassnacht (2004) & Eq 3 of Light (1941)
        //Note: ResistanceAerodynamic_SnowGround_s_p_m uses wind speed measured at station height
        ResistanceAerodynamic_SnowGround_s_p_m = (log(Height_WeatherStationWindSensor_m / RoughnessLength_Snow_m) *
            log(Height_WeatherStationWindSensor_m / RoughnessLength_Snow_m)) / (pow(vonKarman_Constant, 2) * WindSpeed_Station_m_p_s);

        //PotentialEvapotranspiration_TreeCover_m_p_s (m/s), evapotranspiration from canopy and ground, from Eq 13.1 of Chin (2021), Penman Monteith Evaporation method
        //Note: Equation 13.1 of Chin (2021) requires all terms are in fundamental SI units (J, Pa, kg, s, m, C, etc.)
        //Note: Issues with Eq 4.2.27 of Shuttleworth (1993) include not dividing by density of water, using mixed units MJ/day, kPa/C, etc.
        PotentialEvapotranspiration_TreeCover_m_p_s = 1 / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) *
            ((VaporPressureGradient_Pa * (Radiation_Net_Wpm2 - GroundHeatFlux_W_p_m2)) +
                (DensityAir_kg_p_m3 * SpecificHeat_MoistAir_J_p_kg_p_C * VaporPressure_Deficit_Pa / ResistanceAerodynamic_TreeCanopy_s_p_m)) /
            (VaporPressureGradient_Pa + PsychrometricConstant_Pa_p_C *
                (1 + ResistanceSurface_Tree_and_Soil_s_p_m / ResistanceAerodynamic_TreeCanopy_s_p_m));
        //PotentialEvapotranspiration_TreeCover_m (m) converted from (m/s) with SimulationTimeStep_Duration_sec[timeStep]
        PotentialEvapotranspiration_TreeCover_m = PotentialEvapotranspiration_TreeCover_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
        //If PotentialEvapotranspiration_TreeCover_m negative then to zero
        if (PotentialEvapotranspiration_TreeCover_m < 0) { PotentialEvapotranspiration_TreeCover_m = 0.0; }
        //PotentialEvapotranspiration_TreeCover_m (m) stored in vector for each time step
        input->PotentialEvapotranspiration_TreeCover_m.push_back(PotentialEvapotranspiration_TreeCover_m);

        //Note: Adjustment to approximate theory of Jensen and Allen (2015) Chapter 6, Evaporation from Water Surfaces:
        //Note: ... "An important distinction between Rn for a water body and Rn for vegetation or soil is that with soil and
        //Note: ... vegetation, essentially all of the Rn quantity is captured at the “opaque” surface and is immediately available
        //Note: ... solar radiation, Rs, for conversion to λE or H or conduction into the surface as G. 
        //Note: ... With water, however, much of the penetrates to some depth in the water body, depending on the turbidity
        //Note: ... of the water, where it is converted to Qt, heat storage." This approach is modified for smaller waters below.
        //HeatStorage_Water_J_p_m2_p_s (J/m2/s) from Eq 6-14a and 6-14b of Jensen and Allen (2015) modified from large lakes
        //Note: Large lakes used HeatStorage_Coefficient_a = 0.5 and HeatStorage_Coefficient_b = 0.8; reduced for smaller waters
        //Note: Large lakes increased loss from long wave radiation at Julian Day = 180; not simulated for smaller waters
        //Note: Modification includes Heat Storage transferring energy into Ground Heat Flux 
        double HeatStorage_Coefficient_a = 0.25;
        //HeatStorage_Coefficient_b reduced from large lake value of 0.8
        double HeatStorage_Coefficient_b = 0.05;
        //HeatStorage_Water_J_p_m2_p_s (J/m2/s) heat storage of water from Eq 6-14a and 6-14b of Jensen and Allen (2015) 
        HeatStorage_Water_J_p_m2_p_s = HeatStorage_Coefficient_a * (Radiation_Shortwave_Direct_Wpm2 + Radiation_Shortwave_Diffuse_Wpm2) - 
            HeatStorage_Coefficient_b * Radiation_Longwave_Upwelling_Wpm2;

        //PotentialEvaporation_WaterOnGround_m_p_s (m/s), evaporation from water surface on ground, from Eq 13.1 of Chin (2021), Penman Monteith Evaporation method
        //Note: Equation 13.1 of Chin (2021) requires all terms are in fundamental SI units (J, Pa, kg, s, m, C, etc.)
        //Note: Issues with Eq 4.2.27 of Shuttleworth (1993) include not dividing by density of water, using mixed units MJ/day, kPa/C, etc.
        //Note: ResistanceSurface_Water_s_p_m set to zero. 
        PotentialEvaporation_WaterOnGround_m_p_s = 1 / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) *
            ((VaporPressureGradient_Pa * (Radiation_Net_Wpm2 - HeatStorage_Water_J_p_m2_p_s)) +
                (DensityAir_kg_p_m3 * SpecificHeat_MoistAir_J_p_kg_p_C * VaporPressure_Deficit_OverWater_Pa / ResistanceAerodynamic_SurfaceWater_s_p_m)) /
            (VaporPressureGradient_Pa + PsychrometricConstant_Pa_p_C *
                (1 + ResistanceSurface_Water_s_p_m / ResistanceAerodynamic_SurfaceWater_s_p_m));
        //PotentialEvaporation_WaterOnGround_m (m) converted from (m/s) with SimulationTimeStep_Duration_sec[timeStep]
        PotentialEvaporation_WaterOnGround_m = PotentialEvaporation_WaterOnGround_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
        //If PotentialEvaporation_WaterOnGround_m negative then to zero
        if (PotentialEvaporation_WaterOnGround_m < 0) { PotentialEvaporation_WaterOnGround_m = 0.0; }
        //PotentialEvaporation_WaterOnGround_m (m) stored in vector for each time step
        input->PotentialEvaporation_WaterOnGround_m.push_back(PotentialEvaporation_WaterOnGround_m);

        //PotentialEvaporation_WaterOnTree_m (m), evaporation from liquid water in canopy, set equal to PotentialEvaporation_WaterOnGround_m (m)
        PotentialEvaporation_WaterOnTree_m = PotentialEvaporation_WaterOnGround_m;
        //PotentialEvaporation_WaterOnTree_m (m) stored in vector for each time step
        input->PotentialEvaporation_WaterOnTree_m.push_back(PotentialEvaporation_WaterOnTree_m);

        //PotentialSublimation_SnowOnGround_m_p_s (m/s), sublimation from canopy, from Eq 2 of Lundberg et al. (1998), presented as Eq 13.1 of Chin (2021).
        //Note: Eq 2 of Lundberg et al. (1998) has no rs/ra term in the denominator, which Eq 13.1 maintains, given surface resistance rs = 0
        //Note: ResistanceAerodynamic_Snow_to_Rain_ratio term multiplied w to represent Lundberg et al. (1998) 10x factor from Eq 6 to 7
        //Note: Theory: COMET UCAR Snowmelt Processes International Edition states: Snow on vegetation is exposed to ...
        //Note: ... more wind and sun than snow on the ground and has a higher surface area to mass ratio, thereby ...
        //Note: ... making it more prone to sublimation and/or melting.
        //PotentialSublimation_SnowOnGround_m_p_s (m/s), sublimation from ground, from Eq 2 of Lundberg et al. (1998), presented as Eq 13.1 of Chin (2021).
        //Note: Eq 2 of Lundberg et al. (1998) has no rs/ra term in the denominator, which Eq 13.1 maintains, given surface resistance rs = 0
        //Note: ResistanceAerodynamic_Snow_to_Rain_ratio term multiplied w to represent Lundberg et al. (1998) 10x factor from Eq 6 to 7
        PotentialSublimation_SnowOnGround_m_p_s = 1 / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) *
            ((VaporPressureGradient_Pa * (Radiation_Net_Wpm2 - HeatStorage_Water_J_p_m2_p_s)) +
                (DensityAir_kg_p_m3 * SpecificHeat_MoistAir_J_p_kg_p_C * VaporPressure_Deficit_Pa /
                    (ResistanceAerodynamic_SnowGround_s_p_m * ResistanceAerodynamic_Snow_to_Rain_ratio))) /
            (VaporPressureGradient_Pa + PsychrometricConstant_Pa_p_C *
                (1 + ResistanceSurface_Water_s_p_m / ResistanceAerodynamic_SnowGround_s_p_m));
        //PotentialSublimation_SnowOnGround_m (m) converted from (m/s) with SimulationTimeStep_Duration_sec[timeStep]
        PotentialSublimation_SnowOnGround_m = PotentialSublimation_SnowOnGround_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
        //If PotentialSublimation_SnowOnGround_m negative then to zero
        if (PotentialSublimation_SnowOnGround_m < 0) { PotentialSublimation_SnowOnGround_m = 0.0; }
        //PotentialSublimation_SnowOnGround_m (m) stored in vector for each time step
        input->PotentialSublimation_SnowOnGround_m.push_back(PotentialSublimation_SnowOnGround_m);

        //PotentialSublimation_SnowOnTree_m_p_s (m/s), sublimation from canopy, from Eq 2 of Lundberg et al. (1998), presented as Eq 13.1 of Chin (2021).
        PotentialSublimation_SnowOnTree_m_p_s = 1 / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) *
            ((VaporPressureGradient_Pa * (Radiation_Net_Wpm2 - HeatStorage_Water_J_p_m2_p_s)) +
                (DensityAir_kg_p_m3 * SpecificHeat_MoistAir_J_p_kg_p_C * VaporPressure_Deficit_Pa /
                    (ResistanceAerodynamic_SnowTreeCanopy_s_p_m * ResistanceAerodynamic_Snow_to_Rain_ratio))) /
            (VaporPressureGradient_Pa + PsychrometricConstant_Pa_p_C *
                (1 + ResistanceSurface_Water_s_p_m / ResistanceAerodynamic_SnowTreeCanopy_s_p_m));
        //PotentialSublimation_SnowOnTree_m (m) converted from (m/s) with SimulationTimeStep_Duration_sec[timeStep]
        PotentialSublimation_SnowOnTree_m = PotentialSublimation_SnowOnTree_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
        //If PotentialSublimation_SnowOnTree_m negative then to zero
        if (PotentialSublimation_SnowOnTree_m < 0) { PotentialSublimation_SnowOnTree_m = 0.0; }
        //PotentialSublimation_SnowOnTree_m (m) stored in vector for each time step
        input->PotentialSublimation_SnowOnTree_m.push_back(PotentialSublimation_SnowOnTree_m);
    }
}
