eccodes/tools/grib_to_netcdf.cc

/*
 * (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.
 */

#include "grib_api_internal.h"

#ifdef HAVE_NETCDF

#include <cmath>
#include <stdarg.h>
#include <stdlib.h>
#include <ctype.h>

#include <time.h>
#include <float.h>
#include <netcdf.h>

#include "grib_tools.h"
#include "eccodes_windef.h"

#ifdef ECCODES_ON_WINDOWS
#include <stdint.h>
#endif

const char* tool_description =
    "Convert GRIB file(s) to netCDF format."
    "\n\tNote: The GRIB geometry should be a regular lat/lon grid or a regular Gaussian grid"
    "\n\t(the key \"typeOfGrid\" should be \"regular_ll\" or \"regular_gg\")";
const char* tool_name        = "grib_to_netcdf";
const char* tool_online_doc  = "https://confluence.ecmwf.int/display/ECC/grib_to_netcdf";
const char* tool_usage       = "[options] -o output_file grib_file grib_file ... ";
static char argvString[2048] = {0,};

/*=====================================================================*/

static grib_context* ctx           = NULL;
static double global_missing_value = 9.9692099683868690e+36; /* See GRIB-953 */

/*===============================================================================*/
/* request from mars client */
/*===============================================================================*/

#define NUMBER(x) (sizeof(x) / sizeof(x[0]))

typedef int err;
typedef off_t file_offset;

static int files = 0;

struct value
{
    struct value* next;
    char* name;
};
typedef struct value value;

struct request
{
    struct request* next;
    struct parameter* params;
    char* name;
    int order;
};
typedef struct request request;
/*  request part */
/* language part */

struct parameter
{
    struct parameter* next;
    struct value* values;
    char* name;
    int count;
};
typedef struct parameter parameter;

static const char* get_value(const request*, const char* name, int n);
static bool parsedate(const char* name, long* julian, long* second, bool* isjul);

void usage_and_exit(void)
{
    usage(); // this calls exit(1)
}

static bool eq_string(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static bool eq_integer(const char* l, const char* r)
{
    if (l && r)
        return atol(l) == atol(r);
    return false;
}

static bool eq_null(const char* l, const char* r)
{
    return true;
}

static bool eq_coord(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static bool eq_range(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static bool eq_date(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static bool eq_param(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static bool eq_time(const char* l, const char* r)
{
    if (l && r)
        return strcmp(l, r) == 0;
    return false;
}

static value* new_value(const char* name)
{
    value* v = (value*)calloc(1, sizeof(value));
    Assert(v);
    v->name = grib_context_strdup(ctx, name);
    return v;
}

static parameter* find_parameter(const request* r, const char* parname)
{
    if (!parname)
        return 0;
    if (r) {
        parameter* p = r->params;
        while (p) {
            if (strcmp(p->name, parname) == 0)
                return p;
            p = p->next;
        }
    }
    return NULL;
}

static void ecc_reqmerge(parameter* pa, const parameter* pb, request* a)
{
    const value* vb = pb->values;

    if (strcmp(pa->name, pb->name) != 0)
        return;

    while (vb) {
        value* va      = pa->values;
        value* last    = 0;
        const char* nb = vb->name;
        bool add    = true;

        while (va) {
            if (strcmp(va->name, nb) == 0) {
                add = false;
                break;
            }

            last = va;
            va   = va->next;
        }

        if (add) {
            value* v = new_value(nb);
            if (last)
                last->next = v;
            else
                pa->values = v;
            pa->count = 0;
        }

        vb = vb->next;
    }
}

/* Fast version if a && b same */
static bool ecc_reqmerge1(request* a, const request* b)
{
    parameter* pa       = a->params;
    const parameter* pb = b->params;

    while (pa && pb) {
        if (strcmp(pa->name, pb->name) != 0)
            return false;

        ecc_reqmerge(pa, pb, a);

        pa = pa->next;
        pb = pb->next;
    }

    return (pa == NULL && pb == NULL);
}

static void free_one_value(value* p)
{
    grib_context_free(ctx, p->name);
    grib_context_free(ctx, p);
}

static void free_all_values(value* p)
{
    while (p) {
        value* q = p->next;
        free_one_value(p);
        p = q;
    }
}

/* Convert the first part of the string 'p' to a number (x) and set n to its length. */
/* Return the rest of the string */
static const char* parse1(const char* p, int* x, int* n)
{
    *x = *n = 0;
    while (*p && isdigit(*p)) {
        (*x) *= 10;
        (*x) += *p - '0';
        (*n)++;
        p++;
    }
    return p;
}

static bool is_number(const char* name)
{
    const char* p = name;
    int x, n;

    if (p == 0 || *p == 0)
        return false;

    if (*p == '-')
        p++;
    else if (*p == '+')
        p++;

    p = parse1(p, &x, &n);
    if (n == 0 && *p != '.')
        return false;

    if (*p == '.') {
        p++;
        p = parse1(p, &x, &n);
    }

    if (*p == 'e' || *p == 'E') {
        p++;
        if (*p == '-')
            p++;
        else if (*p == '+')
            p++;
        p = parse1(p, &x, &n);
        if (n == 0)
            return false;
    }

    return *p == 0 ? true : false;
}
static parameter* new_parameter(char* name, value* v)
{
    parameter* p = (parameter*)calloc(1, sizeof(parameter));
    Assert(p);
    p->name   = grib_context_strdup(ctx, name);
    p->values = v;
    return p;
}

static void put_value(request* r, const char* parname, const char* valname, bool append, bool unique, bool ascending)
{
    parameter* p;
    value* v;

    if (!r)
        return;

    if ((p = find_parameter(r, parname)) != NULL) {
        if (append) {
            value *a = p->values, *b = NULL, *c = NULL;
            while (a) {
                b = a;
                if (unique) {
                    if (is_number(a->name) && is_number(valname)) {
                        if (atof(a->name) == atof(valname))
                            return;
                    }
                    else if (strcmp(a->name, valname) == 0)
                        return;
                }

                if (ascending) {
                    if (is_number(a->name)) {
                        if (atof(valname) < atof(a->name))
                            break;
                    }
                    else if (strcmp(valname, a->name) < 0)
                        break;
                }
                c = b;
                a = a->next;
            }
            v = new_value(grib_context_strdup(ctx, valname));
            if (ascending) {
                if (c) {
                    if (b && b != c)
                        v->next = b;
                    c->next = v;
                }
                else {
                    if (a)
                        v->next = a;
                    p->values = v;
                }
            }
            else {
                if (b)
                    b->next = v;
                else
                    p->values = v;
            }
            /* p->count++; */
            p->count = 0;
        }
        else {
            if (p->values) {
                free_all_values(p->values->next);
                p->values->next = NULL;
                /* p->count = 1; */
                p->count = 0;
                if (strcmp(p->values->name, valname) == 0)
                    return;
                else {
                    grib_context_free(ctx, p->values->name);
                    p->values->name = grib_context_strdup(ctx, valname);
                }
            }
            else {
                v         = new_value(grib_context_strdup(ctx, valname));
                p->values = v;
                /* p->count = 1; */
                p->count = 0;
            }
        }
    }
    else {
        parameter* q = NULL;
        parameter* s = r->params;
        v            = new_value(grib_context_strdup(ctx, valname));
        p            = new_parameter(grib_context_strdup(ctx, parname), v);
        while (s) {
            q = s;
            s = s->next;
        }
        if (q)
            q->next = p;
        else
            r->params = p;
    }
}

static void add_value(request* r, const char* parname, const char* fmt, ...)
{
    char buffer[1024];
    va_list list;

    va_start(list, fmt);
    vsnprintf(buffer, sizeof(buffer), fmt, list);
    va_end(list);

    put_value(r, parname, buffer, true, false, false);
}

static void ecc_reqmerge2(request* a, const request* b)
{
    const parameter* pb = b->params;

    while (pb) {
        parameter* pa = find_parameter(a, pb->name);

        if (pa == NULL) {
            value* v = pb->values;
            while (v) {
                put_value(a, pb->name, v->name, true, true, false);
                v = v->next;
            }
        }
        else {
            ecc_reqmerge(pa, pb, a);
        }

        pb = pb->next;
    }
}

static void reqmerge(request* a, const request* b)
{
    if (a && b) {
        if (!ecc_reqmerge1(a, b))
            ecc_reqmerge2(a, b);
    }
}

static void save_name(FILE* f, const char* name, int n)
{
    int i = 0, cnt = 0;

    if (name == NULL)
        name = "(null)";

    cnt = fprintf(f, "%s", name);

    for (i = cnt; i < n; i++)
        putc(' ', f);
}

static void save_one_value(FILE* f, value* r)
{
    save_name(f, r->name, 0);
}

static void save_all_values(FILE* f, value* r)
{
    while (r) {
        save_one_value(f, r);
        if (r->next)
            putc('/', f);
        r = r->next;
    }
}

static void save_all_parameters(FILE* f, parameter* r)
{
    while (r) {
        if (r->values) {
            fprintf(f, ",\n    ");
            save_name(f, r->name, 10);
            fprintf(f, " = ");
            save_all_values(f, r->values);
        }

        r = r->next;
    }
    putc('\n', f);
}

static void save_one_request(FILE* f, const request* r)
{
    if (r) {
        save_name(f, r->name, 0);
        save_all_parameters(f, r->params);
        putc('\n', f);
    }
}

static void save_all_requests(FILE* f, const request* r)
{
    while (r) {
        save_one_request(f, r);
        r = r->next;
    }
}

/* Not used for the moment
static void print_one_request(const request *r)
{
    save_one_request(stdout, r);
}
 */

static void print_all_requests(const request* r)
{
    save_all_requests(stdout, r);
}

static void free_one_parameter(parameter* p)
{
    grib_context_free(ctx, p->name);
    free_all_values(p->values);
    /*free_all_requests(p->interface);*/
    grib_context_free(ctx, p);
}

static void free_all_parameters(parameter* p)
{
    while (p) {
        parameter* q = p->next;
        free_one_parameter(p);
        p = q;
    }
}

static void free_one_request(request* r)
{
    grib_context_free(ctx, r->name);
    free_all_parameters(r->params);
    grib_context_free(ctx, r);
}

static void free_all_requests(request* p)
{
    while (p) {
        request* q = p->next;
        free_one_request(p);
        p = q;
    }
}

static void set_value(request* r, const char* parname, const char* fmt, ...)
{
    char buffer[10240];
    va_list list;

    va_start(list, fmt);
    vsnprintf(buffer, sizeof(buffer), fmt, list);
    va_end(list);

    put_value(r, parname, buffer, false, false, false);
}

static err handle_to_request(request* r, grib_handle* g)
{
    grib_keys_iterator* ks;
    char name[256];
    char value[256];
    size_t len = sizeof(value);
    int e      = 0;

    if (!g)
        return -1;

    /* printf("------------\n"); */
    ks = grib_keys_iterator_new(g, GRIB_KEYS_ITERATOR_ALL_KEYS, "mars");

    while (grib_keys_iterator_next(ks)) {
        strcpy(name, grib_keys_iterator_get_name(ks));

        if ((e = grib_keys_iterator_get_string(ks, value, &len)) != GRIB_SUCCESS)
            grib_context_log(ctx, GRIB_LOG_ERROR, "Cannot get %s as string %d (%s)", name, e, grib_get_error_message(e));

        set_value(r, name, "%s", value);
        len = sizeof(value);
    }

    strcpy(name, "stepUnits");
    if ((e = grib_get_string(g, name, value, &len)) == GRIB_SUCCESS) {
        set_value(r, name, "%s", value);
    }
    else {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Cannot get %s as string (%s)", name, grib_get_error_message(e));
    }

    /*
     Assert(grib_get_long(g, "validityDate", &l ) == 0);
     set_value(r, "validityDate", "%ld", l);

     Assert(grib_get_long(g, "validityTime", &l ) == 0);
     set_value(r, "validityTime", "%ld", l);
     */

    len = sizeof(value);
    if (grib_get_string(g, "cfVarName", name, &len) == 0) {
        if (strcmp(name, "unknown") != 0) {
            set_value(r, "param", "%s", name);
        }
        else {
            len = sizeof(value);
            if (grib_get_string(g, "shortName", name, &len) == 0) {
                set_value(r, "param", "%s", name);
            }
        }
    }

    len = sizeof(value);
    if (grib_get_string(g, "name", name, &len) == 0) {
        if (strcmp(name, "unknown") != 0) {
            set_value(r, "_long_name", "%s", name);
        }
    }

    len = sizeof(value);
    if (grib_get_string(g, "units", name, &len) == 0) {
        if (strcmp(name, "unknown") != 0) {
            set_value(r, "_units", "%s", name);
        }
    }

    len = sizeof(value);
    if (grib_get_string(g, "cfName", name, &len) == 0) {
        if (strcmp(name, "unknown") != 0) {
            set_value(r, "_cf_name", "%s", name);
        }
    }

    grib_keys_iterator_delete(ks);
    return e;
}

/*===============================================================================*/
/*===============================================================================*/

typedef bool (*namecmp)(const char*, const char*);

typedef struct hypercube
{
    request* cube;
    request* r;
    request* iterator;
    char* set;
    int count;
    int size;
    int max;
    int* index_cache;
    int index_cache_size;
    namecmp* compare;
} hypercube;

typedef struct axis_t
{
    const char* name;
    namecmp compare;
} axis_t;

/* This should be c++ ... */

typedef enum field_state
{
    unknown,
    packed_mem,
    packed_file,
    expand_mem
} field_state;

typedef struct
{
    int refcnt;
    request* r;
} field_request;

/* One field .. */

typedef struct field
{
    int refcnt;

    field_state shape;

    grib_handle* handle;
    double* values;
    size_t value_count;

    /* if on file */

    file_offset offset;
    size_t length;
    grib_file* file;

    /* missing fields/values */

    /*bool is_missing;*/ /* field is missing */
    bool has_bitmap;     /* field has missing values (= bitmap) */

    field_request* r;

} field;

typedef struct fieldset
{
    int refcnt;

    /* if fields */

    int max;
    int count;

    field** fields;

} fieldset;

/*
 #define MISSING_VALUE(n)      ((n) == missing_value)
 #define MISSING_FIELD(f)      ((f)->missing)
 #define FIELD_HAS_BITMAP(f)   ((f)->bitmap)

 #define FASTNEW(type)         (type*)calloc(1, sizeof(type))
 #define grib_context_free(ctx,x)           grib_context_free(ctx,x)
 */

static field* get_field(fieldset* v, int n, field_state shape);
static hypercube* new_hypercube_from_mars_request(const request* r);
static void release_field(field* g);
static int count_axis(const hypercube* h);
static const char* get_axis(const hypercube* h, int pos);
static const char* get_axis(const hypercube* h, int pos);
static int cube_order(const hypercube* h, const request* r);
static void free_hypercube(hypercube* h);
static int ecc_cube_position(const hypercube* h, const request* r, bool remove_holes);

static value* clone_one_value(const value* p)
{
    value* q = (value*)calloc(1, sizeof(value));
    Assert(q);
    q->next = NULL;
    q->name = grib_context_strdup(ctx, p->name);
    return q;
}

static value* clone_all_values(const value* p)
{
    if (p) {
        value* q = clone_one_value(p);
        q->next  = clone_all_values(p->next);
        /* q->alias = cone_value(p->alias); */
        return q;
    }
    return NULL;
}

static parameter* clone_one_parameter(const parameter* p)
{
    parameter* q = (parameter*)calloc(1, sizeof(parameter));
    Assert(q);
    q->next   = NULL;
    q->name   = grib_context_strdup(ctx, p->name);
    q->values = clone_all_values(p->values);
    return q;
}

static parameter* clone_all_parameters(const parameter* p)
{
    if (p) {
        parameter* q = clone_one_parameter(p);
        q->next      = clone_all_parameters(p->next);
        return q;
    }
    return NULL;
}

static request* clone_one_request(const request* r)
{
    if (r) {
        request* p = (request*)calloc(1, sizeof(request));
        Assert(p);
        p->name = grib_context_strdup(ctx, r->name);

        p->params = clone_all_parameters(r->params);
        p->next   = NULL;
        return p;
    }
    return NULL;
}

static request* new_request(const char* name, parameter* p)
{
    request* r = (request*)calloc(1, sizeof(request));
    Assert(r);
    r->name   = grib_context_strdup(ctx, name);
    r->params = p;
    return r;
}

static request* empty_request(const char* name)
{
    return new_request(name ? name : "", NULL);
}

static field_request* new_field_request(request* r)
{
    field_request* g = (field_request*)grib_context_malloc_clear(ctx, sizeof(field_request));
    g->r             = clone_one_request(r);
    return g;
}

static void free_field_request(field_request* g)
{
    if (!g)
        return;
    g->refcnt--;
    if (g->refcnt <= 0) {
        free_all_requests(g->r);
        grib_context_free(ctx, g);
    }
}

static void free_field(field* g)
{
    if (!g)
        return;
    g->refcnt--;
    if (g->refcnt <= 0) {
        /*free_gribfile(g->file);*/
        free_field_request(g->r);
        if (g->values)
            grib_context_free(ctx, g->values);
        grib_handle_delete(g->handle);
        grib_context_free(ctx, g);
    }
}

static void free_fieldset(fieldset* v)
{
    if (!v)
        return;
    v->refcnt--;
    if (v->refcnt <= 0) {
        int i;
        grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: free_fieldset (%d fields) : ", v->count);

        for (i = 0; i < v->count; i++)
            free_field(v->fields[i]);

        grib_context_free(ctx, v->fields);
        grib_context_free(ctx, v);
    }
}

static field* new_field()
{
    return (field*)grib_context_malloc_clear(ctx, sizeof(field));
}

#define INIT_SIZE 1024
static void grow_fieldset(fieldset* v, int n)
{
    int m = v->count;
    int x = v->max;

    if (n < v->count)
        return;

    v->count = n;

    while (v->count >= v->max)
        if (v->max < INIT_SIZE)
            v->max = INIT_SIZE;
        else
            v->max += v->max / 2 + 1;

    if (v->max != x) {
        int i;
        if (v->fields == NULL) {
            v->fields = (field**)grib_context_malloc(ctx, sizeof(field*) * v->max);
            Assert(v->fields);
        }
        else {
            field** f = (field**)grib_context_malloc(ctx, sizeof(field*) * v->max);
            Assert(f);
            for (i = 0; i < m; i++)
                f[i] = v->fields[i];
            grib_context_free(ctx, v->fields);
            v->fields = f;
        }

        for (i = m; i < v->max; i++)
            v->fields[i] = NULL;
    }
}

static fieldset* new_fieldset(int n)
{
    fieldset* f = (fieldset*)calloc(1, sizeof(fieldset));
    Assert(f);
    grow_fieldset(f, n);
    return f;
}

static field* read_field(grib_file* file, file_offset pos, long length)
{
    field* g = new_field();

    g->file   = file;
    g->offset = pos;
    g->length = length;
    g->shape  = packed_file;

    return g;
}

static err set_field(fieldset* v, field* g, int pos)
{
    if (pos >= 0) {
        field* h;
        grow_fieldset(v, pos + 1);
        h              = v->fields[pos];
        v->fields[pos] = g;
        g->refcnt++;
        if (h)
            free_field(h);
    }
    return 0;
}

static void count_parval(parameter* p)
{
    int n    = 0;
    value* v = p->values;

    while (v) {
        n++;
        v = v->next;
    }

    p->count = n;
}

static int count_values(const request* r, const char* parname)
{
    parameter* p = find_parameter(r, parname);

    if (p == NULL)
        return 0;
    if (p->count)
        return p->count;

    count_parval(p);

    return p->count;
}
static const char* get_value(const request* r, const char* parname, int nth)
{
    parameter* p = find_parameter(r, parname);
    value* v;
    int i = 0;

    if (p == NULL)
        return NULL;

    if (!p->count)
        count_values(r, parname);

    v = p->values;

    while (v) {
        if (nth == i++)
            return v->name;
        v = v->next;
    }

    return NULL;
}

static err to_packed_mem(field* g)
{
    if (g->shape == packed_mem)
        return 0;

    if (g->shape == expand_mem) {
        if (g->values)
            grib_context_free(ctx, g->values);
        g->values      = NULL;
        g->value_count = 0;
        g->shape       = packed_mem;

        return 0;
    }

    if (g->shape == packed_file) {
    }
    return 0;
}

static err to_expand_mem(field* g)
{
    err e = 0;

    Assert(g);
    if (g->shape == expand_mem)
        return 0;

    if (g->shape == packed_file) {
        const void* dummy = NULL;

        grib_file* file = grib_file_open(g->file->name, "r", &e);
        if (!file || !file->handle) {
            grib_context_log(ctx, GRIB_LOG_ERROR | GRIB_LOG_PERROR, "%s", g->file->name);
            return -1;
        }

        fseeko(file->handle, g->offset, SEEK_SET);

        g->handle = grib_handle_new_from_file(ctx, file->handle, &e);
        Assert(g->handle);

        if (g->handle)
            grib_get_message(g->handle, &dummy, &g->length);

        grib_file_close(file->name, 0, &e);
        if (!g->handle)
            return -1;

        if (g->values)
            grib_context_free(ctx, g->values);
        g->values = NULL;
    }

    if (g->values == NULL) {
        size_t count = 0;
        long bitmap  = 0;

        if ((e = grib_get_size(g->handle, "values", &g->value_count))) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get number of values: %s", grib_get_error_message(e));
            return e;
        }

        count = g->value_count;

        if ((e = grib_set_double(g->handle, "missingValue", global_missing_value))) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot set missingValue: %s", grib_get_error_message(e));
            return e;
        }

        g->values = (double*)grib_context_malloc(ctx, sizeof(double) * g->value_count);
        if ((e = grib_get_double_array(g->handle, "values", g->values, &count))) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot decode values: %s", grib_get_error_message(e));
            return e;
        }

        if (count != g->value_count)
            grib_context_log(ctx, GRIB_LOG_FATAL, "ecCodes: value count mismatch %ld %ld", count, g->value_count);

        if ((e = grib_get_long(g->handle, "missingValuesPresent", &bitmap))) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get missingValuesPresent: %s", grib_get_error_message(e));
            return e;
        }

        g->has_bitmap = (bitmap != 0);

#ifdef COMEBACK
        set g->missing
#endif
    }

    g->shape = expand_mem;

    return e;
}

static void set_field_state(field* g, field_state shape)
{
    switch (shape) {
        case expand_mem:
            to_expand_mem(g);
            break;

        case packed_mem:
            to_packed_mem(g);
            break;

        case packed_file:
            release_field(g);
            break;

        default:
            grib_context_log(ctx, GRIB_LOG_FATAL, "Internal error %s %d", __FILE__, __LINE__);
            break;
    }
}

static field* get_field(fieldset* v, int n, field_state shape)
{
    field* g = v->fields[n];
    set_field_state(g, shape);
    return g;
}

static void release_field(field* g)
{
    if (g->file) {
        if (g->values)
            grib_context_free(ctx, g->values);
        g->values = NULL;
        g->shape  = packed_file;
        grib_handle_delete(g->handle);
        g->handle = NULL;
    }
}

static request* field_to_request(field* f)
{
    if (f->r == 0) {
        field_state state = f->shape;
        request* r        = empty_request(
#ifdef COMEBACK
            ((f->ksec1 == NULL) || (f->ksec1[2] != mars.computeflg)) ? "GRIB" : "COMPUTED");
#else
            "GRIB");
#endif

        set_field_state(f, packed_mem);
        handle_to_request(r, f->handle);
        set_field_state(f, state);

        f->r = new_field_request(r);
        free_all_requests(r);
    }
    return f->r->r;
}

static request* fieldset_to_request(fieldset* fs)
{
    int i;
    request* r = empty_request("GRIB");

    if (!fs) {
        free_one_request(r);
        return 0;
    }

    for (i = 0; i < fs->count; i++) {
        request* s = field_to_request(fs->fields[i]);
        reqmerge(r, s);
    }

    return r;
}

/*===============================================================================*/
/* hypercube from mars client */
/*===============================================================================*/

static bool eq_string(const char*, const char*);
static bool eq_integer(const char*, const char*);
static bool eq_range(const char*, const char*);
static bool eq_param(const char*, const char*);
static bool eq_coord(const char*, const char*);
static bool eq_date(const char*, const char*);
static bool eq_time(const char*, const char*);
static bool eq_null(const char*, const char*);

static axis_t global_axis[] = {

    /* From dhsbase.c 'check_grib' */
    {
        "class",
        eq_string,
    },
    {
        "type",
        eq_string,
    },
    {
        "stream",
        eq_string,
    },
    {
        "levtype",
        eq_string,
    },

    {
        "origin",
        eq_string,
    },
    {
        "product",
        eq_string,
    },
    {
        "section",
        eq_string,
    },
    {
        "method",
        eq_integer,
    },
    {
        "system",
        eq_integer,
    },

    /* testing */
    /* {"repres",      eq_null, }, */

    /* from field order */
    {
        "date",
        eq_date,
    },
    {
        "refdate",
        eq_date,
    },
    {
        "hdate",
        eq_date,
    },
    {
        "time",
        eq_time,
    },

    {
        "reference",
        eq_range,
    },
    {
        "step",
        eq_range,
    },
    {
        "fcmonth",
        eq_integer,
    },
    {
        "fcperiod",
        eq_range,
    },

    {
        "leadtime",
        eq_range,
    },
    {
        "opttime",
        eq_range,
    },

    {
        "expver",
        eq_string,
    },
    {
        "domain",
        eq_string,
    },

    {
        "diagnostic",
        eq_integer,
    },
    {
        "iteration",
        eq_integer,
    },

    {
        "quantile",
        eq_range,
    },
    {
        "number",
        eq_integer,
    },

    {
        "levelist",
        eq_coord,
    },
    {
        "latitude",
        eq_coord,
    },
    {
        "longitude",
        eq_coord,
    },
    {
        "range",
        eq_range,
    },

    {
        "param",
        eq_param,
    },

    {
        "ident",
        eq_integer,
    },
    {
        "obstype",
        eq_integer,
    },
    {
        "instrument",
        eq_integer,
    },

    {
        "frequency",
        eq_integer,
    },
    {
        "direction",
        eq_integer,
    },

    {
        "channel",
        eq_integer,
    },
};

static int axisindex(const char* name)
{
    size_t i = 0;
    for (i = 0; i < NUMBER(global_axis); i++) {
        if (STR_EQUAL(name, global_axis[i].name)){
            return (int)i;
        }
    }
    return -1;
}

static namecmp comparator(const char* name)
{
    static char* dontcompare = NULL;
    static bool first     = 1;
    int i                    = 0;

    if (first) {
        dontcompare = getenv("MARS_DONT_CHECK");
        first       = false;
    }

    if (dontcompare != NULL) {
        if (STR_EQUAL(dontcompare, name))
            return eq_null;
    }

    if ((i = axisindex(name)) != -1)
        return global_axis[i].compare;
    grib_context_log(ctx, GRIB_LOG_ERROR, "No comparator for %s", name);
    return eq_string;
}

/********************/
/* index accessors  */
/********************/

static int count_index(const hypercube* h)
{
    int i = 0, n = 0;

    for (i = 0; i < h->size; ++i)
        n += h->set[i];
    return n;
}

static int count_holes(const hypercube* h, int cnt)
{
    int i = 0, n = 0;

    for (i = 0; i < cnt; ++i)
        n += h->set[i];

    return (n == cnt) ? 0 : (cnt - n);
}

static void reset_index(hypercube* h, int v)
{
    memset(h->set, v, h->size);
}

static void set_index(hypercube* h, int index, int value)
{
    if (index < 0 || index >= h->count) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error, bad hypercube index %d", index);
        exit(1);
    }

    if (index >= h->max) {
        int old = h->max;
        while (index >= h->max)
            h->max += 4096;

        h->set = h->set ? (char*)grib_context_realloc(ctx, h->set, h->max) : (char*)grib_context_malloc(ctx, h->max);
        Assert(h->set);
        memset(h->set + old, 0, h->max - old);
    }

    if (index >= h->size)
        h->size = index + 1;

    h->set[index] = value;
}

/**************************/
/* End of index accessors */
/**************************/

/*******************/
/* axis accessors  */
/*******************/

static int count_axis(const hypercube* h)
{
    if (h && h->cube)
        return count_values(h->cube, "axis");

    return -1;
}

static const char* get_axis(const hypercube* h, int pos)
{
    const char* axis = NULL;
    if (pos < count_axis(h)) {
        axis = get_value(h->cube, "axis", pos);
    }
    return axis;
}

static void add_axis(hypercube* h, const char* axis)
{
    add_value(h->cube, "axis", "%s", axis);
}

static void unset_value(request* r, const char* parname)
{
    parameter *p, *q = NULL;
    if (!r)
        return;

    p = r->params;
    while (p) {
        if (STR_EQUAL(parname, p->name)) {
            if (q)
                q->next = p->next;
            else
                r->params = p->next;
            free_one_parameter(p);
            return;
        }
        q = p;
        p = p->next;
    }
}

static void reset_axis(hypercube* h)
{
    unset_value(h->cube, "axis");
}

static void valcpy(request* a, request* b, char* aname, char* bname)
{
    parameter* p;
    if (a && b && (p = find_parameter(b, bname))) {
        bool z = false;
        value* v  = p->values;
        while (v) {
            put_value(a, aname, v->name, z, false, false);
            z = true;
            v = v->next;
        }
    }
}

static void cube_values(hypercube* h, const char* p)
{
    valcpy(h->cube, h->r, (char*)p, (char*)p);
}

static int count_dimensions(const hypercube*, const char*);

static int set_axis(hypercube* h)
{
    int i     = 0;
    int count = (h && h->r) ? 1 : -1;

    reset_axis(h);
    for (i = (NUMBER(global_axis) - 1); i >= 0; --i) {
        int n = count_dimensions(h, global_axis[i].name);
        if (n > 1) {
            add_axis(h, global_axis[i].name);
            cube_values(h, global_axis[i].name);
            count *= n;
        }
    }

    return count;
}

/*************************/
/* End of axis accessors */
/*************************/

/*******************/
/* Cube dimensions */
/*******************/

static int count_dimensions(const hypercube* h, const char* axis)
{
    int dims = -1;
    if (h && h->r)
        dims = count_values(h->r, axis);
    return dims;
}

/**************************/
/* End of cube dimensions */
/**************************/

/**************************/
/* Auxiliary functions    */
/**************************/

static int count_hypercube(const request* r)
{
    int count = 1;
    size_t i = 0;
    for (i = 0; i < NUMBER(global_axis); ++i) {
        int c = count_values(r, global_axis[i].name);
        count *= c ? c : 1;
    }

    return count;
}

static int cube_order(const hypercube* h, const request* r)
{
    return ecc_cube_position(h, r, true);
}

static int cube_position(const hypercube* h, const request* r)
{
    return ecc_cube_position(h, r, false);
}

static void reserve_index_cache(hypercube* h, int size)
{
    if (size <= 0)
        return;

    if (h->index_cache != 0)
        grib_context_free(ctx, h->index_cache);
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: Allocating hypercube index_cache: %d entries", size);
    h->index_cache = (int*)calloc(size, sizeof(int));
    Assert(h->index_cache);
    h->index_cache_size = size;
}

static int ecc_cube_position(const hypercube* h, const request* r, bool remove_holes)
{
    request* cube = h->cube;
    int c         = count_axis(h);
    int index     = 0;
    int i         = 0;
    int n         = 1;
    int ok        = 0;

    if (h->index_cache == 0 || h->index_cache_size != c)
        reserve_index_cache((hypercube*)h, c);

    for (i = 0; i < c; ++i) {
        const char* axis = get_axis(h, i);
        const char* v    = get_value(r, axis, 0);
        const char* w    = NULL;
        int dims         = count_dimensions(h, axis);
        int k            = 0;
        int count        = count_values(cube, axis);
        int last         = h->index_cache[i];

        for (k = 0; k < count; k++) {
            int j = (k + last) % count;
            w     = get_value(cube, axis, j);
            if (h->compare ? h->compare[i](w, v) : (strcmp(w, v) == 0)) {
                index += j * n;
                n *= dims;
                ok++;
                ((hypercube*)h)->index_cache[i] = j;
                break;
            }
            else
                grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: ecc_cube_position, %s, %s != %s [%scompare function available]", axis, w, v, h->compare ? "" : "no ");
        }
    }

    if (remove_holes) {
        int holes = 0;
        if (count_index(h) != h->size)
            holes = count_holes(h, index);
        index -= holes;
    }

    return (ok == c) ? index : -1;
}

static void cube_indexes(
    const hypercube* h, request* r,
    char** times_array, size_t times_array_size,
    int* indexes, int size)
{
    request* cube = h->cube;
    int c         = count_axis(h);
    int i         = 0;
    int index     = 0;
    int n         = 1;

    if (size < c) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error in cube_indexes. size=%d < axis=%d", size, c);
    }

    if (h->index_cache == 0 || h->index_cache_size != c)
        reserve_index_cache((hypercube*)h, c);

    for (i = 0; i < c; ++i) {
        const char* axis           = get_axis(h, i);
        const char* v              = get_value(r, axis, 0);
        const char* w              = NULL;
        int dims                   = count_dimensions(h, axis);
        int j                      = 0;
        int k                      = 0;
        int count                  = count_values(cube, axis);
        int last                   = h->index_cache[i];
        const bool is_time_axis = (STR_EQUAL(axis, "time"));
        if (is_time_axis) {
            Assert(times_array);
            Assert(times_array_size == count);
        }

        for (k = 0; k < count; k++) {
            j = (k + last) % count;
            if (is_time_axis) {
                /* GRIB-792: use fast lookup */
                Assert(j >= 0 && j < times_array_size);
                w = times_array[j];
                /* For testing:
                 * Assert( strcmp(w, get_value(cube, axis, j))==0 );
                 * */
            }
            else {
                /* slow access method */
                w = get_value(cube, axis, j);
            }

            if (h->compare ? h->compare[i](w, v) : (w == v)) {
                index += j * n;
                n *= dims;
                ((hypercube*)h)->index_cache[i] = j;
                break;
            }
        }
        indexes[i] = j;
    }
    (void)index;
}

/*********************************/
/* End of Auxiliary functions    */
/*********************************/

static hypercube* new_hypercube(const request* r)
{
    hypercube* h = (hypercube*)calloc(1, sizeof(hypercube));
    int total = 0, count = 0;
    size_t n        = 0;
    const char* val = 0;
    Assert(h);
    h->r    = clone_one_request(r);
    h->cube = empty_request("CUBE");

    h->count = total = count_hypercube(r);
    count            = set_axis(h);

    h->compare = 0;

    if ((total != count) || (count == 0)) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error while computing hypercube fields");
        grib_context_log(ctx, GRIB_LOG_ERROR, "Number of fields in request %d", total);
        grib_context_log(ctx, GRIB_LOG_ERROR, "Number of fields in hypercube %d", count);
    }
    set_index(h, count - 1, 1);
    memset(h->set, 1, count);

    /* This is expensive, but makes the iterator with only
       those parameters found as axis */
    h->iterator = empty_request(0);
    for (n = 0; n < NUMBER(global_axis); ++n)
        if ((val = get_value(h->r, global_axis[n].name, 0)) != NULL)
            set_value(h->iterator, global_axis[n].name, val);

    return h;
}

static void print_hypercube(const hypercube* h)
{
    Assert(h);
    if (!h) return;
    print_all_requests(h->r);
    print_all_requests(h->cube);
    grib_context_log(ctx, GRIB_LOG_INFO, "%d active out of %d fields described\n", count_index(h), h->size);
}

static void free_hypercube(hypercube* h)
{
    free_all_requests(h->r);
    free_all_requests(h->cube);
    free_all_requests(h->iterator);
    grib_context_free(ctx, h->index_cache);
    grib_context_free(ctx, h->compare);
    grib_context_free(ctx, h->set);
    grib_context_free(ctx, h);
}

struct stuff_1
{
    hypercube* c;
    request* r;
};

static void reqcb_1(const request* r, int count, axis_t* names, char* vals[], void* data)
{
    struct stuff_1* s = (struct stuff_1*)data;
    int i;

    for (i = 0; i < count; i++)
        if (vals[i])
            set_value(s->r, names[i].name, vals[i]);

    set_index(s->c, cube_position(s->c, s->r), 1);
}

typedef void (*loopproc)(const request*, int, axis_t*, char**, void*);

static void names_loop(const request* r, loopproc proc, void* data);

static hypercube* new_hypercube_from_mars_request(const request* r)
{
    int i;
    int n;

    struct stuff_1 s;

//     const request *lang = mars_language_from_request(r);
//     int count = 0;
//     count = init_axis(lang);
//     grib_context_log(ctx,GRIB_LOG_DEBUG,"cube %s",r->kind);
//     /* print_all_requests(mars_language_from_request(r)); */
//     grib_context_log(ctx,GRIB_LOG_INFO,"NUMBER(axis): %d, number axisnew: %d",NUMBER(axis),count);

    s.c = new_hypercube(r);
    s.r = clone_one_request(r);

    reset_index(s.c, 0);
    names_loop(r, reqcb_1, &s);

    free_one_request(s.r);

    /* add single parameters */

    for (i = 0; i < NUMBER(global_axis); i++) {
        int m = count_values(r, global_axis[i].name);
        if (m == 1) {
            add_value(s.c->cube, "axis", global_axis[i].name);
            set_value(s.c->cube, global_axis[i].name, get_value(r, global_axis[i].name, 0));
        }
    }

    n = count_values(s.c->cube, "axis");
    if (n) {
        s.c->compare = (namecmp*)calloc(n, sizeof(namecmp));
        Assert(s.c->compare);
    }

    for (i = 0; i < n; i++)
        s.c->compare[i] = comparator(get_value(s.c->cube, "axis", i));

    return s.c;
}

/* This one doesn't have single parameters in CUBE */
static hypercube* new_simple_hypercube_from_mars_request(const request* r)
{
    int i;
    int n;

    struct stuff_1 s;
    s.c = new_hypercube(r);
    s.r = clone_one_request(r);

    reset_index(s.c, 0);

    names_loop(r, reqcb_1, &s);

    free_one_request(s.r);
    n = count_values(s.c->cube, "axis");
    if (n) {
        s.c->compare = (namecmp*)calloc(n, sizeof(namecmp));
        Assert(s.c->compare);
    }

    for (i = 0; i < n; i++)
        s.c->compare[i] = comparator(get_value(s.c->cube, "axis", i));

    return s.c;
}

/*===========================================================================================*/

/*===========================================================================================*/

/* TODO:
 - Print usage in log file
 - consider FCMONTH and Climatology
 - Build logic to create validationtime when only one of DATE or TIME or STEP have multiple values:
 for example: date=-1, time=12, step=24/48
 - parametrise the type of data for each axis (function fill_netcdf_dimensions & define_netcdf_dimensions)
 Now, all are INT but TIME could be float
 - allow user specified scale_factor
 - insert 'user input request'
 - ORIGIN
 - Monthly means seem not to ignore STEP (data server, era40, eg retrieve 39/142)
 */

typedef struct ncatt
{
    char name[1024];
    char* long_name;
    char* units;
    char* short_name;
    char* standard_name;
    request* metadata;
    nc_type nctype;
} ncatt_t;

typedef struct filter_type
{
    fieldset* fset;
    hypercube* filter;
    int count;
    double scale_factor;
    double add_offset;
    double missing;
    bool bitmap;
    ncatt_t att;
    request* filter_request;
    bool scale;
} dataset_t;

/*
 * typedef struct ncfile {
 *     dataset_t *filters;
 *     int ncid;
 * } ncfile_t;
*/

typedef struct ncoptions
{
    bool usevalidtime; /* Whether to use valid TIME only or not */
    bool auto_refdate; /* Automatic Reference Date */
    long refdate;         /* Reference date */
    const char* version;

    char* title;
    char* history;
    char* unlimited;
    bool checkvalidtime;
    request* mars_description;
    bool mmeans;      /* Whether this dataset is Monthly Means */
    bool climatology; /* Whether this dataset is climatology */
    bool shuffle;
    long deflate;
} ncoptions_t;

ncoptions_t setup;

#define NC_TYPES 7
struct nc_types_values
{
    double nc_type_max;
    double nc_type_min;
    double nc_type_missing;
} nc_type_values[NC_TYPES] = {
    /* In some occasions, SHRT_MIN-2 for the minimum value, makes ncview display
 missing values for -32766, while NC_FILL_SHORT=-32767, and SHRT_MIN=-32768 */
    { 0, 0, 0 },                                  /* NC_NAT,   'Not A Type' (c.f. NaN) */
    { 0x7f, NC_FILL_BYTE + 1, NC_FILL_BYTE },     /* NC_BYTE,   signed 1 byte integer */
    { 0xff, NC_FILL_CHAR + 1, NC_FILL_CHAR },     /* NC_CHAR,   ISO/ASCII character */
    { 0x7fff, NC_FILL_SHORT + 1, NC_FILL_SHORT }, /* NC_SHORT,  signed 2 byte integer */
    { 0x7ffffff, NC_FILL_INT + 1, NC_FILL_INT },  /* NC_INT,    signed 4 byte integer */
    { FLT_MAX, -FLT_MAX, NC_FILL_FLOAT },         /* NC_FLOAT,  single precision floating point number */
    { DBL_MAX, -DBL_MAX, NC_FILL_DOUBLE },        /* NC_DOUBLE, double precision floating point number */
};

static long fcmonth2days(long date, long months)
{
    long julianfrom = grib_date_to_julian(date);
    long years      = (long)(months / 12);
    long leap       = (long)(years / 4) - (long)(years / 100);
    long to         = years * 365 + (months % 12) * 32 + 1 + leap; /* FCMONTH can't be > 28 */
    long julianto   = julianfrom + to;
    long days       = 0;

    long dd       = date % 100;
    long dateto   = grib_julian_to_date(julianto);
    long nextdate = (dateto / 100) * 100 + dd;

    julianto = grib_date_to_julian(nextdate);
    days     = julianto - julianfrom;

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: date: %ld + %ld months = %ld days", date, months, days);

    return days;
}

static long request_fields(request* r)
{
    long cnt     = 1;
    parameter* p = r->params;
    while (p) {
        if (p->name[0] != '_') {
            cnt *= count_values(r, p->name);
        }
        p = p->next;
    }
    return cnt;
}

/* Transform MARS FCMONTHs into number of months from base date.
 For example, FCMONTH=1 is current month == 0 */
static void fcmonth2nbmonths(request* r)
{
    long n = count_values(r, "fcmonth");

    if (n == 0)
        return;
    n = atol(get_value(r, "fcmonth", 0));
    set_value(r, "fcmonth", "%ld", n - 1);
}

static long monthnumber(const char* m)
{
    const char* months[] = { "JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC" };
    int i                = 0;

    while (i < 12)
        if (strcmp(m, months[i++]) == 0)
            return i;
    grib_context_log(ctx, GRIB_LOG_ERROR, "Error. Translation for MONTH='%s' not found", m);
    return -1;
}

static int check_stepUnits(const char* step_units_str)
{
    /* Only hours, minutes and seconds supported */
    if (STR_EQUAL(step_units_str, "h") ||
        STR_EQUAL(step_units_str, "m") ||
        STR_EQUAL(step_units_str, "s")) {
        return GRIB_SUCCESS;
    }
    return GRIB_WRONG_STEP_UNIT;
}

/* The argument represents 1 field */
static void validation_time(request* r)
{
    long date        = 0;
    long time        = 0;
    double step      = 0;
    long fcmonthdays = 0;
    long fcmonth     = 0;
    double v;
    long julian            = 0;
    const char* step_units = NULL;

    long nstep    = count_values(r, "step");
    long ndate    = count_values(r, "date");
    long ntime    = count_values(r, "time");
    long nfcmonth = count_values(r, "fcmonth");

    static long julianrefdate = 0;

    if (nstep > 1 || ndate > 1 || ntime > 1 || nfcmonth > 1) {
        if (nstep > 1)
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error. Field has more than 1 STEP");
        if (ntime > 1)
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error. Field has more than 1 TIME");
        if (ndate > 1)
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error. Field has more than 1 DATE");
        if (nfcmonth > 1)
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error. Field has more than 1 FCMONTH");
        print_all_requests(r);
        exit(1);
    }

    if (nstep)
        step = atof(get_value(r, "step", 0));

    if (ndate) {
        const char* p         = get_value(r, "date", 0);
        const char* marsClass = get_value(r, "class", 0);
        if (eq_string(marsClass, "s2")) {
            /* S2S Data. See GRIB-699 and GRIB-762 */
            const char* hdate = get_value(r, "hdate", 0);
            grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: S2S Data");
            if (hdate) {
                p = hdate; /* This is a hindcast */
            }
        }
        if (is_number(p))
            date = atol(p);
        else {
            long second = 0;
            bool isjul, date_ok;
            date_ok = parsedate(p, &julian, &second, &isjul);
            if (!date_ok)
                grib_context_log(ctx, GRIB_LOG_ERROR, "Failed to parse date: '%s'", p);
            date = grib_julian_to_date(julian);
        }

        /* Climatology */
        if (strlen(p) == 3) {
            if (setup.usevalidtime)
                grib_context_log(ctx, GRIB_LOG_ERROR, "Climatology data. Setting usevalidtime=OFF");
            setup.auto_refdate = false;
            setup.usevalidtime = false;
            setup.climatology  = true;
        }
    }

    if (ntime)
        time = atol(get_value(r, "time", 0));

    if (nfcmonth) {
        fcmonth = atol(get_value(r, "fcmonth", 0));
        /* FCMONTH needs base DATE */
        if (date == 0) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error: FCMONTH needs DATE");
            exit(1);
        }
        fcmonthdays = fcmonth2days(date, fcmonth);
    }

    julian     = grib_date_to_julian(date);
    step_units = get_value(r, "stepUnits", 0);
    if (step_units) {
        if (check_stepUnits(step_units) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR,
                             "Cannot convert stepUnits of '%s'. Only hours, minutes and seconds supported.", step_units);
        }
        if (STR_EQUAL("m", step_units)) {
            step /= 60;
        }
        else if (STR_EQUAL("s", step_units)) {
            step /= 3600;
        }
    }
    v = julian * 24.0 + fcmonthdays * 24.0 + time / 100.0 + step * 1.0;
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: date=%ld, julian=%ld, fcmonthdays=%ld, time=%ld, step=%g, validation=%.3f", date, julian, fcmonthdays, time, step, v);
    set_value(r, "_validation", "%lf", v);
    set_value(r, "_juliandate", "%ld", julian);

    if (!julianrefdate)
        julianrefdate = grib_date_to_julian(setup.refdate);

    set_value(r, "_validationtime", "%lf", v - julianrefdate * 24.0);

    /* Remove minutes from TIME */
    if (ntime)
        set_value(r, "time", "%.1lf", time / 100.0);
}

static void free_nc_options()
{
    if (setup.title)
        grib_context_free(ctx, setup.title);
    if (setup.history)
        grib_context_free(ctx, setup.history);
    if (setup.unlimited)
        grib_context_free(ctx, setup.unlimited);
    if (setup.mars_description)
        free_all_requests(setup.mars_description);
}

static void get_nc_options(const request* user_r)
{
    const char* checkvalidtime_env = NULL;
    const char* validtime          = get_value(user_r, "usevalidtime", 0);
    const char* refdate            = get_value(user_r, "referencedate", 0);
    const char* shuffle            = get_value(user_r, "shuffle", 0);
    const char* deflate            = get_value(user_r, "deflate", 0);

    const char* title     = get_value(user_r, "title", 0);
    const char* history   = get_value(user_r, "history", 0);
    const char* unlimited = get_value(user_r, "unlimited", 0);

    setup.shuffle      = shuffle ? (strcmp(shuffle, "true") == 0) : false;
    setup.deflate      = deflate ? ((strcmp(deflate, "none") == 0) ? -1 : atol(deflate)) : -1;
    setup.usevalidtime = validtime ? (strcmp(validtime, "true") == 0) : false;
    setup.refdate      = refdate ? atol(refdate) : 19000101;
    setup.auto_refdate = refdate ? (strcmp(get_value(user_r, "referencedate", 0), "AUTOMATIC") == 0) : false;

    setup.title     = title ? grib_context_strdup(ctx, (title)) : NULL;
    setup.history   = history ? grib_context_strdup(ctx, (history)) : NULL;
    setup.unlimited = unlimited ? grib_context_strdup(ctx, ((unlimited))) : NULL;

    setup.checkvalidtime   = true;
    checkvalidtime_env     = getenv("GRIB_TO_NETCDF_CHECKVALIDTIME");
    if (checkvalidtime_env) {
        const long v = atol(checkvalidtime_env);
        if (v == 0) setup.checkvalidtime = false;
    }
    setup.mars_description = empty_request("MARS");
}

static nc_type translate_nctype(const char* name)
{
    if (!name)
        return NC_SHORT;

    if (STR_EQUAL(name, "NC_BYTE"))
        return NC_BYTE;

    if (STR_EQUAL(name, "NC_SHORT"))
        return NC_SHORT;

    if (STR_EQUAL(name, "NC_INT"))
        return NC_INT;

    if (STR_EQUAL(name, "NC_FLOAT"))
        return NC_FLOAT;

    if (STR_EQUAL(name, "NC_DOUBLE"))
        return NC_DOUBLE;

    grib_context_log(ctx, GRIB_LOG_ERROR, "Unknown netCDF type '%s'. Using NC_SHORT", name);
    return NC_SHORT;
}

static void check_err(const char* function, const int stat, const int line)
{
    if (stat != NC_NOERR) {
        /* (void) fprintf(stderr, "line %d of %s: %s\n", line, tool_name, nc_strerror(stat)); */
        (void)fprintf(stderr, "\n%s ERROR: line %d, %s: %s\n",
                      tool_name, line, function, nc_strerror(stat));
        if (stat == NC_EVARSIZE) {
            (void)fprintf(stderr,
                          "\nCannot create netCDF classic format, dataset is too large!\n"
                          "Try splitting the input GRIB(s).\n");
        }
        exit(1);
    }
}

#define DIM_ID 1
static int set_dimension(int ncid, const char* name, int n, int xtype, const char* units, const char* long_name)
{
    int var_id = 0;
    int stat   = 0;
    int dim_id = DIM_ID;
    int dim_vec[DIM_ID];

    if ( setup.unlimited && (STR_EQUAL(name, setup.unlimited)) )
        n = NC_UNLIMITED;

    stat = nc_def_dim(ncid, name, n, &dim_id);
    check_err("nc_def_dim", stat, __LINE__);

    dim_vec[0] = dim_id;
    stat       = nc_def_var(ncid, name, (nc_type)xtype, 1, dim_vec, &var_id);
    check_err("nc_def_var", stat, __LINE__);

    if (units != NULL) {
        stat = nc_put_att_text(ncid, var_id, "units", strlen(units), units);
        check_err("nc_put_att_text", stat, __LINE__);
    }

    if (long_name != NULL) {
        stat = nc_put_att_text(ncid, var_id, "long_name", strlen(long_name), long_name);
        check_err("nc_put_att_text", stat, __LINE__);
    }

    return var_id;
}

static int check_grid(field* f)
{
    err e = 0;
    char grid_type[80];
    size_t size = sizeof(grid_type);

    if ((e = grib_get_string(f->handle, "typeOfGrid", grid_type, &size)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get typeOfGrid %s", grib_get_error_message(e));
        return e;
    }

    if ( !STR_EQUAL(grid_type, "regular_ll") && !STR_EQUAL(grid_type, "regular_gg") ) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Grid type = %s", grid_type);
        grib_context_log(ctx, GRIB_LOG_ERROR, "First GRIB is not on a regular lat/lon grid or on a regular Gaussian grid. Exiting.\n");
        return GRIB_GEOCALCULUS_PROBLEM;
    }
    return e;
}

static int get_num_latitudes_longitudes(grib_handle* h, size_t* nlats, size_t* nlons)
{
    err e = 0;
    char grid_type[80];
    size_t size = sizeof(grid_type);

    if (grib_get_string(h, "typeOfGrid", grid_type, &size) == GRIB_SUCCESS &&
        strcmp(grid_type, "regular_ll") == 0) {
        /* Special shortcut for regular lat/on grids */
        long n;
        Assert(!grib_is_missing(h, "Ni", &e));
        if ((e = grib_get_long(h, "Ni", &n)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Ni: %s", grib_get_error_message(e));
            return e;
        }
        *nlons = n;

        if ((e = grib_get_long(h, "Nj", &n)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Nj: %s", grib_get_error_message(e));
            return e;
        }
        *nlats = n;
    }
    else {
        if ((e = grib_get_size(h, "distinctLatitudes", nlats)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get distinctLatitudes: %s", grib_get_error_message(e));
            return e;
        }
        if ((e = grib_get_size(h, "distinctLongitudes", nlons)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get distinctLongitudes: %s", grib_get_error_message(e));
            return e;
        }
    }
    return e;
}

static int def_latlon(int ncid, fieldset* fs)
{
    int n        = 0;
    size_t nlats = 0, nlons = 0;
    err e = 0;

    field* g = get_field(fs, 0, expand_mem);

    DEBUG_ASSERT(check_grid(g) == GRIB_SUCCESS);

    if ((e = get_num_latitudes_longitudes(g->handle, &nlats, &nlons)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get lat/lon info: %s", grib_get_error_message(e));
        return e;
    }

    /* Define longitude */
    n = (int)nlons;
    set_dimension(ncid, "longitude", n, NC_FLOAT, "degrees_east", "longitude");

    /* Define latitude */
    n = nlats;
    set_dimension(ncid, "latitude", n, NC_FLOAT, "degrees_north", "latitude");

    /* g->purge_header = true; */
    release_field(g);

    return e;
}

static int put_latlon(int ncid, fieldset* fs)
{
    int var_id = 0;
    size_t i   = 0;
    size_t n   = 0;
    int stat   = 0;
    err e      = 0;

    field* g = get_field(fs, 0, expand_mem);

    double* dvalues = NULL;
    float* fvalues  = NULL;
    long nv         = 0;

    size_t ni;
    size_t nj;

#if 0
    /* Get info in degrees */
    if((e = grib_get_double(g->handle, "iDirectionIncrementInDegrees", &ew_stride)) != GRIB_SUCCESS)
    {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get iDirectionIncrementInDegrees: %s", grib_get_error_message(e));
        return e;
    }

    if((e = grib_get_double(g->handle, "jDirectionIncrementInDegrees", &ns_stride)) != GRIB_SUCCESS)
    {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get jDirectionIncrementInDegrees %s", grib_get_error_message(e));
        return e;
    }

    /* Define longitude */
    if((e = grib_get_long(g->handle, "Ni", &ni)) != GRIB_SUCCESS)
    {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Ni %s", grib_get_error_message(e));
        return e;
    }
    /* Define latitude */
    if((e = grib_get_long(g->handle, "Nj", &nj)) != GRIB_SUCCESS)
    {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Nj %s", grib_get_error_message(e));
        return e;
    }

#endif
    if ((e = get_num_latitudes_longitudes(g->handle, &nj, &ni)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: put_latlon: cannot get distinctLatitudes: %s", grib_get_error_message(e));
        return e;
    }

    /* Compute max. # values and allocate */
    nv = ni;
    if (nv < nj)
        nv = nj;

    fvalues = (float*)grib_context_malloc(ctx, sizeof(float) * nv);
    dvalues = (double*)grib_context_malloc(ctx, sizeof(double) * nv);

    /* longitude */
    n    = ni;
    stat = nc_inq_varid(ncid, "longitude", &var_id);
    check_err("nc_inq_varid", stat, __LINE__);
    if ((e = grib_get_double_array(g->handle, "distinctLongitudes", dvalues, &n)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: put_latlon: cannot get distinctLongitudes: %s", grib_get_error_message(e));
        return e;
    }
    if (n != ni) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Number of distinctLongitudes is not the same as Ni (%zu!=%zu)",n,ni);
        return GRIB_GEOCALCULUS_PROBLEM;
    }

    for (i = 0; i < n; i++) {
        fvalues[i] = dvalues[i];
    }

    stat = nc_put_var_float(ncid, var_id, fvalues);
    check_err("nc_put_var_float", stat, __LINE__);

    /* latitude */
    n    = nj;
    stat = nc_inq_varid(ncid, "latitude", &var_id);
    check_err("nc_inq_varid", stat, __LINE__);
    if ((e = grib_get_double_array(g->handle, "distinctLatitudes", dvalues, &n)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: put_latlon: cannot get distinctLatitudes: %s", grib_get_error_message(e));
        return e;
    }
    if (n != nj) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Number of distinctLatitudes is not the same as Nj (%zu!=%zu)",n,nj);
        return GRIB_GEOCALCULUS_PROBLEM;
    }

    for (i = 0; i < n; i++) {
        fvalues[i] = dvalues[i];
    }
    stat = nc_put_var_float(ncid, var_id, fvalues);
    check_err("nc_put_var_float", stat, __LINE__);

    /* g->purge_header = true; */
    release_field(g);
    grib_context_free(ctx, fvalues);
    grib_context_free(ctx, dvalues);

    return e;
}

static int compute_scale(dataset_t* subset)
{
    double max            = -DBL_MAX;
    double min            = DBL_MAX;
    double median         = 0;
    long i                = 0;
    size_t j              = 0;
    int64_t scaled_max    = 0;
    int64_t scaled_min    = 0;
    int64_t scaled_median = 0;
    double ao = 0.0, sf = 0.0;
    double x;

    char test_scaled_max    = 0;
    char test_scaled_min    = 0;
    char test_scaled_median = 0;

    err e = 0;

    fieldset* fs = subset->fset;
    int idx      = subset->att.nctype;

    for (i = 0; i < fs->count; i++) {
        field* g = get_field(fs, i, expand_mem);
        size_t len;
        static double* vals    = NULL;
        static size_t vals_len = 0;

        if ((e = grib_get_size(g->handle, "values", &len)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get size of values: %s", grib_get_error_message(e));
            return e;
        }

        if (len > vals_len) {
            if (vals)
                grib_context_free(ctx, vals);
            vals     = (double*)grib_context_malloc(ctx, sizeof(double) * len);
            vals_len = len;
        }
        if ((e = grib_get_double_array(g->handle, "values", vals, &len)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get values: %s", grib_get_error_message(e));
            return e;
        }

        if (g->has_bitmap) {
            subset->bitmap = true;
            for (j = 0; j < len; ++j) {
                if (vals && vals[j] != global_missing_value) {
                    if (vals[j] > max) max = vals[j];
                    if (vals[j] < min) min = vals[j];
                }
            }
        }
        else {
            for (j = 0; j < len; ++j) {
                if (vals) {
                    if (vals[j] > max) max = vals[j];
                    if (vals[j] < min) min = vals[j];
                }
            }
        }
        /* g->purge_header = true; */
        release_field(g);
    }

    median = (max + min) / 2.0;

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: max_int: %lf, min_int: %lf", nc_type_values[idx].nc_type_max, nc_type_values[idx].nc_type_min);

    sf = (double)((max - min) / (double)(nc_type_values[idx].nc_type_max - nc_type_values[idx].nc_type_min));
    ao = ((max + min) - sf * (nc_type_values[idx].nc_type_min + nc_type_values[idx].nc_type_max)) / 2;

    if (min == max) {
        sf = 1.0; /* Prevent divide by zero later. Constant field grib has max == min */
    }

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: idx is: %d", idx);
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: max: %lf, min: %lf, median: %lf, scale factor: %lf, add_offset: %lf", max, min, median, sf, ao);

    x = ((median - ao));
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: x=%lf", x);
    x /= sf;
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: x=%lf", x);

    scaled_max    = rint((max - ao) / sf);
    scaled_min    = rint((min - ao) / sf);
    scaled_median = rint((median - ao) / sf);

    if (scaled_max > nc_type_values[idx].nc_type_max) {
        grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: scaled_max (=%lld) > nc_type_max (=%lf). Set sf to 1.0",
                         (long long)scaled_max, nc_type_values[idx].nc_type_max);
        sf = 1.0; /* ECC-685 */
    }

    test_scaled_max    = (char)scaled_max;
    test_scaled_min    = (char)scaled_min;
    test_scaled_median = (char)scaled_median;

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: scaled_max: %lld, scaled_min: %lld, scaled_median: %lld, x: %lf",
                     (long long)scaled_max, (long long)scaled_min, (long long)scaled_median, x);

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: test_scaled_max: %x, test_scaled_min: %x, test_scaled_median: %x",
                     test_scaled_max, test_scaled_min, test_scaled_median);

    max    = scaled_max * sf + ao;
    min    = scaled_min * sf + ao;
    median = scaled_median * sf + ao;

    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: max: %lf, min: %lf, median: %lf", max, min, median);

    subset->scale_factor = sf;
    subset->add_offset   = ao;

    return 0;
}

static int nc_put_att_type(int ncid, int varid, const char* name, nc_type nctype, int n, double value)
{
    int r = 0;
    switch (nctype) {
        case NC_BYTE: {
            unsigned char val_uchar = (unsigned char)value;
            r                       = nc_put_att_uchar(ncid, varid, name, nctype, n, &val_uchar);
            break;
        }
        case NC_SHORT: {
            short int val_short = (short int)value;
            r                   = nc_put_att_short(ncid, varid, name, nctype, n, &val_short);
            break;
        }
        case NC_INT: {
            int val_int = (int)value;
            r           = nc_put_att_int(ncid, varid, name, nctype, n, &val_int);
            break;
        }
        case NC_FLOAT: {
            float val_flt = (float)value;
            r             = nc_put_att_float(ncid, varid, name, nctype, n, &val_flt);
            break;
        }
        case NC_DOUBLE: {
            double val_dbl = (double)value;
            r              = nc_put_att_double(ncid, varid, name, nctype, n, &val_dbl);
            break;
        }
        default:
            grib_context_log(ctx, GRIB_LOG_ERROR, "nc_put_att_type(...): Unknown netCDF type '%d'", nctype);
            break;
    }
    return r;
}

static int nc_put_vara_type(int ncid, int varid, const size_t start[], const size_t count[], void* valuesp, nc_type nctype)
{
    int r = 0;
    switch (nctype) {
        case NC_BYTE:
            r = nc_put_vara_uchar(ncid, varid, start, count, (unsigned char*)valuesp);
            break;
        case NC_SHORT:
            r = nc_put_vara_short(ncid, varid, start, count, (short int*)valuesp);
            break;
        case NC_INT:
            r = nc_put_vara_int(ncid, varid, start, count, (int*)valuesp);
            break;
        case NC_FLOAT:
            r = nc_put_vara_float(ncid, varid, start, count, (float*)valuesp);
            break;
        case NC_DOUBLE:
            r = nc_put_vara_double(ncid, varid, start, count, (double*)valuesp);
            break;
        default:
            grib_context_log(ctx, GRIB_LOG_ERROR, "nc_put_vara_type(...): Unknown netCDF type '%d'", nctype);
            break;
    }
    return r;
}

static void scale_bitmap(double* vals, long n, void* data, dataset_t* subset)
{
    int i          = 0;
    nc_type nctype = subset->att.nctype;

    /* if(!subset->bitmap) {
        grib_context_log(ctx,GRIB_LOG_DEBUG,"No scale of bitmap required");
        return;
     } */
    if (n > 0 && !vals) {
        Assert(!"scale_bitmap: n > 0 but vals == NULL");
        return;
    }

    switch (nctype) {
        case NC_BYTE: {
            unsigned char* vscaled = (unsigned char*)data;
            for (i = 0; i < n; ++i) {
                if (vals[i] == global_missing_value) {
                    vscaled[i] = (unsigned char)subset->missing;
                }
            }
            break;
        }

        case NC_SHORT: {
            short int* vscaled = (short int*)data;
            for (i = 0; i < n; ++i) {
                if (vals[i] == global_missing_value) {
                    vscaled[i] = (short int)subset->missing;
                }
            }
            break;
        }

        case NC_INT: {
            int* vscaled = (int*)data;
            for (i = 0; i < n; ++i) {
                if (vals[i] == global_missing_value) {
                    vscaled[i] = (int)subset->missing;
                }
            }
            break;
        }

        case NC_FLOAT: {
            float* vscaled = (float*)data;
            for (i = 0; i < n; ++i) {
                if (vals[i] == global_missing_value) {
                    vscaled[i] = (float)subset->missing;
                }
            }
            break;
        }

        case NC_DOUBLE: {
            double* vscaled = (double*)data;
            for (i = 0; i < n; ++i) {
                if (vals[i] == global_missing_value) {
                    vscaled[i] = (double)subset->missing;
                }
            }
            break;
        }

        default:
            grib_context_log(ctx, GRIB_LOG_ERROR, "scale(...): Unknown netCDF type %d", nctype);
            break;
    }
}

static void scale(double* vals, long n, void* data, dataset_t* g)
{
    int i               = 0;
    nc_type nctype      = g->att.nctype;
    double scale_factor = g->scale_factor;
    double add_offset   = g->add_offset;

    /*
    if(!subset->scale)
    {
     grib_context_log(ctx,GRIB_LOG_DEBUG,"No scale required");
     return;
    }
    */
    DEBUG_ASSERT(vals);
    DEBUG_ASSERT(n > 0);
    if (!vals) return;

    switch (nctype) {
        case NC_BYTE: {
            unsigned char* vscaled = (unsigned char*)data;
            for (i = 0; i < n; ++i) {
                if (!g->bitmap || (vals[i] != global_missing_value)) {
                    double d = rint((vals[i] - add_offset) / scale_factor);
                    if (!(d >= nc_type_values[nctype].nc_type_min && d <= nc_type_values[nctype].nc_type_max)) {
                        grib_context_log(ctx, GRIB_LOG_ERROR, "Scaling for type NC_BYTE failed");
                        return;
                    }
                    vscaled[i] = d;
                }
            }
            break;
        }

        case NC_SHORT: {
            short int* vscaled = (short int*)data;
            for (i = 0; i < n; ++i) {
                if (!g->bitmap || (vals[i] != global_missing_value)) {
                    double d = 0;
                    Assert(scale_factor > 0);
                    d = rint((vals[i] - add_offset) / scale_factor);
                    if (!(d >= nc_type_values[nctype].nc_type_min && d <= nc_type_values[nctype].nc_type_max)) {
                        grib_context_log(ctx, GRIB_LOG_ERROR, "Scaling for type NC_SHORT failed");
                        return;
                    }
                    vscaled[i] = d;
                }
            }
            break;
        }

        case NC_INT: {
            int* vscaled = (int*)data;
            for (i = 0; i < n; ++i) {
                if (!g->bitmap || (vals[i] != global_missing_value)) {
                    double d = rint((vals[i] - add_offset) / scale_factor);
                    if (!(d >= nc_type_values[nctype].nc_type_min && d <= nc_type_values[nctype].nc_type_max)) {
                        grib_context_log(ctx, GRIB_LOG_ERROR, "Scaling for type NC_INT failed");
                        return;
                    }
                    vscaled[i] = d;
                }
            }
            break;
        }

        case NC_FLOAT: {
            float* vscaled = (float*)data;
            for (i = 0; i < n; ++i) {
                if (!g->bitmap || (vals[i] != global_missing_value)) {
                    double d = vals[i];
                    if (!(d >= nc_type_values[nctype].nc_type_min && d <= nc_type_values[nctype].nc_type_max)) {
                        grib_context_log(ctx, GRIB_LOG_ERROR, "Scaling for type NC_FLOAT failed");
                        return;
                    }
                    vscaled[i] = d;
                }
            }
            break;
        }

        case NC_DOUBLE: {
            double* vscaled = (double*)data;
            for (i = 0; i < n; ++i) {
                if (!g->bitmap || (vals[i] != global_missing_value)) {
                    double d = vals[i];
                    if(!(d >= nc_type_values[nctype].nc_type_min && d <= nc_type_values[nctype].nc_type_max)) {
                        grib_context_log(ctx, GRIB_LOG_ERROR, "Scaling for type NC_DOUBLE failed");
                        return;
                    }
                    vscaled[i] = d;
                }
            }
            break;
        }

        default:
            grib_context_log(ctx, GRIB_LOG_ERROR, "scale(...): Unknown netCDF type %d", nctype);
            break;
    }
}

/* Return array of strings which are the "time" values */
static char** create_times_array(const request* cube, size_t* size)
{
    char** result         = NULL;
    const char* time_axis = "time"; /* special case */
    parameter* the_param  = find_parameter(cube, time_axis);
    *size                 = 0;
    if (the_param) {
        size_t i = 0, num_values = 0;
        value* va = NULL;

        if (!the_param->count)
            count_values(cube, time_axis);

        /* Go thru all values to count how many there are */
        va = the_param->values;
        while (va) {
            ++num_values;
            va = va->next;
        }
        /* Create and populate array */
        result = (char**)grib_context_malloc(ctx, sizeof(char*) * num_values);
        va     = the_param->values;
        while (va) {
            result[i++] = va->name;
            va          = va->next;
        }
        *size = num_values;
    }
    return result;
}

static int put_data(hypercube* h, int ncid, const char* name, dataset_t* subset)
{
    int i      = 0;
    int stat   = 0;
    int dataid = 0;
    int naxis  = count_axis(h);
    size_t start[NC_MAX_DIMS];
    size_t count[NC_MAX_DIMS];
    char** times_array      = NULL;
    size_t times_array_size = 0;
    fieldset* fs            = subset->fset;
    field* f                = get_field(fs, 0, expand_mem);

    void* vscaled       = NULL;
    size_t vscaled_length = 0;

    long ni;
    long nj;
    err e = 0;

    /* Define longitude */
    if ((e = grib_get_long(f->handle, "Ni", &ni)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Ni: %s", grib_get_error_message(e));
        return e;
    }
    /* Define latitude */
    if ((e = grib_get_long(f->handle, "Nj", &nj)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Nj: %s", grib_get_error_message(e));
        return e;
    }

    /* Start filling dimensions at first value */
    for (i = 0; i < 2 + naxis; ++i)
        start[i] = 0;

    /* Count dimensions per axis */
    for (i = 0; i < naxis; ++i)
        count[naxis - i - 1] = 1;

    count[naxis]     = nj; /* latitude */
    count[naxis + 1] = ni; /* longitude */

    /* f->purge_header = true; */
    release_field(f);

    stat = nc_inq_varid(ncid, name, &dataid);
    check_err("nc_inq_varid", stat, __LINE__);

    /* GRIB-792: Build fast array storing values for the "time" axis. */
    /* This is for performance reasons */
    times_array = create_times_array(h->cube, &times_array_size);

    for (i = 0; i < fs->count; i++) {
        field* g = get_field(fs, i, expand_mem);
        size_t len;
        static double* vals    = NULL;
        static size_t vals_len = 0;
        bool missing        = 0;

        request* r;
        int j = 0;
        int idx[1024];
        int idxsize = 1024;

        if ((e = grib_get_size(g->handle, "values", &len)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get size of values: %s", grib_get_error_message(e));
            return e;
        }

        if (len > vals_len) {
            if (vals)
                grib_context_free(ctx, vals);
            vals     = (double*)grib_context_malloc(ctx, sizeof(double) * len);
            vals_len = len;
        }
        if ((e = grib_get_double_array(g->handle, "values", vals, &len)) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get values: %s", grib_get_error_message(e));
            return e;
        }

        // bool missing = (g->ksec4[0] < 0); /* If negative number of values, field is missing */

        r = field_to_request(g);
        if (!missing) {
            /* Reserved the maximum memory needed */
            /* This should only be done once, as all fields have the same geometry */
            if ((vscaled_length == 0) || (vscaled_length < sizeof(double) * len)) {
                if (vscaled)
                    grib_context_free(ctx, vscaled);
                vscaled        = (void*)grib_context_malloc(ctx, sizeof(double) * len);
                vscaled_length = sizeof(double) * len;
            }

            scale(vals, len, vscaled, subset);
            if (subset->bitmap)
                scale_bitmap(vals, len, vscaled, subset);

            if ((e = grib_get_long(g->handle, "Ni", &ni)) != GRIB_SUCCESS) {
                grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Ni: %s", grib_get_error_message(e));
                return e;
            }
            /* Define latitude */
            if ((e = grib_get_long(g->handle, "Nj", &nj)) != GRIB_SUCCESS) {
                grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Nj: %s", grib_get_error_message(e));
                return e;
            }

            if (nj != count[naxis] || ni != count[naxis + 1]) {
                grib_context_log(ctx, GRIB_LOG_ERROR, "GRIB message %d has different resolution\n", i + 1);
                grib_context_log(ctx, GRIB_LOG_ERROR, "lat=%ld, long=%ld instead of lat=%ld, long=%ld\n", nj, ni, count[naxis], count[naxis + 1]);
                exit(1);
            }

            cube_indexes(h, r, times_array, times_array_size, idx, idxsize);
            for (j = 0; j < naxis; ++j)
                start[naxis - j - 1] = idx[j];

            grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: Put data from field %d", i);

            stat = nc_put_vara_type(ncid, dataid, start, count, vscaled, subset->att.nctype);
            check_err("nc_put_vara_type", stat, __LINE__);
        }

        /* g->purge_header = true; */
        release_field(g);
    }

    grib_context_free(ctx, vscaled);
    grib_context_free(ctx, times_array);
    return 0;
}

static void set_always_a_time(hypercube* h, request* data_r)
{
    if (setup.usevalidtime && count_values(data_r, "time") == 1) {
        set_value(h->cube, "time", "%.2lf", atof(get_value(data_r, "_validationtime", 0)));
        add_value(h->cube, "axis", "time");

        {
            int i = 0;
            int n = count_values(h->cube, "axis");
            if (n) {
                h->compare = (namecmp*)calloc(n, sizeof(namecmp));
                Assert(h->compare);
            }

            for (i = 0; i < n; i++)
                h->compare[i] = comparator(get_value(h->cube, "axis", i));
        }
    }
}

static int define_netcdf_dimensions(hypercube* h, fieldset* fs, int ncid, dataset_t* subsets, int subsetcnt, const request* data_r)
{
    const request* cube = h->cube;
    int naxis           = count_axis(h);
    int i               = 0;
    int stat            = 0;
    int n               = 0;
    int var_id          = 0; /* Variable ID */
    int dims[1024];

    size_t chunks[NC_MAX_DIMS] = {0,}; /* For chunking */
    err e = 0;

    long ni;
    long nj;

    field* f = get_field(fs, 0, expand_mem);

    if ((e = check_grid(f)) != GRIB_SUCCESS) {
        release_field(f);
        return e;
    }

    if ((e = grib_get_long(f->handle, "Ni", &ni)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Ni: %s", grib_get_error_message(e));
        return e;
    }
    if ((e = grib_get_long(f->handle, "Nj", &nj)) != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "ecCodes: cannot get Nj: %s", grib_get_error_message(e));
        return e;
    }
    release_field(f);

    /* Count dimensions per axis */
    for (i = 0; i < naxis; ++i)
        chunks[naxis - i - 1] = 1;

    chunks[naxis]     = nj; /* latitude */
    chunks[naxis + 1] = ni; /* longitude */

    /* START DEFINITIONS */

    /* Define latitude/longitude dimensions */
    e = def_latlon(ncid, fs);
    if (e != GRIB_SUCCESS)
        return e;

    /* Define netCDF dimensions */
    for (i = 0; i < naxis; ++i) {
        int nctype       = NC_INT;
        const char* axis = get_axis(h, i);
        const char* units      = NULL;
        char u[10240];
        const char* lowaxis = (axis);
        const char* longname = (char*)lowaxis;
        n                   = count_values(cube, axis);

        if (count_values(data_r, "levtype") > 1) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Cannot handle fields for different levtypes.\n");
            grib_context_log(ctx, GRIB_LOG_ERROR, "Please split input data into different files. Exiting!\n");
            exit(1);
        }
        if (count_values(data_r, "stepUnits") > 1) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Cannot handle fields for different stepUnits.\n");
            grib_context_log(ctx, GRIB_LOG_ERROR, "Please split input data into different files. Exiting!\n");
            exit(1);
        }

        if (strcmp(axis, "levelist") == 0) {
            const char* levtype = get_value(data_r, "levtype", 0);
            if (levtype) {
                if (strcmp("pl", levtype) == 0) {
                    units    = "millibars";
                    longname = "pressure_level";
                }

                if (strcmp("ml", levtype) == 0) {
                    longname = "model_level_number";
                }
            }

            lowaxis = "level";
        }

        if (strcmp(axis, "number") == 0) {
            longname = "ensemble_member";
        }

        if (strcmp(axis, "date") == 0) {
            snprintf(u, sizeof(u), "days since %ld-%02ld-%02ld 00:00:0.0", setup.refdate / 10000, (setup.refdate % 10000) / 100, (setup.refdate % 100));
            units    = u;
            longname = "Base_date";
            if (setup.climatology) {
                snprintf(u, sizeof(u), "months");
                units = u;
            }
        }

        if (strcmp(axis, "time") == 0) {
            bool onedtime = (count_values(cube, "date") == 0 && count_values(cube, "step") == 0);
            snprintf(u, sizeof(u), "hours since 0000-00-00 00:00:00.0");
            longname = "reference_time";
            if (setup.usevalidtime || onedtime) {
                snprintf(u, sizeof(u), "hours since %ld-%02ld-%02ld 00:00:00.0", setup.refdate / 10000, (setup.refdate % 10000) / 100, (setup.refdate % 100));
                longname = "time";
            }
            if (setup.climatology) {
                snprintf(u, sizeof(u), "hours");
            }
            units = u;
            /* nctype = NC_FLOAT; */
        }

        if (strcmp(axis, "step") == 0) {
            units    = "hours";
            longname = "time_step";
            if (count_values(cube, "date") == 0 && count_values(cube, "time") == 0) {
                const char* d = get_value(data_r, "date", 0);
                const char* t = get_value(data_r, "time", 0);
                long date     = d ? atol(d) : 0;
                long hour     = t ? atol(t) : 0;
                long min      = t ? 60 * (atof(t) - hour) : 0;
                snprintf(u, sizeof(u), "hours since %ld-%02ld-%02ld %02ld:%02ld:00.0", date / 10000, (date % 10000) / 100, (date % 100), hour, min);
                units = u;
            }
        }

        if ((strcmp(axis, "fcmonth") == 0)) {
            const char* date = get_value(data_r, "date", 0);
            char ymd[32]     = "";

            if (date) {
                strncat(ymd, date, 4);
                strcat(ymd, "-");
                strncat(ymd, date + 4, 2);
                strcat(ymd, "-");
                /* udunits is a bit tricky with month being 30.4 days */
                /* ncview doesn't display properly */
                strcat(ymd, "01");
            }
            else {
                snprintf(ymd, sizeof(ymd), "00-00-00");
            }
            snprintf(u, sizeof(u), "months since %s 00:00:00.0", ymd);
            units    = u;
            longname = "time";
        }

        var_id = set_dimension(ncid, lowaxis, n, nctype, units, longname);

        if (strcmp(axis, "time") == 0) {
            if (setup.usevalidtime) {
                const char* cal = "gregorian";

                if (setup.mmeans) {
                    const char* period = "0000-01-00 00:00:00";
                    stat         = nc_put_att_text(ncid, var_id, "avg_period", strlen(period), period);
                    check_err("nc_put_att_text", stat, __LINE__);
                }

                stat = nc_put_att_text(ncid, var_id, "calendar", strlen(cal), cal);
                check_err("nc_put_att_text", stat, __LINE__);
            }
        }
    }

    /* Define data dimension */
    n = 1 + 1 + naxis; /* longitude + latitude + # axis */
    for (i = 0; i < n; ++i)
        dims[i] = n - i - 1;

    for (i = 0; i < subsetcnt; ++i) {
        printf("%s: Defining variable '%s'.\n", tool_name, subsets[i].att.name);

        stat = nc_def_var(ncid, subsets[i].att.name, subsets[i].att.nctype, n, dims, &var_id);
        check_err("nc_def_var", stat, __LINE__);

        if (setup.deflate > -1) {
#ifdef NC_NETCDF4
            stat = nc_def_var_chunking(ncid, var_id, NC_CHUNKED, chunks);
            check_err("nc_def_var_chunking", stat, __LINE__);

            /* Set compression settings for a variable */
            stat = nc_def_var_deflate(ncid, var_id, setup.shuffle, 1, setup.deflate);
            check_err("nc_def_var_deflate", stat, __LINE__);
#else
            (void)chunks;
            grib_context_log(ctx, GRIB_LOG_ERROR, "Deflate option only supported in netCDF4");
#endif
        }
        if (subsets[i].scale) {
            compute_scale(&subsets[i]);
            stat = nc_put_att_double(ncid, var_id, "scale_factor", NC_DOUBLE, 1, &subsets[i].scale_factor);
            check_err("nc_put_att_double", stat, __LINE__);

            stat = nc_put_att_double(ncid, var_id, "add_offset", NC_DOUBLE, 1, &subsets[i].add_offset);
            check_err("nc_put_att_double", stat, __LINE__);
        }

        stat = nc_put_att_type(ncid, var_id, "_FillValue", subsets[i].att.nctype, 1, nc_type_values[subsets[i].att.nctype].nc_type_missing);
        check_err("nc_put_att_type", stat, __LINE__);
        stat = nc_put_att_type(ncid, var_id, "missing_value", subsets[i].att.nctype, 1, nc_type_values[subsets[i].att.nctype].nc_type_missing);
        check_err("nc_put_att_type", stat, __LINE__);

        if (subsets[i].att.units) {
            const char* txt = subsets[i].att.units;
            stat            = nc_put_att_text(ncid, var_id, "units", strlen(txt), txt);
            check_err("nc_put_att_text", stat, __LINE__);
        }

        if (subsets[i].att.long_name) {
            const char* txt = subsets[i].att.long_name;
            stat            = nc_put_att_text(ncid, var_id, "long_name", strlen(txt), txt);
            check_err("nc_put_att_text", stat, __LINE__);
        }

        if (subsets[i].att.short_name) {
            const char* txt = subsets[i].att.short_name;
            stat            = nc_put_att_text(ncid, var_id, "short_name", strlen(txt), txt);
            check_err("nc_put_att_text", stat, __LINE__);
        }

        if (subsets[i].att.standard_name) {
            const char* txt = subsets[i].att.standard_name;
            stat            = nc_put_att_text(ncid, var_id, "standard_name", strlen(txt), txt);
            check_err("nc_put_att_text", stat, __LINE__);
        }

        // if(subsets[i].att.other)
        // {
        //     const char *txt = subsets[i].att.long_name;
        //     stat = nc_put_att_text(ncid, var_id, "other",strlen(txt),txt);
        //     check_err("nc_put_att_text", stat,__LINE__,__FILE__);
        // }

        if (subsets[i].att.metadata) {
            parameter* p = subsets[i].att.metadata->params;

            while (p) {
                const char* txt = p->values->name;
                stat            = nc_put_att_text(ncid, var_id, p->name, strlen(txt), txt);
                check_err("nc_put_att_text", stat, __LINE__);

                p = p->next;
            }
        }
    }

    /* Dimension-less variable for MARS request */
    if ((0)) /* reset when we have proper & fast mars_description */
    {
        /* parameter *p = data_r->params; */
        parameter* p = setup.mars_description->params;
        stat         = nc_def_var(ncid, "MARS", NC_CHAR, 0, 0, &var_id);
        check_err("nc_def_var", stat, __LINE__);

        /* Store request for those parameters with single value */
        while (p) {
            /* if((p->name[0] != '_') && (p->count == 1)) */
            if (p->name[0] != '_') {
                char par[1024];
                char val[1024000] = "";
                snprintf(par, sizeof(par), "%s", (p->name));
#if 0
                value2string(p->values,val);
#else
                snprintf(val, sizeof(val), "%s", (p->values->name));
#endif
                stat = nc_put_att_text(ncid, var_id, par, strlen(val), (val));
                if (stat != NC_NOERR) {
                    printf("Error setting request for %s = %s\n", par, val);
                }
                check_err("nc_put_att_text", stat, __LINE__);
            }
            p = p->next;
        }
    }

    /* Global attributes */
    {
        char timestamp[80];
        time_t now;
        /* char *convention = "MARS;CF"; */
        const char* convention = "CF-1.6";
        char history[10240];
        /* char *institution = "ECMWF Meteorological Archival and Retrieval System"; */

        /* Convention */
        stat = nc_put_att_text(ncid, NC_GLOBAL, "Conventions", strlen(convention), convention);
        check_err("nc_put_att_text", stat, __LINE__);

        /* Use history provided or Timestamp */
        if (setup.history) {
            snprintf(history, sizeof(history), "%s", setup.history);
        }
        else {
            int major    = ECCODES_MAJOR_VERSION;
            int minor    = ECCODES_MINOR_VERSION;
            int revision = ECCODES_REVISION_VERSION;

            time(&now);
            strftime(timestamp, sizeof(timestamp), "%Y-%m-%d %H:%M:%S GMT", gmtime(&now));
            snprintf(history, sizeof(history), "%s by grib_to_netcdf-%d.%d.%d: %s", timestamp, major, minor, revision, argvString);
        }
        stat = nc_put_att_text(ncid, NC_GLOBAL, "history", strlen(history), history);
        check_err("nc_put_att_text", stat, __LINE__);

        //stat = nc_put_att_text(ncid, NC_GLOBAL, "source",strlen(setup.source),setup.source);
        //check_err(stat,__LINE__,__FILE__);
        //stat = nc_put_att_text(ncid, NC_GLOBAL, "institution",strlen(institution),institution);
        //check_err(stat,__LINE__,__FILE__);

        if (setup.title) {
            stat = nc_put_att_text(ncid, NC_GLOBAL, "title", strlen(setup.title), setup.title);
            check_err("nc_put_att_text", stat, __LINE__);
        }
    }
    return e;
}

static size_t string_to_unique_number(const char* axis, const char* str)
{
    size_t result = 0;
    if (strcmp(axis, "type") == 0) {
        /* TODO: not ideal but capture the most common MARS types */
        if (strcmp(str, "an") == 0)
            return 2;
        else if (strcmp(str, "fc") == 0)
            return 9;
        else if (strcmp(str, "cf") == 0)
            return 10;
        else if (strcmp(str, "pf") == 0)
            return 11;
        else if (strcmp(str, "em") == 0)
            return 17;
        else if (strcmp(str, "es") == 0)
            return 18;
        else if (strcmp(str, "ep") == 0)
            return 30;
        else if (strcmp(str, "4i") == 0)
            return 33;
        else if (strcmp(str, "4g") == 0)
            return 8;
        else if (strcmp(str, "ia") == 0)
            return 3;
        else if (strcmp(str, "efi") == 0)
            return 27;
    }
    /* Fallback general case: Use hashing */
    result = 5381;
    while (*str) {
        result = 33 * result ^ (unsigned char)*str++;
    }
    return result;
}

static int fill_netcdf_dimensions(hypercube* h, fieldset* fs, int ncid)
{
    const request* cube = h->cube;
    int naxis           = count_axis(h);
    int i               = 0;
    int var_id          = 0;
    int stat            = 0;

    /* Put latitude/longitude values */
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: Put latitude/longitude values");
    put_latlon(ncid, fs);

    /* Put netCDF axis values */
    grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: Put netcdf axis values");
    for (i = 0; i < naxis; ++i) {
        int j               = 0;
        const char* axis    = get_axis(h, i);
        int n               = count_values(cube, axis);
        int* values         = (int*)grib_context_malloc(ctx, sizeof(int) * n);
        const char* lowaxis = (axis);

        if (!values) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "fill_netcdf_dimensions: cannot allocate %ld bytes", sizeof(int) * n);
            exit(1);
        }

        if (strcmp("levelist", axis) == 0)
            lowaxis = "level";

        if (strcmp("date", axis) == 0) {
            if (setup.climatology)
                for (j = 0; j < n; ++j)
                    values[j] = monthnumber(get_value(cube, axis, j));
            else
                for (j = 0; j < n; ++j)
                    values[j] = grib_date_to_julian(atol(get_value(cube, axis, j))) - grib_date_to_julian(setup.refdate);
        }
        else {
            for (j = 0; j < n; ++j) {
                long lv        = 0;
                const char* sv = get_value(cube, axis, j);
                if (is_number(sv)) {
                    lv = atol(sv); /* Detect error? */
                }
                else {
                    /* ECC-725: Convert string-valued dimension to integer
                     * e.g. mars type or stream */
                    lv = string_to_unique_number(axis, sv);
                }
                values[j] = lv;
            }
        }

        stat = nc_inq_varid(ncid, (lowaxis), &var_id);
        check_err("nc_inq_varid", stat, __LINE__);

        /* if ( strcmp("time", axis) == 0 && setup.unlimited != NULL && strcmp(setup.unlimited, "time") == 0 && setup.usevalidtime) */
        /* GRIB-437, GRIB-625 Special treatment of RECORD (unlimited) dimension */
        /* See "The NetCDF C Interface Guide" Section 6.23 */
        if (setup.unlimited != NULL && strcmp(setup.unlimited, axis) == 0) {
            /* This is tricky. I'm not sure it works when this dimension is not outer dimension */
            size_t start[NC_MAX_DIMS];
            size_t count[NC_MAX_DIMS];
            nc_type dim_type = 0;

            start[0] = 0;
            count[0] = n;

            stat = nc_inq_vartype(ncid, var_id, &dim_type); /* get the type of this dimension */
            check_err("nc_inq_vartype", stat, __LINE__);

            stat = nc_put_vara_type(ncid, var_id, start, count, values, dim_type);
            check_err("nc_put_vara_type", stat, __LINE__);
        }
        else {
            stat = nc_put_var_int(ncid, var_id, values);
            check_err("nc_put_var_int", stat, __LINE__);
        }

        grib_context_free(ctx, values);
    }

    return 0;
}

static void remove_param(request* r, void* data, const char* p)
{
    request* config = (request*)data;
    const char* ignore;
    int i = 0;

    while ((ignore = get_value(config, p, i++)) != NULL) {
        unset_value(r, ignore);
    }
}

static void print_ignored_keys(FILE* f, request* data)
{
    const char* ignore = NULL;
    int i              = 0;
    while ((ignore = get_value(data, "ignore", i)) != NULL) {
        if (i == 0) {
            fprintf(f, "%s: Ignoring key(s): %s", tool_name, ignore);
        }
        else {
            fprintf(f, ", %s", ignore);
        }
        ++i;
    }
    if (i > 0)
        fprintf(f, "\n");
}

#define NO_TABLE -1
#define NO_PARAM 0

static void paramtable(const char* p, long* param, long* table, bool paramIdMode)
{
    const char* q = p;
    int len       = strlen(p);

    *param = atol(p);

    while (p && (*p != '.') && ((p - q) < len))
        ++p;

    if ((*p == '.'))
        *table = atol(++p);

    /* This version is grib_api... It should rely on what grib_api returns,
       either param.table or paramId
     */
    if (paramIdMode) {
        /* Special case for param=228015 => 15.228 */
        if ((*param != NO_PARAM) && (*table == NO_TABLE) && (len == 6)) {
            char tbl[4];
            char par[4];
            p = q;
            strncpy(tbl, p, 3);
            tbl[3] = '\0';
            strncpy(par, p + 3, 3);
            par[3] = '\0';
            *param = atol(par);
            *table = atol(tbl);
        }
    }
}

static void find_nc_attributes(const request* subset_r, const request* user_r, ncatt_t* att, const request* config_r, const request* data_r)
{
    const char* split         = NULL;
    int j                     = 0;
    bool set_param_as_name = true;
    long datatable            = 0; /* = atol(get_value(data_r,"_CODETABLE2",0)); */

    if (count_values(user_r, "split") == 0)
        strcat(att->name, "data");

    while ((split = get_value(user_r, "split", j++)) != NULL) {
        if (strcmp(split, "param") != 0) {
            if (count_values(data_r, split) > 1)
                set_param_as_name = false;
        }
    }

    j = 0;
    while ((split = get_value(user_r, "split", j++)) != NULL) {
        bool found    = false;
        request* cfg     = (request*)config_r;
        bool is_param = strcmp(split, "param") == 0;
        /* Only use this parameter in the name if there is more
           than one value in the original request or if param */
        bool setname = ((count_values(data_r, split) > 1) || (is_param && set_param_as_name));

        while (cfg) {
            const char* cfgname = get_value(cfg, "NAME", 0);
            const char* cfgval  = get_value(cfg, "VALUE", 0);
            const char* dataval = NULL;
            int i               = 0;

            if (strcmp(split, cfgname) == 0) {
                while ((dataval = get_value(subset_r, cfgname, i++)) != NULL) {
                    const char* tablestr = get_value(cfg, "TABLE2", 0);
                    long cfgtable        = (is_param && tablestr != NULL) ? atol(get_value(cfg, "TABLE2", 0)) : -1;

                    long cfgparam  = atol(cfgval);
                    long dataparam = atol(dataval);
                    paramtable(dataval, &dataparam, &datatable, false);

                    /* If it is not param and they're EXACTLY equal or
                       being param, they're the same parameter and table */
                    if ((!is_param && (strcmp(dataval, cfgval) == 0)) || (is_param && (dataparam == cfgparam) && (datatable == cfgtable || (datatable == 0 && (cfgtable == 128))))) {
                        const char* val      = NULL;
                        const char* metafile = NULL;

                        if ((val = get_value(cfg, "accuracy", 0)) != NULL)
                            att->nctype = translate_nctype(val);
                        att->long_name  = grib_context_strdup(ctx, (get_value(cfg, "LONG_NAME", 0)));
                        att->short_name = grib_context_strdup(ctx, (get_value(cfg, "SHORT_NAME", 0)));
                        att->units      = grib_context_strdup(ctx, (get_value(cfg, "UNITS", 0)));

                        /* Check if there is more metadata for this variable */
                        if ((metafile = get_value(cfg, "METADATA", 0)) != NULL) {
                            static const char* metadata_dir = NULL;
                            char metapath[1024];

                            if (!metadata_dir)
                                metadata_dir = getenv("METADATA_DIR");

                            snprintf(metapath, sizeof(metapath), "%s/%s", metadata_dir ? metadata_dir : ".", metafile);
                            att->metadata = 0; /* read_request_file(metapath); */
                        }

                        if (setname) {
                            const char* pname = get_value(cfg, "DATA", 0);
                            if (strlen(att->name))
                                strcat(att->name, "_");
                            strcat(att->name, pname);
                        }
                        found = true;

                        grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: Var. name found: '%s'", att->name);
                    }
                }
            }
            cfg = cfg->next;
        }

        /* We have not found configuration for this attribute */
        if (!found) {
            const char* val = get_value(subset_r, split, 0);
            if (val && (setname)) {
                if (strlen(att->name))
                    strcat(att->name, "_");
                strcat(att->name, val);
            }
        }
    }

    /* NetCDF does not allow variable names to start with a digit */
    if (!isalpha(att->name[0])) {
        char buf[1048];
        const char* val = get_value(subset_r, "param", 0);
        snprintf(buf, sizeof(buf), "%s_%s", (val ? val : "p"), att->name);
        strcpy(att->name, buf);
    }
}

static request* first;
static request* last;

static void reqcb(const request* r, int count, axis_t* names, char* vals[], void* data)
{
    request* w = clone_one_request(r);
    int i;
    /*request **d = (request**)data;*/

    int* n = (int*)data;

    w->order = (*n)++;
    for (i = 0; i < count; i++) {
        if (vals[i])
            put_value(w, names[i].name, vals[i], false, false, false);
    }

    if (first == NULL)
        first = w;
    else
        last->next = w;
    last = w;
}

static bool chk_152(int count, axis_t* names, char* vals[])
{
    return true;
}

static void loop(const request* r, bool ml, int index, int count, axis_t* strings, char* values[], loopproc callback, void* data)
{
    if (index < count) {
        parameter* p = find_parameter(r, strings[index].name);

        (void)count_values(r, strings[index].name); /* force list expansion */

        if (p) {
            value* v = p->values;

            while (v) {
                values[index] = v->name;
                loop(r, ml, index + 1, count, strings, values, callback, data);
                v = v->next;
            }
        }
        else {
            values[index] = NULL;
            loop(r, ml, index + 1, count, strings, values, callback, data);
        }
    }
    else if (!ml || chk_152(count, strings, values))
        callback(r, count, strings, values, data);
}

static void values_loop(const request* r, int count, axis_t* parnames, loopproc callback, void* data)
{
    char** values = (char**)grib_context_malloc(ctx, sizeof(char*) * count);
    const char* p = get_value(r, "levtype", 0);
    bool ml    = (bool)(p && (strcmp(p, "ml") == 0));

    if (ml) {
        p = get_value(r, "expect", 0);
        if (p && atol(p) != 0) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "EXPECT provided, special treatment of LNSP");
            grib_context_log(ctx, GRIB_LOG_ERROR, "and other single level parameters disabled");
            ml = false;
        }
    }

    loop(r, ml, 0, count, parnames, values, callback, data);

    grib_context_free(ctx, values);
}

static void names_loop(const request* r, loopproc proc, void* data)
{
    values_loop(r, NUMBER(global_axis), global_axis, proc, data);
}

static request* unwind_one_request(const request* r)
{
    int n = 0;
    first = last = NULL;
    names_loop(r, reqcb, &n);
    return first;
}

static int split_fieldset(fieldset* fs, request* data_r, dataset_t** subsets, const request* user_r, const request* config_r)
{
    const char* split  = NULL;
    int count          = 1;
    int i              = 0;
    request* s         = NULL;
    request* u         = NULL;
    dataset_t* filters = NULL;
    nc_type nctype     = translate_nctype(get_value(user_r, "accuracy", 0));

    s = empty_request("filter");
    while ((split = get_value(user_r, "split", i++)) != NULL) {
        int cnt = count_values(data_r, split);
        if (cnt) {
            count *= count_values(data_r, split);
            valcpy(s, data_r, (char*)split, (char*)split);
        }
    }

    u = unwind_one_request(s);

    free_all_requests(s);

    filters = (dataset_t*)calloc(count, sizeof(dataset_t));
    Assert(filters);

    s = u;
    for (i = 0; i < count; ++i) {
        Assert(s);
        filters[i].filter         = new_hypercube_from_mars_request(s);
        filters[i].fset           = new_fieldset(1);
        filters[i].count          = 0;
        filters[i].filter_request = clone_one_request(s);
        filters[i].bitmap         = false;
        /* filters[i].mmeans  = false; */

        s = s->next;
    }

    for (i = 0; i < fs->count; ++i) {
        bool ok = false;
        field* f   = get_field(fs, i, packed_mem);
        request* g = field_to_request(f);
        int j      = 0;

        while (!ok && (j < count)) {
            ok = (cube_order(filters[j].filter, g) != -1);
            if (ok) {
                const char* p;
                set_field(filters[j].fset, f, filters[j].count++);
                filters[j].bitmap |= f->has_bitmap;
                if ((p = get_value(f->r->r, "_units", 0)) != NULL) {
                    filters[j].att.units = grib_context_strdup(ctx, p);
                }
                if ((p = get_value(f->r->r, "_long_name", 0)) != NULL) {
                    filters[j].att.long_name = grib_context_strdup(ctx, p);
                }
                if ((p = get_value(f->r->r, "_cf_name", 0)) != NULL) {
                    filters[j].att.standard_name = grib_context_strdup(ctx, p);
                }
            }
            j++;
        }
        if (!ok) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Internal error. Lost field %d while splitting fieldset", i);
            print_all_requests(g);
            grib_context_log(ctx, GRIB_LOG_ERROR, "count is %d", count);
            grib_context_log(ctx, GRIB_LOG_ERROR, "First cube is:");
            print_hypercube(filters[0].filter);
            exit(1);
        }
        /* f->purge_header = true; */
        release_field(f);
    }

    for (i = 0; i < count; ++i) {
        filters[i].att.nctype = nctype;
        filters[i].scale      = true;
        filters[i].missing    = nc_type_values[nctype].nc_type_missing;
        find_nc_attributes(filters[i].filter_request, user_r, &(filters[i].att), config_r, data_r);
        grib_context_log(ctx, GRIB_LOG_DEBUG, "grib_to_netcdf: filter[%d] found.- Var. name '%s', nctype: %d, found nctype: %d", i, filters[i].att.name, nctype, filters[i].att.nctype);

        if (strlen(filters[i].att.name) == 0) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Variable name not found");
            exit(1);
        }

        /* Check if we need scaling */
        switch (filters[i].att.nctype) {
            case NC_FLOAT:
            case NC_DOUBLE:
                filters[i].scale = false;
                break;
            default:
                filters[i].scale = true;
                break;
        }
    }

    free_all_requests(u);
    *subsets = filters;
    return count;
}

static void free_subsets(dataset_t* subsets, int count)
{
    int i = 0;

    for (i = 0; i < count; ++i) {
        grib_context_free(ctx, subsets[i].att.units);
        grib_context_free(ctx, subsets[i].att.long_name);
        grib_context_free(ctx, subsets[i].att.short_name);
        grib_context_free(ctx, subsets[i].att.standard_name);
        free_all_requests(subsets[i].att.metadata);
        free_hypercube(subsets[i].filter);
        free_fieldset(subsets[i].fset);
        free_all_requests(subsets[i].filter_request);
    }
    grib_context_free(ctx, subsets);
}

/* Return the number of the given month or -1 if it fails to match */
static int convert_month(const char* pMonthString)
{
    if (strcmp(pMonthString, "jan") == 0)
        return 1;
    if (strcmp(pMonthString, "feb") == 0)
        return 2;
    if (strcmp(pMonthString, "mar") == 0)
        return 3;
    if (strcmp(pMonthString, "apr") == 0)
        return 4;
    if (strcmp(pMonthString, "may") == 0)
        return 5;
    if (strcmp(pMonthString, "jun") == 0)
        return 6;
    if (strcmp(pMonthString, "jul") == 0)
        return 7;
    if (strcmp(pMonthString, "aug") == 0)
        return 8;
    if (strcmp(pMonthString, "sep") == 0)
        return 9;
    if (strcmp(pMonthString, "oct") == 0)
        return 10;
    if (strcmp(pMonthString, "nov") == 0)
        return 11;
    if (strcmp(pMonthString, "dec") == 0)
        return 12;
    return -1; /*Failed*/
}

static bool parsedate(const char* name, long* julian, long* second, bool* isjul)
{
    const char* p = name;
    int n;
    int y = 0, m = 0, d = 0, H = 0, M = 0, S = 0;

    *julian = *second = 0;
    *isjul            = false;

    if (p == 0 || *p == 0)
        return false;

    /* Special ERA Interim grib1 date format: jul-21, sep-02 etc
     * See GRIB-416
     */
    if (isalpha(*p)) {
        char month[32];
        int day = 0;
        n = sscanf(p, "%[^-]-%d", month, &day);
        /* Matched two items (month and day) and month is 3 letters */
        if (n == 2 && strlen(month) == 3) {
            y = 1900; /* no year specified */
            m = convert_month(month);
            if (m == -1)
                return false;
            *julian = grib_date_to_julian(y * 10000 + m * 100 + day);
            *second = 0;
            return true;
        }
    }

    /* year */
    p = parse1(p, &y, &n);
    if (n != 2 && n != 4) /* year string must be 2 or 4 characters long: 93 or 1993 */
        return false;
    if (*p++ != '-')
        return false;

    /* month */
    p = parse1(p, &m, &n);
    if (n == 2) {
        /* day */
        if (*p++ != '-')
            return false;
        p = parse1(p, &d, &n);
        if (n != 2)
            return false;
    }
    else if (n == 3) {
        long j = grib_date_to_julian(y * 10000 + 101) + m - 1;
        j      = grib_julian_to_date(j);
        /* julian day */;
        d      = j % 100;
        m      = (j % 10000) / 100;
        *isjul = true;
    }
    else
        return false;

    if (m == 0 || m > 12) {
        return false; /* month out of range */
    }

    while (*p && isspace(*p))
        p++;

    /* hour */
    p = parse1(p, &H, &n);
    if (n != 0) {
        if (n != 2)
            return false;

        /* minute */
        if (*p++ != ':')
            return false;

        p = parse1(p, &M, &n);
        if (n != 2)
            return false;

        if (*p != 0) {
            /* second */
            if (*p++ != ':')
                return false;
            p = parse1(p, &S, &n);
            if (n != 2)
                return false;
        }
    }

    *julian = grib_date_to_julian(y * 10000 + m * 100 + d);
    *second = H * 3600 + M * 60 + S;

    return *p == 0 ? true : false;
}

static bool check_dimension_name(const char* dim)
{
    /* Dimension name must begin with an alphabetic character, followed by zero
     * or more alphanumeric characters including the underscore */
    int i = 0, len = 0;
    if (!dim)
        return false;

    len = strlen(dim);
    if (len == 0)
        return false;
    if (!isalpha(dim[0]))
        return false;

    for (i = 1; i < len; ++i) {
        const char c = dim[i];
        const int ok = isalnum(c) || c == '_';
        if (!ok)
            return false;
    }
    return true;
}

static int get_creation_mode(int option_kind)
{
    /* Return the mode flag for nc_create based   */
    /* on the kind of netCDF user wants to create */
    int creation_mode = NC_CLOBBER;
    switch (option_kind) {
        case NC_FORMAT_CLASSIC:
            printf("%s: Creating classic file format.\n", tool_name);
            break;
        case NC_FORMAT_64BIT:
            creation_mode |= NC_64BIT_OFFSET;
            printf("%s: Creating large (64 bit) file format.\n", tool_name);
            break;
#ifdef NC_NETCDF4
        case NC_FORMAT_NETCDF4:
            creation_mode |= NC_NETCDF4;
            printf("%s: Creating netCDF-4/HDF5 format.\n", tool_name);
            break;
        case NC_FORMAT_NETCDF4_CLASSIC:
            creation_mode |= NC_NETCDF4 | NC_CLASSIC_MODEL;
            printf("%s: Creating netCDF-4 classic model file format.\n", tool_name);
            break;
#else
        case NC_FORMAT_NETCDF4:
        case NC_FORMAT_NETCDF4_CLASSIC:
            grib_context_log(ctx, GRIB_LOG_ERROR, "%s not built with netCDF4, cannot create netCDF-4 files.", tool_name);
            exit(1);
            break;
#endif
        default:
            fprintf(stderr, "Bad value (%d) for -k option\n", option_kind);
            exit(1);
            break;
    }
    return creation_mode;
}
/*=====================================================================*/

grib_option grib_options[] = {
    { "I:", "key1,key2,...", "\n\t\tIgnore keys. Default method,type,stream,refdate,hdate\n", 0, 1, "method,type,stream,refdate,hdate" },
    { "S:", "key1,key2,...", "\n\t\tSplit according to keys. Default param,expver\n", 0, 1, "param,expver" },
    { "R:", "date", "\n\t\tReference date in the format YYYYMMDD. Default value 19000101.\n", 0, 1, "19000101" },
    { "D:", "NC_DATATYPE",
      "\n\t\tType of data. Possible values NC_BYTE, NC_SHORT, NC_INT, NC_FLOAT, NC_DOUBLE."
      "\n\t\tDefault NC_SHORT\n",
      0, 1, "NC_SHORT" },
    { "T", 0, "Don't use time of validity.\n", 0, 1, 0 },
    { "f", 0, 0, 0, 1, 0 },
    { "o:", "output_file", "\n\t\tThe name of the netCDF output file.\n", 1, 1, 0 },
    { "V", 0, 0, 0, 1, 0 },
    { "M", 0, 0, 0, 1, 0 },
    { "k:", "kind",
      "\n\t\tSpecifies the kind of file to be created. Possible values are:"
      "\n\t\t1 -> netCDF classic file format"
      "\n\t\t2 -> netCDF 64 bit classic file format (Default)"
      "\n\t\t3 -> netCDF-4 file format"
      "\n\t\t4 -> netCDF-4 classic model file format\n",
      0, 1, "2" },
    { "d:", "level",
      "\n\t\tDeflate data (compression level). Only for netCDF-4 output format."
      "\n\t\tPossible values [0,9]. Default None."
      "\n\t\tChunking strategy based on GRIB message.\n",
      0, 1, "6" },
    { "s", 0, "Shuffle data before deflation compression.\n", 0, 1, 0 },
    { "u:", "dimension", "\n\t\tSet dimension to be an unlimited dimension.\n", 0, 1, "time" },
    { "h", 0, 0, 0, 1, 0 },
};

int grib_options_count    = sizeof(grib_options) / sizeof(grib_option);
static fieldset* fs       = NULL;
static request* data_r    = NULL;
request* user_r           = NULL;
static int option_kind    = 2; /* By default NetCDF3, 64-bit offset */
static int deflate_option = 0;

/* Table of formats for legal -k values. Inspired by nccopy */
struct KindValue
{
    const char* name;
    int kind;
} legalkinds[] = {
    { "1", NC_FORMAT_CLASSIC },
    { "classic", NC_FORMAT_CLASSIC },

    /* The 64-bit offset kind */
    { "2", NC_FORMAT_64BIT },
    { "64-bit-offset", NC_FORMAT_64BIT },

    /* NetCDF-4 HDF5 format */
    { "3", NC_FORMAT_NETCDF4 },
    { "hdf5", NC_FORMAT_NETCDF4 },
    { "netCDF-4", NC_FORMAT_NETCDF4 },
    { "netCDF4", NC_FORMAT_NETCDF4 },
    { "enhanced", NC_FORMAT_NETCDF4 },

    /* NetCDF-4 HDF5 format, but using only nc3 data model */
    { "4", NC_FORMAT_NETCDF4_CLASSIC },
    { "hdf5-nc3", NC_FORMAT_NETCDF4_CLASSIC },
    { "netCDF-4 classic model", NC_FORMAT_NETCDF4_CLASSIC },
    { "netCDF4_classic", NC_FORMAT_NETCDF4_CLASSIC },
    { "enhanced-nc3", NC_FORMAT_NETCDF4_CLASSIC },

    /* null terminate*/
    { NULL, 0 }
};

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

    /* GRIB-413: Collect all program arguments into a string */
    const size_t maxLen = sizeof(argvString);
    size_t currLen      = 0;
    for (i = 0; i < argc; ++i) {
        currLen += strlen(argv[i]);
        if (currLen >= maxLen - 1)
            break;
        strcat(argvString, argv[i]);
        if (i != argc - 1)
            strcat(argvString, " ");
    }
    return grib_tool(argc, argv);
}

int grib_tool_before_getopt(grib_runtime_options* options)
{
    return 0;
}

int grib_tool_init(grib_runtime_options* options)
{
    char* p            = NULL;
    char* list         = NULL;
    ctx                = grib_context_get_default();
    options->onlyfiles = 1;
    fs                 = new_fieldset(0);
    data_r             = empty_request(0);
    user_r             = empty_request(0);

    printf("%s: Version ", tool_name);
    grib_print_api_version(stdout);
    printf("\n");

    if (grib_options_on("D:")) {
        set_value(user_r, "accuracy", grib_options_get_option("D:"));
    }
    else {
        set_value(user_r, "accuracy", "NC_SHORT");
    }

    /* Option -S: Split according to keys */
    if (grib_options_on("S:")) {
        char* lasts         = NULL;
        list = grib_options_get_option("S:");
        p    = strtok_r(list, ",", &lasts);
        set_value(user_r, "split", p);
        p = strtok_r(NULL, ",", &lasts);
        while (p != NULL) {
            add_value(user_r, "split", p);
            p = strtok_r(NULL, ",", &lasts);
        }
    }
    else {
        set_value(user_r, "split", "param");
        add_value(user_r, "split", "expver");
    }

    /* Option -I: Ignore keys */
    if (grib_options_on("I:")) {
        char* lasts = NULL;
        list = grib_options_get_option("I:");
        p    = strtok_r(list, ",", &lasts);
        set_value(user_r, "ignore", p);
        p = strtok_r(NULL, ",", &lasts);
        while (p != NULL) {
            add_value(user_r, "ignore", p);
            p = strtok_r(NULL, ",", &lasts);
        }
    }
    else {
        set_value(user_r, "ignore", "method");
        add_value(user_r, "ignore", "type");
        add_value(user_r, "ignore", "stream");
        add_value(user_r, "ignore", "refdate");
        add_value(user_r, "ignore", "hdate");
    }

    if (grib_options_on("T"))
        set_value(user_r, "usevalidtime", "false");
    else
        set_value(user_r, "usevalidtime", "true");

    if (grib_options_on("k:")) {
        struct KindValue* kvalue = NULL;
        char* kind_name          = grib_options_get_option("k:");
        for (kvalue = legalkinds; kvalue->name; kvalue++) {
            if (strcmp(kind_name, kvalue->name) == 0) {
                option_kind = kvalue->kind; /* Found the right kind */
                break;
            }
        }
        if (kvalue->name == NULL) {
            fprintf(stderr, "Invalid value for -k option: %s\n", kind_name);
            fprintf(stderr, "Please use one of:\n");
            for (kvalue = legalkinds; kvalue->name; kvalue++) {
                if (is_number(kvalue->name))
                    fprintf(stderr, "\t%s\n", kvalue->name);
            }
            usage_and_exit();
        }
    }

    if (grib_options_on("d:")) {
        if (option_kind == 3 || option_kind == 4) { /* netCDF-4 */
            char* theArg = grib_options_get_option("d:");
            if (!is_number(theArg) || atol(theArg) < 0 || atol(theArg) > 9) {
                fprintf(stderr, "Invalid deflate option: %s (must be 0 to 9)\n", theArg);
                usage_and_exit();
            }
            set_value(user_r, "deflate", theArg);
            deflate_option = 1;
        }
        else {
            fprintf(stderr, "Invalid deflate option for non netCDF-4 output formats\n");
            usage_and_exit();
        }
    }
    else {
        set_value(user_r, "deflate", "none");
    }

    if (grib_options_on("s")) {
        if (deflate_option)
            set_value(user_r, "shuffle", "true");
        else {
            fprintf(stderr, "Invalid shuffle option. Deflate option needed.\n");
            usage_and_exit();
        }
    }
    else
        set_value(user_r, "shuffle", "false");

    if (grib_options_on("R:")) {
        char* theArg = grib_options_get_option("R:");
        if (!is_number(theArg)) {
            fprintf(stderr, "Invalid reference date: %s\n", theArg);
            usage_and_exit();
        }
        set_value(user_r, "referencedate", theArg);
    }
    else {
        set_value(user_r, "referencedate", "19000101");
    }

    if (grib_options_on("u:")) {
        char* theDimension = grib_options_get_option("u:");
        if (!check_dimension_name(theDimension)) {
            fprintf(stderr, "Invalid dimension: \"%s\"\n", theDimension);
            exit(1);
        }
        set_value(user_r, "unlimited", theDimension);
    }

    get_nc_options(user_r);

    return 0;
}

int grib_tool_new_filename_action(grib_runtime_options* options, const char* filename)
{
    char buf[1024] = {0,};
    int e           = 0;
    int i           = 0;
    grib_handle* h  = NULL;
    grib_file* file = NULL;

    printf("%s: Processing input file '%s'.\n", tool_name, filename);

    file = grib_file_open(filename, "r", &e);
    if (!file || !file->handle)
        return e;

    fseeko(file->handle, 0, SEEK_SET);

    files++;

    while ((h = grib_handle_new_from_file(ctx, file->handle, &e)) != NULL) {
        long length;
        field* g;
        request* r;

        /* process only GRIB for the moment*/
        size_t size = sizeof(buf);
        Assert(grib_get_string(h, "identifier", buf, &size) == 0);
        if (strcmp(buf, "GRIB")) {
            grib_handle_delete(h);
            continue;
        }
        Assert(grib_get_long(h, "totalLength", &length) == 0);

        g = read_field(file, h->offset, length);

        r = empty_request("");
        if (handle_to_request(r, h) != GRIB_SUCCESS) {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Failed to convert GRIB handle to a request");
        }

        /* Keep full MARS description */
        /* copy = clone_one_request(r); */
        /* reqmerge(setup.mars_description,copy); */
        /* reqmerge(setup.mars_description,r); */

        if (i == 1) {
            const char* mmeans = get_value(r, "_MONTHLY_MEANS", 0);
            setup.mmeans       = mmeans ? (atol(mmeans) == 1) : false;
        }
        fcmonth2nbmonths(r);

        if (!setup.refdate) {
            if (setup.auto_refdate)
                setup.refdate = atol(get_value(r, "date", 0));
            else {
                const char* p = get_value(user_r, "referencedate", 0);
                if (is_number(p))
                    setup.refdate = atol(p);
                else {
                    long julian = 0, second = 0;
                    bool isjul;
                    parsedate(p, &julian, &second, &isjul);
                    setup.refdate = grib_julian_to_date(julian);
                }
            }
        }

        validation_time(r);
        if (setup.usevalidtime) {
            unset_value(r, "date");
            unset_value(r, "time");
            unset_value(r, "step");
            unset_value(r, "fcmonth");

            set_value(r, "time", "%.2lf", atof(get_value(r, "_validationtime", 0)));
        }

        g->r = new_field_request(r);

        set_field(fs, g, i++);

        reqmerge(data_r, r);
        free_all_requests(r);

        grib_handle_delete(h);
    }

    if (e != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "%s (message %d)", grib_get_error_message(e), i);
    }

    int e2 = 0;
    grib_file_close(file->name, 0, &e2);
    if (e2 != GRIB_SUCCESS) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Failed to close file %s (%s)", file->name, grib_get_error_message(e2));
    }

    {
        /* Now do some checks */
        request* temp_data_r = fieldset_to_request(fs);
        if (setup.checkvalidtime) {
            int cnt = request_fields(temp_data_r);
            if (fs->count != i || (cnt < i)) {
                grib_context_log(ctx, GRIB_LOG_ERROR, "Wrong number of fields");
                grib_context_log(ctx, GRIB_LOG_ERROR, "File contains %d GRIBs, %d left in internal description, %d in request", i, fs->count, cnt);
                grib_context_log(ctx, GRIB_LOG_ERROR, "The fields are not considered distinct!\n");
                /*grib_context_log(ctx, GRIB_LOG_ERROR, "MARS description");*/
                /*print_all_requests(setup.mars_description);*/
                if (ctx->debug) {
                    grib_context_log(ctx, GRIB_LOG_ERROR, "Internal description");
                    print_all_requests(temp_data_r);
                }
                if (grib_options_on("T")) {
                    grib_context_log(ctx, GRIB_LOG_ERROR, "Hint: This may be due to several fields having the same date, time and step.");
                }
                else {
                    grib_context_log(ctx, GRIB_LOG_ERROR, "Hint: This may be due to several fields having the same validity time.");
                    grib_context_log(ctx, GRIB_LOG_ERROR, "Try using the -T option (Do not use time of validity)");
                }
                exit(1);
            }
        }
        free_all_requests(temp_data_r);
    }

    return e;
}

int grib_tool_new_file_action(grib_runtime_options* options, grib_tools_file* file)
{
    return 0;
}

int grib_tool_new_handle_action(grib_runtime_options* options, grib_handle* h)
{
    return 0;
}

int grib_tool_skip_handle(grib_runtime_options* options, grib_handle* h)
{
    return 0;
}

int grib_tool_finalise_action(grib_runtime_options* options)
{
    request* config_r  = NULL;
    hypercube* dims    = NULL;
    dataset_t* subsets = NULL;
    int count;
    int i;
    int ncid;
    int stat;
    int err           = 0;
    int creation_mode = NC_CLOBBER;

    if (options->outfile == NULL || options->outfile->name == NULL) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "No output file. Exiting!");
        usage_and_exit();
    }

    if (fs->count == 0) {
        grib_context_log(ctx, GRIB_LOG_ERROR, "Input does not contain any field. Exiting!");
        return -1;
    }

    printf("%s: Found %d GRIB field%s in %d file%s.\n", tool_name, fs->count, fs->count > 1 ? "s" : "", files, files > 1 ? "s" : "");

    if (ctx->debug) {
        grib_context_log(ctx, GRIB_LOG_INFO, "Request representing %d fields ", fs->count);
        print_all_requests(data_r);
    }

    /* Split the SOURCE from request into as many datasets as specified */
    count = split_fieldset(fs, data_r, &subsets, user_r, config_r);
    remove_param(data_r, (void*)user_r, "ignore");
    remove_param(data_r, (void*)user_r, "split");
    print_ignored_keys(stdout, user_r);

    dims = new_simple_hypercube_from_mars_request(data_r);
    if (ctx->debug) {
        grib_context_log(ctx, GRIB_LOG_INFO, "Hypercube");
        print_hypercube(dims);
    }

    /* In case there is only 1 DATE+TIME+STEP, set at least 1 time as axis */
    set_always_a_time(dims, data_r);

    /* Create netCDF file */

    printf("%s: Creating netCDF file '%s'\n", tool_name, options->outfile->name);
    printf("%s: NetCDF library version: %s\n", tool_name, nc_inq_libvers());

    creation_mode = get_creation_mode(option_kind);
    stat          = nc_create(options->outfile->name, creation_mode, &ncid);
    if (stat != NC_NOERR) {
        char msg[1024];
        snprintf(msg, sizeof(msg), "nc_create: '%s'", options->outfile->name);
        check_err(msg, stat, __LINE__);
    }

    /* Define netCDF dataset */
    err = define_netcdf_dimensions(dims, fs, ncid, subsets, count, data_r);
    if (err != GRIB_SUCCESS) {
        stat = nc_close(ncid);
        check_err("nc_close", stat, __LINE__);
        stat = nc_delete(options->outfile->name);
        check_err("nc_delete", stat, __LINE__);

        exit(1);
    }

    /* End definitions */
    stat = nc_enddef(ncid);
    check_err("nc_enddef", stat, __LINE__);

    /* Fill dimensions */
    fill_netcdf_dimensions(dims, fs, ncid);

    /* Put data values */
    for (i = 0; i < count; ++i) {
        if (subsets[i].fset) {
            char dataset[100];
            snprintf(dataset, sizeof(dataset), subsets[i].att.name, i + 1);
            put_data(dims, ncid, dataset, &subsets[i]);
        }
        else {
            grib_context_log(ctx, GRIB_LOG_ERROR, "Fieldset %d is empty!!", i + 1);
        }
    }

    stat = nc_close(ncid);
    check_err("nc_close", stat, __LINE__);

    free_all_requests(data_r);
    free_hypercube(dims);
    free_fieldset(fs);
    free_subsets(subsets, count);

    free_nc_options();
    printf("%s: Done.\n", tool_name);

    return 0;
}

int grib_no_handle_action(grib_runtime_options* options, int err)
{
    fprintf(dump_file, "\t\t\"ERROR: unreadable message\"\n");
    return 0;
}


#else
#include <stdio.h>
int main(int argc, char** argv)
{
    printf("\n");
    printf("grib_to_netcdf:\n");
    printf("\n");
    printf(" ecCodes was not compiled with NETCDF enabled\n");
    printf("\n");
    exit(1);
}
#endif