#!/usr/bin/perl -w

#  Create Colorado groundwater input data files
#  based on well by well acreage assignment for 2016
#  This program is unique to 2016 as a result of the 
#  phasing in of meter records.
#
#  Willem A Schreuder
#  March 30, 2017

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

#
#  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);
}

#  Read county flags
my (@flg) = (0);
open(FLG , "<../data0/cty.flg") || die "Cannot open file ../data0/cty.flg\n";
my $line = <FLG>;
chomp $line;
my @cty = split(' ',$line);
my $wid = (@cty>9) ? 2 : 1;
for(my $i=0;$i<$Ni;$i++)
{
   my $line = <FLG>;
   chomp $line;
   for (my $j=0;$j<$Nj;$j++)
   {
      my $l = substr($line,$wid*$j,$wid);
      $flg[$i*$Nj+$j] = $l;
      ($l<0 || $l>@cty) && die "Invalid county id $l\n";
   }
}
close(FLG);
#  Make offset 1
unshift @cty , 'OUT';

#  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 volumes for 2016
my %CCP = (01=>1950,02=>1840,03=>1240,10=>1240,11=>1900,12=>1970);
#  CCP wells
my %WDID = (6506306=>0,6506993=>0,6506330=>0,6506995=>0,6506800=>0,6506302=>0,6506299=>0,6506797=>0);
my ($CCP,%ccp);
foreach my $vol (values %CCP)
{
   $CCP += $vol;
}

#  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);
}
#  The 8 Colorado counties are numbered 1-8
my @N    = (undef,0,0,0,0,0,0,0,0);
my @Pump = (undef,0,0,0,0,0,0,0,0);
my @Area = (undef,0,0,0,0,0,0,0,0);

#  Well list file
open(STE , ">co/wells.ste") || die "Cannot open file co/wells.ste\n";
print STE "\@1\n";

#  Snarf meter readings
open(DAT , "<../data/co/2016.meter") || die "Cannot open file ../data/co/2016.meter\n";
<DAT>;
while (my $line = <DAT>)
{
   chomp $line;
   my ($wdid,undef,$use,$pump,$Xutm,$Yutm,undef,$a15,undef,$a16,$type) = split(',',$line);
   ($pump>0 && $Xutm) || next;
   my $area = ($a16==0 && $a15) ? $a15 : $a16;
   ($area==0) && ($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;
      printf STE "%7d %9d 3 %s\n" , $x , $y , $wdid;
   }
   #  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;
      }
      my $cty = ($flg[$off]>8) ? 0 : $flg[$off];
      $N[$cty]++;
      $Pump[$cty] += $pump;
      $Area[$cty] += $area;
      my $afa = ($area==0) ? 'NaN' : sprintf '%.3f' , $pump/$area;
      printf STE "%7d %9d 1 %-6s %6s\n" , $x , $y , $wdid , $afa;
   }
   else
   #  M&I pumping
   {
      $mi += $pump;
      $mi[$off] += 0.5*$pump;
      printf STE "%7d %9d 2 %s\n" , $x , $y , $wdid;
   }
}
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 $k (0 .. 8)
{
   $N += $N[$k];
   $Pump += $Pump[$k];
   $Area += $Area[$k];
   printf "%-12s %4d %8d %8d\n" , $cty[$k] , $N[$k] , $Pump[$k] , $Area[$k];
}
printf "%-12s %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/2016.agw",@agw);
#  M&I net pumping
save("co/2016.mi",@mi);
#  Monthly pumping and recharge
foreach my $mo (sort keys %frac)
{
   save("co/2016.$mo.pmp",@{$pmp{$mo}});
   save("co/2016.$mo.rcg",@{$rcg{$mo}});
}
