pro read_3dctm,n,filename,head_all,press,psurf_all,press3d,array_all, $
               lowres=lowres

nhead=75
nlon=144
nlat=89
nlev=28

if lowres ne 0 then begin

nhead=75
nlon=72
nlat=46

endif

head=fltarr(nhead)
psurf=fltarr(nlon,nlat)
array=fltarr(nlon,nlat,28)

head_all=fltarr(nhead,n)
psurf_all=fltarr(nlon,nlat,n)
array_all=fltarr(nlon,nlat,28,n)


get_lun,lun

openr,lun,filename,f77_unformatted

for i=0,n-1 do begin
   readu,lun,head,psurf,array
   head_all(*,i)=head
   psurf_all(*,*,i)=psurf
   array_all(*,*,*,i)=array
endfor

close,lun
free_lun,lun

;** calculate 3d pressure array from specified pressure levels

; model pressure level edges:

prslay=[1000.00, 843.304, 707.036, 595.789, 494.923 $
         , 393.598, 308.007, 247.462, 204.531, 174.100, 150.415 $
         , 130.000, 108.628, 93.8084, 82.3165, 72.4290, 62.3167 $
         , 51.4364, 38.3119, 26.1016, 17.7828, 12.1153, 8.25404 $
         , 5.62341, 3.83119, 2.61016, 1.77828, 1.00000, .430576]

; prslay is used to define the sigma levels in the following two-step
; process.  First, pressure levels are found which are the geometric
; means of the edge pressures defined above:

nlevs=n_elements(prslay)
press=sqrt(prslay(0:nlevs-2)*prslay(1:nlevs-1))

; next, calculate sigma levels corresponding to these levels. 
; In the 3dctm code, ktrop=11, kstrat=17, and levels (the total
; number of levels is ktrop + kstrat = levels = 28. prslay is 
; defined from (0,levels) so prslay(11) is the interface pressure,
; or pint = 130.00

ktrop=7
pint=prslay(ktrop)
ptop=prslay(nlevs-1)
ps=prslay(0)

; In the troposphere the sigma levels range from 0 to 1, where sigma = 1
; is the ground prslay(0) and sigma = 0 is the interface pressure:

sigma=fltarr(nlevs-1)

sigma(0:ktrop-1)=(press(0:ktrop-1)-pint)/(ps-pint)
sigma(ktrop:nlevs-2)=(press(ktrop:nlevs-2)-pint)/(pint-ptop)

; OK, now I can calculate the 3d pressure field:

press3d=fltarr(nlon,nlat,28,n)

for i=0,nlon-1 do for j=0,nlat-1 do for k=0,ktrop-1 do for l=0,n-1 do $
   press3d(i,j,k,l)=sigma(k)*(psurf_all(i,j,l)-pint)+pint

for i=0,nlon-1 do for j=0,nlat-1 do for k=ktrop,nlevs-2 do for l=0,n-1 do $
   press3d(i,j,k,l)=sigma(k)*(pint-ptop)+pint

; Now, pint=130 mbar now, and generally we are above that value
; so the model surfaces are pressure surfaces.  Below that value,
; the model surfaces are not parallel with pressure surfaces.
; Most of the time, I will be looking above 130 mbar, so most
; of the time press3d will not be necessary.

return
end