/********************************************************************
        Priestley-Taylor法による可能蒸発散量の計算

                                                A.Kondoh 1994
 ********************************************************************/

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

#define         BUF_SIZE        4320

#define         RO              (1.96)            /* 太陽定数(ly/min,cal/cm2min) */
#define         SC              (1382.0)          /* 太陽定数(W/m") */
#define         PI              3.141592653
#define         D2R             1.745329E-02      /* 度→ﾗｼﾞｱﾝ */
#define         R2D             57.29577952       /* ﾗｼﾞｱﾝ -> 度 */
#define         A               1.0E03            /* W/m" -> mm/day への変換係数  */

#define         SE              stderr
#define         CRLF            puts("")
#define         PAUSE           do{}while(!isspace(getch())) /* ﾎﾟｰｽﾞ */
#define         FILE_ERR        -1

int             nday[13][32];        /* 各月の日数 */
int             buffer[BUF_SIZE];    /* バッファー */
signed short    ET[BUF_SIZE];        /* 蒸発量計算結果 */
unsigned char   buffer1[BUF_SIZE];   /* １バイト整数バッファ */
short           buffer2[BUF_SIZE];   /* ２バイト整数バッファ */

/* 関数 */

float           PT_method(int,double,double,double,double,double,double,int,int); /* PT法の計算 */
double          sekii();             /* 赤緯 */
double          get_lat();           /* 緯度をラジアンで返す */
double          get_lon();           /* 経度をラジアンで返す */
double          Qa_calc();           /* 大気圏外日射量 */
double          Sun();               /* 可照時間の計算 */
double          cal_delta();         /* Δの計算 */
double          cal_gamma();         /* γの計算 */
double          cal_ramda();         /* 蒸発の潜熱の計算 */
void            cal_nday();          /* １月０日からの日数 */
double          get_temperature();   /* 気温の読み取り */
double          get_albedo();        /* アルベドの読み取り */
double          get_cloudiness();    /* 日照率の読み取り */
double          get_elevation();     /* 標高の読み取り */
double          get_precipitation(); /* 降水量の読み取り */

static int      nmonth[13]={0,31,28,31,30,31,30,31,31,30,31,30,31};

FILE            *fp;                            /* ﾊﾟﾗﾒｰﾀｰﾌｧｲﾙ */
int             fat,fal,fel,fcl,fpr,fo;        /* ﾌｧｲﾙﾃﾞｨｽｸﾘﾌﾟﾀ */

/*-------------------------------------------------------------------------------------------------
        function main()
  -------------------------------------------------------------------------------------------------*/

void    main(argc,argv)
int     argc;
char    *argv[];
{
        char    fname[80];
        int     year,month;                 /* 年･月 */
        int     number_of_lines;
        int     number_of_pixels;
        short   AT_lines,AT_pixels;         /* 気温 */
        short   AL_lines,AL_pixels;         /* アルベド */
        short   CL_lines,CL_pixels;         /* 日照率 */
        short   PR_lines,PR_pixels;         /* 降水量 */
        short   EL_lines,EL_pixels;         /* 標高 */
	int     i,j;
        double  lat,lon;        /* 緯度・経度 */
        double  albedo;         /* アルベド */
	double  air_tmp;        /* 気温 */
	double  sunshine;       /* 日照率 */
	double  elev;           /* 標高 */
	double  precip;         /* 降水量 */

/* ファイルオープン */

     /* パラメータファイル */

	number_of_lines=360;
	number_of_pixels=720;
        year=2001;

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

     /* 月 */

        fscanf(fp,"%d",&month);
        fprintf(SE,"Month : %2d\n",month);

     /* 気温ﾃﾞｰﾀﾌｧｲﾙ */

	fscanf(fp,"%s %d %d",fname,&AT_lines,&AT_pixels);
        if( (fat=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

        fprintf(SE,"Air temperature  : %12s  ( L:%4d P:%4d )\n",fname,AT_lines,AT_pixels);

     /* アルベド */

	fscanf(fp,"%s %d %d",fname,&AL_lines,&AL_pixels);
        if( (fal=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

        fprintf(SE,"Albedo           : %12s  ( L:%4d P:%4d )\n",fname,AL_lines,AL_pixels);

     /* 日照率 */

	fscanf(fp,"%s %d %d",fname,&CL_lines,&CL_pixels);
        if( (fcl=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

        fprintf(SE,"Cloudiness       : %12s  ( L:%4d P:%4d )\n",fname,CL_lines,CL_pixels);

     /* 降水量 */

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

        fprintf(SE,"Precipitation    : %12s  ( L:%4d P:%4d )\n",fname,CL_lines,CL_pixels);

     /* 標高 */

	fscanf(fp,"%s %d %d",fname,&EL_lines,&EL_pixels);
        if( (fel=open(fname,O_BINARY | O_RDONLY))==FILE_ERR){
                fprintf(SE,"Can't open %s.\n",fname);
                exit(-1);
        }

        fprintf(SE,"Elevation        : %12s  ( L:%4d P:%4d )\n",fname,EL_lines,EL_pixels);

     /* 出力ﾃﾞｰﾀﾌｧｲﾙ */

        fscanf(fp,"%s",fname);
        if( (fo=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);
        }

        fprintf(SE,"Output file name : %s\n",fname);

        fclose(fp);

/* 計算準備 */

        cal_nday(year);

/* グリッドごとに計算 */

        for(i=0;i<number_of_lines;i++){        /************ ライン方向 *****************/
                if( !(i%10) ) fprintf(SE,"*");

            /* 緯度の計算 */

                lat=get_lat(i);

		for(j=0;j<number_of_pixels;j++){  /************ ピクセル方向 ************/

                     /* 経度の計算 */

                        lon=get_lon(j);

                     /* 計算パラメータの取得 */
/*
fprintf(SE,"\n");
fprintf(SE,"L:%3d P:%3d ",i,j);
*/
                        air_tmp=get_temperature(fat,AT_lines,AT_pixels,i,j);   /* 気温 ℃ */
/*
fprintf(SE,"AT:%6.1f ",(float)air_tmp);
*/
			albedo=get_albedo(fal,AL_lines,AL_pixels,i,j);         /* アルベド (double) */
/*
fprintf(SE,"AL:%6.1f ",(float)albedo);
*/
                        if(albedo<=0.0){         /* 海は計算しない */
                                ET[j]=-9999;
                                continue;
                        }

			sunshine=get_cloudiness(fcl,CL_lines,CL_pixels,i,j);   /* 日照率 (double) */
/*
fprintf(SE,"CL:%6.1f ",(float)sunshine);
*/
                        if(sunshine>1.0){        /* 海は計算しない */
                                ET[j]=-9999;
                                continue;
                        }

			elev=get_elevation(fel,EL_lines,EL_pixels,i,j);        /* 標高 (m) */
			precip=get_precipitation(fpr,PR_lines,PR_pixels,i,j);  /* 降水量 (mm) */

                     /* 可能蒸発量の計算 0.1mmで出力*/

			ET[j]=(signed short)(PT_method(month,lat,air_tmp,albedo,sunshine,elev,precip,i,j)*10.0+0.5);

                     /* 結果が負になったらゼロ */

                        if(ET[j]<0)   ET[j]=0;
		}

/* 結果の出力 */

		write(fo,ET,number_of_pixels*sizeof(short));

        }/* next line */

/* ファイルクローズ */

        close(fat);
        close(fal);
        close(fcl);
        close(fel);
        close(fpr);
        close(fo);

}/* end of main */

/*-------------------------------------------------------------------------------------------------
        function PT_method
  -------------------------------------------------------------------------------------------------*/

float   PT_method(
int     month,
double  lat,               /* 緯度：ラジアン */
double  air_tmp,
double  albedo,
double  sunshine,
double  elev,
double  precip,
int     i,
int     j)
{
        double   R;          /* extra terestrial radiation */
	double   Rs;
	double   Rn;
        double   Ep;
        double   delta;      /* Δ */
	double   gamma;      /*    */
        double   ramda;      /* λ */
	double   w0;         /* 日出、日没の時角 */
        double   a,b;
        double   dummy;

/* 大気外日射量(ly/day)->(J/m"/day) -> W/m" */

        dummy=Sun(lat,month,&w0,i,j);             /* 戻り値は各月の可照時間、ここでは w0 のみ使用 */

        R=Qa_calc(month,lat,w0)*4.19E4;       /* ly/day -> J/m"/day */
        R=R/86400.0;                          /* J/m"/day -> W/m" (1W/m"=1J/s/m") */

/* 短波放射量 W/m" */

        a = 0.29 * cos(lat);
        b = 0.71 - a;

        Rs = ( a + b*sunshine ) * R;

/* 正味放射量 W/m" */

        Rn=(1.0-albedo)*Rs - 1.11*(0.2+0.8*sunshine)*(100.0-air_tmp);

/* Δ、γ、λ */

        delta=cal_delta(air_tmp);
        gamma=cal_gamma(elev,air_tmp);
        ramda=cal_ramda(air_tmp);

/* Priestley and Taylor method */

	Ep=1.26*delta/(delta+gamma)*Rn * 86400.0 / ramda  / (1000.0) * (double)nmonth[month];
/*                     W/m"           -> J/day/m" -> g/m" ->  mm/day      ->   mm/month           */
/*
if(i==46 && j==450){
        printf("SUB PET 46-450\n");
        printf("dummy:%f\n",dummy);
        printf("month:%d w0:%f\n",month,w0);
        fprintf(SE,"Lat:%5.1f Ep:%8.3f R:%8.1f Rs:%8.1f Rn:%8.1f\n",lat*R2D,Ep,R,Rs,Rn);
        fprintf(SE,"Sunshine:%6.2f  Air T.:%6.2f  Albedo:%6.2f  El:%8.2f\n",sunshine,air_tmp,albedo,elev);
        fprintf(SE,"ramda:%8.2f  delta:%8.2f  gamma:%8.2f\n",ramda,delta,gamma);
        fprintf(SE,"w0:%f\n",w0);
        fprintf(SE,"\n");
}*/

        return( (float)Ep );
}

/*-------------------------------------------------------------------------------------------------
        function Qa_calc()    戻り値に注意 W/m" or ly/day
  -------------------------------------------------------------------------------------------------*/

double  Qa_calc(
int     mon,            /* 月 */
double  lat,            /* 緯度(radian) */
double  w0)             /* 日出、日没の時角(radian) */
{
        int     day=15;
        double  sek,dummy,ret;

        sek=sekii(mon,day);

/* 直達放射ﾋﾞｰﾑに垂直な平面上の大気外放射量(Rsex)の計算 */

        dummy=(double)(360.0*(double)nday[mon][15]/370.0);
/*
fprintf(SE,"Qa_calc\n");
fprintf(SE,"dummy : %f\n",dummy);
*/

/* RO : 太陽定数 (ly/min) */
/* SC : 太陽定数 (1382 W/m") */

     /* 太陽・地球間の距離 */

        dummy=RO*( 1.0 + 0.033 * cos( D2R*dummy ) ); 
/*      dummy=SC*( 1.0 + 0.033 * cos( D2R*dummy ) ); */

     /* 水平面が受け取る日射量 */

        ret=dummy*( w0*sin(lat)*sin(sek)+cos(lat)*cos(sek)*sin(w0) );

        ret=(1440.0/PI)*ret;     /* ly/day */

        return( ret );
}

/*-------------------------------------------------------------------------------------------------
              function sekii() 赤緯
  -------------------------------------------------------------------------------------------------*/

double  sekii(month,day)
int     month,day;
{
        double x,ret;
        x=(double)nday[month][day];   /* １月０日からの日数 */
        ret=0.01698*(x-80.0);      /* ﾙﾝﾃﾞの式 */
        return( 0.4093*sin( (double)ret ) );      /* ﾙﾝﾃﾞの式 */
}

/*-------------------------------------------------------------------------------------------------
        function cal_nday() 1月0日からの日数
  -------------------------------------------------------------------------------------------------*/

void    cal_nday(year)
int     year;
{
        int     i,j,dummy;

        if( (year%4)==0 ){                           /* 閏年の処理 */
                if((year%100==0) && (year%400!=0)){  /* 但、100で割れるが400で */
                        nmonth[2] = 28;              /* 割れない年は平年 */
                }else{
                        nmonth[2] = 29;
                }
        }

        dummy=0;
        for(i=1;i<=12;i++){
                for(j=1;j<=nmonth[i];j++){
                        dummy++;
                        nday[i][j]=dummy;
                }
        }
}/* end of cal_nday */

/*-------------------------------------------------------------------------------------------------
          function cal_delta() �凾ﾌ計算
  -------------------------------------------------------------------------------------------------*/

double  cal_delta(
double  Temp)                 /* 気温 */
{
        double   dummy1,dummy2;

        dummy1= -1.545232E-08*Temp*Temp*Temp+3.54531E-06*Temp*Temp
                     -6.00776E-04*Temp+0.0726612;
        dummy2= -3.86308E-9*Temp*Temp*Temp*Temp+1.18177E-06*Temp*Temp*Temp
                      -3.00388E-4*Temp*Temp+0.0726612*Temp+1.8095;

        return( dummy1*exp( (double)dummy2 ) );
}

/*-------------------------------------------------------------------------------------------------
        γの計算
  -------------------------------------------------------------------------------------------------*/

double  cal_gamma(
double  elev,
double  air_tmp)
{
	double	ret;

	ret=(1005.0*(1013.25-0.11986*elev+0.000005356*elev*elev))/(622.0*(2501.0-2.37*air_tmp));

	return ( ret );
}

/*-------------------------------------------------------------------------------------------------
        function Sun() 可照時間の計算
  -------------------------------------------------------------------------------------------------*/

double  Sun(lat,mon,w0,i,j)
double  lat;
int     mon;
double  *w0;
int     i;
int     j;
{
        int     day=15;
        double  t1,t2,t0,cos_t0,sin_t0;
        double  sek,ret;
        double  tp; /* 可照時間 */ 

        sek=sekii(mon,day);       /* 各月15日の赤緯 */

        cos_t0=-tan( lat )*tan( sek );
        if(cos_t0<-1.0){
                sin_t0=0.0;   /* チェック！*/
        }else{
                sin_t0=sqrt( (double)( 1.0-cos_t0*cos_t0 ) );
        }

        t0=atan2( (double)sin_t0,(double)cos_t0 );

        t1=-3.81972*t0+12.0;     /* 日出 */
        t2= 3.81972*t0+12.0;     /* 日没 */
        *w0=t0;                  /* 日出、日没の時角 */
        tp=t2-t1;
        if(tp>24.0){
                tp=24.0;
        }

        ret=(tp);
/*
if(i==46 && j==450){
printf("SUN 46-450\n");
fprintf(SE,"Sun : \n");
fprintf(SE,"cos_t0 : %f\n",cos_t0);
fprintf(SE,"sin_t0 : %f\n",sin_t0);
fprintf(SE,"    t0 : %f\n",    t0);
}*/

        return( ret );         /* 可照時間(月間総時間) */
}

/*-------------------------------------------------------------------------------------------------
         function ramda() 蒸発の潜熱(J/g)の計算
  -------------------------------------------------------------------------------------------------*/

double  cal_ramda(Ta)        /* Ta : 気温(ﾟC) */
double  Ta;
{
        return( 2501.0-2.37*Ta );
}

/*-------------------------------------------------------------------------------------------------
         緯度  (ラジアンで返す)
  -------------------------------------------------------------------------------------------------*/

double  get_lat(
int     line)
{
	double	ret;

	ret=(90.0-(double)line*0.5 )*D2R;

        return( ret );
}

/*-------------------------------------------------------------------------------------------------
         経度  (ラジアンで返す)
  -------------------------------------------------------------------------------------------------*/

double  get_lon(
int     pixel)
{
	double	ret;

	ret=( (double)pixel*(double)0.5-(double)180.0 )*(double)D2R;

        return( ret );
}

/*-------------------------------------------------------------------------------------------------
        気温
  -------------------------------------------------------------------------------------------------*/

double  get_temperature(
int     fat,
short   AT_lines,
short   AT_pixels,
short   line,
short   pixel)
{
	long    offset;

	offset=(long)line*(long)AT_pixels*(long)2;
	lseek(fat,offset,SEEK_SET);

	read(fat,buffer2,AT_pixels*2);

	return( (double)buffer2[pixel]/(double)10.0 );
}

/*-------------------------------------------------------------------------------------------------
        アルベド
  -------------------------------------------------------------------------------------------------*/

double  get_albedo(
int     fal,
short   AL_lines,
short   AL_pixels,
short   line,
short   pixel)
{
	long	offset;

	offset=(long)line*(long)AL_pixels;
	lseek(fal,offset,SEEK_SET);

	read(fal,buffer1,AL_pixels);

	return( (double)buffer1[pixel]/100.0 );
}

/*-------------------------------------------------------------------------------------------------
        日照率
  -------------------------------------------------------------------------------------------------*/

double  get_cloudiness(
int     fcl,
short   CL_lines,
short   CL_pixels,
short   line,
short   pixel)
{
	long	offset;

	offset=(long)line*(long)CL_pixels;
	lseek(fcl,offset,SEEK_SET);

	read(fcl,buffer1,CL_pixels);

	return( (double)buffer1[pixel]/100.0 );
}

/*-------------------------------------------------------------------------------------------------
        標高
  -------------------------------------------------------------------------------------------------*/

double  get_elevation(
int     fel,
short   EL_lines,
short   EL_pixels,
short   line,
short   pixel)
{
	long	offset;

	offset=(long)line*(long)EL_pixels*(long)2;
	lseek(fel,offset,SEEK_SET);

	read(fel,buffer2,EL_pixels*2);

	return( (double)buffer2[pixel] );
}

/*-------------------------------------------------------------------------------------------------
        降水量
  -------------------------------------------------------------------------------------------------*/

double  get_precipitation(
int     fpr,
short   PR_lines,
short   PR_pixels,
short   line,
short   pixel)
{
	long	offset;

	offset=(long)line*(long)PR_pixels*(long)2;
	lseek(fpr,offset,SEEK_SET);

	read(fpr,buffer2,PR_pixels*2);

	return( (double)buffer2[pixel] );
}