#!/usr/bin/perl -w

use strict;
use lib '/pm/grflib/perl';
use MODFLOW;
use DBF;
use GetOpt;

my %args = GetOpt('d:' , d=>'../data/monet.dbf');
my $file = $args{d};
(@ARGV>0) || die "Usage: $0 <year(s)>\n";

#
#  Number of days in month
#
my @ndm = (undef,31,28,31,30,31,30,31,31,30,31,30,31);
my $Ncol = 326;
my $Nrow = 165;

# Station    Long   Lat    Easting  Northing  Row Column
# Akron    -103:09 40:09    480549  14607776   68    41
# McCook   -100:35 40:14   1198494  14614955   67   177
# RedCloud  -98:17 40:05   1840882  14557195   68   299
my @stn = ('Akron','McCook','RedCloud');
my %col = ('Akron'=>47,'McCook'=>177,'RedCloud'=>299);

my $modflow = MODFLOW->new(Ni=>165,Nj=>326,Nk=>1,X0=>266023,Y0=>14092806,DX=>326*5280,DY=>165*5280);
my $dbf = DBF->open($file);
($dbf) || die "Cannot open DBF file $file\n";
$dbf->index('Month','N');

#
#  Read ET area
#
my $dat = $modflow->open('../data0/etarea.dat','TXT');
defined($dat) || die "Cannot open file ../data0/etarea.dat\n";
$modflow->skip($dat,1);
my (undef,undef,undef,undef,@area) = $modflow->read($dat,'(326F6.1)');
$modflow->close($dat);
#  Convert acres to fraction
@area = map {$_/640} @area;

#
#  Output file
#
my $num=0;
printf "%10d%10d\n" , 1 , 40;
foreach my $ARGV (@ARGV)
{
   my ($yr0,$yr1) = split('-',$ARGV);
   defined($yr1) || ($yr1=$yr0);
   foreach my $yr ($yr0 .. $yr1)
   {
      foreach my $mo (1 .. 12)
      {
         my %in = $dbf->lookup('Month',10000*$yr+100*$mo,'N');
         my $nsec = 24*60*60*(($yr%4==0 && $mo==2) ? 29 : $ndm[$mo]);

         #  Header
         printf "%10d%10d%10d%10d\n" , ($num)?-1:1 , 1 , ($num)?-1:1 , 0;
         #  ET Surface (for first stress period)
         ($num) || printf "OPEN/CLOSE ../static/11.etsurf 1.0 (free) -1\n";

         #
         #  East-west distribution of ET rate by column based on Akron, McCook and RedCloud
         #
         my @rate;
         for (my $j=1;$j<$col{McCook};$j++)
         {
            $rate[$j] = ($j-$col{Akron})/($col{McCook}-$col{Akron})*($in{McCook}-$in{Akron}) + $in{Akron};
         }
         for (my $j=$col{McCook};$j<=$Ncol;$j++)
         {
            $rate[$j] = ($j-$col{McCook})/($col{RedCloud}-$col{McCook})*($in{RedCloud}-$in{McCook}) + $in{McCook};
         }
         #
         #  Sanity check
         #
         for (my $j=1;$j<=$Ncol;$j++)
         {
            ($rate[$j]<0) && ($rate[$j] = 0); #  If end is 0, extrapolation could be negative
            ($rate[$j]>60) && warn "Huge rate $mo/$yr $j $rate[$j]\n";
         }
         #
         #  Output ET rate
         #
         my $in2fps = 1/12/$nsec;  #  Convert inches to feet per second
         print "INTERNAL $in2fps (free) -1  $mo/$yr\n";
         my $k=0;
         for (my $i=1;$i<=$Nrow;$i++)
         {
            for (my $j=1;$j<=$Ncol;$j++,$k++)
            {
               printf "%6.2f" , $area[$k]*$rate[$j];
            }
            print "\n";
         }
         #  Extinction depth (for first stress period)
         ($num) || printf "CONSTANT 10.0\n";
         $num++;
      }
   }
}
