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 ; ; Modified DPS 21-Dec-96: extended available height range to 81 km. IF ((nlayers-0.5)*del_h) GT 81. THEN BEGIN Message,'Density data available only up to 81 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. $ ; 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.] p_mb=[1013.,788.93,614.42,478.51,372.66,290.23,226.03,176.03 $ ; mbar ,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,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] temp=[295.1,283.1,272.1,259.9,247.0,234.6,223.6,214.8 $ ; Kelvin ,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,259.1,253.9 $ ,248.5,243.6,238.9,234.5,230.4,226.4,222.6,219.2 $ ,216.4,214.3,212.5,210.8,209.2,207.8,206.3,204.4] p_Pa=100.*p_mb ; Pascals ngeom=p_Pa/(k*temp) ; m**-3 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 PRO atenuate,debug=debug,ps=ps,eps=eps ; PATHS.PRO - Compute path lengths for NITEOWL s/c viewing ; RLG Oct/96 ; Cosmetic changes added by D P Steele, October 28/96 ; Major changes made by D P Steele, November 6, 1996 nh=40 ; Assume 40 layers, each 2 km thick ozden=DblArr(nh) tpath=DblArr(nh) dpath=DblArr(nh) ozpath=DblArr(nh) ozpath(*)=0.0 ozpaths=ozpath neutpath=ozpath no2path=ozpath no3path=ozpath re=6371.2D0 ; Earth radius, km ; centrht=Dindgen(nh)*2 +1 yn=' ' ; o3den(*)=[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 ; cm**-3 ozdens,centrht,o3den,nlayers=nh,del_h=2.,debug=debug o3den=o3den/1.0E6 ; convert from m**-3 to cm**-3 neutdens,centrht,neutden,nlayers=nh,del_h=2.,debug=debug neutden=neutden/1.0E6 noxprofile,no2iht,no2idens,no3iht,no3idens no2dens=Interpol(no2idens,no2iht,centrht) no3dens=Interpol(no3idens,no3iht,centrht) CASE 1 OF N_Elements(ps): psopen,'c:\usask\893\crossect\atenuate.ps' N_Elements(eps): psopen,'c:\usask\893\crossect\atenuate.eps',/encaps ELSE: set_plot,'win' ; Line added 961028 by D P Steele ENDCASE ; openw,1,'ozpath.txt' ; printf,1, ' Ozone Densities [10(12)] at 1,3,5, .... 49 km ' ; printf,1,o3den o3dens=o3den below10=Where(centrht LT 10.) o3dens(below10)=0.0 ; Assume no tropospheric O3 to do comparison. for th=0,nh-1 do begin ; Compute [O3] weighted path for 25 layers curpath=0.0 ; Start with zero length at tangent point dpath(*)=0.0 & tpath(*)=0.0 ; Clean up for h=th,nh-1 do begin ; Scan from current tangent upwards to top tpath(h)=sqrt((re+(h+1.D0)*2.D0)^2 - (re+th*2.D0)^2)*1.0D5 ; cm dpath(h)=tpath(h)-curpath ; Store layer geometric length curpath=tpath(h) ozpath(th)=ozpath(th)+o3den(h)*dpath(h) ; Strato and tropo O3 ozpaths(th)=ozpaths(th)+o3dens(h)*dpath(h) ; Strato O3 only neutpath(th)=neutpath(th)+neutden(h)*dpath(h) ; Neutral atmosphere no2path(th)=no2path(th)+no2dens(h)*dpath(h) ; NO2 no3path(th)=no3path(th)+no3dens(h)*dpath(h) ; NO3 endfor ; if th eq 0 then begin ; A check print ; printf,1,' Path Lengths (Km) and Deltas for Tangent = ',th*2+1,' Km' ; printf,1,tpath ; printf,1,dpath ; endif endfor ; printf,1,' Column O3 [cm-2] vs Tangent Height With and Without Tropo [O3]' ; printf,1,ozpath ; printf,1,ozpaths IF NOT Keyword_Set(eps) THEN BEGIN plot,ozpath,centrht,linestyle=0 $ ,title=' Column O3 Concentration (cm!e-2!n) vs Tangent Height' $ ,xtitle='Column Ozone Amount (cm!e-2!n)' $ ,ytitle='Tangent Height (km)' oplot,ozpaths,centrht,linestyle=2 ENDIF read,'Press to see the assumed ozone profiles.',yn IF NOT Keyword_Set(eps) THEN BEGIN plot,o3den,centrht,linestyle=0,title=' Assumed Ozone Profiles' $ ,xtitle='Ozone Density (cm!e-3!n)' $ ,ytitle='Height (km)' oplot,o3dens,centrht,linestyle=2 ENDIF ; close,1 Read,'Show absorption/scattering effects on a stellar spectrum? ([Y]/N) ' $ ,yn doabs=(StrUpCase(yn) NE 'N') IF doabs THEN BEGIN Read,'Magnitude of stellar spectrum? ',starmag factor=10.D0^((2.486593-starmag)*0.4) ; Read,'Altitude of tangent ray? (1..79 km) ',alt FOR alt=10,60,10 DO BEGIN ozcol=(Interpol(ozpath,centrht,alt))(0) neutcol=(Interpol(neutpath,centrht,alt))(0) no2col=(Interpol(no2path,centrht,alt))(0) no3col=(Interpol(no3path,centrht,alt))(0) IF alt EQ 10 THEN BEGIN Read,'Ozone absorption? ([Y]/N) ',yn ozabs=(StrUpCase(yn) NE 'N') Read,'Rayleigh scattering? ([Y]/N) ',yn raysc=(StrUpCase(yn) NE 'N') Read,'NO2 absorption? ([Y]/N) ',yn no2abs=(StrUpCase(yn) NE 'N') Read,'NO3 absorption? ([Y]/N) ',yn no3abs=(StrUpCase(yn) NE 'N') IF Keyword_Set(debug) THEN Stop ; Load star spectrum for G-type star OpenR,lun,'\usask\893\crossect\wavelgth.txt',/Get_LUN starwave=DblArr(1084) ReadF,lun,starwave Free_LUN,lun dwl=Median(starwave-Shift(starwave,1)) OpenR,lun,'\usask\893\crossect\g1--2v.txt',/Get_LUN starflux=DblArr(1084) ReadF,lun,starflux Free_LUN,lun starflux=starflux*factor ; Load ozone absorption cross-section from Shettle et al., GRL rdcols,'\usask\893\crossect\ozonecro.dat',no3,o3wl,o3xs o3wl=[Findgen(70)+4000.,o3wl] ; Extrapolate xsection to 4000 A. o3xs=[Replicate(2.0E-24,70),o3xs] no3=no3+70 ; Load Rayleigh scattering cross-section from Jayson Eppler (MODTRAN) rdcols,'\usask\893\crossect\ray2xsec.txt',nrs,rswl,rsxs rswl=10.*rswl ; convert to Angstroms ; Load NO2 cross-section no2xsect,no2wl,no2xs,/angstrom ; Load NO2 density profile no3lores,no3wl,no3xs,/angstrom ; Stay within wavelength range covered by all data ; wllow=Min(starwave)>Min(o3wl)>Min(rswl)>Min(no2wl)>Min(no3wl) wllow=4000. ; wlhigh=Max(starwave)