;@/usr/local/rsi/idl/lib/arrow.pro ; remove this for mac & window versions pro mini_plan,date,latpole=lat00,lonpole=lon00,kiruna=kiruna,takeoff_time=takeoff_time,$ dryden=dryden,scale=scale0,proj=proj,no_haloe=no_haloe,no_sage=no_sage,no_small_plots=no_small_plots,nmc=nmc,$ old_plan=old_plan,color=color,airspeed=airspeed,solar_zenith_cont=solar_zenith_cont,$ standard=standard,ndays=ndays ;+ ; this routine does a sza plot for the location ; clidk on map to find the location ; plot appears in second window ; inputs: ; date= date in 'yymmdd' format ; keyword inputs: ; latpole=latitude map pole ; lonpole = longitude map pole ; kiruna = set the map pole to kiurna ; dryden = set the map pole to dryden ; scale = zoom in parameter (default 4.e7 ) - larger numbers to zoom out, smaller to zoom in ; takeoff_time = local takeoff time in hours and a fraction ; nmc = use nmc forecast instead of default ; solar_zenith_cont = plot contours of solar zenith angle ; ; notes ; ; takeoff_gmt is the takeoff time in gmt ; local_to_gmt is the local time to gmt offset local_time=gmt-local_to_gmt ; proj values: ; 1 = stereographic ; 2 = orthographic ; 3 = conic ; 4 = lambert equal ; 5 = gnomic ; 9 = mercator ; 10 = mollweide ; 11 = sine ; 14 = albers ; ;- ; set defaults theta=460 if n_elements(proj) eq 0 then proj=4 ; use the statements below for X window ;set_plot,'X' ;dir='/app/ftp/pub/solve/met_data/mini_plan/' ; this is the location of the data files ;dir='/home/schom/solve/' ; use the statements below for MAC set_plot,'MAC' dir='Macintosh HD:Applications:RSI:IDL 5.2:my idl stuff:' ; this is the location of the data files if n_elements(airspeed) eq 0 then airspeed=828. ; DC-8 speed in km/hr !quiet=1 if !d.name eq 'X' then !p.charsize=0.5 if n_elements(date) eq 0 then date=today(3) if n_elements(scale0) eq 0 then scale0=4.0e7 if n_elements(ndays) eq 0 then ndays=5 ; delete all open windows while !d.window ne -1 do tvdelete,!d.window if keyword_set(kiruna) then begin lat_home=67.84 lon_home=20.41 city='Kiruna' local_to_gmt=1. if n_elements(takeoff_time) eq 0 then takeoff_time=9. endif if keyword_set(dryden) then begin lat_home=34.92 lon_home=-117.89+360. city='Dryden' local_to_gmt=-8 if n_elements(takeoff_time) eq 0 then takeoff_time=16. endif if n_elements(lat_home) eq 0 then lat_home=67.84 if n_elements(lon_home) eq 0 then lon_home=20.41 if n_elements(city) eq 0 then city='Kiruna' if n_elements(local_to_gmt) eq 0 then local_to_gmt=1. if n_elements(takeoff_time) eq 0 then takeoff_time=9. if n_elements(old_plan) gt 0 then begin lat_pathx=0. lon_pathx=0. restore,dir+old_plan ;,takeoff_time,lat00,lon00,lat_home,lon_home,city,lat_path,lon_path,scale0 endif if n_elements(lat00) eq 0 then lat00=lat_home if n_elements(lon00) eq 0 then lon00=lon_home date_arr=strarr(4) for j=0,3 do date_arr(j)=newdatex(date,j-3) n=where(date_arr eq date) id_date=n(0) lat_path=[lat_home] ; way point locations lon_path=[lon_home] local_to_gmt=float(local_to_gmt) takeoff_gmt=takeoff_time-local_to_gmt lat_pathd=lat_path ; the path with a denser series of points lon_pathd=lon_path segment=[1.] iseg=0 flttime=[0] ; way point flight times flttimed=flttime ; a denser series of flight times yr=fix(strmid(date,0,2)) mon=fix(strmid(date,2,2)) day=fix(strmid(date,4,2)) szz=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime) szzm=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime,/moon) dy19=day1900(yr,mon,day) read_sage2,dd,sage_lat,sage_lon,dir=dir read_haloe,ddh,haloe_lat,haloe_lon,dir=dir read_poam,ddp,poam_lat,poam_lon,aztop=aztop,dir=dir read_forecast,date_arr(0),pv_arr0,t_arr0,m_arr0,lats0,lons0,theta0,err=err0,source=source0,$ trop_hgt=trop_hgt0,nmc=nmc,standard=standard,dir=dir read_forecast,date_arr(1),pv_arr1,t_arr1,m_arr1,lats1,lons1,theta1,err=err1,source=source1,$ trop_hgt=trop_hgt1,nmc=nmc,standard=standard,dir=dir read_forecast,date_arr(2),pv_arr2,t_arr2,m_arr2,lats2,lons2,theta2,err=err2,source=source2,$ trop_hgt=trop_hgt2,nmc=nmc,standard=standard,dir=dir read_forecast,date_arr(3),pv_arr3,t_arr3,m_arr3,lats3,lons3,theta3,err=err3,source=source3,$ trop_hgt=trop_hgt3,nmc=nmc,standard=standard,dir=dir if not keyword_set(nmc) then begin ; load up nmc plots for comparison - unless nmc is selected read_forecast,date_arr(0),pv_arrn0,t_arrn0,m_arrn0,latsn0,lonsn0,thetan0,err=errn0,source=sourcen0,$ trop_hgt=trop_hgtn0,/nmc,dir=dir read_forecast,date_arr(1),pv_arrn1,t_arrn1,m_arrn1,latsn1,lonsn1,thetan1,err=errn1,source=sourcen1,$ trop_hgt=trop_hgtn1,/nmc,dir=dir read_forecast,date_arr(2),pv_arrn2,t_arrn2,m_arrn2,latsn2,lonsn2,thetan2,err=errn2,source=sourcen2,$ trop_hgt=trop_hgtn2,/nmc,dir=dir read_forecast,date_arr(3),pv_arrn3,t_arrn3,m_arrn3,latsn3,lonsn3,thetan3,err=errn3,source=sourcen3,$ trop_hgt=trop_hgtn3,/nmc,dir=dir endif irt=1 if id_date ge 0 then begin if id_date eq 0 then begin source=source0 & data_err=err0 if err0 eq 0 then begin & pv_arr=pv_arr0 & t_arr=t_arr0 & m_arr=m_arr0 & lats=lats0 & lons=lons0 & theta=theta0 & trop_hgt=trop_hgt0 & endif endif if id_date eq 1 then begin data_err=err1 & source=source1 if err1 eq 0 then begin & pv_arr=pv_arr1 & t_arr=t_arr1 & m_arr=m_arr1 & lats=lats1 & lons=lons1 & theta=theta1 & trop_hgt=trop_hgt1 & endif endif if id_date eq 2 then begin if err2 eq 0 then begin & pv_arr=pv_arr2 & t_arr=t_arr2 & m_arr=m_arr2 & lats=lats2 & lons=lons2 & theta=theta2 & trop_hgt=trop_hgt2 & endif data_err=err2 & source=source2 endif if id_date eq 3 then begin data_err=err3 & source=source3 if err3 eq 0 then begin & pv_arr=pv_arr3 & t_arr=t_arr3 & m_arr=m_arr3 & lats=lats3 & lons=lons3 & theta=theta3 & trop_hgt=trop_hgt3 & endif endif endif else begin source='No data' data_err=-1 endelse device,retain=2,/cursor_orig,pseudo=8 win1=!d.window+1 ; main plot window win2=!d.window+2 ; zenith angle window win3=!d.window+3 win4=!d.window+4 win5=!d.window+5 win6=!d.window+6 win7=!d.window+7 ; PV and temperature plot window if not keyword_set(nmc) then begin ; add more little windows for nmc plots for comparison - unless nmc is selected winn3=!d.window+7 winn4=!d.window+8 winn5=!d.window+9 winn6=!d.window+10 endif loadct,39 ; little plot windows if keyword_set(no_small_plots) then goto,sk5 if err0 eq 0 then begin window,win3,xpos=10,ypos=500,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(theta0)+'K PV '+source0+' '+nicedate(date_arr(0)) field_map,pv_arr0,t_arr0,m_arr0,lons0,lats0,theta0,trop_hgt0,lat_home,lon_home map_continents map_grid endif if err1 eq 0 then begin window,win4,xpos=10+220,ypos=500,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(theta1)+'K PV '+source1+' '+nicedate(date_arr(1)) field_map,pv_arr1,t_arr1,m_arr1,lons1,lats1,theta1,trop_hgt1,lat_home,lon_home map_continents map_grid endif if err2 eq 0 then begin window,win5,xpos=10+440,ypos=500,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(theta2)+'K PV '+source2+' '+nicedate(date_arr(2)) field_map,pv_arr2,t_arr2,m_arr2,lons2,lats2,theta2,trop_hgt2,lat_home,lon_home map_continents map_grid endif if err3 eq 0 then begin window,win6,xpos=10+660,ypos=500,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(theta3)+'K PV '+source3+' '+nicedate(date_arr(3)) field_map,pv_arr3,t_arr3,m_arr3,lons3,lats3,theta3,trop_hgt3,lat_home,lon_home map_continents map_grid endif if not keyword_set(nmc) then begin ; add more windows for nmc plots for comparison - unless nmc is selected yoff=500-250 if errn0 eq 0 then begin window,winn3,xpos=10,ypos=yoff,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(thetan0)+'K PV '+sourcen0+' '+nicedate(date_arr(0)) field_map,pv_arrn0,t_arrn0,m_arrn0,lonsn0,latsn0,thetan0,trop_hgtn0,lat_home,lon_home map_continents map_grid endif if errn1 eq 0 then begin window,winn4,xpos=10+220,ypos=yoff,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(thetan1)+'K PV '+sourcen1+' '+nicedate(date_arr(1)) field_map,pv_arrn1,t_arrn1,m_arrn1,lonsn1,latsn1,thetan1,trop_hgtn1,lat_home,lon_home map_continents map_grid endif if errn2 eq 0 then begin window,winn5,xpos=10+440,ypos=yoff,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(thetan2)+'K PV '+sourcen2+' '+nicedate(date_arr(2)) field_map,pv_arrn2,t_arrn2,m_arrn2,lonsn2,latsn2,thetan2,trop_hgtn2,lat_home,lon_home map_continents map_grid endif if errn3 eq 0 then begin window,winn6,xpos=10+660,ypos=yoff,xsize=200,ysize=200 erase map_set,lat00,lon00,proj=proj,scale=1.e8,title=str(thetan3)+'K PV '+sourcen3+' '+nicedate(date_arr(3)) field_map,pv_arrn3,t_arrn3,m_arrn3,lonsn3,latsn3,thetan3,trop_hgtn3,lat_home,lon_home map_continents map_grid endif endif ; NMC sk5: ; open main window window,win1,xsize=450,ysize=450, title='Map Window' erase ; draw the map of the region and overplot fields map_set,lat00,lon00,proj=proj,scale=scale0,title=str(theta)+'K PV '+source+' '+nicedate(date)+$ ' Takeoff GMT '+hrsmins(takeoff_gmt) + ' Local '+hrsmins(takeoff_time) sat_col=255 path_col=200 if data_err eq 0 then begin sat_col=0 path_col=20 field_map,pv_arr,t_arr,m_arr,lons,lats,theta,trop_hgt,lat_home,lon_home endif map_continents,col=230,/countries,/coasts map_grid ; plot 4 hour range ring range_ring,lat_home,lon_home,airspeed*4,360,bearing,latp,lonp oplot,lonp,latp,col=50 ; plot great circle route from dryden to kiruna gr_circ_rte,34.9,360.-117.9,67.84,20.48,30,bearing,dist,del,latp,lonp,rd oplot,lonp,latp,col=50 ; plot 500 km range rings out to 3000 km for j=500,3000,500 do begin ; 500 km range rings range_ring,lat_home,lon_home,j,360,bearing,latp,lonp oplot,lonp,latp,col=75 endfor ; plot location of city oplot,[lon_home],[lat_home],psym=1,col=80 xyouts,lon_home,lat_home,city,col=sat_col symb,10,/fill ; plot poam points plot_sat_points,ddp,dy19,takeoff_gmt,poam_lat,poam_lon,sat_col,local_to_gmt,times=poam_times,sat_char=' ',/print,aztop=aztop ; plot sage 2 points if not keyword_set(no_sage) then $ plot_sat_points,dd,dy19,takeoff_gmt,sage_lat,sage_lon,sat_col,local_to_gmt,times=sage_times,sat_char=' S' ; plot haloe points if not keyword_set(no_haloe) then $ plot_sat_points,ddh,dy19,takeoff_gmt,haloe_lat,haloe_lon,sat_col,local_to_gmt,times=haloe_times,sat_char=' H' ; plot solar zenith contours if keyword_set(solar_zenith_cont) then begin latbase=findgen(60)*2.-30. lonbase=findgen(91)*4. sz2d=fltarr(n_elements(lonbase),n_elements(latbase)) for j=0,n_elements(latbase)-1 do $ for i=0,n_elements(lonbase)-1 do $ sz2d(i,j)=solar_zenith(latbase(j),lonbase(i),mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime) contour,sz2d,lonbase,latbase,c_col=220,/over,/follow,levels=[70.,75.,80.,85.,90.,95.,100.,105.] endif ; zenith angle window window,win2,xsize=400,ysize=300,title='Flight window',xpos=10,ypos=10 zen_plot,date,takeoff_time,takeoff_gmt,poam_times,sage_times,haloe_times print,'Wy Pt Flt Time Local GMT Lat Lon Nom. Flight Len (hrs)' print,str(1),' '+hrsmins(0.),' '+hrsmins(takeoff_time),' '+hrsmins(takeoff_gmt),lat_home,lon_home,hrsmins(0.); ok all done with set up istar=1 ; back to map window xx: ; find a location wset,win1 symb,12 map_set,lat00,lon00,proj=proj,scale=scale0,/noerase if n_elements(old_plan) gt 0 then begin ; skip this section if loading an old plan lat0=lat_pathx(istar) lon0=lon_pathx(istar) istar=istar+1 endif else cursor,lon0,lat0,/data,/down ; plot our tenative point oplot,[lon0],[lat0],psym=8,col=path_col,symsize=1.2 wset,win2 ; go to zenith angle window ns=n_elements(lat_path) oldlat=lat_path(ns-1) oldlon=lon_path(ns-1) gr_circ_rte,oldlat,oldlon,lat0,lon0,40,bearing,dist,del,latp,lonp,rd ; find the distance gr_circ_rte,oldlat,oldlon,lat0,lon0,1,bearing,dist2,del,latp1,lonp1,rd2 ; find the distance lat_path3=[lat_path,latp1] lon_path3=[lon_path,lonp1] lat_path2=[lat_pathd,latp] lon_path2=[lon_pathd,lonp] flttime2=[flttimed,flttimed(n_elements(flttimed)-1)+rd/airspeed] ; this is how long it takes to fly there flttimex=[flttime,flttime(n_elements(flttime)-1)+rd2/airspeed] ; this is how long it takes to fly there szzx=solar_zenith(lat_path3,lon_path3,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttimex) szz2=solar_zenith(lat_path2,lon_path2,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime2) szzm2=solar_zenith(lat_path2,lon_path2,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime2,/moon) zen_plot,date,takeoff_time,takeoff_gmt,poam_times,sage_times,haloe_times oplot,flttime2+takeoff_time,szz2 oplot,flttimex+takeoff_time,szzx,psym=1 ; way point spots oplot,flttime2+takeoff_time,szzm2,linestyle=1 mp=max(flttime2) ; now compute the time it takes to get home gr_circ_rte,lat0,lon0,lat_home,lon_home,20,bearing,dist,del,latp,lonp,rd ; find the distance home_time=dist/airspeed rr=n_elements(flttimex) print,str(rr)-1,' '+hrsmins(mp),' '+hrsmins(mp+takeoff_time),' '+hrsmins(mp+takeoff_gmt),lat0,lon0,hrsmins(mp+home_time) ; now decide if this is a permanent point if n_elements(old_plan) gt 0 then begin xx1=1 if istar gt n_elements(lat_pathx)-1 then xx1=-1 goto,xy1 endif else widget2,xx1 ; draw widget xx1=1 means accept the point, xx1= 0 means no, xx1=-1 means finished if xx1 eq 0 then goto,xx ; we have selected a point xy1: ns=n_elements(lat_path) oldlat=lat_path(ns-1) oldlon=lon_path(ns-1) iseg=iseg+1 ; go to map window and plot the path wset,win1 map_set,lat00,lon00,proj=proj,scale=scale0,/noerase gr_circ_rte,oldlat,oldlon,lat0,lon0,40,bearing,dist2,del,latp,lonp,rdx oplot,lonp,latp,col=path_col,thick=2 lat_path=[lat_path,lat0] lon_path=[lon_path,lon0] segment=[segment,latp*0. + iseg] lat_pathd=[lat_pathd,latp] lon_pathd=[lon_pathd,lonp] gr_circ_rte,oldlat,oldlon,lat0,lon0,2,bearing,dist,del,latp,lonp,rd flttime=[flttime,flttime(ns-1)+dist/airspeed] nss=n_elements(flttimed) flttimed=[flttimed,flttimed(nss-1)+rdx/airspeed] ; switch to flight window and plot zenith angles wset,win2 szz=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime) szzm=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime,/moon) szzd=solar_zenith(lat_pathd,lon_pathd,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttimed) szzmd=solar_zenith(lat_pathd,lon_pathd,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttimed,/moon) zen_plot,date,takeoff_time,takeoff_gmt,poam_times,sage_times,haloe_times oplot,flttime+takeoff_time,szz,psym=1 oplot,flttime+takeoff_time,szzm,psym=1 oplot,flttimed+takeoff_time,szzd oplot,flttimed+takeoff_time,szzmd,linestyle=1 if flttime(n_elements(flttime)-1) gt 10. then goto,quitit if xx1 eq -1 then goto,quitit ; quit plotting points and clean up wset,win1 map_set,lat00,lon00,proj=proj,/noerase,scale=scale0 goto,xx print,'out of range' quitit: ; okay we are all done now - print out results print,'Way Pt. Flt TimeLocal GMT Lat Lon Solar Z Lunar Z' for j=0,n_elements(flttime)-1 do print,j+1,' '+hrsmins(flttime(j)), ' '+hrsmins(flttime(j)+takeoff_time),' '+hrsmins(flttime(j)+takeoff_gmt),lat_path(j),$ lon_path(j),szz(j),szzm(j) lat_pathx=lat_path lon_pathx=lon_path if n_elements(old_plan) eq 0 then begin ; save plan only if this is a new plan name='fltplan_'+date+'-'+strmid(systime(0),11,2)+'-'+strmid(systime(0),14,2) save,file=dir+name,date,takeoff_time,lat00,lon00,lat_home,lon_home,city,lat_pathx,lon_pathx,local_to_gmt,scale0 print,'saved: '+name old_plan=name endif ; now create temperature and pv plot window if data_err eq 0 then begin window,win7,xsize=400,ysize=300,title='MPV and Temperature',xpos=10,ypos=350 s=size(pv_arr) if lons(0) eq 180. then begin pv_arr2=rotate_array(pv_arr,lons,outlon=lons2) t_arr2=rotate_array(t_arr,lons,outlon=lons2) hgt_2=rotate_array(trop_hgt,lons,outlon=lons2) endif else begin pv_arr2=pv_arr(0:s(1)-2,*) t_arr2=t_arr(0:s(1)-2,*) hgt_2=trop_hgt(0:s(1)-2,*) lons2=lons(0:s(1)-2) endelse ; interpolate over the flight path pv=latlonintrp(pv_arr2,lons2,lats,lon_pathd,lat_pathd,/wrap,/reglon,/regular) temper=latlonintrp(t_arr2,lons2,lats,lon_pathd,lat_pathd,/wrap,/regular) hgtt=latlonintrp(hgt_2,lons2,lats,lon_pathd,lat_pathd,/wrap,/regular) yrange=[1.0,3.5] xrange=takeoff_time+[0,10] position=[.15,.1,.78,.85] plot,flttimed+takeoff_time,pv,title='MPV and Temperature '+str(theta)+'K '+nicedate(date),xrange=xrange,yrange=yrange,$ xtitle='Flight Time (local)',ytitle='MPV',xst=1,yst=9,xticks=10.,position=position gmts=takeoff_gmt+findgen(10.) n=where(gmts gt 24.) if n(0) gt 0 then gmts(n)=gmts(n)-24. xinc=takeoff_time+findgen(10) xyouts,xinc,yrange(0)+.05+fltarr(10),string(gmts,format='(f4.1)') xyouts,takeoff_time-1.5,yrange(0)+.05,'GMT',size=0.9 xyouts,xinc,yrange(0)+.1+fltarr(10),string(findgen(10),format='(i2)') xyouts,takeoff_time-1.5,yrange(0)+.1,'FLT',size=0.9 oplot,(takeoff_time+8.)*[1.,1.],yrange legend,['MPV','T','Trop H'],lines=[0,2,1],box=0 axis,yaxis=1,yrange=[180,220],/save,ytitle='Temperature' oplot,flttimed+takeoff_time,temper,linesty=2 tnat=nat_theta(10.e-9,4.e-6,theta,type=1) tice=frost_theta(4.e-6,theta) oplot,takeoff_time+[0,10.],[tnat,tnat],linest=2 xyouts,takeoff_time+.5,tnat,'NAT',/data oplot,takeoff_time+[0,10.],[tice,tice],linest=2 xyouts,takeoff_time+.5,tice,'Ice',/data xrange=takeoff_time+[0,11.5] plot,/noerase,flttimed+takeoff_time,hgtt*1.0e-3,position=[position(0),position(1),position(0)+(position(2)-position(0))*11.5/10.,position(3)],$ xst=5,yst=5,yrange=[5,15],xrange=xrange,linesty=1 axis,yaxis=1,yrange=[5,15],/save,ytitle='Trop Height (km)' ; oplot,flttimed+takeoff_time,hgtt*1.0e-3,linesty=1 endif ; trajectory calculation widget3,rr1,'Run back trajectory?' if rr1 eq 0 then goto,bbr mtraj,startday=date,starttime=12.,ndays=ndays,forward=-1,theta=theta,$ latin=lat_pathd,lonin=lon_pathd,dir=dir,scale=scale0,latpole=lat00,$ lonpole=lon00,segment=segment bbr: widget3,rr1,'Run RDF?' if rr1 eq 0 then goto,bb1 rdf_gen,date,ndays,latpole=latpole,lonpole=lonpole,dir=dir,latdraw=lat_pathd,londraw=lon_pathd bb1: widget3,rr1,'Print plan?' if rr1 eq 0 then goto,endall print_plan,old_plan,nmc=nmc,color=color,airspeed=airspeed,scale0=scale0,solar_zenith_cont=solar_zenith_cont,dir=dir endall: return end