/*============================================================================
 *
 * Name        : LST_Retrieval.cpp
 * Author      : Li Fang
 * Version     : v1.0
 * Description : Generate GOES Imager Continental U.S. (CONUS) and Full Disk (FD)
 *               Land Surface Temperature (LST) products.
 *               LST retrievals will be extracted through Dual Window algorithm
 *               (Sun and Pinker, 2004: J. Applied Meteorology)
 *               The output product include LST retrievals, observation time,
 *               pixel location, data quality flags and metadata
 *
 * Li Fang,      lfang1@gmu.edu  George Mason University
 * Yuling Liu,   yliu16@gmu.edu  George Mason University
 * Donglian Sun, dsun@gmu.edu    George Mason University
 *
 *Created on: Aug 19, 2010
 *
 *---------------------------------------------------------------------------
 * HDF LIBRARY
 * Hierarchical Data Format (HDF) Software Library and Utilities
 * Copyright 2006-2011 by The HDF Group.
 * NCSA Hierarchical Data Format (HDF) Software Library and Utilities
 * Copyright 1988-2006 by the Board of Trustees of the University of Illinois
 * All rights reserved.
 *---------------------------------------------------------------------------
 *NetCDF LIBRARY
 *Copyright 1993-2004 University Corporation for Atmospheric Research/Unidata
 * Portions of this software were developed by the Unidata Program
 * at the University Corporation for Atmospheric Research.
 *---------------------------------------------------------------------------
 *
 *============================================================================
 */

#include "LST_Retrieval.h"
#include <stdlib.h>
#include <stdio.h>
#include <netcdfcpp.h>
#include "ImageSampler.h"
#include <mfhdf.h>
#include "Para_define.h"
#include "grib2include.h"
#include "execinfo.h"
#include <limits.h>
#include <math.h>
#include <memory.h>
#include "GOES_LST_MAIN.h"

using namespace std;

LST_Retrieval::LST_Retrieval()
{
   /*Variables initialization*/
   root_allday = NULL; /*Regression model for day time */
   root_correction = NULL; /*Regression model for night time */
   ImgHeight = 0; /*Initiation of image height*/
   ImgWidth = 0; /*Initiation of image width*/
   img_date = 0; /*Initiation of image date*/
   img_time = 0; /*Initiation of image time*/
   GridWidth = 0;
   GridHeight = 0;

}

LST_Retrieval::~LST_Retrieval()
{
   if (root_correction != NULL)
   {
      delete root_correction;
   }
   if (root_allday != NULL)
   {
      delete root_allday;
   }
}

/*Read coefficients of regression model*/
Condition* LST_Retrieval::parser(string filename, vector<string> columns)
{
   ifstream file;
   file.open(filename.c_str());
   if (file.fail())
   {
      param::instance << getTimeString() << "File " << filename
            << " could not be opened." << endl;
      return NULL; /* error return for file open failure*/
   }
   file.seekg(0, ios_base::end);
   long flLength = file.tellg();
   file.seekg(0, ios_base::beg);
   char *temp = new char[flLength + 1];
   file.read(temp, flLength);
   temp[flLength] = 0;
   string input_coef(temp);
   delete[] temp;
   file.close();
   /* success file open/read return */
   return ConditionParser(input_coef, columns).analysis();
}

/*Read in and preprocess data for GOES 13 data from GSIP products*/
/* Modified by Yuling Liu for GSIP V3 process */
bool LST_Retrieval::GetData()
{

   /*get data from GSIP product in HDF format*/

   string dim_name;

   /*Diagnostic code begin*/
   /*test_file to output internal binary arrays*/
   ofstream test_file;
   string test_file_name;
   /*Diagnostic code end*/

   /* Get data from GSIP L1 product:
    brightness temperature at 3.9_micron
    brightness temperature at 11_micron
    pixel latitude
    pixel longitude
    sensor zenith
    solar zenith
    surface type
    */

   dim_name = param::instance.DIM_NAME_L1_3_9_BT;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, GOES_TP2_Aft_Cal);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   GetImgInfor();
   SetOutputFilename();
   sprintf(suffix_str, "_%d_%04d", img_date, img_time);
   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal brightness temperature 3.9*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_BT_39_"
            + suffix_str;
      ;
      test_file.open(test_file_name.c_str());
      GOES_TP2_Aft_Cal.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_L1_11_BT;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, GOES_TP4_Aft_Cal);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal brightness temperature 11*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_BT_11_"
            + suffix_str;
      test_file.open(test_file_name.c_str());
      GOES_TP4_Aft_Cal.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_PIXEL_LATITUDE;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, lat_degree);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal latitude*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_latitude_"
            + suffix_str;
      //test_file_name = "UT_latitude";
      test_file.open(test_file_name.c_str());
      lat_degree.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_PIXEL_LONGITUDE;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, lon_degree);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal longitude*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_longitude_"
            + suffix_str;
      //test_file_name = "UT_longitude";
      test_file.open(test_file_name.c_str());
      lon_degree.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_PIXEL_SENSOR_ZEN;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, Stz_Img);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal satellite zenith angle*/
      //test_file_name = "UT_satellite_zenith_angle_" + Result_FN;
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_satellite_zenith_angle_" + suffix_str;
      test_file.open(test_file_name.c_str());
      Stz_Img.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_PIXEL_SOLAR_ZEN;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, Soz_Img);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal solar zenith angle*/
      //test_file_name = "UT_solar_zenith_angle_" + Result_FN;
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_solar_zenith_angle_" + suffix_str;
      test_file.open(test_file_name.c_str());
      Soz_Img.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   /*Check the existence of land/sea mask in GSIP L1 product
    *if so, read land/sea mask in HDF format
    *otherwise, read alternative land/sea mask
    */

   /* read alternative land/sea mask*/
   if (param::instance.IS_GSIP_LAND_MASK == 0)
   {
      dim_name = param::instance.DIM_NAME_PIXEL_LAND_MASK;
      Get_Alternative_Landsea_mask();
   } else /* read GSIP land/sea mask*/
   {
      dim_name = param::instance.DIM_NAME_PIXEL_LAND_MASK;
      Get_GSIP_V3_Byte(GSIP_Data_V3_FN, dim_name, LandMask_Img);
   }
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal land/sea mask*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_land_Sea_mask";
      test_file.open(test_file_name.c_str());
      LandMask_Img.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   /*GSIP L2 product:
    cloud mask
    */
   dim_name = param::instance.DIM_NAME_PIXEL_CLOUD_MASK;
   Get_GSIP_V3_Byte(GSIP_Data_V3_FN, dim_name, Cloud_Mask_Img);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal cloud mask*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_cloud_mask";
      test_file.open(test_file_name.c_str());
      Cloud_Mask_Img.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }


   /*GSIP L3 product:
    snow fraction
    water vapor
    aerosol content
    */

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      IS_GSIP_SNOW_TPW = 1;
      /*Diagnostic code end*/
   }
   dim_name = param::instance.DIM_NAME_PIXEL_SNOW_FRAC;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, Snowmask_Img);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal snow mask*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_snow_mask";
      test_file.open(test_file_name.c_str());
      Snowmask_Img.MatrixWrite(test_file);
      SampleLine_L3.MatrixWrite(test_file);
      SamplePixel_L3.MatrixWrite(test_file);
      test_file.close();
      IS_GSIP_SNOW_TPW++;
      /*Diagnostic code end*/
   }

   dim_name = param::instance.DIM_NAME_PIXEL_TPW;
   Get_GSIP_V3_Float(GSIP_Data_V3_FN, dim_name, TPW_Img);
   param::instance << getTimeString() << dim_name << " load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal TPW*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_TPW";
      test_file.open(test_file_name.c_str());
      TPW_Img.MatrixWrite(test_file);
      SampleLine_L3.MatrixWrite(test_file);
      SamplePixel_L3.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   /*Read and preprocess MODIS emissivity DATA*/
   Get_Emissivity_RAW();
   param::instance << getTimeString() << "emissivity file load completion."
         << endl;

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*output internal emissivity*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_emissivity"
            + suffix_str;
      test_file.open(test_file_name.c_str());
      Em_TP2_Resampled.MatrixWrite(test_file);
      Em_TP4_Resampled.MatrixWrite(test_file);
      SampleLine_L3.MatrixWrite(test_file);  // yuling
      SamplePixel_L3.MatrixWrite(test_file);  // yuling
      test_file.close();
      /*Diagnostic code end*/
   }

   /*Read and process coefficient DATA*/
   GetCoefData();
   param::instance << getTimeString()
         << "regression coefficient file load completion." << endl;

   return true; /* success data input reached */
}

/*Read alternative land/sea mask*/
bool LST_Retrieval::Get_Alternative_Landsea_mask()
{
   /*0   Water (and Goode's interrupted space)
    * 1   Evergreen Needleleaf Forest
    * 2   Evergreen Broadleaf Foreset
    * 3   Deciduous Needleleaf Forest
    * 4   Deciduous Broadleaf Forest
    * 5   Mixed Forest
    * 6   Woodland
    * 7   Wooded Grassland
    * 8   Closed Shrubland
    * 9   Open Shrubland
    * 10   Grassland
    * 11   Cropland
    * 12   Bare Ground
    * 13   Urban and Built-up
    */

   /*read global land/sea mask*/
   CMatrix<char> globalMask(param::instance.ALTERNATIVE_LAND_MASK_HEIGHT,
         param::instance.ALTERNATIVE_LAND_MASK_WIDTH);
   ifstream file;
   string landseamask_fn = getEnvDefine("GLST_STATIC_DIR")
         + param::instance.ALTERNATIVE_LAND_MASK_NAME;
   file.open(landseamask_fn.c_str());
   if (file.fail())
   {
      param::instance << "land mask read failure!" << endl;
      exit(2); /* error return */
   }
   globalMask.MatrixsRead(file);
   file.close();

   double SampleLine = 0.0;
   double SamplePixel = 0.0;
   int SampleLine_i = 0;
   int SamplePixel_i = 0;
   char cur_landmask = 0;

   LandMask_Img.Resize(ImgHeight, ImgWidth);

   for (int i = 0; i < ImgHeight; i++)
   {
      for (int j = 0; j < ImgWidth; j++)
      {
         float cur_lon = lon_degree.GetElement(i, j);
         float cur_lat = lat_degree.GetElement(i, j);
         if ((cur_lon == -1) || (cur_lat == -1))
         {
            cur_landmask = -1;
         } else
         {
            SamplePixel = (180 + cur_lon)
                  / param::instance.ALTERNATIVE_LAND_MASK_RESOLUTION;
            SampleLine = (90 - cur_lat)
                  / param::instance.ALTERNATIVE_LAND_MASK_RESOLUTION;
            SamplePixel_i = (int) floor(SamplePixel);
            SampleLine_i = (int) floor(SampleLine);

            if ((SamplePixel < 0) || (SampleLine < 0) || (SamplePixel
                  > param::instance.ALTERNATIVE_LAND_MASK_WIDTH) || (SampleLine
                  > param::instance.ALTERNATIVE_LAND_MASK_HEIGHT))
            {
               cur_landmask = -1;
            } else
            {
               cur_landmask
                     = globalMask.GetElement(SampleLine_i, SamplePixel_i);
            }
         }
         LandMask_Img.SetElement(i, j, cur_landmask);
      }
   }

   if (param::instance.IsTest)
   {
      ofstream outputfile;
      string test_file_name = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_land_sea_mask.dat";
      outputfile.open(test_file_name.c_str());
      LandMask_Img.MatrixWrite(outputfile);
      outputfile.close();
   }

   return true; /* success alternative land/sea mask read reached */

}

/*Get the image information: image size and mode*/
/* Modified by Yuling Liu for GSIP V3 format update */
bool LST_Retrieval::GetImgInfor()
{
   /*Variable assignment for image height and image width*/
   ImgHeight = GOES_TP2_Aft_Cal.GetNumRows();
   ImgWidth = GOES_TP2_Aft_Cal.GetNumColumns();
   if ((ImgHeight == param::instance.DEFAULT_IMG_HEIGHT_CONUS) && (ImgWidth
         == param::instance.DEFAULT_IMG_WIDTH_CONUS))
   {
      GOES_MODE = 1; // CONUS
   } else if ((ImgHeight == param::instance.DEFAULT_IMG_HEIGHT_FD) && (ImgWidth
         == param::instance.DEFAULT_IMG_WIDTH_FD))
   {
      GOES_MODE = 0; // FD
   } else
   {
      param::instance << getTimeString() << "Incorrect GOES image size."
            << endl;
      exit(2); /* error return for incorrect GOES image size*/
   }

   /*variable assignment for image date and image time*/
   int temp_goes_id = 0;
   int status = sscanf(GSIP_Data_V3_FN_ONLY.c_str(),
         param::instance.GSIP_V3_PRODUCT_NAME_FORMAT.c_str(), &temp_goes_id,
         scanmode, &img_date, &img_time);

   if (status == 0)
   {
      param::instance << getTimeString() << "GSIP file name convert fail."
            << endl;
      exit(2); /* error return for incorrect GSIP file interpretation*/
   }
   if (temp_goes_id != GOES_ID)
   {
      param::instance << getTimeString()
            << "The satellite ID of GSIP products doesn't match the assigned GOES ID"
            << endl;
      exit(2); /* error return*/
   }

  
   return true; /* success image information interpretation reached */
}

// modified by Yuling on Aug. 29, 2014 for GSIP V3 change
bool LST_Retrieval::SetOutputFilename()
{
   string temp_result_name = "glstL3_g";
   char buf[255];
   sprintf(buf, "%d_", GOES_ID);
   temp_result_name += buf;
   temp_result_name += scanmode;

   sprintf(buf, "_%d_%04d_v2", img_date, img_time);
   temp_result_name += buf;
   Result_FN = getEnvDefine("GLST_OUTPUT_DIR") + temp_result_name;

   return true; /*success output file name set reached*/
}

/*Read in variable (float) from GSIP data
 IN: GSIP data file name;
 IN: Dim Name;
 IN: Array to read in the variable
 */
bool LST_Retrieval::Get_GSIP_L1_L2_Float(string file_name,
      string DimName,
      FMatrix& ScaledData)
{
   int32 sd_id = Error_HDF(SDstart(file_name.c_str(), DFACC_READ), file_name);
   int32 sds_index =
         Error_HDF(SDnametoindex(sd_id, DimName.c_str()), file_name);
   int32 sds_id = Error_HDF(SDselect(sd_id, sds_index), file_name);
   int32 edges[2] =
   { 0 , 0 };
   int32 start[2] =
   { 0 , 0 };
   char sds_name[255];
   memset(sds_name, 0, sizeof(char) * 255);
   int32 var_rank = 0;
   int32 sd_data_type = 0;
   int32 num_attri = 0;
   Error_HDF(
         SDgetinfo(sds_id, sds_name, &var_rank, edges, &sd_data_type,
               &num_attri), file_name);

   SMatrix data_orig;
   BMatrix data_orig2;
   if ((DimName == param::instance.DIM_NAME_PIXEL_SENSOR_ZEN) || (DimName
         == param::instance.DIM_NAME_PIXEL_SOLAR_ZEN))
   {
      data_orig2.Resize(edges[0], edges[1]);
      Error_HDF(
            SDreaddata(sds_id, start, NULL, edges, (VOIDP) data_orig2.GetData()),
            file_name);
   } else
   {
      data_orig.Resize(edges[0], edges[1]);
      Error_HDF(
            SDreaddata(sds_id, start, NULL, edges, (VOIDP) data_orig.GetData()),
            file_name);
   }
   ScaledData.Resize(edges[0], edges[1]);

   char if_scale = 0;
   float RANGE_MIN = 0.0;
   float RANGE_MAX = 0.0;
   int32 SCALED_MIN = 0;
   int32 SCALED_MAX = 0;
   int32 missingValue = 0;
   short tempvalue = 0;
   float scaledValue = 0.0;
   int32 attr_index = Error_HDF(
         SDfindattr(sds_id, param::instance.HDF_ATTRI_NAME_IS_SCALED.c_str()),
         file_name);
   Error_HDF(SDreadattr(sds_id, attr_index, &if_scale), file_name);
   if (if_scale != 0)
   {
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_RANGE_MIN.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &RANGE_MIN));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_RANGE_MAX.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &RANGE_MAX));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_SCALED_MIN.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &SCALED_MIN));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_SCALED_MAX.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &SCALED_MAX));
   }

   attr_index = Error_HDF(
         SDfindattr(sds_id,
               param::instance.HDF_ATTRI_NAME_SCALED_MISSING.c_str()),
         file_name);
   Error_HDF(SDreadattr(sds_id, attr_index, &missingValue), file_name);

   for (int32 i = 0; i < edges[0]; i++)
   {
      for (int32 j = 0; j < edges[1]; j++)
      {
         if ((DimName == param::instance.DIM_NAME_PIXEL_SENSOR_ZEN) || (DimName
               == param::instance.DIM_NAME_PIXEL_SOLAR_ZEN))
            tempvalue = data_orig2.GetElement(i, j);
         else
            tempvalue = data_orig.GetElement(i, j);
         if (if_scale)
         {
            if (tempvalue != missingValue)
            {
               scaledValue = (((float) tempvalue - SCALED_MIN) / (SCALED_MAX
                     - SCALED_MIN)) * (RANGE_MAX - RANGE_MIN) + RANGE_MIN;
            } else
            {
               scaledValue = param::instance.INVALID_VALUE_DOUBLE;
            }
         } else
         {
            if (tempvalue != missingValue)
            {
               scaledValue = tempvalue;
            } else
            {
               scaledValue = param::instance.INVALID_VALUE_DOUBLE;
            }
         }
         ScaledData.SetElement(i, j, scaledValue);
      }
   }

   /* Terminate access to the array. */
   Error_HDF(SDendaccess(sds_id), file_name);
   /* Terminate access to the SD interface and close the file. */
   Error_HDF(SDend(sd_id), file_name);

   return true; /* success GSIP product read reached */
}

int32 LST_Retrieval::Error_HDF(int32 stat)
{
   if (stat == FAIL)
   {
      param::instance << getTimeString() << " HDF file read fail." << endl;
      exit(2); /* error return for HDF file read failure*/
   }
   return stat; /* status return */
}

int32 LST_Retrieval::Error_HDF(int32 stat, std::string filename)
{
   if (stat == FAIL)
   {
      param::instance << getTimeString() << "  " << filename
            << " HDF file read fail." << endl;
      exit(2); /* error return for HDF file read failure*/
   }
   return stat; /* status return */
}

/*Read in variable (byte) from GSIP data
 IN: GSIP data file name;
 IN: Dim Name;
 IN: Array to read in the variable
 */
bool LST_Retrieval::Get_GSIP_L1_L2_Byte(string file_name,
      string DimName,
      BMatrix& ScaledData)
{
   int32 sd_id = Error_HDF(SDstart(file_name.c_str(), DFACC_READ), file_name);
   int32 sds_index =
         Error_HDF(SDnametoindex(sd_id, DimName.c_str()), file_name);
   int32 sds_id = Error_HDF(SDselect(sd_id, sds_index), file_name);
   int32 edges[2] =
   { 0 , 0 };
   int32 start[2] =
   { 0 , 0 };
   char sds_name[255];
   memset(sds_name, 0, sizeof(char) * 255);
   int32 var_rank = 0;
   int32 sd_data_type = 0;
   int32 num_attri = 0;

   Error_HDF(
         SDgetinfo(sds_id, sds_name, &var_rank, edges, &sd_data_type,
               &num_attri), file_name);

   BMatrix data_orig(edges[0], edges[1]);
   ScaledData.Resize(edges[0], edges[1]);
   Error_HDF(
         SDreaddata(sds_id, start, NULL, edges, (VOIDP) data_orig.GetData()),
         file_name);
   char if_scale;
   float RANGE_MIN = 0;
   float RANGE_MAX = 0;
   int32 SCALED_MIN = 0;
   int32 SCALED_MAX = 0;
   int32 missingValue = 0;
   char temp_value = 0;
   char ScaledValue = static_cast<char> (0.0);
   int32 attr_index = Error_HDF(
         SDfindattr(sds_id, param::instance.HDF_ATTRI_NAME_IS_SCALED.c_str()),
         file_name);
   Error_HDF(SDreadattr(sds_id, attr_index, &if_scale));
   if (if_scale != 0)
   {
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_RANGE_MIN.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &RANGE_MIN));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_RANGE_MAX.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &RANGE_MAX));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_SCALED_MIN.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &SCALED_MIN));
      attr_index
            = Error_HDF(
                  SDfindattr(sds_id,
                        param::instance.HDF_ATTRI_NAME_SCALED_MAX.c_str()),
                  file_name);
      Error_HDF(SDreadattr(sds_id, attr_index, &SCALED_MAX));
   }

   attr_index
         = Error_HDF(
               SDfindattr(sds_id,
                     param::instance.HDF_ATTRI_NAME_RANG_MISSING.c_str()),
               file_name);
   Error_HDF(SDreadattr(sds_id, attr_index, &missingValue), file_name);

   for (int32 i = 0; i < edges[0]; i++)
   {
      for (int32 j = 0; j < edges[1]; j++)
      {
         temp_value = data_orig.GetElement(i, j);
         if (if_scale)
         {
            if (temp_value != missingValue)
            {
               ScaledValue =static_cast<char> ( ((temp_value - SCALED_MIN) / (SCALED_MAX
                     - SCALED_MIN)) * (RANGE_MAX - RANGE_MIN) + RANGE_MIN );
            } else
            {
               ScaledValue = param::instance.INVALID_VALUE_INT;
            }
         } else
         {
            if (temp_value != missingValue)
            {
               ScaledValue = temp_value;
            } else
            {
               ScaledValue = param::instance.INVALID_VALUE_INT;
            }
         }
         ScaledData.SetElement(i, j, ScaledValue);
      }
   }

   /* Terminate access to the array. */
   Error_HDF(SDendaccess(sds_id), file_name);
   /* Terminate access to the SD interface and close the file. */
   Error_HDF(SDend(sd_id), file_name);
   return true; /* success GSIP product read reached */
}

template<class ct>
ct* Error_NetCDF(ct* stat)
{
   if (stat == NULL)
   {
      param::instance << getTimeString() << " NetCDF file read fail." << endl;
      exit(2); /* error return for NetCDF file read failure*/
   }
   return stat; /* status return */
}

bool LST_Retrieval::Error_NetCDF_b(bool stat)
{
   if (stat == FALSE)
   {
      param::instance << getTimeString() << " NETCDF file read fail." << endl;
      exit(2); /* error return for HDF file read failure*/
   }
   return stat; /* status return */
}

/*Read in GSIP L3 product
 IN: GSIP data file name;
 IN: Dim Name;
 IN: Array to read in the variable
 */
bool LST_Retrieval::Get_GSIP_L3(string file_name,
      string DimName,
      FMatrix& ScaledData)
{

   NcFile nc(file_name.c_str(), NcFile::ReadOnly);
   if (nc.is_valid() == FALSE)
   {
      param::instance << getTimeString() << " NetCDF file open fail." << endl;
      exit(2); /* error return */
   }

   NcDim* dim = Error_NetCDF(nc.get_dim("nlon_grid"));
   int32 array_lon = dim->size();
   dim = Error_NetCDF(nc.get_dim("nlat_grid"));
   int32 array_lat = dim->size();
   NcVar* ncvar = Error_NetCDF(nc.get_var(DimName.c_str()));

   NcAtt* att = Error_NetCDF(ncvar->get_att("RANGE_MIN"));
   NcValues* ncvalue = att->values();
   float RANGE_MIN = ncvalue->as_float(0);

   att = Error_NetCDF(ncvar->get_att("RANGE_MAX"));
   ncvalue = att->values();
   float RANGE_MAX = ncvalue->as_float(0);

   att = Error_NetCDF(ncvar->get_att("SCALED_MIN"));
   ncvalue = att->values();
   int32 SCALED_MIN = ncvalue->as_int(0);

   att = Error_NetCDF(ncvar->get_att("SCALED_MAX"));
   ncvalue = att->values();
   int32 SCALED_MAX = ncvalue->as_int(0);

   long array_size = array_lon * array_lat;
   ncbyte ch2_values[array_size];
   float lon_cel[array_size];
   memset(lon_cel, 0, sizeof(float) * array_size);
   float lat_cel[array_size];
   memset(lat_cel, 0, sizeof(float) * array_size);
   long array[] =
   { array_lat , array_lon };
   Error_NetCDF_b(ncvar->get(ch2_values, array));
   ncvar = Error_NetCDF(nc.get_var("lon_cell"));
   Error_NetCDF_b(ncvar->get(lon_cel, array));
   ncvar = Error_NetCDF(nc.get_var("lat_cell"));
   Error_NetCDF_b(ncvar->get(lat_cel, array));

   float delta_lat = 0.0;
   float delta_lon = 0.0;
   float north_lat = 0.0;
   float west_lon = 0.0;
   float south_lat = 0.0;
   float east_lon = 0.0;
   ncvalue = Error_NetCDF(nc.get_var("delta_lat"))->values();
   delta_lat = ncvalue->as_float(0);
   ncvalue = Error_NetCDF(nc.get_var("delta_lon"))->values();
   delta_lon = ncvalue->as_float(0);
   ncvalue = Error_NetCDF(nc.get_var("north_lat"))->values();
   north_lat = ncvalue->as_float(0);
   ncvalue = Error_NetCDF(nc.get_var("south_lat"))->values();
   south_lat = ncvalue->as_float(0);
   ncvalue = Error_NetCDF(nc.get_var("west_lon"))->values();
   west_lon = ncvalue->as_float(0);
   ncvalue = Error_NetCDF(nc.get_var("east_lon"))->values();
   east_lon = ncvalue->as_float(0);

   long lon_boundary =long(((east_lon - west_lon) / delta_lon)) + 1;
   long lat_boundary =long(((north_lat - south_lat) / delta_lat)) + 1;
   FMatrix Scaled_Data_Origin(lat_boundary, lon_boundary);
   for (int32 i = 0; i < lat_boundary; i++)
   {
      for (int32 j = 0; j < lon_boundary; j++)
      {
         Scaled_Data_Origin.SetElement(i, j,
               param::instance.INVALID_VALUE_DOUBLE);
      }
   }
   ScaledData.Resize(ImgHeight, ImgWidth);
   float Unscaled_Value = 0.0;
   int32 lon_samp = 0;
   int32 lat_samp = 0;

   /*Diagnostic code begin*/
   FILE* test_insert_loc_file;
   /*Diagnostic code end*/

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      test_insert_loc_file = fopen(
            "getEnvDefine(\"GLST_OUTPUT_DIR\")+UT_insert_loc.txt", "wt");
      fprintf(test_insert_loc_file, "boundary:\rwest_lon:%f\t\tnorth_lat:%f\r",
            west_lon, north_lat);
      fprintf(test_insert_loc_file,
            "No\tlon\tlat\tvalue_s\tvalue_u\tlon_samp\tlat_samp\r");
      /*Diagnostic code end*/
   }

   for (int32 i = 0; i < array_lon; i++)
   {
      char scaled_value = ch2_values[i];
      Unscaled_Value = (((float) scaled_value - SCALED_MIN) / (SCALED_MAX
            - SCALED_MIN)) * (RANGE_MAX - RANGE_MIN) + RANGE_MIN;
      lon_samp = (int32) floor((lon_cel[i] - west_lon) / delta_lon);
      lat_samp = (int32) floor((north_lat - lat_cel[i]) / delta_lat);
      Scaled_Data_Origin.SetElement((int32) lat_samp, (int32) lon_samp,
            Unscaled_Value);
      if (param::instance.IsTest)
      {
         /*Diagnostic code begin*/
         fprintf(test_insert_loc_file, "%d\t%f\t%f\t%d\t%d\r", i, lon_cel[i],
               lat_cel[i], lon_samp, lat_samp);
         /*Diagnostic code end*/
      }
   }

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      fclose(test_insert_loc_file);
      SampleLine_L3.Resize(ImgHeight, ImgWidth);
      SamplePixel_L3.Resize(ImgHeight, ImgWidth);
      string outputfilename = "";
      if (IS_GSIP_SNOW_TPW == 1)
      {
         outputfilename = getEnvDefine("GLST_OUTPUT_DIR")
               + "UT_snow_original_L3.dat";
      } else if (IS_GSIP_SNOW_TPW == 2)
      {
         outputfilename = getEnvDefine("GLST_OUTPUT_DIR")
               + "UT_tpw_original_L3.dat";
      }
      ofstream test_file;
      test_file.open(outputfilename.c_str());
      Scaled_Data_Origin.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   ImageSampler32f DataSampler(lon_boundary, lat_boundary,
         Scaled_Data_Origin.GetData(), false);
   float lon_cur = 0.0;
   float lat_cur = 0.0;
   for (int32 i = 0; i < ImgHeight; i++)
   {
      for (int32 j = 0; j < ImgWidth; j++)
      {
         lon_cur = lon_degree.GetElement(i, j);
         lat_cur = lat_degree.GetElement(i, j);
         int32 SamplePixel = (int32) floor((lon_cur - west_lon) / delta_lon);
         int32 SampleLine = (int32) floor((north_lat - lat_cur) / delta_lat);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            SampleLine_L3.SetElement(i, j, (float) SampleLine);
            SamplePixel_L3.SetElement(i, j, (float) SamplePixel);
            /*Diagnostic code end*/
         }

         if ((SamplePixel < 0) || (SampleLine < 0) || (SamplePixel
               > lon_boundary) || (SampleLine > lat_boundary))
         {
            Unscaled_Value = param::instance.INVALID_VALUE_DOUBLE;
         } else
         {
            Unscaled_Value = DataSampler.ResampleNN(SamplePixel, SampleLine);
         }
         ScaledData.SetElement(i, j, Unscaled_Value);
      }
   }

   return true; /* success GSIP product read reached */
}

/*Read in global MODIS emissivity data in RAW format, and pre-process the emissivity data*/
bool LST_Retrieval::Get_Emissivity_RAW()
{
   DMatrix Em_TP2_Img;
   DMatrix Em_TP4_Img;

   Get_Emissivity_HDF(param::instance.DIM_NAME_EMIS_22, Em_TP2_Img);
   Get_Emissivity_HDF(param::instance.DIM_NAME_EMIS_31, Em_TP4_Img);

   Em_TP2_Resampled.Resize(ImgHeight, ImgWidth);
   Em_TP4_Resampled.Resize(ImgHeight, ImgWidth);

   ImageSampler64f emSampler_tp2(param::instance.MODIS_EMI_WIDTH,
         param::instance.MODIS_EMI_HEIGHT, Em_TP2_Img.GetData(), false);
   ImageSampler64f emSampler_tp4(param::instance.MODIS_EMI_WIDTH,
         param::instance.MODIS_EMI_HEIGHT, Em_TP4_Img.GetData(), false);

/*  yuling

   ImageSampler64f emSampler_tp2(param::instance.MODIS_EMI_HEIGHT,
         param::instance.MODIS_EMI_WIDTH, Em_TP2_Img.GetData(), false);
   ImageSampler64f emSampler_tp4(param::instance.MODIS_EMI_HEIGHT,
         param::instance.MODIS_EMI_WIDTH, Em_TP4_Img.GetData(), false);

*/
   double SamplePixel = 0.0;
   double SampleLine = 0.0;
   float Em_Value2 = 0.0;
   float Em_Value4 = 0.0;
   /* added by Yuling Liu for GSIP V3 update */
   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      SampleLine_L3.Resize(ImgHeight, ImgWidth);
      SamplePixel_L3.Resize(ImgHeight, ImgWidth);
   }

   for (int32 i = 0; i < ImgHeight; i++)
   {
      for (int32 j = 0; j < ImgWidth; j++)
      {
         float lon_temp = lon_degree.GetElement(i, j);
         float lat_temp = lat_degree.GetElement(i, j);
         SamplePixel = (lon_temp - param::instance.MODIS_LON_LEFT_UPPER_CORNER)
               / param::instance.MODIS_EMI_RESOLUTION;
         SampleLine = (param::instance.MODIS_LAT_LEFT_UPPER_CORNER - lat_temp)
               / param::instance.MODIS_EMI_RESOLUTION;


         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
	     /* modified by Yuling Liu for GSIP V3 update */
            SampleLine_L3.SetElement(i, j, (int) floor(SampleLine));
            SamplePixel_L3.SetElement(i, j, (int) floor(SamplePixel));
            /*Diagnostic code end*/
         }

         if ((SamplePixel < 0) || (SampleLine < 0) || (SamplePixel
               > param::instance.MODIS_EMI_WIDTH) || (SampleLine
               > param::instance.MODIS_EMI_HEIGHT))
         {
            Em_Value2 = param::instance.INVALID_VALUE_DOUBLE;
            Em_Value4 = param::instance.INVALID_VALUE_DOUBLE;
         } else
         {
            Em_Value2 = emSampler_tp2.ResampleNN(SamplePixel, SampleLine);
            Em_Value4 = emSampler_tp4.ResampleNN(SamplePixel, SampleLine);
         }

         Em_TP2_Resampled.SetElement(i, j, Em_Value2);
         Em_TP4_Resampled.SetElement(i, j, Em_Value4);
      }
   }

   return true; /* success emissivity data read reached */
}

/*Read in MODIS emissivity data in HDF format*/
bool LST_Retrieval::Get_Emissivity_HDF(string DimName, DMatrix& ScaledData)
{
   int32 sd_id = Error_HDF(SDstart(Em_Img_FN.c_str(), DFACC_READ), Em_Img_FN);
   int32 sds_index =
         Error_HDF(SDnametoindex(sd_id, DimName.c_str()), Em_Img_FN);
   int32 sds_id = Error_HDF(SDselect(sd_id, sds_index), Em_Img_FN);
   int32 edges[2] =
   { 0 , 0 };
   int32 start[2] =
   { 0 , 0 };
   char sds_name[255];
   memset(sds_name, 0, sizeof(char) * 255);
   int32 var_rank = 0;
   int32 sd_data_type = 0;
   int32 num_attri = 0;

   Error_HDF(
         SDgetinfo(sds_id, sds_name, &var_rank, edges, &sd_data_type,
               &num_attri), Em_Img_FN);

   ScaledData.Resize(edges[0], edges[1]);
   Error_HDF(
         SDreaddata(sds_id, start, NULL, edges, (VOIDP) ScaledData.GetData()),
         Em_Img_FN);

   /* Terminate access to the array. */
   Error_HDF(SDendaccess(sds_id), Em_Img_FN);
   /* Terminate access to the SD interface and close the file. */
   Error_HDF(SDend(sd_id), Em_Img_FN);

   return true; /* success emissivity data read reached */
}

/*get file name for algorithm coefficients*/
bool LST_Retrieval::GetCoefFilename(string& Coeff_Day_FN,
      string& Coeff_Night_FN)
{
   int nYear = 0;
   int DayOfYear = 0;
   int nMonth = 0;
   //int nDay = 0;
   DayOfYear = img_date % 1000;
   nYear = (int) img_date / 1000;

   DayOfYear2MonthDay(nYear, DayOfYear, nMonth, nDay);
   if ((nMonth >= 3) && (nMonth <= 5))
   {
      Coeff_Day_FN = param::instance.COEF_FILENAME_GOES13_SPRING;
      Coeff_Night_FN = param::instance.COEF_FILENAME_GOES13_SPRING_CORRECTION;
   } else if ((nMonth >= 6) && (nMonth <= 8))
   {
      Coeff_Day_FN = param::instance.COEF_FILENAME_GOES13_SUMMER;
      Coeff_Night_FN = param::instance.COEF_FILENAME_GOES13_SUMMER_CORRECTION;
   } else if ((nMonth >= 9) && (nMonth <= 11))
   {
      Coeff_Day_FN = param::instance.COEF_FILENAME_GOES13_FALL;
      Coeff_Night_FN = param::instance.COEF_FILENAME_GOES13_FALL_CORRECTION;
   } else
   {
      Coeff_Day_FN = param::instance.COEF_FILENAME_GOES13_WINTER;
      Coeff_Night_FN = param::instance.COEF_FILENAME_GOES13_WINTER_CORRECTION;
   }

   return true; /*success automatically match coefficients file names*/
}

/*Read in coefficient of regression model from *.txt file*/
bool LST_Retrieval::GetCoefData()
{
   vector<string> columns;
   columns.push_back("tp4");
   columns.push_back("bt_dif");
   columns.push_back("bt_dif2");
   columns.push_back("stz");
   columns.push_back("soz");
   columns.push_back("emi");
   columns.push_back("emi_dif");

   vector<string> columns_correction;
   columns_correction.push_back("goeslst");
   columns_correction.push_back("bt_3.9");
   columns_correction.push_back("bt_11");
   columns_correction.push_back("emi_3.9");
   columns_correction.push_back("emi_11");
   columns_correction.push_back("stz");
   columns_correction.push_back("soz");

   string Coeff_AllDay_FN = "";
   string Coeff_Correction_FN = "";
   GetCoefFilename(Coeff_AllDay_FN, Coeff_Correction_FN);

   root_allday = parser(getEnvDefine("GLST_STATIC_DIR") + Coeff_AllDay_FN,
         columns);
   if (root_allday == NULL)
   {
      param::instance << getTimeString()
            << "daytime training coefficients file read failure." << endl;
      exit(2); /* error return for file open failure*/
   }
   root_correction = parser(
         getEnvDefine("GLST_STATIC_DIR") + Coeff_Correction_FN,
         columns_correction);
   if (root_correction == NULL)
   {
      param::instance << getTimeString()
            << "nighttime training coefficients file read failure." << endl;
      exit(2); /* error return for file open failure*/
   }
   param::instance << getTimeString()
         << "regression coefficients file load completion." << endl;
   return true; /* success  coefficients data read reached */
}

/*Retrieve LST from GSIP products*/
void LST_Retrieval::Get_LST_Dual_Window()
{
   /*get data from GSIP product in HDF format*/
   if (GetData() == FALSE)
   {
      return; /* error return for data read failure*/
   }

   float tp2_cpixel = 0.0;
   float tp4_cpixel = 0.0;
   float stz_cpixel = 0.0;
   float soz_cpixel = 0.0;
   float em_tp2_cpixel = 0.0;
   float em_tp4_cpixel = 0.0;
   char cloud_mask_cpixel = 0;
   int32 IsLand = 0;
   float coef[7];
   memset(coef, 0, sizeof(float) * 7);
   float coef_correction[7];
   memset(coef_correction, 0, sizeof(float) * 7);
   Condition* root_allday_temp = NULL;
   Condition* root_correction_temp = NULL;
   float RetrivalTemp = 0.0;
   short ScaledValue = 0;
   LST_Retrievals.Resize(ImgHeight, ImgWidth);
   LST_before_scale.Resize(ImgHeight, ImgWidth);

   /*Diagnostic code begin*/
   FILE * test_coef_night;
   FILE * test_coef_day;
 
   /*Diagnostic code end*/
   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      string temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_test_coef_"
            + suffix_str + "_nighttime.csv";
      test_coef_night = fopen(temp_fl.c_str(), "wt");
      fprintf(
            test_coef_night,
            "soz,tp4,tp2,stz,em_tp2,em_tp4,a0,a1,a2,a3,a4,a5,a6,coef0,coef1,coef2,coef3,coef4,coef5,coef6,coef7,retrieval,a4_coef4\r");

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_test_coef_" + suffix_str
            + "_daytime.csv";
      test_coef_day = fopen(temp_fl.c_str(), "wt");
      fprintf(
            test_coef_day,
            "soz,tp4,tp2,stz,em_tp2,em_tp4,a0,a1,a2,a3,a4,a5,a6,coef0,coef1,coef2,coef3,coef4,coef5,coef6,coef7,retrieval,a4_coef4\r");
      /*Diagnostic code end*/
   }



   for (int32 iHeight = 0; iHeight < ImgHeight; iHeight++)
   {
      for (int32 iWidth = 0; iWidth < ImgWidth; iWidth++)
      {
         RetrivalTemp = param::instance.INVALID_VALUE_DOUBLE;
         ScaledValue = 0;
         IsLand = LandMask_Img.GetElement(iHeight, iWidth);
         if (IsLand > 0)
         {
            cloud_mask_cpixel = Cloud_Mask_Img.GetElement(iHeight, iWidth);
            if ((cloud_mask_cpixel == param::instance.THRESHOLD_CLOUD_CLEAR)
                  || (cloud_mask_cpixel
                        == param::instance.THRESHOLD_CLOUD_PRO_CLEAR))
            {
               tp2_cpixel = GOES_TP2_Aft_Cal.GetElement(iHeight, iWidth);
               tp4_cpixel = GOES_TP4_Aft_Cal.GetElement(iHeight, iWidth);
               stz_cpixel = Stz_Img.GetElement(iHeight, iWidth);
               soz_cpixel = Soz_Img.GetElement(iHeight, iWidth);

               em_tp2_cpixel = Em_TP2_Resampled.GetElement(iHeight, iWidth);
               em_tp4_cpixel = Em_TP4_Resampled.GetElement(iHeight, iWidth);
               if ((tp2_cpixel != param::instance.INVALID_VALUE_DOUBLE)
                     && (tp4_cpixel != param::instance.INVALID_VALUE_DOUBLE)
                     && (stz_cpixel != param::instance.INVALID_VALUE_DOUBLE)
                     && (stz_cpixel <= 89) && (soz_cpixel
                     != param::instance.INVALID_VALUE_DOUBLE) && (em_tp2_cpixel
                     != param::instance.INVALID_VALUE_DOUBLE) && (em_tp4_cpixel
                     != param::instance.INVALID_VALUE_DOUBLE))
               {
                  if ((tp2_cpixel
                        >= param::instance.THRESHOLD_GSIP_BT_NORMAL_LOW)
                        && (tp2_cpixel
                              <= param::instance.THRESHOLD_GSIP_BT_NORMAL_UPPER)
                        && (tp4_cpixel
                              >= param::instance.THRESHOLD_GSIP_BT_NORMAL_LOW)
                        && (tp4_cpixel
                              <= param::instance.THRESHOLD_GSIP_BT_NORMAL_UPPER))
                  {
                     float bt_diff = tp4_cpixel - tp2_cpixel;

                     coef[0] = tp4_cpixel;
                     coef[1] = bt_diff;
                     coef[2] = 0;
                     //coef[3] = 1 / cos(stz_cpixel * M_PI / 180) - 1;
                     coef[3] = cos(stz_cpixel * M_PI / 180);
                     //coef[3] = 1 / (1 - cos(stz_cpixel * M_PI / 180));
                     coef[4] = tp2_cpixel * cos(soz_cpixel * M_PI / 180);
                     //coef[4] = cos(soz_cpixel * M_PI / 180);
                     //coef[4] = bt_diff*cos(soz_cpixel * M_PI / 180);
                     coef[5] = 1 - 0.5 * (em_tp2_cpixel + em_tp4_cpixel);
                     //coef[5] = 1 - em_tp4_cpixel;
                     coef[6] = 0;
                     root_allday_temp = root_allday;
                     root_correction_temp = root_correction;

                     while (true)
                     {
                        if (root_allday_temp->model_coef == NULL)
                        {
                           if (coef[root_allday_temp->name_index]
                                 <= root_allday_temp->bound_val)
                           {
                              root_allday_temp = root_allday_temp->top;
                           } else
                              root_allday_temp = root_allday_temp->bottom;
                        } else
                        {
                           RetrivalTemp = root_allday_temp->model_coef[0]
                                 * coef[0] + root_allday_temp->model_coef[1]
                                 * coef[1] + root_allday_temp->model_coef[2]
                                 * coef[2] + root_allday_temp->model_coef[3]
                                 * coef[3] + root_allday_temp->model_coef[4]
                                 * coef[4] + root_allday_temp->model_coef[5]
                                 * coef[5] + root_allday_temp->model_coef[6]
                                 * coef[6] + root_allday_temp->model_coef[7];

                           break; /*success coefficients match reached*/
                        }
                     }

                     coef_correction[0] = RetrivalTemp;
                     coef_correction[1] = tp2_cpixel;
                     coef_correction[2] = tp4_cpixel;
                     coef_correction[3] = em_tp2_cpixel;
                     coef_correction[4] = em_tp4_cpixel;
                     coef_correction[5] = stz_cpixel;
                     coef_correction[6] = soz_cpixel;
                     float DT = 0.0;
                     while (true)
                     {
                        if (root_correction_temp->model_coef == NULL)
                        {
                           if (coef[root_correction_temp->name_index]
                                 <= root_correction_temp->bound_val)
                           {
                              root_correction_temp = root_correction_temp->top;
                           } else
                              root_correction_temp
                                    = root_correction_temp->bottom;
                        } else
                        {
                           DT = root_correction_temp->model_coef[0]
                                 * coef_correction[0]
                                 + root_correction_temp->model_coef[1]
                                       * coef_correction[1]
                                 + root_correction_temp->model_coef[2]
                                       * coef_correction[2]
                                 + root_correction_temp->model_coef[3]
                                       * coef_correction[3]
                                 + root_correction_temp->model_coef[4]
                                       * coef_correction[4]
                                 + root_correction_temp->model_coef[5]
                                       * coef_correction[5]
                                 + root_correction_temp->model_coef[6]
                                       * coef_correction[6]
                                 + root_correction_temp->model_coef[7];
                           break; /*success coefficients match reached*/
                        }
                     }
                     /*SYSTEMETIC CORRECTION*/
                     RetrivalTemp = RetrivalTemp - DT;
                     if (param::instance.IsTest)
                     {
                        /*Diagnostic code begin*/
                        if ((soz_cpixel >= 100) && (soz_cpixel < 170))
                        {
                           fprintf(
                                 test_coef_night,
                                 "%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\r",
                                 soz_cpixel, tp4_cpixel, tp2_cpixel,
                                 stz_cpixel, em_tp2_cpixel, em_tp4_cpixel,
                                 coef[0], coef[1], coef[2], coef[3], coef[4],
                                 coef[5], coef[6],
                                 root_allday_temp->model_coef[0],
                                 root_allday_temp->model_coef[1],
                                 root_allday_temp->model_coef[2],
                                 root_allday_temp->model_coef[3],
                                 root_allday_temp->model_coef[4],
                                 root_allday_temp->model_coef[5],
                                 root_allday_temp->model_coef[6],
                                 root_allday_temp->model_coef[7], RetrivalTemp,
                                 coef[4] * root_allday_temp->model_coef[4]);
                        }
                        if ((soz_cpixel >= 10) && (soz_cpixel < 70))
                        {
                           fprintf(
                                 test_coef_day,
                                 "%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\r",
                                 soz_cpixel, tp4_cpixel, tp2_cpixel,
                                 stz_cpixel, em_tp2_cpixel, em_tp4_cpixel,
                                 coef[0], coef[1], coef[2], coef[3], coef[4],
                                 coef[5], coef[6],
                                 root_allday_temp->model_coef[0],
                                 root_allday_temp->model_coef[1],
                                 root_allday_temp->model_coef[2],
                                 root_allday_temp->model_coef[3],
                                 root_allday_temp->model_coef[4],
                                 root_allday_temp->model_coef[5],
                                 root_allday_temp->model_coef[6],
                                 root_allday_temp->model_coef[7], RetrivalTemp,
                                 coef[4] * root_allday_temp->model_coef[4]);
                        }
                        /*Diagnostic code end*/
                     }

                  }/*END IF Normal GSIP BT check*/
               } /*ENDIF valid pixel checking*/
               /*large viewing angle filter*/
               /*if(stz_cpixel > 77)
                {
                RetrivalTemp = param::instance.INVALID_VALUE_DOUBLE;
                }*/
            } /*ENDIF cloud free pixel*/
         }/*ENDIF land pixel checking*/

         if (RetrivalTemp != param::instance.INVALID_VALUE_DOUBLE)
         {
            if (RetrivalTemp < param::instance.THRESHOLD_LST_RANGE_LOW)
            {
               RetrivalTemp = param::instance.THRESHOLD_LST_RANGE_LOW;
            }
            if (RetrivalTemp > param::instance.THRESHOLD_LST_RANGE_UPPER)
            {
               RetrivalTemp = param::instance.THRESHOLD_LST_RANGE_UPPER;
            }
            ScaledValue = static_cast<short> ( (RetrivalTemp
                  - param::instance.METADATA_LST_ATTRI_OFFSET)
                  / param::instance.METADATA_LST_ATTRI_FACTOR );
         }
         LST_Retrievals.SetElement(iHeight, iWidth, ScaledValue);
         LST_before_scale.SetElement(iHeight, iWidth, RetrivalTemp);

      } /*ENDFOR ImgWidth*/
   }/*ENDFOR ImgHeight*/

   param::instance << getTimeString() << " LST retrieval completion." << endl;

   /*Generate QC flag and metadata information*/

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      string test_file_name = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_LST_unscale_" + suffix_str;
      ofstream test_file;
      test_file.open(test_file_name.c_str());
      LST_before_scale.MatrixWrite(test_file);
      test_file.close();
      string temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_LST_scale.dat";
      test_file.open(temp_fl.c_str());
      LST_Retrievals.MatrixWrite(test_file);
      test_file.close();
      fclose(test_coef_night);
      fclose(test_coef_day);
      /*Diagnostic code end*/
   }

   /*Generate QC flag and metadata information*/
   Generate_QC_Metadata();
   param::instance << getTimeString()
         << " QC and Metadata generation completion." << endl;
   /*scale Lon/Lat for output*/
   ScaleLonLat();
   /*output LST in NetCDF format*/
   OutputNetCDF();
   param::instance << getTimeString() << " NetCDF output completion." << endl;
   /*output LST in GRIB2 format*/
   OutputGRIB2();   
   param::instance << getTimeString() << " GRIB2 output completion." << endl;
   param::instance << getTimeString() << "Successful completion" << endl;

   return; /* success LST retrieval reached */
}

bool LST_Retrieval::OutputGRIB2()
{

   GetGridLST_QC();
   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*test_file to output internal binary arrays*/
      ofstream test_file;
      string test_file_name;
      /*output internal grid LST*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_grid_LST";
      test_file.open(test_file_name.c_str());
      LST_grid.MatrixWrite(test_file);
      test_file.close();
      /*output internal grid qc*/
      test_file_name = getEnvDefine("GLST_OUTPUT_DIR") + "UT_grid_qc";
      test_file.open(test_file_name.c_str());
      QC_Matrix_grid.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   /*local parameters*/
   long ngrdpts = GridHeight * GridWidth; /*Number of grid points in the defined grid*/
   long lcgrib = ngrdpts * 4;
   unsigned char* cgrib = new unsigned char[lcgrib];

   /*set up indicator section*/
   /*
    * Discipline-GRIB Master Table Number (see Code Table 0.0)
    * listsec0[0] = 0  // 0 - meteorological products: temperature
    * listsec0[1] = 2  // GRIB Edition Number (currently 2)
    */
   g2int ListSec0[2] =
   { 0 , 2 };
   g2int ListSec1[13];
   memset(ListSec1, 0, sizeof(long) * 13);

   ListSec1[0] = 7; /* Id of originating centre [Common Code Table C-1]*/
   /* 7 - US NWS, NCEP*/
   ListSec1[1] = 2; /* Id of originating sub-centre [Sub-Centres table]*/
   /* 2 NCEPCentral Operations*/
   ListSec1[2] = 2; /* GRIB Master Tables Version Number [Code Table 1.0]*/
/*   ListSec1[3] = 0;  GRIB Local Tables Version Number [Code Table 1.1]*/
   ListSec1[3] = 1; /* Modified by yuling, local table should be set to 1 to be consistent with wgrib standard*/
   ListSec1[4] = 1; /* Significance of Reference Time [Code Table 1.2]*/

   //get year month day hour from imgDate and imgTime
   int DayOfYear = img_date % 1000;
   int nYear = (int) img_date / 1000;
   int nMonth = 0;
   int nDay = 0;
   DayOfYear2MonthDay(nYear, DayOfYear, nMonth, nDay);
   int nMin = img_time % 100;
   int nHour = (int) img_time / 100;
   ListSec1[5] = nYear; /* Reference Time - Year [4 digits]*/
   ListSec1[6] = nMonth; /* Reference Time - Month*/
   ListSec1[7] = nDay; /* Reference Time - Day*/
   ListSec1[8] = nHour; /* Reference Time - Hour*/
   ListSec1[9] = nMin; /* Reference Time - Minute*/
   ListSec1[10] = 0; /* Reference Time - Second*/
   ListSec1[11] = 0; /* Production status of data [Code Table 1.3]*/
   /* 0 - Operational products*/
   ListSec1[12] = 6; /* Type of processed data [Code Table 1.4]*/
   /* 6 - Processed satellite observations */

   /*set up grid definition section*/
   long igds[5];
   memset(igds, 0, sizeof(long) * 5);
   igds[0] = 0; /* Source of grid definition [see Code Table 3.0]*/
   igds[1] = ngrdpts; /* Number of grid points in the defined grid.*/
   igds[2] = 0; /* Number of octets needed for each additional grid points definition*/
   igds[3] = 0; /* Interpretation of list for optional points definition. */
   /* [Code Table 3.11]*/
   /* 0 - There is no appended list*/
   igds[4] = 0; /* Grid Definition Template Number [Code Table 3.1]*/
   /* 0 - Latitude/longitude*/

   long igdstmpl[19];
   memset(igdstmpl, 0, sizeof(long) * 19);
/* 
   long LatFirst = static_cast<long> (lat_max * 1000000);
   long LatEnd = long ((lat_max - param::instance.GRID_RESOLUTION * GridHeight)
         * 1000000 ); 
   long LonFirst = (lon_min + 360) * 1000000;
   long LonEnd = (lon_min - param::instance.GRID_RESOLUTION * GridWidth)
         * 1000000; 
*/
   /* modified by yuling, to make lon in range and consistent with wgrib etc */
   long LatFirst = long(lat_min * 1000000) ; 
   long LatEnd   = long(lat_max * 1000000) ; 
   long LonFirst = long(lon_min + 360) * 1000000 ;
   long LonEnd   = long(lon_max * 1000000) ;
 
   igdstmpl[0] = 6;
   igdstmpl[7] = GridWidth;
   igdstmpl[8] = GridHeight;
   igdstmpl[11] = LatFirst;
   igdstmpl[12] = LonFirst;
   igdstmpl[13] = 48;
   igdstmpl[14] = LatEnd;
   igdstmpl[15] = LonEnd;
   igdstmpl[16] = long (param::instance.GRID_RESOLUTION * 1000000 );
   igdstmpl[17] = long (param::instance.GRID_RESOLUTION * 1000000 );
   igdstmpl[18] = 64;

   /*set up product definition/representation/bit-map/data sections for temperature*/
   int ipdsnum = 0; /*product definition template number*/
   long ipdstmpl[15];
   memset(ipdstmpl, 0, sizeof(long) * 15);
   ipdstmpl[0] = 0; /*Parameter category,see Code table 4.1. 0 - temperature*/
   ipdstmpl[1] = 17; /*Parameter number,see Code table 4.2. 17-skin temperature K*/
   ipdstmpl[2] = 0; /*Type of generating process (see Code table 4.3).0 - Analysis*/
   ipdstmpl[3] = 255; /*Background generating process identifier*/
   ipdstmpl[4] = 255; /*Analysis or forecast generating process identified*/

   ipdstmpl[5] = 0; /*Hours of observational data cutoff after reference time*/
   ipdstmpl[6] = 0; /*Minutes of observational data cutoff after reference time*/
   ipdstmpl[7] = 1; /*Indicator of unit of time range (see Code table 4.4) */
   ipdstmpl[8] = 0; /*Forecast time in units defined by octet 18*/
   ipdstmpl[9] = 255; /*Type of first fixed surface (see Code table 4.5)*/
   ipdstmpl[10] = 0; /*Scale factor of first fixed surface*/
   ipdstmpl[11] = 0; /*Scaled value of first fixed surface*/
   ipdstmpl[12] = 255;/*Type of second fixed surfaced (see Code table 4.5)*/
   ipdstmpl[13] = 0; /*Scale factor of second fixed surface*/
   ipdstmpl[14] = 0; /*Scaled value of second fixed surfaces*/

   /*set up product definition/representation/bit-map/data sections for QC*/
   long ipdstmpl_qc[15];
   memset(ipdstmpl_qc, 0, sizeof(long) * 15);
   ipdstmpl_qc[0] = 255; /*Parameter category,see Code table 4.1.*/
   ipdstmpl_qc[1] = 255; /*Parameter number,see Code table 4.2.*/
   ipdstmpl_qc[2] = 0; /*Type of generating process (see Code table 4.3).0 - Analysis*/
   ipdstmpl_qc[3] = 255; /*Background generating process identifier*/
   ipdstmpl_qc[4] = 255; /*Analysis or forecast generating process identified*/
   ipdstmpl_qc[5] = 0; /*Hours of observational data cutoff after reference time*/
   ipdstmpl_qc[6] = 0; /*Minutes of observational data cutoff after reference time*/
   ipdstmpl_qc[7] = 1; /*Indicator of unit of time range (see Code table 4.4) */
   ipdstmpl_qc[8] = 0; /*Forecast time in units defined by octet 18*/
   ipdstmpl_qc[9] = 255; /*Type of first fixed surface (see Code table 4.5)*/
   ipdstmpl_qc[10] = 0; /*Scale factor of first fixed surface*/
   ipdstmpl_qc[11] = 0; /*Scaled value of first fixed surface*/
   ipdstmpl_qc[12] = 255;/*Type of second fixed surfaced (see Code table 4.5)*/
   ipdstmpl_qc[13] = 0; /*Scale factor of second fixed surface*/
   ipdstmpl_qc[14] = 0; /*Scaled value of second fixed surfaces*/

   float *coordlist = NULL;
   long numcoord = 0;
   /*Data Representation Template Number (see Code Table 5.0)
    *0 - Grid Point Data - Simple Packing */
   int idrsnum = 0; /*Data Representation Template Number*/

   /*Temperature ranges from 210.0K to 350.0K.
    * temperature filled value is -999.0.
    * D value is chosen 1.
    * All residuals are pack into "word" in 14bits long*/
   long idrstmpl[5];
   memset(idrstmpl, 0, sizeof(long) * 5);
   idrstmpl[0] = -999; /*Reference value (R) (IEEE 32-bit floating-point value)*/
   idrstmpl[1] = 0; /*Binary scale factor (E), not used*/
   idrstmpl[2] = 1; /*Decimal scale factor (D)*/
   idrstmpl[3] = 14; /*Number of bits used for each packed value for simple packing*/
   idrstmpl[4] = 0; /*Type of original field values. 0 - float*/

   /*QC ranges from 0 to 65536.
    * D value is chosen 1.
    * All residuals are pack into "word" in 16bits long*/
   long idrstmpl_qc[5];
   memset(idrstmpl_qc, 0, sizeof(long) * 5);
   idrstmpl_qc[0] = 0; /*Reference value (R) (IEEE 32-bit floating-point value)*/
   idrstmpl_qc[1] = 0; /*Binary scale factor (E), not used*/
   idrstmpl_qc[2] = 0; /*Decimal scale factor (D)*/
   idrstmpl_qc[3] = 16; /*Number of bits used for each packed value for simple packing*/
   idrstmpl_qc[4] = 1; /*Type of original field values. 1 - int*/

   int ibmap = 255; /*Bitmap indicator ( see Code Table 6.0 )*/

   /*call g2_create function to encodes Section 0 and Section 1
    * at the beginning of the message*/
   int status = g2_create(cgrib, ListSec0, ListSec1);
   if (status == -1)
   {
      param::instance << getTimeString() << "Write GRIB2 ERROR." << endl;
      exit(2); /* error return for GRIB message create*/
   }

   /*call g2_addgrid function to encodes a grid definition into Section 3*/
   status = g2_addgrid(cgrib, igds, igdstmpl, NULL, 0);
   if (status < 0)
   {
      param::instance << getTimeString() << "Write GRIB2 ERROR." << endl;
      exit(2); /* error return for grid add*/
   }

   g2float *fld = LST_grid.GetData();
   /*call g2_addfield to encode data field*/
   status = g2_addfield(cgrib, ipdsnum, ipdstmpl, coordlist, numcoord, idrsnum,
         idrstmpl, fld, ngrdpts, ibmap, NULL);
   if (status < 0)
   {
      param::instance << getTimeString() << "Write GRIB2 ERROR." << endl;
      exit(2); /* error return for field add*/
   }

   g2float *fld_qc = QC_Matrix_grid.GetData();
   /*call g2_addfield to encode data field*/
   status = g2_addfield(cgrib, ipdsnum, ipdstmpl_qc, coordlist, numcoord,
         idrsnum, idrstmpl_qc, fld_qc, ngrdpts, ibmap, NULL);
   if (status < 0)
   {
      param::instance << getTimeString() << "Write GRIB2 ERROR." << endl;
      exit(2); /* error return for field add*/
   }

   /*call g2_gribend add the final section 8 to the message*/
   status = g2_gribend(cgrib);
   if (status < 0)
   {
      param::instance << getTimeString() << "Write GRIB2 ERROR." << endl;
      exit(2); /* error return for end GRIB message*/
   }

   FILE * GribFile = fopen((Result_FN + ".grb").c_str(), "w");
   fwrite(cgrib, sizeof(unsigned char), status, GribFile);
   fclose(GribFile);

   if (cgrib != NULL)
   {
      delete[] cgrib;
   }

   return true; /*success GRIB format output reached*/
}

bool LST_Retrieval::Generate_QC_Metadata()
{

   /*Diagnostic code begin*/
   CMatrix<char> comb_filter_clearland;
   CMatrix<char> comb_filter_goodinput;
   char comb_filter_clearland_pixel = 0;
   char comb_filter_goodinput_pixel = 0;
   CMatrix<char> byte1_bit_2_3;
   CMatrix<char> byte1_bit_4_5;
   CMatrix<char> byte1_bit_6_7;
   CMatrix<char> byte2_bit_0_1;
   CMatrix<char> byte2_bit_2;
   CMatrix<char> byte2_bit_3;
   CMatrix<char> byte2_bit_4_5;
   CMatrix<char> byte2_bit_6_7;
   FMatrix FilteredBT39;
   FMatrix FilteredBT11;
   FMatrix FilteredSOZ;
   FMatrix FilteredSTZ;
   FMatrix FilteredSnow;
   FMatrix FilteredTPW;
   FMatrix FilteredEmTP2;
   FMatrix FilteredEmTP4;
   CMatrix<char> CombinedFilterMask;
   char CombinedFilterMaskPixel = 0;
   float em_tp2_cpixel = 0.0;
   float em_tp4_cpixel = 0.0;
   /*Diagnostic code end*/

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      comb_filter_clearland.Resize(ImgHeight, ImgWidth);
      comb_filter_goodinput.Resize(ImgHeight, ImgWidth);
      byte1_bit_2_3.Resize(ImgHeight, ImgWidth);
      byte1_bit_4_5.Resize(ImgHeight, ImgWidth);
      byte1_bit_6_7.Resize(ImgHeight, ImgWidth);
      byte2_bit_0_1.Resize(ImgHeight, ImgWidth);
      byte2_bit_2.Resize(ImgHeight, ImgWidth);
      byte2_bit_3.Resize(ImgHeight, ImgWidth);
      byte2_bit_4_5.Resize(ImgHeight, ImgWidth);
      byte2_bit_6_7.Resize(ImgHeight, ImgWidth);
      FilteredBT39.Resize(ImgHeight, ImgWidth);
      FilteredBT11.Resize(ImgHeight, ImgWidth);
      FilteredSOZ.Resize(ImgHeight, ImgWidth);
      FilteredSTZ.Resize(ImgHeight, ImgWidth);
      FilteredSnow.Resize(ImgHeight, ImgWidth);
      FilteredTPW.Resize(ImgHeight, ImgWidth);
      FilteredEmTP2.Resize(ImgHeight, ImgWidth);
      FilteredEmTP4.Resize(ImgHeight, ImgWidth);
      CombinedFilterMask.Resize(ImgHeight, ImgWidth);
      /*Diagnostic code end*/
   }

   /*Quality control flags*/
   float tp2_cpixel = 0.0;
   float tp4_cpixel = 0.0;
   float stz_cpixel = 0.0;
   float soz_cpixel = 0.0;
   float tpw_cpixel = 0.0;
   float snow_cpixel = 0.0;
   float lst_cpixel = 0.0;
   int32 IsLand = 0;
   char cloud_mask_cpixel = 0;
   short qc_bit = 0;
   short qc_value = 0;
   QC_Matrix.Resize(ImgHeight, ImgWidth);
   /*initiation of statistics of LST specific metadata*/
   lst_min = param::instance.THRESHOLD_LST_NORMAL_UPPER;
   lst_max = param::instance.THRESHOLD_LST_NORMAL_LOW;
   lst_sum = 0.0;
   lst_squre_sum = 0.0;
   lst_mean = 0.0;
   lst_std = 0.0;
   good_input_percent = 0.0;
   Num_Land = 0;
   Num_Good_Input = 0;
   Num_Cloud_Free = 0;
   Num_Snow = 0;
   Num_LST_Normal = 0;
   Num_LST_Bad = 0;
   Num_Large_View_Angle = 0;
   Num_Cold_Surface = 0;

   int32 count_pixels = 0;

   for (int32 iHeight = 0; iHeight < ImgHeight; iHeight++)
   {
      for (int32 iWidth = 0; iWidth < ImgWidth; iWidth++)
      {
         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            comb_filter_clearland_pixel = 0;
            comb_filter_goodinput_pixel = 0;
            CombinedFilterMaskPixel = 0;
            /*Diagnostic code end*/
         }

         QC_Matrix.SetElement(iHeight, iWidth,
               param::instance.INVALID_VALUE_INT);
         IsLand = LandMask_Img.GetElement(iHeight, iWidth);
         cloud_mask_cpixel = Cloud_Mask_Img.GetElement(iHeight, iWidth);
         tp2_cpixel = GOES_TP2_Aft_Cal.GetElement(iHeight, iWidth);
         tp4_cpixel = GOES_TP4_Aft_Cal.GetElement(iHeight, iWidth);
         stz_cpixel = Stz_Img.GetElement(iHeight, iWidth);
         soz_cpixel = Soz_Img.GetElement(iHeight, iWidth);
         tpw_cpixel = TPW_Img.GetElement(iHeight, iWidth);
         snow_cpixel = Snowmask_Img.GetElement(iHeight, iWidth);
         lst_cpixel = LST_before_scale.GetElement(iHeight, iWidth);

         /*LST Specific Metadata statistic*/
         if (IsLand != 0)
         {
            Num_Land++;

            if ((cloud_mask_cpixel == param::instance.THRESHOLD_CLOUD_CLEAR)
                  || (cloud_mask_cpixel
                        == param::instance.THRESHOLD_CLOUD_PRO_CLEAR))
            {
               Num_Cloud_Free++;
               if (param::instance.IsTest)
               {
                  /*Diagnostic code begin*/
                  CombinedFilterMaskPixel = 1;
                  FilteredBT39.SetElement(iHeight, iWidth, tp2_cpixel);
                  FilteredBT11.SetElement(iHeight, iWidth, tp4_cpixel);
                  FilteredSOZ.SetElement(iHeight, iWidth, soz_cpixel);
                  FilteredSTZ.SetElement(iHeight, iWidth, stz_cpixel);
                  FilteredSnow.SetElement(iHeight, iWidth, snow_cpixel);
                  FilteredTPW.SetElement(iHeight, iWidth, tpw_cpixel);
                  em_tp2_cpixel = Em_TP2_Resampled.GetElement(iHeight, iWidth);
                  em_tp4_cpixel = Em_TP4_Resampled.GetElement(iHeight, iWidth);
                  FilteredEmTP2.SetElement(iHeight, iWidth, em_tp2_cpixel);
                  FilteredEmTP4.SetElement(iHeight, iWidth, em_tp4_cpixel);
                  CombinedFilterMask.SetElement(iHeight, iWidth,
                        CombinedFilterMaskPixel);
                  /*Diagnostic code end*/
               }

            }

            if (snow_cpixel >= param::instance.THRESHOLD_SNOW)
               Num_Snow++;

            if (((cloud_mask_cpixel == param::instance.THRESHOLD_CLOUD_CLEAR)
                  || (cloud_mask_cpixel
                        == param::instance.THRESHOLD_CLOUD_PRO_CLEAR))
                  && (snow_cpixel < param::instance.THRESHOLD_SNOW))
            {
               if (param::instance.IsTest)
               {
                  /*Diagnostic code begin*/
                  comb_filter_clearland_pixel = 1;
                  comb_filter_clearland.SetElement(iHeight, iWidth,
                        comb_filter_clearland_pixel);
                  /*Diagnostic code end*/
               }
               if ((stz_cpixel < param::instance.THRESHOLD_VIEW_ANGLE))
               {
                  if (param::instance.IsTest)
                  {
                     /*Diagnostic code begin*/
                     comb_filter_goodinput_pixel = 1;
                     comb_filter_goodinput.SetElement(iHeight, iWidth,
                           comb_filter_goodinput_pixel);
                     /*Diagnostic code end*/
                  }

                  Num_Good_Input++;

                  if ((lst_cpixel >= param::instance.THRESHOLD_LST_NORMAL_LOW)
                        && (lst_cpixel
                              <= param::instance.THRESHOLD_LST_NORMAL_UPPER))
                  {
                     Num_LST_Normal++;
                  } else if ((lst_cpixel
                        >= param::instance.THRESHOLD_LST_COLD_SURFACE_LOW)
                        && (lst_cpixel
                              <= param::instance.THRESHOLD_LST_COLD_SURFACE_UPPER))
                  {
                     Num_Cold_Surface++;
                  } else if (lst_cpixel != param::instance.INVALID_VALUE_DOUBLE)
                     Num_LST_Bad++;
               } else
                  Num_Large_View_Angle++;
            }

            if (lst_cpixel != param::instance.INVALID_VALUE_DOUBLE)
            {
               count_pixels++;
               lst_sum += lst_cpixel;
               lst_squre_sum += lst_cpixel * lst_cpixel;
               if (lst_min > lst_cpixel)
               {
                  lst_min = lst_cpixel;
               }
               if (lst_max < lst_cpixel)
               {
                  lst_max = lst_cpixel;
               }
            }

         }

         /*quality control flag at bit level*/
         qc_value = 0;
         /*Availability*/
         if ((tp2_cpixel == param::instance.INVALID_VALUE_DOUBLE)
               || (tp4_cpixel == param::instance.INVALID_VALUE_DOUBLE))
         {
            qc_bit = 2; /*missing data*/
         } else if ((tp2_cpixel >= param::instance.THRESHOLD_GSIP_BT_NORMAL_LOW)
               && (tp2_cpixel <= param::instance.THRESHOLD_GSIP_BT_NORMAL_UPPER)
               && (tp4_cpixel >= param::instance.THRESHOLD_GSIP_BT_NORMAL_LOW)
               && (tp4_cpixel <= param::instance.THRESHOLD_GSIP_BT_NORMAL_UPPER))
         {
            qc_bit = 0; /*normal data*/
         } else
            qc_bit = 1; /*bad data*/
         qc_value = qc_value | (qc_bit << 12);
         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte1_bit_2_3.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*land type*/
         if (IsLand == 0)
         {
            qc_bit = 1; /*not land*/
         } else if (IsLand == -1)
         {
            qc_bit = 2; /*out of space*/
         } else
            qc_bit = 0; /*land*/
         qc_value = qc_value | (qc_bit << 10);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte1_bit_4_5.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*Cloud Index*/
         if (cloud_mask_cpixel == param::instance.THRESHOLD_CLOUD_CLEAR)
         {
            qc_bit = 0; /*cloud free*/
         } else if (cloud_mask_cpixel
               == param::instance.THRESHOLD_CLOUD_PRO_CLEAR)
         {
            qc_bit = 1; /*probably clear*/
         } else if (cloud_mask_cpixel
               == param::instance.THRESHOLD_CLOUD_PRO_CLOUDY)
         {
            qc_bit = 2; /*probably cloudy*/
         } else if (param::instance.THRESHOLD_CLOUD_CLOUDY)
         {
            qc_bit = 3; /*cloudy*/
         }
         qc_value = qc_value | (qc_bit << 8);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte1_bit_6_7.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*Snow Index*/
         if ((snow_cpixel < param::instance.THRESHOLD_SNOW)
               && (snow_cpixel >= 0))
         {
            qc_bit = 0; /*snow free*/
         } else if ((snow_cpixel >= param::instance.THRESHOLD_SNOW)
               && (snow_cpixel <= 1))
         {
            qc_bit = 1; /*snow*/
         } else
            qc_bit = 2; /*filled value*/
         qc_value = qc_value | (qc_bit << 6);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte2_bit_0_1.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*day/night flag*/
         if ((soz_cpixel > 0) && (soz_cpixel
               <= param::instance.THRESHOLD_SOLAR_ZEN))
         {
            qc_bit = 0; /*day time*/
         } else
            qc_bit = 1; /*night time*/
         qc_value = qc_value | (qc_bit << 5);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte2_bit_2.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*View Angle*/
         if (stz_cpixel > param::instance.THRESHOLD_VIEW_ANGLE)
         {
            qc_bit = 1; /*large view angle*/
         } else
            qc_bit = 0; /*normal*/
         qc_value = qc_value | (qc_bit << 4);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte2_bit_3.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*Atmospheric condition*/
         if (tpw_cpixel == param::instance.INVALID_VALUE_DOUBLE)
         {
            qc_bit = 3; /*filled value*/
         } else if (tpw_cpixel <= param::instance.THRESHOLD_TPW_LOW)
         {
            qc_bit = 0; /*dry*/
         } else if ((tpw_cpixel > param::instance.THRESHOLD_TPW_LOW)
               && (tpw_cpixel < param::instance.THRESHOLD_TPW_UPPER))
         {
            qc_bit = 1; /*moist*/
         } else if (tpw_cpixel >= param::instance.THRESHOLD_TPW_UPPER)
         {
            qc_bit = 2; /*very moist*/
         }
         qc_value = qc_value | (qc_bit << 2);

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte2_bit_4_5.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         /*LST Quality*/
         if (lst_cpixel == param::instance.INVALID_VALUE_DOUBLE)
         {
            qc_bit = 3; /*filled value*/
         } else if ((lst_cpixel >= param::instance.THRESHOLD_LST_NORMAL_LOW)
               && (lst_cpixel <= param::instance.THRESHOLD_LST_NORMAL_UPPER))
         {
            qc_bit = 0; /*normal*/
         } else if ((lst_cpixel
               >= param::instance.THRESHOLD_LST_COLD_SURFACE_LOW)
               && (lst_cpixel
                     <= param::instance.THRESHOLD_LST_COLD_SURFACE_UPPER))
         {
            qc_bit = 2; /*cold surface*/
         } else
            qc_bit = 1; /*out of range*/
         qc_value = qc_value | qc_bit;

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            byte2_bit_6_7.SetElement(iHeight, iWidth, qc_bit);
            /*Diagnostic code end*/
         }

         QC_Matrix.SetElement(iHeight, iWidth, qc_value);
      }
   }

   /*LST Specific Metadata statistic*/
   good_input_percent = (float) Num_Good_Input / (float) Num_Land;
   lst_mean = lst_sum / count_pixels;
   lst_std = sqrt(
         lst_squre_sum / count_pixels - lst_sum * lst_sum / count_pixels
               / (count_pixels - 1));

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      /*test_file to output QC flag*/
      ofstream test_file;
      string temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_QC_flag.dat";
      test_file.open(temp_fl.c_str());
      QC_Matrix.MatrixWrite(test_file);
      test_file.close();

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_Comb_filter_clearland.dat";
      test_file.open(temp_fl.c_str());
      comb_filter_clearland.MatrixWrite(test_file);
      test_file.close();

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_comb_filter_goodinput.dat";
      test_file.open(temp_fl.c_str());
      comb_filter_goodinput.MatrixWrite(test_file);
      test_file.close();

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_QC_decompose.dat";
      test_file.open(temp_fl.c_str());
      byte1_bit_2_3.MatrixWrite(test_file);
      byte1_bit_4_5.MatrixWrite(test_file);
      byte1_bit_6_7.MatrixWrite(test_file);
      byte2_bit_0_1.MatrixWrite(test_file);
      byte2_bit_2.MatrixWrite(test_file);
      byte2_bit_3.MatrixWrite(test_file);
      byte2_bit_4_5.MatrixWrite(test_file);
      byte2_bit_6_7.MatrixWrite(test_file);
      test_file.close();

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_Filtered_dataset.dat";
      test_file.open(temp_fl.c_str());
      FilteredBT39.MatrixWrite(test_file);
      FilteredBT11.MatrixWrite(test_file);
      FilteredSOZ.MatrixWrite(test_file);
      FilteredSTZ.MatrixWrite(test_file);
      FilteredSnow.MatrixWrite(test_file);
      FilteredTPW.MatrixWrite(test_file);
      FilteredEmTP2.MatrixWrite(test_file);
      FilteredEmTP4.MatrixWrite(test_file);
      test_file.close();

      temp_fl = getEnvDefine("GLST_OUTPUT_DIR") + "UT_Combined_FilterMask.dat";
      test_file.open(temp_fl.c_str());
      CombinedFilterMask.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/
   }

   return true; /* success metadata generation reached */
}

/*Retrieve LST for GOES 8-11*/
void LST_Retrieval::Get_LST_GOES8_11()
{

   if (GOES_ID == param::instance.GOES_ID_13)
   {
      param::instance << getTimeString() << "It is not an operational call."
            << endl;
      exit(2); /* error return for operational call*/
   }
   /////////////////////////////////////
   /*
    * LST Algorithm Coefficients.
    * Note that the coeffs are stratified in difference atmospheric conditions:
    * Daytime (solar zenith < 85 deg) and nighttime,
    * dry (atmospheric total column water vapor <=2.0 g/cm2)  and moist atmospheric conditions.
    */
   float coef_param[4][5]; /* Row0:Day Dry; Row1: Day Moist; Row2:Night Dry; Row3:Night Moist*/
   memset(coef_param, 0, sizeof(float) * 20);

   if (GOES_ID == param::instance.GOES_ID_8)
   {
      float coef_param_temp[4][5] =
      {
      { param::instance.COEF_GOES_8_0_0 , param::instance.COEF_GOES_8_0_1 ,
            param::instance.COEF_GOES_8_0_2 , param::instance.COEF_GOES_8_0_3 ,
            param::instance.COEF_GOES_8_0_4 } ,
      { param::instance.COEF_GOES_8_1_0 , param::instance.COEF_GOES_8_1_1 ,
            param::instance.COEF_GOES_8_1_2 , param::instance.COEF_GOES_8_1_3 ,
            param::instance.COEF_GOES_8_1_4 } ,
      { param::instance.COEF_GOES_8_2_0 , param::instance.COEF_GOES_8_2_1 ,
            param::instance.COEF_GOES_8_2_2 , param::instance.COEF_GOES_8_2_3 ,
            param::instance.COEF_GOES_8_2_4 } ,
      { param::instance.COEF_GOES_8_3_0 , param::instance.COEF_GOES_8_3_1 ,
            param::instance.COEF_GOES_8_3_2 , param::instance.COEF_GOES_8_3_3 ,
            param::instance.COEF_GOES_8_3_4 } };
      memcpy(coef_param, coef_param_temp, sizeof(float) * 4 * 5);
   } else if (GOES_ID == param::instance.GOES_ID_10)
   {
      float coef_param_temp[4][5] =
      {
            { param::instance.COEF_GOES_10_0_0 ,
                  param::instance.COEF_GOES_10_0_1 ,
                  param::instance.COEF_GOES_10_0_2 ,
                  param::instance.COEF_GOES_10_0_3 ,
                  param::instance.COEF_GOES_10_0_4 } ,
            { param::instance.COEF_GOES_10_1_0 ,
                  param::instance.COEF_GOES_10_1_1 ,
                  param::instance.COEF_GOES_10_1_2 ,
                  param::instance.COEF_GOES_10_1_3 ,
                  param::instance.COEF_GOES_10_1_4 } ,
            { param::instance.COEF_GOES_10_2_0 ,
                  param::instance.COEF_GOES_10_2_1 ,
                  param::instance.COEF_GOES_10_2_2 ,
                  param::instance.COEF_GOES_10_2_3 ,
                  param::instance.COEF_GOES_10_2_4 } ,
            { param::instance.COEF_GOES_10_3_0 ,
                  param::instance.COEF_GOES_10_3_1 ,
                  param::instance.COEF_GOES_10_3_2 ,
                  param::instance.COEF_GOES_10_3_3 ,
                  param::instance.COEF_GOES_10_3_4 } };
      memcpy(coef_param, coef_param_temp, sizeof(float) * 4 * 5);
   } else if (GOES_ID == param::instance.GOES_ID_9)
   {
      float coef_param_temp[4][5] =
      {
      { param::instance.COEF_GOES_9_0_0 , param::instance.COEF_GOES_9_0_1 ,
            param::instance.COEF_GOES_9_0_2 , param::instance.COEF_GOES_9_0_3 ,
            param::instance.COEF_GOES_9_0_4 } ,
      { param::instance.COEF_GOES_9_1_0 , param::instance.COEF_GOES_9_1_1 ,
            param::instance.COEF_GOES_9_1_2 , param::instance.COEF_GOES_9_1_3 ,
            param::instance.COEF_GOES_9_1_4 } ,
      { param::instance.COEF_GOES_9_2_0 , param::instance.COEF_GOES_9_2_1 ,
            param::instance.COEF_GOES_9_2_2 , param::instance.COEF_GOES_9_2_3 ,
            param::instance.COEF_GOES_9_2_4 } ,
      { param::instance.COEF_GOES_9_3_0 , param::instance.COEF_GOES_9_3_1 ,
            param::instance.COEF_GOES_9_3_2 , param::instance.COEF_GOES_9_3_3 ,
            param::instance.COEF_GOES_9_3_4 } };
      memcpy(coef_param, coef_param_temp, sizeof(float) * 4 * 5);
   } else if (GOES_ID == param::instance.GOES_ID_11)
   {
      float coef_param_temp[4][5] =
      {
            { param::instance.COEF_GOES_11_0_0 ,
                  param::instance.COEF_GOES_11_0_1 ,
                  param::instance.COEF_GOES_11_0_2 ,
                  param::instance.COEF_GOES_11_0_3 ,
                  param::instance.COEF_GOES_11_0_4 } ,
            { param::instance.COEF_GOES_11_1_0 ,
                  param::instance.COEF_GOES_11_1_1 ,
                  param::instance.COEF_GOES_11_1_2 ,
                  param::instance.COEF_GOES_11_1_3 ,
                  param::instance.COEF_GOES_11_1_4 } ,
            { param::instance.COEF_GOES_11_2_0 ,
                  param::instance.COEF_GOES_11_2_1 ,
                  param::instance.COEF_GOES_11_2_2 ,
                  param::instance.COEF_GOES_11_2_3 ,
                  param::instance.COEF_GOES_11_2_4 } ,
            { param::instance.COEF_GOES_11_3_0 ,
                  param::instance.COEF_GOES_11_3_1 ,
                  param::instance.COEF_GOES_11_3_2 ,
                  param::instance.COEF_GOES_11_3_3 ,
                  param::instance.COEF_GOES_11_3_4 } };
      memcpy(coef_param, coef_param_temp, sizeof(float) * 4 * 5);
   }

   float tp4_cpixel = 0.0;
   float tp5_cpixel = 0.0;
   float stz_cpixel = 0.0;
   float soz_cpixel = 0.0;
   float em_tp4_cpixel = 0.0;
   float em_tp5_cpixel = 0.0;
   float tpw_pixel = 0.0;
   float snow_pixel = 0.0;

   float coef[5];
   memset(coef, 0, sizeof(float) * 5);
   char cloud_mask_pixel = static_cast<char> (0);

   float RetrivalTemp = 0.0;
   char ScaledValue = static_cast<char> ( 0.0 );
   LST_Retrievals.Resize(ImgHeight, ImgWidth);
   QC_Matrix.Resize(ImgHeight, ImgWidth);
   for (int32 iHeight = 0; iHeight < ImgHeight; iHeight++)
   {
      for (int32 iWidth = 0; iWidth < ImgWidth; iWidth++)
      {
         RetrivalTemp = param::instance.INVALID_VALUE_DOUBLE;
         ScaledValue = param::instance.INVALID_VALUE_INT;
         QC_Matrix.SetElement(iHeight, iWidth,
               param::instance.INVALID_VALUE_INT);
         int32 IsLand = LandMask_Img.GetElement(iHeight, iWidth);
         if (IsLand > 0)
         {
            tp4_cpixel = GOES_TP2_Aft_Cal.GetElement(iHeight, iWidth);
            tp5_cpixel = GOES_TP4_Aft_Cal.GetElement(iHeight, iWidth);
            stz_cpixel = Stz_Img.GetElement(iHeight, iWidth);
            soz_cpixel = Soz_Img.GetElement(iHeight, iWidth);
            cloud_mask_pixel = Cloud_Mask_Img.GetElement(iHeight, iWidth);
            tpw_pixel = TPW_Img.GetElement(iHeight, iWidth);
            snow_pixel = Snowmask_Img.GetElement(iHeight, iWidth);
            em_tp4_cpixel = Em_TP2_Resampled.GetElement(iHeight, iWidth);
            em_tp5_cpixel = Em_TP4_Resampled.GetElement(iHeight, iWidth);

            coef[0] = 1;
            coef[1] = tp4_cpixel;
            coef[2] = tp4_cpixel - tp5_cpixel;
            coef[3] = 0.5 * (em_tp4_cpixel + em_tp5_cpixel);
            coef[4] = (tp4_cpixel - tp5_cpixel) * (1 / (cos(
                  stz_cpixel * M_PI / 180)) - 1);

            if ((soz_cpixel > 0) && (soz_cpixel
                  <= param::instance.THRESHOLD_SOLAR_ZEN) && ((cloud_mask_pixel
                  == param::instance.THRESHOLD_CLOUD_CLEAR)
                  || (cloud_mask_pixel
                        == param::instance.THRESHOLD_CLOUD_PRO_CLEAR)))
            {
               if (tpw_pixel < param::instance.THRESHOLD_TPW_LOW)
               {
                  RetrivalTemp = (coef_param[0][0] * coef[0])
                        + (coef_param[0][1] * coef[1]) + (coef_param[0][2]
                        * coef[2]) + (coef_param[0][3] * coef[3])
                        + (coef_param[0][4] * coef[4]);
               } else
               {
                  RetrivalTemp = (coef_param[1][0] * coef[0])
                        + (coef_param[1][1] * coef[1]) + (coef_param[1][2]
                        * coef[2]) + (coef_param[1][3] * coef[3])
                        + (coef_param[1][4] * coef[4]);
               }
            }
            if (((soz_cpixel < 0) || (soz_cpixel
                  > param::instance.THRESHOLD_SOLAR_ZEN)) && ((cloud_mask_pixel
                  == param::instance.THRESHOLD_CLOUD_CLEAR)
                  || (cloud_mask_pixel
                        == param::instance.THRESHOLD_CLOUD_PRO_CLEAR)))
            {
               if (tpw_pixel < param::instance.THRESHOLD_TPW_LOW)
               {
                  RetrivalTemp = (coef_param[2][0] * coef[0])
                        + (coef_param[2][1] * coef[1]) + (coef_param[2][2]
                        * coef[2]) + (coef_param[2][3] * coef[3])
                        + (coef_param[2][4] * coef[4]);
               } else
               {
                  RetrivalTemp = (coef_param[3][0] * coef[0])
                        + (coef_param[3][1] * coef[1]) + (coef_param[3][2]
                        * coef[2]) + (coef_param[3][3] * coef[3])
                        + (coef_param[3][4] * coef[4]);
               }
            }
            if ((cloud_mask_pixel == 0) && (stz_cpixel < 45) && (tp4_cpixel
                  >= 240))
            {
               QC_Matrix.SetElement(iHeight, iWidth, 0);
            }
            if ((cloud_mask_pixel == 1) || ((stz_cpixel >= 45) && (stz_cpixel
                  <= 70 * M_PI / 180)) || (tp4_cpixel < 240))
            {
               QC_Matrix.SetElement(iHeight, iWidth, 1);
            }
            if ((cloud_mask_pixel > 1) || (stz_cpixel > 70) || (tp4_cpixel
                  < 240))
            {
               QC_Matrix.SetElement(iHeight, iWidth, 2);
            }

         }
         if (RetrivalTemp != param::instance.INVALID_VALUE_DOUBLE)
         {
            if (RetrivalTemp < param::instance.THRESHOLD_LST_COLD_SURFACE_LOW)
               RetrivalTemp = param::instance.THRESHOLD_LST_COLD_SURFACE_LOW;
            if (RetrivalTemp > param::instance.THRESHOLD_LST_NORMAL_UPPER)
               RetrivalTemp = param::instance.THRESHOLD_LST_NORMAL_UPPER;
            ScaledValue = static_cast<char> ( (RetrivalTemp
                  - param::instance.METADATA_LST_ATTRI_OFFSET)
                  / param::instance.METADATA_LST_ATTRI_FACTOR );
         }
         LST_Retrievals.SetElement(iHeight, iWidth, ScaledValue);
      }
   }

   ScaleLonLat();
   OutputNetCDF();
   param::instance << getTimeString() << "Successful completion." << endl;

   return; /* success LST retrieval reached */
}

/*output results in NetCDF format*/
void LST_Retrieval::OutputNetCDF()
{
   NcFile nc((Result_FN + ".nc").c_str(), NcFile::Replace);
   Error_NetCDF(
         nc.add_dim(param::instance.METADATA_DIM_NAME_LON_SIZE.c_str(),
               ImgWidth));
   Error_NetCDF(
         nc.add_dim(param::instance.METADATA_DIM_NAME_LAT_SIZE.c_str(),
               ImgHeight));
   long array_output[2] =
   { ImgHeight , ImgWidth };

   /*write LST layer*/
   NcDim* dimLat = Error_NetCDF(
         nc.get_dim(param::instance.METADATA_DIM_NAME_LAT_SIZE.c_str()));
   NcDim* dimLon = Error_NetCDF(
         nc.get_dim(param::instance.METADATA_DIM_NAME_LON_SIZE.c_str()));
   Error_NetCDF(
         nc.add_var(param::instance.METADATA_DIM_NAME_LST.c_str(), ncShort,
               dimLat, dimLon));
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> put(
         LST_Retrievals.GetData(), array_output);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> add_att(
         "Factor", param::instance.METADATA_LST_ATTRI_FACTOR);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> add_att(
         "Offset", param::instance.METADATA_LST_ATTRI_OFFSET);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> add_att(
         "LST Form", param::instance.METADATA_LST_ATTRI_FORM.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> add_att(
         "FillValue", param::instance.METADATA_LST_ATTRI_FILL_VALUE);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LST.c_str()))-> add_att(
         "NumberType", param::instance.METADATA_LST_ATTRI_NUMBER_TYPE.c_str());

   /*write Longitude & Latitude layer*/
   if (param::instance.IsTest)
   {
      /*diagnostic code for write unscale LON/LAT*/
      Error_NetCDF(nc.add_var("Longitude", ncFloat, dimLat, dimLon));
      Error_NetCDF(nc.get_var("Longitude"))->put(lon_degree.GetData(),
            array_output);
      Error_NetCDF(nc.add_var("Latitude", ncFloat, dimLat, dimLon));
      Error_NetCDF(nc.get_var("Latitude"))->put(lat_degree.GetData(),
            array_output);
      /*diagnostic code for write unscale LON/LAT*/
   }

   Error_NetCDF(
         nc.add_var(param::instance.METADATA_DIM_NAME_LON.c_str(), ncShort,
               dimLat, dimLon));
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->put(
         lon_degree_scaled.GetData(), array_output);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "Factor", param::instance.METADATA_LON_ATTRI_FACTOR);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "Offset", param::instance.METADATA_LON_ATTRI_OFFSET);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "Longitude Form", param::instance.METADATA_LON_ATTRI_FORM.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "FillValue", param::instance.METADATA_LON_ATTRI_FILL_VALUE);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "NumberType",
         param::instance.METADATA_LONLAT_ATTRI_NUMBER_TYPE.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LON.c_str()))->add_att(
         "UNITS", param::instance.METADATA_LONLAT_ATTRI_UNIT.c_str());

   nc.add_var(param::instance.METADATA_DIM_NAME_LAT.c_str(), ncShort, dimLat,
         dimLon);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->put(
         lat_degree_scaled.GetData(), array_output);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "Factor", param::instance.METADATA_LAT_ATTRI_FACTOR);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "Offset", param::instance.METADATA_LAT_ATTRI_OFFSET);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "Latitude Form", param::instance.METADATA_LAT_ATTRI_FORM.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "FillValue", param::instance.METADATA_LAT_ATTRI_FILL_VALUE);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "NumberType",
         param::instance.METADATA_LONLAT_ATTRI_NUMBER_TYPE.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_LAT.c_str()))->add_att(
         "UNITS", param::instance.METADATA_LONLAT_ATTRI_UNIT.c_str());

   /*write quality control layer*/
   nc.add_var(param::instance.METADATA_DIM_NAME_QC.c_str(), ncShort, dimLat,
         dimLon);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_QC.c_str()))->put(
         QC_Matrix.GetData(), array_output);

   /*write common metadata
    a.   Product name
    b.   Satellite: string
    c.   Instrument: string
    d.   Projection: string
    e.   ProductResolution (at nadir): 16-bit integer
    f.   DateTime :   32-bit integer,   units = ccyyddd
    g.   RowNumbers:  16-bit integer
    h.   ColumnNumbers :  16-bit integer
    i.   LatitudeRange
    j.   LongitudeRange
    k.   CentralLatLongInfor
    l.   ByteOrderInfo (leftmost/rightmost)
    m.   Product Version Number
    n.   Data Compression Type
    o.   Ancillary data to produce product
    p.   Production Location
    q.   Contact Information
    */
   char mystr[256];
   memset(mystr, 0, 256);
   sprintf(mystr, "GOES %d", GOES_ID);
   Error_NetCDF_b(
         nc.add_att("ProductName",
               param::instance.METADATA_COMMON_PRODUCT_NAME.c_str()));
   Error_NetCDF_b(nc.add_att("Satellite", mystr));
   Error_NetCDF_b(
         nc.add_att("Instrument",
               param::instance.METADATA_COMMON_INSTRUMENT.c_str()));
   Error_NetCDF_b(
         nc.add_att("Projection",
               param::instance.METADATA_COMMON_PROJECTION.c_str()));
   Error_NetCDF_b(
         nc.add_att("ProductResolution at nadia (km)",
               param::instance.METADATA_COMMON_RESOLUTION));
   Error_NetCDF_b(nc.add_att("Date", img_date));
   Error_NetCDF_b(nc.add_att("Date", img_time));
   Error_NetCDF_b(nc.add_att("RowNumbers", ImgHeight));
   Error_NetCDF_b(nc.add_att("ColumnNumbers", ImgWidth));
   GetBoundingBox(lon_min, lon_max, lat_min, lat_max, lon_central, lat_central);
   Error_NetCDF_b(nc.add_att("longitude min", lon_min));
   Error_NetCDF_b(nc.add_att("longitude max", lon_max));
   Error_NetCDF_b(nc.add_att("latitude min", lat_min));
   Error_NetCDF_b(nc.add_att("latitude max", lat_max));
   Error_NetCDF_b(nc.add_att("central lon", lon_central));
   Error_NetCDF_b(nc.add_att("central lat", lat_central));
   Error_NetCDF_b(
         nc.add_att("ByteOrderInfo",
               param::instance.METADATA_COMMON_BYTE_ORDER.c_str()));
   Error_NetCDF_b(
         nc.add_att("Product Version Number",
               param::instance.METADATA_COMMON_VERSION.c_str()));
   Error_NetCDF_b(
         nc.add_att("Data format",
               param::instance.METADATA_COMMON_DATA_FORMAT.c_str()));
   Error_NetCDF_b(
         nc.add_att("Ancillary data to produce product",
               param::instance.METADATA_COMMON_ANCILLARY.c_str()));
   Error_NetCDF_b(
         nc.add_att("Production Location",
               param::instance.METADATA_COMMON_PRODUCTION_LOCATION.c_str()));
   Error_NetCDF_b(
         nc.add_att("Contact Information",
               param::instance.METADATA_COMMON_CONTACT.c_str()));

   /*write LST specific metadata
    a.   Min LST
    b.   Max LST
    c.   Mean LST
    d.   Standard Deviation
    e.   Percent of pixels with good input data
    f.   Total land pixels
    g.   Total numbers of cloud free pixels (land only)
    h.   Total numbers of snow pixels (land only)
    i.   Total number of good LST retrievals (land, cloud free, snow free, aerosol free, normal view angle, normal LST quality)
    j.   Total numbers of LST retrievals out of range (land, cloud free ,snow free, normal view angle and aerosol free)
    k.   Total numbers of LST retrievals with large view angle ( land, cloud free ,snow free and aerosol free)
    l.   Total numbers of LST retrievals with extreme cold surface ( land, cloud free , snow free and aerosol free)
    */
   Error_NetCDF_b(nc.add_att("Product Unit", "Kelvin"));
   Error_NetCDF_b(nc.add_att("Min LST", lst_min));
   Error_NetCDF_b(nc.add_att("Max LST", lst_max));
   Error_NetCDF_b(nc.add_att("Mean LST", lst_mean));
   Error_NetCDF_b(nc.add_att("Standard Deviation", lst_std));
   Error_NetCDF_b(
         nc.add_att("Percent of pixels with good input data",
               good_input_percent));
   Error_NetCDF_b(nc.add_att("Total land pixels", Num_Land));
   Error_NetCDF_b(
         nc.add_att("Total numbers of cloud free pixels", Num_Cloud_Free));
   Error_NetCDF_b(nc.add_att("Total numbers of snow pixels", Num_Snow));
   Error_NetCDF_b(
         nc.add_att("Total number of good LST retrievals", Num_LST_Normal));
   Error_NetCDF_b(
         nc.add_att("Total numbers of LST retrievals out of range", Num_LST_Bad));
   Error_NetCDF_b(
         nc.add_att("Total numbers of LST retrievals with large view angle",
               Num_Large_View_Angle));
   Error_NetCDF_b(
         nc.add_att(
               "Total numbers of LST retrievals with extreme cold surface",
               Num_Cold_Surface));

   /*image date*/
   nc.add_var(param::instance.METADATA_DIM_NAME_DATE.c_str(), ncInt);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_DATE.c_str()))->put(
         &img_date);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_DATE.c_str()))->add_att(
         "long_name", param::instance.METADATA_DATE_ATTRI_LONG_NAME.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_DATE.c_str()))->add_att(
         "units", param::instance.METADATA_DATE_ATTRI_UNIT.c_str());

   /*image time*/
   nc.add_var(param::instance.METADATA_DIM_NAME_TIME.c_str(), ncInt);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_TIME.c_str()))->put(
         &img_time);
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_TIME.c_str()))->add_att(
         "long_name", param::instance.METADATA_TIME_ATTRI_LONG_NAME.c_str());
   Error_NetCDF(nc.get_var(param::instance.METADATA_DIM_NAME_TIME.c_str()))->add_att(
         "units", param::instance.METADATA_TIME_ATTRI_UNIT.c_str());

   return; /* success NetCDF write reached */
}

/*scale longitude & latitude to 16-bit int32*/
void LST_Retrieval::ScaleLonLat()
{
   lon_degree_scaled.Resize(ImgHeight, ImgWidth);
   lat_degree_scaled.Resize(ImgHeight, ImgWidth);
   float beforeScale = 0.0;
   unsigned short afterScale = 0;
   short afterscale_short = 0;
   for (int32 iHeight = 0; iHeight < ImgHeight; iHeight++)
   {
      for (int32 iWidth = 0; iWidth < ImgWidth; iWidth++)
      {
         beforeScale = lon_degree.GetElement(iHeight, iWidth);
         if (beforeScale == param::instance.INVALID_VALUE_DOUBLE)
         {
            afterScale = param::instance.METADATA_LON_ATTRI_FILL_VALUE;
            afterscale_short = ConvertUStoS(afterScale);
            lon_degree_scaled.SetElement(iHeight, iWidth, afterscale_short);
         } else
         {
            afterScale = (unsigned short) ((beforeScale
                  - param::instance.METADATA_LON_ATTRI_OFFSET)
                  / param::instance.METADATA_LON_ATTRI_FACTOR);
            afterscale_short = ConvertUStoS(afterScale);
            lon_degree_scaled.SetElement(iHeight, iWidth, afterscale_short);
         }
         beforeScale = lat_degree.GetElement(iHeight, iWidth);
         if (beforeScale == param::instance.INVALID_VALUE_DOUBLE)
         {
            afterScale = param::instance.METADATA_LAT_ATTRI_FILL_VALUE;
            afterscale_short = ConvertUStoS(afterScale);
            lat_degree_scaled.SetElement(iHeight, iWidth, afterscale_short);
         } else
         {
            afterScale = (unsigned short) ((beforeScale
                  - param::instance.METADATA_LAT_ATTRI_OFFSET)
                  / param::instance.METADATA_LAT_ATTRI_FACTOR);
            afterscale_short = ConvertUStoS(afterScale);
            lat_degree_scaled.SetElement(iHeight, iWidth, afterscale_short);
         }
      }
   }

   if (param::instance.IsTest)
   {
      /*Diagnostic code begin*/
      ofstream test_file;

      string test_file_name = getEnvDefine("GLST_OUTPUT_DIR")
            + "UT_LonLat_scaled";
      test_file.open(test_file_name.c_str());
      lon_degree_scaled.MatrixWrite(test_file);
      lat_degree_scaled.MatrixWrite(test_file);
      test_file.close();
      /*Diagnostic code end*/

   }

   return; /* success LON/LAT scaling reached */
}

/*convert day of year to month and day*/
void LST_Retrieval::DayOfYear2MonthDay(int year,
      int dayofyear,
      int& nmonth,
      int& nday)
{
   int days_in_prev_months_u[13] =
   { 0 , 31 , 59 , 90 , 120 , 151 , 181 , 212 , 243 , 273 , 304 , 334 , 365 };
   int days_in_prev_months_l[13] =
   { 0 , 31 , 60 , 91 , 121 , 152 , 182 , 213 , 244 , 274 , 305 , 335 , 366 };

   if ((year % 4) == 0)
   {
      for (int im = 0; im < 13; im++)
      {
         if ((dayofyear > days_in_prev_months_l[im] && (dayofyear
               <= days_in_prev_months_l[im + 1])))
         {
            nmonth = im + 1;
            nday = dayofyear - days_in_prev_months_l[im];
         }
      }
   } else
   {
      for (int im = 0; im < 13; im++)
      {
         if ((dayofyear > days_in_prev_months_u[im] && (dayofyear
               <= days_in_prev_months_l[im + 1])))
         {
            nmonth = im + 1;
            nday = dayofyear - days_in_prev_months_u[im];
         }
      }
   }
}

/*get bounding box*/
void LST_Retrieval::GetBoundingBox(float& lon_min,
      float& lon_max,
      float& lat_min,
      float& lat_max,
      float& lon_cen,
      float& lat_cen)
{
   GetMinMax(ImgHeight, ImgWidth, lon_min, lon_max, lon_degree,
         param::instance.METADATA_LON_ATTRI_OFFSET);
   GetMinMax(ImgHeight, ImgWidth, lat_min, lat_max, lat_degree,
         param::instance.METADATA_LAT_ATTRI_OFFSET);
   int central_lon = (int) (ImgHeight / 2);
   int central_lat = (int) (ImgWidth / 2);
   lon_cen = lon_degree.GetElement(central_lon, central_lat);
   lat_cen = lat_degree.GetElement(central_lon, central_lat);
}

void LST_Retrieval::GetMinMax(int height,
      int width,
      float& min,
      float& max,
      FMatrix& matrix,
      float FilledValue)
{
   float temp_value = 0.0;
   min = 10000.0;
   max = -10000.0;
   for (int i = 0; i < height; i++)
   {
      for (int j = 0; j < width; j++)
      {
         temp_value = matrix.GetElement(i, j);
         if ((temp_value != param::instance.INVALID_VALUE_DOUBLE)
               && (temp_value != FilledValue))
         {
            if (temp_value < min)
            {
               min = temp_value;
            }
            if (temp_value > max)
            {
               max = temp_value;
            }
         }
      }
   }

}

/*Get Grid LST*/
float LST_Retrieval::GetGridLSTPixel(FMatrix& matrix,
      int centerx,
      int centery,
      int size)
{

   int x = centerx - size / 2;
   int y = centery - size / 2;
   float sum = 0;
   int count = 0;
   float value = 0;

   int xsize = x + size;
   int ysize = y + size;
   for (int i = x; i < xsize; i++)
   {
      for (int j = y; j < ysize; j++)
      {
         if (i == centerx && j == centery)
         {
            continue;
         }
         if (i < 0 || j < 0 || i >= GridHeight || j >= GridWidth)
         {
            continue;
         }
         value = matrix.GetElement(i, j);
         if (value != param::instance.INVALID_VALUE_DOUBLE)
         {
            sum += value;
            count++;
         }
      }
   }
   if (count > 0)
   {
      return sum / count;
   } else
   {
      return param::instance.INVALID_VALUE_DOUBLE;
   }
}

/*Get Grid QC*/
float LST_Retrieval::GetGridQCPixel(FMatrix& matrix,
      int centerx,
      int centery,
      int size)
{
   int x = centerx - size / 2;
   int y = centery - size / 2;
   float value = 0;
   float minvalue = param::instance.INVALID_VALUE_DOUBLE;
   int distance = size / 2 + 1;
   int xsize = x + size;
   int ysize = y + size;
   for (int i = x; i < xsize; i++)
   {
      for (int j = y; j < ysize; j++)
      {
         if (i == centerx && j == centery)
         {
            continue;
         }
         if (i < 0 || j < 0 || i >= GridHeight || j >= GridWidth)
         {
            continue;
         }
         value = matrix.GetElement(i, j);
         if (value != param::instance.INVALID_VALUE_DOUBLE)
         {
            int dis = min(abs(centerx - i), abs(centery - j));
            if (dis < distance)
            {
               distance = dis;
               minvalue = value;
               if (dis == 1)
                  continue;
            }
         }
      }
   }
   return minvalue;
}

/*Get Grid LST and QC*/
bool LST_Retrieval::GetGridLST_QC()
{
   int32 SampleLon = 0;
   int32 SampleLat = 0;
   GridWidth  = (int32) ((lon_max - lon_min) / param::instance.GRID_RESOLUTION + 1 );
   GridHeight = (int32) ((lat_max - lat_min) / param::instance.GRID_RESOLUTION + 1 );
   float cur_lon = 0.0;
   float cur_lat = 0.0;
   float cur_lst = 0.0;
   float cur_qc = 0;
   LST_grid.Resize(GridHeight, GridWidth);
   FMatrix LST_grid_temp(GridHeight, GridWidth);
   QC_Matrix_grid.Resize(GridHeight, GridWidth);
   FMatrix QC_Matrix_temp(GridHeight, GridWidth);

   for (int i = 0; i < GridHeight; i++)
   {
      for (int j = 0; j < GridWidth; j++)
      {
         LST_grid_temp.SetElement(i, j, param::instance.INVALID_VALUE_DOUBLE);
         LST_grid.SetElement(i, j, param::instance.INVALID_VALUE_DOUBLE);
         QC_Matrix_grid.SetElement(i, j, param::instance.INVALID_VALUE_DOUBLE);
         QC_Matrix_temp.SetElement(i, j, param::instance.INVALID_VALUE_DOUBLE);
      }
   }

   for (int i = 0; i < ImgHeight; i++)
   {
      for (int j = 0; j < ImgWidth; j++)
      {
         cur_lon = lon_degree.GetElement(i, j);
         cur_lat = lat_degree.GetElement(i, j);
         cur_lst = LST_before_scale.GetElement(i, j);
         cur_qc = (float) (QC_Matrix.GetElement(i, j));
         SampleLon = (int32) floor(
               (cur_lon - lon_min) / param::instance.GRID_RESOLUTION);
/* 
         SampleLat = (int32) floor(
              (lat_max - cur_lat) / param::instance.GRID_RESOLUTION);
*/
         /* modified by yuling, to make the lattitude dimension varies from south to north */
         SampleLat = (int32) floor(
                (cur_lat - lat_min) / param::instance.GRID_RESOLUTION); 

         if ((SampleLon >= 0) && (SampleLon < GridWidth) && (SampleLat >= 0)
               && (SampleLat < GridHeight))
         {
            LST_grid_temp.SetElement(SampleLat, SampleLon, cur_lst);
            QC_Matrix_temp.SetElement(SampleLat, SampleLon, cur_qc);
         }
      }
   }

   for (int i = 1; i < GridHeight - 1; i++)
   {
      for (int j = 1; j < GridWidth - 1; j++)
      {
         cur_lst = LST_grid_temp.GetElement(i, j);
         cur_qc = QC_Matrix_temp.GetElement(i, j);
         if (cur_lst != param::instance.INVALID_VALUE_DOUBLE)
         {
            LST_grid.SetElement(i, j, cur_lst);
            QC_Matrix_grid.SetElement(i, j, cur_qc);
         } else
         {
            LST_grid.SetElement(i, j, GetGridLSTPixel(LST_grid_temp, i, j, 5));
            QC_Matrix_grid.SetElement(i, j,
                  GetGridQCPixel(QC_Matrix_temp, i, j, 11));
         }
      }
   }
   return true; /* success grid LST reached */
}

short LST_Retrieval::ConvertUStoS(unsigned short data)
{
   short result;
   result = *((short *) &data);
   return result;
}




/*Read in GSIP L3 product
 IN: GSIP data file name;
 IN: Dim Name;
 IN: Array to read in the variable
     Latitude, Longitude
     Brightness Temperature at 3.9 and 11 micron
     Sensor Zenith Angle
     Solar Zenith Angle
     Total Precipitable Water
     Snow Fraction 
 Modified by Yuling for GSIP V3 update
 */
bool LST_Retrieval::Get_GSIP_V3_Float(string file_name,
      string DimName,
      FMatrix& ScaledData)
{

   NcFile nc(file_name.c_str(), NcFile::ReadOnly);
   short tempvalue = 0;
   float tempvalue2 = 0.0;
   if (nc.is_valid() == FALSE)
   {
      param::instance << getTimeString() << " NetCDF file open fail." << endl;
      exit(2); /* error return */
   }
 
   NcAtt* att= Error_NetCDF(nc.get_att("number_lines"));
   NcValues* ncvalue= att->values();
   ImgHeight = ncvalue->as_int(0);
   att= Error_NetCDF(nc.get_att("number_elements"));
   ncvalue= att->values();
   ImgWidth = ncvalue->as_int(0);

 
   NcVar* ncvar = Error_NetCDF(nc.get_var(DimName.c_str()));
   att = Error_NetCDF(ncvar->get_att("SCALED"));
   ncvalue= att->values();
   int32 if_scale =  ncvalue->as_int(0);

   ScaledData.Resize(ImgHeight, ImgWidth);
   float Unscaled_Value = 0.0;
 
   long array[] =
   { ImgHeight , ImgWidth };
  
   if ((DimName ==  param::instance.DIM_NAME_PIXEL_LATITUDE) || 
       (DimName ==  param::instance.DIM_NAME_PIXEL_LONGITUDE) )
   {
	 FMatrix ch2_values(ImgHeight,ImgWidth);

	 Error_NetCDF_b(ncvar->get(ch2_values.GetData(), array));

        att = Error_NetCDF(ncvar->get_att("_FillValue"));
        ncvalue = att->values();
        float Fill_Value_f = ncvalue->as_float(0);

	 for (int32 i = 0; i < ImgHeight; i++)
	 {
	      for (int32 j = 0; j < ImgWidth; j++)
	      {
		  tempvalue2 = ch2_values.GetElement(i,j);

       	  if (param::instance.IsTest)
	         {
       	     /*Diagnostic code begin*/
	            /*Diagnostic code end*/
       	  }
	         if (tempvalue2  != Fill_Value_f )
		  {
 		       Unscaled_Value=tempvalue2 ;
		  }
	         else
		  {
       	       Unscaled_Value = param::instance.INVALID_VALUE_DOUBLE;
		  }
	         ScaledData.SetElement(i, j, Unscaled_Value);
	      }
	  }
   }
   else
   {	
	 SMatrix ch2_values(ImgHeight,ImgWidth);
	 Error_NetCDF_b(ncvar->get(ch2_values.GetData(), array));

        att = Error_NetCDF(ncvar->get_att("scale_factor"));
        ncvalue = att->values();
        float SCALE_FACTOR = ncvalue->as_float(0);

        att = Error_NetCDF(ncvar->get_att("add_offset"));
        ncvalue = att->values();
        float OFFSET = ncvalue->as_float(0);
         
        att = Error_NetCDF(ncvar->get_att("_FillValue"));
        ncvalue = att->values();
        int32 Fill_Value_i = ncvalue->as_int(0);
   for (int32 i = 0; i < ImgHeight; i++)
   {
      for (int32 j = 0; j < ImgWidth; j++)
      {

         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            /*Diagnostic code end*/
         }
         if ( if_scale == 1 )
         {
	    tempvalue =ch2_values.GetElement(i,j) ;
	    if ( tempvalue  != Fill_Value_i )
 	       Unscaled_Value=float(tempvalue)  * SCALE_FACTOR + OFFSET;
            else
              Unscaled_Value = param::instance.INVALID_VALUE_DOUBLE;
         }
         else
         {
	      Unscaled_Value =float(tempvalue);	
	 }     
         ScaledData.SetElement(i, j, Unscaled_Value);
      }
   }
   }

   return true; /* success GSIP product read reached */
}


/*Read in GSIP V3 product
 IN: GSIP data file name;
 IN: Dim Name;
 IN: Array to read in the variable with Byte type
     Cloud Mask
 Modified by Yuling Liu for GSIP V3 update 
 */
bool LST_Retrieval::Get_GSIP_V3_Byte(string file_name,
      string DimName,
      BMatrix& ScaledData)
{

   NcFile nc(file_name.c_str(), NcFile::ReadOnly);
   if (nc.is_valid() == FALSE)
   {
      param::instance << getTimeString() << " NetCDF file open fail." << endl;
      exit(2); /* error return */
   }
 
   NcVar* ncvar = Error_NetCDF(nc.get_var(DimName.c_str()));
   ScaledData.Resize(ImgHeight, ImgWidth);
   long array[] = { ImgHeight , ImgWidth };

   BMatrix origData(ImgHeight,ImgWidth);
   /* has to convert to ncbyte */
   Error_NetCDF_b(ncvar->get((ncbyte*)origData.GetData(), array));   

   NcAtt* att = Error_NetCDF(ncvar->get_att("_FillValue"));
   NcValues* ncvalue = att->values();
   char Fill_Value_c = ncvalue->as_ncbyte(0);

   char tempvalue = 0;   

   for (int32 i = 0; i < ImgHeight; i++)
   {
      for (int32 j = 0; j < ImgWidth; j++)
      {
	 tempvalue =origData.GetElement(i,j) ;
         if (param::instance.IsTest)
         {
            /*Diagnostic code begin*/
            /*Diagnostic code end*/
         }
         if ( tempvalue  == Fill_Value_c )
         {
 	    tempvalue  = param::instance.INVALID_VALUE_INT;	 
	  }
 
         ScaledData.SetElement(i, j, tempvalue);
      }
   }

   return true; /* success GSIP product read reached */
}