<code perl>
#!/usr/bin/perl -w

# (C) Copyright 2015, met.no
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
# 02110-1301, USA.

# pod included at end of file

use strict;
use Getopt::Long;
use Pod::Usage qw(pod2usage);
use Geo::BUFR;

# Will be used if neither --tablepath nor $ENV{BUFR_TABLES} is set
use constant DEFAULT_TABLE_PATH => '/usr/local/lib/bufrtables';

our @pressure_level_desc = (
    '004086',  # Long time period or displacement [s]
    '008042',  # Extended vertical sounding significance [flag table]
    '007004',  # Pressure [PA]
    '010009',  # Geopotential height [GPM]
    '005015',  # Latitude displacement (high accuracy) [degree]
    '006015',  # Longitude displacement (high accuracy) [degree]
    '012101',  # Temperature/dry-bulb temperature [K]
    '012103',  # Dew-point temperature [K]
    '011001',  # Wind direction [degree true]
    '011002',  # Wind speed [m/s]
    );

# Parse command line options
our %option = ();
GetOptions(
           \%option,
           'help',
           'outfile=s',
           'tablepath=s',
           'verbose=i',
       ) or pod2usage(-verbose => 0);

# User asked for help
pod2usage(-verbose => 1) if $option{help};

# Make sure there is 1 input file
pod2usage(-verbose => 0) unless @ARGV == 1;

our $verbose = $option{verbose} ? $option{verbose} : 0;

# Set BUFR table path
if ($option{tablepath}) {
    # Command line option --tablepath overrides all
    Geo::BUFR->set_tablepath($option{tablepath});
} elsif ($ENV{BUFR_TABLES}) {
    # If no --tablepath option, use the BUFR_TABLES environment variable
    Geo::BUFR->set_tablepath($ENV{BUFR_TABLES});
} else {
    # If all else fails, use the libbufr bufrtables
    Geo::BUFR->set_tablepath(DEFAULT_TABLE_PATH);
}

my $infile = $ARGV[0];
open(my $IN, '<',$infile)
    or die "Cannot open $infile: $!";

# Where to print the altered BUFR message(s)
my $OUT;
if ($option{outfile}) {
    open($OUT, '>', $option{outfile})
        or die "Cannot open $option{outfile} for writing: $!";
} else {
    $OUT = *STDOUT;
}

my $bufr = Geo::BUFR->new();
$bufr->fopen($infile);
my @subset_data; # Will contain data values for subset 1,2...
my @subset_desc; # Will contain the set of descriptors for subset 1,2...
my $new_bufr;

READLOOP:
while (not $bufr->eof()) {
    # Read (and decode) next observation, skipping messages that
    # cannot be decoded
    my ($data, $desc);
        eval {
            ($data, $desc) = $bufr->next_observation();
        };
        if ($@) {
            warn $@;
            # Try to extract message number and ahl of the bulletin
            # where the error occurred
            my $current_message_number = $bufr->get_current_message_number();
            if (defined $current_message_number) {
                my $error_msg = "In message $current_message_number";
                my $current_ahl = $bufr->get_current_ahl();
                $error_msg .= " contained in bulletin with ahl $current_ahl\n"
                    if $current_ahl;
                warn $error_msg if $error_msg;
            }
            next READLOOP;
        }
    next READLOOP if !$data;

    my $isub = $bufr->get_current_subset_number();
    my $nsubsets = $bufr->get_number_of_subsets();
    my $imsg = $bufr->get_current_message_number();
    my $ahl = $bufr->get_current_ahl();
    if ($isub == 1) {
        $new_bufr = Geo::BUFR->new();
        $new_bufr->copy_from($bufr,'metadata');
        @subset_data = ();
        @subset_desc = ();
    }

    my $descriptors_unexpanded = $bufr->get_descriptors_unexpanded();
    if ($descriptors_unexpanded !~ /^30905(2|3)/) {
	if ($verbose > 0 && $isub == 1) {
	    print "Removing message $imsg";
	    print " contained in bulletin with ahl $ahl" if $ahl;
	    print " since unexpanded descriptors did not start with 309052 or 309053\n";
	}
	next;
    } else {
	if ($verbose > 3) {
	    print "Processing subset $isub in message $imsg";
	    print " contained in bulletin with ahl $ahl" if $ahl;
	    print "\n";
	}
        my @plevels; # array of references to the 10 data values of
                     # one pressure level
        my $i_eddrf = ($descriptors_unexpanded =~ /^309052/) ? 29 : 19;
	my $i1 = $i_eddrf;
        while ($desc->[$i1] ne '031001') {
            my @record = @$data[$i1 .. $i1+9];
            if (defined $record[2]) {
                push @plevels, \@record;
            } elsif ($verbose > 2) {
                print "Missing pressure in level ", ($i1 - $i_eddrf)/10 + 1
                    . " in subset $isub in Message number $imsg";
                print " in $ahl" if $ahl;
                print "\n";
            }
            $i1 += 10;
        }
        if (!@plevels || already_sorted(\@plevels)) {
            # Reencode the BUFR message with no changes
            $subset_data[$isub] = $data;
            $subset_desc[$isub] = $desc;
        } else {
            if ($verbose > 0) {
                print "Need to sort (or merge) subset $isub"
                    . " in Message number $imsg";
                print " in $ahl" if $ahl;
                print "\n";
            }
            my ($new_data, $new_desc)
                = get_sorted_data($data,$desc,\@plevels,$i1,$i_eddrf);
            $subset_data[$isub] = $new_data;
            $subset_desc[$isub] = $new_desc;
        }
    }

    if ($isub == $nsubsets) {
        print $OUT $new_bufr->encode_message(\@subset_data, \@subset_desc);
    }
}
$bufr->fclose();

## Will return in ($new_data, $new_desc) the same data as in
## ($data,$desc) except that pressure levels are sorted on decreasing
## pressure values, also merging levels with same pressure (if time
## and or position coordinates are equal in the 2 pressure levels or
## missing)
sub get_sorted_data {
    my ($data,$desc,$plevels,$i1,$i_eddrf) = @_;
    my ($new_data,$new_desc);

    # Section 4 up to where pressure level data starts is unaltered
    for my $ii (0 .. $i_eddrf-2) {
        $new_desc->[$ii] = $desc->[$ii];
        $new_data->[$ii] = $data->[$ii];
    }

    # Then handle pressure levels (and eddrf)
    $new_desc->[$i_eddrf-1] = '031002';
    $new_data->[$i_eddrf-1] = undef; # Will be filled with $eddrf later
    my $eddrf = 0; # extended delayed descriptor replication factor

    # Use a Schwartzian transform to sort on descending pressure levels
    my @sorted_plevels = map { $_->[0] }
                         sort { $b->[1] <=> $a->[1] }
                         map { [$_, $_->[2]] } @{$plevels};

    # Levels with same pressure and same time and position should be
    # coalesced into one
    my @prev_plevel = @{$sorted_plevels[0]};
    my $prev_PP = $prev_plevel[2];
    for my $ii (1 .. (@{$plevels}-1))  {
        my $PP = $sorted_plevels[$ii]->[2];
        if ($PP == $prev_PP
            && same_pos_time(\@prev_plevel,$sorted_plevels[$ii])) {
            # Combine the two levels into 1
            if ($verbose > 2) {
                print "Merging pressure level for pressure = $PP\n";
            }
            my @new_plevel;
            $new_plevel[0] = defined($prev_plevel[0])
                ? $prev_plevel[0] : $sorted_plevels[$ii]->[0];
            my $evss = combine_evss
                ($sorted_plevels[$ii-1]->[1], $sorted_plevels[$ii]->[1]);
            $new_plevel[1] = $evss;
            $new_plevel[2] = $PP;
            for (3 .. 9) {
                $new_plevel[$_] = defined($prev_plevel[$_])
                    ? $prev_plevel[$_] : $sorted_plevels[$ii]->[$_];
            }
            @prev_plevel = @new_plevel;
        } else {
            push @$new_data, @prev_plevel;
            push @$new_desc, @pressure_level_desc;
            @prev_plevel = @{$sorted_plevels[$ii]};
            $eddrf++;
        }
        $prev_PP = $PP;
    }
    # Add last level
    push @$new_data, @prev_plevel;
    push @$new_desc, @pressure_level_desc;
    $eddrf++;

    $new_data->[$i_eddrf-1] = $eddrf;

    # Section 4 after pressure level data is unaltered
    # (assuming windshear levels are sorted already)
    push @$new_data, @$data[$i1 .. @$data-1];
    push @$new_desc, @$desc[$i1 .. @$data-1];

    return ($new_data, $new_desc);
}

## Combine two extended vertical sounding significance (evss) into one
sub combine_evss {
    my ($evss1,$evss2) = @_;
    if (defined($evss1) && defined($evss2)) {
        my $binary1 = pack "N", $evss1; # Packed as 32 bits in big-endian order
        my $binary2 = pack "N", $evss2;
        my $binary = $binary1 | $binary2; # bitwise or
        return unpack "N", $binary;
    } elsif (defined($evss1)) {
        return $evss1;
    } else {
        return $evss2;
    }
}

## Return true (1) if pressure levels are sorted and unique, else
## return false (0)
sub already_sorted {
    my $plevels = shift;

    return 1 if !@{$plevels};
    my @PP = map { $_->[2] } @{$plevels};
    my $prev_PP = 1000000000; # Surely bigger than any real PP
    my $ii = 0;
    foreach my $PP (@PP) {
	if ($PP > $prev_PP) {
	    if ($verbose > 1) {
		print "Not sorted after pressure level $prev_PP";
		print " (next level is $PP)\n";
	    }
	    return 0;
	}
        if ($PP == $prev_PP) {
            if (same_pos_time($plevels->[$ii],$plevels->[$ii-1])) {
		if ($verbose > 1) {
		    print "Duplicate levels at pressure $PP\n";
		}
		return 0;
	    }
	}
        $prev_PP = $PP;
        $ii++;
    }
    return 1;
}

## Return false (0) if the two (references) to pressure levels have
## position or time defined with differing values, else returns true
## (1)
sub same_pos_time {
    my ($p1,$p2) = @_;

    if (defined($p1->[0]) && defined($p2->[0]) && $p1->[0] != $p2->[0]) {
        return 0; # 004086 'Long time period or displacement' differ
    } elsif (defined($p1->[4]) && defined($p2->[4]) && $p1->[4] != $p2->[4]) {
        return 0; # 005015 'Latitude displacement' differ
    } elsif (defined($p1->[5]) && defined($p2->[5]) && $p1->[5] != $p2->[5]) {
        return 0; # 005015 'Latitude displacement' differ
    }
    return 1;
}

=pod

=head1 SYNOPSIS

  sortbufrtemp.pl <bufr file>
      [--outfile <sorted bufr file>]
      [--tablepath <path to BUFR tables>]
      [--verbose n]
      [--help]

=head1 DESCRIPTION

Will extract from <bufr file> all BUFR TEMPs (recognized by unexpanded
descriptors in section 3 starting with 309052 or 309053), sorting the
pressure levels and merging duplicate pressure levels (having same
pressure, time and location coordinates) into one single level,
printing to STDOUT unless C<--outfile> is set. Pressure levels with
missing pressure will be removed.

Execute without arguments for Usage, with option C<--help> for some
additional info.

=head1 OPTIONS


   --outfile <filename>
                    Will print to <filename> instead of STDOUT
   --tablepath <path to BUFR tables>
                    Set path to BUFR tables (overrides $ENV{BUFR_TABLES})
   --verbose n      Set verbose level to n, 0<=n<=3 (default 0).
                    n>=1: prints which BUFR messages were actually altered,
                          and which messages were removed because unexpanded
                          descriptors did not start with 309052 or 309053
                    n>=2: prints at which pressure level sorting
                          breaks, or first duplicate level encountered
                          (whichever of these conditions happened
                          first)
                    n>=3: prints which pressure levels were merged,
                          and which levels were removed due to missing
                          pressure
                    n=4: prints which subset, message number and ahl
                         (if exists) is currently processed
                    Verbose output is sent to STDOUT, so ought to be
                    combined with option --outfile
   --help           Display Usage and explain the options used. Almost
                    the same as consulting perldoc sortbufrtemp.pl

Options may be abbreviated, e.g. C<--h> or C<-h> for C<--help>.

To avoid having to use the C<--tablepath> option, you are adviced to
set the environment variable BUFR_TABLES to the directory where your
BUFR tables are located (unless the default path provided by
sortbufrtemp.pl works for you).

=head1 AUTHOR

Pål Sannes E<lt>pal.sannes@met.noE<gt>

=head1 COPYRIGHT

Copyright (C) 2015 met.no

=cut
</code>