#!/bin/perl

# This perl script decodes GRIB code. Well, to be more precise, it decodes
# enough GRIB to allow me to read MPI model grib fields. But it can be 
# easily extended to read a bit more.
#
# Copyright (C) 1996, 2002 W. M. Connolley
#
# $Author: wmc $:$Date: 2002/07/22 21:07:25 $:$Revision: 1.1 $
#
# 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.
# http://www.gnu.org/copyleft/gpl.html

# Setup options
eval "\$$1=\$2" while $ARGV[0] =~ /^(\w+)=(.*)/ && shift;

if ($#ARGV eq -1) {
  print "Usage: grib.pl options filenames-list\n";
#  print "Output: info to tty and data to input-file-name.out.field-number\n";
  print "Output: info to tty and data to grib-out.field-id.year.month\n";
  print "Options: M=month - only output given month\n";
  print "       : Y=year  - only output given year\n";
  print "       : P=path - put output here (relative or abs)\n";
  exit;
}

# Decode GRIB...
while ( $Filei=shift @ARGV ) {

open(IN,$Filei);
print "Reading from file: $Filei\n";

# Main loop; try to read "GRIB"
$N=0;
READING: while (read(IN,$in,4)) {

# We may have to skip padding 4-octets (why?)
$Skipped=0;
while ($in ne "GRIB") {
  $Skipped=+4;
  if ( !read(IN,$in,4) ) { print "End of file; read $N fields\n"; last READING};
};
print "Read \"GRIB\" OK after skipping $Skipped octets\n"; 
$N++;

# Read PBD
read(IN,$in,3) || die "Can't read PDB block length";
$in=pack("xa3",$in);
$n=vec($in,0,32);
print "Length of PDB (1) block: $n\n";
read(IN,$in,$n-3) || die "Can't read remaining $n-3 bytes of PDB";
# Check we have a GDB
if (vec($in,4,8) & 128) { print "We have a GDB: OK\n"}
else { die "No GDB present: don't understand this grib\n"};
# Check we have a BDB
if (vec($in,4,8) & 64) { die "We have a BDB: fail"}
else { print "BDB not present: OK\n"};
# Get variable type
$Var=vec($in,5,8);
print "Variable code: $Var\n";
# Get year and month
$Year=vec($in,9,8);
# This is MPI-specific!
$Century=vec($in,15,8)*256+vec($in,16,8);	
$Year=$Year+$Century*100;
$Month=vec($in,10,8);
print "Year, Month: $Year, $Month\n";
$Skip=0;
if ($M && !($M==$Month)) { $Skip=1 };
if ($Y && !($Y==$Year)) { $Skip=1 };

# Read GDB
read(IN,$in,3) || die "Can't read GDB block length";
$in=pack("xa3",$in);
$n=vec($in,0,32);
print "Length of GDB (2) block: $n\n";
read(IN,$in,$n-3) || die "Can't read remaining $n-3 bytes of GDB";

# Interpret GDB info
$grid_type=vec($in,2,8);
if ($grid_type == 4) { print "Grid type is gaussian lat/regular long: OK\n"}
else { die "Don't recognize grid type: $grid_type"};
$nlat=vec($in,3,8)*256+vec($in,4,8);
$nlon=vec($in,5,8)*256+vec($in,6,8);
print "Nlat: $nlat, Nlon: $nlon\n";
$zlat=(vec($in,7,8)*256*256+vec($in,8,8)*256+vec($in,9,8))/1000;
$zlon=(vec($in,10,8)*256*256+vec($in,11,8)*256+vec($in,12,8))/1000;
print "Origin (lat, lon) [may not be coded]: $zlat, $zlon\n";

# Read DataB
read(IN,$in,3) || die "Can't read DataB block length";
$in=pack("xa3",$in);
$n=vec($in,0,32);
print "Length of DataB (2) block: $n\n";
read(IN,$in,$n-3) || die "Can't read remaining $n-3 bytes of DataB";

# Decode DataB
$scale_factor=2**(-(vec($in,1,8)*256+vec($in,2,8)-2**15));
print "Scale factor: $scale_factor\n";
$ref_frac=vec($in,4,8)*256*256+vec($in,5,8)*256+vec($in,6,8);
$ref_mant=vec($in,3,8);
if ($ref_mant & 128) {$ref_sign=-1} else {$ref_sign=1};
$ref_mant=($ref_mant & 127) - 64;
$ref_value=$ref_sign*$ref_frac*2**(-24)*16**($ref_mant);
print "Ref value: $ref_value\n";

$bits_per_value=vec($in,7,8);
print "Bits per packed value: $bits_per_value\n";
$spare=(vec($in,0,8) & 15);
if ($n == $nlat*$nlon*$bits_per_value/8+11+$spare/8) {
  print "OK: length of DataB equals lats*longs*bits/8+11+spare/8\n"
} else {
  die "DataB length not consistent with data it is supposed to contain"
};
if (!$Skip) {
  @data=unpack("x8n*",$in);
  foreach $x (@data) { $x=$x*$scale_factor+$ref_value };
  # $Fileo="$Filei.out.$N";
  $Fileo="Grib-out.$Var.$Year.$Month";
  if ($P) { $Fileo=$P."/".$Fileo };
  print "Output to $Fileo\n";
  open(OUT,"> ".$Fileo);
  foreach $x (@data) { print OUT "$x \n" };
  close OUT;
} else {
  if ($M) { print "Month $Month does not match requested month $M\n" };
  if ($Y) { print "Year $Year does not match requested year $Y\n" }
};

# Read "End"
read(IN,$in,4) || die "Can't read End";
if ($in eq "7777") { print "Read End OK\n"}
else { die "Didn't find 7777 at End: found: $in"};

}

print "Read $N fields from file $Filei\n";

}