interface.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	xbgrun
subroutine	watchdog (h)
real *8 function	getfai (x, y)
subroutine	det2web (x, y, xs, ys)

Function/Subroutine Documentation

det2web()

subroutine det2web	(	real*8	x,
		real*8	y,
		real*8	xs,
		real*8	ys
	)

Definition at line 789 of file interface.f.

Referenced by xbgrun().

      implicit none 
#include "Zmaxdef.h"
#include "Zobs.h"
#include "Zprivate.h"
      real*8 x, y ! input.  x,y component of the coordinate
                  !   or direction cos  in the detctor system.
      real*8 xs,ys ! ouput. values transformed to web system
               !  where the x-axis is directed to the incident
               !  xs,ys can be x,y
      
!      Z* = Z * exp(-Fai0) = (x + iy)(cos-i sin)
!                          = x cos  + y sin + i ( y cos -x sin)
!
      real*8 temp
      temp = x*cosrot + y*sinrot
      ys = y*cosrot - x*sinrot
      xs = temp

Here is the caller graph for this function:

getfai()

real*8 function getfai	(	real*8	x,
		real*8	y
	)

Definition at line 782 of file interface.f.

      implicit none
#include "Zglobalc.h"
!        get polar angle of the point (x,y) in deg.
      real*8 x, y  ! input postion or direction-cos components
      getfai = atan2(y, x)*todeg

watchdog()

subroutine watchdog

(

type(histogram1)

h

)

Definition at line 766 of file interface.f.

      implicit none
#include "Zmaxdef.h"
#include "Zobs.h"
#include "Zprivate.h"
#include "Zprivate2.h"
      type(histogram1):: h
      integer i
      write(0,*) ' h%c%init=', h%c%init
      write(0,*) ' h%c%title=', h%c%title
      if( h%c%init .eq. 'initend' ) then
         do i = 1, 20
            write(0,*) 'i=',i, ' h%dnw(i)=',h%dnw(i)
         enddo
      endif

xbgrun()

subroutine xbgrun

(

)

Definition at line 58 of file interface.f.

References b, bin, cdedxinair(), charge, cmintime2websec(), cmudedx(), cos, cqincident(), csortstack(), d, d0, depth, det2web(), false, height, histdep(), howmuch(), kelec, kgetenv2(), kmuon, kphoton, kseblk(), kwhist(), kwhistai(), kwhistc(), kwhistd(), kwhistdir(), kwhistev(), kwhisti(), kwhistid(), kwhistp(), kwhists(), kwhistso(), mass, mpirank, ne, nfai, ng, nmu, nrbin, p, rndc(), subcode, t, tkarspec, tkrtspec, true, w2hl(), x, xyz, and y.

      implicit none
#include "Zmaxdef.h"
#include "Zglobalc.h"
#include "Zmanagerp.h"
#include "Ztrack.h"
#include "Ztrackv.h"
#include "Ztrackp.h"
#include "Zcode.h"
#include "Zheavyp.h"
#include "Zobs.h"
#include "Zobsp.h"
#include "Zobsv.h"
#include  "Zstackv.h"
#include  "Zincidentv.h"
#include "Zprivate.h"
#include "Zprivate1.h"
#include "Zprivate2.h"
#include "Zprivate3.h"     
#include "Zprivate4.h" 
#if defined (DOMPI)
#include "mpif.h"
#include "Zmpi.h"
      real*8  ager(maxmpisize)
      real*8  esizer(maxmpisize) 
      real*8  wage, wsize
      real*8  etime1, etime2 
      character*10 numbr
      integer lengnr, rank, rankc
!         to gatter sum of the  all rank data 
      real nrfaireca(nrbin, nfai, 4, nsites)
      real nrfaialla(nrbin, nfai, 4, nsites)
      real derfaia(nrbin, nfai,  nsites)
      real*8 nga(maxnoofassites)
      real*8 nea(maxnoofassites)
      real*8 nmua(maxnoofassites)
      real*8 nhada(maxnoofassites)
      real*8 sumelossa(maxnoofassites)
#define  REALLIMITg  15000
#define  REALLIMITe  9000
#define  REALLIMITmu 7500
#define  REALLIMITh  7500
#endif
      type(track):: incident
      type(coord):: angleatobs

      logical  heobs            ! if T, currently observing 
      common /zheobs/ heobs     !  particle is the one  obsrved at skeelton making time



      integer id  ! input.  1 ==> aTrack is going out from                      
!                                 outer boundery.
!                           2 ==> reached at an observation level
!                           3 ==> reached at inner boundery.
      type(track):: atrack

      type(track):: inci
      type(coord):: angle
      type(coord):: tetafai
      
      character*128 input
      character*64 dirstr
      real sr, dr
      integer i, j, k, icon, ir, l, m
      integer nn
      integer iij, code, codex
      integer i1, i2, ic, ifai, ji

      real*8 r, eloss, rinmu, cosang
      real*8 webmin(nrbin, nfai, nsites)
      real*8 dedt, rho, dist, disto, dedxmu
      real*8 fai, rminm

      real*8 u
      logical accept

      real*8 wx, wy, wz
      real   za
      real  de, ek, f, molu
      real*8  cvh2den
      integer ldep, ridx, faiidx, depidx
      integer lengenv
      integer  ncpu, mcpu ! no. of smashed skeletons, and actully used skeletonsn
      integer  margin
      real*4 wgt, wwgt
      real*4 enhance0      ! since we use only mcpu, the result must be enhanced 
                          ! by a factor of ncpu/mcpu
      real*4 enhance      ! enhance0 * wgt
      integer binw
      character*9 ptcl(4)
      data ptcl/"Photons", "Electrons","Muons", "hadrons"/
      character*9 ptcl2(3)
      data ptcl2/"Electrons", "Muons","All"/
      real power(4)
      integer nstr
      real e0,  cosz, limit(4), limite(4)
      data power/1.,1.,1.,1./
      real  power2(4)
      data power2/1.,1.,1.,1./
      character*128 title
      character*96 evid(nsites), plotid
      real*8 cog, cog2, sumne,  obstimes, savederg(5)
      real*8 firstz, dd
      real*8 fai0
      real*8 getfai
      logical dosort
      real prob, tmin, dt
!     ***********************
#include "interface2.h"
!     *********************
!     example
!       histdep:  2 5 6 7 10 /
!       depth   1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
!       ansites = 5
!       w2hl: 0 1 0 0 2 3 4 0 0 5  0 0 ...
      


      lengenv = kgetenv2("LIMIT", input) 
      read(input(1:lengenv),*)  limit
      lengenv = kgetenv2("BINW", input) 
      read(input(1:lengenv),*)  binw
      lengenv = kgetenv2("NCPU", input) 
      read(input(1:lengenv),*)  ncpu
      lengenv = kgetenv2("MCPU", input)
      read(input(1:lengenv),*)  mcpu
      lengenv = kgetenv2("MARGIN", input)
      read(input(1:lengenv),*)  margin
      lengenv = kgetenv2("SeeLowdE", input)
      if(lengenv .le. 0)  then
         write(0,*) ' SeeLowdE not given'
         stop
      endif
      seelowde = input(1:lengenv) .eq. "yes"
      lengenv = kgetenv2("KeepWeight", input)
      keepweight = input(1:lengenv) .eq. "yes"

      enhance0 = ncpu
      enhance0 = enhance0/mcpu
      do i = 1, 4
         limite(i) = limit(i) * enhance0 
               !     this must be enhanced from the firs.
      enddo


      write(0,*) ' *** ncpu=',ncpu,
     *            ' mcpu=',mcpu, ' enhance factor=',enhance0
      write(0,*) ' margin=',margin
      do i = 1, noofsites
         w2hl(i) = 0
         w2il(i) = 0
      enddo
      do i = 1, noofsites
         if(histdep(i) .eq. 0) exit
         hnsites = i
         w2hl( histdep(i)) = i
      enddo

      do i = 1, noofsites
         if( indivdep(i) .eq. 0) exit
         ansites = i
         w2il( indivdep(i) ) = i
      enddo

!        now hnsites is the  actual number of sites where statistics
!        by histogram is needed.
!        mapping; depth index to array index; where2loc
!     

 
      call kwhistso( binw )  ! specify binary or ascii write of histogaram
                             !  1--> ascii  2--> binary


      rminm = rmin/10.**(bin/2.0d0)
      r=rmin
      dr = 10.**bin 
      rbin(1) = 0.
      do i = 2, nrbin
         rbin(i) = r
         r = r* dr
      enddo

      return      
!    ******************
      entry xihist
!
!

!     histogram: instanciate
!         rspec (lateral):  
!       

      if(tkarspec) then
         do i = 1, hnsites
            do k = 1, 4
               do j = 1, nfai
                  call kwhisti(rspec(k, j,  i),
     *                 rmin, bin, nrbin,  b'00011' )
                  write(plotid, 
     *            '("Lateral dist. of ",a," at given fai")')
     *                 ptcl(k)
                  call kwhistai(rspec(k, j, i), 
     *            plotid,  "ar", "ptcls", .true.,  power(k), 
     *            "r", "m.u")
               enddo
            enddo
         enddo
      endif
         
      if(tkrtspec) then
          do i = 1, hnsites
             l = histdep(i)
             call cmintime2websec(obssites(l).pos.xyz,
     *       l, i,   webmin)
          enddo

          do i = 1, hnsites
             do j = 1, 4
!                at  center
                call kwhisti( tspec0(j, i),
     *               -5., 0.1, 200, b'00000')
                call kwhistai( tspec0(j, i),
     *          "Arrival time dist. of "//ptcl(j)//" at center",
     *          "t", "ptcls", .false., 0., "time", "ns")

                do ir=1, nrbin
                   do ifai=1, nfai
                      tmin = webmin(ir, ifai,i)
                      dt = 0.01*10.0**(bin*(ir-1))*100.  ! approx core distnace m
                      dt = dt**0.675*1.e9/3.0e8/30.     ! if sqrt 1m-->0.1 ns 10 m-->0.3 ns
                                       ! 100m 1ns 1km 3ns   4km 6ns
                                       ! dt**0.675  makes larger bin at large distance (<=x2)
                      if(j .eq. 4) dt=dt*10.0 ! for delayed hadrons
                      dt= max(dt, 0.1)
                      call kwhisti( tspec(j, ir, ifai, i),
     *                tmin,  dt, 1000,   b'00000')
                      call kwhistai( tspec(j, ir, ifai,  i),
     *              "Arrival time of "//ptcl(j)//" at (r,faI)",
     *              "rt", "ptcls", .false., 0., "time", "ns")
                  enddo
               enddo
            enddo
         enddo
      endif
      if(tkrtspec) then
         do i = 1, hnsites
            do j = 1, 4
               call kwhistc(tspec0(j,i))
               do ir = 1, nrbin
                  do ifai = 1, nfai
                     call kwhistc(tspec(j, ir, ifai, i))
                  enddo
               enddo
            enddo
         enddo
      endif

      if(tkarspec) then
         do i = 1, hnsites
            do j = 1, 4 
               do ifai = 1, nfai
                  call kwhistc(rspec(j, ifai, i))
               enddo
            enddo
         enddo
      endif
      return
!    ******************
      entry xclearnrfai
!
      
      do i = 1, noofsites
         do j= 1, 4
            do k = 1, nfai
               do l = 1,  nrbin
                  nrfaiall(l, k, j, i) = 0.
                  derfai(l, k, i) = 0.   !  no j (for all charged ptcls)  > 500 keV
                  if(seelowde) derfai2(l, k, i) = 0.   !  no j (for all charged ptcls) < 500 keV
               enddo
            enddo
         enddo
      enddo

      do i = 1, ansites
         do j= 1, 4
            do k = 1, nfai
               do l = 1,  nrbin
                  nrfairec(l, k, j, i) = 0.
               enddo
            enddo
         enddo
      enddo


      return
!     *********************************** hook for Beginning of  1 event
!     *  All system-level initialization for 1 event generation has been
!     *  eneded at this moment.
!     *  After this is executed, event generation starts.
!     *
      entry xbgevent
#if defined (DOMPI)
      etime1 = mpi_wtime()
#endif      

      call cqincident(inci, angle)


#if  FNODATDEF > 0
      bufc=0 
#endif

      cosz = -angle.r(3)
      fai0 = getfai( -angle.r(1), -angle.r(2)) ! in deg.
!          to rotate coordinate and direction cos in the
!          Obsplane =1 system (x is directed to magnetic east
!          y is to the magnetic  north) to the one in the 
!          system with the x-axis directed to the incident
!          direction. All outuput from this program is
!          measured in this system.  (web data too).
!          Z* = Z exp(-(Fai0)) ;  Z in the detector system
!          Z* in the web system.
!
      cosrot = cos(fai0*torad)
      sinrot = sin(fai0*torad)
!       
!
      e0 = inci.p.fm.p(4)
      if(inci.p.code .eq. 9) then
         nn= inci.p.subcode
      elseif(inci.p.code .eq. 1) then
         nn=0
      else
         nn=1
      endif
!         next is only available for parallel job.  For normal job,
!         fisrt col.depth is not yet fixed.
      firstz= zfirst.pos.depth*0.1
      write(0,'(a,1pE11.3,a,E11.3,a,E11.3,a)')
     *     ' E0=',e0, ' cosz=',cosz, ' firstz=',
     *    firstz, ' g/cm2' 
!          here we use enhanced limit
      call howmuch(limite, e0, nn, cosz)

!      
      do i = 1, noofsites
         sumeloss(i) = 0.
         ng(i) = 0.
         ne(i) = 0.
         nmu(i) = 0.
         nhad(i) = 0. 
      enddo

      obstimes = 0.

      return
!     ***************
      entry xobs(atrack, id)
!
!     For id =2, you need not output the z value, because it is always
!     0 (within the computational accuracy).
!
!    ------------------------------------------
      if(fai0 .ne. 0. ) then
!         rotate to the web system
         call det2web( atrack.vec.w.r(1), atrack.vec.w.r(2),
     *                 atrack.vec.w.r(1), atrack.vec.w.r(2))
         call det2web( atrack.pos.xyz.x, atrack.pos.xyz.y,
     *       atrack.pos.xyz.x, atrack.pos.xyz.y)
      endif
!    ------------------------------------------
      obstimes = obstimes + 1.d0
      if(mod(obstimes, 500000.d0) .eq. 0. ) then 
         dosort=.false.
         do i = 1, min(4,stack_pos)
            if(stack(i).p.fm.p(4) .ne. savederg(i)) then
               savederg(i)=stack(i).p.fm.p(4) 
               dosort=.true.
            endif
         enddo
         if(dosort) then
            call csortstack
         endif
#if defined (DOMPI)
         write(0, *)'rank=',mpirank,'  obstimes=', obstimes,
     *    ' ptclE=',atrack.p.fm.p(4)
#else
         write(0, *),'  obstimes=', obstimes,
     *    ' ptclE=',atrack.p.fm.p(4)
#endif
         do i = 1, min(4,stack_pos)
            write(0,*)' stack tops=', stack(i).p.fm.p(4)
         enddo
      endif
!     ***************
      code = atrack.p.code
      ldep =  atrack.where

      if(id .eq. 2 .and. code .le. 6) then
         codex=min(code, 4)
         wz = atrack.vec.w.r(3)  ! downgoing < 0
         if(wz .gt. 0) return
         wz = -wz
         if( keepweight ) then
            wgt = atrack.wgt
            enhance= enhance0*wgt
         else
            wgt = 1.
            enhance= enhance0
         endif

         if(code .eq. kphoton) then
            ng(ldep) = ng(ldep) + enhance
         elseif(code .eq. kelec) then
            ne(ldep) = ne(ldep) + enhance
         elseif(code .eq. kmuon) then
            nmu(ldep) = nmu(ldep) + enhance
         else
            nhad(ldep)=nhad(ldep) + enhance
         endif

         r = sqrt( atrack.pos.xyz.x**2 +
     *                 atrack.pos.xyz.y**2 )

         ek = atrack.p.fm.p(4) -atrack.p.mass

!       ------------- compute energy loss rate
         if(atrack.p.charge .ne. 0  ) then

!             -----------------
!                      /|    |
!                     / |   1g/cm2k
!                    /A |    |
!            -------------------
!                  / ptcl direction  
!         get energy loss when aTrack goes some distance
!         of which vertical gramage is 1g/cm2.
!         Gramage the particle travel is 
!         1/cos where cos is the cos of angle (i.e, A if Fig)
!          in the detctor system.
!         1g/cm^2 = 10-3kg/10-4 m^2 =10 kg/m^2.
!         To travel  1 g/cm^2  along shower axis, the ptcl must
!         run dist kg/m^2
            if(abs(wz) .gt. 1.d-2) then
               dist =10./wz    ! in kg/m2/(g/cm2)
            else
!                 for safety
               dist =1000.
            endif
            if( heobs ) then
!                   the ptcls is the one obsrved at skeleton making time
!                   we must compute dE/dt here
               rho = cvh2den(atrack.pos.height)
               call cdedxinair(atrack.p, rho, dedt) ! dedt; GeV/(kg/m2)
               if(atrack.p.code .eq. kmuon ) then
!                 dE/dx due to muon pair, brem, nuc.i
                  call cmudedx(muni, mubr, mupr, atrack.p.fm.p(4),
     *                 dedxmu)  ! dedxmu in GeV/(g/cm2)
                  dedxmu = dedxmu /10. !  GeV/(kg/m2)
                  dedt = dedt + dedxmu
               endif
            else
!               we can use already computed one
               call cqelossrate(dedt) !  loss rate GeV/(kg/m^2)
            endif
!                       energy in 1 g/cm2 of vertical direction
            eloss =min( dedt*dist, dble(ek)) * enhance  !  GeV/(g/cm2)

            sumeloss(ldep)=sumeloss(ldep) + eloss
         else
            eloss=0.
         endif
!          
!          get web sector indexes (r and fai index )
!
         molu = obssites(ldep).mu
         rinmu =r/molu
         sr = rinmu
         if(rinmu .gt. rminm) then
            ridx= log10(rinmu/rminm)/bin + 1
         else
            ridx =0
         endif
!              fai is  in    -15 to 345  (for dfai=30.)
         fai=getfai(atrack.pos.xyz.x, atrack.pos.xyz.y)
         fai= mod(fai + 360.d0,   360.d0)
         if(fai .gt. (360.d0-dfai/2.0d0)) fai= fai-360.d0
         faiidx=(fai+dfai/2.0d0) /dfai + 1
!
         if(ridx .gt. 0 .and. ridx .le. nrbin ) then
!               for all particles
            nrfaiall(ridx, faiidx, codex, ldep) =
     *           nrfaiall(ridx, faiidx, codex, ldep) + enhance
            if(seelowde) then 
!                 separate low energy contribution ; Eloss = 0 for neutral 
               if(ek .gt. 500.d-6) then
                  derfai(ridx, faiidx, ldep) =
     *                 derfai(ridx, faiidx,  ldep) +  eloss
               else
                  derfai2(ridx, faiidx, ldep) =
     *                 derfai2(ridx, faiidx,  ldep) +  eloss
               endif
            else
!                don't separate low E and high E ptcls
               derfai(ridx, faiidx, ldep) =
     *              derfai(ridx, faiidx,  ldep) +  eloss
            endif  
!
!             ==============
!               for  individual particle
            ji = w2il(ldep)     !  indivdep index.
            if(ji .gt.  0 ) then
               prob =  recprob(ridx, codex, ji)
               if(prob .gt. 1.) then
                  wwgt = wgt
                  accept = .true.
               else
!                 Example:  if prob=10^-3 and wgt=2x10^3, 
!                  record it with weight =2 (=prob*wgt)
                  prob = prob*wgt ! wgt =1 or wgt > 1
                  if(prob .gt. 1.) then
                     wwgt = prob
                     accept = .true.
                  else
                     call rndc(u)
                     if(u .lt. prob) then
                        accept =.true.
                        wwgt=1.
                     else
                        accept = .false.
                     endif
                  endif                  
               endif  ! computtion of accept end
               if(accept) then
                  nrfairec(ridx, faiidx, codex, ji) =
     *              nrfairec(ridx, faiidx, codex, ji) + wwgt
#if  FNODATDEF > 0
                 if(bufc .lt. bufsize) then
                    bufc = bufc + 1
                    buf(bufc).ldep=ldep
                    buf(bufc).code=code
                    buf(bufc).subcode=atrack.p.subcode
                    buf(bufc).charge = atrack.p.charge
                    buf(bufc).ridx=ridx
                    buf(bufc).faiidx= faiidx
                    buf(bufc).rinmu = rinmu
                    buf(bufc).fai= fai
                    buf(bufc).ek = ek
                    buf(bufc).t = atrack.t
                    buf(bufc).wx=-atrack.vec.w.r(1)
                    buf(bufc).wy=-atrack.vec.w.r(2)
                    buf(bufc).wz=wz
                    buf(bufc).wgt = wwgt !  not wgt
                 else
#if BUFSIZE > 200000
                    write(0,*) ' too many partciles --> buf'
                    write(0,*) ' you must make BUFSIZE in interface.f'
                    write(0,*) ' to the standard value (<2x10^5) '
                    stop
#else
                    write(fnodat) bufc, buf
                    bufc= 0
#endif
                 endif  
#else
                 write(*,
     *  '(6i3, 1pE11.3, 0p,f6.1,1p2E11.3,0p,2f8.4,f10.6,1pE11.3)'
     *    )
     *          ldep,  code,  atrack.p.subcode,
     *          atrack.p.charge, ridx, faiidx,
     *          rinmu, fai,
     *          ek, atrack.t, 
     *          -atrack.vec.w.r(1), 
     *          -atrack.vec.w.r(2),  wz, wwgt  ! not wgt
#endif
               endif  ! end accepted ptcl treatment
            endif   ! end individual ptcl output treatment
         endif  ! ridx  >0 ended



!        ================ for histograming
         i=w2hl(ldep)
         if(i .gt. 0 ) then
            if(tkarspec) then
               call kwhist( rspec(codex,faiidx, i),  sr, enhance)
            endif

            if( tkrtspec ) then
               if( ridx .eq. 0 ) then
                  call kwhist( tspec0(codex, i),
     *              sngl( atrack.t ), enhance)
               elseif(ridx .le. nrbin) then
                  call kwhist( tspec(codex, ridx, faiidx,  i),
     *                 sngl( atrack.t ), enhance)
               endif
            endif
         endif       
      endif   ! code and id judge
      return
!     **************
      entry xenevent
!     **************
#if FNODATDEF > 0
#if BUFSIZE < 200000 
      if(bufc .gt. 0) then
         write(fnodat)  bufc, buf
         bufc=0
      endif
      close(fnodat)
#endif
#endif


      firstz= zfirst.pos.depth*0.1

#if defined (DOMPI)
      write(0,*) ' rank=',mpirank, ' closed main file'
      etime2 = mpi_wtime()
      write(0,*) ' elapsed time for this event = ',
     * etime2-etime1
      call mpi_barrier(mpi_comm_world, icon)
#include "inc_gatherNrfai.f"
#include "inc_gatherHyb.f"
      if( mpirank .eq. 0) then
         write(0,*) " rank 0 is reading all data and combining now"
#include "inc_readAndput.f"
      endif  
#endif
#if defined (DOMPI)
      if( mpirank .eq.  0) then
#endif
#include "inc_writeHyb.f"
#include "inc_writeNrfai.f"
#if defined (DOMPI)
      endif
#endif

!       fwollowings are not supported in MPI job.
#if ! defined (DOMPI)
      do i = 1, hnsites
         j = histdep(i)
         write( evid(i), 
     *   '(i3, f7.1,  f6.3, f6.3, i5,  i4)')
     *   histdep(i), asdepthlist(j)*0.1,
     *   asobssites(j).age, asdepthlist(j)*0.1/cog,
     *   int(asobssites(j).mu), int(cog)
      enddo

      if(tkarspec) then
         do i = 1, hnsites
            k=histdep(i)
            do j = 1, 4
               write(dirstr,'(a,"/d",i4, "/")') 
     *              ptcl(j), int( depthlist(k)*0.1 )  
               call kseblk(dirstr, "|", nstr)
               do l = 1, nfai
                  call kwhistdir(rspec(j, l,  i),  dirstr)
                  call kwhists(  rspec(j, l, i), 0. )
                  call kwhistev( rspec(j, l, i), eventno)
                  call kwhistid( rspec(j, l,  i), evid(i))
                  call kwhistp( rspec(j, l, i),  fno)
!                        *********** deallocate ********           
                  call kwhistd( rspec(j, l, i) )
               enddo  ! code loop
            enddo ! fai loop
         enddo ! depth loop
      endif

      if( tkrtspec ) then
         do i = 1, hnsites
            do j = 1, 4   
               call kwhists( tspec0(j,i), 0.)
               call kwhistev( tspec0(j,i), eventno)
               call kwhistid( tspec0(j,i), evid(i))
               k=histdep(i)
               dirstr = " "
               write( dirstr,'(a,"/d",i4, "/")')
     *              ptcl(j), int( asdepthlist(k)*0.1 )
               call kseblk( dirstr, "|", nstr)
               call kwhistdir( tspec0(j,i),  dirstr )
               call kwhistp( tspec0(j,i),  fno )
!                 *********** deallocate ********                  
               call kwhistd( tspec0(j,i) )

               do ifai= 1, nfai
                  do ir= 1, nrbin
                     call kwhists( tspec(j,ir, ifai,i), 0.)
                     call kwhistev(tspec(j,ir, ifai,i), eventno)
                     call kwhistid( tspec(j,ir, ifai,i), evid(i))
                     dirstr = " "
                     write(dirstr,'(a,"/d",i4, "/F",i2,"/")')
     *                    ptcl(j), int( depthlist(k)*0.1), ifai
                     call kseblk(dirstr, "|", nstr)
                     call kwhistdir(tspec(j,ir, ifai,i),  dirstr)
                     call kwhistp( tspec(j,ir, ifai,i),  fno)
!                        *********** deallocate ********           
                     call kwhistd( tspec(j, ir, ifai, i) )
                  enddo
               enddo
            enddo  ! code loop
         enddo   ! depth loop
      endif
#endif

Here is the call graph for this function: