토양 증발산량(ETO) 계산 JAVA 버전
ETO - 토양 증발산량
토양 증발산량이란 토양에서 증발되는 계수를 의미한다.
가뭄과 관련된 직종 또는 직업에서 토양 증발산량을 필요로한다.
아래 내용은 토양 증발산량의 대표적인 측정법인
FAO Penman-Monteith 모형을 기준으로 JAVA로 작성한 내용이다.
** 1시간 기준으로 토양 증발산량을 계산한 내용
** j 값은 날짜를 365일변 변경한 값
public static DateTimeFormatter fm_day = DateTimeFormatter.ofPattern("D");
** t 값은 현재 시간을 나타냄
public static DateTimeFormatter fm_hour = DateTimeFormatter.ofPattern("HH");
public class EtoCalcService {
//DEM, LAT, LONG 값은 지정된 장소에서 측정한다면 상수값으로 사용할 수있다.
static final int DEM = 45;
static final float LAT = 36.22f;
static final float LONG = 127.21f;
public double pmet(double j ,double Thr ,double RHhr ,double t ,double Rss ,double u2 ,double Pa)
{
/*
j 오늘 올해에 몇일(day)
DEM 해발(m) 현장 45m
lat 위도(현장 위도36.22)
lon 경도(현장 경도127.21) 공식계산할 때 필요없음, 이미 공식에서 실제 현장 데이터를 입력했음
Thr 한시간에 평균 온도(℃) 센서데이터
RHhr 한시간에 평균 상대습도(%) 센서데이터
t 현제시간(예:새벽2시부터 3시사이의 eto를 계산하면 2.5로 입력,오후14시부터 15시사이면 14.5로 입력)
Rss 일사량(w/m2) 센서데이터
u2 2미터에 풍속(m/s) 센서데이터
Pa 실제 센서데이터（hPa）
*/
//γ,psychrometric constant [kPa °C-1]
// P=101.3*((293-0.0065*DEM)/293)**5.26 #해발 데이터로 계산한 기압값이
//hPa는 Kpa로 환산
double P = Pa*10;
double gama = (double)(0.665/1000*P);
//e0-ea（kPa）
//eo,saturation vapour pressure at air temperature Thr [kPa]
double e0 = (double)(0.6108*Math.exp(17.27*Thr/(237.3+Thr)));
//ea,average hourly actual vapour pressure [kPa],상대습도로 계산법,공식54
double ea = (double)e0*RHhr/100;
double VPD = (double)e0-ea;
//saturation slope vapour pressure curve at ThrΔ(kPa℃-1)
double deta = 4098*0.6108*Math.exp((17.27*Thr)/(Thr+237.3))/Math.pow((Thr+237.3),2);
//Rn，net radiation at the grass surface(MJ m-2/hour)
//dr，extraterrestrial radiation: correction due to distance Earth-Sun,공식23
double dr = 1+0.033*Math.cos(2*Math.PI/365*j);
//δ,공식24
double delta = 0.409*Math.sin(2*Math.PI/365*j-1.39);
// ϕ,conversion from decimal degrees(latitude) to radians,공식22
double fi = Math.PI/180*LAT;
//omegas = math.acos(-math.tan(fi)*math.tan(delta)) #ws,공식25,day,month,year로 계산
//#공시33
double b=2*Math.PI*(j-81)/364;
//공식32
double Sc=0.1645*Math.sin(2*b)-0.1255*Math.cos(b)-0.025*Math.sin(b);
//공식31
double omega=Math.PI/12*((t+0.06667*(120-127.21)+Sc)-12);
//공식29,한시간당,공식중t1은 1로
double omega1=omega-(Math.PI/24);
//공식30,한시간당,공식중t1은 1로
double omega2=omega+(Math.PI/24);
//공식 25
double ws = Math.acos(-Math.tan(fi) * Math.tan(delta));
// N = 24 / Math.PI * ws-daylight hours공식34 실제 일사량 측정할 수 있어서 필요없음
// Extraterrestrial radiation for hourly or shorter periods (Ra)
double Ra = 0.0;
if (omega < -ws) {
Ra = 0;
}else
{
//Ra(MJm-2day-1)공식28
Ra = 12*60/Math.PI*0.082*dr*((omega2-omega1)*Math.sin(fi)*Math.sin(delta)+Math.cos(fi)*Math.cos(delta)*(Math.sin(omega2)-Math.sin(omega1)));
}
//실제 일사량 데이터 w/m-2,MJm-2로 환산
double Rs = Rss*0.0036;
//Clear-sky solar radiation (Rso)(MJm-2/hour) #공식37
double Rso = (0.75+(2*Math.pow(10,(-5)))*DEM)*Ra;
//Clear-sky solar radiation (Rso )(MJm-2/hour)공식38
double Rns = (1-0.23)*Rs;
double aa = 0.0;
if (((ws-0.79)<=omega) && (omega<=(ws-0.52)))
{
//현장 기후는 subhumid climates
aa=0.6;
} else {
aa = (Rs / Rso);
}
//#Rnl(MJm-2/hour),공식39
double Rnl = 0.000000004903*(Math.pow((Thr+273.16),4))*(0.34-0.14*Math.pow(ea,0.5))*(1.35*aa-0.35);
//#공식40
double Rn = Rns-Rnl;
//u2 = u10*4.87/(math.log(67.8*10-5.42)) - (m/s)공식에서 2미터에 측정하는 풍속, 실제현장 센서 거의 2미터에 실치해서 환산 필요없음
// G soil heat flux（MJm-2/hour）, daytime Rn>0, nighttime Rn<0,공식45,46
double G = 0.0;
if (Rn<0){
G=0.5*Rn;
}else
{
G=0.1*Rn;
}
//#단위：mm
double eto = (0.408*deta*(Rn-G)+gama*(37/(Thr+273))*u2*VPD)/(deta+gama+gama*0.34*u2);
System.err.println(deta);
System.err.println(gama);
System.err.println(VPD);
return eto;
}
}
***위의 내용에서 잘못 계산되거나 잘못된 파라미터 값이 있다면 알려주세요
***또는 코드를 사용하시는 분들은 댓글 부탁드립니다.
댓글
댓글 쓰기