ࡱ> E@ dvbjbj +dn:::::::NN:::::dffffff$R-~::::: N1::::d d  0::0:. 0@F60d00|0NN:::::04"7 OcNNy{ NN{Test Program for 1-D Euler Equations PROGRAM testFluid ! ! Program to test numerical methods ! USE Input USE FluidArrays USE Trans USE Matrix USE Location IMPLICIT NONE ! CALL SetIC ! Set initial conditions CALL SetLocAr ! Set arrays to cross reference data structures CALL AllocMat ! Allocate Matrices CALL InitFluid ! Initialize dependent variables CALL Transient ! Run a transient ! STOP ! END Define Basic Data Types MODULE IntrType IMPLICIT NONE ! ! Set parameters that determine the precision of floating point ! arithmetic and number of digits in integers ! INTEGER, PARAMETER :: sdk = selected_real_kind(13,307) INTEGER, PARAMETER :: sik = kind(10000000) END MODULE IntrType Data Structures I find that the best way to start a program is to make a list of all the information that you will need. Determine its scope in the program (Global or Local), and establish data structures to carry this information. I usually begin with scalar data. MODULE ScalarDat USE IntrType ! ! Global scalar data needed by the program ! REAL(sdk) :: dt = 0.1_sdk ! time step REAL(sdk) :: time = 0.0_sdk ! time REAL(sdk) :: endTime = 1.0_sdk ! end time for the transient REAL(sdk) :: eps = 1.e-8_sdk ! convergence criterion REAL(sdk) :: residMax ! maximum factional residual INTEGER(sik) :: it ! current iteration number INTEGER(sik) :: itmax = 10 ! number of iterations INTEGER(sik) :: nsmax = 100000 ! maximum number of time steps INTEGER(sik) :: nstep = 0 ! current time step number INTEGER(sik) :: ncell = 10 ! number of computational volumes INTEGER(sik) :: nvar ! total number of independent variables INTEGER(sik) :: nbnd = 1 ! number of boundary cells INTEGER(sik) :: ivLB ! lower bound of velocity indices INTEGER(sik) :: ivUB ! upper bound of velocity indices CHARACTER(LEN=1) :: lBC, uBC ! Boundary condition flags ! 'p' for pressure boundary condition ! 'v' for velocity boundary condition LOGICAL :: converged ! test on iteration convergence END MODULE ScalarDat Physical Array Data MODULE FluidArrays USE IntrType ! ! Arrays containing state information about the fluid in 1-D flow ! All use SI units ! ! Arrays names ending with "N" are evaluated at the new (n+1) time level ! others are values at the start of the time step. ! REAL(sdk), ALLOCATABLE :: p(:) ! Start of step pressure REAL(sdk), ALLOCATABLE :: pN(:) ! End of step pressure REAL(sdk), ALLOCATABLE :: T(:) ! Start of step temperature REAL(sdk), ALLOCATABLE :: TN(:) ! End of step temperature REAL(sdk), ALLOCATABLE :: e(:) ! Start of step specific internal energy REAL(sdk), ALLOCATABLE :: eN(:) ! End of step specific internal energy REAL(sdk), ALLOCATABLE :: rho(:) ! Start of step density REAL(sdk), ALLOCATABLE :: rhoN(:) ! End of step density REAL(sdk), ALLOCATABLE :: rhoe(:) ! Start of step internal energy per vol REAL(sdk), ALLOCATABLE :: rhoeN(:) ! End of step internal energy per vol REAL(sdk), ALLOCATABLE :: v(:) ! Start of step cell edge velocity REAL(sdk), ALLOCATABLE :: vN(:) ! End of step velocity REAL(sdk), ALLOCATABLE :: q(:) ! power source per volume REAL(sdk), ALLOCATABLE :: fric(:) ! wall friction coefficient per volume ! ! Geometry arrays ! REAL(sdk), ALLOCATABLE :: dx(:) ! length of each computational volume REAL(sdk), ALLOCATABLE :: vol(:) ! fluid volume of each finite volume REAL(sdk), ALLOCATABLE :: fa(:) ! fluid flow area at each volume edge REAL(sdk), ALLOCATABLE :: hd(:) ! hydraulic diameter at each volume edge ! ! Pressure, temperature and velocity are independent variables. The ! TARGET attribute makes it easier to add perturbations from the solution ! arrays at the end of each iteration. ! TARGET :: pN, TN, vN CONTAINS SUBROUTINE StartStep IMPLICIT NONE ! ! Copy end-of-step variables from the last step into start-of-step ! variables for this time step ! rho = rhoN rhoe = rhoeN p = pN v = vN e = eN T = TN END SUBROUTINE StartStep SUBROUTINE InitFluid USE Eos IMPLICIT NONE INTEGER(sik) :: i ! ! Initialize all dependent fluid state variables ! DO i = LBOUND(p,DIM=1), UBOUND(p,DIM=1) rhoN(i) = rhoLiq (TN(i), pN(i)) eN(i) = SpEnergy(TN(i), pN(i)) rhoeN(i) = rhoN(i)*eN(i) ENDDO ! END SUBROUTINE InitFluid SUBROUTINE AllocFluidAr USE ScalarDat ! Allocate state variable arrays, space is reserved at each end for ! storage of boundary information. ALLOCATE(pN(1-nbnd:ncell+nbnd)) ALLOCATE(p(1-nbnd:ncell+nbnd)) ALLOCATE(TN(1-nbnd:ncell+nbnd)) ALLOCATE(T(1-nbnd:ncell+nbnd)) ALLOCATE(rho(1-nbnd:ncell+nbnd)) ALLOCATE(rhoN(1-nbnd:ncell+nbnd)) ALLOCATE(e(1-nbnd:ncell+nbnd)) ALLOCATE(eN(1-nbnd:ncell+nbnd)) ALLOCATE(rhoe(1-nbnd:ncell+nbnd)) ALLOCATE(rhoeN(1-nbnd:ncell+nbnd)) ALLOCATE(v(1-nbnd:ncell+nbnd+1)) ALLOCATE(vN(1-nbnd:ncell+nbnd+1)) ALLOCATE(q(1:ncell)) ALLOCATE(fric(1:ncell+1)) ALLOCATE(dx(1-nbnd:ncell+nbnd)) ALLOCATE(fa(1-nbnd:ncell+nbnd+1)) ALLOCATE(hd(1-nbnd:ncell+nbnd+1)) ALLOCATE(vol(1-nbnd:ncell+nbnd)) END SUBROUTINE AllocFluidAr END MODULE FluidArrays Array Data for Equation Solution MODULE Matrix USE IntrType ! ! Array storage for the linear system used in the Newton iteration ! REAL(sdk), ALLOCATABLE :: a(:,:), x(:), b(:) INTEGER(sik), ALLOCATABLE :: ipvt(:) CONTAINS SUBROUTINE AllocMat ! ! Allocate space for the linear system ! USE ScalarDat ALLOCATE ( a(nvar,nvar), x(nvar), b(nvar), ipvt(nvar) ) END SUBROUTINE AllocMat SUBROUTINE Solve ! Apply a linear solver appropriate to the matrix USE ScalarDat USE LUsolve ! IMPLICIT NONE INTEGER(sik) :: info CALL sgefat(a,nvar,nvar,ipvt,info) CALL sgeslt(a,nvar,nvar,ipvt,b,0) END SUBROUTINE Solve END MODULE Matrix Data Locators MODULE Location ! ! This contains several types of locaters. The first gives locations and ! associated weighting used to perform averages or differentials of ! quantities. The second provides pointers to move information from ! the data structure associated with the full system of equations and ! unknowns, to the data structure associated with physical state variables ! and spatial locations. The third is a set of translation indices to ! provide the full system variable index for each independent state ! variable. The fourth lists independent variables associated with USE Intrtype IMPLICIT NONE !======================================================================= ! Information for averages (and differences) TYPE averageT INTEGER(sik) :: n ! number of points used INTEGER(sik), POINTER :: i(:) ! cell or face index REAL(sdk), POINTER :: wf(:) ! weight factor END TYPE ! TYPE (averageT), ALLOCATABLE, TARGET:: fluxAv(:), dPdx(:), dV(:) ! fluxAv - information used to construct averages in cell edge ! mass and energy fluxes ! dPdx - information used to evaluate the pressure gradient ! in the momentum equation ! dV - information used to evaluate the velocity derivative ! in momentum transport term ! !======================================================================== ! Translation from system to physical variables TYPE varPtrT REAL(sdk), POINTER :: var INTEGER(sik) :: i, itype END TYPE ! %var - pointer to the variable storage in the fluid arrays ! %i - index in the fluid array ! %itype - variable type ! 1 - pressure ! 2 - temperature ! 3 - velocity ! TYPE (varPtrT), ALLOCATABLE :: loc(:) ! !======================================================================= ! Translation for physical to system variables ! INTEGER(sik), ALLOCATABLE :: ivarP(:), ivarT(:), ivarV(:) ! ! ivarP - index in the linear solution array "b" for the change ! in the corresponding element of pN ! ivarT - index in the linear solution array "b" for the change ! in the corresponding element of TN ! ivarV - index in the linear solution array "b" for the change ! in the corresponding element of vN ! !======================================================================== ! List of active varibles in each equation ! Also can be considered as a list of column numbers ! in each matrix row with potential for non-zero elements TYPE varListT INTEGER(sik) :: n, iloc INTEGER(sik), POINTER :: i(:) INTEGER(sik) :: eqnType END TYPE ! n - number of columns with potential for non-zero elements ! iloc - index of cell or face location for the equation ! i - list of column indices ! eqnType - Equation type associated with this row ! 1 - mass ! 2 - energy ! 3 - momentum ! TYPE (varListT), ALLOCATABLE :: varList(:) CONTAINS SUBROUTINE SetLocAr ! ! Allocate and load indices in the locator arrays ! ! Note that I only record potential nonzero contributions to matricies ! Depending on direction of velocity, some coefficients marked as potentially ! non-zero in structures below may actually be zero. ! USE ScalarDat USE FluidArrays IMPLICIT NONE INTEGER(sik) :: i, icol, irow, ivar, iLB, iUB, j, jj, nv ! ! Lower and Upper Bounds of velocity (cell edge) arrays ivLB = 1 ivUB = ncell+1 ! Calculate the total number of unknowns nvar = 3*ncell + 1 ! p & T in each volume, velocity at each edge IF (lBC.EQ.'v') THEN ivLB = 2 ! velocity fixed at left edge, do not evaluate nvar = nvar - 1 ! momentum equation at face 1 ENDIF IF (uBC.EQ.'v') THEN ivUB = ncell ! velocity fixed at right edge, do not evaluate nvar = nvar -1 ! momentum equation at face ncell+1 ENDIF ALLOCATE(loc(nvar),varList(nvar)) ALLOCATE(fluxAv(ncell+1), dPdx(ivLB:ivUB), dV(ivLB:ivUB)) ! ! Load tables used for weight factors ! DO i = 1,ncell+1 ! Number of volumes that might contribute to an edge average fluxAv(i)%n = 2*nbnd ALLOCATE (fluxAv(i)%i(fluxAv(i)%n)) ALLOCATE (fluxAv(i)%wf(fluxAv(i)%n)) ! Load indices of volumes that might contribute to an edge average DO j = -nbnd, nbnd-1 fluxAv(i)%i(j+nbnd+1) = i + j ENDDO ENDDO DO i = ivLB,ivUB ! Number of volumes contributing to the pressure derivative dPdx(i)%n = 2 ALLOCATE (dPdx(i)%i(dPdx(i)%n)) ALLOCATE (dPdx(i)%wf(dPdx(i)%n)) dPdx(i)%i(1) = i - 1 ! Indicies of volumes contributing dPdx(i)%i(2) = i ! to the derivative ! Weighting factors to generate the derivative. ! Note that calculation here requires a fixed mesh. dPdx(i)%wf(2) = 2.0_sdk/(dx(i)+dx(i-1)) dPdx(i)%wf(1) = - dPdx(i)%wf(2) ENDDO DO i = ivLB, ivUB ! Number of edges contributing to the velocity derivative dV(i)%n = 1+ 2*nbnd ALLOCATE (dV(i)%i(dV(i)%n)) ALLOCATE (dV(i)%wf(dV(i)%n)) ! Load indices of edges that might contribute DO j = -nbnd, nbnd dV(i)%i(j+nbnd+1) = i + j ENDDO ENDDO ! ! Translation of state to system variables iLB = LBOUND(p,DIM=1) iUB = UBOUND(p,DIM=1) ALLOCATE ( ivarP(iLB:iUB), ivarT(iLB:iUB)) iLB = LBOUND(v,DIM=1) iUB = UBOUND(v,DIM=1) ALLOCATE ( ivarV(iLB:iUB)) ivarP = 0 ivarT = 0 ivarV = 0 ivar = 0 ! variable index in the full system of equations IF (ivLB.EQ.1) THEN ivar = ivar+1 loc(ivar)%var => vN(1) loc(ivar)%i = 1 loc(ivar)%itype = 3 ivarV(1) = ivar ENDIF DO i = 1,ncell ivar = ivar+1 ! Pressure loc(ivar)%var => pN(i) loc(ivar)%i = i loc(ivar)%itype = 1 ivarP(i) = ivar ivar = ivar+1 ! Temperature loc(ivar)%var => TN(i) loc(ivar)%i = i loc(ivar)%itype = 2 ivarT(i) = ivar ! No system variable at right edge IF(i+1.GT.ivUB) EXIT ! if it has a constant velocity BC ivar = ivar+1 ! Velocity loc(ivar)%var => vN(i+1) loc(ivar)%i = i+1 loc(ivar)%itype = 3 ivarV(i+1) = ivar ENDDO ! Calculate and store indices of all variables used ! in each equation. This is a list of columns in ! each row of the matrix with potential for non-zero ! elements. ! irow = 0 IF(ivarV(1).NE.0) THEN irow = irow + 1 varList(irow)%n = 1 ! Velocity at this edge appears in the equation varList(irow)%iloc = 1 varList(irow)%eqnType = 1 ! ! Velocity at edge 2 may appear in the equation IF( ivarV(i+1).NE.0) varList(irow)%n = varList(irow)%n + 1 varList(irow)%n = varList(irow)%n + 2*nbnd ! Variables in Dp & 1/rho terms ALLOCATE(varList(irow)%i(varList(irow)%n)) j = 1 varList(irow)%i(j) = ivarV(1) ! Local Velocity varList(irow)%i(j+1) = ivarT(1) ! Temp in 1/rho varList(irow)%i(j+2) = ivarP(1) ! press in 1/rho and Grad P IF( ivarV(2).NE.0) varList(irow)%i(j+3) = ivarV(2) ENDIF ! DO i = 1,ncell irow = irow + 1 ! Mass Equation varList(irow)%iloc = i varList(irow)%eqnType = 1 varList(irow)%n = 2*( 2*nbnd+1) ! Temperatures and Pressures IF( ivarV(i).NE.0) varList(irow)%n = varList(irow)%n + 1 IF( ivarV(i+1).NE.0) varList(irow)%n = varList(irow)%n + 1 ALLOCATE(varList(irow)%i(varList(irow)%n)) nv = 0 IF(ivarV(i).NE.0) THEN ! Velocity for left edge flux nv = nv+1 varList(irow)%i(nv) = ivarV(i) ENDIF ! Pressure and Temperatures in time derivative and ! Flux terms DO j = -nbnd,nbnd nv = nv+1 varList(irow)%i(nv) = ivarP(i+j) nv = nv+1 varList(irow)%i(nv) = ivarT(i+j) ENDDO IF(ivarV(i+1).NE.0) THEN ! Velocity for right edge nv = nv+1 varList(irow)%i(nv) = ivarV(i+1) ENDIF irow = irow + 1 ! Energy Equation ! Same footprint as Mass Equation ALLOCATE(varList(irow)%i(varList(irow-1)%n)) varList(irow) = varList(irow-1) varList(irow)%eqnType = 2 ! ! Momentum Equation at the Right cell face IF(ivarV(i+1).EQ.0) CYCLE irow = irow + 1 varList(irow)%iloc = i+1 varList(irow)%eqnType = 3 varList(irow)%n = 1 IF( ivarV(i+2).NE.0) varList(irow)%n = varList(irow)%n + 1 IF( ivarV(i).NE.0) varList(irow)%n = varList(irow)%n + 1 varList(irow)%n = varList(irow)%n + 2*nbnd ! Variables in 1/rho term IF(ivarP(i+1).NE.0) THEN ! Dp and 1/rho contributions to right varList(irow)%n = varList(irow)%n + 2*nbnd ENDIF ALLOCATE(varList(irow)%i(varList(irow)%n)) varList(irow)%i = 0 nv = 0 IF( ivarV(i).NE.0) THEN nv = nv +1 varList(irow)%i(nv) = ivarV(i) ENDIF nv = nv + 1 varList(irow)%i(nv) = ivarV(i+1) ! Local Velocity DO j = -nbnd+1,nbnd ! 1/rho and Grad P terms IF(ivarP(i+j).EQ.0) CYCLE nv = nv + 1 varList(irow)%i(nv) = ivarT(i+j) nv = nv + 1 varList(irow)%i(nv) = ivarP(i+j) ENDDO IF( ivarV(i+2).NE.0) varList(irow)%i(nv+1) = ivarV(i+2) ENDDO END SUBROUTINE SetLocAr END MODULE Location Set Initial Conditions MODULE Input USE Intrtype USE ScalarDat USE FluidArrays ! ! This module sets initial conditions directly that would normally be ! obtained from an input file ! CONTAINS SUBROUTINE SetIC ncell = 12 dt = 0.1_sdk CALL AllocFluidAr DO i = LBOUND(v,DIM=1), UBOUND(v,DIM=1) vN(i) = 5.0_sdk fa(i) = 0.4_sdk ENDDO DO i = LBOUND(v,DIM=1), UBOUND(p,DIM=1) pN(i) = 1.0e5_sdk TN(i) = 300.0_sdk dx(i) = 0.5_sdk vol(i) = dx(i)*fa(i) ENDDO q = 6.98e8_sdk fric = 0.1_sdk lBC = 'v' uBC = 'p' END SUBROUTINE SetIC END MODULE Input Create a Subprogram for Each Flow Equation MODULE FlowEqn ! ! This contains functions needed to evaluate flow equation residuals ! USE IntrType USE Location USE FluidArrays USE ScalarDat Use Eos CONTAINS Start With the Mass Conservation Equation FUNCTION MassEqn(i) ! ! Residual of the mass equation at cell i ! IMPLICIT NONE ! ! Input ! ! i - cell index at which the equation is evaluated ! REAL(sdk) :: MassEqn INTEGER(sik), INTENT(IN) :: i MassEqn = (rhoN(i) - rho(i))*vol(i)/dt + Flux(rhoN,i+1) - & Flux(rhoN,i) END FUNCTION MassEqn ! Energy Conservation FUNCTION EnergyEqn(i) ! ! Residual of the energy equation at cell i ! IMPLICIT NONE ! ! Input ! ! i - cell index at which the equation is evaluated ! REAL(sdk) :: EnergyEqn INTEGER(sik), INTENT(IN) :: i EnergyEqn = (rhoeN(i)- rhoe(i))*vol(i)/dt + & Flux(rhoeN,i+1) - Flux(rhoeN,i) & + pN(i)*(fa(i+1)*vN(i+1)-fa(i)*vN(i)) - q(i)*vol(i) END FUNCTION EnergyEqn Momentum Conservation FUNCTION MomenEqn(i) ! ! Residual of the momentum equation at face i (low edge of cell i) ! IMPLICIT NONE ! ! Input ! ! i - edge index at which the equation is evaluated ! REAL(sdk) :: MomenEqn INTEGER(sik), INTENT(IN) :: i MomenEqn = (vN(i)- v(i))/dt + Vdv(i) + GradP(i)/EdgeAv(rhoN,i) + & fric(i)*vN(i)*ABS(vN(i)) END FUNCTION MomenEqn ! Average to Get Cell Edge Values for State Variables FUNCTION EdgeAv(x,j) ! ! Calculate the value of scalar field x at cell face j based on ! weighting coefficients in fluxAv ! IMPLICIT NONE REAL(sdk) :: EdgeAv ! quantity being fluxed REAL(sdk), INTENT(IN) :: x(LBOUND(rho,1):UBOUND(rho,1)) INTEGER(sik), INTENT(IN) :: j ! edge index TYPE (averageT), POINTER :: a INTEGER(sik) :: i a => fluxAv(j) EdgeAv = 0 DO i = 1,a%n EdgeAv = EdgeAv + a%wf(i)*x(a%i(i)) ENDDO END FUNCTION EdgeAv ! Calculate Flux of Mass or Energy at Cell Edge FUNCTION Flux(x,j) ! ! Calculate the flux of scalar field x at cell face j based on ! weighting coefficients in fluxAv ! IMPLICIT NONE REAL(sdk) :: Flux ! Quantity being fluxed REAL(sdk), INTENT(IN) :: x(LBOUND(rho,1):UBOUND(rho,1)) INTEGER(sik), INTENT(IN) :: j ! edge index TYPE (averageT), POINTER :: a INTEGER(sik) :: i a => fluxAv(j) Flux = 0 ! First calculate edge value of x DO i = 1,a%n Flux = Flux + a%wf(i)*x(a%i(i)) ENDDO Flux = Flux*vN(j)*fa(j) ! Multiply edge value by edge area and velocity END FUNCTION Flux ! Special Function for Momentum Transport FUNCTION Vdv(j) ! ! Calculate the momentum flux term at face j based upon weighting ! information in derived type array dV ! IMPLICIT NONE REAL(sdk) :: Vdv INTEGER(sik), INTENT(IN) :: j ! Face where momentum Equation is evaluated INTEGER(sik) :: i TYPE (averageT), POINTER :: a a => dV(j) Vdv = 0 DO i = 1,a%n Vdv = Vdv + a%wf(i)*vN(a%i(i)) ENDDO Vdv = vN(j)*Vdv END FUNCTION Vdv Special Function for Pressure Gradient FUNCTION GradP(j) ! ! Calculate an approximation to the pressure gradient at face j ! based on information in the dPdx array ! IMPLICIT NONE REAL(sdk) :: GradP TYPE (averageT), POINTER :: a INTEGER(sik), INTENT(IN) :: j ! Face where momentum Equation is evaluated INTEGER(sik) :: i a => dPdx(j) GradP = 0 DO i = 1,a%n GradP = GradP + a%wf(i)*pN(a%i(i)) ENDDO END FUNCTION GradP END MODULE FlowEqn Drive the Transient Solution MODULE Trans USE IntrType ! ! Subprograms used in the flow transient ! CONTAINS SUBROUTINE Transient ! ! Driver for the transient ! USE ScalarDat USE FluidArrays USE Matrix IMPLICIT NONE time = 0.0_sdk DO nstep = 1,nsmax CALL StartStep converged = .FALSE. CALL SetWf CALL Residuals(0) DO it = 1,itmax CALL SetWf CALL Jacobian CALL Solve CALL EvalVar CALL Residuals(1) IF(converged) EXIT ENDDO IF(.NOT.converged) THEN WRITE (*, '(a,f7.2)') 'Iteration failed at time ',time STOP ENDIF time = time + dt IF( time.GT. endtime-0.01*dt) EXIT ENDDO END SUBROUTINE Transient SUBROUTINE EvalVar USE FluidArrays USE Location USE Matrix USE ScalarDat USE Eos IMPLICIT NONE INTEGER(sik) :: i ! ! Update the state variable arrays ! DO i = 1,nvar loc(i)%var = loc(i)%var + b(i) ENDDO DO i = 1,ncell rhoN(i) = rhoLiq (TN(i), pN(i)) eN(i) = SpEnergy(TN(i), pN(i)) rhoeN(i) = rhoN(i)*eN(i) ENDDO END SUBROUTINE EvalVar SUBROUTINE Residuals(mode) ! ! Evaluate Mass, Energy, and Momentum equation residuals ! ! Input ! ! mode - 0 only evaluate residuals ! 1 find maximum scaled residual ! USE Matrix USE FluidArrays USE Location USE FlowEqn USE ScalarDat ! IMPLICIT NONE INTEGER(sik), INTENT(IN) :: mode INTEGER(sik) :: i, ivar ! DO i = 1,nvar SELECT CASE (varList(i)%eqnType) CASE(1) ! Mass Equation b(i) = MassEqn(varList(i)%iloc) CASE(2) ! Energy Equation b(i) = EnergyEqn(varList(i)%iloc) CASE(3) ! Momentum Equation b(i) = MomenEqn(varList(i)%iloc) END SELECT ENDDO ! IF(mode.EQ.0) RETURN ! ivar = 0 residMax = 0.0_sdk IF(ivarV(1).NE.0) THEN ivar = ivar+1 residMax = MAX(residMax,b(ivar)/MAX(0.01, v(1), vN(1))*dt) ENDIF DO i = 1,ncell ivar = ivar+1 residMax = MAX(residMax,b(ivar)/(vol(i)*MAX(0.01, rho(i), rhoN(i)))*dt) ivar = ivar+1 residMax = MAX(residMax,b(ivar)/(vol(i)*MAX(0.01, rhoe(i), rhoeN(i)))*dt) IF(ivarV(i+1).EQ.0) CYCLE ivar = ivar+1 residMax = MAX(residMax,b(ivar)/MAX(0.01, v(i+1), vN(i+1))*dt) ENDDO ! IF(residMax.LT.eps) converged = .TRUE. ! END SUBROUTINE Residuals SUBROUTINE Jacobian ! ! Evaluate the Jacobian Matrix using finite approximations to each ! partial derivative ! USE Matrix USE FluidArrays USE ScalarDat USE FlowEqn USE Eos USE Location ! IMPLICIT NONE INTEGER(sik) :: i, ii, ivar, j REAL(sdk) :: fbase,varBase, rhoBase, rhoeBase, fracDvar=0.001_sdk, dvar REAL(sdk) :: fpert ! First clear the Jacobian Matrix a = 0.0_sdk ! Now the set non-zero Jacobian elements ! Loop over all rows in the matrix DO i = 1,nvar DO j = 1,varList(i)%n ! Loop over all non-zero columns in the row ivar = varList(i)%i(j) IF(ivar.EQ.0) CYCLE varBase = loc(ivar)%var ! hold the old variable value ! ! Perturb the variable ! IF(loc(ivar)%itype.EQ.3) THEN ! velocity IF(loc(ivar)%var.GE.0) THEN dvar = MAX(fracDvar, fracDvar*loc(ivar)%var) loc(ivar)%var = loc(ivar)%var + dvar ELSE dvar = MIN(-fracDvar, fracDvar*loc(ivar)%var) loc(ivar)%var = loc(ivar)%var + dvar ENDIF ELSE ! pressure or temperature dvar = fracDvar*loc(ivar)%var loc(ivar)%var = loc(ivar)%var + dvar ii = loc(ivar)%i rhoBase = rhoN(ii) rhoeBase = rhoeN(ii) rhoN(ii) = RhoLiq(TN(ii),pN(ii)) rhoeN(ii) = rhoN(ii)*SpEnergy(TN(ii),pN(ii)) ENDIF ! Reevaluate the function SELECT CASE (varList(i)%eqnType) CASE(1) ! Mass Equation fpert = MassEqn(varList(i)%iloc) CASE(2) ! Energy Equation fpert = EnergyEqn(varList(i)%iloc) CASE(3) ! Momentum Equation fpert = MomenEqn(varList(i)%iloc) END SELECT ! Jacobian a(i,ivar) = (b(i) - fpert)/dvar ! loc(ivar)%var = varBase ! Restore variables IF(loc(ivar)%itype.NE.3) THEN ii = loc(ivar)%i rhoN(ii) = rhoBase rhoeN(ii) = rhoeBase ENDIF ENDDO ENDDO END SUBROUTINE Jacobian SUBROUTINE SetWf ! ! Set weighting factors for edge averages and differences ! Contents of this subroutine are specific to the 1st order upwind ! method used as a sample ! USE FluidArrays USE ScalarDat USE Location ! IMPLICIT NONE INTEGER(sik) :: i DO i = 1,ncell+1 IF(vN(i).LT.0.0_sdk) THEN fluxAv(i)%wf(1) = 0.0_sdk fluxAv(i)%wf(2) = 1.0_sdk IF(i.LT.ivLB.OR.i.GT.ivUB) CYCLE dV(i)%wf(1) = 0.0_sdk dV(i)%wf(2) = dPdx(i)%wf(1) dV(i)%wf(3) = dPdx(i)%wf(2) ELSE fluxAv(i)%wf(2) = 0.0_sdk fluxAv(i)%wf(1) = 1.0_sdk IF(i.LT.ivLB.OR.i.GT.ivUB) CYCLE dV(i)%wf(1) = dPdx(i)%wf(1) dV(i)%wf(2) = dPdx(i)%wf(2) dV(i)%wf(3) = 0.0_sdk ENDIF ENDDO END SUBROUTINE SetWf END MODULE Trans Equation of State MODULE Eos USE IntrType ! ! Equation of state information and other physical properties ! REAL(sdk) :: mu = 0.00086_sdk ! water viscosity (Pa*s) REAL(sdk) :: cond = 0.62_sdk ! water conductivity (w/m/K) REAL(sdk) :: cv = 4186.0_sdk ! water specific heat (J/kg/K) REAL(sdk) :: dRhoDt =-0.2538_sdk ! derivative of density with ! respect to temperature (kg/m**3/K) CONTAINS FUNCTION RhoLiq(T,p) ! ! Water density ! ! Input ! T - Temperature (K) ! p - pressure (Pa) ! Output ! RhoLiq - density (kg/m**3) ! IMPLICIT NONE REAL(sdk), INTENT(IN) :: T, p REAL(sdk) :: RhoLiq ! RhoLiq = 923.86 + dRhoDt*T ! END FUNCTION RhoLiq FUNCTION SpEnergy(T,p) ! ! Water Specific Internal Energy ! ! Input ! T - Temperature (K) ! p - pressure (Pa) ! Output ! SpEnergy - Specific Internal Energy (J/kg) ! IMPLICIT NONE REAL(sdk), INTENT(IN) :: T, p REAL(sdk) :: SpEnergy ! SpEnergy = cv*T ! END FUNCTION SpEnergy END MODULE Eos %)   K M ] ^ Z [ xz:<]^ ! !!´¦ї{i[iL>h 7h 75CJ0^JaJ0h 75CJOJQJ^JaJhk0hk05CJ(^JaJ(#hk0hk05CJOJQJ^JaJhk0h5CJ0^JaJ0hk0hk05CJ0^JaJ0hk05CJOJQJ^JaJhh5CJ^JaJhh5CJ(^JaJ(h5CJOJQJ^JaJ#hh5CJOJQJ^JaJhCJOJQJ^JaJhh5CJ0^JaJ0%&'()?Agiw L {    @8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd$a$gddv   , B D 2 L ] ^ Z [ p  T S@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8$a$gdgd!o@_`y\T@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gdk0$a$gdk0gd3rKbdxz`Enp@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gdk0 +-<MXcnyzHn @8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gdk0  Syz-Sv 2Ki ;]^@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8$a$gdk0gdk0^p$%23KMvx> P ` b t @8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gdk0 !!!*!,!x!!"M""")#n#o#####$$6$x$$$$$C%%@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7$a$gdgdk0!!&&**IIIIIJLLL_LvLwLxL3M9MbMNNNNPPPPPPPJRLRRRTT𨚌}k}]k}}k}k}}k}}k}hQPhQP5CJ(^JaJ(#hQPhQP5CJOJQJ^JaJhQP5CJOJQJ^JaJhQPh 75CJ0^JaJ0hQPhQP5CJ0^JaJ0hQP5CJ0^JaJ0h 7hk05CJ(^JaJ(h 7h 75CJ(^JaJ(h5CJOJQJ^JaJ#h 7h 75CJOJQJ^JaJh 75CJOJQJ^JaJ$%%%&_&&&&!'2'R'q'~'''(9(`((((((9);)z)|)))D*@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7D*z****F+++,",@,d,,,,$-R----. .9.:.G.H.`.b...@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7..1/h/j/|////////0<0O0000K111111>2H2I2o222@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 72222@3[33334>4J4T4U4j4444#5o555556c666666@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 76717T777777(8B8\8]8888888888 9 9Y9Z9t9999@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 799999 :8:U:k::::::;;0;v;;;; <"<<<T<U<_<<<6=@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 76=V=X=e====>">q>>??5?A?B???@K@U@W@j@@@@*AiAAA@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7AA2BDBkBwBBBBC:CLCuCCCCCDDHDDDD E EVEvEEEE@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7EE&FeFFFG7GCGtGGGGGGHHTHUHHHHHI8IDIIIII@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gd 7IIIIIIIJJJfJJJJJJJJJK$K:KDKpKKKKKK@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8$a$gdQPgd 7KKLLL6LKLwLxLLLLLLLMM&M'McMdM|M~MMMMMM@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8t@8@8@8@8@8@8@8@8gdQP$a$gdQPgd 7MMN N"NDNENNNNNNNNNN O O;O=OOOQO[O]OOOOOO@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8$a$gdQPgdQPO P?PPPPPPPPPPPPQQ1Q3Q=Q?QyQ{QQQQR)R*RDR@8@8@8@8@8@8@8@8Y@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8$a$gdQPgdQPDRFRKRRRRRRRRSSS/SbSSSTTT-T[U[s[t[[[[[[[[[[@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8@8gdzgdQP[ \ \\3\B\C\U\V\i\\\\\\\]]*]?]Y]t]]]]]] ^6^@^gdz@^A^^^_^v^^^^^^^^^^__,_Q_[_\_o_______``!`gdz!`]`_`j`l``````` aa!a3aXatavaaaabXbbbbccccgdzc+c-c:cQclcccccc@dTdddde'e)eTeVeseteweeeeeeegdzef%f5fAfRfTffffff-g>gggg)hHhdhhhhhh$i]iiiijgdzjj[jjjjjk@kykkkk*lXlllm=mPmmmmmn:nWnvnnngdznnnnnnn oPolonooooooooop'pIprpppppq'qPqgdzPqtqqqqqqqqqrr/r1rqrsrrr8szsssssstt t*t$a$gdzgdz*tEtPtvtxttttttttuuu"u$uIuKuUuruuuuuuuvvvgdzv3v5vOvQvdvgdz)1h0:p= /!"#$%H@H Normal CJ$OJQJ_HaJ$mH sH tH DA@D Default Paragraph FontRi@R  Table Normal4 l4a (k@(No Listdn%&'()?AgiwL{,BD2L]^Z[pTS!o@ _ ` y   \  T 3 r K bdxz`Enp +-<MXcnyzHn  Syz-Sv 2Ki ;]^p$%23KMvx>P`bt *,xM)no$6xC_!2Rq~ 9 ` 9!;!z!|!!!D"z""""F###$"$@$d$$$$$%R%%%%& &9&:&G&H&`&b&&&&1'h'j'|'''''''/(<(O((((K))))))>*H*I*o******@+[++++,>,J,T,U,j,,,,#-o----5.c....../1/T//////(0B0\0]0000000000 1 1Y1Z1t11111111 282U2k2222223303v3333 4"4<4T4U4_44465V5X5e55556"6q667757A7B7778K8U8W8j8888*9i99992:D:k:w::::;:;L;u;;;;;<<H<<<< = =V=v=====&>e>>>?7?C?t??????@@T@U@@@@@A8ADAAAAAAAAAAABBBfBBBBBBBBBC$C:CDCpCCCCCCCDDD6DKDwDxDDDDDDDEE&E'EcEdE|E~EEEEEEEF F"FDFEFFFFFFFFFF G G;G=GOGQG[G]GGGGGG H?HHHHHHHHHHHHII1I3I=I?IyI{IIIIJ)J*JDJFJKJJJJJJJJKKK/KbKKKLLL-LSUSsStSSSSSSSSSS T TT3TBTCTUTVTiTTTTTTTUU*U?UYUtUUUUUU V6V@VAV^V_VvVVVVVVVVVVWW,WQW[W\WoWWWWWWWXX!X]X_XjXlXXXXXXX YY!Y3YXYtYvYYYYZXZZZZ[[[[+[-[:[Q[l[[[[[[@\T\\\\]'])]T]V]s]t]w]]]]]]]^%^5^A^R^T^f^^^^-_>____)`H`d``````$a]aaaabb[bbbbbc@cycccc*dXddde=ePeeeeef:fWfvfffffffff gPglgngggggggggh'hIhrhhhhhi'iPitiiiiiiiiijj/j1jqjsjjj8kzkkkkkkll l*lElPlvlxllllllllmmm"m$mImKmUmrmmmmmmmnnn3n5nOnQnfn000p0p00000000000000000000000p00p0p0000000000p0p00p0000000000000 0000 000 00 000 0 0p00000000000000000000000000000 00(0 :0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F:0F000000000000000000000000000000 0 0(0 0p0000000000000000000000000000 0 0 00p00000000000000000000000000000 000000000000 0 0 0 0 0 0 00 000000000000 00 00 00 0 0 0 0 0 00 00 0 0 0 00 000 0 0 0 0 0 0 000000 00 00 0(0(0(0(0(0(0(0(00(0(00(00(0(0(0(0(0(0(000(000(0 0 0(000(0 0(0(0(0(0(0(0(00(0(0(000(0(00000000000000000000000000000000000000000000000000000000000000000000000000000 0 0 00 0 0000000000000000000000000000000808080808080800080 0 00 0 0 0 0808080808080 0080800808008080808080808080800800800 0 0 0 0 0 0808080@0@0@00@0@00@0000@0@0@00p00000000000000000000000000000 000p0@0@0@0@0@00@0@000@0@0@0@0@0@0@0@0@0000@0@0@0@0@0 0 000p0@0p0p0p0@0@0@00@0@0@00@0@0@0H0H00H0H0H0H00p0p0p0H0p0H0H00H0H0H0H0H000H0H0H0H0H00H0p0p0p0p0H0p0H0H00H0H0H0H0H000H0H0H0H0H00H00000p0H0H00 0 0H0(00 0 0 00H0H0H0H0H0H000H0H0H0H000p0p0H0H0H0H00H0H0H0H0H0H0H0H0H0 0P0 00 000p0p0p0P0P0P0P0P0P00P0P0P000P0P0P000P0p0p0p:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:000P0P0P0P0P0P000P00P0P00000000000000000000 0 00:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:0:000X0X0X0X0X0X000`0`000000000000`000000000 0 00(0`0`0`0`0`0`0`000`00`0`00`0`00`0`0`0`0`00`00`00`00`0 0 0 0 0`0h0h0h0h0h000h0h00h0h00h0h0h00 00p0p0h0h0h0h00h00h0h00h0h0h000000h00000000000 0 0 0 0h0h0h0h0h00h0h00h0h00h0h000h,2p _ `  p *AAAAAADD6DDHHHHHHS%S&S=S>SSiiiiijn3n5nOnQnfn:0C:0:0\H:0@0:0:0:0:0@0:0 Ř:0 :0 :0:0@0:0Pt:0:0:0:0@0:0Y:0:0:0:0@0:0d:0:0:0@0@0:0:0:0@0:0'0z:0':0':0':0+:0':0'|:0':0':0:0@0:03{:03:03:0:0@0:09:09:09:0:0:0!Tdv<FV  ^ %D*.2696=AEIKMODRTWY[@^!`cejnPq*tvdv=?@ABCDEGHIJKLMNOPQRSTUWXYZ[\]^_`abdv>8@0(  B S  ?3 3 fnfn8*urn:schemas-microsoft-com:office:smarttagsCity9*urn:schemas-microsoft-com:office:smarttagsplace 5> U]#+ CKfox")]`eh"%*/_bgl-059{~ X[]`o x   ] `   < ? Q T { ~   ! G J T W  2613il~~147;@DGLUW`bkm#$.3?DNRSTX^cdgijktvwx| :>X\vx/:x  BJFOX_w-5DG)/48=?HN!*1=@NQ`chikp G!J!]!b!g!l!q!v!!!!!" """""g#o#$!$0$3$;$?$N$Q$_$`$r$u$z$$$$3%4%a%h%&&.&5&W&_&&&r'{'''''''''''''''''''''3(7(@(D((((((())Q)U)X)\)))))))*** *Z*^*`*g*h*l*|************F+L+M+N+V+Z+k+q+r+s+u+v+w+}+~++++++++++++ ,, ,&,',(,*,+,8,9,\,],`,i,,,,,,,,,,,,,,, ----------)---.-/-1-2-8-9-N-V-u-y-z-{-}-~---;.?.@.A.C.E.T.V.W.X.Z.\.i.m.n.o.q.s.{................ ///!/"/#/%/&/'/)/*/+/A/C/D/E/G/I/J/L/M/N/////////////,0/090>0F0I0S0X0l0q0r0y0|000000000000000000001111z1~11111111111111111111122B2F2H2K2O2Q2R2S2_2c2e2f2i2j2u2y2{2222222222222222222222 3333 3%3&3'3+3/3333333444444,4042474B4G4O4S4\5`5l5q555555555555555666666{66666666666666666666667777!7"7#7*7+7/7H7O7P7T7V7W7]7b777777777777777778!8+828387898:8B8G8^8_8p8t8w8{88888888888888888884999:9;9C9J9K9O9U9\9]9a9s9x9999999999999999999999999::<:L:S:T:X:Z:[:\:^:b:g:h:i:::; ;; ;!;%;';(;);+;/;4;5;8;B;D;T;[;\;`;b;c;d;f;j;o;p;s;;;;;;;;;;;;;;;<<< <<<<<<<<<<<<<<<<<<<<=_=d=|======================>>> >>>>>0>5>6>7>?>F>G>K>Q>X>Y>]>k>r>s>w>}>>>>>>>>>> ?????%?&?*?2?6?R?Y?Z?^?`?a?b?i?j?n?z????????????????????????????@@ @@@ @!@%@'@(@)@+@/@4@@@@@@@@@@@@@@@@@@@@@AAA AAAA#A%A&A'A)A-A2A3A6ANASA_AfAgAkAmAnAwA|AAAAAABBBBBBBBBBBBBBBC CCCCC*C,C-C.CKCLCVC[CgClCvCxCyCzCCCCCCCCCCCCCCCCCCCCCDDDD0D5DDDDDDEEEqExEyEzEEEEEFFF!F.F1FBFCFIFPFTFXFYFZF^FaFbFcFfFiFjFkFmFoFFFFFFGGG9G:GaGbGGGGGGGGGGGGGGGGGGGGGGHHH2H9HQHSHTHUHXHZH`HbHhHjHkHlHnHpHqHrHyHzH|HHHHHHHHHHIIIICIDIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIJJJJJJJJ"J$J%J&J;JCJJJJJJK K#K(K.KkKnKKKKKLLLL#L)L1L7LCLDLSLYL\LbLeLiLjLkLoLrLsLtLLLLLHMNMlMoMMMMM9NAN]N`NeNfNqNwNNNNNNNNNNNOOO OOOP PCPFPKPNP[P^PPPPPPPPPPPPP Q QQQQQQQQQ Q#Q$Q%Q8Q;Q>Q@QDQGQZQ]QQQRR,R/R4R9RDRLRhRkRRRRRRRRRRRRRSSS S SSSSSSSS7SUUU V VVV-V/VnVuV~VVVVVVVV!W"W6W7W9W[F[X[][r[v[[[[[[[[[[[[[[[[\\\\\\\\\*\-\.\/\2\6\7\8\<\>\F\J\Z\b\i\s\t\x\{\~\\\\\\\\\\\\\\\\\\\\\\\]]]]0]?]]]]]^^^$^-^4^r^u^z^{^^^^^^^^^^^^^^^^^^^_%_t_|_____1`5`8`?`@`A`C`D`l`s`z`~`````aa0a4a;aCaEaMaRaVaXa[amaqasava}aaaaaaaaaaaaaaaaaaaaaaabbbgbkbnbvb{bbbbbbbbbbbbbbbbbbbbbbbc ccc#c*c0c8c:cLcQcXc\cacicqcsccccccc6d;d?dFdGdNdOdPdRdVdddddddddddee#e+e,e3e4e5e7e;eeeeeeeeeeeeeeeeef f2f6f8f9fDfHfOfVfafffmfufffffvgggggggggggggg hhhhhh/h5h6h7h9h;hThjhzh|h}h~hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh iiiiii2iHiXiZi[i\i^i`ifijikilinipi|i~iiiiiiiiiiiiiiiiiiijj&j.j|jjjjjjjjjkkkAkDkIkOkkkkkWl]lllllllllllllmmm mmmmm nnnn#n+n.n0nFnNnfn 9<X]|">DW_{%-INs{ %DN & A E P [ ` f   X ] 7 < v {   O T ~diIOt|14@D'.NStw(/X_~ 1:W`z6?OXmvx|  QY-5<Dag~&(HQ!oq8=X` - 8 Q _ x ?!G!!!!!" "W"Y"""""##($0$F$N$j$r$$$$$;%?%a%h%%%%%%%&,&5'='''3(7(@(D((((()))&)Q)U)l)t)))))***$*M*V*s*|***F+M+k+r+++ ,',d,i,,,,, --)-.-u-z---;.@.i.n.../"/A/D///,0/0F0I0j0q000000000001111*121z1~1111111111222>2B2[2_2q2u2222222223 3 3&3|3333333344(4,4B4H44455L5T5\5`5i5l555555566w6z66666 77H7P777777788c8i8p8t88888880939o9r9999999::<:L:T:::; ;;!;B;D;T;\;;;;;;;<<<<<<<<\=_=|==========,>/>k>s>>>>> ??I?R?z??????????@@@!@i@n@@@@@@@AAAAJAMAkBsBBBBBBBCC*C-COCVCvCyCCCCCCCCCDDDD?ECEqEyEEEEE FF&F.FTFYFFFFGG GiGmGGGGGGGH HQHTHHHHHKIOIIIIIIIJJJJJJK KLKTKKKKKKKLLLLHLJLgLjLLL.M7MgMlMMM$N(NDNMNUN]NlNmNNNNNOOOOOP>PCPSP[PPPPPPPPQQQ>QAQQQQQ'R,RORXR`RhRRRRRRR S SZT^TyTTTTTTTTLUVUaUdUUUUUUVVVVV&W+W2W6WhWnWuWzWWWWWXXqXuXXX7Y?Y\YdYYYYYYYYY Z%ZbZdZZZZZ[[1[5[>[F[U[X[r[v[[[[[[[F\J\\\\\\\-]0]]]j^r^^^^^1_5_______1`5`P`S`l`s```````a a0a4aiamaaaaaCbKbgbkbbbbbbbbcc$cLcRccccc6d;dbdgdddddeeeeeeff)f+fDfIfafgfTgZggggggg hh/h6hQhThzh}hhhhhhh ii/i2iXi[i|iiiijjwj|jjjjjjjjjk k?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcefghijkmnopqrstuvwxyz{|}~Root Entry F Data d1TablelWordDocument+SummaryInformation(DocumentSummaryInformation8CompObjj  FMicrosoft Word Document MSWordDocWord.Document.89q