add2hyb.f

Go to the documentation of this file.

      implicit none
!   Suppose  we have n hybrid data
!     h1, h2, h3, ....hn
!   1) cp h1 h0
!   2) h0 + h2 --> h; mv  h h0
!   3) h0 + h3 --> h; mv  h h0
!   ..
!   n) h0 + hn --> h; mv  h h0
! 
!   This program add two hybrid data; h0 + hx--> h
!            environmental variable
!   file h0: HYBFILE0
!   file hx: HYBFILEX
!   file h:  HYBFILET
!

      integer ndepth
      parameter(ndepth= 50)
      real*8  ASdep(ndepth),  muunit(ndepth), sumEsize, sumEsize2
      real*8  Esize0(ndepth),
     *     age0(ndepth),  cogdep0(ndepth),
     *     seloss0(ndepth),
     *     ng0(ndepth), ne0(ndepth), nmu0(ndepth), nhad0(ndepth),
     *     cog0, cog20
      real*8  Esizex(ndepth),
     *     agex(ndepth),  cogdepx(ndepth),
     *     selossx(ndepth), 
     *     ngx(ndepth), nex(ndepth), nmux(ndepth), nhadx(ndepth),
     *     cogx, cog2x

      real*8  Esizet(ndepth),
     *     aget(ndepth),  cogdept(ndepth),
     *     selosst(ndepth),
     *     ngt(ndepth), net(ndepth), nmut(ndepth), nhadt(ndepth),
     *     cogt, cog2t
            
      integer EvNo0, EvNox
      integer fn0, fnx, fnt
      integer kgetenv2
      integer klena
      integer code, subcode, charge
      integer leng, i, j, lengflesh
      integer icon0, iconx, icont
      real dd, firstz, w1, w2, w3, E0
      character*128 hyb0, hybx, hybt, flesher
      character*128 input0, inputx, inputt
      fn0 = 2
      fnx = 3
      fnt = 4
      leng = kgetenv2("HYBFILE0", hyb0)
      call copenfw2(fn0, hyb0, 1, icon0)
      if(icon0 .ne. 1)  then
         write(0,*) hyb0(1:leng)
         if( icon0 .eq. 0) then
            write(0,*) 'not exists'
         else
            write(0,*) ' cannot be opened '
         endif
         write(0,*) ' icon=',icon0
         stop 9999
      else
         write(0,*)  hyb0(1:leng), ' opened'
      endif
      leng = kgetenv2("HYBFILEX", hybx)
      call copenfw2(fnx, hybx, 1, iconx)
      if(iconx .ne. 1)  then
         write(0,*) hybx(1:leng)
         if( iconx .eq. 0) then
            write(0,*) 'not exists'
         else
            write(0,*) ' cannot be opened '
         endif
         write(0,*) ' icon=',iconx
         stop 9999
      else
         write(0,*)  hybx(1:leng), ' opened'
      endif

      
      leng = kgetenv2("HYBFILET", hybt)
      call copenfw2(fnt, hybt, 1, icont)
      if(icont .ne. 0)  then
         write(0,*) hybt(1:leng)
         write(0,*) ' cannot be opened '
         write(0,*) ' icon=',icont
         stop 9999
      else
         write(0,*)  hybt(1:leng), ' opened'
      endif
      lengflesh = kgetenv2("FLESHDIR", flesher)

      i =0
      do while(.true.)
         input0 = ' '
         read( fn0, '(a)',  end=1000 ) input0
         inputx = ' '
         read( fnx, '(a)'  ) inputx 
         if(input0 .ne. " ") then
            if(flesher(1:lengflesh) .eq. "FleshHist") then
               if(input0(1:1) .eq. "h" ) then
                  read(input0(3:klena(input0)),* ) evno0, code,
     *            subcode, charge, e0, w1, w2, w3,firstz, 
     *            cog0, cog20

                  read(inputx(3:klena(inputx)),*) evnox, code,
     *            subcode, charge, e0, w1, w2, w3,firstz,
     *            cogx, cog2x

               else
                  read(input0(3:klena(input0)), *) 
     *            i, asdep(i),  muunit(i), age0(i), cogdep0(i),
     *            ng0(i), ne0(i), nmu0(i), nhad0(i), 
     *            esize0(i),  seloss0(i)

                  read(inputx(3:klena(inputx)), *) 
     *            i, asdep(i),  muunit(i), agex(i), cogdepx(i),
     *            ngx(i), nex(i), nmux(i), nhadx(i), 
     *            esizex(i),  selossx(i)
               endif
            else
               read(input0(3:klena(input0)), *) 
     *            i, asdep(i),  muunit(i), age0(i), cogdep0(i),
     *            ng0(i), ne0(i), nmu0(i), nhad0(i), 
     *            esize0(i)
               selossx(i) =0.

               read(inputx(3:klena(inputx)), *) 
     *            i, asdep(i),  muunit(i), agex(i), cogdepx(i),
     *            ngx(i), nex(i), nmux(i), nhadx(i), 
     *            esizex(i)
               selossx(i) =0.
            endif
         else
!              1 event read
            if(inputx .ne. " ") then
               write(0,*) ' event differ', inputx
               stop 1111
            endif
            do j = 1, i
               esizet(j) = esize0(j) + esizex(j)
               if(esizet(j)  .gt. 0.) then
                  aget(j) = ( esize0(j)*age0(j) + esizex(j)*agex(j))/
     *               esizet(j) 
               else
                  aget(j)= 0.
               endif
               selosst(j) = seloss0(j) + selossx(j)
               ngt(j) = ng0(j) + ngx(j)
               net(j) = ne0(j) + nex(j)
               nmut(j) = nmu0(j) + nmux(j)
               nhadt(j) = nhad0(j) + nhadx(j)
            enddo
            cogt =0.
            cog2t =0.
            sumesize = 0. 
            sumesize2 = 0. 
            do j = 1, i
               if(j .gt. 1 .and. j  .lt. i ) then
                  dd =( asdep(j+1) - asdep(j-1))/2.0
               elseif(j .eq. 1) then
                  dd =(asdep(2) - asdep(1))
               else
                  dd =(asdep(i) - asdep(i-1))
               endif
               cogt = cogt + esizet(j)*asdep(j)*dd
               sumesize = sumesize + esizet(j)*dd
               
               if( aget(j) .gt.  2.0-aget(i)  )  then
                  if(j .gt. 1 .and. j  .lt. i ) then
                     dd =( asdep(j+1) - asdep(j-1))/2.0
                  elseif(j .eq. 1) then
                     dd =(asdep(2) - asdep(1))
                  else
                     dd =(asdep(i) - asdep(i-1))
                  endif
                  cog2t = cog2t + esizet(j)*asdep(j)*dd
                  sumesize2= sumesize2 + esizet(j)*dd
               endif
            enddo
            cogt =  cogt / sumesize
            if(sumesize2 .gt. 0.) then
               cog2t =  cog2t / sumesize2
            else
               cog2t = asdep(i)
            endif
            do j = 1, i
               cogdept(j) = asdep(j)/cog2t
            enddo
            write(fnt, 
     *         '("h ", i4,  3i3, 1pE11.3, 0p 3f11.7, f7.2, 2f7.0)')
     *            evno0, code,
     *            subcode, charge, e0, w1, w2, w3,firstz, 
     *            cogt, cog2t
            do j = 1, i
               if(flesher(1:lengflesh) .eq. "FleshHist") then
                  write(fnt, '("t ", i3, 2f7.1,  2f6.3,
     *                 1p6E11.3)')
     *                 j,
     *                 asdep(j), muunit(j),
     *                 aget(j),  cogdept(j),
     *                 ngt(j), net(j), nmut(j), nhadt(j),
     *                 esizet(j), selosst(j)
               else
                  write(fnt, '("t ", i3, 2f7.1,  2f6.3,
     *                 1p5E11.3)')
     *                 j,
     *                 asdep(j), muunit(j),
     *                 aget(j),  cogdept(j),
     *                 ngt(j), net(j), nmut(j), nhadt(j),
     *                 esizet(j)
               endif
            enddo
            write(fnt, *) 
         endif
      enddo
 1000 continue
      write(0,*) 'all events processed '
      end