WSL/SLF GitLab Repository

ARPSIO.cc 19.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
18
19
#include <meteoio/plugins/ARPSIO.h>
#include <meteoio/IOUtils.h>
20
#include <meteoio/FileUtils.h>
21
#include <meteoio/IOExceptions.h>
22
#include <meteoio/meteoLaws/Meteoconst.h> //for PI
23
24
#include <meteoio/MathOptim.h>

25
26
#include <string.h>
#include <algorithm>
27
#include <fstream>
28

29
30
31
32
33
34
using namespace std;

namespace mio {
/**
 * @page arps ARPSIO
 * @section arps_format Format
35
36
37
38
 * This is for reading grid data in the ARPS grid format (it transparently supports both true ARPS ascii grids and 
 * grids modified by the ARPSGRID utility). DEM reading works well while reading meteo parameters might be a rough ride 
 * (since ARPS files do not always contain a consistent set of meteo fields).
 * 
39
40
 * It is possible to manually print the list of available variables as well as extract a whole section, for a given variable 
 * with the following one line scripts, in order to ease debugging:
41
 * @code
42
 * awk '/^ [a-z]/{print $0}'  {arps_file}
43
44
 * awk '/^ [a-z]/{isVar=0} /^ radsw/{isVar=1;next} {if(isVar) print $0}' {arps_file}
 * @endcode
45
46
 *
 * @section arps_units Units
47
 * All units are assumed to be MKSA.
48
49
50
51
52
53
 *
 * @section arps_keywords Keywords
 * This plugin uses the following keywords:
 * - COORDSYS: coordinate system (see Coords); [Input] and [Output] section
 * - COORDPARAM: extra coordinates parameters (see Coords); [Input] and [Output] section
 * - DEMFILE: path and file containing the DEM; [Input] section
54
55
56
 * - ARPS_XCOORD: x coordinate of the lower left corner of the grids; [Input] section
 * - ARPS_YCOORD: y coordinate of the lower left corner of the grids; [Input] section
 * - GRID2DPATH: path to the input directory where to find the arps files to be read as grids; [Input] section
57
58
 * - GRID3DPATH: path to the input directory where to find the arps 3D files to be read as grids; [Input] section
 * - ARPS_EXT: arps file extension, or <i>none</i> for no file extension (default: .asc)
59
60
61
 */

const double ARPSIO::plugin_nodata = -999.; //plugin specific nodata value
62
const char* ARPSIO::default_ext=".asc"; //filename extension
63

64
ARPSIO::ARPSIO(const std::string& configfile)
65
        : cfg(configfile),
66
          coordin(), coordinparam(), coordout(), coordoutparam(),
67
          grid2dpath_in(), grid3dpath_in(), ext(default_ext), dimx(0), dimy(0), dimz(0), cellsize(0.),
68
          xcoord(IOUtils::nodata), ycoord(IOUtils::nodata), zcoord(), is_true_arps(true)
69
70
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
71
	setOptions();
72
73
}

74
ARPSIO::ARPSIO(const Config& cfgreader)
75
        : cfg(cfgreader),
76
          coordin(), coordinparam(), coordout(), coordoutparam(),
77
          grid2dpath_in(), grid3dpath_in(), ext(default_ext), dimx(0), dimy(0), dimz(0), cellsize(0.),
78
          xcoord(IOUtils::nodata), ycoord(IOUtils::nodata), zcoord(), is_true_arps(true)
79
80
{
	IOUtils::getProjectionParameters(cfg, coordin, coordinparam, coordout, coordoutparam);
81
82
83
	setOptions();
}

84
ARPSIO& ARPSIO::operator=(const ARPSIO& source) {
85
	if (this != &source) {
86
87
88
89
90
		coordin = source.coordin;
		coordinparam = source.coordinparam;
		coordout = source.coordout;
		coordoutparam = source.coordoutparam;
		grid2dpath_in = source. grid2dpath_in;
91
		grid3dpath_in = source.grid3dpath_in;
92
93
94
95
96
97
98
99
100
101
102
103
104
		ext = source.ext;
		dimx = source.dimx;
		dimy = source.dimy;
		dimz = source.dimz;
		cellsize = source.cellsize;
		xcoord = source.xcoord;
		ycoord = source.ycoord;
		zcoord = source.zcoord;
		is_true_arps = source.is_true_arps;
	}
	return *this;
}

105
106
void ARPSIO::setOptions()
{
107
108
	const string grid_in = cfg.get("GRID2D", "Input", IOUtils::nothrow);
	if (grid_in == "ARPS") { //keep it synchronized with IOHandler.cc for plugin mapping!!
109
110
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
	}
111
112
113
114
115
	
	const string grid3d_in = cfg.get("GRID3D", "Input", IOUtils::nothrow);
	if (grid3d_in == "ARPS") { //keep it synchronized with IOHandler.cc for plugin mapping!!
		cfg.getValue("GRID3DPATH", "Input", grid3dpath_in);
	}
116

117
118
	cfg.getValue("ARPS_XCOORD", "Input", xcoord, IOUtils::nothrow);
	cfg.getValue("ARPS_YCOORD", "Input", ycoord, IOUtils::nothrow);
119

120
	//default value has been set in constructor
121
	cfg.getValue("ARPS_EXT", "Input", ext, IOUtils::nothrow);
122
	if (ext=="none") ext.clear();
123
124
}

125
void ARPSIO::listFields(const std::string& i_name)
126
{
127
128
129
	const size_t pos = i_name.find_last_of(":");//a specific parameter can be provided as {filename}:{parameter}
	const std::string filename = (pos!=IOUtils::npos)? grid2dpath_in +"/" + i_name.substr(0, pos) : grid2dpath_in +"/" + i_name;
	
130
131
	FILE *fin;
	openGridFile(fin, filename);
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
	const std::string dem_marker = (is_true_arps)? "zp coordinat" : "zp_coordinat";
	moveToMarker(fin, filename, dem_marker);
	
	char dummy[ARPS_MAX_LINE_LENGTH];
	int nb_elems = 0;
	do {
		nb_elems = fscanf(fin," %[^\t\n] ",dummy); //HACK: possible buffer overflow
		if (dummy[0]>='A' && dummy[0]<='z') { //white spaces have been skipped above
			std::string tmp( dummy );
			IOUtils::trim(tmp);
			std::cout << "Found '" << tmp << "'\n";
		}
	} while (!feof(fin) && nb_elems!=0);
	
	fclose(fin);
}
148

149
150
151
152
153
154
155
156
157
158
159
void ARPSIO::read2DGrid(Grid2DObject& grid_out, const std::string& i_name)
{
	const size_t pos = i_name.find_last_of(":");//a specific parameter can be provided as {filename}:{parameter}
	const std::string filename = (pos!=IOUtils::npos)? grid2dpath_in +"/" + i_name.substr(0, pos) : grid2dpath_in +"/" + i_name;
	if (pos==IOUtils::npos) { //TODO: read by default the first data grid that is found?
		listFields(i_name);
		throw InvalidArgumentException("Please provide the parameter that has to be read!", AT);
	}
	
	const std::string param_str = (pos!=IOUtils::npos)?  i_name.substr(pos+1) : "";
	const size_t param = MeteoGrids::getParameterIndex( param_str );
160
	if (param==IOUtils::npos)
161
162
163
164
165
166
167
168
169
170
171
		throw InvalidArgumentException("Invalid MeteoGrids Parameter requested: '"+param_str+"'", AT);
	
	FILE *fin;
	openGridFile(fin, filename);
	read2DGrid_internal(fin, filename, grid_out, static_cast<MeteoGrids::Parameters>(param));
	if (grid_out.empty()) {
		ostringstream ss;
		ss << "No suitable data found for parameter " << MeteoGrids::getParameterName(param) << " ";
		ss << " in file \"" << filename << "\"";
		throw NoDataException(ss.str(), AT);
	}
172
	fclose(fin);
173
174
}

175
176
177
178
void ARPSIO::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
{
	std::string date_str = date.toString(Date::ISO);
	std::replace( date_str.begin(), date_str.end(), ':', '.');
179
180
181
	const std::string filename = grid2dpath_in + "/" + date_str + ext;
	FILE *fin;
	openGridFile(fin, filename);
182
183
184
185
186
187
188
	read2DGrid_internal(fin, filename, grid_out, parameter);
	if (grid_out.empty()) {
		ostringstream ss;
		ss << "No suitable data found for parameter " << MeteoGrids::getParameterName(parameter) << " ";
		ss << "at time step " << date.toString(Date::ISO) << " in file \"" << filename << "\"";
		throw NoDataException(ss.str(), AT);
	}
189

190
191
192
193
194
195
196
197
198
	fclose(fin);
}

void ARPSIO::read2DGrid_internal(FILE* &fin, const std::string& filename, Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter)
{
	//extra variables: 
	//kmh / kmv: horiz./vert. turbulent mixing coeff. m^2/s
	//tke: turbulent kinetic energy. (m/s)^2
	
199
	//Radiation parameters
200
	if (parameter==MeteoGrids::ISWR) readGridLayer(fin, filename, "radsw", 2, grid_out);
201
	if (parameter==MeteoGrids::RSWR) {
202
		Grid2DObject net;
203
204
		readGridLayer(fin, filename, "radsw", 2, grid_out);
		readGridLayer(fin, filename, "radswnet", 2, net);
205
206
		grid_out.grid2D -= net.grid2D;
	}
207
	if (parameter==MeteoGrids::ILWR) readGridLayer(fin, filename, "radlwin", 2, grid_out);
208
	if (parameter==MeteoGrids::ALB) {
209
		Grid2DObject rswr, iswr;
210
211
		readGridLayer(fin, filename, "radsw", 2, iswr);
		readGridLayer(fin, filename, "radswnet", 2, grid_out); //net radiation
212
213
214
215
216
		rswr.grid2D = iswr.grid2D - grid_out.grid2D;
		grid_out.grid2D = iswr.grid2D/rswr.grid2D;
	}

	//Wind grids
217
218
219
	if (parameter==MeteoGrids::U) readGridLayer(fin, filename, "u", 2, grid_out);
	if (parameter==MeteoGrids::V) readGridLayer(fin, filename, "v", 2, grid_out);
	if (parameter==MeteoGrids::W) readGridLayer(fin, filename, "w", 2, grid_out);
220
	if (parameter==MeteoGrids::VW) {
221
		Grid2DObject V;
222
223
		readGridLayer(fin, filename, "u", 2, grid_out); //U
		readGridLayer(fin, filename, "v", 2, V);
224
225
		for (size_t jj=0; jj<grid_out.getNy(); jj++) {
			for (size_t ii=0; ii<grid_out.getNx(); ii++) {
226
				grid_out(ii,jj) = sqrt( Optim::pow2(grid_out(ii,jj)) + Optim::pow2(V(ii,jj)) );
227
228
229
			}
		}
	}
230
	if (parameter==MeteoGrids::DW) {
231
		Grid2DObject V;
232
233
		readGridLayer(fin, filename, "u", 2, grid_out); //U
		readGridLayer(fin, filename, "v", 2, V);
234
235
		for (size_t jj=0; jj<grid_out.getNy(); jj++) {
			for (size_t ii=0; ii<grid_out.getNx(); ii++) {
236
				grid_out(ii,jj) = fmod( atan2( grid_out(ii,jj), V(ii,jj) ) * Cst::to_deg + 360., 360.); // turn into degrees [0;360)
237
238
239
240
241
			}
		}
	}

	//Basic meteo parameters
242
	if (parameter==MeteoGrids::P) readGridLayer(fin, filename, "p", 2, grid_out);
243
	if (parameter==MeteoGrids::TSG) readGridLayer(fin, filename, "tsoil", 2, grid_out); //or is it tss for us?
244
	/*if (parameter==MeteoGrids::RH) {
245
		//const double epsilon = Cst::gaz_constant_dry_air / Cst::gaz_constant_water_vapor;
246
		readGridLayer(filename, "qv", 2, grid_out); //water vapor mixing ratio
247
		//Atmosphere::vaporSaturationPressure(T);
248
		//TODO: compute relative humidity out of it!
249
		//through potential temperature "pt" -> local temperature?
250
251
252
	}*/

	//Hydrological parameters
253
	if (parameter==MeteoGrids::HS) readGridLayer(fin, filename, "snowdpth", 2, grid_out);
254
	if (parameter==MeteoGrids::PSUM) {
255
		readGridLayer(fin, filename, "prcrate1", 2, grid_out); //in kg/m^2/s
256
257
258
259
		grid_out.grid2D *= 3600.; //we need kg/m^2/h
	}

	//DEM
260
	const std::string dem_marker = (is_true_arps)? "zp coordinat" : "zp_coordinat";
261
	if (parameter==MeteoGrids::DEM) readGridLayer(fin, filename, dem_marker, 2, grid_out);
262
	if (parameter==MeteoGrids::SLOPE) {
263
264
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
265
		readGridLayer(fin, filename, dem_marker, 2, dem);
266
		dem.update();
267
		grid_out.set(dem.cellsize, dem.llcorner, dem.slope);
268
	}
269
	if (parameter==MeteoGrids::AZI) {
270
271
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
272
		readGridLayer(fin, filename, dem_marker, 2, dem);
273
		dem.update();
274
		grid_out.set(dem.cellsize, dem.llcorner, dem.azi);
275
276
	}

277
	rewind(fin);
278
279
}

280
void ARPSIO::read3DGrid(Grid3DObject& grid_out, const std::string& i_name)
281
{
282
	const size_t pos = i_name.find_last_of(":");//a specific parameter can be provided as {filename}:{parameter}
283
284
	const std::string filename = (pos!=std::string::npos)? grid3dpath_in +"/" + i_name.substr(0, pos) : grid3dpath_in +"/" + i_name;
	if (pos==std::string::npos) { //TODO: read by default the first data grid that is found?
285
286
287
288
		listFields(i_name);
		throw InvalidArgumentException("Please provide the parameter that has to be read!", AT);
	}
	
289
	std::string parameter = (pos!=std::string::npos)?  i_name.substr(pos+1) : "";
290
291
292
293
	if (parameter=="DEM") { //this is called damage control... this is so ugly...
		parameter = (is_true_arps)? "zp coordinat" : "zp_coordinat";
	}
	
294
295
	FILE *fin;
	openGridFile(fin, filename);
296
297

	//resize the grid just in case
298
299
300
301
	Coords llcorner(coordin, coordinparam);
	llcorner.setXY(xcoord, ycoord, IOUtils::nodata);
	grid_out.set(dimx, dimy, dimz, cellsize, llcorner);
	grid_out.z = zcoord;
302
303

	// Read until the parameter is found
304
	moveToMarker(fin, filename, parameter);
305
	
306
	//read the data we are interested in
307
	for (size_t iz = 0; iz < dimz; iz++) {
308
		for (size_t iy = 0; iy < dimy; iy++) {
309
			for (size_t ix = 0; ix < dimx; ix++) {
310
				double tmp;
311
				if (fscanf(fin," %16lf%*[\n]",&tmp)==1) {
312
313
					grid_out.grid3D(ix,iy,iz) = tmp;
				} else {
314
315
					fclose(fin);
					throw InvalidFormatException("Failure in reading 3D grid in file "+filename, AT);
316
				}
317
318
319
320
			}
		}
	}

321
	fclose(fin);
322
323
324
325
}

void ARPSIO::readDEM(DEMObject& dem_out)
{
326
327
328
	const std::string filename = cfg.get("DEMFILE", "Input");
	FILE *fin;
	openGridFile(fin, filename);
329
	read2DGrid_internal(fin, filename, dem_out, MeteoGrids::DEM);
330
	fclose(fin);
331
332
}

333
void ARPSIO::initializeGRIDARPS(FILE* &fin, const std::string& filename)
334
335
{
	//go to read the sizes
336
	moveToMarker(fin, filename, "nnx");
337
	//finish reading the line and move to the next one
338
	if (fscanf(fin,"%*[^\n]")!=0) {
339
		fclose(fin);
340
341
		throw InvalidFormatException("Error in file format of file "+filename, AT);
	}
342
	if (fscanf(fin," %u %u %u \n",&dimx,&dimy,&dimz)!=3) {
343
		fclose(fin);
344
345
346
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
347
		fclose(fin);
348
349
350
351
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

	//initializing cell size
352
	moveToMarker(fin, filename, "x_coordinate");
353
	double v1, v2;
354
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
355
		fclose(fin);
356
357
358
		throw InvalidFormatException("Can not read first two x coordinates from file "+filename, AT);
	}
	const double cellsize_x = v2 - v1;
359
	moveToMarker(fin, filename, "y_coordinate");
360
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
361
		fclose(fin);
362
363
364
		throw InvalidFormatException("Can not read first two y coordinates from file "+filename, AT);
	}
	const double cellsize_y = v2 - v1;
365
	if (cellsize_x!=cellsize_y) {
366
		fclose(fin);
367
368
369
		throw InvalidFormatException("Only square cells currently supported! Non compliance in file "+filename, AT);
	}
	cellsize = cellsize_y;
370
371
372

	//HACK: zcoords must be read from zp. But they are NOT constant...

373
374
}

375
void ARPSIO::initializeTrueARPS(FILE* &fin, const std::string& filename, const char curr_line[ARPS_MAX_LINE_LENGTH])
376
377
378
{
	//go to read the sizes
	if (sscanf(curr_line," nx = %u, ny = %u, nz = %u ",&dimx,&dimy,&dimz)!=3) {
379
		fclose(fin);
380
381
382
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
383
		fclose(fin);
384
385
386
387
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

	//initializing cell size
388
	moveToMarker(fin, filename, "x coordinate");
389
	double v1, v2;
390
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
391
		fclose(fin);
392
393
394
		throw InvalidFormatException("Can not read first two x coordinates from file "+filename, AT);
	}
	const double cellsize_x = v2 - v1;
395
	moveToMarker(fin, filename, "y coordinate");
396
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
397
		fclose(fin);
398
399
400
		throw InvalidFormatException("Can not read first two y coordinates from file "+filename, AT);
	}
	const double cellsize_y = v2 - v1;
401
	if (cellsize_x!=cellsize_y) {
402
		fclose(fin);
403
404
405
		throw InvalidFormatException("Only square cells currently supported! Non compliance in file "+filename, AT);
	}
	cellsize = cellsize_y;
406

407
	moveToMarker(fin, filename, "z coordinate");
408
409
410
	while (fscanf(fin,"%lg",&v1)==1) {
		zcoord.push_back( v1 );
	}
411
	if (zcoord.size()!=dimz) {
412
		ostringstream ss;
413
		ss << "Expected " << dimz << " z coordinates in file \""+filename+"\", found " << zcoord.size();
414
		fclose(fin);
415
416
		throw InvalidFormatException(ss.str(), AT);
	}
417
418
}

419
void ARPSIO::openGridFile(FILE* &fin, const std::string& filename)
420
{
421
	if (!FileUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
422
	if ((fin=fopen(filename.c_str(),"r")) == NULL) {
423
		throw AccessException("Can not open file "+filename, AT);
424
425
426
427
	}

	//identify if the file is an original arps file or a file modified by ARPSGRID
	char dummy[ARPS_MAX_LINE_LENGTH];
428
	for (unsigned char j=0; j<5; j++) {
429
		//the first easy difference in the structure happens at line 5
430
		if (fgets(dummy,ARPS_MAX_STRING_LENGTH,fin)==NULL) {
431
			fclose(fin);
432
433
			throw InvalidFormatException("Fail to read header lines of file "+filename, AT);
		}
434
	}
435
	zcoord.clear(); //reset the zcoord
436
	unsigned int v1;
437
438
439
	if (sscanf(dummy," nx = %u, ny = ", &v1)<1) {
		//this is an ASCII file modified by ARPSGRID
		is_true_arps=false;
440
		initializeGRIDARPS(fin, filename);
441
442
	} else {
		//this is a true ARPS file
443
		initializeTrueARPS(fin, filename, dummy);
444
	}
445
446
447
	//set xcoord, ycoord to a default value if not set by the user
	if (xcoord==IOUtils::nodata) xcoord = -cellsize;
	if (ycoord==IOUtils::nodata) ycoord = -cellsize;
448
449
450
451
452

	//come back to the begining of the file
	rewind(fin);
}

453
454
455
456
457
/** @brief Skip the given number of layers from the ARPS file
 * Usually, the first layer for a given parameter is purely numeric and therefore should
 * be skipped.
 * @param fin         file pointer
 * @param filename  file name to use for error messages
458
 * @param layer     Index of the layer to skip to (1 to dimz)
459
*/
460
void ARPSIO::skipToLayer(FILE* &fin, const std::string& filename, const unsigned int& layers) const
461
462
463
464
465
466
467
468
469
470
471
472
473
{
	if (layers>1) {
		double tmp;
		const size_t jmax=dimx*dimy*(layers-1);
		for (size_t j = 0; j < jmax; j++) {
			if (fscanf(fin," %16lf%*[\n]",&tmp)==EOF) {
				fclose(fin);
				throw InvalidFormatException("Fail to skip data layers in file "+filename, AT);
			}
		}
	}
}

474
475
476
477
478
479
480
481
482
483
484
/** @brief Read a specific layer for a given parameter from the ARPS file
 * @param parameter The parameter to extract. This could be any of the following:
 *        - x_coordinate for getting the X coordinates of the mesh
 *        - y_coordinate for getting the Y coordinates of the mesh
 *        - zp_coordinat for getting the Z coordinates of the mesh
 *        - u for getting the u component of the wind field
 *        - v for getting the v component of the wind field
 *        - w for getting the w component of the wind field
 * @param layer     Index of the layer to extract (1 to dimz)
 * @param grid      [out] grid containing the values. The grid will be resized if necessary.
*/
485
void ARPSIO::readGridLayer(FILE* &fin, const std::string& filename, const std::string& parameter, const unsigned int& layer, Grid2DObject& grid)
486
{
487
	if (layer<1 || layer>dimz) {
488
		fclose(fin);
489
		ostringstream tmp;
490
		tmp << "Layer " << layer << " does not exist in ARPS file " << filename << " (nr layers=" << dimz << ")";
491
492
493
494
495
496
497
498
499
		throw IndexOutOfBoundsException(tmp.str(), AT);
	}

	//resize the grid just in case
	Coords llcorner(coordin, coordinparam);
	llcorner.setXY(xcoord, ycoord, IOUtils::nodata);
	grid.set(dimx, dimy, cellsize, llcorner);

	// Read until the parameter is found
500
	moveToMarker(fin, filename, parameter);
501
502

	// move to the begining of the layer of interest
503
	skipToLayer(fin, filename, layer);
504
505

	//read the data we are interested in
506
507
	for (size_t iy = 0; iy < dimy; iy++) {
		for (size_t ix = 0; ix < dimx; ix++) {
508
			double tmp;
509
510
			const int readPts = fscanf(fin," %16lf%*[\n]",&tmp);
			if (readPts==1) {
511
				grid(ix,iy) = tmp;
512
			} else {
513
				char dummy[ARPS_MAX_LINE_LENGTH];
514
				const int status = fscanf(fin,"%s",dummy);
515
				fclose(fin);
516
517
518
519
				string msg = "Fail to read data layer for parameter '"+parameter+"' in file '"+filename+"'";
				if (status>0) msg += ", instead read: '"+string(dummy)+"'";
				
				throw InvalidFormatException(msg, AT);
520
			}
521
522
523
524
		}
	}
}

525
void ARPSIO::moveToMarker(FILE* &fin, const std::string& filename, const std::string& marker)
526
527
528
{
	char dummy[ARPS_MAX_LINE_LENGTH];
	do {
529
530
531
532
533
		const int nb_elems = fscanf(fin," %[^\t\n] ",dummy); //HACK: possible buffer overflow
		if (nb_elems==0) {
			fclose(fin);
			throw InvalidFormatException("Matching failure in file "+filename+" when looking for "+marker, AT);
		}
534
535
536
		if (dummy[0]>='A' && dummy[0]<='z') { //this is a "marker" line
			if (strncmp(dummy,marker.c_str(), marker.size()) == 0) return;
		}
537
	} while (!feof(fin));
538
	
539
540
	fclose(fin);
	throw InvalidFormatException("End of file "+filename+" should NOT have been reached when looking for "+marker, AT);
541
542
543
}

} //namespace