c     Turbulent Extreme Multi-Zone (TEMZ) model for blazar multi-waveband
c       light curves, SEDs, and (synchrotron) polarization 
c
c     Program that calculates the light curve and dynamic SED of a jet with turbulent
c       plasma passing through a conical standing shock
c     This version includes variations on times as short as the time for a cell to
c       cross the shock, as well as longer-term variations according to a given PSD slope
c     Includes synchrotron radiation and IC of seed photons from a dust torus
c       and from a Mach disk
c     In this version, a turbulent velocity field is approximated by nested zones; all
c       except single cells rotate about random axis with rotational velocity following
c       Kolmogorov spectrum
c
c     To compile:
c     gfortran -fopenmp -o temz.out temz.f
c
c     To run: need input file temzinp.txt
c
c     Depending on the system, may need to set environment variables
c       to allow large arrays. For example, in cshell put the following
c       4 lines into the .sh or .csh file that is read when you log in
c    OMP_NUM_THREADS=16   ! For 16 processors, for example
c    KMP_STACKSIZE=unlimited
c    OMP_PROC_BIND=TRUE
c    OMP_WAIT_POLICY=passive
c     Then, on the command line, type:
c    limit stacksize unlimited
c
      use OMP_LIB
c
      dimension itarra(3),spsdx(32768),spsd(32768),spsd1(32768),
     , pq(400,169,68),pu(400,169,68),fpol(400,169,68),
     , flux(400,169,68),flsync(400,169,68),flcomp(400,169,68),
     , flec(400,169,68),flssc(400,169,68),
     , gammin(400,169),gammax(400,169),gmax0(400,169),n0(400,169),
     , egam(400,169,44),enofe(400,169,44),
     , bx(400,169),by(400,169),bz(400,169),delta(400,169),
     , betadz(400,169),betadx(400,169),betady(400,169),betad(400,169),
     , gammad(400,169),zcell(400,169),layer(169),tlf1(169),
     , xcell(169),ycell(169),rcell(169),sinph(169),cosph(169),
     , phcell(169),psi(169),tanpsi(169),sinpsi(169),cospsi(169),
     , imax(169),nouter(169),bperp(169),bfield(169),
     , betaux(169),betauy(169),betauz(169),betau(169),gammau(169),
     , deltmd(169),dmd(169),betamz(169),betamr(169),gamamr(169),
     , ng1x(169),ng1y(169),ng2x(169),ng2y(169),ng3x(169),ng3y(169),
     , ididj1(169),ididj2(169),ididj3(169),
     , ggam(44),edist(44),gcnt(44),igcnt(44),
     , fsynmd(64000,68),fsscmd(64000,68),alphmd(64000,68),
     , absorb(64000,68),alfmdc(64000,68),abexmd(64000,68),
     , bmdx(64000),bmdy(64000),bmdz(64000),bmdtot(64000),mdd(64000),
     , zsmd(64000),bcmd(64000,3),tlf(400),betamx(169),betamy(169),
     , useed(401),seedpk(401),tdelr(401),dustii(401,22),
     , dustf(401,22),dustnu(22),dusti(22),
     , nu(68),snu(68),syseed(68),scseed(68),ssseed(68),fsync(68),
     , fsync2(68),phots(68),phalph(68),tauxmd(68),fmdall(68),
     , ssabsm(68),ssabs(68),sstau(68),srcfn(68),rnu(68),
     , ecflux(68),emisco(68),tauexp(68),taupp(68),
     , ajofnu(68),fsnoab(68),specin(68),sscflx(68),spxec(68),spxssc(68),
     , ididb1(8,8),ididb2(4,4),ididb3(2,2),ididc1(8,8),ididc2(4,4),
     , ididc3(2,2),ididd1(8,8,8),ididd2(4,4,8),ididd3(2,2,8),
     , ididl1(8,8,8),ididl2(4,4,8),ididl3(2,2,8),
     , ididk1(169,8),ididk2(169,8),ididk3(169,8),
     , phi11(8,8),phi21(8,8),thet11(8,8),thet21(8,8),
     , phi12(4,4),phi22(4,4),thet12(4,4),thet22(4,4),
     , phi13(2,2),phi23(2,2),thet13(2,2),thet23(2,2),
     , phv11(8,8),phv12(4,4),phv13(2,2),cosv11(8,8),cosv12(4,4),
     , cosv13(2,2),sinv11(8,8),sinv12(4,4),sinv13(2,2),
     , bfrac1(8,8),bfrac2(4,4),bfrac3(2,2),
     , bux1(8,8,400),buy1(8,8,400),buz1(8,8,400),bux2(4,4,400),
     , buy2(4,4,400),buz2(4,4,400),
     , bux3(2,2,400),buy3(2,2,400),buz3(2,2,400),
     , bu1x1(8,8),bu1y1(8,8),bu1z1(8,8),bu2x1(8,8),bu2y1(8,8),
     , bu2z1(8,8),cu1x1(8,8),cu1y1(8,8),cu1z1(8,8),angrt1(8,8),
     , cosrt1(8,8),cosrt2(4,4),cosrt3(2,2),
     , bu1x2(4,4),bu1y2(4,4),bu1z2(4,4),bu2x2(4,4),bu2y2(4,4),
     , bu2z2(4,4),cu1x2(4,4),cu1y2(4,4),cu1z2(4,4),angrt2(4,4),
     , bu1x3(2,2),bu1y3(2,2),bu1z3(2,2),bu2x3(2,2),bu2y3(2,2),
     , bu2z3(2,2),cu1x3(2,2),cu1y3(2,2),cu1z3(2,2),angrt3(2,2)
      character*10 dumdum
      character*16 filnam,filnm2
      real*8 pol,pqcum,pucum,pq,pu,pmean,polc,ai2,
     ,  pcum,tanv0,cosv0,betaup,gamup,beta,sinz,cosz,phcell,
     ,  cosph,sinph,thlos,betau,gammau,opang,tanop,cosop,sinop,
     ,  zeta,tanz,betad,betadx,betady,betadz,slos,clos,
     ,  eta,tanxi,xi,betacs,gammad,bdx,bdy,bdz,
     ,  dcsth1,dcsth2,dth1,dth2,dsnth1,dsnth2,n0ave,
     ,  betamd,betarl,betamx,betamy,betamz,betamr,gamamr,
     ,  cosmd,tanmd,bup,gmax0,ustob,tlfact,delt,
     ,  gamb,glow,gmrat,gmratl,gmratm,ggam,tloss,t1,t2,tlmin,
     ,  eterm1,eterm2,gammax,glim,t2max,phots,betaux,betauy,
     ,  betauz,psi,cospsi,sinpsi,betat,
     ,  tanpsi,betd,gamd,bd2,gm1,angm,bmparx,bmpary,bmparz,
     ,  bmprpx,bmprpy,bmprpz,dtnth1,dtnth2,dotprd,
     ,  betupx,betupy,betupz,betatx,betaty,betatz,thetat,
     ,  sintht,costht,phit,btprpx,btprpy,btprpz,
     ,  btparx,btpary,btparz,sx,sy,sz,slx,sly,slz,denomi,
     ,  vvx,vvy,vvz,uux,uuy,uuz,denomv,gvv,vvx1,vvy1,vvz1,
     ,  anx,any,anz,betadd,gammdd,csth1,csth2,tcpu1,
     ,  emisco,fsync2,ssabs,sstau,ssabsm,tauexp,ecflux,sscflx,fsnoab,
     ,  srcfn,ssseed,fmdall,tauxmd,ajofnu,absorb,abssss,fsynmd,fsyncc
      data ididb1/ 64*0/, ididb2/ 16*0/, ididb3/ 4*0/,
     ,  bfrac1/ 64*0.0/, bfrac2/ 16*0.0/, bfrac3/ 4*0.0/,
c     ,  angrt1/ 64*0.0/, angrt2/ 16*0.0/, angrt3/ 4*0.0/,
c     ,  cosrt1/ 64*0.0/, cosrt2/ 16*0.0/, cosrt3/ 4*0.0/,
     ,  ididj1/ 169*0/,ididj2/ 169*0/,ididj3/ 169*0/,
     ,  ididc1/ 64*0/,ididc2/ 16*0/,ididc3/ 4*0/,
     ,  ididd1/ 512*0/,ididd2/ 128*0/,ididd3/ 32*0/,
     ,  ididl1/ 512*0/,ididl2/ 128*0/,ididl3/ 32*0/,
     ,  ididk1/ 1352*0/,ididk2/ 1352*0/,ididk3/ 1352*0/,
     ,  bux1/ 25600*0.0/,buy1/ 25600*0.0/,buz1/ 25600*0.0/,
     ,  bux2/ 6400*0.0/,buy2/ 6400*0.0/,buz2/ 6400*0.0/,
     ,  bux3/ 1600*0.0/,buy3/ 1600*0.0/,buz3/ 1600*0.0/,
     ,  bu1x1/ 64*0.0/,bu1y1/ 64*0.0/,bu1z1/ 64*0.0/,bu2x1/ 64*0.0/,
     ,  bu2y1/ 64*0.0/,bu2z1/ 64*0.0/,cu1x1/ 64*0.0/,cu1y1/64*0.0/,
     ,  cu1z1/ 64*0.0/,angrt1/ 64*0.0/,cosrt1/ 64*0.0/,
     ,  cosrt2/ 16*0.0/,cosrt3/ 4*0.0/
      real*4 n0mean,n0,nu,ldust
      integer dstart
      character*6 :: dpath='mapsp/'
c      common/ci/i,j
      common/cvel/bdx,bdy,bdz,gamd,betd
      common/cparm/zred1,bfld,bperpp
      common/cdust/csth1,csth2,dcsth1,dcsth2,dsang,tdust
      common/cssc/snu,ssseed,cscat
      common/cdist/ggam,edist
      common/cseed/dustnu,dusti,csang
      common/cang/cosz,sinz
      common/crite/iwrite
ccc
ccc   Set nomap=0 if you want an image (position, flux, and polarization for each cell)
ccc     In this case, you need to create directory 'mapsp'
ccc
ccc     Otherwise, set nomap=1
      iwrote=0
      nomap=0
c
      open (2,iostat=ios, err=9000, file='temzinp.txt',
     ,   status='old')
      open (3,iostat=ios, err=9000, file='temzspec.txt',
     , status='new')
      open (4,iostat=ios, err=9000, file='temzlc.txt',
     , status='new')
      open (5,iostat=ios, err=9000, file='temzpol.txt',
     , status='new')
      open (6,iostat=ios, err=9000, file='temzcheck.txt',
     , status='new')
      open (9,iostat=ios, err=9000, file='ompcheck.txt',
     , status='new')
      it=0
      if(nomap.eq.1)go to 1111
      write(filnam,"(a7,i4.4,a4)") "temzmap",it,".txt"
      open (7,iostat=ios, err=9000, file=dpath//filnam,
     ,   status='new')
      write(filnm2,"(a7,i4.4,a4)") "temzmpo",it,".txt"
      open (8,iostat=ios, err=9000, file=dpath//filnm2,
     ,   status='new')
      close(7)
      close(8)
 1111 continue
c     Input file format: 1st line characters, then values of parameters,
c       one per line
      read(2,9111)dumdum,nrad,itlast,modul,zred,dgpc,alpha,pn,pb,bave,
     , psdslp,uratio,rsize,gmaxmn,gmrat,gmin,betaup,betat,zeta,thlos,
     , opang,tdust,ldust,dtdist,dtrad,zdist0,vmd,bratio,bpitch
      close(2)
c     openmp initialization
c      time1=omp_get_wtime()
      iprocs = omp_get_num_procs()
      nthread = iprocs
      call omp_set_num_threads(nthread)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 7701 inu=1,40
      idth = omp_get_thread_num()
      fsscmd(1,inu)=0.0
      write(9,9994)idth,inu,fsscmd(1,inu)
 7701 continue
!$omp end parallel do
      call itime(itarra)
      ita=itarra(1)+itarra(2)+itarra(3)
      iseed = rand(ita)
      randstart=rand(iseed)
c     Turn pair-production absorption on (ipprod=1) or off (ipprod=0)
      ipprod=0
      do 6543 jn=1,8
      do 6543 jm=1,8
      angrt1(jm,jn)=0.0
 6543 cosrt1(jm,jn)=0.0
      do 6544 jn=1,4
      do 6544 jm=1,4
      angrt2(jm,jn)=0.0
 6544 cosrt2(jm,jn)=0.0
      do 6545 jn=1,2
      do 6545 jm=1,2
      angrt3(jm,jn)=0.0
 6545 cosrt3(jm,jn)=0.0
c     jcells = numbers of cells in transverse direction (cross-section)
      jcells=3*nrad*(nrad-1)+1
      ancol=2*nrad-1
      rbound=(ancol-1.0)*rsize
c     Scaling of turbulent velocity, density, and magnetic field with scale size
      tvexp=1.0/3.0
      tnexp=2.0/3.0
      tbexp=1.0/3.0
c     Normalization of turbulent density; from Kolmogorov spectrum
      anorm0=1.0/((0.5*ancol)**tnexp+(0.25*ancol)**tnexp+
     ,   (0.125*ancol)**tnexp+1.0)
      anorm1=anorm0*(0.125*ancol)**tnexp
      anorm2=anorm0*(0.25*ancol)**tnexp
      anorm3=anorm0*(0.5*ancol)**tnexp
c     Normalization of turbulent field strengths & velocities; from Kolmogorov spectrum
      bnorm0=1.0/((0.5*ancol)**tbexp+(0.25*ancol)**tbexp+
     ,   (0.125*ancol)**tbexp+1.0)
      bnorm1=bnorm0*(0.125*ancol)**tbexp
      bnorm2=bnorm0*(0.25*ancol)**tbexp
      bnorm3=bnorm0*(0.5*ancol)**tbexp
      vnorm0=1.0/((0.5*ancol)**tvexp+(0.25*ancol)**tvexp+
     ,   (0.125*ancol)**tvexp+1.0)
      vnorm1=vnorm0*(0.125*ancol)**tvexp
      vnorm2=vnorm0*(0.25*ancol)**tvexp
      vnorm3=vnorm0*(0.5*ancol)**tvexp
      write(6,9433)nrad,ancol,bnorm0,bnorm1,bnorm2,bnorm3,anorm0,
     ,  anorm1,anorm2,anorm3,vnorm0,vnorm1,vnorm2,vnorm3
 9433 format(i4,f5.0,1p12e10.4)
c     zrat later defines zsize = zrat*rsize/tanz
c       If zrat is decreased, size of dust arrays idsize needs to be increased
      zrat=0.25
      amppsd=1.0
      ndim=32768
      andim=ndim
      mdmd=500
c     (mdmax-1) should be evenly divisible by 8
      mdmax=32001
c     An SED will be printed out every ispecs time steps
      ispecs=1
      ispec=1
c     Set up frequencies
      do 1 inu=1,68
      phots(inu)=0.0
      nu(inu)=10.0**(10.0+0.25*(inu-1))
      if(inu.eq.2)nu(inu)=1.5e10
      if(inu.eq.3)nu(inu)=4.3e10
      if(inu.eq.4)nu(inu)=8.6e10
      if(inu.eq.6)nu(inu)=2.3e11
      if(inu.eq.19)nu(inu)=3.72e14
      if(inu.eq.20)nu(inu)=4.56e14
      if(inu.eq.21)nu(inu)=6.74e14
      if(inu.eq.22)nu(inu)=8.22e14
      if(inu.eq.23)nu(inu)=1.7783e15
      if(inu.eq.24)nu(inu)=3.1623e15
      if(inu.eq.25)nu(inu)=5.6234e15
      if(inu.eq.26)nu(inu)=1.0000e16
      if(inu.eq.28)nu(inu)=4.8e16
      if(inu.eq.29)nu(inu)=1.21e17
      if(inu.eq.30)nu(inu)=2.42e17
      if(inu.eq.31)nu(inu)=4.84e17
      if(inu.eq.32)nu(inu)=9.68e17
      if(inu.eq.33)nu(inu)=1.936e18
      if(inu.eq.34)nu(inu)=2.42e18
      if(inu.eq.35)nu(inu)=4.84e18
    1 snu(inu)=nu(inu)
cccc  Choose frequencies of radio map and optical maps
      nrmap=6
      nopmap=20
c     alpha = (s-1)/2 is the underlying spectral index, where s = slope of 
c       electron E dist.
c     gmaxmn,gmaxmx are the min/max values of initial gamma_max of
c       electrons in a cell
      gmaxmx=gmrat*gmaxmn
c     gmin is the minimum value of initial gamma of electrons
c     2p is the slope of the volume vs. initial gamma-max law over all cells,
c          V = V0*gamma_max**(-2p)
      pexp=-2.0
c     Set up compilation of distribution of initial gamma_max values
      gmrf=gmaxmx/gmaxmn
      gmrfl=alog10(gmrf)/43.0
      do 2 ie=1,44
      gfl=alog10(gmaxmn)+gmrfl*(ie-1)
      gcnt(ie)=10.0**gfl
      igcnt(ie)=0
    2 continue
      yr=3.16e7
      rad=180.0/3.14159
      pi=3.14159
      pio2=pi*0.5
      sq2=1.41421
      sq3=1.73205
      parsec=3.086e18
      c=3.0e10
      emc2=8.186e-7
      cc2=1.29e-9
c     fgeom is multiplier to take into account that grid of cells occupies less
c       than a full cross-section of radius (2*nrad-1)*rcell
      fgeom=1.33
      emfact=3.74e-23
      c6gam=8.2e-21
      amjy=1.0e-26
      zred1=1.0+zred
      sen=2.0*alpha+1.0
      zeta=zeta/rad
      opang=opang/rad
c     brand is the turbulent component of the B field, bord is the ordered component
      brand=bave*(1.0-bratio)
      bord=bave*bratio
      bpsin=sin(bpitch/rad)
      bpcos=cos(bpitch/rad)
c     uratio is the user-specified ratio of energy density of electrons 
c        to that of mag. field
c     Set the normalization of the electron energy distribution accordingly:
      n0ave=uratio*bave**2*(sen-2.0)/(8.0*pi*emc2)/
     , (gmin**(2.0-sen)-gmaxmn**(2.0-sen))
c     Line of sight
      thlos=thlos/rad
      slos=dsin(thlos)
      clos=dcos(thlos)
      sinz=dsin(zeta)
      cosz=dcos(zeta)
      tanz=sinz/cosz
      gamup=1.0d0/dsqrt(1.0d0-betaup**2)
c     Compression ratio of Mach disk
c       Ultra-relativistic eq. of state assumed, so compression ratio is that
c       given by Hughes, Aller, & Aller (1989, ApJ, 341, 54)
      etac=dsqrt(8.0d0*gamup**4-1.70d1*gamup**2+9.0d0)/gamup
c     Speed downstream of conical shock if turbulent velocity is ignored;
c        for setting cell length zsize and for first estimate of time delay
c      betadd=dsqrt((1.0d0-(betaup*cosz)**2)**2+9.0d0*
c     ,(betaup**2*cosz*sinz)**2)/(3.0*betaup*sinz)
      sx=-sinz
      sy=0.0
      sz=cosz
      call vdcalc(0.0d0,0.0d0,betaup,sx,sy,sz,anx,any,anz,
     ,  betadd,gammdd,eta)
      if(betadd.gt.0.57735.and.betadd.lt.betaup)go to 7891
      write(6,9891)betaup,betadd,sinz,cosz
      go to 9000
 7891 betd=betadd
      gamd=gammdd
c     Length of a cylindrical cell in pc
c     Set so that there are 8 cells parallel to axis from the start of one column
c       to the start of the next column
      zsize=zrat*rsize/tanz
c     zss is used as a diagnostic to make sure that zsize does not change because of memory problems
      zss=1.001*zsize
      volc=pi*rsize**2*zsize
c     Length and volume of cell in plasma proper frame
      zsizep=zsize/gamd
      volcp=volc/gamd
      svmd=sqrt(vmd)
      delobs=1.0d0/(gamd*(1.0d0-betd*clos))
c     Time step in observer's frame in days
      dtfact=(1.0d0-betd*clos)/betd
      dtime=1190.0*zsize*dtfact*zred1
c      ip0=rand(0)*5000
      ip0=100
c     Distance of shock from axis and apex of conical jet
      tanop=dtan(opang)
c     Next line is specific to the selected number of cells per slice
      rshock=2*(nrad-1)*rsize
c     Distance of Mach disk from z value where conical shock intersects radially outermost cells
      zshock=rshock/tanz
c     Maximum number of cells between shock and rarefaction; applies to outermost cells
      icells=2.0*zshock/zsize+0.5
      kmd=zshock/zsize+1.1
      dstart=mdmd+icells/betd
      write(5,9494)kmd
 9494 format('# kmd: ',i5)
c     Distance of Mach disk from from vertex of jet cone
      zsvtex=(rshock+rsize)/tanop+zshock
c      write(3,8888)iseed,alpha
      expon=alpha+1.0
      exp1=0.5*expon-1.0
      anorm=(1.0+alpha)/(alpha+5.0/3.0)
c     Computation of IR seed photon density from dust torus as function of distance down jet
c     Area filling factor of hot dust, from IR luminosity. 1.05E8 = 1E45/(parsec**2)
      filld=1.05e8*ldust/(5.67e-5*2.0*pi**2*dtrad*dtdist*tdust**4)
      icend=icells+1
c      icend=icells+nrad/zrat
      do 333 id=1,icend
      zdist=zdist0+(id-1)*zsize
c      zdist=zdist0+(id-nrad/zrat)*zsize+zshock
c     Calculate min & max angles of dust torus in plasma frame
      dphi1=asin(dtrad/sqrt(zdist**2+dtdist**2))
      dphi2=atan(dtdist/zdist)
      dth1=dphi2-dphi1
      dth2=dphi2+dphi1
      csth1=dcos(dth1)
      csth2=dcos(dth2)
      dcsth1=-(csth1-betadd)/(1.0d0-csth1*betadd)
      dcsth2=-(csth2-betadd)/(1.0d0-csth2*betadd)
c     Doppler factor of dust torus emission in frame of cells, used
c       to estimate frequency of peak intensity in plasma frame
      tdel=gammdd*(1.0d0-betadd*0.5*(csth1+csth2))
      tdelr(id)=tdel
      dsang=dth2-dth1
c     Calculate seed photon field from dust emission in plasma frame
c     Peak frequency of dust thermal emission for part of torus closest to shock
      seedpk(id)=5.88e10*tdust*tdel
c     Use this to set the frequency array of the seed photons
      do 3 inu=1,22
      dustnu(inu)=seedpk(id)*10**(-1.4+(inu-1)*0.1)
      dusti(inu)=filld*seedph(dustnu(inu))
      if(dusti(inu).lt.1.0e-20)dusti(inu)=1.0e-20
      dustf(id,inu)=dustnu(inu)
    3 dustii(id,inu)=dusti(inu)
c      do 3331 inu=1,22
  333 continue
      do 336 id=1,icend
      dflux=0.0
      do 6 inu=2,22
      if(dustii(id,inu-1).le.0.0.or.dustii(id,inu).le.0.0)go to 4
      a=alog10(dustii(id,inu)/dustii(id,inu-1))/
     ,  alog10(dustnu(inu)/dustnu(inu-1))
      dflux=dflux+dustii(id,inu-1)/(a+1.0)*dustnu(inu-1)*
     ,  ((dustnu(inu)/dustnu(inu-1))**(a+1.0)-1.0)
      go to 5
    4 dflux=dflux+0.5*(dustii(id,inu)+dustii(id,inu-1))*
     ,  (dustnu(inu)-dustnu(inu-1))
    5 continue
    6 continue
      useed(id)=4.0*pi*dflux/c
  336 continue
      gmratl=log10(gmaxmx/gmin)/40.0
      write(3,6665)nrad,itlast,dtime,zred,dgpc,alpha,pn,pb,bave,
     , uratio,bratio,bpitch,rsize,gmaxmn,
     , gmrat,gmin,betaup,betat,(zeta*rad),(thlos*rad),(opang*rad),
     , tdust,ldust,dtdist,dtrad,zdist0,useed(5),psdslp,vmd,filld
      write(4,6665)nrad,itlast,dtime,zred,dgpc,alpha,pn,pb,bave,
     , uratio,bratio,bpitch,rsize,gmaxmn,
     , gmrat,gmin,betaup,betat,(zeta*rad),(thlos*rad),(opang*rad),
     , tdust,ldust,dtdist,dtrad,zdist0,useed(5),psdslp,vmd,filld
      write(5,6665)nrad,itlast,dtime,zred,dgpc,alpha,pn,pb,bave,
     , uratio,bratio,bpitch,rsize,gmaxmn,
     , gmrat,gmin,betaup,betat,(zeta*rad),(thlos*rad),(opang*rad),
     , tdust,ldust,dtdist,dtrad,zdist0,useed(5),psdslp,vmd,filld
 6665 format('#Radius in no. of cells: ',i3,4x,'No. of time steps:',i6,
     ,   ' time step: ',f8.5/
     ,  '#redshift: ',f6.3,2x,'Distance in Gpc: ',f7.3/
     ,  '#spectral index: ',f4.2,2x,'standard dev./mean of density, ',
     ,  'field fluctuations: ',2f6.2/
     ,  '#mean unshocked magnetic field: ',f5.3,2x,'ratio of electron ',
     ,  'to mag. energy: ',e8.2/'#ratio of ordered to total B field: ',
     ,  f5.2,2x,'pitch angle of helical field: ',f5.2/
     ,  '#cell radius (pc): ',f6.3/
     ,  '#Min. value of gamma_max: ',f8.1,2x,'ratio of max. to min. ',
     ,  'values of gamma_max: ',f5.1/'#gamma_min: ',f6.1,/
     ,  '#upstream laminar velocity: ',f9.5,'c',2x,
     ,  'upstream turbulent velocity: ',f9.5,'c'/'#shock angle: ',f6.3,
     ,  2x,'viewing angle: ',f6.3,2x,'opening angle: ',f6.3/
     ,  '#Dust temperature: ',f6.0,2x,'dust luminosity ',f5.2,
     ,  'x10**45 erg/s'/'#distance of center of dust ',
     ,  'torus from black hole: ',f3.1,' pc'/'#radius of torus: ',
     ,  f3.1,' pc   ','Distance of shock from central ',
     ,  'engine: ',f5.2,' pc'/'#Energy density of seed photons in ',
     ,  'plasma frame: ',e9.2/'#-Slope of PSD: ',f5.1,
     ,  '     Area of Mach disk relative to other cell',1pe9.2/
     ,  '# Area filling factor of dust emission: ',e9.2)
c     Set up variation of energy density according to PSD; see Done et al.,
c       1992, ApJ, 400, 138, eq. B1
      tinc=dtime
      open (10,iostat=ios, err=9000, file='simdat.txt',
     , status='new')
      call psdsim(ndim,-psdslp,-psdslp,1.0,tinc,spsdx)
      spxmn=0.0
      psdmn=0.0
      spmin=1.0
      spmax=1.01
      spsig=0.0
      spxsig=0.0
      ipulse=dstart+100+ip0
      spexp=1.0/(0.5*expon+1.0)
      do 4995 ip=1,ndim
c     Calculate mean of spsdx
      spxmn=spxmn+spsdx(ip)/andim
 4995 continue
      do 4996 ip=1,ndim
c     Normalize spsd by standard deviation
      spxsig=spxsig+(spsdx(ip)-spxmn)**2/(andim-1.0)
 4996 continue
      spxsig=sqrt(spxsig)
      do 4997 ip=1,ndim
c     Normalize spsd by standard deviation and take exponential of result
c       to get amplitude of flux variation
      spsd1(ip)=exp(amppsd*spsdx(ip)/spxsig)
c     Need to scale n0 and B by a different factor to get the desired amplitude
c     -- But in this version of temz, B is not varied with spsd(ip), so ignore
c      spsd(ip)=spsd(ip)**(spexp)
c     Average of 8 time steps to smooth variations so that discreteness
c       of columns of cells does not cause artificial spikes of flux
      spsd(ip)=spsd1(ip)
      if(ip.gt.7)spsd(ip)=0.125*(spsd1(ip)+spsd1(ip-1)+spsd1(ip-2)+
     ,  spsd1(ip-3)+spsd1(ip-4)+spsd1(ip-5)+spsd1(ip-6)+spsd1(ip-7))
      psdmn=psdmn+spsd(ip)/andim
 4997 continue
      do 4998 ip=1,ndim
      spsig=spsig+(spsd(ip)-psdmn)**2/(andim-1.0)
      if(spsd(ip).lt.spmin)spmin=spsd(ip)
      if(spsd(ip).gt.spmax)spmax=spsd(ip)
 4998 continue
      spsig=sqrt(spsig)
      write(6,9992)nrad,ip0,dstart,zsize,zshock,clos,slos,
     ,   gammdd,delobs,dtfact,spxsig,spxmn,spsig,psdmn,spmin
c      normalize spsd by its minimum value
      spnorm=sqrt(spmin*spmax)
      do 4999 ip=1,ndim
      spsd(ip)=spsd(ip)/spnorm
c     If no modulation at input is desired, set spsd = 1
      if(modul.eq.0)spsd(ip)=1.0
      write(10,9345)ip,spsd(ip),spsdx(ip)
 4999 continue
 1000 continue
 9345 format(i10,1p2e14.5)
      close(10)
c     Set parameters of each cell at initial time
c     There are jcells columns of cells, with icells cells per column
c     The central cell (j=1) is a Mach disk with low velocity; its observed radiation
c      is ignored, but it is an  important source of IC seed photons
      write(4,6667)
      write(5,6668)
      i=1
      j=1
      jzero=2
      do 999 nnn=1,(nrad-1)
      do 998 jcnt=1,(6*nnn)
      j=j+1
      annn=nnn
      layer(j)=nrad-nnn
      angc=60.0/(rad*annn)
      rcell(j)=2.0*rsize*annn
      xcell(j)=rcell(j)*cos(angc*(j-jzero))
      ycell(j)=rcell(j)*sin(angc*(j-jzero))
      cosph(j)=1.0
      sinph(j)=0.0
      zcol=rcell(j)/tanz
      imax(j)=2.0*zcol/zsize+0.01
      if(imax(j).lt.2)imax(j)=2
c     nouter(j) = approx. no. of cells between cell of interest and observer
      nouter(j)=imax(j)
      cosph(j)=xcell(j)/rcell(j)
      sinph(j)=ycell(j)/rcell(j)
      tanpsi(j)=rcell(j)/(zsvtex-zcol)
      psi(j)=datan(tanpsi(j))
      cospsi(j)=dcos(psi(j))
      sinpsi(j)=tanpsi(j)*cospsi(j)
c     Set up turbulent zone numbers
      ng1x(j)=1+2*(xcell(j)/rsize+2*nrad)/nrad
      ng1y(j)=1+2*(ycell(j)/rsize+2*nrad)/nrad
      if(ng1y(j).eq.5.and.ng1x(j).eq.4)ng1y(j)=4
      if(ng1x(j).lt.1)ng1x(j)=1
      if(ng1y(j).lt.1)ng1y(j)=1
      if(ng1x(j).gt.8)ng1x(j)=8
      if(ng1y(j).gt.8)ng1y(j)=8
      ng2x(j)=1+(xcell(j)/rsize+2*nrad)/nrad
      ng2y(j)=1+(ycell(j)/rsize+2*nrad)/nrad
      ranc=rand(0)
      if(ranc.lt.0.74)go to 555
      if(ng2y(j).eq.3.and.ng2x(j).eq.2)ng2y(j)=2
  555 continue
      if(ng2x(j).lt.1)ng2x(j)=1
      if(ng2y(j).lt.1)ng2y(j)=1
      if(ng2x(j).gt.4)ng2x(j)=4
      if(ng2y(j).gt.4)ng2y(j)=4
      ng3x(j)=1
      if(xcell(j).gt.0)ng3x(j)=2
      ng3y(j)=1
      if(ycell(j).gt.0)ng3y(j)=2
  998 continue
      jzero=j+1
  999 continue
      xcell(1)=0.0
      ycell(1)=0.0
      rcell(1)=0.0
      cosph(1)=1.0/sq2
      sinph(1)=1.0/sq2
      ng1x(1)=4
      ng1y(1)=5
      ng2x(1)=3
      ng2y(1)=2
      ng3x(1)=2
      ng3y(1)=1
      layer(1)=nrad
c     Determine indices that indicate which cells are first ones
c       encountered for a given section across the jet in later loop
      do 888 j=jcells,1,-1
      if(ididd1(ng1x(j),ng1y(j),layer(j)).eq.0)ididk1(j,layer(j))=1
      if(ididd2(ng2x(j),ng2y(j),layer(j)).eq.0)ididk2(j,layer(j))=1
      if(ididd3(ng3x(j),ng3y(j),layer(j)).eq.0)ididk3(j,layer(j))=1
      if(ididc1(ng1x(j),ng1y(j)).eq.0)ididj1(j)=1
      if(ididc2(ng2x(j),ng2y(j)).eq.0)ididj2(j)=1
      if(ididc3(ng3x(j),ng3y(j)).eq.0)ididj3(j)=1
      ididc1(ng1x(j),ng1y(j))=1
      ididc2(ng2x(j),ng2y(j))=1
      ididc3(ng3x(j),ng3y(j))=1
      ididd1(ng1x(j),ng1y(j),layer(j))=1
      ididd2(ng2x(j),ng2y(j),layer(j))=1
      ididd3(ng3x(j),ng3y(j),layer(j))=1
      do 887 i=1,imax(j)+1
      zcell(i,j)=zshock+(i-1)*zsize-rcell(j)/tanz
  887 continue
  888 continue
      imax(1)=1
      zcell(1,1)=zshock
      zrf=zshock
      xrf=0.0
c
c     *** Set up Mach disk emission for earlier times ***
c
c     Compute time delay (no. of time steps) between when plasma passes Mach disk
c       and when it passes reference point of conical shock, in plasma frame
      zmd=zrf-zshock
      idelmd=zmd/zsize
c     Velocity parameters of plasma in Mach disk
      betadx(1,1)=0.0
      betady(1,1)=0.0
      betadz(1,1)=1.0d0/3.0d0
      betamd=betadz(1,1)
      gammd=1.0d0/dsqrt(1.0d0-betamd**2)
      betad(1,1)=betamd
      gammad(1,1)=gammd
      dopref=gammdd/gammd
c     Determine B vector of MD assuming random magnetic field orientation
c     Randomly select magnetic field direction for every nth cell, then interpolate
c       inside loop to get direction for intermediate cells, where n = 1 for
c       zone 0 (each cell separately), n = 2 for zone 1, n=4 for zone 2, and
c       n=8 for zone 3
c     First need to set up field direction for the md=2 cell
      i=1
      j=1
      phi21(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      phi22(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      phi23(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet21(ng1x(j),ng1y(j))=acos(costh)
      call randln(1.0,pb,rndb21)
      call randln(1.0,pn,rndn21)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet22(ng2x(j),ng2y(j))=acos(costh)
      call randln(1.0,pb,rndb22)
      call randln(1.0,pn,rndn22)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet23(ng3x(j),ng3y(j))=acos(costh)
      call randln(1.0,pb,rndb23)
      call randln(1.0,pn,rndn23)
ccccc     Loop over mdmax cells that pass through Mach disk
      tcpu1 = omp_get_wtime()
      do 130 md=1,(mdmax-1)
      phi10=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet10=acos(costh)
      sint1=sin(thet10)
      call randln(1.0,pb,rndb10)
      call randln(1.0,pn,rndn10)
      bux0=sint1*cos(phi10)*rndb10
      buy0=sint1*sin(phi10)*rndb10
      buz0=costh*rndb10
      an00=n0ave*rndn10
c
      if(ididb1(ng1x(j),ng1y(j)).eq.1)go to 6140
      phi11(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet11(ng1x(j),ng1y(j))=acos(costh)
      sint1=sin(thet11(ng1x(j),ng1y(j)))
      sint2=sin(thet21(ng1x(j),ng1y(j)))
      call randln(1.0,pb,rndb11)
      call randln(1.0,pn,rndn11)
      bu1x1(ng1x(j),ng1y(j))=sint1*cos(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1y1(ng1x(j),ng1y(j))=sint1*sin(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1z1(ng1x(j),ng1y(j))=costh*rndb11
      bu2x1(ng1x(j),ng1y(j))=sint2*cos(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2y1(ng1x(j),ng1y(j))=sint2*sin(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2z1(ng1x(j),ng1y(j))=cos(thet21(ng1x(j),ng1y(j)))*rndb21
      cosrt1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))+
     ,  bu1y1(ng1x(j),ng1y(j))*bu2y1(ng1x(j),ng1y(j))+
     ,  bu1z1(ng1x(j),ng1y(j))*bu2z1(ng1x(j),ng1y(j))
      cdenom=sqrt((bu1x1(ng1x(j),ng1y(j))**2+bu1y1(ng1x(j),ng1y(j))**2
     ,  +bu1z1(ng1x(j),ng1y(j))**2)*(bu2x1(ng1x(j),ng1y(j))**2+
     ,   bu2y1(ng1x(j),ng1y(j))**2+bu2z1(ng1x(j),ng1y(j))**2))
      cosrt1(ng1x(j),ng1y(j))=cosrt1(ng1x(j),ng1y(j))/cdenom
      if(cosrt1(ng1x(j),ng1y(j)).ge.1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=0.999999
      if(cosrt1(ng1x(j),ng1y(j)).le.-1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=-0.999999
      angrt1(ng1x(j),ng1y(j))=acos(cosrt1(ng1x(j),ng1y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt1(ng1x(j),ng1y(j))=
     ,  angrt1(ng1x(j),ng1y(j))-2.0*pi
      cu1x1(ng1x(j),ng1y(j))=bu1y1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))-bu1z1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))
      cu1y1(ng1x(j),ng1y(j))=-bu1x1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))+bu1z1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      cu1z1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))-bu1y1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      ididb1(ng1x(j),ng1y(j))=1
 6140 continue
      anggrr=bfrac1(ng1x(j),ng1y(j))*angrt1(ng1x(j),ng1y(j))
      call vecrot(bu1x1(ng1x(j),ng1y(j)),bu1y1(ng1x(j),ng1y(j)),
     ,  bu1z1(ng1x(j),ng1y(j)),cu1x1(ng1x(j),ng1y(j)),
     ,  cu1y1(ng1x(j),ng1y(j)),cu1z1(ng1x(j),ng1y(j)),anggrr,
     ,  bux1(ng1x(j),ng1y(j),kmd),buy1(ng1x(j),ng1y(j),kmd),
     ,  buz1(ng1x(j),ng1y(j),kmd))
      an01=n0ave*(rndn11*bfrac1(ng1x(j),ng1y(j))+
     ,         rndn21*(1.0-bfrac1(ng1x(j),ng1y(j))))
      bfrac1(ng1x(j),ng1y(j))=bfrac1(ng1x(j),ng1y(j))+0.9999
      if(bfrac1(ng1x(j),ng1y(j)).lt.1.0)go to 6141
      thet21(ng1x(j),ng1y(j))=thet11(ng1x(j),ng1y(j))
      phi21(ng1x(j),ng1y(j))=phi11(ng1x(j),ng1y(j))
      bfrac1(ng1x(j),ng1y(j))=0.0
      ididb1(ng1x(j),ng1y(j))=0
 6141 continue
      if(ididb2(ng2x(j),ng2y(j)).eq.1)go to 6150
      phi12(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet12(ng2x(j),ng2y(j))=acos(costh)
      sint1=sin(thet12(ng2x(j),ng2y(j)))
      sint2=sin(thet22(ng2x(j),ng2y(j)))
      call randln(1.0,pb,rndb12)
      call randln(1.0,pn,rndn12)
      bu1x2(ng2x(j),ng2y(j))=sint1*cos(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1y2(ng2x(j),ng2y(j))=sint1*sin(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1z2(ng2x(j),ng2y(j))=costh*rndb12
      bu2x2(ng2x(j),ng2y(j))=sint2*cos(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2y2(ng2x(j),ng2y(j))=sint2*sin(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2z2(ng2x(j),ng2y(j))=cos(thet22(ng2x(j),ng2y(j)))*rndb22
      cosrt2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))+
     ,  bu1y2(ng2x(j),ng2y(j))*bu2y2(ng2x(j),ng2y(j))+
     ,  bu1z2(ng2x(j),ng2y(j))*bu2z2(ng2x(j),ng2y(j))
      cdenom=sqrt((bu1x2(ng2x(j),ng2y(j))**2+bu1y2(ng2x(j),ng2y(j))**2
     ,  +bu1z2(ng2x(j),ng2y(j))**2)*(bu2x2(ng2x(j),ng2y(j))**2+
     ,   bu2y2(ng2x(j),ng2y(j))**2+bu2z2(ng2x(j),ng2y(j))**2))
      cosrt2(ng2x(j),ng2y(j))=cosrt2(ng2x(j),ng2y(j))/cdenom
      if(cosrt2(ng2x(j),ng2y(j)).ge.1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=0.999999
      if(cosrt2(ng2x(j),ng2y(j)).le.-1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=-0.999999
      angrt2(ng2x(j),ng2y(j))=acos(cosrt2(ng2x(j),ng2y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt2(ng2x(j),ng2y(j))=
     ,  angrt2(ng2x(j),ng2y(j))-2.0*pi
      cu1x2(ng2x(j),ng2y(j))=bu1y2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))-bu1z2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))
      cu1y2(ng2x(j),ng2y(j))=-bu1x2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))+bu1z2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      cu1z2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))-bu1y2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      ididb2(ng2x(j),ng2y(j))=1
 6150 continue
      anggrr=bfrac2(ng2x(j),ng2y(j))*angrt2(ng2x(j),ng2y(j))
      call vecrot(bu1x2(ng2x(j),ng2y(j)),bu1y2(ng2x(j),ng2y(j)),
     ,  bu1z2(ng2x(j),ng2y(j)),cu1x2(ng2x(j),ng2y(j)),
     ,  cu1y2(ng2x(j),ng2y(j)),cu1z2(ng2x(j),ng2y(j)),anggrr,
     ,  bux2(ng2x(j),ng2y(j),kmd),buy2(ng2x(j),ng2y(j),kmd),
     ,  buz2(ng2x(j),ng2y(j),kmd))
      an02=n0ave*(rndn12*bfrac2(ng2x(j),ng2y(j))+
     ,         rndn22*(1.0-bfrac2(ng2x(j),ng2y(j))))
c      if(md.ge.10000.and.md.le.10020)
c     , write(5,9650)bu1x2(ng2x(j),ng2y(j)),bu2x2(ng2x(j),ng2y(j)),
c     , bux2(ng2x(j),ng2y(j),kmd),bu1y2(ng2x(j),ng2y(j)),
c     , bu2y2(ng2x(j),ng2y(j)),buy2(ng2x(j),ng2y(j),kmd),
c     , bu1z2(ng2x(j),ng2y(j)),bu2z2(ng2x(j),ng2y(j)),
c     , buz2(ng2x(j),ng2y(j),kmd),bfrac2(ng2x(j),ng2y(j)),
c     , angrt2(ng2x(j),ng2y(j))
c 9650 format('zone 2 ',1p11e11.4)
      bfrac2(ng2x(j),ng2y(j))=bfrac2(ng2x(j),ng2y(j))+0.33333
      if(bfrac2(ng2x(j),ng2y(j)).lt.1.0)go to 6151
      thet22(ng2x(j),ng2y(j))=thet12(ng2x(j),ng2y(j))
      phi22(ng2x(j),ng2y(j))=phi12(ng2x(j),ng2y(j))
      bfrac2(ng2x(j),ng2y(j))=0.0
      ididb2(ng2x(j),ng2y(j))=0
 6151 continue
c
      if(ididb3(ng3x(j),ng3y(j)).eq.1)go to 6160
      phi13(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet13(ng3x(j),ng3y(j))=acos(costh)
      sint1=sin(thet13(ng3x(j),ng3y(j)))
      sint2=sin(thet23(ng3x(j),ng3y(j)))
      call randln(1.0,pb,rndb13)
      call randln(1.0,pn,rndn13)
      bu1x3(ng3x(j),ng3y(j))=sint1*cos(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1y3(ng3x(j),ng3y(j))=sint1*sin(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1z3(ng3x(j),ng3y(j))=costh*rndb13
      bu2x3(ng3x(j),ng3y(j))=sint2*cos(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2y3(ng3x(j),ng3y(j))=sint2*sin(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2z3(ng3x(j),ng3y(j))=cos(thet23(ng3x(j),ng3y(j)))*rndb23
      cosrt3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))+
     ,  bu1y3(ng3x(j),ng3y(j))*bu2y3(ng3x(j),ng3y(j))+
     ,  bu1z3(ng3x(j),ng3y(j))*bu2z3(ng3x(j),ng3y(j))
      cdenom=sqrt((bu1x3(ng3x(j),ng3y(j))**2+bu1y3(ng3x(j),ng3y(j))**2
     ,  +bu1z3(ng3x(j),ng3y(j))**2)*(bu2x3(ng3x(j),ng3y(j))**2+
     ,   bu2y3(ng3x(j),ng3y(j))**2+bu2z3(ng3x(j),ng3y(j))**2))
      cosrt3(ng3x(j),ng3y(j))=cosrt3(ng3x(j),ng3y(j))/cdenom
      if(cosrt3(ng3x(j),ng3y(j)).ge.1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=0.999999
      if(cosrt3(ng3x(j),ng3y(j)).le.-1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=-0.999999
      angrt3(ng3x(j),ng3y(j))=acos(cosrt3(ng3x(j),ng3y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt3(ng3x(j),ng3y(j))=
     ,  angrt3(ng3x(j),ng3y(j))-2.0*pi
      cu1x3(ng3x(j),ng3y(j))=bu1y3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))-bu1z3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))
      cu1y3(ng3x(j),ng3y(j))=-bu1x3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))+bu1z3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      cu1z3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))-bu1y3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      ididb3(ng3x(j),ng3y(j))=1
 6160 continue
      anggrr=bfrac3(ng3x(j),ng3y(j))*angrt3(ng3x(j),ng3y(j))
      call vecrot(bu1x3(ng3x(j),ng3y(j)),bu1y3(ng3x(j),ng3y(j)),
     ,  bu1z3(ng3x(j),ng3y(j)),cu1x3(ng3x(j),ng3y(j)),
     ,  cu1y3(ng3x(j),ng3y(j)),cu1z3(ng3x(j),ng3y(j)),anggrr,
     ,  bux3(ng3x(j),ng3y(j),kmd),buy3(ng3x(j),ng3y(j),kmd),
     ,  buz3(ng3x(j),ng3y(j),kmd))
      an03=n0ave*(rndn13*bfrac3(ng3x(j),ng3y(j))+
     ,         rndn23*(1.0-bfrac3(ng3x(j),ng3y(j))))
      bfrac3(ng3x(j),ng3y(j))=bfrac3(ng3x(j),ng3y(j))+0.14285
      if(bfrac3(ng3x(j),ng3y(j)).lt.1.0)go to 6161
      thet23(ng3x(j),ng3y(j))=thet13(ng3x(j),ng3y(j))
      phi23(ng3x(j),ng3y(j))=phi13(ng3x(j),ng3y(j))
      bfrac3(ng3x(j),ng3y(j))=0.0
      ididb3(ng3x(j),ng3y(j))=0
 6161 continue
      bux=bnorm0*bux0+bnorm1*bux1(ng1x(j),ng1y(j),kmd)
     ,  +bnorm2*bux2(ng2x(j),ng2y(j),kmd)
     ,  +bnorm3*bux3(ng3x(j),ng3y(j),kmd)
      buy=bnorm0*buy0+bnorm1*buy1(ng1x(j),ng1y(j),kmd)
     ,  +bnorm2*buy2(ng2x(j),ng2y(j),kmd)
     ,  +bnorm3*buy3(ng3x(j),ng3y(j),kmd)
      buz=bnorm0*buz0+bnorm1*buz1(ng1x(j),ng1y(j),kmd)
     ,  +bnorm2*buz2(ng2x(j),ng2y(j),kmd)
     ,  +bnorm3*buz3(ng3x(j),ng3y(j),kmd)
      bcmd(md,1)=bux
      bcmd(md,2)=buy
      bcmd(md,3)=buz
c     Compute B field components downstream of Mach disk shock
      idelay=dstart-mdmd+md-idelmd
      if(idelay.lt.1.or.idelay.gt.(ndim-ip0))
     ,  write(6,9299)idelay,i,jcells,md,idelmd,ip0,mdmax,dstart,
     ,  zshock,delobs,gamd,betd,clos
      n0mean=spsd(idelay+ip0)*(anorm0*an00+anorm1*an01+anorm2*an02+
     ,  anorm3*an03)
      n0(i,j)=etac*n0mean
      bx(i,j)=(brand*bux-bord*bpsin*sinph(j))*etac
      by(i,j)=(brand*buy+bord*bpsin*cosph(j))*etac
      bz(i,j)=brand*buz+bord*bpcos
      bfield(j)=sqrt(bx(i,j)**2+by(i,j)**2+bz(i,j)**2)
      bfld=bfield(j)
      bmdx(md)=bx(i,j)
      bmdy(md)=by(i,j)
      bmdz(md)=bz(i,j)
      bmdtot(md)=bfld
c     Calculate the initial maximum electron energy in the Mach disk
      ball2=bux*bux+buy*buy+buz*buz
c     Next line relates this to direction of upstream B field relative to shock
c      gmax0(i,j)=gmaxmx*buz*buz/ball2
c      if(gmax0(i,j).lt.gmaxmn)gmax0(i,j)=gmaxmn
c     Next 3 lines assume a power-law distribution of gmax0, unrelated
c         to direction of B field
c     See http://mathworld.wolfram.com/RandomNumber.html for choosing
c      random numbers from a power-law distribution
c      xrand=rand(0)
c      if(pexp.eq.0.0)gmax0(i,j)=gmaxmx
c      if(pexp.lt.0.0)gmax0(i,j)=((gmaxmx**pexp-gmaxmn**pexp)*xrand+
c     ,  gmaxmn**pexp)**(1.0/pexp)
      gmax0(i,j)=gmaxmn
      gminmd=gmin
c     Calculate energy distribution in the Mach disk cell
c     Compute energy density of photons in Mach disk, time delayed by 1 step
c       Ignore photons from dust torus since beaming is small in Mach disk
      mdi=md-1
      if(mdi.le.0)mdi=1
c     Value of SSC photon density for fast cooling case, from Sari and Esen (2001)
      uphmd=bmdtot(mdi)**2/(8.0*pi)*(sqrt(1.0+4.0*uratio)-1.0)/2.0
      ustob=8.0*pi*uphmd/(bfield(j))**2
      tlfact=cc2*bfield(j)**2*(1.0+ustob)
c     Effective length of Mach disk zsmd = loss time of highest-E
c       electrons (to energy glow where radiation is negligible) times flow velocity
      glow=10.0
      delt=(gmax0(i,j)-glow)/(tlfact*glow*gmax0(i,j))
      zsmd(md)=delt/(yr*3.26)*gammd*betamd
c     Use 0.99 instead of 1 to avoid singularity of N(gamma) at gmax0
      gmrat=0.99*gmax0(i,j)/gminmd
      gmratl=dlog10(gmrat)/32.0
      gmratm=dlog10(gminmd/glow)/11.0
      do 123 ie=1,44
      if(ie.ge.12)egam(i,j,ie)=gminmd*10.0**(gmratl*(ie-12))
      if(ie.lt.12)egam(i,j,ie)=glow*10.0**(gmratm*(ie-1))
      ggam(ie)=egam(i,j,ie)
      t2=(gmax0(i,j)-ggam(ie))/(tlfact*ggam(ie)*gmax0(i,j))
      enofe(i,j,ie)=0.0
      eterm1=ggam(ie)/gminmd
      eterm2=1.0d0-ggam(ie)*t2*tlfact
      if(eterm2.lt.0.0)go to 122
      if(ie.ge.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(1.0d0-eterm2**(sen-1.0))
      if(ie.lt.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(eterm1**(sen-1.0)-
     ,  eterm2**(sen-1.0))
c     Divide by delt since integral is over time
      delt=zsize*yr*3.26/(gammad(i,j)*betad(i,j))
      enofe(i,j,ie)=enofe(i,j,ie)/delt
      if(enofe(i,j,ie).lt.1.0e-35)enofe(i,j,ie)=0.0
  122 edist(ie)=enofe(i,j,ie)
  123 continue
      bperpp=bfld
c     Now estimate synchrotron emission seen by other cells
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 125 inu=1,40
      restnu=nu(inu)
      fsscmd(md,inu)=0.0
      sstau(inu)=0.0
c     Synchrotron mean intensity for SSC calculation inside Mach disk
      ssabs(inu)=0.0
      if(inu.lt.16)ssabs(inu)=1.02e4*(sen+2.0)*akapnu(restnu)*bperpp
      abssss=(parsec/nu(inu))/nu(inu)
      if(ssabs(inu).gt.(1.0e-35/abssss))sstau(inu)=ssabs(inu)*
     ,  abssss*rsize*svmd
      fsync(inu)=ajnu(nu(inu))*bperpp*rsize*svmd*(emfact*parsec)
      if(sstau(inu).lt.0.01)go to 124
      srcfn(inu)=fsync(inu)/sstau(inu)
      if(sstau(inu).gt.5.0)fsync(inu)=srcfn(inu)
      if(sstau(inu).le.5.0)fsync(inu)=
     ,  srcfn(inu)*(1.0-exp(-sstau(inu)))
  124 continue
      fsynmd(md,inu)=ajnu(nu(inu)/dopref)*dopref**2*bperpp
      absorb(md,inu)=0.0
      abssss=1.02e4*(sen+2.0)*akapnu(nu(inu)/dopref)*bperpp/nu(inu)
      if(abssss.gt.(1.0e-35*nu(inu)))absorb(md,inu)=(abssss/nu(inu))
      alphmd(md,inu)=10.0
      ssseed(inu)=fsync(inu)
c      if(md.eq.1.or.md.eq.428)
c     ,  write(4,9994)md,inu,restnu,fsync(inu),(restnu/dopref),
c     ,   fsynmd(md,inu),sstau(inu),
c     ,   absorb(md,inu),alphmd(md,inu),bperpp,dopref
  125 continue
!$omp end parallel do
c     Calculate spectral indices
c     Need to add a lower frequency to get spectral index of nu(1)
      fq1=0.98*nu(1)
c      fsyn1=ajnu(fq1/dopref)*dopref**2*bperpp/dtfact
      fsyn1=ajnu(fq1/dopref)*dopref**2*bperpp
      absrb1=0.0
      abssss=1.02e4*(sen+2.0)*akapnu(fq1/dopref)*bperpp
      if(abssss.gt.1.0e-18)absrb1=absss/(fq1/dopref)**2
      do 126 inu=1,40
      restnu=nu(inu)
      alfmdc(md,inu)=10.0
      abexmd(md,inu)=1.7
      if(fsynmd(md,inu).gt.0.0.and.fsyn1.gt.0.0)
     ,alphmd(md,inu)=-dlog10(fsynmd(md,inu)/fsyn1)/alog10(restnu/fq1)
      if(absorb(md,inu).gt.0.0.and.absrb1.gt.0.0)
     ,abexmd(md,inu)=-dlog10(absorb(md,inu)/absrb1)/alog10(restnu/fq1)
      fq1=restnu
      fsyn1=fsynmd(md,inu)
      absrb1=absorb(md,inu)
  126 continue
      do 127 inu=41,68
      alphmd(md,inu)=10.0
      ssseed(inu)=0.0
      fsync(inu)=0.0
      fsynmd(md,inu)=0.0
      absorb(md,inu)=0.0
  127 continue
c     Calculate the SSC emission from the Mach disk for reference
c       relative Doppler factor dopref
      betd=betamd
      gamd=gammd
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 128 inu=17,68
      fsscmd(md,inu)=ssc(nu(inu)/dopref)*dopref**2
  128 continue
!$omp end parallel do
c     Calculate the SSC spectral index
      fq1=0.98*nu(17)
      fssc1=ssc(fq1/dopref)*dopref**2
      do 129 inu=17,68
      restnu=nu(inu)
      alfmdc(md,inu)=10.0
      if(fsscmd(md,inu).gt.0.0.and.fssc1.gt.0.0)
     ,alfmdc(md,inu)=-alog10(fsscmd(md,inu)/fssc1)/alog10(restnu/fq1)
      fq1=restnu
      fssc1=fsscmd(md,inu)
  129 continue
      if((mdmax-md).gt.icells)go to 131
c     Carry field values of zones including Mach disk to later times
      k=kmd+icells
      bux1(ng1x(j),ng1y(j),k)=bux1(ng1x(j),ng1y(j),kmd)
      buy1(ng1x(j),ng1y(j),k)=buy1(ng1x(j),ng1y(j),kmd)
      buz1(ng1x(j),ng1y(j),k)=buz1(ng1x(j),ng1y(j),kmd)
      bux2(ng2x(j),ng2y(j),k)=bux2(ng2x(j),ng2y(j),kmd)
      buy2(ng2x(j),ng2y(j),k)=buy2(ng2x(j),ng2y(j),kmd)
      buz2(ng2x(j),ng2y(j),k)=buz2(ng2x(j),ng2y(j),kmd)
      bux3(ng3x(j),ng3y(j),k)=bux3(ng3x(j),ng3y(j),kmd)
      buy3(ng3x(j),ng3y(j),k)=buy3(ng3x(j),ng3y(j),kmd)
      buz3(ng3x(j),ng3y(j),k)=buz3(ng3x(j),ng3y(j),kmd)
  131 continue
  130 continue
c
c     *** End Mach disk set-up ***
c
c     After initial set-up, time loop starts here
    9 i=1
      ncells=0
      it=it+1
      itt=it-icells
      md=mdmax
      ididj3(jcells)=1
      ididj2(1)=0
      ididj3(1)=0
c
c     Start loop over all cells in first layer to set up physical parameters
c
      do 80 j=jcells,1,-1
c
      do 81 inu=1,68
      fsynmd(md,inu)=0.0
      fsscmd(md,inu)=0.0
      fmdall(inu)=0.0
   81 continue
      k=kmd-(zshock-zcell(i,j))/zsize
      if(k.lt.1)k=1
      idelay=dstart+it-0.5+((zrf-zcell(i,j))*clos
     ,   +(xrf-xcell(j))*slos)*betadd/zsize
      if(idelay.le.0.or.idelay.gt.(ndim-ip0))
     ,  write(5,9222)idelay,dstart,j,i,md,zshock,zrf,zcell(i,j),
     ,  xrf,xcell(j),betadd,zsize
      phcell(j)=datan2(sinph(j),cosph(j))
c     Determine v_turb & B vectors of cell (ignore v_turb for Mach disk)
c     Assume that v_turb follows Kolmogorov spectrum. Note: if v_turb >~ 0.5 should use
c       4-velocity version given in eq. 3 of Zrake & MacFadyen (2013, ApJL, 763, L12).
c     Assume random B orientations in each of 4 zones of different sizes
c     In the systemic fluid frame, each zone rotates at a turbulence velocity about a
c       randomly directed axis; the turbulent component of the velocity decreases with
c       distance from the axis
c
c     For Mach disk, continue from previous calculation of B direction
c     Interpolate inside loop to get direction for cells not on
c       a section boundary
c     First need to set up field direction for the more downstream sections
      if(ididk1(j,layer(j)).eq.0.or.ididb1(ng1x(j),ng1y(j)).eq.1)
     ,  go to 6035
      phi21(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet21(ng1x(j),ng1y(j))=acos(costh)
      call randln(1.0,pb,rndb21)
      call randln(1.0,pn,rndn21)
 6035 if(ididk2(j,layer(j)).eq.0.or.ididb2(ng2x(j),ng2y(j)).eq.1)
     ,  go to 6036
      phi22(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet22(ng2x(j),ng2y(j))=acos(costh)
      call randln(1.0,pb,rndb22)
      call randln(1.0,pn,rndn22)
 6036 if(ididk3(j,layer(j)).eq.0.or.ididb3(ng3x(j),ng3y(j)).eq.1)
     ,  go to 6037
      phi23(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet23(ng3x(j),ng3y(j))=acos(costh)
      call randln(1.0,pb,rndb23)
      call randln(1.0,pn,rndn23)
ccccc   Determine field component for each section of cells
c     Start with single cell-sized section 0
 6037 phi10=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet10=acos(costh)
      sint1=sin(thet10)
      phv10=2.0*pi*rand(0)
      cosv=2.0*(rand(0)-0.5)
      if(cosv.ge.1.0)cosv=0.999999
      if(cosv.le.-1.0)cosv=-0.999999
      thv10=acos(cosv)
      sinv=sin(thv10)
      vtz0=cosv
      vtx0=sinv*cos(phv10)
      vty0=sinv*sin(phv10)
      call randln(1.0,pb,rndb10)
      call randln(1.0,pn,rndn10)
      bux0=sint1*cos(phi10)*rndb10
      buy0=sint1*sin(phi10)*rndb10
      buz0=costh*rndb10
      an00=n0ave*rndn10
c     Next is section 1, equal in size to 2x2=4 cells
      if(ididj1(j).eq.1)go to 6038
      if(ididk1(j,layer(j)).eq.0) go to 6041
 6038 if(ididl1(ng1x(j),ng1y(j),layer(j)).eq.1)go to 6040
      phi11(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet11(ng1x(j),ng1y(j))=acos(costh)
      sint1=sin(thet11(ng1x(j),ng1y(j)))      
      sint2=sin(thet21(ng1x(j),ng1y(j)))
      phv11(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      cosv11(ng1x(j),ng1y(j))=2.0*(rand(0)-0.5)
      if(cosv11(ng1x(j),ng1y(j)).ge.1.0)cosv11(ng1x(j),ng1y(j))=0.999999
      if(cosv11(ng1x(j),ng1y(j)).le.-1.0)cosv11(ng1x(j),ng1y(j))=
     ,  -0.999999
      thv=acos(cosv11(ng1x(j),ng1y(j)))
      sinv11(ng1x(j),ng1y(j))=sin(thv)
      call randln(1.0,pb,rndb11)
      call randln(1.0,pn,rndn11)
      bu1x1(ng1x(j),ng1y(j))=sint1*cos(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1y1(ng1x(j),ng1y(j))=sint1*sin(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1z1(ng1x(j),ng1y(j))=costh*rndb11
      bu2x1(ng1x(j),ng1y(j))=sint2*cos(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2y1(ng1x(j),ng1y(j))=sint2*sin(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2z1(ng1x(j),ng1y(j))=cos(thet21(ng1x(j),ng1y(j)))*rndb21
      cosrt1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))+
     ,  bu1y1(ng1x(j),ng1y(j))*bu2y1(ng1x(j),ng1y(j))+
     ,  bu1z1(ng1x(j),ng1y(j))*bu2z1(ng1x(j),ng1y(j))
      cdenom=sqrt((bu1x1(ng1x(j),ng1y(j))**2+bu1y1(ng1x(j),ng1y(j))**2
     ,  +bu1z1(ng1x(j),ng1y(j))**2)*(bu2x1(ng1x(j),ng1y(j))**2+
     ,   bu2y1(ng1x(j),ng1y(j))**2+bu2z1(ng1x(j),ng1y(j))**2))
      cosrt1(ng1x(j),ng1y(j))=cosrt1(ng1x(j),ng1y(j))/cdenom
      if(cosrt1(ng1x(j),ng1y(j)).ge.1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=0.999999
      if(cosrt1(ng1x(j),ng1y(j)).le.-1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=-0.999999
      angrt1(ng1x(j),ng1y(j))=acos(cosrt1(ng1x(j),ng1y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt1(ng1x(j),ng1y(j))=
     ,  angrt1(ng1x(j),ng1y(j))-2.0*pi
      cu1x1(ng1x(j),ng1y(j))=bu1y1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))-bu1z1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))
      cu1y1(ng1x(j),ng1y(j))=-bu1x1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))+bu1z1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      cu1z1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))-bu1y1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      ididl1(ng1x(j),ng1y(j),layer(j))=1
 6040 continue
      anggrr=bfrac1(ng1x(j),ng1y(j))*angrt1(ng1x(j),ng1y(j))
      call vecrot(bu1x1(ng1x(j),ng1y(j)),bu1y1(ng1x(j),ng1y(j)),
     ,  bu1z1(ng1x(j),ng1y(j)),cu1x1(ng1x(j),ng1y(j)),
     ,  cu1y1(ng1x(j),ng1y(j)),cu1z1(ng1x(j),ng1y(j)),anggrr,
     ,  bux1(ng1x(j),ng1y(j),k),buy1(ng1x(j),ng1y(j),k),
     ,  buz1(ng1x(j),ng1y(j),k))
      an01=n0ave*(rndn11*bfrac1(ng1x(j),ng1y(j))+
     ,         rndn21*(1.0-bfrac1(ng1x(j),ng1y(j))))
      dvx=xcell(j)-(0.5*ng1x(j)-2.25)*nrad*rsize
      dvy=ycell(j)-(0.5*ng1y(j)-2.25)*nrad*rsize
      dvz=(0.5-bfrac1(ng1x(j),ng1y(j)))*zsize
c     If nrad=8, adjust placement of zone center so it is not at a cell center
      if(nrad.eq.8)then
        if(ng1x(j).lt.5)dvx=dvx+rsize
        if(ng1x(j).ge.5)dvx=dvx-rsize
        if(ng1y(j).lt.5)dvy=dvy+rsize
        if(ng1y(j).ge.5)dvy=dvy-rsize
        endif
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv11(ng1x(j),ng1y(j))
      axx=sinv11(ng1x(j),ng1y(j))*cos(phv11(ng1x(j),ng1y(j)))
      axy=sinv11(ng1x(j),ng1y(j))*sin(phv11(ng1x(j),ng1y(j)))
      vtx1=(axy*dvz-axz*dvy)/dvs
      vty1=(axz*dvx-axx*dvz)/dvs
      vtz1=(axx*dvy-axy*dvx)/dvs
      bfrac1(ng1x(j),ng1y(j))=bfrac1(ng1x(j),ng1y(j))+0.99999
      if(bfrac1(ng1x(j),ng1y(j)).lt.1.0)go to 6041
      thet21(ng1x(j),ng1y(j))=thet11(ng1x(j),ng1y(j))
      phi21(ng1x(j),ng1y(j))=phi11(ng1x(j),ng1y(j))
      bfrac1(ng1x(j),ng1y(j))=0.0
      ididl1(ng1x(j),ng1y(j),layer(j))=0
 6041 continue
c     Next is section 2, equal in size to 4x4=16 cells
      if(ididj2(j).eq.1)go to 6042
      if(ididk2(j,layer(j)).eq.0) go to 6051
 6042 if(ididl2(ng2x(j),ng2y(j),layer(j)).eq.1)go to 6050
      phi12(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet12(ng2x(j),ng2y(j))=acos(costh)
      sint1=sin(thet12(ng2x(j),ng2y(j)))
      sint2=sin(thet22(ng2x(j),ng2y(j)))
      phv12(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      cosv12(ng2x(j),ng2y(j))=2.0*(rand(0)-0.5)
      if(cosv12(ng2x(j),ng2y(j)).ge.1.0)cosv12(ng2x(j),ng2y(j))=0.999999
      if(cosv12(ng2x(j),ng2y(j)).le.-1.0)cosv12(ng2x(j),ng2y(j))=
     ,  -0.999999
      thv=acos(cosv12(ng2x(j),ng2y(j)))
      sinv12(ng2x(j),ng2y(j))=sin(thv)
      call randln(1.0,pb,rndb12)
      call randln(1.0,pn,rndn12)
      bu1x2(ng2x(j),ng2y(j))=sint1*cos(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1y2(ng2x(j),ng2y(j))=sint1*sin(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1z2(ng2x(j),ng2y(j))=costh*rndb12
      bu2x2(ng2x(j),ng2y(j))=sint2*cos(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2y2(ng2x(j),ng2y(j))=sint2*sin(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2z2(ng2x(j),ng2y(j))=cos(thet22(ng2x(j),ng2y(j)))*rndb22
      cosrt2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))+
     ,  bu1y2(ng2x(j),ng2y(j))*bu2y2(ng2x(j),ng2y(j))+
     ,  bu1z2(ng2x(j),ng2y(j))*bu2z2(ng2x(j),ng2y(j))
      cdenom=sqrt((bu1x2(ng2x(j),ng2y(j))**2+bu1y2(ng2x(j),ng2y(j))**2
     ,  +bu1z2(ng2x(j),ng2y(j))**2)*(bu2x2(ng2x(j),ng2y(j))**2+
     ,   bu2y2(ng2x(j),ng2y(j))**2+bu2z2(ng2x(j),ng2y(j))**2))
      cosrt2(ng2x(j),ng2y(j))=cosrt2(ng2x(j),ng2y(j))/cdenom
      if(cosrt2(ng2x(j),ng2y(j)).ge.1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=0.999999
      if(cosrt2(ng2x(j),ng2y(j)).le.-1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=-0.999999
      angrt2(ng2x(j),ng2y(j))=acos(cosrt2(ng2x(j),ng2y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt2(ng2x(j),ng2y(j))=
     ,  angrt2(ng2x(j),ng2y(j))-2.0*pi
      cu1x2(ng2x(j),ng2y(j))=bu1y2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))-bu1z2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))
      cu1y2(ng2x(j),ng2y(j))=-bu1x2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))+bu1z2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      cu1z2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))-bu1y2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      ididl2(ng2x(j),ng2y(j),layer(j))=1
 6050 continue
      anggrr=bfrac2(ng2x(j),ng2y(j))*angrt2(ng2x(j),ng2y(j))
      call vecrot(bu1x2(ng2x(j),ng2y(j)),bu1y2(ng2x(j),ng2y(j)),
     ,  bu1z2(ng2x(j),ng2y(j)),cu1x2(ng2x(j),ng2y(j)),
     ,  cu1y2(ng2x(j),ng2y(j)),cu1z2(ng2x(j),ng2y(j)),anggrr,
     ,  bux2(ng2x(j),ng2y(j),k),buy2(ng2x(j),ng2y(j),k),
     ,  buz2(ng2x(j),ng2y(j),k))
      an02=n0ave*(rndn12*bfrac2(ng2x(j),ng2y(j))+
     ,         rndn22*(1.0-bfrac2(ng2x(j),ng2y(j))))
      dvx=xcell(j)-(ng2x(j)-2.5)*nrad*rsize
      dvy=ycell(j)-(ng2y(j)-2.5)*nrad*rsize
      dvz=3.0*(0.5-bfrac2(ng2x(j),ng2y(j)))*zsize
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv12(ng2x(j),ng2y(j))
      axx=sinv12(ng2x(j),ng2y(j))*cos(phv12(ng2x(j),ng2y(j)))
      axy=sinv12(ng2x(j),ng2y(j))*sin(phv12(ng2x(j),ng2y(j)))
      vtx2=(axy*dvz-axz*dvy)/dvs
      vty2=(axz*dvx-axx*dvz)/dvs
      vtz2=(axx*dvy-axy*dvx)/dvs
      bfrac2(ng2x(j),ng2y(j))=bfrac2(ng2x(j),ng2y(j))+0.33333
      if(bfrac2(ng2x(j),ng2y(j)).lt.1.0)go to 6051
      thet22(ng2x(j),ng2y(j))=thet12(ng2x(j),ng2y(j))
      phi22(ng2x(j),ng2y(j))=phi12(ng2x(j),ng2y(j))
      bfrac2(ng2x(j),ng2y(j))=0.0
      ididl2(ng2x(j),ng2y(j),layer(j))=0
 6051 continue
c     Last is section 3, equal in size to 8x8=64 cells
      if(ididj3(j).eq.1)go to 6052
      if(ididk3(j,layer(j)).eq.0) go to 6061
 6052 if(ididl3(ng3x(j),ng3y(j),layer(j)).eq.1)go to 6060
      phi13(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet13(ng3x(j),ng3y(j))=acos(costh)
      sint1=sin(thet13(ng3x(j),ng3y(j)))
      sint2=sin(thet23(ng3x(j),ng3y(j)))
      phv13(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      cosv13(ng3x(j),ng3y(j))=2.0*(rand(0)-0.5)
      if(cosv13(ng3x(j),ng3y(j)).ge.1.0)cosv13(ng3x(j),ng3y(j))=0.999999
      if(cosv13(ng3x(j),ng3y(j)).le.-1.0)cosv13(ng3x(j),ng3y(j))=
     ,  -0.999999
      thv=acos(cosv13(ng3x(j),ng3y(j)))
      sinv13(ng3x(j),ng3y(j))=sin(thv)
      call randln(1.0,pb,rndb13)
      call randln(1.0,pn,rndn13)
      bu1x3(ng3x(j),ng3y(j))=sint1*cos(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1y3(ng3x(j),ng3y(j))=sint1*sin(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1z3(ng3x(j),ng3y(j))=costh*rndb13
      bu2x3(ng3x(j),ng3y(j))=sint2*cos(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2y3(ng3x(j),ng3y(j))=sint2*sin(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2z3(ng3x(j),ng3y(j))=cos(thet23(ng3x(j),ng3y(j)))*rndb23
      cosrt3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))+
     ,  bu1y3(ng3x(j),ng3y(j))*bu2y3(ng3x(j),ng3y(j))+
     ,  bu1z3(ng3x(j),ng3y(j))*bu2z3(ng3x(j),ng3y(j))
      cdenom=sqrt((bu1x3(ng3x(j),ng3y(j))**2+bu1y3(ng3x(j),ng3y(j))**2
     ,  +bu1z3(ng3x(j),ng3y(j))**2)*(bu2x3(ng3x(j),ng3y(j))**2+
     ,   bu2y3(ng3x(j),ng3y(j))**2+bu2z3(ng3x(j),ng3y(j))**2))
      cosrt3(ng3x(j),ng3y(j))=cosrt3(ng3x(j),ng3y(j))/cdenom
      if(cosrt3(ng3x(j),ng3y(j)).ge.1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=0.999999
      if(cosrt3(ng3x(j),ng3y(j)).le.-1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=-0.999999
      angrt3(ng3x(j),ng3y(j))=acos(cosrt3(ng3x(j),ng3y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt3(ng3x(j),ng3y(j))=
     ,  angrt3(ng3x(j),ng3y(j))-2.0*pi
      cu1x3(ng3x(j),ng3y(j))=bu1y3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))-bu1z3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))
      cu1y3(ng3x(j),ng3y(j))=-bu1x3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))+bu1z3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      cu1z3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))-bu1y3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      ididl3(ng3x(j),ng3y(j),layer(j))=1
 6060 continue
      anggrr=bfrac3(ng3x(j),ng3y(j))*angrt3(ng3x(j),ng3y(j))
      call vecrot(bu1x3(ng3x(j),ng3y(j)),bu1y3(ng3x(j),ng3y(j)),
     ,  bu1z3(ng3x(j),ng3y(j)),cu1x3(ng3x(j),ng3y(j)),
     ,  cu1y3(ng3x(j),ng3y(j)),cu1z3(ng3x(j),ng3y(j)),anggrr,
     ,  bux3(ng3x(j),ng3y(j),k),buy3(ng3x(j),ng3y(j),k),
     ,  buz3(ng3x(j),ng3y(j),k))
      an03=n0ave*(rndn13*bfrac3(ng3x(j),ng3y(j))+
     ,         rndn23*(1.0-bfrac3(ng3x(j),ng3y(j))))
      dvx=xcell(j)-(2.0*ng3x(j)-3.0)*sq2*(nrad-1)*rsize
      dvy=ycell(j)-(2.0*ng3y(j)-3.0)*sq2*(nrad-1)*rsize
      dvz=7.0*(0.5-bfrac3(ng3x(j),ng3y(j)))*zsize
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv13(ng3x(j),ng3y(j))
      axx=sinv13(ng3x(j),ng3y(j))*cos(phv13(ng3x(j),ng3y(j)))
      axy=sinv13(ng3x(j),ng3y(j))*sin(phv13(ng3x(j),ng3y(j)))
      vtx3=(axy*dvz-axz*dvy)/dvs
      vty3=(axz*dvx-axx*dvz)/dvs
      vtz3=(axx*dvy-axy*dvx)/dvs
      bfrac3(ng3x(j),ng3y(j))=bfrac3(ng3x(j),ng3y(j))+0.14285
      if(bfrac3(ng3x(j),ng3y(j)).lt.1.0)go to 6061
      thet23(ng3x(j),ng3y(j))=thet13(ng3x(j),ng3y(j))
      phi23(ng3x(j),ng3y(j))=phi13(ng3x(j),ng3y(j))
      bfrac3(ng3x(j),ng3y(j))=0.0
      ididl3(ng3x(j),ng3y(j),layer(j))=0
 6061 continue
c     Now add B components of all sections containing cell j
c     Weight according to section size a la scheme of Jones (1988)
      bux=bnorm0*bux0+bnorm1*bux1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*bux2(ng2x(j),ng2y(j),k)+bnorm3*bux3(ng3x(j),ng3y(j),k)
      buy=bnorm0*buy0+bnorm1*buy1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*buy2(ng2x(j),ng2y(j),k)+bnorm3*buy3(ng3x(j),ng3y(j),k)
      buz=bnorm0*buz0+bnorm1*buz1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*buz2(ng2x(j),ng2y(j),k)+bnorm3*buz3(ng3x(j),ng3y(j),k)
c     
c     Compute B field components downstream of shock in the plasma frame
c       by transforming the shock normal to the upstream plasma
c       frame, and compressing the component of B parallel to the shock
c
      bux=brand*bux-bord*bpsin*sinph(j)
      buy=brand*buy+bord*bpsin*cosph(j)
      buz=brand*buz+bord*bpcos
c     Unit vector of shock normal at current position
      anx=cosz*cosph(j)
      any=cosz*sinph(j)
      anz=sinz
c     Velocity vector of laminar component of pre-shock flow
      betupx=betaup*cosph(j)*sinpsi(j)
      betupy=betaup*sinph(j)*sinpsi(j)
      betupz=betaup*cospsi(j)
c     Velocity vector of the turbulent component of pre-shocked plasma
c     Add velocities of nested zones relativistically
      vvx=betat*vnorm0*vtx0
      vvy=betat*vnorm0*vty0
      vvz=betat*vnorm0*vtz0
      uux=betat*vnorm1*vtx1
      uuy=betat*vnorm1*vty1
      uuz=betat*vnorm1*vtz1
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      vvx1=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      vvy1=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      vvz1=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=1
c      write(3,9059)i,j,numv,vvx1,vvy1,vvz1,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
      vvx=vvx1
      vvy=vvy1
      vvz=vvz1
      uux=betat*vnorm2*vtx2
      uuy=betat*vnorm2*vty2
      uuz=betat*vnorm2*vtz2
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      vvx1=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      vvy1=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      vvz1=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=2
c      write(3,9059)i,j,numv,vvx1,vvy1,vvz1,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
      vvx=vvx1
      vvy=vvy1
      vvz=vvz1
      uux=betat*vnorm3*vtx3
      uuy=betat*vnorm3*vty3
      uuz=betat*vnorm3*vtz3
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      betatx=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      betaty=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      betatz=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=3
c      write(3,9059)i,j,numv,betatx,betaty,betatz,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
c     Velocity vector of the pre-shock plasma including turbulent component
      dotprd=betupx*betatx+betupy*betaty+betupz*betatz
      btparx=dotprd*betupx/betaup**2
      btpary=dotprd*betupy/betaup**2
      btparz=dotprd*betupz/betaup**2
      btprpx=betatx-btparx
      btprpy=betaty-btpary
      btprpz=betatz-btparz
      betaux(j)=(betupx+btparx+btprpx/gamup)/(1.0d0+dotprd)
      betauy(j)=(betupy+btpary+btprpy/gamup)/(1.0d0+dotprd)
      betauz(j)=(betupz+btparz+btprpz/gamup)/(1.0d0+dotprd)
      betau(j)=dsqrt(betaux(j)**2+betauy(j)**2+betauz(j)**2)
      gammau(j)=1.0d0/dsqrt(1.0d0-betau(j)**2)
c     Unit vector of shock front at current position
      sx=-sinz*cosph(j)
      sy=-sinz*sinph(j)
      sz=cosz
c     Velocity vector downstream of shock + compression ratio of shock
      if(j.gt.1)call vdcalc(betaux(j),betauy(j),betauz(j),sx,sy,sz,
     ,  betadx(i,j),betady(i,j),betadz(i,j),betad(i,j),gammad(i,j),eta)
c     Upstream velocity component parallel to shock normal
      velpar=betaux(j)*anx+betauy(j)*any+betauz(j)*anz
c     If Mach disk, compute separately
      if(j.eq.1)go to 12
c     Calculate upstream B field components parallel + perpendicular to shock normal
      call bcalc(betaux(j),betauy(j),betauz(j),betau(j),gammau(j),
     ,  anx,any,anz,bux,buy,buz,bparx,bpary,bparz,bprpx,bprpy,bprpz,
     ,  bpar,bprp)
      bx(i,j)=bparx+eta*bprpx
      by(i,j)=bpary+eta*bprpy
      bz(i,j)=bparz+eta*bprpz
      go to 13
c     Set N0 and field of plasma in central cell, which is a Mach disk
   12 n0mean=spsd(idelay+ip0)*(anorm0*an00+anorm1*an01+anorm2*an02+
     ,  anorm3*an03)
      n0(i,j)=etac*n0mean
      bx(i,j)=bux*etac
      by(i,j)=buy*etac
      bz(i,j)=buz
      bfield(j)=sqrt(bx(i,j)**2+by(i,j)**2+bz(i,j)**2)
      bfld=bfield(j)
      bmdx(md)=bx(i,j)
      bmdy(md)=by(i,j)
      bmdz(md)=bz(i,j)
      bmdtot(md)=bfld
c     Calculate the initial maximum electron energy in the Mach disk
      ball2=bux*bux+buy*buy+buz*buz
c     Next line relates this to direction of upstream B field relative to shock
c      gmax0(i,j)=gmaxmx*buz*buz/ball2
c      if(gmax0(i,j).lt.gmaxmn)gmax0(i,j)=gmaxmn
c      gmax0(i,j)=gmaxmn
c     Next 3 lines assume a power-law distribution of gmax0, unrelated
c         to direction of B field
c     See http://mathworld.wolfram.com/RandomNumber.html for choosing
c      random numbers from a power-law distribution
c      xrand=rand(0)
c      if(pexp.eq.0.0)gmax0(i,j)=gmaxmx
c      if(pexp.lt.0.0)gmax0(i,j)=((gmaxmx**pexp-gmaxmn**pexp)*xrand+
c     ,  gmaxmn**pexp)**(1.0/pexp)
      gmax0(i,j)=gmaxmn
      gminmd=gmin
c     Calculate energy distribution in the Mach disk cell
c     Compute energy density of photons in Mach disk, time delayed by 1 step
c       Ignore photons from dust torus since beaming is small in Mach disk
      mdi=md-1
      if(mdi.le.0)mdi=1
      uphmd=bmdtot(mdi)**2/(8.0*pi)*(sqrt(1.0+4.0*uratio)-1.0)/2.0
      ustob=8.0*pi*uphmd/(bfield(j))**2
      tlfact=cc2*bfield(j)**2*(1.0+ustob)
c     Effective length of Mach disk zsmd = loss time of highest-E
c       electrons (to energy glow where radiation is negligible) times flow velocity
      glow=10.0
      delt=(gmax0(i,j)-glow)/(tlfact*glow*gmax0(i,j))
      zsmd(md)=delt/(yr*3.26)*gammd*betamd
c     Use 0.99 instead of 1 to avoid singularity of N(gamma) at gmax0
      gmrat=0.99*gmax0(i,j)/gminmd
      gmratl=dlog10(gmrat)/32.0
      gmratm=dlog10(gminmd/glow)/11.0
      do 1121 ie=1,44
      if(ie.ge.12)egam(i,j,ie)=gminmd*10.0**(gmratl*(ie-12))
      if(ie.lt.12)egam(i,j,ie)=glow*10.0**(gmratm*(ie-1))
      ggam(ie)=egam(i,j,ie)
      t2=(gmax0(i,j)-ggam(ie))/(tlfact*ggam(ie)*gmax0(i,j))
      enofe(i,j,ie)=0.0
      eterm1=ggam(ie)/gminmd
      eterm2=1.0d0-ggam(ie)*t2*tlfact
      if(eterm2.lt.0.0)go to 5124
      if(ie.ge.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(1.0d0-eterm2**(sen-1.0))
      if(ie.lt.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(eterm1**(sen-1.0)-
     ,  eterm2**(sen-1.0))
c     Divide by delt since integral is over time
      delt=zsize*yr*3.26/(gammad(i,j)*betad(i,j))
      enofe(i,j,ie)=enofe(i,j,ie)/delt
      if(enofe(i,j,ie).lt.1.0e-35)enofe(i,j,ie)=0.0
 5124 edist(ie)=enofe(i,j,ie)
 1121 continue
      bperpp=bfld
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 1122 inu=1,40
      restnu=nu(inu)
      alfmdc(md,inu)=10.0
      fsscmd(md,inu)=0.0
c     Synchrotron mean intensity for SSC calculation inside Mach disk
      abssss=1.02e4*(sen+2.0)*akapnu(nu(inu))*bperpp
      if(abssss.lt.(1.0e-35*nu(inu)))abssss=0.0
      sstau(inu)=abssss*(parsec/nu(inu))*rsize*svmd/nu(inu)
c     fsync(inu) is the optically thin intensity
      fsync(inu)=0.0
      fsyncc=ajnu(restnu)
      if(fsyncc.gt.1.0e-30)fsync(inu)=fsyncc*
     ,  bperpp*rsize*svmd*(emfact*parsec)
      if(sstau(inu).lt.0.01)go to 1123
c     Optically thick source function = emis. coef./abs. coef.
c     Note that the length of the path through the source cancels out
      srcfn(inu)=fsync(inu)/sstau(inu)
      if(sstau(inu).gt.5.0)fsync(inu)=srcfn(inu)
      if(sstau(inu).gt.0.01.and.sstau(inu).le.5.0)fsync(inu)=
     ,  srcfn(inu)*(1.0-exp(-sstau(inu)))
 1123 continue
c     Now estimate synchrotron emission seen by other cells
      fsynmd(md,inu)=ajnu(nu(inu)/dopref)*dopref**2*bperpp
      absorb(md,inu)=0.0
      abssss=1.02e4*(sen+2.0)*akapnu(nu(inu)/dopref)*bperpp
      if(abssss.lt.((1.0e-35*nu(inu))*nu(inu)))abssss=0.0
      absorb(md,inu)=(abssss/nu(inu))/nu(inu)
      alphmd(md,inu)=10.0
      abexmd(md,inu)=1.7
      ssseed(inu)=fsync(inu)
 1122 continue
!$omp end parallel do
c     Need to add a lower frequency to get spectral index of nu(1)
      fq1=0.98*nu(1)
      fsyn1=ajnu(fq1/dopref)*dopref**2*bperpp
      absrb1=1.02e4*(sen+2.0)*akapnu(fq1/dopref)*bperpp/(fq1/dopref)**2
c     Compute spectral indices
      do 1124 inu=1,40
      restnu=nu(inu)
      if(fsynmd(md,inu).gt.0.0.and.fsyn1.gt.0.0)
     ,alphmd(md,inu)=-dlog10(fsynmd(md,inu)/fsyn1)/alog10(restnu/fq1)
      if(absorb(md,inu).gt.0.0.and.absrb1.gt.0.0)
     ,abexmd(md,inu)=-dlog10(absorb(md,inu)/absrb1)/alog10(restnu/fq1)
      fq1=restnu
      fsyn1=fsynmd(md,inu)
      absrb1=absorb(md,inu)
 1124 continue
      do 1125 inu=41,68
      alphmd(md,inu)=10.0
      ssseed(inu)=0.0
      fsync(inu)=0.0
      fsynmd(md,inu)=0.0
      absorb(md,inu)=0.0
 1125 continue
c     Calculate the SSC emission from the Mach disk for reference
c       relative Doppler factor dopref
      betd=betamd
      gamd=gammd
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 1126 inu=9,68
      fsscmd(md,inu)=ssc(nu(inu)/dopref)*dopref**2
      alfmdc(md,inu)=10.0
 1126 continue
!$omp end parallel do
      fq1=0.98*nu(9)
      fssc1=ssc(fq1/dopref)*dopref**2
      do 1127 inu=9,68
      restnu=nu(inu)
      if(fsscmd(md,inu).gt.0.0.and.fssc1.gt.0.0)
     ,alfmdc(md,inu)=-alog10(fsscmd(md,inu)/fssc1)/alog10(restnu/fq1)
      fq1=restnu
      fssc1=fsscmd(md,inu)
 1127 continue
      go to 88
   13 continue
c     Calculate component of magnetic field that is perpendicular to
c       the aberrated line of sight in the plasma frame
c     Line-of-sight unit vector
      slx=slos
      sly=0.0
      slz=clos
      call bcalc(betadx(i,j),betady(i,j),betadz(i,j),betau(j),gammau(j),
     ,  slx,sly,slz,bx(i,j),by(i,j),bz(i,j),dum1,dum2,dum3,dum4,
     ,  dum5,dum6,dum7,bperp(j))
      bfield(j)=sqrt(bx(i,j)**2+by(i,j)**2+bz(i,j)**2)
      n0mean=spsd(idelay+ip0)*(anorm0*an00+anorm1*an01+anorm2*an02+
     ,  anorm3*an03)
      n0(i,j)=eta*n0mean
c     Calculate the initial maximum electron energy in each cell
c       from the ratio of B(perp. to shock) to B_total in shock frame
c      gmax0(i,j)=gmaxmx*(bprp*bprp/(bpar*bpar+bprp*bprp))
c      if(gmax0(i,j).lt.gmaxmn)gmax0(i,j)=gmaxmn
      gmax0(i,j)=gmaxmx
c     Maximum electron energy reached only if shock is subluminal
c       in de Hoffman-Teller frame
      htlim=1.0/(gammau(j)*velpar)
      bfang=bprp/bpar
      if(bfang.gt.htlim)gmax0(i,j)=gmaxmn
c     Next 3 lines assume a power-law distribution of gmax0, unrelated
c         to direction of B field
c     See http://mathworld.wolfram.com/RandomNumber.html for choosing
c      random numbers from a power-law distribution
c      xrand=rand(0)
c      if(pexp.eq.0.0)gmax0(i,j)=gmaxmx
c      if(pexp.lt.0.0)gmax0(i,j)=((gmaxmx**pexp-gmaxmn**pexp)*xrand+
c     ,  gmaxmn**pexp)**(1.0/pexp)
      do 15 ig=1,43
   15 if(gmax0(i,j).le.gcnt(ig+1).and.gmax0(i,j).ge.gcnt(ig))
     , igcnt(ig)=igcnt(ig)+1
   80 continue
c
c     End initial loop over all cells in first layer
c
c     Start another loop over all cells in first layer to calculate the
c       electron energy distribution and the emission
c
   88 continue
      do 100 j=jcells,2,-1
      ncells=ncells+1
      do 101 inu=1,68
      emisco(inu)=0.0
      ecflux(inu)=0.0
  101 continue
c     Determine Doppler factor relative to the observer
      bdx=betadx(i,j)
      bdy=betady(i,j)
      bdz=betadz(i,j)
      betacs=bdx*slos+bdz*clos
      delta(i,j)=1.0d0/(gammad(i,j)*(1.0d0-betacs))
      bperpp=bperp(j)
      bfld=bfield(j)
c     Determine velocity of the cell plasma relative to the MD plasma
      bmparx=betamd*bdx*bdz/betad(i,j)**2
      bmpary=betamd*bdy*bdz/betad(i,j)**2
      bmparz=betamd*bdz*bdz/betad(i,j)**2
      bmprpx=-bmparx
      bmprpy=-bmpary
      bmprpz=betamd-bmparz
      betamx(j)=(bdx-bmparx-bmprpx/gammad(i,j))/(1.0d0-betamd*bdz)
      betamy(j)=(bdy-bmpary-bmprpy/gammad(i,j))/(1.0d0-betamd*bdz)
      betamz(j)=(bdz-bmparz-bmprpz/gammad(i,j))/(1.0d0-betamd*bdz)
      betamr(j)=dsqrt(betamx(j)**2+betamy(j)**2+betamz(j)**2)
      gamamr(j)=1.0d0/dsqrt(1.0d0-betamr(j)**2)
      zcl=zcell(i,j)
      zrel=zshock-0.5*zsize-zcl
      dmd(j)=sqrt(zrel**2+rcell(j)**2)
c     Calculate angle between MD seed photon and scattered l.o.s. photon
      call scatcs(bdx,bdy,bdz,clos,0.0d0,slos,xcell(j),ycell(j),
     ,   -zrel,cscat)
c     Determine Doppler factor of the MD emission in the cell's frame
      dmdp=sqrt(rcell(j)**2+(zrel/gamamr(j))**2)
      cosmd=(betamx(j)*xcell(j)+betamy(j)*ycell(j)+
     ,  betamz(j)*zrel/gamamr(j))/(betamr(j)*dmdp)
      tanmd=dsqrt(1.0d0/cosmd**2-1.0d0)
      angm=datan(tanmd)
      tanmd=dtan(2.0*datan(dtan(0.5d0*angm)/
     ,  (gamamr(j)*(1.0d0+betamr(j)))))
      cosmd=1.0d0/dsqrt(1.0d0+tanmd**2)
      deltmd(j)=1.0d0/(gamamr(j)*(1.0d0-betamr(j)*cosmd))
c     Determine time step of Mach disk seed photons from light-travel delay
      delcor=deltmd(j)/dopref
      md=mdmd-0.5-(((zrf-zcl)*clos+(xrf-xcell(j))*slos)+
     ,   dmd(j))*betadd/zsize
      if(md.lt.1.or.md.gt.64000)write(6,9992)i,j,md,zcell(i,j),xcell(j),
     ,   dmd(j),clos,slos,betadd,zrf,xrf,zshock,zsize
      do 3146 inu=1,68
      fmdall(inu)=0.0
      syseed(inu)=0.0
      scseed(inu)=0.0
      tauxmd(inu)=0.0
 3146 ssseed(inu)=0.0
c     Mach disk's B field in frame of cell plasma transverse to cell's
c       line of sight to Mach disk (for synchrotron calculation)
      sx=rcell(j)*cosph(j)
      sy=rcell(j)*sinph(j)
      sz=zrel
      call bcalc(betamx(j),betamy(j),betamz(j),betamr(j),gamamr(j),
     ,  sx,sy,sz,bmdx(md),bmdy(md),bmdz(md),dum1,dum2,dum3,dum4,
     ,  dum5,dum6,dum7,bmperp)
c     Apply as a correction factor to previous estimate of B_perpendicular
      bpcorr=bmperp/bmdtot(md)
c
c     Calculate emission from various processes
c
c     Calculate synchrotron flux from Mach disk in frame of cell plasma
c
      path=2.0*rsize*svmd/(tanmd*cosmd)
      if(path.gt.zsmd(md))path=zsmd(md)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 146 inu=1,40
      sstau(inu)=0.0
      tauexp(inu)=0.0
      ssabs(inu)=0.0
      ajofnu(inu)=0.0
      if(fsynmd(md,inu).lt.1.0e-30.or.alphmd(md,inu).gt.9.0)
     ,   go to 145
c     Synchrotron mean intensity for inverse Compton calculation
      ssabs(inu)=absorb(md,inu)*bpcorr**(0.5*(sen+2.0))*
     ,  delcor**(abexmd(md,inu))
      sstau(inu)=ssabs(inu)*parsec*rsize*svmd
      ajofnu(inu)=fsynmd(md,inu)*bpcorr**(1.0+alphmd(md,inu))*
     ,  delcor**(2.0+alphmd(md,inu))
      ajofnu(inu)=ajofnu(inu)*zsmd(md)*(emfact*parsec)
      syseed(inu)=ajofnu(inu)
      if(sstau(inu).lt.0.01)go to 143
      srcfn(inu)=ajofnu(inu)/sstau(inu)
      if(sstau(inu).gt.5.0)syseed(inu)=srcfn(inu)
      if(sstau(inu).gt.0.1.and.sstau(inu).le.5.0)syseed(inu)=
     ,  srcfn(inu)*(1.0-exp(-sstau(inu)))
  143 syseed(inu)=syseed(inu)*pi*rsize**2*vmd/(dmd(j)**2)
c     Now estimate exponential attenuation of MD seed photons by
c      synchrotron self-absorption in intervening cells
      tauexp(inu)=1.02e4*(sen+2.0)*akapnu(nu(inu))*
     ,  bperp(j)/nu(inu)*(parsec/nu(inu))*dmd(j)
      tauxmd(inu)=tauexp(inu)
      if(tauexp(inu).gt.15.0)syseed(inu)=0.0
      if(tauexp(inu).le.15.0)syseed(inu)=syseed(inu)/exp(tauexp(inu))
  145 continue
      if(inu.lt.16)fmdall(inu)=fmdall(inu)+syseed(inu)
  146 continue
!$omp end parallel do
c
c     Calculate inverse Compton flux from Mach disk in frame of cell plasma
c
      do 147 inu=16,68
c     Inverse Compton mean intensity from Mach disk for 2nd-order
c       inverse Compton calculation in cells
      scseed(inu)=fsscmd(md,inu)*
     ,  delcor**(2.0+alfmdc(md,inu))
      scseed(inu)=scseed(inu)*zsmd(md)*pi*rsize**2*
     ,  (amjy*parsec)*vmd/dmd(j)**2
      if(tauxmd(inu).gt.15.0)scseed(inu)=0.0
      if(tauxmd(inu).le.15.0)scseed(inu)=scseed(inu)/
     ,exp(tauxmd(inu))
      fmdall(inu)=fmdall(inu)+syseed(inu)+scseed(inu)
  147 continue
      nuhi=2
      do 2148 inu=1,68
      ssseed(inu)=fmdall(inu)
 2148 if(fmdall(inu).gt.0.0)nuhi=inu
c     Calculate seed photon energy density energy loss calculation
      usdmd=0.0
      do 149 inu=2,nuhi
c     Only calculate up to Klein-Nishina limit of gmin
      epslon=6.63e-27*snu(inu)/(gmin*emc2)
      if(epslon.gt.1.0)go to 150
      if(fmdall(inu).le.1.0e-30.or.fmdall(inu-1).le.1.0e-30)go to 148
      aaa=dlog10(fmdall(inu-1)/fmdall(inu))/alog10(snu(inu)/snu(inu-1))
      if(aaa.gt.0.999.and.aaa.lt.1.0001)go to 2147
      usdmd=usdmd+0.5/c*snu(inu-1)/(1.0-aaa)*fmdall(inu-1)*
     ,  ((snu(inu)/snu(inu-1))**(1.0-aaa)-1.0)
      go to 2149
 2147 usdmd=usdmd+0.5*fmdall(inu-1)*snu(inu-1)*
     ,  alog(snu(inu)/snu(inu-1))/c
 2149 continue
      if(usdmd.ne.usdmd)write(6,9994)inu,j,aaa,snu(inu-1),snu(inu),
     ,  fmdall(inu-1),fmdall(inu)
  148 continue
  149 continue
  150 continue
c     Calculate energy distribution in the cell
      id=zcell(i,j)/zsize+1
      if(id.lt.1)id=1
      if(id.gt.(icells))write(6,6992)it,i,j,id,md,k,zcell(i,j),
     ,  zshock,zsize
 6992 format('*x*x*x ',6i6,15f11.5)
      zdist=zdist0+(id-1)*zsize
c     Calculate min & max angles of dust torus in plasma frame
      dphi1=asin(dtrad/sqrt(zdist**2+dtdist**2))
      dphi2=atan(dtdist/zdist)
      dth1=dphi2-dphi1
      dth2=dphi2+dphi1
      csth1=dcos(dth1)
      csth2=dcos(dth2)
      dcsth1=-(csth1-betad(i,j))/(1.0d0-csth1*betad(i,j))
      dcsth2=-(csth2-betad(i,j))/(1.0d0-csth2*betad(i,j))
      csang=0.5*(dcsth1+dcsth2)
c     Doppler factor of dust torus emission in frame of cells, used
c       to estimate frequency of peak intensity in plasma frame
      tdel=gammad(i,j)*(1.0d0-betad(i,j)*0.5d0*(csth1+csth2))
c     Correct energy density of seed photons from dust for Doppler factor
c       that include turbulent component of velocity; approximates that
c       motion is directed parallel to jet axis; also include first-order
c       dependence on scattering angle in Klein-Nishina cross-section
      useedr=useed(id)*(tdel/tdelr(id))**2/(1.0d0+betad(i,j)*csang)
      ustob=8.0*pi*(useedr+usdmd)/(bfield(j))**2
      tlfact=cc2*bfield(j)**2*(1.0+ustob)
      tlf1(j)=tlfact
      delt=zsize*yr*3.26/(gammad(i,j)*betad(i,j))
      gamb=gmax0(i,j)/(1.0+tlfact*gmax0(i,j)*delt)
      glow=gmin/(1.0+tlfact*gmin*delt)
      gmrat=0.99*gmax0(i,j)/gmin
      gmratl=dlog10(gmrat)/32.0
      gmratm=dlog10(gmin/glow)/11.0
      ibreak=0
      do 90 ie=1,44
      egam(i,j,ie)=gmin*10.0**(gmratl*(ie-12))
      if(ie.lt.12)egam(i,j,ie)=glow*10.0**(gmratm*(ie-1))
      if(ie.eq.1)egam(i,j,ie)=glow+0.2*(glow*10.0**(gmratm)-glow)
      if(ibreak.eq.1.or.egam(i,j,ie).lt.gamb)go to 89
      egam(i,j,ie)=gamb
      ibreak=1
   89 ggam(ie)=egam(i,j,ie)
      tloss=(gmax0(i,j)-ggam(ie))/(tlfact*ggam(ie)*gmax0(i,j))
      t2=dmin1(tloss,delt)
      enofe(i,j,ie)=0.0
      eterm1=ggam(ie)/gmin
      eterm2=1.0d0-ggam(ie)*t2*tlfact
      if(eterm2.lt.0.0)go to 5089
      if(ie.ge.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(1.0d0-eterm2**(sen-1.0))
      if(ie.lt.12)enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(eterm1**(sen-1.0)-
     ,  eterm2**(sen-1.0))
c     Divide by delt since integral is over time
      enofe(i,j,ie)=enofe(i,j,ie)/delt
      if(enofe(i,j,ie).lt.1.0e-35)enofe(i,j,ie)=0.0
 5089 edist(ie)=enofe(i,j,ie)
   90 continue
      gammax(i,j)=gmax0(i,j)/(1.0d0+tlfact*delt*gmax0(i,j))
      gammin(i,j)=gmin/(1.0d0+tlfact*delt*gmin)
      emold=0.0
c     Start loop over observed frequencies
      do 91 inu=1,22
      dustnu(inu)=dustf(id,inu)
      hnukt=4.8e-11*dustnu(inu)/tdust
      hnuktr=4.8e-11*dustnu(inu)*(tdel/tdelr(id))/tdust
      dusti(inu)=dustii(id,inu)
      if(hnuktr.gt.60.0)go to 5093
      if(dusti(inu).gt.1.0e-30)dusti(inu)=dusti(inu)*
     ,  (tdel/tdelr(id))**3*(1.0-exp(hnukt))/(1.0-exp(hnuktr))
      go to 5094
 5093 dusti(inu)=0.0
 5094 continue
   91 continue
      betd=betad(i,j)
      gamd=gammad(i,j)
      call polcalc(bfield(j),bx(i,j),by(i,j),bz(i,j),clos,slos,chipol)
      if(chipol.lt.-pio2)chipol=chipol+pi
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 92 inu=1,68
      rnu(inu)=nu(inu)*zred1/delta(i,j)
      specin(inu)=0.0001
      fsync2(inu)=0.0
      srcfn(inu)=0.0
      fsnoab(inu)=0.0
      ssabs(inu)=0.0
      ecflux(inu)=0.0
      sscflx(inu)=0.0
      ssseed(inu)=0.0
      flsync(i,j,inu)=0.0
      flcomp(i,j,inu)=0.0
      flec(i,j,inu)=0.0
      flssc(i,j,inu)=0.0
      flux(i,j,inu)=0.0
      specin(inu)=alpha
      if(inu.gt.44)go to 93
      emisco(inu)=ajnu(rnu(inu))*bperp(j)*delta(i,j)**2
      fsync2(inu)=emisco(inu)*zsize*(emfact*parsec)
      fsnoab(inu)=fsync2(inu)
c     Skip frequencies that are obviously optically thin
      if(inu.gt.15)go to 93
      sssab=akapnu(rnu(inu))
      if(sssab.gt.1.0e-30)ssabs(inu)=1.02e4*(sen+2.0)*sssab*bperp(j)/
     ,   rnu(inu)*(parsec/rnu(inu))/delta(i,j)
      ssabs(inu)=ssabs(inu)*zsize*delta(i,j)
c     Attenuate Mach disk emission from synchrotron absorption on way to cell
      ssabsm(inu)=ssabs(inu)*dmd(j)/zsize
      if(ssabsm(inu).ge.10.0)ssseed(inu)=0.0
      if(ssabsm(inu).lt.10.0.and.ssabsm(inu).gt.0.1)ssseed(inu)=
     ,     ssseed(inu)/exp(ssabsm(inu))
c     Return to absorption within cell
c     Use rsize instead of zsize because of aberration
      ssabs(inu)=ssabs(inu)*rsize/zsize
      fsync2(inu)=fsync2(inu)*rsize/zsize
      if(ssabs(inu).gt.0.01)srcfn(inu)=fsync2(inu)/ssabs(inu)
      if(ssabs(inu).gt.5.0)fsync2(inu)=srcfn(inu)
      if(ssabs(inu).gt.0.1.and.ssabs(inu).le.5.0)fsync2(inu)=
     ,  srcfn(inu)*(1.0-exp(-ssabs(inu)))
      tauexp(inu)=(nouter(j)-i)*ssabs(inu)
      if(rcell(j).gt.(0.98*rbound).and.xcell(j).le.0.0)
     ,         tauexp(inu)=0.0
      if(thlos.eq.0.0)tauexp(inu)=(nouter(j)-i)*ssabs(inu)
      if(tauexp(inu).gt.15.0)fsync2(inu)=0.0
      if(tauexp(inu).gt.0.01.and.tauexp(inu).le.15.0)fsync2(inu)=
     ,   fsync2(inu)/exp(tauexp(inu))
c     Revert to using zsize since next step divides intensity by zsize
      fsync2(inu)=fsync2(inu)*zsize/rsize
      specin(inu)=alpha
   93 continue
      flsync(i,j,inu)=fsync2(inu)*(volc/zsize)*zred1/
     ,  (1.0e18*amjy*dgpc**2)*fgeom
      flux(i,j,inu)=flsync(i,j,inu)
   92 continue
!$omp end parallel do
      fq1=0.98*rnu(1)
      emold=ajnu(fq1)*bperp(j)*delta(i,j)**2*zsize*(emfact*parsec)
      fq1=0.98*nu(1)
      do 94 inu=1,43
      if(emold.gt.0.0.and.fsnoab(inu).gt.0.0)specin(inu)=
     ,  dlog10(emold/fsnoab(inu))/alog10(nu(inu)/fq1)
      if(specin(inu).lt.alpha)specin(inu)=alpha
c     Determine degree of polarization
      poldeg=(specin(inu)+1.0)/(specin(inu)+5.0/3.0)
      if(ssabs(inu).gt.1.0)poldeg=3.0/(12.0*specin(inu)+19.0)
      if(ssabs(inu).gt.1.0)chipol=chipol+pio2
      if(chipol.gt.pi)chipol=chipol-pi
      fpol(i,j,inu)=poldeg*flsync(i,j,inu)
      pq(i,j,inu)=fpol(i,j,inu)*cos(2.0*chipol)
      pu(i,j,inu)=fpol(i,j,inu)*sin(2.0*chipol)
      fq1=nu(inu)
      emold=fsnoab(inu)
   94 continue
      betd=betad(1,1)
      gamd=gammad(1,1)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 95 inu=14,68
      ecflux(inu)=0.0
      sscflx(inu)=0.0
      if(rnu(inu).lt.1.0e14)go to 85
      ecflux(inu)=ecdust(rnu(inu))*3.086*volc*zred1*delta(i,j)**2/
     ,   dgpc**2*fgeom
      sscflx(inu)=ssc(rnu(inu))*3.086*volc*zred1*delta(i,j)**2/
     ,   dgpc**2*fgeom
   85 taupp(inu)=0.0
   95 continue
!$omp end parallel do
      do 86 inu=14,68
      if(nu(inu).lt.1.0e22.or.ipprod.eq.0)go to 98
c     Pair production opacity calculation
c     Expression for anumin includes typical interaction angle
      anumin=(1.24e20/(nu(inu)*zred1))*1.24e20*(2.0*gammad(i,j))**2
      alnumn=alog10(anumin)
      inumin=(alnumn-10.0)*4+1
      if(inumin.gt.40)go to 98
      if(anumin.gt.nu(inumin))anumin=nu(inumin)
      do 96 iinu=inumin,40
      bp=sqrt(1.0-anumin/nu(iinu))
      if(bp.le.0.0)bp=0.001
      xsecp1=1.25e-25*(1.0-bp*bp)*((3.0-bp**4)*
     *  alog((1.0+bp)/(1.0-bp))-2.0*bp*(2.0-bp*bp))
      bp=sqrt(1.0-anumin/nu(iinu+1))
      if(bp.le.0.0)bp=0.001
      xsecp2=1.25e-25*(1.0-bp*bp)*((3.0-bp**4)*
     *  alog((1.0+bp)/(1.0-bp))-2.0*bp*(2.0-bp*bp))
      taupp(inu)=taupp(inu)+0.5*(phots(iinu)*xsecp1+
     ,  phots(iinu+1)*xsecp2)*
     ,  (nu(iinu+1)-nu(iinu))*(nouter(j)-i)*zsize*parsec
   96 continue
      if(taupp(inu).lt.0.1)go to 98
      if(taupp(inu).lt.10.0)go to 97
   87 ecflux(inu)=0.0
      sscflx(inu)=0.0
      go to 98
   97 ecflux(inu)=ecflux(inu)/exp(taupp(inu))
      sscflx(inu)=sscflx(inu)/exp(taupp(inu))
   98 flec(i,j,inu)=ecflux(inu)
      flssc(i,j,inu)=sscflx(inu)
      flcomp(i,j,inu)=ecflux(inu)+sscflx(inu)
      flux(i,j,inu)=flsync(i,j,inu)+flcomp(i,j,inu)
   86 continue
      do 99 inu=15,68
      if(ecflux(inu-1).gt.0.0.and.ecflux(inu).gt.0.0)spxec(inu)=
     ,  dlog10(ecflux(inu-1)/ecflux(inu))/alog10(nu(inu)/nu(inu-1))
      if(sscflx(inu-1).gt.0.0.and.sscflx(inu).gt.0.0)spxssc(inu)=
     ,  dlog10(sscflx(inu-1)/sscflx(inu))/alog10(nu(inu)/nu(inu-1))
   99 continue
ccc    Test
  100 continue
c
c     End loop over first layer of cells
c
c      go to 2200
      if(icells.eq.1)go to 2200
      i=i+1
      istart=i
c     Start loop over downstream cells
      do 200 j=jcells,2,-1
      do 200 i=istart,imax(j)
      ncells=ncells+1
      if(it.gt.1)go to 110
c
c     *** Initial set-up of downstream cells; skip after 1st time step
c
      k=kmd-(zshock-zcell(i,j))/zsize
c     Set up physical parameters of downstream cells at first time step
c     Cells farther downstream contain plasma ejected earlier
      idelay=dstart+it-0.5+((zrf-zcell(i,j))*clos
     ,   +(xrf-xcell(j))*slos)*betadd/zsize
      if(idelay.lt.1.or.idelay.gt.(ndim-ip0))
     ,  write(5,9223)idelay,dstart,j,i,md,zshock,zrf,zcell(i,j),
     ,  xrf,xcell(j),betad(i,j),zsize
c     Velocity vector of laminar component of pre-shock flow
      betupx=betaup*cosph(j)*sinpsi(j)
      betupy=betaup*sinph(j)*sinpsi(j)
      betupz=betaup*cospsi(j)
cccccccccc
c     Determine turbulent velocity component & B vector of cell assuming random magnetic
c       field orientation in each of 4 sections of different sizes
c     For Mach disk, continue from previous calculation of B direction
c     Interpolate inside loop to get direction for cells not on a section boundary
ccccc   Determine field component for each section of cells
c     Start with single cell-sized section 0
      phi10=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet10=acos(costh)
      sint1=sin(thet10)
      call randln(1.0,pb,rndb10)
      call randln(1.0,pn,rndn10)
      bux0=sint1*cos(phi10)*rndb10
      buy0=sint1*sin(phi10)*rndb10
      buz0=costh*rndb10
      an00=n0ave*rndn10
      phv10=2.0*pi*rand(0)
      cosv=2.0*(rand(0)-0.5)
      if(cosv.ge.1.0)cosv=0.999999
      if(cosv.le.-1.0)cosv=-0.999999
      thv10=acos(cosv)
      sinv=sin(thv10)
      vtz0=cosv
      vtx0=sinv*cos(phv10)
      vty0=sinv*sin(phv10)
c     Next is section 1, equal in size to 2x2=4 cells
      if(ididj1(j).eq.1)go to 6238
      if(ididk1(j,layer(j)).eq.0) go to 6241
 6238 if(ididb1(ng1x(j),ng1y(j)).eq.1)go to 6240
      phi11(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet11(ng1x(j),ng1y(j))=acos(costh)
      sint1=sin(thet11(ng1x(j),ng1y(j)))
      sint2=sin(thet21(ng1x(j),ng1y(j)))
      phv11(ng1x(j),ng1y(j))=2.0*pi*rand(0)
      cosv11(ng1x(j),ng1y(j))=2.0*(rand(0)-0.5)
      if(cosv11(ng1x(j),ng1y(j)).ge.1.0)cosv11(ng1x(j),ng1y(j))=0.999999
      if(cosv11(ng1x(j),ng1y(j)).le.-1.0)cosv11(ng1x(j),ng1y(j))=
     ,  -0.999999
      thv=acos(cosv11(ng1x(j),ng1y(j)))
      sinv11(ng1x(j),ng1y(j))=sin(thv)
      call randln(1.0,pb,rndb11)
      call randln(1.0,pn,rndn11)
      bu1x1(ng1x(j),ng1y(j))=sint1*cos(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1y1(ng1x(j),ng1y(j))=sint1*sin(phi11(ng1x(j),ng1y(j)))*rndb11
      bu1z1(ng1x(j),ng1y(j))=costh*rndb11
      bu2x1(ng1x(j),ng1y(j))=sint2*cos(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2y1(ng1x(j),ng1y(j))=sint2*sin(phi21(ng1x(j),ng1y(j)))*rndb21
      bu2z1(ng1x(j),ng1y(j))=cos(thet21(ng1x(j),ng1y(j)))*rndb21
      cosrt1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))+
     ,  bu1y1(ng1x(j),ng1y(j))*bu2y1(ng1x(j),ng1y(j))+
     ,  bu1z1(ng1x(j),ng1y(j))*bu2z1(ng1x(j),ng1y(j))
      cdenom=sqrt((bu1x1(ng1x(j),ng1y(j))**2+bu1y1(ng1x(j),ng1y(j))**2
     ,  +bu1z1(ng1x(j),ng1y(j))**2)*(bu2x1(ng1x(j),ng1y(j))**2+
     ,   bu2y1(ng1x(j),ng1y(j))**2+bu2z1(ng1x(j),ng1y(j))**2))
      cosrt1(ng1x(j),ng1y(j))=cosrt1(ng1x(j),ng1y(j))/cdenom
      if(cosrt1(ng1x(j),ng1y(j)).ge.1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=0.999999
      if(cosrt1(ng1x(j),ng1y(j)).le.-1.0)
     ,  cosrt1(ng1x(j),ng1y(j))=-0.999999
      angrt1(ng1x(j),ng1y(j))=acos(cosrt1(ng1x(j),ng1y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt1(ng1x(j),ng1y(j))=
     ,  angrt1(ng1x(j),ng1y(j))-2.0*pi
      cu1x1(ng1x(j),ng1y(j))=bu1y1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))-bu1z1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))
      cu1y1(ng1x(j),ng1y(j))=-bu1x1(ng1x(j),ng1y(j))*
     ,  bu2z1(ng1x(j),ng1y(j))+bu1z1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      cu1z1(ng1x(j),ng1y(j))=bu1x1(ng1x(j),ng1y(j))*
     ,  bu2y1(ng1x(j),ng1y(j))-bu1y1(ng1x(j),ng1y(j))*
     ,  bu2x1(ng1x(j),ng1y(j))
      ididb1(ng1x(j),ng1y(j))=1
 6240 continue
      anggrr=bfrac1(ng1x(j),ng1y(j))*angrt1(ng1x(j),ng1y(j))
      call vecrot(bu1x1(ng1x(j),ng1y(j)),bu1y1(ng1x(j),ng1y(j)),
     ,  bu1z1(ng1x(j),ng1y(j)),cu1x1(ng1x(j),ng1y(j)),
     ,  cu1y1(ng1x(j),ng1y(j)),cu1z1(ng1x(j),ng1y(j)),anggrr,
     ,  bux1(ng1x(j),ng1y(j),k),buy1(ng1x(j),ng1y(j),k),
     ,  buz1(ng1x(j),ng1y(j),k))
      an01=n0ave*(rndn11*bfrac1(ng1x(j),ng1y(j))+
     ,         rndn21*(1.0-bfrac1(ng1x(j),ng1y(j))))
      dvx=xcell(j)-(0.5*ng1x(j)-2.25)*nrad*rsize
      dvy=ycell(j)-(0.5*ng1y(j)-2.25)*nrad*rsize
      dvz=(0.5-bfrac1(ng1x(j),ng1y(j)))*zsize
c     If nrad=8, adjust placement of zone center so it is not at a cell center
      if(nrad.eq.8)then
        if(ng1x(j).lt.5)dvx=dvx+rsize
        if(ng1x(j).ge.5)dvx=dvx-rsize
        if(ng1y(j).lt.5)dvy=dvy+rsize
        if(ng1y(j).ge.5)dvy=dvy-rsize
        endif
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv11(ng1x(j),ng1y(j))
      axx=sinv11(ng1x(j),ng1y(j))*cos(phv11(ng1x(j),ng1y(j)))
      axy=sinv11(ng1x(j),ng1y(j))*sin(phv11(ng1x(j),ng1y(j)))
      vtx1=(axy*dvz-axz*dvy)/dvs
      vty1=(axz*dvx-axx*dvz)/dvs
      vtz1=(axx*dvy-axy*dvx)/dvs
      bfrac1(ng1x(j),ng1y(j))=bfrac1(ng1x(j),ng1y(j))+0.99999
      if(bfrac1(ng1x(j),ng1y(j)).lt.1.0)go to 6241
      thet21(ng1x(j),ng1y(j))=thet11(ng1x(j),ng1y(j))
      phi21(ng1x(j),ng1y(j))=phi11(ng1x(j),ng1y(j))
      bfrac1(ng1x(j),ng1y(j))=0.0
      ididb1(ng1x(j),ng1y(j))=0
 6241 continue
c     Next is section 2, equal in size to 4x4=16 cells
      if(ididj2(j).eq.1)go to 6242
      if(ididk2(j,layer(j)).eq.0) go to 6251
 6242 if(ididb2(ng2x(j),ng2y(j)).eq.1)go to 6250
      phi12(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet12(ng2x(j),ng2y(j))=acos(costh)
      sint1=sin(thet12(ng2x(j),ng2y(j)))
      sint2=sin(thet22(ng2x(j),ng2y(j)))
      phv12(ng2x(j),ng2y(j))=2.0*pi*rand(0)
      cosv12(ng2x(j),ng2y(j))=2.0*(rand(0)-0.5)
      if(cosv12(ng2x(j),ng2y(j)).ge.1.0)cosv12(ng2x(j),ng2y(j))=0.999999
      if(cosv12(ng2x(j),ng2y(j)).le.-1.0)cosv12(ng2x(j),ng2y(j))=
     ,  -0.999999
      thv=acos(cosv12(ng2x(j),ng2y(j)))
      sinv12(ng2x(j),ng2y(j))=sin(thv)
      call randln(1.0,pb,rndb12)
      call randln(1.0,pn,rndn12)
      bu1x2(ng2x(j),ng2y(j))=sint1*cos(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1y2(ng2x(j),ng2y(j))=sint1*sin(phi12(ng2x(j),ng2y(j)))*rndb12
      bu1z2(ng2x(j),ng2y(j))=costh*rndb12
      bu2x2(ng2x(j),ng2y(j))=sint2*cos(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2y2(ng2x(j),ng2y(j))=sint2*sin(phi22(ng2x(j),ng2y(j)))*rndb22
      bu2z2(ng2x(j),ng2y(j))=cos(thet22(ng2x(j),ng2y(j)))*rndb22
      cosrt2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))+
     ,  bu1y2(ng2x(j),ng2y(j))*bu2y2(ng2x(j),ng2y(j))+
     ,  bu1z2(ng2x(j),ng2y(j))*bu2z2(ng2x(j),ng2y(j))
      cdenom=sqrt((bu1x2(ng2x(j),ng2y(j))**2+bu1y2(ng2x(j),ng2y(j))**2
     ,  +bu1z2(ng2x(j),ng2y(j))**2)*(bu2x2(ng2x(j),ng2y(j))**2+
     ,   bu2y2(ng2x(j),ng2y(j))**2+bu2z2(ng2x(j),ng2y(j))**2))
      cosrt2(ng2x(j),ng2y(j))=cosrt2(ng2x(j),ng2y(j))/cdenom
      if(cosrt2(ng2x(j),ng2y(j)).ge.1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=0.99999
      if(cosrt2(ng2x(j),ng2y(j)).le.-1.0)
     ,  cosrt2(ng2x(j),ng2y(j))=-0.99999
      angrt2(ng2x(j),ng2y(j))=acos(cosrt2(ng2x(j),ng2y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt2(ng2x(j),ng2y(j))=
     ,  angrt2(ng2x(j),ng2y(j))-2.0*pi
      cu1x2(ng2x(j),ng2y(j))=bu1y2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))-bu1z2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))
      cu1y2(ng2x(j),ng2y(j))=-bu1x2(ng2x(j),ng2y(j))*
     ,  bu2z2(ng2x(j),ng2y(j))+bu1z2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      cu1z2(ng2x(j),ng2y(j))=bu1x2(ng2x(j),ng2y(j))*
     ,  bu2y2(ng2x(j),ng2y(j))-bu1y2(ng2x(j),ng2y(j))*
     ,  bu2x2(ng2x(j),ng2y(j))
      ididb2(ng2x(j),ng2y(j))=1
 6250 continue
      anggrr=bfrac2(ng2x(j),ng2y(j))*angrt2(ng2x(j),ng2y(j))
      call vecrot(bu1x2(ng2x(j),ng2y(j)),bu1y2(ng2x(j),ng2y(j)),
     ,  bu1z2(ng2x(j),ng2y(j)),cu1x2(ng2x(j),ng2y(j)),
     ,  cu1y2(ng2x(j),ng2y(j)),cu1z2(ng2x(j),ng2y(j)),anggrr,
     ,  bux2(ng2x(j),ng2y(j),k),buy2(ng2x(j),ng2y(j),k),
     ,  buz2(ng2x(j),ng2y(j),k))
      an02=n0ave*(rndn12*bfrac2(ng2x(j),ng2y(j))+
     ,         rndn22*(1.0-bfrac2(ng2x(j),ng2y(j))))
      dvx=xcell(j)-(ng2x(j)-2.5)*nrad*rsize
      dvy=ycell(j)-(ng2y(j)-2.5)*nrad*rsize
      dvz=3.0*(0.5-bfrac2(ng2x(j),ng2y(j)))*zsize
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv12(ng2x(j),ng2y(j))
      axx=sinv12(ng2x(j),ng2y(j))*cos(phv12(ng2x(j),ng2y(j)))
      axy=sinv12(ng2x(j),ng2y(j))*sin(phv12(ng2x(j),ng2y(j)))
      vtx2=(axy*dvz-axz*dvy)/dvs
      vty2=(axz*dvx-axx*dvz)/dvs
      vtz2=(axx*dvy-axy*dvx)/dvs
      bfrac2(ng2x(j),ng2y(j))=bfrac2(ng2x(j),ng2y(j))+0.33333
      if(bfrac2(ng2x(j),ng2y(j)).lt.1.0)go to 6251
      thet22(ng2x(j),ng2y(j))=thet12(ng2x(j),ng2y(j))
      phi22(ng2x(j),ng2y(j))=phi12(ng2x(j),ng2y(j))
      bfrac2(ng2x(j),ng2y(j))=0.0
      ididb2(ng2x(j),ng2y(j))=0
 6251 continue
c     Last is section 2, equal in size to 8x8=64 cells
      if(ididj3(j).eq.1) go to 6252
      if(ididk3(j,layer(j)).eq.0) go to 6261
 6252 if(ididb3(ng3x(j),ng3y(j)).eq.1)go to 6260
      phi13(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      costh=2.0*(rand(0)-0.5)
      if(costh.ge.1.0)costh=0.999999
      if(costh.le.-1.0)costh=-0.999999
      thet13(ng3x(j),ng3y(j))=acos(costh)
      sint1=sin(thet13(ng3x(j),ng3y(j)))
      sint2=sin(thet23(ng3x(j),ng3y(j)))
      phv13(ng3x(j),ng3y(j))=2.0*pi*rand(0)
      cosv13(ng3x(j),ng3y(j))=2.0*(rand(0)-0.5)
      if(cosv13(ng3x(j),ng3y(j)).ge.1.0)cosv13(ng3x(j),ng3y(j))=0.999999
      if(cosv13(ng3x(j),ng3y(j)).le.-1.0)cosv13(ng3x(j),ng3y(j))=
     ,  -0.999999
      thv=acos(cosv13(ng3x(j),ng3y(j)))
      sinv13(ng3x(j),ng3y(j))=sin(thv)
      call randln(1.0,pb,rndb13)
      call randln(1.0,pn,rndn13)
      bu1x3(ng3x(j),ng3y(j))=sint1*cos(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1y3(ng3x(j),ng3y(j))=sint1*sin(phi13(ng3x(j),ng3y(j)))*rndb13
      bu1z3(ng3x(j),ng3y(j))=costh*rndb13
      bu2x3(ng3x(j),ng3y(j))=sint2*cos(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2y3(ng3x(j),ng3y(j))=sint2*sin(phi23(ng3x(j),ng3y(j)))*rndb23
      bu2z3(ng3x(j),ng3y(j))=cos(thet23(ng3x(j),ng3y(j)))*rndb23
      cosrt3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))+
     ,  bu1y3(ng3x(j),ng3y(j))*bu2y3(ng3x(j),ng3y(j))+
     ,  bu1z3(ng3x(j),ng3y(j))*bu2z3(ng3x(j),ng3y(j))
      cdenom=sqrt((bu1x3(ng3x(j),ng3y(j))**2+bu1y3(ng3x(j),ng3y(j))**2
     ,  +bu1z3(ng3x(j),ng3y(j))**2)*(bu2x3(ng3x(j),ng3y(j))**2+
     ,   bu2y3(ng3x(j),ng3y(j))**2+bu2z3(ng3x(j),ng3y(j))**2))
      cosrt3(ng3x(j),ng3y(j))=cosrt3(ng3x(j),ng3y(j))/cdenom
      if(cosrt3(ng3x(j),ng3y(j)).ge.1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=0.999999
      if(cosrt3(ng3x(j),ng3y(j)).le.-1.0)
     ,  cosrt3(ng3x(j),ng3y(j))=-0.999999
      angrt3(ng3x(j),ng3y(j))=acos(cosrt3(ng3x(j),ng3y(j)))
      xsign=rand(0)-0.5
      if(xsign.ge.0.0)xsign=1.0
      if(xsign.lt.0.0)xsign=-1.0
      if(xsign.lt.0.0)angrt3(ng3x(j),ng3y(j))=
     ,  angrt3(ng3x(j),ng3y(j))-2.0*pi
      cu1x3(ng3x(j),ng3y(j))=bu1y3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))-bu1z3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))
      cu1y3(ng3x(j),ng3y(j))=-bu1x3(ng3x(j),ng3y(j))*
     ,  bu2z3(ng3x(j),ng3y(j))+bu1z3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      cu1z3(ng3x(j),ng3y(j))=bu1x3(ng3x(j),ng3y(j))*
     ,  bu2y3(ng3x(j),ng3y(j))-bu1y3(ng3x(j),ng3y(j))*
     ,  bu2x3(ng3x(j),ng3y(j))
      ididb3(ng3x(j),ng3y(j))=1
 6260 continue
      anggrr=bfrac3(ng3x(j),ng3y(j))*angrt3(ng3x(j),ng3y(j))
      call vecrot(bu1x3(ng3x(j),ng3y(j)),bu1y3(ng3x(j),ng3y(j)),
     ,  bu1z3(ng3x(j),ng3y(j)),cu1x3(ng3x(j),ng3y(j)),
     ,  cu1y3(ng3x(j),ng3y(j)),cu1z3(ng3x(j),ng3y(j)),anggrr,
     ,  bux3(ng3x(j),ng3y(j),k),buy3(ng3x(j),ng3y(j),k),
     ,  buz3(ng3x(j),ng3y(j),k))
      an03=n0ave*(rndn13*bfrac3(ng3x(j),ng3y(j))+
     ,         rndn23*(1.0-bfrac3(ng3x(j),ng3y(j))))
      dvx=xcell(j)-(2.0*ng3x(j)-3.0)*sq2*(nrad-1)*rsize
      dvy=ycell(j)-(2.0*ng3y(j)-3.0)*sq2*(nrad-1)*rsize
      dvz=7.0*(0.5-bfrac3(ng3x(j),ng3y(j)))*zsize
      dvs=sqrt(dvx*dvx+dvy*dvy+dvz*dvz)
      axz=cosv13(ng3x(j),ng3y(j))
      axx=sinv13(ng3x(j),ng3y(j))*cos(phv13(ng3x(j),ng3y(j)))
      axy=sinv13(ng3x(j),ng3y(j))*sin(phv13(ng3x(j),ng3y(j)))
      vtx3=(axy*dvz-axz*dvy)/dvs
      vty3=(axz*dvx-axx*dvz)/dvs
      vtz3=(axx*dvy-axy*dvx)/dvs
      bfrac3(ng3x(j),ng3y(j))=bfrac3(ng3x(j),ng3y(j))+0.14285
      if(bfrac3(ng3x(j),ng3y(j)).lt.1.0)go to 6261
      thet23(ng3x(j),ng3y(j))=thet13(ng3x(j),ng3y(j))
      phi23(ng3x(j),ng3y(j))=phi13(ng3x(j),ng3y(j))
      bfrac3(ng3x(j),ng3y(j))=0.0
      ididb3(ng3x(j),ng3y(j))=0
 6261 continue
c     Now add B components of all sections containing cell j
c     Weight according to section size a la scheme of Jones (1988)
      bux=bnorm0*bux0+bnorm1*bux1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*bux2(ng2x(j),ng2y(j),k)+bnorm3*bux3(ng3x(j),ng3y(j),k)
      buy=bnorm0*buy0+bnorm1*buy1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*buy2(ng2x(j),ng2y(j),k)+bnorm3*buy3(ng3x(j),ng3y(j),k)
      buz=bnorm0*buz0+bnorm1*buz1(ng1x(j),ng1y(j),k)
     ,  +bnorm2*buz2(ng2x(j),ng2y(j),k)+bnorm3*buz3(ng3x(j),ng3y(j),k)
c     Compute B field components downstream of shock in the plasma frame
c       by transforming the shock normal to the upstream plasma
c       frame, then compressing the component of B parallel to the shock
      bux=brand*bux-bord*bpsin*sinph(j)
      buy=brand*buy+bord*bpsin*cosph(j)
      buz=brand*buz+bord*bpcos
      bup=sqrt(bux*bux+buy*buy+buz*buz)
c     Calculate upstream B field components parallel + perpendicular to shock normal
c     Unit vector of shock normal at current position
      anx=cosz*cosph(j)
      any=cosz*sinph(j)
      anz=sinz
c     Velocity vector of laminar component of pre-shock flow
      betupx=betaup*cosph(j)*sinpsi(j)
      betupy=betaup*sinph(j)*sinpsi(j)
      betupz=betaup*cospsi(j)
c     Velocity vector of the turbulent component of pre-shocked plasma
c     Add velocities of nested zones relativistically
      vvx=betat*vnorm0*vtx0
      vvy=betat*vnorm0*vty0
      vvz=betat*vnorm0*vtz0
      uux=betat*vnorm1*vtx1
      uuy=betat*vnorm1*vty1
      uuz=betat*vnorm1*vtz1
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      vvx1=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      vvy1=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      vvz1=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=1
c      write(3,9059)i,j,numv,vvx1,vvy1,vvz1,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
      vvx=vvx1
      vvy=vvy1
      vvz=vvz1
      uux=betat*vnorm2*vtx2
      uuy=betat*vnorm2*vty2
      uuz=betat*vnorm2*vtz2
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      vvx1=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      vvy1=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      vvz1=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=2
c      write(3,9059)i,j,numv,vvx1,vvy1,vvz1,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
      vvx=vvx1
      vvy=vvy1
      vvz=vvz1
      uux=betat*vnorm3*vtx3
      uuy=betat*vnorm3*vty3
      uuz=betat*vnorm3*vtz3
      denomv=1.0d0+vvx*uux+vvy*uuy+vvz*uuz
      gvv=1.0d0/dsqrt(1.0d0-(vvx*vvx+vvy*vvy+vvz*vvz))
      betatx=(vvx+uux/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvx)/denomv
      betaty=(vvy+uuy/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvy)/denomv
      betatz=(vvz+uuz/gvv+gvv/(1.0d0+gvv)*(denomv-1.0d0)*vvz)/denomv
c      numv=3
c      write(3,9059)i,j,numv,betatx,betaty,betatz,vvx,vvy,vvz,uux,
c     ,   uuy,uuz,denomv,gvv
c     Velocity vector of the pre-shock plasma including turbulent component
      dotprd=betupx*betatx+betupy*betaty+betupz*betatz
      btparx=dotprd*betupx/betaup**2
      btpary=dotprd*betupy/betaup**2
      btparz=dotprd*betupz/betaup**2
      btprpx=betatx-btparx
      btprpy=betaty-btpary
      btprpz=betatz-btparz
      betaux(j)=(betupx+btparx+btprpx/gamup)/(1.0d0+dotprd)
      betauy(j)=(betupy+btpary+btprpy/gamup)/(1.0d0+dotprd)
      betauz(j)=(betupz+btparz+btprpz/gamup)/(1.0d0+dotprd)
      betau(j)=dsqrt(betaux(j)**2+betauy(j)**2+betauz(j)**2)
      gammau(j)=1.0d0/dsqrt(1.0d0-betau(j)**2)
c     Unit vector of shock front at current position
      sx=-sinz*cosph(j)
      sy=-sinz*sinph(j)
      sz=cosz
c     Velocity vector downstream of shock + compression ratio of shock
      if(j.gt.1)call vdcalc(betaux(j),betauy(j),betauz(j),sx,sy,sz,
     ,  betadx(i,j),betady(i,j),betadz(i,j),betad(i,j),gammad(i,j),eta)
      betacs=betadx(i,j)*slos+betadz(i,j)*clos
      delta(i,j)=1.0d0/(gammad(i,j)*(1.0d0-betacs))
c     Upstream velocity component parallel to shock normal
      velpar=betaux(j)*anx+betauy(j)*any+betauz(j)*anz
c     Calculate magnetic field components (in plasma frame) downstream of shock
c      parallel and perpendicular to shock normal
      call bcalc(betaux(j),betauy(j),betauz(j),betau(j),gammau(j),
     ,  anx,any,anz,bux,buy,buz,bparx,bpary,bparz,bprpx,bprpy,bprpz,
     ,  bpar,bprp)
      bx(i,j)=bparx+eta*bprpx
      by(i,j)=bpary+eta*bprpy
      bz(i,j)=bparz+eta*bprpz
      n0mean=spsd(idelay+ip0)*(anorm0*an00+anorm1*an01+anorm2*an02+
     ,  anorm3*an03)
      n0(i,j)=eta*n0mean
c     Calculate the initial maximum electron energy in each cell
c     Next line relates this to direction of B field relative to shock
c      gmax0(i,j)=gmaxmx*(bprp*bprp/(bpar*bpar+bprp*bprp))
c      if(gmax0(i,j).lt.gmaxmn)gmax0(i,j)=gmaxmn
      gmax0(i,j)=gmaxmx
c     Maximum electron energy reached only if shock is subluminal
c       in de Hoffman-Teller frame
      htlim=1.0/(gammau(j)*velpar)
c     Ratio of B field transverse and parallel to shock normal
      bfang=bprp/bpar
      if(bfang.gt.htlim)gmax0(i,j)=gmaxmn
c     Next 3 lines assume a power-law distribution of gmax0, unrelated
c         to direction of B field
c     See http://mathworld.wolfram.com/RandomNumber.html for choosing
c      random numbers from a power-law distribution
c      xrand=rand(0)
c      if(pexp.eq.0.0)gmax0(i,j)=gmaxmx
c      if(pexp.lt.0.0)gmax0(i,j)=((gmaxmx**pexp-gmaxmn**pexp)*xrand+
c     ,  gmaxmn**pexp)**(1.0/pexp)
      do 105 ig=1,43
  105 if(gmax0(i,j).le.gcnt(ig+1).and.gmax0(i,j).ge.gcnt(ig))
     , igcnt(ig)=igcnt(ig)+1
c
c     *** Time loop resumes here ***
c
  110 continue
c     Calculate component of magnetic field that is perpendicular to
c       the aberrated line of sight in the plasma frame
c     Line-of-sight vector in plasma frame
      slx=slos
      sly=0.0
      slz=clos
      dum1=0.0
      dum2=0.0
      dum3=0.0
      dum4=0.0
      dum5=0.0
      dum6=0.0
      dum7=0.0
      call bcalc(betadx(i,j),betady(i,j),betadz(i,j),betad(i,j),
     ,  gammad(i,j),slx,sly,slz,bx(i,j),by(i,j),bz(i,j),dum1,dum2,
     ,  dum3,dum4,dum5,dum6,dum7,bperp(j))
      bfield(j)=sqrt(bx(i,j)**2+by(i,j)**2+bz(i,j)**2)
      do 111 inu=1,68
      emisco(inu)=0.0
      ecflux(inu)=0.0
      flcomp(i,j,inu)=0.0
      flsync(i,j,inu)=0.0
      flec(i,j,inu)=0.0
      flssc(i,j,inu)=0.0
  111 flux(i,j,inu)=0.0
      bdx=betadx(i,j)
      bdy=betady(i,j)
      bdz=betadz(i,j)
      delt=(i-1)*zsize*yr*3.26/(gammad(i,j)*betad(i,j))
      bperpp=bperp(j)
      bfld=bfield(j)
c     Determine velocity of the cell plasma relative to the MD plasma
      bmparx=betamd*bdx*bdz/betad(i,j)**2
      bmpary=betamd*bdy*bdz/betad(i,j)**2
      bmparz=betamd*bdz*bdz/betad(i,j)**2
      bmprpx=-bmparx
      bmprpy=-bmpary
      bmprpz=betamd-bmparz
      betamx(j)=(bdx-bmparx-bmprpx/gammad(i,j))/(1.0d0-betamd*bdz)
      betamy(j)=(bdy-bmpary-bmprpy/gammad(i,j))/(1.0d0-betamd*bdz)
      betamz(j)=(bdz-bmparz-bmprpz/gammad(i,j))/(1.0d0-betamd*bdz)
      betamr(j)=dsqrt(betamx(j)**2+betamy(j)**2+betamz(j)**2)
      gamamr(j)=1.0d0/dsqrt(1.0d0-betamr(j)**2)
      zcl=zcell(i,j)
      zrel=zshock-0.5*zsize-zcl
      dmd(j)=sqrt(zrel**2+rcell(j)**2)
c     Calculate angle between MD seed photon and scattered l.o.s. photon
      call scatcs(bdx,bdy,bdz,clos,0.0d0,slos,xcell(j),ycell(j),
     ,   -zrel,cscat)
c     Determine Doppler factor of the MD emission in the cell's frame
      dmdp=sqrt(rcell(j)**2+(zrel/gamamr(j))**2)
      cosmd=(betamx(j)*xcell(j)+betamy(j)*ycell(j)+
     ,  betamz(j)*zrel/gamamr(j))/(betamr(j)*dmdp)
      if((cosmd**2).lt.0.00001)go to 4143
      tanmd=dsqrt(1.0d0/cosmd**2-1.0d0)
      angm=datan(tanmd)
      go to 4144
 4143 angm=pio2
 4144 tanmd=dtan(2.0*datan(dtan(0.5d0*angm)/
     ,  (gamamr(j)*(1.0d0+betamr(j)))))
      cosmd=1.0d0/dsqrt(1.0d0+tanmd**2)
      deltmd(j)=1.0d0/(gamamr(j)*(1.0d0-betamr(j)*cosmd))
c     Determine time step of Mach disk seed photons from light-travel delay
      delcor=deltmd(j)/dopref
      md=mdmd-0.5-(((zrf-zcl)*clos+(xrf-xcell(j))*slos)+
     ,   dmd(j))*betadd/zsize
      do 4145 inu=1,68
      fmdall(inu)=0.0
      syseed(inu)=0.0
      scseed(inu)=0.0
      tauxmd(inu)=0.0
      ajofnu(inu)=0.0
      ssabs(inu)=0.0
      sstau(inu)=0.0
      srcfn(inu)=0.0
 4145 ssseed(inu)=0.0
c     Mach disk's B field in frame of cell plasma transverse to cell's
c       line of sight to Mach disk (for synchrotron calculation)
      sx=rcell(j)*cosph(j)
      sy=rcell(j)*sinph(j)
      sz=zrel
      dum1=1.0
      dum2=1.0
      dum3=1.0
      dum4=1.0
      dum5=1.0
      dum6=1.0
      dum7=1.0
      call bcalc(betamx(j),betamy(j),betamz(j),betamr(j),gamamr(j),
     ,  sx,sy,sz,bmdx(md),bmdy(md),bmdz(md),dum1,dum2,dum3,dum4,
     ,  dum5,dum6,dum7,bmperp)
c     Apply as a correction factor to previous estimate of B_perpendicular
      bpcorr=bmperp/bmdtot(md)
c
c     Calculate emission from various processes
c
c     Calculate synchrotron flux from Mach disk in frame of cell plasma
c
      do 1146 inu=1,40
c     Synchrotron mean intensity for inverse Compton calculation
      ssabs(inu)=0.0
      if(fsynmd(md,inu).eq.0.0.or.alphmd(md,inu).gt.9.0)
     ,   go to 1145
      if(absorb(md,inu).gt.1.0e-35)ssabs(inu)=
     ,  absorb(md,inu)*bpcorr**(0.5*(sen+2.0))*delcor**(abexmd(md,inu))
      path=2.0*rsize*svmd/(tanmd*cosmd)
      if(path.gt.zsmd(md))path=zsmd(md)
      sstau(inu)=ssabs(inu)*parsec*path
      if(fsynmd(md,inu).lt.1.0e-35)go to 1546
      ajofnu(inu)=fsynmd(md,inu)*bpcorr**(1.0+alphmd(md,inu))*
     ,  delcor**(2.0+alphmd(md,inu))
      ajofnu(inu)=ajofnu(inu)*zsmd(md)*(emfact*parsec)
 1546 syseed(inu)=ajofnu(inu)
      srcfn(inu)=ajofnu(inu)
      if(sstau(inu).gt.0.01)srcfn(inu)=ajofnu(inu)/sstau(inu)
      if(sstau(inu).gt.5.0)syseed(inu)=srcfn(inu)
      if(sstau(inu).gt.0.01.and.sstau(inu).le.5.0)syseed(inu)=
     ,  srcfn(inu)*(1.0-exp(-sstau(inu)))
      syseed(inu)=syseed(inu)*pi*rsize**2*vmd/(dmd(j)**2)
      tauexp(inu)=0.0
      abssss=1.02e4*(sen+2.0)*akapnu(nu(inu))*bperp(j)*(parsec/nu(inu))
      if(abssss.gt.(1.0e-35*nu(inu)))tauexp(inu)=abssss/nu(inu)*dmd(j)
      tauxmd(inu)=tauexp(inu)
      if(tauexp(inu).gt.15.0)syseed(inu)=0.0
      if(tauexp(inu).le.15.0)syseed(inu)=syseed(inu)/exp(tauexp(inu))
 1145 scseed(inu)=0.0
      if(inu.lt.16)fmdall(inu)=fmdall(inu)+syseed(inu)
 1146 continue
c
c     Calculate inverse Compton flux from Mach disk in frame of cell plasma
c
      do 1147 inu=16,68
c     Inverse Compton mean intensity from Mach disk for 2nd-order
c       inverse Compton calculation in cells
      scseed(inu)=fsscmd(md,inu)*
     ,  delcor**(2.0+alfmdc(md,inu))
      scseed(inu)=scseed(inu)*zsmd(md)*pi*rsize**2*
     ,  (amjy*parsec)*vmd/dmd(j)**2
      if(tauxmd(inu).gt.15.0)scseed(inu)=0.0
      if(tauxmd(inu).le.15.0)scseed(inu)=scseed(inu)/
     ,exp(tauxmd(inu))
      fmdall(inu)=fmdall(inu)+syseed(inu)+scseed(inu)
 1147 continue
      nuhi=2
      do 4148 inu=1,68
      ssseed(inu)=fmdall(inu)
 4148 if(fmdall(inu).gt.0.0)nuhi=inu
c     Calculate seed photon energy density energy loss calculation
      usdmd=0.0
      do 1149 inu=2,nuhi
c     Only calculate up to Klein-Nishina limit of gmin
      epslon=6.63e-27*snu(inu)/(gmin*emc2)
      if(epslon.gt.1.0)go to 1150
      if(fmdall(inu).le.1.0e-38.or.fmdall(inu-1).le.1.0e-38)go to 1148
      aaa=dlog10(fmdall(inu-1)/fmdall(inu))/alog10(snu(inu)/snu(inu-1))
      if(aaa.gt.0.999.and.aaa.lt.1.0001)go to 4147
      usdmd=usdmd+0.5/c/(1.0-aaa)*fmdall(inu-1)*snu(inu-1)*
     ,  ((snu(inu)/snu(inu-1))**(1.0-aaa)-1.0)
      go to 4149
 4147 usdmd=usdmd+0.5*fmdall(inu-1)*snu(inu-1)*
     ,  alog(snu(inu)/snu(inu-1))/c
 4149 continue
      if(usdmd.ne.usdmd)write(6,9994)i,j,aaa,snu(inu-1),snu(inu),
     ,  fmdall(inu-1),fmdall(inu)
 1148 continue
 1149 continue
 1150 continue
      id=zcell(i,j)/zsize+1
      if(id.lt.1)id=1
c     Skip flux calculation for cell if gamma_max is too low to emit at lowest frequency
c     Ratio of energy density of photons emitted by hot dust + Mach disk to
c       energy density of the magnetic field
      zdist=zdist0+(id-1)*zsize
c     Calculate min & max angles of dust torus in plasma frame
      dphi1=asin(dtrad/sqrt(zdist**2+dtdist**2))
      dphi2=atan(dtdist/zdist)
      dth1=dphi2-dphi1
      dth2=dphi2+dphi1
      csth1=dcos(dth1)
      csth2=dcos(dth2)
      dcsth1=-(csth1-betad(i,j))/(1.0d0-csth1*betad(i,j))
      dcsth2=-(csth2-betad(i,j))/(1.0d0-csth2*betad(i,j))
      csang=0.5*(dcsth1+dcsth2)
c     Doppler factor of dust torus emission in frame of cells, used
c       to estimate frequency of peak intensity in plasma frame
      tdel=gammad(i,j)*(1.0d0-betad(i,j)*0.5d0*(csth1+csth2))
c     Correct energy density of seed photons from dust for Doppler factor
c       that includes turbulent component of velocity; approximates that
c       motion is directed parallel to jet axis; also include first-order
c       dependence on scattering angle in Klein-Nishina cross-section
      useedr=useed(id)*(tdel/tdelr(id))**2/(1.0d0+betad(i,j)*csang)
      ustob=8.0*pi*(useedr+usdmd)/(bfield(j))**2
c     Calculate the maximum electron energy from gmax0 and solution to equation
c       d gamma/dt = - cc2*(b**2+8*pi*useed)*gamma**2
      tlfact=cc2*bfield(j)**2*(1.0+ustob)
      cellno=i
      tlavg=tlf1(j)/cellno
      tlf(i)=tlfact
      do 112 l=istart,i
  112 tlavg=tlavg+tlf(l)/cellno
c     Calculate energy distribution for the cell
      id=zcell(i,j)/zsize+1
      if(id.lt.1)id=1
      if(id.gt.(icells))write(6,9992)i,j,id,zcell(i,j),zshock,
     ,  zsize
      delt=zsize*yr*3.26/(gammad(i,j)*betad(i,j))
      glow=gammin(i-1,j)/(1.0+tlfact*gammin(i-1,j)*delt)
      gmrat=0.99*gammax(i-1,j)/gammin(i-1,j)
      gmratl=dlog10(gmrat)/32.0
      gmratm=dlog10(gammin(i-1,j)/glow)/11.0
      t2max=i*delt
      gamb=gammax(i-1,j)/(1.0+tlfact*gammax(i-1,j)*delt)
      ibreak=0
      if(gamb.le.glow.or.gamb.ge.gammax(i-1,j))ibreak=1
      do 190 ie=1,44
      egam(i,j,ie)=gammin(i-1,j)*10.0**(gmratl*(ie-12))
      if(ie.lt.12)egam(i,j,ie)=glow*10.0**(gmratm*(ie-1))
      if(ie.eq.1)egam(i,j,ie)=glow+0.2*(glow*10.0**(gmratm)-glow)
      if(ibreak.eq.1.or.egam(i,j,ie).lt.gamb)go to 189
      egam(i,j,ie)=gamb
      ibreak=1
  189 ggam(ie)=egam(i,j,ie)
      enofe(i,j,ie)=0.0
      t1=(i-1)*delt
      tloss=(gmax0(i,j)-ggam(ie))/(tlavg*ggam(ie)*gmax0(i,j))
      if(tloss.le.t1.or.enofe((i-1),j,ie).lt.1.0e-35)go to 188
      t2=dmin1(tloss,t2max)
      tlmin=t2max-(gmin-ggam(ie))/(tlavg*gmin*ggam(ie))
      if(ie.lt.12)t1=tlmin
      eterm1=1.0d0-ggam(ie)*t1*tlfact
      eterm2=1.0d0-ggam(ie)*t2*tlfact
      if(eterm1.lt.0.0.or.eterm2.lt.0.0)go to 188
      enofe(i,j,ie)=n0(i,j)/((sen-1.0)*tlavg*ggam(ie)**(sen+1.0))*
     ,  (eterm1**(sen-1.0)-eterm2**(sen-1.0))
c     Divide by cell crossing time since integral is over time
      if(enofe(i,j,ie).lt.1.0e-30)enofe(i,j,ie)=0.0
      enofe(i,j,ie)=enofe(i,j,ie)/delt
      if(enofe(i,j,ie).ne.enofe(i,j,ie))write(6,9994)j,ie,ggam(ie),
     ,  enofe(i,j,ie),tloss,tlfact,t1,t2,t2max,gmax0(i,j),
     ,  useedr,usdmd,bfield(j)
  188 edist(ie)=enofe(i,j,ie)
  190 continue
      gammax(i,j)=gammax(i-1,j)/(1.0d0+tlfact*delt*gammax(i-1,j))
      gammin(i,j)=gammin(i-1,j)/(1.0d0+tlfact*delt*gammin(i-1,j))
      if(gammin(i,j).lt.1.0)gammin(i,j)=1.0
      if(gammax(i,j).lt.2.0)gammax(i,j)=2.0
      do 4191 inu=1,68
      flsync(i,j,inu)=0.0
      flcomp(i,j,inu)=0.0
      flec(i,j,inu)=0.0
      flssc(i,j,inu)=0.0
      flux(i,j,inu)=0.0
      fpol(i,j,inu)=0.0
      pq(i,j,inu)=0.0
      pu(i,j,inu)=0.0
 4191 continue
      if(gammax(i,j).le.2.0)go to 201
      emold=0.0
      emeold=0.0
      emsold=0.0
c     calculate flux in mJy from cell
      do 191 inu=1,22
      dustnu(inu)=dustf(id,inu)
      hnukt=4.8e-11*dustnu(inu)/tdust
      hnuktr=4.8e-11*dustnu(inu)*(tdel/tdelr(id))/tdust
      dusti(inu)=dustii(id,inu)
 9120 format('dust ',3i6,1p7e12.4)
      if(hnuktr.gt.60.0)go to 5193
      if(dusti(inu).gt.1.0e-30)dusti(inu)=dusti(inu)*
     ,  (tdel/tdelr(id))**3*(1.0-exp(hnukt))/(1.0-exp(hnuktr))
      go to 5194
 5193 dusti(inu)=0.0
 5194 continue
  191 continue
      betd=betad(i,j)
      gamd=gammad(i,j)
      call polcalc(bfield(j),bx(i,j),by(i,j),bz(i,j),clos,slos,chipol)
      if(chipol.lt.-pio2)chipol=chipol+pi
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 192 inu=1,68
      rnu(inu)=nu(inu)*zred1/delta(i,j)
      specin(inu)=0.0001
      emisco(inu)=0.0
      fsync2(inu)=0.0
      ssabs(inu)=0.0
      ecflux(inu)=0.0
      sscflx(inu)=0.0
      if(inu.gt.44)go to 193
      emisco(inu)=ajnu(rnu(inu))*bperp(j)*delta(i,j)**2
      fsync2(inu)=emisco(inu)*zsize*(emfact*parsec)
      fsnoab(inu)=fsync2(inu)
c     Skip frequencies that are obviously optically thin
      if(inu.gt.15)go to 193
      abfact=(parsec/rnu(inu))/(rnu(inu)*delta(i,j))
      ssssab=1.02e4*(sen+2.0)*akapnu(rnu(inu))*bperp(j)
      if(ssssab.gt.(1.0e-35/abfact))ssabs(inu)=ssssab*abfact*
     ,  zsize*delta(i,j)
c     Attenuate Mach disk emission from synchrotron absorption on way to cell
      ssabsm(inu)=ssabs(inu)*dmd(j)/zsize
      if(ssabsm(inu).ge.15.0)ssseed(inu)=0.0
      if(ssabsm(inu).lt.15.0.and.ssabsm(inu).gt.0.1)ssseed(inu)=
     ,     ssseed(inu)/exp(ssabsm(inu))
c     Return to absorption within cell
c     Use rsize instead of zsize because of aberration
      ssabs(inu)=ssabs(inu)*rsize/zsize
      fsync2(inu)=fsync2(inu)*rsize/zsize
      if(ssabs(inu).gt.0.001)srcfn(inu)=fsync2(inu)/ssabs(inu)
      if(ssabs(inu).gt.5.0)fsync2(inu)=srcfn(inu)
      if(ssabs(inu).gt.0.01.and.ssabs(inu).le.5.0)fsync2(inu)=
     ,  srcfn(inu)*(1.0-exp(-ssabs(inu)))
      tauexp(inu)=(nouter(j)-i)*ssabs(inu)
      if(rcell(j).gt.(0.98*rbound).and.xcell(j).le.0.0)
     ,         tauexp(inu)=0.0
      if(thlos.eq.0.0)tauexp(inu)=(nouter(j)-i)*ssabs(inu)
      if(tauexp(inu).gt.15.0)fsync2(inu)=0.0
      if(tauexp(inu).le.15.0)fsync2(inu)=fsync2(inu)/exp(tauexp(inu))
c     Revert to using zsize since next step divides intensity by zsize
      fsync2(inu)=fsync2(inu)*zsize/rsize
      specin(inu)=alpha
  193 continue
      flsync(i,j,inu)=0.0
      if(fsync2(inu).gt.1.0e-30)flsync(i,j,inu)=fsync2(inu)*
     ,  (volc/zsize)*zred1/(1.0e18*amjy*dgpc**2)*fgeom
      if(flsync(i,j,inu).lt.0.0)flsync(i,j,inu)=0.0
      flux(i,j,inu)=flsync(i,j,inu)
  192 continue
!$omp end parallel do
      fq1=0.98*rnu(1)
      emold=ajnu(fq1)*bperp(j)*delta(i,j)**2*zsize*(emfact*parsec)
      fq1=0.98*nu(1)
      do 194 inu=1,43
      if(emold.gt.0.0.and.fsnoab(inu).gt.0.0)specin(inu)=
     ,  dlog10(emold/fsnoab(inu))/alog10(nu(inu)/fq1)
      if(specin(inu).lt.alpha)specin(inu)=alpha
c     Determine degree of polarization
      poldeg=(specin(inu)+1.0)/(specin(inu)+5.0/3.0)
      if(ssabs(inu).gt.1.0)poldeg=3.0/(12.0*specin(inu)+19.0)
      if(ssabs(inu).gt.1.0)chipol=chipol+pio2
      if(chipol.gt.pi)chipol=chipol-pi
      fpol(i,j,inu)=poldeg*flsync(i,j,inu)
      pq(i,j,inu)=fpol(i,j,inu)*cos(2.0*chipol)
      pu(i,j,inu)=fpol(i,j,inu)*sin(2.0*chipol)
      fq1=nu(inu)
      emold=fsnoab(inu)
  194 continue
      betd=betad(1,1)
      gamd=gammad(1,1)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 195 inu=14,68
      ecflux(inu)=0.0
      sscflx(inu)=0.0
      if(rnu(inu).lt.1.0e14)go to 185
      ecflux(inu)=ecdust(rnu(inu))*3.086*volc*zred1*delta(i,j)**2/
     ,   dgpc**2*fgeom
      sscflx(inu)=ssc(rnu(inu))*3.086*volc*zred1*delta(i,j)**2/
     ,   dgpc**2*fgeom
  185 taupp(inu)=0.0
  195 continue
!$omp end parallel do
      do 186 inu=9,68
      if(nu(inu).lt.1.0e22.or.ipprod.eq.0)go to 198
c     Pair production opacity calculation
c     Expression for anumin includes typical interaction angle
      anumin=(1.24e20/(nu(inu)*zred1))*1.24e20*(2.0*gammad(i,j))**2
      alnumn=alog10(anumin)
      inumin=(alnumn-10.0)*4+1
      if(inumin.gt.40)go to 198
      if(anumin.gt.nu(inumin))anumin=nu(inumin)
      do 196 iinu=inumin,40
      bp=sqrt(1.0-anumin/nu(iinu))
      if(bp.le.0.0)bp=0.001
      xsecp1=1.25e-25*(1.0-bp*bp)*((3.0-bp**4)*
     *  alog((1.0+bp)/(1.0-bp))-2.0*bp*(2.0-bp*bp))
      bp=sqrt(1.0-anumin/nu(iinu+1))
      if(bp.le.0.0)bp=0.001
      xsecp2=1.25e-25*(1.0-bp*bp)*((3.0-bp**4)*
     *  alog((1.0+bp)/(1.0-bp))-2.0*bp*(2.0-bp*bp))
      taupp(inu)=taupp(inu)+0.5*(phots(iinu)*xsecp1+
     ,  phots(iinu+1)*xsecp2)*
     ,  (nu(iinu+1)-nu(iinu))*(nouter(j)-i)*zsize*parsec
  196 continue
      if(taupp(inu).lt.0.1)go to 198
      if(taupp(inu).lt.10.0)go to 197
  187 ecflux(inu)=0.0
      sscflx(inu)=0.0
      go to 198
  197 ecflux(inu)=ecflux(inu)/exp(taupp(inu))
      sscflx(inu)=sscflx(inu)/exp(taupp(inu))
  198 flec(i,j,inu)=ecflux(inu)
      flssc(i,j,inu)=sscflx(inu)
      flcomp(i,j,inu)=ecflux(inu)+sscflx(inu)
      flux(i,j,inu)=flsync(i,j,inu)+flcomp(i,j,inu)
  186 continue
      do 199 inu=15,68
      if(ecflux(inu-1).gt.0.0.and.ecflux(inu).gt.0.0)spxec(inu)=
     ,  dlog10(ecflux(inu-1)/ecflux(inu))/alog10(nu(inu)/nu(inu-1))
      if(sscflx(inu-1).gt.0.0.and.sscflx(inu).gt.0.0)spxssc(inu)=
     ,  dlog10(sscflx(inu-1)/sscflx(inu))/alog10(nu(inu)/nu(inu-1))
  199 continue
  201 continue
  200 continue
 2200 continue
c     End loop over downstream cells
      tflold=0.0
      do 500 inu=1,68
      tflux=0.0
      tsflux=0.0
      tcflux=0.0
      tecfl=0.0
      tsscfl=0.0
      alph=0.0
      qcum=0.0
      ucum=0.0
      iwp=0
      if(itt.eq.2)iwp=1
      if(itt.eq.25.or.itt.eq.50)iwp=1
      if(itt.eq.75.or.itt.eq.100)iwp=1
      if(itt.eq.125.or.itt.eq.150)iwp=1
      if(itt.eq.175.or.itt.eq.200)iwp=1
      if(itt.eq.225.or.itt.eq.250)iwp=1
      if(itt.eq.275.or.itt.eq.300)iwp=1
      if(itt.eq.325.or.itt.eq.350)iwp=1
      if(itt.eq.375.or.itt.eq.400)iwp=1
      if(itt.eq.425.or.itt.eq.450)iwp=1
      if(itt.eq.475.or.itt.eq.500)iwp=1
      if(itt.eq.525.or.itt.eq.550)iwp=1
      if(itt.eq.575.or.itt.eq.600)iwp=1
      if(itt.eq.625.or.itt.eq.650)iwp=1
      if(itt.eq.675.or.itt.eq.700)iwp=1
      if(itt.eq.725.or.itt.eq.750)iwp=1
      if(itt.eq.775.or.itt.eq.800)iwp=1
      if(itt.eq.825.or.itt.eq.850)iwp=1
      if(itt.eq.875.or.itt.eq.900)iwp=1
      if(itt.eq.925.or.itt.eq.950)iwp=1
      if(itt.eq.975.or.itt.eq.1000)iwp=1
      if(itt.eq.1025.or.itt.eq.1050)iwp=1
      if(itt.eq.1075.or.itt.eq.1100)iwp=1
      if(itt.eq.1125.or.itt.eq.1150)iwp=1
      if(itt.eq.1175.or.itt.eq.1200)iwp=1
      if(itt.eq.1225.or.itt.eq.1250)iwp=1
      if(itt.eq.1275.or.itt.eq.1300)iwp=1
      if(itt.eq.1325.or.itt.eq.1350)iwp=1
      if(itt.eq.1375.or.itt.eq.1400)iwp=1
      if(itt.eq.1425.or.itt.eq.1450)iwp=1
      if(itt.eq.1475.or.itt.eq.1500)iwp=1
      if(itt.eq.1525.or.itt.eq.1550)iwp=1
      if(itt.eq.1575.or.itt.eq.1600)iwp=1
      if(itt.eq.1625.or.itt.eq.1650)iwp=1
      if(itt.eq.1675.or.itt.eq.1700)iwp=1
      if(itt.eq.1725.or.itt.eq.1750)iwp=1
      if(itt.eq.1775.or.itt.eq.1800)iwp=1
      if(itt.eq.1825.or.itt.eq.1850)iwp=1
      if(itt.eq.1875.or.itt.eq.1900)iwp=1
      if(itt.eq.1925.or.itt.eq.1950)iwp=1
      if(itt.eq.1975.or.itt.eq.2000)iwp=1
      if(itt.eq.2025.or.itt.eq.2050)iwp=1
      if(itt.eq.2075.or.itt.eq.2100)iwp=1
      if(itt.eq.2125.or.itt.eq.2150)iwp=1
      if(itt.eq.2175.or.itt.eq.2200)iwp=1
      if(itt.eq.2225.or.itt.eq.2250)iwp=1
      if(itt.eq.2275.or.itt.eq.2300)iwp=1
      if(itt.eq.2325.or.itt.eq.2350)iwp=1
      if(itt.eq.2375.or.itt.eq.2400)iwp=1
      if(itt.eq.2425.or.itt.eq.2450)iwp=1
      if(itt.eq.2475.or.itt.eq.2500)iwp=1
      if(itt.eq.2525.or.itt.eq.2550)iwp=1
      if(itt.eq.2575.or.itt.eq.2600)iwp=1
      if(itt.eq.2625.or.itt.eq.2650)iwp=1
      if(itt.eq.2675.or.itt.eq.2700)iwp=1
      if(itt.eq.2725.or.itt.eq.2750)iwp=1
      if(itt.eq.2775.or.itt.eq.2800)iwp=1
      if(itt.eq.2825.or.itt.eq.2850)iwp=1
      if(itt.eq.2875.or.itt.eq.2900)iwp=1
      if(itt.eq.2925.or.itt.eq.2950)iwp=1
      if(itt.eq.2975.or.itt.eq.3000)iwp=1
      if(itt.eq.3025.or.itt.eq.3050)iwp=1
      if(itt.eq.3075.or.itt.eq.3100)iwp=1
      if(itt.eq.3125.or.itt.eq.3150)iwp=1
      if(itt.eq.3175.or.itt.eq.3200)iwp=1
      if(itt.eq.3225.or.itt.eq.3250)iwp=1
      if(itt.eq.3275.or.itt.eq.3300)iwp=1
      if(itt.eq.3325.or.itt.eq.3350)iwp=1
      if(itt.eq.3375.or.itt.eq.3400)iwp=1
      if(itt.eq.3425.or.itt.eq.3450)iwp=1
      if(itt.eq.3475.or.itt.eq.3500)iwp=1
      if(itt.eq.3525.or.itt.eq.3550)iwp=1
      if(itt.eq.3575.or.itt.eq.3600)iwp=1
      if(itt.eq.3625.or.itt.eq.3650)iwp=1
      if(itt.eq.3675.or.itt.eq.3700)iwp=1
      if(itt.eq.3725.or.itt.eq.3750)iwp=1
      if(itt.eq.3775.or.itt.eq.3800)iwp=1
      if(itt.eq.3825.or.itt.eq.3850)iwp=1
      if(itt.eq.3875.or.itt.eq.3900)iwp=1
      if(itt.eq.3925.or.itt.eq.3950)iwp=1
      if(itt.eq.3975.or.itt.eq.4000)iwp=1
      if(itt.eq.4025.or.itt.eq.4050)iwp=1
      if(itt.eq.4075.or.itt.eq.4100)iwp=1
      if(itt.eq.4125.or.itt.eq.4150)iwp=1
      if(itt.eq.4175.or.itt.eq.4200)iwp=1
      if(itt.eq.4225.or.itt.eq.4250)iwp=1
      if(itt.eq.4275.or.itt.eq.4300)iwp=1
      if(itt.eq.4325.or.itt.eq.4350)iwp=1
      if(itt.eq.4375.or.itt.eq.4400)iwp=1
      if(itt.eq.4425.or.itt.eq.4450)iwp=1
      if(itt.eq.4475.or.itt.eq.4500)iwp=1
      if(itt.eq.4525.or.itt.eq.4550)iwp=1
      if(itt.eq.4575.or.itt.eq.4600)iwp=1
      if(itt.eq.4625.or.itt.eq.4650)iwp=1
      if(itt.eq.4675.or.itt.eq.4700)iwp=1
      if(itt.eq.4725.or.itt.eq.4750)iwp=1
      if(itt.eq.4775.or.itt.eq.4800)iwp=1
      if(itt.eq.4825.or.itt.eq.4850)iwp=1
      if(itt.eq.4875.or.itt.eq.4900)iwp=1
      if(itt.eq.4925.or.itt.eq.4950)iwp=1
      if(itt.eq.4975.or.itt.eq.5000)iwp=1
c     Do not make maps if nomap=1 was specified at the beginning
      if(nomap.eq.1)iwp=0
      if(iwp.ne.1.or.inu.ne.nrmap)go to 7300
      write(filnam,"(a7,i4.4,a4)") "temzmap",itt,".txt"
      open (7,iostat=ios, err=9000, file=dpath//filnam,status='new')
      write(filnm2,"(a7,i4.4,a4)") "temzmpo",itt,".txt"
      open (8,iostat=ios, err=9000, file=dpath//filnm2,status='new')
 7300 continue
      do 300 j=2,jcells
      do 298 i=1,imax(j)
      if(inu.gt.43)go to 296
      qcum=qcum+pq(i,j,inu)
      ucum=ucum+pu(i,j,inu)
  296 continue
      tflux=tflux+flux(i,j,inu)
      tsflux=tsflux+flsync(i,j,inu)
      tcflux=tcflux+flcomp(i,j,inu)
      tecfl=tecfl+flec(i,j,inu)
      tsscfl=tsscfl+flssc(i,j,inu)
c     Position of cell on sky, in parsecs
      xobs=(zcell(i,j)*slos+xcell(j)*clos)
      yobs=ycell(j)
      poldg=0.0
      poldgg=0.0
      pang=0.0
      if(flsync(i,j,inu).le.1.0e-6)go to 297
      if(iwp.ne.1)go to 297
      if(inu.ne.nrmap.and.inu.ne.nopmap)go to 297
      poldg=100.0*fpol(i,j,inu)/flsync(i,j,inu)
      pangg=0.5*atan2(pu(i,j,inu),pq(i,j,inu))*rad
      if(inu.eq.nrmap)write(7,9876)i,j,xobs,yobs,
     ,  (0.001*flsync(i,j,inu)),
     ,  (0.001*pq(i,j,inu)),(0.001*pu(i,j,inu)),poldg,pangg,
     ,  tflux,qcum,ucum,betad(i,j),it
      if(inu.eq.nopmap)write(8,9876)i,j,xobs,yobs,
     ,  (0.001*flsync(i,j,inu)),
     ,  (0.001*pq(i,j,inu)),(0.001*pu(i,j,inu)),poldg,pangg,
     ,  tflux,qcum,ucum,betad(i,j),it
  297 continue
  298 continue
c  299 continue
  300 continue
      phots(inu)=0.0
      phalph(inu)=alph+1.0
      poldeg=0.0
      if(tflold.le.0.0.or.tflux.lt.1.0e-35)go to 390
      if(inu.gt.1)alph=alog10(tflold/tflux)/
     ,alog10(nu(inu)/nu(inu-1))
      if(inu.le.43)poldeg=sqrt(qcum**2+ucum**2)/tflux
      if(inu.le.43)pang=0.5*atan2(ucum,qcum)*rad
      tflux=0.001*nu(inu)*tflux
      tsflux=0.001*nu(inu)*tsflux
      tcflux=0.001*nu(inu)*tcflux
      tecfl=0.001*nu(inu)*tecfl
      tsscfl=0.001*nu(inu)*tsscfl
      tsscf2=0.001*nu(inu)*tsscf2
c     Print out SED for every ispecs-th time step
  390 continue
      timeo=itt*dtime
      if(inu.eq.1.and.itt.gt.0)write(3,9990)timeo
c     Write SED to file
      tfsync=1.0e3*tsflux/nu(inu)
      if(itt.gt.0)write(3,9996)nu(inu),tflux,alph,tsflux,tecfl,tsscfl
  391 if(inu.eq.18)tfl18=tflux
      if(inu.eq.19)tfl19=tflux
      if(inu.eq.20)tfl20=tflux
      if(inu.eq.21)tfl21=tflux
      if(inu.eq.22)tfl22=tflux
      if(inu.eq.23)tfl23=tflux
      if(inu.eq.31)tfl31=tflux
      if(inu.eq.32)tfl32=tflux
      if(inu.eq.33)tfl33=tflux
      if(inu.eq.34)tfl34=tflux
      if(inu.eq.33)tfcomx=1.0e6*tcflux/nu(inu)
      if(inu.eq.33)tfec=tecfl
      if(inu.eq.33)tfssc=tsscfl
      if(inu.eq.52)tfl52=tflux
      if(inu.eq.53)tfl53=tflux
      if(inu.eq.54)tfl54=tflux
      if(inu.eq.53)tfcomp=1.0e9*tcflux/nu(inu)
      if(inu.eq.53)tfec=tecfl
      if(inu.eq.53)tfssc=tsscfl
      if(inu.eq.11)tfl11=tflux
      if(inu.eq.6)tfl6=tflux
      tflold=1000.0*tflux/nu(inu)
      if(inu.eq.3)pdeg3=100.0*poldeg
      if(inu.eq.3)pang3=pang
      if(inu.eq.4)pdeg4=100.0*poldeg
      if(inu.eq.4)pang4=pang
      if(inu.eq.6)pdeg6=100.0*poldeg
      if(inu.eq.6)pang6=pang
      if(inu.eq.8)pdeg8=100.0*poldeg
      if(inu.eq.8)pang8=pang
      if(inu.eq.12)pdeg12=100.0*poldeg
      if(inu.eq.12)pang12=pang
      if(inu.eq.16)pdeg16=100.0*poldeg
      if(inu.eq.16)pang16=pang
      if(inu.eq.19)pdeg19=100.0*poldeg
      if(inu.eq.19)pang19=pang
      if(inu.eq.20)pdeg20=100.0*poldeg
      if(inu.eq.20)pang20=pang
      if(inu.eq.21)pdeg21=100.0*poldeg
      if(inu.eq.21)pang21=pang
      if(inu.eq.22)pdeg22=100.0*poldeg
      if(inu.eq.22)pang22=pang
      if(inu.eq.24)pdeg24=100.0*poldeg
      if(inu.eq.24)pang24=pang
      if(inu.eq.28)pdeg28=100.0*poldeg
      if(inu.eq.28)pang28=pang
      if(inu.eq.29)pdeg29=100.0*poldeg
      if(inu.eq.29)pang29=pang
      if(inu.eq.31)pdeg31=100.0*poldeg
      if(inu.eq.31)pang31=pang
      if(inu.eq.32)pdeg32=100.0*poldeg
      if(inu.eq.32)pang32=pang
      if(inu.eq.33)pdeg33=100.0*poldeg
      if(inu.eq.33)pang33=pang
  500 continue
      close(7)
      close(8)
      alp19a=0.0
      alp19b=0.0
      alp20a=0.0
      alp20b=0.0
      alp21a=0.0
      alp21b=0.0
      alp22a=0.0
      alp22b=0.0
      alp33a=0.0
      alp33b=0.0
      alp53a=0.0
      alp53b=0.0
      if(tfl18.gt.0.0.and.tfl19.gt.0.0)alp19a=alog10(tfl18/tfl19)/
     ,alog10(nu(19)/nu(18))+1.0
      if(tfl19.gt.0.0.and.tfl20.gt.0.0)alp19b=alog10(tfl19/tfl20)/
     ,alog10(nu(20)/nu(19))+1.0
      alp20a=alp19b
      if(tfl20.gt.0.0.and.tfl21.gt.0.0)alp20b=alog10(tfl20/tfl21)/
     ,alog10(nu(21)/nu(20))+1.0
      alp21a=alp20b
      if(tfl21.gt.0.0.and.tfl22.gt.0.0)alp21b=alog10(tfl21/tfl22)/
     ,alog10(nu(22)/nu(21))+1.0
      alp22a=alp21b
      if(tfl22.gt.0.0.and.tfl23.gt.0.0)alp22b=alog10(tfl22/tfl23)/
     ,alog10(nu(23)/nu(22))+1.0
      if(tfl32.gt.0.0.and.tfl33.gt.0.0)alp33a=alog10(tfl32/tfl33)/
     ,alog10(nu(33)/nu(32))+1.0
      if(tfl33.gt.0.0.and.tfl34.gt.0.0)alp33b=alog10(tfl33/tfl34)/
     ,alog10(nu(34)/nu(33))+1.0
      if(tfl52.gt.0.0.and.tfl53.gt.0.0)alp53a=alog10(tfl52/tfl53)/
     ,alog10(nu(53)/nu(52))+1.0
      if(tfl53.gt.0.0.and.tfl54.gt.0.0)alp53b=alog10(tfl53/tfl54)/
     ,alog10(nu(54)/nu(53))+1.0
      alph19=0.5*(alp19a+alp19b)
      alph20=0.5*(alp20a+alp20b)
      alph21=0.5*(alp21a+alp21b)
      alph22=0.5*(alp22a+alp22b)
      alph33=0.5*(alp33a+alp33b)
      alph53=0.5*(alp53a+alp53b)
      timeo=itt*dtime
c     Write light curve points to file
      if(itt.gt.0)write(4,9991)itt,timeo,nu(6),tfl6,nu(20),tfl20,
     ,  alph20,nu(21),tfl21,alph21,nu(22),tfl22,alph22,nu(33),
     ,  tfl33,alph33,nu(53),tfl53,alph53
c     Write selected polarization data to file
      if(itt.gt.0)write(5,9988)itt,timeo,nu(3),pdeg3,pang3,nu(4),pdeg4,
     ,  pang4,nu(6),pdeg6,pang6,nu(19),pdeg19,pang19,nu(20),pdeg20,
     ,  pang20,nu(21),pdeg21,pang21,nu(22),pdeg22,pang22,nu(28),pdeg28,
     ,  pang28,nu(29),pdeg29,pang29,nu(32),pdeg32,pang32,
     ,  nu(33),pdeg33,pang33
c     Set up next time step by shifting physical conditions 1 slice down jet
      if(itt.eq.itlast)go to 9000
      do 598 j=2,jcells
      do 598 i=imax(j)+1,istart,-1
      gmax0(i,j)=gmax0(i-1,j)
      bx(i,j)=bx(i-1,j)
      by(i,j)=by(i-1,j)
      bz(i,j)=bz(i-1,j)
      n0(i,j)=n0(i-1,j)
      betadx(i,j)=betadx(i-1,j)
      betady(i,j)=betady(i-1,j)
      betadz(i,j)=betadz(i-1,j)
      betad(i,j)=betad(i-1,j)
      gammad(i,j)=gammad(i-1,j)
      delta(i,j)=delta(i-1,j)
  598 continue
      do 1995 kk=icells,2,-1
      do 1995 ng1=1,8
      do 1995 ng2=1,8
      bux1(ng1,ng2,kk)=bux1(ng1,ng2,(kk-1))
      buy1(ng1,ng2,kk)=buy1(ng1,ng2,(kk-1))
      buz1(ng1,ng2,kk)=buz1(ng1,ng2,(kk-1))
 1995 continue
      do 1996 kk=icells,2,-1
      do 1996 ng1=1,4
      do 1996 ng2=1,4
      bux2(ng1,ng2,kk)=bux2(ng1,ng2,(kk-1))
      buy2(ng1,ng2,kk)=buy2(ng1,ng2,(kk-1))
      buz2(ng1,ng2,kk)=buz2(ng1,ng2,(kk-1))
 1996 continue
      do 1997 kk=icells,2,-1
      do 1997 ng1=1,2
      do 1997 ng2=1,2
      bux3(ng1,ng2,kk)=bux3(ng1,ng2,(kk-1))
      buy3(ng1,ng2,kk)=buy3(ng1,ng2,(kk-1))
      buz3(ng1,ng2,kk)=buz3(ng1,ng2,(kk-1))
 1997 continue
      ispec=ispec-1
      if(ispec.eq.0)ispec=ispecs+1
c     Move cells in time array to make room for next time step
      do 599 inu=1,68
      do 599 md=1,(mdmax-1)
      alphmd(md,inu)=alphmd(md+1,inu)
      alfmdc(md,inu)=alfmdc(md+1,inu)
      fsynmd(md,inu)=fsynmd(md+1,inu)
      fsscmd(md,inu)=fsscmd(md+1,inu)
      absorb(md,inu)=absorb(md+1,inu)
      abexmd(md,inu)=abexmd(md+1,inu)
  599 continue
c      if(it.eq.1)go to 9000
      go to 9
 6666 format('# freq(Hz)  flux density (Jy Hz) spectral index')
 6667 format('#'/'#   time(d)      freq(Hz)  F(Jy Hz)  freq  ',
     ,   '    F(Jy Hz) ',
     ,   ' alpha     Freq      F(muJy)    alpha  ',
     ,   '  freq(Hz)  F(Jy Hz)   alpha      Freq    F(Jy Hz) ',
     ,      '   alpha    freq   F(Jy Hz)   alpha')
 6668 format('#'/'#     time(d)  freq(Hz)   p(%)   chi(deg) ',
     ,   '  freq      p      chi     freq       p      chi  ',
     ,   '   freq     p      chi     freq        p        chi  ',
     ,   '   freq     p      chi     freq       p      chi  ',
     ,   '   freq     p      chi     freq       p      chi  ',
     ,   '   freq     p      chi     freq       p      chi')
 7777 format(i4,1x,i4,2f12.4,2f8.2,2x,f10.4)
 8888 format(i5,f9.2)
 8889 format(/)
 9059 format(2i5,i3,11f8.4)
 9111 format(a10/i3,i8,i3/f6.3/f7.3/f4.2/f4.2/f4.2/f5.3/f3.1/e7.1/f5.3/
     ,f7.1/f5.1/f6.1/d9.5/d9.5/d6.1/d6.1/d6.2/f6.1/f4.1/f3.1/f3.1/
     ,f3.1/e7.2/f5.2/f5.2)
 9222 format('idelay out of bounds ',2i12,3i6,1p7e11.3)
 9223 format('*idelay out of bounds ',2i12,3i6,1p7e11.3)
 9232 format(2i4,3i6,1p12e10.2)
 9299 format(i12,2i5,2i6,3i10,1p5e12.4)
 9333 format(4i7,1p12e10.2)
 9623 format('* ',3i5,13f8.4)
 9736 format(6i7,1p10e12.3)
 9875 format('#    i     j    x(mas)     y(mas)   flux(Jy)     Q(Jy)',
     ,   '       U(Jy)      P(%)    EVPA(deg)',
     ,   '  tot flux    qcum      ucum     betad      it')
 9876 format(2i5,1x,1p11e11.3,i5)
 9891 format('Program halted: betadd < sound speed ',1p6e12.5)
 9911 format('seed:  ',2i5,1p11e10.3)
 9988 format(i5,f8.2,2x,11(1pe8.2,1x,0pf7.3,1x,f8.3,1x))
 9989 format(i6,2x,i6,2x,i6,1p16e10.2)
 9990 format(///'Time = ',f8.2,' days'/'  freq(Hz))',2x,
     , 'Ftot(Jy Hz)  sp. index    Fsynch',4x,'     F(EC)',
     , 4x,'F(SSC-MD)')
 9991 format(i5,f8.2,2x,1p18e10.2,1x,i7)
 9992 format('** ',3i6,1p20e12.4)
 9993 format(e12.4,1x,i10,1x,2f8.3)
 9994 format(i5,2x,i5,1p16e12.4)
 9995 format('E dist',2x,1p2e11.4)
 9996 format(1p12e12.4)
 9997 format(3e10.2,2x,1p10e12.4)
 9998 format('# mean, std. dev. pol and EVPA',4f12.5/
     , '# mean, std. dev. pol., layers 1-7',2f12.5)
 9999 format(10f10.3)
      close (4, status='keep')
      close (3, status='keep')
 9000 stop
      end
c
c     Subroutine that calculates the seed photon intensity from emission by hot dust
      function seedph(f)
      common/cdust/csth1,csth2,dcsth1,dcsth2,dsang,tdust
c      common/cvel/bdx,bdy,bdz,gammad,betad
      common/cfreq/freq
      real*8 csth1,csth2,dcsth1,dcsth2,cs1,cs2,csm
      external sdgran
      freq=f
      cs1=csth2
      cs2=csth1
      csm=0.5d0*(cs1+cs2)
      call qg5(cs1,csm,sdgran,ans1)
      call qg5(csm,cs2,sdgran,ans2)
c     Multiply by 2pi, but then divide by 4pi to get mean intensity
      seedph=ans1+ans2
      return
      end
c
      function sdgran(cs)
      common/cdust/csth1,csth2,dcsth1,dcsth2,dsang,tdust
      common/cvel/bdx,bdy,bdz,gammad,betad
      common/cfreq/freq
      real*8 cs,csp,bdx,bdy,bdz,gammad,betad,csth1,csth2,
     ,  dcsth1,dcsth2
      hok=4.80e-11
      sdgran=0.0
      csp=-(cs-betad)/(1.0d0-betad*cs)
      tdel=1.0/(gammad*(1.0d0-betad*csp))
      f=freq/tdel
      expon=hok*f/tdust
      if(expon.lt.0.01)go to 10
      if(expon.gt.10.0)go to 20
      val=(1.33e-26*f)*f/((9.0e20/f)*(exp(expon)-1.0))
      go to 25
   10 val=(2.76e-16*f)*tdust*(f/9.0e20)
      go to 25
   20 if(expon.gt.20.0)go to 26
      val=(1.33e-26*f)*f/((9.0e20/f)*exp(expon))
   25 sdgran=val*tdel*tdel
   26 return
      end
c
c     Subroutine to calculate synchrotron emission coefficient
      function ajnu(anu)
      common/cparm/zred1,bfield,b
      common/cdist/gam,edist
      dimension gam(44),edist(44)
      real*8 gam
      ajnu=0.0
      if(b.lt.1.0e-5)go to 1000
      if(edist(2).lt.1.0e-35.and.edist(40).lt.1.0e-35)go to 1000
      xfact=2.38e-7*anu/b
      x=xfact/(gam(1)*gam(1))
      if(x.gt.0.01)go to 10
      rfact=1.5*x**0.33333333
      go to 20
   10 rfact=0.0
      if(x.lt.25.0)rfact=(1.77205*x**0.293627-
     , 6.51657e-4/x**0.50372)/exp(x)
   20 gran1=edist(1)*rfact
      do 100 i=2,44
      x=xfact/(gam(i)*gam(i))
      if(x.gt.0.01)go to 30
      rfact=1.5*x**0.33333333
      go to 40
   30 rfact=0.0
      if(x.lt.25.0)rfact=(1.77205*x**0.293627-
     , 6.51657e-4/x**0.50372)/exp(x)
   40 gran2=edist(i)*rfact
      if(gran1.eq.0.0.or.gran2.eq.0.0)go to 45
      grat=gam(i)/gam(i-1)
      alg=alog10(grat)
      if(alg.eq.0.0)go to 45
      a=1.0+alog10(gran2/gran1)/alg
      if(a.gt.5.0.or.a.lt.-5.0)go to 45
      if(a.lt.0.01)go to 45
      addit=gran1*(grat**a-1.0)*gam(i-1)/a
      go to 46
   45 addit=0.5*(gran1+gran2)*(gam(i)-gam(i-1))
   46 ajnu=ajnu+addit
      gran1=gran2
  100 continue
 1111 format(1p14e11.4)
 1000 return
      end
c
c     Subroutine to calculate synchrotron absorpion coefficient
      function akapnu(anu)
      common/cparm/zred1,bfield,b
      common/cdist/gam,edist
      dimension gam(44),edist(44)
      real*8 gam
      akapnu=0.0
      if(edist(2).lt.1.0e-35.and.edist(40).lt.1.0e-35)go to 1000
      xfact=2.38e-7*anu/b
      x=xfact/(gam(1)*gam(1))
      akapnu=0.0
      f1=0.0
      if(x.gt.25.0)go to 10
      f1=1.8*x**0.3/exp(x)
   10 gran1=f1*edist(1)/gam(1)
      do 100 i=2,44
      x=xfact/(gam(i)*gam(i))
      gran2=0.0
      if(x.gt.25.0)go to 45
      f2=1.8*x**0.3/exp(x)
      gran2=f2*edist(i)/gam(i)
      if(gran1.eq.0.0.or.gran2.eq.0.0)go to 45
      grat=gam(i)/gam(i-1)
      alg=alog10(grat)
      if(alg.eq.0.0)go to 45
      a=1.0+alog10(gran2/gran1)/alg
      if(a.gt.5.0.or.a.lt.-5.0)go to 45
      if(a.lt.0.01)go to 45
      addit=gran1*(grat**a-1.0)*gam(i-1)/a
      go to 46
   45 addit=0.5*(gran1+gran2)*(gam(i)-gam(i-1))
   46 akapnu=akapnu+addit
      gran1=gran2
  100 continue
      if(akapnu.lt.1.0e-25)akapnu=0.0
 1000 return
      end
c     Subroutine to carry out 10-point Gaussian integration
      subroutine qg10(a,b,func,ss)
      real ss,func
      real*8 a,b
      external func
c     Returns as ss the integral of the function func between a and 
c     b, by ten-point Gauss-Legendre integration:
c     the function is evaluated exactly ten
c     times at interior points in the range of integration.
      integer j
c     The abscissas and weights:
      REAL dx,xm,xr,w(5),x(5) 
      SAVE w,x
      DATA w/.2955242247,.2692667193,.2190863625,.1494513491,
     ,  .0666713443/
      DATA x/.1488743389,.4333953941,.6794095682,.8650633666,
     ,  .9739065285/
      xm=0.5*(b+a)
      xr=0.5*(b-a)
      ss=0
c     Will be twice the average value of the function, since the ten
c     weights (five numbers above each used twice) sum to 2. 
      do 11 j=1,5
      dx=xr*x(j)
      ss=ss+w(j)*(func(xm+dx)+func(xm-dx))
   11 continue
c     Scale the answer to the range of integration.
      ss=xr*ss 
      return
      END
c     Subroutine to carry out 5-point Gaussian integration
      subroutine qg5(a,b,func,ss)
      real*4 ss,func
      external func
c     Returns as ss the integral of the function func between a and 
c     b, by 5-point Gauss-Legendre integration:
c     the function is evaluated exactly ten
c     times at interior points in the range of integration.
      integer j
c     The abscissas and weights:
      real*8 a,b,dx,xm,xr,w(5),x(5) 
      save w,x
      data w/.236926885,.478628670,.568888889,.47862867,
     ,  .236926885/
      data x/-.906179846,-.538469310,0.0,.53846931,.906179846/
      ss=0
      xm=0.5*(a+b)
      xr=0.5*(b-a)
      do 11 j=1,5
      ss=ss+w(j)*func(xm+x(j)*xr)
   11 continue
c     Scale the answer to the range of integration.
      ss=xr*ss 
      return
      end
c     Subroutine to carry out 3-point Gaussian integration     
      subroutine qg3(xl,xu,fct,y)
c     Gaussian quadrature integration using 3 points
      external fct
      a=0.5*(xu+xl)
      b=xu-xl
      y=0.3872983*b
      y=0.2777778*(fct(a+y)+fct(a-y))
      y=b*(y+0.4444444*fct(a))
      return
      end
c     The next subroutine produce a log-normal distributed random value
      subroutine randln(mu,sigma,ans)
      real*4 mu
c     Normalize by calculated mean of 2.94 for input mu=1.00
c     Input sigma = 0.4 then gives calculated sigma ~0.5
      ans=exp(randn(mu,sigma))/2.94
      return
      end
c     Produce normally distributed random value
      function randn(mean,stdev)
      real*4 mean,stdev,randm
      real*8 tmp1,tmp2
      tmp1=rand(0)
      tmp2=rand(0)
      if(tmp1.lt.1.0d-11)tmp1=1.0d-11
      r=(-2.0*dlog(tmp1))**0.5
      theta=2.0d0*3.14159*tmp2
      randn=mean+stdev*r*sin(theta)
      return
      end
c     The next 4 subroutines are from Ritaban Chatterjee to create variations
c       according to an input PSD with slopes beta1 and beta2 at
c       variational frequencies below and above a break frequency nu_break
C**************************************************************
      subroutine psdsim(N,beta1,beta2,nu_break,t_incre1,lc_sim)
C**************************************************************
C N=number of data points in the lc, should be an integer power of 2
C t_incre1=increment in time in each step while resampling the simulated data. 
C This must be larger than the smallest interval between successive data points in the input light curve

      implicit none

      real*8 nu(65536),flux_s1(16384),dat(65536),R(65536)
      real*8 dataim_s1(16384),datareal_s1(16384),amp(16384)
      real*8 flux_s2(16384),dataim_s2(16384),datareal_s2(16384)
      real*8 dataim(16384),datareal(16384),flux_s(16384)
      real*8 fac_norm,fac_norm2
      real*4 ann,beta1,beta2,nu_break,t_incre1,lc_sim(32768)
      real*4 randn
      integer N,j,i,ifile,ik,nn,ISEED1,ISEED2,isign

      fac_norm=1./(N*t_incre1)
      fac_norm2=N**2./(2.*N*t_incre1)

        do j=1,N/2  
            R(j)=randn(0.0,1.0)
            nu(j)= j*fac_norm/86400.                      
            if (nu(j) .le. nu_break) then
              flux_s(j)=nu(j)**beta1
              datareal(j)=dsqrt(fac_norm2*0.5*flux_s(j))*R(j)
            else
              flux_s(j)=(nu(j)**beta2)
     >                      *((nu_break**beta1)/(nu_break**beta2))
              datareal(j)=dsqrt(fac_norm2*0.5*flux_s(j))*R(j)
            endif
        end do
        do j=1,N/2
            R(j)=randn(0.0,1.0)
            dataim(j)=dsqrt(0.5*fac_norm2*flux_s(j))*R(j)
        end do

           dat(1)=0.0
           dat(2)=0.0
           dat(N+1)=datareal(N/2)
           dat(N+2)=0 ! this is at +/-1/T and has to obey A+iB = A-iB => B=0          
           j=1

        do i=3,N-1,2
             dat(i)=datareal(j)
             dat(2*(N+1)-i)=dat(i)
             dat(i+1)=dataim(j)
             dat(2*(N+1)-i+1)=-1.0*dat(i+1)
           j=j+1
        end do
        nn=N
        isign=-1
        call four1(dat,nn,isign)
        j=1
           ann=nn
        do i=1,2*nn,2
           lc_sim(j)=dat(i)/ann !taking the real part
           j=j+1
        end do

      return
      end  
C*****************************************************************
       SUBROUTINE four1(data,nn,isign)
C****************************************************************

       INTEGER isign,nn
       REAL*8 data(2*nn)
c  Replaces data(1:2*nn) by its discrete Fourier transform, 
c  if isign is input as 1; or replaces
c  data(1:2*nn) by nn times its inverse discrete Fourier 
c  transform, if isign is input as -1.
c  data is a complex array of length nn or, equivalently, 
c  a real array of length 2*nn.
c  nn MUST be an integer power of 2 (this is not checked!).
       INTEGER i,istep,j,m,mmax,n
       REAL*8 tempi,tempr
       real*8 theta,wi,wpi,wpr,wr,wtemp !Double precision 
                              !   for the trigonometric recurrences. 
       n=2*nn
       j=1
       do i=1,n,2 !This is the bit-reversal section of the routine.
         if(j.gt.i)then
            tempr=data(j) ! Exchange the two complex numbers.
            tempi=data(j+1)
            data(j)=data(i)
            data(j+1)=data(i+1)
            data(i)=tempr
            data(i+1)=tempi
         endif
       m=nn
1      if ((m.ge.2).and.(j.gt.m)) then
         j=j-m
         m=m/2
         goto 1
       endif
         j=j+m
       end do 
       mmax=2 ! Here begins the Danielson-Lanczos section of 
              ! the routine.
2      if (n.gt.mmax) then ! Outer loop executed log2 nn times.
           istep=2*mmax
           theta=6.28318530717959d0/(isign*mmax) !Initialize for 
                                  ! the trigonometric recurrence.
           wpr=-2.d0*sin(0.5d0*theta)**2
           wpi=sin(theta)
           wr=1.d0
           wi=0.d0
           do m=1,mmax,2 ! Here are the two nested inner loops. 

              do i=m,n,istep
                j=i+mmax        !This is the Danielson-Lanczos formula:
                tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
                tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
                data(j)=data(i)-tempr
                data(j+1)=data(i+1)-tempi
                data(i)=data(i)+tempr
                data(i+1)=data(i+1)+tempi
              end do
            wtemp=wr                    !Trigonometric recurrence.
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
          end do
        mmax=istep
        goto 2                          !Not yet done.
      endif                             !All done.
      return
      END
c
c     polcalc computes the polarization angle chi based on relations
c       derived by Lyutikov et al. 2005, MNRAS, 360, 869
c
      Subroutine polcalc(b,bx,by,bz,clos,slos,chi)
      common/cvel/bdx,bdy,bdz,gammad,betad
      real*8 clos,slos,bdx,bdy,bdz,gammad,betad,term,ndq,q2,
     ,  qx,qy,qz,ex,ey,ez,term1,term2,term3,coschi
      bxh=bx/b
      byh=by/b
      bzh=bz/b
      term1=slos*bxh+clos*bzh
      term2=slos*bdx+clos*bdz
      term3=(bdx*bxh+bdy*byh+bdz*bzh)*gammad/(1.0d0+gammad)
      qx=bxh+(term1-term3)*bdx-term2*bxh
      qy=byh+(term1-term3)*bdy-term2*byh
      qz=bzh+(term1-term3)*bdz-term2*bzh
      q2=qx*qx+qy*qy+qz*qz
      ndq=qx*slos+qz*clos
      term1=-qy*clos
      term2=qx*clos-qz*slos
      term3=qy*slos
      term=dsqrt(q2-ndq*ndq)
      ex=term1/term
      ey=term2/term
      ez=term3/term
      coschi=(-clos*ex+slos*ez)
      if(coschi.gt.1.0d0.and.coschi.lt.1.000001d0)coschi=0.99999999d0
      if(coschi.lt.-1.0d0.and.coschi.gt.-1.000001d0)coschi=-0.99999999d0
      if(ey.gt.1.0d0.and.ey.lt.1.000001d0)ey=0.99999999d0
      if(ey.lt.-1.0d0.and.ey.gt.-1.000001d0)ey=-0.99999999d0
      chi=atan2(ey,coschi)
      return
      end
c
c     vdcalc computes downstream velocity vector; Lorentz transformation of unit
c     vector along shock front follows Lyutikov et al. (2003, ApJ, 597, 998)
c
      subroutine vdcalc(vx,vy,vz,sx,sy,sz,vdx,vdy,vdz,vd,gd,eta)
      common/cang/cosz,sinz
      real*8 vx,vy,vz,sx,sy,sz,v2,v,g,spx,spy,spz,s,g2,
     ,  dotprd,vparx,vpary,vparz,vprpx,vprpy,vprpz,vprp2,vprp,uprp,
     ,  vd,gd,vd2,vdx,vdy,vdz,vfact,vdprpx,vdprpy,vdprpz,vdprp2,
     ,  gdprp,eta,cosz,sinz
      v2=vx*vx+vy*vy+vz*vz
      v=dsqrt(v2)
      g2=1.0d0/(1.0d0-v2)
      g=dsqrt(g2)
      dotprd=vx*sx+vy*sy+vz*sz
c     Calculate components of velocity vector parallel and perpendicular to shock front
      vparx=dotprd*sx
      vpary=dotprd*sy
      vparz=dotprd*sz
      vprpx=vx-vparx
      vprpy=vy-vpary
      vprpz=vz-vparz
      vprp2=vprpx*vprpx+vprpy*vprpy+vprpz*vprpz
      vprp=dsqrt(vprp2)
c     Shock jump condition, from Konigl (1980, Phys. Fluids, 23, 1083)
      vfact=1.0d0
      uprp=g*vprp
c     uprp must exceed the proper sound speed, 1/sqrt(2) for a shock
c     Otherwise, it is a sound wave and the velocity does not change significantly
      if(uprp.gt.0.7071e0)vfact=(1.0d0+1.0d0/(g2*vprp2))/(3.0d0)
      vdprpx=vprpx*vfact
      vdprpy=vprpy*vfact
      vdprpz=vprpz*vfact
      vdprp2=vdprpx*vdprpx+vdprpy*vdprpy+vdprpz*vdprpz
c     Shock compression ratio (downstream density/upstream density)
c       From Cawthorne + Cobb (1990, ApJ, 350, 536)
      eta=dsqrt((8.0d0*v2*sinz*sinz-1.0d0/g2)/(1.0d0-v2*cosz*cosz))*
     ,  v*g*sinz
c     eta must be >= 1; if formula gives < 1 or if argument to dsqrt < 0, it is not a shock
      if(eta.lt.1.0.or.eta.ne.eta)eta=1.0
c      eta=dsqrt(vprp2*(1.0d0-vdprp2)/(vdprp2*(1.0d0-vprp2)))
      vdx=vparx+vdprpx
      vdy=vpary+vdprpy
      vdz=vparz+vdprpz
      vd2=vdx*vdx+vdy*vdy+vdz*vdz
      vd=dsqrt(vd2)
      gd=1.0d0/dsqrt(1.0d0-vd2)
      return
      end
c
c     bcalc computes magnetic field component parallel and perpendicular to shock
c     normal or line of sight; follows Lyutikov et al. (2003, ApJ, 597, 998)
c
      subroutine bcalc(vx,vy,vz,v,g,ax,ay,az,bx,by,bz,
     ,  bparx,bpary,bparz,bprpx,bprpy,bprpz,bpar,bprp)
      real*8 vx,vy,vz,ax,ay,az,v2,v,g,g1,gp1,apx,apy,apz
      common/ccell/xcell,ycell,zcell
      common/crite/iwrite
c     v is velocity in units of c, a is a unit vector normal to the shock front
c     b is magnetic field vector
c     Determine unit vector of shock normal in plasma frame
      v2=v*v
      g1=g-1.0d0
      apx=(1.0d0+g1*vx*vx/v2)*ax+(g1*vx*vy/v2)*ay+(g1*vx*vz/v2)*az
      apy=(g1*vx*vy/v2)*ax+(1.0d0+g1*vy*vy/v2)*ay+(g1*vy*vz/v2)*az
      apz=(g1*vx*vz/v2)*ax+(g1*vy*vz/v2)*ay+(1.0d0+g1*vz*vz/v2)*az
      gp1=dsqrt(apx*apx+apy*apy+apz*apz)
      apx=apx/gp1
      apy=apy/gp1
      apz=apz/gp1
      dotprd=bx*apx+by*apy+bz*apz
c     Calculate components of B field parallel and perpendicular to shock normal
      bparx=dotprd*apx
      bpary=dotprd*apy
      bparz=dotprd*apz
      bpar=bparx*bparx+bpary*bpary+bparz*bparz
      bpar=sqrt(bpar)
      bprpx=bx-bparx
      bprpy=by-bpary
      bprpz=bz-bparz
      bprp=bprpx*bprpx+bprpy*bprpy+bprpz*bprpz
      bprp=sqrt(bprp)
c      if(iwrite.eq.1)write(6,9000)ax,ay,az,apx,apy,apz,bx,by,bz,
c     ,  vx,vy,vz,bpar,bprp,xcell,ycell,zcell
c     , bparx,bpary,bparz,bprpx,bprpy,bprpz,bpar,bprp
c 9000 format(21f8.4)
      return
      end
c
c     vecrot rotates unit vector a by an angle psi about unit vector c along
c     a great circle to create new vector v
c
      subroutine vecrot(ax,ay,az,cx,cy,cz,psi,vx,vy,vz)
      cs=cos(psi)
      s=sin(psi)
      cs1=1.0-cs
      vx=ax*(1.0+cs1*(cx*cx-1.0))-ay*(cz*s-cs1*cx*cy)+
     ,  az*(cy*s+cs1*cx*cz)
      vy=ax*(cz*s+cs1*cx*cy)+ay*(1.0+cs1*(cy*cy-1.0))-
     ,  az*(cx*s-cs1*cy*cz)
      vz=ax*(-cy*s+cs1*cx*cz)+ay*(cx*s+cs1*cy*cz)+
     ,  az*(1.0+cs1*(cz*cz-1.0))
      return
      end
c
      subroutine scatcs(vx,vy,vz,sx,sy,sz,x,y,z,cscat)
      real*8 vx,vy,vz,sx,sy,sz,g,g1,g2,slx,sly,slz,
     ,  smx,smy,smz
c     v is velocity in units of c, s is line-of-sight unit vector
c     (x,y,z) is position of cell relative to Mach disk, cscat is
c     cosine of scattering angle in plasma frame
      v2=vx*vx+vy*vy+vz*vz
      g2=1.0/(1.0d0-v2)
      g=dsqrt(g2)
      g1=g-1.0d0
c     Determine unit vector of l.o.s. in plasma frame (slx,sly,slz)
      slx=(1.0+g1*vx*vx/v2)*sx+g1*vx*vy*sy/v2+g1*vx*vz*sz/v2
      sly=g1*vx*vy*sx/v2+(1.0+g1*vy*vy/v2)*sy+g1*vy*vz*sz/v2
      slz=g1*vx*vz*sx/v2+g1*vy*vz*sy/v2+(1.0+g1*vz*vz/v2)*sz
      s=slx*slx+sly*sly+slz*slz
      s=sqrt(s)
      slx=slx/s
      sly=sly/s
      slz=slz/s
c     Determine unit vector from Mach disk to cell in plasma frame (smx,smy,smz)
      smx=(1.0+g1*vx*vx/v2)*x+g1*vx*vy*y/v2+g1*vx*vz*z/v2
      smy=g1*vx*vy*x/v2+(1.0+g1*vy*vy/v2)*y+g1*vy*vz*z/v2
      smz=g1*vx*vz*x/v2+g1*vy*vz*y/v2+(1.0+g1*vz*vz/v2)*z
      s=smx*smx+smy*smy+smz*smz
      s=sqrt(s)
      smx=smx/s
      smy=smy/s
      smz=smz/s
c     Calculate cosine of scattering angle in plasma frame using head-on
c      approximation (photon is scattering in direction of electron motion)
      cscat=slx*smx+sly*smy+slz*smz
      return
      end
c
c  Subroutine to calculate the Klein-Nishina cross-section; from Dermer & Menon
c     (2009, High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays,
c      and Neutrinos, Princeton U. Press), eq. 6.31
c
      subroutine xseckn(epsi,epsf,g,y,xsec)
      real*8 xsec
      z=epsf/(g*epsi*y)
      xsec=(y+1.0/y-2.0*z+z*z)/(g*epsi)
      return
      end
c  Subroutine to calculate inverse Compton emission from external sources of seed photons
      function ecdust(anuf)
      common/cparm/zred1,bfield,bperp
      common/cdust/csth1,csth2,dcsth1,dcsth2,dsang,tdust
      common/cdist/gam,edist
c      common/cvel/bdx,bdy,bdz,gammad,betad
      common/cseed/dnu,di,cscat
      common/crite/iwrite
      dimension gam(44),edist(44),dnu(22),di(22)
      real*8 gam,dcsth1,dcsth2,dsnth1,
     ,  dsnth2,val1,val2,vala,valb,csth1,csth2,xsecc,xsecx
      ecdust=0.0
      if(edist(2).lt.1.0e-35.and.edist(40).lt.1.0e-35)go to 4000
      gran=0.0
      nhi=22
      nuhi=22
c     Set the highest frequency of seed photons for scattering
      do 2011 nn=10,nhi
      if(dnu(nn).lt.anuf.and.di(nn).gt.1.0e-30) go to 2010
      nuhi=nn-1
      go to 2012
 2010 continue
 2011 continue
 2012 continue
c     Flux will be in mJy, so set x-section as (3e26/8)sigt
      s0=25.0
      homc2=8.099e-21
      ie1=1
      ef=homc2*anuf
      ef1=ef+1.0
      if(ef1.le.gam(1))go to 1
      if(ef1.ge.gam(44).or.(gam(44)/gam(1)).lt.1.01)go to 4000
      do 10 ie=2,44
      if(gam(ie).gt.ef1)go to 11
   10 continue
      ie=44
   11 ie1=ie-1
      g1=ef1
      go to 12
c     In this version, approximate that plasma velocity is along jet axis
c       and that Doppler factor of dust torus is the mean over its solid angle
c       as viewed in the plasma frame
    1 g1=gam(1)
   12 if((gam(ie1+1)/gam(ie1)).gt.1.0002)go to 13
      ie1=ie1+1
      if(ie1.ge.43)go to 4000
      go to 12
   13 vala=s0*edist(ie1)/g1
c     Loop to integrate over electron Lorentz factors (gam)
      do 3000 ie=ie1,43
      g2=gam(ie+1)
      valb=s0*edist(ie+1)/g2
      val1=0.0
      val2=0.0
      gran1=0.0
      addit=0.0
c     Set up loop 1 to integrate over incident photon frequency anui
      anumax=amin1(anuf,dnu(nuhi))
      di1=di(1)
      id=2
      anumin=anuf/((2.0d0*g1)*(g1-ef)*
     ,  (1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1))))
      if(anumin.ge.anumax)go to 601
      if(anumin.gt.dnu(1))go to 2
      anumin=dnu(1)
      go to 5
    2 continue
      do 3 id=2,nuhi
      if(anumin.le.dnu(id))go to 4
    3 continue
      id=nuhi
    4 continue
      a=alog10(di(id-1)/di(id))/alog10(dnu(id-1)/dnu(id))
      di1=di(id-1)*(anumin/dnu(id-1))**a
    5 ide=nuhi
      if(anumax.ge.dnu(nuhi))go to 8
      do 6 idd=id,nuhi
      if(anumax.le.dnu(idd))go to 7
    6 continue
      idd=nuhi
    7 a=alog10(di(idd-1)/di(idd))/alog10(dnu(idd-1)/dnu(idd))
      die=di(idd-1)*(anumax/dnu(idd-1))**a
      ide=idd
      go to 9
    8 die=di(nuhi)
    9 continue
      anui1=anumin
      xx=1.0-ef/g1
      rat=homc2*anui1*g1*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1)))
      call xseckn(rat,ef,g1,xx,xsecc)
      xsecc=xsecc*(g1*homc2*anui1)
      val1=xsecc*(1.0e20/anui1)*(anuf/anui1)*di1*vala
      if(val1.lt.1.0d-40)val1=0.0d0
   25 continue
c     Loop 1 to integrate over incoming photon frequency anui for lower gam value
      do 600 nu=id,ide
      anui2=dnu(nu)
      di2=di(nu)
      if(nu.lt.ide)go to 30
      anui2=anumax
   30 rat=homc2*anui2*g1*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1)))
      call xseckn(rat,ef,g1,xx,xsecc)
      xsecc=xsecc*(g1*homc2*anui2)
      val2=xsecc*(1.0e20/anui2)*(anuf/anui2)*di2*vala
      if(val2.lt.1.0d-40)val2=0.0d0
  525 if(val1.eq.0.0d0.or.val2.eq.0.0d0)go to 845
      test=abs((anui1-anui2)/anui1)
      if(test.lt.0.001)go to 847
      ratnu=anui2/anui1
      a=1.0+dlog10(val2/val1)/alog10(ratnu)
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 845
      addit=val1*(ratnu**a-1.0)*anui1/a
      go to 846
  845 addit=0.5*(val1+val2)*(anui2-anui1)
  846 gran1=gran1+addit
  847 anui1=anui2
      val1=val2
      di1=di2
  600 continue
c     End anui loop 1
  601 continue
      ratg=g2/g1
      ratgl=alog10(ratg)
      valb=s0*edist(ie+1)/g2
      val1=0.0
      val2=0.0
      gran2=0.0
c     Set up loop 2 to integrate over incident photon frequency anui
      anumax=amin1(anuf,dnu(nuhi))
      di1=di(1)
      id=2
      anumin=anuf/((2.0d0*g2)*(g2-ef)*
     ,  (1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2))))
      if(anumin.ge.anumax)go to 1601
      if(anumin.gt.dnu(1))go to 1002
      anumin=dnu(1)
      go to 1005
 1002 continue
      do 1003 id=2,nuhi
      if(anumin.le.dnu(id))go to 1004
 1003 continue
      id=nuhi
 1004 continue
      di1=0.0
      addit=0.0
      if(di(id).lt.1.0e-25.or.di(id-1).lt.1.0e-25)go to 1005
      a=alog10(di(id)/di(id-1))/alog10(dnu(id)/dnu(id-1))
      di1=di(id)*(anumin/dnu(id))**a
 1005 ide=nuhi
      if(anumax.ge.dnu(nuhi))go to 1008
      do 1006 idd=id,nuhi
      if(anumax.le.dnu(idd))go to 1007
 1006 continue
      idd=nuhi
 1007 a=alog10(di(idd-1)/di(idd))/alog10(dnu(idd-1)/dnu(idd))
      die=di(idd-1)*(anumax/dnu(idd-1))**a
      ide=idd
      go to 1009
 1008 die=di(nuhi)
 1009 continue
      anui1=anumin
      xx=1.0-ef/g2
      rat=homc2*anui1*g2*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2)))
      call xseckn(rat,ef,g2,xx,xsecc)
      xsecc=xsecc*(g2*homc2*anui1)
      val1=xsecc*(1.0e20/anui1)*(anuf/anui1)*di1*valb
      if(val1.lt.1.0d-40)val1=0.0d0
 1025 gran2=0.0
c     Loop 2 to integrate over incoming photon frequency anui for upper gam value
      if(id.ge.ide)go to 1601
      do 1600 nu=id,ide
      anui2=dnu(nu)
      di2=di(nu)
      if(nu.lt.ide)go to 1030
      anui2=anumax
      di2=die
 1030 rat=homc2*anui2*g2*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2)))
      call xseckn(rat,ef,g2,xx,xsecc)
      xsecc=xsecc*(g2*homc2*anui2)
      val2=xsecc*(1.0e20/anui2)*(anuf/anui2)*di2*valb
      if(val2.lt.1.0d-40)val2=0.0d0
      if(val1.lt.1.0e-38.or.val2.lt.1.0e-38)go to 1845
      test=abs((anui1-anui2)/anui1)
      if(test.lt.0.001)go to 1847
      ratnu=anui2/anui1
      a=1.0+dlog10(val2/val1)/alog10(ratnu)
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 1845
      addit=val1*(ratnu**a-1.0)*anui1/a
      go to 1846
 1845 addit=0.5*(val1+val2)*(anui2-anui1)
 1846 gran2=gran2+addit
 1847 continue
 6599 continue
      anui1=anui2
      val1=val2
      di1=di2
 1600 continue
c     End anui loop 2
 1601 continue
      if(gran1.le.1.0e-28.or.gran2.le.1.0e-28)go to 2845
      if(ratgl.lt.0.0001)go to 2845
      a=1.0+alog10(gran2/gran1)/ratgl
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 2845
      addit=gran1*(ratg**a-1.0)*g1/a
      go to 2846
 2845 addit=0.5*(gran1+gran2)*(g2-g1)
 2846 gran=gran+addit
      g1=g2
      vala=valb
 3000 continue
c     End gam loop
      if(gran.gt.1.0e-20)ecdust=gran*1.0e-20
 6666 format('In ecdust',i5,2x,1p16e12.3)
 9996 format('a ',1p12e11.3)
 9997 format('b ',1p12e11.3)
 9998 format('c ',1p12e11.3)
 9999 format('d ',1p12e11.3)
 4000 return
      end
c
c  Subroutine to calculate synchrotron self-Compton emission
c    of seed photons from other cells
      function ssc(anuf)
      common/cparm/zred1,bfield,bperp
      common/cdist/gam,edist
      common/cssc/dnu,di,cscat
      common/crite/iwrite
      dimension gam(44),edist(44),dnu(68),di(68)
      real*8 gam,addit,gran,ssca,val1,val2,xsecc,xsecx,di
      ssc=0.0
      if(edist(2).lt.1.0e-35.and.edist(40).lt.1.0e-35)go to 4000
      ssca=0.0
      gran=0.0
      nhi=68
      nuhi=68
c     Set the highest frequency of seed photons for scattering
      do 2011 nn=12,nhi
      if(dnu(nn).lt.anuf.and.di(nn).gt.1.0e-30) go to 2010
      nuhi=nn-1
      go to 2012
 2010 continue
 2011 continue
 2012 continue
c     Flux will be in mJy, so set x-section as (3e26/8)sigt
      s0=25.0
      homc2=8.099e-21
      ie1=1
      ef=homc2*anuf
      ef1=ef+1.0
      if(ef1.le.gam(1))go to 1
      if(ef1.ge.gam(44).or.(gam(44)/gam(1)).lt.1.01)go to 4000
      do 10 ie=2,44
      if(gam(ie).gt.ef1)go to 11
   10 continue
      ie=44
   11 ie1=ie-1
      g1=ef1
      go to 12
    1 g1=gam(1)
   12 if((gam(ie1+1)/gam(ie1)).gt.1.0002)go to 13
      ie1=ie1+1
      if(ie1.ge.43)go to 4000
      go to 12
   13 vala=s0*edist(ie1)/g1
c     Loop to integrate over electron Lorentz factors (gam)
      do 3000 ie=ie1,43
      g2=gam(ie+1)
      valb=s0*edist(ie+1)/g2
      val1=0.0
      val2=0.0
      gran1=0.0
      addit=0.0
c     Set up loop 1 to integrate over incident photon frequency anui
      anumax=amin1(anuf,dnu(nuhi))
      di1=di(1)
      id=2
      anumin=anuf/((2.0d0*g1)*(g1-ef)*
     ,  (1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1))))
      if(anumin.ge.anumax)go to 601
      if(anumin.gt.dnu(1))go to 2
      anumin=dnu(1)
      go to 5
    2 continue
      do 3 id=2,nuhi
      if(anumin.le.dnu(id))go to 4
    3 continue
      id=nuhi
    4 continue
      di1=0.0
      if(di(id).lt.1.0e-25.or.di(id-1).lt.1.0e-25)go to 5
      a=dlog10(di(id)/di(id-1))/alog10(dnu(id)/dnu(id-1))
      di1=di(id)*(anumin/dnu(id))**a
    5 ide=nuhi
      if(anumax.ge.dnu(nuhi))go to 8
      do 6 idd=id,nuhi
      if(anumax.le.dnu(idd))go to 7
    6 continue
      idd=nuhi
    7 a=dlog10(di(idd-1)/di(idd))/alog10(dnu(idd-1)/dnu(idd))
      die=di(idd-1)*(anumax/dnu(idd-1))**a
      ide=idd
      go to 9
    8 die=di(nuhi)
    9 continue
      anui1=anumin
      xx=1.0-ef/g1
      rat=homc2*anui1*g1*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1)))
      call xseckn(rat,ef,g1,xx,xsecc)
      xsecc=xsecc*(g1*homc2*anui1)
      val1=xsecc*(1.0e20/anui1)*(anuf/anui1)*di1*vala
      if(val1.lt.1.0d-40)val1=0.0d0
   25 continue
c     Loop 1 to integrate over incoming photon frequency anui for lower gam value
      do 600 nu=id,ide
      anui2=dnu(nu)
      di2=di(nu)
      if(nu.lt.ide)go to 30
      anui2=anumax
   30 rat=homc2*anui2*g1*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g1*g1)))
      call xseckn(rat,ef,g1,xx,xsecc)
      xsecc=xsecc*(g1*homc2*anui2)
      val2=xsecc*(1.0e20/anui2)*(anuf/anui2)*di2*vala
      if(val2.lt.1.0d-40)val2=0.0d0
  525 if(val1.eq.0.0d0.or.val2.eq.0.0d0)go to 845
      test=abs((anui1-anui2)/anui1)
      if(test.lt.0.001)go to 847
      ratnu=anui2/anui1
      a=1.0+dlog10(val2/val1)/alog10(ratnu)
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 845
      addit=val1*(ratnu**a-1.0)*anui1/a
      go to 846
  845 addit=0.5*(val1+val2)*(anui2-anui1)
  846 gran1=gran1+addit
  847 anui1=anui2
      val1=val2
      di1=di2
  600 continue
c     End anui loop 1
  601 continue
      ratg=g2/g1
      ratgl=alog10(ratg)
      valb=s0*edist(ie+1)/g2
      val1=0.0
      val2=0.0
      gran2=0.0
c     Set up loop 2 to integrate over incident photon frequency anui
      anumax=amin1(anuf,dnu(nuhi))
      di1=di(1)
      id=2
      anumin=anuf/((2.0d0*g2)*(g2-ef)*
     ,  (1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2))))
      if(anumin.ge.anumax)go to 1601
      if(anumin.gt.dnu(1))go to 1002
      anumin=dnu(1)
      go to 1005
 1002 continue
      do 1003 id=2,nuhi
      if(anumin.le.dnu(id))go to 1004
 1003 continue
      id=nuhi
 1004 continue
      di1=0.0
      addit=0.0
      if(di(id).lt.1.0e-25.or.di(id-1).lt.1.0e-25)go to 1005
      a=dlog10(di(id)/di(id-1))/alog10(dnu(id)/dnu(id-1))
      di1=di(id)*(anumin/dnu(id))**a
 1005 ide=nuhi
      if(anumax.ge.dnu(nuhi))go to 1008
      do 1006 idd=id,nuhi
      if(anumax.le.dnu(idd))go to 1007
 1006 continue
      idd=nuhi
 1007 a=dlog10(di(idd-1)/di(idd))/alog10(dnu(idd-1)/dnu(idd))
      die=di(idd-1)*(anumax/dnu(idd-1))**a
      ide=idd
      go to 1009
 1008 die=di(nuhi)
 1009 continue
      anui1=anumin
      xx=1.0-ef/g2
      rat=homc2*anui1*g2*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2)))
      call xseckn(rat,ef,g2,xx,xsecc)
      xsecc=xsecc*(g2*homc2*anui1)
      val1=xsecc*(1.0e20/anui1)*(anuf/anui1)*di1*valb
      if(val1.lt.1.0d-40)val1=0.0d0
 1025 gran2=0.0
c     Loop 2 to integrate over incoming photon frequency anui for upper gam value
      if(id.ge.ide)go to 1601
      do 1600 nu=id,ide
      anui2=dnu(nu)
      di2=di(nu)
      if(nu.lt.ide)go to 1030
      anui2=anumax
      di2=die
 1030 rat=homc2*anui2*g2*(1.0d0-cscat*dsqrt(1.0d0-1.0d0/(g2*g2)))
      call xseckn(rat,ef,g2,xx,xsecc)
      xsecc=xsecc*(g2*homc2*anui2)
      val2=xsecc*(1.0e20/anui2)*(anuf/anui2)*di2*valb
      if(val2.lt.1.0d-40)val2=0.0d0
      if(val1.eq.0.0.or.val2.eq.0.0)go to 1845
      test=abs((anui1-anui2)/anui1)
      if(test.lt.0.001)go to 1847
      ratnu=anui2/anui1
      a=1.0+dlog10(val2/val1)/alog10(ratnu)
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 1845
      addit=val1*(ratnu**a-1.0)*anui1/a
      go to 1846
 1845 addit=0.5*(val1+val2)*(anui2-anui1)
 1846 gran2=gran2+addit
 1847 continue
      anui1=anui2
      val1=val2
      di1=di2
 1600 continue
c     End anui loop 2
 1601 continue
      if(gran1.lt.1.0e-28.or.gran2.lt.1.0e-28)go to 2845
      if(ratgl.lt.0.0001)go to 2845
      a=1.0+alog10(gran2/gran1)/ratgl
      if(abs(a).lt.0.01.or.abs(a).gt.5.0)go to 2845
      addit=gran1*(ratg**a-1.0)*g1/a
      go to 2846
 2845 addit=0.5*(gran1+gran2)*(g2-g1)
 2846 gran=gran+addit
      g1=g2
      vala=valb
 3000 continue
c     End gam loop
      if(gran.gt.1.0e-15)ssca=gran*1.0d-20
      ssc=ssca
 6666 format('In ssc',i5,2x,1p16e12.3)
 4000 return
      end
