/****************************************************************************** ＭＯＲＴＯＮ（１９８３）の方法による蒸発散量計算プログラム A.Kondoh 1992 ******************************************************************************/ #include #include #include #include #include #include #include #include #include #include #include #include #define Stefan (5.67032E-08) /* W/m"/K4 */ #define PI 3.141592653 #define D2R (1.745329E-02) /* 度→ﾗｼﾞｱﾝ */ #define R2D 57.29577952 /* ﾗｼﾞｱﾝ→度 */ #define EMISSIVITY 0.92 /* 放射率 */ #define SE stderr #define FILE_ERR -1 #define CRLF puts("") #define CLS printf("\033[2J") #define PAUSE do{}while(!isspace(getch())) /* ﾎﾟｰｽﾞ */ /* ファイルディスクリプタ */ FILE *fi; FILE *fo; /* 入力用変数 */ int year; /* 計算年度 */ float lat; /* 緯度 */ float alt; /* 標高 */ float Pa; /* 年降水量 */ float P_Ps; /* 現地気圧と海面気圧との比 */ float azd; /* 乾季、無雪期、晴天時のアルベドの天頂値 */ int month; /* 月 */ float emissivity; /* 放射率 */ float Tmax,Tmin; /* 最高・最低気温 */ float Ta; /* 平均気温 */ float Rh; /* 相対湿度 */ float sd; /* 日照時間 */ float Temp[13]; /* 最高気温と最低気温の平均（℃）*/ float Td[13]; /* 露点温度（℃）*/ float Sr[13]; /* 日照率 */ /* 出力用変数 */ float Etp[13]; /* 可能蒸発散量 */ float Rtp[13]; /* 平衡温度における土壌−植生表面の正味放射量 */ float Etw[13]; /* 湿潤環境の地域蒸発散量 */ float Et[13]; /* 補間関係による地域蒸発散量 */ static int nmonth[13]={0,31,28,31,30,31,30,31,31,30,31,30,31}; int nday[13][32]; void morton(); void get_input_parameters(); float Esa(); /* 飽和水蒸気圧 */ float max_sunshine(); /* 可照時間 */ float get_Tdp(); /* 露点温度の計算 */ void cal_nday(); void output_results(); float sekii(); /* 太陽の赤緯 */ int Ta_switch; /* 0:Tmax,Tmin / 1:Tmean */ /*----------------------------------------------------------------------------- ﾓｰﾄﾝ法ﾒｲﾝﾌﾟﾛｸﾞﾗﾑ -----------------------------------------------------------------------------*/ void main(void) { /* タイトル */ CLS; fprintf(SE,"====== Areal Evapotranspiration by Morton(1983) ======\n"); /* データ入力 */ get_input_parameters(); /* ｽﾃｰｼｮﾝの緯度（度）・標高（ｍ）・年降水量（ｍｍ）*/ for(month=1;month<=12;month++){ P_Ps=pow( ((288.0-0.0065*alt)/288.0) , 5.256 ); azd=0.26-0.00012*Pa*pow(P_Ps,0.5)*(1.0+fabs(lat/42.0)+pow(lat/42.0,2.0)); if(azd<0.11) azd=0.11; /* ほとんどの地表面は0.11になる */ if(azd>0.17) azd=0.17; /* 砂漠 */ /* ﾓｰﾄﾝ法の計算 */ morton(); }/* next month */ /* 計算結果出力 */ output_results(); }/* end of main */ /*----------------------------------------------------------------------------- MORTON(1983)の方法による地域蒸発散量の計算 -----------------------------------------------------------------------------*/ void morton() { float vd; /* 露点温度における飽和水蒸気圧(mbar) */ float vt; /* 気温における飽和水蒸気圧 */ float delta; /* 気温における飽和水蒸気圧(mbar) */ float theta; /* 太陽の赤緯（度）*/ float cosZ; /* 正午における太陽の天頂角距離 */ float Za; /* 〃 */ float cosW; /* 日の出から正午までの時角 */ float Wd; /* 〃 */ float cosz; /* 太陽の平均天頂角距離 */ float rv; /* 太陽の位置ﾍﾞｸﾄﾙ radius vector */ float Ge; /* 大気圏外における全天日射量 W/m" */ float azz; /* 雪のない晴天時のｱﾙﾍﾞﾄﾞの天頂値 */ float c0; float az; /* 晴天ｱﾙﾍﾞﾄﾞの天頂値 */ float a0; /* 晴天ｱﾙﾍﾞﾄﾞ */ float Wp; /* 可降水量 */ float c1; float tc; /* turbidity coefficient */ float tau; /* 直達日射の晴天時の透過率 */ float tau_a; /* τの吸収による成分 */ float G0; /* 晴天時の全天日射量 */ float Gi; /* 直達日射量 */ float c2; float albedo; /* 平均ｱﾙﾍﾞﾄﾞ */ float rho; /* 雲による大気放射量の比例増加 */ float Rt; /* 気温における土壌−植生表面の正味放射量 */ float Rtc; float fz; float rPs; float b0; float rp; float sf; /* the stability factor の逆数 */ float fT; /* 水蒸気輸送係数 */ float ramda; /* 熱輸送係数 */ float Tp; float vp; float delta_p; float dTp; float b1; float b2; float Lw; float lh; /* 潜熱 */ float alpha,beta; /* Morton(1983) Step 2 */ /* 露点温度における飽和水蒸気圧(mbar) */ vd=6.11*exp( 17.27*Td[month]/(Td[month]+237.3)); /* 気温における飽和水蒸気圧(mbar) */ if(Temp[month]>=0.0){ alpha=17.27; beta=237.3; }else{ alpha=21.88; beta=265.5; } vt=6.11*exp( alpha*Temp[month]/(Temp[month]+beta) ); /* 気温における飽和水蒸気圧曲線の勾配 */ delta=alpha*beta*vt/pow( (Temp[month]+beta),2.0 ); /* Morton(1983) Step 3 大気圏外における日射 Ge(W/m") の計算 */ /* 太陽の赤緯 */ theta=23.2*sin( D2R*( 29.5*(float)month - 94.0 ) ); /* 正午における太陽の天頂角距離 */ cosZ=cos( D2R*( lat-theta ) ); if(cosZ<0.001) cosZ=0.001; /* 太陽の昇らない極域 */ Za= lat - theta; /* degree */ /* 日の出から正午までの時角 W */ cosW=1.0 - cosZ/( cos(D2R*lat)*cos(D2R*theta) ); if(cosW<-1.0) cosW=-1.0; /* 太陽の沈まない極域 */ Wd=R2D*acos(cosW); /* degree */ /* 太陽の平均天頂角距離 */ cosz=cosZ+( (180.0/PI)*sin(D2R*Wd)/Wd-1.0 )*cos(D2R*lat)*cos(D2R*theta); /* 太陽の位置ﾍﾞｸﾄﾙ(radius vector) */ rv=1.0+(1.0/60.0)*sin( D2R*( 29.5*(float)month -106.0 )); /* 大気圏外における全天放射量 (W/m") */ Ge=(1354.0/(rv*rv))*(Wd/180.0)*cosz; /* Morton(1983) step 4 */ /* 雪のない晴天時のｱﾙﾍﾞﾄﾞの天頂値 */ azz=azd; if(azz<0.11) azz=0.11; if(azz>0.5*(0.91-vd/vt)) azz=0.5*(0.91-vd/vt); /* 晴天時のｱﾙﾍﾞﾄﾞの天頂値 */ c0=vt-vd; if(c0<0.0) c0=0.0; if(c0>1.0) c0=1.0; az=azz+(1.0-c0*c0)*(0.34-azz); /* 晴天時のｱﾙﾍﾞﾄﾞ */ a0=az*( exp(1.08)-(2.16*cosZ/PI+sin(D2R*Za) )*exp(0.012*Za)); a0=a0/( 1.473*(1.0-sin(D2R*Za)) ); /* Morton(1983) step 5 */ /* 可降水量 (mm) */ Wp=vd/(0.49+Temp[month]/129.0); /* turbidity coefficient */ c1=21.0-Temp[month]; if(c1<0.0) c1=0.0; if(c1>5.0) c1=5.0; tc=(0.5+2.5*cosz*cosz)*exp(c1*(P_Ps-1.0)); /* Morton(1983) step 6 τ:直達日射の晴天時の放射率 */ tau=exp( -0.089*pow( (P_Ps/cosz),0.75 ) -0.083*pow( (tc/cosz) , 0.90 ) -0.029*pow( (Wp/cosz) , 0.60 ) ); /* Morton(1983) step 7 吸収によるτの成分 */ tau_a=exp( -0.0415*pow( (tc/cosz),0.90 ) -pow(0.0029,0.5)*pow( (Wp/cosz),0.3 )); if(tau_a1.0) c2=1.0; rho=0.18*( (1.0-c2)*pow( (1.0-Sr[month]),2.0 )+c2*pow( (1.0-Sr[month]),0.5 ) )/P_Ps; /* Morton(1983) step 11 気温における土壌−植生表面の長波出放射 */ Lw=emissivity*Stefan*pow( (Temp[month]+273.0),4.0 )*( 1.0-(0.71+0.007*vd*P_Ps)*(1.0+rho) ); if( Lw<0.05*emissivity*Stefan*pow( (Temp[month]+273.0),4.0 ) ){ Lw=0.05*emissivity*Stefan*pow( (Temp[month]+273.0),4.0 ); } /* Morton(1983) step 12 */ /* 気温における土壌−植生表面の正味放射量 */ Rt=(1.0-albedo)*Gi-Lw; Rtc=Rt; if(Rtc<0.0) Rtc=0.0; if(Temp[month]>=0.0){ fz=28.0; /* W/m" */ rPs=0.66; /* mbar/ﾟC */ }else{ fz=28.0*1.15; rPs=0.66/1.15; } b0=1.00; /* CRAE model */ rp=rPs*P_Ps; /* sf=1/ζ*/ sf=0.28*(1.0+vd/vt)+delta*Rtc/( rp*pow( (1.0/P_Ps),0.5)*b0*fz*(vt-vd) ); if(sf>1.0) sf=1.0; fT=pow( (1.0/P_Ps),0.5)*fz*sf; ramda=rp+4.0*emissivity*Stefan*pow( (Temp[month]+273.0),3.0 )/fT; /* Morton(1983) step 13 */ Tp=Temp[month]; vp=vt; delta_p=delta; dTp=99999.9; while(dTp>0.01){ dTp=( Rt/fT+vd-vp+ramda*(Temp[month]-Tp) )/(delta_p+ramda); Tp=Tp+dTp; vp=6.11*exp( alpha*Tp/(Tp+beta) ); delta_p=alpha*beta*vp/pow( (Tp+beta),2.0 ); } /* Morton(1983) step 14 */ Etp[month]=Rt - ramda*fT*(Tp-Temp[month]); Rtp[month]=Etp[month] + rp*fT*(Tp-Temp[month]); b1=14.0; /* W/m" */ b2=1.20; Etw[month]=b1+b2/(1.0+rp/delta_p)*Rtp[month]; if(Etw[month]<0.5*Etp[month]) Etw[month]=0.5*Etp[month]; if(Etw[month]>Etp[month]) Etw[month]=Etp[month]; /* Morton(1983) step 15 補間関係 */ Et[month]=2.0*Etw[month]-Etp[month]; /* Morton(1983) step 16 単位の変換 (W/m"->mm/month) */ if(Temp[month]>=0.0){ lh=28.5; /* 蒸発の潜熱 W*days/kg */ }else{ lh=28.5*1.15; /* 昇華の潜熱 */ } /* 正味放射量 */ Rtp[month] = Rtp[month]/lh*(float)nmonth[month]; /* 可能蒸発散量 */ Etp[month] = Etp[month]/lh * (float)nmonth[month]; /* 湿潤環境の蒸発散量 */ Etw[month] = Etw[month]/lh * (float)nmonth[month]; /* 地域蒸発散量 */ Et[month] = Et[month]/lh * (float)nmonth[month]; }/* end of morton */ /*----------------------------------------------------------------------------- 入力パラメーター -----------------------------------------------------------------------------*/ void get_input_parameters() { int i; char fname[80]; CRLF; fprintf(SE,">> input file name : "); scanf("%s",fname); if((fi=fopen(fname,"rt"))==NULL){ fprintf(SE,"Can't open %s.\n",fname); exit(-1); } /* fprintf(SE,">> emissivity : "); scanf("%f",&emissivity); */ emissivity=EMISSIVITY; fprintf(SE,">> Ta (Tmax,Tmin:0/Tmean:1) : ",Ta_switch); scanf("%d",&Ta_switch); /* ｽﾃｰｼｮﾝﾃﾞｰﾀ入力 */ fscanf(fi,"%d %f %f %f",&year,&lat,&alt,&Pa); CRLF; fprintf(SE,"Year : %6d\n",year); fprintf(SE,"Latitude : %6.1f (degree)\n",lat); fprintf(SE,"Altiitude a.s.l. : %6.1f (m)\n",alt); fprintf(SE,"Annual precipitation : %6.1f (mm)\n",Pa); /* 気象ﾃﾞｰﾀ入力 */ for(i=1;i<=12;i++){ switch(Ta_switch){ case 0: fscanf(fi,"%f %f %f %f",&Tmax,&Tmin,&Rh,&sd); break; case 1: fscanf(fi,"%f %f %f",&Ta,&Rh,&sd); } Rh=Rh/100.0; /* % -> 小数 */ /* 気温は最大と最小気温の平均 */ switch(Ta_switch){ case 0: Temp[i]=(Tmax+Tmin)/2.0; break; case 1: Temp[i]=Ta; } /* 露点温度 */ Td[i]=get_Tdp(i); /* 日照率 */ Sr[i]=sd/max_sunshine(i); }/* next month */ /* 入力パラメータ確認 */ CRLF; fprintf(SE,"------T-------Td-------Sr-----\n"); for(i=1;i<=12;i++){ fprintf(SE," %6.1f %6.1f %6.3f\n",Temp[i],Td[i],Sr[i]); } fprintf(SE,"-----ﾟC-------ﾟC--------------\n"); CRLF; PAUSE; cal_nday(year); }/* end of get_input_parameters */ /*----------------------------------------------------------------------------- 露点温度の計算 (農業気象の測器と測定法 p.74) -----------------------------------------------------------------------------*/ float get_Tdp( int mon) /* 月 */ { float e; /* 水蒸気圧 */ float t1,t2,ft,ftd; float a0,a1,a2,a3,a4; float tol; /* 繰返し計算の許容誤差 */ int count; float At; /* 気温 */ At=Temp[mon]; tol=0.01; a0=1.809378; /* 岡田の式の係数 */ a1=0.07266115; a2=-3.003879E-4; a3=1.181765E-6; a4=-3.863083E-9; e=Rh*Esa(mon); /* 水蒸気圧 mbar */ t1=At; /* 繰返し計算による露点温度の初期値 */ /* 許容誤差がtol未満になるか30回計算したら終了 */ count=0; while( (count++)<30){ ft =exp(a0+a1*t1+a2*t1*t1+a3*t1*t1*t1+a4*t1*t1*t1*t1)-e; ftd=(a1+2.0*a2*t1+3.0*a3*t1*t1+4.0*a4*t1*t1*t1)*(ft+e); t2=t1-ft/ftd; if( fabs(t1-t2)=0.0){ ret=( -3.86308E-09*At + 1.18177E-06)*At; ret=((ret-3.00388E-04)*At+0.0726612)*At+1.8095; ret=exp( ret ); }else{ ret=((1.36277E-06*At-2.99091E-04)*At+0.0823895)*At; ret=exp(ret+1.80938); } return( ret ); } /*------------------------------------------------------------------- 可照時間の計算 -------------------------------------------------------------------*/ float max_sunshine(mon) int mon; { int day=15; float t1,t2,t0,cos_t0,sin_t0; float sek,ret; sek=D2R*sekii(mon,day); /* 各月15日の赤緯 */ cos_t0=-tan( D2R*lat )*tan( (double)sek ); sin_t0=sqrt( 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; /* 日没 */ ret=(t2-t1); return( ret ); /* 可照時間 */ } /*----------------------------------------------------------------------------- 計算結果の表示・出力 -----------------------------------------------------------------------------*/ void output_results() { int i; float Total_Etp; float Total_Rtp; float Total_Etw; float Total_Et; CLS; fprintf(SE,"====== Areal Evapotranspiration by Morton(1983) ======\n"); CRLF; fprintf(SE," Year : %6d\n",year); fprintf(SE," Latitude : %6.1f (degree)\n",lat); fprintf(SE," Altiitude a.s.l. : %6.1f (m)\n",alt); fprintf(SE," Annual precipitation : %6.1f (mm)\n",Pa); CRLF; Total_Etp=0.0; Total_Rtp=0.0; Total_Etw=0.0; Total_Et=0.0; fprintf(SE," Month----Etp------Rtp------Etw-------Et---\n"); for(i=1;i<=12;i++){ fprintf(SE," M:%3d %6.1f %6.1f %6.1f %6.1f\n" ,i,Etp[i],Rtp[i],Etw[i],Et[i]); Total_Etp += Etp[i]; Total_Rtp += Rtp[i]; Total_Etw += Etw[i]; Total_Et += Et[i]; } fprintf(SE," ------------------------------(mm/month)--\n"); fprintf(SE," Total %6.1f %6.1f %6.1f %6.1f\n" ,Total_Etp,Total_Rtp,Total_Etw,Total_Et); fprintf(SE," ------------------------------(mm/month)--\n"); PAUSE; }/* end of output */ /*------------------------------------------------------------------- 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 sekii() 赤緯(度) -------------------------------------------------------------------*/ float sekii(month,day) int month,day; { float x,ret; x=(float)nday[month][day]; /* １月０日からの日数 */ ret=23.45*sin( D2R*( (x-80.0)/370.0*360.0 ) ); /* ﾙﾝﾃﾞの式 */ return( ret ); }