      program readem
*------------------------------------------------------------------------------*
*                                                                              *
* READEM illustrates how to access the FGR-13 dose and risk coefficients       *
* via FORTRAN.  The index file FGR13COF.NDX provides the pointers into the     *
* FGR-13 coefficient files and assembles the decay chain for the specified     *
* radionuclide. The dose coefficient files are:                                *
*   F12TIII1.DAT - FGR 12 Table III.1 for air submersion                       *
*   F12TIII3.DAT - FGR 12 Table III.3 for ground surface contamination         *
*   F12TIII7.DAT - FGR 12 Table III.7 contaminated soil infinite thickness     *
*   FGR13ING.DAT - Age-specific absorbed dose coefficients for ingestion       *
*   FGR13INH.DAT - Age-specific absorbed dose coefficients for inhalation,     *
* the risk coefficient files are:                                              *
*   FGR13EXT.RBS - Details of risk via external exposures                      *
*   FGR13ING.RBS - Details of risk via ingestion intakes                       *
*   FGR13INH.RBS - Details of risk via inhalation intakes                      *
* and a file of chemical form information is                                   *
*   CHEMDATA.DAT - Names of chemical forms for identification of coefficients. *
*                                                                              *
* All coefficient files and FGR13CD.NDX are formatted direct access files      *
* and can be examined using an ANSI editor, includeding Windows WORDPAD.       *
* The file CHEMDATA.DAT is a sequential dat file and includes notes regarding  *
* its format and usuage.                                                       *
* Note that attempts to edit the direcet access files may make them unuseable. *
* The RECL specification in the FORTRAN open statement for the direct access   *
* files is obtained from the file DCFPAK.INI.  The purpose of the INI file is  *
* to provide the full path to these files. See DCFPAK.INI for information      *
* regarding useage with different FORTRAN compliers.                           *
*                                                                              *
* author:  k.f. eckerman                                                       *
* date:    12/2/99                                                             *
*------------------------------------------------------------------------------*
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      include 'batch.cmn'
      character*7 nuke
      logical ipath(3), ok
      include 'iolist.cmn'
*
*     setting dbatch to true results in the DCFPAK_P package running in a
*     silent mode; i.e., not print to the screen.
*      
      dbatch = .false.
*
      call cls
      write(*,'(1x,'' View FGR-13 Dose & Risk Coefficients'')')
      write(*,'(3x,'' K.F. Eckerman and R.W. Leggett'')')
      write(*,'(3x,'' Oak Ridge National Laboratory'')')
      write(*,'(3x,'' Oak Ridge, TN 37831-6480'',/)')
*
*     call epson to compute some machine constants
*
      call epson      
*
*     following do-loop sets pathway flags to true; i.e., DCFPAK_P will return 
*     cofficients for all pathways.
*
      do i = 1, 3
         ipath(i) = .true.
      end do
*
*     open direct access files and write text to screen
*
      call openem
*
*     label 10 is the subject of a backward pointing goto
*
   10 write(*,'(a\)')' Input nuclide (e.g.; Ba-137m or <Enter> to quit)-
     .> '
      read(*,'(bn, a7)') nuke
      if (lentrim(nuke) .ne. 0) then          ! branch to exit program
*
*        call nukeok to determine if nuke in data bases
*
         call nukeok (nuke, ok)
         if (.not. ok) then
            write(*,'(1x, ''Radionuclide '', (a), 
     :           '' is not in index file!'')') nuke(:lentrim(nuke))
            write(*,*)
            write(*,*)
            call pauseit
         else
*
*           lets get some factors
*
            call getcof(nuke, ipath)
*
*           call tablem to print factors to screen
*
           call tablem (ipath)
*
         endif
         call cls
*
*        compute activity of chain members as at user specified times
*        note label 20 is also the subject of a backward pointing goto
*
        write(*,*)
        write(*,*)'Compute Activity and Time-Integrated Activity Ratio' 
        write(*,*)
   20     write(*,'(a\)')' Input time (d) or zero to quit --> '
          read(*,'(bn, e10.0)') t
          if (t .ne. 0.) then
           write(*,'(5x,''Nuclide  T1/2        A(t)/Ao   intA/Ao(d)'')')
*
*           do over the maximum length of the retained chain members
*
           do ispec = 1, max(nint, next)
              call birch(ispec, t, a, aint)
              write(*,'(5x, a7, 1x, a8, a2, 1p2e12.5)') nucnam(ispec), 
     :             thalf(ispec), iu(ispec), a, aint
           end do
           goto 20
          else
*
*          go get another radionuclide
*          
           call cls
          go to 10        
          end if 
      endif
*
*     close all files
*
      call closem
      end
*-----------------------------------------------------------------------
*
      subroutine tablem (ipath)
*
*-----------------------------------------------------------------------
*     routine:  tablem
*     author:   k.f. eckerman
*     data:     12/2/99
*     purpose:  write dose coefficients to screen. 
*
      include 'pakparm.cmn'
      include 'fgr13pak.cmn'
      character*70 head
      character*25 chemform(10)
      character*8 t, ltrim, buffer
      logical ipath(3)
*
      if (ipath(1)) then
         do ispec = 1, nint
            ne = NLET(ispec)
            if (iflag(ispec, 1)) then
              do itype = 1, ntypes(ispec)
               head = 'Inhalation Absorbed Dose Coefficients (Gy/Bq)'
               do iage = 1, mage
                 call cls
                 ip = (80-lentrim(head)) / 2 - 4
                 call curpos(2, ip)
                 write(*,*) head(:lentrim(head)), ': ', namage(iage)
                 write(*,*)
                 t = ltrim(Thalf(ispec))
                 call curright(3)
                 write(buffer, '(1pe8.1)') f1inh(itype,ispec,iage)
                 if (nfinh(ispec) .eq. 0 .or. itype .le. 3) then
                   write(*,*) nucnam(ispec) // ' T1/2 = ' // 
     :                t(:lentrim(t)) // ' ' // iu(ispec) // ' Type: ' 
     :                // type(itype) // ' f_1 ='  // buffer
                 else
                   write(*,*) nucnam(ispec) // ' T1/2 = ' // 
     :                t(:lentrim(t)) // ' ' // iu(ispec) // ' Type: ' 
     :                // type(itype) // ' f_1 ='  // buffer //
     :                ' ' // chemfinh(ispec, itype - 3)
                 end if     
                 if (ne .eq. 1) then
                   write(*,*)
                   if (LETh(ispec, ne) .eq. 'L') then
                     write(*,'(/,t7,''Organ'',t19,''d(low)'',t32,
     :                 ''Organ'', t44,''d(low)'')')
                   else
                     write(*,'(/,t7,''Organ'',t19,''d(high)'',t32,
     :                ''Organ'', t44,''d(high)'')')
                   end if                     
                   write(*,'(2(6x,a9,2x,1pe8.2))') (organ(j), 
     :              (dfinh(itype,ispec,iage,j,ilet),ilet=1,ne),j=1,morg)
                 else
                   write(*,'(/,t6,''Organ'',t18,''d(low)'',t28,
     :              ''d(high)'',t40,''Organ'',t52,''d(low)'',t62,
     :              ''d(high)'')')
                   write(*,'(2(5x,a9,1p2e10.2))') (organ(j), 
     :              (dfinh(itype,ispec,iage,j,ilet),ilet=1,ne), 
     :               j=1,morg-2)
                   write(*,'(5x,a9,1p2e10.2,5x,a9,e10.2))')
     :              organ(morg-1),(dfinh(itype,ispec,iage,morg-1,ilet),
     :              ilet=1,ne),organ(morg), 
     :              dfinh(itype,ispec,iage,morg,1)
                 end if
                 call curpos (23, 62)
                 write(*,'('' Member:'',i2,''/'',i2)') ispec, nspec
                 call curpos (23, 1)
                 call pauseit
               end do                            ! End of age loop
               if (nfinh(ispec) .eq. 0 .or. itype .le. 3) then
                  head = 'Risk Coefficient (/Bq) for Inhalation Intakes'
               else
                  head = 'Risk Coefficient (/Bq) for Inhalation Intakes'
     :                   // ': ' // chemfinh(ispec, itype -3)             
               end if   
               call cls
               ip = (80-lentrim(head)) / 2
               call curpos(2, ip)
               write(*,*) head(:lentrim(head))
               write(*,*)
               t = ltrim(Thalf(ispec))
               call curright(3)
               write(*,*) nucnam(ispec) // ' T1/2 = ' // t(:lentrim(t)) 
     :              // ' ' // iu(ispec) // ' Type: ' // type(itype)
                write(*,*)
                write(*,'(t6,''Cancer   ''t19,''Mortality'',t32,
     :                 ''Morbidity'')')
                do ican = 1, mcan
                   write(*,'(5x,a9, 1p2e13.2)') cancer(ican), 
     :                   (rinh(itype,ispec,ican,k), k = 1, 2)
                end do
                call curpos (23,62)
                write(*,'('' Member:'',i2,''/'',i2)') ispec, nspec
                call curpos (23, 1)
                call pauseit
              end do                      ! end type  
            end if
         end do                                  ! End of species loop
      end if                                     ! End of ipath(1)
*
      if (ipath(2)) then
         do ispec = 1, nint
            ne = NLET(ispec)
            if (iflag(ispec, 2)) then
             do itype = 1, nfings(ispec)           
               head = 'Ingestion Absorbed Dose Coefficients (Gy/Bq)'
               do iage = 1, mage
                 call cls
                 ip = (80-lentrim(head)) / 2 - 4
                 call curpos(2, ip)
                 write(*,*) head(:lentrim(head)), ': ', namage(iage)
                 write(*,*)
                 t = ltrim(Thalf(ispec))
                 call curright(3)
                 write(buffer, '(1pe8.1)') f1ing(itype,ispec,iage)
                 if (nfing(ispec) .eq. 0) then
                   write(*,*) nucnam(ispec) // ' T1/2 = ' // 
     :                t(:lentrim(t)) // ' ' // iu(ispec) // ' f_1 ='
     :                // buffer 
                 else
                   write(*,*) nucnam(ispec) // ' T1/2 = ' // 
     :                t(:lentrim(t)) // ' ' // iu(ispec) // ' f_1 ='
     :                // buffer // ' ' // chemfing(ispec, itype)
                 end if                      
                 if (ne .eq. 1) then
                   write(*,*)
                   if (LETg(ispec, ne) .eq. 'L') then
                     write(*,'(/,t7,''Organ'',t19,''d(low)'',t32,
     :                 ''Organ'', t44,''d(low)'')')
                   else
                     write(*,'(/,t7,''Organ'',t19,''d(high)'',t32,
     :                ''Organ'', t44,''d(high)'')')
                   end if                     
                   write(*,'(2(6x,a9,2x,1pe8.2))') (organ(j), 
     :                (dfing(itype,ispec,iage,j,ilet),ilet=1,ne),
     :                j = 1, morg)
                 else
                   write(*,'(/,t6,''Organ'',t18,''d(low)'',t28,
     :              ''d(high)'',t40,''Organ'',t52,''d(low)'',t62,
     :              ''d(high)'')')
                     write(*,'(2(5x,a9,1p2e10.2))') (organ(j), 
     :                (dfing(itype,ispec,iage,j,ilet),ilet=1,ne),
     :                j = 1, morg-2)
                   write(*,'(5x,a9,1p2e10.2,5x,a9,e10.2))')
     :              organ(morg-1),(dfing(itype,ispec,iage,morg-1,ilet),
     :              ilet=1,ne),organ(morg), 
     :              dfinh(itype,ispec,iage,morg,1)
                 end if
                 call curpos (23, 62)
                 write(*,'('' Member:'',i2,''/'',i2)') ispec, nint
                 call curpos (23, 1)
                 call pauseit
               end do                            ! End of age loop
               if (nfing(ispec) .eq. 0) then
                 head = 'Risk Coefficient (/Bq) for Ingestion Intakes'
               else
                 head = 'Risk Coefficient (/Bq) for Ingestion Intakes'
     :                  // ': ' // chemfing(ispec, itype)
               end if       
               call cls
               ip = (80-lentrim(head)) / 2
               call curpos(2, ip)
               write(*,*) head(:lentrim(head))
               write(*,*)
               t = ltrim(Thalf(ispec))
               call curright(3)
               write(*,*) nucnam(ispec) // ' T1/2 = ' // t(:lentrim(t)) 
     :              // ' ' // iu(ispec)
               write(*,*)
               write(*,'(t19,''--- Drinking Water ---'',
     :                   t48,''--- Dietary Intake ---'')')
               write(*,'(t6,''Cancer   ''t19,''Mortality'',t32,
     :                 ''Morbidity'',t48,''Mortality'',t61,
     :                 ''Morbidity'')')
               do ican = 1, mcan
                  write(*,'(5x,a9,1p2e13.2,3x,2e13.2)') cancer(ican), 
     :                 ((ring(itype,ispec,ican,j,k), k = 1, 2), j=1,2)
               end do
               call curpos (23,62)
               write(*,'('' Member:'',i2,''/'',i2)') ispec, nint
               call curpos(23, 1)
               call pauseit
             end do  
            end if
         end do                                  ! End of species loop
      end if                                     ! End of ipath(1)
*
*     external coefficients
*   
      if (ipath(3)) then
         do ip = 1, mext
            do ispec = 1, next
               if (ip .eq. 1) then
                  head = 'External Absorbed Dose Coefficients:' //
     :                   ' Submersion (Gy per Bq-s/m^3)'
               elseif (ip .eq. 2) then
                  head = 'External Absorbed Dose Coefficients:' //
     :                   ' Ground Surface (Gy per Bq-s/m^2)'
               else
                  head = 'External Absorbed Dose Coefficients:' //
     :                   ' Soil Volume (Gy per Bq-s/m^3)'
               end if
               if (iflag(ispec, 3)) then
                 call cls
                 ipx = (80-lentrim(head)) / 2
                 call curpos(2, ipx)
                 write(*,*) head(:lentrim(head))
                 write(*,*)
                 t = ltrim(Thalf(ispec))
                 call curright(3)
                 write(*,*) nucnam(ispec) //
     :                    ' T1/2 = ' // t(:lentrim(t)) // ' ' //
     :                    iu(ispec) 
                 write(*,*)
                 write(*,'(/,t7,''Organ'',t20,''h_T'',t32,''Organ'',t45,
     :                   ''h_T'')')
                 write(*,'(2(6x,a9,2x,1pe8.2))') (organ(j), 
     :              dfext(ispec, ip, j), j = 1, morg)
                 call curpos (23,62)
                 write(*,'('' Member:'',i2,''/'',i2)') ispec, nint
                 call curpos (23, 1)
                 call pauseit
                 head = 'External Risk Coefficient'
                 if (ip .eq. 1) then
                    head = 'Risk Coefficients for' //
     :                     ' Submersion (/Bq-s/m^3)'
                 elseif (ip .eq. 2) then
                    head = 'Risk Coefficients for' //
     :                     ' Ground Surface (/Bq-s/m^2)'
                 else
                    head = 'Risk Coefficients for' //
     :                     ' Soil Volume (/Bq-s/m^3)'
                 end if
                 call cls
                 ipx = (80-lentrim(head)) / 2
                 call curpos(2, ipx)
                 write(*,*) head(:lentrim(head))
                 write(*,*)
                 t = ltrim(Thalf(ispec))
                 call curright(3)
                 write(*,*) nucnam(ispec) //
     :                    ' T1/2 = ' // t(:lentrim(t)) // ' ' //
     :                    iu(ispec)
                 write(*,*)
                 write(*,'(t6,''Cancer   ''t19,''Mortality'',t32,
     :                 ''Morbidity'')')
                 do ican = 1, mcan
                    write(*,'(5x,a9, 1p2e13.2)') cancer(ican), 
     :                   (rext(ispec, ican, ip, k), k = 1, 2)
                 end do
                 call curpos (23,62)
                 write(*,'('' Member:'',i2,''/'',i2)') ispec, nint
                 call curpos (23, 1)
                 call pauseit
               end if
            end do                               ! End of ispec loop
         end do                                  ! End of ipath loop
      end if                                     ! End of ipath(1)
      return
      end
