TRSTAPT.c = Trajectory STAtion by Particle Travel times coordinates 3/07/10
1. STAPTTCT.c<--->PTTCT.c = trSTApt.c
Due to the new shear explained previously and the new RUC predicted soundings, we have modified
PTTC.c to call it PTTCT.c = Particle Travel Times Coordinates Trajectories .
The sum of the shear coordinates (tilt) is eliminated and the highest heights of these predicted soundings are about 74,000',
a little less than the earlier Storm Machine web site soundings. This should suffice for the Nuclear Cloud tops but this will
be determined after many runs. Coupling PTTCT.c
into this main program is now the subroutine pttct() . We have the back and forth
interchange between these two explained below.
The purpose of TRSTAPT.c is to calculate the horizontal vector distance a standard particle sp travels using the latest pttct.mb[m] file created by the subroutine pttct() using the latest RUC predicted soundings. Each trdis > 48 km or tth > 1 hr traveled calls for a new RUC file (mcs.txt) from a new METAR's STA and/or time for pttct() to convert to a usable sounding and send to pttct.mb[m] for the main() TRSTAPT.c to repeat the cycle again. However, we only use the original RUC (000.txt) for the ground zero STA and list STA name & time in ruco.sta. Otherwise, we use the cde,cdn and ctth for each pttct.mb[m]'s.
Simply, we start to descend particles by their fall velocities. If the particles pass closer to another METAR'S STA coordinate that is closer than the original METAR's STA but > 48 km, we pick up on this new RUC predicted sounding for that grid station and update the pttct.mb[m] file continually until all 5 representative particles of discrete height reach the surface. Thus, we have a pttct.mb[m] file for each of the m=6 representative particles. Besides a particle travel > 48 km, we also use a time criteria. If the accumulated travel time (tth) > 1 hr, we also use a new RUC sounding. With higher winds with altitudes, these representative particles cause the last pttct.mb[m] to become a composite of more STA locations with different times. We have to request these new soundings from the RUC web site using subroutine rucreq(). To date, this is a semi-automatic procedure which requires some interfacing with cut and paste editor commands. We created a file called 'ruc.dir' that contains mcs.txt - STA,yr,mo,day & Zulu hr request name files seen below. We link these to reuse some mcs.txt files again, if possible. The combining of these RUC soundings into the various pttct.* files is displayed and explained best in trSTApt.dia & trSTAptL.dia.
There is also some adjustments to layer travel time tt because of possible different thicknesses of layer levl between the former and later predicted RUC soundings. This is seen an explained better in DELtt.dia. A list of all these stations is compiled in 3-let US METAR Stations. We made this list into a Folder of Placemarks in *.kml format for Google Earth again and used their screen captured tool as seen in US METAR Stations.
B. main() TRSTAPT.c = Trajectory STAtion Particle Travel time :
We find out the computer's date & time by _strdate(date); _strtime(time); to submit to the
RUC website for the time of the initial requested predicted sounding. We read city.dat for
the city's name that is closest to a METAR station in file
3-let US METAR Stations
usually the same as the city's name. We set combined soundings cs=0 and # of levels levl=61 then send to RUC
request rucreq() copying the min STA distance let3m = name. We obtain the initial RUC format
file from the RUC web site .
The # of predicted hours is 1.0 since we only want the first one - which is our default. The time is a few hours after the 00z or 12z sounding or the first forecasted sounding available for the METAR's sta ( 3-let US METAR Stations) alphabetical listing.
(1) RUC Request by html:
We now set up the request for this sounding by subroutine rucreq(). When requesting the original RUC
sounding manually, the address showed in the address window. We highlighted and copied this request address string,
it revealed to be about 285 characters long. This was saved as a generic htm format file
(staTRsta.htm) then copied to a new *.htm file called RUCreq.htm with the new appropriate
fields (hr, 5-min, time, month year & sta) replaced into the generic address string. The
needed section of this 285 character is given below:
(2) Run html RUC Request:
One method is to use Windows Explorer in My Computer, then find RUCreq.htm and just click. Keep this window open for the next request put into RUCreq.htm. The other and best method is to have another DOS window open and type the command line RUCreq.htm. The system(" "); command in C didn't work like it does for *.bat & *.exe. When RUCreq.htm is opened locally off line in Internet Explorer or run by DOS command, when we click on the displayed associated address tag Station Request , the sounding appears as a text file on the screen.
(3) Saving Requested RUC File:
Using MS Notepad, menu is top bar Edit > Select All then Copy. If we use New then pasted, it can
be saved as an older 000.txt. The paste may have more lines below 68 but we only read 62 lines
of data past the headers. The important thing is to copy the last sounding in the list, usually just one,
because that matches the expanded layer travel time tth request time. The main() in TRSTAPT.c becomes suspended
with a 3 line message:
RUC Request for PHL. Run RUCreq.htm from another DOS Window
Click on Station Request
Save as 000.txt in Radfo\PHL
(4) RUC Request by PERL (*.pl):
Returning to RUC web site again, explaining the importance of this program,
soon to retire was a Bill Moninger who informed me about a automate access program that many use. A sample perl script file
was e-mailed to me and have only to replace values. Simply station request and times. However, the time values are in seconds,
rounding off the nearest hour from # of secs from 1/1/1970 - the birth of baby UNIX. This does eliminate year, month, day and hour.
After doing a request manually, the seconds converted in the address queue for the 1st day in April 2010 = 1270080000 secs.
From there, all future requests are calculated. The smoothed generic file used for this is seen in
The latest Window's perl version is 188.8.131.52 as seen by their web site Strawberry Perl . After learning some basics, it is incorporated in this program and run as a system("perl *.pl"); command. What used to take 20 minutes, is now done within 20 seconds automatically and now saved as 000.ruc, etc. in the city name directory (NYC,PHL,CHI). End Update 04/29/10
C. Subroutine pttct()- or PTTCT.c = Particle Travel Times Coordinates Trajectories:
Now we run pttct() with # of combined METAR stations (cs) and last LEVeL (levl) of the last METAR station.
For our initial run, cs=0 levl=61 explained below. Reading the RUC's PHL\000.txt, we extract date-time
and METAR station name (let3), and create our detailed file name 'fon' with the 3 letter suffix is the
station PHL. This created file also goes into the city's directory as shown here as IN: = PHL\000.txt,
OUT: = PHL\09102814.PHL. We used the yr,mo,day & Zulu time and .STA to identify the converted sounding.
The weather parameters are converted to real values and mean layer parameters also in directory PHL as
PHL\09102814.PHL that can now be easily sequenced with the other converts. The # of wind levels is
ln=0->61 levels = 62 or 1->61 layers = 61 (ln). This initial file name must be listed in ruco.sta to be
run in RUCDHM.c for initial blast plume rise and nuclear cap & stem
cloud dimensions (cap.dh) used in the local city dust inflow calculations
INFLOW.c . Then finally, PHL\09102814.PHL will be used somewhat in the
main grid program RADFORUG.c .
A new file (shortened pttc.cor) is now called pttct.mb[m], where the t. means trajectory coordinates. Here, we reverse reading the height level order from up @ pttc.cor to down because before trajectories, we assumed a uniform sounding everywhere and kept the same METAR's station. Now, particles falling from near the top of the sounding pass through the first upper level (levl=61) which also identifies the first layer with mean weather & wind vectors for this 1st layer 61->60. We rewrite pttct.mb[m] down to our last layer, the surface l=1 (1->0) with new discrete layers for each trdis > 48 km or tth > 1). This is best displayed and explained in trSTApt.dia & trSTAptL.dia.
(1) Standard Particle Sizes:
We need to track a dust particle that flow past METAR'S stations. Based on this trajectory of a particle
passing within and closest to another downwind METAR's STA, we need a representative particle size sp[m]
and standard pressure level mb[m] for tracing. The series of bulk tons of dust concentration from nf1.NYC
data determined these standard pressure millibars (mb) levels-heights as seen below:
mb=1000 @ levl= 3 = 102m = 334.6' sp=220u pttct.1000 = pttct.0
mb= 700 @ levl=16 = 2837m = 9307.6' sp=210u pttct.700
mb= 500 @ levl=20 = 5279m = 17319.3' sp=200u pttct.500
mb= 300 @ levl=38 = 8764m = 28752.9' sp=160u pttct.300
mb= 150 @ levl=53 = 13310m = 43667.4' sp= 85u pttct.150
mb= 50 @ levl=60 = 20214m = 66318.1' sp= 60u pttct.50
The levl's may change slightly but it's heights (m & ') will change somewhat with each new RUC sounding. The pressure mb and particle size sp remain constant. All the RUC level counts are the same, I hope, are fixed at 62. It was best to choose a mandatory pressure level near these heights as seen above. In pttct(), we find the new heights each time for mb with new RUC's. However, since they don't change much and we also need these mb 6 heights' to match the hc heights' for the main 1.5 hr+ run final program RADFORUG.c , we will use only the 1st converted sounding from (000.txt) which is listed in ruco.sta as PHL\09102814.PHL. For PPTTCKML.c , we convert all the 62 pressure levl's to one of the m=0->5 standard pressure levl's mb shown in the left column by rounding off. This string of 62 l-m's is put to file levl.mb. This way, we can link to pttct.mb[m] by each level lev without using to evaluate all the levels in PHL\09102814.PHL again.
D. main() TRSTAPT.c = Trajectory STAtion Particle Travel time :
(1) Nearest Station Calculations:
Returning back to the main program, we still have cs=0, but now start with mb=1000, levl=3. With the initial
pttct.000 file written from pttct(), the downwind particle distance trdis, using cde,cdn -> de,dn conversion
to meters for check against a distance greater than 48 km. We keep increasing the layer's de,dn by adding
these next layer's components below it until trdis > 48 km or tth > 1hr, where we then update cdem,cdnm (the cumulated
travel components from the city's name STA). We use cdem,cdnm with the city.dat's utme,utmn and convert this to Long,Lat
by utmll3(). Keeping track of cde,cdn and cdem,cdnm with de,dn is best seen and explained in
trSTApt.dia & trSTAptL.dia again.
The RUC grid is 40 km, and a diagonal would be 56.5 km so we use 48 km. We need the nearest METAR's STA lat & lon nearest the grid point, as we can't request the RUC for the grid point itself. We search 3-let US METAR Stations for the STA's decimal lon/lat (dlon,dlat) and convert it to meters. Since we don't have to be exact, just relative by searching for the closest distance dis, we use a simpler method. The earth's radius R is 6378137 m perpendicular from the polar axis @ the equator (Lat 0), but at particle location, Lat is r=R*cos(Lat). So, arc dis Long (x) is r*(dlon-Long) in meters if dlon-Long are in radians. The S-N distance (y) is METAR'S STA Lat - Lat in deg = dlat-Lat. With 60 n-mi/deg, we now also have linear distance for Lat (y) differences from dlat. After conversion to meters, we have dis=sqrt(x^2 + y^2). We find this minimum dis and now download a new RUC profile for the next METAR's STAtion (let3).
(2) Saving Time & Pressure Level:
Along with the distance, we also extract tth - the expanded layer's travel time in hours. This enables us to
request the time of this new sounding for the METAR's station. The decimal hours is converted to hours
and 5 min intervals. To date, 5 min intervals don't make a difference in the sounding data.
We now have cs = 1 at the level where the particle trdis > 48 km or tth > 1 hr, say levl=47 for pttct() to
use again for overwriting pttct.000. With l=levl=47, we also save this last height
of the levl h[lev] as hl1 for travel time adjustment between the 2 pair of layers later.
(3) Saving More Time Requests:
Since upper air features are usually more pronounced in the lower levels, the higher levels are usually more streamlined. Therefore, we increase the hour interval and the 48 km distance by a factor with height. With the 6 standard particle levels above (mb[m]), m=0->5. For each m, we now have a factor (m+1)/2, which is also the new tested hour travel time criteria. Successively, we have the factor-hours 1,1,2,2,3,3 for m=0->5. This was a real help in reducing the # of requests especially for the smaller particles at the higher levels without sacrificing much accuracy. End Update 04/29/10
E. Subroutine pttct() Again:
(1) Using Saved Parameters:
With cs= 1 & levl= 47, we can now open the next RUC file of mcs.txt as PHL\001.txt just created in RUCreq.htm. We open PHL\001.txt and create the corrected usable weather layers file or IN: = PHL\001.txt, OUT: = PHL\09102815.PHL. The 01 in 001.txt is the cs term. The final cs is the # of combined stations used to rewrite pttct.mb[m] files. However, we only overwrite a new section (1-levl) layers of each file for each cs for cde,cdn,ctth. Layers above levl remain untouched until a higher mb[m] level is started. This is very evident and explained more in trSTApt.dia & trSTAptL.dia.
(2) Different Layer Thicknesses:
There is a slight problem. Say with our 1st levl=53, the new sounding thickness of this layer in 001.txt
might be different from that in the previous 000.txt. This would change the expanded layer travel time tth when
the particles would reach the surface. The interface layer's tt adjustment must be made is why we saved the top of the
1st sounding @ levl as hl1 in the main() of TRSTAPT.c . We have the present tt @ levl as tt=dz/fps=(h[l]-h[l-1])/fps
after converting h's to feet by fps='/sec for particle sp. Now we have a new tt where we adjust dz
because of possible different thicknesses of layer levl between these 2 soundings. This is best explained
and seen in diagram DELTT.dia.
These corrections are made before we write over the pttct.mb[m] file. With this new profile below levl, we have new accumulated cde,cdn & ccth terms to write. We find file position filp=(id-1)*ln*6+l*6 or reading the 1st l, then rewrite these new cde,cdn & ccth terms for each increasing l level for all np particle sizes, not just the standard sp u sizes used in TRSTAPT.c . The main program, ( RADFORUG.c ) needs all these levels & particles vectors and ccth stored in the pttct.mb[m] files.
F. main() TRSTAPT.c = Trajectory STAtion Particle Travel time:
We return to main() with cs= 1 but still have levl= 47. We continue down the l level loop we left from with the rewritten new section of pttct.000 file, down to where trdis > than another 48 km and create the next RUCreq.htm request file for the next RUC sounding station. We now read pttct.000 in subroutine rcortt() for the last level's cumulative terms cde,cdn & ccth. Since this is now where we continue from, we add the new travel time and components from this new initial level levl=47. Therefore, we use these @ levl=47 as cde0,cdn0 & ctth0 as the initial levl for the continuing descent. So then all new vector components and travel times are added from these held new initial values for the next de,dn expansion to calculate distances trdis.
(1) Reuse Previous RUC Sounding in Subroutine rucreq()?:
If the search for the minimum STA in 3-let US METAR Stations is still the same station and the
tth is < 1 hr, then we keep the same sounding and continue down the l level loop. We don't need a new IN: = PHL\0??.txt, as
the last one is still good. How could this happen? Consider a particle traveling with a west wind out to sea past ACY
or HAT. No METARS stations out there, unless automated buoys are considered as METARS STA's. A message is
written to the screen and bug file as OUT: = PHL\09102814.PHL with hr = ..15., ..16. ..17.,etc. These files are
not really needed any more as we used the creation of them to only rewrite the pttct.000 file. However, we still
store them in the \PHL directory but mostly with the new station name like PHL\09102816.WRI.
A problem occurred when we wanted to know this output STA - Date OUT: file before we read it's header. This would mean we downloaded it but we really want to determine if we can reuse an older file before we download. Since we already have the request parameters in the staTRsta.htm address, we can use them as an id. Also, we can use a string greater than a 12345678.STA = 11 chars + a '.'. This string becomes the stored name in a file called 'ruc.dir'. This will link the associated mcs.txt with the STA yr,mo,day & Zulu hr file. Such a sample section is seen below:
We updated the fon name with a format like this seen below:
Besides winds out to sea, suppose we have little directional wind shear aloft. Also, if coupled with this, we have light winds such as during summer. What appears is that the particles from higher mb levels will take nearly the same direction and travel time as the lower mb levels. This indicates that we may be required to request a previous sounding. One of the semi-manual complexities to reduce is the number of requests and collection of these mcs.txt RUC STA's soundings. We can only match date time - STA or the 'fon' file name. So, we read left & right file names and if 'fon' matches our new request fon's yr,mo,day & Zulu hr .STA file, we copy the already used left mcs.txt file name to the new requested one and ignore the request routine. The screen and bug file writes this out as seen in the example below:
IN: PHL\402.txt OUT: PHL\10010903.A04
Previously Used OUT: PHL\10010903.A04 - No New Request Needed
copy PHL\200.txt PHL\402.txt
This system("copy PHL\200.txt PHL\402.txt "); command is created and executed within the program.
G. File Checks:
(1) File Limits:
The complexity of this program, just to prepare the weather input, probably rivals any of these other programs. Besides the subroutines rcortt(int,int),cortt(),utmll3(),rucreq(),reqnam(),reuse(),nameout(),pttct() and midlayer(), we skip in and out of these usually from main() while some are nested within each other. One problem is the number of files allowed to be opened at one time. The compiler requires 5, and the limit is 15 and we initially had more than 10. The first remedy is to open and close the files if they are only used once by the file handle (ft) denoting temporary. Once finished reading, the same file handle (ft) can be used again by another file.
(2) File Open - Closed Balance:
Skipping around this program is dangerous as it requires that a particular file must be opened once and left opened or closed
before opening it again. Conversely, closing the same file twice can also be devastating to reading or writing the other files.
In the bug program trSTApt.bug, many steps and data can be seen written out. If we edit this file
to show only the closed - open handles, then we can do an easier check. Starting with a shorter cycle run, we have some file
count handle errors and these effects can be seen with explanations below:
(3) File Open - Closed Balance Errors:
H. RUC Data Message:
We have some new problems with the predicted RUC profiles. We correct these in subroutine pttct().
A. Dew points prediction > temperature or supersaturation. Here we just simply set dew point = temperature.
B. Predicted significant levels too close to the mandatory levels. Even though the pressure is reported to the nearest .1 mb, a case was found that the height of a significant level was within 1 meter from the height of the mandatory 1000.0 mb level. Close, but we could consider that 1 m thickness, although absurd, must keep the total # of layers = 61. However, the calculated pressure rounded of to the nearest .1 mb was still 1000.0 mb. The problem is now that we use pressures difference between levels to calculate the density which would give us density = 0, or a vacuum. This would play havoc with the fall velocity formulas for particles. For the insignificant 1 m layer, or nearly so, we use the density of the layer immediately above.
C. Another case was found that a predicted significant level was the 100 mb level and the height of 16056m were the same, so obviously were the temps & winds. Worse, was another case near the 200 mb level where the pressure was the same but the higher level's height was lower than the level below by 1 m. Originally thought the best way to handle these cases for any layer is to skip this level and reduce the ln=61, to 60.
D. Another case was found that a predicted significant level pressure of 1000.3mb was different than the next level of 1000.0 mb but their associated height's were 212 & 213 m msl. With a layer 1 m thick, the density was not representative. Also, the travel time tt to the nearest .001 hr can become less than the minimum .001 hr for this layer with larger particle diameters. So a ctth=0 hr is possible at the surface layer and would cause havoc with the rest of the program and is not representative of this layer. In the same sounding, 925.3 mb @ 844 m and the next level was 925.0 mb @ 846 m, a 2 m layer thickness. This caused a 31% jump in the layer's density compared the layer above.
A problem with skipping levels in C. above, is that altering the ln=61 levels can lead to other problems when overlapping different soundings in & out at the interface level l=levl. With these cases, it was best to combine B. & C. above with D. here for a new procedure. We take the average of the pressures below and above this too close level n. So, if the layer <= 4 m, we apply this averaging. We average pressure rather height so to split the weight of the air rather than the height. This makes little difference in the lower levels where it is most likely to occur. We then determine dz, explained below, for the average pressure, then add dz to the lower level. For our 2nd case, p= 925.3 -> 920.3 -> 915.3 where 920.3 replaces the 925.0 mb which causes the layer thickness increase from 2 -> 43.6 m.
For temp & dew point, the average of the n-1 & n+1 levels were also used, just as was done with the pressure levels. However, for the wind levels, we proportioned dz the wind @ n-1 -> n+1 layer thickness for velocity and direction. For direction, the wind shear must be corrected for the 360 degree crossover twice so the shear is no more or less than +,- 180 degrees. This easy case for direction was:
dir=195.0 -> dir=193.0 -> dir=191.0
E. The surface pressure is below a mandatory level. With a deep low pressure area or higher terrain, this is possible. First discovered in the PIT sounding given below:
PIT 12 kt
9709 324 134 124 247 4
10000 99999 99999 99999 99999 99999
9685 345 134 122 240 9
By convention, they report the surface data first here at 970.9 mb at 324 m msl. Therefore, the 1000.0 mb level is above that level or out of order. The first correction step is to reverse the order so the 1000.0 mb level is the first. Next, we fill the missing data (99999) by using the same dew point, wind dir & speed as the 970.9 mb level. What about the height and temperature? The same temp wouldn't be bad if we use the same height, but this creates problems later with a 0' thickness layer.
The first thing we do is use the hydrostatic equation dz=Rg*vt*log(p[n-1]/p[n]). Constant Rg=R/Ga where R=universal gas constant and Ga=acceleration by gravity. Since we don't have the air temperature below ground, we don't have a mean virtual temperature vt. We use the surface temperature for the mean vt which is only a slight error. Here we calculated dz=247.7 m. So the height of 1000.0 mb is 247.7 m below ground level (324 m msl) - 247.7m = 76.3m. Now for air temp, we could use the lapse rate of the layer above and extrapolate down. However, lapse rates near the surfaces can be highly super adiabatic (my thesis) and to extend it below the surface could be quite non-representative of this layer. We simply extend down (76.3m) the average standard lapse rate of 3.3 deg F/1000', adjusted to deg C/m gives us 14.89 deg C @ 76.3 m. The adjusted first 3 layers is then given below:
LVL Pres H(m) Deg C Avg C/100m DP C Avg Den(g/cm^3) Vis Dir Vel
2 968.5 345 13.4 13.40 0.00 12.2 12.30 1.165e-003 1.749e-004 240 10.3
1 970.9 324 13.4 14.14 -0.60 12.4 12.40 1.198e-003 1.752e-004 247 4.6
0 1000.0 76 14.9 247 4.6
I. Condensed Program Summary:
In order to comprehend this program more clearly, it was necessary to condense the program into a summary. The skipping from main() to subroutines and back to main(), transferring in & out of routines and tracking the open/close files problem made this imperative. This is all seen in trSTApt.fld.
2. Return Back To Radiosonde Soundings: