/*
 * (C) Copyright 2005- ECMWF.
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 *
 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
 */

/* cmake config header */
#ifdef HAVE_ECCODES_CONFIG_H
#include "eccodes_config.h"
#endif

/* autoconf config header */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <stdio.h>
#include <stdlib.h>
#include <grib_api.h>
#include <errno.h>
#include <string.h>
#include <math.h>
#include <sys/types.h>

#ifndef ECCODES_ON_WINDOWS
  #include <dirent.h>
#else
  #include <direct.h>
  #include <io.h>

  #ifdef _MSC_VER
  #   define strcasecmp _stricmp
  #endif
#endif

#include <assert.h>
#include "tigge_tools.h"
void validate(const char* path);

#define CHECK(a)  check(#a,a)
#define WARN(a)   warn(#a,a)
#define NUMBER(a) (sizeof(a)/sizeof(a[0]))

typedef struct pair {
    const char* key;
    int         key_type;
    long        value_long;
    char*       value_string;
} pair;

typedef struct parameter parameter;

typedef void (*check_proc)(grib_handle*,const parameter*,double min,double max);

struct parameter {
    const char* name;
    double min1;
    double min2;
    double max1;
    double max2;
    pair   pairs[15];
    check_proc checks[5];
};

static void point_in_time(grib_handle*,const parameter*,double,double);
static void statistical_process(grib_handle*,const parameter*,double,double);
static void six_hourly(grib_handle*,const parameter*,double,double);
static void since_prev_pp(grib_handle*,const parameter*,double,double);
static void three_hourly(grib_handle* h,const parameter* p,double min,double max);
static void from_start(grib_handle*,const parameter*,double,double);
static void daily_average(grib_handle*,const parameter*,double,double);
static void given_level(grib_handle*,const parameter*,double,double);
static void predefined_level(grib_handle*,const parameter*,double,double);
static void predefined_thickness(grib_handle*,const parameter*,double,double);
static void given_thickness(grib_handle*,const parameter*,double,double);
static void has_bitmap(grib_handle*,const parameter*,double,double);
static void has_soil_level(grib_handle*,const parameter*,double,double);
static void has_soil_layer(grib_handle*,const parameter*,double,double);
static void resolution_s2s(grib_handle*,const parameter*,double,double);
static void resolution_s2s_ocean(grib_handle*,const parameter*,double,double);

static void height_level(grib_handle*,const parameter*,double,double);
static void pressure_level(grib_handle*,const parameter*,double,double);
static void potential_temperature_level(grib_handle*,const parameter*,double,double);
static void potential_vorticity_level(grib_handle*,const parameter*,double,double);

/*
TODO:
- Shape of the earth
- Levels
*/

#include "tigge_check.h"


int error = 0;
const char* file = 0;
int field = 0;
int warning = 0;
int valueflg = 0;
const char* param = "unknown";
int warnflg = 0;
int zeroflg = 0;
int is_lam = 0;
int is_s2s = 0;
int is_s2s_refcst = 0;
int is_uerra = 0;
int is_crra = 0;

const char* good = NULL;
const char* bad = NULL;

FILE* fgood = NULL;
FILE* fbad  = NULL;

static void check(const char* name,int a)
{
    if(!a) {
        printf("%s, field %d [%s]: %s failed\n",file,field,param,name);
        error++;
    }
}

/*
static void warn(const char* name,int a)
{
    if(!a) {
        printf("%s, field %d [%s]: %s failed\n",file,field,param,name);
        warning++;
    }
}
*/

static void save(grib_handle* h, const char *name,FILE* f)
{
    size_t size;
    const void *buffer;
    int e;

    if(!f) return;

    if((e = grib_get_message(h,&buffer,&size)) != GRIB_SUCCESS)
    {
        printf("%s, field %d [%s]: cannot get message: %s\n",file,field,param,grib_get_error_message(e));
        exit(1);
    }

    if(fwrite(buffer,1,size,f) != size)
    {
        perror(name);
        exit(1);
    }
}

static long get(grib_handle *h,const char* what)
{
    int e; long val;
    if((e = grib_get_long(h,what,&val)) != GRIB_SUCCESS)
    {
        printf("%s, field %d [%s]: cannot get %s: %s\n",file,field,param,what,grib_get_error_message(e));
        error++;
        val = -1;
    }
    return val;
}

static double dget(grib_handle *h,const char* what)
{
    int e; double val;
    if((e = grib_get_double(h,what,&val)) != GRIB_SUCCESS)
    {
        printf("%s, field %d [%s]: cannot get %s: %s\n",file,field,param,what,grib_get_error_message(e));
        error++;
        val = -1;
    }
    return val;
}

static int missing(grib_handle *h,const char* what)
{
    int err=0;
    return grib_is_missing(h,what,&err);
}

static int eq(grib_handle *h,const char* what,long value)
{
    return get(h,what) == value;
}

static int ne(grib_handle *h,const char* what,long value)
{
    return get(h,what) != value;
}

static int ge(grib_handle *h,const char* what,long value)
{
    return get(h,what) >= value;
}

static int le(grib_handle *h,const char* what,long value)
{
    return get(h,what) <= value;
}

static int DBL_EQUAL(double d1, double d2, double tolerance)
{
    return fabs(d1-d2) <= tolerance;
}

static void gaussian_grid(grib_handle* h)
{
    const double tolerance = 1.0/1000000.0; /* angular tolerance for grib2: micro degrees */
    long n = get(h,"numberOfParallelsBetweenAPoleAndTheEquator"); /* This is the key N */
    static double* values = NULL;
    static long last_n = 0;
    double north = dget(h,"latitudeOfFirstGridPointInDegrees");
    double south = dget(h,"latitudeOfLastGridPointInDegrees");

    double west = dget(h,"longitudeOfFirstGridPointInDegrees");
    double east = dget(h,"longitudeOfLastGridPointInDegrees");

    int e;

    if(n != last_n)
    {
        if(values) free(values);

        values = (double*)malloc(2*sizeof(double)*n);
        if(!values)
        {
            printf("%s, field %d [%s]: failed to allocate %ld bytes\n",file,field,param,2*(long)sizeof(double)*(n));
            error++;
            return;
        }

        if((e = grib_get_gaussian_latitudes(n,values)) != GRIB_SUCCESS)
        {
            printf("%s, field %d [%s]: cannot get gaussian latitudes for N%ld: %s\n",file,field,param,n,grib_get_error_message(e));
            error++;
            free(values);
            last_n = 0;
            return;
        }

        last_n = n;
    }

    if (!values) { assert(0); return; }
    if (values) {
        values[0] = rint(values[0]*1e6)/1e6;
    }

    if ( !DBL_EQUAL(north, values[0], tolerance) || !DBL_EQUAL(south, -values[0], tolerance) )
        printf("N=%ld north=%f south=%f v(=gauss_lat[0])=%f north-v=%0.30f south-v=%0.30f\n",
                n,north,south,values[0],north-values[0],south+values[0]);

    CHECK(DBL_EQUAL(north, values[0], tolerance));
    CHECK(DBL_EQUAL(south, -values[0], tolerance));

    if(missing(h,"numberOfPointsAlongAParallel"))   /* same as key Ni */
    {
        /* If missing, this is a REDUCED gaussian grid */
        const long MAXIMUM_RESOLUTION = 640;
        CHECK(get(h,"PLPresent"));
        CHECK(DBL_EQUAL(west, 0.0, tolerance));
        if (n > MAXIMUM_RESOLUTION) {
            printf("Gaussian number N (=%ld) cannot exceed %ld\n", n, MAXIMUM_RESOLUTION);
            CHECK(n <= MAXIMUM_RESOLUTION);
        }
    }
    else
    {
        /* REGULAR gaussian grid */
        long l_west = get(h,"longitudeOfFirstGridPoint");
        long l_east = get(h,"longitudeOfLastGridPoint");
        long parallel = get(h,"numberOfPointsAlongAParallel");
        long we = get(h,"iDirectionIncrement");
        double dwest  = dget(h,"longitudeOfFirstGridPointInDegrees");
        double deast  = dget(h,"longitudeOfLastGridPointInDegrees");
        double dwe = dget(h,"iDirectionIncrementInDegrees");
        /*printf("parallel=%ld east=%ld west=%ld we=%ld\n",parallel,east,west,we);*/

        CHECK(parallel == (l_east-l_west)/we + 1);
        CHECK(fabs((deast-dwest)/dwe + 1 - parallel) < 1e-10);
        CHECK(!get(h,"PLPresent"));
    }

    CHECK(ne(h,"Nj",0));

    if(get(h,"PLPresent")) {
        size_t count, i, nPl;
        int err_code = grib_get_size(h,"pl",&count);
        long *pl = NULL;
        double expected_lon2 = 0;
        long total, max_pl = 0;
        long numberOfValues = get(h,"numberOfValues");
        long numberOfDataPoints = get(h,"numberOfDataPoints");

        if(err_code) {
            printf("%s, field %d [%s]: cannot get size of pl: %s\n",file,field,param,grib_get_error_message(err_code));
            error++;
            return;
        }

        pl = (long*)malloc(sizeof(long)*count);
        CHECK(pl != NULL);
        if (!pl) return;

        nPl = count;
        if((err_code =  grib_get_long_array(h, "pl", pl, &count)))
        {
            printf("%s, field %d [%s]: cannot get pl: %s\n",file,field,param,grib_get_error_message(err_code));
            free(pl);
            error++;
            return;
        }
        if(nPl != count)
            printf("nPl=%ld count=%ld\n",(long)nPl,(long)count);

        CHECK(nPl == count);
        CHECK(nPl == (size_t)2*n);

        total = 0;
        max_pl = pl[0]; /* max elem of pl array = num points at equator */
        for(i = 0 ; i < count; i++) {
            total += pl[i];
            if (pl[i] > max_pl) max_pl = pl[i];
        }

        free(pl);

        /* Do not assume maximum of pl array is 4N! not true for octahedral */
        expected_lon2 = 360.0 - 360.0/max_pl;
        if (!DBL_EQUAL(expected_lon2, east, tolerance))
            printf("east actual=%g expected=%g diff=%g\n",east, expected_lon2, expected_lon2-east);
        CHECK(DBL_EQUAL(expected_lon2, east, tolerance));

        if(numberOfDataPoints != total)
            printf("GAUSS numberOfValues=%ld numberOfDataPoints=%ld sum(pl)=%ld\n",
                    numberOfValues,
                    numberOfDataPoints,
                    total);

        CHECK(numberOfDataPoints == total);

        CHECK(missing(h,"iDirectionIncrement"));
        CHECK(missing(h,"iDirectionIncrementInDegrees"));

        CHECK(eq(h,"iDirectionIncrementGiven",0));
        CHECK(eq(h,"jDirectionIncrementGiven",1));
    }

    CHECK(eq(h,"resolutionAndComponentFlags1",0));
    CHECK(eq(h,"resolutionAndComponentFlags2",0));
    CHECK(eq(h,"resolutionAndComponentFlags6",0));
    CHECK(eq(h,"resolutionAndComponentFlags7",0));
    CHECK(eq(h,"resolutionAndComponentFlags8",0));
}

static void check_validity_datetime(grib_handle* h)
{
    /* If we just set the stepRange (for non-instantaneous fields) to its
     * current value, then this causes the validity date and validity time
     * keys to be correctly computed.
     * Then we can compare the previous (possibly wrongly coded) value with
     * the newly computed one
     */
    char stepType[15]={0,};
    int err = 0;
    size_t str_len = 100;

    err = grib_get_string(h, "stepType", stepType, &str_len);
    if (err) return;
    if (strcmp(stepType, "instant")!=0) { /* not instantaneous */
        char stepRange[100]={0,};
        long saved_validityDate, saved_validityTime;
        long validityDate, validityTime;

        /* Check only applies to accumulated, max etc. */
        str_len = 100;
        err = grib_get_string(h, "stepRange", stepRange, &str_len);
        if (err) return;

        saved_validityDate = get(h, "validityDate");
        saved_validityTime = get(h, "validityTime");

        err = grib_set_string(h, "stepRange", stepRange, &str_len);
        if (err) return;
        validityDate = get(h, "validityDate");
        validityTime = get(h, "validityTime");

        if (validityDate!=saved_validityDate || validityTime!=saved_validityTime) {
            printf("warning: %s, field %d [%s]: invalid validity Date/Time (Should be %ld and %ld)\n",
                    file,field,param, validityDate, validityTime);
            warning++;
        }
    }
}

static void check_range(grib_handle* h,const parameter* p,double min,double max)
{
    double missing = 0;
    if(!valueflg)
        return;

    missing = dget(h,"missingValue");

    /* See ECC-437 */
    if (!(get(h,"bitMapIndicator") == 0 && min == missing && max == missing)) {
        if(min < p->min1 || min > p->min2)
        {
            printf("warning: %s, field %d [%s]: %s minimum value %g is not in [%g,%g]\n",file,field,param,
                    p->name,
                    min,p->min1,p->min2);
            printf("  => [%g,%g]\n",min < p->min1 ? min : p->min1, min > p->min2 ? min : p->min2);

            warning++;
        }

        if(max < p->max1 || max > p->max2 )
        {
            printf("warning: %s, field %d [%s]: %s maximum value %g is not in [%g,%g]\n",file,field,param,
                    p->name,
                    max,p->max1,p->max2);
            printf("  => [%g,%g]\n",max < p->max1 ? max : p->max1, max > p->max2 ? max : p->max2);
            warning++;
        }
    }
}

static void point_in_time(grib_handle* h,const parameter* p,double min,double max)
{
    switch(get(h,"typeOfProcessedData"))
    {
    case 0: /* Analysis */
        if (is_uerra)
            CHECK(eq(h,"productDefinitionTemplateNumber",0)||eq(h,"productDefinitionTemplateNumber",1));
        if (get(h,"productDefinitionTemplateNumber") == 1){
            CHECK(ne(h,"numberOfForecastsInEnsemble",0));
            CHECK(le(h,"perturbationNumber",get(h,"numberOfForecastsInEnsemble")));
        }
        break;

    case 1: /* Forecast */
        if (is_uerra)
            CHECK(eq(h,"productDefinitionTemplateNumber",0)||eq(h,"productDefinitionTemplateNumber",1));
        if (get(h,"productDefinitionTemplateNumber") == 1){
            CHECK(ne(h,"numberOfForecastsInEnsemble",0));
            CHECK(le(h,"perturbationNumber",get(h,"numberOfForecastsInEnsemble")));
        }
        break;

    case 2: /* Analysis and forecast products */
        CHECK(eq(h,"productDefinitionTemplateNumber",0));
        break;

    case 3: /* Control forecast products */
        CHECK(eq(h,"perturbationNumber",0));
        CHECK(ne(h,"numberOfForecastsInEnsemble",0));
        if (is_s2s_refcst)
            CHECK(eq(h,"productDefinitionTemplateNumber",60));
        else if (is_s2s)
            /*CHECK(eq(h,"productDefinitionTemplateNumber",60)||eq(h,"productDefinitionTemplateNumber",11)||eq(h,"productDefinitionTemplateNumber",1));*/
            CHECK(eq(h,"productDefinitionTemplateNumber",1));
        else
            CHECK(eq(h,"productDefinitionTemplateNumber",1));
        break;

    case 4: /* Perturbed forecast products */
        CHECK(ne(h,"perturbationNumber",0));
        CHECK(ne(h,"numberOfForecastsInEnsemble",0));
        if (is_s2s_refcst)
            CHECK(eq(h,"productDefinitionTemplateNumber",60));
        else if (is_s2s)
            /*CHECK(eq(h,"productDefinitionTemplateNumber",60)||eq(h,"productDefinitionTemplateNumber",11)||eq(h,"productDefinitionTemplateNumber",1));*/
            CHECK(eq(h,"productDefinitionTemplateNumber",1));
        else
            CHECK(eq(h,"productDefinitionTemplateNumber",1));
        if (is_lam) {
            CHECK(le(h,"perturbationNumber",get(h,"numberOfForecastsInEnsemble")));
        } else {
            /* Is there always cf in tigge global datasets?? */
            CHECK(le(h,"perturbationNumber",get(h,"numberOfForecastsInEnsemble")-1));
        }
        break;

    default:
        printf("Unsupported typeOfProcessedData %ld\n",get(h,"typeOfProcessedData"));
        CHECK(0);
        break;
    }

    if (is_lam) {
        if(get(h,"indicatorOfUnitOfTimeRange") == 10 ) /*  three hours */
        {
            /* Three hourly is OK */
            ;
        }
        else
        {
            CHECK(eq(h,"indicatorOfUnitOfTimeRange",1));/* Hours */
            CHECK((get(h,"forecastTime") % 3) == 0);  /* Every three hours */
        }
    }
    else if (is_uerra) {
        if(get(h,"indicatorOfUnitOfTimeRange") == 1) /*  hourly */
        {
            CHECK((eq(h,"forecastTime",1)||eq(h,"forecastTime",2)||eq(h,"forecastTime",4)||eq(h,"forecastTime",5))||(get(h,"forecastTime") % 3) == 0);
        }
    }
    else {
        if(get(h,"indicatorOfUnitOfTimeRange") == 11) /*  six hours */
        {
            /* Six hourly is OK */
            ;
        }
        else
        {
            CHECK(eq(h,"indicatorOfUnitOfTimeRange",1));/* Hours */
            CHECK((get(h,"forecastTime") % 6) == 0);  /* Every six hours */
        }
    }

    check_range(h,p,min,max);
}

static void height_level(grib_handle* h,const parameter* p,double min,double max)
{
    long level = get(h,"level");

    if (is_uerra){
        switch(level)
        {
        case   15:
        case   30:
        case   50:
        case   75:
        case  100:
        case  150:
        case  200:
        case  250:
        case  300:
        case  400:
        case  500:
            break;
        default:
            printf("%s, field %d [%s]: invalid height level %ld\n",file,field,param,level);
            error++;
            break;
        }
    }
}

static void pressure_level(grib_handle* h,const parameter* p,double min,double max)
{
    long level = get(h,"level");

    if (is_uerra && !is_crra){
        switch(level)
        {
        case 1000:
        case  975:
        case  950:
        case  925:
        case  900:
        case  875:
        case  850:
        case  825:
        case  800:
        case  750:
        case  700:
        case  600:
        case  500:
        case  400:
        case  300:
        case  250:
        case  200:
        case  150:
        case  100:
        case   70:
        case   50:
        case   30:
        case   20:
        case   10:
            break;
        default:
            printf("%s, field %d [%s]: invalid pressure level %ld\n",file,field,param,level);
            error++;
            break;
        }
    }
    else if (is_uerra && is_crra){
        switch(level)
        {
        case 1000:
        case  975:
        case  950:
        case  925:
        case  900:
        case  875:
        case  850:
        case  825:
        case  800:
        case  750:
        case  700:
        case  600:
        case  500:
        case  400:
        case  300:
        case  250:
        case  200:
        case  150:
        case  100:
        case   70:
        case   50:
        case   30:
        case   20:
        case   10:
        case    7:
        case    5:
        case    3:
        case    2:
        case    1:
            break;
        default:
            printf("%s, field %d [%s]: invalid pressure level %ld\n",file,field,param,level);
            error++;
            break;
        }
    }
    else if (is_s2s){
        switch(level)
        {
        case 1000:
        case  925:
        case  850:
        case  700:
        case  500:
        case  300:
        case  200:
        case  100:
        case  50:
        case  10:
            break;
        default:
            printf("%s, field %d [%s]: invalid pressure level %ld\n",file,field,param,level);
            error++;
            break;
        }
    } else {
        switch(level)
        {
        case 1000:
        case  200:
        case  250:
        case  300:
        case  500:
        case  700:
        case  850:
        case  925:
        case  50:
            break;
        default:
            printf("%s, field %d [%s]: invalid pressure level %ld\n",file,field,param,level);
            error++;
            break;
        }
    }
}

static void potential_vorticity_level(grib_handle* h,const parameter* p,double min,double max)
{
    long level = get(h,"level");

    switch(level)
    {
    case 2:
        break;
    default:
        printf("%s, field %d [%s]: invalid potential vorticity level %ld\n",file,field,param,level);
        error++;
        break;
    }
}

static void potential_temperature_level(grib_handle* h,const parameter* p,double min,double max)
{
    long level = get(h,"level");

    switch(level)
    {
    case 320:
        break;
    default:
        printf("%s, field %d [%s]: invalid potential temperature level %ld\n",file,field,param,level);
        error++;
        break;
    }
}

static void statistical_process(grib_handle* h,const parameter* p,double min,double max)
{
    switch(get(h,"typeOfProcessedData"))
    {
    case 0: /* Analysis */
        if (is_uerra)
            CHECK(eq(h,"productDefinitionTemplateNumber",8)||eq(h,"productDefinitionTemplateNumber",11));
        break;

    case 1: /* Forecast */
        if (is_uerra)
            CHECK(eq(h,"productDefinitionTemplateNumber",8)||eq(h,"productDefinitionTemplateNumber",11));
        break;

    case 2: /* Analysis and forecast products */
        CHECK(eq(h,"productDefinitionTemplateNumber",8));
        break;

    case 3: /* Control forecast products */
        if (!is_s2s_refcst)
            CHECK(eq(h,"productDefinitionTemplateNumber",11));
        else
            CHECK(eq(h,"productDefinitionTemplateNumber",61));
        break;

    case 4: /* Perturbed forecast products */
        if (!is_s2s_refcst)
            CHECK(eq(h,"productDefinitionTemplateNumber",11));
        else
            CHECK(eq(h,"productDefinitionTemplateNumber",61));
        break;

    default:
        printf("Unsupported typeOfProcessedData %ld\n",get(h,"typeOfProcessedData"));
        error++;
        return;
        break;
    }

    if (is_lam) {
        if(get(h,"indicatorOfUnitOfTimeRange") == 10 ) /*  three hours */
        {
            /* Three hourly is OK */
            ;
        }
        else
        {
            CHECK(eq(h,"indicatorOfUnitOfTimeRange",1));/* Hours */
            CHECK((get(h,"forecastTime") % 3) == 0);  /* Every three hours */
        }
    }
    else if (is_uerra)
    {
/*  forecastTime for uerra might be all steps decreased by 1 i.e 0,1,2,3,4,5,8,11...29 too many... */
        if(get(h,"indicatorOfUnitOfTimeRange") == 1)
        {
            CHECK(le(h,"forecastTime",30));
        }
    }
    else
    {
        if(get(h,"indicatorOfUnitOfTimeRange") == 11) /*  six hours */
        {
            /* Six hourly is OK */
            ;
        }
        else
        {
            CHECK(eq(h,"indicatorOfUnitOfTimeRange",1));/* Hours */
            CHECK((get(h,"forecastTime") % 6) == 0);  /* Every six hours */
        }
    }

    CHECK(eq(h,"numberOfTimeRange",1));
    CHECK(eq(h,"numberOfMissingInStatisticalProcess",0));
    CHECK(eq(h,"typeOfTimeIncrement",2));
    /*CHECK(eq(h,"indicatorOfUnitOfTimeForTheIncrementBetweenTheSuccessiveFieldsUsed",255));*/

    if (is_s2s)
        if(get(h,"typeOfStatisticalProcessing") == 0)
            CHECK(eq(h,"timeIncrementBetweenSuccessiveFields",1)||eq(h,"timeIncrementBetweenSuccessiveFields",4));
        else
            CHECK(eq(h,"timeIncrementBetweenSuccessiveFields",0));
    else
        CHECK(eq(h,"timeIncrementBetweenSuccessiveFields",0));


    CHECK(eq(h,"minuteOfEndOfOverallTimeInterval",0));
    CHECK(eq(h,"secondOfEndOfOverallTimeInterval",0));

    if (is_uerra)
    {
        CHECK((eq(h,"endStep",1)||eq(h,"endStep",2)||eq(h,"endStep",4)||eq(h,"endStep",5))||(get(h,"endStep") % 3) == 0);
    }
    else if (is_lam)
    {
        CHECK((get(h,"endStep") % 3) == 0);  /* Every three hours */
    }
    else
    {
        CHECK((get(h,"endStep") % 6) == 0); /* Every six hours */
    }

    if(get(h,"indicatorOfUnitForTimeRange") == 11)
    {
        /* Six hourly is OK */
        CHECK(get(h,"lengthOfTimeRange")*6 + get(h,"startStep") == get(h,"endStep"));
    }
    else if(get(h,"indicatorOfUnitForTimeRange") == 10)
    {
        /* Three hourly is OK */
        CHECK(get(h,"lengthOfTimeRange")*3 + get(h,"startStep") == get(h,"endStep"));
    }
    else
    {
        CHECK(eq(h,"indicatorOfUnitForTimeRange",1)); /* Hours */
        CHECK(get(h,"lengthOfTimeRange") + get(h,"startStep") == get(h,"endStep"));
    }
}

static void has_bitmap(grib_handle* h,const parameter* p,double min,double max)
{
    /* printf("bitMapIndicator %ld\n",get(h,"bitMapIndicator")); */
    CHECK(eq(h,"bitMapIndicator",0));
}

static void has_soil_level(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(get(h,"topLevel") == get(h,"bottomLevel"));
    CHECK(le(h,"level",14)); /* max in UERRA */
}

static void has_soil_layer(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(get(h,"topLevel") == get(h,"bottomLevel") - 1);
    CHECK(le(h,"level",14)); /* max in UERRA */
}

static void resolution_s2s(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(eq(h,"iDirectionIncrement",1500000));
    CHECK(eq(h,"jDirectionIncrement",1500000));
}

static void resolution_s2s_ocean(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(eq(h,"iDirectionIncrement",1000000));
    CHECK(eq(h,"jDirectionIncrement",1000000));
}

static void six_hourly(grib_handle* h,const parameter* p,double min,double max)
{
    statistical_process(h,p,min,max);

    if(get(h,"indicatorOfUnitForTimeRange") == 11)
        CHECK(eq(h,"lengthOfTimeRange",1));
    else
        CHECK(eq(h,"lengthOfTimeRange",6));

    CHECK(get(h,"endStep") == get(h,"startStep") + 6);
    check_range(h,p,min,max);
}

static void since_prev_pp(grib_handle* h,const parameter* p,double min,double max)
{
    statistical_process(h,p,min,max);

    CHECK(eq(h,"indicatorOfUnitForTimeRange",1));
    CHECK(get(h,"endStep") == get(h,"startStep") + get(h,"lengthOfTimeRange"));
    check_range(h,p,min,max);
}

static void three_hourly(grib_handle* h,const parameter* p,double min,double max)
{
    statistical_process(h,p,min,max);

    if(get(h,"indicatorOfUnitForTimeRange") == 11)
        CHECK(eq(h,"lengthOfTimeRange",1));
    else
        CHECK(eq(h,"lengthOfTimeRange",3));

    CHECK(get(h,"endStep") == get(h,"startStep") + 3);
    check_range(h,p,min,max);
}

static void from_start(grib_handle* h,const parameter* p,double min,double max)
{
    long step = get(h,"endStep");
    statistical_process(h,p,min,max);
    CHECK(eq(h,"startStep",0));

    if(step == 0){
        if(!is_uerra){
            CHECK(min == 0 && max == 0); /* ??? xxx */
        }
    }
    else
    {
        check_range(h,p,min/step,max/step);
    }
}

static void daily_average(grib_handle* h,const parameter* p,double min,double max)
{
    long step = get(h,"endStep");
    CHECK(get(h,"startStep") == get(h,"endStep") - 24);
    statistical_process(h,p,min,max);

    if(step == 0)
        CHECK(min == 0 && max == 0);
    else
    {
        check_range(h,p,min,max);
    }
}

static void given_level(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(ne(h,"typeOfFirstFixedSurface",255));
    CHECK(!missing(h,"scaleFactorOfFirstFixedSurface"));
    CHECK(!missing(h,"scaledValueOfFirstFixedSurface"));

    CHECK(eq(h,"typeOfSecondFixedSurface",255));
    CHECK(missing(h,"scaleFactorOfSecondFixedSurface"));
    CHECK(missing(h,"scaledValueOfSecondFixedSurface"));
}

static void predefined_level(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(ne(h,"typeOfFirstFixedSurface",255));
    CHECK(missing(h,"scaleFactorOfFirstFixedSurface"));
    CHECK(missing(h,"scaledValueOfFirstFixedSurface"));

    CHECK(eq(h,"typeOfSecondFixedSurface",255));
    CHECK(missing(h,"scaleFactorOfSecondFixedSurface"));
    CHECK(missing(h,"scaledValueOfSecondFixedSurface"));
}

static void predefined_thickness(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(ne(h,"typeOfFirstFixedSurface",255));
    CHECK(missing(h,"scaleFactorOfFirstFixedSurface"));
    CHECK(missing(h,"scaledValueOfFirstFixedSurface"));

    CHECK(ne(h,"typeOfSecondFixedSurface",255));
    CHECK(missing(h,"scaleFactorOfSecondFixedSurface"));
    CHECK(missing(h,"scaledValueOfSecondFixedSurface"));
}

static void given_thickness(grib_handle* h,const parameter* p,double min,double max)
{
    CHECK(ne(h,"typeOfFirstFixedSurface",255));
    CHECK(!missing(h,"scaleFactorOfFirstFixedSurface"));
    CHECK(!missing(h,"scaledValueOfFirstFixedSurface"));

    CHECK(ne(h,"typeOfSecondFixedSurface",255));
    CHECK(!missing(h,"scaleFactorOfSecondFixedSurface"));
    CHECK(!missing(h,"scaledValueOfSecondFixedSurface"));
}

static void latlon_grid(grib_handle* h)
{
    const double tolerance = 1.0/1000000.0; /* angular tolerance for grib2: micro degrees */
    long data_points = get(h,"numberOfDataPoints");
    long meridian    = get(h,"numberOfPointsAlongAMeridian");
    long parallel    = get(h,"numberOfPointsAlongAParallel");

    long north = get(h,"latitudeOfFirstGridPoint");
    long south = get(h,"latitudeOfLastGridPoint");
    long west  = get(h,"longitudeOfFirstGridPoint");
    long east  = get(h,"longitudeOfLastGridPoint");

    long ns = get(h,"jDirectionIncrement");
    long we = get(h,"iDirectionIncrement");

    double dnorth = dget(h,"latitudeOfFirstGridPointInDegrees");
    double dsouth = dget(h,"latitudeOfLastGridPointInDegrees");
    double dwest  = dget(h,"longitudeOfFirstGridPointInDegrees");
    double deast  = dget(h,"longitudeOfLastGridPointInDegrees");

    double dns = dget(h,"jDirectionIncrementInDegrees");
    double dwe = dget(h,"iDirectionIncrementInDegrees");

    if(eq(h,"basicAngleOfTheInitialProductionDomain",0))
    {
        CHECK(missing(h,"subdivisionsOfBasicAngle"));
    }
    else
    {
        /* long basic    = get(h,"basicAngleOfTheInitialProductionDomain"); */
        /* long division = get(h,"subdivisionsOfBasicAngle"); */
        CHECK(!missing(h,"subdivisionsOfBasicAngle"));
        CHECK(!eq(h,"subdivisionsOfBasicAngle",0));
    }

    if(missing(h,"subdivisionsOfBasicAngle"))
    {
        CHECK(eq(h,"basicAngleOfTheInitialProductionDomain",0));
    }

    CHECK(meridian*parallel == data_points);

    CHECK(eq(h,"resolutionAndComponentFlags1",0));
    CHECK(eq(h,"resolutionAndComponentFlags2",0));
    CHECK(eq(h,"resolutionAndComponentFlags6",0));
    CHECK(eq(h,"resolutionAndComponentFlags7",0));
    CHECK(eq(h,"resolutionAndComponentFlags8",0));

    CHECK(eq(h,"iDirectionIncrementGiven",1));
    CHECK(eq(h,"jDirectionIncrementGiven",1));

    CHECK(eq(h,"numberOfOctectsForNumberOfPoints",0));
    CHECK(eq(h,"interpretationOfNumberOfPoints",0));

    if(get(h,"iScansNegatively"))
    {
        long tmp    = east;
        double dtmp = deast;

        east = west;
        west = tmp;

        deast = dwest;
        dwest = dtmp;
    }

    if(get(h,"jScansPositively"))
    {
        long   tmp  = north;
        double dtmp = dnorth;

        north = south;
        south = tmp;

        dnorth = dsouth;
        dsouth = dtmp;
    }

    if (!(is_lam || is_uerra))
    {
        double area, globe;
        CHECK(north > south);
        CHECK(east  > west);

        /* Check that the grid is symmetrical */
        CHECK(north == -south);
        CHECK( DBL_EQUAL(dnorth, -dsouth, tolerance) );
        CHECK(parallel == (east-west)/we + 1);
        CHECK(fabs((deast-dwest)/dwe + 1 - parallel) < 1e-10);
        CHECK(meridian == (north-south)/ns + 1);
        CHECK(fabs((dnorth-dsouth)/dns + 1 - meridian) < 1e-10 );

        /* Check that the field is global */
        area  = (dnorth-dsouth) * (deast-dwest);
        globe = 360.0*180.0;
        CHECK(area <= globe);
        CHECK(area >= globe*0.95);
    }

    /* GRIB2 requires longitudes are always positive */
    CHECK(east >= 0);
    CHECK(west >= 0);

    /*
      printf("meridian=%ld north=%ld south=%ld ns=%ld \n",meridian,north,south,ns);
      printf("meridian=%ld north=%f south=%f ns=%f \n",meridian,dnorth,dsouth,dns);
      printf("parallel=%ld east=%ld west=%ld we=%ld \n",parallel,east,west,we);
      printf("parallel=%ld east=%f west=%f we=%f \n",parallel,deast,dwest,dwe);
     */

}

#define X(x) printf("%s=%ld ",#x,get(h,#x))

static void check_parameter(grib_handle* h,double min,double max)
{
    int err;
    int best = -1;
    int match = -1;
    size_t i = 0;

    for(i = 0; i < NUMBER(parameters); i++)
    {
        int j = 0;
        int matches = 0;
        while(parameters[i].pairs[j].key != NULL)
        {
            long val = -1;
            const int ktype = parameters[i].pairs[j].key_type;
            if (ktype == GRIB_TYPE_LONG) {
                if(grib_get_long(h,parameters[i].pairs[j].key,&val) == GRIB_SUCCESS) {
                    if(parameters[i].pairs[j].value_long == val) {
                        matches++;
                    }
                }
            }
            else if (ktype == GRIB_TYPE_STRING) {
                char strval[256]={0,};
                size_t len = 256;
                if (is_uerra && strcasecmp(parameters[i].pairs[j].key,"model")==0) {
                    /* printf("Skipping model keyword for UERRA class\n"); */
                    matches++; /*xxx hack to pretend that model key was matched.. */
                } else {
                    if (strcasecmp(parameters[i].pairs[j].value_string,"MISSING")==0) {
                        int is_miss = grib_is_missing(h, parameters[i].pairs[j].key, &err);
                        if (err == GRIB_SUCCESS && is_miss) {
                            matches++;
                        }
                    }
                    else if(grib_get_string(h,parameters[i].pairs[j].key,strval,&len) == GRIB_SUCCESS) {
                        if(strcmp(parameters[i].pairs[j].value_string, strval) == 0) {
                            matches++;
                        }
                    }
                }
            }
            else {
                assert(!"Unknown key type");
            }
            j++;
#if 0
            printf("%s %s %ld val -> %d %d %d\n",
                    parameters[i].pairs[j].key,
                    parameters[i].pairs[j].value_string,
                    val,
                    matches,
                    j,
                    best);
#endif
        }

        if(matches == j && matches > best)
        {
            best = matches;
            match = i;
        }
    }

    if(match >= 0)
    {

        /*int j = 0;*/
        param = parameters[match].name;
        i = 0;
        while(parameters[match].checks[i])
            (*parameters[match].checks[i++])(h,&parameters[match],min,max);
        /*
                printf("=========================\n");
                printf("%s -> %d %d\n",param, match, best);
                while(parameters[match].pairs[j].key != NULL)
                {
                     printf("%s val -> %ld %d\n",parameters[match].pairs[j].key,parameters[match].pairs[j].value,j);
                     j++;
                }
                printf("matched parameter: %s\n", param);
         */
    }
    else
    {
        printf("%s, field %d [%s]: cannot match parameter\n",file,field,param);
        X(origin);
        X(discipline);
        X(parameterCategory);
        X(parameterNumber);
        X(typeOfFirstFixedSurface);
        X(scaleFactorOfFirstFixedSurface);
        X(scaledValueOfFirstFixedSurface);
        X(typeOfSecondFixedSurface);
        X(scaleFactorOfSecondFixedSurface);
        X(scaledValueOfSecondFixedSurface);
        printf("\n");
        error++;
    }
}

static void check_packing(grib_handle* h)
{
    /* ECC-1009: Warn if not using simple packing */
    int err = 0;
    char packingType[254] = {0,};
    size_t len = sizeof(packingType);
    const char* expected_packingType = "grid_simple";

    err = grib_get_string(h, "packingType", packingType, &len);
    if (err) return;
    if (strcmp(packingType, expected_packingType)!=0) {
        printf("warning: %s, field %d [%s]: invalid packingType %s (Should be %s)\n",
               file, field, param, packingType, expected_packingType);
        warning++;
    }
}

static void verify(grib_handle* h)
{
    double min = 0,max = 0;

    CHECK(eq(h,"editionNumber",2));
    CHECK(missing(h,"reserved") || eq(h,"reserved",0));

    if (valueflg)
    {
        size_t count, n;
        int e = grib_get_size(h,"values",&count);
        double *values;
        size_t i;
        int bitmap = !eq(h,"bitMapIndicator",255);

        if(e) {
            printf("%s, field %d [%s]: cannot number of values: %s\n",file,field,param,grib_get_error_message(e));
            error++;
            return;
        }

        values = (double*)malloc(sizeof(double)*(count));
        if(!values)
        {
            printf("%s, field %d [%s]: failed to allocate %ld bytes\n",file,field,param,(long)(sizeof(double)*count));
            error++;
            return;
        }

        CHECK(eq(h,"numberOfDataPoints",count));

        n = count;
        if((e =  grib_get_double_array(h,"values",values,&count)))
        {
            printf("%s, field %d [%s]: cannot get values: %s\n",file,field,param,grib_get_error_message(e));
            free(values);
            error++;
            return;
        }

        if(n != count)
        {
            printf("%s, field %d [%s]: value count changed %ld -> %ld\n",file,field,param,(long)count,(long)n);
            free(values);
            error++;
            return;
        }

        if(bitmap)
        {
            double missing = dget(h,"missingValue");

            min = max = missing;

            for(i = 0; i < count ; i++)
            {
                if((min == missing) || ((values[i] != missing) && (min>values[i]))) min = values[i];
                if((max == missing) || ((values[i] != missing) && (max<values[i]))) max = values[i];
            }
        }
        else
        {
            min = max = values[0];
            for(i = 1; i < count ; i++)
            {
                if(min>values[i]) min = values[i];
                if(max<values[i]) max = values[i];
            }
        }
        free(values);
    }

    check_parameter(h,min,max);

    check_packing(h);

    /* Section 1 */

    CHECK(ge(h,"gribMasterTablesVersionNumber",4));
    CHECK(eq(h,"versionNumberOfGribLocalTables",0)); /* Local tables not used */

    CHECK(eq(h,"significanceOfReferenceTime",1)); /* Start of forecast */

    if (!is_s2s){
        /* todo check for how many years back the reforecast is done? Is it coded in the grib??? */
        /* Check if the date is OK */
        {
            long date = get(h,"date");
            /* CHECK(date > 20060101); */

            CHECK( (date / 10000)         == get(h,"year"));
            CHECK( ((date % 10000) / 100) == get(h,"month"));
            CHECK( ((date % 100))         == get(h,"day"));
        }
    }

    if (is_uerra){
        CHECK(le(h,"hour",24));
    }
    else if (is_lam){
        CHECK(eq(h,"hour",0) || eq(h,"hour",3) || eq(h,"hour",6) || eq(h,"hour",9) || eq(h,"hour",12) || eq(h,"hour",15) || eq(h,"hour",18) || eq(h,"hour",21));
    }
    else
    {
    /* Only 00, 06 12 and 18 Cycle OK */
        CHECK(eq(h,"hour",0) || eq(h,"hour",6) || eq(h,"hour",12) || eq(h,"hour",18));
    }
    CHECK(eq(h,"minute",0));
    CHECK(eq(h,"second",0));
    CHECK(ge(h,"startStep",0));

    if (is_s2s){
        CHECK(eq(h,"productionStatusOfProcessedData",6)||eq(h,"productionStatusOfProcessedData",7)); /*  S2S prod||test */
        CHECK(le(h,"endStep",100*24));
    }
    else if (!is_uerra)
    {
        CHECK(eq(h,"productionStatusOfProcessedData",4)||eq(h,"productionStatusOfProcessedData",5)); /*  TIGGE prod||test */
        CHECK(le(h,"endStep",30*24));
    }

    if (is_uerra){
        CHECK((eq(h,"step",1)||eq(h,"step",2)||eq(h,"step",4)||eq(h,"step",5))||(get(h,"step") % 3) == 0);
    }
    else if (is_lam) {
        CHECK((get(h,"step") % 3) == 0);
    }
    else
    {
        CHECK((get(h,"step") % 6) == 0);
    }

    if (is_uerra){
        if (is_crra){
            CHECK(eq(h,"productionStatusOfProcessedData",10)||eq(h,"productionStatusOfProcessedData",11)); /*  CRRA prod||test */
        }
        else
        {
            CHECK(eq(h,"productionStatusOfProcessedData",8)||eq(h,"productionStatusOfProcessedData",9)); /*  UERRA prod||test */
        }
        CHECK(le(h,"endStep",30));
        /* 0 = analysis , 1 = forecast */
        CHECK(eq(h,"typeOfProcessedData",0)||eq(h,"typeOfProcessedData",1));
        if (get(h,"typeOfProcessedData") == 0)
            CHECK(eq(h,"step",0));
        else
            CHECK((eq(h,"step",1)||eq(h,"step",2)||eq(h,"step",4)||eq(h,"step",5))||(get(h,"step") % 3) == 0);
    }
    else
    {
        /* 2 = analysis or forecast , 3 = control forecast, 4 = perturbed forecast */
        CHECK(eq(h,"typeOfProcessedData",2)||eq(h,"typeOfProcessedData",3)||eq(h,"typeOfProcessedData",4));
    }

    /* TODO: validate local usage. Empty for now xxx*/
    /* CHECK(eq(h,"section2.sectionLength",5)); */

    /* Section 3 */

    CHECK(eq(h,"sourceOfGridDefinition",0)); /* Specified in Code table 3.1 */

    switch(get(h,"gridDefinitionTemplateNumber"))
    {
    case 0:
    case 1: /*rotated latlon*/
        latlon_grid(h);
        break;

    case 30: /*Lambert conformal*/
        /*lambert_grid(h); # TODO xxx
        printf("warning: Lambert grid - geometry checking not implemented yet!\n"); */
       /*CHECK(eq(h,"scanningMode",64));*/ /* M-F data used to have it wrong.. but it might depends on other projection set up as well!*/
        break;

    case 40: /* gaussian grid (regular or reduced) */
        gaussian_grid(h);
        break;

    default:
        printf("%s, field %d [%s]: Unsupported gridDefinitionTemplateNumber %ld\n",
                file,field,param,
                get(h,"gridDefinitionTemplateNumber"));
        error++;
        return;
        break;
    }

    /* If there is no bitmap, this should be true */
    /* CHECK(eq(h,"bitMapIndicator",255));*/

    if(eq(h,"bitMapIndicator",255))
        CHECK(get(h,"numberOfValues") == get(h,"numberOfDataPoints"));
    else
        CHECK(get(h,"numberOfValues") <= get(h,"numberOfDataPoints"));

    /* Check values */
    CHECK(eq(h,"typeOfOriginalFieldValues",0)); /* Floating point */

    check_validity_datetime(h);

    /* do not store empty values e.g. fluxes at step 0
        todo ?? now it's allowed in the code here!
        if(!missing(h,"typeOfStatisticalProcessing"))
          CHECK(ne(h,"stepRange",0));*/
}

void validate(const char* path)
{
    FILE *f = fopen(path,"rb");
    grib_handle *h = 0;
    int err;
    int count = 0;

    file  = path;
    field = 0;

    if(!f) {
        printf("%s: %s\n",path,strerror(errno));
        error++;
        return;
    }

    while( (h= grib_handle_new_from_file(0,f,&err)) != NULL)
    {
        int last_error   = error;
        int last_warning = warning;

        ++field;
        verify(h);

        if((last_error != error) || (warnflg && (last_warning != warning)))
            save(h,bad,fbad);
        else
            save(h,good,fgood);


        grib_handle_delete(h);
        count++;
        param = "unknown";
    }
    fclose(f);

    if(err) {
        printf("%s: grib_handle_new_from_file: %s\n",path,grib_get_error_message(err));
        error++;
        return;
    }

    if(count == 0) {
        printf("%s does not contain any GRIBs\n",path);
        error++;
        return;
    }
}

static void usage()
{
    printf("tigge_check [options] grib_file grib_file ...\n");
    printf("   -l: check local area model fields\n");
    printf("   -v: check value ranges\n");
    printf("   -w: warnings are treated as errors\n");
    printf("   -g: write good gribs\n");
    printf("   -b: write bad  gribs\n");
    printf("   -z: return 0 to calling shell\n");
    printf("   -s: check s2s fields\n");
    printf("   -r: check s2s reforecast fields\n");
    printf("   -u: check uerra fields\n");
    printf("   -c: check crra fields (-u must be also used in this case)\n");
    exit(1);
}

int main(int argc, char** argv)
{
    int i;
    int err = 0;

    for(i = 1; i < argc; i++)
    {
        if(argv[i][0] == '-')
        {
            switch(argv[i][1])
            {
            case 'w':
                warnflg++;
                break;

            case 'z':
                zeroflg++;
                break;

            case 'v':
                valueflg++;
                break;

            case 'g':
                if(++i == argc) usage();
                good = argv[i];
                fgood = fopen(good,"w");
                if(!fgood) { perror(good); exit(1); }
                break;

            case 'b':
                if(++i == argc) usage();
                bad = argv[i];
                fbad = fopen(bad,"w");
                if(!fbad) { perror(bad); exit(1); }
                break;

            case 'l':
                is_lam=1;
                break;

            case 's':
                is_s2s=1;
                break;

            case 'r':
                is_s2s_refcst=1;
                break;

            case 'u':
                is_uerra=1;
                break;

            case 'c':
                is_crra=1;
                break;


            default:
                usage();
                break;
            }
        }
        else
        {
            break;
        }
    }

    if(i == argc)
        usage();

    for(; i < argc; i++)
    {
        error = 0;
        scan(argv[i]);
        if(error) err = 1;
        if(warning && warnflg) err = 1;
    }

    if(fgood && fclose(fgood)) { perror(good); exit(1); }
    if(fbad  && fclose(fbad )) { perror(bad ); exit(1); }

    return zeroflg ? 0: err;
}