#! /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;