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

   file = '/users/mhervig/pmcs/clim/haloe_vpmc_nh_seasonal-series_all-years_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
      
;- set up the plot

   set_plot,'ps'
   
   plotfile = '/users/mhervig/pmcs/clim/haloe_vpmc_clim_plot.ps'
   plotfile = 'haloe_vpmc_clim_plot.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

   k = where( h2o_gridz(*,ila,iz) gt 0)
   
   plot,doygrid(k),h2o_gridz(k,ila,iz),position=[0.2,0.7,0.6,0.9],title=lat_tit,$
     xrange=[50,280],xtitle='day of year',$
     yrange=[0,12],ytitle='H!d2!nO (ppmv)'
   
     oplot,doygrid(k),h2o_gridz(k,ila,iz)+h2o_sd_gridz(k,ila,iz),line=2
     oplot,doygrid(k),h2o_gridz(k,ila,iz)-h2o_sd_gridz(k,ila,iz),line=2
     
     xyouts,55,10, 'a) ' + string(z1(iz),format='(f4.1)') + ' km'
     
;- plot a profile

   k = where(h2o_gridz(idy,ila,*) gt 0)
   
   plot,h2o_gridz(idy,ila,k),z1(k),position=[0.2,0.44,0.6,0.64],/noerase,$
     xrange=[0,12],xtitle='H!d2!nO (ppmv)', $
     yrange=[50,90],ytitle='altitude (km)'
   
     oplot,h2o_gridz(idy,ila,k)+h2o_sd_gridz(idy,ila,k),z1(k),line=2
     oplot,h2o_gridz(idy,ila,k)-h2o_sd_gridz(idy,ila,k),z1(k),line=2
   
     xyouts,0.5,55,'b) day ' + string(doygrid(idy),format='(i3)')
     
;- done

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