WSL/SLF GitLab Repository

ARPSIO.cc 19.8 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
#include <meteoio/meteoLaws/Meteoconst.h> //for PI
19
20
#include <meteoio/MathOptim.h>

21
22
23
#include <string.h>
#include <algorithm>

24
#include <meteoio/plugins/ARPSIO.h>
25
26
27
28
29
30
31

using namespace std;

namespace mio {
/**
 * @page arps ARPSIO
 * @section arps_format Format
32
33
34
35
 * 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).
 * 
36
37
 * 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:
38
 * @code
39
 * awk '/^ [a-z]/{print $0}'  {arps_file}
40
41
 * awk '/^ [a-z]/{isVar=0} /^ radsw/{isVar=1;next} {if(isVar) print $0}' {arps_file}
 * @endcode
42
43
 *
 * @section arps_units Units
44
 * All units are assumed to be MKSA.
45
46
47
48
49
50
 *
 * @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
51
52
53
 * - 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
54
55
 * - 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)
56
57
58
 */

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

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

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

81
ARPSIO& ARPSIO::operator=(const ARPSIO& source) {
82
	if (this != &source) {
83
84
85
86
87
		coordin = source.coordin;
		coordinparam = source.coordinparam;
		coordout = source.coordout;
		coordoutparam = source.coordoutparam;
		grid2dpath_in = source. grid2dpath_in;
88
		grid3dpath_in = source.grid3dpath_in;
89
90
91
92
93
94
95
96
97
98
99
100
101
		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;
}

102
103
void ARPSIO::setOptions()
{
104
105
	const string grid_in = cfg.get("GRID2D", "Input", IOUtils::nothrow);
	if (grid_in == "ARPS") { //keep it synchronized with IOHandler.cc for plugin mapping!!
106
107
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
	}
108
109
110
111
112
	
	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);
	}
113

114
115
	cfg.getValue("ARPS_XCOORD", "Input", xcoord, IOUtils::nothrow);
	cfg.getValue("ARPS_YCOORD", "Input", ycoord, IOUtils::nothrow);
116

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

122
void ARPSIO::listFields(const std::string& i_name)
123
{
124
125
126
	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;
	
127
128
	FILE *fin;
	openGridFile(fin, filename);
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
	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);
}
145

146
147
148
149
150
151
152
153
154
155
156
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 );
157
	if (param==IOUtils::npos)
158
159
160
161
162
163
164
165
166
167
168
		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);
	}
169
	fclose(fin);
170
171
}

172
173
174
175
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(), ':', '.');
176
177
178
	const std::string filename = grid2dpath_in + "/" + date_str + ext;
	FILE *fin;
	openGridFile(fin, filename);
179
180
181
182
183
184
185
	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);
	}
186

187
188
189
190
191
192
193
194
195
	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
	
196
	//Radiation parameters
197
	if (parameter==MeteoGrids::ISWR) readGridLayer(fin, filename, "radsw", 2, grid_out);
198
	if (parameter==MeteoGrids::RSWR) {
199
		Grid2DObject net;
200
201
		readGridLayer(fin, filename, "radsw", 2, grid_out);
		readGridLayer(fin, filename, "radswnet", 2, net);
202
203
		grid_out.grid2D -= net.grid2D;
	}
204
	if (parameter==MeteoGrids::ILWR) readGridLayer(fin, filename, "radlwin", 2, grid_out);
205
	if (parameter==MeteoGrids::ALB) {
206
		Grid2DObject rswr, iswr;
207
208
		readGridLayer(fin, filename, "radsw", 2, iswr);
		readGridLayer(fin, filename, "radswnet", 2, grid_out); //net radiation
209
210
211
212
213
		rswr.grid2D = iswr.grid2D - grid_out.grid2D;
		grid_out.grid2D = iswr.grid2D/rswr.grid2D;
	}

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

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

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

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

274
	rewind(fin);
275
276
}

277
void ARPSIO::read3DGrid(Grid3DObject& grid_out, const std::string& i_name)
278
{
279
	const size_t pos = i_name.find_last_of(":");//a specific parameter can be provided as {filename}:{parameter}
280
	const std::string filename = (pos!=IOUtils::npos)? grid3dpath_in +"/" + i_name.substr(0, pos) : grid3dpath_in +"/" + i_name;
281
282
283
284
285
286
287
288
289
290
	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);
	}
	
	std::string parameter = (pos!=IOUtils::npos)?  i_name.substr(pos+1) : "";
	if (parameter=="DEM") { //this is called damage control... this is so ugly...
		parameter = (is_true_arps)? "zp coordinat" : "zp_coordinat";
	}
	
291
292
	FILE *fin;
	openGridFile(fin, filename);
293
294

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

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

318
	fclose(fin);
319
320
321
322
}

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

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

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

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

370
371
}

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

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

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

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

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

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

451
452
453
454
455
/** @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
456
 * @param layer     Index of the layer to skip to (1 to dimz)
457
*/
458
void ARPSIO::skipToLayer(FILE* &fin, const std::string& filename, const unsigned int& layers) const
459
460
461
462
463
464
465
466
467
468
469
470
471
{
	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);
			}
		}
	}
}

472
473
474
475
476
477
478
479
480
481
482
/** @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.
*/
483
void ARPSIO::readGridLayer(FILE* &fin, const std::string& filename, const std::string& parameter, const unsigned int& layer, Grid2DObject& grid)
484
{
485
	if (layer<1 || layer>dimz) {
486
		fclose(fin);
487
		ostringstream tmp;
488
		tmp << "Layer " << layer << " does not exist in ARPS file " << filename << " (nr layers=" << dimz << ")";
489
490
491
492
493
494
495
496
497
		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
498
	moveToMarker(fin, filename, parameter);
499
500

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

	//read the data we are interested in
504
505
	for (size_t iy = 0; iy < dimy; iy++) {
		for (size_t ix = 0; ix < dimx; ix++) {
506
			double tmp;
507
508
			const int readPts = fscanf(fin," %16lf%*[\n]",&tmp);
			if (readPts==1) {
509
				grid(ix,iy) = tmp;
510
			} else {
511
				char dummy[ARPS_MAX_LINE_LENGTH];
512
				const int status = fscanf(fin,"%s",dummy);
513
				fclose(fin);
514
515
516
517
				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);
518
			}
519
520
521
522
		}
	}
}

523
void ARPSIO::moveToMarker(FILE* &fin, const std::string& filename, const std::string& marker)
524
525
526
{
	char dummy[ARPS_MAX_LINE_LENGTH];
	do {
527
528
529
530
531
		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);
		}
532
533
534
		if (dummy[0]>='A' && dummy[0]<='z') { //this is a "marker" line
			if (strncmp(dummy,marker.c_str(), marker.size()) == 0) return;
		}
535
	} while (!feof(fin));
536
	
537
538
	fclose(fin);
	throw InvalidFormatException("End of file "+filename+" should NOT have been reached when looking for "+marker, AT);
539
540
541
}

} //namespace