December, 1982 Department of Agronomy and Soils Departmental Series No. 79 Alabama Agricultural Experiment Station Auburn U~niversity Gale A. Buchanan, Director Auburn University, Alabama 2 CONTENTS Page SUMMARY ...... .. ** ***** * S * 4 5 6 INTIRODUCTION................... * 4 NUMERICAL IMPLEMENTATION. .......... * 0 9 * 4 CALCULATION OF THE TIME STEP.... CALCULATION OF THE SPACE STEP*... RESULTS ........... ......... ..... 13 14 17 17 19 20 21 * 0 * 4 * le 1. Drainage from a Uniform Soil Profile..... 2. Steady State Infiltration into a Sandy Profile 3'. Evaporation and, Drainage in- a Uniform Soil Profi"I wi th Fi~ ed Water Table .. .. .. .. 0 4 e **0 0 DISCUSSION...............** ** LITERATU-RE CITED..............* * S 0 23 I 0 0 APPENDIX I: Approximation of Space Steps APPENDIX .'1: Executi on of WAFLOW * .... . 0 I S I I 24 29 0 0 0 S a 9 9 0 0 3 LIST OF TABLES AND FIGURES Page Table 1. Comparison of the Described Model WAFLOW with UNSATI for a 1.40-m Deep Drainage Profile Consisting of 29 Gridpoints . . . . . . . Table 2. Effect of the Number of Space Steps on Computer Time and Mass Balance during 8 hr. of Infiltration in a 70-cm Deep Profile (Mass Balance Criterion s is 0.001) . . . . . . . . . . . . . . . . . Table 3. Effect of Mass Balance Criterion on Computing Time and Mass Balance 65 hr. after Drainage Started . . . . . . . . . . . . . . . . 31 32 33 Fig. 1. Diagram Showing the Finite Difference Grid Superimposed on . . . . . . . . . .. the Depth-time Region of a Soil Profile Depth-time Region under Consideration about the General Fig. 2. Gridpoint (i,j) Showing the Identification of the Gridpoint Values of . 34 Pressure Head and Flux . . . . . . . .. . . . . . . . . . * . . . . 35 Fig. 3. Matrix Form for the System of Linear Equations as Presented . . . . . . . . . . . . . . . 36 . . . .0.. in Eq. [12] . . . . .e.0.. 2 Fig. 4. Water Content y and as a Function of Depth During 37 Infiltration (41 Gridpoints) . . . . . . . . . . . . . . . . . . . . . Water Content Profiles Calculated with the Models UNSAT1 Fig. 5. and WAFLOW. Number of Gridpoints: 29. Time from Start of Drainage = 14.4 hr. The Individual Data Points Indicate the Spacing in the z-Direction . . . . . . . . . . . . . . . . . . . . . . 38 4 LIST OF TABLES AND FIGURES (continued) Fig. 6. Water Content Profiles during Infiltration, 0.2019 and 0.4019 hr. since Start of Simulation. The Individual Data Points Indicate the Spacing in the. z-direction . . . . . . . . . . . . . Fig. 7. Water Content Profile after 4.0106 hr., when Evaporation Rate is 1/24 cm/hr. The Individual Data Points Indicate the . . . . . . . . . . . Spacing in the z-direction . . . . . . . . . Page 39 40 Appendix Table 1. Definition of the Main Program Variables in WAFLOW . . Appendix Table 2. Required Input Data and a Listing of the Actual . . . . . . Input Data for Example Problem 2 (Infiltration Profile) Appendix Table 3. Description and Listing of Output. Appendix Table 4. Listing of WAFLOW . 41 45 50 58 . . . . . . . . . . . . . . . . . . . . . ... . . . . Appendix Fig. 1. Illustrations for Selection of zi . . . . . . . . . . . 72 Appendix Fig. 2. Flowchart of WAFLOW. ....... . . . ... .... 73 5 SUMMARY The pressure head form of the general flow equation for water in a porous medium was numerically solved using a scheme that allowed both the time step (At) and the space increment (Az) to be changed during the flow process. mass balance equation -was used to check the accuracy of the simulation. The If a predefined accuracy criterion was not met, the time step as well as the space increments between a fixed number of gridpoints were changed. The model concentrates a large number of gridpoints in regions of large changes in pressure head gradients. As the number of gridpoints is fixed, fewer gridThe method is points are distributed in regions where small changes occur. demonstrated for an infiltration, drainage, and evaporation process. 6 AN ADAPTIVE SIMULATION TECHNIQUE FOR THE ONE-DIMENSIONAL WATER FLOW EQUATION J. H. Dane, J. W.L.M. Hopmans, and F. H. Mathis INTRODUCTION Many authors have numerically solved the one-dimensional pressure head form of the general flow equation for water in porous media. A comparison of In most the more popular procedures was presented by Haverkamp et al. (1977). cases presented in the literature, a fixed az value is used in the spatial domain. An adaptive numerical scheme to solve the pressure head form of the The general flow equation was recently published by Dane and Mathis (1981). difference between their model and most other numerical procedures is that both time and space increments are allowed to change during the simulation. This allows the use of a fine grid in those portions of the flow region known to contain large changes in pressure head gradients. The fine grid spacing moves up and down the soil profile, corresponding with those portions where large changes in pressure head gradient occur. Furthermore, their proposed In this report, is scheme is applicable to non-steady flow conditions as well. the source text, together with a description and discussion of the model, presented. THEORY The general transport equation for water movement in the soil which was used in the simulation described in this report can be formulated by combining *Associate Professor and Graduate Assistant, Department of Agronomy and Soils, Auburn Universityand Associate Professor, Department of Mathematics, Baylor University, Texas, respectively. 7 the Darcy equation and the mass balance equation of soil water. dimensional form the Darcy equation may be written as: In one v. -K(h) 1 where K(h) is the hydraulic conductivity function, h is soil water pressure head, H is hydraulic head, z is distance (z=O at reference level and z > 0 above reference level), and j- is the hydraulic gradient. The algebraic sign of the flux v indicates the direction of the flow, i.e., v is upward if positive and downward if negative. Since the hydraulic head is the sum of the may be written as: pressure and gravitational head, Eq. C(1 v = -K(h)fL + ah , 1 [2] where is the pressure head gradient. For a certain set of assumptions the mass balance equation can be written as a volume balance equation: a at 3av az [3] , of a This equation relates the time rate of change of water content, -differential volume element of soil to the difference of inflow and outflow across that element, ~z- , which is called the divergence of the flux. Com- bining Eq. [2] with Eq. [3] yields the general equation for vertical flow of water in soil: at az az [41 8 To obtain an equation for water flow in one dependent variable, relation between e and h is required. another Introducing the water capacity of the the time soil as C(h) =de (i.e., the slope of the water retention function) dh derivative term in Eq. (41 can be transformed to: a de ah ah [5] t: ~'at= C(h) at Substituting Eq. [5] in Eq. [43 yields: ah a= a C(h) K(h( +1) [61 which is called the pressure head form of the one-dimensional general flow equation. Equation [6] can be expressed as a system of two first order non-linear differential equations with h(z,t) and v(z,t) being the unknowns: v(z,t) = -K(h)(ah.+ 1) [7a] and C.(h)-- = --- z [7b] for -L - z 5 0 and t > 0, where v(z,t) is the Darcy flux, and L is the depth of the profile under consideration. An adaptive numerical method to solve this set of equations, subjected to certain initial and boundary conditions, will be presented. NUMERICAL IMPLEMENTATION The solution of Eq. [7] by a finite difference technique requires that a grid be superimposed upon the time (t)-depth (z) region of the flow system under consideration. The independent variables t and z will be subscripted with j and i, respectively. For-L _z .0 and 0 - t - T, the gridpoints (zit j ) are defined for i=O,1,...,N and j=O,1,...,M such that z0=0, 1). tjtj ZN=-L, t0=0, and tM=T (see figure i=z i - zi 1 The grid spacing may be arbitrary and is denoted by - tj , and i.e., both time and space increments may be of nonuniform size. However, the number of gridpoints in the depth direction is fixed. Space and time derivatives are approximated at time t. (unknown time level) and at the midpoint of (z i , zi- 1 ) using backward differences. Both K and C are evaluated at (h i j + hi, ) and vi, such that: For i=1,...,N and j=1,...,M we wish to find h i -K1 + = (v j + [8a] and C [hi j 2At where: + h i 1, - hi i,j-L -1,j h i i j -- . i i 1j [8b] K= K 2hi -1 , 3) and 10 C. 2 This method is implicit for which, at each time step, a system of (2N+2) nonlinear equations must be solved. To solve Eq. [8] one needs to define one set of initial conditions and two sets of boundary conditions. The initial condition for a given problem is specified by assigning a pressure head value to each node at t=O: h(z,O) = g(z) -L z 0 . [9a] This initial condition yields, according to Eq. [7a], the initial flux values: v(z,) = -KEg ( 3 z)+ . [9b] Boundary conditions are defined for the top and bottom nodal point for t,> 0 and are be expressed by: 1 h(Ot) + c2 v(0,t) = g (t) 0 z=O z=-L . [10a] [10b] a. h(-L,t) + a2 v(-L,t) = g (t) The advantage of expressing the boundary conditions in this manner is the ease with which one can switch from a pressure head to a flux boundary condition and vice versa by assigning the values 0 or 1 to c and 8. It should also be noted that the functions g0 and gl are time dependent. gridpoints, Eq. [9] and Eq. [10] can be written as: Expressed in terms of 11 hi so 0 g(zi) V and 0 OOK~g~zi)J L z.) 1i1a] ~h 1 02 +a.V~ g(~ 1lh $+2 vNj 1 t11.Cb) (g Eq. £8) can 'be solved by the foll1owi ng iterati ve procedure- (see fi gure 2 also). Sup-pose hi.i. and vij. are known for iO1.,N. If H12 and V2 are defined as the unknown -h and. v at time j, respectively, Eq. [8a] can then be wri tten as: --K2' 11 EV2(i) +V2,(i-1)] Equation [8b] Let 11(i), i=Q,...,N be the known pressure heads at time j-l. can then -be approximated by: 2At EH2(i) Gruig4 h H2(P'1) -wHl~ i li) Hl(i1)J ] 2 ils + [2i v(~) nkonA2zn and 12 -jC;R*H2(i-1) + V2(i-1) Az - CR-H2(i) - V2(i) =-C*R*HI* [12b] where R.: - At and H1* = Hl(i) + HI(i-1) . During the first iteration of each time step j, K and C are calculated from the known pressure head values HI at time j-l, while during subsequent iterations (all during the same time step), K and C are determined from the latest values of H2. If for all depths both H2(i) and V2(i), i=0,...,N differ less than a certain small amount from their corresponding values during the previous iteration, then the latest values for H2(i) and V2(i) are accepted as the true values. If .the difference between the previous and the latest computed pressure heads and fluxes is larger than that small amount, the solution process is repeated as before. Consequently, Eq. [123 represents a set of linear algebraic equations in the unknowns H2(i) and V2(i) for i=O,...,N with H2(0), H2(N), V2(0), and V2(N) satisfying the boundary conditions given in Eq. [11b]. Keller (1974) showed that H2(i) and V2(i) will converge to the solution of Eq. [8], provided that the pressure head and flux values over-a time step are suitably close. -At a certain time j, this system of linear equations can be put in the matrix form Ax=B, as shown in figure 3. In the simulation program, an IMSL (International Mathematical and Statistical Library) subroutine is used to solve this linear system. However, this subroutine requires that the coefficient matrix A is stored in band storage mode to minimize memory requi rements. A n by n matrix with k lower codiagonals and 1 upper codiagonals, stored in band storage mode, is reduced to a matrix of dimensions n by (k+l+1). The matrix is stored row wise so that the zero-elements are compressed out of the . matrix and the main diagonal elements fall in column (k+1). The 2N+2 by 2N+2 13 coefficient matrix A, with two upper and lower codiagonals, can be stored in band storage mode as shown below: 5 columns 0 0 2 K 0 -K a IAzl 0 iAz 1 -iC-R (2N+2) rows 0 -C*R 0 1 -C.R -1 0 -KCAz 1 1 N K -1 0 iAzN 0 0 -4C'R 2 The solution process uses the initial condition (Eq. [11]) to solve for the This set of pressure heads and fluxes at the end of the first time step#. solutions is then used to solve for the next time step. thus marches forward in time by increments at. The solution process CALCULATION OF THE TIME STEP The mass balance equation was used to control the time step size. large time step will result in an inaccurate approximation of h (see Eq. A too and v [8]), which in turn will influence the outcome of the mass balance equation. At time t., the mass balance MB. is defined as: MB =LP [e(z,tj) - e(z,tj- 1 )]dz -J" 3 [v(O,t) - v(-L,t)]dt , [13] where e(z,t) is a function of pressure head and can thus be calculated. The increase in water in a profile (first integral form in Eq. [13]) was estimated 14 by applying the trapezoidal rule. The second integral is simply a subtraction of the fluxes at the top and bottom of the profile under consideration and integrated over atj. hi ,j If MBj was larger than a specific value c, the values of and v i , j were rejected,, the time step decreased, and new values for hi 1 Very small values for MB, e.g. 0.1e, resulted in and vi, were calculated. an increase of the time step. A relative mass balance was also calculated after each time step. relative mass balance (percent) The is defined as MB. x 100 divided by the total amount of water that enters (or leaves) the profile across the bottom and top boundary over a given time step, while the overall relative mass balance (percent) is defined similarly but over the total time of simulation. The relative mass balance may be a better tool to check accuracy, especially if the second term of the right hand side of Eq. [13] is small relative to MB. CALCULATION OF THE SPACE STEP Space steps between the various gridpoints are nonuniform and may be changed during the simulation. following way. The space derivative _ihis approximated exactly if h were a piecewise az linear function of z. A function F(z) is piecewise linear with respect to the grid z i if F(z) is continous for -L<:z <0 < The space steps Azi are estimated in the and linear over the interval z 0 The results of this example problem were compared with results obtained from a simulation model published by Van Genuchten (1978), finite element scheme. Table 1 shows that no substantial differences occurred between the two methods. However, there is an important distinction between the two, as is When theta is plotted against depth, the described model, who used a Hermitian shown in figure 5. WAFLOW, concentrates its gridpoints in regions where the largest changes uses a constant space step occur, whereas UNSATI (Van Genuchten's model) between all gridpoints. 19 2. Steady State Infiltration into a Sandy Profile. In'this example (Haverkamp et al., 1977) water is allowed to infiltrate into a 70-cm deep homogeneous sandy soil profile having the following hydraulic properties: K K A A+Ihi K =34 cm/hr. 6 A = 1.175 * 10 S= 4.74 and +Ih+eres + r r = 0.287 " = 0.075 . = 1.611 * 106 S=-3.96 The initial and boundary conditions for infiltration of water in sand were: h = -61.5cm v = -13.69 cm/hr. h = -61.5 cm 0 :5Ls70cm z=0 z=-70 cm t=O t .0 t >0 , where v denotes a constant downward flux and the bottom boundary condition a constant pressure head. Figure 6 shows the water content profiles after 0.2 and 0.4 hr. simulation time (using 71 gridpoints) and also the redistribution of the gridpoints as the wetting front moves down. The smallest space steps The for are at the wetting front where the largest values of 32h/az2 occur. results compared favorably with those published by Haverkamp et al. (1977) the same soil characteristics. 20 Table 2 reflects the effect of the number of space steps on computer time and accuracy. It is interesting that the accuracy, expressed in absolute and relative-mass balance, is markedly reduced only when the number of space steps decreases below 20. However, the required computer time is greatly reduced. 3. Evaporation and Drainage in a Uniform Soil Profile with Fixed Water Table As was shown by Dane and Mathis (1981), evaporation can also be simulated with the model. Their example problem was published by Klute and Heermann (1978), who simulated evaporation and drainage in uniform, coarse, uranium mill tailings with a fixed water table at 3 m depth. The hydraulic properties of the tailings were represented by the following empirical functions: = e(h) =e(h)cosh B - T osh' + -7320~h 0 cm (h)a e(h) = ( h < -7320 cm and eoer where e, e , ho , b, and d are parameters which determine the function of the r initial drainage curve. functions were given: The following values for the parameters of the 21 o = 0.433 d -0.0171 e r= 0.076 h = -131 b =-0.468 a = -0.1964 A = 0.225*108 = 54.29 The simulation of drainage with evaporation was subjected to the conditions: h 1-l cm t=0 t>O t>0 -300z=0 z=-300 cm zs 0 h > -15000 cm q = 1/24 cm/hr h = 0 cm At the surface, a constant flux condition was changed to a constant pressure head condition when the pressure head reached a value of -15,000 cm water. Figure 7 shows an evaporation profile. Again, the data points not only The indicate the water content, but also the space step distribution. obtained results were in close agreement with those published by Klute and Heermann. DISCUSSION The aspect of the number of gridpoints required to simulate a particular flow problem has not been considered in this communication, but it will largely depend on the type of problem and the depth of the soil profile. shown in table 2, reducing the number of gridpoints will greatly reduce computing time. It might also be possible to let the model itself determine This can probably be done by assigning a Values of G above that criterion would As the number of gridpoints required. criterion value to G in Eq. [14]. then result in an increase of gridpoints and vice versa. Values of pressure heads and fluxes could then be found by piecewise linear interpolation between the old gridpoint depths. Table I shows a relatively large computing time for the model presented in this report. However, by increasing the mass balance criterion e, which 22 was originally set to 0.001 cm, considerably. the amount of computing time decreased Different values for e, with their respective computing times and mass balance values, are listed in table 3 for the first example problem. As expected, the accuracy of the solution decreased with larger values for e, Furthermore, the convergence criterion requires relatively small time steps, as the differences in values for the flux and pressure head at two consecutive time steps must be small. The simulation model presented in this report is applicable only to uniform soils. However, it should be possible to extend the model to flow problems in layered soils. The uniqueness of the model is its ability to redistribute the space steps whenever this is necessary, i.e., when large changes in pressure head gradients occur. It, therefore, eliminates manual redistribution of gridpoints and makes computer use more efficient. 23 LITERATURE CITED (1) Dane, J. H. 1980. Comparison of Field and Laboratory Determined Soil Sci. Soc. Am. J. 44:228-231. Hydraulic Conductivity Values. (2) Dane, J. H., and F. H. Mathis. 1981. An Adaptive Finite Difference Soil Sci. Soc. Am. Scheme for the One-dimensional Water Flow Equation. J. 45:1048-1054. (3) De Boor, C. York, N. Y. 1978. A Practical Guide to Splines. Springer-Verlag, New (4) Haverkamp, R., M. Vauclin, J. Touma, P. J. Wierenga, and G. Vachaud. 1977. A Comparison of Numerical Simulation Models for One-dimensional Soil Sci. Soc. Am. J. 41:285-294. Infiltration. (5) Keller, H. B. 1974. Accurate Difference Methods for Nonlinear Twopoint Siam J. Numer. Anal. 11:305-320. Boundary Value Problems. (6) Klute, A., and D. F. Heermann. Tailings Profiles. 1978. Water Movement in Uranium Mill U. S. Environmental Nev. 89114. Technical note ORP LV-78-8. Las Vegas, Protection Agency, P. 0. Box 15027. (7) Van Genuchten, M.Th. 1978. Numerical Solutions of the One-dimensional Research Report NR 78-WR-09. Saturated-unsaturated Flow Equation. Princeton Univ., Princeton, N. J. 24 APPENDIX I APPROXIMATION OF SPACE STEPS For a choice of Az , the following procedure was used. i First note that the difference equations approximate the space derivatives exactly if h and v (Eq. [8) are piecewise linear functions of z. For t and N fixed, let h(z,tj) be the exact (unknown) solution to the boundary value problem given in Eq. [71]. We will assume that h(z,t.) has a continuous second partial derivative with respect to z. Then the zi 's are selected such that, if F(z) is the piecewise linear function with F(zi)=h(zi. j), for i=0,1,...,N, the error, maxjF(z)-h(z,t.) I, -L 5 z 5 0, is as small as possible. That is, the 's we wish to use are those which will uniquely determine the best possible piecewise interpolation of the pressure head at t using N+1 gridpoints. is accomplished if 2h(z,t)2 amax a 2Azi This zi-1 s z S z i is constant for each i=1,2,3,....,N. Although this is, in general, a difficult problem in nonlinear programming, we may obtain a reasonable approximation of the zi ' ssimilarly as in the algorithm NEWNOT (page 184, de Boor, 1978). Let g(z) and = 1 ,z Ng(s)ds = N zo i=1 at-1 zig(s)ds . 25 Since zi.f' z M g(s)ds is approximately equal to Smax zi..1 5s :z z we may select the g(z) AZ i so that zi- 1 r Iig(s)ds /(g(s)ds o To accomplish this, define z G(z) z0 g(s)ds and assume that g(s)$O. G' Since Then G is monotone increasing and so has an inverse G(zN) N = zo zNg(s)ds = N.zi/zig(s)ds 1 , G(z) */zTg(s)ds 1 =-f-G(zN) N) U Hence: Si=1,,...,N. =3 If a is zi = The following example illustrates the selection of z., Suppose h(z)z 0.z 51, and N=2 so that z0=O, z2=1, and zl=a is to be determined. fixed and a piecewise linear function F(z) is constructed, it will look like 26 the dashed line (Appendix figure 1). Approximation theory tells us that the error of.approximating h(z) with F(z) is less than max Osza h"(z) 8 2 = 6a a 2 8Z. a in the first interval and less than max. a8zs1 h"(6 8 2 2 6 8 1 2 in the second interval. Hence, the overall error will be less than max (1-a) The best choice of a is the one that makes 6a max -- '-6 (-) 2 as small as-possible. are equal. That is We claim that this is accomplished if the two errors 6a = (1-a)2 or a = 0.56984. To see that this is true, let a be the solution of ga = (1-a) and 27 max -a6 -a )(-a -(1 )2 Now if a max 6a 6 (1a)2 8· = (l-a.) 2 and again the error is larger. using a.! The smallest error is therefore obtained by In general, however, we do not know h"(z) exactly, but suppose we have approximations, for example h"(z) =1.5 h"(z) =4.5 Now: 0 < z . 0.5 0.5 < z 5 1. fh"(zl) 8 A2 Z 1 Ih"(z2) 8 " AZ2 2 or !nh ti rI Az, =11 hr, -z 2 Therefore, the problem of finding a* may be approximated by finding a*, Assuming a >0.5 then: so that the two areas in Appendix figure 1 are equal. 28 A +(a~o.5)"~~ 1 = Fii..5 and A2 (i.a5 = Solving Al =A2 will yield a= 0.6057 as an approximation to a*,o 29 APPENDIX II EXECUTION OF WAFLOW Appendix figure 2 shows the flow chart of the simulation model WAFLOW. The program consists of two subroutines (UNSTE and REZEE) and four functions (FK, FC, FTH, and UIN), INPUT1.DATA. The input data are read from the data file The functions FK, FC, and FTH define the hydraulic functions K(h), C(h), and s(h), respectively, for the soil type under consideration. The function UIN provides the initial conditions, expressed in pressure heads as a.function of depth. Upon execution of the program, a listing is printed The subroutine of the soil's hydraulic properties and the initial conditions. UNSTE allows changes in the top and bottom boundary conditions with respect to time. To solve the linear system of differential equations, an IMSL Depending on the desired accuracy this external LEQT2B uses an iterative improvement subroutine was used. subroutine is either LEQT1B or LEQT2B. to obtain a more accurate solution. If the solution does not satisfy the criterion e of the mass balance equation, new space steps are generated in the subroutine REZEE. time step is then decreased and the iteration process is repeated. Also the To save computer time, the maximum number of iterations in one time step may not exceed a predefined value IT. The simulation, on the other hand, proceeds in time if both the convergence and mass balance criteria are met. The time step size will increasefor subsequent calculations if the mass balance is less than 1/10 of the imposed criterion. The simulation proceeds in time until either the maximum simulation time is reached (TIM >TEND), or a maximum number of time steps is executed (KOUNT >NT). 30 Appendix table 1 gives a list of the most significant variables used in WAFLOW.. Instructions for preparing the input data, together with a listing of actual input data, are given in Appendix table 2. These input data refer to example problem 2, the infiltration profile. A description and listing of the output of this example are given in Appendix table 3, while the listing of the program itself is given in Appendix table 4. -ow Table 1. Comparison of the Des"cribed Model*WAIL0W with UNSAT1ifo'r a 1.40-rn eepDrai rig of 29 Gridpoints Simulated time__(hr.,) 14.4-2 WF[OW UNSAT1 1.348 1*387 0.51 WAFLOW- UNSATI Drainage rate 5.01 26.44 WA-FLOW 0.*157. INSAT 0.160 ..W.M.M.do- 50,60 W1AFLOW 0.065 UNSATI 00066 (cm/hr. ). 3*346 343-51 0.366 0.379 Overall -absolute mass balance (cm.),, 0.0004 090002 04029 0.037 29s 566 0.055 23.*081 0.033 23.050 0.003 '0.0005 0.08 17*.776 2.41 0*005 Moisture in profile (cm).. 40.o285 Total CPU 40*287 29.543 20.*208 20.159 17,711 1,37 La) E-j Table 2. Effect of the Number of Space Steps on Computer Time and Mass Balance during 0.8 hr. of Infiltration in a 70-cm Deep Profile (Mass Balance Criterion E is 0.001) 3aSS Number of space steps NZ CPU-time (min.) Relative mass balance (pct.) at time (hr.) 04 0.8 1.14 0.5 0.82 1.14 1.18 Overall absolute mass balance (cm) at t 08 hr. 0.13 0.15 Pverall relative aIss balance (pct.) at t 0.8 r. 1.21 1.33 100......... 4.72 70 ......... 3.31 60 ......... 2.86 50 .g...... 1.53 1.53 0.98 1.61 1-27 1.50 0.85 1.42 1.72 0.14 0.14 1.27 ).i A 2.37 1.29 1.20 1.19 1.09 40 ......... 2.08 30 ......... 1.64 20 ......... 1.15 15 ......... -- 1.50 1.54 0.23 0.13 0.13 0.12 - 1.12 1.04 0.91 1.21 error criterion not met WA 33 Table 3. Effect of Mass Balance Criterion on Computing Time and Mass Balance 65 hr. after Orainage Started . ... ... . . . .. ' ' !!l . 1 . I. . .h . " ,,i, .,,J. .... (cm) CM CPU (min.) _NNW" Overall absolute mass balance (cm) at t = 65 hr. -. . .. Overall relative mass balance (pet.) at at t = 65 hr. 0.34 0.78 1.07 .... . ,... 0.001 ..... 2.40 0.005 0. 1.44 ..... 1.15 0.0855 0.1957 0.2703 .. -..;. - 0.01 -: ----.r-.:.--.- . :-. 34 -. · iI (time) 0 I i-I 1 2 M-I M i (depth) N FIG. 1. Diagram showing the finite difference grid superimposed on the depth-time region of a soil profile. 35 ... .0 j1 j (time) M- I 1J'A i-I' (- - i(depth) i=N x 0 I initial and boundary conditions known values for flux and pressure head unknown values for flux and pressure head estimated values for flux and pressure head FIG. 2. Depth-time region under consideration about the general gridpoint (i j) showing the identification of the .gridpoi-nt values of pressure head and flux. a2 0 0 AZ 1 0 0 -1 0 0 0 0 00 0 0 0 0 0 0 .0 0 0 0 0 0 0 O0 O0 0 0 0 O~ -K AZ K V2(0) -.KAz 4, *aH 1 0 0 0 0 0 0 0 0 0 0 -K -ICIR AZ2 K 1 -jC OR 10 0 0 0 0 0 0 0 Sz H2(N-). 0 0 0 0 0 5.1 -K IAZ l 0 0 K 0 1-IC *R -1 0 0 0 0 0 0 0 0 0 o o -K IA ZN K JA ZN V2 (N--1) V2(N BAD C* AZ N 1 112 (N) 0 0 0 11 02 Fig. 3. Matrix Form for the System of Linear Equations as Presented in Eq. [12]. The coefficient matrix A has (2-N+i2) rows and (2N+2) -columns. all other symbols are explained in the text. B.C. refers to boundary condition, while (A) C~ 37 -2.50 FIRST OR SECONO D.ER. 0.00 2450 5.00 1"Wom ... . I . 1400 7.-50 THETA 0 .10 0 0 0 0 0 20 00 (. 0 0 40s 0 0 50 o TH ETA .0 -FIRST DE.R. xSECOND DER 60, 70 FIG. 4. Water Content, z. and as a Function of Depth During infiltration (41 gridpoirits).- THETA 0.00 0 0. 10 0 0 0 0.15 0.20 "MOM 0.10 Wr 0.15. .'.." _ i _ : -. . : ' 0.20 . i''' : I PA ,1p, ,, ,. 20 .h 0 0 0 0. 0 0 O 0 0 0 0 40 '0 0 O 0 0 0 0 0 0 o 0 0 0 0 2 Ia. wJ 60 801- UNSAT. I 9 0 0 WAFLOW 0 0 o 100 .-0 0 0 0 0 0 O 0 O 0 120 1 0 0 o 0 140 FIG. 5. 0 Water content profiles calculated with the models UNSATI and WAFLOW. Number of gridpoints: 29. Time from start of drainage = 14.4 hr. The individual data points indicate the space in the z-direction. 39 0.00 a . 0.05 0.10 THETA 0.15 -- O.20 --- 0.25 - " . 0 0 0 0 O " 0 0 0 0 0 T= O.2019 hrs. SO 0 0 0 30 s CL 0 Ia 0 40 T =0.4019 hrs. 50 0 0 60 0 0 70 FIG. 6. 0 Water content profiles during infiltration, 0.2019 and 0.4019 hr. since start of simulation. The individual data points indicate the spacing in the z-direction. 40 THETA 0.00 0. I0 0.20 0.30 0 0 0 0.40 0.50 500 0 0 0 0 0 0 0 a. 250- 0 0 0 0 300 o 0 Appendix Table 1. Definition of the Main Program Variables in WAFLOW. If the variable represents an array, the maximum dimension of that array is specified. i 1 .1.0 1 Variable Definition Coefficient matrix of Ax=B. system. Input matrix for solution of linear A(310,5) Al (310,7) ALPHI and ALPH2 Working space for IMSL subroutine LEQT2B. Specification whether the top boundary condition is a pressure head or a flux: ALPH1 0 and ALPH2 = 1 + -constant flux. constant pressure head. ALPH1 - 1 and ALPH2 = 0 BETHI and BETH2 Specification whether the bottom boundary condition is a pressure head or a flux: BETHI = 0 and BETH2 = 1 BETH2 = 0 + constant flux. constant pressure head. BETH2 = 1 and CTNEWS DEN - Increase in stored moisture in one time step (cm). Substitute for U2 and V2, used in calculations to check if U2 and V2 meet the criterion for convergence. DT Time step (hr.) Space step (cm) Criterion for the mass balance (cm). Function, which computes water capacity from pressure head value DZ EPS FC (cm-1). FK Function, which computes hydraulic conductivity from pressure head value (cm/hr.). 11 1 LYTIIY-~ 111 III~11~1 Il~-C1 Appendix Table 1. (continued) M I Variable FTH IA IER IJOB Definition Function, which computes theta from pressure head value. Row dimension of A. Error message from IMSL-subroutine. Optional parameter in argument list of IMSL-subroutine. Maximum iterations allowed for one time step. Parameter which keeps track of the times when output-has to be printed. IT ITIM IFLAG Parameter which is initially set equal to 0, indicating that REZEE must be called in the first time step. IFLOW Parameter defining the main direction of flow, i.e. infiltration or evaporation: IFLOW = 0 IFLOW = 1 + Evaporation model. Infiltration model. ISWIT Parameter initially set to 0, indicating that REZEE must be called during the time step when the top boundary condition is changed from constant flux to a constant pressure head. KK KOUNT LEQT1(2)B Number of iterations executed in current time step. Number of time steps executed. IMSL-subroutine. Number of columns in RH. NC NL NSOIL NT Maximum number of space steps allowed. Number of lower codiagonals in coefficient matrix A. Parameter which refers to the soil type. Total number of time steps allowed in one simulation. -r I 1 43 Appendix Table 1. (continued) Variable NTOUT NZ NZPI NZP4 R REZEE RH(310) Definition Number of times that output must be printed. Number of space steps. Number of gridpoints, i.e. NZ+1. Last row of coefficient matrix A, i.e. 2*(NZ-1) + 4. Ratio of DZ and DT. Subroutine that calculates new space steps between gridpoints. Input-output matrix of Ax=B. On input RH contains B, on output X, the unknowns of the linear system. RMAXK SAVDEN SAVEI MAX(X1, X2, RMAXK). Summation of SAVE1 since start of simulation (cm). Increase in stored water over current time step (cm) as calculated from 2 consecutive water content profiles. SAVE2 Intake of moisture by profile over current time step as calculated from fluxes at the boundaries (cm). SAVMAS TEND TH(150) TIM TM(10) Summation of SAVE2 since start of simulation (cm). End of simulation (hr.). Theta, volumetric water content. Time since start of simulation (hr.). Array containing the times that output is desired. Arrays containing previous, stored, and current values of UO(150), U!(150), pressure head for all gridpoints. U2(150) I\ ' I~~C rIIyy Il yI IIIy IIVYIY C~I 11i311 1 IY 44 Appendix Table 1. (continued) .: , :.., .. . t : . . . : : " . . 1 . . . . " . ... . ; Variable UBOT UIN Definition Boundary condition at bottom of profile. Function, containing the initial pressure head values as a function of depth. UMAX UNSTE Convergence criterion for iterative solution process (0.0005). Subroutine which may account for transient top or bottom boundary condi tions. US(150) Array containing the pressure heads at the first iteration of the current'time step. UTOP Boundary condition at the top of profile under consideration. Arrays containing previous, stored, and current flux values for all gridpoints. vo(150), Vi (150), V2(150) VS(150) Array containing the flux values at the first iteration of the current time step. WIN X1 = SAVE2 Absolute value of maximum difference between flux values of first and last iteration in current time step (cm/hr.). X2 Absolute value of maximum difference between pressure head values of first and last iteration in current time step (cm). XL(2100) XMB XMB1 Work space for IMSL-subroutine. Absolute mass balance at current time step (cm). Relative mass balance at current time step (%). Array containing the depths of all gridpoints. Depth of profile (cm). ~~IT~ ~I 45 Appendix Table 2. Required Input Data and a Listing of Actual Input Data for Example Problem 2 (Infiltration Profile) 1.Input- data file*: INPUTi .DATA VARIABLE DESCRIPTION FORMAT COLUMN IFLOW 0: evaporation at top boundary. =1: i"nfiltratio-n or zero flux -at top boun-dary. Number -of space step's. NZ 15 Maximum number of time steps. NT NTOUT Number of times output is printed. 16-20 NSOIL Soil type. is 21-25 11W10 F1095 21-30 31-40 41-50 F1O95 F1O.5 ZBOT Depth of profile under consideration (cm). -UTOP Value for top boundary condition: Negative pressure head-negative-value (cm). Positive press-ure head positive value (cm). Downward flux - negative value (cm/hr.). Upward flux - positive 'value (cm/hr.). UBOT Value for bottom boundary condition. OT Initial time step (hr.)., TEND Maximum simulation time (hr.). 1005 6.-10 11-15 16-20 F5.2 F5. 2 F5*2 F5.,2 ALPHI ALPH2 BETHI1 BETH2 See table 4. 1-50 10F52 TM(10) Array containing the times (hr.) that output must be printed. 46 Appendix Table 2. (continued) 2. Initial conditions. The initial conditions, used for the example run, are listed in the function UIN (Z,NSOIL). Only pressure head values can be assigned to the initial conditions. So I properties. Analytical expressions for e(h), K(h),and C(h) are defined in the functions FTH (U, NSOIL), FK (U, NSOIL), and FC (U, NSOIL), respectively. Transient boundary conditions can be defined in the subroutine UNSTE (UTOP, UBOT, Ul, Vl, NZP1). 3. 4. 47 INITIALIZATTONS AND 5OUNCARY CONCITION-S OF INFI-MCUEL NR. OF SPFACESTEPS NR* OP TIYESSTE?-S ...... ...... 35 600 TIMES OUTPUT IS PRINTEC 000 SOIL TYPE 2' ALPHA2 -o.e0. BETHA20.oO 0.0.*. 0 -046.0 *...- 70.00000o ooeo* DEPTH Of PPCF'ILE (CM) TOP BOUNOARY CONDITI-ON 6.06V mmt3.6900O AL PHA 1=0. 0 BOTTOM BOUNDARY CO3NDITION * -61.50000 6ESTHA1= 1 *0 INITIAL TIMESTE.P (HOUR-S) ..... .. 0.00100 MODEL STOPS AT 6.30000 OUTPUT IS -PRINTED AT: 1.000 0 .00400.eo00.0 Note: 00 0.0 0.0 00 HCUR S Printout was- sp ecified for 1.00 hours. The remiing 9 options -were not use 40 THE PHYSICAL. PROPERTIES -OF SCILN.R PRESSURE HEAO cm THETA 0.2 86 Ce270 0 *22.2 0.164 0.124 0.102 00091 0.084 010081 0*.079 0*078 0.076 0 .076 0.0O76 0.076 0.*076 0.075 0.075 0o075 0*075 0*075 D0*75 0.0O75 0.075' 0*075 0 *075 0.075 0.8075 0.075. 0*075 0.o 075 0.075 0*075 0.075 0*075 0*075 0*075 0.75 0oO75 0*075 09075 0*075 0.075 HYDRAUJLIC CCOUCT'IVITY C M/1HR 0.3248 1E+02. 0.w 15 112E +02 0.35635E+01 0.98843 E+00 0a3-4-98 78+0-0 0. 14832E+000 4715.88 E-b0 1 0.2 1784E-01 0. 1.3224E-01RO 0. 84 18OE-02 0.a557.3 5Et--02. 0*38140-E-0-2 0.26 844E-02 0*.193's56E--02 0.142.55E-M02 0.1069 5E-02 0. 815 6 SE-w03 0.*6312.8E-0O3 0 .495C3.Em03 0*.7243 8E-m04 18 0.* 525E-04 0*64329E-05 0 2' 271.0 78-05 0. 13 054E-m'05 0.6 9322E-06 0.396-6 58-06 0924072E-0.6 0.153 228-0O6 .1.01.44E-06 0.6541 1E-07 0 *488518--07 0.352258-0O7 0.e3 1288 -08 0 *63480O-0 0. 192888-0O9 0.e 7450 88E-10 0.337c53E8-10 0.*17 129 E-- 10 0.94.639E-11l 0.m55-86 lE-l1 00 10'l6OE-l1 0.22355-11. 0. 15211E-11. 0,9 105628-1.1 0.a 7502 73e- 12 WATERC.APAC I TY CM '1 0.v 4 6 9 93Ew03 0. 30 5 E-C2 123 0.e 3931 9Ew0 2 0.e 5118 5E-C2 0*2988 1E--92 0.e 15.5 85E C2 0*81842E-C 0.b4487 7Em03't 0 2 5877 E-C3 0 .1564 8E-'020 0. 98712E-oC4 0.64603E-o-04 0.43655'E- C4 0.3*0333E-C04 0 *2,159 5E-w0 4 0.a157*0 7 Em- 040. 116-4 .2E-oC 4 O.87769-E-w5 o * 67172 E-0 5 O *5211.3 EwpCi o .:69,88-6E-m 06 0*167 8 ZE-oC6 0..548 9E-C7 0. 22464E-07 0. 10'.58E-07C 0.17829 E-m 0 0 -a11111.3E-Oe 0.*72174-E-CS 0.*48 5 26 E-09 C 0.a 335951E-09D 0.a23 8 63.20E--01 0. *1893 8E--10 0.w 3 56190E-11 0.1026IE-ll 0 .37925E-1l2 0.*163006 IE- 12 0.* 814 3 6 E- 1.3 0. 43774E--m1.3 0.2 52 13E13 0. 15347E-13 0. 97739E-14 -30*.0 -40.a0 -60.*0 qm790 ml 30 * 0 i19.0. *0 -300 -l600.0 -o700 -800.0 -I90l00 -1000.0 -1100.*0 -1200.0 o-13000a -14000 -1000 -120000 -35400.0m-4 500.0*! -7500. 0 -83500. 0 -105500 a 0 -1500.0C-127500.90 -1350000 -14500.0 -1.5500. C Note: Above tables were actually generated by the hydraulic relationships specified in the program. VALUES FO0' THETAPRSSO2 HEAD AN TOPBCUN:O*CCN*ITtOCN BO;TTCP4BGUNDCCCNCITIO'N GRIOPOINT 3 4 OEPTH-(CM) 0.00 1.3 .6 0 C 1 FLUX FOR THE SUCCSStIVE <1500 GCRIDP01NTS AT TIM RETA 0.99 8531 -01o 0.9958 5 -12E-01 H (C.14 I -O.-61500E+02 -0.6 15002+02 -0.61.5002+.02 -0.6.15.002E+0(2 -0.,6 1500E+032 -0.615;002+02 FLUX( CM/KR) -0 * 132002+00 -. mC13200 2+00 -. 0. 13200E+00 0 -0. 132002+-00 -C 1320.02e+0 -C4 132002E+00 -0C .132002E+00 -00*13200E2+00C 5 6 7 ;-0.400002+01. 0 19 q a851 E-&0 1 -0. +600002*01+C 0.9585 asIE-Ol0 -0.O 00.00GE+0 1 ,8 O * 998 51 E-001 -0. 140002E+02 -0.*160302+02 -0.v6 15002E+02 -0.V6 1500-E+Q2. -0. 15 002,+02 *6 -0.v6150 02E+-02 -0.*6 15 00 E+.CZ -0.*61.5002E+02 -0.C6 150 02E+02 0 * 9 98 512E-01-I 0.9'9851-0w. 0.99851E-01O -0.O200002E+02 10 -0.G 000 E+ 02 22 0*~ 13200E+(0 C. 13 2002f+00 -C 13200.E+00 O* -0. 132002E+00 13 16 19 o0o.240002+02 G0. 2 80002-+02. -0.30000E+02 -0.O-340002.+02 -0.o360002+02 -0.* 8000 E+42 3 0.,998512-01. -0.6 15009.,+,C2 m0.6 1500-E+02 -0.615 0 0 r202 0o 99g8512-01l 0.99851.2-001 019,98 51 E-01 0.99851E-01l 0.9985.2,-01.l 0. 95985 12- 0 1 00 958512,-01' 0.958512-01 0o 9'98 512E-0 1 0.*998512E--01 -0.*615002+02 -0.6 15002+02 -0.*6 1500.E+0.2 -=0.6153002+0.2 -0.6 1500..+02 -w0.615002+02 -0.-6 150 0 E40 2 -0.*615002+02 -0. 13 20 025+0 0 -0. 13200E+00 -13200E+00 -. 20 02E+00 13 20 21 23 24 25 26 27 -0.o40002+0C2 -0.*4,20002+02 -0.*440002E+02 -C * 13200E+00 -0.*13.2002+00 -0 * 132002+00 -C. 132002+00 -C.* 13 200E+ 00C -0C .132002+00 -. 1l32002+00 -0.* 13200E+C0 -0. 1l3200E+00 -C. 132002+00 '.o-o13 2002E+0:0 -0.13200E+00c -0. o13200.E+00 0. *132002+00 -0.9460002E+02 -0O.480 002I+02 -0.o500002-+02 -0.615002+-02 -0.*6 15002E+02 -0.615002'+02 -004615002+02 .- 61500E+02 0~. -0.6 15002+02 -0. 15002E+02 v6 -0.*61500 2+02 -0. 6 15002+02 -0.O6 15002+.02 -0.61.5002+02Z -0 *6 15002E+02 28 29 0.520:010E+02 -0.~ 5400 0E0-2 -0.*5600.02+02 -W0580002+02 -0 *60CCO +0C2 0.98512-01 0.99851E-01l 0.o 958 5 1E-01 0. 98531 E-*0 1 0.998-51 E-01.O 30 32 -*0.62000E+.0*2 -0. ,64000E+02 S3 34 35 36 -m0.660002+02 -0.O680002E+0 2 -0. 700002+02 WATER IN PROFILE AT TIME 0.0 Note: : 6.98956 CM 'Above data were -generated by assigning values to the parameters that specify the initial and boundary conditions. 50 Appendix Table 3. Description and Listing of Output The following variables are printed for the initial time (TIM=O) and for the selected times [TM(ITIM)1. TIM KOUNT UTOP UBOT SAVDEN - Time since start of simulation (hr.). - Time step number. - Top boundary condition - pressure head or flux. - Bottom boundary condition - pressure head or flux. - Increase in stored moisture over current time step as calculated from two consecutive water content profiles (cm). SAVMAS - Moisture intake by profile over current timestep as calculated from fluxes at the boundaries (cm). TM(ITIM+1) - Next time that output will be printed (hr.). Z TM U2 V2 - Depths of all gridpoints at time TM(ITIM). - Theta's of all gridpoints at time TM(ITIM). - Pressure heads of all gridpoints at time TM(ITIM). - Fluxes of all gridpoints at time TM(ITIM). The following variables are listed at all time steps: TIM DT XMB XMB1 RMAXK - Time since start of simulation (hr.). - Time step size (hr.). - Absolute mass balance in current time step (cm). - Relative mass balance in current time step (pct.). - The maximum difference in pressure head and flux .over all gridpoints, which occurred between the first and last iteration of the current time step (cm or cm/hr.). IT - Number of iterations which were needed to solve the linear system in current time step. TI ME T.1ME S TF.P HIR S NRS 0. 10,000E.W02 0.0 0.0 0.-500-008-03 0.0 0 2. 0 O0E-74Th 0.*0 C.l250E&03 0.00012 0* 1250Oe-0G 0.00025 C-.2500E~-03 0.00 037 0. 15000E-PO3 0.000O6 2 C*2 5000OF.-&3G 0.*00087 0.25QQ-OE-03 0.00112 Cm0.250 OfE-"3 Ov.0137 0 *500,E-o0a 0.00187, 0.5000OE6-03 0.0023-7 0.500005-03 0.00 287 O 0-OOOSoE-03 0.*003,37 G0.50000:8-03. 00.00387 0.10.000E8-0,2 00048 7 O'sIO0QOOE-0 2 0.0-0687 0. 100,0 : 0.2 00.0787 0. 1 0000-*0.2 0.00887 C *10000CE-w0 2. 0.*00987 0.100008E-0-2 0o.01087 0.o100 00r-0G2. O001187 0Q.l10000E- 0Z2 0.01287 0.10000E--0.2 0*01387 0. 10 0.0-0E-0 2 0.01-487 0.10000E-0C2. 0.0158*7 0.100008*-0O2 0.01687 a0.100008--02 0.01787 .0.e100.008I-0 2 001887 0o 00 8O-02 0.01G9.87 0.o1000 06--02. 0.*02087 0 a1 000El!'02 0.002187 0.e100008O-02 0.02287 0.100008-02 0.0*2387 0.100008E-02 0 w02'4'87 0. 1000 08-0 2 0. 025 87- 0.o100008-02 0.02687 00.10000-0w2 0.02787 0.000-02 0.02q87 0*.200008-02 0403187 G0.200008-,-02 0003387 0.-a20000E-02- aBs *-MASS 8AL. C-m 0*65587]E-PCI R ELa .MASSS5AL. 0.E R065.04 0.468083 2P+0 2 0. 8116458-02,j 0.o 62680 E-C3 0. 18388-C3 0o.146 958-0.3 0O.79749-0-*4 0.o291580 0.37T9098E-03 0 i. 24 11 S.-w3 0. 11287r:-03 0.7 480.48E-CS 0-.73971E+~02 0.-a478,948+01l 0.7-410 8O.+01 .*54-6 33E+ 0 1 0. 440668E+0G1 .242078.01 0*43769E8OQ1 00-5591.9E+01 0. 35572E+01 0.1.6650E-01 0*.l134E-00 0.s551608E+01~ 0. 413818+01 0.421128+01l .371-818+01 0.*33*8248E+01 0.3 L373E*+01 0. 244 18 8+01 0.22952E+01 0.212578+01 0. *203.0 l38.01 0. 184818+01l 041 72238+01 0. 161788+01 0.151248+01. 0. 13704E+01 0. 127538+01 0. 11734E+01 0.ollQ04E+01 0.,927058+00 0.859708+0-0 0*700058+00 0.5945E+01 0.249508+01 0 *226 90E+01 0.2 1649CE+01 ER RCRCCNTRCL CM C M/HR c. 43166 7 E+01I 0.a434-93.-01I C*.3C8543E+O0I 0.17 831 -8+C1I 0.44 4448E+0 C. 747058+00 0 * 43755068E+0.. 0.74-79 2E C 3 0. 6425CEIO3 0. 5710-0-m3 0. 50421E-0.3 0.458708-0C3 0.42545 8-03 0.o 3 87328- 3 0. 36.0408-*0.3 0.33.1148-0Q3 0.3!2882 8C3 027*5348-0C3 0 e 2506JE-03 0 *233 518-C3 0 .2 1940-0 03 0.-205 108-.03 0.o 185858.-03 0a.17'295--w03 00 159 1i8-0C3 0. 125728-03 0. 1165,9E-03 0. 10195-0w3 0. 9493 5E-C4 0. 703708-03 04'6 766 58E-03 0.*64457 8-mC3 0.-615288-0C3 0.4 131.1E+00. 74 7 3 84E8+00 4.6 7138E+00 Ceo 3507 88E+00a 0.5 '237'98+00 0. 493-878+00 C.627608+00 0.*60-7548 +00. 0C.*57 .7398+0 .5403.08+00 0.w5 14 106 +;0 0.48870E+00 .466.658+00 0.e446258+00 0.423948+00 06*41223E+-00 0,*318 20GE+0 0 0.3 84298+00 0.372698+00 0.3W60858+00 0 .I51128E+00 0.340928+00 0.*332 6.58E+0 0 0424098+00 Oo 0.309378+00 0 e 3028 5E-CC 0.296528+00 Co29061E+00 0.285188+00 0. 507758+00 0.o490918+00 .*478578+00 0. 46 5438 +00 0.a45-434F+00 10 a 5 5' 1.2 6 5 6 7 6 6 5 6 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5' 5, 5 5 0*05787 0.20000-C2w 0.05987 0.200008E-0(2 0.06187- 0.200008-.02 6.028 -03 0. 4655 88E-0..3 0-o466388-0C3 0.4 0.1709 LE+0l 0.171752+01 0. 17 204E8+0 1 C. 369948+00 0* 306576E+C0 C*1*6078E"OC 0.-0638.7 0.ZOO.0OOE-02 0.606587 020000E-0O2 0.067870v0.000Eaft 02 0.06987 0.2 0 000E--02 00Z 0.07,137 0 *2O0 0.07387 0.200002-M02, 0 0.07587be 2 0-0O02 2 0.2000-0f" 2 0*07787 0.-07 9 87 0*2ZO0O0-OZ* 0,908S187 G0.200008ft-02 0.08387 0.*200008E-002 0.0-8587 C*0.OOOE;-0G2 0 4 4580.2E-03 0 *4416'72E-3 0 *3920 A52,-C3 0.41689 42E+0 1 0.o L6 291 E+0I 0 14460E+01 13 0.. 5378+0 t 0.1187701 0. 9106 14E+01I 0.78 1438E+00 0.03256 28e+-00 5 5 S 0o.3670 28 -0.3 322 2-0C3 0 18 0.21.878 E-C3 0 111* 25456-03 00 145862-03C 0. 11.852-03 3.343202E+0 0 Cl.o. *463:5E+00 0 * 3Z4268+00 5 5 5 5 0.669978+00 0.3o-37522.+00 0.437 08 E+oo 0.O39705S+-00 0.-32931E+100 0. 31 34-CC0 0.33 1308+00 03-0 5-4662+0 5 5 5 0.-08787 0.200008C-02~o 2 O008'987 0 .2O00O:o0 "c 0.09187 C 4O0 00OE-E 02 002-m0 2 0~.09187' 0a.200 0 0.*0 938 7'0400 0.*09517'C2008o 0.09787 0*.20000E8-02 0.09987 0.*200D0-0.2 0.10187 0.20008-02 0*10387 0.2 00004-0O2 0.105-87 0.2000080-02O 0010787 C.20000-0m2 0*10987 0.o2 000.0 -02, Tel I18 7 C 0.00 002E-02 0*11387 0.20:0002---02 0*11587 C0.200002-02. 0.11787 0.920000e.-02 0.119-87 i0.20000EIQ-02 0.1 2187. o00O0OOE-0O2 0.*12387 0.200002E--02. 0 1'2587 0*.200 00 *-02 5 5 5 5 0.102766-03 0. 270170e8-03 0 o .2972-C3 0. 317.5 3-0C3 0.32685Eq'03 0.226922+01l 0. 105352+01 0. 107972+01 0.*11 1308E+0O1 0.*3016938E+ 00 0.3C9172a+OC 0 2,30548E+0 0 .3-7548-+0:0 0.32903E41+00 0259058E+00 0.289588+00 C0.872E+0 0 6 5 00'11543E+01 0.117852+01 0. 9120558E+01 00121.098+01 0.e 126 0 48-+0I 0.12-695E+01 0.*12173.E+01. 0o.122912E+01 0. 12142E+01 0. 1 1580E2+0 1 0. 114358+01 0.10O882+01 0.104162+01 0.*9 95 84E2+.00 0.9473-52+00 0.*879702E+00 0 *80966E+00 0.756402'+00 0.681748+00 0.*6 0870E-+00 0.a5 13032E+00 0.448452+00 0.3 8264E+QO 0.2 72+5 2+ 0 0.189'332E+01L 0 .84741E+00 0.m8g8492+00 5 5 5 5 5 5 5 0.3417 2,-&C3 0.3441 -0-C3 0*33003E--C3 0.* 3 2.5E--C3 33 1004E-*C 3 0*1.3 5 .5 5 5 5 5 .5 5 5 0.29706E+00 0. 29603E+00 0O~.26403+00 0.289142-+00 0.2.87172+00 0.280678+00 0.:2840623+00 0.284012-03 0.0919E 03 2684-C.*3 0.* 0. 2 368508-C3 0 2 19 5 1E-O3 0.a 205072-w03. 0.w 18 4838-0 3 0.e 16530.3-10 3 0 . .1.390S52-003 0. 121582-0C3 0. 1037 4E-C3 0.7 *39025--04 0. 1026g8-02 0. 22'97 IE-C3 0n*24364 2ma- 0*12787 C*2000 0E402 0.12987 0.1-3187 0*13387 0*13587 0.13787 0.13987 0. 14187 0.14387 0.200008E-02 0.200002-.02 0.200008-02 C0.00008-02 0.20000E-0.2 C0.00002-02 C. 200 0 0E-02, C.2 000CE-w0 2 0.a14587T0 .20000E2-02 0.*14787 O0.400002'-022 0.147*87 -0.200002-02 0.1498q7 -0?.20000202 0.-*2+83116E+00 0.2800362+00 0.2799388+00 0*27535E+00 09274616+00 C0o273482-+0 0 0.a2 7 2762E+00 0.5043.72+00 0026937E+00 0426883E+00 .5 5 5 5 5 4 4 4 6 4 4 0.17387 C.200008--02 0*17587 0.200008-02 0.17787 C.20000E-02 0.26441.-03 0.o 25383.3 E-C 3 0.2471 IE-C3 0.975 122+00 0.9527028+00 0.*9,11,32E+00 0.4260982+*00 4 4 4 C.*259139E+00 0.9 258t324E+C0 0.*17987 0,200Q2f-02 0.183870.0022 0*1.85870200-2 0.1.8787 C0.200002-02 0.18987 Cw000E&0-2 0*1.9187 0.002*OE-02 0.v1913 87 Q0.2002D -0a2 0.,19587 0i.,200002-02 04.19787 0 *2'000OEft0 2 0.19987 0.400002-02 0.20387 0.*400002-w0 2 0o.0787 0.400020 0.21187 0.*400002-02) 0.215-870.C400002E-0O2 0.*21987* 0.4010002o-0-2 0*22387 0.*400002E-02 0o.22 78.7 0.*40000O2G-0 2 0.e235222E-03 0. 229'54f-*C3 1.2- mC3S 0. 212 0t.05 52a-03 0.1860,5.5-0oC3 0. 1692 12 0. 153522"C,3 0.o 13'02 942E- 03 0. *114-7 I-,-C3 0.*9 1.5232-0-C4 .G4771.2-03 .a,. 9 8 09 72E-3 0.86746EW+0-0 0. -a78334E+00 0.*687982+00 0.s566 17S+00 o * 2 093.2E+0CC C025-941.2+00 2 C.* 5- 59S 0 0 .25547892.+00G OC4 C*.2-56*7 1E +0O 0.255782+00. 4 4 4 0. 98589-0*3 0.e 9895-6e-C3 0948042-0.3 0.0 839 052E-0-3 0.*80 21 62-0C3 0.754862-0n3 00 484442-03 0 *.608622 0 3 0.512572-0C3 0.38 7072E-0:3 0.539422-0C4 0. 7549S2-0C3 0.77471-0-m3 0. 827792-0w3 0.*869182-o03 0.9 07362-0C3 0,02894-E-03 0.941432-03 0. 9-3696E-C3 0.924142-03 0.903792-03 0.867312-03 0.83 13 1E- 03 . 78544E2-03 0*7430',6E-03 0.69350 E-C3 0. 640182E--03 0.33954 2-C3 0.423052E+00 0.337532.+00 *174762+ 0 1 0. 0.175872a E+0 1 -+0 0.o 18012 1 0.s179 292E+0 1476E+01 04 0 * 168.32E+0 1 0.183722e+01 0.o1,954692E+0 1 0.14789E+01 0 o139171E+01 0. 126 182E+0.1 0.1 1220E+01 0.9449 1.2+00 0.7 13532+00 0.416132+00 0 .994302-01.I 0.3 1497E+01 0.139192+01 0.0142852+01 0.a1,52642E*0 1 0.* 16022E+01I 0o.167312+01 0a171292E+01 0.*173592E+01 0. 172772+01l 0.170402+01l 0.166652+01 0.159942+01 0. 153292*01 001.44832;+ 01 0.137022+01 0. 12788E2*01 0.1.18052+01 0.108622+01 00*25493E2+00 0.*23187C40022 0.400002,-02 0.400002-02 0'40000-0.2 0.400002-02 0.400002-m-02 04400-00E-W02 0*25587 0.25987 :0.400 E-0 2 0.26387 0.400002E-02 0.26787* 0.400002,-02 0.27187 Cs 8000 E-0G2 0.2'1187 0.*400002E-02 0.o27T587' 0 .0000-02 0.27987 0.e400 002-0 2 0*28381 .0.400002-0a2 0.28787 O0.400002-0Q2 0.o29187 0.400002.-02 0.29587 0.40000-0O2 0.2-9987 0.400002G-0O2 0.3-0387 0.s400002--02 0.-30787 0.400002E-02 0.31187 0.400002-02 0.w31.5 87 0.400002E-02 0.31987 C*470000E-02 0*32387 0.400002-0O2 0.32787 0s400002E-02 0.33187 C.4G000E-020.23-587 0.23987 0.243 87 0.24787 0.25187 2 0.* 55 142+00 0.475032+00 0.47215E+00 0.947 106E+00 0.* 46'96 52e+0 0 0.467322+00 O .466762+0.0 0.46 5992+00 0.464302E+00 0.459412+400 0.4597*52+00 0o-.5 9712E+0 0 0.4567.22+00 0.456402E+00 0.*4 55782S+0 0 0.45-5012+00 0* 45417E+C0 0.8 13322E+0 0 C. 4530232E+0 0 0.447882+00 0* 447672+00 0*.445 26E + 0O 0.444992+00 0*444242+00 0 *4437 82E+00 0.o444652+00 0.444142+00 0.43938E+00 0.,442492+00 0.443902+00' 0.44081.2+0.0 0.442872+00 C .44.20 12+00 0.44 1972+00 ..44 Cr 16 621:+00 4 4 4 4 4 4 4 5 5 5 S 5 S 5 5 5 g C C C C C C A-03358Q7 i0.400002n-02- 0.37987 0.400002--02 0*.383.87 0.*400002-02 0938787 0.400002-02 0.8025 52-03 0.82 1252-0C3 0.841.022-03 0.147992+01 0.1.51452+01 0. 15509 2+01 0.4222262E+00 0. 433 182+00 0 .4316CO+C.0 54 0*.391'87 0.400008-02 0.39587 0*.400008-02 0,3998-7 0.40000]E-02 0.40387.0*.400008E.-02 0.407'8710o.400006-02 0.*41187 0.400008-02 0.41587 C.40000E-02 0.41987 0.400008-02 0.42387 O0.400008-02 0.42787 C.40000-02 0.43187 0.400008-0O2 0.43587 10-a40000E-02 0.43 98$1 C.400008--02 0.443 87 0.4 0000OE-ft02 0.44787 C.40000-I8-02 0.45187 C*40000E-0Z 0*4558 7 C.400008E-02 0-45987 C.40000E-02 0.4638,7 0.o4000--;02 0.46787 0.800008-.02 0*46781 0.400008-02 0*47187 C0.400008E-02 0.*475:87 C,40000&*02. 0.47987 0.400008-02 0.48387 0.400008-02 40008t-02 0488771 0.'49187 0*.400008-"02 0*49587 C940000E-02 0.49987 C.4000-00 0.503-87 0.a4:040 0 E-0 2 0 0.50787 C.*4000 -- 2 0.l51187 0.400008.-.02 0.515.87 0.400008E-02 0.51987 0.*400008-02 0.52387 0.400008-02 0.52787 C.400008E-02 0.53187 0.400008E-02 0.53587 0.400008E-02 0.539-87 0.40000-E-02 0*54387 0.*400008-02 0.547 87'0.40000E-02' 0.55187 0.400008-02 0.55587 0940000E-02 0.55987 0.800008-02 0.559-87 0.-*4 0E-02 000 0.56387 G0.400008,-02 0.56787 0.400008-02 0.57187 .0.40000E-02 0.57587 0.400008-02' 0*57987 0.40000E-02 0.*58387 0.40000-0W2 0.58'787 0.*4000-02 0.59187 C0.40000E-02 0.59587 0*.400008-02 0.'59987 0.400008-02 0.60381 0.400008-02 0.60787 0.400008-02 0.61187C0.400OOE-02 0.61587 0.400008-02 0 849408E-C3 0.* 5580 -E-C3 8 0~. 853 0 3 8-03 843798-.03 0.o 0. 83 157 E-*C3 0.o8 16 17e-03 0.e79113E-03. 0.o764 58--03 0. 737 16E-C3 0.71076 8-C3 0.681768E-03 0.6426.20E-03 0. 6 0433 -C3 0. 54285E-03 0.,47502E--030.38132E-0:3 0.25252E-:3 0. 1L0297E--03 0.o 98765E-04 .32558E-02 0. 665198-03 0 * 707548- C3 0.74533E-0C3 0.7863*.38-03 0. 8a985E-C.3 0.* 83E-03 864 0. 8866S88-C3 0. 89505E-03 0. 88584E-C3 0.87166-E-03 0.8505 3E-C3 0.8259 18-03 797668-0C3 0.s 0. 7740 CE-C3 0.74.1428-03 0. 7050 CE-03 0.65241.4E03 0. 586218-03 0.50324E-03 0.3959 28-03 04.24035E-03 0. 103258-03 0. 815098-04 0.3 17508-C2 0.75164E-03 0.76240E-C3 0.-7860 4E-03 0. 810308--C3 0.8308 6E-C3 0.83852E-C3 0.84198E-03 0.83756E-C3 0*8 143 88-C3 0.7429 7E-03 0.6 7884E-0 3 0.5908 3E-03 0.463.34E-C3 0 *15664C+01 00 a43 2908E+00 0*1 5718+01 0.-153'.348-+01 0.150508+01 0. 145888+01 0.,140998+01 0.135938+01 0.1131078+01 0.*125 72E+01I 0.s1-18 51E8+01 0.111458.01L 0.410012]E+01 0.m876 108+00 0.703.3 1+00 0.465748E+00 0.189918+00 0.132 14E+0.0 0.*3001I88-+01 0o012.268.+0 1 0*130468+01 0.a1374 18+01I 0.14496E+01 0 *152988+01l 0.159438+01 0.163458+01 0.a165 008E+0 1 00163308+01 0.160698+01 0 .156808+01 0*152278+01 C.1470 78+01" 0.142718+01 0.136718+01 *0.129998+01 04111030E+01 0. £08 118+01 0.927966+00 C.730078+00 0 .4800 7E+0 0 0.190468+0.0 0. 0150298e+00 0 .292158-+01 0.138608+01 0.136858-+01 0.140618+01 0.144998+01 0.149 508+01 0.153338+01 0.154788+01 0.155458+01 0.15468E+01 0.159.0448+01 0.145 54E+01 0.137368+01l 0.91255 78+01 0.1093-68+01 0.858308C+00 0.43 210E+00 0.443 3178+00 431698.0 0.4 0.4331938+00 0*.433'898+00 0.433 738+00 0.*43 8E+0 0 342 0.43346t+00 0.43415E+0C 0,* 43435E+00 Co.4 33908+0 0 0.434118-+00 0.434768.+00 0.43467E+00 0.43454E+00 0.v27 868+0.0 4 0.-427388+00 0.*42 8986+0 0 0.427988+00 0.v42856E+00 .28 0C* 138E+O0 0.428098+00 0.a42 85 68e+0 0 0.43 0628+00 0.430188+00 004259 18+00 0.*427788+00 0.410538+00 0.42 8878 +00 Go.430 37E8+ 00 0.4320388+00 0.429578+00 0.43 094E+00 .o43 1038+00 C.43 1208+00 0.43 1598+00 0.v432-098+00 0.77877E+OC 0.*428-56E+00 0,,42824.r+00 0. 426008 'O0 0.42644E+00 C. 425 798+00 0.o42 5338+00 0.425958+.00 0.m424658+00 0.426578+00 0. 425048 +00 0. 427 678+00 0.426478+00 0.424818+00 0. 42 792 8+00 0.426818+00 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 7 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5. S 5 7 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 55 0.623$7 0.0 OO0 E0Z 0.*62787 C*.80O08,OZ 0.63187 0.63587 0.w631987 0.,o643807 0*663787 0*65187 0.65587 0.w67947 0.968S7 0.o66787 0.67187 0.67587 0.67987 0.48387 0.67387 0.69187 0.69587 C.400-O08'-02 0.400008-O2 0.400038E-02 40.400008D6 -2 0.40000-QZ 0.40-000-4m32. 0.400008OJ-002 0.400008-02 0.*4000002o 0.4,000068-02 0.400008-0(2 0.400008E--02 (09.400008I-C2 0.400008-02 O0.40000E.-02 0.400008E-0,2 0.400008I-02 0.*30607E-04 Qv.475 32E-Q3 Q0. 7Z61 Cv-3 0. 78-2)18-m3 0.s568 16E-01O 0.250166+21 00.878828+00o 0*.1-33 608*0 1 0.13937801O 0o.144396+0 1 0.148688+01 0o 146468+0 1 0*.140978+101 0.132478+01 C*412 1618+.~ C. 1033-0 1 0.932308E+00 a0.7700.Q5 00 00.54622E+00 0.399971E+00 0.234338E+00 0.964288E-01 0.296668+01 0.135528.01l 0.1287 18+01 0.e13-4468E+0 1 0.*14.111840 1 0.147298.01 0. 157 098E+0 1 0.15864e+01 0013 14r.+01 0. 147528+01 0 * 39458+01 C. 128308.01 0.*114 758E+01 0.*95346E.00 0.7 1624E+00 0.3943 18+00 0.323348-01L 0.439*108+01 0.263628E+0.1 0. 17662F+0 1 0.655918+00 0. 63798-+00 09.6 86168E+00 0.711708+00 0.739'868+00 0 *760058+00 0.77171E+.00 0.775 168.00 0 *777408+00 0.*769638+00 0.74698+00 0.722218+00 C*691238+00 0.6005 8+00 0*.28788+00 Cs.77 220OE+00C 0**4242648+00 42.7E'+:CQ S* 0.422278+00 0.423208.00C 0.422258.,00 0.422328+0'0 C.4 246.+00 4 0 4 21908+0C 04 42Z3340 0.424148+00, 0 a 4238.00+ 0. 424398+00 0o4423400C 0.76554C-00~ 0.41783e+.C0 40.4195386#00 0.4 19a98+CC C.420-038+00 4.42.1198+00 0.420288+00 0.42268E+00 0 4 4244 1 E+Q 0 0e.421-73E+00 0.424178+00 0442,5778+0 0 0.423348+00 0.4 26 25E8+00 0.421538+00 0444'0666E+120 C0*771115E.-a0 G.424258+00 0.22417E+CC 0.225028+00O 00.2 25348+00 0.22-5308+00 0.2 2474840 0.22 5098+00 0. 2248 18+00 0.224-68-CC 0.224828+00 C*22442.E00 0.224858E+00 0.224918++00 0.24678+00 0.22 548 E+OC A.'2563E+0 k0 24L3800 C.222548+00 0 *796478-03 0. 7665,5-E-03 0. 720268-03 0.s661618.SIWO3 00.59 1108-F-03' 0.o506.828-0**C3 0*418606-03 0 .1868-*03 0.217448-03 0.* 12740E-03 0.46998-0-4' 0. 3Z>222-02 0.7 361.28E-3 0. 765988E-03 0 *79,966E-03 .008 28598 0.8467 7E-C3 0.-*85S33-6 E-03 0.85106E-C 3 0.*8"32208G-w03 o0. 18 3*-C3 0,*758918-603 0.66976 1 E-3 0.o240 0.3. 6 0.518508-42C3 0.a3 894'W8-03 0 *46932E-02 0. 1411S8-.-02 0.*47335E-03 0. 1735!-'*C3 0.018748E-03 0. 194518-03 0. 19936E-(:3 0. 20185E8-C3 0.202118-03 0.201948-03 0. 19229E-03 0.*184858-C3 0. 173878-03 0. 166358--03 0. 1589 18-03 0.*149 218E-23 0.69987 C 4t00-0 002. 0m.72938 7' 0400008S-0 2 0.70387 0.40000e-"02 04.730787 0.400008-02 0.731187 0.400008-02 0. 7158'7 0*40000E-02 0.s71987* C.4'0000-02 0.723787 0.400-00-02 0*727187 0.40000;-0(2* .731587 0.400008-0C2 0.7587 G0.40000E-02 0.3987 0.400006-02 0.74387 0.4000E-02' 07787 Go 400008-02!0 0*77187 C.40000E-02 Os075587.020 008E-02 0.759387 0.400008-02 0.76387 (0*20000E-02 0.7787 0.4008C-02 0.*7*7987 O0.800008-02 0*77187 O0.400008.-02 .77387 0o78587 0.77787 0*77987. 0.78187 0.78387 0.78587 0*7787 O0.00008-02 0.*200-00E-02 0.200008c-02 0.200008-02 0.200008-02 0.200008E-02 0.200008E-02 0.200008S-02 6 .7 5 5 5 S 5 5 5 5 5 5 5 5 5 5 5 5 5 5 7 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 .5 5 5 5 4 4 4 4 4 4 4 4 4 0.793.87- 0Co2000E0.2 0.80187 0.200008S-02 0.80387 0.200008-02 4 4 4 4 4 4 4 4 56 0*80987 0 *200002-02 00.2-00-002-02 0.8 1387 Cis 200-02E-G 0*.a.1587 C.2000a2---02 0.81787 0*81987 C. 00002O-60:2 0982187 0.200002-02 0.82587 0.20000-02 0.200002-0O2 0.*82 787 0.,82986' 0.20000E-02 0983186 0*83386 0.*83586 0.200002-02 0.6-837186 0.w 00006-02 2 0.8'3986 C .200002I-02 0.a84186 C.200002-02. 0.e84386 0.84586 0.8a4786 0.8-4986. C. 200:00-02 0085186 Coo0.00aem02 0*85386 0.85 786 0.859,86 0*86186 0.o863.8'6 0. 865 86 0.*86986 0.87386 0487386 0.87586 0. 87786 0.87986 0.* 88.18 6 0.883 86 0. 885 86 0.*88 786 0.89186 0.89586 0*89986 0 *9078,6 0.91186 0.01586 0 *91986 0.'92386 0,93186 0.93586 0.93986 0*94786 0*95586 0.9 6386 0*97186 0*97986 C.* 0002j-m02 20 C.200002Eii--02 C. *200002---02. C.*2 00002:-02 0c 200002-02 0.*400002E-02 0940000E-02 0.20000E-.02 C.20 00 02E-02 0. 20002 C .200002E-02 C.200002--02 0.20000E2-02 C.200002-02 C. 20000E-02 C.200002-02 C.e40000E-02 C.400002-02 0.4-00002E-02 0.400002-0.2 C*400 0 02E-02 C .4000-0O2 0*400008-02 0. 400002E-02 0.*400002E-02 0.40000E2-02 C. 40000E-02 09 80000E-02 C *800002E-02 C. 800002--02 C.080000E-02 C.* 800002E-0 2 0. 800002E-02 0. 12182-OO3 0.116 502'-p03 0. 112 18E-03 0.411 200-0-C3 0. 11407t"-03 0. 11510E--k'0 0.1Z&26E-C3 0. 1322CE-C3 0*.1 582-03 0.14465E-03 0.*148972-0C3 0.m 115606E-0 0. 155642,-C3 0016136E-03 0.1611-9E-03 0.16366E-*03' 0.157782-0O3 0.1563 62E-0 3 0.15.284E-G3 0.14869-0C3 0. 1428 5F..03 0.13689--3 0.13408E2-03 0. L24.8gE~-C3 0.*12240E-03 0. 11223E-u^3 0~o.111 252-C 3 0.10242E-C3 0.98299E-C4 06*800675-G3 0.12-030. 282322-C3 0.1763 7E-C03 0. 141812--03 0. 1255*92-03 0. 1145 1E-C3 0o'105'8 1-C3 0. 100802-03 0.954422-04 0 *2542 5EO-C3 0.30736E-03 0.272162-03 0.241542-03 00 .268 8E-0a3 0*19110E-03 0. 163602-03 0. 1528CE-C3 0. 133572-03 0.118092-03 0. 104752-C3 0.*93 1432-C 4 0. 112492E-03 0 *2470 92-03 0.*172162-03 0.145012-03 0. 123191E-03 A0^n9-8786e 0o.06 662E+0 0 0.4gl22E+00 0.w4,80112+00 0. 487342+00 045051-9E+00 0o-.20152+00 0.540 72E +0.0 0 5 5969 3 e+00 0.641142 +00 0.671112*00. 0.-74266e+00 0.*789446+00 0.i 855'45E+00 0. 884312E+00 0.9'52272+00 0.990 15E+*00 0. l04,852E+01 0. "105642.01l 0.010961E+01 0.w 11237E+01l 0.114862+01 . 1.1613E+01I 0.*i17 29E+01I 0. 12124 2*01 0.0119 34E+01 0.123 74E+01I 0.120 172+01 0. 126282+01 0. 1233 52+0 1 0. 12572E+01 0.542432+01 0. 1582178+02 0 61 11 8,9E+02 0 5 19 282+01 0.35601E.01 0.3073 22+01 0.2'90132+01 0. 28133 2+01 0. 276262E+01 0.*2795852+01 0. 28 12 12+01 0.402912I+0O1 0.5376-42+01 0.532632+01 0.528632+01 0.555952+01 0.523842+01 0.4999 12+01 0.519862*01l 0* .05462+0O1 0.496602+01 0 .4 890 12+01 0.482222+01 0.3 18 132+01 0 * 7845 12+01 0 * 75563 2+0 1 0.754692+01 0.745892+01 0.747072+01 C *22345 2+00 0.6224522.+00 0.2248.32+00 0.2 24452E+00 0.2233-82+00 0.*222872+00 0.222632E+00 0.22 1642+00 C.220422+00 0.2 19342E+00 0. 217482+00. 0*21593E+00 0.21372.00 0.211345+00 0.20867E+C 0.*205292+00 0.202 28*2+OC 0.198512+00 0.s194052+00 0.189372+00 0.18486E-+00 0. 417971E-+00 0.174012+.00 0 * 67922+00 0.*16 147E+00 0*.54792+C0 0.1479,92+00 0v.L44+00 0.134-282+00 0.*128012E+0,0 C. 890462+00 0.2 1053 2+00 0. 149622 +000.8-2170E-01 0.697242-0 1 0.646702-01 0. 60 8092-01. 0457330E-01 0c 503728'-0 1 0.478-592-01 0.*789902-01, 0.702142-01 0*624802-Cl 0.55518E-01 0.53672-01 0.447882-Cl 0.105593E-01 0 0 '317332-01 0.283372E-01 0.253262-01 0 *226312-01 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 4 9 7 5 5 4 4 4 4 4 4 5 S 5 5 4 4 4 4 I' 4 4 4 0.3la4 76 82E-01 0~'.2562-010 *250502-0-1 0.2 13 292-01 00*18 155-C-01 0.154142-01 4 4 4 4 4 4 VALUES FGA. THFETAtP!RESSUR-E HEAD ANO FLUX FOR THE SUCCESSIVe TCFBOUNO.aCCNoITICN -3.00 SOTTOtSOUNOvCGAOI T IOGN -.6.1.5sWOO GR I OPIO IN T 3 4 5 6 7 9 12 13 14 16 18 19 21 -2 -23 24 25 2627 26 29 30 33 34 35 36 OEPT-H C4) -. *3519 7eC02 -0.a454-12 C-4C.2 TNET'A 0o.6741E4'0O 0.267 ZZE*+OO 0.266.5 9.E+0Oa 0 026359se*12 O.Z62C1EiOO 0-00 6 02 89E002 -0.*62017f+-02 -0. 640122 0.02 go6 5 4 E +-0Q 0*Z5302E'OO 0.289 4E+000.t344IE+00' GRIO1POtNT$ AT TIME:?.:e 020 72690+02 .0082217396+02 m-0.o2311E0*02-0. 222200+02 -0.o23538+02 -0.w25 650.02 -0.3289e+02 -0. 270570+02 -0.e 20870IE+02 -04905E+02 -0. 51455E+ 02 o3 -0. 36494E+02 -0.*37839E0.02 -0.* 4215 E+02 3 0-.4206 00E+ 02 "*Q43526E0+CZ -. -0045022E+02 -0.O4654604-02 -wc*13690.02 -- * 13664+02 C ,-43 * 13 6491E02 -0., 13622E+02 -C* 136 1IIf +,02 -C. 136016-E 02 -oOw-135926+02 -C.Z 13584E402 -C.4135701 -oC e 13570E+"02 -Ca *13364E+%02 -C *13 558Ee:CZ 0-* 13 5543E+ C2 -. *13546e+,02 -0.135460+02 Z+0 01*.13 5400 2 08Ct 13538.02 Col13533E+02 -0.6 1I33 +0 2. -1.135330+02 -C. 135320+02 -0. 1353.2t+02 4C.13532E+02 -C *13 52E+CZ 1.3 -0I~ 5 3 E+-a2 '-C. 13532E-+02 -0. 135320+02 -C. 13532E+02 -0.67275IE+C2 0-*.67698F42 -0*.6805 0+012 -12. 683490E+02 -"0.6859,90-+2 -0.w68991E+02.C -0.*6945+1-2 -0.69.2750.02 -0.w693880+02 -w0.694S4E0402 -12.6 9638E+02. -0.*69000E*Q2 -0.69?54E+C2 0069801E+02 .0.69,842.E+02 -0.*698770E+012 -0.699090+02 "0.69936E.02 -0*699600+1-2 -0.*69980E+22 -.*7ooooEi'02 0.Z20632E*00. 0. 198-8 1 EeCO 0o-.l9LLS8E+Oa a*I3.c 37593E+ 00 0.1 L6126E+00 0.a 14143-E+00. C. 135559E+00* 0. 13011E.10 0. 12507 E+00 0. 12041E+QO 09'116 13 E+0 0.e 10541 E+00 0. 102 50E+00 04.998-5 1E-01 -0.5292 1E+02 -0. 54583,E+02 -0. 5 6 27 1E+ 0.2 -0.579330+02 -0.o 59730E+02 -0.6 l50OE+02 CHANGE IN STORAGE IN THE WHOLE PRO3FILE TOTAL -FLUX 0IFFPERENCE BETIWEEN TOiP A40 BOTTC?4 TILL NCh COMPUTATIONS WILL PROCEED TILL. C.1l4227Ei-Q2 CM 0.1156470+02 CM 6.o HOR -HGURS 58 C Appendix Table 4. Listing of WAFLOW ONE DIMENSTONA L SIMULATION MOD-EL CW A F LOaw ' C*A C COmMON TFL0WvNSnILTIM 0Imf NSION U01(1501 fWAI 50),U2.l1501,9VOt 50),Vlftt5O) DIMEFNSION US(15-0),-VS( 1501,vTH( 5) Z(1501),T#~i 1). INTEGER IFt /IAAA4'/,J1IEVAP'/J21'TNFI'/ c Cr***THE DATA ARE READ FRIO' DISK ************** RSAOI3,051 ITFLOWNZNTNTr,-UTNSOIlL 25 30 'FORMATI'5T5) FORMAT( 5F10.5 9 It4F5*2 /9,OF5*21 FL=Jl1 TlF(IFtOW.o.O*) 1Ff IFLOW.EQ.I) F-L-J2 C *** 40 * 45 THFSE DATA ARE WRITTEN TO UNIT 6 ************ W*ITE(6,45) FL ,NZNTNTO'jTNSfTtZBOT,UTOP,&LPHt, ALPH2,UPIOTP 1 BETHI, BETH2,O0TTFNDTM FOR7MAT!' INTrIALIZATIDNS AND 13OUNDARY CONDITTINS 0F',?Xf/ 2' MR* ~'A0NQ. ,5 nF SPACESTEPS ..... OF TTM4.STEPIS .IP1...,59/0 41TTIMES OUTPUT IS PRINTED ***vj,5,/, 59 SOIL TYPE.........%5/ 6' DEPTH 'OF PROFIL.E (CM) 0000019FIO951/9 7' TOP BOUNDARY CONITTION .... ,FlO.105, ALPHAt=',F3.1,' A1P1A2=', 11 9VTTOM BOUNIARY CONDTION *.',10.5,' BETHAl=fF3.1.,' INITIAL TIMESTE-P (HOUPS) .. ',FlO.5,f, MOnFEL STOPS AT.......'FO5'H'JS,() 2' OUTPUT IS PRINTE0 AT vf,2Xl0F5*2,' HOURS,3(/) Of 1' BFTHA2=0# 59 c c C c C ! FLnw*o TV LOW* I "Sm.!Lm t NSIOILw=2 ....... ...... .. ....... E&PMT!N4OL NP!L.TIATVrN -14OfEL e i...... LUTh '400E1 C NI fLw3oooo ANO ALPHAZSI ANO AtL04A 2=0* AND0 SI1A2=1 ANO 8ETN&A2*0 C TOP ROUNDARY CONOT!ON: >AL+44tuO CON4STANT FLUX r. ALOH*9&lsl c Ct1NS1TANT VISSSURFE4-A0 C BOTTOMN BOUNDARY CtnN0!T!ON:* CnNSTANT FLUX -)A"H&1-ml c RETH41t c CONSTANT PQESURE HAO c 'r C o C.0- *0ARCY POSI'TIVE FLUX NEGATTVE FLUX - UPWAPV) FLOW* PLOW -mvqo>OWNWARO CSI 10 46 1*1092.00vIO THIET a=IFTHf.UNSOTL) COMO = FK(UNSctL) CA P = FC(uvNSOIL) WITE(6,56) U,,THWFT,CONDCAP COWTNUS O0l 47 130091400t0O -T)4ET 46 = FTH!UNSOIL) CO~nO *FK-(U9NSOILI AP=FC(UN5fL)1 W~tTS(6*56) U9TH4ETCONDC.&P 47 Cn-NTYNUF n0 49!=5055000 U..i .*!~ THET = FT4(UNSOTJL) CVOND CAP =FKfUvNSOIL) 2 FC(UJNIS1J!L) 56 41* WRTTF(6,56)' u#TH$ETv,oNn),CAP FlM"?,~I7,538,1.,XF25 CfNTTNUS 60 Cl C SOMS INTIALIZATICONS SAVMAS =uO SAVOEN = 0.a FPS .001 11250 T'T M = 0. oz -Zr4OT/NZ NC 2150 N7+1 NZPI TFLAG =0 C INITTAL. VALtJFS OF ZIJV AMD TM AT TIME5 ZERO 50 00 60 !=1,NZ'PI ZI!) 2 FLtAT(!-wt)*0Z 60 65 U-011)=UTN( Z t t) NSOTL) Till Y)=FTH'(UO(l I ,NSOIL) CONTIlNUF 00 65 Tvr7,NPt CIN= FK( .5*(U0(-t)-+UO 1-11lNSVIV. VfTI-zCOlN*f((JOE!)-U0(!-1l)/OZ) *CON CONTYINU F MEAD AND FLUX PESP.* C **L!ST T"S INTIAL VALUES OF '0EPTH,TMFTAvPRESS-URE WQ!TE(6 921TIKOUNTtUTOPUEO0T 00 77 Il,NZII U(I WRIT'E(6,751 I,(I T(1 CON'TNIJF CTNEW = 0.0+00 orl 80 1ItNZ CTNEC-WTNEW-(TI(Y CCNT NUE CTNFrWS =2 CTNEWf2.WRTE(6,85) 7114CTNEWS WRTT(6, p-60) V(I T7 +THII+1) )*.t! 80 C STrPF VALUES FEW NP4 =2*(MZ-1)+4 NINZ-l n0 qo J=1,N71ol UlTJV=Uo(J) 'ATRTX SFTUP Vii J)=Vfl(J) A% CONalrV NUSI 61 c C 'T 4ATN LOIP: NMTTS IS'4UNT OF TIME STtPS &LLOWET) TO RI&Fl #ENO# Sn ON:IN * TRA S t4T TOP 0 O TON04S Rr? T TOM R fq1 1ifA R YCONODIT CALL UNST-EIUTOPUiDOTptJ1,vV1,ZP1.) C 1o0, t00 r*,Nzo1 ItO UOtI) u Ul) TH(IO FTHft-lfI),.SOT.L) -V*Itl= VtII) CORTITNUE IF(!LOWEO.)GOTO 12() C C C 14CONSTANT FLUX TO A CnNSTANT PRESSUPE HFAO ROIJNOARY IN"* CHANGE rOM ( -)ALPHAIwI CASIE OF EVAP.'40EL ANO IF H4< -o.5000 CM tVFUOEI).GT.-15000.1 (qfl 120tz 1Ff TSWT.EQwf.) CALL PE7E.r.t ZU2,V~tU2,V2#NtN.C) PEZPE TS ALWAYS CALLF0 FOR THE FPSTr TIMIE STEUP AFTER CHANGE TO) CT7NSA:NT POE$SUP5 HEAO.0 ALPH41 a 1. AL0142 a 0* UTOnP = tOWi C r t2' C C CO)NTINJF. C C 'TT" IS THF' MAXs MUJ'4EP Of ITEPATIONS ALLOWED IN ONE TIME STEP* THE IMPUT MATRIX a IS rrOtpqFSSED) IN BAND ST1ORAGE MCDE. IT HAS 5 COLUMNs AMD t?*NZ +2) ROWS. no0 700 KKIt,tT =0.0+00 0. M,+00 A(194) At 1,51 P14( I) =AIPHI =ALPH2 fl 0+D00 =UTI2P 62 on 130 1=19NZ j =2*1 At Jl)= 0*040'0 4(J3) +.5*O)Z A(J14l A(J'2 A(J#5) %*D RH(J.)3A(.*D A(J*1) A(4,2) fJ,3) A (J,4) 41 A(J,5) = .5*R*FC (.5* (UlII1+UI T+tfINSIL) +to 0.0+00 =AtJ,1,) RHfJ) 3A(Jvl)*(ti0(I)NJ~fT+11) 10 CONTINUE A(N2P4vl) a0.0+00 A( M2P4,2) = SETHI A(N294,9) STH? A(N?04,4) A (-M2P4, 51 RH(NZP4) =0.0+30 0.0+00 UROT M-mI I JOB=.O SUROITTRE LIMT2 c c IS MapE ACCURATE THAN LEOT1B A,!JOB1,Xt.,!ER) CALL L OT1P(AN2-P4,NLNL,!AP.RH'4 CALL LPQT2R-(ANP4,NLNMLIARH,'4,!A,!JO36At,IAXL, IERI ItFflSRm.Q.29-) WRITE(6,*)TIM C r * EP9.1 MESSAGE: ** IEvk=12Q ERPI?0 -. MATR!X IS SINGULAR MATRIX A- IS TOO ILL CONO)ITIONED, ODOES NOT C U2(1) 00 140 RH(1) 1,1NZ =RHf 4.l) UZ(T+11 140 CONTINUE VNP.LVRH(N2P41 **THE FIRST COMPUTATION !F(KK.GTOl)GO2TO 160 00 150 I=107NPI USII) = '2tI) VSIT) =V2(T1 FO.R EACH TIME STEP IS MEM4ORITZEO:l US ANO VIS* I1r 160 C * r * CONTI NUIS NO MATTER HOW LARGE TH15 ERROR TNJ THE 14ASS BALANCE TS, TAKE THE LAST COMPUTED VALUES IF KK SOUALS IT's IF(K*5Q*tTIGOTO 113. C'0-4PUTFl AND) THE RSEVTCU.S -MFWLY C **CHP-CK FOR CCNVlERGE.NC. -T0 THE C * VALUFS DIFFPE0 TOO) 4JC,TR.Y AG;AIN WITH THE, NEW VALUES AS tNPUT* O EION * U4.AX 4UST 8f. SMALLEP THAN Q,,O005 C *~THEC0. UM A X 00 17T=i 2,?4P1t I-EN UZ fI)I -'UfI1))1/OEN.') IIO c UMAX* AM 1 f UMAX, A 3S ((U2(fI CONTIrLNUE .I U4AX.GT...00051 GOV] 600 CONVERGIENGE CRITEIRMN IS 14ET KL -a(K 00 190 !w-,NZPI 1807 WAIT'S VAI.UE IF PRESSURF HEAD ITSPOSITIVE AND CALCULATE THE4F C * C **CORRESPONDI1NG CH-AN-GE IN TH4ETA FROM WATEAR ETriNTION D.ATA,. I ,-U2( I) 119 U2(T ).GS.O.000001 I WRITE.(6tS~v CO3NT INUEF 190 c c * THE FOLLOWING STATF.4ENTS CHEC'K THE MRASS 9ALANCSE EOUATION. IS RELATIVE M4ASS tALANCER-*ELATTVlE Tfl TH FOIFFERENCE C ** XMl~ C.** IN FLUX SETWEEN THE TOP ANI T HE 8OTT11M OF THE PROFILE. c CTN-W 10.0*+00 oil. 200 IT=10? CTNEW 2 CTNEW +T(7TI41((VW!4) 200 CONTINUE CTNfIWS 2CTNFW/?2. SAVEI 2 CTNFW$ , V0( ' 0 V (Nt. W,N -V2(tfV SAVS2 mWIN XMR - ASICTNEWS-WIN) XM8L =(X48fA8S(WIN))*lOO.O c c C THE MAX. PRFSS5UIR HEAD OR FLUX DIFFECRENCF, OCCURINMG OVER THEL WHOLFE C q~rOFtLE BETWEEN THEF.rIRST ANI3 TH4E LAST- ITERATION TN THE CURRENT TI'4 IE C' STIE- IS PRINTE-D AS AN FRROR CONTOOL. RMAXK = 0.0+00 Om 210 I12tNZPI =AtStUIIl-4SU)) 2I 211 X2 =Ai SfV2Z(Tl-VSfl) R.MAXK a'AKAXI(X1X2.RMAXK) CO.NTINUF PMAXK =.5*PMAXK 64 SCUT5 MASS BALANCE * XM8.1 REgLECTS THE * tp 'XmB' IS LESTHAN A PREDEFTNEO VALUP C **COMPLET-R3. C c C c c IN T14E CUJRRNT TXME STEP-* *EFPS' THE :TIMESTEP [S '4XM85> 0'EQSlCAL EZEE. WHT-CI4RECZO41'ITES THE IF H0Wf'VER rv (RIDPOfTTS fSrA-NCElS, AL'SO THE TIEME STEP WILL THEN SBE H.ALVEO. WRQITEF(6,730) TlMQTXl8,K'Ml~RRAXK,$KK REZEE IS ALWAYS CALLED -IN THE FIRST TIREST'eP: IF-(TfFLAG.FEQ.O)XM-0lS1 IFLAG =I IF(XMR.LfwEPS )GOTC) ?50 CALL PEZE5(ZvUtllO,10VONZvNC) 011 240 TIV1,NZPl Ulf T)MUO( 1) V-11 T )=vat I) CONTINUF OT=OT/2.* GoTO 100 P******** EADY FOR THE 240 c c C 250 c C * C * NEXT TI4E STEP ********* TIM a TIM4 + OT CU4UtATIVIE CON4PlNENTS OF THE MASS SAVMA.S AND SAVIEN ARE THE TWOC BALANCE EQUATION .URfNfG THE COM-PUTATIO'NS. SAVV'AS sSAVMAS + SAV'9l SAVOEN SAVO.EN + SAV'E2 PRINT VALUES IF *TI.1' EQUALS ANY 'TM(ITIMJ' !F(TIM.LT.TM(14o.R.TIP*.E0.(NTDUT4I)) GOTO 500 TIP43 TTM+1 CALL RE!EFtZvU2,V2vU2,V2,NZvNC) C C c C C * * 300 C C * COMPUTE MOISTURE CONTENTS FROM WATER.RETENTION DATA 00 .300 1NP TH(T) =FTH(U'?(!)*NSOIL). CONTMNE ALL THE COMPUTED VALUES ARE NOW WRITTEN TO DISK WR!TFf2, '10)Tl.,D0T WRT (2,2q01zI(rlT=lqNZPl) WQITF(2,.2'0) ((T),=,~ WRITrEf2,7Po) 1V2ff) ,131,ZPL) WRITFf6,71") TTM9K0UNTUTOPUBOT c 00 150 1=1.NZPI 65 ** 5T1JESIPf) CVPUTArTCNS ARIF FTN!SHFO ECUAL-S TEND l.THF TO t .TM TV T14f 'ASS 8AI.ACr.-tS MtJCI4 SALLER -THAN 'THq p#s0eFlNc-r) 'cPSOVALJE Nt)THE f AN E C PNW GRlfl "s0!NTS DISTANCES ARIE CCUL t C * TIf4SSTE'P MAY OF LAP(vR FC07PSUeSF.0UJ.NT TIM-E STE:0S. O T!8P20 l'fTlw*(TI?4.GT.TEND) 1540 ftPS ) GOTa 600 T9XM9G C C CAL c s0n C *T1e CONTIN4UE 'CC-IMPUTED~ VAtUSFOP FLUX ANO PREtSSURE H1CA0 ARE NOW RIEST01ED. 00 650 I*INZPt* U1111 =U2(11 VI =V24I) CONTINUE IF(KK.EQ.!IT.OR.U'4AX*.LT..00iO05) GO TO 800 650 c 700 I-On 820 C 1515 CtNT INUE CONTINUEF WR!TS(6,840)1 TEND FtDPMAT(111 THE P14YSICAL PROPEPTIES OF SOILNIP * lt3,04(f),' PPS-SURE HFAD',0X9 2THFTA',!X,lMV%)PAULC CONDOUCTTYVTTY',3X.'WATERCAPACITYt, lFWlQAT(1HI'VAlLUFES FOR THETAP'RESSURF. HIEA' AND FLUX FOR THE'1, ct0*5,* TT45STEPNR. '405o 1.' SUCCESSIVE ORIOPOYNTS AT TIME:-, OTT04BOUNO.CDNDIT ION *B15,, ?(f)tf TOlP8UN0.CO0NOTIOfN tOHC)1XFL(CfR) 45X,*0EP1H(C4)q7%qTHETA FrDR~lAT25v4Xq4(3XElZ.5) E,74,:,tD5'CMI#3(/i) FMR'4AT(f,' WATFR TN PROFYLE AT FfDP'ATfT75,'PI'S. PRESSUPr- IHAD ATGPDIT:,3lZ5 FOOM&T(Imi2xv I' M'F3,4XIIMSTIXA8S~m&SS8AL*' 7? 75 l5 18,5 760 913X, 'PEPC' ,5Y, 'V149,4fHR' ,SX,' TTEQATIONS') -;7X,'HR*S',tt3X, 'C", 2 %.30 FnPMAT($ 27 51 l 04l FOrQ4AT('f3'Cn*PUTATVnms aqF EmnFi3, TTME IS 'tg7.4,1 HRS') NT-FWHOLE PROFILEI I15XE15w6, TRG RSO~fR'Artr/,'CHANGE T iw C-4',*/,tyvlTtDTAt FLUX 0lFrFREFIC'F !rTWEfrN TOP ANC) BOTTOM TTLL', 7' NcWtE15.6q' CM',?(/I,' COMPUTATTflNS WILL 00OCFEE TILL tp lF7*4,' HOURS') STOP END C 66 C SURMtTINE PEzEE(zDtnfVUVNNC? c C c DIMENSION Z(NC) tJ4NC) ,V(NC) ,flV(NC) ,Dt(NC) DIMENSION Ut'EWI 150) ,VREWIf 15*),ZNE.Wt 1501 ,VMNEWH 150) DIMENSION TH(150),CO)EPG(2,300PO),DV1(150),oulUSO0) c THI 1)nFTH(fljU( 1),NSft) TM( KV=FTHIDfUI(K) ,NSO1L) PUpposE: 0TV RECOMPUTE NF'W SPACE STF-PS Ct"4MMON !FLO'INSTLpTIM O.K=(FK(D-UII41)mNtrhLJ-*FK(rnUII),PN.S01')I/DZ VNFWtIUV= CONTINE (DVlfl)-OUIHf)*0K--YK)/0DIVK I r C c c !F(TI.*LT**041 GO0 TC 8 Z---> DEPTH Z VNFW --- > SECON1 0ER!VATIVE OF M WfTH Z 9Ml -> FT*RST 0EQ~[VATIV'E Of H WITH Z OVt --- > FIRST OERIVATYVE OF V WITH Z WRT(1890) TtM WRTEII,?O0) (7IEI i ,)N T= WP!ITEVfI,7-00-) Ot ( I -p I -P N) Tz WR TTF ( t?001 (DV () T,IflN) CONTINU(E DVMYN v .010 TTFUF1OLWEO*t) rOV!ATN=0. CO:F(1,1) =0. 01.110 T=1IVN IF (VNEWC!).LT.oV'41NIVNFW(I)=IVMTNf COSEFG ( 2,91T) zSQRT ( VNF-W ( T CloFIVltT+1) = CEG1!).OF(,)l(I)Z STEP = COFG(l,N+t)/N IF (STEP*LE.0s1 GOTO 90 1J- I C c 8 L0 +) 67 42.1 i) ZN~FWC!I1vwIZtj-)+Z(J~t) )/20 CONTINUE GVT0O10 00 9S vt IN, Z!NEW (T I vZ(I I-wLO AT (TAW I*STI' GITO t127 101. 2, M0 110 TV(ZN4EWfTI.GE.Z(-j4tl) GO+TO-101 c c 127 1111 -L!NIAI1? TNTEMPOtATION NEW~ COM4PUTED GRID POINTS. OT r-R,4YIN0PRO-04 T4C P H AND V ARE E ?fZNEW(1V-Zf.J) )/(zfj~lzfiltJ UNE+W( T)=tJJ)4+(U:(J~1)-lJlJ) i*R CONTI NUE TU(T )UNEW(IJ~lF6,l~)TU( IN RZCE#v'5S20*4) FVVMATftX, Z(T)=*ZMFW(T) VTTlM*LTe..C'4) GO TO3 25n 12' 1211 C 1I> Ig0 210 C 25O TI~f r )wTH(-Uf II ,NSO!L) CONTINUJE WRPTFEtlt00) (ZH(lYhlvLK) 7s. Fn A f -FrQMAT 6F.2o 51 Q-FTUR N .68 c C * c t, 2 C 5 4C C IO PUNCTYnN FTH(HtNSOTL) COM-PUTFS WATF-R CONTIENT FROM4 WATER RETSMNTIO REAL MITF..pI YF(NSO.EC*.31 '00Ttl 10 fF(NS0TL*F~ot GIO TO 5 SOTLN-R. 2 !FfH*.GT.*-*Ol GOTO 2 VMf.T.-7120.) GnlTO I BEv t3t./ASS(Hfl***46A rvAM =f.4-13-..76))ff.43~l+.O76O) FTI4 .443*(COSH(BE)-GAMI)f(,CSS~4EGAM4, 1R FTUP N FTH = (t*tF-2A3SH.*,t6 RIETURN PTH = .413 RETUJRN CUPQVFS. FTH=I.611E+O6*.2217/t'I.&tI'Eo06.AB.S(H)**3*96)..0075 IQET UR N SOTLNqo 3 T.F(HGT*-.OO0t) GO TO t5 4=3 s 57168 TE--f 1.0/f l.O+f&LP*&BS-(H) )**N) )**M 9TM4=O. '98.*TE+0. 069 P ETUR N 15 FTH=1.365 RETURN 69 C * FUNCT1.0'N FKL)'4NSVtt)HY:OR&AJIC CONDUCTIVITY VAtUES FROM4PRESSURE !P(NSOIL.S-0*I) GO TTO20 0TO Ii. O~.~O1) ov TP( TP-(H.LT.-T32.0.) GOT') I HEAO OATA 4 rAM a(*413-*760)/(*433+*076O) T r3T70 10 171s2AR()**l6 T T 2to T - .3 ' CONTTNUE FX=5.4E.0--w0*vEXPf4.2q*T)/24e 15 17E+6fI,@ t7E+06 ASS(H I**4.e74) -FK=4.*to 2t) lF-H.GT,-m.0OOIj) TS=l0/fj(2 TE~fI(Q1-1 QETURN TFt GO TO0 30 P*8StM I*AfSH))*) (A*mvi * 30 70. C ** VALUES FROM4 P~fSShJRE MEAf) DATA* WATER CAPOACT ~AK N.?q4T.1* T~fNSfltL*.FQw3) 170TO 1I) GO TO I lr-(NS0,Tt.El SCtLN~. 2 tffH.GT.-.Ol) ME RETURN RFTIJRN C~ GOVI 2 WI*46 /ASSl.f*13lAt H**. QETUPN C 5 *61o6AS()*A Q C' C SO!LNRo 3 TV(HeGTs-o00O1) GO TO 15 10 Wm*Iw(+1f ALP* A8S tH I t**N I S PC aO. V2,46*4M*HH*N*t ALP** N) (4AR (4)**T I Q ETUR N~ ~C1C0 RESTUR N Fmo 71 c FUNCTION UTf(ZNSf-T L) C **1HP TNtTTAL CONOITtON~qt EX0RfSSEOI)IN C **-PUNCrION OP DEPTH* !P~tN.SCI1Tt.E-o.1) GO Tr 5 ' VPfStL-..0,3)(GO Tfl 10 UTM . -t 5 to Ulm a -61.3 RETUR-N UIN 0-26-000 PETURN ENO~ "PESSURE 4644f) VALUES AS A SUBROUTINE C tNSTFUTOPtUlTUIV1,NZPI) PURPOSE: TO ACCOUNT FOR TRAN'SYEVT TOP OR BOTTOM BOUNDARY CONO[TTONS. CIV4R'ON !FLOWtNSOVITIMOP4E.NSTON Ultl5O0YVt(15OIVp (TTM*EQ.0O) GO TO 20 UB*T AU ( ZPI) -UROT *-FK(U80TAN4SOIL) USOs(+**I-*~I UTOP UTOp RETUR N END c c 20 F~IO **~*A JrSB p751 AYL59JH5 6201 AE 1 t7 PAGES T 8 YN 2 124 Mo O0.l 10*050(it AM 0! ******* 0 i 010 81IN 1 0 72 / I I A2 AZ AZI 05 for selection APPENDIX FIGo lo APPEDIXI~IG 1.Illustrations of z.i. 7 -. 3 APPENDIX FIG. 2. Flowchart of WAFLOW.