ecearth3/sources/lpjg/cru/guessio/cruncep_1901_2015.h

//////////////////////////////////////////////////////////////////////////////////////
// CRUNCEP_1901_2015.H
// Header file for input from a fast data archive
// Created automatically by FastArchive on Mon Sep 12 18:08:07 2016
//
// The following #includes should appear in your source code file:
//
//   #include <stdio.h>
//   #include <stdlib.h>
//   #include <string.h>
//   #include "/data/cruncep_v7/to_fast_archive_conversion/nc_monthly_files_verJuly1//cruncep_1901_2015.h"
//
// Functionality to retrieve data from the archive is provided by class Cruncep_1901_2015Archive.
// The following public functions are provided:
//
// bool open(char* filename)
//   Attempts to open the specified file as a fast data archive. The format must be
//   exactly compatible with this version of cruncep_1901_2015.h (normally the archive and
//   header file should have been produced together by the same program using class
//   CFastArchive). Returns false if the file could not be opened or had format
//   errors. open() with no argument is equivalent to open("cruncep_1901_2015.bin").
//
// void close()
//   Closes the archive (if open).
//
// bool rewind()
//   Sets the file pointer to the first record in the archive file. Returns false if
//   no archive file is currently open.
//
// bool getnext(Cruncep_1901_2015& obj)
//   Retrieves the next record in the archive file and advances the file pointer to
//   the next record. Data are written to the member variables of obj. Returns false if
//   no archive file is currently open or if the file pointer is beyond the last
//   record. Use rewind() and getnext() to retrieve data sequentially from the archive.
//
// bool getindex(Cruncep_1901_2015& obj)
//   Searches the archive for a record matching the values specified for the index
//   items (lon and lat) in obj. If a matching record is found, the data are
//   written to the member variables of obj. Returns true if the archive was open and
//   a matching record was found, otherwise false. The search is iterative and fast.
//
// Sample program:
//
//   Cruncep_1901_2015Archive ark;
//   Cruncep_1901_2015 data;
//   bool success,flag;
//
//   // Retrieve all records in sequence and print values of lon and lat:
//
//   success=ark.open("cruncep_1901_2015.bin");
//   if (success) {
//      flag=ark.rewind();
//      while (flag) {
//         flag=ark.getnext(data);
//         if (flag)
//            printf("Loaded record: lon=%g, lat=%g\n",data.lon,data.lat);
//      }
//   }
//   
//   // Look for a record with lon=-180, lat=-90:
//
//   data.lon=-180;
//   data.lat=-90;
//   success=ark.getindex(data);
//   if (success) printf("Found it!\n");
//   else printf("Not found\n");
//
//   ark.close();


struct Cruncep_1901_2015 {

	// Index part

	double lon;
	double lat;

	// Data part

	double mtemp[1380];
	double mprec[1380];
	double mswrad[1380];
	double soilcode[1];
};


const long CRUNCEP_1901_2015_NRECORD=59191;
const int CRUNCEP_1901_2015_DATA_LENGTH=12422;
const int CRUNCEP_1901_2015_INDEX_LENGTH=7;
const int CRUNCEP_1901_2015_HEADERSIZE=659;
unsigned char CRUNCEP_1901_2015_HEADER[CRUNCEP_1901_2015_HEADERSIZE-4]={
	0x01,0x02,0x93,0x00,0x00,0x00,0x07,0x00,0x02,0x04,0x6C,0x6F,0x6E,0x00,0x2D,0x31,0x38,0x30,0x00,0x00,
	0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,
	0xD3,0x9D,0x42,0x7F,0x00,0x00,0x31,0x38,0x30,0x00,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,
	0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x30,0x2E,
	0x32,0x35,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,
	0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x0B,0x04,0x6C,0x61,0x74,0x00,0x2D,0x39,0x30,0x00,
	0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,
	0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x39,0x30,0x00,0x00,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,
	0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,
	0x30,0x2E,0x32,0x35,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,
	0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x0A,0x00,0x00,0x30,0x86,0x00,0x04,0x06,
	0x6D,0x74,0x65,0x6D,0x70,0x00,0x2D,0x31,0x30,0x30,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,
	0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x31,0x30,
	0x30,0x00,0x00,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,
	0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x30,0x2E,0x30,0x30,0x31,0x00,0x00,0x00,0xA0,0xB5,
	0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,
	0x00,0x00,0x12,0x00,0x00,0x05,0x64,0x06,0x6D,0x70,0x72,0x65,0x63,0x00,0x30,0x00,0x05,0x64,0x31,0x00,
	0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,
	0xD3,0x9D,0x42,0x7F,0x00,0x00,0x31,0x30,0x30,0x30,0x30,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,
	0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x30,0x2E,
	0x30,0x30,0x30,0x31,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,
	0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x1B,0x00,0x00,0x05,0x64,0x07,0x6D,0x73,0x77,0x72,
	0x61,0x64,0x00,0x30,0x00,0x05,0x64,0x30,0x31,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,
	0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x31,0x30,0x30,0x30,0x30,
	0x30,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,
	0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x30,0x2E,0x30,0x30,0x31,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,
	0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x1B,
	0x00,0x00,0x05,0x64,0x09,0x73,0x6F,0x69,0x6C,0x63,0x6F,0x64,0x65,0x00,0x30,0x00,0x05,0x64,0x31,0x00,
	0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,
	0xD3,0x9D,0x42,0x7F,0x00,0x00,0x35,0x30,0x00,0x64,0x31,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,
	0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x30,0x2E,
	0x31,0x00,0x31,0x00,0x00,0x00,0xA0,0xB5,0xBF,0x01,0x00,0x00,0x00,0x00,0xD2,0x00,0x41,0x00,0x00,0x00,
	0x00,0x00,0x15,0x95,0xD3,0x9D,0x42,0x7F,0x00,0x00,0x09,0x00,0x00,0x00,0x01};


class Cruncep_1901_2015Archive {

private:

	FILE* pfile;
	long recno;
	long datano;
	unsigned char pindex[CRUNCEP_1901_2015_INDEX_LENGTH];
	unsigned char pdata[CRUNCEP_1901_2015_DATA_LENGTH];
	bool iseof;

	long readbin(int nbyte) {

		unsigned char buf[4];
		long mult[4]={0x1000000,0x10000,0x100,1},val;
		int i;

		fread(buf,nbyte,1,pfile);

		val=0;
		for (i=0;i<nbyte;i++) {
			val+=buf[i]*mult[4-nbyte+i];
		}

		return val;
	}

	void getindex(long n) {

		fseek(pfile,CRUNCEP_1901_2015_INDEX_LENGTH*(n-recno),SEEK_CUR);
		fread(pindex,CRUNCEP_1901_2015_INDEX_LENGTH,1,pfile);
		datano=pindex[CRUNCEP_1901_2015_INDEX_LENGTH-4]*0x1000000+pindex[CRUNCEP_1901_2015_INDEX_LENGTH-3]*0x10000+
			pindex[CRUNCEP_1901_2015_INDEX_LENGTH-2]*0x100+pindex[CRUNCEP_1901_2015_INDEX_LENGTH-1];
		recno=n+1;
		iseof=(recno==CRUNCEP_1901_2015_NRECORD);
	}

	void getdata() {

		fseek(pfile,CRUNCEP_1901_2015_INDEX_LENGTH*-recno+(datano-CRUNCEP_1901_2015_NRECORD)*CRUNCEP_1901_2015_DATA_LENGTH,SEEK_CUR);
		fread(pdata,CRUNCEP_1901_2015_DATA_LENGTH,1,pfile);
		fseek(pfile,CRUNCEP_1901_2015_DATA_LENGTH*(CRUNCEP_1901_2015_NRECORD-datano-1)+CRUNCEP_1901_2015_INDEX_LENGTH*recno,SEEK_CUR);
	}

	double popreal(unsigned char* bits,int nbyte,int nbit,double scalar,double offset) {

		unsigned char buf;
		int nb=nbit/8,i;
		double rval=0.0;
		long mult[4]={1,0x100,0x10000,0x1000000};

		for (i=0;i<4;i++) {
			if (i<nb) rval+=bits[nbyte-i-1]*mult[i];
			else if (i==nb) {
				buf=bits[nbyte-i-1]<<(8-nbit%8);
				buf>>=8-nbit%8;
				rval+=buf*mult[i];
			}
		}

		for (i=nbyte-1;i>=0;i--) {
			if (i>=nb)
				bits[i]=bits[i-nb];
			else
				bits[i]=0;
		}

		nb=nbit%8;

		for (i=nbyte-1;i>=0;i--) {
			bits[i]>>=nb;
			if (i>0) {
				buf=bits[i-1];
				buf<<=8-nb;
				bits[i]|=buf;
			}
		}

		rval=rval*scalar+offset;

		return rval;
	}

	void bitify(unsigned char buf[4],double fval,double offset,double scalar) {

		long ival = (long)((fval-offset)/scalar + 0.5);
		buf[0]=(unsigned char)(ival/0x1000000);
		ival-=buf[0]*0x1000000;
		buf[1]=(unsigned char)(ival/0x10000);
		ival-=buf[1]*0x10000;
		buf[2]=(unsigned char)(ival/0x100);
		ival-=buf[2]*0x100;
		buf[3]=(unsigned char)(ival);
	}

	void merge(unsigned char ptarget[3],unsigned char buf[4],int bits) {

		int nb=bits/8;
		int i;
		unsigned char nib;
		for (i=0;i<3;i++) {

			if (i<3-nb)
				ptarget[i]=ptarget[i+nb];
			else
				ptarget[i]=0;
		}
		nb=bits%8;
		for (i=0;i<3;i++) {
			ptarget[i]<<=nb;
			if (i<3-1) {
				nib=ptarget[i+1]>>(8-nb);
				ptarget[i]|=nib;
			}
		}

		nb=bits/8;
		if (bits%8) nb++;
		for (i=1;i<=nb;i++)
			ptarget[3-i]|=buf[4-i];
	}

	int compare_index(unsigned char* a,unsigned char* b) {

		int i;
		for (i=0;i<3;i++) {
			if (a[i]<b[i]) return -1;
			else if (a[i]>b[i]) return +1;
		}

		return 0;
	}

	bool initialise(const char* filename) {

		int i;
		unsigned char* pheader;

		if (pfile) fclose(pfile);
		pfile=fopen(filename,"rb");
		if (!pfile) {
			printf("Could not open %s for input\n",filename);
			return false;
		}

		pheader=new unsigned char[CRUNCEP_1901_2015_HEADERSIZE-4];
		if (!pheader) {
			printf("Out of memory\n");
			fclose(pfile);
			pfile=NULL;
			return false;
		}
		::rewind(pfile);
		fread(pheader,CRUNCEP_1901_2015_HEADERSIZE-4,1,pfile);
		for (i=0;i<CRUNCEP_1901_2015_HEADERSIZE-4;i++) {
			if (pheader[i]!=CRUNCEP_1901_2015_HEADER[i]) {
				printf("Format of %s incompatible with this version of cruncep_1901_2015.h\n",filename);
				fclose(pfile);
				pfile=NULL;
				delete pheader;
				return false;
			}
		}
		delete[] pheader;

		::rewind(pfile);
		fseek(pfile,CRUNCEP_1901_2015_HEADERSIZE+CRUNCEP_1901_2015_DATA_LENGTH*CRUNCEP_1901_2015_NRECORD,SEEK_CUR);
		recno=0;
		iseof=false;

		return true;
	}

public:

	Cruncep_1901_2015Archive() {
		pfile=NULL;
	}

	~Cruncep_1901_2015Archive() {
		if (pfile) fclose(pfile);
	}

	bool open(const char* filename) {
		return initialise(filename);
	}

	bool open() {
		return open("cruncep_1901_2015.bin");
	}

	void close() {
		if (pfile) {
			fclose(pfile);
			pfile=NULL;
		}
	}

	bool rewind() {

		if (!pfile) return false;

		::rewind(pfile);
		fseek(pfile,CRUNCEP_1901_2015_HEADERSIZE+CRUNCEP_1901_2015_DATA_LENGTH*CRUNCEP_1901_2015_NRECORD,SEEK_CUR);
		recno=0;
		iseof=false;

		return true;
	}

	bool getnext(Cruncep_1901_2015& obj) {

		if (!pfile || iseof) return false;

		int i;

		getindex(recno);
		getdata();

		obj.lat=popreal(pindex,3,10,0.25,-90);
		obj.lon=popreal(pindex,3,11,0.25,-180);

		for (i=0;i>=0;i--) obj.soilcode[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,9,0.1,0);
		for (i=1379;i>=0;i--) obj.mswrad[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,27,0.001,0);
		for (i=1379;i>=0;i--) obj.mprec[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,27,0.0001,0);
		for (i=1379;i>=0;i--) obj.mtemp[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,18,0.001,-100);

		return true;
	}

	bool getindex(Cruncep_1901_2015& obj) {

		if (!CRUNCEP_1901_2015_NRECORD || !pfile) return false;

		// else

		unsigned char ptarget[3]={0,0,0};
		unsigned char buf[4];
		bitify(buf,obj.lon,-180,0.25);
		merge(ptarget,buf,11);
		bitify(buf,obj.lat,-90,0.25);
		merge(ptarget,buf,10);

		long begin = 0;
		long end = CRUNCEP_1901_2015_NRECORD;

		while (begin < end) {
			long middle = (begin+end)/2;

			getindex(middle);
			getdata();

			int c = compare_index(pindex, ptarget);

			if (c < 0) {
				begin = middle + 1;
			}
			else if (c > 0) {
				end = middle;
			}
			else {

				for (int i=0;i>=0;i--) obj.soilcode[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,9,0.1,0);
				for (int i=1379;i>=0;i--) obj.mswrad[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,27,0.001,0);
				for (int i=1379;i>=0;i--) obj.mprec[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,27,0.0001,0);
				for (int i=1379;i>=0;i--) obj.mtemp[i]=popreal(pdata,CRUNCEP_1901_2015_DATA_LENGTH,18,0.001,-100);

				return true;
			}
		}

		return false;
	}
};