pro plot_ver7_1z,fdat=fdat,filename=filename,flashes=flashes,plotarea=plotarea,$
      cross=cross,locations=locations,varlists=varlists,checkcross=checkcross,noshow=noshow,$
      nocross=nocross,ifterrain=ifterrain,outfile=outfile
;    
;   f:   the structure containing  the 1Z data   
;   filename: the filename have to be read to get the 1Z data 
;   flashes:  if you want to overplot the flash locations
;   plotarea: the area select to plot, if not set, will have to select inter actively
;   cross:    the cross section locations, multiple crosses is ok
;   locations: the variable holds the choice of area and cross and return to the calling program
;   varlists: maximum 9 variables may be called from following
;       'NEARSURFZ','NEARSURFRAIN','PCT85','PCT37','V19','CH4'
;       '2A12RAIN','STRATCONV','2A12CLDICE','2A12CLDWATER'
;       '2A12PRECIPWATER','2A12PRECIPICE','STORMHT'
;   checkcross: check the cross section location and see the details
;   noshow:  knob to show the ps files in ggv
;          
;   this code plots the variables in 1Z file
;   with horizontal and vertical sections.
;            Chuntao Liu
;                    10. 2005  
;
;

;goto,jump3
set_plot,'x'
!p.multi=0
;
; read in the data 
;
   if(n_elements(fdat) le 0) then begin
      if(n_elements(filename) le 0) then begin
        filename=dialog_pickfile(filter='1Z*.HDF*')
      endif    
      if(strpos(filename,'gz') ge 0) then begin
        spawn,'cp '+filename+' .'
        ftmp=strsplit(filename,'/',/extract) 
        ftmp=ftmp[n_elements(ftmp)]
        spawn,'gunzip '+ftmp 
        ftmp=strsplit(ftmp,'.',/extract)
        filename=ftmp[0]+'.'+ftmp[1]+'.'+ftmp[2]+'.6.HDF'
      endif
      read_pf_level1_v7_hdf,filename,fdat
   endif 
pr_lon=fdat.pr.lon
tmi_lon=fdat.tmi.lonhi
virs_lon=fdat.virs.lon
lis_lon=fdat.lis.flhlon
if(n_elements(plotarea) eq 4) then begin
if (max(plotarea[[1,3]]) ge 180 ) then begin
 id=where(pr_lon lt 0)
 pr_lon[id]=360+pr_lon[id]
 id=where(tmi_lon lt 0)
 tmi_lon[id]=360+tmi_lon[id]
 id=where(virs_lon lt 0)
 virs_lon[id]=360+virs_lon[id]
 id=where(lis_lon lt 0)
 lis_lon[id]=360+lis_lon[id]
endif
endif

   f=fdat
   if(n_elements(varlists) le 0) then begin
     varlists=['PCT85','PCT37','CH4','NEARSURFZ','NEARSURFRAIN','2A12RAIN'] 
   endif
   ez_color,1,colors,ncolor=60 
   orbit=stringit(long(f.pr.orbit))
;
; define the plotting area
;
 answer=''
 if(n_elements(plotarea) lt 4) then begin
  plotarea=[-36.,-180.,36,180]
  zoomin:
   id_pr=where(pr_lon ge plotarea[1] and pr_lon le plotarea[3] and $
               f.pr.lat ge plotarea[0] and f.pr.lat le plotarea[2])
   jutc=julday(f.pr.month[id_pr/49],f.pr.day[id_pr/49],f.pr.year[id_pr/49],f.pr.hour[id_pr/49],$
      f.pr.minute[id_pr/49],f.pr.second[id_pr/49])
   tjutc=jutc[sort(jutc)]   
   caldat,tjutc,tmonth,tday,tyear,thour,tminute,tsecond 
   date=stringit(long(mean(tyear[0])))+'-'+stringit(long(mean(tmonth[0])))+'-'+$
        stringit(long(mean(tday[0])))+' '+stringit(long(mean(thour[0])))+':'+$
        stringit(long(mean(tminute[0])))+$
        ' UTC '
   print,plotarea
   map_set,limit=plotarea,/hires,/label,ymargin=[6,2],title=date+orbit+' 2A25 near surface reflectivity '
;    plot_pixels,f.pr.nearsurfz,pr_lon,f.pr.lat
    contour,f.pr.nearsurfz,pr_lon,f.pr.lat,levels=indgen(10)*6,$
            c_color=colors[indgen(10)*6],/overplot,c_label=intarr(10)+1,/fill
    contour,f.tmi.pct85,tmi_lon,f.tmi.lathi,$
        levels=[250],/overplot      
   map_set,limit=plotarea,/hires,/label,ymargin=[6,2],title=date+orbit+' 2A25 near surface reflectivity ',/noerase
   colorbar,indgen(11)*6,colors[indgen(10)*6],/ivestyle,/noframe,format='(i3)',ysize=0.01
   read,'Zoom in? (y/n) ',answer
    if answer eq 'y' then begin &$
        print,'Click on lower left corner of desired zoom'
        cursor,x1,y1
        print,'x1, y1 =',x1,y1
        read,'Press return',answer
        print,'Click on upper right corner of desired zoom'
        cursor,x2,y2
        print,'x2, y2 =',x2,y2
        plotarea[0]=min([y1,y2])
        plotarea[1]=min([x1,x2])
        plotarea[2]=max([y1,y2])
        plotarea[3]=max([x1,x2])
        goto,zoomin
    endif
 endif
;
; define the cross section location
;
 if(n_elements(nocross) gt 0) then cross=0
 if(n_elements(cross) lt 4 and n_elements(nocross) eq 0) then begin
  cross=0
  crossin:
  read,'want to drow a cross section? (y/n)',answer
    if answer eq 'y' then begin &$
        print,'Click on the first point'
        cursor,x1,y1
        oplot,[x1],[y1],psym=1,symsize=3.0,color=55
;        print,'x1, y1 =',x1,y1
        read,'Press return',answer
        print,'Click on the second point'
        cursor,x2,y2
        oplot,[x2],[y2],psym=1,symsize=3.0,color=55
        oplot,[x1,x2],[y1,y2],thick=3
;        print,'x2, y2 =',x2,y2
        if(n_elements(cross) ge 4) then cross=[[cross],[x1,y1,x2,y2]]
        if(n_elements(cross) lt 4) then cross=[x1,y1,x2,y2]
        goto,crossin
    endif
 endif 
   id_pr=where(pr_lon ge plotarea[1] and pr_lon le plotarea[3] and $
               f.pr.lat ge plotarea[0] and f.pr.lat le plotarea[2])
   jutc=julday(f.pr.month[id_pr/49],f.pr.day[id_pr/49],f.pr.year[id_pr/49],f.pr.hour[id_pr/49],$
               f.pr.minute[id_pr/49],f.pr.second[id_pr/49])
   tjutc=jutc[sort(jutc)]   
   caldat,tjutc,tmonth,tday,tyear,thour,tminute,tsecond 
   date=stringit(long(mean(tyear[0])))+'-'+stringit(long(mean(tmonth[0])))+'-'+$
        stringit(long(mean(tday[0])))+' '+stringit(long(mean(thour[0])))+':'+$
        stringit(long(mean(tminute[0])))+':'+stringit(long(mean(tsecond[0])))+$
        ' UTC '
 locations={plotarea:plotarea,cross:cross}
;
;  start to plot out the results
;
  ;
  ; cross section 
  ;
jump3:
  if(n_elements(cross) ge 4) then begin

     for k=0,n_elements(cross[0,*])-1 do begin
;     window,k
       set_plot,'ps'
       !p.font=6
        ez_color,1,colors,ncolor=60
       if(n_elements(outfile) eq 0) then outfile=orbit
       outcross=outfile+'_cross_section_'+stringit(k)+'.ps'
       device,/color,filename=outcross,xsize=25,ysize=15
       id_pr=where(pr_lon ge min(cross[[0,2],k])-0.02 and pr_lon le max(cross[[0,2],k])+0.02 and $
               f.pr.lat ge  min(cross[[1,3],k])-0.02 and f.pr.lat le  max(cross[[1,3],k])+0.02)
       distance=cos((cross[1,k]+cross[3,k])/2.0*!pi/180.0)*111.0*$
              sqrt((cross[3,k]-cross[1,k])^2.0+(cross[2,k]-cross[0,k])^2.0) 
        print,distance
        dlatlon=cos((cross[1,k]+cross[3,k])/2.0*!pi/180.0)*111.0
        dis=abs((cross[2,k]-cross[0,k])*(cross[1,k]-f.pr.lat[id_pr])-$
                 (cross[3,k]-cross[1,k])*(cross[0,k]-pr_lon[id_pr])) $
           /sqrt((cross[0,k]-cross[2,k])^2.0+(cross[1,k]-cross[3,k])^2.0)*dlatlon
         idc=where(dis lt 3.5)               ; choose the points close to the cross line
         vdis=dis[idc]
         id=id_pr[idc]
           dis1=sqrt((cross[0,k]-pr_lon[id])^2.0+(cross[1,k]-f.pr.lat[id])^2.0$
               -(vdis/dlatlon)^2.0)*dlatlon  
          id1=sort(dis1) 
          id=id[id1]
          dis1=dis1[id1]                    ; distance from the first point

       print,'min dis1:',min(dis1)

       if(n_elements(checkcross) gt 0) then begin             ; check if the cross section is in the right order
        !p.multi=[0,1,2]
         map_set,limit=[min(cross[[1,3],k])-0.1,min(cross[[0,2],k])-0.1,$
             max(cross[[1,3],k])+.1,max(cross[[0,2],k])+.1],/label
          id_pr=where(pr_lon[0,*] ge plotarea[1]-3 and pr_lon[0,*] le plotarea[3]+3)
        lon=pr_lon[*,min(id_pr):max(id_pr)]
        lat=f.pr.lat[*,min(id_pr):max(id_pr)]
        data=f.pr.nearsurfz[*,min(id_pr):max(id_pr)]
        plot_pixels,data,lon,lat

         oplot,cross[[0,2],k],cross[[1,3],k],psym=1,thick=2
         oplot,cross[[0,2],k],cross[[1,3],k],thick=3
         oplot,pr_lon[id],f.pr.lat[id],psym=2
         oplot,pr_lon[id],f.pr.lat[id],line=0
         !p.multi=[0,1,2]
         map_set,limit=[min(cross[[1,3],k])-0.1,min(cross[[0,2],k])-0.1,$
             max(cross[[1,3],k])+.1,max(cross[[0,2],k])+.1],/label,/noerase
       endif
                                                             ;  start to plot cross section
        plot,[0,distance],[0,20],psym=3,xtitle='distance(km)',ytitle='Altitude (km), rainrate (10 mm/hr)',$
             ymargin=[8,2],xstyle=1,title=orbit+' '+date,xmargin=[10,10]
        xpicsize=(dis1-[0,dis1[0:n_elements(id)-2]])>0
        pct85=fltarr(n_elements(id))
        pct37=fltarr(n_elements(id))
        ch4=fltarr(n_elements(id))
        dbz_cross=fltarr(n_elements(id),80)
        rain_2a25=fltarr(n_elements(id))
        rain_2a12=fltarr(n_elements(id)) 
        xlon=fltarr(n_elements(id)) 
        xlat=fltarr(n_elements(id)) 
        elev=fltarr(n_elements(id)) 
        read_sds,'/uufs/chpc.utah.edu/common/home/ez_grp/data/trmm_ana/pf_v6.1/soft/terrain.HDF',terr

        for i=0,n_elements(id)-1 do begin
          rain_2a25[i]=f.pr.nearsurfrain[id[i]]/10.0
          xlon[i]=pr_lon[id[i]]
          xlat[i]=f.pr.lat[id[i]]
          elev[i]=terr.terr[(xlon[i]+180.0)/0.16666666,n_elements(terr.terr[0,*])-(xlat[i]+90)/0.16666666]/1000.0
;          print,xlon[i],xlat[i],terr.lon[(xlon[i]+180.0)/0.16666666,(xlat[i]+90)/0.16666666],terr.lat[(xlon[i]+180.0)/0.16666666,1080-(xlat[i]+90)/0.16666666],elev[i]
          ch4[i]=f.virs.ch4[id[i]]
          rain_2a12[i]=f.tmi.rain[f.colo.hi[id[i]]]/10.0
          pct85[i]=f.tmi.pct85[f.colo.hi[id[i]]]
          pct37[i]=f.tmi.pct37[f.colo.hi[id[i]]/2]
          did=where(id[i] eq f.ray+f.scan*49)
;          print,xpicsize[i]
          if(xpicsize[i] gt 0 and did[0] ge 0 and xpicsize[i] lt 7.0) then begin
            dbz=reform(f.pr_dbz[did[0],*])/100.0
            dbz_cross[i,*]=dbz
            z=(80-indgen(80))/4.0+0.125
            x=intarr(80)+dis1[i]-xpicsize[i]*0.5
            did=where(dbz gt 0)
            if(did[0] ge 0) then plot_pixels_color,dbz[did],x[did],z[did],$       ; reflectivity
                 levels=indgen(61),color=indgen(60),pix_size=[xpicsize[i]/2.0,0.125]
          endif
        endfor
        if(n_elements(ifterrain) gt 0) then oplot,dis1-xpicsize*0.5,elev,line=1,thick=16,color=colors[10]
        oplot,dis1-xpicsize*0.5,rain_2a25,line=0,thick=2
        oplot,dis1-xpicsize*0.5,rain_2a12,line=2,thick=2
        xyouts,0,-1.5,'('+strcompress(string(cross[0,k],format='(f8.2)'))+$
               ','+strcompress(string(cross[1,k],format='(f7.2)'))+')' 
        xyouts,distance*0.85,-1.5,'('+strcompress(string(cross[2,k],format='(f8.2)'))+$
               ','+strcompress(string(cross[3,k],format='(f7.2)'))+')' 
        axis,yaxis=1,yrange=[0,300],ytitle='Brightness Temperature (K)',$
        ytickv=[75,150,225],yticks=3,ticklen=1,/save
        oplot,dis1-xpicsize*0.5,pct85,thick=5
        oplot,dis1-xpicsize*0.5,pct37,thick=5,line=2
        colorbar,indgen(61),indgen(60),levelind=indgen(10)*6,format='(i2)',/ivestyle,/noframe,ysize=0.01
       device,/close
       if(n_elements(noshow) eq 0) then spawn,'ggv '+outcross+'&'
     endfor
  endif
;
; panels
;
 set_plot,'ps'
 ez_color,1,colors,ncolor=60
 if(n_elements(outfile) eq 0) then outfile=orbit
 device,/color,filename=outfile+'_panel.ps',xsize=30,ysize=20
 !p.font=6
 !p.charsize=1.0
 
; varlists=['PCT85','PCT37','V19','2A12RAIN','NEARSURFRAIN','NEARSURFZ',$
;       'CH4','STRATCONV','STORMHT']
; varlists=['PCT85','PCT37','V19','2A12RAIN','NEARSURFRAIN','NEARSURFZ']
; varlists=[ '2A12RAIN','STRATCONV','2A12CLDICE','2A12CLDWATER']
; varlists=['PCT85', 'PCT37','V19']
; varlists=['PCT85']

 ix=1
 iy=1
 xsize=15 & ysize=15 & xmargin=8
 kw=[0]
 if(n_elements(varlists) gt 1 and n_elements(varlists) le 4) then begin      ; 4 panels
  ix=2 & iy=2 & kw=[0,3,2,1] 
   xsize=25 & ysize=25 & xmargin=8
 endif 
 if(n_elements(varlists) gt 4 and n_elements(varlists) le 6) then begin      ; 6 panels
  ix=3 & iy=2 & kw=[0,5,4,3,2,1] 
   xsize=30 & ysize=20 & xmargin=16
 endif 
 if(n_elements(varlists) gt 6 and n_elements(varlists) le 9) then begin      ; 9 panels
  ix=3 & iy=3 & kw=[0,8,7,6,5,4,3,2,1] 
   xsize=30 & ysize=30 & xmargin=16
 endif 
 set_plot,'ps'
 ez_color,1,colors,ncolor=60
 if(n_elements(outfile) eq 0) then outfile=orbit
 device,/color,filename=outfile+'_panel.ps',xsize=xsize,ysize=ysize
 !p.font=6

 
   for k=0,n_elements(varlists)-1 do begin
     if(varlists[k] eq 'PCT85') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       data=(300-f.tmi.pct85[*,min(id_tmi):max(id_tmi)])/2.5
       levels=300-indgen(61)*2.5
       colors=indgen(60)
       levelind=indgen(10)*6
       tic=' 85GHz PCT (K)'
       format='(i3)'
     endif     
    if(varlists[k] eq 'PCT37') then begin
       id_tmi=where(f.tmi.lonlo[0,*] ge plotarea[1]-5 and f.tmi.lonlo[0,*] le plotarea[3]+5)
       lon=f.tmi.lonlo[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.latlo[*,min(id_tmi):max(id_tmi)]
       data=(315-f.tmi.pct37[*,min(id_tmi):max(id_tmi)])/1.5
       levels=315-indgen(61)*1.5
       colors=indgen(60)
       levelind=indgen(10)*6
       tic=' 37GHz PCT (K)'
       format='(i3)'
     endif
    if(varlists[k] eq 'V19') then begin
       id_tmi=where(f.tmi.lonlo[0,*] ge plotarea[1]-5 and f.tmi.lonlo[0,*] le plotarea[3]+5)
       lon=f.tmi.lonlo[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.latlo[*,min(id_tmi):max(id_tmi)]
       data=61-(324-f.tmi.V19[*,min(id_tmi):max(id_tmi)])/2.2
       levels=324-indgen(61)*2.2
       colors=reverse(indgen(60))
       levelind=indgen(10)*6
       tic='19GHz Vertical (K)'
       format='(i3)'
     endif
    if(varlists[k] eq 'H19') then begin
       id_tmi=where(f.tmi.lonlo[0,*] ge plotarea[1]-5 and f.tmi.lonlo[0,*] le plotarea[3]+5)
       lon=f.tmi.lonlo[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.latlo[*,min(id_tmi):max(id_tmi)]
       data=61-(324-f.tmi.H19[*,min(id_tmi):max(id_tmi)])/2.2
       levels=324-indgen(61)*2.2
       colors=reverse(indgen(60))
       levelind=indgen(10)*6
       tic='19GHz Horizontal (K)'
       format='(i3)'
     endif


    if(varlists[k] eq 'CH4') then begin
       id_virs=where(virs_lon[0,*] ge plotarea[1]-5 and virs_lon[0,*] le plotarea[3]+5,vcount)
      if(vcount gt 0) then begin
       lon=virs_lon[*,min(id_virs):max(id_virs)]
       lat=f.virs.lat[*,min(id_virs):max(id_virs)]
       data=(304-f.virs.ch4_all[*,min(id_virs):max(id_virs)])/2
       levels=304-indgen(61)*2
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='CH4 T!dB!n (K)'
       format='(i3)'
      endif
     endif
    if(varlists[k] eq 'NEARSURFZ') then begin
       id_pr=where(pr_lon[0,*] ge plotarea[1]-5 and pr_lon[0,*] le plotarea[3]+5)
       lon=pr_lon[*,min(id_pr):max(id_pr)]
       lat=f.pr.lat[*,min(id_pr):max(id_pr)]
       data=f.pr.nearsurfz[*,min(id_pr):max(id_pr)]
       levels=indgen(61)
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='PR near surface reflectivity (dBZ)'
       format='(i3)'
     endif
    if(varlists[k] eq 'NEARSURFRAIN') then begin
       id_pr=where(pr_lon[0,*] ge plotarea[1]-5 and pr_lon[0,*] le plotarea[3]+5)
       lon=pr_lon[*,min(id_pr):max(id_pr)]
       lat=f.pr.lat[*,min(id_pr):max(id_pr)]
       data=f.pr.nearsurfrain[*,min(id_pr):max(id_pr)]
       dtmp=data
       cid=where(data gt 0,ccc)
       if(ccc gt 0) then data[where(data gt 0)]=dtmp[where(data gt 0)]/60.0*45.0+15
       levels=indgen(61)
;       colors=indgen(60)
       colors=[0,indgen(59)/60.0*45+15]
       levelind=indgen(10)*6
       tic='2A25 near surface rainrate (mm/hr)'
       format='(i3)'
     endif
    if(varlists[k] eq 'STORMHT') then begin
       id_pr=where(pr_lon[0,*] ge plotarea[1]-5 and pr_lon[0,*] le plotarea[3]+5)
       lon=pr_lon[*,min(id_pr):max(id_pr)]
       lat=f.pr.lat[*,min(id_pr):max(id_pr)]
       data=(f.pr_2a23.stormh[*,min(id_pr):max(id_pr)]/1000.0*3) >0
       levels=indgen(61)/3.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='2A23 storm height (km)'
       format='(f5.1)'
     endif

    if(varlists[k] eq '2A12RAIN') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       data=f.tmi.rain[*,min(id_tmi):max(id_tmi)]
       dtmp=data
       cid=where(data gt 0,ccc)
       if(ccc gt 0) then data[where(data gt 0)]=dtmp[where(data gt 0)]/60.0*45.0+15
       levels=indgen(61)
;       colors=indgen(60)
       colors=[0,indgen(59)/60.0*45+15]
       levelind=indgen(10)*6
       tic='2A12 rainrate (mm/hr)'
       format='(i3)'
     endif
    if(varlists[k] eq 'STRATCONV') then begin
       id_pr=where(pr_lon[0,*] ge plotarea[1]-5 and pr_lon[0,*] le plotarea[3]+5)
       lon=pr_lon[*,min(id_pr):max(id_pr)]
       lat=f.pr.lat[*,min(id_pr):max(id_pr)]
       dat=f.pr_2a23.raintype2a23[*,min(id_pr):max(id_pr)]
       data=dat
       data[*,*]=0
       id=where(dat ge 100 and dat le 153)
       if(id[0] ge 0) then data[id]=30
       id=where(dat ge 200 and dat le 293)
       if(id[0] ge 0) then data[id]=45
       levels=[0,29,43,51]
       colors=[0,30,45]
       levelind=[0]
       tic='2A23 rain types '
       format='(i3)'
     endif
    if(varlists[k] eq '2A12PRECIPICE') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       dat=f.tmi.rain
       dat[*,*]=0
       dat[f.tmi.profile_ray,f.tmi.profile_scan]=total(f.tmi.precip_ice,2)/1000.0
       data=dat[*,min(id_tmi):max(id_tmi)]*4
       levels=indgen(61)/4.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='2A12 total colume !c precipitation ice (g/m!u2!n)'
       format='(f4.1)'
     endif
    if(varlists[k] eq '2A12PRECIPWATER') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       dat=f.tmi.rain
       dat[*,*]=0
       dat[f.tmi.profile_ray,f.tmi.profile_scan]=total(f.tmi.precip_water,2)/1000.0
       data=dat[*,min(id_tmi):max(id_tmi)]*4
       levels=indgen(61)/4.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='2A12 total colume !c precipitation water (g/m!u2!n)'
       format='(f4.1)'
     endif
    if(varlists[k] eq '2A12CLDWATER') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       dat=f.tmi.rain
       dat[*,*]=0
       dat[f.tmi.profile_ray,f.tmi.profile_scan]=total(f.tmi.cld_water,2)/1000.0
       data=dat[*,min(id_tmi):max(id_tmi)]*15
       levels=indgen(61)/15.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='2A12 total colume !c cloud water (g/m!u2!n)'
       format='(f4.1)'
     endif
    if(varlists[k] eq '2A12CLDICE') then begin
       id_tmi=where(tmi_lon[0,*] ge plotarea[1]-5 and tmi_lon[0,*] le plotarea[3]+5)
       lon=tmi_lon[*,min(id_tmi):max(id_tmi)]
       lat=f.tmi.lathi[*,min(id_tmi):max(id_tmi)]
       dat=f.tmi.rain
       dat[*,*]=0
       dat[f.tmi.profile_ray,f.tmi.profile_scan]=total(f.tmi.cld_ice,2)/1000.0
       data=dat[*,min(id_tmi):max(id_tmi)]*15
       levels=indgen(61)/15.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='2A12 total colume !c cloud ice (g/m!u2!n)'
       format='(f4.1)'
     endif
;
;    start to plot
;
;  print,plotarea
    
    !p.multi=[kw[k],ix,iy]
    if(k eq 0) then title=orbit+' '+date
    if(k gt 0) then title=''
     if(max(plotarea[[1,3]]) lt 180) then map_set,limit=plotarea,/hires,/noerase,xmargin=[2,xmargin],/noborder,ymargin=[4,4],$
          title=title,charsize=1.0,/cont
     if(max(plotarea[[1,3]]) ge 180) then map_set,0,180,limit=plotarea,/hires,/noerase,xmargin=[2,xmargin],/noborder,ymargin=[4,4],$
          title=title,charsize=1.0,/cont
     map_grid,/label,charsize=0.5
;      plot_pixels,data,lon,lat
      plot_orbit_pixels_color,data,lon,lat,levels=indgen(61),color=indgen(60)
    !p.multi=[kw[k],ix,iy]
     if(max(plotarea[[1,3]]) lt 180) then map_set,limit=plotarea,/hires,/noerase,xmargin=[2,xmargin],/noborder,ymargin=[4,4],$
         title=title,charsize=1.0,/cont
     if(max(plotarea[[1,3]]) ge 180) then map_set,0,180,limit=plotarea,/hires,/noerase,xmargin=[2,xmargin],/noborder,ymargin=[4,4],$
         title=title,charsize=1.0,/cont
     map_grid,/label,charsize=0.5
     xyouts,plotarea[1]+(plotarea[3]-plotarea[1])*0.05,$
            plotarea[2]-(plotarea[2]-plotarea[0])*0.05,tic,charsize=1.0
     if(n_elements(cross) ge 4) then begin                  ;draw cross section line
       for kc=0,n_elements(cross[0,*])-1 do begin
        oplot,cross[[0,2],kc],cross[[1,3],kc],psym=1,thick=2
        oplot,cross[[0,2],kc],cross[[1,3],kc],thick=3
       endfor
     endif
     if(n_elements(flashes) gt 0) then begin                ;draw flashes 
      id_flash=where(lis_lon ge plotarea[1]-3 and lis_lon le plotarea[3]+3)   
      if(id_flash[0] ge 0) then $
        oplot,lis_lon[id_flash],f.lis.flhlat[id_flash],psym=1,thick=4,symsize=0.5
     endif
      id_pr=where(pr_lon[0,*] ge plotarea[1]-3 and pr_lon[0,*] le plotarea[3]+3)
      oplot,pr_lon[0,id_pr],f.pr.lat[0,id_pr],line=2
      oplot,pr_lon[48,id_pr],f.pr.lat[48,id_pr],line=2
     if(varlists[k] ne 'STRATCONV') then colorbar,levels,colors,format=format,/ivestyle,/col,$
          xsize=0.005,/noframe,levelind=levelind,charsize=0.8
     if(varlists[k] eq 'STRATCONV') then colorbar,levels,colors,format=format,/ivestyle,/col,$
          xsize=0.005,levelind=[1,2,3],charsize=0.8,tics=['Norain','N/A','Strat','Conv']
   endfor 
 device,/close
 if(n_elements(noshow) eq 0) then spawn,'ggv '+outfile+'_panel.ps &'
end