************************************************************************ * * PROGRAM IDCF * * ************************************************************************ * * * * Steve Fetter * Version 1.0 for MSDOS * * Subcontractor to EG&G, Idaho * Date: 06/12/90 * * Idaho Falls, ID * Copyright (c) 1990 * * * * ************************************************************************ * * * Program IDCF calculates the dose, delivered over a variety of * * time intervals, to various organs resulting from the one-time * * inhalation or ingestion of one curie of a radionuclide. IDCF * * uses the models presented in ICRP Publication No. 30 to model * * the retention and distribution of radionuclides in the lung, GI * * tract, bone, and other organs. This code is documented in * * Steve Fetter, "Internal Dose Conversion Factors for 19 Target * * Organs and 9 Irradiation Times and External Dose-Rate Conversion * * Factors for 21 Target Organs for 259 Radionuclides Produced in * * Potential Fusion-Reactor Materials," EGG-FSP-8036 (Idaho Falls: * * EG&G, March 1988). There have been only cosmetic changes to the * * code since March 1988. * * * * There are three special cases: inert gases, alkaline earths, * * and highly soluble gases. As argued by the ICRP, the internal * * dose from a lung-full of inert gas will be small compared with * * the dose from external radiation. The only exception is Ar-37, * * which emits very-low-energy electrons and photons; such cases * * are ignored here. Alkaline-earth metabolism is described by a * * power-function model that is very difficult to treat * * analytically, so the power function has been approximated by a * * sum of exponentials and is treated the same as other compounds. * * Highly soluble gases (H2O, CO2) are assumed to be instantaneously * * distributed uniformly throughout the body, without passing * * through the lung or GI. * * * * The ICRP model has also been implemented at Oak Ridge; this * * model is described in ORNL/NUREG/TM-190. There is one major * * difference between that implementation and the present program * * related to the biological behavior of radioactive daughters: the * * ORNL model assumes that daughters behave biologically according * * to their atomic number; the ICRP assumes that they behave * * according to the atomic number of the parent radionuclide (except * * in the case of radionuclides decaying to an isotope of iodine). * * Both the ORNL and ICRP models assume that inert gases produced by * * radioactive decay in the body escape with a half-life of 2 hours. * * In my model, daughters behave like the parent without exception. * * This will lead to small differences in the dose computed for * * Te-129m, Te-129, Te-131m, Te-131, Te-132, I-133, and I-135. * * Since none of these radionuclides are important in fusion * * applications, these differences are relatively unimportant. * * * * Dose to the lung: the ICRP neglects the dose to the * * nasopharyngeal of the lung because "for most particle sizes it is * * usually small in comparison with the dose received by other * * regions." This is not true in the case of radionuclides with * * short half-lives; in fact, the total dose can be underestimated * * by as much as 50%. I have followed the ORNL model in this case * * and included the dose to the N-P region. * * * * And now, a word from our sponsor: * * * * "This code was developed with funds from the U.S. Department of * * Energy. Neither DOE nor persons acting on behalf of DOE makes * * any warrenty or other representation, express or implied, with * * respect to the accuracy, completeness, or usefulness of any * * information furnished hereunder; that the use of any such infor- * * mation will not infringe privately owned rights; that the * * services, materials, or information furnished hereunder will not * * result in injury or damage when used for any such purpose; nor * * that the services, materials or information furnished hereunder * * will accomplish the intended results or are safe for any purpose * * including the intended purpose." * * * * Input Files * * ----- ----- * * decaylib: decay-scheme data file; attached to unit 1; * * biolib: biological data file; attached to unit 2; * * saflib: specific absorbed fractions file; attached to unit 3. * * * * Output Files * * ------ ----- * * IDCF.1,2,3,4: the results in text file format; * * IDCF-1.DAT, IDCF-2.DAT: results in data file format. * * * ************************************************************************ logical surface integer dk, bio, afu, key, con, out(6), nuclide, + intake, time, target, source, daughter, ndaughters, z, + nsources, ntargets, ntimes, nboxes(99), organ(15,99), + zai(0:14), nintakes, nnuclides, t, box, organ_p(15), + start, end, i, j, nbranches, ndaughter(3), + daughter_zais(3,5) character*1 lung_class(99), bone_class(99) character*4 time_name(9) character*8 nuclide_name, target_name(22), old_name/'nonamyet'/ character*12 outfile(6) real branch, branch_sum, dcf(22,9,2), see(19,19,0:14), + lambda_rad(0:14), U(0:14,19,9,2), delta_t, f1(99), + bio_halflife(15,99), fraction(15,99), int_time(9), + Er(12), saf(12,19,19), Bqxgxrad_per_CixMeV, + bio_halflife_p(15), fraction_p(15), u_wb, sec_per_day, + sec_per_year, branchfrc(3), decay_c parameter (dk = 1, bio = 2, afu = 3, key = 5, con = 6, + Bqxgxrad_per_CixMeV = 592.81004, + nsources = 19, ntargets = 19, nintakes = 2, ntimes = 9, + sec_per_day = 24.0 * 3600.0, sec_per_year = 3.156E+07) data (out(i), i = 1, 6) + / 7, 8, 9, 10, 11, 12 / data (outfile(i), i = 1, 6) + / 'IDCF.1', 'IDCF.2', 'IDCF.3', 'IDCF.4', + 'IDCF-1.DAT', 'IDCF-2.DAT' / data (int_time(t), t = 1, ntimes) + / 2.0, 7.0, 30.0, 1.0, 10.0, 20.0, 30.0, 40.0, 50.0 / data (time_name(t), t = 1, ntimes) + / ' 2 d', ' 7 d', '30 d', + ' 1 y', '10 y', '20 y', '30 y', '40 y', '50 y' / call under0('true') call ovefl('true') * ** Print sign on message: write(con,1000) * ** Get integration times: do t = 1, 3 int_time(t) = int_time(t) * sec_per_day end do do t = 4, ntimes int_time(t) = int_time(t) * sec_per_year end do * ** Read in data and initialize data file: call INITIALIZE (dk, bio, afu, key, con, out, outfile, + ntargets, nsources, f1, lung_class, + bone_class, nboxes, organ, bio_halflife, + fraction, er, saf, target_name) * ** Start timer: call TIMER(start) * ** Begin processing radionuclide data: write (con,1001) read (dk,2000) nnuclides write (out(5),4002) nnuclides write (out(6),4002) nnuclides do nuclide = 1, nnuclides * ** Initialize sums: branch_sum = 0.0 do intake = 1, nintakes do time = 1, ntimes do target = 1, (ntargets+3) dcf(target,time,intake) = 0.0 end do end do end do * ** Begin main loop: do while (branch_sum .lt. 0.99) read (dk,2001) nuclide_name, branch, ndaughters write(con,1002) nuclide_name branch_sum = branch + branch_sum if (nuclide_name .eq. old_name) then nbranches = nbranches + 1 else old_name = nuclide_name nbranches = 1 ndaughter(2) = 0 ndaughter(3) = 0 branchfrc(2) = 0.0 branchfrc(3) = 0.0 do 1, i = 1,3 do 2, j = 1,5 daughter_zais(i,j) = 0 2 continue 1 continue end if ndaughter(nbranches) = ndaughters branchfrc(nbranches) = branch/10 read (dk,2002) decay_c backspace (dk) * ** Calculate the specific effective energy (MeV/g/dis)x(rem/rad) * ** absorbed in each target organ from decays of each daughter in * ** source organ: call SEE_CALC (dk, ndaughters, bone_class, nbranches, + Er, saf, ntargets, nsources, + see, lambda_rad, z, zai, surface, + daughter_zais) * ** Only continue of the parent is not inert: if ((z .ne. 2) .and. (z .ne. 10) .and. + (z .ne. 18) .and. (z .ne. 36) .and. + (z .ne. 54) .and. (z .ne. 86)) then * ** Calculate the number of decays of each daughter in each source * ** organ over each integration time, for intake by inhalation and * ** by ingestion (dis/Bq): do box = 1, nboxes(z) organ_p(box) = organ(box,z) bio_halflife_p(box) = bio_halflife(box,z) fraction_p(box) = fraction(box,z) end do call U_CALC (ndaughters, lambda_rad, f1(z), + lung_class(z), nboxes(z), organ_p, + bio_halflife_p, fraction_p, surface, + nsources, int_time, ntimes, nintakes, U) * ** Now sum over all sources organ and daughters to get the dose * ** conversion factors for each target organ, integration time, * ** and intake route (rem/Ci): do intake = 1, nintakes do time = 1, ntimes do target = 1, ntargets do source = 1, nsources do daughter = 0, ndaughters dcf(target,time,intake) = + branch * + Bqxgxrad_per_CixMeV * + U(daughter,source,time,intake) * + see(target,source,daughter) + + dcf(target,time,intake) end do ! daughter end do ! source end do ! target do daughter = 0, ndaughters u_wb = 0.0 do source = 1, nsources u_wb=U(daughter,source,time,intake)+ + u_wb end do dcf(20,time,intake) = branch * + Bqxgxrad_per_CixMeV * + u_wb * see(19,19,daughter) + + dcf(20,time,intake) end do end do ! time end do ! intake end if ! inert end do ! branch call EDE_CALC (ntargets, ntimes, nintakes, dcf) call OUTPUT (out, nuclide_name, zai, decay_c, nbranches, + branchfrc, ndaughter, daughter_zais, time_name, + dcf, ntimes, nintakes, target_name, nuclide) end do ! nuclides call TIMER(end) delta_t = 0.01 * REAL(end - start) / 60.0 write (con,1003) delta_t close (unit = dk, status = 'keep') do i = 1, 6 close (unit = out(i), status = 'keep') end do stop 'adios' 1000 format (//16X,'*************************************************' + /16X,'* IDCF *' + /16X,'* A Program by Steve Fetter to Calculate *' + /16X,'* Internal Dose Conversion Factors *' + /16X,'*************************************************') 1001 format ( /16X, 'Calculating Dose Factors...') 1002 format ('+', A8) 1003 format ( /16X, 'IDCF ran in ', F5.1, ' minutes.', + /16X, 'The results are in IDCF.1,.2,.3,.4, ' + 'IDCF-1.DAT, and IDCF-2.DAT'/) 2000 format (6X, I3) 2001 format (/ A8, 2X, F6.4, 2X, I1) 2002 format (8X, E9.3E2) 4002 format (I3, ' = number of radionuclides in file') end ! idcf *************************** LEVEL 2 ROUTINES *************************** * * ************************************************************************ * * SUBROUTINE INITIALIZE (dk, bio, afu, key, con, + out, outfile, ntargets, nsources, + f1, lung_class, bone_class, nboxes, + organ, bio_halflife, fraction, er, + saf, target_name) * * ************************************************************************ * * * This subroutine reads the metabolic data file and the specific * * absorbed fraction data file. * * * * Input * * ----- * * dk, bio, afu, key, con, out: I/O units for radionuclide decay * * file, metabolic data file, specific absorbed fraction file, * * console input, console output, and output files. * * ntargets, nsources: total number of target and source organs. * * * * Output * * ------ * * f1, lung_class, bone_class, nboxes: the fraction transferred from * * the GI to the blood, the lung clearance class, the bone dep- * * osition class, and the number of compartments in the model * * for each element; * * organ, bio_halflife, fraction: for each element and each compart- * * ment, the organ index, biological halflife, and fraction * * transferred from the blood to the compartment; * * er, saf: the reference energies, and the specific absorbed frac- * * tion for each energy, target organ, and source organ; * * outfile: file name of output files. * * * ************************************************************************ logical exist integer dk, bio, afu, key, con, out(6), z, organ(15,99), + nelements, box, nboxes(99), source, nsources, + target, ntargets, energy, element, i character*1 lung_class(99), bone_class(99), char character*8 target_name(22), dk_version, dy_mo_yr, hr_min_sec character*12 decaylib, biolib, outfile(6) real f1(99), bio_halflife(15,99), fraction(15,99), + Er(12), saf(12,ntargets,nsources) * ** Open the decay data file: * ** ************************ write (con,1000) read (key,1001) decaylib inquire (file = decaylib, exist = exist) if (exist) then open (unit = dk, file = decaylib, status = 'old', + action = 'read') else write (con,1002) decaylib stop 'Sorry' end if * ** Swallow comments: read (dk,2000) dk_version char = '*' do while (char .eq. '*') read (dk,2001) char end do * ** Open the metabolic data file: * ** **************************** write (con,1003) read (key,1001) biolib inquire (file = biolib, exist = exist) if (exist) then open (unit = bio, file = biolib, status = 'old', + action = 'read') else write (con,1002) biolib stop 'Sorry' end if * ** Swallow the comments: char = '*' do while (char .eq. '*') read (bio,2001) char end do * ** Now read the data in: read (bio,2002) nelements do element = 1, nelements read (bio,2003) z, f1(z), lung_class(z), bone_class(z) read (bio,2002) nboxes(z) do box = 1, nboxes(z) read (bio,2004) organ(box,z), bio_halflife(box,z), + fraction(box,z) end do ! boxes end do ! elements close (unit = bio, status = 'keep') * ** Open the specific absorbed fraction data file: * ** ********************************************* inquire (file = 'saflib.dat', exist = exist) if (exist) then open (unit = afu, file = 'saflib.dat', status = 'old', + action = 'read') else write (con,1002) 'saflib.dat' stop 'Sorry' end if * ** Swallow comments char = '*' do while (char .eq. '*') read (afu,2001) char end do * ** Now read the reference energies: read (afu,3000) (Er(energy), energy = 1, 12) do energy = 1, 12 Er(energy) = log(Er(energy)) end do * ** Now read the specific absorbed fractions do source = 1, nsources read (afu,3001) target_name(1), + (saf(energy,1,source), energy = 1, 12) do target = 2, ntargets read (afu,3002) target_name(target), + (saf(energy,target,source), energy = 1, 12) end do end do target_name(19) = 'others' target_name(20) = 'whl.body' target_name(21) = 'tot.body' target_name(22) = ' EDE' do source = 1, nsources do target = 1, ntargets do energy = 1, 12 if (saf(energy,target,source) .gt. 1.0E-35) then saf(energy,target,source) = + log(saf(energy,target,source)/1000.0) else saf(energy,target,source) = -88.74 end if end do end do end do close (unit = afu, status = 'keep') * ** Now open the output files: * ** ************************** exist = .false. do i = 1, 6 inquire (file = outfile(i), exist = exist) end do if (exist) then write (con,1004) read (key,1005) char if ((char .eq. 'Y') .or. (char .eq. 'y')) then call SYSTEM ('del idcf.?') call SYSTEM ('del idcf-?.dat') else stop 'rename files and restart program' end if end if do i = 1, 6 open (unit = out(i), file = outfile(i), status = 'new') end do call DATE (dy_mo_yr) call TIME (hr_min_sec) write (out(5),4000) dy_mo_yr, hr_min_sec, dk_version, + (target_name(target), target = 3, 18), + (target_name(target), target = 21, 22) write (out(6),4001) dy_mo_yr, hr_min_sec, dk_version, + (target_name(target), target = 3, 18), + (target_name(target), target = 21, 22) return 1000 format (/16X,'What is the name of the radionuclide file? ') 1001 format (A12) 1002 format (/16X,'Error: file ', A12, ' doesn''t exist') 1003 format (/16X,'What is the name of the metabolic data file? ') 1004 format (/16X,'IDCF output files already exist.' + /16X,' Do you want to erase them (Y/N)? ') 1005 format (A1) 2000 format (//63X, A5) 2001 format (A1) 2002 format (I2) 2003 format (/ I2, 4X, E8.2, X, A1, X, A1) 2004 format (I2, E9.2, E9.2) 3000 format (/ 14X, 12(4X, F5.3)) 3001 format (//////// A8, 7X, 12(2X, E7.1)) 3002 format (A8, 7X, 12(2X, E7.1)) 4000 format (A8, X, A8, ' DECAYLIB version ', A8 / 162('*') + / '* This file contains internal dose ' + 'conversion factors (IDCFs) for inhalation, ' + 'for use by the DOSE code. IDCFs are given for ' + 'the following 18 target organs: *' + / '*', X, A7, 16(X, A8), X, A7, '*' / 162('*')) 4001 format (A8, X, A8, ' DECAYLIB version ', A8 / 162('*') + / '* This file contains internal dose ' + 'conversion factors (IDCFs) for ingestion, ' + 'for use by the DOSE code. IDCFs are given for ' + 'the following 18 target organs: *' + / '*', X, A7, 16(X, A8), X, A7, '*' / 162('*')) end ! intialize ************************************************************************ * * SUBROUTINE SEE_CALC (dk, ndaughters, bone_class, nbranches, + Er, saf, ntargets, nsources, + see, lambda_rad, z, zai, surf, + daughter_zais) * * ************************************************************************ * * * Subroutine SEE_CALC computes the specific effective energy * * (MeV/g/dis) absorbed in each target organ from radiations of * * each daughter in each source organ, using the rules and given * * in sections 4, 5, 6, and 7 of ICRP Publication No. 30, Part 1. * * All decay energies except protons, neutrons, fission fragments, * * and neutrinos are considered, because no radionuclides in the * * decay library exhibited these decay modes. * * * * Input * * ----- * * dk: the logical unit associated with DECAYLIB. * * ndaughters: the number of radioactive daughter products in the * * decay chain; * * bone_class: bone dosimetry class; * * Er: the reference energies for the array saf; * * saf: the array of specific absorbed energies, saf(e,t,s), for * * each source organ s, target organ t, and energy e; * * ntargets, nsources: the number of target and source organs. * * * * Output * * ------ * * see: the array of specific effective energies, see(t,s,j), * * absorbed in each target organ t per disintegration of * * daughter j in source organ s (MeV/g/dis); * * lambda_rad: the array of radioactive decay constants of each * * daughter (1/s); * * z, zai: the atomic number and identifier of the parent. * * * ************************************************************************ logical surf integer dk, ndaughters, ntargets, nsources, z, i, j, k, + nalphas, nelectrons, nphotons, t, s, zai(0:14), + daughter_zais(3,5), nbranches character*1 bone_class(99) real lambda_rad(0:ndaughters), Er(12), saf(12,ntargets,nsources), + see(ntargets,nsources,0:ndaughters), Q_alpha, Q_recoil, + atomic_weight, e, f, AF, afp(12) parameter (Q_alpha = 20.0, Q_recoil = 20.0) * ** Begin main loop... do j = 0, ndaughters read (dk,1000) zai(j), lambda_rad(j) if (j .eq. 0) then z = zai(0) / 10000 else daughter_zais(nbranches,j) = zai(j) end if atomic_weight = REAL (INT((zai(j) - 10000 * + (zai(j)/10000)) / 10)) * ** Determine the bone class; if bone_class='N' and the halflife is * ** less than 15 days, then it's surface deposited: if ((bone_class(z) .eq. 'S') .or. + ((bone_class(z) .eq. 'N') .and. + (lambda_rad(j) .ge. 5.348E-07))) then surf = .true. else surf = .false. end if * ** Zero sums: do s = 1, nsources do t = 1, ntargets see(t,s,j) = 0.0 end do end do * ** Compute the alpha and atomic recoil energies... read (dk,2000) nalphas do i = 1, nalphas read (dk,2001) e, f do s = 1, nsources do t = 1, ntargets see(t,s,j) = e * f * (Q_alpha * + AF(t,s,'a',e,surf,Er,afp) + + Q_recoil * (4.0/atomic_weight) * + AF(t,s,'r',e,surf,Er,afp)) + + see(t,s,j) end do ! targets end do ! sources end do ! alphas * ** Compute the electron energy: do k = 1, 3 read (dk,2000) nelectrons do i = 1, nelectrons if (k .eq. 3) then read (dk,2001) e, f else read (dk,2002) f, e end if do s = 1, nsources do t = 1, ntargets see(t,s,j) = e * f * + AF(t,s,'e',e,surf,Er,afp) + + see(t,s,j) end do ! targets end do ! sources end do ! electrons end do ! k * ** Compute the photon energy: read (dk,2000) nphotons do i = 1, nphotons read (dk,2001) e, f do s = 1, nsources do t = 1, ntargets do k = 1, 12 afp(k) = saf(k,t,s) end do see(t,s,j) = e * f * AF(t,s,'p',e,surf,Er,afp) + + see(t,s,j) end do ! targets end do ! sources end do ! photons end do ! daughter return 1000 format (I7, X, E9.3) 2000 format (I2) 2001 format (5X, F9.7, 2X, F9.7) 2002 format (5X, 9X, 2X, F9.7 / 5X, F9.7) end ! see_calc ************************************************************************ * * SUBROUTINE U_CALC (ndaughters, lambda_rad, f1, + lung_class, nboxes, organ, bio_halflife, + fraction, surface, nsources, time, ntimes, + nintakes, U) * * ************************************************************************ * * * Subroutine U_CALC calculates the number of disintegrations of each * * radioactive daughter in each source organ for each integration * * time, using a superposition of linear chains. There are 15 * * possible chains, of which 2 are ingestion routes, and 13 are * * inhalation routes. The solution as a function of integration * * time is exact, and is much less time consuming than eigenvalue * * methods. As far as I know, this method is original. * * * * Input * * ----- * * ndaughters: the number of radioactive daughters; * * lambda_rad: the array of their radioactive decay constants; * * f1, lung_class: the fraction transferred from the GI to the blood * * and the lung clearance class of the parent; * * nboxes: the number of organ compartments in the metabolic model; * * organ: the organ index of each compartment; * * bio_halflife: the biological halflife of each compartment; * * fraction: the fraction of the parent transferred from the blood * * to the compartment; * * nsources, time, ntimes, nintakes: the number of source organs, * * the array of integration times and the number of times, and * * the number of intake modes. * * * * Output * * ------ * * U: the array of disintegrations of each daughter in each source * * organ for each integration time and intake mode. * * * * * * To facilitate computation, the removal times, etc. have been * * subscripted, so that compartments a,b,c...j are now 1,2,3...10. * * Similarily, the lung clearance class is also a subscript, with * * D,W,Y becoming 1,2,3. * * * * ICRP lung model: * * ================================================================ * * Class * * ----------------------------------------- * * Compart- D=1 W=2 Y=3 * * Region ment T F T F T F * * ---------------------------------------------------------------- * * N-P a=1 0.01 0.5 0.01 0.1 0.01 0.01 * * (Dnp=0.30) b=2 0.01 0.5 0.40 0.9 0.40 0.99 * * ---------------------------------------------------------------- * * T-B c=3 0.01 0.95 0.01 0.5 0.01 0.01 * * (Dtb=0.08) d=4 0.2 0.05 0.2 0.5 0.2 0.99 * * ---------------------------------------------------------------- * * P e=5 0.5 0.8 50 0.15 500 0.05 * * (Dp=0.25) f=6 n.a. n.a. 1.0 0.4 1.0 0.4 * * g=7 n.a. n.a. 50 0.4 500 0.4 * * h=8 0.5 0.2 50 0.05 500 0.15 * * ---------------------------------------------------------------- * * L i=9 0.5 1.0 50 1.0 1000 0.9 * * j=10 n.a. n.a. n.a. n.a. infin. 0.1 * * ================================================================ * * Removal half-times T measured in days. * * n.a. = not applicable * * * * Special organs: transfer compartment (20), lung (9), ST (5), * * SI (4), ULI (6), LLI (3), and whole body (19). * * * ************************************************************************ logical unused(19), surface integer ndaughters, nboxes, organ(nboxes), nsources, ntimes, + nintakes, i, j, s, t, lng_class, box, org character*1 lung_class real lambda_rad(0:ndaughters), f1, bio_halflife(nboxes), + fraction(nboxes), U(0:14,nsources,ntimes,nintakes), + time(ntimes), fraction_lung(10,4), lambda_lung(10,4), + lambda_gi(4), lambda_tc, lambda_b, V(0:14,6), sec_per_day, + lambda_bio(6), transfer(6), sec_per_hour, D_np, D_tb, D_p, + tau(4), u_tc(0:14,2), mass(19), + lambda_organ, unused_mass parameter (sec_per_hour = 3600.0) parameter (sec_per_day = sec_per_hour * 24.0) * ** The fraction of activity of an aerosol which is deposited in * ** the N-P, T-B, and P regions for an Activity Median Aerodynamic * ** Diameter (AMAD) of 1.0E-6 m, from ORNL/NUREG/TM-190: parameter (D_np = 0.31, D_tb = 0.08, D_p = 0.249) data ((fraction_lung(i,j), i = 1, 10), j = 1, 4) + / 0.5, 0.5, 0.95, 0.05, 0.8, 0.0, 0.0, 0.2, 1.0, 0.0, + 0.1, 0.9, 0.5, 0.5, 0.15, 0.4, 0.4, 0.05, 1.0, 0.0, + 0.01, 0.99, 0.01, 0.99, 0.05, 0.4, 0.4, 0.15, 0.9, 0.1, +3.22580645, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 / data ((lambda_lung(i,j), i = 1, 10), j = 1, 4) + / 8.02E-04, 8.02E-04, 8.02E-04, 4.01E-05, 1.60E-05, + 8.02E-06, 1.60E-05, 1.60E-05, 1.60E-05, 0.0, + 8.02E-04, 2.01E-05, 8.02E-04, 4.01E-05, 1.60E-07, + 8.02E-06, 1.60E-07, 1.60E-07, 1.60E-07, 0.0, + 8.02E-04, 2.01E-05, 8.02E-04, 4.01E-05, 1.60E-08, + 8.02E-06, 1.60E-08, 1.60E-08, 8.02E-09, 0.0, + 1.00E-03, 1.00E-03, 1.00E-03, 1.00E-03, 1.00E-03, + 1.00E-03, 1.00E-03, 1.00E-03, 1.00E-03, 1.00E-03 / data (tau(i), i = 1, 4) + / 1.0, 4.0, 13.0, 24.0 / data (mass(s), s = 1, 19) + / 14.0, 1500.0, 135.0, 400.0, 250.0, 220.0, 310.0, + 1800.0, 1000.0, 28000.0, 11.0, 100.0, 1500.0, 5000.0, + 180.0, 35.0, 20.0, 200.0, 29325.0 / * ** Set everything equal to zero to start: if (surface) then mass(14) = 120.0 mass(19) = 34205.0 else mass(14) = 5000.0 mass(19) = 29325.0 end if do i = 1, nintakes do t = 1, ntimes do s = 1, nsources do j = 0, ndaughters U(j,s,t,i) = 0.0 end do end do end do end do * ** Get the lung clearance class index: if (lung_class .eq. 'D') then lng_class = 1 else if (lung_class. eq. 'W') then lng_class = 2 else if (lung_class. eq. 'Y') then lng_class = 3 else if (lung_class. eq. 'N') then lng_class = 4 end if * ** Get the decay constants for the GI tract: do i = 1, 4 if (lung_class .eq. 'N') then lambda_gi(i) = 1.00E-03 else lambda_gi(i) = 1.0 / (sec_per_hour * tau(i)) end if end do * ** Get the decay constant of the transfer compartment: lambda_tc = log(2.0) / (sec_per_day * 0.25) do box = 1, nboxes if (organ(box) .eq. 20) then lambda_tc = log(2.0) / (sec_per_day * bio_halflife(box)) end if end do if (lung_class .eq. 'N') then lambda_tc = 1.00E-03 end if * ** Get the decay constant for transfer from the SI to the TC: if (f1 .eq. 1.0) then f1 = 0.95 end if lambda_b = f1 * lambda_gi(2) / (1.0 - f1) * ** Figure out what organs aren't explicitly used in the model: do s = 1, 19 unused(s) = .true. end do do box = 1, nboxes if (organ(box) .le. 18) then unused(organ(box)) = .false. end if end do unused_mass = 0.0 do s = 1, 19 if (unused(s)) then unused_mass = mass(s) + unused_mass end if end do * ** Begin main loop: * ** *************** do t = 1, ntimes do i = 1, nintakes do j = 0, ndaughters u_tc(j,i) = 0.0 end do end do * ** First do all pathways leading to the transfer compartment: do box = 1, nboxes * ** If it's the transfer compartment, then give it a very small * ** halflife so that it's similar to direct excretion: org = organ(box) if (org .eq. 20) then lambda_organ = log(2.0) / (sec_per_day * 0.01) else lambda_organ = log(2.0) / + (sec_per_day * bio_halflife(box)) end if if ((org .eq. 20) .or. (org .eq. 99)) then org = 19 end if * ** The lung(a)-->tc-->organ route: lambda_bio(1) = lambda_lung(1,lng_class) lambda_bio(2) = lambda_tc lambda_bio(3) = lambda_organ transfer(1) = D_np * fraction_lung(1,lng_class) transfer(2) = 1.0 transfer(3) = fraction(box) call CHAIN (3,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + U(j,9,t,1) u_tc(j,1) = V(j,2) + u_tc(j,1) end if U(j,org,t,1) = V(j,3) + U(j,org,t,1) end do * ** The lung(c)-->tc-->organ route: lambda_bio(1) = lambda_lung(3,lng_class) lambda_bio(2) = lambda_tc lambda_bio(3) = lambda_organ transfer(1) = D_tb * fraction_lung(3,lng_class) transfer(2) = 1.0 transfer(3) = fraction(box) call CHAIN (3,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + U(j,9,t,1) u_tc(j,1) = V(j,2) + u_tc(j,1) end if U(j,org,t,1) = V(j,3) + U(j,org,t,1) end do * ** The lung(e)-->tc-->organ route: lambda_bio(1) = lambda_lung(5,lng_class) lambda_bio(2) = lambda_tc lambda_bio(3) = lambda_organ transfer(1) = D_p * fraction_lung(5,lng_class) transfer(2) = 1.0 transfer(3) = fraction(box) call CHAIN (3,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + U(j,9,t,1) u_tc(j,1) = V(j,2) + u_tc(j,1) end if U(j,org,t,1) = V(j,3) + U(j,org,t,1) end do * ** The lung(h)-->lung(i)-->tc-->organ route: lambda_bio(1) = lambda_lung(8,lng_class) lambda_bio(2) = lambda_lung(9,lng_class) lambda_bio(3) = lambda_tc lambda_bio(4) = lambda_organ transfer(1) = D_p * fraction_lung(8,lng_class) transfer(2) = fraction_lung(9,lng_class) transfer(3) = 1.0 transfer(4) = fraction(box) call CHAIN (4,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + V(j,2) + U(j,9,t,1) u_tc(j,1) = V(j,3) + u_tc(j,1) end if U(j,org,t,1) = V(j,4) + U(j,org,t,1) end do * ** The lung(b)-->stomach-->small intestine-->tc-->organ route: lambda_bio(1) = lambda_lung(2,lng_class) lambda_bio(2) = lambda_gi(1) lambda_bio(3) = lambda_b + lambda_gi(2) lambda_bio(4) = lambda_tc lambda_bio(5) = lambda_organ transfer(1) = D_np * fraction_lung(2,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = f1 transfer(5) = fraction(box) call CHAIN (5,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + U(j,9,t,1) U(j,5,t,1) = V(j,2) + U(j,5,t,1) U(j,4,t,1) = V(j,3) + U(j,4,t,1) u_tc(j,1) = V(j,4) + u_tc(j,1) end if U(j,org,t,1) = V(j,5) + U(j,org,t,1) end do * ** The lung(d)-->stomach-->small intestine-->tc-->organ route: lambda_bio(1) = lambda_lung(4,lng_class) lambda_bio(2) = lambda_gi(1) lambda_bio(3) = lambda_b + lambda_gi(2) lambda_bio(4) = lambda_tc lambda_bio(5) = lambda_organ transfer(1) = D_tb * fraction_lung(4,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = f1 transfer(5) = fraction(box) call CHAIN (5,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + U(j,9,t,1) U(j,5,t,1) = V(j,2) + U(j,5,t,1) U(j,4,t,1) = V(j,3) + U(j,4,t,1) u_tc(j,1) = V(j,4) + u_tc(j,1) end if U(j,org,t,1) = V(j,5) + U(j,org,t,1) end do * ** The lung(f)-->lung(d)-->stomach-->SI-->tc-->organ route: lambda_bio(1) = lambda_lung(6,lng_class) lambda_bio(2) = lambda_lung(4,lng_class) lambda_bio(3) = lambda_gi(1) lambda_bio(4) = lambda_b + lambda_gi(2) lambda_bio(5) = lambda_tc lambda_bio(6) = lambda_organ transfer(1) = D_p * fraction_lung(6,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 transfer(5) = f1 transfer(6) = fraction(box) call CHAIN (6,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + V(j,2) + U(j,9,t,1) U(j,5,t,1) = V(j,3) + U(j,5,t,1) U(j,4,t,1) = V(j,4) + U(j,4,t,1) u_tc(j,1) = V(j,5) + u_tc(j,1) end if U(j,org,t,1) = V(j,6) + U(j,org,t,1) end do * ** The lung(g)-->lung(d)-->stomach-->SI-->tc-->organ route: lambda_bio(1) = lambda_lung(7,lng_class) lambda_bio(2) = lambda_lung(4,lng_class) lambda_bio(3) = lambda_gi(1) lambda_bio(4) = lambda_b + lambda_gi(2) lambda_bio(5) = lambda_tc lambda_bio(6) = lambda_organ transfer(1) = D_p * fraction_lung(7,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 transfer(5) = f1 transfer(6) = fraction(box) call CHAIN (6,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,9,t,1) = V(j,1) + V(j,2) + U(j,9,t,1) U(j,5,t,1) = V(j,3) + U(j,5,t,1) U(j,4,t,1) = V(j,4) + U(j,4,t,1) u_tc(j,1) = V(j,5) + u_tc(j,1) end if U(j,org,t,1) = V(j,6) + U(j,org,t,1) end do * ** The stomach-->small intestine-->tc-->organ ingestion route: lambda_bio(1) = lambda_gi(1) lambda_bio(2) = lambda_b + lambda_gi(2) lambda_bio(3) = lambda_tc lambda_bio(4) = lambda_organ transfer(1) = 1.0 transfer(2) = 1.0 transfer(3) = f1 transfer(4) = fraction(box) call CHAIN (4,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters if (box .eq. 1) then U(j,5,t,2) = V(j,1) + U(j,5,t,2) U(j,4,t,2) = V(j,2) + U(j,4,t,2) u_tc(j,2) = V(j,3) + u_tc(j,2) end if U(j,org,t,2) = V(j,4) + U(j,org,t,2) end do end do ! box * ** Now do the lower GI pathways: * ** **************************** * ** The lung(b)-->stomach-->small intestine-->uli-->lli route: lambda_bio(1) = lambda_lung(2,lng_class) lambda_bio(2) = lambda_gi(1) lambda_bio(3) = lambda_gi(2) + lambda_b lambda_bio(4) = lambda_gi(3) lambda_bio(5) = lambda_gi(4) transfer(1) = D_np * fraction_lung(2,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 - f1 transfer(5) = 1.0 call CHAIN (5,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,6,t,1) = V(j,4) + U(j,6,t,1) U(j,3,t,1) = V(j,5) + U(j,3,t,1) end do * ** The lung(d)-->stomach-->small intestine-->uli-->lli route: lambda_bio(1) = lambda_lung(4,lng_class) lambda_bio(2) = lambda_gi(1) lambda_bio(3) = lambda_gi(2) + lambda_b lambda_bio(4) = lambda_gi(3) lambda_bio(5) = lambda_gi(4) transfer(1) = D_tb * fraction_lung(4,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 - f1 transfer(5) = 1.0 call CHAIN (5,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,6,t,1) = V(j,4) + U(j,6,t,1) U(j,3,t,1) = V(j,5) + U(j,3,t,1) end do * ** The lung(f)-->lung(d)-->stomach-->SI-->ULI-->LLI route: lambda_bio(1) = lambda_lung(6,lng_class) lambda_bio(2) = lambda_lung(4,lng_class) lambda_bio(3) = lambda_gi(1) lambda_bio(4) = lambda_gi(2) + lambda_b lambda_bio(5) = lambda_gi(3) lambda_bio(6) = lambda_gi(4) transfer(1) = D_p * fraction_lung(6,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 transfer(5) = 1.0 - f1 transfer(6) = 1.0 call CHAIN (6,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,6,t,1) = V(j,5) + U(j,6,t,1) U(j,3,t,1) = V(j,6) + U(j,3,t,1) end do * ** The lung(g)-->lung(d)-->stomach-->SI-->ULI-->LLI route: lambda_bio(1) = lambda_lung(7,lng_class) lambda_bio(2) = lambda_lung(4,lng_class) lambda_bio(3) = lambda_gi(1) lambda_bio(4) = lambda_gi(2) + lambda_b lambda_bio(5) = lambda_gi(3) lambda_bio(6) = lambda_gi(4) transfer(1) = D_p * fraction_lung(7,lng_class) transfer(2) = 1.0 transfer(3) = 1.0 transfer(4) = 1.0 transfer(5) = 1.0 - f1 transfer(6) = 1.0 call CHAIN (6,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,6,t,1) = V(j,5) + U(j,6,t,1) U(j,3,t,1) = V(j,6) + U(j,3,t,1) end do * ** The stomach-->small intestine-->uli-->lli ingestion route: lambda_bio(1) = lambda_gi(1) lambda_bio(2) = lambda_gi(2) + lambda_b lambda_bio(3) = lambda_gi(3) lambda_bio(4) = lambda_gi(4) transfer(1) = 1.0 transfer(2) = 1.0 transfer(3) = 1.0 - f1 transfer(4) = 1.0 call CHAIN (4,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,6,t,2) = V(j,3) + U(j,6,t,2) U(j,3,t,2) = V(j,4) + U(j,3,t,2) end do * ** The lung(h)-->lung(j) route: lambda_bio(1) = lambda_lung( 8,lng_class) lambda_bio(2) = lambda_lung(10,lng_class) transfer(1) = D_p * fraction_lung( 8,lng_class) transfer(2) = fraction_lung(10,lng_class) call CHAIN (2,lambda_bio,ndaughters,lambda_rad, + transfer,time(t),V) do j = 0, ndaughters U(j,9,t,1) = V(j,2) + U(j,9,t,1) end do * ** Distribute the transformations in the "whole body" or in * ** "all other" among all the tissues not explicitly listed in * ** the metabolic model: do i = 1, nintakes do s = 1, 19 if (unused(s)) then do j = 0, ndaughters if (s .ne. 19) then U(j,s,t,i) = U(j,s,t,i) + U(j,19,t,i) * + mass(s) / unused_mass else U(j,s,t,i) = U(j,19,t,i) * + mass(s) / unused_mass end if ! 19 end do ! daughters end if ! unused end do ! s end do ! intakes * ** Add the transformations in the transfer compartment to those * ** in the whole body: do i = 1, nintakes do j = 0, ndaughters U(j,19,t,i) = U(j,19,t,i) + u_tc(j,i) end do end do end do ! time return end ! u_calc ************************************************************************ * * SUBROUTINE EDE_CALC (ntargets, ntimes, nintakes, + dcf) * * ************************************************************************ * * * Subroutine EDE_CALC calculations the whole-body dose and the * * effective dose equivalent from the weighted sum of the dcfs for * * individual organs. * * * ************************************************************************ integer target, ntargets, time, ntimes, intake, nintakes, i, j real dcf(22,ntimes,nintakes), dummy, array(11), mass(19) data (mass(target), target = 1, 19) + / 14.0, 1500.0, 160.0, 640.0, 150.0, 210.0, 310.0, + 1800.0, 1000.0, 28000.0, 11.0, 100.0, 1500.0, 120.0, + 180.0, 35.0, 20.0, 45.0, 70000.0 / do intake = 1, nintakes do time = 1, ntimes * ** Calculate the whole-body dcf by using the mass-weighted sum * ** of the dcfs for all organs: do target = 1, ntargets dcf(21,time,intake) = dcf(target,time,intake) * + mass(target) / mass(19) + dcf(21,time,intake) end do ! target * ** Next, calculate the "effective dose equivalent": if (dcf(11,time,intake) .gt. dcf(16,time,intake)) then dcf(22,time,intake) = 0.25 * dcf(11,time,intake) else dcf(22,time,intake) = 0.25 * dcf(16,time,intake) end if dcf(22,time,intake) = 0.15 * dcf(10,time,intake) + + 0.12 * dcf(13,time,intake) + + 0.12 * dcf( 9,time,intake) + + 0.03 * dcf(17,time,intake) + + 0.03 * dcf(14,time,intake) + + dcf(22,time,intake) do i = 1, 8 array(i) = dcf(i,time,intake) end do array( 9) = dcf(12,time,intake) array(10) = dcf(15,time,intake) array(11) = dcf(18,time,intake) do i = 1, 10 do j = (i+1), 11 if (array(i) .le. array(j)) then dummy = array(i) array(i) = array(j) array(j) = dummy end if end do end do do i = 1, 5 dcf(22,time,intake) = 0.06 * array(i) + + dcf(22,time,intake) end do end do ! time end do ! intake return end ! ede_calc ************************************************************************ * * SUBROUTINE OUTPUT (out, nuclide_name, zai, decay_c, nbranches, + branchfrc, ndaughter, daughter_zais, + time, dcf, ntimes, nintakes, + target_name, nuclide) * * ************************************************************************ * * * Subroutine OUTPUT prints the results in a text file and in a data * * file for later use. * * * * Input * * ----- * * out: logical units associated with the output files; * * nuclide_name: name of the radionuclide; * * zai: nuclide identifier; * * time: the array of integration times; * * dcf: array of dose conversion factors. * * ntargets, ntimes, nintakes: number of target organs, integration * * times, and intake modes. * * * ************************************************************************ integer out(6), zai(0:14), ntimes, nintakes, t, nuclide, s integer u, v, w, x, nbranches, ndaughter(3), daughter_zais(3,5) character*4 time(ntimes) character*8 nuclide_name, target_name(22) character*9 dy_mo_yr real dcf(22,ntimes,nintakes), branchfrc(3), decay_c ****** Write text output file: if (mod((nuclide-1),5) .eq. 0) then call DATE (dy_mo_yr) write(out(1),2000) dy_mo_yr, 'INHALED', + (target_name(s), s = 1, 11) write(out(2),2000) dy_mo_yr, 'INHALED', + (target_name(s), s = 12, 22) write(out(3),2000) dy_mo_yr, 'INGESTED', + (target_name(s), s = 1, 11) write(out(4),2000) dy_mo_yr, 'INGESTED', + (target_name(s), s = 12, 22) end if write(out(1),1000) nuclide_name, (time(t), (dcf(s,t,1), + s = 1, 11), t = 1, ntimes) write(out(2),1000) nuclide_name, (time(t), (dcf(s,t,1), + s = 12, 22), t = 1, ntimes) write(out(3),1000) nuclide_name, (time(t), (dcf(s,t,2), + s = 1, 11), t = 1, ntimes) write(out(4),1000) nuclide_name, (time(t), (dcf(s,t,2), + s = 12, 22), t = 1, ntimes) ****** Write data output file: write(out(5),3000) nuclide_name, zai(0), decay_c, nbranches, + (branchfrc(u), u = 1,3), (ndaughter(v), v = 1,3), + ((daughter_zais(w,x), x = 1,5) , w = 1,3), + ((dcf(s,t,1), s = 3, 18), + (dcf(s,t,1), s = 21, 22), + t = 1, ntimes) write(out(6),3000) nuclide_name, zai(0), decay_c, nbranches, + (branchfrc(u), u = 1,3), (ndaughter(v), v = 1,3), + ((daughter_zais(w,x), x = 1,5) , w = 1,3), + ((dcf(s,t,2), s = 3, 18), + (dcf(s,t,2), s = 21, 22), + t = 1, ntimes) return ****** NOTE: the "8" in the following format statement is (ntimes - 1) 1000 format (3X, A8, 3X, A4, 3X, 11(2X,1PE8.2), + 8(/3X, 8X, 3X, A4, 3X, 11(2X,1PE8.2)) + / 3X, 128('-')) 2000 format (2X, A9, X, 8X, 21X, + 'INTERNAL DOSE CONVERSION FACTORS: REM/CI ', A8 + // 3X, 72X, 'Organ' + / 3X, 20X, 108('-') + / 3X, 'Nuclide Time ', 11(2X, A8) + / 3X, 128('=')) 3000 format (A8, X, I6, 2X, 1PE9.3, 2X, I1, 3(X, F5.3), 3(X, I1), + 15(X, I6), 9(/ 18(X, 1PE8.2))) end ! output ************************** LEVEL 3 ROUTINES **************************** * * ************************************************************************ * * REAL FUNCTION AF (target, source, type, energy, + surface, Er, afp) * * ************************************************************************ * * * Function AF calculates the fraction of energy absorbed in the * * target organ per gram from radiations of a given type in the * * source organ. * * * * Input * * ----- * * target, source: the target and source organ identifiers; * * type: the particle type (a=alpha, r=recoil, e=electron, p=photon); * * energy: the particle energy (MeV); * * surface: is the radionuclide deposited primarily on bone surfaces? * * Er, afp: reference energies and absorbed fractions for photons * * emitted in source organ and absorbed in target organ; * * * ************************************************************************ logical surface, GI, bone, marrow, bladder integer target, source character type*1 real energy, Er(12), afp(12), mass(19), path(19), INTERPOLATE, mu data (mass(t), t = 1, 19) + / 14.0, 1500.0, 160.0, 640.0, 150.0, 210.0, 310.0, + 1800.0, 1000.0, 28000.0, 11.0, 100.0, 1500.0, 120.0, + 180.0, 35.0, 20.0, 45.0, 70000.0 / data (path(t), t = 1, 19) + / 0.25, 1.0, 0.0, 0.0, 0.0, 0.0, 0.67, + 0.51, 0.49, 0.16, 0.38, 0.66, 0.021, 0.0, + 0.56, 0.34, 0.35, 0.0, 0.49 / * ** Figure out if any of three special conditions hold: * ** GI <-- GI, bone <-- bone, marrow <-- bone GI = .false. bone = .false. marrow = .false. bladder = .false. if ((source .ge. 3) .and. (source .le. 6) .and. + (target .eq. source)) then GI = .true. else if ((source .eq. 14) .and. (target .eq. 14)) then bone = .true. else if ((source .eq. 14) .and. (target .eq. 13)) then marrow = .true. else if ((source .eq. 18) .and. (target .eq. 18)) then bladder = .true. end if if (type .eq. 'a') then ! ALPHAS ! ****** if ((GI) .or. (bladder)) then AF = 0.005 else if (bone) then if (surface) then AF = 0.25 ! 0.5*0.25 + 0.5*0.25 else AF = 0.013 ! 0.2*0.025 + 0.8*0.01 end if else if (marrow) then if (surface) then AF = 0.25 ! 0.5*0.5 + 0.5*0.0 else AF = 0.010 ! 0.2*0.05 + 0.8*0.0 end if else if (target .eq. source) then AF = 1.0 else AF = 0.0 end if AF = AF / mass(target) else if (type .eq. 'r') then ! RECOILS ! ******* if ((GI) .or. (bone) .or. (marrow) .or. (bladder)) then AF = 0.0 else if (target .eq. source) then AF = 1.0 else AF = 0.0 end if AF = AF / mass(target) else if (type .eq. 'e') then ! ELECTRONS ! ********* if ((GI) .or. (bladder)) then AF = 0.50 else if (bone) then if (surface) then if (energy .lt. 0.2) then AF = 0.25 ! 0.5*0.25 + 0.5*0.25 else AF = 0.020 ! 0.5*0.025 + 0.5*0.015 end if else AF = 0.017 ! 0.2*0.025 + 0.8*0.015 end if else if (marrow) then if (surface) then AF = 0.25 ! 0.5*0.5 + 0.5*0.0 else AF = 0.070 ! 0.2*0.35 + 0.8*0.0 end if else if (target .eq. source) then AF = 1.0 else AF = 0.0 end if AF = AF / mass(target) else if (type .eq. 'p') then ! PHOTONS ! ******* if (energy .lt. 0.01) then mu = 1.06 * 2.027E-06 * energy**(-3.178) if ((GI) .or. (bladder)) then AF = exp(INTERPOLATE (log(energy), Er, afp, 1, 12)) else if (bone) then if (surface) then AF = (1.0 - exp(-0.001*mu)) / mass(target) else AF = exp(INTERPOLATE(log(energy), Er, afp,1,12)) end if else if (marrow) then if (surface) then AF = 0.25 * exp(-0.001*mu) / mass(target) else AF = exp(INTERPOLATE(log(energy), Er, afp,1,12)) end if else if (target .eq. source) then AF = (1.0 - exp(-mu*path(target))) / mass(target) else AF = 0.0 end if else AF = exp(INTERPOLATE (log(energy), Er, afp, 2, 12)) end if end if ! type return end ! AF ************************************************************************ * * SUBROUTINE CHAIN (norgans, lambda_bio, + ndaughters, lambda_rad, transfer, time, V) * * ************************************************************************ * * * Subroutine CHAIN solves for the integral of the activity of each * * daughter in each compartment of a linear chain. * * * * Input: * * ----- * * norgans: the number of organ compartments in the linear chain; * * lambda_bio: the biological decay constants of each compartment; * * ndaughters: the number of radioactive daughters; * * lambda_rad: the radioactive decay constants of each daughter; * * transfer: the transfer fractions to each compartment; * * time: the integration time. * * * * Ouput: * * ----- * * V: the results for each daughter in each compartment. * * * ************************************************************************ logical simple integer norgans, ndaughters, i, j, k, s real lambda_bio(6), lambda_rad(0:ndaughters), + transfer(norgans), time, V(0:14,6), box_transfer, + daughter_transfer, EXPINT, PRODUCT, PRODUCT_MINUS, sum real*8 DPRODUCT, DPRODUCT_MINUS, dsum * ** If (lambda_bio(i) + lambda_rad(k))*time >> 1 for all i and k, * ** then the calculation is greatly simplified: simple = .true. do s = 1, norgans do j = 0, ndaughters if (((lambda_bio(s)+lambda_rad(j))*time) .lt. 13.) then simple = .false. end if end do end do if (simple) then do s = 1, norgans do j = 0, ndaughters if (s .eq. 1) then if (j .eq. 0) then box_transfer = transfer(1) else box_transfer = 0.0 end if else box_transfer = transfer(s) * V(j,(s-1)) * + lambda_bio(s-1) end if if (j .eq. 0) then daughter_transfer = 0.0 else daughter_transfer = V((j-1),s) * lambda_rad(j) end if V(j,s) = (box_transfer + daughter_transfer) / + (lambda_bio(s) + lambda_rad(j)) end do end do else * ** Remove redundancies in lambda arrays: do s = 1, (norgans-1) do i = (s+1), norgans if (lambda_bio(s) .eq. lambda_bio(i)) then lambda_bio(i) = (1.000 + 0.010*(0.5 - RRAND()))* + lambda_bio(i) end if end do end do do j = 0, (ndaughters-1) do i = (j+1), ndaughters if (lambda_rad(j) .eq. lambda_rad(i)) then lambda_rad(i) = (1.000 + 0.010*(0.5 - RRAND()))* + lambda_rad(i) end if end do end do * ** Now do the main calculation: do s = 1, norgans do j = 0, ndaughters if (PRODUCT(lambda_rad,0,ndaughters,0,j) .gt. 0.0) then sum = 0.0 do i = 1, s do k = 0, j sum = EXPINT((lambda_bio(i)+ + lambda_rad(k)),time) / + (PRODUCT_MINUS(lambda_bio,1,s,i) * 1.0E+30 * + PRODUCT_MINUS(lambda_rad,0,j,k)) + sum end do ! k end do ! i V(j,s) = PRODUCT(transfer,1,norgans,1,s) * + PRODUCT(lambda_bio,1,norgans,1,(s-1))*1.E30* + PRODUCT(lambda_rad,0,ndaughters,1,j) * sum else dsum = 0.0D+00 do i = 1, s do k = 0, j dsum = DBLE(EXPINT((lambda_bio(i) + + lambda_rad(k)),time)) / + (DPRODUCT_MINUS(lambda_bio,1,s,i) * + DPRODUCT_MINUS(lambda_rad,0,j,k)) + dsum end do ! k end do ! i V(j,s) = REAL(DPRODUCT(transfer,1,norgans,1,s) * + DPRODUCT(lambda_bio,1,norgans,1,(s-1)) * + DPRODUCT(lambda_rad,0,ndaughters,1,j) * dsum) end if end do ! j end do ! s end if ! simple return end ! chain ************************** LEVEL 4 ROUTINES **************************** * * ************************************************************************ * * REAL FUNCTION INTERPOLATE (x, xr, yr, degree, + npts) * * ************************************************************************ * * * Function INTERPOLATE performs a Lagrangian interpolation of the * * dependent variable y for a given value of the independent * * variable x. Reference: James & James, Mathematics Dictionary, * * 3rd edition, Van Nostrand. * * * * Input * * ----- * * x: the value of the absicca for which the interpolation is * * desired. * * xr: the array of distinct, monotonically increasing abcissas. * * yr: the array of corresponding ordinates. * * degree: the degree of the interpolation polynomial; * * 1 = linear, 2 = quadratic, 3 = cubic, etc. * * npts: the number of elements in the xr and yr arrays. * * * ************************************************************************ integer degree, npts, i, j, k, firstpt, lastpt real x, xr(npts), yr(npts), y, sum, product * ** Locate the position of x in the array xr... i = 1 do while (((x - xr(i)) .gt. 0.0) .and. (i .lt. npts)) i = i + 1 end do * ** If x is equal to one of the points in xr, then no interpolation * ** is necessary... if ((x - xr(i)) .eq. 0.0) then y = yr(i) * ** Otherwise, use interpolation. First, make sure that the inputs * ** make sense... else if (degree .gt. (npts - 1)) then degree = npts - 1 else if (degree .eq. 0) then stop 'degree of polynomial=0' end if * ** Now calculate the first and last points of the polynomial. * ** If x < xr(2) or x > xr(npts-1), then the job is simple... if ((i .eq. 1) .or. (i .eq. 2)) then firstpt = 1 else if (i .eq. npts) then firstpt = npts - degree * ** If xr(2) < x < xr(npts-1), and x is closer to * ** xr(i-1) than xr(i), then prefer xr(i-2) over xr(i+1)... else if (x .lt. ((xr(i-1) + xr(i)) / 2.0)) then firstpt = i - (degree + 2) / 2 * ** On the other hand, if it's closer to xr(i), then... else firstpt = i - (degree + 1) / 2 end if if (firstpt .lt. 1) then firstpt = 1 end if lastpt = firstpt + degree if (lastpt .gt. npts) then lastpt = npts end if * ** Now let's interpolate... sum = 0.0 do k = firstpt, lastpt product = 1.0 do j = firstpt, lastpt if (k .ne. j) then product = product * (x - xr(j)) / + (xr(k) - xr(j)) end if end do sum = sum + yr(k) * product end do y = sum end if INTERPOLATE = y return end ! interpolate ************************************************************************ * * REAL FUNCTION EXPINT (decay_constant, + integration_time) * * ************************************************************************ * * * Function EXPINT calculates the integral of exp(-kt) between * * t = 0 and t = integration_time for k = decay_constant, allowing * * for the special case of very small (or zero) decay_constant. * * * ************************************************************************ real decay_constant, integration_time, product product = decay_constant * integration_time if (product .lt. 1.0E-04) then EXPINT = integration_time * (1.0 - product / 2.0) else EXPINT = (1.0 - EXP(-product)) / decay_constant end if return end ! expint ************************************************************************ * * REAL FUNCTION PRODUCT (lambda, l, u, lower, upper) * * ************************************************************************ * * * Function PRODUCT calculates the product of the elements of an * * array between the lower and upper subscripts. * * * ************************************************************************ integer l, u, lower, upper, i real lambda(l:u), prod prod = 1.0 do i = lower, upper prod = lambda(i) * prod end do PRODUCT = prod return end ************************************************************************ * * REAL FUNCTION PRODUCT_MINUS (lambda, lower, + upper, index) * * ************************************************************************ * * * Function PRODUCT calculates the product of all elements of an * * array minus a certain element. It checks to make sure that there * * are no redundancies in the lambda array, and removes any that * * are found. * * * ************************************************************************ integer lower, upper, i, index real lambda(lower:upper), prod prod = 1.0 do i = lower, upper if (i .ne. index) then prod = (lambda(i) - lambda(index)) * prod end if end do if (prod .eq. 0.0) then stop 'redundancy in lambda' else PRODUCT_MINUS = prod end if return end ************************************************************************ * * REAL*8 FUNCTION DPRODUCT (lambda, l, u, + lower, upper) * * ************************************************************************ * * * Function PRODUCT calculates the product of the elements of an * * array between the lower and upper subscripts. * * * ************************************************************************ integer l, u, lower, upper, i real*4 lambda(l:u) real*8 prod prod = 1.0D+00 do i = lower, upper prod = DBLE(lambda(i)) * prod end do DPRODUCT = prod return end ************************************************************************ * * REAL*8 FUNCTION DPRODUCT_MINUS (lambda, lower, + upper, index) * * ************************************************************************ * * * Function PRODUCT calculates the product of all elements of an * * array minus a certain element. It checks to make sure that there * * are no redundancies in the lambda array, and removes any that * * are found. * * * ************************************************************************ integer lower, upper, i, index real*4 lambda(lower:upper) real*8 prod prod = 1.0D+00 do i = lower, upper if (i .ne. index) then prod = DBLE((lambda(i) - lambda(index))) * prod end if end do if (prod .eq. 0.0D+00) then stop 'redundancy in lambda' else DPRODUCT_MINUS = prod end if return end