pro plot_gpm_figure1

;******************
;  Input
;******************
gpm_dir='/uufs/chpc.utah.edu/common/home/mace-group4/gpm/gpm1.gesdisc.eosdis.nasa.gov/data/'

fdir_ku='/GPM_L2/GPM_2AKu.05/2018/001/'
fname_ku='2A.GPM.Ku.V7-20170308.20180101-S010928-E024202.021833.V05A.HDF5'

radar='ku'
radar='ka'
if radar eq 'ku' then begin
  radarname='KuPR'
  radarfreq='13.6 GHz'
  radarswath='245 km'
endif else if radar eq 'ka' then begin
  radarname='KaPR'
  radarfreq='35.5 GHz'
  radarswath='120 km'
endif

;  Set up the map
;ulat=-10.0 & llat=-70.0
;llon=70 & rlon=180+45  ;0 to 360 degrees

;  Set up the map
ulat=-45.0 & llat=-70.0
llon=120.0 & rlon=180.0  ;0 to 360 degrees

;  All data
;ulat=90.0 & llat=-90.0
;llon=0 & rlon=360  ;0 to 360 degrees

;Swath S1 has channels 1-9: 9channels
;10V 10H 19V 19H 23V 37V 37H 89V 89H.
;Swath S2 has channels 10-13: 4channels
;166V 166H 183+/-3V 183+/-8V.
;chan=7  ;channel
;channels=['10V','10H','19V','19H','23V','37V','37H','89V','89H',$
;          '166V','166H','183+/-3V','183+/-8V']
          
;******************
;  Read data from hdf5 file
;******************

;presult=h5_parse(fdir+fname)

;  Open an existing hdf5 file
file_id=h5f_open(gpm_dir+fdir_ku+fname_ku)
ns_id=h5g_open(file_id,'NS')  
dset_id=h5d_open(ns_id,'Latitude')
lat=h5d_read(dset_id)
h5d_close,dset_id

dset_id=h5d_open(ns_id,'Longitude')
lon=h5d_read(dset_id)
h5d_close,dset_id

slv_id=h5g_open(ns_id,'SLV')
dset_id=h5d_open(slv_id,'zFactorCorrected')
dbz=h5d_read(dset_id)
h5d_close,dset_id
h5g_close,slv_id

time_id=h5g_open(ns_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,ns_id
h5f_close,file_id

;****************
;  nbin,nray,nscan
;****************

s=size(dbz,/dimensions)
nbin=s[0]  ;176
nray=s[1]  ;49
nscan=s[2] ;7937

;*****************
;  Calculate time
;*****************

msecond=msecond*1e-3
second=second+msecond
julian_day=julday(month,day,year,hour,minute,second)

;  loop through the range bins to see where the data is
for i=0,nbin-1 do begin
  dbz1=reform(dbz[i,*,*])
  result=where(dbz1 gt -9999,count)
  if count le 0 then begin
    print,i,count
  endif else if count gt 0 then begin
    print,i,count,min(dbz1[result]),max(dbz1[result])
  endif
endfor

;*****************************
;  Subset the data in the region of interest box          
;*****************************

;if rlon le 180 then begin
;  idx=where(lon ge llon-1.0 and lon le rlon+1.0 and $
;            lat ge llat-1.0 and lat le ulat+1.0,count)
;  lon=lon[idx]
;  lat=lat[idx]
;  julian_day=julian_day[idx]
;  stop
;endif else if rlon gt 180 then begin
  ;  Calculate longitude as 0 to 360
  colon=lon
  result=where(colon lt 0,count)
  if count gt 0 then colon[result]=360.0+colon[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
  ;  Subset data within colon bounds
  idx=where(colon ge lcolon-1.0 and colon le rcolon+1.0 and $
            lat ge llat-1.0 and lat le ulat+1.0,cidx)           
  julian_day=julian_day[idx]
  colon=colon[idx]
  lon=lon[idx]
  lat=lat[idx]
  
  dbz_sub=make_array(nbin,cidx)
  ;  loop through the range bins to see where the data is
  for i=0,nbin-1 do begin
    dbz1=reform(dbz[i,*,*])
    dbz_sub[i,*]=dbz1[idx]
    dbz2=dbz_sub[i,*]
    result=where(dbz2 gt -9999,count)
    if count le 0 then begin
      print,i,count
    endif else if count gt 0 then begin
      print,i,count,min(dbz2[result]),max(dbz2[result])
    endif       
  endfor  

 
  ;convert these back to longitude
  rlon=rcolon
  if rlon gt 180 then rlon=rlon-360.0
  if llon gt 180 then llon=llon-360.0
;endif
print,llon,rlon
print,lcolon,rcolon

;*********************
;  Subset out subsatellite point
;*********************

lat1=lat[24,*]
lon1=lon[24,*]
dbz1=dbz_sub[*,24,*]
stop
;*************************************************
;  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_s1[*,0,0]
  result=where(var gt -9999.90,count)  
  dmin=min(var[result])
  dmax=max(var[result])
  p2=plot(julian_day_s1,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_1bgmi.png'
  p0.save,imagename,height=pydim
endif  ;end of plot to check

;***************************
;  Make a plot of this subset
;***************************

;  Set up plot size
pxdim=900 & pydim=600
;  Position the plots
xl=0.07 & xr=0.80
yb=0.10 & yt=0.80
sx=0.07
sy=0.05
numplots_x=1
numplots_y=1
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

;********************************
;  My 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])]
mycbtable=mytable[0:top_color,*]
revcbtable=revmytable[0:top_color,*]

;*******************************************
;  Set up the color tables
;*******************************************

dmax=300;340.0
dmin=200;180;160.0
mp=220.0  ; 243 cirrus start here
xp=fix(256*((dmin-mp)/(dmin-dmax)))

common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr

R_new=indgen(256)
G_new=indgen(256)
B_new=indgen(256)

loadct, 15, ncolors=xp
R_new(0:xp-1)=R_curr(0:xp-1)
G_new(0:xp-1)=G_curr(0:xp-1)
B_new(0:xp-1)=B_curr(0:xp-1)

loadct, 0, ncolors=255-xp+1
gamma_ct, 1.4, /intensity
R_new(xp:255)=reverse(R_curr(0:255-xp))
G_new(xp:255)=reverse(G_curr(0:255-xp))
B_new(xp:255)=reverse(B_curr(0:255-xp))

;  IR table
ir_table=make_array(/integer,3,256)
ir_table[0,*]=R_new
ir_table[1,*]=G_new
ir_table[2,*]=B_new

;  16 Color tablle
;  Top is the last color to scale 256 colors, 0-255
top_color=252
;  Colortable  0-252  253=white
ref_table=colortable(38,ncolors=254)
;254=hot pink               ;gray=255
ref_table=[ref_table,transpose([238,18,137]),transpose([150,150,150])]
ref_cbtable=ref_table[0:top_color,*]

;********Color table end********

;  Strings
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])


;  Calculate regular grid array
delta=0.1
ynum=(ulat-llat)/delta
grid_lat=findgen(ynum)*delta+llat
;if rlon-llon lt 0 then begin
  xnum=(rcolon-lcolon)/delta
  grid_colon=findgen(xnum)*delta+lcolon  
  grid_lon=grid_colon
  result=where(grid_lon gt 180,count)
  if count gt 0 then grid_lon[result]=grid_lon[result]-360.0
;endif else if rlon-llon gt 0 then begin  
;  xnum=(rlon-llon)/delta
;  grid_lon=findgen(xnum)*delta+llon
;endif

;  Triangulate the data
print,'start qhull'
;qhull,lon,lat,triangles,/delaunay;,sphere=s
qhull,colon,lat,triangles,/delaunay,sphere=s
print,'end qhull, start gridding'

var=dbz_sub[166,*]
result=where(var lt -9999,count)
var[result]=-7777
;grid_tb=griddata(lon,lat,tb,xout=grid_lon,yout=grid_lat,$
grid_var=griddata(colon,lat,var,xout=grid_colon,yout=grid_lat,$  
  /grid,/degrees,/sphere,$
  ;/KRIGING,min_points=10,sectors=8,empty_sectors=4,$
  ;method='NearestNeighbor',$
  /inverse_distance,min_points=10,sectors=8,empty_sectors=3,$;max_per_sector=10,$  
  ;/inverse_distance,empty_sectors=1,max_per_sector=1,sectors=1,min_points=1,$
  missing=-9999,$
  triangles=triangles);,$
;missing=!values.f_nan)
print,'end gridding'

;  Bytscal the data
result=where(grid_var gt 0)
dmin=min(grid_var[result])
dmax=max(grid_var[result])
data_image=bytscl(grid_var,top=top_color,min=dmin,max=dmax)
result=where(grid_var gt -8000 and grid_var le -7000,count)
if count gt 0 then data_image[result]=255
result=where(grid_var gt -10000 and grid_var le -9000,count)
if count gt 0 then data_image[result]=253

p0=image(data_image,grid_colon,grid_lat,limit=[llat,lcolon,ulat,rcolon],$  
  /buffer,dimensions=[pxdim,pydim],grid_units='degrees',$
  center_longitude=(rcolon-lcolon)/2.0,$
  rgb_table=mytable,$;,min_value=dmin,max_value=dmax,$
  map_projection='Mercator',$  
  position=pos[pnum,*])
  
p1=mapcontinents(color='black')
  p0.mapgrid.box_axes=1  
  p0.mapgrid.box_color='black'
  p0.mapgrid.color='black'
  p0.mapgrid.label_color='black'
  p0.mapgrid.linestyle='dot'  
  p0.mapgrid.label_position=0

  c0=colorbar($
    rgb_table=mycbtable,$
    range=[dmin,dmax],$
    orientation=1,position= reform(cbpos[pnum,*]),font_size=10,$
    tickdir=1,border_on=1,title='DbZ')
print,dmin,dmax 
stime=min(julian_day)
etime=max(julian_day)
caldat,stime,smm,sdd,syy,shh,smi,sss
caldat,etime,emm,edd,eyy,ehh,emi,ess

sdate_str=string(syy,format='(I4)')+string(smm,format='(I02)')+string(sdd,format='(I02)')
stime_str=string(shh,format='(I02)')+string(smi,format='(I02)')+string(sss,format='(I02)')  
edate_str=string(eyy,format='(I4)')+string(emm,format='(I02)')+string(edd,format='(I02)')
etime_str=string(ehh,format='(I02)')+string(emi,format='(I02)')+string(ess,format='(I02)')  

stime_cstr=string(shh,format='(I02)')+':'+string(smi,format='(I02)')+':'+string(sss,format='(I02)') 
etime_cstr=string(ehh,format='(I02)')+':'+string(emi,format='(I02)')+':'+string(ess,format='(I02)')
datastream='GPM_2ADPR.05'
;t3=text(0.1,0.90,'Channel: '+channels[chan])
t1=text(0.1,0.96,datastream+' '+sdate_str+' '+stime_cstr+' - '+edate_str+' '+etime_cstr)
t2=text(0.1,0.93,'Dual-Frequency Precipitation Radar (DPR)')
t3=text(0.1,0.90,'Reflectivity Factor with attenuation correction (dBZ) at height bin 166') 
imagename='map.'+datastream+'.dbz.'+$
  sdate_str+stime_str+'.'+edate_str+etime_str+'.NS1.png'
    
p0.save,imagename,height=pydim
stop
i0=image(grid_tb,grid_lon,grid_lat,/buffer,$
         dimensions=[pxdim,pydim],grid_units='degrees',$
         limit=[llat,llon,ulat,rlon],$
         rgb_table=mytable,$
         map_projection='Mercator',$
         position=pos[pnum,*])
i1=mapcontinents(color='white')
i0.mapgrid.box_axes=1
       ;;i0.mapgrid.label_position=0
       i0.mapgrid.box_color='black'
       i0.mapgrid.color='black'
       i0.mapgrid.label_color='black'
       i0.mapgrid.linestyle='dot'

p0.save,'map_tb.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