#!/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.

# First we set the period of analysis.
$MaxYear=2000;
$MinYear=1900;

# This sets the minimum number of years that need to be 
#  in the rainfall record of a station for it to be used
$MinCount=1*($MaxYear-$MinYear)/2;

# Here we read in the list of data files
open(FILELIST,"ls prcp|");

# This loop processes the data files, parsing the name, 
#  and then reading in the data
$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 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 10 mm. 
#
# This should be almost enough information to start calculating...

# This variable sets the minimum standard deviation for rainfall 
#  anomalies.
$MinDevRain=10; # 10 tenths of a millimeter

# In this loop we calculate the standard deviation of monthly rainfall
#  for each station for each month.  We store the standard deviation 
#  rainfall, and if this is below the minimum value, we set it to the 
#  minimum 
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;
    }
  }
 }
}

# In this set of loops, we calculate the monthly anomaly index 
#  for each station, for each month and year, and we calculate the average 
#  anomaly index for all stations for each month and year. We weight 
#  the average so that stations with longer records have more weight 
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];
      $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;
   }
   printf "%5.2f ",$AvgAnom1[$year][$mon];
  }
}

# Give a new name to the RainAnom array.
@AnomTot=@RainAnom;


# The next step is to factor out the first order spacial anomaly mode, 
#  and to calculate the higher order rainfall anomalies.  
#
# How do we do this?  For each station and each month, we have a series
#  of annual anomalies, and the question is to what extent this series
#  of anomalies matches the regional pattern. So we compare the station 
#  series to the regional series, and calculate that component of the 
#  station's anomaly that is explained by the first order regional 
#  pattern. 
#
# We then subtract the first order anomaly component of each station 
#  from the full rainfall anomaly, and this gets us the first order 
#  anomaly residual.
#
# As a reporting measure, we should report the RMS anomaly 
#  amplitude, for each month and each year, the RMS of the first order
#  anomaly, and the RMS of the residual.
#
# Now given this residual, we need to calculate that anomaly pattern
#  that explains the most of the remaining residual. 
#
# One method for doing this is to start with a random spatial and 
#  temporal pattern, and then to project it onto the Anomaly 
#  residual pattern, renomalize the result, and then project again.
# 
# This is one way of coming up with a 'fastest growing eigenvector'
#  or the pattern explains the most of the anomaly residual pattern.
#
# So now armed with this computational methodology, lets just do it.

# First note that we have the following variables already calculated:
#
# RainAnom[$StnIndex][$year][$mon] = The rainfall anomaly
# RainStdDev[$StnIndex][$mon] = The rainfall standard deviation 
# RainAvgCnt[$StnIndex][$mon] = The number of available data points
# DataAvail[$StnIndex][$year][$mon] = Index of data availability
# AvgAnom1[$year][$mon] = The average first order anomaly index
# AnomTot[$StnIndex][$year][$mon] = Total normalized anomaly
# 
# We next have to calculate the following:
#
# AnomTotRMS[$year][$mon] = Total RMS rainfall anomaly index
# Anom1Amp[$StnIndex][$mon] = The amplitude of the first order anomaly.
# Anom1RMS[$year][$mon] = Anomaly 1 RMS.
# Resid2[$StnIndex][$year][[$mon] = Second order anomaly residual
# Resid2RMS[$year][$mon] = Anomaly  2 RMS.
# AvgAnom2[$year][$mon] = The average second order anomaly index
# Anom2Amp[$StnIndex][$mon] = The amplitude of the second order anomaly.
# Anom2RMS[$year][$mon] = The amplitude of the second order anomaly.
# Resid3[$StnIndex][$year][[$mon] = Third order anomaly residual
# Resid3RMS[$year][$mon] = Anomaly 3 RMS.
#


# Next, we calculate the projection of the first order anomaly onto 
#  the anomalies of the individual rainfall stations and subtract 
#  that out of the total anomaly to calculate the residual
#
# The inputs of this calculation are:
#
#  1) $AvgAnomNP1[$year][$mon]
#  2) $RainCnt[$StnIndex][$mon]
#  3) $AnomN[$StnIndex][$year][$mon]
#
# The outputs of the calculation are:
#
#  1) $AnomNP1[$StnIndex][$year][$mon]

# The following lines are also available in a different 
#  form in the subroutine NextOrderAnom

print "\n\nPrinting First Order Anomaly Amplitudes:\n";
for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
 if($StnCount{$StnIndex} >= $MinCount){
  print "\nStn= $StnIndex";
  for($mon=1;$mon<=12;$mon++){
   $Anom1VecAmp=0;
   $Anom1VecAvg=0;
   for($year=$MinYear;$year<=$MaxYear;$year++){
     $Anom1Vec[$year]=$AvgAnom1[$year][$mon]
		*$DataAvail[$StnIndex][$year][$mon];
     $Anom1VecAvg=$Anom1VecAvg+$Anom1Vec[$year];
   }
   if ($RainCnt[$StnIndex][$mon] >= 1){
    $Anom1VecAvg=$Anom1VecAvg/$RainCnt[$StnIndex][$mon];
   }else{
    $Anom1VecAvg=0.;
   }
   for($year=$MinYear;$year<=$MaxYear;$year++){
     # Make sure anomaly vector has zero mean
     $Anom1Vec[$year]=$Anom1Vec[$year]-$Anom1VecAvg;
     $Anom1VecAmp=$Anom1VecAmp+$Anom1Vec[$year]**2;
   }
   $Anom1VecAmp=sqrt($Anom1VecAmp);
   $Anom1Amp[$StnIndex][$mon]=0;
   if($Anom1VecAmp > 0){
    for($year=$MinYear;$year<=$MaxYear;$year++){
     $Anom1Vec[$year]=$Anom1Vec[$year]/$Anom1VecAmp;
     $Anom1Amp[$StnIndex][$mon]=$Anom1Amp[$StnIndex][$mon]+
	$AnomTot[$StnIndex][$year][$mon]*$Anom1Vec[$year];
    }
   }else{
    $Anom1Amp[$StnIndex][$mon]=0;
   }
   for($year=$MinYear;$year<=$MaxYear;$year++){
    $Anom1[$StnIndex][$year][$mon]=$Anom1Amp[$StnIndex][$mon]*$Anom1Vec[$year];
    $Resid2[$StnIndex][$year][$mon]=$AnomTot[$StnIndex][$year][$mon]
		-$Anom1[$StnIndex][$year][$mon];
   }
   printf "%5.2f ",$Anom1Amp[$StnIndex][$mon];
  }
 }
}
#  The above section is also available in the section NextOrderAnom
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

# First we calculate the AnomTotRMS array and output the result. 
print "\n\nPrinting Total Anomaly Index RMS:\n";
for($year=$MinYear;$year<=$MaxYear;$year++){
  for($mon=1;$mon<=12;$mon++){
   $AnomTotRMS[$year][$mon]=0;
   $Anom1RMS[$year][$mon]=0;
   $Resid2RMS[$year][$mon]=0;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
     if ($DataAvail[$StnIndex][$year][$mon] == 1){
      $AnomTotRMS[$year][$mon]=$AnomTotRMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($AnomTot[$StnIndex][$year][$mon])**2;
      $Anom1RMS[$year][$mon]=$Anom1RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Anom1[$StnIndex][$year][$mon])**2;
      $Resid2RMS[$year][$mon]=$Resid2RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Resid2[$StnIndex][$year][$mon])**2;
      $AvgAnomCount[$year][$mon]=$AvgAnomCount[$year][$mon]+
	+$RainCnt[$StnIndex][$mon];
     }
    }
   }
   if ($AvgAnomCount[$year][$mon] >= 1){
    $AnomTotRMS[$year][$mon]=sqrt($AnomTotRMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Anom1RMS[$year][$mon]=sqrt($Anom1RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Resid2RMS[$year][$mon]=sqrt($Resid2RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
   }else{
    $AnomTotRMS[$year][$mon]=0;
    $Anom1RMS[$year][$mon]=0;
    $Resid2RMS[$year][$mon]=0;
   }
  }
  printf "\n%4d AnT",$year;
  for($mon=1;$mon<=12;$mon++){
   printf "%5.2f ",$AnomTotRMS[$year][$mon];
  }
  printf "\n%4d An1",$year;
  for($mon=1;$mon<=12;$mon++){
   printf "%5.2f ",$Anom1RMS[$year][$mon];
  }
  printf "\n%4d Rs2",$year;
  for($mon=1;$mon<=12;$mon++){
   printf "%5.2f ",$Resid2RMS[$year][$mon];
  }
}

# Next, we wish to decompose the second order residual into the next
#  anomaly mode and a third order residual. So lets try.
#
# The method we do is we create an initial Anom2Amp[$SntIndex][$mon]
#  and an initial AvgAnom2[$year][$mon].  We then multiply  
#  $Resid2[$StnIndex][$year][$mon]*$AvgAnom2[$year][$mon]*
#  *$Anom2Amp[$StnIndex][$mon], and average over the stations to 
#  get the new AvgAnom2, and average over the years to get the new
#  Anom2Amp.  We then renormalize and iterate to convergence. 
#
#  So lets try:

for ($order=2;$order<=4;$order++){
 $itmax=5;
 $bitmax=5;
 if ($order == 2){
  @Resid=@Resid2;
 }elsif($order == 3){
  @Resid=@Resid3;
 }elsif($order == 4){
  @Resid=@Resid4;
 }
 #@Resid=@AnomTot;
 
 # Initialize with random values
 for($mon=1;$mon<=12;$mon++){
  for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    $StnAmp[$StnIndex][$mon]=(rand 2)-1.;
  }
 }
 
 $innoc=0.01;
 for($bit=1;$bit<=$bitmax;$bit++){
  # Innoculate with random values
  for($mon=1;$mon<=12;$mon++){
   $StnAmpTot=0.;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
     $StnAmp[$StnIndex][$mon]=$StnAmp[$StnIndex][$mon]+$innoc*((rand 2)-1.);
     $StnAmpTot=$StnAmpTot+$StnAmp[$StnIndex][$mon];
   }
   if ($StnAmpTot < 0){
    for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
     $StnAmp[$StnIndex][$mon]=-$StnAmp[$StnIndex][$mon];
    }
   }
  }
  for($it=1;$it<=$itmax;$it++){
  # So I first take a random station pattern
  # Multiply times the station residuals to calculate the 
  # Annual anomaly for that pattern. Then given the resulting 
  # Anomaly, one re-calculates the spacial pattern. 
  # We do this in a subrouting IterateAnomCalc, which takes the 
  # $StnAmp[$StnIndex][$mon], $DataAvail[$StnIndex][$year][$mon],
  # $Resid[$StnIndex][$year][$mon], and returns
  # $AnnAnom[$year][$mon], $StnAmpIt[$StnIndex][$mon]
   &IterateAnomCalc;
   @StnAmp=@StnAmpIt;
 
   # Print out the new station amplitudes:
   print "Station Amplitudes Average Iteration = $it";
   $StnAvg=0;
   $StnAvgCnt=0;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
 #    print "\nStn=$StnIndex";
     for($mon=1;$mon<=12;$mon++){
 #     printf "%5.2f ",$StnAmpIt[$StnIndex][$mon];
      $StnAvg=$StnAvg+abs($StnAmpIt[$StnIndex][$mon]);
      $StnAvgCnt++;
     }
    }
   }
   $StnAvg=$StnAvg/$StnAvgCnt;
   printf " %6.4f\n",$StnAvg;
  }
 }
 
 # Now that we have AnnAnom[$year][$mon] we can calculate 
 #  the next order anomaly. 
 @AvgAnomN=@AnnAnom;
 @AnomNM1=@Resid;
 &NextOrderAnom;
 
 if ($order==2){
  @Anom2=@AnomN;
  @Resid3=@ResidN;
  @AnomN=();
  @ResidN=();
 }elsif($order==3){
  @Anom3=@AnomN;
  @Resid4=@ResidN;
  @AnomN=();
  @ResidN=();
 }elsif($order==4){
  @Anom4=@AnomN;
  @Resid5=@ResidN;
  @AnomN=();
  @ResidN=();
 }

} #End Order loop


# First we calculate the AnomTotRMS array and output the result. 
open(CC,">calcD.txt");
print CC "\n\nPrinting Total Anomaly Index RMS:\n";
for($year=$MinYear;$year<=$MaxYear;$year++){
  for($mon=1;$mon<=12;$mon++){
   $AnomTotRMS[$year][$mon]=0;
   $Anom1RMS[$year][$mon]=0;
   $Resid2RMS[$year][$mon]=0;
   $Anom2RMS[$year][$mon]=0;
   $Resid3RMS[$year][$mon]=0;
   $Anom3RMS[$year][$mon]=0;
   $Resid4RMS[$year][$mon]=0;
   $Anom4RMS[$year][$mon]=0;
   $Resid5RMS[$year][$mon]=0;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
     if ($DataAvail[$StnIndex][$year][$mon] == 1){
      $AnomTotRMS[$year][$mon]=$AnomTotRMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($AnomTot[$StnIndex][$year][$mon])**2;
      $Anom1RMS[$year][$mon]=$Anom1RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Anom1[$StnIndex][$year][$mon])**2;
      $Resid2RMS[$year][$mon]=$Resid2RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Resid2[$StnIndex][$year][$mon])**2;
      $Anom2RMS[$year][$mon]=$Anom2RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Anom2[$StnIndex][$year][$mon])**2;
      $Resid3RMS[$year][$mon]=$Resid3RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Resid3[$StnIndex][$year][$mon])**2;
      $Anom3RMS[$year][$mon]=$Anom3RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Anom3[$StnIndex][$year][$mon])**2;
      $Resid4RMS[$year][$mon]=$Resid4RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Resid4[$StnIndex][$year][$mon])**2;
      $Anom4RMS[$year][$mon]=$Anom4RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Anom4[$StnIndex][$year][$mon])**2;
      $Resid5RMS[$year][$mon]=$Resid5RMS[$year][$mon]+
 	+$RainCnt[$StnIndex][$mon]*($Resid5[$StnIndex][$year][$mon])**2;
      $AvgAnomCount[$year][$mon]=$AvgAnomCount[$year][$mon]+
	+$RainCnt[$StnIndex][$mon];
     }
    }
   }
   if ($AvgAnomCount[$year][$mon] >= 1){
    $AnomTotRMS[$year][$mon]=sqrt($AnomTotRMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Anom1RMS[$year][$mon]=sqrt($Anom1RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Resid2RMS[$year][$mon]=sqrt($Resid2RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Anom2RMS[$year][$mon]=sqrt($Anom2RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Resid3RMS[$year][$mon]=sqrt($Resid3RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Anom3RMS[$year][$mon]=sqrt($Anom3RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Resid4RMS[$year][$mon]=sqrt($Resid4RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Anom4RMS[$year][$mon]=sqrt($Anom4RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
    $Resid5RMS[$year][$mon]=sqrt($Resid5RMS[$year][$mon]/
		$AvgAnomCount[$year][$mon]);
   }else{
    $AnomTotRMS[$year][$mon]=0;
    $Anom1RMS[$year][$mon]=0;
    $Resid2RMS[$year][$mon]=0;
    $Anom2RMS[$year][$mon]=0;
    $Resid3RMS[$year][$mon]=0;
    $Anom3RMS[$year][$mon]=0;
    $Resid4RMS[$year][$mon]=0;
    $Anom4RMS[$year][$mon]=0;
    $Resid5RMS[$year][$mon]=0;
   }
  }
  printf CC "\n%4d AnT",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$AnomTotRMS[$year][$mon];
  }
  printf CC "\n%4d An1",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Anom1RMS[$year][$mon];
  }
  printf CC "\n%4d Rs2",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Resid2RMS[$year][$mon];
  }
  printf CC "\n%4d An2",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Anom2RMS[$year][$mon];
  }
  printf CC "\n%4d Rs3",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Resid3RMS[$year][$mon];
  }
  printf CC "\n%4d An3",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Anom3RMS[$year][$mon];
  }
  printf CC "\n%4d Rs4",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Resid4RMS[$year][$mon];
  }
  printf CC "\n%4d An4",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Anom4RMS[$year][$mon];
  }
  printf CC "\n%4d Rs5",$year;
  for($mon=1;$mon<=12;$mon++){
   printf CC "%5.2f ",$Resid5RMS[$year][$mon];
  }
}
close(CC);

exit 0;

#------------Subroutine NextOrderAnom--------------------#

sub NextOrderAnom{

 # Next, we calculate the projection of the first order anomaly onto 
 #  the anomalies of the individual rainfall stations and subtract 
 #  that out of the total anomaly to calculate the residual
 #
 # The inputs of this calculation are:
 #
 #  1) $AvgAnomN[$year][$mon]
 #  2) $RainCnt[$StnIndex][$mon]
 #  3) $AnomNM1[$StnIndex][$year][$mon]
 #  4) $StnAmp[$StnIndex][$mon]
 #
 # The outputs of the calculation are:
 #
 #  1) $AnomN[$StnIndex][$year][$mon]
 #  2) $ResidN[$StnIndex][$year][$mon]
 
 #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
 #            WE NEED TO TURN THIS INTO A SUBROUTINE
 
 print "\n\nPrinting First Order Anomaly Amplitudes:\n";
 for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
  if($StnCount{$StnIndex} >= $MinCount){
   print "\nStn= $StnIndex";
   for($mon=1;$mon<=12;$mon++){
    $AnomNVecAmp=0;
    $AnomNVecAvg=0;
    for($year=$MinYear;$year<=$MaxYear;$year++){
      $AnomNVec[$year]=$AvgAnomN[$year][$mon]
 		*$DataAvail[$StnIndex][$year][$mon];
      $AnomNVecAvg=$AnomNVecAvg+$AnomNVec[$year];
    }
    if ($RainCnt[$StnIndex][$mon] >= 1){
     $AnomNVecAvg=$AnomNVecAvg/$RainCnt[$StnIndex][$mon];
    }else{
     $AnomNVecAvg=0.;
    }
    for($year=$MinYear;$year<=$MaxYear;$year++){
      # Make sure anomaly vector has zero mean
      $AnomNVec[$year]=$AnomNVec[$year]-$AnomNVecAvg;
      $AnomNVecAmp=$AnomNVecAmp+$AnomNVec[$year]**2;
    }
    $AnomNVecAmp=sqrt($AnomNVecAmp);
    $AnomNAmp[$StnIndex][$mon]=0;
    if($AnomNVecAmp > 0){
     for($year=$MinYear;$year<=$MaxYear;$year++){
      $AnomNVec[$year]=$AnomNVec[$year]/$AnomNVecAmp;
      $AnomNAmp[$StnIndex][$mon]=$AnomNAmp[$StnIndex][$mon]+
  	$AnomNM1[$StnIndex][$year][$mon]*$AnomNVec[$year];
     }
    }else{
     $AnomNAmp[$StnIndex][$mon]=0;
    }
    for($year=$MinYear;$year<=$MaxYear;$year++){
     $AnomN[$StnIndex][$year][$mon]=$AnomNAmp[$StnIndex][$mon]*$AnomNVec[$year];
     $ResidN[$StnIndex][$year][$mon]=$AnomNM1[$StnIndex][$year][$mon]
 		-$AnomN[$StnIndex][$year][$mon];
    }
    printf "%5.2f ",$AnomNAmp[$StnIndex][$mon];
   }
  }
 }
 #              WE NEED TO TURN THE ABOVE INTO A SUBROUTINE 
 #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
}

#-------------Subroutine IterateAnomCalc------------------------#
 # We do this in a subrouting IterateAnomCalc, which takes the 
 # $StnAmp[$StnIndex][$mon], $DataAvail[$StnIndex][$year][$mon],
 # $Resid[$StnIndex][$year][$mon], and returns
 # $AnnAnom[$year][$mon], $StnAmpIt[$StnIndex][$mon]

sub IterateAnomCalc{

  #First we normalize the $StnAmp
  for($mon=1;$mon<=12;$mon++){
   $StnAmpNorm=0.;
   $StnCnt=0;
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
     $StnAmpNorm=$StnAmpNorm+$StnAmp[$StnIndex][$mon]**2;
     $StnCnt++;
    }
   }
   if($StnCnt >=1 ){
    $StnAmpNorm=sqrt($StnAmpNorm/$StnCnt);
    for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
     if($StnCount{$StnIndex} >= $MinCount){
      if($StnAmpNorm > 0){
       $StnAmp[$StnIndex][$mon]=$StnAmp[$StnIndex][$mon]/$StnAmpNorm;
      }else{
       $StnAmp[$StnIndex][$mon]=0.;
      }
     }else{
      $StnAmp[$StnIndex][$mon]=0.;
     }
    }
   }else{
    for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
      $StnAmp[$StnIndex][$mon]=0.;
    }
   }
  }
# Comment out the prints?
#  print "Station Amplitudes, Iteration = $it\n";
#  for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
#   if($StnCount{$StnIndex} >= $MinCount){
#    print "\nStn=$StnIndex";
#    for($mon=1;$mon<=12;$mon++){
#    printf "%5.2f ",$StnAmp[$StnIndex][$mon];
#    }
#   }
#  }

  # Second we calculate the AnnAnom from the StnAmp 
  for($mon=1;$mon<=12;$mon++){
   for($year=$MinYear;$year<=$MaxYear;$year++){
    $AnnAnom[$year][$mon]=0.;
    $StnCnt=0;
    for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
     if($StnCount{$StnIndex} >= $MinCount){
      if ($DataAvail[$StnIndex][$year][$mon] == 1){
       $AnnAnom[$year][$mon]=$AnnAnom[$year][$mon]+
		$StnAmp[$StnIndex][$mon]*$Resid[$StnIndex][$year][$mon];
       $StnCnt++;
      }
     }
    }
    if ($StnCnt >= 1){
     $AnnAnom[$year][$mon]=$AnnAnom[$year][$mon]/$StnCnt;
    }else{
     $AnnAnom[$year][$mon]=0.;
    }
   }
  }

  # Next from the AnnAnom we calculate the StnAmpIt
  for($mon=1;$mon<=12;$mon++){
   for($StnIndex=1;$StnIndex<=$StnIndexMax;$StnIndex++){
    if($StnCount{$StnIndex} >= $MinCount){
     $StnAmpIt[$StnIndex][$mon]=0.;
     $AnomCnt=0;
     for($year=$MinYear;$year<=$MaxYear;$year++){
      if ($DataAvail[$StnIndex][$year][$mon] == 1){
       $StnAmpIt[$StnIndex][$mon]=$StnAmpIt[$StnIndex][$mon]+
		$AnnAnom[$year][$mon]*$Resid[$StnIndex][$year][$mon];
       $AnomCnt++;
      }
     }
     if ($AnomCnt >= 1){
      $StnAmpIt[$StnIndex][$mon]=$StnAmpIt[$StnIndex][$mon]/$AnomCnt;
     }else{
      $StnAmpIt[$StnIndex][$mon]=0;
     }
    }else{
     $StnAmpIt[$StnIndex][$mon]=0;
    }
   }
  }
}