; Calculates averaging kernels from a MOPITT retrieved covariance matrix
;   on 7-levels, using apriori profile and covariance matrix
;   for V3
;
; INPUT: 
;   covmat_ret: Retrieved MOPITT covariance matrix (7x7)
;   psurf:  Surface pressure, from MOPITT retrieval
;
; OUTPUT:
;   avgker: Averaging kernels (nlev x nlev)  
;   nlev: Number of levels in averaging kernels
;
; USES:
;   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,  13 June 2002
;   
;-------------------------------------------------------------------------

pro calc_avgker_v3, covmat_ret, psurf, avgker, nlev, retlev

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

; check for surface pressure being above lowest retrieval level (Psurf<850hPa)
; 'retlev' contains the pressures of the surface and the retrieval levels
;  with lower pressures.

if (psurf lt 850.) then begin
  levind = where(mopittlev lt psurf,nlev)
  retlev = [psurf, mopittlev[levind]]
  nlev = nlev+1
endif else begin
  retlev = [psurf, mopittlev[1:6]]
  nlev = 7
endelse


; Read the a priori data file (35 levels)

read_apriori_dat_v3, co_ap, ch4_ap, pres_ap, covmat_ap


; Make an NLEV (nominally 7) level a priori CO profile and cov. matrix, 
; that runs in same direction as MOPITT levels (1000 to 150 mb).
; Pick out the levels that correspond to MOPITT levels

ind7 = [31, 25, 22, 18, 15, 13, 11] 

retind = ind7[(7-nlev):6]
psind = where(pres_ap ge psurf)
retind[0] = psind[0]

co_ap_7 = co_ap[retind]
ctmp = covmat_ap[*,retind]
covmat_ap_7 = ctmp[retind,*]


; Invert the a priori covariance matrix (nlev)

Ca_1_7=invert(covmat_ap_7, status)

; status = 0: successful
; status = 1: singular array
; status = 2: warning that a small pivot element was used and that
;             significant accuracy was probably lost.
; Status is usually 2
;print,'Inversion status: ',status
if status eq 1 then begin
  print,'*** Using Identity Matrix instead of a priori cov matrix ***'
  Ca_1_7 = Identity(nlev)
endif

; Remove bad levels from retrieved covariance matrix (if psurf lt 850)

if (nlev lt 7) then begin
  c_ret = fltarr(nlev,nlev)
  for i = 7-nlev, 6 do begin
     j = i - (7-nlev)
     c_ret[j,*] = [reform(covmat_ret[i,(7-nlev):6])]
  endfor
  c_ret[0,0] = covmat_ret[0,0]
endif else $
  c_ret = covmat_ret

; Calculate averaging kernel matrix using matrix multiplication
; A = I - Cx Ca^-1

avgker = Identity(nlev) - c_ret ## Ca_1_7


end