radixsnd2arl.cc 25.1 KB
Newer Older
1
2
3
4
5
/*
 * Example utility to convert a vertical profile of meteorological
 * data to the ARL format
 */

6
#include <algorithm>
7
#include <cmath>
8
#include <ctime>
9
10
11
12
#include <iostream>
#include <string>
#include <vector>

13
14
#include "radixcommand/commandline.hh"
#include "radixcore/stringfunctions.hh"
15
#include "radixmath/constants.hh"
16
#include "radixmath/interpolate.hh"
17
#include "radixmath/util.hh"
18

19
#include "radixio/arldatastream.hh"
20
21
22
#include "radixio/csvfile.hh"

using namespace radix;
23

24
25
26
27
28
29
30
31
32
void addHour(struct tm *time, int hours)
{
  int seconds = hours * 60 * 60;

  time_t date_seconds = mktime(time) + seconds;

  *time = *localtime(&date_seconds);
}

33
34
int main(int argc, char **argv)
{
35
36
37
38
39
40
41
42
43
44
45
  std::cout << "************************" << std::endl;
  std::cout << "***** radixsnd2arl *****" << std::endl;
  std::cout << "************************" << std::endl;

  // Set up command line options
  CommandLine commandLine(argc, argv);
  commandLine.addOption("i", "Input csv file containing met data", true);
  commandLine.addOption("clat", "Centre latitude of output ARL file (degrees)",
                        true);
  commandLine.addOption("clon", "Centre longitude of output ARL file (degrees)",
                        true);
46
47
  commandLine.addOption("e", "Extent of output ARL file (km) [500]", false);
  commandLine.addOption("r", "Resolution of output ARL file (km) [10]", false);
48
49
  commandLine.addOption("t", "Time of data start (YYYYMMDDHH) [1951111917]",
                        false);
50
51
  commandLine.addOption("n", "Number of one hour timesteps to output [1]",
                        false);
52
53
  commandLine.addOption("g", "Add ground level elevation to height values",
                        false);
54
  commandLine.addOption("o", "Output ARL file", false);
55
56
57
  commandLine.addOption(
      "rs", "Resample input data to uniform pressure levels (value=interval)",
      false);
58
59
60
61
62
63
64
  commandLine.addOption("rsstart",
                        "Resampling start pressure level - "
                        "Defaults to use input pressure extent",
                        false);
  commandLine.addOption("rsend",
                        "Resampling end pressure level - "
                        "Defaults to use input pressure extent",
65
                        false);
66
67
68
69

  // Ensure required options present
  std::vector<std::string> commandErrors;
  if (!commandLine.validate(commandErrors))
70
  {
71
72
73
74
75
76
    std::cout << "Error in arguments..." << std::endl;
    for (std::string error : commandErrors)
    {
      std::cout << "\t" << error << std::endl;
    }
    std::cout << std::endl;
Purves, Murray's avatar
Purves, Murray committed
77
    commandLine.printParsedLine(std::cout);
78

79
80
    return -1;
  }
81
82
83
84
85

  // Get command line options
  std::string inputCsvPath = commandLine.get<std::string>("i");
  std::string outputArlPath =
      commandLine.get<std::string>("o", inputCsvPath + ".bin");
86
87
88
89
90
91
92
93
  float extent           = commandLine.get<float>("e", 500.0);
  float resolution       = commandLine.get<float>("r", 10.0);
  float centreLat        = commandLine.get<float>("clat");
  float centreLon        = commandLine.get<float>("clon");
  int numberTimesteps    = commandLine.get<int>("n", 1);
  float groundElevation  = commandLine.get<float>("g", 0.f);
  std::string startTime  = commandLine.get<std::string>("t", "1951111917");
  float resampleInterval = commandLine.get<float>("rs", -1.f);
94
95
  float resampleStart    = commandLine.get<float>("rsstart", -1.f);
  float resampleEnd      = commandLine.get<float>("rsend", -1.f);
96
97
98
99
100
101
102
103
104
105
106
107
  // Get the grid size
  int numberGridCells = (int)(extent / resolution);
  if (!numberGridCells % 2 == 1)
  {
    numberGridCells++;
  }
  // Parse the start time
  int year  = from_string(startTime.substr(0, 4), 1951);
  int month = from_string(startTime.substr(4, 2), 11);
  int day   = from_string(startTime.substr(6, 2), 19);
  int hour  = from_string(startTime.substr(8, 2), 17);

108
109
110
111
112
113
  struct tm metTime = {0, 0, 0};
  metTime.tm_year   = year - 1900;
  metTime.tm_mon    = month - 1;
  metTime.tm_mday   = day;
  metTime.tm_hour   = hour;

114
  std::cout << "Creating ARL-formatted met file with parameters:" << std::endl;
115
  std::cout << "  " << extent << " by " << extent << " km grid centred on ("
116
117
118
119
            << centreLat << "," << centreLon << ") with resolution "
            << resolution << " (" << numberGridCells
            << " cells in each direction)" << std::endl;
  std::cout << "  " << numberTimesteps
120
            << " timesteps (of 1 hour each) starting from " << metTime.tm_year
Purves, Murray's avatar
Purves, Murray committed
121
            << "-" << metTime.tm_mon + 1 << "-" << metTime.tm_mday << " at "
122
            << metTime.tm_hour << "00" << std::endl
123
124
125
126
127
128
            << std::endl;

  // If we have correct options, parse input file
  std::cout << "Reading input csv file: " << inputCsvPath << std::endl;

  // Header strings to denote fields
129
  std::string pressureString = "pres", tempString = "temp", relHumString = "rh",
130
131
              dewPtString = "tdew", wDirString = "wdir", wSpdString = "wspd",
              heightString = "zhgt";
132
133
134
  int pressureIndex = -1, tempIndex = -1, relHumIndex = -1, dewPtIndex = -1,
      wDirIndex = -1, wSpdIndex = -1, heightIndex = -1;
  bool usingRelHum = false, usingDewPt = false;
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

  // Open and read the csv
  CSVFile inputCsv(inputCsvPath);
  std::vector<std::vector<std::string>> inputData;
  bool inputReadSuccess = inputCsv.data(inputData);
  if (!inputReadSuccess)
  {
    std::cout << "Failed to read input csv!" << std::endl;
    std::cout
        << "Please ensure the file exists and has the appropriate permissions."
        << std::endl;
    return -2;
  }
  // Ensure there is data in here
  if (inputData.size() < 5)
  {
    std::cout << "Not enough data in input csv!" << std::endl;
    std::cout
        << "There are only " << inputData.size()
        << " lines in this file - at least 5 are required (4 headers + data)"
        << std::endl;
    return -3;
  }

  // Work out which fields are which data points - should be 2nd line (1)
  std::cout << "Finding data fields..." << std::endl;
LEFEBVREJP email's avatar
LEFEBVREJP email committed
161
  for (size_t entry = 0; entry < inputData[1].size(); ++entry)
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  {
    if (inputData[1][entry] == pressureString)
    {
      std::cout << "  Found pressure at index " << entry << std::endl;
      pressureIndex = entry;
    }
    else if (inputData[1][entry] == tempString)
    {
      std::cout << "  Found temperature at index " << entry << std::endl;
      tempIndex = entry;
    }
    else if (inputData[1][entry] == relHumString)
    {
      std::cout << "  Found relative humidity at index " << entry << std::endl;
      relHumIndex = entry;
177
178
179
180
181
182
183
      usingRelHum = true;
    }
    else if (inputData[1][entry] == dewPtString)
    {
      std::cout << "  Found dew point at index " << entry << std::endl;
      dewPtIndex = entry;
      usingDewPt = true;
184
185
186
187
188
189
190
191
192
193
194
    }
    else if (inputData[1][entry] == wDirString)
    {
      std::cout << "  Found wind direction at index " << entry << std::endl;
      wDirIndex = entry;
    }
    else if (inputData[1][entry] == wSpdString)
    {
      std::cout << "  Found wind speed at index " << entry << std::endl;
      wSpdIndex = entry;
    }
195
196
197
198
199
    else if (inputData[1][entry] == heightString)
    {
      std::cout << "  Found height at index " << entry << std::endl;
      heightIndex = entry;
    }
200
201
  }
  // If we are missing one, we can't continue
202
203
204
  if ((pressureIndex == -1) || (tempIndex == -1) ||
      (relHumIndex == -1 && dewPtIndex == -1) || (wDirIndex == -1) ||
      (wSpdIndex == -1))
205
  {
Purves, Murray's avatar
Purves, Murray committed
206
    std::cout << "Missing data fields in input!" << std::endl;
207
208
209
210
211
212
213
214
215
216
    if (pressureIndex == -1)
    {
      std::cout << "  Missing pressure field (denoted by " << pressureString
                << ")";
    }
    if (tempIndex == -1)
    {
      std::cout << "  Missing temperature field (denoted by " << tempString
                << ")";
    }
217
    if (relHumIndex == -1 && dewPtIndex == -1)
218
    {
219
220
221
      std::cout << "  Missing both relative humidity field (denoted by "
                << relHumString << ") and dew point field (denoted by "
                << dewPtString << ") - need at least one of these";
222
223
224
225
226
227
228
229
230
231
232
    }
    if (wDirIndex == -1)
    {
      std::cout << "  Missing wind direction field (denoted by " << wDirString
                << ")";
    }
    if (wSpdIndex == -1)
    {
      std::cout << "  Missing wind speed field (denoted by " << wSpdString
                << ")";
    }
233
234
235
236
    if (heightIndex == -1)
    {
      std::cout << "  Missing height field (denoted by " << heightString << ")";
    }
237
238
239
240
    std::cout << "Please ensure these fields are present and rerun.";
    return -4;
  }

241
242
  // Ensure we aren't using both relative humidity and dew point (duplicate
  // data) Prefer use of relative humidity
243
244
245
246
247
  if (usingRelHum)
  {
    usingDewPt = false;
  }

248
  // Read the data
249
250
  std::vector<float> inputPressures, inputTemps, inputRelHums, inputDewPts,
      inputWDirs, inputWSpds, inputHeights;
251
252
  std::cout << "Data fields found - reading data..." << std::endl;
  float missingValue = -9999.f;
LEFEBVREJP email's avatar
LEFEBVREJP email committed
253
  for (size_t row = 4; row < inputData.size(); ++row)
254
255
256
257
258
259
260
261
  {
    if (inputData[row].size() < inputData[1].size())
    {
      std::cout << "  Warning: this row (" << row + 1
                << ") has less entries than the field names row" << std::endl;
      ;
    }

262
    bool foundPressure = false, foundTemp = false, foundRelHum = false,
263
264
         foundDewPt = false, foundWDir = false, foundWSpd = false,
         foundHeight = false;
LEFEBVREJP email's avatar
LEFEBVREJP email committed
265
    for (size_t entry = 0; entry < inputData[row].size(); ++entry)
266
267
268
269
270
271
272
273
274
275
276
277
278
    {
      float thisValue = from_string(inputData[row][entry], missingValue);

      if (entry == pressureIndex)
      {
        foundPressure = true;
        inputPressures.push_back(thisValue);
      }
      else if (entry == tempIndex)
      {
        foundTemp = true;
        inputTemps.push_back(thisValue);
      }
279
      else if (usingRelHum && entry == relHumIndex)
280
281
282
283
      {
        foundRelHum = true;
        inputRelHums.push_back(thisValue);
      }
284
285
286
287
288
      else if (usingDewPt && entry == dewPtIndex)
      {
        foundDewPt = true;
        inputDewPts.push_back(thisValue);
      }
289
290
291
292
293
294
295
296
297
298
      else if (entry == wSpdIndex)
      {
        foundWSpd = true;
        inputWSpds.push_back(thisValue);
      }
      else if (entry == wDirIndex)
      {
        foundWDir = true;
        inputWDirs.push_back(thisValue);
      }
299
300
301
302
303
      else if (entry == heightIndex)
      {
        foundHeight = true;
        inputHeights.push_back(thisValue);
      }
304
305
306
307
308
309
310
311
312
313
314
315
316
317
    }
    // Add missing data if any of the above weren't found
    if (!foundPressure)
    {
      std::cout << "  Warning: couldn't find pressure in row " << row
                << std::endl;
      inputPressures.push_back(missingValue);
    }
    if (!foundTemp)
    {
      std::cout << "  Warning: couldn't find temperature in row " << row
                << std::endl;
      inputTemps.push_back(missingValue);
    }
318
    if (usingRelHum && !foundRelHum)
319
320
321
322
323
    {
      std::cout << "  Warning: couldn't find relative humidity in row " << row
                << std::endl;
      inputRelHums.push_back(missingValue);
    }
324
325
326
327
328
329
    if (usingDewPt && !foundDewPt)
    {
      std::cout << "  Warning: couldn't find dew point in row " << row
                << std::endl;
      inputDewPts.push_back(missingValue);
    }
330
331
332
333
334
335
336
337
338
339
340
341
    if (!foundWSpd)
    {
      std::cout << "  Warning: couldn't find wind speed in row " << row
                << std::endl;
      inputWSpds.push_back(missingValue);
    }
    if (!foundWDir)
    {
      std::cout << "  Warning: couldn't find wind direction in row " << row
                << std::endl;
      inputWDirs.push_back(missingValue);
    }
342
343
344
345
346
347
    if (!foundHeight)
    {
      std::cout << "  Warning: couldn't find height in row " << row
                << std::endl;
      inputHeights.push_back(missingValue);
    }
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368

    // If pressure value in this row is -9999, remove the row
    //  - pressure entry is required for interpolation & by HYSPLIT
    if (inputPressures.back() == missingValue)
    {
      std::cout << "  Warning: pressure entry in row " << row
                << " is -9999; removing row" << std::endl;
      inputPressures.pop_back();
      inputTemps.pop_back();
      if (usingRelHum)
      {
        inputRelHums.pop_back();
      }
      if (usingDewPt)
      {
        inputDewPts.pop_back();
      }
      inputWSpds.pop_back();
      inputWDirs.pop_back();
      inputHeights.pop_back();
    }
369
370
371
372
  }
  std::cout << "Data read complete: " << inputPressures.size()
            << " entries read." << std::endl;

373
374
375
376
377
378
379
380
381
382
383
384
  // We will need the ground level pressure later
  // Assume it is the first entry in the input
  float groundPressure = inputPressures[0];

  // Interpolate the data to remove missing values
  std::cout << "Interpolating values to remove missing data..." << std::endl;
  // Resample the pressures if required
  if (resampleInterval > 0)
  {
    std::cout << "  Resampling input data using sampling interval: "
              << resampleInterval << " mb" << std::endl;
    std::vector<float> resampledPressures;
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405

    // Calculate the upper and lower bounds of the resampled pressure array
    radix_line("    input start = " << resampleStart
                                    << "; end = " << resampleEnd);
    if (resampleStart < 0)
    {
      resampleStart =
          (float)roundUpInt((int)*std::max_element(std::begin(inputPressures),
                                                   std::end(inputPressures)),
                            resampleInterval);
    }
    if (resampleEnd < 0)
    {
      resampleEnd =
          (float)roundDownInt((int)*std::min_element(std::begin(inputPressures),
                                                     std::end(inputPressures)),
                              resampleInterval);
    }
    std::cout << "    Using start = " << resampleStart
              << "; end = " << resampleEnd << std::endl;

406
407
    // Calculate the resampled pressure values using the interval given
    std::cout << "    Calculating resampled pressure values..." << std::endl;
Purves, Murray's avatar
Purves, Murray committed
408
    float thisPressure;
409
    for (thisPressure = resampleStart; thisPressure > resampleEnd;
410
411
         thisPressure -= resampleInterval)
    {
Purves, Murray's avatar
Purves, Murray committed
412
413
414
415
416
417
      resampledPressures.push_back(thisPressure);
    }
    // Add another entry below the min. pressure
    // - but only if that pressure won't be 0
    if (thisPressure > 0)
    {
418
419
420
421
      resampledPressures.push_back(thisPressure);
    }

    std::cout << "  Temperature:" << std::endl;
422
423
    inputTemps = interpolateToOtherBaseValues(
        inputPressures, resampledPressures, inputTemps, true);
424
425
426
427
    if (usingRelHum)
    {
      std::cout << "  Relative humidity:" << std::endl;
      inputRelHums = interpolateToOtherBaseValues(
428
          inputPressures, resampledPressures, inputRelHums, true);
429
430
431
432
433
    }
    if (usingDewPt)
    {
      std::cout << "  Dew point:" << std::endl;
      inputDewPts = interpolateToOtherBaseValues(
434
          inputPressures, resampledPressures, inputDewPts, true);
435
436
    }
    std::cout << "  Wind speed:" << std::endl;
437
    inputWSpds = interpolateToOtherBaseValues(
438
        inputPressures, resampledPressures, inputWSpds, true);
439
440
    std::cout << "  Wind direction:" << std::endl;
    inputWDirs = interpolateToOtherBaseValues(
441
        inputPressures, resampledPressures, inputWDirs, true, true);
442
443
    std::cout << "  Height:" << std::endl;
    inputHeights = interpolateToOtherBaseValues(
444
        inputPressures, resampledPressures, inputHeights, true);
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459

    std::cout << "  Replacing pressure data with resampled array" << std::endl;
    radix_block(
        radix("    Initial pressure array: ");
        for (float p
             : inputPressures) { radix(" " << p); } radix_line("");
        radix("    Resampled pressure array: ");
        for (float p
             : resampledPressures) { radix(" " << p); } radix_line(""););
    inputPressures = resampledPressures;
  }
  else
  {
    std::cout << "  Not resampling data" << std::endl;
    std::cout << "  Temperature:" << std::endl;
460
    inputTemps = interpolateValues(inputPressures, inputTemps, true);
461
462
463
    if (usingRelHum)
    {
      std::cout << "  Relative humidity:" << std::endl;
464
      inputRelHums = interpolateValues(inputPressures, inputRelHums, true);
465
466
467
468
    }
    if (usingDewPt)
    {
      std::cout << "  Dew point:" << std::endl;
469
      inputDewPts = interpolateValues(inputPressures, inputDewPts, true);
470
471
    }
    std::cout << "  Wind speed:" << std::endl;
472
    inputWSpds = interpolateValues(inputPressures, inputWSpds, false);
473
    std::cout << "  Wind direction:" << std::endl;
474
    inputWDirs = interpolateValues(inputPressures, inputWDirs, false, true);
475
    std::cout << "  Height:" << std::endl;
476
    inputHeights = interpolateValues(inputPressures, inputHeights, true);
477
478
479
  }
  std::cout << "  Interpolation complete." << std::endl;

480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
  // Check if the grid is of appropriate size - with too small grid/too many
  // vertical levels, the index header can become oversized leading to errors
  int indexBase = 108, levelBase = 48;
  int indexSize  = indexBase + (inputPressures.size() * levelBase);
  int recordSize = numberGridCells * numberGridCells;
  if (indexSize > recordSize)
  {
    // The index header will be too big - instruct user to either reduce number
    // of vertical levels or increase number of grid cells
    std::cout << "Error: index header will be too large!" << std::endl
              << "  The size of the index header in the ARL-formatted output ("
              << indexBase << " + number of vertical levels * " << levelBase
              << " = " << indexSize
              << ") will be larger than the size of each record (number of "
                 "grid cells ^ 2 = "
              << recordSize << ")." << std::endl;
    std::cout << "  Please either increase the number of horizontal grid cells "
                 "(increase extent/resolution) or reduce the number of input "
                 "vertical levels."
              << std::endl;
    throw std::exception();
  }

503
504
505
506
507
508
509
510
511
  // Convert the data into ARL format
  std::cout << "Converting data to ARL format..." << std::endl;
  ARLDataStream outputStream(outputArlPath, std::ios::out);
  for (int timestep = 0; timestep < numberTimesteps; ++timestep)
  {
    // Write index header section
    std::cout << "    Writing index headers for timestep " << timestep << "..."
              << std::endl;
    ARLRecordHeader thisRecordHeader;
512
513
514
515
    thisRecordHeader.year  = metTime.tm_year;
    thisRecordHeader.month = metTime.tm_mon + 1;
    thisRecordHeader.day   = metTime.tm_mday;
    thisRecordHeader.hour  = metTime.tm_hour;
516
517
518
519
520
521
522
523
    thisRecordHeader.ic    = 0;
    thisRecordHeader.il    = 0;
    thisRecordHeader.cgrid = "99";
    thisRecordHeader.kvar  = "INDX";
    thisRecordHeader.nexp  = 0;
    thisRecordHeader.prec  = 0.f;
    thisRecordHeader.var1  = 0.f;
    ARLIndexHeader thisIndexHeader;
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
    thisIndexHeader.model_id = "NFDB";
    thisIndexHeader.icx      = 0;
    thisIndexHeader.mn       = 0;
    thisIndexHeader.pole_lat = centreLat;
    thisIndexHeader.pole_lon = centreLon;
    thisIndexHeader.ref_lat  = centreLat;
    thisIndexHeader.ref_lon  = centreLon;
    thisIndexHeader.size     = resolution;
    thisIndexHeader.orient   = 0.f;
    thisIndexHeader.tang_lat = centreLat;
    thisIndexHeader.sync_xp  = (numberGridCells + 1) / 2;
    thisIndexHeader.sync_yp  = (numberGridCells + 1) / 2;
    thisIndexHeader.sync_lat = centreLat;
    thisIndexHeader.sync_lon = centreLon;
    thisIndexHeader.dummy    = 0.f;
    thisIndexHeader.nx       = numberGridCells;
    thisIndexHeader.ny       = numberGridCells;
LEFEBVREJP email's avatar
LEFEBVREJP email committed
541
    thisIndexHeader.nz       = int(inputPressures.size());
542
543
544
545
546
    thisIndexHeader.z_flag   = 2;
    thisIndexHeader.lenh     = numberGridCells * numberGridCells;
    thisIndexHeader.levels   = inputPressures;
    thisIndexHeader.num_vars_at_levels =
        std::vector<int>(inputPressures.size(), 5);
547
548
549
550
    std::vector<std::string> surfaceVarNames = {"PRSS", "TEMP", "RELH", "UWND",
                                                "VWND"};
    std::vector<std::string> varNames        = {"HGTS", "TEMP", "RELH", "UWND",
                                         "VWND"};
551

552
    thisIndexHeader.var_names.push_back(surfaceVarNames);
553
554
    for (size_t level = 1; level < inputPressures.size(); ++level)
    {
555
      thisIndexHeader.var_names.push_back(varNames);
556
557
558
559
    }
    std::vector<int> checkSums = std::vector<int>(5, 0);
    for (size_t level = 0; level < inputPressures.size(); ++level)
    {
560
      thisIndexHeader.check_sums.push_back(checkSums);
561
    }
562

563
564
565
566
567
568
569
570
571
572
573
    // Write the headers
    outputStream.write_record_header(thisRecordHeader);
    outputStream.write_index_header(thisRecordHeader, thisIndexHeader);
    std::cout << "    Headers written" << std::endl;

    // Write meteorological variables for each level
    for (int level = 0; level < inputPressures.size(); ++level)
    {
      std::cout << "    Writing meteorological data for timestep " << timestep
                << ", level " << level << "..." << std::endl;

574
575
576
      thisRecordHeader.il = level;

      // Write pressure/height variables
577
      {
578
579
580
        if (level == 0)
        {
          thisRecordHeader.kvar = "PRSS";
581
          thisRecordHeader.var1 = groundPressure;
582
          std::vector<std::vector<float>> thisData(
583
584
              numberGridCells,
              std::vector<float>(numberGridCells, groundPressure));
585
586
587
588
589
          std::cout << "      pressure...";
          outputStream.write_record_header(thisRecordHeader);
          outputStream.write_record(thisRecordHeader, thisIndexHeader,
                                    thisData);
          std::cout << "written" << std::endl;
590
591
592
        }
        else
        {
593
          thisRecordHeader.kvar = "HGTS";
594
595
          float thisHeight      = inputHeights[level] + groundElevation;
          thisRecordHeader.var1 = thisHeight;
596
          std::vector<std::vector<float>> thisData(
597
              numberGridCells, std::vector<float>(numberGridCells, thisHeight));
598
599
600
601
602
          std::cout << "      height...";
          outputStream.write_record_header(thisRecordHeader);
          outputStream.write_record(thisRecordHeader, thisIndexHeader,
                                    thisData);
          std::cout << "written" << std::endl;
603
        }
604
605
606
      }

      // Write temperature levels
607
608
      // First convert temperature from Celsius to Kelvin
      inputTemps[level] = inputTemps[level] - ABS_ZERO_CELSIUS;
609
610
611
612
613
614
615
616
617
618
619
620
      {
        thisRecordHeader.kvar = "TEMP";
        thisRecordHeader.var1 = inputTemps[level];
        std::vector<std::vector<float>> thisData(
            numberGridCells,
            std::vector<float>(numberGridCells, inputTemps[level]));
        std::cout << "      temperature...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }

621
      // Write relative humidity levels
622
      {
623
624
625
626
627
        double thisVar = 0.0;
        if (usingRelHum)
        {
          thisVar = inputRelHums[level];
        }
628
629
        else if (usingDewPt)
        {
630
631
          thisVar =
              dewPointToRelativeHumidity(inputDewPts[level], inputTemps[level]);
632
        }
633

634
635
        thisRecordHeader.kvar = "RELH";
        thisRecordHeader.var1 = thisVar;
636
        std::vector<std::vector<float>> thisData(
637
            numberGridCells, std::vector<float>(numberGridCells, thisVar));
638
639
640
        std::cout << "      relative humidity...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
641
642
643
644
645
646
        std::cout << "written" << std::endl;
      }

      // Write wind data
      // Need to calculate wind u and v components from direction/speed
      float thisWindU =
Purves, Murray's avatar
Purves, Murray committed
647
648
          (inputWSpds[level] * 0.5144444444) *
          sin(toRadians(fmod((inputWDirs[level] + 180.0), 360.0)));
649
      float thisWindV =
Purves, Murray's avatar
Purves, Murray committed
650
651
          (inputWSpds[level] * 0.5144444444) *
          cos(toRadians(fmod((inputWDirs[level] + 180.0), 360.0)));
652
      std::cout << "      Initial wspd = " << inputWSpds[level]
653
654
                << " knots, wdir: " << inputWDirs[level] << " degrees"
                << std::endl;
655
      std::cout << "      Converted wind components: u = " << thisWindU
656
                << " m/s, v = " << thisWindV << "m/s" << std::endl;
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
      {
        thisRecordHeader.kvar = "UWND";
        thisRecordHeader.var1 = thisWindU;
        std::vector<std::vector<float>> thisData(
            numberGridCells, std::vector<float>(numberGridCells, thisWindU));
        std::cout << "      wind (u component)...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }
      {
        thisRecordHeader.kvar = "VWND";
        thisRecordHeader.var1 = thisWindV;
        std::vector<std::vector<float>> thisData(
            numberGridCells, std::vector<float>(numberGridCells, thisWindV));
        std::cout << "      wind (v component)...";
        outputStream.write_record_header(thisRecordHeader);
        outputStream.write_record(thisRecordHeader, thisIndexHeader, thisData);
        std::cout << "written" << std::endl;
      }
    }
    std::cout << " Written timestep " << timestep << std::endl;
679
680

    // Add an hour to time for next timestep
681
682
683
#ifdef _WIN32
    metTime.tm_year = metTime.tm_year + 400;
#endif  // _WIN32
684
    addHour(&metTime, 1);
LEFEBVREJP email's avatar
LEFEBVREJP email committed
685
686
    printf("%d%02d%02d:%02d\n", metTime.tm_year, (metTime.tm_mon + 1),
           metTime.tm_mday, metTime.tm_hour);
687
688
689
#ifdef _WIN32
    metTime.tm_year = metTime.tm_year - 400;
#endif
690
691
692
693
  }

  std::cout << "File write complete!" << std::endl;

694
695
  return 0;
}