; FILE: butterfly.ivp.idl [IDL source code] ;=========================================================== ; NONLINEAR SPHERICAL MEAN-FIELD DYNAMO ;----------------------------------------------------------- ; APAS7500 Fall 1997 Problem II.X.Y ; see Section X.Y of notes, and Appendix X ;----------------------------------------------------------- ; This codes produces butterfly diagrams of B_r and ; B_phi at depths rbr and rbphi, as set by the user (see below). ; The solution is reconstructed from a solution ; produced by the fortran code sphdyn_eig.f ; ; Solution is plotted as B_phi isocontours and poloidal ; fiedlines in meridional quadrant, over half a dynamo ; cycle at phase intervals of pi/6 ;----------------------------------------------------------- ; HARDCOPY: when producing a postscript hardcopy, before ; running the code (.run medplane.idl) make sure to set ; ; device,/portrait,yoffset=5.,xsize=18,ysize=18 ;----------------------------------------------------------- ; INPUT: this code requires an input file, the name of which ; being set immediately below ;=========================================================== nnx=1200 nny=1200 window,15,retain=2,xpos=50,ypos=200,xs=nnx,ys=nny solfile='saasfe/BL_ivp.3030.i3e' parity=-1. parity=0. rbphi=0.7 ; depth of toroidal field cut rbr =1.0 loadct,0 tvlct,rr,gg,bb,/get rr(1)=255 & gg(1)=0 & bb(1)=0 ; red tvlct,rr,gg,bb a=2.*!PI*findgen(16)/15. usersym,sin(a),cos(a),/fill nskip=10 ;***** you shouldn't have to touch anything below this point ***** !ORDER=0 close,1 & openr,1,/f77_unformatted,solfile title0=string(' ',format='(a80)') & ctime=1.e0 readu,1,title0 print,'Using Input file: ',solfile print,'Run title: ',title0 idum=lonarr(6) vdum=fltarr(8) readu,1,idum numnp =idum(0) numel =idum(1) nelx =idum(2) nely =idum(3) nt0 =idum(4) nsteps=idum(5) print,idum readu,1,comega,calpha,reynolds,etadf,ybl,rzc,to ; & etadf=0.001 print,comega,calpha,reynolds,etadf,ybl,rzc,to ; & etadf=0.001 nnodes=numnp/2 nt=nt0*nsteps xx =fltarr(nelx+1) yy =fltarr(nely+1) b_r =fltarr(nt+1,nelx+1) b_y =fltarr(nt+1,nelx+1) omg =fltarr(nelx+1,nely+1) eta =fltarr(nelx+1,nely+1) ur =fltarr(nelx+1,nely+1) uth =fltarr(nelx+1,nely+1) aquad=fltarr(nelx+1) daqdm=fltarr(nelx+1) bquad=fltarr(nelx+1,nely+1) bphi=fltarr(nelx+1,nely+1) ;to=to*nt/nt0 time=findgen(nt+1)/(nt)*to ;*(1.*nt/nt0) emag=fltarr(nt+1) readu,1,xx readu,1,yy readu,1,omg readu,1,ur readu,1,uth readu,1,eta xxa=acos(xx) ;---------- find boundaries of layer and top iblay=0 ibp=0 for k=1,nely-1 do if yy(k+1) gt rbphi and yy(k) le rbphi then ibp=k for k=1,nely-1 do if yy(k+1) gt rbr and yy(k) le rbr then ibr=k print, 'Bphi Butterfly diagram constructed at r/R=',yy(ibp) print, 'Br Butterfly diagram constructed at r/R=',yy(ibr) ;----------------------------------------------------------------------- ; select contour interval lev=2.*findgen(16)/15.-1. clst=[1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0] lev=[-1.0,-0.5,-0.2,-0.1,-0.05,-0.02,-0.01,-0.005,-0.002,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.5,1.0] lev2=[-0.9,-0.7,-0.5,-0.3,-0.1,0.1,0.3,0.5,0.7,0.9] clst2=[1,1,1,1,1,0,0,0,0,0] lev3=[-0.9,-0.7,-0.5,-0.3,-0.1,-0.01,-0.001,0.001,0.01,0.1,0.3,0.5,0.7,0.9] clst3=[1,1,1,1,1,1,1,0,0,0,0,0,0,0] ;----------------------------------------------------------------------- ; plot grid pos1 =[0.15,0.425,0.95,0.675] pos2 =[0.15,0.125,0.95,0.375] pos3 =[0.15,0.725,0.95,0.975] ;----------------------------------------------------------------------------- ; loop over timesteps idum=lonarr(2) ;----------------------------------------------------------------------------- ; loop over timesteps istep=lonarr(1) vdum=fltarr(5) ifr=-1 for it=0,nt do begin ; print,it readu,1,istep readu,1,vdum & emag(it)=vdum(2) ; time(it)=vdum(0) & emag(it)=vdum(2) ; print,(iskip+1),vdum readu,1,bquad readu,1,bphi ;---------- Extract toroidal/poloidal fields at desired depths for j=0,nelx do begin ; rsin=yy(ibr)*sqrt(1.-xx(j)^2) b_y(it,j) =bphi(j,ibp) ; aquad(j) =rsin*bquad(j,ibr) aquad(j) =bquad(j,ibr) endfor ;---------- Compute vector potential angular derivative (to get B_r) ; this computes d/dmu (A sin theta) ( mu=cos(theta) ) for j=1,nelx-1 do begin ; sinthp=sin(acos(xx(j+1))) ; sinthm=sin(acos(xx(j-1))) sinthp=sin(xxa(j+1)) sinthm=sin(xxa(j-1)) ; delth=xxa(j+1)-xxa(j-1) ; daqdm(j)=(aquad(j+1)*sinthp-aquad(j-1)*sinthm)/delth delmu=xx(j+1)-xx(j-1) daqdm(j)=(aquad(j+1)*sinthp-aquad(j-1)*sinthm)/delmu endfor ; patch for pole sinthp=sin(acos(xx(1))) sinthm=sin(acos(xx(0))) delmu=xx(1)-xx(0) daqdm(0)=(aquad(1)*sinthp-aquad(0)*sinthm)/delmu ; patch for equator sinthp=sin(acos(xx(nelx))) sinthm=sin(acos(xx(nelx-1))) delmu=xx(nelx)-xx(nelx-1) daqdm(nelx)=(aquad(nelx)*sinthp-aquad(nelx-1)*sinthm)/delmu ;---------- Now finally calculate B_r at surface for j=0,nelx do begin ; rsin=1.0*sin(acos(xx(j))) ; b_r(it,j)=daqdm(j)/rsin b_r(it,j)=-daqdm(j) endfor ; thmid=0.5*(xxa(1)+xxa(0)) ; rsinthmid=1.0*sin(thmid) ; b_r(it,0)=daqdm(0) ; thmid=0.5*(xxa(nelx)+xxa(nelx-1)) ; rsinthmid=1.0*sin(thmid) ; b_r(it,nelx)=daqdm(nelx)/rsinthmid if parity lt 0. then b_r(it,0)=0. ; if parity gt 0. then b_r(it,0)=b_r(it,nelx-1) endfor ;---------- Now normalize fields to [-1.,1.] maxbr=max(b_r(nt/2:nt,*)) maxby=max(b_y(nt/2:nt,*)) b_r=b_r/maxbr b_y=b_y/maxby title1='B!9P!6(r/R='+string(yy(ibp),format='(f4.2)')+'), max(B!9P!6)='$ +string(maxby,format='(e8.2)') title2='Br(r/R='+string(yy(ibr),format='(f4.2)')+'), max(Br)='$ +string(maxbr,format='(e8.2)') ;title3='!17Babcock-Leighton solar-cycle model' ;title3='!17Stochastic !7a!17-quenched solar-like dynamo model'' title3='!17'+title0 ; plot butterfly diagrams ;------------------------------------------------------------------------------ !X.STYLE=1 & !X.THICK=2 !Y.STYLE=1 & !Y.THICK=2 ti=0. xxan=-(2.*xxa/!PI-1.) ;to=0.05 ;---------- 1. plot color-filled contours for B-phi contour,pos=pos1,levels=lev2,$ c_linestyle=clst2,xcharsize=0.0000000001,charsize=1.25,$ xrange=[ti,to],yrange=[0.,1.],xticklen=0.04,$ ytitle='!6(latitude)/(!7p!6/2)',$ ; ytitle='!17cos!7(h)!17',$ title='!6 ',$ b_y,time,xxan ; b_y,time,xx if parity ne 0. then begin contour,/noerase,pos=pos1,levels=lev2,$ c_linestyle=clst2,xcharsize=0.0000000001,charsize=1.25,$ xrange=[ti,to],yrange=[0.,1.],xticklen=0.04,$ ytitle='!6(latitude)/(!7p!6/2)',$ ; ytitle='!17cos!7(h)!17',$ title='!6 ',$ parity*b_y,time,-xxan ; parity*b_y,time,-xx endif ; xyouts,0.,1.5,title0 xyouts,ti+0.,1.1,title1 for ii=0,nt-1 do begin dum=max(abs(b_y(ii:ii,0:nelx/2)),imax) oplot,psym=8,symsize=0.5,[time(ii),time(ii)],[xxan(imax),xxan(imax)],color=1 ; dum=max(abs(b_y(ii:ii,nelx/2:nelx)),imax) ; oplot,psym=8,symsize=0.5,[time(ii),time(ii)],[xxan(imax),xxan(imax)],color=1 endfor contour,/noerase,pos=pos2,levels=lev3,$ c_linestyle=clst3,charsize=1.25,$ xrange=[ti,to],yrange=[0.,1.],xticklen=0.04,$ ytitle='!6(latitude)/(!7p!6/2)',$ ; ytitle='!17cos!7(h)!17',$ xtitle='!18t/!7s!17',$ b_r,time,xxan ; b_r,time,xx if parity ne 0. then begin contour,/noerase,pos=pos2,levels=lev3,$ c_linestyle=clst3,charsize=1.25,$ xrange=[ti,to],yrange=[0.,1.],xticklen=0.04,$ ytitle='!6(latitude)/(!7p!6/2)',$ ; ytitle='!17cos!7(h)!17',$ xtitle='!18t/!7s!17',$ parity*b_r,time,-xxan ; parity*b_r,time,-xx endif reylab='!6Rm='+string(reynolds,format='(e8.1)') wlab='!18w/R!17='+string(ybl,format='(f5.3)') etalab='!7g!D!17c!n!17/!7g!D!17e!n!17='+string(etadf,format='(e8.1)') colab='!6C!7X!6='+string(comega,format='(e8.1)') calab='!6C!7a!6='+string(calpha,format='(f5.1)') xyouts,ti+0.,1.1,title2 xyouts,ti+0.,-1.45,calab xyouts,ti+0.,-1.65,colab xyouts,ti+0.,-1.85,reylab xyouts,ti+2.*(to-ti)/3.,-1.45,wlab xyouts,ti+2.*(to-ti)/3.,-1.65,etalab if nt lt 100 then ssym=4 if nt ge 100 then ssym=3 ssym=-4 ;yr=[0.,max(emag(nt/2:nt-2))] yr=[0.,max(emag(nt/2:nt-2))] yr=[min(emag(0:nt-2)),max(emag(0:nt-2))] plot,/noerase,pos=pos3,time,emag,psym=ssym,yrange=yr,$ xcharsize=0.00000001,xrange=[ti,to],xticklen=0.04,$ ytitle='!18E!DB!b!17',title=title3 ;endfor end