Line data Source code
1 : /*******************************************************************************
2 : * ALMA - Atacama Large Millimiter Array
3 : * (c) CSIC - Instituto de Física Fundamental (Spain), 2024
4 : *
5 : * This library is free software; you can redistribute it and/or
6 : * modify it under the terms of the GNU Lesser General Public
7 : * License as published by the Free Software Foundation; either
8 : * version 2.1 of the License, or (at your option) any later version.
9 : *
10 : * This library is distributed in the hope that it will be useful,
11 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 : * Lesser General Public License for more details.
14 : *
15 : * You should have received a copy of the GNU Lesser General Public
16 : * License along with this library; if not, write to the Free Software
17 : * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 : *
19 : * "@(#) $Id: ATMProfile.cpp Exp $"
20 : *
21 : * who when what
22 : * -------- -------- ----------------------------------------------
23 : * pardo 24/03/09 created
24 : */
25 :
26 : #include "ATMProfile.h"
27 : #include "ATMException.h"
28 : #include <iostream>
29 : #include <math.h>
30 : #include <sstream>
31 :
32 :
33 :
34 : ATM_NAMESPACE_BEGIN
35 :
36 0 : AtmProfile::AtmProfile(unsigned int n)
37 : {
38 0 : numLayer_ = n;
39 0 : initBasicAtmosphericParameterThresholds();
40 0 : for(unsigned int i = 0; i < numLayer_; ++i) {
41 0 : v_layerO3_.push_back(0.0);
42 0 : v_layerCO_.push_back(0.0);
43 0 : v_layerN2O_.push_back(0.0);
44 0 : v_layerNO2_.push_back(0.0);
45 0 : v_layerSO2_.push_back(0.0);
46 0 : v_layerHCl_.push_back(0.0);
47 0 : v_layerHCN_.push_back(0.0);
48 0 : v_layerThickness_.push_back(0.0);
49 0 : v_layerTemperature_.push_back(0.0);
50 0 : v_layerTemperature0_.push_back(0.0);
51 0 : v_layerTemperature1_.push_back(0.0);
52 0 : v_layerPressure_.push_back(0.0);
53 0 : v_layerPressure0_.push_back(0.0);
54 0 : v_layerPressure1_.push_back(0.0);
55 0 : v_layerWaterVapor_.push_back(0.0);
56 0 : v_layerWaterVapor0_.push_back(0.0);
57 0 : v_layerWaterVapor1_.push_back(0.0);
58 : }
59 0 : }
60 :
61 :
62 :
63 5 : AtmProfile::AtmProfile(const Length &altitude,
64 : const Pressure &groundPressure,
65 : const Temperature &groundTemperature,
66 : double tropoLapseRate,
67 : const Humidity &relativeHumidity,
68 : const Length &wvScaleHeight,
69 : const Pressure &pressureStep,
70 : double pressureStepFactor,
71 : const Length &topAtmProfile,
72 5 : unsigned int atmType) :
73 : // Atmospheretype atmType):
74 5 : typeAtm_(atmType), groundTemperature_(groundTemperature),
75 5 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
76 5 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
77 5 : pressureStep_(pressureStep), pressureStepFactor_(1.075), // pressureStepFactor_(pressureStepFactor), 05JUN2024
78 10 : altitude_(altitude), topAtmProfile_(topAtmProfile)
79 : {
80 :
81 5 : if(pressureStep_.get("mb")<5.0 || pressureStep_.get("mb")>12.5){pressureStep_=Pressure(10,"mb");}
82 : // std::cout << "Warning pressureStep out of limits. Fixed to 10 mb" << std::endl; } //pressureStep has to be betwee 5 and 12.5 mb 05JUN2024
83 5 : numLayer_ = 0;
84 5 : numLayer_ = mkAtmProfile();
85 : // std::cout << "ad numLayer_=" << numLayer_ << std::endl;
86 5 : initBasicAtmosphericParameterThresholds();
87 5 : newBasicParam_ = true;
88 :
89 :
90 5 : }
91 :
92 0 : AtmProfile::AtmProfile(const Length &altitude,
93 : const Pressure &groundPressure,
94 : const Temperature &groundTemperature,
95 : double tropoLapseRate,
96 : const Humidity &relativeHumidity,
97 : const Length &wvScaleHeight,
98 : const Pressure &pressureStep,
99 : double pressureStepFactor,
100 : const Length &topAtmProfile,
101 : unsigned int atmType,
102 : const vector<Length> &v_layerBoundaries,
103 0 : const vector<Temperature> &v_layerTemperature):
104 0 : typeAtm_(atmType), groundTemperature_(groundTemperature),
105 0 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
106 0 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
107 0 : pressureStep_(pressureStep), pressureStepFactor_(1.075), // pressureStepFactor_(pressureStepFactor), 05JUN2024
108 0 : altitude_(altitude), topAtmProfile_(topAtmProfile)
109 : {
110 0 : if(pressureStep_.get("mb")<5.0 || pressureStep_.get("mb")>12.5){pressureStep_=Pressure(10,"mb");}
111 : // std::cout << "Warning pressureStep out of limits. Fixed to 10 mb" << std::endl; } //pressureStep has to be between 5 and 12.5 mb 05JUN2024
112 0 : numLayer_ = 0;
113 0 : numLayer_ = mkAtmProfile();
114 0 : initBasicAtmosphericParameterThresholds();
115 0 : newBasicParam_ = true;
116 0 : unsigned int nL1 = v_layerBoundaries.size();
117 0 : unsigned int nL2 = v_layerTemperature.size();
118 0 : if(nL1 == nL2 ) {
119 0 : double h=altitude_.get(Length::UnitMeter);
120 : double h0;
121 : double h1;
122 : double counter;
123 : double avT;
124 0 : for(unsigned int n = 0; n < numLayer_; n++) {
125 0 : h0=h;
126 0 : h1 = h0 + v_layerThickness_[n];
127 0 : counter = 0.0;
128 0 : avT = 0.0;
129 0 : for(unsigned int m = 0; m < nL1; m++) {
130 0 : if( h0 <= v_layerBoundaries[m].get(Length::UnitMeter) && h1 >= v_layerBoundaries[m].get(Length::UnitMeter) ){
131 : // std::cout << "n=" << n << " h0=" << h0 << " h1=" << h1 << " v_layerBoundaries[" << m << "]=" << v_layerBoundaries[m].get(Length::UnitMeter) << std::endl;
132 0 : avT = avT + v_layerTemperature[m].get(Temperature::UnitKelvin);
133 0 : counter = counter + 1.0;
134 : }
135 : }
136 0 : if(avT > 0.0){
137 : // std::cout << "layer" << n << "old average temperature:" << v_layerTemperature_[n] << std::endl;
138 0 : v_layerTemperature_[n] = avT/counter;
139 : // std::cout << "layer" << n << "new average temperature:" << v_layerTemperature_[n] << std::endl;
140 : }
141 0 : h=h1;
142 : }
143 : }
144 0 : }
145 :
146 0 : AtmProfile::AtmProfile(const Length &altitude,
147 : const Pressure &groundPressure,
148 : const Temperature &groundTemperature,
149 : double tropoLapseRate,
150 : const Humidity &relativeHumidity,
151 : const Length &wvScaleHeight,
152 0 : unsigned int atmType) :
153 : // Atmospheretype atmType):
154 0 : typeAtm_(atmType), groundTemperature_(groundTemperature),
155 0 : tropoLapseRate_(tropoLapseRate), groundPressure_(groundPressure),
156 0 : relativeHumidity_(relativeHumidity), wvScaleHeight_(wvScaleHeight),
157 0 : pressureStep_(10.0, Pressure::UnitMilliBar), pressureStepFactor_(1.075), altitude_(altitude), //pressureStepFactor fixed to 1.1 05JUN2024
158 0 : topAtmProfile_(48.0, Length::UnitKiloMeter)
159 : {
160 0 : numLayer_ = 0;
161 0 : numLayer_ = mkAtmProfile();
162 0 : initBasicAtmosphericParameterThresholds();
163 0 : newBasicParam_ = true;
164 0 : }
165 :
166 0 : AtmProfile::AtmProfile(const Length &altitude,
167 : const vector<Length> &v_layerThickness,
168 : const vector<Pressure> &v_layerPressure,
169 : const vector<Temperature> &v_layerTemperature,
170 0 : const vector<MassDensity> &v_layerWaterVapor)
171 : {
172 0 : newBasicParam_ = true;
173 0 : unsigned int nL1 = v_layerThickness.size();
174 0 : unsigned int nL2 = v_layerPressure.size();
175 0 : unsigned int nL3 = v_layerTemperature.size();
176 0 : unsigned int nL4 = v_layerWaterVapor.size();
177 :
178 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4) {
179 0 : numLayer_ = nL1;
180 0 : for(unsigned int n = 0; n < numLayer_; n++) {
181 0 : v_layerO3_.push_back(0.0);
182 0 : v_layerCO_.push_back(0.0);
183 0 : v_layerN2O_.push_back(0.0);
184 0 : v_layerNO2_.push_back(0.0);
185 0 : v_layerSO2_.push_back(0.0);
186 0 : v_layerHCl_.push_back(0.0);
187 0 : v_layerHCN_.push_back(0.0);
188 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
189 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
190 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
191 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
192 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
193 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
194 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
195 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
196 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
197 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
198 : }
199 0 : } else {
200 0 : numLayer_ = 0;
201 : }
202 0 : initBasicAtmosphericParameterThresholds();
203 0 : }
204 :
205 0 : AtmProfile::AtmProfile(const vector<Length> &v_layerBoundaries,
206 : const vector<Pressure> &v_layerPressure,
207 : const vector<Temperature> &v_layerTemperature,
208 0 : const vector<MassDensity> &v_layerWaterVapor)
209 : {
210 0 : newBasicParam_ = true;
211 0 : unsigned int nL1 = v_layerBoundaries.size();
212 0 : unsigned int nL2 = v_layerPressure.size();
213 0 : unsigned int nL3 = v_layerTemperature.size();
214 0 : unsigned int nL4 = v_layerWaterVapor.size();
215 :
216 0 : if(nL1 == nL2 && nL2 == nL3 && nL3 == nL4) {
217 0 : numLayer_ = nL1 - 1;
218 0 : altitude_ = v_layerBoundaries[0];
219 0 : for(unsigned int n = 0; n < numLayer_; n++) {
220 0 : v_layerO3_.push_back(0.0);
221 0 : v_layerCO_.push_back(0.0);
222 0 : v_layerN2O_.push_back(0.0);
223 0 : v_layerNO2_.push_back(0.0);
224 0 : v_layerSO2_.push_back(0.0);
225 0 : v_layerHCl_.push_back(0.0);
226 0 : v_layerHCN_.push_back(0.0);
227 0 : v_layerThickness_.push_back(v_layerBoundaries[n+1].get(Length::UnitMeter) - v_layerBoundaries[n].get(Length::UnitMeter));
228 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
229 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
230 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
231 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
232 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
233 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
234 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
235 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
236 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
237 : }
238 0 : } else {
239 0 : numLayer_ = 0;
240 : }
241 0 : initBasicAtmosphericParameterThresholds();
242 0 : }
243 :
244 0 : AtmProfile::AtmProfile(const Length &altitude,
245 : const vector<Length> &v_layerThickness,
246 : const vector<Pressure> &v_layerPressure,
247 : const vector<Temperature> &v_layerTemperature,
248 0 : const vector<NumberDensity> &v_layerWaterVapor)
249 : {
250 0 : newBasicParam_ = true;
251 0 : unsigned int nL1 = v_layerThickness.size();
252 0 : unsigned int nL2 = v_layerPressure.size();
253 0 : unsigned int nL3 = v_layerTemperature.size();
254 0 : unsigned int nL4 = v_layerWaterVapor.size();
255 :
256 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4) {
257 0 : numLayer_ = nL1;
258 0 : for(unsigned int n = 0; n < numLayer_; n++) {
259 0 : v_layerO3_.push_back(0.0);
260 0 : v_layerCO_.push_back(0.0);
261 0 : v_layerN2O_.push_back(0.0);
262 0 : v_layerNO2_.push_back(0.0);
263 0 : v_layerSO2_.push_back(0.0);
264 0 : v_layerHCl_.push_back(0.0);
265 0 : v_layerHCN_.push_back(0.0);
266 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
267 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
268 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
269 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
270 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
271 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
272 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
273 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
274 0 : / (1000.0 * 6.023e23));
275 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
276 0 : / (1000.0 * 6.023e23));
277 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
278 0 : / (1000.0 * 6.023e23));
279 : }
280 0 : } else {
281 0 : numLayer_ = 0;
282 : }
283 0 : initBasicAtmosphericParameterThresholds();
284 0 : }
285 :
286 0 : AtmProfile::AtmProfile(const vector<Length> &v_layerBoundaries,
287 : const vector<Pressure> &v_layerPressure,
288 : const vector<Temperature> &v_layerTemperature,
289 0 : const vector<NumberDensity> &v_layerWaterVapor)
290 : {
291 0 : newBasicParam_ = true;
292 0 : unsigned int nL1 = v_layerBoundaries.size();
293 0 : unsigned int nL2 = v_layerPressure.size();
294 0 : unsigned int nL3 = v_layerTemperature.size();
295 0 : unsigned int nL4 = v_layerWaterVapor.size();
296 :
297 0 : if(nL1 == nL2 && nL2 == nL3 && nL3 == nL4) {
298 0 : numLayer_ = nL1 - 1;
299 0 : altitude_ = v_layerBoundaries[0];
300 0 : for(unsigned int n = 0; n < numLayer_; n++) {
301 0 : v_layerO3_.push_back(0.0);
302 0 : v_layerCO_.push_back(0.0);
303 0 : v_layerN2O_.push_back(0.0);
304 0 : v_layerNO2_.push_back(0.0);
305 0 : v_layerSO2_.push_back(0.0);
306 0 : v_layerHCl_.push_back(0.0);
307 0 : v_layerHCN_.push_back(0.0);
308 0 : v_layerThickness_.push_back(v_layerBoundaries[n+1].get(Length::UnitMeter) - v_layerBoundaries[n].get(Length::UnitMeter));
309 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
310 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
311 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
312 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
313 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
314 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
315 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
316 0 : / (1000.0 * 6.023e23));
317 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
318 0 : / (1000.0 * 6.023e23));
319 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
320 0 : / (1000.0 * 6.023e23));
321 : }
322 0 : } else {
323 0 : numLayer_ = 0;
324 : }
325 0 : initBasicAtmosphericParameterThresholds();
326 0 : }
327 :
328 0 : AtmProfile::AtmProfile(const Length &altitude,
329 : const vector<Length> &v_layerThickness,
330 : const vector<Pressure> &v_layerPressure,
331 : const vector<Temperature> &v_layerTemperature,
332 : const vector<MassDensity> &v_layerWaterVapor,
333 0 : const vector<NumberDensity> &v_layerO3)
334 : {
335 0 : newBasicParam_ = true;
336 0 : unsigned int nL1 = v_layerThickness.size();
337 0 : unsigned int nL2 = v_layerPressure.size();
338 0 : unsigned int nL3 = v_layerTemperature.size();
339 0 : unsigned int nL4 = v_layerWaterVapor.size();
340 0 : unsigned int nL5 = v_layerO3.size();
341 :
342 0 : if(nL1 * 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5) {
343 0 : numLayer_ = nL1;
344 0 : for(unsigned int n = 0; n < numLayer_; n++) {
345 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
346 0 : v_layerCO_.push_back(0.0);
347 0 : v_layerN2O_.push_back(0.0);
348 0 : v_layerNO2_.push_back(0.0);
349 0 : v_layerSO2_.push_back(0.0);
350 0 : v_layerHCl_.push_back(0.0);
351 0 : v_layerHCN_.push_back(0.0);
352 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
353 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
354 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
355 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
356 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
357 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
358 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
359 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
360 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
361 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
362 : }
363 0 : } else {
364 0 : numLayer_ = 0;
365 : }
366 0 : initBasicAtmosphericParameterThresholds();
367 0 : }
368 :
369 0 : AtmProfile::AtmProfile(const Length &altitude,
370 : const vector<Length> &v_layerThickness,
371 : const vector<Pressure> &v_layerPressure,
372 : const vector<Temperature> &v_layerTemperature,
373 : const vector<NumberDensity> &v_layerWaterVapor,
374 0 : const vector<NumberDensity> &v_layerO3)
375 : {
376 0 : newBasicParam_ = true;
377 0 : unsigned int nL1 = v_layerThickness.size();
378 0 : unsigned int nL2 = v_layerPressure.size();
379 0 : unsigned int nL3 = v_layerTemperature.size();
380 0 : unsigned int nL4 = v_layerWaterVapor.size();
381 0 : unsigned int nL5 = v_layerO3.size();
382 :
383 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5) {
384 0 : numLayer_ = nL1;
385 0 : for(unsigned int n = 0; n < numLayer_; n++) {
386 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
387 0 : v_layerCO_.push_back(0.0);
388 0 : v_layerN2O_.push_back(0.0);
389 0 : v_layerNO2_.push_back(0.0);
390 0 : v_layerSO2_.push_back(0.0);
391 0 : v_layerHCl_.push_back(0.0);
392 0 : v_layerHCN_.push_back(0.0);
393 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
394 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
395 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
396 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
397 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
398 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
399 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
400 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
401 0 : / (1000.0 * 6.023e23));
402 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
403 0 : / (1000.0 * 6.023e23));
404 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
405 0 : / (1000.0 * 6.023e23));
406 : }
407 0 : } else {
408 0 : numLayer_ = 0;
409 : }
410 0 : initBasicAtmosphericParameterThresholds();
411 0 : }
412 :
413 0 : AtmProfile::AtmProfile(const Length &altitude,
414 : const vector<Length> &v_layerThickness,
415 : const vector<Pressure> &v_layerPressure,
416 : const vector<Temperature> &v_layerTemperature,
417 : const vector<MassDensity> &v_layerWaterVapor,
418 : const vector<NumberDensity> &v_layerO3,
419 : const vector<NumberDensity> &v_layerCO,
420 : const vector<NumberDensity> &v_layerN2O,
421 : const vector<NumberDensity> &v_layerNO2,
422 : const vector<NumberDensity> &v_layerSO2,
423 : const vector<NumberDensity> &v_layerHCl,
424 0 : const vector<NumberDensity> &v_layerHCN)
425 : {
426 0 : newBasicParam_ = true;
427 0 : unsigned int nL1 = v_layerThickness.size();
428 0 : unsigned int nL2 = v_layerPressure.size();
429 0 : unsigned int nL3 = v_layerTemperature.size();
430 0 : unsigned int nL4 = v_layerWaterVapor.size();
431 0 : unsigned int nL5 = v_layerO3.size();
432 0 : unsigned int nL6 = v_layerCO.size();
433 0 : unsigned int nL7 = v_layerN2O.size();
434 0 : unsigned int nL8 = v_layerNO2.size();
435 0 : unsigned int nL9 = v_layerSO2.size();
436 0 : unsigned int nL10 = v_layerHCl.size();
437 0 : unsigned int nL11 = v_layerHCN.size();
438 :
439 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5 && nL5 == nL6 && nL6
440 0 : == nL7 && nL7 == nL8 && nL8 == nL9 && nL9 == nL10 && nL10 == nL11) {
441 0 : numLayer_ = nL1;
442 0 : for(unsigned int n = 0; n < numLayer_; n++) {
443 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
444 0 : v_layerCO_.push_back(v_layerCO[n].get(NumberDensity::UnitInverseCubicMeter));
445 0 : v_layerN2O_.push_back(v_layerN2O[n].get(NumberDensity::UnitInverseCubicMeter));
446 0 : v_layerNO2_.push_back(v_layerNO2[n].get(NumberDensity::UnitInverseCubicMeter));
447 0 : v_layerSO2_.push_back(v_layerSO2[n].get(NumberDensity::UnitInverseCubicMeter));
448 0 : v_layerHCl_.push_back(v_layerHCl[n].get(NumberDensity::UnitInverseCubicMeter));
449 0 : v_layerHCN_.push_back(v_layerHCN[n].get(NumberDensity::UnitInverseCubicMeter));
450 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
451 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
452 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
453 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
454 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
455 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
456 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
457 0 : v_layerWaterVapor_.push_back(exp((log(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter))+ log(v_layerWaterVapor[n + 1].get(MassDensity::UnitKiloGramPerCubicMeter)))/2.0));
458 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(MassDensity::UnitKiloGramPerCubicMeter));
459 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(MassDensity::UnitKiloGramPerCubicMeter));
460 : }
461 0 : } else {
462 0 : numLayer_ = 0;
463 : }
464 0 : initBasicAtmosphericParameterThresholds();
465 0 : }
466 :
467 0 : AtmProfile::AtmProfile(const Length &altitude,
468 : const vector<Length> &v_layerThickness,
469 : const vector<Pressure> &v_layerPressure,
470 : const vector<Temperature> &v_layerTemperature,
471 : const vector<NumberDensity> &v_layerWaterVapor,
472 : const vector<NumberDensity> &v_layerO3,
473 : const vector<NumberDensity> &v_layerCO,
474 : const vector<NumberDensity> &v_layerN2O,
475 : const vector<NumberDensity> &v_layerNO2,
476 : const vector<NumberDensity> &v_layerSO2,
477 : const vector<NumberDensity> &v_layerHCl,
478 0 : const vector<NumberDensity> &v_layerHCN)
479 : {
480 0 : newBasicParam_ = true;
481 0 : unsigned int nL1 = v_layerThickness.size();
482 0 : unsigned int nL2 = v_layerPressure.size();
483 0 : unsigned int nL3 = v_layerTemperature.size();
484 0 : unsigned int nL4 = v_layerWaterVapor.size();
485 0 : unsigned int nL5 = v_layerO3.size();
486 0 : unsigned int nL6 = v_layerCO.size();
487 0 : unsigned int nL7 = v_layerN2O.size();
488 0 : unsigned int nL8 = v_layerNO2.size();
489 0 : unsigned int nL9 = v_layerSO2.size();
490 0 : unsigned int nL10 = v_layerHCl.size();
491 0 : unsigned int nL11 = v_layerHCN.size();
492 :
493 0 : if(nL1 + 1 == nL2 && nL2 == nL3 && nL3 == nL4 && nL4 == nL5 && nL5 == nL6 && nL6
494 0 : == nL7 && nL7 == nL8 && nL8 == nL9 && nL9 == nL10 && nL10 == nL11) {
495 0 : numLayer_ = nL1;
496 0 : for(unsigned int n = 0; n < numLayer_; n++) {
497 0 : v_layerO3_.push_back(v_layerO3[n].get(NumberDensity::UnitInverseCubicMeter));
498 0 : v_layerCO_.push_back(v_layerCO[n].get(NumberDensity::UnitInverseCubicMeter));
499 0 : v_layerN2O_.push_back(v_layerN2O[n].get(NumberDensity::UnitInverseCubicMeter));
500 0 : v_layerNO2_.push_back(v_layerNO2[n].get(NumberDensity::UnitInverseCubicMeter));
501 0 : v_layerSO2_.push_back(v_layerSO2[n].get(NumberDensity::UnitInverseCubicMeter));
502 0 : v_layerHCl_.push_back(v_layerHCl[n].get(NumberDensity::UnitInverseCubicMeter));
503 0 : v_layerHCN_.push_back(v_layerHCN[n].get(NumberDensity::UnitInverseCubicMeter));
504 0 : v_layerThickness_.push_back(v_layerThickness[n].get(Length::UnitMeter));
505 0 : v_layerTemperature_.push_back((v_layerTemperature[n].get(Temperature::UnitKelvin) + v_layerTemperature[n+1].get(Temperature::UnitKelvin))/2.0);
506 0 : v_layerTemperature0_.push_back(v_layerTemperature[n].get(Temperature::UnitKelvin));
507 0 : v_layerTemperature1_.push_back(v_layerTemperature[n+1].get(Temperature::UnitKelvin));
508 0 : v_layerPressure_.push_back(exp((log(v_layerPressure[n].get(Pressure::UnitMilliBar))+ log(v_layerPressure[n + 1].get(Pressure::UnitMilliBar)))/2.0));
509 0 : v_layerPressure0_.push_back(v_layerPressure[n].get(Pressure::UnitMilliBar));
510 0 : v_layerPressure1_.push_back(v_layerPressure[n+1].get(Pressure::UnitMilliBar));
511 0 : v_layerWaterVapor_.push_back( (exp((log(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter))+ log(v_layerWaterVapor[n + 1].get(NumberDensity::UnitInverseCubicMeter)))/2.0))* 18.0
512 0 : / (1000.0 * 6.023e23));
513 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor[n].get(NumberDensity::UnitInverseCubicMeter) * 18.0
514 0 : / (1000.0 * 6.023e23));
515 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor[n+1].get(NumberDensity::UnitInverseCubicMeter) * 18.0
516 0 : / (1000.0 * 6.023e23));
517 : }
518 0 : } else {
519 0 : numLayer_ = 0;
520 : }
521 0 : initBasicAtmosphericParameterThresholds();
522 0 : }
523 :
524 202 : AtmProfile::AtmProfile(const AtmProfile &a)
525 : { //:AtmType(a.type_){
526 : // std::cout<<"AtmProfile copy constructor"<<endl; COMMENTED OUT BY JUAN MAY/16/2005
527 202 : typeAtm_ = a.typeAtm_;
528 202 : groundTemperature_ = a.groundTemperature_;
529 202 : tropoLapseRate_ = a.tropoLapseRate_;
530 202 : groundPressure_ = a.groundPressure_;
531 202 : relativeHumidity_ = a.relativeHumidity_;
532 202 : wvScaleHeight_ = a.wvScaleHeight_;
533 202 : pressureStep_ = a.pressureStep_;
534 202 : pressureStepFactor_ = a.pressureStepFactor_;
535 202 : altitude_ = a.altitude_;
536 202 : topAtmProfile_ = a.topAtmProfile_;
537 202 : numLayer_ = a.numLayer_;
538 202 : newBasicParam_ = a.newBasicParam_;
539 202 : v_layerThickness_.reserve(numLayer_);
540 202 : v_layerPressure_.reserve(numLayer_);
541 202 : v_layerPressure0_.reserve(numLayer_);
542 202 : v_layerPressure1_.reserve(numLayer_);
543 202 : v_layerTemperature_.reserve(numLayer_);
544 202 : v_layerTemperature0_.reserve(numLayer_);
545 202 : v_layerTemperature1_.reserve(numLayer_);
546 202 : v_layerWaterVapor_.reserve(numLayer_);
547 202 : v_layerWaterVapor0_.reserve(numLayer_);
548 202 : v_layerWaterVapor1_.reserve(numLayer_);
549 202 : v_layerCO_.reserve(numLayer_);
550 202 : v_layerO3_.reserve(numLayer_);
551 202 : v_layerN2O_.reserve(numLayer_);
552 202 : v_layerNO2_.reserve(numLayer_);
553 202 : v_layerSO2_.reserve(numLayer_);
554 202 : v_layerHCl_.reserve(numLayer_);
555 202 : v_layerHCN_.reserve(numLayer_);
556 : // std::cout << "numLayer_=" << numLayer_ << std::endl; COMMENTED OUT BY JUAN MAY/16/2005
557 8642 : for(unsigned int n = 0; n < numLayer_; n++) {
558 8440 : v_layerThickness_.push_back(a.v_layerThickness_[n]);
559 8440 : v_layerTemperature_.push_back(a.v_layerTemperature_[n]);
560 8440 : v_layerTemperature0_.push_back(a.v_layerTemperature0_[n]);
561 8440 : v_layerTemperature1_.push_back(a.v_layerTemperature1_[n]);
562 : //cout << "n=" << n << std::endl;
563 8440 : v_layerWaterVapor_.push_back(a.v_layerWaterVapor_[n]);
564 8440 : v_layerWaterVapor0_.push_back(a.v_layerWaterVapor0_[n]);
565 8440 : v_layerWaterVapor1_.push_back(a.v_layerWaterVapor1_[n]);
566 8440 : v_layerPressure_.push_back(a.v_layerPressure_[n]);
567 8440 : v_layerPressure0_.push_back(a.v_layerPressure1_[n]);
568 8440 : v_layerPressure1_.push_back(a.v_layerPressure1_[n]);
569 8440 : v_layerCO_.push_back(a.v_layerCO_[n]);
570 8440 : v_layerO3_.push_back(a.v_layerO3_[n]);
571 8440 : v_layerN2O_.push_back(a.v_layerN2O_[n]);
572 8440 : v_layerNO2_.push_back(a.v_layerNO2_[n]);
573 8440 : v_layerSO2_.push_back(a.v_layerSO2_[n]);
574 8440 : v_layerHCl_.push_back(a.v_layerHCl_[n]);
575 8440 : v_layerHCN_.push_back(a.v_layerHCN_[n]);
576 : }
577 202 : altitudeThreshold_ = a.altitudeThreshold_;
578 202 : groundPressureThreshold_ = a.groundPressureThreshold_;
579 202 : groundTemperatureThreshold_ = a.groundTemperatureThreshold_;
580 202 : tropoLapseRateThreshold_ = a.tropoLapseRateThreshold_;
581 202 : relativeHumidityThreshold_ = a.relativeHumidityThreshold_;
582 202 : wvScaleHeightThreshold_ = a.wvScaleHeightThreshold_;
583 202 : }
584 0 : void AtmProfile::setBasicAtmosphericParameterThresholds(const Length &altitudeThreshold,
585 : const Pressure &groundPressureThreshold,
586 : const Temperature &groundTemperatureThreshold,
587 : double tropoLapseRateThreshold,
588 : const Humidity &relativeHumidityThreshold,
589 : const Length &wvScaleHeightThreshold)
590 : {
591 0 : altitudeThreshold_ = altitudeThreshold;
592 0 : groundPressureThreshold_ = groundPressureThreshold;
593 0 : groundTemperatureThreshold_ = groundTemperatureThreshold;
594 0 : tropoLapseRateThreshold_ = tropoLapseRateThreshold;
595 0 : relativeHumidityThreshold_ = relativeHumidityThreshold;
596 0 : wvScaleHeightThreshold_ =wvScaleHeightThreshold;
597 0 : }
598 :
599 5 : void AtmProfile::initBasicAtmosphericParameterThresholds()
600 : {
601 5 : altitudeThreshold_ = Length(1.0,Length::UnitMeter);
602 5 : groundPressureThreshold_ = Pressure(99.,Pressure::UnitPascal); // DB 2014-03-06 : Choose 99 Pascal instead of 100 Pascal because the threshold must be lower than the value of delta_pressure used in getAverageDispersiveDryPathLength_GroundPressureDerivative()
603 5 : groundTemperatureThreshold_ = Temperature(0.3,Temperature::UnitKelvin);
604 5 : tropoLapseRateThreshold_ = 0.01;
605 5 : relativeHumidityThreshold_ = Humidity(100.,Percent::UnitPercent);
606 5 : wvScaleHeightThreshold_ = Length(20.,Length::UnitMeter);
607 5 : }
608 :
609 120 : bool AtmProfile::updateAtmProfile(const Length &altitude,
610 : const Pressure &groundPressure,
611 : const Temperature &groundTemperature,
612 : double tropoLapseRate,
613 : const Humidity &relativeHumidity,
614 : const Length &wvScaleHeight)
615 : {
616 :
617 : /* TODO A faire: pour decider s'il faut recalculer le profile on devrait plutot donner des seuils, eg
618 : if(fabs(altitude_.get()-altitude.get())>0.1)mkNewProfile=true;
619 : */
620 :
621 : unsigned int numLayer;
622 120 : bool mkNewProfile = false;
623 :
624 : // if(altitude_.get() != altitude.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"altitude has changed" <<endl;
625 : // if(groundPressure_.get() != groundPressure.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"ground pressure has changed" <<endl;
626 : // if(groundTemperature_.get() != groundTemperature.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"ground temperature has changed"<<endl;
627 : // if(wvScaleHeight_.get() != wvScaleHeight.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"wv scale height has changed" <<endl;
628 : // if(tropoLapseRate_ != tropoLapseRate) mkNewProfile = true; //if(mkNewProfile)cout<<"tropo lapse rate has changed" <<endl;
629 : // if(relativeHumidity_.get() != relativeHumidity.get()) mkNewProfile = true; //if(mkNewProfile)cout<<"relative humidity has changed" <<endl;
630 : //! all these thresholds shoudl be parametrized.
631 120 : if(fabs(altitude_.get()-altitude.get())> altitudeThreshold_.get()) { // Length(0.1,Length::UnitMeter).get()) {
632 0 : mkNewProfile = true;
633 : //cout<<"altitude has changed" <<endl;
634 : }
635 120 : if(fabs(groundPressure_.get()-groundPressure.get())> groundPressureThreshold_.get()) { // Pressure(100.,Pressure::UnitPascal).get()) {
636 120 : mkNewProfile = true;
637 : //cout<<"ground pressure has changed" <<endl;
638 : }
639 120 : if(fabs(groundTemperature_.get()-groundTemperature.get()) > groundTemperatureThreshold_.get()) { // Temperature(0.3,Temperature::UnitKelvin).get()) {
640 0 : mkNewProfile = true;
641 : //cout<<"ground temperature has changed"<<endl;
642 : }
643 120 : if(fabs(wvScaleHeight_.get()-wvScaleHeight.get()) > wvScaleHeightThreshold_.get()) { // Length(20.,Length::UnitMeter).get() ) {
644 0 : mkNewProfile = true;
645 : //cout<<"wv scale height has changed" <<endl;
646 : }
647 120 : if(fabs(tropoLapseRate_-tropoLapseRate) > tropoLapseRateThreshold_) { // 0.01) {
648 0 : mkNewProfile = true;
649 : //cout<<"tropo lapse rate has changed" <<endl;
650 : }
651 : // we use WVR...
652 120 : if(fabs(relativeHumidity_.get()-relativeHumidity.get())> relativeHumidityThreshold_.get()) { // Humidity(100.,Percent::UnitPercent).get()) {
653 0 : mkNewProfile = true;
654 : //cout<<"relative humidity has changed" <<endl;
655 : }
656 120 : if(mkNewProfile) {
657 120 : newBasicParam_ = true;
658 120 : altitude_ = altitude;
659 120 : groundPressure_ = groundPressure;
660 120 : groundTemperature_ = groundTemperature;
661 120 : tropoLapseRate_ = tropoLapseRate;
662 120 : relativeHumidity_ = relativeHumidity;
663 120 : wvScaleHeight_ = wvScaleHeight;
664 120 : numLayer = mkAtmProfile();
665 120 : numLayer_ = numLayer;
666 : // std::cout << "There are new basic parameters, with " << numLayer_ << " layers " << std::endl;
667 : } else {
668 0 : numLayer = getNumLayer();
669 0 : numLayer_ = numLayer;
670 : }
671 :
672 120 : return mkNewProfile;
673 : }
674 :
675 : // Note that this setBasicAtmosphericParameters will be overrided by the subclasses.
676 0 : bool AtmProfile::setBasicAtmosphericParameters(const Length &altitude,
677 : const Pressure &groundPressure,
678 : const Temperature &groundTemperature,
679 : double tropoLapseRate,
680 : const Humidity &relativeHumidity,
681 : const Length &wvScaleHeight)
682 : {
683 0 : return updateAtmProfile(altitude,
684 : groundPressure,
685 : groundTemperature,
686 : tropoLapseRate,
687 : relativeHumidity,
688 0 : wvScaleHeight);
689 : }
690 :
691 3 : string AtmProfile::getAtmosphereType() const
692 : {
693 3 : return getAtmosphereType(typeAtm_);
694 : }
695 :
696 3 : string AtmProfile::getAtmosphereType(unsigned int typeAtm)
697 : {
698 3 : string type;
699 : string typeNames[] = { "TROPICAL", "MIDLATSUMMER", "MIDLATWINTER",
700 24 : "SUBARTSUMMER", "SUBARTWINTER","US_ST76"};
701 :
702 3 : if(typeAtm < typeATM_end) {
703 3 : type = typeNames[typeAtm-1];
704 : } else {
705 0 : type = "DEFAULT";
706 : }
707 6 : return type;
708 21 : }
709 :
710 1 : vector<Temperature> AtmProfile::getTemperatureProfile() const
711 : {
712 1 : vector<Temperature> t;
713 1 : t.reserve(v_layerTemperature_.size());
714 43 : for(unsigned int i = 0; i < v_layerTemperature_.size(); i++) {
715 42 : Temperature tt(v_layerTemperature_[i], Temperature::UnitKelvin);
716 42 : t.push_back(tt);
717 42 : }
718 1 : return t;
719 0 : }
720 :
721 7452 : Temperature AtmProfile::getLayerTemperature(unsigned int i) const
722 : {
723 : /*if(i > v_layerTemperature_.size() - 1) {
724 : Temperature t(-999.0, Temperature::UnitKelvin);
725 : return t;
726 : } else {
727 : Temperature t(v_layerTemperature_[i], Temperature::UnitKelvin);
728 : return t;
729 : }*/
730 7452 : if(i > v_layerTemperature_.size() - 1) {
731 0 : std::ostringstream oss;
732 0 : oss << "Not a valid layer: " << i;
733 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
734 0 : }
735 7452 : return Temperature(v_layerTemperature_[i], Temperature::UnitKelvin);
736 : }
737 :
738 0 : Temperature AtmProfile::getLayerBottomTemperature(unsigned int i) const
739 : {
740 0 : if(i > v_layerTemperature0_.size() - 1) {
741 0 : std::ostringstream oss;
742 0 : oss << "Not a valid layer: " << i;
743 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
744 0 : }
745 0 : return Temperature(v_layerTemperature0_[i], Temperature::UnitKelvin);
746 : }
747 :
748 0 : Temperature AtmProfile::getLayerTopTemperature(unsigned int i) const
749 : {
750 0 : if(i > v_layerTemperature1_.size() - 1) {
751 0 : std::ostringstream oss;
752 0 : oss << "Not a valid layer: " << i;
753 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
754 0 : }
755 0 : return Temperature(v_layerTemperature1_[i], Temperature::UnitKelvin);
756 : }
757 :
758 :
759 0 : void AtmProfile::setAltitude(const Length &groundaltitude)
760 : {
761 :
762 0 : if (groundaltitude <= altitude_){
763 :
764 : // std::cout << "extrapolar para abajo" << std::endl;
765 :
766 0 : int nextralayers = int( 0.50001+(altitude_-groundaltitude).get(Length::UnitMeter)/v_layerThickness_[0]);
767 0 : if (nextralayers == 0) {nextralayers=1;}
768 : // std::cout << "aaa=" << ( 0.50001+(altitude_-groundaltitude).get(Length::UnitMeter)/v_layerThickness_[0]) << std::endl;
769 0 : double newThickness = (altitude_-groundaltitude).get(Length::UnitMeter) / nextralayers;
770 :
771 : // std::cout << "Number of extra layers: " << nextralayers << " of thickness = " << newThickness << " m" << std::endl;
772 : // std::cout << "nextralayers-1 " << nextralayers-1 << std::endl;
773 :
774 0 : for(int i=nextralayers; i>0; i--){
775 0 : v_layerThickness_.push_back(newThickness);
776 0 : v_layerTemperature_.push_back(v_layerTemperature_[0]);
777 0 : v_layerTemperature0_.push_back(v_layerTemperature0_[0]);
778 0 : v_layerTemperature1_.push_back(v_layerTemperature1_[0]);
779 0 : v_layerWaterVapor_.push_back(v_layerWaterVapor_[0]);
780 0 : v_layerWaterVapor0_.push_back(v_layerWaterVapor0_[0]);
781 0 : v_layerWaterVapor1_.push_back(v_layerWaterVapor1_[0]);
782 0 : v_layerPressure_.push_back(v_layerPressure_[0]);
783 0 : v_layerPressure0_.push_back(v_layerPressure0_[0]);
784 0 : v_layerPressure1_.push_back(v_layerPressure1_[0]);
785 0 : v_layerCO_.push_back(v_layerCO_[0]);
786 0 : v_layerO3_.push_back(v_layerO3_[0]);
787 0 : v_layerN2O_.push_back(v_layerN2O_[0]);
788 0 : v_layerNO2_.push_back(v_layerNO2_[0]);
789 0 : v_layerSO2_.push_back(v_layerSO2_[0]);
790 0 : v_layerHCl_.push_back(v_layerHCl_[0]);
791 0 : v_layerHCN_.push_back(v_layerHCN_[0]);
792 : }
793 :
794 0 : for(int i=v_layerThickness_.size()-1; i>nextralayers-1; i--){
795 0 : v_layerThickness_[i] = v_layerThickness_[i-nextralayers];
796 0 : v_layerTemperature_[i] = v_layerTemperature_[i-nextralayers];
797 0 : v_layerTemperature0_[i] = v_layerTemperature0_[i-nextralayers];
798 0 : v_layerTemperature1_[i] = v_layerTemperature1_[i-nextralayers];
799 0 : v_layerWaterVapor_[i] = v_layerWaterVapor_[i-nextralayers];
800 0 : v_layerWaterVapor0_[i] = v_layerWaterVapor0_[i-nextralayers];
801 0 : v_layerWaterVapor1_[i] = v_layerWaterVapor1_[i-nextralayers];
802 0 : v_layerPressure_[i] = v_layerPressure_[i-nextralayers];
803 0 : v_layerPressure0_[i] = v_layerPressure0_[i-nextralayers];
804 0 : v_layerPressure1_[i] = v_layerPressure1_[i-nextralayers];
805 0 : v_layerCO_[i] = v_layerCO_[i-nextralayers];
806 0 : v_layerO3_[i] = v_layerO3_[i-nextralayers];
807 0 : v_layerN2O_[i] = v_layerN2O_[i-nextralayers];
808 0 : v_layerNO2_[i] = v_layerNO2_[i-nextralayers];
809 0 : v_layerSO2_[i] = v_layerSO2_[i-nextralayers];
810 0 : v_layerHCl_[i] = v_layerHCl_[i-nextralayers];
811 0 : v_layerHCN_[i] = v_layerHCN_[i-nextralayers];
812 : }
813 :
814 : // std::cout << "nextralayers=" << nextralayers << std::endl;
815 :
816 0 : for(int i=nextralayers; i>0 ; i--){
817 :
818 0 : v_layerThickness_[i-1] = newThickness;
819 0 : v_layerTemperature1_[i-1] = v_layerTemperature0_[i];
820 0 : v_layerTemperature0_[i-1] = v_layerTemperature0_[i]-tropoLapseRate_*0.001*v_layerThickness_[i-1];
821 0 : v_layerTemperature_[i-1] = (v_layerTemperature0_[i]+v_layerTemperature1_[i])/2.0;
822 :
823 0 : v_layerPressure1_[i-1] = v_layerPressure0_[i];
824 0 : v_layerPressure0_[i-1] = v_layerPressure0_[i]* exp(-0.0341695 *
825 0 : (-newThickness)/ v_layerTemperature_[i-1] );
826 0 : v_layerPressure_[i-1] = exp((log(v_layerPressure0_[i-1])+ log(v_layerPressure1_[i-1]))/2.0);
827 :
828 0 : v_layerWaterVapor1_[i-1] = v_layerWaterVapor0_[i];
829 0 : v_layerWaterVapor0_[i-1] = v_layerWaterVapor0_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
830 0 : v_layerWaterVapor_[i-1] = exp((log(v_layerWaterVapor0_[i-1])+ log(v_layerWaterVapor1_[i-1]))/2.0);
831 :
832 0 : v_layerCO_[i-1] = v_layerCO_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
833 0 : v_layerO3_[i-1] = v_layerO3_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
834 0 : v_layerN2O_[i-1] = v_layerN2O_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
835 0 : v_layerNO2_[i-1] = v_layerNO2_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
836 0 : v_layerSO2_[i-1] = v_layerSO2_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
837 0 : v_layerHCl_[i-1] = v_layerHCl_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
838 0 : v_layerHCN_[i-1] = v_layerHCN_[i]*(v_layerPressure0_[i-1]/v_layerPressure0_[i]);
839 : }
840 :
841 0 : tropoLayer_ = tropoLayer_ + (v_layerThickness_.size()-numLayer_);
842 0 : numLayer_=v_layerThickness_.size();
843 :
844 : }else{
845 : // std::cout << "extrapolar para arriba" << std::endl;
846 0 : double extraheight = fabs(groundaltitude.get(Length::UnitMeter)-altitude_.get(Length::UnitMeter));
847 0 : double cumulheight = 0.0;
848 :
849 : // std::cout << "before v_layerTemperature_.size()=" << v_layerTemperature_.size() << std::endl;
850 :
851 0 : for(unsigned int i=0; i<numLayer_; i++){
852 0 : cumulheight = cumulheight + v_layerThickness_[i];
853 : // std::cout << i << " cumulheight = " << cumulheight << " extraheight= " << extraheight << std::endl;
854 0 : if (cumulheight>=extraheight) {
855 :
856 0 : double remainingheight = fabs(cumulheight-extraheight);
857 : // std::cout << "remainingheight = " << remainingheight << " m" << std::endl;
858 0 : v_layerThickness_.erase(v_layerThickness_.begin(),v_layerThickness_.begin()+i);
859 0 : v_layerThickness_[0] = remainingheight;
860 :
861 :
862 0 : v_layerTemperature0_.erase(v_layerTemperature0_.begin(),v_layerTemperature0_.begin()+i);
863 0 : v_layerTemperature1_.erase(v_layerTemperature1_.begin(),v_layerTemperature1_.begin()+i);
864 0 : v_layerTemperature_.erase(v_layerTemperature_.begin(),v_layerTemperature_.begin()+i);
865 0 : v_layerTemperature0_[0] = v_layerTemperature0_[1]-tropoLapseRate_*0.001*v_layerThickness_[0];
866 0 : v_layerTemperature_[0] = (v_layerTemperature0_[0]+v_layerTemperature1_[0])/2.0;
867 :
868 0 : v_layerPressure0_.erase(v_layerPressure0_.begin(),v_layerPressure0_.begin()+i);
869 : // std::cout << "v_layerPressure1_[0]=" << v_layerPressure1_[0] << " remainingheight=" << remainingheight << std::endl;
870 0 : v_layerPressure1_.erase(v_layerPressure1_.begin(),v_layerPressure1_.begin()+i);
871 0 : v_layerPressure_.erase(v_layerPressure_.begin(),v_layerPressure_.begin()+i);
872 0 : v_layerPressure0_[0] = v_layerPressure1_[0]* exp(-0.0341695 *
873 0 : (-remainingheight)/ v_layerTemperature_[0] );
874 0 : v_layerPressure_[0] = exp((log(v_layerPressure0_[0])+ log(v_layerPressure1_[0]))/2.0);
875 :
876 0 : v_layerWaterVapor0_.erase(v_layerWaterVapor0_.begin(),v_layerWaterVapor0_.begin()+i);
877 0 : v_layerWaterVapor1_.erase(v_layerWaterVapor1_.begin(),v_layerWaterVapor1_.begin()+i);
878 0 : v_layerWaterVapor_.erase(v_layerWaterVapor_.begin(),v_layerWaterVapor_.begin()+i);
879 0 : v_layerWaterVapor0_[0] = v_layerWaterVapor1_[0]*(v_layerPressure0_[0]/v_layerPressure1_[0]);
880 0 : v_layerWaterVapor_[0] = exp((log(v_layerWaterVapor0_[0])+ log(v_layerWaterVapor1_[0]))/2.0);
881 :
882 0 : v_layerCO_.erase(v_layerCO_.begin(),v_layerCO_.begin()+i);
883 0 : v_layerO3_.erase(v_layerO3_.begin(),v_layerO3_.begin()+i);
884 0 : v_layerN2O_.erase(v_layerN2O_.begin(),v_layerN2O_.begin()+i);
885 0 : v_layerNO2_.erase(v_layerNO2_.begin(),v_layerNO2_.begin()+i);
886 0 : v_layerSO2_.erase(v_layerSO2_.begin(),v_layerSO2_.begin()+i);
887 0 : v_layerHCl_.erase(v_layerHCl_.begin(),v_layerHCl_.begin()+i);
888 0 : v_layerHCN_.erase(v_layerHCN_.begin(),v_layerHCN_.begin()+i);
889 :
890 0 : tropoLayer_ = tropoLayer_ - i;
891 :
892 : // std::cout << "after v_layerTemperature_.size()=" << v_layerTemperature_.size() << std::endl;
893 0 : break;
894 : }
895 : }
896 :
897 0 : numLayer_ = v_layerThickness_.size();
898 :
899 : }
900 :
901 0 : altitude_ = groundaltitude;
902 0 : groundTemperature_ = v_layerTemperature0_[0];
903 :
904 : // std::cout << "v_layerPressure0_[0]=" << v_layerPressure0_[0] << std::endl;
905 :
906 0 : groundPressure_ = v_layerPressure0_[0]*100.0; // default units for Pressure are Pascals
907 :
908 : /*
909 : const Temperature Tdif(((true_antenna_altitude.get(Length::UnitKiloMeter)-altitude.get(Length::UnitKiloMeter))*tropoLapseRate),Temperature::UnitKelvin);
910 : groundTemperature_ = groundTemperature+Tdif;
911 : groundPressure_ = groundPressure * exp(-0.0341695 * pow((6.371/(6.371+altitude.get(Length::UnitKiloMeter))),2.0)
912 : * (true_antenna_altitude.get(Length::UnitMeter)-altitude.get(Length::UnitMeter))
913 : / ((groundTemperature.get(Temperature::UnitKelvin)+Tdif.get(Temperature::UnitKelvin))/2.0) );
914 : */
915 :
916 :
917 :
918 0 : }
919 :
920 0 : void AtmProfile::setLayerTemperature(unsigned int i, const Temperature &layerTemperature)
921 : {
922 0 : if(i < v_layerTemperature_.size()) {
923 0 : v_layerTemperature_[i] = layerTemperature.get(Temperature::UnitKelvin);
924 : }
925 0 : }
926 :
927 0 : vector<Length> AtmProfile::getThicknessProfile() const
928 : {
929 0 : vector<Length> l;
930 0 : l.reserve(v_layerThickness_.size());
931 0 : for(unsigned int i = 0; i < v_layerThickness_.size(); i++) {
932 0 : Length ll(v_layerThickness_[i], Length::UnitMeter);
933 0 : l.push_back(ll);
934 0 : }
935 0 : return l;
936 0 : }
937 :
938 4270 : Length AtmProfile::getLayerThickness(unsigned int i) const
939 : {
940 : /*if(i > v_layerThickness_.size() - 1) {
941 : Length l(-999.0, Length::UnitMeter);
942 : return l;
943 : } else {
944 : Length l(v_layerThickness_[i], Length::UnitMeter);
945 : return l;
946 : }*/
947 4270 : if(i > v_layerThickness_.size() - 1) {
948 0 : std::ostringstream oss;
949 0 : oss << "Not a valid layer: " << i;
950 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
951 0 : }
952 4270 : return Length(v_layerThickness_[i], Length::UnitMeter);
953 : }
954 :
955 0 : Length AtmProfile::getLayerBottomHeightAboveGround(unsigned int i) const
956 : {
957 0 : if(i > v_layerThickness_.size() - 1) {
958 0 : std::ostringstream oss;
959 0 : oss << "Not a valid layer: " << i;
960 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
961 0 : }
962 0 : double h=0.0;
963 0 : for(unsigned int j = 0; j < i; j++) {
964 0 : h = h + v_layerThickness_[j];
965 : }
966 0 : return Length(h, Length::UnitMeter);
967 : }
968 :
969 0 : Length AtmProfile::getLayerTopHeightAboveGround(unsigned int i) const
970 : {
971 0 : if(i > v_layerThickness_.size() - 1) {
972 0 : std::ostringstream oss;
973 0 : oss << "Not a valid layer: " << i;
974 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
975 0 : }
976 0 : double h=0.0;
977 0 : for(unsigned int j = 0; j < i+1; j++) {
978 0 : h = h + v_layerThickness_[j];
979 : }
980 0 : return Length(h, Length::UnitMeter);
981 : }
982 :
983 0 : Length AtmProfile::getLayerBottomHeightAboveSeaLevel(unsigned int i) const
984 : {
985 0 : if(i > v_layerThickness_.size() - 1) {
986 0 : std::ostringstream oss;
987 0 : oss << "Not a valid layer: " << i;
988 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
989 0 : }
990 0 : double h=altitude_.get(Length::UnitMeter);
991 0 : for(unsigned int j = 0; j < i; j++) {
992 0 : h = h + v_layerThickness_[j];
993 : }
994 0 : return Length(h, Length::UnitMeter);
995 : }
996 :
997 0 : Length AtmProfile::getLayerTopHeightAboveSeaLevel(unsigned int i) const
998 : {
999 0 : if(i > v_layerThickness_.size() - 1) {
1000 0 : std::ostringstream oss;
1001 0 : oss << "Not a valid layer: " << i;
1002 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1003 0 : }
1004 0 : double h=altitude_.get(Length::UnitMeter);
1005 0 : for(unsigned int j = 0; j < i+1; j++) {
1006 0 : h = h + v_layerThickness_[j];
1007 : }
1008 0 : return Length(h, Length::UnitMeter);
1009 : }
1010 :
1011 :
1012 0 : void AtmProfile::setLayerThickness(unsigned int i, const Length &layerThickness)
1013 : {
1014 0 : if(i < v_layerThickness_.size()) {
1015 0 : v_layerThickness_[i] = layerThickness.get(Length::UnitMeter);
1016 : }
1017 0 : }
1018 :
1019 112 : MassDensity AtmProfile::getLayerWaterVaporMassDensity(unsigned int i) const
1020 : {
1021 : /*if(i > v_layerWaterVapor_.size() - 1) {
1022 : MassDensity m(-999.0, MassDensity::UnitKiloGramPerCubicMeter);
1023 : return m;
1024 : } else {
1025 : MassDensity m(v_layerWaterVapor_[i], MassDensity::UnitKiloGramPerCubicMeter);
1026 : return m;
1027 : }*/
1028 112 : if(i > v_layerWaterVapor_.size() - 1) {
1029 0 : std::ostringstream oss;
1030 0 : oss << "Not a valid layer: " << i;
1031 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1032 0 : }
1033 112 : return MassDensity(v_layerWaterVapor_[i], MassDensity::UnitKiloGramPerCubicMeter);
1034 : }
1035 0 : MassDensity AtmProfile::getLayerBottomWaterVaporMassDensity(unsigned int i) const
1036 : {
1037 0 : if(i > v_layerWaterVapor0_.size() - 1) {
1038 0 : std::ostringstream oss;
1039 0 : oss << "Not a valid layer: " << i;
1040 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1041 0 : }
1042 0 : return MassDensity(v_layerWaterVapor0_[i], MassDensity::UnitKiloGramPerCubicMeter);
1043 : }
1044 0 : MassDensity AtmProfile::getLayerTopWaterVaporMassDensity(unsigned int i) const
1045 : {
1046 0 : if(i > v_layerWaterVapor1_.size() - 1) {
1047 0 : std::ostringstream oss;
1048 0 : oss << "Not a valid layer: " << i;
1049 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1050 0 : }
1051 0 : return MassDensity(v_layerWaterVapor1_[i], MassDensity::UnitKiloGramPerCubicMeter);
1052 : }
1053 :
1054 112 : NumberDensity AtmProfile::getLayerWaterVaporNumberDensity(unsigned int i) const
1055 : {
1056 : /*if(i > v_layerWaterVapor_.size() - 1) {
1057 : NumberDensity m(-999.0, NumberDensity::UnitInverseCubicMeter);
1058 : return m;
1059 : } else {
1060 : NumberDensity
1061 : m(v_layerWaterVapor_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1062 : return m;
1063 : }*/
1064 112 : if(i > v_layerWaterVapor_.size() - 1) {
1065 0 : std::ostringstream oss;
1066 0 : oss << "Not a valid layer: " << i;
1067 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1068 0 : }
1069 112 : return NumberDensity(v_layerWaterVapor_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1070 : }
1071 0 : NumberDensity AtmProfile::getLayerBottomWaterVaporNumberDensity(unsigned int i) const
1072 : {
1073 0 : if(i > v_layerWaterVapor0_.size() - 1) {
1074 0 : std::ostringstream oss;
1075 0 : oss << "Not a valid layer: " << i;
1076 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1077 0 : }
1078 0 : return NumberDensity(v_layerWaterVapor0_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1079 : }
1080 0 : NumberDensity AtmProfile::getLayerTopWaterVaporNumberDensity(unsigned int i) const
1081 : {
1082 0 : if(i > v_layerWaterVapor1_.size() - 1) {
1083 0 : std::ostringstream oss;
1084 0 : oss << "Not a valid layer: " << i;
1085 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1086 0 : }
1087 0 : return NumberDensity(v_layerWaterVapor1_[i] * 6.023e23 * 1000.0 / 18.0, NumberDensity::UnitInverseCubicMeter);
1088 : }
1089 :
1090 0 : void AtmProfile::setLayerWaterVaporMassDensity(unsigned int i, const MassDensity &layerWaterVapor)
1091 : {
1092 0 : if(i <= v_layerWaterVapor_.size() - 1) {
1093 0 : v_layerWaterVapor_[i] = layerWaterVapor.get(MassDensity::UnitKiloGramPerCubicMeter);
1094 : }
1095 0 : }
1096 :
1097 0 : void AtmProfile::setLayerWaterVaporNumberDensity(unsigned int i, const NumberDensity &layerWaterVapor)
1098 : {
1099 0 : if(i <= v_layerWaterVapor_.size() - 1) {
1100 0 : v_layerWaterVapor_[i] = layerWaterVapor.get(NumberDensity::UnitInverseCubicMeter) * 18.0 / (6.023e23 * 1000.0);
1101 : }
1102 0 : }
1103 :
1104 0 : vector<Pressure> AtmProfile::getPressureProfile() const
1105 : {
1106 0 : vector<Pressure> p;
1107 0 : p.reserve(v_layerPressure_.size());
1108 0 : for(unsigned int i = 0; i < v_layerPressure_.size(); i++) {
1109 0 : Pressure pp(v_layerPressure_[i], Pressure::UnitMilliBar);
1110 0 : p.push_back(pp);
1111 0 : }
1112 0 : return p;
1113 0 : }
1114 :
1115 3336 : Pressure AtmProfile::getLayerPressure(unsigned int i) const
1116 : {
1117 : /*if(i > v_layerPressure_.size() - 1) {
1118 : Pressure p(-999.0, Pressure::UnitMilliBar);
1119 : return p;
1120 : } else {
1121 : Pressure p(v_layerPressure_[i], Pressure::UnitMilliBar);
1122 : return p;
1123 : }*/
1124 3336 : if(i > v_layerPressure_.size() - 1) {
1125 0 : std::ostringstream oss;
1126 0 : oss << "Not a valid layer: " << i;
1127 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1128 0 : }
1129 3336 : return Pressure(v_layerPressure_[i], Pressure::UnitMilliBar);
1130 : }
1131 :
1132 0 : Pressure AtmProfile::getLayerBottomPressure(unsigned int i) const
1133 : {
1134 0 : if(i > v_layerPressure0_.size() - 1) {
1135 0 : std::ostringstream oss;
1136 0 : oss << "Not a valid layer: " << i;
1137 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1138 0 : }
1139 0 : return Pressure(v_layerPressure0_[i], Pressure::UnitMilliBar);
1140 : }
1141 :
1142 0 : Pressure AtmProfile::getLayerTopPressure(unsigned int i) const
1143 : {
1144 0 : if(i > v_layerPressure1_.size() - 1) {
1145 0 : std::ostringstream oss;
1146 0 : oss << "Not a valid layer: " << i;
1147 0 : throw AtmException(ATM_EXCEPTION_ARGS(oss.str().c_str()));
1148 0 : }
1149 0 : return Pressure(v_layerPressure1_[i], Pressure::UnitMilliBar);
1150 : }
1151 :
1152 :
1153 11507 : Length AtmProfile::getGroundWH2O() const
1154 : {
1155 11507 : double wm = 0;
1156 347535 : for(unsigned int j = 0; j < numLayer_; j++) {
1157 336028 : wm = wm + v_layerWaterVapor_[j] * v_layerThickness_[j]; // kg/m**2 or mm (from m*kg/m**3 IS units)
1158 : }
1159 11507 : wm = wm * 1e-3; // (pasar de mm a m)
1160 11507 : Length wh2o(wm);
1161 11507 : return wh2o;
1162 : }
1163 :
1164 125 : MassDensity AtmProfile::rwat(const Temperature &tt, const Humidity &rh, const Pressure &pp) const
1165 : {
1166 125 : double t = tt.get(Temperature::UnitKelvin);
1167 125 : double p = pp.get(Pressure::UnitMilliBar);
1168 125 : double u = rh.get(Percent::UnitPercent);
1169 : double e, es, rwat0;
1170 :
1171 125 : if(p <= 0 || t <= 0 || u <= 0) {
1172 0 : return MassDensity(0.0, MassDensity::UnitGramPerCubicMeter);
1173 : } else {
1174 125 : es = 6.105 * exp(25.22 / t * (t - 273.0) - 5.31 * log(t / 273.0));
1175 125 : e = 1.0 - (1.0 - u / 100.0) * es / p;
1176 125 : e = es * u / 100.0 / e;
1177 125 : rwat0 = e * 216.502 / t; //(en g/m*3)
1178 : }
1179 125 : return MassDensity(rwat0, MassDensity::UnitGramPerCubicMeter);
1180 : }
1181 :
1182 0 : Humidity AtmProfile::rwat_inv(const Temperature &tt, const MassDensity &dd, const Pressure &pp) const
1183 : {
1184 0 : double p = pp.get(Pressure::UnitMilliBar);
1185 0 : double t = tt.get(Temperature::UnitKelvin);
1186 0 : double r = dd.get(MassDensity::UnitGramPerCubicMeter);
1187 : double es, e, rinv;
1188 :
1189 0 : if(p <= 0 || t <= 0 || r <= 0) {
1190 0 : rinv = 0.0;
1191 : } else {
1192 0 : es = 6.105 * exp(25.22 / t * (t - 273.0) - 5.31 * log(t / 273.0));
1193 0 : e = r * t / 216.502;
1194 0 : rinv = 100 * (e * (p - es) / (es * (p - e)));
1195 0 : if(rinv < 0 && p < 3) rinv = 0.0;
1196 : }
1197 0 : return Humidity(rinv, Percent::UnitPercent);
1198 : }
1199 :
1200 3448 : vector<NumberDensity> AtmProfile::st76(const Length &h, unsigned int tip, unsigned int seasonal_diurnal) const
1201 : {
1202 : unsigned int i1, i2, i3, i_layer;
1203 : double x1, x2, x3, d;
1204 3448 : vector<NumberDensity> minorden;
1205 3448 : NumberDensity o3den, n2oden, coden, no2den, so2den, hclden, hcnden;
1206 : static const double avogad = 6.022045E+23;
1207 : static const double airmwt = 28.964;
1208 : // static const double h2omwt=18.015;
1209 3448 : double ha = h.get(Length::UnitKiloMeter);
1210 : // std::cout << "tip = " << tip << std::endl;
1211 : // if (tip > 10) {seasonal_diurnal=tip-10; tip=1;}else{seasonal_diurnal=0;}
1212 : // std::cout << "in st76 seasonal_diurnal,tip = " << seasonal_diurnal << "," << tip << std::endl;
1213 : static const double
1214 : alt[50] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
1215 : 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0,
1216 : 24.0, 25.0, 27.5, 30.0, 32.5, 35.0, 37.5, 40.0, 42.5, 45.0, 47.5,
1217 : 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0,
1218 : 105.0, 110.0, 115.0, 120.0 }; // REFERENCE LEVELS IN KM
1219 :
1220 : static const double
1221 : ozone[6][50] = { { 2.869E-02, 3.150E-02, 3.342E-02, 3.504E-02, 3.561E-02, 3.767E-02, 3.989E-02, 4.223E-02,
1222 : 4.471E-02, 5.000E-02, 5.595E-02, 6.613E-02, 7.815E-02, 9.289E-02, 1.050E-01, 1.256E-01,
1223 : 1.444E-01, 2.500E-01, 5.000E-01, 9.500E-01, 1.400E+00, 1.800E+00, 2.400E+00, 3.400E+00,
1224 : 4.300E+00, 5.400E+00, 7.800E+00, 9.300E+00, 9.850E+00, 9.700E+00, 8.800E+00, 7.500E+00,
1225 : 5.900E+00, 4.500E+00, 3.450E+00, 2.800E+00, 1.800E+00, 1.100E+00, 6.500E-01, 3.000E-01,
1226 : 1.800E-01, 3.300E-01, 5.000E-01, 5.200E-01, 5.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1227 : 5.000E-03, 5.000E-04 },
1228 : { 3.017E-02, 3.337E-02, 3.694E-02, 4.222E-02, 4.821E-02, 5.512E-02, 6.408E-02, 7.764E-02,
1229 : 9.126E-02, 1.111E-01, 1.304E-01, 1.793E-01, 2.230E-01, 3.000E-01, 4.400E-01, 5.000E-01,
1230 : 6.000E-01, 7.000E-01, 1.000E+00, 1.500E+00, 2.000E+00, 2.400E+00, 2.900E+00, 3.400E+00,
1231 : 4.000E+00, 4.800E+00, 6.000E+00, 7.000E+00, 8.100E+00, 8.900E+00, 8.700E+00, 7.550E+00,
1232 : 5.900E+00, 4.500E+00, 3.500E+00, 2.800E+00, 1.800E+00, 1.300E+00, 8.000E-01, 4.000E-01,
1233 : 1.900E-01, 2.000E-01, 5.700E-01, 7.500E-01, 7.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1234 : 5.000E-03, 5.000E-04 },
1235 : { 2.778E-02, 2.800E-02, 2.849E-02, 3.200E-02, 3.567E-02, 4.720E-02, 5.837E-02, 7.891E-02,
1236 : 1.039E-01, 1.567E-01, 2.370E-01, 3.624E-01, 5.232E-01, 7.036E-01, 8.000E-01, 9.000E-01,
1237 : 1.100E+00, 1.400E+00, 1.800E+00, 2.300E+00, 2.900E+00, 3.500E+00, 3.900E+00, 4.300E+00,
1238 : 4.700E+00, 5.100E+00, 5.600E+00, 6.100E+00, 6.800E+00, 7.100E+00, 7.200E+00, 6.900E+00,
1239 : 5.900E+00, 4.600E+00, 3.700E+00, 2.750E+00, 1.700E+00, 1.000E-00, 5.500E-01, 3.200E-01,
1240 : 2.500E-01, 2.300E-01, 5.500E-01, 8.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1241 : 5.000E-03, 5.000E-04 },
1242 : { 2.412E-02, 2.940E-02, 3.379E-02, 3.887E-02, 4.478E-02, 5.328E-02, 6.564E-02, 7.738E-02,
1243 : 9.114E-02, 1.420E-01, 1.890E-01, 3.050E-01, 4.100E-01, 5.000E-01, 6.000E-01, 7.000E-01,
1244 : 8.500E-01, 1.000E+00, 1.300E+00, 1.700E+00, 2.100E+00, 2.700E+00, 3.300E+00, 3.700E+00,
1245 : 4.200E+00, 4.500E+00, 5.300E+00, 5.700E+00, 6.900E+00, 7.700E+00, 7.800E+00, 7.000E+00,
1246 : 5.400E+00, 4.200E+00, 3.200E+00, 2.500E+00, 1.700E+00, 1.200E+00, 8.000E-01, 4.000E-01,
1247 : 2.000E-01, 1.800E-01, 6.500E-01, 9.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1248 : 5.000E-03, 5.000E-04 },
1249 : { 1.802E-02, 2.072E-02, 2.336E-02, 2.767E-02, 3.253E-02, 3.801E-02, 4.446E-02, 7.252E-02,
1250 : 1.040E-01, 2.100E-01, 3.000E-01, 3.500E-01, 4.000E-01, 6.500E-01, 9.000E-01, 1.200E+00,
1251 : 1.500E+00, 1.900E+00, 2.450E+00, 3.100E+00, 3.700E+00, 4.000E+00, 4.200E+00, 4.500E+00,
1252 : 4.600E+00, 4.700E+00, 4.900E+00, 5.400E+00, 5.900E+00, 6.200E+00, 6.250E+00, 5.900E+00,
1253 : 5.100E+00, 4.100E+00, 3.000E+00, 2.600E+00, 1.600E+00, 9.500E-01, 6.500E-01, 5.000E-01,
1254 : 3.300E-01, 1.300E-01, 7.500E-01, 8.000E-01, 8.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1255 : 5.000E-03, 5.000E-04 },
1256 : { 2.660E-02, 2.931E-02, 3.237E-02, 3.318E-02, 3.387E-02, 3.768E-02, 4.112E-02, 5.009E-02,
1257 : 5.966E-02, 9.168E-02, 1.313E-01, 2.149E-01, 3.095E-01, 3.846E-01, 5.030E-01, 6.505E-01,
1258 : 8.701E-01, 1.187E+00, 1.587E+00, 2.030E+00, 2.579E+00, 3.028E+00, 3.647E+00, 4.168E+00,
1259 : 4.627E+00, 5.118E+00, 5.803E+00, 6.553E+00, 7.373E+00, 7.837E+00, 7.800E+00, 7.300E+00,
1260 : 6.200E+00, 5.250E+00, 4.100E+00, 3.100E+00, 1.800E+00, 1.100E+00, 7.000E-01, 3.000E-01,
1261 : 2.500E-01, 3.000E-01, 5.000E-01, 7.000E-01, 7.000E-01, 4.000E-01, 2.000E-01, 5.000E-02,
1262 : 5.000E-03, 5.000E-04 } };
1263 :
1264 : static const double
1265 : n2o[6][50] = { { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1266 : 3.200E-01, 3.195E-01, 3.179E-01, 3.140E-01, 3.095E-01, 3.048E-01, 2.999E-01, 2.944E-01,
1267 : 2.877E-01, 2.783E-01, 2.671E-01, 2.527E-01, 2.365E-01, 2.194E-01, 2.051E-01, 1.967E-01,
1268 : 1.875E-01, 1.756E-01, 1.588E-01, 1.416E-01, 1.165E-01, 9.275E-02, 6.693E-02, 4.513E-02,
1269 : 2.751E-02, 1.591E-02, 9.378E-03, 4.752E-03, 3.000E-03, 2.065E-03, 1.507E-03, 1.149E-03,
1270 : 8.890E-04, 7.056E-04, 5.716E-04, 4.708E-04, 3.932E-04, 3.323E-04, 2.837E-04, 2.443E-04,
1271 : 2.120E-04, 1.851E-04 },
1272 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1273 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1274 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.607E-01, 1.323E-01, 1.146E-01, 1.035E-01,
1275 : 9.622E-02, 8.958E-02, 8.006E-02, 6.698E-02, 4.958E-02, 3.695E-02, 2.519E-02, 1.736E-02,
1276 : 1.158E-02, 7.665E-03, 5.321E-03, 3.215E-03, 2.030E-03, 1.397E-03, 1.020E-03, 7.772E-04,
1277 : 6.257E-04, 5.166E-04, 4.352E-04, 3.727E-04, 3.237E-04, 2.844E-04, 2.524E-04, 2.260E-04,
1278 : 2.039E-04, 1.851E-04 },
1279 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1280 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1281 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.621E-01, 1.362E-01, 1.230E-01, 1.124E-01,
1282 : 1.048E-01, 9.661E-02, 8.693E-02, 7.524E-02, 6.126E-02, 5.116E-02, 3.968E-02, 2.995E-02,
1283 : 2.080E-02, 1.311E-02, 8.071E-03, 4.164E-03, 2.629E-03, 1.809E-03, 1.321E-03, 1.007E-03,
1284 : 7.883E-04, 6.333E-04, 5.194E-04, 4.333E-04, 3.666E-04, 3.140E-04, 2.717E-04, 2.373E-04,
1285 : 2.089E-04, 1.851E-04 },
1286 : { 3.100E-01, 3.100E-01, 3.100E-01, 3.100E-01, 3.079E-01, 3.024E-01, 2.906E-01, 2.822E-01,
1287 : 2.759E-01, 2.703E-01, 2.651E-01, 2.600E-01, 2.549E-01, 2.494E-01, 2.433E-01, 2.355E-01,
1288 : 2.282E-01, 2.179E-01, 2.035E-01, 1.817E-01, 1.567E-01, 1.350E-01, 1.218E-01, 1.102E-01,
1289 : 9.893E-02, 8.775E-02, 7.327E-02, 5.941E-02, 4.154E-02, 3.032E-02, 1.949E-02, 1.274E-02,
1290 : 9.001E-03, 6.286E-03, 4.558E-03, 2.795E-03, 1.765E-03, 1.214E-03, 8.866E-04, 6.756E-04,
1291 : 5.538E-04, 4.649E-04, 3.979E-04, 3.459E-04, 3.047E-04, 2.713E-04, 2.439E-04, 2.210E-04,
1292 : 2.017E-04, 1.851E-04 },
1293 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1294 : 3.195E-01, 3.163E-01, 3.096E-01, 2.989E-01, 2.936E-01, 2.860E-01, 2.800E-01, 2.724E-01,
1295 : 2.611E-01, 2.421E-01, 2.174E-01, 1.843E-01, 1.621E-01, 1.362E-01, 1.230E-01, 1.122E-01,
1296 : 1.043E-01, 9.570E-02, 8.598E-02, 7.314E-02, 5.710E-02, 4.670E-02, 3.439E-02, 2.471E-02,
1297 : 1.631E-02, 1.066E-02, 7.064E-03, 3.972E-03, 2.508E-03, 1.726E-03, 1.260E-03, 9.602E-04,
1298 : 7.554E-04, 6.097E-04, 5.024E-04, 4.210E-04, 3.579E-04, 3.080E-04, 2.678E-04, 2.350E-04,
1299 : 2.079E-04, 1.851E-04 },
1300 : { 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01, 3.200E-01,
1301 : 3.200E-01, 3.195E-01, 3.179E-01, 3.140E-01, 3.095E-01, 3.048E-01, 2.999E-01, 2.944E-01,
1302 : 2.877E-01, 2.783E-01, 2.671E-01, 2.527E-01, 2.365E-01, 2.194E-01, 2.051E-01, 1.967E-01,
1303 : 1.875E-01, 1.756E-01, 1.588E-01, 1.416E-01, 1.165E-01, 9.275E-02, 6.693E-02, 4.513E-02,
1304 : 2.751E-02, 1.591E-02, 9.378E-03, 4.752E-03, 3.000E-03, 2.065E-03, 1.507E-03, 1.149E-03,
1305 : 8.890E-04, 7.056E-04, 5.716E-04, 4.708E-04, 3.932E-04, 3.323E-04, 2.837E-04, 2.443E-04,
1306 : 2.120E-04, 1.851E-04 } };
1307 :
1308 : static const double
1309 : co[6][50] = { { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1310 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1311 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1312 : 1.400E-02, 1.521E-02, 1.722E-02, 1.995E-02, 2.266E-02, 2.487E-02, 2.738E-02, 3.098E-02,
1313 : 3.510E-02, 3.987E-02, 4.482E-02, 5.092E-02, 5.985E-02, 6.960E-02, 9.188E-02, 1.938E-01,
1314 : 5.688E-01, 1.549E+00, 3.849E+00, 6.590E+00, 1.044E+01, 1.705E+01, 2.471E+01, 3.358E+01,
1315 : 4.148E+01, 5.000E+01 },
1316 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1317 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1318 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1319 : 1.400E-02, 1.521E-02, 1.722E-02, 1.995E-02, 2.266E-02, 2.487E-02, 2.716E-02, 2.962E-02,
1320 : 3.138E-02, 3.307E-02, 3.487E-02, 3.645E-02, 3.923E-02, 4.673E-02, 6.404E-02, 1.177E-01,
1321 : 2.935E-01, 6.815E-01, 1.465E+00, 2.849E+00, 5.166E+00, 1.008E+01, 1.865E+01, 2.863E+01,
1322 : 3.890E+01, 5.000E+01 },
1323 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1324 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1325 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1326 : 1.400E-02, 1.498E-02, 1.598E-02, 1.710E-02, 1.850E-02, 1.997E-02, 2.147E-02, 2.331E-02,
1327 : 2.622E-02, 3.057E-02, 3.803E-02, 6.245E-02, 1.480E-01, 2.926E-01, 5.586E-01, 1.078E+00,
1328 : 1.897E+00, 2.960E+00, 4.526E+00, 6.862E+00, 1.054E+01, 1.709E+01, 2.473E+01, 3.359E+01,
1329 : 4.149E+01, 5.000E+01 },
1330 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1331 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1332 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1333 : 1.400E-02, 1.510E-02, 1.649E-02, 1.808E-02, 1.997E-02, 2.183E-02, 2.343E-02, 2.496E-02,
1334 : 2.647E-02, 2.809E-02, 2.999E-02, 3.220E-02, 3.650E-02, 4.589E-02, 6.375E-02, 1.176E-01,
1335 : 3.033E-01, 7.894E-01, 1.823E+00, 3.402E+00, 5.916E+00, 1.043E+01, 1.881E+01, 2.869E+01,
1336 : 3.892E+01, 5.000E+01 },
1337 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1338 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1339 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1340 : 1.400E-02, 1.521E-02, 1.722E-02, 2.037E-02, 2.486E-02, 3.168E-02, 4.429E-02, 6.472E-02,
1341 : 1.041E-01, 1.507E-01, 2.163E-01, 3.141E-01, 4.842E-01, 7.147E-01, 1.067E+00, 1.516E+00,
1342 : 2.166E+00, 3.060E+00, 4.564E+00, 6.877E+00, 1.055E+01, 1.710E+01, 2.473E+01, 3.359E+01,
1343 : 4.149E+01, 5.000E+01 },
1344 : { 1.500E-01, 1.450E-01, 1.399E-01, 1.349E-01, 1.312E-01, 1.303E-01, 1.288E-01, 1.247E-01,
1345 : 1.185E-01, 1.094E-01, 9.962E-02, 8.964E-02, 7.814E-02, 6.374E-02, 5.025E-02, 3.941E-02,
1346 : 3.069E-02, 2.489E-02, 1.966E-02, 1.549E-02, 1.331E-02, 1.232E-02, 1.232E-02, 1.307E-02,
1347 : 1.400E-02, 1.498E-02, 1.598E-02, 1.710E-02, 1.850E-02, 2.009E-02, 2.220E-02, 2.497E-02,
1348 : 2.824E-02, 3.241E-02, 3.717E-02, 4.597E-02, 6.639E-02, 1.073E-01, 1.862E-01, 3.059E-01,
1349 : 6.375E-01, 1.497E+00, 3.239E+00, 5.843E+00, 1.013E+01, 1.692E+01, 2.467E+01, 3.356E+01,
1350 : 4.148E+01, 5.000E+01 } };
1351 :
1352 : static const double
1353 : den[6][50] = { { 2.450E+19, 2.231E+19, 2.028E+19, 1.827E+19, 1.656E+19, 1.499E+19, 1.353E+19, 1.218E+19,
1354 : 1.095E+19, 9.789E+18, 8.747E+18, 7.780E+18, 6.904E+18, 6.079E+18, 5.377E+18, 4.697E+18,
1355 : 4.084E+18, 3.486E+18, 2.877E+18, 2.381E+18, 1.981E+18, 1.651E+18, 1.381E+18, 1.169E+18,
1356 : 9.920E+17, 8.413E+17, 5.629E+17, 3.807E+17, 2.598E+17, 1.789E+17, 1.243E+17, 8.703E+16,
1357 : 6.147E+16, 4.352E+16, 3.119E+16, 2.291E+16, 1.255E+16, 6.844E+15, 3.716E+15, 1.920E+15,
1358 : 9.338E+14, 4.314E+14, 1.801E+14, 7.043E+13, 2.706E+13, 1.098E+13, 4.445E+12, 1.941E+12,
1359 : 8.706E+11, 4.225E+11 },
1360 : { .496E+19, 2.257E+19, 2.038E+19, 1.843E+19, 1.666E+19, 1.503E+19, 1.351E+19, 1.212E+19,
1361 : 1.086E+19, 9.716E+18, 8.656E+18, 7.698E+18, 6.814E+18, 6.012E+18, 5.141E+18, 4.368E+18,
1362 : 3.730E+18, 3.192E+18, 2.715E+18, 2.312E+18, 1.967E+18, 1.677E+18, 1.429E+18, 1.223E+18,
1363 : 1.042E+18, 8.919E+17, 6.050E+17, 4.094E+17, 2.820E+17, 1.927E+17, 1.338E+17, 9.373E+16,
1364 : 6.624E+16, 4.726E+16, 3.398E+16, 2.500E+16, 1.386E+16, 7.668E+15, 4.196E+15, 2.227E+15,
1365 : 1.109E+15, 4.996E+14, 1.967E+14, 7.204E+13, 2.541E+13, 9.816E+12, 3.816E+12, 1.688E+12,
1366 : 8.145E+11, 4.330E+11 },
1367 : { .711E+19, 2.420E+19, 2.158E+19, 1.922E+19, 1.724E+19, 1.542E+19, 1.376E+19, 1.225E+19,
1368 : 1.086E+19, 9.612E+18, 8.472E+18, 7.271E+18, 6.237E+18, 5.351E+18, 4.588E+18, 3.931E+18,
1369 : 3.368E+18, 2.886E+18, 2.473E+18, 2.115E+18, 1.809E+18, 1.543E+18, 1.317E+18, 1.125E+18,
1370 : 9.633E+17, 8.218E+17, 5.536E+17, 3.701E+17, 2.486E+17, 1.647E+17, 1.108E+17, 7.540E+16,
1371 : 5.202E+16, 3.617E+16, 2.570E+16, 1.863E+16, 1.007E+16, 5.433E+15, 2.858E+15, 1.477E+15,
1372 : 7.301E+14, 3.553E+14, 1.654E+14, 7.194E+13, 3.052E+13, 1.351E+13, 6.114E+12, 2.952E+12,
1373 : 1.479E+12, 7.836E+11 },
1374 : { .549E+19, 2.305E+19, 2.080E+19, 1.873E+19, 1.682E+19, 1.508E+19, 1.357E+19, 1.216E+19,
1375 : 1.088E+19, 9.701E+18, 8.616E+18, 7.402E+18, 6.363E+18, 5.471E+18, 4.699E+18, 4.055E+18,
1376 : 3.476E+18, 2.987E+18, 2.568E+18, 2.208E+18, 1.899E+18, 1.632E+18, 1.403E+18, 1.207E+18,
1377 : 1.033E+18, 8.834E+17, 6.034E+17, 4.131E+17, 2.839E+17, 1.938E+17, 1.344E+17, 9.402E+16,
1378 : 6.670E+16, 4.821E+16, 3.516E+16, 2.581E+16, 1.421E+16, 7.946E+15, 4.445E+15, 2.376E+15,
1379 : 1.198E+15, 5.311E+14, 2.022E+14, 7.221E+13, 2.484E+13, 9.441E+12, 3.624E+12, 1.610E+12,
1380 : 7.951E+11, 4.311E+11 },
1381 : { .855E+19, 2.484E+19, 2.202E+19, 1.950E+19, 1.736E+19, 1.552E+19, 1.383E+19, 1.229E+19,
1382 : 1.087E+19, 9.440E+18, 8.069E+18, 6.898E+18, 5.893E+18, 5.039E+18, 4.308E+18, 3.681E+18,
1383 : 3.156E+18, 2.704E+18, 2.316E+18, 1.982E+18, 1.697E+18, 1.451E+18, 1.241E+18, 1.061E+18,
1384 : 9.065E+17, 7.742E+17, 5.134E+17, 3.423E+17, 2.292E+17, 1.533E+17, 1.025E+17, 6.927E+16,
1385 : 4.726E+16, 3.266E+16, 2.261E+16, 1.599E+16, 8.364E+15, 4.478E+15, 2.305E+15, 1.181E+15,
1386 : 6.176E+14, 3.127E+14, 1.531E+14, 7.244E+13, 3.116E+13, 1.403E+13, 6.412E+12, 3.099E+12,
1387 : 1.507E+12, 7.814E+11 },
1388 : { .548E+19, 2.313E+19, 2.094E+19, 1.891E+19, 1.704E+19, 1.532E+19, 1.373E+19, 1.228E+19,
1389 : 1.094E+19, 9.719E+18, 8.602E+18, 7.589E+18, 6.489E+18, 5.546E+18, 4.739E+18, 4.050E+18,
1390 : 3.462E+18, 2.960E+18, 2.530E+18, 2.163E+18, 1.849E+18, 1.575E+18, 1.342E+18, 1.144E+18,
1391 : 9.765E+17, 8.337E+17, 5.640E+17, 3.830E+17, 2.524E+17, 1.761E+17, 1.238E+17, 8.310E+16,
1392 : 5.803E+16, 4.090E+16, 2.920E+16, 2.136E+16, 1.181E+16, 6.426E+15, 3.386E+15, 1.723E+15,
1393 : 8.347E+14, 3.832E+14, 1.711E+14, 7.136E+13, 2.924E+13, 1.189E+13, 5.033E+12, 2.144E+12,
1394 : 9.688E+11, 5.114E+11 } };
1395 :
1396 : static const double so2[50] = { 3.00E-04, 2.74E-04, 2.36E-04, 1.90E-04, 1.46E-04,
1397 : 1.18E-04, 9.71E-05, 8.30E-05, 7.21E-05, 6.56E-05,
1398 : 6.08E-05, 5.79E-05, 5.60E-05, 5.59E-05, 5.64E-05,
1399 : 5.75E-05, 5.75E-05, 5.37E-05, 4.78E-05, 3.97E-05,
1400 : 3.19E-05, 2.67E-05, 2.28E-05, 2.07E-05, 1.90E-05,
1401 : 1.75E-05, 1.54E-05, 1.34E-05, 1.21E-05, 1.16E-05,
1402 : 1.21E-05, 1.36E-05, 1.65E-05, 2.10E-05, 2.77E-05,
1403 : 3.56E-05, 4.59E-05, 5.15E-05, 5.11E-05, 4.32E-05,
1404 : 2.83E-05, 1.33E-05, 5.56E-06, 2.24E-06, 8.96E-07,
1405 : 3.58E-07, 1.43E-07, 5.73E-08, 2.29E-08, 9.17E-09};
1406 :
1407 : static const double no2[50] = { 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05,
1408 : 2.30E-05, 2.30E-05, 2.30E-05, 2.30E-05, 2.32E-05,
1409 : 2.38E-05, 2.62E-05, 3.15E-05, 4.45E-05, 7.48E-05,
1410 : 1.71E-04, 3.19E-04, 5.19E-04, 7.71E-04, 1.06E-03,
1411 : 1.39E-03, 1.76E-03, 2.16E-03, 2.58E-03, 3.06E-03,
1412 : 3.74E-03, 4.81E-03, 6.16E-03, 7.21E-03, 7.28E-03,
1413 : 6.26E-03, 4.03E-03, 2.17E-03, 1.15E-03, 6.66E-04,
1414 : 4.43E-04, 3.39E-04, 2.85E-04, 2.53E-04, 2.31E-04,
1415 : 2.15E-04, 2.02E-04, 1.92E-04, 1.83E-04, 1.76E-04,
1416 : 1.70E-04, 1.64E-04, 1.59E-04, 1.55E-04, 1.51E-04};
1417 :
1418 : static const double hcl[50] = { 1.00E-03, 7.49E-04, 5.61E-04, 4.22E-04, 3.19E-04,
1419 : 2.39E-04, 1.79E-04, 1.32E-04, 9.96E-05, 7.48E-05,
1420 : 5.68E-05, 4.59E-05, 4.36E-05, 6.51E-05, 1.01E-04,
1421 : 1.63E-04, 2.37E-04, 3.13E-04, 3.85E-04, 4.42E-04,
1422 : 4.89E-04, 5.22E-04, 5.49E-04, 5.75E-04, 6.04E-04,
1423 : 6.51E-04, 7.51E-04, 9.88E-04, 1.28E-03, 1.57E-03,
1424 : 1.69E-03, 1.74E-03, 1.76E-03, 1.79E-03, 1.80E-03,
1425 : 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03,
1426 : 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03,
1427 : 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03, 1.80E-03};
1428 :
1429 : static const double hcn[50] = { 1.70E-04, 1.65E-04, 1.63E-04, 1.61E-04, 1.60E-04,
1430 : 1.60E-04, 1.60E-04, 1.60E-04, 1.60E-04, 1.60E-04,
1431 : 1.60E-04, 1.60E-04, 1.60E-04, 1.59E-04, 1.57E-04,
1432 : 1.55E-04, 1.52E-04, 1.49E-04, 1.45E-04, 1.41E-04,
1433 : 1.37E-04, 1.34E-04, 1.30E-04, 1.25E-04, 1.19E-04,
1434 : 1.13E-04, 1.05E-04, 9.73E-05, 9.04E-05, 8.46E-05,
1435 : 8.02E-05, 7.63E-05, 7.30E-05, 7.00E-05, 6.70E-05,
1436 : 6.43E-05, 6.21E-05, 6.02E-05, 5.88E-05, 5.75E-05,
1437 : 5.62E-05, 5.50E-05, 5.37E-05, 5.25E-05, 5.12E-05};
1438 :
1439 :
1440 3448 : bool v2018 = false;
1441 :
1442 : // --------------------------------------------[ altitude for interpolation ]
1443 :
1444 3448 : if((ha < 0.0) || (ha > 120.0)) {
1445 :
1446 0 : o3den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1447 0 : n2oden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1448 0 : coden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1449 0 : no2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1450 0 : so2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1451 0 : hclden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1452 0 : hcnden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1453 :
1454 : } else {
1455 :
1456 3448 : i1 = 0; i2 = 0; i3 = 0;
1457 3448 : x1 = 0.0; x2 = 0.0; x3 = 0.0;
1458 :
1459 48761 : for(i_layer = 0; i_layer < 50; i_layer++) {
1460 :
1461 48760 : if(ha < alt[i_layer]) {
1462 :
1463 3447 : if(i_layer == 49) {
1464 :
1465 1 : i1 = i_layer - 2; i2 = i_layer - 1; i3 = i_layer;
1466 :
1467 : } else {
1468 :
1469 3446 : if(i_layer == 0) {
1470 :
1471 0 : i1 = i_layer; i2 = i_layer + 1; i3 = i_layer + 2;
1472 :
1473 : } else {
1474 :
1475 3446 : i1 = i_layer - 1; i2 = i_layer; i3 = i_layer + 1;
1476 :
1477 : }
1478 :
1479 : }
1480 :
1481 3447 : x1 = alt[i1]; x2 = alt[i2]; x3 = alt[i3];
1482 3447 : goto calc;
1483 : }
1484 : }
1485 :
1486 1 : calc:
1487 :
1488 3448 : if(x1 == 0.0 && x2 == 0.0 && x3 == 0.0) {
1489 :
1490 1 : o3den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1491 1 : n2oden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1492 1 : coden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1493 1 : no2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1494 1 : so2den = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1495 1 : hclden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1496 1 : hcnden = NumberDensity(0.0, NumberDensity::UnitInverseCubicMeter);
1497 :
1498 : } else {
1499 :
1500 : // ----- ---------------------------------------[ pressure ]
1501 :
1502 : // p=poli2(ha,x1,x2,x3,p[tip][i1],p[tip][i2],p[tip][i3]);
1503 :
1504 : // --------------------------------------------[ TEMPERATURE ]
1505 :
1506 : // t=poli2(ha,x1,x2,x3,t[tip][i1],t[tip][i2],t[tip][i3]);
1507 :
1508 : // --------------------------------------------[ DENSITY ]
1509 :
1510 3447 : d = poli2(ha,x1,x2,x3,den[tip - 1][i1],den[tip - 1][i2],den[tip - 1][i3])
1511 3447 : * airmwt * 1e6 / avogad; // en g/m**3
1512 :
1513 : // D=A+B*HA+C*HA2 en cm**-3, pero dividido por 1E19 D=D*airmwt*1e6/avogad // en g/m**3
1514 :
1515 : // --------------------------------------------[ H2O ]
1516 :
1517 : // rr=poli2(ha,x1,x2,x3,agua[tip][i1],agua[tip][i2],agua[tip][i3]); // RR=A+B*HA+C*HA2 en ppmv
1518 : // rr=rr*1e-6*(h2omwt/airmwt)*d; // en g/m**3
1519 :
1520 : // --------------------------------------------[ OZONE ]
1521 :
1522 6894 : o3den= NumberDensity(poli2(ha, x1, x2, x3, ozone[tip - 1][i1], ozone[tip- 1][i2], ozone[tip - 1][i3])
1523 3447 : * 1e-12 * d * avogad / airmwt,NumberDensity::UnitInverseCubicCentiMeter);
1524 :
1525 3447 : if(v2018){
1526 :
1527 : }else{
1528 3447 : if (seasonal_diurnal == 1) {
1529 : // std::cout << "MODEL Winter&Day (W/D): 20210624 20220511 20220825 20220828 20220901 20220903 20220905" << std::endl;
1530 3398 : if(i2 < 29){
1531 3141 : if(i2 < 24){o3den=o3den*0.650;}else{o3den=o3den*0.860;} //0.770 0.830 1.35 1.65 2.40
1532 : }else{
1533 257 : if(i2 < 32){o3den=o3den*1.35;}else{if(i2 < 35){o3den=o3den*1.60;}else{o3den=o3den*2.45;} } // 1.45 1.50 2.40 on APRIL 24, 2022
1534 : }
1535 49 : }else if(seasonal_diurnal == 2){
1536 : // std::cout << "MODEL Winter&Night (W/N): 20220827" << std::endl;
1537 0 : if(i2 < 29){
1538 0 : if(i2 < 24){o3den=o3den*0.650;}else{o3den=o3den*0.890;}
1539 : }else{
1540 0 : if(i2 < 32){o3den=o3den*1.45;}else{if(i2 < 35){o3den=o3den*1.6;}else{o3den=o3den*2.55;}}
1541 : }
1542 49 : }else if(seasonal_diurnal == 3){
1543 : // std::cout << "MODEL Summer&Day (S/D): 20201206 end of observing run" << std::endl;
1544 0 : if(i2 < 30){
1545 0 : if(i2 < 24){o3den=o3den*0.500;}else{o3den=o3den*0.800;} //PROBAR 0.600 0.850
1546 : }else{
1547 0 : if(i2 < 32){o3den=o3den*1.65;}else{if(i2 < 35){o3den=o3den*1.70;}else{o3den=o3den*2.40;}} //PROBAR 1.67 2.25 2.95
1548 : }
1549 49 : }else if(seasonal_diurnal == 4){
1550 : // std::cout << "MODEL Summer&Night (S/N): 20201021 20201205 20201206 20211031" << std::endl;
1551 0 : if(i2 < 30){
1552 0 : if(i2 < 24){o3den=o3den*0.500;}else{o3den=o3den*0.800;} //PROBAR 0.600 0.850
1553 : }else{
1554 0 : if(i2 < 32){o3den=o3den*1.65;}else{if(i2 < 35){o3den=o3den*1.70;}else{o3den=o3den*2.50;}} //PROBAR 1.67 2.25 2.95
1555 : }
1556 : }else{
1557 : // MODEL St76 Default
1558 : }
1559 : }
1560 :
1561 : // std::cout << "ha o3den " << ha << " " << o3den.get() << std::endl;
1562 :
1563 : // OZONO=A+B*HA+C*HA2 // en ppmv
1564 :
1565 : // --------------------------------------------[ N2O ]
1566 :
1567 6894 : n2oden = NumberDensity(poli2(ha, x1, x2, x3, n2o[tip - 1][i1], n2o[tip - 1][i2], n2o[tip - 1][i3])
1568 3447 : * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1569 3447 : if(v2018){
1570 :
1571 : }else{
1572 3447 : if(ha>36.0){n2oden=n2oden*1.5; }else if(ha>30){ n2oden=n2oden*2.5; }else{ }
1573 : }
1574 :
1575 :
1576 : // 12FEB2024 increasing N2O above 30m to match observed lines at Chajnantor
1577 : // N2O=A+B*HA+C*HA2 // en ppmv
1578 :
1579 : // --------------------------------------------[ NO2 ]
1580 :
1581 6894 : no2den = NumberDensity(poli2(ha, x1, x2, x3, no2[i1], no2[i2], no2[i3])
1582 3447 : * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1583 :
1584 : // no2den = no2den * 100.0; // USADO EN PRUEBAS 20240117
1585 : // N2O=A+B*HA+C*HA2 // en ppmv
1586 :
1587 : // --------------------------------------------[ HCl ]
1588 :
1589 3447 : hclden = NumberDensity(poli2(ha, x1, x2, x3, hcl[i1], hcl[i2], hcl[i3]) * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1590 :
1591 3447 : if(v2018){
1592 0 : hclden=0.0;
1593 : }else{
1594 :
1595 3447 : if(i2 < 30){
1596 3170 : hclden=hclden*1.0; // FIX 15/1/2024
1597 : }else{
1598 277 : hclden=hclden*3.5;
1599 : }
1600 :
1601 : }
1602 :
1603 : // --------------------------------------------[ HCN ]
1604 :
1605 3447 : hcnden = NumberDensity(poli2(ha, x1, x2, x3, hcn[i1], hcn[i2], hcn[i3]) * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1606 :
1607 3447 : if(v2018){
1608 0 : hcnden=0.0;
1609 : }else{
1610 :
1611 3447 : if(i2 < 25){
1612 2915 : hcnden=hcnden*0.9; // FIX 15/1/2024
1613 : }else{
1614 532 : hcnden=hcnden*2.5;
1615 : }
1616 :
1617 : }
1618 :
1619 :
1620 : // --------------------------------------------[ SO2 ]
1621 :
1622 3447 : so2den = NumberDensity(poli2(ha, x1, x2, x3, so2[i1], so2[i2], so2[i3]) * 1e-12 * d * avogad / airmwt, NumberDensity::UnitInverseCubicCentiMeter);
1623 :
1624 : // --------------------------------------------[ CO ]
1625 :
1626 6894 : coden = NumberDensity(poli2(ha,x1,x2,x3,co[tip-1][i1],co[tip-1][i2],co[tip-1][i3])
1627 3447 : * 1e-12 * d * avogad / airmwt,NumberDensity::UnitInverseCubicCentiMeter);
1628 :
1629 :
1630 :
1631 : }
1632 :
1633 : }
1634 :
1635 3448 : minorden.push_back(o3den);
1636 3448 : minorden.push_back(n2oden);
1637 3448 : minorden.push_back(coden);
1638 3448 : minorden.push_back(no2den);
1639 3448 : minorden.push_back(so2den);
1640 3448 : minorden.push_back(hclden);
1641 3448 : minorden.push_back(hcnden);
1642 :
1643 6896 : return minorden;
1644 3448 : }
1645 :
1646 27576 : double AtmProfile::poli2(double ha,
1647 : double x1,
1648 : double x2,
1649 : double x3,
1650 : double y1,
1651 : double y2,
1652 : double y3) const
1653 : {
1654 : double a, b, c;
1655 :
1656 27576 : c = (y3 - y2) * (x2 - x1) - (y2 - y1) * (x3 - x2);
1657 27576 : b = (x2 - x1) * (x3 * x3 - x2 * x2) - (x2 * x2 - x1 * x1) * (x3 - x2);
1658 27576 : c = c / b;
1659 27576 : b = (y2 - y1) - c * (x2 * x2 - x1 * x1);
1660 27576 : b = b / (x2 - x1);
1661 27576 : a = y1 - c * x1 * x1 - b * x1;
1662 :
1663 27576 : return a + b * ha + c * pow(ha, 2);
1664 : }
1665 :
1666 125 : unsigned int AtmProfile::mkAtmProfile()
1667 : {
1668 125 : const int MAXNHX = 20;
1669 : static const double
1670 : hx[MAXNHX] = { 9.225, 10.225, 11.225, 12.850, 14.850, 16.850, 18.850, 22.600, 26.600, 30.600,
1671 : 34.850, 40.850, 46.850, 52.850, 58.850, 65.100, 73.100, 81.100, 89.100, 95.600 };
1672 :
1673 : static const double
1674 : px[6][MAXNHX] = { { 0.3190E+03, 0.2768E+03, 0.2391E+03, 0.1864E+03, 0.1354E+03, 0.9613E+02, 0.6833E+02,
1675 : 0.3726E+02, 0.2023E+02, 0.1121E+02, 0.6142E+01, 0.2732E+01, 0.1260E+01, 0.6042E+00,
1676 : 0.2798E+00, 0.1202E+00, 0.3600E-01, 0.9162E-02, 0.2076E-02, 0.6374E-03 },
1677 : { 0.3139E+03, 0.2721E+03, 0.2350E+03, 0.1833E+03, 0.1332E+03, 0.9726E+02, 0.7115E+02,
1678 : 0.3992E+02, 0.2185E+02, 0.1216E+02, 0.6680E+01, 0.2985E+01, 0.1400E+01, 0.6780E+00,
1679 : 0.3178E+00, 0.1380E+00, 0.4163E-01, 0.9881E-02, 0.2010E-02, 0.5804E-03 },
1680 : { 0.2892E+03, 0.2480E+03, 0.2124E+03, 0.1649E+03, 0.1206E+03, 0.8816E+02, 0.6433E+02,
1681 : 0.3558E+02, 0.1901E+02, 0.1014E+02, 0.5316E+01, 0.2255E+01, 0.1022E+01, 0.4814E+00,
1682 : 0.2206E+00, 0.9455E-01, 0.3000E-01, 0.8729E-02, 0.2332E-02, 0.8164E-03 },
1683 : { 0.3006E+03, 0.2587E+03, 0.2223E+03, 0.1739E+03, 0.1288E+03, 0.9495E+02, 0.7018E+02,
1684 : 0.3983E+02, 0.2199E+02, 0.1232E+02, 0.6771E+01, 0.3056E+01, 0.1452E+01, 0.7051E+00,
1685 : 0.3353E+00, 0.1459E+00, 0.4431E-01, 0.1024E-01, 0.1985E-02, 0.5627E-03 },
1686 : { 0.2731E+03, 0.2335E+03, 0.1995E+03, 0.1546E+03, 0.1130E+03, 0.8252E+02, 0.6017E+02,
1687 : 0.3314E+02, 0.1751E+02, 0.9306E+01, 0.4826E+01, 0.1988E+01, 0.8645E+00, 0.4000E+00,
1688 : 0.1819E+00, 0.7866E-01, 0.2639E-01, 0.8264E-02, 0.2364E-02, 0.8439E-03 },
1689 : { 0.2978E+03, 0.2559E+03, 0.2191E+03, 0.1698E+03, 0.1240E+03, 0.9060E+02, 0.6621E+02,
1690 : 0.3688E+02, 0.1999E+02, 0.1087E+02, 0.5862E+01, 0.2565E+01, 0.1182E+01, 0.5572E+00,
1691 : 0.2551E+00, 0.1074E+00, 0.3224E-01, 0.8697E-02, 0.2158E-02, 0.6851E-03 } };
1692 :
1693 : static const double
1694 : tx[6][MAXNHX] = { { 0.2421E+03, 0.2354E+03, 0.2286E+03, 0.2180E+03, 0.2046E+03, 0.1951E+03, 0.2021E+03,
1695 : 0.2160E+03, 0.2250E+03, 0.2336E+03, 0.2428E+03, 0.2558E+03, 0.2686E+03, 0.2667E+03,
1696 : 0.2560E+03, 0.2357E+03, 0.2083E+03, 0.1826E+03, 0.1767E+03, 0.1841E+03 },
1697 : { 0.2403E+03, 0.2338E+03, 0.2273E+03, 0.2167E+03, 0.2157E+03, 0.2157E+03, 0.2177E+03,
1698 : 0.2223E+03, 0.2271E+03, 0.2349E+03, 0.2448E+03, 0.2596E+03, 0.2740E+03, 0.2727E+03,
1699 : 0.2603E+03, 0.2394E+03, 0.2045E+03, 0.1715E+03, 0.1644E+03, 0.1785E+03 },
1700 : { 0.2239E+03, 0.2196E+03, 0.2191E+03, 0.2183E+03, 0.2173E+03, 0.2163E+03, 0.2153E+03,
1701 : 0.2152E+03, 0.2152E+03, 0.2178E+03, 0.2275E+03, 0.2458E+03, 0.2637E+03, 0.2633E+03,
1702 : 0.2531E+03, 0.2407E+03, 0.2243E+03, 0.2072E+03, 0.1991E+03, 0.2091E+03 },
1703 : { 0.2301E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03, 0.2252E+03,
1704 : 0.2252E+03, 0.2298E+03, 0.2361E+03, 0.2468E+03, 0.2648E+03, 0.2756E+03, 0.2763E+03,
1705 : 0.2662E+03, 0.2392E+03, 0.2023E+03, 0.1681E+03, 0.1609E+03, 0.1770E+03 },
1706 : { 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2172E+03, 0.2161E+03, 0.2149E+03,
1707 : 0.2126E+03, 0.2127E+03, 0.2165E+03, 0.2222E+03, 0.2368E+03, 0.2516E+03, 0.2601E+03,
1708 : 0.2523E+03, 0.2486E+03, 0.2388E+03, 0.2215E+03, 0.2033E+03, 0.2113E+03 },
1709 : { 228.260, 221.838, 216.778, 216.700, 216.700, 216.700, 216.700, 219.200, 223.136,
1710 : 227.340, 236.110, 252.746, 268.936, 265.057, 250.174, 233.026, 212.656, 196.466,
1711 : 187.260, 189.204} };
1712 :
1713 125 : double T_ground = groundTemperature_.get(Temperature::UnitKelvin); // std::cout<<"T_ground: " << T_ground <<"K"<<endl;
1714 125 : double P_ground = groundPressure_.get(Pressure::UnitMilliBar); // std::cout<<"P_ground: " << P_ground <<"mb"<<endl;
1715 125 : double rh = relativeHumidity_.get(Percent::UnitPercent); // std::cout<<"rh: " << rh <<"%"<<endl;
1716 125 : double h0 = wvScaleHeight_.get(Length::UnitKiloMeter); // std::cout<<"h0: " << h0 <<"km"<<endl;
1717 125 : double dp = pressureStep_.get(Pressure::UnitMilliBar); // std::cout<<"dp: " << dp <<"mb"<<endl;
1718 125 : double alti = altitude_.get(Length::UnitKiloMeter); // std::cout<<"alti: " << alti <<"km"<<endl;
1719 125 : double atmh = topAtmProfile_.get(Length::UnitKiloMeter); // std::cout<<"atmh: " << atmh <<"km"<<endl;
1720 125 : double dp1 = pressureStepFactor_; // std::cout<<"dp1: " << dp1 <<" "<<endl;
1721 125 : double dt = tropoLapseRate_; // TODO implementer des unites (K/km) ici localement
1722 : double prLimit;
1723 : double minmin;
1724 125 : unsigned int seasonal_diurnal = 0;
1725 :
1726 125 : if(dp<5.0 || dp>12.5){dp = 10.0; std::cout << "ATMProfile WARNING: Primary Pressure step fixed to 10 mb for numerical reasons" << std::endl; }
1727 125 : if(dp1<1.074 || dp1 > 1.076 ){dp = 1.075; std::cout << "ATMProfile WARNING: Pressure Step factor fixed to 1.075 for numerical reasons" << std::endl; }
1728 :
1729 125 : if(typeAtm_ == 1) { // TROPICAL US St. 76
1730 125 : prLimit = 125; //100.0; // 100.0 230.2; TEMPORAL FIX 20220606
1731 125 : seasonal_diurnal=1;
1732 0 : } else if(typeAtm_ == 11) { // CHAJNANTOR Winter + Daytime
1733 0 : typeAtm_ = 1;
1734 0 : seasonal_diurnal=1;
1735 0 : prLimit = 125.0; //100.0; // 100.0 230.2; TEMPORAL FIX 20220606
1736 0 : } else if(typeAtm_ == 12) { // CHAJNANTOR Winter + Nighttime
1737 0 : typeAtm_ = 1;
1738 0 : seasonal_diurnal=2;
1739 0 : prLimit = 125.0; //100.0; // 100.0 230.2; TEMPORAL FIX 20220606
1740 0 : } else if(typeAtm_ == 13) { // CHAJNANTOR Summer + Daytime
1741 0 : typeAtm_ = 1;
1742 0 : seasonal_diurnal=3;
1743 0 : prLimit = 125.0; //100.0; // 100.0 230.2; TEMPORAL FIX 20220606
1744 0 : } else if(typeAtm_ == 14) { // CHAJNANTOR Summer + Nighttime
1745 0 : typeAtm_ = 1;
1746 0 : seasonal_diurnal=4;
1747 0 : prLimit = 125.0;
1748 0 : } else if(typeAtm_ == 2) {
1749 0 : prLimit = 100.0; // 198.0;
1750 0 : seasonal_diurnal=1;
1751 0 : } else if(typeAtm_ == 3) {
1752 0 : prLimit = 100.0; // 300.0
1753 0 : seasonal_diurnal=1;
1754 0 : } else if(typeAtm_ == 4) {
1755 0 : prLimit = 100.0; // 311.0
1756 0 : seasonal_diurnal=1;
1757 0 : } else if(typeAtm_ == 5) {
1758 0 : prLimit = 100.0; // 332.0
1759 0 : seasonal_diurnal=1;
1760 0 : } else if(typeAtm_ == 6) {
1761 0 : prLimit = 100.0; // 332.0
1762 0 : seasonal_diurnal=1;
1763 : } else {
1764 0 : typeAtm_ = 6; // Default: US St76 Atm.
1765 0 : prLimit = 100.0; // 250.0;
1766 0 : seasonal_diurnal=1;
1767 : }
1768 :
1769 125 : unsigned int npp = 0; // number of layers initialized
1770 :
1771 125 : double rt = 6371.2E+0; // Earth radius in km
1772 125 : double g0 = 9.80665E+0; // Earth gravity at the surface (m/s**2)
1773 :
1774 : // static bool first = true; // [-Wunused_but_set_variable]
1775 :
1776 : unsigned int i;
1777 : unsigned int i0;
1778 : unsigned int j;
1779 : unsigned int k;
1780 : // double wh2o, wh2o0; // [-Wunused_but_set_variable]
1781 : double wgr0;
1782 125 : double g = g0;
1783 : double www, altura;
1784 : // double abun_ozono, abun_n2o, abun_co;
1785 125 : NumberDensity ozono, n2o, co;
1786 : double wgr, dh;
1787 : // double humrel; // [-Wunused_but_set_variable]
1788 : // bool errcode;
1789 :
1790 : //int nmaxLayers=40; // FV peut etre devrions nos avoir un garde-fou au cas ou le nb de couches serait stupidement trop grand
1791 :
1792 125 : vector<double> v_layerPressure;
1793 125 : vector<double> v_layerPressure0;
1794 125 : vector<double> v_layerPressure1;
1795 125 : vector<double> v_layerTemperature;
1796 125 : vector<double> v_layerTemperature0;
1797 125 : vector<double> v_layerTemperature1;
1798 125 : vector<double> v_layerThickness;
1799 125 : vector<double> v_layerWaterVapor;
1800 125 : vector<double> v_layerWaterVapor0;
1801 125 : vector<double> v_layerWaterVapor1;
1802 125 : vector<double> v_layerO3;
1803 125 : vector<double> v_layerCO;
1804 125 : vector<double> v_layerN2O;
1805 125 : vector<double> v_layerNO2;
1806 125 : vector<double> v_layerSO2;
1807 125 : vector<double> v_layerHCl;
1808 125 : vector<double> v_layerHCN;
1809 :
1810 125 : vector<double> v_layerO3_aux;
1811 125 : vector<double> v_layerCO_aux;
1812 125 : vector<double> v_layerN2O_aux;
1813 125 : vector<double> v_layerNO2_aux;
1814 125 : vector<double> v_layerSO2_aux;
1815 125 : vector<double> v_layerHCl_aux;
1816 125 : vector<double> v_layerHCN_aux;
1817 125 : vector<double> v_layerPressure_aux;
1818 125 : vector<double> v_layerPressure0_aux;
1819 125 : vector<double> v_layerPressure1_aux;
1820 125 : vector<double> v_layerTemperature_aux;
1821 125 : vector<double> v_layerTemperature0_aux;
1822 125 : vector<double> v_layerTemperature1_aux;
1823 125 : vector<double> v_layerThickness_aux;
1824 125 : vector<double> v_layerWaterVapor_aux;
1825 125 : vector<double> v_layerWaterVapor0_aux;
1826 125 : vector<double> v_layerWaterVapor1_aux;
1827 125 : vector<NumberDensity> minorden;
1828 :
1829 125 : v_layerPressure.push_back(P_ground);
1830 125 : v_layerThickness.push_back(alti * 1000);
1831 125 : v_layerTemperature.push_back(T_ground);
1832 :
1833 250 : wgr = rwat(Temperature(v_layerTemperature[0], Temperature::UnitKelvin),
1834 250 : Humidity(rh, Percent::UnitPercent),
1835 375 : Pressure(v_layerPressure[0], Pressure::UnitMilliBar)).get(MassDensity::UnitGramPerCubicMeter);
1836 :
1837 : // Absolute Humidity in gr/m**3 at altitude alti
1838 125 : wgr0 = wgr * exp(alti / h0); // h0 in km ==> wgr0 would be the absolute humidity (gr/m**3) at sea level
1839 : // wh2o0 = wgr0 * h0; // wh2o0 is the zenith column of water above sea level (kgr/m**2 or mm) // [-Wunused_but_set_variable]
1840 : // wh2o = wh2o0 * exp(-alti / h0); // wh2o is the zenith column of water above alti (kgr/m**2 or mm) // [-Wunused_but_set_variable]
1841 :
1842 125 : v_layerWaterVapor.push_back(wgr); // in gr/m**3
1843 :
1844 125 : i0 = 0;
1845 125 : i = 0;
1846 125 : j = 0;
1847 125 : npp = 0;
1848 :
1849 : // std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1850 :
1851 125 : bool control = true;
1852 : // bool control2 = true;
1853 : while(true) {
1854 3523 : i++;
1855 :
1856 : // std::cout << "layer: " << i << " " << v_layerPressure[i - 1] - dp * pow(dp1, i - 1) << " " << prLimit << std::endl;
1857 :
1858 3523 : if( (v_layerPressure[i - 1] - dp * pow(dp1, i - 1) <= prLimit ) ) { //&& (v_layerTemperature[i - 1] <= tx[typeAtm_ - 1][0]) ) { // &&(fabs(prLimit-v_layerPressure[i - 1]-dp*pow(dp1,i-1))>=20.0) ) {
1859 :
1860 1085 : if(control) {
1861 :
1862 : //std::cout << "capa " << i-1 << " temperature=" << v_layerTemperature[i - 1] << " K " << " tx[0]=" << tx[typeAtm_ - 1][0] << std::endl;
1863 :
1864 : // std::cout << "prLimit=" << prLimit << std::endl;
1865 : // std::cout << "aaa " << v_layerPressure[i - 1] - dp * pow(dp1, i - 1) << std::endl;
1866 :
1867 : /* for(k = 0; k < 20; k++) {
1868 : if(control2) {
1869 : if(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) >= px[typeAtm_ - 1][k]) {
1870 : j = k;
1871 : std::cout << "P=" << v_layerPressure[i - 1] - dp * pow(dp1, i - 1)
1872 : << " prLimit=" << prLimit << " px[" << typeAtm_ - 1 << "][" << k << "]="
1873 : << px[typeAtm_ - 1][k] << std::endl;
1874 : control2 = false;
1875 : }
1876 : }
1877 : } */
1878 :
1879 125 : minmin = 20000.0;
1880 :
1881 2625 : for(k = 0; k < MAXNHX; k++) {
1882 :
1883 4328 : if( (fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) > 1.05*px[typeAtm_ - 1][k])) &&
1884 1828 : (fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) - px[typeAtm_ - 1][k])) <= minmin ) {
1885 :
1886 125 : j = k;
1887 :
1888 : /* std::cout << "P=" << v_layerPressure[i - 1] - dp * pow(dp1, i - 1)
1889 : << " prLimit=" << prLimit << " px[" << typeAtm_ - 1 << "][" << k << "]="
1890 : << px[typeAtm_ - 1][k] << std::endl; FIXED 20240606 */
1891 125 : minmin = fabs(v_layerPressure[i - 1] - dp * pow(dp1, i - 1) - px[typeAtm_ - 1][k]);
1892 : // std::cout << " minmin=" << minmin << std::endl; FIXED 20240606
1893 :
1894 : }
1895 :
1896 : }
1897 :
1898 :
1899 : }
1900 :
1901 : // std::cout << "i,j,v_layerPressure.size()-1=" << i << "," << j << "," << v_layerPressure.size() - 1 << std::endl;
1902 :
1903 1085 : if(i0 == 0) i0 = i - 1;
1904 :
1905 : // std::cout << "i0,j " << i0 << "," << j << std::endl;
1906 :
1907 :
1908 1085 : if(i < v_layerPressure.size() - 1) {
1909 0 : if(control) {
1910 0 : int j0 = (j>0) ? j-1:0;
1911 0 : v_layerPressure[i] = px[typeAtm_ - 1][j0];
1912 0 : v_layerTemperature[i] = tx[typeAtm_ - 1][j0];
1913 :
1914 0 : std::cout << "v_layerTemperature[i]=" << v_layerTemperature[i] << std::endl;
1915 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]; COMMENTED OUT 31/5/2017
1916 :
1917 0 : www = v_layerWaterVapor[i - 1] / 1000.0;
1918 0 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1919 0 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1920 0 : v_layerThickness[i] = v_layerThickness[i - 1] + dh;
1921 : } else {
1922 0 : j++;
1923 0 : v_layerPressure[i] = px[typeAtm_ - 1][j - 1];
1924 0 : v_layerTemperature[i] = tx[typeAtm_ - 1][j - 1];
1925 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]; COMMENTED OUT 31/5/2017
1926 0 : if(j < MAXNHX) {
1927 0 : v_layerThickness[i] = (hx[j] - hx[j - 1]) * 1000.0 + v_layerThickness[i - 1];
1928 : }
1929 : }
1930 : // std::cout << "layer " << i << " j=" << j << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1931 0 : v_layerWaterVapor[i] = wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0));
1932 : } else {
1933 1085 : if(control) {
1934 125 : int j0 = (j>0) ? j-1:0;
1935 125 : v_layerPressure.push_back(px[typeAtm_ - 1][j0]);
1936 125 : v_layerTemperature.push_back(tx[typeAtm_ - 1][j0]);
1937 : //- tx[typeAtm_ - 1][0] + v_layerTemperature[i0]); COMMENTED OUT 31/5/2017
1938 :
1939 125 : www = v_layerWaterVapor[i - 1] / 1000.0;
1940 125 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1941 125 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1942 125 : v_layerThickness.push_back(v_layerThickness[i - 1] + dh);
1943 : } else {
1944 960 : j++;
1945 960 : v_layerPressure.push_back(px[typeAtm_ - 1][j - 1]);
1946 960 : v_layerTemperature.push_back(tx[typeAtm_ - 1][j - 1]);
1947 : // - tx[typeAtm_ - 1][0] + v_layerTemperature[i0]); COMMENTED OUT 31/5/2017
1948 960 : if(j < MAXNHX) {
1949 959 : v_layerThickness.push_back((hx[j] - hx[j - 1]) * 1000.0 + v_layerThickness[i - 1]);
1950 : }
1951 :
1952 : }
1953 1085 : v_layerWaterVapor.push_back(wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0)));
1954 : }
1955 : // std::cout << "layer " << i << " j=" << j << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1956 1085 : if(control) {
1957 125 : tropoLayer_ = i - 1;
1958 : // std::cout << "tropoLayer_=" << tropoLayer_ << std::endl;
1959 125 : control = false;
1960 : }
1961 : // std::cout << "2 layer " << i << " alt at end of layer = " << v_layerThickness[i] << " atmh=" << atmh*1000 << std::endl;
1962 :
1963 : // besides the thickness threshold, if j goes >=MAXNHX this loop cannot iterate further
1964 1085 : if((j >= MAXNHX) || (v_layerThickness[i] > (atmh * 1000.0)))
1965 125 : break;
1966 : } else {
1967 :
1968 : // std::cout << "i,j,v_layerPressure.size()-1=" << i << "," << j << "," << v_layerPressure.size() - 1 << std::endl;
1969 :
1970 2438 : if(i > v_layerPressure.size() - 1) {
1971 2438 : v_layerPressure.push_back(v_layerPressure[i - 1] - dp * pow(dp1, i - 1));
1972 :
1973 2438 : www = v_layerWaterVapor[i - 1] / 1000.0; // in kg/m**3
1974 2438 : g = g0 * pow(1. + ((v_layerThickness[i - 1] / 1000.0)) / rt, -2); // gravity corrected for the height
1975 2438 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1976 2438 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1977 2438 : v_layerThickness.push_back(v_layerThickness[i - 1] + dh);
1978 : //std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1979 2438 : v_layerTemperature.push_back(v_layerTemperature[i - 1] + dt * dh / 1000.0);
1980 2438 : v_layerWaterVapor.push_back(wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0))); //r[i] in kgr/(m**2*1000m) [gr/m**3]
1981 : } else {
1982 0 : v_layerPressure[i] = v_layerPressure[i - 1] - dp * pow(dp1, i - 1);
1983 0 : www = v_layerWaterVapor[i - 1] / 1000.0; // in kg/m**3
1984 0 : g = g0 * pow(1. + ((v_layerThickness[i - 1] / 1000.0)) / rt, -2); // gravity corrected for the height
1985 0 : dh = 288.6948 * v_layerTemperature[i - 1] * (1.0 + 0.61 * www / 1000.0)
1986 0 : * log(v_layerPressure[i - 1] / v_layerPressure[i]) / g;
1987 0 : v_layerThickness[i] = v_layerThickness[i - 1] + dh;
1988 : //std::cout << "layer " << i << " v_layerThickness[" << i << "]=" << v_layerThickness[i] << " v_layerPressure[" << i << "]=" << v_layerPressure[i] << std::endl;
1989 0 : v_layerTemperature[i] = v_layerTemperature[i - 1] + dt * dh / 1000.0;
1990 0 : v_layerWaterVapor[i] = wgr0 * exp(-v_layerThickness[i] / (1000.0 * h0)); //r[i] in kgr/(m**2*1000m) [gr/m**3]
1991 : }
1992 :
1993 : //std::cout << "dh=" << dh << std::endl;
1994 2438 : if( v_layerTemperature[i] <= tx[typeAtm_ - 1][4] && v_layerPressure[i] <= px[typeAtm_ - 1][4] ) // [0]
1995 : {
1996 : // std::cout << "v_layerThickness[i]=" << v_layerThickness[i] << " " << v_layerThickness[i-1] <<std::endl;
1997 : // std::cout << "dh=" << dh << std::endl;
1998 : // std::cout << "old prLimit =" << prLimit << std::endl;
1999 11 : prLimit=v_layerPressure[i]; //14NOV2017
2000 : // std::cout << "new prLimit =" << prLimit << std::endl;
2001 : }
2002 :
2003 : // humrel = rwat_inv(Temperature(v_layerTemperature[i], Temperature::UnitKelvin),
2004 : // MassDensity(v_layerWaterVapor[i], MassDensity::UnitGramPerCubicMeter),
2005 : // Pressure(v_layerPressure[i], Pressure::UnitMilliBar)).get(Percent::UnitPercent);
2006 :
2007 :
2008 : /* cout << "layer " << i
2009 : << " P " << v_layerPressure[i] << " prLimit " << prLimit
2010 : << " T " << v_layerTemperature[i]
2011 : << " Alti (layer top in m) " << v_layerThickness[i]
2012 : << " WaterVapor " << v_layerWaterVapor[i] << std::endl; */
2013 : // std::cout << "1 layer " << i << " alt at end of layer = " << v_layerThickness[i] << " atmh=" << atmh*1000 << std::endl;
2014 2438 : if(v_layerThickness[i] > (atmh * 1000.0)) break;
2015 : }
2016 : // std::cout << "v_layerPressure[" << i << "]" << v_layerPressure[i] << " v_layerTemperature[" << i << "] " << v_layerTemperature[i] << std::endl;
2017 3398 : }
2018 :
2019 125 : npp = i - 1;
2020 :
2021 : // std::cout << "npp=" << npp << std::endl;
2022 : // std::cout << "tropoLayer=" << tropoLayer_ << std::endl;
2023 125 : tropoTemperature_ = Temperature(v_layerTemperature[tropoLayer_], Temperature::UnitKelvin);
2024 125 : tropoAltitude_ = Length(v_layerThickness[tropoLayer_], Length::UnitMeter);
2025 : // std::cout << "tropoAltitude=" << tropoAltitude_.get(Length::UnitKiloMeter) << " km" << std::endl;
2026 : // std::cout << "tropoTemperature=" << tropoTemperature_.get(Temperature::UnitKelvin) << " K" << std::endl;
2027 : // std::cout << "ground Altitude=" << altitude_.get(Length::UnitKiloMeter) << " km" << std::endl;
2028 : // std::cout << "ground Temperature=" << groundTemperature_.get(Temperature::UnitKelvin) << " K" << std::endl;
2029 : // std::cout << "Calculated Lapse Rate=" << (tropoTemperature_.get(Temperature::UnitKelvin)-groundTemperature_.get(Temperature::UnitKelvin))/(tropoAltitude_.get(Length::UnitKiloMeter)-altitude_.get(Length::UnitKiloMeter)) << " K/km" << std::endl;
2030 :
2031 :
2032 125 : altura = alti;
2033 :
2034 : /*
2035 : if(first){
2036 : v_layerO3_.reserve(npp);
2037 : v_layerCO_.reserve(npp);
2038 : v_layerN2O_.reserve(npp);
2039 : v_layerNO2_.reserve(npp);
2040 : v_layerSO2_.reserve(npp);
2041 : v_layerHCl_.reserve(npp);
2042 : v_layerHCN_.reserve(npp);
2043 : }
2044 : */
2045 :
2046 :
2047 3523 : for(unsigned int jj = 0; jj < npp; jj++) {
2048 :
2049 3398 : v_layerTemperature0.push_back(v_layerTemperature[jj]);
2050 3398 : v_layerTemperature1.push_back(v_layerTemperature[jj + 1]);
2051 3398 : v_layerPressure0.push_back(v_layerPressure[jj]);
2052 3398 : v_layerPressure1.push_back(v_layerPressure[jj + 1]);
2053 3398 : v_layerWaterVapor0.push_back(1.0E-3*v_layerWaterVapor[jj]);
2054 3398 : v_layerWaterVapor1.push_back(1.0E-3*v_layerWaterVapor[jj + 1]);
2055 :
2056 : }
2057 :
2058 3523 : for(j = 0; j < npp; j++) {
2059 :
2060 3398 : v_layerThickness[j] = (v_layerThickness[j + 1] - v_layerThickness[j]); // in m
2061 3398 : altura = altura + v_layerThickness[j] / 2.0E3; // in km
2062 3398 : v_layerTemperature[j] = (v_layerTemperature[j + 1] + v_layerTemperature[j]) // AVERAGING MAY/30/2017
2063 3398 : / 2.; // in K
2064 3398 : v_layerPressure[j] = exp((log(v_layerPressure[j + 1])
2065 3398 : + log(v_layerPressure[j])) / 2.0); // in mb
2066 3398 : if(altura < 40.0 ){
2067 6532 : v_layerWaterVapor[j] = 1.0E-3 * exp((log(v_layerWaterVapor[j + 1]) + // units:
2068 3266 : log(v_layerWaterVapor[j])) / 2.0); // kg/m**3
2069 : }else{
2070 132 : v_layerWaterVapor[j] = v_layerWaterVapor[j-1]; // 12FEB2024 TO KEEP MESOSPHERIC WATER CONSTANT
2071 : }
2072 : // std::cout << "type_=" << type_ << std::endl;
2073 : // std::cout << "typeAtm_=" << typeAtm_ << std::endl;
2074 :
2075 3398 : unsigned int atmType = typeAtm_; // conversion in int
2076 :
2077 : // std::cout << "going to minorden with atmType=" << atmType << std::endl;
2078 :
2079 : // std::cout << "Going to st76 seasonal_diurnal=" << seasonal_diurnal << std::endl;
2080 3398 : minorden = st76(Length(altura, Length::UnitKiloMeter), atmType, seasonal_diurnal);
2081 :
2082 : // std::cout << "Ozone: " << abun_ozono << " " << ozono.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
2083 : // std::cout << "N2O : " << abun_n2o << " " << n2o.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
2084 : // std::cout << "CO : " << abun_co << " " << co.get(NumberDensity::UnitInverseCubicCentiMeter) << std::endl;
2085 :
2086 :
2087 : /* std::cout << j << " " << v_layerO3_.size()-1 << std::endl; */
2088 :
2089 : /* if(j>v_layerO3_.size()-1){
2090 : v_layerO3_.push_back( 1.E6*abun_ozono); // in m**-3
2091 : v_layerCO_.push_back( 1.E6*abun_co ); // in m**-3
2092 : v_layerN2O_.push_back( 1.E6*abun_n2o ); // in m**-3
2093 : v_layerNO2_.push_back( 1.E6*abun_no2 ); // in m**-3
2094 : v_layerSO2_.push_back( 1.E6*abun_so2 ); // in m**-3
2095 : v_layerHCl_.push_back( 1.E6*abun_hcl ); // in m**-3
2096 : v_layerHCN_.push_back( 1.E6*abun_hcn ); // in m**-3
2097 : std::cout << "uno" << std::endl;
2098 : }else{ */
2099 :
2100 : /* v_layerO3_[j] = 1.E6*abun_ozono; // in m**-3
2101 : v_layerCO_[j] = 1.E6*abun_co; // in m**-3
2102 : v_layerN2O_[j] = 1.E6*abun_n2o; // in m**-3
2103 : v_layerNO2_[j] = 1.E6*abun_no2; // in m**-3
2104 : v_layerSO2_[j] = 1.E6*abun_so2; // in m**-3
2105 : v_layerHCl_[j] = 1.E6*abun_hcl; // in m**-3
2106 : v_layerHCN_[j] = 1.E6*abun_hcn; // in m**-3 */
2107 :
2108 : // v_layerO3.push_back(1.E6*abun_ozono); // in m**-3
2109 : // v_layerCO.push_back(1.E6*abun_co); // in m**-3
2110 : // v_layerN2O.push_back(1.E6*abun_n2o); // in m**-3
2111 : // v_layerNO2.push_back(1.E6*abun_no2); // in m**-3
2112 : // v_layerSO2.push_back(1.E6*abun_so2); // in m**-3
2113 : // v_layerHCl.push_back(1.E6*abun_hcl); // in m**-3
2114 : // v_layerHCN.push_back(1.E6*abun_hcn); // in m**-3
2115 :
2116 3398 : v_layerO3.push_back(1.E6 * minorden[0].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2117 3398 : v_layerN2O.push_back(1.E6 * minorden[1].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2118 3398 : v_layerCO.push_back(1.E6 * minorden[2].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2119 : //std::cout << "minorden[3] " << minorden[3].get("m**-3") << std::endl;
2120 3398 : v_layerNO2.push_back(1.E6 * minorden[3].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2121 3398 : v_layerSO2.push_back(1.E6 * minorden[4].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2122 3398 : v_layerHCl.push_back(1.E6 * minorden[5].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2123 3398 : v_layerHCN.push_back(1.E6 * minorden[6].get(NumberDensity::UnitInverseCubicCentiMeter)); // in m**-3
2124 :
2125 :
2126 :
2127 : /* } */
2128 :
2129 3398 : altura = altura + v_layerThickness[j] / 2.0E3;
2130 :
2131 : /* std::cout << "j=" << j << "v_layerThickness.size()=" << v_layerThickness.size() << std::endl; */
2132 :
2133 : }
2134 :
2135 : /* std::cout << "j=" << j << " v_layerThickness_aux.size()= " << v_layerThickness_aux.size() << std::endl; */
2136 :
2137 : /* if(v_layerThickness.size()>npp){
2138 : for(j=npp; j<v_layerThickness.size(); j++){
2139 : delete v_layerThickness[j];
2140 : }
2141 : } */
2142 :
2143 3523 : for(j = 0; j < npp; j++) {
2144 :
2145 3398 : v_layerPressure_aux.push_back(v_layerPressure[j]);
2146 3398 : v_layerPressure0_aux.push_back(v_layerPressure0[j]);
2147 3398 : v_layerPressure1_aux.push_back(v_layerPressure1[j]);
2148 3398 : v_layerTemperature_aux.push_back(v_layerTemperature[j]);
2149 3398 : v_layerTemperature0_aux.push_back(v_layerTemperature0[j]);
2150 3398 : v_layerTemperature1_aux.push_back(v_layerTemperature1[j]);
2151 3398 : v_layerThickness_aux.push_back(v_layerThickness[j]);
2152 3398 : v_layerWaterVapor_aux.push_back(v_layerWaterVapor[j]);
2153 3398 : v_layerWaterVapor0_aux.push_back(v_layerWaterVapor0[j]);
2154 3398 : v_layerWaterVapor1_aux.push_back(v_layerWaterVapor1[j]);
2155 3398 : v_layerO3_aux.push_back(v_layerO3[j]);
2156 3398 : v_layerCO_aux.push_back(v_layerCO[j]);
2157 3398 : v_layerN2O_aux.push_back(v_layerN2O[j]);
2158 3398 : v_layerNO2_aux.push_back(v_layerNO2[j]);
2159 3398 : v_layerSO2_aux.push_back(v_layerSO2[j]);
2160 3398 : v_layerHCl_aux.push_back(v_layerHCl[j]);
2161 3398 : v_layerHCN_aux.push_back(v_layerHCN[j]);
2162 :
2163 : }
2164 :
2165 125 : if(j > v_layerPressure_.size() - 1) { // ?????
2166 119 : v_layerPressure_.reserve(npp);
2167 119 : v_layerPressure0_.reserve(npp);
2168 119 : v_layerPressure1_.reserve(npp);
2169 119 : v_layerTemperature_.reserve(npp);
2170 119 : v_layerTemperature0_.reserve(npp);
2171 119 : v_layerTemperature1_.reserve(npp);
2172 119 : v_layerThickness_.reserve(npp);
2173 119 : v_layerWaterVapor_.reserve(npp);
2174 119 : v_layerWaterVapor0_.reserve(npp);
2175 119 : v_layerWaterVapor1_.reserve(npp);
2176 119 : v_layerO3_.reserve(npp);
2177 119 : v_layerCO_.reserve(npp);
2178 119 : v_layerN2O_.reserve(npp);
2179 119 : v_layerNO2_.reserve(npp);
2180 119 : v_layerSO2_.reserve(npp);
2181 119 : v_layerHCl_.reserve(npp);
2182 119 : v_layerHCN_.reserve(npp);
2183 : }
2184 :
2185 125 : v_layerPressure_ = v_layerPressure_aux;
2186 125 : v_layerPressure0_ = v_layerPressure0_aux;
2187 125 : v_layerPressure1_ = v_layerPressure1_aux;
2188 125 : v_layerTemperature_ = v_layerTemperature_aux;
2189 125 : v_layerTemperature0_ = v_layerTemperature0_aux;
2190 125 : v_layerTemperature1_ = v_layerTemperature1_aux;
2191 125 : v_layerThickness_ = v_layerThickness_aux;
2192 125 : v_layerWaterVapor_ = v_layerWaterVapor_aux;
2193 125 : v_layerWaterVapor0_ = v_layerWaterVapor0_aux;
2194 125 : v_layerWaterVapor1_ = v_layerWaterVapor1_aux;
2195 125 : v_layerO3_ = v_layerO3_aux;
2196 125 : v_layerCO_ = v_layerCO_aux;
2197 125 : v_layerN2O_ = v_layerN2O_aux;
2198 125 : v_layerNO2_ = v_layerNO2_aux;
2199 125 : v_layerSO2_ = v_layerSO2_aux;
2200 125 : v_layerHCl_ = v_layerHCl_aux;
2201 125 : v_layerHCN_ = v_layerHCN_aux;
2202 :
2203 : // first = false; // ????? [-Wunused_but_set_variable]
2204 :
2205 125 : return npp;
2206 125 : }
2207 :
2208 : ATM_NAMESPACE_END
2209 :
|
