WSL/SLF GitLab Repository

ARPSIO.cc 20.3 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
36
37
38
39
 * 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).
 * 
 * It is possible to manually extract a whole section, for a given variable with the following script, in order to ease debugging:
 * @code
 * awk '/^ [a-z]/{isVar=0} /^ radsw/{isVar=1;next} {if(isVar) print $0}' {arps_file}
 * @endcode
40
41
 *
 * @section arps_units Units
42
 * All units are assumed to be MKSA.
43
44
45
46
47
48
 *
 * @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
49
50
51
52
 * - 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
 * - GRID2DEXT: arps file extension, or <i>none</i> for no file extension (default: .asc)
53
54
55
 */

const double ARPSIO::plugin_nodata = -999.; //plugin specific nodata value
56
const std::string ARPSIO::default_ext=".asc"; //filename extension
57

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

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

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

98
99
void ARPSIO::setOptions()
{
100
101
	const string grid_in = cfg.get("GRID2D", "Input", IOUtils::nothrow);
	if (grid_in == "ARPS") { //keep it synchronized with IOHandler.cc for plugin mapping!!
102
103
104
		cfg.getValue("GRID2DPATH", "Input", grid2dpath_in);
	}

105
106
	cfg.getValue("ARPS_XCOORD", "Input", xcoord, IOUtils::nothrow);
	cfg.getValue("ARPS_YCOORD", "Input", ycoord, IOUtils::nothrow);
107

108
	//default value has been set in constructor
109
	cfg.getValue("GRID2DEXT", "Input", ext, IOUtils::nothrow);
110
	if (ext=="none") ext.clear();
111
112
}

113
void ARPSIO::listFields(const std::string& i_name)
114
{
115
116
117
	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;
	
118
119
	FILE *fin;
	openGridFile(fin, filename);
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
	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);
}
136

137
138
139
140
141
142
143
144
145
146
147
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 );
148
	if (param==IOUtils::npos)
149
150
151
152
153
154
155
156
157
158
159
		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);
	}
160
	fclose(fin);
161
162
}

163
164
165
166
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(), ':', '.');
167
168
169
	const std::string filename = grid2dpath_in + "/" + date_str + ext;
	FILE *fin;
	openGridFile(fin, filename);
170
171
172
173
174
175
176
	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);
	}
177

178
179
180
181
182
183
184
185
186
	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
	
187
	//Radiation parameters
188
	if (parameter==MeteoGrids::ISWR) readGridLayer(fin, filename, "radsw", 1, grid_out);
189
	if (parameter==MeteoGrids::RSWR) {
190
		Grid2DObject net;
191
192
		readGridLayer(fin, filename, "radsw", 1, grid_out);
		readGridLayer(fin, filename, "radswnet", 1, net);
193
194
		grid_out.grid2D -= net.grid2D;
	}
195
	if (parameter==MeteoGrids::ILWR) readGridLayer(fin, filename, "radlwin", 1, grid_out);
196
	if (parameter==MeteoGrids::ALB) {
197
		Grid2DObject rswr, iswr;
198
199
		readGridLayer(fin, filename, "radsw", 1, iswr);
		readGridLayer(fin, filename, "radswnet", 1, grid_out); //net radiation
200
201
202
203
204
		rswr.grid2D = iswr.grid2D - grid_out.grid2D;
		grid_out.grid2D = iswr.grid2D/rswr.grid2D;
	}

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

	//Basic meteo parameters
230
	if (parameter==MeteoGrids::P) readGridLayer(fin, filename, "p", 2, grid_out);
231
	if (parameter==MeteoGrids::TSG) readGridLayer(fin, filename, "tsoil", 1, grid_out); //or is it tss for us?
232
	/*if (parameter==MeteoGrids::RH) {
233
		//const double epsilon = Cst::gaz_constant_dry_air / Cst::gaz_constant_water_vapor;
234
		readGridLayer(filename, "qv", 2, grid_out); //water vapor mixing ratio
235
		//Atmosphere::waterSaturationPressure(T);
236
		//TODO: compute relative humidity out of it!
237
		//through potential temperature "pt" -> local temperature?
238
239
240
	}*/

	//Hydrological parameters
241
	if (parameter==MeteoGrids::HS) readGridLayer(fin, filename, "snowdpth", 1, grid_out);
242
	if (parameter==MeteoGrids::PSUM) {
243
		readGridLayer(fin, filename, "prcrate1", 1, grid_out); //in kg/m^2/s
244
245
246
247
		grid_out.grid2D *= 3600.; //we need kg/m^2/h
	}

	//DEM
248
249
	const std::string dem_marker = (is_true_arps)? "zp coordinat" : "zp_coordinat";
	if (parameter==MeteoGrids::DEM) readGridLayer(fin, filename, dem_marker, 1, grid_out);
250
	if (parameter==MeteoGrids::SLOPE) {
251
252
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
253
		readGridLayer(fin, filename, dem_marker, 1, dem);
254
		dem.update();
255
		grid_out.set(dem.cellsize, dem.llcorner, dem.slope);
256
	}
257
	if (parameter==MeteoGrids::AZI) {
258
259
		DEMObject dem;
		dem.setUpdatePpt(DEMObject::SLOPE);
260
		readGridLayer(fin, filename, dem_marker, 1, dem);
261
		dem.update();
262
		grid_out.set(dem.cellsize, dem.llcorner, dem.azi);
263
264
	}

265
	rewind(fin);
266
267
}

268
void ARPSIO::read3DGrid(Grid3DObject& grid_out, const std::string& i_name)
269
{
270
271
272
273
274
275
276
277
278
279
280
281
	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);
	}
	
	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";
	}
	
282
283
	FILE *fin;
	openGridFile(fin, filename);
284
285

	//resize the grid just in case
286
287
288
289
	Coords llcorner(coordin, coordinparam);
	llcorner.setXY(xcoord, ycoord, IOUtils::nodata);
	grid_out.set(dimx, dimy, dimz, cellsize, llcorner);
	grid_out.z = zcoord;
290
291

	// Read until the parameter is found
292
	moveToMarker(fin, filename, parameter);
293
294

	//read the data we are interested in
295
296
297
	for (size_t ix = 0; ix < dimx; ix++) {
		for (size_t iy = 0; iy < dimy; iy++) {
			for (size_t iz = 0; iz < dimz; iz++) {
298
				double tmp;
299
				if (fscanf(fin," %16lf%*[\n]",&tmp)==1) {
300
301
					grid_out.grid3D(ix,iy,iz) = tmp;
				} else {
302
303
					fclose(fin);
					throw InvalidFormatException("Failure in reading 3D grid in file "+filename, AT);
304
				}
305
306
307
308
			}
		}
	}

309
	fclose(fin);
310
311
312
313
}

void ARPSIO::readDEM(DEMObject& dem_out)
{
314
315
316
	const std::string filename = cfg.get("DEMFILE", "Input");
	FILE *fin;
	openGridFile(fin, filename);
317
	read2DGrid_internal(fin, filename, dem_out, MeteoGrids::DEM);
318
	fclose(fin);
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
}

void ARPSIO::readLanduse(Grid2DObject& /*landuse_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::readAssimilationData(const Date& /*date_in*/, Grid2DObject& /*da_out*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::readStationData(const Date&, std::vector<StationData>& /*vecStation*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::readMeteoData(const Date& /*dateStart*/, const Date& /*dateEnd*/,
Mathias Bavay's avatar
Mathias Bavay committed
340
                           std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
341
                           const size_t&)
342
343
344
345
346
347
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::writeMeteoData(const std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
Mathias Bavay's avatar
Mathias Bavay committed
348
                            const std::string&)
349
350
351
352
353
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

354
void ARPSIO::readPOI(std::vector<Coords>&)
355
356
357
358
359
360
361
362
363
364
365
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

void ARPSIO::write2DGrid(const Grid2DObject& /*grid_in*/, const std::string& /*name*/)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

366
367
368
369
370
371
void ARPSIO::write2DGrid(const Grid2DObject&, const MeteoGrids::Parameters&, const Date&)
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

372
void ARPSIO::initializeGRIDARPS(FILE* &fin, const std::string& filename)
373
374
375
376
{
	double v1, v2;

	//go to read the sizes
377
	moveToMarker(fin, filename, "nnx");
378
	//finish reading the line and move to the next one
379
	if (fscanf(fin,"%*[^\n]")!=0) {
380
		fclose(fin);
381
382
		throw InvalidFormatException("Error in file format of file "+filename, AT);
	}
383
	if (fscanf(fin," %u %u %u \n",&dimx,&dimy,&dimz)!=3) {
384
		fclose(fin);
385
386
387
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
388
		fclose(fin);
389
390
391
392
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

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

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

413
414
}

415
void ARPSIO::initializeTrueARPS(FILE* &fin, const std::string& filename, const char curr_line[ARPS_MAX_LINE_LENGTH])
416
417
418
419
420
{
	double v1, v2;

	//go to read the sizes
	if (sscanf(curr_line," nx = %u, ny = %u, nz = %u ",&dimx,&dimy,&dimz)!=3) {
421
		fclose(fin);
422
423
424
		throw InvalidFormatException("Can not read dimx, dimy, dimz from file "+filename, AT);
	}
	if (dimx==0 || dimy==0 || dimz==0) {
425
		fclose(fin);
426
427
428
429
		throw IndexOutOfBoundsException("Invalid dimx, dimy, dimz from file "+filename, AT);
	}

	//initializing cell size
430
	moveToMarker(fin, filename, "x coordinate");
431
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
432
		fclose(fin);
433
434
435
		throw InvalidFormatException("Can not read first two x coordinates from file "+filename, AT);
	}
	const double cellsize_x = v2 - v1;
436
	moveToMarker(fin, filename, "y coordinate");
437
	if (fscanf(fin,"%lg %lg",&v1,&v2)!=2) {
438
		fclose(fin);
439
440
441
		throw InvalidFormatException("Can not read first two y coordinates from file "+filename, AT);
	}
	const double cellsize_y = v2 - v1;
442
	if (cellsize_x!=cellsize_y) {
443
		fclose(fin);
444
445
446
		throw InvalidFormatException("Only square cells currently supported! Non compliance in file "+filename, AT);
	}
	cellsize = cellsize_y;
447

448
	moveToMarker(fin, filename, "z coordinate");
449
450
451
	while (fscanf(fin,"%lg",&v1)==1) {
		zcoord.push_back( v1 );
	}
452
	if (zcoord.size()!=dimz) {
453
		ostringstream ss;
454
		ss << "Expected " << dimz << " z coordinates in file \""+filename+"\", found " << zcoord.size();
455
		fclose(fin);
456
457
		throw InvalidFormatException(ss.str(), AT);
	}
458
459
}

460
void ARPSIO::openGridFile(FILE* &fin, const std::string& filename)
461
462
463
{
	unsigned int v1;

464
	if (!IOUtils::fileExists(filename)) throw AccessException(filename, AT); //prevent invalid filenames
465
	if ((fin=fopen(filename.c_str(),"r")) == NULL) {
466
		fclose(fin);
467
		throw AccessException("Can not open file "+filename, AT);
468
469
470
471
	}

	//identify if the file is an original arps file or a file modified by ARPSGRID
	char dummy[ARPS_MAX_LINE_LENGTH];
472
	for (unsigned char j=0; j<5; j++) {
473
		//the first easy difference in the structure happens at line 5
474
		if (fgets(dummy,ARPS_MAX_STRING_LENGTH,fin)==NULL) {
475
			fclose(fin);
476
477
			throw InvalidFormatException("Fail to read header lines of file "+filename, AT);
		}
478
	}
479
	zcoord.clear(); //reset the zcoord
480
481
482
	if (sscanf(dummy," nx = %u, ny = ", &v1)<1) {
		//this is an ASCII file modified by ARPSGRID
		is_true_arps=false;
483
		initializeGRIDARPS(fin, filename);
484
485
	} else {
		//this is a true ARPS file
486
		initializeTrueARPS(fin, filename, dummy);
487
	}
488
489
490
	//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;
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506

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

/** @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.
*/
507
void ARPSIO::readGridLayer(FILE* &fin, const std::string& filename, const std::string& parameter, const unsigned int& layer, Grid2DObject& grid)
508
{
509
	if (layer<1 || layer>dimz) {
510
		fclose(fin);
511
		ostringstream tmp;
512
		tmp << "Layer " << layer << " does not exist in ARPS file " << filename << " (nr layers=" << dimz << ")";
513
514
515
516
517
518
519
520
521
		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
522
	moveToMarker(fin, filename, parameter);
523
524

	// move to the begining of the layer of interest
525
	if (layer>1) {
526
		double tmp;
527
		const size_t jmax=dimx*dimy*(layer-1);
528
		for (size_t j = 0; j < jmax; j++) {
529
			if (fscanf(fin," %16lf%*[\n]",&tmp)==EOF) {
530
				fclose(fin);
531
532
				throw InvalidFormatException("Fail to skip data layers in file "+filename, AT);
			}
533
		}
534
535
536
	}

	//read the data we are interested in
537
538
	for (size_t iy = 0; iy < dimy; iy++) {
		for (size_t ix = 0; ix < dimx; ix++) {
539
			double tmp;
540
541
			const int readPts = fscanf(fin," %16lf%*[\n]",&tmp);
			if (readPts==1) {
542
				grid(ix,iy) = tmp;
543
			} else {
544
				char dummy[ARPS_MAX_LINE_LENGTH];
545
				(void)fscanf(fin,"%s",dummy);
546
				fclose(fin);
547
				throw InvalidFormatException("Fail to read data layer for parameter '"+parameter+"' in file '"+filename+"', instead read: '"+string(dummy)+"'", AT);
548
			}
549
550
551
552
		}
	}
}

553
void ARPSIO::moveToMarker(FILE* &fin, const std::string& filename, const std::string& marker)
554
555
{
	char dummy[ARPS_MAX_LINE_LENGTH];
556
	int nb_elems = 0;
557
	do {
558
559
560
561
562
		nb_elems = fscanf(fin," %[^\t\n] ",dummy); //HACK: possible buffer overflow
		if (dummy[0]>='A' && dummy[0]<='z') { //this is a "marker" line
			if (strncmp(dummy,marker.c_str(), marker.size()) == 0) return;
		}
	} while (!feof(fin) && nb_elems!=0);
563
	
564
	if (feof(fin)) {
565
		fclose(fin);
566
567
568
		const std::string message = "End of file "+filename+" should NOT have been reached when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
569
	if (nb_elems==0) { //there was a line that did not match the format
570
		fclose(fin);
571
572
573
		const std::string message = "Matching failure in file "+filename+" when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
574
575
576
}

} //namespace