#!/usr/bin/perl -w
#   makes a GrADS control file for grib files
#
#   requires wgrib and Perl5
#
#   usage: grib2ctl [options] [grib file] [optional index file] >[control file]
#
#   note: this script does not make the index file .. you have to run gribmap
#
#   Analyses: (using initial time)
#
#      $ grib2ctl.pl example.grib >example.ctl
#      $ gribmap -i example.ctl -0
#
#   Forecasts: (using verifiation time)
#
#      $ grib2ctl.pl -verf example.grib >example.ctl
#      $ gribmap -i example.ctl
#
#   bugs:
#         many
#	  will fail under number of situations
#         finite number of NCEP grids are supported
#
# requires wgrib v1.8.0.12 or higher
# 
# documentation:  http://www.cpc.ncep.noaa.gov/products/wesley/grib2ctl.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
#
#
# added -analysis option
#             December 2007, Davide Sacchetti davide.sacchetti.AT.arpal.org
#
# rotated lat-lon grid:  Graziano Giuliani
#
$version="0.9.15";
use POSIX;
use Math::Trig qw(deg2rad rad2deg);

# ***** if wgrib is not on path, add it here
# $wgrib='/u/wx51we/home/bin/wgrib';
$wgrib='wgrib';

# **** directory of interpolation files
$pdef_dir='/usr/local/lib/grads';
#$pdef_dir=$ENV{"HOME"}/grads;

$wflag="";
$file="";
$index="";
$pdef="";
$prs="prs";
$suffix="";
$model="GFS";
$calendar="";
$lc="";
$global_map="";
$timestep="";
$rr="";
$nearest_neighbor="";
$kludge="";
$no235z="";
$template="";
$ec_gbl="";

$pdef_offset = 0;

$raw="no";
$pdef_nearest=1;
$analysis="";

# vertical levels allowed, index = grib level number
# 100 = mb, 103=m above msl 105-m above ground 107=sigma 109=hybrid level 111=cm down 113=K
# 115=mb above ground 125=cm above ground 160=m below sea level
$allow_profile[20] = '+';
$allow_profile[100] = '-';
$allow_profile[103] = '+';
$allow_profile[105] = '+';
$allow_profile[107] = '-';
$allow_profile[109] = '+';
$allow_profile[111] = '-';
$allow_profile[113] = '+';
$allow_profile[115] = '+';
$allow_profile[125] = '-';
$allow_profile[160] = '-';
$allow_profile[235] = '-';

for ($i = 0; $i <= $#ARGV; $i++) {
   $_ = $ARGV[$i];
   SWITCH: {
      /^-verf/ && do { $wflag="$wflag -verf" ; last SWITCH; };
      /^-ncep_opn/ && do { $wflag="$wflag -ncep_opn" ; last SWITCH; };
      /^-ncep_rean/ && do { $wflag="$wflag -ncep_rean" ; last SWITCH; };
      /^-no_prs/ && do { $prs="" ; last SWITCH; };
      /^-no_suffix/ && do { $suffix="no" ; last SWITCH; };
      /^-rev_z/ && do { last SWITCH; };
      /^-365/ && do { $calendar="365"; last SWITCH; };
      /^-ts(\d+\w+)/ && do { $timestep=$1; last SWITCH; };
      /^-lc$/ && do { $lc="on"; last SWITCH; };
      /^-rr_nn$/ && do { $rr="on"; $nearest_neighbor="on"; $model="ETA"; last SWITCH; };
      /^-rr$/ && do { $rr="on"; $model="ETA"; last SWITCH; };
      /^-eta$/ && do { $model="ETA"; last SWITCH; };
      /^-mrf$/ && do { $model="GFS"; last SWITCH; };
      /^-gfs$/ && do { $model="GFS"; last SWITCH; };
      /^-ruc$/ && do { $model="RUC"; last SWITCH; };
      /^-global$/ && do { $global_map = "on"; last SWITCH; };
      /^-noah$/ && do { $wflag="$wflag -ncep_opn"; last SWITCH; };
      /^-osu/ && do { last SWITCH; };
      /^-ec_gbl/ && do { $ec_gbl = 1; last SWITCH; };
      /^-kludge/ && do { $kludge="on"; last SWITCH; };
      /^-no_235z/ && do { $no235z="on"; last SWITCH; };

      /^-iso_profile/ && do { undef @allow_profile; $allow_profile[235]='-'; last SWITCH };
      /^-prs_profile/ && do { undef @allow_profile; $allow_profile[100]='-'; last SWITCH };
      /^-m_profile/ && do { undef @allow_profile; $allow_profile[103]='+'; 
                                  $allow_profile[105]='+' ; last SWITCH };
      /^-no_profile/ && do { undef @allow_profile; last SWITCH };

      /^-raw/ && do { $raw="yes" ; last SWITCH };
      /^-pdef_linear/ && do { $pdef_nearest=0 ; last SWITCH };
      /^-analysis/ && do { $analysis=$ARGV[++$i]; last SWITCH; };
      /^-/ && do { print STDERR "unknown option: $_\n"; exit 8; };
      if ($file eq "") {
         $file = $_;
      }
      elsif ($index eq "") {
         $index = $_;
      }
      else {
         $pdef = $_;
     }
   }
}

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 control file for grib files\n";
   print STDERR " usage: $0 [options] [grib_file] [optional index file] [optional pdef file] >[ctl file]\n";
   print STDERR " -ncep_opn            .. use NCEP opn grib table for T62 NCEP fields\n";
   print STDERR " -ncep_rean           .. use NCEP reanalysis grib table for T62 NCEP fields\n";
   print STDERR " -verf                .. use forecast verification times\n";
   print STDERR " -no_prs              .. no prs suffix on variable name\n";
   print STDERR " -no_suffix           .. no suffix on variable name\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 " -iso_profile         .. set z coordinate to ocean isotherms\n";
   print STDERR " -prs_profile         .. set z coordinate to pressure (mb)\n";
   print STDERR " -m_profile           .. set z coordinate to meters above ground/msl\n";
   print STDERR " -no_profile          .. no z coordinates\n";
   print STDERR " -raw                 .. raw grid\n";
   print STDERR " -pdef_linear         .. use linear interpolation for thinned Gaussian grids\n";
   print STDERR " -iso_profile         .. make profile using subsurface isoterms\n";
   print STDERR " -analysis YYYYMMDDHH .. set initial time (e.g. -analysis 1999123100)\n";
   print STDERR "                         initial time must have correct time spacing\n";
   print STDERR "\n";
   print STDERR "Note 1: the index file will be generated by the gribmap program, default: grib_file.idx\n";
   print STDERR "Note 2: the pdef file is only generated for thinned lat-lon/Gaussian grids, default: [grib_file].pdef\n";
   exit 8;
}

$_ = $file;
if (/%y4/ || /%y2/ || /%m2/ || /%m1/ || /%d2/ || /%d1/ || /%h2/ ||
    /%h1/ || /%f2/ || /%f3/) { $template='on'; }

if (-d "c:\\") {
   $ListA="c:\\g$$.tmp";
   $TmpFile="c:\\h$$.tmp";
   unlink ($ListA, $TmpFile);
   $sys="win";
}
else {
   $ListA="/tmp/g$$.tmp";
   $TmpFile="/dev/null";
   unlink $ListA;
   $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 ($sys eq 'win') {
      $gfile =~ s=\\=/=g;
   }
   $gfile =~ s/%y4/\\d{4}/g;
   $gfile =~ s/%y2/\\d{2}/g;
   $gfile =~ s/%m2/\\d{2}/g;
   $gfile =~ s/%m1/\\d{1,2}/g;
   $gfile =~ s/%d2/\\d{2}/g;
   $gfile =~ s/%d1/\\d{1,2}/g;
   $gfile =~ s/%h2/\\d{2}/g;
   $gfile =~ s/%h1/\\d{1,2}/g;
   $gfile =~ s/%h3/\\d{3}/g;
   $gfile =~ s/%f2/\\d{2,8}/g;
   $gfile =~ s/%f3/\\d{3,8}/g;

   $dir=$gfile;
   $dir =~ s=(/*)[^/]*$=$1=;
   $gfile =~ s/^$dir//;

   $head=$gfile;
   $head =~ s=\\d\{.*==;
   $tail=$gfile;
   $tail =~ s=.*\\d\{.*\}==;

   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;
   }

# allfiles has the name of all the files
# need to find times: t0, t1, .. t-las

   for ($i = 0; $i <= $#allfiles; $i++) {
      $_ = $allfiles[$i];
      $_ =~ s/$head//;
      $_ =~ s/$tail//;
#     now $_ has the date code / forecast hour

      if ($i == 0) {
         $min_tval = $_;
         $min_t2val = $_;
         $max_tval = $_;
      }
      elsif ($min_tval eq $min_t2val && $min_tval eq $max_tval) {
          if ($_ > $min_tval) {
              $min_t2val = $_;
              $max_tval = $_;
          }
          else {
             $min_tval = $_;
          }
      }
      else {
         if ($_ > $max_tval) { $max_tval = $_; }
         elsif ($_ < $min_tval) {
            $min_t2val = $min_tval;
            $min_tval = $_;
         }
         elsif ($_ < $min_t2val && $_ > $min_tval) {
            $min_t2val = $_;
         }
      }
   }
   $file="$dir$head$min_tval$tail";
   if ($sys eq 'win') {
      $file =~ s=/=\\=g;
   }

#  make inventory of first two files and last file
#  need to get dt and last date

   system "$wgrib $wflag -v \"$dir$head$min_tval$tail\" >$ListA";
   if ($#allfiles >= 1) {
      system "$wgrib $wflag -v \"$dir$head$min_t2val$tail\" >>$ListA";
   }
   if ($#allfiles >= 2) {
      system "$wgrib $wflag -v \"$dir$head$max_tval$tail\" >>$ListA";
   }
}
else {
   system "$wgrib $wflag -v \"$file\" >$ListA";
}

if ( ! -s $ListA ) {
    print STDERR "Big problem:  Possible causes:\n";
    print STDERR "    $file does not exist (bad name? wrong directory?)\n";
    print STDERR "    $file is not a grib file\n";
    print STDERR "    cannot write a temporary file $ListA (permissions? full disk?)\n";
    print STDERR "    wgrib was not found ($wgrib not on path, wrong directory?)\n";
    print STDERR "    $wgrib was compiled for a different machine or operating system\n";
    unlink $ListA;
    exit 8;
}

# make table of dates and variables

open (FileDate, "<$ListA");
while (defined($_ = <FileDate>)) {

   # date table
   $_ =~ s/^.*D=//;
   $d=substr($_, 0, 10);
   $dates{$d}="";

   # variable/level list
   @Fld = split(':', $_, 99);
   $kpds=substr($Fld[3],5);
   ($kpds5,$kpds6,$kpds7) = split(/,/,$kpds);
   $varname = "$Fld[1]:$kpds6";
   if (defined $flevels{$varname}) {
      if (!($flevels{$varname} =~ / $kpds7 /)) {
         $flevels{$varname} .= "$kpds7 ";
      }
   }
   else {
      $flevels{$varname} = " $kpds7 ";
      $fcomments{$varname} = "$kpds5:$Fld[$#Fld]";
   }
}
close (FileDate);
if ($sys eq "win") {
   unlink $TmpFile;
}
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);

if ("$analysis" ne "") {
   $yeara = substr($analysis,0,4);
   $moa = substr($analysis,4,2);
   $daya = substr($analysis,6,2);
   $houra = substr($analysis,8,2);
}
else {
   $yeara = $year;
   $moa = $mo;
   $daya = $day;
   $houra = $hour;
}

if ($mo < 0 || $mo > 12) {
   print "illegal date code $time\n";
   exit 8;
}

if ($ntime == 1 && "$analysis" ne "") {
   if ($analysis <= $time) {
      $year1 = $year;
      $mo1 = $mo;
      $day1 = $day;
      $hour1 = $hour;

      $year_last = $year;
      $mo_last = $mo;
      $day_last = $day;
      $hour_last = $hour;

      $year = $yeara;
      $mo = $moa;
      $day = $daya;
      $hour = $houra;
   } else {
      $year1 = $yeara;
      $mo1 = $moa;
      $day1 = $daya;
      $hour1 = $houra;

      $year_last = $yeara;
      $mo_last = $moa;
      $day_last = $daya;
      $hour_last = $houra;

      $yeara = $year;
      $moa = $mo;
      $daya = $day;
      $houra = $hour;
   }
   $ntime = 2;
} 
elsif ($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);

    $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);
}

# ---------------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 grib2ctl v$version\n";
print "* command line options: @ARGV\n";

# ------------------- grid -----------------------
$griddef = `$wgrib $wflag -V "$file" -d 1 -o $TmpFile`;
$_=$griddef;

/ winds\((\S*)\)/;
$winds=$1;
/ center (\S*) /;
$center=$1;
/ subcenter (\S*) /;
$subcenter=$1;
/ grid=(\S*) /;
$grid=$1;
$LCC='lcc';
if ( $winds eq 'grid') { $LCC='lccr'; }
if ($center == 7 && $subcenter == 15) { $LCC='lcc'; }


# ------------ test for bad ecmwf reduced gaussian grids ------------
if ($center == 98) {
   if (/  thinned gaussian:/) {
      / long (\S*) to (\S*), (\S*) grid pts   \(-1 x (\S*)\) /;
      $lon0=$1;
      $lon1=$2;
      $nxny = $3;
      $ny = $4;

      if ($lon0 != $lon1) {
         s/\n//g;
         / bdsgrid \S* (.*) min/;
         $list=$1;

         ($maxlons) = split(' ' , $list);
         $minlons = $maxlons;
         foreach $t (split(' ',$list)) {
            if ($t > $maxlons) { $maxlons = $t; }
            if ($t < $minlons) { $minlons = $t; }
         }
         if ($lon1 < $lon0) { $lon1 += 360.0; }

         if (abs(360 - 360/$maxlons - $lon1) < 0.001) {
             print "* -ec_gbl set\n";
             $ec_gbl='1';
         }
         if (abs(360 - 360/$minlons - $lon1) < 0.001) {
             print "* -ec_gbl set.\n";
             $ec_gbl='1';
         }
      }
   }
}


print "dtype grib $grid\n";
if ($template eq "on") {
   print "options template\n";
}

if ($raw eq 'yes') {
   / nx (\S*) ny (\S*) /;
   $nx=$1;
   $ny=$2;
   print "xdef $nx linear 0 0.1\n";
   print "ydef $ny linear 0 0.1\n";
}
elsif (/ thinned latlon: /) {
   / lat *(\S*) to (\S*) by (\S*) /;
   $lat0 = $1;
   $lat1 = $2;
   $dy = $3;
   / long (\S*) to (\S*), (\S*) grid pts   \(-1 x (\S*)\) /;
   $lon0=$1;
   $lon1=$2;
   $nxny = $3;
   $ny = $4;
   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/\n//g;
   / bdsgrid \S* (.*) min/;
   $list=$1;
   $i = 1;
   foreach $t (split(' ',$list)) {
     $nlon[$i++] = $t;
   }
   $nx=$nlon[1];
   if ($nx < $nlon[$ny]) { $nx = $nlon[$ny]; }
   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");
   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);
}
elsif (/  latlon: /) {
   / lat  (\S*) to (\S*) by (\S*) /;
   $lat0=$1;
   $lat1=$2;
   $dlat=$3;

   / long (\S*) to (\S*) by (\S*), \((\S*) x (\S*)\)/;
   $lon0=$1;
   # $lon1=$2;
   $dlon=$3;
   $nx =$4;
   $ny =$5;

   if ($lat0 > $lat1) {
      print "options yrev\n";
      print "ydef $ny linear $lat1 ", abs($dlat), "\n"
   }
   else {
      print "ydef $ny linear $lat0 ", abs($dlat), "\n"
   }
   print "xdef $nx linear $lon0 $dlon\n";
}
elsif ($grid == 5 && $center == 7) {
   print "pdef 53 57 nps 27 49 -105 190.5\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 6 && $center == 7) {
   print "pdef 53 45 nps 27 49 -105 190.5\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 87 && $center == 7) {
   print "pdef 81 62 nps 31.9 112.53 -105 68.513\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 96 && $kludge eq "on" && $center == 7) {
   # quick and dirty pdef for eta 12 km
   # use -kludge option
   print "pdef 606 1067 eta.u -111 50 0.17520661  0.075046904\n";
   print "xdef 1440 linear -200 0.125\n";
   print "ydef 721 linear 0 0.125\n";
}
elsif ($grid == 101 && $center == 7) {
   print "pdef 113 91 nps 58.5 92.5 -105 91.452\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 104 && $center == 7) {
   print "pdef 147 110 nps 75.5 109.5 -105 90.75464\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 105 && $center == 7) {
   print "pdef 83 83 nps 40.5 88.5  -105 90.75464\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 106 && $center == 7) {
   print "pdef 165 117 nps 80 176 -105 45.37732\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 107 && $center == 7) {
   print "pdef 120 92 nps 46 167  -105 45.37732\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 192 && $rr ne "" && $center == 7) {
   # quick and dirty pdef for RR egrid
   # need -rr option to use
   print "pdef 237 387 eta.u -111 50 0.4491525  0.2072539\n";
   print "xdef 360 linear -180 0.5\n";
   print "ydef 181 linear 0 0.5\n";
}
elsif ($grid == 201 && $center == 7) {
   print "pdef 65 65 nps 33 33 -105 381\n";
   print "xdef 180 linear -180 2\n";
   print "ydef 51 linear -10 2\n";
}
elsif ($grid == 202 && $center == 7) {
   print "pdef 65 43 nps 33 45 -105 190.5\n";
   print "xdef 91 linear -200 2\n";
   print "ydef 41 linear 10 2\n";
}
elsif ($grid == 203 && $center == 7) {
   print "pdef  45 39 nps 27 37 -150 190.5\n";
   print "xdef 103 linear -250 2\n";
   print "ydef 33 linear 26 2\n";
}
elsif ($grid == 205 && $center == 7) {
   print "pdef  45 39 nps 27 57 -60 190.5\n";
   print "xdef 50 linear -120 2\n";
   print "ydef 46 linear 0 2\n";
}
elsif ($grid == 207 && $center == 7) {
   print "pdef 49 35 nps 25 51 -150 95.25\n";
   print "xdef 51 linear -200 2\n";
   print "ydef 30 linear 45 1\n";
}
elsif ($grid == 211 && $center == 7) {
   # awips lambert conformal
   print "pdef 93 65 $LCC 12.19 -133.459 1 1 25 25 -95 81270.5 81270.5\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
elsif ($grid == 212 && $center == 7) {
   # awips lambert conformal
   print "pdef 185 129 $LCC 35.0 -95.0 105 49 25 25 -95 40635 40635\n";
   print "xdef 181 linear -140 0.5\n";
   print "ydef  91 linear  15 0.5\n";
}
elsif ($grid == 213 && $center == 7) {
   print "pdef  129 85 nps 65 89 -105 95.25\n";
   print "xdef 170 linear -190 1\n";
   print "ydef 81 linear 10 1\n";
}
elsif ($grid == 215 && $center == 7) {
   # lambert conformal
   print "pdef 369 257 $LCC 12.19 -133.46 1 1 25 25 -95 20318 20318\n";
   print "xdef 289 linear -136 0.25\n";
   print "ydef  157 linear 18 0.25\n";
}
elsif ($grid == 216 && $center == 7) {
   print "pdef 147 110 nps 75.5 109.5 -105 91.452\n";
   print "xdef 181 linear -180 1\n";
   print "ydef 91 linear 0 1\n";
}
elsif ($grid == 221 && $nearest_neighbor eq "on" && $center == 7) {
   # awips lambert conformal nearest neighbor
   print "pdef 96673 1 file 1 sequential binary-little $pdef_dir/grib221nn.pdef\n";
   print "xdef 1111 linear -250 0.333333\n";
   print "ydef 247 linear 8 0.333333\n";
}
elsif ($grid == 221 && $center == 7) {
   # awips lambert conformal
   print "pdef 349 277 $LCC 1 -145.5 1 1 50 50 -107 32463 32463\n";
   print "xdef 1111 linear -250 0.333333\n";
   print "ydef 247 linear 8 0.333333\n";
}
# old grid 240
#elsif ($grid == 240) {
#   # nps usa
#   print "pdef 1160 880 nps 441 1601 255 4.763\n";
#   print "xdef 801 linear -130 0.1\n";
#   print "ydef 401 linear 20 0.1\n";
#}
# new grid 240
elsif ($grid == 240 && $center == 7) {
   # nps usa
   print "pdef 1121 881 nps 401 1601 255 4.7625\n";
   print "xdef 1601 linear -130 0.05\n";
   print "ydef 801 linear 20 0.05\n";
}
elsif ($grid == 241 && $center == 7) {
   print "pdef 386 293 nps 147.315 534.0 -105 14.2875\n";
   print "xdef 161 linear -140 0.5\n";
   print "ydef 81 linear 20 0.5\n";
}
else {
   # unknown grid

   $_ = $griddef;
   GRD: {

      / polar stereo: Lat1 16.125000 Long1 234.983000 Orient -100.0/ && do {
         print "options yrev\n";
         print "pdef 129 86 nps 64 136 -100 60\n";
         print "xdef 720 linear  0 0.5\n";
         print "ydef 148 linear 16 0.5\n";
	 last GRD; };

      / polar stereo: Lat1 -4.860000 Long1 -122.614000 Orient -80.000000/ && do {
         print "pdef 49 51 nps 24 26 -80 381\n";
         print "xdef 144 linear 0 2.5\n";
         print "ydef 45 linear -20 2.5\n";
         last GRD; };

      / Lambert Conf:.* Lov 265.*\(151 x 113\)/s && do {
         print "options yrev\n";
         print "pdef 151 113 $LCC 16.281 233.8622 1 1 25 25 265 40635 40635\n";
         print "xdef 141 linear -130 0.5\n";
         print "ydef 71 linear 20 0.5\n";
         last GRD; };

      # beta: mercator 
      # scan modes .. assumes west to east
      / Mercator: / && do {
	/lat  *(\S*) to (\S*) /;
        $lat1 = $1;
        $lat2 = $2;
	/long (\S*) to (\S*) /;
        $lon1 = $1;
        $lon2 = $2;
	/ nx (\S*) ny (\S*) /;
	$nx = $1;
	$ny = $2;

	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"; }

	$dlon = $lon2 - $lon1;
	if ($dlon < 0) {
	    $dlon = $dlon + 360;
	}
	$dlon = $dlon / ($nx - 1);
	print "xdef $nx linear $lon1 $dlon\n";

	last GRD; };


      # beta: generalized lambert conformal pdef/xdef/ydef
      #  not very good .. needs to calculate the all
      #  vertices for better xdef and ydef
      #  for improvements .. pull out code from lcgrib

      / Lambert Conf: / && do {
	/ Lat1 (\S*) Lon1 (\S*) Lov (\S*)/;
	$lat1 = $1;
	$lon1 = $2;
	$lov = $3;

	/Latin1 (\S*) Latin2 (\S*) /;
	$latin1 = $1;
	$latin2 = $2;

	/Pole \((\S*) x (\S*)\) Dx (\S*) Dy (\S*) /;
	$nx = $1;
	$ny = $2;
	$dx = 1000*$3;
	$dy = 1000*$4;

#       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; }
        if ($lon1 < -180) { $lon1 += 360; $lov += 360; }
        if ($lon1 >= 180) { $lon1 -= 360; $lov -= 360; }

	print "pdef $nx $ny $LCC $lat1 $lon1 1 1 $latin1 $latin2 $lov $dx $dy\n";
        if ($global_map eq "") {
            $dx = $dx / (110000.0 * cos($lat1*3.141592654/180.0));
	    $dy = $dy / 110000.0;
	    if ($lon1 > 180) {
		$lon1 = $lon1 - 360.0;
	    }
	    print "xdef $nx linear $lon1 $dx\n";
	    print "ydef $ny linear $lat1 $dy\n";
	}
	else {
	    print "xdef 360 linear 0 1\n";
	    print "ydef 181 linear -90 1\n";
	}

	last GRD; };

      /  thinned gaussian:/ && do {

        / 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);

#       get gaussian latitudes

        $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; }

	if ($ec_gbl eq "") {
           $dx = ($lon1 - $lon0) / ($nx - 1);
	   $t = 1;
        }
        else {
           $dx = 360.0 / $nx;
	   $t = 0;
        }
        print "xdef $nx linear $lon0 $dx\n";

#  now to create the pdef file

        open (PDEF, ">$pdef");
        binmode PDEF;

        if ($pdef_nearest == 0) {
           gau_pdf_linear ();
        }
        else {
           gau_pdf_near ();
        }

        close(PDEF);
        last GRD; };

      /  gaussian:/ && do {
         / lat  (\S*) to (\S*)/;
         $lat0=$1;
         $lat1=$2;

         / long (\S*) to (\S*) by (\S*), \((\S*) x (\S*)\)/;
         $lon0=$1;
         # $lon1=$2;
         $dlon=$3;
         $nx =$4;
         $ny =$5;
         $dlon = 360 / $nx;

         if ($lat0 > $lat1) {
            print "options yrev\n";
         }
         print "xdef $nx linear $lon0 $dlon\n";
         print "ydef $ny levels\n";

	$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];

         last GRD; };
#     rotated LatLon grid
      / rotated LatLon grid/ && do {
         / LatLon grid  lat (\S*) to (\S*)  lon (\S*) to (\S*)/;
         $lat0 = $1;
         $lat1 = $2;
         $lon0 = $3;
         $lon1 = $4;
         /nxny \S+  \((\S*) x (\S*)\)/;
         $nx = $1;
         $ny = $2;
         / south pole lat (\S*) lon (\S*)  rot angle (\S*)/;
         $lat_sp = $1;
         $lon_sp = $2;
         $rot_angle = $3;
         print "* Rotated LatLon grid: South pole lat $lat_sp lon $lon_sp",
               "  rot angle $rot_angle\n";

         $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 eq '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";
        last GRD; };

#     polar stereographic
     /  polar stereo: / && do {
        / Lat1 (\S*) Long1 (\S*) Orient (\S*)/;
        $lat1=$1;
        $lon1=$2;
        $orient=$3;

        / (\S*) pole \((\S*) x (\S*)\) Dx (\S*) Dy (\S*) scan (\S*)/;
        $pole=$1;
        $nx=$2;
        $ny=$3;
        $dx=$4;
        $dy=$5;
        $scan=$6;

        # 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;
        }
        $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";

        # need to do a better job here
        # need to find lat/lon of end points to right domain
        # for now, just do simple stuff
        $dx=abs($dx)*1000;
        $nx = int(2 * $rearth * 3.14 / $dx + 1);
        $dx = 360/$nx;
        printf "xdef $nx linear 0 $dx\n";
        if ($proj eq 'sps') {
           $ny=int(($lat1+90)/$dx)+1;
           printf "ydef $ny linear -90 $dx\n",
        }
        else {
           $ny=int((90-$lat1)/$dx)+1;
           $l=90-($ny-1)*$dx;
           printf "ydef $ny linear $l $dx\n",
        }
        last GRD; };


      print STDERR "*** script needs to be modified ***\n";
      print STDERR "unknown user-defined grid\n";
   }
}


# make the tdef statement

&tdef;

# ------------------var-------------------------------------;

%tails =(
   '1' => 'sfc',
   '2' => 'clb',
   '3' => 'clt',
   '4' => 'zdg',
   '5' => 'adcl',
   '6' => 'mwl',
   '7' => 'trp',
   '8' => 'toa',
   '9' => 'bos',
   '10' => 'clm',
   '12' => 'lcb',
   '13' => 'lct',
   '14' => 'loc',
   '20' => 'tmpl',
   '22' => 'mcb',
   '23' => 'mct',
   '24' => 'mdc',
   '32' => 'hcb',
   '33' => 'hct',
   '34' => 'hic',
   '100' => 'prs',
   '101' => 'plr',
   '102' => 'msl',
   '103' => 'hml',
   '104' => 'zlr',
   '105' => 'hag',
   '106' => 'hlr',
   '107' => 'sig',
   '108' => 'slr',
   '109' => 'hbl',
   '110' => 'blr',
   '111' => 'dpl',
   '112' => 'dlr',
   '113' => 'tht',
   '114' => 'tlr',
   '114' => 'ceil',
   '116' => 'plg',
   '117' => 'pv',
   '119' => 'eta',
   '121' => 'plr',
   '126' => 'pa',
   '128' => 'slr',
   '141' => 'plr',
   '160' => 'dsl',
   '200' => 'clm',
   '201' => 'oclm',
   '204' => 'htfl',
   '206' => 'gcbl',
   '207' => 'gctl',
   '209' => 'bcb',
   '210' => 'bct',
   '211' => 'bcl',
   '212' => 'lcb',
   '213' => 'lct',
   '214' => 'lcl',
   '220' => 'pbl',
   '222' => 'mcb',
   '223' => 'mct',
   '224' => 'mcl',
   '232' => 'hcb',
   '233' => 'hct',
   '234' => 'hcl',
   '235' => 'tmp',
   '237' => 'bmxl',
   '238' => 'bitl',
   '240' => 'omxl',
   '242' => 'cvb',
   '243' => 'cvt',
   '244' => 'cvl',
   '248' => 'sccb',
   '249' => 'scct',
   '251' => 'dccb',
   '252' => 'dcct',
   );

#  Find the vertical profile

$nlevelmax=0;
$levelsmax=0;
$profile_kpds6=-1;

foreach $fname (sort keys(%flevels)) {
   ($name, $kpds6) = split(/:/, $fname);
   $_=$flevels{$fname};
   $nlev = (tr/ / /) - 1;
   if ($nlev > $nlevelmax) {
      if (defined($allow_profile[$kpds6])) {
         $nlevelmax = $nlev;
         $levelsmax = $flevels{$fname};
         $profile_kpds6 = $kpds6;
      } 
   }
}

if ($nlevelmax <= 1) {
   print "zdef 1 linear 1 1\n";
}
else {
   print "*  z has $nlevelmax levels, for $tails{$profile_kpds6}\n";
   print "zdef $nlevelmax levels\n";
   ($_ = $levelsmax) =~ s/.//;
   chop($_);
   if ($allow_profile[$profile_kpds6] eq '+') {
      print join ' ', sort {$a <=> $b} split(/ /,$_);
   }
   else { 
      print join ' ', sort {$b <=> $a} split(/ /,$_);
   }
   print "\n";
}
if ($nlevelmax == 1) { $profile_kpds6 = -1; }

$tails{'100'} = "$prs";
$nvar = 0;

foreach $fname (sort keys(%flevels)) {
   ($name, $kpds6) = split(/:/, $fname);
   ($kpds5, $comment) = split(/:/, $fcomments{$fname});
   $comment = substr($comment,1);

   $_=$flevels{$fname};
   $nlev = (tr/ / /) - 1;
   $kpds7s = $_;

   # fix names to be grads compatible
   # eliminate dashes, underscores, blanks and put no in front of leading digits

   $_ = $name;
   $_ =~ tr/_\- //d;
   if ($lc) { $_ =~ s/(.*)/\L$1/gi; }
   if ( /^\d/ ) { $_ = "no$_"; }
   $name = $_;
   
   $tail = $suffix eq 'no' ? "" : $tails{$kpds6};
   if (! defined $tail) { $tail="l$kpds6"; }

   # convert profile variable

   if ($profile_kpds6 == $kpds6 && $nlev != 1) {
      $out[$nvar++] = "$name$tail $nlev $kpds5,$kpds6,0 ** (profile) $comment";
   }

   # convert profiles to individual layers
   elsif ($kpds6 == 100) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}${ll}mb  0 $kpds5,$kpds6,$ll ** $ll mb $comment";
      }
   }
   elsif ($kpds6 == 101) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $out[$nvar++] = "${name}${lev1}0_${lev2}0mb  0 $kpds5,$kpds6,$ll ** ${lev1}0-${lev2}0 mb $comment";
      }
   }
   elsif ($kpds6 == 103) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}${ll}m  0 $kpds5,$kpds6,$ll ** $ll m above msl $comment";
      }
   }
   elsif ($kpds6 == 105) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}${ll}m  0 $kpds5,$kpds6,$ll ** $ll m above ground $comment";
      }
   }
   elsif ($kpds6 == 106) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $lev1f=100*$lev1;
         $lev2f=100*$lev2;
         $out[$nvar++] = "${name}${lev2f}_${lev1f}m   0 $kpds5,$kpds6,$ll ** ${lev2f}-${lev1f} m above ground $comment";
      }
   }
   elsif ($kpds6 == 107) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         if ($ll < 10) {$ll="0$ll";}
         if ($ll < 100) {$ll="0$ll";}
         if ($ll < 1000) {$ll="0$ll";}
         $lev1 = $ll;
         $lev1 =~ s/0+$/ /;
         $out[$nvar++] = "${name}sig$lev1  0 $kpds5,$kpds6,$ll ** sigma=.${lev1} $comment";
      }
   }
   elsif ($kpds6 == 108) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $lev1f = $lev1 / 100;
         $lev2f = $lev2 / 100;
         $out[$nvar++] = "${name}sg${lev1}_${lev2}   0 $kpds5,$kpds6,$ll ** sigma=${lev1f}-${lev2f} layer $comment";
      }
   }
   elsif ($kpds6 == 109) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}hlev$ll  0 $kpds5,$kpds6,$ll ** hybrid level $ll $comment";
      }
   }
   elsif ($kpds6 == 110) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $out[$nvar++] = "${name}hg${lev1}_${lev2}   0 $kpds5,$kpds6,$ll ** hybrid=${lev1}-${lev2} layer $comment";
      }
   }
   elsif ($kpds6 == 111) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}_${ll}cm  0 $kpds5,$kpds6,$ll ** $ll cm underground $comment";
      }
   }
   elsif ($kpds6 == 112) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $out[$nvar++] = "${name}${lev1}_${lev2}cm  0 $kpds5,$kpds6,$ll ** $lev1-$lev2 cm underground $comment";
      }
   }
   elsif ($kpds6 == 20) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $flt = $ll / 100.0;
	 # cK = centi-Kelvins
         $out[$nvar++] = "${name}${ll}cK  0 $kpds5,$kpds6,$ll ** ${flt}K level $comment";
      }
   }
   elsif ($kpds6 == 113) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}${ll}K  0 $kpds5,$kpds6,$ll ** ${ll}K level $comment";
      }
   }
   elsif ($kpds6 == 116) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256);
         $lev2 = $ll % 256;
         $out[$nvar++] = "${name}${lev1}_${lev2}mb  0 $kpds5,$kpds6,$ll ** ${lev1}-${lev2} mb above gnd $comment";
      }
   }
   elsif ($kpds6 == 117) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1=$ll;
         if ($lev1 > 32668) { $lev1 = 32768 - $lev1; }
         $lev2 = $lev1;
         if ($lev1 < 0) { $lev1 = -$lev1; $lev1 = "neg$lev1"; }
         $out[$nvar++] = "${name}pv${lev1}  0 $kpds5,$kpds6,$ll ** pot vorticity = ${lev2} units level $comment";
      }
   }
   elsif ($kpds6 == 119) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $lev1 = $ll / 10000.0;
         $lev2 = $lev1;
         $lev2 =~ s/\./p/;
         $out[$nvar++] = "${name}eta${lev2}  0 $kpds5,$kpds6,$ll ** Eta=${lev1} level $comment";
      }
   }
   elsif ($kpds6 == 160) {
      foreach $ll (sort {$b <=> $a} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}m${ll}  0 $kpds5,$kpds6,$ll ** ${ll} m below sea level $comment";
      }
   }




   elsif ($kpds6 == 235) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $val = $ll;
         if ($val > 32768) { $val = 32768 - $val; }
         $iso=$val/10;
         $isovar=$iso;
         $isovar =~ s/-/neg/;
         $isovar =~ s/[.]/p/;
         $out[$nvar++] = "${name}t${isovar}c  0 $kpds5,$kpds6,$ll ** ${iso}C isotherm $comment";
      }
   } 
   elsif ($kpds6 == 236) {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $lev1 = int($ll / 256)*10;
         $lev2 = ($ll % 256)*10;
         $out[$nvar++] = "${name}${lev1}_${lev2}m  0 $kpds5,$kpds6,$ll ** $lev1-$lev2 m under water $comment";
      }
   }
   elsif ($kpds6 == 1) {
      $kpds7s =~ s/^ //;
      $out[$nvar++] = "${name}sfc  0 $kpds5,$kpds6,$kpds7s ** surface $comment";
   }
   elsif ($kpds6 == 2) {
      $out[$nvar++] = "${name}clb  0 $kpds5,$kpds6,0 ** cloud base $comment";
   }
   elsif ($kpds6 == 3) {
      $out[$nvar++] = "${name}clt  0 $kpds5,$kpds6,0 ** cloud top $comment";
   }
   elsif ($kpds6 == 4) {
      $out[$nvar++] = "${name}0deg  0 $kpds5,$kpds6,0 ** 0C isotherm level $comment";
   }
   elsif ($kpds6 == 5) {
      $out[$nvar++] = "${name}adcl  0 $kpds5,$kpds6,0 ** adiabatic lifting condensation level $comment";
   }
   elsif ($kpds6 == 6) {
      $out[$nvar++] = "${name}mwl  0 $kpds5,$kpds6,0 ** max wind level $comment";
   }
   elsif ($kpds6 == 7) {
      $out[$nvar++] = "${name}trp  0 $kpds5,$kpds6,0 ** tropopause $comment";
   }
   elsif ($kpds6 == 8) {
      $out[$nvar++] = "${name}toa  0 $kpds5,$kpds6,0 ** top of atmos $comment";
   }
   elsif ($kpds6 == 102) {
      $out[$nvar++] = "${name}msl  0 $kpds5,$kpds6,0 ** mean-sea level $comment";
   }
   elsif ($kpds6 == 200) {
      $out[$nvar++] = "${name}clm  0 $kpds5,$kpds6,0 ** atmos column $comment";
   }
   elsif ($kpds6 == 201) {
      $out[$nvar++] = "${name}oclm  0 $kpds5,$kpds6,0 ** ocean column $comment";
   }
   elsif ($kpds6 == 204) {
      $out[$nvar++] = "${name}htfl  0 $kpds5,$kpds6,0 ** highest trop freezing level $comment";             
   }
   elsif ($kpds6 == 206) {
      $out[$nvar++] = "${name}gcbl  0 $kpds5,$kpds6,0 ** grid-scale cloud bottom level $comment";             
   }
   elsif ($kpds6 == 207) {
      $out[$nvar++] = "${name}gctl  0 $kpds5,$kpds6,0 ** grid-scale cloud top level $comment";             
   }
   elsif ($kpds6 == 209) {
      $out[$nvar++] = "${name}bcb  0 $kpds5,$kpds6,0 ** boundary cld base $comment";
   }
   elsif ($kpds6 == 210) {
      $out[$nvar++] = "${name}bct  0 $kpds5,$kpds6,0 ** boundary cld top $comment";
   }
   elsif ($kpds6 == 211) {
      $out[$nvar++] = "${name}bcl  0 $kpds5,$kpds6,0 ** boundary cld layer $comment";
   }
   elsif ($kpds6 == 212) {
      $out[$nvar++] = "${name}lcb  0 $kpds5,$kpds6,0 ** low cloud base $comment";
   }
   elsif ($kpds6 == 213) {
      $out[$nvar++] = "${name}lct  0 $kpds5,$kpds6,0 ** low cloud top $comment";
   }
   elsif ($kpds6 == 214) {
      $out[$nvar++] = "${name}lcl  0 $kpds5,$kpds6,0 ** low cloud level $comment";
   }
   elsif ($kpds6 == 222) {
      $out[$nvar++] = "${name}mcb  0 $kpds5,$kpds6,0 ** mid-cloud base $comment";
   }
   elsif ($kpds6 == 223) {
      $out[$nvar++] = "${name}mct  0 $kpds5,$kpds6,0 ** mid-cloud top $comment";
   }
   elsif ($kpds6 == 224) {
      $out[$nvar++] = "${name}mcl  0 $kpds5,$kpds6,0 ** mid-cloud level $comment";
   }
   elsif ($kpds6 == 232) {
      $out[$nvar++] = "${name}hcb  0 $kpds5,$kpds6,0 ** high cloud base $comment";
   }
   elsif ($kpds6 == 233) {
      $out[$nvar++] = "${name}hct  0 $kpds5,$kpds6,0 ** high cloud top $comment";
   }
   elsif ($kpds6 == 234) {
      $out[$nvar++] = "${name}hcl  0 $kpds5,$kpds6,0 ** high cloud level $comment";
   }
   elsif ($kpds6 == 237) {
      $out[$nvar++] = "${name}${tail}  0 $kpds5,$kpds6,0 ** bottom ocn mixed layer $comment";
   }
   elsif ($kpds6 == 238) {
      $out[$nvar++] = "${name}${tail}  0 $kpds5,$kpds6,0 ** bottom ocn isothermal layer $comment";
   }
   elsif ($kpds6 == 242) {
      $out[$nvar++] = "${name}cvb  0 $kpds5,$kpds6,0 ** convective cld base $comment";
   }
   elsif ($kpds6 == 243) {
      $out[$nvar++] = "${name}cvt  0 $kpds5,$kpds6,0 ** convective cld top $comment";
   }
   elsif ($kpds6 == 244) {
      $out[$nvar++] = "${name}cvl  0 $kpds5,$kpds6,0 ** convective cld layer $comment";
   }
   elsif ($kpds6 == 246) {
      $out[$nvar++] = "${name}mept  0 $kpds5,$kpds6,0 ** max e-pot-t level $comment";
   }
   elsif ($kpds6 == 248) {
      $out[$nvar++] = "${name}sccb  0 $kpds5,$kpds6,0 ** shallow convective cloud base $comment";
   }
   elsif ($kpds6 == 249) {
      $out[$nvar++] = "${name}scct  0 $kpds5,$kpds6,0 ** shallow convective cloud top $comment";
   }
   elsif ($kpds6 == 251) {
      $out[$nvar++] = "${name}dccb  0 $kpds5,$kpds6,0 ** deep convective cloud base $comment";
   }
   elsif ($kpds6 == 252) {
      $out[$nvar++] = "${name}dcct  0 $kpds5,$kpds6,0 ** deep convective cloud top $comment";
   }
   elsif ($nlev == 1) {
      # unknown single level
      $kpds7s =~ s/^ //;
      $out[$nvar++] = "${name}$tail  0 $kpds5,$kpds6,$kpds7s ** unknown level $comment";
   }
   else {
      foreach $ll (sort {$a <=> $b} split(' ',$kpds7s)) {
         $out[$nvar++] = "${name}$tail$ll  0 $kpds5,$kpds6,$ll **level=$ll $comment";
      }
   }
}

print "vars $nvar\n";
for ($i = 0; $i < $nvar; $i++) {
   print $out[$i];
}
print "ENDVARS\n";

#if ($sys eq "win") {
#   unlink $TmpFile;
#}
# unlink $ListA;
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;

   $month=substr("janfebmaraprmayjunjulaugsepoctnovdec",$moa*3-3,3);

   if ($timestep) { $dt=$timestep  }
   else { 
      if ($ntime == 1) {
         if ($timestep) { $dt=$timestep }
         else { $dt="1mo"; }
      }
      elsif ($hour != $hour1) {
         $tmp= (&jday($year1,$mo1,$day1) - &jday($year,$mo,$day)) * 24 + $hour1 - $hour;
         $dt="${tmp}hr";
         $n = (&jday($year_last,$mo_last,$day_last) - &jday($yeara,$moa,$daya)) * 24 + $hour_last - $houra;
         $n = int($n / $tmp) + 1;
      }
      elsif ($day != $day1) {
         # assume that dt < 365 days
         $tmp = &jday($year1,$mo1,$day1) - &jday($year,$mo,$day);
         $dt="${tmp}dy";
         $n=int((&jday($year_last,$mo_last,$day_last) - &jday($yeara,$moa,$daya))/$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 - $yeara*12 - $moa) / $tmp) + 1;
      }
      else {
         $tmp = $year1 - $year;
         $dt="${tmp}yr";
         $n = int(($year_last - $yeara) / $tmp) + 1;
      }
   }
   if ($calendar eq "365") {
       print "options 365_day_calendar\n";
   }
   print "tdef $n linear ${houra}Z$daya$month$yeara $dt\n";
}

sub gau_pdf_near {
   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 - $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 = $nxny;
           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 gau_pdf_linear {

#  linear interpolation along latitude

   my @n1 = ();  # left locations
   my @n2 = ();  # right locations
   my @w1 = ();  # left weights
   my @w2 = ();  # right weights

   print "pdef $nxny 1 file 2 stream binary $caret3$pdef\n";
   print "* nx $nx ny $ny nxny $nxny\n";

   if ($yrev == 0) {
           $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 = $nxny;
      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;
   }
}

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);
}