-
Notifications
You must be signed in to change notification settings - Fork 46
/
temperature_component.cpp
248 lines (205 loc) · 10.7 KB
/
temperature_component.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
/*
* temp_component.cpp
* hector
*
* Created by Ben on 05/25/2011.
*
*/
#include "components/temperature_component.hpp"
#include "core/core.hpp"
#include "core/dependency_finder.hpp"
#include "h_util.hpp"
#include "visitors/avisitor.hpp"
namespace Hector {
using namespace std;
//------------------------------------------------------------------------------
/*! \brief Constructor
*/
TemperatureComponent::TemperatureComponent() : internal_Ftot(0.0), last_Ftot(0.0) {
}
//------------------------------------------------------------------------------
/*! \brief Destructor
*/
TemperatureComponent::~TemperatureComponent() {
}
//------------------------------------------------------------------------------
// documentation is inherited
string TemperatureComponent::getComponentName() const
{
const string name = TEMPERATURE_COMPONENT_NAME;
return name;
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::init( Core* coreptr ) {
logger.open( getComponentName(), false, Logger::DEBUG );
H_LOG( logger, Logger::DEBUG ) << "hello " << getComponentName() << std::endl;
tgaveq.set( 0.0, U_DEGC, 0.0 );
tgav.set( 0.0, U_DEGC, 0.0 );
core = coreptr;
tgav_constrain.allowInterp( true );
tgav_constrain.name = D_TGAV_CONSTRAIN;
FCO2_record.allowInterp( true );
// Register the data we can provide
core->registerCapability( D_GLOBAL_TEMP, getComponentName() );
core->registerCapability( D_GLOBAL_TEMPEQ, getComponentName() );
// Register our dependencies
core->registerDependency( D_RF_TOTAL, getComponentName() );
}
//------------------------------------------------------------------------------
// documentation is inherited
unitval TemperatureComponent::sendMessage( const std::string& message,
const std::string& datum,
const message_data info ) throw ( h_exception )
{
unitval returnval;
if( message==M_GETDATA ) { //! Caller is requesting data
return getData( datum, info.date );
} else if( message==M_SETDATA ) { //! Caller is requesting to set data
H_THROW("Temperature: sendMessage not yet implemented for message=M_SETDATA.");
//TODO: call setData below
//TODO: change core so that parsing is routed through sendMessage
//TODO: make setData private
} else { //! We don't handle any other messages
H_THROW( "Caller sent unknown message: "+message );
}
return returnval;
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::setData( const string& varName,
const message_data& data ) throw ( h_exception )
{
H_LOG( logger, Logger::DEBUG ) << "Setting " << varName << "[" << data.date << "]=" << data.value_str << std::endl;
try {
if( varName == D_ECS ) {
H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" );
S = unitval::parse_unitval( data.value_str, data.units_str, U_DEGC );
} else if( varName == D_TGAV_CONSTRAIN ) {
H_ASSERT( data.date != Core::undefinedIndex(), "date required" );
tgav_constrain.set( data.date, unitval::parse_unitval( data.value_str, data.units_str, U_DEGC ) );
} else if( varName == D_SO2I_B ) {
H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" );
so2i_b = unitval::parse_unitval( data.value_str, data.units_str, U_UNITLESS );
} else if( varName == D_SO2D_B ) {
H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" );
so2d_b = unitval::parse_unitval( data.value_str, data.units_str, U_UNITLESS );
} else if( varName == D_OC_B ) {
H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" );
oc_b = unitval::parse_unitval( data.value_str, data.units_str, U_UNITLESS );
} else if( varName == D_BC_B ) {
H_ASSERT( data.date == Core::undefinedIndex() , "date not allowed" );
bc_b = unitval::parse_unitval( data.value_str, data.units_str, U_UNITLESS );
} else {
H_THROW( "Unknown variable name while parsing " + getComponentName() + ": "
+ varName );
}
} catch( h_exception& parseException ) {
H_RETHROW( parseException, "Could not parse var: "+varName );
}
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::prepareToRun() throw ( h_exception ) {
H_LOG( logger, Logger::DEBUG ) << "prepareToRun " << std::endl;
if( tgav_constrain.size() ) {
Logger& glog = Logger::getGlobalLogger();
H_LOG( glog, Logger::WARNING ) << "Temperature will be overwritten by user-supplied values!" << std::endl;
}
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::run( const double runToDate ) throw ( h_exception ) {
// We track total radiative forcing using internal variable `internal_Ftot`
// Need to do this because if we're subject to a user constraint (being forced
// to match a temperature record), need to track the Ftot that *would* have
// produced the observed temperature record. This way there's a smooth
// transition when we exit the constraint period, after which internal_Ftot
// will rise in parallel with the value reported by ForcingComponent.
// If we never had any temperature constraint, `internal_Ftot` will match `Ftot`.
// If the user has supplied temperature data, use that (except if past its end)
if( tgav_constrain.size() && runToDate <= tgav_constrain.lastdate() ) {
H_LOG( logger, Logger::NOTICE ) << "** Using user-supplied temperature" << std::endl;
tgav = tgav_constrain.get( runToDate );
tgaveq = tgav;
// Update last_Ftot and internal_Ftot for when we exit constraint
last_Ftot = core->sendMessage( M_GETDATA, D_RF_TOTAL ).value( U_W_M2 );
const double heatflux = core->sendMessage( M_GETDATA, D_HEAT_FLUX ).value( U_W_M2 );
internal_Ftot = 3.7 * tgav / S + heatflux; // inverting normal T=S/3.7*F eqn
} else { // not constrained
// Not subject to a constraint. Get total forcing and update our
// internal tracking variables `internal_Ftot` (parallels Ftot, may
// be identical if never had a constraint) and `last_Ftot`.
double Ftot = core->sendMessage( M_GETDATA, D_RF_TOTAL ).value( U_W_M2 );
// A few forcing agents occur disproportionately over land, and thus exert
// a stronger effect than their global mean RF would suggest. In the future,
// we want to compute RF separately for land and ocean; for now, add a 'bonus'
// for these forcing agents.
// NOTE: this is commented out and not used for now.
/*
const double fbc = core->sendMessage( M_GETDATA, D_RF_BC ).value( U_W_M2 );
Ftot += fbc * bc_b; // % extra
const double foc = core->sendMessage( M_GETDATA, D_RF_OC ).value( U_W_M2 );
Ftot += foc * oc_b; // % extra
const double fso2d = core->sendMessage( M_GETDATA, D_RF_SO2d ).value( U_W_M2 );
Ftot += fso2d * so2d_b; // % extra
const double fso2i = core->sendMessage( M_GETDATA, D_RF_SO2i ).value( U_W_M2 );
Ftot += fso2i * so2i_b; // % extra
*/
// CO2 is subject to a lag effect, coded empirically here because a simple model
// like Hector doesn't account for all buffers and processes in the earth system
// We use a window (size CO2_WINDOW) CO2_LAG years in the past
const double FCO2_current = core->sendMessage( M_GETDATA, D_RF_CO2 ).value( U_W_M2 );
FCO2_record.set( runToDate, FCO2_current );
#define CO2_LAG 20
#define CO2_WINDOW 10
double FCO2_lagged_mean = 0.0; /* window mean of FCO2_past */
if( runToDate > core->getStartDate() + CO2_LAG ) {
for( int i=runToDate-CO2_LAG-CO2_WINDOW; i<runToDate-CO2_LAG; i++ ) {
FCO2_lagged_mean += FCO2_record.get( i );
}
FCO2_lagged_mean /= CO2_WINDOW;
}
// Subtract off the current FCO2 computed by ForcingComponent, and use the lagged
// window mean computed above instead
Ftot = Ftot - FCO2_current + FCO2_lagged_mean;
// At this point, adjust out internal Ftot tracking variable by the delta between
// Ftot and last_Ftot
internal_Ftot += ( Ftot - last_Ftot );
last_Ftot = Ftot;
// S is the equilibrium climate sensitivity for a doubling of CO2
// this corresponds, according to Knutti and Hegerl (2008), to ~3.7 W/m2
// Use this to compute a temperature resulting from total forcing
tgaveq = S / 3.7 * internal_Ftot;
const double heatflux = core->sendMessage( M_GETDATA, D_HEAT_FLUX ).value( U_W_M2 );
tgav = S / 3.7 * ( internal_Ftot - heatflux );
}
H_LOG( logger, Logger::DEBUG ) << " internal_Ftot=" << internal_Ftot << " tgav=" << tgav << " in " << runToDate << std::endl;
}
//------------------------------------------------------------------------------
// documentation is inherited
unitval TemperatureComponent::getData( const std::string& varName,
const double date ) throw ( h_exception ) {
unitval returnval;
H_ASSERT( date == Core::undefinedIndex(), "Only current temperatures provided" );
if( varName == D_GLOBAL_TEMP ) {
returnval = tgav;
} else if( varName == D_GLOBAL_TEMPEQ ) {
returnval = tgaveq;
} else {
H_THROW( "Caller requested unknown datum: " + varName );
}
return returnval;
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::shutDown() {
H_LOG( logger, Logger::DEBUG ) << "goodbye " << getComponentName() << std::endl;
logger.close();
}
//------------------------------------------------------------------------------
// documentation is inherited
void TemperatureComponent::accept( AVisitor* visitor ) {
visitor->visit( this );
}
}