get_Reynolds_SST.c

/*******************************************************************************
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);

}