[gradsusr] How to combine streamline wind and precipitation into one file and make the code less messy?

Tue Jan 28 20:34:59 EST 2014

Hi, so what I am trying to do here is to plot Streamline plots and
Precipitation plots on the same plot...but when I try to create a loop for
that I get lots of error messages that it doesn't recognize my variable
'precip' that's with my precip netcdf file.  So, then I tried appending the
precip file onto the end and just ended up with even messier code.

Is there a lot shorter, easier way I could do this code and get both
precipitation and streamline winds on the same plot.  Note that my years in
this code come from text file 'NovJanPosQBO.txt'  I also include that code
at the end.

'reinit'
'set display color white';'c'
'set grid off'

*- - -
 xl=0.4 ; yt=7.5
 dx=2.3 ; dy=2.4
*- - -

size.1=10
size.2=42

file='NovJanPosQBO.txt'
ret=read(file)

n=1;while(n<=100)
ret=read(file)
chk=sublin(ret,1);if(chk>0);break;endif
ret=sublin(ret,2)

ip=1;while(ip<=2)
  if(n<=size.ip)
*    year.ip.n=substr(ret,ip)
     year.1.n=substr(ret,1,4)
     year.2.n=substr(ret,6,4)
  endif
ip=ip+1;endwhile

n=n+1;endwhile

*= = = = = = = = =

'sdfopen /export1/Databank/NCEP_Month/pressure/uwnd.mon.mean.nc'
'sdfopen /export1/Databank/NCEP_Month/pressure/vwnd.mon.mean.nc'

 lv=1 ; while(lv<=2)
'set t 17 775'
 if(lv=1); 'set lev 30' ;endif
 if(lv=2); 'set lev 200' ;endif
*if(lv=3); 'set lev 700' ;endif
*if(lv=4); 'set lev 1000' ;endif

'define pu1=ave(uwnd.1,t-3,t+3)'
'define pu2=ave(uwnd.1,t-16,t+16)'
'define pufilter=pu1-pu2'
'define pv1=ave(vwnd.2,t-3,t+3)'
'define pv2=ave(vwnd.2,t-16,t+16)'
'define pvfilter=pv1-pv2'

 ip=1 ; while(ip<=2)

'define ucomp'ip'=const(pufilter,0)'
'define vcomp'ip'=const(pvfilter,0)'

 n=1 ; while(n<=size.ip)

   'set time apr'year.ip.n
   'ucomp'ip'=ucomp'ip'+ave(pufilter,t-0,t+2)/'size.ip
   'vcomp'ip'=vcomp'ip'+ave(pvfilter,t-0,t+2)/'size.ip

 n=n+1 ; endwhile

 ip=ip+1 ; endwhile
* = = = = = = = = =
*'set gxout shaded'
'set gxout vector'
'umean=(ucomp1 + ucomp2)/2'
'vmean=(vcomp1 + vcomp2)/2'

 ip=1 ; while(ip<=2)

'define utot'1'= ucomp'ip' - umean'
'define vtot'1'= vcomp'ip' - vmean'
ip=ip+1 ; endwhile

* plotting

 ip=1 ; while(ip<=1)

 xr=xl+dx
 yb=yt-dy
'set parea 'xl' 'xr' 'yb' 'yt
'set lat -90 90'
'set lon 0 360'
'set lev 30'

'set mproj scaled'
'set mpdset mres'
'set ccolor 1'
'set map 11'
'set annot 1 2'
'set xlopts 15 1 0.08'
'set ylopts 15 1 0.08'

len=0.1
scale=3

 'color -10 10 -kind green->white->orangered'
 'd utot'1

*'set gxout vector'
*'set arrscl 'len' 'scale
*'set arrlab off'
*'set arrowhead -0.2'
*'d skip(utot'ip',6);skip(vtot'ip',6)'

 'set gxout stream'
'set cthick 1' ; 'set ccolor 1'
 'set strmden 1'
 'd utot'1';vtot'1

'set string 1 c 5'
'set strsiz 0.1 0.1'
'draw string 1.5 8.1 (1) Apr-Jun 30mb'
'draw string 1.5 7.9 Wind Anomaly for'
'draw string 1.5 7.7 QBO Positive Phase in Nov-Jan'

'set string 1 c 5'
'set strsiz 0.1 0.1'
'draw string 1.5 3.7 (1) Apr-Jun 200mb'
'draw string 1.5 3.5 Wind Anomaly for'
'draw string 1.5 3.3 QBO Positive Phase in Nov-Jan'

 xl=xr+0.4
 ip=ip+1 ; endwhile

 xl=0.4
 yt=yb-2.0
 lv=lv+1 ; endwhile
**********************QBO Positive Phase************************************

x0=5.5;y0=8.5
xl=x0; yt=y0
dx=3.0;dy=2.0

yb=yt-dy;xr=xl+dx

'sdfopen /export1/Databank/PRECIP/precl.mon.mean.1x1.nc'
*'sdfopen /export1/Databank/PRECIP/precip.1x1.mon.mean.nc'
'set t 17 770'
'define p1=ave(precip,t-3,t+3)'
'define p2=ave(precip,t-16,t+16)'
'define pfilter=p1-p2'
yy.1=1964
yy.2=1967
yy.3=1976
yy.4=1983
yy.5=1986
yy.6=1988
yy.7=1991
yy.8=1993
yy.9=2009
yy.10=2011

month.1='apr'
month.2='may'
month.3='jun'

yr=1; while(yr<=10)

mm=1; while(mm<=3)
if(yr=1);'define sea'mm'= const(pfilter,0)';endif
year=yy.yr;mon=month.mm
*say mm'/'year
'define sea'mm'= sea'mm'+ pfilter(time= 'mon''year')/10'
mm=mm+1;endwhile
yr=yr+1;endwhile

'set parea 'xl' 'xr' 'yb' 'yt
*say 'parea 'xl' 'xr' 'yb' 'yt
'define SeaTotal=(sea1+sea2+sea3)/3'

'set gxout shaded'
'color -0.3 0.3 0.03 -kind
maroon->red->orangered->orange->yellow->white->white->lime->forestgreen->blue->midnightblue->darkviolet'

'set lat 20 50'
'set lon 190 300'
'set t 17'
'define avera1= ave(pfilter,time=01apr1950,time=01apr2011,12)/3'
'define avera2= ave(pfilter,time=01may1950,time=01may2011,12)/3'
'define avera3= ave(pfilter,time=01jun1950,time=01jun2011,12)/3'

'define TotalAvera= avera1+avera2+avera3'
'set mpdset hires'
'd (Seatotal-TotalAvera)'
'set string 1 tc 5'
'set strsiz 0.1 0.1'
'draw string 7.0 8.2 Apr-Jun Precipitation Anomaly'
'draw string 7.0 8.0 for QBO Nov-Jan Positive Years'
'set csmooth on'
cbar(xr+0.14,yt-2.4,ik)
function cbar(x,y,k)
 'query shades'
  shdinfo = result
  if (subwrd(shdinfo,1)='None')
    say 'Cannot plot color bar: No shading information'
    return
  endif
*
  cnum = subwrd(shdinfo,5)
*
*  Plot colorbar
*
   dx=0.07
   dy=0.12
    x=x
   y0=y

  num = 0
  while (num<cnum)

    rec = sublin(shdinfo,num+2)
    col = subwrd(rec,1)
     hi = subwrd(rec,3)
    'set line 'col
    'draw recf 'x-dx' 'y' 'x' 'y+dy
    'set line 1  1  2'
    'draw rec  'x-dx' 'y' 'x' 'y+dy
     if(num<cnum-1)
      'set string 1 l 2 0' ; 'set strsiz 0.07 0.08'
      'draw string 'x+0.02' 'y+dy' `0'hi
     endif
    num = num + 1
    y=y+dy
  endwhile

if(ik=1);'draw string 'x+0.02' 'y+dy'(km)';endif

NovJanPosQBO.txt text file code:

years others
1964 1960
1967 1961
1976 1962
1983 1963
1986 1965
1988 1966
1991 1968
1993 1969
2009 1970
2011 1971
     1972
     1973
     1974
     1975
     1977
     1978
     1979
     1980
     1981
     1982
     1984
     1985
     1987
     1989
     1990
     1992
     1994
     1995
     1996
     1997
     1998
     1999
     2000
     2001
     2002
     2003
     2004
     2005
     2006
     2007
     2008
     2010
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20140128/ee0ef9de/attachment-0001.html 