using Newtonsoft.Json.Linq;
using NHibernate.Criterion;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Reflection;
using System.Runtime.InteropServices;
using System.Runtime.Remoting.Messaging;
using System.Text;
using System.Text.RegularExpressions;
using System.Web.UI.WebControls;
using System.Web.UI.WebControls.WebParts;


namespace WeatherPrep
{
    /// <summary>
    /// Several functions to read raw surface weather data and compute weather variables
    /// </summary>
    public enum VEGTYPE { TREE, SHRUB, GRASS, CORN };

    public class SurfaceWeather
    {
        #region class_variables
        const double ConversionFactor_273pt15K_to_0C = 273.15;    ///<value>Conversino factor for Celsius to Kelvin</value>
        const double XLIM = 0.000000001;
        const double ABS_PRS = 1013.25;                     ///<value>mBar</value>
        const double ABS_TEMP = 288.15;                     ///<value>Kelvin</value>
        const int SOL_CONST = 1367;
        const double ES_CONS = 0.6108;                      ///<value>saturated vapor pressure constant (kPa)</value>
        const double ATRN = 0.957;
        const double K1 = 0.1;
        const double BA = 0.84;
        const double PA_CONS = 3.486;                       ///<value>Density for moist air sensible heat flux constant</value>
        const double CHI_CONS = 2.501;                      ///<value>latent heaf for evaoration of water, adjusted by temperature (used in meteo and evap routines)</value>
        const double GAMMA_C = 0.0016286;                   ///<value>psychrometric constant, adjusted by vapor pressure (used in meteo and evap routines)</value>
        const double Stefan_Boltzmann_sigma_Wpm2pK4 = 0.0000000567;                  ///<value>Stefan-Boltzmann constant (MJ m^-2 K^-4 day^-1)</value>
        const double Emissivity_surface_frac = 0.97;                             ///<value>surface emissivity for longwave radiation</value>
        const double Emissivity_surface_98 = 0.98;                          ///<value>surface emissivity (M.C. Llasat, R.L. Snyder, Ag and Forest Meteorology, 1998)</value>
        const double CP = 1013;                             ///<value>specific heat of moist air (J/Kg C)</value>
        const int ZCONS = 2;                                ///<value>height of weather station (m)</value>
        const double rd = 0.00137;                          ///<value>roughness height for water (m)</value>
        const double rDS = 0.005;                           ///<value>roughness height for snow (m)</value>
        const double rDT = 0.95;                            ///<value>roughness height for tree (m)</value>
        const double rTreeHeight = 5;                       ///<value>tree height (m)</value>
        const double rGroundHeight_m = 0.1;                  ///<value>ground height (m)</value>
        const double rZom = 0.123;                          ///<value>roughness length coefficient controlling momentum transfer defined Eq 13.3 Chin (2021)</value>
        const double rZov = 0.0123;                         ///<value>roughness length coefficient controlling heat and vapor transfer defined Eq 13.4 and following text Chin (2021)</value>
        const double rZu = ZCONS;                           ///<value>wind measurement height (m) from constants</value>
        const double rZe = ZCONS;                           ///<value>vapor measurement height (m) from constants</value>
        const double rZuT = rTreeHeight + rZu;              ///<value>wind estimate height in tree (m)</value>
        const double rZuG = 0.1 + rd;                       ///<value>wind estimate height for ground (m)</value>
        const double rZuG_m = rGroundHeight_m + rZu;                       ///<value>wind estimate height for ground (m)</value>

        public const double LEAF_STORAGE_M = 0.0002;               ///<value>specific leaf storage of water = 0.2 mm = 0.0002 m</value>
        const double IMPERV_STORAGE_M = 0.0015;             ///<value>maximum impervious cover storage of water = 1.5 mm = 0.0015 m</value>
        const double PERV_STORAGE_M = 0.001;                ///<value>maximum pervious cover storage of water = 1 mm = 0.001 m</value>                                    

        public const double NODATA = -999;
        const int NOWBAN = 99999;
        public const string SFC_TABLE = "SurfaceWeather";
        const string HRMIX_TABLE = "HourlyMixHt";

        const string NEW_INTERNATIONAL = "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW ZZ ZZ ZZ W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD";
        const string NEW_US_CANADA = "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB MW MW MW MW AW AW AW AW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD";
        const string ISD_old = "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD";
        const string NARCCAP = "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01      NETRAD PCPXX SD";
        // Field number for the old surface weather format
        const int USAF = 0;
        const int WBAN = 1;
        const int DATETIME = 2;
        const int DIR = 3;
        const int SPD = 4;
        const int GUS = 5;
        const int CLG = 6;
        const int SKC = 7;
        const int TEMP = 16;
        const int DEWP = 17;
        const int SLP = 18;
        const int ALT = 19;
        const int STP = 20;
        const int MAX = 21;
        const int MIN = 22;
        const int PCP01 = 23;
        const int NETRAD = 24;
        const int SD = 27;

        public DateTime TimeStamp;
        public int Jday;
        public int Year;
        public int Month;
        public int Day;
        public int Hour;
        public int Minute;
        public double Radiation_Longwave_Downwelling_Wpm2;
        public double Radiation_Longwave_Upwelling_Wpm2;
        public double LatRad;
        public double HourAngleRad;
        public double DecAngleRad;
        public double AirDens_kg_p_m3;
        public double airMass_relativeOptical_kg_p_m2;
        public double Altimeter_inHg;
        public double Altimeter_kPa;
        public double Tair_2m_C;
        public double Tair_2m_F;
        public double Tair_2m_K;
        public double Tdew_2m_C;
        public double Tdew_2m_F;
        public double Radiation_Shortwave_Diffuse_Wpm2;
        public double Radiation_Shortwave_Direct_Wpm2;
        public double Radiation_Shortwave_Total_Wpm2;
        public double Radiation_Net_Wpm2;
        public double AtmP_hundredths_inHg;
        public double AtmP_kPa;
        public double AtmP_mbar;
        public double PeGrMh;
        public double PeTrMh;
        public double PeSnGrMh;
        public double PeSnTrMh;
        public double PtTrMh;
        public double VegEvMh;                  ///<value>Evaporation from vegetations (m/h)</value>
        public double VegStMh;                  ///<value>Storage of precipitation by vegetation (m/h)</value>
        public double VegIntcptMh;              ///<value>Precipitation interception by vegetation (m/h)</value>
        public double UnderCanThrufallMh;       ///<value>Under canopy throughfall of precipitation (m/h)</value>
        public double UnderCanPervEvMh;         ///<value>Evaporation from under canopy pervious cover (m/h)</value>
        public double UnderCanPervStMh;         ///<value>Storage of precipitation by under canopy pervious cover (m/h)</value>                                       
        public double UnderCanPervInfilMh;      ///<value>Infiltration from under canopy pervious cover (m/h)</value>
        public double UnderCanImpervEvMh;       ///<value>Evaporation from under canopy impervious cover (m/h)</value>
        public double UnderCanImpervStMh;       ///<value>Storage of precipitation by under canopy impervious cover (m/h)</value>                                       
        public double UnderCanImpervRunoffMh;   ///<value>Runoff from under canopy impervious cover (m/h)</value>
        public double NoCanPervEvMh;            ///<value>Evaporation from no canopy pervious cover (m/h)</value>
        public double NoCanPervStMh;            ///<value>Storage of precipitation by no canopy pervious cover (m/h)</value>                                       
        public double NoCanPervInfilMh;         ///<value>Infiltration from no canopy pervious cover (m/h)</value>
        public double NoCanImpervEvMh;          ///<value>Evaporation from no canopy impervious cover (m/h)</value>
        public double NoCanImpervStMh;          ///<value>Storage of precipitation by no canopy impervious cover (m/h)</value>                                       
        public double NoCanImpervRunoffMh;      ///<value>Runoff from no canopy impervious cover (m/h)</value>
        public double PARuEm2s;
        public double PARWm2;
        public double Ppt_1hr_inches;
        public double Ppt_6hr_inches;
        public double Ppt_1hr_m;
        public double Ppt_6hr_m;
        public double Ppt_24hr_m;
        public double RelHum;
        public double VapPres_Sat_kPa;
        public double SnowDepth_inches;
        public double SnowDepth_m;
        public double SolZenAgl_deg;
        public double CloudCover_Total_Tenth;
        public double CloudCover_Opaque_tenths;
        public double CloudCover_Translucent_tenths;
        public double CeilingHeight_hundredths_ft;
        public double VapPres_Act_kPa;
        public double Wdir_deg;
        public double Wspd_kt;
        public double Wspd_mph;
        public double Wspd_mps;

        public double CloudCover_1_Tenth;      // 0–10, NODATA if missing
        public double CloudCover_2_Tenth;      // 0–10
        public double CloudCover_3_Tenth;      // 0–10

        public double CloudBaseHeight_1_m;    // meters, NODATA if missing
        public double CloudBaseHeight_2_m;
        public double CloudBaseHeight_3_m;
        public double CloudBaseHeight_m;

        public bool Flag_SkyData_hasMultipleLayers;

        public enum WeatherDataFormat { GHCNh_csv, ISD_old, ISD_intl, ISD_us_can, NARCCAP };
        public static WeatherDataFormat WeatherData_Format;
        #endregion

        //ProcessSurfaceWeatherData receives variables from CreateOutput.cs
        public static void ProcessSurfaceWeatherData(string MeteorologicalDataInputFile_string, string PrecipitationDataInputFile_string, LocationData LocationData_Object,string LeafAreaIndex_dataBase, VEGTYPE VegetationType_class, int DateYear_Start_int, int DateYear_Stop_int, string Meteo_DatabasePath_string, double Height_WeatherStationWindSensor_m, double Height_TreeCanopyTop_m, ref List<SurfaceWeather> SurfaceWeather_List_final)
        {
            List<SurfaceWeather> SurfaceWeather_List_raw = new List<SurfaceWeather>();
            List<SurfaceWeather> SurfaceWeather_List_Hr = new List<SurfaceWeather>();
            List<LeafAreaIndex> LeafAreaIndex_list = new List<LeafAreaIndex>();
            int Variable_TimeSeries_Hr_Length = 0;

            try {
                //WeatherData_Format = CheckWeatherDataFormat(MeteorologicalDataInputFile_string), as ISD or GHCNh
                WeatherData_Format = CheckWeatherDataFormat(MeteorologicalDataInputFile_string);

                //If WeatherData_Format == WeatherDataFormat.GHCNh_csv is true then call ReadSurfaceWeatherData_GHCNh
                if (WeatherData_Format == WeatherDataFormat.GHCNh_csv) {
                    //Note: GHCNh header and data are formatted as: YR--MODAHRMN, temperature_C, dew_point_temperature_C, station_level_pressure_hPa, sea_level_pressure_hPa, altimeter_hPa, wind_direction_deg, wind_speed_mps, wind_gust_mps, precipitation_mm, precipitation_6_hour_mm, precipitation_24_hour_mm, sky_cover_baseht_1_m, sky_cover_1_code, sky_cover_baseht_2_m, sky_cover_2_code, sky_cover_baseht_3_m, sky_cover_3_code
                    ReadSurfaceWeatherData_GHCNh(MeteorologicalDataInputFile_string, ref SurfaceWeather_List_raw, DateYear_Start_int);
                }
                //Else call ReadSurfaceWeatherData_ISD for legacy formats
                else {
                    ReadSurfaceWeatherData_ISD(MeteorologicalDataInputFile_string, ref SurfaceWeather_List_raw, DateYear_Start_int);
                }

                //AdjustTimeStamp will convert observations from off hour to on hour when that hour is missing, e.g., 3:54 to 4:00
                AdjustTimeStamp(SurfaceWeather_List_raw);

                //RawData_Index_Start to RawData_Index_Stop will check for actual length of record
                int RawData_Index_Start = -1, RawData_Index_Stop = -1;
                //For loop through SurfaceWeather_List_raw.Count to find index of first and last record
                for (int i = 0; i < SurfaceWeather_List_raw.Count; i++) {
                    //If SurfaceWeather_List_raw[i].Tair_2m_C is a number (or not Not a Number) then enter
                    if (!double.IsNaN(SurfaceWeather_List_raw[i].Tair_2m_C) || !double.IsNaN(SurfaceWeather_List_raw[i].Tdew_2m_C)) {
                        //If RawData_Index_Start == -1 then at start of list, convert to i
                        if (RawData_Index_Start == -1) RawData_Index_Start = i;
                        //RawData_Index_Stop updatd to i repeatedly, untile last instance of Tair_2m_C is a number
                        RawData_Index_Stop = i;
                    }
                }

                DateTime DateYear_Start = new DateTime(DateYear_Start_int, 1, 1, 0, 0, 0); 
                DateTime DateYear_Stop = new DateTime(DateYear_Stop_int, 12, 31, 23, 0, 0);

                //RawData_DateTime_Start and RawData_DateTime_Stop created as fallbacks if nothing found
                DateTime RawData_DateTime_Start = (RawData_Index_Start == -1) ? new DateTime(DateYear_Start_int, 1, 1, 0, 0, 0)
                                                      : SurfaceWeather_List_raw[RawData_Index_Start].TimeStamp.Date;
                DateTime RawData_DateTime_Stop = (RawData_Index_Stop == -1) ? new DateTime(DateYear_Stop_int, 12, 31, 23, 0, 0)
                                                      : SurfaceWeather_List_raw[RawData_Index_Stop].TimeStamp.Date.AddHours(23);

                //RawData_DateTime_Start clamped to start to at least DateYear_Start
                if (RawData_DateTime_Start < DateYear_Start) {RawData_DateTime_Start = DateYear_Start; }
                //RawData_DateTime_Stop clamped to start to at least DateYear_Stop
                if (RawData_DateTime_Stop > DateYear_Stop) { RawData_DateTime_Stop = DateYear_Stop; }

                //CreateHourlyRecords converts SurfaceWeather_List_raw to SurfaceWeather_List_Hr, and gets length of list
                CreateHourlyRecords(SurfaceWeather_List_raw, ref SurfaceWeather_List_Hr, RawData_DateTime_Start, RawData_DateTime_Stop);

                //FillExtrapolate will extrapolate values to the end of the observation period
                FillExtrapolate(SurfaceWeather_List_Hr);
                //FillInterpolate will interpolate values for missing hours
                FillInterpolate(SurfaceWeather_List_Hr);

                //Variable_TimeSeries_Hr_Length is SurfaceWeather_List_Hr.count
                Variable_TimeSeries_Hr_Length = SurfaceWeather_List_Hr.Count;

                //AdjustTimeZone converts from UTC to the GMTOffset of the station location, using i-Tree Database
                AdjustTimeZone(SurfaceWeather_List_Hr, LocationData_Object.GMTOffset, Variable_TimeSeries_Hr_Length);

                //FillTimeZoneAdjustmentDataGap function fills missing hours caused by the time zone adjustment
                FillTimeZoneAdjustmentDataGap(SurfaceWeather_List_Hr, Variable_TimeSeries_Hr_Length, ref SurfaceWeather_List_final);

                //Variable_TimeSeries_Hr_Length is SurfaceWeather_List_final.count
                Variable_TimeSeries_Hr_Length = SurfaceWeather_List_final.Count;

                //If WeatherPrepConfig.xml has entry for PrecipitationDataCsv then update Ppt_1hr_m variable in SurfaceWeather_List_final 
                if (!string.IsNullOrEmpty(PrecipitationDataInputFile_string)) {
                    //ReplaceRainData function sent the PrecipitationDataInputFile_string file from WeatherPrepConfig.xml
                    ReplaceRainData(PrecipitationDataInputFile_string, ref SurfaceWeather_List_final);
                }

                CalcAirDensity(SurfaceWeather_List_final, Variable_TimeSeries_Hr_Length);
                CalcSolarZenithAngle(SurfaceWeather_List_final, LocationData_Object, Variable_TimeSeries_Hr_Length);
                CalcSolarRadiation(SurfaceWeather_List_final, LocationData_Object, Variable_TimeSeries_Hr_Length);
                LeafAreaIndex.ReadLAIPartialRecords(LeafAreaIndex_dataBase, ref LeafAreaIndex_list, SurfaceWeather_List_final[0].TimeStamp, SurfaceWeather_List_final[Variable_TimeSeries_Hr_Length - 1].TimeStamp);
                CalcET(SurfaceWeather_List_final, LeafAreaIndex_list, Variable_TimeSeries_Hr_Length, Height_WeatherStationWindSensor_m, Height_TreeCanopyTop_m);
                CalcPrecipInterceptByCanopy(SurfaceWeather_List_final, LeafAreaIndex_list, Variable_TimeSeries_Hr_Length, VegetationType_class);
                CalcPrecipInterceptByUnderCanopyCover(SurfaceWeather_List_final, Variable_TimeSeries_Hr_Length);
                CalcPrecipInterceptByNoCanopyCover(SurfaceWeather_List_final, Variable_TimeSeries_Hr_Length);
            }
            catch (Exception) {
                throw;
            }
        }
        //ReplaceRainData function uses XML-config provided rainFile to replace raw NOAA precipitation record
        //ReplaceRainData function uses XML-config provided rainFile to merge/replace raw NOAA precipitation record
        static void ReplaceRainData(string rainFile, ref List<SurfaceWeather> list)
        {
            List<SurfaceWeather> rList = new List<SurfaceWeather>();

            // Call ReadPrecip to populate rList from PrecipitationDataCsv
            ReadPrecip(rainFile, ref rList);

            if (rList == null || rList.Count == 0) {
                Console.WriteLine("Warning: PrecipitationDataFile appears to be empty. " +
                                  "SurfaceWeather precipitation (Ppt_1hr_m) was not modified.");
                return;
            }

            // Sort precipitation list by time, just to be safe
            rList = rList.OrderBy(r => r.TimeStamp).ToList();

            // Build a lookup: TimeStamp -> Ppt_1hr_m (m/h)
            Dictionary<DateTime, double> precipByTime = new Dictionary<DateTime, double>();
            DateTime firstPrecipTs = rList[0].TimeStamp;
            DateTime lastPrecipTs = rList[rList.Count - 1].TimeStamp;

            foreach (var rec in rList) {
                // If duplicates exist, the last one wins
                precipByTime[rec.TimeStamp] = rec.Ppt_1hr_m;
            }

            int nMatch = 0;
            int nNoMatch = 0;
            HashSet<DateTime> usedPrecipTimes = new HashSet<DateTime>();

            // Merge by timestamp: replace Ppt_1hr_m when a Giovanni/PrecipitationDataCsv value exists
            foreach (var wx in list) {
                if (precipByTime.TryGetValue(wx.TimeStamp, out double ppt_m)) {
                    wx.Ppt_1hr_m = ppt_m;
                    nMatch++;
                    usedPrecipTimes.Add(wx.TimeStamp);
                }
                else {
                    // Keep whatever was in list[i].Ppt_1hr_m (often NOAA-derived or zero)
                    nNoMatch++;
                }
            }

            int nUniquePrecipTimes = precipByTime.Count;
            int nUnusedPrecipTimes = nUniquePrecipTimes - usedPrecipTimes.Count;

            Console.WriteLine();
            Console.WriteLine("------------------------------------------------------------");
            Console.WriteLine("PrecipitationDataFile merge summary:");
            Console.WriteLine("  Weather records total     : {0}", list.Count);
            Console.WriteLine("  Weather records updated   : {0}", nMatch);
            Console.WriteLine("  Weather records unchanged : {0}", nNoMatch);

            // Report the precip file temporal extent
            Console.WriteLine();
            Console.WriteLine("  First precipitation time  : {0}", firstPrecipTs.ToString("yyyy-MM-dd HH:mm:ss"));
            Console.WriteLine("  Last  precipitation time  : {0}", lastPrecipTs.ToString("yyyy-MM-dd HH:mm:ss"));

            if (nUnusedPrecipTimes > 0) {
                Console.WriteLine("  Precipitation times not used (no matching weather record): {0}", nUnusedPrecipTimes);
            }

            // Strong warning when not all weather rows were covered
            if (nNoMatch > 0) {
                Console.WriteLine("************************************************************");
                Console.WriteLine(" ***  WARNING: INCOMPLETE PRECIPITATION COVERAGE  ***");
                Console.WriteLine("************************************************************");
                Console.WriteLine("{0} WeatherPrep records did NOT have matching timestamps in", nNoMatch);
                Console.WriteLine("the PrecipitationDataFile (PrecipitationDataCsv).");
                Console.WriteLine();
                Console.WriteLine("These records retain their original precipitation values,");
                Console.WriteLine("which may be zeros or placeholders, and can cause under-");
                Console.WriteLine("estimation of total precipitation for the affected period.");
                Console.WriteLine();
                Console.WriteLine("Check that the WeatherPrepConfig.xml date range aligns with");
                Console.WriteLine("the precipitation file coverage, or obtain additional");
                Console.WriteLine("precipitation data if needed.");
                Console.WriteLine("************************************************************");
            }
            else {
                Console.WriteLine("All WeatherPrep records had matching precipitation timestamps.");
                Console.WriteLine("Precipitation timeseries from PrecipitationDataFile fully");
                Console.WriteLine("replaced the SurfaceWeatherDataFile precipitation.");
            }
            Console.WriteLine("------------------------------------------------------------");
            Console.WriteLine();
        }
        //ReadPrecip is used to parse Date, Time and Precipitation data from the optional Precip.csv input file for replacing NOAA raw weather values
        //JK 20211029: Renamed from "ReadCsv" to be more precise, and to avoid confusion with the functions handling the new CSV files for LAI & weather inputs
        static void ReadPrecip(string sFile, ref List<SurfaceWeather> rList)
        {

            // ------------------------------------------------------------
            // NEW: Check that the file exists before trying to read it
            // ------------------------------------------------------------
            if (string.IsNullOrWhiteSpace(sFile) || !File.Exists(sFile)) {
                Console.WriteLine();
                Console.WriteLine("************************************************************");
                Console.WriteLine(" ***  WARNING: PRECIPITATION DATA FILE NOT FOUND  ***");
                Console.WriteLine("************************************************************");
                Console.WriteLine("WeatherPrepConfig.xml specifies a PrecipitationDataCsv file:");
                Console.WriteLine("    {0}", sFile);
                Console.WriteLine();
                Console.WriteLine("But this file does NOT exist. SurfaceWeather precipitation");
                Console.WriteLine("fields (Ppt_1hr_m) will NOT be replaced.");
                Console.WriteLine();
                Console.WriteLine("Please verify the file path in <PrecipitationDataCsv>");
                Console.WriteLine("and ensure the file is accessible to WeatherPrep.");
                Console.WriteLine("************************************************************");
                Console.WriteLine();
                return;
            }
            string line;
            SurfaceWeather data;
            Regex reg;
            Match m;

            // Updated regex to handle scientific notation in precipitation values
            reg = new Regex(@"^\s*(?<dt>\d{8}),\s*(?<hr>\d{1,2}):(?<mn>\d{2}):(?<sc>\d{2}),\s*(?<p>[0-9.]+(?:[eE][+-]?[0-9]+)?)\s*$");

            // sr created as StreamReader type of sFile passed to ReadPrecip
            using (StreamReader sr = new StreamReader(sFile))
            {
                try {
                    // Read the header line
                    sr.ReadLine();

                    while ((line = sr.ReadLine()) != null) {
                        // Match the line using the updated regex
                        m = reg.Match(line);

                        // If the match is successful, continue processing
                        if (m.Success) {
                            // Create a new SurfaceWeather object to hold the data
                            data = new SurfaceWeather();

                            // Parse date and time information
                            data.Year = int.Parse(m.Groups["dt"].ToString().Substring(0, 4));
                            data.Month = int.Parse(m.Groups["dt"].ToString().Substring(4, 2));
                            data.Day = int.Parse(m.Groups["dt"].ToString().Substring(6, 2));
                            data.Hour = int.Parse(m.Groups["hr"].ToString());
                            data.Minute = int.Parse(m.Groups["mn"].ToString());

                            // Create a timestamp from the parsed date and time
                            data.TimeStamp = DateTime.Parse($"{data.Month}/{data.Day}/{data.Year} {data.Hour}:{data.Minute}");

                            // Parse the precipitation value, handling both decimal and scientific notation
                            data.Ppt_1hr_m = double.Parse(m.Groups["p"].ToString());

                            // Add the parsed data to the result list
                            rList.Add(data);
                        }
                        else {
                            // Handle the case where the line doesn't match the expected pattern
                            Console.WriteLine($"Warning: Line did not match the expected format: {line}");
                        }
                    }
                }
                catch (Exception ex) {
                    // Catch and re-throw any exceptions, or log them
                    Console.WriteLine("Error processing file: " + ex.Message);
                    throw;
                }
            }
        }

        /// <summary>
        /// Read WeatherPrep-ready CSV produced from GHCN-Hourly.
        /// Expected columns (order can vary):
        ///   YR--MODAHRMN, temperature_C, dew_point_temperature_C, station_level_pressure_hPa, altimeter_hPa, wind_direction_deg, wind_speed_mps, precipitation_mm, precipitation_6_hour_mm, precipitation_24_hour_mm, snow_depth_mm, sky_cover_baseht_1_m, sky_cover_1_code, sky_cover_baseht_2_m, sky_cover_2_code, sky_cover_baseht_3_m, sky_cover_3_code
        /// </summary>
        static void ReadSurfaceWeatherData_GHCNh(string sFile, ref List<SurfaceWeather> WeatherVariables_list, int StartYear_input)
        {
            using (StreamReader sr = new StreamReader(sFile)) {
                string header = sr.ReadLine();
                if (string.IsNullOrWhiteSpace(header)) {
                    Console.WriteLine("Error: Empty header in GHCN CSV file: " + sFile);
                    return;
                }

                string[] headers = header.Split(',');

                // Helper: get column index by name (case-insensitive)
                int idx(string name) {
                    for (int i = 0; i < headers.Length; i++) {
                        if (string.Equals(headers[i].Trim(), name, StringComparison.OrdinalIgnoreCase))
                            return i;
                    }
                    return -1;
                }

                // Datetime (YYYYMMDDHHMM)
                int idx_DateTime_YRMM = idx("YR--MODAHRMN");

                // Temperatures (°C)
                int idx_Temperature_C = idx("temperature_C");
                int idx_DewPoint_C = idx("dew_point_temperature_C");

                // Pressures (hPa)
                int idx_StationPressure_hPa = idx("station_level_pressure_hPa");
                int idx_SeaLevelPressure_hPa = idx("sea_level_pressure_hPa");
                int idx_Altimeter_hPa = idx("altimeter_hPa");

                // Wind (deg, m/s)
                int idx_WindDirection_deg = idx("wind_direction_deg");
                int idx_WindSpeed_mps = idx("wind_speed_mps");
                int idx_WindGust_mps = idx("wind_gust_mps");

                // Precipitation (mm)
                int idx_Precip_1hr_mm = idx("precipitation_mm");
                int idx_Precip_6hr_mm = idx("precipitation_6_hour_mm");
                int idx_Precip_24hr_mm = idx("precipitation_24_hour_mm");
                int idx_SnowDepth_mm = idx("snow_depth_mm");

                // Cloud base heights (m)
                int idx_SkyCoverBaseHt1_m = idx("sky_cover_baseht_1_m");
                int idx_SkyCoverBaseHt2_m = idx("sky_cover_baseht_2_m");
                int idx_SkyCoverBaseHt3_m = idx("sky_cover_baseht_3_m");

                // Raw sky cover codes (alphanumeric)
                int idx_SkyCoverCode1 = idx("sky_cover_1_code");
                int idx_SkyCoverCode2 = idx("sky_cover_2_code");
                int idx_SkyCoverCode3 = idx("sky_cover_3_code");

                if (idx_DateTime_YRMM < 0)
                {
                    Console.WriteLine("Error: GHCN CSV is missing YR--MODAHRMN column: " + sFile);
                    return;
                }

                string line;
                int counter = 0;

                // Helper: safe string access
                string field(string[] f, int index)
                {
                    if (index < 0 || index >= f.Length) return string.Empty;
                    return (f[index] ?? string.Empty).Trim();
                }

                //parseOrNoData Helper for safe double parsing with NODATA default
                double parseOrNoData(string[] f, int index)
                {
                    if (index < 0 || index >= f.Length) return NODATA;
                    string s = (f[index] ?? string.Empty).Trim();
                    if (string.IsNullOrEmpty(s)) return NODATA;
                    if (double.TryParse(s, out double val)) return val;
                    return NODATA;
                }

                int debugRow = 0;
                while ((line = sr.ReadLine()) != null) {

                    debugRow++;

                    if (string.IsNullOrWhiteSpace(line)) {
                        continue;
                    }

                    string[] f = line.Split(',');
                    if (f.Length <= idx_DateTime_YRMM) {
                        continue;
                    }

                    SurfaceWeather data = new SurfaceWeather();

                    string dt = f[idx_DateTime_YRMM].Trim(); // YYYYMMDDHHMM
                    if (dt.Length < 10) {
                        continue;
                    }

                    data.Year = int.Parse(dt.Substring(0, 4));
                    data.Month = int.Parse(dt.Substring(4, 2));
                    data.Day = int.Parse(dt.Substring(6, 2));
                    data.Hour = int.Parse(dt.Substring(8, 2));
                    data.Minute = (dt.Length >= 12) ? int.Parse(dt.Substring(10, 2)) : 0;

                    data.TimeStamp = new DateTime(data.Year, data.Month, data.Day, data.Hour, data.Minute, 0);

                    // First-record consistency check with StartYear_input, similar to legacy reader
                    counter++;
                    if (counter == 1 && StartYear_input != data.Year)
                    {
                        Console.WriteLine("\nWarning: The WeatherPrepConfig.xml StartYear is " + StartYear_input +
                                          " and not the same as the input weather file year of " + data.Year + ".");
                        Console.WriteLine("Exiting: WeatherPrep cannot run when the desired year in WeatherPrepConfig.xml is different from the actual data.");
                        Console.WriteLine("Suggestion: Make sure the StartYear in WeatherPrepConfig.xml matches the year of your input data.");
                        Environment.Exit(1);
                    }

                    // --- Temperatures (convert °C -> °F) ---
                    double Tair_C = parseOrNoData(f, idx_Temperature_C);
                    double Tdew_C = parseOrNoData(f, idx_DewPoint_C);

                    //Tair_2m_C air temperature at 2 m height, degrees C
                    if (Tair_C != NODATA) {
                        data.Tair_2m_C = Tair_C;
                        //Tair_2m_K air temperature at 2 m height, degrees K
                        data.Tair_2m_K = Tair_C + 273.15;
                        //Tair_2m_F air temperature at 2 m height, degrees F
                        data.Tair_2m_F = Tair_C * 9.0 / 5.0 + 32.0;
                    }
                    else {
                        data.Tair_2m_C = NODATA;
                        data.Tair_2m_K = NODATA;
                        data.Tair_2m_F = NODATA;
                    }
                    if (Tdew_C != NODATA)
                    {
                        data.Tdew_2m_C = Tdew_C;
                        data.Tdew_2m_F = Tdew_C * 9.0 / 5.0 + 32.0;
                    }
                    else
                    {
                        data.Tdew_2m_C = NODATA;
                        data.Tdew_2m_F = NODATA;
                    }
                    // --- Pressure (hPa ≈ mbar) ---
                    double stnP_hPa = parseOrNoData(f, idx_StationPressure_hPa);
                    double slp_hPa = parseOrNoData(f, idx_SeaLevelPressure_hPa);

                    // Prefer station pressure; fall back to sea-level if needed
                    double AtmP_hPa = (stnP_hPa != NODATA) ? stnP_hPa : slp_hPa;
                    data.AtmP_kPa = (AtmP_hPa != NODATA) ? AtmP_hPa / 10.0 : NODATA;
                    data.AtmP_mbar = (AtmP_hPa != NODATA) ? AtmP_hPa : NODATA;
                    //AtmP_hundredths_inHg  (hundredths of inHg) is (Atm_mbar/33.8639) * 100
                    data.AtmP_hundredths_inHg = (AtmP_hPa != NODATA) ? (AtmP_hPa / 33.8639) * 100.0 : NODATA;

                    // --- Altimeter_inHg (hPa -> kPa) ---
                    double alt_hPa = parseOrNoData(f, idx_Altimeter_hPa);
                    if (alt_hPa != NODATA) {
                        // 1 hPa ≈ 0.1 kPa
                        data.Altimeter_kPa = alt_hPa / 10.0;
                    }
                    else {
                        data.Altimeter_kPa = NODATA;
                    }

                    // --- Wind (deg, m/s -> mph) ---
                    data.Wdir_deg = parseOrNoData(f, idx_WindDirection_deg);

                    double ws_mps = parseOrNoData(f, idx_WindSpeed_mps);
                    //If ws_mps not NODATA then 
                    if (ws_mps != NODATA) {
                        //Wspd_mps (meters per second)
                        data.Wspd_mps = ws_mps;
                        //Wspd_mph (miles per hour)
                        data.Wspd_mph = ws_mps * 2.23694;
                        //Wspd_kt (knot)
                        data.Wspd_kt = ws_mps * 1.94384;
                    }
                    //Else NODATA
                    else {
                        data.Wspd_mps = NODATA;
                        data.Wspd_mph = NODATA;
                        data.Wspd_kt = NODATA;
                    }

                    // Gust is currently unused in SurfaceWeather, but we could parse it if needed:
                    // double gust_mps = parseOrNoData(f, idx_WindGust_mps);

                    // --- Precipitation (mm -> inches) ---
                    double ppt1_mm = parseOrNoData(f, idx_Precip_1hr_mm);
                    double ppt6_mm = parseOrNoData(f, idx_Precip_6hr_mm);
                    double ppt24_mm = parseOrNoData(f, idx_Precip_24hr_mm);

                    //Ppt_1hr_m Precipitation (meters, per hour), converting mm to m
                    data.Ppt_1hr_m = (ppt1_mm != NODATA) ? ppt1_mm / 1000.0 : NODATA;
                    //Ppt_6hr_m Precipitation (meters, per hour), converting mm to m
                    data.Ppt_6hr_m = (ppt6_mm != NODATA) ? ppt6_mm / 1000.0 : NODATA;
                    //Ppt_24hr_m Precipitation (meters, per hour), converting mm to m
                    data.Ppt_24hr_m = (ppt24_mm != NODATA) ? ppt24_mm / 1000.0 : NODATA;

                    // --- Sky cover: use last non-empty layer (3 → 2 → 1) ---
                    // --- Sky cover: 3 layers --- 
                    string code1 = field(f, idx_SkyCoverCode1);
                    string code2 = field(f, idx_SkyCoverCode2);
                    string code3 = field(f, idx_SkyCoverCode3);

                    double baseHt1_m = parseOrNoData(f, idx_SkyCoverBaseHt1_m);
                    double baseHt2_m = parseOrNoData(f, idx_SkyCoverBaseHt2_m);
                    double baseHt3_m = parseOrNoData(f, idx_SkyCoverBaseHt3_m);

                    // Map to okta 0–8 (NOT tenths) so we preserve raw coverage
                    data.CloudCover_1_Tenth = MapGHCNSkyCodeToTenths(code1, NODATA);
                    data.CloudCover_2_Tenth = MapGHCNSkyCodeToTenths(code2, NODATA);
                    data.CloudCover_3_Tenth = MapGHCNSkyCodeToTenths(code3, NODATA);

                    data.CloudBaseHeight_1_m = baseHt1_m;
                    data.CloudBaseHeight_2_m = baseHt2_m;
                    data.CloudBaseHeight_3_m = baseHt3_m;

                    // Map to tenths for each layer
                    data.CloudCover_1_Tenth = MapGHCNSkyCodeToTenths(code1, NODATA);
                    data.CloudCover_2_Tenth = MapGHCNSkyCodeToTenths(code2, NODATA);
                    data.CloudCover_3_Tenth = MapGHCNSkyCodeToTenths(code3, NODATA);

                    data.CloudBaseHeight_1_m = baseHt1_m;
                    data.CloudBaseHeight_2_m = baseHt2_m;
                    data.CloudBaseHeight_3_m = baseHt3_m;

                    // Determine if we have multiple layers
                    bool hasLayer1 = (data.CloudCover_1_Tenth != NODATA && data.CloudBaseHeight_1_m != NODATA);
                    bool hasLayer2 = (data.CloudCover_2_Tenth != NODATA && data.CloudBaseHeight_2_m != NODATA);
                    bool hasLayer3 = (data.CloudCover_3_Tenth != NODATA && data.CloudBaseHeight_3_m != NODATA);

                    // Store Flag_SkyData_hasMultipleLayers for later use (e.g., CalcNetRadWm2, CalcCloudCover)
                    data.Flag_SkyData_hasMultipleLayers = ((hasLayer1 && hasLayer2) ||
                                                            (hasLayer1 && hasLayer3) ||
                                                            (hasLayer2 && hasLayer3));

                    // GHCNh guidance: total state of sky is best determined by the *last* non-empty layer
                    string Sky_Code_effective = string.Empty;
                    double CloudBaseHt_effective_m = NODATA;

                    if (!string.IsNullOrEmpty(code3)) {
                        Sky_Code_effective = code3;
                        CloudBaseHt_effective_m = baseHt3_m;
                    }
                    else if (!string.IsNullOrEmpty(code2)) {
                        Sky_Code_effective = code2;
                        CloudBaseHt_effective_m = baseHt2_m;
                    }
                    else if (!string.IsNullOrEmpty(code1)) {
                        Sky_Code_effective = code1;
                        CloudBaseHt_effective_m = baseHt1_m;
                    }

                    // Total cover in tenths is from that effective code
                    if (!string.IsNullOrEmpty(Sky_Code_effective))
                        data.CloudCover_Total_Tenth = MapGHCNSkyCodeToTenths(Sky_Code_effective, NODATA);
                    else
                        data.CloudCover_Total_Tenth = NODATA;

                    // CeilingHeight_hundredths_ft uses the same effective base height (ISD & GHCNh consistent)
                    if (CloudBaseHt_effective_m != NODATA) {
                        double effectiveBaseHt_ft = CloudBaseHt_effective_m * 3.28084;
                        data.CeilingHeight_hundredths_ft = effectiveBaseHt_ft / 100.0;
                    }
                    else {
                        data.CeilingHeight_hundredths_ft = NODATA;
                    }

                    // --- Snow depth (mm in GHCNh; convert to meters)
                    double snow_mm = parseOrNoData(f, idx_SnowDepth_mm);

                    if (snow_mm != NODATA) {
                        // convert mm → m
                        data.SnowDepth_m = snow_mm / 1000.0;
                    }
                    else {
                        data.SnowDepth_m = NODATA;
                    }

                    // Preserve SnowDepth_inches for compatibility with ISD
                    data.SnowDepth_inches = (snow_mm != NODATA)
                                            ? (snow_mm / 25.4)  // mm → inches
                                            : 0.0;

                    WeatherVariables_list.Add(data);
                }
            }
        }

        //MapGHCNSkyCodeToTenths uses the GHCNh sky_cover_1 variables to return cloud cover values of tenths 
        //Note: GHCNh reports: Since up to 3 cloud layers can be reported, the full state of the sky can best be determined by the last layer's value. In other words if three layers are reported and the third layer uses BKN then the total state of sky is BKN which is similar in definition to “mostly cloudy.” OVC is similar to “cloudy” or overcast and FEW or SCT is similar to “partly cloudy.” In cases where there are more than 3 cloud layers, the highest layers will not be reported.
        //Values in oktas:
        //CLR:00 None, SKC or CLR
        //FEW:01 One okta - 1/10 or less but not zero
        //FEW:02 Two oktas - 2/10 - 3/10, or FEW
        //SCT:03 Three oktas - 4/10
        //SCT:04 Four oktas - 5/10, or SCT
        //BKN:05 Five oktas - 6/10
        //BKN:06 Six oktas - 7/10 - 8/10
        //BKN:07 Seven oktas - 9/10 or more but not 10/10, or BKN
        //OVC:08 Eight oktas - 10/10, or OVC
        //VV:09 Sky obscured, or cloud amount cannot be estimated
        //X:10 Partial obscuration
        //Note: The reported oktas for each level are presumed to be the cloud coverage at that level, even if the view was blocked by lower cloud cover and the higher cloud cover fraction could be less.
        //Note: The obscured condition is interpreted as full cloud cover. 
        static double MapGHCNSkyCodeToTenths(string code, double NODATA)
        {
            if (string.IsNullOrWhiteSpace(code))
                return NODATA;

            code = code.Trim().ToUpperInvariant();

            // NEW: handle combined codes like "OVC:08" or "BKN:07"
            int colonIndex = code.IndexOf(':');
            if (colonIndex >= 0)
            {
                // e.g. "OVC"
                string left = code.Substring(0, colonIndex).Trim();
                // e.g. "08"
                string right = code.Substring(colonIndex + 1).Trim();

                // Prefer the numeric piece if it’s present; otherwise use the CloudBase_normalized_frac.
                if (!string.IsNullOrWhiteSpace(right))
                    code = right;
                else
                    code = left;
            }

            // Now 'code' is like "OVC" or "08" or "BKN" or "07"

            switch (code)
            {
                case "CLR":
                case "SKC":
                case "00":
                    return 0.0;

                case "FEW":
                case "01":
                    return 1.0;

                case "02":
                    return 2.5;

                case "SCT":
                case "03":
                    return 4.0;

                case "04":
                    return 5.0;

                case "05":
                    return 6.0;

                case "BKN":
                case "06":
                    return 7.5;

                case "07":
                    return 8.5;

                case "OVC":
                case "08":
                    return 10.0;

                case "VV":
                case "09":
                    return 10.0;

                case "X":
                case "10":
                    return 10.0;

                default:
                    return NODATA;
            }
        }

        //CompressCloudLayers_Tenths will compress the 3 layer cloud heights into an effective single height
        //Note: CloudBase_effective_km is a weighted mean of reported cloud heights (no clamping)
        //Note: CloudCover_frac is the maximum cloud cover fraction across layers, per GHCNh guidance
        public static void CompressCloudLayers_Tenths(
            double CloudCover_1_tenths, double CloudBase_1_m,
            double CloudCover_2_tenths, double CloudBase_2_m,
            double CloudCover_3_tenths, double CloudBase_3_m,
            double NODATA,
            out double CloudCover_frac,
            out double CloudBase_effective_km)
        {
            double[] CloudCover_vec_frac = { CloudCover_1_tenths, CloudCover_2_tenths, CloudCover_3_tenths };
            double[] CloudBase_vec_m = { CloudBase_1_m, CloudBase_2_m, CloudBase_3_m };

            double CloudCover_sum_frac = 0.0;
            double CloudBase_CoverWeighted_sum_km = 0.0;

            double CloudCover_max_frac = 0.0;
            bool Flag_CloudCover_Present = false;

            for (int i = 0; i < 3; i++)
            {
                if (CloudCover_vec_frac[i] == NODATA || CloudBase_vec_m[i] == NODATA)
                {
                    continue;
                }

                double CloudCover_i_frac = CloudCover_vec_frac[i] / 10.0;  // 0–1 fraction
                if (CloudCover_i_frac <= 0.0)
                {
                    continue;
                }

                Flag_CloudCover_Present = true;

                // GHCNh recommended total cover is max across layers
                if (CloudCover_i_frac > CloudCover_max_frac)
                {
                    CloudCover_max_frac = CloudCover_i_frac;
                }

                double CloudBase_i_km = CloudBase_vec_m[i] / 1000.0;     // km, NO clamping here
                CloudCover_sum_frac += CloudCover_i_frac;
                CloudBase_CoverWeighted_sum_km += CloudCover_i_frac * CloudBase_i_km;
            }

            if (!Flag_CloudCover_Present)
            {
                CloudCover_frac = 0.0;
                CloudBase_effective_km = NODATA;
                return;
            }

            CloudCover_frac = CloudCover_max_frac;

            if (CloudCover_sum_frac > 0.0)
            {
                CloudBase_effective_km = CloudBase_CoverWeighted_sum_km / CloudCover_sum_frac;
            }
            else
            {
                CloudBase_effective_km = NODATA;
            }
        }

        //WeatherDataFormat CheckWeatherDataFormat defines weather input data format
        //Note: GHCNh simply distinguished by extension .csv, while ISD arrives as .txt in multiple formats
        public static WeatherDataFormat CheckWeatherDataFormat(string MeteorologicalDataInputFile_string)
        {
            string line;
            WeatherDataFormat rc = WeatherDataFormat.ISD_old;

            //Note: GHCNh header and data are formatted as: YR--MODAHRMN, temperature_C, dew_point_temperature_C, station_level_pressure_hPa, sea_level_pressure_hPa, altimeter_hPa, wind_direction_deg, wind_speed_mps, wind_gust_mps, precipitation_mm, precipitation_6_hour_mm, precipitation_24_hour_mm, sky_cover_baseht_1_m, sky_cover_1_code, sky_cover_baseht_2_m, sky_cover_2_code, sky_cover_baseht_3_m, sky_cover_3_code
            if (Path.GetExtension(MeteorologicalDataInputFile_string).Equals(".csv", StringComparison.OrdinalIgnoreCase)) {
                return WeatherDataFormat.GHCNh_csv;
            }

            using (StreamReader sr = new StreamReader(MeteorologicalDataInputFile_string)) {
                try {
                    //read the header line
                    line = sr.ReadLine();
                    //Console.WriteLine("\nMeteorological Input Data Header: " + line + "\n");

                    if (Regex.IsMatch(line, "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW ZZ ZZ ZZ W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD"))
                    {
                        rc = WeatherDataFormat.ISD_intl;
                    }
                    else if (Regex.IsMatch(line, "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB MW MW MW MW AW AW AW AW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD"))
                    {
                        rc = WeatherDataFormat.ISD_us_can;
                    }
                    else if (Regex.IsMatch(line, "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01 PCP06 PCP24 PCPXX SD"))
                    {
                        rc = WeatherDataFormat.ISD_old;
                    }
                    else if (Regex.IsMatch(line, "  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB WW WW WW W TEMP DEWP    SLP   ALT    STP MAX MIN PCP01      NETRAD PCPXX SD"))
                    {
                        rc = WeatherDataFormat.NARCCAP;
                    }
                    else
                    {
                        //throw new InvalidOperationException("\nWarning: WeatherPrepConfig.xml seems to point to a faulty SurfaceWeatherDataFile.");
                        Console.WriteLine("\nWarning: The WeatherPrepConfig.xml file contains a faulty value for SurfaceWeatherDataFile.");
                        Console.WriteLine("Check: " + MeteorologicalDataInputFile_string);
                        Console.WriteLine("Confirm: The header should contain  USAF  WBAN YR--MODAHRMN DIR SPD GUS CLG SKC L M H  VSB ...");
                        Console.WriteLine("Suggestion: Files downloaded from NOAA often need to be processed by ishapp2.exe.");
                    }
                    return rc;
                }
                catch (Exception)
                {
                    throw;
                }
            }
        }

        //ReadSurfaceWeatherData_ISD function reads in the raw meteorological data
        //Note: Presumes data arrives in BG units, and is then converted to SI units to match default GHCNh units
        public static void ReadSurfaceWeatherData_ISD(string sFile, ref List<SurfaceWeather> WeatherVariables_list, int StartYear_input)
        {
            string line;
            SurfaceWeather data;
            Regex reg;

            using (StreamReader sr = new StreamReader(sFile)) {
                try {
                    switch (WeatherData_Format) {
                        case WeatherDataFormat.ISD_intl:
                            reg = new Regex(@"^(?<usaf>.{6}) (?<wban>.{5}) (?<dt>.{12}) (?<dir>.{3}) (?<spd>.{3}) .{3} (?<clg>.{3}) (?<skc>.{3}) .{30} (?<temp>.{4}) (?<dewp>.{4}) .{6} (?<alt>.{5}) (?<stp>.{6}) .{7} (?<pcp01>.{5}) (?<pcp06>.{5}) .{11} (?<sd>.{2})\s*$");
                            break;
                        case WeatherDataFormat.ISD_us_can:
                            reg = new Regex(@"^(?<usaf>.{6}) (?<wban>.{5}) (?<dt>.{12}) (?<dir>.{3}) (?<spd>.{3}) .{3} (?<clg>.{3}) (?<skc>.{3}) .{36} (?<temp>.{4}) (?<dewp>.{4}) .{6} (?<alt>.{5}) (?<stp>.{6}) .{7} (?<pcp01>.{5}) (?<pcp06>.{5}) .{11} (?<sd>.{2})\s*$");
                            break;
                        case WeatherDataFormat.ISD_old:
                            reg = new Regex(@"^(?<usaf>.{6}) (?<wban>.{5}) (?<dt>.{12}) (?<dir>.{3}) (?<spd>.{3}) .{3} (?<clg>.{3}) (?<skc>.{3}) .{21} (?<temp>.{4}) (?<dewp>.{4}) .{6} (?<alt>.{5}) (?<stp>.{6}) .{7} (?<pcp01>.{5}) (?<pcp06>.{5}) .{11} (?<sd>.{2})\s*$");
                            break;
                        case WeatherDataFormat.NARCCAP:
                            reg = new Regex(@"^(?<usaf>.{6}) (?<wban>.{5}) (?<dt>.{12}) (?<dir>.{3}) (?<spd>.{3}) .{3} (?<clg>.{3}) (?<skc>.{3}) .{21} (?<temp>.{4}) (?<dewp>.{4}) .{6} (?<alt>.{5}) (?<stp>.{6}) .{7} (?<pcp>.{5}) (?<netrad>.{11}) .{5} (?<sd>.{2})\s*$");
                            break;
                        default:
                            reg = new Regex(@"^(?<usaf>.{6}) (?<wban>.{5}) (?<dt>.{12}) (?<dir>.{3}) (?<spd>.{3}) .{3} (?<clg>.{3}) (?<skc>.{3}) .{21} (?<temp>.{4}) (?<dewp>.{4}) .{6} (?<alt>.{5}) (?<stp>.{6}) .{7} (?<pcp01>.{5}) (?<pcp06>.{5}) .{11}  (?<sd>.{2})\s*$");
                            break;
                    }
                    Match m;
                    //read the header line
                    sr.ReadLine();

                    //CounterWhileLoop_int initialized to zero for error handling below, when CounterWhileLoop_int == 1
                    int CounterWhileLoop_int = 0;

                    while ((line = sr.ReadLine()) != null)
                    {
                        //CounterWhileLoop_int is advanced, and when CounterWhileLoop_int == 1 it will be used for error handling
                        CounterWhileLoop_int = CounterWhileLoop_int + 1;

                        line = Regex.Replace(line, "([0-9.]+)T", @"$1 ");   // remove trace (T) for PCP01
                        m = reg.Match(line);

                        data = new SurfaceWeather();

                        data.Year = (m.Groups["dt"].ToString().Substring(0, 4) == "****" || m.Groups["dt"].ToString().Substring(0, 4) == "    ") ? (int)NODATA : int.Parse(m.Groups["dt"].ToString().Substring(0, 4));

                        //If CounterWhileLoop_int equals 1 and the StartYear_input is not equal to data.Year then send error message
                        if (CounterWhileLoop_int == 1 & StartYear_input != data.Year)
                        {
                            //Notify user that the input instructions in the WeatherPrepConfigFile and weather data are problematic
                            Console.WriteLine("\nWarning: The WeatherPrepConfig.xml StartYear is " + StartYear_input + " and not the same as the input weather file year of " + data.Year + ".");
                            Console.WriteLine("Exiting: WeatherPrep cannot run when the desired year in WeatherPrepConfig.xml is different from the actual data.");
                            Console.WriteLine("Suggestion: Make sure the StartYear in WeatherPrepConfig.xml matches the year of your input data.");
                            //exit program when this error emerges
                            System.Environment.Exit(1);
                        }

                        data.Month = (m.Groups["dt"].ToString().Substring(4, 2) == "**" || m.Groups["dt"].ToString().Substring(4, 2) == "  ") ? (int)NODATA : int.Parse(m.Groups["dt"].ToString().Substring(4, 2));
                        data.Day = (m.Groups["dt"].ToString().Substring(6, 2) == "**" || m.Groups["dt"].ToString().Substring(6, 2) == "  ") ? (int)NODATA : int.Parse(m.Groups["dt"].ToString().Substring(6, 2));
                        data.Hour = (m.Groups["dt"].ToString().Substring(8, 2) == "**" || m.Groups["dt"].ToString().Substring(8, 2) == "  ") ? (int)NODATA : int.Parse(m.Groups["dt"].ToString().Substring(8, 2));
                        data.Minute = (m.Groups["dt"].ToString().Substring(10, 2) == "**" || m.Groups["dt"].ToString().Substring(10, 2) == "  ") ? (int)NODATA : int.Parse(m.Groups["dt"].ToString().Substring(10, 2));
                        data.TimeStamp = DateTime.Parse(data.Month + "/" +
                                                        data.Day + "/" +
                                                        data.Year + " " +
                                                        data.Hour + ":" +
                                                        data.Minute);
                        data.Wdir_deg = (m.Groups["dir"].ToString() == "***" || m.Groups["dir"].ToString() == "   ") ? NODATA : double.Parse(m.Groups["dir"].ToString());
                        //read in wind spead (mph)
                        data.Wspd_mph = (m.Groups["spd"].ToString() == "***" || m.Groups["spd"].ToString() == "   ") ? NODATA : double.Parse(m.Groups["spd"].ToString());
                        //read in cloud ceiling (lowest opaque layer w/ 5/8 or greater coverage (hundreds ft); 722 = unlimited
                        data.CeilingHeight_hundredths_ft = (m.Groups["clg"].ToString() == "***" || m.Groups["clg"].ToString() == "   ") ? NODATA : double.Parse(m.Groups["clg"].ToString());
                        //read in sky cover (CLR-CLEAR, SCT-SCATTERED-1/8 TO 4/8, BKN-BROKEN-5/8 TO 7/8, OVC-OVERCAST, OBS-OBSCURED, POB-PARTIAL OBSCURATION
                        switch (m.Groups["skc"].ToString())
                        {
                            case "***":
                            case "   ":
                                data.CloudCover_Total_Tenth = NODATA;
                                break;
                            case "CLR":
                                data.CloudCover_Total_Tenth = 0;
                                break;
                            case "SCT":
                                //Note: CloudCover_Total_Tenth = 3.125 is average of SKC = SCT = average(1/8, 2/8, 3/8,4/8) = 2.5/8 (not 3.75)
                                data.CloudCover_Total_Tenth = 3.125;
                                break;
                            case "BKN":
                                data.CloudCover_Total_Tenth = 7.5;
                                break;
                            case "OVC":
                                data.CloudCover_Total_Tenth = 10;
                                break;
                            default:
                                data.CloudCover_Total_Tenth = 10;
                                break;
                        }
                        //read in temperature (F)
                        data.Tair_2m_F = (m.Groups["temp"].ToString() == "****" || m.Groups["temp"].ToString() == "    ") ? NODATA : double.Parse(m.Groups["temp"].ToString());
                        //read in dew point (F)
                        data.Tdew_2m_F = (m.Groups["dewp"].ToString() == "****" || m.Groups["dewp"].ToString() == "    ") ? NODATA : double.Parse(m.Groups["dewp"].ToString());
                        //read altimeter setting (inches Hg)
                        data.Altimeter_inHg = (m.Groups["alt"].ToString() == "*****" || m.Groups["alt"].ToString() == "     ") ? NODATA : double.Parse(m.Groups["alt"].ToString());
                        //read in station atmospheric pressure (mbar)
                        data.AtmP_mbar = (m.Groups["stp"].ToString() == "******" || m.Groups["stp"].ToString() == "      ") ? NODATA : double.Parse(m.Groups["stp"].ToString());
                        //read in liquid precipitation for 1-hr (inches)
                        data.Ppt_1hr_inches = (m.Groups["pcp01"].ToString() == "*****" || m.Groups["pcp01"].ToString() == "     ") ? 0 : double.Parse(m.Groups["pcp01"].ToString());
                        if (WeatherData_Format == WeatherDataFormat.NARCCAP) {
                            //read in net radiation (W/m2)
                            data.Radiation_Net_Wpm2 = (m.Groups["netrad"].ToString() == "***********" || m.Groups["netrad"].ToString() == "           ") ? NODATA : double.Parse(m.Groups["netrad"].ToString());
                        }
                        else {
                            //read in liquid precipitation for 6-hr (inches)
                            data.Ppt_6hr_inches = (m.Groups["pcp06"].ToString() == "*****" || m.Groups["pcp06"].ToString() == "     ") ? 0 : double.Parse(m.Groups["pcp06"].ToString());
                        }
                        //read in snow depth (inches)
                        data.SnowDepth_inches = (m.Groups["sd"].ToString() == "**" || m.Groups["sd"].ToString() == "  ") ? 0 : double.Parse(m.Groups["sd"].ToString());

                        //Conversion to SI
                        //If data.Tair_2m_F != NODATA) then convert
                        if (data.Tair_2m_F != NODATA) {
                            data.Tair_2m_C = (data.Tair_2m_F - 32.0) * (5.0 / 9.0);
                            data.Tair_2m_K = data.Tair_2m_C + ConversionFactor_273pt15K_to_0C;
                        }
                        //Else assign to NODATA
                        else {
                            data.Tair_2m_C = NODATA;
                            data.Tair_2m_K = NODATA;
                        }

                        //If data.Tdew_2m_F != NODATA) then convert
                        if (data.Tdew_2m_F != NODATA) {
                            data.Tdew_2m_C = (data.Tdew_2m_F - 32.0) * (5.0 / 9.0);
                        }
                        else {
                            data.Tdew_2m_C = NODATA;
                        }

                        // ---- Pressure: mbar → kPa + legacy inHg ----
                        if (data.AtmP_mbar != NODATA) {
                            data.AtmP_kPa = data.AtmP_mbar / 10.0;
                            //AtmP_hundredths_inHg  (hundredths of inHg) is (Atm_mbar/33.8639) * 100
                            data.AtmP_hundredths_inHg = (data.AtmP_mbar / 33.8639) * 100.0;
                        }
                        else {
                            data.AtmP_kPa = NODATA;
                            data.AtmP_hundredths_inHg = NODATA;
                        }

                        if (data.Altimeter_inHg != NODATA) {
                            // Convert inHg → kPa
                            data.Altimeter_kPa = data.Altimeter_inHg * 3.386389;
                        }
                        else {
                            data.Altimeter_kPa = NODATA;
                        }
                        
                        // ---- Wind: mph → m/s (+ kt) ----
                        if (data.Wspd_mph != NODATA) {
                            //Wspd_mps (meters per second) from Wspd_mph (miles per hour) conversion
                            data.Wspd_mps = data.Wspd_mph * 0.44704;
                            //Wspd_kt (knot) converted from Wpsd_mps
                            data.Wspd_kt = data.Wspd_mps * 1.94384;
                        }
                        else {
                            data.Wspd_mps = NODATA;
                            data.Wspd_kt = NODATA;
                        }

                        // ---- Precipitation: inches → m ----
                        // Precip: inches → m (0 or positive)
                        data.Ppt_1hr_m = data.Ppt_1hr_inches * 0.0254;

                        // Snow depth: inches → m (0 or positive)
                        data.SnowDepth_m = data.SnowDepth_inches * 0.0254;

                        // ---- Cloud height: hundreds of ft → m ----
                        if (data.CeilingHeight_hundredths_ft != NODATA) {
                            // CLG is hundreds of ft; you stored hundredths-of-ft, so:
                            //Note: This converts the ISD code 722 (indicating unlimited) to 22,000 m, which is unlimited
                            double base_ft = data.CeilingHeight_hundredths_ft * 100.0;
                            data.CloudBaseHeight_1_m = base_ft * 0.3048;
                            //CloudBaseHeight_2_m and CloudBaseHeight_3_m are NODATA for ISD
                            data.CloudBaseHeight_2_m = NODATA;
                            data.CloudBaseHeight_3_m = NODATA;
                        }
                        else {
                            data.CloudBaseHeight_1_m = NODATA;
                            //CloudBaseHeight_2_m and CloudBaseHeight_3_m are NODATA for ISD
                            data.CloudBaseHeight_2_m = NODATA;
                            data.CloudBaseHeight_3_m = NODATA;
                        }

                        // ---- Cloud cover: tenths, 1 layer to 3 layers ----
                        if (data.CloudCover_Total_Tenth != NODATA) {
                            data.CloudCover_1_Tenth = data.CloudCover_Total_Tenth;
                            //CloudCover_2_Tenth and CloudCover_3_Tenth are NODATA for ISD
                            data.CloudCover_2_Tenth = NODATA;
                            data.CloudCover_3_Tenth = NODATA;
                        }
                        else {
                            data.CloudCover_1_Tenth = NODATA;
                            //CloudCover_2_Tenth and CloudCover_3_Tenth are NODATA for ISD
                            data.CloudCover_2_Tenth = NODATA;
                            data.CloudCover_3_Tenth = NODATA;
                        }
                        WeatherVariables_list.Add(data);
                    }
                }
                catch (Exception)
                {
                    throw;
                }
            }
        }

        public static void FillExtrapolate(List<SurfaceWeather> SurfaceWeather_List_Hr)
        {
            if (SurfaceWeather_List_Hr == null || SurfaceWeather_List_Hr.Count == 0)
                return;

            // Station pressure (kPa)
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.AtmP_kPa,
                (sw, v) => sw.AtmP_kPa = v,
                varLabel: "AtmP_kPa");

            // Air temperature (C)
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.Tair_2m_C,
                (sw, v) => sw.Tair_2m_C = v,
                varLabel: "Tair_2m_C");

            // Dewpoint temperature (C)
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.Tdew_2m_C,
                (sw, v) => sw.Tdew_2m_C = v,
                varLabel: "Tdew_2m_C");

            // Wind direction
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.Wdir_deg,
                (sw, v) => sw.Wdir_deg = v,
                varLabel: "Wdir_deg");

            // Wind speed
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.Wspd_mps,
                (sw, v) => sw.Wspd_mps = v,
                varLabel: "Wspd_mps");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_1_m,
                (sw, v) => sw.CloudBaseHeight_1_m = v,
                varLabel: "CloudBaseHeight_1_m");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_2_m,
                (sw, v) => sw.CloudBaseHeight_2_m = v,
                varLabel: "CloudBaseHeight_2_m");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_3_m,
                (sw, v) => sw.CloudBaseHeight_3_m = v,
                varLabel: "CloudBaseHeight_3_m");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_1_Tenth,
                (sw, v) => sw.CloudCover_1_Tenth = v,
                varLabel: "CloudCover_1_Tenth");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_2_Tenth,
                (sw, v) => sw.CloudCover_2_Tenth = v,
                varLabel: "CloudCover_2_Tenth");

            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_3_Tenth,
                (sw, v) => sw.CloudCover_3_Tenth = v,
                varLabel: "CloudCover_3_Tenth");

            // Total cloud cover
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_Total_Tenth,
                (sw, v) => sw.CloudCover_Total_Tenth = v,
                varLabel: "SKC",
                allNoDataFill: 3.75);

            // Ceiling height
            ExtrapolateEnds(
                SurfaceWeather_List_Hr,
                sw => sw.CeilingHeight_hundredths_ft,
                (sw, v) => sw.CeilingHeight_hundredths_ft = v,
                varLabel: "CLG",
                allNoDataFill: 100.0);

            // Net radiation (NARCCAP only)
            if (WeatherData_Format == WeatherDataFormat.NARCCAP)
            {
                ExtrapolateEnds(
                    SurfaceWeather_List_Hr,
                    sw => sw.Radiation_Net_Wpm2,
                    (sw, v) => sw.Radiation_Net_Wpm2 = v,
                    varLabel: "Radiation_Net_Wpm2");
            }
        }

        private static void ExtrapolateEnds(
            List<SurfaceWeather> SurfaceWeather_List_Hr,
            Func<SurfaceWeather, double> getter,
            Action<SurfaceWeather, double> setter,
            string varLabel,
            double? allNoDataFill = null)
        {
            if (SurfaceWeather_List_Hr == null || SurfaceWeather_List_Hr.Count == 0)
                return;

            // Find first non-NODATA
            int first = -1;
            for (int i = 0; i < SurfaceWeather_List_Hr.Count; i++)
            {
                if (getter(SurfaceWeather_List_Hr[i]) != NODATA)
                {
                    first = i;
                    break;
                }
            }

            // No valid values at all
            if (first == -1)
            {
                if (allNoDataFill.HasValue)
                {
                    Console.WriteLine(
                        $"Warning: SurfaceWeatherDataFile only has NODATA values for {varLabel}. " +
                        $"Converting all values to {allNoDataFill.Value}.");
                    Console.WriteLine(
                        $"Suggestion: Replace SurfaceWeatherDataFile with one containing observations for {varLabel}.");

                    for (int i = 0; i < SurfaceWeather_List_Hr.Count; i++)
                    {
                        setter(SurfaceWeather_List_Hr[i], allNoDataFill.Value);
                    }
                }
                return;
            }

            // Fill leading NODATA up to first non-NODATA
            double firstVal = getter(SurfaceWeather_List_Hr[first]);
            for (int i = 0; i < first; i++)
            {
                setter(SurfaceWeather_List_Hr[i], firstVal);
            }

            // Find last non-NODATA
            int last = -1;
            for (int i = SurfaceWeather_List_Hr.Count - 1; i >= 0; i--)
            {
                if (getter(SurfaceWeather_List_Hr[i]) != NODATA)
                {
                    last = i;
                    break;
                }
            }

            // Fill trailing NODATA after last non-NODATA
            double lastVal = getter(SurfaceWeather_List_Hr[last]);
            for (int i = last + 1; i < SurfaceWeather_List_Hr.Count; i++)
            {
                setter(SurfaceWeather_List_Hr[i], lastVal);
            }
        }

        public static void FillInterpolate(List<SurfaceWeather> SurfaceWeather_List_Hr)
        {
            if (SurfaceWeather_List_Hr == null || SurfaceWeather_List_Hr.Count == 0)
                return;

            // Linear fillers
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.AtmP_kPa,
                (sw, v) => sw.AtmP_kPa = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.Tair_2m_C,
                (sw, v) => sw.Tair_2m_C = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.Tdew_2m_C,
                (sw, v) => sw.Tdew_2m_C = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.Wdir_deg,
                (sw, v) => sw.Wdir_deg = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.Wspd_mps,
                (sw, v) => sw.Wspd_mps = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_1_m,
                (sw, v) => sw.CloudBaseHeight_1_m = v);
           
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_2_m,
                (sw, v) => sw.CloudBaseHeight_2_m = v);
            
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudBaseHeight_3_m,
                (sw, v) => sw.CloudBaseHeight_3_m = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_1_Tenth,
                (sw, v) => sw.CloudCover_1_Tenth = v);
            
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_2_Tenth,
                (sw, v) => sw.CloudCover_2_Tenth = v);
            
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_3_Tenth,
                (sw, v) => sw.CloudCover_3_Tenth = v);
            
            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CeilingHeight_hundredths_ft,
                (sw, v) => sw.CeilingHeight_hundredths_ft = v);

            InterpolateLinear(
                SurfaceWeather_List_Hr,
                sw => sw.CloudCover_Total_Tenth,
                (sw, v) => sw.CloudCover_Total_Tenth = v);

            if (WeatherData_Format == WeatherDataFormat.NARCCAP)
            {
                InterpolateLinear(
                    SurfaceWeather_List_Hr,
                    sw => sw.Radiation_Net_Wpm2,
                    (sw, v) => sw.Radiation_Net_Wpm2 = v);
            }
        }

        private static void InterpolateLinear(
            List<SurfaceWeather> SurfaceWeather_List_Hr,
            Func<SurfaceWeather, double> getter,
            Action<SurfaceWeather, double> setter)
        {
            if (SurfaceWeather_List_Hr == null || SurfaceWeather_List_Hr.Count == 0)
                return;

            int i = 0;
            while (i < SurfaceWeather_List_Hr.Count)
            {
                // Look for a run of NODATA starting at i
                if (getter(SurfaceWeather_List_Hr[i]) == NODATA)
                {
                    int bfIdx = i - 1;
                    if (bfIdx < 0)
                    {
                        // Can't interpolate if we don't have a "before" value;
                        // edge gaps are handled by ExtrapolateEnds.
                        i++;
                        continue;
                    }

                    int j = i + 1;
                    // Find next non-NODATA
                    while (j < SurfaceWeather_List_Hr.Count && getter(SurfaceWeather_List_Hr[j]) == NODATA)
                        j++;

                    if (j >= SurfaceWeather_List_Hr.Count)
                    {
                        // No "after" value to anchor interpolation; edges are
                        // handled by ExtrapolateEnds, so we stop here.
                        break;
                    }

                    int afIdx = j;

                    // Now interpolate between bfIdx and afIdx
                    DateTime tBefore = SurfaceWeather_List_Hr[bfIdx].TimeStamp;
                    DateTime tAfter = SurfaceWeather_List_Hr[afIdx].TimeStamp;
                    double vBefore = getter(SurfaceWeather_List_Hr[bfIdx]);
                    double vAfter = getter(SurfaceWeather_List_Hr[afIdx]);

                    double totalMinutes = (tAfter - tBefore).TotalMinutes;
                    if (totalMinutes <= 0)
                    {
                        // Pathological case; just copy vBefore
                        for (int k = bfIdx + 1; k < afIdx; k++)
                            setter(SurfaceWeather_List_Hr[k], vBefore);
                    }
                    else
                    {
                        for (int k = bfIdx + 1; k < afIdx; k++)
                        {
                            double minutesSinceBefore =
                                (SurfaceWeather_List_Hr[k].TimeStamp - tBefore).TotalMinutes;
                            double frac = minutesSinceBefore / totalMinutes;
                            double val = vBefore + frac * (vAfter - vBefore);
                            setter(SurfaceWeather_List_Hr[k], val);
                        }
                    }

                    i = afIdx; // continue from after the filled segment
                }
                else
                {
                    i++;
                }
            }
        }

        public static void FillTimeZoneAdjustmentDataGap(List<SurfaceWeather> SurfaceWeather_List_Hr, int Variable_TimeSeries_Hr_Length, ref List<SurfaceWeather> SurfaceWeather_List_final)
        {
            if (SurfaceWeather_List_Hr == null || SurfaceWeather_List_Hr.Count == 0 || Variable_TimeSeries_Hr_Length <= 0) {
                SurfaceWeather_List_final = new List<SurfaceWeather>();
                return;
            }

            int startIdx = 0;
            int endIdx = Math.Min(Variable_TimeSeries_Hr_Length, SurfaceWeather_List_Hr.Count) - 1;

            // Enforce first 00:00 and last 23:00 if needed
            if (SurfaceWeather_List_Hr[startIdx].Hour != 0) {
                for (int i = 0; i < Variable_TimeSeries_Hr_Length; i++) {
                    if (SurfaceWeather_List_Hr[i].Hour == 0) { startIdx = i; break; }
                }
            }

            if (SurfaceWeather_List_Hr[endIdx].Hour != 23) {
                for (int i = endIdx; i >= startIdx; i--) {
                    if (SurfaceWeather_List_Hr[i].Hour == 23) { endIdx = i; break; }
                }
            }

            // Optional: clamp trailing all-NODATA hours
            int lastDataIdx = endIdx;
            for (int i = endIdx; i >= startIdx; i--) {
                if (HasAnyObservation(SurfaceWeather_List_Hr[i], SurfaceWeather.NODATA)) {
                    lastDataIdx = i;
                    break;
                }
            }
            endIdx = Math.Max(lastDataIdx, startIdx);

            // Slice but reuse original objects
            int count = endIdx - startIdx + 1;
            SurfaceWeather_List_final = SurfaceWeather_List_Hr.GetRange(startIdx, count);

            //Pad missing Jan 1 / Dec 31 caused by GMT offset
            PadYearEdgesIfMissing(SurfaceWeather_List_final);
        }

        private static void PadYearEdgesIfMissing(List<SurfaceWeather> list)
        {
            if (list == null || list.Count < 48)
                return;

            var (jan1Missing, dec31Missing) = ComputeEdgePaddingFlags(list);

            // ---- Pad January 1 from January 2 ----
            if (jan1Missing) {
                Console.WriteLine("Notice: Due to the time zone for your Place being west of UTC, ...");
                Console.WriteLine("... January 1 data are incomplete and January 2 data were used.");

                DateTime jan1Date = list[0].TimeStamp.Date.AddDays(-1);

                // Assume we have full 24 hours for Jan 2 starting at index 0
                var jan1Block = new List<SurfaceWeather>(24);

                for (int hour = 0; hour < 24; hour++) {
                    var src = list[hour];
                    var clone = CloneSurfaceWeather(src);

                    // Set the new timestamp to Jan 1 local, hour = 0..23
                    clone.TimeStamp = jan1Date.AddHours(hour);
                    clone.Year = clone.TimeStamp.Year;
                    clone.Month = clone.TimeStamp.Month;
                    clone.Day = clone.TimeStamp.Day;
                    clone.Hour = clone.TimeStamp.Hour;

                    // Precip should be 0 so we don't over-count the year
                    clone.Ppt_1hr_m = 0.0;

                    // Snow depth kept from donor day (as in your current logic)
                    // clone.SnowDepth_m already copied from src

                    jan1Block.Add(clone);
                }

                // Insert Jan 1 at the start
                list.InsertRange(0, jan1Block);
            }

            // ---- Pad December 31 from December 30 ----
            if (dec31Missing && list.Count >= 24) {
                DateTime dec31Date = list[list.Count - 1].TimeStamp.Date.AddDays(1);

                var dec31Block = new List<SurfaceWeather>(24);
                int donorStart = list.Count - 24; // last full day (= Dec 30) before we append

                for (int hour = 0; hour < 24; hour++) {
                    var src = list[donorStart + hour];
                    var clone = CloneSurfaceWeather(src);

                    clone.TimeStamp = dec31Date.AddHours(hour);
                    clone.Year = clone.TimeStamp.Year;
                    clone.Month = clone.TimeStamp.Month;
                    clone.Day = clone.TimeStamp.Day;
                    clone.Hour = clone.TimeStamp.Hour;

                    // Again, no extra precip on the synthetic last day
                    clone.Ppt_1hr_m = 0.0;

                    dec31Block.Add(clone);
                }

                // Append Dec 31
                list.AddRange(dec31Block);
            }
        }

        private static (bool jan1Missing, bool dec31Missing)
        ComputeEdgePaddingFlags(List<SurfaceWeather> list)
        {
            bool enoughForEdgeCheck = list != null && list.Count >= 48;
            if (!enoughForEdgeCheck)
                return (false, false);

            var firstDate = list[0].TimeStamp.Date;
            var lastDate = list[list.Count - 1].TimeStamp.Date;

            bool jan1Missing =
                firstDate.Month == 1 &&
                firstDate.Day == 2;

            bool dec31Missing =
                lastDate.Month == 12 &&
                lastDate.Day == 30;

            return (jan1Missing, dec31Missing);
        }

        private static SurfaceWeather CloneSurfaceWeather(SurfaceWeather src)
        {
            return new SurfaceWeather
            {
                TimeStamp = src.TimeStamp,
                Jday = src.Jday,
                Year = src.Year,
                Month = src.Month,
                Day = src.Day,
                Hour = src.Hour,

                Tair_2m_C = src.Tair_2m_C,
                Tdew_2m_C = src.Tdew_2m_C,
                Tair_2m_F = src.Tair_2m_F,
                Tdew_2m_F = src.Tdew_2m_F,
                Tair_2m_K = src.Tair_2m_K,
                Wdir_deg = src.Wdir_deg,
                Wspd_mps = src.Wspd_mps,
                Wspd_kt = src.Wspd_kt,
                AtmP_kPa = src.AtmP_kPa,
                AtmP_hundredths_inHg = src.AtmP_hundredths_inHg,
                Ppt_1hr_m = src.Ppt_1hr_m,
                SnowDepth_m = src.SnowDepth_m,

                Radiation_Net_Wpm2 = src.Radiation_Net_Wpm2,
                CloudCover_Total_Tenth = src.CloudCover_Total_Tenth,
                CloudBaseHeight_1_m = src.CloudBaseHeight_1_m,
                CloudBaseHeight_2_m = src.CloudBaseHeight_2_m,
                CloudBaseHeight_3_m = src.CloudBaseHeight_3_m,
                CloudCover_1_Tenth = src.CloudCover_1_Tenth,
                CloudCover_2_Tenth = src.CloudCover_2_Tenth,
                CloudCover_3_Tenth = src.CloudCover_3_Tenth,

                // plus any other fields you rely on downstream
            };
        }


        /// <summary>
        /// If an observation on the hour is missing, round observation **up** to the next top-of-hour
        /// 
        /// Examples:
        ///   00:54 → 01:00  
        ///   06:00 → 06:00
        /// 
        /// Notes:
        /// - GHCNh and ISD reports often occur at off-hour times, such as :54, :07, :23, etc.
        /// - After rounding:
        ///      • PCP01: we retain the **maximum** PCP01 among obs that rounded into the same hour.
        ///      • All other variables: we retain the **last** observation within the rounded hour.
        /// </summary>
        public static void AdjustTimeStamp(List<SurfaceWeather> SurfaceWeather_List_raw)
        {
            foreach (var SurfaceWeather_record in SurfaceWeather_List_raw) {
                int TS_minute = SurfaceWeather_record.TimeStamp.Minute;

                if (TS_minute != 0) {
                    // Round UP to next hour
                    int TS_minute_gap = 60 - TS_minute;
                    SurfaceWeather_record.TimeStamp = SurfaceWeather_record.TimeStamp.AddMinutes(TS_minute_gap);
                }
            }
        }

        //CreateHourlyRecords converts SurfaceWeather_List_raw to SurfaceWeather_List_Hr, and gets length of list
        //Note: Will convert multiple observations per hour to a single observation; takes maximum precipitation
        //Note: Will trim calendar to raw observed dates, to avoid filling year with constant or empty values
        public static void CreateHourlyRecords(List<SurfaceWeather> SurfaceWeather_List_raw, ref List<SurfaceWeather> SurfaceWeather_List_Hr, DateTime RawData_DateTime_Start, DateTime RawData_DateTime_Stop)
        {
            if (SurfaceWeather_List_raw == null || SurfaceWeather_List_raw.Count == 0)
                return;

            if (RawData_DateTime_Stop < RawData_DateTime_Start)
                return;

            //Variable_TimeSeries_Hr_Length is total hours in timespan from start to stop date
            int Variable_TimeSeries_Hr_Length = (int)(RawData_DateTime_Stop - RawData_DateTime_Start).TotalHours + 1;

            // Initialize the hourly list
            SurfaceWeather_List_Hr = new List<SurfaceWeather>(Variable_TimeSeries_Hr_Length);
            for (int i = 0; i < Variable_TimeSeries_Hr_Length; i++) {
                var data = new SurfaceWeather();
                data.TimeStamp = RawData_DateTime_Start.AddHours(i);

                // Initialize key fields
                data.Wdir_deg = data.Wspd_mps = data.CloudBaseHeight_1_m = data.CloudBaseHeight_2_m = data.CloudBaseHeight_3_m = data.CloudCover_1_Tenth = data.CloudCover_2_Tenth = data.CloudCover_3_Tenth = data.CloudCover_Total_Tenth = data.Tair_2m_C = data.Tdew_2m_C = data.AtmP_kPa = data.Radiation_Net_Wpm2 = NODATA;

                data.Ppt_1hr_m = 0.0;
                data.SnowDepth_m = 0.0;

                // (Optional but cleaner: also init the related fields)
                data.Tair_2m_F = data.Tair_2m_K = NODATA;
                data.Tdew_2m_F = NODATA;

                SurfaceWeather_List_Hr.Add(data);
            }

            // Merge each SurfaceWeather_Variable_i_raw record into the appropriate hourly slot
            foreach (var SurfaceWeather_Variable_i_raw in SurfaceWeather_List_raw) {
                // Map timestamp to hourly index
                int idx = (int)(SurfaceWeather_Variable_i_raw.TimeStamp - RawData_DateTime_Start).TotalHours;
                if (idx < 0 || idx >= Variable_TimeSeries_Hr_Length)
                    continue;   // outside range, ignore

                var SurfaceWeather_Variable_i_Hr = SurfaceWeather_List_Hr[idx];

                MergeSurfaceWeatherRecord(SurfaceWeather_Variable_i_raw, SurfaceWeather_Variable_i_Hr);
            }
        }

        private static void MergeSurfaceWeatherRecord(SurfaceWeather SurfaceWeather_Variable_i_raw, SurfaceWeather SurfaceWeather_Variable_i_Hr)
        {
            // Precip: keep max over all reports in the hour
            if (SurfaceWeather_Variable_i_raw.Ppt_1hr_m >= SurfaceWeather_Variable_i_Hr.Ppt_1hr_m)
                SurfaceWeather_Variable_i_Hr.Ppt_1hr_m = SurfaceWeather_Variable_i_raw.Ppt_1hr_m;

            if (SurfaceWeather_Variable_i_raw.SnowDepth_m >= SurfaceWeather_Variable_i_Hr.SnowDepth_m)
                SurfaceWeather_Variable_i_Hr.SnowDepth_m = SurfaceWeather_Variable_i_raw.SnowDepth_m;

            // Wind
            if (SurfaceWeather_Variable_i_raw.Wdir_deg != NODATA)
                SurfaceWeather_Variable_i_Hr.Wdir_deg = SurfaceWeather_Variable_i_raw.Wdir_deg;
            if (SurfaceWeather_Variable_i_raw.Wspd_mps != NODATA)
                SurfaceWeather_Variable_i_Hr.Wspd_mps = SurfaceWeather_Variable_i_raw.Wspd_mps;

            // Cloud height 3 layers
            if (SurfaceWeather_Variable_i_raw.CloudBaseHeight_1_m != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudBaseHeight_1_m = SurfaceWeather_Variable_i_raw.CloudBaseHeight_1_m;
            if (SurfaceWeather_Variable_i_raw.CloudBaseHeight_2_m != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudBaseHeight_2_m = SurfaceWeather_Variable_i_raw.CloudBaseHeight_2_m;
            if (SurfaceWeather_Variable_i_raw.CloudBaseHeight_3_m != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudBaseHeight_3_m = SurfaceWeather_Variable_i_raw.CloudBaseHeight_3_m;
            // Cloud cover 3 layers
            if (SurfaceWeather_Variable_i_raw.CloudCover_1_Tenth != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudCover_1_Tenth = SurfaceWeather_Variable_i_raw.CloudCover_1_Tenth;
            if (SurfaceWeather_Variable_i_raw.CloudCover_2_Tenth != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudCover_2_Tenth = SurfaceWeather_Variable_i_raw.CloudCover_2_Tenth;
            if (SurfaceWeather_Variable_i_raw.CloudCover_3_Tenth != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudCover_3_Tenth = SurfaceWeather_Variable_i_raw.CloudCover_3_Tenth;
            // Cloud cover effective layer
            if (SurfaceWeather_Variable_i_raw.CloudCover_Total_Tenth != NODATA)
                SurfaceWeather_Variable_i_Hr.CloudCover_Total_Tenth = SurfaceWeather_Variable_i_raw.CloudCover_Total_Tenth;

            // Temperatures (now SI-first)
            if (SurfaceWeather_Variable_i_raw.Tair_2m_C != NODATA) {
                SurfaceWeather_Variable_i_Hr.Tair_2m_C = SurfaceWeather_Variable_i_raw.Tair_2m_C;
                SurfaceWeather_Variable_i_Hr.Tair_2m_K = SurfaceWeather_Variable_i_raw.Tair_2m_K;  // already computed from C in reader
                SurfaceWeather_Variable_i_Hr.Tair_2m_F = SurfaceWeather_Variable_i_raw.Tair_2m_F;
            }
            if (SurfaceWeather_Variable_i_raw.Tdew_2m_C != NODATA) {
                SurfaceWeather_Variable_i_Hr.Tdew_2m_C = SurfaceWeather_Variable_i_raw.Tdew_2m_C;
                SurfaceWeather_Variable_i_Hr.Tdew_2m_F = SurfaceWeather_Variable_i_raw.Tdew_2m_F;
            }

            // Pressure
            if (SurfaceWeather_Variable_i_raw.AtmP_kPa != NODATA)
                SurfaceWeather_Variable_i_Hr.AtmP_kPa = SurfaceWeather_Variable_i_raw.AtmP_kPa;

            // Radiation (if applicable)
            if (SurfaceWeather_Variable_i_raw.Radiation_Net_Wpm2 != NODATA)
                SurfaceWeather_Variable_i_Hr.Radiation_Net_Wpm2 = SurfaceWeather_Variable_i_raw.Radiation_Net_Wpm2;

            // If you later add more SI fields (e.g., VapPrsKpa, RelHum, etc),
            // this is the *only* place you need to extend.
        }

        /// <summary>
        /// Adjust Time Zone from UTC to station locaton
        public static void AdjustTimeZone(List<SurfaceWeather> SurfaceWeather_List_Hr, double timeZone, int Variable_TimeSeries_Hr_Length)
        {
            DateTime dateTimeLST;
            int i;

            try
            {
                for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
                {
                    dateTimeLST = SurfaceWeather_List_Hr[i].TimeStamp.AddHours((int)timeZone);    //GMT to LST

                    SurfaceWeather_List_Hr[i].TimeStamp = dateTimeLST;
                    SurfaceWeather_List_Hr[i].Year = dateTimeLST.Year;
                    SurfaceWeather_List_Hr[i].Month = dateTimeLST.Month;
                    SurfaceWeather_List_Hr[i].Day = dateTimeLST.Day;
                    SurfaceWeather_List_Hr[i].Hour = dateTimeLST.Hour;
                    SurfaceWeather_List_Hr[i].Jday = dateTimeLST.DayOfYear;

                }
            }
            catch (Exception)
            {
                throw;
            }
        }

        /// <summary>
        /// Returns true if this hourly record contains any real observation signal.
        /// Adjust fields if your project uses different members.
        /// </summary>
        private static bool HasAnyObservation(SurfaceWeather dataFile, double NODATA)
        {
            // Treat precipitation/snow > 0 as signal.
            // Treat typical met drivers != NODATA as signal (values can be zero legitimately).
            return
                dataFile.Tair_2m_C != NODATA ||
                dataFile.Tdew_2m_C != NODATA ||
                dataFile.AtmP_kPa != NODATA ||
                dataFile.Radiation_Net_Wpm2 != NODATA ||
                dataFile.Wspd_mps != NODATA ||
                dataFile.Wdir_deg != NODATA ||
                dataFile.CloudCover_Total_Tenth != NODATA ||
                dataFile.Ppt_1hr_m > 0.0 ||
                dataFile.SnowDepth_m > 0.0;
        }

        /// <summary>
        /// Calculate air density.
        /// </summary>
        /// <param name="WeatherVariables_list">List holding hourly surface weather data</param>
        /// <param name="Variable_TimeSeries_Hr_Length">Record count in WeatherVariables_list</param>
        public static void CalcAirDensity(List<SurfaceWeather> WeatherVariables_list, int Variable_TimeSeries_Hr_Length)
        {
            int i;

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                //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: Given equation uses Tair_2m_K, the 273.15 in the denominator is removed
                //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
                WeatherVariables_list[i].AirDens_kg_p_m3 = 3.486 * (WeatherVariables_list[i].AtmP_kPa / WeatherVariables_list[i].Tair_2m_K);
            }
        }
        /// <summary>
        /// Calculate solar zenith angle.
        /// </summary>
        /// <param name="WeatherVariables_list">List holding hourly surface weather data</param>
        /// <param name="LocationData_Object">Location data</param>
        /// <param name="Variable_TimeSeries_Hr_Length">Record count in WeatherVariables_list</param>
        public static void CalcSolarZenithAngle(List<SurfaceWeather> WeatherVariables_list, LocationData LocationData_Object, int Variable_TimeSeries_Hr_Length)
        {
            double dayAngle;  //day angle
            double solDecl; //solar declination
            double eqnTime;   //equation of time
            double stdMer;    //standard meridian
            double Latitude_rad;       //local latitude in radians
            double trSolTime; //true solar time in hours
            double hourAngle; //hour angle
            double cosZ;      //cosine of the zenith angle
            double exAtmZenAgl_deg;  //exoatmospheric solar zenith angle
            double exAtmElvAgl_deg; //exoatmospheric solar elevation angle
            double refCorr_deg;   //refraction correction
            int i;

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                //dayAngle (radians)
                dayAngle = Deg2Rad(360 * (WeatherVariables_list[i].Jday - 1) / 365.0);
                //eqnTime (decimal minutes), equation of time
                eqnTime = (0.000075 + 0.001868 * Math.Cos(dayAngle) - 0.032077 * Math.Sin(dayAngle)
                        - 0.014615 * Math.Cos(2 * dayAngle) - 0.040849 * Math.Sin(2 * dayAngle)) * 229.18;
                //solDecl (radians), solar declination angle
                solDecl = (0.006918 - 0.399912 * Math.Cos(dayAngle) + 0.070257 * Math.Sin(dayAngle)
                            - 0.006758 * Math.Cos(2 * dayAngle) + 0.000907 * Math.Sin(2 * dayAngle)
                            - 0.002697 * Math.Cos(3 * dayAngle) + 0.00148 * Math.Sin(3 * dayAngle));
                //WeatherVariables_list stores solDecl (radians)
                WeatherVariables_list[i].DecAngleRad = solDecl;
                //To handle both Northern/Southern hemisphere, accept negative latitude (for Southern)
                //Latitude_rad (radians), converted from LocationData_Object.Latitude (degrees)
                Latitude_rad = Deg2Rad(LocationData_Object.Latitude);
                //WeatherVariables_list stores latitude
                WeatherVariables_list[i].LatRad = Latitude_rad;
                //stdMer (degrees) takes 15 deg for each GMT offset hour
                stdMer = 15.0 * LocationData_Object.GMTOffset;
                //trSolTime (hour with decimal minutes), true solar time
                trSolTime = WeatherVariables_list[i].Hour + (4 * (Math.Abs(stdMer) - Math.Abs(LocationData_Object.Longitude)) + eqnTime) / 60.0;
                //hourAngle (degrees)
                hourAngle = 15.0 * trSolTime - 180;
                //hourAngle (radians)
                hourAngle = Deg2Rad(hourAngle);
                //WeatherVariables_list stores hourAngle
                WeatherVariables_list[i].HourAngleRad = hourAngle;
                //cosZ (radians) cosine Zenith angle
                cosZ = Math.Sin(solDecl) * Math.Sin(Latitude_rad) + Math.Cos(solDecl) * Math.Cos(Latitude_rad) * Math.Cos(hourAngle);
                if (cosZ > 1)
                {
                    cosZ = 1.0;
                }
                else if (cosZ < -1)
                {
                    cosZ = -1.0;
                }

                //calculate atmosphere correction for the solar zenith angle
                exAtmZenAgl_deg = Rad2Deg(Math.Acos(cosZ));    //in degrees
                exAtmElvAgl_deg = 90.0 - exAtmZenAgl_deg;           //in degrees

                refCorr_deg = 0.0;
                if (-1 <= exAtmElvAgl_deg && exAtmElvAgl_deg < 15)
                {
                    refCorr_deg = ABS_PRS / ABS_TEMP
                            * (0.1594 + 0.0196 * exAtmElvAgl_deg + 0.00002 * Math.Pow(exAtmElvAgl_deg, 2))
                            / (1 + 0.505 * exAtmElvAgl_deg + 0.0845 * Math.Pow(exAtmElvAgl_deg, 2));
                }
                else if (15 <= exAtmElvAgl_deg && exAtmElvAgl_deg < 90)
                {
                    refCorr_deg = 0.00452 * ABS_PRS / ABS_TEMP / Math.Tan(Deg2Rad(exAtmElvAgl_deg));
                }
                WeatherVariables_list[i].SolZenAgl_deg = exAtmZenAgl_deg - refCorr_deg;
            }
        }
        /// <summary>
        /// Convert unit of angle from degrees to radians.
        /// </summary>
        /// <param name="angle">Angle in degrees</param>
        /// <returns>Angle in radians</returns>
        static double Deg2Rad(double angle)
        {
            return Math.PI * angle / 180.0;
        }
        /// <summary>
        /// Convert unit of angle from radians to degrees.
        /// </summary>
        /// <param name="angle">Angle in radians</param>
        /// <returns>Angle in degrees</returns>
        static double Rad2Deg(double angle)
        {
            return 180.0 * angle / Math.PI;
        }
        /// <summary>
        /// Calculate solar radiation, aka shortwave radiation.
        /// </summary>
        /// <param name="WeatherVariables_list">List holding hourly surface weather data</param>
        /// <param name="LocationData_Object">Location data</param>
        /// <param name="Variable_TimeSeries_Hr_Length">Record count in WeatherVariables_list</param>
        public static void CalcSolarRadiation(List<SurfaceWeather> WeatherVariables_list, LocationData LocationData_Object, int Variable_TimeSeries_Hr_Length)
        {
            double ETR;       //ETR:   extraterrestrial radiation
            double airMass_relativeOpticalPressueCorrected_kg_p_m2; //M:     pressure-corrected relative optical air mass
            double CloudCover_Opaque_DirectNormalTransmittance_tenths;      //OPQN:  adjusted opaque cloud cover for diffuse radiation calculation
            double CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths;      //OPQD:  adjusted opaque cloud cover for direct radiation calculation
            double Kn;        //Kn:    transmission value of extraterrestrial direct normal radiation
            double Tr;        //TR:    Rayleigh scattering transmittance
            double Toz;       //TO:    ozone absorption transmittance
            double Tum;       //TUM:   uniformly mixed gas transmittance
            double Ta;        //TA:    aerosol absorption and scattering transmittance
            double Kd;
            int i;

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                CalcExtraterrestrialRadiation(WeatherVariables_list[i].Jday, WeatherVariables_list[i].SolZenAgl_deg, out ETR);

                CalcAirMass(WeatherVariables_list[i].SolZenAgl_deg, ref WeatherVariables_list[i].airMass_relativeOptical_kg_p_m2, WeatherVariables_list[i].AtmP_kPa, out airMass_relativeOpticalPressueCorrected_kg_p_m2);

                CalcCloudCover(LocationData_Object.Latitude, WeatherVariables_list[i].airMass_relativeOptical_kg_p_m2, WeatherVariables_list[i].CloudCover_Total_Tenth, WeatherVariables_list[i].Flag_SkyData_hasMultipleLayers, WeatherVariables_list[i].CloudCover_1_Tenth, WeatherVariables_list[i].CloudBaseHeight_1_m, WeatherVariables_list[i].CloudCover_2_Tenth, WeatherVariables_list[i].CloudBaseHeight_2_m, WeatherVariables_list[i].CloudCover_3_Tenth, WeatherVariables_list[i].CloudBaseHeight_3_m, ref WeatherVariables_list[i].CloudCover_Opaque_tenths, ref WeatherVariables_list[i].CloudCover_Translucent_tenths, out CloudCover_Opaque_DirectNormalTransmittance_tenths, out CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths);

                //CalcAtmosphericTransmission_Kn calculates Eq 15 to obtain Kn, direct normal insolation reaching Earth
                CalcAtmosphericTransmission_Kn(WeatherVariables_list[i].Day, WeatherVariables_list[i].airMass_relativeOptical_kg_p_m2, WeatherVariables_list[i].Tair_2m_C, WeatherVariables_list[i].Tdew_2m_C, WeatherVariables_list[i].AtmP_mbar, ref WeatherVariables_list[i].VapPres_Sat_kPa, ref WeatherVariables_list[i].VapPres_Act_kPa, ref WeatherVariables_list[i].RelHum, LocationData_Object.ozone[WeatherVariables_list[i].Month - 1], LocationData_Object.CoeffA, LocationData_Object.CoeffPHI, LocationData_Object.CoeffC, LocationData_Object.Elevation, airMass_relativeOpticalPressueCorrected_kg_p_m2, WeatherVariables_list[i].CloudCover_Opaque_tenths, WeatherVariables_list[i].CloudCover_Translucent_tenths, CloudCover_Opaque_DirectNormalTransmittance_tenths, out Tr, out Toz, out Tum, out Ta, out Kn);

                //CalcAtmosphericTransmission_Kd calculates Eq 32 to obtain Kd, for diffuse horizontal insolation reaching Earth
                CalcAtmosphericTransmission_Kd(WeatherVariables_list[i].airMass_relativeOptical_kg_p_m2, WeatherVariables_list[i].CloudCover_Opaque_tenths, WeatherVariables_list[i].CloudCover_Translucent_tenths, WeatherVariables_list[i].Ppt_1hr_m, LocationData_Object.AlbedoVal[WeatherVariables_list[i].Month - 1], Tr, Toz, Tum, Ta, CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, Kn, out Kd);

                //Calculate Earth's Surface Direct Shortwave Radiation with Eq 14
                //Radiation_Shortwave_Direct_Wpm2 = Kn * ETR
                //Note: Kn is Atmospheric Transmission Direct Normal, ETR is extraterreistrial insolation
                //Note: Radiation_Shortwave_Direct_Wpm2 from Eq 14 of Hirabayashi and Endreny (2025)
                WeatherVariables_list[i].Radiation_Shortwave_Direct_Wpm2 = Kn * ETR;
                //Calculate Earth's Surface Diffuse Shortwave Radiation
                //Radiation_Shortwave_Diffuse_Wpm2 = Kd * ETR
                //Note: Kd is Atmospheric Transmission Diffuse Horizontal, ETR is extraterreistrial insolation
                //Note: Radiation_Shortwave_Diffuse_Wpm2 from Eq 31 of Hirabayashi and Endreny (2025)
                WeatherVariables_list[i].Radiation_Shortwave_Diffuse_Wpm2 = Kd * ETR;
                //Calculate Earth's Surface Total Shortwave Radiation
                //Radiation_Shortwave_Total_Wpm2 = (Kn + Kd) * ETR
                //Note: Radiation_Shortwave_Total_Wpm2 from combined Eq 14 and 31 of Hirabayashi and Endreny (2025)
                WeatherVariables_list[i].Radiation_Shortwave_Total_Wpm2 = (Kn + Kd) * ETR;
                WeatherVariables_list[i].PARWm2 = WeatherVariables_list[i].Radiation_Shortwave_Total_Wpm2 * 0.46;
                WeatherVariables_list[i].PARuEm2s = WeatherVariables_list[i].PARWm2 * 4.6;

                if (WeatherData_Format != WeatherDataFormat.NARCCAP) {

                    double hourOfDay = WeatherVariables_list[i].TimeStamp.Hour + WeatherVariables_list[i].TimeStamp.Minute / 60.0;

                    //CalcNetRadWm2 calculates Radiation_Net_Wpm2
                    CalcNetRadWm2(WeatherVariables_list[i].Tair_2m_K, WeatherVariables_list[i].Tdew_2m_C, WeatherVariables_list[i].CloudCover_Total_Tenth, WeatherVariables_list[i].Flag_SkyData_hasMultipleLayers, WeatherVariables_list[i].CloudCover_1_Tenth, WeatherVariables_list[i].CloudBaseHeight_1_m, WeatherVariables_list[i].CloudCover_2_Tenth, WeatherVariables_list[i].CloudBaseHeight_2_m, WeatherVariables_list[i].CloudCover_3_Tenth, WeatherVariables_list[i].CloudBaseHeight_3_m, WeatherVariables_list[i].Radiation_Shortwave_Total_Wpm2, WeatherVariables_list[i].AtmP_mbar, hourOfDay, ref WeatherVariables_list[i].Radiation_Longwave_Downwelling_Wpm2, ref WeatherVariables_list[i].Radiation_Longwave_Upwelling_Wpm2, ref WeatherVariables_list[i].Radiation_Net_Wpm2, LocationData_Object.AlbedoVal[WeatherVariables_list[i].Month - 1]);
                }
            }
        }
        /// <summary>
        /// Calculate 
        /// </summary>
        /// <param name="jDay"></param>
        /// <param name="zenAgl_deg"></param>
        /// <param name="ETR"></param>
        /// <param name="ETRN"></param>
        static void CalcExtraterrestrialRadiation(int jDay, double zenAgl_deg, out double ETR)
        {
            double dayAngle;  //psi:   day angle
            double eccent;    //eo:    eccentricity of the earth's orbit

            dayAngle = Deg2Rad(360 * (jDay - 1) / 365.0);                        //in radians
            eccent = 1.00011 + 0.034221 * Math.Cos(dayAngle) + 0.00128 * Math.Sin(dayAngle)
                    + 0.000719 * Math.Cos(2 * dayAngle) + 0.000077 * Math.Sin(2 * dayAngle);

            //ETR is extraterrestrial radiation on a curved earth, adjusted for zenAgl_deg, from Eq 11 and Eq 12 combined
            ETR = eccent * SOL_CONST * Math.Cos(Deg2Rad(zenAgl_deg));

            //If ETR < 0 then set to 0
            if (ETR < 0)
            {
                ETR = 0.0;
            }
        }
        static void CalcAirMass(double zenAgl_deg, ref double airMass_relativeOptical_kg_p_m2, double prsKpa, out double airMass_relativeOpticalPressueCorrected_kg_p_m2)
        {
            double solElvAgl_deg; //CloudBase_normalized_frac: solar elevation angle
            double prsMbar = prsKpa * 10.0;

            solElvAgl_deg = 90.0 - zenAgl_deg;   //in degrees
            if (zenAgl_deg <= 90)
            {
                //airMass_relativeOptical_kg_p_m2 is 1 / (Math.Sin(Deg2Rad(solElvAgl_deg)) + 0.50572 * Math.Pow((solElvAgl_deg + 6.07995), -1.6364))
                //Note: airMass_relativeOptical_kg_p_m2 from Eq 18 of Hirabayashi and Endreny (2025) and Eq 3-1 of NREL (1995)
                //Note: airMass_relativeOptical_kg_p_m2 should range from 1 at solar noon, to about 30 or 40 at sunrise and sunset
                airMass_relativeOptical_kg_p_m2 = 1 / (Math.Sin(Deg2Rad(solElvAgl_deg)) + 0.50572 * Math.Pow((solElvAgl_deg + 6.07995), -1.6364));
            }
            else
            {
                airMass_relativeOptical_kg_p_m2 = 99.0;
            }
            if (airMass_relativeOptical_kg_p_m2 >= 30)
            {
                airMass_relativeOptical_kg_p_m2 = 30.0;
            }
            //airMass_relativeOpticalPressueCorrected_kg_p_m2 is airMass_relativeOptical_kg_p_m2 * (prsMbar) / 1013.0
            //Note: airMass_relativeOpticalPressueCorrected_kg_p_m2 from Eq 17 of Hirabayashi and Endreny (2025) and Table 3-4 of NREL (1995)
            airMass_relativeOpticalPressueCorrected_kg_p_m2 = airMass_relativeOptical_kg_p_m2 * (prsMbar) / 1013.0;
            //Warning: This seems to be an error, and prsKpa / 10 should be prsKpa * 10 to convert to prsMbar
            //Note: 2024 and prior error used: airMass_relativeOpticalPressueCorrected_kg_p_m2 = airMass_relativeOptical_kg_p_m2 * (prsKpa / 10.0) / 1013.0;
        }

        public static double HeightWeight(double CloudBase_i_m, double Latitude_rad)
        {
            double CloudBase_min_m = 30.0;
            //Cloud maximum height at equator is 18 to 20 km, at pole 6 to 9 km 
            const double CloudBase_equator_m = 18000.0;
            const double CloudBase_pole_m = 7000.0;
            double CloudBase_max_m = CloudBase_pole_m + (CloudBase_equator_m - CloudBase_pole_m) * Math.Cos(Latitude_rad);
            //CloudBase_max_m clamaped in case of any rounding weirdness, clamp to [CloudBase_pole_m, CloudBase_equator_m]
            if (CloudBase_max_m < CloudBase_pole_m) CloudBase_max_m = CloudBase_pole_m;
            if (CloudBase_max_m > CloudBase_equator_m) CloudBase_max_m = CloudBase_equator_m;

            //Opacity_Cloud_High_min_frac set to 0.1 (fraction), the minimum on scale of 0 to 1 (fraction)
            double Opacity_Cloud_High_min_frac = 0.1;

            //If CloudBase_i_m == NODATA then return
            if (CloudBase_i_m == NODATA) {
                return 0.0;
            }

            //If CloudBase_i_m <= 0.0 then a very low cloud
            if (CloudBase_i_m <= 0.0) {
                CloudBase_i_m = CloudBase_min_m;
            }

            //Clamp CloudBase_i_m between min and max
            double z_clamped = Math.Max(CloudBase_min_m, Math.Min(CloudBase_max_m, CloudBase_i_m));

            //Take log of CloudBase_i_m to evenly scale from 30 to 18000 
            double log_CloudBase_min = Math.Log(CloudBase_min_m);
            double log_CloudBase_max = Math.Log(CloudBase_max_m);
            double log_CloudBase = Math.Log(z_clamped);
            //CloudBase_normalized_frac = (log_CloudBase - log_CloudBase_min) / (log_CloudBase_max - log_CloudBase_min)
            double CloudBase_normalized_frac = (log_CloudBase - log_CloudBase_min) / (log_CloudBase_max - log_CloudBase_min);

            //HeightWeight_frac = Opacity_Cloud_High_min_frac + (1.0 - Opacity_Cloud_High_min_frac) * (1.0 - CloudBase_normalized_frac)
            //Note: 1.0 at low clouds, Opacity_Cloud_High_min_frac at high clouds
            double HeightWeight_frac = Opacity_Cloud_High_min_frac + (1.0 - Opacity_Cloud_High_min_frac) * (1.0 - CloudBase_normalized_frac);

            return Math.Max(Opacity_Cloud_High_min_frac, Math.Min(1.0, HeightWeight_frac));
        }


        //CalcCloudCover function computes CloudCover_Opaque_DirectNormalTransmittance_tenths and CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths
        static void CalcCloudCover(
            double Latitude_deg,
            double airMass_relativeOptical_kg_p_m2,
            double CloudCover_Total_Tenth,
            bool Flag_SkyData_hasMultipleLayers,
            double CloudCover_1_Tenth,
            double CloudBaseHeight_1_m,
            double CloudCover_2_Tenth,
            double CloudBaseHeight_2_m,
            double CloudCover_3_Tenth,
            double CloudBaseHeight_3_m,
            ref double CloudCover_Opaque_tenths,
            ref double CloudCover_Translucent_tenths,
            out double CloudCover_Opaque_DirectNormalTransmittance_tenths,
            out double CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths
            )

        {
            double a1;
            double b1;
            double corrN;
            double corrD;
            //Latitude_rad = Deg2Rad(Math.Abs(Latitude_deg)) is absolute latitude (south was negative degrees)
            double Latitude_rad = Deg2Rad(Math.Abs(Latitude_deg));

            //Algorithm is loosely based on Atkinson et al. (1992) and NREL (1995) NSRDB partitioning total cloud cover to opaque and translucent cover
            //Reference:
            //Atkinson, D. and R. F. Lee (1992). "Procedures for Substituting Values for Missing NWS Meteorological Data for Use in Regulatory Air Quality Models." Environmental Protection Agency (EPA).
            //USEPA(1997).Analysis of the Affect of ASOS - Derived Meteorological Data on Refined Modeling.Research Triangle Park, NC, US Environmental Protection Agency.
            /* 
            //Original SAS Code to assign opaque cloud cover based on ceiling height (CLG) and total cloud cover (SKC)
            //Presumed OPQCC is missing and OPQCC_T is CLG
            // IF OPQCC has a value then total cloud cover is OPQCC
            IF OPQCC GE .Z THEN OPQCC_T = OPQCC;
            // If OPQCC is missing and TOTCC is not missing
            IF OPQCC = . AND TOTCC NE . THEN OPQCC = TOTCC;
            // Else if OPQCC is missing and TOTCC is missing ...
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING > 70 THEN OPQCC = 0;
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING LE 70 AND CEILING GT 60 THEN OPQCC = 7;
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING LE 60 AND CEILING GT 50 THEN OPQCC = 8;
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING LE 50 AND CEILING GT 30 THEN OPQCC = 9;
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING LE 30 THEN OPQCC = 10;
            ELSE IF OPQCC = . AND TOTCC =  . AND CEILING = .   THEN OPQCC = OPQCC_T;
            IF TOTCC  = . AND OPQCC NE . THEN TOTCC = OPQCC;
            */

            // Cloud cover in fraction (0–1)
            double CloudCover_Total_frac = (CloudCover_Total_Tenth != NODATA)
                ? CloudCover_Total_Tenth / 10.0
                : 0.0;

            //UnlimitedCloudBase_m defined at 18000 m, below the 22000 m ~ 722 hundredths ft, ISD code for unlimited sky 
            const double UnlimitedCloudBase_m = 18000.0;

            bool noCloudsOrUnlimited = CloudCover_Total_Tenth == NODATA || CloudCover_Total_frac <= 0.0 || CloudBaseHeight_1_m == NODATA || CloudBaseHeight_1_m >= UnlimitedCloudBase_m;

            if (noCloudsOrUnlimited) {
                CloudCover_Opaque_tenths = 0.0;
            }

            //Else if Flag_SkyData_hasMultipleLayers false then a single layer is processed
            else if (!Flag_SkyData_hasMultipleLayers)
            {
                // ---------- SINGLE-LAYER CASE (ISD or 1-layer GHCNh) ----------

                double CloudCover_frac = CloudCover_Total_Tenth / 10.0;
                if (CloudCover_frac <= 0.0) {
                    CloudCover_Opaque_tenths = 0.0;
                }
                else {
                    double CloudBase_min_m = 30.0;

                    //If ceiling missing or <=0, treat as low cloud at CloudBase_min_m
                    //Note: GHCNh and ISD parameterize CloudBaseHeight_1_m
                    double CloudBase_i_m;
                    if (CloudBaseHeight_1_m == NODATA || CloudBaseHeight_1_m <= 0.0) {
                        CloudBase_i_m = CloudBase_min_m;
                    }
                    else {
                        CloudBase_i_m = CloudBaseHeight_1_m;
                    }

                    double w = HeightWeight(CloudBase_i_m, Latitude_rad);
                    double Cloud_Opaque_frac = CloudCover_frac * w;

                    Cloud_Opaque_frac = Math.Max(0.0, Math.Min(1.0, Cloud_Opaque_frac));
                    CloudCover_Opaque_tenths = Cloud_Opaque_frac * 10.0;
                }
            }


            //Else if Flag_SkyData_hasMultipleLayers is true then multiple cloud heights
            else
            {
                // ---------- MULTI-LAYER CASE (GHCNh with 2–3 layers) ----------

                // Cumulative cover in tenths per layer (could be NODATA)
                CloudCover_1_Tenth = (CloudCover_1_Tenth != NODATA) ? CloudCover_1_Tenth : 0.0;
                CloudCover_2_Tenth = (CloudCover_2_Tenth != NODATA) ? CloudCover_2_Tenth : 0.0;
                CloudCover_3_Tenth = (CloudCover_3_Tenth != NODATA) ? CloudCover_3_Tenth : 0.0;

                // Convert to cumulative fractions
                double CloudCover_1_frac = CloudCover_1_Tenth / 10.0;
                double CloudCover_2_frac = CloudCover_2_Tenth / 10.0;
                double CloudCover_3_frac = CloudCover_3_Tenth / 10.0;

                // Incremental fractions by layer:
                // layer 1: c1
                // layer 2: max(c2 - c1, 0)
                // layer 3: max(c3 - c2, 0)
                CloudCover_1_frac = Math.Max(0.0, CloudCover_1_frac);
                double CloudCover_2_1_frac = Math.Max(0.0, CloudCover_2_frac - CloudCover_1_frac);
                double CloudCover_3_2_frac = Math.Max(0.0, CloudCover_3_frac - CloudCover_2_frac);

                // Use reported base heights for each layer
                double CloudBase_1_m = CloudBaseHeight_1_m;
                double CloudBase_2_m = CloudBaseHeight_2_m;
                double CloudBase_3_m = CloudBaseHeight_3_m;

                double CloudBase_weight_1 = HeightWeight(CloudBase_1_m, Latitude_rad);
                double CloudBase_weight_2 = HeightWeight(CloudBase_2_m, Latitude_rad);
                double CloudBase_weight_3 = HeightWeight(CloudBase_3_m, Latitude_rad);

                // Sum height-weighted opaque fractions
                double Cloud_Opaque_frac =
                    CloudCover_1_frac * CloudBase_weight_1 +
                    CloudCover_2_1_frac * CloudBase_weight_2 +
                    CloudCover_3_2_frac * CloudBase_weight_3;

                // Cap at total cover fraction and 1.0
                Cloud_Opaque_frac = Math.Max(0.0,
                                      Math.Min(CloudCover_Total_frac,
                                      Math.Min(1.0, Cloud_Opaque_frac)));

                CloudCover_Opaque_tenths = Cloud_Opaque_frac * 10.0;
            }

            // Translucent is the remainder of total
            CloudCover_Translucent_tenths = CloudCover_Total_Tenth - CloudCover_Opaque_tenths;
            if (CloudCover_Translucent_tenths < 0.0)
                CloudCover_Translucent_tenths = 0.0;


            // Translucent cloud is the remainder of total cover
            // (Eq 44 Hirabayashi & Endreny, Eq 3-24 NREL (1995))
            CloudCover_Translucent_tenths = CloudCover_Total_Tenth - CloudCover_Opaque_tenths;

            // Ensure non-negative
            if (CloudCover_Translucent_tenths < 0.0)
            {
                CloudCover_Translucent_tenths = 0.0;
            }

            /*
            // ----------------------------------------------------------
            // Optional debug print for FIRST DAY of simulation only
            // ----------------------------------------------------------
            try
            {
                Console.WriteLine(
                    
                    $"Height={CeilingHeight_hundredths_ft:F1} hundredths ft, " +
                    $"Total={CloudCover_Total_Tenth:F1} tenths, " +
                    $"Opaque={CloudCover_Opaque_tenths:F1} tenths, " +
                    $"Translucent={CloudCover_Translucent_tenths:F1} tenths");
            }
            catch { }
            */

            // -------- end new opaque/translucent block --------

            //OPQD is Effective opaque cloud cover for diffuse horizontal calculations (tenths)
            //OPQN is Effective opaque cloud cover for direct normal transmittance (tenths)

            //If (1 <= CloudCover_Opaque_tenths && CloudCover_Opaque_tenths <= 9)
            //Note: Criteria that OPQ is 1 to 9 explained in NREL (1995) under Eq 3-18 ...
            //...Furthermore, under clear (OPQ = 0) or overcast (OPQ = 10) skies, no adjustments were made to the observed values, i.e., both N and D were equal to zero
            //Note: From NREL (1995) Algorithm around Eq 3-16
            if (1 <= CloudCover_Opaque_tenths && CloudCover_Opaque_tenths <= 9) {
                //a1 = 4.955 * (1 - Math.Exp(-0.454 * airMass_relativeOptical_kg_p_m2)) - 3.4
                //Note: a1 from Eq 28 of Hirabayashi and Endreny (2025) and Eq 3-16 of NREL (1995)
                a1 = 4.955 * (1 - Math.Exp(-0.454 * airMass_relativeOptical_kg_p_m2)) - 3.4;
                //if (a1 <= 0)
                if (a1 <= 0) {
                    //b1 = -0.2 * a1
                    //Note: b1 from conditional Eq 29 of Hirabayashi and Endreny (2025) and Eq 3-17 of NREL (1995)
                    b1 = -0.2 * a1;
                }
                else {
                    //Note: b1 from conditional Eq 29 of Hirabayashi and Endreny (2025) and Eq 3-17 of NREL (1995)
                    b1 = 0.1 * a1;
                }
                //corrN = a1 * Math.Sin(Math.PI / 10 * CloudCover_Total_Tenth) + b1 * Math.Sin(2 * Math.PI / 10 * CloudCover_Total_Tenth)
                //Note: corrN from Eq 27 of Hirabayashi and Endreny (2025) and Eq 3-15 of NREL (1995)
                //Note: Original equation used sin(18 * OPQ) and sin(36 * OPQ) which is equivalent of sin(180 * OPQ/10) and sin(360 * OPQ/10) 
                //Note: ... where OPQ units tenths (1 to 10) and sin was operating on degrees not radians
                corrN = a1 * Math.Sin(Math.PI / 10 * CloudCover_Total_Tenth) + b1 * Math.Sin(2 * Math.PI / 10 * CloudCover_Total_Tenth);
                //corrD = a1 * Math.Sin(Math.PI / 10 * CloudCover_Total_Tenth) * 0.5
                //Note: corrD from Eq 42 of Hirabayashi and Endreny (2025) and Eq 3-18 of NREL (1995)
                //Note: Original equation used sin(18 * OPQ) which is equivalent of sin(180 * OPQ/10) 
                //Note: ... where OPQ units tenths (1 to 10) and sin was operating on degrees not radians
                corrD = a1 * Math.Sin(Math.PI / 10 * CloudCover_Total_Tenth) * 0.5;
            }
            //else Algorithm from below Eq 3-18 in NREL (1995)
            //Note: Furthermore, under clear (OPQ = 0) or overcast (OPQ = 10) skies, no adjustments were made to the observed values, i.e., both N and D were equal to zero
            else
            {
                //corrN = 0.0 when OPQ = 0 or 10, algorithm from below Eq 3-18 in NREL (1995)
                corrN = 0.0;
                //corrD = 0.0 when OPQ = 0 or 10, algorithm from below Eq 3-18 in NREL (1995)
                corrD = 0.0;
            }
            //CloudCover_Opaque_DirectNormalTransmittance_tenths = CloudCover_Opaque_tenths + corrN;
            //OPQN is Effective Opaque Cloud Cover for Direct Normal Transmittance (tenths)
            //Note: CloudCover_Opaque_DirectNormalTransmittance_tenths within Eq 26 of Hirabayashi and Endreny (2025) and Eq 3-20 of NREL (1995)
            CloudCover_Opaque_DirectNormalTransmittance_tenths = CloudCover_Opaque_tenths + corrN;

            //CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths = CloudCover_Opaque_tenths + corrD
            //OPQD is Effective opaque cloud cover for diffuse horizontal calculations (tenths)
            //Note: CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths from Eq 41 of Hirabayashi and Endreny (2025) and below Eq 3-23 of NREL (1995)
            CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths = CloudCover_Opaque_tenths + corrD;
        }

        //CalcAtmosphericTransmission_Kn function computes direct normal transmittance
        //Note: Kn = direct beam transmittance, from Eq 2-1 of NREL (1995), Kn = In / I0
        //Note: Kt = clearness index or effective global horizontal transmittance, from Eq 2-2 of NREL (1995), Kt = It / (I0 cos Z)
        //Note: Kd = effective diffuse horizontal transmittance, from Eq 2-3 of NREL (1995),  Kd = Id / (I0 cos Z)
        //Note: Kt from Eq 2-4 of NREL (1995),  Kt = Kn + Kd
        //Note: I0 = extraterrestrial direct normal irradiance
        //Note: In = direct normal irradiance at the Earth's surface
        //Note: 1t = total global horizontal irradiance at the Earth's surface
        //Note: Id = diffuse horizontal irradiance at the Earth's surface
        //Note: I0 cos z = extraterrestrial irradiance on a surface parallel to the Earth's surface
        //Note: Z = solar zenith angle
        static void CalcAtmosphericTransmission_Kn(int day, double airMass_relativeOptical_kg_p_m2, double Tair_2m_C, double Tdew_2m_C, double AtmPres_Station_Mb, ref double VapPres_Sat_kPa, ref double VapPres_Act_kPa, ref double RelHum, double ozone, double a, double phi, double C, double ElevationStation_m, double airMass_relativeOpticalPressueCorrected_kg_p_m2, double CloudCover_Opaque_tenths, double CloudCover_Translucent_tenths, double CloudCover_Opaque_DirectNormalTransmittance_tenths, out double Tr, out double Toz, out double Tum, out double Ta, out double Kn)
        {
            double Tw;        //Tw:    water vapor absorption transmittance
            double Ttrn;      //TTRN:  translucent cloud transmittance
            double Topq;      //TOPQ:  opaque cloud cover transmittance
            double Xo;        //XO:    total amount of ozone in cm in a slanted path
            double tau;       //tauA:  broadband aerosol optical depth
            double vPrsMb;    //vapor pressure in mBar
            double vPrs_StationPressureCorrected_Mb; //pressure adjusted vapor pressure in mBar
            double PrecipitableWater_mm;        //PrecipitableWater_mm:    total precipitable water
            double Xw;        //Xw:    precipitable water vapor in cm in a slant path
            double ATRN_directNormalTransmittanceIntercept = 1.0;
            double BTRN_directNormalTransmittanceSlope = 0.0;
            double M_PrecipitableWater_slope;
            double B_PrecipitableWater_intercept;

            //Tr is Raleigh scattering transmittance
            //Tr = Math.Exp((-0.0903 * Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 0.84)) * (1 + airMass_relativeOpticalPressueCorrected_kg_p_m2 - Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 1.01)))
            //Note: Tr from Eq 16 of Hirabayashi and Endreny (2025) and Eq 3-4 of NREL (1995) 
            Tr = Math.Exp((-0.0903 * Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 0.84)) * (1 + airMass_relativeOpticalPressueCorrected_kg_p_m2 - Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 1.01)));

            //Xo = ozone * airMass_relativeOptical_kg_p_m2
            //Note: Xo from Eq 20 of Hirabayashi and Endreny (2025) and Table (of symbols) 3-4 of NREL (1995) 
            Xo = ozone * airMass_relativeOptical_kg_p_m2;

            //Toz is ozone absorption transmittance
            //Toz = 1 - 0.1611 * Xo * Math.Pow((1 + 139.45 * Xo), -0.3035) - 0.002715 * Xo / (1 + 0.044 * Xo + 0.0003 * Math.Pow(Xo, 2))
            //Note: Toz from Eq 19 of Hirabayashi and Endreny (2025) and Eq 3-5 of NREL (1995) 
            Toz = 1 - 0.1611 * Xo * Math.Pow((1 + 139.45 * Xo), -0.3035) - 0.002715 * Xo / (1 + 0.044 * Xo + 0.0003 * Math.Pow(Xo, 2));

            //Tum is uniformly mixed gas absorption transmittance
            //Tum = Math.Exp(-0.0127 * Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 0.26))
            //Note: Tum from Eq 21 of Hirabayashi and Endreny (2025) and Eq 3-6 of NREL (1995) 
            Tum = Math.Exp(-0.0127 * Math.Pow(airMass_relativeOpticalPressueCorrected_kg_p_m2, 0.26));

            //tau = a * Math.Sin(360 * (double)day / 365 - phi) + C
            //Note: tau from Eq 25 of Hirabayashi and Endreny (2025) and Eq 6-5 of NREL (1995) and Eq 4 of Maxwell (1998)
            //Note: Algorithm and Coefficients a, phi, C from NREL (1995) Section 6.0, Daily Estimates of Aerosol Optical Depths
            tau = a * Math.Sin(360 * (double)day / 365 - phi) + C;

            //Ta is aerosol absorption and scattering transmittance
            //Ta = Math.Exp(-tau * airMass_relativeOptical_kg_p_m2)
            //Note: Ta from Eq 24 of Hirabayashi and Endreny (2025) and Eq 3-8 of NREL (1995)
            Ta = Math.Exp(-tau * airMass_relativeOptical_kg_p_m2);

            //VapPres_Sat_kPa = ES_CONS * Math.Exp((17.27 * Tair_2m_C) / (237.3 + Tair_2m_C))
            //Note: VapPres_Sat_kPa from Eq 66 of Hirabayashi and Endreny (2025)
            VapPres_Sat_kPa = ES_CONS * Math.Exp((17.27 * Tair_2m_C) / (237.3 + Tair_2m_C));
            //VapPres_Act_kPa = ES_CONS * Math.Exp((17.27 * Tdew_2m_C) / (237.3 + Tdew_2m_C))
            //Note: VapPres_Act_kPa from Eq 67 of Hirabayashi and Endreny (2025)
            VapPres_Act_kPa = ES_CONS * Math.Exp((17.27 * Tdew_2m_C) / (237.3 + Tdew_2m_C));
            //RelHum = VapPres_Act_kPa / VapPres_Sat_kPa
            //Note: Standard relation
            RelHum = VapPres_Act_kPa / VapPres_Sat_kPa;
            //vPrsMb = VapPres_Act_kPa * 10.0
            //Note: Standard conversion
            vPrsMb = VapPres_Act_kPa * 10.0;

            //vPrs_StationPressureCorrected_Mb = vPrsMb * AtmPres_Station_Mb / 1013.25
            //Note: vPrs_StationPressureCorrected_Mb from Eq (missing) of Hirabayashi and Endreny (2025) and Eq 5-8 of NREL (1995)
            //Note: Algorithm and Coefficients M and B from NREL (1995) Section 5.0, Hourly Estimates of Precipitable Water
            vPrs_StationPressureCorrected_Mb = vPrsMb * AtmPres_Station_Mb / 1013.25;

            //M_PrecipitableWater_slope = (0.0004 * ElevationStation_m + 1.1)
            //Note: M_PrecipitableWater_slope from Eq (missing) of Hirabayashi and Endreny (2025) and Eq 5-10 of NREL (1995) 
            M_PrecipitableWater_slope = (0.0004 * ElevationStation_m + 1.1);
            //B_PrecipitableWater_intercept is 1 
            //Warning: B_PrecipitableWater_intercept should be 1 or 2 based on humidity and Table 5-2 of NREL (1995)
            B_PrecipitableWater_intercept = 1;

            //PrecipitableWater_mm = vPrs_StationPressureCorrected_Mb * M_PrecipitableWater_slope + B_PrecipitableWater_intercept
            //Note: vPrs_StationPressureCorrected_Mb from Eq (missing) of Hirabayashi and Endreny (2025) and Eq 5-9 of NREL (1995) 
            //Warning: Coefficient B set to 1, but can be 1 or 2 as shown in Table 5-2 of NREL (1995)
            PrecipitableWater_mm = vPrs_StationPressureCorrected_Mb * M_PrecipitableWater_slope + B_PrecipitableWater_intercept;

            //PrecipitableWater_mm = PrecipitableWater_mm / 10.0
            //Note: PrecipitableWater_cm should have max value = 6 cm, stated below Eq 3-7 and as known in meteorology
            //Note: Standard conversion
            double PrecipitableWater_cm = PrecipitableWater_mm / 10.0;

            //Note: Xw = PrecipitableWater_cm from Eq (missing) of Hirabayashi and Endreny (2025) and Eq 5-9 and Eq 5-10 of NREL (1995)
            //Warning: PrecipitableWater_mm was erroneously multipled by airMass_relativeOptical_kg_p_m2; confused m air mass and slope symbols 
            Xw = PrecipitableWater_cm;

            //Tw is water vapor transmittance
            //Tw = 1 - 1.668 * Xw / (Math.Pow((1 + 54.6 * Xw), 0.637) + 4.042 * Xw)
            //Note: Tw from Eq 22 of Hirabayashi and Endreny (2025) and Eq 3-7 of NREL (1995) and Eq 3 of Maxwell (1998)
            //Note: Xw (cm) is the precipitable water vapor in a slant path through the atmosphere
            Tw = 1 - 1.668 * Xw / (Math.Pow((1 + 54.6 * Xw), 0.637) + 4.042 * Xw);

            //Translucent transmittance
            //Ttrn = 1, initialized to unlimitted; developed from 3.4.3 Translucent Cloud Algorithms of NREL (1995)
            //Note: Ttrn from Eq (missing) of Hirabayashi and Endreny (2025) and Eq 3-25 of NREL (1995)
            Ttrn = 1;
            //If CloudCover_Translucent_tenths >= 0.5 then enter, which is an integer value of 1
            //Note: The algorithm is only active for CloudCover_Translucent_tenths integer values of 1 or higher
            if (CloudCover_Translucent_tenths > 0.5) {
                //BTRN_0 and Coeff_c initialized
                double BTRN_0 = 0.0;
                double Coeff_c = 0.0;
                //CloudCover_Translucent_tenths_int is CloudCover_Translucent_tenths rounded to nearest integer
                int CloudCover_Translucent_tenths_int = (int)Math.Round(CloudCover_Translucent_tenths);
                //GetTranslucentCloudCoefficients takes CloudCover_Translucent_tenths_int and returns BTRN_0 and Coeff_c
                //Note: GetTranslucentCloudCoefficients accesses values from Table 3-6 in NREL (1995)
                GetTranslucentCloudCoefficients(CloudCover_Translucent_tenths_int, out BTRN_0, out Coeff_c);

                //BTRN_directNormalTransmittanceSlope = BTRN_0 + Coeff_c * CloudCover_Opaque_tenths
                //Note: BTRN_directNormalTransmittanceSlope from Eq 3-28 in NREL (1995)
                BTRN_directNormalTransmittanceSlope = BTRN_0 + Coeff_c * CloudCover_Opaque_tenths;

                //if (CloudCover_Opaque_tenths + CloudCover_Translucent_tenths > 10)
                //Note: This condition should not happen, given 10 is the maximum sky cloud cover
                //Note: Conditioanl algorithm developed with DeepSeek to satisfy all values of OPQ and avoid using Table 3-7
                //Note: Eq 3-27 of NREL (1995) does not fit the data values in Table 3-7; it is valid only when m = 1.0 and OPQ = 0.
                if (CloudCover_Opaque_tenths + CloudCover_Translucent_tenths > 10) {
                    ATRN_directNormalTransmittanceIntercept = 0.0;
                }
                //else if (CloudCover_Opaque_tenths == 0 || CloudCover_Translucent_tenths_int == 1)
                //Note: This condition is follows both the top row and left column of Table 3-7 in NREL (1995)
                else if (CloudCover_Opaque_tenths == 0 || CloudCover_Translucent_tenths_int == 1) {
                    //ATRN_directNormalTransmittanceIntercept = 0.995 * Math.Pow(0.9958, CloudCover_Translucent_tenths - 1) * Math.Pow(0.9975, CloudCover_Opaque_tenths)
                    //Note: Equation was developed in DeepSeek to improve on retrieving values from Table 3-25 in NREL (1995)
                    ATRN_directNormalTransmittanceIntercept = 0.995 * Math.Pow(0.9958, CloudCover_Translucent_tenths - 1) * Math.Pow(0.9975, CloudCover_Opaque_tenths);
                }
                //else
                //Note: This condition is approximates all cells but top row and left column of Table 3-7 in NREL (1995)
                else {
                    //ATRN_directNormalTransmittanceIntercept = 0.995 * Math.Pow(0.993, CloudCover_Translucent_tenths - 1) * Math.Pow(0.994, CloudCover_Opaque_tenths)
                    //Note: Equation was developed in DeepSeek to improve on retrieving values from Table 3-25 in NREL (1995)
                    ATRN_directNormalTransmittanceIntercept = 0.995 * Math.Pow(0.993, CloudCover_Translucent_tenths - 1) * Math.Pow(0.994, CloudCover_Opaque_tenths);
                }
                //Ttrn = ATRN_directNormalTransmittanceIntercept - BTRN_directNormalTransmittanceSlope * airMass_relativeOptical_kg_p_m2
                //Note: Ttrn from Eq 30 of Hirabayashi and Endreny (2025) and Eq 3-25 of NREL (1995)
                Ttrn = ATRN_directNormalTransmittanceIntercept - BTRN_directNormalTransmittanceSlope * airMass_relativeOptical_kg_p_m2;
                //If (Ttrn < 0 ) then set to 0
                if (Ttrn < 0 ) {
                    Ttrn = 0;
                }
            }

            //Opaque cloud
            //Topq = (10 - CloudCover_Opaque_DirectNormalTransmittance_tenths) / 10.0
            //Note: Topq from Eq 26 of Hirabayashi and Endreny (2025) and Eq 3-19 of NREL (1995)
            //Note: CloudCover_Opaque_DirectNormalTransmittance_tenths combines opq + n from 
            Topq = (10 - CloudCover_Opaque_DirectNormalTransmittance_tenths) / 10.0;

            //Kn is direct beam transmittance
            //Kn = 0.9751 * Tr * Toz * Tum * Tw * Ta * Topq * Ttrn
            //Note: Kn from Eq 15 of Hirabayashi and Endreny (2025) and Eq 5 of Maxwell (1998)
            Kn = 0.9751 * Tr * Toz * Tum * Tw * Ta * Topq * Ttrn;
              
        }

        //GetTranslucentCloudCoefficients accesses values from Table 3-6 in NREL (1995)
        static void GetTranslucentCloudCoefficients(int CloudCover_Translucent_tenths, out double BTRN_0, out double Coeff_c)
        {
            BTRN_0 = 0.0;
            Coeff_c = 0.0;

            //switch contains values from Table 3-6 in NREL (1995)
            switch (CloudCover_Translucent_tenths)
            {
                case 1: BTRN_0 = 0.004; Coeff_c = 0.003; break;
                case 2: BTRN_0 = 0.006; Coeff_c = 0.004; break;
                case 3: BTRN_0 = 0.009; Coeff_c = 0.005; break;
                case 4: BTRN_0 = 0.012; Coeff_c = 0.006; break;
                case 5: BTRN_0 = 0.016; Coeff_c = 0.006; break;
                case 6: BTRN_0 = 0.020; Coeff_c = 0.006; break;
                case 7: BTRN_0 = 0.026; Coeff_c = 0.006; break;
                case 8: BTRN_0 = 0.032; Coeff_c = 0.006; break;
                case 9: BTRN_0 = 0.040; Coeff_c = 0.006; break;
                case 10: BTRN_0 = 0.050; Coeff_c = 0.006; break;
                default:
                    Console.WriteLine("Invalid translucent cloud cover value.");
                    break;
            }
        }

        //CalcAtmosphericTransmission_Kd function computes effective diffuse horizontal transmittance 
        //Note: Kn = direct beam transmittance, from Eq 2-1 of NREL (1995), Kn = In / I0
        //Note: Kt = clearness index or effective global horizontal transmittance, from Eq 2-2 of NREL (1995), Kt = It / (I0 cos Z)
        //Note: Kd = effective diffuse horizontal transmittance, from Eq 2-3 of NREL (1995),  Kd = Id / (I0 cos Z)
        //Note: Kt from Eq 2-4 of NREL (1995),  Kt = Kn + Kd
        //Note: I0 = extraterrestrial direct normal irradiance
        //Note: In = direct normal irradiance at the Earth's surface
        //Note: 1t = total global horizontal irradiance at the Earth's surface
        //Note: Id = diffuse horizontal irradiance at the Earth's surface
        //Note: I0 cos z = extraterrestrial irradiance on a surface parallel to the Earth's surface
        //Note: Z = solar zenith angle
        static void CalcAtmosphericTransmission_Kd(double airMass_relativeOptical_kg_p_m2, double CloudCover_Opaque_tenths, double Trans, double rain, double Albedo_surface_frac, double Tr, double Toz, double Tum, double Ta, double CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, double Kn, out double Kd)
        {
            double Ksr;       //KSR:   diffuse radiation from Rayleigh scattering
            double Tas;       //TAS:   trensmittance of aerosol scattering
            double Taa;       //TAA:   transmittance of aerosol absorptance
            double Ksa;       //KSA:   diffuse radiation from aerosol scattering
            double fM;        //f(M):  empirical air mass function
            double Ksopq;     //KSOPQ:
            double Kstrn;     //KSTRN:
            double Ksgrf;     //KSGRF: diffuse radiation from ground reflectance
            double Ksgrf1;    //KSGRF1:
            double Ksgrf2;    //KSGRF2:

            double b2;
            double c2;
            double PSW_precipitableWaterSwitch;       //PSW:
            double Kd0;
            double Rcld;    //RCLD:  broadband (solar) cloud reflectance
            double Ratm;      //RATM:  broadband (solar) atmospheric reflectance

            //Taa = 1 - K1 * (1 - airMass_relativeOptical_kg_p_m2 + Math.Pow(airMass_relativeOptical_kg_p_m2, 1.06)) * (1 - Ta)
            //Note: Taa from Eq 36 of Hirabayashi and Endreny (2025) and Eq 3-9 of NREL (1995)
            Taa = 1 - K1 * (1 - airMass_relativeOptical_kg_p_m2 + Math.Pow(airMass_relativeOptical_kg_p_m2, 1.06)) * (1 - Ta);
            //Ksr = 0.5 * (1 - Tr) * Toz * Tum * Taa
            //Note: Ksr from Eq 35 of Hirabayashi and Endreny (2025) and Eq 3-11 of NREL (1995)
            Ksr = 0.5 * (1 - Tr) * Toz * Tum * Taa;
            //Ksa = BA * (1 - Ta) * Toz * Tum * Taa
            //Note: Ksa from Eq 37 of Hirabayashi and Endreny (2025) and Eq 3-12 of NREL (1995)
            Ksa = BA * (1 - Ta) * Toz * Tum * Taa;
            //fM = 0.38 + 0.925 * Math.Exp(-0.851 * airMass_relativeOptical_kg_p_m2)
            //Note: fM from Eq 34 of Hirabayashi and Endreny (2025) and Eq 3-14 of NREL (1995)
            fM = 0.38 + 0.925 * Math.Exp(-0.851 * airMass_relativeOptical_kg_p_m2);

            //    if CloudCover_Translucent_tenths > 0 {
            //        Kstrn = -0.00235 + 0.00689 * CloudCover_Translucent_tenths + 0.000203 * CloudCover_Translucent_tenths ^ 2
            //    }
            Kstrn = 0.0;

            Ksopq = 0.0;
            if (CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths >= 5)
            {
                //b2 = 0.0953 + 0.137 * CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths - 0.0409 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 2) + 0.00579 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 3) - 0.000328 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 4)
                //Note: b2 from Eq 39 of Hirabayashi and Endreny (2025) and Eq 3-22 of NREL (1995)
                b2 = 0.0953 + 0.137 * CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths - 0.0409 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 2) + 0.00579 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 3) - 0.000328 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 4);
                
                //c2 = -0.109 - 0.02 * CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths + 0.011 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 2) - 0.00156 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 3) + 0.000121 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 4)
                //Note: c2 from Eq 40 of Hirabayashi and Endreny (2025) and Eq 3-23 of NREL (1995)
                c2 = -0.109 - 0.02 * CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths + 0.011 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 2) - 0.00156 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 3) + 0.000121 * Math.Pow(CloudCover_Opaque_DiffuseHorizontalTransmittance_tenths, 4);
                //Ksopq = -0.06 + b2 * Ta + c2 * Math.Pow(Ta, 2)
                //Note: Ksopq from Eq 38 of Hirabayashi and Endreny (2025) and Eq 3-21 of NREL (1995)
                Ksopq = -0.06 + b2 * Ta + c2 * Math.Pow(Ta, 2);

                if (Ksopq <= 0)
                {
                    Ksopq = 0.0;
                }
            }

            //Note: PSW algorithm from Algorithm 3.5.2 The Precipitation Switch of NREL (1995) 
            //Note: PSW is used as a multiplier to modify the diffuse horizontal value
            if (CloudCover_Opaque_tenths >= 8) {
                if (rain != 0) {
                    PSW_precipitableWaterSwitch = 0.6;
                }
                else {
                    PSW_precipitableWaterSwitch = 1.0;
                }
            }
            else {
                PSW_precipitableWaterSwitch = 1.0;
            }

            //Kd0 = (fM * (Ksr + Ksa) + Ksopq + Kstrn) * PSW_precipitableWaterSwitch
            //Note: Kd0 from Eq 33 of Hirabayashi and Endreny (2025) and within Eq 6 of Maxwell (1998)
            //Note: Kd0 is computed to avoid an implicit function, separating out Ksgrf given Ksgrf is a function of Kd
            Kd0 = (fM * (Ksr + Ksa) + Ksopq + Kstrn) * PSW_precipitableWaterSwitch;

            //Rcld = 0.06 * CloudCover_Opaque_tenths + 0.02 * Trans
            //Note: Rcld from Eq 48 of Hirabayashi and Endreny (2025) and Eq 3-31 of NREL (1995)
            Rcld = 0.06 * CloudCover_Opaque_tenths + 0.02 * Trans;
            //Note: Rcld from Eq 46 of Hirabayashi and Endreny (2025) and Eq 3-33 of NREL (1995)
            Ksgrf1 = (Kn + Kd0) * Rcld * (Albedo_surface_frac - 0.2);

            if (Ksgrf1 <= 0) {
                Ksgrf1 = 0.01;
            }
            //Tas = Ta / Taa
            //Note: Tas from Eq 50 of Hirabayashi and Endreny (2025) and Eq 3-10 of NREL (1995)
            Tas = Ta / Taa;
            //Ratm = (0.0685 + 0.16 * (1 - Tas)) * (10 - CloudCover_Opaque_tenths) / 10.0
            //Note: Ratm from Eq 49 of Hirabayashi and Endreny (2025) and Eq 3-32 of NREL (1995)
            Ratm = (0.0685 + 0.16 * (1 - Tas)) * (10 - CloudCover_Opaque_tenths) / 10.0;
            //Ksgrf2 = (Kn + Kd0) * Ratm * Albedo_surface_frac
            //Note: Ksgrf2 from Eq 46 of Hirabayashi and Endreny (2025) and Eq 3-34 of NREL (1995)
            Ksgrf2 = (Kn + Kd0) * Ratm * Albedo_surface_frac;

            //Ksgrf = Ksgrf1 + Ksgrf2
            //Note: Ksgrf from Eq 45 of Hirabayashi and Endreny (2025) and Eq 3-35 of NREL (1995)
            //Note: Ksgrf is not adjusted for PSW in Eq 6 of Maxwell (1998)
            Ksgrf = Ksgrf1 + Ksgrf2;

            //Kd is effective diffuse horizontal transmittance
            //Kd = Kd0 + Ksgrf
            //Note: Kd from Eq 32 of Hirabayashi and Endreny (2025) and within Eq 6 of Maxwell (1998)
            Kd = Kd0 + Ksgrf;
        }


        static void CalcNetRadWm2(
            double Tair_2m_K,
            double Tdew_2m_C,
            double CloudCover_Total_Tenth,
            bool Flag_SkyData_hasMultipleLayers,
            double CloudCover_1_Tenth,
            double CloudBaseHeight_1_m,
            double CloudCover_2_Tenth,
            double CloudBaseHeight_2_m,
            double CloudCover_3_Tenth,
            double CloudBaseHeight_3_m,
            double Radiation_Shortwave_Total_Wpm2,
            double AtmP_mbar,
            double hourOfDay_local,
            ref double Radiation_Longwave_Downwelling_Wpm2,
            ref double Radiation_Longwave_Upwelling_Wpm2,
            ref double Radiation_Net_Wpm2,
            double Albedo_surface_frac)

        {
            //References: 
            //Hirabayashi, S., & Endreny, T. (2025). Surface and Upper Weather Pre-processor for i-Tree Eco and Hydro. i-Tree Tools White Paper, Version 1.2. 
            //Martin, M., & Berdahl, P. (1984). Characteristics of infrared sky radiation in the United States. Solar Energy, 33(3), 321-336. https://www.sciencedirect.com/science/article/pii/0038092X84901622
            //Stull, R.B. (2017).Practical Meteorology: An Algebra-based Survey of Atmospheric Science. Version 1.02b(Version 1.02b ed.).Univ.of British Columbia: Creative Commons.

            double Emissivity_sky_frac;
            double Radiation_Longwave_Surface_Upwelling_Wpm2;    //Longwave radiation from the surface (W/m^2)
            double Radiation_Longwave_Cloud_Downwelling_Wpm2;    //Longwave radiaton from cloud cover w/ Emissivity set to 1 (W/m^2)
            double Radiation_Longwave_Sky_Downwelling_Wpm2;    //Longwave radiation from clear sky with Emissivity set to Ea*(1-cc) (W/m^2)
            double Radiation_Shortwave_Net_Wpm2;     //Total Albedo_surface_frac corrected incoming shortwave radiation (W/m^2)
            double Radiation_Longwave_Net_Wpm2;    //Total longwave radiation Emissivity_sky_frac corrected (incoming sky - outgoing surface) (W/m^2)
            //Emissivity_cloud_1_frac set to 1 from Martin & Berdahl (1984) Eq 4a when C is 1
            double Emissivity_cloud_1_frac = 1.0;

            double Tair_2m_C = Tair_2m_K - ConversionFactor_273pt15K_to_0C;
            //DewPointDepression_C = Tair_C - Tdew_C
            double DewPointDepression_C = Tair_2m_C - Tdew_2m_C;

            //Note: Sky_height_km of ~ 1.5 to 4 km and EnvironmentalLapseRate_Kpkm of 6.5 K/km create a temperature delta between ...
            //... Sky_Temperature_K Tair_2m_K that generates approximately 98.5 W/m2 net longwave radiation, which is the value ...
            //... Stull (2017) provides as global average net infrared radiation in Eq 2.39 and Fig 2.10.
            //Note: The net longwave radiation results from longwave upwelling radiation - longwave downwelling radiation ...
            //... where it is recognized WeatherPrep cannot estimate ground surface temperature for longwave upwelling radiation ...
            //... and WeatherPrep will use Tair_2m_K as the estimated Tair_surface_K for computing longwave upwelling radiation 

            // LCL in km (approx)
            double z_LCL_km = 0.125 * DewPointDepression_C;

            // Now define clear-sky effective height as:
            // base + factor * (dewpointDep – reference), then clamp.
            double dewRef_C = 10.0;           // a reference dewpoint depression (tunable)
            double baseClear_km = 2.5;        // base effective height (tunable)
            double slopeClear_km_perC = 0.05; // how much height moves per °C (tunable)

            double Sky_height_clear_km =baseClear_km + slopeClear_km_perC * (DewPointDepression_C - dewRef_C);

            // Clamp to a reasonable range, e.g. 1.5–4.0 km
            double Sky_height_min_km = 1.5;
            double Sky_height_max_km = 4.0;

            Sky_height_clear_km = Math.Max(Sky_height_min_km,
                                    Math.Min(Sky_height_max_km, Sky_height_clear_km));
            //EnvironmentalLapseRate_Kpkm set to standard of 6.5 K/km or C/km
            double EnvironmentalLapseRate_Kpkm = 6.5;

            //Sky_Temperature_K 
            double Sky_Temperature_K = Tair_2m_K - (EnvironmentalLapseRate_Kpkm * Sky_height_clear_km);

            double CloudCover_frac = CloudCover_Total_Tenth / 10.0;

            // Convert reported ceiling if present (>0), else use LCL-based estimate
            double Cloud_height_raw_km;

            // If Flag_SkyData_hasMultipleLayers true there are multiple cloud layers (GHCNh), use compressed effective height
            if (Flag_SkyData_hasMultipleLayers)
            {
                double effCover_frac;
                double effHeight_km;

                SurfaceWeather.CompressCloudLayers_Tenths(
                    CloudCover_1_Tenth, CloudBaseHeight_1_m,
                    CloudCover_2_Tenth, CloudBaseHeight_2_m,
                    CloudCover_3_Tenth, CloudBaseHeight_3_m,
                    NODATA,
                    out effCover_frac,
                    out effHeight_km);

                if (effHeight_km != NODATA) {
                    Cloud_height_raw_km = effHeight_km;
                }
                //If CloudBaseHeight_1_m > 0 then clouds present
                else if (CloudBaseHeight_1_m > 0.0) {
                    // Fall back to CLG-based height
                    double Ratio_km_per_m = 0.001;

                    Cloud_height_raw_km = CloudBaseHeight_1_m * Ratio_km_per_m;
                }
                else {
                    // Fallback to LCL
                    Cloud_height_raw_km = z_LCL_km;
                }
            }
            //Else Flag_SkyData_hasMultipleLayers is false, a single cloud layer
            else {
                // Single-layer (ISD or 1-layer GHCNh): use CLG if available, else LCL
                if (CloudBaseHeight_1_m > 0.0) {
                    double Ratio_km_per_m = 0.001;

                    Cloud_height_raw_km = CloudBaseHeight_1_m * Ratio_km_per_m;
                }
                else {
                    Cloud_height_raw_km = z_LCL_km;
                }
            }

            // Clamp cloud height to 1.5–4.0 km (only here, not in CompressCloudLayers_Tenths)
            double Cloud_height_min_km = 1.5;
            double Cloud_height_max_km = 4.0;

            double Cloud_height_km = Math.Max(Cloud_height_min_km,
                                       Math.Min(Cloud_height_max_km, Cloud_height_raw_km));


            //Cloud_Temperature_K is Tair_2m_K - EnvironmentalLapseRate_Kpkm * Cloud_height_km
            //Note: This crude approximation avoids the difficulty of estimating cloud temperatures with cloud physics
            double Cloud_Temperature_K = Tair_2m_K - EnvironmentalLapseRate_Kpkm * Cloud_height_km;

            //Emissivity via Martin & Berdahl (1984) Eq 1 to 3 clear-sky emissivity ----------
            //Note: Emissivity_sky_frac from Eq 58, 58a, 58b, 58c Hirabayashi and Endreny (2025)
            //2025: Removed as undocumented: Emissivity_sky_frac = 0.741 + (0.0062 * Tdew_2m_C);
            //Note: Martin & Berdahl (1984) emissivity as function of dew point (Td in °C) ...
            // ... ε0 = 0.711 + 0.56*(Td/100) + 0.73*(Td/100)^2
            double Emissivity_sky_monthly_frac = 0.711 + 0.56 * (Tdew_2m_C / 100.0) + 0.73 * (Tdew_2m_C / 100.0) * (Tdew_2m_C / 100.0);

            // Diurnal adjustment: Δ_t = 0.013 * cos(2π t / 24)
            double Emissivity_delta_t = 0.013 * Math.Cos(2.0 * Math.PI * hourOfDay_local / 24.0);

            // Pressure adjustment (P in mbar): Δ_P = 0.00012 * (P - 1000)
            double Emissivity_delta_p = 0.00012 * (AtmP_mbar - 1000.0);

            //Emissivity_sky_frac = Emissivity_sky_monthly_frac + Emissivity_delta_t + Emissivity_delta_p;
            Emissivity_sky_frac = Emissivity_sky_monthly_frac + Emissivity_delta_t + Emissivity_delta_p;
            //if (Emissivity_sky_frac > 1 then set to 1
            Emissivity_sky_frac = Math.Max(0.0, Math.Min(1.0, Emissivity_sky_frac));

            //Longwave radiation upwelling from the surface (W/m^2)
            //Note: Radiation_Longwave_Surface_Upwelling_Wpm2 from Eq 59 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Surface_Upwelling_Wpm2 = Emissivity_surface_frac * Stefan_Boltzmann_sigma_Wpm2pK4 * Math.Pow(Tair_2m_K, 4);
            //Longwave radiaton downwelling from cloud cover (W/m2) w/ Emissivity set to 1 (W/m^2)
            //Note: Radiation_Longwave_Cloud_Downwelling_Wpm2 from Eq 57 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Cloud_Downwelling_Wpm2 = Emissivity_cloud_1_frac * (CloudCover_Total_Tenth / 10.0) * Stefan_Boltzmann_sigma_Wpm2pK4 * Math.Pow(Cloud_Temperature_K, 4);
            //Longwave radiation from clear sky (W/m2) with Emissivity set to Ea*(1-cc) (W/m^2)
            //Note: Radiation_Longwave_Sky_Downwelling_Wpm2 from Eq 56 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Sky_Downwelling_Wpm2 = Emissivity_sky_frac * (1 - (CloudCover_Total_Tenth / 10.0)) * Stefan_Boltzmann_sigma_Wpm2pK4 * Math.Pow(Sky_Temperature_K, 4);

            //Total Albedo_surface_frac corrected incoming shortwave radiation (W/m^2)
            //Note: Radiation_Shortwave_Net_Wpm2 from Eq 54 Hirabayashi and Endreny (2025)
            Radiation_Shortwave_Net_Wpm2 = (1 - Albedo_surface_frac) * Radiation_Shortwave_Total_Wpm2;
            //Radiation_Longwave_Net_Wpm2 = (Radiation_Longwave_Sky_Downwelling_Wpm2 + Radiation_Longwave_Cloud_Downwelling_Wpm2) - Radiation_Longwave_Surface_Upwelling_Wpm2
            //Note: Radiation_Longwave_Net_Wpm2 from Eq 55 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Net_Wpm2 = (Radiation_Longwave_Sky_Downwelling_Wpm2 + Radiation_Longwave_Cloud_Downwelling_Wpm2) - Radiation_Longwave_Surface_Upwelling_Wpm2;

            //Total net radation (downwelling shortwave and longwave, Radiation_Net_Wpm2) (W/m^2)
            //Note: Radiation_Net_Wpm2 from Eq 53 Hirabayashi and Endreny (2025)
            Radiation_Net_Wpm2 = Radiation_Shortwave_Net_Wpm2 + Radiation_Longwave_Net_Wpm2;

            //Radiation_Longwave_Downwelling_Wpm2 and Radiation_Longwave_Upwelling_Wpm2 sent back to function calling CalcNetRadWm2 for use in Radiation.csv
            //Note: part of Eq 55 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Downwelling_Wpm2 = (Radiation_Longwave_Sky_Downwelling_Wpm2 + Radiation_Longwave_Cloud_Downwelling_Wpm2);
            //Total longwave radiation upward (W/m2)
            //Note: part of Eq 55 Hirabayashi and Endreny (2025)
            Radiation_Longwave_Upwelling_Wpm2 = Radiation_Longwave_Surface_Upwelling_Wpm2;

        }

        //Function CalcET using descriptive variable names with units
        public static void CalcET(List<SurfaceWeather> weatherData, List<LeafAreaIndex> LeafAreaIndex_list, int Variable_TimeSeries_Hr_Length, double Height_WeatherStationWindSensor_m, double Height_TreeCanopyTop_m)
        {
            //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, Second Edition (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.
            //UCAR (2011).Snowmelt Processes: International Edition. University Corporation for Atmospheric Research COMET MetEd Program. https://download.comet.ucar.edu/memory-stick/hydro/basic_int/snowmelt/index.htm

            int i;
            double VaporPressure_Actual_kPa;
            double VaporPressure_Saturated_kPa;
            double VaporPressure_Deficit_kPa;
            double VaporPressure_Deficit_Pa;
            double VaporPressure_Deficit_OverWater_Pa;
            double WindSpeed_Station_m_p_s;
            double WindSpeed_TreeCanopy_m_p_s;
            double DensityWater_kg_p_m3;
            double DensityAir_kg_p_m3;
            double Radiation_Net_Wpm2;
            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 vonKarman_Constant = 0.41;
            double Ratio_sec_to_hr = 3600.0;
            double Ratio_Pa_to_kPa = 1000.0;
            double Ratio_J_to_MJ = 1000000.0;
            double Ratio_MJ_to_J = 1.0 / 1000000.0;
            //Roughness Lengths from Davenport-Wieringa Zo (m) to Classification and Landscape Table 4.1 in Stull (2000)
            //Note: Classificaiton smooth = Snow; Classificaiton Open = Airport; Note UCAR Community Land Model uses 0.024 
            double RoughnessLength_Snow_m = 0.005;
            double RoughnessLength_Airport_m = 0.03;
            //Roughness Lengths for open water from Shuttleworth (1993) page 4.15, above Eq 4.2.30, citing Thom and Oliver (1977)
            double RoughnessLength_Water_m = 0.00137;
            //ZeroPlaneDisplacementHeight_frac (fraction) from Eq 13.3 Chin (2021) adjusts cover height to get zero plane displacement height
            double ZeroPlaneDisplacementHeight_frac = 0.67;
            //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_TreeCanopyTop_m;
            double Height_TreeZeroPlaneDisplacement_m = Height_TreeCanopyTop_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 LAI_plus_BAI_plus_GAI_Tree_to_SVeg_ratio;
            double LAI_plus_BAI_plus_GAI_Ground_to_Tree_ratio;
            //BarkAreaIndex_Tree_m2_p_m2 for trees is estimated at 1.7 from HydroPlus test cases
            double BarkAreaIndex_Tree_m2_p_m2 = 1.7;
            //LeafAreaIndex_SVeg_m2_p_m2 for shrubs is estimated at 2.3 from HydroPlus test cases
            double LeafAreaIndex_SVeg_m2_p_m2 = 2.2;
            //BarkAreaIndex_SVeg_m2_p_m2 for shrubs is estimated at 0.5 from HydroPlus test cases
            double BarkAreaIndex_SVeg_m2_p_m2 = 0.5;
            //GroundAreaIndex_m2_p_m2 for trees, shrubs, and ground is estimated at 1
            double GroundAreaIndex_m2_p_m2 = 1;

            //SpecificHeat_MoistAir_J_p_kg_p_C is specific heat of moist air at constant pressure, defined by Chin (2021) w Eq 13.37
            double SpecificHeat_MoistAir_J_p_kg_p_C = 1013;
            double LatentHeatOfVaporizaton_J_p_kg;
            double LatentHeatOfVaporizaton_MJ_p_kg;
            double PsychrometricConstant_Pa_p_C;
            double PsychrometricConstant_kPa_p_C;
            double VaporPressureGradient_Pa;
            double VaporPressureGradient_kPa;
            double Ratio_MolecularWeightVapor_to_MolecularWeightDryAir = 0.622;
            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;


            try {
                //For loop through Variable_TimeSeries_Hr_Length, the length of the input data
                for (i = 0; i < Variable_TimeSeries_Hr_Length; i++) {

                    //VaporPressure_Actual_kPa is actual vapor pressure (kPa)
                    VaporPressure_Actual_kPa = 0.6108 * Math.Exp((17.27 * weatherData[i].Tdew_2m_C) / (237.3 + weatherData[i].Tdew_2m_C));
                    VaporPressure_Saturated_kPa = 0.6108 * Math.Exp((17.27 * weatherData[i].Tair_2m_C) / (237.3 + weatherData[i].Tair_2m_C));
                    //VaporPressure_Deficit_kPa is vapor pressure deficit (kPa), using previously computed saturated vapor pressure VapPres_Sat_kPa (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 Tair_2m_C > freezing then reduce vapor pressure deficit due to high humidity of air by evaporation; otherwise no adjustment
                    /*
                    if (weatherData[i].Tair_2m_C > 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].Tair_2m_C / 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;
                    }
                    */

                    //WindSpeed_Station_m_p_s is wind speed at station (m/s)
                    WindSpeed_Station_m_p_s = weatherData[i].Wspd_mps;
                    //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) using Eq 4.14b of Stull (2000)
                    WindSpeed_TreeCanopy_m_p_s = WindSpeed_Station_m_p_s * (Math.Log(Height_TreeWindSensor_m / RoughnessLength_Airport_m) / Math.Log(Height_WeatherStationWindSensor_m / RoughnessLength_Airport_m));

                    //DensityWater_kg_p_m3 is density of water (kg/m^3) from Equation in Figure 11.1.1 McCutcheon (1993)
                    DensityWater_kg_p_m3 = 1000 * (1 - ((weatherData[i].Tair_2m_C + 288.9414) / (508929.2 * (weatherData[i].Tair_2m_C + 68.12963))) *
                        Math.Pow((weatherData[i].Tair_2m_C - 3.9863), 2));

                    //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 * (weatherData[i].AtmP_kPa / (273.15 + weatherData[i].Tair_2m_C));

                    //LatentHeatOfVaporizaton_MJ_p_kg (MJ/kg) from Eq 13.36 Chin (2021) is latent heat of vaporization
                    LatentHeatOfVaporizaton_MJ_p_kg = 2.501 - 0.002361 * weatherData[i].Tair_2m_C;
                    //LatentHeatOfVaporizaton_J_p_kg (J/kg) is converted from MJ to K
                    LatentHeatOfVaporizaton_J_p_kg = LatentHeatOfVaporizaton_MJ_p_kg * Ratio_J_to_MJ;

                    //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 * weatherData[i].AtmP_kPa) /
                        (Ratio_MolecularWeightVapor_to_MolecularWeightDryAir * LatentHeatOfVaporizaton_MJ_p_kg);
                    //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 * Math.Exp((17.27 * weatherData[i].Tair_2m_C) / (weatherData[i].Tair_2m_C + 237.3))) / Math.Pow((weatherData[i].Tair_2m_C + 237.3), 2);
                    //VaporPressureGradient_Pa (Pa) is converted from kPa to Pa
                    VaporPressureGradient_Pa = VaporPressureGradient_kPa * Ratio_Pa_to_kPa;

                    //Radiation_Net_Wpm2 (W/m2) is energy from downwelling shortwave and longwave radiation (W/m^2)
                    Radiation_Net_Wpm2 = weatherData[i].Radiation_Net_Wpm2;

                    //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
                    //Mote: 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 = LeafAreaIndex_list[i].Lai * 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 evoled 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 = Math.Pow((1 - VaporPressure_Deficit_kPa / weatherData[i].VapPres_Sat_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 (LeafAreaIndex_list[i].Lai < 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;

                    //Note: These terms used in HydroPlus to compute evaporation for different cover types
                    //LAI_plus_BAI_plus_GAI_Tree_to_SVeg_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_Tree_to_SVeg_ratio = (LeafAreaIndex_list[i].Lai + BarkAreaIndex_Tree_m2_p_m2 + GroundAreaIndex_m2_p_m2) /
                        (LeafAreaIndex_SVeg_m2_p_m2 + BarkAreaIndex_SVeg_m2_p_m2 + GroundAreaIndex_m2_p_m2);
                    //LAI_plus_BAI_plus_GAI_Ground_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_Ground_to_Tree_ratio = (0 + 0 + GroundAreaIndex_m2_p_m2) /
                        (LeafAreaIndex_list[i].Lai + BarkAreaIndex_Tree_m2_p_m2 + GroundAreaIndex_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_WeatherStationWindSensor_m is used in place of Height_WeatherStationWindSensor_m to ensure real values from log
                    ResistanceAerodynamic_TreeCanopy_s_p_m = (Math.Log((Height_TreeWindSensor_m - Height_TreeZeroPlaneDisplacement_m) /
                        (rZom * Height_TreeCanopyTop_m)) * Math.Log((Height_TreeWindSensor_m - Height_TreeZeroPlaneDisplacement_m) /
                        (rZov * Height_TreeCanopyTop_m))) / (Math.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 * Math.Pow((Math.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 & Eq 3 of Light (1941)
                    ResistanceAerodynamic_SnowTreeCanopy_s_p_m = (Math.Log(Height_TreeWindSensor_m / RoughnessLength_Snow_m) *
                        Math.Log(Height_TreeWindSensor_m / RoughnessLength_Snow_m)) / (Math.Pow(vonKarman_Constant, 2) * WindSpeed_TreeCanopy_m_p_s);

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

                    //PtTrMh (m/h), evapotranspiration from canopy and ground, from Eq 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, usng mixed units MJ/day, kPa/C, etc.
                    weatherData[i].PtTrMh = 1 / (LatentHeatOfVaporizaton_J_p_kg * DensityWater_kg_p_m3) *
                        ((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));
                    //PtTrMh (m/h) converted from (m/s) with Ratio_sec_to_hr
                    weatherData[i].PtTrMh = weatherData[i].PtTrMh * Ratio_sec_to_hr;
                    //If negative pet to zero
                    if (weatherData[i].PtTrMh < 0) { weatherData[i].PtTrMh = 0.0; }

                    //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, (aka shortwave radiation) 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 longwave radiation at Julian Day = 180; not simulated for smaller waters
                    //Note: Modification includes Heat Storage transfering 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 * (weatherData[i].Radiation_Shortwave_Direct_Wpm2 + weatherData[i].Radiation_Shortwave_Diffuse_Wpm2) -
                        HeatStorage_Coefficient_b * weatherData[i].Radiation_Longwave_Upwelling_Wpm2;

                    //PeGrMh (m/h), evaporation from water surface on ground, from Eq 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, usng mixed units MJ/day, kPa/C, etc.
                    //Note: ResistanceSurface_Water_s_p_m set to zero. 
                    weatherData[i].PeGrMh = 1 / (LatentHeatOfVaporizaton_J_p_kg * DensityWater_kg_p_m3) *
                        ((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));
                    //PeGrMh (m/h) converted from (m/s) with Ratio_sec_to_hr
                    weatherData[i].PeGrMh = weatherData[i].PeGrMh * Ratio_sec_to_hr;
                    //If negative pet to zero
                    if (weatherData[i].PeGrMh < 0) { weatherData[i].PeGrMh = 0.0; }

                    //PeTrMh (m/s), evaporation from canopy, set equal to PeGrMh  (m/s)
                    //Note: Variable name contains Mh but it was previously converted to m/s
                    weatherData[i].PeTrMh = weatherData[i].PeGrMh;

                    //PeSnGrMh (m/h), 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
                    weatherData[i].PeSnGrMh = 1 / (LatentHeatOfVaporizaton_J_p_kg * DensityWater_kg_p_m3) *
                        ((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));
                    //PeSnGrMh (m/h) converted from (m/s) with Ratio_sec_to_hr
                    weatherData[i].PeSnGrMh = weatherData[i].PeSnGrMh * Ratio_sec_to_hr;
                    //If negative pet to zero
                    if (weatherData[i].PeSnGrMh < 0) { weatherData[i].PeSnGrMh = 0.0; }

                    //PeSnTrMh (m/h), 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.
                    weatherData[i].PeSnTrMh = 1 / (LatentHeatOfVaporizaton_J_p_kg * DensityWater_kg_p_m3) *
                        ((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));
                    //PeSnTrMh (m/h) converted from (m/s) with Ratio_sec_to_hr
                    weatherData[i].PeSnTrMh = weatherData[i].PeSnTrMh * Ratio_sec_to_hr;
                    //If negative pet to zero
                    if (weatherData[i].PeSnTrMh < 0) { weatherData[i].PeSnTrMh = 0.0; }
                }
            }
            catch (Exception) {
                throw;
            }
        }
        public static void CalcPrecipInterceptByCanopy(List<SurfaceWeather> wList, List<LeafAreaIndex> LeafAreaIndex_list, int Variable_TimeSeries_Hr_Length, VEGTYPE vegType)
        {
            int i;
            double maxSt;       // maximum canopy storage of precipitation
            double canopyCover_frac;           // canopy cover fraction
            double potEvMh;     // potential evaporation for tree or ground (for vegType=Tree and Shrub, respectively)
            double canPrcp;     // precipitation fallen and contacted to canopy at the time step
            double freeThf;     // precipitation fallen through canopy without contact (free throughfall) at the time step
            double canDrip;     // drip of water from canopy in the 2nd stage
            double preSt;       // canopy storage of precipitation at the previous time step
            double preEv;       // evaporation from canopy at the previous time step

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                // previous hour condition 
                preSt = (i == 0) ? 0 : wList[i - 1].VegStMh;
                preEv = (i == 0) ? 0 : wList[i - 1].VegEvMh;

                potEvMh = (vegType == VEGTYPE.TREE) ? wList[i].PeTrMh : wList[i].PeGrMh;
                canopyCover_frac = (vegType == VEGTYPE.TREE) ? 1 - Math.Exp(-0.7 * LeafAreaIndex_list[i].Lai) : 1 - Math.Exp(-0.3 * LeafAreaIndex_list[i].Lai);

                maxSt = LEAF_STORAGE_M * LeafAreaIndex_list[i].Lai;

                freeThf = wList[i].Ppt_1hr_m * (1 - canopyCover_frac);
                canPrcp = wList[i].Ppt_1hr_m - freeThf;

                // storage
                wList[i].VegStMh = preSt + canPrcp - preEv;
                if (wList[i].VegStMh > maxSt)
                {
                    wList[i].VegStMh = maxSt;
                }
                else if (wList[i].VegStMh < 0)
                {
                    wList[i].VegStMh = 0;
                }

                // evaporation
                if (maxSt == 0) // when evergreen % = 0, LAI = 0 in the winter and thus maxSt = 0
                {
                    wList[i].VegEvMh = 0;
                }
                else
                {
                    wList[i].VegEvMh = Math.Pow(wList[i].VegStMh / maxSt, 2.0 / 3.0) * potEvMh;
                    // cap evaporation with storage
                    if (wList[i].VegEvMh > wList[i].VegStMh)
                    {
                        wList[i].VegEvMh = wList[i].VegStMh;
                    }
                }

                // drip
                canDrip = 0;
                if (wList[i].VegStMh < maxSt)                     // 1st stage
                {
                    //                    wList[i].ThrufallMh = freeThf - wList[i].VegEvMh;   // should be this.
                    wList[i].UnderCanThrufallMh = freeThf;                       // but for compatibility with Hydro, use this
                }
                else
                {
                    if (preSt < maxSt)                          // 2nd stage for the first time
                    {
                        canDrip = canPrcp - (maxSt - preSt) - preEv;
                        if (canDrip < 0)
                        {
                            canDrip = 0;
                        }
                    }
                    else                                       // 2nd stage
                    {
                        canDrip = canPrcp - preEv;
                        if (canDrip < 0)
                        {
                            canDrip = 0;
                        }
                    }
                    wList[i].UnderCanThrufallMh = canDrip + freeThf;
                }

                // interception
                // While it's raining hourly interception is calculated from input - output of water
                // (i.e., Rain on canopy - Under canopy throughfall
                if (wList[i].Ppt_1hr_inches > 0)
                {
                    wList[i].VegIntcptMh = wList[i].Ppt_1hr_m - wList[i].UnderCanThrufallMh;
                }
                else
                {
                    wList[i].VegIntcptMh = 0;
                }
            }
        }
        public static void CalcPrecipInterceptByUnderCanopyCover(List<SurfaceWeather> wList, int Variable_TimeSeries_Hr_Length)
        {
            int i;
            double maxPervSt;       // maximum pervious cover storage of precipitation
            double maxImpervSt;     // maximum impervious cover storage of precipitation
            double potEvMh;         // potential evaporation for ground
            double prePervSt;       // pervious storage of precipitation at the previous time step
            double preImpervSt;     // impervious storage of precipitation at the previous time step
            double prePervEv;       // pervious evaporation at the previous time step
            double preImpervEv;     // impervious evaporation at the previous time step

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                // previous hour condition 
                prePervSt = (i == 0) ? 0 : wList[i - 1].UnderCanPervStMh;
                preImpervSt = (i == 0) ? 0 : wList[i - 1].UnderCanImpervStMh;
                prePervEv = (i == 0) ? 0 : wList[i - 1].UnderCanPervEvMh;
                preImpervEv = (i == 0) ? 0 : wList[i - 1].UnderCanImpervEvMh;

                potEvMh = wList[i].PeGrMh;

                maxPervSt = PERV_STORAGE_M;
                maxImpervSt = IMPERV_STORAGE_M;

                // storage
                //  pervious cover
                wList[i].UnderCanPervStMh = prePervSt + wList[i].UnderCanThrufallMh - prePervEv;
                if (wList[i].UnderCanPervStMh > maxPervSt)
                {
                    wList[i].UnderCanPervStMh = maxPervSt;
                }
                else if (wList[i].UnderCanPervStMh < 0)
                {
                    wList[i].UnderCanPervStMh = 0;
                }
                //  impervious cover
                wList[i].UnderCanImpervStMh = preImpervSt + wList[i].UnderCanThrufallMh - preImpervEv;
                if (wList[i].UnderCanImpervStMh > maxImpervSt)
                {
                    wList[i].UnderCanImpervStMh = maxImpervSt;
                }
                else if (wList[i].UnderCanImpervStMh < 0)
                {
                    wList[i].UnderCanImpervStMh = 0;
                }

                // evaporation
                //  pervious cover
                wList[i].UnderCanPervEvMh = (wList[i].UnderCanPervStMh / maxPervSt) * potEvMh;
                //  cap evaporation with storage
                if (wList[i].UnderCanPervEvMh > wList[i].UnderCanPervStMh)
                {
                    wList[i].UnderCanPervEvMh = wList[i].UnderCanPervStMh;
                }
                //  impervious cover
                wList[i].UnderCanImpervEvMh = (wList[i].UnderCanImpervStMh / maxImpervSt) * potEvMh;
                //  cap evaporation with storage
                if (wList[i].UnderCanImpervEvMh > wList[i].UnderCanImpervStMh)
                {
                    wList[i].UnderCanImpervEvMh = wList[i].UnderCanImpervStMh;
                }

                // runoff
                wList[i].UnderCanImpervRunoffMh = 0;
                if (preImpervSt < maxImpervSt)
                {
                    wList[i].UnderCanImpervRunoffMh = wList[i].UnderCanThrufallMh - (maxImpervSt - preImpervSt) - wList[i].UnderCanImpervEvMh;
                }
                else if (preImpervSt == maxImpervSt)
                {
                    wList[i].UnderCanImpervRunoffMh = wList[i].UnderCanThrufallMh - wList[i].UnderCanImpervEvMh;
                }
                if (wList[i].UnderCanImpervRunoffMh < 0)
                {
                    wList[i].UnderCanImpervRunoffMh = 0;
                }

                // infiltration
                wList[i].UnderCanPervInfilMh = 0;
                if (prePervSt < maxPervSt)
                {
                    wList[i].UnderCanPervInfilMh = wList[i].UnderCanThrufallMh - (maxPervSt - prePervSt) - wList[i].UnderCanPervEvMh;
                }
                else if (prePervSt == maxPervSt)
                {
                    wList[i].UnderCanPervInfilMh = wList[i].UnderCanThrufallMh - wList[i].UnderCanPervEvMh;
                }
                if (wList[i].UnderCanPervInfilMh < 0)
                {
                    wList[i].UnderCanPervInfilMh = 0;
                }
            }
        }
        public static void CalcPrecipInterceptByNoCanopyCover(List<SurfaceWeather> wList, int Variable_TimeSeries_Hr_Length)
        {
            int i;
            double maxPervSt;       // maximum pervious cover storage of precipitation
            double maxImpervSt;     // maximum impervious cover storage of precipitation
            double potEvMh;         // potential evaporation for ground
            double prePervSt;       // pervious storage of precipitation at the previous time step
            double preImpervSt;     // impervious storage of precipitation at the previous time step
            double prePervEv;       // pervious evaporation at the previous time step
            double preImpervEv;     // impervious evaporation at the previous time step

            for (i = 0; i < Variable_TimeSeries_Hr_Length; i++)
            {
                // previous hour condition 
                prePervSt = (i == 0) ? 0 : wList[i - 1].NoCanPervStMh;
                preImpervSt = (i == 0) ? 0 : wList[i - 1].NoCanImpervStMh;
                prePervEv = (i == 0) ? 0 : wList[i - 1].NoCanPervEvMh;
                preImpervEv = (i == 0) ? 0 : wList[i - 1].NoCanImpervEvMh;

                potEvMh = wList[i].PeGrMh;

                maxPervSt = PERV_STORAGE_M;
                maxImpervSt = IMPERV_STORAGE_M;

                // storage
                //  pervious cover
                wList[i].NoCanPervStMh = prePervSt + wList[i].Ppt_1hr_m - prePervEv;
                if (wList[i].NoCanPervStMh > maxPervSt)
                {
                    wList[i].NoCanPervStMh = maxPervSt;
                }
                else if (wList[i].NoCanPervStMh < 0)
                {
                    wList[i].NoCanPervStMh = 0;
                }
                //  impervious cover
                wList[i].NoCanImpervStMh = preImpervSt + wList[i].Ppt_1hr_m - preImpervEv;
                if (wList[i].NoCanImpervStMh > maxImpervSt)
                {
                    wList[i].NoCanImpervStMh = maxImpervSt;
                }
                else if (wList[i].NoCanImpervStMh < 0)
                {
                    wList[i].NoCanImpervStMh = 0;
                }

                // evaporation
                //  pervious cover
                wList[i].NoCanPervEvMh = (wList[i].NoCanPervStMh / maxPervSt) * potEvMh;
                //  cap evaporation with storage
                if (wList[i].NoCanPervEvMh > wList[i].NoCanPervStMh)
                {
                    wList[i].NoCanPervEvMh = wList[i].NoCanPervStMh;
                }
                //  impervious cover
                wList[i].NoCanImpervEvMh = (wList[i].NoCanImpervStMh / maxImpervSt) * potEvMh;
                //  cap evaporation with storage
                if (wList[i].NoCanImpervEvMh > wList[i].NoCanImpervStMh)
                {
                    wList[i].NoCanImpervEvMh = wList[i].NoCanImpervStMh;
                }

                // runoff
                wList[i].NoCanImpervRunoffMh = 0;
                if (preImpervSt < maxImpervSt)
                {
                    wList[i].NoCanImpervRunoffMh = wList[i].Ppt_1hr_m - (maxImpervSt - preImpervSt) - wList[i].NoCanImpervEvMh;
                }
                else if (preImpervSt == maxImpervSt)
                {
                    wList[i].NoCanImpervRunoffMh = wList[i].Ppt_1hr_m - wList[i].NoCanImpervEvMh;
                }
                if (wList[i].NoCanImpervRunoffMh < 0)
                {
                    wList[i].NoCanImpervRunoffMh = 0;
                }

                // infiltration
                wList[i].NoCanPervInfilMh = 0;
                if (prePervSt < maxPervSt)
                {
                    wList[i].NoCanPervInfilMh = wList[i].Ppt_1hr_m - (maxPervSt - prePervSt) - wList[i].NoCanPervEvMh;
                }
                else if (prePervSt == maxPervSt)
                {
                    wList[i].NoCanPervInfilMh = wList[i].Ppt_1hr_m - wList[i].NoCanPervEvMh;
                }
                if (wList[i].NoCanPervInfilMh < 0)
                {
                    wList[i].NoCanPervInfilMh = 0;
                }
            }
        }
        public static void CopySurfaceWeatherData(string inFile, string outFile, bool blCopy)
        {
            StreamReader sr;
            StreamWriter sw;
            string lineOrg = "";
            string line = "";
            WeatherDataFormat WeatherData_Format;

            using (sr = new StreamReader(inFile))
            {
                WeatherData_Format = CheckWeatherDataFormat(inFile);
                try
                {
                    //First file to be copied: copy the entire input file to the output file.
                    if (blCopy)
                    {
                        using (sw = new StreamWriter(outFile, false))
                        {
                            while ((lineOrg = sr.ReadLine()) != null)
                            {
                                switch (WeatherData_Format)
                                {
                                    case WeatherDataFormat.ISD_intl:
                                        line = lineOrg.Substring(0, 66) + lineOrg.Substring(75, 66);
                                        break;
                                    case WeatherDataFormat.ISD_us_can:
                                        line = lineOrg.Substring(0, 66) + lineOrg.Substring(81, 66);
                                        break;
                                    case WeatherDataFormat.ISD_old:
                                        line = lineOrg;
                                        break;
                                    case WeatherDataFormat.NARCCAP:
                                        line = lineOrg;
                                        break;
                                    default:
                                        break;
                                }
                                sw.WriteLine(line);
                            }
                        }
                    }
                    //Append the file by excluding the first header line.
                    else
                    {
                        using (sw = new StreamWriter(outFile, true))
                        {
                            //Read the header line
                            lineOrg = sr.ReadLine();
                            while ((lineOrg = sr.ReadLine()) != null)
                            {
                                switch (WeatherData_Format)
                                {
                                    case WeatherDataFormat.ISD_intl:
                                        line = lineOrg.Substring(0, 66) + lineOrg.Substring(75, 66);
                                        break;
                                    case WeatherDataFormat.ISD_us_can:
                                        line = lineOrg.Substring(0, 66) + lineOrg.Substring(81, 66);
                                        break;
                                    case WeatherDataFormat.ISD_old:
                                        line = lineOrg;
                                        break;
                                    default:
                                        break;
                                }
                                sw.WriteLine(line);
                            }
                        }
                    }
                }
                catch (Exception)
                {
                    throw;
                }
            }
        }
    }
}
