/******************************************************************************
       WbudgetRas.c -Water Budget Raster Version-       A.Kondoh 2001-2002
 ******************************************************************************/

#include        <stdio.h>
#include        <stdlib.h>
#include        <fcntl.h>
#include        <sys\types.h>
#include        <sys\stat.h>
#include        <process.h>
#include        <dos.h>
#include        <io.h>
#include        <share.h>
#include        <math.h>
#include        <string.h>
#include        <ctype.h>
#include        <conio.h>
#include        <errno.h>

#define         FILE_ERR        -1
#define         SE              stderr

#define		TOLERANCE	(0.1)	/* allowance of iteration error */
#define         WHC_FACTOR      (3.0)

unsigned char	buffer[720];
short		buffer2[720];

/* input */

short	Pr[13][360][720];   /* monthly precipitation */
short	Ep[13][360][720];   /* monthly potential evaporation */

/* output */

short	Ea[13][360][720];         /* Actual Evapotranspiration */
short   Sm[13][360][720];   /* Soil water storage */
short   Surplus[13][360][720];    /* Water surplus */
short   Deficit[13][360][720];    /* Water deficit */

/* working */

float	mPr[13];	/* precipitation */
float	mEp[13];	/* potential evaporation */
float	mEa[13];	/* actual evaporation */
float	mSurplus[13];   /* water surplus */
float	mDeficit[13];   /* water deficit */
float	mSm[13];        /* Soil water */
float	whc;		/* water holding capacity */

void	WaterBudget();
void	get_monthly_data();

/*-----------------------------------------------------------------------------
        Main
  -----------------------------------------------------------------------------*/

void    main(argc,argv)
int     argc;
char    *argv[];
{
        int     month;
        int     i,j;
	FILE    *fi;              /* parameter file */
	int	fpr,fep,fhc;
	int	fea,fsm,fsp,fdf;
	char	fname[80];

/* open parameter file */

        if( (fi=fopen(argv[1],"rt"))==NULL){
                fprintf(SE,"Can't open %s.\n",argv[1]);
                exit(-1);
        }

/* open input file */

     /* 1. Precipitation Jan. to Dec. BSQ (short) */

	fscanf(fi,"%s",fname);
        if( (fpr=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

     /* 2. Potential Evaporation Jan. to Dec. BSQ (short) */

	fscanf(fi,"%s",fname);
        if( (fep=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

     /* 3. Soil Water Holding Capacity */

	fscanf(fi,"%s",fname);
        if( (fhc=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

/* open output file */

     /* 4. Actual evaporation */

	fscanf(fi,"%s",fname);
        if( (fea=open(fname,O_RDWR | O_TRUNC | O_CREAT | O_BINARY,
                        S_IREAD | S_IWRITE))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

     /* 5. Soil water storage */

	fscanf(fi,"%s",fname);
        if( (fsm=open(fname,O_RDWR | O_TRUNC | O_CREAT | O_BINARY,
                        S_IREAD | S_IWRITE))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

     /* 6. Water surplus */

	fscanf(fi,"%s",fname);
        if( (fsp=open(fname,O_RDWR | O_TRUNC | O_CREAT | O_BINARY,
                        S_IREAD | S_IWRITE))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

     /* 7. Water deficit */

	fscanf(fi,"%s",fname);
        if( (fdf=open(fname,O_RDWR | O_TRUNC | O_CREAT | O_BINARY,
                        S_IREAD | S_IWRITE))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

/* calculation */

	for(i=0;i<360;i++){
                if( !(i%10) ){printf("*");}

		for(j=0;j<720;j++){

	             /* set input variables */

		        get_monthly_data(fpr,fep,fhc,i,j);

                        whc=whc*WHC_FACTOR;

                        if(mEp[1]>=0.0){   /* not ocean */

                             /* get water budget on the land */

        	        	WaterBudget(i,j);

                             /* layout to output array (0.1mm) */

	                	for(month=1;month<=12;month++){
                                        if(mEa[month]<0.0){
                                                Ea[month][i][j]=0;
                                        }else{
        		                	Ea[month][i][j]=(short)(mEa[month]*10.0+0.5);   /* 0.1mm */
                                        }
			                Surplus[month][i][j]=(short)(mSurplus[month]*10.0+0.5);
			                Deficit[month][i][j]=(short)(mDeficit[month]*10.0+0.5);
			                Sm[month][i][j]=(short)(mSm[month]*10.0+0.5);
		                }
                        }else{ /* ocean */

	                	for(month=1;month<=12;month++){
		                	Ea[month][i][j]=-999;
			                Surplus[month][i][j]=-999;
			                Deficit[month][i][j]=-999;
			                Sm[month][i][j]=-999;
		                }
                        }
                }
	}

/* output */

        printf("\nWriting!\n");
	for(month=1;month<=12;month++){
		for(i=0;i<360;i++){
			write(fea,Ea[month][i],720*sizeof(short));
		}
	}

	for(month=1;month<=12;month++){
		for(i=0;i<360;i++){
			write(fsp,Surplus[month][i],720*sizeof(short));
		}
	}

	for(month=1;month<=12;month++){
		for(i=0;i<360;i++){
			write(fdf,Deficit[month][i],720*sizeof(short));
		}
	}

	for(month=1;month<=12;month++){
		for(i=0;i<360;i++){
			write(fsm,Sm[month][i],720*sizeof(short));
		}
	}

}/* end of main*/

/*-----------------------------------------------------------------------------
        Water budget calculation 
  -----------------------------------------------------------------------------*/

void	WaterBudget(
int     i,
int     j)
{
        int	month, previous_month;
        int	count;
        float	Previous_Sm[13]; /* Soil Moisture of Previous Iteration */
        float	mAWL[13];        /* acumulated potential water loss */
        float   sum;

/* initialize */

	for(month=1;month<=12;month++){
		Previous_Sm[month]=3000.0; /* Soil Moisture of Previous Iteration */
		mSm[month]=0.0;		  /* Soil Mositure */
		mAWL[month]=0.0;	  /* Accumulated Potential Water Loss */
	}

/* iteration */

        count=0;
	while(1){

		for(month=1;month<=12;month++){

		     /* adjust month */

			if(month==1){
				previous_month=12;
			}else{
				previous_month=month-1;
			}

		     /* soil water storage */

			if(mPr[month]<mEp[month]){
				mAWL[month]=mAWL[previous_month]+mEp[month]-mPr[month];
				mSm[month]=whc*exp(-mAWL[month]/whc);
			}else{
				mSm[month]=mSm[previous_month]+(mPr[month]-mEp[month]);
			}
                        if(mSm[month]>whc){
                                mSm[month]=whc;
                        }

     	     	/* actual evaporation */

			if(mPr[month]<mEp[month]){
				mEa[month]=mPr[month]-(mSm[month]-mSm[previous_month]);
			}else{
				mEa[month]=mEp[month];
			}

                        if(mEa[month]<0.0){
                                mEa[month]=0.0;
                        }

     	     	/* water surplus, deficit */

			if(mSm[month]>=whc){
				mSurplus[month]=mPr[month]-mEa[month]-(mSm[month]-mSm[previous_month]);
				mDeficit[month]=0.0;
			}else{
				mSurplus[month]=0.0;
				mDeficit[month]=mEp[month]-mEa[month];
			}

		}/* next month */

	      /* judge */

		sum=0.0;
		for(month=1;month<=12;month++){
			sum += Previous_Sm[month]-mSm[month];
		}

		if(fabs(sum)<TOLERANCE){
			break;
		}else{
			for(month=1;month<=12;month++){
				Previous_Sm[month]=mSm[month];
			}
		}

	}/* wend */

}/* end of WaterBudget */

/*-----------------------------------------------------------------------------
	monthly time series of the pixel
  -----------------------------------------------------------------------------*/

void	get_monthly_data(
int	fpr,
int	fep,
int	fhc,
int	line,
int	pixel)
{
	int	month,l,p;
	long	offset, buffer_size;

	buffer_size=(long)720*(long)sizeof(short);

	for(month=1;month<=12;month++){
		offset=(long)(month-1)*buffer_size*(long)360
                      +(long)line*buffer_size;

		lseek(fpr,offset,SEEK_SET);
		lseek(fep,offset,SEEK_SET);

		read(fpr,buffer2,720*sizeof(short));
		mPr[month]=(float)buffer2[pixel];

		read(fep,buffer2,720*sizeof(short));
		mEp[month]=(float)buffer2[pixel]/10.0;
	}

	offset=(long)line*(long)720;
	lseek(fhc,offset,SEEK_SET);
	read(fhc,buffer,720);
	whc=(float)buffer[pixel];

        if(whc>=250.0){

                for(l=line-1;l<=line+1;l++){
                	offset=(long)l*(long)720;
	                lseek(fhc,offset,SEEK_SET);
	                read(fhc,buffer,720);
                        for(p=pixel-1;p<=pixel+1;p++){
                                whc=(float)buffer[p];
                                if(whc<250.0){
                                        goto end_whc;
                                }
                        }
                }
        }

end_whc:;

}/* end of get_monthly_data */