PRO ozdens,hkm,o3dens,nlayers=nlayers,del_h=del_h,debug=debug

; This little routine returns a nominal ozone profile (CIRA-86, average
; ozone for December, 50 deg N) interpolated to a specified regular grid
; in altitude. It is modelled on NeutDens.PRO.
;
; NOTE: the ozone density returned is molecules per cubic meter, not per
; cubic centimeter.

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

k=1.380658D-23  ; Boltzmann's constant

; Pressure level (mb)
tpl=[20.,15.,10.,7.0,5.0,3.0,2.0,1.5,1.0,0.7,0.5 $
    ,0.3,0.2,0.15,0.1,0.07,0.05,0.03,0.02,0.015,0.01,0.007,0.005,0.003]

; Average ozone (ppmv) for December at 60 deg N, from p. (12)344 of CIRA-86,
to3mr=[4.77,5.10,5.42,5.69,5.94,6.10,5.79,5.10,3.74,2.68,1.94 $
; extended upward using data from 50 deg N since no data for 60 deg
      ,1.78,1.25,0.99,0.70,0.51,0.39,0.31,0.27,0.24,0.21,0.33,0.57,0.96]

; Reference pressure levels (mb) for tabulation of geopotential height,
; from p. (12)405 of CIRA-86
plref=[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.764,6.826,5.316,4.140,3.224 $
      ,2.511,1.956,1.523,1.186,0.9237,0.7194,0.5603,0.4363 $
      ,0.3398,0.2647,0.2061,0.1605,0.1250,9.736E-2,7.582E-2,5.905E-2 $
     ,4.599E-2,3.582E-2,2.789E-2,2.172E-2,1.692E-2,1.318E-2,1.026E-2 $
     ,7.992E-3,6.22E-3,4.85E-3,3.78E-3,2.94E-3]

o3mr=Interpol(to3mr,tpl,plref)>1E-3

; CIRA-86 grand mean geometric height (m), from p. (12)9 of CIRA-86
zgeom=[57.,2166.,4197.,6148.,8011.,9771.,11450.,13055. $        ; meters
      ,14604.,16116.,17625.,19151.,20697.,22267.,23860.,25478. $
      ,27120.,28788.,30479.,32198.,33948.,35737.,37563.,39425. $
      ,41325.,43265.,45231.,47213.,49198.,51170.,53116.,55029. $
      ,56905.,58739.,60539.,62308.,64048.,65754.,67433.,69087. $
      ,70720.,72333.,73932.,75521.,77097.,78662.,80217.,81766. $
      ,83298.,84804.,86284.,87733.]

; December zonal mean temperature (K) for 60 deg N,
; from p. (12)381 of CIRA-86
temp=[268.7,257.7,249.1,237.8,226.5,219.4,217.6,217.0 $
     ,216.9,216.6,216.0,215.3,214.2,212.4,210.9,209.3 $
     ,208.6,209.1,210.4,212.1,214.1,217.0,220.7,225.0 $
     ,229.7,234.8,239.6,244.5,249.0,252.5,254.1,253.6 $
     ,251.9,249.8,247.6,245.4,243.2,241.3,239.2,237.0 $
     ,235.0,233.2,231.4,229.5,227.7,226.0,224.6,223.6 $
     ,222.6,221.6,220.0,217.1]

; Determine ozone partial pressures (nb) at pressure levels 'pl'
o3pp=o3mr*plref
o3geom=1E-9*o3pp*1E5/(k*temp)   ; nb -> b -> Pa

; Extend downward with Dick Gattinger's regular tabulation.
dgoz=[1.0,1.2,1.3,1.0,.5,.3,.3,.5,1.3,2.3,3.8,4.4,4.8 $
     ,4.5,3.7,3.0,2.2,1.7,1.2,.8,.4,.3,.2,.15,.1]*1.0D18 ; /m^3
hkm=Findgen(nlayers)*del_h+del_h/2.

o3dens=Interpol(o3geom,zgeom,hkm*1000.)

; Use Dick's numbers where they are greater than the CIRA-86 ones, below
; the CIRA peak.
wo3pk=Where(o3dens EQ Max(o3dens))
o3dens(0:wo3pk(0))=o3dens(0:wo3pk(0))>dgoz(0:wo3pk(0))

Return
END