#!/bin/perl -w

# wunconv - re-insert modified pp-fields into unified model ancillary fields
#
# use:	  - wunconv.pl ancil-file pp-file output-ancil-file
#
# diags to: stdout
#
# method:
#	1. read in all of ancillary file (first argument)
#	2. read in pp-field(s) (second arg)
#	3. diagnostic info to wunconv.info
#	4. interpret header (mostly for info + to make sure we know what we're doing)
#	5. write pp-headers into lookup tables and data into data block 
#          (check that the headers are the same; warn if not)
#	6. write out result to stdout
#
# reference:
#	UMDP f3 (format of files), R. S. Bell, v 22 27/1/95

$N_ins=0;
$MAX_warn=5;
$Verbose=1;
$FIXHD10=-1; %FIXHD10=qw(0 single  1 timeseries  2 repeating);
$DATE=-1;
$ENDDATE=-1;
$MakeAncilOcean=0;				# Set H[1] to 2

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

if (scalar @ARGV != 3) {
  die "use: wunconv.pl ancil-file pp-file output-ancil-file"
};

# 1. readin ancil
#
undef $/;
$ANCIL=shift;
open(IN,$ANCIL) or die "Failed to open ancillary file <$ANCIL> (first arg)\n"; $IN=<IN>; close(IN);
$NL=length($IN);

# 2. readin pp(s)
#
$PPF=shift;
open(IN,$PPF) or die "Failed to open pp file <$PPF> (second arg)\n"; $PIN=<IN>; close(IN);
$NLP=length($PIN);

# 3. Output
#
$OUTFILE=shift;
if (!$OUTFILE) { die "No output file specified" };
open(OUT,"> $OUTFILE") or die "Can't open output to $OUTFILE";

# 4. Info
#
@H=unpack("i256",substr($IN,0,256*4));
if ($Verbose) {

  if ($H[1] == 1) { $m="Atmos" } 
  elsif ($H[1] == 2) { $m="Ocean" }
  else { $m="Unknown" };

  print "Ancil model type: $m\n";

  print "\n\n";

};

if ($Verbose > 1) { for ($I=0; $I<256; $I++) { 
  if ($H[$I] == -32768) {
    print sprintf("%4d%7s",$I+1,"null")
  } else {
    print sprintf("%4d%7d",$I+1,$H[$I])
  }; 
  if (($I+1) % 8 == 0) { print "\n" }; 
} };
if (($H[0] == -32768) or ($H[0] == 15)) { 
  if ($Verbose > 3) { print "H0=15 or null; OK\n" }
} else { 
  wwarn("H0 <> 15: $H[0]\n") 
};
if ($MakeAncilOcean) { $H[1]=2 };
@INTS=unpack("i*",substr($IN,$H[99]*4-4,$H[100]*4));
$NX=$INTS[5]; $NY=$INTS[6];		# Ancil specific
$NF=$INTS[14]*$H[151];			# Ancil specific

if ($Verbose > 0) {

  print "Ty: ",$H[1],"\n";
  print "Vn: ",$H[11]/100.,"\n";
  print "t0 (start) Y/M/D:H/M/S:D $H[20]/$H[21]/$H[22]:$H[23]/$H[24]/$H[25]:$H[26]\n";
  print "t1 (end)   Y/M/D:H/M/S:D $H[27]/$H[28]/$H[29]:$H[30]/$H[31]/$H[32]:$H[33]\n";
  print "t2 (skip)  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";
  print "Start, len of ints: $H[99],$H[100]\n";
  print "We seem to have ${NX}x${NY}x$NF\n";
  print "Start (StartI+LenI), len of real: $H[104] (",$H[99]+$H[100],"),$H[105]\n";
  print "Start (StartR+lenR), len of Lookup: $H[149] (",$H[104]+$H[105],"), \[$H[150],$H[151]\]\n";

};

if ($Verbose > 1) { print "  ints: ",join(",",@INTS),"\n" };

# Check out various constants...
if ($Verbose > 1) {
  
  print "Constants: start d1 (d2)\n";
  print "  Integer: $H[99] $H[100]\n";
  print "  Real   : $H[104] $H[105]\n";
  print "  LevDepC: $H[109] $H[110] ($H[111])\n";
  print "  RowDepC: $H[114] $H[115] ($H[116])\n";
  print "  ColDepC: $H[119] $H[120] ($H[121])\n";
  print "  FieldsC: $H[124] $H[125] ($H[126])\n";
  print "  Extra  : $H[129] $H[130]\n\n";
 
};

if ($NF == $H[151]) { 
  if ($Verbose > 2) { print "N Field types matches 2nd dim of lookup; good ($NF)\n" }
} else { 
  wwarn("N Field types ($NF) does not match 2nd dim of lookup ($H[151])\n");
};

# Start of data is officially item 160, with dimension 161
# But there may be a connection to start_of_lookup+length_of_lookup (150+151*152)
if ($Verbose > 0) { print "Start (StartLookup+lenL), len of data: $H[159] (",$H[149]+$H[150]*$H[151],"),$H[160]\n" };
if ($H[149]+$H[151]*$H[150] < $H[159]) { 
  if ($Verbose > 2) { print "Data starts ($H[159]) after end of lookup ($H[149]+$H[151]*$H[150]=",$H[149]+$H[151]*$H[150],") ... good!\n" } 
} else {
  wwarn("Data starts before end of lookup... bad! (probably fatal)\n")
};

if ($H[160] == $NX*$NY*$NF) { 
  if ($Verbose > 2) { print "Length of data matches NX*NY*NF; good\n" } 
} else { 
  warn "Length of data ($H[160]) != NX*NY*NF (",$NX*$NY*$NF,")\n" 
};

# Grid info (PP)
$pbdx=unpack("f",substr($PIN,(46+16)*4,4));
$pbzx=unpack("f",substr($PIN,(46+15)*4,4));
$pbdy=unpack("f",substr($PIN,(46+14)*4,4));
$pbzy=unpack("f",substr($PIN,(46+13)*4,4));
if ($Verbose > 1) { print "PP bdx,y: [$pbdx,$pbdy] bzx,y: [$pbzx,$pbzy]\n" };
# Grid info (ancil)



# Check expected pp-length
$LREC=$NX*$NY;
$PPL=4*(4+64+$LREC);
$LREC1=unpack("i",substr($PIN,15*4,4));
$NY1=unpack("i",substr($PIN,18*4,4));
$NX1=unpack("i",substr($PIN,19*4,4));
$PPL1=($LREC1 +64 +4)*4;
if ($PPL != $PPL1) { 
  print "  these PP-fields are not the same as the ancil. Ah well. $PPL. $PPL1\n" 
} else {
  if ($Verbose > 3) { print "  these PP-fields are the same as the ancil: good: $PPL\n" }
};
$PPL=$PPL1; $LREC=$LREC1;
$NX=$NX1; $NY=$NY1;
$NPP=$NLP/$PPL;
if ($Verbose > 1) { print "Expected pp-field length is: $PPL; that would mean we read in ",$NPP," fields\n" };

# 5. Insert fields
#
if ($N_ins) { $N=$N_ins } else { $N=$H[151] };
if ($H[151] != $NPP or $H[151] != $N_ins) { 

  if ($Verbose > 1) { print " H151,2: $H[150],$H[151] NPP: $NPP N_ins: $N_ins \n" };
  wwarn("I think I read in $NPP fields but I'm inserting $N and the original ancil had $H[151]\n");
  if ($Verbose > 0) { print "I'm going to change H[152(151)] (currently $H[151]) to $N (or try to...)\n" };
  substr($IN,151*4,4)=pack "i",$N;
  if ($Verbose > 0) { print "I think I'll change NX ($NX1) and NY ($NY1) too... (to ${NX}x${NY})\n" };
  substr($IN,$H[99]*4-4+5*4,4)=pack "i",$NX;
  substr($IN,$H[99]*4-4+6*4,4)=pack "i",$NY;
# Start-of-data seems to pad from end-of-lookup so its at 16385 (and 2^14=16384, by chance).
# So rather than just putting our stuff immeadiately after lookup, lets pad too.
  $OurDataStart=(int(($H[149]+$N*$H[150])/16384)+1)*16384+1;
  if ($Verbose > 0) { print "I think I'll change 160 (Start of data: $H[159]) and 161 (Dim of data: $H[160]) too... (to: ".($OurDataStart).", ".($N*$LREC).")\n" };
  substr($IN,159*4,4)=pack "i",$OurDataStart;
  substr($IN,160*4,4)=pack "i",$N*$LREC;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));
};
if ($FIXHD10 != -1) {
  print "Changing Fixhd(10 (ftn indexing)) from $H[9] ($FIXHD10{$H[9]}) to $FIXHD10 ($FIXHD10{$FIXHD10})\n";
  substr($IN,9*4,4)=pack "i",$FIXHD10;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));
};
if ($DATE != -1) {
  print "changing Fixhd(21 (ftn indexing)) from $H[20] to $DATE \n";
  substr($IN,20*4,4)=pack "i",$DATE;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));  
};
if ($ENDDATE != -1) {
  print "change Fixhd(28 (ftn indexing)) from $H[27] to $ENDDATE \n";
  substr($IN,27*4,4)=pack "i",$ENDDATE;
# Redo this with the new values.
  @H=unpack("i256",substr($IN,0,256*4));  
};  
for ($F=0; $F<$N; $F++) {
  if ($F >= $NPP) {
    print "Run out of pp-fields to insert; stopping\n"; last
  };
  if ($Verbose > 1 ) { print "Inserting field number $F / Header from: ",$H[149]+$F*$H[150]," to ",$H[149]+(1+$F)*$H[150]-1, " / Data from  : ",$H[159]+$F*$LREC," to ",$H[159]+(1+$F)*$LREC-1,"\n" };
  $PP=substr($PIN,$PPL*$F,$PPL);
  $PPH=substr($PP,4,64*4); $PPD=substr($PP,4*(3+64),4*$LREC);
  if ($PPH ne substr($IN,($H[149]+$F*$H[150]-1)*4,$H[150]*4)) {
    wwarn("header inserting for field $F does not match the one in the ancil\n")
  } else {
    if ($Verbose > 2) { print "  header to insert matches that in ancil; good\n" };
  };
# Insert pp header into lookup
# Remembering to change LBEGIN
  substr($IN,($H[149]+$F*$H[150]-1)*4,$H[150]*4)=$PPH;
  substr($IN,($H[149]+$F*$H[150]-1+29-1)*4,4)=pack("i",($H[159]+$F*$LREC-1));
# Insert data into data
  substr($IN,($H[159]+$F*$LREC-1)*4,$LREC*4)=$PPD;
  @data=unpack "f*",$PPD;
  if ($Verbose > 3) { print "data range: ".makerange(\@data,unpack("f",substr($PPH,(46+16)*4,4))).", [0,0]: $data[0]\n" };

};
if ($Verbose > 1) { print "(Length (words) of file read in: $NL (",$NL/4,"))\n" };

# 6. Write out result
print OUT $IN;



sub wwarn {

  $MAX_warn--;
  if ($MAX_warn >= 0) {
    print " *** Warning: @_"
  };
  if ($MAX_warn == 0) { print " ... and this is the last warning! (set MAX_warn=nn on the command line if you want more)\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)."] (mdi: $bmdi)";

}