c     Turbulent Extreme Multi-Zone (TEMZ) model for blazar multi-waveband
c       light curves, SEDs, and (synchrotron) polarization 
c
c     Author: Alan Marscher, Dept. of Astronomy & Institute for Astrophysical Reearch, Boston University
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 allows an ordered helical magnetic field in addition to a turbulent field
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 the synchrotron
c       radiation plus a dust torus
c
c     To compile:
c     gfortran -o temzssc.out temzssc.f
c
c     To run: need input file temzsscinp.txt in working directory
c     If maps are recorded (nomap=0 in temzsscinp.txt file), need to create empty sub-directory 'mapsp'
c       where map files will be written
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),
     , bfij(125,112,169),bxfij(125,112,169),byfij(125,112,169),
     , bzfij(125,112,169),bpp(125,112,169),ebk(125,112,169),
     , ecofij(125,112,169,44),acofij(125,112,169,16),
     , alphsy(125,112,169,44),
     , ecij(125,112,169,48),alphec(125,112,169,48),
     , sscij(125,112,169,48),alphsc(125,112,169,48),
     , pq(112,169,64),pu(112,169,64),fpol(112,169,64),
     , flux(112,169,64),flsync(112,169,64),flcomp(112,169,64),
     , flec(112,169,64),flssc(112,169,64),prevsd(112,169,64),
     , gmax0(125,112,169),
     , egam(112,169,44),enofe(112,169,44),delta(125,112,169),
     , betadz(125,112,169),betadx(125,112,169),betady(125,112,169),
     , betad(125,112,169),gammad(125,112,169),ussold(125,112,169),
     , n0(125,112,169),betupx(169),betupy(169),betupz(169),
     , zcell(112,169),layer(169),xmaxx(112,169),angc(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),
     , 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),omptst(125,64),
     , useed(113),seedpk(113),tdelr(113),dustii(113,22),
     , dustf(113,22),dustnu(22),dusti(22),kd(112,169),
     , nu(64),snu(64),syseed(64),scseed(64),ssseed(64),
     , fsync2(64),phots(64),fmdall(64),
     , ssabs(64),sstau(64),srcfn(64),rnu(64),
     , ecflux(64),emisco(64),tauexp(64),taupp(64),
     , ajofnu(64),fsnoab(64),specin(64),sscflx(64),
     , 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,112),buy1(8,8,112),buz1(8,8,112),bux2(4,4,112),
     , buy2(4,4,112),buz2(4,4,112),
     , bux3(2,2,112),buy3(2,2,112),buz3(2,2,112),
     , 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 pq,pu,gammmx,gammmn,
     ,  betaup,gamup,sinz,cosz,phcell,
     ,  cosph,sinph,thlos,betau,gammau,opang,tanop,
     ,  zeta,tanz,betad,betadx,betady,betadz,slos,clos,
     ,  eta,betacs,gammad,bdx,bdy,bdz,
     ,  dcsth1,dcsth2,dth1,dth2,n0ave,
     ,  betamx,betamy,betamz,betamr,gamamr,
     ,  cosmd,tanmd,gmax0,ustob,tlfact,delt,
     ,  gamb,glow,gmrat,gmratl,gmratm,ggam,tloss,t2,
     ,  tlavg,tlavgo,
     ,  eterm1,eterm2,gammax,phots,betaux,betauy,
     ,  betauz,psi,cospsi,sinpsi,betat,
     ,  tanpsi,betd,gamd,angm,bmparx,bmpary,bmparz,
     ,  bmprpx,bmprpy,bmprpz,dotprd,
     ,  betupx,betupy,betupz,betatx,betaty,betatz,thetat,
     ,  sintht,costht,phit,btprpx,btprpy,btprpz,
     ,  btparx,btpary,btparz,sx,sy,sz,slx,sly,slz,
     ,  anx,any,anz,betadd,gammdd,csth1,csth2,
     ,  emisco,fsync2,ssabs,sstau,tauexp,ecflux,sscflx,fsnoab,
     ,  srcfn,ssseed,syseed,scseed,fmdall,ajofnu,
     ,  denm,xcl,ycl,zcl,xcls,ycls,zcls,xdf,ydf,zdf,ddf,xydf
      data ididb1/ 64*0/, ididb2/ 16*0/, ididb3/ 4*0/,
     ,  bfrac1/ 64*0.0/, bfrac2/ 16*0.0/, bfrac3/ 4*0.0/,
     ,  angrt1/ 64*0.0/, angrt2/ 16*0.0/, angrt3/ 4*0.0/,
     ,  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/ 7168*0.0/,buy1/ 7168*0.0/,buz1/ 7168*0.0/,
     ,  bux2/ 1792*0.0/,buy2/ 1792*0.0/,buz2/ 1792*0.0/,
     ,  bux3/ 448*0.0/,buy3/ 448*0.0/,buz3/ 448*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/
      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
      iwrote=0
c
      open (2,iostat=ios, err=9000, file='temzsscinp.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 (1,iostat=ios, err=9000, file='temzcheck.txt',
     , status='new')
      open (9,iostat=ios, err=9000, file='ompcheck.txt',
     , status='new')
      it=0
c     Input file format: 1st line characters, then values of parameters,
c       one per line
      read(2,9111)dumdum,nrad,itlast,modul,nomap,zred,dgpc,alpha,pn,pb,
     , bave,psdslp,uratio,rsize,gmaxmn,gmrat,gmin,betaup,betat,zeta,
     , thlos,opang,tdust,ldust,dtdist,dtrad,zdist0,bratio,bpitch
      close(2)
      if(nomap.eq.1)go to 1111
      write(filnam,"(a8,i4.4,a4)") "temzmap",it,".txt"
      open (7,iostat=ios, err=9000, file=dpath//filnam,
     ,   status='new')
      write(filnm2,"(a8,i4.4,a4)") "temzmpo",it,".txt"
      open (8,iostat=ios, err=9000, file=dpath//filnm2,
     ,   status='new')
      close(7)
      close(8)
 1111 continue
      iccnt=0
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,idth )
!$omp& schedule( dynamic, 1 )
      do 7701 inu=1,40
      idth = omp_get_thread_num()
      omptst(1,inu)=0.0
      write(9,9994)idth,inu,omptst(1,inu)
 7701 continue
!$omp end parallel do
      call itime(itarra)
      ita=itarra(1)+itarra(2)+itarra(3)
      iseed = rand(ita)
      randstart=rand(iseed)
      istop=0
c     Turn pair-production absorption on (ipprod=1) or off (ipprod=0)
      ipprod=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**2, density, and magnetic field with scale size
      tnexp=2.0/3.0
      tbexp=1.0/3.0
      tvexp=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; from Kolmogorov spectrum
c     Note: this assumes that the turbulent velocity does not exceec c/sqrt(3), otherwise
c       need to use 4-velocity 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(1,9433)nrad,ancol,bnorm0,bnorm1,bnorm2,bnorm3,anorm0,
     ,  anorm1,anorm2,anorm3,vnorm0,vnorm1,vnorm2,vnorm3
 9433 format(i4,f5.0,1p12e10.3)
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
c     An SED will be printed out every ispecs time steps
      ispecs=1
      ispec=1
c     Set up frequencies
      do 1 inu=1,64
      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
      if(inu.eq.40)nu(inu)=1.0e20
      if(inu.eq.41)nu(inu)=2.3713e20
      if(inu.eq.42)nu(inu)=5.6234e20
      if(inu.eq.43)nu(inu)=1.3335e21
      if(inu.eq.44)nu(inu)=3.1623e21
      if(inu.eq.45)nu(inu)=1.0e22
      if(inu.eq.46)nu(inu)=2.42e22
      if(inu.gt.46)nu(inu)=2.42e22*10.0**(-12.25+0.25*(inu+3))
      snu(inu)=nu(inu)
    1 continue
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)
c     This is an option for determining probability distribution of gamma_max
      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 mean energy density of electrons 
c        to that of mag. field in upstream (unshocked) plasma
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
c     Bulk Lorentz factor of laminar component of flow upstream of shock
      gamup=1.0d0/dsqrt(1.0d0-betaup**2)
c     Speed downstream of conical shock if turbulent velocity is ignored;
c        for setting cell length zsize and for first estimate of time delay
      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(1,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 (8 is arbitrary, but reasonable)
      zsize=zrat*rsize/tanz
c     seff is the integral over the volume of a cell divided by the square of the
c       distance from its center; used to calculate the seed photon field from within a cell
      seff=0.0
      rzrat=rsize/zsize
      ahi=0.5*pi
      do 7 iv=1,100
      rz1=(iv-1)*0.01*rzrat
      rz2=iv*0.01*rzrat
      alo=ahi
      ahi=atan2((0.5/rz2),1.0)
      seff=seff+0.5*(rz2-rz1)*(alo+ahi)
    7 continue
      write(1,9521)rzrat,seff
 9521 format('rzrat, seff ',2f10.4)
c     anglim is used later when path length through a cell is needed
      anglim=atan2(2.0*rsize,zsize)
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
      delobs=1.0d0/(gamd*(1.0d0-betd*clos))
c     Time step in observer's frame in days: crossing time of a single cell at
c       laminar velocity
      dtfact=(1.0d0-betd*clos)/betd
      dtime=1190.0*zsize*dtfact*zred1
c     Set arbitrary reference point in array of master time variations from random
c       fluctuations of energy density inject into jet upstream of the shock
      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 shock convergence point from z value where conical shock 
c        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=ip0+icells/betd
      ndmax=sqrt(4.0*(rshock/rsize)**2+icells*icells)
      write(5,9494)kmd,ndmax
 9494 format('# kmd: ',i5,'    ndmax: ',i5)
c
c     Time delay is measured relative to time when plasma passes
c       upstream end of standing shock, located distance zrf from
c       upstream start of shock
      zrf=0.0
      xrf=0.0
c     Distance of downstream end of standing shock from vertex of jet cone
c     Need for calculation of seed photon field from dusty torus
      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)
      dustii(id,inu)=dusti(inu)
    3 continue
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=alog10(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,filld,zdist0,useed(5),psdslp
      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,filld,zdist0,useed(5),psdslp
      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,filld,zdist0,useed(5),psdslp
 6665 format('#Radius in no. of cells: ',i3,4x,'No. of time steps:',i6,
     ,   ' time step: ',f8.5/
     ,  '#redshift: ',f5.3,2x,'Distance in Gpc: ',f5.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): ',f7.5/
     ,  '#Min. value of gamma_max: ',f8.1,2x,'ratio of max. to min. ',
     ,  'values of gamma_max: ',f5.1/'#gamma_min: ',f8.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   ','# Area filling factor of dust emission: ',e9.2/
     ,  '#Distance of shock from central engine: ',f5.2,' pc'/
     ,  '#Energy density of seed photons in plasma frame: ',e9.2/
     ,  '# -Slope of PSD: ',f5.1,'    SSC, no Mach disk')
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
      spsig=0.0
      spxsig=0.0
c     If a pulse of extra energy input is desired - usually for testing purposes -
c       ipulse sets a starting point in the variation array; if desired, code
c       needs to be added to insert the pulse over a range of times
c      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)
 4998 continue
      spsig=sqrt(spsig)
      write(1,9992)nrad,ip0,dstart,zsize,zshock,clos,slos,
     ,   gammdd,delobs,dtfact,spxsig,spxmn,spsig,psdmn,spmin
c      normalize spsd by its minimum value
      spnorm=spmin
      do 4999 ip=1,ndim
      spsd(ip)=spsd(ip)/spnorm
c     For testing purposes or if no fluctuations at input are desired
c      spsd(ip)=1.0
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
 9345 format(i10,1p2e14.5)
      close(10)
c     Set parameters of each cell at initial time
c     There are jcells columns of cells, with imax (<= icells) cells per column
c     The central cell (j=1) is ignored in this version of TEMZ
      write(4,6667)
      write(5,6668)
      i=1
      j=1
      ncells=0
      jzero=2
      do 999 nnn=1,(nrad-1)
      do 998 jcnt=1,(6*nnn)
      j=j+1
      annn=nnn
      layer(j)=nrad-nnn
      angc(j)=60.0/(rad*annn)
      rcell(j)=2.0*rsize*annn
      xcell(j)=rcell(j)*cos(angc(j)*(j-jzero))
      ycell(j)=rcell(j)*sin(angc(j)*(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
      nouter(j)=imax(j)
      cosph(j)=xcell(j)/rcell(j)
      sinph(j)=ycell(j)/rcell(j)
      phcell(j)=datan2(sinph(j),cosph(j))
      tanpsi(j)=rcell(j)/(zsvtex-zcol)
      psi(j)=datan(tanpsi(j))
      cospsi(j)=dcos(psi(j))
      sinpsi(j)=tanpsi(j)*cospsi(j)
c     Velocity vector of laminar component of pre-shock flow
      betupx(j)=betaup*cosph(j)*sinpsi(j)
      betupy(j)=betaup*sinph(j)*sinpsi(j)
      betupz(j)=betaup*cospsi(j)
c     Set up magnetic field cell 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
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,2,-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
      write(1,9192)j,imax(j),
     ,ididd1(ng1x(j),ng1y(j),layer(j)),ididd2(ng2x(j),ng2y(j),layer(j)),
     , ididd3(ng3x(j),ng3y(j),layer(j)),ididc1(ng1x(j),ng1y(j)),
     , ididc2(ng2x(j),ng2y(j)),ididc3(ng3x(j),ng3y(j)),ididj1(j),
     , ididj2(j),ididj3(j),ididk1(j,layer(j)),ididk2(j,layer(j)),
     , ididk3(j,layer(j))
      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)
      ncells=ncells+1
      zcell(i,j)=zshock+(i-1)*zsize-rcell(j)/tanz
      xmaxx(i,j)=rbound*abs(cos(angc(j)*(j-jzero)))-xcell(j)
      rhole=(zshock-zcell(i,j)+(i-1)*zsize)*tanz
      xhole=rhole*abs(cos(angc(j)*(j-jzero)))-xcell(j)
      xmaxx(i,j)=xmaxx(i,j)-xhole
      write(1,9191)j,i,imax(j),ncells,xcell(j),ycell(j),zcell(i,j),
     ,ididd1(ng1x(j),ng1y(j),layer(j)),ididd2(ng2x(j),ng2y(j),layer(j)),
     , ididd3(ng3x(j),ng3y(j),layer(j)),ididc1(ng1x(j),ng1y(j)),
     , ididc2(ng2x(j),ng2y(j)),ididc3(ng3x(j),ng3y(j)),ididj1(j),
     , ididj2(j),ididj3(j),ididk1(j,layer(j)),ididk2(j,layer(j)),
     , ididk3(j,layer(j)),xhole,rhole,xmaxx(i,j)
  887 continue
  888 continue
 9191 format(4i5,3f11.5,12i3,3f11.5)
 9192 format(2i5,12i3)
      write(*,*)'no. of cells = ',ncells
      imax(1)=1
      zcell(1,1)=zshock
c     kd(i,j) is the z distance from the most upstream cell in units of zsize
      do 33 j=1,jcells
      do 3333 i=1,imax(j)
      kd(i,j)=(zcell(i,j)-zrf)*betadd/zsize
      if(kd(i,j).gt.icells)kd(i,j)=icells
      if(kd(i,j).lt.1)kd(i,j)=1
 3333 continue
   33 continue
c
c     After initial set-up, time loop starts here
    9 i=1
c     it is the time step; itt is the time step after the initial icells
c       time steps; results of code are recorded starting when itt >= 1
c       so that transients are no longer present
      it=it+1
      itt=it-icells-2
      ididj3(1)=0
c
c     Start loop over all cells in first layer to set up physical parameters
c
c      write(*,*)'start loop over first layer, it =',it
      do 80 j=jcells,2,-1
c
c     Pre-shock energy density of cell is obtained from power spectrum array
c       point (idelay+ip0); starting point is somewhat arbitrary
c     Plasma farther downstream was injected at an earlier time
c     nd is the time delay of a cell at (x,y,z) relative to position (0,0,zrf)
c     Lower nd corresponds to earlier time, farther from reference point
      nd=ndmax
      k=kd(i,j)
      idelay=dstart+it-0.5-((zcell(i,j)-zrf)*clos+
     ,   (xcell(j)-xrf)*slos)*betadd/zsize
      if(idelay.le.0.or.idelay.gt.(ndim-ip0))write(1,9222)idelay,
     ,  dstart,j,i,it,imax(j),ndim,zshock,zrf,zcell(i,j),rcell(j),
     ,  xrf,xcell(j),betadd,zsize,clos
      if(idelay.le.0.or.idelay.gt.(ndim-ip0))go to 9000
c     Determine B vector of cell assuming random magnetic field orientation
c        in each of 4 sections of different sizes
c     Interpolate inside loop to get direction for cells not on
c       a section boundary
      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
      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))=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))
      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 transverse to the shock normal
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(j)=betaup*cosph(j)*sinpsi(j)
      betupy(j)=betaup*sinph(j)*sinpsi(j)
      betupz(j)=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     Velocity vector of the pre-shock plasma including turbulent component
      dotprd=betupx(j)*betatx+betupy(j)*betaty+betupz(j)*betatz
      btparx=dotprd*betupx(j)/betaup**2
      btpary=dotprd*betupy(j)/betaup**2
      btparz=dotprd*betupz(j)/betaup**2
      btprpx=betatx-btparx
      btprpy=betaty-btpary
      btprpz=betatz-btparz
      betaux=(betupx(j)+btparx+btprpx/gamup)/(1.0d0+dotprd)
      betauy=(betupy(j)+btpary+btprpy/gamup)/(1.0d0+dotprd)
      betauz=(betupz(j)+btparz+btprpz/gamup)/(1.0d0+dotprd)
      betau=dsqrt(betaux**2+betauy**2+betauz**2)
      gammau=1.0d0/dsqrt(1.0d0-betau**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,betauy,betauz,sx,sy,sz,
     ,  betadx(nd,i,j),betady(nd,i,j),betadz(nd,i,j),betad(nd,i,j),
     ,  gammad(nd,i,j),eta)
c     Calculate the Doppler factor
      bdx=betadx(nd,i,j)
      bdy=betady(nd,i,j)
      bdz=betadz(nd,i,j)
      betacs=bdx*slos+bdz*clos
      delta(nd,i,j)=1.0d0/(gammad(nd,i,j)*(1.0d0-betacs))
c     Upstream velocity component parallel to shock normal
      velpar=betaux*anx+betauy*any+betauz*anz
c     Calculate upstream B field components parallel + perpendicular to shock front
      call bcalc(betaux,betauy,betauz,betau,gammau,
     ,  anx,any,anz,bux,buy,buz,bparx,bpary,bparz,bprpx,bprpy,bprpz,
     ,  bpar,bprp)
      bx=bparx+eta*bprpx
      by=bpary+eta*bprpy
      bz=bparz+eta*bprpz
c     Calculate downstream magnetic field, including component bperp that  
c       is perpendicular to 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(nd,i,j),betady(nd,i,j),betadz(nd,i,j),
     , betau,gammau,slx,sly,slz,bx,by,bz,dum1,dum2,dum3,dum4,
     ,  dum5,dum6,dum7,bperp)
      bfij(nd,i,j)=sqrt(bx**2+by**2+bz**2)
      bxfij(nd,i,j)=bx
      byfij(nd,i,j)=by
      bzfij(nd,i,j)=bz
      bpp(nd,i,j)=bperp
c     Set reference B strength arbitrarily as that of cell (1,jcells) at it=1
      if(it.eq.1.and.j.eq.jcells)then
        bpref=bperp
        blref=bfij(nd,i,j)
      else
        continue
      end if
      n0mean=spsd(idelay+ip0)*(anorm0*an00+anorm1*an01+anorm2*an02+
     ,  anorm3*an03)
      n0(nd,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
      if(it.eq.1)then
        ustob=1.0
      else
        ustob=8.0*pi*ussold(nd,i,j)/(bfij(nd,i,j))**2
      end if
      htlim=1.0/(gammau*velpar)
      bfang=bprp/bpar
c     For HT-frame 'superluminal' shock, no particle acceleration
      if(bfang.gt.htlim)then
        gmax0(nd,i,j)=gmaxmn
      else
c     For HT-frame subluminal shock, gamma_max,0 is from Baring et al. (2017),
c       equating cooling time with acceleration time with mean-free path
c       of pitch-angle scattering proportional to gamma
        gmax0(nd,i,j)=gmaxmx*((bfij(nd,i,j)/blref)/
     ,   (bfij(nd,i,j)*bfij(nd,i,j)*(1.0+ustob)/(blref*blref)))**0.3333
         if(gmax0(nd,i,j).lt.gmaxmn)gmax0(nd,i,j)=gmaxmn
      end if
c     Other schemes tried in the past:
c      gmax0(nd,i,j)=gmaxmx-(gmaxmx-gmaxmn)*(bfang/htlim)
c      gmax0(nd,i,j)=gmaxmx*(bprp*bprp/(bpar*bpar+bprp*bprp))
c      if(gmax0(nd,i,j).lt.gmaxmn)gmax0(nd,i,j)=gmaxmn
c      if(itt.gt.0)write(1,9059)i,j,betatx,betaty,
c     ,  betatz,bdx,bdy,bdz,delta(nd,i,j),bx,by,bz,htlim,bfang,
c     ,  gmax0(nd,i,j)
      if(itt.gt.0)write(1,9059)i,j,nd,bfij(nd,i,j),
     ,  bprp,bpar,bdx,bdy,ussold(nd,i,j),ussold(nd-1,i,j),htlim,bfang,
     ,  ustob,blref,gmaxmx,gmax0(nd,i,j)
 9059 format('+++ ',3i5,1p14e11.4)
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(nd,i,j)=gmaxmx
c      if(pexp.lt.0.0)gmax0(nd,i,j)=((gmaxmx**pexp-gmaxmn**pexp)*xrand+
c     ,  gmaxmn**pexp)**(1.0/pexp)
      do 15 ig=1,43
      if(gmax0(nd,i,j).le.gcnt(ig+1).and.gmax0(nd,i,j).ge.gcnt(ig))
     , igcnt(ig)=igcnt(ig)+1
   15 continue
c
c     Initialize radiation field via approximations, since calculations have not
c       yet been carried out for earlier times
c
c     Cosine of angle of SSC scattering; use =0.7 for initialization
      cscat=0.7
c     Compute seed photon field from hot dust in torus
      id=zcell(i,j)/zsize+1
      if(id.lt.1)id=1
      if(id.gt.(icells+1))write(1,6992)it,i,j,id,md,k,zcell(i,j),
     ,  zshock,zsize
      if(id.gt.icells)id=icells
      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(nd,i,j))/(1.0d0-csth1*betad(nd,i,j))
      dcsth2=-(csth2-betad(nd,i,j))/(1.0d0-csth2*betad(nd,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(nd,i,j)*(1.0d0-betad(nd,i,j)*0.5d0*(csth1+csth2))
      do 73 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 7073
      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 7074
 7073 dusti(inu)=0.0
 7074 continue
   73 continue
c     Use radiation density from previous time step to calculate energy loss
      if(it.eq.1)then
        ustob=0.0
      else
        if(nd.le.1)ustob=8.0*pi*ussold(nd,i,j)/(bfij(nd,i,j))**2
        if(nd.gt.1)ustob=8.0*pi*ussold(nd-1,i,j)/(bfij(nd,i,j))**2
      end if
      tlfact=cc2*bfij(nd,i,j)**2*(1.0+ustob)
      delt=zsize*yr*3.26/(gammad(nd,i,j)*betad(nd,i,j))
      gamb=gmax0(nd,i,j)/(1.0+tlfact*gmax0(nd,i,j)*delt)
      glow=gmin/(1.0+tlfact*gmin*delt)
      gmrat=0.99*gmax0(nd,i,j)/gmin
      gmratl=dlog10(gmrat)/38.0
      gmratm=dlog10(gmin/glow)/5.0
      ibreak=0
      do 74 ie=1,44
      egam(i,j,ie)=gmin*10.0**(gmratl*(ie-6))
      if(ie.lt.6)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)
      tloss=(gmax0(nd,i,j)-ggam(ie))/(tlfact*ggam(ie)*gmax0(nd,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 1089
      if(ie.ge.6)enofe(i,j,ie)=n0(nd,i,j)/((sen-1.0)*tlfact*
     ,  ggam(ie)**(sen+1.0))*(1.0d0-eterm2**(sen-1.0))
      if(ie.lt.6)enofe(i,j,ie)=n0(nd,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
 1089 edist(ie)=enofe(i,j,ie)
   74 continue
      ebk(nd,i,j)=egam(i,j,6)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 773 inu=1,48
      ecij(nd,i,j,inu)=ecdust(nu(inu+16))
  773 continue
!$omp end parallel do
      bperpp=bpref
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 774 inu=1,64
      ssseed(inu)=0.0
      if(inu.gt.44)go to 7741
      ecofij(nd,i,j,inu)=ajnu(nu(inu))*bperpp
      if(inu.lt.17)acofij(nd,i,j,inu)=akapnu(nu(inu))
      if(ecofij(nd,i,j,inu).ne.ecofij(nd,i,j,inu))then
      write(5,6993)it,nd,i,j,inu,inu,inu,
     ,  ecofij(nd,i,j,inu),bpref,
     ,  ggam(1),edist(1),ggam(44),edist(44),n0mean
      stop
      else
      continue
      endif
c     Very rough approximation of synchrotron seed photon field
 7741 continue
      if(it.eq.1) then
      if(inu.lt.45)syseed(inu)=ecofij(nd,i,j,inu)*seff*fgeom*
     ,   (emfact*parsec)*nrad**2
      if(inu.lt.45)ssseed(inu)=syseed(inu)
      else
      ssseed(inu)=prevsd(i,j,inu)
      endif
  774 continue
!$omp end parallel do
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 775 inu=1,48
      sscij(nd,i,j,inu)=ssc(nu(inu+16))
c      if(i.eq.1.and.j.eq.169)write(3,9996)snu(inu),ssseed(inu),
c     ,  ecij(nd,i,j,inu),sscij(nd,i,j,inu),prevsd(i,j,inu)
  775 continue
!$omp end parallel do
      do 776 inu=2,44
      alphsy(nd,i,j,inu)=0.0
      if(ecofij(nd,i,j,(inu-1)).gt.1.0e-30.and.
     ,   ecofij(nd,i,j,inu).gt.1.0e-30)alphsy(nd,i,j,inu)=
     ,   alog10(ecofij(nd,i,j,(inu-1))/ecofij(nd,i,j,inu))/
     ,   alog10(nu(inu)/nu(inu-1))
  776 continue
      do 777 inu=2,48
      alphec(nd,i,j,inu)=0.0
      if(ecij(nd,i,j,(inu-1)).gt.1.0e-30.and.
     ,   ecij(nd,i,j,inu).gt.1.0e-30)alphec(nd,i,j,inu)=
     ,   alog10(ecij(nd,i,j,(inu-1))/ecij(nd,i,j,inu))/
     ,   alog10(nu(inu+16)/nu(inu+15))
      alphsc(nd,i,j,inu)=0.0
      if(sscij(nd,i,j,(inu-1)).gt.1.0e-30.and.
     ,   sscij(nd,i,j,inu).gt.1.0e-30)alphsc(nd,i,j,inu)=
     ,   alog10(sscij(nd,i,j,(inu-1))/sscij(nd,i,j,inu))/
     ,   alog10(nu(inu+16)/nu(inu+15))
  777 continue
c     Approximate spectral index of lowest frequency to be same as next-lowest
      alphsy(nd,i,j,1)=alphsy(nd,i,j,2)
      alphec(nd,i,j,1)=alphec(nd,i,j,2)
      alphsc(nd,i,j,1)=alphsc(nd,i,j,2)
      if(it.ne.1)go to 780
c
c     Copy initial values of physical parameters to all cells for retarded
c       times earlier than the start of the computation
c     Code will not record emission results until after these initial values
c       are replaced by computed values (i.e., after icells time steps)
      do 79 ii=1,imax(j)
      do 78 ndn=1,ndmax
      bxfij(ndn,ii,j)=bxfij(ndmax,i,j)
      byfij(ndn,ii,j)=byfij(ndmax,i,j)
      bzfij(ndn,ii,j)=bzfij(ndmax,i,j)
      bfij(ndn,ii,j)=bfij(ndmax,i,j)
      bpp(ndn,ii,j)=bpp(ndmax,i,j)
      betad(ndn,ii,j)=betad(ndmax,i,j)
      betadx(ndn,ii,j)=betadx(ndmax,i,j)
      betady(ndn,ii,j)=betady(ndmax,i,j)
      betadz(ndn,ii,j)=betadz(ndmax,i,j)
      gammad(ndn,ii,j)=gammad(ndmax,i,j)
      delta(ndn,ii,j)=delta(ndmax,i,j)
      ebk(ndn,ii,j)=ebk(ndmax,i,j)
      gmax0(ndn,ii,j)=gmax0(ndmax,i,j)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 76 inu=1,48
      if(inu.le.44)then
       ecofij(ndn,ii,j,inu)=ecofij(ndmax,i,j,inu)
       alphsy(ndn,ii,j,inu)=alphsy(ndmax,i,j,inu)
       if(inu.le.16)acofij(ndn,ii,j,inu)=acofij(ndmax,i,j,inu)
      endif
c     Seed photon field from dusty torus emission + 1st-order SSC
      ecij(ndn,ii,j,inu)=ecij(ndmax,i,j,inu)
      sscij(ndn,ii,j,inu)=sscij(ndmax,i,j,inu)
      alphec(ndn,ii,j,inu)=alphec(ndmax,i,j,inu)
      alphsc(ndn,ii,j,inu)=alphsc(ndmax,i,j,inu)
   76 continue
!$omp end parallel do
   78 continue
   79 continue
  780 continue
c      if(i.eq.1.or.i.eq.2)write(1,9992)nd,i,j,ecofij(2,i,j,14),
c     , alphsy(2,i,j,14),acofij(2,i,j,4),ecij(2,i,j,30),alphec(2,i,j,30),
c     , sscij(2,i,j,30),alphsc(2,i,j,30),ssseed(14),ssseed(30),
c     , bperpp,egam(i,j,6)
   80 continue
c      call cpu_time(cput)
c      write(*,*)'end loop over first layer, it =',it
c
c     End initial loop over all cells in first layer
c
c     Start another loop over all cells to calculate the
c       electron energy distribution and the emission
c
c      write(3,9666)ecofij(1,2,169,6),acofij(1,2,169,6)
c 9666 format('after first loop  ',1p2e12.3)
      do 100 j=jcells,2,-1
      tlavg=0.0d0
      do 1001 i=1,imax(j)
c     nd is the time index, = ndmax for current time, lower in past
      nd=ndmax
      do 101 inu=1,64
      emisco(inu)=0.0
      ecflux(inu)=0.0
      fmdall(inu)=0.0
      ssseed(inu)=0.0
  101 continue
      bdx=betadx(nd,i,j)
      bdy=betady(nd,i,j)
      bdz=betadz(nd,i,j)
c      betacs=bdx*slos+bdz*clos
c      delta(nd,i,j)=1.0d0/(gammad(nd,i,j)*(1.0d0-betacs))
      bperp=bpp(nd,i,j)
      bfld=bfij(nd,i,j)
c     Calculate seed photon (SP) field from all other cells in
c     cell (nd,i,j) frame; include time delays and relative Doppler factors
      do 183 jj=2,jcells
      do 83 ii=1,imax(jj)
c     Determine the position vector between the two cells
      xcl=xcell(j)
      ycl=ycell(j)
      zcl=zcell(i,j)
      xcls=xcell(jj)
      ycls=ycell(jj)
      zcls=zcell(ii,jj)
      xdf=xcl-xcls
      ydf=ycl-ycls
      zdf=zcl-zcls
      xydf=dsqrt(xdf*xdf+ydf*ydf)
      ddf=dsqrt(xdf*xdf+ydf*ydf+zdf*zdf)
c     Light-travel time delay between cells in units of time increments
      ndl=nd-ddf*betadd/zsize
      if(ndl.lt.1)ndl=1
c      if(ndl.gt.ndmax)write(1,9144)ndl,nd,ddf,betadd,zsize
c 9144 format(2i6,3f11.5)
      bdsx=betadx(ndl,ii,jj)
      bdsy=betady(ndl,ii,jj)
      bdsz=betadz(ndl,ii,jj)
      denm=1.0d0-(bdx*bdsx+bdy*bdsy+bdz*bdsz)
c     Determine velocity of the cell (nd,i,j) plasma relative to the
c       plasma of SP cell
      bmparx=(bdsx*bdx+bdsy*bdy+bdsz*bdz)*bdx/
     ,  (betad(nd,i,j)*betad(nd,i,j))
      bmpary=bmparx*bdy/bdx
      bmparz=bmparx*bdz/bdx
      bmprpx=bdx-bmparx
      bmprpy=bdy-bmpary
      bmprpz=bdz-bmparz
      betamx=(bdx-bmparx-bmprpx/gammad(nd,i,j))/denm
      betamy=(bdy-bmpary-bmprpy/gammad(nd,i,j))/denm
      betamz=(bdz-bmparz-bmprpz/gammad(nd,i,j))/denm
      betamr=dsqrt(betamx*betamx+betamy*betamy+betamz*betamz)
      gamamr=1.0d0/dsqrt(1.0d0-betamr*betamr)
c     Calculate angle between SP and scattered l.o.s. photon; need single precision x,y,z
      xxdf=xdf
      yydf=ydf
      zzdf=zdf
      call scatcs(bdx,bdy,bdz,clos,0.0d0,slos,xxdf,yydf,
     ,   -zzdf,cscat)
c     Determine relative Doppler factor
      deltmd=1.0
      bpcorr=1.0
      if(betamr.lt.1.0d-3.or.ddf.lt.1.0e-6)go to 1144
      cosmd=(betamx*xdf+betamy*ydf+
     ,  betamz*zdf)/(betamr*ddf)
      tanmd=dsqrt(1.0d0/(cosmd*cosmd)-1.0d0)
      angm=datan(tanmd)
      tanmd=dtan(2.0*datan(dtan(0.5d0*angm)/
     ,  (gamamr*(1.0d0+betamr))))
      cosmd=1.0d0/dsqrt(1.0d0+tanmd*tanmd)
      deltmd=1.0d0/(gamamr*(1.0d0-betamr*cosmd))
c     SP cell's B field in frame of cell plasma transverse to cell's
c       line of sight to SP cell (for synchrotron calculation)
      call bcalc(betamx,betamy,betamz,betamr,
     ,  gamamr,xdf,ydf,zdf,bxfij(ndl,ii,jj),byfij(ndl,ii,jj),
     ,  bzfij(ndl,ii,jj),dum1,dum2,dum3,dum4,dum5,dum6,dum7,bmperp)
      bpcorr=bmperp/bpref
 1144 continue
c
c     Calculate emission from various processes
c
c     Calculate synchrotron flux of SP in frame of cell plasma
c
c     Effective path across SP cell in direction of emitting cell
      path=seff
      if(ii.eq.i.and.jj.eq.j)go to 144
      angl=acos(abs(zdf/ddf))
      if(angl.gt.anglim)path=2.0*rsize*ddf/abs(xydf)
      if(angl.le.anglim)path=zsize*ddf/abs(zdf)
  144 continue
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu,abreak,absexp,sexp )
!$omp& schedule( dynamic, 1 )
      do 146 inu=1,44
      syseed(inu)=0.0
      scseed(inu)=0.0
      if(inu.le.20)syseed(inu+44)=0.0
      if(inu.le.20)scseed(inu+44)=0.0
      sstau(inu)=0.0
      tauexp(inu)=0.0
      ssabs(inu)=0.0
      ajofnu(inu)=0.0
      if(ecofij(ndl,ii,jj,inu).lt.1.0e-30.or.alphsy(ndl,ii,jj,inu).
     ,   gt.9.0)go to 145
c     Critical frequency of electrons at low-energy break of energy 
c        distribution (energy number 6)
      abreak=2.8e6*(bpref*bpcorr)*ebk(ndl,ii,jj)*ebk(ndl,ii,jj)
c     Synchrotron intensity for inverse Compton calculation
      if(inu.gt.16)go to 142
c     Correct previously calculated absorption coefficient for actual value
c     of B_perp; use different exponent below abreak
      if ((nu(inu)/deltmd).lt.abreak) then
      absexp=1.6667
      sexp=0.66667
      else
      absexp=0.5*(sen+4.0)
      sexp=0.5*(sen+2.0)
      endif
      ssabs(inu)=acofij(ndl,ii,jj,inu)*1.02e4*(sen+2.0)*bpref*
     ,  (parsec/nu(inu))/nu(inu)
      ssabs(inu)=ssabs(inu)*bpcorr**sexp*deltmd**absexp
  142 sstau(inu)=ssabs(inu)*path
      ajofnu(inu)=ecofij(ndl,ii,jj,inu)*bpcorr**(1.0+
     ,  alphsy(ndl,ii,jj,inu))*deltmd**(2.0+alphsy(ndl,ii,jj,inu))*path*
     ,  (parsec*emfact)
      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 if(ii.ne.i.or.jj.ne.j)then
      syseed(inu)=syseed(inu)*0.25*rsize**2*fgeom/(ddf*ddf)
      else
      syseed(inu)=syseed(inu)*seff*fgeom
      endif
c     Now estimate exponential attenuation of seed photons by
c      synchrotron self-absorption in intervening cells
c     Approximate that the absorption coefficient is that of SP cell
      tauexp(inu)=sstau(inu)*ddf/path
      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
      fmdall(inu)=fmdall(inu)+syseed(inu)
  146 continue
!$omp end parallel do
c
c     Calculate SP inverse Compton flux in frame of cell plasma
c
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 147 inu=1,64
c     Inverse Compton mean SP intensity for 2nd-order
c       inverse Compton calculation
      if(inu.le.16)go to 1147
      scseed(inu)=sscij(ndl,ii,jj,inu-16)*deltmd**
     ,  (2.0+alphsc(ndl,ii,jj,inu-16))
      if(ii.ne.i.or.jj.ne.j)then
      scseed(inu)=scseed(inu)*zsize*0.25*rsize*rsize*
     ,  (amjy*parsec)/(ddf*ddf)
      else
      scseed(inu)=scseed(inu)*seff*fgeom*(parsec*amjy)
      endif
      fmdall(inu)=fmdall(inu)+scseed(inu)
 1147 continue
  147 continue
!$omp end parallel do
   83 continue
  183 continue
c      call cpu_time(cput)
c      write(*,*)'end loop over all other cells, i,j = ',i,j
c      if(i.gt.1)go to 2146
      do 2145 inu=1,64
      ssseed(inu)=fmdall(inu)
      prevsd(i,j,inu)=ssseed(inu)
 2145 continue
c      go to 2147
c 2146 continue
c 2147 continue
c
c     End calculation of seed photons from all other cells
c
c     Calculate energy distribution and emission from cell (i,j)
      nuhi=2
      do 2148 inu=17,64
      if(fmdall(inu).gt.0.0)nuhi=inu
 2148 continue
c     Estimate seed photon field from same cell & from other cells
c     Calculate seed photon energy density for energy loss calculation
      usp=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 2149
      usp=usp+0.5/c*snu(inu-1)/(1.0-aaa)*fmdall(inu-1)*
     ,  ((snu(inu)/snu(inu-1))**(1.0-aaa)-1.0)
      go to 2150
 2149 usp=usp+0.5*(fmdall(inu-1)*snu(inu-1)/c)*
     ,  alog(snu(inu)/snu(inu-1))
 2150 continue
      if(usp.ne.usp.or.usp.lt.1.0e-30)write(1,9992)it,i,j,aaa,
     ,  snu(inu-1),snu(inu),fmdall(inu-1),fmdall(inu),usp
  148 continue
  149 continue
      ussold(nd,i,j)=usp
  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+1))write(1,6992)it,i,j,id,md,k,zcell(i,j),
     ,  zshock,zsize
 6992 format('*x*x*x ',6i6,15f11.5)
      if(id.gt.icells)id=icells
      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(nd,i,j))/(1.0d0-csth1*betad(nd,i,j))
      dcsth2=-(csth2-betad(nd,i,j))/(1.0d0-csth2*betad(nd,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(nd,i,j)*(1.0d0-betad(nd,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
      useedd=useed(id)*(tdel/tdelr(id))**2/(1.0d0+betad(nd,i,j)*csang)
 6101 continue
      if(i.eq.1)ussold(nd,i,j)=useedd+usp
      ustob=8.0*pi*(useedd+usp)/(bfij(nd,i,j))**2
c     Calculate energy distribution for the cell
c     Derive the maximum electron energy from gmax0 and solution to equation
c       d gamma/dt = - cc2*(b**2+8*pi*useed)*gamma**2
      tlfact=cc2*bfij(nd,i,j)**2*(1.0+ustob)
      cellno=i
      tlavgo=tlavg
      if(i.eq.1)tlavgo=tlfact
      tlavg=(tlavg*(cellno-1.0)+tlfact)/cellno
      delt=zsize*yr*3.26/(gammad(nd,i,j)*betad(nd,i,j))
      t1=(i-1)*delt
      t2max=i*delt
      gammmn=gmin
      if(i.eq.1)then
        glow=gmin/(1.0+tlfact*gmin*delt)
        gammmx=gmax0(nd,i,j)
c     Break in energy distribution occurs at gamma_max evaluated at downstream end of cell
        gamb=0.9999*gammmx/(1.0+tlfact*gammmx*delt)
      else
        glow=gmin/(1.0+tlavg*gmin*t2max)
        gammmx=gmax0(nd,i,j)/(1.0+tlavgo*gmax0(nd,i,j)*t1)
        gamb=0.9999*gammmx/(1.0+tlavg*gammmx*t2max)
      endif
        gmratm=dlog10(gmin/glow)/5.0
c     Need 0.9999 factor since gmax is a singularity
      gmrat=0.9999*gammmx/gammmn
      gmratl=dlog10(gmrat)/38.0
      ibreak=0
      if(gamb.le.glow)ibreak=1
c     Set up energies so that breaks in E distribution occur at one of the array values
      do 84 ie=1,44
      egam(i,j,ie)=gammmn*10.0**(gmratl*(ie-6))
      if(ie.lt.6)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)
   84 continue
!$omp parallel do
!$omp& default( shared )
!$omp& private( ie,tloss,t2,eterm1,eterm2,efill )
!$omp& schedule( dynamic, 1 )
      do 90 ie=1,44
      enofe(i,j,ie)=0.0
      tloss=(gmax0(nd,i,j)-ggam(ie))/(tlavgo*ggam(ie)*gmax0(nd,i,j))
c      exxx=0.0
c      exx2=0.0
c      if(i.gt.1.and.tloss.le.t1)go to 188
      if(i.gt.1.and.enofe((i-1),j,ie).lt.1.0e-35)go to 188
      tloss=(gmax0(nd,i,j)-ggam(ie))/(tlavg*ggam(ie)*gmax0(nd,i,j))
      t2=dmin1(tloss,t2max)
      efill=1.0
      eterm1=1.0
      if(i.gt.1)eterm1=1.0d0-tlavgo*t1*ggam(ie)
c     Volume filling factor of e-s with ggam < gmin
      if(ie.lt.6)efill=(t2-(gmin-ggam(ie))/(tlavg*gmin*ggam(ie)))/delt
      eterm2=1.0d0-tlavg*t2*ggam(ie)
      if(eterm1.lt.0.0)eterm1=0.0
      if(eterm2.lt.0.0)eterm2=0.0
      enofe(i,j,ie)=n0(nd,i,j)/((sen-1.0)*tlavg*ggam(ie)**(sen+1.0))*
     ,  (eterm1**(sen-1.0)-eterm2**(sen-1.0))*efill
c      exxx=enofe(i,j,ie)
      if(enofe(i,j,ie).lt.1.0e-30)enofe(i,j,ie)=0.0
c     Divide by delt since integral is over time
      enofe(i,j,ie)=enofe(i,j,ie)/delt
c      exx2=enofe(i,j,ie)
  188 continue
      edist(ie)=enofe(i,j,ie)
c      if(itt.eq.itlast.and.i.lt.3)then
c        write(5,9987)ie,i,j,nd,it,ggam(ie),exxx,exx2, 
c     ,  eterm1,eterm2,delt,tloss,tlmin,tlavg,tlavgo,
c     ,  t1,t2,gmax0(nd,i,j),gmin,glow,gamb,n0(nd,i,j),bfij(nd,i,j),
c     ,  useedd,usp
c      end if
   90 continue
!$omp end parallel do
      ebk(nd,i,j)=egam(i,j,6)
      emold=0.0
c      if(itt.eq.itlast.and.i.lt.3)then
c      aeee=1.0
c        do 288 ie=1,44
c        if(ie.gt.1)aeee=alog10(enofe(i,j,ie)/enofe(i,j,ie-1))/
c     *    dlog10(ggam(ie)/ggam(ie-1))
c        write(5,9987)ie,i,j,nd,it,ggam(ie),enofe(i,j,ie),aeee
c  288    continue
c      end if
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(nd,i,j)
      gamd=gammad(nd,i,j)
      call polcalc(bfij(nd,i,j),bxfij(nd,i,j),byfij(nd,i,j),
     ,     bzfij(nd,i,j),clos,slos,chipol)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu,bperpp,sssab )
!$omp& schedule( dynamic, 1 )
      do 92 inu=1,64
      rnu(inu)=nu(inu)*zred1/delta(nd,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
      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
      bperpp=bperp
      emisco(inu)=ajnu(rnu(inu))*bperp*delta(nd,i,j)**2
c     Calculate emission coefficient with no beaming and reference B_perp
c     Replaces crude initial calculation
c      if(emisco(inu).ne.emisco(inu)) then
c      iccnt=iccnt+1
c
c     For test purposes:
c      if(inu.ne.17)go to 9516
c      if(i.eq.1.and.j.eq.20)write(5,9515)i,j,nd,delta(nd,i,j),
c     bperp,bfij(nd,i,j),bxfij(nd,i,j),byfij(nd,i,j),
c     , bzfij(nd,i,j),betadx(nd,i,j),betady(nd,i,j),betadz(nd,i,j),
c     ,  gammad(nd,i,j),delt,tlfact,gmax0(nd,i,j),gamb,ggam(36),
c     ,  edist(36),ggam(42),edist(42),ggam(44),edist(44)
c    ,(chipol*57.3),emisco(inu)
c 9516 continue
      bperpp=bpref
      ecofij(nd,i,j,inu)=ajnu(nu(inu))*bpref
      bperpp=bperp
      fsync2(inu)=emisco(inu)*zsize*(emfact*parsec)
      fsnoab(inu)=fsync2(inu)
c     Skip frequencies that are obviously optically thin
      if(inu.gt.16)go to 93
      sssab=akapnu(rnu(inu))
      if(sssab.gt.1.0e-30)ssabs(inu)=1.02e4*(sen+2.0)*sssab*bperp/
     ,   rnu(inu)*(parsec/rnu(inu))/delta(nd,i,j)
      ssabs(inu)=ssabs(inu)*zsize*delta(nd,i,j)
      bperpp=bpref
      acofij(nd,i,j,inu)=akapnu(nu(inu))
      bperpp=bperp
c     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)=xmaxx(i,j)*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
c      if(iccnt.gt.50)go to 9000
      fq1=0.98*rnu(1)
      emold=ajnu(fq1)*bperp*delta(nd,i,j)**2*zsize*(emfact*parsec)
      fq1=0.98*nu(1)
      do 94 inu=1,44
      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)then
        poldeg=3.0/(12.0*specin(inu)+19.0)
        chipol=chipol+pio2
      end if
      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(nd,i,j)
      gamd=gammad(nd,i,j)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 95 inu=17,64
      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(nd,i,j)**2/
     ,   dgpc**2*fgeom
      sscflx(inu)=ssc(rnu(inu))*3.086*volc*zred1*delta(nd,i,j)**2/
     ,   dgpc**2*fgeom
      ecij(nd,i,j,inu-16)=ecdust(nu(inu))
      sscij(nd,i,j,inu-16)=ssc(nu(inu))
      go to 88
   85 taupp(inu)=0.0
      ecij(nd,i,j,inu-16)=0.0
      sscij(nd,i,j,inu-16)=0.0
   88 continue
   95 continue
!$omp end parallel do
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu,nnu )
!$omp& schedule( dynamic, 1 )
      do 7795 inu=17,64
      nnu=inu-16
      alphec(nd,i,j,nnu)=0.0
      alphsc(nd,i,j,nnu)=0.0
      if(inu.eq.64) go to 7785
      if(ecij(nd,i,j,nnu).lt.1.0e-30.or.
     ,   ecij(nd,i,j,(nnu+1)).lt.1.0e-30) go to 7784
      alphec(nd,i,j,nnu)=-alog10(ecij(nd,i,j,(nnu+1))/
     ,   ecij(nd,i,j,(nnu)))/alog10(nu(inu+1)/nu(inu))
 7784 if(sscij(nd,i,j,nnu).lt.1.0e-30.or.
     ,   sscij(nd,i,j,(nnu+1)).lt.1.0e-30) go to 7785
      alphsc(nd,i,j,nnu)=-alog10(sscij(nd,i,j,(nnu+1))/
     ,   sscij(nd,i,j,nnu))/alog10(nu(inu+1)/nu(inu))
 7785 continue
 7795 continue
!$omp end parallel do
c     Approximate that highest frequency IC spectral index is same as for
c       next-highest frequency
      alphec(nd,i,j,48)=alphec(nd,i,j,47)
      alphsc(nd,i,j,48)=alphsc(nd,i,j,47)
      do 86 inu=14,64
      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(nd,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
c      write(*,*)'end calculation for it,i,j =',it,i,j
 1001 continue
  100 continue
c      write(*,*)'end loop over all cells, it =',it
c
c     End loop over all cells
c
      if(itt.ge.(itlast-1))then
        nd=ndmax
        do 5858 j=1,jcells
        do 5857 i=1,imax(j)
        write(1,9994)i,j,bfij(nd,i,j),n0(nd,i,j),enofe(i,j,35),
     ,    egam(i,j,43),nu(19),flsync(i,j,19),nu(30),flsync(i,j,30),
     ,    sscij(nd,i,j,14),alphsc(nd,i,j,14),ecij(nd,i,j,14),nu(46),
     ,    sscij(nd,i,j,30),alphsc(nd,i,j,30),ecij(nd,i,j,30)
 5857 continue
 5858 continue
      end if

      tflold=0.0
c     iwp determines whether map of emission will be written to a file
cccccc Set iwp=1 for recording map at all time steps; set iwp=0 otherwise
      iwp=1
      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
      do 500 inu=1,64
      tflux=0.0
      tsflux=0.0
      tcflux=0.0
      tecfl=0.0
      tsscfl=0.0
      alph=0.0
      qcum=0.0
      ucum=0.0
c     Do not make maps if nomap=1 was specified at the beginning
c      if(nomap.eq.1)iwp=0
      if(iwp.ne.1.or.inu.ne.nrmap)go to 7300
      if(nomap.eq.1.or.itt.le.0)go to 7300
      write(filnam,"(a8,i4.4,a4)") "temzmap",itt,".txt"
      open (7,iostat=ios, err=9000, file=dpath//filnam,status='new')
      write(filnm2,"(a8,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)
      poldgg=0.0
      pang=0.0
      if(flsync(i,j,inu).le.1.0e-6)go to 297
      if(iwp.ne.1.or.itt.le.0)go to 297
      if(inu.ne.nrmap.and.inu.ne.nopmap)go to 297
      if(nomap.eq.1)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(nd,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(nd,i,j),it
  297 continue
  298 continue
c  299 continue
  300 continue
      phots(inu)=0.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(tflux.gt.0.0.and.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 spectral energy distribution (SED) to file
      tfsync=1.0e3*tsflux/nu(inu)
c      write(3,9996)nu(inu),tflux,alph,tsflux,tecfl,tsscfl
      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.48)tfl48=tflux
      if(inu.eq.49)tfl49=tflux
      if(inu.eq.50)tfl50=tflux
      if(inu.eq.51)tfl51=tflux
      if(inu.eq.49)tfcomp=1.0e9*tcflux/nu(inu)
      if(inu.eq.49)tfec=tecfl
      if(inu.eq.49)tfssc=tsscfl
      if(inu.eq.11)tfl11=tflux
      if(inu.eq.6)tfl6=tflux
      tflold=1000.0*tflux/nu(inu)
c     Record linear polarization in perent and position angle in degrees at selected frequencies
      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)
c     Compute spectral index at selected frequencies
      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(tfl48.gt.0.0.and.tfl49.gt.0.0)alp49a=alog10(tfl48/tfl49)/
     ,alog10(nu(49)/nu(48))+1.0
      if(tfl49.gt.0.0.and.tfl50.gt.0.0)alp49b=alog10(tfl49/tfl50)/
     ,alog10(nu(50)/nu(49))+1.0
      if(tfl50.gt.0.0.and.tfl51.gt.0.0)alp50a=alog10(tfl50/tfl51)/
     ,alog10(nu(50)/nu(49))+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)
      alph49=0.5*(alp49a+alp49b)
      alph50=0.5*(alp49b+alp50a)
      timeo=itt*dtime
c     Write light curve points to file
c      write(4,9991)itt,timeo,nu(6),tfl6,nu(20),tfl20,
      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(50),tfl50,alph50
c     Write selected polarization data to file
c      write(5,9988)itt,timeo,nu(3),pdeg3,pang3,nu(4),pdeg4,
      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
      close(1)
      open (1,iostat=ios, err=9000, file='temzcheck.txt',
     , status='old')
      close(3)
      open (3,iostat=ios, err=9000, file='temzspec.txt',
     , status='old')
      close(4)
      open (4,iostat=ios, err=9000, file='temzlc.txt',
     , status='old')
      close(5)
      open (5,iostat=ios, err=9000, file='temzpol.txt',
     , status='old')
 1997 continue
      ispec=ispec-1
      if(ispec.eq.0)ispec=ispecs+1
c     For testing purposes:
c      ndhalf=0.5*ndmax
c      j=20
c      do 7599 i=imax(j),2,-1
c      nd=1
c      write(1,9517)itt,nd,i,bfij(nd,i,j),
c     ,  bxfij(nd,i,j),
c     ,  byfij(nd,i,j),bzfij(nd,i,j),bpp(nd,i,j),gammad(nd,i,j),
c     ,  betad(nd,i,j),delta(nd,i,j),gmax0(nd,i,j)
c      nd=ndhalf
c      write(1,9517)itt,nd,i,
c     ,  bfij(nd,i,j),bxfij(nd,i,j),
c     ,  byfij(nd,i,j),bzfij(nd,i,j),bpp(nd,i,j),gammad(nd,i,j),
c     ,  betad(nd,i,j),delta(nd,i,j),gmax0(nd,i,j)
c      nd=ndmax-1
c      write(1,9517)itt,nd,i,bfij(nd,i,j),
c     ,  bxfij(nd,i,j),
c     ,  byfij(nd,i,j),bzfij(nd,i,j),bpp(nd,i,j),gammad(nd,i,j),
c     ,  betad(nd,i,j),delta(nd,i,j),gmax0(nd,i,j)
c 7599 continue
c     Move cells in time and space arrays to make room for next time step
      do 599 j=2,jcells
      do 1599 i=imax(j),1,-1
      do 2599 nd=1,ndmax
      if(i.gt.1)then
      ussold(nd,i,j)=ussold(nd,i-1,j)
      gmax0(nd,i,j)=gmax0(nd,i-1,j)
      bfij(nd,i,j)=bfij(nd,i-1,j)
      bxfij(nd,i,j)=bxfij(nd,i-1,j)
      byfij(nd,i,j)=byfij(nd,i-1,j)
      bzfij(nd,i,j)=bzfij(nd,i-1,j)
      bpp(nd,i,j)=bpp(nd,i-1,j)
      n0(nd,i,j)=n0(nd,i-1,j)
      betad(nd,i,j)=betad(nd,i-1,j)
      betadx(nd,i,j)=betadx(nd,i-1,j)
      betady(nd,i,j)=betady(nd,i-1,j)
      betadz(nd,i,j)=betadz(nd,i-1,j)
      gammad(nd,i,j)=gammad(nd,i-1,j)
      delta(nd,i,j)=delta(nd,i-1,j)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 4599 inu=1,48
      if(inu.le.44)then
      ecofij(nd,i,j,inu)=ecofij(nd,i-1,j,inu)
      if(inu.le.16)acofij(nd,i,j,inu)=acofij(nd,i-1,j,inu)
      alphsy(nd,i,j,inu)=alphsy(nd,i-1,j,inu)
      endif
      ecij(nd,i,j,inu)=ecij(nd,i-1,j,inu)
      alphec(nd,i,j,inu)=alphec(nd,i-1,j,inu)
      sscij(nd,i,j,inu)=sscij(nd,i-1,j,inu)
      alphsc(nd,i,j,inu)=alphsc(nd,i-1,j,inu)
 4599 continue
!$omp end parallel do
      else
      if(nd.eq.ndmax)go to 2598
      ussold(nd,1,j)=ussold(nd+1,1,j)
      gmax0(nd,1,j)=gmax0(nd+1,1,j)
      bfij(nd,1,j)=bfij(nd+1,1,j)
      bxfij(nd,1,j)=bxfij(nd+1,1,j)
      byfij(nd,1,j)=byfij(nd+1,1,j)
      bzfij(nd,1,j)=bzfij(nd+1,1,j)
      bpp(nd,1,j)=bpp(nd+1,1,j)
      n0(nd,1,j)=n0(nd+1,1,j)
      betad(nd,1,j)=betad(nd+1,1,j)
      betadx(nd,1,j)=betadx(nd+1,1,j)
      betady(nd,1,j)=betady(nd+1,1,j)
      betadz(nd,1,j)=betadz(nd+1,1,j)
      gammad(nd,1,j)=gammad(nd+1,1,j)
      delta(nd,1,j)=delta(nd+1,1,j)
!$omp parallel do
!$omp& default( shared )
!$omp& private( inu )
!$omp& schedule( dynamic, 1 )
      do 5599 inu=1,48
      if(inu.le.44)then
      ecofij(nd,i,j,inu)=ecofij(nd+1,1,j,inu)
      if(inu.le.16)acofij(nd,i,j,inu)=acofij(nd+1,1,j,inu)
      alphsy(nd,i,j,inu)=alphsy(nd+1,1,j,inu)
      endif
      ecij(nd,i,j,inu)=ecij(nd+1,1,j,inu)
      alphec(nd,i,j,inu)=alphec(nd+1,1,j,inu)
      sscij(nd,i,j,inu)=sscij(nd+1,1,j,inu)
      alphsc(nd,i,j,inu)=alphsc(nd+1,1,j,inu)
 5599 continue
!$omp end parallel do
      endif
 2598 continue
 2599 continue
 1599 continue
  599 continue
c      for testing purposes:
c      j=20
c      do 4598 i=1,5
c      write(1,9517)itt,j,i,bpp(ndmax,i,j),bpp(ndmax-1,i,j),
c     ,  bpp(ndmax-2,i,j),bpp(ndmax-3,i,j),bpp(ndmax-4,i,j)
c 4598 continue
      write(*,*)'end of time step it.itt =',it,itt
      if(itt.eq.itlast)go to 9000
c      if(it.eq.2)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')
 6993 format('* ',7i5,1p15e10.3)
 7777 format(i4,1x,i4,2f12.4,2f8.2,2x,f10.4)
 8888 format(i5,f9.2)
 8889 format(/)
 9111 format(a10/i3,i8,i2,i1/f6.4/f7.3/f4.2/f4.2/f4.2/f5.3/f3.1/e7.1/
     ,f7.5/f7.1/f5.1/f7.1/d10.6/d9.5/d6.1/d6.1/d6.2/f6.1/f6.3/f4.2/
     ,f4.2/f4.2/f5.2/f5.2)
 9123 format(i5,9f10.5)
 9222 format('idelay out of bounds ',2i9,4i6,i10,1p11e11.3)
 9223 format('*idelay out of bounds ',2i12,3i6,1p7e11.3)
 9224 format('idelay ',2i9,4i6,i10,1p11e11.3)
 9232 format(2i4,3i6,1p12e10.2)
 9299 format(i12,2i5,2i6,3i10,1p5e12.4)
 9333 format(4i7,1p12e10.2)
 9515 format(' i j nd delta bp b(4) beta(3) chi tlfact gam N j ',3i4,
     ,       2f8.3,1p10e11.4)
 9517 format('itt i B(4) BP gam beta delta gmax0 ',3i5,10f11.3)
 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)
 9987 format('N ',5i5,1p20e12.4)
 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)')
 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.77205*x**0.293627-6.51657e-4/x**0.50372)/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.77205*x**0.293627-6.51657e-4/x**0.50372)/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,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,randn
      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),dat(65536),R(65536)
      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,nn,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,g2,
     ,  dotprd,vparx,vpary,vparz,vprpx,vprpy,vprpz,vprp2,vprp,uprp,
     ,  vd,gd,vd2,vdx,vdy,vdz,vfact,vdprpx,vdprpy,vdprpz,vdprp2,
     ,  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, it is not a shock
      if(eta.lt.1.0)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     front 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/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
      if(v.lt.1.0e-3)go to 10
      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
      go to 20
 10   dotprd=bx*ax+by*ay+bz*az
      apx=ax
      apy=ay
      apz=az
c     Calculate components of B field parallel and perpendicular to shock normal
 20   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(1,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/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,val1,val2,vala,valb,csth1,csth2,xsecc
      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
c   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
c 6666 format('In ecdust',i5,2x,1p18e12.3)
c 9996 format('a ',1p12e11.3)
c 9997 format('b ',1p12e11.3)
c 9998 format('c ',1p12e11.3)
c 9999 format('d ',1p12e11.3)
c 4000 write(10,6666)ie1,anuf,gam(1),edist(1),gam(44),edist(44),dnu(3),
c     , di(3),dnu(20),di(20),cscat,xsecc,addit,val2,valb,a,gran,ecdust
 4000 return
      end
c
c  Subroutine to calculate synchrotron self-Compton emission
c    of seed photons from other cells
      function ssc(anuf)
      common/cdist/gam,edist
      common/cssc/dnu,di,cscat
      common/crite/iwrite
      dimension gam(44),edist(44),dnu(64),di(64)
      real*8 gam,addit,gran,ssca,val1,val2,xsecc,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=64
      nuhi=64
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
c 4000 write(4,6666)ie1,anuf,gam(1),edist(1),gam(44),edist(44),dnu(14),
c     , di(14),dnu(15),di(15),cscat,xsecc,addit,val2,valb,a,gran,ssca
c 6666 format('In ssc',i5,1p18e11.3)
 4000 return
      end
