Climate.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_CLIMATE_HPP
#define OPENSWMM_CLIMATE_HPP
#include "Radiation.hpp"
 
namespace openswmm {
 
struct SimulationContext;
 
namespace climate {
 
// ============================================================================
// Constants
// ============================================================================
 
constexpr double MM_PER_INCH = 25.40;
 
// ============================================================================
// Evaporation method
// ============================================================================
 
enum class EvapMethod : int {
    CONSTANT = 0,
    MONTHLY = 1,
    TIMESERIES = 2,
    TEMPERATURE = 3, 
    PAN = 4
};
 
// ============================================================================
// 7-day temperature moving average (for Hargreaves ET)
// ============================================================================
 
struct MovingAvg7 {
    double ta[7] = {}; 
    double tr[7] = {}; 
    int front = 0;     
    int count = 0;     

    void push(double t_avg, double t_range);

    double avg_temp() const;

    double avg_range() const;
};
 
// ============================================================================
// Daily climate state (scalar — broadcast to all subcatchments)
// ============================================================================
 
struct ClimateState {
    RadiationState radiation;
    double temperature = 70.0; 
    double temp_range = 0.0;   
    double evap_rate = 0.0;    
    double wind_speed = 0.0;   
    double humidity = 50.0;    
 
    // Source/default bases (pre-forcing). The per-step source lookup writes
    // these; forcing then resolves effective values on top. Keeping them
    // separate lets a one-shot/cleared forcing revert to the data source (or
    // default) even when no source overwrites the value each step.
    double temperature_src = 70.0; 
    double file_tmin = 70.0;       
    double file_tmax = 70.0;       
    double wind_speed_src = 0.0;   
 
    // Derived values
    double gamma = 0.0; 
    double ea = 0.0;    
 
    // Hargreaves parameters
    double latitude = 0.0; 
 
    // Site elevation for psychrometric constant (matching legacy Temp.elev)
    double elev = 0.0;       
    double snow_divt = 34.0; 
 
    // Monthly evaporation table (for MONTHLY method)
    double monthly_evap[12] = {};
 
    // EVAPRATE unit-conversion factor (in/day US or mm/day SI → ft/sec) for the
    // project's unit system. Set at init; MONTHLY uses it so SI models convert
    // mm/day correctly (defaults to the US factor for back-compatibility).
    double evaprate_ucf = 1036800.0;
 
    // Monthly adjustment factors
    double adjust_evap[12] = {
        0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0}; 
    // Legacy climate.c NextEvapDate / NextEvapRate: the date the evaporation
    // source next changes (every type — a year ahead for CONSTANT, the first
    // of next month for MONTHLY, the next series entry for TIMESERIES, the
    // next day for a climate file) and the rate that entry carries. A
    // time-series evaporation is STEP-WISE: the rate holds until the date
    // is reached, and the runoff step never crosses it (runoff_getTimeStep).
    double next_evap_date = 0.0;
    double next_evap_rate = 0.0; 
    int evap_ts_pos =
        -1; 
    double adjust_temp[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    double adjust_rain[12] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
    double adjust_hydcon[12] = {1, 1, 1, 1, 1, 1,
                                1, 1, 1, 1, 1, 1}; 
 
    EvapMethod evap_method = EvapMethod::CONSTANT;
 
    // Monthly adjustment for infiltration capacity (scales f0/fmin/Ks)
    // @see Legacy: InfilFactor (set from hydcon adjustment pattern)
    double infil_factor = 1.0;
 
    // 7-day moving average for Hargreaves ET (matching legacy TMovAve Tma)
    MovingAvg7 temp_ma;
 
    // Recovery factor for infiltration models
    // @see Legacy: Evap.recoveryFactor (set from recovery adjustment pattern)
    double recovery_factor = 1.0;
 
    // Resolved timeseries/pattern indices (set at init, -1 = none)
    int temp_ts_index = -1;      
    int evap_ts_index = -1;      
    int humidity_ts_index = -1;  
    int recovery_pat_index = -1; 
 
    // Sub-daily sinusoidal temperature interpolation (Gap #9)
    // Valid only when temperature comes from a daily min/max source.
    // @see Legacy: setTemp() / updateTempTimes() in climate.c
    double tmin_daily = 0.0; 
    double tmax_daily = 0.0; 
    double prev_tmax = 0.0;  
    double hrsr = 6.0;       
    double hrss = 15.0;      
    double hrday = 10.5;     
    double dhrdy = -9.0;     
    double dydif = 15.0;     
    double dtlong = 0.0;     
    bool has_minmax = false; 
    int last_temp_doy = -1;  
};
 
// ============================================================================
// Functions
// ============================================================================

double hargreaves(double latitude, int day_of_year, double t_avg, double t_range);

void updateDailyClimate(ClimateState& state, int day_of_year, int month);

void updateTempTimes(ClimateState& state, int doy);

double getSubdailyTemp(const ClimateState& state, double hour);

void batchDistributeEvap(double evap_rate, const double* ponded_depth, double* evap_out, int n,
                         double dt);
 
} // namespace climate
} // namespace openswmm
 
#endif // OPENSWMM_CLIMATE_HPP