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
 
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 {
    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 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;   
 
    // 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]   = {1,1,1,1,1,1,1,1,1,1,1,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