<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.6000.16788" name=GENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=#ffffff>
<DIV><FONT face=Arial size=2>
<DIV><FONT face=Arial size=2>Hello all.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>I'm trying to plot ensemble ECMWF data for
wind.</FONT></DIV>
<DIV><FONT face=Arial size=2>It is not a much difficult job to plot the time
series of standard deviation and mean for intensities.</FONT></DIV>
<DIV><FONT face=Arial size=2>The problem comes with directions. Sorry, it could
be just a mistake of mine..... but I cannot identify it!</FONT></DIV>
<DIV><FONT face=Arial size=2>I'm working with version 1.8.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>I start with 51 grib files (one for each ensemble
member), where latitude, longitude, level and time vary.</FONT></DIV>
<DIV><FONT face=Arial size=2>I want to have a diagram for each grid point and
level, so I let only time to vary (in 40 steps) when I run the
plotting job.</FONT></DIV>
<DIV><FONT face=Arial size=2>An example of the resulting diagram is attached as
a png file.</FONT></DIV>
<DIV><FONT face=Arial size=2><STRONG>I can have a scatter plot for the
(u,v) couples, but it seems that I must not set time, even if such a
setting covers the whole time range. Otherwise, the plot is empty. That is, I
must avoid any 'set t .....' command in the part of the script ruling the
scatter plot, anyway I obtain too many points in the plot (maybe all
members x all times = 2040 points in the u-v
plane!!!!!).</STRONG></FONT></DIV>
<DIV><STRONG><FONT face=Arial size=2>How can I have one plot for each time step,
instead of a plot for the whole time range?</FONT></STRONG></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>This is one of the ctl files ("p2.ctl", the grib
file is attached):</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2>*************************************************************************************************</FONT></DIV>
<DIV><SPAN lang=IT>
<P>dset ^p2.grib</P>
<P>index ^p2.grib.idx</P>
<P>undef 9.999E+20</P>
<P>title p2.grib</P>
<P>* produced by grib2ctl v0.9.12.5p41</P>
<P>dtype grib 255</P>
<P>options yrev</P>
<P>ydef 4 linear 45.000000 0.5</P>
<P>xdef 6 linear 10.500000 0.500000</P>
<P>tdef 40 linear 18Z23feb2009 6hr</P>
<P>* z has 3 levels, for prs</P>
<P>zdef 3 levels</P>
<P>1000 925 850</P>
<P>vars 12</P>
<P>no10Usfc 0 165,1,0 ** surface 10 metre u wind component m s**-1</P>
<P>no10Vsfc 0 166,1,0 ** surface 10 metre v wind component m s**-1</P>
<P>no2Tsfc 0 167,1,0 ** surface 2 metre temperature K</P>
<P>SFsfc 0 144,1,0 ** surface Snow fall m(of water equivalent)</P>
<P>Tprs 3 130,100,0 ** (profile) Temperature K</P>
<P>TPsfc 0 228,1,0 ** surface Total precipitation m</P>
<P>U850mb 0 131,100,850 ** 850 mb U-velocity m s**-1</P>
<P>V850mb 0 132,100,850 ** 850 mb V-velocity m s**-1</P>
<P>Z850mb 0 129,100,850 ** 850 mb Geopotential (at the surface = orography) m**2
s**-2</P>
<P>var121sfc 0 121,1,0 ** surface undefined</P>
<P>var122sfc 0 122,1,0 ** surface undefined</P>
<P>var123sfc 0 123,1,0 ** surface undefined</P>
<P>ENDVARS</P>
<P>*************************************************************************+</P></SPAN></DIV>
<DIV><FONT face=Arial size=2>Then, here is the grads job to plot the
data:</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2>*********************************************************************************************************************<SPAN
lang=IT>
<P>function rose(args)</P>
<P>'reinit'</P>
<P>llon=subwrd(args,1)</P>
<P>llat=subwrd(args,2)</P>
<P>llev=subwrd(args,3)</P>
<P>*Scelgo il titolo con la localita' corrispondente alle coordinate</P>
<P>'open c.ctl'</P>
<P>'set lat 'llat</P>
<P>'set lon 'llon</P>
<P>nomepunto=''</P>
<P>coordinate=llat' 'llon</P>
<P>if (coordinate='46.00 12.00') ; nomepunto='FELTRE'; endif</P>
<P>if (coordinate='46.50 12.00') ; nomepunto='LIVINALLONGO'; endif</P>
<P>if (coordinate='46.50 12.25') ; nomepunto='MISURINA'; endif</P>
<P>if (coordinate='46.50 12.50') ; nomepunto='AURONZO'; endif</P>
<P>if (coordinate='46.25 12.00') ; nomepunto='AGORDO'; endif</P>
<P>if (coordinate='46.25 12.25') ; nomepunto='LONGARONE'; endif</P>
<P>if (coordinate='46.25 12.50') ; nomepunto='CIMOLAIS'; endif</P>
<P>if (coordinate='46.00 11.50') ; nomepunto='MARCESINA'; endif</P>
<P>if (coordinate='46.00 11.75') ; nomepunto='ARSIE'; endif</P>
<P>if (coordinate='46.00 12.00') ; nomepunto='FELTRE'; endif</P>
<P>if (coordinate='46.00 12.25') ; nomepunto='VITTORIO VENETO'; endif</P>
<P>if (coordinate='45.75 11.25') ; nomepunto='VALLI DEL PASUBIO'; endif</P>
<P>if (coordinate='45.75 11.50') ; nomepunto='CARRE VI'; endif</P>
<P>if (coordinate='45.75 11.75') ; nomepunto='BASSANO'; endif</P>
<P>if (coordinate='45.75 12.00') ; nomepunto='MONTEBELLUNA'; endif</P>
<P>if (coordinate='45.75 12.25') ; nomepunto='VILLORBA'; endif</P>
<P>if (coordinate='45.75 12.50') ; nomepunto='ODERZO'; endif</P>
<P>if (coordinate='45.75 12.75') ; nomepunto='PORTOGRUARO'; endif</P>
<P>if (coordinate='45.75 13.00') ; nomepunto='SAN MICHELE AL TAGLIAMENTO';
endif</P>
<P>if (coordinate='45.50 10.75') ; nomepunto='LAZISE'; endif</P>
<P>if (coordinate='45.50 11.00') ; nomepunto='VR TORRICELLE'; endif</P>
<P>if (coordinate='45.50 11.25') ; nomepunto='MONTECCHIA DI CROSARA'; endif</P>
<P>if (coordinate='45.50 11.50') ; nomepunto='ALTAVILLA VIC.NA'; endif</P>
<P>if (coordinate='45.50 11.75') ; nomepunto='GRISIGNANO DI ZOCCO'; endif</P>
<P>if (coordinate='45.50 12.00') ; nomepunto='PIANIGA'; endif</P>
<P>if (coordinate='45.50 12.25') ; nomepunto='MESTRE'; endif</P>
<P>if (coordinate='45.50 12.50') ; nomepunto='CAVALLINO'; endif</P>
<P>if (coordinate='45.50 12.75') ; nomepunto='LARGO DI JESOLO'; endif</P>
<P>if (coordinate='45.25 11.00') ; nomepunto='ISOLA DELLA SCALA'; endif</P>
<P>if (coordinate='45.25 11.25') ; nomepunto='ROVERCHIARA'; endif</P>
<P>if (coordinate='45.25 11.50') ; nomepunto='MONTAGNANA'; endif</P>
<P>if (coordinate='45.25 11.75') ; nomepunto='MONSELICE'; endif</P>
<P>if (coordinate='45.25 12.00') ; nomepunto='BOVOLENTA'; endif</P>
<P>if (coordinate='45.25 12.25') ; nomepunto='LAGUNA SUD'; endif</P>
<P>if (coordinate='45.00 11.25') ; nomepunto='CARBONARA DI PO'; endif</P>
<P>if (coordinate='45.00 11.50') ; nomepunto='BAGNOLO DI PO'; endif</P>
<P>if (coordinate='45.00 11.75') ; nomepunto='ARQUA POLESINE'; endif</P>
<P>if (coordinate='45.00 12.00') ; nomepunto='ADRIA'; endif</P>
<P>if (coordinate='45.00 12.25') ; nomepunto='PORTO VIRO'; endif</P>
<P>if (coordinate='45.00 12.50') ; nomepunto='PILA'; endif</P>
<P>if (coordinate='45.75 10.75') ; nomepunto='GARDA NORD'; endif</P>
<P>if (coordinate='45.75 11.00') ; nomepunto='ALA (TN)'; endif</P>
<P>*Epsigramma in basso, con la statistica del modulo del vento</P>
<P>*'set parea 1 8 0.5 5'</P>
<P>'set parea 1 11 0.5 3.75'</P>
<P>*Definisco le variabili da plottare nel diagramma in basso:</P>
<P>*media e varianza del modulo del vento</P>
<P>'set lev 'llev</P>
<P>'set t 1 40'</P>
<P>if llev=1000</P>
<P>'define campoc=sqrt(no10Usfc*no10Usfc+no10Vsfc*no10Vsfc)'</P>
<P>'define campocu=no10Usfc'</P>
<P>'define campocv=no10Vsfc'</P>
<P>else</P>
<P>'define campoc=sqrt(U850mb*U850mb+V850mb*V850mb)'</P>
<P>'define campocu=U850mb'</P>
<P>'define campocv=v850mb'</P>
<P>endif</P>
<P>i=1</P>
<P>'define ps=0.0'</P>
<P>'define pv=0.0'</P>
<P>while(i<51)</P>
<P>'open p'i'.ctl'</P>
<P>'set dfile 2'</P>
<P>'set lat 'llat</P>
<P>'set lon 'llon</P>
<P>'set lev 'llev</P>
<P>'set t 1 40'</P>
<P>if llev=1000</P>
<P>'define campop'i'=sqrt(no10Usfc*no10Usfc+no10Vsfc*no10Vsfc)'</P>
<P>'define campop'i'u=no10Usfc'</P>
<P>'define campop'i'v=no10Vsfc'</P>
<P>else</P>
<P>'define campop'i'=sqrt(U850mb*U850mb+V850mb*V850mb)'</P>
<P>'define campop'i'u=U850mb'</P>
<P>'define campop'i'v=V850mb'</P>
<P>endif</P>
<P>'define ps=ps+campop'i</P>
<P>'close 2'</P>
<P>i=i+1</P>
<P>endwhile</P>
<P>'define ps=ps+'campoc</P>
<P>'define pm=ps/51.0'</P>
<P>*Plot delle medie nel diagramma in basso</P>
<P>'set grads off'</P>
<P>'set gxout line'</P>
<P>i=1</P>
<P>while(i<51)</P>
<P>'set vrange 0 20'</P>
<P>'set ccolor 1'</P>
<P>'set cmark 3'</P>
<P>'set cthick 1'</P>
<P>'set csmooth on'</P>
<P>'set missconn on'</P>
<P>'d campop'i</P>
<P>i=i+1</P>
<P>endwhile</P>
<P>i=1</P>
<P>while(i<51)</P>
<P>'define pvd=pm-campop'i</P>
<P>'define pv=pv+pvd*pvd'</P>
<P>i=i+1</P>
<P>endwhile</P>
<P>'define pvd=pm-'campoc</P>
<P>'define pv=pv+pvd*pvd'</P>
<P>'define pstd=sqrt(pv/50.0)'</P>
<P>'set ccolor 2'</P>
<P>'set cmark 3'</P>
<P>'set cthick 3'</P>
<P>'set csmooth on'</P>
<P>'set missconn on'</P>
<P>'set vrange 0 20'</P>
<P>'d campoc'</P>
<P>'set ccolor 4'</P>
<P>'set cmark 1'</P>
<P>'set cthick 3'</P>
<P>'set cstyle 0'</P>
<P>'set missconn off'</P>
<P>'d pm'</P>
<P>*Barre di deviazione standard nel diagramma in basso</P>
<P>'set gxout bar'</P>
<P>'set bargap 50'</P>
<P>'set baropts outline'</P>
<P>'set barbase 0'</P>
<P>'set ccolor 4'</P>
<P>'set cthick 10'</P>
<P>'define pp1=pm-pstd/2.0'</P>
<P>'define pp2=pm+pstd/2.0'</P>
<P>'set vrange 0 20'</P>
<P>'d pp1;pp2'</P>
<P>*Diagrammi in alto: scatter plots con le punte dei vettori vento</P>
<P>*'set parea 1 8 5.5 10.5'</P>
<P>'set parea 1 11 4.25 7.5'</P>
<P>*Plot nel diagramma in alto</P>
<P>'set grads off'</P>
<P>'set grid off'</P>
<P>'set xlab off'</P>
<P>'set ylab off'</P>
<P>'set gxout scatter'</P>
<P>'set vrange -20 20'</P>
<P>'set vrange2 -20 20'</P>
<P>i=1</P>
<P>while(i<51)</P>
<P>'set ccolor 1'</P>
<P>'set digsize 0.05'</P>
<P>'set cmark 3'</P>
<P>'d campop'i'u;campop'i'v'</P>
<P>i=i+1</P>
<P>endwhile</P>
<P>'set ccolor 2'</P>
<P>'set digsize 0.05'</P>
<P>'set cmark 3'</P>
<P>'d campocu;campocv'</P>
<P>'draw title 'nomepunto</P>
<P>'enable print rose.gmf'</P>
<P>'print'</P>
<P>'disable print'</P>
<P>'quit'</P>
<P>*******************************************************************************************</P>
<P><EM>Antonino Claudio Bonan</EM></P>
<P><EM>Meteorologo</EM></P>
<P><EM>Centro Meteorologico di Teolo - ARPAV</EM></P>
<P><EM>via Marconi 55</EM></P>
<P><EM>35037 Teolo (PD) -
Italy</EM></P></SPAN></FONT></DIV></FONT></DIV><p></p><p><a href="http://www.spambrave.com/">I'm protected by SpamBrave</a></p></BODY></HTML>