pro get_tgcm_on2_day,histfile,slt,on2_tmp

; some constants

missing=-999.00

; read tie-gcm model output(secondary history file)
read_netcdf,histfile,d,attr

lon=d.lon             ; longitude array
lat=d.lat             ; latitude array
lev=d.lev             ; vertical model level array 
time=d.time

nlon=n_elements(lon)
nlat=n_elements(lat)
nlev=n_elements(lev)
ntime=n_elements(time)

; get mass mixing ratio of O1,O2,neutral temperature TN, and geopotential height Z
; all the four fields are 4 dimensional array, e.g., TN(nlon,nlat,nlev,ntime)
; note: O1,O2, and TN are on midpoint, geopotential height gz is on interface
; in tiegcm1.6 history file, only geopotential height zint is on the model interface, all 
; the other variables(i.e. O1,O2,etc) are on half model levels (midpoint), we
; we will do all our calculation on half model levels.

o1=d.o1
o2=d.o2
tn=d.tn
gzint=d.z

; get total number density, total mass density, mean molecular weight
; and geopotential height at model mid-point(half level) 

mmr2nd,o1,o2,tn,gzint,lat,lev,nlat,nlon,ntime,nt,rhot,mbar,gzmid

; calculate number density of O and N2
no1=fltarr(nlon,nlat,nlev,ntime)
nn2=fltarr(nlon,nlat,nlev,ntime)
for i=0,nlon-1 do begin
   for j=0,nlat-1 do begin
      for k=0,nlev-1 do begin
         no1[i,j,k,*]=o1[i,j,k,*]*nt[i,j,k,*]*mbar[i,j,k,*]/16.
         nn2[i,j,k,*]=(1-o1[i,j,k,*]-o2[i,j,k,*])*nt[i,j,k,*]*mbar[i,j,k,*]/28.
      endfor
   endfor
endfor

; TIEGCM history was saved every hour for each day starting UT=0
ut=intarr(ntime)
ut[*]=fix(findgen(ntime))

; define output 
on2_tmp=fltarr(nlat)

;define tmp array
on2_tgcm=fltarr(nlat,ntime)
ilon=intarr(nlat)
the_lon=fltarr(nlat)

;start time loop
for m=0,ntime-1 do begin

   ; initializing:set on2 to missing
   on2_tgcm[*,m]=missing

   ;find ilon
   ; handle missing SLT situation. slt(nlat)
   ; here we assume TIMED 15 orbits on a same day have same local time 
   ; for each latitude zone
   ;find ilon for each latitude
   index=where(slt eq missing,count)
   if count gt 0 then begin 
       ilon[index]=missing
       the_lon[index]=missing
   endif
   index=where(slt ne missing,count)
   if count gt 0 then the_lon[index]=15.*(slt[index]-ut[m])
   index=where((the_lon lt -180.) and (the_lon ne missing),count)
   if count gt 0 then the_lon[index]=the_lon[index]+360.
   index=where(the_lon gt 175.,count)
   if count gt 0 then the_lon[index]=the_lon[index]-360.
   ;
   for j=0,nlat -1 do begin
      if the_lon[j] ne missing then begin
         lon_tmp=abs(lon-the_lon[j])
         index=where(lon_tmp eq min(lon_tmp))
         ilon[j]=index[0]
      endif
   endfor

   ; integrated column number density at each latiude for all ilon
   ; referenced at where N2 integrated column number density is 1.e17cm^-2
   no1_1=fltarr(nlat)
   nn2_1=fltarr(nlat)           ; nn2_1[*]=1.e17 cm^-2
   for j=0,nlat-1 do begin      ; zbottom and ztop are on2 dimensional array with number of elements of nilon
      if ilon[j] ne missing then begin
      ; make vertical grid that has 1 km resolution
      zbottom=fix(gzmid[ilon[j],j,0,m])
      ztop=fix(gzmid[ilon[j],j,nlev-2,m])   ; top level nlev-1 has gzmid=0.
      nz=(ztop-zbottom)+1
      z=findgen(nz)+zbottom
      ; interpolate (log) O and N2 number density onto 1 km resolution
      result_o1=interpol(alog(no1[ilon[j],j,0:nlev-2,m]),gzmid[ilon[j],j,0:nlev-2,m],z)
      result_o1=exp(result_o1)
      result_n2=interpol(alog(nn2[ilon[j],j,0:nlev-2,m]),gzmid[ilon[j],j,0:nlev-2,m],z)
      result_n2=exp(result_n2)

      ; calculate column number density for each model level. Column number
      ; density is calculated from top of the atmosphere doward. The column
      ;depth of each model level is from the top of the model to that model level
      no1_2=fltarr(nz)
      nn2_2=fltarr(nz)
      no1_2[nz-1]=result_o1[nz-1]*1.e5
      nn2_2[nz-1]=result_n2[nz-1]*1.e5
      for k=nz-2,0,-1 do begin
         no1_2[k]=no1_2[k+1]+result_o1[k]*1.e5
         nn2_2[k]=nn2_2[k+1]+result_n2[k]*1.e5
      endfor
      ; find the altitude where N2 column number density is 1.e17
      index=where(nn2_2 le 1.e17)
      id1=index[0]-1
      id2=index[0]
   ;   print,'id2=',id2,'z=',z(id2)
   ;   print,'id1=',id1,'z=',z(id1)
      zx=z[id2]+ $
         (alog(1.e17)-alog(nn2_2[id2]))/  $
         (alog(nn2_2[id1])-alog(nn2_2[id2]))* $
         (z[id1]-z[id2])
      print,'zx=',zx
      ; obtain O column number density at this altitude
      no1_1[j]=alog(no1_2[id2])+(alog(no1_2[id1])-alog(no1_2[id2]))/ $
                   (z[id1]-z[id2])*(zx-z[id2])
      no1_1[j]=exp(no1_1[j])  
      nn2_1[j]=1.e17
      on2_tgcm[j,m]=no1_1[j]/nn2_1[j]
      endif else on2_tgcm[j,m]=missing
   endfor   ; latitude loop
endfor  ; time loop
; calculate daily average for each latitude (to make latitude-dayofyear map of ON2)
for j=0,nlat-1 do begin
   index=where(on2_tgcm[j,*] ne missing,count)
   if count gt 0 then on2_tmp[j]=mean(on2_tgcm[j,index]) else on2_tmp[j]=missing
endfor
return
end