[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: New Ticket - [netCDFDecoders !YUU-262958]: ua2nc bug
- Subject: Re: New Ticket - [netCDFDecoders !YUU-262958]: ua2nc bug
- Date: Wed, 13 Sep 2006 09:52:57 -0600 (MDT)
Ryan,
the problem was the bulletin came in on same day so it was filed into the
0th hour for that day. i added some code to file the bulletin into the
next day 0th hour. the ua2nc code is attached, could you test it out
for any problems.
robb...
On Thu, 7 Sep 2006, Ryan Torn wrote:
> New Ticket: ua2nc bug
>
> I think I have found a fairly major bug in the decoder ua2nc. The
> report time on Chinese rawinsonde observations at 00UTC are typically
> marked as 23 UTC on the previous day. For example, the 00 UTC
> rawinsonde observation for the 7th has a report time of 23 UTC on the
> 6th. When this report is fed through the ua2nc decoder, the program
> thinks the time is 00 UTC on the 6th because the program rounds up the
> hour, but cannot round up the day. This error caused me to assimilate
> observations from the wrong time over China. Could you please come up
> with a solution to this problem? I have enclosed a Chinese rawinsonde
> observation as an example
>
> Ryan
>
> TTAA 03231 52203 99928 16662 00000 00103 ///// ///// 92768 20071 000^M^M
> 00 85501 20880 11505 70140 07663 25007 50580 12390 29512 40746 25128 ^M^M
> 27015 30949 41559 28022 25069 52358 28025 20213 53159 28025 15398 ^M^M
> 54963 29035 10653 599// 31024 88228 54558 28527 88109 601// 30026 ^M^M
> 77156 28535=^M^M
>
> --
> Ryan Torn
> Graduate Research Assistant
> Dept. of Atmospheric Sciences
> University of Washington, Box 351640
> Seattle, WA 98195
> Phone 206-685-1851
>
>
>
> Ticket Details
> ===================
> Ticket ID: YUU-262958
> Department: Support netCDF Decoders
> Priority: Normal
> Status: Open
> Link:
> http://www.unidata.ucar.edu/esupport/staff/index.php?_m=tickets&_a=viewticket&ticketid=1855
>
===============================================================================
Robb Kambic Unidata Program Center
Software Engineer III Univ. Corp for Atmospheric Research
address@hidden WWW: http://www.unidata.ucar.edu/
===============================================================================
#!/usr/bin/perl
#
# Decodes Upperair (WMO TEMP) reports according to:
# WMO Manual on codes volume 1 No. 306
# FM 35-IX Ext. TEMP - Upper-level pressure, temperture, humidity, and
# wind reports from a fixed land station.
#
# usage with input from STDIN: ua2nc [-v] [-h] [-n (old|new)] cdlfile
[datadir] [yyyymm]
#
# -v flag prints out the reports in the uaLog.log file, very verbose
#
use NetCDF ;
use Time::Local ;
#no encoding;
# process command line switches
$upperair = "new";
while( $_ = $ARGV[0], /^-/ ) {
shift ;
last if /^--$/ ;
/^(-v)/ && $verbose++ ;
/^(-h)/ && $hourly++; #create hourly files instead of daily
/^(-n)/ && ( $upperair = shift );
}
# new netCDF file naming covention is default
if( $upperair eq "new" ) {
if( $hourly ) {
$upperair = "00";
} else {
$upperair = "_0000";
}
$upperairPrefix = "Upperair_";
} else {
$upperair = "_upperair";
$upperairPrefix = "";
}
# process input parameters
if( $#ARGV == 0 ) {
$cdlfile = $ARGV[ 0 ] ;
} elsif( $#ARGV == 1 ) {
$cdlfile = $ARGV[ 0 ] ;
if( $ARGV[ 1 ] =~ /^\d/ ) {
$yyyymm = $ARGV[ 1 ] ;
} else {
$datadir = $ARGV[ 1 ] ;
}
} elsif( $#ARGV == 2 ) {
$cdlfile = $ARGV[ 0 ] ;
$datadir = $ARGV[ 1 ] ;
$yyyymm = $ARGV[ 2 ] ;
} else {
die "Usage with input on STDIN: ua2nc [-v] [-h] [-n (old|new)] cdlfile
[datadir] [yyyymm]\n" ;
}
# check for cdl and netCDF ncgen
die "Missing cdlfile parameter: $!\n" unless -e $cdlfile ;
if( -e "util/ncgen" ) {
$ncgen = "util/ncgen" ;
} elsif( -e "/usr/local/ldm/util/ncgen" ) {
$ncgen = "/usr/local/ldm/util/ncgen" ;
} elsif( -e "/upc/netcdf/bin/ncgen" ) {
$ncgen = "/upc/netcdf/bin/ncgen" ;
} elsif( -e "./ncgen" ) {
$ncgen = "./ncgen" ;
} else {
open( NCGEN, "which ncgen |" ) ;
$ncgen = <NCGEN> ;
close( NCGEN ) ;
if( $ncgen =~ /no ncgen/ ) {
die "Can't find NetCDF utility 'ncgen' in PATH, util/ncgen
/usr/local/ldm/util/ncgen, /upc/netcdf/bin/ncgen, or ./ncgen : $!\n" ;
} else {
$ncgen = "ncgen" ;
}
}
# the data directory
$datadir = "." if( ! $datadir ) ;
system( "mkdir -p $datadir" ) if( ! -e $datadir ) ;
if( -e "$datadir/uaLog.log.2" ) {
`rm -f $datadir/uaLog.log.3`;
`mv $datadir/uaLog.log.2 $datadir/uaLog.log.3`;
}
if( -e "$datadir/uaLog.log.1" ) {
`mv $datadir/uaLog.log.1 $datadir/uaLog.log.2`;
}
if( -e "$datadir/uaLog.log" ) {
`mv $datadir/uaLog.log $datadir/uaLog.log.1`;
}
# redirect STDOUT and STDERR
open( STDOUT, ">$datadir/uaLog.log" ) ||
die "could not open $datadir/uaLog.log: $!\n" ;
open( STDERR, ">&STDOUT" ) ||
die "could not dup stdout: $!\n" ;
select( STDERR ) ; $| = 1 ;
select( STDOUT ) ; $| = 1 ;
# year and month
if( ! $yyyymm ) {
$theyear = (gmtime())[ 5 ] ;
$theyear = ( $theyear < 100 ? $theyear : $theyear - 100 ) ;
$thedecade = sprintf( "%02d", $theyear ) ;
$theyear = "20" . sprintf( "%02d", $theyear ) ;
$themonth = (gmtime())[ 4 ] ;
$themonth++ ;
$themonth = sprintf( "%02d", $themonth ) ;
$yyyymm = "$theyear$themonth" ;
} else {
die "yyyymm must be 6 in length: $!\n" if( length( $yyyymm ) != 6 ) ;
$theyear = substr( $yyyymm, 0, 4 ) ;
$themonth = substr( $yyyymm, 4 ) ;
$thedecade = substr( $yyyymm, 2, 2 ) ;
}
%bin = ( "23", "00", "00", "00", "01", "00", "02", "03", "03", "03",
"04", "03", "05", "06", "06", "06", "07", "06", "08", "09",
"09", "09", "10", "09", "11", "12", "12", "12", "13", "12",
"14", "15", "15", "15", "16", "15", "17", "18", "18", "18",
"19", "18", "20", "21", "21", "21", "22", "21" ) ;
# set error handling to verbose only
$status = NetCDF::opts( VERBOSE ) ;
# set interrupt handler
$SIG{ 'INT' } = 'atexit' ;
$SIG{ 'KILL' } = 'atexit' ;
$SIG{ 'TERM' } = 'atexit' ;
$SIG{ 'QUIT' } = 'atexit' ;
# set defaults
%mand_ttaa = ("99", 99999, "00", 1000, "92", 925, "85", 850, "70", 700,
"50", 500, "40", 400, "30", 300, "25", 250, "20", 200,
"15", 150, "10", 100 ) ;
%mand_ttcc = ("70", 1, "50", 1, "30", 1, "20", 1, "10", 1, "07", 1, "05", 1,
"03", 1, "02", 1, "01", 1 ) ;
%mand_pres = ( "1000", 1, "925", 2, "850", 3, "700", 4, "500", 5, "400", 6,
"300", 7, "250", 8, "200", 9, "150", 10, "100", 11 , "70", 12,
"50", 13, "30", 14, "20", 15, "10", 16, "07", 17, "05", 18, "03", 19,
"02", 20, "01", 21 ) ;
%TTAA_top_wind_level = ( "/", 99999, "0", 1000, "9", 925, "8", 850, "7", 700,
"5", 500, "4", 400, "3", 300, "2", 200, "1", 100 ) ;
%TTCC_top_wind_level = ( "/", 100, "0", 0, "7", 70, "5", 50, "3", 30,
"2", 20, "1", 0) ;
# defaults
$manLevel = 22 ;
$sigTLevel = 65 ;
$sigWLevel = 50 ;
$mWndNum = 4 ;
$mTropNum = 4 ;
$staNameLen = 6 ;
# open cdl and get dimensions for variables
open( CDL, "$cdlfile" ) || die "could not open $cdlfile: $!\n" ;
$i = 0 ;
while( <CDL> ) {
if( s#^\s*manLevel\s*=\s*(\d{1,5})## ) {
$manLevel = $1 ;
} elsif( s#^\s*sigTLevel\s*=\s*(\d{1,5})## ) {
$sigTLevel = $1 ;
} elsif( s#^\s*sigWLevel\s*=\s*(\d{1,5})## ) {
$sigWLevel = $1 ;
} elsif( s#^\s*mWndNum\s*=\s*(\d{1,5})## ) {
$mWndNum = $1 ;
} elsif( s#^\s*mTropNum\s*=\s*(\d{1,5})## ) {
$mTropNum = $1 ;
} elsif( s#^\s*staNameLen\s*=\s*(\d{1,5})## ) {
$staNameLen = $1 ;
} elsif( s#^\s*variables## ) {
last ;
}
}
close CDL ;
# initialize record buffer and variables
makedataref();
# read in station data
if( -e "etc/snstns.tbl" ) {
$sfile = "etc/snstns.tbl" ;
} elsif( -e "./snstns.tbl" ) {
$sfile = "./snstns.tbl" ;
} else {
die "Can't find snstns.tbl station file.: $!\n" ;
}
open( STATION, "$sfile" ) || die "could not open $sfile: $!\n" ;
while( <STATION> ) {
s#^(\w{3,6})?\s+(\d{4,5}).{40}## ;
$id = $1 ;
$wmo_id = $2 ;
$wmo_id = "0" . $wmo_id if( length( $wmo_id ) == 4 ) ;
( $lat, $lon, $elev ) = split ;
$lat = sprintf( "%7.2f", $lat / 100 ) ;
$lon = sprintf( "%7.2f", $lon / 100) ;
# set these vars ( $lat, $lon, $elev, $id )
$STATIONS{ "$wmo_id" } = "$lat $lon $elev $id" ;
}
close STATION ;
# read in list of already processed reports if it exists
# open ua.lst, list of reports processed in the last 24 hours.
if( -e "$datadir/ua.lst" ) {
open( LST, "$datadir/ua.lst" ) ||
die "could not open $datadir/ua.lst: $!\n" ;
while( <LST> ) {
( $stn, $yyyymmddhh, $record, @R ) = split ;
$rpt_hash{ "$stn $yyyymmddhh" } = "$record @R" ;
}
close LST ;
}
# Now begin parsing file and decoding observations breaking on cntrl C
$/ = "\cC" ;
# set select processing here from STDIN
START:
while( 1 ) {
open( STDIN, '-' ) ;
vec( $rin,fileno( STDIN ),1 ) = 1 ;
$timeout = 1200 ; # 20 minutes
$nfound = select( $rout = $rin, undef, undef, $timeout ) ;
# timed out
if( ! $nfound ) {
print "Shut down, time out $timeout seconds on select\n" ;
atexit() ;
}
atexit( "eof" ) if( eof( STDIN )) ;
# Process each line of ua bulletins, header first
$_ = <STDIN> ;
next unless( /TT|PP/ ) ;
s#\cC## ;
s#\cM##g ;
s#\cA\n## ;
s#\c^##g ;
s#^\n##g ;
s#\d\d\d \n## ;
s#\w{4}\d{1,2} \w{4} (\d{2})(\d{2})(\d{2}).*\n## ;
$bday = $1 ;
$bhour = $2 ;
$bhour = "23" if( $bhour eq "24" ) ;
$bmin = $3 ;
next unless ( $bday && defined( $bhour ) && defined( $bmin ) ) ;
next if( $bmin > 59 || $bhour > 23 || $bday > 31 ) ;
# check for valid transmission times against current time
$cday = (gmtime())[ 3 ] ;
$chour = (gmtime())[ 2 ] ;
# skip bulletins over 24 hours old or in the future
if( $bday == $cday ) {
next if( $bhour > $chour ) ;
} else { # $bday != $cday, skip over day old reports
next if( $bday < $cday -1 ) ;
if( $bday > $cday ) {
next if( $cday != 1 || $bday < 28) ;
}
next if( $bhour < $chour ) ;
}
# reset year and month
$theyear = (gmtime())[ 5 ] ;
$theyear = ( $theyear < 100 ? $theyear : $theyear - 100 ) ;
$thedecade = sprintf( "%02d", $theyear ) ;
$theyear = "20" . sprintf( "%02d", $theyear ) ;
$themonth = (gmtime())[ 4 ] ;
$themonth++ ;
$themonth = sprintf( "%02d", $themonth ) ;
$yyyymm = "$theyear$themonth" ;
# Separate bulletins into reports
if( /=\n/ ) {
s#=\s+\n#=\n#g ;
} else {
s#\n# #g ;
}
@reports = split( /=\n/ ) ;
for( @reports ) { # Process each report in the bulletin
undef( $stn_name ) ;
next if( /NIL/ ) ;
s#\n# #g ;
s#^\d{5}.*((TT)(AA|BB|CC|DD)|(PP)(BB|DD))#$1# ;
next unless (
s# ?((TT)(AA|BB|CC|DD)|(PP)(BB|DD)) {1,2}(\d\d)(\d\d)([\d\/]) (\d{5}) ## ) ;
#print "$&\n" ;
$type = $1 ;
$rday = $6 ;
$rhour = $7 ;
$indicator = $8 ;
$stn = $9 ;
if( $rday =~ /^[5678]/ ) {
$rday -= 50 ;
$knots = 1 ;
} else {
$knots = 0 ;
}
# check for valid times
next unless ($rday && defined( $rhour )) ;
next if( $rhour > 23 || $rday > 31 ) ;
# skip reports over 24 hours old
$tmpyyyymm = $yyyymm ;
if( $rday == $cday ) {
next if( $rhour > $chour ) ;
# increment rday & check if last day of month
if( $rhour == 23 ) {
incrementRDAY() ;
$tmpyyyymm = $theyear . $themonth ;
}
} else { # $rday != $cday, skip over day old reports
next if( $rday < ( $cday -1 ) ) ;
if( $rday > $cday ) {
next if( $cday != 1 || $rday < 28) ;
# cday = 1, reset month and year
$tmpmonth = sprintf( "%02d", $themonth -1 ) ;
if( $tmpmonth == 0 ) {
$tmpmonth = "12" ;
$tmpyear =
sprintf( "%04d", $theyear -1 );
} else {
$tmpyear = $theyear ;
}
$tmpyyyymm = $tmpyear . $tmpmonth ;
}
next if( $rhour < $chour ) ;
}
$yyyymmddhh = $tmpyyyymm . sprintf( "%02d", $rday ) .
$bin{ "$rhour" } ;
if( defined( $STATIONS{ "$stn" } )) {
# station's lat, lon, elev, id.
( $lat, $lon, $elev, $id ) =
split( ' ', $STATIONS{ "$stn" } ) ;
$id = $AS6 unless( defined( $id )) ;
} else {
print "WMO Id $stn not in station table\n" ;
next ;
}
if( defined( $rpt_hash{ "$stn $yyyymmddhh" } )) {
( $record, @R ) = split( ' ',
$rpt_hash{ "$stn $yyyymmddhh" } ) ;
$rpt_hash{ "$stn $yyyymmddhh" } .= " $type"
if( index( "@R", $type ) == -1 ) ;
} else {
$record = -1 ;
}
# clean data
s#\s+# #g ;
$_ .= ' ' ;
$j = 0 ;
# process reports
SWITCH: {
ttaa( $_ ), last SWITCH if( $type eq "TTAA" ) ;
ttbb( $_ ), last SWITCH if( $type eq "TTBB" ) ;
ttcc( $_ ), last SWITCH if( $type eq "TTCC" ) ;
ttdd( $_ ), last SWITCH if( $type eq "TTDD" ) ;
ppbbdd( $_ ), last SWITCH
if( $type eq "PPBB" || $type eq "PPDD" ) ;
}
# no data extracted
next unless $j ;
# covert knots to meters/second if necessary
if( $knots && $#wspd != -1 ) {
for( $i = 0; $i < $j; $i++ ) {
$wspd[ $i ] *= 0.5144 if( $wspd[ $i ] != $F ) ;
}
}
# we have a legal report, open a Netcdf file
$status = doNet() ;
# print output
printvars() if( $verbose ) ;
# output the NetCDF data here
makedataref();
# record already exists for $stn
if( $record >= 0 ) {
$status = NetCDF::recget( $ncid, $record, \@dataref ) ;
print "NetCDF::recget status = $status\n"
if( $status ) ;
} else { # initialize new record
$record = $recnum ;
$rpt_hash{ "$stn $yyyymmddhh" } = "$record $type" ;
$recnum++ ;
$dataref[ 0 ] = \$stn ;
if( defined( $id )) {
$id = padstr( $id, $staNameLen ) ;
$dataref[ 1 ] = \$id ;
$dataref[ 2 ] = \$lat ;
$dataref[ 3 ] = \$lon ;
$dataref[ 4 ] = \$elev ;
}
$dataref[ 5 ] = \theTime( "synoptic" ) ;
$dataref[ 11 ] = \theTime( "release" ) ; # release time
}
setvars() ;
$status = NetCDF::recput( $ncid, $record, \@dataref ) ;
if( $status ) { # failure
print "Wmo =$stn, Type =$type, Rpt: $rday$rhour, ",
"Trans: $bday$bhour, recput =$status\n" ;
} else {
$status = NetCDF::sync( $ncid ) ;
#print "Syncing $ncfile with ncid $ncid\n" ;
}
atexit( "eof" ) if( eof( STDIN )) ;
if( $type eq "TTAA" || $type eq "TTCC" ) {
undef( $numMand ) ;
undef( @prMan ) ;
undef( @htMan ) ;
undef( @tpMan ) ;
undef( @tdMan ) ;
undef( @wdMan ) ;
undef( @wsMan ) ;
undef( $numTrop ) ;
undef( @prTrop ) ;
undef( @tpTrop ) ;
undef( @tdTrop ) ;
undef( @wdTrop ) ;
undef( @wsTrop ) ;
undef( $numMwnd ) ;
undef( @prMaxW ) ;
undef( @wdMaxW ) ;
undef( @wsMaxW ) ;
} elsif( $type eq "TTBB" || $type eq "TTDD" ) {
undef( $numSigT ) ;
undef( @prSigT ) ;
undef( @tpSigT ) ;
undef( @tdSigT ) ;
undef( %prSigT ) ;
} elsif( $type eq "PPBB" || $type eq "PPDD" ) {
undef( $numSigW ) ;
undef( @htSigW ) ;
undef( @wsSigW ) ;
undef( @wdSigW ) ;
undef( %htSigW ) ;
}
undef( $i ) ;
undef( $j ) ;
undef( $k ) ;
undef( @pres ) ;
undef( @ht ) ;
undef( @t ) ;
undef( @td ) ;
undef( @wdir ) ;
undef( @wspd ) ;
undef( $trop_pres ) ;
undef( $trop_t ) ;
undef( $trop_td_dep ) ;
undef( $trop_wdir ) ;
undef( $trop_wspd ) ;
undef( $max_pres ) ;
undef( $max_wdir ) ;
undef( $max_wspd ) ;
} # end foreach report
} # end while( 1 )
atexit( "eof" ) ;
exit( 0 ) ; #should never get here
# setvars into record
sub setvars
{
# enter report data into proper record location
SWITCH1: {
( $type eq "TTAA" || $type eq "TTCC" ) && do {
$numMand = ${$dataref[ 6 ]} ; # number of mandatory
# data already inserted into record, append/overwrite
if( $numMand != $F ) {
@prMan = @{ $dataref[ 13 ] } ;
@htMan = @{ $dataref[ 14 ] } ;
@tpMan = @{ $dataref[ 15 ] } ;
@tdMan = @{ $dataref[ 16 ] } ;
@wdMan = @{ $dataref[ 17 ] } ;
@wsMan = @{ $dataref[ 18 ] } ;
} else { # first data write
@prMan = @ML ;
@htMan = @ML ;
@tpMan = @ML ;
@tdMan = @ML ;
@wdMan = @ML ;
@wsMan = @ML ;
}
if( $pres[ 1 ] == 1000 ) {
$prMan[ 0 ] = $pres[ 0 ] ;
$htMan[ 0 ] = $ht[ 0 ] ;
$tpMan[ 0 ] = $t[ 0 ] ;
$tdMan[ 0 ] = $td[ 0 ] ;
$wdMan[ 0 ] = $wdir[ 0 ] ;
$wsMan[ 0 ] = $wspd[ 0 ] ;
$j = 1 ;
} else {
$j = 0 ;
}
for( $i = $j; $i < $manLevel; $i++ ) {
$k = $mand_pres{ "$pres[ $i ]" } ;
next unless( defined( $k )) ;
$prMan[ $k ] = $pres[ $i ] ;
$htMan[ $k ] = $ht[ $i ] ;
$tpMan[ $k ] = $t[ $i ] ;
$tdMan[ $k ] = $td[ $i ] ;
$wdMan[ $k ] = $wdir[ $i ] ;
$wsMan[ $k ] = $wspd[ $i ] ;
}
$numMand = 0 ;
for( $i = 0 ; $i < $manLevel; $i++ ) {
$numMand++ if( $prMan[ $i ] != $F ) ;
}
$dataref[ 6 ] = \$numMand ;
$dataref[ 13 ] = \@prMan ;
$dataref[ 14 ] = \@htMan ;
$dataref[ 15 ] = \@tpMan ;
$dataref[ 16 ] = \@tdMan ;
$dataref[ 17 ] = \@wdMan ;
$dataref[ 18 ] = \@wsMan ;
# enter topopause data
$numTrop = ${$dataref[ 10 ]} ; # number of tropopause
# data already inserted into record, append/overwrite
if( $numTrop != $F ) {
@prTrop = @{ $dataref[ 25 ] } ;
@tpTrop = @{ $dataref[ 26 ] } ;
@tdTrop = @{ $dataref[ 27 ] } ;
@wdTrop = @{ $dataref[ 28 ] } ;
@wsTrop = @{ $dataref[ 29 ] } ;
} else { # first data write
@prTrop = @TL ;
@tpTrop = @TL ;
@tdTrop = @TL ;
@wdTrop = @TL ;
@wsTrop = @TL ;
}
# enter new tropopause data
if( defined( $trop_pres ) && $numTrop == $F ) {
$prTrop[ 0 ] = $trop_pres ;
$tpTrop[ 0 ] = $trop_t ;
$tdTrop[ 0 ] = $trop_td_dep ;
$wdTrop[ 0 ] = $trop_wdir ;
$wsTrop[ 0 ] = $trop_wspd ;
$numTrop = 1;
} elsif( defined( $trop_pres ) ) {
for( $i = 0; $i < $mTropNum; $i++ ) {
if( $prTrop[ $i ] == $trop_pres ) {
$index = $i ;
last ;
} elsif( $prTrop[ $i ] == $F ) {
$numTrop++ ;
$index = $i ;
last ;
}
}
$prTrop[ $index ] = $trop_pres ;
$tpTrop[ $index ] = $trop_t ;
$tdTrop[ $index ] = $trop_td_dep ;
$wdTrop[ $index ] = $trop_wdir ;
$wsTrop[ $index ] = $trop_wspd ;
}
$dataref[ 25 ] = \@prTrop ;
$dataref[ 26 ] = \@tpTrop ;
$dataref[ 27 ] = \@tdTrop ;
$dataref[ 28 ] = \@wdTrop ;
$dataref[ 29 ] = \@wsTrop ;
$dataref[ 10 ] = \$numTrop ;
# enter maximum wind data
$numMwnd = ${$dataref[ 9 ]} ; # number of max winds
# data already inserted into record, append/overwrite
if( $numMwnd != $F ) {
@prMaxW = @{ $dataref[ 30 ] } ;
@wdMaxW = @{ $dataref[ 31 ] } ;
@wsMaxW = @{ $dataref[ 32 ] } ;
} else { # first data write
@prMaxW = @MWL ;
@wdMaxW = @MWL ;
@wsMaxW = @MWL ;
}
# enter new max wind data
if( defined( $max_pres ) && $numMwnd == $F ) {
$prMaxW[ 0 ] = $max_pres ;
$wdMaxW[ 0 ] = $max_wdir ;
$wsMaxW[ 0 ] = $max_wspd ;
$numMwnd = 1;
} elsif( defined( $max_pres ) ) {
for( $i = 0; $i < $mWndNum; $i++ ) {
if( $prMaxW[ $i ] == $max_pres ) {
$index = $i ;
last ;
} elsif( $prMaxW[ $i ] == $F ) {
$numMwnd++ ;
$index = $i ;
last ;
}
}
$prMaxW[ $index ] = $max_pres ;
$wdMaxW[ $index ] = $max_wdir ;
$wsMaxW[ $index ] = $max_wspd ;
}
$dataref[ 30 ] = \@prMaxW ;
$dataref[ 31 ] = \@wdMaxW ;
$dataref[ 32 ] = \@wsMaxW ;
$dataref[ 9 ] = \$numMwnd ;
last SWITCH1 ;
} ;
( $type eq "TTBB" || $type eq "TTDD" ) && do {
$k = ${$dataref[ 7 ]} ; # number of Sig Temps
# data already inserted into record, overwrite or
# append
if( $k != $F ) {
@prSigT = @{ $dataref[ 19 ] } ;
@tpSigT = @{ $dataref[ 20 ] } ;
@tdSigT = @{ $dataref[ 21 ] } ;
for( $i = 0; $i < $k; $i++ ) {
$prSigT[ $i ] = sprintf( "%3.1f",
$prSigT[ $i ] ) ;
$prSigT{ "$prSigT[ $i ]" } = $i ;
}
for( $i = 0; $i <= $#pres; $i++ ) {
if( defined( $prSigT{ "$pres[ $i ]"} )) {
$j = $prSigT{ "$pres[ $i ]" } ;
$tpSigT[ $j ] = $t[ $i ] ;
$tdSigT[ $j ] = $td[ $i ] ;
} else { # append
$prSigT[ $k ] = $pres[ $i ] ;
$tpSigT[ $k ] = $t[ $i ] ;
$tdSigT[ $k ] = $td[ $i ] ;
$k++ ;
}
}
while( $k > $sigTLevel ) {
$k-- ;
pop( @prSigT ) ;
pop( @tpSigT ) ;
pop( @tdSigT ) ;
print "High sigTLevel $stn $type $k\n" ;
}
$dataref[ 7 ] = \$k ;
$dataref[ 19 ] = \@prSigT ;
$dataref[ 20 ] = \@tpSigT ;
$dataref[ 21 ] = \@tdSigT ;
} else { # first data entry
$numSigT = $#pres +1 ;
$dataref[ 7 ] = \$numSigT ;
for( $i = $numSigT ; $i < $sigTLevel; $i++ ) {
$pres[ $i ] = $F ;
$t[ $i ] = $F ;
$td[ $i ] = $F ;
}
$dataref[ 19 ] = \@pres ;
$dataref[ 20 ] = \@t ;
$dataref[ 21 ] = \@td ;
}
$dataref[ 12 ] = \$indicator if( $type eq "TTBB" ) ;
last SWITCH1 ;
} ;
( $type eq "PPBB" || $type eq "PPDD" ) && do {
$k = ${$dataref[ 8 ]} ; # number of Sig Winds
# data already inserted into record, overwrite or
# append
if( $k != $F ) {
@htSigW = @{ $dataref[ 22 ] } ;
@wdSigW = @{ $dataref[ 23 ] } ;
@wsSigW = @{ $dataref[ 24 ] } ;
for( $i = 0; $i < $k; $i++ ) {
$htSigW[ $i ] = sprintf( "%3.1f",
$htSigW[ $i ] ) ;
$htSigW{ "$htSigW[ $i ]" } = $i ;
}
for( $i = 0; $i <= $#ht; $i++ ) {
if( defined( $htSigW{ "$ht[ $i ]" } )) {
$j = $htSigW{ "$ht[ $i ]" } ;
$wsSigW[ $j ] = $wspd[ $i ] ;
$wdSigW[ $j ] = $wdir[ $i ] ;
} else { # append
$htSigW[ $k ] = $ht[ $i ] ;
$wsSigW[ $k ] = $wspd[ $i ] ;
$wdSigW[ $k ] = $wdir[ $i ] ;
$k++ ;
}
}
while( $k > $sigWLevel ) {
$k-- ;
pop( @htSigW ) ;
pop( @wsSigW ) ;
pop( @wdSigW ) ;
print "High sigWLevel $stn $type $k\n" ;
}
$dataref[ 8 ] = \$k ;
$dataref[ 22 ] = \@htSigW ;
$dataref[ 23 ] = \@wdSigW ;
$dataref[ 24 ] = \@wsSigW ;
} else { # first data entry
$numSigW = $#wdir +1 ;
$dataref[ 8 ] = \$numSigW ;
for( $i = $numSigW; $i < $sigWLevel; $i++ ) {
$ht[ $i ] = $F ;
$wdir[ $i ] = $F ;
$wspd[ $i ] = $F ;
}
$dataref[ 22 ] = \@ht ;
$dataref[ 23 ] = \@wdir ;
$dataref[ 24 ] = \@wspd ;
}
last SWITCH1 ;
} ;
}
} # end setvars
# create a netcdf file or reopen a existing one
sub doNet
{
my( $Ncfile, $Id, $Num, $Time, $baseTime, $offset, $rpt, $stn ) ;
#$ncfile = $datadir . "/" . $yyyymmddhh . "_perl_ua.nc" ;
$ncfile = $datadir . "/" . $upperairPrefix . substr( $yyyymmddhh, 0, 8 )
. $upperair . ".nc" ;
if( $hourly ) {
$ncfile = $datadir . "/" . $upperairPrefix . substr( $yyyymmddhh, 0, 8
) .
"_" . substr( $yyyymmddhh, 8 ) . $upperair . ".nc" ;
}
# writing to same file
return 1 if( $ncfile eq $lastNc ) ;
# current time
$thetime = time() ;
# save current file info
$Nets{ $lastNc } = "$ncid $recnum $thetime" if( $lastNc ) ;
# File is open, get ncfile id and recnum and reset the time
if( defined( $Nets{ $ncfile } )) { # already open for writes
( $ncid, $recnum, $ncTime ) = split( " ", $Nets{ $ncfile } ) ;
$ncTime = $thetime ;
$lastNc = $ncfile ;
return 1 ;
}
# close files with no activity for 20 minutes
foreach $Ncfile ( keys %Nets ) {
( $Id, $Num, $Time ) = split( " ", $Nets{ $Ncfile } ) ;
if( $thetime - $Time > 1200 ) {
print "Closing $Ncfile with ncid $Id, No write for > 20
Minutes\n" ;
$status = NetCDF::close( $Id ) ;
delete( $Nets{ $Ncfile } ) ;
# } elsif( $Ncfile eq $lastNc ) {
# print "Syncing $lastNc with ncid $ncid\n" ;
# $status = NetCDF::sync( $ncid ) ;
}
}
# remove rpt entries older than 24 hours
$baseTime = $yyyymm . $cday . $chour ;
foreach $rpt ( keys %rpt_hash ) {
( $stn, $Time ) = split( " ", $rpt ) ;
$offset = $baseTime - $Time ;
next if( $offset < 100 ) ; # same day ok & previous day ok > $chour
delete( $rpt_hash{ $rpt } ) ;
}
# open or create ncfiles
if( -e $ncfile ) {
$ncid = NetCDF::open( "$ncfile", WRITE ) ;
return 0 if( $ncid == -1 ) ;
$recNum_id = NetCDF::dimid( $ncid, "recNum" ) ;
$name_id = "xxxxxxxx" ;
$recnum = -1 ;
# get current value of recnum
NetCDF::diminq( $ncid, $recNum_id, $name_id, $recnum ) ;
} else {
system( "$ncgen -o $ncfile $cdlfile" ) ;
$ncid = NetCDF::open( "$ncfile", WRITE ) ;
return 0 if( $ncid == -1 ) ;
# NetCDF record counter
$recnum = 0 ;
}
$Nets{ $ncfile } = "$ncid $recnum $thetime" ;
$lastNc = $ncfile ;
print "Opening $ncfile with ncid $ncid\n" ;
return 1 ;
} # end doNet
# print vars
sub printvars
{
my( @rec ) ;
for( $i=0; $i<$j; $i++ ) {
$rec[$i] = sprintf "%6s%8s%8s%10s%10s%13s",
$pres[$i], $ht[$i], $t[$i], $td[$i], $wdir[$i], $wspd[$i] ;
}
# Output decoded sounding
print "Station: $stn, report type: $type, j = $j\n" ;
print "Time: $yyyymmddhh\n" ;
print "---------------------------------------------------------\n" ;
print "Pres Ht Temp DewPt Dir Spd\n" ;
print " mb m C C deg m/s \n" ;
print "---------------------------------------------------------\n" ;
for( $i = 0; $i < $j; $i++ ) { print "$rec[$i]\n" ; }
print "\n" ;
if( $trop_t ) {
print "Tropopause values\n" ;
print "---------------------------------------------------------\n" ;
# trop_t, trop_td_dep trop_wdir, trop_wspd
$trop = sprintf "%6s%15s%10s%10s%13s",
$trop_pres, $trop_t, $trop_td_dep, $trop_wdir, $trop_wspd ;
print "$trop\n\n" ;
}
if( $max_wdir ) {
# max_pres, max_wdir, max_wspd
print "Maximum values\n" ;
print "---------------------------------------------------------\n" ;
$max = sprintf "%6s%35s%13s", $max_pres, $max_wdir, $max_wspd ;
print "$max\n\n" ;
}
} # end printvars
# execute at exit
sub atexit
{
my( $sig ) = @_ ;
if( $sig eq "eof" ) {
print "eof on STDIN --shutting down\n" ;
} elsif( defined( $sig )) {
print "Caught SIG$sig --shutting down\n" ;
}
# open ua.lst, list of reports processed in the last 8 hours.
open( LST, ">$datadir/ua.lst" ) ||
die "could not open $datadir/ua.lst: $!\n" ;
select( LST ) ;
# remove rpt entries older than 24 hours
$yyyymmddhh = $yyyymm . $cday . $chour ;
foreach $rpt ( keys %rpt_hash ) {
( $stn, $Time ) = split( " ", $rpt ) ;
$offset = $yyyymmddhh - $Time ;
next unless( $offset < 100 ) ; # same day ok & previous day ok > $chour
print "$rpt $rpt_hash{ $rpt }\n" ;
}
close LST ;
foreach $Ncfile ( keys %Nets ) {
( $ncid, $recnum, $nctime ) = split( " ", $Nets{ $Ncfile } ) ;
print STDOUT "Closing $Ncfile with ncid $ncid\n" ;
$status = NetCDF::close( $ncid ) ;
print "STDOUT NetCDF::close status = $status\n" if( $status ) ;
}
close( STDOUT ) ;
close( STDERR ) ;
exit( 0 ) ;
} # end atexit
# pad str to correct length
sub padstr
{
( $str, $len ) = @_ ;
my( $size, $i ) ;
$size = length( $str ) ;
for( $i = $size; $i < $len; $i++ ) {
$str .= "\0" ;
#print "$str,\n" ;
}
if( $size > $len ) {
print STDOUT "String length is over $len chars long:\n $str\n"
if( $verbose ) ;
$str = substr( $str, 0, $len ) ;
}
return $str ;
} # end padstr
sub makedataref
{
undef( @dataref ) ;
$F = 99999 ;
$F0 = 99999 ;
$F1 = 99999 ;
$F2 = 99999 ;
$F3 = 99999 ;
$F4 = 99999 ;
$F5 = 99999 ;
$F6 = 99999 ;
$F7 = 99999 ;
$F8 = 99999 ;
$F85 = 99999 ;
$F9 = 99999 ;
$F10 = 99999 ;
$A = \$F ;
$S6 = "\0" x $staNameLen ;
$AS6 = \$S6 ;
@ML = ( 99999 ) x ( $manLevel ) ;
@ML0 = ( 99999 ) x ( $manLevel ) ;
@ML1 = ( 99999 ) x ( $manLevel ) ;
@ML2 = ( 99999 ) x ( $manLevel ) ;
@ML3 = ( 99999 ) x ( $manLevel ) ;
@ML4 = ( 99999 ) x ( $manLevel ) ;
@ML5 = ( 99999 ) x ( $manLevel ) ;
@STL = ( 99999 ) x ( $sigTLevel ) ;
@STL0 = ( 99999 ) x ( $sigTLevel ) ;
@STL1 = ( 99999 ) x ( $sigTLevel ) ;
@STL2 = ( 99999 ) x ( $sigTLevel ) ;
@SWL = ( 99999 ) x ( $sigWLevel ) ;
@SWL0 = ( 99999 ) x ( $sigWLevel ) ;
@SWL1 = ( 99999 ) x ( $sigWLevel ) ;
@SWL2 = ( 99999 ) x ( $sigWLevel ) ;
@TL = ( 99999 ) x ( $mTropNum ) ;
@TL0 = ( 99999 ) x ( $mTropNum ) ;
@TL1 = ( 99999 ) x ( $mTropNum ) ;
@TL2 = ( 99999 ) x ( $mTropNum ) ;
@TL3 = ( 99999 ) x ( $mTropNum ) ;
@TL4 = ( 99999 ) x ( $mTropNum ) ;
@MWL = ( 99999 ) x ( $mWndNum ) ;
@MWL0 = ( 99999 ) x ( $mWndNum ) ;
@MWL1 = ( 99999 ) x ( $mWndNum ) ;
@MWL2 = ( 99999 ) x ( $mWndNum ) ;
# default netCDF record structure, contains all vars for the UA reports
@dataref = ( \$F0, \$S6, \$F1, \$F2, \$F3, \$F4, \$F5, \$F6, \$F7, \$F8,
\$F85, \$F9, \$F10, \@ML0, \@ML1, \@ML2, \@ML3, \@ML4, \@ML5, \@STL0, \@STL1,
\@STL2, \@SWL0, \@SWL1, \@SWL2, \@TL0, \@TL1, \@TL2, \@TL3, \@TL4, \@MWL0,
\@MWL1, \@MWL2 ) ;
}
# TTAA type of reports, mandatory levels: surface to 100 mb
sub ttaa
{
my( $TTAA, $stride, $toplevel, $offset, $thermo, $wind ) ;
( $TTAA ) = @_ ;
$stride = 18 ; # when winds present in report
$toplevel = $TTAA_top_wind_level{ "$indicator" } ;
# mandatory levels, sfc thru 100 mb
for( $offset = 0 ; $offset < length( $TTAA ) ; $offset += $stride, $j++ ) {
$pres[ $j ] = $ht[ $j ] = $t[ $j ] = $td[ $j ] = $wdir[ $j ] =
$wspd[ $j ] = $F ;
$pres[ $j ] = substr( $TTAA, $offset, 2 ) ;
last unless( defined( $mand_ttaa{ "$pres[ $j ]" } ) ) ;
$pres[ $j ] = $mand_ttaa{ "$pres[ $j ]" } ;
# $stride is set to 12 when no winds are reported
$stride = 12 if( $pres[ $j ] < $toplevel ) ;
if( substr( $TTAA, $offset+2, 10 ) =~ m#([\d/]{3}) ([\d/]{5}) # ) {
$ht[ $j ] = $1 ;
dec_thermo( $2 ) ;
} else {
$ht[ $j ] = "///" ;
}
if(( $stride == 18 || $pres[ $j ] == 99999 ) &&
substr( $TTAA, $offset + 12, 6 ) =~ m#([\d/]{5}) # ) {
dec_wind( $1 ) ;
}
if( $pres[ $j ] == 99999 ) { # surface
$sfc_pres = $pres[ $j ] =
( $ht[ $j ] < 100 ) ? $ht[ $j ] + 1000 : $ht[ $j ] ;
$ht[ $j ] = $elev ;
}
$ht[ $j ] .= "0" if( $pres[ $j ] <= 500 ) ; # 500 mb and above
if( $ht[ $j ] =~ m#/# ) {
$ht[ $j ] = $F ;
# 1000 mb special case for below sea level
} elsif( $pres[ $j ] == 1000 && $pres[ $j ] > $sfc_pres ) {
if( $ht[ $j ] > 500 ) {
$ht[ $j ] -= 500 ;
$ht[ $j ] *= -1 ;
}
# 925 mb special case for below sea level
} elsif( $pres[ $j ] == 925 && $pres[ $j ] > $sfc_pres ) {
if( $ht[ $j ] > 500 ) {
$ht[ $j ] -= 500 ;
$ht[ $j ] *= -1 ;
}
} elsif( $pres[ $j ] == 925 && $ht[ $j ] < 250 ) { # 925 mb
$ht[ $j ] = "1" . $ht[ $j ] ;
} elsif( $pres[ $j ] == 850 ) { # 850 mb
$ht[ $j ] = "1" . $ht[ $j ] ;
} elsif( $pres[ $j ] == 700 ) { # 700 mb
if( $ht[ $j ] < 350 ) {
$ht[ $j ] = "3" . $ht[ $j ] ;
} else {
$ht[ $j ] = "2" . $ht[ $j ] ;
}
} elsif( $pres[ $j ] == 250 && $ht[ $j ] < 5000 ) { # 250 mb
$ht[ $j ] = "1" . $ht[ $j ] ;
} elsif( $pres[ $j ] =~ /200|150|100$/ ) { # 200, 150, &100 mb
$ht[ $j ] = "1" . $ht[ $j ] ;
}
} # End of loop thru TTAA mandatory levels
pop( @pres ) unless( $mand_pres{ "$pres[ $j ]" } ) ;
# Tropopause near end of TTAA
if( $TTAA !~ / 88999/ &&
$TTAA =~ m# 88(\d\d\d) ([\d/]{5}) ([\d/]{5}) # ) {
$trop_pres = $1 ;
$thermo = $2 ;
$wind = $3 ;
# sets trop_t, trop_td_dep
trop_dec_thermo( $thermo ) ;
# sets trop_wdir, trop_wspd
trop_dec_wind( $wind ) ;
} # End of tropopause level in TTAA group
# Maximum Wind
if( $TTAA !~ / 77999/ &&
$TTAA =~ m# (66|77)(\d\d\d) ([\d/]{5}) # ) {
$max_pres = $2 ;
# sets max_wdir, max_wspd
max_dec_wind( $3 ) ;
} # End of maximum wind level in TTAA group
} #end ttaa
# TTBB reports, significant levels with respect to temp: surface to 100 mb
sub ttbb
{
my( $TTBB, $offset ) ;
( $TTBB ) = @_ ;
$indicator = $F if( $indicator eq '/' ) ;
$offset = 0 ;
while( $offset < length( $TTBB )) {
last unless substr( $TTBB, $offset, 2 ) =~ /00|11|22|33|44|55|66|77|88/;
$pres[ $j ] = $ht[ $j ] = $t[ $j ] = $td[ $j ] = $wdir[ $j ] =
$wspd[ $j ] = $F ;
if( substr( $TTBB, $offset+2, 10 ) =~ /(\d{3}) (\d{5}) / ) {
$pres[ $j ] = $1 ;
dec_thermo( $2 ) ;
$pres[ $j ] = "1" . $pres[ $j ] if( $pres[ $j ] < 100 ) ;
if( defined( $mand_pres{ "$pres[ $j ]" } )) {
$offset += 12 ;
next ;
}
$pres[ $j ] = sprintf( "%3.1f", $pres[ $j ] ) ;
$j++ ;
}
$offset += 12 ;
} # End of TTBB
} #end ttbb
# TTCC type of reports, mandatory levels: under 100 mb
sub ttcc
{
my( $TTCC, $stride, $toplevel, $offset, $thermo, $wind ) ;
( $TTCC ) = @_ ;
$stride = 18 ; # when winds present in report
$toplevel = $TTCC_top_wind_level{ "$indicator" } ;
if( ! defined( $toplevel )) {
print "Undefined Wind termination level: $level, $stn $type ",
"Rpt: $rday$rhour Trans: $bday$bhour\n" ;
return ;
}
# mandatory levels, 70 mb thru 1 mb
for( $offset = 0 ; $offset < length( $TTCC ) ; $offset += $stride, $j++ ) {
$pres[ $j ] = $ht[ $j ] = $t[ $j ] = $td[ $j ] = $wdir[ $j ] =
$wspd[ $j ] = $F ;
$pres[ $j ] = substr( $TTCC, $offset, 2 ) ;
last unless( defined( $mand_ttcc{ "$pres[ $j ]" } ) ) ;
# $stride is set to 12 when no winds are reported
$stride = 12 if( $pres[ $j ] < $toplevel ) ;
if( substr( $TTCC, $offset+2, 10 ) =~ m#([\d/]{3}) ([\d/]{5}) # ) {
$ht[ $j ] = $1 ;
dec_thermo( $2 ) ;
} else {
$ht[ $j ] = "///" ;
}
if( $stride == 18 &&
substr( $TTCC, $offset +12, 6 ) =~ m#([0-3/][\d/]{4}) # ) {
dec_wind( $1 ) ;
}
if( $ht[ $j ] =~ m#/# ) {
$ht[ $j ] = $F ;
} else {
$ht[ $j ] .= "0" ;
}
if( $pres[ $j ] == 70 ) { # 70 mb
$ht[ $j ] = "1" . $ht[ $j ] ;
} elsif( $pres[ $j ] =~ /50|30|20/ ) { # 50,30,20 mb
$ht[ $j ] = "2" . $ht[ $j ] ;
} else { # 10, 7, 5, 3, 2, 1 mb
$ht[ $j ] = "3" . $ht[ $j ] ;
}
} # end of loop thru TTCC mandatory levels
pop( @pres ) unless( $mand_pres{ "$pres[ $j ]" } ) ;
# Tropopause near end of TTCC if not found in TTAA group
if( $TTCC !~ / 88999/ && $TTCC =~ m# 88(\d\d)(\d) ([\d/]{5}) ([\d/]{5}) # ) {
$trop_pres = $1 . "." . $2 ;
$thermo = $3 ;
$wind = $4 ;
$trop_pres += 0.5 ;
$trop_pres =~ s#\.\d*## ;
# sets trop_t, trop_td_dep
trop_dec_thermo( $thermo ) ;
# sets trop_wdir, trop_wspd
trop_dec_wind( $wind ) ;
} # End of tropopause level in TTCC group
# Maximum Wind
if( $TTCC !~ / 77999/ && $TTCC =~ m# (66|77)(\d\d)(\d) ([\d/]{5}) # ) {
$max_pres = $2 . "." . $3 ;
max_dec_wind( $4 ) ;
$max_pres += 0.5 ;
$max_pres =~ s#\.\d*## ;
}# End of maximum wind level in TTCC group
} #end ttcc
# TTDD reports, significant levels with respect to temperature: under 100 mb
sub ttdd
{
my( $TTDD, $offset ) ;
( $TTDD ) = @_ ;
$offset = 0 ;
while( $offset < length( $TTDD )) {
last unless substr( $TTDD, $offset, 2 ) =~ /00|11|22|33|44|55|66|77|88/;
$pres[ $j ] = $ht[ $j ] = $t[ $j ] = $td[ $j ] = $wdir[ $j ] =
$wspd[ $j ] = $F ;
if( substr( $TTDD, $offset +2, 10 ) =~ /(\d{3}) (\d{5}) / ) {
$pres[ $j ] = sprintf( "%3.1f", $1 / 10 ) ;
dec_thermo( $2 ) ;
$j++ ;
}
$offset += 12 ;
}
} #end ttdd
# PP(BB|DD) reports, significant levels with respect to winds: all levels
sub ppbbdd
{
my( $PP, $offset, $i, $ten_thousand, $thousand ) ;
( $PP ) = @_ ;
$offset = 0 ;
while( $offset < length( $PP )) {
last unless substr( $PP, $offset, 1 ) eq "9" ;
$ten_thousand = substr( $PP, $offset +1, 1 ) ;
for( $i = 0; $i < 3; $i++ ) {
last if(( $thousand =
substr( $PP, $offset +2 + $i, 1 )) eq "/" ) ;
$pres[ $j ] = $ht[ $j ] = $t[ $j ] = $td[ $j ] = $wdir[ $j ] =
$wspd[ $j ] = $F ;
$ht[ $j ] = ( $ten_thousand . $thousand . "000" ) * 0.3048 ;
next unless $ht[ $j ] ;
if( substr( $PP, $offset +6 + $i * 6, 6 ) =~ /(\d{5}) / ) {
dec_wind( $1 ) ;
}
next if( $wdir[ $j ] eq $F || $wspd[ $j ] eq $F ) ;
$ht[ $j ] = sprintf( "%3.1f", $ht[ $j ] ) ;
$j++ ;
}
$offset += 6 + $i * 6 ;
}
} # end ppbbdd
sub dec_thermo
{
my( $thermo ) = @_ ;
my( $td_dep ) ;
$t[ $j ] = substr( $thermo, 0, 2 ) . "." . substr( $thermo, 2, 1 ) ;
$td_dep = substr( $thermo, 3, 2 ) ;
if( $t[ $j ] =~ m#/# ) {
$t[ $j ] = $F ;
} else {
$t[ $j ] =~ s/^0// ;
if( $t[ $j ] =~ m#[13579]$# ) { $t[ $j ] = "-" . $t[ $j ] ; }
}
if( $td_dep =~ m#/# ) {
$td[ $j ] = $F ;
} else {
# change to report Dew-point depression instead of Dew-point
if( $td_dep <= 50 ) {
#$td[ $j ] = $t[ $j ] - $td_dep / 10.0 ;
$td[ $j ] = $td_dep / 10.0 ;
} else {
#$td[ $j ] = $t[ $j ] - ( $td_dep -50.0 ) ;
$td[ $j ] = $td_dep -50.0 ;
}
$td[ $j ] *= 1.0 ;
}
# Convert to kelvin
$t[ $j ] = celsius2kelvin( $t[ $j ] ) ;
return ;
} # end dec_thermo
sub trop_dec_thermo
{
my( $thermo ) = @_ ;
my( $td_dep ) ;
$trop_t = substr( $thermo, 0, 2 ) . "." . substr( $thermo, 2, 1 ) ;
$td_dep = substr( $thermo, 3, 2 ) ;
if( $trop_t =~ m#/# ) {
$trop_t = $F ;
} else {
$trop_t =~ s/^0// ;
$trop_t = "-" . $trop_t if( $trop_t =~ m#[13579]$# ) ;
}
if( $td_dep =~ m#/# ) {
$trop_td_dep = $F ;
} else {
# change to report Dew-point depression instead of Dew-point
if( $td_dep <= 50 ) {
#$trop_td_dep = $trop_t - $td_dep/10.0 ;
$trop_td_dep = $td_dep/10.0 ;
} else {
#$trop_td_dep = $trop_t - ( $td_dep -50.0 ) ;
$trop_td_dep = $td_dep -50.0 ;
}
$trop_td_dep *= 1.0 ;
}
# Convert to kelvin
$trop_t = celsius2kelvin( $trop_t ) ;
} # end trop_dec_thermo
#
sub dec_wind
{
my( $wind ) = @_ ;
$wdir[ $j ] = substr( $wind, 0, 2 ) ;
$wspd[ $j ] = substr( $wind, 2, 3 ) ;
if( $wdir[ $j ] =~ m#/# ) {
$wdir[ $j ] = $wspd[ $j ] = $F ;
} else {
$wdir[ $j ] =~ s/^0// ;
$wdir[ $j ] .= "0" ;
if( $wspd[ $j ] > 500 ) {
$wdir[ $j ] += 5 ;
$wspd[ $j ] -= 500 ;
}
}
} # end dec_wind
#
sub trop_dec_wind
{
my( $wind ) = @_ ;
$trop_wdir = substr( $wind, 0, 2 ) ;
$trop_wspd = substr( $wind, 2, 3 ) ;
if( $trop_wdir =~ m#/# ) {
$trop_wdir = $trop_wspd = $F ;
} else {
$trop_wdir =~ s/^0// ;
$trop_wdir .= "0" ;
if( $trop_wspd > 500 ) {
$trop_wdir += 5 ;
$trop_wspd -= 500 ;
}
$trop_wspd *= 0.5144 if( $knots ) ;
}
} # end trop_dec_wind
# calculate max winds
sub max_dec_wind
{
my( $wind ) = @_ ;
$max_wdir = substr( $wind, 0, 2 ) ;
$max_wspd = substr( $wind, 2, 3 ) ;
if( $max_wdir =~ m#/# ) {
$max_wdir = $max_wspd = $F ;
} else {
$max_wdir =~ s/^0// ;
$max_wdir .= "0" ;
if( $max_wspd > 500 ) {
$max_wdir += 5 ;
$max_wspd -= 500 ;
}
$max_wspd *= 0.5144 if( $knots ) ;
}
} # end max_dec_wind
# returns kelvin given celsius temperature
sub celsius2kelvin
{
my( $celsius ) = @_ ;
if( $celsius == $F ) {
return $F ;
} else {
return $celsius + 273.15 ;
}
} # end celsius2kelvin
sub incrementRDAY
{
if( $themonth =~ /02/ ) {
if( $rday < 28 ) {
$rday++;
return;
}
} elsif( $themonth =~ /04|06|09|11/ ) {
if( $rday < 30 ) {
$rday++;
return;
}
} elsif( $themonth =~ /01|03|05|07|08|10|12/ ) {
if( $rday < 31 ) {
$rday++;
return;
}
}
# have to calculate first day of month
$thetime = timegm(0, 0, $rhour, $rday, $themonth -1, $thedecade, 0,0,0);
$thetime += 3600;
( $ss, $mm, $rhour, $rday, $themonth, $theyear, $wday, $yday, $isdst ) =
gmtime( $thetime ) ;
$theyear = ( $theyear < 100 ? $theyear : $theyear - 100 ) ;
$theyear = "20" . sprintf( "%02d", $theyear ) ;
$themonth++ ;
$themonth = sprintf( "%02d", $themonth ) ;
$rday = sprintf( "%02d", $rday ) ;
$rhour = sprintf( "%02d", $rhour ) ;
} # end incrementRDAY
sub theTime
{
my( $ss, $mm, $hh, $mday, $mon, $fmon, $year, $fyear, $wday, $yday, $isdst ) ;
my( $when ) = @_ ;
$mm = 0 ;
if( $when eq "synoptic" ) {
$mday = substr( $yyyymmddhh, 6, 2 ) ;
$hh = substr( $yyyymmddhh, 8, 2 ) ;
$fmon = substr( $yyyymmddhh, 4, 2 ) ;
$fyear = substr( $yyyymmddhh, 2, 2 ) ;
} elsif( $when eq "release" ) {
$mday = $rday ;
$hh = $rhour ;
$fmon = $themonth ;
$fyear = $thedecade ;
}
$time = timegm( 0, $mm, $hh, $mday, $fmon -1, $fyear, 0, 0, 0 ) ;
return $time ;
} # end theTime