#!/bin/perl 

# Decode ECMWF-style GRIB.
#
# Copyright (C) 1996, 2002 W. M. Connolley
#
# $Author: wmc $:$Date: 2002/07/22 21:10:03 $:$Revision: 1.9 $           
#
# 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
#
# Call as "wmc_grib.pl" for options list
#
# ECMWF seem to write EDITION=1 grib. I don't have the full specs for this...
# so I've had to guess.
# You might find https://wms.ecmwf.int/documents/manuals/libraries/gribex/
# helpful (you need a certificate to visit those pages)
#
# From the EMCWF manual:
#
#  Content   Grib-0  Grib-1
#  GRIB         1-4     1-4
#  Total l        *     5-7
#  Edition #      *       8
#  Sect1 l      5-7    9-11
#  Reserved (0)   8       *
#  Ver of ct2     *      12
#  Id of center   9      13
#  ...and from here on, grib-1 is at grib-0 +4 octets.
#  So, to tell if we're grib-0 or grib-1, we can
#  look at octet 8, I guess: this should be 0 in grib-0 and
#  1 in grib-1.
#
# If you encounter problems with this, you could look in "MARS user guide (for
# data retrieval) [also ECMWF computer bulletin B6.7/2] revision 11 - sept/95,
# appendix E. Also, try setting EDITION=1
#
# Even worse, there are *BUGS* in their complex packing of spherical harmonics
# (I think)
#
# To just extract specific fields and output just those unchanged, call as
#    wmc_grib.pl OUTFORM=grib selection-options filesnames,
# eg
#    wmc_grib.pl OUTFORM=grib V=129 Y=1991 M=1 LEV=500 filenames
# will just select 500 hPa height in january 1991.
#
# Notice: usual usage restrictions - no commmercial use, free academic
# use, this notice must be left intact and a similar condition imposed
# on any third party. No warranty, bugs very likely exist.
# W. M. Connolley August 1996 - Sept 2000
# wmc@bas.ac.uk / http://www.wmc.care4free.net/

# Restrictions: 
#  Only works on ECMWF "grib". 
#  Only decodes "ll" and "gp" files. 
#  HEADER stuff is UKMO pp-field specific, ignore it.

# Known bugs:
#  10 hPa geopotential in llpr93010100 causes floating exception.
#  The reading of the "scale factor" may be dodgy... search 
#    the code for string "scale_factor".

# Fixes:
#  WMC 27/08/97 - add PDBGUESS mode. Discover ECMWF accumulated-ppn fields require PDBGUESS=66
#  ditto        - fix (I think!) bug in guessing scale factor if exponent negative
#  WMC 28/08/97 - add EDITION=1 mode
#  WMC 11/09/97 - add option to include level type in filename
#  WMC 15/09/00 - add OUTMODE option of ">>". Rewrite the above header a bit to be less
#                 nasty to ECMWF who are nice really.
#  WMC 21/09/01 - Add option "VCP" to allow an offsete in the grid definition section when
#                 ECMWF add in the vertical coords when transferring from model to sfc grid.
#                 Fortunately the extra 122 reals are included in the length of the grid section
#                 so there is no need to do adjustments elsewhere.
#  WMC 01/10/01 - Note that VCP=416 for ERA-15 (I think)

$O1="/dev/null";		# Can be overridden as desired
$OUTFORM="txt";			# Default output in text format
$PDBGUESS="56";			# Default guess if not set to 24 in the grib. Could also try 66!
$OUTMODE=">";			# Default is to overwrite existing (non-avg) files.
				# Set this to append; useful if data from more than one
                                # timestep is present
$USETIME="";			# Set to "T1" to use forecast time
$FAKEAVG="";			# Set to, eg, ".avg"
$VCP=0;				# If you need to use this, a good guess is 488=122*4;
                                # (ref: mail from Sakari to GJM).  
                                # Or 416 (or 256?!?) (for ERA-15)

# UKMO Note - if "bin" is selected as the OUTFORM, then the resulting file
# could be used as the data block of a pp-file, ie it just needs the header
# cat'd on. This is actually done if option "HEADER" is set, whereupon
# that header is added to the output, with just lbyr, mon, yrd and mond
# modified from the header.

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

if ($#ARGV eq -1) {
  print "Usage: wmc_grib.pl options filenames-list\n";
  print "Output: info: tty (see O1), data: Grib-out.field-id.level.year.month.day.hour\n";
  print "Options: M=month - only output given month\n";
  print "       : Y=year  - only output given year\n";
  print "       : YG=year - only output ge given year\n";
  print "       : V=nnn - only output if variable is nnn (ECMWF table 2)\n";
  print "       : P=path - put output here (relative or abs)\n";
  print "       : O1=filename - mostly junk text out, defaults /dev/null\n";
  print "       : NTW=nn - only output the nth field in the file\n";
  print "       : OUTFORM=bin - output binary 4-byte floats, default txt\n";
  print "       :        =grib - output grib unchanged (but to files and can select)\n";
  print "       :        =bin2 - output binary data as in the data block";
  print "       : TRI=1 - output as lat lon value, default just value\n";
  print "       : AVG=1 - also output average of all written fields\n";
  print "       : AVG=2 - output average but not individual fields\n";
  print "       : LEV=x - only output level with value x. For llpr files, x is in hPa\n";
  print "       : HEADER=filename - prefix bin output with pp-header in filename\n";
  print "       : [removed: use EDITION=1] PDBGUESS=g or =nn - guess proper PDB block length, or use supplied value\n";
  print "       : EDITION=1 - set code edition to 1, appropriate for some ERA grib\n";
  print "       : LEVT_IN_OUTNAME=1 - include level type code in the filename output\n";
  print "       : OUTMODE=\">>\" - set to append\n";
  print "       : USETIME=T1 - use the forecast time T1 as a delta for the output file timestamp\n";
  print "       : USE_CALC_DB=1 - use the calculated length of the data block, not the specified one (see wmc's ERA web pages under 'lwd' if you end up needing this...)\n";
  print "       : DEBUG=1\n";
  exit;
}

if ($DEBUG) { $|=1 };

use Date::Manip;
Date_Init("DateFormat=nonUS");

# Read in the HEADER if required. Only read in the first 64*4 bytes - 
# this allows us to read from an actual pp-field. Neat, eh?
if ($OUTFORM eq "bin" && -r $HEADER) {
  open(HEAD,$HEADER) || die "Can't open header file";
  read(HEAD,$Header,64*4) || die "Can't read from header file";
  close(HEAD);
}

#
# Decode GRIB..., looping over files in arg list
#
while ( $Filei=shift @ARGV ) {

open(IN,$Filei);
print "$Filei -> ";

if (!open(OUT1,"> $O1")) {
 print "Can't open info output file - trying /dev/null"; 
 $O1="/dev/null";
 open(OUT1,"> $O1") || die "Can't open junk out to $O1" ;
};
print OUT1 "Reading from file: $Filei\n";

#
# Inner-Main loop; try to read "GRIB"
#
# Count input and output fields
$N=0; $NW=0;
READING: while (read(IN,$in,4)) {

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

#
# Read PBD
# This becomes a bit complex for ERA code which sets "code edition 1" and
# expects us to understand that. Ho hum.
#
$Save=$in;
if ($EDITION eq "1") {read(IN,$in1,4) or die "Can't read junk in PDB"; $Save.=$in1};
read(IN,$in,3) or die "Can't read PDB block length";
$Save.=$in;
$in=pack("xa3",$in);
$n=vec($in,0,32);
print OUT1 "Length of PDB (1) block: $n\n";
if ($EDITION ne "1" and $n > 500) { 
  die "Warning: length of PDB is $n which is rather large; try EDITION=1, perhaps" 
};

# OK, have length, read in the rest of the PDB
read(IN,$in,$n-3) || die "Can't read remaining $n-3 bytes of PDB";
$Save.=$in;
$in=$in1.$in;

# Check we have a GDB
if (vec($in,4,8) & 128) { print OUT1 "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 OUT1 "BDB not present: OK\n"; if ($DEBUG) { print "BDB not present: OK\n" } };

# Get variable type, level type and level value from PDB
$Var=vec($in,9,8);
$LevT=vec($in,10,8);
$LevV=vec($in,11,8)*256+vec($in,12,8);
print OUT1 "Variable code, Level type, value: $Var $LevT $LevV\n";

# Get year and month, etc
$Year=vec($in,13,8);
# This [may be] MPI-specific!
$Century=vec($in,25,8);
$Year=$Year+($Century-1)*100;
$Month=vec($in,14,8);
$Day=vec($in,15,8);
$Hour=vec($in,16,8);

# Get forecast stuff
$TimeRangeUnitInd=vec($in,18,8);
$Time1=vec($in,19,8);
$Time2=vec($in,20,8);
$TimeRangeInd=vec($in,21,8);
$Navg=vec($in,22,8)*256+vec($in,23,8);
$CT5{"113"}="Avg of $Navg forecasts/analyses";
$CT4{"1"}="minute"; $CT4{"1"}="hour"; $CT4{"2"}="day"; $CT4{"3"}="month"; $CT4{"4"}="year";
if (!($CT5=$CT5{$TimeRangeInd})) { $CT5=$TimeRangeInd };
if (!($CT4=$CT4{$TimeRangeUnitInd})) { $CT4=$TimeRangeUnitInd };
print OUT1 "T1: $Time1 T2: $Time2 TRI: $CT5 ($CT4; $TimeRangeInd) [$Year/$Month/$Day/$Hour]\n";

# If asked, we can reset the time used on the filename to be the forecast time...
if ($USETIME eq "T1") {
  if ($DEBUG) { print "Not using std time $Year/$Month/$Day $Hour but using forecast time " };
  $d=ParseDate("$Year/$Month/$Day $Hour:00") or die "Failed to ParseDate [$Year,$Month,$Day,$Hour]"; 
  $d1=DateCalc($d,"+".$Time1." $CT4") or die "@_";
  ($Year,$Month,$Day,$Hour)=split("/",&UnixDate($d1,"%Y/%m/%d/%H"));
  if ($DEBUG) { print "$Year/$Month/$Day $Hour\n" };
};

# Read some EMCWF-sepcific stuff
# See: https://wms.ecmwf.int/documents/manuals/libraries/gribex/localDefinition8.html
my $o41 = vec($in,41,8);
my $o42 = vec($in,42,8);
my $o43 = vec($in,43,8);
my $o44_45 = vec($in,44,16);
print OUT1 "[octets 41, 42, 43, 44-45] (8=era): $o41; (3=era): $o42; (2=an, 9=fc, 19=fa): $o43; (stream: 1043=mon mn; 1071=mon mn day mn): $o44_45\n";

# Decide whether we want this field or not...
$Skip=0;
if ($M && $M!=$Month) { $Skip=1; print OUT1 "Skipping due to month\n" };
if ($Y && $Y!=$Year) { $Skip=1; print OUT1 "Skipping due to year\n" };
if ($YG && $YG>$Year) { $Skip=1; print OUT1 "Skipping due to yearG\n" };
if ($NTW && $N!=$NTW) { $Skip=1; print OUT1 "Skipping due to # in file\n" };
if ($V && $V!=$Var) { $Skip=1; print OUT1 "Skipping due to var $Var ne $V\n"};
if ($LEV && $LevV!=$LEV) { $Skip=1; print OUT1 "Skipping due to level $LevV ne $LEV\n" };

$Text="V Y M D H: $Var, $Year, $Month, $Day, $Hour ";
if ($Skip) { print OUT1 "$Text ...skipping\n"; if ($DEBUG) { print "$Text ...skipping\n" } } 
      else { print "$Text -> "; $NW++ };

#
# Read GDB
#
read(IN,$in,3) || die "Can't read GDB block length";
$Save.=$in;
$in=pack("xa3",$in);
$n=vec($in,0,32);
print OUT1 "Length of GDB (2) block: $n\n"; if ($DEBUG) { print "Length of GDB (2) block: $n\n" };

read(IN,$in,$n-3) || die "Can't read remaining $n-3 bytes of GDB";
$Save.=$in;

#
# Interpret GDB info
#
$grid_type=vec($in,2,8);
if ($DEBUG) { print "Grid type: $grid_type\n" };
if ($grid_type == 0) { 
#
# Regular lat-lon grid
#
  print OUT1 "Grid type is regular lat/regular long: OK\n";
  $nlon=vec($in,3,8)*256+vec($in,4,8);
  $nlat=vec($in,5,8)*256+vec($in,6,8);
  $totalpoints=$nlat*$nlon;
  print OUT1 "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 OUT1 "Origin (lat, lon) [may not be coded]: $zlat, $zlon\n";
  $dlat=-(vec($in,20,8)*256+vec($in,21,8))/1000;
  $dlon=(vec($in,22,8)*256+vec($in,23,8))/1000;
  print OUT1 "Delta (lat, lon) [may not be coded]: $dlat, $dlon\n"; 

} elsif ($grid_type == 4) {
#
# Gaussian lat-lon grid.
# Complication: may be regular (eg MPI) or reduced (eg ECMWF)
# For ECMWF-ERA, look at nlat to decide: if its 65535 (-1 by another name)
# then we have a reduced grid.
#
  $nlat=vec($in,3,8)*256+vec($in,4,8);
  $reduced=($nlat==65535);
  $nlon=vec($in,5,8)*256+vec($in,6,8);
  print OUT1 "Grid type is Gaussian lat/";
  undef @NLONS;
  if ($reduced) {
# Reduced grid. Complicated.
    print OUT1 "reduced lon; OK\n"; 
    $npoletoeq=vec($in,23,8);
    print OUT1 "Points from pole to equator: $npoletoeq\n"; if ($DEBUG) { print "reduced lon; pte: $npoletoeq\n" };
    print OUT1 "Points along latitudes:\n";
    $totalpoints=0;
    for ($i=0; $i<$npoletoeq; $i++) {
      $j=vec($in,29+$i*2+$VCP,8)*256+vec($in,30+$i*2+$VCP,8);
      print OUT1 "$j "; 
      $NLONS[$i]=$j;
      $totalpoints+=$j*2;
    }; print OUT1 "\n"; if ($DEBUG) { print " NL: ",@NLONS+0,"\n" };
    @NLONS=(@NLONS,reverse(@NLONS));
  } else {
# Regular grid. Simple.
    print OUT1 "regular lon; OK\n";
    $totalpoints=$nlat*$nlon;
  };

  print OUT1 "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 OUT1 "Origin (lat, lon) [may not be coded]: $zlat, $zlon\n"; if ($DEBUG) { print "$zlat, $zlon\n" };

} elsif ($grid_type == 50) {
#
# Spherical harmonic grid. Ohh-err.
#
  print OUT1 "Spherical harmonic coefficients grid (grid type 50)\n";
  $in=" ".$in;		# Add leading byte to allow easy use of vec
  $Jpar=vec($in,2,16);
  $Kpar=vec($in,3,16);
  $Mpar=vec($in,4,16);
  print OUT1 "J, K, M: $Jpar, $Kpar, $Mpar\n";
  if ($Jpar==$Kpar && $Jpar==$Mpar) {
    print OUT1 "Triangular representation: OK\n";
# Total points is 0...J triangle, *2 for real+complex
# Note that this is modified by "complex packing" in DataB
    $totalpoints=($Jpar+1)*($Jpar+2)
  } else {
    die "Spherical harmonics not triangular. Dying"
  };
  $RepType=vec($in,10,8); $RepMode=vec($in,11,8);
  print OUT1 "Rep type and mode: $RepType, $RepMode\n";
  $Lat=(vec($in,29,8)*256+vec($in,30,8))*256+vec($in,31,8);
  print OUT1 "S Pole Lat: $Lat\n";

} else { 
#
# OK, we don't understand the grid type. This is OK if skip is set, or if
# OUTFORM is "grib" in which case we don't have to understand it!
#
  if ($OUTFORM eq "grib" || $Skip) {
    print OUT1 "Don't recognize grid type: $grid_type - OK because not interpreting\n"
  } else {
    die "Don't recognize grid type: $grid_type" 
  }
};

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

# Decode DataB
$flag=vec($in,0,8);
# Decode scale factor. Check this if weird things happen in the output.
# The first form of getting the scale factor may be machine independent. The
# second form needs "n" for alpha-like machines and "s" for sun-like ones.
# Appears to fail for negative exponents: if (substr($in,1,1)&0x80) {$pm=-1} else {$pm=1};
if (vec($in,1,8)&0x80) {$pm=-1} else {$pm=1};
$scale_factor=2**($pm*( (vec($in,1,8)&0x7F)*256+vec($in,2,8) ));
# $scale_factor=2**(unpack("n",substr($in,1,2)));
print OUT1 "Scale factor: $scale_factor\n";
print OUT1 "          is: ",($pm*( (vec($in,1,8)&0x7F)*256+vec($in,2,8) )),"\n";
print OUT1 "          is: ",vec($in,1,8)," ",vec($in,2,8)," (sign: $pm)\n";
$ref_value=&get_ref_val(substr($in,3,4));
print OUT1 "Ref value: $ref_value\n";

$bits_per_value=vec($in,7,8);
print OUT1 "Bits per packed value: $bits_per_value";
if ($bits_per_value % 8 != 0) { 
  print "\n\n***argh! bits_per_value = $bits_per_value\n\n***I'm going to fail, later...\n\n" ;
  print OUT1 " ...this is **not OK**\n"
} else {
  print OUT1 " ...OK\n"
};

$spare=(vec($in,0,8) & 15);
if ($grid_type != 50) {
  $n1=$totalpoints*$bits_per_value/8+11+$spare/8
} else {
  if ($flag & 64) {
    print OUT1 "Spectral packing, complicated: OK\n"
  } else {
    die "Spectral packing simple mode: not understood"
  };
  $xN=vec($in,4,16); print OUT1 "N: $xN\n";
  $J1=vec($in,12,8); print OUT1 "J1: $J1\n";
  $K1=vec($in,13,8); print OUT1 "K1: $K1\n";
  $M1=vec($in,14,8); print OUT1 "M1: $M1\n";
# Check that N is consistent with J1, etc
# This seems to come out wrong at the moment - odd
  if ($xN-18-($J1+1)*($J1+2)*2*2 == 0) {
    print OUT1 "Number of complicated packed values OK\n"
  } else {
    print OUT1 "Wrong number? of complicated packed values ($xN, ",($J1+1)*($J1+2)*2*2+18,")" 
  };
  $xP=vec($in,5,16); print OUT1 "P: $xP\n";
  $n1=$totalpoints*$bits_per_value/8+18+($J1+1)*($J1+2)*2
};
$db_length="*"; 
if ($n == $n1) {
  print OUT1 "Length of DataB is OK\n"
} else {
  unless ($OUTFORM eq "grib" || $Skip) {
    if ($USE_CALC_DB) {
# If we are here, its because we have set USE_CALC_DB, so we probably don't want
# to see an error message. Put it to the O1 file instead.
      print OUT1 "DataB length ($n) not consistent with data it is supposed to contain ($n1) ...Going with calculated version ($n1)";
      $db_length=$totalpoints;
    } else {
      warn "DataB length ($n) not consistent with data it is supposed to contain ($n1) ...Going with DataB version ($n)";
    }
  }
};

# Decode and output data, unless "Skip" set
if (!$Skip) {
# Don't decode if "grib" output
  if ($OUTFORM ne "grib") {

    if ($grid_type eq 50) {
# Harmonics are complicated to unpack

# Define data array length
      $big=($Jpar+1)*($Jpar+2)/2*2;
      $small=($J1+1)*($J1+2)/2;
      print OUT1 "Big, Small: $big, $small\n";
      $#data=$big;
 
# First we have to unpack (J1+1)*(J1+2)/2 real-imag pairs coded like the
# reference value (ie, unpacked). 
# Due to a bug in ECMWF encoding, we have to scale the last in each row
      $#data=$big-1;
      $pt1=0; $pt=0;
      for ($I=0;$I<=$J1*2+1;$I+=2) {
        for ($J=$I;$J<=$J1*2+1;$J++) {
          $data[$pt+$J]=&get_ref_val(substr($in,15+4*$pt1++,4));
        };
        $pt+=$Jpar*2-$I;
      };

# Now the rest of the truncated triangle, packed
      $pt0=$pt1; $pt=0;
      @data1=unpack("n*",substr($in,15+$small*8,($big-$small)*4));
      for ($I=0;$I<=$Jpar*2+1;$I+=2) {
        for ($J=0;$J<=$Jpar*2+1-$I;$J++) {
          if ($J>=($J1*2+2-$I)) {
            if (!$J || $J==($J1*2+2-$I)) { $Ss=1 } else {
              $Z=int(($J+$I)/2);
              $Ss=1.0/($Z*($Z+1))**($xP/1000.)
            };
            $data[$pt+$J]=($data1[$pt1++ -$pt0]*$scale_factor+$ref_value)*$Ss;
          };
        };
        $pt+=$Jpar*2+2-$I;
      };

# Now correct ECMWFs mis-coding
      $pt=0;
      for ($I=0;$I<=$J1;$I++) {
        if ($I<=$J1) {$Z=$J1} else {$Z=$I};
        $Ss=1.0/($Z*($Z+1))**($xP/1000.); 
        $data[2*($pt+10)]*=$Ss;
        $data[2*($pt+10)+1]*=$Ss;
        if ($I<=$J1) {$Z=$J1+1} else {$Z=$I};
        $Ss=1.0/($Z*($Z+1))**($xP/1000.);
        $data[2*($pt+10)+2]*=$Ss;
        $pt+=($Jpar-$I);
      };
      for ($I=$J1+1;$I<=$Jpar;$I++) {
        $Ss=1.0/($I*($I+1))**($xP/1000.);
        $data[2*($pt+10)+2]*=$Ss;
        $pt+=($Jpar-$I+1);
      };

# Check that $pt1 now equals $big... I hope so!
      if ($pt1 == $big) {
        print OUT1 "pt equals big: counted points OK\n"
      } else {
        die "Counted wrong: pt1 ne big: $pt1, $big"
      } 

    } else {
# Normal type (not spectral) packing. Easy

# This was the old code. It is valid, as long as bits-per-packed-value is 16.
# Sadly, this is not always true: esp after "compute" at ECMWF. So... lets make it 
# a bit more general
#      @data=unpack("x8n$db_length",$in);
#      foreach $x (@data) { $x=$x*$scale_factor+$ref_value };
# OK, this is more general, but probably hideously inefficient...
      my $nb=$bits_per_value/8.;
      for (my $i=0; $i<$totalpoints; $i++) {
        $data[$i]=0;
        my $s = 256;
        for (my $j=0; $j<$nb; $j++) {
          $data[$i]=($data[$i]*$s)+vec($in,$nb*$i+$j+8,8);
        };
        $data[$i]=$data[$i]*$scale_factor+$ref_value;
      }; 

    };

    if ($AVG) { 
      if ($#data_avg+1) { $i=0; foreach $x (@data) { $data_avg[$i++]+=$x } }
      else { @data_avg=@data }
    }
  };
  if ($LEVT_IN_OUTNAME) {
    $Fileo="Grib-out.$Var.$LevV.$LevT.$Year.$Month.$Day.$Hour$FAKEAVG";
  } else {
    $Fileo="Grib-out.$Var.$LevV.$Year.$Month.$Day.$Hour$FAKEAVG";
  };
  if ($P) { $Fileo=$P."/".$Fileo };
  print "$Fileo";
  print OUT1 "First/last data value: $data[0]/$data[scalar(@data)-1] ".(scalar(@data))."\n";
# If AVG=2, skip data output for individual fields.
  if ($AVG!=2) { 
# OUTMODE is usually ">" but can be set to ">>"
    open(OUT,$OUTMODE.$Fileo);
# grib output - output records unchanged. Useful to split up files
# and to select only certain fields.
    if ($OUTFORM eq "grib" ) {
      print OUT $Save."7777";
    } elsif ($OUTFORM eq "txt") {
# Text output format one value per line
      $i=-1;
      if ($reduced) { $i0=0; $i1=0 };
      foreach $x (@data) { 
        $i++;
        if ($TRI) { 
          if (!$reduced) {
            print OUT $zlat+$dlat*int($i/$nlon)," ",$zlon+($i % $nlon)*$dlon," $x\n" 
          } else {        
            if ($i == $NLONS[$i0]) { $i0++; $i1+=$NLONS[$i0] };
            print OUT "$zlat+$i0, ",($i-$i1)*360/$NLONS[$i0],", $x\n";
          }
        } else { print OUT "$x\n" }
      };
    } else {
# Binary output. But we allow 2 formats - 4-byte and 2-byte.
      $n=$#data+1;
# Maybe with header - in which case fake fortran-unformatted records to
# make it a pp-file.
      if ($HEADER) {
        print OUT pack("i",64*4),$Header,pack("i2",64*4,$n);
      } 
# Output data values. Either 2 or 4 byte
      if ($OUTFORM =~ /bi(n)?2/) {
# 2-byte format. Output number of pts, scale factor & ref value as 4-bytes,
        print OUT pack("if2",$totalpoints,$scale_factor,$ref_value);
# Then the data itself unchanged
        print OUT substr($in,8);
      } else {
# 4-byte format (default) - just the values.
        print OUT pack("f$n",@data)
      }; 
# Fake closing bit of f-unformatted too (only if HEADER set).
      if ($HEADER) { print OUT pack("i",$n) };
    };
    close OUT;
  };
};

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

}

print OUT1 "Read $N fields from file $Filei and wrote $NW\n";
print " ($NW/$N)\n";

# Keep track of total number written
$NWT+=$NW;

};

print "Total number written: $NWT\n";

# Write average, if requested. Txt or binary.
# Little check to make sure at least some have been written
# Output to last-written-file-name plus ".avg"
if ($AVG and $NWT) {
  foreach $x (@data_avg) {$x/=$NWT};
  open(OUT,"> $Fileo.avg");
  if ($OUTFORM eq "txt") { 
    foreach $x (@data_avg) { print OUT "$x\n" }
  } else {
    $n=$#data+1;
    print OUT pack("f$n",@data_avg);
  }
};

# Subroutine to unpack the reference values (and ref-value like
# objects) from 4 bytes

sub get_ref_val {

local ($in)=@_;

$ref_frac=vec($in,1,8)*256*256+vec($in,2,8)*256+vec($in,3,8);
$ref_mant=vec($in,0,8);
if ($ref_mant & 128) {$ref_sign=-1} else {$ref_sign=1};
$ref_mant=($ref_mant & 127) - 64;
$ref_sign*$ref_frac*2**(-24)*16**($ref_mant);

}