/******************************************************************************

Description:

  Integer function read_GEOS.f
  Reads desired variables from input file (in calling argument list).

Input arguments:
  geos_file             input GEOS file name
  geos_time_index       indicates earlier (=1) or later (=2) GEOS file time

Output arguments:
  none

Function output:
  int return_code       successful completion is zero, otherwise non-zero

  Output variables in ancillary.h

Revision History:
  06/2021     R. Frey   Original version

Calls:
  none

*******************************************************************************/

// Includes

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>
#include <float.h>

//MVCCM includes
//#include "granule_geolocation.h"
//#include "granule.h"
#include "ancillary.h"

//netcdf includes
#include <netcdf.h>

/******************************************************************************/

int read_GEOS(char *geos_file, int geos_time_index)

{

/******************************************************************************/

//  Declarations and initializations.

/******************************************************************************/

    size_t lat_len, lon_len, start[3], edge[3];
    int status, idf, latln_id, lonln_id, nx, ny, nxny, varid;
    int return_code = 0;

/******************************************************************************/

    start[0] = 0;
    start[1] = 0;
    start[2] = 0;
    edge[0] = 0;
    edge[1] = 0;
    edge[2] = 0;

/******************************************************************************/

//  Open file.
//  printf("Opening GEOS file: %s\n", geos_file);
    status = nc_open(geos_file, NC_NOWRITE, &idf);
    if (status != NC_NOERR) {
      printf("Could not open GEOS file: %s %d %d\n", geos_file, idf, status);
      return (-1);
    }

/******************************************************************************/

//  Get data array dimensions.  All meteorologial variables are the same size.
//  First dimension (time) is 1.

    status = nc_inq_dimid(idf, "lat", &latln_id);
    if(status < 0) {
      printf("Could not acquire lat-dimension data: %s\n", geos_file);
      return_code = -2;
      return (return_code);
    }
    status = nc_inq_dimlen(idf, latln_id, &lat_len);

    status = nc_inq_dimid(idf, "lon", &lonln_id);
    if(status < 0) {
      printf("Could not acquire lon-dimension data: %s\n", geos_file);
      return_code = -3;
      return (return_code);
    }
    status = nc_inq_dimlen(idf, lonln_id, &lon_len);

/******************************************************************************/

    nx = lon_len;
    ny = lat_len;
    nxny = lat_len * lon_len;
    nwp_cells = nxny;
    nwp_cols = nx;
    nwp_rows = ny;

//  printf("dims: %d %d %d %d %d %d\n", idf, latln_id, lonln_id, lat_len, lon_len, nxny);
//  printf("time index: %d\n", geos_time_index);

    if(geos_time_index == 1) {

      if ((g_sfct1 = (float *) malloc(sizeof(float) * nxny)) == NULL) {
        printf("Could not allocate memory for g_sfct1\n");
        return(-4);
      }
      if ((g_tpw1 = (float *) malloc(sizeof(float) * nxny)) == NULL) {
        printf("Could not allocate memory for g_tpw1\n");
        return(-5);
      }

    }
    else if (geos_time_index == 2) {

      if ((g_sfct2 = (float *) malloc(sizeof(float) * nxny)) == NULL) {
        printf("Could not allocate memory for g_sfct2\n");
        return(-6);
      }
      if ((g_tpw2 = (float *) malloc(sizeof(float) * nxny)) == NULL) {
        printf("Could not allocate memory for g_tpw2\n");
        return(-7);
      }

    }

/******************************************************************************/

    edge[0] = 1;
    edge[1] = lat_len;
    edge[2] = lon_len;

//  Read skin temperatures.

    status = nc_inq_varid (idf, "TS", &varid);
    if(status != NC_NOERR) {
      printf("Could not find TS variable in %s\n", geos_file);
      return_code = -8;
      return (return_code);
    }

    if(geos_time_index == 1) {
      status = nc_get_vara_float(idf, varid, start, edge, g_sfct1);
    }
    else if (geos_time_index == 2) {
      status = nc_get_vara_float(idf, varid, start, edge, g_sfct2);
    }
    if (status != NC_NOERR) {
      printf("could not get TS data in %s\n", geos_file);
      printf("status: %d\n", status);
      return_code = -9;
      return (return_code);
    }

//  Read total precipitable water.

    status = nc_inq_varid (idf, "TQV", &varid);
    if(status != NC_NOERR) {
      printf("Could not find TQV variable in %s\n", geos_file);
      return_code = -10;
      return (return_code);
    }

    if(geos_time_index == 1) {
      status = nc_get_vara_float(idf, varid, start, edge, g_tpw1);
    }
    else if (geos_time_index == 2) {
      status = nc_get_vara_float(idf, varid, start, edge, g_tpw2);
    }
    if (status != NC_NOERR) {
      printf("could not get TQV data in %s\n", geos_file);
      printf("status: %d\n", status);
      return_code = -11;
      return (return_code);
    }

/******************************************************************************/

//  if(geos_time_index == 1) {
//    printf("GEOS1 data: %f %f %f\n", g_sfct1[0], g_sfct1[9], g_sfct1[nxny-1]);
//    printf("GEOS1 data: %f %f %f\n", g_tpw1[0], g_tpw1[9], g_tpw1[nxny-1]);
//    printf("GEOS1 data: %f %f %f %f %f %f %f %f %f\n",
//           g_sfct1[0], g_lst[0], g_tpw1[0],
//           g_sfct1[32580], g_lst[32580], g_tpw1[32580],
//           g_sfct1[65159], g_lst[65159], g_tpw1[65159]);
//  }
//  if(geos_time_index == 2) {
//    printf("GEOS2 data: %f %f %f\n", g_sfct2[0], g_sfct2[9], g_sfct2[nxny-1]);
//    printf("GEOS2 data: %f %f %f\n", g_tpw2[0], g_tpw2[9], g_tpw2[nxny-1]);
//    printf("GEOS2 data: %f %f %f %f %f %f\n",
//           g_sfct2[0], g_tpw2[0],
//           g_sfct2[32580], g_tpw2[32580],
//           g_sfct2[65159], g_tpw2[65159]);
//  }
//  fflush;

/******************************************************************************/

//  Close file.
    status = nc_close(idf);

/******************************************************************************/

    return (return_code);

/******************************************************************************/

}