#!/usr/bin/perl -w
# 2012 Public Domain Wesley Ebisuzaki
#
# makes a GrADS control file for grib2 files
#
# requires wgrib2 and Perl5
#
# usage: alt_g2ctl [options] [grib file] [optional index file] >[control file]
#
# note: this script does not make the index file .. you have to run alt_gmp
#
# Analyses: (using initial time)
#
# $ alt_g2ctl -0 example.grib >example.ctl
# $ alt_gmp example.ctl
#
# Forecasts: (using verifiation time)
#
# $ alt_g2ctl example.grib >example.ctl
# $ alt_gmp example.ctl
#
# bugs:
# many
#
# requires newer wgrib2
#
# documentation: http://www.cpc.ncep.noaa.gov/products/wesley/alt_g2ctl_gmp.html
#
#
# added output for rotated LatLon grids
# Helmut P. Frank, Helmut.Frank@dwd.de
# Fri Sep 14 13:54:00 GMT 2001
# -ts, -lc options: Ag Stephens, BADC 3/2003
#
# requires grads 2.0a3+ - thinned grid update
#
$version="1.0.8";
use POSIX;
use Math::Trig qw(deg2rad rad2deg);
# ***** if wgrib2 is not on path, add it here
$wgrib2='wgrib2';
$wgrib2_flags='-npts -set_ext_name 1 -end_FT -ext_name -lev';
$wgrib2_inv=".invd01";
$file="";
$index="";
$pdef="";
$calendar="";
$timestep="";
# $nearest_neighbor="";
$kludge="";
$template="";
$pdef_nearest=1;
$raw="no";
$profile="mb";
$profile_sort='reverse';
$update='0';
$big='0';
$nthreads=4;
$short_varname=0;
$match="";
# set pdef_offset = 0 for old code, 1 for new code
$pdef_offset=0;
for ($i = 0; $i <= $#ARGV; $i++) {
$_ = $ARGV[$i];
SWITCH: {
/^-verf$/ && do { last SWITCH; };
/^-0$/ && do { $wgrib2_flags = '-npts -set_ext_name 1 -T -ext_name -lev'; $wgrib2_inv = ".invd03"; last SWITCH };
/^-0t$/ && do { $wgrib2_flags = '-npts -set_ext_name 1 -T -ext_name -ftime -lev'; $wgrib2_inv = ".invd02"; last SWITCH };
/^-b$/ && do { $wgrib2_flags = '-npts -set_ext_name 1 -start_FT -ext_name -lev'; $wgrib2_inv = ".invd04"; last SWITCH };
/^-bt$/ && do { $wgrib2_flags = '-npts -set_ext_name 1 -start_FT -ext_name -lev'; $wgrib2_inv = ".invd05"; last SWITCH };
/^-e$/ && do { $wgrib2_flags='-npts -set_ext_name 1 -end_FT -ext_name -lev'; $wgrib2_inv=".invd01"; last SWITCH };
/^-et$/ && do { $wgrib2_flags='-npts -set_ext_name 1 -end_FT -ext_name -ftime -lev'; $wgrib2_inv=".invd01"; last SWITCH };
/^-update$/ && do { $update="1"; last SWITCH };
/^-(prs|mb)$/ && do { $profile = "mb"; $profile_sort='reverse'; last SWITCH };
/^-no_profile$/ && do { $profile = ""; $profile_sort=''; last SWITCH };
/^-iso$/ && do { $profile = "iso"; $profile_sort='' ; last SWITCH };
/^-bsl$/ && do { $profile = "bsl"; $profile_sort='reverse' ; last SWITCH };
/^-dsl$/ && do { $profile = "dsl"; $profile_sort='reverse' ; last SWITCH };
/^-hyp$/ && do { $profile = "hyp"; $profile_sort='' ; last SWITCH };
/^-hply$/ && do { $profile = "hply"; $profile_sort='' ; last SWITCH };
/^-sig$/ && do { $profile = "sig"; $profile_sort='' ; $profile_sort='reverse'; last SWITCH };
/^-365$/ && do { $calendar="365"; last SWITCH; };
/^-ts(\d+\w+)/ && do { $timestep=$1; last SWITCH; };
/^-kludge$/ && do { $kludge="on"; last SWITCH; };
# /^-sub$/ && do { push @from, $ARGV[$i+1]; push @to, $ARGV[$i+2]; $i+=2; last SWITCH; };
/^-del$/ && do { push @from, $ARGV[$i+1]; push @to,""; $i+=1; last SWITCH; };
/^-raw$/ && do { $raw="yes" ; last SWITCH };
/^-pdef_linear$/ && do { $pdef_nearest=0 ; last SWITCH };
/^-nthreads$/ && do { $nthreads=$ARGV[++$i]; ; last SWITCH };
/^-wgrib2$/ && do { $wgrib2=$ARGV[++$i]; last SWITCH };
/^-short$/ && do { $short_varname = 1; last SWITCH };
/^-match$/ && do { $match = $match . " -match \"$ARGV[++$i]\""; last SWITCH };
/^-not$/ && do { $match = $match . " -not \"$ARGV[++$i]\""; last SWITCH };
# /^-big$/ && do { $big='1'; ; last SWITCH };
/^-/ && do { print STDERR "unknown option: $_\n"; exit 8; };
if ($file eq "") {
$file = $_;
}
elsif ($index eq "") {
$index = $_;
}
else {
$pdef = $_;
}
}
}
# print STDERR "file=$file index=$index\n";
if ($file eq "") {
if ($#ARGV >= 0) {
print STDERR "*** missing grib file ***\n\n\n";
}
print STDERR "$0 $version wesley ebisuzaki\n";
print STDERR " makes a Grads alt control file for grib2 files\n";
print STDERR " usage: $0 [options] [grib_file] [optional index file] [optional pdef file] >[ctl file]\n";
print STDERR " -0 .. use analysis time\n";
print STDERR " -0t .. use analysis time + fhour\n";
print STDERR " -b .. use use start of ave/acc period or fcst time\n";
print STDERR " -bt .. use use start of ave/acc period or fcst time + fhour\n";
print STDERR " -e .. use use end of ave/acc period or fcst time (default\n";
print STDERR " -et .. use use end of ave/acc period or fcst time + fhour\n";
# print STDERR " -big .. for files > 2GB (only works on 64-bit machines/codes\n";
print STDERR " -update .. alt_gmp will be in fast update mode\n";
print STDERR " -nthreads N .. number of threads used by alt_gmp (default=1)\n";
print STDERR " -wgrib2 EXE .. replace wgrib2 by EXE\n";
print STDERR " -short .. remove comments and shorten variables names\n";
print STDERR " -match X .. only match X (regex) from wgrib2 inv, can be repeated\n";
print STDERR " -not X .. not X (regex) from wgrib2 inv, can be repeated\n";
# print STDERR " -sub \"A\" \"B\" .. substiute variables names, A -> B, use regular expressions\n";
# print STDERR " .. any number of -sub options can be used\n";
print STDERR " -prs .. pressure (mb) vertical coordinates (default)\n";
print STDERR " -iso .. pot temp (K) vertical coordinates\n";
print STDERR " -dsl .. below sea level (m) vertical coordinates\n";
print STDERR " -bsl .. below sea level (m) vertical coordinates (obsolete)\n";
print STDERR " -hyp .. hybrid pressure level\n";
print STDERR " -hply .. hybrid pressure layer\n";
print STDERR " -sig .. sigma (0..1) vertical coordinates\n";
print STDERR " -no_profile .. no vertical coordinates\n";
print STDERR " -365 .. 365 day calendar\n";
print STDERR " -ts[timestep] .. set timestep for individual time files (e.g. -ts6hr)\n";
print STDERR " -lc .. set lowercase option for parameter names\n";
print STDERR " -pdef_linear .. linear interpolation for thinned grids\n";
print STDERR " -raw .. raw grid\n";
print STDERR "\n";
print STDERR "Note 1: the index file will be generated by the alt_gmp program, default: grib_file.idx\n";
print STDERR "Note 2: the pdef file is only generated for thinned lat-lon grids, default: grib_file.pdef\n";
print STDERR "Note 3: template options supported: %y4 %y2 %m2 %m1 %d2 %d1 %h2 %h1 %n2 %f2 %f3 %e %j3\n";
print STDERR " %iy4 %iy2 %im2 %im1 %id2 %id1 %ih2 %ih1 %in2 %fn2 %fhn %fdhn\n";
exit 8;
}
$_ = $file;
if (/%y4/ || /%y2/ || /%m2/ || /%m1/ || /%d2/ || /%d1/ || /%h2/ || /%h1/ || /%n2/
|| /%f2/ || /%j3/ || /%f3/ || /%e/ || /%iy2/ || /%iy4/ || /%im1/ || /%im2/ ||
/%id1/ || /%id2/ || /%ih1/ || /%ih2/ || /%in2/ || /%fn2/ || /%fhn/ || /%fdhn/ ) { $template='on'; }
if ($^O eq 'MSWin32') {
$ListA="$ENV{TEMP}\\h$$.tmp";
$TmpFile="NUL";
unlink $ListA;
$sys="win";
}
else {
$TmpFile="/dev/null";
$j=0;
while (-f "/tmp/g.$$.$j.tmp") {
$j++;
}
$ListA="/tmp/g.$$.$j.tmp";
$sys="unix";
}
# ctlfilename = name used by control file (different for template option(
# file = file name (of first matching file)
$ctlfilename=$file;
# inventory of All records
if ($template eq "on") {
$gfile=$file;
# if windows change backslashes to regular slashes
if ($sys eq 'win') { $gfile =~ s=\\=/=g; }
# change GrADS templates to "capture" regular expressions
$gfile =~ s/%y4/(\?<year>\\d{4})/g;
$gfile =~ s/%iy4/(\?<year>\\d{4})/g;
$gfile =~ s/%y2/(\?<year>\\d{2})/g;
$gfile =~ s/%iy2/(\?<year>\\d{2})/g;
$gfile =~ s/%m2/(\?<month>\\d{2})/g;
$gfile =~ s/%im2/(\?<month>\\d{2})/g;
$gfile =~ s/%m1/(\?<month>\\d{1,2})/g;
$gfile =~ s/%im1/(\?<month>\\d{1,2})/g;
$gfile =~ s/%d2/(\?<day>\\d{2})/g;
$gfile =~ s/%id2/(\?<day>\\d{2})/g;
$gfile =~ s/%d1/(\?<day>\\d{1,2})/g;
$gfile =~ s/%j3/(\?<day>\\d{3})/g;
$gfile =~ s/%id1/(\?<day>\\d{1,2})/g;
$gfile =~ s/%h2/(\?<hour>\\d{2})/g;
$gfile =~ s/%ih2/(\?<hour>\\d{2})/g;
$gfile =~ s/%h1/(\?<hour>\\d{1,2})/g;
$gfile =~ s/%ih1/(\?<hour>\\d{1,2})/g;
$gfile =~ s/%n2/(\?<min>\\d{1,2})/g;
$gfile =~ s/%in2/(\?<min>\\d{1,2})/g;
$gfile =~ s/%h3/(\?<hour>\\d{3})/g;
$gfile =~ s/%f2/(\?<fhr>\\d{2,8})/g;
$gfile =~ s/%f3/(\?<fhr>\\d{3,8})/g;
$gfile =~ s/%fn2/(\?<fmin>\\d{2})/g;
$gfile =~ s/%fhn/(\?<fhn>\\d{4})/g;
$gfile =~ s/%fdhn/(\?<fdhn>\\d{6})/g;
$dir=$gfile;
$dir =~ s=(/*)[^/]*$=$1=;
$gfile =~ s/^$dir//;
if (($gfile =~ s/%/%/g) != 0) {
print STDERR "unrecognized GrADS template";
exit 8;
}
if ($dir eq "") {
opendir(DIR,'.');
}
else {
opendir(DIR,$dir);
}
@allfiles = grep /^$gfile$/, readdir DIR;
closedir DIR;
if ($#allfiles <= -1 ) {
print STDERR "\nError: could not find any files in directory: $dir\n";
exit 8;
}
# find 1st, 2nd and last files in chronological order
for ($i = 0; $i <= $#allfiles; $i++) {
$tfile = $allfiles[$i];
$tfile =~ m/$gfile/;
# find time_index ..
$t_index = 0;
if (defined($+{min})) { $t_index += $+{min}; }
if (defined($+{fhn})) { $t_index += $+{fhn}; }
if (defined($+{fdhn})) { $t_index += $+{fdhn}; }
if (defined($+{hour})) { $t_index += 100*$+{hour}; }
if (defined($+{fhr})) { $t_index += 100*$+{fhr}; }
if (defined($+{day})) { $t_index += 10000*$+{day}; }
if (defined($+{month})) { $t_index += 1000000*$+{month}; }
if (defined($+{year})) { $t_index += 100000000*$+{year}; }
if (!defined($file1)) { # 1st time through
$file1=$tfile;
$file2=$tfile;
$filelast=$tfile;
$index1 = $t_index;
$index2 = $t_index;
$indexlast = $t_index;
}
elsif ( $t_index < $index1) {
$index2 = $index1;
$file2 = $file1;
$index1 = $t_index;
$file1 = $tfile;
}
elsif ($t_index > $indexlast) {
if ($index1 == $index2) {
$index2 = $t_index;
$file2 = $tfile;
}
$indexlast = $t_index;
$filelast = $tfile;
}
elsif ($t_index < $index2 && $t_index > $index1) {
$index2 = $t_index;
$file2 = $tfile;
}
}
# used in grid section
$file = "$dir$file1";
# print STDERR "1,2,n= $file1 $file2 $filelast\n";
# make inventory of times 1, 2 and last
system "$wgrib2 $match -lev0 $wgrib2_flags \"$dir$file1\" >$ListA";
if ($file1 ne $file2) {
system "$wgrib2 $match -lev0 $wgrib2_flags \"$dir$file2\" >>$ListA";
}
if ($file2 ne $filelast) {
system "$wgrib2 $match -lev0 $wgrib2_flags \"$dir$filelast\" >>$ListA";
}
}
else {
system "$wgrib2 $match -lev0 $wgrib2_flags \"$file\" >$ListA";
}
if ( ! -s $ListA ) {
if ($match ne '') {
print STDERR "Control file was not generated, variables not found\n";
print STDERR " $match used\n";
}
else {
print STDERR "Big problem:\n";
print STDERR " either $file is missing or not a grib file\n";
print STDERR " or wgrib2 is not on your path or wgrib2 is too old\n";
print STDERR " or can not write to $ListA (full disk or permissions)\n";
}
unlink $ListA;
exit 8;
}
# make table of dates and variables
scan_ListA();
unlink $ListA;
@sdates=sort keys(%dates);
# number of time 1 or greater
$ntime=$#sdates + 1;
$time=$sdates[0];
$year = substr($time,0,4);
$mo = substr($time,4,2);
$day = substr($time,6,2);
$hour = substr($time,8,2);
$minute = substr($time,10,2);
if ($mo < 0 || $mo > 12) {
print "illegal date code $time\n";
exit 8;
}
$month=substr("janfebmaraprmayjunjulaugsepoctnovdec",$mo*3-3,3);
if ($ntime > 1) {
$year1 = substr($sdates[1],0,4);
$mo1 = substr($sdates[1],4,2);
$day1 = substr($sdates[1],6,2);
$hour1 = substr($sdates[1],8,2);
$minute1 = substr($sdates[1],10,2);
$year_last = substr($sdates[$#sdates],0,4);
$mo_last = substr($sdates[$#sdates],4,2);
$day_last = substr($sdates[$#sdates],6,2);
$hour_last = substr($sdates[$#sdates],8,2);
$minute_last = substr($sdates[$#sdates],10,2);
}
# ---------------intro------------------------------------
if ("$index" eq "") {$index="$file.idx";}
if ("$pdef" eq "") { $pdef = "$file.pdef";}
if ($sys eq "unix") {
$caret1 = (substr($file,0,1) eq "/") ? "" : '^';
$caret2 = (substr($index,0,1) eq "/") ? "" : '^';
$caret3 = (substr($pdef,0,1) eq "/") ? "" : '^';
}
else {
$caret1 = (substr($file,1,1) eq ":") ? "" : '^';
$caret2 = (substr($index,1,1) eq ":") ? "" : '^';
$caret3 = (substr($pdef,1,1) eq ":") ? "" : '^';
}
print "dset $caret1$ctlfilename\nindex $caret2$index\n";
print "undef 9.999E+20\ntitle $file\n* produced by alt_g2ctl v$version, use alt_gmp to make idx file\n";
print "* command line options: @ARGV\n";
print "* alt_gmp options: update=$update\n";
print "* alt_gmp options: nthreads=$nthreads\n";
print "* alt_gmp options: big=$big\n";
print "* alt_gmp options: match=$match\n";
print "* wgrib2 inventory flags: $wgrib2_flags\n";
print "* wgrib2 inv suffix: $wgrib2_inv\n";
# ------------------- grid -----------------------
$griddef = `$wgrib2 "$file" -one_line -d 1 -nxny -grid -vector_dir`;
$_=$griddef;
$_ = $griddef;
$_ =~ s/^[^(]*//;
$_ =~ s/:.*//;
/\((\S*) x -*(\d*)\)/;
$nx=$1;
$ny=$2;
$_=$griddef;
$t = substr($griddef,0,240);
print "* griddef=$t\n";
print "dtype grib2\n";
if ($template eq "on") { print "options template\n"; }
if ($raw eq 'yes') {
print "xdef $nx linear 0 0.1\n";
print "ydef $ny linear 0 0.1\n";
}
elsif (/: Gaussian grid:/) {
/ lon (\S*) to (\S*)/;
$lon0=$1;
$lon1=$2;
if ($lon1 <= $lon0) { $lon1 += 360.0; }
$dlon = ($lon1 - $lon0) / ($nx - 1);
/ lat (\S*) to (\S*)/;
if ($1 < $2) {
$lat0=$1;
$lat1=$2;
}
else {
$lat0=$2;
$lat1=$1;
}
/number of latitudes between pole-equator=(\S*)/;
$nGauLats=2*$1;
print "* ny=$ny nlat=$nGauLats lat0=$lat0 lat1=$lat1\n";
print "xdef $nx linear $lon0 $dlon\n";
print "ydef $ny levels\n";
&print_range_gau_lats;
}
elsif (/ thinned global Gaussian grid:/) {
/ input (\S*)/;
$scan = $1;
if ($scan ne 'WE:NS' && $scan ne 'WE:SN') { print "\n* **** unsupported scan mode $scan\n"; }
/#grid points by latitude: (.*)/;
$list = $1;
/#points=(\d*)/;
$npnts=$1;
$i = 1; $nx = 0;
foreach $t (split(' ',$list)) {
$t =~ s/:.*//;
$nlon[$i++] = $t;
if ($nx < $t) { $nx = $t; }
}
$dx = 360.0 / $nx;
$t = 0;
printf "xdef $nx linear 0 $dx\n";
print "ydef $ny levels\n";
&print_gau_lats;
# now to create the pdef file
open (PDEF, ">$pdef") or die "Could not open pdef file for write: $pdef";
binmode PDEF;
$regional_thinned_grid=0;
if ($pdef_nearest == 0) {
pdef_linear();
} else {
pdef_near();
}
close(PDEF);
}
elsif (/ thinned (global|regional) lat-lon grid:/) {
/ lat *(\S*) to (\S*) by (\S*) /;
$lat0 = $1;
$lat1 = $2;
$dy = $3;
/ lon (\S*) to (\S*) with variable spacing/;
$lon0=$1;
$lon1=$2;
/\(-1 x (\d*)\)/;
$ny = $1;
/ #points=(\d*) /;
$npnts = $1;
/ input (\S*)/;
$scan = $1;
if ($scan ne 'WE:NS' && $scan ne 'WE:SN') { print "\n* **** unsupported scan mode $scan\n"; }
if ($lat0 > $lat1) {
$yrev = 1;
print "ydef $ny linear $lat1 ", abs($dy), "\n"
}
else {
$yrev = 0;
print "ydef $ny linear $lat0 ", abs($dy), "\n"
}
$t=$_;
$_ =~ s/^.*points by latitude: //;
$_ =~ s/:.*//;
$list=$_;
$i = 1;
foreach $t (split(' ',$list)) {
$nlon[$i++] = $t;
}
if ($lon1 <= $lon0) { $lon1 += 360.0; }
#
# check if global
#
$nx = $nlon[$ny];
$dx = ($lon1 - $lon0) / ($nx - 1);
$dx=$dx*$nx;
$regional_thinned_grid = 1;
if (abs($dx-360.0) < 0.01) {$regional_thinned_grid=0; }
#
# figure out dx and nx for the projection grid
#
$nx=$nlon[1];
if ($nx < $nlon[$ny]) { $nx = $nlon[$ny]; }
$dx = ($lon1 - $lon0) / ($nx - 1);
print "xdef $nx linear $lon0 $dx\n";
open (PDEF, ">$pdef") or die "Could not open pdef file for write: $pdef";
binmode PDEF;
$regional_thinned_grid=1;
if ($pdef_nearest == 0) {
pdef_linear();
} else {
pdef_near();
}
close(PDEF);
}
elsif (/: lat-lon grid:/) {
/ lat (\S*) to (\S*) by (\S*)/;
$lat0=$1;
$lat1=$2;
$dlat=$3;
/ lon (\S*) to (\S*) by (\S*)/;
$lon0=$1;
# $lon1=$2;
$dlon=$3;
if ($lat0 > $lat1) {
print "ydef $ny linear $lat1 ", abs($dlat), "\n"
}
else {
print "ydef $ny linear $lat0 ", abs($dlat), "\n"
}
if ($lon0 + ($nx-1) * $dlon > 360.0) { $lon0 -= 360.0; }
print "xdef $nx linear $lon0 $dlon\n";
}
elsif (/ Mercator grid: /) {
# beta: mercator
# scan modes .. assumes west to east
/lat *(\S*) to (\S*) /;
$lat1 = $1;
$lat2 = $2;
/lon (\S*) to (\S*) /;
$lon1 = $1;
$lon2 = $2;
/grid: \((\d*) x (\d*)\) /;
$nx = $1;
$ny = $2;
$dlon = $lon2 - $lon1;
if ($dlon <= 0) {
$dlon = $dlon + 360;
}
$dlon = $dlon / ($nx - 1);
if ($lon1 + ($nx-1) * $dlon > 360.0) { $lon1 -= 360.0 }
print "xdef $nx linear $lon1 $dlon\n";
if ($lat1 > $lat2) {
# print "options yrev\n";
$t = $lat2;
$lat2 = $lat1;
$lat1 = $t;
}
print "ydef $ny levels\n";
$i = 0;
$n1 = log(tan((45+$lat1/2)*3.1415927/180));
$n2 = log(tan((45+$lat2/2)*3.1415927/180));
$dy = ($n2 - $n1) / ($ny - 1);
while ($i < $ny) {
$nn = $n1 + $dy * $i;
$lat = (atan(exp($nn))*180/3.1415927-45)*2;
printf ("%9.4f ", $lat);
$i++;
if ($i % 7 == 0) { print "\n"; }
}
if ($i % 7 != 0) { print "\n"; }
}
elsif (/ Lambert Conformal: /) {
/ Lat1 (\S*) Lon1 (\S*) LoV (\S*) LatD \S* Latin1 (\S*) Latin2 (\S*) /;
$lat1 = $1;
$lon1 = $2;
$lov = $3;
$latin1 = $4;
$latin2 = $5;
/Pole \((\S*) x (\S*)\) Dx (\S*) m Dy (\S*) m/;
$nx = $1;
$ny = $2;
$dx = $3;
$dy = $4;
$t='lcc';
if (/:winds\(grid\)/) { $t='lccr' ; }
if ($dx == 0 || $dy == 0) {
print STDERR "Problem: DX or DY is missing! Cannot make a control file\n";
exit;
}
&domain;
# make sure than lon1 is within 180 degrees of lov
# lov, lon1 should be between 0 and 360
if ($lon1 < $lov - 180.0) { $lon1 += 360.0; }
if ($lon1 > $lov + 180.0) { $lon1 -= 360.0; }
# make sure that lov is within 180 degrees of west
# lov is between 0 and 360 and west is between -180 and 180
if ($lov < ($west+$east)/2 - 180.0) { $lov += 360.0 ; $lon1 += 360.0; };
if ($lov > ($west+$east)/2 + 180.0) { $lov -= 360.0 ; $lon1 -= 360.0; };
print "pdef $nx $ny $t $lat1 $lon1 1 1 $latin1 $latin2 $lov $dx $dy\n";
# if ($lon1 - $west >= 360.0) {
# $lon1 -= 360.0;
# $lov -= 360.0;
# }
# if ($lov < $lon1 - 180.0) { $lov += 360.0; }
# if ($lov > $lon1 + 180.0) { $lov -= 360.0; }
# print "pdef $nx $ny $t $lat1 $lon1 1 1 $latin1 $latin2 $lov $dx $dy\n";
$dx = $dx / (110000.0 * cos($lat1*3.141592654/180.0));
$dy = $dy / 110000.0;
$tmp = int(($east - $west) / $dx + 1);
if (($west == -180) && ($east == 180)) {
$tmp -= 1;
}
print "xdef $tmp linear $west $dx\n";
$tmp = int(($north - $south) / $dy + 1);
print "ydef $tmp linear $south $dy\n";
}
elsif (/ polar stereographic grid: /) {
/ (\S*) pole lat1 (\S*) lon1 (\S*) latD (\S*) lonV (\S*) dx (\S*) m dy (\S*) m/;
$pole = $1;
$lat1 = $2;
$lon1 = $3;
$orient = $5;
$dx = $6;
$dy = $7;
/ \((\S*) x (\S*)\) input (\S*) /;
$nx=$1;
$ny=$2;
$scan=64;
# $scan=-1;
# if ($3 eq 'WE:NS') { $scan = 0; }
# if ($3 eq 'WE:SN') { $scan = 64; }
# if ($scan == -1) { printf STDERR "HELP polar stereographic code needs to be extended!!\n"; }
# probably only works for scan=64
$dpr=3.14159265358979/180.0;
$rearth=6.3712e6;
$h=1;
$proj="nps";
if ($pole eq "South") {
$h=-1;
$orient -= 180.0;
$proj="sps";
}
$hi=1;
$hj=-1;
if (($scan/128 % 2) == 1) {
$hi=-1;
}
if (($scan/64 % 2) == 1) {
$hj=1;
}
if ($dx == 0 || $dy == 0) {
print STDERR "Problem: DX or DY is missing! Cannot make a control file\n";
exit;
}
$dxs=$dx*$hi;
$dys=$dy*$hj;
$de=(1+sin(60*$dpr))*$rearth;
$dr=$de*cos($lat1*$dpr)/(1+$h*sin($lat1*$dpr));
$xp=1-$h*sin(($lon1-$orient)*$dpr)*$dr/$dxs;
$yp=1+cos(($lon1-$orient)*$dpr)*$dr/$dys;
$dx=$h*$dx/1000;
printf "pdef $nx $ny $proj $xp $yp $orient $dx\n";
&domain;
$dx = 1000.0 * abs($dx);
$nx = int(2 * $rearth * 3.14 * ($east-$west)/360.0 / $dx + 1);
$dx = ($east-$west) /($nx-1);
# make dx a faction of 360 so grads can id as a global field
$nx = int(360.0 / $dx) + 1.0;
$dx = 360.0 / $nx;
$nx = int( ($east-$west)/$dx ) + 1;
$dy = abs($dy);
$ny = int($rearth * 3.14 * ($north-$south)/180.0 / $dy + 1);
$dy = ($north-$south)/($ny-1);
$tmp = int(($east - $west) / $dx + 1);
if (($west == -180) && ($east == 180)) {
$tmp -= 1;
}
print "xdef $tmp linear $west $dx\n";
$tmp = int(($north - $south) / $dy + 1);
print "ydef $tmp linear $south $dy\n";
# $dx = 1000.0 * abs($dx);
# $nx = int($rearth * 3.14 / $dx + 1);
# $dx = 360/$nx;
# printf "xdef $nx linear 0 $dx\n";
# if ($proj eq 'sps') {
# $ny=int(($north+90)/$dx)+1;
# printf "ydef $ny linear -90 $dx\n",
# }
# else {
# $ny=int((90-$south)/$dx)+1;
# $l=90-($ny-1)*$dx;
# printf "ydef $ny linear $l $dx\n",
# }
}
elsif ( / thinned gaussian:/) {
/ lat *(\S*) to (\S*)/;
$lat0 = $1;
$lat1 = $2;
/ long (\S*) to (\S*), (\S*) grid pts \(-1 x (\S*)\) /;
$lon0=$1;
$lon1=$2;
$nxny = $3;
$ny = $4;
print "ydef $ny levels\n";
$yrev = 0;
if ($lat0 > $lat1) { $yrev = 1; }
$eps = 3e-14;
$m=int(($ny+1)/2);
$i=1;
while ($i <= $m) {
$z=cos(3.141592654*($i-0.25)/($ny+0.5));
do {
$p1 = 1;
$p2 = 0;
$j = 1;
while ($j <= $ny) {
$p3 = $p2;
$p2 = $p1;
$p1=((2*$j-1)*$z*$p2-($j-1)*$p3)/$j;
$j++;
}
$pp = $ny*($z*$p1-$p2)/($z*$z-1);
$z1 = $z;
$z = $z1 - $p1/$pp;
} until abs($z-$z1) < $eps;
$x[$i] = -atan2($z,sqrt(1-$z*$z))*180/3.141592654;
$x[$ny+1-$i] = -$x[$i];
$i++;
}
$i = 1;
while ($i < $ny) {
printf " %7.3f", $x[$i];
if (($i % 10) == 0) { print "\n"; }
$i++;
}
printf " %7.3f\n", $x[$ny];
$t=$_;
s/\n//g;
/ bdsgrid \S* (.*) min/;
$list=$1;
$i = 1;
$nx = 0;
foreach $t (split(' ',$list)) {
$nlon[$i++] = $t;
if ($nx < $t) { $nx = $t; }
}
if ($lon1 <= $lon0) { $lon1 += 360.0; }
$dx = ($lon1 - $lon0) / ($nx - 1);
print "xdef $nx linear $lon0 $dx\n";
# now to create the pdef file
open (PDEF, ">$pdef") or die "Could not open pdef file for write: $pdef";
binmode PDEF;
print "pdef $nxny 1 file 1 stream binary $caret3$pdef\n";
if ($yrev == 0) {
$offset = 0;
for ($j = 1; $j <= $ny; $j++) {
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - 1.0) * ($nlon[$j] - 1);
$x = floor($x + 0.5) + $offset + $pdef_offset;
print PDEF pack("L", $x);
}
$offset += $nlon[$j];
}
}
else {
$offset = $nxny;
for ($j = $ny; $j >= 1; $j--) {
$offset -= $nlon[$j];
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - 1.0) * ($nlon[$j] - 1);
$x = floor($x + 0.5) + $offset + $pdef_offset;
print PDEF pack("L", $x);
}
}
}
# print weights
$x = pack("f", 1.0);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
# print wind rotation
$x = pack("L", -999);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
close(PDEF);
}
# rotated LatLon grid
elsif (/ rotated lat-lon grid/) {
/ lat (\S*) to (\S*) by (\S*)/;
$lat0 = $1;
$lat1 = $2;
$dlat = $3;
/ lon (\S*) to (\S*) by (\S*)/;
$lon0 = $1;
$lon1 = $2;
$dlon = $3;
if ($lon0 > 180.0) { $lon0 -= 360.0; }
/rotated lat-lon grid:\((\S*) x (\S*)\)/;
$nx = $1;
$ny = $2;
/ south pole lat=(\S*) lon=(\S*) angle of rot=(\S*):/;
$lat_sp = $1;
$lon_sp = $2;
$rot_angle = $3;
print "* Rotated xxxxLatLon grid: South pole lat $lat_sp lon $lon_sp",
" rot angle $rot_angle\n";
if ($dlon < $dlat) { $dlat = $dlon ;}
$dlon = $dlat;
# $dlon = ($lon1-$lon0)/($nx-1);
# $dlat = ($lat1-$lat0)/($ny-1);
$t0 = $lon_sp+180.0;
if ($t0 > 360.0) { $t0 -= 360.0; }
$t1 = -$lat_sp;
$pdef='rotll';
if (/:winds\(grid\)/) { $pdef='rotllr'; }
print "pdef $nx $ny $pdef $t0 $t1 $dlon $dlat $lon0 $lat0\n";
($swlat, $swlon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat1,$lon0);
if ($lat0 > $lat1) {
($swlat, $swlon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat1,$lon0);
($nwlat, $nwlon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat0,$lon0);
($nelat, $nelon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat0,$lon1);
($selat, $selon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat1,$lon1);
$northlat=($nwlat>$nelat)?$nwlat:$nelat;
$southlat=($swlat<$selat)?$swlat:$selat;
$goodny=($northlat-$southlat)/$dlat;
$nlat=$northlat-fmod($northlat,abs($dlat));
print "options yrev\n";
print "ydef $goodny linear $nlat ", abs($dlat), "\n"
}
else {
($swlat, $swlon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat0,$lon0);
($nwlat, $nwlon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat1,$lon0);
($nelat, $nelon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat1,$lon1);
($selat, $selon) = rr2ll($lat_sp,$lon_sp,$rot_angle,$lat0,$lon1);
$northlat=($nwlat>$nelat)?$nwlat:$nelat;
$southlat=($swlat<$selat)?$swlat:$selat;
$goodny=floor(($northlat-$southlat)/$dlat)+1;
$slat=$southlat-fmod($southlat,abs($dlat));
print "ydef $goodny linear $slat ", abs($dlat), "\n"
}
$eastlon=($nelon>$selon)?$nelon:$selon;
$westlon=($swlon<$nwlon)?$swlon:$nwlon;
$goodnx=floor(($eastlon-$westlon)/$dlon)+1;
$wlon=$westlon-fmod($westlon,abs($dlon));
print "xdef $goodnx linear $wlon $dlon\n";
}
elsif (/: Irregular Grid:/) {
/: Irregular Grid:\((\S*) x/;
$nx = $1;
print "xdef $nx linear 0 1\n";
print "ydef 1 linear 0 1\n";
}
elsif (/:grid_template=204:/) {
print "pdef 1200 1684 bilin sequential binary-big pdef_ncep_12.atl\n";
print "xdef 1300 linear -105.04 0.1325\n";
print "ydef 1500 linear -25.2 0.0678\n";
}
else {
print STDERR "*** script needs to be modified ***\n";
print STDERR "*** unknown grid ***\n";
exit 8;
}
# make the tdef statement
&tdef;
# ------------------var-------------------------------------;
# Find the vertical profile
foreach $var (keys(%levels)) {
$var_old=$var;
$vert_lev = '';
if ($profile eq 'mb') {
# $var =~ s/:([\d.e+-]*) mb$/prs/ && do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/prs$/:%s mb/; };
# $var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) mb$/prs/ && do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/prs$/:%s mb/; };
$var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) mb$/:prs/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:prs$/:%s mb/; };
}
if ($profile eq 'iso') {
$var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) K isentropic level$/:iso/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:iso$/:%s K isentropic level/; };
}
if ($profile eq 'bsl') {
$var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) m below sea level$/:bsl/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:bsl$/:%s m below sea level/; };
}
if ($profile eq 'dsl') {
$var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) m below sea level$/:dsl/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:dsl$/:%s m below sea level/; };
}
if ($profile eq 'sig') {
$var =~ s/:([0-9]*\.?[0-9]+([eE][+-]?[0-9]+)?) sigma level$/:sig/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:sig$/:%s sigma level/; };
}
if ($profile eq 'hyp') {
$var =~ s/:([0-9]*) hybrid pressure level$/:hyp/ &&
do { $vert_lev = $1 ; $var_id = $var; $var_id =~ s/:hyp$/:%s hybrid pressure level/; };
}
if ($profile eq 'hply') {
$var =~ /:([0-9]*)-([0-9]*) hybrid pressure layer$/ && do {
if ($1 + 1 == $2) {
$var =~ s/:([0-9]*)-([0-9]*) hybrid pressure layer$/:hply/;
$vert_lev = $1 ; $var_id = $var; $var_id =~ s/:hply$/:%s-%s hybrid pressure layer/;
}
}
}
if ($vert_lev ne '') {
$levels{$var_old} = $var;
if (defined($prof_lev{$var})) {
$prof_lev{$var} = $prof_lev{$var} . ' ' . $vert_lev;
}
else {
$prof_lev{$var} = $vert_lev;
$prof_id{$var} = $var_id;
}
}
}
# find variable with most levels
# sort the levels
$nmax = 1;
$max_levels='1';
foreach $var (keys(%prof_lev)) {
if ($profile_sort eq 'reverse') { $levels = join ' ', sort {$b <=> $a} split(/ /,$prof_lev{$var}); }
else { $levels = join ' ', sort {$a <=> $b} split(/ /,$prof_lev{$var}); }
$prof_lev{$var} = $levels;
$n=($levels =~ s/ //g)+1;
if ($n <= 3) { $prof_lev{$var} = ''; }
else {
if ($n > $nmax) { $max_levels = $prof_lev{$var}; $nmax = $n;}
}
}
print "zdef $nmax levels $max_levels\n";
$max_levels = ' ' . $max_levels . ' ';
# to use a profile, the levels must be a subset of the max_levels
# generate ctl lines for profiles
$iout=0;
foreach $var (keys(%prof_lev)) {
if ($prof_lev{$var} ne '') {
$levels = $prof_lev{$var};
if ($max_levels =~ /$levels/) {
$_ = $';
# WNE??? 9/2013
$n=$nmax - (s/ //g) + 1;
$name = $var;
shorten_name();
if ($short_varname == 0) {
$out[$iout++] = "$name $n 0 \"$prof_id{$var}\" * profile $prof_id{$var}\n";
}
else {
$out[$iout++] = "$name $n 0 \"$prof_id{$var}\"\n";
}
}
else {
$prof_lev{$var} = '';
}
}
}
# generate the ctl line for single levels
foreach $var (sort keys(%levels)) {
$lev = $levels{$var};
$multi_level = 0;
if ($lev ne '0') {
if ($prof_lev{$lev} ne '') { $multi_level = 1; }
}
if ($multi_level == 0) {
$name = $var;
$name =~ s/:[^:]*$/:$lev0{$var}/;
shorten_name();
if ($short_varname == 0) {
$out[$iout++] = "$name 0 0 \"$var\" * $var\n";
}
else {
$out[$iout++] = "$name 0 0 \"$var\"\n";
}
}
}
print "vars $iout\n";
if ($short_varname == 0) {
foreach $line (sort @out) {
print $line;
}
}
else {
$i = 1;
foreach $line (sort @out) {
$_ = $line;
$_ =~ s/^[^ ]*/v$i/;
print $_;
$i++;
}
}
print "endvars\n";
exit 0;
#------------------ jday --------------------
# jday(year,mo,day) return the julian day relative to jan 0
# mo=1..12
#
sub jday {
local($n);
local($nleap);
local($year1);
$n=substr(" 000 031 059 090 120 151 181 212 243 273 304 334",($_[1]-1)*4,4);
$n = $n + $_[2];
$year1 = $_[0] - 1905;
if ($calendar eq '365') {
$n += $year1 * 365;
}
else {
if ($_[1] > 2 && $_[0] % 4 == 0) {
if ($_[0] % 400 == 0 || $_[0] % 100 != 0) {
$n++;
}
}
$nleap = int($year1 / 4);
$n = $n + $nleap + $year1 * 365;
}
$n;
}
#------------------ write tdef statement ------------------
# still not great but better than before
sub tdef {
local($tmp);
local($n);
$n=$ntime;
if ($timestep) { $dt=$timestep }
else {
if ($ntime == 1) {
if ($timestep) { $dt=$timestep }
else { $dt="1mo"; }
}
elsif ( ($minute != $minute1) || ($hour != $hour1)) {
$tmp= (&jday($year1,$mo1,$day1) - &jday($year,$mo,$day)) * 24 * 60 + ($hour1 - $hour)*60
+ $minute1 - $minute;
if (($tmp % 60) == 0) {
$tmp = $tmp / 60;
$dt="${tmp}hr";
$n = (&jday($year_last,$mo_last,$day_last) - &jday($year,$mo,$day)) * 24 + $hour_last - $hour;
$n = int($n / $tmp) + 1;
}
else {
$dt="${tmp}mn";
$n= (&jday($year_last,$mo_last,$day_last) - &jday($year,$mo,$day)) * 24 * 60
+ ($hour_last - $hour)*60 + ($minute_last - $minute);
$n = int($n / $tmp) + 1;
}
}
elsif ($day != $day1) {
$tmp = &jday($year1,$mo1,$day1) - &jday($year,$mo,$day);
$dt="${tmp}dy";
$n=int((&jday($year_last,$mo_last,$day_last) - &jday($year,$mo,$day))/$tmp)+1;
}
elsif ($mo != $mo1) {
# assume that dt < 12 months
$tmp = $year1*12+$mo1 - $year*12-$mo;
$dt="${tmp}mo";
$n = int(($year_last*12+$mo_last - $year*12 - $mo) / $tmp) + 1;
}
else {
$tmp = $year1 - $year;
$dt="${tmp}yr";
$n = int(($year_last - $year) / $tmp) + 1;
}
}
if ($calendar eq "365") {
print "options 365_day_calendar\n";
}
if ($minute == 0) { print "tdef $n linear ${hour}Z$day$month$year $dt\n"; }
else { print "tdef $n linear ${hour}:${minute}Z$day$month$year $dt\n"; }
}
#----------------- find size of domain for regular grids
sub domain {
$_ = `$wgrib2 -d 1 "$file" -domain`;
/:N=(\S*) S=(\S*) W=(\S*) E=(\S*)/;
$north=$1;
$south=$2;
$west=$3;
$east=$4;
}
#----------------- print Gaussian latitudes from lat1..lat2
sub print_range_gau_lats {
local($i);
local($eps);
local($m);
local($z);
local($i);
local($j);
local($p1);
local($p2);
local($p3);
local($z1);
local($x);
$eps = 3e-14;
$m=int(($nGauLats+1)/2);
$i=1;
while ($i <= $m) {
$z=cos(3.141592654*($i-0.25)/($nGauLats+0.5));
do {
$p1 = 1;
$p2 = 0;
$j = 1;
while ($j <= $nGauLats) {
$p3 = $p2;
$p2 = $p1;
$p1=((2*$j-1)*$z*$p2-($j-1)*$p3)/$j;
$j++;
}
$pp = $nGauLats*($z*$p1-$p2)/($z*$z-1);
$z1 = $z;
$z = $z1 - $p1/$pp;
} until abs($z-$z1) < $eps;
$x[$i] = -atan2($z,sqrt(1-$z*$z))*180/3.141592654;
$x[$nGauLats+1-$i] = -$x[$i];
$i++;
}
$i = $j = 1;
while ($i <= $nGauLats) {
if ($x[$i] > $lat0 - 0.002 && $x[$i] < $lat1 + 0.002) {
printf " %7.3f", $x[$i];
if (($j % 14) == 0) { print "\n"; }
$j++;
}
$i++;
}
if ((($j-1) % 14) != 0) { print "\n"; }
}
#----------------- print Gaussian latitudes
sub print_gau_lats {
local($i);
local($eps);
local($m);
local($z);
local($i);
local($j);
local($p1);
local($p2);
local($p3);
local($z1);
local($x);
$eps = 3e-14;
$m=int(($ny+1)/2);
$i=1;
while ($i <= $m) {
$z=cos(3.141592654*($i-0.25)/($ny+0.5));
do {
$p1 = 1;
$p2 = 0;
$j = 1;
while ($j <= $ny) {
$p3 = $p2;
$p2 = $p1;
$p1=((2*$j-1)*$z*$p2-($j-1)*$p3)/$j;
$j++;
}
$pp = $ny*($z*$p1-$p2)/($z*$z-1);
$z1 = $z;
$z = $z1 - $p1/$pp;
} until abs($z-$z1) < $eps;
$x[$i] = -atan2($z,sqrt(1-$z*$z))*180/3.141592654;
$x[$ny+1-$i] = -$x[$i];
$i++;
}
$i = 1;
while ($i < $ny) {
printf " %7.3f", $x[$i];
if (($i % 14) == 0) { print "\n"; }
$i++;
}
printf " %7.3f\n", $x[$ny];
}
sub gau_pdf_near {
#
# writes pdef file for thinned gaussian grid
#
# nx, ny = dimension of output grid
#
print "pdef $npnts 1 general 1 stream binary $caret3$pdef\n";
if ($scan eq 'WE:NS') {
$offset = $npnts;
for ($j = $ny; $j >= 1; $j--) {
$offset -= $nlon[$j];
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
$x = floor($x + 0.5);
if ($x >= $nlon[$j]) { $x -= $nlon[$j] };
$x += $offset + $pdef_offset;
print PDEF pack("L", $x);
}
}
}
else {
$offset = 0;
for ($j = 1; $j <= $ny; $j++) {
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
$x = floor($x + 0.5);
if ($x >= $nlon[$j]) { $x -= $nlon[$j]; }
$x += $offset + $pdef_offset;
print PDEF pack("L", $x);
}
$offset += $nlon[$j];
}
}
# print weights
$x = pack("f", 1.0);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
# print wind rotation
$x = pack("L", -999);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
}
sub gau_pdf_linear {
#
# not working
# linear interpolation along latitude
my @n1 = (); # left locations
my @n2 = (); # right locations
my @w1 = (); # left weights
my @w2 = (); # right weights
print "pdef $npnts 1 general 2 stream binary $caret3$pdef\n";
$offset = 0;
for ($j = 1; $j <= $ny; $j++) {
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
# grid points and weights
$nl = floor($x);
$wl = 1 - ($x - $nl);
$nr = $nl + 1;
$wr = $x - $nl;
if ($nl >= $nlon[$j]) { $nl -= $nlon[$j]; }
if ($nr >= $nlon[$j]) { $nr -= $nlon[$j]; }
$nr += $offset + $pdef_offset;
$nl += $offset + $pdef_offset;
push @n1, $nl;
push @w1, $wl;
push @n2, $nr;
push @w2, $wr;
}
$offset += $nlon[$j];
}
# write locations and weights to PDEF file
$rem =($nx*$ny) % 4;
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("L", $n1[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("LLLL", $n1[$i], $n1[$i+1], $n1[$i+2], $n1[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("f", $w1[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("ffff", $w1[$i], $w1[$i+1], $w1[$i+2], $w1[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("L", $n2[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("LLLL", $n2[$i], $n2[$i+1], $n2[$i+2], $n2[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("f", $w2[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("ffff", $w2[$i], $w2[$i+1], $w2[$i+2], $w2[$i+3]);
}
# wind rotation
$x = pack("L", -999);
for ($i = 0; $i < $rem; $i++) {
print PDEF $x;
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF $x,$x,$x,$x;
}
}
sub rr2ll{
my $polelat=shift;
my $polelon=shift;
my $rotation=shift;
my $rlat=shift;
my $rlon=shift;
my $a = deg2rad(90.0+$polelat);
my $b = deg2rad($polelon);
my $r = deg2rad($rotation);
my $pr = deg2rad($rlat);
my $gr = -deg2rad($rlon);
my $pm = asin(cos($pr)*cos($gr));
my $gm = atan2(cos($pr)*sin($gr),-sin($pr));
my $lat = rad2deg(asin(sin($a)*sin($pm)-cos($a)*cos($pm)*cos($gm-$r)));
my $lon = -rad2deg((-$b+atan2(cos($pm)*sin($gm-$r),
sin($a)*cos($pm)*cos($gm-$r)+cos($a)*sin($pm))));
return ($lat, $lon);
}
#
# reads inventory
# picks out variables and levels
#
sub scan_ListA {
my ($ens, $rec_no, $rec_loc, $npts, $time, $var, $lev00);
open (FileDate, "<$ListA") or die "Could not open wgrib2 inventory: $ListA";
while (<FileDate>) {
chomp;
($rec_no,$rec_loc,$lev00, $npts, $time, $var) = split(/:/,$_,6);
$npts =~ s/npts=//;
$time =~ s/.*=//;
$dates{$time} = 0;
$levels{$var} = 0;
$lev0{$var} = $lev00;
}
close (FileDate);
}
sub gdate{
my ($year, $mo, $day, $hour, $month, $minute);
$_ = $_[0];
$year = substr($_,0,4);
$mo = substr($_,4,2);
$day = substr($_,6,2);
$hour = substr($_,8,2);
$minute = substr($_,10,2);
if ($mo < 0 || $mo > 12) {
print "illegal date code $time\n";
exit 8;
}
$month=substr("janfebmaraprmayjunjulaugsepoctnovdec",$mo*3-3,3);
if ($minute != 0) { return " ${hour}:${minute}Z$day$month$year"; }
return " ${hour}Z$day$month$year";
}
sub pdef_near {
# now to create the pdef file
# regional means grid points are on left and right boundaries
#
# $regional_thinned_grid = 1
#
print "pdef $npnts 1 file 1 stream binary $caret3$pdef\n";
$t=0;
if ($regional_thinned_grid == 1) { $t = 1; }
if ($scan eq 'WE:SN') {
$offset = 0;
for ($j = 1; $j <= $ny; $j++) {
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
$x = floor($x + 0.5);
if ($x >= $nlon[$j]) { $x -= $nlon[$j]; }
$x += $offset + $pdef_offset;
print PDEF pack("L", $x);
}
$offset += $nlon[$j];
}
}
else {
$offset = $npnts;
for ($j = $ny; $j >= 1; $j--) {
$offset -= $nlon[$j];
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
$x = floor($x + 0.5);
if ($x >= $nlon[$j]) { $x -= $nlon[$j] };
$x += $offset + $pdef_offset;
print PDEF pack("L", $x);
}
}
}
# print weights
$x = pack("f", 1.0);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
# print wind rotation
$x = pack("L", -999);
for ($i = 0; $i < $nx*$ny; $i++) {
print PDEF $x;
}
}
sub pdef_linear {
# linear interpolation along latitude
my @n1 = (); # left locations
my @n2 = (); # right locations
my @w1 = (); # left weights
my @w2 = (); # right weights
print "pdef $npnts 1 file 2 stream binary $caret3$pdef\n";
$t=0;
if ($regional_thinned_grid == 1) { $t = 1; }
if ($scan eq 'WE:SN') {
$offset = 0;
for ($j = 1; $j <= $ny; $j++) {
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
# grid points and weights
$nl = floor($x);
$wl = 1 - ($x - $nl);
$nr = $nl + 1;
$wr = $x - $nl;
if ($nl >= $nlon[$j]) { $nl -= $nlon[$j]; }
if ($nr >= $nlon[$j]) { $nr -= $nlon[$j]; }
$nr += $offset + $pdef_offset;
$nl += $offset + $pdef_offset;
push @n1, $nl;
push @w1, $wl;
push @n2, $nr;
push @w2, $wr;
}
$offset += $nlon[$j];
}
}
else {
$offset = $npnts;
for ($j = $ny; $j >= 1; $j--) {
$offset -= $nlon[$j];
for ($i = 0; $i < $nx; $i++) {
$x = $i / ($nx - $t) * ($nlon[$j] - $t);
# grid points and weights
$nl = floor($x);
$wl = 1.0 - ($x - $nl);
$nr = $nl + 1;
$wr = 1.0 - $wl;
if ($nl >= $nlon[$j]) { $nl -= $nlon[$j]; }
if ($nr >= $nlon[$j]) { $nr -= $nlon[$j]; }
$nr+= $offset + $pdef_offset;
$nl+= $offset + $pdef_offset;
push @n1, $nl;
push @w1, $wl;
push @n2, $nr;
push @w2, $wr;
}
}
}
# write locations and weights to PDEF file
$rem =($nx*$ny) % 4;
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("L", $n1[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("LLLL", $n1[$i], $n1[$i+1], $n1[$i+2], $n1[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("f", $w1[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("ffff", $w1[$i], $w1[$i+1], $w1[$i+2], $w1[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("L", $n2[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("LLLL", $n2[$i], $n2[$i+1], $n2[$i+2], $n2[$i+3]);
}
for ($i = 0; $i < $rem; $i++) {
print PDEF pack("f", $w2[$i]);
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF pack("ffff", $w2[$i], $w2[$i+1], $w2[$i+2], $w2[$i+3]);
}
# wind rotation
$x = pack("L", -999);
for ($i = 0; $i < $rem; $i++) {
print PDEF $x;
}
for ($i = $rem; $i < $nx*$ny; $i += 4) {
print PDEF $x,$x,$x,$x;
}
}
# this routine makes shorter names
sub shorten_name {
$_ = $name;
s/,missing=\d*([.:])/$1/;
s/\.std_dev([.:])/sDEV$1/;
s/\.ENS=\+(\d*)([.:])/ENS$1$2/;
s/\.ENS=-(\d*)([.:])/ENSm$1$2/;
s/\.Confidence_Indicator/qc/;
s/\.spatial_max/spmax/;
s/\.spatial_ave/spave/;
s/\.spatial_min/spmin/;
s/\.(\d*)[%]_level/$1PC/;
tr/ :-//d;
tr/\./d/;
if (/^\d/) { $_ = "no" . $_; }
$name = $_;
return;
# old code
for ($i = 0; $i <= $#from; $i++) {
s/$from[$i]/$to[$i]/;
}
$name = $_;
}