C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ c *dvlp/t2cg2.f* C********************************************************************* C C C PROGRAM TOUGH2 C C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C @ @ C @ TOUGH2, MODULE T2CG2, VERSION 2.0, October 1999 @ C @ @ C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C C-----PROGRAM MAIN IS A HIGH-LEVEL EXECUTIVE ROUTINE. C C MAIN CALLS SEVERAL LOWER-LEVEL EXECUTIVE ROUTINES, WHICH EXECUTE C THE ACTUAL CALCULATIONS. C C ALL LARGE ARRAYS ARE DIMENSIONED IN PARAMETER STATEMENTS IN AN C INCLUDE FILE 'T2'. IF MODIFICATIONS IN ARRAY DIMENSIONS ARE C DESIRED, THESE NEED TO BE DONE ONLY IN FILE 'T2'. C C======================================================================= C C GLOSSARY OF MAJOR ARRAYS, IN ORDER OF APPEARANCE BELOW C C======================================================================= C C ELEMENTS: C *ELEM* HOLDS ELEMENT CODE NAMES (IDENTIFIERS) C *MATX* HOLDS MATERIAL DOMAIN IDENTIFIER C *EVOL* HOLDS VOLUME C *PHI* IS POROSITY C *P* IS PRESSURE AT LAST CONVERGED TIME STEP C *T* IS TEMPERATURE AT LAST CONVERGED TIME STEP C *AI* IS FOR HEAT EXCHANGE WITH CONFINING LAYERS C *AHT* ARE HEAT TRANSFER AREAS C C C PRIMARY VARIABLES: C *X* LATEST CONVERGED VALUES OF PRIMARY VARIABLES C *DX* LATEST INCREMENTS OF PRIMARY VARIABLES C *DELX* SMALL INCREMENTS OF PRIMARY VARIABLES C FOR NUMERICALLY COMPUTING DERIVATIVES C *R* LATEST RESIDUALS C *DOLD* ACCUMULATION TERMS AT BEGINNING OF TIME STEP C *ROW* SCALING FACTORS FOR MATRIX ROWS C *COL* SCALING FACTORS FOR MATRIX COLUMNS C C C CONNECTIONS (INTERFACES): C *NEX1* INDICES OF FIRST ELEMENTS IN A CONNECTION C *NEX2* INDICES OF SECOND ELEMENTS IN A CONNECTION C *DEL1* FIRST ELEMENT NODAL DISTANCES C *DEL2* SECOND ELEMENT NODAL DISTANCES C *AREA* INTERFACE AREAS C *BETA* COSINE OF ANGLE BETWEEN NODAL LINE AND C THE VERTICAL C *ISOX* ISOTROPY INDICES C *GLO* HEAT FLOW RATES C *ELEM1* CODE NAMES OF FIRST ELEMENT IN A CONNECTION C *ELEM2* CODE NAMES OF SECOND ELEMENT C C C LINEAR EQUATION SOLVER: C *IRN* ROW INDICES OF NON-ZERO MATRIX ELEMENTS C *ICN* COLUMN INDICES OF NON-ZERO MATRIX ELEMENTS C *WKAREA* WORKSPACE USED IN MA28 C *IKEEP* INTERNAL STORAGE SPACE FOR MA28 C *JVECT* COLUMN INDICES OF NON-ZERO MATRIX ELEMENTS C *IW* INTERNAL STORAGE SPACE FOR MA28 C *W19A* WORKSPACE FOR SCALING ROUTINE (MC19A) C C C SINKS AND SOURCES: C *F1* TABLE OF GENERATION TIMES C *F2* TABLE OF GENERATION RATES C *F3* TABLE OF INJECTION ENTHALPIES C C C OTHER LARGE ARRAYS: C *PAR* SECONDARY (THERMOPHYSICAL) PARAMETERS FOR C ALL ELEMENTS C *FLO* RATES OF GAS AND LIQUID PHASE FLOWS ACROSS C ALL INTERFACES C *VEL* PORE VELOCITIES OF GAS AND LIQUID PHASE FLOW C C C======================================================================= C C$$$$$$$$$ ASSIGN PARAMETERS FOR FLEXIBLE DIMENSIONING OF LARGE ARRAYS $ C C*********************************************************************** C* * C* THE FILE 'T2' IS INCLUDED * C* * C*********************************************************************** C C INCLUDE 'T2' C COMMON/MADIM/M1,M2,M3,M4,M5,M6,M7,M8 C C======================================================================= C C C$$$$$$$$$ COMMON BLOCKS FOR ELEMENTS $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C THESE BLOCKS HAVE A LENGTH OF NEL (= NUMBER OF ELEMENTS) C COMMON/E1/ELEM(MNEL) COMMON/E2/MATX(MNEL) COMMON/E3/EVOL(MNEL) COMMON/E4/PHI(MNEL) COMMON/E5/P(MNEL) COMMON/E6/T(MNEL) common/e7/pm(mnel) COMMON/VINWES/AI(MNEL) COMMON/AHTRAN/AHT(MNEL) c.....7-20-93: define coordinate arrays common/xyz1/x1(mnel) common/xyz2/x2(mnel) common/xyz3/x3(mnel) C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C$$$$$$$$$ COMMON BLOCKS FOR PRIMARY VARIABLES $$$$$$$$$$$$$$$$$$$$$$$$$ C C DATA BLOCKS /P1/, /P2/, /P3/ HAVE A LENGTH OF (NK+1)*NEL. C NK IS THE NUMBER OF MASS BALANCE EQUATIONS PER GRID BLOCK. C DATA BLOCKS /P4/ THROUGH /P7/ HAVE A LENGTH OF NEQ*NEL. C NEQ IS THE NUMBER OF BALANCE EQUATIONS PER GRID BLOCK. C COMMON/P1/X((MNK+1)*MNEL) COMMON/P2/DX((MNK+1)*MNEL) COMMON/P3/DELX((MNK+1)*MNEL) COMMON/P4/R(MNEQ*MNEL+1) COMMON/P5/DOLD(MNEQ*MNEL) COMMON/P6/ROW(MNEQ*MNEL) COMMON/P7/COL(MNEQ*MNEL) C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C C$$$$$$$$$ COMMON BLOCKS FOR CONNECTIONS $$$$$$$$$$$$$$$$$$$$$$$$$$$C C C THESE BLOCKS HAVE A LENGTH OF NCON (=NUMBER OF CONNECTIONS). C COMMON/C1/NEX1(MNCON) COMMON/C2/NEX2(MNCON) COMMON/C3/DEL1(MNCON) COMMON/C4/DEL2(MNCON) COMMON/C5/AREA(MNCON) COMMON/C6/BETA(MNCON) COMMON/C7/ISOX(MNCON) COMMON/C8/GLO(MNCON) COMMON/C9/ELEM1(MNCON) COMMON/C10/ELEM2(MNCON) COMMON/C11/FVD(MNCON) c.....6-8-95: new array SIG(MNCON) to hold coefficients for c radiative heat transfer. common/c12/sig(mncon) c.....6-4-93: append two common blocks used in the T2VOC version c of TOUGH2. common/flovp1/flovg(mnph*mncon) common/flovp2/flovw(mnph*mncon) c.....diffusive fluxes COMMON/FMOLDIF/FDIF(MNCON*MNPH*MNK) COMMON/VOLBC/VBC(MNEL) c.....6-09-99: COMMON blocks for T2DM C*** BEGIN ADDITION FOR T2DM C COMMON/NUMGBXYZ/NGB(3) COMMON/NDUPPAR/NZMULTI,NZDISF,NZDUP COMMON/DARVELM/DVELM((2*MNEQ+1)*MNPH*MNCON) COMMON/FMECDIS/FDIS(MNCON*MNPH*MNK) C PARAMETER (NDIMM = 2) PARAMETER (NNZERO = (4*NDIMM+1)*MNEQ*NREDM) C COMMON/SD1/N2Z(NNZERO) COMMON/SD2/IDLI(NNZERO) COMMON/SD3/IDLII(NNZERO) COMMON/SD4/IDLIII(NNZERO) COMMON/SD5/IFLAGG(NNZERO) COMMON/SD6/NFLAGG(NNZERO) COMMON/SD7/N3ZMRK,ID C C*** END ADDITION FOR T2DM C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C$$$$$$$$$ COMMON BLOCKS FOR LINEAR EQUATIONS $$$$$$$$$$$$$$$$$$$$$$$$$ C c arrays used by both MA28 and the conjugate gradient package. COMMON/L1/IRN(mnz) COMMON/L2/ICN(mnz) COMMON/L3/CO(mnz+1) COMMON/L4/WKAREA(MNEQ*MNEL+10) C WKAREA HAS A LENGTH OF NEQ*NEL+10. C c JVECT is given a slightly larger dimension than needed for c MA28, to accommodate the conjugate gradients. COMMON/L7/JVECT(niwork) C JVECT HAS A LENGTH OF NZ. c c arrays used by MA28 only. COMMON/L5/IKEEP(5*MNEQ*MNEL) C IKEEP HAS A LENGTH OF 5*NEQ*NEL. COMMON/L6/IW(8*MNEQ*MNEL) C IW HAS A LENGTH OF 8*NEQ*NEL. COMMON/L8/W19A(5*MNEQ*MNEL) C W19A HAS A LENGTH OF 5*NEQ*NEL. COMMON/AMMIS/MA,IPIV,U,IAB,NZ C c array used by conjugate gradient solvers only. c COMMON/CGARA6/RWORK(NRWORK) c COMMON/SOLVER/MATSLV,NMAXIT,ICLOSR,CLOSUR,ISYM,IUNIT,NVECTR,seed common/soll/lenw,leniw C C arrays used by luband only. COMMON/lub1/AB(nrwork) COMMON/lub3/NSUPDI,NSUBDI,mnzp1,mnetp1,mnelp1,nnnbig COMMON/lub4/matord,nsubdg,nsupdg,ntotd C C********* END modification. C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C C$$$$$ COMMON BLOCKS FOR SINKS/SOURCES $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C COMMON/G1/F1(MGTAB) COMMON/G2/F2(MGTAB) COMMON/G3/F3(MGTAB) common/g3a/pw(mgtab) C ARRAYS F1, F2, AND F3 HOLD, RESPECTIVELY, TABLES OF TIMES, RATES, C AND ENTHALPIES OF TIME-DEPENDENT SINKS OR SOURCES. THEIR DIMENSIONS C MUST NOT BE SMALLER THAN THE TOTAL NUMBER OF SUCH DATA. C COMMON/G4/ELEG(MNOGN) COMMON/G5/SOURCE(MNOGN) COMMON/G6/LTABG(MNOGN) COMMON/G7/G(MNOGN) COMMON/G8/EG(MNOGN) COMMON/G9/NEXG(MNOGN) COMMON/G10/ITABG(MNOGN) COMMON/G11/NGIND(MNOGN) COMMON/G12/LCOM(MNOGN) COMMON/G13/PI(MNOGN) COMMON/G14/PWB(MNOGN) COMMON/G15/HG(MNOGN) COMMON/G16/GPO(MNOGN) COMMON/G17/SDENS(MNOGN) COMMON/G18/SSAT(MNOGN) COMMON/G19/GVOL(MNOGN) COMMON/G20/HL(MNOGN) COMMON/G21/HS(MNOGN) COMMON/G22/QVGC(MNOGN) COMMON/G23/QVWC(MNOGN) COMMON/G24/QVOC(MNOGN) COMMON/G25/GRAD(MNOGN) C THE DIMENSION OF ALL ARRAYS IN COMMON G4 THROUGH G25 MUST BE C EQUAL TO OR LARGER THAN THE TOTAL NUMBER OF SINKS AND SOURCES. C COMMON/G26/FF(MNPH*MNOGN) C THE DIMENSION OF ARRAY FF MUST BE EQUAL TO OR LARGER THAN THE C PRODUCT OF NUMBER OF PHASES AND TOTAL NUMBER OF SINKS AND SOURCES. c common/g27/fnam(mnogn) common/g28/nftab(mnogn) common/g29/iftit(mnogn) common/g30/jftit(mnogn) common/g31/ijf(mnogn) C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C COMMON/SECPAR/PAR((MNEQ+1)*(MNPH*(MNB+MNK)+2)*MNEL) C ARRAY PAR HAS A LENGTH OF (NEQ+1)*NSEC*NEL. C NSEC = NPH*(NB+NK)+2 IS THE NUMBER OF SECONDARY PARAMETERS. C NPH IS THE NUMBER OF PHASES. C NB IS THE NUMBER OF THERMOPHYSICAL PARAMETERS (USUALLY 6). C NK IS THE NUMBER OF COMPONENTS (SPECIES). C COMMON/COMPO/FLO(MNPH*MNCON) COMMON/PORVEL/VEL(MNPH*MNCON) C COMMON/TITLE/TITLE CHARACTER*80 TITLE CHARACTER*5 ELEM,ELEM1,ELEM2,ELEG,SOURCE,VV LOGICAL EX COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/DMN/INUM,IPRINT,MCYC,MCYPR,MSEC,TZERO,TIMP1 COMMON/SVZ/NOITE,MOP(24) COMMON/LDIM/LICN,LIRN COMMON/BC/NELA COMMON/V/IV,IS common/ran/iran common/ff/h1 character*1 h1 COMMON/SOLVR1/matslv,nmaxit,nnvvcc,iiuunn,iissoo,nactdi COMMON/SOLVR2/ritmax,closur COMMON/SOLVR3/ordrng,oprocs,zprocs,coord CHARACTER*2 ordrng,oprocs,zprocs CHARACTER*5 coord C C DIMENSION VV(26) C CRAY-SPECIFIC: CONNECT TO FILE "INPUT" AS DEFAULT INPUT, FILE "OUTPUT" C AS DEFAULT OUTPUT. C CALL LINK('UNIT5=(INPUT,OPEN,TEXT),READ5, C X UNIT6=(OUTPUT,CREATE,HC),PRINT6 //') C C mnelp1 = mnel+1 mnzp1 = mnz+1 mnetp1 = nredm+1 lenw=nrwork leniw=niwork C h1=char(12) PRINT 1,h1 1 FORMAT(A1/7X,'@@@@@ @@ @ @ @@@ @ @ @@ @@@ @ ', A ' @ @ ', 1 '@ @ @ @@ @@@@@ @ @@ @ @ @@@ @ @ @ @'/ 2 7X,' @ @ @ @ @ @ @ @ @ @ @ @ @@ @@', B ' @ @', 3 ' @ @ @ @ @ @ @ @@ @ @ @ @ @ @@ @'/ 4 7X,' @ @ @ @ @ @ @@ @@@@ @ @@ @ @ @ @', C ' @ @', 5 ' @ @@@@ @ @ @ @ @ @ @ @@@ @ @ @ @ @'/ 6 7X,' @ @ @ @ @ @ @ @ @ @ @ @ @ @', D ' @ @', 71X, ' @ @ @ @ @ @ @ @ @@ @ @ @ @ @ @@'/ 8 7X,' @ @@ @@ @@@ @ @ @@@@ @@@ @ @ @', E ' @@ ', 9 ' @@@@ @ @ @ @ @@ @ @ @ @ @@ @ @'//) C PRINT 4 4 FORMAT(20X,'TOUGH2 IS A PROGRAM FOR MULTIPHASE MULTICOMPONENT', X' FLOW IN PERMEABLE MEDIA, INCLUDING HEAT FLOW.'/17X, X' IT IS A MEMBER OF THE MULKOM FAMILY OF CODES, DEVELOPED', X' AT LAWRENCE BERKELEY NATIONAL LABORATORY.') PRINT 5 5 FORMAT(/26X,80('*')/26X,19('*'),41X,20('*')/ X26X,19('*'),' TOUGH2 - VERSION 2.0 (OCTOBER 1999) ',20('*')/ X26X,19('*'),' T2CG2 Solver Package ',20('*')/ X26X,19('*'),41X,20('*')/26X,80('*')/) print 23 23 format(8X,'Copyright 1999 by The Regents of the University of ', x'California (subject to approval by the U.S. Department of ', x'Energy).') print 24 24 format(/ x' NOTICE: This software was developed under funding from the ', x'U.S. Department of Energy and the U.S. Government consequently ', x'retains'/' certain rights as follows: the U.S. Government has ', x'been granted for itself and others acting on its behalf a ', x'paid-up, nonexclusive,'/' irrevocable, worldwide license in ', x'the software to reproduce, prepare derivative works, and ', x'perform publicly and display publicly.'/' Beginning five (5) ', x'years after the date permission to assert copyright is obtained', x' from the U.S. Department of Energy, and subject'/' to any ', x'subsequent five (5) year renewals, the U.S. Government is ', x'granted for itself and others acting on its behalf a paid-up,'/ x' nonexclusive, irrevocable, worldwide license in the software ', x'to reproduce, prepare derivative works, distribute copies to ', x'the'/' public, perform publicly and display publicly, and to ', x'permit others to do so.') c print 22 22 format( x5X,'To engage radiative heat transfer, enter a coefficient "S', x'IGX" in the last field (columns 71-80) of a connection record.'/ x5X,'SIGX is the (relative) "radiant emittance".'/ x5X,'Entering SIGX < 0 will cause conductive heat transfer to be ', x'suppressed; will use abs(SIGX) as radiative transfer ', x'coefficient.'/ x5X,'Internally, SIGX will be multiplied with the Stefan-', x'Boltzmann constant, which is 5.669e-8 J/(m^2 K^4 s).') C M1=MNEL M2=MNCON M3=MNEQ M4=MNPH M5=MNB+MNK M6=MNOGN M7=MGTAB M8=MNK c revise length assignment of MA28 arrays LIRN=MNZ LICN=MNZ C MPRIM=(MNK+1)*MNEL MSEC=(MNEQ+1)*(MNPH*(MNB+MNK)+2)*MNEL PRINT 6,MNEL,MNCON,MNEQ,MNK,MNPH,MNB,MNOGN,MGTAB 6 FORMAT(//' PARAMETERS FOR FLEXIBLE DIMENSIONING OF MAJOR ARRAYS ', X'(MAIN PROGRAM) ARE AS FOLLOWS'// X' MNEL =',I7,' MNCON =',I7,' MNEQ =',I3,' MNK =',I3, X' MNPH =',I3, X' MNB =',I3,' MNOGN =',I5,' MGTAB =',I6/' ',131('=')) C PRINT 7,MNEL,MNCON,MPRIM,MNOGN 7 FORMAT(/' MAXIMUM NUMBER OF VOLUME ELEMENTS (GRID BLOCKS):',12X, X'MNEL =',I8/ X' MAXIMUM NUMBER OF CONNECTIONS (INTERFACES):',17X,'MNCON =',I8/ X' MAXIMUM LENGTH OF PRIMARY VARIABLE ARRAYS:',18X,'MPRIM =',I8/ X' MAXIMUM NUMBER OF GENERATION ITEMS (SINKS/SOURCES):',9X, X'MNOGN =',I8) PRINT 9,MGTAB,MSEC,MNZ,LIRN,LICN 9 FORMAT( X' MAXIMUM NUMBER OF TABULAR (TIME-DEPENDENT) GENERATION DATA', X': MGTAB =',I8/ X' LENGTH OF SECONDARY PARAMETER ARRAY:',24X,'MSEC =',I8/ X' MAXIMUM NUMBER OF JACOBIAN MATRIX ELEMENTS:',17X,'MNZ =',I8/ X/' LARGE LINEAR EQUATION ARRAYS: LENGTH OF IRN IS',14X, X'LIRN =',I8/ X' LENGTH OF ICN AND CO IS',7X, X'LICN =',I8//' ',131('=')) C print 21 21 format(/' array dimensioning is made according to the needs of', x' the conjugate gradient solvers'/ x' when using MA28 or LUBAND, only a smaller-size problem can be', x' accommodated'/' restriction with MA28 is: ', x'{number of elements} + 2 * {number of connections} < {MNEL', x' + 2* MNCON}/4'//' ',131('=')) c---*----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8 c CALL IO CALL SECOND(TZERO) C WRITE(11,899) c 8 FORMAT(/' TOUGH2 1.11 CG 28 January 1994',6X, c 8 FORMAT(/' TOUGH2 1.12 CG 12 February 1997',6X, c 8 FORMAT(/' TOUGH2 1.20 CG 26 September 1997',6X, c 899 FORMAT(/' TOUGH2 1.20 CG 12 February 1998',6X, c 899 FORMAT(/' TOUGH2 1.20 CG 11 March 1998',6X, c 899 FORMAT(/' TOUGH2 2.00 19 May 1999',6X, 899 FORMAT(/' TOUGH2 2.00 8 June 1999',6X, X'MAIN PROGRAM'/ x47X,'special version for conjugate gradient package T2CG2'/ x47X,'includes definition of coordinate arrays', x' and radiative heat transfer capability') C READ 2,TITLE 2 FORMAT(A80) PRINT 3,TITLE 3 FORMAT(/1X,131('=')//5X,'PROBLEM TITLE: ',A80//) C C-----FETCH INPUT DATA. C CALL INPUT C EVALUATE FLOATING POINT PROCESSOR. CALL FLOP C IF(IS.NE.0) GOTO 100 C CALL RFILE PRINT 10,NEL,NELA,NCON,NOGN 10 FORMAT(/' MESH HAS',I7,' ELEMENTS (',I7,' ACTIVE) AND',I7, A' CONNECTIONS (INTERFACES)', A' BETWEEN THEM'/' GENER HAS',I6,' SINKS/SOURCES') if(iran.ne.0) print 112 112 format(' BLOCK-BY-BLOCK PERMEABILITY MODIFICATION IS IN EFFECT') C c c.....set solver type and make informative printout call sinsub c IF(MOP(7).NE.0) CALL INDATA C CALL SECOND(ELT1) ELT1=ELT1-TZERO PRINT 11,ELT1 11 FORMAT(//' END OF TOUGH2 INPUT JOB --- ELAPSED TIME = ',F8.4,' SEC AONDS'/) C PRINT 79,h1 79 FORMAT(A1/' ',131('*')) PRINT 80 80 FORMAT(' * ARRAY *MOP*', X' ALLOWS TO GENERATE MORE PRINTOUT IN VARIOUS SUBROUTINES, AND', X' TO MAKE SOME CALCULATIONAL CHOICES. *') PRINT 78 78 FORMAT(' ',131('*')/) PRINT 81,MOP(1) 81 FORMAT(' ', X' MOP(1) =',I2,' *** ALLOWS TO GENERATE A SHORT PRINTOUT FOR' X,' EACH NEWTON-RAPHSON ITERATION'/ X' = 0, 1, OR 2: GENERATE 0, 1, OR 2 LINES OF PRINTOUT'// X12X,'MORE PRINTOUT IS GENERATED FOR MOP(I) > 0 IN THE FOLLOWING', X' SUBROUTINES (THE LARGER MOP IS, THE MORE WILL BE PRINTED).'/) PRINT 82,MOP(2),MOP(3),MOP(4),MOP(5),MOP(6),mop(8) 82 FORMAT(' ', X' MOP(2) =',I2,' *** CYCIT ', X' MOP(3) =',I2,' *** MULTI ', X' MOP(4) =',I2,' *** QU ', X' MOP(5) =',I2,' *** EOS ', X' MOP(6) =',I2,' *** LINEQ '/ X' MOP(8) =',I2,' *** DISF (T2DM ONLY)'/) PRINT 83,MOP(7) 83 FORMAT(' ', X' MOP(7) =',I2,' *** IF UNEQUAL ZERO, WILL GENERATE A PRINTOUT', X' OF INPUT DATA'/) C PRINT 84,MOP(9),MOP(10) 84 FORMAT(' ', X8X,' CALCULATIONAL CHOICES OFFERED BY MOP ARE AS FOLLOWS:'// X' MOP(9) =',I2,' *** CHOOSES FLUID COMPOSITION ON WITHDRAWAL', X' (PRODUCTION).'/ X' = 0: ACCORDING TO RELATIVE MOBILITIES.'/ X' = 1: ACCORDING TO COMPOSITION IN PRODUCING ELEMENT.'// X' MOP(10) =',I2,' *** CHOOSES INTERPOLATION FORMULA FOR DEPEN', X'DENCE OF THERMAL CONDUCTIVITY ON LIQUID SATURATION (SL).'/ X' = 0: K = KDRY + SQRT(SL)*(KWET-KDRY)'/ X' = 1: K = KDRY + SL*(KWET-KDRY)'/) PRINT 85,MOP(11) 85 FORMAT(' ', X' MOP(11) =',I2,' *** CHOOSES EVALUATION OF MOBILITY', X' AND ABSOLUTE PERMEABILITY AT INTERFACES.'/ X' = 0: MOBILITIES ARE UPSTREAM WEIGHTED WITH WUP.', X' (DEFAULT IS WUP = 1.0). PERMEABILITY IS UPSTREAM WEIGHTED.'/ X' = 1: MOBILITIES ARE AVERAGED BETWEEN', X' ADJACENT ELEMENTS. PERMEABILITY IS UPSTREAM WEIGHTED.'/ X' = 2: MOBILITIES ARE UPSTREAM WEIGHTED WITH WUP.', X' (DEFAULT IS WUP = 1.0). PERMEABILITY IS HARMONIC WEIGHTED.'/ X' = 3: MOBILITIES ARE AVERAGED BETWEEN', X' ADJACENT ELEMENTS. PERMEABILITY IS HARMONIC WEIGHTED.'/ X' = 4: MOBILITY * PERMEABILITY PRODUCT IS HARMONIC', X' WEIGHTED.'/) PRINT 86,MOP(12),mop(13),MOP(14),MOP(15) 86 FORMAT( X' MOP(12) =',I2,' *** CHOOSES PROCEDURE FOR INTERPOLATING ', X'GENERATION RATES FROM A TIME TABLE.'/ X' = 0: TRIPLE LINEAR INTERPOLATION.'/ X' = 1: "STEP FUNCTION" OPTION.'// X' MOP(13) =',I2,' *** T2DM ONLY. SPECIFIES ASSIGNMENT OF', X' COMPONENTS OF BOUNDARY VELOCITY AND', X' CONCENTRATION GRADIENT VECTORS.'/ X' AFFECTS ONLY THE INTERPOLATED ', X' COMPONENTS. DIRECT COMPONENTS ARE USED WHERE AVAILABLE.'/ X' = 0: VELOCITY AND GRADIENT AT BOUNDARY ARE ZERO.'/ X' = 1: VELOCITY IS ZERO; GRADIENT AT BOUNDARY IS', X' NEAREST-NEIGHBOR.'/ X' = 2: VELOCITY IS NEAREST NEIGHBOR; GRADIENT AT', X' BOUNDARY IS ZERO.'/ X' = 3: VELOCITY AND GRADIENT AT BOUNDARY', X' ARE NEAREST-NEIGHBOR.'// X' MOP(14) =',I2,' *** SPECIFIES THE HANDLING OF PIVOT FAILURES', X' IN THE LINEAR EQUATION SOLUTION (MA28 only)'/ X' = 0: PERFORM NEW MATRIX DECOMPOSITION AFTER PIVOT', X' FAILURE'/ X' > 0: IGNORE PIVOT FAILURE AND PROCEED',// X' MOP(15) =',I2,' *** ALLOWS TO SELECT A SEMI-ANALYTICAL HEAT', X' EXCHANGE CALCULATION WITH CONFINING BEDS.'/ X' = 0: NO SEMI-ANALYTICAL HEAT EXCHANGE'/ X' > 0: SEMI-ANALYTICAL HEAT EXCHANGE ENGAGED (WHEN A', X' SPECIAL SUBROUTINE MODULE *QLOSS* IS PRESENT)'/) PRINT 87,MOP(16),MOP(17),MOP(18) 87 FORMAT(' ', X' MOP(16) =',I2,' *** PERMITS TO CHOOSE TIME STEP SELECTION ', X'OPTION'/ X' = 0: USE TIME STEPS EXPLICITLY PROVIDED AS INPUT.'/ X' > 0: INCREASE TIME STEP BY AT LEAST A FACTOR 2, IF', X' CONVERGENCE OCCURS IN .LE. MOP(16) ITERATIONS.'// X' MOP(17) =',I2,' *** HANDLES SCALING OPTIONS.'/ X' = 0: NO SCALING.'/ X' = 7: SCALING.'// X' MOP(18) =',I2,' *** ALLOWS TO SELECT HANDLING OF INTERFACE ', X'DENSITY.'/ X' = 0: PERFORM UPSTREAM WEIGHTING FOR INTERFACE DENSITY. X'/' > 0: COMPUTE INTERFACE DENSITY AS AVERAGE OF', X' THE TWO GRID BLOCK DENSITIES.'/ X' HOWEVER, WHEN ONE OF THE TWO PHASE SATURATIONS', X' IS ZERO, DO UPSTREAM WEIGHTING.'/) c print 88,mop(19),mop(20),mop(21) 88 format(' MOP(19) =',I2,' *** is used in some EOS-modules for', X' selecting different options'// x' MOP(20) =',I2,' *** is used in some EOS-modules for', X' selecting different options'// x' MOP(21) =',I2,' *** PERMITS TO SELECT LINEAR', x' EQUATION SOLVER'/ X' = 0: DEFAULTS TO MOP(21) = 3'/ x' = 1: DIRECT SOLVER - MA28'/ x' = 2: SUBROUTINE DSLUBC: BI-CONJUGATE GRADIENT', x' SOLVER; PRECONDITIONER: INCOMPLETE LU FACTORIZATION'/ x' = 3: SUBROUTINE DSLUCS: BI-CONJUGATE GRADIENT', x' SOLVER - LANCZOS TYPE; PRECONDITIONER: INCOMPLETE LU', x' FACTORIZATION'/ x' = 4: SUBROUTINE DSLUGM: GENERALIZED MINIMUM', x' RESIDUAL CONJUGATE GRADIENTS; PRECONDITIONER: INCOMPLETE LU', x' FACTORIZATION'/ x' = 5: SUBROUTINE DLUSTB: Stabilized bi-conjugate', x' gradient solver; PRECONDITIONER: INCOMPLETE LU', x' FACTORIZATION'/ x' = 6: SUBROUTINE LUBAND: Direct solver using', x' LU decomposition'//) c PRINT 89,MOP(22),MOP(23) 89 FORMAT(' MOP(22) =',I2,' *** T2DM ONLY. SPECIFIES METHOD OF', X' SUMMATION OF DUPLICATE ELEMENTS IN THE JACOBIAN. ***',/, X' = 0: USE SUMDUP2.',/, X' = 1: USE SUMDUP.',/, X' = 2: IF MATSLV = 1, USE MA28 INTERNAL SUMMATION.'//, X' MOP(23) =',I2,' *** T2DM ONLY. HANDLES EFFECTS OF', X' NON-CONNECTED NEIGHBOR GRID BLOCKS ON DISPERSIVE FLUX IN ', X' DISF. ***',/, X' = 0: INCLUDE INFLUENCE OF NEIGHBOR GRID BLOCKS.', X' (MORE ACCURATE JACOBIAN, SLOWER LINEAR EQUATION SOLUTION).'/, X' = 1: NEGLECT INFLUENCE OF NEIGHBOR GRID BLOCKS.', X' (LESS ACCURATE JACOBIAN, FASTER LINEAR EQUATION SOLUTION).'//) c print 90,mop(24) 90 format(' MOP(24) =',I2,' *** PERMITS TO SELECT HANDLING OF', X' MULTIPHASE DIFFUSIVE FLUXES AT INTERFACES'/ X' = 0: HARMONIC WEIGHTING OF FULLY-COUPLED EFFECTIVE ', X'MULTIPHASE DIFFUSIVITY'/ x' = 1: SEPARATE HARMONIC WEIGHTING OF GAS AND LIQUID ', X'PHASE DIFFUSIVITIES.'// X' ',131('*')) C C-----SOLVE MASS- AND ENERGY-TRANSPORT EQUATIONS. IF(IS.EQ.0) CALL CYCIT C IF IS.NE.0 (KEYWORD "ENDFI" IN INPUT FILE, OR PROBLEM IN RFILE), C NO FLOW SIMULATION WILL BE PERFORMED. C 100 CONTINUE IF(IV.EQ.0) GOTO 15 PRINT 20,h1 20 FORMAT(A1/ X' ',131('*')/' *',129X,'*'/' *',49X,'SUMMARY OF PROGRAM UNITS', X' USED',51X,'*'/' *',129X,'*'/' ',131('*')// X' UNIT VERSION DATE COMMENTS'/ X1X,131('_')/) C ENDFILE 11 REWIND 11 C 16 CONTINUE DO17 I=1,26 17 VV(I)=' ' C READ(11,14,END=15) (VV(I),I=1,26) 14 FORMAT(26A5) PRINT 14,(VV(I),I=1,26) GOTO 16 C 15 CONTINUE PRINT 19 19 FORMAT(//1X,131('*')) CALL SECOND(ELT) ELT=ELT-TZERO ELTC=ELT-ELT1 PRINT 12,ELT,ELTC,ELT1,h1 12 FORMAT(//' END OF TOUGH2 SIMULATION RUN --- ELAPSED TIME = ', AF9.3,' SEC-- CALCULATION TIME = ',F9.3,' SEC-- DATA INPUT TIME = ' B,F8.3,' SEC'/A1) STOP END C C C C********************************************************************* C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C********************************************************************* C C C BLOCK DATA VV COMMON/V/IV,IS DATA IV/1/ END C C C C********************************************************************* C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C********************************************************************* C C C SUBROUTINE IO COMMON/SVZ/NOITE,MOP(24) LOGICAL EX C SAVE ICALL DATA ICALL/0/ ICALL=ICALL+1 PRINT 1 1 FORMAT(//' SUMMARY OF DISK FILES'/) C INQUIRE(FILE='VERS',EXIST=EX) IF(EX) GOTO 62 PRINT 63 63 FORMAT(' FILE *VERS* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(11,FILE='VERS',STATUS='NEW') ENDFILE 11 GOTO 70 C 62 PRINT 64 64 FORMAT(' FILE *VERS* EXISTS --- OPEN AS AN OLD FILE') OPEN(11,FILE='VERS',STATUS='OLD') C 70 CONTINUE REWIND 11 C IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'IO 1.0 15 APRIL 1991',6X,'OPEN', X' FILES *VERS*, *MESH*, *INCON*, *GENER*, *SAVE*, *LINEQ*,', X' AND *TABLE*') C INQUIRE(FILE='MESH',EXIST=EX) IF(EX) GOTO 2 PRINT 3 3 FORMAT(' FILE *MESH* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(4,FILE='MESH',STATUS='NEW') GOTO 10 C 2 PRINT 4 4 FORMAT(' FILE *MESH* EXISTS --- OPEN AS AN OLD FILE') OPEN(4,FILE='MESH',STATUS='OLD') C 10 INQUIRE(FILE='INCON',EXIST=EX) IF(EX) GOTO 12 PRINT 13 13 FORMAT(' FILE *INCON* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(1,FILE='INCON',STATUS='NEW') ENDFILE 1 GOTO 20 C 12 PRINT 14 14 FORMAT(' FILE *INCON* EXISTS --- OPEN AS AN OLD FILE') OPEN(1,FILE='INCON',STATUS='OLD') C 20 INQUIRE(FILE='GENER',EXIST=EX) IF(EX) GOTO 22 PRINT 23 23 FORMAT(' FILE *GENER* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(3,FILE='GENER',STATUS='NEW') ENDFILE 3 GOTO 30 C 22 PRINT 24 24 FORMAT(' FILE *GENER* EXISTS --- OPEN AS AN OLD FILE') OPEN(3,FILE='GENER',STATUS='OLD') C 30 INQUIRE(FILE='SAVE',EXIST=EX) IF(EX) GOTO 32 PRINT 33 33 FORMAT(' FILE *SAVE* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(7,FILE='SAVE',STATUS='NEW') GOTO 40 C 32 PRINT 34 34 FORMAT(' FILE *SAVE* EXISTS --- OPEN AS AN OLD FILE') OPEN(7,FILE='SAVE',STATUS='OLD') C 40 INQUIRE(FILE='LINEQ',EXIST=EX) IF(EX) GOTO 42 PRINT 43 43 FORMAT(' FILE *LINEQ* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(15,FILE='LINEQ',STATUS='NEW') GOTO 50 C 42 PRINT 44 44 FORMAT(' FILE *LINEQ* EXISTS --- OPEN AS AN OLD FILE') OPEN(15,FILE='LINEQ',STATUS='OLD') REWIND 15 C 50 CONTINUE INQUIRE(FILE='TABLE',EXIST=EX) IF(EX) GOTO 52 PRINT 53 53 FORMAT(' FILE *TABLE* DOES NOT EXIST --- OPEN AS A NEW FILE') OPEN(8,FILE='TABLE',STATUS='NEW') ENDFILE 8 GOTO 60 C 52 PRINT 54 54 FORMAT(' FILE *TABLE* EXISTS --- OPEN AS AN OLD FILE') OPEN(8,FILE='TABLE',STATUS='OLD') C 60 RETURN END C C C C********************************************************************* C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C********************************************************************* C C C SUBROUTINE FLOP C C-----CALCULATE NUMBER OF SIGNIFICANT DIGITS FOR FLOATING POINT C PROCESSING. ASSIGN DEFAULT FOR DFAC, AND PRINT APPROPRIATE C WARNING WHEN MACHINE ACCURACY IS INSUFFICIENT. C ! IMPLICIT REAL*8(A-H,O-Z) C COMMON/CONTST/RE1,RE2,RERM,NER,KER,DFAC SAVE ICALL DATA ICALL/0/ ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'FLOP 1.0 11 APRIL 1991',6X, X'CALCULATE NUMBER OF SIGNIFICANT DIGITS FOR FLOATING POINT', X' ARITHMETIC') C A=SQRT(.99) B=A DO1 N=1,260 B=B/2. C=A+B D=C-A IF(D.EQ.0.) GOTO 2 1 CONTINUE C 2 B2=B*2. N10=-INT(LOG10(B2)) DF=SQRT(B2) C PRINT 3,N10,DF 3 FORMAT(/24X,84('*')/24X,'*',23X,'EVALUATE FLOATING POINT', X' ARITHMETIC',25X,'*'/24X,84('*')/24X,'*',82X,'*'/24X,'*', X' FLOATING POINT PROCESSOR HAS APPROXIMATELY',I3, X' SIGNIFICANT DIGITS',17X,'*'/24X,'*',82X,'*'/24X,'*', X' DEFAULT VALUE OF INCREMENT FACTOR FOR NUMERICAL DERIVATIVES', X' IS DFAC = ',E10.4,' *') C IF(DFAC.EQ.0.) THEN DFAC=DF PRINT 10 10 FORMAT(24X,'* DEFAULT VALUE FOR DFAC WILL BE USED',46X,'*') ELSE PRINT 11,DFAC 11 FORMAT(24X,'*', X' USER-SPECIFIED VALUE DFAC = ',E10.4,' WILL BE USED',30X,'*') ENDIF PRINT 6 6 FORMAT(24X,'*',82X,'*'/24X,84('*')/) C IF(N10.LE.12.AND.N10.GT.8) PRINT 4 4 FORMAT(' WWWWWWWWWW WARNING WWWWWWWWWW: NUMBER OF SIGNIFICANT', X' DIGITS IS MARGINAL; EXPECT DETERIORATED CONVERGENCE BEHAVIOR'/) IF(N10.LE.8) PRINT 5 5 FORMAT(' WWWWWWWWWW WARNING WWWWWWWWWW: NUMBER OF SIGNIFICANT', X' DIGITS IS INSUFFICIENT; CONVERGENCE WILL BE POOR OR FAIL'/ X' WWWWWWWWWWWWWWWWWWWWWWWWWWWWW: CODE SHOULD BE RUN IN DOUBLE' X,' PRECISION!'/) C RETURN END C C C C********************************************************************* C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C********************************************************************* C C C subroutine second(time) c timing function for IBM RS/6000 C ! IMPLICIT REAL*8(A-H,O-Z) C save icall data icall/0/ icall=icall+1 if(icall.eq.1) write(11,899) 899 FORMAT(6X,'SECOND 1.0 6 September 1994',6X, x'CPU timing function for IBM RS/6000') time=float(mclock())/100. c return end C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C subroutine sinsub C ! IMPLICIT REAL*8(A-H,O-Z) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/SVZ/NOITE,MOP(24) COMMON/BC/NELA common/ff/h1 character*1 h1 COMMON/SOLVR1/matslv,nmaxit,nnvvcc,iiuunn,iissoo,nactdi COMMON/SOLVR2/ritmax,closur COMMON/SOLVR3/ordrng,oprocs,zprocs,coord CHARACTER*2 ordrng,oprocs,zprocs CHARACTER*5 coord c SAVE ICALL DATA ICALL/0/ C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) c 899 FORMAT(6X,'SIN 1.00 26 May 1999',6X, 899 FORMAT(6X,'SINSUB 1.00 1 October 1999',6X, x'initialize parameters for the solver package, and generate', x' informative printout') C---*----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8 C print 1,h1 1 format(A1) c.....imported from INPUT of t2fa.f (GM) C C*********************************************************************** C* * C* SETTING SOLVER TYPE AND CORRESPONDING PARAMETERS * C* * C*********************************************************************** C C ordrng = 'ST' IF(iissoo.eq.0) THEN matslv = mop(21) zprocs = 'Z1' if(matslv.eq.1.or.matslv.eq.6) zprocs='Z0' oprocs = 'O0' ritmax = 1.0e-1 closur = 1.0e-6 END IF C PRINT 6015 IF(iissoo.eq.0) THEN PRINT 6020 ELSE PRINT 6021 END IF C IF(matslv.eq.0.or.matslv.gt.6) THEN PRINT 6025, matslv matslv = 3 ENDIF PRINT 6026, matslv C IF(matslv.EQ.1.or.matslv.eq.6) GO TO 3085 IF(ritmax.EQ.0.0e0.or.ritmax.GT.1.0e0) ritmax = 1.0e-1 IF((abs(closur)).gt.1.0e-6) closur = 1.0e-6 IF((abs(closur)).lt.1.0e-12) closur = 1.0e-12 riter = NELA*neq*ritmax nmaxit = MAX(20,INT(riter)) PRINT 6030, ritmax,closur,nmaxit 3085 continue c if(zprocs.ne.'Z0') then if(matslv.eq.1.or.matslv.eq.6) then print 6101, matslv,zprocs zprocs='Z0' endif endif c if(oprocs.ne.'O0') then if(matslv.eq.1.or.matslv.eq.6) then print 6036, matslv,oprocs oprocs='O0' endif endif c IF(zprocs.NE.'Z0'.AND.zprocs.NE.'Z1'.AND. & zprocs.NE.'Z2'.AND.zprocs.NE.'Z3'.AND. & zprocs.NE.'Z4') THEN print 6102, zprocs zprocs = 'Z1' END IF c IF(oprocs.NE.'O0'.AND.oprocs.NE.'O1'.AND. & oprocs.NE.'O2'.AND.oprocs.NE.'O3'.AND. & oprocs.NE.'O4') THEN print 6037, oprocs oprocs = 'O0' END IF c if(neq.eq.1) then c.....no preprocessing for NEQ=1 print 6200 zprocs='Z0' oprocs='O0' endif c print 6103, zprocs print 6038, oprocs C 6015 FORMAT(130('*'),/,'*',T130,'*',/,'*', & T29,'M A T R I X S O L V E R A N D ', & 'R E L E V A N T I N F O R M A T I O N',T130,'*', & /,'*',T130,'*',/,130('*'),/) 6020 FORMAT(T5,'The solver is determined from MOP(21)') 6021 FORMAT(T5,'The solver is determined from the SOLVR ', & 'data block - MOP(21) IS OVERRIDEN !!!',/) 6025 FORMAT(T5,'The solution method indicator MATSLV = ',i2, & 5x,'- reset internally to the default, MATSLV = 3') 6026 FORMAT(/T5,'The solution method indicator MATSLV = ',i2,/, & T10,'MATSLV = 1: SUBROUTINE MA28 - ', & 'Direct solver using LU decomposition ',/, & T10,'MATSLV = 2: SUBROUTINE DSLUBC - ', & 'BI-CONJUGATE GRADIENT SOLVER',/, & T42,'Incomplete LU factorization preconditioning',/, & T10,'MATSLV = 3: SUBROUTINE DSLUCS (DEFAULT) - ', & 'Lanczos-type Conjugate Gradient Squared solver ',/, & T52,'Incomplete LU factorization preconditioning',/, & T10,'MATSLV = 4: SUBROUTINE DSLUGM - ', & 'Generalized Minimum Residual Conjugate Gradient ', & 'solver',/, & T42,'Incomplete LU factorization preconditioning',/, & T10,'MATSLV = 5: SUBROUTINE DLUSTB - ', & 'STABILIZED BI-CONJUGATE GRADIENT SOLVER ',/, & T42,'Incomplete LU factorization preconditioning',/ & T10,'MATSLV = 6: SUBROUTINE LUBAND - ', & 'Direct solver using LU decomposition ',/) 6030 FORMAT(T5,'RITMAX: Maximum # of CG iterations as fraction of ', & 'the total number of equations ',T95,'= ',1PE12.5,/, & T10,'(0.0 < RITMAX <= 1.0, Default = 0.1)',/, & T5,'CLOSUR: Convergence criterion for the CG ', & 'iterations ',T95,'= ',1PE12.5,/, & T10,'(1.0e-12 <= CLOSUR <= 1.0e-6, Default = 1.0e-6)',/, & T5,'NMAXIT: Maximum # of CG iterations - not to exceed ', & 'the total number of equations NELA*NEQ', & T95,'= ',I5,/T10,'(20 < NMAXIT <= NREDM)',/) C---*----1----*----2----*----3----*----4----*----5----*----6----*----7----*----8 6101 FORMAT(/,' ',15('WARNING-'),//T14,'For MATSLV = ',I1,', no Z-', &'preprocessing can be used.'/T14,'Action taken: reset ZPROCS = ', &A2,' to *Z0*; continue execution') 6036 FORMAT(/,' ',15('WARNING-'),//T14,'For MATSLV = ',I1,', no O-', &'preprocessing can be used.'/T14,'Action taken: reset OPROCS = ', &A2,' to *O0*; continue execution') 6102 FORMAT(/,' ',15('WARNING-'),//T14,'Unknown matrix preprocessing', x' option ZPROCS = ',a2/T14,'Action taken: reset ZPROCS to', x' *Z1*; continue execution') 6037 FORMAT(/,' ',15('WARNING-'),//T14,'Unknown matrix preprocessing', x' option OPROCS = ',a2/T14,'Action taken: reset OPROCS to', x' *O0*; continue execution') 6200 format(' !!!!! NEQ=1; do not perform any matrix preprocessing') 6103 FORMAT(/T5,'The matrix Z-preprocessing system is ZPROCS = ',a2,// & T10,'ZPROCS = Z0: No Z-preprocessing; ', & 'default for NEQ = 1 and for MATSLV = 1, 6'/ & T10,'ZPROCS = Z1: Replacement of zeros on ', & 'the main-diagonal by a small number; '/ & T24,'default for NEQ > 1 and for 1 < MATSLV < 6'/ & T10,'ZPROCS = Z2: Linear combination of equations ', & 'in each element to produce non-zero ', & 'main diagonal entries',/, & T10,'ZPROCS = Z3: Normalization of equations, ', & 'followed by Z2',/, & T10,'ZPROCS = Z4: Same as in OPROCS = O4',/) 6038 FORMAT(/T5,'The matrix O-preprocessing system is OPROCS = ',a2,// & T10,'OPROCS = O0: No O-preprocessing; ', & 'default for NEQ = 1 and for MATSLV = 1, 6'/ & T10,'OPROCS = O1: Elimination of lower half of ', & 'the main-diagonal submatrix with center ', & 'pivoting',/, & T10,'OPROCS = O2: O1+Elimination of upper half of ', & 'the main-diagonal submatrix with center ', & 'pivoting',/, & T10,'OPROCS = O3: O2+Normalization - Results in ', & 'unit main-diagonal submatrices ',/, & T10,'OPROCS = O4: pre-processing which results ', & 'in unit main-diagonal submatrices without center ', & 'pivoting',/) C c return end C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE LINEQ1 C C C ! IMPLICIT REAL*8(A-H,O-Z) C C-----THIS SUBROUTINE CALLS THE LINEAR EQUATION SOLVER MA28. C IT HAS LOGIC TO HANDLE FAILURES IN LINEAR EQUATION SOLUTION. C C AFTER SOLUTION, LATEST UPDATED ITERATES ARE OBTAINED FOR C ALL PRIMARY DEPENDENT VARIABLES. C C C C$$$$$$$$$ COMMON BLOCKS FOR LINEAR EQUATIONS $$$$$$$$$$$$$$$$$$$$$$$$$ C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) COMMON/L4/WKAREA(1) COMMON/L5/IKEEP(1) COMMON/L6/IW(1) COMMON/L7/JVECT(1) COMMON/L8/W19A(1) COMMON/AMMIS/MA,IPIV,U,IAB,NZ C C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C COMMON/E1/ELEM(1) COMMON/P2/DX(1) COMMON/P4/R(1) COMMON/P6/ROW(1) COMMON/P7/COL(1) COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/SVZ/NOITE,MOP(24) COMMON/DM/DELTEN,DELTEX,FOR,FORD COMMON/KONIT/KON,DELT,IGOOD COMMON/DG/WUP,WNR COMMON/MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK, AABORT1,ABORT2 COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT COMMON/LDIM/LICN,LIRN COMMON/BC/NELA CHARACTER*5 ELEM C SAVE ICALL,N DATA ICALL/0/ ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'LINEQ1 1.0 22 JANUARY 1990',6X, X'INTERFACE FOR THE LINEAR EQUATION SOLVER MA28') C ISCALE=0 IF(MOP(17).GE.7) ISCALE=1 IF(MOP(17).GE.5.AND.ITER.EQ.1) ISCALE=1 MA=MA+1 IF(MA.EQ.1) N=NEQ*NELA IF(ISCALE.EQ.0) GOTO 300 C 128 CONTINUE C-----COME HERE TO COMPUTE SCALING FACTORS. IF(IAB.EQ.0) CALL MC19A(N,NZ,CO,IRN,ICN,ROW,COL,W19A) IF(IAB.NE.0) CALL MC19A(N,NZ,CO,IRN,JVECT,ROW,COL,W19A) DO301 I=1,N ROW(I)=EXP(ROW(I)) COL(I)=EXP(COL(I)) 301 CONTINUE 300 CONTINUE C C C-----FOR IAB.EQ.0 CALL MA28A, OTHERWISE CALL MA28B. C INUM=0 IGOOD=0 C-----MA COUNTS CALLS TO LINEQ1. C IF(MOP(6).NE.0) PRINT 109,MA,N,NZ 109 FORMAT(/,' !!!! L I N E Q 1 !!!! !!!!! !!!!! MA =', & I4,' N =',I4,' NZ =',I6) C 117 CONTINUE C IF(ISCALE.EQ.0) GOTO 200 C C-----COME HERE TO SCALE MATRIX. DO 201 II=1,NZ I=IRN(II) IF(IAB.EQ.0) J=ICN(II) IF(IAB.NE.0) J=JVECT(II) 201 CO(II)=CO(II)*ROW(I)*COL(J) 200 CONTINUE C C IF(MOP(6).LT.7) GOTO99 PRINT 102 102 FORMAT(/23H MATRIX OF COEFFICIENTS/) IF(IAB.NE.0) PRINT 100,(IRN(NN),JVECT(NN),CO(NN),NN=1,NZ) IF(IAB.EQ.0) PRINT 100,(IRN(NN),ICN(NN),CO(NN),NN=1,NZ) 100 FORMAT(5(1X,2I5,E14.6)) 99 CONTINUE C C C IF(MOP(6).NE.0) CALL THYME(0,TS,TT) C IF(IAB.EQ.0) ACALL MA28A(N,NZ,CO,LICN,IRN,LIRN,ICN,U,IKEEP,IW,WKAREA,IFLAG) C C IF(IAB.NE.0) ACALL MA28B(N,NZ,CO,LICN,IRN,JVECT,ICN,IKEEP,IW,WKAREA,IFLAG) C IAB=1 C IF(IFLAG.EQ.-2) GOTO 2 INUM=0 C-----IF NO NUMERICAL SINGULARITY IS ENCOUNTERED, SET INUM=0. IF(IFLAG.LE.0) GOTO115 C C=====PIVOT FAILURE=====PIVOT FAILURE=====PIVOT FAILURE===== C IPIV=IPIV+1 C IF(MOP(17).EQ.2.OR. MOP(17).EQ.3) GOTO 137 GOTO 138 137 ISCALE=1 IF(MOP(6).NE.0) XPRINT 140,IFLAG,IPIV 140 FORMAT(21H PIVOT FAILURE IN ROW,I4,19H --- RESCALE MATRIX,13H == A= IPIV =,I5) GOTO 128 138 CONTINUE C IF(MOP(14).EQ.0) GOTO1 C-----COME HERE TO PROCEED AFTER PIVOT FAILURE WITHOUT MAKING NEW C DECOMPOSITION. IF(MOP(6).NE.0) XPRINT 118,IFLAG,IPIV 118 FORMAT(24H VERY SMALL PIVOT IN ROW,I4,13H === IPIV =,I5) C FOR MOP(14) > 0 IGNORE PIVOT PROBLEM AND PROCEED. GOTO 115 C C C-----COME HERE FOR PIVOT FAILURE IN MA28B, AND PROCEED TO CALL C MA28A TO PERFORM A NEW DECOMPOSITION. C 1 IF(MOP(6).NE.0) PRINT 116,IFLAG,IPIV 116 FORMAT(24H VERY SMALL PIVOT IN ROW,I4,18H --- PERFORM NEW,14H DE ACOMPOSITION,12H === IPIV =,I5) GOTO 3 C 2 CONTINUE IF(MOP(17).EQ.1.OR. MOP(17).EQ.3) GOTO 127 GOTO 20 127 ISCALE=1 IF(MOP(6).NE.0) PRINT 141 141 FORMAT(74H MATRIX NUMERICALLY SINGULAR --------------------------- A PERFORM SCALING) GOTO 128 C 20 PRINT 4 4 FORMAT(84H MATRIX NUMERICALLY SINGULAR --------------------------- A PERFORM NEW DECOMPOSITION) C INUM=INUM+1 C-----COUNT THE NUMBER OF SINGULARITIES THAT ARISE WITHOUT C INTERRUPTION BY A NON-SINGULAR DECOMPOSITION. C IAB=0 IF(INUM.EQ.2) GOTO 107 C-----IF TWO SUBSEQUENT DECOMPOSITIONS ENCOUNTER A SINGUL@RITY, C REDUCE THE TIME STEP. C 3 IAB=0 C C-----RECOMPUTE MATRIX---------- CALL MULTI C IF(IGOOD.EQ.0) GOTO117 PRINT 119 119 FORMAT(54H ++++++++++++++++++++ MULTI FAILS ++++++++++++++++++++) RETURN C 115 CONTINUE C C IF(MOP(6).EQ.0) GOTO110 CALL THYME(1,TD,TT) PRINT 111,TD,IFLAG 111 FORMAT(22H DECOMPOSITION TIME = ,E12.4,9H SECONDS,21H ********* A* IFLAG =,I4) 110 CONTINUE C C IF(IFLAG.LT.0.AND.IFLAG.GE.-12) GOTO107 C C IF(MOP(6).NE.0) CALL THYME(0,TS,TT) C IF(ISCALE.EQ.0) GOTO 210 C C-----COME HERE TO SCALE RIGHT HAND SIDE VECTOR. DO211 II=1,N 211 R(II)=R(II)*ROW(II) 210 CONTINUE C CALL MA28C(N,CO,LICN,ICN,IKEEP,R,WKAREA,1) C IF(ISCALE.EQ.0) GOTO 220 DO 221 II=1,N 221 R(II)=R(II)*COL(II) 220 CONTINUE C IF(MOP(6).EQ.0) GOTO112 CALL THYME(1,TSS,TT) PRINT 113,TSS 113 FORMAT(22H SOLUTION TIME = ,E12.4,9H SECONDS) 112 CONTINUE C IF(MOP(6).GE.2) PRINT 114,EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN, AMINICN,IRANK 114 FORMAT(7H EPS = ,E10.4,9H RMIN = ,E10.4,10H RESID = ,E10.4,9H I ARNCP =,I6,9H ICNCP =,I6,10H MINIRN =,I6,10H MINICN =,I6,9H IRA BNK =,I4) C GOTO103 C 107 IGOOD=1 IAB=0 108 CONTINUE C-----COME HERE FOR FAILURE IN SOLVING LINEAR EQUATION. PRINT 104,KCYC,ITER,IFLAG IF(IFLAG.GE.0) RETURN 104 FORMAT(' LINEQ1 FAILS AT (',I4,1H,,I3,') IFLAG =,',I4, x'; WILL REDUCE DELT') C IFL=-IFLAG GOTO(402,202,203,204,205,206,207,208,209,410,411,212),IFL 402 PRINT 401 RETURN 202 PRINT 302 RETURN 203 PRINT 303,MINIRN RETURN 204 PRINT 304 RETURN 205 PRINT 305,MINICN RETURN 206 PRINT 306,MINIRN,MINICN RETURN 207 PRINT 307 RETURN 208 PRINT 308,NZ RETURN 209 PRINT 309,NZ RETURN 410 PRINT 310,NZ RETURN 411 PRINT 311 RETURN 212 PRINT 312 RETURN 401 FORMAT(32H MATRIX IS STRUCTURALLY SINGULAR) 302 FORMAT(31H MATRIX IS NUMERICALLY SINGULAR) 303 FORMAT(45H INCREASE DIMENSION OF ARRAY IRN TO AT LEAST ,I6,9H ELEM AENTS) 304 FORMAT(51H DIMENSIONS OF ARRAYS ICN AND CO ARE MUCH TOO SMALL) 305 FORMAT(53H INCREASE DIMENSION OF ARRAYS ICN AND CO TO AT LEAST ,I7 A,9H ELEMENTS) 306 FORMAT(47H INSUFFICIENT DIMENSIONS; NEED AT LEAST IRN - ,I6,27H E ALEMENTS; ICN AND CO - ,I7,9H ELEMENTS) 307 FORMAT(33H ERROR IN BLOCK TRIANGULARIZATION) 308 FORMAT(28H MAKE LIRN LARGER THAN NZ = ,I7) 309 FORMAT(28H MAKE LICN LARGER THAN NZ = ,I7) 310 FORMAT(41H NUMBER OF NON-ZERO MATRIX ELEMENTS NZ = ,I7,25H MUST BE A LARGER THAN ZERO) 311 FORMAT(44H ORDER OF MATRIX MUST BE BETWEEN 1 AND 32767) 312 FORMAT(46H ROW OR COLUMN INDEX IS OUT OF RANGE, INCREASE,17H ARRAY A DIMENSIONS) C 103 CONTINUE C C-----UPDATE CHANGES. C IF(MOP(6).GE.5) PRINT 31 31 FORMAT(/' ===== INCREMENTS ==== MASS BALANCES FIRST,', X' ENERGY BALANCE LAST'/) DO5 NN=1,NELA NLOC=(NN-1)*NEQ NLOCP=(NN-1)*NK1 C IF(MOP(6).GE.5) PRINT 30,ELEM(NN),(R(NLOC+K),K=1,NEQ) 30 FORMAT(' AT ELEMENT *',A5,'* ',8(1X,E12.6)) C DO5 K=1,NEQ DX(NLOCP+K)=DX(NLOCP+K)+WNR*R(NLOC+K) 5 CONTINUE C C ISCALE=0 RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE LINEQ C C C ! IMPLICIT REAL*8(A-H,O-Z) C C-----THIS SUBROUTINE CALLS THE LINEAR EQUATION SOLVER T2CG2. C IT HAS LOGIC TO HANDLE FAILURES IN LINEAR EQUATION SOLUTION. C C AFTER SOLUTION, LATEST UPDATED ITERATES ARE OBTAINED FOR C ALL PRIMARY DEPENDENT VARIABLES. C C C$$$$$$$$$ COMMON BLOCKS FOR LINEAR EQUATIONS $$$$$$$$$$$$$$$$$$$$$$$$$ C PARAMETER(seed = 1.0e-25) PARAMETER(iunit = 0, nvectr = 30) C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) COMMON/L4/WKAREA(1) COMMON/L7/JVECT(1) C COMMON/AMMIS/MA,IPIV,U,IAB,NZ C common/soll/lenw,leniw C C arrays used by luband only. COMMON/lub1/AB(1) COMMON/lub3/NSUPDI,NSUBDI,mnzp1,mnetp1,mnelp1,nnnbig COMMON/lub4/matord,nsubdg,nsupdg,ntotd C COMMON/E1/ELEM(1) COMMON/P2/DX(1) COMMON/P4/R(1) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/SVZ/NOITE,MOP(24) COMMON/DM/DELTEN,DELTEX,FOR,FORD COMMON/KONIT/KON,DELT,IGOOD COMMON/DG/WUP,WNR COMMON/MA28F/EPS,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN,IRANK, & ABORT1,ABORT2 COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT COMMON/BC/NELA C CHARACTER*5 ELEM CHARACTER*2 ordrng,oprocs,zprocs CHARACTER*5 coord C COMMON/SOLVR1/matslv,nmaxit,nnvvcc,iiuunn,iissoo,nactdi COMMON/SOLVR2/ritmax,closur COMMON/SOLVR3/ordrng,oprocs,zprocs,coord C INTEGER BNDWTH COMMON/D4COM6/BNDWTH,NCRTM,NEQRM,NEQUH,IDIM,IRELSZ C SAVE ICALL,N,iteruc,iprpro,izero0 DATA ICALL,iteruc,iprpro,izero0/0,0,0,0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of LINEQ C C C********************************************************************** c.....for MATSLV=1, call LINEQ of t2m.f (1991-version of TOUGH2), c which here is renamed LINEQ1; it will call MA28. if(matslv.eq.1) then call lineq1 return endif C********************************************************************** C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) c 899 FORMAT(6X,'LINEQ 0.91 CG 31 January 1994',6X, c 899 FORMAT(6X,'LINEQ 2.00 11 January 1997',6X, c 899 FORMAT(6X,'LINEQ 2.00 21 May 1999',6X, c 899 FORMAT(6X,'LINEQ 2.00 26 May 1999',6X, c 899 FORMAT(6X,'LINEQ 2.00 1 October 1999',6X, c 899 FORMAT(6X,'LINEQ 2.00 4 October 1999',6X, 899 FORMAT(6X,'LINEQ 2.00 2 May 2006',6X, & 'Interface for linear equation solvers T2CG2'/ & 47X,'Can call a direct solver or a package of ', & 'conjugate gradient solvers') C MA = MA+1 IF(MA.EQ.1) N=NEQ*NELA C C C********************************************************************* C* * C* ACCOUNTING FOR ZEROs ON THE MAIN DIAGONAL * C* * C********************************************************************* C C izerod = 0 iprpro = 0 IF(MATSLV.NE.6) THEN NNF = N*NEQ DO 20 I=1,NNF IF(IRN(I).EQ.ICN(I).AND.ABS(CO(I)).EQ.0.0e0) izerod = izerod+1 20 CONTINUE END IF C if(izerod.ne.izero0) then WRITE(15,6003) kcyc,iter,izerod,zprocs,oprocs izero0 = izerod endif C zertio=1.0d2*izerod/(neq*nela) IF(izerod.GT.0) THEN IF(zprocs.EQ.'Z1') iprpro = 1 IF(zprocs.EQ.'Z2') iprpro = 2 IF(zprocs.EQ.'Z3') iprpro = 3 IF(zprocs.EQ.'Z4') iprpro = 4 IF(oprocs.NE.'O0') THEN IF(zprocs.EQ.'Z0'.OR.zprocs.EQ.'Z1') THEN IF(zertio.GT.2.0e+1) THEN PRINT 6004, zertio,oprocs oprocs='O0' END IF END IF END IF END IF C 6003 FORMAT(1X,'At KCYC=',i5,' and ITER=',i5,', IZEROD=',I5, & ', ZPROCS = ',a2,' and OPROCS = ',a2) 6004 FORMAT(//,' ',15('WARNING-'),//T14,F5.2,'% of the elements on ', & 'the main diagonal of the Jacobian matrix are zeros'/ & T14,'The matrix preprocessor OPROCS = ',A2,' cannot be used'/ & T14,'Action taken: reset OPROCS to *O0* (no O-preprocessing);', & ' continue execution.') C C INUM = 0 IGOOD = 0 C C-----MA COUNTS CALLS TO LINEQ. C IF(MOP(6).NE.0) PRINT 6005,kcyc,iter,icall,N,NZ,izerod,zertio 6005 FORMAT(' !!!!! L I N E Q !!!!! at [KCYC, ITER] = [',I5,',',I3, &']',' ICALL =',I5,' N =',I4,' NZ =',I8,' IZEROD =',I7, &' or',F6.2,' % zeros') C C C IF(MOP(6).GE.7) THEN PRINT 6015 PRINT 6010,(IRN(NN),ICN(NN),CO(NN),NN=1,NZ) END IF C 6010 FORMAT(5(1X,2I5,E14.6)) 6015 FORMAT(/,' MATRIX OF COEFFICIENTS',/) C C********************************************************************* C* * C* DETERMINATION OF THE BANDWIDTHS OF THE U AND L MATRICES * C* OR PLACEMENT OF THE ELEMENTS INTO THE CG SOLUTION ARRAY * C* * C********************************************************************* C IF(icall.EQ.1) THEN co(mnzp1) = 0.0e0 r(mnetp1) = 0.0e0 END IF C C C IF(MATSLV.EQ.6.AND.icall.EQ.1) THEN iddfup = 0 iddfdn = 0 DO 400 j=1,nz iddffu = icn(j)-irn(j) iddffd = irn(j)-icn(j) IF(iddffu.GE.iddfup) iddfup=iddffu IF(iddffd.GE.iddfdn) iddfdn=iddffd 400 CONTINUE C matord = neq*nela nsupdg = iddfup nsubdg = iddfdn ntotd = 2*nsubdg+nsupdg+1 navdia = lenw/matord C IF(ntotd.GT.navdia) THEN PRINT 6020, ntotd,navdia,matord*ntotd STOP END IF END IF C C C IF(icall.EQ.1.AND.NEQ.GT.1.AND.MATSLV.NE.6.AND.MATSLV.NE.1) THEN IF(oprocs.NE.'O0'.OR.(izerod.NE.0.AND.iprpro.GT.1)) THEN CALL ELINDX CALL REASSN END IF END IF C C C*********************************************************************** C* * C* MATRIX SOLUTION * C* * C*********************************************************************** C IF(MATSLV.EQ.6) THEN CALL LUBAND(matord,nz,nsubdg,nsupdg,ntotd,ab, & matord,r,JVECT,info) C ELSE C IF(izerod.NE.0.AND.NEQ.GT.1.AND.zprocs.NE.'Z0') THEN CALL MTRXIN(iprpro) END IF c IF(oprocs.NE.'O0'.AND.NEQ.EQ.1) GO TO 415 IF(oprocs.eq.'O0') THEN GO TO 415 ELSE IF(oprocs.eq.'O1') THEN CALL MTRXPR(1) ELSE IF(oprocs.eq.'O2') THEN CALL MTRXPR(2) ELSE IF(oprocs.eq.'O3') THEN CALL MTRXPR(3) ELSE IF(oprocs.eq.'O4') THEN CALL MTRXIN(4) END IF C 415 CONTINUE C DO 440 i=1,n wkarea(i) = 0.0e0 440 CONTINUE C IF(MOP(6).NE.0) CALL THYME(0,TS,TT) C IF(MATSLV.EQ.2) THEN CALL DSLUBC(N,r,wkarea,NZ,irn,icn,co, & CLOSUR,NMAXIT,ITERU,ERR,IERR,IUNIT, & AB,LENW,jvect,LENIW) ELSE IF(MATSLV.EQ.3) THEN CALL DSLUCS(N,r,wkarea,NZ,irn,icn,co, & CLOSUR,NMAXIT,ITERU,ERR,IERR,IUNIT, & AB,LENW,jvect,LENIW) ELSE IF(MATSLV.EQ.4) THEN CALL DSLUGM(N,r,wkarea,NZ,irn,icn,co,NVECTR, & CLOSUR,NMAXIT,ITERU,ERR, & IERR,IUNIT,AB,LENW,jvect,LENIW) ELSE IF(MATSLV.EQ.5) THEN CALL DLUSTB(N,r,wkarea,NZ,irn,icn,co, & CLOSUR,NMAXIT,ITERU,ERR,IERR,IUNIT, & AB,LENW,jvect,LENIW) END IF C iteruc = iteruc+iteru C WRITE(15,6045) kcyc,iter,deltex,ierr,err,iteru,iteruc igood=ierr if(igood.ne.0) print 6021 6021 format(' failure in linear equation solution') c DO 450 I=1,N R(I) = wkarea(I) 450 CONTINUE C IF(MOP(6).NE.0) THEN CALL THYME(1,TSS,TT) PRINT 6040,TSS END IF END IF C C 6020 FORMAT(//,20('ERROR-'),//,T33, & ' S I M U L A T I O N A B O R T E D', & /,T40,'DECLARED BANDWIDTH SMALLER THAN NEEDED', & /,T33,' PLEASE CORRECT AND TRY AGAIN', & //,T20,'THE NUMBER OF NEEDED DIAGONALS IS', I4, & ' WHILE THE AVAILABLE NUMBER IS ', I4, & //,T27,'THE PARAMETER nrwork MUST BE INCREASED TO ', & 'AT LEAST ',I10, & //,20('ERROR-')) 6030 FORMAT(' EPS = ',E10.4, ' RMIN = ',E10.4, & ' RESID = ',E10.4,' IRNCP = ',I6, & ' MINIRN = ',I6, ' MINICN = ',I6,'IRANK = ',I4) C 6040 FORMAT(' SOLUTION TIME = ',1PE12.4,' SECONDS') 6045 FORMAT(T5,' At [',I4,',',I3,']',' DELT=',E12.6, & ' IERR=',I1,'& ERR=',1PE12.6,' IT=',I5,' ITC=',I10) C C C********************************************************************* C* * C* UPDATE CHANGES * C* * C********************************************************************* C 455 IF(MOP(6).GE.5) PRINT 6050 6050 FORMAT(/,' ===== INCREMENTS ==== in order of primary', & ' variables',/) C DO 500 NN=1,NELA NLOC = (NN-1)*NEQ NLOCP = (NN-1)*NK1 C IF(MOP(6).GE.5) PRINT 6055,ELEM(NN),(R(NLOC+K),K=1,NEQ) C DO 480 K=1,NEQ DX(NLOCP+K) = DX(NLOCP+K)+WNR*R(NLOC+K) 480 CONTINUE 500 CONTINUE C 6055 FORMAT(' AT ELEMENT *',A5,'* ',8(1X,E12.6)) C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of LINEQ C C RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE LUBAND(N,NZ,KL,KU,LDAB,AB,LDB,B,IPIV8,INFO) C C c IMPLICIT REAL*8 (A-H,O-Z) C C*********************************************************************** C* * C*! C O M M O N D E C L A R A T I O N S * C* * C*********************************************************************** C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) C COMMON/SVZ/NOITE,MOP(24) C CHARACTER*2 ordrng,oprocs,zprocs CHARACTER*5 coord COMMON/SOLVR3/ordrng,oprocs,zprocs,coord C INTEGER IPIV8(*) c REAL*8 AB(LDAB,*),B(LDB,*) Dimension AB(LDAB,*),B(LDB,1) C SAVE ICALL DATA ICALL/0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of LUBAND C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'LUBAND 1.0 12 January 1997',6X, & 'Direct banded matrix solver using LU decomposition') C C C********************************************************************* C* * C* I N I T I A L I Z A T I O N * C* * C********************************************************************* C info = 0 DO 502 j=1,n DO 502 i=1,LDAB AB(i,j) = 0.0e0 502 CONTINUE C C********************************************************************* C* * C* PLACEMENT OF MATRIX ELEMENTS INTO THE AB ARRAY * C* * C********************************************************************* C DO 503 i=1,NZ nfrst = KL+KU+1+irn(i)-icn(i) AB(nfrst,icn(i)) = co(i) 503 CONTINUE C IF(MOP(6).NE.0) CALL THYME(0,TS,TT) C C********************************************************************* C* * C* LU DECOMPOSITION * C* * C********************************************************************* C CALL DGBTRF(N,N,KL,KU,AB,LDAB,IPIV8,INFO) C IF(MOP(6).NE.0) THEN CALL THYME(1,TD,TT) PRINT 6001,TD,IFLAG END IF C 6001 FORMAT(' LU DECOMPOSITION TIME = ',1PE12.4, & ' SECONDS',' ********** INFO =',I4) C C********************************************************************* C* * C* SOLUTION USING THE LU FACTORIZATION COMPUTED BY DGBTRF * C* * C********************************************************************* C IF(info.eq.0) THEN C IF(MOP(6).NE.0) CALL THYME(0,TS,TT) C CALL DGBTRS(N,KL,KU,AB,LDAB,IPIV8,B,LDB) C IF(MOP(6).NE.0) THEN CALL THYME(1,TSS,TT) PRINT 6002, TSS END IF C ELSE PRINT 6003, info STOP END IF C 6002 FORMAT(' SOLUTION TIME = ',1PE12.4,' SECONDS') 6003 FORMAT(//,20('ERROR-'),//,T33, & ' S I M U L A T I O N A B O R T E D', & /,T24,'THE UPPER TRIANGULAR MATRIX ELEMENT U(I,I) ', & 'IS EXACTLY ZERO AT I = ',I5, & /,T38,'THE FACTORIZATION IS COMPLETED BUT DIVISION', & /,T38,'BY ZERO WILL OCCUR IF SOLUTION IS ATTEMPTED', & //,20('ERROR-')) C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of LUBAND C C RETURN C END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE MTRXPR(ilevel) C C C*********************************************************************** C* * C* THE FILE 'T2' IS INCLUDED * C* * C*********************************************************************** C C INCLUDE 'T2' C C*********************************************************************** C* * C*! C O M M O N D E C L A R A T I O N S * C* * C*********************************************************************** C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) COMMON/P4/R(1) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/BC/NELA COMMON/AMMIS/MA,IPIV,U,IAB,NZ C CHARACTER*2 ordrng,oprocs,zprocs CHARACTER*5 coord COMMON/SOLVR3/ordrng,oprocs,zprocs,coord COMMON/SOLVR1/matslv,nmaxit,nnvvcc,iiuunn,iissoo,nactdi C SAVE ICALL DATA ICALL/0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of MTRXPR C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(6X,'MTRXPR 1.0 19 January 1997',6X, & 'Routine for O-preprocessing ', & 'of the Jacobian') C C C MXMYMZ = NELA IPHASE = NEQ 9801 format(t5,'=>=>=> MTRXPR Flag # ',i2) C C C IF(IPHASE.EQ.1) GO TO 50 C C*********************************************************************** C* * C* ELIMINATION OF LEFT-OFF-M.BAND ELEMENTS * C* * C*********************************************************************** C DO 45 IJK=1,MXMYMZ DO 15 LB=1,IPHASE-1 LB1=LB+1 DO 12 LA=LB1,IPHASE nb1 = NB0(LB,LA,IJK) nb2 = NB0(LB,LB,IJK) c QUOT = CO(nb1)/CO(nb2) R(NX0(LA,IJK)) = R(NX0(LA,IJK))-QUOT*R(NX0(LB,IJK)) C DO 4 L=LB1,IPHASE CO(NB0(L,LA,IJK)) = CO(NB0(L,LA,IJK)) & -QUOT*CO(NB0(L,LB,IJK)) 4 CONTINUE C DO 10 L=1,IPHASE IF(IJKMM(ijk).LT.1) GO TO 6 DO 5 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,L,LA)) = CO(NA0(NA0LOC,L,LA)) & -CO(NA0(NA0LOC,L,LB))*QUOT 5 CONTINUE c 6 IF(IJKPP(ijk).LT.1) GO TO 10 DO 8 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,L,LA)) = CO(NC0(NC0LOC,L,LA)) & -CO(NC0(NC0LOC,L,LB))*QUOT 8 CONTINUE 10 CONTINUE 12 CONTINUE C DO 14 L=LB1,IPHASE CO(NB0(LB,L,IJK)) = 0.0e0 14 CONTINUE 15 CONTINUE IF(ilevel.EQ.1) GO TO 45 C C*********************************************************************** C* * C* ELIMINATION OF RIGHT-OFF-M.BAND ELEMENTS * C* * C*********************************************************************** C DO 25 LB=IPHASE,2,-1 LB1=LB-1 DO 23 LA=LB1,1,-1 QUOT = CO(NB0(LB,LA,IJK))/CO(NB0(LB,LB,IJK)) R(NX0(LA,IJK)) = R(NX0(LA,IJK))-QUOT*R(NX0(LB,IJK)) C DO 22 L=1,IPHASE IF(IJKMM(ijk).LT.1) GO TO 20 DO 18 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,L,LA)) = CO(NA0(NA0LOC,L,LA)) & -CO(NA0(NA0LOC,L,LB))*QUOT 18 CONTINUE c 20 IF(IJKPP(ijk).LT.1) GO TO 22 DO 21 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,L,LA)) = CO(NC0(NC0LOC,L,LA)) & -CO(NC0(NC0LOC,L,LB))*QUOT 21 CONTINUE 22 CONTINUE 23 CONTINUE C DO 24 L=LB1,1,-1 CO(NB0(LB,L,IJK)) = 0.0e0 24 CONTINUE 25 CONTINUE IF(ilevel.EQ.2) GO TO 45 C C*********************************************************************** C* * C* N O R M A L I Z A T I O N * C* * C*********************************************************************** C DO 30 LB=1,IPHASE QQQ =-1.0e0/CO(NB0(LB,LB,IJK)) CO(NB0(LB,LB,IJK)) =-1.0e0 R(NX0(LB,IJK)) = R(NX0(LB,IJK))*QQQ C DO 29 L=1,IPHASE IF(IJKMM(ijk).LT.1) GO TO 27 DO 26 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,L,LB)) = CO(NA0(NA0LOC,L,LB))*QQQ 26 CONTINUE c 27 IF(IJKPP(ijk).LT.1) GO TO 29 DO 28 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,L,LB)) = CO(NC0(NC0LOC,L,LB))*QQQ 28 CONTINUE 29 CONTINUE 30 CONTINUE 45 CONTINUE GO TO 999 C C*********************************************************************** C* * C* THE NEQ=1 CASE * C* * C*********************************************************************** C 50 DO 100 IJK=1,MXMYMZ QQQ =-1.0e0/CO(NB0(1,1,IJK)) c CO(NB0(1,1,IJK)) =-1.0e0 R(NX0(1,IJK)) = R(NX0(1,IJK))*QQQ C IF(IJKMM(ijk).LT.1) GO TO 56 DO 55 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,1,1)) = CO(NA0(NA0LOC,1,1))*QQQ 55 CONTINUE c 56 IF(IJKPP(ijk).LT.1) GO TO 100 DO 60 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,1,1)) = CO(NC0(NC0LOC,1,1))*QQQ 60 CONTINUE c 100 CONTINUE C 999 CONTINUE C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of MTRXPR C C RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE REASSN C C C*********************************************************************** C* * C* THE FILE 'T2' IS INCLUDED * C* * C*********************************************************************** C C INCLUDE 'T2' C C*********************************************************************** C* * C*! C O M M O N D E C L A R A T I O N S * C* * C*********************************************************************** C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) COMMON/L7/JVECT(1) COMMON/P4/R(1) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/BC/NELA C COMMON/AMMIS/MA,IPIV,U,IAB,NZ COMMON/lub3/NSUPDI,NSUBDI,mnzp1,mnetp1,mnelp1,nnnbig common/soll/lenw,leniw C SAVE ICALL DATA ICALL/0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of REASSN C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) ! 899 FORMAT(6X,'REASSN 1.0 17 January 1997',6X, 899 FORMAT(6X,'REASSN 1.1 09 September 2000',6X, & 'Establish the sparsity pattern ', & 'of the Jacobian for the Z preprocessors') C C C C*********************************************************************** C* * C* NUMBERING TO MAP CO TO THE A0,C0,B0,X0 SYSTEM * C* * C*********************************************************************** C C DO 100 iee=1,neq*nela IJKR = iee/neq IMODR = MOD(iee,neq) IF(IMODR.EQ.0) THEN IJKMEL = IJKR IPH1 = neq ELSE IJKMEL = IJKR+1 IPH1 = IMODR END IF NX0(IPH1,IJKMEL) = iee 100 CONTINUE C C C C*********************************************************************** C* * C* INITIALIZATION * C* * C*********************************************************************** C C ISUMPP(0) = 0 ISUMMM(0) = 0 DO 136 ijk=1,MNEL+1 IJKPP(ijk) = 0 ISUMPP(ijk) = 0 c IJKMM(ijk) = 0 ISUMMM(ijk) = 0 DO 136 i2=1,neq DO 136 i1=1,neq NB0(i1,i2,ijk) = mnzp1 136 CONTINUE C C C DO 138 i2=1,neq DO 138 i1=1,neq DO 138 ic8=1,NCONUP NA0(ic8,i1,i2) = mnzp1 138 CONTINUE C DO 140 i2=1,neq DO 140 i1=1,neq DO 140 ic8=1,NCONDN NC0(ic8,i1,i2) = mnzp1 140 CONTINUE C C C C*********************************************************************** C* * C* DETERMINING NEIGBORING ELEMENT SPECIFICS ALONG SAME ROW * C* * C*********************************************************************** C C DO 200 n=1,nela c IF(no(iijjkk(n)+1).GT.n) THEN IJKPP(n) = iijjkk(n+1)-iijjkk(n)-1 jvect(n) = iijjkk(n)+1 jvect(nela+n) = iijjkk(n+1)-1 c IJKMM(n) = 0 jvect(2*nela+n) = 0 jvect(3*nela+n) = 0 GO TO 200 END IF c IF(no(iijjkk(n+1)-1).LT.n) THEN IJKPP(n) = 0 jvect(n) = 0 jvect(nela+n) = 0 c IJKMM(n) = iijjkk(n+1)-iijjkk(n)-1 jvect(2*nela+n) = iijjkk(n)+1 jvect(3*nela+n) = iijjkk(n+1)-1 GO TO 200 END IF c DO 150 index = iijjkk(n)+1, iijjkk(n+1)-1 IF(n.GT.no(index).AND.n.LT.no(index+1)) THEN IJKPP(n) = iijjkk(n+1)-1-index jvect(n) = index+1 jvect(nela+n) = iijjkk(n+1)-1 c IJKMM(n) = index-iijjkk(n) jvect(2*nela+n) = iijjkk(n)+1 jvect(3*nela+n) = index GO TO 200 END IF 150 CONTINUE 200 CONTINUE C C C DO 220 n=1,nela ISUMPP(n) = ISUMPP(n-1)+IJKPP(n) ISUMMM(n) = ISUMMM(n-1)+IJKMM(n) 220 CONTINUE C C C IF(ISUMPP(nela).GT.NCONUP) THEN PRINT 6001, ISUMPP(nela) STOP END IF C IF(ISUMMM(nela).GT.NCONDN) THEN PRINT 6002, ISUMMM(nela) STOP END IF C C C C*********************************************************************** C* * C* MAPPING THE GLOBAL MATRIX ELEMENT NUMBER TO CELL-SPECIFIC POINTERS * C* * C*********************************************************************** C C DO 400 izz=1,nz IJKR = irn(izz)/neq IMODR = MOD(irn(izz),neq) IF(IMODR.EQ.0) THEN IJKMEL = IJKR IPH2 = neq ELSE IJKMEL = IJKR+1 IPH2 = IMODR END IF i8 = ijkmel C C C IJKC = icn(izz)/neq IMODC = MOD(icn(izz),neq) IF(IMODC.EQ.0) THEN IJKCON = IJKC IPH1 = neq ELSE IJKCON = IJKC+1 IPH1 = IMODC END IF C C C IF(IJKMEL.EQ.IJKCON) THEN NB0(IPH1,IPH2,i8) = izz ELSE c IF(i8.LT.IJKCON) THEN C IF(jvect(i8).EQ.0) GO TO 400 nabov = 0 DO 340 index = jvect(i8),jvect(nela+i8) id = no(index) nabov = nabov+1 IF(id.EQ.ijkcon) THEN jjj = nabov GO TO 342 END IF 340 CONTINUE 342 NC0(jjj+ISUMPP(i8-1),IPH1,IPH2) = izz ELSE C IF(jvect(2*nela+i8).EQ.0) GO TO 400 nbelo = 0 DO 380 index = jvect(2*nela+i8),jvect(3*nela+i8) id = no(index) nbelo = nbelo+1 IF(id.EQ.ijkcon) THEN jjj = nbelo GO TO 382 END IF 380 CONTINUE 382 NA0(jjj+ISUMMM(i8-1),IPH1,IPH2) = izz END IF END IF C C C 400 CONTINUE C C C DO 500 ii = 1,leniw jvect(ii) = 0 500 CONTINUE C C C 6001 FORMAT(//,20('ERROR-'),//,T33, & ' S I M U L A T I O N A B O R T E D', & /,T24,'THE FIRST DIMENSION OF THE isumpp ARRAY ', & 'IS INSUFFICIENT FOR THIS PROBLEM', & /,T12,'THE PARAMETER nconup MUST BE AT LEAST ',I10, & //,T33,' CORRECT AND TRY AGAIN', & //,20('ERROR-')) 6002 FORMAT(//,20('ERROR-'),//,T33, & ' S I M U L A T I O N A B O R T E D', & /,T24,'THE FIRST DIMENSION OF THE isummm ARRAY ', & 'IS INSUFFICIENT FOR THIS PROBLEM', & /,T12,'THE PARAMETER ncondn MUST BE AT LEAST ',I10, & //,T33,' CORRECT AND TRY AGAIN', & //,20('ERROR-')) C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of REASSN C C 999 RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE MTRXIN(IOPT) C C C*********************************************************************** C* * C* THE FILE 'T2' IS INCLUDED * C* * C*********************************************************************** C C INCLUDE 'T2' C PARAMETER (SEED = 1.0e-25) C** C*********************************************************************** C* * C*! C O M M O N D E C L A R A T I O N S * C* * C*********************************************************************** C COMMON/L1/IRN(1) COMMON/L2/ICN(1) COMMON/L3/CO(1) COMMON/P4/R(1) C DIMENSION B0INV(5,5),TA0(5,5),TC0(5,5),TR0(5) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/BC/NELA COMMON/AMMIS/MA,IPIV,U,IAB,NZ C SAVE ICALL DATA ICALL/0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of MTRXIN C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) ! 899 FORMAT(6X,'MTRXIN 1.0 24 January 1997',6X, 899 FORMAT(6X,'MTRXIN 1.1 10 September 2000',6X, & 'Routine for Z-preprocessing ', & 'of the Jacobian') C C C MXMYMZ = NELA IPHASE = NEQ C C C*********************************************************************** C* * C* REPLACEMENT OF 0 BY SEED ON THE MAIN DIAGONAL * C* * C*********************************************************************** C C C IF(IOPT.GT.1) GO TO 61 c DO 20 I=1,NELA*NEQ DO 20 I=1,NZ IF(IRN(I).EQ.ICN(I).AND.ABS(CO(I)).EQ.0.0e0) THEN CO(I) = SEED END IF 20 CONTINUE GO TO 999 C C C*********************************************************************** C* * C* ADDITION OF EQUATIONS TO PREVENT ZEROS ON THE MAIN DIAGONAL * C* * C*********************************************************************** C C 61 IF(IOPT.GT.2) GO TO 101 DO 100 IJK=1,MXMYMZ DO 90 LA=1,IPHASE C IF(CO(NB0(LA,LA,IJK)).EQ.0.0e0) THEN DO 80 LB=1,IPHASE IF(LB.EQ.LA) GO TO 80 C IF(CO(NB0(LA,LB,IJK)).NE.0.0e0) THEN C R(NX0(LA,IJK)) = R(NX0(LA,IJK))+1.75e0*R(NX0(LB,IJK)) DO 70 L1=1,IPHASE CO(NB0(L1,LA,IJK)) & = CO(NB0(L1,LA,IJK)) & +1.75e0*CO(NB0(L1,LB,IJK)) 70 CONTINUE DO 75 L1=1,IPHASE IF(IJKMM(ijk).LT.1) GO TO 73 DO 72 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,L1,LA)) & = CO(NA0(NA0LOC,L1,LA)) & +1.75e0*CO(NA0(NA0LOC,L1,LB)) 72 CONTINUE 73 IF(IJKPP(ijk).LT.1) GO TO 75 DO 74 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,L1,LA)) & = CO(NC0(NC0LOC,L1,LA)) & +1.75e0*CO(NC0(NC0LOC,L1,LB)) 74 CONTINUE 75 CONTINUE GO TO 90 END IF 80 CONTINUE C END IF C 90 CONTINUE 100 CONTINUE GO TO 999 C C C*********************************************************************** C* * C* NORMALIZATION + ADDITION OF EQUATIONS * C* * C*********************************************************************** C C 101 IF(IOPT.GT.4) GO TO 301 DO 200 IJK=1,MXMYMZ C ROWMAX = 0.0e0 DO 150 LB=1,IPHASE DO 115 LA=1,IPHASE ABB0 = ABS(CO(NB0(LA,LB,IJK))) IF(ABB0.GT.ROWMAX) THEN ROWMAX = ABB0 ROMX = SIGN(ABB0,CO(NB0(LA,LB,IJK))) END IF IF(iopt.EQ.3) GO TO 115 c IF(IJKMM(ijk).LT.1) GO TO 113 DO 112 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) ABA0 = ABS(CO(NA0(NA0LOC,LA,LB))) IF(ABA0.GT.ROWMAX) THEN ROWMAX = ABA0 ROMX = SIGN(ABA0,CO(NA0(NA0LOC,LA,LB))) END IF 112 CONTINUE c 113 IF(IJKPP(ijk).LT.1) GO TO 115 DO 114 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) ABC0 = ABS(CO(NC0(NC0LOC,LA,LB))) IF(ABC0.GT.ROWMAX) THEN ROWMAX = ABC0 ROMX = SIGN(ABC0,CO(NC0(NC0LOC,LA,LB))) END IF 114 CONTINUE c 115 CONTINUE C QUOT = 1.0e0/ROMX C R(NX0(LB,IJK)) = R(NX0(LB,IJK))*QUOT DO 120 LA=1,IPHASE CO(NB0(LA,LB,IJK)) = CO(NB0(LA,LB,IJK))*QUOT IF(IJKMM(ijk).LT.1) GO TO 117 DO 116 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,LA,LB)) = CO(NA0(NA0LOC,LA,LB))*QUOT 116 CONTINUE c 117 IF(IJKPP(ijk).LT.1) GO TO 120 DO 118 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,LA,LB)) = CO(NC0(NC0LOC,LA,LB))*QUOT 118 CONTINUE c 120 CONTINUE 150 CONTINUE 200 CONTINUE C C C 255 DO 300 IJK=1,MXMYMZ DO 290 LA=1,IPHASE C IF(CO(NB0(LA,LA,IJK)).EQ.0.0e0) THEN DO 280 LB=1,IPHASE IF(LB.EQ.LA) GO TO 280 C IF(CO(NB0(LA,LB,IJK)).NE.0.0e0) THEN C R(NX0(LA,IJK)) = R(NX0(LA,IJK))+1.75e0*R(NX0(LB,IJK)) DO 270 L1=1,IPHASE CO(NB0(L1,LA,IJK)) = CO(NB0(L1,LA,IJK)) & +1.75e0*CO(NB0(L1,LB,IJK)) 270 CONTINUE DO 275 L1=1,IPHASE c IF(IJKMM(ijk).LT.1) GO TO 273 DO 272 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) CO(NA0(NA0LOC,L1,LA)) & = CO(NA0(NA0LOC,L1,LA)) & +1.75e0*CO(NA0(NA0LOC,L1,LB)) 272 CONTINUE c 273 IF(IJKPP(ijk).LT.1) GO TO 275 DO 274 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) CO(NC0(NC0LOC,L1,LA)) & = CO(NC0(NC0LOC,L1,LA)) & +1.75e0*CO(NC0(NC0LOC,L1,LB)) 274 CONTINUE c 275 CONTINUE GO TO 290 END IF 280 CONTINUE C END IF C 290 CONTINUE 300 CONTINUE GO TO 999 C C C*********************************************************************** C* * C* MULTIPLICATION BY THE INVERSE OF THE MAIN-DIAGONAL SUBMATRIX * C* * C*********************************************************************** C C 301 IF(NEQ.EQ.3) GO TO 401 C C C DO 400 IJK=1,MXMYMZ C DETR = CO(NB0(1,1,IJK))*CO(NB0(2,2,IJK))- & CO(NB0(1,2,IJK))*CO(NB0(2,1,IJK)) IF(ABS(DETR).LE.SEED) THEN PRINT 6001 STOP END IF DETI = 1.0e0/DETR C B0INV(1,1) = CO(NB0(2,2,IJK))*DETI B0INV(1,2) =-CO(NB0(1,2,IJK))*DETI B0INV(2,1) =-CO(NB0(2,1,IJK))*DETI B0INV(2,2) = CO(NB0(1,1,IJK))*DETI C CO(NB0(1,1,IJK)) = 1.0e0 CO(NB0(2,2,IJK)) = 1.0e0 CO(NB0(1,2,IJK)) = 0.0e0 CO(NB0(2,1,IJK)) = 0.0e0 C DO 390 LB=1,IPHASE TR0(LB) = 0.0e0 DO 355 L1=1,IPHASE TR0(LB) = TR0(LB)+B0INV(L1,LB)*R(NX0(L1,IJK)) 355 CONTINUE R(NX0(LB,IJK)) = TR0(LB) C DO 380 LA=1,IPHASE TA0(LA,LB) = 0.0e0 TC0(LA,LB) = 0.0e0 c IF(IJKMM(ijk).LT.1) GO TO 371 DO 370 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) DO 360 L1=1,IPHASE TA0(LA,LB) = TA0(LA,LB) & +CO(NA0(NA0LOC,L1,LB))*B0INV(LA,L1) 360 CONTINUE CO(NA0(NA0LOC,LA,LB)) = TA0(LA,LB) 370 CONTINUE c 371 IF(IJKPP(ijk).LT.1) GO TO 380 DO 375 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) DO 372 L1=1,IPHASE TC0(LA,LB) = TC0(LA,LB) & +CO(NC0(NC0LOC,L1,LB))*B0INV(LA,L1) 372 CONTINUE CO(NC0(NC0LOC,LA,LB)) = TC0(LA,LB) 375 CONTINUE c 380 CONTINUE 390 CONTINUE 400 CONTINUE GO TO 999 C C C 401 DO 500 IJK=1,MXMYMZ C DET1 = CO(NB0(2,2,IJK))*CO(NB0(3,3,IJK)) & -CO(NB0(2,3,IJK))*CO(NB0(3,2,IJK)) DET2 = CO(NB0(2,1,IJK))*CO(NB0(3,3,IJK)) & -CO(NB0(2,3,IJK))*CO(NB0(3,1,IJK)) DET3 = CO(NB0(2,1,IJK))*CO(NB0(3,2,IJK)) & -CO(NB0(2,2,IJK))*CO(NB0(3,1,IJK)) DETR = CO(NB0(1,1,IJK))*DET1 & +CO(NB0(1,2,IJK))*DET2 & +CO(NB0(1,3,IJK))*DET3 C IF(ABS(DETR).LE.SEED) THEN PRINT 6001 STOP END IF DETI = 1.0e0/DETR C B0INV(1,1) = DET1*DETI B0INV(1,2) =(CO(NB0(1,3,IJK))*CO(NB0(3,2,IJK)) & -CO(NB0(1,2,IJK))*CO(NB0(3,3,IJK)))*DETI B0INV(1,3) =(CO(NB0(1,2,IJK))*CO(NB0(2,3,IJK)) & -CO(NB0(1,3,IJK))*CO(NB0(2,2,IJK)))*DETI C B0INV(2,1) =-DET2*DETI B0INV(2,2) =(CO(NB0(1,1,IJK))*CO(NB0(3,3,IJK)) & -CO(NB0(1,3,IJK))*CO(NB0(3,1,IJK)))*DETI B0INV(2,3) =(CO(NB0(1,3,IJK))*CO(NB0(2,1,IJK)) & -CO(NB0(1,1,IJK))*CO(NB0(2,3,IJK)))*DETI C B0INV(3,1) = DET3*DETI B0INV(3,2) =(CO(NB0(1,2,IJK))*CO(NB0(3,1,IJK)) & -CO(NB0(1,1,IJK))*CO(NB0(3,2,IJK)))*DETI B0INV(3,3) =(CO(NB0(1,1,IJK))*CO(NB0(2,2,IJK)) & -CO(NB0(1,2,IJK))*CO(NB0(2,1,IJK)))*DETI C CO(NB0(1,1,IJK)) =-1.0e0 CO(NB0(2,2,IJK)) =-1.0e0 CO(NB0(3,3,IJK)) =-1.0e0 CO(NB0(1,2,IJK)) = 0.0e0 CO(NB0(1,3,IJK)) = 0.0e0 CO(NB0(2,1,IJK)) = 0.0e0 CO(NB0(2,3,IJK)) = 0.0e0 CO(NB0(3,1,IJK)) = 0.0e0 CO(NB0(3,2,IJK)) = 0.0e0 C DO 440 LB=1,IPHASE TR0(LB) = 0.0e0 DO 405 L1=1,IPHASE TR0(LB) = TR0(LB)+B0INV(L1,LB)*R(NX0(L1,IJK)) 405 CONTINUE R(NX0(LB,IJK)) =-TR0(LB) C DO 430 LA=1,IPHASE TA0(LA,LB) = 0.0e0 TC0(LA,LB) = 0.0e0 c IF(IJKMM(ijk).LT.1) GO TO 421 DO 420 NDM=1,IJKMM(ijk) NA0LOC = NDM+ISUMMM(ijk-1) DO 410 L1=1,IPHASE TA0(LA,LB) = TA0(LA,LB) & +CO(NA0(NA0LOC,L1,LB))*B0INV(LA,L1) 410 CONTINUE CO(NA0(NA0LOC,LA,LB)) =-TA0(LA,LB) 420 CONTINUE c 421 IF(IJKPP(ijk).LT.1) GO TO 430 DO 425 NDM=1,IJKPP(ijk) NC0LOC = NDM+ISUMPP(ijk-1) DO 422 L1=1,IPHASE TC0(LA,LB) = TC0(LA,LB) & +CO(NC0(NC0LOC,L1,LB))*B0INV(LA,L1) 422 CONTINUE CO(NC0(NC0LOC,LA,LB)) =-TC0(LA,LB) 425 CONTINUE c 430 CONTINUE 440 CONTINUE 500 CONTINUE C 6001 FORMAT(//,20('ERROR-'),//,T40, & ' S I M U L A T I O N H A L T E D', & /,T35,'THE DETERMINANT OF THE MAIN DIAGONAL SUBMATRIX ', & 'IS ZERO',/, & T31,'Preprocessors ZPROCS = Z4 and OPROCS = O4 cannot ', & 'be used! ',/, & T40,' PLEASE CORRECT AND TRY AGAIN',/, & //,20('ERROR-')) C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of MTRXIN C C 999 RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE ELINDX C C********************************************************************* C* * C*! DETERMINE THE POSITIONING ARRAYS ia AND ja * C*! WHICH IDENTIFY THE LOCATIONS OF THE NEIGHBORS * C*! Version 1.00, January 16, 1998 * C* * C********************************************************************* C INCLUDE 'T2' C C********************************************************************* C* * C*! C O M M O N D E C L A R A T I O N S * C* * C********************************************************************* C COMMON/C1/NEX1(1) COMMON/C2/NEX2(1) C COMMON/L7/JVECT(1) C COMMON/NN/NEL,NCON,NOGN,NK,NEQ,NPH,NB,NK1,NEQ1,NBK,NSEC,NFLUX COMMON/BC/NELA COMMON/soll/lenw,leniw C SAVE ICALL DATA ICALL/0/ C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of ELINDX C C ICALL=ICALL+1 IF(ICALL.EQ.1) WRITE(11,899) 899 FORMAT(/6X,'ELINDX 1.00 16 January 1998',6X, & 'Routine for neighbor element indexing') C C********************************************************************* C* * C*! I N I T I A L I Z A T I O N S * C* * C********************************************************************* C DO 100 ii = 1,nel iijjkk(ii) = 1 100 CONTINUE C DO 110 ii = 1,mncon+mnel+1 no(ii) = 0 110 CONTINUE C DO 120 ii = 1,leniw jvect(ii) = 0 120 CONTINUE C C********************************************************************* C* * C*! DETERMINE NUMBER OF NEIGHBORING ELEMENTS AND STORE IN iw * C* * C********************************************************************* C DO 200 n = 1,ncon n1 = nex1(n) n2 = nex2(n) c IF(n1.EQ.0.OR.n2.EQ.0) GO TO 200 c IF(n1.LE.nela) iijjkk(n1) = iijjkk(n1)+1 IF(n2.LE.nela) iijjkk(n2) = iijjkk(n2)+1 c 200 CONTINUE C C********************************************************************* C* * C*! STORE TEMPORARILY THE # OF NEIGHBORS PER ELEMENT ii IN jvect * C*! (FIRST nela ELEMENTS OF jvect) * C* * C********************************************************************* C DO 250 ii = 1,nel jvect(ii) = iijjkk(ii) 250 CONTINUE C C********************************************************************* C* * C*! DETERMINE THE CUMULATIVE NUMBER OF NEIGBORING ELEMENTS * C*! STORE IN jvect - FROM jvect(mshift) TO jvect(mshift+nel) * C* * C********************************************************************* C mshift = nel+1 no(mshift) = 0 C isum = 0 DO 300 ii = 1,nel isum = isum+iijjkk(ii) jvect(mshift+ii) = isum 300 CONTINUE C C********************************************************************* C* * C*! DETERMINE THE iw AND ja ARRAYS WITH THE NEIGHBOR INFORMATION * C* * C********************************************************************* C DO 320 ii = 1,nel iijjkk(ii) = 1 no(jvect(mshift+ii-1)+1) = ii 320 CONTINUE C DO 350 n = 1,ncon n1 = nex1(n) n2 = nex2(n) c IF(n1.GT.nela.AND.n2.GT.nela) GO TO 350 IF(n1.EQ.0.OR.n2.EQ.0) GO TO 350 c iijjkk(n1) = iijjkk(n1)+1 inx1 = jvect(mshift+n1-1)+iijjkk(n1) no(inx1) = n2 c iijjkk(n2) = iijjkk(n2)+1 inx2 = jvect(mshift+n2-1)+iijjkk(n2) no(inx2) = n1 c 350 CONTINUE C C! C iijjkk(1) = 1 no(1) = 1 DO 400 ii = 1,nel c itemp = 0 iijjkk(ii+1) = iijjkk(1)+jvect(mshift+ii) no(iijjkk(ii+1)) = ii+1 c c ------------ c ...... Sort the subarrays ja(ii), ii=iijjkk(ii),...,iijjkk(ii+1)-1 c ------------ c istart = iijjkk(ii) numsrt = iijjkk(ii+1)-iijjkk(ii) c CALL SORT(no(istart),numsrt) c c ------------ c ...... Put the diagonal element first, c ...... swap with element in sorted subarray c ------------ c DO 380 i5 = iijjkk(ii),iijjkk(ii+1)-1 IF(no(i5).EQ.ii) THEN itemp = i5 GO TO 390 END IF 380 CONTINUE c 390 no(itemp) = no(iijjkk(ii)) no(iijjkk(ii)) = ii c 400 CONTINUE C DO 420 ii = 1,leniw jvect(ii) = 0 420 CONTINUE C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of ELINDX C C 9999 RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C*********************************************************************** C C C SUBROUTINE SORT(iar,n) C C********************************************************************* C* * C*! SORTING THE VECTOR iar IN ASCENDING ORDER * C*! Version 1.00, January 14, 1998 * C* * C********************************************************************* C c IMPLICIT REAL*8 (A-H,O-Z) C C********************************************************************* C* * C*! L O C A L A R R A Y S * C* * C********************************************************************* C DIMENSION iar(n) C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> Main body of SORT C C m = n c DO 1000 i5 = 1,n m = m/2 IF(m.EQ.0) GO TO 9999 c max = n-m DO 100 j = 1,max c DO 50 k = j,1,-m iert = iar(k+m) IF(iert.GE.iar(k)) go to 100 c itemp = iar(k+m) iar(k+m) = iar(k) iar(k) = itemp 50 CONTINUE c 100 CONTINUE c 1000 CONTINUE C C C! =>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=> End of SORT C C 9999 RETURN END C C C C*********************************************************************** C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ C***********************************************************************