/*************************************************************************** * Routines to read grib and other files for the GDS into structures used * by the other parts of the code. Output information is * * Output into structure modelVals from PGRB NCEP files : * * Atmospheric profile data (NCEP): * Pressure * Height * Temperature * Relative Humidity * Ozone * Surface Pressure * Surface Temperature * Land Flag * Ice fraction * Wind Speed (U and V) * 2m Air temperature * * Output into structure modelVals2 from SFLUXGRB NCEP files : * (Note a different structure due to difference spatial resolution of * the sflux files) * * SSI (Surface Solar Irradiance) from an sflux grib file (NCEP) * * Output into structure sstStr from rtg HR sst file * * Previous days SST from RTG High resolution SST data * * Output into structure aerosolStr from aerosol file * * Aerosol data from aerosol file * * Version History : * * 30 March 2006 - Original written : Jonathan Mittaz * 17 April 2006 - Tidy up code : Jonathan Mittaz * 17 April 2006 - Fix bug in aerosol latitude (-70,70) -> (-90,90) not mapped * in output. makes no difference to output level2p since * aerosol data is not interpolated, but changed for sake of * clarity : Jonathan Mittaz * 19 April 2006 - Fix bug in addNewHours which had wrong number of days for * year end. Also removed superfluous call to * getMonthFromDayNo : Jonathan Mittaz */ /* * The exact desicion tree as to which files to use for the different outputs * can be quite complex and determined by what is available at the time. * * For pgrb files simply the files that bracket the current time * This may involve using the previous 6 hour files ie. * * Current time (hours) : 13.00Z * Could use: * * gfs.12z.pgrbf00 and gfs.12z.pgrbf03 * * OR * * gfs.06z.pgrbf06 and gfs.06z.pgrbf09 * * depending on what is available * * For the insolation (Surface Solar Irradiance) it is a bit more complicated * * Here we are dependent on the f03 and f06 files. The f03 is a 3 hour total * and for f06 is a 6 hour total. So we have * * Current time (hours) : 13.00Z * Could use * * gfs.12z.sfluxgrbf03 if available * * OR * * gfs.06zsfluxgrbf09 * which covers the same time period and is a 3 hour average * * In the case of current time (hours) : 15.00Z may need to go back one main * time level for the 2 files * * After 12 hours has been written can get * * gfs.12z.sfluxgrbf03 and gfs.12z.sfluxgrbf06 * * If these are not there then have to get * * gfs.06z.sfluxgrbf09 and gfs.06z.sfluxgrbf12 * * Note that when reading the grib files account must be taken of the fact that * the writing order of entries can change so an inventory is made and parsed * for the required entries (getInventory and getInventory2). Also a * consistency check on array sizes is made just in case. * * Final file is the previous day SST grib file * */ /* * Routine names : * * read_in_aerosol_data : Read in aerosol files into struct aerosolStr * read_in_previous_sst : Read in previous days SST data into structure * void getFilenameSst : Get filename of best previous SST data file * get_sst_file : Get filename of best previous SST GRIB file * returns 1 on success, 0 on failure * convert_sst_file : Converts previous days SST GRIB file to binary * read_in_model_data2 : Read in the SFLUX file for the SSI data * getFilename2 : Get bect SFLUX file(s) for SSI * convert_sflux_file : convert GRIB sflux file to simple binary * getInventory2 : make inventory file for GRIB conversion to * select the required parameters from an SFLUX * file * get_sflux_file : get GRIB SFLUX file * read_in_model_data : read in the PGRB model file (atmospheric profiles * etc) * convert_pgrb_file : Convert a GRIB PGRB file to simple binary * get_pgrb_file : Get a PGRB GRIB model file * getFilename : get best file for model interpolation based on the * current date * getInventory : Get input file to extract required parameters from * the PGRB GRIB file * addNewHours : Add hours to the current time to get new * hour/month/day/year * subtractNewHours : Subtract hours to the current time to get new * hour/month/day/year * getDayNoFromMonthDay : Get the day number from a calender date * getMonthFromDayNo : Get month/day from a day number */ /* Include headers */ #include #include #include #include /* Note need to include read_grib.h as some of the structures/prototypes * are defined in this header file */ #include #include /* Aerosol structures */ /* These will help in parsing the aerosol file */ /* Note more than one structure due to the fact that the aerosol data is * stored in 2 files and there are also data and rowID lines in the data * file */ struct aerosolRowIdStr { int row; int spare; int spare2; unsigned char marker; char spare3[3]; int time; int day; int year; }; struct aerosolDataStr { short optThickness; short avGradient; short gradientXp; short gradientXm; short gradientYp; short gradientYm; char physDescriptor; char spare; unsigned char nObs; unsigned char age; short weight; short coverage; char covarianceXp; char covarianceXm; char covarianceYp; char covarianceYm; short temp; short spare2; }; /* Header record taken from documentation */ struct aerosolDocStr { int start_record; float min_lat; float max_lat; float min_long; float max_long; float resolution; int spare[26]; int num_rows; int num_cols; int spare2[115]; int young_year; int young_month; int young_day; int young_hour; int old_year; int old_month; int old_day; int old_hour; int spare3; char rest_of_record[9476]; }; /* Taken from IDL code */ struct aerosolDocStrOld { short spare1; /* 0 */ short spare2; /* 1 */ short old_year; /* 2 */ short old_month; /* 3 */ short old_day; /* 4 */ short old_hour; /* 5 */ short young_year; /* 6 */ short young_month; /* 7 */ short young_day; /* 8 */ short young_hour; /* 9 */ short spare3; /* 10 */ short spare4; /* 11 */ short spare5; /* 12 */ short spare6; /* 13 */ short spare7; /* 14 */ short spare8; /* 15 */ short num_cols; /* 16 */ short num_rows; /* 17 */ short max_lat; /* 18 */ short max_long; /* 19 */ short min_lat; /* 20 */ short min_long; /* 21 */ short resolution; /* 22 */ short spare9; /* 23 */ short offset; /* 24 */ short spare10; /* 25 */ short spare11; /* 26 */ short spare12; /* 27 */ }; /* Prototypes */ void getFilename( char *pFilename, int Year, int Month, int Day, int Hours, int Type ); void getFilename2( char *pFilename, char *pFilename2, int Year, int Month, int Day, int Hours ); void getFilename2old( char *pFilename, char *pFilename2, int Year, int Month, int Day, int Hours ); void getFilenameSst( char *pFilename, int Year, int Month, int Day ); void addNewHours( int *pYearOut, int *pMonthOut, int *pDayOut, int *pNewHours, int hourStep ); void getInventory( void ); void getInventory2( void ); int get_pgrb_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pCheckFilename ); int get_sflux_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours ); void convert_sflux_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pFilename, char *pFilenameOut ); int get_sst_file( int *pNewYear, int *pNewMonth, int *pNewDay, char *pFilename ); void convert_sst_file( int newYear, int newMonth, int newDay, char *pFilename ); void convert_pgrb_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pFilename, char *pFilenameOut ); int find_filename2_back( int *newYear, int *newMonth, int *newDay, int *mainHours, int suppHours, char *pFilename, int shift_back ); /***************************************************************************** * * AEROSOL * * Routines to read in aerosol data * */ /* Routine to read in the aerosol file if present */ /* Note that the current file has no date in the filename - it is stored * internally. Therefore, if the latest file has not been copied * the last file copies will always be there to be read */ void read_in_aerosol_data( int Year, int Month, int Day, float Hours, struct aerosolStr *pAerosol ) { int i = 0; int j = 0; /* Structures to help read the file */ struct aerosolDocStr doc; struct aerosolDataStr data; struct aerosolRowIdStr rowId; int dayNo = 0; int dayNo2 = 0; int dayDifference = 0; int hourDifference = 0; int nLongitude = 0; int nLatitude = 0; int nLongitudeNew = 0; int nLatitudeNew = 0; int offset = 0; int shift_offset = 0.; float *pAerosolData = NULL; float *pAerosolDate = NULL; FILE *fpScript = NULL; char filename1[MAX_STRING_LENGTH]; char filename2[MAX_STRING_LENGTH]; unsigned char test_hdr[10108]; unsigned char test_data[28]; /* Get Day No for sst data */ getDayNoFromMonthDay( Year, Month, Day, &dayNo ); /* Get the aerosol data file */ strcpy(filename1,DATA_DIR); strcat(filename1,"aerosol/aerosoldev.field"); /* Now have data - read it */ if( NULL == (fpScript = fopen(filename1,"r")) ) printError("L2P","Cannot get aerosol file",NULL,NULL); /* Read document header file */ if( 1 != fread(&doc,sizeof(struct aerosolDocStr),1,fpScript) ) printError("L2P","Aerosol read error",NULL,NULL); /* Check byte swap etc. */ if( 2 != byteSwapInt(doc.start_record) ){ fprintf(stderr,"start record : %d\n",doc.start_record); printError("L2P","Aerosol data header byteswap error",NULL,NULL); } /* Get positional information */ pAerosol->minLatitude = byteSwapFlt(doc.min_lat); pAerosol->maxLatitude = byteSwapFlt(doc.max_lat); pAerosol->stepLatitude = byteSwapFlt(doc.resolution); pAerosol->minLongitude = byteSwapFlt(doc.min_long); pAerosol->maxLongitude = byteSwapFlt(doc.max_long); pAerosol->stepLongitude = byteSwapFlt(doc.resolution); /* Get size of stored array */ nLatitude = byteSwapInt(doc.num_rows); /* Note ncols includes an ID column */ nLongitude = byteSwapInt(doc.num_cols)-1; /* Because the aerosol data is defined from -70 to +70 (lat) redefine things * to force the aerosol data to match the other data needed in level 2p * which goes from -90 to +90 (lat) */ nLongitudeNew = 360/pAerosol->stepLongitude; nLatitudeNew = 181/pAerosol->stepLatitude; pAerosol->nDims1 = nLongitudeNew; pAerosol->nDims2 = nLatitudeNew; #if 0 /* Offset in data file to start of data */ offset = byteSwapShort(doc.offset); #endif /* One of the outputs of level 2p is the time difference between the * aerosol measurement and the SST time */ /* Time of last used data */ pAerosol->Year = byteSwapInt(doc.young_year); pAerosol->Month = byteSwapInt(doc.young_month); pAerosol->Day = byteSwapInt(doc.young_day); pAerosol->Hour = byteSwapInt(doc.young_hour); /* Now read data section */ /* Allocate output structure */ if( NULL == (pAerosol->pAerosol = (float *) calloc(nLongitudeNew*nLatitudeNew,sizeof(float))) ) printError("L2P","Cannot allocate pAerosol",NULL,NULL); if( NULL == (pAerosol->pAerosolDate = (float *) calloc(nLongitudeNew*nLatitudeNew,sizeof(float))) ) printError("L2P","Cannot allocate pAerosol",NULL,NULL); memset(pAerosol->pAerosol,'\0',nLongitudeNew*nLatitudeNew*sizeof(float)); /* Set a flag for no data */ for(i=0;ipAerosolDate+i) = -1; #if 0 memset(pAerosol->pAerosolDate,'\0',nLongitudeNew*nLatitudeNew*sizeof(float)); #endif /* Allocate array for temp storage - this keep the data as stored in the file which is only -70 to +70 in latitude */ if( NULL == (pAerosolData = (float *) calloc(nLongitude*nLatitude,sizeof(float))) ) printError("L2P","Cannot allocate pAerosol",NULL,NULL); if( NULL == (pAerosolDate = (float *) calloc(nLongitude*nLatitude,sizeof(float))) ) printError("L2P","Cannot allocate pAerosol",NULL,NULL); /* Get Day No in year for youngest time */ getDayNoFromMonthDay( pAerosol->Year, pAerosol->Month, pAerosol->Day, &dayNo2 ); /* Not get difference between SST date and last Aerosol date in days */ /* Differences in individual locations are referenced to this time */ if( Year == pAerosol->Year ) dayDifference = dayNo - dayNo2; else { if( (0 == pAerosol->Year%4 && 0 != pAerosol->Year%100) || 0 == pAerosol->Year%400 ) dayDifference = dayNo + (366 - dayNo2); else dayDifference = dayNo + (365 - dayNo2); } /* Convert to hours */ dayDifference *= 24; /* Hour difference */ hourDifference = ((int) Hours) - pAerosol->Hour; /* Read in data */ for(i=0;i *(pAerosolData+j+i*nLongitude)) *(pAerosolData+j+i*nLongitude) = 0.; /* date.age is time in hours from minimun time. Need to work out * how long this is from observation */ *(pAerosolDate+j+i*nLongitude) = (float) (data.age + dayDifference + hourDifference); } /* After each set of rows - there is a rowID block - here used as a * check as to if we're reading the file correctly */ /* Check for rowId */ if( 1 != fread(&rowId,sizeof(struct aerosolRowIdStr),1,fpScript) ) printError("L2P","Aerosol rowId error",NULL,NULL); /* Check marker */ /* Copy data over */ if( 255 != rowId.marker ){ printf("Marker : %d\n",rowId.marker); printError("L2P","Incorrect read of aerosol data",NULL,NULL); } } /* Now read in all the data - close file */ fclose(fpScript); /* Calculate shift in output array to take into account the difference * between -90,90 (output) and -70,70 (input) latitude ranges * Note at this time pAerosol->minLatitude contains the file value (-70) * rather that the output value (90) */ shift_offset = ((90 + pAerosol->minLatitude) / pAerosol->stepLongitude) * nLongitude; /* Reset min/max values */ /* This was not in code before April 17 2005 */ pAerosol->minLatitude = -90.; pAerosol->maxLatitude = +90; /* Now re-jigg data so it matches model ranges ie. from 0 - 360 rather than -180,180 longitude */ offset = (0 - pAerosol->minLongitude) / pAerosol->stepLongitude; /* Reset min/max values */ pAerosol->minLongitude += offset * pAerosol->stepLongitude; pAerosol->maxLongitude += offset * pAerosol->stepLongitude; for(i=0;ipAerosol+j+offset+i*nLongitude + shift_offset) = *(pAerosolData+j+i*nLongitude); *(pAerosol->pAerosolDate+j+offset+i*nLongitude + shift_offset) = *(pAerosolDate+j+i*nLongitude); } else { *(pAerosol->pAerosol+j-offset+i*nLongitude + shift_offset) = *(pAerosolData+j+i*nLongitude); *(pAerosol->pAerosolDate+j-offset+i*nLongitude + shift_offset) = *(pAerosolDate+j+i*nLongitude); } } } #ifdef TEST printf("Date on file : %d %d %d %d\n",pAerosol->Year,pAerosol->Month, pAerosol->Day,pAerosol->Hour); fpScript = fopen("test.dat","w"); fwrite(&nLongitude,sizeof(nLongitude),1,fpScript); fwrite(&nLatitude,sizeof(nLatitude),1,fpScript); fwrite(pAerosolData,sizeof(float),nLongitude*nLatitude, fpScript); fclose(fpScript); fpScript = fopen("test2.dat","w"); fwrite(&nLongitudeNew,sizeof(nLongitude),1,fpScript); fwrite(&nLatitudeNew,sizeof(nLatitude),1,fpScript); fwrite(pAerosol->pAerosol,sizeof(float),nLongitudeNew*nLatitudeNew, fpScript); fclose(fpScript); #endif /* Free memory for orignal data (opt depth and date) */ free(pAerosolData); free(pAerosolDate); } /***************************************************************************** * * PREVIOUS DAYS SST * * Routines to read in the previous days high resolution SST data product * */ /* Routine to read in the previous days SST from the RTG_HR_SST file * Note in this case there will be times when the previous calander day * SST hasn't be made yet, so then we have to go back to the day before * that */ #define NSST_LONG 4320 #define NSST_LAT 2160 void read_in_previous_sst( int Year, int Month, int Day, int *pNx, int *pNy, struct sstStr *pSstVals ) { FILE *fp = NULL; int ndata = 0; char filename[MAX_STRING_LENGTH]; char command[MAX_STRING_LENGTH]; /* Clear old version of SST of there */ if( 1 == pSstVals->Filled ){ free(pSstVals->pSST); pSstVals->pSST = NULL; pSstVals->Filled = 0; } /* Setup min/max's */ /* There are hardwired because the output of wgrib2 as a binary file does * not have header information in it */ pSstVals->minLongitude = 0.04166; pSstVals->maxLongitude = 359.95833; pSstVals->stepLongitude = 0.08333; pSstVals->minLatitude = 89.95833; pSstVals->maxLatitude = -89.95833; pSstVals->stepLatitude = -0.08333; /* Allocate arrays */ pSstVals->nDims1 = NSST_LONG; pSstVals->nDims2 = NSST_LAT; /* Now for this we have to work out if we need to read in 2 files or 1 */ /* SST */ if( NULL == (pSstVals->pSST = (float *) calloc(NSST_LONG*NSST_LAT,sizeof(float))) ) printError("L2P","Cannot allocate pSST",NULL,NULL); /* Set date of true SST level2p data for this file */ pSstVals->Year = Year; pSstVals->Month = Month; pSstVals->Day = Day; /* Now read binary file */ /* Get filename - this routines does the checking for previous versions etc */ getFilenameSst( filename, Year, Month, Day ); if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","Cannot open %s",NULL,filename); if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Current data file is currupted - needs to be recopied so remove it and stop processing of level 2p */ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSST ",NULL,filename); } if( ndata != NSST_LONG*NSST_LAT*sizeof(float) ) printError("L2P","sst data header out of range",NULL,NULL); if( NSST_LONG*NSST_LAT != fread(pSstVals->pSST,sizeof(float), NSST_LONG*NSST_LAT,fp)){ /* Remove bad file again */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pSST",NULL,NULL); } /* close file */ fclose(fp); pSstVals->Filled = 1; } /* Routine to get best previous day SST file and if needed convert * the GRIB format to simple binary to be read in */ void getFilenameSst( char *pFilename, int Year, int Month, int Day ) { int newYear = 0; int newMonth = 0; int newDay = 0; int dayNo = 0; char checkFilename[MAX_STRING_LENGTH]; FILE *fp = NULL; newYear = Year; newMonth = Month; newDay = Day; /* Get day no */ getDayNoFromMonthDay( newYear, newMonth, newDay, &dayNo ); /* Now go back one day */ dayNo -= 1; /* Re calculate the Year, Month, Day for this day no */ /* If gone back to previous year, then go back to 31 Dec */ if( 0 == dayNo ){ newYear -= 1; newMonth = 12; newDay = 31; } else getMonthFromDayNo( newYear, dayNo, &newMonth, &newDay ); /* Check to see if this file is there in converted format */ sprintf(checkFilename,"rtg/rtg_sst_grb_ht_0.083.%4.4d%2.2d%2.2d.dat", newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Copy over grib file and create dat (converted from GRIB) file */ /* Note routine returns a 1 on successfully finding the GRIB file 0 on failure */ if( 0 == get_sst_file(&newYear,&newMonth,&newDay, pFilename) ){ /* File not there - exit */ printError("L2P","Cannot access previous day SST",NULL,NULL); } /* Make readable binary version */ convert_sst_file( newYear, newMonth, newDay, pFilename ); } else { /* File is there - just need to output the filename and file will be * read by calling routine */ fclose(fp); } } /* Routine to find and return name of best GRIB file available */ /* Returns a 1 on successfully finding a GRIB file to use * Returns a 0 on no valid GRIB file found */ int get_sst_file( int *pNewYear, int *pNewMonth, int *pNewDay, char *pFilename ) { char checkFilename[MAX_STRING_LENGTH]; FILE *fpScript = NULL; int DayNo = 0; int newYear2 = 0; int newMonth2 = 0; int newDay2 = 0; int newYear = 0; int newMonth = 0; int newDay = 0; newYear = *pNewYear; newMonth = *pNewMonth; newDay = *pNewDay; /* File should be there */ sprintf(checkFilename,"rtg/rtg_sst_grb_hr_0.083.%4.4d%2.2d%2.2d", newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fpScript = fopen(pFilename,"r")) ){ /* Because of when the files are copied it's possible we have to * use the previous previous day */ getDayNoFromMonthDay( newYear, newMonth, newDay, &DayNo ); newYear2 = newYear; /* Step back an extra day */ DayNo -= 1; if( 0 == DayNo ){ newYear2 -= 1; newMonth2 = 12; newDay2 = 31; } else getMonthFromDayNo( newYear2, DayNo, &newMonth2, &newDay2 ); sprintf(checkFilename,"rtg/rtg_sst_grb_hr_0.083.%4.4d%2.2d%2.2d", newYear2,newMonth2,newDay2); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); *pNewYear = newYear2; *pNewMonth = newMonth2; *pNewDay = newDay2; if( NULL == (fpScript = fopen(pFilename,"r")) ){ /* No GRIB file available either 1 or 2 days old */ return 0; } else { /* Successfully found 2 day onld GRIB file */ fclose(fpScript); return 1; } } else { /* Successfully found 1 day onld GRIB file */ fclose(fpScript); return 1; } } /* Routine to convert intput previous days SST GRIB file to simple * binary format */ void convert_sst_file( int newYear, int newMonth, int newDay, char *pFilename ) { char command[MAX_STRING_LENGTH]; char runCommand[MAX_STRING_LENGTH]; char checkFilename[MAX_STRING_LENGTH]; char checkFilename2[MAX_STRING_LENGTH]; char invDir[MAX_STRING_LENGTH]; FILE *fpScript = NULL; /* Get executable directory */ strcpy(runCommand,WGRIB2_EXEC); /* Get inventory */ sprintf(checkFilename,"rtg/rtg_sst_grb_hr_0.083.%4.4d%2.2d%2.2d", newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); strcpy(invDir,DATA_DIR); strcat(invDir,"work/new-inventory.txt"); sprintf(command,"%s -s %s > %s",runCommand,checkFilename2,invDir); if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); sprintf(checkFilename,"rtg/rtg_sst_grb_hr_0.083.%4.4d%2.2d%2.2d.dat", newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* * Note extract data in same sense as defined by grid values */ sprintf(command, "%s -i -lola 0.04166:%d:0.08333 89.95833:%d:-0.08333 %s bin %s < %s >& /dev/null", runCommand,NSST_LONG,NSST_LAT,pFilename,checkFilename2,invDir); if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); if( NULL == (fpScript = fopen(pFilename,"r")) ) printError("L2P","Cannot open file %s",NULL,pFilename); fclose(fpScript); } /***************************************************************************** * * SURFACE SOLAR IRRADIANCE * * Routines to read in the SSI data from the SFLUX GRIB file * */ /* Reads in the correct set of SFLUX files to get a 3 hour average SSI * * This means that if we are in the 0-3 hour forecast range, use the f03 data * but if we are inthe 3-6 hour forecast range, we need to use the f06 and * f03 data to get the correct 3 hour average since the f06 file contains * the 6 hour average, not the 3 hour average. * * The routine also has to deal with the cases where the f03/f06 files * for the current time have not been copied yet, and we have to go back * to the previous forcast and use f09/f12 instead. */ void read_in_model_data2( int Year, int Month, int Day, int Hour, int *pNx, int *pNy, struct modelVals2 *pModelVals ) { int i = 0; FILE *fp = NULL; int newHours = 0; int ndata = 0; float *pSSI1 = NULL; float *pSSI2 = NULL; char filename[MAX_STRING_LENGTH]; char filename2[MAX_STRING_LENGTH]; char command[MAX_STRING_LENGTH]; /* If data already filled in structure - remove and reset */ if( 1 == pModelVals->Filled ){ free(pModelVals->pSurfaceSolarIrradiance); pModelVals->pSurfaceSolarIrradiance = NULL; pModelVals->Filled = 0; } /* Setup min/max's */ /* Hard wired because simple binary version of GRIB does not have a header */ /* pModelVals->minLongitude = 0; pModelVals->maxLongitude = 359.687; pModelVals->stepLongitude = 0.313; pModelVals->minLatitude = 89.761; pModelVals->maxLatitude = -89.761; pModelVals->stepLatitude = 0.313; */ pModelVals->minLongitude = 0; pModelVals->maxLongitude = 360; pModelVals->stepLongitude = 1.; pModelVals->minLatitude = -90.; pModelVals->maxLatitude = 90.; pModelVals->stepLatitude = 1.; /* Allocate arrays */ pModelVals->nDims1 = NLONG2; pModelVals->nDims2 = NLAT2; /* Now for this we have to work out if we need to read in 2 files or 1 */ /* Allocate structure memory */ /* SSIe */ if( NULL == (pModelVals->pSurfaceSolarIrradiance = (float *) calloc(NLONG2*NLAT2,sizeof(float))) ) printError("L2P","Cannot allocate pSurfaceSolarIrradiance", NULL,NULL); /* Get filename */ /* Need to know which 6 hour bin we're in, and which 3 hour - remember that * NCEP model files are labelled with a 6 hour bin z time and a forecase * window f0,f3,f6 etc. */ pModelVals->Year = Year; pModelVals->Month = Month; pModelVals->Day = Day; /* Get hour to nearest 3 hour (minimum needed bin) size */ newHours = (Hour/3)*3; pModelVals->Hour = newHours; /* Now read binary file */ /* Note getFilename2 returns 2 filenames. If filename2 = "IGNORE" then * we only need to use the f03 file alond. If it has a valid filename * then we are in the 3-6 hour regime and need to recalculate the SSI */ getFilename2( filename, filename2, Year, Month, Day, Hour ); /* Open f03 version of the file */ /* Note that if the file is corrupt, it is removed and will be recreated * in the next run */ if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","Cannot open %s",NULL,filename); /* SSI1 */ if( NULL == (pSSI1 = (float *) calloc(NLONG2*NLAT2,sizeof(float))) ) printError("L2P","Cannot allocate pSSI1",NULL,NULL); if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSSI",NULL,NULL); } if( ndata != NLONG2*NLAT2*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","sflux data header out of range",NULL,NULL); } if( NLONG2*NLAT2 != fread(pSSI1,sizeof(float),NLONG2*NLAT2,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pSSI1",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSSI",NULL,NULL); } /* close file */ fclose(fp); /* If we need to read in the f06 file, read it in. This is determined by * the string IGNORE being present or not */ if( NULL == strstr(filename2,"IGNORE") ){ if( NULL == (fp = fopen(filename2,"r")) ) printError("L2P","Cannot open %s",NULL,filename2); /* SSI2 */ if( NULL == (pSSI2 = (float *) calloc(NLONG2*NLAT2,sizeof(float))) ) printError("L2P","Cannot allocate pSSI2",NULL,NULL); if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename2); system(command); printError("L2P","Error reading header pSSI",NULL,NULL); } if( ndata != NLONG2*NLAT2*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename2); system(command); printError("L2P","sflux data header out of range",NULL,NULL); } if( NLONG2*NLAT2 != fread(pSSI2,sizeof(float),NLONG2*NLAT2,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename2); system(command); printError("L2P","Error reading pSSI2",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename2); system(command); printError("L2P","Error reading header pSSI",NULL,NULL); } /* close file */ fclose(fp); /* Because we have a 6 hour average data *pSST2 and 3 hour average data * *pSST1, we can get the correct 3 hour average for the 3-6 hour regime * by removing the f03 3 hour average from the data an recomputing * * 3-6 hourly av = f06 * 6 hours - f03 * 3 hours * * 3-6 3 hour av = (f06 * 6 hours - f03 * 3 hours)/3 hours * * or * * 3-6 3 hour av = (f06 * 2 - f03 * 1) */ /* Create correct 3 hour average */ for(i=0;ipSurfaceSolarIrradiance+i) = *(pSSI2+i)*2 - *(pSSI1+i); /* Free memory */ free(pSSI2); } else { /* Second file not needed - just use the original f03 three hour average */ /* Copy data */ for(i=0;ipSurfaceSolarIrradiance+i) = *(pSSI1+i); } free(pSSI1); } /* Set filled flag */ pModelVals->Filled = 1; } /* Routine to get the files needed to calculate the SSI from the SFLUX * Files */ void getFilename2old( char *pFilename, char *pFilename2, int Year, int Month, int Day, int Hours ) { /* Stuff to get filenames given year/dayno/time */ int mainHours = 0; int suppHours = 0; int newHours = 0; int newYear = 0; int newMonth = 0; int newDay = 0; int dayNo = 0; char checkFilename[MAX_STRING_LENGTH]; char checkFilename2[MAX_STRING_LENGTH]; int returnVal = 0; FILE *fp = NULL; /* Setup INGORE filenames */ strcpy(pFilename,"IGNORE"); strcpy(pFilename2,"IGNORE"); /* Create filenames */ newHours = (Hours/3)*3; /* Dependent on if before 03 or 06 */ newYear = Year; newMonth = Month; newDay = Day; /* Which 6 hour window are we in */ mainHours = (newHours/6)*6; /* Tells us which regime (f00, f03, f06 etc.) we are in */ suppHours = newHours - mainHours; if( 0 == suppHours ){ /* Within 0-3 so can only access f03 if possible */ /* If not latest f03, then get f09 from previous 6 hour model */ returnVal = 0; suppHours = 3; /* First try and find most recent (but may not be there yet) */ /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); fp = NULL; if( NULL == (fp = fopen(pFilename,"r")) ){ /* Converted file not there - see if GRIB file is there and convert */ returnVal = get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours); if( 0 == returnVal ){ /* 03 not there so have to get 03 from previous 6 hour slot */ /* Shift back main time by 6 hours and look for f06 file */ mainHours -= 6; /* Check to see if we've gone into the previous day */ if( 0 > mainHours ){ /* Re-set main hours to 24-6 */ mainHours = 18; /* SHift back ucrrent day to yesterday */ getDayNoFromMonthDay( newYear, newMonth, newDay, &dayNo ); dayNo -= 1; if( 0 == dayNo ){ newYear -= 1; newMonth = 12; newDay = 31; } else { getMonthFromDayNo( newYear, dayNo, &newMonth, &newDay ); } } /* Now set forecast hours to f09 */ suppHours = 9; /* Check to see if there's already a converted binary file version */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Copy over grib file and create dat file */ /* If zero returned than GRIB file not there so exit */ returnVal = get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours); /* 2nd time - look for 15 hour forecast */ if( 0 == returnVal ){ /* 03 not there so have to get 03 from previous 6 hour slot */ /* Shift back main time by 6 hours and look for f06 file */ mainHours -= 6; /* Check to see if we've gone into the previous day */ if( 0 > mainHours ){ /* Re-set main hours to 24-6 */ mainHours = 18; /* SHift back ucrrent day to yesterday */ getDayNoFromMonthDay( newYear, newMonth, newDay, &dayNo ); dayNo -= 1; if( 0 == dayNo ){ newYear -= 1; newMonth = 12; newDay = 31; } else { getMonthFromDayNo( newYear, dayNo, &newMonth, &newDay ); } } /* Not set forecast hours to f15 */ suppHours = 15; /* Check to see if there's already a converted binary file version */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Copy over grib file and create dat file */ /* If zero returned than GRIB file not there so exit */ if( 0==get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } } } /* If just copied over file need to convert it */ if( NULL == fp ){ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2,pFilename ); } else fclose(fp); /* Note that the 09 one contains the 3 hour average so it's the same as 03 from that point of view */ } else { fclose(fp); } } else { sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename ); } } } else if( 3 == suppHours ){ /* Within 3-6 so use f03 and f06 if there*/ /* If not latest f03,f06, then get f09,f12 from previous 6 hour model */ suppHours = 3; /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Main if - check for f03 .dat * simple binary file */ /* Copy over grib file and create dat file */ if( 0 == get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ /* Now must get f09 and f12 files */ mainHours -= 6; if( 0 > mainHours ){ mainHours = 18; getDayNoFromMonthDay( newYear, newMonth, newDay, &dayNo ); dayNo -= 1; if( 0 == dayNo ){ newYear -= 1; newMonth = 12; newDay = 31; } else { getMonthFromDayNo( newYear, dayNo, &newMonth, &newDay ); } } suppHours = 9; sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Copy over grib file and create dat file */ if( 0==get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2,pFilename ); } else fclose(fp); /* Now get f12 */ suppHours = 12; sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); if( NULL == (fp = fopen(pFilename2,"r")) ){ /* Copy over grib file and create dat file */ if( 0==get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2,pFilename2 ); } else fclose(fp); } else { /* f03 GRIB file is there - convert to simple binary format */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename ); /* Now the f06 */ suppHours = 6; /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); /* Note that if it is there, the already in correct output string var*/ strcat(pFilename2,checkFilename); if( NULL == (fp = fopen(pFilename2,"r")) ){ /* Copy over grib file and create dat file */ if( 0 == get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename2 ); } } } else { /* main If - f03 .dat (simple bianry file) is there - look for * f06 file */ suppHours = 6; /* Check to see if this file is there in simple binary format */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); if( NULL == (fp = fopen(pFilename2,"r")) ){ /* f06 .dat file not there - get and convert GRIB version */ if( 0==get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } /* Get input and output filenames for conversion */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2,pFilename2 ); } } } } /* Routine to get the files needed to calculate the SSI from the SFLUX * Files */ void getFilename2( char *pFilename, char *pFilename2, int Year, int Month, int Day, int Hours ) { /* Stuff to get filenames given year/dayno/time */ int mainHours = 0; int suppHours = 0; int newHours = 0; int newYear = 0; int newMonth = 0; int newDay = 0; int dayNo = 0; char checkFilename[MAX_STRING_LENGTH]; char checkFilename2[MAX_STRING_LENGTH]; int returnVal = 0; FILE *fp = NULL; /* Setup INGORE filenames */ strcpy(pFilename,"IGNORE"); strcpy(pFilename2,"IGNORE"); /* Create filenames */ newHours = (Hours/3)*3; /* Dependent on if before 03 or 06 */ newYear = Year; newMonth = Month; newDay = Day; /* Which 6 hour window are we in */ mainHours = (newHours/6)*6; /* Tells us which regime (f00, f03, f06 etc.) we are in */ suppHours = newHours - mainHours; if( 0 == suppHours ){ /* Within 0-3 so can only access f03 if possible */ /* If not latest f03, then get f09 from previous 6 hour model */ returnVal = 0; suppHours = 3; /* First try and find most recent (but may not be there yet) */ /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); fp = NULL; if( NULL == (fp = fopen(pFilename,"r")) ){ /* Converted file not there - see if GRIB file is there and convert */ returnVal = get_sflux_file(newYear,newMonth,newDay,mainHours, suppHours); if( 0 == returnVal ){ /* Now set forecast hours to f09 */ suppHours = 9; /* Shift back 6 hours */ returnVal = find_filename2_back( &newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename, 1 ); if( 0 == returnVal ){ /* Now set forecast hours to f15 */ suppHours = 15; /* Shift back 6 hours */ returnVal = find_filename2_back( &newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename, 1 ); if( 0 == returnVal ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } } } else { sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename ); } } else fclose(fp); } else if( 3 == suppHours ){ /* Within 3-6 so use f03 and f06 if there*/ /* If not latest f03,f06, then get f09,f12 from previous 6 hour model */ suppHours = 3; /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ /* Main if - check for f03 .dat * simple binary file */ /* Copy over grib file and create dat file */ if( 0 == get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ suppHours = 9; /* Shift back 6 hours */ returnVal = find_filename2_back(&newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename, 1); if( 0 == returnVal ){ suppHours = 15; /* Shift back 6 hours */ returnVal = find_filename2_back(&newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename, 1); if( 0 == returnVal ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } suppHours = 18; /* Do not shift back 6 hours */ returnVal = find_filename2_back(&newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename2, 0); if( 0 == returnVal ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } } else { suppHours = 12; /* Do not shift back 6 hours */ returnVal = find_filename2_back(&newYear, &newMonth, &newDay, &mainHours, suppHours, pFilename2, 0); if( 0 == returnVal ){ sprintf(checkFilename, "gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } } } else { /* f03 GRIB file is there - convert to simple binary format */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename ); /* Now the f06 */ suppHours = 6; /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); /* Note that if it is there, the already in correct output string var*/ strcat(pFilename2,checkFilename); if( NULL == (fp = fopen(pFilename2,"r")) ){ /* Copy over grib file and create dat file */ if( 0 == get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2, pFilename2 ); } } } else { /* main If - f03 .dat (simple bianry file) is there - look for * f06 file */ suppHours = 6; /* Check to see if this file is there in simple binary format */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); if( NULL == (fp = fopen(pFilename2,"r")) ){ /* f06 .dat file not there - get and convert GRIB version */ if( 0==get_sflux_file(newYear,newMonth,newDay,mainHours,suppHours) ){ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); printError("L2P","Cannot open file %s",NULL,checkFilename); } /* Get input and output filenames for conversion */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename2,DATA_DIR); strcat(pFilename2,checkFilename); /* Convert to binary grib file */ convert_sflux_file( newYear,newMonth,newDay,mainHours,suppHours, checkFilename2,pFilename2 ); } } } } int find_filename2_back( int *newYear, int *newMonth, int *newDay, int *mainHours, int suppHours, char *pFilename, int shift_back ) { int dayNo = 0; int returnVal = 0; char checkFilename[MAX_STRING_LENGTH]; char checkFilename2[MAX_STRING_LENGTH]; FILE *fp = NULL; /* 03 not there so have to get 03 from previous 6 hour slot */ /* Shift back main time by 6 hours and look for f06 file */ returnVal = 1; if( 1 == shift_back ){ *mainHours -= 6; } /* Check to see if we've gone into the previous day */ if( 0 > *mainHours ){ /* Re-set main hours to 24-6 */ *mainHours = 18; /* SHift back ucrrent day to yesterday */ getDayNoFromMonthDay( *newYear, *newMonth, *newDay, &dayNo ); dayNo -= 1; if( 0 == dayNo ){ *newYear -= 1; *newMonth = 12; *newDay = 31; } else { getMonthFromDayNo( *newYear, dayNo, newMonth, newDay ); } } /* Check to see if there's already a converted binary file version */ sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", *mainHours,suppHours,*newYear,*newMonth,*newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); #ifdef DEBUG printf("Looking for %s\n",checkFilename); #endif if( NULL == (fp = fopen(pFilename,"r")) ){ /* Copy over grib file and create dat file */ /* If zero returned than GRIB file not there so exit */ returnVal = get_sflux_file(*newYear,*newMonth,*newDay,*mainHours, suppHours); if( 0 != returnVal ){ sprintf(checkFilename, "ncep/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.grib2", *mainHours,suppHours,*newYear,*newMonth,*newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,checkFilename); sprintf(checkFilename, "work/gfs.t%2.2dz.sfluxgrbf%2.2d.%4.4d%2.2d%2.2d.dat", *mainHours,suppHours,*newYear,*newMonth,*newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Convert to binary grib file */ convert_sflux_file( *newYear,*newMonth,*newDay,*mainHours,suppHours, checkFilename2, pFilename ); } } else { #ifdef DEBUG printf("Found %s\n",checkFilename); #endif fclose(fp); } return(returnVal); } /* Converts sfile GRIB file (pFilename) to putput file simple binary * format stored in pFilenameOut */ void convert_sflux_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pFilename, char *pFilenameOut ) { char command[MAX_STRING_LENGTH]; char runCommand[MAX_STRING_LENGTH]; char invDir[MAX_STRING_LENGTH]; FILE *fpScript = NULL; /* If we want this code to do the GRIB2->.dat conversion then when compiling set -DCONVERT_GRIB. If not (default) then request to convert will result in an error since the conversion should have already been done elsewhere*/ #ifndef CONVERT_GRIB printError("L2P","convert_sflux_file: Attempt to convert GRIB2 to .dat not allowed in this version",NULL,NULL); #endif strcpy(runCommand,WGRIB2_EXEC); /* Get inventory */ strcpy(invDir,DATA_DIR); strcat(invDir,"work/inventory.txt"); sprintf(command, "%s -s %s > %s",runCommand,pFilename,invDir); if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); /* This routine parses the output inventory file and picks out the required * entry - must do it this way because GRIB files can be written in any * order * Output inventory list is written to new-invenotry.txt */ getInventory2(); strcpy(invDir,DATA_DIR); strcat(invDir,"work/new-inventory.txt"); sprintf(command,"%s -i -lola 0:360:1 -90:181:1 %s bin %s < %s >& /dev/null", runCommand,pFilenameOut,pFilename,invDir); #ifdef DEBUG printf("COMMAND: %s\n",command); #endif /* Run wgrib2.x to convert GRIB file to a simple binary format */ if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); /* Check to see if wgrib2 run above has worked */ if( NULL == (fpScript = fopen(pFilenameOut,"r")) ) printError("L2P","Cannot open file %s",NULL,pFilenameOut); fclose(fpScript); } /* Routine to parse inventory list of a GRIB file and make a new input * list to be used to strip out the relevant entries into a simple * binary format to be read in */ /* max lines is used by both getInventory and getInventory2 */ #define MAX_LINES 500 /* Define ehich entries to use - in this case just the SSI */ #define NLINES_DEFINED2 1 char chosenLines2[NLINES_DEFINED2][MAX_STRING_LENGTH] = {"DSWRF:surface"}; char chosenLines2_second[NLINES_DEFINED2][MAX_STRING_LENGTH] = {"hour ave fcst"}; void getInventory2( void ) { int i = 0; int j = 0; int nlines = 0; FILE *fp = NULL; char **pLines = NULL; char invDir[MAX_STRING_LENGTH]; /* Make lines array to read in inventory file */ if( NULL == (pLines = (char **) calloc(MAX_LINES,sizeof(char *))) ) printError("L2P","getInventory2: Cannot allocate pLines",NULL,NULL); for(i=0;iFilled ){ for(i=0;inLevels;i++){ free(*(pModelVals->ppHGT+i)); free(*(pModelVals->ppTemperature+i)); free(*(pModelVals->ppRelativeHumidity+i)); free(*(pModelVals->ppOzone+i)); } free(pModelVals->ppHGT); pModelVals->ppHGT = NULL; free(pModelVals->ppTemperature); pModelVals->ppTemperature = NULL; free(pModelVals->ppRelativeHumidity); pModelVals->ppRelativeHumidity = NULL; free(pModelVals->ppOzone); pModelVals->ppOzone = NULL; free(pModelVals->pPressure); pModelVals->pPressure = NULL; free(pModelVals->pSurfacePressure); pModelVals->pSurfacePressure = NULL; free(pModelVals->pSurfaceTemperature); pModelVals->pSurfaceTemperature = NULL; free(pModelVals->pLand); pModelVals->pLand = NULL; free(pModelVals->pIce); pModelVals->pIce = NULL; free(pModelVals->pTemperature2m); pModelVals->pTemperature2m = NULL; free(pModelVals->pWindSpeedU); pModelVals->pWindSpeedU = NULL; free(pModelVals->pWindSpeedV); pModelVals->pWindSpeedV = NULL; free(pModelVals->pTau4); pModelVals->pTau4 = NULL; pModelVals->Filled = 0; } /* Setup min/max's */ pModelVals->minLongitude = 0.; pModelVals->maxLongitude = 359.; pModelVals->stepLongitude = 1.; /* Need to check if it's this way up */ pModelVals->minLatitude = -90.; pModelVals->maxLatitude = 90.; pModelVals->stepLatitude = 1.; /* Allocate arrays */ pModelVals->nDims1 = NLONG; pModelVals->nDims2 = NLAT; pModelVals->nLevels = NUMBER_OF_LEVELS; /* First the profile data */ /* Pressure */ if( NULL == (pModelVals->pPressure = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pPressure",NULL,NULL); for(i=0;ipPressure+i) = pressureArray[i]; /* Height */ if( NULL == (pModelVals->ppHGT = (float **) calloc(NUMBER_OF_LEVELS, sizeof(float *))) ) printError("L2P","Cannot allocate ppHGT",NULL,NULL); for(i=0;ippHGT+i) = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate *(ppHGT+i)",NULL,NULL); } /* Temperature */ if( NULL == (pModelVals->ppTemperature = (float **) calloc(NUMBER_OF_LEVELS, sizeof(float *))) ) printError("L2P","Cannot allocate ppTemp",NULL,NULL); for(i=0;ippTemperature+i) = (float *) calloc(NLONG*NLAT,sizeof(float))) ) printError("L2P","Cannot allocate *(ppTemp+i)",NULL,NULL); } /* Relative Humidity */ /* Note not available for all levels - not the top (10-100 mb) */ if( NULL == (pModelVals->ppRelativeHumidity = (float **) calloc(NUMBER_OF_LEVELS,sizeof(float *))) ) printError("L2P","Cannot allocate ppRH",NULL,NULL); for(i=MIN_RELATIVE_HUMIDITY_LEVEL;ippRelativeHumidity+i) = (float *) calloc(NLONG*NLAT,sizeof(float))) ) printError("L2P","Cannot allocate *(ppRH+i)",NULL,NULL); } for(i=0;ippRelativeHumidity+i) = NULL; /* Ozone */ /* Note that not all levels are included */ if( NULL == (pModelVals->ppOzone = (float **) calloc(NUMBER_OF_LEVELS, sizeof(float *))) ) printError("L2P","Cannot allocate *pOzone",NULL,NULL); for(i=0;ippOzone+i) = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate *(pOzone+i)",NULL,NULL); } for(i=MAX_OZONE_LEVEL;ippOzone+i) = NULL; /* Now the single point location data */ /* Surface Pressure */ if( NULL == (pModelVals->pSurfacePressure = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pSurfacePressure",NULL,NULL); /* Surface Temp */ if( NULL == (pModelVals->pSurfaceTemperature = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pSurfaceTemp",NULL,NULL); /* Land */ if( NULL == (pModelVals->pLand = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pLand",NULL,NULL); /* Ice */ if( NULL == (pModelVals->pIce = (float *) calloc(NLONG*NLAT,sizeof(float))) ) printError("L2P","Cannot allocate pIce",NULL,NULL); /* 2m Temp */ if( NULL == (pModelVals->pTemperature2m = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pTemperature2m",NULL,NULL); /* windSpeedU */ if( NULL == (pModelVals->pWindSpeedU = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pWindSpeedU",NULL,NULL); /* windSpeedV */ if( NULL == (pModelVals->pWindSpeedV = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pWindSpeedV",NULL,NULL); /* Optical depth */ if( NULL == (pModelVals->pTau4 = (float *) calloc(NLONG*NLAT, sizeof(float))) ) printError("L2P","Cannot allocate pTau4",NULL,NULL); memset(pModelVals->pTau4,'\0',NLONG*NLAT*sizeof(float)); /* First for time get filename */ /* This routines does all the checking and getting/converting of GRIB files * etc */ getFilename( filename, Year, Month, Day, Hour, Type ); /* Set time in structure */ pModelVals->Year = Year; pModelVals->Month = Month; pModelVals->Day = Day; newHours = (Hour/3)*3; switch(Type){ case 1: /* First file */ pModelVals->Hour = newHours; break; case 2: /* Second file */ pModelVals->Hour = newHours+3; break; default: printError("L2P","Read model: Type of out range",NULL,NULL); break; } /* Now read binary file */ if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","Cannot open %s",NULL,filename); /* Note that compared to HDF file levels are reversed */ /* Since code was written for HDF files originally we must reverse * positions */ /* Not also that if there is any problem with reading the file it is * deleted and hopefully a corrected version will be remade next time */ /* This code could be cleaned up by writing a routine to do the read as * a generic function */ /* First read in height */ for(i=0;ippHGT+NUMBER_OF_LEVELS-i-1), sizeof(float),NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading ppHGT",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header ppHGT",NULL,NULL); } } /* Next read in temperature profile */ for(i=0;ippTemperature+NUMBER_OF_LEVELS-i-1), sizeof(float),NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading ppTemperature",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header ppTemperature",NULL,NULL); } } /* Next read in relative humidity */ for(i=0;ippRelativeHumidity+pos),sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading ppRelativeHumidity",NULL, NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header ppRelativeHumidity", NULL,NULL); } } /* Now read in Ozone */ for(i=0;ippOzone+i),sizeof(float),NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading ppOzone",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header ppOzone",NULL,NULL); } } /* Now read in the surface parameters - there are single nos on a long/lat * grid */ /* Surface pressure */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSurfacePressure",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pSurfacePressure", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pSurfacePressure,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pSurfacePressure",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSurfacePressure",NULL,NULL); } /* Convert to correct units for RTM model */ for(i=0;ipSurfacePressure+i) /= 100.; /* Surface temperature */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSurfaceTemperature",NULL, NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pSurfaceTemperature", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pSurfaceTemperature,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pSurfaceTemperature",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pSurfaceTemperature",NULL, NULL); } /* Land Flag */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pLand",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pLand", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pLand,sizeof(float),NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pLand",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pLand",NULL,NULL); } /* Ice Fraction */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pIce",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pIce", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pIce,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pIce",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pIce",NULL,NULL); } /* 2m temperature */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pTemperature2m",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pTemperature2m", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pTemperature2m,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pTemperature2m",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pTemperature2m",NULL,NULL); } /* Wind speed U */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pWindSpeedU",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pWindSpeedU", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pWindSpeedU,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pWindSpeedU",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pWindSpeedU",NULL,NULL); } /* Wind speed V */ if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pWindSpeedV",NULL,NULL); } if( ndata != NLONG*NLAT*sizeof(float) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P", "Number of entries (wgrib2) incorrect : pWindSpeedV", NULL,NULL); } if( NLONG*NLAT != fread(pModelVals->pWindSpeedV,sizeof(float), NLONG*NLAT,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading pWindSpeedV",NULL,NULL); } if( 1 != fread(&ndata,sizeof(int),1,fp) ){ /* Remove bad file */ sprintf(command,"rm -rf %s",filename); system(command); printError("L2P","Error reading header pWindSpeedV",NULL,NULL); } /* Set flag to say that structure has been filled */ pModelVals->Filled = 1; /* close file */ fclose(fp); } /* Convert to GRIB model file to a simple binary format using wgrib2.x */ void convert_pgrb_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pFilename, char *pFilenameOut ) { char command[MAX_STRING_LENGTH]; char runCommand[MAX_STRING_LENGTH]; char checkFilename[MAX_STRING_LENGTH]; char invDir[MAX_STRING_LENGTH]; FILE *fpScript = NULL; /* If we want this code to do the GRIB2->.dat conversion then when compiling set -DCONVERT_GRIB. If not (default) then request to convert will result in an error since the conversion should have already been done elsewhere*/ #ifndef CONVERT_GRIB printError("L2P","convert_pgrb_file: Attempt to convert GRIB2 to .dat not allowed in this version",NULL,NULL); #endif /* Check file is there */ if( NULL == (fpScript = fopen(pFilename,"r")) ) printError("L2P","Cannot find file %s",NULL,pFilename); /* Get executable directory */ strcpy(runCommand,WGRIB2_EXEC); /* Get inventory list */ strcpy(invDir,DATA_DIR); strcat(invDir,"work/inventory.txt"); sprintf(command,"%s -s %s > %s",runCommand,pFilename,invDir); if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); /* getInventory parses the inventory list and makes a list of only * those parameters we need - have to do it this way because the * entries in the grib file can be in any order */ getInventory(); /* Output filename */ sprintf(checkFilename, "work/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilenameOut,DATA_DIR); strcat(pFilenameOut,checkFilename); /* Use new inventory list to extract only those parameters we need */ strcpy(invDir,DATA_DIR); strcat(invDir,"work/new-inventory.txt"); sprintf(command,"%s -i -bin %s %s < %s >& /dev/null", runCommand,pFilenameOut,pFilename,invDir); #ifdef DEBUG printf("COMMAND: %s\n",command); #endif if( -1 == system(command) ) printError("L2P","wgrib2 failed",NULL,NULL); /* Check that output filename has been created */ if( NULL == (fpScript = fopen(pFilenameOut,"r")) ) printError("L2P","Cannot open file %s",NULL,pFilenameOut); fclose(fpScript); } /* Get the PGRB file */ int get_pgrb_file( int newYear, int newMonth, int newDay, int mainHours, int suppHours, char *pCheckFilename ) { char checkFilename2[MAX_STRING_LENGTH]; FILE *fpScript = NULL; /* File should be copied already */ /* Make filename to search for */ sprintf(pCheckFilename, "ncep/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.1p0deg", mainHours,suppHours,newYear,newMonth,newDay); strcpy(checkFilename2,DATA_DIR); strcat(checkFilename2,pCheckFilename); #ifdef DEBUG fprintf(stderr,"Attempting to read %s\n",pCheckFilename); #endif /* Check to see if file exists * * If file exists the return a 1 * If file does not exist, return a 0 */ if( NULL == (fpScript = fopen(checkFilename2,"r")) ){ #ifdef DEBUG printf("NOT Found %s\n",pCheckFilename); #endif return 0; } else { fclose(fpScript); return 1; } } /* Get the name of the best PGRB file (in converted binary format) to use */ /* Note this routine does all the getting and converting of the GRIB files * as part of it */ void getFilename( char *pFilename, int Year, int Month, int Day, int Hours, int Type ) { /* Stuff to get filenames given year/dayno/time */ int mainHours = 0; int suppHours = 0; int newHours = 0; int newYear = 0; int newMonth = 0; int newDay = 0; char checkFilename[MAX_STRING_LENGTH]; FILE *fp = NULL; FILE *fpScript = NULL; /* Create filenames */ /* Get hours at a resolution of 3 hours - smalled time step of models */ newHours = (Hours/3)*3; /* Dependent on type - before or after time */ /* Type = 1 means get the filejust before the current time * Type = 2 means get the file just after the current time */ switch( Type ){ case 1: newYear = Year; newMonth = Month; newDay = Day; mainHours = (newHours/6)*6; suppHours = newHours - mainHours; break; case 2: /* Only difference from case 1 is the inclusion of addNewHours which adds * 3 hours to time taking into accout day changes etc */ newYear = Year; newMonth = Month; newDay = Day; addNewHours( &newYear, &newMonth, &newDay, &newHours, 3 ); newHours = (newHours/3)*3; mainHours = (newHours/6)*6; suppHours = newHours - mainHours; break; default: printError("L2P","getFilename Type out of range",NULL,NULL); } /* Check to see if this file is there */ sprintf(checkFilename, "work/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); #ifdef DEBUG fprintf(stderr,"Attempting to read %s\n",checkFilename); #endif strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fp = fopen(pFilename,"r")) ){ #ifdef DEBUG fprintf(stderr,"Failed to read %s\n",checkFilename); fprintf(stderr,"Getting GRIB File\n"); #endif /* If converted file does not exist Copy over grib file and create dat * file */ if( 0 == get_pgrb_file(newYear,newMonth,newDay,mainHours,suppHours, checkFilename) ){ #ifdef DEBUG fprintf(stderr,"Failed to get %s\n",checkFilename); #endif /* File not there - back off to previous 6 hourly model */ newHours = mainHours; /* Subtract 6 hours from current time to get previous model */ subtractNewHours( &newYear, &newMonth, &newDay, &newHours, 6 ); /* To get back the the data time must add 6 the the f%% value */ suppHours += 6; mainHours = newHours; /* Check to see if we've copied the previous file before */ sprintf(checkFilename, "work/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); #ifdef DEBUG fprintf(stderr,"Attempting to read %s\n",checkFilename); #endif strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fpScript = fopen(pFilename,"r")) ){ /* Have not copied this file before - create from GRIB file */ /* Copy over grib file and create dat file */ if( 0 == get_pgrb_file(newYear,newMonth,newDay,mainHours,suppHours, checkFilename ) ){ /* File not there - back off to previous 6 hourly model */ newHours = mainHours; /* Subtract 6 hours from current time to get previous model */ subtractNewHours( &newYear, &newMonth, &newDay, &newHours, 6 ); /* To get back the the data time must add 6 the the f%% value */ suppHours += 6; mainHours = newHours; /* Check to see if we've copied the previous file before */ sprintf(checkFilename, "work/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.dat", mainHours,suppHours,newYear,newMonth,newDay); #ifdef DEBUG fprintf(stderr,"Attempting to read %s\n",checkFilename); #endif strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); if( NULL == (fpScript = fopen(pFilename,"r")) ){ /* Have not copied this file before - create from GRIB file */ /* Copy over grib file and create dat file */ if( 0 == get_pgrb_file(newYear,newMonth,newDay,mainHours,suppHours, checkFilename ) ) printError("L2P","Cannot find correct GRIB file %s",NULL, checkFilename); /* Make readable binary version */ /* Input filename */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.1p0deg", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Get new binary version of pgrb file */ convert_pgrb_file( newYear, newMonth, newDay, mainHours, suppHours, pFilename,checkFilename ); /* Output filename is stored in checkFilename - copy to output */ /* This should be tidied up as it's a bit silly with swapping * checkFilename/pFilename around */ strcpy(pFilename,checkFilename); } } else { /* Make readable binary version */ /* Input filename */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.1p0deg", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Get new binary version of pgrb file */ #ifdef DEBUG printf("Converting %s\n",checkFilename); #endif convert_pgrb_file( newYear, newMonth, newDay, mainHours, suppHours, pFilename,checkFilename ); /* Output filename is stored in checkFilename - copy to output */ /* This should be tidied up as it's a bit silly with swapping * checkFilename/pFilename around */ strcpy(pFilename,checkFilename); } } } else { /* Previous GRIB file exists - convert it to simple binary */ sprintf(checkFilename, "ncep/gfs.t%2.2dz.pgrb2f%2.2d.%4.4d%2.2d%2.2d.1p0deg", mainHours,suppHours,newYear,newMonth,newDay); strcpy(pFilename,DATA_DIR); strcat(pFilename,checkFilename); /* Make readable binary version */ convert_pgrb_file( newYear, newMonth, newDay, mainHours, suppHours, pFilename, checkFilename ); /* CheckFIlename contains output name - copy to output variable */ /* This should be tidied up as it's a bit silly with swapping * checkFilename/pFilename around */ strcpy(pFilename,checkFilename); } } else { /* Converted binary file already exists - just return */ fclose(fp); } } /* getInventory parses the inventory list and makes a list of only * those parameters we need - have to do it this way because the * entries in the grib file can be in any order */ /* These are the entries we need from the GRIB file */ #define NLINES_DEFINED 86 char chosenLines[NLINES_DEFINED][MAX_STRING_LENGTH] = {"HGT:1000 mb", /* Height profile */ "HGT:975 mb", "HGT:950 mb", "HGT:925 mb", "HGT:900 mb", "HGT:850 mb", "HGT:800 mb", "HGT:750 mb", "HGT:700 mb", "HGT:650 mb", "HGT:600 mb", "HGT:550 mb", "HGT:500 mb", "HGT:450 mb", "HGT:400 mb", "HGT:350 mb", "HGT:300 mb", "HGT:250 mb", "HGT:200 mb", "HGT:150 mb", "HGT:100 mb", "HGT:70 mb", "HGT:50 mb", "HGT:30 mb", "HGT:20 mb", "HGT:10 mb", "TMP:1000 mb", /* Temperature profile */ "TMP:975 mb", "TMP:950 mb", "TMP:925 mb", "TMP:900 mb", "TMP:850 mb", "TMP:800 mb", "TMP:750 mb", "TMP:700 mb", "TMP:650 mb", "TMP:600 mb", "TMP:550 mb", "TMP:500 mb", "TMP:450 mb", "TMP:400 mb", "TMP:350 mb", "TMP:300 mb", "TMP:250 mb", "TMP:200 mb", "TMP:150 mb", "TMP:100 mb", "TMP:70 mb", "TMP:50 mb", "TMP:30 mb", "TMP:20 mb", "TMP:10 mb", "RH:1000 mb", /* Relative Humidity profile */ "RH:975 mb", "RH:950 mb", "RH:925 mb", "RH:900 mb", "RH:850 mb", "RH:800 mb", "RH:750 mb", "RH:700 mb", "RH:650 mb", "RH:600 mb", "RH:550 mb", "RH:500 mb", "RH:450 mb", "RH:400 mb", "RH:350 mb", "RH:300 mb", "RH:250 mb", "RH:200 mb", "RH:150 mb", "RH:100 mb", "O3MR:100 mb", /* Ozone profile */ "O3MR:70 mb", "O3MR:50 mb", "O3MR:30 mb", "O3MR:20 mb", "O3MR:10 mb", "PRES:surface", /* Surface Pressure */ "TMP:surface", /* Surfece Tempereature */ "LAND:surface", /* Land flag */ "ICEC:surface", /* Ice fraction */ "TMP:2 m above ground", /* Temperature 2m above ground */ "UGRD:10 m above ground", /* Wind speed */ "VGRD:10 m above ground" }; void getInventory( void ) { int i = 0; int j = 0; int nlines = 0; FILE *fp = NULL; char **pLines = NULL; char invDir[MAX_STRING_LENGTH]; /* Make lines array */ if( NULL == (pLines = (char **) calloc(MAX_LINES,sizeof(char *))) ) printError("L2P","getInventory: Cannot allocate pLines",NULL,NULL); for(i=0;i *pNewHours ){ /* Get hours in previous day */ *pNewHours += 24.; /* Get day number of current time */ getDayNoFromMonthDay( *pYearOut, *pMonthOut, *pDayOut, &DayNo ); /* Go back a day */ DayNo -= 1; /* If dayNo equals zero need to go back a year */ if( 0 == DayNo ){ *pYearOut -= 1; *pMonthOut = 12; *pDayOut = 31; } else { /* Get new day/month/year combination */ getMonthFromDayNo( *pYearOut, DayNo, pMonthOut, pDayOut ); } } } /* Deal with dates etc */ /* Get day number from a date */ /* Note have to take into account leap years etc */ void getDayNoFromMonthDay( int Year, int Month, int Day, int *pDayNo ) { int DayToMonth[13] = {0,31,59,90,120,151,181,212,243,273,304,334,365}; int DayToMonthLeap[13] = {0,31,60,91,121,152,182,213,244,274,305,335,366}; if( (0 == Year%4 && 0 != Year%100) || 0 == Year%400 ) *pDayNo = DayToMonthLeap[Month-1] + Day; else *pDayNo = DayToMonth[Month-1] + Day; } /* Get the month in the year dependent on the day number */ /* Note have to take into account leap years etc */ void getMonthFromDayNo( int Year, int DayNo, int *pMonth, int *pDay ) { int i = 0; int DayToMonth[13] = {0,31,59,90,120,151,181,212,243,273,304,334,365}; int DayToMonthLeap[13] = {0,31,60,91,121,152,182,213,244,274,305,335,366}; if( (0 == Year%4 && 0 != Year%100) || 0 == Year%400 ){ /* If we are in a leap year - use appropriate array */ for(i=1;i<13;i++){ if( DayNo < DayToMonthLeap[i] + 1 ){ *pMonth = i; *pDay = DayNo - DayToMonthLeap[i-1]; break; } } } else { /* Not a leap year */ for(i=1;i<13;i++){ if( DayNo < DayToMonth[i] + 1 ){ *pMonth = i; *pDay = DayNo - DayToMonth[i-1]; break; } } } /* If leap year loop */ }