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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
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 );
	if (param==IOUtils::npos) {
		throw InvalidArgumentException("Invalid MeteoGrids Parameter requested: '"+param_str+"'", AT);
		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);
	}
162
	fclose(fin);
163
164
}

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

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

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

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

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

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

267
	rewind(fin);
268
269
}

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

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

	// Read until the parameter is found
294
	moveToMarker(fin, filename, parameter);
295
296

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

311
	fclose(fin);
312
313
314
315
}

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

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
342
                           std::vector< std::vector<MeteoData> >& /*vecMeteo*/,
343
                           const size_t&)
344
345
346
347
348
349
{
	//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
350
                            const std::string&)
351
352
353
354
355
{
	//Nothing so far
	throw IOException("Nothing implemented here", AT);
}

356
void ARPSIO::readPOI(std::vector<Coords>&)
357
358
359
360
361
362
363
364
365
366
367
{
	//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);
}

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

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

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

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

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

415
416
}

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

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

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

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

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

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

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

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

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

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

555
void ARPSIO::moveToMarker(FILE* &fin, const std::string& filename, const std::string& marker)
556
557
{
	char dummy[ARPS_MAX_LINE_LENGTH];
558
	int nb_elems = 0;
559
	do {
560
561
562
563
564
		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);
565
	
566
	if (feof(fin)) {
567
		fclose(fin);
568
569
570
		const std::string message = "End of file "+filename+" should NOT have been reached when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
571
	if (nb_elems==0) { //there was a line that did not match the format
572
		fclose(fin);
573
574
575
		const std::string message = "Matching failure in file "+filename+" when looking for "+marker;
		throw InvalidFormatException(message, AT);
	}
576
577
578
}

} //namespace