#!/bin/perl -w

#
# wconv - take UKMO unified model ***4 byte*** ancillary fields and convert them to pp-fields
#
# Use: wconv.pl [options] infile [outfile]
#
# If not set, outfile defaults to infile.pp
#
# method:
#	1. read in all of file (first argument or stdin)
#	2. diagnostic info to wconv.info
#	3. interpret header (mostly for info + to make sure we know what we're doing)
#	4. read pp-headers from lookup tables and data from data block and write to stdout
#
# options:
#	arbitrary_skip=nnn - set padding between data records to value nnn
#	arbitrary_skip=auto - set padding so that records start on 2048-boundaries
#	arbitrary_skip=lbegin - set padding from LBEGIN from LOOKUP table
#
#	Verbose=n - set verbose output mode. Can be useful if things are going wrong...
#
# Notes:
#	The format of the actual data is something of a mess. Usuable ancils sometimes
#	have padding between data, sometimes not. Hence the arbitrary_skip stuff.
#
# reference:
#	UMDP f3 (format of files), R. S. Bell, v 22 27/1/95

# 0. Setup
$Verbose=1;
$DefOut=0;
$arbitrary_skip=0;
# Command-line setup
eval "\$$1=\$2" while $ARGV[0] =~ /^(\w+)=(.*)/ && shift;
if (scalar @ARGV != 2) { 

  if (scalar @ARGV == 1) {

    $DefOut=1;

  } else {

    die "use: wconv.pl [options] infile [outfile]" 

  };

};

# 1. readin
#
$INFILE=shift;
if ($DefOut) { $OUTFILE="$INFILE.pp" } else { $OUTFILE=shift };
undef $/;
open IN,$INFILE or die "Can't open $INFILE"; $IN=<IN>; close IN;
$NL=length($IN);

# 2. Diags to stdout, output to OUTFILE
#
if (!$OUTFILE) { die "No output specified" };
open OUT,"> $OUTFILE" or die "Can't open $OUTFILE";

# 3. Info
#
@H=unpack("i256",substr($IN,0,256*4));
if ($Verbose > 1) { for ($I=0; $I<256; $I++) { print sprintf("%4d%7d",$I+1,$H[$I]); if (($I+1) % 8 == 0) { print "\n" }; }; };
if ($H[0] == 15 or $H[0] == -32768) { print "H0=$H[0] (= 15 or -32768); OK\n" } else { warn "H0 <> 15 or -32768: is: $H[0]\n" };
print "Vn: ",$H[11]/100.,"\n";
print "(t0) Y/M/D:H/M/S:D $H[20]/$H[21]/$H[22]:$H[23]/$H[24]/$H[25]:$H[26]\n";
print "(t1) Y/M/D:H/M/S:D $H[27]/$H[28]/$H[29]:$H[30]/$H[31]/$H[32]:$H[33]\n";
print "(t2) Y/M/D:H/M/S:D $H[34]/$H[35]/$H[36]:$H[37]/$H[38]/$H[39]:$H[40]\n";
print "Global (0): $H[3], Type (4=ancil): $H[4]\n";
@INTS=unpack("i*",substr($IN,$H[99]*4-4,$H[100]*4));
if ($Verbose > 1) { print "  ints: ",join(",",@INTS),"\n" };
$NX=$INTS[5]; $NY=$INTS[6];		# Ancil specific
$NF=$INTS[14]*$H[151];				# Ancil specific
# Probably should be *H[151]... see wunconv.pl
print "Data array: ${NX}x${NY}x$NF\n";
print "Start, len of integer consts: $H[99], $H[100]\n";
print "Start (StartI+LenI), len of real consts: $H[104] (",$H[99]+$H[100],"),$H[105]\n";
print "Start (StartR+lenR), len, [dim1,dim2] of lookup table: $H[149] (",$H[104]+$H[105],"), \[$H[150],$H[151]\]\n";
print "Start (StartL+lenL), len of data: $H[159] (",$H[149]+$H[150]*$H[151],"),$H[160]\n";
if ($H[160] == $NX*$NY*$NF) { 
  print "Length of data matches NX*NY*NF; good\n" 
} else { 
  warn "Length of data ($H[160]) != NX*NY*NF (",$NX*$NY*$NF,")\n" 
};

# 4. Output fields
#
$Tot=0;
for ($F=0; $F<$H[151]; $F++) {
  print "F$F";
  print " Head: ",$H[149]+$F*$H[150],"-",$H[149]+(1+$F)*$H[150]-1;
  print " Data: ",$H[159]+$F*$NX*$NY,"-",$H[159]+(1+$F)*$NX*$NY-1;
  $h0=($H[149]+$F*$H[150]-1)*4;
  @lbuser=unpack("i7",substr($IN,$h0+38*4,7*4));
  print " Lbuser(1-7): ", join ",",@lbuser," (".$F*$NX*$NY.")\n";
  print OUT pack "i", $H[150]*4 or die "Failed to write Header-1"; $Tot+=4;
  print OUT substr($IN,($H[149]+$F*$H[150]-1)*4,$H[150]*4) or die "Failed to write Header"; $Tot+=$H[150]*4;
  print OUT pack "i", $H[150]*4 or die "Failed to write Header-2"; $Tot+=4;
  print OUT pack "i", $NX*$NY*4 or die "Failed to write data header"; $Tot+=4;

# Unpack some info from LOOKUP table 
# We don't need year or mon - thats just to check we are where we think we are.
# We use LBEGIN if arbitrary_skip is set to LBEGIN
  $LBYR=unpack("i",substr($IN,$h0+(1-1)*4,4));
  $LBMON=unpack("i",substr($IN,$h0+(2-1)*4,4));
  $LBEGIN=unpack("i",substr($IN,$h0+(29-1)*4,4));
#  print "\n\nLBYR,LBMON LBEGIN, *4: $LBYR,$LBMON $LBEGIN, ",$LBEGIN*4,"\n\n";

# From tlc's IDL: length_padding=(2048-(size_data mod 2048)) mod 2048
# This is only true sometimes!
  $lpad=(2048-($NX*$NY % 2048)) % 2048;
  if ($arbitrary_skip eq "auto") { 
    $lpad1=$lpad ;
# Or, we might decide to trust LBEGIN from the LOOKUP table
  } elsif ($arbitrary_skip eq "lbegin") { 
    $lpad1=($LBEGIN-($H[159]+$F*($NX*$NY)-1)); if ($F > 0) { $lpad1=$lpad1/$F };
  } else { 
    $lpad1=$arbitrary_skip 
  };
  $p=($H[159]+$F*($NX*$NY+$lpad1)-1)*4;
  $l=$NX*$NY*4;
# lbuser(0) tells us if data is real (1) integer (2) or logical (3; treat as 2)
  if ($lbuser[0] == 2 or $lbuser[0] == 3) {
    print "Note: data type is $lbuser[0]\n";
    @data=unpack "i*",substr($IN,$p,$l)
  } else {
    @data=unpack "f*",substr($IN,$p,$l)
  }
  print "data range: ".makerange(\@data,unpack("f",substr($IN,$h0+(46+16)*4,4))).", [0,0]: $data[0] (Lpad, arbitrary skip, data_from: $lpad1, $arbitrary_skip, $p)\n";
  print OUT substr($IN,$p,$l) or die "Failed to write data"; $Tot+=$NX*$NY*4;
  if (length($IN) < $p+$l-1) {
    print "\n***Warning! Not enough data available...\n***You asked for $p to ".($p+$l)." and I only had up to ".length($IN)."\n\n";
  };
  print OUT pack "i", $NX*$NY*4 or die "Failed to write tail"; $Tot+=4;
};
print "(Length (words) of file read in: $NL (",$NL/4,")). Written: $Tot of $H[151] fields\n";


sub makerange {

  ($rdata,$bmdi)=@_;

  $begun=0;
  for (@$rdata) {
    if (abs($_-$bmdi) > abs($bmdi)/100) {
      if (!$begun) { $begun=1; $mn=$_; $mx=$_ };
      if ($_>$mx) { $mx=$_ };
      if ($_<$mn) { $mn=$_ }
    }
  };
  return "[".sprintf("%.2f",$mn).",".sprintf("%.2f",$mx)."]";

}