;-------------------------------------------------------------------------------------
;
; Routine shows examples of how to use the HALOE Vpmc climatology data files.
; This version is for the 14-year average seasonal time series data.
;
; Source: Mark Hervig
;
;-------------------------------------------------------------------------------------


;- get the climatology data file

   ans = ' '
   read,'Read the climatology file (y/n)? ',ans
   if (ans eq 'y') then begin
   
   file = '/users/mhervig/pmcs/clim/haloe_vpmc_clim_nh_time-series_all-obs.idlsav'   
   file = dialog_pickfile(title='Choose the HALOE PMC climatology file')
   
   read_haloe_vpmc_climatology,file,$
        yeargrid,doygrid,latgrid,z1,$
        nob_grid,fpmc_grid,$
       pmcz_grid,pmcz_sd_grid,$
     pmcemx_grid,pmcemx_sd_grid,$
       
         t_gridz,    t_sd_gridz,$
         p_gridz,    p_sd_gridz,$
       h2o_gridz,  h2o_sd_gridz,$
        o3_gridz,   o3_sd_gridz,$
        no_gridz,   no_sd_gridz,$
      a245_gridz, a245_sd_gridz,$
      a340_gridz, a340_sd_gridz,$
      a346_gridz, a346_sd_gridz,$
      a526_gridz, a526_sd_gridz,$
      a626_gridz, a626_sd_gridz
      
   endif
   
;- set up the plot

   set_plot,'ps'
   
   plotfile = '/users/mhervig/pmcs/clim/haloe_vpmc_clim_plot2.ps'
   plotfile = 'haloe_vpmc_clim_plot2.ps'
   
   device,xoffset=2.0,xsize=17,yoffset=2.0,ysize=22,filename=plotfile

   !y.style     = 1
   !x.style     = 1
   
;- get array indices for plotting 
   
   zin = 83.  ; Z, km
   lin = 67.  ; latitude, deg N
   din = 197. ; day of year
   
   d = abs(z1      - zin) & k=where(d eq min(d)) & iz  = k(0)  ; altitude index
   d = abs(latgrid - lin) & k=where(d eq min(d)) & ila = k(0)  ; latitude index   
   d = abs(doygrid - din) & k=where(d eq min(d)) & idy = k(0)
   
   lat_tit = 'HALOE Vpmc Climatology, ' + string(latgrid(ila),format='(f6.2)') + ' deg. latitude'
   
;- plot a seasonal time series   
   
   plot,doygrid,h2o_gridz(*,ila,iz),psym=4,position=[0.1,0.7,0.9,0.9],title=lat_tit ,$
     xrange=[50,280],xtitle='day of year',$
     yrange=[0,11],ytitle='H!d2!nO (ppmv)'
     
     xyouts,70,9, 'a) all years, ' + string(z1(iz),format='(f4.1)') + ' km'
     
;- plot a profile

   k = where(h2o_gridz(*,ila,iz) gt 0)
   
   plot,yeargrid(*),h2o_gridz(*,ila,iz), position=[0.1,0.44,0.9,0.64],/noerase,$
     xrange=[1992,2006],xtitle='year', $
     yrange=[0,11],ytitle='H!d2!nO (ppmv)'
   
     xyouts,1993,9,'b) ' + string(z1(iz),format='(f4.1)') + ' km'
     
;- done

   device,/close
   set_plot,'x'
   print,'done'
   end