; Example program to illustrate reading covariance matrix and calculating 
;   averaging kernels.
;
; OUTPUT:
;   plot of averaging kernels
;
; SUBROUTINES CALLED:
;   calc_avger_v3.pro
;   read_apriori_dat_v3: reads the MOPITT apriori data file
;
; NOTE: apriori arrays and matrix run from low pressure to surface,
;   while MOPITT output runs from surface upwards
; 
; L. Emmons, NCAR,  9 July 2002
;   
;-------------------------------------------------------------------------

@calc_avgker_v3

pro example_avgker_calc

mopittlev=[1000., 850., 700., 500., 350., 250., 150.]

;-------------------------------------------------------------------------
; Read a priori profile (convert to ppbv)
;-------------------------------------------------------------------------
print,'Reading V3 apriori profile'
read_apriori_dat_v3,apriori_co,apriori_ch4,apriori_pres,covmat_ap
apriori_co = apriori_co * 1.e9  ; convert to ppbv

;-------------------------------------------------------------------------
; Example aircraft data (already binned to standard MOPITT levels
;-------------------------------------------------------------------------
ac_pres = [1000., 850., 700., 500., 350., 250., 150.]
ac_co = [ 75, 75, 75, 80, 85, 75, 60]

;-------------------------------------------------------------------------
; Example MOPITT data: 
;   retrieved profile, errors, covariance matrix, surface pres.
;   commands for reading a MOPITT HDF file are:
;    co_surf_mixing_ratio = get_sd(mopfile, 'Retrieval Bottom CO Mixing Ratio')
;    co_mixing_ratio = get_sd(mopfile, 'CO Mixing Ratio')
;    covmatrix = get_sd(mopfile, 'Retrieval Error Covariance Matrix')
;    psurf = get_sd(mopfile, 'Retrieval Bottom Pressure')
;   which will return all pixels, not the single pixel illustrated here.
;   See the data file specification for details.
;-------------------------------------------------------------------------

; from file, MOP02-20001101-L2V5.3.1.hdf, pixel 141707, at (26.10S, 142.45W)

 psurf = 1019.10
 tsurf = 293.815
 co_surf_mixing_ratio = [ 79.9311,   70.4192]
 co_mixing_ratio = fltarr(2,6)
 co_mixing_ratio[0,*] = [71.8540, 72.5114, 80.7527, 83.1405, 76.5784, 60.0378]
 co_mixing_ratio[1,*] = [30.3219, 21.2398, 16.3888, 12.4246, 11.4563, 14.0957]
 covmatrix = [6.52006e-16, -3.08134e-16, -3.68012e-16, -2.05053e-16,  $
  1.79254e-17,  6.88825e-17, 1.72149e-17, -2.33923e-16, -1.94607e-16, $
  -3.66491e-17 , 6.19993e-17, -1.33302e-16, -1.09088e-16, -6.29910e-17, $
  -4.98229e-17,  6.78655e-17, -6.35391e-17, -1.00982e-16, 2.13220e-17, $
  -2.95665e-17,  9.61127e-17]

help,co_mixing_ratio, covmatrix

  co_surf_err = co_surf_mixing_ratio[1]
  co_err = co_mixing_ratio[1,*]
  mopittlev[0] = psurf
  
;-------------------------------------------------------------------------
; Construct full retrieved covariance matrix (7x7 array)
;-------------------------------------------------------------------------
  cov_retr = dblarr(7,7)
  cov_retr[0,0] = (co_surf_err*1.e-9)^2
  for i=1,6 do cov_retr[i,i] = (co_err[i-1]*1.e-9)^2
  cov_retr[0,1:6] = covmatrix[0:5]
  cov_retr[1,2:6] = covmatrix[6:10]
  cov_retr[2,3:6] = covmatrix[11:14]
  cov_retr[3,4:6] = covmatrix[15:17]
  cov_retr[4,5:6] = covmatrix[18:19]
  cov_retr[5,6] = covmatrix[20]
  for j=0,5 do $
    for i=1,6 do cov_retr[i,j]=cov_retr[j,i]

  ;-------------------------------------------------------------------------
  ; Calculate the averaging kernels
  ;   retlev contains the nlev levels that avgker is for (psurf and 
  ;   the levels above)
  ;-------------------------------------------------------------------------
  calc_avgker_v3, cov_retr, psurf, avgker, nlev, retlev

  ;-------------------------------------------------------------------------
  ; Interpolate aircraft ("true") profile to MOPITT retrieval levels
  ;-------------------------------------------------------------------------
  xac_7 = Interpol(ac_co, ac_pres, retlev)

  ;-------------------------------------------------------------------------
  ; Interpolate apriori profile to MOPITT retrieval levels
  ;-------------------------------------------------------------------------
  xa = Interpol(apriori_co, apriori_pres, retlev)

  ;-------------------------------------------------------------------------
  ; Calculate the "retrieved" in situ value
  ;-------------------------------------------------------------------------
  xr = xa + avgker ## (xac_7 - xa)


  ; Plot profiles and averaging kernels
  set_plot,'ps'
  device,/color,/landscape,file = 'example_avgker_calc.ps'
  !p.font=0

  !p.multi=[0,2,1]
  !p.charsize=1.2
  loadct, 12
  xmax = Max([xr,ac_co,xa])

  plot,ac_co,ac_pres, psym=-4,thick=3,/ylog,yrange=[1060,100],/ystyle, $
    xrange=[0,xmax],xtitle='CO (ppbv)', ytitle='Pres (hPa)' 
  oplot,xr,retlev, thick=3,linest=2,psym=-2, color=180
  oplot,apriori_co, apriori_pres,linest=1, thick=3, color=70

oplot, [2, 20], [110, 110], thick=3,linest=0,psym=-4
xyouts, 20, 110, ' True CO', size=1.2
oplot, [2, 20], [120, 120], thick=3,linest=2,psym=-2,color=180
xyouts, 20, 120, ' Retrieved CO', size=1.2
oplot, [2, 20], [133, 133], thick=3, linestyle=1, color=70
xyouts, 20, 133, ' Apriori CO', size=1.2

plot,reform(avgker[*,0]),retlev,/ylog,yrange=[1060,100],/ystyle,thick=3, $
    xrange=[-0.2,0.7],/xstyle, xtitle='Averaging Kernels', $
    title = String(psurf,tsurf,format='("Psurf: ",f6.1," mb, Tsurf: ",f5.1," K")')
oplot, reform(avgker[*,1]), retlev, linestyle=0, psym=-1, color=70
oplot, reform(avgker[*,2]), retlev, linestyle=0, psym=-2, color=90
oplot, reform(avgker[*,3]), retlev, linestyle=4, psym=-3, color=100
oplot, reform(avgker[*,4]), retlev, linestyle=6, psym=-4, color=120
oplot, reform(avgker[*,5]), retlev, linestyle=6, psym=-5, color=150
oplot, reform(avgker[*,6]), retlev, linestyle=6 , psym=-6, color=200
 
xbar1 = !x.window[0] + 0.05*(!x.window[1]-!x.window[0])
xbar2 = !x.window[0] + 0.15*(!x.window[1]-!x.window[0])
ybar = 0.95*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm, thick=3
xyouts,1.05*xbar2,ybar,'Surface',/norm, charsize=1.
ybar = 0.90*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm, psym=-1, color=70
xyouts,1.05*xbar2,ybar,'850 hPa',/norm, charsize=1.
ybar = 0.85*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm, psym=-2, color=90
xyouts,1.05*xbar2,ybar,'700 hPa',/norm, charsize=1.
xbar1 = !x.window[0] + 0.65*(!x.window[1]-!x.window[0])
xbar2 = !x.window[0] + 0.77*(!x.window[1]-!x.window[0])
ybar = 0.95*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm,psym=-3,linestyle=4,color=100
xyouts,1.02*xbar2,ybar,'500 hPa',/norm, charsize=1.
ybar = 0.90*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm,psym=-4,linestyle=6,color=120
xyouts,1.02*xbar2,ybar,'350 hPa',/norm,  charsize=1.
ybar = 0.85*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm,psym=-5,linestyle=6, color=150
xyouts,1.02*xbar2,ybar,'250 hPa',/norm, charsize=1.
ybar = 0.8*(!y.window[1]-!y.window[0]) + !y.window[0]
plots, [xbar1,xbar2],[ybar,ybar],/norm,psym=-6,linestyle=6, color=200
xyouts,1.02*xbar2,ybar,'150 hPa',/norm, charsize=1.


device,/close

end