pro plot_gpm_l2 fdir='/uufs/chpc.utah.edu/common/home/mace-group4/gpm/gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L2/GPM_2BCMB.05/2018/001/' fname='2B.GPM.DPRGMI.CORRA2016.20180101-S071950-E085225.021837.V05A.HDF5' presult=h5_parse(fdir+fname) ; Open an existing hdf5 file file_id=h5f_open(fdir+fname) ; Open the top level (root) root_id=h5g_open(file_id,'/') ; Root / what (object, version, date, time, source) group_id=h5g_open(root_id,'MS') dset_id=h5d_open(group_id,'Latitude') lat=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(group_id,'Longitude') lon=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(group_id,'cloudLiqWaterCont') lwc=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(group_id,'airPressure') pres=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(group_id,'correctedReflectFactor') reflect=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(group_id,'simulatedBrightTemp') tb=h5d_read(dset_id) h5d_close,dset_id time_id=h5g_open(group_id,'ScanTime') dset_id=h5d_open(time_id,'DayOfYear') doy=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Year') year=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Month') month=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'DayOfMonth') day=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Hour') hour=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Minute') minute=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'Second') second=h5d_read(dset_id) h5d_close,dset_id dset_id=h5d_open(time_id,'MilliSecond') msecond=h5d_read(dset_id) h5d_close,dset_id h5g_close,time_id h5g_close,group_id h5g_close,root_id lat=transpose(lat) lon=transpose(lon) lwc=transpose(lwc) pres=transpose(pres) reflect=transpose(reflect) tb=transpose(tb) msecond=msecond*1e-3 second=second+msecond julian_day=julday(month,day,year,hour,minute,second) ;************************************************* ; Plot lat-lon track and timeseries of variables ;************************************************* do_plot='yes' if do_plot eq 'yes' then begin pxdim=900 pydim=700 dummy=label_date(date_format=['%Z/%M/%D!C%H:%I']) ; Set up the map coordinates p0=map('cylindrical',dimensions=[pxdim,pydim],/buffer,$ position=[0.1,0.4,0.9,0.9]) m1=mapcontinents(fill_color="pale goldenrod") p1=plot(lon[*,0],lat[*,0],color='blue',thick=3,/overplot) ;geoprof ;p3=plot(mlons_cs,mlats_cs,color='red',thick=1,/overplot) var=tb[*,0,0] result=where(var gt -9999.90,count) dmin=min(var[result]) dmax=max(var[result]) p2=plot(julian_day,var,/current,$ position=[0.1,0.1,0.9,0.4],xtickformat='label_date',$ /xstyle,/ystyle,color='red',yrange=[dmin,dmax]) ;imagename='/uufs/chpc.utah.edu/common/home/u0079358/cs_tracks_joint/images/'+$ ; 'modis_opt_thickness.'+orbit+'.png' imagename='map.png' p0.save,imagename,height=pydim endif ;end of plot to check stop ;*************************** ; Make a plot of this subset ;*************************** ; Set up the plot for overpass pxdim=900 & pydim=900 ; Position the plots xl=0.07 & xr=0.80 yb=0.05 & yt=0.60 sx=0.07 sy=0.05 numplots_x=2 numplots_y=2 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; Colorbar position cbpos=pos cbpos[*,0]=pos[*,2]+0.1 cbpos[*,2]=cbpos[*,0]+0.01 ; Color table ; Top is the last color to scale 256 colors, 0-255 top_color=252 ; Colortable 0-252 253=white mytable=colortable(39,ncolors=254) revmytable=reverse(mytable) ;254=hot pink ;gray=255 mytable=[mytable,transpose([238,18,137]),transpose([150,150,150])] revmytable=[revmytable,transpose([238,18,137]),transpose([150,150,150])] cbtable=mytable[0:top_color,*] revcbtable=revmytable[0:top_color,*] dx=0.01 & dy=0.01 font_size=12 ;Clear plot space pnum=0 p0=plot([0,1],[0,1],position=pos[pnum,*],/buffer,dimensions=[pxdim,pydim],$ axis_style=4,/nodata) d=p0.convertcoord([pos[pnum,0],pos[pnum,2]],[pos[pnum,1],pos[pnum,3]],/normal,/to_device) isx=(d[0,1,0]-d[0,0,0]) isy=(d[1,1,0]-d[1,0,0]) ; Set up the map ulat=-10.0 & llat=-70.0 llon=70 & rlon=180+45 m0=map('mercator',box_axes=1,label_position=0,linestyle='dotted',$ position=pos[pnum,*],/current,font_size=font_size);,limit=[llat,llon,ulat,rlon]) m1=mapcontinents(fill_color="pale goldenrod") p1=plot(lon[*,0],lat[*,0],color='blue',thick=3,/overplot) p0.save,'map.png',height=pydim stop ; Plot latitude pnum=0 dmin=min(lat) dmax=max(lat) data_image=bytscl(lat,top=top_color,min=dmin,max=dmax) p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],$ rgb_table=mytable) ;c1=contour(data_image,lat,height_regrid[hidx],ytitle='height (km)',$;xtitle='latitude',$ ; xstyle=1,ystyle=1,font_size=font_size,/nodata,/current,position=pos[pnum,*],$ ; irregular=0,xtickdir=1,ytickdir=1,xtickformat='(A1)',axis_style=2) c2=colorbar(orientation=1,rgb_table=cbtable,range=[dmin,dmax],title='degree',$ position=reform(cbpos[pnum,*]),font_size=font_size,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+dy,'Latitude') ; Plot longitude pnum=1 dmin=min(lon) dmax=max(lon) data_image=bytscl(lon,top=top_color,min=dmin,max=dmax) p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],$ rgb_table=mytable) ;c1=contour(data_image,lat,height_regrid[hidx],ytitle='height (km)',$;xtitle='latitude',$ ; xstyle=1,ystyle=1,font_size=font_size,/nodata,/current,position=pos[pnum,*],$ ; irregular=0,xtickdir=1,ytickdir=1,xtickformat='(A1)',axis_style=2) c2=colorbar(orientation=1,rgb_table=cbtable,range=[dmin,dmax],title='degree',$ position=reform(cbpos[pnum,*]),font_size=font_size,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+dy,'Longitude') ; Histogram latitude pnum=1 data_var=lat start_bin=dmin end_bin=dmax dbin=(end_bin-start_bin)/100. hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1 p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos_hist[pnum,*],$ ytitle='frequency',xtitle='longitude') ; Plot lwc pnum=2 data_var=pres dmin=min(data_var) dmax=max(data_var) data_image=bytscl(data_var,top=top_color,min=dmin,max=dmax) p1=image(data_image,/current,image_dimensions=[isx,isy],position=pos[pnum,*],$ rgb_table=mytable) ;c1=contour(data_image,lat,height_regrid[hidx],ytitle='height (km)',$;xtitle='latitude',$ ; xstyle=1,ystyle=1,font_size=font_size,/nodata,/current,position=pos[pnum,*],$ ; irregular=0,xtickdir=1,ytickdir=1,xtickformat='(A1)',axis_style=2) c2=colorbar(orientation=1,rgb_table=cbtable,range=[dmin,dmax],title='degree',$ position=reform(cbpos[pnum,*]),font_size=font_size,tickdir=1,border_on=1) t1=text(pos[pnum,0]+dx,pos[pnum,3]+dy,'LWC') ; Histogram latitude pnum=2 ;data_var=lwc start_bin=dmin end_bin=dmax dbin=(end_bin-start_bin)/100. hist_generic,data_var,start_bin,end_bin,dbin,bins1,data_freq1,data_counts1 p1=plot(bins1,data_freq1,/current,color='black',/hist,position=pos_hist[pnum,*],$ ytitle='frequency',xtitle='lwc') p0.save,'gpm.png',height=pydim stop ;date='20180104' ; ;syear=strmid(date,0,4) ;smonth=strmid(date,4,2) ;sday=strmid(date,6,2) ;julian_date,syear,smonth,sday,doy ;sdoy=string(doy,format='(I03)') ; ;; Build the first part of the url 1hr time-averaged ;;url = 'http://sbenson:R1a3i5n*@goldsmr4.sci.gsfc.nasa.gov:80/opendap/MERRA2/M2T1NXSLV.5.12.4/'+$ ;; syear+'/'+smonth+'/MERRA2_'+asm_stream+'.tavg1_2d_slv_Nx.'+syear+smonth+sday+'.nc4?' ;url='https://gpm1.gesdisc.eosdis.nasa.gov:443/opendap/GPM_L2/GPM_2BCMB.05/'+syear+'/'+$ ; sdoy+'/2B.GPM.DPRGMI.CORRA2016.20180104-S000846-E014121.021879.V05A.HDF5?' ; ;; Second part of the url that selects out the time, lat and lon array ;url2='MS_ScanTime_SecondOfDay[0:1:7935]';,lat[0:1:360],lon[0:1:575]' ;print,'getting time,lat,lon' ;file_id=ncdf_open(url+url2) ; ;; Get the time variable from the merra file ;; minutes since file date 00:00:00 ;var_id=ncdf_varid(file_id,'MS_ScanTime_SecondOfDay') ;ncdf_varget,file_id,var_id,minutes ; ;ncdf_close,file_id ; stop ;; convert merra time to julian day ;julian_day_merra=julday(smonth,sday,syear,0,minutes,0) ;caldat,julian_day_merra,mm,dd,yy,hh,mi,ss ;;print,mm,dd,yy,hh,mi,ss ; ; ; Get latitude array ; var_id=ncdf_varid(file_id,'lat') ; ncdf_varget,file_id,var_id,all_lat_merra ; ; ; Get longitude array ; var_id=ncdf_varid(file_id,'lon') ; ncdf_varget,file_id,var_id,all_lon_merra ; ; ncdf_close,file_id ; Find the closest merra time to synoptic closest=min(abs(julian_day_synoptic[i]-julian_day_merra),time_idx) caldat,julian_day_merra[time_idx],mm_merra,dd_merra,yy_merra,hh_merra,mi_merra,ss_merra print,mm_merra,dd_merra,yy_merra,hh_merra,mi_merra,ss_merra ; Find the latitude range lat_idx=where(all_lat_merra le ulat+1.0 and all_lat_merra ge llat-1.0,count) lat_merra=all_lat_merra[lat_idx] ; Strings of the indices to use for lat and lon slat0=string(lat_idx[0],format='(I03)') slat1=string(lat_idx[-1],format='(I03)') ; Select lat indexes url2='H250['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'U250['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'V250['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'H500['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'U500['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'V500['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'Q850['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'T850['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'U850['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'V850['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'H850['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'H1000['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'U2M['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575],'+$ 'V2M['+string(time_idx,format='(I03)')+']['+slat0+':'+slat1+'][0:1:575]' print,'getting profiles' file_id=ncdf_open(url+url2) var_id=ncdf_varid(file_id,'H250') ncdf_varget,file_id,var_id,h250_merra var_id=ncdf_varid(file_id,'U250') ncdf_varget,file_id,var_id,u250_merra var_id=ncdf_varid(file_id,'V250') ncdf_varget,file_id,var_id,v250_merra var_id=ncdf_varid(file_id,'H500') ncdf_varget,file_id,var_id,h500_merra var_id=ncdf_varid(file_id,'U500') ncdf_varget,file_id,var_id,u500_merra var_id=ncdf_varid(file_id,'V500') ncdf_varget,file_id,var_id,v500_merra var_id=ncdf_varid(file_id,'Q850') ncdf_varget,file_id,var_id,q850_merra var_id=ncdf_varid(file_id,'T850') ncdf_varget,file_id,var_id,t850_merra var_id=ncdf_varid(file_id,'U850') ncdf_varget,file_id,var_id,u850_merra var_id=ncdf_varid(file_id,'V850') ncdf_varget,file_id,var_id,v850_merra var_id=ncdf_varid(file_id,'H850') ncdf_varget,file_id,var_id,h850_merra var_id=ncdf_varid(file_id,'H1000') ncdf_varget,file_id,var_id,h1000_merra var_id=ncdf_varid(file_id,'U2M') ncdf_varget,file_id,var_id,u2m_merra var_id=ncdf_varid(file_id,'V2M') ncdf_varget,file_id,var_id,v2m_merra ncdf_close,file_id print,'done getting profiles' ;***************************** ; Subset out the longitudes ;***************************** if rlon le 180 then begin lon_idx=where(all_lon_merra le rlon+1.0 and all_lon_merra ge llon-1.0,count) lon_merra=all_lon_merra[lon_idx] h250_merra=h250_merra[lon_idx,*] u250_merra=u250_merra[lon_idx,*] v250_merra=v250_merra[lon_idx,*] h500_merra=h500_merra[lon_idx,*] u500_merra=u500_merra[lon_idx,*] v500_merra=v500_merra[lon_idx,*] q850_merra=q850_merra[lon_idx,*] t850_merra=t850_merra[lon_idx,*] u850_merra=u850_merra[lon_idx,*] v850_merra=v850_merra[lon_idx,*] h850_merra=h850_merra[lon_idx,*] h1000_merra=h1000_merra[lon_idx,*] u2m_merra=u2m_merra[lon_idx,*] v2m_merra=v2m_merra[lon_idx,*] endif else if rlon gt 180 then begin ; Calculate longitude as 0 to 360 colon_merra=all_lon_merra result=where(colon_merra lt 0,count) if count gt 0 then colon_merra[result]=360.0+colon_merra[result] ; Calculate map boundaries in colongitude if rlon lt 0 then rcolon=360.0+rlon else rcolon=rlon if llon lt 0 then lcolon=360.0+llon else lcolon=llon ; Reorder the colon array idx=sort(colon_merra) ; Reorder the data before subsetting colon_merra=colon_merra[idx] h250_merra=h250_merra[idx,*] u250_merra=u250_merra[idx,*] v250_merra=v250_merra[idx,*] h500_merra=h500_merra[idx,*] u500_merra=u500_merra[idx,*] v500_merra=v500_merra[idx,*] q850_merra=q850_merra[idx,*] t850_merra=t850_merra[idx,*] u850_merra=u850_merra[idx,*] v850_merra=v850_merra[idx,*] h850_merra=h850_merra[idx,*] h1000_merra=h1000_merra[idx,*] u2m_merra=u2m_merra[idx,*] v2m_merra=v2m_merra[idx,*] ; Subset data within colon bounds lon_idx=where(colon_merra le rcolon+1.0 and colon_merra ge lcolon-1.0,count) colon_merra=colon_merra[lon_idx] h250_merra=h250_merra[lon_idx,*] u250_merra=u250_merra[lon_idx,*] v250_merra=v250_merra[lon_idx,*] h500_merra=h500_merra[lon_idx,*] u500_merra=u500_merra[lon_idx,*] v500_merra=v500_merra[lon_idx,*] q850_merra=q850_merra[lon_idx,*] t850_merra=t850_merra[lon_idx,*] u850_merra=u850_merra[lon_idx,*] v850_merra=v850_merra[lon_idx,*] h850_merra=h850_merra[lon_idx,*] h1000_merra=h1000_merra[lon_idx,*] u2m_merra=u2m_merra[lon_idx,*] v2m_merra=v2m_merra[lon_idx,*] rlon=rcolon llon=lcolon lon_merra=colon_merra ; convert colon back to longitude lon_merra=colon_merra result=where(colon_merra lt 0,count) if count gt 0 then lon_merra[result]=colon_merra[result]-360.0 endif ;*************** ; Replace missing values with nan. ; Nan is recognized by IDL as the missing value ;*************** result=where(u850_merra eq 1e15,count) if count gt 0 then u850_merra[result]=!values.f_NAN result=where(v850_merra eq 1e15,count) if count gt 0 then v850_merra[result]=!values.f_NAN result=where(h250_merra eq 1e15,count) if count gt 0 then h250_merra[result]=!values.f_NAN result=where(u250_merra eq 1e15,count) if count gt 0 then u250_merra[result]=!values.f_NAN result=where(v250_merra eq 1e15,count) if count gt 0 then v250_merra[result]=!values.f_NAN result=where(h500_merra eq 1e15,count) if count gt 0 then h500_merra[result]=!values.f_NAN result=where(u500_merra eq 1e15,count) if count gt 0 then u500_merra[result]=!values.f_NAN result=where(v500_merra eq 1e15,count) if count gt 0 then v500_merra[result]=!values.f_NAN result=where(q850_merra eq 1e15,count) if count gt 0 then q850_merra[result]=!values.f_NAN result=where(t850_merra eq 1e15,count) if count gt 0 then t850_merra[result]=!values.f_NAN result=where(u850_merra eq 1e15,count) if count gt 0 then u850_merra[result]=!values.f_NAN result=where(v850_merra eq 1e15,count) if count gt 0 then v850_merra[result]=!values.f_NAN result=where(h850_merra eq 1e15,count) if count gt 0 then h850_merra[result]=!values.f_NAN result=where(h1000_merra eq 1e15,count) if count gt 0 then h1000_merra[result]=!values.f_NAN result=where(u2m_merra eq 1e15,count) if count gt 0 then u2m_merra[result]=!values.f_NAN result=where(v2m_merra eq 1e15,count) if count gt 0 then v2m_merra[result]=!values.f_NAN ;*************** ; Start synoptic 4 panel ;*************** ; int_flag=0 automatically determines the contour intervals ; int_flag=1 use a standard set of contour intervals if delta_degree gt 10.0 then begin int_flag=1 endif else if delta_degree le 10.0 then begin int_flag=0 endif ; Plot dimensions pxdim=1000 pydim=1000 ; Position plots xl=0.05 & xr=0.95 yb=0.1 & yt=0.90 sx=0.05 sy=0.15 numplots_x=2 numplots_y=2 position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos ; Colorbar position cbpos=pos cbpos[*,1]=pos[*,1]-0.04 ;ybottom cbpos[*,3]=pos[*,1]-0.03 ;ytop cbpos[*,0]=pos[*,0]+0.05 cbpos[*,2]=pos[*,2]-0.05 tdy=0.01 p0=plot([0,1],[0,1],/buffer,dimensions=[pxdim,pydim],/nodata,axis_style=4) ; Define the color table levels = 8 ;ctable = COLORTABLE([[200,180,255],[140,170,255],[20,170,255],$ ; [40,255,200],[50,255,50],[255,255,70],[255,165,40],[255,70,70]],ncolors=levels,/TRANSPOSE) ctable = COLORTABLE(['lavender','thistle','light sky blue',$ 'aquamarine','light green','khaki','sandy brown','light coral'],ncolors=levels,/TRANSPOSE) fontsize1=12 map_fontsize=12 continent_color='green' ;*********************************************************** ; Panel 1 - 300 mb Height and 300 mb wind ;*********************************************************** ; Plot location pnum=0 ; Panel title '250 mb Geopotential Height' t0=text( pos[pnum,0]+(pos[pnum,2]-pos[pnum,0])/2.0,$ ;xcoor pos[pnum,3]+tdy,alignment=0.5,$ ;ycoor '250 mb Geopotential Height',font_size=fontsize1) ; Set up the map m0=map('mercator',limit=[llat,llon,ulat,rlon],$ box_axes=1,label_position=0,linestyle='dotted',$ position=pos[pnum,*],/current,font_size=map_fontsize) ; Change the lat/lons separately. m0['Latitudes'].LABEL_ANGLE = 90 m0['Longitudes'].LABEL_ANGLE = 0 ; Compute the wind speed wspd=sqrt(u250_merra^2.+v250_merra^2.) ;meters/second ; Compute the wind direction ;wind_direction = atan(-1.0*u250_merra,-1.0*v250_merra) * (180.0/!pi) ; Levels for color contours of windspeed if int_flag eq 0 then begin dmax=max(wspd) dmin=min(wspd) endif else begin dmax=70.0;180.0 dmin=5.0;9.0 endelse levels2=(indgen(8)*ceil(ceil(dmax-dmin)/7.0))+floor(dmin) ; Filled color contours of windspeed c1=contour(wspd,lon_merra,lat_merra,/overplot,$ c_value=levels2,c_color=ctable,/fill,grid_units='degrees',$ xstyle=1,ystyle=1,font_size=fontsize1) ; Plot the 250mb windbarbs ; Plot a wind barb of u and v components v0=vector(u250_merra,v250_merra,lon_merra,lat_merra,/overplot,$ ;/auto_subsample,$ x_subsample=ceil(delta_degree/2.0),y_subsample=ceil(delta_degree/2.0),$ color='blue',grid_units='degrees',vector_style=1) ; Continents mc=mapcontinents(/countries,color=continent_color,thick=1.5) mc2=mapcontinents(/USA,color=continent_color) ; Contour 300mb geopotential height (120m interval) c0=contour(H250_merra,lon_merra,lat_merra,/overplot,font_size=fontsize1,$ color='black',c_thick=1.5,grid_units='degrees',position=pos[pnum,*]) ; Symbol at track location s0=symbol(lon1,lat1,'circle',sym_color='red',/data,/current,target=m0,sym_size=2) ; colorbar cb0=colorbar(target=c1,position=reform(cbpos[pnum,*]),/border,$ title='250 mb Wind Speed (m/s)',font_size=fontsize1) ;************************************************************ ; Panel 2 - 500 mb Height and Vorticity ;************************************************************ ; Plot location pnum=1 ; Panel title '500 mb Geopotential Height' t0=text( pos[pnum,0]+(pos[pnum,2]-pos[pnum,0])/2.0,$ ;xcoor pos[pnum,3]+tdy,alignment=0.5,$ ;ycoor '500 mb Geopotential Height',color='black',/normal,font_size=fontsize1) ; Set up the map m0=map('mercator',limit=[llat,llon,ulat,rlon],/current,font_size=map_fontsize,$ box_axes=1,label_position=0,linestyle='dotted',position=pos[pnum,*]) ; Change the lat/lons separately. m0['Latitudes'].LABEL_ANGLE = 90 m0['Longitudes'].LABEL_ANGLE = 0 ; Calculate relative vorticity relative_vorticity,lat_merra,lon_merra,u500_merra,v500_merra,theta ; Calculate earth vorticity ; EARTH VORTICITY = 2 x (Rate of Rotation of Earth) x sin (Latitude) earth_vorticity=reverse(2.0*7.27e-5 * sin(lat_merra*!DtoR)) ;radians/sec for j=0,n_elements(lon_merra)-1 do begin if j eq 0 then begin earth_vorticity2d=transpose(earth_vorticity) endif else begin earth_vorticity2d=[earth_vorticity2d,transpose(earth_vorticity)] endelse endfor ; Absolute vorticity is relative vorticity + earth vorticity abs_vort=theta+earth_vorticity2d ; Convert abs_vort from /s to 10^-5/s abs_vort5=abs_vort*100000.0 ; Levels for color contours of vorticity if int_flag eq 0 then begin dmax=max(abs_vort5) dmin=min(abs_vort5) endif else begin dmax=2.0;5.0 dmin=-20.0;-25.0 endelse levels2=indgen(8)*(floor((dmax-dmin)/7.0))+floor(dmin) ; Filled color contours of absolute vorticity c1=contour(abs_vort5,lon_merra,lat_merra,/overplot,xstyle=1,ystyle=1,$ c_value=levels2,c_color=ctable,/fill,grid_units='degrees') ; Plot the 500mb (i=5) windbarbs ; Plot a wind barb of u and v components v0=vector(u500_merra,v500_merra,lon_merra,lat_merra,vector_style=1,/overplot,$ x_subsample=ceil(delta_degree/2.0),y_subsample=ceil(delta_degree/2.0),$ color='blue',grid_units='degrees') ; Continents mc=mapcontinents(/countries,color=continent_color,thick=1.5) mc2=mapcontinents(/USA,color=continent_color) ; Contour 500mb geopotential height c0=contour(H500_merra,lon_merra,lat_merra,/overplot,$ color='black',c_thick=1.5,grid_units='degrees',position=pos[pnum,*]) ; track location s0=symbol(lon1,lat1,'circle',sym_color='red',/data,/current,target=m0,sym_size=2) ; colorbar cb0=colorbar(target=c1,position=reform(cbpos[pnum,*]),/border,$ title='Absolute Vorticity (10^-5/s)',font_size=fontsize1) ;****************************************************** ; Panel 3 - 850 mb Data ;****************************************************** ; Plot location pnum=2 ; Panel title '850 mb Geopotential Height' t0=text( pos[pnum,0]+(pos[pnum,2]-pos[pnum,0])/2.0,$ ;xcoor pos[pnum,3]+tdy,alignment=0.5,$ ;ycoor '850 mb Geopotential Height',color='black',/normal,font_size=fontsize1) ; Set up the map m0=map('mercator',limit=[llat,llon,ulat,rlon],font_size=map_fontsize,$ box_axes=1,label_position=0,linestyle='dotted',position=pos[pnum,*],/current) ; Change the lat/lons separately. m0['Latitudes'].LABEL_ANGLE = 90 m0['Longitudes'].LABEL_ANGLE = 0 ; ; Compute relative humidity ; t=t850_merra ;K ; t0=273.15 ;K ; q=q850_merra ; p=850.0*100.0 ;Pa ; s=size(t,/dimensions) ; rh=make_array(/float,s[0],s[1],value=!values.f_NAN) ; ;result=where( finite(t,/nan) or finite(q,/nan) ,countnan) ; result=where( finite(t) and finite(q),count) ; if count gt 0 then begin ; rh[result]=0.263*p*q[result]*(exp( (17.67*(t[result]-t0)) / (t[result]-29.65) ) )^(-1) ; endif t850_merra=t850_merra-273.15 ; Levels for color contours of temperature if int_flag eq 0 then begin result=where(finite(t850_merra)) dmax=max(t850_merra[result]) dmin=min(t850_merra[result]) levels2=(indgen(8)*round(round(dmax-dmin)/7.0))+floor(dmin) endif else begin dmax=24.0 dmin=-18.0 levels2=(indgen(8)*ceil(ceil(dmax-dmin)/7.0))+floor(dmin) endelse ; Filled color contours of temp c1=contour(t850_merra,lon_merra,lat_merra,/overplot,xstyle=1,ystyle=1,$ c_value=levels2,c_color=ctable,/fill,grid_units='degrees',font_size=fontsize1) ; Plot the 850mb windbarbs ; Plot a wind barb of u and v components v0=vector(u850_merra,v850_merra,lon_merra,lat_merra,vector_style=1,/overplot,$ x_subsample=ceil(delta_degree/2.0),y_subsample=ceil(delta_degree/2.0),$ color='blue',grid_units='degrees') ; Continents mc=mapcontinents(/countries,color=continent_color,thick=1.5) mc2=mapcontinents(/USA,color=continent_color) ; Contour 850mb geopotential height c0=contour(H850_merra,lon_merra,lat_merra,/overplot,font_size=fontsize1,$ color='black',c_thick=1.5,grid_units='degrees',position=pos[pnum,*]) ; Contour temp lines ;t850_merra=t850_merra-273.15 ;convert from kelvin to celcius ;c2=contour(t850_merra,lon_merra,lat_merra,/overplot,c_linestyle=2,$ ; color='black',c_thick=2,grid_units='degrees',position=pos[pnum,*]) ; track location s0=symbol(lon1,lat1,'circle',sym_color='red',/data,/current,target=m0,sym_size=2) ; Colorbar cb0=colorbar(target=c1,position=reform(cbpos[pnum,*]),/border,$ title='850 mb Temperature (C)',font_size=fontsize1) ;************************************************************* ; Panel 4 - 1000 mb pressure and winds ;************************************************************* ; Plot location pnum=3 ; Panel title '1000 mb Geopotential Height' t0=text( pos[pnum,0]+(pos[pnum,2]-pos[pnum,0])/2.0,$ ;xcoor pos[pnum,3]+tdy,alignment=0.5,$ ;ycoor '1000 mb Geopotential Height',color='black',/normal,font_size=fontsize1) ; Set up the map m0=map('mercator',limit=[llat,llon,ulat,rlon],/current,font_size=map_fontsize,$ box_axes=1,label_position=0,linestyle='dotted',position=pos[pnum,*]) ; Change the lat/lons separately. m0['Latitudes'].LABEL_ANGLE = 90 m0['Longitudes'].LABEL_ANGLE = 0 ; Contour lines delta_h=(h500_merra-h1000_merra)/1000.0 ; Levels for color contours of deltah if int_flag eq 0 then begin dmax=max(delta_h) dmax=string(dmax,format='(F3.1)') dmax=float(dmax) dmin=min(delta_h) dmin=string(dmin,format='(F3.1)') dmin=float(dmin) xx=(dmax-dmin)/7.0 xx=string(xx,format='(F4.2)') xx=float(xx) levels2=findgen(8)*xx+dmin endif else begin dmax=5.8 dmin=5.2 levels2=indgen(8)*((dmax-dmin)/7.0)+floor(dmin) levels2=[5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9] endelse ; Filled color contours thickness c1=contour(delta_h,lon_merra,lat_merra,/overplot,c_value=levels2,$ c_color=ctable,/fill,grid_units='degrees',xstyle=1,ystyle=1,font_size=fontsize1) ; Plot the 2meter windbarbs ; Plot a wind barb of u and v components v0=vector(u2m_merra,v2m_merra,lon_merra,lat_merra,vector_style=1,/overplot,$ x_subsample=ceil(delta_degree/2.0),y_subsample=ceil(delta_degree/2.0),$ color='blue',grid_units='degrees') ; Continents mc=mapcontinents(/countries,color=continent_color,thick=1.5) mc2=mapcontinents(/USA,color=continent_color) ; Contour lines - 1000mb c0=contour(h1000_merra,lon_merra,lat_merra,/overplot,$ color='black',c_thick=1.5,grid_units='degrees',position=pos[pnum,*]) ; track location s0=symbol(lon1,lat1,'circle',sym_color='red',/data,/current,target=m0,sym_size=2) ; Colorbar cb0=colorbar(target=c1,position=reform(cbpos[pnum,*]),/border,$ title='500mb-1000mb Thickness (km)',font_size=fontsize1) ; Titles if program eq 'socrates' then cap_program='SOCRATES' t5=text(0.09,0.95,cap_program,font_size=fontsize1) t0=text(0.2,0.95,string(yy_track,format='(I4)')+$ string(mm_track,format='(I02)')+string(dd_track,format='(I02)')+' '+$ string(hh_track,format='(I02)')+':'+string(mi_track,format='(I02)')+':'+$ string(ss_track,format='(I02)'),color='black',font_size=fontsize1) t1=text(0.09,0.97,'MERRA-2',font_size=fontsize1,font_name='Helvetica') t1=text(0.2,0.97,string(yy_merra,format='(I4)')+$ string(mm_merra,format='(I02)')+string(dd_merra,format='(I02)')+' '+$ string(hh_merra,format='(I02)')+':'+string(mi_merra,format='(I02)')+':'+$ string(ss_merra,format='(I02)'),color='black',/normal,font_size=fontsize1) t2=text(0.38,0.95,string(lat1,format='(F6.2)')+' lat '+string(lon1,format='(F6.2)')+' lon',$ font_size=fontsize1) ;t3=text(0.45,0.95,string(lat_merra[lat_idx],format='(F9.4)')+'lat '+$ ; string(lon_merra[lon_idx],format='(F9.4)')+'lon',$ ; color='black',/normal) ;merra_hour_str=string(hh_merra,format='(I02)')+'h' imagedir='/uufs/chpc.utah.edu/common/home/mace-group2/'+program+'/synoptic/' imagename='synoptic.'+program+'.'+date+'.'+shour+'h.'+$ string(delta_degree,format='(I03)')+'deg.png' m0.save,imagedir+imagename,height=pydim print,imagename ;endif ;close enough ;endfor ;loop through airplane track stop ; URLS ;url = 'http://eosdap.hdfgroup.uiuc.edu:8080/opendap/data/NASAFILES/hdf4/AIRS.2008.10.27.L3.RetStd001.v5.2.2.0.G08303124144.hdf' ;url = 'http://goldsmr5.sci.gsfc.nasa.gov:80/opendap/MERRA2/M2T3NVASM.5.12.4/'+$ ; '2014/05/MERRA2_400.tavg3_3d_asm_Nv.20140501.nc4?'+$ ; 'T[0:1:7][0:1:71][0:1:360][0:1:575],lat[0:1:360],lon[0:1:575]' ;url = 'http://goldsmr5.sci.gsfc.nasa.gov:80/opendap/MERRA2/M2T3NVASM.5.12.4/'+$ ; syear+'/'+smonth+'/MERRA2_400.tavg3_3d_asm_Nv.'+syear+smonth+sday+'.nc4?'+$ ; 'T[0:1:7][0:1:71][0:1:360][0:1:575],lat[0:1:360],lon[0:1:575],time[0:1:7]' end