# Perl module providing methods for producing normals for marine data
# This module only provides class methods
# Hacked from the MDS2 version to add sea-ice normals (PB, 2008-06-18).

package Normals2;
use strict;
use Carp;
use PP;
use Exporter;
@Normals2::ISA    = qw(Exporter);
@Normals2::EXPORT = qw(daily_normal get_normals superob_normal);

sub get_normals {
    my $Parameter = shift;
    unless ( $Parameter =~ /SST|DAT|NAT|AST|PRE|LAT|ICE/ ) {
        return;
    }    # Use zero normals where not available
    my $Filename = get_normals_file_name($Parameter);
    my @Normals;
    open( DIN_N, "$Filename" ) or die "Unable to open normals file $Filename";
    for ( my $i = 0 ; $i < 73 ; $i++ ) {
        $Normals[$i] = pp_read( \*DIN_N ) or die "Failed to read in normals";
    }
    close(DIN_N);
    return \@Normals;
}

# Get normals file name
sub get_normals_file_name {
    my $Parameter = shift;
    unless ( $Parameter =~ /SST|DAT|NAT|AST|PRE|LAT|ICE/ ) {
        croak
          "Argument to get_default_normals_file must be one of SST, MAT, AST, LAT, ICE and pressure";
    }

 #   unless(defined($ENV{MDS2})) { croak "Environment variable MDS2 not set."; }
    my $Normals_file_name;
    if ( $Parameter =~ /SST/ ) {
        $Normals_file_name =
          "/ibackup/cr1/hadobs/OBS/marine/HadSST2/norms/HadSST2_pn1d6190.pp";
    }
    elsif ( $Parameter =~ /DAT/ ) {
        $Normals_file_name =
          "/ibackup/cr1/hadobs/OBS/marine/MOHMAT/norms/MOHMATD41_pn1dg6190.pp";
    }
    elsif ( $Parameter =~ /NAT/ ) {

#      $Normals_file_name = "/ibackup/cr1/hadobs/OBS/marine/MOHMAT/norms/MOHMATN4_pn1dg6190.pp";
        $Normals_file_name =
          "/ibackup/cr1/hadobs/OBS/marine/HadNAT2/norms/HadNAT2_climatology_pn1d.pp";
    }
    elsif ( $Parameter =~ /AST/ ) {
        $Normals_file_name =
          "/ibackup/cr1/hadobs/OBS/marine/MOHAST/norms/MOHAST_pn1dg6190.pp";
    }
    elsif ( $Parameter =~ /PRE/ ) {
        $Normals_file_name =
          "/ibackup/cr1/hadobs/OBS/pressure/HADSLP2/norms/gmslp3_norms_6190_90n90s_121_pentad_1x1.pp";
    }
    elsif ( $Parameter =~ /LAT/ ) {
        $Normals_file_name =
          "/data/cr2/hadpb/CRU_climatology/data/temp/pentad_1x1_tmp.pp";
    }
    elsif ( $Parameter =~ /ICE/ ) {
        $Normals_file_name =
          "/data/cr2/hadpb/ICE_climatology/ice_clim_hadISST_pentadInt_1870.pp";
    }
    else {
        croak "Unsupported value for parameter";
    }
    return $Normals_file_name;
}

# Return the normal for an ob interpolated to daily resolution
sub daily_normal {
    my ( $Normals, $longitude, $latitude, $Pentad, $Penday ) = @_;
    unless ( defined($Normals) ) { return 0; }    # Normals default to zero
    my ( $Pen1, $Pen2, $Weight1, $Weight2 );
    if ( $Penday < 3 ) {
        $Pen1 = $Pentad - 1;
        $Pen2 = $Pentad;
    }
    elsif ( $Pentad > 3 ) {
        $Pen1 = $Pentad;
        $Pen2 = $Pentad + 1;
    }
    else {
        $Pen1 = $Pentad;
        $Pen2 = $Pentad;
    }
    if ( $Pen1 < 1 )  { $Pen1 = 73; }
    if ( $Pen2 > 73 ) { $Pen2 = 1; }
    $Weight1 = ( 2, 1, 1, 4, 3 )[ $Penday - 1 ];
    $Weight2 = ( 3, 4, 1, 1, 2 )[ $Penday - 1 ];
    my $Lat_i   = lat_index($latitude);
    my $Lon_i   = lon_index($longitude);
    my $Normal1 = ${ $$Normals[ $Pen1 - 1 ]->{data} }[$Lat_i][$Lon_i];
    my $Pcount  = 0;

    while ( $Normal1 == $$Normals[ $Pen1 - 1 ]->{bmdi} && $Pcount < 8 ) {
        ( $Lat_i, $Lon_i ) =
          llPerturb( $Normals->[0], $Pcount, $Lat_i, $Lon_i );
        $Normal1 = ${ $$Normals[ $Pen1 - 1 ]->{data} }[$Lat_i][$Lon_i];
        $Pcount++;
#        warn "$Lat_i $Lon_i";
    }
    my $Normal2 =
      ${ $$Normals[ $Pen2 - 1 ]->{data} }[ $Lat_i ]
      [ $Lon_i ];
    if (   $Normal1 == $$Normals[ $Pen1 - 1 ]->{bmdi}
        || $Normal2 == $$Normals[ $Pen2 - 1 ]->{bmdi} )
    {
        return;
    }
    else {
        return ( $Normal1 * $Weight1 + $Normal2 * $Weight2 ) /
          ( $Weight1 + $Weight2 );
    }
}

# Perturb the latitude and longitude indices to give a
#  normal for an adjacent point if possible
sub llPerturb {
    my $PP     = shift;
    my $Pcount = shift;
    my $Lat_i  = shift;
    my $Lon_i  = shift;
    if ( $Pcount == 0 ) {
        $Lon_i += 1;
    }
    elsif ( $Pcount == 1 ) {
        $Lon_i -= 2;
    }
    elsif ( $Pcount == 2 ) {
        $Lon_i += 1;
        $Lat_i += 1;
    }
    elsif ( $Pcount == 3 ) {
        $Lat_i -= 2;
    }
    elsif ( $Pcount == 4 ) {
        $Lat_i += 1;
        $Lon_i += 2;
    }
    elsif ( $Pcount == 5 ) {
        $Lon_i -= 4;
    }
    elsif ( $Pcount == 6 ) {
        $Lon_i += 2;
        $Lat_i += 2;
    }
    elsif ( $Pcount == 7 ) {
        $Lat_i -= 4;
    }
    else {
        die "Exceeded pcount";
    }
    if ( $Lon_i < 0 ) { $Lon_i += $PP->{lbnpt}; }
    if ( $Lon_i >= $PP->{lbnpt} ) { $Lon_i -= $PP->{lbnpt}; }
    if ( $Lat_i < 0 )             { $Lat_i = 0; }
    if ( $Lat_i >= $PP->{lbrow} ) { $Lat_i = $PP->{lbrow}; }
    return ( $Lat_i, $Lon_i );
}

# Return the normal for a given time and place, uninterpolated.
sub superob_normal {
    my ( $Normals, $longitude, $latitude, $Pentad ) = @_;
    unless ( defined($Normals) ) { return 0; }    # Normals default to zero
    if ( ${ $$Normals[ $Pentad - 1 ]->{data} }[ lat_index($latitude) ]
        [ lon_index($longitude) ] == $$Normals[ $Pentad - 1 ]->{bmdi} )
    {
        return;
    }
    return ${ $$Normals[ $Pentad - 1 ]->{data} }[ lat_index($latitude) ]
      [ lon_index($longitude) ];
}

# Convert a longitude into an index for normal data
sub lon_index {
    my $longitude  = shift;
    my $long_local = $longitude;
    if ( $long_local == -180 ) { $long_local += 0.01; }
    if ( $long_local == 180 ) { $long_local -= 0.01; }
    if   ( $long_local > 0.0 ) { return int($long_local) + 180; }
    else                       { return int($long_local) + 179; }
}

# Convert a latitude into an index for normal data
sub lat_index {
    my $latitude  = shift;
    my $lat_local = $latitude;
    if ( $lat_local == -90 ) { $lat_local += 0.01; }
    if ( $lat_local == 90 ) { $lat_local -= 0.01; }
    if   ( $latitude > 0.0 ) { return 89 - int($lat_local); }
    else                     { return 90 - int($lat_local); }
}