<?php
// Global Constants

$PI = 3.141593;

function PENMAN_Kimb82($tmax, $tmin, $tprv, $XJDAY, $VAPOR, $S_RAD, $WINDT, $LAT, $ELEV, $ANEMH1)
{
    
    $PI = 3.141593;    
    /*
    This subroutine calculates evapotranspiration rate using 
    the 1982 Penman Combination method as presented in 
    ASCE No. 70, Evapotranspiration and Irrigation Water Requirements, 1990.

    Input:
	$tmax - 24 hr max temp (degF)
	$tmin - 24 hr min temp (degF)
	$tprv - Mean temp for prev. 3 days (degF)
	$XJDAY - Day of year
	$VAPOR - Vapor pressure (mb)
	$S_RAD - Solar radiation (langleys)
	$WINDT - Wind run (miles)
	$LAT   - Latitude (deg)
	$ELEV  - Elevation (ft)
	$ANEMH1- Anemometer height (ft)

    */

    // print "$tmax, $tmin, $tprv, $XJDAY, $VAPOR, $S_RAD, $WINDT, $LAT,$ELEV,$ANEMH1\n";
    
    $TC_MX = f2c($tmax);
    $TC_MN = f2c($tmin);
    $TC_AV = ($TC_MX + $TC_MN) / 2;

    //
    //  Put in reasonable values if these are missing.
    //
    if($LAT    == 0)	$LAT    = 40;	// Latitude
    if($ELEV   == 0)	$ELEV   = 5000;	// Elevation
    if($ANEMH1 == 0)	$ANEMH1 = 6.56;

    $REFCROP = 'ALFA';

    // Calculate saturation vapor pressures in millibars for daily
    // temperature in degrees C, then convert to kPa.

    $SVAP_MX = esat_poly($TC_MX);
    $SVAP_MX /= 10;
  
    $SVAP_MN = esat_poly($TC_MN);
    $SVAP_MN = $SVAP_MN / 10;

    $SVAP_AV = esat_poly($TC_AV);
    $SVAP_AV /= 10;

    $SVAP_DW = $VAPOR;
    $SVAP_DW /= 10;

    // Calculate by method 3 the average daily saturation vapor pressure 
    // SVAP_AV in millibars.
  
    $SVAP_AVE = .5 * ($SVAP_MN + $SVAP_MX);

    // Calculate the slope of the vapor pressure-temperature curve DELTA 
    // in millibars/C.
  
    $DELTA = delta_poly($TC_AV);
    $DELTA /= 10;

    // Calculate the latent heat of vaporization XLAT in MJ/kg.
  
    $XLAT = 2.501 - 2.361E-03 * $TC_AV;
    
    // Calculate the average station atmospheric pressure APRESS in kPa.
    
    $APRESS = 101.325 * pow(((288.16 - 0.0065 * ($ELEV / 3.281)) / 288.16),(9.80665 / 0.0065 / 287.04));

    // Calculate  psychrometer constant GAMMA in kPa/C.

    $GAMMA = 1.013E-03 * $APRESS / 0.622 / $XLAT;

    // Calculate the coefficients C1, C2 and C3.

    $C1 = $DELTA / ($DELTA + $GAMMA);
    $C2 = $GAMMA / ($DELTA + $GAMMA);

    // C3 = 6.43 is common estimate derived from equation 6.70

    $C3 = 0.35 * $XLAT / 0.13332;

    // Calculate the net outgoing long wave radiation RBO on a clear day.
    
    $A1 = .26 + .1 * exp(- pow( 0.0154 * ($XJDAY - 176),2));
    $B1 = -0.139;
    $RBO = ($A1 + $B1 * pow($SVAP_DW,0.5) ) * 4.903E-09 * (pow($TC_MX + 273.16,4) + pow($TC_MN + 273.16,4)) / 2;

    // Calculate clear day solar radiation RSO in cal/cm2.

    // DECL = Declination, radians
    $DECL = 0.4093 * sin( 2 * $PI * ( 284 + $XJDAY ) / 365 );
    
    // DR = Distance frm earth to sun
    $DR = 1 + 0.033 * cos( 2 * $PI * $XJDAY / 365);

    // XLATITUDR = Latitude, radians
    $XLATITUDR = $LAT * 2 * $PI / 360;

    // WS = Sunset hour angle, radians
    $WS = acos(-tan($XLATITUDR) * tan($DECL));
    
    // GSC = Solar constant, MJ/m2/min
    $GSC = 0.08202;

    // RA = Extraterrestrial radiation, MJ/m2/day
    $RA = (24.*60./$PI)*$GSC*$DR*(($WS*sin($XLATITUDR)*sin($DECL))+ (cos($XLATITUDR)*cos($DECL)*sin($WS)));
    $RSO = 0.75 * $RA;

    // Calculate net outgoing longwave radiation RB.
  
    $SOL = $S_RAD * 4.1855E-02 / $RSO;

    if($SOL > 1) $SOL = 1; 

    if($SOL > 0.7)
    {
	$AA = 1.126;
	$BB = -0.07;
    }
    if ($SOL <= 0.7)
    {
	$AA = 1.017;
	$BB = -0.06;
    }

    $RB = ($AA * $SOL + $BB) * $RBO;

    // Set coefficient for ALBEDO or reflected short wave radiation.
    
    // ALBEDO  = .29 + .06 * DSIND( 30. * ( XMO + 0.0333 * XDAY + 2.25 ))
    // $ALBEDO = .29 + .06 * sin(30 * ($XMO + 0.0333 * $XDAY + 2.25)
    $ALBEDO = .29 + .06 * sin(($XJDAY + 97.5) * $PI / 180.0);

    // Calculate the net daily radiation RN in cal/cm2.
    $RN = (1 - $ALBEDO ) * $S_RAD * 4.1855E-02 - $RB;

    // Find soil heat flux G towards the surface in cal/cm2.
  
    $PRE_AV = &f2c($tprv);
    $G = 5.0 * ($TC_AV - $PRE_AV ) * 1.8 * 4.1855E-02;

    // Calculate daily wind travel at standard 2 meter anemometer height 
    // WIND2 in kilometers.

    $ANEMH = $ANEMH1;
    $ANEMH *= 30.48;
    $HC = 50;
    $Z0 = 0.122626 * $HC;
    $DH = (2 / 3) * $HC;
  
    $WIND2 = $WINDT * (5280.*.3048/24./3600.)* log((200 - $DH ) / $Z0) / log(($ANEMH - $DH) / $Z0);

    // WR = Ratio of daytime to nighttime wind
    $WR = 2;
    $WINDD2 = 2 * $WIND2 * $WR / (1 + $WR);
  
    $RR  = .622 * $VAPOR / $APRESS;
    $TKV = ($TC_AV + 273.16) * ((1 + $RR / .622)/(1 + $RR));
    $DENSITY = $APRESS / .28704 / $TKV;

    // Set the wind function coefficients for alfalfa and method 3 of 
    // calculating the mean daily vapor pressure deficit.
  
    $AW = 0.4 + 1.4 * exp(-( pow(($XJDAY - 173) / 58,2)));
    $BW = 0.605 + 0.345 * exp(-( pow(($XJDAY - 243) / 80,2)));

    // Solve Penman Combination Equation for reference 
    // evapotranspiration ETR_PM for alfalfa in inches of H2O/day.

    $ETR_PM = $C1 * ($RN - $G) + $C2 * $C3 * ($AW + $BW * $WIND2) * ($SVAP_AVE - $SVAP_DW);
    $ETR_PM = $ETR_PM / 25.4 / $XLAT;

    if($REFCROP == 'TURF')	$ETR_PM *= .8;
    if($ETR_PM < 0)		$ETR_PM = 0;

    return $ETR_PM;
}


function PAN()
{
    // This subroutine calculates Class A Pan evaporation using 
    // the FAO 24 procedure as presented in Appendix E, ASCE No. 70,
    // Evapotranspiration and Irrigation Water Requirements, 1990.

    // COMMON/BLANK1/TF_MX,TF_MN,TF_AV,TF_DW,VAPOR,RHUM,WINDT,WIND2,
    //+  S_RAD,ET_PAN,ETR_PM,T_RH_HT,ANEMH1,ANEMH2_DATE,ANEMH2,DN_RATIO,
    //+  XLATITUDE,XLONGITUDE,ELEV,ALFA_HT,TURF_HT,VEG_HT,PAN_FETCH,
    //+  P_MX(3),P_MN(3),P_AV(3)

    $XRHUM = $RHUM;
    if($XRHUM < 30)	$XRHUM = 30;
    if($XRHUM > 84)	$XRHUM = 84;
    
    // Calculate daily wind travel at standard 2 meter anemometer height 
    // WIND2 in kilometers.
    
    $ANEMH = $ANEMH1;
    $ANEMH *= 30.48;
    $HC = 50;
    $Z0 = 0.122626 * $HC;
    $DH = (2 / 3) * $HC;

    // $U2 is adjusted windrun  
    $U2 = $WINDT * (5280.*.3048/24./3600.)* log((200 - $DH ) / $Z0) / log(($ANEMH - $DH) / $Z0);
    if($U2 < 84)  $U2 = 84; 
    if($U2 > 700) $U2 = 700;
    $FETCH = 1000;
    $XKP = 0.108 - 0.000331*$U2 + 0.0422*log($FETCH) + 0.1434*log($XRHUM) - 0.000631*( pow(log($FETCH),2))*log($XRHUM);
    $ET_PAN = $ETR_PM / $XKP / 1.25;
}

function esat_poly($TC)
{
    // Polynomial coefficients for saturation vapor pressure in
    // millibars vs. temperature in degrees C.
    $C = array(6.10780E+00,4.43652E-01,1.42895E-02,2.65065E-04, 3.03124E-06,2.03408E-08,6.13682E-11);
    $esat = $C[0] + $TC * ($C[1] + $TC * ($C[2] + $TC * ($C[3] + $TC * ($C[4] + $TC * ($C[5] + $TC * $C[6])))));
    
    if($esat < 0.0) $esat = 0.0; 
    return $esat;
}

function delta_poly($TC)
{
    $C = array(6.10780E+00,4.43652E-01,1.42895E-02,2.65065E-04, 3.03124E-06,2.03408E-08,6.13682E-11);

    $delta = $C[1] + $TC * (2 * $C[2] + $TC * (3 * $C[3] 
		   + $TC * (4 * $C[4] + $TC * (5 * $C[5] + $TC * 6 * $C[6]))));
    if($delta < 0.0) $delta = 0;
    return $delta;
}

function e_sat_bosen($T)
{
    // Calculates Saturation Vapor Pressure from Temperature
    // using the Bosen Equation.
    // Pass:    Temperature (C)
    // Return:  Saturation Vapor Pressure (mb)
    
    $esat = 33.8639 * (pow(0.00738*$T+0.8072,8)-(0.000019*(1.8*$T+48))+0.001316);
    return $esat;
}

function et ($t, $tprv, $doy, $vp, $sr, $wr)
{
    // $t - mean temperature, deg F
    // $tprv - previous 3 days mean temperature, deg F
    // $doy - day of year
    // $vp - vapor pressure, mb
    // $sr - solar radiation, Langleys
    // $wr - wind run, miles

  $C2 = 0.959 - 0.0125 * $t + 0.00004534 * $t * $t;
  $C1 = 1 - $C2;
  $Es = -0.6959 + (0.2946 * $t) - (0.005195 * pow($t,2)) + (0.000089 * pow($t,3));
  $Rbo= (0.37 - 0.044 * sqrt($vp)) * 0.0000001171 * pow(f2c($t)+273.16,4);
  $Rso= 760 * exp(- pow(($doy - 166)/157, 2));
  $Rn = 0.77 * $sr - (0.9 * $sr / $Rso + 0.1) * $Rbo;
  $G  = 5 * ($t - $tprv);
  $ETr= 0.000673 * ($C1 * ($Rn - $G) + 15.36 * $C2 * (1.1 + 0.017 * $wr) * ( $Es - $vp ));

  if($ETr < 0) $ETr = 0;

  return $ETr;
}

function penman ($t, $tprv, $doy, $vp, $sr, $wr)
{
    // $t - mean temperature, deg F
    // $tprv - previous 3 days mean temperature, deg F
    // $doy - day of year
    // $vp - vapor pressure, mb
    // $sr - solar radiation, Langleys
    // $wr - wind run, miles

    $PI = 3.1415926;

    // printf("penman 2 - %3d %6.1f %6.1f %6.3f %6.1f %6.0f", $doy, $t, $tprv, $vp, $sr, $wr);

    $P = 850;			// Mean station pressure
    $z = 2;			// WS Height
    $LAT = 40.433;		// Latitude
    $E = 5120 * .304048;	// Elevation

    //  Convert to metric
    $T_avg = f2c($t);
    $tprv = f2c($tprv);
    $Uz = $wr * 1.609344;
    $Rs = $sr;

    // Latent Heat of Vaporization
    $L = 597.3 - (0.564 * $T_avg);


    // Psychrometric Weighting Term
    $gamma = (0.386 * $P)/$L;
    $delta = 2.01 * ( pow(0.00738 * $T_avg + 0.8072,7)) - 0.00116;
    $C1    = $delta/($delta + $gamma);

    // Saturation Vapor Pressures
    $VP_sat_Tavg = e_sat_bosen($T_avg);
    $VP_sat_TDP  = $vp;

    // Vapor Pressure Deficit
    $VP_deficit = $VP_sat_Tavg - $VP_sat_TDP;

    // Net Solar Radiation
    $A_prime = 31.54 - 0.2734 * $LAT + 0.0007813 * $E;
    $B_prime = -0.2986 + 0.2678 * $LAT + 0.0004102 * $E;
    $C_prime = 2.92;
    $Rso = 23.89 * ($A_prime + $B_prime * cos(2*$PI*$doy/365 - $C_prime));
    $Rbo = 0.325 - 0.044 * sqrt($VP_sat_TDP) * 11.71E-8 * pow($T_avg + 273.15,4);
    $Rb  = $Rbo * (1.22 * $Rs/$Rso - 0.18);
    $Rn  = 0.77 * $Rs - $Rb;
	  
    // Heat Flux
    $G = 9 * ($tprv - $T_avg);

    // Wind Function 
    $U2 = $Uz * pow(2/$z,0.2);
    $Wf = 0.75 + 0.0115 * $U2;

    // Penman Reference Evapotranspiration 
    $ETr = (10/$L) * (($C1*($Rn + $G)) + (15.36 * ( 1 - $C1) * $Wf * $VP_deficit));
    $ETr /= 25.4;
    if($ETr < 0) $ETr = 0; 

    //printf(" %6.3f\n", $ETr / 25.4);

    return $ETr;
}
/*
function get_alphalpha()
{
	$arr = array();
	$arr["crop"] = "Alfalfa";
	$arr["coef_2"] = {0.33748,-0.0024345,0.00026031,-0.0000016861 };
	$arr["coef_3"] = { 1.0,0.0,0.0,0.0 };
	//$arr[0.35,1.00,1.0,1.0
	//4, 24
	//0.0,6.0,6.0
}
*/
?>