;	read in O3 loss rates from qqjk from 4 rate-limiting step loss reactions and 
;	the production rate for 1 year of MERRA. qk reactions are 2,24,49,72. They are
;	loss of O from O3, HO2, NO2, and ClO.

;	Compare for new merra and old merra

qkspec=[2,24,49,72]
;
dir1='/misc/gmibig1/public/gmic/MERRA/2005/'
dir2='/misc/gmibig1/public/gmic/MERRA_noclouds/2005/'
expt1='gmic_MERRA_2005_'
expt2='gmic_MERRA_noclouds_2005_'

mlab=['jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec']
dtype='.qqjk.nc'

file1=strarr(12) & file2=file1
for m=0,11 do file1(m)=dir1+expt1+mlab(m)+dtype
for m=0,11 do file2(m)=dir2+expt2+mlab(m)+dtype

mo3loss=fltarr(91,72,12,4)
omo3loss=fltarr(91,72,12,4)
mprod=fltarr(91,72,12) & omprod=mprod

;	Convert data into ppb/month

for m=0,11 do begin $
  fname=file1(m) & $
  for s=0,3 do begin $
    sp=qkspec(s) & $
    read_qqjk,'qqk',fname,sp,144,91,72,data,labels & $
    mo3loss(*,*,m,s)=avg(data,0)*8.64e4*30.5*1e9 & $
  endfor & $
  read_qqjk,'qqj',fname,1,144,91,72,data,labels & $
  mprod(*,*,m)=avg(data,0)*8.64e4*30.5*1e9 & $
;
  fname=file2(m) & $
  for s=0,3 do begin $
    sp=qkspec(s) & $
    read_qqjk,'qqk',fname,sp,144,91,72,data,labels & $
    omo3loss(*,*,m,s)=avg(data,0)*8.64e4*30.5*1e9 & $
  endfor & $
  read_qqjk,'qqj',fname,1,144,91,72,data,labels & $
  omprod(*,*,m)=avg(data,0)*8.64e4*30.5*1e9 & $
endfor