************************************************************************ * * PROGRAM EDRCF * * ************************************************************************ * * * * Steve Fetter * Version 1.0 for MSDOS * * Subcontractor for EG&G, Idaho * Date: 06/12/90 * * Idaho Falls, ID * Copyright (c) 1990 * * * * ************************************************************************ * * * Program EDRCF calculates photon dose-rate conversion factors for * * external exposure to a semi-infinite cloud, an infinite plane, * * and a semi-infinite slab of soil uniformly contaminated with an * * isotope. The clound and plane EDRCFs are calculated using * * the equations in D.C. Kocher, "Dose-Rate Conversion Factors for * * External Exposure to Photon and Electron Radiation from Radionu- * * clides Occuring in Routine Releases from Nuclear Fuel Cycle Facil- * * ities", ORNL/NUREG/TM-283 (1979). The photon DRCFs are * * estimated based on the calculations of H. Beck and G. de Planque, * * "The Radiation Field in Air Due to Distributed Gamma-Ray Sources * * in the Ground," HASL-195 (U.S. AEC, May 1968). * * * * To facilitate waste disposal calculations, the shielding factor * * provided by 1 meter of dry earth has been calculated. This is * * estimated here by a slab of 100% SiO2 at 1.6 g/cm**3. The * * buildup factor is estimated by the factor (1.0 + kur), where * * k = (u - uen)/uen. Although this procedure is at least as good * * as that used in 10 CFR 61, this formula is known to underestimate * * the buildup factor somewhat for large values of ur. For the * * gamma emission of Cs-137, the result is 0.70 that given by 10 CFR * * 61, so in order to be conservative, the buildup factor given by * * the above equation is multiplied by 1.50. * * * * 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." * * * * This version compiled with Lahey Fortran and these options: * * * * /0/N2/N3/N7/A/A1/B/NC/ND/E/NF/NH/NI/NK/L/NO/P/R/NS/NT/NU/W/NX/Z1 * * * * Input files * * ----- ----- * * decaylib.dat: decay-scheme data file; * * * * Output files * * ------ ----- * * edrcf.1,2,3,4,5,6,7: the DRCF text output files; * * edrcf.dat: the DRCF data output file. * * * ************************************************************************ logical exist integer nuclide, nnuclides, dk, key, con, dat, out(7), zai, + t1, t2, org, i character*8 organ_name(22), nuclide_name, dy_mo_yr, hr_min_sec, + dk_version character char*1, decaylib*12, outfile*12(8) real drcf(22,3), shield_factor, deltat parameter (dk = 1, key = 5, con = 6, dat = 14) data (out(i), i = 1, 7) / 7, 8, 9, 10, 11, 12, 13 / data (organ_name(org), org = 1, 22) + / ' Skin ', ' Thymus ', 'GI: LLI ', ' GI: SI ', ' GI: ST ', + 'GI: ULI ', 'Kidneys ', ' Liver ', ' Lungs ', ' Muscle ', + 'Ovaries ', 'Pancreas', ' Marrow ', 'Skeleton', ' Spleen ', + ' Testes ', 'Thyroid ', 'Bladder ', ' Uterus ', 'Y.Marrow', + 'Tot.Body', ' EDE ' / data (outfile(i), i = 1, 8) + / 'EDRCF.1', 'EDRCF.2', 'EDRCF.3', 'EDRCF.4', 'EDRCF.5', + 'EDRCF.6', 'EDRCF.7', 'EDRCF.DAT' / call ovefl('true') call under0('true') ****** Print sign on message: write (con,1000) ****** Start timer... call TIMER(t1) call DATE (dy_mo_yr) call TIME (hr_min_sec) ****** Open the data and output files... write (con,1001) read (key,1002) decaylib inquire (file = decaylib, exist = exist) if (exist) then open (unit = dk, file = decaylib, status = 'old', + action = 'read') else write (con,1003) decaylib stop 'Sorry' end if exist = .false. do i = 1, 8 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 edrcf.?') call SYSTEM ('del edrcf.dat') else stop 'rename files and restart program' end if end if do i = 1, 7 open (unit = out(i), file = outfile(i), status = 'new') end do open (unit = dat, file = 'EDRCF.DAT', status = 'new') ****** Begin main loop... read (dk,2000) dk_version char = '*' do while (char .eq. '*') read (dk,2001) char end do write (dat,4000) dy_mo_yr, hr_min_sec, dk_version, + (organ_name(org), org = 3, 18), + (organ_name(org), org = 21, 22) write (con,1008) read (dk, 2002) nnuclides write (dat,4001) nnuclides do nuclide = 1, nnuclides call DRCF_CALC (dk, nuclide, nuclide_name, zai, drcf, + shield_factor) write (con,1009) nuclide_name call EDE_CALC (drcf) call OUTPUT (out, dat, organ_name, nuclide, nuclide_name, + zai, drcf, shield_factor) end do ! nuclide call TIMER(t2) deltat = 0.01 * REAL(t2 - t1) / 60.0 write (con,1010) deltat close (unit = dk, status = 'keep') do i = 1, 7 close (unit = out(i), status = 'keep') end do close (unit = dat, status = 'keep') stop 'adios' 1000 format (//16X,'*************************************************' + /16X,'* EDRCF *' + /16X,'* A Program by Steve Fetter to Calculate *' + /16X,'* External Dose-Rate Conversion Factors *' + /16X,'*************************************************') 1001 format ( /16X, 'What is the name of the radionuclide file? ') 1002 format (A12) 1003 format ( /16X, 'Error: file ', A12, ' doesn''t exist!') 1004 format ( /16X, 'EDRCF output files already exist.' + /16X, ' Do you want to erase them (Y/N)? ') 1005 format (A1) 1008 format ( /16X, 'Calculating Dose-Rate Conversion Factors...') 1009 format ('+', A8) 1010 format (// 16X, 'Job ran in', F4.1, ' minutes.' + / 16X, 'The results are in EDRCF.1,2,3,4,5,6,7 ' + 'and EDRCF.DAT.' /) 2000 format (//63X, A5) 2001 format (A1) 2002 format (3X, 3X, I3) 4000 format (A8, X, A8, ' DECAYLIB version ', A8 / 162('*') + / '* This data file contains external dose-rate ' + 'conversion factors (EDRCFs) for exposure to an ' + 'infinite, uniformly contaminated cloud, plane, ' + 'and slab, for use by *' + / '* the DOSE code. EDRCFs are given for the ' + 'following 18 target organs:', 91(' '), '*' + / '*', X, A7, 16(X, A8), X, A7, '*' / 162('*')) 4001 format (I3, ' = number of radionuclides in file') end ! edrcf ************************** LEVEL 2 ROUTINES ************************** * * ************************************************************************ * * BLOCK DATA DRCF_DATA * * ************************************************************************ integer i, org real er(25), muar(25), muear(25), muetr(25), car(25), dar(25), + musr(25), csr(25), eAFr(12), AFr(12,21), + evolr(15), rvolr(15) common /drcf_data/ er, muar, muear, muetr, car, dar, + musr, csr, eAFr, AFr, evolr, rvolr ****** Photon attenuation data taken from Kocher, ORNL/NUREG/TM-283, ****** and the Radiological Health Handbook, HEW. data (er(i), i = 1, 25) + /0.010,0.015,0.020,0.030,0.040,0.050,0.060, +0.080,0.100,0.150,0.200,0.300,0.400,0.500,0.600, + 0.800,1.000,1.500,2.000,3.000,4.000,5.000,6.000, + 8.000, 10.00 / ****** Photon mass attenuation coefficient in air (cm2/g): data (muar(i), i = 1, 25) + / 4.99,1.55,0.752,0.349,0.248,0.208,0.188, + 0.167,0.154,0.136,0.123,0.107,0.0954,0.0870,0.0805, + 0.0707,0.0636, 0.0518, 0.0445, 0.0358, 0.0308, 0.0275, 0.0252, + 0.0223, 0.0204 / ****** Photon mass energy-absorption coefficient in air (cm2/g): data (muear(i), i = 1, 25) + / 4.61, 1.27, 0.511, 0.148, 0.0668, 0.0406, 0.0305, + 0.0243, 0.0234, 0.0250, 0.0268, 0.0287, 0.0295, 0.0296, 0.0295, + 0.0289, 0.0278, 0.0254, 0.0234, 0.0205, 0.0186, 0.0174, 0.0164, + 0.0152, 0.0145 / ****** Photon mass energy-absorption coefficient in tissue (cm2/g): data (muetr(i), i = 1, 25) + / 4.87,1.32,0.533,0.154,0.0701, 0.0431, 0.0328, + 0.0264, 0.0256, 0.0275, 0.0294, 0.0317, 0.0325, 0.0328, 0.0325, + 0.0318, 0.0308, 0.0282, 0.0259, 0.0226, 0.0203, 0.0188, 0.0178, + 0.0163, 0.0154 / ****** Photon mass attenuation coefficient in SiO2 (cm2/g): data (musr(i), i = 1, 25) + / 19.0, 5.73, 2.49, 0.859, 0.463, 0.318, 0.252, + 0.194, 0.169, 0.140, 0.126, 0.108, 0.0959, 0.0874, 0.0808, + 0.0707, 0.0636, 0.0518, 0.0447, 0.0363, 0.0317, 0.0287, 0.0266, + 0.0241, 0.0226 / ****** These C coefficients were calculated by hand for SiO2: data (csr(i), i = 1, 25) + / 0.0306, 0.0762, 0.142, 0.384, 0.783, 1.34, 1.98, + 3.22, 3.96, 4.01, 3.51, 2.71, 2.23, 1.93, 1.73, + 1.45, 1.29, 1.04, 0.910, 0.746, 0.653, 0.589, 0.533, + 0.458, 0.407 / ****** Berger C and D coefficients for buildup factors in air: data (car(i), i = 1, 25) + / 0.0330,0.132,0.330,1.12,2.31,3.46,4.20, +4.46,4.15,3.25,2.70,2.19,1.87,1.68,1.53, + 1.34,1.21,0.995,0.838,0.687,0.602,0.546,0.508, + 0.457,0.426 / data (dar(i), i = 1, 25) + / -0.0680,-0.0693,-0.0607,-0.0230, 0.0284,0.0743, 0.108, + 0.151,0.168,0.178,0.172,0.142,0.125,0.106,0.0946, + 0.0745,0.0604,0.0373,0.0242,0.0100, 0.0032,-0.0012,-0.0029, + -0.0035,-0.0033 / ****** Ratio of organ dose rate to body-surface dose rate given ****** by Poston and Snyder, Health Physics, Vol. 26, p. 287: data (eAFr(i), i = 1, 12) + / 0.010, 0.015, 0.020, 0.030, 0.050, 0.100, + 0.200, 0.500, 1.000, 1.500, 2.000, 4.000 / data (AFr(i,1), i = 1, 12) + / 0.170, 0.310, 0.391, 0.514, 0.659, 0.762, + 0.803, 0.755, 0.815, 0.757, 0.830, 0.800 / data (AFr(i,2), i = 1, 12) + / 1.10e-07, 1.78e-04, 5.80e-03, 0.119, 0.325, 0.340, + 0.544, 0.385, 0.424, 0.640, 0.716, 1.23 / data (AFr(i,3), i = 1, 12) + / 3.08e-04, 9.10e-03, 8.59e-03, 2.76e-02, 0.201, 0.346, + 0.372, 0.463, 0.399, 0.601, 0.578, 0.643 / data (AFr(i,4), i = 1, 12) + / 1.45e-06, 1.13e-04, 1.79e-03, 3.95e-02, 0.221, 0.426, + 0.465, 0.461, 0.483, 0.534, 0.548, 0.693 / data (AFr(i,5), i = 1, 12) + / 1.87e-03, 1.79e-03, 1.04e-02, 8.99e-02, 0.317, 0.426, + 0.437, 0.628, 0.529, 0.648, 0.624, 0.480 / data (AFr(i,6), i = 1, 12) + / 2.65e-06, 2.84e-05, 3.46e-03, 5.59e-02, 0.245, 0.405, + 0.488, 0.513, 0.563, 0.712, 0.830, 0.800 / data (AFr(i,7), i = 1, 12) + / 1.08e-06, 3.28e-04, 2.25e-02, 0.127, 0.317, 0.477, + 0.486, 0.547, 0.592, 0.489, 0.525, 0.585 / data (AFr(i,8), i = 1, 12) + / 7.93e-06, 3.58e-04, 7.70e-03, 8.96e-02, 0.298, 0.498, + 0.514, 0.545, 0.529, 0.584, 0.616, 0.693 / data (AFr(i,9), i = 1, 12) + / 7.44e-04, 3.85e-03, 1.13e-02, 0.112, 0.373, 0.562, + 0.597, 0.590, 0.588, 0.659, 0.660, 0.691 / data (AFr(i,10), i = 1, 12) + / 3.08e-03, 2.54e-02, 7.37e-02, 0.216, 0.437, 0.574, + 0.605, 0.631, 0.626, 0.696, 0.683, 0.742 / data (AFr(i,11), i = 1, 12) + / 3.97e-04, 1.79e-03, 1.13e-02, 8.24e-02, 0.305, 0.276, + 0.365, 0.231, 0.379, 0.645, 0.729, 0.749 / data (Afr(i,12), i = 1, 12) + / 5.33e-04, 2.83e-03, 1.67e-02, 2.83e-02, 0.221, 0.339, + 0.395, 0.345, 0.634, 0.598, 0.714, 1.33 / data (AFr(i,13), i = 1, 12) + / 1.45e-03, 2.10e-02, 8.24e-02, 0.374, 1.04, 1.19, + 0.998, 0.772, 0.702, 0.712, 0.687, 0.695 / data (AFr(i,14), i = 1, 12) + / 1.36e-03, 1.93e-02, 7.70e-02, 0.343, 0.974, 1.19, + 0.994, 0.766, 0.693, 0.721, 0.704, 0.747 / data (AFr(i,15), i = 1, 12) + / 1.87e-04, 1.79e-03, 6.70e-03, 7.48e-02, 0.289, 0.405, + 0.541, 0.672, 0.622, 0.788, 0.431, 0.690 / data (AFr(i,16), i = 1, 12) + / 6.61e-03, 2.24e-02, 0.112, 0.311, 0.439, 0.617, + 0.866, 0.713, 0.567, 0.517, 0.347, 0.448 / data (AFr(i,17), i = 1, 12) + / 3.74e-04, 1.49e-03, 2.88e-02, 0.237, 0.524, 0.817, + 0.641, 0.525, 0.458, 0.701, 0.637, 1.33 / data (AFr(i,18), i = 1, 12) + / 1.30e-10, 3.73e-05, 2.68e-03, 7.48e-02, 0.232, 0.464, + 0.420, 0.502, 0.563, 0.640, 0.483, 1.33 / data (AFr(i,19), i = 1, 12) + / 7.40e-08, 2.54e-06, 7.37e-04, 2.91e-02, 8.32e-02, 0.276, + 0.399, 0.423, 0.437, 0.469, 0.544, 0.587 / data (AFr(i,20), i = 1, 12) + / 1.70e-03, 2.39e-02, 9.49e-02, 0.432, 1.16, 1.28, + 1.08, 0.808, 0.731, 0.743, 0.723, 0.746 / data (AFr(i,21), i = 1, 12) + / 9.03e-03, 3.28e-02, 7.81e-02, 0.222, 0.528, 0.660, + 0.669, 0.632, 0.622, 0.696, 0.683, 0.746 / ****** Dose rate 1 m above a uniformly contaminated infinite slab of ****** earth, from Table 2 of HASL-195 (MeV/g-s)/(dis/cm3-s): data (evolr(i), i = 1, 15) + / 0.0, 0.25, 0.364, 0.50, 0.662, 0.75, 1.00, + 1.25, 1.46, 1.76, 2.00, 2.25, 2.50, 2.62, 2.75 / data (rvolr(i), i = 1, 15) + / 0.0, 0.0599, 0.0968, 0.137, 0.186, 0.214, 0.288, + 0.364, 0.428, 0.518, 0.590, 0.659, 0.732, 0.769, 0.806 / end ! block data ************************************************************************ * * SUBROUTINE DRCF_CALC (dk, nuclide, nuclide_name, + zai, drcf, shield_factor) * * ************************************************************************ * * * Subroutine DRCF_CALC computes the dose factors for external * * radiation from photons. * * * * Input * * ----- * * dk: the logical unit associated with DECAYLIB; * * * * Output * * ------ * * nuclide_name: the name of the parent radionuclide; * * zai: radionuclide identifier number; * * drcf: array of DRCFs for cloud, plane, and soil contamination; * * shield_factor: array of shielding factors. * * * ************************************************************************ integer nalphas, nbetas, nelectrons, ngammas, org, i, nuclide, + dk, zai, z, j, k, daughter, ndaughters, order character nuclide_name*8 real er(25), muar(25), muear(25), muetr(25), car(25), dar(25), + musr(25), csr(25), eAFr(12), leAFr(12), AFr(12,21), + lAFr(12,21), evolr(15), rvolr(15), e, f, mua, muea, muet, + ca, da, mus, must, cs, ds, dummy(12), AF(21), rvol, + INTERPOLATE, ratio, branch, branch_sum, density_of_air, + height, thickness, density_of_dirt, drcf_sum(21,3), + sum_shield_factor, drcf(22,3), shield_factor, + erg_per_gxrad, dis_per_Cixs, erg_per_MeV, + disxgxrad_per_CixsxMeV, cm3_per_m3, cm2_per_m2, E1 common /drcf_data/ er, muar, muear, muetr, car, dar, + musr, csr, eAFr, AFr, evolr, rvolr ****** Density of air, height above the ground, and particle-medium ****** correction factors are taken from Kocher; other constants ****** taken from "Table of Isotopes"... parameter (density_of_air = 1.189E-03, height = 100.0) parameter (density_of_dirt = 1.60, thickness = 100.0) parameter (erg_per_gxrad = 100.0) parameter (dis_per_Cixs = 3.700E+10) parameter (erg_per_MeV = 1.6021892E-06) parameter (disxgxrad_per_CixsxMeV = + dis_per_Cixs * erg_per_MeV / erg_per_gxrad) parameter (cm3_per_m3 = 100.0**3, cm2_per_m2 = 100.0**2) if (nuclide .eq. 1) then ****** Divide dose-rate by energy for more accurate interpolation: do i = 2, 15 rvolr(i) = rvolr(i) / evolr(i) end do ****** Log-log interpolation is more accurate in some cases... do i = 1, 12 leAFr(i) = LOG(eAFr(i)) end do do org = 1, 21 do i = 1, 12 lAFr(i,org) = LOG(AFr(i,org)) end do end do end if ****** Set sums equal to zero... do i = 1, 3 do org = 1, 21 drcf_sum(org,i) = 0.0 end do end do sum_shield_factor = 0.0 ****** Start main loop... branch_sum = 0.0 do while (branch_sum .lt. 1.00) read (dk,1000) nuclide_name, branch, ndaughters branch_sum = branch + branch_sum read (dk,1001) zai z = zai / 10000 ****** Read the alpha, beta(-), beta(+), and electron data and ****** throw it away... read (dk,2000) nalphas do j = 1, nalphas read (dk,2001) e, f end do do k = 1, 2 read (dk,2000) nbetas do j = 1, nbetas read (dk,2001) e, f read (dk,2002) e end do ! betas end do ! k read (dk,2000) nelectrons do j = 1, nelectrons read (dk,2001) e, f end do ! electrons ****** Compute the gamma energy... read (dk,2000) ngammas do j = 1, ngammas read (dk,2001) e, f if (e .lt. 0.01) then mua = (3.022E-06 * e**(-3.109)) * density_of_air mus = (1.150E-05 * e**(-3.109)) * density_of_dirt muea = (2.027E-06 * e**(-3.178)) muet = 1.06 * muea ca = (e / 0.01) * car(1) cs = (e / 0.01) * csr(1) do org = 1, 21 AF(org) = (e / 0.01) * AFr(1,org) end do else mua = INTERPOLATE (e, er, muar, 2, 25) * + density_of_air mus = INTERPOLATE (e, er, musr, 2, 25) * + density_of_dirt muea = INTERPOLATE (e, er, muear, 2, 25) muet = INTERPOLATE (e, er, muetr, 2, 25) ca = INTERPOLATE (e, er, car, 2, 25) cs = INTERPOLATE (e, er, csr, 2, 25) do org = 1, 21 do i = 1, 12 dummy(i) = lAFr(i,org) end do AF(org) = EXP(INTERPOLATE (LOG(e), leAFr, + dummy, 1, 12)) end do end if ! e if (e .lt. 0.4) then order = 2 else order = 1 end if rvol = e * INTERPOLATE (e, evolr, rvolr, order, 15) ratio = muet / muea da = INTERPOLATE (e, er, dar, 2, 25) ds = da do org = 1, 21 drcf_sum(org,1) = e * f * AF(org) * branch * + ratio + drcf_sum(org,1) if (mua .lt. 0.50) then drcf_sum(org,2) = e * f * AF(org) * + branch * muet * + (E1 (DBLE(mua), DBLE(height)) - + (ca / (da - 1.0)) * + EXP((da-1.0) * mua * height)) + + drcf_sum(org,2) end if drcf_sum(org,3) = rvol * f * AF(org) * branch * + ratio + drcf_sum(org,3) end do must = mus * thickness sum_shield_factor = rvol * f * AF(21) * branch * ratio * + (1.0 + cs * must * EXP(ds*must)) * + EXP(-must) + sum_shield_factor end do ! gamma ****** Now gobble up the data for the radioactive daughters... do daughter = 1, ndaughters read (dk,1001) i read (dk,2000) nalphas do j = 1, nalphas read (dk,2001) e, f end do do k = 1, 2 read (dk,2000) nbetas do j = 1, nbetas read (dk,2001) e, f read (dk,2002) e end do end do do k = 1, 2 read (dk,2000) nelectrons do j = 1, nelectrons read (dk,2001) e, f end do end do end do ! daughter end do ! branch_sum ****** Now compute the DRCFs for contaminated air, surface and soil: do org = 1, 21 drcf(org,1) = disxgxrad_per_CixsxMeV * drcf_sum(org,1) / + (2.0 * cm3_per_m3 * density_of_air) if ((z .eq. 2) .or. (z .eq. 10) .or. (z .eq. 18) .or. + (z .eq. 36) .or. (z .eq. 54) .or. (z .eq. 86)) then drcf(org,2) = 0.0 else drcf(org,2) = disxgxrad_per_CixsxMeV * drcf_sum(org,2) / + (2.0 * cm2_per_m2) end if drcf(org,3) = disxgxrad_per_CixsxMeV * drcf_sum(org,3) / + cm3_per_m3 end do ! org ****** Finally, the shielding factor: if (drcf_sum(21,3) .eq. 0.0) then shield_factor = 0.0 else shield_factor = sum_shield_factor / drcf_sum(21,3) end if return 1000 format (/ A8, 2X, F6.4, 2X, I1) 1001 format (I7) 2000 format (I2) 2001 format (5X, F10.8, X, F11.9) 2002 format (5X, F9.7) end ! edrcf_calc ************************************************************************ * * SUBROUTINE EDE_CALC (drcf) * * ************************************************************************ * * * Subroutine EDE_CALC computes the effective dose equivalent. * * * * Input/Ouput * * ----------- * * drcf: the array of DRCFs for this nuclide. * * * ************************************************************************ integer i, j, k real drcf(22,3), array(12), dummy do i = 1, 3 ! for each exposure mode ****** the gonad dose is the greater of the ovaries and testes dose: if (drcf(11,i) .gt. drcf(16,i)) then drcf(22,i) = 0.25 * drcf(11,i) else drcf(22,i) = 0.25 * drcf(16,i) end if ****** now add breasts, red marrow, lung, thyroid, and bone surfaces: drcf(22,i) = 0.15 * drcf(10,i) + 0.12 * drcf(13,i) + + 0.12 * drcf( 9,i) + 0.03 * drcf(17,i) + + 0.03 * drcf(14,i) + drcf(22,i) ****** Finally, add the five other organs (except the skin) that have ****** the highest doses: do j = 1, 7 ! put the remaining array(j) = drcf((j+1),i) ! organs except skin end do ! into an array... array( 8) = drcf(12,i) array( 9) = drcf(15,i) array(10) = drcf(18,i) array(11) = drcf(19,i) array(12) = drcf(20,i) do j = 1, 11 ! and sort the array do k = (j+1), 12 if (array(j) .le. array(k)) then dummy = array(j) array(j) = array(k) array(k) = dummy end if end do end do do j = 1, 5 drcf(22,i) = 0.06 * array(j) + drcf(22,i) end do end do ! i return end ! ede_calc ************************************************************************ * * SUBROUTINE OUTPUT (out, dat, organ_name, nuclide, + nuclide_name, zai, drcf, shield_factor) * * ************************************************************************ * * * Subroutine OUTPUT prints the results in six different files. * * * * Input * * ----- * * out: array of logical units for output; * * dat: logical unit associated with the data output file; * * organ_name: the array of organ names; * * nuclide: the number of nuclides processed to this point; * * nuclide_name: the name of this nuclide; * * zai: radionuclide identifier number; * * drcf: the array of DRCFs for this nuclide; * * shield_factor: the shielding factor. * * * ************************************************************************ integer out(7), dat, nuclide, zai, org character*8 nuclide_name, organ_name(22) character*9 dy_mo_yr real drcf(22,3), shield_factor ****** Write text output file: ****** Every 48 data lines print a heading: if (mod((nuclide-1),48) .eq. 0) then call DATE (dy_mo_yr) write (out(1),3001) dy_mo_yr, 'CLOUD', '3', + (organ_name(org), org = 1, 11) write (out(2),3001) dy_mo_yr, 'CLOUD', '3', + (organ_name(org), org = 12, 22) write (out(3),3001) dy_mo_yr, 'PLANE', '2', + (organ_name(org), org = 1, 11) write (out(4),3001) dy_mo_yr, 'PLANE', '2', + (organ_name(org), org = 12, 22) write (out(5),3001) dy_mo_yr, 'SOIL ', '3', + (organ_name(org), org = 1, 11) write (out(6),3001) dy_mo_yr, 'SOIL ', '3', + (organ_name(org), org = 12, 22) write (out(7),3003) dy_mo_yr end if write (out(1),3000) nuclide_name, + (drcf(org,1), org = 1, 11) write (out(2),3000) nuclide_name, + (drcf(org,1), org = 12, 22) write (out(3),3000) nuclide_name, + (drcf(org,2), org = 1, 11) write (out(4),3000) nuclide_name, + (drcf(org,2), org = 12, 22) write (out(5),3000) nuclide_name, + (drcf(org,3), org = 1, 11) write (out(6),3000) nuclide_name, + (drcf(org,3), org = 12, 22) write (out(7),3002) nuclide_name, shield_factor ****** Write data output file: write (dat,4002) nuclide_name, zai, shield_factor, + (drcf(org,1), org = 3, 18), (drcf(org,1), org = 21, 22), + (drcf(org,2), org = 3, 18), (drcf(org,2), org = 21, 22), + (drcf(org,3), org = 3, 18), (drcf(org,3), org = 21, 22) return 3000 format ( X, 6X, A8, 2X, 11(2X, 1PE8.2)) 3001 format (//A9, 17X, 'EXTERNAL DOSE-RATE CONVERSION FACTORS', + ' FOR EXPOSURE TO A CONTAMINATED ', A5, + ' (REM/SEC)/(CI/M**', A1, ')' + // X, 6X, 'Nuclide ', 11(2X, A8) + / X, 6X, 120('=')) 3002 format ( X, 6X, A8, 2X, 2X, 1PE8.2) 3003 format (//A9, 10X, 'SHIELDING FACTOR FOR 1 M OF SOIL', + // X, 6X, 'Nuclide ', 2X, ' SF' + / X, 6X, 20('=')) 4002 format (A8, X, I6, 4X, 1PE8.2, 3(/ 18(X, 1PE8.2))) end ! output ************************** LEVEL 3 ROUTINES ************************** * * ************************************************************************ * * REAL FUNCTION E1 (mua, height) * * ************************************************************************ * * * Function E1 calculates the exponential integral using polynomial * * and rational approximations; see Abramowitz and Stegun, "Handbook * * of Mathematical Functions", Dover Press (1965), pg 231. This * * code is taken directly from Kocher. * * * ************************************************************************ double precision arg, al0, al1, al2, al3, al4, al5, ag1, ag2, ag3 double precision ag4, b1, b2, b3, b4, dnum, dnom, mua, height parameter (al0 = -0.57721566D0, al1 = 0.99999193D0, + al2 = -0.24991055D0, al3 = 0.05519968D0, + al4 = -0.00976004D0, al5 = 0.00107857D0) parameter (ag1 = 8.5733287401D0, ag2 = 18.059016973D0, + ag3 = 8.6347608925D0, ag4 = 0.2677737343D0) parameter (b1 = 9.5733223454D0, b2 = 25.6329561486D0, + b3 = 21.0996530827D0, b4 = 3.9584969228D0) arg = mua * height if (arg .le. 1.0D0) then E1 = al0 + al1 * arg + al2 * arg**2 + al3 * arg**3 + + al4 * arg**4 + al5 * arg**5 - DLOG(arg) else dnum = arg**4 + ag1 * arg**3 + ag2 * arg**2 + ag3 * arg +ag4 dnom = arg**4 + b1 * arg**3 + b2 * arg**2 + b3 * arg + b4 E1 = dnum / (dnom * arg * DEXP(arg)) end if return end ! E1 ************************************************************************ * * 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