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

;
;here lev is model interface. in older version,the variable is "lev", in current version, 
; it is 'ilev'. So the main program that calls mmr2nd needs to pass in 'ilev'.
;o1,o2,tn are arrays of (nlon,nlat,nlev,ntime)
; some constants (in cgs units)
boltz= 1.38e-16           ; Boltzman constant
avo=6.023e23              ; Avogadro's number
;mp=1.67e-24               ; proton rest mass
p0=5.e-4                  ; reference pressure in tiegcm1.6
missing=1.e36
dgtr=1.74533e-2

;calculate effective earth radius and gravitation acceralation at earch's surface
r0=fltarr(nlat)
g0=fltarr(nlat)
for j=0,nlat-1 do begin
   c2 = cos(2.*dgtr*lat(j))
   g0[j]=980.616*(1.-.0026373*c2)
   r0[j] = 2.*g0[j]/(3.085462e-6 + 2.27e-9*c2)
endfor

;n2 mmr
n2=1-o1-o2
; 
nlev=n_elements(lev)
dz=(lev[nlev-1]-lev[0])/(nlev-1)
g=fltarr(nlev)                        ; gravitational accelaration
zmid=fltarr(nlev)                     ; pressure level at midpoint

; define output array
nt=fltarr(nlon,nlat,nlev,ntime)       ; total number density
rhot=fltarr(nlon,nlat,nlev,ntime)     ; total mass density
gzmid=fltarr(nlon,nlat,nlev,ntime)    ; altitude at midpoint
mbar=fltarr(nlon,nlat,nlev,ntime)     ; mean molecular weight

; 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.


; gzint read from tiegcm1.6 history file is calculated based on
; constant g, therefore, we calculate gzint here with g varying 
; with height
; calculate altitude at model half level

; get half model level array 
for k=0,nlev-2 do begin
    zmid(k)=0.5*(lev(k)+lev(k+1))
endfor
; pressure at model half levels
pmid=fltarr(nlev)
pmid=p0*exp(-zmid)

for m=0,ntime-1 do begin
  for i=0, nlon-1 do begin
    for j=0,nlat-1 do begin
       ; geopotential height at model half level,gzmid(nlev)

       g(0)=g0[j]*(r0[j]/(r0[j]+0.5*(gzint[i,j,0,m]+ $
             gzint[i,j,1,m])))^2                      
       ; mean molecular weight, mbar=fltarr(nlev)
       mbar(i,j,*,m)=1./( o1[i,j,*,m]/16. +  $
          o2[i,j,*,m]/32. +  $
          n2[i,j,*,m]/28.)
       ; total number density, nt=fltarr(nlev)
       nt[i,j,*,m]=pmid/(boltz*tn[i,j,*,m])
       rhot[i,j,*,m]=nt[i,j,*,m]*mbar[i,j,*,m]/avo

       mbar_1=mbar(i,j,*,m)/avo
       for k=1,nlev-1 do begin
          gzint[i,j,k,m]=gzint[i,j,k-1,m] + $
                boltz*tn[i,j,k-1,m]*dz/                    $
               (mbar_1[k-1]*g[k-1])
           if k ne nlev-1 then g[k]=g0[j]*(r0[j]/(r0[j]+0.5*(gzint[i,j,k,m]+ $
                   gzint[i,j,k+1,m])))^2
        endfor
        ; now calculate the geometric height at midpoint
        for k=0,nlev-2 do begin
             gzmid(i,j,k,m)=0.5*(gzint[i,j,k,m]+ $
                 gzint[i,j,k+1,m])
             gzmid(i,j,k,m)=gzmid(i,j,k,m)*1.e-5    ; cm to km
        endfor
     endfor
  endfor
endfor

return
end