*-----------------------------------------------------------------------------*
*     FGR13PAK: Dose/Risk Coefficient File Package                            *
*                                                                             *
*     K.F. Eckerman and R.W. Leggett                                          *
*     Oak Ridge National Laboratory                                           *
*                                                                             *
*     Purpose: Provide a FORTRAN module which accesses the dose and risk      *
*              coefficients assembled during the preparation of FGR-13.       *
*                                                                             *
*  FRG13PAK is a collection of FORTRAN routines to assemble the dose and risk *
*  coefficients for the specified radionuclide and its chain members.  The    *
*  package is invoked by the user's code via calls to:                        *
*          OPENEM  <--  to open the data files,                               *
*  calls to the following routines for the specified radionuclide             *
*          NUKEOK  <--  to check the the nuclide name is valid.               *
*          DOSECOF <--  to assemble the decay chain                           *
*  and to CLOSEM to close all open files                                      *
*          CLOSEM  <--  to close the data files.                              *
*  The user's code should type the variable NUKE as character*7 and the       *
*  variables OK and IPATH(3) as logical variables.  IPATH indicates to        *
*  FGR13PAK which coefficients of interest as:                                *
*     IPATH(1) = .true./.false.: Return/not return inhalation coefficients    *
*     IPATH(2) = .true./.false.: Return/not return ingestion coefficients     *
*     IPATH(3) = .true./.false.: Return/not return external coefficients      *
*  The call to NUKEOK is                                                      *
*     CALL NUKEOK(nuke, OK)                                                   *
*  checks for proper notation of the nuclide name and returns OK as .true. if *
*  the nuclide is present in the database, other wise .false.  If OK is true  *
*  then the coefficients are retreived by calling DOSECOF as                  *
*     CALL GETCOF(nuke, ipath)                                                *
*  the output is through the common blocks in the include file DCFPAK.COM.    *
*                                                                             *
*    common /dfacts/chemfinh(mspec,10),chemfing(mspec,10),namage(mage),       *
*   :       cancer(mcan),organ(morg), dfinh(minh,mspec,mage,morg,mlet),       *
*   :       LETh(mspec,mlet), dfing(ming,mspec,mage,morg,mlet),               *
*   :       LETg(mspec,mlet), f1inh(minh,mspec,mage),                         *
*   :       f1ing(ming,mspec,mage), dfext(mspec, mext, morg),                 *
*   :       rinh(minh,mspec,mcan,2), ring(ming,mspec,mcan,2,2),               *
*   :       rext(mspec,mcan,mext,2), NLET(mspec),Type(minh),                  *
*   :       iflag(mspec, mfact), ntypes(mspec), nfings(mspec),                *
*   :       nfinh(mspec), nfing(mspec), nint, next                            *
*                                                                             *
*       chemfinh   character*25 array of chemical form names for inhlation    *
*       chemfing   character*25 array of chemical form names for ingestion    *
*       nameage    character*9 array of names for the mage ages               *
*       cancer     character*9 array of the names of the mcan cancer sites    *
*       organ      character*9 array of the morg organ names                  *
*       dfinh      inhalation absorbed dose coefficient array by type, chain  *
*                  member, age, organ, and let                                *
*       LETh       character array of LET designation (L or H) for ispec for  *
*                  each LET type for inhalation                               *
*       dfinh      ingestion absorbed dose coefficient array by chemical form,* 
*                  chain member, age, organ, and let                          *
*       LETg       character array of LET designation (L or H) for ispec for  *
*                  each LET type for ingestion                                *
*       f1inh      real array of f_1 values for inhalation by type, chain     *
*                  member, and age                                            *
*       f1ing      real array of f_1 values for ingestion by chemical form,   *
*                  chain member and age                                       *
*       dfext      external absorbed dose coefficient array by chain member,  *
*                  external pathway, and organ                                *
*       rinh       inhalation risk coefficient by type, chain member, cancer  *
*                  site, and type of risk (mortality and morbidity)           *
*       ring       ingestion risk coefficient by chemical form, chain member, *
*                  cancer site, pathway, and type of risk (mortality and      *
*                  morbidity)                                                 *
*       rext       external risk coefficient by chain member, cancer site,    *
*                  exposure pathway, and risk (mortality and morbidity)       *
*       NLET       array of the number of LET types for each chain member     *
*       type       character*1 array of the Type designation for each chain   *
*                  member                                                     *
*       iflag      logical array indicating the presence of dose and risk     *
*                  coefficients for each chain member and exposure route      *
*       ntypes     array containing the number of inhalation types for each   *
*                  chain member                                               *
*       nfings     array containing the number of ingestion chemical forms    *
*                  of each chain member                                       *
*       iflag      logical flags for dose factor by chain member and pathway  *
*       nfinh      array of number of special chemical forms for inhalation   *
*       nfing      array of number of special chemical forms for ingestion    *
*       nint       length of chain for internal factors                       *
*       next       length of chain for external factors                       *
*                                                                             *
*     common /radat/ thalf(mspec), iu(mspec), nucnam(mspec),                  *
*    :               branch(mspec, mspec), lmr(mspec), ibr(mspec,mspec),      *
*    :               nbr(mspec), nspec                                        *
*                                                                             *
*       thalf      half-life of chain members                                 *
*       iu         units of the half-lives                                    *
*       nucnam     names of the chain members                                 *
*       branch     branching fraction                                         *
*       lmr        decay constant in 1/d                                      *
*       ibr        pointer to branching fraction ispec --> ispec+1            *
*       nbr        number of branches for chain members                       *
*       nspec      length of chain                                            *
*                                                                             *
* The file FGR13COF.NDX contains the following information: The first record  *
* gives the record numbers of the first (2) and last data record (839),       *
* format (2i4).  The format of records 2 through 839 is                       *
* (a7,e8.0,a2,a8,i7,i5,i6,i4,3(i4,e11.0),e7.0,2e8.0,3i4,i5,i4,i3,2i6,i4,3i6). *
*                                                                             *
* Variable         Description                  Format                        *
* Nuke             Nuclide Name                  a7                           *
* T                Half-life                     e8.0                         *
* Tu               Units of T1/2                 a2                           *
* Mode             Decay mode                    a8                           *
* mdec             Pointer in ICRP38.RAD         i7                           *
* ndec             No radiation records          i5                           *
* mbet             Beta spectrum in ICRP38.BET   i6                           *
* nbet             No points in spectrum         i4                           *
* id1              1st daughter in index         i4                           *
* Y1               Yield of 1st daughter         e11.0                        *
* id2              2nd daughter in index         i4                           *
* Y2               Yield of 2nd daughter         e11.0                        *
* id3              2nd daughter in index         i4                           *
* Y3               Yield of 2nd daughter         e11.0                        *
* Ea               Alpha Energy MeV/nt           e7.0                         *
* Ee               Electron Energy MeV/nt        e8.0                         *
* Ep               Photon Energy MeV/nt          e8.0                         *
* Nplt10           # photon < 10 keV             i4                           *
* Npgt10           # photon >= 10 keV            i4                           *
* Nbeta            # beta particles              i4                           *
* Nelec            # monoenergetic electrons     i5                           *
* Nalpha           # alpha particles             i4                           *
* Nspont           Flag for spf                  i3                           *
* Ing              Ingestion dose coefficient    i6                           *
* Inh              Inhalation dose coefficient   i6                           *
* Iext             External dose coefficient     i4                           *
* Ingr             Ingestion risk                i6                           *
* Inhr             Inhalation risk               i6                           *
* Iextr            External risk                 i6                           *
*                                                                             *
*                                                                             *
*     The following is FGR13PAK's source code.  The routines are presented    *
*     in alphabetical order by function; the order is:                        *
*        1. computational subroutines,                                        *
*        2. screen routines,                                                  *
*        3. functions routines.                                               *
*     The code was written in FORTRAN 77 using the WATCOM FORTRAN Compiler.   *
*     K.F. Eckerman  Dec 2, 1999.                                             *
*-----------------------------------------------------------------------------*
*                                                                             *
*     1. Computational Routines                                               *
*-----------------------------------------------------------------------------*
*
      block data
*
*-----------------------------------------------------------------------
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      data organ /'Adrenals ', 'Bld Wall ', 'B Surface', 'Brain    ',
     :            'Breasts  ', 'Esophagus', 'St Wall  ', 'SI Wall  ',
     :            'ULI Wall ', 'LLI Wall ', 'Kidneys  ', 'Liver    ',
     :            'Lungs    ', 'Muscle   ', 'Ovaries  ', 'Pancreas ',
     :            'R Marrow ', 'Skin     ', 'Spleen   ', 'Testes   ',
     :            'Thymus   ', 'Thyroid  ', 'Uterus   ', '  E(Sv)  ' /
      data indorg / 1,  2,  3,  4,  5, 27,  6,  7,  8,  9, 10, 11,
     :             31, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 33 /
      data iextorg/ 2, 23, 3, 4, 5, 7, 8, 9, 10, 11, 13, 14, 15, 
     :             25, 16, 17, 1, 18, 19, 20, 21, 22, 24, 27 /
      data cancer/'Esophagus', 'Stomach  ', 'Colon    ', 'Liver    ',
     :            'Lung     ', 'Bone     ', 'Skin     ', 'Breast   ',
     :            'Ovary    ', 'Bladder  ', 'Kidneys  ', 'Thyroid  ',
     :            'Leukemia ', 'Residual ', 'Total    '/
      data namage/'Infant   ', '1 yr old ', '5 yr old ', '10 yr old',
     :            '15 yr old', 'Adult    '/
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine batman(b0, zk, zkt, an1, an2, t, n)
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  batman
*   author:   k. f. eckerman
*   date:     01/14/92
*   purpose:  bateman equations.
*
      include 'pakparm.cmn'
*
*     call variables.
      double precision b0, zk, zkt, an1, an2, zero
      integer n
      dimension b0(mspec), zkt(mspec), zk(mspec)
*     local variables.
      double precision s1, s2, ss1, ss2, prod, expfun, expf1
      integer i, j, k
      parameter (zero=0.0d0)
      include 'iolist.cmn'
      an1 = zero
      an2 = zero
      do i = 1, n
        if (b0(i) .ne. zero) then
           s1 = zero
           s2 = zero
           ss1 = zero
           ss2 = zero
           do j = i, n
              prod = zkt(n) / zk(n) * zk(j) / zkt(i)
              do k = i, n
                 if (k .ne. j) prod = prod * zk(k) / (zkt(k) - zkt(j))
              end do
              if (prod .lt. zero) then
                 s1 = s1 + dabs(prod) * expfun(-zkt(j) * dble(t))
                 ss1 = ss1 + dabs(prod) * expf1(zkt(j), dble(t))
              else
                 s2 = s2 + prod * expfun(-zkt(j) * dble(t))
                 ss2 = ss2 + prod * expf1(zkt(j), dble(t))
              end if
           end do
*
*          only positive values are retained; negatives are zero
*
           if (s2 .gt. s1)   an1 = an1 + b0(i) * (s2 - s1)
           if (ss2 .gt. ss1) an2 = an2 + b0(i) * (ss2 - ss1)
        end if
      end do
      return
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine birch(imem, t, rx1, rx2)
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  birch
*   author:   k. f. eckerman
*   date:     01/14/92
*   purpose:  set ipath matrix
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      integer mpath, max 
      common/calcul/ max(mspec), mpath(mspec, mspec) 
*     call variables.
      integer imem
*     local variables.
      integer i, j, mark, jpath, ipath, nmem, m
      double precision zkt, zk, b, b0, an1, an2, x1, x2, zerod
      dimension b(mspec), b0(mspec), zkt(mspec), zk(mspec), mark(mspec), 
     :          jpath(mspec), ipath(mspec)
      parameter(zero = 0.0, zerod = 0.0d0)
*
*     Trace the pathway backwards from Imem to decide which elements
*     of the Mpath matrix to choose.
*
      rx1 = zero
      rx2 = zero
      x1 = zerod
      x2 = zerod
      do i = 1, nspec
        mark(i) = 1
        b(i) = dble(branch(i, i))
      end do
   31 nmem = 1
      jpath(1) = imem
*
      if (max(imem) .eq. 0) goto 35
   33   imem = mpath(mark(imem), imem)
        nmem = nmem + 1
        jpath(nmem) = imem
        if (max(imem) .gt. 0) goto 33
   35 do i = 1, nmem
       ipath(i) = jpath(nmem - i + 1)
      end do
      imem = ipath(nmem)
      do i = 1, nmem
         b0(i) = b(ipath(i))    
         zkt(i) = lmr(ipath(i))
         if (i .lt. nmem) then 
            zk(i) =  dble(branch(ipath(i), ipath(i + 1))) * zkt(i)
         else
            zk(i) = zkt(i)
         end if
      end do
      call batman(b0, zk, zkt, an1, an2, t, nmem)
      x1 = x1 + an1
      x2 = x2 + an2
   60 do 80 i = 1, nmem
        b(ipath(i)) = zerod
        if (i .gt. 1) then
          if (mark(ipath(i)) .ne. max(ipath(i))) then
            m = ipath(i)
            mark(m) = mark(m) + 1
            do j = 1, m - 1
              mark(j) = 1
              b(j) = dble(branch(j, j))
            end do
            goto 31
          end if
        end if
   80 continue
      imem = ipath(nmem)
      rx1 = sngl(x1)
      rx2 = sngl(x2)
      return
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine chain(nuke)
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  chain
*   author:   k. f. eckerman
*   date:     04/06/89 : 09/25/91 : 08/06/92: 02/28/95 :06/04/96
*   purpose:  assemble decay chain.  the name of the parent nuclide should
*             be passed to the routine as namen(1) in the common block
*             /chains/. upon return the chain members will be contained
*             in namen, and ndau is the length of the chain.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      include 'batch.cmn'
* common block /chaind/.
      character*7 named, nuke
      real fhold
      integer iptb, iparb, ibrch, ipar, ipt, n
      logical eob, pob
      common/chaind/named(mspec), fhold(mspec), iptb(mspec),
     :              iparb(mspec), ipt, ibrch, ipar, eob, pob
      common/energy/ealpha(mspec), ebeta(mspec), egamma(mspec)
      dimension eat(mspec), ebt(mspec), egt(mspec)
*     local variables.
      character*78 text
      character*8 t12
      character*2 t12u
      double precision zln2, timest
      parameter(zln2=0.693147181d0, zero=0.0)
      include 'iolist.cmn'
*
*     initialize chain parameters
*
      ibrch = 0
      ipar = 1
      nspec = 1
      eob = .true.
      pob = .false.
      nucnam(1) = nuke
      do i = 1, mspec
         do j = 1, mspec
            branch(i, j) = zero
         end do
      end do
*
*     assign one unit of activity to the parent, rest are zero
*
      branch(1, 1) = 1.0
      lmr(1) = 0.0d0
*
   20 call frward
      if (nspec .le. 0) then
        i = lentrim(nucnam(1))
        write(*,'(1x,(a),'' is not in data base!'')') nucnam(1)(:i)
        return
      endif
      call recver
      if (.not. eob) goto 20
      nspec = nspec - 1
      if (pob) call order
      do i = 1, nspec
         if (nucnam(i)(:2) .eq. 'Sf') then
            lmr(i) = 0.0D0
         else
            lmr(i) = zln2 / timest(thalf(i), iu(i))
         end if
      end do
      t12 = thalf(1)
      t12u = iu(1)
      text = nucnam(1)(:lentrim(nucnam(1))) // ' Decay Chain:'
      text = text(:lentrim(text)) // ' Half-lives and Branching'
      text = text(:lentrim(text)) // ' Fractions'
      if(.not. dbatch) then
         write(*,*) text(:lentrim(text))
         write(olog, '(/a)') text(:lentrim(text))
         call printm
         if (nspec .gt. 5) call cls
      end if
      call path
      timess = 36525.0
*
      if(.not. dbatch) then
      text = ': Activity, Transformations, & Cumulative Energies '
     :       // '(MeV/Bq) at 100y'
      write(*,*) nucnam(1)(:lentrim(nucnam(1))), text(:lentrim(text))
      write(olog,*)nucnam(1)(:lentrim(nucnam(1))),text(:lentrim(text))
      write(*,'(''     Nuclide  T1/2        A(t)/Ao   intA/Ao(d)   Ealph
     .a    Ebeta    Egamma'')')
      write(olog,'(''     Nuclide  T1/2        A(t)/Ao   intA/Ao(d)   Ea
     .lpha    Ebeta    Egamma'')')
      end if
*
      ea = zero
      eb = zero
      eg = zero
      do 50 ispec = 1, nspec
        if (nucnam(ispec)(:2) .eq. 'Sf') goto 50
        call birch(ispec, timess, rx1, rx2)
        ea = ea + 8.64E+04 * rx2 * ealpha(ispec)
        eb = eb + 8.64E+04 * rx2 * ebeta(ispec)
        eg = eg + 8.64E+04 * rx2 * egamma(ispec)
*
        if (.not. dbatch) then
          write(*,'(i4, 1x, a7, 1x, a8, a2, 1p2e12.5,3e9.2)') ispec,
     :     nucnam(ispec), thalf(ispec), iu(ispec), rx1, rx2, ea, eb, eg
          write(olog,'(i4, 1x, a7, 1x, a8, a2, 1p2e12.5,3e9.2)') ispec,
     :     nucnam(ispec), thalf(ispec), iu(ispec), rx1, rx2, ea, eb, eg
        end if
*
        eat(ispec) = ea
        ebt(ispec) = eb
        egt(ispec) = eg
   50 continue
      if (.not. dbatch) call pauseit
      if (nspec .eq. 1) then
         nint = 1
         next = 1
      else
        if (nucnam(nspec)(:2) .eq. 'Sf') then
           nint = icutoff(eat, ebt, egt, .false., nspec-1)
           next = icutoff(eat, ebt, egt, .true., nspec-1)
        else
           nint = icutoff(eat, ebt, egt, .false., nspec)
           next = icutoff(eat, ebt, egt, .true., nspec)
        end if
      end if  
      do jspec = 2, nspec
         n = 0
         do ispec = 1, nspec
            if (nucnam(ispec)(:2) .ne. 'Sf') then
               if (branch(ispec, jspec) .gt. 1.0E-6) then
                  n = n + 1
                  ibr(n, jspec) = ispec
               end if
            end if
         end do
         nbr(jspec) = n
      end do
*
      return
      end
*-----------------------------------------------------------------------
*
      subroutine chekab( nuke )
*
*-----------------------------------------------------------------------
*     routine:  chekab
*     author:   M. Cristy
*     date:     05/05/93
*     purpose:  if nuclide not found, checks whether "a" & "b" isomers
*               exist
*
      character*7 nuke, nukeab, nukea, nukeb, check
      character*6 thalf1, thalf2
      parameter (nnuke = 9)
      dimension nukeab(nnuke), nukea(nnuke), nukeb(nnuke),
     :          thalf1(nnuke), thalf2(nnuke)
      data nukeab /'Eu-150 ', 'In-110 ', 'Ir-186 ', 'Nb-89  ',
     . 'Np-236 ', 'Re-182 ', 'Sb-120 ', 'Sb-128 ', 'Ta-178 '/
      data nukea /'Eu-150a', 'In-110a', 'Ir-186a', 'Nb-89a ',
     . 'Np-236a', 'Re-182a', 'Sb-120a', 'Sb-128a', 'Ta-178a'/
      data nukeb /'Eu-150b', 'In-110b', 'Ir-186b', 'Nb-89b ',
     . 'Np-236b', 'Re-182b', 'Sb-120b', 'Sb-128b', 'Ta-178b'/
      data thalf1/'12.62h', '69.1m', '15.8h', '66m',  '115E3y',
     . '12.7h', '15.89m', '10.4m', '9.31m'/
      data thalf2/'34.2y',  '4.9h',  '1.75h', '122m', '22.5h',
     . '64.0h', '5.76d',  '9.01h', '2.2h'/
      do i = 1, nnuke
        if (nuke .eq. nukeab(i)) then
          write(*,9110) nuke, nukea(i), thalf1(i),nukeb(i),thalf2(i)
          write(*,'(a\)')' Input nuclide or <Enter> to quit)-> '
          read(*,'(bn, a7)') nuke
          nuke = check( nuke )
          if (lentrim(nuke) .eq. 0) stop
          return
        endif
      end do
      return
 9110 format( 4x,'Nuclide ',a,'has 2 isomers:'/20x,a,' with halflife ',a
     ./16x,'and ',a,' with halflife ',a/4x,
     .'Re-input entire name with appropriate "a" or "b" designation',/)
      end
*-----------------------------------------------------------------------
*
      subroutine chekmn( nuke )
*
*-----------------------------------------------------------------------
*     routine: chekmn
*     author:  M. Cristy
*     date:    05/05/93, revised 8/16/93
*     purpose: if "m" isomer is requested, checks whether "n" isomer
*              also exists.
*
      character*7 nuke, nukem, nuken
      character*6 thalfm, thalfn
      character*1 meta
      dimension nukem(3), nuken(3), thalfm(3), thalfn(3)
      data nukem /'Ir-190m', 'Sb-124m', 'Tb-156m'/
      data nuken /'Ir-190n', 'Sb-124n', 'Tb-156n'/
      data thalfm/'1.2h', '93s',   '24.4h'/
      data thalfn/'3.1h', '20.2m', '5.0h'/
      data nnuke /3/
      do i = 1, nnuke
        if (nuke .eq. nukem(i)) then
          write(*, 9110) nuke, nukem(i), thalfm(i), nuken(i), thalfn(i)
   10     write(*, 9120) nukem(i), nuken(i)
          read(*,'(a1)') meta
          if (meta.eq.' ' .or. meta.eq.'m' .or. meta.eq.'M') then
              return
            elseif (meta.eq.'n' .or. meta.eq.'N') then
              nuke = nuke(1:6) // 'n'
              return
            else
              goto 10
          endif
        endif
      end do
      return
 9110 format(4x,'Nuclide ',a,' has 2 metastable isomers:'/20x,
     :a,' with halflife ',a/16x,'and ', a,' with halflife ',a)
 9120 format(4x,'Input <Enter> to accept ',a,', or input "n" for ',a,
     :': ',/)
      end
*-----------------------------------------------------------------------
*
      subroutine chemdat(ipath, ispec)
      
*-----------------------------------------------------------------------
*   routine:  chemdat
*   author:   k. f. eckerman
*   date:     12/10/99
*   purpose:  obtain chemical form infor for element chem
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      character*28 recrd
      character*25 form(10)
      character*2 chem
      integer ipath
      include 'iolist.cmn'
      rewind(i38)
      chem = nucnam(ispec)(:2)
      loop
         read(i38,'(a28)') recrd        ! header record
         if (recrd(:5) .eq. 'START') quit
      end loop   
      if (ipath .eq. 1) then            ! inhalation data
        n = 0
        loop
          read(i38,'(a28)') recrd
          if (recrd(:3) .eq. 'END' .or. recrd(:3) .eq. 'EOF') quit
          if (recrd(:2) .eq. chem) then
             n = n + 1
             chemfinh(ispec, n) = recrd(4:)
          end if
        end loop
        if(n .gt. 0) nfinh(ispec) = n
      else if (ipath .eq. 2) then       ! get ingestion data
         n = 0
        loop                            ! read past the inhalation data
          read(i38,'(a28)') recrd
          if (recrd(:3) .eq. 'END') quit
        end loop
        loop                         ! now at ingestion info
          read(i38,'(a28)')recrd
          loop
            read(i38,'(a28)') recrd
            if (recrd(:3) .eq. 'END' .or. recrd(:3) .eq. 'EOF') quit
            if (recrd(:2) .eq. chem) then
               n = n + 1
               chemfing(ispec, n) = recrd(4:)
            end if
          end loop
          if(n .gt. 0) nfing(ispec) = n           
          if (recrd(:3) .eq. 'END' .or. recrd(:3) .eq. 'EOF') quit
        end loop
      else
        write(*,'(''*** ERROR CHEMDATA.DAT records for '',a2,
     :            '' and pathway '',i2, '' are not valid.'')')
     :   chem, ipath
      end if
      return
      end
*-----------------------------------------------------------------------
*
      subroutine closem
*
*-----------------------------------------------------------------------
*   routine:  closem
*   author:   k. f. eckerman
*   date:     12/3/99
*   purpose:  close files
*
      include 'iolist.cmn'
      close(idex)
      close(i30)
      close(i31)
      close(i32)
      close(i33)
      close(i34)
      close(i35)
      close(i36)
      close(i37)
      close(i38)
      close(olog)
      return
      end
*-----------------------------------------------------------------------
*
      subroutine getcof(nuke, ipath)
*
*-----------------------------------------------------------------------
*     routine:  getcof
*     author:   k.f. eckerman
*     data:     12/01/99
*     purpose:  assemble the decay chain and extract the dose coefficients
*               for the chain members.
*-----------------------------------------------------------------------
*     input variables
*       nuke       parent of the chain in standard notation; e.g., Cs-137
*       ipath(i)   pathway logical flag of size 9. 1 inhalation, 2 ingestion,
*                  3 submersion, 4 ground surface, 5 infinite, 
*
*     output is through the common blocks in the include file fgr13pak.cmn.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      include 'batch.cmn'
      character*7 nuke, nuclide
      character*1 class, let
      integer ipt, ibinry
      logical ipath(*)
      dimension dx(33), ifi(8)
      include 'iolist.cmn'
      data ifi /i30, i31, i32, i33, i34, i35, i36, i37/
*
      if (.not. dbatch) call cls
*
*     zero out the dose factor array
*
      call zerom
*
*     assemble the decay chain
*
      call chain (nuke)
*
*     first do external coefficients for the next chain members
*
      if (ipath(3)) then
        do ispec = 1, next
          iflag(ispec, 3) = .false.
          nuke = nucnam(ispec)
          ipt = ibinry(nuke)
          if (ipt .ne. 0) then
            read(idex,'(a7,a8,a2,a8,126x,i4,12x,i6)', rec=ipt)
     :        nuke, t, ix, mode, iex, iex1
*            write(*,*)'Nuclide - ',Nuke,' pointer - ',ipt,iex,iex1
            if (iex .ne. 0) then
               iflag(ispec, 3) = .true.
               do ip = 1, mext
                  read(ifi(ip+2),'(a7,27e9.0)', rec=iex+3) nuclide,
     :                (dx(j), j=1, 27)
                  do ir = 1, 24     
                     dfext(ispec,ip, ir)= dx(iextorg(ir))
                  end do   
               end do
               jrec = iex1 + 2      ! 3 head records
               do ip = 1, mext
                  do ican = 1, mcan
                     jrec = jrec + 1
                     read(i35,'(a7,57x,e9.0,36x,e9.0)',rec=jrec) 
     :                   nuclide, (rext(ispec,ican,ip,k), k=1,2) 
                  end do
               end do
            end if
          end if
        end do                                   ! End species loop
      end if
*     call pauseit
*
*     now do the inhalation and ingestion coefficients
*
      do 200 ispec = 1, nint
         nuke = nucnam(ispec)
         ipt = ibinry(nuke)
         read(idex,'(a7,132x,2i6,i4,3i6)', rec = ipt) 
     :        nuke, ing, inh, iext, ingr, inhr, iextr
c
c        inhalation
c
         iflag(ispec, 1) = .false.
         if (ipath(1).and. inh .ne. 0) then
           iflag(ispec, 1) = .true.
           read(i31,'(a7,22x,i2,1x,a1)',rec=inh+2)nuclide,nlets,let
           if (nlets.eq. 1) then
              LETh(ispec,1) = let
           else
              LETh(ispec,1) = 'L'
              LETh(ispec,2) = 'H'
           end if
           NLET(ispec) = nlets
           itype = 0
           irec = inh + 1        ! offset 2 for header
           loop                  !  loop over types
             itype = itype + 1
             do iage = 1, 6
               irec = irec + 1
               read(i31,'(a7,i5,f5.0,1x,a1,2x,e8.0,4x,33e10.0)',
     :            rec=irec) nuclide, kage, amad, class,
     :            f1inh(itype,ispec,iage),(dx(j), j = 1, 33)
               do ir = 1, 24
                  dfinh(itype,ispec,iage,ir,1) = dx(indorg(ir))
               end do
               type(itype) = class
               if (nlets .eq. 2) then
                  irec = irec + 1
                  read(i31,'(33x,33e10.0)',rec=irec)(dx(j), j=1,31)
                  do ir = 1, 23
                     dfinh(itype,ispec,iage,ir,2) = dx(indorg(ir))
                  end do 
               end if
             end do
             read(i31,'(a7)',rec=irec+1)nuclide
             if (nuclide .ne. nuke) quit
           end loop
           ntypes(ispec) = itype
           irec = inhr + 2       ! offset 3 for header records
           do itype = 1, ntypes(ispec)
             do ican = 1, mcan
                irec = irec + 1
                read (i37,'(a7,63x,e9.0,36x,e9.0)', rec=irec) nuclide,
     :               (rinh(itype,ispec,ican,k),k=1,2)
             end do
           end do
           call chemdat(1, ispec)
         end if                               ! End of inhalation reads
c
c        ingestion
c
         iflag(ispec, 2) = .false.
         if (ipath(2) .and. ing .ne. 0) then
           iflag(ispec, 2) = .true.
           read(i30,'(a7,13x,i2,1x,a1)',rec=ing+2)nuclide,nlets,let
           if (nlets.eq. 1) then
              LETg(ispec,1) = let
           else
              LETg(ispec,1) = 'L'
              LETg(ispec,2) = 'H'
           end if
           NLET(ispec) = nlets
           itype = 0
           irec = ing + 1  ! offset 2 for header records
           loop            !  loop over types
             itype = itype + 1
             do iage = 1, 6
               irec = irec + 1 
               read(i30,'(a7,i5,e8.0,4x,33e10.0)',
     :            rec=irec) nuclide, kage, f1ing(itype,ispec,iage),
     :            (dx(j), j = 1, 33)
               do ir = 1, 24
                 dfing(itype,ispec,iage,ir,1) = dx(indorg(ir))
               end do   
               if (nlets .eq. 2) then
                  irec = irec + 1
                  read(i30,'(24x,33e10.0)',rec=irec)(dx(j), j = 1, 31)
                  do ir = 1, 23
                    dfing(itype,ispec,iage,ir,2) = dx(indorg(ir))
                 end do   
               end if
             end do
             read(i30,'(a7)',rec=irec+1)nuclide
             if (nuclide .ne. nuke) quit
           end loop
           nfings(ispec) = itype
           irec = ingr + 2     ! offset 3 for header records
           do itype = 1, nfings(ispec)
             do ip = 1, 2
               do ican = 1, mcan
                irec = irec + 1
                   read (i36,'(a7,62x,e9.0,36x,e9.0)', rec=irec) 
     :                   nuclide, (ring(itype,ispec,ican,ip,k),
     :                   k = 1, 2)
               end do
            end do
           end do  
           call chemdat(2, ispec)
         end if  
  200 continue
      return
      end
*------------------------------------------------------------------------------*
*
      subroutine epson
*
*------------------------------------------------------------------------------*
      double precision one, two, eps
      parameter (one = 1.0d0, two = 2.0d0)
      common/epsons/eps
      eps = one
      loop
        eps = sps/two
        if (eps + one .eq. one) quit
      end loop
      eps = dsqrt(eps)
      return
      end  
*------------------------------------------------------------------------------*
*
      subroutine fileini(fnpath, target, nlen, nulog, nuini, program)
*
*------------------------------------------------------------------------------*
*   routine:  fileini
*   author:   k. f. eckerman
*   date:     01/14/92
*   purpose:  opens the file 'program'.ini to find the file name, including path,
*             of target file and its record length nlen.  See the ini file for
*             additional comments.
*   input variable:
*             target = file name passed to this subprogram (8.3 format)
*             program = name of ini file
*   output variable:
*             fnpath = file name including path
*             nlen = recl for direct access file
*
      character*(*) fnpath, target, program
      character*12 fnini, fnstd, lcase
      integer nulog
      logical test
      parameter (maxfil = 40)
      fnini = program(1:lentrim(program)) // '.ini'
      inquire (file = fnini, exist = test)
      if (.not. test) then
         write(*,'(''*** FATAL ERROR: Unable to find INI file '', a12)')
     :    fnini    
         write(nulog,'(''*** FATAL ERROR: Unable to find INI file.'')')
         close(nuIni)
         stop 1
      end if   
      open(nuIni, file=fnini, status='old')
      target = lcase(target)
      do i = 1, maxfil
        read(nuIni, *, end=15) fnstd, fnpath, nlen
        if (fnstd(:3) .eq. 'EOF') then
           write(nuLog,9110) target(:lentrim(target)), 
     :                       fnini(:lentrim(fnini))
           write( *   ,9110) target(:lentrim(target)), 
     :                       fnini(:lentrim(fnini))
           close(nuIni)
           stop 1
        end if   
        if (lcase(fnstd) .eq. target) then
           close(nuIni)
           return
        end if
      end do
   15 write(*,'(''*** FATAL ERROR in function FileIni: file,'',
     : a,a,'' has not EOF record.'')') fnini(:lentrim(fnini))
      write(nulog,'(''*** FATAL ERROR in function FileIni: file,'',
     : a,a,'' has not EOF record.'')') fnini(:lentrim(fnini))
       close(nuIni)
       stop 1
 9110 format(' **** FATAL ERROR in function FileIni: Unable to find the
     :file ', a, ';'/6x,'check ',a,' for proper assignments ****')
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine frward
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  frward
*   author:   k. f. eckerman
*   date:     01/14/92
*   purpose:  read down a branch of a decay chain.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      character*7 named
      real fhold
      integer iptb, iparb, ibrch, ipar, ipt
      logical eob, pob
      common/chaind/named(mspec), fhold(mspec), iptb(mspec), 
     :              iparb(mspec), ipt, ibrch, ipar, eob, pob
      common/energy/ealpha(mspec), ebeta(mspec), egamma(mspec)
*
*     functions referenced.
*
      integer ibinry
      character*8 t, mode
      character*7 nuke, d1
      character*2 ix
      integer j
      include 'iolist.cmn'
*
*     get parent record.
*
      if (ipar .eq. 1) then
         nuke = nucnam(ipar)
         ipt = ibinry(nuke)
         if (ipt .eq. 0) then
            nspec = 0
            return
         endif
      endif
   10 if (ipt .lt. 999) then
         read(idex,'(a7,a8,a2,a8,22x,3(i4,e11.0),e7.0,2f8.0)',
     :      rec=ipt) nuke, t, ix, mode, id1, f1, id2, f2, id3, f3, 
     :      ea, eb, eg
      else
        id1 = 0
        f1 = 0.0
        id2 = 0
        f2 = 0.0
        id3 = 0
        f3 = 0.0
        nuke = 'Sf'
        ea = 0.
        eb = 0.
        eg = 0.
        t = ' '
        ix = ' '
      end if
*
*     ids = 999 denotes "sf" which is not a daughter product, thus set
*     the ids to zero if "sf".
*
      if (id1 .eq. 999 .and. id2 .ne. 0) then
         if (id3 .eq. 0) then
            fhld = f1
            ihld = id1
            f1 = f2
            id1 = id2
            f2 = fhld
            id2 = ihld
         else
            fhld = f1
            ihld = id1
            f1 = f2
            id1 = id2
            f2 = f3
            id2 = id3
            f3 = fhld
            id3 = ihld
         end if
      end if
      if (id2 .eq. 999 .and. id3 .ne. 0) then
         ihld = id2
         fhld = f2
         id2 = id3
         f2 = f3
         id3 = ihld
         f3 = fhld
      end if
*
*     if processing a branch then check to see if d1 has already
*     been included in nucnam, if so only fix up branch(ipar,past) and
*     terminate chain, i.e., chain has converged.
*
      if (pob) then
         if (id1 .gt. 0 .and. id1 .lt. 999) then
            read(idex, '(a7)', rec = id1) d1
         elseif (id1 .eq. 999) then
            d1 = 'Sf'
         end if   
         do j = 1, nspec - 1
            if (d1 .eq. nucnam(j)) goto 16
         end do
         goto 17
*
*        have already handled this daughter; chain has converged.
*        set end of chain and return.
*
   16    branch(ipar, j) = f1
         nucnam(ipar) = nuke
         thalf(ipar) = t
         ealpha(ipar) = ea
         ebeta(ipar) = eb
         egamma(ipar) = eg
         iu(ipar) = ix
         nspec = nspec + 1
         return
      end if
*
*     need to treat this chain member.
*
   17 nucnam(ipar) = nuke
      thalf(ipar) = t
      ealpha(ipar) = ea
      ebeta(ipar) = eb
      egamma(ipar) = eg
      iu(ipar) = ix
      nspec = nspec + 1
      branch(ipar, nspec) = f1
*
*     no further daughters in chain - set end of chain
*
      if (id1 .ne. 0) then
*       further daughters, treat id1 and check for possible branches.
        ipt = id1
*
        if (id2 .ne. 0 ) then
*
*         set end of branch to false, increment branch counter, store 
*         pointer of parent, and record number of second or third daughter
*         while following current chain.  routine recver will direct 
*         recovery of branches.
*
          eob = .false.
          ibrch = ibrch + 1
          iptb(ibrch) = id2
          fhold(ibrch) = f2
          iparb(ibrch) = ipar
          if (id2 .ne. 999) then
             read(idex, '(a7)', rec = id2) named(ibrch)
          else
             named(ibrch) ='Sf'
          end if
        endif
*       third daughter, branch info held as above.
        if (id3 .ne. 0) then
          eob = .false.
          ibrch = ibrch + 1
          iptb(ibrch) = id3
          fhold(ibrch) = f3
          iparb(ibrch) = ipar
          if (id3 .ne. 999) then
             read(idex, '(a7)', rec = id3) named(ibrch)
          else
             named(ibrch) = 'Sf'
          end if
        endif
        ipar = nspec
        if (ipt .ne. 999) goto 10
      endif
      return
*
      end
*-----------------------------------------------------------------------
*
      subroutine nukeok (nuke, ok)
*
*-----------------------------------------------------------------------
*   routine:  nukeok
*   author:   k. f. eckerman
*   date:     01/14/92
*   purpose:  check for valid format of nuclide name and that it exists
*             in the index file
*
      character*7 nuke, check
      logical ok
*
*     function check ensure proper format of user input
*
      nuke = check( nuke )
      call chekmn( nuke )
*
*     find the nuclide in the index file
*
      ipt = ibinry( nuke )
      if (ipt .eq. 0) then
         call chekab( nuke )
         ipt = ibinry( nuke )
      end if
*
*     if pointer ipt is zero, the nuclide is not in the index file
*     else we have a valid nuclide to process
*
      if (ipt .eq. 0) then
         ok = .false.
      else
         ok = .true.
      end if
      return
      end
*-----------------------------------------------------------------------
*
      subroutine openem
*
*-----------------------------------------------------------------------
*     routine:  openem
*     author:   k.f. eckerman
*     date:     12/05/99
*     purpose:  open index and dose coefficient files.
*
      character*64 fpath
      character*12 target
      character*8 prog
      include 'iolist.cmn'
      open(olog, file = 'df_read.log')
      prog = 'fgr13pak'
*
      target = 'index.ndx'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(idex, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'ingestsf.dfs'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i30, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'inhalesf.dfs'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i31, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'submrsin.dfs'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i32, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target='grsurf00.dfs'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i33, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'grvolinf.dfs'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i34, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'external.rsk'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i35, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'ingest.rsk'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i36, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'inhale.rsk'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i37, file=fpath, access='direct', recl=nlen, form=
     :    'formatted', status='old')
*
      target = 'chemdata.dat'
      call fileini(fpath, target, nlen, olog, 41, prog)
      open(i38, file=fpath, status='old')
      return
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine order
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  order
*   author:   k. f. eckerman
*   date:     04/06/89
*   purpose:  order the chain so daughter index > parents.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      common/energy/ealpha(mspec), ebeta(mspec), egamma(mspec)
      character*8 thold
      character*7 nuke
      character*2 ix
      real rsave, csave
      integer i, j, ip, jp, ipass, move
      dimension rsave(mspec), csave(mspec)
      include 'iolist.cmn'
*
*     move # of elements to move
*
      ipass = 0
      loop                          ! ye old 100 label
      move = 0
      ipass = ipass + 1
      if (ipass .gt. 4*nspec) then
         write(olog,'('' Failure in routine order: greater than'',i3,
     :   '' passes for '', a7, ''.'')') ipass, nucnam(1)
         write(*,'('' Failure in routine order: greater than'',i3,
     :   '' passes for '', a7, ''.'')') ipass, nucnam(1)
         stop 1
      endif
*
      do i = 1, nspec
         do j = 1, i-1
            if (branch(i, j) .ne. 0.) then
               ip = i
               jp = j
               move = 1
               go to 15
            endif
         end do
      end do
*
*     if no elements to move then return
*
  15  if (move .eq. 0) quit     ! ye old return statement
      nuke = nucnam(ip)
      thold = thalf(ip)
      ea = ealpha(ip)
      eb = ebeta(ip)
      eg = egamma(ip)
      ix = iu(ip)
      do j = 1, nspec
         rsave(j) = branch(ip, j)
      end do
      do i = ip - 1, jp, -1
         nucnam(i + 1) = nucnam(i)
         thalf(i + 1) = thalf(i)
         ealpha(i + 1) = ealpha(i)
         ebeta(i + 1) = ebeta(i)
         egamma(i + 1) = egamma(i)
         iu(i + 1) = iu(i)
         do j = 1, nspec
            branch(i + 1, j) = branch(i, j)
         end do
      end do
      nucnam(jp) = nuke
      thalf(jp) = thold
      iu(jp) = ix
      ealpha(jp) = ea
      ebeta(jp) = eb
      egamma(jp) = eg
      do j = 1, nspec
         branch(jp, j) = rsave(j)
      end do
      do i = 1, nspec
         csave(i) = branch(i, ip)
      end do
      do j = ip - 1, jp, -1
         do i = 1, nspec
            branch(i, j + 1) = branch(i, j)
         end do
      end do
      do i = 1, nspec
         branch(i, jp) = csave(i)
      end do
      end loop       ! ye old goto 100 statement
      return
*
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine path
*                                                                      *
*----------------------------------------------------------------------*
*     author:  K.F. Eckerman
*     date:    04/06/89
*     purpose: initialize mpath and max matrices
*              Adopted from A. Birchall, Health Phys. 50, 3, 389-397, 1986.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      integer max, mpath
      common/calcul/ max(mspec), mpath(mspec, mspec)
*     initializes mpath and max matrics.
      do i = 1, nspec
        max(i) = 0
        do j = 1, nspec
           mpath(i,j) = 0
        end do
      end do
      do j = 2, nspec
         do i = 1, j - 1
           if (branch(i, j) .ne. 0.) then
              max(j) = max(j) + 1
              mpath(max(j), j) = i
           end if
         end do
      end do
      return
      end
*-----------------------------------------------------------------------
*
      subroutine pauseit
*
*-----------------------------------------------------------------------
*     routine: pause
*     author:  K.F. Eckerman
*     date:    10/23/93
*     purpose: pause w/o line feed
*
      character*1 dumy
      write(*,'(a\)')' Hit <Enter> to continue.'
      read(*,'(bn, a1)') dumy
      return
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine printm
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  printm
*   author:   k. f. eckerman with modifications by j. c. ryman
*   date:     01/15/92
*   purpose:  print the decay chain
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
*
*     local variables.
      character*7 nuke, rlist
      integer i, j, k, nlist, jlist
      dimension jlist(mspec), nuke(mspec), rlist(mspec)
      include 'iolist.cmn'
      write(*, 8001)
      write(olog, 8001)
      if (nspec .eq. 1) then
         write(*, 8002) nspec, nucnam(1), thalf(1), iu(1)
         write(olog, 8002) nspec, nucnam(1), thalf(1), iu(1)
      else
         do 20 i = 1, nspec
            if(nucnam(i)(:2) .eq. 'Sf') goto 20
            nlist = 0
            do j = 1, nspec
               if (i .ne. j .and. branch(i,j) .ne. 0.0) then
                  nlist = nlist + 1
                  jlist(nlist) = j
                  nuke(nlist) = nucnam(j)
                  write (rlist(nlist), '(1pe7.1)') branch(i,j)
                  do k = 4, 6
                     rlist(nlist)(k:k) = rlist(nlist)(k+1:k+1)
                  end do
                  rlist(nlist)(7:7) = ' '
               endif
            end do
            write(*, 8002) i, nucnam(i), thalf(i), iu(i),
     :                    (rlist(j), jlist(j), nuke(j), j = 1, nlist)
            write(olog, 8002) i, nucnam(i), thalf(i), iu(i),
     :                    (rlist(j), jlist(j), nuke(j), j = 1, nlist)

   20    continue
         if (nspec .gt. 5) pause
      endif
*
*     clear screen and return
*
      if (nspec .gt. 5) call cls
      return
*
*     formats
*
 8001 format(' ',3x,'Nuclide  Halflife    f1',7x,'Nuclide   f2',7x,
     : 'Nuclide   f3',7x,'Nuclide')
 8002 format(' ',i2,1x,a7,1x,a8,a2,:,3(1x,a6,'->',i2,1x,a7))
      end
*----------------------------------------------------------------------*
*                                                                      *
      subroutine recver
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  recver
*   author:   k. f. eckerman
*   date:     04/06/89
*   purpose:  recover info on branches in the chain that were detected
*             by frward and direct the reading of the new branch.
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      character*7 named
      real fhold
      integer iptb, iparb, ibrch, ipar, ipt
      logical eob, pob
      common/chaind/named(mspec), fhold(mspec), iptb(mspec),
     :              iparb(mspec), ipt, ibrch, ipar, eob, pob
* local variables.
      character*7 nuke
      integer i
      include 'iolist.cmn'
*
*     no branches to treat, set end of branch to true and return.
*
    1 if (ibrch .eq. 0) then
         eob = .true.
*      elseif (iptb(ibrch) .eq. 999) then
*         eob = .true.
      else
*
*        consider remaining branches. recover parent's
*        index at branch (ipar) and daughter's record number (ipt).  
*        decrement branch counter and return.
*
         pob = .true.
         ipar = iparb(ibrch)
         ipt = iptb(ibrch)
         nuke = named(ibrch)
*
*        need to check to see of the daughter of the branch has already
*        a member of the chain.
*
         do i = 1, nspec - 1
            if (nuke .eq. nucnam(i)) goto 15
         end do
         nucnam(nspec) = nuke
         branch(ipar, nspec) = fhold(ibrch)
         ibrch = ibrch - 1
         ipar = nspec
         return
*
*        if already a chain member set r, decrement the branch counter 
*        and look for another branch to process.
*
  15     branch(ipar,i) = fhold(ibrch)
*
         ibrch = ibrch - 1
         go to 1
*
      endif
*
      return
      end
*-----------------------------------------------------------------------
*
      subroutine zerom
*
*-----------------------------------------------------------------------
*     routine:  zerom
*     author:   k.f. eckerman
*     date:     12/10/99
*     purpose:  set the dose coefficient arrays to zero
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      do ispec = 1, mspec
         nfinh(ispec) = 0
         nfing(ispec) = 0
         do ifact = 1, mfact
            iflag(ispec, ifact) = .false.
         end do
      end do
      do i = 1, minh
        do j = 1, mspec
          do k = 1, mage
            do l = 1, morg
              do m = 1, mlet
                dfinh(i, j, k, l, m) = 0.0
              end do  
            end do
          end do  
        end do
      end do
      do i = 1, ming
        do j = 1, mspec
          do k = 1, mage
            do l = 1, morg
              do m = 1, mlet
                dfing(i, j, k, l, m) = 0.0
              end do  
            end do
          end do  
        end do
      end do
      do i = 1, mspec
         do j = 1, mext
            do k = 1, morg
               dfext(i, j, k) = 0.0
            end do
         end do
      end do
      do j = 1, mspec
        do k = 1, mcan
          do l = 1, mlet
            do i = 1, minh
              rinh(i, j, k, l) = 0.0
            end do
            do i = 1, ming
              do m = 1, 2
                ring(i, j, k, l, m) = 0.0
              end do  
            end do  
          end do
        end do  
      end do
      do i = 1, mspec
         do j = 1, mcan
            do k = 1, mext
               do l = 1, 2
                  rext(i, j, k, l) = 0.0
               end do
            end do
         end do
      end do
      return
      end

*-----------------------------------------------------------------------------*
*
*     2. Screen Routines
*
*-----------------------------------------------------------------------------*
*                                                                             *
      subroutine cls 
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     12/08/93
*     routine to clear screen 
      write(*,*) char(27),'[2J' 
      return 
      end 
*-----------------------------------------------------------------------------*
*                                                                             *
       subroutine curright(icol)
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     12/08/93
*     move the cursor icol columns right on the display.
*
      character*2  col
      write(col, '(i2.2)') icol
      write(*,'(a\)') ' ' // char(27) // '[' // col //'C'
      return
      end
*-----------------------------------------------------------------------------*
*                                                                             *
       subroutine curpos(irow, icol)
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     12/08/93
*     move cursor to the indicated row and column.
*
      character*2  row, col
      write(row, '(i2.2)') irow
      if (icol .ne. 0) then
         write(col, '(i2.2)') icol
      else
         col ='01'
      end if
      write(*,'(a\)') ' ' // char(27) // '[' // row // ';' // col // 'H'
      return
      end
*-----------------------------------------------------------------------------*
*
*     3. Function Routines
*
*-----------------------------------------------------------------------
*
      character*(*) function check(nuke)
*
*-----------------------------------------------------------------------
*     function: check
*     author:   r.j. westfall
*     date:     07/20/89
*     purpose:  convert chemical symbol in nuclide name to proper
*               notation, e.g.; Kr-85m, etc.
*
      character*(*) nuke
      character*7 ltrim
*
*     remove any leading blanks from nuke
*
      nuke = ltrim( nuke )
*
*     ensure first character is upper case.
*
      if (nuke(:1) .ge. 'a' .and. nuke(:1) .le. 'z')
     :    nuke = char(ichar(nuke(:1)) - 32) // nuke(2:7)
*
*     ensure second character, if present, is lower case.
*
      if (nuke(2:2) .ge. 'A' .and. nuke(2:2) .le. 'Z')
     :    nuke = nuke(:1) // char(ichar(nuke(2:2)) + 32) // nuke(3:7)
*
*     ensure metastable notation, if present, is lower case.
*
      do j = 4, 7
         if (nuke(j:j) .ge. 'A' .and. nuke(j:j) .le. 'Z')
     :   nuke = nuke(:j-1) // char(ichar(nuke(j:j)) + 32) // nuke(j+1:)
      end do
      check = nuke
      end
*-----------------------------------------------------------------------------* 
*                                                                             *
      double precision function expf1 (lm, t) 
*                                                                             *
*-----------------------------------------------------------------------------* 
*     author:   k. f. eckerman
*     date:     12/06/99
*     purpose:  routine to compute [1.0 - exp(-lm * t)] / lm. 
*
      double precision lm, t, lmt, one, two, expfun, eps
      parameter(one = 1.0d0, two = 2.0d0)
      common/epsons/eps
      lmt = lm * t
      if (lmt .lt. eps) then
         expf1 = t * (one - lmt / two)
      else
         expf1 = (one - expfun(-lmt)) / lm
      end if 
      return
      end 
*-----------------------------------------------------------------------------*
*                                                                             *
      double precision function expfun (t) 
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     10/04/94
*     purpose:  routine to compute exp (t). 
*  
      double precision t, zero, upval 
      parameter(zero = 0.0d0, upval = 180.d0)
      if (t .lt. -upval) then
         expfun = zero 
      else
         expfun = dexp(t)
      end if
      return 
      end 
*-----------------------------------------------------------------------
*
      integer function ibinry ( target )
*
*-----------------------------------------------------------------------
*     function: ibinry
*     author:   k.f. eckerman
*     date:     06/20/88
*     purpose:  locate record sort by target key
*
      integer left, right, try
      character*7 target, a1
      include 'iolist.cmn'
*
*     initialization.
*
      read (idex, '(2i4)', rec = 1) left, right
*
*     begin attempts to find target.
*
   10 try = int((left + right) / 2)
      read (idex, '(a7)', rec = try) a1
      if (a1 .lt. target) then
         left = try + 1
      elseif (a1 .gt. target) then
         right = try - 1
      else
         ibinry = try
         return
      end if
*
*     continue search unless left > right then set ibinry to zero and
*     let calling deal with the unidentified target.
*
      if (left .lt. right + 1) then
         goto 10
      else
         ibinry = 0
         return
      end if
      end
*----------------------------------------------------------------------*
*                                                                      *
      integer function icutoff(eat, ebt, egt, eflag, nspec)
*                                                                      *
*----------------------------------------------------------------------*
*   routine:  icutoff
*   author:   k. f. eckerman
*   date:     04/06/93
*   purpose:  cut of the decay chain at member for which at least 99%
*             of the alpha, electron, and photon energy has been 
*             released during a 100-y period.
*
      logical eflag
      dimension eat(*), ebt(*), egt(*)
      ea = eat(nspec)
      eb = ebt(nspec)
      eg = egt(nspec)
      if (ea .gt. 0.0) then
         do i = nspec-1, 1, -1
            if (eat(i) .lt. 0.99 * ea) then
               ia = i + 1
               goto 15
            end if
         end do
         ia = 1
      else
         ia = 0
      end if
   15 if (eb .gt. 0.0) then
         do i = nspec-1, 1, -1
            if (ebt(i) .lt. 0.99 * eb) then
               ib = i + 1
               goto 25
            end if
         end do
         ib = 1
      else
         ib = 0
      end if
   25 if (eg .gt. 0.0) then
         do i = nspec-1, 1, -1
            if (egt(i) .lt. 0.99 * eg) then
               ig = i + 1
               goto 35
            end if
         end do
         ig = 1
      else
         ig = 0
      end if
   35 if (eflag) then
        icutoff = max0(ig, ib)
      else
        icutoff = max0(ig, max0(ia, ib))
      end if  
      return 
      end
*-----------------------------------------------------------------------------*
*                                                                             *
      character*(*) function lcase(a)
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     10/04/94
*     purpose:  convert character variable a to lower case.
*
      character*(*) a
      lcase = a
      do i = 1, lentrim(lcase)
         ix = ichar(lcase(i:i))
         if (ix .gt. 64 .and. ix .lt. 91) then 
            lcase(i:i) = char(ix + 32)
         end if
      end do
      return
      end
*-----------------------------------------------------------------------
*
      character*(*) function ltrim(a)
*
*-----------------------------------------------------------------------
*     function: ltrim
*     author:  K.F. Eckerman
*     date:    10/23/93
*     purpose: trim leading blanks from string a
*
      character*(*) a
      logical ok
      ok = .false.
   10 if (a(1:1) .ne. ' ') ok = .true.
      if (.not. ok) then
         n = lentrim(a) - 1
         do i = 1, n
            a(i:i) = a(i+1:i+1)
         end do
         a(n+1:n+1) = ' '
         goto 10
      end if
      ltrim = a
      return
      end
*-----------------------------------------------------------------------------*
*                                                                             *
      double precision function timest(t, ix)
*                                                                             *
*-----------------------------------------------------------------------------* 
*     author:   k. f. eckerman
*     date:     10/04/94
*     purpose:  function returns time in days given time string t and its
*               units ix.
*     function arguments.
      character*2 ix
      character*8 t
      double precision tp
*
      read(t,'(E8.0)') tp
      if (ix .eq. 'us') then
         tp = tp / 8.64d+10
      elseif (ix .eq. 'ms') then
         tp = tp / 8.64d+07
      elseif (ix .eq. 's ') then
         tp = tp / 8.64d+04
      elseif (ix .eq. 'm ') then
         tp = tp / 1.44d+03
      elseif (ix .eq. 'h ') then
         tp = tp / 24.d0
      elseif (ix .eq. 'y ') then
         tp = tp * 365.25d0
      endif
      timest = tp
      return
      end 
*-----------------------------------------------------------------------------*
*                                                                             *
      character*(*) function ucase(a)
*                                                                             *
*-----------------------------------------------------------------------------*
*     author:   k. f. eckerman
*     date:     10/04/94
*     purpose:  convert character variable a to upper case.
*
      character*(*) a
      ucase = a
      do i = 1, lentrim(ucase)
         ix = ichar(ucase(i:i))
         if (ix .gt. 96 .and. ix .lt. 123) then
            ucase(i:i) = char(ix - 32)
         end if
      end do
      return
      end
c-----------------------------------------------------------------------
c
c      integer function len_trim(a)
c
c-----------------------------------------------------------------------
c     author:   k.f. eckerman
c     date:     06/25/96
c     purpose:  determine the trim length of a character variable a
c               in the manner of Microsoft len_trim function. if the
c               len_trim function is not supported by the FORTRAN
c               compiler then active this routine by removal of the
c               'c's on the statements below and in the above name.
c      character*(*) a
c      n = len(a)
c      do i = n, 1, -1
c         if (a(i:i) .ne. ' ' .and. ichar(a(i:i)) .ne. 0) then
c            len_trim = i
c            return
c         end if
c      end do
c      len_trim = 0
c      end
