; file: makefigbetasstage.pro = main to produce departure coefficient figure ; init: Feb 13 2012 ; last: Nov 11 2016 Rob Rutten Deil ; note: first runrh on the desired XXX.atmos and XXX.atom RH setup ; adapt RH run dir doinition.pro ; then in new IDL session first @doinitrh, @doinition or program ; then .r makefigbetasstage ; parameter-free except xrange, yrange axis, label placement tuning etc ; ========== parameter business (model dependent) ; axis stuff xrange=[-100.,2300.] ; standard ALC7 for ssx xrange=[-20,1480] xtickinterval=200 xminor=2 yrange=[-2,+2] ; ALC7 yrange=[-5,+2] ; VALIIIC ; get populations, compute departure coefficients sizelevelselec=size(levelselection) nlevels=sizelevelselec[1] n=fltarr(nh,nlevels) nstar=fltarr(nh,nlevels) b=fltarr(nh,nlevels) for ilevel=0,nlevels-1 do begin level=levelselection[ilevel] n[*,ilevel]=(*thisatom.n_ptr)[*,level+ionstage*nelemicont] nstar[*,ilevel]=(*thisatom.nstar_ptr)[*,level+ionstage*nelemicont] b[*,ilevel]=n[*,ilevel]/nstar[*,ilevel] endfor ; compute Menzel format if (menzel eq 1) then begin if (ionstage eq 0) then begin ncont=(*thisatom.n_ptr)[*,nelemicont] nstarcont=(*thisatom.nstar_ptr)[*,nelemicont] endif if (ionstage eq 1) then begin ncont=(*thisatom.n_ptr)[*,nelemiicont] nstarcont=(*thisatom.nstar_ptr)[*,nelemiicont] endif for ilevel=0,nlevels-1 do b[*,ilevel]=b[*,ilevel]/(ncont/nstarcont) endif ; make betas figure ; ----------------- if (opengv) then spawn,'pkill -s 0 gv' if (ionstage eq 0) then stagename='I' if (ionstage eq 1) then stagename='II' psfilename='fig-betas-'+modelname+'-'+ionelementname+stagename+'.ps' openpsplot,psfilename,thick=3 xtitle='height [km]' if (menzel eq 1) then ytitle='log b (Menzel)' else $ ytitle='log b' logb=alog10(b) ;; xrange=[min(heightkm),max(heightkm)] ;; xrange=[-100,2100] ;; if (n_elements(maxplotheight) eq 0) then maxplotheight=max(heightkm) ;; if (maxplotheight eq 0) then maxplotheight=max(heightkm) ;; xrange=[min(heightkm),maxplotheight] ;; ymin=min(logb) ;; ymax=max(logb) ;; yrange=[ymin-0.1*(ymax-ymin),ymax+0.1*(ymax-ymin)] ;; yrange=[-3.5,+3.5] ; same as popdeps ms ;; yrange=[-0.75,+0.3] ; same as popdeps ms ;; yrange=[-2,+4] plot,heightkm,logb[*,0],$ position=[0.2,0.2,0.95,0.95],/normal,$ ; set margins around plot xticklen=0.03,yticklen=0.03/1.6,$ ; same-length ticks xtitle=xtitle,ytitle=ytitle,$ xtickinterval=xtickinterval,xminor=xminor,$ xrange=xrange,xstyle=1, yrange=yrange,ystyle=1,yminor=1 ; overplot others for ilevel=1,nlevels-1 do oplot,heightkm,logb[*,ilevel] ; add level n (starting at 1) ; ####### adjust labelind=nh/2 ;; labelind=95 ;; for ilevel=0,nlevels-2 do $ ;; xyouts,heightkm[labelind],logb[labelind,ilevel]+0.02,$ ;; string(levelselection[ilevel]+1,format='(i3)') xyouts,heightkm[labelind],logb[labelind,nlevels-1]+0.1,'cont' xyouts,heightkm[labelind],logb[labelind,0]+0.04,'1' xyouts,heightkm[labelind],logb[labelind,1]-0.23,'2' ; add model name xyouts,0.25,0.85,/norm,modelname ; add element name xyouts,0.25,0.78,/norm,ionelementname+' '+stagename ; add run dir if (addrundir eq 1) then begin cd,c=thedir lastdir=strsplit(thedir,'/',/extract) sizelastdir=size(lastdir) xyouts,0.01,0.01,/norm,lastdir[sizelastdir[1]-1] endif ; close plot closepsplot,psfilename,opengv=opengv end