/***************************************************************************
 *   Copyright (C) 2006 by Nicolas PASCAL   *
 *   pascal@icare-pc12   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "ecmwffiledata.h"
#include "caliopfiledata.h"

int main(int argc, char *argv[]) {
    /* Open and create the access to the file  */
//     string fname = "/DATA/LIENS/ECMWF/operational/2009/2009_01_01/ECMWF_200901010000_v1.0.0_1deg_AN.grb";
//     string fname = "/DATA/LIENS/ECMWF/operational.v1.2.0/2011/2011_05_17/ECMWF_201105170000_v1.2.0_1deg_AN.grb";
    string fname = "/DATA/LIENS/ECMWF/operational.v1.2.0/2011/2011_05_20/ECMWF_201105200000_v1.2.0_1deg_AN.grb";
    ECMWFFileData f (fname);

    /* print out some informations */
    printf ("--- %s ---\n", f.get_name().c_str());
    printf ("Product Type %s\n", f.get_s_product_id ().c_str());
    printf ("Start Acq. Time %s\n", f.get_date().to_string().c_str());
    printf ("Time coverage (sec) %f\n", f.get_time_coverage());
    printf ("- Grid -\n%s", f.get_grid_desc().c_str());

    vector <int> pressure_levels = f.get_pressure_levels("z");
    size_t nb_pressure_levels = pressure_levels.size();
    cout<<"Nb pressure levels "<<nb_pressure_levels<<endl;
    cout << "Isobaric Pressure levels : ";
    for (size_t ilvl = 0 ; ilvl < nb_pressure_levels ; ++ilvl)
        cout << pressure_levels [ilvl] << " ";
    cout << endl;

    /* pressure altitudes interpolation */
    const float32 * v_outalt =  CALIOPFileData::meteo_level_altitude;
    int sz_outalt = CALIOPFileData::nb_meteo_level;

    float32 v_geop_isobaric [] = {-551.095, 3510.92, 5605.53, 7746.42, 12183, 16843.4, 26935.3, 38263.4, 51209.3, 66342.5, 85050.8, 96733.9, 111039, 129511, 155535, 178568, 200329, 234182, 261483, 309208, 334397, 358841, 397801, 430071, 485990, 0.687012};
    const unsigned short * v_pres_isobaric = ECMWFFileData::isobaric_level_pressures;
    int sz_pres_isobaric = ECMWFFileData::nb_max_isobaric_level;
    double v_outpres [sz_outalt];
    f.get_v_pres_on_alt (v_outpres, v_outalt, sz_outalt,
                         v_geop_isobaric, v_pres_isobaric, sz_pres_isobaric);

    cout << "Output pressures : ";
    MyTools::print_vec ((double*)(v_outpres), sz_outalt);
    cout << "Output altitudes : ";
    MyTools::print_vec ((float32*)(v_outalt), sz_outalt);
    cout << "Input pressures : ";
    MyTools::print_vec ((const unsigned short*)(v_pres_isobaric), sz_pres_isobaric);
    cout << "Input Geopotential : ";
    MyTools::print_vec ((float32*)(v_geop_isobaric), sz_pres_isobaric);

    /* (lat,lon) Coincidence search */
    int ipix [] = {-1, -1};
    float lat = -90., lon = -180.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = -179.5;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = -179.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = -178.5;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = 178.5;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = 179.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = 179.5;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = -90., lon = -90.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = 0., lon = 0.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    lat = 90., lon = 180.;
    f.get_index (lat, lon, ipix);
    printf ("(%.1f, %.1f) -> [%d, %d]\n", lat, lon, ipix[0], ipix[1]);

    /* data extraction */
    f.open_data_file();

    // --- 2D
    int edges_2d [] = {1, 1};
    int rank;
    string var_name ("");
    double * data_2d = NULL;

    // unique grid cell
    printf ("--- read 2D dataset ---\n");
    var_name = "skt";
    rank = 2;

    lat = 0., lon = 0.;
    f.get_index (lat, lon, ipix);
    data_2d = (double*)f.read_data(data_2d, var_name.c_str(), ipix, NULL, edges_2d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = %.2f\n", var_name.c_str(), lat, lon, ipix[0], ipix[1], *data_2d);

    lat = -90., lon = -180.;
    f.get_index (lat, lon, ipix);
    data_2d = (double*)f.read_data(data_2d, var_name.c_str(), ipix, NULL, edges_2d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = %.2f\n", var_name.c_str(), lat, lon, ipix[0], ipix[1], *data_2d);

    lat = 90., lon = -180.;
    f.get_index (lat, lon, ipix);
    data_2d = (double*)f.read_data(data_2d, var_name.c_str(), ipix, NULL, edges_2d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = %.2f\n", var_name.c_str(), lat, lon, ipix[0], ipix[1], *data_2d);

    lat = 90., lon = 180.;
    f.get_index (lat, lon, ipix);
    data_2d = (double*)f.read_data(data_2d, var_name.c_str(), ipix, NULL, edges_2d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = %.2f\n", var_name.c_str(), lat, lon, ipix[0], ipix[1], *data_2d);

    lat = 17.3, lon = 106.6;
    f.get_index (lat, lon, ipix);
    data_2d = (double*)f.read_data(data_2d, var_name.c_str(), ipix, NULL, edges_2d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = %.2f\n", var_name.c_str(), lat, lon, ipix[0], ipix[1], *data_2d);

    delete[] data_2d, data_2d = NULL;

    // --- 3D
    int start_3d [] = {0, -1, -1};
    int edges_3d [] = {0, 1, 1}; // get all pressure profiles
    double * data_3d = NULL;

    // unique grid cell
    printf ("--- read 3D dataset ---\n");
    var_name = "u";
    rank = 3;
    edges_3d [0] = f.get_nb_pressure_levels(var_name);

    lat = 0., lon = 0.;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    lat = -90., lon = -180.;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    lat = 90., lon = -180.;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    lat = 90., lon = 180.;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    lat = -90., lon = 180.;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    lat = 17.3, lon = 106.6;
    f.get_index (lat, lon, start_3d + 1);
    data_3d = (double*)f.read_data(data_3d, var_name.c_str(), start_3d, NULL, edges_3d, rank);
    printf ("%s (%.1f, %.1f) [%d, %d] = ", var_name.c_str(), lat, lon, start_3d[1], start_3d[2]);
    for (int i = 0 ; i < edges_3d [0] ; ++i)
        printf ("%.2f ", data_3d[i]);
    printf ("\n");

    delete[] data_3d, data_3d = NULL;


    f.close_data_file();

    return 0;
}