PRO neutdens,hkm,ndens,nlayers=nlayers,del_h=del_h,debug=debug

; This little routine returns neutral density (MKS units) vs. geometric
; height over the height range 0 - 50 km.  The data are taken from the
; CIRA-86 grand mean (which is nominally for 30 deg N).
;
; DPSteele, Nov. 6, 1996

IF ((nlayers-0.5)*del_h) GT 51. THEN BEGIN
    Message,'Density data available only up to 51 km at present' $
           ,/Informational
    Return
    ENDIF
    
IF Keyword_Set(debug) THEN Stop

k=1.380658D-23  ; Boltzmann's constant

hgeom=[57.,2166.,4197.,6148.,8011.,9771.,11450.,13055. $
      ,14604.,16116.,17625.,19151.,20697.,22267.,23860.,25478. $
      ,27120.,28788.,30479.,32198.,33948.,35737.,37563.,39425. $
      ,41325.,43265.,45231.,47213.,49198.,51170.]

p_mb=[1013.,788.93,614.42,478.51,372.66,290.23,226.03,176.03 $
     ,137.09,106.77,83.15,64.76,50.43,39.28,30.59,23.82 $
     ,18.55,14.45,11.25,8.7642,6.8255,5.3157,4.1399,3.2242 $
     ,2.5110,1.9556,1.5230,1.1861,0.9237,0.7194]

temp=[295.1,283.1,272.1,259.9,247.0,234.6,223.6,214.8 $
     ,208.5,205.8,207.1,210.0,212.8,215.6,218.3,221.0 $
     ,223.9,226.9,230.3,234.4,238.9,243.7,248.7,253.8 $
     ,258.8,263.0,266.0,267.1,266.5,263.5]

p_Pa=100.*p_mb
ngeom=p_Pa/(k*temp)

IF N_Elements(nlayers) EQ 0 THEN nlayers=25
IF N_Elements(del_h) EQ 0 THEN del_h=2.

hkm=Findgen(nlayers)*del_h+del_h/2.
ndens=Interpol(ngeom,hgeom,hkm*1000.)

Return
END