Pergamon Com/Wers % Geosciences Vol. 22, No. 10, pp. 1083-1096, 1996 Copyright 0 1996 Elsevier Science Ltd PII: SOO98-30@4(%)00047-7 Printed in Great Britain. All rights rewwd 0098-3004/96 515.00 + 0.00 AN IMPROVED 2-D FINITE-DIFFERENCE CIRCULATION MODEL FOR TIDE- AND WIND-INDUCED FLOWS FERNANDO J. CAVIGLIA’ and WALTER C. DRAGANF Facultad de Ciencias Naturales y Museo, UNLP, Casilla de Correo 45 (1900), La Plata, Argentina, and *Departamento de Oceanografia, Servicio de Hidrografia Naval, Av. Montes de Cka 2124 (1271), Buenos Aires, Argentina (e-mail: dragani@oceanar.mil.ar) (Received 23 August 199.7; revised 20 March 1996) Abstract-An RM/FORTRAN computer program for the numerical solution of tide- and wind-driven currents is presented. The model solves the nonlinear two-dimensional hydrodynamic equations by a finite differencing method on a constructed zigzag boundary. The alternating-direction implicit (ADI) algorithm is implemented, and the model requires only the solution of block-tridiagonal systems of lin- ear algebraic equations. In contrast to the conventional vertically integrated model, the bed shear stress is determined from the bottom currents, which are computed using the linear three-dimensional hydro- dynamic equations and the Galerkin technique in the vertical. The program includes an option to run the model using the conventional representation of the bed shear stress based on the depth mean cur- rents. Applications under idealized and field conditions are described also. Copyright 0 1996 Elsevier Science Ltd Key Words: Numerical hydrodynamic model, Tide- and wind-driven circulation, Alternating-direction implicit algorithm, Galerkin technique, Continuous vertical current profile, FORTRAN. INTRODUCTION The prediction of tide- and wind-driven circulation is of the utmost importance in ocean and coastal engineering. Tide level, tidal currents, and wind-gen- erated currents usually are studied using well- established, two-dimensional, vertically-integrated models. Since these models are vertically integrated, the bed stress is computed from the depth mean current. However, this is not physically correct, and a more realistic parameterization is obtained by using the bottom current. Davies (1988) presented a model which consists in solving a nonlinear two- dimensional, vertically integrated model where the bottom stress is no longer related to the average depth currents but is computed from the bottom currents. This transformation can be performed if one supposes that the current velocity is composed of the depth mean current and its deviation from the depth mean value. Then, the linear three-dimen- sional hydrodynamic equations are split into two systems where only one of them is depth dependent. The application of the Galerkin method to the depth-dependent system permits a continuous verti- cal representation of the magnitude of the current deviation from its depth mean value at the calcu- lation points of the two-dimensional model. The re- lation between the vertically integrated model and the three-dimensional one is given through the bot- tom shear stress, which does not differ significantly in magnitude from that obtained with a full three- dimensional model (Davies, 1988). It should be noted that the model is linear only in the sense of neglecting nonlinear terms involving current struc- ture. This paper describes the computer implemen- tation of a two-dimensional model in which the bot- tom stress is calculated from the bottom currents. The computational requirement of the model is not significantly different from that of the conventional two-dimensional model and is small compared with the requirement of a full three-dimensional model. To run the model efficiently on a desktop computer some simplifications were made which also are described here. HYDRODYNAMIC MODEL Vertically integrated equations Mathematical characterization of tide- and wind- driven currents requires the simultaneous solution of the dynamic equations of motion and the unsteady continuity equations. If it is assumed that CAGE0 22/l@ B 1083 1084 F. J. Caviglia and W. C. Dragani vertical accelerations are negligible, pressure is hydrostatic over the depth, and the fluid density is homogeneous, the two-dimensional depth-averaged nonlinear equations can be written as follows: all a ,+,(h+~)U+-&(h+~)V=O (1) all 1 apa (k - Tbxbx) =-gax-pax+ p(h+q) +&V2U (2) g+ug+vg+/Z! aq 1 apa by - rby) =-gay-pay+ @I+?/) + AHV2V (3) where t is the time, U, V are mean vertically aver- aged velocities in the x and y directions respectively, g the acceleration due to gravity, q the free surface elevation above mean water level, h the depth below mean water level, f the Coriolis parameter, p the water density, P, the atmospheric pressure, 7sx and rSY the x and y components of surface wind stress, rbX and rby the x and y components of bottom shear stress, AH the horizontal eddy viscosity coeffi- cient and V2 the horizontal Laplacian operator. The expressions for the x and y components of surface wind stress are given by rs, = k’aCdaiVlOiVl0, 7s~ = L’aCda~~lO~VlO, (4) where lVlol is the wind velocity at 10 m above water level while Vlo, and VQ are the x and y com- ponents of wind velocity, pa the density of air, and Cd, the wind drag coefficient that in accordance with the WAMDI Group (1988) can be calculated as follows: Cda = ( 0.0012875 IVi01<7.5 m/s 0.0008 + 0.000065lVt0l I VIO( 1 7.5 m/s (5) In the traditional two-dimensional vertically inte- grated model, the bottom stress is determined from the depth mean current using either a linear or quadratic form of bed stress. In a three-dimensional model the bed stress is computed from the bottom currents. In order to enhance the two-dimensional model, the bed stress is determined from the bottom currents using a quadratic law (see Owen, 1980; Davies, 1988; Davies and Jones, 1993): 7b.x = pkub ub + “b 7by = &Vb ub + vb d= d= (6) where r& and vb are the bottom currents in the x- and y directions respectively, and k a dimensionless coefficient of bottom friction. The bottom currents and bed stress are computed without formulating a full three-dimensional model. The method developed by Davies (1988), and Davies and Jones (1993) is used. This method con- siders a point model to determine the vertical cur- rent profiles and hence the bed current and bed stress. Vertical projles of the currents In many tide- and wind-driven circulation pro- blems of practical importance, the vertical velocity component and its gradient are considerably smaller than horizontal velocities and gradients, and can be neglected. Thus, flows induced by wind and tidal forces can be assumed to approximate horizontal flows. Neglecting nonlinear terms and the horizontal shear, the hydrodynamic equation of motion at a point can be written in complex notation as $+i$W=-g(g+i$) where w = u + i v, u and v are the velocities in the x and y directions respectively, d = h + q, A, is the coefficient of vertical eddy viscosity, and 0 = (‘1 - z)/(h + u), a convenient nondimensional vertical coordinate, z being the vertical coordinate measured from the mean water level and positive upwards. Following the development of Davies (1988), and Davies and Jones (1993), the current velocity is written as w = W + w’, where W = U + i V is the depth-averaged current and w’ = u’ + i v’ the cur- rent deviation from its depth mean value. Then from Equation (7) follows $+ifiW=-g($+iz) (8) where T, = 7sx + i r,,j and ?-b = rbX + i Thy. Equation (8) gives the depth mean current while Equation (9) defines the current deviation from its depth mean value and gives no net transport. In order to obtain the current profile and its as- sociated bed stress, Equation (9) is solved using the Galerkin method in the vertical using a finite ele- ment representation. In this regard, the current vel- ocity w’ is expanded in terms of m complex coefficients A,(t) and depth-dependent real, continu- ous functions g,(a) (the basis functions) as follows: An improved 2-D finite-difference circulation model 1085 (10) The coefficients A,(t) are obtained by substituting Equation (10) into Equation (9). In general, there will be a residual error R associated with this substi- tution because the basis functions (10) are not the exact solution of the Equation (9). The Galerkin method consists of multiplying the residual R by the basis functions and integrating over the region q 2 z 2 - h or O< u 11. Integration of Equation (9) yields (11) with k = 1, 2, 3,...m. Integrating the term involving A, by parts gives grgk da + U-T-4 J 1 grgk do r=I 0 where A, was considered independent of the depth coordinate, although still varying with time, and the boundary conditions given by -P Avg _ [ 1 =T, -p A,? Tb (13) O-O [ 1 a0 o=, = were used. It should be noted that A, can depend on the depth coordinate in the general Galerkin method. For the particular example in which A, is independent of the vertical coordinate the basis functions are eigenfunctions of the eddy viscosity profile (for more details of the integration procedure see Davies, 1980a, 1980b). Many authors have used the Galerkin method with a basis set of cosine functions (Davies and Furnes, 1980; Davies and Heaps, 1980; Pearce and Cooper, 1981; Davies, 1986): g, = cos((Y,u) (14) Since the boundary conditions in Equation (13) are included explicitly in Equation (12) it is not necessary to choose u, such that these conditions are satisfied by Equation (14) (Davies, 1980a, 1980b). Then, a suitable choice of CI, is ar = (r - 1)K (15) Due to this choice for cl*, the basis functions (14) produce a mode (r = 1) which is constant through the vertical. However, by definition w’ cannot con- tain a mean flow component which implies that A, must be zero. Taking into account Equations (14) and (15), Equation ( 12) becomes f$= _[if+!!$] A, +$(T, - T&OS&) (16) with r =2, 3,...m. In order to solve Equation (16) it is necessary to specify the form of A,. Hess (1976) has presented an expression for A, as a function of tidal currents. From this result, it is assumed that A, = r(U ’ + V 2)“2(h + q) (17) where r is a dimensionless coefficient of order 0.005-0.0005. The simple eddy viscosity formulation given by Equation (17) appears to be supported by more complex turbulence energy models (Davies, 1991). By solving Equation (16), the current velocity wb’ at the sea bed is determined from m wb’ = W + xA,cosa, r=2 (18) and the bed stress computed from Equation (6). Davies (1988) showed that using three or four terms in Equation (lo), the computed bed stress was equivalent to that determined with a full three- dimensional model. To minimize the computer memory requires the coefficients A, to be held, therefore only three terms were retained in this paper. NUMERICAL MODEL Finite difference equations The numerical scheme describes the current vel- ocities, the undisturbed water depth, and the free surface elevation at different grid points as shown in Figure 1. The alternating-direction implicit (ADI) algorithm was used to solve the partial differential equations (see Douglas and Gunn, 1964). The ADI method is known widely because it avoids the matrix inversion problem related to the implicit schemes. However, the method is not unconditionally stable when it is applied to the complete set of Equations (l)-(3) in an area with real bathymetry and irregular bound- aries. Thus, a limitation on the time-step is necess- ary to ensure numerical stability. Calculations usually become unstable for Courant numbers, C,, greater than approximately 3-5 (Weare, 1979), 1086 F. J. Caviglia and W. C. Dragani where At is the finite-difference time-step, Ax and Ay are the horizontal grid spacings in the x and y directions respectively (Fig. l), and h,,, the maximum water depth of the integrational area. Before replacing the derivatives with an appropri- ate set of finite-difference analogs, the continuity equation (1) is used to eliminate the terms $-$ from the momentum equation (2) for the x component, and $$ from the momentum equation (3) for the y compcnent: au uarl u2ah uvad ---------- at d at d ax d ay 1 a& + (Gx - tbx) + AHV2U --- P ax dh + 4 (20) av vaq v2ah uvad --_--__--_ at d at d ay d ax 1 apa + (% - Thy) + AHV2 v --- P ay Ah + 4 (21) The finite-difference equations corresponding to differential Equations (l), (20), and (21) are expressed in an AD1 form, with the time-step split into two halves. For the first time-step, from time VAt to (n + i)At, the terms involving U and n of Equations (1) and (20) are expressed implicitly in the x direction while the terms involving V are expressed explicitly. Once the values of U and rl are computed, Equation (3) is expressed implicitly in the y direction and V is determined, at the same time-step, using the values of U and ‘1 previously computed. At the next step, from time (n + $At to (n + l)At, the procedure is reversed so that V and 9 are computed from Equations (1) and (21) which are expressed implicitly in the y direction with all terms involving U expressed explicitly. Using the computed values of V and g, the new values of U are determined from Equation (2) expressed im- plicitly in the x direction. i+l i i -1 i -2 J-2 J-l Jil J/2 -Y + water surface elevation, i odd and j even n mean depth, i even and j odd I x direction velocity, i even and j even B y direction velocity, i odd and j odd Figure 1. Computing grid. An improved 2-D finite-difference circulation model 1087 The equations are discretized in such a way that the continuity equation is evaluated at nodes corresponding to r) values and the momentum equations for the x and y components are evalu- ated at nodes corresponding to U and V values, respectively. Consider an M x N matrix, with M and N being the number of grid points in the x and y directions. At the interior nodes (x, y) and denoted by the sub- scripts i, j, where x = (i - 1)Ax and y = (j - l)Ay, the Equations (1) and (20) are expressed implicitly in the x direction using the following notation: f $ = f (iAx, jAy, nAt) .f n+@ = &f y + (1 - p)jynj rx = ;ti+Lj +A&) r= ;u;+,,,, +f;-1j+1 +h-L-1 +h+lJ-l) as rl n+‘/* = nn _ $ [(JY + $)u];+‘“/* -$-[(I;‘+ cy)V]lf at i,j ++q”-& g- [ ( u *y+w* - @Y + (fX)“> 1 s+p’2 + (GX - r;:““, + At A u ,, p(hY + (ij”)“) ix7 H .Xx At + 8Ay2 -&U;yatif l,j (22) The approximation is central in space and the coefficient p takes any value between 0 and 1. A value of p > 4 permits theoretically the use of higher Courant numbers. It should be noted that Equations (22) and (23) form a system of nonlinear algebraic equations. In order to avoid iterative procedures, which are time consuming, and noting that the simulated flow con- stitutes a slow transient phenomena, a linearization procedure has been adopted to reduce these equations to a system of linear algebraic equations from which a direct solution can be obtained. Each term Q involving a nonlinear function of U”+‘/’ and/or $‘+ ‘I2 is linearized by means of a local Taylor expansion, Q n+’ = Q” + B”(U TZ+’ _ U “) +c n(?f+’ - f) + O(At*) (24) where B” = [ $$I” and C n = [PI”. Furthermore, at each time-step, the total water lbepth is expressed as &’ = h + $‘- Ii2 Applying the linearization procedure to Equations (22) and (23) we obtain the following system of linear algebraic equations r);J+r/* + G T,, U ;;,‘,y + G ;,, U ;-:‘;’ = G ;,d (25) At each time-step, Equations (25) and (26) lead to a tridiagonal system of linear algebraic equations which is solved by a forward sweep and a backward substi- tution. The algorithm is simply Gaussian elimination. Once the values of U”+1’2 and r)“+“’ are com- puted throughout the computational domain, the values of V”+li2 are determined from Equation (3) which is expressed implicitly in the y direction, v n+t/* = v n _ & u”+p’/7 vx”,: At At + 8Ax* -z&V;_~ +- AH V n+p/2 at i, j. 8Ay2 J’Y (27) Applying the linearization procedure to each term involving a nonlinear function of V”+“‘, we obtain the following tridiagonal system of linear algebraic equations F n ‘u V ;J’f + F ;,J V ;J+“* + F ;,, V $4’ = F ;, (28) Once the values of V”+“2 are computed, the time-step is incremented from (n + i)At to (n + 1)At 1088 F. J. Caviglia and W. C. Dragani and the procedure is reversed. Equations (l), (21), and (2) are written in discrete form as tl n+l = p/2 - & [(hY + :x)u];+“2 _ & [(p + +jY) ~]~+(‘+1)‘2 at i, j (29) vn+l - v n+l’2 + v n+(r+‘W - (I;x + @Jy+“2) x [(tjy)n+’ [ _ (qJ)n+v2] + e v n+ol+lvy~); *?lt +&u) - n+1’2(hx + py+W)] +&P n+(P+w(fiY)~+l’2 - [(U)(tix)x]n+“Z] +5+-)“+1’2 2 n+bL+lV2 g - (;:+)($)“+“‘) 1 n+(ti+lW ‘y tTsy - T:y+(p+1)‘2) At AH v n+1,2 +&x + (?y)n+l’2) + g&2 xx At +- 8Ay2 AH V zz”2 at i, j + 1 u n+l _ _ un+iP - & [q p)J+(l*+1)‘2 At -- a*~ v”+(~+‘v2(~Y);+“2 (30) + (‘y - T::(p+‘)‘2) + *t p(hY + (ij”)“+“q is AHU;;(lr+,),2 (31) where the following notation has been used f y+‘)‘2 = pf ;:‘j” + (1 - &??-“2 1J At each time-step, the vertically integrated model is coupled with the three-dimensional model through the bottom shear stress. In order to avoid an iterative scheme between the two models, the nu- merical solution of Equation (16) is carried out using the bed stress and eddy viscosity computed at the previous time-step. Thus, the complex coeffi- cients can be determined directly for each individual grid element as follows An+112 r 1 -$(l -/.L)M: +$[T’- T;cos(Y,] = 1 [1++4:] (32) where 44: = [if+$q. Initial and boundary conditions The initial conditions are taken as still-water values, that is U, V, and rl are zero everywhere. However, the model accepts as initial conditions for surface fluctuation and current velocities nonzero values which must be dynamically consistent. Along closed boundaries the normal component of current is set to zero. For tide-induced circulation within the compu- tational domain, it is necessary to specify the tide ele- vation and current velocities along open boundaries. Since, generally, the current velocities are not known accurately enough to prescribe them for the model, the values of U and V on the boundary were set equal to their values at the first interior point. For flows induced by the atmospheric forcing the variables at the open boundaries are calculated using a partially clamped, explicit, radiation con- dition given by where 4 = rl, U, or V, 1 = x or y, and Tf is the friction time scale, and whose value should be about 100 At (Chapman, 1985). The radiation con- dition (33) is applied only to those nodes with depths less than 200 m. The other nodes are ap- proximated by clamped boundary conditions in which the open boundary variables are fixed at zero throughout the computations. The model is spun up by increasing the tidal and atmospheric forcing from zero linearly over one inertial period. Nonlinear aliasing effects The model includes the horizontal eddy viscosity AH which filters out the spurious growth of short waves. However, unreal values of AH sometimes have to be used to avoid the large aliasing error introduced by the occurrence of nonlinear instabil- ities. For these cases the nine point filter of Shapiro (1970) is introduced. It consists of the periodic ap- An improved 2-D finite-difference circulation model 1089 18 17 16 16 14 13 12 11 I@ 1 9 B 7 6 6 4 3 2 1 i i i i i i i i i iT ~--~---L--~--~---L~-*--*--~~~~*~~*~ [ i i i i i i i i i __f__j___LL_i -‘+‘-*---~--~“+-“~--) _-_*__ap i i i ,a_. _+__i___p_$_ __?___~_‘l___~___~__I___~___~__-_i___t___~__~__~___~__~__ 111111111I I II 11 11 11 1 2 3 4 6 6 7 6 9 10 11 12 13 14 16 16 17 18 19 m 21 J Figure 2. Idealized domain with computing grid superimposed. plication of an operator which spatially averages the computed velocities and water surface displace- ment at the end of a fixed number of time-steps by using a predetermined weighting factor as follows Yjdf’ Yjj+;(* -s) x("i+2,J + yi-2,j + \Yi,j+2 + *jj-2 - 4*,ij) +T(oiz.l+Z + \vi+2J+2 + *i+2,j_2 + *;-2j_2 - 4*j,,) (34) where Y’ is the averaged value of the variable Y which stands for U, V, or r~, and s is the weighting factor which can take any value between 0 and 0.5. The adopted procedure has the effect of simulat- ing the transfer of turbulent energy to scales smaller than the double of the grid resolution. It should be noted that the user must specify the number of time-steps at which the calculated values of U, V, and rl are filtered. where the nodes with even values of i and j corre- spond to U, the nodes with odd values of i and j correspond to V, the nodes with odd values of i and even values ofj correspond to n, and the nodes with even values of i and odd values ofj correspond to h (Fig. 1). The input variables and parameters of the program are prescribed in five files. In order to describe these files better, an application of the model to an idealized domain, as shown in Figure 2, is presented where the computational domain and the boundary are superimposed. The first input file is called PARAMETR.DAT and contains 35 data in I3 lines. Figure 3 shows the corresponding file for the example given in Figure 2. The data contained in the PARAMETR.DAT file are explained in sequential form as follows: NIV number of columns to be swept in the x direction. The first three records of the configurational file given by Figure 4 are indicating the initial and final nodes of each column. In the example, from j = 9 to j = 15, six columns must be considered because the island is not included in the computational COMPUTER IMPLEMENTATION domain; Input speciJications We first describe the input specifications and then the various subroutines of the program. The compu- tational domain is contained in an M x N matrix MM: number of rows to be swept in the y direction. The fourth, fifth and sixth records of the configurational file given by Figure 4 are indicating the initial and final nodes of each row. In the example, from i =8 to i 1090 F. J. Caviglia and W. C. Dragani 11 7 21 18 22 I5 11 8 1.0 0.6666666 0.0 500000.0 500000.0 180.0 48 981.0 1.0 2000000.0 0.005 7200.0 0.0050 1 12 43.00 1 0 BOUND.DAT CONDINIT.DAT METEOROL.DAT 5 ETAVEL.DAT 17 27 37 47 NNMMNM NNN MMM NNNN MMMU & & it KEND g P A, k T/ r IBAND IFILT c$ ICORI IMETEO File name 1 File name 2 File name 3 NE File name 4 TX1 TX2 TX3 TX4 = 10, two rows must be considered as a consequence of the island; N: columns of the matrix (odd); M: rows of the matrix (even); NNN: columns of the computational domain. In the example, from j = 10 to j = 14, ten columns must be considered as a consequence of the island; MMM: rows of the computational domain. In the example, at i =9, two rows must be considered as a consequence of the island; NNNN: columns between the columns to be swept in the x direction; MMMM: rows between the rows to be swept in the y direction; cc this parameter takes the value 0 or 1. If c( = 0 the nonlinear terms of the acceleration are not considered; Figure 3. View of the PARAMETR.DAT file for example given in Figure 2. j.~: implication coefficient; y: this parameter takes the value 0 or 1. If y =0 the atmospheric pressure gradients are not considered; Ax: horizontal grid spacing in the x direction, in cm; AJJ: horizontal grid spacing in the y direction, in cm; At: time-step in set; KEND: total simulation, time, in h; g: acceleration due to gravity (cm/set*); p: water density (gr/cm3); AH: horizontal eddy viscosity coefficient (cm*/sec) and whose value should be about AH u U L, where U is a typical velocity and L is horizontal grid spacing (Ax or Ay); k: coefficient of bottom friction; Tr: friction time scale; 6 4 3 3 10 3 10 610 6 8 10 12 14 8 14 8 14 8 16 16 14 4 6 8 lo 10 12 12 14 14 16 18 5 7 9 9 11 13 15 5 3 315 5 7 13 13 17 9 19 19 19 17 5533222292929 2 9 5 9 5 5 5 7 7 11 11 13 13 15 15 15 9 15 9 15 9 15 9 17 9 17 17 17 17 15 15 3 4 5 6 7 8 9 lo lo 11 11 12 12 13 13 14 14 15 16 17 18 19 3456789 9 10 11 12 13 14 15 16 i 6 4 4 2 2 2 214 2 4 4 6 61212 14 14 14 18 18 20 10 20 20 20 20 20 20 18 18 3 3 3 16 16 8 10 12 14 16 -1 -1 -1 1 1 6 4 3 10 3 10 6 10 610 8 10 12 8 14 8 14 a 14 8 16 14 5 7 9 9 11 11 13 13 15 15 17 4 6 8 8 10 10 12 14 7 5 3 15 5 15 7 13 13 13 9 17 9 19 19 17 Figure 4. View of configuration file, for example, given in Figure 2. An improved 2-D finite-difference circulation model 1091 0 5 IO Km I? dimensionless coefficient of proportionality of the vertical eddy viscosity; ZBAND: this parameter takes the value 0 or 1. If ZBAND =0 the velocity vertical profile is not computed and the bed shear stress is calculated using the depth average currents; IFILT: number of time-steps at which the calculated values of U, I’, and 9 are filtered; r#~: mean latitude of modelled region, in degrees. Positive in the Northern hemisphere; ICORI: this parameter takes the value 0 or 1. If ICORI =0 the Coriolis force is not considered; IMETEO: this parameter takes the value 0 or 1. If IMETEO =0 the atmospheric forcing is not considered; File name 1: name of the file that contains the positions i, j of the computational domain boundaries; File name 2: name of the file that contains the M x N matrix with the initial conditions and depths; File name 3: name of the file that contains the M x N matrix with the atmospheric forcing; NE: node numbers of the open boundaries where the U, V, or r] values must be specified as boundary conditions or calculated with the gravity radiation condition. When the computational domain does not have open boundaries NE must be set equal to 1. File name 4: name of the file that contains the U, V, and/or r] values which are used as forcing in the open boundaries for tide-induced circulation; TX, TH, TX3, and TX4: times at which model outputs are desired, in h. In general, the computational domain will have an irregular geometry so that the nodes i, j of the boundary must be given in a configuration file whose name is given in the PARAMETR.DAT file. The file corresponding to the example given in Figure 2 is shown in Figure 4. This configuration file consists of 21 lines which are explained in sequential form as follows: I-initial i of each one of the NN columns to be swept in the x direction; II--final i of each one of the NN columns to be swept in the x direction; III-j of the NN columns given by I and II. For the example given in Figure 2, the first column to be swept is given by j = 4 from i =6 to i = 10 (see also Fig. 4); IV-i of the MM rows given below by V and VI; V-initial j of each one of the MM rows to be swept in the y direction; VI-final j of each one of the MM rows to be swept in the y direction. For the example given in Figure 2, the first row to be swept is given by i =5 from j =5 to j = 13 (see also Figure 5. Resulting surface water elevation in cm and depth of mean currents for example of Figure 2 1092 F. J. Caviglia and W. C. Dragani Fig. 4); VII-i of the previous nodes to the initial nodes of each of the NNN columns of the computational domain; VIII-i of the subsequent nodes to the final nodes of each of the NNN columns of the computational domain; IX-j of the NNN columns given by VII and VIII. For the example given in Figure 2, the first column is given by j =3 from i = 5 to i = 11 (see also Fig. 4); X-i of the MMM rows given below by XI and XII; XI-j of the previous nodes to the initial nodes of each of the MMM rows of the computational domain; XII-j of the subsequent nodes to the final nodes of each of the MMM rows of the computational domain. For the example given in Figure 2, the first row is given by i=3 from j=6 to j= 14 (see also Fig. 4); XIII-i of each one of the NE nodes of the open boundaries. If open boundaries do not exist then i= 1; XIV-j of each one of the NE nodes of the open boundaries. If open boundaries do not exist then j= 1; XV-the nodes of the open boundaries can be placed on the right, left, top and/or bottom margins of the computational domain, and they are recognized using the numbers 2, -2, 1 and -1, respectively. This is necessary for solving the radiation condition. Zero must be used if there are no open boundaries. For the example given in Figure 2, the first node is given by i = 3 and j = 8, and the corresponding number is -1 because it is placed in a bottom margin (see also Fig. 4); XVI-initial i of the NNNN columns between columns to be swept in the x direction; XVII-final i of the NNNN columns between columns to be swept in the x direction; XVIII-j of the NNNN columns given by XVI and XVII. For the example given in Figure 2, the first column is given by j = 5 from i = 6 to i = 10 (see also Fig. 4); XIX-initial j of the MMMM rows between rows to be swept in the y direction; XX--final j of the MMMM rows between rows to be swept in the y direction; XXI-i of the MMMM rows given by XIX and XX. For the example given in Figure 2, the first row is given by i=4 from j=7 to j= 13 (see also Fig. 4). The bathymetry and the initial conditions are specified in a file whose name is given in the PARAMETR.DAT file. This file contains M lines with N data each. The h and q values must be given in centimeters, and the U and V values in cm/set. It should be noted that if the initial conditions are taken as still-water values, then U, I’, and r~ are zero everywhere. Table 1 shows the corresponding file for the example of Figure 2. The atmospheric forcing is specified in a file whose name is given in the PARAMETR.DAT file. The wind velocity (in knots) in the x and y direc- tions must be given in the nodes with even and odd ij values, respectively. The atmospheric pressure (in HPa) is given in the nodes with odd i and even j. The open boundary conditions are specified in a file whose name is given in the PARAMETR.DAT file. This file contains L lines, where ‘= KEND 3600 (At/2) +I with NE data each. Table 1. Example of bathymetry and initial condition data file for computational domain given in Figure 2 18 000000000000000000000 17 000000000000000000000 16 0 0 0 0 0 0 0 0 0 0 0 0 200 0 500 0 200 0 0 0 0 15 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 14 0 0 0 0 0 0 200 0 200 0 200 0 200 0 500 0 200 0 0 0 0 13 000000000000000000000 12 0 0 0 0 200 0 200 0 300 0 300 0 300 0 300 0 200 0 200 0 0 11 000000000000000000000 10 0 0 200 0 200 0 300 0 100 0 100 0 100 0 100 0 300 0 200 0 0 9 000000000000000000000 8 0 0 200 0 500 0 300 0 100 0 100 0 100 0 100 0 300 0 200 0 0 7 000000000000000000000 6 0 0 200 0 200 0 300 0 300 0 300 0 200 0 200 0 200 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 200 0 500 0 500 0 500 0 200 0 0 0 0 0 0 0 0 3 000000000000000000000 2 000000000000000000000 1 000000000000000000000 123456789 10 11 12 13 14 15 16 17 18 19 20 21 An improved 2-D finite-difference circulation model 1093 Table 2. Program outputs at specified times Times Free surface elevation Depth averaged Surface currents Half depth currents Bottom currents currents (0 = 0) (u = 0.5) (u = 1) TX1 ETAOUTOI .GRD VELOAOI .DAT VELOSO 1 .DAT VELOMOI .DAT VELOBOI .DAT TX2 ETAOUT02,GRD VELOA02.DAT VELOSOZ.DAT VELOM02.DAT VELOB02.DAT TX3 ETAOUT03,GRD VELOA03.DAT VELOS03.DAT VELOM03.DAT VELOB03.DAT TX4 ETAOUT04,GRD VELOA04,DAT VELOS04.DAT VELOMO4.DAT VELOBO4.DAT Main program The main program TAWIC reads the first data from the subroutines PARAM, INITI, and METE0 and after some preliminary calculations calls the subroutines BOUND, SWEEPX, SWEEPY, FILTER, ENERGY, and OUTPUT. The program then loops through the main subrou- tines SWEEPX and SWEEPY which are concerned with the AD1 algorithm for each time-step. Before calling these subroutines all the boundary con- ditions of the computational domain are updated on the subroutine BOUND. After a predetermined number of time-steps the subroutine FILTER is called and the calculated values of U, IJ’, and n are filtered. At the end of each loop the total mechan- ical energy of the system is computed in the subrou- tine ENERGY and stored in the ENERGY.DAT file. Optionally, the user can create a file with the time history of U, V, or q at any desired node. This output must be achieved by modifying the program just below where the subroutine ENERGY is called. Subroutines SWEEPX (SWEEPY): this subroutine performs the bulk of the work involved in applying the AD1 algorithm. The U (v) and 9 values are computed by forward sweep and backward substitution in the x direction 0, direction). GEGEX (GEGEY): this subroutine called from SWEEPX (SWEEPY) computes the elements of the tridiagonal matrix when the computational domain is sweeped in the x direction b direction). CALV (CALU): this subroutine called from SWEEPX (SWEEPY) computes the V (U) values by forward sweep and backward substitution in the y direction (x direction) after the U (V) and r) values were computed in SWEEPX (SWEEPY) at the same time-step, BOTTOM: this subroutine (called from GEGEX, GEGEY, CALV, and CALU) performs a point cal- culation of the bottom currents applying the Galerkin technique in the vertical. This subroutine is called only when the depths are less than 200 m. If the depths exceed 200 m the bed stress is com- puted from rbx = pkU,/m rby = pkVdm (36) OUTPUT: this subroutine performs the outputs of the program at those times specified in TXl, IOo_Cm/s 0 50 IOOKm N t 15Cm/s 0 50 IOOKm N t Figure 6. Steady-state surface-water elevation in cm and surface (A) and bottom (B) currents produced by a south- erly wind of 12.5 m/set in a closed rectangular basin. 1094 F. J. Caviglia and W. C. Dragani 0 IO 20Km. - lsOcllv% (A) . . . < A,, _- I..... I.... (Cl 5 (B) (D) Figure 7. Tidal amplitude in cm and tidal currents for semidiurnal lunar M2 tide at cotidal hours 3 (A), 6 (B), 9 (C), and 12 (D) corresponding times of flood, high, ebb, and low tide at Golfo Nuevo (Lat: 435, Long: 64”W), Argentina. TX2, TX3, and TX4. The outputs consist of five files (Table 2) in a format readable by the Golden Sofware SURFER version 4. The files with the GRD extension contain the free surface elevation for each specified time. The other files are SURFER post files with the depth average currents, and the surface (a = 0), half depth (a = OS), and bottom (a = 1) currents when depths on the computational domain are less than 200 m. EXAMPLES OF OUTPUT Examples of output are provided to demonstrate the program capabilities. All simulations were made on a PC 486 DX2 running at 66 MHz. Figure 5 shows a result of the model implementation for the idealized domain of Figure 2. The free surface el- evation and depth mean currents are shown after 37 h of simulation using a time-step of At = 180 set and a resolution of Ax = AJJ = 5 km (Fig. 3). The 1095 Figure 8. Steady-state sea level in cm and currents produced by southeasterly wind of 25 mjsec in Rio de la Plata (Lat: 37”S, Long: 55”W), Argentina. open boundary conditions were r~ = 30 cos wt for the nodes i = 3, j = 8, 10, and 12, and U = 10 sin wt for the nodes i = 16, j = 14 and 16, with w = 27r/l2 h. Figure 6A,B shows a closed rectangular domain with the free surface elevation, and the surface (A) and bottom (B) currents, respectively, which result from forcing the model with a constant wind of 12.5 mjsec from the south. The basin has dimension of 800 km (north-south) and 400 km (west-east), with 65 m depth. This simulation was made without considering the Coriolis force. The bottom currents have south directions as is expected. In order to apply the model under field con- ditions it was implemented for Golfo Nuevo (Lat: 43”S, Long: 64”W) and the Rio de la Plata (Lat: 37”S, Long: 55”W), both on the Argentine coast, to simulate tide-driven currents and wind-induced cir- culation, respectively. Figure 7A,B,C,D shows the sea level and currents for the cotidal hours 3, 6, 9, and 12, respectively, and corresponds to the times of flood, high, ebb, and low tide within Golfo Nuevo. This run was performed using a time-step of At = 180 set and a horizontal grid spacing of Ax = Ay = 2.75 km. Figure 8 shows the steady-state free-surface el- evation and the depth mean currents for the Rio de la Plata after two days of wind-induced circulation simulation for a southeasterly wind of 25 mjsec. This run was made using a time-step of At = 80 s with AX = 7.94 km and Ay = 10.87 km. The domain has depth values of 1 m in the northwest corner and greater than 5.000 in the southeast cor- ner. CONCLUSIONS Using the two-dimensional depth-averaged non- linear equations and three-dimensional linear equations, a quasi-three-dimensional numerical model was devised by coupling two modules com- puting the depth averaged current and vertical dis- tribution of the velocities. In contrast to conventional two-dimensional models, the bottom stress is no longer related to the average depth cur- rents but is computed from the bottom currents. The two modules are coupled through the bottom shear stress. The program includes the option of running the model using the conventional represen- tation of the bed shear stress based on the depth mean currents, which requires less computational time. The model described yields reliable knowledge of the two- and three-dimensional characteristics of the flow field for a wide range of shallow-water flows, and can be adapted easily to any domain of irregu- lar boundaries. The program and test datasets are available by anonymous FTP from the server IAMG.ORG. REFERENCES Chapman, D. C., 1985, Numerical treatment of cross-shelf open boundaries in a barotropic coastal ocean model: Jour. Physical Oceanography, v. 15, no. 8, p. 106& 1075. Davies, A. M., 1980a, On formulating a three-dimensional hydrodynamic sea model with an arbitrary variation of vertical eddy viscosity: Computer Method. Appl. Mech. Eng., v. 22, no. 2, p. 187-211. 1096 F. J. Caviglia and W. C. Dragani Davies, A. M., 1980b, Application of the Galerkin method Douglas, J. Jr., and Gunn, J. E., 1964, A general formu- to the formulation of a three-dimensional nonlinear hy- lation of alternating-direction methods: Numer. Math., v. 6, no. 6, p. 4288453. drodynamic numerical sea model: Appl. Math. Modelling, v. 4, no. 4, p. 245-256. Davies, A. M., 1986, A three-dimensional model of the northwest european continental shelf, with application to the M4 tide: Jour. Physical Oceanography, v. 16, no. 5, p. 797-813. Davies, A. M., 1988, On formulating two-dimensional ver- tically integrated hydrodynamic numerical models with an enhanced representation of bed stress: Jour. Geophys. Res., v. 93, no. C2, p. 1241-1263. Davies, A. M., 1991, On using turbulence energy models to develop spectral viscosity models: Continental Shelf Research, v. 11, no. 11, p. 1313-1353. Davies, A. M., and Furnes, G. K., 1980, Observed and computed M2 tidal currents in the North Sea: Jour. Physical Oceanography, v. 10, no. 2, p. 2377257. Davies, A. M., and Heaps, N. S., 1980, Influence of the Norwegian trench on the wind-driven circulation of the North Sea: Tellus, v. 32, no. 21, p. 164175. Davies, A. M., and Jones, J. E., 1993, On improving the bed stress formulation in storm surge models: Jour. Geophysical Res., v. 98, no. C4, p. 7023-7038. Hess, K. W., 1976, A three-dimensional numerical model of the estuary circulation and salinity in Narragansett Bay: Estuarine Coastal Marine Science, v. 4, no. 3, p. 3255338. Owen, A., 1980, A three-dimensional model of the Bristol Channel: Jour. Physical Oceanography, v. 10, no. 8, p. 1290-1302. Pearce, B. R., and Cooper, C. K., 1981, Numerical circula- tion model for wind induced flow: Jour. Hydraulics Division, Am. Sot. Civil Engineers, v. 107, no. HY3, p. 286302. Shapiro, R., 1970, Smoothing, filtering, and boundary effects: Rev. Geophys. Space Phys., v. 8, no. 2, p. 359- 387. WAMDI Group, 1988, The WAM model-A third gener- ation ocean wave prediction model: Jour. Physical Oceanoaraohv, v. 18. no. 12. P. 1775-1810. Weare, T. j.,‘1979, Errors arising from irregular bound- aries in AD1 solutions of the shallow-water equations: Intern. Jour. Numerical Methods Engineering, v. 14, no. 6, p. 921-931.