Skip to content

Commit

Permalink
Added binary I/O for atmospheric data.
Browse files Browse the repository at this point in the history
  • Loading branch information
@lars2015
lars2015 committed Jun 29, 2018
1 parent cae7a3b commit d1b357a
Show file tree
Hide file tree
Showing 4 changed files with 130 additions and 59 deletions.
Binary file modified doc/refman.pdf
View file
Binary file not shown.
130 changes: 93 additions & 37 deletions src/libtrac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copright (C) 2013-2015 Forschungszentrum Juelich GmbH
Copright (C) 2013-2018 Forschungszentrum Juelich GmbH
*/

/*!
Expand Down Expand Up @@ -426,23 +426,50 @@ void read_atm(
if (!(in = fopen(filename, "r")))
ERRMSG("Cannot open file!");

/* Read line... */
while (fgets(line, LEN, in)) {

/* Read data... */
TOK(line, tok, "%lg", atm->time[atm->np]);
TOK(NULL, tok, "%lg", atm->p[atm->np]);
TOK(NULL, tok, "%lg", atm->lon[atm->np]);
TOK(NULL, tok, "%lg", atm->lat[atm->np]);
/* Read binary data... */
if (ctl->atm_bin && ctl->atm_gpfile[0] == '-') {
FREAD(&atm->np, int,
1,
in);
FREAD(atm->time, double,
(size_t) atm->np,
in);
FREAD(atm->p, double,
(size_t) atm->np,
in);
FREAD(atm->lon, double,
(size_t) atm->np,
in);
FREAD(atm->lat, double,
(size_t) atm->np,
in);
for (iq = 0; iq < ctl->nq; iq++)
TOK(NULL, tok, "%lg", atm->q[iq][atm->np]);
FREAD(atm->q[iq], double,
(size_t) atm->np,
in);
}

/* Read ASCII data... */
else {

/* Convert altitude to pressure... */
atm->p[atm->np] = P(atm->p[atm->np]);
/* Read line... */
while (fgets(line, LEN, in)) {

/* Increment data point counter... */
if ((++atm->np) > NP)
ERRMSG("Too many data points!");
/* Read data... */
TOK(line, tok, "%lg", atm->time[atm->np]);
TOK(NULL, tok, "%lg", atm->p[atm->np]);
TOK(NULL, tok, "%lg", atm->lon[atm->np]);
TOK(NULL, tok, "%lg", atm->lat[atm->np]);
for (iq = 0; iq < ctl->nq; iq++)
TOK(NULL, tok, "%lg", atm->q[iq][atm->np]);

/* Convert altitude to pressure... */
atm->p[atm->np] = P(atm->p[atm->np]);

/* Increment data point counter... */
if ((++atm->np) > NP)
ERRMSG("Too many data points!");
}
}

/* Close file... */
Expand Down Expand Up @@ -618,6 +645,8 @@ void read_ctl(
scan_ctl(filename, argc, argv, "ATM_DT_OUT", -1, "86400", NULL);
ctl->atm_filter =
(int) scan_ctl(filename, argc, argv, "ATM_FILTER", -1, "0", NULL);
ctl->atm_bin =
(int) scan_ctl(filename, argc, argv, "ATM_BIN", -1, "0", NULL);

/* Output of CSI data... */
scan_ctl(filename, argc, argv, "CSI_BASENAME", -1, "-", ctl->csi_basename);
Expand Down Expand Up @@ -1118,7 +1147,7 @@ void timer(

/* Print timer... */
else if (mode == 3)
printf("%s = %g s\n", name, runtime[id]);
printf("%s = %.3f s\n", name, runtime[id]);
}

/*****************************************************************************/
Expand Down Expand Up @@ -1320,31 +1349,58 @@ void write_atm(
ERRMSG("Cannot create file!");
}

/* Write header... */
fprintf(out,
"# $1 = time [s]\n"
"# $2 = altitude [km]\n"
"# $3 = longitude [deg]\n" "# $4 = latitude [deg]\n");
for (iq = 0; iq < ctl->nq; iq++)
fprintf(out, "# $%i = %s [%s]\n", iq + 5, ctl->qnt_name[iq],
ctl->qnt_unit[iq]);
fprintf(out, "\n");
/* Write binary data... */
if (ctl->atm_bin && ctl->atm_gpfile[0] == '-') {
FWRITE(&atm->np, int,
1,
out);
FWRITE(atm->time, double,
(size_t) atm->np,
out);
FWRITE(atm->p, double,
(size_t) atm->np,
out);
FWRITE(atm->lon, double,
(size_t) atm->np,
out);
FWRITE(atm->lat, double,
(size_t) atm->np,
out);
for (iq = 0; iq < ctl->nq; iq++)
FWRITE(atm->q[iq], double,
(size_t) atm->np,
out);
}

/* Write data... */
for (ip = 0; ip < atm->np; ip++) {
/* Write ASCII data... */
else {

/* Check time... */
if (ctl->atm_filter && (atm->time[ip] < t0 || atm->time[ip] > t1))
continue;
/* Write header... */
fprintf(out,
"# $1 = time [s]\n"
"# $2 = altitude [km]\n"
"# $3 = longitude [deg]\n" "# $4 = latitude [deg]\n");
for (iq = 0; iq < ctl->nq; iq++)
fprintf(out, "# $%i = %s [%s]\n", iq + 5, ctl->qnt_name[iq],
ctl->qnt_unit[iq]);
fprintf(out, "\n");

/* Write output... */
fprintf(out, "%.2f %g %g %g", atm->time[ip], Z(atm->p[ip]),
atm->lon[ip], atm->lat[ip]);
for (iq = 0; iq < ctl->nq; iq++) {
fprintf(out, " ");
fprintf(out, ctl->qnt_format[iq], atm->q[iq][ip]);
/* Write data... */
for (ip = 0; ip < atm->np; ip++) {

/* Check time... */
if (ctl->atm_filter && (atm->time[ip] < t0 || atm->time[ip] > t1))
continue;

/* Write output... */
fprintf(out, "%.2f %g %g %g", atm->time[ip], Z(atm->p[ip]),
atm->lon[ip], atm->lat[ip]);
for (iq = 0; iq < ctl->nq; iq++) {
fprintf(out, " ");
fprintf(out, ctl->qnt_format[iq], atm->q[iq][ip]);
}
fprintf(out, "\n");
}
fprintf(out, "\n");
}

/* Close file... */
Expand Down
25 changes: 20 additions & 5 deletions src/libtrac.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
You should have received a copy of the GNU General Public License
along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
Copright (C) 2013-2015 Forschungszentrum Juelich GmbH
Copright (C) 2013-2018 Forschungszentrum Juelich GmbH
*/

/*!
Expand Down Expand Up @@ -77,7 +77,7 @@
#define NQ 12

/*! Maximum number of pressure levels for meteorological data. */
#define EP 137
#define EP 111

/*! Maximum number of longitudes for meteorological data. */
#define EX 1201
Expand Down Expand Up @@ -126,6 +126,18 @@
exit(EXIT_FAILURE); \
}

/*! Read binary data. */
#define FREAD(ptr, type, size, out) { \
if(fread(ptr, sizeof(type), size, out)!=size) \
ERRMSG("Error while reading!"); \
}

/*! Write binary data. */
#define FWRITE(ptr, type, size, out) { \
if(fwrite(ptr, sizeof(type), size, out)!=size) \
ERRMSG("Error while writing!"); \
}

/*! Compute linear interpolation. */
#define LIN(x0, y0, x1, y1, x) \
((y0)+((y1)-(y0))/((x1)-(x0))*((x)-(x0)))
Expand Down Expand Up @@ -370,6 +382,9 @@ typedef struct {
/*! Time filter for atmospheric data output (0=no, 1=yes). */
int atm_filter;

/*! Binary I/O of atmospheric data (0=no, 1=yes). */
int atm_bin;

/*! Basename of CSI data files. */
char csi_basename[LEN];

Expand Down Expand Up @@ -523,13 +538,13 @@ typedef struct {
double q[NQ][NP];

/*! Zonal wind perturbation [m/s]. */
double up[NP];
float up[NP];

/*! Meridional wind perturbation [m/s]. */
double vp[NP];
float vp[NP];

/*! Vertical velocity perturbation [hPa/s]. */
double wp[NP];
float wp[NP];

} atm_t;

Expand Down
34 changes: 17 additions & 17 deletions src/trac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,14 +333,14 @@ int main(
PRINT_TIMER(TIMER_SEDI);

/* Report memory usage... */
printf("MEMORY_ATM = %g MByte\n", 2. * sizeof(atm_t) / 1024. / 1024.);
printf("MEMORY_ATM = %g MByte\n", sizeof(atm_t) / 1024. / 1024.);
printf("MEMORY_METEO = %g MByte\n", 2. * sizeof(met_t) / 1024. / 1024.);
printf("MEMORY_DYNAMIC = %g MByte\n",
NP * sizeof(double) / 1024. / 1024.);
4 * NP * sizeof(double) / 1024. / 1024.);
printf("MEMORY_STATIC = %g MByte\n",
(((EX + EY) + (2 + NQ) * GX * GY * GZ) * sizeof(double)
((3 * GX * GY + 4 * GX * GY * GZ) * sizeof(double)
+ (EX * EY + EX * EY * EP) * sizeof(float)
+ (2 * GX * GY * GZ) * sizeof(int)) / 1024. / 1024.);
+ (GX * GY + GX * GY * GZ) * sizeof(int)) / 1024. / 1024.);

/* Report problem size... */
printf("SIZE_NP = %d\n", atm->np);
Expand Down Expand Up @@ -559,15 +559,15 @@ void module_diffusion_meso(
rs = sqrt(1 - r * r);

/* Calculate mesoscale wind fluctuations... */
atm->up[ip] =
r * atm->up[ip] + rs * gsl_ran_gaussian_ziggurat(rng,
ctl->turb_meso * usig);
atm->vp[ip] =
r * atm->vp[ip] + rs * gsl_ran_gaussian_ziggurat(rng,
ctl->turb_meso * vsig);
atm->wp[ip] =
r * atm->wp[ip] + rs * gsl_ran_gaussian_ziggurat(rng,
ctl->turb_meso * wsig);
atm->up[ip] = (float)
(r * atm->up[ip]
+ rs * gsl_ran_gaussian_ziggurat(rng, ctl->turb_meso * usig));
atm->vp[ip] = (float)
(r * atm->vp[ip]
+ rs * gsl_ran_gaussian_ziggurat(rng, ctl->turb_meso * vsig));
atm->wp[ip] = (float)
(r * atm->wp[ip]
+ rs * gsl_ran_gaussian_ziggurat(rng, ctl->turb_meso * wsig));

/* Calculate air parcel displacement... */
atm->lon[ip] += dx2deg(atm->up[ip] * dt / 1000., atm->lat[ip]);
Expand Down Expand Up @@ -632,7 +632,7 @@ void module_isosurf(

static double *iso, *ps, t, *ts;

static int idx, ip2, n, nb = 100000;
static int idx, ip2, n;

FILE *in;

Expand All @@ -649,9 +649,9 @@ void module_isosurf(
ALLOC(iso, double,
NP);
ALLOC(ps, double,
nb);
NP);
ALLOC(ts, double,
nb);
NP);

/* Save pressure... */
if (ctl->isosurf == 1)
Expand Down Expand Up @@ -689,7 +689,7 @@ void module_isosurf(
/* Read pressure time series... */
while (fgets(line, LEN, in))
if (sscanf(line, "%lg %lg", &ts[n], &ps[n]) == 2)
if ((++n) > 100000)
if ((++n) > NP)
ERRMSG("Too many data points!");

/* Check number of points... */
Expand Down

0 comments on commit d1b357a

Please sign in to comment.