<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 File::Slurp qw(write_file);
use Geo::BUFR;

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

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

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

# Make sure there are 2 input files
pod2usage(-verbose => 0) unless @ARGV == 2;

my $infile1 = $ARGV[0];
my $infile2 = $ARGV[1];

# 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 $bufr1 = Geo::BUFR->new();
$bufr1->fopen($infile1);
my $bufr2 = Geo::BUFR->new();
$bufr2->fopen($infile2);

my ($data1, $desc1) = $bufr1->next_observation();
my ($data2, $desc2) = $bufr2->next_observation();

die "Unexpanded descriptors in section 3 for $infile1 must be 309052, is "
    . $bufr1->get_descriptors_unexpanded()
    if $bufr1->get_descriptors_unexpanded() ne '309052';
die "Unexpanded descriptors in section 3 for $infile2 must be 309052, is "
    . $bufr2->get_descriptors_unexpanded()
    if $bufr2->get_descriptors_unexpanded() ne '309052';
die "Can only handle single subset BUFR messages, but $infile1 contains "
    . $bufr1->get_number_of_subsets() . " subsets"
    if $bufr1->get_number_of_subsets() != 1;
die "Can only handle single subset BUFR messages, but $infile2 contains "
    . $bufr2->get_number_of_subsets() . " subsets"
    if $bufr2->get_number_of_subsets() != 1;

my $bufr3 = Geo::BUFR->new();
$bufr3->copy_from($bufr1,'metadata');

my @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]
    );
my @windshear_level_desc = (
    '004086',  # Long time period or displacement [s]
    '008042',  # Extended vertical sounding significance [flag table]
    '007004',  # Pressure [PA]
    '005015',  # Latitude displacement (high accuracy) [degree]
    '006015',  # Longitude displacement (high accuracy) [degree]
    '011061',  # Absolute wind shear in 1 km layer below [m/s]
    '011062',  # Absolute wind shear in 1 km layer above [m/s]
    );

my ($data, $desc);

# First handle section 4 up to where pressure level data starts
for my $ii (0 .. 27) {
    $desc->[$ii] = $desc1->[$ii];
    if (defined($data1->[$ii])) {
        $data->[$ii] = $data1->[$ii];
    } else {
        $data->[$ii] = $data2->[$ii];
    }
}

# Then handle pressure levels (and eddrf)
$desc->[28] = '031002';
$data->[28] = undef; # Will be filled with $eddrf later
my $eddrf = 0; # extended delayed descriptor replication factor
my ($ii,$i1,$i2) = (29,29,29);
while ($desc1->[$i1] ne '031001' || $desc2->[$i2] ne '031001') {
    push(@$desc, @pressure_level_desc);
    if ($desc1->[$i1] eq '031001') {
        # No more pressure levels in message 1, so fetch from message 2
        push(@$data, @$data2[$i2 .. $i2+9]);
        $i2 += 10;
    } elsif ($desc2->[$i2] eq '031001') {
        # No more pressure levels in message 2, so fetch from message 1
        push (@$data, @$data1[$i1 .. $i1+9]);
        $i1 += 10;
    } elsif ($data1->[$i1+2] eq $data2->[$i2+2]) {
        # Same pressure level, should be merged into one
        # evss = extended vertical sounding significance
        my $evss1 = $data1->[$i1+1];
        my $evss2 = $data2->[$i2+1];
        my $evss = combine_evss($evss1,$evss2);
        push(@$data, defined($data1->[$i1]) ? $data1->[$i1] : $data2->[$i2]);
        push(@$data, $evss);
        for (2 .. 9) {
            push(@$data, defined($data1->[$i1 + $_]) ? $data1->[$i1 + $_]
                 : $data2->[$i2 + $_]);
        }
        $i1 += 10;
        $i2 += 10;
    } elsif ($data1->[$i1+2] > $data2->[$i2+2]) {
        # pick the highest pressure, i.e. data1
        push(@$data, @$data1[$i1 .. $i1+9]);
        $i1 += 10;
    } else {
        # pick the highest pressure, i.e. data2
        push(@$data, @$data2[$i2 .. $i2+9]);
        $i2 += 10;
    }
    $eddrf++;
    $ii += 10;
}
$data->[28] = $eddrf;

# Handle windshear levels
my $iws = $ii;
$desc->[$iws] = '031001';
$data->[$iws] = undef; # Will be filled later
if ($data1->[$i1] == 0 && $data2->[$i2] == 0) {
    $data->[$iws] = 0; # No wind shear levels
} else {
    my $nws = 0; # Number of wind shear levels, to be calculated
    my $nws1 = $data1->[$i1++];
    my $nws2 = $data2->[$i2++];
    $data->[$ii++] = 0;
    while ($nws1 > 0 || $nws2 > 0) {
        push(@$desc, @windshear_level_desc);
        if ($nws1 == 0) {
            # Fetch from message 2
            push(@$data, @$data2[$i2 .. $i2+6]);
            $i2 += 7;
            $nws2--;
        } elsif ($nws2 == 0) {
            # Fetch from message 1
            push(@$data, @$data1[$i1 .. $i1+6]);
            $i1 += 7;
            $nws1--;
        } elsif ($data1->[$i1+2] eq $data2->[$i2+2]) {
            # see no reason why data1 should differ from data2 here,
            # so we use data1
            push(@$data, @$data1[$i1 .. $i1+6]);
            $i1 += 7;
            $i2 += 7;
            $nws2--;
        } elsif ($data1->[$i1+2] > $data2->[$i2+2]) {
            # pick the highest pressure, i.e. data1
            push(@$data, @$data1[$i1 .. $i1+6]);
            $i1 += 7;
            $nws1--;
        } else {
            # pick the highest pressure, i.e. data2
            push(@$data, @$data2[$i2 .. $i2+6]);
            $i2 += 7;
            $nws2--;
        }
        $nws++;
        $ii += 7;
    }
    $data->[$iws] = $nws;
}

my ($data_refs, $desc_refs);
# One subset only
$data_refs->[1] = $data;
$desc_refs->[1] = $desc;

# Print the encoded BUFR message
my $msg = $bufr3->encode_message($data_refs, $desc_refs);
if ($option{outfile}) {
    write_file($option{outfile}, {binmode => ':raw'}, $msg);
} else {
    print $msg;
}


## 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;
    }
}

=pod

=head1 SYNOPSIS

  mergebufrtemp.pl <bufr file 1> <bufr file 2>
      [--outfile <file>]
      [--tablepath <path to BUFR tables>]
      [--help]

=head1 DESCRIPTION

Will merge (the first) BUFR message in <bufr file 1> with (the first)
BUFR message in <bufr file 2> into one single BUFR message, printing
to STDOUT unless C<--outfile> is set. The 2 input BUFR messages are
assumed to be single subset BUFR TEMPs utilizing TM305092 (this is
checked) containing data for the same station and termin (this is not
checked). Metadata are fetched from input file 1, the same applies for
conflicting data in section 4 (e.g. if BUFR launch time differs)
except that input file 2 is used for data missing in input file 1.

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})
   --help           Display Usage and explain the options used. Almost
                    the same as consulting perldoc mergebufrtemp.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
mergebufrtemp.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>