/********************************************************************
 * METEX SOAP service for NCEP reanalysis-1 dataset
 ********************************************************************/

load("soap.dll", "matrix.dll", "netcdf.dll");

SP = soap();

// NCEP netCDF objects;
T = T0 = H = P0 = U = U0 = V = V0 = W = W0 = null;

// vertical wind
WND = matrix("float");

// year of the current time
year = 0;

// path to NCEP data
path = "m:\\ncep1\\";

// these matrices must be global
dims = matrix("int", 4);
grid = matrix("float");

// start service
while(1) {
    try {
        SP.serve();
    }
    catch (e) {
        csv(e);
    }
}


///////////////// SOAP functions ////////////////////////////////////

function get_DIMS()
{
    global:dims[0] = 6;    // hours between datasets
    global:dims[1] = 17;   // levels
    global:dims[2] = 73;   // rows (lat grid)
    global:dims[3] = 144;  // columns (lon grid)
    [ptr, n, e] = dims.ptr();
    return [ptr, n*e];
}

function get_LEV()
{
    d = matrix(1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10);
    grid.clone(d);
    [ptr, n, e] = grid.ptr();
    WND.resize(n, 1);
    WND.fill(0, 0);
    return [ptr, n*e];
}

function get_LAT()
{
    grid.resize(73, 1);
    grid.fill(90.0, -2.5);
    [ptr, n, e] = grid.ptr();
    return [ptr, n*e];
}

function get_LON()
{
    grid.resize(144, 1);
    grid.fill(0.0, 2.5);
    [ptr, n, e] = grid.ptr();
    return [ptr, n*e];
}

function get_Z(yy, mm, dd, hh, i, j)
{
    if (isnull(H)) nc_data(path, yy);
    ret = H.readIJ(yy, mm, dd, hh, i, j);
    [ptr, n, e] = ret.ptr();
    return [ptr, n*e];
}

function get_P0(yy, mm, dd, hh, i, j)
{
    if (isnull(P0)) nc_data(path, yy);
    // Pa to hPa
    return 0.01*P0.readIJ(yy, mm, dd, hh, i, j);
}

function get_U(yy, mm, dd, hh, i, j)
{
    if (isnull(U)) nc_data(path, yy);
    ret = U.readIJ(yy, mm, dd, hh, i, j);
    [ptr, n, e] = ret.ptr();
    return [ptr, n*e];
}

function get_U0(yy, mm, dd, hh, i, j)
{
    if (isnull(U0)) nc_data(path, yy);
    return U0.readIJ(yy, mm, dd, hh, i, j);
}

function get_V(yy, mm, dd, hh, i, j)
{
    if (isnull(V)) nc_data(path, yy);
    ret = V.readIJ(yy, mm, dd, hh, i, j);
    [ptr, n, e] = ret.ptr();
    return [ptr, n*e];
}

function get_V0(yy, mm, dd, hh, i, j)
{
    if (isnull(V0)) nc_data(path, yy);
    return V0.readIJ(yy, mm, dd, hh, i, j);
}

function get_W(yy, mm, dd, hh, i, j)
{
    if (isnull(W)) nc_data(path, yy);
    ret = W.readIJ(yy, mm, dd, hh, i, j);
    // Pa/s to hPa/s
    ret *= 0.01;
    [ptr, n, e] = ret.ptr();
    for (i = 0; i < n; i++) WND[i] = ret[i];
    [ptr, n, e] = WND.ptr();
    return [ptr, n*e];
}

function get_W0(yy, mm, dd, hh, i, j)
{
    if (isnull(W0)) nc_data(path, yy);
    return 0.01*W0.readIJ(yy, mm, dd, hh, i, j);
}

function get_T(yy, mm, dd, hh, i, j)
{
    if (isnull(T)) nc_data(path, yy);
    ret = T.readIJ(yy, mm, dd, hh, i, j);
    [ptr, n, e] = ret.ptr();
    return [ptr, n*e];
}

function get_T0(yy, mm, dd, hh, i, j)
{
    if (isnull(T0)) nc_data(path, yy);
    return T0.readIJ(yy, mm, dd, hh, i, j);
}


function SurfaceU(yy, mm, dd, hh)
{
    if (isnull(U0)) nc_data(path, yy);
    return U0.readSurface(yy, mm, dd, hh);
}

function SurfaceV(yy, mm, dd, hh)
{
    if (isnull(V0)) nc_data(path, yy);
    return V0.readSurface(yy, mm, dd, hh);
}

/////////////////////////////////////////////////////////////////////

function nc_data(path, year)
{
    global:year = year;

    if (isnull(T)) global:T = new ncCDF;
    T.open(path + "air." + year + ".nc", "air");

    if (isnull(T0)) global:T0 = new ncCDF;
    T0.open(path + "air.sig995." + year + ".nc", "air");

    if (isnull(H)) global:H = new ncCDF;
    H.open(path + "hgt." + year + ".nc", "hgt");

    if (isnull(P0)) global:P0 = new ncCDF;
    P0.open(path + "pres.sfc." + year + ".nc", "pres");

    if (isnull(U)) global:U = new ncCDF;
    U.open(path + "uwnd." + year + ".nc", "uwnd");

    if (isnull(U0)) global:U0 = new ncCDF;
    U0.open(path + "uwnd.sig995." + year + ".nc", "uwnd");

    if (isnull(V)) global:V = new ncCDF;
    V.open(path + "vwnd." + year + ".nc", "vwnd");

    if (isnull(V0)) global:V0 = new ncCDF;
    V0.open(path + "vwnd.sig995." + year + ".nc", "vwnd");

    if (isnull(W)) global:W = new ncCDF;
    W.open(path + "omega." + year + ".nc", "omega");

    if (isnull(W0)) global:W0 = new ncCDF;
    W0.open(path + "omega.sig995." + year + ".nc", "omega");
}

class ncCDF
{
    HR = 6;         // time resolution of dataset

    nc = nlev = scale = offset = null;
    
    short = matrix("short");
    float = matrix("float");
    
    function open(fname, vname)
    {
        global:nc = netcdf(fname);


        if (nc.variable("level")) {
            // prepare matrix for reading column data
            [global:nlev] = nc.dims();
            short.resize(nlev, 1);
        }
        else {
            global:nlev = 1;
        }

        nc.variable(vname);
        [n, e, name] = nc.size();        
        if (name != "short") error("unexpected netcdf variable type");
        global:scale = nc.scale_factor;
        global:offset = nc.add_offset;
    }

    function readIJ(yy, mm, dd, hh, i, j)
    {
        if (yy != year) nc_data(path, yy);

        hh /= HR;
        k = integer(cal2jul(yy, mm, dd) - cal2jul(yy, 1, 1));
        k = k*24/HR + hh;

        if (nlev == 1) {
            value = scale*nc[k,i,j] + offset;
            return value;
        }

        short.import(nc[k,*,i,j]);
        float.clone(short);
        float *= scale;
        float += offset;
        return float;
    }

    function readSurface(yy, mm, dd, hh)
    {
        if (yy != year) {
            csv(yy, mm, dd, hh, U0, V0);
            nc_data(path, yy);
        }

        hh /= HR;
        k = integer(cal2jul(yy, mm, dd) - cal2jul(yy, 1, 1));
        k = k*24/HR + hh;

        [m, n] = short.size();
        if (m == 0) {
            // prepare matrix for reading layer data
            [m, m, n] = nc.dims();
            short.resize(m, n);
        }

        short.import(nc[k,*,*]);
        float.clone(short);
        float *= scale;
        float += offset;
    
        [ptr, n, e] = float.ptr();

        return [ptr, n*e];
    }
}