#!/usr/bin/perl

# This script is meant to go through a series of WMO 
#  precipitation stations, rank the stations according 
#  to the completeness of record between 1979 to 2000 
#  inclusive, and then calculate the different rainfall 
#  anomalies for this period. 
#
# We are interested in calculating anomalies for different
#  times of year, and we are interested in decomposing the 
#  anomalies into different spacial patterns. 
#
# We hope to calculate an annual scalar anomaly index 
#  for different spacial modes of anomaly.

$MaxYear=2000;
$MinYear=1940;
$MinCount=2*($MaxYear-$MinYear)/3;

open(FILELIST,"ls prcp|");

$StnIndex=0;
while(<FILELIST>){
  chop;
  $FileName=$_;
  $StnName=$FileName;
  $StnName=~ s/\.prcp//;
  $StnIndex++;
  open(DATAFILE,"prcp/".$FileName);
  $count=0;
  while(<DATAFILE>){
   $year= substr $_,12,4;
   if(($year<=$MaxYear)&&($year>=MinYear)){
     $count++;
     for($mon=1;$mon<=12;$mon++){
      $offset=16+($mon-1)*5;
      $raindata=substr $_,$offset,5;
      $Rain[$StnIndex][$year][$mon]=0;
      $DataAvail[$StnIndex][$year][$mon]=0;
      if($raindata>=0){
       $Rain[$StnIndex][$year][$mon]=$raindata;
       $DataAvail[$StnIndex][$year][$mon]=1;
      }
     }
   }
  }
  close(DATAFILE);
  $StnName{$StnIndex}=$StnName;
  $StnCount{$StnIndex}=$count;
  printf "%2d  $StnName \n",$count;
}
$StnIndexMax=$StnIndex;
close(FILELIST);

# Now that we have all of the data loaded, the first thing we 
#  do is calculate all of the averages for all of the relevant 
#  stations. 

for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
 if($StnCount{$StnIndex} >= $MinCount){
  printf "%19s ",$StnName{$StnIndex};
  for($mon=1;$mon<=12;$mon++){
   $RainAvg[$StnIndex][$mon]=0;
   $RainCnt[$StnIndex][$mon]=0;
   for($year=$MinYear;$year<=$MaxYear;$year++){
    $RainAvg[$StnIndex][$mon]=$RainAvg[$StnIndex][$mon]+
			$Rain[$StnIndex][$year][$mon];
    $RainCnt[$StnIndex][$mon]=$RainCnt[$StnIndex][$mon]+
			$DataAvail[$StnIndex][$year][$mon];
   }
   if ($RainCnt[$StnIndex][$mon]>0){
    $RainAvg[$StnIndex][$mon]=$RainAvg[$StnIndex][$mon]/
		$RainCnt[$StnIndex][$mon];
   }
   printf "%4d ",$RainAvg[$StnIndex][$mon];
  }
  print "\n";
 }
}

# Now that we have most of the station monthly normals calculated, 
#  our next task is to calculate the different modes of rainfall
#  anomaly. 
#
# What should our method be? Well, this will depend on the assumed
#  mathematical model for the anomalies.  We will assume a linear
#  model where we classify things into different anomaly spacial 
#  modes, and each mode has a single scalar value for each 
#  year. Each mode also has a fixed spacial distribution that 
#  is independent of time. 
#
# In determining the wieghts for calculating different modes, we 
#  need to also consider how much we expect anomaly amplitudes to 
#  vary with average monthly rainfall and data record length.  
#
# Note if a one station has a record that is twice as complete as 
#  another station, we want to count its anomaly data with 
#  twice the weight. 
#  
# Meanwhile, we expect the relative size of the monthly anomaly 
#  to scale with the standard deviation of the rainfall anomaly   
#  for that month and that station.  
#
# So that arid stations don't have an undue influence on the 
#  anomaly calculation, we will set the minimum standard 
#  deviation to a minimum value called MinDevRain. 
#
# This should be almost enough information to start calculating...

$MinRain=100; # 100 tenths of a millimeter
$MinDevRain=10; # 10 tenths of a millimeter

for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
 if($StnCount{$StnIndex} >= $MinCount){
  for($mon=1;$mon<=12;$mon++){
    $RainAvgCnt[$StnIndex][$mon]=0;
    $RainStdDev[$StnIndex][$mon]=0;
    for($year=$MinYear;$year<=$MaxYear;$year++){
     if ($DataAvail[$StnIndex][$year][$mon] == 1){
      $diff=$Rain[$StnIndex][$year][$mon]-$RainAvg[$StnIndex][$mon];
      $RainAvgCnt[$StnIndex][$mon]++;
      $RainStdDev[$StnIndex][$mon]=$RainStdDev[$StnIndex][$mon]+$diff*$diff;
     }
    }
    if($RainAvgCnt[$StnIndex][$mon]>0){
     $RainStdDev[$StnIndex][$mon]=sqrt($RainStdDev[$StnIndex][$mon]/
 	$RainAvgCnt[$StnIndex][$mon]);
     if($RainStdDev[$StnIndex][$mon]<$MinDevRain){
      $RainStdDev[$StnIndex][$mon]=$MinDevRain;
     }
    }else{
     $RainStdDev[$StnIndex][$mon]=$MinDevRain;
    }
  }
 }
}


for($year=$MinYear;$year<=$MaxYear;$year++){
  printf "\n%4d ",$year;
  for($mon=1;$mon<=12;$mon++){
   $AvgAnom1[$year][$mon]=0;
   $AvgAnomCount[$year][$mon]=0;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
     if ($DataAvail[$StnIndex][$year][$mon] == 1){
      $anom=$Rain[$StnIndex][$year][$mon]-$RainAvg[$StnIndex][$mon];
#      $anomnorm=$RainAvg[$StnIndex][$mon]+$MinRain;
#      if($anom>=0){$sign=1;
#      }else{$sign=-1;}
#      $RainAnom[$StnIndex][$year][$mon]=$sign*sqrt(abs($anom/$anomnorm));
      $RainAnom[$StnIndex][$year][$mon]=$anom/$RainStdDev[$StnIndex][$mon];
      $AvgAnom1[$year][$mon]=$AvgAnom1[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*$RainAnom[$StnIndex][$year][$mon];
      $AvgAnomCount[$year][$mon]=$AvgAnomCount[$year][$mon]+
	+$RainCnt[$StnIndex][$mon];
     }
    }
   }
   if ($AvgAnomCount[$year][$mon] >= 1){
    $AvgAnom1[$year][$mon]=$AvgAnom1[$year][$mon]/$AvgAnomCount[$year][$mon];
   }else{
    $AvgAnom1[$year][$mon]=0;
   }
   $dyear=$year+($mon-0.5)/12;
   printf "%5.2f ",$AvgAnom1[$year][$mon];
#   printf "%6.2f, %5.2f \n",$dyear,$AvgAnom1[$year][$mon];
  }
}


exit 0;