writeHyb.f

Go to the documentation of this file.

      if(observeas) then
#if defined (DOMPI)
         msg = dir(1:lengdir)//"/"//execid(1:lengid)//".hyb"
#define NgXXX NgA
#define NeXXX NeA
#define NmuXXX NmuA
#define NhadXXX NhadA
#define SumElossXXX SumElossA
#else         
         msg = dir(1:lengdir)//"/"//execid(1:lengid)//
     *        "-@."//numb(1:lengn)//".hyb"
#define NgXXX Ng
#define NeXXX Ne
#define NmuXXX Nmu
#define NhadXXX Nhad
#define SumElossXXX SumEloss
#endif
         call copenfw2(fnob, msg, 1, icon)
         if(icon .gt. 1) then
            write(0,*) ' icon=', icon
            call cerrormsg(msg, 1)
            call cerrormsg('could not be opened', 0)
         endif
         cog = 0.
         sumne = 0.

         do i = 1, noofassites

            asobssites(i).esize = asobssites(i).esize* enhance

            if(i .gt. 1 .and. i  .lt. noofassites ) then
               dd =(asdepthlist(i+1) - asdepthlist(i-1))/2.0
            elseif(i .eq. 1) then
               dd =(asdepthlist(2) - asdepthlist(1))
            else
               dd =(asdepthlist(noofassites) -
     *              asdepthlist(noofassites-1))
            endif
            cog = cog + asobssites(i).esize*dd*asdepthlist(i)
            sumne= sumne +asobssites(i).esize*dd
         enddo
!          0.1 is for g/cm2
         cog = cog*0.1/sumne

         cog2 = 0.
         sumne = 0.
         do i = 1, noofassites
            if( asobssites(i).age .gt.
     *          (2.0-asobssites(noofassites).age))  then
               if(i .gt. 1 .and. i  .lt. noofassites ) then
                  dd =( asdepthlist(i+1) - asdepthlist(i-1))/2.0
               elseif(i .eq. 1) then
                  dd =(asdepthlist(2) - asdepthlist(1))
               else
                  dd =(asdepthlist(noofassites) -
     *              asdepthlist(noofassites-1))
               endif
               dd = dd
               cog2 = cog2 + asobssites(i).esize*asdepthlist(i)*dd
               sumne= sumne +asobssites(i).esize*dd
            endif
         enddo
         if(sumne .gt. 0.) then
            cog2 = cog2*0.1/sumne
         else
!              too deep penetration
            cog2 = asdepthlist(noofassites)*0.1
         endif


         if(fnob .ge. 0 )  then
            write(fnob,
     *      '("h ", i4,  3i3, 1pE11.3, 0p 3f11.7, f7.2, 2f7.0)')
     *      eventno,  inci.p.code,
     *      inci.p.subcode, inci.p.charge,
     *      inci.p.fm.e, -angle.r(1), -angle.r(2), -angle.r(3),
     *      firstz, cog, cog2
         else
            write(*,
     *      '("h ", i4,  3i3, 1pE11.3, 0p 3f11.7, f7.2, 2f7.0)')
     *      eventno,  inci.p.code,
     *      inci.p.subcode, inci.p.charge,
     *      inci.p.fm.e, -angle.r(1), -angle.r(2), -angle.r(3),
     *      firstz, cog, cog2
         endif

         do i = 1, noofassites 
            if(fnob .ge. 0) then
               write(fnob, '("t ", i3, 2f7.1,  2f6.3, 1p6E11.3)')
     *          i, 
     *          asdepthlist(i)*0.1,  asobssites(i).mu,
     *          asobssites(i).age,   asdepthlist(i)*0.1/cog2, 
     *          ngxxx(i), nexxx(i),  nmuxxx(i), nhadxxx(i),
     *          asobssites(i).esize, sumelossxxx(i)  
            else
               write(*, '("t ", i3, 2f7.1,  2f6.3, 1p6E11.3)')
     *          i, 
     *          asdepthlist(i)*0.1,  asobssites(i).mu,
     *          asobssites(i).age,   asdepthlist(i)*0.1/cog2, 
     *          ngxxx(i), nexxx(i), nmuxxx(i), nhadxxx(i),
     *          asobssites(i).esize, sumelossxxx(i)  
            endif
         enddo
         if(fnob .gt. 0 ) then
            write(fnob,*)
            close(fnob)
         else
            write(*,*)
         endif
      endif