*DECK WILLIE PROGRAM WILLIE C-Title : WILLIE Program C-Purpose: Operations on WIMS cross section library C-Author : A.Trkov, Inst."Jozef Stefan", Ljubljana, Slovenia, C-A A.Holubar, Nucl.Research Inst., Rez, Czech Republic. C-Version: WARNING - Test version under development C-Version: WILLIE original code derived from WILIT (workname: WILIT2). C-V 93/05 log-lin interpolation over last RI interval in INTERI. C-V 94/08 - input interpolation flag in PTAB option, C-V - new RINTER routine for Res.Integral interpolation. C-V 94/10 Fix filemark search in Res.Int. tables in COMPMT. C-V 94/11 Allow data insertion BEFORE a specified MAT in INCLMT. C-V 96/07 Interpret different file marks in the binary files. C-V 96/10/03 Define FormFeed with CHAR function. C-V 96/10/07 Define CHCK as Hollerith for compatibility. C-V Define " EXPORT" version as default. C-V 97/10/07 Open new files in "EXPORT" version with 'UNKNOWN' status. C-V 98/11/10 Implement option DELM to delete materials C-V Implement option CHCK to check the library integrity. C-V 98/12/22 Interpret blank keywords as comment lines. C-V 99/01/06 Allow fissile materials redefined as fission products. C-V 99/10/20 Fix undefined LER in COMPMT (found by O.Cabellos). C-V 00/03/20 Add NEWL option. C-V 01/05/22 MIXM group structure specification and description. C-M C-M Manual for Program WILLIE C-M ========================= C-M C-M 1. INTRODUCTION C-M --------------- C-M WILLIE is a service program for various operations on the WIMS-D C-M cross sections data library. It is a direct descendant of the C-M WILIT code by A.Holubar. Major modifications and extensions were C-M introduced by A.Trkov to the input formats to make them more C-M user-friendly when running in the interactive mode. Some C-M modifications were also introduced in the cross sections input C-M file formats for better compatibility with external data C-M processing codes (such as NJOY). Fairly extensive changes were C-M also introduced for clearer, commented coding. In the present C-M version of WILLIE some of the old WILIT options are not supported, C-M while some others are added, such as the capability for producing C-M cross section data files suitable as input to graphics plotting C-M programs (for example, PLOTTAB). C-M C-M 2. PURPOSE AND CHARACTERISTICS OF THE CODE C-M ------------------------------------------ C-M Program WILLIE was designed to offer the possibility of making C-M various changes in a WIMS-D library. In particular, it enables C-M to include new data file generated by a processing code such C-M as NJOY into an existing WIMS library. It is also provides C-M for extracting data for a requested material onto a new file C-M and re-inserting the data (after editing the file or otherwise) C-M back into the library with another material identification C-M number. Comparison of two data sets extracted from a WIMS C-M library by WILLIE can be made. WILLIE mostly deals with the C-M card-image form of the library but it is also used to perform C-M the transformations between the binary and the coded form of C-M the library. C-M C-M A detailed description of the binary and formatted forms of the C-M WIMS library are given in the Appendix, since this information C-M is considered useful for a better understanding of the WILLIE C-M code application and its further improvements. C-M C-M WILLIE is designed to be used from a terminal. Brief instructions C-M for defining and formatting the input parameters are given on C-M the screen. At present, WILLIE provides sixteen functions by C-M selecting the appropriate four-letter keywords (some of which C-M are incative). C-M C-M 3. INPUT DESCRIPTION C-M -------------------- C-M Every service request from WILLIE begins by entering the C-M appropriate keyword on the input line. Only the first four C-M characters are interpreted and they must be uppercase. All C-M required data files are opened explicitly by requesting the C-M appropriate filenames from input. Up to 40 characters are C-M interpreted to define filenames. Up to 20 characters are read C-M when entering numeric values. Blanks are interpreted as null, C-M so the input entries need not be right-justified. A decimal C-M point for floating point values is not mandatory. Currently C-M available services identified by keywords and additional C-M input parameters (if any) are listed below. C-M When running WILLIE in batch mode (i.e. feeding the C-M input directives from a file) it is useful to know that C-M input statements with a "blank" keyword are simply copied C-M to the default output unit. This feature can be used to C-M enter comments into the input file before activating a C-M particular servece request from WILLIE. C-M C-M KWRD. Description C-M """"""""""""""""" C-M FOBI Formatted (card-image) library is transformed into binary C-M library. Additional input parameters are: C-M - source formatted library filename, C-M - target binary library filename. C-M BIFO A binary library is transformed into a formatted library. C-M Additional input parameters are: C-M - source binary library filename, C-M - target formatted library filename. C-M COPY Data for a selected material are extracted from the C-M formatted WIMS library and stored on an output file. C-M The original library remains unchanged. Additional C-M parameters are C-M - source formatted library filename, C-M - target cross section filename. C-M - AMT1 identifier of the nuclide (including the decimal C-M resonance integral identifier), the data of which is C-M to be extracted. C-M COP1 Data for a selected P1 scattering matrix are extracted C-M from the formatted WIMS library and stored on the output C-M file. The original library remains unchanged. Additional C-M parameters are C-M - source formatted library filename, C-M - consecutive number NP1 of the P1 matrix to be extracted. C-M Zero or blank is not a valid entry. It causes a list of C-M materials in the library to be displayed. Entering a C-M negative number aborts the service request. C-M - target P1 matrix filename. C-M NOTE: At present only the standard format of the source C-M WIMS library with four P1 scattering matrices in C-M expanded format is accepted. In future this will be C-M extended to include any number of condensed or C-M uncondensed P1 scattering matrices to be processed, C-M and optionally to convert between the expanded and C-M condensed format of the matrices. In addition, it C-M will be possible to append the P1 matrix to a material C-M cross section file, rather than unconditionally open C-M a new file. C-M NEWL A starter file for a new library is generated. Additional C-M input parameters are: C-M - Filename for the new library. C-M - Total number of groups. C-M - Number of fast neutron energy groups. C-M - Number of resonance energy groups. C-M NOTE: Make sure that the energy boundaries and the spectrum C-M are entered with the INCL command before the library C-M is used. C-M INCL Data for a material extracted from a WIMS library either by C-M the WILLIE option COPY or otherwise are read and added to the C-M WIMS formatted library. If present, the burnup data, nuclide C-M cross section data, resonance integrals tables and P1 matrix C-M are processed. Additional input parameters are: C-M - source formatted library filename, C-M - source cross section filename, C-M - target formatted library filename, C-M - identification number MAT of the material to be inserted C-M where the decimal point labels the resonance tables. C-M SPECIAL CASES: C-M * If a zero (or blank) is entered, the default material C-M number on the cross section file is assumed. C-M * If the selected MAT number is already in the C-M library, an additional request for the MAT number is C-M issued. The following options are available: C-M . blank to replace the data in the WIMS library, C-M . new MAT number to insert the new data BEFORE the C-M data for the previously defined MAT number. C-M * In case of an input error, the user may wish to abort C-M the current service by entering a negative MAT number. C-M SMTR A cross section data file, such as produced by the COPY C-M service option is processed. A correction to the self- C-M scattering terms of the scattering matrices is applied such C-M that the absorption cross section and the sum of the C-M scattering matrix element for a particular group add up to C-M the transport cross section. Resonance integrals are C-M copied, but P1-matrices (if any) are skipped. Additional C-M input parameters are: C-M - source cross section filename, C-M - output cross section filename. C-M CHCK Check the library format and consistency. The following C-M tests are made: C-M * NEL,NG,NG0,NG1,NG2,NG3 (number of materials, number of C-M energy groups, fission spectrum groups, number of fast, C-M resonance and thermal groups, respectively) are checked C-M implicitly. In case of inconsistency an abnormal program C-M termination occurs due to library format errors. This C-M type of inconsistency can not be corrected. C-M * NNFD number of fissile nuclei implies all nuclides C-M with the burnup data vector greater than eight. C-M If in error, the value in the library is corrected. C-M * NNFPD number of fission product nuclei implies all C-M nuclides with the burnup data vector length of eight. C-M If in error, the value in the library is corrected. C-M * Nuclide identification numbers. A check for duplicate C-M entries is made. C-M C-M In addition, the burnup and fission yield vectors of C-M all materials are checked to make sure that the C-M addressed nuclides are burnable. A warning is printed, C-M but no corrective action can be taken. C-M Additional input parameters are: C-M - source library filename, C-M - corrected library filename. C-M COMP Two data files, such as produced by the COPY service option C-M or otherwise can be compared in tabular form. Cross sections, C-M P0 scattering matrices and resonance integrals can be C-M processed. Consistency between the scattering matrix and the C-M transport cross section is checked. A separate tolerance C-M criterion is entered for cross sections, transport C-M consistency and resonance integrals. Only the groups in C-M which the difference exceeds the tolerance limit are printed. C-M If zero is entered, all cross sections are printed. Entering C-M a negative value suppresses the groups printout, so that only C-M the maximum and the minimum difference are printed. Up to 60 C-M lines are printed per page. If a reaction type would not fit C-M onto the same page, a form feed is printed. The full output C-M is written on a file defined from input. A brief summary of C-M differences exceeding the specified threshold are written C-M on the COMPMT.LOG file. Additional input parameters are: C-M - reference cross section filename, C-M - compared cross section filename, C-M - output report filename, C-M - cross section printout tolerance [%], C-M - transport consistency tolerance [%], C-M - resonance integral tolerance [%]. C-M NOTE: Resonance integrals are tabulated and compared on a C-M grid of the Bondarenko Sig0 values. A conversion is C-M made from the WIMS background cross section SigB to C-M the Bondarenko Sig0 by subtracting Lambda*SigP from C-M SigB. On the Sig0 scale the resonance integrals are C-M comparable better since the dependence of the absorber C-M potential cross section is removed from the definition C-M to a large extent. C-M From the Sig0 grid of the reference and the compared C-M data, WILLIE first constructs a union grid. It then C-M eliminates the endpoints to avoid extrapolation for C-M either of the data sets. Since no more than 9 values C-M can be printed per line, it reduces the union grid by C-M eliminating the points which lie close together, so C-M that the final grid is as uniform as possible. C-M Resonance integrals are interpolated to the union Sig0 C-M grid using Segev interpolation, except for very low C-M Sig0 values where Root(Sig0) interpolation is applied. C-M PTAB Cross sections from a WILLIE file produced by the COPY option C-M or otherwise are tabulated in two-column format that can be C-M read by a graphics program such as PLOTTAB. Several input C-M files can be specified so that the data sets can be sorted C-M by reaction (i.e.: same cross section from different files) C-M or by material (i.e.: different cross sections from one or C-M more files). Resonance integrals are tabulated as a function C-M of the internally defined mesh of the WIMS background cross C-M section SigB for a user-specified group or all groups. C-M Additional input parameters are: C-M - reference cross section filename (up to 10 files can be C-M specified, terminated by blank), C-M - data descriptor for each data set, which appears as the C-M first record on the output data file for a data set C-M (default descriptor is the filename), C-M - reaction ident form the following list: C-M 1 - potential cross section C-M 2 - slowing down power C-M 3 - transport cross section C-M 4 - absorption cross section C-M 5 - neutron fission yield C-M 6 - fission cross section C-M 7 - outscattering cross section C-M 8 - total scattering cross section C-M 9 - group scattering matrix C-M 10 - absorption resonance integrals C-M 11 - fission resonance integrals C-M -1 - to finish C-M - output filename; entering blank implies adding data to the C-M default filename (PTABMT.CUR defined initially). By C-M entering a new filename the old output file (if open) is C-M closed and the new file opened as the default output file, C-M - required temperature [K], C-M NOTE: only when more than one temperature is present in C-M the thermal cross section or the resonance integrals C-M tables. C-M - required group C-M NOTE: only for the resonance integrals. If zero is entered C-M the resonance integrals for all resonance groups are C-M tabulated. C-M - interpolation flag C-M NOTE: only for the resonance integrals. An integer entry C-M greater than zero forces Root(SigB) interpolation C-M (except in the last interval where 1/SigB rule is C-M specified). This is equivalent to the original C-M WIMS-D/4 resonance integral interpolation rule. C-M The request for the reaction type and subsequent data entries C-M is repeated until a negative entry is encountered. C-M MIXM Data for two or more material extracted from a WIMS library C-M either by the WILLIE option COPY or otherwise are summed C-M with a weighting factor entered from input in per-cent. C-M Normally the weighting factor is the natural abundance, C-M when reconstructing elemental data from the isotopic. Note C-M that the entered abundances are not normalised to 100%. C-M Additional input parameters are: C-M a. Output (mixed) cross section filename (default MIXMAT.XSW) C-M b. WIMS ID number for the new material. C-M c. Source cross section filename. C-M d. Weighting factor for this set [%]. C-M e. Number of fast groups - requested only when the source C-M cross section set does not contain the data delimiter C-M record with appropriate specifications (default 14). C-M f. Number of resonance groups - same condition as above. C-M (default 13). C-M g. Number of thermal groups - same condition as above. C-M (default 42). C-M n. Input item pair c./d. can be repeated any number of C-M times. The list is terminated by specifying zero C-M weighting factor. C-M C-D Options of the old WILIT which are NOT implemented in WILLIE: C-D incb data for material compiled by libcom-wimlib (using fedgroup-r C-D results) - file wimmat - are read from file number n3 and C-D included into another WIMS formatted library read from the C-D file number n1, the resulting formatted library is written on C-D file number n2. burnup data for material being included are C-D copied from indicated material from the original library. C-D inp1 superseeded by INCL. C-D test consistency of transport cross sections, p0 and p1 scattering C-D matrices is tested for a material from formatted library on C-D file number n1. p1 matrix copied from WIMS library in a C-D previous run is supposed on file number n3. results of the C-D test are printed on file nout1=7. C-D redu contents of WIMS library is reduced - only the specified C-D materials are copied into the new formatted library on file C-D number n2 C-D fisp superseeded by INCL. C-D comw superseeded by COMP. C-D coms superseeded by COMP. C-D comg superseeded by COMP. C-D rtab superseeded by PTAB. C-D reca tables of resonance integrals for one material are C-D recalculated for a new mesh of sigmap values specified in C-D input using the Segev's scheme for interpolation. new tables C-D replace the old ones in the new WIMS library written on C-D file number n2. value of the new mesh points will be read C-D from input. C-M C-M 4. Brief description of the main modules in WILLIE C-M -------------------------------------------------- C-M Program WILLIE is composed of a number of modules (subroutines) C-M which perform various functions which the code offers. Some C-M additional auxilliary routines are called from the modules which C-M are described separately within the modules. C-M C-M Main modules are the following: C-M C-M BINFOR - module for the WILLIE service option BIFO. C-M FORBIN - module for the WILLIE service option FOBI. C-M COPYMT - module for the WILLIE service option COPY. C-M COP1MT - module for the WILLIE service option COP1. C-M INCLMT - module for the WILLIE service option INCL. C-M SMTXTR - module for the WILLIE service option SMTR. C-M COMPMT - module for the WILLIE service option COMP. C-M PTABMT - module for the WILLIE service option PTAB. C-M C-M 5. APPENDIX - WIMS library format C-M --------------------------------- C-M Information is gathered from various manuals on WIMS and from the C-M description of WILMA code. Note that in addition to the standard C-M WIMS/D-4 format there exists a CRNL (Canada) modification and also C-M the WIMS-E format. These are not supported by WILLIE at present. C-M C-M The binary form of WIMS library is composed of a number of records, C-M each being created by one Fortran WRITE statement. In the card- C-M image form of library the data from one record are stored on an C-M appropriate number of cards. These can be written by more than one C-M fortran WRITE starement. The records can be grouped into: C-M a. general and burnup data, C-M b. cross-section data (nuclide files), C-M c. resonance tables, C-M d. P1-scattering matrices. C-M A brief description and format (applicable to the formatted C-M library) is given below. C-M C-M a. General and burnup data C-M -------------------------- C-M The first four records define the general parameters of the library C-M and are common for all nuclides. Note that this includes the C-M fission spectrum, which can not be entered for individual fissile C-M species separately. Then follow the burnup data for each nuclide. C-M Rec. Format Data: C-M 1 (5I15) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NNFPD C-M 2 (5I15) (ID(i),i=1,NEL) C-M 3 (1P,5E15.8) (EG(i),i=1,NG+1) C-M 4 (1P,5E15.8) (CHI(i),i=1,NG0) C-M 4+i (1P,I15/(3(E15.8,I6)) NDAT1,(MTF(j),YLD(j),j=1,NDAT1/2) C-M for i = 1 to NEL C-M 5+NEL (I15) filemark C-M C-M Parameter description: C-M NEL - total number of nuclides in the library, C-M NG - total number of groups, C-M NG0 - number of groups in which fission neutrons appear, C-M NG1 - number of groups in the fast energy range, C-M NG2 - number of groups in the resonance range, C-M NG3 - number of groups in the thermal range, C-M NNFD - number of fissile nuclei in the library (NDAT1>8), C-M NNFPD - number of fission products in the library (NDAT1=8), C-M ID(i) - material identification number (integer), C-M EG(i) - energy group boundaries [eV] in descending order, C-M CHI(i) - fraction of fission neutrons born in group (i), C-M NDAT1 - burnup data length for nuclide (i) C-M = 2 for non-burnable and structural materials, C-M = 8 for burnable materials and fission products, C-M > 8 for fissile materials with fission product C-M yields. C-M MTF(1) - ident of nuclide (i), equal to ID(i), C-M MTF(2) - ident of nuclide (j) formed by capture, C-M MTF(3) - ident of nuclide (j) formed by decay, C-M MTF(4) - a flag that can have one of the following values : C-M -2 fission product with resonance tables C-M -1 fission product without resonance tables C-M 0 not fissile, no resonance tables C-M 1 not fissile with resonance tables C-M 2 fissile with resonance absorption tables C-M 3 fissile with resonance absorption, fission C-M tables C-M 4 fissile without resonance tables C-M MTF(k+4) - ident of the fission product nuclide (k=1,(NDAT1-8)/2) C-M YLD(1) - dummy (normally zero) C-M YLD(2) - fraction of captures producing nuclide MTF(2) C-M YLD(3) - decay constant of nuclide ID(i) [1/s] C-M YLD(4) - energy released per fission [J/mole*1E-24] C-M YLD(5) - fraction of nuclide MTF(2) formed by decay, C-M YLD(k+5) - fraction of nuclide MTF(k+4) produced per fission. C-M C-M C-M b. Cross-section data of the nuclide files C-M ------------------------------------------ C-M The cross sections are given first for the epithermal range and C-M are assumed temperature independent. Temperature and resonance C-M self-shielding are accounted for separately in the resonance C-M integr tables for the nuclides for which these effects are C-M important. Data in the thermal range can be tabulated for one C-M or more temperatures. It is assumed that the energy mesh in the C-M thermal groups is fine enough to account for the self-shielding C-M effects in the low-lying resonances explicitly. C-M C-M Rec. Format Data: C-M 1 (1P,I6,E15.8,6I6) ID,AW,IZ,NF,NT,NR C-M 2 (1P,5E15.8) (XP(NG1+i),i=1,NG2),(XX(NG1+i),i=1,NG2), C-M (XT(i),i=1,NG1+NG2),(XA(i),i=1,NG1+NG2), C-M (XQ(NG1+i),i=1,NG2),(GL(NG1+i),i=1,NG2) C-M 3 (1P,I6,E15.8,6I6) (XF(i),i=1,NG1+NG2),(VF(i),i=1,NG1+NG2) C-M (only if NF > 1) C-M 4 (1P,I15/(5E15.8)) NDAT2,(XSCF(i),i=1,NDAT2) C-M 5 (1P,5E15.8) (TMP(i),i=1,NT) C-M 6 (1P,5E15.8) (XT(NG12+i),i=1,NG3),(XA(NG12+i),i=1,NG3) C-M 7 (1P,5E15.8) (XF(NG12+i),i=1,NG3),(VF(NG12+i),i=1,NG3) C-M (only if NF > 1) C-M 8 (1P,I15/(5E15.8)) NDAT3,(XSCF(i),i=1,NDAT3) C-M NOTE: Records 6,7,8 are repeated for C-M NT-times for each temperature. C-M 6+3NT (I15) filemark C-M C-M Parameter description: C-M ID - material identification number (integer), C-M AW - atomic weight of the nuclide [amu], C-M IZ - atomic number of the nuclide, C-M NF - flag C-M NT - number of thermal range temperature tabulatins, C-M NR - number of resonance integral data sets, C-M NG12 = NG1+NG2 C-M XP(i) - potential scattering cross section in group (i), C-M XX(i) - slowing down power per lethargy width in group (i), C-M XT(i) - transport cross section in group (i), C-M XA(i) - absorption cross section in group (i), C-M XQ(i) - dummy (obsolete) entry, C-M GL(i) - Goldstein-Cohen parameters for group (i), C-M XF(i) - fission cross section in group (i), C-M VF(i) - neutron fission yield (nu*Sig.Fiss.) in group (i), C-M NDAT2 - length of the epithermal scattering matrix vector, C-M XSCF(i) - epithermal scattering matrix in condensed format, C-M TMP(i) - temperature list for the thermal data [K] C-M NDAT3 - length of the thermal scattering matrix vector, C-M XSCT(i) - thermal scattering matrix in condensed format. C-M C-M c. Resonance tables C-M ------------------- C-M There are up to two sets of tables for each of the NR2 resonance C-M group (depending on the NF flag). Resonance absorption integrals C-M appear in record-1 for materials with NF=1,2. Resonance neutron C-D fission yield integrals appear in record-2 if NF=3 only. C-M C-M Rec. Format Data: C-M 1 (I15/1PE15.8,2I6/ M1*M2, RID,M1,M2 C-M (1P5E15.8)) (Tres(i,i=1,M1),(SigB(j,j=1,M2) C-M ((RIa(j,i),j=1,M2),i=1,M1) C-M 2 (I15/1PE15.8,2I6/ M1*M2, RID,M1,M2 C-M (1P5E15.8)) (Tres(i),i=1,M1),(SigB(j),j=1,M2) C-M ((RIa(j,i),j=1,M2),i=1,M1) C-M C-M NOTE: In the binary library the parameter M1*M2 is not written. C-M C-M Resonance integral tables for a group are terminated by an entry C-M with (M1=1,M2=1,RID=0. and RIa=0.), followed by a filemark. C-M C-M Parameter description: C-M M1 - Number of temperatures for the resonance integrals, C-M M2 - Number of SigB values for the resonance integrals, C-M Tres(i) - Temperature at which the resonance integrals are C-M tabulated [K], C-M SigB(j) - Background cross section at which the resonance C-M integrals are tabulated [barns], C-M RIa(j,i) - Absorption resonance integral per unit lethargy at C-M temperature Tres(i) and SigB(j) [barns]. C-M RIf(j,i) - Neutron fission yield resonance integral per unit C-M lethargy at temperature Tres(i) and SigB(j) [barns]. C-M C-M C-M d. P1 scattering matrices C-M ------------------------- C-M There are four P1 scattering matrices given at the end of the WIMS C-M library. Usually they are assigned to hydrogen bound in water, C-M deuterium in heavy water, oxygen and carbon. They are given in C-M the full (expanded) form as NGxNG matrix. C-M A special feature of WILLIE is that it can read the P1 scattering C-M matrices in condensed format when the INCLMT option is used to C-M update the WIMS library. The matrix is expanded and written into C-M the library in the usual way. C-M C-M Rec. Format Data: C-M i (1P,5E15.8) (SCP1(j,i),j=1,NG) C-M C-M The matrix elements are given for all groups (i). Four P1 matrices C-M are given one after the other, without filemarks in between. The C-M file is terminated with the filemark. C-M C-M Parameter description: C-M SCP1(j,i)- P1-scattering cross section for transfer from group C-M (i) into group (j) C-M nres total number of resonance tables (0.le.nres.le.25) C-M it is the same for each resonance group and is equal C-M to the sum (over the all resonance isotopes) of C-M of the products nver(i)*ntab(i), where C-M nver(i) is the number of different versions of C-M resonance tables for i-th isotope (for example C-M the isotope with id=235 has three versions of C-M resonance tables, identified with the rid values C-M 235.2, 235.3 and 235.4) C-M and ntab(i) is the number of reactions for which C-M resonance tables for i-th isotope are given C-M (ntab=1 if only absorption, ntab=2 for absorption C-M and fission) C-M C- C* Reserve space for the real and integer work field PARAMETER (MXRW=80 000, MXIW=5 000) C* To change No.of keywords, alter parameter NKW and DATA KWRD,DSCR PARAMETER (NKW =16) CHARACTER*4 KWRD(16) ,KEY, ALDO1,ALDO2 CHARACTER*40 DSCR(16) ,DMMY(2) DIMENSION RWO(MXRW),IWO(MXIW) COMMON /DATFIL/ RDAT COMMON /TERWOR/ ITERM,ALDO(2) C* Logical filemark DATA JTERM/999999999/,ALDO1/'ALDO'/,ALDO2/'aldo'/ C* Keywords for selecting service options DATA KWRD / 'BIFO', 'FOBI', 'COPY', 'COP1', 'INCL' 1 , 'NEWL', 'DELM', 'SMTR', 'CHCK', 'INCB' 2 , 'COMP', 'PTAB', 'MIXM', 'RTAB', 'RECA' 3 , 'END '/ C* Service options description for output DATA DSCR 1 / ' - Binary to formatted conversion ' 2 , ' - Formatted to binary conversion ' 3 , ' - Extract nuclide data from library ' 4 , ' - Extract P1 matrix from library ' 5 , ' - Insert new data into library ' 6 , ' - Make a starter file for a new library' 7 , ' - Delete a material from the library ' 8 , ' - Transport-correct scattering matrix ' 9 , ' - Check library integrity ' * , ' - inactive INCB' 1 , ' - Compare two x-sect data sets ' 2 , ' - Tabulate data for plotting ' 3 , ' - Mix x-sect to make a new material ' 4 , ' - inactive RTAB' 5 , ' - inactive RECA' * , ' - Finish all processing. '/ C* File units DATA LKB,LTT / 5, 6 / ITERM=JTERM READ(ALDO1,501) ALDO(1) READ(ALDO2,501) ALDO(2) C* Select the WILLIE service option 10 WRITE(LTT,603) (KWRD(I),DSCR(I), I=1,NKW) 11 WRITE(LTT,691) WRITE(LTT,691) '$Enter the selected option : ' 12 READ (LKB,691) DMMY KEY=DMMY(1)(1:4) IF(KEY.EQ.' ') THEN C* Copy input records with blank keyword to output as comments WRITE(LTT,691) DMMY GO TO 12 END IF C* Identify the selected option IK=NKW+1 DO 14 I=1,NKW IF(KEY.EQ.KWRD(I)) IK=I 14 CONTINUE IF(IK.GE.1 .AND. IK.LE.NKW) THEN WRITE(LTT,691) WRITE(LTT,604) KEY,DSCR(IK) WRITE(LTT,691) GO TO ( 18, 20, 30, 40, 50, 60, 70, 80, 90,100 1 ,110,120,130,140,150 ,200), IK END IF C* Print warning on Illegal entry - keyword not defined WRITE(LTT,609) KEY GO TO 10 C* Binary to formatted conversion 18 CALL BINFOR(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Formatted to binary conversion 20 CALL FORBIN(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Extract nuclide data from library 30 CALL COPYMT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Extract a P1 scattering matrix from library 40 CALL COP1MT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Add or replace nuclide data in the library 50 CALL INCLMT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Make a starter file for a new library 60 CALL NEWLIB(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Delete a material from the library 70 CALL DELEMT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Transport-correct self-scattering term of the scattering matrix 80 CALL SMTXTR(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Check the library integrity 90 CALL CHKLIB(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 100 CONTINUE C CALL INCLMB GO TO 150 C* Make comparison tables for x-sect and RI for two WILIT files 110 CALL COMPMT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Tabulate nuclide data for plotting purposes 120 CALL PTABMT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 C* Mix cross sections to make a new material 130 CALL MIXMAT(LKB,LTT,RWO,IWO,MXRW,MXIW) GO TO 11 140 CONTINUE 150 CONTINUE C* All processing completed 200 STOP 'WILLIE Completed' 501 FORMAT(A4) 603 FORMAT(/' WILLIE - WIMS Library Utility'/ ! ' ============================='/ 1 ' The following service options are available:'/ 1 (3X,A4,A40)) 604 FORMAT(' Run ',A4,A40) 609 FORMAT(' WILLIE Error: unknown keyword "',A4,'"') 691 FORMAT(2A40) END SUBROUTINE ERRWAR(LTT,IER) C-Title : ERRWAR Subroutine C-Purpose: Print error warning on incorrect library reading WRITE(LTT,2001) IER STOP 2001 FORMAT(6X,'ERROR DURING READING LIBRARY IN PLACE NUMBER ',I2) END SUBROUTINE XSMSUM(NG,ISS,XSM,XSO) C-Title : C-Purpose: Sum elements of the scattering matrix C-Description: C-D Scattering matrix elements packed in XSM in condensed WIMS C-D format are summed by rows and stored in XSO. If ISS is C-D greater than zero, the self-scattering component is excluded C-D from the sum. C-D DIMENSION XSM(*),XSO(*) L =1 DO 40 I=1,NG S =0. NS=XSM(L )+0.1 NE=XSM(L+1)+0.1 DO 20 J=1,NE IF(J.NE.NS .OR. ISS.EQ.0) S =S + XSM(L+1+J) 20 CONTINUE XSO(I)=S L =L+2+NE 40 CONTINUE RETURN END SUBROUTINE PROCXS(RWO,IWO,LF1,LF2,NG1,NG2,NG3,NF,NT) C-Title : PROCXS Subroutine C-Purpose: Copy(LF2>0) or skip(LF2=0) the WIMS cross section data DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) NG12=NG1+NG2 KLN =2*NG12+4*NG2 C* Epithermal cross sections READ(LF1,702) (RWO(I),I=1,KLN) IF(LF2.GT.0) WRITE(LF2,702) (RWO(I),I=1,KLN) IF(ABS(NF).LE.1) GO TO 40 C* Fission yield and fission cross section KLN=2*NG12 READ(LF1,702) (RWO(I),I=1,KLN) IF(LF2.GT.0 .AND. 1 NF .GT.1) WRITE(LF2,702) (RWO(I),I=1,KLN) C* Epithermal scattering matrix 40 READ(LF1,705) NDAT2,(RWO(I),I=1,NDAT2) IF(LF2.GT.0) WRITE(LF2,705) NDAT2,(RWO(I),I=1,NDAT2) C* Temperature list for thermal cross sections READ(LF1,702) (RWO(I),I=1,NT) IF(LF2.GT.0) WRITE(LF2,702) (RWO(I),I=1,NT) DO 60 K=1,NT KLN=2*NG3 READ(LF1,702) (RWO(I),I=1,KLN) IF(LF2.GT.0) WRITE(LF2,702) (RWO(I),I=1,KLN) IF(ABS(NF).LE.1) GO TO 50 READ(LF1,702) (RWO(I),I=1,KLN) IF(LF2.GT.0 .AND. 1 NF .GT.1) WRITE(LF2,702) (RWO(I),I=1,KLN) 50 READ(LF1,705) NDAT3,(RWO(I),I=1,NDAT3) IF(LF2.GT.0) WRITE(LF2,705) NDAT3,(RWO(I),I=1,NDAT3) 60 CONTINUE RETURN 701 FORMAT(5I15) 702 FORMAT(1P5E15.8) 705 FORMAT(1P,I15/(5E15.8)) END *DECK BINFOR SUBROUTINE BINFOR(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : BINFOR Subroutine C-Purpose: Convert binary form of WIMS library to formatted C- CHARACTER*40 FLWI,FLWO DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Logical file units for the binary (source) and formatted (target) lib DATA LF1,LF2 / 1, 2 / C* Preset dummy file delimiter DATA CHDM/.777777/ C* Define the filenames 12 WRITE(LTT,691) '$Enter source binary library filename : ' READ (LKB,691) FLWI C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=12 C 1 ,FORM='UNFORMATTED',READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=12,FORM='UNFORMATTED') C***** EXPORT ***** 14 WRITE(LTT,691) '$Enter formatted library filename : ' READ (LKB,691) FLWO C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLWO,STATUS='NEW',ERR=14 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** IF(FLWI.EQ.FLWO) GO TO 14 OPEN (UNIT=LF2,FILE=FLWO,STATUS='UNKNOWN',ERR=14) C***** EXPORT ***** C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' Source WIMS-D binary library : ',FLWI WRITE(LTT,691) ' Output WIMS-D formatted library : ',FLWO WRITE(LTT,691) C* C* Begin processing the files - write the header indices READ (LF1) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD WRITE(LF2,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD C* List of nuclides READ (LF1) (IWO(I),I=1,NEL) WRITE(LF2,701) (IWO(I),I=1,NEL) C* Group boundaries KLN=NG+1 READ (LF1) (RWO(I),I=1,KLN) WRITE(LF2,702) (RWO(I),I=1,KLN) C* Fission spectrum READ (LF1) (RWO(I),I=1,NG0) WRITE(LF2,702) (RWO(I),I=1,NG0) C* Burnup data for each nuclide DO 20 J=1,NEL READ (LF1) NDAT1,IWO(1),(RWO(I),IWO(I),I=2,NDAT1/2) RWO(1)=0. WRITE(LF2,701) NDAT1 20 WRITE(LF2,703) (RWO(I/2),IWO(I/2),I=2,NDAT1,2) C* Begin processing the nuclide files - check for logical file mark IER=1 CHCK=CHDM READ (LF1,ERR=21) CHCK 21 IF(CHCK.NE.ALDO(1).AND.CHCK.NE.ALDO(2)) GO TO 190 WRITE(LF2,701) ITERM C* Loop over all nuclides DO 60 J=1,NEL C* Identifier READ (LF1) ID,AW,IZ,NF,NT,NR WRITE(LF2,704) ID,AW,IZ,NF,NT,NR KLN=4*NG2+2*(NG1+NG2) C* Absorption, transport slow.dwn.pwr, Lambda, Ksi, sg.pot READ (LF1) (RWO(I),I=1,KLN) WRITE(LF2,702) (RWO(I),I=1,KLN) IF(NF.LE.1) GO TO 42 C* If nuclide is fissile, process fission yield & fission cross sect KLN=2*(NG1+NG2) READ (LF1) (RWO(I),I=1,KLN) WRITE(LF2,702) (RWO(I),I=1,KLN) C* Epithermal scattering matrix 42 READ (LF1) NDAT2,(RWO(I),I=1,NDAT2) WRITE(LF2,705) NDAT2,(RWO(I),I=1,NDAT2) C* Thermal range temperatures READ (LF1) (RWO(I),I=1,NT) WRITE(LF2,702) (RWO(I),I=1,NT) C* Loop over all temperatures in the thermal range DO 50 K=1,NT KLN=2*(NG-NG1-NG2) C* Thermal absorption and transport READ (LF1) (RWO(I),I=1,KLN) WRITE(LF2,702) (RWO(I),I=1,KLN) IF(NF.LE.1) GO TO 44 C* If nuclide is fissile, process fission yield & fission cross sect READ (LF1) (RWO(I),I=1,KLN) WRITE(LF2,702) (RWO(I),I=1,KLN) C* Thermal scattering matrix 44 READ (LF1) NDAT3,(RWO(I),I=1,NDAT3) 50 WRITE(LF2,705) NDAT3,(RWO(I),I=1,NDAT3) IER=2 CHCK=CHDM READ (LF1,ERR=51) CHCK 51 IF(CHCK.NE.ALDO(1).AND.CHCK.NE.ALDO(2)) GO TO 190 WRITE(LF2,701) ITERM 60 CONTINUE C* Begin processing the resonance files - loop over all res. groups DO 90 J=1,NG2 C* Implied loop over all res.tables in a group 70 READ (LF1) RID,M1,M2,(RWO(I),I=1,M1), * (RWO(I+1000),I=1,M2), * ((RWO(M1*(II-1)+I+2000),I=1,M1),II=1,M2) M12=M1*M2 WRITE(LF2,706) M12,RID,M1,M2,(RWO(I),I=1,M1), * (RWO(I+1000),I=1,M2), * (RWO(I+2000),I=1,M12) IF(RID.EQ.0.) GO TO 80 GO TO 70 C* Grout resonance tables processed - check for file mark 80 IER=3 CHCK=CHDM READ (LF1,ERR=81) CHCK 81 IF(CHCK.NE.ALDO(1).AND.CHCK.NE.ALDO(2)) GO TO 190 WRITE(LF2,701) ITERM 90 CONTINUE C* Begin processing the P1 scattering matrices DO 100 J=1,4 DO 100 K=1,NG READ (LF1) (RWO(I),I=1,NG) 100 WRITE(LF2,702) (RWO(I),I=1,NG) IER=4 CHCK=CHDM READ (LF1,ERR=101) CHCK 101 IF(CHCK.NE.ALDO(1).AND.CHCK.NE.ALDO(2)) GO TO 190 WRITE(LF2,701) ITERM C* Processing terminated successfully CLOSE(UNIT=LF1) CLOSE(UNIT=LF2) RETURN C* Library read error 190 CALL ERRWAR(LTT,IER) RETURN 601 FORMAT(' Formatted WIMS library written on : ',A40) 691 FORMAT(2A40) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(3(1P,E15.8,I6)) 704 FORMAT(I6,1P,E15.8,4I6) 705 FORMAT(I15/(1P,5E15.8)) 706 FORMAT(I15/1P,E15.8,2I6/(1P,5E15.8)) END *DECK FORBIN SUBROUTINE FORBIN(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : FORBIN Subroutine C-Purpose: Convert formatted form of wims library to binary C- CHARACTER*80 REC CHARACTER*40 FLWI,FLWO DIMENSION RWO(*),IWO(*),NP1(4) COMMON /TERWOR/ ITERM,ALDO(2) C* Logical file units for the binary (source) and formatted (target) lib DATA LF1,LF2 / 1, 2 / C* Define the filenames 12 WRITE(LTT,691) '$Enter source formatted library filenm: ' READ (LKB,691) FLWI C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=12,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=12) C***** EXPORT ***** 14 WRITE(LTT,691) '$Enter binary library filename : ' READ (LKB,691) FLWO C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLWO,STATUS='NEW',ERR=14,FORM='UNFORMATTED') C***** VAX ***** C***** EXPORT ***** IF(FLWI.EQ.FLWO) GO TO 14 OPEN (UNIT=LF2,FILE=FLWO,STATUS='UNKNOWN',ERR=14 1 ,FORM='UNFORMATTED') C***** EXPORT ***** C* Define P1 scattering matrix selection ISRT=0 JP1 =1 WRITE(LTT,691) ' First P1 matrix (default first 4) : ' READ (LKB,892,END=16) NP1(1) IF(NP1(1).GT.0) GO TO 17 16 NP1(1)=1 NP1(2)=2 NP1(3)=3 NP1(4)=4 GO TO 19 C* Explicitly defined positions 17 JP1=JP1+1 18 WRITE(LTT,691) ' Next P1 matrix : ' READ (LKB,892,ERR=18) NP1(JP1) IF(NP1(JP1).EQ.0 .OR. NP1(JP1).EQ.NP1(JP1-1)) THEN WRITE(LTT,691) ' Illegal entry - redo from start ' GO TO 18 END IF IF(NP1(JP1).LT.NP1(JP1-1)) ISRT=1 IF(JP1.LT.4) GO TO 17 19 CONTINUE C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' Source WIMS-D formatted library : ',FLWI WRITE(LTT,691) ' Output WIMS-D binary library : ',FLWO WRITE(LTT,691) C* C* Begin processing the files - write the header indices READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD WRITE(LF2) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD C* List of nuclides READ (LF1,701) (IWO(I),I=1,NEL) WRITE(LF2) (IWO(I),I=1,NEL) C* Group boundaries KLN=NG+1 READ (LF1,702) (RWO(I),I=1,KLN) WRITE(LF2) (RWO(I),I=1,KLN) C* Fission spectrum READ (LF1,702) (RWO(I),I=1,NG0) WRITE(LF2) (RWO(I),I=1,NG0) C* Burnup data for each nuclide DO 20 J=1,NEL READ (LF1,701) NDAT1 READ (LF1,703) (RWO(I),IWO(I),I=1,NDAT1/2) 20 WRITE(LF2) NDAT1,IWO(1),(RWO(I),IWO(I),I=2,NDAT1/2) C* Begin processing the nuclide files - check for logical file mark IER=1 READ(LF1,701,ERR=190) ICHK IF(ICHK.NE.ITERM) GO TO 190 WRITE(LF2) ALDO(1) C* Loop over all nuclides DO 60 J=1,NEL C* Identifier READ (LF1,704) ID,AW,IZ,NF,NT,NR WRITE(LF2) ID,AW,IZ,NF,NT,NR KLN=4*NG2+2*(NG1+NG2) C* Absorption, transport slow.dwn.pwr, Lambda, Ksi, sg.pot READ (LF1,702) (RWO(I),I=1,KLN) WRITE(LF2) (RWO(I),I=1,KLN) IF(NF.LE.1) GO TO 42 C* If nuclide is fissile, process fission yield & fission cross sect KLN=2*(NG1+NG2) READ (LF1,702) (RWO(I),I=1,KLN) WRITE(LF2) (RWO(I),I=1,KLN) C* Epithermal scattering matrix 42 READ (LF1,705) NDAT2,(RWO(I),I=1,NDAT2) WRITE(LF2) NDAT2,(RWO(I),I=1,NDAT2) C* Thermal range temperatures READ (LF1,702) (RWO(I),I=1,NT) WRITE(LF2) (RWO(I),I=1,NT) C* Loop over all temperatures in the thermal range DO 50 K=1,NT KLN=2*(NG-NG1-NG2) C* Thermal absorption and transport READ (LF1,702) (RWO(I),I=1,KLN) WRITE(LF2) (RWO(I),I=1,KLN) IF(NF.LE.1) GO TO 44 C* If nuclide is fissile, process fission yield & fission cross sect READ (LF1,702) (RWO(I),I=1,KLN) WRITE(LF2) (RWO(I),I=1,KLN) C* Thermal scattering matrix 44 READ (LF1,705) NDAT3,(RWO(I),I=1,NDAT3) WRITE(LF2) NDAT3,(RWO(I),I=1,NDAT3) 50 CONTINUE IER=2 READ (LF1,701,ERR=190) ICHK IF(ICHK.NE.ITERM) GO TO 190 WRITE(LF2) ALDO(1) 60 CONTINUE C* Begin processing the resonance files - loop over all res. groups DO 90 J=1,NG2 C* Implied loop over all res.tables in a group 70 READ (LF1,706) M12,RID,M1,M2,(RWO(I),I=1,M1), * (RWO(I+1000),I=1,M2), * (RWO(I+2000),I=1,M12) WRITE(LF2) RID,M1,M2,(RWO(I),I=1,M1), * (RWO(I+1000),I=1,M2), * ((RWO(M1*(II-1)+I+2000),I=1,M1),II=1,M2) IF(RID.EQ.0.) GO TO 80 GO TO 70 C* Group resonance tables processed - check for file mark 80 IER=3 READ (LF1,701,ERR=190) ICHK IF(ICHK.NE.ITERM) GO TO 190 WRITE(LF2) ALDO(1) 90 CONTINUE C* C* Begin processing the P1 scattering matrices IF(ISRT.NE.0) THEN ISRT=LF2 LF2 =31 OPEN (UNIT=LF2,FILE='FOBI.TMP',STATUS='UNKNOWN' 1 ,FORM='UNFORMATTED') END IF IP1=0 JP1=1 100 IP1=IP1+1 C* Identify the format on the original library (full or condensed) C* Next card may be the last filemark, matrix length or data IER=56 101 READ (LF1,692,END=190) REC READ (REC,701,ERR=150) KL1 C* Filemark found - read next data set IF(KL1.EQ.ITERM) GO TO 101 C* Case: Matrix length on record - read condensed P1 matrix L1 =NG*NG+1 DO 141 J=1,L1 RWO(J)=0 141 CONTINUE READ (LF1,702) (RWO(L1-1+J),J=1,KL1) C* Expand the matrix JG=0 145 JG=JG+1 NU =RWO(L1 )+0.1 NS =RWO(L1+1)+0.1 DO 148 J=1,NS RWO((JG-1)*NG + JG-NU-1+J)=RWO(L1+1+J) 148 CONTINUE L1=L1+2+NS IF(JG.LT.NG) GO TO 145 GO TO 160 C* Case: read expanded P1 matriX 150 READ (REC,702) (RWO(J),J=1,5) READ (LF1,702) (RWO(J),J=6,NG) DO 152 K=2,NG READ (LF1,702) (RWO((K-1)*NG+J),J=1,NG) 152 CONTINUE C* Write the P1 scattering matrix, if selected 160 DO 163 J=1,4 IF(NP1(J).EQ.IP1) THEN DO 162 K=1,NG WRITE(LF2) (RWO((K-1)*NG+I),I=1,NG) 162 CONTINUE IF(ISRT.EQ.0) 1 WRITE(LTT,693) ' Written P1-matrix from position ',IP1 JP1=JP1+1 END IF 163 CONTINUE IF(JP1.LE.4) GO TO 100 C* Sort the P1 matrices if necessary IF(ISRT.NE.0) THEN LSC=LF2 LF2=ISRT DO 168 JP1=1,4 REWIND LSC C* Find the next P1 matrix to write IP1=1000 DO 164 J=1,4 IF(NP1(J).LT.IP1) IP1=J 164 CONTINUE C* Skip preceeding P1 matrices and copy IP1-th one DO 166 J=1,IP1 DO 165 K=1,NG READ (LSC) (RWO((K-1)*NG+I),I=1,NG) IF(J.EQ.IP1) WRITE(LF2) (RWO((K-1)*NG+I),I=1,NG) 165 CONTINUE 166 CONTINUE WRITE(LTT,693) ' Written P1-matrix from position ' 1 ,NP1(IP1) NP1(IP1)=NP1(IP1)+1000 168 CONTINUE END IF C* Check the last filemark C... IER=4 C... READ (LF1,701,ERR=190) ICHK C... IF(ICHK.NE.ITERM) GO TO 190 WRITE(LF2) ALDO(1) C* Processing terminated successfully CLOSE(UNIT=LF1) CLOSE(UNIT=LF2) RETURN C* Library read error 190 CALL ERRWAR(LTT,IER) RETURN 691 FORMAT(2A40) 692 FORMAT(A80) 693 FORMAT(A40,I4) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,I15/E15.8,2I6/(5E15.8)) 892 FORMAT(BN,I10) END *DECK COPYMT SUBROUTINE COPYMT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : COPYMT Subroutine C-Purpose: Extract data for selected material from wims library C- CHARACTER*40 FLNM DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Logical file units for source library and output file DATA LF1,LF2 / 1, 2 / C* Define the filenames 12 WRITE(LTT,691) '$Enter source formatted library filenm: ' READ (LKB,691) FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12) C***** EXPORT ***** 14 WRITE(LTT,691) '$Enter nuclide data output filename : ' READ (LKB,691) FLNM C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLNM,STATUS='NEW',ERR=14 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF2,FILE=FLNM,STATUS='UNKNOWN',ERR=14) C***** EXPORT ***** C* Read the header indices READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD C* Read the list of nuclides READ (LF1,701) (IWO(I),I=1,NEL) C* Define the data set to be extracted 15 WRITE(LTT,691) '$ Enter the selected material ident. : ' READ (LKB,501,ERR=17) RMT IF(RMT.EQ.0) GO TO 100 C* Check if material is on the list MAT=RMT+0.001 DO 16 I=1,NEL IF(IWO(I).EQ.MAT) GO TO 18 16 CONTINUE WRITE(LTT,601) MAT 17 WRITE(LTT,602) (IWO(J),J=1,NEL) GO TO 15 18 CONTINUE C* Skip group boundaries KLN=NG+1 READ (LF1,702) (RWO(I),I=1,KLN) C* Skip fission spectrum READ (LF1,702) (RWO(I),I=1,NG0) C* Extract burnup data for selected nuclide WRITE(LF2,701) ITERM, 3 DO 20 J=1,NEL RWO(1)=0. READ (LF1,701) NDAT1 READ (LF1,703) (RWO(I),IWO(I),I=1,NDAT1,2) IF(IWO(1).NE.MAT) GO TO 20 WRITE(LF2,701) NDAT1 WRITE(LF2,703) (RWO(I),IWO(I),I=1,NDAT1,2) 20 CONTINUE C* Beginning of the nuclide files - check for logical file mark IER=1 READ (LF1,701,ERR=190) JTERM IF(JTERM.NE.ITERM) GO TO 190 C* Loop over all nuclides DO 60 J=1,NEL LFF=0 C* Read the nuclide header card READ (LF1,704) ID,AW,IZ,NF,NT,NR C* Only one set of resonance parameters is copied IF(NR.GT.1) NR =1 C* If selected nuclide then write to output, else skip to next file mark IF(ID.NE.MAT) GO TO 42 LFF=LF2 WRITE(LF2,701) ITERM, 4,NG1,NG2,NG3 WRITE(LF2,704) ID,AW,IZ,NF,NT,NR 42 CALL PROCXS(RWO,IWO,LF1,LFF,NG1,NG2,NG3,NF,NT) 45 IER=110 READ (LF1,701) JTERM IF(JTERM.NE.ITERM) GO TO 190 60 CONTINUE C* Begin processing the resonance files - loop over all res. groups DO 90 J=1,NG2 NRM=0 C* Implied loop over all res.tables in a group 70 READ (LF1,706) M12,RID,M1,M2,(RWO(I ),I=1,M1), * (RWO(I+1000),I=1,M2), * (RWO(I+2000),I=1,M12) IF( RID .EQ.0. ) GO TO 80 IF(ABS(RID-RMT).GT.0.01) GO TO 70 NRM=NRM+1 IF(NRM.EQ.1) WRITE(LF2,701) ITERM, 5,NG1+J,NG1,NG2 WRITE(LF2,706) M12,RID,M1,M2,(RWO(I ),I=1,M1), * (RWO(I+1000),I=1,M2), * (RWO(I+2000),I=1,M12) GO TO 70 C* Group resonance tables processed - check the file mark 80 READ (LF1,701) JTERM IF(JTERM.NE.ITERM) GO TO 190 90 CONTINUE C* P1-scattering matrices are ignored WRITE(LF2,701) ITERM REWIND LF1 GO TO 15 C* All material processed 100 CLOSE(UNIT=LF1) CLOSE(UNIT=LF2) RETURN C* Library read error 190 CALL ERRWAR(LTT,IER) RETURN 501 FORMAT(BN,F20.0) 601 FORMAT(' WILLIE ERROR - Material',I6,' not found') 602 FORMAT(' Select from the list below:'/ (6X,10I6) ) 603 FORMAT(' Processing material',F8.1,' completed'/) 691 FORMAT(2A40) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,I15/E15.8,2I6/(5E15.8)) END *DECK COP1MT SUBROUTINE COP1MT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : COP1MT Subroutine C-Purpose: Extract P1 matrix for a selected material from WIMS library C- CHARACTER*40 FLNM DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Logical file units for source library and output file DATA LF1,LF2 / 1, 2 / C* Define the filenames 12 WRITE(LTT,691) '$Enter source formatted library filenm: ' READ (LKB,691) FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12) C***** EXPORT ***** C* Read the header indices READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD C* Read the list of nuclides READ (LF1,701) (IWO(I),I=1,NEL) C* Define the data set to be extracted 15 WRITE(LTT,691) '$ Enter consecutive P1-matrix number : ' READ (LKB,501,ERR=17) NP1 IF(NP1.LT.0) GO TO 100 IF(NP1.GT.0) GO TO 20 C* Invalid entry - display materials and inquire again 17 WRITE(LTT,602) (IWO(J),J=1,NEL) WRITE(LTT,691) ' Invalid entry - input -1 to quit, or ' GO TO 15 20 WRITE(LTT,691) '$Enter the data output filename : ' READ (LKB,691) FLNM C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLNM,STATUS='NEW',ERR=20 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF2,FILE=FLNM,STATUS='UNKNOWN',ERR=20) C***** EXPORT ***** C* Skip the nuclide files boundaries CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) DO 22 I=1,NEL CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) 22 CONTINUE C* Skip the resonance files DO 24 I=1,NG2 CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) 24 CONTINUE C* Copy the selected P1 matrix to output file JP1=0 C* At present it works only for standard (expanded) library format 30 JP1=JP1+1 DO 32 IG=1,NG READ (LF1,702,ERR=190) (RWO(J),J=1,NG) IF(JP1.EQ.NP1) WRITE(LF2,702) (RWO(J),J=1,NG) 32 CONTINUE IF(JP1.LT.4) GO TO 30 C* All material processed 100 WRITE(LF2,701) ITERM CLOSE(UNIT=LF1) CLOSE(UNIT=LF2) 120 RETURN C* Library read error 190 WRITE(LTT,691) ' ERROR - only standard format supported ' RETURN 501 FORMAT(BN,I20) 602 FORMAT(' The library contains the following materials:'/(6X,10I6)) 691 FORMAT(2A40) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) END *DECK COMPMT SUBROUTINE COMPMT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : COMPMT Subroutine C-Purpose: Compare two WIMS library cross section data sets C-Description: C-D A cross section file such as produced by WILLIE option COPY C-D is processed. It generously allows labelled, unlabelled C-D as well as no-filemark formats. Cross sections are compared C-D group-by-group, their values together with relative and C-D absolute differences are printed. Maximum and minimum C-D relative differences are recorded explicitly. C-D Comparison of temperature dependent data is limited in the C-D sense that: C-D - no interpolation in temperature is provided, C-D - data are compared if temperatures differ by < 10 deg.C, C-D - no provision is made for the case when the reference C-D data set is given on a fine temperature mesh. C-D CHARACTER*1 FF CHARACTER*10 RITY(3) CHARACTER*40 FLN1,FLN2,FLNO,FLNL,FLNM, COMM DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Default group structure (fast/resonance/thermal) DATA NG1,NG2,NG3/ 14, 13, 42 / C* Alternative group structure DATA MG1,MG2,MG3/ 45, 47, 80 / C* Default output filename DATA FLNO /'COMPMT.OUT'/ 1 FLNL /'COMPMT.LOG'/ C* Local file units DATA LF1,LF2,LOU,LER / 1, 2, 3, 7 / C* Resonance integral types DATA RITY /'ABSORPTION','FISSION ','??????????'/ C* Define output page size (lines) and current page position DATA NP,IP /62,0/ C* Form feed character FF =CHAR(12) C* Define the reference cross section filename LFF=LF1 WRITE(LTT,691) ' REFERENCE CROSS SECTION DATA: ' GO TO 12 11 WRITE(LTT,691) ' ERROR reading file : ',FLN2 12 WRITE(LTT,691) '$Enter the cross section data filename: ' READ (LKB,691) FLN2 C***** VAX ***** C OPEN(UNIT=LFF,FILE=FLN2,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN(UNIT=LFF,FILE=FLN2,STATUS='OLD',ERR=11) C***** EXPORT ***** C* Try to interpret the filemark record READ (LFF,701,ERR=26) JTERM, NTYP,JG1,JG2,JG3 IF(NTYP.EQ.0) GO TO 32 20 IF(NTYP - 4 ) 30,32,38 C* No labelled filemark present 26 REWIND LFF GO TO 40 C* Inappropriate labelled filemark found 30 READ (LFF,701,ERR=30) JTERM, NTYP,JG1,JG2,JG3 GO TO 20 C* Labelled filemark for cross section data found 32 IF(LFF.NE.LF2) GO TO 36 C* - Reference data file IF(JG1.GT.0) NG1=JG1 IF(JG2.GT.0) NG2=JG2 IF(JG3.GT.0) NG3=JG3 GO TO 40 C* - Compared data file 36 IF(JG1.GT.0 .AND. JG1.NE.NG1) GO TO 37 IF(JG2.GT.0 .AND. JG2.NE.NG2) GO TO 37 IF(JG3.GT.0 .AND. JG3.NE.NG3) GO TO 37 GO TO 40 C* Error condition on group structure 37 WRITE(LTT,691) ' Incompatible group structure on file : ',FLN2 GO TO 12 C* No cross section data found 38 WRITE(LTT,691) ' No cross section data on file : ',FLN2 GO TO 12 C* Define the compared cross section filename 40 IF(LFF.EQ.LF2) GO TO 42 C* Check the group structure IF(JG1.EQ.0) THEN NGR=NG1+NG2+NG3 WRITE(LTT,693) ' Assumed number of groups ',NGR WRITE(LTT,693) '$ Enter number of groups to redefine : ' READ (LKB,695) JGR IF (JGR.EQ.0 .OR. JGR.EQ. 69) THEN C* 69-groups default ELSE IF(JGR.EQ.172) THEN C* 172-grous alternative group structure NG1=MG1 NG2=MG2 NG3=MG3 NGR=JGR ELSE C* Arbitrary group structure WRITE(LTT,693) '$ Enter number of fast groups : ' READ (LKB,695) NG1 WRITE(LTT,693) '$ Enter number of resonance groups : ' READ (LKB,695) NG2 NGR=JGR NG3=NGR-NG1-NG2 END IF END IF FLN1=FLN2 LFF =LF2 WRITE(LTT,691) ' COMPARED CROSS SECTION DATA: ' GO TO 12 C* Both data files have been initialized, open the output 42 WRITE(LTT,691) '0Default output report filename : ',FLNO WRITE(LTT,691) '$ Enter new filename to redefine : ' READ (LKB,691) FLNM IF(FLNM.NE.' ') FLNO=FLNM C***** VAX ***** C OPEN (UNIT=LOU,FILE=FLNO,STATUS='NEW',CARRIAGECONTROL='LIST') C OPEN (UNIT=LER,FILE=FLNL,STATUS='NEW',CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LOU,FILE=FLNO,STATUS='UNKNOWN') OPEN (UNIT=LER,FILE=FLNL,STATUS='UNKNOWN') C***** EXPORT ***** C* Enter the threshold for comparison output WRITE(LTT,691) '$Enter threshold diff.for x-sect. [%]: ' READ (LKB,694) EPS WRITE(LTT,691) '$Enter threshold diff.for trans. [%]: ' READ (LKB,694) EPZ WRITE(LTT,691) '$Enter threshold diff.for R.I. [%]: ' READ (LKB,694) EPR NG12=NG1 +NG2 NG = NG12+NG3 C* Begin processing the cross sections 400 READ(LF1,704) MT1,AW1,IZ1,NF1,NT1,NR1 READ(LF2,704) MT2,AW2,IZ2,NF2,NT2,NR2 C* Label the output file WRITE(LOU,607) WRITE(LOU,608) ' Reference',MT1,AW1,IZ1,NF1,NT1,NR1,FLN1 WRITE(LOU,608) ' Compared ',MT2,AW2,IZ2,NF2,NT2,NR2,FLN2 WRITE(LER,607) WRITE(LER,608) ' Reference',MT1,AW1,IZ1,NF1,NT1,NR1,FLN1 WRITE(LER,608) ' Compared ',MT2,AW2,IZ2,NF2,NT2,NR2,FLN2 IP =8 IF(EPS.LE.0) GO TO 402 IF(EPS.LE.999.9) WRITE(LOU,606) EPS IF(EPS.LE.999.9) WRITE(LER,606) EPS WRITE(LER,691) IP =IP+2 C* C* Epithermal cross sections processing C* 402 KLN=4*NG2+2*NG12 KLF=2*NG12 LS =1 C* Read the reference epithermal data L1 =LS READ(LF1,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF1.GE.2) READ(LF1,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ(LF1,705) KL1,(RWO(LS-1+I),I=1,KL1) LS =LS+KL1+NG12 C* Read the compared epithermal data L2=LS READ(LF2,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF2.GE.2) READ(LF2,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ(LF2,705) KL2,(RWO(LS-1+I),I=1,KL2) LS =LS+KL2+NG12 C* Perform comparison for epithermal cross sections LX1=L1 LX2=L2 COMM = ' POTENTIAL SCATTERING CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG1,NG2,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG2 LX2=LX2+NG2 COMM = ' SLOWING DOWN POWER ' CALL XSDIFF(LOU,LER 1 ,NG1,NG2,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG2 LX2=LX2+NG2 COMM = ' TRANSPORT CROSS SECTION ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG12 LX2=LX2+NG12 COMM = ' ABSORPTION CROSS SECTION ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG12 + 2*NG2 LX2=LX2+NG12 + 2*NG2 C* Epithermal fission yield and fission cross section IF(NF1.LT.2 .OR. NF2.LT.2) GO TO 412 COMM = ' NEUTRON FISSION YIELD ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) 412 LX1=LX1+NG12 LX2=LX2+NG12 IF(NF1.LT.2 .OR. NF2.LT.2) GO TO 414 COMM = ' FISSION CROSS SECTION ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) 414 LX1=LX1+NG12 LX2=LX2+NG12 C* Epithermal scattering matrix (compare outscatter only) LO1=LX1+KL1 LO2=LX2+KL2 CALL XSMSUM(NG12,1,RWO(LX1),RWO(LO1)) CALL XSMSUM(NG12,1,RWO(LX2),RWO(LO2)) COMM = ' OUTSCATTERING CROSS SECTION ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LO1),RWO(LO2),RWO(LS),EPS,NP,IP) CALL XSMSUM(NG12,0,RWO(LX1),RWO(LO1)) CALL XSMSUM(NG12,0,RWO(LX2),RWO(LO2)) COMM = ' TOTAL SCATTERING CROSS SECTION ' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LO1),RWO(LO2),RWO(LS),EPS,NP,IP) C* Epithermal transport consistency check LZ1= L1+2*NG2 LZ2= L2+2*NG2 DO 418 I=1,NG12 RWO(LO1-1+I)=RWO(LO1-1+I)+RWO(LZ1-1+NG12+I) RWO(LO2-1+I)=RWO(LO2-1+I)+RWO(LZ2-1+NG12+I) 418 CONTINUE COMM = ' REF : TRANSP. VS. SCATT.MTX CONSISTENCY' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LZ1),RWO(LO1),RWO(LS),EPZ,NP,IP) COMM = ' COMP: TRANSP. VS. SCATT.MTX CONSISTENCY' CALL XSDIFF(LOU,LER 1 , 0,NG12,COMM,RWO(LZ2),RWO(LO2),RWO(LS),EPZ,NP,IP) C* Save Lambda*SigP for RI intercomparison LG1=L1 LG2=LG1+NG2 DO 420 I=1,NG2 RWO(LG1-1+I)=RWO(L1-1+I)*RWO(L1-1+NG2*3+NG12*2+I) RWO(LG2-1+I)=RWO(L2-1+I)*RWO(L2-1+NG2*3+NG12*2+I) 420 CONTINUE LS =LG2+NG2 LS0=LS C* C* Thermal cross sections processing C* 430 KLN=2*NG3 KLF=2*NG3 LS =LS0 C* Read temperature list for thermal cross section tables READ(LF1,702) (RWO(LS-1+I),I=1,NT1) LT1=LS LS =LS+NT1 READ(LF2,702) (RWO(LS-1+I),I=1,NT2) LT2=LS LS =LT2+NT2 L1 =LS C* Read the reference thermal data for one temperature IT2=0 IT1=0 440 IF(IT1.GE.NT1) GO TO 450 IT1=IT1+1 LS =L1 READ(LF1,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF1.GE.2) READ(LF1,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ(LF1,705) KL1,(RWO(LS-1+I),I=1,KL1) LS =LS+KL1+NG3 L2 =LS C* Read the compared thermal data for one temperature 450 IF(IT2.LT.NT2) GO TO 452 IF(IT1.LT.NT1) GO TO 440 GO TO 500 452 IT2=IT2+1 LS =L2 READ(LF2,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF2.GE.2) READ(LF2,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ(LF2,705) KL2,(RWO(LS-1+I),I=1,KL2) LS =LS+KL2+NG3 C* Compare data only if temp. differ by less than 10 degrees 460 IF(ABS( RWO(LT1-1+IT1) - RWO(LT2-1+IT2) ) .LT. 10.) GO TO 470 IF( RWO(LT1-1+IT1).GT.RWO(LT2-1+IT2) ) GO TO 450 C* - No provision is made for fine temperature mesh reference data GO TO 440 C* Begin comparison of thermal cross sections 470 IF(EPS.LT.0) GO TO 440 WRITE(LOU,699) FF WRITE(LOU,609) RWO(LT1-1+IT1) IP =4 LX1=L1 LX2=L2 COMM = ' TRANSPORT CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG3 LX2=LX2+NG3 COMM = ' ABSORPTION CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) LX1=LX1+NG3 LX2=LX2+NG3 C* Neutron fission yield and fission cross section IF(NF1.LT.2 .OR. NF2.LT.2) GO TO 472 COMM = ' NEUTRON FISSION YIELD ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) 472 LX1=LX1+NG3 LX2=LX2+NG3 IF(NF1.LT.2 .OR. NF2.LT.2) GO TO 474 COMM = ' FISSION CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LX1),RWO(LX2),RWO(LS),EPS,NP,IP) 474 LX1=LX1+NG3 LX2=LX2+NG3 C* Thermal scattering matrix (compare outscatter & check transport) LO1=LX1+KL1 LO2=LX2+KL2 CALL XSMSUM(NG3,1,RWO(LX1),RWO(LO1)) CALL XSMSUM(NG3,1,RWO(LX2),RWO(LO2)) COMM = ' OUTSCATTERING CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LO1),RWO(LO2),RWO(LS),EPS,NP,IP) CALL XSMSUM(NG3,0,RWO(LX1),RWO(LO1)) CALL XSMSUM(NG3,0,RWO(LX2),RWO(LO2)) COMM = ' TOTAL SCATTERING CROSS SECTION ' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(LO1),RWO(LO2),RWO(LS),EPS,NP,IP) C* Thermal transport consistency check DO 478 I=1,NG3 RWO(LO1-1+I)=RWO(LO1-1+I)+RWO(L1 -1+NG3+I) RWO(LO2-1+I)=RWO(LO2-1+I)+RWO(L2 -1+NG3+I) 478 CONTINUE COMM = ' REF : TRANSP. VS. SCATT.MTX CONSISTENCY' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(L1 ),RWO(LO1),RWO(LS),EPZ,NP,IP) COMM = ' COMP: TRANSP. VS. SCATT.MTX CONSISTENCY' CALL XSDIFF(LOU,LER 1 ,NG12,NG3,COMM,RWO(L2 ),RWO(LO2),RWO(LS),EPZ,NP,IP) GO TO 440 C* Rrocess resonance integral data if any 500 IF(NF1.LT.1 .OR. NF2.LT.1) GO TO 600 IF(NF1.NE.NF2) GO TO 596 C* Print the header to output IF(EPR.LT.0) GO TO 502 WRITE(LOU,614) FF IP =4 502 JRG=NG1 NR =0 510 LFF=LF1 LS =LS0 C* Read the filemark and check for labels 512 READ (LFF,701,ERR=512,END=592) JTERM,N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 514 IF(LFF.EQ.LF1) JRG=JRG+1 IF(N1.NE.0 .AND. N1.NE.5) GO TO 594 IF(N2.NE.0) JRG=N2 IF(N3.NE.0) NG1=N3 IF(N4.NE.0) NG2=N4 READ (LFF,701,ERR=591,END=592) NST IR =0 C* Read the data 514 READ (LFF,706) AM2,NT2,NS2 LL =NT2+NS2+NT2*NS2 READ (LFF,702) (RWO(LS-1+J),J=1,LL) L2 =LS LS =LS +LL IF(LFF.EQ.LF2) GO TO 520 C* Prepare to read the second data set LFF=LF2 AM1=AM2 NT1=NT2 NS1=NS2 L1 =L2 GO TO 512 520 CONTINUE C* For a group both data sets read - do the comparison IR =MIN(IR+1,3) NR =MAX(IR,NR) LRG=JRG-NG1 WRITE(COMM,612) RITY(IR),AM1,AM2 CALL RIDIFF(LOU,LER,COMM,JRG 1 ,NS1,NT1,RWO(L1),RWO(L1+NT1),RWO(L1+NT1+NS1),RWO(LG1-1+LRG) 2 ,NS2,NT2,RWO(L2),RWO(L2+NT2),RWO(L2+NT2+NS2),RWO(LG2-1+LRG) 3 ,RWO(LS),EPR,NP,IP) C* Check if more groups to be processed, else finish IF(JRG.LT.NG1+NG2 .OR. IR.LT.NR) GO TO 510 GO TO 600 C* Error traps 591 WRITE(LTT,611) ' COMPMT Error reading res. data on unit ',LFF GO TO 600 592 WRITE(LTT,611) ' COMPMT E-O-F reading res. data on unit ',LFF GO TO 600 594 WRITE(LTT,611) ' COMPMT incorrect file mark label found ',N1 GO TO 600 596 WRITE(LTT,611) ' COMPMT incompatible NF for RI compar. ',NF1,NF2 C* All data compared 600 CLOSE (UNIT=LF1) CLOSE (UNIT=LF2) CLOSE (UNIT=LOU) RETURN 606 FORMAT(/' Only differences in excess of',F6.1,' % are displayed') 607 FORMAT(// ' WILLIE - DATA COMPARISON UTILITY'/ 1 ' ================================'// 1 10X,' Mat.ID At.Wt. Z Nf Nt Nr Filename') 608 FORMAT( A10,':', I6, F9.4, I4,I3,I3,I3,1X,A40) 609 FORMAT(' TEMPERATURE',F8.1,' K'/' =====================') 611 FORMAT(A40,2I2) 612 FORMAT(1X,A10,' RI MAT.',F8.1,' AND',F8.1) 614 FORMAT(1X,A1//' RESONANCE INTEGRALS TABLES'/ 1 ,' ==========================') 691 FORMAT(2A40) 693 FORMAT(A40,I6) 694 FORMAT(BN,F20.0) 695 FORMAT(BN,I10) 699 FORMAT(1X,A1//) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,E15.8,2I6) END SUBROUTINE XSDIFF(LOU,LER,NG0,NGR,COM,XS1,XS2,RWO,EPS,NP,IP) C-Title : XSDIFF Subroutine C-Purpose: Compare cross sections CHARACTER*1 FF CHARACTER*40 COM DIMENSION XS1(*),XS2(*),RWO(*) C* Rel.diff.threshold below which data are considered identical DATA ERMX/1.E-5/ C* Form feed character FF =CHAR(12) C* Check if printout is requested IF(EPS.LT.0) RETURN C* Initialize ERM1=0. EAM1=0. ERM2=0. EAM2=0. JP =0 C* Calculate differences DO 60 I=1,NGR IG=NG0+I EA=XS2(I)-XS1(I) ER=0. IF(XS1(I).GT.0) ER= 100.*EA/XS1(I) IF(ER.GT. 999.) ER= 999. IF(ER.LT.-999.) ER=-999. IF(ER.LE.ERM1) GO TO 51 IRM1=IG ERM1=ER 51 IF(EA.LE.EAM1) GO TO 52 IAM1=IG EAM1=EA 52 IF(ER.GE.ERM2) GO TO 53 IRM2=IG ERM2=ER 53 IF(EA.GE.EAM2) GO TO 54 IAM2=IG EAM2=EA 54 RWO(I )=ER RWO(I+NGR)=EA IF(ABS(ER).GE.EPS) JP =JP+1 60 CONTINUE C* Check for space on output page KP =4 IF(ERM1.GT.ERMX .OR. -EAM2.GT.ERMX) KP=KP+3 IF(JP.GT.0) KP=KP+2+JP IF(IP+KP.LT.NP) GO TO 62 WRITE(LOU,99) FF IP =0 62 IP =IP+KP C* Write header WRITE(LOU,90) COM,NG0+1,NG0+NGR C* Check if there ara differences to be printed IF(JP.LT.1) GO TO 82 WRITE(LOU,91) DO 80 I=1,NGR IG=NG0+I ER=RWO(I) EA=RWO(I+NGR) X1=XS1(I) X2=XS2(I) IF(X1.LT.1E-9) X1=0. IF(X2.LT.1E-9) X2=0. IF(ABS(ER).GE.EPS) WRITE(LOU,92) IG,X1,X2,ER,EA 80 CONTINUE C* Print the extreme defferences 82 IF(ERM1.GT.ERMX .OR. -ERM2.GT.ERMX) GO TO 84 WRITE(LOU,96) C WRITE(LER,98) RETURN 84 WRITE(LOU,93) WRITE(LOU,94) ' Max.',IRM1,ERM1,IAM1,EAM1 WRITE(LOU,94) ' Min.',IRM2,ERM2,IAM2,EAM2 IRM=IRM1 IF(-ERM2.GT.ERM1) IRM=IRM2 ERM=ERM1 IF(-ERM2.GT.ERM1) ERM=-ERM2 C WRITE(LER,97) COM,IRM,ERM IF(ABS(ERM).GT.EPS) WRITE(LER,97) COM,IRM,ERM RETURN 90 FORMAT(/A40,' Group',I4,' to',I4/ 1 ,' ========================================================') 91 FORMAT(' Gr. Ref.[b] Comp.[b] Erel.[%] Eabs.[b]'/ 4 ,' --------------------------------------------') 92 FORMAT(1P,I5,2E10.3E1,0P,2F10.3) 93 FORMAT(' Extrema: Gr. Erel.[%] Gr. Eabs.[b]'/ 3 ,' ----------------------------------') 94 FORMAT(6X,A5,I4,F10.3,I6,F10.4) 96 FORMAT(' Data sets identical') 97 FORMAT(A40,' Gr',I3,' Emx',F10.3,' %') 98 FORMAT(A40,' Data identical') 99 FORMAT(1X,A1//) END SUBROUTINE RIDIFF(LOU,LER,COM,JRG 1 ,NS1,NT1,TM1,SG1,RI1,PL1 2 ,NS2,NT2,TM2,SG2,RI2,PL2 ,RWO,EPR,NP,IP) C-Title : RIDIFF Subroutine C-Purpose: Compare Resonance Integral data sets CHARACTER*40 COM CHARACTER*8 DASH, CH8, REC(10,6) CHARACTER*1 FF DIMENSION TM1(*),SG1(*),RI1(NS1,NT1) 1 ,TM2(*),SG2(*),RI2(NS2,NT2) ,RWO(*) DATA DASH /'--------'/ 1 ,REC(1,1) /' Sig0 '/ 1 ,REC(1,2) /' -------'/ 1 ,REC(1,3) /' Ref. '/ 1 ,REC(1,4) /' Comp. '/ 1 ,REC(1,5) /' Er [%] '/ 1 ,REC(1,6) /' Ea [b] '/ C* Rel.diff.threshold below which data are considered identical DATA ERMX/1.E-5/ C* Form feed character FF =CHAR(12) C* Check if printout is requested IF(EPR.LT.0) GO TO 600 EPS=0.01 C* Convert WIMS background x-sect to Bondarenko Sig0 DO 110 I=1,NS1 SG1(I)=SG1(I)-PL1 110 CONTINUE DO 120 I=1,NS2 SG2(I)=SG2(I)-PL2 120 CONTINUE C* Loop over temperatures - select temperature indices JT1=0 130 JT1=JT1+1 IF(JT1.GT.NT1) GO TO 600 JT2=0 131 JT2=JT2+1 JTI=MIN(NT2,JT2+1) IF(TM2(JTI).LT.TM1(JT1) .AND. JTI.LT.NT2) GO TO 131 C* Check for matching data at this temperature IF(ABS(TM1(JT1)-TM2(JT2)).GT.10. .AND. JT2.EQ.JTI) GO TO 130 C* Define Sig0 mesh LSG=1 KSG=NS1+NS2 LR1=LSG+KSG LR2=LR1+KSG LRI=LR2+KSG LBL=LRI+KSG C* Sig0 mesh as the union of both data sets (if possible) 140 JNT=1 SGB=0. NSG=0 C* Skip the smallest Sig0 to avoid extrapolation I1 =0 142 I1 =I1+1 IF(SG1(I1) .LT. SG2(1)) GO TO 142 I2 =0 144 I2 =I2+1 IF(SG2(I2) .LT. SG1(1)) GO TO 144 C* Identify the data set from which the next entry is taken 160 IF(SG1(I1) .GT. SG2(I2)) GO TO 166 C* Case: Next Sig0 mesh point from reference set 162 IF(I1.GT.NS1) GO TO 166 IF(ABS(SG1(I1)-SGB).LT.EPS*SGB) GO TO 164 NSG=NSG+1 SGB=SG1(I1) RWO(LSG-1+NSG)=SGB 164 I1 = I1+1 IF(I1.LE.NS1) GO TO 160 C* Case: Next Sig0 mesh point from compared set 166 IF(I2.GT.NS2) GO TO 168 IF(ABS(SG2(I2)-SGB).LT.EPS*SGB) GO TO 167 NSG=NSG+1 SGB=SG2(I2) RWO(LSG-1+NSG)=SGB 167 I2 = I2+1 IF(I2.LE.NS2) GO TO 160 GO TO 162 C* Last entry is 1.E10 - eliminate higher entries 168 IF(RWO(LSG-1+NSG).LE.1.E10) GO TO 170 RWO(LSG-1+NSG)=1.E10 IF(RWO(LSG-2+NSG).LT.RWO(LSG-1+NSG)) GO TO 170 NSG=NSG-1 GO TO 168 C* Check if Sig0 mesh needs to be reduced 170 JM =NSG RM =1.E30 C* Locate the densest Sig0 mesh spacing DO 172 I=2,NSG RR =ALOG(RWO(LSG+I-1)/RWO(LSG+I-2)) IF(RR.GT.RM) GO TO 172 RM =RR JM =I 172 CONTINUE IF(NSG.LE.9 .AND. (RR-1.).GT.EPS) GO TO 200 C* Reduce the Sig0 mesh by eliminating smallest Sig0 intervals J1 =MAX( 1,JM-2) J2 =MIN(NSG,JM+1) R1 =ALOG(RWO(LSG-1+JM)/RWO(LSG-1+J1 )) R2 =ALOG(RWO(LSG-1+J2)/RWO(LSG-1+JM-1)) C* Case: Eliminate lower point of the narrow Sig0 interval IF(R1.LT.R2 .AND. JM.NE. 2) JD=JM C* Case: Eliminate upper point of the narrow Sig0 interval IF(R1.GE.R2 .AND. JM.NE.NSG) JD=JM+1 C* Shift the Sig0 values DO 178 I=JD,NSG RWO(LSG+I-2)=RWO(LSG+I-1) 178 CONTINUE NSG=NSG-1 GO TO 170 C* Interpolate resonance integrals to a common grid 200 IF(JNT.LE.0) GO TO 300 C* Interpolate in the SigB domain (rather than Sig0) LS1=LBL DO 220 I=1,NSG RWO(LS1-1+I)=RWO(LSG-1+I)+PL1 220 CONTINUE CALL INTERI(NS1,NSG,SG1,RWO(LS1),RI1(1,JT1),RWO(LR1),INT) DO 240 I=1,NSG RWO(LS1-1+I)=RWO(LSG-1+I)+PL2 240 CONTINUE CALL INTERI(NS2,NSG,SG2,RWO(LS1),RI2(1,JT2),RWO(LR2),INT) CALL INTERI(NS2,NSG,SG2,RWO(LS1),RI2(1,JTI),RWO(LRI),INT) 300 IF(ABS(TM1(JT1)-TM2(JT2)).LT.10.) GO TO 400 C* Interpolate second data set for temperature 334 GG= (SQRT(TM1(JT1))-SQRT(TM2(JT2))) 1 / (SQRT(TM2(JTI))-SQRT(TM2(JT2))) DO 336 I=1,NSG RWO(LR2-1+I)=RWO(LR2-1+I) + GG*(RWO(LRI-1+I)-RWO(LR2-1+I)) 336 CONTINUE TM2(JT2)=TM1(JT1) C* Calculate the differences 400 CONTINUE IRM1=1 IAM1=1 IRM2=1 IAM2=1 ERM1=0. EAM1=0. ERM2=0. EAM2=0. DO 460 I=1,NSG EA = RWO(LR2-1+I)-RWO(LR1-1+I) ER = 0. IF(RWO(LR1-1+I).GT.0) ER = 100.*EA/RWO(LR1-1+I) IF(ER.GT. 999.) ER= 999. IF(ER.LT.-999.) ER=-999. IF(ER.LE.ERM1) GO TO 451 IRM1=I ERM1=ER 451 IF(EA.LE.EAM1) GO TO 452 IAM1=I EAM1=EA 452 IF(ER.GE.ERM2) GO TO 453 IRM2=I ERM2=ER 453 IF(EA.GE.EAM2) GO TO 454 IAM2=I EAM2=EA 454 CALL CH8PCK(REC(I+1,1),RWO(LSG-1+I)) REC(I+1,2) = DASH CALL CH8PCK(REC(I+1,3),RWO(LR1-1+I)) CALL CH8PCK(REC(I+1,4),RWO(LR2-1+I)) WRITE(REC(I+1,5),692) ER WRITE(REC(I+1,6),694) EA 460 CONTINUE C* Print the results JP=3 + 4 IF(ERM1.GT.EPR .OR. -ERM2.GT.EPR) JP=JP+7 IF(IP+JP.LT.NP) GO TO 461 WRITE(LOU,609) FF IP=0 461 IP=IP+JP WRITE(LOU,601) COM,TM1(JT1),JRG IF(ERM1.GT.ERMX .OR. -ERM2.GT.ERMX) GO TO 464 WRITE(LOU,696) C WRITE(LER,698) GO TO 130 464 IF(ERM1.LE.EPR .AND. -ERM2.LE.EPR) GO TO 470 NSG1=NSG+1 DO 466 I=1,6 WRITE(LOU,690) (REC(J,I),J=1,NSG1) 466 CONTINUE WRITE(LOU,690) 470 WRITE(LOU,603) WRITE(LOU,604) ' Max.',REC(IRM1+1,1),ERM1,REC(IAM1+1,1),EAM1 WRITE(LOU,604) ' Min.',REC(IRM2+1,1),ERM2,REC(IAM2+1,1),EAM2 CH8=REC(IRM1+1,1) IF(-ERM2.GT.ERM1) CH8=REC(IRM2+1,1) ERM=ERM1 IF(-ERM2.GT.ERM1) ERM=-ERM2 C WRITE(LER,697) COM,TM1(JT1),JRG,CH8,ERM IF(ABS(ERM).GT.EPR) WRITE(LER,697) COM,TM1(JT1),JRG,CH8,ERM GO TO 130 600 RETURN 601 FORMAT(/A40,' TEMP.',F6.1,' K GROUP',I3/ 1 ' ================================', 1 '================================') 603 FORMAT(' Extrema: Sig0 Erel.[%] Sig0 Eabs.[b]'/ 3 ,' ------------------------------------------') 604 FORMAT(6X,A5,A8,F10.3,2X,A8,F10.4) 609 FORMAT(1X,A1//) 690 FORMAT(10A8) 692 FORMAT(F8.2) 694 FORMAT(F8.4) 696 FORMAT(' Data sets identical') 697 FORMAT(A14,' Temp',F6.1,' K Gr',I3,' Sg0',A8,' Emx',F10.3,' %') 698 FORMAT(A14,' Temp',F6.1,' K Gr',I3,' Data identical') END SUBROUTINE INTERI(NA,NB,SA,SB,RA,RB,INT) C-Title : INTERI Subroutine C-Purpose: Interpolate a table of resonance integrals C-Description: C-D A set of NA input resonance integrals RA on Sig0 mesh SA C-D in ascending order is given. The last point is assumed at C-D infinite dilution. Resonance integrals are interpolated to C-D NB values RB at meshpoints SB. Flag INT>0 forces Root(Sig0) C-D interpolation, which is equivalent to the original WIMS-D/4 C-D method. C-D DIMENSION SA(*),SB(*),RA(*),RB(*) C* Define the RI at infinite dilution RI=RA(NA) C* Implied loop over output mesh IB=0 10 IB=IB+1 IF(IB.GT.NB) RETURN C* Locate the correct interval on the input mesh J2=1 20 J1=J2 J2=J1+1 IF(SA(J2).LT.SB(IB) .AND. J2.LT.NA-1) GO TO 20 C* Perform interpolation 30 CALL RINTER(SA(J1),SA(J2),SB(IB),RI,RA(J1),RA(J2),RB(IB),INT) GO TO 10 END SUBROUTINE RINTER(S1,S2,SS,RI,R1,R2,RR,INT) C-Title : RINTER Subroutine C-Purpose: Interpolate resonance integrals using Segev's formula C-Description: C-D Given the resonance integral at infinite dilition RI and C-D two resonance integrals R1,R2 at the corresponding Sig0 C-D values S1,S2, calculate the interpolated resonance integral C-D RR at Sig0 value SS. C-D If INT>0 then the Root(Sig0) interpolation is used, C-D except in the last interval extending to infinity where C-D the 1/Sig0 interpolation scheme is applied. This is C-D equivalent to the original WIMS-D/4 interpolation. This C-D interpolation rule is also forced automatically if the C-D resonance integrals are not monotonic ascending or when C-D the upper interval boundary S2 liev above SG0. Otherwise C-D the Segeev interpolation is used, which requires the C-D root WW of a transcendental equation in an interval W1,W2 C-D which encloses the root. Compared to the original Segev's C-D work, the differences in the present algorithm are a better C-D protection against numerical failure and a new root-finding C-D algorithm which is a combination of the bisection method C-D (HH=1) and the method of false position (HH=0), controlled C-D according to the success of the root estimate in the C-D previous iteration. C-D NOTE: S2 must be smaller than the assumed Sig0 at inf.dil. C-D (R20 80 IF(SS.LE.SB) GO TO 82 C* Special treatment for 1/Sig0 interpolation towards infinite dilution RA=RB RB=RI SA=SB SB=DIL IF(SB/SA.LT.100.) SB=SA*100. QA=1/SA QB=1/SB QQ=1/SS GO TO 84 C* Root(Sig0) interpolation 82 QA=SQRT(SA) QB=SQRT(SB) QQ=SQRT(SS) 84 RR=RA + (RB-RA)*(QQ-QA)/(QB-QA) RETURN END SUBROUTINE RESPOL(S1,S2,S,RI,R1,R2,R) C-Title : RESPOL Subroutine C-Purpose: Interpolate resonance integrals using Segev's formula C-Description: C-D Given the resonance integral at infinite dilition RI and C-D two resonance integrals R1,R2 at the corresponding Sig0 C-D values S1,S2, interpolate the resonance integral R to a C-D Sig0 value S C-D C-Author : Ales Holubar, Nucl.Res.Inst., Rez, Czechoslovakia (1990) C-Refern.: Segev, Nucl.Sci.Eng. 79, 113 (1981). C-Version: 93/2 A.Trkov, use Log interpolation if R1.ge.R2 C-V 94/8 A.Trkov, guard against small WH/WL interval. C-V SUPERSEEDED BY RINTER C- DATA EPS/1.0E-4/ IF(S1.LT.S2) GO TO 1 C* Swap entries if not in ascending order S0=S2 R0=R2 S2=S1 R2=R1 S1=S0 R1=R0 C* Test if RI are monotonic to avoid numerical problems 1 IF((R2-R1)/RI.LT.EPS) GO TO 6 IF((RI-R2)/RI.LT.EPS) GO TO 6 C* Initialize bounds for the root W G1=RI/R1 G2=RI/R2 RW=S1/S2 A1=ALOG(G1) A2=ALOG(G2) T =ALOG(G1)/ALOG(G2) D =ALOG(G2/G1) WL=ALOG(RW*T)/D WH=ALOG(RW)/D C* Iterate for the loop (stop on small diff. or function value 2 W =(WL+WH)*0.5 IF(ABS((WH-WL)/W).LT.1.0E-6) GO TO 4 G = G2**W-1.-RW*(G1**W-1.) IF(ABS(G) .LT. 0.01) GO TO 4 C* Reduce the bounds by bisection and re-iterate IF(G.LT.0.) GO TO 3 WL=W GO TO 2 3 WH=W GO TO 2 C* Iterations complete - evaluate the resonance integral 4 GW=(G1**W-1.) ET= GW*S1 P = 1./W R = RI*(S/(S+ET))**P RETURN C* Log interpolation if non-monotonic res.int. tables are given 6 R =R1+(R2-R1)*ALOG(S/S1)/ALOG(S2/S1) RETURN END SUBROUTINE CH8PCK(CH,R) C-Title : CH8PCK subroutine C-Purpose: Write a real number R into CH string, 8 characters long C-Author : A.Trkov, Institute Jozef Stefan, Ljubljana, Slovenia (1990) CHARACTER*8 CH IF( R .NE. 0) GO TO 10 CH =' 0' RETURN 10 AR =ABS(R) IF( AR .GE. 1.E-4 .AND. AR .LT. 1.E6) GO TO 20 WRITE(CH,11) R RETURN 20 IF( AR .GT. 9.99999) GO TO 30 WRITE(CH,21) R IF( CH(8:8) .NE. '0' ) RETURN 30 IF( AR .GT. 99.9999) GO TO 40 WRITE(CH,31) R IF( CH(8:8) .NE. '0' ) RETURN 40 IF( AR .GT. 999.999) GO TO 50 WRITE(CH,41) R IF( CH(8:8) .NE. '0' ) RETURN 50 IF( AR .GT. 9999.99) GO TO 60 WRITE(CH,51) R IF( CH(8:8) .NE. '0' ) RETURN 60 IF( AR .GT. 99999.9) GO TO 70 WRITE(CH,61) R IF( CH(8:8) .NE. '0' ) RETURN 70 IR=R+0.5 WRITE(CH,71) IR RETURN 11 FORMAT(1PE8.1) 21 FORMAT(F8.5) 31 FORMAT(F8.4) 41 FORMAT(F8.3) 51 FORMAT(F8.2) 61 FORMAT(F8.1) 71 FORMAT(I8) END *DECK PTABMT SUBROUTINE PTABMT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : PTABMT Subroutine C-Purpose: Tabulate nuclide data for plotting purposes C-Description: C-D To use high quality plotting software such as PLOTTAB, the C-D cross sections and resonance integrals must be tabulated on C-D separate files. C-D First the filename containing the group energy boundaries C-D is requested. It may be a WIMS library, a WILLIE file with C-D labelled filemarks or another list file where 70 (default) C-D group boundaries are read, 5 per line, 15 colums each (as C-D in the original WIMS library). C-D Next the list of files to be processed is entered, terminated C-D by a blank entry. The files may be WILLIE files produced by C-D the COPY option or otherwise (for example, from the WIMSR C-D module of the NJOY code). C-D The possibility is offered to define special labels which C-D appear at the beginning of each output set. The labes is C-D specific file containing the data set and may be up to 40 C-D characters long. The reaction type label is appended C-D automatically after column 40. C-D Reaction cross section to be plotted is selected by entering C-D the appropriate key number. A negative entry implies end of C-D processing. The data for the selected reaction will be written C-D on the output file from all input files. The procedure may be C-D repeated for as many reactions as desired. C-D The default output filename is displayed. It may be redefined C-D by entering a new filename, otherwise the data continue to be C-D written on the same file. C-D C- PARAMETER (MXFL=10) CHARACTER*10 IDNT(12) CHARACTER*40 BLNK,FLNM,FLNO,FLNI(MXFL) CHARACTER*80 COMM(MXFL) DIMENSION RWO(*),IWO(*), RMT(MXFL) COMMON /TERWOR/ ITERM,ALDO(2) SAVE FLNO DATA BLNK /' '/ 1 ,FLNO /'PTABMT.CUR'/ C* Number and selectable reaction types DATA MSEL,IDNT / 12 1,'POTENTIAL ','SLOW.DN.PW','TRANSPORT ','ABSORPTION','FISS.YIELD' 1,'FISSION ','OUTSCATT. ','SCATTERING','DIFF.SCAT.','ABSORPT.RI' 2,'FISSION RI','FISS.SPECT'/ DATA RMT/ MXFL*0.0 / C* Default group system (No.of fast,resonance,thermal groups) DATA NG1,NG2,NG3/ 14, 13, 43 / C* Logical file units for source library and output file DATA LF1,LF2 / 1, -2 / C* Default temperatures [K] for cross sect. and res.int. TMX=300. TMR=300. LBL=1 DO 11 I=1,MXFL RMT(I)=0 11 CONTINUE C* Define the filenames GO TO 14 12 WRITE(LTT,691) ' WIMS library of WILLIE file with labell' 1 ,'ed filemarks must be entered. ' 14 WRITE(LTT,691) '$Enter group boundaries filename : ' READ (LKB,691) FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLNM,STATUS='OLD',ERR=12) C***** EXPORT ***** C* Interpret the header READ (LF1,701,ERR=16) JTERM,N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 18 IF(N1.NE.0 .AND. N1.NE.1) GO TO 12 C* Case: WILLIE file with labelled filemarks NG1=N2 NG2=N3 NG3=N4 GO TO 20 C* Case: file without headers 16 REWIND LF1 GO TO 20 C* Case: WIMS library 18 NEL=JTERM NG =N1 NG0=N2 NG1=N3 NG2=N4 READ (LF1,701,ERR=12) NG3,NNFD,NFPD READ (LF1,701,ERR=12) (IWO(I),I=1,NEL) C* Read the group boundaries (in descending energy) 20 NG =NG1+NG2+NG3 NGB=NG+1 READ (LF1,702,ERR=12) (RWO(LBL-1+J),J=1,NGB) LGB=LBL LGC=LGB+NGB LBL=LBL+NGB+NG2 CLOSE(UNIT=LF1) C* Define the data filenames 30 WRITE(LTT,691) ' Define the cross section data files : ' NFL=0 32 WRITE(LTT,691) '$ : ' READ (LKB,691) FLNM IF(FLNM.EQ.BLNK) GO TO 34 NFL=NFL+1 FLNI(NFL)=FLNM GO TO 32 C* Enter the labels for each data set 34 FLNM= ' LABEL ' DO 36 JFL=1,NFL C* One blank entry implies default setting for the remaining files IF(FLNM.EQ.BLNK) GO TO 35 WRITE(LTT,691) ' File : ' 1,FLNI(JFL) WRITE(LTT,691) '$Enter the label (default="filename") : ' READ (LKB,691) FLNM COMM(JFL)(1:40)=FLNM IF(FLNM.NE.BLNK) GO TO 36 35 COMM(JFL)(1:40)=FLNI(JFL) 36 CONTINUE C* Define the parameter to be tabulated 40 WRITE(LTT,691) ' Select from the following reactions : ' WRITE(LTT,691) ' 1 - potential cross section ' WRITE(LTT,691) ' 2 - slowing down power ' WRITE(LTT,691) ' 3 - transport cross section ' WRITE(LTT,691) ' 4 - absorption cross section ' WRITE(LTT,691) ' 5 - neutron fission yield ' WRITE(LTT,691) ' 6 - fission cross section ' WRITE(LTT,691) ' 7 - outscattering cross section ' WRITE(LTT,691) ' 8 - total scattering cross section ' WRITE(LTT,691) ' 9 - group scattering matrix ' WRITE(LTT,691) ' 10 - absorption resonance integrals ' WRITE(LTT,691) ' 11 - fission resonance integrals ' WRITE(LTT,691) ' 12 - fission spectrum ' WRITE(LTT,691) '$ -1 - to finish : ' READ (LKB,501) ISEL IF(ISEL.LT.0) GO TO 500 IF(ISEL.LT.1 .OR. ISEL.GT.MSEL) GO TO 40 C* Define the output filename WRITE(LTT,691) ' Default output table filename : ',FLNO WRITE(LTT,691) '$ Enter new name to redefine : ' READ (LKB,691) FLNM IF(FLNM.EQ.BLNK) FLNM=FLNO IF(FLNM.EQ.FLNO .AND. LF2.GT.0) GO TO 48 FLNO=FLNM IF(LF2.GT.0) CLOSE(UNIT=LF2) LF2=IABS(LF2) C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLNM,STATUS='NEW',CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF2,FILE=FLNM,STATUS='UNKNOWN') C***** EXPORT ***** C* C* Begin processing the data files 48 JFL=1 GO TO 100 C* Error condition reading a file 90 WRITE(LTT,691) ' Error reading file ' 1 , FLNI(JFL) WRITE(LTT,691) '$Enter new filename or blank to skip : ' READ (LKB,691) FLNI(JFL) C* Loop over all the files 100 JEL=0 NEL=1 MAT=RMT(JFL)+0.01 IF(FLNI(JFL) .EQ. BLNK) GO TO 400 C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLNI(JFL),STATUS='OLD',ERR=90,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLNI(JFL),STATUS='OLD',ERR=90) C***** EXPORT ***** C* Identify the file format READ (LF1,701,ERR=160) JTERM,N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 180 110 IF(N1.EQ.0) GO TO 112 IF(N1.EQ.4) GO TO 150 IF(N1.EQ.2) GO TO 120 112 CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) IF(JTERM.LT.0) GO TO 90 GO TO 110 C* Case: WILLIE fission spectrum file with labelled filemarks 120 IF(ISEL.NE.12) GO TO 112 IF(N2.GT.0) NG0=N2 122 IF(NFL.GT.1) THEN COMM(JFL)=IDNT(ISEL)//' ' ELSE COMM(JFL)( 1:10)=IDNT(ISEL) END IF CALL PTABSP(RWO(LBL),COMM(JFL),RWO(LGB),LF1,LF2,NG0) GO TO 400 C* Case: WILLIE cross section file with labelled filemarks 150 IF(N2.GT.0 .AND. N2.NE.NG1) GO TO 90 IF(N3.GT.0 .AND. N3.NE.NG2) GO TO 90 IF(N4.GT.0 .AND. N4.NE.NG3) GO TO 90 GO TO 200 C* Case: file without headers 160 REWIND LF1 GO TO 200 C* Case: WIMS library 180 NEL=JTERM IF(N1.NE.NG )GO TO 90 IF(N3.NE.NG1)GO TO 90 IF(N4.NE.NG2)GO TO 90 READ (LF1,701,ERR=90) NG3,NNFD,NFPD READ (LF1,701,ERR=90) (IWO(I),I=1,NEL) C* Blank-read the group boundaries READ (LF1,702) (RWO(LBL-1+I),I=1,NG) C* Check if fission spectrum is to be processed IF(ISEL.EQ.12) GO TO 122 C* Select the material to be processed 185 IF(RMT(JFL).GT.0) GO TO 188 WRITE(LTT,691) '$ Enter the selected material ident. : ' READ (LKB,502,ERR=187) RMT(JFL) IF(RMT(JFL).EQ.0) GO TO 400 C* Check if material is on the list MAT=RMT(JFL)+0.001 DO 186 I=1,NEL IF(IWO(I).EQ.MAT) GO TO 188 186 CONTINUE WRITE(LTT,601) MAT 187 WRITE(LTT,602) (IWO(J),J=1,NEL) GO TO 185 188 CONTINUE C* Skip to the beginning of the nuclide files CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) C* C* Read the nuclide header card 200 READ (LF1,704) ID,AW,IZ,NF,NT,NR IF(MAT.EQ.0 .OR. MAT.EQ.ID) GO TO 220 C* Skip to the next file CALL SKIPWF(LF1,ITERM,JTERM,N1,N2,N3,N4) IER=109 IF(JTERM.LT.0) GO TO 590 GO TO 240 C* Select the temperature for thermal cross section data 220 IF(ISEL.GE.10) GO TO 230 221 WRITE(LTT,691) '$ Enter temperature for cross sect. : ' READ (LKB,502,ERR=221) TMX C* Tabulate the cross sections 230 IF(NFL.GT.1) THEN COMM(JFL)(41:80)=IDNT(ISEL)//' ' ELSE COMM(JFL)(41:80)=COMM(JFL)( 1:40) COMM(JFL)( 1:40)=IDNT(ISEL)//' ' END IF WRITE(COMM(JFL)(61:68),605) ID CALL PTABXS(RWO(LBL),IWO,COMM(JFL),RWO(LGB),RWO(LGC) 1 ,LF1,LF2,NG1,NG2,NG3,NF,NT,TMX,ISEL) C* Check the filemark IER=110 READ (LF1,701,ERR=590) JTERM IF(JTERM.NE.ITERM) GO TO 590 240 JEL=JEL+1 IF(JEL.LT.NEL) GO TO 200 C* C* Begin processing the resonance files 300 IF(ISEL.LT.10) GO TO 400 IF(JFL .GT.1) GO TO 320 302 WRITE(LTT,691) '$ Enter temperature for R.I. : ' READ (LKB,502,ERR=302) TM IF(TM.NE.0) TMR=TM WRITE(LTT,691) '$ Enter R.I. interp.( >0=SQRT(Sig0) ): ' READ (LKB,502,ERR=302) TRP INT=0 IF(TRP.GT.0) INT=1 WRITE(LTT,691) '$ Enter group (default=all groups) : ' READ (LKB,502,ERR=302) GR JGR=GR+0.1 C* Tabulate the resonance integrals 320 IF(NFL.GT.1) THEN COMM(JFL)(41:50)=IDNT(ISEL) ELSE COMM(JFL)( 1:10)=IDNT(ISEL) END IF CALL PTABRI(RWO(LBL),IWO,COMM(JFL),LF1,LF2 1 ,NG1,NG2,NG3,JGR,RWO(LGC),RMT(JFL),TMR,ISEL,INT) 400 CLOSE(UNIT=LF1) JFL=JFL+1 IF(JFL.LE.NFL) GO TO 100 GO TO 40 C* Processing completed (leave output unit open to allow append) 500 CLOSE(UNIT=LF1) RETURN C* Library read error 590 CALL ERRWAR(LTT,IER) CLOSE(UNIT=LF1) CLOSE(UNIT=LF2) RETURN 501 FORMAT(BN,I20) 502 FORMAT(BN,F20.0) 601 FORMAT(' WILLIE ERROR - Material',I6,' not found') 602 FORMAT(' Select from the list below:'/ (6X,10I6) ) 603 FORMAT(' Processing material',F8.1,' completed'/) 605 FORMAT('MAT',I5) 691 FORMAT(2A40) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,I15/E15.8,2I6/(5E15.8)) END SUBROUTINE SKIPWF(LUN,ITERM,JTERM,N1,N2,N3,N4) C-Title : SKIPWF Subroutine C-Purpose: Skip WIMS or WILLIE file till the next filemark 10 READ (LUN,92,ERR=10,END=90) JTERM,N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 10 RETURN C* EOF encountered before filemark 90 JTERM=-1 RETURN 92 FORMAT(5I15) END SUBROUTINE PTABSP(RWO,COM,GEB,LF1,LF2,NG0) C-Title : PTABSP Subroutine C-Purpose: Tabulate fission spectrum in PLOTTAB histogram format CHARACTER*80 COM DIMENSION RWO(*),GEB(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Read the spectrum READ (LF1,702) (RWO(I),I=1,NG0) C* Write the complete data set to output in ascending energy 80 WRITE(LF2,800) COM DO 82 I=1,NG0 WRITE(LF2,801) GEB(NG0+2-I),RWO(NG0+1-I) WRITE(LF2,801) GEB(NG0+1-I),RWO(NG0+1-I) 82 CONTINUE WRITE(LF2,800) RETURN 702 FORMAT(1P5E15.8) 800 FORMAT(A80) 801 FORMAT(1P,2E11.4) END SUBROUTINE PTABXS(RWO,IWO,COM,GEB,PLG,LF1,LF2 1 ,NG1,NG2,NG3,NF,NT,TMP,ISEL) C-Title : PTABXS Subroutine C-Purpose: Tabulate cross sections in PLOTTAB histogram format CHARACTER*80 COM DIMENSION RWO(*),IWO(*),GEB(*),PLG(*) COMMON /TERWOR/ ITERM,ALDO(2) NG12=NG1+NG2 KLN =2*NG12+4*NG2 LBL =1 L1 =LBL C* Read Epithermal cross sections READ (LF1,702) (RWO(LBL-1+I),I=1,KLN) C* Save SigP*Lambda DO 20 I=1,NG2 PLG(I)=RWO(LBL-1+I)*RWO(LBL+2*NG12+3*NG2-1+I) 20 CONTINUE LBL =LBL+KLN C* Read fission yield and fission cross section KLN =2*NG12 IF(NF.GT.1) READ (LF1,702) (RWO(LBL-1+I),I=1,KLN) LBL =LBL+KLN C* Read epithermal scattering matrix 40 READ (LF1,705) NDAT2,(RWO(LBL-1+I),I=1,NDAT2) LBL =LBL+NDAT2 C* Read temperature list for thermal cross sections READ (LF1,702) (RWO(LBL-1+I),I=1,NT) LT =LBL LBL =LBL+NT L2 =LBL C* Read thermal cross sections for all temperatures DO 60 K=1,NT IF(RWO(LT-1+K).LE.TMP) L2=LBL KLN =2*NG3 READ (LF1,702) (RWO(LBL-1+I),I=1,KLN) LBL =LBL+KLN IF(NF.GT.1) READ (LF1,702) (RWO(LBL-1+I),I=1,KLN) LBL =LBL+KLN 50 READ (LF1,705) NDAT3,(RWO(LBL-1+I),I=1,NDAT3) LBL =LBL+NDAT3 60 CONTINUE IF(ISEL.GT.9) RETURN C* Shift the selected epithermal data to beginning of work field NX =NG12 LB =NG1+NG2+NG3+1 IF(ISEL.LT.3) NX=NG2 IF(ISEL.LT.3) LB=NG12+1 IF(ISEL.EQ.1) LX=L1 IF(ISEL.EQ.2) LX=L1+NG2 IF(ISEL.EQ.3) LX=L1+NG2*2 IF(ISEL.EQ.4) LX=L1+NG2*2+NG12 IF(ISEL.EQ.5) LX=L1+NG2*4+NG12*2 IF(ISEL.EQ.6) LX=L1+NG2*4+NG12*3 IF(ISEL.GT.6) GO TO 72 C* Case: process epithermal cross sections DO 70 I=1,NX RWO(I)=RWO(LX-1+I) 70 CONTINUE GO TO 75 C* Case: process epithermal scattering matrix 72 LX =L1+NG2*4+NG12*4 ISS =8-ISEL IF(ISS.LT.0) STOP ' PTABMT MATRIX PRINTOUT OPTION NOT AVAILABLE' CALL XSMSUM(NG12,ISS,RWO(LX),RWO) C* Shift the selected thermal data to beginning of work field 75 IF(ISEL.LT.3) GO TO 80 IF(ISEL.EQ.3) LX=L2 IF(ISEL.EQ.4) LX=L2+NG3 IF(ISEL.EQ.5) LX=L2+NG3*2 IF(ISEL.EQ.6) LX=L2+NG3*3 IF(ISEL.GT.6) GO TO 77 C* Case: process thermal cross sections DO 76 I=1,NG3 RWO(I+NX)=RWO(LX-1+I) 76 CONTINUE GO TO 78 C* Case: process thermal scattering matrix 77 LX=L2+NG3*4 CALL XSMSUM(NG3,ISS,RWO(LX),RWO(NX+1)) 78 NX =NX+NG3 C* Write the complete data set to output in ascending energy 80 WRITE(LF2,800) COM IF(NF.LT.2 .AND. (ISEL.EQ.5 .OR. ISEL.EQ.6)) GO TO 84 IF(NF.NE.3 .AND. ISEL.EQ.11 ) GO TO 84 DO 82 I=1,NX WRITE(LF2,801) GEB(LB+1-I),RWO(NX+1-I) WRITE(LF2,801) GEB(LB -I),RWO(NX+1-I) 82 CONTINUE 84 WRITE(LF2,800) RETURN 701 FORMAT(5I15) 702 FORMAT(1P5E15.8) 705 FORMAT(1P,I15/(5E15.8)) 800 FORMAT(A80) 801 FORMAT(1P,2E11.4) END SUBROUTINE PTABRI(RWO,IWO,COM,LF1,LF2,NG1,NG2,NG3,JG 1 ,PLG,RMT,TMP,ISEL,INT) C-Title : PTABRI Subroutine C-Purpose: Tabulate Resonance Integrals in PLOTTAB format CHARACTER*80 COM DIMENSION RWO(*),IWO(*),PLG(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Define the number of Sig0 points to be generated DATA NSG /181/ LTM=1 LS0=1000+1 LR0=2000+1 C* Loop over resonance groups DO 90 KG=1,NG2 IG =NG1+KG NRT=0 C* Implied loop over all res.tables in a group 70 READ (LF1,701) M12 IF(M12.EQ.ITERM) GO TO 90 READ (LF1,706) RID,M1,M2,(RWO(LTM-1+I),I=1,M1), * (RWO(LS0-1+I),I=1,M2), * (RWO(LR0-1+I),I=1,M12) IF(RID.EQ.0.) GO TO 80 IF( JG.GT.0 .AND. JG.NE.IG ) GO TO 70 IF(RMT.NE.0.0 .AND. ABS(RID-RMT).GT.0.01) GO TO 70 NRT=NRT+1 IF(ISEL.EQ.10 .AND. NRT.NE.1) GO TO 70 IF(ISEL.EQ.11 .AND. NRT.NE.2) GO TO 70 C* Determine the SigB interpolation table LBL=LR0 +M12 LS1=LBL LR1=LS1+NSG DSG=EXP( ALOG(RWO(LS0+M2-1)/RWO(LS0)) / (NSG-1) ) SG0=RWO(LS0) DO 72 I=1,NSG RWO(LS1-1+I)=SG0 SG0 = SG0*DSG 72 CONTINUE CALL INTERI(M2,NSG,RWO(LS0),RWO(LS1),RWO(LR0),RWO(LR1),INT) C* Printout, convert WIMS SigB to Bondarenko Sif0 WRITE(COM(44:50),699) IG WRITE(LF2,692) COM DO 74 I=1,NSG WRITE(LF2,601) RWO(LS1-1+I)-PLG(KG),RWO(LR1-1+I) 74 CONTINUE WRITE(LF2,692) GO TO 70 C* Group resonance tables processed - check the file mark 80 READ (LF1,701,ERR=98) JTERM IF(JTERM.NE.ITERM) GO TO 98 90 CONTINUE RETURN C* Library read error 98 STOP 'PTABRI Read error' 601 FORMAT(1P,2E11.4) 692 FORMAT(A80) 699 FORMAT('.GR.',I3) 701 FORMAT(I15) 706 FORMAT(1P,E15.8,2I6/(5E15.8)) END *DECK SMTXTR SUBROUTINE SMTXTR(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : SMTXTR Subroutine C-Purpose: Correct the self-scattering term of the scatt.matrix C-Description: C-D A cross section file such as produced by WILLIE option COPY C-D is processed. It generously allows labelled, unlabelled C-D as well as no-filemark formats. Cross sections and scattering C-D matrices are read and checked, a correction to the scattering C-D matrix is performed and a new WILLIE cross section file for C-D a material is written. Any P1 scattering matrices are ignored. C-D CHARACTER*40 FLN1,FLN2,FLNM CHARACTER*80 REC DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) C* Default group structure (fast/resonance/thermal) DATA NG1,NG2,NG3/ 14, 13, 42 / C* Local file units (unit 2 flagged -ve before opening) DATA LF1,LF2,LOU / 1,-2, 3 / C* Define the source cross section filename WRITE(LTT,691) ' SOURCE CROSS SECTION DATA: ' GO TO 12 11 WRITE(LTT,691) ' ERROR reading file : ',FLN2 12 WRITE(LTT,691) '$ Enter the filename: ' READ (LKB,691) FLN2 C***** VAX ***** C OPEN(UNIT=LF1,FILE=FLN2,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN(UNIT=LF1,FILE=FLN2,STATUS='OLD',ERR=11) C***** EXPORT ***** C* Define the output cross section filename IF(LF2.GT.0) GO TO 18 LF2=-LF2 WRITE(LTT,691) ' MODIFIED CROSS SECTION DATA: ' WRITE(LTT,691) '$ Enter the filename: ' READ (LKB,691) FLNM IF(FLNM.NE.' ') FLN2=FLNM C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLN2,STATUS='NEW',CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF2,FILE=FLN2,STATUS='UNKNOWN') C***** EXPORT ***** C* Try to interpret the filemark record 18 REWIND LF2 READ (LF1,701,ERR=26) JTERM, NTYP,JG1,JG2,JG3 WRITE(LF2,701 ) JTERM, NTYP,JG1,JG2,JG3 IF(NTYP.EQ.0) GO TO 32 20 IF(NTYP - 4 ) 30,32,28 C* No cross section data found 28 WRITE(LTT,691) ' No cross section data on file : ',FLN2 GO TO 12 C* No labelled filemark present 26 REWIND LF1 GO TO 400 C* Inappropriate labelled filemark found - skip the data 30 READ (LF1,693) REC WRITE(LF2,693) REC IF(REC(1:15).NE.' 999999999') GO TO 30 READ (REC,701) JTERM, NTYP,JG1,JG2,JG3 GO TO 20 C* Labelled filemark for cross section data found 32 IF(JG1.GT.0) NG1=JG1 IF(JG2.GT.0) NG2=JG2 IF(JG3.GT.0) NG3=JG3 NG12=NG1 +NG2 NG = NG12+NG3 C* Begin processing the cross sections 400 READ (LF1,704) MT1,AW1,IZ1,NF1,NT1,NR1 WRITE(LF2,704) MT1,AW1,IZ1,NF1,NT1,NR1 C* C* Epithermal cross sections processing C* 402 KLN=4*NG2+2*NG12 KLF=2*NG12 LS =1 C* Read the source epithermal data L1 =LS READ (LF1,702) (RWO(LS-1+I),I=1,KLN) WRITE(LF2,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF1.GE.2) READ (LF1,702) (RWO(LS-1+I),I=1,KLF) IF(NF1.GE.2) WRITE(LF2,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ (LF1,705) KL1,(RWO(LS-1+I),I=1,KL1) LX1=LS LS =LS+KL1+NG12 LO1=LS C* Sum the scattering matrix elements CALL XSMSUM(NG12,0,RWO(LX1),RWO(LO1)) C* Calculate the necessary correction to the self-scattering term LZ1= L1+2*NG2 DO 418 I=1,NG12 RWO(LO1-1+I)=RWO(LO1-1+I)+RWO(LZ1-1+NG12+I) - RWO(LZ1-1+I) 418 CONTINUE C* Subtract the difference from the self-scattering term LSS=LX1 DO 420 I=1,NG12 NU =RWO(LSS )+0.1 NS =RWO(LSS+1)+0.1 RWO(LSS+1+NU)=RWO(LSS+1+NU)-RWO(LO1-1+I) LSS=LSS+NS+2 420 CONTINUE WRITE(LF2,705) KL1,(RWO(LX1-1+I),I=1,KL1) C* C* Thermal cross sections processing C* 430 KLN=2*NG3 KLF=2*NG3 LS =1 C* Read temperature list for thermal cross section tables READ (LF1,702) (RWO(LS-1+I),I=1,NT1) WRITE(LF2,702) (RWO(LS-1+I),I=1,NT1) LT1=LS LS =LS+NT1 L1 =LS C* Read the reference thermal data for one temperature IT1=0 440 IT1=IT1+1 IF(IT1.GT.NT1) GO TO 500 LS =L1 READ (LF1,702) (RWO(LS-1+I),I=1,KLN) WRITE(LF2,702) (RWO(LS-1+I),I=1,KLN) LS =LS+KLN IF(NF1.GE.2) READ (LF1,702) (RWO(LS-1+I),I=1,KLF) IF(NF1.GE.2) WRITE(LF2,702) (RWO(LS-1+I),I=1,KLF) LS =LS+KLF READ (LF1,705) KL1,(RWO(LS-1+I),I=1,KL1) LX1=LS LS =LS+KL1 LO1=LS C* Sum the scattering matrix elements CALL XSMSUM(NG3 ,0,RWO(LX1),RWO(LO1)) C* Calculate the necessary correction to the self-scattering term LZ1= L1 DO 448 I=1,NG3 RWO(LO1-1+I)=RWO(LO1-1+I)+RWO(LZ1-1+NG3 +I) - RWO(LZ1-1+I) 448 CONTINUE C* Subtract the difference from the self-scattering term LSS=LX1 DO 450 I=1,NG3 NU =RWO(LSS )+0.1 NS =RWO(LSS+1)+0.1 RWO(LSS+1+NU)=RWO(LSS+1+NU)-RWO(LO1-1+I) LSS=LSS+NS+2 450 CONTINUE WRITE(LF2,705) KL1,(RWO(LX1-1+I),I=1,KL1) GO TO 440 C* Rrocess resonance integral data if any 500 IF(NF1.LT.1 .OR. NF1.GT.3) GO TO 600 C* Read the filemark and check for labels JRG=NR1 510 READ (LF1,701,ERR=591,END=592) JTERM,N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 514 WRITE(LF2,701 ) JTERM,N1,N2,N3,N4 IF(N1.NE.0 .AND. N1.NE.5) GO TO 594 IF(N2.NE.0) JRG=N2 IF(N3.NE.0) NG1=N3 IF(N4.NE.0) NG2=N4 READ (LF1,701,ERR=591,END=592) JTERM JRG=JRG+1 C* Read the resonance data 514 WRITE(LF2,701 ) JTERM READ (LF1,706) AM1,NT1,NS1 READ (LF2,706) AM1,NT1,NS1 LL =NT1+NS1+NT1*NS1 READ (LF1,702) (RWO(LS-1+J),J=1,LL) WRITE(LF2,702) (RWO(LS-1+J),J=1,LL) C* Check if more groups to be processed, else finish IF(JRG.LT.NG1+NG2) GO TO 510 GO TO 600 C* Error traps 591 WRITE(LTT,611) ' COMPMT Error reading res. data on unit ',LF1 GO TO 600 592 WRITE(LTT,611) ' COMPMT E-O-F reading res. data on unit ',LF1 GO TO 600 594 WRITE(LTT,611) ' COMPMT incorrect file mark label found ',N1 C* All data compared 600 WRITE(LF2,701 ) ITERM CLOSE (UNIT=LF1) CLOSE (UNIT=LF2) LF2=-LF2 RETURN 607 FORMAT(// ' WILLIE - SCATTERING MATRIX TRANSPORT CORRECTION'/ 1 ' ==============================================='/) 611 FORMAT(A40,2I2) 691 FORMAT(2A40) 693 FORMAT(A80) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,E15.8,2I6) END *DECK INCLMT SUBROUTINE INCLMT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : INCLMT Subroutine C-Purpose: Add data for a material into the WIMS library C-Description: C-D A deta set extracted previously by the COPY option or C-D prepared otherwise can be inserted into the WIMS library C-D with a different nuclide identification number. C-D CHARACTER*75 REC , REN CHARACTER*40 BLNK, FLNM, FLWI, FLWO, FLNI DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) DATA BLNK /' '/ 1 ,FLWI /'WIMS.LIB'/ 1 ,FLWO /'WIMS.LIB'/ C* Default file units DATA LF1,LF2,LF3/ 1, 2, 3 / C* NDAT1=0 NFBU =0 KR =0 KF =0 KG0 =0 C* Define the filenames GO TO 12 11 WRITE(LTT,691) ' ERROR opening file : ',FLWI 12 WRITE(LTT,691) ' Default source formatted lib.filename: ',FLWI WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691,END=684) FLNM IF(FLNM.NE.BLNK) FLWI=FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11) C***** EXPORT ***** GO TO 14 13 WRITE(LTT,691) ' ERROR opening file : ',FLNI 14 WRITE(LTT,691) '$Enter nuclide data input filename : ' READ (LKB,691) FLNI IF(FLNI.EQ.BLNK) GO TO 683 C***** VAX ***** C OPEN (UNIT=LF2,FILE=FLNI,STATUS='OLD',ERR=13,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF2,FILE=FLNI,STATUS='OLD',ERR=13) C***** EXPORT ***** GO TO 16 15 WRITE(LTT,691) ' ERROR opening file : ',FLWO 16 WRITE(LTT,691) ' Default new formatted lib. filename: ',FLWO WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691,END=682) FLNM IF(FLNM.NE.BLNK) FLWO=FLNM C***** VAX ***** C OPEN (UNIT=LF3,FILE=FLWO,STATUS='NEW',ERR=15 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** IF(FLWO.EQ.FLWI) GO TO 15 OPEN (UNIT=LF3,FILE=FLWO,STATUS='UNKNOWN',ERR=15) C***** EXPORT ***** C* C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' Source WIMS-D library filename : ',FLWI WRITE(LTT,691) ' New WIMS-D library filename : ',FLWO WRITE(LTT,691) ' Cross sections filename : ',FLNI WRITE(LTT,691) C* Check the data to be inserted MAT=-1 LGB=0 LSP=0 LBU=0 LXS=0 LRI=0 LP1=0 LSC=1000 READ (LF2,701,ERR=102) JTERM, N1,N2,N3,N4 GO TO 110 C* If file mark not found, assume only cross sections given 102 REWIND LF2 N1=0 N2=0 N3=0 N4=0 GO TO 140 110 IF(N1.EQ.0 .OR. N1.GT.1) GO TO 120 C* - Group boundaries to be inserted IF(N2.GT.0) NG1=N2 IF(N3.GT.0) NG2=N3 IF(N4.GT.0) NG3=N4 NG =NG1+NG2+NG3 LGB=LSC LSC=LSC+NG+1 READ (LF2,702) (RWO(LGB+I),I=1,NG+1) READ (LF2,701) JTERM, N1,N2,N3,N4 120 IF(N1.EQ.0 .OR. N1.GT.2) GO TO 130 C* - Spectrum to be inserted IF(N2.GT.0) KG0=N2 LSP=LSC LSC=LSC+KG0+1 READ (LF2,702) (RWO(LSP+I),I=1,KG0) READ (LF2,701,END=141) JTERM, N1,N2,N3,N4 130 IF(N1.EQ.0 .OR. N1.GT.3) GO TO 140 C* - Burnup data to be inserted READ (LF2,701) NDAT1 LBU=LSC LSC=LSC+NDAT1/2 READ (LF2,703) (RWO(LBU+I),IWO(LBU+I),I=1,NDAT1/2) READ (LF2,701,END=141) JTERM, N1,N2,N3,N4 C* Save the NF flag for consistency checks IF(NDAT1.EQ.8) NFBU=IWO(LBU+4) C* - Cross section data to be inserted 140 IF(N1.NE.0 .AND. N1.NE.4) GO TO 141 READ (LF2,704,END=141) MAT,WW,KZ,KF,KT,KR C* Check NF in burnup data and cross sections consistency IF(KF.GT.1 .AND. NFBU.LT.0) THEN IF(KF.EQ.2 .AND. NFBU.NE.-2) THEN WRITE(LTT,691) ' WARNING - Absorption RI data present - ' 1 ,'Parameter NF should be -2 ' END IF IF(KF.EQ.3) THEN WRITE(LTT,691) ' ERROR - Fissile material with fission ' 1 ,'RI tables is used as a fission product ' STOP 'WILLIE ERROR in burnup data' END IF WRITE(LTT,691) ' WARNING - Fissile material is used as a' 1 ,' fission product ' WRITE(LTT,691) ' Fission cross section deleted' NFBU= 0 KF =-KF ELSE NFBU= KF END IF C* - Check for the material identifier WRITE(LTT,693) ' Default material ident. on file is : ',MAT WRITE(LTT,691) ' Enter blank to accept , ' WRITE(LTT,691) ' negative to quit , ' WRITE(LTT,691) '$ new value to redefine : ' READ (LKB,696) RMT C* - Desired material is MAT IF(RMT.GT.0) MAT=RMT+0.01 IF(RMT.LT.0) GO TO 681 MAT1=MAT C* Read the source library header and list of nuclides 141 READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD READ (LF1,701) (IWO(I),I=1,NEL) KEL =NEL IF (KG0.LE.0) KG0=NG0 KNFD=NNFD KFPD=NFPD C* - Case: no new material data to be processed IF(MAT.LE.0) GO TO 146 DO 144 I=1,NEL IF(IWO(I).NE.MAT) GO TO 144 C* - Case: material already exists - Define MAT1 to insert before MAT WRITE(LTT,693) ' Material ident. already exists MAT: ',MAT WRITE(LTT,691) ' Enter blank to replace MAT ' WRITE(LTT,691) ' negative to quit ' WRITE(LTT,691) '$ new ident to insert before MAT: ' READ (LKB,696) RMT WRITE(LTT,691) IF(RMT.GT.0) MAT1=RMT+0.01 IF(RMT.LT.0) GO TO 681 C* -- Material MAT is to be replaced MAT=MAT1 IF(MAT.EQ.MAT1) GO TO 146 C* -- Material MAT1 is to be inserted before material MAT MT2=IWO(I) IWO(I)=MAT1 DO 142 J=I,NEL MT1=IWO(J+1) IWO(J+1)=MT2 MT2=MT1 142 CONTINUE GO TO 145 144 CONTINUE C* - Case: material is to be inserted IWO(NEL+1)=MAT C* Add burnup data if missing 145 IF(LBU.LE.0) THEN NDAT1=2 LBU=LSC LSC=LSC+NDAT1/2 RWO(LBU+1)=0. END IF IWO(LBU+1)=MAT1 KEL =NEL+1 KNFD=NNFD KFPD=NFPD IF(NDAT1.GT.8) THEN KNFD=KNFD+1 ELSE IF(NDAT1.EQ.8) THEN IF (IWO(LBU+4).LT.2) THEN KFPD=KFPD+1 ELSE KNFD=KNFD+1 ENDIF ENDIF C* Write the library header 146 WRITE(LF3,701) KEL,NG,KG0,NG1,NG2,NG3,KNFD,KFPD WRITE(LF3,701) (IWO(I),I=1,KEL) C* Write the group boundaries KLEN=NG+1 READ (LF1,702) (RWO( I),I=1,KLEN) WRITE(LF3,702) (RWO(LGB+I),I=1,KLEN) C* Write the fission spectrum READ (LF1,702) (RWO( I),I=1,NG0) WRITE(LF3,702) (RWO(LSP+I),I=1,KG0) C* Process the burnup chains data DO 324 J=1,NEL READ (LF1,701) JDAT1 READ (LF1,703) (RWO( I),IWO( I),I=1,JDAT1/2) L =0 IF(IWO(1).NE.MAT .OR. LBU.EQ.0) GO TO 322 IF(MAT.EQ.MAT1) GO TO 321 C* - Insert the burnup data before MAT WRITE(LF3,701) NDAT1 WRITE(LF3,703) (RWO(LBU+I),IWO(LBU+I),I=1,NDAT1/2) WRITE(LTT,693) ' Burnup data added for material ',MAT1 GO TO 322 C* - Replace the burnup data if same material is processed 321 L =LBU IF (JDAT1.EQ.2 .AND. NDAT1.GT.2) THEN WRITE(LTT,691) ' WARNING - A non burnable material repla' 1 ,'ced by a burnable one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ELSE IF (JDAT1.EQ.8) THEN IF (NDAT1.LT.8) THEN WRITE(LTT,691) ' WARNING - A burnable material replaced ' 1 ,'by a non-burnable one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ELSE IF ((IWO(4).LT.2) .AND. (IWO(L+4).GT.2)) THEN WRITE(LTT,691) ' WARNING - A burnable material replaced ' 1 ,'by a fissile one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ELSE IF ((IWO(4).GT.2) .AND. (IWO(L+4).LT.2)) THEN WRITE(LTT,691) ' WARNING - A fissile material replaced b' 1 ,'y a burnable one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ENDIF ELSE IF (NDAT1.LT.8) THEN WRITE(LTT,691) ' WARNING - A fissile material replaced b' 1 ,'y a non-burnable one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ELSE IF ((NDAT1.EQ.8) .AND. (IWO(L+4).LT.2)) THEN WRITE(LTT,691) ' WARNING - A fissile material replaced b' 1 ,'y a burnable one. ' WRITE(LTT,691) ' Run WILLIE CHCK Option... ' ENDIF ENDIF JDAT1=NDAT1 WRITE(LTT,693) ' Burnup data replaced for material ',MAT 322 WRITE(LF3,701) JDAT1 WRITE(LF3,703) (RWO(L+I),IWO(L+I),I=1,JDAT1/2) 324 CONTINUE IF(KEL.EQ.NEL .OR. MAT.NE.MAT1) GO TO 326 C* - Add the burnup data for the new material IF(LBU.GT.0) GO TO 325 NDAT1=2 RWO(1)=0. IWO(1)=MAT 325 WRITE(LF3,701) NDAT1 WRITE(LF3,703) (RWO(LBU+I),IWO(LBU+I),I=1,NDAT1/2) WRITE(LTT,693) ' Burnup data added for material ',MAT C* Check for the file mark at the end of burnup data 326 IER=51 READ (LF1,701) JTERM IF(JTERM.NE.ITERM) GO TO 690 WRITE(LF3,701) ITERM C* Copy the nuclide cross section data from the old library DO 434 J=1,NEL LFO=LF3 READ (LF1,704) ID,AW,IZ,NF,NT,NR IF(ID.NE.MAT) GO TO 328 C* Copy the new cross section data WRITE(LF3,704) MAT1,WW,KZ,NFBU,KT,KR CALL PROCXS(RWO,IWO,LF2,LF3,NG1,NG2,NG3,KF,KT) C* Mark that resonance data need to be inserted IF(KR.GT.0) LRI=1 IF(MAT.EQ.MAT1) THEN C* If replaced material has resonance data, overwrite IF(NR.GT.0) LRI=0 WRITE(LTT,693) ' Cross sections replaced for material ',MAT ELSE WRITE(LTT,693) ' Cross sections inserted for material ',MAT1 END IF IER=51 READ (LF2,701) JTERM, N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 690 WRITE(LF3,701) ITERM C* To replace cross section data, skip file in the orig. library IF(MAT.EQ.MAT1) LFO=0 C* Copy or skip the cross section data from the original file 328 IF(LFO.GT.0) THEN WRITE(LF3,704) ID,AW,IZ,NF,NT,NR C* Mark the MAT number before which to insert resonance data IF(LRI.GT.0 .AND. NR.GT.0) THEN MAT=ID LRI=0 END IF END IF CALL PROCXS(RWO,IWO,LF1,LFO,NG1,NG2,NG3,NF,NT) IER=51 READ (LF1,701) JTERM, N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 690 IF(LFO.GT.0) WRITE(LF3,701) ITERM 434 CONTINUE IF(KEL.EQ.NEL .OR. MAT.NE.MAT1) GO TO 500 C* Add the nuclide cross section data for the new material at the end WRITE(LF3,704) MAT,WW,KZ,NFBU,KT,KR CALL PROCXS(RWO,IWO,LF2,LF3,NG1,NG2,NG3,KF,KT) WRITE(LTT,693) ' Cross sections added for material ',MAT IER=52 READ (LF2,701) JTERM, N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 690 WRITE(LF3,701) ITERM C* Loop over all groups for the resonance integrals 500 DO 550 J=1,NG2 C* Determine the number of new resonance tables NRT=KR*( (KF+1)/2 ) IF(KF.EQ.4) NRT=0 C* Read the resonance integrals from the old library 532 L1=0 READ (LF1,706) M12,RID,M1,M2 READ (LF1,702) (RWO(L1+I ),I=1,M1), * (RWO(L1+I+100),I=1,M2), * (RWO(L1+I+200),I=1,M12) L2=L1+200+M12 ID=RID+0.001 IF(ID.GT.0 .AND. ID.NE.MAT) GO TO 536 C* Write the new resonance integrals IF(NRT.LE.0) GO TO 535 DO 534 K=1,NRT READ (LF2,706) J12,RJD,J1,J2 READ (LF2,702) (RWO(L2+I ),I=1,J1), * (RWO(L2+I+100),I=1,J2), * (RWO(L2+I+200),I=1,J12) IF(RMT.GT.0) RJD=RMT WRITE(LF3,706) J12,RJD,J1,J2 WRITE(LF3,702) (RWO(L2+I ),I=1,J1), * (RWO(L2+I+100),I=1,J2), * (RWO(L2+I+200),I=1,J12) 534 CONTINUE NRT=0 535 IF(ID.EQ.MAT .AND. MAT.EQ.MAT1) GO TO 532 C* Copy the old resonance integrals 536 WRITE(LF3,706) M12,RID,M1,M2 WRITE(LF3,702) (RWO(L1+I ),I=1,M1), * (RWO(L1+I+100),I=1,M2), * (RWO(L1+I+200),I=1,M12) C* Group resonance integral table terminated by dummy RID=0 IF(RID.NE. 0. ) GO TO 532 C* Check the file marks IER=53 READ (LF1,701) JTERM, K1,K2,K3,K4 IF(JTERM.NE.ITERM) GO TO 690 IER=54 IF(KR.GT.0) READ (LF2,701) JTERM, N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 690 WRITE(LF3,701) ITERM 550 CONTINUE IF(KR.GT.0) 1WRITE(LTT,694) ' Resonance integrals written for mat. ',RJD C* Check for any P1 scattering data JP1=-1 LBL=1 L2 =LBL IC2=0 610 IF(N1.EQ.0 .OR. N1.EQ.6) GO TO 612 CALL SKIPWF(LF2,ITERM,JTERM,N1,N2,N3,N4) IF(JTERM.LT.0 .OR. N2.LT.0) GO TO 630 GO TO 610 612 READ (LF2,692,END=630) REN C* P1 matrix to be added - inquire about the position WRITE(LTT,691) ' Enter the P1-scattering matrix position' WRITE(LTT,691) '$ "n" to replace, blank to append : ' READ (LKB,695) JP1 C* Identify the format of the data to be added READ (REN,701,ERR=630) KL2 READ (LF2,702) (RWO(L2-1+J),J=1,KL2) IC2=1 LBL=LBL+KL2 C* Copy the P1 scattering matrices from the old library 630 IC1=0 NP1=0 640 NP1=NP1+1 C* Identify the format on the original library (full or condensed) C* Next card may be the last filemark, matrix length or data IER=56 READ (LF1,692,END=690) REC READ (REC,701,ERR=650) KL1 IF(KL1.NE.ITERM) GO TO 642 C* Filemark found - Check if new data had already been written IF(JP1.EQ.0 .OR. JP1.GE.NP1) GO TO 660 GO TO 680 C* Case-source lib.: copy condensed P1 matrices 642 L1 =LBL READ (LF1,702) (RWO(L1-1+J),J=1,KL1) IF(JP1.EQ.NP1) GO TO 660 WRITE(LTT,693) ' Copy condensed P1-matrix ',NP1 WRITE(LF3,701) KL1 WRITE(LF3,702) (RWO(L1-1+J),J=1,KL1) GO TO 640 C* Case-source lib.: copy uncondensed P1 matrices (as char.strings) 650 NR =(NG+4)/5 JG =0 IF(JP1.NE.NP1) 1WRITE(LTT,693) ' Copy expanded P1-matrix ',NP1 651 JG =JG+1 IF(JP1.NE.NP1) WRITE(LF3,692) REC DO 652 K=2,NR READ (LF1,692) REC IF(JP1.NE.NP1) WRITE(LF3,692) REC 652 CONTINUE IF(JG.GE.NG) GO TO 654 READ (LF1,692) REC GO TO 651 654 IF(JP1.NE.NP1) GO TO 640 C* Check how the new data are to be written 660 IF(IC2.EQ. 0) GO TO 670 IF(IC1.EQ. 0) GO TO 664 C* Case-target lib.: New data may be written in condensed format WRITE(LTT,693) ' Write condensed P1-matrix ',NP1 WRITE(LF3,701) KL2 WRITE(LF3,702) (RWO(L2-1+J),J=1,KL2) GO TO 640 C* Case-target lib.: Expand new data to uncondensed format 664 JG=0 LJ=L2 LS=LBL WRITE(LTT,693) ' Expand new P1-matrix ',NP1 665 JG=JG+1 DO 667 J=1,NG RWO(LS-1+J)=0. 667 CONTINUE NU =RWO(LJ )+0.1 NS =RWO(LJ+1)+0.1 DO 668 J=1,NS RWO(LS+JG-NU-1+J)=RWO(LJ+1+J) 668 CONTINUE WRITE(LF3,702) (RWO(LS-1+J),J=1,NG) LJ=LJ+2+NS IF(JG.LT.NG) GO TO 665 GO TO 640 C* Case-target lib.: New data in uncondensed format - copy as characters 670 JG =0 WRITE(LTT,693) ' Write expanded P1-matrix ',NP1 671 JG =JG+1 WRITE(LF3,692) REN IER=58 NR =(NG+4)/5 DO 672 K=2,NR READ (LF2,692,ERR=690) REN WRITE(LF3,692) REN 672 CONTINUE IF(JG.GE.NG) GO TO 640 C* Check if any more matrices are to be processed READ (LF2,692) REN GO TO 671 C* All P1 scattering matrices processed - finish 680 WRITE(LF3,701) ITERM C* All processing completed 681 CLOSE(UNIT=LF3) 682 CLOSE(UNIT=LF2) 683 CLOSE(UNIT=LF1) 684 RETURN C* Error condition reading the library 690 CALL ERRWAR(LTT,IER) RETURN 691 FORMAT(2A40) 692 FORMAT(A75) 693 FORMAT(A40,I8) 694 FORMAT(A40,F8.1) 695 FORMAT(BN,I10) 696 FORMAT(BN,F10.0) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) 706 FORMAT(1P,I15/E15.8,2I6) END *DECK MIXMAT SUBROUTINE MIXMAT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : MIXMAT Subroutine C-Purpose: Mix WIMS library cross section data sets C-Description: C-D The module allows wimsing of cross sections of fission spectra C-D to make a new composite data set C-D Cross section files such as produced by WILLIE option COPY C-D are processed. The module generously allows labelled, unlabelled C-D as well as no-filemark formats. If fissin spectra are processed, C-D the source file must be labelled (NTYPE=2). Cross sections from C-D several files are read, each scaled by the weighting factor C-D defined from input. The result is a new material which is written C-D to output in the WIMS-D format so that it can be processed further C-D by other modules if desired. One of the purposes of this module C-D is to generate element cross section files from isotopic data. C-D The WIMS-D/4 code allows resonance integrals to be specified C-D for materials appearing in the fuel only. Such materials are C-D usually given by isotopes. Mixing of the resonance integral data C-D is not supported by this module. C-D CHARACTER*40 FLNI,FLNO,FLNL, BLNK DIMENSION RWO(*),IWO(*),YD(4),ID(4) COMMON /TERWOR/ ITERM,ALDO(2) C* Default group structure (fast/resonance/thermal) DATA NG1,NG2,NG3/ 14, 13, 42 / C* Alternative group structure DATA MG1,MG2,MG3/ 45, 47, 80 / C* Default output filename DATA BLNK /' '/ 7 FLNL /'MIXMAT.LOG'/ C* Local file units DATA LIN,LOU,LFL / 1, 2, 7 / C* Limit for the number of thermal temperatures MXT=6 C* Initialize parameters AB0=0. AW0=0. ZZ0=0. NF0=0 NT0=0 NR0=0 NMT=0 NOD=0 C* Reserve and initialize the output cross section field C* L1 - pot., sl.dwn.pwr, epith.transp., absorpt., Xi and Lambda C* L5 - epithermal neutron fission yield and fission cross section C* L7 - epithermal scattering matrix C* L8 - thermal temperature list C* L9 - thermal transport and absorption cross section C* L10- thermal neutron fission yield and fission cross section C* L12- thermal scattering matrix L1 = 1 LBL= 1 C... N12= NG1+NG2 C... NGR= N12+NG3 C... NSP= N12 C... L1 = 1 C... L5 = L1 +NG2*4 +N12*2 C... L7 = L5 +N12*2 C... L8 = L7 +NGR*N12 C... L9 = L8 +MXT C... L10= L9 +MXT*NG3*2 C... L12= L10+MXT*NG3*2 C... LBL= L12+MXT*NG3*NG3 C... IF(LBL.GT.MXRW) STOP 'WILLIE MIXMAT - Array capacity exceeded 1' C... DO 10 I=L1,LBL C... RWO(I)=0. C...10 CONTINUE C* Open the log file for diagnostic messages C***** VAX ***** C OPEN(UNIT=LFL,FILE=FLNL,STATUS='NEW',CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN(UNIT=LFL,FILE=FLNL,STATUS='UNKNOWN') C***** EXPORT ***** C* Write the header to logfile and terminal WRITE(LFL,606) WRITE(LFL,607) GO TO 6 5 WRITE(LTT,691) ' ERROR opening file : ',FLNO 6 WRITE(LTT,691) ' Enter the OUTPUT cross section filename' WRITE(LTT,691) '$ BLANK = Default = MIXMAT.XSW : ' READ (LKB,691) FLNO IF(FLNO.EQ.BLNK)FLNO='MIXMAT.XSW' C* Open the output cross section filename C***** VAX ***** C OPEN(UNIT=LOU,FILE=FLNO,STATUS='NEW',CARRIAGECONTROL='LIST',ERR=5) C***** VAX ***** C***** EXPORT ***** OPEN(UNIT=LOU,FILE=FLNO,STATUS='UNKNOWN',ERR=5) C***** EXPORT ***** WRITE(LTT,691) ' Enter WIMS ID. for the new material ' WRITE(LTT,691) '$ ID.LE.0 = Default = 1000*Z+A : ' READ (LKB,695) MWID GO TO 12 C* Define the source cross section filename 11 WRITE(LTT,691) ' ERROR reading file : ',FLNI 12 WRITE(LTT,691) WRITE(LTT,691) ' Enter INPUT cross section data filename' WRITE(LTT,691) '$ BLANK = To Finish : ' READ (LKB,691) FLNI IF(FLNI.EQ.BLNK) GO TO 500 C***** VAX ***** C OPEN(UNIT=LIN,FILE=FLNI,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN(UNIT=LIN,FILE=FLNI,STATUS='OLD',ERR=11) C***** EXPORT ***** C* Define the abundance WRITE(LTT,691) ' Enter abundance ( 0 to skip) [at.%] : ' READ (LKB,694) AB1 IF(AB1.GT. 0 ) GO TO 18 CLOSE(UNIT=LIN) GO TO 12 C* Try to interpret the filemark record 18 READ (LIN,701,ERR=22) JTERM, NTYP,JG1,JG2,JG3 20 IF(NTYP.EQ.0 .OR. NTYP.EQ.4) GO TO 32 IF(NTYP.EQ.2) GO TO 40 IF(NTYP.EQ.3) GO TO 25 IF(NTYP.LT.4) GO TO 30 GO TO 38 C* No labelled filemark present - define group structure 22 REWIND LIN IF(NMT.GT.0) GO TO 400 C... WRITE(LTT,691) C... WRITE(LTT,691) ' WARNING - No labelled filemarks on file' C... WRITE(LTT,691) '$ Enter number of fast groups : ' C... READ (LKB,695) JG1 C... WRITE(LTT,691) '$ Enter number of resonance groups : ' C... READ (LKB,695) JG2 C... WRITE(LTT,691) '$ Enter number of thermal groups : ' C... READ (LKB,695) JG3 GO TO 32 C* Burnup data record 25 READ (LIN,701) NNN NNNN=NNN/2 IF (NNNN.GT.4) THEN CLOSE(LIN) CLOSE(LOU) STOP ' WILLIE MIXMAT: Burnup Record Array capacity excedeed 2 ' ENDIF IF (NOD.EQ.0) THEN NOD=NNN ELSE IF(NNN.NE.NOD) THEN WRITE(*,*)' Different burnup record length found, last assumed' ENDIF ENDIF READ (LIN,707) (YD(I),ID(I),I=1,NNNN) C* Inappropriate filemark label found - skip to the next filemark 30 READ (LIN,701,ERR=30) JTERM, NTYP,JG1,JG2,JG3 GO TO 20 C* Labelled filemark for cross section data found 32 IF(NMT.GT.0) GO TO 36 C* First labelled filemarks file may redefine the groups IF(JG1.EQ. 0 ) JG1=NG1 IF(JG2.EQ. 0 ) JG2=NG2 IF(JG3.EQ. 0 ) JG3=NG3 C* Redefine the group structure and re-initialize the work field NG1=JG1 NG2=JG2 NG3=JG3 N12= NG1+NG2 NGR= N12+NG3 L1 = 1 L5 = L1 +NG2*4 +N12*2 L7 = L5 +N12*2 L8 = L7 +NGR*N12 L9 = L8 +MXT L10= L9 +MXT*NG3*2 L12= L10+MXT+NG3*2 LBL= L12+MXT*NG3*NG3 IF(LBL.GT.MXRW) STOP 'MIXMAT ERROR - MXRW limit exceeded' DO 34 I=L1,LBL RWO(I)=0. 34 CONTINUE GO TO 400 C* Subsequent files must have the same group structure 36 IF(JG1.GT.0 .AND. JG1.NE.NG1) GO TO 37 IF(JG2.GT.0 .AND. JG2.NE.NG2) GO TO 37 IF(JG3.GT.0 .AND. JG3.NE.NG3) GO TO 37 GO TO 400 C* Error condition on group structure 37 WRITE(LTT,691) ' Incompatible group structure on file : ',FLNI GO TO 12 C* No cross section data found 38 WRITE(LTT,691) ' No cross section data on file : ',FLNI GO TO 12 C* Labelled filemark for fission spectrum data found 40 WRITE(LTT,691) ' Fission spectrum data encountered ' IF(NMT.GT.0) GO TO 100 C* First labelled filemarks file may redefine the spectral groups IF(JG1.EQ. 0 ) JG1=NSP NSP=JG1 IF(LBL.LE.L1) LBL=L1+NSP DO 42 I=1,NSP RWO(L1-1+I)=0. 42 CONTINUE C* C* Begin processing the fission spectra 100 WRITE(LFL,608) MT0,AW0,IZ0,NF0,NT0,NR0,AB1,FLNI NMT= NMT+ 1 AB0= AB0+ AB1 ABF= AB1* 0.01 READ(LIN,702) (RWO(LBL-1+I),I=1,NSP) DO 120 I=1,NSP 120 RWO(L1-1+I)=RWO(L1-1+I) + ABF*RWO(LBL-1+I) C* Fission spectrum processing completed - any additional data ignored CLOSE (UNIT=LIN) GO TO 12 C* C* Begin processing the cross sections 400 READ (LIN,704) MT1,AW1,IZ1,NF1,NT1,NR1 C* Request additional information for the current material WRITE(LTT,607) WRITE(LTT,608) MT1,AW1,IZ1,NF1,NT1,NR1,AB1,FLNI C* Check the parameter for the thermal temperature list 401 WRITE(LFL,608) MT1,AW1,IZ1,NF1,NT1,NR1,AB1,FLNI IF(NT0.EQ.0) NT0=NT1 IF(NT1.EQ.1 .AND. NT0.GT.1) WRITE(LTT,691) 1 ' MIXMAT WARNING - No temperature depende' 2,'nce for this material ' IF(NT1.GT.1 .AND. NT0.NE.NT1) THEN WRITE(LTT,691) ' MIXMAT ERROR - Incompatible temperature' 1 ,' list - data abandoned ' CLOSE(UNIT=LIN) GO TO 12 END IF C* Check the fission flag parameter IF(NF1.GT.0 .AND. NF1.NE.4) WRITE(LTT,691) 1 ' MIXMAT WARNING - Res. tables ignored ' IF(NF1.GE.2) NF0=4 C* Check the resonance tables parameter IF(NR1.GT.0) WRITE(LTT,691) 1 ' MIXMAT WARNING - Res. tables ignored ' C* Calculate averages NMT= NMT+ 1 AB0= AB0+ AB1 ABF= AB1* 0.01 AW0= AW0+ ABF*AW1 ZZ0= ZZ0+ ABF*FLOAT(IZ1) C* C* Epithermal cross sections processing KLN=4*NG2+2*N12 KLF=2*N12 LSC=MAX(LBL+N12*NGR, LBL+KLN, LBL+KLF) IF(LSC.GT.MXRW) STOP 'WILLIE MIXMAT - Array capacity exceeded 2' C* Process the epithermal cross sections READ(LIN,702,ERR=490) (RWO(LBL-1+I),I=1,KLN) DO 414 I=1,KLN RWO(L1-1+I)=RWO(L1-1+I) + ABF*RWO(LBL-1+I) 414 CONTINUE C* Process the epithermal fission yield and x-set. (if present) IF(NF1.LT.2) GO TO 420 READ(LIN,702) (RWO(LBL-1+I),I=1,KLF) DO 416 I=1,KLF RWO(L5-1+I)=RWO(L5-1+I) + ABF*RWO(LBL-1+I) 416 CONTINUE C* Process the epithermal scattering matrix 420 READ(LIN,705,ERR=490) KL1,(RWO(LBL-1+I),I=1,KL1) LB1=LBL DO 422 I=1,N12 LBB=LB1 LBS=RWO(LBB) +0.1 LBT=RWO(LBB+1)+0.1 LB1=LBB+LBT+2 J12=L7+(I-1)*NGR + I-LBS DO 422 J=1,LBT RWO(J12-1+J)=RWO(J12-1+J) + ABF*RWO(LBB+1+J) 422 CONTINUE C* C* Thermal cross sections processing KLN=2*NG3 KLF=2*NG3 K9 =LBL K10=K9 +KLN K12=K10+KLF LSC=K12+NG3*NG3 IF(LSC.GT.MXRW) STOP 'WILLIE MIXMAT - Array capacity exceeded 3' C* Read and check the thermal cross sections temperature list READ(LIN,702) (RWO(LBL-1+I),I=1,NT1) DO 432 I=1,NT1 IF(NMT.EQ.1) RWO(L8-1+I)=RWO(LBL-1+I) IF(ABS(RWO(L8-1+I)-RWO(LBL-1+I)).GT. 0.5) WRITE(LFL,691) 1 ' MIXMAT WARNING - Incompatible temperatu' 2,'re list encountered ' 432 CONTINUE C* Loop over thermal temperatures DO 480 IT0=1,NT0 M9 = L9+(IT0-1)*KLN M10=L10+(IT0-1)*KLF M12=L12+(IT0-1)*NG3*NG3 C* Read all thermal data for one temperature IF(IT0.GT.NT1) GO TO 450 READ(LIN,702) (RWO(K9-1+I),I=1,KLN) IF(NF1.GE.2) READ(LIN,702) (RWO(K10-1+I),I=1,KLF) READ(LIN,705) KL1,(RWO(K12-1+I),I=1,KL1) C* Process the thermal cross sections 450 DO 454 I=1,KLN RWO( M9-1+I)=RWO( M9-1+I) + ABF*RWO( K9-1+I) 454 CONTINUE C* Process the thermal fission yield and x-set. (if present) IF(NF1.LT.2) GO TO 460 DO 456 I=1,KLF RWO(M10-1+I)=RWO(M10-1+I) + ABF*RWO(K10-1+I) 456 CONTINUE C* Process the thermal scattering matrix 460 LB1=K12 DO 462 I=1,NG3 LBB=LB1 LBS=RWO(LBB) +0.1 LBT=RWO(LBB+1)+0.1 LB1=LBB+LBT+2 J12=M12+(I-1)*NG3 DO 462 J=1,LBT C* Suppress upscattering into the epithermal groups JJ =MAX(0, J-1+I-LBS) RWO(J12+JJ)=RWO(J12+JJ) + ABF*RWO(LBB+1+J) 462 CONTINUE 480 CONTINUE C* Cross section data processing completed - any additional data ignored CLOSE (UNIT=LIN) GO TO 12 C* On read error try for 172-group library 490 IF(NGR.EQ.172 .OR. NMT.GT.1) 1 STOP 'MIXMAT ERROR - Unknown group structure' NMT=NMT-1 NG1=MG1 NG2=MG2 NG3=MG3 JG1=0 JG2=0 JG3=0 NGR=NG1+NG2+NG3 REWIND LIN WRITE(LTT,691) WRITE(LTT,691) ' MIXMAT WARNING - Undefined group struct' 1 ,'ure - Try 172-group scheme ' GO TO 18 C* C* Write the new WIMS-D library file 500 XAB0=ABS(1.0-AB0/100.0) IF (XAB0.GT.1.0E-5) THEN AW0=100.0*AW0/AB0 ZZ0=100.0*ZZ0/AB0 IL0=3*NG2+2*N12 DO IL=1,NG2 RWO( L1-1+IL0+IL)=100.0*RWO(L1-1+IL0+IL)/AB0 END DO ENDIF IZ0=ZZ0+0.5 IA0=AW0+0.5 MT0=IZ0*1000+IA0 IF(MWID.GT.0) MT0=MWID WRITE(LFL,691) '========================================' WRITE(LFL,608) MT0,AW0,IZ0,NF0,NT0,NR0,AB0,FLNO IF(NTYP.NE.2) GO TO 510 C* C* Write the fission spectrum WRITE(LOU,701) ITERM, NTYP,NSP WRITE(LOU,702) (RWO(L1 -1+I),I=1,NSP) GO TO 600 C* C* Write burnup information, if any 510 IF(NOD.GT.0) THEN NTYP=3 WRITE(LOU,701) ITERM, NTYP NNNN=NNN/2 ID(1)=MT0 WRITE(LOU,701) NNN WRITE(LOU,707)(YD(I),ID(I),I=1,NNNN) NTYP=4 WRITE(LOU,701) ITERM, NTYP ENDIF C* C* Write the nuclide header information WRITE(LOU,704) MT0,AW0,IZ0,NF0,NT0,NR0 C* Write the epithermal cross sections KLN=4*NG2+2*N12 KLF=2*N12 WRITE(LOU,702) (RWO( L1-1+I),I=1,KLN) C* Write the epithermal fission yield and x-set. (if present) IF(NF1.LT.2) GO TO 520 WRITE(LOU,702) (RWO( L5-1+I),I=1,KLF) C* Process the epithermal scattering matrix 520 LBB=LBL LSC=LBL+N12*NGR IF(LSC.GT.MXRW) STOP 'WILLIE MIXMAT - Array capacity exceeded 4' DO 540 I=1,N12 M7 =L7+(I-1)*NGR C* Find the self-scattering position in the reduced scatt.matrix LBS=I 522 LB1=I-LBS IF(RWO(M7+LB1).GT.0 .OR. LBS.LE.1) GO TO 524 LBS=LBS-1 GO TO 522 C* Find the current row length of the reduced scattering matrix 524 LBT=NGR-LB1 525 IF(RWO(M7-1+LBT+LB1).GT.0 .OR. LBT.LE.LBS) GO TO 526 LBT=LBT-1 GO TO 525 C* Move the current row to a scratch array 526 RWO(LBB )=LBS RWO(LBB+1)=LBT DO 528 J=1,LBT 528 RWO(LBB+1+J)=RWO(M7+LB1-1+J) LBB=LBB+LBT+2 540 CONTINUE KL1=LBB-LBL WRITE(LOU,705) KL1,(RWO(LBL-1+I),I=1,KL1) C* C* Write the temperature list WRITE(LOU,702) (RWO( L8-1+I),I=1,NT0) KLN=2*NG3 KLF=2*NG3 C* Loop over all thermal temperatures DO 580 IT0=1,NT0 M9 = L9+(IT0-1)*KLN M10=L10+(IT0-1)*KLF M12=L12+(IT0-1)*NG3*NG3 C* Write the thermal cross sections WRITE(LOU,702) (RWO( M9-1+I),I=1,KLN) C* Write the epithermal fission yield and x-set. (if present) IF(NF1.LT.2) GO TO 560 WRITE(LOU,702) (RWO(M10-1+I),I=1,KLF) C* Process the epithermal scattering matrix 560 LBB=LBL LSC=LBL+NG3*NG3 IF(LSC.GT.MXRW) STOP 'WILLIE MIXMAT - Array capacity exceeded 4' DO 570 I=1,NG3 J12=M12+(I-1)*NG3 C* Find the self-scattering position in the reduced scatt.matrix LBS=I 562 LB1=I-LBS IF(RWO(J12+LB1).GT.0 .OR. LBS.LE.1) GO TO 564 LBS=LBS-1 GO TO 562 C* Find the current row length of the reduced scattering matrix 564 LBT=NG3-LB1 565 IF(RWO(J12-1+LBT+LB1).GT.0 .OR. LBT.LE.LBS) GO TO 566 LBT=LBT-1 GO TO 565 C* Move the current row to a scratch array 566 RWO(LBB )=LBS RWO(LBB+1)=LBT DO 568 J=1,LBT 568 RWO(LBB+1+J)=RWO(J12+LB1-1+J) LBB=LBB+LBT+2 570 CONTINUE KL1=LBB-LBL WRITE(LOU,705) KL1,(RWO(LBL-1+I),I=1,KL1) 580 CONTINUE C* New file written 600 WRITE(LOU,701)ITERM CLOSE (UNIT=LOU) CLOSE (UNIT=LFL) RETURN 606 FORMAT(// ' WILLIE - X-SECT. MIXING UTILITY '/ 1 ' ================================'/) 607 FORMAT(/' Mat.ID At.Wt. Z Nf Nt Nr Abund. Filename') 608 FORMAT( I8, F9.4, I4,I3,I3,I3, F9.4,1X,A40) 691 FORMAT(2A40) 694 FORMAT(BN,F20.0) 695 FORMAT(BN,I10) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) C 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 705 FORMAT(1P,I15/(5E15.8)) C 706 FORMAT(1P,E15.8,2I6) 707 FORMAT(1P,(3(E15.8,I6))) END *DECK DELEMT SUBROUTINE DELEMT(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : Subroutine DELEMT C-Purpose: Delete a material from the WIMS-D library CHARACTER*75 REC CHARACTER*40 BLNK, FLNM, FLWI,FLWO DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) DATA BLNK /' '/ 1 ,FLWI /'WIMS.LIB'/ 1 ,FLWO /'WIMS.LIB'/ C* Default file units DATA LF1,LF3/ 1, 3 / C* Define the filenames GO TO 12 11 WRITE(LTT,691) ' ERROR opening file : ',FLWI 12 WRITE(LTT,691) ' Default source formatted lib.filename: ',FLWI WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691,END=684) FLNM IF(FLNM.NE.BLNK) FLWI=FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11) C***** EXPORT ***** GO TO 16 15 WRITE(LTT,691) ' ERROR opening file : ',FLWO 16 WRITE(LTT,691) ' Default new formatted lib. filename: ',FLWO WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691,END=683) FLNM IF(FLNM.NE.BLNK) FLWO=FLNM C***** EXPORT ***** IF(FLWI.EQ.FLWO) GO TO 15 C***** EXPORT ***** C***** VAX ***** C OPEN (UNIT=LF3,FILE=FLWO,STATUS='NEW',ERR=15 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF3,FILE=FLWO,STATUS='UNKNOWN',ERR=15) C***** EXPORT ***** C* Read the MAT number of the material to be deleted WRITE(LTT,691) '$Enter MAT No.of material to remove : ' READ (LKB,696) RMT IF(RMT.GT.0) MAT=RMT+0.01 C* C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' Source WIMS-D library filename : ',FLWI WRITE(LTT,691) ' New WIMS-D library filename : ',FLWO WRITE(LTT,694) ' Removing material MAT No.: ',RMT WRITE(LTT,691) C* C* Read the source library header and list of nuclides 141 READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD READ (LF1,701) (IWO(I),I=1,NEL) KNFD=NNFD KFPD=NFPD C* Check the MAT number validity DO 142 I=1,NEL IF(IWO(I).EQ.MAT) GO TO 144 142 CONTINUE WRITE(LTT,693) ' DELM ERROR - Can not find material MAT ',MAT STOP C* Check the burnup data 144 KLEN=NG+1 READ (LF1,702) (RWO( I),I=1,KLEN) READ (LF1,702) (RWO( I),I=1,NG0) DO 154 J=1,NEL READ (LF1,701) JDAT1 READ (LF1,703) (RWO( I),IWO( I),I=1,JDAT1/2) IF(IWO(1).EQ.MAT) GO TO 162 154 CONTINUE WRITE(LTT,693) ' DELM ERROR - Error in burnup data MAT ',MAT STOP 162 IF(JDAT1.GT.8) THEN KNFD =KNFD-1 WRITE(LTT,693) ' Fissile material deleted ' ELSE IF(JDAT1.EQ.8) THEN IF(IWO(4).GE.2) THEN KNFD =KNFD-1 WRITE(LTT,693) ' Fissile material deleted ' ELSE KFPD=KFPD-1 WRITE(LTT,693) ' Burnable material deleted ' END IF END IF C* C* Process the library, deleting material MAT REWIND LF1 READ (LF1,701) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NFPD READ (LF1,701) (IWO(I),I=1,NEL) KEL=0 DO 164 I=1,NEL IWO(KEL+1)=IWO(I) IF(IWO(I).NE.MAT) KEL=KEL+1 164 CONTINUE C* Write the library header WRITE(LF3,701) KEL,NG,NG0,NG1,NG2,NG3,KNFD,KFPD WRITE(LF3,701) (IWO(I),I=1,KEL) C* Write the group boundaries KLEN=NG+1 READ (LF1,702) (RWO(I),I=1,KLEN) WRITE(LF3,702) (RWO(I),I=1,KLEN) C* Write the fission spectrum READ (LF1,702) (RWO(I),I=1,NG0) WRITE(LF3,702) (RWO(I),I=1,NG0) C* Process the burnup chains data DO 324 J=1,NEL READ (LF1,701) JDAT1 READ (LF1,703) (RWO(I),IWO(I),I=1,JDAT1/2) IF(IWO(1).EQ.MAT) GO TO 324 WRITE(LF3,701) JDAT1 WRITE(LF3,703) (RWO(I),IWO(I),I=1,JDAT1/2) 324 CONTINUE C* Check for the file mark at the end of burnup data 326 IER=51 READ (LF1,701) JTERM IF(JTERM.NE.ITERM) GO TO 690 WRITE(LF3,701) ITERM C* Copy the nuclide cross section data from the old library DO 434 J=1,NEL LFO=LF3 READ (LF1,704) ID,AW,IZ,NF,NT,NR IF(ID.EQ.MAT) LFO=0 C* Copy or skip the cross section data from the original file IF(LFO.GT.0) WRITE(LF3,704) ID,AW,IZ,NF,NT,NR CALL PROCXS(RWO,IWO,LF1,LFO,NG1,NG2,NG3,NF,NT) IER=51 READ (LF1,701) JTERM, N1,N2,N3,N4 IF(JTERM.NE.ITERM) GO TO 690 IF(LFO.GT.0) WRITE(LF3,701) ITERM 434 CONTINUE C* Loop over all groups for the resonance integrals 500 DO 550 J=1,NG2 C* Read the resonance integrals from the old library 532 L1=0 READ (LF1,706) M12,RID,M1,M2 READ (LF1,702) (RWO(L1+I ),I=1,M1), * (RWO(L1+I+100),I=1,M2), * (RWO(L1+I+200),I=1,M12) ID=RID+0.001 IF(ID.EQ.MAT) GO TO 537 C* Copy the old resonance integrals WRITE(LF3,706) M12,RID,M1,M2 WRITE(LF3,702) (RWO(L1+I ),I=1,M1), * (RWO(L1+I+100),I=1,M2), * (RWO(L1+I+200),I=1,M12) C* Group resonance integral table terminated by dummy RID=0 537 IF(RID.NE. 0. ) GO TO 532 C* Check the file marks IER=53 READ (LF1,701) JTERM, K1,K2,K3,K4 IF(JTERM.NE.ITERM) GO TO 690 IER=54 WRITE(LF3,701) ITERM 550 CONTINUE C* Copy the P1 scattering matrices from the old library NP1=0 640 NP1=NP1+1 C* Identify the format on the original library (full or condensed) C* Next card may be the last filemark, matrix length or data IER=56 READ (LF1,692,END=690) REC C* Expanded matrix form if read error occurs READ (REC,701,ERR=650) KL1 C* Finish if Filemark found IF(KL1.EQ.ITERM) GO TO 680 C* Case-source lib.: copy condensed P1 matrices 642 L1 =1 READ (LF1,702) (RWO(L1-1+J),J=1,KL1) WRITE(LTT,693) ' Copy condensed P1-matrix ',NP1 WRITE(LF3,701) KL1 WRITE(LF3,702) (RWO(L1-1+J),J=1,KL1) GO TO 640 C* Case-source lib.: copy uncondensed P1 matrices (as char.strings) 650 NR =(NG+4)/5 JG =0 WRITE(LTT,693) ' Copy expanded P1-matrix ',NP1 651 JG =JG+1 WRITE(LF3,692) REC DO 652 K=2,NR READ (LF1,692) REC WRITE(LF3,692) REC 652 CONTINUE IF(JG.GE.NG) GO TO 640 READ (LF1,692) REC GO TO 651 C* All P1 scattering matrices processed - finish 680 WRITE(LF3,701) ITERM C* All processing completed CLOSE(UNIT=LF3) 683 CLOSE(UNIT=LF1) 684 RETURN C* Error condition reading the library 690 CALL ERRWAR(LTT,IER) RETURN 691 FORMAT(2A40) 692 FORMAT(A75) 693 FORMAT(A40,I8) 694 FORMAT(A40,F8.1) 695 FORMAT(BN,I10) 696 FORMAT(BN,F10.0) 701 FORMAT(5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 706 FORMAT(1P,I15/E15.8,2I6) END *DECK CHKLIB SUBROUTINE CHKLIB(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : Subroutine CHKLIB C-Purpose: Check WIMS-D library integrity CHARACTER*75 REC CHARACTER*40 BLNK, FLNM, FLWI,FLWO DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) DATA BLNK /' '/ 1 ,FLWI /'WIMS.LIB'/ 2 ,FLWO /'WIMS.LIB'/ C* Default file units DATA LF1,LF3/ 1, 3 / C* Define the filenames GO TO 12 11 WRITE(LTT,691) ' ERROR opening file : ',FLWI 12 WRITE(LTT,691) ' Default source formatted lib.filename: ',FLWI WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691) FLNM IF(FLNM.NE.BLNK) FLWI=FLNM C***** VAX ***** C OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11,READONLY) C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF1,FILE=FLWI,STATUS='OLD',ERR=11) C***** EXPORT ***** GO TO 16 15 WRITE(LTT,691) ' ERROR opening file : ',FLWO 16 WRITE(LTT,691) ' Default new formatted lib. filename: ',FLWO WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691) FLNM IF(FLWO.NE.BLNK) FLWO=FLNM C***** EXPORT ***** IF(FLWI.EQ.FLWO) GO TO 15 C***** EXPORT ***** C* C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' Source WIMS-D library filename : ',FLWI WRITE(LTT,691) ' New WIMS-D library filename : ',FLWO WRITE(LTT,691) C* C* Begin reading the library READ (LF1,695,ERR=810) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NNFPD IF(NG1+NG2+NG3 .NE. NG) THEN WRITE(LTT,691) ' Fast/Resonance/Thermal groups mismatch' GO TO 810 END IF C* C* Read MAT numbers and check for duplicate entries READ (LF1,695,ERR=810) (IWO(J),J=1,NEL) IER=0 DO 24 I=2,NEL MAT=IWO(I-1) DO 22 J=I,NEL IF(MAT.EQ.IWO(J)) THEN WRITE(LTT,693) ' Duplicate MAT number ',MAT IER=1 END IF 22 CONTINUE 24 CONTINUE IF(IER.NE.0) GO TO 810 C* C* Read the energy group boundaries and check they are monotonic NGB=NG+1 READ (LF1,696,ERR=810) (RWO(J),J=1,NGB) IER=0 DO 32 J=1,NG IF(RWO(J).LE.RWO(J+1)) THEN WRITE(LTT,693) ' Non-monotonic energy in group',J IER=1 END IF 32 CONTINUE IF(IER.NE.0) GO TO 810 C* C* Read the fission spectrum fraction READ (LF1,696,ERR=810) (RWO(J),J=1,NG0) SP=0.0 DO 42 J=1,NG0 SP=SP+RWO(J) 42 CONTINUE IF(ABS(SP-1.).GT.0.001) 1 WRITE(LTT,694) ' CHKLIB WARNING - Unnormalised spectrum ',SP C* C* Process the burnup data L1= 1+NEL L2=L1+NEL KNFD =0 KNFPD=0 DO 120 I=1,NEL READ (LF1,695) NDT ND2=NDT/2 IF( ((ND2*2).NE.NDT) .OR. 1 ((NDT.LT.8) .AND. (NDT.NE.2)) ) THEN WRITE(LTT,693) ' CHKLIB WARNING - Illegal Bu.length MAT',IWO(I) END IF READ (LF1,697,ERR=810) (RWO(J),J=1,NDT) MAT=NINT(RWO(2)) IF(MAT.NE.IWO(I)) THEN WRITE(LTT,693) ' CHKLIB ERROR - Burnup data MAT found : ',MAT WRITE(LTT,693) ' Expected : ',IWO(I) END IF IWO(L1-1+I)=2*ND2 IF(NDT.GT.8) THEN KNFD =KNFD +1 ELSE IF(NDT.EQ.8) THEN JF=NINT(RWO(8)) IF(JF.GE.2) THEN KNFD =KNFD +1 ELSE KNFPD=KNFPD+1 END IF END IF C* 120 CONTINUE IF(KNFD .NE.NNFD ) THEN WRITE(LTT,693) ' CHKLIB WARNING - NNFD Found : ',NNFD WRITE(LTT,693) ' Expected : ',KNFD END IF IF(KNFPD.NE.NNFPD) THEN WRITE(LTT,693) ' CHKLIB WARNING - NNFPD Found : ',NNFPD WRITE(LTT,693) ' Expected : ',KNFPD END IF C* Check that the next record is the library delimiter READ (LF1,695,ERR=810) JTERM IF(JTERM.NE.ITERM) THEN WRITE(LTT,691) ' Library delimiter record expected ' GO TO 810 END IF C* C* Rewrite the general information and the burnup data C***** VAX ***** C OPEN (UNIT=LF3,FILE=FLWO,STATUS='NEW',ERR=15 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LF3,FILE=FLWO,STATUS='UNKNOWN',ERR=15) C***** EXPORT ***** C* 200 REWIND LF1 C* Indices and MAT numbers READ (LF1,695) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NNFPD WRITE(LF3,695) NEL,NG,NG0,NG1,NG2,NG3,KNFD,KNFPD C* READ (LF1,695) (IWO(J),J=1,NEL) WRITE(LF3,695) (IWO(J),J=1,NEL) C* Energy group boundaries READ (LF1,696) (RWO(J),J=1,NGB) WRITE(LF3,702) (RWO(J),J=1,NGB) C* READ (LF1,696) (RWO(J),J=1,NG0) WRITE(LF3,702) (RWO(J),J=1,NG0) C* C* Nuclide decay and burnup data DO 220 I=1,NEL ND1=IWO(L1-1+I) DO 202 J=1,ND1 RWO(J)=0. 202 CONTINUE READ (LF1,695) NDT READ (LF1,697) (RWO(J),J=1,NDT) MAT=NINT(RWO(2)) C* Check the MAT numbers DO 218 J=4,ND1,2 IF(J.EQ.8) GO TO 218 MT1=NINT(RWO(J)) IF(MT1.EQ.0) GO TO 218 DO 212 K=1,NEL IF(MT1.NE.IWO(K)) GO TO 212 IF(IWO(L1-1+K).GE.8) GO TO 218 WRITE(LTT,693) ' CHKLIB WARNING - Burnup data err.MAT : ',MAT WRITE(LTT,693) ' Non-burnable product specified : ',MT1 IWO(L1-1+K)=8 IF(K.LT.I) THEN REWIND LF3 GO TO 200 END IF 212 CONTINUE WRITE(LTT,693) ' CHKLIB ERROR - Burnup data error MAT : ',MAT WRITE(LTT,693) ' Non-existent burnup product : ',MT1 218 CONTINUE WRITE(LF3,695) ND1 WRITE(LF3,703) (RWO(J),NINT(RWO(J+1)),J=1,ND1,2) 220 CONTINUE C* Next record is the file delimiter READ (LF1,695,ERR=810) JTERM WRITE(LF3,695) ITERM IF(JTERM.NE.ITERM) THEN WRITE(LTT,691) ' CHKLIB ERROR - File delimiter expected ' GO TO 810 END IF C* C* Process the nuclide cross section files JEL=0 NP1=4*NG2+2*(NG1+NG2) NP2= 2*(NG1+NG2) NP3= 2* NG3 C* Read the nuclide header record 300 READ (LF1,698,ERR=810) MAT,AA,IA,NF,NT,NR WRITE(LF3,704) MAT,AA,IA,NF,NT,NR JEL=JEL+1 MT1=IWO(JEL) IF(MAT.NE.MT1) THEN WRITE(LTT,693) ' CHKLIB ERROR - x.s. data MAT found : ',MAT WRITE(LTT,693) ' Expected : ',MT1 END IF C* Read the fast cross section data READ (LF1,696,ERR=810) (RWO(J),J=1,NP1) WRITE(LF3,702) (RWO(J),J=1,NP1) IF(NF.GE.2) THEN READ (LF1,696,ERR=810) (RWO(J),J=1,NP2) WRITE(LF3,702) (RWO(J),J=1,NP2) END IF C* Read the fast scattering matrix READ (LF1,695,ERR=810) LMX READ (LF1,696,ERR=810) (RWO(J),J=1,LMX) C* Analyse the fast scattering matrix LL=1 IT=1 NG12=NG1+NG2 IER=0 DO 322 IG=1,NG12 JER=0 LSS=NINT(RWO(LL)) LST=NINT(RWO(LL+1)) IF(LSS.LT.1 .OR. LSS.GT.LST) THEN C* WRITE(LTT,693) ' CHKLIB ERROR - Incorrect sct. matrix ',MAT WRITE(LTT,693) ' Fast group/self-sc/length ' 1 ,IG ,LSS ,LST GO TO 810 END IF C* No upscattering in the fast range is allowed C***** NOUPSCATTERING C IF(LSS.GT.1) THEN C IF(IER.EQ.0) C 1 WRITE(LTT,693) ' CHKLIB WARNING - Incorrect sct. matrix ',MAT C WRITE(LTT,693) ' Fast group/self-sc/length ' C 1 ,IG ,LSS ,LST C IER=1 C JER=1 C J1=LL+1+LSS C J2=LSS-1 C DO 312 J=J1,LMX C RWO(J-J2)=RWO(J) C 312 CONTINUE C LST=LST-LSS+1 C LSS=1 C END IF C***** NOUPSCATTERING C* Check that downscattering is within the total number of groups IF(IG+LST-1.GT.NG) THEN IF(IER.EQ.0) 1 WRITE(LTT,693) ' CHKLIB WARNING - Incorrect sct. matrix ',MAT IF(JER.EQ.0) 1 WRITE(LTT,693) ' Fast group/self-sc/length ' 1 ,IG ,LSS ,LST J1=LL+2+LST J2=IG+LST-1-NG DO 314 J=J1,LMX RWO(J-J2)=RWO(J) 314 CONTINUE LST=LST-J2 END IF RWO(LL )=LSS RWO(LL+1)=LST 321 LL=LL+LST+2 322 CONTINUE LMX=LL-1 WRITE(LF3,695) LMX WRITE(LF3,702) (RWO(J),J=1,LMX) C* C* Read the thermal data READ (LF1,696) (RWO(J),J=1,NT) WRITE(LF3,702) (RWO(J),J=1,NT) DO 348 IT=1,NT READ (LF1,696) (RWO(J),J=1,NP3) WRITE(LF3,702) (RWO(J),J=1,NP3) IF(NF.GE.2) THEN READ (LF1,696) (RWO(J),J=1,NP3) WRITE(LF3,702) (RWO(J),J=1,NP3) END IF READ (LF1,695) LMX READ (LF1,696) (RWO(J),J=1,LMX) C* Analyse the fast scattering matrix LL=1 NG12=NG1+NG2 DO 342 I=1,NG12 IG =NG1+NG2+I LSS=NINT(RWO(LL)) LST=NINT(RWO(LL+1)) C IF(LSS.LT.1 .OR. LSS.GT.LST) WRITE(LF3,99) MAT,IT,IG,LSS,LST C IF(LST.LT.1) WRITE(LF3,99) MAT,IT,IG,LSS,LST C IF(LST+1-LSS.GT.NG+1-IG) WRITE(LF3,99) MAT,IT,IG,LSS,LST C IF( IG+1-LSS.LT.NG-NG3 ) WRITE(LF3,99) MAT,IT,IG,LSS,LST LL=LL+LST+2 342 CONTINUE WRITE(LF3,695) LMX WRITE(LF3,702) (RWO(J),J=1,LMX) 348 CONTINUE C* Next record is the file delimiter READ (LF1,695,ERR=810) JTERM WRITE(LF3,695) ITERM IF(JTERM.NE.ITERM) THEN WRITE(LTT,691) ' CHKLIB ERROR - File delimiter expected ' GO TO 810 END IF IF(JEL.LT.NEL) GO TO 300 C* C* Process the resonance files IG=NG1 C* Next group 400 IG=IG+1 C* Next material 410 READ (LF1,695,ERR=810) NDT READ (LF1,699,ERR=810) RMT,NT,NS IF(NT*NS.NE.NDT) THEN WRITE(LTT,694) ' Incorrect resonance table for MAT : ',RMT WRITE(LTT,693) ' Resonance group : ',IG GO TO 810 END IF WRITE(LF3,695) NDT WRITE(LF3,706) RMT,NT,NS NP=NT+NS+NT*NS READ (LF1,696,ERR=810) (RWO(J),J=1,NP) WRITE(LF3,702) (RWO(J),J=1,NP) IF(RMT.NE.0) GO TO 410 C* Resonance group processed - next record is file delimiter READ (LF1,695,ERR=810) JTERM WRITE(LF3,695) ITERM IF(JTERM.NE.ITERM) THEN WRITE(LTT,691) ' CHKLIB ERROR - File delimiter expected ' GO TO 810 END IF IF(IG.LT.NG1+NG2) GO TO 400 C* C* process the P1 matrices NP1=0 640 NP1=NP1+1 C* Identify the format on the original library (full or condensed) C* Next card may be the last filemark, matrix length or data IER=56 READ (LF1,692,END=810) REC C* Expanded matrix form if read error occurs READ (REC,701,ERR=650) KL1 C* Finish if Filemark found IF(KL1.EQ.ITERM) GO TO 680 C* Case-source lib.: copy condensed P1 matrices 642 L1 =1 READ (LF1,702) (RWO(L1-1+J),J=1,KL1) WRITE(LTT,693) ' Copy condensed P1-matrix ',NP1 WRITE(LF3,701) KL1 WRITE(LF3,702) (RWO(L1-1+J),J=1,KL1) GO TO 640 C* Case-source lib.: copy uncondensed P1 matrices (as char.strings) 650 NR =(NG+4)/5 JG =0 WRITE(LTT,693) ' Copy expanded P1-matrix ',NP1 651 JG =JG+1 WRITE(LF3,692) REC DO 652 K=2,NR READ (LF1,692) REC WRITE(LF3,692) REC 652 CONTINUE IF(JG.GE.NG) GO TO 640 READ (LF1,692) REC GO TO 651 C* All P1 scattering matrices processed - finish 680 WRITE(LF3,701) ITERM C* All processing completed CLOSE(UNIT=LF3) CLOSE(UNIT=LF1) RETURN C* C* Error traps 810 WRITE(LTT,691) ' CHKLIB ERROR - Incorrect library format' 890 STOP 'CHKLIB ERROR' C* 691 FORMAT(2A40) 692 FORMAT(A75) 693 FORMAT(A40,3I8) 694 FORMAT(A40,F8.1) 695 FORMAT(BN,5I15) 696 FORMAT(BN,5F15.0) 697 FORMAT(3(F15.0,F6.0)) 698 FORMAT(I6,F15.0,4I6) 699 FORMAT(F15.0,2I6) 701 FORMAT( 5I15) 702 FORMAT(1P,5E15.8) 703 FORMAT(1P,3(E15.8,I6)) 704 FORMAT(1P,I6,E15.8,4I6) 706 FORMAT(1P,E15.8,2I6) END *DECK NEWLIB SUBROUTINE NEWLIB(LKB,LTT,RWO,IWO,MXRW,MXIW) C-Title : Subroutine NEWLIB C-Purpose: Make a starter file for a new library C-Description: C-D The starter file for a new library is written. Make sure that C-D the energy boundaries and the spectrum are entered with the C-D INCL command before the library is used. C- CHARACTER*40 BLNK, FLNM, FLWO REAL*8 ZERO DIMENSION RWO(*),IWO(*) COMMON /TERWOR/ ITERM,ALDO(2) DATA BLNK /' '/ 2 ,FLWO /'NEW.LIB'/ C* Default file units DATA LFO/ 2 / C* Initialise DATA ZERO,NG1,NG2,NG3,NG/ 0.D0, 14, 13, 42, 69 / C* Define the new library filename GO TO 16 15 WRITE(LTT,691) ' ERROR - Opening file : ',FLWO 16 WRITE(LTT,691) ' Default new formatted lib. filename: ',FLWO WRITE(LTT,691) '$ Enter filename to change: ' READ (LKB,691) FLNM IF(FLNM.NE.BLNK) FLWO=FLNM C***** VAX ***** C OPEN (UNIT=LFO,FILE=FLWO,STATUS='NEW',ERR=15 C 1 ,CARRIAGECONTROL='LIST') C***** VAX ***** C***** EXPORT ***** OPEN (UNIT=LFO,FILE=FLWO,STATUS='UNKNOWN',ERR=15) C***** EXPORT ***** C* GO TO 18 17 WRITE(LTT,691) ' ERROR - defining energy groups - redo ' 18 WRITE(LTT,691) '$ Enter total number of groups : ' READ (LKB,801) NGI IF(NGI.NE.0) NG =NGI WRITE(LTT,691) '$ Enter number of fast neutron groups : ' READ (LKB,801) NGI IF(NGI.NE.0) NG1=NGI WRITE(LTT,691) '$ Number of resonance groups : ' READ (LKB,801) NGI IF(NGI.NE.0) NG2=NGI NG3=NG-NG1-NG2 IF(NG1.LE.0 .OR. 1 NG2.LE.0 .OR. 2 NG3.LE.0) GO TO 17 C* Verification printout WRITE(LTT,691) WRITE(LTT,691) ' Processing: ' WRITE(LTT,691) ' New WIMS-D library filename : ',FLWO WRITE(LTT,692) ' Number of fast neutron energy groups : ',NG1 WRITE(LTT,692) ' Number of resonance energy groups : ',NG2 WRITE(LTT,692) ' Number of thermal energy groups : ',NG3 WRITE(LTT,692) ' Total number energy groups : ',NG WRITE(LTT,691) C* C* Begin writing the library starter file NEL=0 NG0=NG1 NNFD=0 NNFPD=0 C* Library header index record WRITE(LFO,695) NEL,NG,NG0,NG1,NG2,NG3,NNFD,NNFPD C* No MAT numbers C WRITE(LFO,695) C* Energy group boundaries NGB=NG+1 WRITE(LFO,696) (ZERO,J=1,NGB) C* Fission spectrum fraction WRITE(LFO,696) (ZERO,J=1,NG0) C* File marks before the nuclide files WRITE(LFO,695) ITERM C* File marks before the resonance files DO 30 I=1,NG2 WRITE(LFO,695) ITERM WRITE(LFO,695) 1 WRITE(LFO,697) ZERO,1,1 WRITE(LFO,696) ZERO,ZERO,ZERO 30 CONTINUE C* File mark before the P1 matrix WRITE(LFO,695) ITERM C* Dummy P1 matrices DO 40 I=1,4 DO 38 IG=1,NG WRITE(LFO,696) (ZERO,JG=1,NG) 38 CONTINUE 40 CONTINUE C* End-of-library filemark WRITE(LFO,695) ITERM CLOSE(UNIT=LFO) C* 691 FORMAT(2A40) 692 FORMAT(A40,I5) 695 FORMAT(5I15) 696 FORMAT(5F15.12) 697 FORMAT(F15.0,2I6) 801 FORMAT(BN,I10) END