#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <math.h>
#include <float.h>
#include "gribwlib.h"
#include "pds4.h"
#define HOUR 1
#define DAY 2
#define HOURS3 10
#define HOURS6 11
#define HOURS12 12
/* #define IN_ORDER */
/*
* take a 6-hourly time series, and make a monthly mean
*
* v0.1 10/99
* v0.2 4/00 handles averages of averages
* v0.3 9/00 fixed time_range = 113 -> time_range == 113
* v0.4 12/01 added time_range = 4 -> time_range == 113A
* v0.5 1/02 minor stuff .. overwriting data not found?
* v0.6 12/03 fixup remove monthly mean hard coding
*/
#define N 100000
void update_data(float *sum, int *count,float *data, int ndata);
int write_mean(float *sum, int *count, unsigned char *pds,
unsigned char *gds, int ndata, int dhr,
int NumAve, int NumMissing, int scale2, FILE *output);
int main(int argc, char **argv) {
long int len_grib, pos = 0;
unsigned char *pds, *gds, *old_pds, *old_gds;
FILE *input, *output;
float *sum;
int *count;
int array_size;
unsigned int field, old_field;
int time;
int new_mean, hour, dhr, eof, i;
int scale10, scale2, old_scale2;
char record[400];
float *data;
int ndata, old_ndata, NumAve, NumMissing;
/* preliminaries .. open up all files */
if (argc != 4) {
fprintf(stderr, "%s [in gribfile] [out gribfile] dhr\n", argv[0]);
exit(8);
}
if ((input = fopen(argv[1],"rb")) == NULL) {
fprintf(stderr,"could not open file: %s\n", argv[1]);
exit(7);
}
if ((output = fopen(argv[2],"wb+")) == NULL) {
fprintf(stderr,"could not open file: %s\n", argv[2]);
exit(7);
}
dhr = atoi(argv[3]);
if (dhr > 24) {
fprintf(stderr,"dhr must be less than or equal to 24 (hours)\n");
exit(8);
}
if (24 % dhr != 0) {
fprintf(stderr,"dhr hours must be multiple of 24\n");
exit(8);
}
count = (int *) malloc(N * sizeof(int));
sum = (float *) malloc(N * sizeof(float));
if (count == NULL || sum == NULL) {
fprintf(stderr,"memeory allocation failed\n");
exit(8);
}
array_size = N;
for (i = 0; i < array_size; i++) {
count[i] = 0;
sum[i] = 0.0;
}
pos = 0;
new_mean = 1;
old_field = -1;
NumAve = NumMissing = 0;
old_ndata = 0;
eof = 0;
old_pds = old_gds = NULL;
for(;;) {
#ifdef IN_ORDER
len_grib = rd_grib_rec2(input, pos, &pds, &gds, &data, &ndata,
&scale10, &scale2);
if (len_grib <= 0) eof = 1;
#else
if (fgets(record, sizeof(record), stdin) == NULL) eof = 1;
if (sscanf(record, "%d:%ld:", &i, &pos) != 2) eof = 1;
if (eof == 0) {
len_grib = rd_grib_rec2(input, pos, &pds, &gds, &data, &ndata,
&scale10, &scale2);
if (len_grib <= 0) eof = 1;
}
#endif
/* check to see if we have a new variable/month */
if (eof == 1) new_mean = 1;
field = PDS_Field(pds);
if (old_field != field) new_mean = 1;
time = get_InitYYYYMMDDHH(pds);
if (new_mean == 1 && NumAve > 0) {
/* printf("field %d nave=%d\n",field, NumAve); */
write_mean(sum,count,old_pds,old_gds,
old_ndata,dhr,NumAve,NumMissing,old_scale2,output);
if (old_pds == NULL) printf("old_pds == NULL\n");
if (old_gds == NULL) printf("old_gds == NULL\n");
if (old_pds) free(old_pds);
if (old_gds) free(old_gds);
old_pds = old_gds = NULL;
}
if (eof) break;
if (new_mean) {
/* initialize for new mean */
if (ndata > array_size) {
count = (int *) realloc(count, ndata * sizeof(int));
sum = (float *) realloc(sum, ndata * sizeof(float));
array_size = ndata;
if (count == NULL || sum == NULL) {
fprintf(stderr,"memeory allocation failed\n");
exit(8);
}
}
for (i = 0; i < ndata; i++) {
count[i] = 0;
sum[i] = 0.0;
}
old_field = field;
NumAve = NumMissing = new_mean = 0;
old_ndata = ndata;
old_pds = cpGRIBsec(pds);
old_gds = cpGRIBsec(gds);
old_scale2 = scale2;
}
if (ndata != old_ndata) {
fprintf(stderr,"fatal error .. size of data array changed!"
" %d %d\n",ndata,old_ndata);
exit(8);
}
hour = PDS_Hour(pds);
update_data(sum,count,data,ndata);
NumAve += (PDS_NumAve(pds) > 1 ? PDS_NumAve(pds) : 1);
NumMissing += PDS_NumMissing(pds);
pos += len_grib;
}
fclose(input);
fclose(output);
exit (0);
}
void update_data(float *sum, int *count, float *data, int ndata)
{
int i;
for (i = 0; i < ndata; i++) {
if (! UNDEFINED_VAL(data[i])) {
sum[i] += data[i];
count[i]++;
}
}
}
static int days_in_month(int yyyymm) {
static int days[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
int year, mo = yyyymm % 100;
if (mo <= 0 || mo > 12) exit(9);
if (mo != 2) return days[mo-1];
year = yyyymm / 100;
if (year % 4 != 0) return 28;
if (year % 100 != 0) return 29;
if (year % 400 != 0) return 28;
return 29;
}
int write_mean(float *sum, int *count, unsigned char *pds,
unsigned char *gds, int ndata, int dhr,
int NumAve, int NumMissing, int scale2, FILE *output) {
int dec_scale, i, count_max, p1, p2, time_range, time_units, new_time_range, nhour;
static int m_count = 0;
if (NumAve == 0) return 0;
time_range = PDS_TimeRange(pds);
for (count_max = i = 0; i < ndata; i++) {
count_max = count_max < count[i] ? count[i] : count_max;
if (count[i] > 0) {
sum[i] /= count[i];
}
else {
sum[i] = UNDEFINED;
}
}
printf("%d: writing ave out %d points time_range %d nave=%d\n", ++m_count, ndata, time_range,count_max);
/* number of available and missing obs */
pds = PDStool(pds,P_used(NumAve),P_missing(NumMissing),P_end);
/* increase precision by 1 digit */
/* if (NumAve > 10) dec_scale = PDS_DecimalScale(pds) + 1; */
/* pds = PDStool(pds, P_dec_scale(dec_scale), P_end); */
time_range = PDS_TimeRange(pds);
p1 = PDS_P1(pds);
p2 = PDS_P2(pds);
time_units = PDS_ForecastTimeUnit(pds);
if (time_units != HOUR) {
fprintf(stderr,"help .. Forecast Time Unit not hours\n");
fprintf(stderr,"extend code time_unit = %d\n",time_units);
exit(8);
}
if (time_range == 0) {
new_time_range = 113;
p2 = dhr;
}
else if (time_range == 1) {
new_time_range = 113;
p1 = 0;
p2 = dhr;
}
else if (time_range == 3 && dhr == 24) {
new_time_range = 130;
}
else if (time_range == 3 && dhr == p2-p1) {
new_time_range = 131;
}
else if (time_range == 4 && dhr == 24) {
new_time_range = 128;
}
else if (time_range == 4 && dhr == p2-p1) {
new_time_range = 129;
}
else if (time_range == 10) {
new_time_range = 113;
if (p1 != 0) {
fprintf(stderr,"forecast hour too big 255 used\n");
p1 = 255;
}
else {
p1 = p2;
}
p2 = dhr;
}
else if (time_range == 113) {
new_time_range = 113;
p2 = dhr < p2 ? dhr : p2;
}
else if (time_range == 128 && dhr == p2 - p1) {
new_time_range = 129;
}
else if (time_range == 130 && dhr == p2 - p1) {
new_time_range = 131;
}
else {
fprintf(stderr,"help .. mean with time range %d p1 %d p2 %d\n",time_range,p1,p2);
fprintf(stderr,"record not written\n");
return 0;
}
pds = PDStool(pds, P_time_range(new_time_range),P_p1(p1),P_p2(p2),P_end);
/* eta style output - works with mrf */
set_BDSMinBits(-1);
set_def_power2(-scale2);
/* dont take averages of certain fields */
if (time_range != 113 || (PDS_KPDS6(pds) != 2 && PDS_KPDS6(pds) != 3)) {
wrt_grib_rec(pds, gds, sum, ndata, output);
}
return 0;
}