PRO ozdens,hgeom,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.

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)
pl=[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,
o3mr=[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]

; December zonal mean geopotential height (m) for 60 deg N,
; from p. (12)405 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.]

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

; Tabulation of geopotential and geometric heights (at 60 deg lat) from 
; Houghton, The Physics of Atmospheres, p. 227.
; geopm=[1E4,2E4,3E4,4E4,5E4,1E5,2E5,3E5,4E5,5E5,6E5]
; geomm=[9996.,20024.,30083.,40174.,50297.,101395.,206071.,314191. $
;       ,425927.,541465.,661000.]

; Get the geometric heights corresponding to the geopotential heights at
; which the temperature is tabulated.
; zgeom=Interpol(geomm,geopm,zgeop)

; Calculate the temperature gradient over the given height range
; tempgrad=Deriv(zgeom/1E3,temp)

; Determine ozone partial pressures (nb) at pressure levels 'pl'
o3pp=o3mr*pl
o3dens=1E-9*o3pp*1E5/(k*temp)   ; nb -> b -> Pa
o3dens=o3dens/1E6               ; m^-3 -> cm^-3

; Determine geometric heights at pressure levels 'pl'
o3gh=Interpol(zgeom,plref,pl)

; Plot the results
IF Keyword_Set(ps) THEN psopen,'\usask\893\o3dist.ps',/Long $
                   ELSE Window,0,XSize=500,YSize=750
; !P.Multi=[0,0,2,0,0]
; Plot,tempgrad,zgeom/1E3 $
;     ,XTitle='Temperature Gradient (K/km)' $
;     ,YRange=[0,50],YTitle='Geometric Height (km)' $
;     ,Title='CIRA-86 Vertical Temperature Gradient  December, 60 Deg N'
; TimeStmp

Plot,o3dens,o3gh/1E3 $
    ,XRange=[0,250],XTitle='Ozone Number Density (cm!e-3!n)' $
    ,YRange=[0,80],YTitle='Geometric Height (km)' $
    ,Title='CIRA-86 Ozone Distribution  December, 60 Deg N'

; Label CIRA-86 ozone profile
XYOuts,40,35,'CIRA-86'

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.0D12
hkm=Findgen(nlayers)*del_h+del_h/2.
; ndens=Interpol(ngeom,hgeom,hkm*1000.)

OPlot,dgoz,hkm,Psym=1

;; Add data from Halley Bay, Antarctica, Aug. 15, 1987
;; (Gardiner, Antarctic Ozone Depletion Measured by Balloon Sondes at
;; Halley Bay, pp. 1-9 of _Ozone in the Atmosphere_, Bojkov and Fabian
;; (eds.), A. Deepak, 1989.
;hbpl=[20.,30.,40.,50.,60.,70.,80.,90.,100.,120.,150.,200.]
;hbo3pp=[69.,93.,104.,149.,164.,157.,160.,131.,112.,86.,67.,36.]
;
;; Determine geometric heights for Halley Bay data
;hbgm=Interpol(zgeom,plref,hbpl)
;
;; Overplot Halley Bay ozone partial pressures (multiplied by 1.4 to
;; normalize to CIRA-86) for comparison with CIRA-86 values
;OPlot,1.4*hbo3pp,hbgm/1E3,LineStyle=1
;
;; Label Halley Bay data
;XYOuts,100,10,'Halley Bay (76!9%!XS), Aug. 15, 1987'
;XYOuts,100,7,'(!9X!X1.4 to match CIRA-86)'
;
TimeStmp

IF Keyword_Set(ps) THEN psclose

Return
END