; some constants
missing=1.e36

histfile='pyr079_min_less.nc'
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           ; model time in minutes

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

o1=d.o1
o2=d.o2
tn=d.tn
gzint=d.z
no=d.no
n4s=d.n4s
n2d=d.n2d

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

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

ind=where(lat eq 2.5)
ilat=ind[0]
ind=where(lon eq 165.)
ilon=ind[0]
itime=10
; calculate number density
o_nd=o1[ilon,ilat,*,itime]*nt[ilon,ilat,*,itime]*mbar[ilon,ilat,*,itime]/16
o2_nd=o2[ilon,ilat,*,itime]*nt[ilon,ilat,*,itime]*mbar[ilon,ilat,*,itime]/32
n2_nd=(1.-o1[ilon,ilat,*,itime]-o2[ilon,ilat,*,itime])*nt[ilon,ilat,*,itime]* $
            mbar[ilon,ilat,*,itime]/28
no_nd=no[ilon,ilat,*,itime]*nt[ilon,ilat,*,itime]*mbar[ilon,ilat,*,itime]/30
n4s_nd=n4s[ilon,ilat,*,itime]*nt[ilon,ilat,*,itime]*mbar[ilon,ilat,*,itime]/14
n2d_nd=n2d[ilon,ilat,*,itime]*nt[ilon,ilat,*,itime]*mbar[ilon,ilat,*,itime]/14

close,1
openw,1,'no_min_new.dat'
fmt='$(f7.2,6e11.3,f8.2)'
for k=0,nlev-1 do begin
    printf,1,fmt,gzmid[ilon,ilat,k,itime],o_nd[k],o2_nd[k],n2_nd[k], $
                n4s_nd[k], n2d_nd[k],no_nd[k],$
                tn[ilon,ilat,k,itime]
endfor
close,1

end