/******************************************************************** ペンマン法による月蒸発散量の計算 A.Kondoh 1988-1991 ********************************************************************/ #include #include #include #include #include #include #include #include #include #include #include #define Stefan 8.132E-11 /* ly･min-1･ﾟK-4 */ #define Hg 0.7500616 /* mb -> mmHg */ #define RO 1.96 /* 太陽定数(ly/min,cal/cm2min) */ #define PI 3.141592653 #define D2P 1.745329E-02 /* 度→ﾗｼﾞｱﾝ */ #define SE stderr #define CRLF puts("") #define PAUSE do{}while(!isspace(getch())) /* ﾎﾟｰｽﾞ */ float ET[13]; /* Evapotranspiration (mm/day) */ float Ea[13]; /* Aero dynamic term */ float Ta[13]; /* Air Temperature (ﾟC) */ float U2[13]; /* Wind Speed at 2m (m/s) */ float HM[13]; /* Relative Humidity (%) */ float SD[13]; /* Sunshine Duration (hour) */ float Short_R[13]; /* Short Radiation */ float Long_R[13]; /* Long Radiation */ float Rn[13]; /* Net Radiation(J/m"/day) */ float Qn[13]; /* Effective Energy (mm/day) */ int nday[13][32]; float lat; /* 緯度 */ float albedo; /* アルベド */ int loc_num; /* 地点番号 */ char loc_name[10]; /* 地点名 */ char loc_name1[10]; /* 漢字 */ char loc_name2[10]; /* ローマ字 */ char loc_name3[10]; /* カタカナ */ float lat; /* 緯度 */ float lon; /* 経度 */ float alt; /* 標高 */ float factor; /* 実蒸発散量への変換係数 */ float sekii(); /* 赤緯 */ static int nmonth[13]={0,31,28,31,30,31,30,31,31,30,31,30,31}; int lack[13]; /* 0:ﾃﾞｰﾀあり、1:欠測 */ FILE *fi,*fo,*fd; /* ﾌｧｲﾙﾃﾞｨｽｸﾘﾌﾟﾀ */ /*------------------------------------------------------------------- function main() -------------------------------------------------------------------*/ void main(argc,argv) int argc; char *argv[]; { int mon; int method_switch; void read_data(),write_data(); /* ﾃﾞｰﾀ入出力 */ float Penman(); /* ﾍﾟﾝﾏﾝ法 空力項はDunne & Leopold(1978) */ float Fitz_Stern(); /* Fitzpatrick & Stern(1965) */ int loc_num; /* 地点番号 */ char loc_name3[9]; /* 地点名（カナ)*/ int year,month; /* 年･月 */ /* スタート */ CRLF; fprintf(SE,"<<< Calculation of Evapotranspiration by PENMAN Method >>>\n"); CRLF; fprintf(SE," USAGE: PENMAN [input file] [output file] [index file] [method] [albedo] [f]\n"); CRLF; fprintf(SE," METHOD: Penman ... 1\n"); fprintf(SE," Fitzpatrick & Stern(1965) ... 2\n"); CRLF; fprintf(SE," Albedo : 0-100 integer(percent)\n"); fprintf(SE," f : 0-100 integer(percent)\n"); CRLF; if(argc!=7) exit(-1); /* パラメータ処理 */ method_switch=atoi(argv[4]); /* 手法選択ｽｲｯﾁ */ albedo=(float)atoi(argv[5])/100.0; /* ｱﾙﾍﾞﾄﾞ */ factor=(float)atoi(argv[6])/100.0; /* 実蒸発散量への変換係数 */ CRLF; switch(method_switch){ case 1: fprintf(SE,"Penman method selected\n"); break; case 2: fprintf(SE,"Fitzpatrick & Stern(1965) is selected for Rn.\n"); } fprintf(SE,"Albedo : %f\n",albedo); fprintf(SE,"f : %f\n",factor); /* ファイルオープン */ /* 入力ﾃﾞｰﾀﾌｧｲﾙ */ if( (fi=fopen(argv[1],"rt"))==NULL){ fprintf(SE,"Can't open %s.\n",argv[1]); exit(-1); } /* 出力ﾃﾞｰﾀﾌｧｲﾙ */ if( (fo=fopen(argv[2],"wt"))==NULL){ fprintf(SE,"Can't open %s.\n",argv[2]); exit(-1); } /* ｲﾝﾃﾞｯｸｽﾃﾞｰﾀﾌｧｲﾙ （fd=NULL の場合は計算値のみ出力）*/ if( (fd=fopen(argv[3],"rt"))==NULL){ fprintf(SE,"Can't open %s.\n",argv[3]); } /* 計算 */ while(!feof(fi)){ /* パラメータ入力 */ read_data(&loc_num,loc_name3,&year,&month); fprintf(SE,"Location No: %3d [%8s] Y:%4d\n",loc_num,loc_name3,year); /* 月ごとに計算 */ for(mon=1;mon<=12;mon++){ /* ﾃﾞｰﾀ欠測の場合は計算をﾊﾟｽ */ if(lack[mon]==1){ ET[mon]=-999.9; continue; } /* 可能蒸発散量の計算 */ switch(method_switch){ case 1: ET[mon]=Penman(mon); break; case 2: ET[mon]=Fitz_Stern(mon); break; default: fprintf(SE,"ERR!!\n"); exit(-1); } }/* next month */ /* 結果の出力 */ write_data(loc_num,loc_name3,year,month); }/* wend */ fclose(fi); fclose(fo); fclose(fd); }/* end of main */ /*------------------------------------------------------------------- function Penman() -------------------------------------------------------------------*/ float Penman(mon) int mon; { float sd; /* sunshine duration */ float esa,ea; /* vapor pressure(mmHg) */ float Qa; /* extra terestrial radiation */ float s; /* 日照率 */ float ramda; /* latent heat for evaporation */ float Eo,dl,gm,T4,w0; float Sun(),Esa(),Qa_calc(); float cal_delta(),cal_ramda(); /* estimation of net radiation */ sd=Sun(mon,&w0); /* 各月の可照時間 */ esa=Esa(mon); /* 飽和水蒸気圧(mb) */ ea=esa*HM[mon]; s=SD[mon]/sd; /* 日照率 */ Qa=Qa_calc(mon,w0)*4.19E4; /* 大気外日射量(ly/day)->(J/m"/day) */ ramda=cal_ramda(mon); T4=pow( (double)(273.15+Ta[mon]),(double)4); /* Penman I method */ Short_R[mon]=(1.0-albedo)*Qa*( 0.18+0.55*s );/* J/m"/day */ Long_R[mon]=Stefan*T4*( 0.56-0.092*sqrt( ea*Hg ))*(0.1+0.9*s); Long_R[mon]=Long_R[mon]*6.034E7; /* ly/min -> J/m"/day */ Rn[mon]=Short_R[mon]-Long_R[mon]; /* 単位 J/m"/day */ dl=cal_delta(mon); gm=0.66; /* 1000mb,20ﾟCの時の乾湿計定数 */ Ea[mon]=(0.013+0.00016*U2[mon])*(esa-ea); /* Dunne & Leopold,1978),U2:km/day */ Ea[mon]=Ea[mon]*10.0; /* cm/day -> mm/day */ Qn[mon]=Rn[mon]/(1.0E03*ramda); /* mm/day */ Eo=(Qn[mon]*dl+gm*Ea[mon])/(dl+gm); return( Eo ); } /*------------------------------------------------------------------- function Qa_calc() -------------------------------------------------------------------*/ float Qa_calc(mon,w0) int mon; float w0; { int day=15; float sek,dummy,ret; sek=sekii(mon,day); /* 直達放射ﾋﾞｰﾑに垂直な平面上の大気外放射量(Rsex)の計算 */ dummy=(double)(360.0*(float)nday[mon][15]/370.0); dummy=RO*( 1.0 + 0.033 * cos( D2P*dummy ) ); /* RO : 太陽定数 (ly/min) */ ret=dummy*( w0*sin(lat)*sin(sek)+cos(lat)*cos(sek)*sin(w0) ); ret=(1440.0/PI)*ret; return( ret ); } /*------------------------------------------------------------------- function Sun() 可照時間の計算 -------------------------------------------------------------------*/ float Sun(mon,w0) int mon; float *w0; { int day=15; float t1,t2,t0,cos_t0,sin_t0; float sek,ret; sek=sekii(mon,day); /* 各月15日の赤緯 */ cos_t0=-tan( (double)lat )*tan( (double)sek ); 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; ret=(t2-t1); return( ret ); /* 可照時間(月間総時間) */ } /*------------------------------------------------------------------- function sekii() 赤緯 -------------------------------------------------------------------*/ float sekii(month,day) int month,day; { float x,ret; x=(float)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 Esa() 気温に対する飽和水蒸気圧 -------------------------------------------------------------------*/ float Esa(mon) int mon; { float Temp,ret; Temp=Ta[mon]; ret=(-3.86308E-09*Temp+1.18177E-06)*Temp; ret=((ret-3.00388E-04)*Temp+0.0726612)*Temp+1.8095; return( (float)exp( (double)ret ) ); } /*------------------------------------------------------------------- function cal_delta() 凾ﾌ計算 -------------------------------------------------------------------*/ float cal_delta(mon) int mon; { float Temp,dummy1,dummy2; Temp=Ta[mon]; 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 ) ); } /*------------------------------------------------------------------- function cal_ramda() 蒸発の潜熱(J/g)の計算 -------------------------------------------------------------------*/ float cal_ramda(mon) /* Ta : 気温(ﾟC) */ int mon; { return( 2501.0-2.37*Ta[mon] ); } /*------------------------------------------------------------------- function Fitz_Stern(mon) -------------------------------------------------------------------*/ float Fitz_Stern(mon) int mon; { float sd; /* sunshine duration */ float esa,ea; /* vapor pressure(mmHg) */ float Qa; /* extra terestrial radiation */ float s; /* 日照率 */ float ramda; /* latent heat for evaporation */ float Eo,dl,gm,T4,w0; float Sun(),Esa(),Qa_calc(); float cal_delta(),cal_ramda(); /* estimation of net radiation */ sd=Sun(mon,&w0); /* 各月の可照時間 */ esa=Esa(mon); /* 飽和水蒸気圧(mb) */ ea=esa*HM[mon]; s=SD[mon]/sd; Qa=Qa_calc(mon,w0)*4.19E4; /* 大気外日射量(ly/day)->(J/m"/day) */ ramda=cal_ramda(mon); T4=pow( (double)(273.15+Ta[mon]),(double)4); /* Fitzpatrick & Stern(1965) */ /* short radiation J/m"/day */ Short_R[mon]=(1.0-albedo)*Qa*( (0.385*s+0.373)-0.0042/(s+0.0154) ); /* long radiation */ Long_R[mon]=Stefan*T4*( 0.352-0.049*sqrt( ea*Hg ))*(0.3+0.7*s); Long_R[mon]=Long_R[mon]*6.034E7; /* ly/min -> J/m"/day */ /* net radiation */ Rn[mon]=Short_R[mon]-Long_R[mon]; /* 単位 J/m"/day */ dl=cal_delta(mon); gm=0.66; /* 1000mb,20ﾟCの時の乾湿計定数 */ /* Aerodynamic term */ Ea[mon]=(0.013+0.00016*U2[mon])*(esa-ea); /* Dunne & Leopold,1978) */ Ea[mon]=Ea[mon]*10.0; /* cm/day -> mm/day */ Qn[mon]=Rn[mon]/(1.0E03*ramda); /* mm/day */ Eo=(Qn[mon]*dl+gm*Ea[mon])/(dl+gm); return( Eo ); } /*-------------------------------------------------------------------- function read_data() --------------------------------------------------------------------*/ void read_data( int *loc_num, char loc_name3[], int *year, int *month) { int i; int cur_loc_num; void cal_nday(); float dm; /* ダミー変数 */ /* SDPデータの入力 */ for(i=1;i<=12;i++){ fscanf(fi,"%d %s %s %s %d %f %f %f", loc_num,loc_name1,loc_name2,loc_name3,month,&lat,&lon,&alt); /* 入力ﾃﾞ-ﾀは12ヵ月づつ並んでいることを仮定 */ if(i!=1 && (*loc_num)!=cur_loc_num){ fprintf(SE,"ERROR reading data!\n"); exit(-1); } cur_loc_num=(*loc_num); /* データ入力 */ fscanf(fi,"%f %f %f %f %f %f %f %f %f %f" ,&dm,&dm,&dm,&Ta[i],&dm,&dm,&dm,&HM[i],&dm,&U2[i]); fscanf(fi,"%f %f %f %f %f %f %f %f %f %f" ,&dm,&dm,&dm,&dm,&dm,&SD[i],&dm,&dm,&dm,&dm); fscanf(fi,"%f %f %f %f %f %f %f %f %f %f %f", ,&dm,&dm,&dm,&dm,&dm,&dm,&dm,&dm,&dm,&dm,&dm); /* 欠測処理 */ if(Ta[i]<-900.0 || U2[i]<-900.0 || HM[i]<-900.0 || SD[i]<-900.0){ lack[i]=1; }else{ lack[i]=0; } /* 計算のため、単位変更 */ U2[i]=U2[i]*86.4; /* m/s -> km/day */ HM[i]=HM[i]/100.0; /* % -> rate */ SD[i]=SD[i]*(float)nmonth[i]; /* hour/day -> hout/month */ } /* cal_nday(*year);*/ /* 月の日数の計算 */ } /*------------------------------------------------------------------- function write_data() -------------------------------------------------------------------*/ void write_data( int loc_num, char loc_name3[9], int year, int month) { int i; float total=0.0; int loc_num_idx; float etm; /* 月蒸発散量 */ float lat,lon; float alt; char loc_name1[13]; char loc_name2[19]; int lack_on; /* 欠測有無判定スイッチ */ if(fd!=NULL){ rewind(fd); while(1){ fscanf(fd,"%d %s %s %f %f %f", &loc_num_idx,loc_name1,loc_name2,&lat,&lon,&alt); if(loc_num_idx==loc_num){ break; } if(feof(fd)){ fprintf(SE,"ERROR: END OF INDEX FILE\n"); lat=-999.9; lon=-999.9; alt=-999.9; break; } }/* wend */ }else{ lat=-999.9; lon=-999.9; alt=-999.9; }/* end if */ fprintf(fo,"%3d %8s %10.5f %10.5f %6.1f %4d ",loc_num,loc_name3,lat,lon,alt,year); total=0.0; lack_on=0; for(i=1;i<=12;i++){ if(lack[i]==0){ etm=factor*ET[i]*(float)nmonth[i]; /* 日蒸発散量→月蒸発散量 */ }else{ etm=-999.9; /* 欠測の場合 */ lack_on=1; } fprintf(fo,"%6.1f ",etm); total+=etm; } if(lack_on==0){ fprintf(fo,"%6.1f\n",total); }else{ /* 欠測の場合は年値を出力しない */ total=-999.9; fprintf(fo,"%6.1f\n",total); } }