/*----------------------------------------------------------------------*\
 | Routines to read sea-ice files for use in Fortran programs.          |
 |                                                                      |
 | Public functions:                                                    |
 |  iceopen     open an sst file (contents go to a static array)        |
 |  icexy       retrieve temperature at some x,y position               |
 |  icell       retrieve ice cover at some latitude and longitude       |
 |  iceclose    free space used for sst data                            |
 |                                                                      |
 | These procedures call procedures in read_ice.c                       |
 |                                                                      |
 | Porting instructions:                                                |
 |      Fortran compilers differ in how they expect Fortran-callable    |
 |      C procedures to be named.  I have coded these routines so that  |
 |      the C name is followed by an underscore.  This convention is    |
 |      used by the Green Hills Fortran compiler and may be used by     |
 |      some others.  If you are running UNIX, try f77 -S on a dummy    |
 |      procedure, then look carefully at the name given in the         |
 |      assembly-language output.  This, combined with and compared     |
 |      with the results of the same exercise in C, will help to        |
 |      identify the proper name modification needed to make these      |
 |      routines visible to Fortran.                                    |
 |                                                                      |
 | Fortran declarations:                                                |
 |      INTEGER ICEOPEN                                                 |
 |      INTEGER ICEXY                                                   |
 |      INTEGER ICELL                                                   |
 |      none needed for iceclose, since this is a procedure             |
 |                                                                      |
 | Fortran usage:                                                       |
 |      OK = ICEOPEN (FILENAME)                                         |
 |      OK = ICEXY (IX,IY,T)                                            |
 |      CALL ICECLOSE                                                   |
 |      above assumes that OK, IX, and IY are declared INTEGER, T is    |
 |      declared as REAL*8, and FILENAME is declared CHARACTER*256      |
 |      (FILENAME can be smaller than 256)                              |
 |                                                                      |
 |      ICEXY stores a result in T, and returns status.  If the point   |
 |      given by IX and IY is not in the image, status is 0.  If it is  |
 |      but the pixel at that point indicates missing data or land,     |
 |      the status returned is 1 but the value stored in T is -99.0.    |
 |      So a value in T of -99.0 always means no datum was available,   |
 |      and the value returned (0 or 1) indicates whether this is       |
 |      because the point position is invalid or because the image has  |
 |      no ice value for that pixel.                                    |
 |                                                                      |
 |      ICELL does what ICEXY does, but the first two arguments are     |
 |      double precision latitude and longitude rather than pixel       |
 |      coordinates.                                                    |
 |                                                                      |
 | Peter N. Schweitzer (U.S. Geological Survey, Reston, VA 22092)       |
\*----------------------------------------------------------------------*/

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

extern int DFR8getdims (char *name, long *width, long *height, short *haspalette);
extern int DFR8getimage (char *name, unsigned char *dst, long w, long h, unsigned char *pal);

#define MISSING		157
#define LAND		168

static unsigned char *ice = NULL;
static long width = 0, height = 0;
static char hemisphere = 0;

static int nint (double x) {
	return ((x - floor(x) < 0.5) ? (int) x : (int)(x + 0.5));
	}

/*----------------------------------------------------------------------*\
 | Code to retrieve coordinates given latitude and longitude			|
\*----------------------------------------------------------------------*/

extern void mapll (double latitude, double longitude, double *x, double *y);
extern void mapxy (double x, double y, char hemisphere, double *latitude, double *longitude);

static void xytoll (int x, int y, double *lat, double *lon) {
	double ax,ay;

	ax = ((double) x) * 25.0 - (((hemisphere == 'N') ? 3850.0 : 3950.0) - 12.5);
	ay = ((double) height - y) * 25.0 - (((hemisphere == 'N') ? 5350.0 : 3950.0) - 12.5);
	mapxy (ax,ay,hemisphere,lat,lon);
	}

void xytoll_ (int *x, int *y, double *lat, double *lon) {
	xytoll (*x,*y,lat,lon);
	}

static void lltoxy (double lat, double lon, int *x, int *y) {
	double ax,ay;

	lon += (hemisphere == 'N') ? 45.0 : 180.0;
	if (lon > 360.0) lon -= 360.0;
	mapll (lat,lon,&ax,&ay);
	*x = nint ((ax + ((hemisphere == 'N') ? 3850.0 : 3950.0) - 12.5)/25.0);
	*y = nint ((ay + ((hemisphere == 'N') ? 5350.0 : 3950.0) - 12.5)/25.0);
	*y = height - *y;
	}

void lltoxy_ (double *lat, double *lon, int *x, int *y) {
	lltoxy (*lat,*lon,x,y);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/

char *iceopen_ (char *name, int name_len) {
	int i;
	long w,h;
	short haspalette;

	if (name_len > 0) {
		for (i=0; i < name_len; i++)
			if (isspace (name[i])) {
				name[i] = 0;
				break;
				}
		if (DFR8getdims (name,&w,&h,&haspalette) == 0) {
			if (w != width || h != height) {
				width = w;
				height = h;
				if (ice) free (ice);
				if (!(ice = (unsigned char *) malloc (width*height))) {
					fprintf (stderr,"Error: could not allocate ice array to size %ld\n",width*height);
					return (NULL);
					}
				}
			if (DFR8getimage (name,ice,w,h,NULL) != 0) {
				fprintf (stderr,"Error: could not read ice data from file %s\n",name);
				return (NULL);
				}
			if (height == 448) hemisphere = 'N';
			else
				if (height == 332) hemisphere = 'S';
				else
					fprintf (stderr,"Warning: expected image height of 448 (N) or 332 (S), got %ld\n",height);
			}
		else {
			fprintf (stderr,"Error: could not get image dimensions from file %s\n",name);
			return (NULL);
			}
		return (ice);
		}
	else
		return (NULL);
	}

void iceclose_ (void) {
	if (ice) {
		free (ice);
		ice = NULL;
		}
	}

int icexy_ (int *x, int *y, double *value) {
	unsigned pixel;

	*value = -99.0;
	if (*x < 0 || *x >= width || *y < 0 || *y >= height) return (0);
	pixel = ice [*y * width + *x];
	if (pixel == MISSING || pixel == LAND) return (1);
	*value = (double) pixel;
	return (1);
	}

int icell_ (double *latitude, double *longitude, double *value) {
	int x,y;
	unsigned pixel;

	lltoxy (*latitude,*longitude,&x,&y);
	*value = -99.0;
	if (x < 0 || x >= width || y < 0 || y >= height) return (0);
	pixel = ice [y * width + x];
	if (pixel == MISSING || pixel == LAND) return (1);
	*value = (double) pixel;
	return (1);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/