;read_output.pro
; reads chem1d output- this reads the slightly smaller output file
;	which doesn't have all the production and loss terms
; modified 5/05 to default to 66 level version and read ozone at 54 km
;	reads the apr16.OUT file made by ajk on 4/07
; modified 11/18 to save key fields to look at diurnal clo and related species
; reads the hawaii files which are unformatted- unlike most other "OUT" files
; Nov 2023, used to read the supernova NOx runs
;
nalt=66
;	read,'enter nalt (e.g. 260)',nalt

nmol=89
nreact=254
z=fltarr(nalt)
conc=fltarr(nmol,nalt)
mix=conc & rprod = conc & rfreq = conc
rconst=fltarr(nreact,nalt)
rate=fltarr(nreact,nalt)
vec1=fltarr(nalt) & vec2=vec1 & vec3 = vec2 & vec4 = vec2
vec5=fltarr(nmol) & vec6 = fltarr(nmol)
outfile=''
read,'enter outfile prefix (e.g. betterclo, smiles)',outfile
ntime=975
ntime=892
ntime=999
ntime=550
if outfile eq 'betterhcl' or outfile eq 'acehclch4' or $
   outfile eq 'scaleupch4' or outfile eq 'ace_hih2o' then ntime = 800
time=fltarr(ntime)	;just time of day
sza_all=fltarr(ntime)
xtime = time		;net time of simulation
; outfile='forcedclo'
;  ntime=546		;stimpfle file is a tiny bit shorter for some reason
outfile = outfile + '.OUT'
close,1
openr,1,outfile,/f77_unformatted
;
; loop over time steps
;
ox56=fltarr(ntime) & o3 = ox56
o3 = fltarr(ntime,66) & den =  o3
ox = o3 & h = o3 & oh = o3 & o2 = o3
noden = o3 & no2den = o3
o3mix = o3
clomix = o3
clo_den = o3
o384 = o3 & o84 = o3 & h84 = o3 & oh84 = o3 & o84 = o3
o378 = o3 & o78 = o3 & h78 = o3 & oh78 = o3 & o78 = o3
o372 = o3 & o72 = o3 & h72 = o3 & oh72 = o3 & o72 = o3
no2mix = o3 & nomix = o3 & h2omix = o3 & h2mix = o3
hoclmix = o3 & ohmix = o3 & n2o5mix = o3 & hno3mix = o3
clono2mix = o3 & hclmix = o3 & ch4mix = o3 & ho2mix = o3
clmix = o3
ndays = 0.
for itime=0,ntime-1 do begin
; for itime = 0,60 do begin
readu,1,z
readu,1,thour,sza
;print,itime,thour,sza
if (itime gt 0 and thour lt time(itime-1)) then begin
  ndays = ndays + 1
  print,'new day',ndays,itime
  endif
time(itime)=thour
sza_all(itime)=sza
xtime(itime) = time(itime) + ndays*24.
;
readu,1,srtop,sstop,srbot,ssbot
for imol=0,nmol-1 do begin
 readu,1,vec1
 readu,1,vec2
 readu,1,vec3
 readu,1,vec4
 conc(imol,0)=transpose(vec1)
 mix(imol,0)=transpose(vec2)
 rprod(imol,0)=transpose(vec3)
 rfreq(imol,0)=transpose(vec4)
 endfor
;
for ireac=0,nreact-1 do begin
; readu,1,vec5		;for iform and idest, see chemout.f
; readu,1,vec6
 readu,1,vec1		;altitude length
 readu,1,vec2
 rconst(ireac,0)=transpose(vec1)
 rate(ireac,0)=transpose(vec2)
 endfor

ox56(itime)=conc(0,28) + conc(3,28)
ox(itime,0) = conc(0,*)
o3(itime,0)=conc(3,*)
noden(itime,0) = conc(7,*)
no2den(itime,0) = conc(8,*)
o3mix(itime,0) = mix(3,*)
nomix(itime,0) = mix(7,*)
no2mix(itime,0) = mix(8,*)
h2omix(itime,0) = mix(16,*)
h2mix(itime,0) = mix(15,*)
ho2mix(itime,0) = mix(18,*)
n2o5mix(itime,0) = mix(10,*)
hno3mix(itime,0) = mix(12,*)
clono2mix(itime,0) = mix(39,*)
hoclmix(itime,0) = mix(42,*)
clmix(itime,0) = mix(31,*)
oh(itime,0) = conc(17,*)
ohmix(itime,0) = mix(17,*)
clo_den(itime,0) = conc(32,*)
clomix(itime,0) = mix(32,*)
hclmix(itime,0) = mix(43,*)
ch4mix(itime,0) = mix(20,*)
den(itime,0) = conc(88,*)
o2(itime,0) =  conc(2,*)
;	mix56(itime)=mix(0,28) + mix(3,28)
o384(itime)=conc(3,42)
o84(itime) = conc(0,42)
h84(itime)=conc(14,42)
oh84(itime) = conc(17,42)
;
;o378(itime)=conc(3,39)
;o78(itime) = conc(0,39)
;h78(itime)=conc(14,39)
;oh78(itime) = conc(17,39)
;
;print,'itime,thour =',itime,thour
endfor
;time = time(214:271)		;This is specific to Apr16.OUT
z=findgen(nalt)*2.
o3den=o3
stop,'hit .con to save fields'
fileout = ''
read,'enter rootname of output save file',fileout
filesav = fileout + '.sav'
save,file=filesav,nomix,noden,no2mix,no2den,z,xtime,o3mix,clomix,clo_den,hoclmix, $
  ho2mix,clono2mix,hclmix,ohmix,n2o5mix,hno3mix,h2omix,h2mix,clmix,ch4mix,o3den,$
  time,sza_all
end