From d3252d8f4917f5361c59cda42d2b5b7ae33b5552 Mon Sep 17 00:00:00 2001
From: Paolo Veglio <paolo.veglio@ssec.wisc.edu>
Date: Wed, 15 Jun 2022 19:13:08 +0000
Subject: [PATCH] temp commit 2
---
src/assign_geos_vals.c | 268 ++++++++++++++++
src/bilinearInterpSST.c | 112 +++++++
src/get_GEOS.c | 201 ++++++++++++
src/get_NDVI_background.c | 184 +++++++++++
src/get_Olson_eco.c | 141 +++++++++
src/get_Reynolds_SST.c | 162 ++++++++++
src/get_geos_times.c | 81 +++++
src/get_granule_times.c | 76 +++++
src/get_ti_vars.c | 91 ++++++
src/get_ti_vars_geos.c | 91 ++++++
src/get_ti_weights.c | 69 ++++
src/get_ti_weights_geos.c | 69 ++++
src/getcoord.c | 650 ++++++++++++++++++++++++++++++++++++++
src/read_GEOS.c | 219 +++++++++++++
src/read_GEOS_constants.c | 200 ++++++++++++
src/read_GEOS_lndocn.c | 187 +++++++++++
src/swap_bytes.c | 126 ++++++++
17 files changed, 2927 insertions(+)
create mode 100644 src/assign_geos_vals.c
create mode 100644 src/bilinearInterpSST.c
create mode 100644 src/get_GEOS.c
create mode 100644 src/get_NDVI_background.c
create mode 100644 src/get_Olson_eco.c
create mode 100644 src/get_Reynolds_SST.c
create mode 100644 src/get_geos_times.c
create mode 100644 src/get_granule_times.c
create mode 100644 src/get_ti_vars.c
create mode 100644 src/get_ti_vars_geos.c
create mode 100644 src/get_ti_weights.c
create mode 100644 src/get_ti_weights_geos.c
create mode 100644 src/getcoord.c
create mode 100644 src/read_GEOS.c
create mode 100644 src/read_GEOS_constants.c
create mode 100644 src/read_GEOS_lndocn.c
create mode 100644 src/swap_bytes.c
diff --git a/src/assign_geos_vals.c b/src/assign_geos_vals.c
new file mode 100644
index 0000000..2e4787e
--- /dev/null
+++ b/src/assign_geos_vals.c
@@ -0,0 +1,268 @@
+/*******************************************************************************
+Description:
+
+ Integer function assign_geos_vals.c
+ Assigns GEOS variable values to each pixel location in the input L1b granule.
+ Also spatially interpolates GEOS values to pixel locations if desired.
+ Called from get_GEOS.c
+
+Input arguments:
+ none
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Pointers to inputs are stored in ancillary.h.
+ Pointers to output variables are stored in ancillary.h.
+
+Revision History:
+ 10/2012 R. Frey Original version
+ 01/2021 R. Frey Made changes to accomodate sector processing.
+ Code cleanup.
+ 06/2021 R. Frey Modified from assign_gdas_vals.c
+ 06/2021 R. Frey Added 'snowfr' and 'icefr'.
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+#include "ancillary.h"
+//#include "granule_geolocation.h"
+//#include "granule.h"
+#include "mask_processing_constants.h"
+
+/******************************************************************************/
+
+int assign_geos_vals(float *lat, float *lon)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ extern double fabs(double);
+
+ const int eles=3200;
+ const int lines=3232;
+
+ int row, col;
+ int i, j, ij, i1, i2, j1, j2, ic, jc;
+ int j1c, j2c, i1c, i2c;
+ int line1, line2, ele1, ele2;
+ long int n_pixels = eles * lines;
+ long int nwp_index, nwp_index_c;
+ long int nwp_index_c1, nwp_index_c2, nwp_index_c3, nwp_index_c4;
+ long int indx00, indx01, indx10, indx11;
+ long int idx, idxloc;
+ float lst_sum;
+ double x0 = -180.0;
+ double dx = 0.625;
+ double dxc = 0.3125;
+ double y0 = -90.0;
+ double dy = 0.50;
+ double dyc = 0.25;
+ double x, lat1, lat2, lat00, lat01, lat10, lat11, dyi, pyi, py;
+ double y, lon1, lon2, lon00, lon01, lon10, lon11, dxi, pxi, px;
+ double t, tx, ty, t00, t01, t10, t11;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Allocate memory for assigned data. Structure variables are stored in
+// ancillary.h.
+ if ((grn_anc.sfct = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for sfct\n");
+ return(-1);
+ }
+ if ((grn_anc.tpw = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for tpw\n");
+ return(-2);
+ }
+ if ((grn_anc.snowfr = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for snowfr\n");
+ return(-2);
+ }
+ if ((grn_anc.icefr = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for icefr\n");
+ return(-2);
+ }
+ if ((grn_anc.geos_ocnfr = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for geos_ocnfr\n");
+ return(-2);
+ }
+ if ((grn_anc.landicefr = (float *) malloc(sizeof(float) * n_pixels)) == NULL) {
+ printf("Could not allocate memory for landicefr\n");
+ return(-2);
+ }
+
+/******************************************************************************/
+
+// Get TPW, snow fraction, ice fraction, and bilinearly interpolated surface
+// temperature.
+
+/******************************************************************************/
+
+// Use 0-based line and element numbers.
+ idx = 0;
+ idxloc = 0;
+ ij = 0;
+
+ for ( row=0; row<lines; ++row) {
+ for ( col=0; col<eles; ++col) {
+
+/******************************************************************************/
+
+// Find GEOS 0.5-degree grid box corresponding to current location.
+ y = min(max(lat[idx], -90.00), 90.00);
+ j = (int) ((y - y0 + dy) / dy);
+
+ x = min(max(lon[idx], -180.00), 179.375);
+ i = (int) ((x - x0 + dx) / dx);
+
+ nwp_index = (j-1)*nwp_cols + (i-1);
+ if(nwp_index < 0 || nwp_index > nwp_cells) {
+ printf("NWP grid box out of range: %ld\n", nwp_index);
+ return (-3);
+ }
+
+// Find GEOS 0.25-degree grid box corresponding to current location.
+ jc = (int) ((y - y0 + dyc) / dyc);
+ ic = (int) ((x - x0 + dxc) / dxc);
+ nwp_index_c = ((jc-1)*nwp_cols*2) + (ic-1);
+
+/******************************************************************************/
+// Don't spatially interpolate the following variables.
+ grn_anc.tpw[idxloc] = *(ti_tpw + nwp_index);
+ grn_anc.snowfr[idxloc] = *(g_snow + nwp_index);
+ grn_anc.icefr[idxloc] = *(g_ice + nwp_index);
+ grn_anc.geos_ocnfr[idxloc] = (*(g_frocn + nwp_index_c));
+ grn_anc.landicefr[idxloc] = (*(g_frlandice + nwp_index_c));
+
+// Get non-interpolated value in case interpolation is not performed.
+ grn_anc.sfct[idxloc] = *(ti_sfct + nwp_index);
+
+/******************************************************************************/
+
+// Find four GEOS values that surround the pixel of interest. Begin by
+// finding the two surrounding latitude bands and computing the fractional
+// distance between the pixel and the closest latitude band.
+ lat1 = y0 + ((j-1) * dy);
+ lat2 = y0 + ((j+1) * dy);
+ if( fabs(y-lat1) <= fabs(y-lat2) ) {
+ j1 = j-1;
+ j2 = j;
+ }else {
+ j1 = j;
+ j2 = j+1;
+ }
+ lat00 = y0 + (j1 * dy);
+ lat10 = y0 + (j2 * dy);
+ dyi = lat00 - lat10;
+ pyi = lat00 - y;
+ if(pyi > 0.0) {
+ py = pyi / dyi;
+ }
+ else {
+ py = dy;
+ }
+ if(j1 < 0 || fabs(dyi) < 0.000001 || py > 1.00 || py < 0.00) {
+ printf("problem with j index in spatial interpolation\n");
+ printf("%d %d %ld %d %d %d %f %f %f %f %f %f\n", row,col,idx,j,j1,j2,lat00,lat10,dyi,y,pyi,py);
+ }
+
+// Find surrounding longitudes and compute the fractional distance between
+// the pixel of interest and the closest longitude box.
+ lon1 = x0 + ((i-1) * dx);
+ lon2 = x0 + ((i+1) * dx);
+// if(lon1 < 0.0) lon1 = lon1 + 360.0;
+ if( fabs(x-lon1) < fabs(x-lon2) ) {
+ i1 = i-1;
+ i2 = i;
+ }else {
+ i1 = i;
+ i2 = i+1;
+ }
+
+ lon00 = x0 + (i1 * dx);
+ lon10 = x0 + (i2 * dx);
+ dxi = lon00 - lon10;
+ pxi = lon00 - x;
+ if(pxi > 0.0) {
+ px = pxi/dxi;
+ }
+ else {
+ px = dx;
+ }
+ if(fabs(dxi) < 0.000001 || px > 1.00 || px < 0.00) {
+ printf("problem with i index in spatial interplation\n");
+ printf("%d %d %ld %d %d %d %f %f %f %f %f %f\n", row,col,idx,i,i1,i2,lon00,lon10,dxi,x,pxi,px);
+ }
+
+// Following are indices of surrounding NWP gridboxes.
+ indx00 = (j1 * nwp_cols) + i1;
+ indx01 = (j1 * nwp_cols) + i2;
+ indx10 = (j2 * nwp_cols) + i1;
+ indx11 = (j2 * nwp_cols) + i2;
+
+// Interpolate only if region contains only land or water.
+// Get ocean fraction for the four regions.
+ j1c = (j1*2)-1;
+ j2c = j2*2;
+ i1c = (i1*2)-1;
+ i2c = i2*2;
+ if(j1c < 0) j1c = 1;
+ if(i1c < 0) i1c = 1;
+ nwp_index_c1 = (j1c-1)*nwp_cols*2 + (i1c-1);
+ nwp_index_c2 = (j1c-1)*nwp_cols*2 + (i2c-1);
+ nwp_index_c3 = (j2c-1)*nwp_cols*2 + (i1c-1);
+ nwp_index_c4 = (j2c-1)*nwp_cols*2 + (i2c-1);
+ lst_sum = *(g_frocn+nwp_index_c1) + *(g_frocn+nwp_index_c2) + *(g_frocn+nwp_index_c3) + *(g_frocn+nwp_index_c4);
+// printf("\nlat,lon: %f %f %f %f %f %f\n", grn_geo.lat[idx], grn_geo.lon[idx],lat00,lat10,lon00,lon10);
+// printf("nfrocean1: %d %d %d %d %d %d %d %d\n", j1,j2,i1,i2,j1c,j2c,i1c,i2c);
+// printf("frocean2: %f %f %f %f %f\n", *(g_frocn+nwp_index_c1),*(g_frocn+nwp_index_c2),*(g_frocn+nwp_index_c3),
+// *(g_frocn+nwp_index_c4),lst_sum);
+
+ if(lst_sum == 4.0 || lst_sum == 0.0) {
+
+ t00 = ti_sfct[indx00];
+ t01 = ti_sfct[indx01];
+ tx = ((1.0 - px) * t00) + (px * t01);
+
+ t10 = ti_sfct[indx10];
+ t11 = ti_sfct[indx11];
+ ty = ((1.0 - px) * t10) + (px * t11);
+ t = ((1.0 - py) * tx) + (py * ty);
+
+ grn_anc.sfct[idxloc] = t;
+
+// if(ij < 10)
+// printf("interp: %f %f %f %f %f %f %f %f\n", lst_sum, t00, t01, tx, t10, t11, ty, t);
+// ij++;
+
+ }
+
+ idx++;
+ idxloc++;
+
+ }
+ }
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/bilinearInterpSST.c b/src/bilinearInterpSST.c
new file mode 100644
index 0000000..cea548d
--- /dev/null
+++ b/src/bilinearInterpSST.c
@@ -0,0 +1,112 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "sst.h"
+
+int bilinearInterpSST(int, int, float, float, const SST *, float *);
+
+/* ===== bilinearInterpSST ================================================== */
+/*
+*| Name:
+*| bilinearInterpSST - Bi-linear interpolation scheme used with Reynolds
+*| blended SST data
+*|
+*| Interface:
+*| int
+*| bilinearInterpSST (int xindx, int yindx, float lat, float lon, const SST *sstGrid float returnValue)
+*|
+*| Input:
+*| i Grid index in the x direction
+*| j Grid index in the y direction
+*| lat Pixel latitude value
+*| lon Pixel longitude value
+*| sstGrid The SST grid from the Reynold's SST file
+*|
+*|
+*| Input and Output:
+*| none
+*|
+*| Output:
+*| t Interpolated SST Value
+*|
+*| Return values:
+*| 0 - Success.
+*| -1 - Unsuccessful
+*|
+*|
+*/
+
+int
+ bilinearInterpSST (int i, int j, float lat, float lon, const SST *sstGrid, float * t)
+{
+
+ extern float fabsf(float);
+
+ int npoints_x = 360;
+ int npoints_y = 180;
+ float dx=1.0;
+ float dy=1.0;
+ int doYInterp = 1;
+ int i1, j1;
+ float p, ip0, ip1, jp0, jp1, t00, t01, t0, t10, t11, t1, ldi, ldj;
+
+ if ( j < 0 || j > (npoints_y -1))
+ return(-1);
+ if ( i < 0 || i > (npoints_x -1))
+ return(-1);
+
+ ldi = lon - i;
+ if (ldi >= 0.5) {
+ i1 = i + 1;
+ if ( i == (npoints_x - 1))
+ i1 = 0;
+ }
+ else {
+ i1 = i - 1;
+ if ( i == 0)
+ i1 = npoints_x - 1;
+ }
+
+ ldj = lat - (j - 90.0);
+
+ if (ldj >= 0.5) {
+ j1 = j + 1;
+ if ( j == (npoints_y - 1))
+ doYInterp = 0;
+ }
+ else {
+ j1 = j - 1;
+ if ( j == 0)
+ doYInterp = 0;
+ }
+
+ p = lon - (i1 + 0.5);
+ if ( p < -1.0 )
+ p = p + 360.0;
+ if ( p > 1.0)
+ p = p - 360.0;
+ ip0 = fabsf( p / dx);
+ ip1 = dx - ip0;
+
+ t00 = sstGrid->values[j][i];
+ t01 = sstGrid->values[j][i1];
+ t0 = (ip0 * t00) + (ip1 * t01);
+
+ if( doYInterp ) {
+ p = lat - (j1 - 90.0 + 0.5);
+ jp0 = fabsf( p / dy );
+ jp1 = dy - jp0;
+
+ t10 = sstGrid->values[j1][i];
+ t11 = sstGrid->values[j1][i1];
+ t1 = (ip0 * t10) + (ip1 * t11);
+
+ *t = (jp0 * t0) + (jp1 * t1);
+ }
+ else
+ *t = t0;
+
+ return (0);
+}
+
diff --git a/src/get_GEOS.c b/src/get_GEOS.c
new file mode 100644
index 0000000..a379f1a
--- /dev/null
+++ b/src/get_GEOS.c
@@ -0,0 +1,201 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_GEOS.c
+ Opens and reads selected meteorological variables from two GEOS files.
+ The first file ('geos1') is for a time < 3 hours earlier than the target
+ (granule) time. The second file ('geos2') is for a time < 3 hours later.
+ Variables are temporally interpolated to the target time.
+ Read snow and ice cover values from GEOS land and ocean files, respectively.
+ Surface temperature is then spatially interpolated.
+
+ Called from collect_inputs.c
+
+Input arguments:
+ geos1 file name of earlier GEOS data time
+ geos2 file name of later GEOS data time
+ geos_lnd file name of GEOS land properties file (for snow cover)
+ geos_ocn file name of GEOS ocean properties file (for ice cover)
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Output parameters stored in header file ancillary.h
+
+Revision History:
+ 06/2021 R. Frey Original version
+
+Calls:
+ read_GEOS.c
+ get_granule_times.c
+ get_geos_times.c
+ get_ti_weights.c
+ assign_geos_vals.c
+
+Comments:
+ Variables from the earlier GEOS time have names ending in "1" (e.g., 'sfct1').
+ Variables from the later GEOS time have names ending in "2" (e.g., 'sfct2').
+ These are stored in ancillary.h.
+
+
+*******************************************************************************/
+
+// Includes
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+//#include "granule_geolocation.h"
+//#include "granule.h"
+#include "ancillary.h"
+
+/******************************************************************************/
+
+void get_GEOS(float *lat, float *lon, char *granule_start_time, char *anc_dir, char *geos1, char *geos2, char *geos_lnd,
+ char *geos_ocn, char *geos_cnst, float *tpw, float *snowfr, float *icefr, float *geos_ocnfr,
+ float *landicefr, float *sfct)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ extern int get_granule_times(char[64], float*);
+ extern int get_geos_times(char[64], char[64], int*, int*);
+ extern int get_ti_weights(float, int, int, float*);
+ extern int get_ti_vars(float);
+ extern int assign_geos_vals(float *, float *);
+ extern int read_GEOS(char[256], int);
+ extern int read_GEOS_lndocn(char[256], char[256]);
+ extern int read_GEOS_constants(char[256]);
+
+ char geos1_wpath[256];
+ char geos2_wpath[256];
+ char geosland_wpath[256];
+ char geosocean_wpath[256];
+ char geosconstants_wpath[256];
+ char s1[64], s2[64];
+ char search[2];
+ char *token;
+ int irt;
+ int geoshr1;
+ int geoshr2;
+ int GEOS_time_index, ihr, imin;
+ float targhr;
+ float wt;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+/* Get complete file names. */
+ strcpy(geos1_wpath,anc_dir);strcat(geos1_wpath,"/");strcat(geos1_wpath,geos1);
+ printf("\nGEOS file1: %s\n", geos1_wpath);
+ strcpy(geos2_wpath,anc_dir);strcat(geos2_wpath,"/");strcat(geos2_wpath,geos2);
+ printf("GEOS file2: %s\n", geos2_wpath);
+ strcpy(geosland_wpath,anc_dir);strcat(geosland_wpath,"/");strcat(geosland_wpath,geos_lnd);
+ printf("GEOS LAND file: %s\n", geosland_wpath);
+ strcpy(geosocean_wpath,anc_dir);strcat(geosocean_wpath,"/");strcat(geosocean_wpath,geos_ocn);
+ printf("GEOS OCEAN file: %s\n", geosocean_wpath);
+ strcpy(geosconstants_wpath,anc_dir);strcat(geosconstants_wpath,"/");strcat(geosconstants_wpath,geos_cnst);
+ printf("GEOS CONSTANTS file: %s\n\n", geosconstants_wpath);
+
+/******************************************************************************/
+
+// Read GEOS time-independent (constants) file.
+// Get land/sea tags (fraction ocean) and surface elevation (geopotential ht.).
+ irt = read_GEOS_constants(geosconstants_wpath);
+
+// Read snow and ice fractions from GEOS land and ocean property files.
+// Variables listed in read_GEOS_lndocn.c
+// Pointers to GEOS snow and ice fractions are stored in ancillary.h.
+ irt = read_GEOS_lndocn(geosland_wpath, geosocean_wpath);
+
+// Read selected variables from GEOS 3-hourly analysis files.
+// Variables listed in read_GEOS.c
+// Pointers to GEOS variables are stored in ancillary.h.
+
+ GEOS_time_index = 1;
+ irt = read_GEOS(geos1_wpath, GEOS_time_index);
+ if(irt != 0) {
+ printf("read_GEOS unsuccessful: %s\n", geos1_wpath);
+ // return (-1);
+ }
+
+ GEOS_time_index = 2;
+ irt = read_GEOS(geos2_wpath, GEOS_time_index);
+ if(irt != 0) {
+ printf("read_GEOS unsuccessful: %s\n", geos2_wpath);
+ // return (-2);
+ }
+
+/******************************************************************************/
+
+// Time-interpolate GEOS grids (pointers to data in ancillary.h).
+
+/******************************************************************************/
+
+// Get granule start date, hour, and minute.
+ printf("GEOS interpolation target time from L1b: %s\n", granule_start_time);
+
+ (void) strcpy(search, " ");
+ token = strtok(granule_start_time, search);
+ (void) strcpy(search, ":");
+ token = strtok(NULL, search);
+ ihr = atoi(token);
+ token = strtok(NULL, search);
+ imin = atoi(token);
+ targhr = ihr + (imin / 60.0);
+// printf("hour, min: %d %d %f\n", ihr, imin, targhr);
+
+/******************************************************************************/
+
+// Get GEOS times from file names, at least for now.
+ (void) strcpy(s1, geos1);
+ (void) strcpy(s2, geos2);
+ irt = get_geos_times(s1, s2, &geoshr1, &geoshr2);
+
+/******************************************************************************/
+
+// Calculate time-interpolation weights.
+ irt = get_ti_weights(targhr, geoshr1, geoshr2, &wt);
+ printf("Time-interpolation weights: %f %d %d %f\n", targhr, geoshr1, geoshr2, wt);
+
+/******************************************************************************/
+
+// Perform time interplation. Pointers to inputs and outputs are stored in
+// ancillary.h.
+ irt = get_ti_vars(wt);
+
+/******************************************************************************/
+
+// Assign GEOS variable values to each pixel location in the input L1b granule.
+// Also spatially interpolates GEOS values to pixel locations if desired.
+// Pointers to inputs and outputs are stored in ancillary.h.
+ irt = assign_geos_vals(lat, lon);
+
+/******************************************************************************/
+
+ tpw = grn_anc.tpw;
+ snowfr = grn_anc.snowfr;
+ icefr = grn_anc.icefr;
+ geos_ocnfr = grn_anc.geos_ocnfr;
+ landicefr = grn_anc.landicefr;
+ sfct = grn_anc.sfct;
+
+ int i;
+
+ for (i=0; i<10; i++) {
+ printf("i: %d\ttpw: %f\tsnowfr: %f\ticefr: %f\tgeos_ocnfr: %f\tlandicefr: %f\tsfct: %f\n",
+ i, tpw[i], snowfr[i], icefr[i], geos_ocnfr[i], landicefr[i], sfct[i]);
+ }
+ // return ( return_code);
+
+}
+
diff --git a/src/get_NDVI_background.c b/src/get_NDVI_background.c
new file mode 100644
index 0000000..faa9109
--- /dev/null
+++ b/src/get_NDVI_background.c
@@ -0,0 +1,184 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_NDVI_background.c
+ Opens NDVI background file and reads the data.
+ Called from collect_inputs.c.
+
+Input arguments:
+ Ancillary data directory
+ Input NDVI file name
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Output parameters in stored in ancillary.h.
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//HDF includes
+#include <hdfi.h>
+#include <hdf.h>
+#include <mfhdf.h>
+#include <HdfEosDef.h>
+
+//MVCCM includes
+//#include "granule_geolocation.h"
+//#include "granule.h"
+//#include "ancillary.h"
+
+/******************************************************************************/
+
+void get_NDVI_background(float *lat, float *lon, char *anc_dir, char *NDVI_file, float *g_ndvi_background)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+/******************************************************************************/
+
+ char NDVI_file_wpath[256];
+ char NDVI_data[23][80];
+ char SDS_name[10];
+ char * ndvi_sub="/NDVI/";
+ float map_res, wm_lon, nm_lat;
+ //float *g_ndvi_background;
+ long int nncols = 21600;
+ long int nnrows = 10800;
+ long int xindex, yindex, indx, idx;
+ int32 start[2], stride[2], edge[2], sds_index, sds_id, irt, id;
+ int jd, xerr, yerr, row, col, numberRows;
+ int ii, jj;
+ int msgcounta = 0, msgcounta_lim = 100;
+ int msgcountb = 0, msgcountb_lim = 100;
+ int16 *ndvi_arr;
+ int fnday[23] = {1, 17, 33, 49, 65, 81, 97, 113, 129, 145, 161, 177, 193, 209, 225, 241, 257, 273, 289, 305, 321, 337, 353};
+
+ const int eles=3200;
+ const int lines=3232;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Construct file name.
+
+ (void) strcpy(NDVI_file_wpath, anc_dir);
+ (void) strcat(NDVI_file_wpath, "/");
+ (void) strcat(NDVI_file_wpath, NDVI_file);
+// printf("NDVI file: %s\n", NDVI_file_wpath);
+
+/******************************************************************************/
+
+// Check filename here.
+
+/******************************************************************************/
+
+// Open file. */
+ id = SDstart(NDVI_file_wpath, DFACC_READ);
+ if (id <= 0) {
+ printf("Could not open NDVI file %s\n",NDVI_file_wpath);
+ //return (-1);
+ }
+
+// Select proper SDS and read the data */
+
+ start[0] = 0;
+ start[1] = 0;
+ stride[0] = 1;
+ stride[1] = 1;
+ edge[1] = nncols;
+ edge[0] = nnrows;
+ (void) strcpy(SDS_name, "NDVI");
+
+ ndvi_arr = (int16 *) malloc(21600 * 10800 * 2);
+ sds_index = SDnametoindex(id, SDS_name);
+ sds_id = SDselect(id, sds_index);
+ irt = SDreaddata(sds_id, start, stride, edge, ndvi_arr);
+ if (irt < 0) {
+ printf("Could not read NDVI file %s\n",NDVI_file_wpath);
+ //return (-2);
+ }
+
+// Close NDVI file. */
+ id = SDend(id);
+
+ map_res = 360.0 / (float)(nncols);
+ wm_lon = (map_res / 2.0) - 180.0;
+ nm_lat = 90.0 - (map_res / 2.0);
+
+// Loop through current data and assign NDVI value for each pixel. */
+
+ //g_ndvi_background = (float *)malloc(eles*lines*sizeof(float));
+
+ idx = 0;
+ numberRows = lines;
+ for ( row=0;row<(numberRows );++row) {
+ for ( col=0;col<eles;++col) {
+
+// Initialize.
+ xerr = 0;
+ yerr = 0;
+ g_ndvi_background[idx] = 32.767;
+
+ xindex = (long int)((lon[idx] - wm_lon) / map_res);
+ yindex = (long int)((nm_lat - lat[idx]) / map_res);
+
+ if(xindex < 0 || xindex> nncols) {
+ xerr = 1;
+ if(msgcounta <= msgcounta_lim){
+ printf("xindex out of bounds on NDVI map! \n");
+ msgcounta += 1;
+ }
+ return_code = 1;
+ }
+ if(yindex < 0 || yindex> nnrows) {
+ yerr = 1;
+ if(msgcountb <= msgcountb_lim){
+ printf("yindex out of bounds on NDVI map! \n");
+ msgcountb += 1;
+ }
+ return_code = 2;
+ }
+
+ if(xerr == 0 && yerr == 0) {
+ indx = (yindex*nncols) + xindex;
+ g_ndvi_background[idx] = *(ndvi_arr + indx) * 0.001;
+
+// if (idx < 10) {
+// printf("ndvi %d %d %ld %ld %ld %f %f %d %f\n", row, col, xindex, yindex, indx, grn_geo.lat[idx], grn_geo.lon[idx], *(ndvi_arr + indx), g_ndvi_background[idx]);
+// }
+
+ }
+
+ idx++;
+
+ }
+ }
+
+ //grn_anc.ndvibk = g_ndvi_background;
+
+ free (ndvi_arr);
+
+/******************************************************************************/
+
+// return (return_code);
+
+}
diff --git a/src/get_Olson_eco.c b/src/get_Olson_eco.c
new file mode 100644
index 0000000..c16949f
--- /dev/null
+++ b/src/get_Olson_eco.c
@@ -0,0 +1,141 @@
+/*******************************************************************************
+Description:
+
+ integer function get_Olson_eco.c
+ Opens and reads Olson ecosystem file.
+ Called from collect_inputs.c.
+
+Input arguments:
+ anc_dir ancillary data directory
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Input file name (static) defined here.
+
+ Output stored in header file ancillary.h.
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+
+/******************************************************************************/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+//#include "granule_geolocation.h"
+//#include "granule.h"
+//#include "ancillary.h"
+
+/******************************************************************************/
+
+int get_Olson_eco(float *lat, float *lon, char *anc_dir, unsigned char *granule_ecotype)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+//
+/******************************************************************************/
+
+ extern long int lround(double);
+ extern void getcoord_(double*, double*, double*, double*);
+
+ char eco_file[256];
+ unsigned char eco_index;
+ //unsigned char *granule_ecotype;
+ double lt, ln, ltt, lnn;
+ double drow;
+ double dcol;
+ long int idx, irow, icol, byte_loc;
+ long int max_row = 17347;
+ long int max_col = 40031;
+ int row, col, numberRows;
+
+ const int eles=3200;
+ const int lines=3232;
+
+ size_t bytes_read;
+ FILE *eco_file_ptr;
+
+// Define static ecosystem file name.
+ char eco_file_name[] = "goge1_2_img.v1";
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+/* Get file name */
+ (void) strcpy(eco_file, anc_dir);
+ (void) strcat(eco_file, "/");
+ (void) strcat(eco_file, eco_file_name);
+// printf("ecosystem file: %s \n", eco_file);
+
+/* Open file */
+ eco_file_ptr = fopen(eco_file, "r");
+ if(eco_file_ptr == NULL) {
+ printf("Cannot open ecosystem file %s \n", eco_file);
+ //return (-1);
+ }
+
+ //granule_ecotype = (unsigned char *)malloc(eles*lines*sizeof(unsigned char));
+ numberRows = lines;
+ idx = 0;
+ for ( row=0; row<(numberRows); ++row) {
+ for ( col=0; col<eles ;++col) {
+
+// Get lat/lon
+ lt = (double) lat[idx];
+ ln = (double) lon[idx];
+ ltt = lt;
+ lnn = ln;
+
+// Get row/col of data grid corresponding to lat/lon.
+ (void) getcoord_(<, &ln, &drow, &dcol);
+
+ irow = lround(drow);
+ if(irow > max_row) irow = max_row;
+ if(irow < 1) irow = 1;
+ icol = lround(dcol);
+ if(icol > max_col) icol = max_col;
+ if(icol < 1) icol = 1;
+
+ byte_loc = ( ((irow-1) * max_col) + icol) - 1;
+
+// Read file
+ fseek(eco_file_ptr, byte_loc, SEEK_SET);
+ bytes_read = fread((char *) &eco_index, sizeof(eco_index), 1, eco_file_ptr);
+
+ granule_ecotype[idx] = eco_index;
+
+// if(idx < 10)
+// printf("Olson ecosystem index: %d %d\n", idx, eco_index);
+
+ ++idx;
+
+ }
+ }
+
+/* Close file */
+ (void) fclose(eco_file_ptr);
+
+ //grn_anc.eco = granule_ecotype;
+
+/******************************************************************************/
+
+ //return (return_code);
+
+}
diff --git a/src/get_Reynolds_SST.c b/src/get_Reynolds_SST.c
new file mode 100644
index 0000000..7cef9e7
--- /dev/null
+++ b/src/get_Reynolds_SST.c
@@ -0,0 +1,162 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_Reynolds_SST
+ Opens and reads SSTs from Reynolds SST files.
+ Called from collect_inputs.c.
+
+Input arguments:
+ ancillary data directory
+ Reynolds SST file name
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Output parameters stored in ancillary.h
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ bilinearInterpSST
+ swapbyte4
+
+*******************************************************************************/
+
+// Includes
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+//#include "granule.h"
+#include "sst.h"
+
+/******************************************************************************/
+
+void get_Reynolds_SST(float *lat, float *lon, char *anc_dir, char *SST_file, float *sstInterp)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+/******************************************************************************/
+
+ extern int bilinearInterpSST(int, int, float, float, const SST *, float *);
+ extern int swapbyte4(void *, int);
+ extern float fmaxf(float, float);
+ extern float fminf(float, float);
+
+ /* number of points around a latitude circle */
+ long int npoints_x = 360;
+ /* number of grid points along a meridian */
+ long int npoints_y = 180;
+ /* the 2 possible sizes of an SST file (formatted/unformatted) */
+ long int unfSize = 583264;
+ char sst_file_wpath[256];
+
+ char fileFormat[4] = "unk";
+ float xlon, ytemp;
+ //float * sstInterp;
+ float tempVal;
+ long int fileSize;
+ long int indx;
+ int i, j, xindx, yindx;
+ int ok;
+
+ const int eles=3200;
+ const int lines=3232;
+
+ FILE * sst_file_ptr;
+ size_t bytes_read;
+ SST sstGrid;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+ strcpy(sst_file_wpath, anc_dir);
+ strcat(sst_file_wpath, "/");
+ strcat(sst_file_wpath, SST_file);
+// printf("SST file: %s \n", sst_file_wpath);
+
+ /* Open file */
+ sst_file_ptr = fopen(sst_file_wpath, "rb");
+ if(sst_file_ptr == NULL) {
+ printf("Cannot open Reynolds SST file %s\n", sst_file_wpath);
+ // return (-1);
+ }
+ /* Get the file size */
+ fseek(sst_file_ptr, 0, SEEK_END);
+ fileSize = ftell(sst_file_ptr);
+ fseek(sst_file_ptr, 0, SEEK_SET);
+
+ /* read if file is unformatted only */
+ if ( fileSize == unfSize ) {
+ strcpy(fileFormat,"unf");
+ bytes_read = fread(&sstGrid, sizeof(sstGrid), 1, sst_file_ptr);
+ (void) swapbyte4(&sstGrid, 12 + npoints_x * npoints_y);
+ }
+ else {
+ /* Close file */
+ (void) fclose(sst_file_ptr);
+ // return (-1);
+ }
+ /* Close file */
+ (void) fclose(sst_file_ptr);
+
+ // to make it work with cython I don't need to allocate memory... I guess it's because cython
+ // already does it?!?
+ //sstInterp = (float *) malloc(eles*lines*sizeof(float));
+
+
+ indx = 0;
+ for ( i=0; i < lines; ++i) {
+ for ( j=0; j < eles; ++j) {
+
+ xindx = 0;
+ yindx = 0;
+ ok = -1;
+ sstInterp[indx] = 0.0;
+
+// if(imgr1->landmask[indx] != 1) {
+
+ if (lon[indx] < 0.0)
+ xlon = lon[indx] + 360.0;
+ else
+ xlon = lon[indx];
+
+ xindx = (int) fminf(fmaxf(xlon,0.0),359.99);
+ ytemp = fminf(fmaxf(lat[indx],-89.99), 89.99);
+ yindx = (int) (ytemp + 90.0);
+
+ ok = bilinearInterpSST(xindx,yindx,lat[indx],xlon,&sstGrid,&tempVal);
+
+ if (ok == 0) {
+ sstInterp[indx] = tempVal + 273.15;
+ }
+
+ if(ok != 0)
+ printf("Bad SST interpolation: %d %d %d %d %d %d %f\n", indx, i, j, xindx, yindx, ok, sstInterp[indx]);
+
+// }
+
+// if(indx < 10)
+// printf("SSTs: %d %d %d %d %f\n", indx, xindx, yindx, ok, sstInterp[indx]);
+
+ indx++;
+ }
+ }
+ // grn_anc.reynSST = sstInterp;
+
+/******************************************************************************/
+
+ //return (return_code);
+
+}
diff --git a/src/get_geos_times.c b/src/get_geos_times.c
new file mode 100644
index 0000000..447fc5b
--- /dev/null
+++ b/src/get_geos_times.c
@@ -0,0 +1,81 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_geos_times.c
+ Parses GEOS hour from filename.
+ Called from get_GEOS.c
+
+Input arguments:
+ s1 earlier GEOS file name
+ s2 later GEOS file name
+
+Output arguments:
+ geoshr1 integer time for earlier GEOS file
+ geoshr2 integer time for later GEOS file
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+/******************************************************************************/
+
+int get_geos_times(char s1[64], char s2[64], int *geoshr1, int *geoshr2)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ char search[4];
+ char *geos_time = "0";
+ char *token;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Get GEOS start hours.
+
+// printf("get time 1: %s\n", s1);
+// printf("get time 2: %s\n", s2);
+
+ (void) strcpy(search, "_");
+ token = strtok(s1, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ (void) strcpy(search, "00.");
+ geos_time = strtok(token, search);
+ *(geoshr1) = atoi(geos_time);
+
+ (void) strcpy(search, "_");
+ token = strtok(s2, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ token = strtok(NULL, search);
+ (void) strcpy(search, "00.");
+ geos_time = strtok(token, search);
+ *(geoshr2) = atoi(geos_time);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/get_granule_times.c b/src/get_granule_times.c
new file mode 100644
index 0000000..ef97765
--- /dev/null
+++ b/src/get_granule_times.c
@@ -0,0 +1,76 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_granule_times.c
+ Parses granule start date, hour, and minute from IFF attribute.
+ Called from get_GDAS.c
+
+Input arguments:
+ s1 granule start time string from IFF attribute
+
+Output arguments:
+ targhr decimal granule start hour
+ (target hour for time interplation of GDAS
+ variables)
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+/******************************************************************************/
+
+int get_granule_times(char s1[64], float *targhr)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ char s2[64];
+ char search[2];
+ char *start_time = "0";
+ char *start_date = "0";
+ char *start_hour = "0";
+ char *start_min = "0";
+ int imin;
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Get granule start date, hour, and minute.
+
+ (void) strcpy(search, " ");
+ start_date = strtok(s1, search);
+ start_time = strtok(NULL, search);
+
+ (void) strcpy(s2, start_time);
+ (void) strcpy(search, ":");
+ start_hour = strtok(s2, search);
+ start_min = strtok(NULL, search);
+ printf("\nstart date, hour and minute: %s %s %s\n", start_date, start_hour, start_min);
+
+// Set decimal target hour for time interpolation.
+ imin = atoi(start_min);
+ *(targhr) = (float)atoi(start_hour) + ((float)imin / 60.0);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/get_ti_vars.c b/src/get_ti_vars.c
new file mode 100644
index 0000000..b48470f
--- /dev/null
+++ b/src/get_ti_vars.c
@@ -0,0 +1,91 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_ti_vars_geos.c
+ Calculates time-interpolated variables from the two GEOS times.
+ Called from get_GEOS.c
+
+Input arguments:
+ wt time interpolation weight for later GEOS time
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Inputs are taken from pointers stored in ancillary.h.
+ Pointers to time-interpolated variables are stored in ancillary.h.
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+#include "ancillary.h"
+
+/******************************************************************************/
+
+int get_ti_vars(float wt)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ long int n;
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Allocate memory for time-interpolated data.
+ if ((ti_sfct = (float *) malloc(sizeof(float) * nwp_cells)) == NULL) {
+ printf("Could not allocate memory for ti_sfct\n");
+ return(-2);
+ }
+ if ((ti_tpw = (float *) malloc(sizeof(float) * nwp_cells)) == NULL) {
+ printf("Could not allocate memory for ti_tpw\n");
+ return(-4);
+ }
+
+/******************************************************************************/
+
+ for (n=0; n<nwp_cells; n++) {
+
+ ti_sfct[n] = (g_sfct1[n] >= 0.0 && g_sfct2[n] >= 0.0) ?
+ (wt*(g_sfct2[n]) + (1.0-wt)*(g_sfct1[n])) : -999.99;
+// if(n >= 1000 && n <= 1009)
+// printf("ti sfc t: %ld %f %f %f %f\n", n, wt, g_sfct1[n], g_sfct2[n], ti_sfct[n]);
+
+ ti_tpw[n] = (g_tpw1[n] >= 0.0 && g_tpw2[n] >= 0.0) ?
+ (wt*(g_tpw2[n]/10.0) + (1.0-wt)*(g_tpw1[n]/10.0)) : -999.99;
+// if(n >= 1000 && n <= 1009)
+// printf("ti tpw: %ld %f %f %f %f\n", n, wt, g_tpw1[n], g_tpw2[n], ti_tpw[n]);
+
+ }
+
+/******************************************************************************/
+
+// Free unneeded memory.
+ free(g_sfct1);
+ free(g_sfct2);
+ free(g_tpw1);
+ free(g_tpw2);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/get_ti_vars_geos.c b/src/get_ti_vars_geos.c
new file mode 100644
index 0000000..b48470f
--- /dev/null
+++ b/src/get_ti_vars_geos.c
@@ -0,0 +1,91 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_ti_vars_geos.c
+ Calculates time-interpolated variables from the two GEOS times.
+ Called from get_GEOS.c
+
+Input arguments:
+ wt time interpolation weight for later GEOS time
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Inputs are taken from pointers stored in ancillary.h.
+ Pointers to time-interpolated variables are stored in ancillary.h.
+
+Revision History:
+ 10/2012 R. Frey Original version
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+//MVCCM includes
+#include "ancillary.h"
+
+/******************************************************************************/
+
+int get_ti_vars(float wt)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ long int n;
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Allocate memory for time-interpolated data.
+ if ((ti_sfct = (float *) malloc(sizeof(float) * nwp_cells)) == NULL) {
+ printf("Could not allocate memory for ti_sfct\n");
+ return(-2);
+ }
+ if ((ti_tpw = (float *) malloc(sizeof(float) * nwp_cells)) == NULL) {
+ printf("Could not allocate memory for ti_tpw\n");
+ return(-4);
+ }
+
+/******************************************************************************/
+
+ for (n=0; n<nwp_cells; n++) {
+
+ ti_sfct[n] = (g_sfct1[n] >= 0.0 && g_sfct2[n] >= 0.0) ?
+ (wt*(g_sfct2[n]) + (1.0-wt)*(g_sfct1[n])) : -999.99;
+// if(n >= 1000 && n <= 1009)
+// printf("ti sfc t: %ld %f %f %f %f\n", n, wt, g_sfct1[n], g_sfct2[n], ti_sfct[n]);
+
+ ti_tpw[n] = (g_tpw1[n] >= 0.0 && g_tpw2[n] >= 0.0) ?
+ (wt*(g_tpw2[n]/10.0) + (1.0-wt)*(g_tpw1[n]/10.0)) : -999.99;
+// if(n >= 1000 && n <= 1009)
+// printf("ti tpw: %ld %f %f %f %f\n", n, wt, g_tpw1[n], g_tpw2[n], ti_tpw[n]);
+
+ }
+
+/******************************************************************************/
+
+// Free unneeded memory.
+ free(g_sfct1);
+ free(g_sfct2);
+ free(g_tpw1);
+ free(g_tpw2);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/get_ti_weights.c b/src/get_ti_weights.c
new file mode 100644
index 0000000..607489d
--- /dev/null
+++ b/src/get_ti_weights.c
@@ -0,0 +1,69 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_ti_weights_geos.c
+ Calculates time-interpolation weights from L1b start time and the two GEOS
+ file times.
+ Computed weight is for the later of two GEOS grids. The weight of the earlier
+ one is simply (1.0 - weight).
+ Called from get_GEOS.c
+
+Input arguments:
+ targhr decimal start time of input L1b data
+ gdashr1 integer hour of earlier GEOS file
+ gdashr2 integer hour of later GEOS file
+
+Output arguments:
+ wt time-interpolation weight for later GEOS time
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+Revision History:
+ 10/2012 R. Frey Original version
+ 06/2021 R. Frey Modified to use GMAO GEOS files.
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+/******************************************************************************/
+
+int get_ti_weights(float targhr, int geoshr1, int geoshr2, float *wt)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Calculate time-interpolation weight.
+
+ if(geoshr2 > geoshr1) *(wt) = (targhr-(float)geoshr1) / ((float)geoshr2-(float)geoshr1);
+ if(geoshr2 == geoshr1) *(wt) = 0.0;
+ if(geoshr2 < geoshr1) {
+ geoshr2 += 24;
+ *(wt) = (targhr-(float)geoshr1) / ((float)geoshr2-(float)geoshr1);
+ if( *(wt) < 0.0) {
+ printf("GEOS times are incorrect! T-I weight: %f %d %d %f\n", targhr, geoshr1, geoshr2, *(wt));
+ return_code = -1;
+ }
+ }
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/get_ti_weights_geos.c b/src/get_ti_weights_geos.c
new file mode 100644
index 0000000..607489d
--- /dev/null
+++ b/src/get_ti_weights_geos.c
@@ -0,0 +1,69 @@
+/*******************************************************************************
+Description:
+
+ Integer function get_ti_weights_geos.c
+ Calculates time-interpolation weights from L1b start time and the two GEOS
+ file times.
+ Computed weight is for the later of two GEOS grids. The weight of the earlier
+ one is simply (1.0 - weight).
+ Called from get_GEOS.c
+
+Input arguments:
+ targhr decimal start time of input L1b data
+ gdashr1 integer hour of earlier GEOS file
+ gdashr2 integer hour of later GEOS file
+
+Output arguments:
+ wt time-interpolation weight for later GEOS time
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+Revision History:
+ 10/2012 R. Frey Original version
+ 06/2021 R. Frey Modified to use GMAO GEOS files.
+
+Calls:
+ none
+
+*******************************************************************************/
+
+// Includes
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+/******************************************************************************/
+
+int get_ti_weights(float targhr, int geoshr1, int geoshr2, float *wt)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+ int return_code = 0;
+
+/******************************************************************************/
+
+// Calculate time-interpolation weight.
+
+ if(geoshr2 > geoshr1) *(wt) = (targhr-(float)geoshr1) / ((float)geoshr2-(float)geoshr1);
+ if(geoshr2 == geoshr1) *(wt) = 0.0;
+ if(geoshr2 < geoshr1) {
+ geoshr2 += 24;
+ *(wt) = (targhr-(float)geoshr1) / ((float)geoshr2-(float)geoshr1);
+ if( *(wt) < 0.0) {
+ printf("GEOS times are incorrect! T-I weight: %f %d %d %f\n", targhr, geoshr1, geoshr2, *(wt));
+ return_code = -1;
+ }
+ }
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/getcoord.c b/src/getcoord.c
new file mode 100644
index 0000000..8b5872b
--- /dev/null
+++ b/src/getcoord.c
@@ -0,0 +1,650 @@
+/*
+!C **********************************************************************
+!Description:
+ GIHLS211 transforms latitude, longitude coordinates to line, sample
+ Renamed GIHLS211 to GETCOORD !!
+ coordinates for an image in the Goode's Interrupted Homolosine
+ projection. This routine was complied and run using the C compiler on
+ SunOS 4.2 (UNIX). Results were accurate at the time of testing, but they
+ are by no means guaranteed!
+
+!Input Parameters:
+lt latitude (double precision)
+ln longitude (double precision)
+
+!Output Parameters:
+lne line number of 1km mask for given latitude/longitude
+smp element number of 1km maks for given latitude/longitude
+
+!Revision History:
+
+!Team-unique Header:
+
+!References and Credits:
+
+1. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
+ U.S. Geological Survey Professional Paper 1453 , United State Government
+ Printing Office, Washington D.C., 1989.
+
+2. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
+ Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
+ State Government Printing Office, Washington D.C., 1987.
+
+3. Goode, J.P., 1925, The Homolosine projection: a new device for
+ portraying the Earth's surface entire: Assoc. Am. Geographers, Annals,
+ v. 15, p. 119-125
+
+4. Steinwand, Daniel R., "Mapping Raster Imagery to the Interrupted
+ Goode Homolosine Projection", 1993. In press--IJRS.
+
+!Design Notes:
+!END************************************************************************/
+#include <stdio.h>
+#include <math.h>
+
+
+/* Modified by JC Guu 08/27/96 */
+/* function prototype is given here for the */
+/* functions in this module. */
+
+extern void getcoord_(double *lt, double *ln, double *lne, double *smp);
+/*void getcoord_(double *lt, double *ln, double *lne, double *smp);*/
+void bad_input_parms();
+void goode_init(double r);
+int goode_forward(double lon, double lat, double* x, double* y);
+/*
+void ptitle(char* A);
+void radius(double A);
+*/
+void giherror(char* what, char* where);
+void getcoord_sincos(double val, double* sin_val, double* cos_val);
+int getcoord_sign(double x);
+double getcoord_adjust_lon(double x);
+double getcoord_e0fn(double x);
+double getcoord_e1fn(double x);
+double getcoord_e2fn(double x);
+double getcoord_mlfn(double e0,double e1,double e2,double phi);
+
+
+#define PI 3.141592653589793238
+#define HALF_PI PI*0.5
+#define TWO_PI PI*2.0
+#define EPSLN 1.0e-10
+#define R2D 57.2957795131
+#define D2R 0.0174532925199
+
+#define OK 1
+#define ERROR -1
+#define IN_BREAK -2
+
+/* Variables common to all subroutines in this code file
+ -----------------------------------------------------*/
+static double R; /* Radius of the earth (sphere) */
+static double lon_center[12]; /* Central meridians, one for each region */
+static double feast[12]; /* False easting, one for each region */
+
+/* Transformation routine
+ -------------------------------*/
+void getcoord_(double *lt, double *ln, double *lne, double *smp)
+{
+float pixsiz;
+double x,y;
+double lat, lon, line, samp;
+/*
+int nl,ns;
+*/
+double ul_x, ul_y;
+
+/* Modified by JC Guu 08/30/96 */
+/* The variable status is added here to */
+/* process the return code of */
+/* goode_forward(). */
+
+int status;
+
+/* Modified by JC Guu 08/26/96 */
+/* The type of constant is specified. */
+
+pixsiz = 1000.0F;
+
+/* Report parameters and image geometric characteristics to the user
+ -----------------------------------------------------------------*/
+/*printf("Converting latitude, longitude to line, sample coordinates\n\n");
+printf("Pixel size is %f km\n\n", pixsiz/1000.0); */
+
+/* Modified by JC Guu 08/26/96 */
+/* The equality comparision */
+/* "pixsiz == 1000.0" is changed to an */
+/* inequality comparision */
+/* "fabs(pixsiz-1000.0F) < 0.000001F". */
+
+
+if (fabs(pixsiz-1000.0F) < 0.000001F) {
+ ul_x = -20015000.0;
+ ul_y = 8673000.0;
+/*
+ nl = 17347;
+ ns = 40031;
+*/
+ }
+else {
+ ul_x = -20011500.0;
+ ul_y = 8669500.0;
+/*
+ nl = 2168;
+ ns = 5004;
+*/
+ }
+/*printf("Image size is %d lines by %d samples with an upper left\n",nl,ns);
+printf("corner of UL_X = %lf and UL_Y = %lf meters.\n", ul_x, ul_y); */
+
+/* Initialize the Interrupted Goode Homolosine projection
+ ------------------------------------------------------*/
+goode_init(6370997.0);
+
+/* Process point
+ ------------------*/
+ lat = *lt;
+ lon = *ln;
+/* printf("%12.6f %12.6f\n",lat, lon);*/
+ lat *= D2R;
+ lon *= D2R;
+
+/* Modified by JC Guu 08/30/96 */
+/* The variable status is used to */
+/* capture the return value form */
+/* goode_forward(). */
+/* Valid return code to process the */
+/* return status could be added in the */
+/* furture. */
+
+ status=goode_forward(lon, lat, &x, &y);
+ if(status==ERROR)
+ /*printf("Newton-Raphson method failed to converge\n");*/
+ ;
+
+ line = (ul_y - y) / pixsiz + 1.0;
+ samp = (x - ul_x) / pixsiz + 1.0;
+/* printf("%12.6f %12.6f\n",lat, lon);
+ printf("%12.6f %12.6f\n",line, samp); */
+
+ *lne = line;
+ *smp = samp;
+ *lt = lat;
+ *ln = lon;
+
+ return;
+}
+
+
+
+/* Modified by JC Guu 08/26/96 */
+/* The return type void is given. */
+
+void bad_input_parms() {
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to report bad parameters.
+!Input Parameters: none
+!Output Parameters: none
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 09/03/96 */
+/* The prohibited printf() code is */
+/* commented out. */
+
+/* printf("Syntax: gihll2ls pixsize\n");
+printf(" pixsize in km = 1 or 8\n\n");*/
+}
+
+
+
+/* Modified by JC Guu 08/26/96 */
+/* The old K&R style function definition */
+/* is changed to the modern form and the */
+/* return type is given as void. */
+
+void goode_init(double r)
+{
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Initialize the Goode`s Homolosine projection.
+ Place parameters in static storage for common use.
+!Input Parameters: r Radius of the earth (sphere)
+!Output Parameters: lon_center[12] Central meridians, one for
+ each region
+ feast[12]; False easting, one for each
+ region
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Place parameters in static storage for common use
+ -------------------------------------------------*/
+
+R = r;
+
+/* Initialize central meridians for each of the 12 regions
+ -------------------------------------------------------*/
+lon_center[0] = -1.74532925199; /* -100.0 degrees */
+lon_center[1] = -1.74532925199; /* -100.0 degrees */
+lon_center[2] = 0.523598775598; /* 30.0 degrees */
+lon_center[3] = 0.523598775598; /* 30.0 degrees */
+lon_center[4] = -2.79252680319; /* -160.0 degrees */
+lon_center[5] = -1.0471975512; /* -60.0 degrees */
+lon_center[6] = -2.79252680319; /* -160.0 degrees */
+lon_center[7] = -1.0471975512; /* -60.0 degrees */
+lon_center[8] = 0.349065850399; /* 20.0 degrees */
+lon_center[9] = 2.44346095279; /* 140.0 degrees */
+lon_center[10] = 0.349065850399; /* 20.0 degrees */
+lon_center[11] = 2.44346095279; /* 140.0 degrees */
+
+/* Initialize false eastings for each of the 12 regions
+ ----------------------------------------------------*/
+feast[0] = R * -1.74532925199;
+feast[1] = R * -1.74532925199;
+feast[2] = R * 0.523598775598;
+feast[3] = R * 0.523598775598;
+feast[4] = R * -2.79252680319;
+feast[5] = R * -1.0471975512;
+feast[6] = R * -2.79252680319;
+feast[7] = R * -1.0471975512;
+feast[8] = R * 0.349065850399;
+feast[9] = R * 2.44346095279;
+feast[10] = R * 0.349065850399;
+feast[11] = R * 2.44346095279;
+
+/* Report parameters to the user
+ -----------------------------*/
+/*ptitle("Goode`s Homolosine Equal Area");
+radius(r): */
+
+/* Modified by JC Guu 08/26/96 */
+/* The function should not return anything */
+/* here. */
+
+}
+
+
+
+/* Modified by JC Guu 08/26/96 */
+/* The old K&R style function definition */
+/* is changed to the modern form and the */
+/* return type is given as void. */
+
+int goode_forward(double lon, double lat, double* x, double* y)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description:Goode`s Homolosine forward equations--mapping
+ lat,long to x,y
+!Input Parameters: lon Longitude
+ lat Latitude
+!Output Parameters: x X projection coordinate
+ y Y projection coordinate
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+{
+
+/* Modified by JC Guu 08/30/96 */
+/* Function prototype is provided for */
+/* getcoord_adjust_lon() therefore the following */
+/* declaration is not needed. */
+
+/* double getcoord_adjust_lon();*/ /* Function to adjust longitude to -180 - 180 */
+double delta_lon; /* Delta longitude (Given longitude - center */
+double theta;
+double delta_theta;
+double constant;
+int i;
+int region;
+
+/* Forward equations
+ -----------------*/
+if (lat >= 0.710987989993) /* if on or above 40 44' 11.8" */
+ {
+ if (lon <= -0.698131700798) region = 0; /* If to the left of -40 */
+ else region = 2;
+ }
+else if (lat >= 0.0) /* Between 0.0 and 40 44' 11.8" */
+ {
+ if (lon <= -0.698131700798) region = 1; /* If to the left of -40 */
+ else region = 3;
+ }
+else if (lat >= -0.710987989993) /* Between 0.0 & -40 44' 11.8" */
+ {
+ if (lon <= -1.74532925199) region = 4; /* If between -180 and -100 */
+ else if (lon <= -0.349065850399) region = 5; /* If between -100 and -20 */
+ else if (lon <= 1.3962634016) region = 8; /* If between -20 and 80 */
+ else region = 9; /* If between 80 and 180 */
+ }
+else /* Below -40 44' */
+ {
+ if (lon <= -1.74532925199) region = 6; /* If between -180 and -100 */
+ else if (lon <= -0.349065850399) region = 7; /* If between -100 and -20 */
+ else if (lon <= 1.3962634016) region = 10; /* If between -20 and 80 */
+ else region = 11; /* If between 80 and 180 */
+ }
+
+if (region==1||region==3||region==4||region==5||region==8||region==9)
+ {
+ delta_lon = getcoord_adjust_lon(lon - lon_center[region]);
+ *x = feast[region] + R * delta_lon * cos(lat);
+ *y = R * lat;
+ }
+else
+ {
+ delta_lon = getcoord_adjust_lon(lon - lon_center[region]);
+ theta = lat;
+ constant = PI * sin(lat);
+
+/* Iterate using the Newton-Raphson method to find theta
+ -----------------------------------------------------*/
+ for (i=0;;i++)
+ {
+ delta_theta = -(theta + sin(theta) - constant) / (1.0 + cos(theta));
+ theta += delta_theta;
+ if (fabs(delta_theta) < EPSLN) break;
+ if (i >= 30)
+ {giherror("Iteration failed to converge","Goode-forward");return(ERROR);}
+ }
+ theta /= 2.0;
+ *x = feast[region] + 0.900316316158 * R * delta_lon * cos(theta);
+
+/* Modified by JC Guu */
+/* The explicit cast is provided. */
+
+ *y = R * (1.4142135623731 * sin(theta) - 0.0528035274542 * ((double) getcoord_sign(lat)));
+ }
+
+
+return(OK);
+}
+
+
+
+/* Modified by JC Guu 08/27/96 */
+/* The function definition style is changed */
+/* the K&R style to the modern form. */
+
+/*
+void ptitle(char* A)
+*/
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Functions to report projection parameter
+!Input Parameters: none
+!Output Parameters: none
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 09/03/96 */
+/* printf() is a prohibited function. */
+
+/* {printf("\n%s Projection Parameters:\n\n",A);} */
+
+
+
+/*
+void radius(double A)
+*/
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description:
+!Input Parameters: none
+!Output Parameters: none
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 08/27/96 */
+/* The non-standard format specification */
+/* %lf is changed to %f. */
+
+/* Modified by JC Guu 09/03/96 */
+/* The prohibited function ptintf() is */
+/* commented out. */
+
+/*{ printf(" Radius of Sphere: %f meters\n",A); } */
+
+
+
+/* Modified by JC Guu 08/26/96 */
+/* The old K&R style function definition */
+/* is changed to the modern form and the */
+/* return type is given as void. */
+
+void giherror(char* what, char* where)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to report errors
+!Input Parameters: what cause of the error
+ where where the error occurs
+!Output Parameters: none
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 09/03/96 */
+/* printf() is a prohibited function. */
+/* commented out by fhliang on 12/08/97 */
+{
+/* printf("[%s] %s\n",where,what); */
+}
+
+
+#ifndef sun
+
+
+
+/* Modified by JC Guu 08/26/96 */
+/* The old K&R style function definition */
+/* is changed to the modern form and the */
+/* return type is given as void. */
+
+void getcoord_sincos(double val, double* sin_val, double* cos_val)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to calculate the sine and cosine in
+ one call. Some computer systems have implemented
+ this function, resulting in a faster implementation
+ than calling each function separately. It is
+ provided here for those computer systems which
+ don`t implement this function.
+
+!Input Parameters: val the value to be operated on
+!Output Parameters: sin_val the sine of val
+ cos_val the cosine of val
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 08/26/96 */
+/* A return statement is not needed here. */
+
+{ *sin_val = sin(val); *cos_val = cos(val); }
+#endif
+
+
+
+/* Modified by JC Guu 08/27/96 */
+/* The old K&R style function definition */
+/* is changed to the modern form and the */
+/* return type is given as int. */
+
+int getcoord_sign(double x)
+{ if (x < 0.0) return(-1); else return(1);}
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to return the sign of an argument
+!Input Parameters: x a floating point number
+!Output Parameters: the sign of x
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+
+
+/* Modified by JC Guu 08/27/96 */
+/* The old K&R style function definition */
+/* is changed to the modern. */
+
+double getcoord_adjust_lon(double x)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to adjust longitude to -180 to 180
+!Input Parameters: x
+!Output Parameters:
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+/* Modified by JC Guu 08/27/96 */
+/* The explicit cast is provided. */
+
+{x=(fabs(x)<PI)?x:(x-(((double) getcoord_sign(x))*TWO_PI));return(x);}
+
+
+
+/* Modified by JC Guu 08/27/96 */
+/* The style of function definition is changed */
+/* from the K&R style to the modern form for */
+/* the following functions. */
+
+double getcoord_e0fn(double x)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to compute e0
+!Input Parameters: x
+!Output Parameters:
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+{return(1.0-0.25*x*(1.0+x/16.0*(3.0+1.25*x)));}
+
+
+
+double getcoord_e1fn(double x)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to compute e1
+!Input Parameters: x
+!Output Parameters:
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+{return(0.375*x*(1.0+0.25*x*(1.0+0.46875*x)));}
+
+
+
+double getcoord_e2fn(double x)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description: Function to compute e2
+!Input Parameters: x
+!Output Parameters:
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+{return(0.05859375*x*x*(1.0+0.75*x));}
+
+
+
+double getcoord_mlfn(double e0,double e1,double e2,double phi)
+
+/* Modified by JC Guu 08/30/96 */
+/* The required prolog components are added */
+/* here. */
+
+/*
+!C
+!Description:
+!Input Parameters: e0, e1, e2 and phi
+!Output Parameters:
+!Revision History:
+!Team-unique Header:
+!END
+*/
+
+{
+return(e0*phi-e1*sin(2.0*phi)+e2*sin(4.0*phi));
+}
diff --git a/src/read_GEOS.c b/src/read_GEOS.c
new file mode 100644
index 0000000..5e94a41
--- /dev/null
+++ b/src/read_GEOS.c
@@ -0,0 +1,219 @@
+/******************************************************************************
+
+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);
+
+/******************************************************************************/
+
+}
diff --git a/src/read_GEOS_constants.c b/src/read_GEOS_constants.c
new file mode 100644
index 0000000..9b8d06d
--- /dev/null
+++ b/src/read_GEOS_constants.c
@@ -0,0 +1,200 @@
+/******************************************************************************
+
+Description:
+
+ Integer function read_GEOS_constants.c
+ Reads desired variables from time-invariant (constant) file (calling argument).
+
+Input arguments:
+ geos_file_cnst input GEOS constants file name
+
+Output arguments:
+ none
+
+Function output:
+ int return_code successful completion is zero, otherwise non-zero
+
+ Output variables in ancillary.h
+
+Revision History:
+ 07/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_constants(char *geos_file_cnst)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+/******************************************************************************/
+
+ size_t lat_len, lon_len, start[3], edge[3];
+ int status, idf, latln_id, lonln_id, nxc, nyc, nxnyc, 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 constants file: %s\n", geos_file_cnst);
+ status = nc_open(geos_file_cnst, NC_NOWRITE, &idf);
+ if (status != NC_NOERR) {
+ printf("Could not open GEOS constants file: %s %d %d\n", geos_file_cnst,
+ 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_cnst);
+ 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_cnst);
+ return_code = -3;
+ return (return_code);
+ }
+ status = nc_inq_dimlen(idf, lonln_id, &lon_len);
+
+/******************************************************************************/
+
+ nxc = lon_len;
+ nyc = lat_len;
+ nxnyc = lat_len * lon_len;
+
+ if ((g_frocn = (float *) malloc(sizeof(float) * nxnyc)) == NULL) {
+ printf("Could not allocate memory for g_frocn\n");
+ return(-4);
+ }
+ if ((g_frlandice = (float *) malloc(sizeof(float) * nxnyc)) == NULL) {
+ printf("Could not allocate memory for g_frlandice\n");
+ return(-4);
+ }
+ if ((g_phis = (float *) malloc(sizeof(float) * nxnyc)) == NULL) {
+ printf("Could not allocate memory for g_glst\n");
+ return(-5);
+ }
+
+/******************************************************************************/
+
+ edge[0] = 1;
+ edge[1] = lat_len;
+ edge[2] = lon_len;
+
+/******************************************************************************/
+
+// Read ocean fractions.
+
+ status = nc_inq_varid (idf, "FROCEAN", &varid);
+ if(status != NC_NOERR) {
+ printf("Could not find FROCEAN variable in %s\n", geos_file_cnst);
+ return_code = -8;
+ return (return_code);
+ }
+
+ status = nc_get_vara_float(idf, varid, start, edge, g_frocn);
+ if (status != NC_NOERR) {
+ printf("could not get FROCEAN data in %s\n", geos_file_cnst);
+ printf("status: %d\n", status);
+ return_code = -9;
+ return (return_code);
+ }
+
+/******************************************************************************/
+
+// Read permanent land ice fractions.
+
+ status = nc_inq_varid (idf, "FRLANDICE", &varid);
+ if(status != NC_NOERR) {
+ printf("Could not find FRLANDICE variable in %s\n", geos_file_cnst);
+ return_code = -8;
+ return (return_code);
+ }
+
+ status = nc_get_vara_float(idf, varid, start, edge, g_frlandice);
+ if (status != NC_NOERR) {
+ printf("could not get FROCEAN data in %s\n", geos_file_cnst);
+ printf("status: %d\n", status);
+ return_code = -9;
+ return (return_code);
+ }
+
+/******************************************************************************/
+
+// Read surface geopotential heights.
+
+ status = nc_inq_varid (idf, "PHIS", &varid);
+ if(status != NC_NOERR) {
+ printf("Could not find PHIS variable in %s\n", geos_file_cnst);
+ return_code = -10;
+ return (return_code);
+ }
+
+ status = nc_get_vara_float(idf, varid, start, edge, g_phis);
+ if (status != NC_NOERR) {
+ printf("could not get PHIS data in %s\n", geos_file_cnst);
+ printf("status: %d\n", status);
+ return_code = -11;
+ return (return_code);
+ }
+
+/******************************************************************************/
+
+// printf("GEOS frc ocn data: %f %f %f\n", g_frocn[0], g_frocn[9], g_frocn[nxnyc-1]);
+// printf("GEOS PHIS data: %f %f %f\n", g_phis[0], g_phis[9], g_phis[nxnyc-1]);
+// fflush;
+
+/******************************************************************************/
+
+// Close files.
+ status = nc_close(idf);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/read_GEOS_lndocn.c b/src/read_GEOS_lndocn.c
new file mode 100644
index 0000000..88dc7ba
--- /dev/null
+++ b/src/read_GEOS_lndocn.c
@@ -0,0 +1,187 @@
+/******************************************************************************
+
+Description:
+
+ Integer function read_GEOS_lndocn.c
+ Reads desired variables from input files (in calling argument list).
+
+Input arguments:
+ geos_file_lnd input GEOS land property file name
+ geos_file_ocn input GEOS ocean property file name
+
+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_lndocn(char *geos_file_lnd, char *geos_file_ocn)
+
+{
+
+/******************************************************************************/
+
+// Declarations and initializations.
+
+/******************************************************************************/
+
+ size_t lat_len, lon_len, start[3], edge[3];
+ int status, idfl, idfo, 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 files.
+// printf("Opening GEOS land file: %s\n", geos_file_lnd);
+ status = nc_open(geos_file_lnd, NC_NOWRITE, &idfl);
+ if (status != NC_NOERR) {
+ printf("Could not open GEOS land property file: %s %d %d\n", geos_file_lnd,
+ idfl, status);
+ return (-1);
+ }
+// printf("Opening GEOS ocean file: %s\n", geos_file_ocn);
+ status = nc_open(geos_file_ocn, NC_NOWRITE, &idfo);
+ if (status != NC_NOERR) {
+ printf("Could not open GEOS ocean property file: %s %d %d\n", geos_file_ocn,
+ idfo, status);
+ return (-1);
+ }
+
+/******************************************************************************/
+
+// Get data array dimensions. All meteorologial variables are the same size.
+// First dimension (time) is 1.
+
+ status = nc_inq_dimid(idfl, "lat", &latln_id);
+ if(status < 0) {
+ printf("Could not acquire lat-dimension data: %s\n", geos_file_lnd);
+ return_code = -2;
+ return (return_code);
+ }
+ status = nc_inq_dimlen(idfl, latln_id, &lat_len);
+
+ status = nc_inq_dimid(idfl, "lon", &lonln_id);
+ if(status < 0) {
+ printf("Could not acquire lon-dimension data: %s\n", geos_file_lnd);
+ return_code = -3;
+ return (return_code);
+ }
+ status = nc_inq_dimlen(idfl, 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", idfl, latln_id, lonln_id, lat_len, lon_len, nxny);
+
+ if ((g_snow = (float *) malloc(sizeof(float) * nxny)) == NULL) {
+ printf("Could not allocate memory for g_snow\n");
+ return(-4);
+ }
+ if ((g_ice = (float *) malloc(sizeof(float) * nxny)) == NULL) {
+ printf("Could not allocate memory for g_ice\n");
+ return(-5);
+ }
+
+/******************************************************************************/
+
+ edge[0] = 1;
+ edge[1] = lat_len;
+ edge[2] = lon_len;
+
+// Read snow fractions.
+
+ status = nc_inq_varid (idfl, "FRSNO", &varid);
+ if(status != NC_NOERR) {
+ printf("Could not find FRSNO variable in %s\n", geos_file_lnd);
+ return_code = -8;
+ return (return_code);
+ }
+
+ status = nc_get_vara_float(idfl, varid, start, edge, g_snow);
+ if (status != NC_NOERR) {
+ printf("could not get FRSNO data in %s\n", geos_file_lnd);
+ printf("status: %d\n", status);
+ return_code = -9;
+ return (return_code);
+ }
+
+// Read ice fractions.
+
+ status = nc_inq_varid (idfo, "FRSEAICE", &varid);
+ if(status != NC_NOERR) {
+ printf("Could not find FRSEAICE variable in %s\n", geos_file_ocn);
+ return_code = -10;
+ return (return_code);
+ }
+
+ status = nc_get_vara_float(idfo, varid, start, edge, g_ice);
+ if (status != NC_NOERR) {
+ printf("could not get FRSEAICE data in %s\n", geos_file_ocn);
+ printf("status: %d\n", status);
+ return_code = -11;
+ return (return_code);
+ }
+
+/******************************************************************************/
+
+// printf("GEOS snow data: %f %f %f\n", g_snow[0], g_snow[9], g_snow[nxny-1]);
+// printf("GEOS ice data: %f %f %f\n", g_ice[0], g_ice[9], g_ice[nxny-1]);
+// fflush;
+
+/******************************************************************************/
+
+// Close files.
+ status = nc_close(idfl);
+ status = nc_close(idfo);
+
+/******************************************************************************/
+
+ return (return_code);
+
+/******************************************************************************/
+
+}
diff --git a/src/swap_bytes.c b/src/swap_bytes.c
new file mode 100644
index 0000000..47fd705
--- /dev/null
+++ b/src/swap_bytes.c
@@ -0,0 +1,126 @@
+# include <stdlib.h>
+
+int swapbyte4(void *, int);
+void swbyt4_(void *, int *);
+void fbyte4_(void *, int *);
+
+/* ===== swapbyte4 ================================================== */
+/*
+*| Name:
+*| swapbyte4 - Flips 4 byte words into network byte order.
+*|
+*| Interface:
+*| int
+*| swapbyte4 (void *buf, int n_words)
+*|
+*| Input:
+*| buf - Buffer to be flipped.
+*| n_words - Number of words to flip.
+*|
+*| Input and Output:
+*| none
+*|
+*| Output:
+*| none
+*|
+*| Return values:
+*| 0 - Success.
+*| -1 - buf pointer is NULL.
+*|
+*| Remarks:
+*| Network byte order is big endian.
+*|
+*| Categories:
+*| converter
+*| utility
+*/
+
+int
+ swapbyte4 (void *buf, int n_words)
+{
+ int *b_ptr = (int *) buf;
+ int nw = (int) n_words;
+
+ if (b_ptr == (int *) NULL)
+ {
+ return (-1);
+ }
+
+ if (n_words > 0)
+ {
+ swbyt4_ (b_ptr, &nw);
+ }
+ return (0);
+}
+
+/*
+*$ Name:
+*$ swbyt4 - Re-orders bytes if internal representation is not
+*$ big-endian.
+*$
+*$ Interface:
+*$ subroutine
+*$ swbyt4(integer buf(*), integer n)
+*$
+*$ Input:
+*$ n - Number of 4 byte swaps to be made if order is not
+*$ big-endian.
+*$
+*$ Input and Output:
+*$ buf - Array containing bytes to be swapped.
+*$
+*$ Output:
+*$ none
+*$
+*$ Return values:
+*$ none
+*$
+*$ Remarks:
+*$ none
+*$
+*$ Categories:
+*$ utility
+*/
+
+void
+swbyt4_(void *buf, int *n)
+{
+ /* determine byte order from first principles */
+ union
+ {
+ char bytes[sizeof(int)];
+ int word;
+ } q;
+
+ q.word = 1;
+ if (q.bytes[3] == 1)
+ return;
+
+ /* swap bytes */
+ fbyte4_(buf, n);
+}
+
+ /* buffer - INTEGER*4 array to switch bytes */
+ /* n - INTEGER number of 4 byte switches */
+
+void
+fbyte4_(void *buffer, int *num)
+{
+ char *cbuf = (char *)buffer;
+ int i, n;
+
+ n = *num;
+ for (i = 0; i < n; i++)
+ {
+ char b;
+
+ b = cbuf[0];
+ cbuf[0] = cbuf[3];
+ cbuf[3] = b;
+ b = cbuf[1];
+ cbuf[1] = cbuf[2];
+ cbuf[2] = b;
+ cbuf += 4;
+ }
+}
+
--
GitLab