#!/usr/bin/perl -w

#  Program to calculate monthly maximum ET for phreatophytes
#  
#  Willem A. Schreuder
#  August 17, 2004

use strict;
use POSIX;

my $debug = 0;
while (@ARGV>0 && $ARGV[0] =~ /^-/)
{
   my $arg = shift @ARGV;
   if ($arg eq '-d')
   {
      $debug++;
   }
   else
   {
      die "Unknown parameter $arg\n";
   }
}
(@ARGV>0) || die "Usage: $0 [-d] <year>\n";

#
#  Constants
#
my $pi = 3.141592653589793238462643;
#  Lower and Upper temperature base for GDD
my $T0 = 40;
my $T1 = 77;
#  Crop coefficients by GDD
my %K0 = (
  Akron    => [{GDD=>473,CC=>0.36},{GDD=> 833,CC=>1.08},{GDD=>3575,CC=>1.08},{GDD=>4200,CC=>0.8},{GDD=>4543,CC=>0.36}] ,
  McCook   => [{GDD=>456,CC=>0.36},{GDD=>1045,CC=>1.08},{GDD=>4075,CC=>1.08},{GDD=>4806,CC=>0.8},{GDD=>5288,CC=>0.36}] ,
  RedCloud => [{GDD=>391,CC=>0.36},{GDD=>1012,CC=>1.08},{GDD=>4087,CC=>1.08},{GDD=>4868,CC=>0.8},{GDD=>5344,CC=>0.36}] ,
);
#  ETr/ETo ratio by month
my %ratio = (
  Akron    => [undef,2.34,1.97,1.80,1.56,1.31,1.28,1.33,1.39,1.63,2.01,2.34,2.56] , 
  McCook   => [undef,2.20,1.88,1.70,1.53,1.30,1.31,1.26,1.23,1.47,1.73,2.04,2.31] ,
  RedCloud => [undef,2.19,1.89,1.71,1.49,1.22,1.23,1.17,1.16,1.34,1.73,2.00,2.17] );
# Station Data
my %stn = (
  Akron    =>  {ID=>'050109' , LAT =>  40.15},
  McCook   =>  {ID=>'255310' , LAT =>  40.20},
  RedCloud =>  {ID=>'257070' , LAT =>  40.10});
my @stn = sort keys %stn;

#
#  Calculate monthly ET
#
sub ETmon
{
   my ($stn,$lat,$yr,@data) = @_;
   my (@PET,@PPT);
   #  Convert latitude to radians
   $lat *= $pi/180;
   #  Restart Growing degree days
   my $gdd = 0;
   #  Loop over days of the year
   ($debug>1) && print "\n$stn\n\nDate      Tmin Tmax   Ppt      GDD    ETo  Ratio    ETr    K0    PET\n";
   for (my $J=0;$J<@data;$J++)
   {
      my $mo = $data[$J]{mo};

      #
      #  Calculate Ra (Extraterrestrial Radiation) for 24-hour periods (mm/day)
      #  Source: REF-ET for Windows ver. 2.0 Appendix 2 (pp. 62-63)
      #
      #  Convert temperatures to centigrade
      my $Tmin = ($data[$J]{Tmin}-32)*5/9;
      my $Tmax = ($data[$J]{Tmax}-32)*5/9;
      #  Inverse relative Earth-Sun distance
      my $Dr = 1 + 0.033*cos(2*$pi/365*$J);
      #  Solar Declination
      my $d = 0.4093*sin(2*$pi/365*$J-1.39);
      #  Sunset hour angle
      my $Ws = acos(-tan($lat)*tan($d));
      # Solar constant
      my $Gsc = 0.0820;
      #  Extraterrestrial Radiation
      my $Ra = 0.408*24*60/$pi*$Gsc*$Dr*($Ws*sin($lat)*sin($d)+cos($lat)*cos($d)*sin($Ws));

      #
      #  Calculate ETo by 1985 Hargreaves Method (mm/day)
      #  Source: REF-ET for Windows ver. 2.0 Appendix 2 (p. 58)
      #
      my $ETo = 0.0023*sqrt($Tmax-$Tmin)*(($Tmin+$Tmax)/2+17.8)*$Ra;

      #  Calculate ETr by station using calibrated ratio (inches/day)
      my $ETr = $ratio{$stn}[$mo]*$ETo / 25.4;

      #
      #  Calculate Potential ET from growing degree days
      #
      #  Accumulate growing degree days (F)
      my $t0 = ($data[$J]{Tmin}<$T0) ? $T0 : $data[$J]{Tmin};
      my $t1 = ($data[$J]{Tmax}>$T1) ? $T1 : $data[$J]{Tmax};
      ($t1>$T0) && ($gdd += ($t0+$t1)/2-$T0);
      #  Select crop coefficient based on GDD by interpolation
      my $K0;
      if ($gdd<=$K0{$stn}[0]{GDD})
      {
         $K0 = $K0{$stn}[0]{CC};
      }
      elsif ($gdd>=$K0{$stn}[-1]{GDD})
      {
         $K0 = $K0{$stn}[-1]{CC};
      }
      else
      {
         for (my $k=0;$k<@{$K0{$stn}} && !defined($K0);$k++)
         {
            ($K0{$stn}[$k]{GDD}<=$gdd && $gdd<$K0{$stn}[$k+1]{GDD}) &&
              ($K0 = ($gdd-$K0{$stn}[$k]{GDD})/($K0{$stn}[$k+1]{GDD}-$K0{$stn}[$k]{GDD})*($K0{$stn}[$k+1]{CC}-$K0{$stn}[$k]{CC})+$K0{$stn}[$k]{CC});
         }
         defined($K0) || die "WTF: K0\n";
      }
      #  Potential ET
      my $PET = $K0*$ETr;
      ($debug>1) && printf "$mo/$data[$J]{dy}/$yr %3d %3d %6.2f %8d %6.3f %6.3f %6.3f %6.3f %6.3f\n" ,
        $data[$J]{Tmin} , $data[$J]{Tmax}, $data[$J]{ppt} , $gdd , $ETo/25.4 , $ratio{$stn}[$mo] , $ETr , $K0 , $PET;

      # Accumulate monthly PET and precip
      $PET[$mo] += $PET;
      $PPT[$mo] += $data[$J]{ppt};
   }
   #
   #  Calculate monthly Crop Irrigation Requirement
   #  Source:  Irrigation Water Requirements SCS TR-21, Sep 1970 (p. 27)
   #
   my @CIR;
   ($debug>1) && print "\n# Month   PET    PPTa   PPTe   CIR\n";
   foreach my $mo (1 .. 12)
   {
      #  Application depth
      my $D = 1;
      my $f = (0.531747+0.295164*$D-0.057697*$D**2+0.003804*$D**3);
      #  Effective Precipitation (Re)
      my $PPTe = (0.70917*$PPT[$mo]**0.82416 - 0.1156)*10**(0.02426*$PET[$mo])*$f;
      #  Sanity checks
      ($PPTe<0) && ($PPTe = 0);
      ($PPTe>$PPT[$mo]) && ($PPTe = $PPT[$mo]);
      ($PPTe>$PET[$mo]) && ($PPTe = $PET[$mo]);
      #  CIR for the month (ETmax)
      $CIR[$mo] = $PET[$mo] - $PPTe;
      ($debug>1) && printf "%2d/%4d %6.3f %6.3f %6.3f %6.3f\n" , $mo , $yr ,$PET[$mo] , $PPT[$mo] , $PPTe , $CIR[$mo];
   }
   return @CIR;
}

#
#  Requested years
#
foreach my $yr (@ARGV)
{
   my %ET;
   #  Loop over ET stations
   foreach my $stn (@stn)
   {
      #  Read temperature and precip for this year
      open(TXT , "<$stn{$stn}{ID}.txt") || die "Cannot open file $stn{$stn}{ID}.txt\n";
      <TXT>;
      my @data;
      while (my $line = <TXT>)
      {
         chomp $line;
         my ($date,$Tmin,$Tmax,$ppt) = split(' ',$line);
         my ($m,$d,$y) = split('/',$date);
         ($y==$yr) || next;
         push @data , {Tmin=>$Tmin,Tmax=>$Tmax,ppt=>$ppt,mo=>$m,dy=>$d};
      }
      (@data == ($yr%4?365:366)) || die "Wrong count in met data $stn $yr\n";
      #  Calculate ETmax for this station
      @{$ET{$stn}} = ETmon($stn,$stn{$stn}{LAT},$yr,@data);
   }
   #  Print result
   ($debug>0) && printf "Month     Akron McCook RedCloud\n";
   foreach my $mo (1 .. 12)
   {
      printf "%4d%.2d00 %6.2f %6.2f %6.2f\n" , $yr , $mo , $ET{Akron}[$mo] , $ET{McCook}[$mo] , $ET{RedCloud}[$mo];
   }
}
