/home/pascal/depot/filedata/src/ecmwffiledata.h

/***************************************************************************
 *   Copyright (C) 2005 by Nicolas PASCAL   *
 *   nicolas.pascal@icare.univ-lille1.fr   *
 *                                                                         *
 *   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.             *
 ***************************************************************************/
#ifndef ECMWFFILEDATA_H
#define ECMWFFILEDATA_H

#include "gribfiledata.h"
#include "meteofiledata.h"
#include "interpolation.h"
#include "earth_geometry.h"

class ECMWFFileData : public MeteoFileData, public GribFileData, public FileDataReader {
protected:
    typedef enum {
        UNDEFINED = 0,
        ANALYSIS,
        FORECAST
    } ECMWFProductType;
    ECMWFProductType product_id;
    string version;
    int sz_x;
    int sz_y;
    float lat0;
    float lon0;

    void set_time_coverage();
    void parse_filename(const string & short_filename);
    bool check_filename(const string & short_filename) const;
    void init_grid ();
    virtual void load_time_level() {;};
public:
    static const unsigned short isobaric_level_pressures[];
    static const int nb_max_isobaric_level;
    static const double a_hybrid_91[92];
    static const double b_hybrid_91[92];
    static const int nb_level_hybrid_91;
    ECMWFFileData(const string &name, const string &mode="r");
    ~ECMWFFileData();
    void free_geolocation_data() {;};
    inline const bool get_index(const float &lat, const float& lon, const double &time, int * indexes);
    inline const bool get_index(const float &lat, const float& lon, int * indexes);
    inline const bool get_index(const float &lat, const float& lon, const float &p_lvl, int * indexes) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
    };
    inline const int get_lat_index(const float &lat) const;
    inline const int get_lon_index(const float &lon) const;

    inline const float32 get_lat (const int &i_y) const;
    inline const float32 get_lon (const int &i_x) const;
     void * read_data(void * data, const char * sds_name, int * start=NULL, int * stride=NULL, int *edges=NULL, int rank=-1) {
        return this->GribFileData::read_data(data, sds_name, start, stride, edges, rank);
    };
    ECMWFProductType get_product_id() const {
        return product_id;
    };
    const string get_s_product_id () const {
        string s_prod_id ("");
        switch (product_id) {
            case UNDEFINED:
                s_prod_id = "UNDEFINED";
                break;
            case ANALYSIS:
                s_prod_id = "ANALYSIS";
                break;
            case FORECAST:
                s_prod_id = "FORECAST";
                break;
            default :
                g_exception e(__FILE__ , __LINE__ , "Invalid product type ID " + s_prod_id);
                throw e;
        };
        return s_prod_id;
    };
    int get_sz_x() const {
        return sz_x;
    }
    int get_sz_y() const {
        return sz_y;
    }
    const string get_grid_desc ();
    const vector <int> get_pressure_levels(const string & var_name);
    const int get_nb_pressure_levels (const string & var_name);
    void get_hybrid_level_pressures(const double &p_surf,
        double * v_p_hybrid) const {
        for (int i = 0 ; i < nb_level_hybrid_91 ; ++i)
            v_p_hybrid[i] = get_hybrid_level_pressure (i, p_surf);
    };
    const double get_hybrid_level_pressure (const int i_level, const double &p_surf) const {
        if (i_level < 0 || i_level > nb_level_hybrid_91 ) {
            g_exception e(__FILE__ , __LINE__ , "Invalid hybrid level number " + MyTools::to_string(i_level) + " Must in range [0, 91]");
            throw e;
        }
// cout << "i_level=" << i_level << " a=[" << a_hybrid_91 [i_level] << "," << a_hybrid_91 [i_level+1] << "] b=[" << b_hybrid_91 [i_level] << "," << b_hybrid_91 [i_level+1] << "] psurf=" << p_surf << " -> " << .5 * ((a_hybrid_91 [i_level] + a_hybrid_91 [i_level+1]) * 0.01 + (b_hybrid_91 [i_level] + b_hybrid_91 [i_level+1]) * p_surf) << endl;
        return .5 * ((a_hybrid_91 [i_level] + a_hybrid_91 [i_level+1]) * 0.01 + (b_hybrid_91 [i_level] + b_hybrid_91 [i_level+1]) * p_surf);
    };
    static const unsigned short* get_isobaric_level_pressures() {
        return isobaric_level_pressures;
    }
    template <typename Txout, typename Tyout, typename Txin, typename Tyin>
    void get_v_pres_on_alt (
        Tyout * v_out_pres, const Txout * v_alt, const int sz_alt,
        const Txin * v_geop, const Tyin * v_inpres, const int n_pres,
        const string interp_type = "linear");
    template <typename Txout, typename Txin, typename Tyin>
    double get_pres_on_alt (
        const Txout alt,
        const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres,
        const string interp_type = "linear");
    double alt_to_geop (const double & alt) {
        return alt * Earth::g0;
    };
    template <typename Txout, typename Tyout>
    void get_abscissa_iso (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out_pres);

    template <typename Txout, typename Tyout, typename Txin>
    void get_abscissa_hybrid (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out_pres,
                              const Txin * v_in_pres, const int sz_in_pres);

    template <typename Txout, typename Tyout, typename Txin, typename Tyin>
    void interp_on_pres (Tyout * v_out_var, const Txout * v_out_abs, const int sz_out,
                         const Tyin * v_in_var, const Txin * v_in_pres, const int sz_in);

    void close_data_file();
    void open_data_file();
//     vector<int> get_dataset_dimension(const string &sds_name);
//     void get_dataset_fill_value(const string &name,void *fillvalue);
    int get_n_dataset();
    string get_dataset_name(int i);
    int get_dataset_data_type(string sds_name);
    string get_values_attr_dataset(string sds_name,string attr);
    bool has_attr_dataset(string sds_name,string attr);
};

const int ECMWFFileData::get_lat_index(const float & lat) const {
    if (lat < -90. || lat > 90.) {
        g_exception e(__FILE__ , __LINE__ , "Invalid latitude " + MyTools::to_string(lat));
        throw e;
    }
    return MyTools::round<int>(fabs(lat - lat0) / delta_lat);
}

const float32 ECMWFFileData::get_lat(const int & i_y) const {
    return lat0 - i_y * delta_lat;
}
const float32 ECMWFFileData::get_lon(const int & i_x) const {
    float32 lon = lon0 + i_x * delta_lon;
    if (lon > 180.)
        lon -= 360.;
    return lon;
}
const int ECMWFFileData::get_lon_index(const float & _lon) const {
    if (_lon < -180. || _lon > 180.) {
        g_exception e(__FILE__ , __LINE__ , "Invalid longitude " + MyTools::to_string(_lon));
        throw e;
    }
    // longitudes as [-180, 180] -> [0, 360[
    float lon (_lon);
    if (lon < 0.)
        lon = 360. + lon;
    int ilon = MyTools::round<int>((lon - lon0) / delta_lon);
    // treat extreme case of longitudes in range [-delta_lon/2, 0]
    if (ilon == sz_x)
        ilon = 0;
    return ilon;
}

const bool ECMWFFileData::get_index(const float & lat, const float & lon, int * indexes) {
    indexes [0] = get_lat_index(lat);
    indexes [1] = get_lon_index(lon);
    return true;
}

const bool ECMWFFileData::get_index(const float & lat, const float & lon, const double & time, int * indexes) {
    return get_index (lat, lon, indexes);
}

template <typename Txout, typename Tyout, typename Txin, typename Tyin>
void ECMWFFileData::get_v_pres_on_alt(
    Tyout * v_outpres, const Txout * v_alt, const int sz_alt,
    const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres,
    const string interp_type) {
    // extrapolation forbidden -> NaN returned if out of range
    double nan = numeric_limits<double>::signaling_NaN();

    if (interp_type == "linear" || interp_type == "linear_log") {
        for (int i = 0 ; i < sz_alt ; ++i)
            v_outpres [i] = get_pres_on_alt (v_alt[i], v_ingeop, v_inpres, sz_inpres, interp_type) ;
    } else if (interp_type == "spline") {
        // computations made in double. Last geopotential is ignored because set on a virtual pressure of 0 hPa. This is the reason of the "- 1" in size of the input pressures vector
        int sz_in = sz_inpres - 1;
        double x_in  [sz_in];
        copy (v_ingeop, v_ingeop + sz_in, x_in);
        double y_in  [sz_in];
        copy (v_inpres, v_inpres + sz_in, y_in);
        int sz_out = sz_alt;
        double x_out [sz_out];
        for (int i = 0 ; i < sz_out ; ++i)
            x_out[i] = alt_to_geop (v_alt[i] * 1000.); // * 1000 : km -> m;
        double y_out [sz_out];
        // compute spline interp
        // -------------------------------------------------------
        spline_interp (y_out, x_out, sz_out,
                       y_in,  x_in, sz_in, nan);

        copy (y_out, y_out + sz_out, v_outpres);

    } else {
        g_exception e(__FILE__ , __LINE__ , "Invalid interpolation type " + interp_type);
        throw e;
    }
}

template <typename Txout, typename Txin, typename Tyin>
double ECMWFFileData::get_pres_on_alt (
    const Txout alt,
    const Txin * v_ingeop, const Tyin * v_inpres, const int sz_inpres,
    const string interp_type) {
    // geopotential of searched altitude
    double outgeop = alt_to_geop (alt * 1000.); // * 1000 : km -> m;
    // extrapolation forbidden -> NaN returned if out of range
    double nan = numeric_limits<double>::signaling_NaN();
    // last geopotential is ignored because set on a virtual pressure of 0 hPa -> -2
    if (outgeop < v_ingeop[0] || outgeop > v_ingeop [sz_inpres - 2])
        return nan;
    // process to linear interpolation of pressure
    int i = 0;
    while ((i < (sz_inpres - 2)) && (outgeop > v_ingeop[i]))
        ++i;
//     if (i >= (sz_inpres - 2))
    if (i >= (sz_inpres - 1))
        return nan;
    else {
        double outpres = nan;

        if (interp_type == "linear") {
            outpres = v_inpres[i - 1] + (outgeop - v_ingeop[i - 1]) * (v_inpres[i] - v_inpres[i - 1])/(v_ingeop[i] - v_ingeop[i - 1]);
        } else if (interp_type == "linear_log") {
            outpres = log(v_inpres[i - 1]) + (log(outgeop) - log(v_ingeop[i - 1])) * (log(v_inpres[i]) - log(v_inpres[i - 1]))/(log(v_ingeop[i]) - log(v_ingeop[i - 1]));
            outpres = exp(outpres);
        } else {
            g_exception e(__FILE__ , __LINE__ , "Invalid interpolation type " + interp_type);
            throw e;
        }
        if (outpres < 0)
            outpres = nan;
// cout << alt << " " << outgeop << " v_ingeop=[" << v_ingeop[i - 1] << "," << v_ingeop[i] << "] v_inpres=[" << v_inpres[i - 1] << "," << v_inpres[i] << "] -> " << outpres << endl;
        return outpres;
    }
}

template <typename Txout, typename Tyout>
void ECMWFFileData::get_abscissa_iso(Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out) {
    MyTools::get_log_abscissa_sequential (v_out_abs, v_out_pres, sz_out, isobaric_level_pressures, nb_max_isobaric_level - 1);
}

template <typename Txout, typename Tyout, typename Txin>
void ECMWFFileData::get_abscissa_hybrid (Tyout * v_out_abs, const Txout * v_out_pres, const int sz_out, const Txin * v_in_pres, const int sz_in) {
    MyTools::get_log_abscissa_sequential (v_out_abs, v_out_pres, sz_out, v_in_pres, sz_in);
}

template <typename Txout, typename Tyout, typename Txin, typename Tyin>
void ECMWFFileData::interp_on_pres(Tyout * v_out_var, const Txout * v_out_abs, const int sz_out, const Tyin * v_in_var, const Txin * v_in_pres, const int sz_in) {
    // extrapolation forbidden -> NaN returned if out of range
    Tyout nan = numeric_limits<Tyout>::signaling_NaN();
    int i1, i2;
    double y1, y2;
    double abs;
    for (int i = 0 ; i < sz_out ; ++i) {
        abs = v_out_abs [i];
        if (isnan(abs))
            v_out_var [i] = nan;
        else {
            i1 = long(floor(abs));
            i2 = i1 + 1;
            y1 = v_in_var [i1];
            y2 = v_in_var [i2];
            v_out_var [i] = y1 + (abs - i1) * (y2 - y1);
        }
    }
}

#endif