#! /usr/bin/perl -w

package PrecisPPPackage;

use strict;
use Config;
use Exporter;
use lib "$ENV{UPATH}/perl"; # This is where modules PPPackage and PrecisGeneral live
use Carp;
use PrecisGeneral qw(abs_path
		     acos
		     asin
		     balanced_parentheses
		     decimal
		     file_handle
		     hash_to_f90
		     integer
		     make_x_to_360
		     mkdir_p
		     %month_to_integer
		     natural
		     read_namelist
		     read_STASHmaster
		     real
                     reg2rot
		     replace_directory
		     rm
		     rot2reg
		     same_file
		     shift_180
		     shift_360
		     %Sm
		     time_difference
		     translate_time
		     );

my ($word_size, $marker_size);
my ($machine_bigendian, $machine_int, $bs_machine_int, $machine_float);

BEGIN {
  # Can not assign variables to constants, since "use" pragmas get run
  # at compile time. So need to define these variables in a BEGIN
  # block so that they have a value when the 'use constants' are
  # processed.

  # Set the word size to 4 bytes. This is assuming that even on 64-bit
  # operating systems we are going to be dealing with 32-bit files (Note that
  # the word size of the operating system can be found with: $word_size=length
  # $Config{byteorder};)
  $word_size = 4;

  # The record marker is more tricky - that's right out of Perl's
  # jurisdiction. This is dodgy, but I can not think of another option.
  #
  # It is assumed that 32-bit PP files have 32-bit block control words and
  # anything else has no block control words. This is not a very good
  # assumption!
  $marker_size = ($word_size == 4) ? $word_size : 0;

  $machine_bigendian = 1 if ($Config{byteorder} =~ /4321$/);

  if ($word_size == 4) {
      $machine_float  = 'f';
      $machine_int    = 'i';
      $bs_machine_int = $machine_bigendian ? 'V' : 'N';
  }
}

our $VERSION   = 1.7.1;

our @ISA       = qw(Exporter);
our @EXPORT    = ();
our @EXPORT_OK = qw(@pp_keys
		    @pp_keys_header		   
		    %pp_keys
		    change_ppfile_suffix
		    check_ppfiles
		    create_global_pp
		    create_lam_containing_pp
		    decode_exp
		    encode_exp
		    find_pp4lbc	    
		    global
		    initialize_pp 
		    input_file_names
		    output_file_names
		    p_to_uv
		    parse_pp_header_relationships
		    pp_grid_info
		    pp_gridbox_weights	
		    pp_header_match
		    pp_inside_ppref 
		    pp_set
		    pp_time1
		    pp_time2
		    ppdata_2darray
		    pp_latitude
		    pp_longitude
		    rcm_gridbox_edge_extrema_in_reg_coords
		    read_pp
		    read_pp_header 
		    read_pp_data
		    short_to_long_pp_header_items
		    skip_pp_data
		    subregion_ppref
		    write_pp 
		    wrong_pp_endian
		    );

# Set constants
use constant WORD_SIZE   => $word_size;         # No. of bytes in a word
use constant MARKER_SIZE => $marker_size;       # No. of bytes in a record marker (can be 0)
use constant WORD        => 'a' x WORD_SIZE;    # Unpack for number of bytes in a word 
use constant NLHD        => 45;                 # No. of words in integer header
use constant NBHD        => 19;                 # No. of words in real    header
use constant HEADER_SIZE => 64 * WORD_SIZE;     # No. of bytes in the PP header
use constant MACHINE_BE  => $machine_bigendian; # Defined if machine is big-endian
use constant INTEGER     => $machine_int;       # Unpack for machine integer
use constant BS_INTEGER  => $bs_machine_int;    # Unpack for byte-swapped machine int
use constant FLOAT       => $machine_float;     # Unpack for machine float

my $pi              = atan2(1,1) * 4;
my $degrees2radians = $pi/180.0;

$ENV{VN} = $ENV{UMVERSION} if (!defined $ENV{VN} or $ENV{VN} eq '');
confess "ERROR: Neeed to set enviroment variable VN or UMVERSION\n"
    if (!defined $ENV{VN} or $ENV{VN} eq '');

# Define all of the PP hash keys in order
our @pp_keys_header = qw(lbyr      
			 lbmon     
			 lbdat     
			 lbhr      
			 lbmin     
			 lbday     
			 lbyrd     
			 lbmond    
			 lbdatd    
			 lbhrd     
			 lbmind    
			 lbdayd     
			 lbtim     
			 lbft      
			 lblrec    
			 lbcode    
			 lbhem     
			 lbrow     
			 lbnpt     
			 lbext     
			 lbpack    
			 lbrel     
			 lbfc      
			 lbcfc     
			 lbproc    
			 lbvc      
			 lbrvc     
			 lbexp     
			 lbegin    
			 lbnrec    
			 lbproj    
			 lbtyp     
			 lblev     
			 lbrsvd1    
			 lbrsvd2    
			 lbrsvd3   
			 lbrsvd4    
			 lbsrce    
			 lbuser1 
			 lbuser2 
			 lbuser3 
			 lbuser4 
			 lbuser5 
			 lbuser6 
			 lbuser7 
			 brsvd1     
			 brsvd2     
			 brsvd3     
			 brsvd4     
			 bdatum    
			 bacc      
			 blev      
			 brlev     
			 bhlev     
			 bhrlev    
			 bplat     
			 bplon     
			 bgor      
			 bzy       
			 bdy       
			 bzx       
			 bdx       
			 bmdi      
			 bmks
			 );

our @pp_keys =  ('bcw1',@pp_keys_header,'bcw2','bcw3','data','bcw4');

# Define a hash of the PP keys. The 64 words of the PP header are
# given numbers 1 to 64 in order.
our %pp_keys;
my $n = 0;
foreach my $key (@pp_keys) {
    $pp_keys{$key} = $n++;
}

sub initialize_pp {
    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    my %pp;    
    for my $n (0..$#_) {
	$pp{$pp_keys_header[$n]} = $_[$n];
    }
    $pp{data} = [ (0.0) ];
    return %pp;
}

sub pp_set {
    #---------------------------------------------------------------------------
    #---------------------------------------------------------------------------
    my %pp;    
    
    if (@_) { 
	confess "ERROR: Can only initialize a PP hash from a 64 word array\n" 
	    if (@_ != 64);
	for my $n (0..63) {
	    $pp{$pp_keys_header[$n]} = $_[$n];
	}
    } else {
	for my $key (@pp_keys_header) {
	    $pp{$key} = 0;
	}
	$pp{lbrel} = 2;    
    } 

    foreach my $key (qw(bcw1 bcw2)) {
	$pp{$key} = 256;
    }
    foreach my $key (qw(bcw3 bcw4)) {
	$pp{$key} = 0;
    }
    
    $pp{data}  = [ (0.0) ];
    
    return %pp;
}

sub read_pp {
    #---------------------------------------------------------------------------
    # Read an entire PP field (header and data and, if present, block control
    # words and extra data) into a hash.
    #
    # Call: %pp = read_pp(\*FILEHANDLE);
    #     : %pp = read_pp($filename);
    #---------------------------------------------------------------------------
    my $arg = shift; # A filename or open filehandle

    my $FH = file_handle($arg, 'read', 'binary');

    # Read PP header
    my %pp = read_pp_header(\*$FH);

    # Return nothing if the PP header not defined
    return () if (!%pp);

    # Read the PP data
    %pp = read_pp_data(\*$FH, \%pp);

    # Return the PP hash	    
    return %pp;
}

sub read_pp_header {
    #---------------------------------------------------------------------------
    # Read a PP header into a hash. Reads integer and real header and, if
    # present, the first three block control words and if present any extra data
    #
    # Call: %pp = read_pp_header(\*FILEHANDLE);
    #     : %pp = read_pp_header($filename);
    #---------------------------------------------------------------------------
    my $arg = shift; # A filename or open filehandle

    my $FH = file_handle($arg, 'read', 'binary');

    my %pp;
    my $packed;

    # Read first block control word (bcw1), if it exists
    my $bcw1;        
    if (MARKER_SIZE > 0) {
	my $record_length = MARKER_SIZE; 
	my $bytes_io = sysread $FH, $packed, $record_length;
	return () unless (defined $bytes_io && $bytes_io == $record_length);
	$bcw1 = unpack INTEGER, $packed;
    }

    # Read integer and real header (and bcw2, bcw3 if they exist)
    my $record_length = HEADER_SIZE + 2 * MARKER_SIZE;
    my $bytes_io = sysread $FH, $packed, $record_length;
    return () unless (defined $bytes_io && $bytes_io == $record_length);

    # Unpack header (and bcw[23] if they exist)
    my @header;

    if (defined $bcw1) {
	if (wrong_pp_endian($bcw1)) {

	    #-----------------------------------------------------------------
	    # PP field has BCWs and is NOT machine endian => Byte-swap numbers
	    #-----------------------------------------------------------------

	    # Unpack integer header into native integers (did it this
	    # way to avoid problems with -ve integers)
	    my $offset = 0;
	    foreach (1..NLHD) {
		push @header, unpack(INTEGER,
				     pack(WORD,
					  reverse
					  unpack WORD,
					  substr($packed, $offset, WORD_SIZE)
					  )
				     );
		$offset += WORD_SIZE;
	    }

	    # Unpack real header into native reals (thanks to Emma
	    # Hibling for this bit)
	    foreach (1..NBHD) {
		push @header, unpack(FLOAT,
				     pack(WORD,
					  reverse
					  unpack WORD,
					  substr($packed, $offset, WORD_SIZE)
					  )
				     );
		$offset += WORD_SIZE;
	    }	    

	    # Unpack bcw2 into a non-machine endian integer (bcw1 is
	    # non-machine endian, too. Note that if this field ever gets
	    # written out with subroutine write_pp then bcw1 and bcw2 will be
	    # converted to machine endian)
	    push @header, unpack INTEGER, substr($packed, $offset, MARKER_SIZE);

	    # Unpack bcw3 into a native integer
	    $offset += MARKER_SIZE;
	    push @header, unpack(INTEGER,
				 pack(WORD,
				      reverse
				      unpack WORD,
				      substr($packed, $offset, MARKER_SIZE)
				      )
				 );
	} else {	

	    #----------------------------------------
	    # PP field has BCWs and is machine endian
	    #----------------------------------------

	    # Unpack lhd, bhd, bcw2, bcw3 into native numbers
	    @header = unpack INTEGER.NLHD.FLOAT.NBHD.INTEGER.'2', $packed;
	}	

	# Enter block control words into PP hash
	$pp{bcw1} = $bcw1;
	$pp{bcw3} = pop @header;
	$pp{bcw2} = pop @header;

    } else {

	#-------------------------------------------
	# PP field has no BCWs and is machine endian
	#-------------------------------------------

	# Unpack lhd, bhd into native numbers
	@header = unpack INTEGER.NLHD.FLOAT.NBHD, $packed;
    }

    # Enter PP header into PP hash    
    foreach (0..$#pp_keys_header) {
	$pp{$pp_keys_header[$_]} = $header[$_];
    }

    # Return PP hash (with no data key)
    return %pp;
}

sub read_pp_data {
    #---------------------------------------------------------------------------
    # Read the data array of a PP field from into a PP hash
    #
    # Call: %pp = read_pp_data(\*FILEHANDLE, \%pp);
    #---------------------------------------------------------------------------
    my $FH = shift; # An open filehandle
    my $pp = shift; # HASH reference: PP hash
   
    #----------------------------------------------------
    # Read data as a packed string, and bcw4 if it exists
    #----------------------------------------------------
    my $record_length = $pp->{lblrec} * WORD_SIZE + MARKER_SIZE;
    my $packed;
    my @data;

    my $bytes_io = sysread $FH, $packed, $record_length;
    return () unless (defined $bytes_io && $bytes_io == $record_length);
    
    # Unpack data 
    if (MARKER_SIZE == 0) {

	#-------------------------------------------
	# PP field has no BCWs and is machine endian
	#-------------------------------------------

        # Unpack data into into native reals
	@data = unpack FLOAT.$pp->{lblrec}, $packed;	

    } else {
	if (wrong_pp_endian($pp->{bcw1})) {

	    #--------------------------------------------
	    # PP field has BCWs and is NOT machine endian
	    #--------------------------------------------

	    # Unpack data into native reals (thanks to Emma Hibling for this bit)
	    my $offset = 0;
	    foreach (0..$pp->{lblrec}-1) {
		push @data, unpack FLOAT, pack(WORD,
					       reverse
					       unpack WORD,
					       substr($packed, $offset, WORD_SIZE)
					       );
		$offset += WORD_SIZE;
	    }

	    # Unpack bcw4 into a native integer
	    $pp->{bcw4} = unpack(INTEGER,
				 pack(WORD,
				      reverse
				      unpack WORD,
				      substr($packed, $offset, MARKER_SIZE)
				      )
				 );

	} else {

	    #----------------------------------------
	    # PP field has BCWs and is machine endian
	    #----------------------------------------

	    # Unpack data and bcw4 into native numbers
	    @data = unpack FLOAT.$pp->{lblrec}.INTEGER, $packed;
	    $pp->{bcw4} = pop @data;

	    # Throw away extra data, for now .....
	    if ($pp->{lbext} > 0) {
		$pp->{lblrec} -= $pp->{lbext};
		$pp->{lbext}   = 0;
		@_ = @data[ 0..$pp->{lblrec}-1 ];
		@data = @_;
		print "size data=",(scalar @data),"\n";
	    }
	}
    }

    # Enter data into PP hash
    $pp->{data} = [ @data ];  

    # Change logical data to 0s and 1s
    if ($pp->{lbuser1} == 3) {
	foreach (@{ $pp->{data} }) {
	    $_ = 1.0 if ($_ != 0);
	}
    }    

    # Return PP hash	    
    return %$pp;
}

sub write_pp {
    #---------------------------------------------------------------------------
    # Write a PP hash to a PP field on disk
    #
    # Call: $bytes = write_pp(\*FILEHANDLE, \%pp);
    #     : $bytes = write_pp($filename, \%pp);
    #---------------------------------------------------------------------------
    my $arg = shift; # A filename or filehandle
    my $pp  = shift; # HASH reference: PP hash

    my $FH = file_handle($arg, 'write', 'binary');

    # Check that all of the keys have been defined, but no need to check BCWs
    foreach my $key ((@pp_keys_header,'data')) {
	confess "ERROR: PP hash key \'$key\' is not defined" if (!exists $pp->{$key});
    }

    my $template;
    my @pack_list;
    my $field_size = HEADER_SIZE + $pp->{lblrec} * WORD_SIZE;  # Bytes in PP field (no BCWs)

    if (MARKER_SIZE > 0) {
	#------------------------
	# Use block control words
	#------------------------
	# Set BCWs around data
	$pp->{bcw1} = 256;
	$pp->{bcw2} = $pp->{bcw1};
	$pp->{bcw3} = $pp->{lblrec} * WORD_SIZE;
	$pp->{bcw4} = $pp->{bcw3};

	# Create a pack template for the whole PP field
	$template = INTEGER.INTEGER.NLHD.FLOAT.NBHD.INTEGER.'2'.FLOAT.$pp->{lblrec}.INTEGER;

	# Create list to pack
	push @pack_list, $pp->{bcw1};
	foreach my $key (@pp_keys_header) {
	    push @pack_list, $pp->{$key};
	}
	push @pack_list, ($pp->{bcw2}, $pp->{bcw3}, @{ $pp->{data} }, $pp->{bcw4});
	$field_size += 4 * MARKER_SIZE;      # Add four BCWs to size of field
    } else {
        #-----------------------
	# No block control words
        #-----------------------
	# Create a pack template for the whole PP field
	$template = INTEGER.NLHD.FLOAT.NBHD.FLOAT.$pp->{lblrec};

	# Create list to pack
	foreach my $key (@pp_keys_header) {
	    push @pack_list, $pp->{$key};
	}
	push @pack_list, @{ $pp->{data} };
  }

    # Pack the list and write it to disk
    my $bytes_io = syswrite $FH, (pack $template, @pack_list);
    confess "ERROR: Problem writing PP field.\n".
	    "       Only $bytes_io bytes of an expected $field_size written\n" 
	if ($bytes_io != $field_size);    
    confess "ERROR: Problem writing PP field\n" if (!defined $bytes_io);  

    return $bytes_io;
}

sub rcm_extrema_in_reg_coords {
# Return:
#
# west  = most westerly , in reg lat/lon coords, of RCM grid box centres
# south = most southerly, in reg lat/lon coords, of RCM grid box centres
# east  = most easterly , in reg lat/lon coords, of RCM grid box centres
# north = most northerly, in reg lat/lon coords, of RCM grid box centres
#
# Requires a PP header hash as input (see sub find_pp_header)
    my %pp = @_;  # RCM PP header

    my ($west,$south,$east,$north);
    
    my @lons = $pp{bzx} + $pp{bdx}; # rotated grid longitudes
    my @lats = $pp{bzy} + $pp{bdy}; # rotated grid latitudes
    foreach (2..$pp{lbnpt}) {
	push @lons, $lons[-1] + $pp{bdx};
    }
    foreach (2..$pp{lbrow}) { 
	push @lats, $lats[-1] + $pp{bdy}; 
    }

    # Convert to radians
    foreach (@lons) {
	# Convert longitudes are between -pi and pi
	$_ -= 360.0  if abs $_ > 180.0;
	$_ *= $degrees2radians;
    }
    foreach (@lats) {
	$_ *= $degrees2radians;
    }
    $pp{bplat} *= $degrees2radians;
    $pp{bplon} *= $degrees2radians;
    
    # Rotated grid extrema of rotated coords
    my ($northlat,$westlon,$southlat,$eastlon) = ($lats[0],$lons[0],$lats[-1],$lons[-1]);   
    my ($sock,$cpart,$t1,$t2,$ires,$x); # Variables for rotation calculation
    my (@rlatt,@rlong);                 # Variables for rotation calculation
    
    
# Define SOCK
    if ($pp{bplon} != 0) {
	$sock = $pp{bplon} - $pi;
    } else {
	$sock = 0;
    }
    
#-----------------------
# Find $north and $south
#-----------------------
    foreach my $lat ($northlat,$southlat) {
	foreach (0..@lons-1) {
	    $cpart = (cos $lons[$_]) * cos $lat;	    
	    $x = ((cos $pp{bplat})*$cpart) + (sin $pp{bplat}) * sin $lat;
	    $x =  1.0 if $x >=  1.0;
	    $x = -1.0 if $x <= -1.0;	    
	    $rlatt[$_] = &asin($x);	    
	    # Convert back to degrees
	    $rlatt[$_] = $rlatt[$_] / $degrees2radians;	    
	    # Find the extrema
	    if ($lat == $northlat && $_ == 0) {
		# First time throught loops
		$north = $rlatt[0];
		$south = $rlatt[0];
	    } else {
		# Subsequent time throught loops
		$north = $rlatt[$_] if $rlatt[$_] > $north;
		$south = $rlatt[$_] if $rlatt[$_] < $south;
	    }
	}
    }

#---------------------
# Find $east and $west
#---------------------
    foreach my $lon ($westlon,$eastlon) {
	foreach (0..@lats-1) {
	    $cpart = (cos $lon) * cos $lats[$_];	    
	    $x = ((cos $pp{bplat}) * $cpart)+(sin $pp{bplat}) * sin $lats[$_];
	    $x =  1.0 if $x >=  1.0;
	    $x = -1.0 if $x <= -1.0;	  
	    $rlatt[$_] = &asin($x);	    
	    $t1 = (-1 * cos $pp{bplat}) * sin $lats[$_];
	    $t2 = (sin $pp{bplat}) * $cpart;	    
	    $x = ($t1 + $t2)/cos $rlatt[$_];
	    $x =  1.0 if $x >=  1.0;
	    $x = -1.0 if $x <= -1.0;	    
	    $rlong[$_] = -1 * &acos($x);	    
	    $rlong[$_]=-1*$rlong[$_] if ($lon >= 0.0 && $lon <= $pi);
	    $rlong[$_] = $sock + $rlong[$_];	    
	    # Convert back to degrees
	    $rlong[$_] = $rlong[$_] / $degrees2radians;	    
	    # Find the extrema
	    if ($lon == $westlon && $_ == 0) {
		# First time throught loops
		$east = $rlong[0];
		$west = $rlong[0];
	    } else {
		# Subsequent time throught loops
		$east = $rlong[$_] if $rlong[$_] > $east;
		$west = $rlong[$_] if $rlong[$_] < $west;
	    }
	}
    }

    return ($west,$south,$east,$north);
}

{ 
    # Wrap this subroutine and use caching to make it faster.
    my %decode_cache;
    sub decode_exp {
	my $coded_exp = $_[0];
	if (exists $decode_cache{$coded_exp}) {
	    return $decode_cache{$coded_exp};
	}

	# TODO AFAIK, not portable to T3E!
	my $runid;
	if ($coded_exp < 0) {
	    $runid = '-----';
	    return $runid;
	}

	my @char = qw(a b c d e f g h i j k l m n o p q r
		      s t u v w x y z 0 1 2 3 4 5 6 7 8 9);
	
	# First convert to bit string.
	my $bits = unpack 'B32', (pack 'N', $coded_exp);
	
	# Step through 6 bits at a time (missing off first 2),
	# converting each 6 bit chunk into a decimal. The decimal
	# is the index into the character lookup array. 
	my ($i, $chunk, $j);
	for ($i=2; $i < 32; $i+=6) {
	    $chunk = substr $bits, $i, 6;
	    $j = unpack 'N', pack('B32', substr('0' x 32 . $chunk, -32));
	    $runid .= $char[$j];
	}
	
	$decode_cache{$coded_exp} = $runid;
	return $runid;
    }  # end of wrapped subroutine
}

{ 

    # Wrap this subroutine and use caching to make it faster. I don't
    # understand what this means, but it's what Emma does.
    my %encode_cache;
    sub encode_exp {
	my $runid = shift;

	if (exists $encode_cache{$runid}) {
	    return $encode_cache{$runid};
	}

	if ($runid eq '-----') {
	    return -32768;
	}
                                                                           
	my @char = qw(a b c d e f g h i j k l m n o p q r
		      s t u v w x y z 0 1 2 3 4 5 6 7 8 9);
	my $nchar = @char;

	my $coded_exp = 0;
	for my $i (0..length($runid)-1) {
	    my $letter = substr($runid, -1-$i, 1);
	    my $newnum;
	    for my $j (0..$#char) {
		if ($letter eq $char[$j]) {
		    $newnum = $j*(2**($i*6));                                        
		    $coded_exp += $newnum;
		    last;                                       
		}
	    }
	    confess "ERROR: Incorrect letter in runid: $letter\n" if (!defined $newnum);
	}

	return $coded_exp;

    }  # end of wrapped subroutine
}

sub pp_gridbox_weights {
    #---------------------------------------------------------------------------
    # Calculate the weights of each grid box (in km**2) of a PP field. Only
    # works for regular lat-lon grids, i.e. lbcode = 1,2,101,102
    #
    # Call: @weights = pp_gridbox_weights(\%pp);
    #---------------------------------------------------------------------------
    my $pp = shift; # HASH reference: PP hash

    #--------------------------------------------
    # Return empty lists in certain circumstances
    #--------------------------------------------
    # Return empty list if the number rows or columns is not set
    return () if ($pp->{lbrow} < 1 or $pp->{lbnpt} < 1);

    # Return empty list for LBC fields
    return () if ($pp->{lbhem} == 99);

    # Return empty list for LBC rim data
    return () if ($pp->{lblev} == 7777);

    # Return empty list if the field is packed
    return () if (($pp->{lbpack} % 10) > 0);

    # Return empty list for non regular lat-lon grids
    my $code = $pp->{lbcode};
    if (! ($code == 1 or $code == 2 or $code == 101 or $code == 102)) {
	return ();
    }

    # Return list of 1s for UK national grid data
    return (1.0) x $pp->{lblrec} if ($pp->{lbhem} == 4444);

    use constant PI          => 3.14159265358979323846;
    use constant PI_OVER_180 => PI / 180.0;
    use constant PI_OVER_2   => PI / 2.0;
    use constant RADIUS      => 6371229.0; # Radius of the earth in metres

    # Define grid constants from the field header
    my $zy  = $pp->{bzy} * PI_OVER_180; 
    my $dy  = $pp->{bdy} * PI_OVER_180; 
    my $dy2 = $dy / 2.0;
    my $dy4 = $dy / 4.0;
    my $dx  = abs($pp->{bdx});
  
    # Initialise source grid latitudes (grid box edges)
    my @ycoord;
    $ycoord[0] = $zy + $dy2;
    foreach my $j (1..$pp->{lbrow}) {
	push @ycoord, $ycoord[$j-1] + $dy;
    }

    # Find area (divided by 2*pi*r^2) of each dome defined by the north
    # and south edges of each row.
    # Note that if rows 0 and/or $pp->{lbrow}-1 are polar rows then $dome[0]
    # and/or dome[$pp->{lbrow}] will contain rubbish
    my @dome;
    foreach my $j (0..$pp->{lbrow}) {
	$dome[$j] = 1.0 - sin($ycoord[$j]);
    }
    
    # Find the area (divided by 2*pi*r^2) of each global latitudinal
    # gridbox band
    my @band; 
    if ($dy < 0) {
	foreach my $j (0..$pp->{lbrow}-1) {
	    $band[$j] = $dome[$j+1] - $dome[$j];
	}
    } else {
	foreach my $j (0..$pp->{lbrow}-1) {
	    $band[$j] = $dome[$j] - $dome[$j+1];
	}
    }   
    
    # Correctly set box area (divided by dxrad_r2) for polar rows. 
    # Note that 2 = 1-sin(-pi/2) = surface area of earth (divided by
    # 2*pi*r^2)
    if ($dy < 0) {
	if ($ycoord[0] > PI_OVER_2 - $dy4) {
	    # Weights in north polar row have been divided by $pp->{lbnpt}
	    $band[0] = $dome[1] / $pp->{lbnpt};
	}
	if ($ycoord[$pp->{lbrow}] < -1*PI_OVER_2 + $dy4) {
	    # Weights in south polar row have been divided by $pp->{lbnpt}
	    $band[$pp->{lbrow}-1] = (2.0 - $dome[$pp->{lbrow}-1]) / $pp->{lbnpt};	 
	}
    } else {
	if ($ycoord[0] < -1*PI_OVER_2 - $dy4) {
            # Weights in south polar row have been divided by $pp->{lbnpt}
	    $band[0] = (2.0 - $dome[1]) / $pp->{lbnpt};
	}
	if ($ycoord[$pp->{lbrow}] > PI_OVER_2 + $dy4) {
	    # Weights in north polar row have been divided by $pp->{lbnpt}
	    $band[$pp->{lbrow}-1] = $dome[$pp->{lbrow}-1] / $pp->{lbnpt};
	}
    }
    
    # Find the area (divided by 2*pi*r^2) of a ficticious global
    # latitudinal gridbox band centred on the equator
    my $equator_band = abs(sin(-$dy2) - sin($dy2));
    
    # Find the relative area-weighting for each gridbox in the domain
    # (all gridboxes on a row will have the same weighting)
    my @weights;
    foreach my $j (0..$pp->{lbrow}-1) {
	push @weights, ($band[$j]/$equator_band) x $pp->{lbnpt};
    }
	
    return @weights;
}

sub p_to_uv {

    # Functionally the same as UM vn4.5 deck PTOUV1A.dk
    
    my %pp = @_;   # Incoming P-grid PP field
    my %pp_u = (); # returned UV-grid PP field
    
    if ($pp{lbhem} == 0 or $pp{lbhem} == 3 or $pp{lbhem} == 4) {
	%pp_u = %pp;
	$pp_u{lbrow}--;
	$pp_u{lblrec} = $pp_u{lbnpt} * $pp_u{lbrow};
	$pp_u{bzx} += $pp_u{bdx}/2.0;
	$pp_u{bzy} += $pp_u{bdy}/2.0;
	$pp_u{data}   = [ ($pp_u{bmdi}) * $pp_u{lblrec} ];
	my $i;

	foreach $i (0..$pp_u{lblrec}-2) {
	    $pp_u{data}[$i] = 0.25 * 
		($pp{data}[$i] + 
		 $pp{data}[$i+1] + 
		 $pp{data}[$i+$pp{lbnpt}] + 
		 $pp{data}[$i+1+$pp{lbnpt}]);
	}
	
	# Set end points
	if ($pp{lbhem} == 3) {
	    # No cyclic wrap around: Set last values on row
	    for ($i = $pp{lbnpt}-1; $i < $pp_u{lblrec} ; $i += $pp{lbnpt}) {
		$pp_u{data}[$i]=$pp_u{data}[$i-1];
	    }
	} else {
	    # Cyclic wrap around                                                       	
	    for ($i = $pp{lbnpt}-1; $i < $pp_u{lblrec} ; $i += $pp{lbnpt}) {
		$pp_u{data}[$i] = 0.25 *
		    ($pp{data}[$i] + 
		     $pp{data}[$i+1-$pp{lbnpt}] + 
		     $pp{data}[$i+$pp{lbnpt}] + 
		     $pp{data}[$i+1]);
	    }
	} 
    } # if ($pp{lbhem} == 0 or $pp{lbhem} == 3 or $pp{lbhem} == 4)

    return %pp_u;
}

sub wrong_pp_endian {
    #---------------------------------------------------------------------------
    # If both the endian of both the machine and the PP file can be determined
    # and they are the same, then returns undef, otherwise returns 1.
    #
    # Assumes that if a file has BCWs then it is a 32-bit file.
    #
    # Call: $ok = wrong_pp_endian($bcw1);
    #---------------------------------------------------------------------------
    my $bcw1 = shift;

    if (defined $bcw1) {
	if ($bcw1 == 256) {
	    # PP field is machine endian
	    return;
	} elsif ($bcw1 == 65536) {
	    # PP field is not machine endian
	    return 1; 
	} else {
	    confess "ERROR: Can not determine byte order of PP field: bcw1 = $bcw1\n";
	} 
    }

    # Can not determine byte order of PP field, so assume that it is
    # machine endian
    return;
}

sub ppdata_2darray {
    my $pp  = shift;
    my $tlc = shift;

    my @data;
    my $rs = 0; # Index of start of row
    for (my $row = 0; $row < $pp->{lbrow}; $row++) {
	$data[$row] = [ @{ $pp->{data} }[($rs..$rs+$pp->{lbnpt}-1)] ];
	$rs += $pp->{lbnpt};
    }	
    
    if (defined $tlc) {

	# Reverse columns? 
	if ($pp->{bdx} < 0.0) {
	    my $n = 0;
	    foreach my $row (@data) { 
		my @row = @$row;
		@row = reverse @row;
		$data[$n] = [ @row ];
		$n++;
	    }
	    # Adjust bdx and bzx
	    $pp->{bzx} += $pp->{bdx}*($pp->{lbnpt}+1);
	    $pp->{bdx} *= -1;
	}
	
	# Reverse rows?    
	if ($pp->{bdy} > 0.0) {
	    @data = reverse @data ;
	    # Adjust bdy and bzy
	    $pp->{bzy} += $pp->{bdy}*($pp->{lbrow}+1);
	    $pp->{bdy} *= -1;
	}
    }

    $pp->{data} = [ @data ];

    return %$pp;
}

sub input_file_names {
    #-----------------------------------------------------------------------
    # Call  $n    = input_file_names(\@files ,$file_or_FH
    #                                [,$omit_regex] [,$incl_regex] [,$all]);
    #       @list = input_file_names(\@files, undef
    #                                [,$omit_regex] [,$incl_regex] [,$all]);
    #-----------------------------------------------------------------------
  
    my $input_files = shift; # Reference to an array of files and/or directories
    my $list_in     = shift;
    my $omit_regex  = shift; # Regular expression
    my $incl_regex  = shift; # Regular expression
    my $all         = shift;

    my $n_in = 0; # OUT if  defined $list_in
    my @files;    # OUT if !defined $list_in

    my $FH_IN = file_handle($list_in, 'write') if (defined $list_in);

    $omit_regex = qr{$omit_regex} if (defined $omit_regex);
    $incl_regex = qr{$incl_regex} if (defined $incl_regex);
    
    foreach my $file (@$input_files) {
	if (-d $file) {
	    # Argument is a directory. Finds all non-empty, binary files in tree
	    if (defined $list_in) {
		$n_in += files_in_all_subdirs_($file,\*$FH_IN,$omit_regex,$incl_regex,$all);
	    } else {
		push @files, files_in_all_subdirs_($file,undef,$omit_regex,$incl_regex,$all);
	    }
	    
	} elsif (! -s $file or ! -B _) {

	    die "ERROR: Nonexistent, empty or non-binary input file: ".
		abs_path($file)."\n";

	} else {

	    # Argument is a binary file
	    $file = abs_path($file);
	    # Omit any files from the list
	    next if (defined $omit_regex and $file =~ /$omit_regex/);
	    next if (defined $incl_regex and $file !~ /$incl_regex/);
	    if (defined $list_in) {
		# Write file name to disk
		print $FH_IN "$file\n"; 
		$n_in++;
	    } else {
		# Add file name to array of file names
		push @files, $file;
	    }
	}
    }
    
    if (defined $list_in) {
	close $FH_IN;
	if (-s $list_in) {
	    return $n_in; 
	} else {
	    die "ERROR: No input files found\n";
	}
    } else {
	if (@files) {
	    if (wantarray()) {
		return @files;
	    } elsif (defined wantarray()) {
		return $files[0];
	    } else {
		return;
	    }
	} else {
	    die "ERROR: No input files found\n";
	} 
    }
    
    sub files_in_all_subdirs_ {
	# Recursively called subroutine to find all non-empty, binary files in
	# all subdirectories
	my $file       = shift;
	my $FH_IN      = shift; # Open file handle
	my $omit_regex = shift;
	my $incl_regex = shift;
	my $all        = shift;
	
	my $n = 0; # OUT if $FH_IN is defined: Number of valid files in tree

	if (defined $FH_IN) {
	    foreach my $file (glob("$file/*")) {
		if (-d $file) {
		    $n += files_in_all_subdirs_($file,\*$FH_IN,$omit_regex,$incl_regex,$all);
		} elsif (-s $file and -B _ and ($all or $file =~ /\.pp$/)) {
		    $file = abs_path($file);		    
                    # Omit any files from the list whose names match a regular
                    # expression
		    next if (defined $omit_regex and $file =~ /$omit_regex/);
		    next if (defined $incl_regex and $file !~ /$incl_regex/);
		    # Write file name to disk
		    print $FH_IN "$file\n";		    
		    $n++;
		}
	    }
            # Return the number of files
	    return $n;

	} else {

	    my @files;
	    
	    foreach my $file (glob("$file/*")) {
		if (-d $file) {
		    push @files, files_in_all_subdirs_($file,undef,$omit_regex,$incl_regex,$all);
		} elsif (-s $file and -B _ and ($all or $file =~ /\.pp$/)) {
		    $file = abs_path($file);		    
		    # Omit any files from the list whose names match a regular
		    # expression
		    next if (defined $omit_regex and $file =~ /$omit_regex/);		    
		    next if (defined $incl_regex and $file !~ /$incl_regex/);
		    # Add file name to array of file names
		    push @files, $file;		    
		}
                # Return the list of file names
		return @files;
	    }
	}
    }
} # sub input_file_names

sub output_file_names {

    #---------------------------------------------------------------------------
    # Call 1: output_file_names(\$out,$delete,$list_in);
    #	    
    #         * If $out is defined:
    #           - Replaces $out in the calling routine with its full path
    #           - Checks that $out is different to all files listed in $list_in
    #           - Removes $out
    #           - Returns 1
    #	    
    #         * If $out is undefined:
    #           - Dies
    #	    
    # Call 2: output_file_names(\$out,$delete,$list_in,$list_out);
    #	    
    #         * If $out is defined:
    #           - Replaces $out in the calling routine with its full path
    #           - Checks that $out different to all files listed in $list_in
    #           - Writes $out to new file $list_out
    #           - Removes $out
    #           - Returns 1
    #	    
    #         * If $out is undefined:
    #           - See Call 3
    #	    
    # Call 3: output_file_names(\$out,$delete,$list_in,$list_out,$suffix,\$directory);
    #	    
    #         * If $out is defined:           
    #           - $list_out is ignored
    #           - $suffix   is ignored
    #           - If $out has an absolute or relative path or $directory is undefined:
    #             - Replaces $out in the calling routine with its full path
    #           - If $directory is defined:
    #             - Replaces $directory in the calling routine with its full path
    #             - If $out has neither an absolute nor a relative path:
    #               - Changes the path of $out in the calling routine with to $directory
    #           - Checks that $out different to all files listed in $list_in 
    #           - Removes $out
    #           - Returns 1
    #	    
    #         * If $out is undefined:
    #           - Creates an output file name in file $list_out for every input file
    #             name in file $list_in
    #           - Modifies the suffix and/or path of the output file names if $suffix
    #             and/or $directory are defined 
    #           - If $directory is defined:
    #             - Replaces $directory in the calling routine with its full path
    #           - Removes all files listed in file $list_out
    #           - Returns the number of output file names in file $list_out
    #
    #---------------------------------------------------------------------------
    
    my $out       = shift; # SCALAR reference: File name. Modifies input. May be undef
    my $delete    = shift; #
    my $list_in   = shift; # Must be defined.
    my $list_out  = shift; # May be undef. 
    my $suffix    = shift; # May be undef. 
    my $directory = shift; # SCALAR reference: Directory name. Modifies input. May be undef.

    my $n_out = 0;  # Returned number of output files
    
    confess "ERROR: \$list_in is not defined.\n" if (!defined $list_in);
    open my $FH_IN, "< $list_in" or die "ERROR: Can not open $list_in: $!"; 
    
    $$directory = mkdir_p($$directory) if (defined $$directory);

    if (defined $$out) {

        # $out is defined
        # In this case, $suffix is ignored

	if (defined $$directory and $$out !~ /\//) {
	    $$out = replace_directory($$directory,$$out);
	} else {	    
	    $$out = abs_path($$out);
	}

	while (local $_ = <$FH_IN>) {
	    chomp;
	    die "ERROR: Identical input and output files: $_\n" if (same_file($$out,$_));
	}
	if (defined $list_out) {
	    open my $FH_OUT, "> $list_out" or die "ERROR: Can not open $list_out: $!"; 
	    print $FH_OUT $$out."\n";
	    close $FH_OUT;
	}

	# Output file is different to all input files. Delete if requested.
	rm($$out) if ($delete);

	$n_out = 1;

    } elsif (defined $list_out) {

	open my $FH_OUT, "> $list_out" or die "ERROR: Can not open $list_out: $!"; 
	while (local $_ = <$FH_IN>) {
	    chomp;
	    my $out = $_;
	    $out = change_ppfile_suffix($suffix, $_)    if (defined $suffix);
	    $out = replace_directory($$directory, $out) if (defined $$directory);
	    die "ERROR: Identical input and output files: $_\n" if (same_file($out,$_));
	    print $FH_OUT "$out\n";

	    # Output file is different to all input files. Delete if requested.
	    rm($out) if ($delete);

	    $n_out++;
	}
	close $FH_OUT;

    } else {
	confess "ERROR: Neither \$out nor \$list_out are defined.\n";
    }

    close $FH_IN;

    return $n_out;
} #  sub output_file_names

sub check_ppfiles {

    #----------------------------------------------------------------
    # Call: @ppfiles = check_ppfiles(@files);
    #     : $ppfile  = check_ppfiles($file);
    #----------------------------------------------------------------

    local @_ = @_;

    if (@_) {
	foreach (@_) {
	    $_ = abs_path($_);
	    if (-d $_) {
		die "ERROR: $_ is a directory\n";
	    } elsif (! -s _ or ! -B _) {
		die "ERROR: Nonexistent, empty or non-binary file: $_\n";
	    }
	}
    }

    # Return
    if (wantarray()) {
	return @_;
    } elsif (defined wantarray()) {
	return $_[0];
    } else {
	return;
    }
}

sub change_ppfile_suffix { 
    
    my $suffix = shift; 
    
    local @_ = @_;
    foreach (@_) {
	$_ .= $suffix if (! s/\.pp$/$suffix/);
    }

    if (wantarray()) {
        return @_;
    } elsif (defined wantarray()) {
        return $_[0];
    } else {
        return;
    }
}

sub parse_pp_header_relationships {
    #---------------------------------------------------------------------------
    #
    #
    # Call: %h = parse_pp_header_relationships($header, undef); 
    #     : %h = parse_pp_header_relationships($header, $f90, $ignore_lhs); 
    #---------------------------------------------------------------------------
    my $expression = shift; # String containing expression
    my $f90        = shift; # T = create F90 namelist for testing expression
    my $type       = shift; # 'ppwhere': Parse 'x>67 and not y<87.4' where the 
                            #            right hand operands are ignored
    
    return () if (!defined $expression);
    
    my %h;         # The returned hash
    
    $h{type} = (defined $type) ? $type : '';
    
    # Check parentheses
    balanced_parentheses($expression) or die "ERROR: $@\n";

    # Reformat expression
    $expression = reformat_expression_($expression, \%h);

    # Replace operands
    $expression = operands_($expression, \%h);

    # Reduce expression to irreducible partial expressions and create a nice
    # description
    reduce_expression_($expression, \%h);

    # Add a to the hash a string contining all the hash information an F90
    # namelist
    create_f90_namelist_(\%h) if ($f90);
    
    return %h;
    
    #---------------------------------------------------------------------------
    # End of subroutine. Subroutines start here.
    #---------------------------------------------------------------------------
    
    sub reformat_expression_ {

	local $_ = shift;
	my $h    = shift;
	
	my $len = length($_);
	
	$h->{info} .= sprintf "  %-${len}s : original expression\n", $_;
	
	# Spaces	
	s/\s+//g;
	s/(and|or)/ $1 /g; 
	s/not/not /g;
	$h->{info} .= sprintf "= %-${len}s : space\n", $_;
	
	# Change to lower case
	$h->{info} .= sprintf "= %-${len}s : lower case\n", $_
	    if (/[A-Z]/ and $_ = lc);

	# Change all brackets to round ones
	$h->{info} .= sprintf "= %-${len}s : [] -> ()\n",   $_ if (tr/[]/()/); 
	$h->{info} .= sprintf "= %-${len}s : {} -> ()\n",   $_ if (tr/{}/()/);
	
	# Define a month regular expression.
	# Matches a 1 or two digit integer
	# Matches a three letter month abbreviation, case insentsitive
	my $month_regex = qr{ jan(?:uary)?       |
				  feb(?:ruary)?  |
				  mar(?:ch)?     |
				  apr(?:il)?     |
				  may            |
				  jun(?:e)?      |
				  jul(?:y)?      |
				  aug(?:ust)?    |
				  sep(?:tember)? | 
				  oct(?:ober)?   |
				  nov(?:ember)?  |
				  dec(?:ember)?
			      }ix;

	# Define a season regular expression.
	my $ssn_regex = qr{ djf | mam | jja | son }ix;

	#----------------------------
	# Check the expression syntax
	#----------------------------
	{

	    # Very important to localize $_ here, as it's about to get cleared
	    # but it is needed later
	    local $_ = $_;	    	   

            # Define a time regular expression.	    
	    my $time_regex = qr{ (?:\d+[,-\._]){0,5}\d+[,-\._]? | # Delimiter separated
				 [a-z\d]{5}                     # UM datestamp
			       }x;

	    # General check for disallowed characters
	    die "ERROR: Incorrect character in expression: $1\n" 
		if ( /([^a-z\d,\._><=!\/\s\-\+\(\)])/);

	    # Strip out valid elements of expression, dying if there is
	    # anything left at the end

            # Define a real number regular expression
	    my $real = qr{ [-\+]? (?: \d+\.?\d* | \.\d+) }x; 
	    
	    while (length) {

		# Find number of characters in expression
		my $len = length($_);

		# 'not' and open brackets
		while (s/\A \(* not \s \(*//x) {}; 

		# Yet more open brackets
		s/\A \(* //x; 

		# Header item (fairly loose check)
		die "ERROR: Incorrect expression syntax. Bad item: $_\n"
		    if (! s/\A ([a-z]{1,6}\d?) (?= [><=!\/])//x);
		
		my $item = $1;
		
		# Relational operator
		die "ERROR: Incorrect expression syntax. ".
		    "Bad relational operator: $_\n"
		    if (! s/\A ( >=?|<=?|==?|!=|\/= ) (?= [a-z\.\-\+\d\(])//x);

		my $operator = $1;

		# Value of item
	      ITEM: {
		  if ($item =~ /^(?:vtime|dtime|period)$/) {
		      die "ERROR: Incorrect expression syntax. ".
			  "Bad value for item \'$item\': $item$operator$_\n"
			  if (! s/\A $time_regex (?! [a-z\d]) //x);
		      last ITEM;
		  }
		  if ($item =~ /^(?:lbexp|exp)$/) {
		      die "ERROR: Incorrect expression syntax. ".
			  "Bad value for item \'$item\': $item$operator$_\n"
			  if (! (s/\A $real      (?! [a-z]  ) //x or 
				 s/\A [a-z\d]{5} (?! [a-z\d]) //x));
		      last ITEM;
		  }
		  if ($item =~ /^(lb)?mon$/) {
		      die "ERROR: Incorrect expression syntax. ".
			  "Bad value for item \'$item\': $item$operator$_\n"
			  if (! (s/\A $real        //x or 
			         s/\A $month_regex //x or 
			         s/\A $ssn_regex   //x));
		      last ITEM;
		  }
 		  if ($item =~ /^(lb)?mond?$/) {
		      die "ERROR: Incorrect expression syntax. ".
			  "Bad value for item \'$item\': $item$operator$_\n"
			  if (! (s/\A $real        //x or 
			         s/\A $month_regex //x));
		      last ITEM;
		  }
                  # Still here?		  
		  die "ERROR: Incorrect expression syntax. ".
		      "Bad value for PP header item \'$item\': $item$operator$_\n"
		      if (! s/\A $real (?! [a-z])//x);
	      }
		
		# Close brackets
		s/\A\)*//; 

		# ' and ' or ' or ' followed by yet more open brackets
		s/\A \s (and|or) \s \(* //x; 

		# Die if the nothing has been stripped out of the expression
		die "ERROR: Incorrect expression syntax (6): $_\n" 
		    if ($len == length);
	    }
	}	    

	# Change times to time integers
        # E.g. 'vtime=1960,12,3,4,17,0' to 'vtime=19601203041700'
        # E.g. 'vtime=g0c3h'            to 'vtime=19601203041700'
        #-------------------------------------------------------
        # C. Wang (10/0/2011) made following changes:
        # Change expression to ((\d{1,4}[,\._]){0,5}\d{1,4}[,\._]? (?! [a-z\d]))
        # ((\d{1,4}[,-\._]*?){0,5}\d{1,4}[,-\._]? (?! [a-z\d])) so that
        # vtime etc can be given 8 digits or with /,- etc. e.g
        # vtime=19600101 or vtime=1960-01-01
 
	s/
	    (vtime|dtime|period)
	    ([><=!\/]{1,2})
    	    ((\d{1,4}[,-\._]*?){0,5}\d{1,4}[,-\._]? (?! [a-z\d]))
	 /
	    { $1.$2.translate_time(\$3, 'integer') }
	 /egx;
	
	# Change lettered months to integers
        # E.g. 'mon=apr' to 'mon=4`
	s/
	    ((?:lb)?mond?)
	    ([><=!\/]{1,2})
	    ($month_regex)
	 /
	    { $1.$2.$month_to_integer{lc($3)} }
	 /egx;

	# Change lettered seasons to month values
        # E.g. 'mon=jja' to '(mon=6 or mon=7 or mon=8)`
	if (s/
	    ((?:lb)?mon?)
	     ([><=!\/]{1,2})
	    ($ssn_regex)
      	     /
  	    { season_to_months__(lc($3)) }
	     /egx) {
	
	    my $op = $2;
	
    	    die "ERROR: Incorrect expression syntax. ".
	        "Bad relational operator for seasons: $2\n"
	        if ($op !~ /^==?$/);
        }

	return $_;

        sub season_to_months__  {
	    my $ssn = shift;
	  CASE: { 
	      if ($ssn  eq 'djf') { $_ = "(mon=12 or mon=1 or mon=2)" ; last CASE };  
	      if ($ssn  eq 'mam') { $_ = "(mon=3 or mon=4 or mon=5)"  ; last CASE };
	      if ($ssn  eq 'jja') { $_ = "(mon=6 or mon=7 or mon=8)"  ; last CASE };
	      if ($ssn  eq 'son') { $_ = "(mon=9 or mon=10 or mon=11)"; last CASE };
	    }
	    return $_;
        }

    } # sub reformat_expression_
    
    
    sub operands_ {
	#-----------------------------------------------------------------------
	# Replace operands to the 'and' and 'or' operators with a variable
	# name (p0, p1, ...) . Each operand is an arithmetic expression
	# relating a PP header item to a number. E.g. lbtim=13, dx>0.1

        # Also replace comma separated dates with integers
	#-----------------------------------------------------------------------
	local $_ = shift;
	my $h    = shift;

	my $len = length($_);

	# $real matches 5, .5, 0.5, -5, -.5, -0.5, +5, +.5, +0.5
	my $real = qr{ [-\+]? (?: \d+\.?\d* | \.\d+) }x; 

	# $regexp matches 'bdx>-0.67' or 'exp==addfa' or 'vtime<=1960,1,1'
	my $regexp = qr{      ([a-z]{1,6}\d?) \s*
			      ([><=!\/]{1,2}) \s* 
			      ($real                              (?! [a-z]  ) |
			       [a-z\d]{5}                         (?! [a-z\d]) | 
			       \d{14}                             (?! [a-z\d]) | 
			       (\d{1,4}[,-\._]){0,5}\d{1,4}[,-\._]? (?! [a-z\d]))
			  }x; 

	# Add to expression hash    
	my $p = 0;
	while (s/$regexp/p$p/) {

	    my $value;

	    push @{ $h->{item} }, $1;
	    push @{ $h->{rel}  }, $2;
	    push @{ $h->{c}    }, $3;

	    if ($h->{type} eq 'ppwhere') {
		die "ERROR: Need a better error message here!\n" 
		    if ($h->{c}[-1] =~ /[a-z]/);
	    } else {
		
		$h->{item}[-1] = short_to_long_pp_header_items(1, $h->{item}[-1]) ;
		
		if ($h->{item}[-1] =~ /^(vtime|dtime|period)$/) {
		    $value      = translate_time(\$h->{c}[-1], 'comma');
		    $h->{c}[-1] = decimal(translate_time(\$h->{c}[-1], 'integer'));
		} elsif ($h->{item}[-1] eq 'lbexp' and $h->{c}[-1] =~ /[a-z]/) {
		    $h->{c}[-1] = encode_exp($h->{c}[-1]);
		}
		
		if ($h->{item}[-1] =~ /^lb/) {
		    die "ERROR: Incorrect expression syntax. ".
   "Non-integer value of integer item: $h->{item}[-1]$h->{rel}[-1]$h->{c}[-1]\n"
			if ($h->{c}[-1] !~ /^\d+$/);
		}
    
	    }
	    
	    $h->{c}[-1] = decimal($h->{c}[-1]);

	    $value = $h->{c}[-1] if (!defined $value);
	    
	    $h->{rel}[-1] =~ s/\/=/!=/;           # Change /= to !=
	    $h->{rel}[-1] =~ s/^=$/==/;           # Change =  to ==

	    $h->{info} .= sprintf "= %-${len}s : %s\n", 
	                  $_, "p$p = $h->{item}[-1] $h->{rel}[-1] $value";

	    $p++;
	}

	# Number of operands
	$h->{n_var} = $p;

	die "ERROR: Could not idenitfy any operands\n" if ($h->{n_var} == 0);

	return $_;
    }  # sub operands_

    sub reduce_expression_ {
	#-----------------------------------------------------------------------
	# Reduce expression to a sequence irreducible partial expressions
	# which, when evaluated by subroutine pp_header_match, will return a
	# true or false answer to the question "does this PP field satisfy the
	# PP header relationships described in the original expression?"
	#
        # Call: reduce_expression_($expression, \%h);
	#-----------------------------------------------------------------------
	local $_ = shift;
	my $h    = shift;

	my @temp;

	my $len = length($_);
	my $t = 0;

	# Reduce until we end up with something like 't5' or 'not t5'
	while (!/\A (?:not\s)? [a-z]\d+ \Z/x) { 

	    foreach my $op (qw(and or)) {

		my $regexp = qr{ ( (?:not\s)? [pt]\d+ \s $op \s (?:not\s)? [pt]\d+ ) }x;

		while (s/$regexp/t$t/) { 

		    $temp[$t] = $1;

		    $h->{info} .= sprintf "= %-${len}s : %s\n", $_, "t$t = $temp[$t]";
		    add_to_partial_expression_hash_($1, $h, $t);
		    $t++;

		    while (s/ \( ( (?:not\s)? [pt]\d+ ) \) /$1/gx) {
			$h->{info} .= sprintf "= %-${len}s : removed parentheses\n", $_;
		    }

		    # Only do one 'or' iteration at a time (because 'and' has
		    # higher precedence)
		    last if ($op eq 'or');
		}    
	    }
	}

	$h->{not} = 1 if (/\Anot/);

	$h->{n_temp} = $t;

	# Build description
	for (my $t=$h->{n_temp}-1; $t>=0; $t--) {
	    s/t$t/($temp[$t])/;
	}
	for (my $p = 0; $p < $h->{n_var}; $p++) {
	    my $value = $h->{c}[$p];
	    if ($h->{item}[$p] =~ /^(?:vtime|dtime|period)$/) {
		$value = substr("00000000000000$value",-14);	
		$value = translate_time(\$value, 'comma');
	    }

	    s/p$p/$h->{item}[$p]$h->{rel}[$p]$value/;	   
	}
	s/(\A\(|\)\Z)//g if (! /\Anot/); # Remove null outer brackets

	$h->{description} = $_;

    } # sub reduce_expression_ 

    sub add_to_partial_expression_hash_ {
	#-----------------------------------------------------------------------
	# 
	#
	# Call: add_to_partial_expression_hash_($parital_expression, \%h,  $t);
	#-----------------------------------------------------------------------
	local $_ = shift;
	my $h    = shift;
	my $t    = shift;

	# Initialize to undef;
	($h->{t1}[$t], 
	 $h->{t2}[$t],
	 $h->{p1}[$t], 
	 $h->{p2}[$t],
	 $h->{not1}[$t], 
	 $h->{not2}[$t]) = ();

	my @operand;

	$h->{not1}[$t] = 1 if (s/^not //);
	$h->{not2}[$t] = 1 if (s/not // );

	/^ ([pt]\d+) \s* (and|or) \s* ([pt]\d+) $/x;  # E.g. p3 and p5
	($operand[1], $h->{op}[$t], $operand[2]) = ($1, $2, $3);

	foreach my $n ((1,2)) {
	    $operand[$n] =~ / ([pt]) (\d+) /x;	    
	    $h->{$1.$n}[$t] = $2;                  # Set, e.g., $h->{p1} = 4
	}
    }    

    sub create_f90_namelist_ {
	#-----------------------------------------------------------------------
	# Call: create_f90_namelist_(\%h);
	#-----------------------------------------------------------------------
	my $h = shift;

	my %h2;

	foreach my $key (qw(t1 t2 p1 p2)) {	
	    foreach (@{ $h->{$key} }) {
		push @{ $h2{$key}      }, (defined $_) ? ($_+1) :  0 ; 
		push @{ $h2{'l_'.$key} }, (defined $_) ? 'T'    : 'F';
	    }
	}	    

	foreach my $key (qw(not1 not2)) {	
	    foreach (@{ $h->{$key} }) {
		push @{ $h2{$key} }, (defined $_) ? 'T' : 'F';
	    }
	}

	if ($h->{type} eq 'ppwhere') {
	    $h2{item} = [ (0) x @{ $h->{item} } ];
	} else {
	    $h2{item} = [ @{ $h->{item} } ];

	    my %pp_keys2 = %pp_keys;

	    # Set negative numbers for derived items
	    $pp_keys2{vtime}  = -1; # Validity/start time
	    $pp_keys2{dtime}  = -2; # Data/end time
	    $pp_keys2{period} = -3; # Period (End minus Start times)

	    foreach (@{ $h2{item} }) { 
		$_ = $pp_keys2{$_}
	    }
	}

	$h2{rel} = [ @{ $h->{rel} } ];
	foreach (@{ $h2{rel} }) {
	    s/!=/\/=/;
	    $_ = sprintf '\'%-2s\'', $_;
	}

	$h2{c} = [ @{ $h->{c} } ];

	if (exists $h->{op}) {
	    $h2{op} = [ @{ $h->{op} } ];
	    foreach (@{ $h2{op} }) { 
		$_ = sprintf '\'%-3s\'', $_;
	    }
	}

	foreach my $key (qw(n_temp)) {
	    $h2{$key} = $h->{$key};
	}

	$h2{test} = ($h->{n_var} > 0) ? 'T' : 'F';
	$h2{not}  = ($h->{not}      ) ? 'T' : 'F';

	my $namelist_name  = 'HEADERCRITERIA';
	my $structure_name = 'h';
	if ($h->{type} eq 'edit') {
	    $structure_name = 'e';
    	    $namelist_name  = 'HEADEREDITS';
    	}
	if ($h->{type} eq 'ppwhere') {
	    $structure_name = 'w';
	    $namelist_name  = 'EXPRESSIONDEFINITION';
	}

	$h->{f90_namelist} .= hash_to_f90(\%h2, $structure_name, $namelist_name);
    } # sub create_f90_namelist_ 

} # sub parse_pp_header_relationships

sub pp_header_match {
    #---------------------------------------------------------------------------
    #
    # Call: $match = pp_header_match(\%h, \%pp); 
    #---------------------------------------------------------------------------
    my $h  = shift; # HASH reference: Expression
    my $pp = shift; # HASH reference: PP

    #------------------------------------------
    # Always match if no criteria have been set
    #------------------------------------------
    return 1 if (! %$h);   

    #------------------------------------------------------------------------
    # Evaluate and return when only one relational expression (e.g. lbhr > 3)
    #------------------------------------------------------------------------
    if ($h->{n_temp} == 0) {
	return ((eval pp_header_value_($pp, $h->{item}[0])." $h->{rel}[0] $h->{c}[0]") == 1) ? 1 : undef;
    }

    #---------------------------------------------
    # Evaluate more than one relational expression
    # (e.g. lbhr > 3 and lbmon < 4)
    #---------------------------------------------
    my @temp;
    
    for (my $n=0; $n < $h->{n_temp}; $n++) {

	my $A;  # The result (true or false) of the LHS relational expression
	my $B;  # The result (true or false) of the RHS relational expression 
	
        # Determine if A is true or false
	if (defined $h->{t1}[$n]) {	    

	    # $A has already been determined in a previous pass through the
	    # loop
	    $A = $temp[$h->{t1}[$n]];
	    
	} elsif (defined $h->{p1}[$n]) {
	    
	    # Evaluate the LHS boolean expression
	    $A = ((eval pp_header_value_($pp, $h->{item}[$h->{p1}[$n]])."
                        $h->{rel}[$h->{p1}[$n]] 
                        $h->{c}[$h->{p1}[$n]]") == 1) ? 1 : 0;	    
	} 
	
	# Determine if B is true or false
	if (defined $h->{t2}[$n]) {
	    
	    # $B has already been determined in a previous pass through the
	    # loop
	    $B = $temp[$h->{t2}[$n]];
	    
	} elsif (defined $h->{p2}[$n]) {

	    # Evaluate the RHS boolean expression
	    $B = ((eval pp_header_value_($pp, $h->{item}[$h->{p2}[$n]])."
                        $h->{rel}[$h->{p2}[$n]] 
                        $h->{c}[$h->{p2}[$n]]") == 1) ? 1 : 0; 

	}

        # Negate A and/or B if required
	$A = (! $A) if ($h->{not1}[$n]);
	$B = (! $B) if ($h->{not2}[$n]);
	
	# Evaluate the boolean expression 'A operator B'
	$temp[$n] = ($A && $B and $h->{op}[$n] eq 'and' or 
		     $A || $B and $h->{op}[$n] eq 'or');
    }

    # Negate the final result if required
    $temp[-1] = (! $temp[-1]) if ($h->{not});
    
    # Return result
    return ($temp[-1]) ? 1 : undef;
    
    sub pp_header_value_ {
	#-----------------------------------------------------------------------
        # Call: $value = pp_header_value_(\%pp, $item);
	#-----------------------------------------------------------------------
	my $pp   = shift; # HASH reference: PP field
	my $item = shift; # A PP header item (e.g. lbuser4)
                          # or a derived item (e.g. vtime)

        # Item is a standard PP header item name, so just return its value
	return $pp->{$item} if (defined $pp->{$item});

	# Item is a derived from more than one standard PP header item
      ITEM: {
	  if ($item eq 'vtime') {
	      # Validity/start time
	      my %vtime = pp_vtime($pp);
	      return translate_time(\%vtime, 'integer') * 1;
	  }
	  if ($item eq 'dtime') {
	      # Data/end time
	      my %dtime = pp_dtime($pp);
	      return translate_time(\%dtime, 'integer') * 1;
	  }
	  if ($item eq 'period') {
	      # Period ('end time' minus 'start time')
	      my %period = pp_period($pp);
	      return translate_time(\%period, 'integer') * 1;
	  }
      } # ITEM:

    } # sub pp_header_value_

} # sub pp_header_match 

sub short_to_long_pp_header_items {

    #---------------------------------------------------------------------------
    # Expand PP header item strings without leading 'lb' or 'b' to their
    # full names, e.g. dx   -> bdx
    #                  npt  -> lbnpt
    #                  bgor -> bgor
    #                  lbyr -> lbyr
    #
    # Call: $item  = short_to_long_pp_header_items($allow_derived_items, $x);
    #     : @items = short_to_long_pp_header_items($allow_derived_items, @x);
    #---------------------------------------------------------------------------
    my $allow_derived_items = shift;

    return if (! @_);

    local @_ = @_;

  ITEM: foreach (@_) {
      
      # Skip over derived items
      next ITEM if ($allow_derived_items and /^(?:vtime|dtime|period)$/);	
      
      die "ERROR: Ambiguous PP header item: $_\n" if ($_ eq 'lev');
      
      $_ = 'lbuser4' if ($_ eq 'stash' or $_ eq 'st' or $_ eq 'sc');  	    
      
      if (! /^l?b/) {
	  if (exists $pp_keys{"b$_"}) {
	      $_ = "b$_";                 # Convert to real    header item
	  } elsif (exists $pp_keys{"lb$_"}) {
	      $_ = "lb$_";                # Convert to integer header item
	  } else {
	      die "ERROR: Incorrect PP header item: $_\n";
	  }
      } elsif (! exists $pp_keys{$_}) {
	  die "ERROR: Incorrect PP header item: $_\n";
      }
  } # ITEM:
    
    if (wantarray()) {
	return @_;
    } elsif (defined wantarray()) {
	return $_[0];
    } else {
	return;
    }
}

sub find_pp4lbc {
    #-------------------------------------------------------------------------
    # Call: @pp4lbcs = find_pp4lbc(source          => $ENV{source},
    # 				   umtime          => $current_time,
    # 				   aerosol         => $SULPHUR,
    # 				   pp4lbc_prefix1  => $pp4lbc_prefix1,
    # 				   pp4lbc_prefix2  => $pp4lbc_prefix2,
    # 				   JOBDIR          => $ENV{JOBDIR},
    # 				   PP4LBC_DIR      => $ENV{PP4LBC_DIR},
    # 				   LBC_FINDPP4LBC  => $LBC_FINDPP4LBC,
    # 				   LBC_CHECKPP4LBC => $LBC_CHECKPP4LBC,
    # 				   );
    #
    # The arguments are all optional and read into a hash and their so
    # order is irrelevant. Unset arguments default to undef, which might
    # cause problems.
    #
    #-------------------------------------------------------------------------

    my %a = @_; # Subroutine arguments

    my %rcm = pp_grid_info($a{JOBDIR});   # Get RCM grid details 

    my %filenames; # Names of returned standard and aerosol pp4lbc files
    my %streams;   # 


    # Create filenames and return if they are not to be found
    if (! $a{LBC_FINDPP4LBC}) {

	$filenames{standard}    = "$a{PP4LBC_DIR}/$a{pp4lbc_prefix1}$a{umtime}";	
	if ($a{aerosol}) {
	    $filenames{aerosol} = "$a{PP4LBC_DIR}/$a{pp4lbc_prefix2}$a{umtime}";
	}

	# Check that files which have already been named are ok
	foreach my $key (keys %filenames) {
	    my $pp4lbc = $filenames{$key};
	    die "ERROR: ".ucfirst($key)." PP4LBC file $pp4lbc nonexistent or empty.\n" 
		if (! -s $pp4lbc);
	    # Check that the PP4LBC spans the RCM domain
	    if ($a{LBC_CHECKPP4LBC}) {
		my %pp4lbc = read_pp2_header($pp4lbc);
		die "ERROR: ".ucfirst($key)." PP4LBC file $pp4lbc does not span RCM domain\n"
		    if (! pp_inside_ppref(\%rcm,\%pp4lbc));
	    }    
	}

	# Return now 
	return ($filenames{standard},$filenames{aerosol});
    }
 
    # Still here? Then find for some PP4LBC files

    $streams{standard} = 'pi';
    if ($a{aerosol}) {
	$streams{aerosol} = 'pg';
    }

    #opendir DIR, $a{PP4LBC_DIR} or 
#	die "ERROR: Can not opendir $a{PP4LBC_DIR}: $!";

    # Assume for now that LBCs are always for the atmosphere
    my $submodel = 'a';
    
    foreach my $stream (keys %streams) {
        
#        my @filelist = glob("$a{PP4LBC_DIR}/$a{source}$submodel.*$streams{standard}$a{umtime}*");
        my @filelist = glob("$a{PP4LBC_DIR}/$a{source}$submodel.*$streams{$stream}$a{umtime}*");
	
	die "ERROR: No $a{PP4LBC_DIR} $stream PP4LBC files for time $a{umtime} in $a{PP4LBC_DIR}\n" 
	    if (! @filelist);

	if ($a{LBC_CHECKPP4LBC}) {
	    #------------------------------------------------------
	    # Find the first PP4LBC file which spans the RCM domain
	    #------------------------------------------------------

	    foreach my $ppfile (@filelist) {
		
		next if (! -s $ppfile);
		
		# Get PP4LBC grid details
		my %pp4lbc = read_pp_header($ppfile);
		
		# Stop searching if RCM domain lies wholly inside this PP4LBC
		if (pp_inside_ppref(\%rcm,\%pp4lbc)) {
		    $filenames{$stream} = $ppfile;
		    last;
		}	    
	    }
	    if (!defined $filenames{$stream}) {
		die "ERROR: No $filenames{$stream} $stream PP4LBC files for time $a{umtime} in $a{PP4LBC_DIR}\n".
		    "       Files exist but none span RCM domain\n";
	    }

	} else {
	    
            #-------------------------------------
	    # Find the first non-empty PP4LBC file
	    #-------------------------------------

	    foreach my $ppfile (@filelist) {	       
		next if (! -s $ppfile);
		$filenames{$stream} = $ppfile;
		last;
	    }	
	    if (!defined $filenames{$stream}) {
		die "ERROR: No $filenames{$stream}  $stream PP4LBC files for time $a{umtime} in $a{PP4LBC_DIR}\n"
	    }
	}       		

    } # foreach $stream (keys %streams) 

    return ($filenames{standard},$filenames{aerosol});
}

sub rcm_gridbox_edge_extrema_in_reg_coords {
# Return:
#
# west  = most westerly , in reg lat/lon coords, of RCM grid box edges
# south = most southerly, in reg lat/lon coords, of RCM grid box edges
# east  = most easterly , in reg lat/lon coords, of RCM grid box edges
# north = most northerly, in reg lat/lon coords, of RCM grid box edges
#
# Requires a PP hash as input
#    print "UPATH=$ENV{UPATH}\n";
    my $pp = shift; # HASH reference: PP hash

    my ($reglon,$reglat);
    my @reglon;
    my @reglat;

    my @longitude = $pp->{bzx} + $pp->{bdx}/2.0;
    foreach my $i (1..$pp->{lbnpt}) {
        push @longitude, $longitude[0] + $i*$pp->{bdx};
    }

    my @latitude  = $pp->{bzy} + $pp->{bdy}/2.0;
    foreach my $i (1..$pp->{lbrow}) {
        push @latitude, $latitude[0] + $i*$pp->{bdy};
    }

    foreach my $lon (($longitude[0],$longitude[-1])) {
        foreach my $lat (@latitude) {
            ($reglon,$reglat) = rot2reg($lon,$lat,$pp->{bplon},$pp->{bplat});
            push @reglon, $reglon;
            push @reglat, $reglat;
        }
    }

    foreach my $lat (($latitude[0],$latitude[-1])) {
	foreach my $lon (@longitude) {
            ($reglon,$reglat) = rot2reg($lon,$lat,$pp->{bplon},$pp->{bplat});
            push @reglon, $reglon;
            push @reglat, $reglat;
        }
    }
    # At this moment, @reglon is in the  range [0,360)

    #----------------------------------------------------------------
    # @reglon values are in [0,360), but if the unrotated Greenwich
    # meridian passes through the domain, we nned to shift the values
    # to be in the range [-180,180)
    #----------------------------------------------------------------
    my $rcm_max_lat;
    my $rcm_min_lat;
    if ($pp->{bdy} > 0) {
	$rcm_max_lat = $latitude[-1];
	$rcm_min_lat = $latitude[0];
    } else {
	$rcm_max_lat = $latitude[0];
  	$rcm_min_lat = $latitude[-1];
    }

    my $rcm_max_lon;
    my $rcm_min_lon;
    if ($pp->{bdx} > 0) {
	$rcm_max_lon = $longitude[-1];
	$rcm_min_lon = $longitude[0];
    } else {	  
	$rcm_max_lon = $longitude[0];
  	$rcm_min_lon = $longitude[-1];
    }

    foreach my $reglat (@reglat) {

        # Find rotated location of unrotated greenwhich meridian at each
        # unrotated latitude in the RCM domain ($rot_g_lon,$rot_g_lat)
	my ($rot_g_lon,$rot_g_lat) = reg2rot(0,$reglat,$pp->{bplon},$pp->{bplat});

	# Adjust @reglon if ($rot_g_lon,$rot_g_lat) is inside the RCM domain
	$rot_g_lon = make_x_to_360($rcm_min_lon, $rot_g_lon);
	if ($rot_g_lon <= $rcm_max_lon and 
	    $rot_g_lat <= $rcm_max_lat and
	    $rot_g_lat >= $rcm_min_lat) {
	    @reglon = shift_180(@reglon);
	    last;
	}
    }
    
    # Sort arrays numerically
    @reglon = sort { $a <=> $b } @reglon;
    @reglat = sort { $a <=> $b } @reglat;

    my $west  = $reglon[0];
    my $east  = $reglon[-1];
    my $south = $reglat[0];
    my $north = $reglat[-1];

    $east -= 360.0 if ($east > $west + 360.0);

    return ($west,$south,$east,$north);
}

sub pp_inside_ppref {

    # Test if gridbox edges of PP lie inside gridbox centres of REF

    my $pp  = shift; # HASH reference: PP hash
    my $ref = shift; # HASH reference: PP hash

    # Find extrema of PP gridbox edges in unrotated coordinates
    my ($west,$south,$east,$north) = rcm_gridbox_edge_extrema_in_reg_coords(\%$pp);
    
    # Find extrema of reference PP gridbox centres
    my ($bottom,$top) =
	sort { $a <=> $b } ($ref->{bzy} + $ref->{bdy}, $ref->{bzy} + $ref->{bdy}*$ref->{lbrow});

    my ($left,$right) = 
	sort { $a <=> $b } ($ref->{bzx} + $ref->{bdx}, $ref->{bzx} + $ref->{bdx}*$ref->{lbnpt});

    $left=$left%360	;
    $west=$west%360	;
    while ($right < $left) { $right += 360 }        ;
    while ($right >=  $left+360) { $right -= 360 }  ;
    while ($east < $west) { $east += 360 }          ;
    while ($east >=  $west+360) { $east -= 360 }    ;
    
    my $pp_not_inside_ppref;

    # Is $north in ref domain?
    $pp_not_inside_ppref = 1 if ($north < $bottom or $north > $top);
    # Is $south in ref domain?
    $pp_not_inside_ppref = 1 if ($south < $bottom or $south > $top);

    # Only test east and west if reference PP field does not have 360
    # degree wrap around
    if ($right + abs($ref->{bdx}) < $left + 360) {
	# Is $east in reference domain?
	$pp_not_inside_ppref = 1 if ($east > $right);
	# Is $west in reference domain?
	$pp_not_inside_ppref = 1 if ($west < $left);
    }
    
    ($pp_not_inside_ppref) ? return : (return 1);
}

sub pp_grid_info {
    #---------------------------------------------------------------------------
    # Get grid details from SIZES job file.
    # Returns a hash with the following keys:
    #   bdx
    #   bdy
    #   lbnpt
    #   lbrow
    #   bzy
    #   bzx
    #   bplat
    #   bplon
    #
    # Call: %pp = pp_grid_info($job_files_directory);
    #     : %pp = pp_grid_info($grid_namelist);
    #---------------------------------------------------------------------------
    my $JOBDIR = shift;

    my %pp;
    if (-f $JOBDIR) {          
	# Get details from namelist
	%pp = read_namelist($JOBDIR);
	$pp{bdx}   = $pp{DELTA_LAMBDA_TARG};
	$pp{bdy}   = $pp{DELTA_PHI_TARG};
	$pp{lbnpt} = $pp{POINTS_LAMBDA_TARG}; 
	$pp{lbrow} = $pp{POINTS_PHI_TARG};
	$pp{bzy}   = $pp{PHI_ORIGIN_TARG};
	$pp{bzx}   = $pp{LAMBDA_ORIGIN_TARG};
	$pp{bplat} = $pp{PHI_POLE}; 
	$pp{bplon} = $pp{LAMBDA_POLE};
    } elsif (-d $JOBDIR) {
	# Get details from SIZES job file
	%pp = read_namelist("$JOBDIR/SIZES");	
	$pp{bdx}   = $pp{EWSPACEA};
	$pp{bdy}   = $pp{NSSPACEA};
	$pp{lbnpt} = $pp{ROW_LENGTH}; 
	$pp{lbrow} = $pp{P_ROWS}; 
	$pp{bzy}   = $pp{FRSTLATA};
	$pp{bzx}   = $pp{FRSTLONA};
	$pp{bplat} = $pp{POLELATA};
	$pp{bplon} = $pp{POLELONA};
    }
    $pp{bdy} *= -1 if ($ENV{VN} < 5);
    $pp{bzy} -= $pp{bdy};
    $pp{bzx} -= $pp{bdx};
    
    return %pp;
}

sub create_lam_containing_pp {
    #---------------------------------------------------------------------------
    # Create a LAM non-rotated pole PP field which contains an RCM field
    # modified by C.WANG to create LAM for any rectangular area in non-rotated
    # grid by using buffer with latitude/longitude by setting -r w,s,e,n
    # in callin ppregrid. C wang (6/12/2011)
    # 
    # Call: 
    #---------------------------------------------------------------------------
    my $dx     = shift;
    my $dy     = shift;
    my $lat    = shift;
    my $lon    = shift;
    my $buffer = shift;
    my $point_or_ll = shift; 
    my $ppin   = shift;  # HASH reference: PP hash

    # Default buffer size ('w,s,e,n')
    $buffer = 1    if (!defined $buffer);

    # Default lat-lon anchor
    $lat    = 90.0 if (!defined $lat);
    $lon    =  0.0 if (!defined $lon);
    my %this_ppin = %{$ppin};

    die "ERROR: Must define dx and dy for the target grid.\n" 
	if (!(defined $dx or defined $dy));

    my %ppref = &pp_set;

    # Absolute values of grid spacings
    my $adx = abs($dx);
    my $ady = abs($dy);

    my %buffer;
    my @buffer = split /\s*,\s*/, $buffer;
    my $west;
    my $south;
    my $east;
    my $north; 
    if ($point_or_ll >0 ) {   # When buffer is the no. of grid points each side of lam domain
      # Set buffer hash
      if (@buffer == 1) {
	@buffer = ($buffer[0]) x 4;
      } elsif (@buffer == 2 or @buffer == 3 or @buffer > 4) { 
 	die "ERROR: Buffer does not have 1 or 4 values: buffer=$buffer\n";
      }
      %buffer = ('w' => $buffer[0],   # size of mdi buffer to west  of ppin
        	 's' => $buffer[1],   # size of mdi buffer to south of ppin
		 'e' => $buffer[2],   # size of mdi buffer to east  of ppin
		 'n' => $buffer[3]);  # size of mdi buffer to north of ppin	
       
      foreach my $key (keys %buffer) {
	unless (natural($buffer{$key})) {
	  die "ERROR: Buffer entry \'$key\' is not a non-negative integer: $key=$buffer{$key}\n"; 
	}	    
      }
      ($west,$south,$east,$north) = rcm_gridbox_edge_extrema_in_reg_coords(\%$ppin);
      $west =$west  -($adx *$buffer{w});
      $south=$south -($ady *$buffer{s});
      $east =$east  +($adx *$buffer{e});
      $north=$north +($ady *$buffer{n});
    }else{ #buffer is w,s,e,n
      if (@buffer != 4) {
        die "ERROR:w,s,e,n Buffer does not have  4 values: buffer=$buffer\n";
      }
      %buffer = ('w' => $buffer[0],   # 
                 's' => $buffer[1],   # 
	         'e' => $buffer[2],   # 
	         'n' => $buffer[3]);  #    
      $west  = $buffer{w};
      $south = $buffer{s};
      $east  = $buffer{e};
      $north = $buffer{n}; 
    }	    
   
    my $w = $lon;
    while ($w <  $west) { $w += $adx };
    while ($w >= $west) { $w -= $adx };

    my $e = $w;
    while ($e <= $east) { $e += $adx };

    my $s = $lat;
    while ($s <  $south) { $s += $ady };
    while ($s >= $south) { $s -= $ady };

    my $n = $s;
    while ($n <= $north) { $n += $ady };


    # Make sure that $north and $south are between -90 and 90
    while ($s < -90.0) {
	$s += $ady;
	$buffer{s}--;
    }
    while ($n > 90.0) {
	$n -= $ady;
	$buffer{n}--;
    }


    # Check for 360 degree wraparound
    my $wraparound360;
    my $global;
    if ($e - $w + $adx >= 360.0) {
	if (360/$adx != int(360/$adx)) {
	    die "ERROR: Limited area target grid has 360 degree wrap around\n".
		"       but 360 is not an integer multiple of dx:\n".
		"       dx=$dx, 360/dx=",360/$dx,"\n";
	}
	$buffer{w} = 'n/a';
	$buffer{e} = 'n/a';
	$w = $lon;
	$e = $w + 360 - $adx;

	$wraparound360 = 1;
	$global = 1 if (abs($n - $s) >= 180.0);
    }

    # Check that $e > $w
    die "ERROR: e <= w: e = $e, w = $w\n" if ($e <= $w);    
    # Check that $n > $s
    die "ERROR: n <= s: n = $n, s = $s\n" if ($n <= $s);

    # Set lbnpt lbrow lblrec lbhem lbcode
    $ppref{lbnpt}  = int(($e - $w)/$adx) + 1;
    $ppref{lbrow}  = int(($n - $s)/$ady) + 1;
    $ppref{lblrec} = $ppref{lbrow} * $ppref{lbnpt};
    $ppref{lbnrec} = $ppref{lblrec};
    $ppref{lbcode} = 1;
    $ppref{lbhem}  = $global        ? 0 : 
	$wraparound360 ? 4 : 3;

    # Set lbfc lbuser4 lbrel
    $ppref{lbfc}    = 0;
    $ppref{lbuser4} = 0;

    # Set bzx bzy bdx bdy bplat bplon
    if ($ENV{VN} < 5) {
	$ppref{bdx} = $adx;  
	$ppref{bdy} = -1 * $ady;
	$ppref{bzy} = $n - $ppref{bdy};
	$ppref{bzx} = $w - $ppref{bdx};
    } else {
	$ppref{bdx} = $adx;  
	$ppref{bdy} = $ady;  
	$ppref{bzy} = $s - $ppref{bdy};
	$ppref{bzx} = $w - $ppref{bdx};
    }
    $ppref{bzx}    = shift_180($ppref{bzx});
    $ppref{bplat}  = 90.0;
    $ppref{bplon}  = 0.0;

    # Set data to zero
    $ppref{data} = [ (0.0) x $ppref{lblrec} ];

    return %ppref;
}

sub create_global_pp {

    # Create a global PP file

    my $dx   = shift; # Sets $nx if $nx=undef
    my $dy   = shift; # Sets $nx if $ny=undef
    my $lat  = shift; # 
    my $lon  = shift; # 
    my $nx   = shift; # Sets $dx if $dx=undef
    my $ny   = shift; # Sets $dy if $dy=undef
    my $grid = shift; # 
    my $ppin = shift; # Ignored  if $grid is set.

    # Set grid type if one has not been provided
    if (!defined $grid) {
	if (defined $ppin) {
	    if ($ppin->{lbuser4} <= 0) {
		$grid = 1;
	    } else {
		my %STASHmaster = read_STASHmaster();	
		$grid = $STASHmaster{$ppin->{lbuser4}}[$Sm{Grid}] || 1;
		$grid =~ s/\s+//g;
	    }
	} else {
	    confess "ERROR: Must define \$grid or \%ppin\n";    
	}
    }
    
    # Adjust $lat and $lon for the grid type
  CASE: {
      if ($grid eq '11') { $lat += $dy/2;  $lon += $dx/2; last CASE2; } # Arakawa-B   wind grid
      if ($grid eq '12') { $lat += $dy/2;  $lon += $dx/2; last CASE2; } # Arakawa-B   wind grid
      if ($grid eq '13') { $lat += $dy/2;  $lon += $dx/2; last CASE2; } # Arakawa-B   wind grid
      if ($grid eq '18') { $lat += $dy/2;                 last CASE2; } # Arakawa-C u-wind grid
      if ($grid eq '19') {                 $lon += $dx/2; last CASE2; } # Arakawa-C v-wind grid
  } # CASE:
    
    # Set $dx, $dy, $nx, $ny
    if (defined $dx and defined $dy) {
  	$nx = 360/abs($dx);
	$ny = 180/abs($dy) + 1;
    } else {
 	$dx = 360/$nx;
	$dy = 180/($ny-1);
        $dy *= -1 if (defined $ppin and $ppin->{bdy} < 0);
    } 

    # Adjust $lat to be on first row of grid
    if ($dy < 0) {
	while ($lat <  90.0) { $lat -= $dy };
	while ($lat >  90.0) { $lat += $dy };
    } else {			     	
	while ($lat > -90.0) { $lat -= $dy };
	while ($lat < -90.0) { $lat += $dy };
    }

    # Set PP header and data
    my %ppref = &pp_set;
    $ppref{bdx}   = $dx;
    $ppref{bdy}   = $dy;
    
    $ppref{lbnpt} = $nx;
    $ppref{lbrow} = (abs($lat) == 90.0) ? $ny : $ny-1;

    $ppref{lbrel}   = 2;
    $ppref{bplat}   = 90.0;               
    $ppref{bplon}   = 0.0;                
    $ppref{lbhem}   = 0;
    $ppref{lbcode}  = 1;
    $ppref{lbfc}    = 0;
    $ppref{lbuser4} = 0;
    $ppref{lblrec}  = $ppref{lbnpt} * $ppref{lbrow};

    $ppref{bzx}     = $lon - $ppref{bdx};
    $ppref{bzy}     = $lat - $ppref{bdy};
    $ppref{bmdi}    = -1073741800;

    # Set data array to zero
    $ppref{data} = [ ($ppref{bmdi}) x $ppref{lblrec} ];
 
    # Error check
    if ($ppref{lbnpt} != int($ppref{lbnpt})) {
	die "ERROR: 360 is not an integer multiple of dx: dx=$ppref{bdx}, 360/dx=".360/$ppref{bdx}."\n";
    }
    if ($ppref{lbrow} != int($ppref{lbrow})) {
	die "ERROR: 180 is not an integer multiple of dy: dy=$ppref{bdy}, 180/dy=".180/$ppref{bdy}."\n";
    }    

    # Add an extra useful field
    $ppref{grid} = $grid;

    return %ppref;
}

sub skip_pp_data {

    #--------------------------------------------------------------------------
    # Having read a PP header (up to and including the 3rd block control word),
    # skip the data and the 4th block control word and be ready to read the
    # next PP field.
    #
    # Call: %pp = skip_pp_data(\*FILEHANDLE, \%pp);
    #--------------------------------------------------------------------------
    my $FH   = shift;  # An open filehandle
    my $pp   = shift;  # HASH reference: PP hash

    # Skip to beginning of next PP field
    sysseek $FH, $pp->{lblrec} * WORD_SIZE + MARKER_SIZE, 1;
}

#sub machine_endian {
#    my $bo = $Config{byteorder};
#    if ($bo =~ /^1234/) {
#	return 'big' 
#	    if (defined $ENV{F_UFMTENDIAN} and $ENV{F_UFMTENDIAN} eq 'big');
#	return 'little';
#    } elsif ($bo =~ /4321$/) {
#	return 'little' 
#	    if (defined $ENV{F_UFMTENDIAN} and $ENV{F_UFMTENDIAN} eq 'little');
#	return 'big';
#    }  
#    # Can not determine byte order
#    return;
#}

sub subregion_ppref {
    #--------------------------------------------------------------------------
    # 
    #
    # Call: %region = subregion_ppref(\%ppref);
    #     : %region = subregion_ppref(\*FILEHANDLE);
    #     : %region = subregion_ppref($filename);
    #--------------------------------------------------------------------------
    my $arg = shift; # One of: 1) HASH reference: PP hash
                     #         2) An open filehandle
                     #         3) A PP file name

    my %ppref;
    if (ref($arg) eq 'HASH') {
	%ppref = %$arg;
    } else {
	%ppref = read_pp($arg);
    }

    my %region;

    $region{pole_lat} = $ppref{bplat};
    $region{pole_lon} = $ppref{bplon};

    $region{west}  = $ppref{bzx} + 0.5*$ppref{bdx};
    $region{east}  = $ppref{bzx} + ($ppref{lbnpt}+0.5)*$ppref{bdx};
    $region{north} = $ppref{bzy} + 0.5*$ppref{bdy};
    $region{south} = $ppref{bzy} + ($ppref{lbrow}+0.5)*$ppref{bdy};
    
    if ($ppref{bdx} < 0) {
	my $temp       = $region{west};
	$region{west}  = $region{east};
	$region{east}  = $temp;
    }
    if ($ppref{bdy} > 0) {
	my $temp       = $region{south};
	$region{south} = $region{north};
	$region{north} = $temp;
    }

    $region{west}  = shift_360($region{west});
    $region{east}  = make_x_to_360($region{west},$region{east});

    $region{south} = -90.0 if ($region{south} < -90.0);
    $region{north} =  90.0 if ($region{north} >  90.0);
    
    return %region;
}  # sub subregion_ppref


sub global {
    #--------------------------------------------------------------------------
    # Determine if a PP field is global
    # 
    # Call: $global = global(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift; # HASH reference: PP hash

    # Return 1 if field is global
    return 1 if ($pp->{lbhem} == 0);    

    # Return undef if field is not global
    return;
}

sub pp_longitude {
    #--------------------------------------------------------------------------
    # Find the unique longitudes of gridbox centres of a PP field
    #
    # Call: @longitude = pp_longitude(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift; # HASH reference: PP hash

    my @longitude;
    foreach my $x (1..$pp->{lbnpt}) {
	push @longitude, $pp->{bzx} + $x*$pp->{bdx};  
    }    

    return shift_360(@longitude);
}

sub pp_latitude {
    #--------------------------------------------------------------------------
    # Find the unique latidudes of gridbox centres of a PP field
    #
    # Call: @latitude = pp_latitude(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift; # HASH reference: PP hash

    my @latitude;
    foreach my $y (1..$pp->{lbrow}) {
	push @latitude, $pp->{bzy} + $y*$pp->{bdy};              
    }

    return @latitude;
}

sub pp_time1 {
    #--------------------------------------------------------------------------
    # Return a time hash for the time of the first field in a PP file
    #
    # Call: %time = pp_time1(\*FILEHANDLE);
    #     : %time = pp_time1($filename);
    #--------------------------------------------------------------------------
    my $arg = shift; # An open filehandle or a PP file name

    my $FH = file_handle($arg, 'read', 'binary');

    # Read PP header
    my %pp;
    while (1) {	
    	my %pp0 = read_pp_header(\*$FH);
        %pp = %pp0;
	last if ($pp0{lbuser4} == 10 || $pp0{lbuser4} == 2 ||
                 $pp0{lbuser4} == 3 || $pp0{lbuser4} == 30111 ||
                 $pp0{lbuser4} == 4 || $pp0{lbuser4} == 16004 ||
                 $pp0{lbuser4} == 409 || $pp0{lbuser4} == 1 || 
                 $pp0{lbuser4} == 30417 );
													    	
    	skip_pp_data(\*$FH, \%pp);

    }


    close $FH;

    return (year   => $pp{lbyr},
	    month  => $pp{lbmon},
	    day    => $pp{lbdat},
	    hour   => $pp{lbhr},
	    minute => $pp{lbmin},
	    second => 0);
}

sub pp_time2 {
    #--------------------------------------------------------------------------
    # Return a time hash for the time of the last field in a PP file
    #
    # Call: %time = pp_time2(\*FILEHANDLE);
    #     : %time = pp_time2($filename);
    #--------------------------------------------------------------------------
    my $arg = shift; # An open filehandle or a PP file name

    my $FH = file_handle($arg, 'read', 'binary');

    my %pp;

    # Read PP header
    while (1) {	
	my %pp0 = read_pp_header(\*$FH);
	last if (! %pp0);
	%pp = %pp0;
	skip_pp_data(\*$FH, \%pp)
    }
    close $FH;

    return (year   => $pp{lbyr},
	    month  => $pp{lbmon},
	    day    => $pp{lbdat},
	    hour   => $pp{lbhr},
	    minute => $pp{lbmin},
	    second => 0);
}

sub pp_vtime {
    #--------------------------------------------------------------------------
    # 
    #
    # Call: %vtime = pp_vtime(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift;

    return (year   => $pp->{lbyr},
	    month  => $pp->{lbmon},
	    day    => $pp->{lbdat},
	    hour   => $pp->{lbhr},
	    minute => $pp->{lbmin},
	    second => 0);
}

sub pp_dtime {
    #--------------------------------------------------------------------------
    # 
    #
    # Call: %dtime = pp_dtime(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift;

    return (year   => $pp->{lbyrd},
	    month  => $pp->{lbmond},
	    day    => $pp->{lbdatd},
	    hour   => $pp->{lbhrd},
	    minute => $pp->{lbmind},
	    second => 0);
}

sub pp_period {
    #--------------------------------------------------------------------------
    # 
    #
    # Call: %period = pp_period(\%pp);
    #--------------------------------------------------------------------------
    my $pp = shift;
    # Period ('end time' minus 'start time')
    
    # Retun 0 for instantaneous data, for which 'end time' could well
    # be 'data time', etc.
    return 0 if ($pp->{lbproc} <= 0);
	      
    my %vtime = pp_vtime($pp);
    my %dtime = pp_dtime($pp);
    my $LCAL360 = (($pp->{lbtim} % 10) == 2) ? 1 : undef;
    return time_difference(\%vtime, \%dtime, $LCAL360);
}
1;