*/
*/     **************************************************************
*/     Patch file Version 1.0 to create version MCNP_CEAR from MCNP4b
*/                          by Lianyan Liu
*/         Center of Engineering Applications of Radioistopes
*/                North Carolina State University
*/                    Raleigh, North Carilina
*/                           May 1997
*/     **************************************************************
*/
*id mcliu
*i,mc.34
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                                                         
*d,mc.230
      if(iovr.ne.1) goto 195 
      call imcn        
      if(idum(1).eq.0) goto 195
      inquire(file='imp.in',exist=map_exist)
      if(map_exist) goto 180 
      write(*,*)
      write(iuo,*)
      write(*,*)  'Warning: input importance map does not exist'
      write(iuo,*)'Warning: input importance map does not exist'
      if(idum(1).eq.1) stop 'WW Map file is needed for this option'
      write(*,*)  'Message: start without fine-meshed WW'
      write(iuo,*)'Message: start without fine-meshed WW'         
  180 continue
      if(.not.map_exist) call newmesh
      if(.not.map_exist) goto 185
      call rdmap
      call get_wtexp
  185 continue
      if((.not.map_exist).or.(idum(1).eq.3)) call resetmap
      if(idum(1).ge.2.and.igww.eq.0) stop 'Need WWG card'                                     
  195 continue           
*d,mc.237,mc.246
      write(jtty,200)hovr                                               
  200 format(1x,a6,8h is done)                                          
      if(nwer.eq.1)write(iuo,205)                                       
  205 format(/10x,25h1 warning message so far.)                         
      if(nwer.gt.1)write(iuo,210)nwer                                   
  210 format(/i10,25h warning messages so far.)                         
      if(nfer.eq.1)write(iuo,215)                                       
  215 format(/10x,21h1 fatal error so far.)                             
      if(nfer.gt.1)write(iuo,220)nfer                                   
  220 format(/i10,21h fatal errors so far.)              
*/     Modifications to subroutine HSTORY
*id hsliu
*i,hs.4
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*i,hs.25
      if(idum(1).eq.0) goto 35
      if(uuu.ne.0.0) uinv=1./uuu
      if(vvv.ne.0.0) vinv=1./vvv
      if(www.ne.0.0) winv=1./www
   35 b=-log(rang())
c
*d,hs.33
   40 continue
      if(idum(1).ne.0) jsu=0
      if(lca(llca+icl).lt.0) call chkcel(icl,3,j)
*d,hs.77
      if(ple.eq.0.)go to 80   
*d,hs.80
      if(qpl.le.0.)go to 80  
*i,hs.89
c
   80 dx=huge
      dy=huge
      dz=huge
      if(idum(1).eq.0) goto 88
      if(ix.lt.1     .and.uuu.lt.0.0) goto 82
      if(ix.gt.nxmesh.and.uuu.gt.0.0) goto 82 
      if(uuu.gt.0.0) dx=(xmesh(ix+1)-xxx)*uinv
      if(uuu.lt.0.0) dx=(xmesh(ix)  -xxx)*uinv
   82 if(iy.lt.1     .and.vvv.lt.0.0) goto 84
      if(iy.gt.nymesh.and.vvv.gt.0.0) goto 84
      if(vvv.gt.0.0) dy=(ymesh(iy+1)-yyy)*vinv
      if(vvv.lt.0.0) dy=(ymesh(iy  )-yyy)*vinv
   84 if(iz.lt.1     .and.www.lt.0.0) goto 88
      if(iz.gt.nzmesh.and.www.gt.0.0) goto 88
      if(www.gt.0.0) dz=(zmesh(iz+1)-zzz)*winv
      if(www.lt.0.0) dz=(zmesh(iz  )-zzz)*winv
   88 dmesh=min(dx,dy,dz)    
*d,hs.91
      if(ple.ne.0.0.and.qpl.gt.0.0) pmf=b*gs 
*d,hs4b.15
   90 d=min(pmf,dls,dxl,dtc,deb,dw,dmesh)
*i,hs.124
      b=b-d*qpl
*d,hs4b.25
  170 if(d.ne.dw) goto 171
*i,hs4b.29
c
  171 if(idum(1).eq.0.or.d.ne.dmesh) goto 175
      if(d.eq.dx.and.uuu.gt.0.0) ix=ix+1
      if(d.eq.dx.and.uuu.lt.0.0) ix=ix-1
      if(d.eq.dy.and.vvv.gt.0.0) iy=iy+1
      if(d.eq.dy.and.vvv.lt.0.0) iy=iy-1
      if(d.eq.dz.and.www.gt.0.0) iz=iz+1
      if(d.eq.dz.and.www.lt.0.0) iz=iz-1
      dls=dls-d
      dxl=dxl-d
      dtc=dtc-d
      deb=deb-d
      call binary_search(ntmesh,tmesh,tme,it)
      if(idum(1).eq.1) goto 172
      if(ix.lt.1.or.ix.gt.nxmesh) goto 172
      if(iy.lt.1.or.iy.gt.nymesh) goto 172
      if(iz.lt.1.or.iz.gt.nzmesh) goto 172 
      if(ie.lt.1.or.ie.gt.(nemesh(1)+nemesh(2))) goto 172 
      if(it.lt.1.or.it.gt.ntmesh) goto 172
      nstop=nstop+1
      if(nstop.gt.maxstop) stop 'nstop reaches the maximum'
      nxtrack(nstop)=ix
      nytrack(nstop)=iy
      nztrack(nstop)=iz
      netrack(nstop)=ie
      nttrack(nstop)=it
      rflux(ix,iy,iz,ie,it)=rflux(ix,iy,iz,ie,it)+wgt
  172 if(d.eq.dls) goto 178
      if(.not.map_exist) goto 80
      call wtwndo(4,ww)
      if(nter.ne.0) goto 260
      goto 80
c                              
*d,hs4b.32
  175 if(d.ne.dls)go to 180                                             
  178 continue
*i,hs.178
      jsup=0
*i,hs.188
      if(idum(1).eq.0) goto 195
      call binary_search(ntmesh,tmesh,tme,it)
      ie=0
      if(nemesh(ipt).eq.0) goto 195      
      call binary_search(nemesh(ipt),emesh(1,ipt),erg,ie)
      if(ipt.eq.1.and.ie.ne.0) ie=ie+nemesh(2)
  195 continue
*d,hs4b.35,hs4b.36
  200 jsup=0
      if(abs(wwp(ipt,4)).eq.1..and.(idx.eq.0.or.wwp(ipt,5).lt.0.))      
     1 call wtwndo(4,ww)
      if(idum(1).ne.0.and.map_exist) call wtwndo(4,ww)   
*id spliu
*i,sp.4
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*d,sp.113
  220 if(idum(1).eq.0) goto 225 
      nstop=0 
      frmsrc=frmsrc+wgt
      call binary_search(nxmesh,xmesh,xxx,ix)
      call binary_search(nymesh,ymesh,yyy,iy)
      call binary_search(nzmesh,zmesh,zzz,iz)
      ie=0
      if(nemesh(ipt).eq.0) goto 225      
      call binary_search(nemesh(ipt),emesh(1,ipt),erg,ie)
      if(ipt.eq.1.and.ie.ne.0) ie=ie+nemesh(2)
  225 paxtc(2,1,ipt)=paxtc(2,1,ipt)+wgt          
*id wwliu
*i,ww.5
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*i,ww.6
      if(idum(1).ne.0) goto 25
c
*i,ww.17
c   
   25 inside=1
      if(ix.lt.1.or.ix.gt.nxmesh) inside=0
      if(iy.lt.1.or.iy.gt.nymesh) inside=0
      if(iz.lt.1.or.iz.gt.nzmesh) inside=0
      if(ie.lt.1.or.ie.gt.(nemesh(1)+nemesh(2))) inside=0
      if(it.lt.1.or.it.gt.ntmesh) inside=0
      if(inside.eq.0) goto 50
c
*d,ww4b.2
   30 if(idum(1).eq.0) ww=wwval(ipt,icl,ib,1)
      if(idum(1).gt.0) ww=0.50*wtexp(ix,iy,iz,ie,it) 
*d,ww.22
      if(idum(1).eq.0) t1=wwp(ipt,1)*ww
      if(idum(1).gt.0) t1=4.0*ww 
*d,ww.26
      if(idum(1).eq.0) npa=min(wwp(ipt,3)-1.,wgt/t1)  
      if(idum(1).gt.0) npa=min(5.0d0,wgt/(2.*ww)-1.0)   
*d,ww.34
   40 if(idum(1).eq.0) t2=min(wgt*wwp(ipt,3),ww*wwp(ipt,2))  
      if(idum(1).gt.0) t2=ww*2.0  
*d,ww.40,ww.41
      if(idum(1).eq.0.or.inside.eq.0) t3= fiml9(1,1)/fiml(1)       
      if(idum(1).eq.0.or.inside.eq.0) goto 55
      t4=geimp(ix,iy,iz,ie)
      if(t4.eq.0.0) t4=gimp(ix,iy,iz)
      if(t4.eq.0.0) stop 'Zero geometry importance in WTWNDO'    
      t3= 1./t4    
   55 if(wgt.gt.t3*wcs2tc(ipt)) return                   
      t2=wcs1tc(ipt)*t3    
*id cnliu
*i,cn.4
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*i,cn4b.15
      if(idum(1).ne.0.and.map_exist) goto 160
*id cpliu
*i,cp.4
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*i,cp4b.17
      if(idum(1).ne.0.and.map_exist) goto 70
*id wgliu
*i,wg.5
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
*i,wg.19
c
   50 if(idum(1).eq.0.or.m.eq.1) return 
      rspns=rspns+wt                                                              
      do 80 n=1,nstop
      i1=nxtrack(n)
      i2=nytrack(n)
      i3=nztrack(n)
      i4=netrack(n)
      i5=nttrack(n)
      cflux(i1,i2,i3,i4,i5)=cflux(i1,i2,i3,i4,i5)+wt
  80  continue 
c
*id ouliu
*d,ou.46
      if(idum(1).le.1) goto 65  
      call get_wtexp
      call wtmap
      if(idum(1).eq.3) call resetmap                                               
   65 tlc=cp1                           
*af,ou
*deck bsliu
      subroutine binary_search(ndim,xarray,x,index)
c
*if def,cheap,1  
      implicit double precision (a-h,o-z)   
      dimension xarray(1)
c
      index=0
      if(x.lt.xarray(1).or.x.gt.xarray(ndim+1)) return
c
      ileft=1
      iright=ndim
  10  imid=(ileft+iright)/2
      if(imid.ge.1.and.imid.le.ndim) goto 15 
      write(*,*)'The array=',(xarray(i),i=1,ndim)
      write(*,*)'ileft=',ileft,'iright=',iright,'imid=',imid
      write(*,*)'x value=',x 
      stop 'Error in computing IMID'
  15  continue
      if(x.ge.xarray(imid).and.x.le.xarray(imid+1)) goto 20
      if(x.lt.xarray(imid)) iright=imid-1
      if(x.gt.xarray(imid+1)) ileft=imid+1
      goto 10
c
  20  index=imid 
      return
      end
*deck gwliu
      subroutine get_wtexp
c
*call cm
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
c
      dimension cfluxt(maxnx+maxny+maxnz),rfluxt(maxnx+maxny+maxnz)
c
      if(frmsrc.eq.0.0) stop 'No source particles are simulated.'
      if(rspns.eq.0.0)  stop 'No particles ever score. Bad WW map'
c
      if(.not.x_symm) goto 50
      do 40 i5=1,ntmesh
      do 40 i4=1,(nemesh(1)+nemesh(2))
      do 40 i3=1,nzmesh
      do 40 i2=1,nymesh
      do 20 i1=1,nxmesh
      cfluxt(i1)=cflux(i1,i2,i3,i4,i5)
      rfluxt(i1)=rflux(i1,i2,i3,i4,i5)
  20  continue
      do 30 i1=1,nxmesh
      cflux(i1,i2,i3,i4,i5)=0.5*(cfluxt(i1)+cfluxt(nxmesh+1-i1))
      rflux(i1,i2,i3,i4,i5)=0.5*(rfluxt(i1)+rfluxt(nxmesh+1-i1))
  30  continue
  40  continue
c
  50  if(.not.y_symm) goto 90
      do 80 i5=1,ntmesh 
      do 80 i4=1,(nemesh(1)+nemesh(2))
      do 80 i3=1,nzmesh
      do 80 i1=1,nxmesh
      do 60 i2=1,nymesh
      cfluxt(i2)=cflux(i1,i2,i3,i4,i5)
      rfluxt(i2)=rflux(i1,i2,i3,i4,i5)
  60  continue
      do 70 i2=1,nymesh
      cflux(i1,i2,i3,i4,i5)=0.5*(cfluxt(i2)+cfluxt(nymesh+1-i2))
      rflux(i1,i2,i3,i4,i5)=0.5*(rfluxt(i2)+rfluxt(nymesh+1-i2))
  70  continue
  80  continue
c
 90   if(.not.z_symm) goto 130
      do 120 i5=1,ntmesh
      do 120 i4=1,(nemesh(1)+nemesh(2))
      do 120 i2=1,nymesh
      do 120 i1=1,nxmesh
      do 100 i3=1,nzmesh
      cfluxt(i3)=cflux(i1,i2,i3,i4,i5)
      rfluxt(i3)=rflux(i1,i2,i3,i4,i5)
  100 continue
      do 110 i3=1,nzmesh
      cflux(i1,i2,i3,i4,i5)=0.5*(cfluxt(i3)+cfluxt(nzmesh+1-i3))
      rflux(i1,i2,i3,i4,i5)=0.5*(rfluxt(i3)+rfluxt(nzmesh+1-i3))
  110 continue
  120 continue
c
  130 reference=frmsrc/rspns
c
      do 160  i5=1,ntmesh
      do 160  i4=1,(nemesh(1)+nemesh(2))
      do 160  i3=1,nzmesh
      do 160  i2=1,nymesh
      do 160  i1=1,nxmesh
      t=cflux(i1,i2,i3,i4,i5)
      wtexp(i1,i2,i3,i4,i5)=0.0
      if(t.gt.0.0) wtexp(i1,i2,i3,i4,i5)=rflux(i1,i2,i3,i4,i5) 
     $                                   /(t*reference)
 160  continue
c
      do 180 i3=1,nzmesh
      do 180 i2=1,nymesh
      do 180 i1=1,nxmesh
      t1=0.0
      t2=0.0
      do 170 i4=1,(nemesh(1)+nemesh(2))
      t3=0.0
      t4=0.0
      do 165 i5=1,ntmesh
      t1=t1+cflux(i1,i2,i3,i4,i5)
      t2=t2+rflux(i1,i2,i3,i4,i5)
      t3=t3+cflux(i1,i2,i3,i4,i5)
      t4=t4+rflux(i1,i2,i3,i4,i5)
 165  continue
      geimp(i1,i2,i3,i4)=1.0
      if(t3.gt.0.0.and.t4.gt.0.0) geimp(i1,i2,i3,i4)=t3/t4*reference
 170  continue
      gimp(i1,i2,i3)=1.0
      if(t1.gt.0.0.and.t2.gt.0.0) gimp(i1,i2,i3)=t1/t2*reference
 180  continue
c
      nzero=0
      do 195 n=1,nxmesh+nymesh+nzmesh
      do 190 i3=1,nzmesh
      do 190 i2=1,nymesh
      do 190 i1=1,nxmesh
      if(gimp(i1,i2,i3).gt.0) goto 190
      nzero=nzero+1
      t1=0.0
      if(i1.gt.1)      t1=t1+gimp(i1-1,i2,i3)
      if(i1.lt.nxmesh) t1=t1+gimp(i1+1,i2,i3)
      if(i2.gt.1)      t1=t1+gimp(i1,i2-1,i3)
      if(i2.lt.nymesh) t1=t1+gimp(i1,i2+1,i3)
      if(i3.gt.1)      t1=t1+gimp(i1,i2,i3-1)
      if(i3.lt.nzmesh) t1=t1+gimp(i1,i2,i3+1)
      if(t1.eq.0.0)  gimp(i1,i2,i3)=1.0
      if(t1.gt.0.0)  gimp(i1,i2,i3)=t1/6.0
c      if(t1.gt.0.0)   gimp(i1,i2,i3)=t1/60.0
 190  continue 
      if(nzero.eq.0) goto 220
      nzero=0
 195  continue
c
      do 200 i3=1,nzmesh
      do 200 i2=1,nymesh
      do 200 i1=1,nxmesh
      if(gimp(i1,i2,i3).le.0.0) 
     $      stop 'Zero geometric importace in get_wtexp'
  200 continue
c
 220  map_exist=.true.
c
      return
      end
*deck nmliu
      subroutine newmesh
c
*call cm
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c
c     get the number of uniform regions for position and time
      nxreg=   idum(2)
      nyreg=   idum(3+abs(nxreg))
      nzreg=   idum(4+abs(nxreg)+abs(nyreg))
      ntreg=   idum(5+abs(nxreg)+abs(nyreg)+abs(nzreg))
c
      x_symm=.false.
      y_symm=.false.
      z_symm=.false.
      if(nxreg.lt.0) x_symm = .true.
      if(nyreg.lt.0) y_symm = .true.
      if(nzreg.lt.0) z_symm = .true.
      if(x_symm) nxreg=-nxreg
      if(y_symm) nyreg=-nyreg
      if(z_symm) nzreg=-nzreg
c
      nxmesh=0
      do 10 n=1,nxreg
      nxmesh=nxmesh      +idum(2+n)
   10 continue
      nymesh=0
      do 15 n=1,nyreg
      nymesh=nymesh      +idum(3+nxreg+n)
   15 continue
      nzmesh=0
      do 20 n=1,nzreg
      nzmesh=nzmesh      +idum(4+nxreg+nyreg+n)
   20 continue
      ntmesh=0
      do 25 n=1,ntreg
      ntmesh=ntmesh      +idum(5+nxreg+nyreg+nzreg+n)
   25 continue
c
      i=1
      do 40 n=1,nxreg
      xleft=rdum(n)
      xright=rdum(n+1)
      imesh=idum(2+n)
      dx=(xright-xleft)/float(imesh)
      do 40 m=1,imesh
      xmesh(i)=xleft+dx*float(m-1)
      i=i+1
   40 continue
      xmesh(i)=xmesh(i-1)+dx
c
      if(.not.x_symm) goto 50
      nxmesh=2*nxmesh
      xmid=xmesh(1)
      do i=2,nxmesh/2+1
      xmesh(i+nxmesh/2)=xmesh(i)
      enddo
      do i=1,nxmesh/2
      d=xmesh(nxmesh+2-i)-xmid
      xmesh(i)=xmid-d
      enddo
      xmesh(nxmesh/2+1)=xmid
c
   50 i=1
      do 60 n=1,nyreg
      yleft =rdum(n+nxreg+1)
      yright=rdum(n+nxreg+2)
      imesh=idum(3+nxreg+n)
      dy=(yright-yleft)/float(imesh)
      do 60 m=1,imesh
      ymesh(i)=yleft+dy*float(m-1)
      i=i+1
   60 continue
      ymesh(i)=ymesh(i-1)+dy
c
      if(.not.y_symm) goto 65
      nymesh=2*nymesh
      ymid=ymesh(1)
      do i=2,nymesh/2+1
      ymesh(i+nymesh/2)=ymesh(i)
      enddo
      do i=1,nymesh/2
      d=ymesh(nymesh+2-i)-ymid
      ymesh(i)=ymid-d
      enddo
      ymesh(nymesh/2+1)=ymid
c
   65 i=1
      do 70 n=1,nzreg
      zleft =rdum(n+nxreg+nyreg+2)
      zright=rdum(n+nxreg+nyreg+3)
      imesh=idum(4+nxreg+nyreg+n)
      dz=(zright-zleft)/float(imesh)
      do 70 m=1,imesh
      zmesh(i)=zleft+dz*float(m-1)
      i=i+1
   70 continue 
      zmesh(i)=zmesh(i-1)+dz
c
      if(.not.z_symm) goto 75
      nzmesh=2*nzmesh
      zmid=zmesh(1)
      do i=2,nzmesh/2+1
      zmesh(i+nzmesh/2)=zmesh(i)
      enddo
      do i=1,nzmesh/2
      d=zmesh(nzmesh+2-i)-zmid
      zmesh(i)=zmid-d
      enddo
      zmesh(nzmesh/2+1)=zmid
c
   75 if(ntmesh.eq.0) goto 85
      i=1
      do 80 n=1,ntreg
      tleft =rdum(n+nxreg+nyreg+nzreg+3)
      tright=rdum(n+nxreg+nyreg+nzreg+4)
      imesh=idum(5+nxreg+nyreg+nzreg+n)
      dt=(tright-tleft)/float(imesh)
      do 80 m=1,imesh
      tmesh(i)=tleft+dt*float(m-1)
      i=i+1
   80 continue
      tmesh(i)=tmesh(i-1)+dt
c
  85  do 90 i=1,2
      nemesh(i)=ngww(i)
  90  continue
      do 110 i=1,2
      emesh(1,i)=0.0
      do 100 n=1,nemesh(i)
      emesh(n+1,i)=ewwg(lewg+n+mgww(i))
  100 continue
  110 continue
c   
      write(iuo,120)
  120 format(/6x,42h The new mesh structure for importance map)
      write(iuo,130)nxmesh
  130 format(i4, 9h X meshes)
      write(iuo,140)(xmesh(i),i=1,nxmesh+1)
  140 format(5(1x,e12.4,1x))
      write(iuo,150)nymesh
  150 format(i4, 9h Y meshes)
      write(iuo,140)(ymesh(j),j=1,nymesh+1)
      write(iuo,160)nzmesh
  160 format(i4, 9h Z meshes)
      write(iuo,140)(zmesh(k),k=1,nzmesh+1)
      write(iuo,170)nemesh(1)
  170 format(i4, 20h E meshes of neutron)
      write(iuo,140)(emesh(l,1),l=1,nemesh(1)+1)
      write(iuo,180)nemesh(2)
  180 format(i4, 22h E meshes of gamma-ray)
      write(iuo,140)(emesh(l,2),l=1,nemesh(2)+1)
      write(iuo,190)ntmesh
  190 format(i4, 12h time meshes)
      write(iuo,140)(tmesh(m),m=1,ntmesh+1)
      write(iuo,*)
c
      if(nxmesh               .gt.maxnx) stop 
     $        'Mesh # of X exceeds maximum allowed'
      if(nymesh               .gt.maxny) stop 
     $        'Mesh # of Y exceeds maximum allowed'
      if(nzmesh               .gt.maxnz) stop 
     $        'Mesh # of Z exceeds maximum allowed'
      if((nemesh(1)+nemesh(2)).gt.maxne) stop 
     $        'Mesh # of E exceeds maximum allowed'
      if(ntmesh               .gt.maxnt) stop 
     $        'Mesh # of T exceeds maximum allowed'
      if(nemesh(1).eq.0.and.nemesh(2).eq.0) stop
     $        'Error -- both energy groups are zero'
c      
      do i=1,nxmesh
      if(xmesh(i).gt.xmesh(i+1)) stop 'Strange X-mesh'
      enddo
      do i=1,nymesh
      if(ymesh(i).gt.ymesh(i+1)) stop 'Strange Y-mesh'
      enddo  
      do i=1,nzmesh
      if(zmesh(i).gt.zmesh(i+1)) stop 'Strange Z-mesh'
      enddo
      do i=1,nemesh(1)
      if(emesh(i,1).gt.emesh(i+1,1)) stop 'Strange E-mesh of neutron'
      enddo
      do i=1,nemesh(2)
      if(emesh(i,2).gt.emesh(i+1,2)) stop 'Strange E-mesh of gamma-ray'
      enddo   
      do i=1,ntmesh
      if(tmesh(i).gt.tmesh(i+1)) stop 'Strange time mesh'
      enddo
c
      return
      end
*deck rmliu
      subroutine rdmap
c
*call cm
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
c
      open(1,file='imp.in',status='old',form='unformatted')
      read(1)x_symm,y_symm,z_symm
      read(1)nxmesh,nymesh,nzmesh,nemesh(1),nemesh(2),ntmesh
      read(1)(xmesh(i1),  i1=1,nxmesh+1)
      read(1)(ymesh(i2),  i2=1,nymesh+1)
      read(1)(zmesh(i3),  i3=1,nzmesh+1)
      read(1)(emesh(i4,1),i4=1,nemesh(1)+1)
      read(1)(emesh(i4,2),i4=1,nemesh(2)+1)  
      read(1)(tmesh(i5),  i5=1,ntmesh+1)
      read(1)frmsrc,rspns 
      read(1)(((((rflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $         i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)
      read(1)(((((cflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $         i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)
      close(1)
c      
      write(iuo,85)
  85  format(/4x,46h The read-in mesh structure for importance map)
      write(iuo,90)nxmesh
  90  format(i4, 9h X meshes)
      write(iuo,100)(xmesh(i1),i1=1,nxmesh+1)
 100  format(5(1x,e12.4,1x))
      write(iuo,110)nymesh
 110  format(i4, 9h Y meshes)
      write(iuo,100)(ymesh(i2),i2=1,nymesh+1)
      write(iuo,120)nzmesh
 120  format(i4, 9h Z meshes)
      write(iuo,100)(zmesh(i3),i3=1,nzmesh+1)
      write(iuo,130)nemesh(1)
 130  format(i4, 20h E meshes of neutron)
      write(iuo,100)(emesh(i4,1),i4=1,nemesh(1)+1)
      write(iuo,140)nemesh(2)
 140  format(i4, 22h E meshes of gamma-ray)
      write(iuo,100)(emesh(i4,2),i4=1,nemesh(2)+1)
      write(iuo,150)ntmesh
 150  format(i4, 12h Time meshes)
      write(iuo,100)(tmesh(i5),i5=1,ntmesh+1)
      write(iuo,*)
c
      return
      end
*deck wmliu
      subroutine wtmap
c 
*call cm
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
c
      open(1,file='imp.out',form='unformatted')
      rewind(1)
c
      write(1)x_symm,y_symm,z_symm
      write(1)nxmesh,nymesh,nzmesh,nemesh(1),nemesh(2),ntmesh
      write(1)(xmesh(i1),i1=1,nxmesh+1)
      write(1)(ymesh(i2),i2=1,nymesh+1)
      write(1)(zmesh(i3),i3=1,nzmesh+1)
      write(1)(emesh(i4,1),i4=1,nemesh(1)+1)
      write(1)(emesh(i4,2),i4=1,nemesh(2)+1)
      write(1)(tmesh(i5),i5=1,ntmesh+1)
      write(1)frmsrc,rspns 
      write(1)(((((rflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $      i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)
      write(1)(((((cflux(i1,i2,i3,i4,i5),i1=1,nxmesh),i2=1,nymesh),
     $      i3=1,nzmesh),i4=1,nemesh(1)+nemesh(2)),i5=1,ntmesh)
      close(1)
      return
      end
*deck rpliu
      subroutine resetmap
*if def,cheap,1  
      implicit double precision (a-h,o-z)  
c
      logical map_exist,x_symm,y_symm,z_symm
      dimension iwwmesh(6)
      equivalence (spare,iwwmesh)
      equivalence (iwwmesh(1),ix)
      equivalence (iwwmesh(2),iy)
      equivalence (iwwmesh(3),iz)
      equivalence (iwwmesh(4),ie)
      equivalence (iwwmesh(5),it)
      equivalence (iwwmesh(6),nstop)
      parameter (maxne=6,maxnx=30,maxny=30,maxnz=30,
     $           maxnt=4,maxstop=6000)
      common /liumap/rflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               cflux(maxnx,maxny,maxnz,maxne,maxnt),
     $               wtexp(maxnx,maxny,maxnz,maxne,maxnt),
     $               geimp(maxnx,maxny,maxnz,maxne),
     $               gimp (maxnx,maxny,maxnz),
     $               emesh(maxne+1,2),xmesh(maxnx+1),ymesh(maxny+1),
     $               zmesh(maxnz+1),  tmesh(maxnt+1),
     $               nemesh(2),nxmesh,nymesh,nzmesh,ntmesh,frmsrc,rspns,
     $               x_symm,y_symm,z_symm
      common /liutrack/nxtrack(maxstop),nytrack(maxstop),
     $                 nztrack(maxstop),netrack(maxstop),
     $                 nttrack(maxstop)
      common /inputmap/map_exist 
c                            
c
      frmsrc=0.0
      rspns=0.0
      do 100 i5=1,ntmesh
      do 100 i4=1,(nemesh(1)+nemesh(2))
      do 100 i3=1,nzmesh
      do 100 i2=1,nymesh
      do 100 i1=1,nxmesh
      rflux(i1,i2,i3,i4,i5)=0.0
      cflux(i1,i2,i3,i4,i5)=0.0
 100  continue  
c
      return
      end