ifort 에서 32bit processor 에서 64bit 연산을 하게 해주려면 컴파일 옵션을 어떻게 설정해야 하죠?

[Paradism]

1. 컴퓨터 OS  : linux (및 Window)

2. 사용 포트란 버전 : intel fortran (및 Compaq Visual Fortran)

3. 코딩 소스(첨부파일 불가)  :

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)
      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
      NSUPDI = NSUPD
      NSUBDI = NSUBD
      lenw=nrwork
      leniw=niwork
      nnnbig = NBAN*nredm
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     EVAL‎UATE 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 sin
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 EVAL‎UATION 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
      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,'EVAL‎UATE 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
      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
      subroutine sin
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,'SIN      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
      SUBROUTINE LINEQ
C
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,
     AABORT1,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.....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=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>=>   Main body of LINEQ
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,
  899 FORMAT(6X,'LINEQ    2.00      4 October   1999',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
      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
         izero0 = izerod
      endif
C
         zertio=1.0e2*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,
     &              ' and ZPROCS = ',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
      co(mnzp1) = 0.0e0
      r(mnetp1) = 0.0e0
C
C
C
      IF(icall.EQ.1.AND.NEQ.GT.1) THEN
         IF(oprocs.NE.'O0') THEN
            CALL ELINDX
            CALL REASSN
            GO TO 395
         END IF
         IF(izerod.NE.0.and.zprocs.NE.'Z0') THEN
            CALL ELINDX
            CALL REASSN
            GO TO 395
         END IF
      END IF
  395 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 = nnnbig/matord
C
         IF(ntotd.GT.navdia) THEN
            PRINT 6020, ntotd,navdia
            STOP
         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)
            GO TO 415
         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
 
         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
C
         iteruc = iteruc+iteru
C                           
         WRITE(15,6045) kcyc,iter,deltex,ierr,err,iteru,iteruc
      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,
     &       //,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
      RETURN
      END
C
C
      SUBROUTINE LINEQ1
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(/51H !!!! L I N E Q 1 !!!!     !!!!!     !!!!!     MA =,I4,
     A'   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
      SUBROUTINE LUBAND(N,NZ,KL,KU,LDAB,AB,LDB,B,IPIV8,INFO)
C
C
C
c      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/L1/IRN(1)
      COMMON/L2/ICN(1)
      COMMON/L3/CO(1)
C
      COMMON/D4ARAY0/AABB(1)
      COMMON/D4ARAY1/NROW(1)
      COMMON/D4ARAY2/NCOL(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,*)
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
      IF(ordrng.EQ.'ST') THEN
         DO 503 i=1,NZ
            nfrst = KL+KU+1+irn(i)-icn(i)
            AB(nfrst,icn(i)) = co(i)
  503    CONTINUE
      ELSE
         DO 504 i=1,NZ
            nfrst = KL+KU+1+NROW(i)-NCOL(i)
            AB(nfrst,NCOL(i)) = AABB(i)
  504    CONTINUE
      END IF
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
      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
      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(1,ijk).LT.1) GO TO 6
                  DO 5 NDM=1,IJKMM(1,ijk)
                    CO(NA0(NDM,L,LA,IJK)) = CO(NA0(NDM,L,LA,IJK))
     &                                     -CO(NA0(NDM,L,LB,IJK))*QUOT
    5             CONTINUE
c
    6             IF(IJKPP(1,ijk).LT.1) GO TO 10
                  DO 8 NDM=1,IJKPP(1,ijk)
                    CO(NC0(NDM,L,LA,IJK)) = CO(NC0(NDM,L,LA,IJK))
     &                                     -CO(NC0(NDM,L,LB,IJK))*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(1,ijk).LT.1) GO TO 20
                 DO 18 NDM=1,IJKMM(1,ijk)
                    CO(NA0(NDM,L,LA,IJK)) = CO(NA0(NDM,L,LA,IJK))
     &                                     -CO(NA0(NDM,L,LB,IJK))*QUOT
   18            CONTINUE
c
   20            IF(IJKPP(1,ijk).LT.1) GO TO 22
                 DO 21 NDM=1,IJKPP(1,ijk)
                    CO(NC0(NDM,L,LA,IJK)) = CO(NC0(NDM,L,LA,IJK))
     &                                     -CO(NC0(NDM,L,LB,IJK))*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(1,ijk).LT.1) GO TO 27
               DO 26 NDM=1,IJKMM(1,ijk)
                 CO(NA0(NDM,L,LB,IJK)) = CO(NA0(NDM,L,LB,IJK))*QQQ
   26          CONTINUE
c
   27          IF(IJKPP(1,ijk).LT.1) GO TO 29
               DO 28 NDM=1,IJKPP(1,ijk)
                 CO(NC0(NDM,L,LB,IJK)) = CO(NC0(NDM,L,LB,IJK))*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(1,ijk).LT.1) GO TO 56
         DO 55 NDM=1,IJKMM(1,ijk)
            CO(NA0(NDM,1,1,IJK)) = CO(NA0(NDM,1,1,IJK))*QQQ
   55    CONTINUE
c
   56    IF(IJKPP(1,ijk).LT.1) GO TO 100
         DO 60 NDM=1,IJKPP(1,ijk)
            CO(NC0(NDM,1,1,IJK)) = CO(NC0(NDM,1,1,IJK))*QQQ
   60    CONTINUE
c
  100 CONTINUE
C
  999 CONTINUE
C
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
      COMMON/L3/CO(1)
      COMMON/P4/R(1)
      COMMON/L1/IRN(1)
      COMMON/L2/ICN(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
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,
     &          'Establish the sparsity pattern ',
     &          'of the Jacobian for the LUBAND direct solver')
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
      DO 136 ijk=1,nela
         DO 136 i2=1,neq
            DO 136 i1=1,neq
               NB0(i1,i2,ijk) = mnzp1
               DO 136 jjj=1,ndim
                  NC0(jjj,i1,i2,IJK) = mnzp1
                  NA0(jjj,i1,i2,IJK) = mnzp1
  136 CONTINUE
C
C
C
      DO 200 n=1,nela     
c
         IF(no(iijjkk(n)+1).GT.n) THEN
            IJKPP(1,n) = iijjkk(n+1)-iijjkk(n)-1
            IJKPP(2,n) = iijjkk(n)+1
            IJKPP(3,n) = iijjkk(n+1)-1
            IJKMM(1,n) = 0
            IJKMM(2,n) = 0
            IJKMM(3,n) = 0
            GO TO 200
         END IF
c
         IF(no(iijjkk(n+1)-1).LT.n) THEN
            IJKPP(1,n) = 0
            IJKPP(2,n) = 0
            IJKPP(3,n) = 0
            IJKMM(1,n) = iijjkk(n+1)-iijjkk(n)-1
            IJKMM(2,n) = iijjkk(n)+1
            IJKMM(3,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(1,n) = iijjkk(n+1)-1-index
               IJKPP(2,n) = index+1
               IJKPP(3,n) = iijjkk(n+1)-1
               IJKMM(1,n) = index-iijjkk(n)
               IJKMM(2,n) = iijjkk(n)+1
               IJKMM(3,n) = index
               GO TO 200
            END IF
  150    CONTINUE
  200 CONTINUE
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,IJKMEL) = izz
         ELSE
c
            IF(i8.LT.IJKCON) THEN
C
               IF(IJKPP(2,i8).EQ.0) GO TO 400
               nabov = 0
               DO 340 index = IJKPP(2,i8),IJKPP(3,i8)
                  id    = no(index)
                  nabov = nabov+1
                  IF(id.EQ.ijkcon) THEN
                     jjj = nabov
                     GO TO 342
                  END IF
  340          CONTINUE
  342          NC0(jjj,IPH1,IPH2,IJKMEL) = izz
            ELSE
C
               IF(IJKMM(2,i8).EQ.0) GO TO 400
               nbelo = 0
               DO 380 index = IJKMM(2,i8),IJKMM(3,i8)
                  id    = no(index)
                  nbelo = nbelo+1
                  IF(id.EQ.ijkcon) THEN
                     jjj = nbelo
                     GO TO 382
                  END IF
  380          CONTINUE
  382          NA0(jjj,IPH1,IPH2,IJKMEL) = izz
            END IF
         END IF
C
C
C
  400 CONTINUE
C
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
      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,
     &          '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=2,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(1,ijk).LT.1) GO TO 73
                    DO 72 NDM=1,IJKMM(1,ijk)
                      CO(NA0(NDM,L1,LA,IJK)) = CO(NA0(NDM,L1,LA,IJK))
     &                       +1.75e0*CO(NA0(NDM,L1,LB,IJK))
   72               CONTINUE
   73               IF(IJKPP(1,ijk).LT.1) GO TO 75
                    DO 74 NDM=1,IJKPP(1,ijk)
                      CO(NC0(NDM,L1,LA,IJK)) = CO(NC0(NDM,L1,LA,IJK))
     &                       +1.75e0*CO(NC0(NDM,L1,LB,IJK))
   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.3) 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
c
            IF(IJKMM(1,ijk).LT.1) GO TO 113
            DO 112 NDM=1,IJKMM(1,ijk)
              ABA0 = ABS(CO(NA0(NDM,LA,LB,IJK)))
              IF(ABA0.GT.ROWMAX) THEN
                ROWMAX = ABA0
                ROMX   = SIGN(ABA0,CO(NA0(NDM,LA,LB,IJK)))
              END IF
  112       CONTINUE
c
  113       IF(IJKPP(1,ijk).LT.1) GO TO 115
            DO 114 NDM=1,IJKPP(1,ijk)
              ABC0 = ABS(CO(NC0(NDM,LA,LB,IJK)))
              IF(ABC0.GT.ROWMAX) THEN
                ROWMAX = ABC0
                ROMX   = SIGN(ABC0,CO(NC0(NDM,LA,LB,IJK)))
              END IF
  114       CONTINUE
c
  115     CONTINUE
C
          QUOT           = 1.0e0/ROMX
          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(1,ijk).LT.1) GO TO 117
            DO 116 NDM=1,IJKMM(1,ijk)
              CO(NA0(NDM,LA,LB,IJK)) = CO(NA0(NDM,LA,LB,IJK))*QUOT
  116       CONTINUE
c
  117       IF(IJKPP(1,ijk).LT.1) GO TO 120
            DO 118 NDM=1,IJKPP(1,ijk)
              CO(NC0(NDM,LA,LB,IJK)) = CO(NC0(NDM,LA,LB,IJK))*QUOT
  118       CONTINUE
c
  120     CONTINUE
  150   CONTINUE
  200 CONTINUE
C
C
C
      DO 300 IJK=1,MXMYMZ
        DO 290 LA=2,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(1,ijk).LT.1) GO TO 273
                    DO 272 NDM=1,IJKMM(1,ijk)
                      CO(NA0(NDM,L1,LA,IJK)) = CO(NA0(NDM,L1,LA,IJK))
     &                       +1.75e0*CO(NA0(NDM,L1,LB,IJK))
  272               CONTINUE
c
  273               IF(IJKPP(1,ijk).LT.1) GO TO 275
                    DO 274 NDM=1,IJKPP(1,ijk)
                      CO(NC0(NDM,L1,LA,IJK)) = CO(NC0(NDM,L1,LA,IJK))
     &                       +1.75e0*CO(NC0(NDM,L1,LB,IJK))
  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(1,ijk).LT.1) GO TO 371
            DO 370 NDM=1,IJKMM(1,ijk)
              DO 360 L1=1,IPHASE
                 TA0(LA,LB) = TA0(LA,LB)
     &                       +CO(NA0(NDM,L1,LB,IJK))*B0INV(LA,L1)
  360         CONTINUE
              CO(NA0(NDM,LA,LB,IJK)) = TA0(LA,LB)        
  370       CONTINUE
c
  371       IF(IJKPP(1,ijk).LT.1) GO TO 380
            DO 375 NDM=1,IJKPP(1,ijk)
              DO 372 L1=1,IPHASE
                 TC0(LA,LB) = TC0(LA,LB)
     &                       +CO(NC0(NDM,L1,LB,IJK))*B0INV(LA,L1)
  372         CONTINUE
              CO(NC0(NDM,LA,LB,IJK)) = 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(1,ijk).LT.1) GO TO 421
            DO 420 NDM=1,IJKMM(1,ijk)
              DO 410 L1=1,IPHASE
                 TA0(LA,LB) = TA0(LA,LB)
     &                       +CO(NA0(NDM,L1,LB,IJK))*B0INV(LA,L1)
  410         CONTINUE
              CO(NA0(NDM,LA,LB,IJK)) =-TA0(LA,LB)        
  420       CONTINUE
c
  421       IF(IJKPP(1,ijk).LT.1) GO TO 430
            DO 425 NDM=1,IJKPP(1,ijk)
              DO 422 L1=1,IPHASE
                 TC0(LA,LB) = TC0(LA,LB)
     &                       +CO(NC0(NDM,L1,LB,IJK))*B0INV(LA,L1)
  422         CONTINUE
              CO(NC0(NDM,LA,LB,IJK)) =-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
  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***********************************************************************

4. 질문내용(핵심내용)

       32 bit 프로세서에서 64 bit 연산을 해야 할 때, 

       그러니까 f77 compiler 에서의 -qautodbl=dblpad 옵션과 같은 기능을 주려고 하면

       ifort 에서 그냥 -autodouble 로 하면 되나요?  아님 먼가 다른 것이 필요한 건지..

5. (에러발생시) 에러 메세지 : 컴파일시 이런 워닝이 쭈욱~ 뜹니다. 일단 코드는 돌아가는 것으로 보이는데,  경고메시지를 보니 형 변환의 문제가 생기는 것 같아 불안해서요.

t2cg2.f(1115): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [DELT]
      COMMON/KONIT/KON,DELT,IGOOD
-----------------------^
t2cg2.f(1118): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [ABORT1]
     AABORT1,ABORT2
------^
t2cg2.f(1118): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [ABORT2]
     AABORT1,ABORT2
-------------^
t2cg2.f(1119): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [TIMIN]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
---------------------------------^
t2cg2.f(1119): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [SUMTIM]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
---------------------------------------^
t2cg2.f(1119): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [GF]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
----------------------------------------------^
t2cg2.f(1119): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [TIMOUT]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
-------------------------------------------------^
t2cg2.f(1679): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
   30 FORMAT('       AT ELEMENT *',A5,'*   ',8(1X,E12.6))
---------------------------------------------------^
t2cg2.f(1608): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
  114 FORMAT(7H EPS = ,E10.4,9H  RMIN = ,E10.4,10H  RESID = ,E10.4,9H  I
------------------------^
t2cg2.f(1608): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
  114 FORMAT(7H EPS = ,E10.4,9H  RMIN = ,E10.4,10H  RESID = ,E10.4,9H  I
------------------------------------------^
t2cg2.f(1608): remark #8291: Recommended relationship between field width 'W' and the number of fractional digits 'D' in this edit descriptor is 'W>=D+7'.
  114 FORMAT(7H EPS = ,E10.4,9H  RMIN = ,E10.4,10H  RESID = ,E10.4,9H  I
--------------------------------------------------------------^
t2cg2.f(1423): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [DELT]
      COMMON/KONIT/KON,DELT,IGOOD
-----------------------^
t2cg2.f(1426): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [ABORT1]
     AABORT1,ABORT2
------^
t2cg2.f(1426): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [ABORT2]
     AABORT1,ABORT2
-------------^
t2cg2.f(1427): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [TIMIN]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
---------------------------------^
t2cg2.f(1427): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [SUMTIM]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
---------------------------------------^
t2cg2.f(1427): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [GF]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT
----------------------------------------------^
t2cg2.f(1427): warning #6375: Because of COMMON, the alignment of object is inconsistent with its type   [TIMOUT]
      COMMON/CYC/KCYC,ITER,ITERC,TIMIN,SUMTIM,GF,TIMOUT

그리고 32bit 윈도우의 CVF 에서도, 64 bit 연산을 하기 위해 어떤 설정을 해야 하는지 알려주시면 감사하겠습니다.

[kernelpanic]

software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/fortran/lin/compiler_f/copts/common_options/option_m32_m64.htm원하시는게 이건지.
소스는 T2없다며 에러나네요

[Paradism]

감사합니다만, 제가 원하는 것과는 다른 것 같습니다. 제가 질문을 잘못해서 그런건지도 모르겠네요.^^ -r8 -i8 옵션 추가로 일단 linux ifort 에서는 문제가 해결되었습니다. 답변 주셔서 감사합니다.

[Paradism]

혹시 -r8 -i8 과 같은 옵션을 Compaq Visual Fortran (윈도우) 에서 설정할 수 있나요?? visual 창에서 옵션 내부 메뉴를 뒤져봤는데 찾질 못하고 있습니다. 아시는 분 도움 좀 부탁 드립니다. ^^