120 #define VGROUPCLASS "PlanetaryGrid"
122 #define REGISTRATION CENTER
125 #define MAX_NORTH 89.5
126 #define MAX_SOUTH -89.5
127 #define MAX_WEST -179.375
128 #define MAX_EAST 179.375
129 #define TOMSLATSZ 180
130 #define TOMSLONSZ 288
133 int main(
int argc,
char *argv[]) {
163 float32 stdlimit1, stdlimit2, stdlimit3;
170 float flocnt = 0., fhicnt = 0., coslat;
172 int min_lin = (-20), max_lin = (
TOMSLATSZ + 20), pts_used;
173 float min_lat = -90., max_lat = 90., mmoffset = 10., zen,
lat;
177 int lo_cnts[5] = {0, 0, 0, 0, 0}, hi_cnts[5] = {0, 0, 0, 0, 0};
178 int z_st[4] = {-1, -1, -1, -1}, z_en[4] = {-1, -1, -1, -1};
179 int lo_vals[5] = {50, 100, 150, 200, 250},
180 hi_vals[5] = {400, 450, 500, 550, 600};
181 int in_z_run = 0, zcur = 0, zline, ipx, iv;
188 int32 gridid, sdsid, geomid;
196 int16 *int_SDSdataRT;
197 int16 *int_SDSdataCLMAVG;
198 int16 *int_SDSdataCLMSTD;
199 float32 *float_SDSdataQC;
213 strcpy(infileRT, argv[1]);
214 strcpy(infileCLM, argv[2]);
216 strcpy(monthstr, argv[4]);
218 if (!strcmp(
lowcase(argv[5]),
"ozone"))
219 strcpy(vgroupname,
"Geophysical Data");
220 else if (!strcmp(
lowcase(argv[5]),
"uwind"))
221 pexit(
"Only OZONE parameter supported by this program.");
222 else if (!strcmp(
lowcase(argv[5]),
"vwind"))
223 pexit(
"Only OZONE parameter supported by this program.");
224 else if (!strcmp(
lowcase(argv[5]),
"pres"))
225 pexit(
"Only OZONE parameter supported by this program.");
226 else if (!strcmp(
lowcase(argv[5]),
"rhum"))
227 pexit(
"Only OZONE parameter supported by this program.");
228 else pexit(
"Invalid Parameter string specified");
230 sscanf(argv[6],
"%f", &stdlimit1);
231 sscanf(argv[7],
"%f", &stdlimit2);
232 sscanf(argv[8],
"%f", &stdlimit3);
233 sscanf(argv[9],
"%d", &maxmissing);
234 if (argc >= 11) sscanf(argv[10],
"%d", &loval);
235 if (argc >= 12) sscanf(argv[11],
"%d", &hival);
236 if (argc >= 13) sscanf(argv[12],
"%d", &outmax);
237 if (argc >= 14) sscanf(argv[13],
"%f", &
min_lat);
238 if (argc >= 15) sscanf(argv[14],
"%f", &max_lat);
241 strncpy(gridflag, argv[15], 1);
244 if (argc == 17) ifileflag = atoi(argv[16]);
253 array_size = shape[0] * shape[1];
258 pexit(
"calloc int_SDSdataRT");
260 if ((int_SDSdataCLMAVG =
263 pexit(
"calloc int_SDSdataCLMAVG");
265 if ((int_SDSdataCLMSTD =
268 pexit(
"calloc int_SDSdataCLMSTD");
270 if ((float_SDSdataQC =
272 (float32 *) calloc(array_size,
sizeof (float32))) ==
NULL)
273 pexit(
"calloc float_SDSdataQC");
294 strcpy(sdsname,
"ozone_mean");
296 strcpy(vgname, vgroupname);
304 pexit(
"in getting day, month");
309 if (argc >= 15) mmoffset = max_lat;
310 zen = -23. * cos((
jday + 10.) * 2 *
pi / 365.);
311 min_lat = -90. + zen + mmoffset;
313 max_lat = 90. + zen - mmoffset;
314 if (max_lat > 90.) max_lat = 90.;
315 printf(
"day dependent min, max latitude chosen:\n");
316 printf(
"day = %d, mmoffset = %f, sun lat = %f\n",
317 jday, mmoffset, zen);
321 if ((
rdsds(infileRT, vgname, sdsname, inShape,
322 int_SDSdataRT)) != 0)
pexit(
"rdsds RT");
324 if ((inShape[0] != shape[0]) || (inShape[1] != shape[1]))
325 pexit(
"real-time dimensions not matching expected");
333 strcpy(sdsname,
"ozone_mean");
335 if ((
rdsds(infileCLM, vgname, sdsname, inShape,
336 int_SDSdataCLMAVG)) != 0)
pexit(
"rdsds CLMAVG");
338 if ((inShape[0] != shape[0]) || (inShape[1] != shape[1]))
339 pexit(
"clim avg dimensions not matching expected");
341 strcpy(sdsname,
"ozone_std_dev");
343 if ((
rdsds(infileCLM, vgname, sdsname, inShape,
344 int_SDSdataCLMSTD)) != 0)
pexit(
"rdsds CLMSTD");
346 if ((inShape[0] != shape[0]) || (inShape[1] != shape[1]))
347 pexit(
"clim std dimensions not matching expected");
386 if (int_SDSdataRT[l + ipx] > 0) {
416 if (
i > min_lin) min_lin =
i;
417 if (
i < max_lin) max_lin =
i;
420 float_SDSdataQC[l] = 0.0;
424 if ((int_SDSdataRT[l] <= 0) || (int_SDSdataCLMAVG[l] <= 0)) {
426 if (int_SDSdataRT[l] <= 0) rtmiss++;
427 if (int_SDSdataRT[l] < minval) minval = int_SDSdataRT[l];
428 if (int_SDSdataRT[l] < loval) locnt++;
432 for (iv = 0; iv < 5; iv++) {
433 if (int_SDSdataRT[l] < lo_vals[iv]) lo_cnts[iv]++;
437 if ((!strcmp(gridflag,
"d")) || (!strcmp(gridflag,
"D"))) {
439 (float32) int_SDSdataRT[l] - int_SDSdataCLMAVG[l];
441 if (int_SDSdataCLMSTD[l] > 0)
442 float_SDSdataQC[l] = (float32)
443 (int_SDSdataRT[l] - int_SDSdataCLMAVG[l]) /
444 int_SDSdataCLMSTD[l];
448 if (int_SDSdataRT[l] < minval) minval = int_SDSdataRT[l];
449 if (int_SDSdataRT[l] > maxval) maxval = int_SDSdataRT[l];
452 if (int_SDSdataRT[l] < loval) flocnt += coslat;
453 if (int_SDSdataRT[l] > hival) fhicnt += coslat;
456 for (iv = 0; iv < 5; iv++) {
457 if (int_SDSdataRT[l] < lo_vals[iv]) lo_cnts[iv]++;
458 if (int_SDSdataRT[l] > hi_vals[iv]) hi_cnts[iv]++;
461 if (float_SDSdataQC[l] < -1.0) neg1cnt++;
462 if (float_SDSdataQC[l] < -2.0) neg2cnt++;
463 if (float_SDSdataQC[l] < -3.0) neg3cnt++;
465 if (float_SDSdataQC[l] > 1.0) pos1cnt++;
466 if (float_SDSdataQC[l] > 2.0) pos2cnt++;
467 if (float_SDSdataQC[l] > 3.0) pos3cnt++;
488 totalpts = pts_used - missing;
491 free(int_SDSdataCLMAVG);
492 free(int_SDSdataCLMSTD);
502 if ((
startHDF(outfile, &sdfid, &fid, DFACC_CREATE)) != 0)
503 pexit(
"Fatal error starting HDF file");
510 gridid =
setupGrid(fid,
"QC difference data");
519 pexit(
"Fatal error writing geometry");
525 strcpy(datalabel,
"QC array");
528 dataunit =
"std dev diff";
531 if ((sdsid =
wrtsds(sdfid,
rank, shape, datatype,
533 float_SDSdataQC)) < 0)
pexit(
"main wrtsds");
535 free(float_SDSdataQC);
541 if ((
result =
addAttr(sdsid, dataattr, DFNT_CHAR, dataunit)) != 0)
566 locnt = (
int) flocnt;
567 hicnt = (
int) fhicnt;
573 if (!strcmp(gridflag,
"s")) {
576 "---------------------------------------------------------------\n");
577 printf(
"Results of comparison of real-time and climatological files:\n");
578 printf(
"%s\t%s\n", infileRT, infileCLM);
579 printf(
"Month: %s\tParameter: %s\n", monthstr,
"Ozone");
580 printf(
"Thresholds: %f %f %f Max Missings: %d\n",
581 stdlimit1, stdlimit2, stdlimit3, maxmissing);
582 printf(
"Lines considered range from %d (lat %f) to %d (lat %f)\n",
583 max_lin, max_lat, min_lin,
min_lat);
586 printf(
"Minimum value: %d Maximum: %d\n", minval, maxval);
587 printf(
"Total # non-missing values: %6d / considered: %d\n",
589 printf(
"Total # values <%d: %d >%d: %d allowed: %d",
590 loval, locnt, hival, hicnt, outmax);
591 if (locnt >= outmax || hicnt >= outmax) {
601 printf(
"Total # points/percentage < -1 STD: %6d (%5.2f percent)\n",
602 neg1cnt, (
float) neg1cnt / totalpts * 100.0);
603 printf(
"Total # points/percentage +/- 1 STD: %6d (%5.2f percent)",
604 totalpts - (neg1cnt + pos1cnt),
605 (
float) (totalpts - (pos1cnt + neg1cnt)) / totalpts * 100.0);
607 if (((
float) (totalpts - (pos1cnt + neg1cnt)) / totalpts * 100.0) <
615 printf(
"Total # points/percentage > 1 STD: %6d (%5.2f percent)\n",
616 pos1cnt, (
float) pos1cnt / totalpts * 100.0);
619 printf(
"Total # points/percentage < -2 STD: %6d (%5.2f percent)\n",
620 neg2cnt, (
float) neg2cnt / totalpts * 100.0);
621 printf(
"Total # points/percentage +/- 2 STD: %6d (%5.2f percent)",
622 totalpts - (neg2cnt + pos2cnt),
623 (
float) (totalpts - (pos2cnt + neg2cnt)) / totalpts * 100.0);
625 if (((
float) (totalpts - (pos2cnt + neg2cnt)) / totalpts * 100.0) <
633 printf(
"Total # points/percentage > 2 STD: %6d (%5.2f percent)\n",
634 pos2cnt, (
float) pos2cnt / totalpts * 100.0);
637 printf(
"Total # points/percentage < -3 STD: %6d (%5.2f percent)\n",
638 neg3cnt, (
float) neg3cnt / totalpts * 100.0);
639 printf(
"Total # points/percentage +/- 3 STD: %6d (%5.2f percent)",
640 totalpts - (neg3cnt + pos3cnt),
641 (
float) (totalpts - (pos3cnt + neg3cnt)) / totalpts * 100.0);
643 if (((
float) (totalpts - (pos3cnt + neg3cnt)) / totalpts * 100.0) <
651 printf(
"Total # points/percentage > 3 STD: %6d (%5.2f percent)\n",
652 pos3cnt, (
float) pos3cnt / totalpts * 100.0);
655 printf(
"Total # missing values: %6d\n", missing);
656 printf(
"Total # missing real-time: %6d, max allowed: %6d",
658 if (rtmiss > maxmissing) {
665 "---------------------------------------------------------------\n");
669 printf(
"++%6d+%6d+%6d+%6d+%6d+%6d+%6d+%6d+%6d+%6d+\n", lo_cnts[0],
670 lo_cnts[1], lo_cnts[2], lo_cnts[3], lo_cnts[4], hi_cnts[0],
671 hi_cnts[1], hi_cnts[2], hi_cnts[3], hi_cnts[4]);
672 printf(
"++%5.2f+%5.2f+%5.2f+%6d+%5d+%5d+%5d+%5d+%5d+%5d+%5d+%5d+\n",
673 (
float) (totalpts - (pos1cnt + neg1cnt)) / totalpts * 100.0,
674 (
float) (totalpts - (pos2cnt + neg2cnt)) / totalpts * 100.0,
675 (
float) (totalpts - (pos3cnt + neg3cnt)) / totalpts * 100.0,
676 rtmiss, z_st[0], z_en[0], z_st[1], z_en[1], z_st[2], z_en[2],
678 printf(
"++%5.2f+%5.2f+%5.2f+%5.2f+%5.2f+%5.2f+%5d+%5d+%6.2f+%6.2f+\n",
679 (
float) neg1cnt / totalpts * 100.0, (
float) neg2cnt / totalpts * 100.0,
680 (
float) neg3cnt / totalpts * 100.0, (
float) pos1cnt / totalpts * 100.0,
681 (
float) pos2cnt / totalpts * 100.0, (
float) pos3cnt / totalpts * 100.0,
682 max_lin, min_lin, max_lat,
min_lat);
684 "---------------------------------------------------------------\n");
692 printf(
"Threshold Status: FAILED\n");
695 printf(
"Threshold Status: SUCCESS\n");
709 printf(
"\n\nUsage:\n");
710 printf(
"\t%s <nrtfile><clmfile><outfile><monthstr>\n", argv[0]);
711 printf(
"\t <param><t1><t2><t3><maxmiss><loval><hival><outmax>\n");
712 printf(
"\t <min_lat><max_lat>[diff][type]\n");
713 printf(
"\nWhere:\n");
714 printf(
"\tnrtfile: Real-time file to process\n");
715 printf(
"\tclmfile: Climatology file\n");
716 printf(
"\toutfile: output QC file\n");
717 printf(
"\tmonthstr: Month in string form (i.e., January)\n");
718 printf(
"\t (only used for examining climatology file,\n\t otherwise derived from the input file name\n");
719 printf(
"\tparam: Parameter to run: OZONE ...only\n");
721 "\tt1: Threshold percent of pts w/in 1 STD DEV of climatology\n");
723 "\tt2: Threshold percent of pts w/in 2 STD DEV of climatology\n");
725 "\tt3: Threshold percent of pts w/in 3 STD DEV of climatology\n");
726 printf(
"\tmaxmiss: Max num of missing NRT points permitted.\n");
727 printf(
"\tloval: lowest acceptable value in NRT file\n");
728 printf(
"\thival: highest acceptable value in NRT file\n");
729 printf(
"\toutmax: # of points, weighted by cos( latitude )\n");
730 printf(
"\t that can be outside lo/hival before error set\n");
731 printf(
"\tmin_lat: minimum latitude to consider (default -90.).\n");
732 printf(
"\tmax_lat: maximum latitude to consider (default 90.).\n");
733 printf(
"\t NOTE that if min_lat = 95, the limits will be\n");
734 printf(
"\t determined using the julian day in the equation:\n");
735 printf(
"\t min_lat = -90. + zen + max_lat\n");
736 printf(
"\t max_lat = 90. - zen - max_lat\n");
737 printf(
"\t where max_lat = limit adjustment, default, 10.\n");
738 printf(
"\t and zen = -23. * cos( ( jday + 10 ) * 360 / 365 )\n");
739 printf(
"\t ** [optional] parameters follow. Preceeding args \n");
741 "\t are required for subsequently ones used (fill spaces).\n");
743 "\tdiff: [optional] Enter 'd' to see a simple difference output\n");
744 printf(
"\t grid rather than the default STD variance 's'.\n");
745 printf(
"\t Difference calculations do not product reports.\n");
747 "\ttype [optional] Enter 1 if NRT file is actually a CLM file.\n");
748 printf(
"\t This will allow the program to read the CLM\n");
749 printf(
"\t as of test of CLM to itself.\n");
751 printf(
"Example: \n\n");
753 printf(
"\t o3qc $SDSDEMO/S199407100_TOVS.OZONE.hdf \n");
754 printf(
"\t $SDSDATA/S19891991_TOMS.OZONE.hdf output.hdf \n");
755 printf(
"\t March OZONE 25 50 75 50\n\n");