#!/usr/bin/perl -w

####################################################################
#British Crown Copyright (c) 2009, the Met Office All rights reserved.
#
#Redistribution and use in source and binary forms, with or without 
#modification, are permitted provided that the following conditions 
#are met:
#
#Redistributions of source code must retain the above copyright notice, 
#this list of conditions and the following disclaimer. 
#
#Redistributions in binary form must reproduce the above copyright 
#notice, this list of conditions and the following disclaimer in the 
#documentation and/or other materials provided with the distribution. 
#
#Neither the name of the Met Office nor the names of its contributors 
#may be used to endorse or promote products derived from this software 
#without specific prior written permission. 
#
#THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
#"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
#LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
#A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
#HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
#SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
#LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
#DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
#THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
#(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
#OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#########################################################

# Make a time-series from gridded input file
# Global average calculated as (NH+SH)/2

use strict;
use warnings;

# Extract the best-estimate area averaged time-series for the selected area
my @Date;
my @Best_estimate;
my @Best_estimate_field;
while ( my $Field = read_month( \*STDIN ) ) {
    if ( $Field->{month} == 1 ) { push @Date, $Field->{year}; }
    push @Best_estimate_field, $Field;
}

#average monthly fields to annual
@Best_estimate_field = field_to_annual(@Best_estimate_field);

#calculate area average for each year
for ( my $i = 0 ; $i < scalar(@Best_estimate_field) ; $i++ ) {
    $Best_estimate[$i] = area_average( $Best_estimate_field[$i] );
}

# Output the results
for ( my $i = 0 ; $i < scalar(@Date) ; $i++ ) {
    unless ( defined( $Best_estimate[$i] ) ) { print "\n"; next; }
    printf "%s %6.3f\n", $Date[$i], $Best_estimate[$i];
}

# Convert a monthly series of fields to an annual series
sub field_to_annual {

    my $Year;
    my @Count;
    my @Result;
    for ( my $i = 0 ; $i < scalar(@_) ; $i++ ) {
        my $Field = $_[$i];
        if ( $Field->{month} != 1 && !defined($Year) ) { next; }
        if ( $Field->{month} == 1 ) {
            $Year  = dclone($Field);
            @Count = undef;
            for ( my $i = 0 ; $i < $Year->{rows} ; $i++ ) {
                for ( my $j = 0 ; $j < $Year->{columns} ; $j++ ) {
                    if ( ${ $Year->{data} }[$i][$j] == $Year->{missing} ) {
                        next;
                    }
                    $Count[$i][$j] = 1;
                }
            }
        }
        else {
            for ( my $i = 0 ; $i < $Field->{rows} ; $i++ ) {
                for ( my $j = 0 ; $j < $Field->{columns} ; $j++ ) {
                    if ( ${ $Field->{data} }[$i][$j] == $Field->{missing} ) {
                        next;
                    }
                    if ( ${ $Year->{data} }[$i][$j] == $Year->{missing} ) {
                        ${ $Year->{data} }[$i][$j] =
                          ${ $Field->{data} }[$i][$j];
                        $Count[$i][$j] = 1;
                    }
                    else {
                        ${ $Year->{data} }[$i][$j] +=
                          ${ $Field->{data} }[$i][$j];
                        $Count[$i][$j]++;
                    }
                }
            }
            if ( $Field->{month} == 12 || $i == ( scalar(@_) - 1 ) ) {

                #  Normalise the data
                for ( my $i = 0 ; $i < $Year->{rows} ; $i++ ) {
                    for ( my $j = 0 ; $j < $Year->{columns} ; $j++ ) {
                        if ( ${ $Year->{data} }[$i][$j] == $Year->{missing} ) {
                            next;
                        }
                        ${ $Year->{data} }[$i][$j] /= $Count[$i][$j];
                    }
                }
                push @Result, $Year;
            }
        }

    }

    return @Result;
}

# Global average of field
sub area_average {

    my $Field = shift;

    # Calculate the area average
    my $W_tot_NH = 0;
    my $A_tot_NH = 0;
    my $W_tot_SH = 0;
    my $A_tot_SH = 0;
    for ( my $i = 0 ; $i < $Field->{rows} ; $i++ ) {

        #Grid-box area is proportional to cosine of latitude
        my $Lat = 90 - ( $i + 0.5 ) * 180 / $Field->{rows};
        if ( $Lat > 90 || $Lat < -90 ) { next; }
        my $Weight = cos( $Lat * 3.141592 / 180 );
        for ( my $j = 0 ; $j < $Field->{columns} ; $j++ ) {
            if ( $Field->{data}->[$i][$j] == $Field->{missing} ) { next; }
            if ( $Lat > 0 ) {
                $W_tot_NH += $Weight;
                $A_tot_NH += $Field->{data}->[$i][$j] * $Weight;
            }
            else {
                $W_tot_SH += $Weight;
                $A_tot_SH += $Field->{data}->[$i][$j] * $Weight;
            }

        }
    }

#if northern hemisphere or southern hemisphere average is undefined return nothing
    unless ( $W_tot_NH > 0 && $W_tot_SH > 0 ) { return; }

    #otherwise return the mean of the northern and southern hemisphere averages
    return ( $A_tot_NH / $W_tot_NH + $A_tot_SH / $W_tot_SH ) / 2;

}

# Read in the data for a month
sub read_month {

    my $fh = shift;    # FileHandle to read data from

    # Data structure to store gridded field
    my %data = (
        year    => undef,
        month   => undef,
        rows    => undef,    # Number of rows in the data grid
        columns => undef,    # Number of columns in the data grid
        missing => undef,    # What data value corresponds to missing data
        data    => undef,    # 2D array for gridded data
    );

    # One header line
    my $Header = <$fh>;
    unless ( defined($Header) ) {
        return;              # No more fields
    }
    unless ( $Header =~
        /^\s+(\d+)\s+(\d+)\s+\d+\s+(\d+) rows\s+(\d+) columns[\.\s]+Missing=([e\d\+\-\.]+)/
      )
    {
        die "Bad file format - unexpected header $Header";
    }
    $data{year}    = $1;
    $data{month}   = $2;
    $data{rows}    = $3;
    $data{columns} = $4;
    $data{missing} = $5;

    # Read in the data
    for ( my $i = 0 ; $i < $data{rows} ; $i++ ) {
        my $Row = <$fh>;
        unless ( defined($Row) ) { die "Unexpected end of data file"; }
        $Row =~ s/^\s+//;    # strip leading spaces - simplifies the split
        my @Columns = split /\s+/, $Row;
        unless ( scalar(@Columns) == $data{columns} ) {
            die
              "Unexpected data input - expected $data{columns} of data, got $Row";
        }
        $data{data}[$i] = [@Columns];
    }
    return \%data;
}

# Make a copy of a field
sub dclone {
    my $From = shift;
    my %To;
    foreach (qw(year month rows columns missing)) {
        $To{$_} = $From->{$_};
    }
    for ( my $i = 0 ; $i < $To{rows} ; $i++ ) {
        for ( my $j = 0 ; $j < $To{columns} ; $j++ ) {
            $To{data}[$i][$j] = $From->{data}->[$i][$j];
        }
    }
    return \%To;
}