/* * Routine to make the level2p dataset for the GODAE project * * Reads in SST image, clear sky probability and gets and read * nearest in time NCEP models, aerosol dataset and previous day SST * files. Gets optical depth of 11 micron filter from radiaitive * transfer based on NCEP models and calculates biases * * * Written by Jonathan Mittaz * * Version History * * 30 March 2006 - Original version : Jonathan Mittaz * 19 April 2006 - Tidied up code : Jonathan Mittaz * 19 April 2006 - Fixed error in non-interpolated surface temperature value - * was incorrectly using the 2m temperature * 20 June 2007 - only output parameters when there is clear sku ocean - * should help compression * 02 March 2009 - Fixed leap year problem * 17 March 2009 - Added in code to deal with MTSAT and MSG (names/headers etc) * waiting for update to SSESs * 26 Jan 2012 - Added in code to deal with GOES-15 and made modifications * in preparation to the transition to GHRSST GDS 2.0 * 13 Apr 2012 - Added in new SSES code using the MTLS error flag as input * 9 Jul 2012 - Added in diurnal warming * 26 Sep 2013 - Modifications for RTG in GRIB2 * */ /* include files */ #include #include #include #include #include #include #include /* Include files from other parts of code */ #include #include #include /* Define some macros */ #define myMIN(a,b) ((a < b) ? a : b) #define myMAX(a,b) ((a > b) ? a : b) /* Define success and fail for gentemann diurnal model */ #define SUCCESS 0 #define FAIL -1 #define BADVAL -999. /* Contain model at times before and after current time */ /* contains the model data from before the corrent time */ struct modelVals model1; /* contains the model data from after the corrent time */ struct modelVals model2; /* contains the model data (SFLUX) for the SSI */ struct modelVals2 model1s; /* Contains the aerosol data */ struct aerosolStr aerosol; /* Contains the previous days SST data */ struct sstStr previousSst; /* Diurnal warming model from Gentemann et al. 2003 - for each model (1/2) */ float *pDiurnal1 = NULL; float *pDiurnal2 = NULL; /* Prototypes */ int main( int argc, char *argv[] ); void read_in_all_model_data( int satId, int Year, int Month, int Day, float utc, const float longitude, const float latitude, float *pWindSpeed, float *pIce, int *pLandCode, float *pTau4, float *pTemp2m, float *pSSI, float *pTimeDiff, float *pAerOptDepth, float *pAerOptDepthTime, float *pPreviousSST, char *pFile1, char *pFile1s, char *pFile2, char *pFile2s, float minLong, float maxLong, float minLat, float maxLat, float maxZenith, float *pNcepSST, float *pIce_DTime, float *pDiurnal ); float bilinearInterpolation( float x1, float val1, float x2, float val2, float y1, float val3, float y2, float val4, float xin, float yin ); int getTimeValue( int Year, int Month, int Day, float utc ); int getJulianDay( int mm, int id, int iyyy ); void getBiasMTLS( int satId, int dayTime, float MTLSQual, float pClr, float *pOutBias, float *pOutStdev ); void getBias( int satId, int dayTime, float tauVal, float tempDiff, float pClr, float *pBias, float *pStdev ); void getBiasOld( int satId, int dayTime, float tauVal, float tempDiff, float *pBias, float *pStdev ); void getBiasNew( int satId, float pClr, int dayTime, float dBTdSST, float tempDiff, float *pBias, float *pStdev ); void getBiasValue( int dayTime, float tauVal, float tempDiff, float *pTauArrayDay, float *pTauArrayNight, float *pConstDay, float *pSlopeDay, float *pConstNight, float *pSlopeNight, float *pStdevDay, float *pStdevNight, float *pBias, float *pStdev ); #ifdef GDS16 void read_in_sst_data( char *pFileStem, int *pSatId, int *pYear, int *pMonth, int *pDay, int *pHours, int *pMins, int *pSeconds, float *pUtc, int *pNele, int *pNlin, float **ppLat, float **ppLon, float **ppSatZen, float **ppSolZen, unsigned char **ppSST, float **ppPclr, float **ppBT1, float **ppBT2, float **ppDBTdSST, char *pSatName, char *pNorthSouth, float *pPixelStep, float *pLineTime, char *pSpatial_Resolution, float *pSpatial_lat_res, float *pSpatial_lon_res, char *pSensor ); #else void read_in_sst_data( char *pFileStem, int *pSatId, int *pYear, int *pMonth, int *pDay, int *pHours, int *pMins, int *pSeconds, float *pUtc, int *pNele, int *pNlin, float **ppLat, float **ppLon, float **ppSatZen, float **ppSolZen, signed short **ppSST, float **ppPclr, float **ppBT1, float **ppBT2, float **ppDBTdSST, float **ppMTLS_Qual, char *pSatName, char *pNorthSouth, float *pPixelStep, float *pLineTime, char *pSpatial_Resolution, float *pSpatial_lat_res, float *pSpatial_lon_res, char *pSensor ); #endif void getCellPosition( float newLongitude, float newLatitude, int nDims1, int nDims2, float_fp_kind minLongitude, float_fp_kind maxLongitude, float_fp_kind stepLongitude, float_fp_kind minLatitude, float_fp_kind maxLatitude, float_fp_kind stepLatitude, int *pOffset1, int *pOffset2, int *pOffset3, int *pOffset4, float *pMinLongitudePos, float *pMaxLongitudePos, float *pMinLatitudePos, float *pMaxLatitudePos ); #ifdef RUN_CRTM void runCRTM_Model( int satId, struct modelVals *pModel, float minLong, float maxLong, float minLat, float maxLat, float maxZenith ); void copyModel( struct modelVals *pModel1, struct modelVals *pModel2 ); void copyModel2( struct modelVals2 *pModel1, struct modelVals2 *pModel2 ); float_fp_kind getZenithAnglePos( float_fp_kind sateliteLong, float_fp_kind sateliteLat, float_fp_kind longitude, float_fp_kind latitude, float_fp_kind sateliteHeight ); void getSunAngle( int Year, int Month, int Day, float_fp_kind Hours, float_fp_kind longitude, float_fp_kind latitude, float_fp_kind *pSunAngle ); #endif /* Function call to pCRTM model (F90) */ /* Note this is not currently used in the current version of the code * It is here as a placeholder in case it is needed */ extern void f90_derivemodel_MP_runmodel( /* Coded instrument type */ int *pSatId, /* Number of atmos layers */ int *pNLayers, /* Number of filters (3) */ int *pNFilters, /* Channels to output */ int *pChannelArray, /* satelite zenith angle */ float_fp_kind *pSateliteAngle, /* solar zenith angle */ float_fp_kind *pSolarAngle, /* height array */ float_fp_kind *pHeightArray, /* pressure array */ float_fp_kind *pPressureArray, /* temperaeture array */ float_fp_kind *pTemperatureArray, /* Ozone array */ float_fp_kind *pOzoneArray, /* RH array */ float_fp_kind *pWaterArray, /* surface pressure */ float_fp_kind *pSurfacePressure, /* surface temperature */ float_fp_kind *pSurfaceTemperature, /* wind speed */ float_fp_kind *pWindSpeed, /* OUTPUT: BT(nFilters) */ float_fp_kind *pOutBt, /* Radiance(nFilters) */ float_fp_kind *pOutRad, /* Transmittence(nFilters) */ float_fp_kind *pOutTau, /* dBT/dSST */ float_fp_kind *p_dBT_dSST, /* Error code */ int *pErrorCode ); /* Main rroutine */ int main( int argc, char *argv[] ) { int i = 0; int j = 0; int nele = 0; int nlin = 0; #ifdef GDS16 unsigned char *pUnMaskedSST = NULL; #else signed short *pUnMaskedSST = NULL; #endif float *pFltSST = NULL; float *pCloudProb = NULL; float *pLatitude = NULL; float *pLongitude = NULL; float *pZenith = NULL; float *pSolarAngle = NULL; float *pBT1 = NULL; float *pBT2 = NULL; float *pDBTdSST = NULL; float *pMTLSQual = NULL; float *pDiurnalModel = NULL; float latVal = 0.; float longVal = 0.; int Year = 0; int Month = 0; int Day = 0; int Hours = 0; int Mins = 0; int Seconds = 0; float utc = 0; int time = 0; float windSpeed = 0.; float ice = 0.; int landCode = 0.; float tau4 = 0.; float ssi = 0.; float timeDiff = 0.; float aerOptDepth = 0.; float aerOptDepthTime = 0.; float temp2m = 0.; float ncep_sst = 0.; float Diurnal = 0.; char file1[132]; char file1s[132]; char file2[132]; char file2s[132]; int satId = 78; /* Output arrays */ signed short *pSSTArr = NULL; short *pSST_DTimeArr = NULL; signed char *pSSES_Bias = NULL; signed char *pSSES_StdDev = NULL; signed char *pSSIArr = NULL; signed char *pSSI_DTimeArr = NULL; signed char *pSSISrc = NULL; signed char *pWindSpeedArr = NULL; signed char *pWindSpeed_DTimeArr = NULL; signed char *pWindSpeedSrc = NULL; signed char *pAerosolArr = NULL; signed char *pAerosol_DTimeArr = NULL; signed char *pAerosolSrc = NULL; signed char *pADIArr = NULL; signed char *pADI_DTimeArr = NULL; signed char *pIceArr = NULL; signed char *pIceSrc = NULL; signed char *pZenithArr = NULL; signed char *pSolarZenithArr = NULL; #ifdef GDS16 signed char *pRejection = NULL; #else short *pRejection = NULL; #endif signed char *pConfidence = NULL; signed char *pProximity = NULL; signed char *pClearSkyProb = NULL; signed char *pDiurnalEstimate = NULL; signed char *pDT = NULL; float minLong = 0.; float maxLong = 0.; float minLong2 = 0.; float maxLong2 = 0.; float minLong3 = 0.; float maxLong3 = 0.; float minLat = 0.; float maxLat = 0.; float maxZenith = 0.; float tempDifference = 0.; float biasVal = 0.; float stdevVal = 0.; float pixelTime = 0.; float pixelStep = 0.; float lineTime = 0; float prevSST = 0.; char startTime[MAX_STRING_LENGTH]; char stopTime[MAX_STRING_LENGTH]; char startDate[MAX_STRING_LENGTH]; char stopDate[MAX_STRING_LENGTH]; char inFile[MAX_STRING_LENGTH]; char satelliteName[MAX_STRING_LENGTH]; char northSouth[MAX_STRING_LENGTH]; char spatial_resolution[MAX_STRING_LENGTH]; char sensor[MAX_STRING_LENGTH]; int newDay = 0; int newMonth = 0; int newYear = 0; int newHours = 0; int newMins = 0; int newSeconds = 0; int dayNo = 0; float newUtc = 0.; float temp_time = 0.; float tVal = 0.; float clearSkyProbability = 0.; int validpos = 0; int nmean = 0; float spatial_lat_res = 0.; float spatial_lon_res = 0.; float Ice_DTime = 0.; /* These variables are for a test */ if( 2 != argc ){ fprintf(stderr, "USAGE: makeGODAE-level2p FileStem\n"); exit(-1); } #ifdef GDS16 printf(" Creating NOAA Level 2P data (GDS version 1.6)\n"); #else printf(" Creating NOAA Level 2P data (GDS version 2.0)\n"); #endif /* First get environment variables */ if( NULL == getenv("L2P_DATA") ) printError("L2P","L2P_DATA environment variable not defined",NULL,NULL); strcpy(DATA_DIR,getenv("L2P_DATA")); strcat(DATA_DIR,"/"); if( NULL == getenv("HOME") ) printError("L2P","HOME environment variable not defined",NULL,NULL); strcpy(WGRIB_EXEC,getenv("HOME")); strcat(WGRIB_EXEC,"/mcidas/bin/wgrib.x"); strcpy(WGRIB2_EXEC,getenv("HOME")); strcat(WGRIB2_EXEC,"/mcidas/bin/wgrib2.x"); /* Read in data from SST processing */ #ifdef GDS16 read_in_sst_data( argv[1], &satId, &Year, &Month, &Day, &Hours, &Mins, &Seconds, &utc, &nele, &nlin, &pLatitude, &pLongitude, &pZenith, &pSolarAngle, &pUnMaskedSST, &pCloudProb, &pBT1, &pBT2, &pDBTdSST, satelliteName, northSouth, &pixelTime, &lineTime, spatial_resolution, &spatial_lat_res, &spatial_lon_res, sensor ); #else read_in_sst_data( argv[1], &satId, &Year, &Month, &Day, &Hours, &Mins, &Seconds, &utc, &nele, &nlin, &pLatitude, &pLongitude, &pZenith, &pSolarAngle, &pUnMaskedSST, &pCloudProb, &pBT1, &pBT2, &pDBTdSST, &pMTLSQual, satelliteName, northSouth, &pixelTime, &lineTime, spatial_resolution, &spatial_lat_res, &spatial_lon_res, sensor ); #endif /* Get min/max lat/longitude */ /* Note by definition this has to be from -180 to +180 in longitude */ /* but output longitude is 0,360 so have to switch it */ #if 0 /* This isn't general enough... */ minLong = 1e8; maxLong = -1e8; minLong2 = 1e8; maxLong2 = -1e8; minLat = 1e8; maxLat = -1e8; for(i=0;i 180 ) longVal = longVal - 360.; /* Get Min/Max */ if( longVal > 0 ){ minLong = myMIN(minLong,longVal); maxLong = myMAX(maxLong,longVal); } else { minLong2 = myMIN(minLong2,longVal); maxLong2 = myMAX(maxLong2,longVal); } minLat = myMIN(minLat,latVal); maxLat = myMAX(maxLat,latVal); } } } /* Take into account crossing date line */ if( minLong != 1e8 && minLong2 != 1e8 ){ if( minLong < 10. && maxLong2 > -10. ){ minLong3 = minLong2; maxLong3 = maxLong; } else { minLong3 = maxLong2; maxLong3 = minLong; } minLong = minLong3; maxLong = maxLong3; } #else /* Scan down left and right edges and top bottom to get square */ /* Find position of LHS which has valid longitude data in it */ validpos = -1; for(i=0;i=0;i--){ for(j=0;j=0;i--){ for(j=0;j maxLat ){ latVal = minLat; minLat = maxLat; maxLat = latVal; } if( minLong > 180. ) minLong -= 360; if( maxLong > 180. ) maxLong -= 360; #endif /* Get the time (long integer from 1 Jan 1980) from time of image */ time = getTimeValue( Year, Month, Day, utc ); /* Get start time in correct output format */ #ifdef GDS16 sprintf(startTime,"%2.2d:%2.2d:%2.2d UTC", Hours,Mins,Seconds); sprintf(startDate,"%4.4d-%2.2d-%2.2d", Year,Month,Day); #else sprintf(startTime,"%4.4d%2.2d%2.2dT%2.2d%2.2d%2.2dZ", Year,Month,Day,Hours,Mins,Seconds); #endif /* Get end time from the start time and number of lines times the time it * takes to get each time (100 lines / minute) */ /* New UTC in seconds */ newUtc = utc+nlin*lineTime; /* If gone over the day - get new day/month/year */ if( 86400 < newUtc ){ /* This routine is in read_grib.c */ getDayNoFromMonthDay( Year, Month, Day, &dayNo ); /* Have we crossed over to a new year? */ if( ((0 == Year%4 && 0 != Year%100) || 0 == Year%400) && 366 == dayNo){ /* leap year */ newYear = Year + 1; newMonth = 1; newDay = 1; } else if( 365 == dayNo ){ /* Not a leap year */ newYear = Year + 1; newMonth = 1; newDay = 1; } else { /* Just get new day/month/year */ dayNo += 1; getMonthFromDayNo(Year,dayNo,&newMonth,&newDay); newYear = Year; } /* Get correct utc for this day */ newUtc -= 86400; newHours = (int) (newUtc/3600.); temp_time = ((newUtc - newHours*3600.)/60.); newMins = (int) temp_time; newSeconds = (int) ((temp_time - (int)temp_time)*60.); } else { /* Not crossed over the day boundary */ newYear = Year; newMonth = Month; newDay = Day; newHours = (int) (newUtc/3600.); temp_time = ((newUtc - newHours*3600.)/60.); newMins = (int) temp_time; newSeconds = (int) ((temp_time - (int)temp_time)*60.); } /* Reformat end time */ #ifdef GDS16 sprintf(stopTime,"%2.2d:%2.2d:%2.2d UTC", newHours,newMins,newSeconds); sprintf(stopDate,"%4.4d-%2.2d-%2.2d", newYear,newMonth,newDay); #else sprintf(stopTime,"%4.4d%2.2d%2.2dT%2.2d%2.2d%2.2dZ", newYear,newMonth,newDay,newHours,newMins,newSeconds); #endif /* Now loop round positions and make arrays for wind speed etc */ /* First Allocate arrays of output */ /* THis is a work array for the proper SST (comes in the form of a byte * array */ if( NULL == (pFltSST = (float *) calloc( nele*nlin, sizeof(float) )) ) printError("L2P","Cannot allocate pSSTArr",NULL,NULL); /* Output SST value - as a scaled short */ if( NULL == (pSSTArr = (short *) calloc( nele*nlin, sizeof(short) )) ) printError("L2P","Cannot allocate pSSTArr",NULL,NULL); /* Time from start time of image of each SST value (seconds) */ if( NULL == (pSST_DTimeArr = (short *) calloc( nele*nlin, sizeof(short) )) ) printError("L2P","Cannot allocate pSST_DTimeArr",NULL,NULL); /* DT analysis - difference from now and previous days analysis */ if( NULL == (pDT = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pDT",NULL,NULL); /* SSES */ if( NULL == (pSSES_Bias = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pSSES_Bias",NULL,NULL); if( NULL == (pSSES_StdDev = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pSSES_StdDev",NULL,NULL); /* SSI */ if( NULL == (pSSIArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pSSIArr",NULL,NULL); /* Time difference from current time to estimated SSI time */ if( NULL == (pSSI_DTimeArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pSSI_DTimeArr",NULL,NULL); /* The source of the SSI (in our case the NCEP models */ if( NULL == (pSSISrc = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pSSIArr",NULL,NULL); /* Wind Speed */ if( NULL == (pWindSpeedArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pWindSpeedArr",NULL,NULL); if( NULL == (pWindSpeed_DTimeArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pWindSpeed_DTimeArr",NULL,NULL); /* SOurce of wind speed data (NCEP models) */ if( NULL == (pWindSpeedSrc = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pWindSpeedSrc",NULL,NULL); /* Aerosol data */ if( NULL == (pAerosolArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pWindSpeedArr",NULL,NULL); /* Time difference of aerosol observation to current time */ if( NULL == (pAerosol_DTimeArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pWindSpeed_DTimeArr",NULL,NULL); /* SOurce of aerosol data */ if( NULL == (pAerosolSrc = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pAerosolSrc",NULL,NULL); /* Aerosol dynamic indicator */ if( NULL == (pADIArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pADIArr",NULL,NULL); if( NULL == (pADI_DTimeArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pADI_DTime",NULL,NULL); /* Ice fraction */ if( NULL == (pIceArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pIceArr",NULL,NULL); /* Source of ice fraction data - in our case from NCEP */ if( NULL == (pIceSrc = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pIceSrc",NULL,NULL); /* Zenith Angle - from bob's code */ if( NULL == (pZenithArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pZenithArr",NULL,NULL); /* Solar Zenith Angle - from bob's code */ if( NULL == (pSolarZenithArr = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pZenithArr",NULL,NULL); /* Now a set of flags defined in the GDS */ /* Rejection flag - should a pixel be rejected and why */ /* For GDS 2.0 meaning has changed to just an l2p_flag */ #ifdef GDS16 if( NULL == (pRejection = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pRejection",NULL,NULL); /* Confidence flag */ if( NULL == (pConfidence = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pConfidence",NULL,NULL); #else if( NULL == (pRejection = (short *) calloc( nele*nlin, sizeof(short) )) ) printError("L2P","Cannot allocate pRejection",NULL,NULL); #endif /* Confidence flag */ /* Does not exist in GDS 2.0 */ /* Proximity confidence flag - values based on clear sky probability */ /* For GDS 2.0 changed to quality_level */ if( NULL == (pProximity = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pProximity",NULL,NULL); /* Experimental field - bayesian clear sky proabability */ if( NULL == (pClearSkyProb = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pClearSkyProb",NULL,NULL); /* Experimental field - diurnal warming estimate */ if( NULL == (pDiurnalEstimate = (signed char *) calloc( nele*nlin, sizeof(signed char) )) ) printError("L2P","Cannot allocate pDiurnalEstimate",NULL,NULL); /* This is so we can check the arrays inpendently as floats */ /* Loop round input arrays and get values for all the required parameters * The routine read_in_all_model_data will be linear interpolation in time * and space if necessary */ for(i=0;i *(pCloudProb+j+i*nele) && 0 == *(pRejection+j+i*nele) ) *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_CLOUDY; #ifndef GDS16 /* Set a day/night bits */ if( *(pSolarAngle+j+i*nele) > 90. ){ *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_NIGHT; } if( (*(pSolarAngle+j+i*nele) <= 95.) && (*(pSolarAngle+j+i*nele) > 80.) ){ *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_TWILIGHT; } if(*(pSolarAngle+j+i*nele) <= 80.){ *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_DAY; } if( SST_GOES_SCATTERING_FLAG == *(pUnMaskedSST+j+i*nele) ){ *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_GOES_SCAT; } #endif #ifdef GDS16 /* Check wind speed */ if( WIND_SPEED_MAX <= windSpeed ) *(pConfidence+j+i*nele) += CONFIDENCE_WIND; #endif /* If SST value valid do more checks */ if( 7 < *(pUnMaskedSST+j+i*nele) ){ /* Convert SST to K - to floating point number */ #ifdef GDS16 *(pFltSST+j+i*nele) = 270. + 0.15 * *(pUnMaskedSST+j+i*nele); /* printf("%d %d : %f\n",j,i,*(pFltSST+j+i*nele));*/ #else *(pFltSST+j+i*nele) = 271.35 + 0.01 * *(pUnMaskedSST+j+i*nele); #endif /* Check against previous SST */ tempDifference = *(pFltSST+j+i*nele) - prevSST; /* If out of range set bad */ if( MIN_ALLOWED_SST > *(pFltSST+j+i*nele) || MAX_ALLOWED_SST < *(pFltSST+j+i*nele) ){ /* SST outside allowed range - set rejection flags etc */ *(pSSTArr+j+i*nele) = SST_NODATA; *(pDT+j+i*nele) = DT_ANAL_NODATA; *(pSSES_Bias+j+i*nele) = SSES_BIAS_NODATA; *(pSSES_StdDev+j+i*nele) = SSES_STDEV_NODATA; *(pProximity+j+i*nele) = PROXIMITY_UNPROCESSED; *(pClearSkyProb+j+i*nele) = CLEAR_SKY_PROB_NODATA; *(pDiurnalEstimate+j+i*nele) = DIURNAL_EST_NODATA; } else { /* Convert SST/delta Time to netCDF format */ *(pSSTArr+j+i*nele) = (SST_TYPE) ((*(pFltSST+j+i*nele) - SST_OFFSET)/SST_FACTOR); /* *(pDT+j+i*nele) = (DT_ANAL_TYPE) ((tempDifference - DT_ANAL_OFFSET)/DT_ANAL_FACTOR); */ tVal = ((tempDifference - DT_ANAL_OFFSET)/DT_ANAL_FACTOR); if( tVal < DT_ANAL_MIN_VAL ){ tVal = DT_ANAL_MIN_VAL; } else if( tVal > DT_ANAL_MAX_VAL ){ tVal = DT_ANAL_MAX_VAL; } *(pDT+j+i*nele) = (DT_ANAL_TYPE) tVal; /* Get the biases from the opt depths */ /* First need temperature difference between SST and 2m temp, both * from the NCEP models */ tempDifference = temp2m - ncep_sst; /* Use lookup tables to return bias/stdev */ clearSkyProbability = *(pCloudProb+j+i*nele); /* Use solar angle to determing day/night */ #ifdef GDS16 if( 90 > *(pSolarAngle+j+i*nele) ) getBias( satId, 1, *(pDBTdSST+j+i*nele), tempDifference, clearSkyProbability, &biasVal, &stdevVal ); else getBias( satId, 0, *(pDBTdSST+j+i*nele), tempDifference, clearSkyProbability, &biasVal, &stdevVal ); #else if( 90 > *(pSolarAngle+j+i*nele) ) getBiasMTLS( satId, 1, *(pMTLSQual+j+i*nele), clearSkyProbability, &biasVal, &stdevVal ); else getBiasMTLS( satId, 0, *(pMTLSQual+j+i*nele), clearSkyProbability, &biasVal, &stdevVal ); #endif /* If OK bias value - put in netCDF format */ if( -998 < biasVal && -998 < stdevVal ){ /* Make sure bias is clamped (if needed) to range -2.54,2.54 */ if( SSES_BIAS_MIN_VAL*SSES_BIAS_FACTOR > biasVal-SSES_BIAS_OFFSET ){ biasVal = SSES_BIAS_MIN_VAL*SSES_BIAS_FACTOR; } if( SSES_BIAS_MAX_VAL*SSES_BIAS_FACTOR < biasVal-SSES_BIAS_OFFSET ){ biasVal = SSES_BIAS_MAX_VAL*SSES_BIAS_FACTOR; } *(pSSES_Bias+j+i*nele) = (SSES_BIAS_TYPE) ((biasVal - SSES_BIAS_OFFSET)/SSES_BIAS_FACTOR); *(pSSES_StdDev+j+i*nele) = (SSES_STDEV_TYPE) ((stdevVal - SSES_STDEV_OFFSET)/SSES_STDEV_FACTOR); } else { *(pSSES_Bias+j+i*nele) = SSES_BIAS_NODATA; *(pSSES_StdDev+j+i*nele) = SSES_STDEV_NODATA; } /* Set proximity flag confidence - this is based on clear * sky probability */ /* Proximity confidence as defined */ /* * 0 = Unprocessed * 1 = Cloudy ( < 0.8 probability ) * 2 = Bad ( 0.8 <= probability < 0.95 ) * 3 = Suspect ( 0.95 <= probability < 0.99 ) * 4 = Acceptable ( 0.99 <= probability < 0.999 ) * 5 = Excellent ( 0.999 <= probability ) */ /* Decision tree as to which level we are at */ /* Note min/max values set in writeNetCDF.h */ if( MIN_PROXIMITY_CODE_1 <= clearSkyProbability && MAX_PROXIMITY_CODE_1 > clearSkyProbability ) *(pProximity+j+i*nele) = PROXIMITY_CLOUDY; else if( MIN_PROXIMITY_CODE_2 <= clearSkyProbability && MAX_PROXIMITY_CODE_2 > clearSkyProbability ) *(pProximity+j+i*nele) = PROXIMITY_BAD; else if( MIN_PROXIMITY_CODE_3 <= clearSkyProbability && MAX_PROXIMITY_CODE_3 > clearSkyProbability ) *(pProximity+j+i*nele) = PROXIMITY_SUSPECT; else if( MIN_PROXIMITY_CODE_4 <= clearSkyProbability && MAX_PROXIMITY_CODE_4 > clearSkyProbability ) *(pProximity+j+i*nele) = PROXIMITY_ACCEPTABLE; else if( MIN_PROXIMITY_CODE_5 <= clearSkyProbability && MAX_PROXIMITY_CODE_5 > clearSkyProbability ) *(pProximity+j+i*nele) = PROXIMITY_EXCELLENT; /* Clear sky probability - Experimental field */ *(pClearSkyProb+j+i*nele) = (CLEAR_SKY_PROB_TYPE) ((*(pCloudProb+j+i*nele) - CLEAR_SKY_PROB_OFFSET)/CLEAR_SKY_PROB_FACTOR); } /* Now add rest of data - only needed if there is SST data */ /* SSI checks - no present in GDS2.0 (optional field) */ #ifdef GDS16 if( -999 != ssi ){ /* If SSI data is OK then check that its within range */ if( ssi > 1000. ) ssi = 1000.; if( ssi < 0. ) ssi = 0.; /* Convert to single byte form */ *(pSSIArr+j+i*nele) = (SSI_TYPE) ((ssi - SSI_OFFSET)/ SSI_FACTOR+0.5); *(pSSI_DTimeArr+j+i*nele) = (SSI_TYPE) (((timeDiff-SSI_DTIME_OFFSET)/SSI_DTIME_FACTOR)+0.5); /* At this point still need discussion about some new codes to * include possible land contamination flags - but not yet in * the GDS */ #ifdef NOT_GDS_COMPLIENT_NOT_USE if( 1 == landCode ) *(pSSISrc+j+i*nele) = SSI_LAND_CODE; else *(pSSISrc+j+i*nele) = SSI_CODE; #else *(pSSISrc+j+i*nele) = SSI_CODE; #endif } else { /* Fill value with invalid data value */ *(pSSIArr+j+i*nele) = SSI_NODATA; *(pSSI_DTimeArr+j+i*nele) = SSI_DTIME_NODATA; #ifdef NOT_GDS_COMPLIENT_NOT_USE if( 1 == landCode ) *(pSSISrc+j+i*nele) = SSI_LAND_CODE; else *(pSSISrc+j+i*nele) = SSI_CODE; #else *(pSSISrc+j+i*nele) = SSI_CODE; #endif } #endif /* Wind speed */ /* Check to see if wind speed there */ if( -999 != windSpeed ){ /* Convert to single byte value */ *(pWindSpeedArr+j+i*nele) = (WIND_TYPE) (((windSpeed-WIND_OFFSET)/WIND_FACTOR)+0.5); /* Note wind speed time is 0 since interpolated in time */ *(pWindSpeed_DTimeArr+j+i*nele) = (WIND_DTIME_TYPE) 0.; if( 1 == landCode ) *(pWindSpeedSrc+j+i*nele) = WIND_LAND_CODE; else *(pWindSpeedSrc+j+i*nele) = WIND_CODE; } else { *(pWindSpeedArr+j+i*nele) = WIND_NODATA; *(pWindSpeed_DTimeArr+j+i*nele) = WIND_DTIME_NODATA; if( 1 == landCode ) *(pWindSpeedSrc+j+i*nele) = WIND_LAND_CODE; else *(pWindSpeedSrc+j+i*nele) = WIND_CODE; } /* Output zenith angle */ /* Check to see if zenith angle there */ if( -99 != *(pZenith+j+i*nele) ){ /* Convert to single byte getting NINT value */ if( 0 < *(pZenith+j+i*nele) ) *(pZenithArr+j+i*nele) = (ZENITH_TYPE) (((*(pZenith+j+i*nele)-ZENITH_OFFSET)/ZENITH_FACTOR)+0.5); else *(pZenithArr+j+i*nele) = (ZENITH_TYPE) (((*(pZenith+j+i*nele)-ZENITH_OFFSET)/ZENITH_FACTOR)-0.5); } else { /* No zenith angle data there */ *(pZenithArr+j+i*nele) = ZENITH_NODATA; } /* Check to see if solar zenith angle there */ if( -99 != *(pSolarAngle+j+i*nele) ){ /* Convert to single byte getting NINT value */ if( 0 < *(pSolarAngle+j+i*nele) ) *(pSolarZenithArr+j+i*nele) = (SOLAR_ZENITH_TYPE) (((*(pSolarAngle+j+i*nele)-SOLAR_ZENITH_OFFSET)/ SOLAR_ZENITH_FACTOR)+0.5); else *(pSolarZenithArr+j+i*nele) = (SOLAR_ZENITH_TYPE) (((*(pSolarAngle+j+i*nele)-SOLAR_ZENITH_OFFSET)/ SOLAR_ZENITH_FACTOR)-0.5); } else { /* No zenith angle data there */ *(pSolarZenithArr+j+i*nele) = SOLAR_ZENITH_NODATA; } /* Check to see if ice fraction data there */ if( -999 != ice ){ /* Convert to byte value - this may change */ *(pIceArr+j+i*nele) = (ICE_TYPE) (ice-ICE_OFFSET)/ICE_FACTOR; /* If Ice then set ICE flag */ if( 0. < ice ) *(pRejection+j+i*nele) = *(pRejection+j+i*nele) | REJECTION_ICE; } else { /* No ice data available */ *(pIceArr+j+i*nele) = ICE_NODATA; } *(pIceSrc+j+i*nele) = ICE_CODE; /* Aerosol */ /* Check to see if data there */ if( -999 != aerOptDepth ){ /* Set to max limit */ if( 2.54 <= aerOptDepth) aerOptDepth = 2.54; /* Convert to single byte value */ *(pAerosolArr+j+i*nele) = (AEROSOL_TYPE) (((aerOptDepth-AEROSOL_OFFSET)/AEROSOL_FACTOR)+0.5); /* Limit time to maximum allowed by GDS */ #if 0 if( 25 <= aerOptDepthTime || 0 == aerOptDepthTime ) aerOptDepthTime = 25; #else if( 254 <= aerOptDepthTime || -1 == aerOptDepthTime){ /* No data available */ *(pAerosolArr+j+i*nele) = AEROSOL_NODATA; *(pAerosol_DTimeArr+j+i*nele) = AEROSOL_DTIME_NODATA; *(pAerosolSrc+j+i*nele) = (SRC_AEROSOL_TYPE) AEROSOL_CODE; } else { /* Convert time to byte value */ *(pAerosol_DTimeArr+j+i*nele) = (AEROSOL_DTIME_TYPE) (((aerOptDepthTime-AEROSOL_DTIME_OFFSET)/ AEROSOL_DTIME_FACTOR)+0.5); *(pAerosolSrc+j+i*nele) = (SRC_AEROSOL_TYPE) AEROSOL_CODE; } #endif } else { /* No data available */ *(pAerosolArr+j+i*nele) = AEROSOL_NODATA; *(pAerosol_DTimeArr+j+i*nele) = AEROSOL_DTIME_NODATA; *(pAerosolSrc+j+i*nele) = (SRC_AEROSOL_TYPE) AEROSOL_CODE; } /* Aerosol Dynamic Indicator - needs info from new SST */ /* For moment use Aerosol measurements */ *(pADIArr+j+i*nele) = *(pAerosolArr+j+i*nele); *(pADI_DTimeArr+j+i*nele) = *(pAerosol_DTimeArr+j+i*nele); /* Diurnal warming estimate - Experimental field */ if( -999 != Diurnal ){ *(pDiurnalEstimate+j+i*nele) = (DIURNAL_EST_TYPE) ((Diurnal - DIURNAL_EST_OFFSET)/DIURNAL_EST_FACTOR); } else { *(pDiurnalEstimate+j+i*nele) = DIURNAL_EST_NODATA; } } else { /* No SST data - set flags/data appropriately */ *(pFltSST+j+i*nele) = 0.; #ifdef GDS16 *(pRejection+j+i*nele) += REJECTION_SST_RANGE; #endif *(pSSTArr+j+i*nele) = SST_NODATA; *(pDT+j+i*nele) = DT_ANAL_NODATA; *(pSSES_Bias+j+i*nele) = SSES_BIAS_NODATA; *(pSSES_StdDev+j+i*nele) = SSES_STDEV_NODATA; *(pProximity+j+i*nele) = PROXIMITY_UNPROCESSED; *(pClearSkyProb+j+i*nele) = CLEAR_SKY_PROB_NODATA; *(pDiurnalEstimate+j+i*nele) = DIURNAL_EST_NODATA; /* Rest of data set to no data */ #ifdef GDS16 *(pSST_DTimeArr+j+i*nele) = SST_DTIME_NODATA; *(pSSIArr+j+i*nele) = SSI_NODATA; *(pSSI_DTimeArr+j+i*nele) = SSI_DTIME_NODATA; *(pWindSpeedArr+j+i*nele) = WIND_NODATA; *(pWindSpeed_DTimeArr+j+i*nele) = WIND_DTIME_NODATA; *(pAerosolArr+j+i*nele) = AEROSOL_NODATA; *(pAerosol_DTimeArr+j+i*nele) = AEROSOL_DTIME_NODATA; *(pAerosolSrc+j+i*nele) = AEROSOL_DTIME_NODATA; *(pIceArr+j+i*nele) = ICE_NODATA; *(pIceSrc+j+i*nele) = SRC_ICE_NODATA; *(pSSISrc+j+i*nele) = SRC_SSI_NODATA; *(pWindSpeedSrc+j+i*nele) = SRC_WIND_NODATA; *(pZenithArr+j+i*nele) = ZENITH_NODATA; #else *(pSST_DTimeArr+j+i*nele) = SST_DTIME_NODATA; *(pWindSpeedArr+j+i*nele) = WIND_NODATA; *(pWindSpeed_DTimeArr+j+i*nele) = WIND_DTIME_NODATA; *(pAerosolArr+j+i*nele) = AEROSOL_NODATA; *(pAerosol_DTimeArr+j+i*nele) = AEROSOL_DTIME_NODATA; *(pAerosolSrc+j+i*nele) = AEROSOL_DTIME_NODATA; *(pIceArr+j+i*nele) = ICE_NODATA; *(pIceSrc+j+i*nele) = SRC_ICE_NODATA; *(pWindSpeedSrc+j+i*nele) = SRC_WIND_NODATA; *(pZenithArr+j+i*nele) = ZENITH_NODATA; *(pSolarZenithArr+j+i*nele) = SOLAR_ZENITH_NODATA; #endif } } else { /* No data at all - set flags */ #ifdef GDS16 *(pRejection+j+i*nele) = REJECTION_SST_RANGE; *(pConfidence+j+i*nele) = 0; #endif *(pProximity+j+i*nele) = PROXIMITY_UNPROCESSED; *(pSSTArr+j+i*nele) = SST_NODATA; *(pSST_DTimeArr+j+i*nele) = SST_DTIME_NODATA; *(pSSES_Bias+j+i*nele) = SSES_BIAS_NODATA; *(pSSES_StdDev+j+i*nele) = SSES_STDEV_NODATA; *(pDT+j+i*nele) = DT_ANAL_NODATA; *(pWindSpeedArr+j+i*nele) = WIND_NODATA; *(pWindSpeed_DTimeArr+j+i*nele) = WIND_DTIME_NODATA; *(pAerosolArr+j+i*nele) = AEROSOL_NODATA; *(pAerosol_DTimeArr+j+i*nele) = AEROSOL_DTIME_NODATA; *(pAerosolSrc+j+i*nele) = AEROSOL_DTIME_NODATA; *(pIceArr+j+i*nele) = ICE_NODATA; *(pIceSrc+j+i*nele) = SRC_ICE_NODATA; *(pWindSpeedSrc+j+i*nele) = SRC_WIND_NODATA; *(pZenithArr+j+i*nele) = ZENITH_NODATA; *(pSolarZenithArr+j+i*nele) = SOLAR_ZENITH_NODATA; *(pClearSkyProb+j+i*nele) = CLEAR_SKY_PROB_NODATA; *(pDiurnalEstimate+j+i*nele) = DIURNAL_EST_NODATA; prevSST = 265; *(pFltSST+j+i*nele) = 0.; } } } /* write data to netCDF file */ writeNetCDF( Year,Month,Day,Hours,Mins,Seconds, satelliteName, northSouth, 2, inFile, startTime, stopTime, startDate, stopDate, minLat, maxLat, minLong, maxLong, nele, nlin, time, pLongitude, pLatitude, pSSTArr, pSST_DTimeArr, pSSES_Bias, pSSES_StdDev, pSSIArr, pSSI_DTimeArr, pSSISrc, pDT, pWindSpeedArr, pWindSpeed_DTimeArr, pWindSpeedSrc, pAerosolArr, pAerosol_DTimeArr, pAerosolSrc, pIceArr, pIceSrc, pZenithArr, pRejection, pConfidence, pProximity, pClearSkyProb, Year, Month, Day, Hours, Mins, Seconds, newYear, newMonth, newDay, newHours, newMins, newSeconds, spatial_resolution, spatial_lat_res, spatial_lon_res, sensor, &Ice_DTime, pADIArr, pADI_DTimeArr,pSolarZenithArr,pDiurnalEstimate); /* Free memory */ free(pSSTArr); free(pSST_DTimeArr); free(pSSES_Bias); free(pSSES_StdDev); free(pDT); free(pSSIArr); free(pSSI_DTimeArr); free(pWindSpeedArr); free(pWindSpeed_DTimeArr); free(pAerosolArr); free(pAerosol_DTimeArr); free(pADIArr); free(pADI_DTimeArr); free(pIceArr); free(pIceSrc); free(pZenithArr); free(pSolarZenithArr); free(pWindSpeedSrc); free(pAerosolSrc); free(pSSISrc); free(pRejection); #ifdef GDS16 free(pConfidence); #endif free(pProximity); free(pClearSkyProb); free(pDiurnalEstimate); free(pUnMaskedSST); free(pFltSST); free(pCloudProb); free(pLatitude); free(pLongitude); free(pZenith); free(pSolarAngle); free(pBT1); free(pBT2); free(pDBTdSST); free(pDiurnal1); free(pDiurnal2); /* End of program */ return 0; } /* Routine to read in data from Bob's output */ /* Read in SST and other data */ #define MCIDAS_HEADER_SIZE 64 #ifdef GDS16 void read_in_sst_data( char *pFileStem, int *pSatId, int *pYear, int *pMonth, int *pDay, int *pHours, int *pMins, int *pSeconds, float *pUtc, int *pNele, int *pNlin, float **ppLat, float **ppLon, float **ppSatZen, float **ppSolZen, unsigned char **ppSST, float **ppPclr, float **ppBT1, float **ppBT2, float **ppDBTdSST, char *pSatName, char *pNorthSouth, float *pPixelStep, float *pLineTime, char *pSpatial_Resolution, float *pSpatial_lat_res, float *pSpatial_lon_res, char *pSensor ) #else void read_in_sst_data( char *pFileStem, int *pSatId, int *pYear, int *pMonth, int *pDay, int *pHours, int *pMins, int *pSeconds, float *pUtc, int *pNele, int *pNlin, float **ppLat, float **ppLon, float **ppSatZen, float **ppSolZen, signed short **ppSST, float **ppPclr, float **ppBT1, float **ppBT2, float **ppDBTdSST, float **ppMTLS_Qual, char *pSatName, char *pNorthSouth, float *pPixelStep, float *pLineTime, char *pSpatial_Resolution, float *pSpatial_lat_res, float *pSpatial_lon_res, char *pSensor ) #endif { int mcIdasHeader[MCIDAS_HEADER_SIZE]; int offset = 0; int dayNo = 0; int i = 0; int Hour = 0; int Min = 0; float temp = 0.; int Seconds = 0; char filename[MAX_STRING_LENGTH]; char errorString[MAX_STRING_LENGTH]; FILE *fp = NULL; int idNum = 0; float value = 0.; int nread = 0; int useGoes = 0; #ifdef GDS16 unsigned char *pPclr = NULL; unsigned char *pBT1 = NULL; unsigned char *pBT2 = NULL; unsigned char *pDBTdSST = NULL; #endif float pixelStep_GOES = 0.47; float lineTime_GOES = 0.47; float pixelStep_MTSAT = 0.6; float lineTime_MTSAT = 0.6; float pixelStep_MSG = 0.2; float lineTime_MSG = 0.2; /* Get filename from data input directory */ /* Read in header file with dimensions/dates etc. */ strcpy(filename,DATA_DIR); strcat(filename,"/netcdf_input/"); strcat(filename,pFileStem); strcat(filename,".header"); if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","File not found : %s",NULL,filename); if( MCIDAS_HEADER_SIZE != fread(mcIdasHeader,sizeof(int),MCIDAS_HEADER_SIZE,fp) ) printError("L2P","File reading error : %s",NULL,filename); fclose(fp); /* Get information from McIdas header */ /* Satelite ID */ *pSatId = mcIdasHeader[2]; /* Date */ *pYear = (int) (mcIdasHeader[3]/1000); dayNo = mcIdasHeader[3] - *pYear*1000; *pYear = 1900 + *pYear; getMonthFromDayNo(*pYear,dayNo,pMonth,pDay); /* Time */ Hour = (int) (mcIdasHeader[4]/10000); temp = (mcIdasHeader[4] - Hour * 10000.); Min = (int) (temp/100); Seconds = (int) (temp - Min*100); *pUtc = Hour*3600. + Min*60. + Seconds; *pHours = Hour; *pMins = Min; *pSeconds = Seconds; /* Size of data array */ *pNele = mcIdasHeader[9]; *pNlin = mcIdasHeader[8]; /* From satelite ID get which satelite we're using */ switch( *pSatId ){ case 72: strcpy(pSatName,"GOES9"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 74: strcpy(pSatName,"GOES10"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 76: strcpy(pSatName,"GOES11"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 78: strcpy(pSatName,"GOES12"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 180: strcpy(pSatName,"GOES13"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 184: strcpy(pSatName,"GOES15"); useGoes = 1; *pPixelStep = pixelStep_GOES; *pLineTime = lineTime_GOES; strcpy(pSpatial_Resolution,"2.5 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.023; strcpy(pSensor,"Imager"); break; case 84: strcpy(pSatName,"MTSAT1R"); *pPixelStep = pixelStep_MTSAT; *pLineTime = lineTime_MTSAT; strcpy(pSpatial_Resolution,"4 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.036; strcpy(pSensor,"Imager"); break; case 85: strcpy(pSatName,"MTSAT2"); *pPixelStep = pixelStep_MTSAT; *pLineTime = lineTime_MTSAT; strcpy(pSpatial_Resolution,"4 x 4 km at nadir"); *pSpatial_lat_res = 0.036; *pSpatial_lon_res = 0.036; strcpy(pSensor,"Imager"); break; case 52: strcpy(pSatName,"MSG02"); *pPixelStep = pixelStep_MSG; *pLineTime = lineTime_MSG; strcpy(pSpatial_Resolution,"3 x 3 km at nadir"); *pSpatial_lat_res = 0.027; *pSpatial_lon_res = 0.027; strcpy(pSensor,"SEVIRI"); break; case 53: strcpy(pSatName,"MSG03"); *pPixelStep = pixelStep_MSG; *pLineTime = lineTime_MSG; strcpy(pSpatial_Resolution,"3 x 3 km at nadir"); *pSpatial_lat_res = 0.027; *pSpatial_lon_res = 0.027; strcpy(pSensor,"SEVIRI"); break; default: printf("satID = %d\n",*pSatId); printError("L2P","Satellite ID not allowed",NULL,NULL); } if( 1 == useGoes ){ /* Strip out number from file stem */ /* This number telss us if we're north/south */ sscanf((pFileStem+4),"%4d",&idNum); switch(idNum){ case 142: strcpy(pNorthSouth,"North"); break; case 242: strcpy(pNorthSouth,"North"); break; case 342: strcpy(pNorthSouth,"South"); break; case 442: strcpy(pNorthSouth,"South"); break; case 542: strcpy(pNorthSouth,"North"); break; case 642: strcpy(pNorthSouth,"North"); break; case 742: strcpy(pNorthSouth,"South"); break; case 842: strcpy(pNorthSouth,"South"); break; } } else strcpy(pNorthSouth," "); /* Some informational stuff written to the log file */ if( 1 == useGoes ){ printf("L2P: %s %s dataset\n",pSatName,pNorthSouth); } else { printf("L2P: %s dataset\n",pSatName); } printf("L2P: File Date : %4.4d-%2.2d-%2.2d:%2.2d:%2.2d\n", *pYear,*pMonth,*pDay,Hour,Min); printf("L2P: File Size : %d %d\n",*pNele,*pNlin); /* Open Nav file - contains long/lat/zenith/solar angles */ strcpy(filename,DATA_DIR); strcat(filename,"/netcdf_input/"); strcat(filename,pFileStem); strcat(filename,".nav"); if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","File not found : %s",NULL,filename); printf("Image size : %d %d\n",*pNele,*pNlin); /* Allocate output arrays */ if( NULL == (*ppLat = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppLat",NULL,NULL); if( NULL == (*ppLon = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppLon",NULL,NULL); if( NULL == (*ppSatZen = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppSatZen",NULL,NULL); if( NULL == (*ppSolZen = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppSatZen",NULL,NULL); /* Read in file */ #ifdef GDS16 for(i=0;i<*pNlin;i++){ offset = (i * *pNele); if( *pNele != fread((*ppLat+offset),sizeof(float),*pNele,fp) ) printError("L2P","Reading *ppLat",NULL,NULL); if( *pNele != fread((*ppLon+offset),sizeof(float),*pNele,fp) ) printError("L2P","Reading *ppLon",NULL,NULL); if( *pNele != fread((*ppSatZen+offset),sizeof(float),*pNele,fp) ) printError("L2P","Reading *ppSatZen",NULL,NULL); if( *pNele != fread((*ppSolZen+offset),sizeof(float),*pNele,fp) ) printError("L2P","Reading *ppSolZen",NULL,NULL); } #else if( *pNele * *pNlin != fread(*ppLat,sizeof(float),*pNele * *pNlin,fp) ) printError("L2P","Reading *ppLat",NULL,NULL); if( *pNele * *pNlin != fread(*ppLon,sizeof(float),*pNele * *pNlin,fp) ) printError("L2P","Reading *ppLon",NULL,NULL); if( *pNele * *pNlin != fread(*ppSatZen,sizeof(float),*pNele * *pNlin,fp) ) printError("L2P","Reading *ppSatZen",NULL,NULL); if( *pNele * *pNlin != fread(*ppSolZen,sizeof(float),*pNele * *pNlin,fp) ) printError("L2P","Reading *ppSolZen",NULL,NULL); #endif fclose(fp); /* Open Data file - contains the rest of the data */ strcpy(filename,DATA_DIR); strcat(filename,"/netcdf_input/"); strcat(filename,pFileStem); strcat(filename,".data"); if( NULL == (fp = fopen(filename,"r")) ) printError("L2P","File not found : %s",NULL,filename); #ifdef GDS16 /* Arrays for reading in */ if( NULL == (*ppSST = (unsigned char *) calloc( *pNele * *pNlin, sizeof(unsigned char))) ) printError("L2P","File allocating ppSST",NULL,NULL); if( NULL == (pPclr = (unsigned char *) calloc( *pNele * *pNlin, sizeof(unsigned char))) ) printError("L2P","File allocating pPclr",NULL,NULL); if( NULL == (pBT1 = (unsigned char *) calloc( *pNele * *pNlin, sizeof(unsigned char))) ) printError("L2P","File allocating pBT1",NULL,NULL); if( NULL == (pBT2 = (unsigned char *) calloc( *pNele * *pNlin, sizeof(unsigned char))) ) printError("L2P","File allocating pBT1",NULL,NULL); if( NULL == (pDBTdSST = (unsigned char *) calloc( *pNele * *pNlin, sizeof(unsigned char))) ) printError("L2P","File allocating pDBTdSST",NULL,NULL); #else /* Allocate output arrays - keep SST with flags in */ if( NULL == (*ppSST = (signed short *) calloc( *pNele * *pNlin, sizeof(signed short))) ) printError("L2P","File allocating ppSST",NULL,NULL); #endif if( NULL == (*ppPclr = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppPclr",NULL,NULL); if( NULL == (*ppBT1 = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppBT1",NULL,NULL); if( NULL == (*ppBT2 = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppBT1",NULL,NULL); if( NULL == (*ppDBTdSST = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppDBTdSST",NULL,NULL); #ifndef GDS16 if( NULL == (*ppMTLS_Qual = (float *) calloc( *pNele * *pNlin, sizeof(float))) ) printError("L2P","File allocating ppDBTdSST",NULL,NULL); #endif /* Arrays for reading in */ /* Read in file */ #ifdef GDS16 nread = fread(*ppSST,sizeof(unsigned char), *pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString,"Reading SST from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(pPclr,sizeof(unsigned char),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading prob. clear sky from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(pBT1,sizeof(unsigned char),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading BT1 from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(pBT2,sizeof(unsigned char),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading BT2 from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(pDBTdSST,sizeof(unsigned char), *pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading dBT/dSST from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } #else nread = fread(*ppSST,sizeof(unsigned short), *pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString,"Reading SST from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(*ppPclr,sizeof(float),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading prob. clear sky from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(*ppBT1,sizeof(float),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading BT1 from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(*ppBT2,sizeof(float),*pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading BT2 from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(*ppDBTdSST,sizeof(float), *pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading dBT/dSST from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } nread = fread(*ppMTLS_Qual,sizeof(float), *pNele * *pNlin,fp); if( *pNele * *pNlin != nread ){ sprintf(errorString, "Reading MTLS_Qual from data file (read in %d from %d)", nread, *pNele * *pNlin); printError("L2P",errorString,NULL,NULL); } #endif fclose(fp); /* Convert to floating point */ /* Keep SST as char due to flags in data */ for(i=0;i<*pNele * *pNlin;i++){ if( 0 != *(*ppSST+i) ){ #ifdef GDS16 /* Convert clear sky to real probability */ if( 2 <= *(pPclr+i) ) *(*ppPclr+i) = (1.002 - exp((2. - *(pPclr+i))/40.546121)); else *(*ppPclr+i) = 0.; /* Convert dBT/dSST from byte version to real value */ *(*ppDBTdSST+i) = ((*(pDBTdSST+i)-7.)/247./1.25)+0.1; /* Convert satelite zenith angle to real value */ #ifdef OLD_SST_CODE if( -99 != *(*ppSatZen+i) ){ value = 1./((fabs(*(*ppSatZen+i))/1000.)+1.); *(*ppSatZen+i) = acos(value)*180./M_PI; } #endif #endif /* Convert longitude from -180,180 to 0,360 for consistency with * NCEP models */ /* Note that data is stored degrees west but we want degrees east */ if( 0 > *(*ppLon+i) ){ if( -200 > *(*ppLon+i) ){ fprintf(stderr,"SST = %d\n",*(*ppSST+i)); fprintf(stderr,"ERROR: Longitude = -999\n"); exit(-1); } *(*ppLon+i) = fabs(*(*ppLon+i)); } else *(*ppLon+i) = 360. - *(*ppLon+i); } else { /* No good data - set to FLT_MAX as a flag of this */ *(*ppLat+i) = -FLT_MAX; *(*ppLon+i) = -FLT_MAX; *(*ppSatZen+i) = -FLT_MAX; *(*ppSolZen+i) = -FLT_MAX; *(*ppPclr+i) = -FLT_MAX; } } #ifdef GSD16 /* Free memory */ free(pPclr); free(pBT1); free(pBT2); free(pDBTdSST); #endif } /**************************************************************************** * * The next set of routines impliments Prabhat's method for determining the * bias and standard deviation based on MTLS Quality flag value */ #define GOES13_NDAY 4 float GOES13_DAY_MEAN_ERR[GOES13_NDAY] = {0.12333246,0.30582483,0.49816777,0.81220942}; float GOES13_DAY_BIAS[GOES13_NDAY] = {-0.16714452,-0.17823462,-0.31524014,-0.17634621}; float GOES13_DAY_STDEV[GOES13_NDAY] = {0.40610877,0.53710604,0.70698994,1.0536650}; #define GOES13_NNIGHT 4 float GOES13_NIGHT_MEAN_ERR[GOES13_NNIGHT] = {0.11938851,0.30421204,0.50042965,0.64737414}; float GOES13_NIGHT_BIAS[GOES13_NNIGHT] = {-0.10529313,-0.19734071,-0.50719410,-0.72369719}; float GOES13_NIGHT_STDEV[GOES13_NNIGHT] = {0.39621937,0.44018725,0.67060512,0.73292112}; #define GOES15_NDAY 4 float GOES15_DAY_MEAN_ERR[GOES15_NDAY] = {0.12833993,0.30482618,0.49309807,0.95089391}; float GOES15_DAY_BIAS[GOES15_NDAY] = {-0.17521824,-0.10895059,-0.072773516,-0.11720801}; float GOES15_DAY_STDEV[GOES15_NDAY] = {0.45205945,0.55458385,0.79229355,0.94084549}; #define GOES15_NNIGHT 4 float GOES15_NIGHT_MEAN_ERR[GOES15_NNIGHT] = {0.12080848,0.3036231,0.49403447,1.1285480}; float GOES15_NIGHT_BIAS[GOES15_NNIGHT] = {-0.10380731,-0.16360891,-0.37223655,-0.72134912}; float GOES15_NIGHT_STDEV[GOES15_NNIGHT] = {0.39552018,0.55945760,0.83352172,1.0038176}; #define MTSAT1_NDAY 4 float MTSAT1_DAY_MEAN_ERR[MTSAT1_NDAY] = {0.18517378,0.46828303,0.72531721,1.7187938}; float MTSAT1_DAY_BIAS[MTSAT1_NDAY] = {-0.16046284,-0.16493239,-0.19177689,-0.43815386}; float MTSAT1_DAY_STDEV[MTSAT1_NDAY] = {0.64536619,0.69978702,0.73633587,0.92903745}; #define MTSAT2_NDAY 4 float MTSAT2_DAY_MEAN_ERR[MTSAT2_NDAY] = {0.18517378,0.46828303,0.72531721,1.7187938}; float MTSAT2_DAY_BIAS[MTSAT2_NDAY] = {-0.16046284,-0.16493239,-0.19177689,-0.43815386}; float MTSAT2_DAY_STDEV[MTSAT2_NDAY] = {0.64536619,0.69978702,0.73633587,0.92903745}; #define MTSAT1_NNIGHT 4 float MTSAT1_NIGHT_MEAN_ERR[MTSAT1_NNIGHT] = {0.12043308,0.30219896,0.49609574,1.0430975}; float MTSAT1_NIGHT_BIAS[MTSAT1_NNIGHT] = {-0.059461016,-0.043822456,-0.10529868,-0.40427330}; float MTSAT1_NIGHT_STDEV[MTSAT1_NNIGHT] = {0.46123141,0.47195938,0.58006620,0.68581349}; #define MTSAT2_NNIGHT 4 float MTSAT2_NIGHT_MEAN_ERR[MTSAT2_NNIGHT] = {0.12043308,0.30219896,0.49609574,1.0430975}; float MTSAT2_NIGHT_BIAS[MTSAT2_NNIGHT] = {-0.059461016,-0.043822456,-0.10529868,-0.40427330}; float MTSAT2_NIGHT_STDEV[MTSAT2_NNIGHT] = {0.46123141,0.47195938,0.58006620,0.68581349}; #define MSG2_NDAY 5 float MSG2_DAY_MEAN_ERR[MSG2_NDAY] = {0.12826566,0.30539998,0.50723634,0.68193927,1.0801146}; float MSG2_DAY_BIAS[MSG2_NDAY] = {-0.096471168,-0.11493006,-0.18234481,-0.26645958,-0.51300478}; float MSG2_DAY_STDEV[MSG2_NDAY] = {0.45666140,0.49145469,0.60757136,0.81120104,1.0727277}; #define MSG2_NNIGHT 5 float MSG2_NIGHT_MEAN_ERR[MSG2_NNIGHT] = {0.12417339,0.30663668,0.50377016,0.50377016 }; float MSG2_NIGHT_BIAS[MSG2_NNIGHT] = {0.016464286,-0.019051878,-0.19631855,-0.56603110}; float MSG2_NIGHT_STDEV[MSG2_NNIGHT] = {0.32890087,0.32890087,0.64472902,1.0895201}; #define MSG3_NDAY 6 float MSG3_DAY_MEAN_ERR[MSG3_NDAY] = {0.11996787,0.29320928,0.45214627,0.56021188,0.71262558,1.2584690}; float MSG3_DAY_BIAS[MSG3_NDAY] = {-0.12103125,-0.11725302,-0.14648621,-0.15402645,-0.27195335,-0.89151692 }; float MSG3_DAY_STDEV[MSG3_NDAY] = { 0.43307534, 0.49469094, 0.57786118, 0.69555655, 0.97821244, 1.3449739 }; #define MSG3_NNIGHT 6 float MSG3_NIGHT_MEAN_ERR[MSG3_NNIGHT] = {0.11536383,0.29309112,0.45332946,0.55717735,0.69315567,1.2669497}; float MSG3_NIGHT_BIAS[MSG3_NNIGHT] = {-0.056063475,-0.047242705,-0.058804317,-0.23487946,-0.93923971,-2.1975732}; float MSG3_NIGHT_STDEV[MSG3_NNIGHT] = {0.38485646,0.46549090,0.62359092,0.88330263,1.2260118,1.4877082}; /* Routine to get the bias/stdev */ void getBiasMTLS( int satId, int dayTime, float MTLSQual, float pClr, float *pOutBias, float *pOutStdev ) { int ndata = 0; int i = 0; float *pMean = NULL; float *pBias = NULL; float *pStdev = NULL; float xx = 0.; float x1 = 0.; float x2 = 0.; switch(satId){ case 180: /* GOES-13 */ if( 1 == dayTime ){ ndata = GOES13_NDAY; pMean = &GOES13_DAY_MEAN_ERR[0]; pBias = &GOES13_DAY_BIAS[0]; pStdev = &GOES13_DAY_STDEV[0]; } else { ndata = GOES13_NNIGHT; pMean = &GOES13_NIGHT_MEAN_ERR[0]; pBias = &GOES13_NIGHT_BIAS[0]; pStdev = &GOES13_NIGHT_STDEV[0]; } break; case 184: /* GOES-15 */ if( 1 == dayTime ){ ndata = GOES13_NDAY; pMean = &GOES13_DAY_MEAN_ERR[0]; pBias = &GOES13_DAY_BIAS[0]; pStdev = &GOES13_DAY_STDEV[0]; } else { ndata = GOES13_NNIGHT; pMean = &GOES13_NIGHT_MEAN_ERR[0]; pBias = &GOES13_NIGHT_BIAS[0]; pStdev = &GOES13_NIGHT_STDEV[0]; } break; case 84: /* MTSAT 1 */ if( 1 == dayTime ){ ndata = MTSAT1_NDAY; pMean = &MTSAT1_DAY_MEAN_ERR[0]; pBias = &MTSAT1_DAY_BIAS[0]; pStdev = &MTSAT1_DAY_STDEV[0]; } else { ndata = MTSAT1_NNIGHT; pMean = &MTSAT1_NIGHT_MEAN_ERR[0]; pBias = &MTSAT1_NIGHT_BIAS[0]; pStdev = &MTSAT1_NIGHT_STDEV[0]; } case 85: /* MTSAT 2 */ if( 1 == dayTime ){ ndata = MTSAT2_NDAY; pMean = &MTSAT2_DAY_MEAN_ERR[0]; pBias = &MTSAT2_DAY_BIAS[0]; pStdev = &MTSAT2_DAY_STDEV[0]; } else { ndata = MTSAT2_NNIGHT; pMean = &MTSAT2_NIGHT_MEAN_ERR[0]; pBias = &MTSAT2_NIGHT_BIAS[0]; pStdev = &MTSAT2_NIGHT_STDEV[0]; } break; case 52: /* MSG2 */ if( 1 == dayTime ){ ndata = MSG2_NDAY; pMean = &MSG2_DAY_MEAN_ERR[0]; pBias = &MSG2_DAY_BIAS[0]; pStdev = &MSG2_DAY_STDEV[0]; } else { ndata = MSG2_NNIGHT; pMean = &MSG2_NIGHT_MEAN_ERR[0]; pBias = &MSG2_NIGHT_BIAS[0]; pStdev = &MSG2_NIGHT_STDEV[0]; } case 53: /* MSG3 */ if( 1 == dayTime ){ ndata = MSG3_NDAY; pMean = &MSG3_DAY_MEAN_ERR[0]; pBias = &MSG3_DAY_BIAS[0]; pStdev = &MSG3_DAY_STDEV[0]; } else { ndata = MSG3_NNIGHT; pMean = &MSG3_NIGHT_MEAN_ERR[0]; pBias = &MSG3_NIGHT_BIAS[0]; pStdev = &MSG3_NIGHT_STDEV[0]; } break; default: fprintf(stderr,"Satellite ID = %d\n",satId); printError("L2P","Satellite ID not allowed (MTLS Bias)",NULL,NULL); } /* Get bias and standard deviation assuming linear interpolation */ if( MTLSQual < *pMean ){ *pOutBias = *pBias; *pOutStdev = *pStdev; } else if( MTLSQual > *(pMean+ndata-1) ){ *pOutBias = *(pBias+ndata-1); *pOutStdev = *(pStdev+ndata-1); } else { /* Find location */ *pOutBias = -1.; *pOutStdev = -1.; for(i=0;i= *(pMean+i)) && (MTLSQual <= *(pMean+i+1)) ){ /* Do linear interpolation */ x1 = *(pMean+i); x2 = *(pMean+i+1); xx = (x2-x1); *pOutBias = ((x2-MTLSQual) * *(pBias+i)/xx) + ((MTLSQual-x1) * *(pBias+i+1)/xx); *pOutStdev = ((x2-MTLSQual) * *(pStdev+i)/xx) + ((MTLSQual-x1) * *(pStdev+i+1)/xx); break; } } /*if( -1 == *pOutStdev ){ */ if( (-1 == *pOutBias) || (-1 == *pOutStdev) ){ /* printError("L2P","MTLSQual Error not found in interpolation range",NULL,NULL);*/ printError("L2P","MTLSQual Error not found in interpolation range ",&MTLSQual,NULL); } } } /**************************************************************************** * * The next set of routines impliments Andy's method for determining the * bias and standard deviation based on the atmosphere's transmittence * determined from the dBT/dSST value, and the 2m/surface temperature * differential */ /* First come the lookup tables for GOES 9,10,12 */ #define NBIAS 7 #define NBIAS_PROB 4 /* GOES-9 */ int new_coeffs_goes9 = 0; float tauArrayDay_G9[NBIAS] = {0.361363,0.432442,0.500742,0.587665,0.680806, 0.779026,0.851376}; float tauArrayNight_G9[NBIAS] = {0.352912,0.428582,0.501456,0.586004,0.678332, 0.766490,0.841433}; float constDay_G9[NBIAS] = {0.976278,0.838050,0.628772,0.411598,0.375371, 0.258988,-0.148509}; float slopeDay_G9[NBIAS] = {0.399097,0.240611,0.0952944,-0.0549361,-0.0841420, -0.0986134,-0.144536}; float constNight_G9[NBIAS] = {0.127817,0.0365662,0.00886419,0.0200763,0.0848512, 0.0857580,0.0235307}; float slopeNight_G9[NBIAS] = {0.334651,0.271246,0.215081,0.133458,0.133007, 0.124041,0.0664313}; float stdevDay_G9[NBIAS] = {0.737592,0.728649,0.723216,0.701481,0.762940, 0.931213,1.00757}; float stdevNight_G9[NBIAS] = {0.472973,0.430902,0.440563,0.454246,0.486789, 0.590965,0.664016}; /* GOES-10 */ int new_coeffs_goes10 = 0; float tauArrayDay_G10[NBIAS] = {0.356464,0.433784,0.518690,0.583048, 0.680983,0.811744,0.856426}; float tauArrayNight_G10[NBIAS] = {0.419046,0.416488,0.392434,0.384046, 0.507905,0.607350,0.561964}; float constDay_G10[NBIAS] = {-0.677188,-0.488914,-0.254824,-0.121324, 0.128654,0.426907,0.464489}; float slopeDay_G10[NBIAS] = {0.392443,0.324112,0.224065,0.0894843,0.0781492, 0.0807758,0.0497196}; float constNight_G10[NBIAS] = {-0.221139,-0.159550,-0.0400801,0.0616851, 0.232522,0.487704,0.530731}; float slopeNight_G10[NBIAS] = {0.0914795,0.0492509,0.0735530,0.109503, 0.162550,0.172795,0.140105}; float stdevDay_G10[NBIAS] = {0.691345,0.632672,0.577976,0.556813,0.605931, 0.609013,0.589098}; float stdevNight_G10[NBIAS] = {0.416660,0.415853,0.390858,0.379280, 0.490739,0.582398,0.543434}; /* GOES-11 */ int new_coeffs_goes11 = 0; float tauArrayDay_G11[NBIAS] = {0.352357,0.438551,0.515894,0.570034, 0.703946,0.785589,0.838659}; float tauArrayNight_G11[NBIAS] = {0.355200,0.434927,0.511409,0.569991, 0.672099,0.786754,0.844136}; float constDay_G11[NBIAS] = {-0.219083,-0.0330325,0.157057,0.322864, 0.668739,0.735443,0.668854}; float slopeDay_G11[NBIAS] = {0.564336,0.496418,0.355044,0.293577, 0.150152,0.0860100,0.0920904}; float constNight_G11[NBIAS] = {-0.341596,-0.312718,-0.116543,0.0236273, 0.306915,0.610072,0.686187}; float slopeNight_G11[NBIAS] = {0.233417,0.130715,0.253366,0.300645, 0.345915,0.291719,0.238149}; float stdevDay_G11[NBIAS] = {0.665526,0.596661,0.561807,0.568816, 0.610876,0.610882,0.683165}; float stdevNight_G11[NBIAS] = {0.416660,0.415853,0.390858,0.379280, 0.490739,0.582398,0.543434}; /* GOES-12 */ int new_coeffs_goes12 = 0; float tauArrayDay_G12[NBIAS] = {0.364964,0.429084,0.500204,0.598232, 0.689273,0.782148,0.856534}; float tauArrayNight_G12[NBIAS] = {0.363280,0.430231,0.499296,0.597691, 0.690673,0.794692,0.872125}; float constDay_G12[NBIAS] = {-0.374868,-0.364560,-0.337227,-0.356830, -0.391169,-0.389511,-0.339609}; float slopeDay_G12[NBIAS] = {0.0941708,0.145087,0.173560,0.150618, 0.123108,0.0952590,0.0888319}; float constNight_G12[NBIAS] = {-0.576031,-0.547570,-0.496679,-0.369300, -0.183615,-0.0265943,0.0403807}; float slopeNight_G12[NBIAS] = {0.211028,0.250964,0.232512,0.220121,0.223894, 0.157141,0.130205}; float stdevDay_G12[NBIAS] = {0.832308,0.737539,0.724082,0.749948, 0.808908,0.839862,0.768678}; float stdevNight_G12[NBIAS] = {0.635881,0.515487,0.482030,0.520677,0.613761, 0.651150,0.579933}; /* GOES-13 */ int new_coeffs_goes13 = 0; #ifdef OLD_G13_COEFFS float tauArrayDay_G13[NBIAS] = {0.359431,0.459431,0.538767,0.635820,0.713233,0.811898,0.854530}; float tauArrayNight_G13[NBIAS] = {0.356488,0.456488,0.537534,0.635989,0.703221,0.807931,0.856714}; float constDay_G13[NBIAS] = {-0.475932,-0.475932,-0.396601,-0.528897,-0.534135,-0.493162,-0.479308}; float slopeDay_G13[NBIAS] = {0.756163,0.756163,0.529927,0.399378,0.218199,0.0876433,0.0365910}; float constNight_G13[NBIAS] = {-1.17031,-1.17031,-0.815593,-0.754129,-0.719599,-0.479679,-0.322826}; float slopeNight_G13[NBIAS] = {0.133541,0.133541,0.425010,0.424939,0.416138,0.425239,0.404952}; float stdevDay_G13[NBIAS] = {0.841057,0.841057,1.14205,1.04865,1.05287,1.48009,1.70472}; float stdevNight_G13[NBIAS] = {1.27082,1.27082,1.38483,0.910733,0.761091,0.859492,0.887184}; #else float tauArrayDay_G13[NBIAS] = {0.547514,0.547514,0.547514,0.623922,0.716145,0.845566,0.869016}; float tauArrayNight_G13[NBIAS] = {0.540129,0.540129,0.540129,0.607597,0.696764,0.837278,0.867803}; float constDay_G13[NBIAS] = {0.239007,0.239007,0.239007,0.471898,0.429443,0.214712,0.191647}; float slopeDay_G13[NBIAS] = {0.173595,0.173595,0.173595,0.366129,0.363242,0.135607,0.0909835}; float constNight_G13[NBIAS] = {-0.254484,-0.254484,-0.254484,-0.401498,-0.309745,0.121954,0.236334}; float slopeNight_G13[NBIAS] = {0.0968533,0.0968533,0.0968533,0.0251188,0.215387,0.306057,0.310285}; float stdevDay_G13[NBIAS] = {0.941066,0.941066,0.941066,1.39632,1.35675,1.24420,1.24153}; float stdevNight_G13[NBIAS] = {0.654864,0.654864,0.654864,0.652658,0.695500,0.764373,0.752513}; #endif int new_coeffs_goes15 = 1; /* MTSAT-2 coefficients */ float probArrayDay_G15[NBIAS_PROB][NBIAS] = {{0.896588,0.894435,0.885810,0.885930,0.888507,0.888507,0.888507},{0.972444,0.972745,0.974541,0.975333,0.975272,0.975272,0.975272},{0.996222,0.996123,0.995878,0.996160,0.996612,0.996612,0.996612},{0.999965,0.999964,0.999951,0.999940,0.999916,0.999877,0.999877}}; float probArrayNight_G15[NBIAS_PROB][NBIAS] = {{0.891824,0.891398,0.892363,0.893531,0.890941,0.885290,0.885290},{0.975940,0.975789,0.974561,0.973037,0.972191,0.972191,0.972191},{0.996051,0.995898,0.995324,0.995143,0.994693,0.994693,0.994693},{0.999591,0.999585,0.999520,0.999412,0.999412,0.999412,0.999412}}; float tauArrayDay_G15_P[NBIAS_PROB][NBIAS] = {{0.890362,0.856450,0.712721,0.630086,0.549932,0.549932,0.549932},{0.893722,0.867467,0.709302,0.621743,0.542218,0.542218,0.542218},{0.890185,0.851848,0.709383,0.618418,0.538023,0.538023,0.538023},{0.890353,0.858257,0.709091,0.629263,0.545803,0.460555,0.460555}}; float tauArrayNight_G15_P[NBIAS_PROB][NBIAS] = {{0.879811,0.825196,0.700839,0.624555,0.544037,0.469854,0.469854},{0.884171,0.840928,0.704908,0.630074,0.557852,0.557852,0.557852},{0.891087,0.858898,0.716555,0.638908,0.554656,0.554656,0.554656},{0.892720,0.866502,0.730067,0.658797,0.658797,0.658797,0.658797}}; float constDay_G15_P[NBIAS_PROB][NBIAS] = {{-0.221114,-0.325329,-0.747427,-0.820809,-0.895599,-0.895599,-0.895599},{-0.197292,-0.266412,-0.484008,-0.766115,-0.871642,-0.871642,-0.871642},{-0.727778,-0.810440,-0.939935,-0.935802,-1.13404,-1.13404,-1.13404},{-0.749239,-0.730689,-0.665249,-0.828150,-1.15148,-1.20408,-1.20408}}; float slopeDay_G15_P[NBIAS_PROB][NBIAS] = {{0.339988,0.354462,0.121336,-0.298964,-0.126991,-0.126991,-0.126991},{0.260105,0.313045,0.497260,0.150863,0.406974,0.406974,0.406974},{0.255123,0.285817,0.482822,0.311248,-0.160640,-0.160640,-0.160640},{0.318711,0.322499,0.382440,0.354253,0.267806,0.633179,0.633179}}; float constNight_G15_P[NBIAS_PROB][NBIAS] = {{-1.12837,-1.18794,-1.45191,-1.41693,-1.45131,-1.62287,-1.62287},{-0.949134,-0.989608,-1.25312,-1.28018,-1.22049,-1.22049,-1.22049},{-0.841238,-0.864615,-1.13045,-1.13282,-1.15686,-1.15686,-1.15686},{-0.749875,-0.762817,-1.08174,-1.11575,-1.11575,-1.11575,-1.11575}}; float slopeNight_G15_P[NBIAS_PROB][NBIAS] = {{0.455270,0.424668,0.157852,0.193718,0.0702591,-0.307456,-0.307456},{0.460148,0.453024,0.237260,0.205590,0.282312,0.282312,0.282312},{0.442902,0.438264,0.220385,0.270241,0.279394,0.279394,0.279394},{0.435952,0.430120,0.111134,0.111008,0.111008,0.111008,0.111008}}; float stdevDay_G15_P[NBIAS_PROB][NBIAS] = {{1.20603,1.35068,1.50251,1.07656,1.13287,1.13287,1.13287},{1.27125,1.25877,1.25269,1.40464,1.78233,1.78233,1.78233},{1.82127,1.74895,1.48013,1.22130,0.885423,0.885423,0.885423},{0.897125,0.864811,0.717533,0.879430,1.20912,1.41051,1.41051}}; float stdevNight_G15_P[NBIAS_PROB][NBIAS] = {{0.816385,0.705052,0.436915,0.385035,0.365440,0.456892,0.456892},{0.796105,0.696302,0.371221,0.377136,0.399099,0.399099,0.399099},{0.820281,0.747948,0.362686,0.351159,0.351117,0.351117,0.351117},{0.765090,0.712925,0.335574,0.271424,0.271424,0.271424,0.271424}}; /* MTSAT - using new (prob variable) coefficients */ #ifdef OLD_MTSAT_COEFFS int new_coeffs_mtsat = 0; #else int new_coeffs_mtsat = 1; #endif /* MTSAT */ float tauArrayDay_MTSAT[NBIAS] = {0.346828,0.420803,0.485222,0.584194, 0.695556,0.788623,0.849561}; float tauArrayNight_MTSAT[NBIAS] = {0.344928,0.417574,0.481000,0.581872, 0.696414,0.789087,0.851854}; float constDay_MTSAT[NBIAS] = {0.0790447,0.0858210,0.0224725,-0.0916170, -0.232819,-0.298037,-0.529132}; float slopeDay_MTSAT[NBIAS] = {0.130454,0.0817362,0.0216837,-0.0758651, -0.127496,-0.0485068,-0.0505161}; float constNight_MTSAT[NBIAS] = {-0.00136553,-0.0213738,-0.0633295, -0.0802552,-0.111762,-0.131259,-0.144858}; float slopeNight_MTSAT[NBIAS] = {-0.00672808,-0.00528396,-0.0118487, 0.0102000,0.00739864,0.00697242,0.00972383}; float stdevDay_MTSAT[NBIAS] = {1.03901,0.947220,0.906716,0.884218,0.908570, 1.04282,1.14006}; float stdevNight_MTSAT[NBIAS] = {0.583723,0.504723,0.480892,0.510696, 0.518432,0.518372,0.550007}; /* New probability thresholded MTSAT SSES - divided into 0.8-0.95, 0.95-0.99, 0.99-0.999, 0.999-1. */ float probArrayDay_MTSAT[NBIAS_PROB][NBIAS] = {{0.889179,0.890775,0.891210,0.890845,0.889409,0.887290,0.886926},{0.975048,0.975026,0.975585,0.975582,0.975060,0.975664,0.976231},{0.996077,0.995981,0.995952,0.995822,0.995736,0.995799,0.995677},{0.999821,0.999853,0.999864,0.999830,0.999832,0.999830,0.999825}}; float probArrayNight_MTSAT[NBIAS_PROB][NBIAS] = {{0.884610,0.883293,0.882374,0.882334,0.884472,0.887748,0.889586},{0.973846,0.974271,0.974400,0.974441,0.974752,0.974907,0.975349},{0.995606,0.995923,0.996118,0.996183,0.996101,0.995952,0.995859},{0.999569,0.999557,0.999563,0.999590,0.999605,0.999582,0.999537}}; float tauArrayDay_MTSAT_P[NBIAS_PROB][NBIAS] = {{0.347357,0.420283,0.481578,0.584187,0.698688,0.787518,0.850059},{0.347853,0.420987,0.486089,0.582888,0.694718,0.799169,0.853507},{0.346054,0.416559,0.482527,0.589437,0.701205,0.793640,0.854242},{0.346721,0.421843,0.485608,0.583330,0.694236,0.783547,0.846210}}; float tauArrayNight_MTSAT_P[NBIAS_PROB][NBIAS] = {{0.345450,0.414990,0.477176,0.580058,0.693734,0.789910,0.851720},{0.343993,0.415309,0.480019,0.582322,0.697486,0.791321,0.852822},{0.343347,0.419697,0.484863,0.582984,0.698361,0.790712,0.852549},{0.347651,0.418088,0.478009,0.580021,0.692403,0.782891,0.848644}}; float constDay_MTSAT_P[NBIAS_PROB][NBIAS] = {{0.0123020,0.00650834,-0.0655521,-0.256660,-0.384332,-0.495415,-0.728498},{0.0142619,0.0244465,-0.0354739,-0.173433,-0.342488,-0.523396,-0.776290},{-0.00722106,-0.00596656,-0.0884868,-0.194201,-0.301124,-0.445977,-0.885645},{0.127217,0.126743,0.0644956,-0.0416545,-0.185276,-0.191306,-0.288040}}; float slopeDay_MTSAT_P[NBIAS_PROB][NBIAS] = {{0.238596,0.130463,0.106938,-0.0157086,-0.0852790,-0.0760244,-0.0976610},{0.165915,0.132710,0.0612735,-0.0559162,-0.124086,-0.101425,-0.112346},{0.140287,0.129494,0.0330198,-0.112738,-0.120817,-0.0257806,-0.0538562},{0.123434,0.0553246,0.00804604,-0.0690094,-0.132023,-0.0453742,-0.0205610}}; float constNight_MTSAT_P[NBIAS_PROB][NBIAS] = {{-0.363156,-0.302308,-0.286420,-0.320111,-0.320804,-0.229592,-0.191487},{-0.141410,-0.131377,-0.156824,-0.152016,-0.157577,-0.154939,-0.105294},{-0.0208766,-0.0494505,-0.0887082,-0.0930200,-0.109194,-0.123622,-0.162606},{0.211254,0.142457,0.0716707,0.0155405,-0.0576362,-0.110193,-0.182968}}; float slopeNight_MTSAT_P[NBIAS_PROB][NBIAS] = {{-0.132155,-0.124245,-0.0789476,-0.0489172,-0.0404837,0.00309738,0.0148036},{-0.0491752,-0.0258771,-0.0124117,0.0148066,-0.000722130,0.00671586,0.0236623},{-0.0299099,-0.0386174,-0.0419123,-0.00503642,0.00845810,0.00747703,0.00250814},{0.0896556,0.0657436,0.0290058,0.0210694,0.0162488,0.00523284,-0.00352706}}; float stdevDay_MTSAT_P[NBIAS_PROB][NBIAS] = {{1.15701,1.01172,0.950885,1.02172,1.14542,1.18152,1.15867},{1.09701,1.00986,0.986253,0.992851,1.08222,1.20255,1.15982},{1.09711,1.08183,1.03609,0.936853,0.961405,1.26167,1.54443},{0.999386,0.884163,0.844688,0.834571,0.833876,0.883943,0.906573}}; float stdevNight_MTSAT_P[NBIAS_PROB][NBIAS] = {{0.687858,0.597498,0.570805,0.631921,0.680512,0.711373,0.768826},{0.607443,0.535354,0.519298,0.576690,0.589785,0.572738,0.583195},{0.558662,0.503155,0.478305,0.488844,0.520584,0.515023,0.545341},{0.535194,0.440279,0.411780,0.436511,0.402486,0.427695,0.492234}}; int new_coeffs_mtsat1 = 1; /* MTSAT-1 coefficients coppied from MTSAT-2 bellow */ float probArrayDay_MTSAT1[NBIAS_PROB][NBIAS] = {{0.892879,0.892879,0.887009,0.887832,0.890917,0.890561,0.892410},{0.975035,0.975035,0.975192,0.974857,0.974620,0.975088,0.974133},{0.995964,0.995964,0.995970,0.996004,0.996067,0.996058,0.996017},{0.999964,0.999951,0.999955,0.999959,0.999963,0.999950,0.999941}}; float probArrayNight_MTSAT1[NBIAS_PROB][NBIAS] = {{0.890223,0.887806,0.888529,0.889825,0.885404,0.883992},{0.975096,0.975081,0.975441,0.975559,0.975528,0.974167,0.973505},{0.996300,0.996030,0.996165,0.996143,0.996111,0.996201,0.996171},{0.999647,0.999732,0.999826,0.999849,0.999845,0.999765,0.999589}}; float tauArrayDay_MTSAT1_P[NBIAS_PROB][NBIAS] = {{0.458827,0.458827,0.541249,0.588340,0.670681,0.811232,0.867259},{0.454585,0.454585,0.540652,0.583120,0.666728,0.810672,0.862679},{0.462296,0.462296,0.541680,0.590024,0.664197,0.807784,0.864974},{0.341719,0.465823,0.543701,0.589003,0.672332,0.809173,0.864444}}; float tauArrayNight_MTSAT1_P[NBIAS_PROB][NBIAS] = {{0.358736,0.459332,0.540628,0.584461,0.672363,0.828470,0.867223},{0.352382,0.459387,0.539532,0.586153,0.670754,0.821739,0.868671},{0.346478,0.461809,0.543784,0.586972,0.665000,0.821307,0.867142},{0.864083,0.793702,0.658009,0.585399,0.546372,0.469044,0.350140}}; float constDay_MTSAT1_P[NBIAS_PROB][NBIAS] = {{-0.580545,-0.580545,-0.402047,-0.302657,-0.276859,-0.0684326,-0.0470084},{-0.423700,-0.423700,-0.349808,-0.208088,-0.170566,-0.255334,-0.145164},{-0.736908,-0.736908,-0.677904,-0.414111,-0.218973,-0.348371,-0.356945},{-0.385661,-0.396672,-0.0864515,0.0849145,0.326573,0.286142,0.164524}}; float slopeDay_MTSAT1_P[NBIAS_PROB][NBIAS] = {{0.675647,0.675647,0.551073,0.453368,0.136949,0.0513524,0.00440242},{1.03407,1.03407,0.482739,0.431806,0.121002,-0.0622398,0.0225920},{0.755943,0.755943,0.212591,0.143481,-0.0450738,-0.127220,-0.0773517},{0.665743,0.742486,0.424467,0.302553,0.144884,-0.000618217,-0.0494517}}; float constNight_MTSAT1_P[NBIAS_PROB][NBIAS] = {{-0.811800,-0.301759,-0.299763,-0.310866,-0.279436,-0.177536,-0.122248},{0.0464397,-0.174931,-0.182700,-0.178423,-0.0864035,0.0228643,0.0419102},{0.151174,0.0543383,-0.0182403,0.00333367,0.0846459,0.331485,0.353164},{0.471861,0.443253,0.335218,0.269131,0.258112,0.322784,0.440594}}; float slopeNight_MTSAT1_P[NBIAS_PROB][NBIAS] = {{0.462971,0.232270,0.167608,0.0917647,-0.00138536,-0.0660901,-0.0686854},{0.640030,0.220168,0.146280,0.100344,0.0824346,-0.0313875,-0.0609932},{0.571864,0.293562,0.160348,0.136413,0.0823744,0.0457680,0.0119413},{0.0421150,0.0796281,0.138145,0.133618,0.148019,0.306183,0.556525}}; float stdevDay_MTSAT1_P[NBIAS_PROB][NBIAS] = {{1.62043,1.62043,1.31298,1.24343,1.17212,1.02650,1.08200},{1.21293,1.21293,1.27200,1.22742,1.15772,1.32427,1.35860},{1.63378,1.63378,1.31541,1.20411,1.05623,1.15059,1.20743},{1.12985,1.17375,1.01956,0.963307,0.875261,0.917421,0.941849}}; float stdevNight_MTSAT1_P[NBIAS_PROB][NBIAS] = {{1.78862,1.05432,0.822099,0.800474,0.778860,0.787783,0.775995},{1.19181,0.833355,0.674268,0.632681,0.619918,0.678051,0.690960},{0.886786,0.729767,0.589372,0.569834,0.567384,0.518483,0.481134},{0.390577,0.415907,0.391895,0.384846,0.400397,0.504624,0.500785}}; int new_coeffs_mtsat2 = 1; /* MTSAT-2 coefficients */ float probArrayDay_MTSAT2[NBIAS_PROB][NBIAS] = {{0.892879,0.892879,0.887009,0.887832,0.890917,0.890561,0.892410},{0.975035,0.975035,0.975192,0.974857,0.974620,0.975088,0.974133},{0.995964,0.995964,0.995970,0.996004,0.996067,0.996058,0.996017},{0.999964,0.999951,0.999955,0.999959,0.999963,0.999950,0.999941}}; float probArrayNight_MTSAT2[NBIAS_PROB][NBIAS] = {{0.890223,0.887806,0.888529,0.889825,0.885404,0.883992},{0.975096,0.975081,0.975441,0.975559,0.975528,0.974167,0.973505},{0.996300,0.996030,0.996165,0.996143,0.996111,0.996201,0.996171},{0.999647,0.999732,0.999826,0.999849,0.999845,0.999765,0.999589}}; float tauArrayDay_MTSAT2_P[NBIAS_PROB][NBIAS] = {{0.458827,0.458827,0.541249,0.588340,0.670681,0.811232,0.867259},{0.454585,0.454585,0.540652,0.583120,0.666728,0.810672,0.862679},{0.462296,0.462296,0.541680,0.590024,0.664197,0.807784,0.864974},{0.341719,0.465823,0.543701,0.589003,0.672332,0.809173,0.864444}}; float tauArrayNight_MTSAT2_P[NBIAS_PROB][NBIAS] = {{0.358736,0.459332,0.540628,0.584461,0.672363,0.828470,0.867223},{0.352382,0.459387,0.539532,0.586153,0.670754,0.821739,0.868671},{0.346478,0.461809,0.543784,0.586972,0.665000,0.821307,0.867142},{0.864083,0.793702,0.658009,0.585399,0.546372,0.469044,0.350140}}; float constDay_MTSAT2_P[NBIAS_PROB][NBIAS] = {{-0.580545,-0.580545,-0.402047,-0.302657,-0.276859,-0.0684326,-0.0470084},{-0.423700,-0.423700,-0.349808,-0.208088,-0.170566,-0.255334,-0.145164},{-0.736908,-0.736908,-0.677904,-0.414111,-0.218973,-0.348371,-0.356945},{-0.385661,-0.396672,-0.0864515,0.0849145,0.326573,0.286142,0.164524}}; float slopeDay_MTSAT2_P[NBIAS_PROB][NBIAS] = {{0.675647,0.675647,0.551073,0.453368,0.136949,0.0513524,0.00440242},{1.03407,1.03407,0.482739,0.431806,0.121002,-0.0622398,0.0225920},{0.755943,0.755943,0.212591,0.143481,-0.0450738,-0.127220,-0.0773517},{0.665743,0.742486,0.424467,0.302553,0.144884,-0.000618217,-0.0494517}}; float constNight_MTSAT2_P[NBIAS_PROB][NBIAS] = {{-0.811800,-0.301759,-0.299763,-0.310866,-0.279436,-0.177536,-0.122248},{0.0464397,-0.174931,-0.182700,-0.178423,-0.0864035,0.0228643,0.0419102},{0.151174,0.0543383,-0.0182403,0.00333367,0.0846459,0.331485,0.353164},{0.471861,0.443253,0.335218,0.269131,0.258112,0.322784,0.440594}}; float slopeNight_MTSAT2_P[NBIAS_PROB][NBIAS] = {{0.462971,0.232270,0.167608,0.0917647,-0.00138536,-0.0660901,-0.0686854},{0.640030,0.220168,0.146280,0.100344,0.0824346,-0.0313875,-0.0609932},{0.571864,0.293562,0.160348,0.136413,0.0823744,0.0457680,0.0119413},{0.0421150,0.0796281,0.138145,0.133618,0.148019,0.306183,0.556525}}; float stdevDay_MTSAT2_P[NBIAS_PROB][NBIAS] = {{1.62043,1.62043,1.31298,1.24343,1.17212,1.02650,1.08200},{1.21293,1.21293,1.27200,1.22742,1.15772,1.32427,1.35860},{1.63378,1.63378,1.31541,1.20411,1.05623,1.15059,1.20743},{1.12985,1.17375,1.01956,0.963307,0.875261,0.917421,0.941849}}; float stdevNight_MTSAT2_P[NBIAS_PROB][NBIAS] = {{1.78862,1.05432,0.822099,0.800474,0.778860,0.787783,0.775995},{1.19181,0.833355,0.674268,0.632681,0.619918,0.678051,0.690960},{0.886786,0.729767,0.589372,0.569834,0.567384,0.518483,0.481134},{0.390577,0.415907,0.391895,0.384846,0.400397,0.504624,0.500785}}; /* Use MSG-2 new coefficients */ #ifdef OLD_MSG_COEFFS int new_coeffs_msg2 = 0; #else int new_coeffs_msg2 = 1; #endif /* MSG-2 */ float tauArrayDay_MSG2[NBIAS] = {0.314550,0.420008,0.511455,0.606707, 0.696563,0.780534,0.837080}; float tauArrayNight_MSG2[NBIAS] = {0.307901,0.412722,0.510859,0.603048, 0.704556,0.787009,0.840462}; float constDay_MSG2[NBIAS] = {-0.130997,-0.0318350,-0.0793432,-0.115677, -0.179805,-0.255787,-0.266595}; float slopeDay_MSG2[NBIAS] = {0.332341,0.178034,0.0159455,-0.0605023, -0.145127,-0.160100,-0.0935942}; float constNight_MSG2[NBIAS] = {-0.0973243,-0.189603,-0.230322,-0.209304, -0.110773,0.0419614,0.179236}; float slopeNight_MSG2[NBIAS] = {0.346574,0.248452,0.171780,0.124843,0.103473, 0.123867,0.136059}; float stdevDay_MSG2[NBIAS] = {0.904453,0.806779,0.831525,0.825025,0.834243, 0.884183,0.923727}; float stdevNight_MSG2[NBIAS] = {0.939830,0.817919,0.784054,0.809669,0.768025, 0.744377,0.754268}; /* New coefficients for MSG2 */ #ifdef MSG_COEF_FIRST_TRY float probArrayDay_MSG2[NBIAS_PROB][NBIAS] = {{0.887547,0.890851,0.894018,0.891564,0.891154,0.895382,0.891276},{0.974942,0.974332,0.975602,0.976632,0.975292,0.974994,0.977377},{0.996254,0.995624,0.995672,0.996232,0.996227,0.996147,0.996463},{0.999869,0.999877,0.999862,0.999853,0.999850,0.999851,0.999859}}; float probArrayNight_MSG2[NBIAS_PROB][NBIAS] = {{0.891997,0.888326,0.888659,0.890372,0.891661,0.891779,0.891470},{0.975270,0.975418,0.975187,0.975278,0.974764,0.974526,0.975098},{0.995814,0.995845,0.995788,0.995700,0.995726,0.995960,0.996246},{0.999702,0.999712,0.999726,0.999735,0.999723,0.999699,0.999685}}; float tauArrayDay_MSG2_P[NBIAS_PROB][NBIAS] = {{0.312200,0.417926,0.512627,0.601298,0.704924,0.781200,0.839406},{0.312975,0.421255,0.508218,0.603634,0.700305,0.780423,0.832313},{0.323691,0.418927,0.512772,0.604864,0.701746,0.786448,0.841841},{0.313788,0.419946,0.511812,0.607382,0.695314,0.779678,0.837293}}; float tauArrayNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.304274,0.415066,0.505979,0.606817,0.709584,0.785154,0.839399},{0.304605,0.413666,0.507693,0.606070,0.707762,0.783781,0.838657},{0.307831,0.410343,0.511305,0.602437,0.703892,0.788212,0.841592},{0.310591,0.414093,0.512985,0.601204,0.702227,0.788886,0.840916}}; float constDay_MSG2_P[NBIAS_PROB][NBIAS] = {{-0.229574,-0.191087,-0.445556,-0.478594,-0.375241,-0.471156,-0.744909},{0.232056,-0.0306725,-0.353730,-0.417469,-0.565764,-0.697771,-0.620365},{-0.579094,-0.570791,-0.639485,-0.552455,-0.537639,-0.752596,-0.840732},{-0.115317,0.0461835,0.0461862,-0.0158708,-0.0729978,-0.100741,-0.119188}}; float slopeDay_MSG2_P[NBIAS_PROB][NBIAS] = {{0.442143,0.335897,-0.143489,-0.182639,-0.0920372,-0.163668,-0.260161},{0.760138,0.383800,0.161447,0.0705928,-0.0758799,-0.205991,-0.183782},{0.369690,-0.0565904,-0.0161108,-0.113601,-0.235330,-0.257873,-0.122688},{0.302185,0.192920,0.0137354,-0.0687899,-0.141896,-0.137585,-0.0795861}}; float constNight_MSG2_P[NBIAS_PROB][NBIAS] = {{-0.618860,-0.774110,-0.865725,-0.843158,-0.730721,-0.503794,-0.177249},{-0.460647,-0.470139,-0.539639,-0.534815,-0.403478,-0.234594,-0.0437389},{-0.331605,-0.279752,-0.241671,-0.238321,-0.107793,0.0643218,0.179566},{0.335418,0.107525,0.00716294,0.0611509,0.136486,0.240464,0.329549}}; float slopeNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.205630,0.0677421,0.00502583,0.0307784,0.0154946,0.0793491,0.176670},{0.269025,0.188822,0.105296,0.0783564,0.0609888,0.0778087,0.111832},{0.133473,0.211388,0.225469,0.130550,0.0997079,0.127676,0.145280},{0.501781,0.311292,0.155819,0.146592,0.126530,0.124623,0.114873}}; float stdevDay_MSG2_P[NBIAS_PROB][NBIAS] = {{1.19929,1.25186,1.12093,1.11477,1.10204,1.17468,1.29512},{0.833210,0.840721,0.839270,0.926651,1.04692,1.10575,1.17091},{1.30481,0.798897,0.962638,1.04807,1.01886,1.05145,1.05825},{0.844057,0.780791,0.765368,0.742130,0.733442,0.758978,0.795076}}; float stdevNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.907566,0.909181,0.913858,0.974854,0.971938,0.944861,0.879330},{0.968956,0.840370,0.809461,0.893143,0.872679,0.843599,0.861292},{0.876440,0.834279,0.807832,0.775907,0.712778,0.719318,0.732442},{0.823643,0.694670,0.653640,0.673136,0.621754,0.591061,0.629726}}; #else float probArrayDay_MSG2[NBIAS_PROB][NBIAS] = {{0.888547,0.888477,0.892456,0.891365,0.891988,0.895103,0.896105},{0.974117,0.974698,0.974066,0.975302,0.976344,0.975129,0.974579},{0.996220,0.996081,0.995747,0.995693,0.996075,0.996200,0.996062},{0.999983,0.999984,0.999982,0.999984,0.999984,0.999985,0.999988}}; float probArrayNight_MSG2[NBIAS_PROB][NBIAS] = {{0.889686,0.887512,0.890094,0.890914,0.892078,0.891430,0.888492},{0.974865,0.975656,0.975417,0.975776,0.975234,0.974593,0.974844},{0.995830,0.995879,0.995992,0.995953,0.995951,0.996078,0.996194},{0.999857,0.999857,0.999878,0.999893,0.999881,0.999853,0.999839}}; float tauArrayDay_MSG2_P[NBIAS_PROB][NBIAS] = {{0.324602,0.416618,0.507928,0.605330,0.706155,0.785592,0.838530},{0.306979,0.421951,0.503156,0.605900,0.700293,0.783621,0.836193},{0.317432,0.418441,0.510837,0.603036,0.701600,0.783315,0.835164},{0.318364,0.423997,0.514034,0.609710,0.703721,0.786207,0.841911}}; float tauArrayNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.309682,0.414954,0.505949,0.605670,0.707616,0.786718,0.836898},{0.308795,0.413410,0.507603,0.604706,0.708649,0.786161,0.837455},{0.306782,0.414150,0.508090,0.602326,0.702251,0.788569,0.840009},{0.310733,0.415762,0.513589,0.602053,0.697009,0.790520,0.843508}}; float constDay_MSG2_P[NBIAS_PROB][NBIAS] = {{-0.376744,-0.301169,-0.401568,-0.435039,-0.593038,-0.585056,-0.595208},{0.0482780,-0.125431,-0.402471,-0.640007,-0.720588,-0.707597,-0.693210},{-0.540416,-0.528290,-0.386481,-0.441960,-0.555032,-0.665315,-0.736960},{ 0.0164408,-0.0344198, 0.0000155, 0.0106857,-0.0179222, 0.0265621, 0.0840021}}; float slopeDay_MSG2_P[NBIAS_PROB][NBIAS] = {{0.506721, 0.205913, -0.199585,-0.0389948,-0.0848973, -0.130856, -0.125885},{1.00180, 0.650059, 0.381212,0.0933715,-0.142577,-0.235566,-0.239068},{0.291239, 0.225446, 0.188209, 0.0869040,-0.0578114, -0.172481, -0.167732},{0.323330, 0.158842, 0.0489850,-0.0678940, -0.135649, -0.129956, -0.116706}}; float constNight_MSG2_P[NBIAS_PROB][NBIAS] = {{-0.728557,-0.802580,-0.845953,-0.869644,-0.758978,-0.524964,-0.287198},{ -0.578636, -0.534382, -0.545936, -0.574475, -0.425940, -0.217384,-0.0754272},{-0.324537,-0.352071,-0.290668,-0.237832,-0.116301,0.0836247, 0.247517},{0.362035,0.193630,0.144031,0.139718,0.167273,0.290755,0.395634}}; float slopeNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.144721,0.0983967,0.0662757,0.0337547,0.0281491,0.0918672, 0.147579},{0.109641, 0.138527, 0.115475,0.0750003,0.0748850,0.0935038,0.0921714},{0.133076,0.172025,0.222471,0.177562,0.144870,0.139214,0.131063},{0.507777, 0.411701, 0.349673, 0.223007, 0.140696, 0.120895,0.0822841}}; float stdevDay_MSG2_P[NBIAS_PROB][NBIAS] = {{1.32742,1.25394,1.17877,1.26660,1.22669,1.15129,1.25907},{1.44941,1.18646,1.05226,1.23384,1.15502,1.13482,1.39846},{1.43436,1.39360,1.05844,1.16297,1.11688,1.07288,1.29480},{ 1.02060,0.972632,0.863572,0.762154,0.715097,0.727109,0.752480}}; float stdevNight_MSG2_P[NBIAS_PROB][NBIAS] = {{0.962255,0.942267,0.928625,0.964525,0.977554,0.980117,0.969031},{1.01054,0.925202,0.866742,0.924138,0.889855,0.846424,0.865206},{0.962126,0.866818,0.810073,0.767430,0.714048,0.712792,0.718088},{0.938827,0.846861,0.771563,0.663386,0.562015,0.564683,0.592477}}; #endif /* Use MSG-3 new coefficients *//*COPIED FROM MSG - RIP,GSR 2013apr29 */ #ifdef OLD_MSG_COEFFS int new_coeffs_msg3 = 0; #else int new_coeffs_msg3 = 1; #endif /* MSG-3 */ float tauArrayDay_MSG3[NBIAS] = {0.314550,0.420008,0.511455,0.606707, 0.696563,0.780534,0.837080}; float tauArrayNight_MSG3[NBIAS] = {0.307901,0.412722,0.510859,0.603048, 0.704556,0.787009,0.840462}; float constDay_MSG3[NBIAS] = {-0.130997,-0.0318350,-0.0793432,-0.115677, -0.179805,-0.255787,-0.266595}; float slopeDay_MSG3[NBIAS] = {0.332341,0.178034,0.0159455,-0.0605023, -0.145127,-0.160100,-0.0935942}; float constNight_MSG3[NBIAS] = {-0.0973243,-0.189603,-0.230322,-0.209304, -0.110773,0.0419614,0.179236}; float slopeNight_MSG3[NBIAS] = {0.346574,0.248452,0.171780,0.124843,0.103473, 0.123867,0.136059}; float stdevDay_MSG3[NBIAS] = {0.904453,0.806779,0.831525,0.825025,0.834243, 0.884183,0.923727}; float stdevNight_MSG3[NBIAS] = {0.939830,0.817919,0.784054,0.809669,0.768025, 0.744377,0.754268}; /* New coefficients for MSG3 */ #ifdef MSG_COEF_FIRST_TRY float probArrayDay_MSG3[NBIAS_PROB][NBIAS] = {{0.887547,0.890851,0.894018,0.891564,0.891154,0.895382,0.891276},{0.974942,0.974332,0.975602,0.976632,0.975292,0.974994,0.977377},{0.996254,0.995624,0.995672,0.996232,0.996227,0.996147,0.996463},{0.999869,0.999877,0.999862,0.999853,0.999850,0.999851,0.999859}}; float probArrayNight_MSG3[NBIAS_PROB][NBIAS] = {{0.891997,0.888326,0.888659,0.890372,0.891661,0.891779,0.891470},{0.975270,0.975418,0.975187,0.975278,0.974764,0.974526,0.975098},{0.995814,0.995845,0.995788,0.995700,0.995726,0.995960,0.996246},{0.999702,0.999712,0.999726,0.999735,0.999723,0.999699,0.999685}}; float tauArrayDay_MSG3_P[NBIAS_PROB][NBIAS] = {{0.312200,0.417926,0.512627,0.601298,0.704924,0.781200,0.839406},{0.312975,0.421255,0.508218,0.603634,0.700305,0.780423,0.832313},{0.323691,0.418927,0.512772,0.604864,0.701746,0.786448,0.841841},{0.313788,0.419946,0.511812,0.607382,0.695314,0.779678,0.837293}}; float tauArrayNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.304274,0.415066,0.505979,0.606817,0.709584,0.785154,0.839399},{0.304605,0.413666,0.507693,0.606070,0.707762,0.783781,0.838657},{0.307831,0.410343,0.511305,0.602437,0.703892,0.788212,0.841592},{0.310591,0.414093,0.512985,0.601204,0.702227,0.788886,0.840916}}; float constDay_MSG3_P[NBIAS_PROB][NBIAS] = {{-0.229574,-0.191087,-0.445556,-0.478594,-0.375241,-0.471156,-0.744909},{0.232056,-0.0306725,-0.353730,-0.417469,-0.565764,-0.697771,-0.620365},{-0.579094,-0.570791,-0.639485,-0.552455,-0.537639,-0.752596,-0.840732},{-0.115317,0.0461835,0.0461862,-0.0158708,-0.0729978,-0.100741,-0.119188}}; float slopeDay_MSG3_P[NBIAS_PROB][NBIAS] = {{0.442143,0.335897,-0.143489,-0.182639,-0.0920372,-0.163668,-0.260161},{0.760138,0.383800,0.161447,0.0705928,-0.0758799,-0.205991,-0.183782},{0.369690,-0.0565904,-0.0161108,-0.113601,-0.235330,-0.257873,-0.122688},{0.302185,0.192920,0.0137354,-0.0687899,-0.141896,-0.137585,-0.0795861}}; float constNight_MSG3_P[NBIAS_PROB][NBIAS] = {{-0.618860,-0.774110,-0.865725,-0.843158,-0.730721,-0.503794,-0.177249},{-0.460647,-0.470139,-0.539639,-0.534815,-0.403478,-0.234594,-0.0437389},{-0.331605,-0.279752,-0.241671,-0.238321,-0.107793,0.0643218,0.179566},{0.335418,0.107525,0.00716294,0.0611509,0.136486,0.240464,0.329549}}; float slopeNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.205630,0.0677421,0.00502583,0.0307784,0.0154946,0.0793491,0.176670},{0.269025,0.188822,0.105296,0.0783564,0.0609888,0.0778087,0.111832},{0.133473,0.211388,0.225469,0.130550,0.0997079,0.127676,0.145280},{0.501781,0.311292,0.155819,0.146592,0.126530,0.124623,0.114873}}; float stdevDay_MSG3_P[NBIAS_PROB][NBIAS] = {{1.19929,1.25186,1.12093,1.11477,1.10204,1.17468,1.29512},{0.833210,0.840721,0.839270,0.926651,1.04692,1.10575,1.17091},{1.30481,0.798897,0.962638,1.04807,1.01886,1.05145,1.05825},{0.844057,0.780791,0.765368,0.742130,0.733442,0.758978,0.795076}}; float stdevNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.907566,0.909181,0.913858,0.974854,0.971938,0.944861,0.879330},{0.968956,0.840370,0.809461,0.893143,0.872679,0.843599,0.861292},{0.876440,0.834279,0.807832,0.775907,0.712778,0.719318,0.732442},{0.823643,0.694670,0.653640,0.673136,0.621754,0.591061,0.629726}}; #else float probArrayDay_MSG3[NBIAS_PROB][NBIAS] = {{0.888547,0.888477,0.892456,0.891365,0.891988,0.895103,0.896105},{0.974117,0.974698,0.974066,0.975302,0.976344,0.975129,0.974579},{0.996220,0.996081,0.995747,0.995693,0.996075,0.996200,0.996062},{0.999983,0.999984,0.999982,0.999984,0.999984,0.999985,0.999988}}; float probArrayNight_MSG3[NBIAS_PROB][NBIAS] = {{0.889686,0.887512,0.890094,0.890914,0.892078,0.891430,0.888492},{0.974865,0.975656,0.975417,0.975776,0.975234,0.974593,0.974844},{0.995830,0.995879,0.995992,0.995953,0.995951,0.996078,0.996194},{0.999857,0.999857,0.999878,0.999893,0.999881,0.999853,0.999839}}; float tauArrayDay_MSG3_P[NBIAS_PROB][NBIAS] = {{0.324602,0.416618,0.507928,0.605330,0.706155,0.785592,0.838530},{0.306979,0.421951,0.503156,0.605900,0.700293,0.783621,0.836193},{0.317432,0.418441,0.510837,0.603036,0.701600,0.783315,0.835164},{0.318364,0.423997,0.514034,0.609710,0.703721,0.786207,0.841911}}; float tauArrayNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.309682,0.414954,0.505949,0.605670,0.707616,0.786718,0.836898},{0.308795,0.413410,0.507603,0.604706,0.708649,0.786161,0.837455},{0.306782,0.414150,0.508090,0.602326,0.702251,0.788569,0.840009},{0.310733,0.415762,0.513589,0.602053,0.697009,0.790520,0.843508}}; float constDay_MSG3_P[NBIAS_PROB][NBIAS] = {{-0.376744,-0.301169,-0.401568,-0.435039,-0.593038,-0.585056,-0.595208},{0.0482780,-0.125431,-0.402471,-0.640007,-0.720588,-0.707597,-0.693210},{-0.540416,-0.528290,-0.386481,-0.441960,-0.555032,-0.665315,-0.736960},{ 0.0164408,-0.0344198, 0.0000155, 0.0106857,-0.0179222, 0.0265621, 0.0840021}}; float slopeDay_MSG3_P[NBIAS_PROB][NBIAS] = {{0.506721, 0.205913, -0.199585,-0.0389948,-0.0848973, -0.130856, -0.125885},{1.00180, 0.650059, 0.381212,0.0933715,-0.142577,-0.235566,-0.239068},{0.291239, 0.225446, 0.188209, 0.0869040,-0.0578114, -0.172481, -0.167732},{0.323330, 0.158842, 0.0489850,-0.0678940, -0.135649, -0.129956, -0.116706}}; float constNight_MSG3_P[NBIAS_PROB][NBIAS] = {{-0.728557,-0.802580,-0.845953,-0.869644,-0.758978,-0.524964,-0.287198},{ -0.578636, -0.534382, -0.545936, -0.574475, -0.425940, -0.217384,-0.0754272},{-0.324537,-0.352071,-0.290668,-0.237832,-0.116301,0.0836247, 0.247517},{0.362035,0.193630,0.144031,0.139718,0.167273,0.290755,0.395634}}; float slopeNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.144721,0.0983967,0.0662757,0.0337547,0.0281491,0.0918672, 0.147579},{0.109641, 0.138527, 0.115475,0.0750003,0.0748850,0.0935038,0.0921714},{0.133076,0.172025,0.222471,0.177562,0.144870,0.139214,0.131063},{0.507777, 0.411701, 0.349673, 0.223007, 0.140696, 0.120895,0.0822841}}; float stdevDay_MSG3_P[NBIAS_PROB][NBIAS] = {{1.32742,1.25394,1.17877,1.26660,1.22669,1.15129,1.25907},{1.44941,1.18646,1.05226,1.23384,1.15502,1.13482,1.39846},{1.43436,1.39360,1.05844,1.16297,1.11688,1.07288,1.29480},{ 1.02060,0.972632,0.863572,0.762154,0.715097,0.727109,0.752480}}; float stdevNight_MSG3_P[NBIAS_PROB][NBIAS] = {{0.962255,0.942267,0.928625,0.964525,0.977554,0.980117,0.969031},{1.01054,0.925202,0.866742,0.924138,0.889855,0.846424,0.865206},{0.962126,0.866818,0.810073,0.767430,0.714048,0.712792,0.718088},{0.938827,0.846861,0.771563,0.663386,0.562015,0.564683,0.592477}}; #endif /* Routine to return the bias/std dev. */ void getBias( int satId, int dayTime, float dBTdSST, float tempDiff, float pClr, float *pBias, float *pStdev ) { int use_new_coeffs = 0; switch(satId){ case 72: /* GOES-9 */ use_new_coeffs = new_coeffs_goes9; break; case 74: /* GOES-10 */ use_new_coeffs = new_coeffs_goes10; break; case 76: /* GOES-11 */ use_new_coeffs = new_coeffs_goes11; break; case 78: /* GOES-12 */ use_new_coeffs = new_coeffs_goes12; break; case 180: /* GOES-13 */ use_new_coeffs = new_coeffs_goes13; break; case 184: /* GOES-15 */ use_new_coeffs = new_coeffs_goes15; break; case 84: /* MTSAT */ use_new_coeffs = new_coeffs_mtsat; break; case 85: /* MTSAT 2 */ use_new_coeffs = new_coeffs_mtsat2; break; case 52: /* MSG2 */ use_new_coeffs = new_coeffs_msg2; break; case 53: /* MSG3 */ use_new_coeffs = new_coeffs_msg3; break; default: fprintf(stderr,"Satellite ID = %d\n",satId); printError("L2P","Satellite ID not allowed (2)",NULL,NULL); } if( 1 == use_new_coeffs ) getBiasNew( satId, pClr, dayTime, dBTdSST, tempDiff, pBias, pStdev ); else getBiasOld( satId, dayTime, dBTdSST, tempDiff, pBias, pStdev ); } void getBiasNew( int satId, float pClr, int dayTime, float dBTdSST, float tempDiff, float *pBias, float *pStdev ) { int pos = 0; float tauVal = 0.; float *pTauArrayDay = NULL; float *pTauArrayNight = NULL; float *pConstDay = NULL; float *pSlopeDay = NULL; float *pConstNight = NULL; float *pSlopeNight = NULL; float *pStdevDay = NULL; float *pStdevNight = NULL; /* Determined which set of lookup tables to use */ /* select first from probability */ if( pClr >= 0.8 && pClr < 0.95 ) pos = 0; else if( pClr >= 0.95 && pClr < 0.99 ) pos = 1; else if( pClr >= 0.99 && pClr < 0.999 ) pos = 2; else if( pClr >= 0.999 ) pos = 3; else { /* Cannot calculate Bias/Stdev */ *pBias = -999.; *pStdev = -999.; return; } /* Set pointers to correct tables */ switch(satId){ case 84: /* MTSAT */ /* For MTSAT coefficients were done in tangent linear space */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MTSAT_P[pos][0]; pTauArrayNight = &tauArrayNight_MTSAT_P[pos][0]; pConstDay = &constDay_MTSAT_P[pos][0]; pSlopeDay = &slopeDay_MTSAT_P[pos][0]; pConstNight = &constNight_MTSAT_P[pos][0]; pSlopeNight = &slopeNight_MTSAT_P[pos][0]; pStdevDay = &stdevDay_MTSAT_P[pos][0]; pStdevNight = &stdevNight_MTSAT_P[pos][0]; break; case 85: /* MTSAT 2 */ /* For MTSAT coefficients were done in tangent linear space */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MTSAT2_P[pos][0]; pTauArrayNight = &tauArrayNight_MTSAT2_P[pos][0]; pConstDay = &constDay_MTSAT2_P[pos][0]; pSlopeDay = &slopeDay_MTSAT2_P[pos][0]; pConstNight = &constNight_MTSAT2_P[pos][0]; pSlopeNight = &slopeNight_MTSAT2_P[pos][0]; pStdevDay = &stdevDay_MTSAT2_P[pos][0]; pStdevNight = &stdevNight_MTSAT2_P[pos][0]; break; case 52: /* MSG2 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MSG2_P[pos][0]; pTauArrayNight = &tauArrayNight_MSG2_P[pos][0]; pConstDay = &constDay_MSG2_P[pos][0]; pSlopeDay = &slopeDay_MSG2_P[pos][0]; pConstNight = &constNight_MSG2_P[pos][0]; pSlopeNight = &slopeNight_MSG2_P[pos][0]; pStdevDay = &stdevDay_MSG2_P[pos][0]; pStdevNight = &stdevNight_MSG2_P[pos][0]; break; case 53: /* MSG3 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MSG3_P[pos][0]; pTauArrayNight = &tauArrayNight_MSG3_P[pos][0]; pConstDay = &constDay_MSG3_P[pos][0]; pSlopeDay = &slopeDay_MSG3_P[pos][0]; pConstNight = &constNight_MSG3_P[pos][0]; pSlopeNight = &slopeNight_MSG3_P[pos][0]; pStdevDay = &stdevDay_MSG3_P[pos][0]; pStdevNight = &stdevNight_MSG3_P[pos][0]; break; case 184: /* GOES15 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G15_P[pos][0]; pTauArrayNight = &tauArrayNight_G15_P[pos][0]; pConstDay = &constDay_G15_P[pos][0]; pSlopeDay = &slopeDay_G15_P[pos][0]; pConstNight = &constNight_G15_P[pos][0]; pSlopeNight = &slopeNight_G15_P[pos][0]; pStdevDay = &stdevDay_G15_P[pos][0]; pStdevNight = &stdevNight_G15_P[pos][0]; break; default: fprintf(stderr,"Satellite ID = %d\n",satId); printError("L2P","Satellite ID not allowed (3)",NULL,NULL); } getBiasValue( dayTime, tauVal, tempDiff, pTauArrayDay, pTauArrayNight, pConstDay, pSlopeDay, pConstNight, pSlopeNight, pStdevDay, pStdevNight, pBias, pStdev ); } void getBiasOld( int satId, int dayTime, float dBTdSST, float tempDiff, float *pBias, float *pStdev ) { float tauVal = 0.; float *pTauArrayDay = NULL; float *pTauArrayNight = NULL; float *pConstDay = NULL; float *pSlopeDay = NULL; float *pConstNight = NULL; float *pSlopeNight = NULL; float *pStdevDay = NULL; float *pStdevNight = NULL; /* Determined which set of lookup tables to use */ /* Set pointers to correct tables */ switch(satId){ case 72: /* GOES-9 */ /* transmittence */ tauVal = 1.03026 * dBTdSST - 0.0327601; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G9[0]; pTauArrayNight = &tauArrayNight_G9[0]; pConstDay = &constDay_G9[0]; pSlopeDay = &slopeDay_G9[0]; pConstNight = &constNight_G9[0]; pSlopeNight = &slopeNight_G9[0]; pStdevDay = &stdevDay_G9[0]; pStdevNight = &stdevNight_G9[0]; break; case 74: /* GOES-10 */ /* transmittence */ tauVal = 1.04062 * dBTdSST - 0.0383824; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G10[0]; pTauArrayNight = &tauArrayNight_G10[0]; pConstDay = &constDay_G10[0]; pSlopeDay = &slopeDay_G10[0]; pConstNight = &constNight_G10[0]; pSlopeNight = &slopeNight_G10[0]; pStdevDay = &stdevDay_G10[0]; pStdevNight = &stdevNight_G10[0]; break; case 76: /* GOES-11 */ /* transmittence */ tauVal = 1.04062 * dBTdSST - 0.0383824; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G11[0]; pTauArrayNight = &tauArrayNight_G11[0]; pConstDay = &constDay_G11[0]; pSlopeDay = &slopeDay_G11[0]; pConstNight = &constNight_G11[0]; pSlopeNight = &slopeNight_G11[0]; pStdevDay = &stdevDay_G11[0]; pStdevNight = &stdevNight_G11[0]; break; case 78: /* GOES-12 */ /* transmittence */ tauVal = 1.03854 * dBTdSST - 0.0361785; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G12[0]; pTauArrayNight = &tauArrayNight_G12[0]; pConstDay = &constDay_G12[0]; pSlopeDay = &slopeDay_G12[0]; pConstNight = &constNight_G12[0]; pSlopeNight = &slopeNight_G12[0]; pStdevDay = &stdevDay_G12[0]; pStdevNight = &stdevNight_G12[0]; break; case 180: /* GOES-13 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_G13[0]; pTauArrayNight = &tauArrayNight_G13[0]; pConstDay = &constDay_G13[0]; pSlopeDay = &slopeDay_G13[0]; pConstNight = &constNight_G13[0]; pSlopeNight = &slopeNight_G13[0]; pStdevDay = &stdevDay_G13[0]; pStdevNight = &stdevNight_G13[0]; break; case 84: /* MTSAT */ /* For MTSAT coefficients were done in tangent linear space */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MTSAT[0]; pTauArrayNight = &tauArrayNight_MTSAT[0]; pConstDay = &constDay_MTSAT[0]; pSlopeDay = &slopeDay_MTSAT[0]; pConstNight = &constNight_MTSAT[0]; pSlopeNight = &slopeNight_MTSAT[0]; pStdevDay = &stdevDay_MTSAT[0]; pStdevNight = &stdevNight_MTSAT[0]; break; case 52: /* MSG2 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MSG2[0]; pTauArrayNight = &tauArrayNight_MSG2[0]; pConstDay = &constDay_MSG2[0]; pSlopeDay = &slopeDay_MSG2[0]; pConstNight = &constNight_MSG2[0]; pSlopeNight = &slopeNight_MSG2[0]; pStdevDay = &stdevDay_MSG2[0]; pStdevNight = &stdevNight_MSG2[0]; break; case 53: /* MSG3 */ /* transmittence */ tauVal = dBTdSST; /* Pointers to the lookup tables */ pTauArrayDay = &tauArrayDay_MSG3[0]; pTauArrayNight = &tauArrayNight_MSG3[0]; pConstDay = &constDay_MSG3[0]; pSlopeDay = &slopeDay_MSG3[0]; pConstNight = &constNight_MSG3[0]; pSlopeNight = &slopeNight_MSG3[0]; pStdevDay = &stdevDay_MSG3[0]; pStdevNight = &stdevNight_MSG3[0]; break; default: fprintf(stderr,"Satellite ID = %d\n",satId); printError("L2P","Satellite ID not allowed (2)",NULL,NULL); } getBiasValue( dayTime, tauVal, tempDiff, pTauArrayDay, pTauArrayNight, pConstDay, pSlopeDay, pConstNight, pSlopeNight, pStdevDay, pStdevNight, pBias, pStdev ); } void getBiasValue( int dayTime, float tauVal, float tempDiff, float *pTauArrayDay, float *pTauArrayNight, float *pConstDay, float *pSlopeDay, float *pConstNight, float *pSlopeNight, float *pStdevDay, float *pStdevNight, float *pBias, float *pStdev ) { int i = 0; int pos = 0; float constant = 0.; float slope = 0.; float stdev = 0.; /* If no values assigned - return setting error flags */ if( -1 == tauVal ){ *pBias = -999.; *pStdev = -999.; return; } /* Algorithm depends on day/night */ switch(dayTime){ case 0: /* Night time */ /* Find out where our tau value fits in with lookup tables */ pos = -1; for(i=0;i windU || 25 < windU ){ windU = 25.; } windV = *(pWindSpeedV++); if( 0 > windV || 25 < windV ){ windV = 25.; } *(pWindSpeedPtr++) = sqrt(windU*windU+windV*windV); } } /* Allocate output array */ if( NULL == (pDiurnal1 = (float *) calloc(model1.nDims1*model1.nDims2, sizeof(float))) ){ printError("L2P","Error allocating pWindSpeedArr",NULL,NULL); } if( FAIL == diurnal_gentemann( Year, Month, Day, newHours, mins, pWindSpeedArr,model1.nDims2,model1.nDims1, pDiurnal1) ){ printError("L2P","Diurnal model failed (model1)",NULL,NULL); } /* Set any bad values to 0 diurnal warming until better interpoilation software which can deal * with bad corners */ for(i=0;i windU || 25 < windU ){ windU = 25.; } windV = *(pWindSpeedV++); if( 0 > windV || 25 < windV ){ windV = 25.; } *(pWindSpeedPtr++) = sqrt(windU*windU+windV*windV); } } /* New UTC in seconds - add three hours */ newHours = newHours + 3; /* If gone over the day - get new day/month/year */ if( 24 <= newHours ){ /* This routine is in read_grib.c */ getDayNoFromMonthDay( Year, Month, Day, &dayNo ); /* Have we crossed over to a new year? */ if( ((0 == Year%4 && 0 != Year%100) || 0 == Year%400) && 366 == dayNo){ /* leap year */ newYear = Year + 1; newMonth = 1; newDay = 1; } else if( 365 == dayNo ){ /* Not a leap year */ newYear = Year + 1; newMonth = 1; newDay = 1; } else { /* Just get new day/month/year */ dayNo += 1; getMonthFromDayNo(Year,dayNo,&newMonth,&newDay); newYear = Year; } /* Get correct utc for this day */ newHours =- 24; } else { /* Not crossed over the day boundary */ newYear = Year; newMonth = Month; newDay = Day; } /* Allocate output array */ if( NULL == (pDiurnal2 = (float *) calloc(model1.nDims1*model1.nDims2, sizeof(float))) ){ printError("L2P","Error allocating pWindSpeedArr",NULL,NULL); } if( FAIL == diurnal_gentemann( newYear, newMonth, newDay, newHours, mins, pWindSpeedArr,model2.nDims2, model2.nDims1, pDiurnal2 ) ){ printError("L2P","Diurnal model failed (model2)",NULL,NULL); } /* Set any bad values to 0 diurnal warming until better interpoilation software which can deal * with bad corners */ for(i=0;i 70 */ /* Note that this is nearest neighbour */ /* Mapped aerosol to 1 degree grid */ if( 1 != aerosol.stepLongitude || 1 != aerosol.stepLatitude ) printError("L2P","Aerosol file not on 1 degree grid",NULL,NULL); offset1 = (int) (longitude+0.5); if( 360 == offset1 ) offset1 = 0; if( latitude >= 0 ) offset2 = (int) ((90+latitude)+0.5); else offset2 = (int) ((90+latitude)-0.5); *pAerOptDepth = *(aerosol.pAerosol+offset1+offset2*aerosol.nDims1); *pAerOptDepthTime = *(aerosol.pAerosolDate+offset1+offset2*aerosol.nDims1); /* Get previous SST */ getCellPosition( longitude, latitude, previousSst.nDims1, previousSst.nDims2, previousSst.minLongitude, previousSst.maxLongitude, previousSst.stepLongitude, previousSst.minLatitude, previousSst.maxLatitude, previousSst.stepLatitude, &offset1, &offset2, &offset3, &offset4, &minLongitude, &maxLongitude, &minLatitude, &maxLatitude); /* No need for time interpolation - only a single value */ if( -1 == offset2 && -1 == offset3 && -1 == offset4 ) if( -1 == offset1 ) *pPreviousSST = -999; else *pPreviousSST = *(previousSst.pSST+offset1); else *pPreviousSST = bilinearInterpolation( minLongitude, *(previousSst.pSST+offset1), maxLongitude, *(previousSst.pSST+offset2), minLatitude, *(previousSst.pSST+offset3), maxLatitude, *(previousSst.pSST+offset4), longitude, latitude ); } /* Bi-linear interpolation routine */ float bilinearInterpolation( float x1, float val1, float x2, float val2, float y1, float val3, float y2, float val4, float xin, float yin ) { float t = 0.; float u = 0.; float value = 0.; /* Do bi-linear interpolation */ t = (xin - x1)/(x2 - x1); u = (yin - y1)/(y2 - y1); value = (1-t)*(1-u)*val1 + t*(1-u)*val2 + t*u*val4 + (1-t)*u*val3; return value; } /* Get the indexes of a cell and the min/max values so we can do the * bilinear interpolation. If the position is on a grid point, the * index is returned of offset1 and the other offsets are set to -1 */ void getCellPosition( float newLongitude, float newLatitude, int nDims1, int nDims2, float_fp_kind minLongitude, float_fp_kind maxLongitude, float_fp_kind stepLongitude, float_fp_kind minLatitude, float_fp_kind maxLatitude, float_fp_kind stepLatitude, int *pOffset1, int *pOffset2, int *pOffset3, int *pOffset4, float *pMinLongitudePos, float *pMaxLongitudePos, float *pMinLatitudePos, float *pMaxLatitudePos ) { int cellXposition = 0; int cellYposition = 0; float fltCellXPosition = 0.; float fltCellYPosition = 0.; float checkLongitude = 0.; float checkLatitude = 0.; /* First check that we can access data */ /* Needed for aerosol data */ /* Note that for model is backward due to way data stored */ /* Check is we're looking at a model grib file */ if( minLatitude < 0 && maxLatitude > 0 ){ if( newLatitude < minLatitude || newLatitude > maxLatitude ){ *pOffset1 = -1; *pOffset2 = -1; *pOffset3 = -1; *pOffset4 = -1; return; } } /* Check to see if we need to interpolate or not */ fltCellXPosition = ((newLongitude - minLongitude)/stepLongitude); fltCellYPosition = ((newLatitude - minLatitude)/stepLatitude); cellXposition = (int) fltCellXPosition; cellYposition = (int) fltCellYPosition; /* Do check due to some flt precision issues */ /* Not a number above the first node is OK, below is not */ if( (float) cellXposition == fltCellXPosition ){ checkLongitude = fltCellXPosition*stepLongitude + minLongitude; if( newLongitude < checkLongitude ){ fltCellXPosition -= 1.; cellXposition -= 1; } } if( (float) cellYposition == fltCellYPosition ){ checkLatitude = fltCellYPosition*stepLatitude + minLatitude; if( newLatitude < checkLatitude ){ fltCellYPosition -= 1.; cellYposition -= 1; } } /* Check if we need to interpolate at all */ if( (float) cellXposition == fltCellXPosition && (float) cellYposition == fltCellYPosition ){ /* We do not need to interpolate */ *pOffset1 = cellXposition + nDims1 * cellYposition; *pOffset2 = -1; *pOffset3 = -1; *pOffset4 = -1; } else { /* Note have to take into account going > 359 to wrap around correctly */ if( newLongitude >= maxLongitude || newLongitude < minLongitude){ /* Have to use last and first entries to cross the zero line */ *pOffset1 = nDims1-1 + nDims1 * cellYposition; *pOffset2 = nDims1 * cellYposition; if( nDims2 == (cellYposition+1) ){ *pOffset3 = nDims1-1; *pOffset4 = 0; } else { *pOffset3 = nDims1-1 + nDims1*(cellYposition+1); *pOffset4 = nDims1*(cellYposition+1); } *pMinLongitudePos = maxLongitude; *pMaxLongitudePos = maxLongitude+stepLongitude; /* Check the *pOffsets */ if( 0 > *pOffset1 || nDims1*nDims2 <= *pOffset1 ) printError( "L2P", "cell position 1 out of range(1)", NULL,NULL); if( 0 > *pOffset2 || nDims1*nDims2 <= *pOffset2 ) printError( "L2P", "cell position 2 out of range(1)", NULL,NULL); if( 0 > *pOffset3 || nDims1*nDims2 <= *pOffset3 ) printError( "L2P", "cell position 3 out of range(1)", NULL,NULL); if( 0 > *pOffset4 || nDims1*nDims2 <= *pOffset4 ) printError( "L2P", "cell position 4 out of range(1)", NULL,NULL); } else { /* Get *pOffsets to the square we are at */ /* Note this wont work at the poles - but no sea at the poles (yet) */ *pOffset1 = cellXposition + nDims1 * cellYposition; *pOffset2 = cellXposition+1 + nDims1 * cellYposition; /* Deal with wrap over the poles */ if( nDims2 == (cellYposition+1) ){ *pOffset3 = cellXposition; *pOffset4 = cellXposition+1; } else { *pOffset3 = cellXposition+nDims1*(cellYposition+1); *pOffset4 = cellXposition+1+nDims1*(cellYposition+1); } /* Check the *pOffsets */ if( 0 > *pOffset1 || nDims1*nDims2 <= *pOffset1 ){ printError( "L2P", "cell position 1 out of range(2)", NULL,NULL); } if( 0 > *pOffset2 || nDims1*nDims2 <= *pOffset2 ){ printError( "L2P", "cell position 2 out of range(2)", NULL,NULL); } if( 0 > *pOffset3 || nDims1*nDims2 <= *pOffset3 ){ printError( "L2P", "cell position 3 out of range(2)", NULL,NULL); } if( 0 > *pOffset4 || nDims1*nDims2 <= *pOffset4 ){ printError( "L2P", "cell position 4 out of range(2)", NULL,NULL); } /* Get positions of cell points - Longitude */ *pMinLongitudePos = minLongitude+cellXposition*stepLongitude; *pMaxLongitudePos = minLongitude+(cellXposition+1)* stepLongitude; } *pMinLatitudePos = minLatitude+cellYposition*stepLatitude; *pMaxLatitudePos = minLatitude+(cellYposition+1)*stepLatitude; } } /* Get the time of observations in seconds from 1 Jan 1980 */ /* changing from 2444606.5 to 2444605.5, per reality */ #define START_JULDAY 2444605.5 int getTimeValue( int Year, int Month, int Day, float utc ) { int YYYY = 0; double julDay = 0; double totDays = 0; if( Year < 100 ){ if( Year < 50 ) YYYY = 2000+Year; else YYYY = 1900+Year; } else YYYY = Year; julDay = ((double) getJulianDay( Month, Day, YYYY )) - 0.5; totDays = (julDay - START_JULDAY)*86400. + utc; return (int) totDays; } /* Get the Julian date to be used in getTimeValue */ #define IGREG (15+31*(10+12*1582)) int getJulianDay( int mm, int id, int iyyy ) { int jy = 0; int jm = 0; int ja = 0; int julday = 0; jy = iyyy; if( mm > 2 ) jm = mm+1; else { jy = jy-1; jm = mm+13; } julday = ((int) (365.25*jy)) + ((int) (30.6001*jm)) + id + 1720995; if( id+31*(mm+12*iyyy) >= IGREG ){ ja = ((int) (0.01*jy)); julday += 2 - ja + ((int) (0.25*ja)); } return julday; } /***************************************************************************/ /* Note that this code is not used at the moment - here as a place setted * just in case it is needed */ #ifdef RUN_CRTM_MODEL /* Run the CRTM model */ #define NLEVELS 26 #define NCHANNELS 3 void runCRTM_Model( int satId, struct modelVals *pModel, float minLong, float maxLong, float minLat, float maxLat, float maxZenith ) { int i = 0; int longCntr = 0; int latCntr = 0; int offset = 0; int nlevels = NLEVELS; int nchannels = NCHANNELS; int channelArray[NCHANNELS] = {1,3,4}; float_fp_kind zenithAngle = 0.; float_fp_kind solarAngle = 0.; float_fp_kind height[NLEVELS]; float_fp_kind ozone[NLEVELS]; float_fp_kind surfacePressure = 0.; float_fp_kind RH[NLEVELS]; float_fp_kind surfaceTemperature = 0.; float_fp_kind temperatures[NLEVELS]; float_fp_kind windSpeed = 0.; float_fp_kind BT[NCHANNELS]; float_fp_kind Rad[NCHANNELS]; float_fp_kind Tau[NCHANNELS]; float_fp_kind dBTdSST[NCHANNELS]; float_fp_kind pressure[NLEVELS]; float_fp_kind latitude = 0.; float_fp_kind longitude = 0.; float_fp_kind Hours = 0.; int errorCode = 0; float_fp_kind val1 = 0.; float_fp_kind val2 = 0.; float_fp_kind sateliteLat = 0.; float_fp_kind sateliteLong = 0.; float_fp_kind sateliteHeight = 0.; /* Setup for given GOES satelite */ switch( satId ){ case 72: sateliteLong = 154.5; sateliteLat = 0.; sateliteHeight = 42164.; break; case 74: sateliteLong = 360.-135.; sateliteLat = 0.; sateliteHeight = 42164.; break; case 78: sateliteLong = 360.-75.; sateliteLat = 0.; sateliteHeight = 42164.; break; case 180: sateliteLong = 360.-75.; sateliteLat = 0.; sateliteHeight = 42164.; break; default: printError("L2P","satID out of range",NULL,NULL); } /* Loop round arrays */ for(latCntr=0;latCntrnDims2;latCntr++){ latitude = pModel->minLatitude + latCntr*pModel->stepLatitude; for(longCntr=0;longCntrnDims1;longCntr++){ offset = longCntr + latCntr*pModel->nDims1; longitude = pModel->minLongitude + longCntr*pModel->stepLongitude; /* Get Angles */ zenithAngle = getZenithAnglePos( sateliteLong, sateliteLat, latitude, longitude, sateliteHeight ); #if 0 printf("%f %f %f %f %f = %f\n",sateliteLong, sateliteLat, latitude, longitude, sateliteHeight,zenithAngle); #endif #if 0 if( maxZenith+2 >= zenithAngle ){ #endif Hours = (float_fp_kind ) pModel->Hour; getSunAngle( pModel->Year, pModel->Month, pModel->Day, Hours, longitude, latitude, &solarAngle ); /* Copy across data for float_fp_kind arrays */ for(i=0;ippHGT+i)+offset)); temperatures[i] = (float_fp_kind) (*(*(pModel->ppTemperature+i)+offset)); if( NULL != *(pModel->ppOzone+i) ) ozone[i] = (float_fp_kind) (*(*(pModel->ppOzone+i)+offset)); else ozone[i] = 0.; if( NULL != *(pModel->ppRelativeHumidity+i) ){ RH[i] = (float_fp_kind) (*(*(pModel->ppRelativeHumidity+i)+offset)); } else RH[i] = 0.; pressure[i] = (float_fp_kind) *(pModel->pPressure+i); } surfacePressure = (float_fp_kind) *(pModel->pSurfacePressure+offset); surfaceTemperature = (float_fp_kind) *(pModel->pSurfaceTemperature+offset); val1 = *(pModel->pWindSpeedU+offset); val2 = *(pModel->pWindSpeedV+offset); windSpeed = sqrt(val1*val1 + val2*val2); #if 0 printf("===============================\n"); for(i=0;ipTau4+offset) = Tau[1]; #if 0 printf("Tau val : %f\n",Tau[1]); #endif #if 0 } else *(pModel->pTau4+offset) = -999; #endif } } } /* Copy one model to another */ void copyModel( struct modelVals *pModel1, struct modelVals *pModel2 ) { int i = 0; pModel1->Year = pModel2->Year; pModel1->Month = pModel2->Month; pModel1->Day = pModel2->Day; pModel1->Hour = pModel2->Hour; pModel1->nDims1 = pModel2->nDims1; pModel1->nDims2 = pModel2->nDims2; pModel1->minLatitude = pModel2->minLatitude; pModel1->maxLatitude = pModel2->maxLatitude; pModel1->stepLatitude = pModel2->stepLatitude; pModel1->minLongitude = pModel2->minLongitude; pModel1->maxLongitude = pModel2->maxLongitude; pModel1->stepLongitude = pModel2->stepLongitude; pModel1->nLevels = pModel2->nLevels; /* Arrays */ free(pModel1->pPressure); pModel1->pPressure = pModel2->pPressure; pModel2->pPressure = NULL; free(pModel1->pSurfacePressure); pModel1->pSurfacePressure = pModel2->pSurfacePressure; pModel2->pSurfacePressure = NULL; free(pModel1->pSurfaceTemperature); pModel1->pSurfaceTemperature = pModel2->pSurfaceTemperature; pModel2->pSurfaceTemperature = NULL; free(pModel1->pTemperature2m); pModel1->pTemperature2m = pModel2->pTemperature2m; pModel2->pTemperature2m = NULL; free(pModel1->pLand); pModel1->pLand = pModel2->pLand; pModel2->pLand = NULL; free(pModel1->pIce); pModel1->pIce = pModel2->pIce; pModel2->pIce = NULL; free(pModel1->pWindSpeedU); pModel1->pWindSpeedU = pModel2->pWindSpeedU; pModel2->pWindSpeedU = NULL; free(pModel1->pWindSpeedV); pModel1->pWindSpeedV = pModel2->pWindSpeedV; pModel2->pWindSpeedV = NULL; free(pModel1->ppHGT); /* 2d models */ pModel1->ppHGT = pModel2->ppHGT; pModel2->ppHGT = NULL; free(pModel1->ppTemperature); pModel1->ppTemperature = pModel2->ppTemperature; pModel2->ppTemperature = NULL; free(pModel1->ppOzone); pModel1->ppOzone = pModel2->ppOzone; pModel2->ppOzone = NULL; free(pModel1->ppRelativeHumidity); pModel1->ppRelativeHumidity = pModel2->ppRelativeHumidity; pModel2->ppRelativeHumidity = NULL; for(i=0;inLevels;i++){ free(*(pModel1->ppHGT+i)); *(pModel1->ppHGT+i) = *(pModel2->ppHGT+i); *(pModel2->ppHGT+i) = NULL; free(*(pModel1->ppTemperature+i)); *(pModel1->ppTemperature+i) = *(pModel2->ppTemperature+i); *(pModel2->ppTemperature+i) = NULL; free(*(pModel1->ppOzone+i)); *(pModel1->ppOzone+i) = *(pModel2->ppOzone+i); *(pModel2->ppOzone+i) = NULL; free(*(pModel1->ppRelativeHumidity+i)); *(pModel1->ppRelativeHumidity+i) = *(pModel2->ppRelativeHumidity+i); *(pModel2->ppRelativeHumidity+i) = NULL; } pModel2->Filled = 0; } void copyModel2( struct modelVals2 *pModel1, struct modelVals2 *pModel2 ) { int i = 0; pModel1->Year = pModel2->Year; pModel1->Month = pModel2->Month; pModel1->Day = pModel2->Day; pModel1->Hour = pModel2->Hour; pModel1->nDims1 = pModel2->nDims1; pModel1->nDims2 = pModel2->nDims2; pModel1->minLatitude = pModel2->minLatitude; pModel1->maxLatitude = pModel2->maxLatitude; pModel1->stepLatitude = pModel2->stepLatitude; pModel1->minLongitude = pModel2->minLongitude; pModel1->maxLongitude = pModel2->maxLongitude; pModel1->stepLongitude = pModel2->stepLongitude; /* Arrays */ free(pModel1->pSurfaceSolarIrradiance); pModel1->pSurfaceSolarIrradiance = pModel2->pSurfaceSolarIrradiance; pModel2->pSurfaceSolarIrradiance = NULL; pModel2->Filled = 0; } #define DEG_TO_RADIANS (M_PI/180.) /* Earth radius (km) */ #define EARTH_RADIUS 6378. float_fp_kind getZenithAnglePos( float_fp_kind sateliteLong, float_fp_kind sateliteLat, float_fp_kind longitude, float_fp_kind latitude, float_fp_kind sateliteHeight ) { float_fp_kind c_angle = 0.; float_fp_kind scan_angle = 0.; float_fp_kind angle1 = 0.; float_fp_kind angle2 = 0.; float_fp_kind angle3 = 0.; if( -1. == sateliteHeight ) printError("L2P","-1 == sateliteHeight - should not be here", NULL,NULL); /* First get the zenith angle */ /* Convert to radians and from the north pole */ angle1 = (90.-sateliteLat)*DEG_TO_RADIANS; angle2 = (90.-latitude)*DEG_TO_RADIANS; angle3 = (longitude - sateliteLong)*DEG_TO_RADIANS; /* Get angular distance between footpoint position and satelite */ c_angle = cos(angle1)*cos(angle2) + sin(angle1)*sin(angle2)*cos(angle3); c_angle = acos(c_angle); scan_angle = EARTH_RADIUS*sin(c_angle)/ ((sateliteHeight+EARTH_RADIUS)-EARTH_RADIUS*cos(c_angle)); scan_angle = atan(scan_angle); /* return angle model needs in degrees */ return (c_angle + scan_angle)/DEG_TO_RADIANS; } /* Taken from the web page http://www.astro.uio.no/~bgranslo/aares/calculate.html */ /* Hours in UT */ void getSunAngle( int Year, int Month, int Day, float_fp_kind Hours, float_fp_kind longitude, float_fp_kind latitude, float_fp_kind *pSunAngle ) { int numPI = 0; int nHours = 0; float_fp_kind a = 0.; float_fp_kind b = 0.; float_fp_kind c = 0.; float_fp_kind e = 0.; float_fp_kind f = 0.; float_fp_kind yr = 0.; float_fp_kind mnth = 0.; float_fp_kind numberOfDays = 0.; float_fp_kind numberOfDaysUT = 0.; float_fp_kind numberJulianCenturies = 0.; float_fp_kind numberCenturies = 0.; float_fp_kind siderialTime0Greenwich = 0.; float_fp_kind siderialTime = 0.; float_fp_kind localSiderialTime = 0.; float_fp_kind solarMeanLong = 0.; float_fp_kind solarMeanAnomoly = 0.; float_fp_kind solarCentre = 0.; float_fp_kind solarTrueLong = 0.; float_fp_kind solarObliquity = 0.; float_fp_kind sunRA = 0.; float_fp_kind sunDec = 0.; float_fp_kind hourAngle = 0.; float_fp_kind altitude = 0.; /* Number of days and centuries from epoch J2000.0 */ /* Number of days from J2000.0 Jan 1st */ /* Get JD for this date */ if( 1 == Month || 2 == Month ){ yr = Year - 1; mnth = Month+12; } else { yr = Year; mnth = Month; } a = ((int) (yr/100.))*1.; b = ((int) (a/4.))*1.; c = 2.-a+b; e = ((int) (365.25*(yr+4716)))*1.; f = ((int) (30.6001*(mnth+1.)))*1.; /* This is number of days since 2000 1 jan 0h UT */ numberOfDays = c + e + f + Day - 1524.5 - 2451545.5; numberJulianCenturies = numberOfDays / 36525.; /* Siderial time */ /* Taken from http://aa.usno.navy.mil/faq/docs/GAST.html */ siderialTime0Greenwich = 6.697374558 + 0.06570982441908*numberOfDays + 1.00273790935*Hours; nHours = (int) (siderialTime0Greenwich/24.); /* Get to range 0 - 24 */ siderialTime0Greenwich -= nHours*24.; localSiderialTime = siderialTime0Greenwich + longitude/15.; if( 0 > localSiderialTime ) localSiderialTime += 24.; else if( 24. < localSiderialTime ) localSiderialTime -= 24.; /* Solar coordinates */ numberOfDaysUT = numberOfDays + Hours/24.; numberCenturies = numberOfDaysUT / 36525.; solarMeanLong = 280.466 + 36000.770 * numberCenturies; solarMeanAnomoly = 357.529 + 35999.050 * numberCenturies; solarCentre = (1.915 - 0.005*numberCenturies)* sin(solarMeanAnomoly*DEG_TO_RADIANS) + 0.020*sin(2*solarMeanAnomoly*DEG_TO_RADIANS); /* Get true ecliptic longitude */ solarTrueLong = solarMeanLong + solarCentre; /* Make sure we're in 0 - 360. */ numPI = (int) (solarTrueLong / 360.); if( 0 <= solarTrueLong ){ solarTrueLong -= numPI * 360.; } else { solarTrueLong = (numPI+1) * 360. - solarTrueLong; } solarTrueLong *= DEG_TO_RADIANS; solarObliquity = (23.439 - 0.013 * numberCenturies)*DEG_TO_RADIANS; sunRA = tan(solarTrueLong) * cos(solarObliquity); /* Make sure we're in the right quadrant */ if( 0. <= solarTrueLong && M_PI_2 >= solarTrueLong ){ sunRA = atan(sunRA); } else if ( M_PI_2 < solarTrueLong && 1.5*M_PI >= solarTrueLong ){ sunRA = M_PI + atan(sunRA); } else if ( 1.5*M_PI <= solarTrueLong && 2.*M_PI >= solarTrueLong ){ sunRA = 2*M_PI + atan(sunRA); } sunDec = sin(sunRA) * sin(solarObliquity); hourAngle = localSiderialTime*15. - sunRA/DEG_TO_RADIANS; altitude = sin(latitude*DEG_TO_RADIANS)*sunDec + cos(latitude*DEG_TO_RADIANS)*cos(asin(sunDec))* cos(hourAngle*DEG_TO_RADIANS); /* Sun angle in degrees */ *pSunAngle = asin(altitude)/DEG_TO_RADIANS; } #endif