#!/usr/bin/perl -w

#  Create Colorado groundwater input data files
#  based on well by well acreage assignment.
#
#  Willem A Schreuder
#  March 19, 2018

use strict;
use lib '/pm/grflib/perl';
use GetOpt;
use Transform;
my $xfm = Transform->new('PROJECT UTM83-N13 INV|PROJECT NAD83 INV|PROJECT USGSRR');

my $yr = shift @ARGV;
defined($yr) && ($yr=~/^20\d\d$/) && ($yr>=2017) || die "Usage: $0 <year>\n";

#  Monthly pumping/recharge distribution
my %frac = ('01'=>0,'02'=>0,'03'=>0,'04'=>0,'05'=>0,'06'=>0.147,'07'=>0.320,'08'=>0.328,'09'=>0.205,'10'=>0,'11'=>0,'12'=>0);
#  CCP wells
my %WDID = (6506306=>0,6506993=>0,6506330=>0,6506995=>0,6506800=>0,6506302=>0,6506299=>0,6506797=>0);
#  CCP volumes
my ($CCP,%ccp,%CCP);
open(CCP , "grep ^$yr ../data/co/CCP.dat |") || die "Cannot grep ../data/co/CCP.dat\n";
while (my $line = <CCP>)
{
   my ($mo,$vol) = split(' ',$line);
   $CCP += $CCP{substr($mo,4)} = $vol;
}

#
#  Model domain
#
my ($Ni,$Nj) = (165,326);
my ($X0,$Y0) = (266023,14092806);
my ($DX,$DY) = (5280,5280);
my ($X1,$Y1) = ($X0+$Nj*$DX,$Y0+$Ni*$DY);
sub off {my ($i,$j) = @_; return ($i-1)*$Nj+($j-1);}
my @ibound;
open(DAT , "<../static/02.ibound") || die "Cannot open file ../static/02.ibound\n";
foreach my $i (1 .. $Ni)
{
   my $line = <DAT>; chomp $line;
   foreach my $j (1 .. $Nj)
   {
      $ibound[off($i,$j)] = substr($line,2*($j-1),2);
   }
}
close(DAT);

#
#  Get well location
#
sub welloc
{
   my ($Xutm,$Yutm,$line) = @_;
   $Xutm && $Yutm || (warn "Invalid coordinates $line\n") && return;
   #  Location
   my ($x,$y) = $xfm->forward($Xutm,$Yutm);
   #  Discard wells outside model domain
   ($x<$X0 || $x>$X1 || $y<$Y0 || $y>$Y1) && return;
   #  Discard wells in dead cells
   my ($i,$j) = ($Ni-int(($y-$Y0)/$DY),int(($x-$X0)/$DX)+1);
   my $off = off($i,$j);
   ($ibound[$off] == 0) && return;
   return ($off,$x,$y);
}

#  Initialize arrays
my (%pmp,%rcg,$mi,$ccp);
my @agw = (0) x ($Ni*$Nj);
my @mi  = (0) x ($Ni*$Nj);
foreach my $mo (sort keys %frac)
{
   @{$pmp{$mo}} = (0) x ($Ni*$Nj);
   @{$rcg{$mo}} = (0) x ($Ni*$Nj);
}

#  Snarf meter readings
my (%N,%Pump,%Area);
open(DAT , "<../data/co/$yr.meter") || die "Cannot open file ../data/co/$yr.meter\n";
<DAT>;
while (my $line = <DAT>)
{
   chomp $line;
   my ($use,$cty,$wdid,$pump,undef,$area,$Xutm,$Yutm) = split(',',$line);
   ($pump eq '') && next;
   ($pump>0) || next;
   ($area) || ($area = $pump/1.2);  #  Use 1.2af/acre if no acres are shown
   my ($off,$x,$y) = welloc($Xutm,$Yutm,$line);
   defined($off) || next;
   #  CCP pumping
   if (exists($WDID{$wdid}))
   {
      $WDID{$wdid}++;
      $ccp += $pump;
      $ccp{$off} += $pump;
   }
   #  Agricultural pumping
   elsif ($use eq '1')
   {
      my $rech = 0.17*$pump;
      $agw[$off] += $area;
      while (my ($mo,$f) = each %frac)
      {
         $pmp{$mo}[$off] += $f*$pump;
         $rcg{$mo}[$off] += $f*$rech;
      }
      $N{$cty}++;
      $Pump{$cty} += $pump;
      $Area{$cty} += $area;
   }
   else
   #  M&I pumping
   {
      $mi += $pump;
      $mi[$off] += 0.5*$pump;
   }
}
close(DAT);

#  Add CCP pumping to cells by distributing monthly volumes
#  proportional to annual pumping by well
while (my ($off,$pmp) = each %ccp)
{
   my $f = $pmp/$ccp;
   while (my ($mo,$vol) = each %CCP)
   {
      $pmp{$mo}[$off] += $f*$vol;
   }
}

#  Print totals
printf "M&I    %6d\n" , $mi;
printf "CCPmet %6d\n" , $ccp;
printf "CCPadj %6d\n" , $CCP;
my ($N,$Pump,$Area);
foreach my $cty (sort keys %N)
{
   $N += $N{$cty};
   $Pump += $Pump{$cty};
   $Area += $Area{$cty};
   printf "%-15s %4d %8d %8d\n" , $cty , $N{$cty} , $Pump{$cty} , $Area{$cty};
}
printf "%-15s %4d %8d %8d\n" , 'Ag Total' , $N , $Pump , $Area;

#
#  Save to file
#
sub save
{
   my ($file,@x) = @_;
   #  Skip empty files
   my $sum=0;
   foreach my $x (@x)
   {
      $sum += abs($x);
   }
   $sum || return;
   #  Dump values to file
   open(DAT , ">$file") || die "Cannot open file $file\n";
   print DAT join("\n",@x)."\n";
   close(DAT);
}

#  Annual acres
save("co/$yr.agw",@agw);
#  M&I net pumping
save("co/$yr.mi",@mi);
#  Monthly pumping and recharge
foreach my $mo (sort keys %frac)
{
   save("co/$yr.$mo.pmp",@{$pmp{$mo}});
   save("co/$yr.$mo.rcg",@{$rcg{$mo}});
}
