/*************************************************************************
  Program:     	beps.c (Boreal Ecosystem Productivity Simulator)
  --------
 
  Description:
  -----------
  Calculate key components in carbon and water cycles for a specified area. The key components include:
  NPP
  Plant respiration
  Evaporation
  Transpiration
 ***************************************************************************
  Written originally by:   J. Liu/D. Fraser
  Modified for application to China by:  Gang Mo
  Last update:	Spet. 2008
*****************************************************************************/

#include "beps.h"
#include <string.h>

int main()
{
    char cntrl_fn[50];
	int long j=1;			/* Counters: j for output counter */
	int long line,pix,lin;		/* Counters */
	int long jday;			/* Julian day */
//	int short lastlc;		/* Last landcover value */
	int short rcode;		/* Error code */
	double b[SIZEB];		/* Biological parameters */
	int short *lc;			/* Land cover for a line in a given study area */
	double *lai;			/* LAI for a line in a given study area */
	double *awc;   			/* Av soil water hold cap for a line in a given study area */	
	double *lon,*lat;		/* Long, lat values for a line in study area */
	double *x;				/* X variables */
	double *npp;			/* NPP */
	float *output;	/* For output images */

	struct xvalue *xx;		/* X variables for a line in a given study area */
	struct climatedata *sdat;	/* climate data */
	
/* 	Read Control parameters */
	strcpy(cntrl_fn,"controlfile.txt");
	readconf(cntrl_fn,&rcode);
	if (rcode == ERROR)
          exit(0);

/*	Allocate memory */
	lc=(int short*)malloc(npixels*sizeof(int short));
	awc=(double*)malloc(npixels*sizeof(double));	
	lai=(double*)malloc(npixels*sizeof(double));
	lon=(double*)malloc(npixels*sizeof(double));
	lat=(double*)malloc(npixels*sizeof(double));
	npp=(double*)malloc(npixels*sizeof(double));
	output=(float*)malloc(npixels*sizeof(float));

/*	Daily Model Parameters  (Currently there are 12) */
	x=(double*)malloc(SIZEX*sizeof(double));

/*	Memory for model parameters which change daily */
	xx=(struct xvalue*)malloc(npixels*sizeof(struct xvalue));
		
	sdat=(struct climatedata*)malloc(npixels*sizeof(struct climatedata));

	z= (double *)malloc(SIZEZ*sizeof(double));	/* z variables */


/*	Read default biological paramenters */
	readb_init(b);

/* LLLLLLLLLL	For each line in Study Area 	LLLLLLLLL*/
	for (lin = 0; lin <nlines; lin++) {

/*	  Calculate true line in study area (referenced to 0 origin)*/
	  line = lin_offset + lin;

	  printf("Simulating line %d...\n", line+1);	

/*	  Read longitude and latitude values */
		readlonlat(line,lon,lat);

/*	  Clean memory */
	  for(pix=0; pix<npixels; pix++) {
		npp[pix]=0.0;
		zeroxx(pix,x,xx);
		output[pix]= 0;
	  }

/* 	  Read landcover values (for each pixel in the line) */
	  readlc(line,lc,&rcode); 
	  if (rcode == ERROR) {
	    exit(0);
	  }
			
/*	  Read Available Water Capacity data */
	  readawc(line,awc,&rcode); 	/* in m */

/*	  Read Snow pack data 
	  readsnow(line,xx,rcode);
	  if (rcode == ERROR) {
	    exit(0);
	  }
*/
/*	  Read biomass data */
	  readbm(line,lc,xx,&rcode); 
	  if (rcode == ERROR) {
	    exit(0);
	  }
    
/* DDDDDDDDDD	For each day  DDDDDDDDDD */
	  
 	  for (jday=jday_start; jday<=jday_end; jday++) {  
		            
           /*Read Leaf Area Index*/
	       readlai(line,lai,jday,&rcode); 
	       if (rcode == ERROR)
	       exit(0);

           /*read meteo. data */
	       readclim(line,jday,sdat,&rcode); 
	       if (rcode == ERROR) {
	          exit(0);
		   } 
	       else ; 

           /* PPPPPPPPPP	For each PIXEL in A Line PPPPPPPPPP */
	       for (pix=0; pix<npixels; pix++) 
		   {
			   
            /*If lack of data */
			if (sdat[pix].srad<-1 || sdat[pix].tmax<-69 || sdat[pix].tmin<-89 || sdat[pix].precip<-1 || sdat[pix].dewp<-1)
				lc[pix]=0;

            /*If landcover is a value for simulating */
	        if(lc[pix]>0) 
				{
                  /* Read parameters according to cover type */
	                readb(b,lc[pix]);		
				
				   /* Model calculation for each pixel */
                    model(jday,pix,jday_start,lat[pix],lai[pix],awc[pix],lc[pix],x,b,sdat,xx);					
				}		   					
			}	
           
		   /* PPPPPPPPPP	End of PIXEL cycle PPPPPPPPPP */
//	       if(jday % cycle_display == 0) {
			   if(output_GPP=='y' || output_GPP=='Y') {                
				   for(pix=0; pix<npixels; pix++) {
					  output[pix]=(float)(xx[pix].x11*1000);	
				   }
		           sprintf(outf_name,"%s/gpp%d.img",outp_path,jday);
	               display(outf_name,output,&rcode);
	               if (rcode == ERROR) exit(0);
			   }
	     
			   if(output_NPP=='y' || output_NPP=='Y') {                
				   for(pix=0; pix<npixels; pix++) {
					  output[pix]=(float)(xx[pix].x6*1000);	
				   }
		           sprintf(outf_name,"%s/npp%d.img",outp_path,jday);
	               display(outf_name,output,&rcode);
	               if (rcode == ERROR) exit(0);
			   }

			   if(output_NEP=='y' || output_NEP=='Y') {                
				   for(pix=0; pix<npixels; pix++) {
					  output[pix]=(float)(xx[pix].x6*1000-xx[pix].x12*1000);	
				   }
		           sprintf(outf_name,"%s/nep%d.img",outp_path,jday);
	               display(outf_name,output,&rcode);
	               if (rcode == ERROR) exit(0);
			   }

			   if(output_ET=='y' || output_ET=='Y') {
				   for(pix=0; pix<npixels; pix++) {
					output[pix]=(float)(xx[pix].x4*1000+xx[pix].x5*1000);
				   } 
		           sprintf(outf_name,"%s/et%d.img",outp_path,jday);      
	               display(outf_name,output,&rcode);
	               if (rcode == ERROR) exit(0);
			   }  

//               j++;
//      	    }	/* End of Display */

	  }
/* DDDDDDDDDD	End of daily cycle DDDDDDDDDD*/			    	  	
//	  j=1;
	        
	}	
/* LLLLLLLLLL	End of lin_cycle LLLLLLLLLL*/
       
/*	Free Memory */
    free(lc);
	free(awc);	
	free(lai);
	free(lon);
	free(lat);
	free(npp);
	free(output);
	free(x);
	free(xx);
	free(z);

	printf("\n ***End of BEPS***\n");

	return 1;
}	/* end of main */

