U.S. DEPARTMENT OF COMMERCE / National Bureau of Standards

emiconductor Measurement 1 echnology:

DISTRIB 1,1

An Impurity Redistribution]

Computer Program!

100

U57
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards' was established by an act of Congress March 3, 1 90 1 . The

Bureau's overall goal is to strengthen and advance the Nation's science and technology and

facilitate their effective application for public benefit. To this end, the Bureau conducts

research and provides: (1) a basis for the Nation's physical measurement system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau's technical work is performed by the National Measurement Laboratory, the National Engineering Laboratory,

and the Institute for Computer Sciences and Technology.

THE NATIONAL MEASUREMENT LABORATORY provides the national system of

physical and chemical and materials measurement; coordinates the system with measurement

systems of other nations and furnishes essential services leading to accurate and uniform

physical and chemical measurement throughout the Nation's scientific community, industry,

and commerce; conducts materials research leading to improved methods of measurement,

standards, and data on the properties of materials needed by industry, commerce, educational

institutions, and Government; provides advisory and research services to other Government

Agencies; develops, produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

Absolute Physical Quantities' — Radiation Research — Thermodynamics and

Molecular Science — Analytical Chemistry — Materials Science. THE NATIONAL ENGINEERING LABORATORY provides technology and technical

services to users in the public and private sectors to address national needs and to solve national problems in the public interest; conducts research in engineering and applied science

in support of objectives in these efforts; builds and maintains competence in the necessary

disciplines required to carry out this research and technical service; develops engineering data and measurement capabilities; provides engineering measurement traceability services; develops test methods and proposes engineering standards and code changes; develops and

proposes new engineering practices; and develops and improves mechanisms to transfer

results of its research to the utlimate user. The Laboratory consists of the following centers:

Applied Mathematics — Electronics and Electrical Engineering^ — Mechanical

Engineering and Process Technology^ — Building Technology — Fire Research —

Consumer Product Technology — Field Methods.

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts

research and provides scientific and technical services to aid Federal Agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and

economy in Government operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing the Federal Information Processing Standards Program, developing Federal ADP standards

guidelines, and managing Federal participation in ADP voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal Agencies; and

provides the technical foundation for computer-related policies of the Federal Government.

The Institute consists of the following divisions:

Systems and Software — Computer Systems Engineering — Information Technology.

'Headquarters and Laboratories at Gaithersburg, Maryland, unless otherwise noted; mailing address Washington, D.C. 20234. ^Some divisions within the center are located at Boulder, Colorado, 80303.

The National Bureau of Standards was reorganized, effective April 9, 1978.
MAR 1 5 i979

Semiconductor Measurement Technology: s^c \oo

DISTRIB I, An Impurity Redistribution Computer Program

^^'^

\Qn^

David Gilsinn

Institute for Computer Sciences

and Technology

National Bureau of Standards

Washington, D.C. 20234

and

Richard Kraft

Center for Applied Mathematics

National Engineering Laboratory

National Bureau of Standards

Washington, D.C. 20234

This activity was supported by

The Defense Advanced Research Projects Agency

1400 Wilson Boulevard

Arlington, Va 22209

U.S. DEPARTMENT OF COMMERCE, Juanita M. Kreps, Secretary

Jordan J. Baruch, Assistant Secretary for Science and Technology

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Issued February 1979
Library of Congress Catalog Card Number: 79-600002

National Bureau of Standards Special Publication 400-57

Nat. Bur. Stand. (U.S.), Spec. Publ. 400-57, 130 pages (Feb. 1979) CODEN. XNBSAV

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1979

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402

Stock No 003-003-02030-3 Price $3

(Add 25 percent additional for other than U.S. mailing).
PREFACE

This work was conducted as part of the Semiconductor Technology Program

at the National Bureau of Standards and was funded by the Defense Advanced Research Projects Agency (Order No. 2397, Program Code 5D10)

.

The authors wish to thank Sally Peavy for permission to include in Appen

dix 2 the documentation of her plot routine, W. R. Thurber for the exper

imental data shown in figure 11, and also acknowledge the efforts of M.

G. Buehler in the formulation of the redistribution problem; finally,

the authors wish to thank Bill Keery for his editorial aid and Jane Walters for the typing and layout work.

iii

I

TABLE OF CONTENTS

Page

Preface iii

Abstract 1

1. Introduction 1 1.1. Main Objectives of the Documentation 1

1.2. Tips on Where to Locate the Key Parts of This

Documentation 1

2. Description of the Problem 4

2.1. Physical and Mathematical Description of the

Problem 4

2.2. The Partial Differential Equations in Their

Equivalent Integral Form 9

3. The Computational Procedure 11

3.1. General Description of the Computational Procedure . 11

3.2, Short Guide to the Following Chapters 14

4. Description of Auxiliary Computational Procedures 15

4.1. A Detailed Description of the Finite Difference

Grid 15

4.2. The Description of the Moving Boundary 16

4.3. Dimensionalization 18

4.4. The Description of the Way Arrays CI and C2 Are

Plotted 18

5. Derivation of the Discrete Algebraic Equations 20

5.1. The General Prescription for Obtaining the Tridiagonal Equation Associated with a Fixed Grid

Point 20

5.2. The Index of the Z-Grid Point Ranges from 2 to

NL - 1 When NL > 2 23

5.3. The Index of the Y-Grid Point is Greater than NR . . 26

5.4. The Index of the Z-Grid Point is NL and NL > 1 ... 28

5.5. The Index of the Y-Grid Point is NR 30

5.6. The Grid Points Lie at the Moving Boundary and

NL > 1 33

5.7. The Grid Point Is at the Moving Boundary and NL = 1 . 38

5.8. The Index of the Z-Grid Point is 1 and NL > 1 . . . . 40

5.9. The Index of the Z-Grid Point is 1 and NL =^ 1 . . . . 42

6. Method of Solving the Algebraic Equations 45

6.1. The Derivation of the "Double Sweep" "P, Q" Coefficients Associated with the Grid Points 45

6.2. The First Z-Grid Point VThen NL = 1 and tjyiM = TTMM =0 46

6.3. The First Z-Grid Point When NL = 1 and t^^ = TTMI4 > 0 45

6.4. The Index of the Z-Grid Point is 1 and NL > 1 . . . . 47

6.5. The Index of the Z-Grid Point is Between 1 and NL - 1 and NL > 2 47

6.6. The Index of the Z-Grid Point is NL and NL > 1 ... 48

V
Page

6.7. The Index of the Y-Grid Point is N2 = NYO + 1 . . . . 49

6.8. The Index of the Y-Grid Point is Between NR and N2 . 49

6.9. The Y-Grid Point has Index NR 50

6.10. The Determination of CB(1, TTMMPl) in Terms of

Known Quantities 51

7. Detailed Description of the Program 53

7.1. Block-by-Block Detailed Description of the Calculational Procedure in the Main Program 53

7.2. Block-by-Block Description of Subroutine CENTER ... 62

7.3. Block-by-Block Description of Subroutine TRIDNL ... 63

7.4. Description of Subroutine INTER 64

8. Flow Diagrams 65

8.1. Flow Diagram of Main Program in DISTRIE 65

8.2. Flow Diagram for Subroutine CENTER 78

8.3. Flow Diagram for Subroutine TRIDNL 80

9. How to Use DISTRIB 81

9.1. A Brief Description of DISTRIB: VERSION 1 81

9.2. Definitions and Other Information Related to Input

and Output 81

9.2.1. Real Input 82

9.2.2. Integer Input 84

9.2.3. The Real Output of DISTRIB 85

9.2.4. The Integer Output 86

9.3. A Simple Example 86

10. Verification of Program 94

10.1. Validation of the Program ..... 94

10.2. Limitations of the Program 94

11. Program Listings

.11.1. DISTRIB: VERSION 1 Main Program 95

11.2. Subroutine CENTER 105

11.3. Subroutine TRIDNL 107

11.4. Subroutine INTER , 108

12. Glossary 109

13. References 114

13.1. Bibliography 114

13.2. Supplemental Bibliography ]_]_^ Appendix 1 116

Appendix 2 118

vi
LIST OF FIGURES

Page

1. The problem geometry 5

2. A hypothetical redistribution process 6

3. The moving z-coordinate frame 7

4. Finite difference grids with a = 1/3 and NZO =3 12

5. Motion of the moving boundary in the z- region 16

6. Motion of the moving boundary in the y-region 15

7. Flux cells and flux boundaries in the y-region for fixed

grid points 20

8. Flux cells for fixed grid points in the z-region when

NL > 2 21

9. Flux cells for fixed grid points in the z-region when

NL = 2 22

10. Schematic description of the flux cell associated with the

discretization of the conservation equation across the moving boundary 23

11. Plot of example concentrations in oxide when TTMMPl ~ TFINAL 91

12. Plot of example concentrations in silicon when TTMMPl ~ TFINAL 92

13. Comparison of three calculated impurity concentration distributions in silicon with experimentally determined con- centration data 93

vii
LIST OF TABLES

Page

1. FORTRAN Version of Symbols in Eqs (3.2) to (3.4) in the

Y- and Z-Regions 12

2. FORTRAN Version of Symbols Used in the Moving Boundary

Equation and Finite Difference Grid 17

3. Guide to the Sections Where the Tridiagonal Algebraic

Equations Associated with the Grid Points Are Described

and Derived . 23

4. FORTRAN Version of Quantities in Eq (5.5) 25

5. Intermediate Quantities Used in Eq (5.7) 25

6. FORTRAN Version of Quantities in Eq (5.8) 26

7. Intermediate Quantities Used in Eq (5.10) 27

8. FORTRAN Version of Quantities in Eq (5.16) 30

9. FORTRAN Version of Quantities in Eq (5.24) 32

10. FORTRAN Symbols Used in Eq (5.28) . . 34

11. FORTRAN Version of Symbols in Eq (5.33) 37

12. Additional FORTRAN Symbols for Eq (5.38) 40

13. FORTRi^N Version of Quantities in Eq (5.41) 42

14. FORTRAN Version of Quantities in Eq (5.46) 44

15. Guide to the Subsections Where the "P, Q" Coefficients

Associated With Each of the Fixed Grid Points Are Derived

and Described * 45

16. Brief Definitions of the Real Input 82

17. Real Input 83

18. Brief Definition of Integer Input ... 84

19. Integer Input Information 84

20. Constants Used in Program But Not Made Part of Input .... 85

21. Brief Definitions of the Real Output 85

22. Real Output 86

23. Definitions of Integer Output 86

24. Example of Output of Concentration in Oxide and Silicon at

End of Iteration RM = (1 + L*KPRINT) in Loop 78 for L = 0 . 90

25. Example Concentrations in Oxide and Silicon When TTMMPl = TFINAL 90

26. Listing of Input Physical Data and Important Computational

Data 90

viii
Semiconductor Measurement Technology

:

DISTRIB I, An Impurity Redistribution Computer Progrcun

by

David Gilsinn and Richard Kraft

This report provides docxjmentation of a computer program

which calculates the redistribution of impurities in

silicon during a single oxidation step. The dociomentation provides: (1) a physical and mathematical description of the redistribution process, (2) a detailed description of the discretization of the appropriate partial differential equations, and (3) a complete description of the FORTRAN program for computing the solution.

Key Words: Diffusion; electronic technology; impurity

distribution; material transport; segregation; semiconductor technology.

1 . INTRODUCTION

1.1. Main Objectives of the Documentation

The principal objectives of this documentation are:

(1) to give a detailed derivation of the algebraic equations that

are used in the computer program,

(2) to give in detail, a step-by-step description of the way the

computation in each computational block of the progreun is carried^ out, and

(3) to provide a simple example illustrating how to use the program.

The format of this documentation is modeled after the partial differential

equation computer prograjn format description in [9] in the FORTRAN IV language .

In fabrication processes redistribution of impurities takes place in multiple oxidation and diffusion steps (each step being characterized by con- stant process parameters) . DISTRIB is aimed at simulating this complete

process; however, due to the complexity of the total problem, it was decided to only document a fundamental oxidation step here. This ejcplains

the appellation DISTRIB VERSION 1.

The next section introduces the user to DISTRIB VERSION 1.

1.2. Tips on Where to Locate the Key Parts of This Documentation

1
Question 1:

Answer 1:

Question 2

Answer 2:

Question 3:

Answer 3

Question 4:

Answer 4:

Question 5:

Answer 5

Question 6:

Answer 6:

Question 7:

Answer 7:

Question 8:

Answer 8:

Question 9:

Answer 9:

Question 10;

Answer 10:

Question 11:

Answer 11:

Question 12:

Answer 12:

Question 13:

What should I do if I have a question about this program

that is not answered in the documentation?

Call Richard Kraft, NBS, Applied Mathematics Division or Electron Devices Division, (301) 921-3621.

Where is the physical problem explained?

In section 2.1.

Where are the partial differential equations listed in sum- mary form?

Equations (2.13) to (2.19).

Where are the computational procedures outlined?

In section 3.1.

Where are the discretized equations derived?

In chapters 5 and 6

What machine has the program run on?

UNIVAC 1108.

Where are the flow diagrams?

In chapter 8.

Where are the detailed explanations of the computational

steps?

In chapter 7.

Where are the input and output and a description of how to

work a simple problem?

In chapter 9.

Where is the listing of the computational blocks?

In chapter 11.

Suppose I want to know where to go in the documentation to

learn about some quantity in the listing?

Go to the glossary in chapter 12.

What is the basis for the program's validity and what are

its limitations?

See chapter 10.

What about the "deck cards" and the statements in the listings that are not dociomented in VERSION 1?
Answer 13: Completely ignore them. Everything needed to understand

VERSION 1 is explained in this documentation and is in the

computational blocks (see listing, chapter 11)

.

3
2. DESCRIPTION OF THE PROBLEM

2.1. Physical and Mathematical Description of the Problem

The problem considered here [4] deals with the redistribution of dopant

impurities found in silicon as the silicon is oxidized and converted

into silicon dioxide. For the sake of specificity, we consider the impurity to be boron and have studied the following ideal situation.

Suppose the entire half space lying to the right of the x-z plane, see

figiire la, to be occupied by silicon while the region on the left-hand

side of this plane is filled with oxygen at some elevated temperature.

Also, suppose that the silicon is doped with boron and that the density

distribution of this impurity is constant in planes parallel to the x-z plane so the distribution is essentially only a function of the ydimension.

The time evolution of an initially uniform boron distribution (for appropriate conditions at the oxygen-oxide interface) is illustrated in figure

lb. When the oxygen contacts the silicon, it begins to react to form

silicon dioxide, and a sharply defined moving boundary which separates

the pure silicon from that previously converted into silicon dioxide begins to move into the pure silicon. The position of this moving boundary

is known empirically, see [1] and section 4.2, and is denoted by

where t is time.

Inside the silicon and oxide, the transport of boron is governed by or- dinary Fickian diffusion but with different diffusion coefficients Dj

and D2 in the oxide and silicon, respectively.

The more complicated physical processes take place at the moving oxidesilicon interface. The boron has a preference to be dissolved in the

oxide rather than the silicon. Hence, as the moving boundary moves into

the silicon, the boron's preference to be in the oxide depletes the boron density in the vicinity of the moving boundary. The preference of

the boron to be in the oxide rather than the silicon is quantitatively

expressed by requiring that the moving boundary segregation condition be

satisfied.

y = y (t) , for t > 0 (2.1)

C2(y (t) ,t) =

o m Ci (y^(t) ,t) (2.2)

where m is a given constant (the segregation coefficient) and eq (2.2)

must be satisfied for all time t > 0. In eq (2.2) Ci(y,t) and C2(y/t)

represent the boron densities in the oxide and silicon, respectively.
X

MOVING BOUNDARY

0- SiO:

I 1 Y

Yo(t)

Figure 1. The problem geometry.

In addition to the condition eq (2.2) holding at the moving boundary, it

is also assumed that there is conservation of boron atoms. Mathematically, this condition is expressed by the requirement that

i y^ + (a-^-l)y Ci

y=yjt)

= DoS Co

y

^ y=y^(t)

(2.3)

y=y^(t)

In order to understand the origin of the convective flux (a~ '-1) yQCj , it

is necessary first to understand a second complicated physical effect occurring at the moving oxide-silicon boundary. When silicon is converted

into oxide, its volume expands by a factor a~^, where a 0.45. In the

model chosen here [4] , it is assumed that all the oxidation of the silicon takes place entirely at the moving oxide-silicon interface. In this

case, the expansion tends to push the previously formed oxide backwards

at a certain velocity into the oxygen, creating another moving boundary

between the oxygen and oxide, see figure 2. The velocity with which the

oxide is being convected backwards is assumed to set up a convective

drift transport in the oxide's boron atoms. The magnitude of the velocity of this convective flux of boron atoms is determined in the following

manner

.

5
BORON REDISTRIBUTION DURING OXIDATION

TIME

0

OXYGEN SILICON

+ 1

+ 2

OXYGEN

+ 3

OX DE SILICON

GOVERNS ELECTRONIC

PROPERTIES OF SILICON

DEVICES

Figure 2. A hypothetical redistribution process,

Consider a thin slab with width Ay of pure silicon located exactly at

the moving oxide-silicon interface at some time just before it is con- verted into oxide. After some time At, the moving boundary has traveled from one side of this thin slab to the other

Ay = y (t + At) - y (t ) , o o o o

(2.4)

converting it into a thin slab of pure oxide of width a~-^Ay. The difference between this new larger width and the original width is assumed to

be the distance that the previously formed oxide (i.e., prior to t^) is

moved backward. Therefore, this rate is

lim g-^Ay - Ay

At->0 At

_

^ ' ^

dt

(2.5)

t=t

and the backward velocity of the oxide is -(a ^-Dy^ where y^ = dt

Thus, the full equation describing the transport of the boron in the oxide is

6
) Ci(y,t) = Di9"Ci + (a-l-l)y (t)8 Ci , (2.6)

where the last term in eq (2.6) represents the convective boron transport.

By introducing a coordinate system (see figure 3) , moving with the origin at the oxygen-oxide interface,

z = y + (a-^-l)y (t) o

or (2.7)

z = y + z^(t) - y^(t) ,

where

z (t) = a-^y (t) , o o

(2.8)

eq (2.6) simplifies to

3 Ci = D.a^Ci for 0 < z < z (t) . t-^ -^z-^ — — o

(2.9)

Moreover, in the z-frame there is no longer a moving boundary at the

oxygen-oxide interface to contend with.

The moving boundary conservation condition in the z- and y-frame is

Di3 Ci(z ,t) + z Ci(z ,t) = D29 C2(y ,t) + y C2 (y ,t) . (2.10)

^ z ^ o 0^0 y o 0 0

In this model we assume the transport of boron across the oxygen-oxide

interface to be governed by

Di 9 Ci = C (Ci (0,t) - C J ,

i-~> out

(2.11)

OXYGEN-^ OXIDE- •SILICON

Figure 3. The moving z-coordinate frame.

7
where Cp and Cout assumed to be known constants. Finally, at y = oo we assume the concentration to be stationary and equal to a "bulk" value

B

C2(~.t) = Cg . (2.11')

Since we intend to solve the above equations on the computer, it is mandatory to truncate the infinite y-domain by selecting some point y = BDRY

far from the y-origin and requiring there that

C2(BDRY,t) = C . B

(2.12)

Finally, at time t = 0 it is assumed that the silicon slab is completely

unoxidized [y (0) = 0] and that the concentration has a step function distributxon

C2(y,0) = C

C2(y/0) = C

max y <_ NH

NH < y < BDRY

(2.13)

where C , NH > 0 are parameters,

max

In summary, we will solve eq (2.13) and

a Ci = Di3 ^Ci for 0 < z < z (t) , t ^

i 2 ^ — — o

) Co = Do9 ^Co for y (t) < y < BDRY , t^ Y o — —

Di3 Ci + z Ci

^ z ^ o

^ = DoB Co + y Co

z=z (t) o y=y (t) o

(2.14)

(2.15)

(2.16)

C2 (y /t) = mCi (z ,t) , o o

Di3 Ci - C (Ci - C ) ,z = 0 , and

(2.17)

(2.18)

(2.19)

z=0

C2(BDRY,t) = ,

where the moving boundary yr,(t) or z (t) = a~V (t) is a specified monotonic curve, see section 4.2.

8
2.2. The Partial Differential Equations in Their Equivalent Integral

Form

The usual procedure for the numerical solution of partial differential

equations like eqs (2.13) to (2.19) comes through solving the system of

algebraic equations that results from replacing the partial derivatives

in these equations by their approximating difference quotients. When

this standard procedure was tried (see [10-13] for other approaches) , it

did not lead to a very accurate solution. The discretization of the following system of equations equivalent to the partial differential equations and which involves integrals in place of partial derivatives does

lead to a system of algebraic equations that can be solved readily to

yield an accurate solution of the original system.

The system of equivalent equations is

:

Ci (z,t)dz = -Di3^Ci + Di9 Ci

^ z

^

z=z. = 2B

(2.20)

for every pair of points z^ < Zg such that 0 <_ z^, Zg <_ z^{t) ,

d_

dt /

C2(y,t)dy = -029^02 + Co

^ y

^

y=yA y=yB

(2.21)

for every pair of points y^^ < yg such that yo(t) f_ y^/ Y-q — BDRY,

Ci(z,t)dz t) + // C2(y,t)dy = -Di9^Ci

yo(t)

+ DoS Co

y

Z=Z7

, (2.22)

y=yB

C2(y ,t) = m Ci (z ,t) , o o

DiS Ci - C (Ci (0,t) - C ^) = 0 , z p out

C2(BDRY,t) = Cg , and

(2.23)

(2.24)

(2.25)

C2(y/0) = C^(y) (2.26)

9
where

C^(y) =

'max y <_ NH

:^ NH < y < BDRY .

To see that eq (2.14) is equivalent to eq (2.20) , first cariry out the

differentiation on the left-hand side of eq (2.20) under the integral

sign and substitute for 9^Ci by employing eq (2.14) . One gets

d_

dt I

ci (z,t)dz = / D^a^^c^ dz (2.27)

and

/

DiS^ Cldz = -Di9^Ci

z=z. z=z.

(2.28)

yielding the right-hand side of eq (2.20).

By reversing these steps starting with the right-hand side of eq (2.20), one finds

/

B

O^Ci - Di9 2c,)dz = 0

t ^ ^ z

^ (2.29)

for every 0 1. < Zg <_ ZQ(t) . Therefore, the integral of eq (2.29) must

vanish (it is assumed to be continuous) and eq (2.14) results. By using

trivial modifications of these same techniques, the remaining equations

can be shown to be equivalent to eqs (2.13) to (2.19).

10
3. THE COMPUTATIONAL PROCEDURE

3.1. General Description of the Computational Procedure

The numerical approach to the solution of the partial differential equations in the integral form, eqs (2.20) to (2.26), replaces the problem of

finding the continuous functions Ci(z,t), C2(y,t) with the problem of determining the functions Ci(z, tj^) , 02(7, tMM) which are only defined at

discrete points on a finite difference grid. The time coordinates of

these grid points are specified through the definition of a sequence of

discrete times for MM = 1 , MME where MME is some given positive

integer and the first discrete time, t^, is zero. The spatial coordinates

of the "fixed" grid points (see sec. 4.1 and fig. 4 for further details)

are

z = (J - l)Az for J = 1,2,... , (3.1)

y = (I - l)Ay for I = 1,...,N2 ,

where N2 is some given integer index. The mesh widths Az, Ay are also

given quantities whose magnitudes like those of MME and N2 need not be

considered now. In addition to the fixed grid points in eq (3.1) , there

are two grid points located at the moving boundary (see sec. 2.1). One

is at the moving boundary bordering the z-region and the second is at

the moving boundary bordering the y-region , see figure 4. The space

and time coordinates of these two moving boundary grid points are, re- spectively, (Zq ( tMM) ' ^mm) (yo(tMM) '^MM) • addition to seeking the

concentrations at the fixed grid points C^ (J,tMM) / C2(I,tMM)/ we also desire to find the concentrations at the two moving boundary grid points,

^1 (^^^^^mm) '^MM^ ' C2(y^(t ),t ) forMM=l,..., MME.

o MM MM o MM MM

Because of the initial condition eq (2.26), these concentrations are known at time t^ =0.

Let us assume that we have already determined the concentrations at the

fixed and moving grid points up to some arbitrary time t^n; then we will

proceed to describe how these same concentrations at the next time t^M + At = t^M+i (where At = tj^^j. ~ ''^MM^ ^® found by performing the operations in the following three steps.

Step #1: Discretize the equations.

By discretizing the equivalent description of the partial differential

equations in eqs (2.20) to (2.26) , we derive a system of linear algebraic tridiagonal equations that couple the concentrations at time tj^

with those at t^M+l* Specifically, with each of the fixed grid points

in the z-region, we derive and associate an equation of the form

ZiCl(J - 1,TTMMP1) + Z2CI (J,TTMMP1) + ZsCKJ + 1,TTMMP1) = Z14 (3.2)

11
Yaxis

Z axis

1/3 and

Table 1: FORTRAN Version of Symbols in Eqs (3.2) to (3.4) in the Yand Z-Regions.

Cl((P - l)Az,t^_^^) = CI (P, TTMMPl)

= C2(P, TTMMPl)

t = TTMMPl

MM+1

o MM+1 MM+1

= CB(1, TTMMPl)

'^2(yo^Wi^'Wi^ E CB (2, TTMMPl)

POSITION OF MOVING BOUNDARY

-H AyhI 2 3456 789 N2

Y = 0

Z = 0

MOVING BOUNDARY

GRID POINTS

-t 1- 2 3 4 5 6 7 8 9 10

Figure 4. Finite difference grids with a NZO = 3.

for J = 1,..., NL and where the FORTRAN notation is defined in table 1.

Similarly, with each of the fixed grid points in the y-region, we derive

and associate an equation of the form

YlC2(I - 1, TTMMPl) + Y2C2 (I, TTMMPl) + Y3C2(I + 1, TTMMPl) = Y^ (3.3)

where I = NR, . . . , N2 and table 1 contains the FORTRAN symbol definitions.

The coefficients Zi,...,Zn, Yi , . . . ,Y^. are composed of the mesh width

quantities At, Az, Ay and hence are known. In addition, the quantities

Zi^,Yi^ contain, respectively, (z , tj^^jyj) and C2(y,tj^^) and therefore these

equations couple the concentrations at the two successive time steps

t and t ,

MM MM+1

At the moving boundary grid point in the z-region, there is an additional equation called the preliminary equation (for reasons explained in

step #3) which expresses the value of (zq (tf^+j) , tj^+l) in terms of

the concentrations at the neighboring grid points, at time t^+i.

12
DM1 CB{1,TTMMP1) = DM2 CI (NL,TTMMP1)

+ DM3 C2 (NR,TTMMP1) + DM4 CI (NL - 1,TTMMP1) (3.4)

+ DM5 C2 (NR + 1,TTMMP1) + DM6

where the FORTRAN symbols are defined in table 1. In this equation

DM1,..., DM6 are composed of known mesh quantities and DM6 contains con- centrations, at time tj^^ (at the grid points coincident and adjacent to

the moving boundary grid points) , which are of course assumed to be

known

.

Step #2: Put the tridiagonal equations into "P,Q" form.

By algebraically manipulating eq (3.2), it is possible to put it in

the f orm

C1(J,TTMMP1) = P1(J) C1(J + 1,TTMMP1) + Ql (J) (3.5)

for J = 1,..., NL - 1 and

CI (NL,TTMMP1) = PA(1) CB(1,TTMMP1) + QA(1) . (3.6)

By similar manipulations, eq (3.3) can be put into the form

C2(I,TTMMP1) = P2(I) C2(I - 1,TTMMP1) + Q2 (I) (3.7)

for I = NR + 1,...,N2 and

C2 (NR,TTMMP1) = PA(2) CB(2,TTMMPl) + QA(2) . (3.8)

The "P, Q" coefficients in these equations, i.e., PI ( J) , Ql ( J) , P2(I),

Q2(I), PA(L), QA(L) for L = 1,2 are algebraic combinations of the quantities Y^,..., Yi^, Z^,..., Zi^ and are hence known.

Step #3: Determination of the concentrations at time TTMMPl = tMM+1

The first important use of eqs (3.5) to (3.8) is that it is possible

using these equations with J = NL - 1 in eq (3.5) and I = NR + 1 in

eq (3.7) to eliminate with the help of eq (2.23) all the unknown concentrations on the right-hand side of the preliminary eq (3.4). Thus, we

obtain the value of the concentration at the moving boundary in the z- region in terms of known quantities,

CB(1, TTMMPl) = "Some function of (DMl , , , . ,DM6 , Zi,...,Yi^)" . (3.9)

Since eq (3.4) leads directly in this way to what will turn out to be

the most important concentration at time TTMMPl in the sense that with

this value of CB(1, TTMMPl) all the other concentrations are simply calculated, we call eq (3.4) the preliminary equation for distinction pur13
poses. From this result and the segregation boundary condition eq

(2.23) with SEG = m, we can find DB (2 ,TTMMP1)

CB(2,TTMMP1) = SEG • CB(l,TTMMPl) . (3.10)

By using these two concentrations at the moving boundary and back substituting in eqs (3.6) and (3.8), we can find C1(NL,TTMMP1) and

02 (NR,TTMMPL)

Finally, by using these two quantities and by back substituting in eq

(3.5) starting with J = ML - 1 and by back substituting in eq (3.7) beginning with I = NR + 1 , we can successively generate all the concentrations at the fixed grid points at time t . . TTMMPl

MM+i

By repeating steps one through three starting at time tj^ = 0 when the

concentrations are known in the y-region (there is no z-region since

the moving boundary starts at z = 0 at time t - 0) , all the concentrations at all the grid points at the discrete times t , MM = 1,...,MME MM

can be found.

3.2. Short Guide to the Following Chapters

In this documentation, the derivation of the tridiagonal linear eqs

(3.2), (3.3), and (3.4) is carried out in Chapter 5. The manipulation

of these equations into the form eqs (3.5) to (3.9) is carried through

in Chapter 5. The computational procedure is described in detail in

Chapter 7, and Chapter 8 contains flow diagrams of the listing in Chapter 11. Chapter 9 contains a worked example and a compact and precise

description of the input and output for the program. Chapter 10 describes the validation and limitations of the program, and Chapters 12

and 13 contain a glossary for FORTRAN symbols and references. The appendices collect finite difference approximations used in Chapter 5 and

an "in-house" plot documentation.

The next chapter collects material needed in Chapter 5 and gives a detailed description of the finite difference grid, moving boundary motion, equation dimensionalization , and plotting description.

14
4. DESCRIPTION OF AUXILIARY COMPUTATIONAL PROCEDURES

4.1. A Detailed Description of the Finite Difference Grid

In this section we describe more fully the finite difference grid already introduced in section 3.1, see eq (3.1). We (see fig. 4) begin

by describing how the magnitudes of the mesh widths Ay, Az, At and the

number of grid points, N2, on the y-axis are determined.

The mesh width Ay is calculated in DISTRIB by dividing the program

input quantity BDRY [see above eq (2.12)] by an integer input NYO, i.e.,

BDRY

Ay = . (4.1) ^ NYO ^ '

Consequently, placing the grid point index 1 at y = 0 and calling the

index of the y-grid point at y = BDRY, N2 , implies that

N2 = NYO + 1 . (4.2)

In order to understand the motivation for the way we define the magnitude of the mesh width Az, imagine the width of an interval of size Ay

of silicon after it has been converted into oxide. It swells to a size

[see sec. 2,1] a~-^Ay. Since in DISTRIB a is an unrestricted input and

is possibly very small, this latter quantity, a~-'^Ay, could be very large.

We want to choose the magnitude of Az in such a way that there are enough Az intervals in this width to describe any concentration variation in it. Therefore, we introduce another program input integer NZO

and use it to define

Az = a" ^ Ay/NZO (4,3)

in such a way that a~VNZO % 1 so that Az and Ay are approximately the

same magnitude. Another very important advantage of this choice of Az

is that it implies that when the moving boundary sweeps through a mesh

width of magnitude Ay in the y-region NZO mesh widths of magnitude Az are swept out in the z-region.

In order to explain the reasoning behind the definition of At, we re- mark that the computer program's numerical stability seems to be affected if the moving boimdary moves too rapidly between successive grid

points. Therefore, At is defined in such a way that the moving boundary

does not travel more than a fraction WFRC (real program input) of a Az

mesh width between the times t = 0 and t = At. Thus, we solve the

equation

WFRC • Az = z (At) - z (0) (4.4) o o

for At and then define At as this solution. This is how At for the

first iteration (see sec. 3.1) when t^ = 0 and Zq(0) = 0 is chosen. In

computational blocks #16 and #20 (see sec. 7.1), it is explained how At

is increased during later iterations of the main loop 78 in DISTRIB.

15
time= t^iM + i ^

time = tMM •

NLMI NL

SM(I)AZ

MOVING BOUNDARY

GRID POINTS

NLMI NL

Figure 5 : Motion of the moving boundary in the z-region

MOVING BOUNDARY

GRID POINTS

NR NRPI

S(2)Ay

SM(2)Ay

N2

time = twMf I

N2

time = tMM

Figure 6. Motion of the moving boundary in the y-region.

Before leaving the description of the finite difference grid, it is important and convenient to introduce special names for certain indices

of special grid points. Figure 5 shows the motion of the moving boundary starting at time t|yjj,,j and ending at time t]\^+i in the z-region. We

denote the index of the z-grid point which is at or immediately to the

left of the moving boundary at time by NL and the next grid point

to its left by NLMI (= NL - 1) , Similarly (see fig. 6) , the index of

the first grid point in the y-region greater than the position of the

moving boundary at time tj^^ is denoted by NR, and the index of the next

point to its right is NKPl (= NR + 1) . Of course, the values of NL and

NR change with t,„ . MM

4.2. The Description of the Moving Boundary

In DISTRIB the user can choose between two different formulas for the

moving boundary. These are the two motions related to the wet and dry

oxidations in [1]. In both cases there are three parameters A, B, and

T (each with values that are temperature dependent) in the equations

that describe the motion of the moving boundary.

Case I: Wet oxidation, t = 0.

The equation is [see 1, eq (13)]

A r 4tBll/^

- 1 for t > 0 (4.5)

16
The equation expressing t in terms of is

t = (z^^ + Az^)/B . (4.6)

Table 1 in [1] gives the values of A and B for different temperatures.

Case II: Dry oxidation, t > 0.

For times t > x [1]

/2

= - I 1 + ipr-\ ~ ^ ^ >. - ('i-'^)

When t < T , we use [1]

A r 4tBl^'

where

ZCRIT ^ z = — • t , 4.8

O T

jl/2

ZCRIT = CNSTX = ^ 11 + - 1 . (4.9)

A r 4tBT

CNSTX is a duplicate name for ZCRIT. Table II in [1] gives the values

of A, B, and T for different temperatures.

The FORTRAN symbols for the constants in the previous equations are given in table 2.

Table 2: FORTRAN Version of Symbols Used in the Moving Boundairy Equation and Finite Difference Grid.

A = CNSTA

B = CNSTB

. 5A = CNA

T = CNTAU = TLIN

ZCRIT = CNSTX

Az = DZ

Ay = DY

At = DT

a = ALPHA

17
4.3. Dimensionalization

Although the physical input in DISTRIB is required to be in units of

micrometers and hours, it is computationally more efficient to solve

the partial differential equations with the physical quantities ex- pressed in other units. For the kinds of problems of interest in impurity redistribution, a convenient basic length is BLTH = 0,02 ym.

Given the basic length unit, a basic time is frequently introduced in

diffusion theory by defining

2 BTM = BLTH_

j^Q^j^g (4.10)

D

o

where Dq is some diffusion coefficient in the problem. In our program,

Dq is chosen to be D2, the diffusion coefficient in silicon. When these new units are used to dimensionalize the diffusion equation in the silicon, one finds the new equation has a diffusion constant of unity. That

is , making a change of independent variables

/ _ t / _ y , / ^ z ^ - BTM' ^ ~ BLTH' ^ BLTH

transforms eq (2.14) into

8 , C2 = C 2 - ° / ^2

t y

and eq (2.15) into

\' ^1 = 'I' ^1 •

The rest of the dimensionalizations carried out in block #4 of section

7.1 are standard conversions from micrometers and hours to the new

units BLTH and BTM.

4.4. The Description of the Way Arrays CI and C2 Are Plotted

The plotting is really a peripheral item and not described in as much

detail as substantive parts of the program. All the quantities un- defined are defined or referenced in the glossary and their definition

in the glossary makes their function clear. First, the way CI is

plotted is described. Let NC be the integer part of (NL - 1)/NM0D.

When NC > 0, we plot the concentrations CI and z-coordinate (in units

of ZOPl) at every grid point K = V • NMOD + 1 for V = 1,... , NC. When

NC = 0, we plot the concentrations at z = 0 and at the moving boundary

Zq. The quantities XlS, YlS determine the largest abscissa and ordinate on the graph of CI versus z.

18
The graph of array C2 is plotted in the following way. All abscissa

are in units of UN2 (which is essentially SCALE2, since y is still in

units of BLTH) . The first C2 concentration plotted is CB(2) at the

moving boundary. The next concentration is C2 (NR) with abscissa S(2)DY

Next, we plot ND - 2 concentrations with the indices KS = NR + (V - 2)

NMODl for V = 3,..., ND. Next, we pass these grid points through a filter which eliminates those with an abscissa greater than X2S/2 = X2SCT. The quantities X2S and Y2S determine the largest abscissa and

ordinate that appear on the C2 versus y plot. Note: The specific

scales used as the ordinates and abscissa for these plots were chosen

in order to make the resultant graphs identical in size to those obtained from the equipment which experimentally determined the concen- tration data. This also accounts for the elimination of points beyond

X2S/2. Both conditions may be changed to accommodate the user.

The smallest abscissa and ordinate are determined in both graphs by

(0,0)

.

The plot subroutine is an NBS "house program" that is documented in Appendix 2.

This plotting procedure has worked well in all cases, the only adjustments needing to be made being in NMOD and NMODl (see table 20 in

sec. 9.22).

19
5. DERIVATION OF THE DISCRETE ALGEBRAIC EQUATIONS

5.1. The General Prescription for Obtaining the Tridiagonal Equation

Associated with a Fixed Grid Point

This section, which should be read casually, gives a rough description

of the precise derivations in the following sections.

Each particular tridiagonal equation in eqs (3.2), (3.3) [and (3.4)

also] associated with some particular grid point is derived by discretizing some integral form of the conservation of number eqs (2.20) to

(2.22) over the "flux cell" (to be described below) that contains that

particular grid point. The flux cells are intervals that arise by partitioning the z- and y-regions in a certain way. That is, both the z- region and the y-region are broken up into nonover lapping intervals

called flux cells or concentration cells in such a way that the intervals completely cover the regions [0,Zo(t)] and [yo(t),BDRY] at every

time t; and each grid point is contained in at most one cell except for

the grid points with indices NL and NR. Each of these grid points is

contained in both the irregular sized cell with which it is associated

and the "moving boundary cell"; these statements will be explained more completely later. The two boundaries of each of the concentration cells

are called the flux boundaries of the cell. We interpret the concentration, C1(J,TTMMP1) = C1((J - l)Az,ti4M+i), at some grid point J to represent the average number density in the cell at time t^^+j.

The flux boundaries for the cells in the y-region are shown in figure 7 by the vertical lines located halfway between successive grid points.

Notice that all the cells in the y-region that contain fixed grid points

(see sec. 3.1) are of width Ay except the irregular one of width 0.5 Ay

that contains and leads to the tridiagonal algebraic equation associated

with the grid point with index NR at (NR - l)Ay. This cell has one flux

boundary at (NR - l)Ay and the other is at (NR - 1/2) Ay. Similar re- marks apply to figure 8.

The significance of the flux cell containing some grid point, with index I in the y-region (see fig. 7) is that it leads to a linear algebraic equation containing the concentration C2(I,TTMMP1) as an unknown

when the conservation eq (2.21) defined on this cell is discretized.

NR NRPI I-l lime = t^M + i

time = tMM

Itl N2 •—Trrr 1 ' 1

1

1 N2

f 1 1 1 1 1 1

IRREGULAR FLUX CELL

t

REGULAR FLUX CELL

Figure 7: Flux cells and flux boundaries in the y-region for fixed

grid points.

20
time = tMM+i

Figure 8. Flux cells for fixed grid points in the z-region when NL >

In particular, let yj^ and yg in eq (2.21) be at (I - 3/2) Ay and

(I - 1/2) Ay, respectively. Make now the following approximations in

eq (2.21) (with t = t^^. ) : MM+1

d

dtf

C2(y,t)dy ^ C2(y,t )]dy/At MM

^2(y'Wl^^^ '^^ C2((I

/

C2(y't^)dy C2((I - l)Ay,t^)Ay ,

Do8 Co

y

D2[C2((I - l)Ay,t^^^ ) - C2((I - 2)Ay,t^^^ )] , and

y=yA

C2 Z D2[C2(lAy, t^^^) - C2((I - DAy, t^^^^ ) ] . MM+1

y=yB

Then eq (2.21) becomes, up to terms of higher order, the linear tridiagonal equation in eq (3.4) associated with the grid point with index I

C2((I - 2)Ay,t^^^) - |2 -H C2((I - l)Ay,t^^i)

+ Co(lAy,t ) = - C2((I - l)Ay,t ) MM+1 n^A-i- D2At ^ MM

21
The important fact this example illustrates is that, by evaluating

the fluxes in eq (2.21) at the flux boundaries of the flux cell and

evaluating the integrals by a quadrature approximation with a "node"

located at the grid point in the flux cell, one ends up with a tridiagonal (i.e., an equation involving the concentrations at three con- secutive grid points) equation for the unknown concentration at time

''^MM+1 that grid point. Thus, with the help of these flux cells one can in a natural way associate with each unknown at a fixed grid point

an algebraic equation.

The flux cells of the points in the z-region depend on the values of NL.

There are three cases depending on whether NL = 1, NL = 2, or NL > 2.

The cell configuration for the case NL > 2 can be visualized from figure

8. In the case NL = 2 (see fig. 9) , there are just two irregular cells

each of width 0.5 Az. The flux cell associated with the grid point with

index NL has its right side flux boundary located coincident with the

grid point. The flux cell associated with the first z-grid point with

index 1 has its left-hand flux boundary coincident with the first grid

point. These facts also hold when NL > 1, but in this case there are regular grid cells, i.e., with width Az in between the irregular ones

at the outer boundaries. The case NL = 1 is special because there are

no grid points in the z-region between the moving and fixed boundaries.

Finally, we point out that for all values of NL the conservation formula

eq (2.22) (conservation across the moving boundary) is discretized to

derive the preliminary eq (3. 4), and in this derivation the flux cell is as shown in figure 10 with the flux boundaries at z^ = (NL - l)Az and

yg = (NR - l)Ay. The trapezoidal rule is used to evaluate the integrals

in eq (2.22) with the nodes at NL, NR and the moving boiindary grid

points at S(l)Az + (NL - l)Az and (NR - l)Ay - S(2)Ay, see figures 5 and 6.

These different cases, which arise because we have flux cells with different widths and with different orientation of grid points relative to

the flux boundaries for different values of NL, explain the slight variations, classified in tables 3 and 15, in the way the tridiagonal algebraic equations and their "P,Q" coefficients (see sec. 3.1) are derived

in the following sections.

NL=2

1^ \ f

time = +1

time = t MM

Figure 9. Flux cells for fixed grid points in the z-region when NL = 2.

22
Figure 10. Schematic description of the flux cell associated with the

discretization of the conservation equation across the moving

boundary

.

Table 3: Guide to the Sections Where the Tridiagonal Algebraic Equations

Associated With the Grid Points are Described and Derived.

Case I: The moving boundary has not reached the second z-grid point or NL < 2.

(a) The index of the z-grid point is 1 (see sec. 5.9)

(b) The index of the y-grid point is greater than NR (see sec.

5.3)

(c) The index of the y-grid point is NR (see sec. 5.5)

(d) The grid point is at the moving boundary (see sec. 5.7)

Case II: The moving boundary has reached the second z-grid point or NL ^ 2.

(a) The index of the z-grid point is 1 (see sec. 5.8)

(b) The index of the z-grid point ranges from 2 to NL - 1 (see sec. 5.2)

(c) The index of the z-grid point is NL (see sec. 5.4)

(d) The grid points are at the moving boundary (see sec. 5.6)

(e) The index of the y-grid point is greater than NR (see sec.

5.3)

(f) The index of the y-grid point is NR (see sec. 5.5)

5.2. The Index of the Z-Grid Point Ranges from 2 to NL - 1 When NL > 2.

In this case we suppose the index, J + 1, is such that 2<J+1<NL-1,

and begin with the formula, see eq (2.20)

23
t)dz = - Di 3 Ci

z

^ (5.1)

For the z-grid point, J + 1, we assume zp^ and Zg "the flux boundaries"

to be at (J - 1/2) Az, (J + 1/2) Az, respectively. We speak of these grid

points as regular grid points because the distance betv/een the cell flux

boundaries is Az.

The left-hand side of eq (5.1) is discretized by employing the backward

time difference formula (A7) in Appendix 1 between the times t

= and t = tiyu4 + At. Here, At is some small positive time increment that

sometimes depends on the index MM. After the derivative has been re- placed by the backward difference quotient, the integrals are approximated by the midpoint quadrature rule (A2)

Finally, after further algebraic manipulation we have

• z^ Ci(z,t)dz = II [Ci(JAz,t^ + At)

z.

- Ci(JAz,tj^)] + 0(Az2 + At2)/At . (5.2)

Since all the expressions in eq (5.1) are assumed to be evaluated at

t = t]yiM + At, the approximation of the derivatives on the right-hand

side of eq (5.1) gives

9 Ci

z

and

9 Ci

z

^

= [Ci (JAz,t + At) - Ci((J - l)Az,t,^^ + At)]/Az + 0(Az) (5.3)

J- MM MM

z=ZA

= [Ci((J + l)Az,t,,, + At) - Ci(JAz,t^„ + At)]/Az + 0(Az) .(5.4) MM MM

Z=Zt

By substituting eqs (5.2) to (5.4) into eq (5.1) and further rearranging,

we obtain up to higher order terms the discrete version of eq (5.1) in

tridiagonal form,

DiCi((J - l)Az,t^ + At) - (2 • D + Az2/At) Ci(JAz,tj^ + At)

+ DiCi ( (J + l)Az,t^ + At) = ^-^^^ Ci (JAz,t^)

24
After introducing the FORTRAN quantities in table 4, eq (5.5) becomes

DFl • C1(J,TTMMP1) - (2 • DFl + DZ**2/DT) • CI ( J + 1, TTMMPl)

(5.6)

- (DZ) **2

+ DFl • C1(J + 2), TTMMPl) • CI ( J + 1,TTMM)

Table 4: FORTRAN Version of Quantities in Eq (5.5).

Di = DFl

Cl( , ) = Cl( , ) = TTMM

MM

t + At = TTMMPl

MM

At = DT

Az = DZ

Ci(JAz,t.„) = CI (J + 1,TTMM) MM

Ci (JAz, t + At) = C1(J + 1, TTMMPl) MM

In the FORTFIAN program, the following intermediate quantities in table 5 have been introduced for convenience and computational efficiency.

Table 5: Intermediate Quantities Used in Eq (5.7).

DFl = A2L

-(2 • DFL + (RAl)"^) = A2M

DFl ^ A2R

- (RA1)~^ = AIM

AF(J + 1) = AIM • C1(J + 1,TTMM)

In terms of these intermediate quantities, eq (5.6) can be written in

the compact tridiagonal FORTRAN form

A2L • CI (J, TTMMPl) + A2M • Cl ( J + 1, TTMMPl) + A2R • CI ( J + 2, TTMMPl)

(5.7) = AF (J + 1) .

As has been explained in more detail in section 3.1, the quantities

A2L,..., AF(J) are used in computing the values of the arrays PI, Ql,

and Dl. Specifically, the quantities A2L, . . . , A2R, AIM are computed in

25
subroutine INTER. The subroutine INTER is called by subroutine TRIDNL

when this latter subroutine is calculating the arrays PI, Ql. The array

AF(J) is computed in block #25 of the main program, see section 7.1.

SUMMARY

Schematic summary of the derivation of the tridiagonal FORTRAN equation

for z-grid points with indices between 1 and NL:

(5.1) -> (5.5) (5.6) (5.7) .

Summary notes:

(1) Quantities in table 5 except for AF are computed in subroutine

INTER.

(2) The array AF (J) is computed in block #25 of the main program.

5.3. The Index of the Y-Grid Point is Greater than NR

The derivation of the tridiagonal FORTRAN form of the algebraic equations associated with these grid points is almost entirely analogous

to the derivation of the equations in siibsection 5.2 once the substitutions z for y, for D^, for C^, etc., are made. However, there is one exception; in the y- region we have divided the equivalent of eq

(5.5) by D2 so that the equivalent of that equation is

C2((J - l)Ay,t _ + At) MM

2 + Ay2_]

D2AtJ D2At

+ C2((J + l)Ay,t

C2(JAy,t„„ + At)

MM

MM

+ At) = 5f^C2(JAy,t^)

(5.8)

Table 6 is the equivalent of table 4.

Table 6: FORTRAN Version of Quantities in Eq (5.8).

D2 = DF2

Ay = DY

C2(-,-) = C2(-,-)

Az = DZ

C2(JAy,t,,) = C2(J + 1,TTMM) MM

C2(JAy,t + At) = C2(J + 1,TTMMP1)

26
When the quantities in table 6 are introduced into eq (5.8), it becomes

[DY**2 "I 2 + ^ff~r~^ \ C2(J + 1,TTMMP1) + C2(J + 2,TTMMP1)

(5.9)

DY**2 = C2 (J + 1,TTMM) . DF2 • DT

This form of the equation is further simplified by introducing the intermediate quantities in table 7.

Table 7: Intermediate Quantities Used in Eq (5.10).

DT ^ = RA2

DY2

1 = B2L

-(2 + RA2~VdF2) = B2M

1 = B2R

-RA2"VdF2 = BIM

BF(J + 1) = BIM • C2 ( J + 1,TTMMP1)

In terms of the intermediate quantities in table 7, eq (5.9) assumes

the form

B2L C2(J, TTJ4MP1) + B2M C2 (J + 1, TTMMPl) + B2R C2 (J + 2, TTMMPl)

(5.10)

= BF(J + 1) .

The quantities RA2 , B2L,...,B1M are computed in subroutine INTER. This

subroutine is called by subroutine TRIDNL and in this routine the quantities RA2,...,B1M are used to compute arrays P2 , D2. The quantity

BE (J) is computed in block #25 of the main program.

SUMMARY

Schematic siammary of the derivation of the tridiagonal FORTRAN equations

for y-grid points with indices greater than NR:

(5.8) —> (5.9) —> (5.10) .

Sammary notes:

(1) Quantities in table 7 except for BE are calculated in subroutine INTER.

(2) The array BE (J) is compuced in block #26 of the main program.

27
5.4. The Index of the Z-Grid Point is NL and NL > 1

In this subsection we derive the FORTRAN tridiagonal form of the equation associated with the z-grid point with index NL which is located at

z = (NL - l)Az.

The derivation begins with the integral conservation equation

d

dt/

-— / Ci(z,t)dz -

A

-Di9 Ci

^ z

^ Z=Z7

Di9 Ci (5.11)

Z=Zt

In this equation, Zg = (NL - l)Az, Zji^ = (NL - 3/2) Az are the "flux

boundaries"; notice the flux boundary at Zg coincides with the location

of the grid point associated with the cell (see fig. 9). Also, the

cell width, the distance between the flux boundaries, is Az/2 , and

therefore this grid point is called an irregular grid point, see sec- tion 5.1.

By using the backward difference quotient approximation in eq (A7) to

approximate the time derivative on the left-hand side of eq (5.11) between the times tjyuyi and ty^j^ + At and the right end point quadrature approximation in eq (A3) to approximate the resulting integrals, the

left-hand side of eq (5.11) becomes equivalent to

d_

dt /

(NL-1) Az

Ci(z,t)dz = ff^ ]ci[(NL - DAz, + At]

(NL-3/2) Az (5.12)

- Ci[(NL - l)Az, t^] + 0(At + AzVAt)| .

The spatial derivative on the right-hand side of eq (5.11) at z = z^ is

approximated by the central difference approximation

'd Ci z=z

|^|ci[(NL -DAz, t^ + At] - Ci[(NL - 2)Az, t^ + At]|

A (5.13) + 0(Az) .

To approximate the derivative at z = zg on the right-hand side of eq

(5.11), we use the three point derivative approximation given in eq

(A9) . Setting in that formula

28
XI = z (t + At) o MM

X = (NL - l)Az = z X2 = (NL - 2)Az

h - Az

B

we obtain the sought-after result

1 ( 9 Ci z

^ S (1)

Al )- 1 + S(l) ^1^^^^ - 2)^-' + At] z=zB

+ [(S(l) - 1)/S{1)] Ci[NL - DAz, t^^^ + At] MM (5.14)

where

+ —7T-^

rr^ TTTT Ci [z (t + At) , t + At] } + 0(Az2)

S(l) [1 + S(l)] ^ o MM MM I

S(l) = [z (t^^^ + At) - (NL - l)Az]/Az

o MM (5.15)

When eqs (5.12), (5.13), and (5.14) have been substituted into eq (5.11)

and the expression simplified, it becomes up to terms of higher order

DiS(l)

-

1 + S(l) ^1^^^^ - 2)^^' + At]

r^i • s(i)

R(l) Ci [ (NL - l)Az, t + At]

MM

(5.16)

(1 + S(l)) ^if^o^^MM ^MM

+

S(l) Ci [(NL - DAz, t^^] MM

R(l)

where

R(l) = 2At

Az2

(5.16')

After introducing the FORTRAN notation in table 8 into eq (5.15), it

becomes

29
DFl DFl DFl C1(NL - 1, TTMMPl) + ,\ + ,\ CI (NL,TTMMP1)

1 + S(l) R(l) S(l)

(5.17)

DFl

CB(1) - CKNL, TTMM)

S{1) (1 + S(l)) R(l)

Table 8: FORTRAN Version of Quantities in Eq (5.16).

Di = DFl

t = TTMM

MM

t. . + At = TTMMPl

MM

At = DT

Az = DZ

Z • DT

R(l) - DZ**2

Ci((NL - 2)Az, t + At) = CKNL - 1, TTMMPl)

^ MM

Ci((NL - l)Az, t + At) - CKNL, TTMMPl) L MM

Ci((NL - l)Az, t^) = CKNL, TTMM)

Cl(z L (t + At), t + At) = CB(1)

o MM MM

Equation (5.17) will be used to obtain the "P, Q" quantities PA(1),

QA(1) for z-grid point NL, i.e., quantities such that

CKNL, TTMMPl) = PA(1) • CB(1) + QA(1) , (5.18)

in section 6.6.

SUMMARY

Schematic summary of the derivation of the tridiagonal FORTRAN equation

for the z-grid point with index NL when NL >_ 2

,

(5.11) —^ (5.16) —> (5.17) .

Summary note

Equation (5.17) is used in section 6.6 to obtain quantities PA(1)

,

QA(1) .

5.5. The Index of the Y-Grid Point is NR

The derivation of the tridiagonal FORTRAN formula for this grid point

is analogous to the derivation in section 5.4 for the irregular z-grid

30
point at NL so the user should refer to that section if more details

concerning the following derivation are desired.

Again, we start with the integral conservation equation.

•4B

C2(y,t)dy = -D23 C2

y

+ Do9 Co

y

y=yA

(5.19)

y^yB

where y = {NR-l)Ay and y = (NR - i)Ay.

By using the same approximation procedure that led from eq (5.11) to eq

(5.12) (except for using the left end point quadrature approximation),

we find

dt / ^ C2(y^t)dy = 1^ jC2[(NR - l)Ay, t^ . At]

- C2 [ (NR - 1) Ay, t ] ^ MM

(5.20)

The spatial derivative on the right-hand side of eq (5.19) at y = yB is

approximated by the central difference approximation (see eq (A8)

,

with h = Ay/2)

3 Co

Y

^ Ay

C2(NRAy, t^^ + At) - C2 [ (NR - l)Ay, t^^ + At])-+ 0(Ay) MM

y=yB ;5.2i)

The spatial derivative at y = y^ on the right-hand side of eq (5.19) is

approximated by using the three point formula eq (A9) with:

We find

8 <Z.

XI = NRAy

X = (NR - l)Ay

X2 = y (t^„ + At) . o MM

T- It-^-^^W C2((NR)Ay, t + At) AyJl + S(2) ^ ^ MM

y=yA

. -(2)'^^ C2((NR- DAy, t^ . At) (5.22)

S(2)(l ^ S(2)) ^^(^o^^MM ^ ^MM ^ ° ^'^^ '

31
where

S(2) = -[yo(t^ + At) - (NR - l)Ay]/Ay . (5.23)

When eqs (5.20), (5.21), and (5.22) have been substituted into eq (5.19)

and the resulting expression is simplified, it becomes

-

1 /s(2) ^2[y^(t^^ At), t^ -f At]

(S—(2) \ + ij C2[(NR - DAy, t^ + At] (5.24)

2(2) S (2) -

1 + S(2) C2(NRAy, t^ + At) = ^

C2[.(NR - l)Ay, t^ + At] ,

where

2 D2At

After introducing the FORTRAN notation in table 9 into eq (5.24), it

becomes

1 + S(2) ^^^^^ TTMMPl) +(^^Jjj-+

Ij C2(NR, TTMMPl)

(5.26)

S ( 2 ) ^ ( 7)

-, , '(^s C2(NR + 1, TTMMPl) = C2 (NR, TTMM) .

Table 9: FORTRAN Version of Quantities in Eq (5.24).

D2 = DF2

t = TTMM

MM

t.„, + At = TTMMPl

MM

At = DT

Ay = DY

2 • DF2 • 03? R(2) = DY**

C2(Y (t„^ + At), t + At) = CB(2)

o MM MM

C9((NR - l)Ay, t + At) = C2 (NR, TTMMPl)

^ MM

C2((NR - l)Ay, t^^) = C2 (NR, TTMM)

32
Equation (5.26) will be utilized later in section 6.9 to find the quantities PA(2), QA(2) such that

C2(NR, TTMMPl) = PA(2) CB(2) + QA(2) (5.25')

Schematic summary of the derivation of the tridiagonal FORTRAN equation

for the y-grid point NR,

(5.19) —> (5.24) —> (5.26)

Summary note

Equation (5.26) is used in section 6.9 to determine the quantities

PA(2), QA(2) associated with grid point NR.

5.6. The Grid Points Lie at the Moving Boundary and NL > 1 (The Preliminary Equation)

The main equation associated with the two unknowns Ci [z^iti^ + At) , t^M

+ At], C2[yQ(tf4j4 + At), tf42^^ + At] located at the grid points which form

the interior boundary of the z- and y-regions, respectively, is derived

from the integral form of the conservation of number boundary condition

eq (2. 22) , namely

.z.(t,,,,+At) ^y, ^ ' — - C2(y,t)dy

/o MM r^B

Ci(z,t)dz + /

^A "^^o^W^t)

(5.27)

- Di3 Ci + Co

y

where z = (NL - l)Az and y = (NR - l)Ay.

Furthermore, it is assumed that the moving boundary z^{t) is such that

(NL - l)Az ^ Zq (t) <_ NLAz, for all times t between and tj^jyi + At.

Similarly, for times in this same interval the moving boundary yQ(t) is

assumed to satisfy

(NR - 2)Ay <_ y^(t) <_ (NR - l)Ay .

With these restrictions, eq (5.27) is equivalent to the differential

boundary condition eq (2.16).

Before proceeding with the discretization of eq (5.27), it is convenient

to introduce immediately the FORTRAN notation in table 10.

33
Table 10: FORTRAN Symbols Used in Eq (5.28).

Ci((NL - l)Az, t... + At) = C1(NL, TTMMPl)

Ci({NL - l)Az, t ) = CKNL, TTMM) J- MM

C2((NR - l)Ay, t..,. + At) = C2(NR, TTMMPl) MM

C2((NR - l)Ay, t„) = C2 (NR, TTMM)

^ MM

Ci (z (t + At) , t + At) = CB(1, TTMMPl)

^ o MM MM

Co(y (t + At) , t + At) = CB(2, TTMMPl)

o MM MM

^l^^o^^MM^' ^MM^ = TM^

By using the backward difference quotient approximation in eq (A7) to

approximate the time derivatives on the left-hand side of eq (5.27) and

using the trapezoidal rule to approximate the integrals in the resulting expression, we find (using the notation introduced in table 10) that

the left-hand side of eq (5.27) is equivalent to

S (1) Az S (2) Ay C1(NL, TTMMPl) + CB(1, TTMMPl) + CB(2, TTMMPl)

SM (1) Az

+ C2(NR, TTMMPl) - Cl(NL, TTMM)

(5.28) + CB(1, TTMM) - ^^2At^'^ CB(2, TTMM) + C2 (NR, TTMM)

+ higher order terms in Ay, Az, At ,

where

S(l) - iz (t^, + At) - (NL - l)Az]/Az

o MM

SM(1) = [z (t ) - (NL - l)Az1/Az

o MM

S(2) - [(NR - DAy - Y^it^ + At)]/y (5.29)

SM(2) = [(NR - DAy - y (t,,J]/Ay . o MM

These are consistent with previous definitions of S(l) and S(2),

When eqs (5.14) and (5.22) (re-expressed in the FORTRAN notation of table 10) are employed to approximate the spatial derivatives on the

right-hand side of eq (5.27) and these approximations along with eq

(5.28) are substituted into eq (5.27), there results after simplification and up to terms of higher order

34
^-^^C1(NL, TTMMPl) + CB(1, TTMMPl) + ^——^-^CB(2, TTMMPl)

R(l) R(2)

+ C2(NR, TUVIMPl) - ^^,1]^ CKNL, TTMM) + CB(1, TTMM

)

R ( 1

)

- ^ ' CB(2, TTMM) + C2 (NR, TTMM)

R ( 2

)

^ -Dl

I

-

I Vs\^ ) ~ 1' TTMMPl)

+ (S(l) - 1) C1(NL, TTMMPl) (5.30)

where

S(1)(A s(i)) TM""'|

* - S(2)(l t S(2))

'^^''-

^^^^ C2(NR, TTMMPl)

S { 2

)

S (2) ) C2 (NR + 1, TTMMPDJ. , 1 + S(2)

RM ^ = m

RM = ^

(5.31)

Ay

1 DR = — . D2

In order to eliminate the division by S(2) (which could have the value

zero) in the last term on the right-hand side of eq (5.30), we eliminate

that term by substituting from the relation

1 1 - S (2) - S(2)(l + S(2)) ^^^2' ^ S(2) ^2^^^'

S (2) + -. + g(2) C2(NR, TTMMPl)

(5.32)

= C2 (NR + 1, TTMMPl) - C2 (NR, TTMMPl)

1 R(2) [-C2(NR, TTMMPl) + C2 (NR, TTMM)]

35
This latter equation is just the formula [equivalent before simplification to eq (5.24)] that results when eqs (5.20), (5.21), and (5.22) are substituted into eq (5.19) and the FORTRAN notation of table 9 introduced .

If the terms in the brackets on the right-hand side of eq (5.30) are eliminated by substitution from eq (5.32) and if the resulting equation

is multiplied through by S(l) (to prevent possible division by zero),

then we obtain [after CB(2, TTMMPl) has been replaced by (m = SEG) • CB(1, TTMMPl) using boundary condition eq (2.23)]:

S(l)'-

R(l)

Dl SEG S(l) • S(2)

1 + S(l)

S(l)2

R(2)

+ (1 - S(l))

CB(1, TTMMPl)

R(l)

R • S (1) • S (2)

R(2)

C1(NL, TTMMPl)

- R • S(l) C2(NR, TTMMPl)

S(l)^ • Dl

1 + S(l)

S(l) • SM(1)

R(l)

R • SM(2) • SM(1)

C1(NL - 1, TTMMPl) + R • S (1) C2 (NR + 1, TTMMPl)

R(2)

C1(NL, TTMM) + CB(1, TTMM)

CB(2, TTMM)

(5.33)

+ S(l)

R • SM(2) R \ ^ , + —-T— |C2(NR, TTMM)

R(2) R(2) / Further FORTRAN notation in table 11 can be introduced into eq (5.33)

and by using intermediate variables, it can be more compactly rewritten:

DM1 • CB(1, TTMMPl) = DM2 • C1(NL, TTMMPl) + DM3 • C2 (NR, TTMMPl)

(5.34) + DM4

where

C1(NL - 1, TTMMPl) + DM5 • C2 (NR + 1, TTMMPl) + DM6 ,

DM1 =

DM2

DM3 =^ -

S(l) DFl SEG RT S(l) • S(2)

R(l) 1 + S(l)

S(l)**2

R(2)

R(l)

RT • S(l) • S(2)

DM4 = DFl

R(2)

S (1) **2

+ DFl • (1 - S(l)) ,

- RT • S(l) 1 +

(2) )

1 + S(l)

(5.34')

36
DM5 = RT • S (1) , and

S(l) • SM(1)

DM6 = —r- C1(NL, TTMM) + CB (1 , TTMM) R ( 1

;

RT • SM(2) • S(l) + TT^^ — CB(2, TTMM)

+ S(l) .

• -

. C2(NR, TTMM) , R ( 2

)

and where

DF2 • DZ

R = RT =

R(l) =

R(2) =

DY

2 • DT

. DZ**2 '

2 • DF2 • DT

DY**2

S(l) = (ZOPl - (NL - 1) • DZ)/DZ , (5.35)

S(2) = ( (NR - 1)DY - Y0P1)/DY ,

SM(2) = ( (NR - DDY - y„(t ))/DY , and ^ MM

SM(1) = (z (t ) - (NL - 1)DZ)/DZ . o MM

Table 11: FORTRAN Version of Symbols in Eq (5.33).

At = DT

Ay = DY

Az = DZ

z (t +

o MM

At) = ZOPl

y (t +

-^o MM

At) - YOPl

Di = DFl

D2 = DF2

SUMMARY

Schematic summary of the derivation of the FORTRAN equation for the z- grid point at the moving boundary, i.e., at (^^(tj^ + At), when NL > 1;

(5.27) —> (5.33) —^ (5.34) and (5.34')

37
Summary notes

(1) All the concentrations on the right-hand side of eq (5.34),

i.e., C1(NL, TTMMPl) , . . . ,C2 (NR + 1, TTMMPl) will be eliminated

in terms of known quantities and thereby CB(1, TTMMPl) will be

found. This is carried through in section 6.10.

(2) The concentrations in DM6 are known because they are evaluated

at time TTMM.

5.7. The Grid Point Is at the Moving Boundary and NL = 1

The derivation of the preliminary equations for C 1 ( Zq ( tj^j + At), tf^ +

At) is very similar to the derivation of the preliminary equation for

the same quantity in the previous section, so if the user desires more

details than are given below, he should refer to that section.

We again begin with the conservation eq (5.27). However, instead of

approximating the spatial derivative at z = by eq (5.14), we use boundary condition eq (2.24). Because of this change, we find in place

of eq (5.30)

S (1) —TT-T- [C1(NL, TTMMPl) + CB(1, TTMI4P1) ] R ( 1

)

R • S(2)

R(2)

SM(1)

R(l)

R • SM(2)

[CB(2, TTMMPl) + C2(NR, TTMMPl)]

[C1(NL, TTMM) + CB(1, NR, TTMM)]

[CB(2, TTMM) + C2(NR, TTMM)] (5.36) R(2)

= - [+ C Az(Cl(l, TTMMPl) - C )]

p out

+ R

1 1 - S (2) CB(2, TTMi4Pl) + C2(NR, TTMMPl)

S(2) (1 + S(2) ) ^ ' '

' S(2)

S(2)

1 + S(2)

C2 (NR + 1, TTMMPl)

When the term in the second parenthesis on the right-hand side of eq

(5.35) is eliminated by using eq (5.32), there results

S (1)• [C1(NL, TTMMPl) + CB(1, TTMMPl)]

R ( 1

)

+ ^ • ^^^^ [CB(2, TTMMPl) + C2(NR, TTMMPl)]

R(2)

38
^^777 [C1(NL, TTMM) + CB(1, TTMM)

]

R ( i

)

R • SM(2)

R(2) [CB(2, TTMM) + C2(NR, TTMM)] (5.37)

= - [C Az(Cl(l, TTMMPl) - C J] p out

+ R [C2(NR + 1, TTMMPl) - C2 (NR, TTMMPl)]

R

R(2) [-C2(NR, TTMMPl) + C2(NR, TTMM)]

By rearranging terms and introducing the FORTRAN symbols in eq (5.35)

into eq (5.37), we get

DM1 • CB(1, TTMMPl) = DM2 • Cl (NL, TTMMPl) + DM3 • C2 (NR, TTMMPl)

+ DM4 • C1(NL - 1, TTMMPl) (5.38)

+ DM5 • C2(NR + 1, TTMMPl) + DM6 ,

where

S(l) RT '3(2) • SEG

DMl = + , R(l) R(2)

DM2 = - ^-TTT" + CONPRO • DZ , R(l)

= - iif - ^(^ ' ^) ' ^'-'^'^

DM4 = 0 ,

DM5 = RT , and

DM6 = ^^11^ [C1(NL, TTMM) + CB(1, TTMM)]

R ( X

)

• SM(21fc3(2, TTMM)]

R(2)

RT • SM(2) - RT

[C2 (NR, TTMM)

]

R(2)

+ DZ • CONPRO • CONOUT ,

and where the FORTRAN symbols above are the same as in eq (5.35) with

the addition of the ones in table 12.

39
Table 12: Additional FORTRAN Symbols for Eq (5.38).

C = CONPRO

P

C ^ - CONOUT

out

Notice that the quantities DM3, DM5, and DM6 [with DZ • CONPRO • CONOUT

added to it in eq (5.34')] are equal to their counterparts in eq (5.38')

when S (1) - 1. Hence, by setting S(l) to 1 when we are computing these

quantities, i.e., in the case NL = 1, we can use the same FORTRAN expres

sion as is used to compute these quantities in eq (5.34') when NL > 1.

Of course in computing in this way we must leave the S(l) in DMl and

DM2 of eq (5.38') unchanged. This can be done, when NL = 1, by introducing the quantity SAVE which is equal to S(l) and replacing S(l) where

it appears in DMl and DM2 by the symbol SAVE.

SUMMARY

Schematic summary for obtaining the preliminary equation for the z-grid

point at the moving boundary when NL 1.

(5.27) —> (5.37) --^ (5.38) and (5.38') .

Summary notes

(1) Since the coefficients of the concentrations on the righthand side of eq (5.38) are the same as those on the righthand side of eq (5-34), the concentrations from both of these

preliminary equations can be eliminated by the same procedure

which will be described in section G.IO.

(2) The actual computation of the quantities in eq (5.38') are carried through as described above in computing block #4 of

subroutine CENTER, see section 7.2.

5.8. The Index of the Z-Grid Point is 1 and NL > 1

The derivation of the tridiagonal algebraic equation that is associated

with the z-grid point 1 (when the moving boundary has reached the sec- ond z-grid point) begins with the conservation of number eq (2.20)

(5.39)

where z^^ = 0 and Zg = Az/2 denote the positions of the flux boundaries.

d

dt/

Ci (z, t)dz = -Di3 Ci

z

^ + Di9 Ci z=z.A

40
After approximating the time derivative on the left-hand side of eq

(5.38) by the back difference quotient approximation eq (A7) and the

integrals by quadrature formulas with a node at the left-hand end point,

we find

d

dt

2 (t +At)

o MM

Ci(z,t)dz = —— [Ci(0 • Az, t + At) - Ci (0 • Az, t )]

^ 2At ^ MM ^ MM

Az

(5.40) + 0(Az^ + At^/At) .

By substituting from eq (2.24),

C (C1(0 • Az, t.„, + At) - C ) p MM out

for Di 9 Ci

z on the right-hand side of eq (5.39) and then

z=z.

D, (Ci (Az, t + At) - Ci (0 • Az, t + At)

)

1 1 MM ^ MM

for DiS Ci z on the right-hand side of eq (5.39), that equation

z=z

B

yields together with eq (5.40) after simplification

Az^

+ Az • C + Di

P

Ci(0 • Az, t,,, + At) = DiCi(Az, t^,^ + At)

^ MM Mil

+ : + Az • C • C ^ 2 • At p out

Az2

(5.41)

After introducing into eq (5.41) the FORTRAN notation in table 13, that

equation becomes

/ 2dt \

^ \DZ2j + DZ • CONPRO + DFl

Did, TTMM)

CKy, TTMMPl) = DFl • C1(2, TTMMPl)

(5.42)

2DT

DZ**2

+ DZ • CONPRO • CONOUT

41
Table 13: FORTRAN Version of Quantities in Eq (5.41).

At = DT

Az = DZ

C = CONPRO

P

C = CONOUT

out

Di = DFl

Ci(0 • Az, t „ + At) = Cl(l, TTMMPl) MM

Ci(0 J- • Az, t,.) - Cl(l, TTMM) MM

Co(Az, t + At) = C2{2, TTMMPl)

^ MM

Equation (5.42) is the FORTRAN tridiagonal equation which will give rise

to the "P, Q" quantities associated with the z-grid point 1 and unknown

Cl(l, TTMMPl). This is carried through in section 6.4.

SUMMARY

Schematic summary of the derivation of the FORTRAN tridiagonal equation

for the z-grid point 1 when the moving boundary is beyond z-grid point 2

(5.39) —> (5.41) —> (5.42).

Summary note

The FORTRAN eq (5.42) leads to values of Pl(l), Ql(l) that are com- puted in computational block #24 of the main program. See section

5.4 for the derivation of these "P, Q" quantities.

5.9. The Index of the Z-Grid Point is 1 and NL = 1

This equation is different from all the others in that it is not developed from an integral conservation equation. It is basically derived

from combining eq (2.14) with a Taylor series expansion.

From a spatial Taylor's formula expansion of Ci (z, tMM + At), we have

Ci (z (t + At) ) - Ci (0, t + At) + 9 Ci

^ o MM ^ MM z

^ z (t + At) o MM

z=0 (5.43)

+ 1/2 aci (z (t^^ + At))2 + 0(Az3) . o MM

z=0

42
By defining

S(l) = z (t + At)/Az

o MM

and substituting this together with (from eq (2.14))

into eq (5,43), it becomes after rearranging and up to higher order terms

DnCi i^oMM (z (t + At) , t + At) = DnCi (0, t + At) + Di 9 Ci (S(l)Az) MM ^^ MM -^z^

z=0

+ 1/2 3^Ci (S(l) • Az)^

(5.44)

z=0

When the backward difference approximation,

3^Ci = (Ci(0, t^^ + At) - Ci(0, t,^))/At + 0(At)

t . MM MM

and the expression from eq (2.24),

Dia^ci = C (Ci (0, t^^ + At) - C , ) , p MM out

(5.45)

z=0

are substituted into eq (5.44), it becomes, after rearrangement (ignoring

higher order terms)

,

At • 2 p

Ci(0, t + At) MM

= DiCi(z (t^ + At), t^^ + At) o MM MM

(5.46)

+ C, (0, t ) ^ MM

(S(l) Az)

2At

+ S (1) • Az • C • C

p out

By transferring the symbols in eq (5.46) through the correspondences in

table 14, that equation becomes

(S(l) • DZ)**2 DFl + ^ 2 ' DT ^^"^^ *

* ' TTMMPl)

= DFl CB(1, TTMMPl) + Cl(l, TTMM)

+ S(l) • DZ • CONPRO • CONOUT .

(S(l) • DZ)**2

2DT (5.47)

43
Table 14: FORTRAN Version of Quantities in Eq (5.46).

Di = DFl

Az = DZ

At = DT

Cl(0, + At) = Cl(l, TTMMPl) MM

Cl(0, tj = Cl(l, TTMM)

MM

Ci (z (t + At) , t + At) = CB(1, TTMMPl)

^ o MM MM

C = CONPRO

P

C ^ = CONOUT

out

FORTRAN tridiagonal eq (5.47) will lead to the quantities Pl(l), QKD

in section 6.3. These are the "P, Q" associated with z-grid point 1 when

the moving boundary has not yet reached the second z-grid point.

Notice that when t^j^^ = 0 the quantity Cl(l,0) is unknown, i.e., it is not

given as data in the problem. Thus, the equation we use in place of eq

(5.45) for the first iteration MM = 1 when tj,^ = 0 is

(DFl + S(l) • DZ • CONPRO) • Cl(l, TTMMPl) = DFl • CB(1, TTMMPl)

(5.48)

+ S(l) • DZ • CONPRO • CONOUT .

This is the equation we would have arrived at if we had truncated the

Taylor's formula in eq (5.43) to first order terms and proceeded just as above, and eq (5.48) is consistent with, i.e., the same as, the formula

we would get by approximating eq (5.45) directly.

SUMMARY

Schematic summary of the derivation of the FORTRAN tridiagonal equation

associated with the z-grid point 1 when NL = 1,

(5.43) —> (5.45) —^ (5.47) or (5.48) (use (5.48) when t^ = 0).

Summary note

The "P, Q" terms that result from eq (5.47) or eq (5.48) are derived in sections 6.2 and 6.3.

44
6. METHOD OF SOLVING THE ALGEBRAIC EQUATIONS

6.1. The Derivation of the "Double Sweep" "P, Q" Coefficients Associated

with the Grid Points

In this section we put each of the tridiagonal FORTRAN equations associated with the fixed grid points in "double sweep" form. Specifically, we'

manipulate these tridiagonal FORTRAN eqs (3.2) and (3.3) into the forms

C1(J, TTMMPl) - P1(J) CKJ + 1, TTMMPl) + Ql ( J) (6.1)

for J = 1,...,NL - 1, and

CKNL, TTMMPl) = PA(1) CB(1, TTMMPl) + QA(1) . (6.2)

Also, at the fixed grid points in the y-region,

C2(J, TTMMPl) = P2(J) CI (J + 1, TTMMPl) + Q2 ( J) (6.3)

for J = NR + 1,...,N2, and

C2(NR, TTMMPl) = PA(2) CB(2, TTMMPl) + QA(2) . (6.4)

Once again the specific composition of the double sweep "P, Q" factors

associated with a grid point depends on the location of the grid point

and the position of the moving boundary at time TTMMPl. The following

table 15 completely classifies all the cases that arise. Of course, be- cause the "P, Q" coefficients for a grid point are obtained from the tridiagonal equation associated with that grid point, this table is the same as table 3 except for Case I (e) and Case II (g) which were not included in table 3 (since they are so trivial)

.

Table 15: Guide to the Subsections Where the "P, Q" Coefficients Associated With Each of the Fixed Grid Points Are Derived and

Described.

Case I : The index NL equals 1

:

(a) The first z-grid point, see sections 6.2 and 6.3.

(b) The index of the y-grid point is greater than NR, see sec- tion 6.8.

(c) The index of the y-grid point is NR, see section 5.9.

(d) The grid point lies at the moving boundary, see section

6. 10.

(e) The largest y-grid point, which has index N2 and is located at y = BDRY = NYO • DY, see section 6.7.

45
Case II: The index NL is greater than 1.

(a) The index of the z-grid point is 1, see section 6.4.

(b) The index of the z-grid points ranges from 2 to NL - 1, see section 6.8.

(c) The index of the z-grid point is NL, see section 6.6.

(d) The grid point is at the moving boundary, see section 6.10.

(e) The index of the y-grid point is greater than NR, see sec- tion 6.8.

(f) The index of the y-grid point is NR, see section 6.9.

(g) The index of the y-grid point is N2 and is at y = BDRY, see section 6.7.

6.2. The First Z-Grid Point When NL = 1 and t = TTMM = 0

MM

After dividing both sides of eq (5.48) by (DFl + S(l) • DZ • CONRPO)

,

we get

Cl(l, TTMMPl) = Pl(l) • CB(1, TTMMPl) + Ql(l) (6.5)

where

Pl(l) = DM • DFl

Ql(l) = DM • (Cl(l) • DM3 + DM1 • CONPRO • CONOUT)

DM1 = S(l) • DZ

DM2 = 0

DM3 = DM2

DM = (DFl + DM2 + DMl • CONPRO)"^ .

The introduction of the zero quantities DM2, DM3 into the above equations enables us to use below in section 6.3 the same expression DM but

with DM2, DM3 ^ 0.

SUMMARY

Equations (6.5') are calculated in computational block #23 of the main

program (see sec. 7.1).

6.3. The First Z-Grid Point When NL = 1 and t = TTMM > 0

MM

f S (1) • DZ) **2

When eq (5.47) is divided through by DFl + ^ ^ 2 DT

^ ^^^^ * ' CONPRO, there results

Cl(l, TTMMPl) = Pl(l) CB(1, TTMMPl) + Ql (1)

with

(6.5')

(6.6)

46
(6.6')

Pl(l) - DM • DFl

QUI) = DM • (Cl(l, TTMM) • DM3 + DMl • CONPRO • CONOUT)

DM = (DFl + DM2 + DMl • CONPRO)

DMl = S(l) • DZ

DM2 = (DMl) 2/2 • DT

DM3 = DM2 .

SUMMARY

Note that all quantities in eq (6.6') are the same as in eq (6,5') except

for DM2 and DM3. The calculation of the quantities in eq (6.6') takes

place in computational block #23 of the main program (see sec. 7.1).

6.4. The Index of the Z-Grid Point is 1 and NL > 1

(DZ) **2

After dividing eq (5.42) by + DZ • CONPRO + DFl), it becomes ^ ^ 2 • DT

Cl(l, TTMMPl) = Pl(l) • Cl(2, TTMMPl) + Ql (1) (6.7)

where

Pl(l) = DFl • DM

DZ**2

01(1) = DM Cl(l, TTMM) • + DZ • CONPRO • CONOUT (6.7')

2 • DT

DM = ^

£_ + DZ • CONPRO + DFl . 2 • DT

SUMMARY

The quantities in eq (5.7') are computed in computational block #24 of

the main program (see sec. 7.1).

6.5, The Index of the Z-Grid Point is Between 1 and NL - 1 and NL > 2

We begin with eq (5.7) which holds for the indices J - 2,...,NL - 1. Assume that for one of the indices in this range

C1(J - 1, TTMMPl) = P1(J - 1) • CKJ, TTMMPl) + Ql ( J - 1) . (6.8)

After substituting this value of Cl ( J - 1, TTMMPl) into the left-hand

I side of eq (5.7) and rearranging the terms, there results

CKJ, TTMMPl) = PI (J) • C1(J + 1, TTMMPl) + Ql (J) (6.9)

where

47
PI

Ql

Dl

(J)

(J)

(J)

= -A2R/D1(J)

= (AF(J) - QKJ - 1)A2L)/D1(J)

= A2M = PI (J - 1) • A2L .

(6.9')

Notice that eq (6.8) with J - 2 holds because of eq (6.7). Hence, by induction we get eqs (6.9) and (6.9') for J = 2,...,NL - 1.

Summary notes

(1) The arrays Pi (J) , Dl (J) are computed in computational block

#3 of subroutine TRIDNL.

(2) The quantities A2L,...,A1M (see table 5) are computed in subroutine INTER which is called by subroutine TRIDNL when Pi (J) , Dl(J) are being computed.

(3) The arrays AF(J), Ql (J) are calculated in computational block

#25 of the main program.

6.6 The Index of the Z-Grid Point is NL and NL > 1

When eq (5.17) is multiplied through by S(l) (since S(l) may be zero,

this prevents possible division by zero) and the terms are transposed,

there results

SUMMARY

DM1 C1(NL, TTMMPl) = DM2 CB(1, TTMMPl)

+ DM3 C1(NL - 1, TTMMPl) + DM4 (6.10)

with

DM = 1 + S(l)

DM = DM2 DF1/(1 + S(l)) = DFl/DM (6.10'

)

DM4

DM3 DFl • S (1)/(1 + S (1) ) = DFl

S (1)

;., : Cl (NL, TTMM) .

DM2

We also have from eq (6.9) with J = NL - 1,

C1(NL - 1, TTMMPl) = P(l) C1(NL, TTMMPl) + Q(l) (6.11)

where we have used the abbreviations

P(l) = P1(NL - 1)

Q(l) = QKNL - 1) 6.11')

48
When eq (6.11) is substituted into the right-hand side of eq (6.10) and

the result is simplified, we find

C1(NL, TTMMPl) = PA(1) CB(1, TTMMPl) + QA(1) (6.12)

where

PA(1) = DM2/ (DM1 - DM3 • P(l))

(6.12')

QA(1) = (DM4 + Q(l) • DM3) /(DM1 - DM3 • P(l)) .

SUMMARY

The quantities in eqs (5.10') and (6.12') are evaluated in computational

I block #3 of subroutine CENTER.

:i 6.7. The Index of the Y-Grid Point is N2 = NYO + 1 and is located at

y = BDRY

I The boundary condition for C2 ( , ) at the boundary point y = BDRY is

j (see eq (2. 12)

)

C2(N2, TTMMPl) = CBULK . (6.13)

j We rewrite this simple equation for use later on in the more convenient

l' form

C2(N2, TTMMPl) ^ P2 (N2) C2 (N2 - 1, TTMMPl) = Q2(N2) , (6.14)

if we assume (and we do)

P2(N2) = 0

Q2(N2) = CBULK .

SUMMARY

|| The quantities in eq (6.14^) are calculated in computational block #13

of the main program. These quantities are used in block #26 of the main

:.rogram (see sec. 7.1) and block #4 of siabroutine INTER (see sec. 7.3).

6.8. The Index of the Y-Grid Point is Between NR and N2

We start with tridiagonal eq (5.10) which holds for J = NR + 1,...,

112 - 1.

Assume for one of the indices in this range that

C2(J + 1, TTMMPl) = P2(J + 1) C2(J, TTMMPl) + Q2(J + 1) . (6.15)

After substituting this value of C2 (J + 1, TTMMPl) into the left-hand

side of eq (5.10), with J in place of J - 1, and rearranging terms, it

becomes for J = NR + 1, . .

. , N2 - 1 49
C2{J, TTMMPl) = P2(J) C2(J - 1, TTMMPl) + Q2(J) (6.16)

where

1

P2(J) = - D2 (J)

Q2(J) = (BF{J) - Q2(J + 1) B2R)/D2(J) (6.16')

D2(J) = B2M + P2(J + 1)B2R

or, since B2R = 1 from table 7,

D2 (J) = B2M + P2 (J + 1) .

Notice that eq (6.15) with J = N2 - 1 holds because of eq (6.14). Hence,

by induction we get eq (6.16) and eq (6.16') for J = NR + 1,...,N2 - 1.

SUMMARY

Summary notes

(1) The arrays P2 (J) , D2(J) are calculated in computational block

#4 of subroutine TRIDNL.

(2) The quantities B2M, BlM (the other B's equal 1) are computed in

subroutine INTER which is called by subroutine TRIDNL when

P2 (J) , D2 (J) are computed.

(3) The arrays Q2(J), BF(J) are calculated in computational block

#26 of the main program.

6.9. The Y-Grid Point has Index NR

After transposing terms in the tridiagonal FORTRAN eq (5.26), it takes

the form

DM1 C2(NR, TTMMPl) = DM2 CB(2, TTMMPl)

(6.17) + DM3 C2(NR + 1, TTMMPl) + DM4

where

DMl = 1 +

DM2

R(2)

1

^ + '^'2> (6.17')

DM3 = S (2)/(l + S (2) ) S (2) DM4 = ^l^j

C2 (NR, TTMM) .

We also have from eq (6.15) with J = NR + 1

50
C2 (NR + 1, TTMMPl) = P(2) • C2 (NR, TTMMPl) + Q(2) (6.18)

where we have used the abbreviations

P(2)=P2WR+1)

Q(2) = Q2(NR + 1) .

When eq (5.18) is substituted into the right-hand side of eq (6.17) and

the result is simplified, we find

C2(NR, TTMMPl) = PA(2) CB(2, TTMMPl) + QA(2) (6.19)

where

PA (2) = DM2/ (DM1 - DM3 • P(2))

QA(2) = (DM4 + Q(2) • DM3) /(DMl - DM3 • P(2))

SUMMARY

(6.19')

The quantities eqs (6.1?') and (6.19') are evaluated in computational

block #3 of subroutine CENTER.

6.10. The Determination of CB(1, TTMMPl) in Terms of Known Quantities

The values of the "P, Q" coefficients in the eqs (6.1) through (6.4) are determined by the coefficients in the tridiagonal eqs (3.2) and (3.3)

and these in turn are composed of known mesh parameters, material parameters and concentrations evaluated at time TTMM. Hence, by using eqs (6.1)

through (6.4), we can eliminate the concentrations evaluated at time

TTMMPl that occur on the right-hand side of the "preliminary" eq (5.34)

and thereby determine the value of CB(1, TTMMPl) entirely in known quantities .

More precisely, substituting eqs (6.11) and (6.18) into the right-hand

side of the "preliminary equation," i.e., eq (5.34) for NL > 1 (and eq

(5.38) for NL = 1), we get

(6.20)

DMl CB(1, TTMMPl) = DM8 • Cl (NL, TTMMPl)

+ DM9 C2(NR, TTMMPl) + DMlO

where

DM8 = DM2 + DM4 • P(l)

DM9 = DM3 + DM5 • P(2) (6. 20'

)

DMlO = DM6 + Q(l) • DM4 + Q(2) • DM5 .

Finally eliminating the concentrations at the grid point NL and NR on the

right-hand side of eq (6.20) by use of eqs (6.12) and (6.19) leads to the

51
expression for CB(1, TTMMPl) in terms of known quantities

DMll CB{1, TTMMPl) = , (6.21)

with

DM12 = DM1 - DM8 • PA(1) - SEG • DM9 • PA (2)

DMll = DM10 + QA(1) • DM8 + QA{2) • DM9 , and where

DM1

DM2

DM3

DM4

DM5

DM6

. eq (5.39) if NL = 1 are given m /c o/in • ,.tt i eq (5. 34) if NL > 1

and furthermore.

PA(2)

QA(2) are given in eq (6.19')

^^2) given in eq (6.18 )

PA(1) . . ,„/. are given m eq (6.12 ) QA(1)

P(l)

Q(l) are given in eq (6.1l')

(6.21')

(6.22)

R(l) is given in eq (5.16 ) R(2) is given in eq (5.35)

RT = R is given in eq (5.31)

SUMMARY

The quantities in eqs (6.21), (6.2l') and all the quantities listed below

eq (5.2l') are computed in computational block #4 of subroutine CENTER to

find the value of CB(1, TTMMPl).

52
7. DETAILED DESCRIPTION OF THE PROGRAM

7,1. Block-by-Block Detailed Description of the Calculational Procedure

in the Main Program

As is shown in the program listing (Chapter 11) , the overall calculations

in DISTRIB: VERSION 1 are broken down into a set of consecutive blocks.

The computations in each of these blocks have a unified purpose which is

described with pertinent comments for each of these blocks #1 through

#36 in this chapter. In Chapter 8, the logic in these blocks is flow

charted. The precise calculations carried through in each of these

blocks are shown in the program listing in Chapter 11, and the material

in the following blocks is intended to illuminate the contents by organizing and elaborating on the material in the listing rather than to re- describe the contents in a different language.

IMPORTANT : Quantities with names which appear in the program listing v but do not appear in the glossary are not needed or used in DISTRIB:

VERSION 1 and therefore should be ignored, i.e., all material in the

listing that is not in a block or glossary.

Block #1: Input Declaration, Type Statements, Dimensionalize Arrays

and Common Blocks

The arrays SRL{ ) and YRL( ) are declared real because of certain con- straints in the plotting routines where these arrays are employed. The

fact that they are dimensionalized by 200 means that only 200 data points

can be plotted in a single call of the plot routine.

Arrays that have an element associated with each grid point, e.g., C2 ( J)

(see glossary) are dimensionalized by 2000. This means the maximum number of grid points is 2000 and the maximum number of mesh widths in either region must be less than 2000.

Since the number of cells of the array Cl ( ) actually used is, at time

TTMMPl, equal to NR • NZO, the input quantities NZO and TFINAL (which

together with the equation of the moving boundary determine the largest

possible value NR takes) must be constrained so that the product NR • NZO is always less than the dimension of Cl ( ).

Many of the arrays are dimensionalized in common blocks because they are also used in subroutines.

Block #2: Input the Proper Problem and Material Parameters with Accompanying Format Declarations

The definition of each of the input quantities is to be found in the

glossary, where references are given to places in the main body of this

documentation which contain more elaborate descriptions of these quantities. The list of input and output is given in section 9.2 in a systematic fashion and in classified form.

53
Real input is read under format 4 and real output printed under format

5. Integer input and output are under format 1.

A table listing the input data is found in section 9.2 along with ranges

of values for the input parameters that have worked well in trial runs. The units for the input data are micrometers, hours. The problem with

CONOUT = 0 is homogeneous in the concentrations, and therefore the input concentrations are normalized to CBULK =1. By homogeneous is meant

the value of C. (x,t) for CBULK = A is equal to C. (x,t) [for CBULK = 12

1 1 txmes A.

Block #3: Store Input in Terms of Their Original Dimension Units

In the next block certain input parameters are dimensionalized (rescaled)

Consequently, these input quantities are stored in their original units

in the array PARAM so they can be printed out later in block #34 in their

original units. Also, before its defining quantities are changed, the

quantity SCALE2 (in micrometers which is used to scale the distaTice along

the y-coordinate system, i.e., in the silicon, in the plots that are printed in block #34 is calculated and stored for later output. The rea- son SCALE2 does not appear in the plot routine is that its role there is

played by UN2 which represents the same physical quantity as SCALE2 but

is in units of BLTH rather than micrometers. The fact that DY, YOPl,

etc., when they appear in the plot routines are in BLTH units requires

the use of UN2 instead of SCALE2. A detailed description of the dimensionalization procedure is given in section 4.3.

Block #4: Redimension the Input Data by Introducing Computationally

Efficient Units

The dimensionalization carried out here is described in section 4.3 in

detail. Note that, in Version 1, DFC is set equal to 1 in a declare

card. Thus, during the calculations, DF2 will have the value 1 and DFl

will be the ratio DF1/DF2. The computational formula for BTM is given

in eq (4. 10)

Block #5: Compute the Values of the Grid Mesh Widths

This is the step described in section 4.1 where DY and DZ, the mesh

widths in the y- and z-regions, respectively, are calculated. The com- putational formula for DY is eq (4.1) and eq (4.3) for DZ.

Block #6: Compute the Value of ZCRIT = CNSTX.

In this step CNSTX, which is another name for ZCRIT, is calculated, see section 4.2. The quantities CNA and CNB are just two intermediate variables used in computing DNSTX and ZOPl, Z0P2 later on.

54
Block #7: Compute the Time Increment DT for the First Iteration of the

Main Loop 78

In this block the main purpose is to compute the initial value of At

(= DT) . Recall, from section 4.2, that in the "wet oxidation" Case I the time TLIN = CNTAU = x = 0, so the motion of the moving boundary is

always parabolic. In this case DT is defined as the time it takes the

moving boundary to move the distance DM = WFRC • DZ, having started at

Z = 0. The number WFRC must be positive, real, and less than 1, for

the program to be stable. In the case of a "dry oxidation," Case II in

section 4.2 which starts at Z = 0 (in VERSION 1) TLIN ='CNTAU = t > 0,

the boundary moves linearly for the time period CNTAU with a speed

SLOPE = ZCRIT/CNTAU. Hence, the time DT the moving boundary needs to

transverse the distance DM (where DM is as above) is DM/SLOPE. Values

of WFRC that have been successfully tested range between 0.01 and 0.33.

We have required TLIN to be less than or equal to 10~^ rather than zero

in order to distinguish between the wet and dry oxidation cases. This

does not seem now to be necessary, but we have left this as it is be- cause in all practical cases TLIN is always much greater than 10~^ whenever it is greater than zero, so no harm is done.

The quantity DTS is used to store the quantity DT for recall at a later

step, see computational block #17.

When CNTAU = t = 0, the quantities ZCRIT = CNSTX and SLOPE which are not

used but still appear in later formulas are assigned the values 0 and 1,

respectively, which do not foul up these formulas, for example by causing division by zero.

In case x > 0 (i.e. , dry oxidation) , WFRC must certainly be chosen so that DM < ZCRIT for the above prescription for DT to make sense. In

all the cases we have run, this has been naturally satisfied when data

in the literature are satisfied.

Block #8: Computation and Storage of Initial Value DT in Hours

The quantity DT just calculated in the previous block is multiplied by

BTM (the unit in which time has been rescaled) to determine its magnitude in units of hours. It is stored in the array PARAM for printout

in block #34.

Block #9: Calculation of Initial Values of NL, NR and Value of N2 for

the First Iteration of the Main Loop 78

The quantity Nl is another name for the index, NL, of the first grid

point with a z-coordinate less than or equal to the z-coordinate of the

moving boundary at the beginning time of the main iteration loop 78.

Because YOl , ZOl are zero in Version 1, Nl has the value 1, as it should

since at this point in the program the moving boundary is at z = 0 (recall we are computing here NL for the first iteration of loop 78) NR

55
(which is the index of the. first grid point that has a y-coordinate

greater than the y-coordinate of the moving boundary at the start of

the iterations in loop 78) has the value 2,

The quantity N2 is the index of the last grid point on the y-axis and

since NYO is the number of grid mesh widths in the interval (0, BDRY)

(see below eq (2.11) for the definition of BDRY) N2 = NYO + 1 (see com- putational formula eq (4.2)).

Block #10: Calculation of Initial Values of C2 ( )

The initial values of array C2 ( ) at all the grid points are determined

from the prescription in eq (2.13) of C2(y,0). Because of this prescription, the grid points from NR through NR + NH (NH is an input integer) are assigned the value CMAX (input value) while all the following grid

points are assigned the value CBULK (input value)

Since in VERSION 1 no oxide is initially present, there is no necessity

to assign array Cl ( )

Notice that C2 (1) is essentially assigned through the assignment of

CB(2). This is because CB(2) denotes by definition the concentration

in the silicon at the position of the moving boundary and initially the

moving boundary is located at the position of the first y-grid point.

Block #11: Calculation of the Value of TTMM for MM = 1

This is the step in which TTMM, the time at the beginning of the MMth

iteration of loop 78, for the first iteration is calculated. Since ZOl

in VERSION 1 has been assigned the value zero, in a- declaration card,

this quantity obviously has the correct value which is zero.

Block #12: Calculation of the Initial Values of SM( )

The quantities SM(1), SM(2) are, respectively, the distance between grid

point NL and ZOPl (the position of the moving boundary in the z- coordinate system at time TTMMPl) and the grid point NR (for more details see figs. 5 and 6)

The quantity FNL represents the niamber of mesh widths between z-grid

point 1 and z-grid point NL, and FNR is the number of mesh widths between y-grid point NR and y-grid point N2 (the last grid point on the

y-axis in the domain of the silicon)

Since ZOl and YOl are declared zero in VERSION 1, SM(1) = 0 and SM(2) =

1 as they should be when we are calculating these quantities for the

initial iteration of loop 78.

56
Block #13: Calculation of the "P, Q" Coefficients for Grid Point N2

As has been explained in section 6.7, the satisfaction of the boundary

condition C2(N2) = CBULK (see eq (2.12)) at y = BDRY is equivalent (in

the way this problem is being solved) to assuming P2(N2), Q2(N2) are assigned in such a way that C2(N2) = P2(N2) C2 (N2 - 1) + Q2(N2) [see

eqs (6.14) and (6.14') ]

.

The values of P2(N2), Q2(N2) are used in subroutine TRIDNL.

Block #14: Beginning of the Main Loop which Calculates Concentrations

at Time = TTMMPl

At the beginning of this loop, 78, the quantities in the arrays CI ( ),

C2 ( ) representing the concentrations in the z- and y-regions at each

grid point have their values at time TTMM. At the end of the MMth itera

tion, these same quantities will have their proper values at time TTMMPl

Before the main loop begins, certain starting values for quantities used

in the loop are set. The definitions and use of these quantities ZOPS,

TJ, etc., are explained in later steps.

The overall manner in which the new concentrations are calculated is

explained in section 3.

Block #15: Setting the Print-out Control Parameter

The integer KW controls print-out in block #30. This data will be

printed out after the first iteration, then skip KPRINT - 1 iterations

of loop and print out the data again, and so on.

Block #16: Increasing the Distance that the Moving Boundary Travels

in an Iteration of Loop 78

During the MMth iteration of loop 78, the quantity ZOACEL measures that

fraction of the mesh width DZ that the moving boundary transversed during the previous iteration of loop 78. (If MM - 1, ZOACEL = 0 since

ZOPl = 0. and ZOPS = 0. In this case, see block #14.) If ZOACEL is

less than 0.2 and KPOW (defined below) is greater than 100, then the

magnitude to DT is increased by 0.05 in order to accelerate the distance

the boundary moves during the MMth iteration of loop 78. This acceleration is put into the algorithm in order to decrease the number of iterations that are required in order to compute the solution to a problem in

which the moving boundary moves a certain specified distance.

Of course, subroutine TRIDNL, which computes and stores quantities that

depend on DT, has to be recomputed and the new value of DT is placed in

DTS. When leaving the loop in which DT has been increased, the value of

the control parameter KPOW is increased by 100 so that DT will not be

increased again for at least another 100 iterations of loop 78.

57
The quantities 0.2 and 100 mentioned above are, of course, somewhat

arbitrary, but the above values have worked satisfactorily in a great

number of trial runs.

Block #17: Adjustments Required When Moving Boundary Reached New Z-Grid

Point on Previous Iteration of Loop 78.

The quantity TJ will be computed in a later step, block #20, in such a way that when and only when the moving boundary reached a new z-grid

point during the previous step will TJ have a value greater than or equal to 1. Hence, assuming TJ >_ 1 , NL, which is the first grid point

less than or equal to the position of the moving boundary at the time

when MMth iteration of loop 78 starts, has to be increased in value by

one. The quantity FNL (the number of mesh widths to the left of NL)

also must be increased by 1 and NLMl = NL - 1, NLM2 = NL - 2, increased

accordingly.

Also, CI at the new grid point, NL, must be defined.

Furthermore, since DT has been adjusted (see block #20) during the previous iteration of loop 78, this quantity is reset to the value it had

during the (MM - 2)th iteration, specifically DTS. Also, TRIDNL, the

subroutine computing the arrays PI, Ql, etc., depends on DT (besides a new PI, Dl being added in the y-region) and hence needs to be recomputed.

The value SM(1) (see sees. 5.6 and 5.7 for its meaning) must be set

equal to zero because during the (MM - l)th iteration the moving boundary reached a new grid point in the z-region.

All these changes are effected when TJ > 1. Of course, in the contrary

case the implication is that the moving boundary had not yet reached a new z-grid point and the above changes are not required.

Block #18: Computing the Time and Position of the Moving Boundaries at

the End of the MMth Iteration of Loop 78

In this step we are simply computing what the value of the time will be

at the end of the MMth iteration of loop 78. Since at the beginning of

the loop, t has the value TTMM by definition and since time is made to

increase by DT during the loop, the time at the end of the loop TTMMPl

is given by TTMM + DT. After this value of TTMMPl is computed, the location of the moving boundary in the z- and y-regions, respectively,

ZOPl, YOPl can be computed through the use of formulas in section 4.2

and eq (2.8). Remember, ZW and TW equal zero in VERSION 1.

Block #19- Adjustments Required When Moving Boundary Reached New YGrid Point on Previous Iteration of Loop 78

The quantity SES is the fraction which the signed distance between NR

and YOPl is of the mesh width DY = AY. If YOPl is beyond the y-grid

point NR (during the MMth iteration of loop 78) , then SES is negative

58
and NR, FNR, and NRPl = NR + 1 have to be increased by 1. When SES is

negative at the beginning of the MMth loop, it means (see the next block)

the moving boundary is located at time TTMM on a grid point, and hence

SM(2) should be set equal to 1.

Block #20: If the Moving Boundary Moves for a Time 2DT, Will It Reach

a New Z-Grid Point?

During the MMth iteration of loop 78, we look at where ( ) is at time

TTMMP2 = TTMMPl + DT. This position is denoted as Z0P2. Then we com- pute TJ which is that fraction the signed distance between grid point

NL and ZOPl is of the mesh width AZ = DZ . If this fraction is less

than one, it means that during the (MM + l)th iteration of loop 78 the

moving boundary will not reach a new z-grid point, and we proceed with

the MMth iteration. On the other hand, if TJ is greater than or equal

to 1, we know that the moving boundary will reach or cross grid point

NL + 1 on the next iteration (MI'I + 1) . Therefore, we decide to change

the value of DT so that while still in the MMth iteration of loop 78,

Zq travels all the way to the next grid point. At the time TTr4MPl when

this happens, i.e., Zo(TTr4MPl) = (FNL +1), DZ is calculated (as shown

in the listing) for the two cases CNTAU = 0 or CNTAU > 0. (Remember

that ZW = TW = 0 in VERSION 1.) Finally, we calculate the new value of

DT that will take the moving boundary Zq ( ), starting at time TTMM, to

its new position. .Lastly, we recompute TRIDNL because the arrays Pi,

P2, etc., computed there depend on the value of DT.

Block #21: Calculation of the Parameters S(l), S(2) Needed for Subroutine CENTER

The quantities FNL, FNR, etc., have now been reset (if it was necessary).

These quantities S.(l) and S(2) (see sees. 5.4 to 5.7 for details) are used in computing the coefficients of the unknowns in the algebraic

equations that will be solved later in subroutine CENTER.

Block #22: For the First Iteration of Loop 78, Call TRIDNL

In the first iteration of loop 78, the arrays P2, D2 that are needed

later and which are calculated in subroutine TRIDNL have not been calculated before in the program.

Notice that whenever a new NL occurs, and hence the quantities in TRIDNL

need to be updated, this is done in step #20.

Block #23: Calculating of Pl(l), QKD When the Moving Boundary Zq

Has not Yet Reached Grid Point Two

In the case NL = 1, the moving boundary at time TTMMPl has still not

reached the, second z-grid point. Therefore, Pl(l), Ql(l) must be calculated using formulas in sections 6.2 and 5.3, specifically eqs (6.5')

when TTim. = 0 and (6.6') when TTMM > 0.

Of course, DM, DMl, etc., are just intermediate variables.

59
For the computational significance of Pl(l), Ql (1) , see again sections

5.2 and 6.3. Pl(l) is not calculated in subroutine TRIDNL like the rest

of its elements, and Ql(l) is calculated here rather than with the other

elements of Ql in block #25.

Block #24: Computation of Pl(l), Ql(l) When Zq Has Moved Beyond the

First Grid Mesh Widths

In this case (NL > 1) , the moving boundary has traveled through the

first z-mesh width and Pl(l), Ql(l) are to be calculated by the formulas

in section 6.4, specifically eq (6.7'). Pl(l) is not calculated in subroutine TRIDNL like the rest of its elements, and Ql(l) is calculated

here rather than with the other elements of Ql in block #25.

Block #25: Calculation of Arrays AF ( ), Ql ( ) VJhen There are Regular

Points in the Z~Region

When NL > 2, there are grid points with cell widths (see sec. 5.1) of

magnitude DZ in the z-region, and the arrays AF and Ql associated with

these points are therefore computed. [This amounts to saying that the

arrays AF ( ), Ql ( ) that are associated with these points are calculated

by using the formulas in table 5.2 and eq (6.9'). ]

E*or the significance of the array Ql in the computational procedure,

see chapter 3. The array AF plays the role of the quantity Zi+ in eq

(3.2); see also eq (5.7). Array AF is used to compute array Dl in subroutine TRIDNL. The array Ql, which is composed of AF and Dl, is used

in block #29.

Block #26: Calculation of Arrays BF and Q2 Associated with the Grid

Points in the Y-Region

The array BF is used to compute Q2 later in this block.

For further elucidation of the computational significance of these arrays, see the references to BF in section 5.3, table 7, and eq (5.10).

See section 6.8 and eqs (6.16) and (6.16') for the significance and

definition of Q2(J). Notice P2 (N2) , Q2(N2) are calculated in block

#13. Arrays P2, Q2 are used in block #28.

Block #27: Calculation of CB(1), CB(2), CI (NL) , C2 (NR) at Time TTMMPl

The quantities above at the end of the MMth iteration of loop 78 (i.e.,

at time TTMMPl) are computed in subroutine CENTER which is documented

in detailed block form in section 7.2.

Block #28: Computation of C2 ( ) at Time TTMMPl

By making use of C2 (NR) that was computed in subroutine CENTER in block

#27, we can carry out the backward sweep; see steps #2 and #3 in chapter

3 that generate C2(J) for J = NR + 1 to N2 - 1.

60
Block #29: Computation of Array CI ( ) at Time TTMMPl

In this case NL = 1, the value of Cl(l), which is the only grid point

besides CB(1) in the y-region, is already known, since it has been com- puted in subroutine CENTER as CI (NL) . However, when NL is greater than

one, the calculation of CI (J) for J = 1, NL - 1 must be completed by a backward sweep, as described in section 3.1, steps #2 and #3.

Block #30: Print-out Concentrations When KW = 0

The integer KW is zero when (MM - 1) is divisible by KPRINT,

The concentrations are printed under a time TTMMPl which is in units of

BTM (see sec. 4.3). The user is advised not to print out every concen- tration in each region because there are so many, especially in the yregion (i.e., the silicon).

Also printed out are the values of the concentrations at the boundary

and CI (NL) . (Note that C2 (NR) will always be printed with the other con- centrations in the y-region.

)

Block #31: Reset Quantities in Preparation for Next Iteration of Loop

78

Unless a new grid point is reached by the moving boundary in the next

iteration, the new values of SM(1), SM(2), TTMM are as given (see figs.

5 and 6) . If in the next iteration of loop 78 a new grid point is

reached, the values are corrected appropriately at that time, in blocks

#17 and #19.

Block #32: Test to Determine Whether TTMMPl Has Reached TFINAL

The quantity TFINAL (INPUT) is the time we want to stop the problem

and TIMEND is the value of TFINAL in units of BTM which is the unit of

TTMMPl. Hence, as long as TIMEND is greater than TTMMPl, we proceed to

the next iteration of loop 78. If MME is not large enough to reach

TFINAL, a diagnostic to this effect is printed before control shifts to

statement 79 and then the end card.

Block #33: Shift Program to Statement 800 in VERSION 1

The time is greater than or equal to TFINAL and in VERSION 1 the value

of JSTOP is set equal to 1 in a data card. Hence, control sends the

program to statement 800 where printing and plotting occur prior to the

end of the program.

Block #34: Print-out of Data and Plot of Concentrations at the Final

Value of TTMMPl

The detailed description of the plotting subroutine is given in Appendix

2 and section 4.4.

61
Block #35: Redimensionalize Parameters to Original Units

Parameters such as ZOPl and TTMMPl (at the last iteration of loop 78

that occurs) are changed into micrometer and hour units and printed out

with other pertinent parameters.

Block #36: Control in VERSION 1 Sends Program to Its End

The integers JSTOP and JSAVE have the values one and zero (see below

data cards) in VERSION 1; consequently, in this case the program is sent

to the end.

If TIMEND has not been reached by TTI-lMPl before the number of iterations MME of loop 78 is completed, a diagnostic message to this effect

is printed out.

Block #37: This Block Contains Statement 79 and Ends Program.

7.2. Block-by-Block Description of Subroutine CENTER

This subroutine, which is called in computational block #27 of the main

program, is where the quantities CB(1), CI (NL) , CB(2), C2 (NR) are evaluated at time TTMMPl.

Block #0. Insert Input Through COMMON

The following input for the subroutine passes in through COMMON: S( ),

SM( ), N2M1, N2, NLMl, NL, NR, CI ( ), C2 ( ), CB ( ), DFl, DF2 , CONPRO,

CONOUT, DT, DZ, DY, PI, P2, Ql, Q2

The output is fed to the main program through COMMON.

Block #1: Computation of the Quantities P(2), Q(2)

These quantities defined in eqs (5.18) and (6.18') are used to determine PA(2), QA(2) in computational block #3.

Block #2: Compute P(l), Q(l)

These quantities [see eqs (6.11) and (6.11^)] are used when NL > 1 to

determine PA(1), QA(1). When NL 1, the role played by PA(1), QA(1)

in linking the value of CB(1) to the value of Cl(l) is played by Pl(l),

Ql(l), which are computed in block #23 of the main program. Therefore,

at the end of block #3, we set PA(1) = Pl(l), QA(1) = QKD when NL - 1.

Block #3: Computation of PA(L), QA(L) for L = 1,2

The computation formula for these quantities when L = 1 is described

in eqs (6.10'), (6.1l'), and (6.12'). [See eqs (5.16') and (5.25) for

the definition of R(l) and R(2); also see tables 7 and 8.] The computation formulas of the above quantities in the case L = 2 are described

62
in eqs (6.17'), (6.18'), and (6.19'). The significance of PA(L), QA(L)

is seen from eqs (5.26') and (5.18) or eqs (6.12) and (6.19).

Block #4: Computation of CB(1), CB(2) at Time TTMMPl in Terms of Knovm

Quantities

The concentration at the moving boundary in the z-region CB(1, TTMMPl)

CB(1) is computed by using formula eqs (6.21) and (6.2l') and the aux- iliary eq (6.22). To find CB(2) = CB(2, TTMMPl), we just use the boundary condition in eq C2.23). For the significance and explanation of

how SAVE is used, see the end of section 5,7, before SUMMARY.

Block #5: Calculation of CI (NL) , C2 (NR) at Time TTMMPl

We introduce the abbreviations

C(l) = C1(NL, TTMMPl) = CI (NL)

C(2) = C2(NR, TTMMPl) = C2 (NR)

and compute these desired quantities from CB(1), CB(2), PA(L), QA(L)

known from blocks #1 to #4 by employing eqs (6.12) and (6.19).

7.3. Block-by-Block Description of Subroutine TRIDNL

The input for this subroutine is the output of subroutine INTER, DT,

DZ, DY, DFl, DF2, and CONPRO. The arrays Pi, Dl, P2, D2 in contrast

to the Ql, Q2 arrays do not depend on concentrations. The output of

this subroutine is used both in the main program blocks #25, #26, #28,

and #29 and subroutine CENTER (see blocks #1 and #2) . The output of

the subroutine is PI (J), Dl(J) for J = 2, NL - 1 and P2 (J) , D2(J) for

J = NR + 1, . . . , N2 - 1.

Block #1: Call Subroutine INTER

By calling this subroutine, we input A2L, AIM, B2L, BIM

through COI-IMON. These quantities are used in blocks #3 and #4 here.

Block #2: Calculation of Pl(l) when NL > 1

When NL > 2, this value of Pl(l) is used in the next block. Calculation formulas of the above quantities are found in eqs (6.7) and (6.7')

(This is a duplication of Pl(l) calculated in block #24 of the main pro

gram and is not necessary in VERSION 1.)

Block #3: Calculation of Arrays Pi (J) , Ql(J)

The calculation formulas for these quantities when J = 2, NL - 1 is

eq (6.9'). These quantities are used in blocks #25 and #29 of the main

program.

63
Block #4: Calculation of Arrays P2 (J) , D2(J)

The calculation formulas for these arrays when J = NR + 1, . . .

, N2 - 1 are given in eq (6.16'). Notice that P2 (N2) is calculated in block #13

of the main program and used here.

7.4. Description of Subroutine INTER

This subroutine is so trivial we only describe it briefly.

The routine uses the value of DT, DZ, DY, DFl, DF2 to compute the value

of A2L, A2M, A2R, B2L, B2M, B2R, AIM, BlM for use in subroutine TRIDNL

(see blocks #3 and #4)

Since the values of A2L, BlM depend on DT, every time this quantity changes these quantities need to be recomputed. This explains why

subroutine TRIDNL is called whenever DT changes in the main program,

i.e., because TRIDNL calls INTER (see block #1 in TRIDNL) and thereby

corrects the values of A2L, BlM for the new DT and, consequently,

simultaneously corrects the values of the arrays PI, Dl, P2, D2 for the

values of the new DT.

64
8. FLOW DIAGRAMS

8.1. Flow Diagram of Main Program in DISTRIB

#1

#2

#3

#4

#5

^ START ^

r DECLARATIONS

,

DIMENSION, AND

COM>'ON BLOCK- STATEMENTS

r INPUT DATA

AND FORI^iAT STATEMENTS

f

STORE INPUT

IN ORIGINAL

PHYSICAL UNITS

f

REDIMENSION

THE INPUT

DATA

COMPUTE THE MAGNITUDES OF

THE GRID MESH

V;iDTHS

A

65
0

#6

#7

#8

COMPUTE THE

VALUE OF CNSTX

= ZCRIT

COMPUTE DT FOR

THE FIRST

ITERATION OF THE

MAIN LOOP 78

COIJVERT DT

TO HOUR UNIT

AND STORE

CALCULATION OF THE

VALUES OF GBZD

INDICES NL, NR,

AND. N2 FOR THE

FIRST ITERATION OF

THE MAIN LOOP 78
©

#10

f

SET

NH +NH = NR

f

DO 96,

J --- KR, NH

C2 (J) = CMAX

DO 97, J = NH

C2(J) =+ 1, N2

CBULK

r

SET CB(2) = CMAX

1 f

67
CALCUTJ^TE TTMM

FOR FIRST ITERATION OF THE

MAIN LOOP 78

#12

CALCULATION OF SM(1)

SM(2) FOR THE FIRST

ITERATION OF THE

MAIN LOOP 78

#13

#14

CALCULATION OF

THE "P,Q" COEFFICIENTS FOR

GRID POINT N2

BEGIN MAIN LOOP

78 TO FIND CONCENTRATIONS AT

TIME TTMMPl,

SET CONSTANTS IN LOOP

i 68
#15

#16

SET PRINa-OUT

PAFii^METER ia-J T

COMPUTE DT

TIME /vCCELERATION SIGNAL

ZOACEL

YES

NO

r

UPVALUE SIGNAL

KPOV;, KPOW = KPCW + 1

<^ KPOW > 100

NO

YES

INCREASE DT. CALL SUBROUTINE TRIDNL.

DL-VALUE KPOW BY

100.

RESET 20PS -4

69
#17

UPVALUE NL AKD; FNL BY 1. SET

CI (ML) = CB(1)

#18

CAiyjUI^TE

NLMl NI, - 1 4

NLH2 = NL - 2

NO

YES

SET DT = DTt CALL TRIDNL

SET SM(1) =

CALCUIJ^TE TTMIIPI FOR tmth ITERATION

OF xMAIN LOOP 78 AND

.POSITION OF nOVING

BOUNDARIES AT TTMMPl

V 0

70
0

#19 CALCULATE

THE VALUE OF

SES FOR t-lMth ITERATION OF

LOOP 78

YES

SET NEV? VALUES

OF NR, JNR AND

SMi2)

SET VALUES OF

NR + 1, NR +

2 AND NRMl

#20 CALCUIJ^TE

QUA-NTITIES

NEEDED FOR

EVALUATING

TJ
#21

1

CALCUIiATx-- VALUE OF

TJ

YES

•CALCULATE NEW

VALUES OF ZOPl,

YOPl, DT AND

CALL TRIDNL

CALCULATE VAL- UES OF S (1) , S (2) FOR tmth

ITERATION OF , MAIN LOOP 78

'

72
©

#22 CALL TRIDNL

FOR USE IN THE FIRST

ITEPJ\TION OF

LOOP 78

#23

#24

<r NL > 1 ^

NO

f

CALCULATE

P] (NL) , gl (ML)

WHEN NL = 1

<^ NL > 1

YES

CALCULATE

PI (NL) , Ql (NL)

WHEN NL > 1

YES

NO

73
0

#25

#26

#27

CALCULATE

ARRAYS AF(J)

Ql(J)

CALCULATE

ARRAYS BF(J)

AND Q2 (J)

CB(1) , CB(2) , CI (NL) , C2 (NR) ARE CALCULATED

IN SUBROUTINE

CENTER, SEE

SECTION 7.3

YES

#28 CALCULATE

C2(J) J = NR

+ 1, TO J = N2 - 1 AT

TTIMPl

74
#29

#30

YES

1

NO

r CALCULATE Cl (J> J = 1 tc J = ML - 1 AT TIM£

ITIlMPl

r

YES

f PRINT Cl (J)

,

C2 (J) , CB (1)

,

CB (2) , KL, NR,

MM, DTS , TTMMP].

r

NO

75
#31

#32

#33

PRINT AND PLOT

CONCENTRATIONS

AT PRESENT VAL- UE OF TTMMPl

8 <^ MM. > KME

YES

r

DIAGNOSTIC

INCREASE

WME

UPVALUE MM AND

PERFORM NEXT

ITERATION OF

i'lAIN LOOP 78

76
0

#35

r

PRINTOUT

ARRAY

PARAM (J)

#36

#37

GO TO 73 X©—

f END

j

77
8.2. Flow Diagram for Subroutine CENTER

^ START ^

#0

r

INSERT INPUT

THROUGH

COMt-ION

#1

COMPUTE

P(2), Q(2)

#3

COMPUTE <

R(l) , R(2)

r

78
0

DO 850, L = 1,2

YES

#4

COMPUTE

PA{L)

QA(L)

YES

i

SET FA(1) = Pl(l), pA(l)

#5

CALCULATE

Cl(NL)

,

C2 (NR)

r ^ RETURN ^

79
8.3. Flow Diagram for Subroutine TRIDNL.

START

#1

CALL

INTER

#4

CALCULATE

P2(J), C2(J)

<

80
9. HOW TO USE DISTRIB

9.1. A Brief Description of DISTRIB: VERSION 1

The intent of this section is to give a quick general introduction to

DISTRIB: VERSION 1.

The meaning of the undefined terms in this chapter is explained in the

glossary, chapter 12. A more detailed description of the physical as- pects of the redistribution problem is given in chapter 2. Chapter 3 contains a more detailed description of the overall computational procedure.

DISTRIB: VERSION 1 calculates a redistributed initial impurity profile

in a crystal wafer that has undergone an oxidation at a constant temperature for an input time t = TFINAL. The impurity concentration in

the oxide Ci (2,t) is calculated at time t = TFINAL. Here, z denotes

(fig. 3) the distance measured into the oxide from the interface between the ambient oxygen and the oxide at time t = TFINAL. (The z- coordinate system changes with time.) Also calculated at time = TFINAL

are the impurity concentrations in the silicon C2(y/ t) where y represents (fig. 3) distance into the oxide-silicon composite measured from

the position of the oxygen-silicon interface at time t = 0. At this

initial time, the silicon is completely unoxidized. The principal output quantities are the concentrations CI (L) = C^ ( (L - 1)DZ, TFINAL) for

L = 1, NL where NL represents the index of the first grid point to

the left of the position of the moving oxide-silicon interfacial boundary, Zo(t), at t = TFINAL (fig. 3). Also computed are C2(K) = C2 ( (K - 1)

DY, TFINAL) for K = NR + 1 , . . .

, N2 where NR is the index of the first

grid point to the right of the moving boundary, yo(t) , at time t - TFINAL

(fig. 3) and N2 represents the number of grid points in the original silicon wafer. The quantities DZ and DY above are uniform grid widths in the

z- and y-coordinate systems. Besides the concentrations CI ( ) and C2 ( ),

DISTRIB computes CB(1) = Cj (Zq (TFINAL) ) and CB(2) = C2 (Yq (TFINAL) ) , the

impurity concentrations in the oxide and silicon at the moving oxidesilicon interface at time TFINAL. All of these concentrations are also

plotted

.

9.2. Definitions and Other Information Related to Input and Output

In section 9.21 the real input quantities are briefly defined. A table

follows the definitions giving further information about the previously

defined input, such as: the computational block number in the listing

that contains the read and format cards, format information, data

sources, and numerical ranges for which the parameters have been suc- cessfully tested, etc. This same procedure is used to describe the integer input in section 9.22. The same procedure is also used to describe

the real and integer output in sections 9.23 and 9.24. Consult the

glossary for information concerning undefined quantities and places in

the documentation where more thorough definitions are given. All physical quantities must be entered in units of hours and micrometers.

81
9.2.1. Real Input

The real input is defined in table 16 in the same order they are read

(see computational block #2 in the listing, chapter 11, for the read

and format card)

Table 16: Brief Definitions of the Real Input

Card 1 CBULK The program handles automatically piece-wise constant initial

concentrations in the silicon. The grid points in the silicon

beyond the input grid point NR are assigned the value CBULK

representing the concentration in the bulk of the silicon.

CMAX The grid points 1 through NH are assigned initial concentrations

of magnitude CMAX.

ALPHA A given voliime of silicon swells to a volume (ALPHA) times

(the original amount of silicon at the oxide-oxygen interface)

TFINAL The total time that the wafer is to undergo redistribution.

SEG The redistribution coefficient.

Card 2 DFl The diffusion coefficient of boron in the oxide.

DF2 The diffusion coefficient of boron in the silicon.

WFRC The fraction of a mesh width distance, DZ , that the moving

boundary moves in the first iteration of loop 78, i.e., in the

first time step.

BDRY A positive number representing the width of the wafer.

CONPRO A proportionality constant representing the evaporation rate of

boron from the oxide into the oxygen ambient.

Card 3 CONOUT A constant representing the concentration of boron in the oxygen

at the oxygen-oxide interface. (In our testing CONOUT was always set equal to zero.)

CNSTA

CNSTB

[

CNTAU

)

Constants appearing in the defining equation of the moving boundary, section 4.2.

YIS Same definition as Y2S, etc., below.

Card 4

Y2S

I Quantities that are used to scale the concentrations that are

X2S > outputted in the plotting routines, section 4.4, and computaXIS) tional block #34.

82
Runs of DISTRIB have been successfully completed using data in the ranges

shown in table 17. Sources for the data and formal information are also

given.

Table 17: Real Input (4 F0RI4AT (5 F 10.5)) (see block #2 in listing,

chapter 11) . Physical quantities are assumed to be in units

of micrometers

.

VARIABLE RANGE (Note 4)

(1) CBULK 1. (Note 4)

(2) CMAX 1. (Not tested with other values)

(3) ALPHA 0.44

(4) TFINAL 0.1 to 4.

(5) SEG

(6) DFl Qop> [41 rfii r7i [Ri

(7) DF2

(8) ^'JFRC \J * \J J- L-SJ W.JO ^i>J^^Uti

(9) BDRY o . uo o •

(10) CONPRO n 1 t-r-i 1 non

(11) CONOUT AllL J-vnnQ Lillo Vi^^7(^ LLCLVKIZ }~if^^n XJ^dl wiWXl^ii t'h ^K^d-\J 'Zf^T'O Vif^f^^ncio JOC^^QLiow

of our experimental conditions.

(12) CNSTA See tables in [1]

(13) CNSTB ^ : (14) CNTAU II

(15)

1 YIS 5.

(16) Y2S 2.

(17)

_ 1 X2S 4.

(18) XIS 2.

Notes : (1) The ranges in the subsections are values from successful runs

and are not to be interpreted as precise limitations.

(2) This relationship is cited in [4]

.

(3) WFRC must be chosen small enough so that DM = WRRC*DZ < ZCRIT. For an explanation of this restriction, see the

comment at the end of block #7 of the main program in chapter 7.

(4) The problem with CONOUT =0. is homogeneous in the concen- trations and therefore all runs have been made with CBULKl.

83
9.2.2. Integer Input

Table 18 contains brief definitions of DISTRIB's integer input.

Table 18: Brief Definition of Integer Input

Card 5 NZO A parameter used to determine the magnitude of DZ, eq (4.3).

MME The concentrations are calculated at times ti, tjyuyu^. MME

must be large enough for tjvjME > TFINAL or the program will stop

before calculating the concentrations at TFINAL. (See table and

sec. 9.3 for guidance in choice of MME.)

ND This integer determines how many grid points in the silicon adjacent to the moving boundary will be plotted.

NYO The number of mesh widths of length DY that equal BDRY.

NH The grid point on the y-axis at which the initial boron concentration changes from CMAX to CBULK, eq (2.10) and section 7.2,

block #10.

KPRINT The concentrations and a few other integer parameters such as

MM, NL, NR will be printed out, whenever MM = 0 modulo (KPRINT).

The range of values of integer data successfully calculated in DISTRIB

is given in table 19.

Table 19: Integer Input Information

Integer Input (1 FORMAT (616)) (see block #2 of listing)

VARIABLE RMvIGE

(1) NZO 2. to 8.

(2) MME up to 4000

(3) ND 20. to 30.

(4) NYO 100 to 800

(5) NH 10. to 20.

Certain constants that can be used to control various aspects of how

DISTRIB runs but do not have to be changed from run to run (and consequently, are not entered formally) are listed in the next table.

84
Table 20: Constants Used in Program But Not Made Part of Input.

(1) 0.2 and 100 used in block #16 in the listing of the main program

(see sec. 7.1., block #16).

(2) NMQD and NMODl (sec. 4.4.) and block #34 of listing of the main

program. These quantities need to be changed frequently to get

plotted points nicely spaced. We recommend starting with NMOD - 8,

NMODl = 4.

(3) BLTH (sec. 4.3). The value 0.02 micrometers has worked all the

time

.

(4) In printing out CI, C2 in blocks #30 and #34, certain indices are skipped. Of course, this is something one changes from run to run routinely, in order to find the values at the grid points that interests one.

9.2.3. The Real Output of DISTRIB

Table 21: Brief Definitions of the Real Output. See Tables 24 and 25.

DISTRIB calculates concentrations in the oxide and silicon for

a sequence of iterations of loop 78 with indices MI^l = 1,

I'lME . The time at the end of the MMth iteration of this loop

is TTi-UlPl (in units of BTM) . The concentrations of boron in the oxide and silicon at the

moving boundary

.

These are the concentrations in the oxide and silicon at the

Jth grid point (for the specific values of J see below) at

time TTMMPl.

Cl(l) The concentration at the oxide-oxygen interface.

See table 25.

PARAM is an array that contains material constants (diffusion

constants, etc.) in the original units of micrometers and

hours. SEG is the redistribution (= segregation) coefficient.

The swelling constant ALPHA represents the reciprocal of the

factor that a given volume of silicon swells when converted to

oxide at the moving boundary. TFINAL is the length of time of

oxidation. ZOCR is the position of the moving boundary in the

z-coordinate frame in micrometers. SCALE2 is a scale factor

used in the plots.

Further specific information about the previously defined output is

given in table 22.

PARAM

SEG

ALPHA

TFINAL

ZOCR

SCALE

2

85
Table 22: Real Output (5 FORMAT (7F14.8)) (see blocks #30 and #34).

Whenever MM = 1 modulo (KPRINT) , we print out

(1) C1(J) for J = 1, NL, 15.

(2) C2(J) for J = NR, N2 , 20.

(3) TTMMPl (in units of BTM)

(4) Cl(l). CB(1), CB(2), and DTS (in BTM units).

In addition, when TIMEND < TTMMPl, we print out

(5) Array PARAM, SEG, ALPHA, TFINAL, ZOCR, and SCALE2 (all in micrometers and hours)

9.2.4. The Integer Output

Table 23: Definitions of Integer Output. See tables 24 and 25.

NL is the index of the grid point in the oxide that is nearest

NL, NR the moving boundary grid point. NR is the grid point in the

silicon that is nearest the grid point at the moving boundary.

MM

The index of loop 78 that calculates the solution at time t

TTMMPl from the known solution at time t = TTMM.

See table 25.

NZO

MME

ND

NYO

NH

KPRINTJ

These quantities were defined in table 11

9.3. A Simple Example

In this section we illustrate how DISTRIB: VERSION 1 is used to com- pute the redistribution of an initially uniform distribution of boron

under wet oxidation conditions . The material constants used in this

sample calculation were chosen to obtain the best fit to the experimental data shown in figure 13. For further comments about the proper values for these material constants, see the end of this section.

The first card in computational block #2 reads the real input under format four. (This format card is also in computational block #2.) Since

the initial boron distribution is uniform, the value of CMAX which is

the concentration at the first NH grid points in the silicon (y-region)

is equal to CBULK which represents the initial concentration in the

86
bulk of the silicon and is the value assigned in the program to the

grid points NH + 1, NH + 2, . . . , N2 , Because the partial differential

equations are homogeneous in the concentrations when CONOUT (sec. 2.1)

is zero, as in this example, we normalize the input quantities CMAX,

CBULK by the value of CBULK (atoms per cubic micrometer); therefore,

these two input quantities have the value 1 and the output concentrations of these quantities should be multiplied through by CBULK to obtain the true concentrations. Once again, the glossary is for answering

these questions. The next parameter read is the swelling ratio ALPHA

(sec. 2.1) and it is assigned the value 0.44 [4].

We desire to find the concentration of the redistributed profile after

TFINAL = 0.3 h. When TTMMPl (T^mpi) is greater than TFINAL, the program's

control is shifted to the output printing blocks #34 and #35. The segregation coefficient SEG (dimensionless) is chosen to be 0.2, in order to

obtain the fit shown in figure 13. The diffusion coefficients DFl, DF2

for the boron in the oxide and silicon, respectively, were specifically

assigned the values 0.001 and 0.13 3 in order to obtain the fit shown in

figure 13, and these values are within the range of the reported values

for these quantities (see the end of this section) . (Note that the ratios of the diffusion constants used in the program, see [7] , were re- laxed in order to obtain the fit in this example.) The meaning of the

real grid input parameter WFRC = 0.5 is explained in section 7.1 in the

material concerned with block #7.

The meaning of the next parameters read, CONPRO (micrometers per hour)

and CONOUT (atoms per cubic micrometer) , are related to the flux of boron

between the ambient and oxide (sec. 2.1, eq (2.11)). As mentioned already, CONOUT = 0 in this example and CONPRO is chosen to have the value

1; other calculations show that results are not sensitive to the choice

of CONPRO including the value zero. The next three parameters read,

CNSTA (micrometers) , CNSTB (square micrometers) , CNTAU (hours) , relate

to the moving boundary motion (sec. 4.2, eq (4.5)) . The wet oxidation

is assumed to take place at 1100°C and therefore from table 1 in [1]

CNTAU(t) = 0.0, CNSTA(A) = 0.11, and CNSTB(B) = 0.51. The real input

read next, YIS , Y2S, X2S, and XlS (sec. 4.4.), are parameters that assign maximum abscissa and ordinates that are plotted in regions one and

two, i.e., the oxide and silicon regions, respectively, in the graphical

output in computational block #34. The values YIS = 5, Y2S = 2, X2S = 4,

and XlS = 2 have always yielded good graphical output, and we recommend

their use.

The integer input for DISTRIB is read following the real input in com- putational block #2 under format number 1 (also in computational block

#2). The integer NZO = 3, whose significance is explained in section

4.1 below, eq (4.2), is read first. The quantity MME = 1500 represents

the number of iterations we a priori assign to the loop 78 (computational

block #14). In this example, only 595 iterations of loop 78 occur before

"^I^MPl ^ TFINAL = 0.3 h and the program stops after printing the output

in computational blocks #34 and #35. Integer ND represents the number

of grid points in silicon to the right of the moving boundary point in

87
the silicon that are plotted in computational block #34 (see DO LOOP 234)

The integer NYO - 350 is read next and its function in defining DY(Ay)

is explained by eq (4.1). The integer NH = 12 (sec. 2.1, eq (2.13))

signifies that the initial concentration value CMAX is assigned to the

first twelve grid points in the silicon region at time zero. Since in

this example CMAX = CBULK, it makes no difference what the value of NH

is except NH should be less than NYO in order for the problem to make

sense. The print-out parameter KPRINT = 1000 signals (in block #30) the

program to print out CI ( ) , C2 ( ), the concentrations in the oxide and

the silicon, respectively, at the end of the iterations (1 + L*KPRINT)

of loop 78 where L = 0, 1, etc. Since the largest value of MM is,

for this example, 595, the only values of arrays CI ( ), C2 ( ) printed

out are those at the end of the first iteration of loop 78.

These are certain constant parameters that are not read into the program.

These are described in section 9.2 (table 20).

For various literature references for material constants , we recommend

[1,3,4,5,5,7,9,12,13] in the supplemental bibliography in chapter 13.

The output of DISTRIB is listed and described in section 9.2. The printout in block #30 is effected after every MM = (1 + L'KPRINT)th iteration

(for L=0, 1, 2, ...) of loop 78, which is the loop computing the con- centrations at the successive times t^M. Since KPRINT - 1000 for this

example and only 595 iterations of loop 78 are required to compute the

concentrations at time = TFINAL = 0.3 h in this example, the only

concentrations printed, under block #30, are those after the first iteration and are shown in table 24. From this table and computational block

#30, we see printed out the values of

CI (J) , J = 1, NL, 15

C2 (J) , J = NR, N2, 20.

Since, at time t^^ = 0.0389 (BTM units) , there is only one fixed grid

point in the oxide, the only value in the oxide that is printed is the

grid point at the oxide-oxygen interface. However, between the values

of Cl ( ) and C2 ( ) that we have printed out in table 24, there appear

the values of boron concentration at the moving boundary grid points in

the oxide and silicon CB(1), CB(2). Also printed in the same row with

these quantities are DTS (in BTM units) and the mesh width lengths in

the oxide and silicon, respectively, DZ and DY in units of BLTH = 0.02

ym. In the final print-out (table 26), BTM = 0.0030 h is printed. Also,

table 24 shows the last index MM of loop 78 and index NL of the first

fixed grid point to the left of the moving boundary in the oxide. Finally, in the same line appears the index NR of the first grid point to the

right of the moving boundary in the silicon.

The printout from block #34 is shown in table 25. The main output are

the concentrations

Cl (J) , J = 1, NL, 15

C2 (J) , J = NR, N2, 20. :. 88
These are the concentrations at the first time that TTMMPl (= 99.756).

is greater than TFINAL (in BTM imits) . It is seen from table 25 that

this occurs at the end of the MMth = 595 iteration. The meaning of the

other output in table 25 is similar to the corresponding output in table

24 that has already been described (see also sec. 9.23).

The concentrations at time TTMMPl = 99.756 (in BTM vmits) or TTMMPl ~ TFINAL = 0.3 h in the oxide are plotted in figure 11. The distances

from the oxide-oxygen interface (scaled in ZOPl units which is the

length of the oxide region at time TTMMPl ~ TFINAL) are plotted on the

X-axis. The concentrations at the positions of the first and (1 + L*NMOD)th grid points (for L = 1, 2, ...) are plotted on the y-axis, in

this example, NMOD = 8. The value of ZOPl (in micrometers) is computed

in the last table (26) . At the same time that the concentrations in the

oxide are plotted in figiare 11, the concentrations in the silicon are plotted in figure 12. The x-axis in this figure represents distance in

the unit UN2 = 2SQRT (DF2 'TTI-^MPl) (where TTMMPl is very nearly equal to

TFINAL) from the oxide-silicon moving boundary. The concentrations at

the moving boundary (in the silicon) and at the location of the grid

point with index NR are plotted. After this, the concentrations at the

grid points with the successive locations (NR - 1) DZ + (J - 2) 'NMODl'DZ,

for J = 3, . . . , ND (with NMODl =^ 4 in this particular example) are plotted. The lanit length of the abscissa, UN2 (micrometers), is com- puted in table 26 from computational block #35. In addition, in table

26 are shown all the input quantities in units of micrometers and hours

as printed out in computational block. #35. Besides the input data that

are printed out, the quantities ZO (TFINAL) (this is the same as ZOPl,

when TFINAL = TTMMPl) and 2 'SQRT (DF2 'TFINAL) in micrometers which, as we have already explained, are used as quantities to scale the x-axis

in the oxide and silicon concentration plots in figures 11 and 12. Also

printed is the value of BTM in hours; hence, the times in tables 24 and

25 which are in units of BTM, can be converted to hours, if so desired.

89
Table 24. Example of output of concentration in oxide and silicon at end of

iteration MM = (1 + L*KPRINT) in loop 78 for L = 0.

C I < . ) AT T IME

2. 52280046

NL

1 1 CI ( I ) 2 .522eC046

) AT TIME

. 9a4e'*473

1 .OOCCCOOG

1 .oocooooo

MM

C2( ,

.038^98095

CB ( I ) 3. 68795223

.033996095

1 . 00 cooooo

I. cooooooo 1.00000000

CB( 2) .77759045

1 .00000000

1 .00000000

1. ocoooooo

DTS

. 03899809

1 .00000000

1 . 00000000

I .00000000

DZ .541 12554

1 .00000000

1 .00000000

DV .71428571

1 .00000000

1 .00000000

I . 00000000

1 . 00000000

Table 25. Example concentrations in oxide and silicon when TTMMPl = TFINAL.

C l( . ) AT T IME .03302989

UH NL 595 32

C 1 ( 1 ) .03302989

C2< . ) AT TIME

. 189 53686

1 .00000000

1 .00000000

99.756364064

.8754271

NR

12 CB( 1 ) .813071 12 99.756364064

. 844 37386

1 .00000000

1. OCOOOOOO

, 83500469

CB( 2) .16261422

. 98733067

1 .00000000

I .00000000

. 82 160333

DTS

. 28899809

. 99958049

1 .00000000

.81678325

OZ

,541 12554 .99999438

1 .00000000

81453516

DY .71428571

.99999997

1. OOOOOOOO

.81331455

1 .00000000

1 . OOOOOOOO

Table 26.

ITA FOH STEP

OF 1 . 00 1 0000 0 SEG .20000000

CB< 1 ) .81 3071 1

CMAX

1 .OOOOOOOO

Nzc fue

3 1500

Listing of input physical data (micrometers,

computational data.

liours) and important

NO

30

1 DF2

. 1 33000 00 ALPHA .44000000

CB( 2 ) . 1626 1422 CBULK

1.00000000

NYO

350

CONPRQ

1 .OOOOOOOO

WFRC .05000000

FI RST+DT

. 0008691

VIS 5.00000000

KPR I NT 1000

CONOUT

.OOOOOOOO TFINAL

. 30000000

REAL TIME

. 30001914

y2S

2.00000000

CN5TA

. I 1 000000

ZO ( TF I NAL

)

. 34001 236 BDRV

5 .OOOOOOOO

X2S

4 .OOOOOOOO

CNSTB

. 5 1 000000

2*SQRT{ DF2*TF INAL

)

.39949969

CNTAU

, OOOOOOOO

XI S 2. OOOOOOOO

BTM .00300752

90
I I I I I I I t I I I I I I I I + I I I I I I I I I I I I I I I I + I I I I I I I I 4-

I I I I I I I I I I I I -f I 11111(1

91
I I I I I I I I I I 4- I I I I I I I • I I I I I 'f I I I I I I I I 'f

O I O I « I • III O t + I O I o * o I w I n I • IIIIII)II HCc0o

•H o 11 -1 1 • < 1 z r1 IL

V- 1 * i U. ! a 11 K o o 1 o 1 + « 1 o 1 o + o N. o 1 <0 o 1 • > 1 1 1 > 111 V) 1 >

•f1 IS 1 w 1 s 1 u 1 zo 1 o u 1 o : * o -1 oo in

wc

•H

cnCo

•H

+J

nsUPC

0)oaoo

(D

i-<&H

92
Figure 13 . Comparison of three calculated impurity concentration

distributions in silicon with experimentally determined concentration data.

93
10. VERIFICATION OF PROGRAM

10.1. Validation of the Program

A solution calculated by DISTRIB for a problem with the specified con- centration C2^(0,t) = Cq at the oxygen-oxide interface has been compared

with an analytic similarity solution given in [4] . The calculated boundary concentrations C2^(zQ(t), t) agree with the similarity solution [5]

to within a few percent for a variety of problems with differing segregation constants. The calculated solutions also have qualitative properties possessed by the similarity solution such as moving boundary con- centrations that are constant with time. There is also an analytic solution to special problems with zero flux at the oxygen-oxide interface.

In this problem, the concentration in the oxide or z-region is constant

and DISTRIB 's calculated solution is constant also for these conditions.

10.2. Limitations of the Program

The previous remarks must all be qualified for very short times, up to

5At [5] . The difficulties in obtaining the solutions for the first few

iterations are similar to those described in [2] for analogous problems.

The problem with the inaccuracy of the first few iterations can be some- what mollified, however, by choosing smaller and smaller values of WFRC.

Our recommendation for obtaining a solution at small times is to calculate solutions using a series of decreasing values of WFRC and then use those solutions to extrapolate the solutions to time zero.

94
11. PROGRAM LISTINGS

1. DISTRIB: VERSION 1 MAIN PROGRAM

1

5

41

1+38

9

26

27

36

32

39

4b

*+0 38

28

29

30

31

33

35

20

21

22

23

24

25

19

lb

START COVPUTATIOM BLOCK f^UMBEP 1 IMPLICIT POUBLE PRECISION (A-H»0-Z)

REAL XPL(200) »YPL(200)

DIMENSION TT(2000)»ZO(20no)»AF(200n)rRF(2ono)»PARAM(21)

DIMENSION BT(20) » DFA(2C) »nFB(?0) »CONPR(20) fCNST] {2Q) »CMSt?(2o) » 1 CNTA(2C) .TFIfNlN(20) »ZONEW(?n) rTNEW(20)

COMMOM/RL?0/DT»nY»DZ»THETA

C0r^M0N/PL21/DFl »DF2

C0'^M0N/eL22/DZP rPYP

C0MMGN/BL23/ 82M » RiiR » P2L » A2L » A2fw , A2R ' A I L » A 1 M » ftlP»PlL»RlMrPlP

C0MM0N/PL2'+/C2L»C2M»C2R»CIL»C1'^»C1R

C0MM0N/BL25/P(2) »0(2)»PA(2)»QA(2)»C(2)fCP(2)»P(?)»S(2)»SM{P)

C0MM0M/FL26/P1 (2000) »P2 (2000) fOl (2000) »Q2(?00n)»Dl (2000) 'D9(2nno)

C0N'M0N/PL2B/N2M1 »NLM1 rM2 f NRPl»NL»NR

C0MM0rVBL29/Cl (2000) C2 (2000)

COMMON /BL30 / CO f CK » TZERO » JEE » ALPHA f C AP A » SeG » CWY

COMMON /BL31 / PIST(100)» SOL (100)

C0MM0N/PL32/C0MPRO»C0N0UT

END COMPUTATION BLOCK ^lUMBER 1 START COVPUTATIOM~BLOCK NUMBER 2 ~ READ ( 5r4) CRULK »CMAX» ALPHA »TFlNALfSEG»DFl »DF2»wFPC»nDRY»C0NnR0»

1C0N0UT»CNSTA»CNSTR»CNTAU» Y1S» Y2S»X2S»X.1S

PEA0(5» 1 )NZ0»MME»MD»NY0»NH»KPRIr|T

FORMAT (616)

FORMAT

FORMAT

FORMAT

FORMAT

FORMAT (16H Cl(.) AT TIME »1F14.9)

FORMAT (16H C2(.) AT TIME »1F14.9)

FOPMAT( IHl ) FORMATCiOH TLIN CMTAU TT^^'v SLOPeZCRIT CnSTX )

(5F10.5)

(7F14.8)

( 4 OH MM ML NR

(40H MME MUST BE INCREASED

( 16H CI ( . ) AT TIME

) AT TIME

FORMAT

FORMAT

FORMAT

FORMAT

( 45H PROFILE IN OxIPE C/CB VS Z/ZO— STEP NO flI5 ) (IHl )

( 40H TIMENO TNEW(J) BTM CNTAU TLIN BT(J) ) (lOlH CKl) CP(1) CR{2)

DZ OY

COMPUTATION BLOCK NUMBER 2

( 52H Z0Pl»0F2»UN2»nYfnZfFNC» (NL-1) fTTMMpi,c;cALE2»XPSCUT

( 19H DATA FOP STEP NO »1I5)

33H NZO MME ND NYO NH JSTOP

43H S(l) ST STM El E2 E3 ERROR

)

)

)

) CONPRO CONOUT

DTS

) END

FORMAT

FORMAT

FORMAT(

FORMAT

(

FORMAT( 33H LIST OF Z0( ,

)

FORMAT( 33H LIST OF TT(.)

FORMAT (lOOH DFl DF2

1 CMSTA CNSTB CNTAU ) FORMAT{ 67H PROFILE IN SILICON C/CP VS ( Y-YO ) / ( 2*SQRT ( DF2*TF INaL)

1) — STEP NO »1I5 ) FORMAT (21H LIST AF ( . ) ) FORMAT (21H LIST PF(.) ) FORMAT ( 30H COXIDE CBULKZOl ALPHA TFINAL ) FORMAT (103H SEG PRINT DFl DF2 WFRC BDPY CMA^ CONPRO CONOUT

1 CNSTA CNSTB CNTAU YIS Y2S X2S ) FORMAT ( 33H LIST N2 NH Nl NR NL ) FORMAT( 33H LIST OZ DY OT TT(1) SM(1) 5M2) ) FORMAT (21H LIST PK.) 02(.) > FORMAT ( 33H LIST tJ2L B2P P2M PlL RlR ^IM )

95
17 FORMAT ( 33H LIST A2L A2M AlL AIR AIM ) 15 FORMAT( 20H LIST OF PlOfQK.) ) It^ FORMAT (60H VALUES OF RM R7 RY B C F CMVl CMV2 PIS P?5 AT m=

1 »1I9 ) 13 FORMAT (61H VALUES OF FNL 07P DYP TES ZOPl FN" TT(MMPl) SFS AT

1 MM= fllB)

12 FORMAT (28H DIVIDE FAULT AT MOV. RDRY. ) 11 FORMAT (28H DIVIDE FAULT AT LOOP 107 ) 10 FORMAT (2PH DIVIDE FAULT AT LOOP 105 ) 34 FORMAT (lOlH SEG ALPHA WFRC TFlNAL

1 ZO(TFINAL) 2*S0RT(nF2*TFINAL) ) 6 FORMAT ( 7F10.5 ) 37 FORMAT( 80H CBd) CB(2) FIRST+DT

1 BDRY

41 FORMAT (lOlH CMAX CBULK YIS

1 X2S

42 FORMAT ( 50H NH KPPINT

REAL TIME

) Y2S

)

)

CMAX CBULK

XIS BTM

NZO MME ND NYO

C DATA FIXED IN VERSION 1 DATA DFACl) »DFP(1) » CONPR ( 1 ) » CNSTl ( 1 ) » CNST2 ( D ^ CNTA ( 1 ) » TF IMN { 1)

1 /.0001D0» .OOOIDO^ .OlDOf . 027D0 » . IDO » . iDO » . 0 IDO/

TW=0,

ZW=0.

JSAVE=0

JSTOP =1

ZOl=0.

C0XIDE=1.

C JSTOP TELLS NUMBER OF DRIVE IN STEPS

JSTM1=JST0P-1

C START COMPUTATION BLOCK NUMBER 3

5CALE2 = 2.* SQRT(DF2*TFINAL)

C THESE MATERIAL PARAMETERS ARE RESCALED THUS WE HAVE TO STO^E

C THESE GIVEN PARAMETERS IN UNITS OF MICRONS ANP HOURS

PARAM(1)=DF1

PARAM(2)=nF2

PARAM{3)= CONPRO

PARAM(a)=CONOUT

PARAM(5)=CNSTA

PARAM(6)=CNSTB

PARAM(7)=CNTAU

C END COMPUTATION BLOCK NUMBER 3

C START COMPUTATION BLOCK NUMBER 4

BLTH= .02

C DEC MODIFIES PROGRAM TO DIMENEnTIONALISE ABOUT DFC+DF2

DFC= 1.

DF2=DF2*DFC

BTM=(RLTH**2)/DF2

C0NPR0=C0NPR0*(BTM/BLTH)

CNSTA=CNSTA/BLTH

CNSTB=CNSTR/DF2

CNTAU=CNTAU/BTM

BDRY=RDRY*(1./BLTH)

Z01=Z01* (1 ./BLTH)

DF1=DF1/DF2

DF2=1 ./DFC

TIMEND=TFINAL*(1./BTM)

C END COMPUTATION BLOCK NUMBER 4

C ZOl IS POSITION ON Z-AXIS WHERE THE NON DEGENERATE {Z0{2ER0)

C NOT EQUAL TO ZERO ) MOVING BOUNDARY BEGINS

C ZCRIT IS The position on the Z-AXIIS where the moving BOUNDARY

C CHANGES FROM LINEAR TO PARABOLIC

C CALCULATE DY»DZ»DT AND AUXILLARY QUANTITIES

Y01=ALPHA*Z01

96
C START COMPUTATION BLOCK MUI^BER 5

DY=BDRY/NYO

FNZO=FLOAT(NZO)

DZ=DY*(1 ./( ALPHA*FNZO)

)

C END COMPUTATION BLOCK NUMBER 5

C ================r==========r=r==r====r==========r=========r===r===

C START COMPUTATION BLOCK NUMBER 6

CNA=,5*CN5TA

cnb=cnstb

c take out following card if cnstx is known

CNSTX=SORT (CNA**2+CNB*CnTAU) -CNA

C END COMPUTATION BLOCK NUMBER 6

CNTaU=CNSTX* (CNSTX+CNSTA) /CNSTb

C START COMPUTATION BLOCK NUMBER 7

DM=DZ*WFRC

C INSERT IF MOVING BDRY IS LINEAR AND THEN PARABOLIC

TLlNrCNTAU

C IF MOVING BOUNDARY STARTS AT Y»0. THESE CARDS TAKE CARE OF STTUaTIO

C N THAT LINEAR PART OF MOVING BOUNDARY IS MOT PPESFNT

IF( TLIN .LE. .00000001 ) DTr ( DM* ( DM+CNSTA )) /CNSTB

IF( TLIN .LE. .00000001 )nTS=DT

IF{ TLIN .LE. .00000001 )ZCRIT=0.

IF( TLIN .LE. .00000001 )SL0PE=1.

IF(TLIN .LF. .00000001 ) GO TO 216

ZCRITrCNSTX

SLOPE=ZCRIT/CNTAU

DT=DM/SLOPE

DTS=DT

216 CONTINUE

C END COMPUTATION BLOCK NUMBER 7

C START COMPUTATION BLOCK NUMBER B PARAM( 10)=DT*BTM

C END COMPUTATION BLOCK NUMBER 8

C START COMPUTATION BLOCK NUMBER Q

C CALCULATE INITIAL INITIAL DISTRIBUTIONS

DUM=Y01/DY

NR=INT(DUM)+2

N1=INT{Z01/DZ)+1

NL=N1

N2=NY0+1

N2M1=N2-1

C END COMPUTATION BLOCK NUMBER 9

DO 98 J=1»N1

CI ( J)=COXIDE

98 CONTINUE

CB{1)=C0XIDE

C ===r===============r==r===r=====r=r===r========r=======rr===r===r=

C START COMPUTATION BLOCK NUMBER 10

N1P1=M1+1

NH=NH +rJR NHP1=NH+1

DO 96 J=NR »NH

C2( J)=CMAX

96 CONTINUE

DO 97 J=NHP1»N2

C2( J)=CBULK

97 CONTINUE

CB(2)=CMAX

C END COMPUTATION BLOCK NUMBER 10

97
START COMPUTATION BLOCK MUMBER 11

IF(Z01 .LF. ZCRIT) TTMM = 701/5L0PE

IF{Z01 .GT. ZCRIT) TTMM- ( ZOl* ( ZOl+CMSTA ) ) /CNSTR

END COI^PUTATION BLOCK ^)UlV!BER 11

- ------ START COMPUTATION BLOCK NUMBER

^^37 — — — — — r^ — — — — — — ^ — — ~ — FNL=FLOAT (ML) -1

FNR=FL0AT(NR)-1.

5M{1)=(Z01-FNL*DZ)/DZ

SM(2)=(FNR*DY-Y01 ) /DY

END COMPUTATION BLOCK NUMBER 12

START COMPUTATION BLOCK NUMBER 13

MAIN LOOP

P2(N2)=n.

Q2(N2)=C2(rJ2)

END COMPUTATION BLOCK NUMBER 13 ================================

START COMPUTATION BLOCK NUMBER

— — — — — — — — — — ~ — — ~ — Z. — —~— — — — — — — ~ — — — — — —

ZOPS=0,

ZOP1=0.

TJ=0.

N0RUN=1

KPOW=0

JSAVE=0

DO 78 MM=1»MME

END COMPUTATION BLOCK NUMBER 1*+ START COMPUTATION BLOCK NUMBER 15

KW=MOD(MMrKPRINT)

END COMPUTATION BLOCK NUMBER 15

^ —ZZ~^ — ~ — — — ~ — — ~~~~ — ~ — — — — — ~ — — — — — — — START COMPUTATION BLOCK NUMBER

— — ———— ~~——~~—~— — — — ———————————————

16

SPEEDS UP MOVING BOUNDARY

ZOACEL=(ZOP1-ZOPS)/DZ

IF(ZOACEL .GE. .2 ) GO TO 88

KP0W=KP0W+1

IF(KPOW .GT. 100 )DTS=DTS+,05

IF(KPOW .GT. 100 )DT=DTS

IF(KPOW .GT, 100 )CALL TRTDNL

IF(KPOW .GT. 100 ) KPOW=KPOW-100

CONTINUE

ZOPS IS ZO(TTMM)

ZOPS=ZOP1

END COMPUTATION BLOCK NUMBER 16 ================================

START COMPUTATION BLOCK NUMBER 17

IF(TJ .GE. 1.) NL=NL+1

IF(TJ .GE. l.)Cl(NL)=CP(l)

IF(TJ .GE. 1. )FNL=FNL+1

NLM1=NL-1

NLM2=NL-2

IF(TJ .GE. l.)DT=DTS

IF(TJ .GE. 1.) CALL TRIDNL

IFCTJ .GE.l. ) SM(1)=0.

END COMPUTATION BLOCK NUMBER 17

START COMPUTATION BLOCK NUMBER 18

MMP1=MM+1

TTMMP1=TTMM+DT

IF( TTMMPl .LE. TLIN ) ZOPl = ( TTMMPl - TW)* SLOPE +ZW

IF(TTMMP1 .GT. TL IN ) Z0P1=SQRT ( ( CNA**2 ) +CNB* ( TTMMP l-TW ) )-CNfl+ZW

Y0P1=ALPHA*Z0P1

END COMPUTATION BLOCK ^IUMBEP 18

98
C START COMPUTATION BLOCK NUMBER 19

C CHECK IF NEW CELLS ARE ADDED

SES= ( FNR+DY-YOPl) /DY

IF(SES .LT. 0.) NR=NR+1

IF(SES .LT. 0.) FNR=FNR+1,

IF(SES .LT. 0.)SM(2)=1.

NRP1=NR+1

NRP2=NR+2

NRM1=NR-1

C END COMPUTATION BLOCK NUMBER 19

C START COMPUTATION BLOCK NUMBER 20

MMP2=MMP1+1

TTMMP2=TTMMP1+DT

IF( TTMMP2 ,LE. TLIN ) Z0P2=(TTMMP2 -TW )* SLOPE + ZW

IF( TTMMP2 ,6T, TLIN ) zOP2=SQRT { { CNA**2 ) +CNB* ( TTMMP2-TW )

) -CNA+zW

TJ= (Z0P2-FNL*DZ) /DZ

IF(TJ .LT. 1.) GO TO 219

Z0P1=(FNL+1. )*DZ

IF( ZOPl .LE. ZCRIT)TTMMP1=(Z0P1-ZW)/5L0PE + TW

IF( ZOPl .GT. ZCRIT)TTMMP1=( (ZOPl-ZW)*( ZOPl-ZW+CMSTA ) ) /CNSTB +TW

DT=TTMMP1-TTMM

CALL TRIDNL

Y0P1=ALPHA*Z0P1

219 CONTINUE

C END COMPUTATION BLOCK NUMBER 20

C CONSERVATIOfJ OF NUMBER CALCULATION

S1M=C1 (rJL)*DZ*.5+(Cl(NL)+CB(l) )*SM{1)*.5*DZ

IF(NL .EQ. 1) S1M=(C1<1>+ CB ( 1 H *SM ( 1 ) * . 5*0Z

IF(NL .GE. 2)51rv'=SlM+.5*nz*Cl (1)

IF(NL .LE. 2) GO TO 222

DO 111 J=2»NLM1

S1M=S1M+DZ*C1 ( J)

111 COrJTINUE

222 52M=C2(NR) *DY* . 5+ ( C2 ( NR ) +CB ( 2 ) ) * SM ( 2 ) * . 5*DY

DO 112 J=NRP1»N2M1

S2M=S?M+C2( J)*DY

112 CONTINUE

C ======r==========r====r================rr=======r=====r===r==r==r=

C START COMPUTATION BLOCK NUMBER 21

5(1 )=(Z0P1-FNL*DZ)/DZ

5{2)=(FNR*nY-Y0Pl )/DY

C END COMPUTATION BLOCK NUMBER 21

C =r=======r=====r==z==r===========r======r===r==rr==r==============

C START COMPUTATION BLOCK NUMBER 22

IF(MM .EQ. 1) CALL TRIDNL

C END COMPUTATION BLOCK NUMBER 22

C START COMPUTATION BLOCK NUMBER 23

IF(NL.GT, 1) GO TO 223

IF(NL .EG. 1)DM1=S(1)*DZ

IF(NL .EQ. l)DM2rDMl*DMl/(2.*DT)

IF(NL .EQ. 1 .AND. MM .EG. 1) DM2=0.

IF(NL .EG. 1)DM3=DM2

IF(NL .EG. l)DM=l./( DF1+0M2+DM1*C0NPR0)

IFCNL .eg. 1 )P1(1)=DM*DF1

IF(NL .EG. 1)Q1{1)=DM*(C1(1)*DM3+DM1*CONPRO*CONOUT)

C END COMPUTATION BLOCK NUMBER 23

99
C START COMPUTATION BLOCK NUMBER 24

223 IF(NL .6T. 1 ) DM=1 . / {DZ*DZ/ (2

*DT ) +DZ* CONPRO+DFl ) IF{NL .GT. 1 )P1(1)=DF1*0M

IF(NL .GT. 1 )Q1(1)=(+C1(1)*(DZ**2)/(2.*Ot)+DZ*CONPRO*COMOUT1*Dm

C ^ END COMPUTATION BLOCK NUMBER 24

C START COMPUTATION BLOCK NUMBER 25

IF(NL .LE. 2) GO TO 213

DO 203 J=2»NLM1

AF( J)=A1M*C1( J)

203 CONTINUE

DO 205 J=2»NLM1

JM1=J-1

Q1U) = (AF( J)-Q1( JM1)*A2L)/D1 ( J)

205 CONTINUE

C END COMPUTATION BLOCK NUMBER 25

IF(MM .EQ. KD) WRITE( 6»15)

IF(MM .EQ.KD) WR ITE ( 6 » 5) ( Pi ( J) » J=l f NLMl » 1

)

IF(MM .EQ.KD) WRITE(6»5) (Ql ( J) » J=l rNLMl » 1

)

IF( MM .EQ. KD ) WRITE(6,20)

IF(MM .EQ. KD) WRITEC6»5) (AF(J) »J=2»NLMlfl)

C ======================r==============r======r========r==r===r=====

C START COMPUTATION BLOCK NUMBER 26

213 CONTINUE

DO 204 d=NRPlfN2Ml

BF( J)=B1M*C2( J)

204 CONTINUE

NE=N2-NRP1

DO 207 J=1»NE

K=N2-J

KP1=K+1

Q2(K)=(BF(K)-Q2(KP1) )/D2(K)

207 CONTINUE

C END COMPUTATION BLOCK NUMBER 26

C ==============================r====r========r===========r===rrr=r=

C START COMPUTATION BLOCK NUMBER 27

CALL CENTER

C END COMPUTATION BLOCK NUMBER 27

C START COMPUTATION BLOCK NUMBER 28

DO 209 J=NRP1»N2M1

JM1=J-1

C2( J)=P2( J)*C2( JMl)-«-G2( J)

209 CONTINUE

C END COMPUTATION BLOCK NUMBER 28

C START COMPUTATION BLOCK NUMBER 29

IF(NL .LE. 1) 60 TO 215

DO 210 K=1»NLM1

J=NL-K

JP1=J+1

C1(J)=P1{J)*C1(JP1)+Q1(J)

210 CONTINUE

215 CONTINUE

C END COMPUTATION BLOCK NUMBER 29

C START COMPUTATION BLOCK NUMBER 30

IF(KW .EQ. 0 .OR. KW .GT. 1) GO TO 799

WRITE(6»32)

IF(KW .EQ. 1) WRITE(6»8) TTMMPl

IF(KW .EQ. 1 )WRITE(6»5) (CKj) »J=lrNL»l5 ) WRITE(6»41)

IF(KW .EQ. 1) WRITE(6»1) MM»NL»NR

100
WRITE(6»45)

WRITE (6» 5) Cl(l) »C6(1) »CR(2) »DT5»DZ»DY

IF(KW .EQ. 1) WRITE{6»9) TTMMPl

IF(KW .EQ. 1)WRITE(6»5) (C2(J)»J=NR »M2'?0)

799 CONTINUE

C END COMPUTATION BLOCK NUMBER 30

C START COMPUTATION BLOCK NUMBER 31

SM(1)=S(1)

5M(2)=S(2)

TTMMrTTMMPi

C END COMPUTATION BLOCK NUMBER 31

C SHIFT TO NEW MOVirJG BOUNH-VRY AND CALC NEW PA^AMETFRS

C START COMPUTATION BLOCK NUMBER 32

IF ( TIMEND .GT. TTMMPl ) GO TO 78

C END COMPUTATION BLOCK NUMBER 32

C ==r====r==r=r=====:r=r===rr==============r===r=rr======rz=r======

C START COViPUTATION BLOCK NUMBER 33

IF(JSTOP .EQ. 1 ) GO TO 900

C END COMPUTATION BLOCK NUMBER 33

MTIME=MTIME+1

TT(MTIME)=TTMMP1

Z0(MTIME)=70P1

C CALCULATE CQNeRVED QUANTITIES AND GRAPH Op STANDARD S DIMENSION

IF(KW .LT. 1) GO TO 800

AF(1)=Y0P1

DM=S(2)

Af(2)=AF(1)+DM DY

BF(1)=CP(2)

BF(2)=C(2)

DO 8O3 J=3»ND

AF(J)=AF(2)+( J-2. )*DY*1G.

K=NR +(J-2)*10

PF( J)=C2 (K

)

803 CONTINUE

DO 250 J=lrND

XRL( J)=AF( J)

YPL( J)=HF( J)

250 CONTirjUE

WRITE(6»26)

CALL PLOT (ND»XPL»YRL)

C CONSERVATION OF NUMBER CALCULATION

Si =C1(NL)*DZ*.5+(C1(NL)+CR(1) )* S(l)*.5*nz

IF(NL .EQ. 1) SI =(C1(1)+CP(1) )*S(1)*.5 *DZ

IF(NL .GE. 2)S1 =S1 +.5*DZ*C1(1)

IF(NL .LE. 2) GO TO 221

DO 113 J=2»NLM1

51 =S1 +DZ*C1(J)

113 COUTirJUE

221 S2 =C2(riR) *DY*.5+(C2(NR)+CB(2) )* S(?)*.5*DY

DO 11*+ J=MRP1»N2M1

52 =52 +C2(J)*DY

114 CONTINUE

5TM= S1M+S2M

ST=S1+S2

E1=ST-STM

RY={DF2*DT) / (DY+*2)

E2=-DT*(C0rjPR0) * (CI (l)-CO^lOUT)

E3=RY*{C2(N2)-C2(N2Ml) )*DY

ERR0R=E1-(F2+E3)

101
C START COMPUTATION BLOCK NUMBER 31

800 CONTINUE

WRlTE(6r26)

WRITE(6f8) TTMMPl

WRITE(6»5) (C1(J) »J=1»NL» 5

WRITE(6»41)

WRITE(6fl) MMfNLrNR

WRITECeft+S)

WRITE (6» 5) CI (1) »CB(1) »CB(2) »DTSrDZ»DY

WRITE(6»9) TTMMPl

WRITE(6»5) (C2(J)»J=NR »N2»20)

NM0D=8

NM0D1=U

UN1=1./Z0P1

UN2=2.*SQRT(DF2*TTMMP1)

UN2=1./UN2

EZ=DZ*UN1

NLP1=NL+1

NLP2=NLP1+1

NLP3=NLP2+1

AF{1)=0.

BF(1)=C1(1)

NC=INT( (NL-D/NMOO)

FNC=FLOAT(NC)

IF(NC .LT. 1 ) GO TO 231

DO 230 J=1»NC

K=J*NMOD +1

JP1=J+1

BF( JP1)=C1(K)

AF( JP1)=(K-1)*EZ

230 CONTINUE

231 NCPl=NC+2

NCP2=NC+3

NCP3=NC+i+

AF(NCP1)=Z0P1*UN1

BF(NCP1)=CB(1)

BF(NCP2)=Y1S

AF(NCP2)=X1S

AF(NCP3)=0,

BF{NCP3)=0.

DO 252 J=1»NCP3

XRL( J)=AF( J)

YRL(J)=BF(J)

252 CONTINUE

WRlTE(6r26)

CALL PLOT (NCP3rXRL»YRL)

WRITE (6» 36) NORUN

EY=DY*UN2

AF(NLP2)=0.

BF(NLP2)=CB(2)

BF( NLP3)=C2(NR)

AF(NLP3)=AF(NLP2)+S(2)*EY

DO 232 J=3»ND

K=NLP3+J-2

KS=NR+( J-2)*NM0ni

BF(K)=C2(KS)

AF (K ) =AF (NLP3) + ( J-2) *EY*NM0D1

232 CONTINUE

DO 234 J=1»ND

K=J+NLP1

AF( J)=AF(K)

BF(J)=BF(K)

102
23** CONTINUE

X2SCT=X2S/2.

NDE=0

DO 235 J=lrND

IF(AF(J) .LT. X2SCT) NDE=NDE+1

235 CONTINUE

NDP1=NDE+1

AF(NDP1)=X2S

BF(NDP1)=Y2S

NDP2=NDP1+1

AF(NDP2)= 0.

BF{NDP2)=0.

DO 251 J=1»NDP2

XRL( J)=AF( J)

YRL(J)=BF(J)

251 CONTINUE

WRITE(6»26)

CALL PLOT (N0P2»XRL»YRL)

C END COMPUTATION BLOCK NUMBER

C START COMPUTATION BLOCK NUMBER 35

WRITE (6f35) NORUN

PARAM{8)=CP(1)

PARAM(9)=CB(2)

PARAM(10)=DT*BTM

PARAM (11) =TTMMP1*BTM

PARAM( 12)=BDRY*BLTH

Z0CR=Z0P1*BLTH

WRITE(6»32)

WRITE(6»38) NORUN

WRITE(6»33)

WRITE(6»5) (PARAM( J)

»

J=1.7 ) WRITE (6»34)

WRITE(6»5) SEG» ALPHA»WFRC» TFINaL»Z0CR »scale?

WRITE(6»37)

WRITE(6»5) {PARAM(J) » J=fl»12 ) PARAM{12)=CMAX

PARAM( 13)=CBULK

PARAM(14)=Y1S

PARAM(15)=Y2S

PARAM(16)=X25

PARAM(17)=X1S

PARAM(18)=BTM

WRlTE(6r4'+)

WRITE(6»5) (PARAM(J)» J=12»18 ) WRlTE(6»'+2)

WRITE(6»1) NZO»MME »ND »NYO»NH tKPRINT

C END COMPUTATION BLOCK NUMBER 35

C EtaRT COMPUTATIOn'bLOCK NUMBER 36

J5AVE =JSAVE +1

IF( JSTOP .EG). JSaVE ) GO TO 79

C END COMPUTATION BLOCK NUMBER 36

J = JSAVE

DF2= DFB(J)*DFC

BT (J) =( BLTH**2)/ DF2

WRlTE(6»5)nFB{l) »DFC»DF2»3LTH»BT(

D

TNEW(J)= TTMMPl* BTM/BT (J)

CONPRO = CONPR (J)* ( BT (J)/ BLTH)

CNSTA = CNST1(J)/BLTH

CNSTB = CNST2(J)/DF2

CNTAU = CNTA (J)/BT (J)

DFl = DFA(J)/DF2

DF2 = l./DFC

103
WRlTE(6»5)nFB{l ) » DFC r DF?

rLTH » RT ( 1 ) »DFl,nFA(l)

CNA = .5 * CNSTA

CNB = CMSTB

CN5TX = SORT ( CNA**2 +CmP*CnTAU )-CNA

TLIN = CNTAU + TNEW(J)

: CALCULATE DT SLOPE ZCRIT

DM=DZ*WFRC

IF( CNTAU .LE. .00000001)DT= (DM*(Dm+ CnSTA) )/CNSTR

IF( CNTAU .LE. . 0000000 1 ) DT5 =DT

IF( CNTAU .LE. .00000001) ZCRIT = ZOPl

IF( CNTAU .LE. .00000001) SLOPE =1.

IF( CNTAU .LE. .00000001) GO TO 316

ZCRIT = ZOPl + CN5TX

SLOPE = CNSTX/CNTAU

DT = DM/SLOPE

DTS=DT

316 CONTINUE

ZONEW(J) = ZOPl

TTMM= TNEW(J)

N0RUN=N0RUN+1

WRITE{6»27 ) WRITE(6»5)TLIN»CNTAU»TTMM, SLO^E »ZCRIT» CNSTV

TW =TNEW(J)

TIMEND = TFINN(J) * l./BT (J)

WRlTE(6r39)

WRITE(6f5) TIMEMD»TNEW( J) » PTM

CnTAU » TLT N » BT ( J)

ZW= ZONEW(J)

SCALE2 = 2.* SQRT(DFB( J)*TFINN(J)

)

BTM=BT( J)

TFINAL=TFINN( J)

CALL TRIDNL

PARAM(1)=DFA( J)

PARAM(2)=DFB{ J)

PARAM(3)= CONPR(J)

PARAM(U)=CONOUT

PARAM(5)=CNST1 (J)

PARAM(6)=CNST2( J)

PARAM(7)=CNTA( J)

C ===r=====r===z=r===rrrr===r===rr=======-==r====r:r=r===:

C START COMPUTATION BLOCK NUMBER 37

78 CONTINUE

IF(TIMEND .GT. TTMMPl) WRITE(6»43)

79 CONTINUE

STOP

END

104
11.2. Subroutine CENTER

C START COMPUTATION BLOCK NUMBER 0

IMPLICIT DOUBLE PRECISION (A-H»0-Z)

COMMON/BL20/DT»DYf DZrTHETA

C0MM0N/BL21/DF1 »0F2

C0MM0N/BL25/P(2) » Q ( 2 ) » PA ( 2 ) » QA ( 2 ) » C ( 2 ) » CB ( ? ) » R ( 2 ) f S ( 2 ) » SM ( ?

)

C0MM0N/BL26/P1 (2000) »P2(2000) »Q1 (2000)Q2 ( 2000 ) r DK 2000 )» D? ( 2000

)

C0MM0N/BL28/N2M1»NLM1»N2 f NRPl»NL»MR

COMMON/BL29/C1(2000) »C2(2000)

COMMON /BL30 / CO » CK » TZERO » JEE » ALPHA » CAPA » SE6 » CMVY

C0MM0N/RL32/C0NPR0»C0N0UT

C END COMPUTATION BLOCK NUMBER 0

5 FORMAT (7Flt|.5)

30 FORMAT ( '+3H S(l) SU SUM FL FR ERR ) C CONSERVE CALCULATION

SUM= DZ* (CI (NL)+CB(1) )*.5 *5M(1) +DY*

1SM(2) *(C2(NR) +CB(2) ) * . 5+C2 ( NR ) DY* .

5

C rr=====r===============r======r=r======r===r===============r=====

C START COMPUTATION BLOCK NUMBER 1 P(2)=P2(NPP1)

0(2)=Q2(NRP1)

C END COMPUTATION BLOCK NUMBER 1 C ==rr=========:=r=====r=i:=============rr==rr===rr======r===rr===r=

C START COMPUTATION BLOCK NUMBER 2 IF(NL .EQ. 1) GO TO 84°

P(1)=P1 (NLMl)

Q{1)=Q1(NLM1)

C END COMPUTATION BLOCK NUMBER 2

C ===r======rr===rr=====r=r=rr======r=====r========r====~=========

C START COMPUTATION BLOCK NUMBER 3

8U9 R(l)= 2.*DT/(DZ**2)

R (2)=2.*DF?*DT/ (DY**2)

DO 850 L=l»2

IF(NL .EG. 1 .AND. L .EQ. 1 ) GO TO 850

C COMPUTE P AND Q FOR HALF CELL

DM=1. +5(L)

IF(L .EO. 2 ) GO TO 847

DM1=DF1+S(L)/R(L)

DM2=DF1/DM

DM3=DM2*S(L)

IF(L .EO. l)DMt|=(S(L)/R(L) )*C1(NL)

IF(L .EG. 1) GO TO 848

8i+7 DM1 = 1.+S(L)/R(L)

DM2=1 ./DM

DM3=S(L)/DM

IF(L .EQ. 2)DM4=(S(L)/R(L) )*C2(NR)

848 PA(L)=DM2/(DM1-DM3*P(L)

)

GA(L)=(DM4+Q(L)*DM3)/(DM1-DM3*P(L)

)

850 CONTINUE

IF(NL .EQ. 1)PA(1)=P1(1)

IF(NL .EQ. 1)QA(1)=Q1(1)

C END COMPUTATION BLOCK NUMBER 3

C START COMPUTATIOn'bLOCK NUMBER 4

C CALCULATION OF BOUNDARY DENSITIES

IF(NL .EQ. 1)SAVE= S(l)

IF(NL .EQ. 1)S(1)=1.

105
DM=1.+S(1)

D11=S(1)**2

D12=S{1)*S(2)

RT={DZ*DF2)/(DY)

DMl=Dll/R{ 1)+ (SEG*RT*D12)/R(2) +^^f^/n^

IF(NL .EG. 1)DM1=SAVE/R(1) + (SEG*RT*S(2) )/R(2.)

DM2=-Dll/R(l)-»-(l.-S(l) )*DF1

IF(NL .EQ. 1)DM2=-(SAVE/R(1)+DZ*C0NPR0)

DM3=-RT*D12/R(2)-RT*S(1)*(1,+1,/R(2)

)

DM4=D11*DF1/DM

IF(NL .EQ. 1)DM4=0.

DM5=RT*S(1)

DM6=S(l)*SM(l)*(Cl (NL)+C3(1) ) /R ( 1 ) +RT*SM ( 2 ) S ( 1 ) *cB ( 2 ) /R ( 2

)

1 +S(1)*(RT*SM(2)/R(2)+RT/R(2) )*C2(NR)

IF(NL .EQ. 1)DM6=DN'6+DZ*C0NPR0*C0N0UT

IF(NL .EQ. 1)S(1)=SAVE

DM8=DM2+DMt+*P(l)

DM9=DM3+DM5*P(2)

DM10=DM6+G{1 ) *DMt++DM5*Q(2)

DM11=DM10+QA(1)*DM8 +QA(2)*DM9

DM12=DM1-DM8*PA ( 1 ) -SEG*DMq*P A ( 2

)

CB(1)=DM11/DM12

CB(2)=SEG*CB(1

)

C ^ END COMPUTATION BLOCK NUMBER U

C START COMPUTATION BLOCK NUMBER 5

DO 855 L=i»2

C(L)= PA(L)*CB(L)+QA(L)

855 CONTINUE

Cl(NL)=C(l)

C2(NR)=C(2)

C END COMPUTATION BLOCK NUMBER 5

C CONSERVE CALCULATION

SU= DZ*S(1)*(C1(NL)+CR(1) )*.5 + .5* DY*

1S(2)*(C2(NR)+CB(2) ) +C2 ( NR ) *DY* .

IF{NL .GT.l) FL=-R(1)*(-(S{1)**2)* Cl ( NLMl ) / ( 5 ( D +1 . ) +Cl ( ML)

1*{S(1)-1. )+CB(l)/(S(l)+l. ) )*DZ *.5*DF1

IF(NL .EQ.1)FL=-DT *CONPRO* (C 1 ( 1

) -CONOUT

)

IF(NL .EQ. 1)FR=R(2)* ( C2 (NRPl ) -C2 (NR ) ) *DY* .

IF(NL .EQ. 1)ERR=SU-SUM-{FL+FR)

IF(NL .EQ. 1)60 TO 860

FR=R(2)*S(1)*{C2(NRP1)-C2(NR) ) *DY *.5

ERR=S ( 1 ) * ( SU-SUM ) - ( FL+FR

)

860 CONTINUE

RETURN

END

106
.3. Subroutine TRIDNL

OFORrlS TRIDNL

SUBROUTINE TRIDNL

IMPLICIT DOUBLE PRECISION (A-H»0-Z)

COMMON/RL20/DT » DY » DZ » THETa

C0MM0N/BL21/DF1 »DF2

C0MM0N/BL23/ B2M » B2R

»

B2L » A2L » A2M » A2R » AlL » A 1 M f A IP f B1 L ' RlM r R1

R

C0MM0N/BL25/P(2) »Q(2)»PA(2)»QA(2)»C(2)»CB(2)»o(2)»S(2)»5M(2)

C0MM0N/BL26/P1 (2000) »P2(2000) »Q1(2000)

»

q2 ( ?00n )» Dl ( 2000 )» D? ( 2000

)

C0MM0N/PL28/N2M1 f NLMl »N2 » NRPl»NLfMR

COMMON/BL29/C1(20 00) »C2(2000)

C0MM0N/RL32/C0NPR0»C0N0UT

5 FORMAT (TFll.S)

C CALC. P»S AND DUMMIES

C START COMPUTATION BLOCK NUMBER 1 CALL INTER

C END COMPUTATION BLOCK NUMBER 1 IF(NL .EQ. 1)0M1=S(1)*DZ

IF(NL .EG. 1)OM2=Dm1*DM1/(2,*DT)

IF(NL .EQ. l)DM=i./( DF1+DM2+0M1*C0MPR0)

IFCNL .EQ. 1 )P1 (1)=DM*D'='1

C START COMPUTATION BLOCK NUMBER 2 IF(NL .GT. 1 )DM=1./(DZ*DZ/(2.*DT) +D7* CONPRO+DFl ) IFCNL .GT. 1 )P1 (DrDFl+n-^

C r============r=rr=====r======r=r========r=r==r=========rr=r=======

C END COMPUTATION BLOCK NUMBER 2 C START COMPUTATIOrj BLOCK NUMBER 3

IF(NL .LE. 2 )G0 TO 106

DO lOS J= 2fNLMl

JM1=J-1

DKJ) =A2M+P1 ( JM1)*A2L

PI ( J)=-A2R/D1 ( J)

105 CONTIrjUE

106 CONTINUE

C END COMPUTATION BLOCK NUMBER 3

C START COMPUTATION BLOCK NUMBER 4

DO 101 J= 1»N2M1

K=rJ2-J

KP1=K+1

D2(K)=R2M+P2(KP1)

P2(K)=-1./D2(K)

101 CONTINUE

C END COMPUTATION BLOCK NUMBER 4

RETURN

END

107
Subroutine INTER

SUBROUTINE INTER

IMPLICIT DOUBLE PRECISION (A-H»0-Z)

IMPLICIT DOUBLE PRECISION (A-H»0-Z)

COMMON/BL20/DT»DY»DZ»THETA

C0MM0N/BL21/DF1»DF2

C0MM0N/BL23/ B2M r B2R r B2L f A2L » A2M » A2R r A IL » AIM » AlR » BIL » RIM » BIR

5 FORMAT (TFU.S)

RA1=DT/(DZ**2)

RA2=DT/(DY**2)

A2L= DFl

A2R=A2L

A2M=-(2.*DF1+1./RA1)

A1L=1.

A1R= AIL

A1M=-1 ,/RAl

B2L=1.

B2R=B2L

B2M=-(2.+l./(DF2*RA2) ) B1L=1.

BIR =B1L

BIM = -l./(DF2*RA2)

RETURN

END

108
12 . GLOSSARY

See eqs (5.6) and (5.7) and table 5.

An array that stores the inhomogeneous terms in the tridiagonal equations for the z-region (table 5)

Swelling ratio a [eq (2.5)].

See eqs (5.10) and table 7.

The point representing infinity in the y-region [eq (2.11^)].

The unit being used to dimensionalize lengths (sec. 4.3).

The unit being used to dimensionalize time (sec. 4.3).

An array that stores the inhomogeneous terms in tridiagonal

equations in the y-region [eq (5.10) and table 7].

Material parameters [eqs (2.11), (2.12), and (2.13)].

Arrays representing concentrations at the grid points in the

z- and y-regions, respectively. At the beginning of an iteration of loop 78, CI ( ) represents Ci ( , TTMM), and at the end

of an iteration, CI ( ) represents C^ ( , TTMMPl) . Likewise,

C2( ) . These represent C^ (zq (TTMMPl) , TTMMPl), i = 1, 2, i.e., the

concentrations in the oxide and silicon at the moving boundary.

The concentration in the bulk of the silicon and at y = 0

[eqs (2.13) and (2.19), Cg = CBULK]

A subroutine for calculating CI (NL) , C2 (NR) , CB(1), CB(2)

(sees. 7.2 and 8.2).

The value of C2(y,0) for y £ NH in eq (2.13).

By definition, A/2.

Another name for CNSTB.

The constants A and B, respectively, in eqs (4.5) to (4.9)

(table 2)

Defined in eq (4.9) and FORTRAN name for ZCRIT.

The constant t in eqs (4.8) and (4.9).

109
CONOUT

CONPRO

Dl

D2

DFl

DF2

DFC

DM1

I through > DM4 ' DM1 1 through > DM10 J DT

DTS

DY

DZ

EY I EZ ' FNL

FNR

FNZO

JSAVE

JSTOP

KPOW

KPRINT

LOOP 78

m

- A constant in eq (2.18) representing the boron concentration

in the oxygen at the oxygen-oxide interface (sec. 5.7);

CONOUT = Cout-

- A proportionality constant in eq (2.18) representing the

evaporation rate of boron from the oxide into the oxygen

ambient (sec. 5.7); CONPRO = Cp.

- An intermediate array used to compute the arrays Ql , PI

[eq (6.9') and sec. 7.3, computational block #3].

- An intermediate array used to compute arrays P2, Q2 [sec.

7.3, block #4, and eq (6.16')].

- The diffusion coefficient Di in the oxide or z-region.

- The diffusion coefficient D2 in the silicon or y-region.

- A normalizing factor set equal to one in VERSION 1.

- For these quantities used in computational block #3 of siibroutine CENTER, see eqs (6.I0M and (6.17').

- For these quantities used in computational block #4 of subroutine CENTER, see eqs (6.20), (6.21), and (6.22).

- The time mesh width At (sec. 4.4)

- A value of DT that is stored before the moving boundary zq crosses a grid point (sec. 7.1, computational block #20).

- The mesh width Ay in the y-region (sec. 4.1)

- The mesh width Az in the z-region (sec. 4.1).

- Scaled value of DY, DZ used in plot output.

- The real value of NL - 1. - The real value of NR - 1. - The quantity NZO floated, i.e., made real.

- A parameter used to shift control in the program (computational block #36, sec. 7.1). JSAVE = 0 in VERSION 1 to

begin with.

- In VERSION 1, JSTOP equals 1. - A signal to increase the value of DT when KPOW = 100 and

ZOACEL <_ 0.2 (computational block #16, sec. 7.1).

- The concentrations and certain other integer parameters are printed out when MM = 0 modulo (KPRINT)

- Each iteration of this loop calculates the concentrations at one time. At. (Later see computational block #14, sec. 7.1).

- Another symbol for segregation coefficient SEG [eq (2.2)].

110
This is the "DO" index in loop 78.

This integer specifies the number of iterations in loop 78

(computational block #14, sec. 7.1).

The value of NL at time t = 0 in VERSION 1 (computational

block #9, sec. 7.1)

The number of fixed grid points in the y-region, N2 = NYO +

1 [computational block #9, sec. 7.1, and eq (4.2)].

Tells how many y-grid points to the right of NR will be

plotted in output graph (block #34, sec. 7.1).

A loop limit index in block #26, sec. 7.1.

The point on the y-axis where the initial concentrations in

the silicon changes from CMAX to CBULK [eq (2.13) and sec.

7.1, computational block #10].

See end of section 4.1.

The z-grid points are plotted modulo NMOD (block #34, sec.

7.1) . The y-grid points from NR to the right are plotted modulo

NMODl (block #34, sec. 7.1).

In VERSION 1. NORUN = 1.

See end of section 4.1.

The number of mesh widths of length DY that cover (0, BDRY)

(sec. 4.1) . A parameter used in the definition of DZ [eq (4.3)].

These are the arrays whose functions are described in eqs

(3.5) and (3.7). The quantities with the exception of Pl(l) are computed in subroutine TRIDNL. See eqs (6.5'), (6.5'),

and (6.7') for Pl(l) and (6.9') for P1(J), J > 1. See com- putational blocks #23 and #24 for computation of Pl(l).

See eq (6.16') for the derivation of (6.16).

Alternate names for Pl(NLMl), P2 (NRPl) [see eq (6.18')].

Coefficient linking CI (NL) with CB(1) and C2 (NR) with CB(2)

calculated in subroutine CENTER, block #3 [eqs (6.19') and

(5.12')]. Function of this array is explained in section

3.1, steps #2 and #3.

An array containing various material input constants and

computer program constants of possible interest that are program output (computational blocks #2 and #35, sec. 7.1.

An NBS "in-house" plotting program available from NBS Applied

Mathematics Division. Documented in Appendix 2.

Everything to do with these arrays is explained in the same places as Pi and P2

Ill
Q(l)l

Q(2)»

R

R(l))

R(2)>

RAl)

RA2'

RT

S(l)i

S(2)!

SAVE

SCALE2

SEG

SES

SLOPE

SM(1))

SM(2)

TFINAL

TIMEND

TJ

TLIN

TRIDNL

TTMM

TTMMPl

- Alternate names for Ql(NLMl), Q2 (NRPl) [eq (6.18')].

- Another name for RT.

- These are quantities that are defined in eqs (5.35) and

calculated in computational block #3 of subroutine CENTER.

[See also eqs (5.25) and (5.16, 5.16')].

- Combinations of mesh widths used in subroutine INTER (tables

5 and 7)

- This combination of mesh quantities is defined in eq (5.35)

and appears also in eqs (5.38) and (5.34) (computational

block #4 of subroutine CENTER)

- See eqs (5.35), (5.29), (5.22), (5.15) and figures 5 and 6.

These quantities are computed in computational block #21,

section 7.1.

- See computational block #4 of subroutine CENTER; see remarks

at the end of section 5.9 before the SUMMARY.

- Essentially the same as UN2 , a scale factor used in plotting.

(See sec. 4.4)

- The FORTRAN symbol for m, the segregation coefficient [eq

(2.2) ] . - A signal which when negative implies the moving boundary yo

reached a new y-grid point during the previous iteration of

loop 78 (computational block #19, sec. 7.1).

- The speed of the moving boundary in dry oxidations [eq (4.8)]

- See the same equations and figures cited in description of

S(l), S (2) . These quantities are computed in computational

block #12 and #31, section 7.1.

- This is the time in hours that is inputted into the program.

- When TTMMPl is greater than TIMEND (which is TFINAL in BTM

units) , the program stops and control is shifted to printing

output (computational blocks #32 and #4, sec. 7.1).

- A signal which, when its value is greater than or equal to

1 , indicates that the moving boundary Zq has reached a new z-grid point during the previous iteration of loop 78-.

- The quantity T in dry oxidations (table 4.1); also called

CNTAU.

- The subroutine where arrays PI, Dl are calculated (sees.

7.3 and 8.3)

- The time at the beginning of the MMth iteration of loop 78

in units of BTM (computational block #18, sec. 7.1).

- The time at the end of the MMth iteration of loop 78 (in BTM

units)

112
This quantity is always zero in VERSION 1.

Scaling factors used in the plot output (computational block

#34, sec. 7.1, and sec. 4.4).

The fraction of DZ that the moving boundary Zq moves in the

first iteration of loop 78 [eq (4.4) and computational block

#7, sec. 7.1] . Quantities used in setting scales for plots of concentrations

(computational block #34, sec. 7.1, and sec. 4.4).

A cutoff for points with large abscissas (sec. 4.4).

An array for abscissas used in plotting (sec. 3.4).

The position at time t in the y-coordinate frame, calculated

by using eq (2.8).

Quantities used in setting scales for plots of concentrations (computational block #34, sec. 7,1, and sec. 4.4).

Always zero in VERSION 1.

Position in y-region of moving boundary at time, TTMMPl (computational block #18, sec. 7.1).

An array for ordinates used in plotting (sec. 3.4).

The position of the moving boundary in the z-coordinate

frame at time t (sec. 4.2).

Same as CNSTX.

Always zero in VERSION 1.

When the moving boundary is moving slowly, this parameter is

small and the value of DT is increased in order to speed up

distance the boundary moves in a given iteration of loop 78

(computational block #16, sec. 7.1).

ZOPl in units of micrometers (computational block #34, sec.

7.1) . Position in z-region of moving boundary at time, TTMMPl (computational block #18, sec. 7.1).

See computational block #20, sec. 7.1.

See computational block #16, sec. 7.1.

Always zero in VERSION 1.

113
13 . REFERENCES

13,1. Bibliography

1. Deal, B. E., and Grove, A. G. , General Relationship for Thermal

Oxidation of Silicon, J. Appt. Phys. 36_, 3770-3778 (1965) .

2. Fox, L., What Are the Best Numerical Methods in "Moving Boundaries

in Heat Flow and Diffusion," J. R. Ockendon, Ed, pp. 210-241 (Clare

don Press, Oxford, 1975).

3. Gelfand, I., and Lokutsievski , O. V., in Appendix 2 of Theory of

Finite Difference Schemes (authored by Godunov, S. R. , and Ryabenki

V. S.) (North Holland, 1964).

4. Grove, A. S., Leistiko, 0., Jr., and Sah, C. T. , Redistribution of

Acceptor and Donor Impurities During Thermal Oxidation of Silicon,

J, Appl. Phys. 35.' 2695-2701 (1961) .

5. Kraft, R. , Finite Difference Techniques for Diffusion and Redistribution Problems with Segregation Type Boundary Conditions, PvoG.

AICA Int. Sym, , Bethlehem, Pennsylvania, June 17-19, 1975, pp. 328-

333.

6. Kurtz, A., and Lee, R. , Diffusion of Boron in Silicon, J. Appl.

Phys. 31, 303-305 (1960)

7. Horiuchi, S., and Yamaguchi , S., Diffusion of Boron in Silicon

Through Oxide Layer, Jap. J. Appl. Phys. 314-323 " (1962) .

8. Prince, J. L. , and Schwettman, F. N. , Diffusion of Boron from Implanted Sources Under Oxidizing Conditions, J. Electroohem. Soo

121, 705-710 (1974) .

9. Ralston, A. V., and Wilf, H., Mathematical Methods for Digital Computers (John Wiley & Sons, New York, 1960).

10. Semiconductor Measurement Technology, W. M. Bullis, Ed, NBS Spec.

Publ. 400-1, pp. 9-12 (March 1974).

11. Semiconductor Measurement Technology, W. M. Bullis, Ed, NBS Spec.

Pi±>l. 400-4, p. 8 (November 1974).

12. Semiconductor Measurement Technology, W. M. Bullis, Ed, NBS Spec.

Publ. 400-17, pp. 11-13 (November 1975)

13. Semiconductor Measurement Technology, W. M. Bullis, Ed, NBS Spec.

Publ. 400-19, p. 25 (April 1976).

114
13.2. Supplemental Bibliography

1. Allen, W. , and Atkinson, P., Comparison of Models for Redistribution of Dopants in Silicon During Thermal Oxidation, Sotid-State

Electronics 16, 1283-1287 (1973) .

2. Atalla, M. , and Tamenbaum, E., Impurity Redistribution and Junction

Formation in Silicon by Thermal Oxidation, Bell System Tech. J. 39 , 933-946 (1960)

3. Fuller, C. S., and Ditzenberger , J. A., Diffusion of Donor and Acceptor Elements in Silicon, J. Appl . Phys. 21_, 544-547 (1956) .

4. Ghezzo, M. , and Brown, D. , Diffusivity Nummary of B, Ga, P, As and

St, in Si02, J. Electrochem. Soc. 120 , 146-148 (1973)

5. Ghoshtoyore, R. N. , Intrinsic Diffusion of Boron and Phosphorus in

Silicon Free from Surface Effects, Phys. Rev. B. 2_, 389-403 (1971).

6. Kato , T. , and Nishi, Y. , Redistribution of Diffused Boron in Silicon by Thermal Oxidation, Jap. J. Appl. Phys. 2^, 377-383 (1964) .

7. Kxirtz, A., and Yee , R. , Diffusion of Boron in Silicon, J. Appl.

Phys. 303-305 (1960)

8. Morgalit, S., Neugroschel, A., and Bar Lev, A., Redistribution of

Boron in Silicon After Two Oxidation Steps in MOST Fabrication,

IEEE Trans. Electron Devices ED-19 , 861-868 (1972).

9. Prince, J. L. , and Schwettmann, F. N., Diffusion of Boron from Implanted Sources Under Oxidizing Conditions, J. Electrochem. Soc.

121, 705-710 (1974)

10. Rodoni, M. , Buneman, 0., and Dutton, R. , Boron Redistribution After

Oxidation, Abstract #333 (Electrochemical Society, Fall 1976)

11. Shimakura, K. , Suzuki, T., and Yadoviva, Y. , Boron and Phosphorus

Diffusion Through an Si02 Layer from a Doped Polycrystalline Si

Source Under Various Drive-in Ambients, Solid-State Electronics 18_,

991-997 (1975) .

12. Williams, E. L. , Boron Diffusion in Silicon, J. Electrochem. Soc.

108, 795-798 (1961)

13. Wilson, P. R. , The Diffusion of Boron in the Si-Si02 System, SolidState Electronics 15, 961-970 (1972).

115
Appendix 1

FORMULAS USED IN DERIVATION OF DISCRETE ALGEBRAIC EQUATIONS

We collect here for convenience some standard finite difference approximations used in the discretizations effected in chapter five. The deriva

tions of these formulas can be found in a text book on numerical analysis

1. TAYLOR'S SERIES

f (x+Ax) = f (x) + f ' (x) Ax + Ax2 +. . .+ Ax" +. . . (Al)

2 nl

2. QUADRATURE APPROXIMATIONS

a. Midpoint rule

f (x)dx = f(^) (b-a) + 0[(b-a)3] (A2)

b. Right end point

f(x)dx = f(b)Ax + 0[(b-a)2] (A3)

c. Left end point

/ f(x)dx = f(a) (b-a) + 0[(b-a)2] (A4) a

d. Trapezoid rule

f(x)dx = ^^^^^^^^^ (b-a) + 0[(b-a)3] (a5)

3. DERIVATIVE APPROXIMATIONS

a. Forward differences

df (x) ^ f(x+Ax) - f(x)

dx Ax

+ 0[Ax] (A6)

116
Backward differences

df(x) f(x) - f(x-Ax)

dx Ax

Central differences

+ 0[Ax] (A?)

—

df (x) f (x+Ax)-f (x-Ax) ^ 2i :;

-, - —_ . + 0[Ax^] - . (A8) dx 2Ax

4. THREE-POINT APPROXIMATIONS

Let X2 , X, X]^ be three points such that xj > x > X2 ; then a three-point

approximation to the derivative of t(x) at point x is

,^ Ax. f(x-Ax„) (Ax -Ax^)f(x) df _ _ 1 2 _^ 1 2 ^ dx Ax^ (Ax^+Ax^) Ax^ Ax^

(A9)

Ax^ f (x-Ax^)

Ax^ (Ax^+Ax^)

where Ax^ = x^^ - x and Ax^ = x - x^

117
APPENDIX 2

Plotting vuh routines

B. L. Joiner and S. T. Peavy

Five FORTRAN subroutines for producing plots similar to the one shoun in Figjre 1 are now available on the NBS UNIVAC 1108

FASTRAND file named PLOTS.

Questions related to the operation of these subroutines may- be directed to:

Brian L. Joiner or Sally T. Peavy

A337 Administration Bldg.

National Bureau of Standards

Washington, D. C. 20234

921-2315

Task numbers for the IINIVAC 1108 may be obtained from the

Ccmputer Services Division

A238 Adninistration Bldg.

National Bureau of Standards

Washington, D. C. 20234

921-3364

General Remarks

Except as specified below, the suLx-outines automatically

figure out limits for the plots based on the smallest and largest

data points. A new page is not called by any of these

subroutines. This allows the user to label the top of the page

before calling for tl^e plot. These subroutines do not change

any of the values of tlie arguments. Plotting is done by

repetitively seardung the arrays rather than by sorting them.

The resolution of the plots is 51 characters high by 101

characters v.ide and each plot consumes 54 lines counting borders

and scale labeling. The length and v,ldth of the actual plotting

area are Sh inches and 10 inches respectively, and the overall

dimensions including borders and scale labeling are 9 inches

and 11'4 inches respectively.

118
These five subroutines may be called in from FASTRAND by inserting

the follov/ing pair of cards after the RUN card:

e XQT CUR

INF PLOTS

viiere @ means a 7 and 8 punched in card column one.

Subroutines

I. PLOT (N, X, Y)

Plots the data in one column versus that in another column.

N : -The number of points to be plotted.

X : A vector (one dimensional array) containing the

values of the abscissa.

Y : A vector containing the values of the ordinate.

Example of main program.

DIMENSION X(10), Y(10)

DO 2 I = 1,5

X(I) = I 2 Y(I) = SIN(X(I))

N = 5 CALL PLOT CN, X, Y)

STOP

II. PLOTS G^^GS, X, Y, NRMX, NROW)

Plots the data in up to 5 pairs of columns. The symbols

used are . * + , - respectively for the 5 curves.

NARGS : The number of curves to be plotted.

X : A matrix (t\co dimensional array) having up to

5 coluTjis and containing the values of the

abscissa.

Y : A matrix having the sa-ne dimensions as X and

containing the values of the ordinate. The

first colunm of Y is plotted versus the first

column of X, etc.

NRMX : A vector (one dimensional array) containing the

number of points to be plotted in each of the

(5 or less) column pairs.

NROW : The number of rows specified for X and Y in their

dimension statement (s) in the main program.

(An example is given on the following page.)

119
Example of main program.

DINENSION XC85,4), YC85,4), NRNKC4)

DO 5 I = 1,20

X(I,1) = I-IO

-5 Y(I,1) = ABS(X(I,1))

DO 10 I = 1,37

X(I,2) = 1-20

10 YCI,2) = 20. - X(I,2)**2

NARGS = 2 NR^IX (1) = 20

NRMX (2) = 37

NROW =85

CALL PLOTS OIARGS, X, Y, NRMX, NRO'.V) STOP

III. PLOTL (N, X, Y, IT)

Plots the data in one column versus that in another column

using the symbols specified in the IT column.

N : The number of points to be plotted.

X : A vector (one dimensional array) containing

the values of the abscissa.

Y : A vector containing the values of the ordinate.

IT : A vector containing numbers bet\veen 1 and 26

which dictate the letter of the alphabet to

be used as a plotting symbol.

number: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

letter: ABCDEFGHI J K L M N 0

16 17 18 19 20 21 22 23 24 25 26 PQRSTUVWXYZ

Example of main program.

DIMENSION X(IOO), Y(IOO), ITClOO)

DO 20 I = 1,5

X(I) = I Y(I) = SIN(XCI))

20 IT(I) = 19

DO 30 I = 6,10

X(I) = 1-5

Y(I) = coscxci))

30 ITCD = 3 CALL PLOTL (10, X, Y, IT)

STOP

120
IV. PLOTIA CN, X, Y, IT, YMINS, YMAXS)

This subroutine is the same as PLOTL except that the

limits for the Y axis (ordinate) are specified in the call

Statement. If any points fall outside the specified limits,

the limits are stretched so as to include all points.

Example of main program is same as for PLOTL except for call.

CALL PLOTLA (10, X, Y, IT, -1., 1.)

V. PLOTLF (N, X, Y, IT, XMIN, XMOX, YMIN, YMAX)

This subroutine is also similar to PLOTL except that here

the limits for both X and Y are specified in the call statement.

Any points falling outside. the specified limits are omitted but

a talley is kept and the number of offending points is printed

below the plot.

Exaiiple of main program is same as for PLOTL except for call.

CALL PLOTLF (10, X, Y, IT, 0., 5., -1., 1.)
o

« •

• *

* •

• «

o

oo

•

+o

o

• * *

•

oo

+o

« o

I

o

CO

•f I I t I oo

I t I I I fill o

o

s

Illlll-fllllXlllt-l-llll I I I I • I I I I I I I I I

o

I

in

I

o o o o

•4- -0 o vn

o> o

• •

•-4 •-<

1 1 Figure 1.

122
NBS.114A REV. 9-78)

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO.

Spec. Publ. 400-57

a.Gov't Accession NoJ 3. ftacifitfeftt's Accession No.

I

4, TITLE AND SUBTITLE Semioonduatov Measurement Teahnology: DISTPIB I, An

Impurity Redistribution Computer Program

5. Publication Date February 1979

%. Perfafmimg Or^nization Cod«

7. AUTHOR'

David Gilsinn and Richard Kraft

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, DC 20234

1ft, Pfojpcf/Task/Work Unit No.

11. Contract/Grant No. ARPA Order No. 2397

12. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS fSfreet, City, State, ZIP) ARPA, 1400 Wilson Boulevard, Arlington, VA 22209

13. Type of Report & Period Covered Final

July 1973 - Dec. 1975

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTESLibrary of Congress Catalog Card Number: 79-600002

I

' Document descrit)es a computer program; SF-185, FIPS Software Summary, is attached.

16. ABSTRACT M 200-word or leas tactual summary of most si gni licar^t information. If document includes a significant bibliography or

literature survey, mention it here,)

This report provides documentation of a computer program which calculates the re- distribution of impurities in silicon during a single oxidation step. The documentation provides: (1) a physical and mathematical description of the redistribution process, (2) a detailed description of the discretization of the appropriate

partial differential equations, and (3) a complete description of t±ie FORTRAN program for computing the solution.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of the first key word unless a proper name;

separated by semicolons) Diffusion; electronic technology; impurity distribution; material transport; segregation; semiconductor technology.

18. AVAILABILITY |_}^Unlimited

For Official Distribution. Do Not Release to NTIS

1 1 Order From Sup. of Doc, U.S. Government Printing Office, Washington, DC

20402, SD Stock No. SN003-003-0 2 030-3

(3 Older From National Technical Information Service (NTIS), Springfield, VA. 22161

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

21. NO. OF

PRINTED PAGES

130

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

22. Price $3.00

*U.S. GOVERNMENT PRINTING OFFICE : 1979 0-281-067/6

NBS TECHNICAL PUBLICATIONS

PERIODICALS

ifOURNAL OF RESEARCH—The Journal of Research

|>f the National Bureau of Standards reports NBS research development in those disciplines of the physical and

Engineering sciences in which the Bureau is active. These

Include physics, chemistry, engineering, mathematics, and

|:omputer sciences. Papers cover a broad range of subjects,

(With major emphasis on measurement methodology, and

lie basic technology underlying standardization. Also in- j:Iuded from time to time are survey articles on topics closely

(-elated to the Bureau's technical and scientific programs. As

^ special service to subscribers each issue contains complete

|dtations to all recent NBS publications in NBS and non- NBS media. Issued six times a year. Annual subscription:

{domestic $17.00; foreign $21.25. Single copy, $3.00 domestic;

53.75 foreign. Note: The Journal was formerly published in two sections: Section A "Physics and Chemistry" and Section B "Mathematical Sciences."

>DI MENSIGNS/ NBS

This monthly magazine is published to inform scientists, engineers, businessmen, industry, teachers, students, and

jconsumers of the latest advances in science and technology,

iwith primary emphasis on the work at NBS. The magazine

highlights and reviews such issues as energy research, fire iprotection, building technology, metric conversion, pollution

labatement, health and safety, and consumer product per- iformance. In addition, it reports the results of Bureau pro- grams in measurement standards and techniques, properties

.of matter and materials, engineering standards and services, linstrumentation, and automatic data processing.

Annual subscription: Domestic, $1 1.00; Foreign $13.75

' NONPERIODICALS

(Monographs—Major contributions to the technical liter- ature on various subjects related to the Bureau's scientific and technical activities. iHandbooks—Recommended codes of engineering and industri:il practice (including safety codes) developed in coopera1 with interested industries, professional organizations,

J regulatory bodies.

Special Publications—Include proceedings of conferences

^p>)nso^ed by NBS, NBS annual reports, and other special

publications appropriate to this grouping such as wall charts, pocket cards, and bibliographies.

Applied Mathematics Series—Mathematical tables, manuals, and studies of special interest to physicists, engineers,

chemists, biologists, mathematicians, computer programmers,

and others engaged in scientific and technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical properties of

materials, compiled from the world's literature and critically evaluated. Developed under a world-wide program coordinated by NBS. Program under authority of National

Standard Data Act (Public Law 90-396).

NOTE: At present the principal publication outlet for these data is the Journal of Physical and Chemical Reference Data (JPCRD) published quarterly for NBS by the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints, and supplements

available from ACS, 1155 Sixteenth St. N.W., Wash., D.C.

20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials, components,

systems, and whole structures. The series presents research

results, test methods, and performance criteria related to the

structural and environmental functions and the durability and safety characteristics of building elements and systems. Technical Notes—Studies or reports which are complete in themselves but restrictive in their treatment of a subject.

Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject area. Often

serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures

published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations. The purpose

of the standards is to establish nationally recognized require- ments for products, and to provide all concerned interests with a basis for common understanding of the characteristics of the products. NBS administers this program as a supple- ment to the activities of the private sector standardizing

organizations.

Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the consumer. Easily understandable language and

illustrations provide useful background knowledge for shopping in today's technological marketplace.

Order above NBS publications from: Superintendent of Documents, Government Printing Office, Washington, D.C.

20402. Order following NBS publications—NBSlR's and FIPS from

the National Technical Information Services, Springfield, Va. 22161.

Federal Information Processing Standards Publications (FIPS PUB)—Publications in this series collectively constitute the Federal Information Processing Standards Register.

Register serves as the official source of information in the Federal Government regarding standards issued by NBS

pursuant to the Federal Property and Administrative Services Act of 1949 as amended. Public Law 89-306 (79 Stat. 1127), and as implemented by Executive Order 11717

(38 FR 12315, dated May 11, 1973) and Part 6 of TiUe 15 CFR (Code of Federal Regulations).

NBS Interagency Reports (NBSIR)—A special series of

interim or final reports on work performed by NBS for outside sponsors (both government and non-government).

In general, initial distribution is handled by the sponsor;

public distribution is by the National Technical Information

Services (Springfield, Va. 22161) in paper copy or microfiche form.

BIBLIOGRAPHIC SUBSCRIPTION SERVICES

The following current-awareness and literature-survey bibli- ographies are issued periodically by the Bureau:

Cryogenic Data Center Current Awareness Service. A litera- ture survey issued biweekly. Annual subscription: Domestic, $25.00; Foreign, $30.00.

Liquified Natural Gas. A literature survey issued quarterly. Annual subscription: $20.00.

Superconducting Devices and Materials. A literature sirrvey

issued quarterly. Annual subscription: $30.00. Send subscription orders and remittances for the preceding bibliographic

services to National Bureau of Standards, Cryogenic Data

Center (275.02) Boulder, Colorado 80302.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington. D.C. 20234

OFFICIAL BUSINESS

Penalty for Private Use, $300

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE COM-21S

SPECIAL FOURTH-CLASS RATE

BOOK