pro plot_gpm_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,level=level,outdata=outdata

;    
;   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  
;
;


;  hdf_2_struct,'../out/1Z.GPM.PF.20130930-S163408-E180630.090427.V00H.HDF',fdat
;  hdf_2_struct,'1Z.GPM.PF.20141201-S112043-E125315.004307.V01A.HDF',fdat
;plotarea=0
;cross=0

;goto,jump3
;set_plot,'x'
;!p.multi=0
;
; read in the data 
;
   hdf_2_struct,'/disk1/v7/soft/terrain.hdf',terr
   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
      hdf_2_struct,filename,fdat
   endif 
f=fdat
ku_lon=fdat.ku.lon
gmi_lon=fdat.gmi.lon
ka_lon=fdat.ka.lon
;dku_lon=fdat.dpr.lon
if(n_elements(plotarea) eq 4) then begin
if (max(plotarea[[1,3]]) ge 180 ) then begin
 id=where(ku_lon lt 0)
 ku_lon[id]=360+ku_lon[id]
 id=where(ka_lon lt 0)
 ka_lon[id]=360+ka_lon[id]
; id=where(dku_lon lt 0)
; dku_lon[id]=360+dku_lon[id]
 id=where(gmi_lon lt 0)
 gmi_lon[id]=360+gmi_lon[id]
endif
endif
 tics=strsplit(fdat.filename,'.',/extract)
 orbit=tics[n_elements(tics)-3]
   f=fdat
   if(n_elements(varlists) le 0) then begin
     varlists=['PCT85','PCT37','TB1837','KUNSZ','KUNSPRECIP','GMIPRECIP'] 
   endif
   orbit=orbit
;
; define the plotting area
;
 answer=''
 if(n_elements(plotarea) lt 4) then begin
   set_plot,'x'
   !p.multi=0
   ez_color,1,colors,ncolor=60 
  plotarea=[-66.,-180.,66,180]
  zoomin:
   id_ku=where(ku_lon ge plotarea[1] and ku_lon le plotarea[3] and $
               f.ku.lat ge plotarea[0] and f.ku.lat le plotarea[2])
   jutc=julday(f.ku.month[id_ku/49],f.ku.day[id_ku/49],f.ku.year[id_ku/49],f.ku.hour[id_ku/49],$
      f.ku.minute[id_ku/49],f.ku.second[id_ku/49])
   tjutc=jutc[sort(jutc)]   
   caldat,tjutc,tmonth,tday,tyear,thour,gminute,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(gminute[0])))+$
        ' UTC '
   print,plotarea
   map_set,limit=plotarea,/hires,/label,ymargin=[6,2],title=date+orbit+' Ku near surface reflectivity '
;    plot_pixels,f.ku.nearsurfz,ku_lon,f.ku.lat
    contour,f.ku.nsz,ku_lon,f.ku.lat,levels=indgen(10)*6,$
            c_color=colors[indgen(10)*6],/overplot,c_label=intarr(10)+1,/fill
    contour,f.gmi.pct85,gmi_lon,f.gmi.lat,$
        levels=[250],/overplot      
    oplot,f.ka.lon[0,*],f.ka.lat[0,*],line=1
    oplot,f.ka.lon[23,*],f.ka.lat[23,*],line=1
   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,/data
        print,'x1, y1 =',x1,y1
        read,'Press return',answer
        print,'Click on upper right corner of desired zoom'
        cursor,x2,y2,/data
        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,/data
        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,/data
        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_ku=where(ku_lon ge plotarea[1] and ku_lon le plotarea[3] and $
               f.ku.lat ge plotarea[0] and f.ku.lat le plotarea[2])
   jutc=julday(f.ku.month[id_ku/49],f.ku.day[id_ku/49],f.ku.year[id_ku/49],f.ku.hour[id_ku/49],$
               f.ku.minute[id_ku/49],f.ku.second[id_ku/49])
   tjutc=jutc[sort(jutc)]   
   caldat,tjutc,tmonth,tday,tyear,thour,gminute,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(gminute[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'
       device,/color,filename=outcross+'.ps',xsize=25,ysize=25
       !p.multi=[0,1,2] 
       id_ku=where(ku_lon ge min(cross[[0,2],k])-0.10 and ku_lon le max(cross[[0,2],k])+0.10 and $
               f.ku.lat ge  min(cross[[1,3],k])-0.10 and f.ku.lat le  max(cross[[1,3],k])+0.10)
       dlatlon=cos((cross[1,k]+cross[3,k])/2.0*!pi/180.0)*111.0
       distance=sqrt((dlatlon*$
              (cross[2,k]-cross[0,k]))^2.0+((cross[3,k]-cross[1,k])*dlatlon)^2.0) 
        print,distance
        dis=abs((cross[2,k]-cross[0,k])*(cross[1,k]-f.ku.lat[id_ku])-$
                 (cross[3,k]-cross[1,k])*(cross[0,k]-ku_lon[id_ku])) $
           /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_ku[idc]
          dis1=sqrt(((cross[0,k]-ku_lon[id])*dlatlon)^2.0+((cross[1,k]-f.ku.lat[id])*111.)^2.0$
               -(vdis)^2.0)
          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_ku=where(ku_lon[0,*] ge plotarea[1]-3 and ku_lon[0,*] le plotarea[3]+3)
        lon=ku_lon[*,min(id_ku):max(id_ku)]
        lat=f.ku.lat[*,min(id_ku):max(id_ku)]
        data=f.ku.nearsurfz[*,min(id_ku):max(id_ku)]
        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,ku_lon[id],f.ku.lat[id],psym=2
         oplot,ku_lon[id],f.ku.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))
        tb1837=fltarr(n_elements(id))
        hilon=fltarr(n_elements(id))
        hilat=fltarr(n_elements(id))
        dbz_cross=fltarr(n_elements(id),176)
        dbz_cross=fltarr(n_elements(id),176)
        rain_ku=fltarr(n_elements(id))
        rain_gmi=fltarr(n_elements(id))
        xlon=fltarr(n_elements(id)) 
        xlat=fltarr(n_elements(id)) 
        elev=fltarr(n_elements(id)) 

        for i=0,n_elements(id)-1 do begin
          rain_ku[i]=f.ku.nsprecip[id[i]]/10.0
          xlon[i]=ku_lon[id[i]]
          xlat[i]=f.ku.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]
          rain_gmi[i]=f.gmi.precip[f.colo.pku[id[i]]]/10.0
          pct85[i]=f.gmi.pct85[f.colo.pku[id[i]]]
          pct37[i]=f.gmi.pct37[f.colo.pku[id[i]]]
          tb1837[i]=(reform(f.gmi.hitb[3,*,*]))[f.colohi.pku[id[i]]]
          hilon[i]=(f.gmi.hilon[*,*])[f.colohi.pku[id[i]]]
          hilat[i]=(f.gmi.hilat[*,*])[f.colohi.pku[id[i]]]
;          did=where(id[i] eq f.ray+f.scan*49)
           scan=long(id[i]/49)
           ray=id[i]-scan*49
;          print,xpicsize[i]
          if(xpicsize[i] gt 0 and xpicsize[i] lt 7.0) then begin
            dbz=reform(f.ku.dbz[*,ray,scan])
            kudbz=reform(f.ku.dbz[*,ray,scan])
            dbz_cross[i,*]=dbz
            z=(176-indgen(176))/8.0+0.125
            x=intarr(80)+dis1[i]-xpicsize[i]*0.5
            did=where(dbz ge 13 and z lt 20)
            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_ku,line=0,thick=2
        oplot,dis1-xpicsize*0.5,rain_gmi,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
        oplot,dis1-xpicsize*0.5,tb1837,thick=5,line=1
        colorbar,indgen(61),indgen(60),levelind=indgen(10)*6,format='(i2)',/ivestyle,/noframe,ysize=0.01

       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]
       matchedka=f.ku.dbz
       matchedka[*]=0
       matchedka[*,12:36,*]=f.dpr.dbz
       for i=0,n_elements(id)-1 do begin
           scan=long(id[i]/49)
           ray=id[i]-scan*49
          if(xpicsize[i] gt 0 and xpicsize[i] lt 7.0) then begin
            dbz=reform(matchedka[*,ray,scan])
            dbz_cross[i,*]=dbz
            z=(176-indgen(176))/8.0+0.125
            x=intarr(80)+dis1[i]-xpicsize[i]*0.5
            did=where(dbz ge 12 and z lt 20)
            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
        colorbar,indgen(61),indgen(60),levelind=indgen(10)*6,format='(i2)',/ivestyle,/noframe,ysize=0.01
       device,/close
       spawn,'convert -density 100 '+outcross+'.ps '+outcross+'.png'
     endfor
  endif
; stop
;
; panels
;

 set_plot,'ps'
 ez_color,1,colors,ncolor=60
 if(n_elements(outfile) eq 0) then outfile=orbit
 device,/color,filename=outfile+'.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+'.ps',xsize=xsize,ysize=ysize
 !p.font=6
 
   for k=0,n_elements(varlists)-1 do begin
     print,varlists[k]
     if(varlists[k] eq 'PCT85') then begin
       id_gmi=where(gmi_lon[0,*] ge plotarea[1]-15 and gmi_lon[0,*] le plotarea[3]+15)
       lon=gmi_lon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.lat[*,min(id_gmi):max(id_gmi)]
;       data=(300-f.gmi.pct85[*,min(id_gmi):max(id_gmi)])/2.5
       data=f.gmi.pct85[*,min(id_gmi):max(id_gmi)]
       levels=reverse(300-indgen(61)*2.5)
       colors=reverse(indgen(60))
       levelind=indgen(10)*6
       tic=' 89GHz PCT (K)'
       format='(i3)'
     endif     
     if(varlists[k] eq 'TB1837') then begin
       id_gmi=where(f.gmi.hilon[0,*] ge plotarea[1]-15 and f.gmi.hilon[0,*] le plotarea[3]+15)
       lon=f.gmi.hilon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.hilat[*,min(id_gmi):max(id_gmi)]
;       data=(290-reform(f.gmi.hitb[3,*,min(id_gmi):max(id_gmi)]))/3.0
        data=reform(f.gmi.hitb[3,*,min(id_gmi):max(id_gmi)])
       levels=reverse(290-indgen(61)*3.0)
       colors=reverse(indgen(60))
       levelind=indgen(10)*6
       tic=' TB 183+7 (K)'
       format='(i3)'
     endif

    if(varlists[k] eq 'PCT37') then begin
       id_gmi=where(f.gmi.lon[0,*] ge plotarea[1]-15 and f.gmi.lon[0,*] le plotarea[3]+15)
       lon=f.gmi.lon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.lat[*,min(id_gmi):max(id_gmi)]
;       tmpdata=reform((0.545*f.gmi.tb[6,*,*]-f.gmi.tb[5,*,*])/(0.545-1))
;       data=(315-tmpdata[*,min(id_gmi):max(id_gmi)])/1.5
       data=f.gmi.pct37[*,min(id_gmi):max(id_gmi)]
;       data=tmpdata[*,min(id_gmi):max(id_gmi)]
       levels=reverse(315-indgen(61)*1.5)
       colors=reverse(indgen(60))
       levelind=indgen(10)*6
       tic=' 37GHz PCT (K)'
       format='(i3)'
     endif

    if(varlists[k] eq 'V19') then begin
       id_gmi=where(f.gmi.lon[0,*] ge plotarea[1]-5 and f.gmi.lon[0,*] le plotarea[3]+5)
       lon=f.gmi.lon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.lat[*,min(id_gmi):max(id_gmi)]
;       data=61-(324-reform(f.gmi.tb[2,*,min(id_gmi):max(id_gmi)]))/2.2
       data=reform(f.gmi.tb[2,*,min(id_gmi):max(id_gmi)])
       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_gmi=where(f.gmi.lon[0,*] ge plotarea[1]-5 and f.gmi.lon[0,*] le plotarea[3]+5)
       lon=f.gmi.lon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.lat[*,min(id_gmi):max(id_gmi)]
;       data=61-(324-reform(f.gmi.tb[3,*,min(id_gmi):max(id_gmi)]))/2.2
       data=reform(f.gmi.tb[3,*,min(id_gmi):max(id_gmi)])
       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 'KUNSZ') then begin
       id_ku=where(ku_lon[0,*] ge plotarea[1]-5 and ku_lon[0,*] le plotarea[3]+5)
       lon=ku_lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       data=f.ku.nsz[*,min(id_ku):max(id_ku)]
       levels=indgen(61)
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='Ku near surface reflectivity (dBZ)'
       format='(i3)'
     endif

   if(varlists[k] eq 'DBZCROSS') then begin
       id_ku=where(ku_lon[0,*] ge plotarea[1]-5 and ku_lon[0,*] le plotarea[3]+5)
       lon=ku_lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       if(n_elements(level) eq 0) then level=1
       z=22-indgen(176)/8.0
       zid=where(abs(z-level) eq min(abs(z-level)))
       data=reform(f.ku.dbz[zid[0],*,min(id_ku):max(id_ku)])
       levels=indgen(61)
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='Ku reflectivity at '+string(z[zid[0]],format='(f6.2)')+' km (dBZ)'
       format='(i3)'
     endif

    if(varlists[k] eq 'KUNSPRECIP') then begin
       id_ku=where(ku_lon[0,*] ge plotarea[1]-5 and ku_lon[0,*] le plotarea[3]+5)
       lon=ku_lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       data=f.ku.nsprecip[*,min(id_ku):max(id_ku)]
       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=10^(indgen(61)/60.0*3-1)
       colors=indgen(60)
;       colors=[0,indgen(59)/60.0*45+15]
       levelind=indgen(10)*6
       tic='Ku near surface preciprate (mm/hr)'
       format='(f5.1)'
     endif

    if(varlists[k] eq 'KANSZ') then begin
       id_ka=where(ka_lon[0,*] ge plotarea[1]-5 and ka_lon[0,*] le plotarea[3]+5)
       lon=ka_lon[*,min(id_ka):max(id_ka)]
       lat=f.ka.lat[*,min(id_ka):max(id_ka)]
       data=f.ka.nsz[*,min(id_ka):max(id_ka)]
       levels=indgen(61)
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='Ka near surface reflectivity (dBZ)'
       format='(i3)'
     endif
    if(varlists[k] eq 'KANSPRECIP') then begin
       id_ka=where(ka_lon[0,*] ge plotarea[1]-5 and ka_lon[0,*] le plotarea[3]+5)
       lon=ka_lon[*,min(id_ka):max(id_ka)]
       lat=f.ka.lat[*,min(id_ka):max(id_ka)]
       data=f.ka.nsprecip[*,min(id_ka):max(id_ka)]
       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)
       levels=10^(indgen(61)/60.0*3-1)
       colors=indgen(60)
;       colors=[0,indgen(59)/60.0*45+15]
       levelind=indgen(10)*6
       tic='Ka near surface preciprate (mm/hr)'
       format='(f5.1)'
     endif


    if(varlists[k] eq 'STORMHT') then begin
       id_ku=where(ku_lon[0,*] ge plotarea[1]-5 and ku_lon[0,*] le plotarea[3]+5)
       lon=ku_lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       data=(f.ku.stormht[*,min(id_ku):max(id_ku)]/1000.0) >0
       levels=indgen(61)/3.0
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='Ku storm height (km)'
       format='(f5.1)'
     endif

    if(varlists[k] eq 'VIL') then begin
       id_ku=where(f.ku.lon[0,*] ge plotarea[1]-5 and f.ku.lon[0,*] le plotarea[3]+5)
       lon=f.ku.lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       ku_z=10^((f.ku.dbz<56)/10.0)
;       ku_z=f.ku.dbz<56
       data=fltarr(49,n_elements(id_ku))
       z=(176-indgen(176))/8.0+0.125
       for i=0L,n_elements(id_ku)-1 do begin
        for j=0,48 do begin
          ntop=(22.125-f.ku.stormht[j,id_ku[i]]/1000.)*8
          nbot=(22.125-terr.terr[(f.ku.lon[j,id_ku[i]]+180.0)/0.16666666,(f.ku.lat[j,id_ku[i]]+90.0)/0.16666666]/1000.0-0.5)*8
;          if(ntop lt nbot) then print,z[ntop],z[nbot],f.ku.stormht[j,id_ku[i]]/1000.,terr.terr[(f.ku.lon[j,id_ku[i]]+180.0)/0.16666666,(f.ku.lat[j,id_ku[i]]+90.0)/0.16666666]/1000.0
          for n=ntop,nbot-1 do begin
           data[j,i]=data[j,i]+((ku_z[n,j,id_ku[i]]+ku_z[n+1,j,id_ku[i]])/2.0)^0.571429*125*3.44/1000000.
          endfor
        endfor
       endfor
;       print,max(data)
;       stop
       levels=10^(indgen(61)/60.0*2.5-1)
       colors=indgen(60)
       levelind=indgen(10)*6
       tic='VIL (kg/m!u2!n)'
       format='(f4.1)'
     endif

    if(varlists[k] eq 'GMIPRECIP') then begin
       id_gmi=where(gmi_lon[0,*] ge plotarea[1]-15 and gmi_lon[0,*] le plotarea[3]+15)
       lon=gmi_lon[*,min(id_gmi):max(id_gmi)]
       lat=f.gmi.lat[*,min(id_gmi):max(id_gmi)]
       data=f.gmi.precip[*,min(id_gmi):max(id_gmi)]
       cid=where(data lt 0.1,ccc)
       if(ccc gt 0) then data[cid]=0
       dtmp=data
       cid=where(data gt 0,ccc)
;       if(ccc gt 0) then data[cid]=dtmp[cid]/60.0*45.0+15
       levels=10^(indgen(61)/60.0*3-1)
       colors=indgen(60)
;       colors=[0,indgen(59)/60.0*45+15]
       levelind=indgen(10)*6
       tic='GMI preciprate (mm/hr)'
       format='(f5.1)'
     endif

    if(varlists[k] eq 'STRATCONV') then begin
       id_ku=where(ku_lon[0,*] ge plotarea[1]-5 and ku_lon[0,*] le plotarea[3]+5)
       lon=ku_lon[*,min(id_ku):max(id_ku)]
       lat=f.ku.lat[*,min(id_ku):max(id_ku)]
       dat=f.ku.preciptype[*,min(id_ku):max(id_ku)]
       data=dat
       data[*,*]=0
       id=where(dat eq 1)
       if(id[0] ge 0) then data[id]=30
       id=where(dat eq 2)
       if(id[0] ge 0) then data[id]=45
       levels=[0,29,43,51]
       colors=[0,30,45]
       levelind=[0]
       tic='precipitation types '
       format='(i3)'
     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)
      plot_orbit_pixels_color,data,lon,lat,levels=levels,color=colors
      outdata={lon:lon,lat:lat,data:data}
    !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
      id_ku=where(ku_lon[0,*] ge plotarea[1]-3 and ku_lon[0,*] le plotarea[3]+3)
      oplot,ku_lon[0,id_ku],f.ku.lat[0,id_ku],line=2
      oplot,ku_lon[48,id_ku],f.ku.lat[48,id_ku],line=2
      oplot,ka_lon[0,id_ku],f.ka.lat[0,id_ku],line=1
      oplot,ka_lon[23,id_ku],f.ka.lat[23,id_ku],line=1
;
;      oplot,hilon,hilat,line=1,thick=3
     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
 spawn,'convert -density 100 '+outfile+'.ps '+outfile+'.png'
end