program edlu_bndf c*************************************************** c This program simulated flare energy release as c avalanches in a 3D cellular automaton system c slowly driven to self-organized criticality. c c The code can run two versions of the model, c respectively taken from c c Lu, E.T., & Hamilton, R.J., ApJL 380, L89-L92 (1991) [LH91] c c Lu, E.T., Hamilton, R.J., McTiernan, J.M., & Bromund, K.R., c ApJ 412, 841-852 (1993) [LHMB93] c c and is currently set up to run the LHMB93 model c==================================================== c Names and dimensions of main variables: c c nl : dimension of lattice c niter : number of time-like iteration in run c bnm1(i,x,y,z) : i th field component at spatial node (x,y,z) c bn(i,x,y,z) : storage array for field increments at node (x,y,z) c angcrit : critical field gradient [ B_c in LH91 ] c en(niter) : time series of energy release c elatt(niter) : time series of total lattice energy c ecount(niter) : time series of number of avalanching nodes c irst : I/O flag; =0 --> start from zero field; c =1 --> start from previous run c ioutput : screen output flag; c =1 --> prints energy release info at each iteration c xfrc : forcing factor on x directions at upper/lower boundaries. c*************************************************** c Feb 2000 version; last modified 7 June 2000 c*************************************************** c This version incorporates forcing on upper and lower boundaries. c June 20, 2000. c************************************************** implicit none integer nl,ni parameter(nl=16,ni=500000) c parameter(nl=16,ni=1998791) real bn,bnm1,delb,bnn dimension bn(3,nl,nl,nl),bnm1(3,nl,nl,nl),delb(3),bnn(3) real en,estep,angcrit,angc2,ang2,ang,delb7,elatt,sum1, + bnxyz,dfac,sumb1,sumb2,sumoff,sume,efall,xfrc,diffac integer ecount,i,ix,iy,iz,iter,niter,iexcd,inn,seed,irst,nl0, + niter0,ix1,ix2,iy1,iy2,iz1,iz2,ioutput,iiter dimension en(ni),ecount(ni),elatt(0:ni),efall(0:ni) character*50 infile,outfile,outfilei real urand external urand niter=ni dfac=1. diffac = 1. angcrit=7. angc2=angcrit*angcrit xfrc = 0.1 c I/O: c infile: used only if run starts from a previous run (irst=1) c outfile: time series of energy release and lattice energy irst=0 infile='avrun16.bndf.i3e' outfile='avbf16.i3e' outfilei='avbf16_i.i3e' ioutput=1 close(unit=1) close(unit=2) close(unit=11) c---------- initialize random number generator c seed=54321 call urand_init(seed) c---------- initialize fields if(irst.eq.0)then write(*,11) c --> begin from scratch: start from zero field do ix=1,nl do iy=1,nl do iz=1,nl do i=1,3 bnm1(i,ix,iy,iz)=0. bn(i,ix,iy,iz)=0. enddo enddo enddo enddo elatt(0)=0. else c --> use field configuration from previous run write(*,12) open(2,file=infile,form='unformatted') read(2) niter0,nl0 if(nl0.ne.nl)then write(*,13) nl0,nl stop ' *ERR*' endif read(2) (en(iter),iter=1,niter0) read(2) (ecount(iter),iter=1,niter0) read(2) (elatt(iter),iter=1,niter0) elatt(0)=elatt(niter0) do i=1,3 do ix=1,nl0 do iy=1,nl0 read(2) (bnm1(i,ix,iy,iz),iz=1,nl0) do iz=1,nl0 bn(i,ix,iy,iz)=0. enddo enddo enddo enddo endif 11 format(1x,'Starting from uniform field, B_x=B_y=B_z=1') 12 format(1x,'Starting from previous solution') 13 format(1x,'ERROR: Input array incompatible; nl0= ',i3,' nl= ',i3) c---------- set x,y,z range of grid to be scanned; c on first pass (now!) scan whole grid; later, this range c will change to scan whole grid only if avalanche is in c progress, otherwise only nodes neighbouring perturbed nodes c will be scanned (for efficiency); more on this further below. ix1=2 iy1=2 iz1=2 ix2=nl-1 iy2=nl-1 iz2=nl-1 c********* This is the main time-like iteration loop do iter=1,niter c estep: energy released during this iteration c ( later stored in array en(ni) ) c iexcd: number of nodes having exceeded instability threshold c ( later stored in array ecount(ni) ) c sume: total magnetic energy of lattice c ( later stored in array elatt(ni) ) c sumoff: magnetic energy zeroed out at boundaries c ( later stored in array efall(ni) ) estep =0. iexcd =0 sume =0. sumoff =0. c========== first loop over all lattice sites to check c for stability and compute energy transfer and losses c stability is first checked at all nodes, and if threshold c exceeded the corresponding field increments at local and c neighbouring nodes are icomputed and accumulated in array bn. c Then the lattice is updated simultaneously at all nodes c in a second lattice loop. do ix=ix1,ix2 do iy=iy1,iy2 do iz=iz1,iz2 c---------- Compute deltaB c first compute average of neighbours and gradient; c and accumulate contributions to lattice energy (sume) on the fly do i=1,3 bnn(i)=0. bnxyz=bnm1(i,ix,iy,iz) do inn=-1,1,2 bnn(i)=bnn(i) + +bnm1(i,ix+inn,iy,iz) + +bnm1(i,ix,iy+inn,iz) + +bnm1(i,ix,iy,iz+inn) enddo c [ Following is eq (1) in LH91 ; also eq (1) in LHMB93 if w_j=1/6 in eq (3)] delb(i)=bnxyz-bnn(i)/6. enddo c---------- now compute critical gradient [ |dB_i| in LH91 ] c (what is computed here is the square of critical angle, c to save on a sqrt operation) ang2=0. do i=1,3 ang2=ang2+delb(i)**2 enddo c---------- if in excess of critical angle, compute contribution c to neighbours, zero angle at site, accumulate energy released c field increments are a accumulated in array bn(3,nl,nl,nl) if(ang2.gt.angc2)then c [ Use following line only for LHMB93 redistribution rule ] ang=sqrt(ang2) do i=1,3 c [ Use following line for LH91 redistribution rule ] c delb7=delb(i)/7. c [ Use following line for LHMB93 redistribution rule ] delb7=delb(i)/7. *angcrit/ang*diffac do inn=-1,1,2 bn(i,ix+inn,iy,iz)=bn(i,ix+inn,iy,iz)+delb7 bn(i,ix,iy+inn,iz)=bn(i,ix,iy+inn,iz)+delb7 bn(i,ix,iy,iz+inn)=bn(i,ix,iy,iz+inn)+delb7 enddo bn(i,ix,iy,iz)=bn(i,ix,iy,iz)-6.*delb7 enddo c---------- contribution of this node to energy release iexcd=iexcd+1 c [ Use following line for LH91 redistribution rule ] c estep=estep+angc2 c [ Use two following lines for LHMB93 redistribution rule ] estep=estep+(2.*ang/angcrit-1.) endif c========== end of first lattice loop enddo enddo enddo if(iexcd.gt.0)then c========== second lattice loop: if avalanche has occurred (iexcd>0), c update field at all lattice nodes and zero out increment c array for next iteration c Since an avalanche is in progress, we'll update all nodes c at this iteration and scan all nodes at the next iteration: ix1=2 iy1=2 iz1=1 ix2=nl-1 iy2=nl-1 iz2=nl do ix=ix1,ix2 do iy=iy1,iy2 do iz=iz1,iz2 c---------- update field components; compute magnetic energy do i=1,3 bnm1(i,ix,iy,iz)=bnm1(i,ix,iy,iz)+bn(i,ix,iy,iz) bn(i,ix,iy,iz)=0. sume=sume+bnm1(i,ix,iy,iz)**2 enddo enddo enddo enddo c---------- enforce boundary conditions and compute falloff energy, c defined as magnetic energy zeroed out at boundary do i=1,3 do iy=2,nl-1 do iz=2,nl-1 sumoff=sumoff+bn(i,1,iy,iz)**2+bn(i,nl,iy,iz)**2 bn (i,1 ,iy,iz)=0. bn (i,nl,iy,iz)=0. enddo enddo do ix=2,nl-1 do iz=2,nl-1 sumoff=sumoff+bn(i,ix,1,iz)**2+bn(i,ix,nl,iz)**2 bn (i,ix,1 ,iz)=0. bn (i,ix,nl,iz)=0. enddo enddo do ix=2,nl-1 do iy=2,nl-1 bn (i,ix,iy,1 )=0. bn (i,ix,iy,nl)=0. enddo enddo enddo c---------- compute and store time series quantities c bring back numerical factors in expressions for energy release (estep) c [ Use following line for LH91 redistribution rule ] c estep=estep*6./7. c [ Use following line for LHMB93 redistribution rule ] estep=estep*angc2*6./7. c all nodes have been scanned at this iteration, thus sume c is the total lattice energy elatt(iter)=sume efall(iter)=sumoff en(iter)=estep ecount(iter)=iexcd else c========== avalanche has NOT occurred in this iteration; therefore c add perturbation to random site in lattice lattice (excluding c boundaries) c---------- select random location, excluding boundary nodes ix=2+nint( urand()*(nl-3.) ) iy=2+nint( urand()*(nl-3.) ) iz=1+nint( urand()*(nl-1.) ) c---------- add random perturbation to field at that location; c compute energy before addition (sumb1) and after (sumb2) c to calculate increment in lattice energy; forcing c has non-zero mean (see LH91 and LHMB93) sumb1=0. sumb2=0. c******Add upper and lower boundary forcing******** if(iz.eq.1) then sumb1=sumb1+bnm1(1,ix,iy,iz)**2 bnm1(1,ix,iy,iz)=bnm1(1,ix,iy,iz) + +dfac*(urand()-xfrc) sumb2=sumb2+bnm1(1,ix,iy,iz)**2 do i=2,3 sumb1=sumb1+bnm1(i,ix,iy,iz)**2 bnm1(i,ix,iy,iz)=bnm1(i,ix,iy,iz) sumb2=sumb2+bnm1(i,ix,iy,iz)**2 enddo else if(iz.eq.nl) then sumb1=sumb1+bnm1(1,ix,iy,iz)**2 bnm1(1,ix,iy,iz)=bnm1(1,ix,iy,iz) + +dfac*(urand()+xfrc) sumb2=sumb2+bnm1(1,ix,iy,iz)**2 do i=2,3 sumb1=sumb1+bnm1(i,ix,iy,iz)**2 bnm1(i,ix,iy,iz)=bnm1(i,ix,iy,iz) sumb2=sumb2+bnm1(i,ix,iy,iz)**2 enddo else do i=1,3 sumb1=sumb1+bnm1(i,ix,iy,iz)**2 bnm1(i,ix,iy,iz)=bnm1(i,ix,iy,iz) + +dfac*(urand()) sumb2=sumb2+bnm1(i,ix,iy,iz)**2 enddo endif c---------- compute and store time series quantities elatt(iter)=elatt(iter-1)+(sumb2-sumb1) efall(iter)=0. en(iter)=0. ecount(iter)=0 c At next iteration, we only need to check stability one node c away in each direction from the newly perturbed node, hence: ix1=max(2,ix-1) ix2=min(nl-1,ix+1) iy1=max(2,iy-1) iy2=min(nl-1,iy+1) iz1=max(2,iz-1) iz2=min(nl-1,iz+1) endif c========== output to screen results at this step (if ioutput=1) if(ioutput.eq.1.and.iexcd.ge.10)then write(*,10) iter,iexcd,estep,elatt(iter),efall(iter) endif 10 format(1x,'Step:',i8,1x,'Av:',i5,1x, + 'E rel.:',e10.2,' Latt.E:',e12.4,' E off:',e10.2) c========== safety check : terminate run in case of blowup if(estep.ge.1.e20) stop ' BLOWUP' if(iter.eq.498558) then c open(11,file=outfilei,form='formatted') write(11,*) iter,nl c write(11,*) (en(iiter) ,iiter=1,iter) c write(11,*) (ecount(iiter),iiter=1,iter) c write(11,*) (elatt(iiter) ,iiter=1,iter) c write(11,*) (efall(iiter) ,iiter=1,iter) c do i=1,3 c do ix=1,nl c do iy=1,nl c write(11,*) (bnm1(i,ix,iy,iz),iz=1,nl) c enddo c enddo c enddo endif c********** end of timelike iteration loop enddo c---------- store time series to file open(1,file=outfile,form='unformatted') write(1) niter,nl write(1) (en(iter) ,iter=1,niter) write(1) (ecount(iter),iter=1,niter) write(1) (elatt(iter) ,iter=1,niter) write(1) (efall(iter) ,iter=1,niter) c---------- store current field configuration, to permit restart do i=1,3 do ix=1,nl do iy=1,nl write(1) (bnm1(i,ix,iy,iz),iz=1,nl) enddo enddo enddo close(unit=1) close(unit=2) close(unit=11) end c******************************************************************** function urand() c c Return the next pseudo-random deviate from a sequence which is c uniformly distributed in the interval [0,1] c c Uses the function ran2, from Press et al, Numerical Recipes, c 2nd ed., Cambridge Univ. Press, 1992. c implicit none c c Input - none c c Output real urand c c Local integer iseed real ran2 external ran2 c c Common block to make iseed visible to urand_init (and to save c it between calls) common /urand_seed/ iseed c urand = ran2( iseed ) return end c******************************************************************* subroutine urand_init( urand_seed ) c c Initialize random number generator urand with given seed c implicit none c c Input integer urand_seed c c Output - none c c Local integer iseed c c Common block to communicate with urand common /urand_seed/ iseed c c Seet the seed value iseed = urand_seed return end c****************************************************************** FUNCTION ran2(idum) INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV REAL ran2,AM,EPS,RNMX PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) INTEGER idum2,j,k,iv(NTAB),iy SAVE iv,iy,idum2 DATA idum2/123456789/, iv/NTAB*0/, iy/0/ if (idum.le.0) then idum=max(-idum,1) idum2=idum do 11 j=NTAB+8,1,-1 k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif k=idum/IQ1 idum=IA1*(idum-k*IQ1)-k*IR1 if (idum.lt.0) idum=idum+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idum if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return END