      program origens_3g_lib_prod
c
c program to update origen-s library based on maxs file
c
      include 'array.inc'
      call flux_read
      call maxs_cs_read
      call maxs_isomer_read 
      call process_cs
      call process_isom
      call update_origen_lib
      stop
      end
c
c
      subroutine flux_read
c
c flux data read
c
      include 'array.inc'
      open (50,file='flux',status='old',err=990)
      read(50,*) (e,flux(i),i=1,igroup)
      close (50)
      go to 999
  990 write(6,*) 'flux file open error!'
      stop
  999 return
      end
c
      subroutine maxs_cs_read
c
c maxs cross section data read
c
      include 'array.inc'
      parameter (maxmt=100)
      dimension cs(maxmt,igroup+1)
      character*12 c12
      character*3 cmt(maxmt)
      character*7 nuclid 
c
      do i=1,ireact
        imt(i)=0
        do j=1,igroup
          csmt(i,j)=0.0
        enddo
      enddo
      open (50,file='maxs_xs',status='old',err=990)
   10 read(50,1000) c12
 1000 format(1x,a12)
      if (c12.ne.' Nuclide ID ') go to 10
      read(50,1010) nuclid 
 1010 format(a6)
      read(nuclid,*) idnucl
   20 read(50,1000) c12
      if (c12.ne.'Number of Re') go to 20
      read(50,1100) nomt
 1100 format(i2)
      read(50,1000) c12
      read(50,1000) c12
      read(50,1200) (cmt(i),i=1,nomt)
 1200 format(20x,50(a3,11x))
      read(50,1000) c12
      do j=1,igroup+1
        read(50,1300) energy(j),(cs(i,j),i=1,nomt)
 1300   format(5x,51(1x,e13.6))
      enddo
      close (50)
      do i=1,nomt
        if (cmt(i).eq.'16 ') then
          imt(1)=1
          do j=1,igroup
            csmt(1,j)=cs(i,j)
          enddo
        endif
        if (cmt(i).eq.'17 ') then    
          imt(2)=1
          do j=1,igroup
            csmt(2,j)=cs(i,j)
          enddo
        endif
        if (cmt(i).eq.'18 ') then    
          imt(3)=1
          do j=1,igroup
            csmt(3,j)=cs(i,j)
          enddo
        endif
        if (cmt(i).eq.'102') then    
          imt(4)=1
          do j=1,igroup
            csmt(4,j)=cs(i,j)
          enddo
        endif
        if (cmt(i).eq.'103') then    
          imt(5)=1
          do j=1,igroup
            csmt(5,j)=cs(i,j)
          enddo
        endif
        if (cmt(i).eq.'107') then    
          imt(6)=1
          do j=1,igroup
            csmt(6,j)=cs(i,j)
          enddo
        endif
      enddo
      go to 999
  990 write(6,*) 'MAXS X-sec file open error!'
      stop
  999 return
      end
c
c
      subroutine maxs_isomer_read
c
c maxs isomer ratio data read
c
      include 'array.inc'
      character*24 c24
      character*5 c5
      character*3 cmt
c
      do i=1,ireact
        isom(i)=0
      enddo
      open (50,file='maxs_isom',status='old',err=990)
   10 read(50,1000) c24
 1000 format(1x,a24)
      if (c24.ne.'Number of Reaction Types') go to 10
      read(50,1100) nomt
 1100 format(i2)
      do i=1,nomt
   20   read(50,1200) c5
 1200   format(1x,a5)
        if (c5.ne.'MT-ID') go to 20
        read(50,1300) cmt
 1300   format(a3)
        read(50,1200) c5
        read(50,1400) noisom
 1400   format(i1)
        read(50,1200) c5
        if (noisom.gt.1) then
          if (cmt.eq.'16 ') then
            isom(1)=noisom 
            do j=1,igroup
              read(50,1500) (risom(1,k,j),k=1,noisom)
 1500         format(4x,10(1x,f8.5))
            enddo
          endif
          if (cmt.eq.'17 ') then
            isom(2)=noisom 
            do j=1,igroup
              read(50,1500) (risom(2,k,j),k=1,noisom)
            enddo
          endif
          if (cmt.eq.'18 ') then
            isom(3)=noisom 
            do j=1,igroup
              read(50,1500) (risom(3,k,j),k=1,noisom)
            enddo
          endif
          if (cmt.eq.'102') then
            isom(4)=noisom 
            do j=1,igroup
              read(50,1500) (risom(4,k,j),k=1,noisom)
            enddo
          endif
          if (cmt.eq.'103') then
            isom(5)=noisom 
            do j=1,igroup
              read(50,1500) (risom(5,k,j),k=1,noisom)
            enddo
          endif
          if (cmt.eq.'107') then
            isom(6)=noisom 
            do j=1,igroup
              read(50,1500) (risom(6,k,j),k=1,noisom)
            enddo
          endif
        endif
      enddo
      close (50)
      go to 999
  990 write(6,*) 'MAXS isomer ratio file open error!'
      stop
  999 return
      end
c
      subroutine process_cs
c
c
c
      include 'array.inc'
      dimension sumflux(3),ig(4)
c
      ig(1)=0
      ig(4)=igroup
      do i=1,igroup+1
        if (energy(i).lt.fastlow) then 
          ig(2)=i-2
          goto 10
        endif
      enddo
   10 do i=ig(2)+1,igroup+1
        if (energy(i).lt.thermup) then 
          ig(3)=i-2
          goto 20
        endif
      enddo
   20 do i=ig(3)+1,igroup+1
        if (energy(i).lt.e0) then 
          ie0=i-1
          goto 30
        endif
      enddo
   30 do i=1,3
        sumflux(i)=0.0
        do j=ig(i)+1,ig(i+1)
          sumflux(i)=sumflux(i)+flux(j)
        enddo
      enddo
      therm=0.0
      do i=ig(3)+1,igroup
        roote=sqrt((energy(i)+energy(i+1))*0.5)
        therm=therm+flux(i)/roote
      enddo
      therm=therm*sqrt((energy(ie0)+energy(ie0+1))*0.5)/sumflux(3)
c     write(6,*) 'therm=',therm
      do i=1,ireact
        do j=1,3
          cs3g(i,j)=0.0
          do k=ig(j)+1,ig(j+1)
            cs3g(i,j)=cs3g(i,j)+flux(k)*csmt(i,k)
          enddo
          if (j.eq.3) then
            cs3g(i,j)=cs3g(i,j)/sumflux(j)/therm
          else
            cs3g(i,j)=cs3g(i,j)/sumflux(j)
          endif
        enddo
      enddo
      do i=1,ireact
        cs0(i)=csmt(i,ie0)
      enddo
      return
      end
c
c
      subroutine process_isom
c
c process isomer data
c
      include 'array.inc'
c
      do i=1,ireact
        f1(i)=0.0
        if (isom(i).le.1) go to 10
        sumfsig=0.0
        do j=1,igroup
          sumfsig=sumfsig+flux(j)*csmt(i,j)
        enddo
        do j=1,igroup
          sumrisom=0.0
C         do k=2,isom(i)
          do k=2,2
            sumrisom=sumrisom+risom(i,k,j)
          enddo
          f1(i)=f1(i)+flux(j)*csmt(i,j)*sumrisom
        enddo
        if (sumfsig.ne.0.0) then
          f1(i)=f1(i)/sumfsig
        else
          f1(i)=0.0
        endif
   10 enddo
      return
      end
c
c
      subroutine update_origen_lib
c
c update origen-s library
c
      include 'array.inc'
      character*85 title,c85
      character*5 name
      character*11 mod,modr
c
      modr='modified   '
      open (50,file='origen.lib',status='old',err=990)
      open (60,file='update.lib',status='new',err=990)
      read(50,1000) title
      write(60,1000) title
 1000 format(a85)
   10 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
 1100 format(1x,i6,4x,5i11,1x,a5,2x,a11)
      if (nucl.eq.idnucl) go to 20
      if (nucl.eq.-1) go to 100
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod 
      do i=1,n1d
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 10
   20 call write_xs
   30 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      if (nucl.eq.-1) go to 100
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      do i=1,n1d
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 30
  100 write(60,1100) nucl
c   actinides
      read(50,1000) title
      write(60,1000) title
  110 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      if (nucl.eq.idnucl) go to 120
      if (nucl.eq.-1) go to 200
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      do i=1,n1d
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 110
  120 call write_xs
  130 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      if (nucl.eq.-1) go to 200
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      do i=1,n1d
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 130
  200 write(60,1100) nucl
c   fission products
      read(50,1000) title
      write(60,1000) title
  210 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      if (nucl.eq.idnucl) go to 220
      if (nucl.eq.-1) go to 300
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod 
      do i=1,n1d+nfy/5
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 210
  220 call write_xs
      do i=1,nfy/5
      read(50,1000) c85
      write(60,1000) c85
      enddo
  230 read(50,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      if (nucl.eq.-1) go to 300
      write(60,1100) nucl,n1d,kfis,nfy,iact,no,name,mod
      do i=1,n1d+nfy/5
        read(50,1000) c85
        write(60,1000) c85
      enddo
      go to 230
  300 write(60,1100) nucl
      close (50)
      close (60)
      go to 999
  990 write(6,*) 'origen or update library open error!'
      stop
  999 return
      end
c
      subroutine write_xs
c
c write x-sec data
c
      include 'array.inc'
      character*5 name
      character*11 mod,modr
      dimension rmt(ireact),cs3g0(3)
c
      modr='modified   '
      rmt(1)=16.0
      rmt(2)=17.0
      rmt(3)=18.0
      rmt(4)=102.0
      rmt(5)=103.0
      rmt(6)=107.0
      backspace (50)
      read(50,1000) nucl,n1d,kfis,nfy,iact,no,name,mod
 1000 format(1x,i6,4x,5i11,1x,a5,2x,a11)
      write(60,1000) nucl,n1d,kfis,nfy,iact,no,name,modr
      do i=1,n1d
        read(50,1100) nucl,rmtno,(cs3g0(k),k=3,1,-1),f10,
     1    name,mod
 1100   format(1x,i6,4x,5e11.4,1x,a5,2x,a11)
        do j=1,ireact
          if ((imt(j).eq.1).and.(rmtno.eq.rmt(j))) then
            if (f10.ne.0.0) then
              write(60,1200) nucl,rmt(j),(cs3g(j,k),k=3,1,-1),f1(j),
     1          name,modr
            else
              write(60,1200) nucl,rmt(j),(cs3g(j,k),k=3,1,-1),f10,
     1          name,modr
            endif
 1200       format(1x,i6,4x,1p5e11.4,1x,a5,2x,a11)
          endif
          if ((imt(j).eq.0).and.(rmtno.eq.rmt(j))) then
            write(60,1200) nucl,rmtno,(cs3g0(k),k=3,1,-1),f10,
     1        name,mod
          endif
        enddo
      enddo
      return
      end
