c!MS$OPTIMIZE:'ON:x'
c/optimize
c$LARGE
c$DEBUG
c$STORAGE:2
c      USE DFLIB
C  *****  PROGRAM AC-UNSAF2D *****************************************
C   ACADEMIC DISTRIBUTION VERSION
C     THIS IS THE PROGRAM UNSAF FOR PC/WIN95 AC-UNSAF2D(acunsf2d) 
C
C                               BY DR. M.NISHIGAKI ( OKAYAMA  UNIV. )      
C                                      Y.TAKESHITA ( OKAYAMA  UNIV. )
C
C      OCTOBER,1986  ;  JULY,1987 
C                                                                           
C   ----- FILE INFORMATION -----                                            
C                                                                           
C           * = INPUT FROM KEYBOARD & OUTPUT TO DISPLAY                    
C      FILE 2 = A(MBAND,NUMNP) & B(NUMNP) MATRIX
C      FILE 3 = CALCULATED PRESSURE HEAD DATA  P(NUMNP) OR C(NUMNP)
C      FILE 4 = FOR RECAL DATA  P,P1,Q,KODE
C      FILE 5 = ALL INPUT DATA FROM DISK FILE
C      FILE 6 = WRITE ALL RESULTS AT EACH TIME STEP                      
C      FILE 10= NODAL POINTS & ELEMENTS DATA FOR MESH CHECK AND RESULT
C      FILE 98= WORK & PLOT FILE LIST
C   -----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER*4 MTOTBS,MTOTBD,MTOTIB,ISEGM,ILENR,ILENI,
     &          N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,
     &          N11,N12,N13,N14,N15,N16,NQ16,N17,N18,N19,
     &          L1,L2,L3,L4,L5,L6,L7,LQ7,LQQ7,LQQQ7,L8,L9,L10,
     &          L11,L12,L13,L14,
     &          M1,M2,M3,M4,M5
C
      CHARACTER*4  HED
      CHARACTER*20 FNAME2,FNAME3,FNAME4,FNAME5,FNAME6,FNAM10,DDDD
C 
c      DIMENSION BS(150000),BD(250000),IB(150000)
C
      REAL*8 , ALLOCATABLE::BS(:)  !　配列の大きさを変更できるような宣言
      REAL*8 , ALLOCATABLE::BD(:) 
      INTEGER , ALLOCATABLE::IB(:)
C   
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /FIX/ HED(20),INIT,KTYPE,MAXIT,KAT,INTEG,NSEEP,NCUL,NPRINT,
     &             MAXITO,NTYPE,JMES
      COMMON/WELL1/NB(20),NE(20),NBW,NPT,NPT1,NPB,NPB1,LW,NEP,NSW,KCON
	COMMON/WELL2/AL,AL1,QP,QA,RW,RT,ALFA
      COMMON /MESHZU/ MESH,INCOR,IVFLAG,NT,NUMEP,ITQ,NBOUN,NBNQ
      COMMON /CRTL1/LA,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL2/DT,DT1
	COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C -------------------------------------------------------------------     
      DATA MTOTBS/1500000/,MTOTBD/2500000/,MTOTIB/1500000/
	DATA  ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN/98,2,3,4,5,8,10,20/
cc      result=SETEXITQQ(QWIN$EXITNOPERSIST)
      ISEGM = 65536
C
C--------  COMMENT ABOUT VERSION ------------------------------------
C
      WRITE(*,8000)
 8000 FORMAT(//' THIS PROGRAM WAS DEVELOPED BY'
     & /'    DR. NISHIGAKI   ( ORIGINAL, MAIN-FRAME )'
     & /'    DR. TAKESHITA   ( PC VERSION)   1986    '
     & /'                             (V.2.1 ) 1987.07'   
     & /'                             (V.2.11) 1988.02'        
     & /'                             (V.2.11) 1988.06'            
     & /'                       INTEGERLIZE(1990.11.25)' 
     & /'   BY Y.SHINSHI     COMPLETELY REVISED VERSION 1996.01'
     & /'                    ACADEMIC VERSION 1998.09' 
     & /'   BY T.IMAI                 VERSION 2005.10')
C 
C ================================================== READ OPEN FILE NAME 
C ------------------------------- READ CONTOROL DATA FROM No.5----------

      ILENR = ISEGM/8
      ILENI = ISEGM/2
      WRITE (*,100)
      READ (*,20) FNAME5
C-------------------------------- ALL RESULTS WRITE TO FILE No.6 --------
      WRITE (*,101)
      READ (*,20) FNAME6
C--------------------------------------------------- READ OPEN FILE NAME  
      OPEN (ND98,FILE='WORK.LST',STATUS='OLD',FORM='FORMATTED')
      READ(ND98,20) FNAME2
      READ(ND98,20) FNAME3
      READ(ND98,20) FNAME4
      READ(ND98,20) FNAM10
C
      MTOTBS2  = 0
      MTOTBD2  = 0
      MTOTIB2  = 0
C
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
      READ (ND98,20,END=999) DDDD
C
      READ (ND98,*,END=999) MTOTBS2
      READ (ND98,*,END=999) MTOTBD2 
      READ (ND98,*,END=999) MTOTIB2 
  999 CONTINUE
	IF(MTOTBS2  .LE. MTOTBS)  MTOTBS2 =MTOTBS 
	IF(MTOTBD2  .LE. MTOTBD)  MTOTBD2 =MTOTBD 
	IF(MTOTIB2  .LE. MTOTIB)  MTOTIB2 =MTOTIB 
	 MTOTBS =MTOTBS2 
	 MTOTBD =MTOTBD2 
	 MTOTIB =MTOTIB2  
  997 CONTINUE
      CLOSE(ND98)
      ALLOCATE( BS(MTOTBS))         !　配列の宣言
      ALLOCATE( BD(MTOTBD)) 
      ALLOCATE( IB(MTOTIB))
C=================================================== FILE OPEN =========
      OPEN (ND2, FILE=FNAME2, STATUS='unknown',FORM='UNFORMATTED')
      OPEN (ND3, FILE=FNAME3, STATUS='unknown',FORM='UNFORMATTED')
      OPEN (ND4, FILE=FNAME4, STATUS='unknown',FORM='UNFORMATTED')
      OPEN (ND5, FILE=FNAME5, STATUS='OLD',FORM='FORMATTED')
      OPEN (ND6, FILE=FNAME6, STATUS='unknown',FORM='FORMATTED')
      OPEN (ND10,FILE=FNAM10, STATUS='unknown',FORM='UNFORMATTED')
C-------------------------------------------------------------------------
      READ(ND5, 25)  HED
C      
      READ(ND5, 30 )  NUMNP,NUMEL,KAT,KCON,NUMMAT,INIT,KTYPE,MAXIT,
     &             INTEG,NSEEP,MAXS,NBOUN,MNBN,MIHT
C     
      READ(ND5, 30 )  NBW,NPB,NPT,NEP,LW
C                                        
      READ(ND5, 30 )  NUMEP,ITQ,NBNQ,MNBQ,MIQT                                                
C
      READ(ND5, 30)  NCUL,NT,MAXITO,NTYPE,NSEEP1,MAXS1,IVFLAG,
     &             MESH,INCOR,JMES,IRELM,IRNDE,NSOLV,NDSPL
C     
      IF( INCOR .EQ. 0 .AND. ILENR.LT. NUMNP ) THEN
       WRITE(*,*)' YOU CANNOT USE OUTER DISK SYSTEM'
       WRITE(*,*)' YOU CAN HAVE 2 WAYS AS FOLLOW'
       WRITE(*,*)' WAY-1 TO SET [INCOR=1]'
       WRITE(*,*)' WAY-2 TO ENLARGE WRITING SEGMENT CAPACITY'
       WRITE(*,*)'                     NOW IT IS',ISEGM,'B.'
       STOP
      END IF
C------------------------------------------------------------------------------
      IF(NDSPL.EQ.0) OPEN (UNIT=NDCN,FILE='CON')
      IF(NDSPL.NE.0) OPEN (UNIT=NDCN,FILE='NUL')
C------------------------------------------------------------------------------
      IF (KTYPE.NE.-1) GO TO 70
        MK=1
        MP=1
        GO TO 75
   70 READ(ND5, 30) (NUMK(J),NUMP(J),J=1,NUMMAT)                              
        MK=NUMK(1)                                                            
        MP=NUMP(1)                                                            
       IF (NUMMAT .NE. 1) THEN                                            
      DO J=2,NUMMAT                                                      
       IF (NUMK(J) .GT. MK) MK=NUMK(J)                                       
       IF (NUMP(J) .GT. MP) MP=NUMP(J)                                       
      END DO                                                              
       ENDIF
C --------------------------------------------------------
   75 NT1 = NT                                                              
       IF(NT1 .LT. 1) NT1=1                                                  
       IF(NCUL .LT.1) NCUL=3                                                 
       IF(MAXITO .LT.1) MAXITO=20                                           
      NUMDP=NUMEP                                                           
      ITQ1=ITQ                                                              
       IF(NUMDP .LT. 1) NUMDP=1                                             
       IF(ITQ .LT. 1) ITQ1=1                                                 
      NBOUN1=NBOUN                                                          
      MNBN1=MNBN                                                            
      MIHT1=MIHT                                                            
       IF(NBOUN .LT.1) NBOUN1=1                                              
       IF(MNBN .LT. 1) MNBN1=1                                              
       IF(MIHT .LT. 1) MIHT1=1  
C
       NBNQ1= NBNQ
	 MNBQ1= MNBQ
	 MIQT1= MIQT
	 IF(NBNQ1.LT.1) NBNQ1=1
	 IF(MNBQ1.LT.1) MNBQ1=1
	 IF(MIQT1.LT.1) MIQT1=1

C	                                              
         MAXSP=1                                                               
         NDIM=1                                                                
       IF (NSEEP .GE. 1) THEN                                            
         NDIM=NSEEP                                                            
         MAXSP=MAXS
	 END IF                                                            
   40  IF(NSEEP1 .GE.  1)  THEN                                         
       IF(NSEEP1 .GT. NSEEP) NDIM=NSEEP1                                     
       IF(MAXS1 .GT. MAXS) MAXSP=MAXS1
      END IF 
C
      INUMNP=1
       IF( INCOR .EQ. 1 ) INUMNP=NUMNP
      NUMNPZ=1
C       IF( KAT.LE.0) NUMNPZ=NUMNP                                
C --------------------------------   DYNAMIC  ALLOCATION --------------
      N1  = 1                                                           Q
      N2  = N1  + NUMNP                                                 P
      N3  = N2  + NUMNP                                                 P1
      N4  = N3  + NUMNP                                                 TDCRC
      N5  = N4  + NUMNP*3                                               SCANG 
      N6  = N5  + NUMEL*2                                               XY   
      N7  = N6  + NUMNP*2                                               VXY
      N8  = N7  + NUMEL*2                                               H
      N9  = N8  + NUMNP                                                 CPI
      N10 = N9  + INUMNP*5                                              XYK
      N11 = N10 + MK*NUMMAT*2                                           XYP
      N12 = N11 + MP*NUMMAT*2                                           XYC
      N13 = N12 + MP*NUMMAT*2                                           WIDBTP
      N14 = N13 + NUMDP*2                                               XYQ
      N15 = N14 + ITQ1*2                                                TM
      N16 = N15 + NT1                                                   XYH
      NQ16 = N16 + MIHT1*NBOUN1*2                                        XYQT 
      N17 = NQ16 + MIQT1*NBNQ1*2                                         C
      N18 = N17 + NUMNP                                                 Z
      N19B= N18 + NUMNPZ  
      N19 = N19B+ NDIM                                                  ! SEEPVAL
C-------------------------------------------------------
      L1  = 1                                                           KODE
      L2  = L1  + NUMNP                                                 KX
      L3  = L2  + NUMEL*5                                               NP
      L4  = L3  + NDIM*MAXSP                                            NSPKDS
      L5  = L4  + NDIM*2                                                NBN
      L6  = L5  + NBOUN1                                                NBOU 
      L7  = L6  + NBOUN1*MNBN1                                          NIHT
      LQ7 = L7  + NBOUN1                                                NBQ
      LQQ7= LQ7 + NBNQ1                                                 NBOUQ
      LQQQ7=LQQ7 + NBNQ1*MNBQ1                                          NIQT 
      L8  = LQQQ7  + NBNQ1                                              NSIC
      L9  = L8  + NUMDP                                                 IC
      L10 = L9  + NUMEL                                                 KPM
      L11 = L10 + NUMNP                                                 KODEW
	L12 = L11 + INUMNP                                               ! IS
C-------------------------------------------------------
C
      M1=1  
       WRITE (NDCN,80) N19-1,M1-1,L12-1
       WRITE (ND6,80) N19-1,M1-1,L12-1
      IF( N19-1 .GT. MTOTBS .OR. M1-1 .GT. MTOTBD 
     &    .OR. L12-1 .GT. MTOTIB ) THEN       
       WRITE (NDCN,90)
       STOP
      ELSE
C**********************************************************************
       CALL MAFILL ( IB(L2),BS(N5),BS(N5+NUMEL),NUMEL,MBAND,IRELM ) 
      END IF
C**********************************************************************      
      NUMNP1=1
      NWK = MBAND*NUMNP
      NWK1 = 1
       WRITE(NDCN,*) '  NSOLV = ', NSOLV
       IF( NSOLV .EQ. 1 ) NUMNP1=NUMNP+1
       L13 = L12 + NUMNP1                                               ISP 
       L14 = L13 +NWK1
       M1 = 1
       M2 = M1 + NUMEL*2 
       WRITE (NDCN,80) N19-1,M2-1,L13-1
       WRITE (ND6,80) N19-1,M2-1,L13-1
C
       IF( N19-1 .GT. MTOTBS .OR. M2-1 .GT. MTOTBD 
     &    .OR. L13-1 .GT. MTOTIB ) THEN       
        WRITE (NDCN,90)
        STOP
       END IF
       IF( NSOLV .EQ. 1 ) THEN
C----------------------------------------------------------------------      
          CALL COLHT ( IB(L2),IB(L12),BD(M1),NUMNP,NUMNP1,NUMEL,NWK )
        NWK1 = NWK
        L14 = L13+NWK1 
        WRITE (NDCN,80) N19-1,M2-1,L14-1
        WRITE (ND6,80) N19-1,M2-1,L14-1
       IF( N19-1 .GT. MTOTBS .OR. M2-1 .GT. MTOTBD 
     &    .OR. L14-1 .GT. MTOTIB ) THEN       
          WRITE (NDCN,90)
          STOP
       ELSE              
        CALL ADRSK(IB(L2),IB(L12),IB(L13),BD(M1),NUMNP,NUMNP1,NUMEL,NWK)
       END IF
       END IF
C      
      M1  = 1                                                           A
      M2  = M1  + NWK                                                   B
      M3  = M2  + NUMNP                                                 AW     
C------------------------------------------- FOR IN CORE
      IMBAND=1
      INWK=1
       IF( INCOR .EQ. 1 .AND. NSOLV .EQ. 0 ) INWK=MBAND*NUMNP
       IF( INCOR .EQ. 1 .AND. NSOLV .EQ. 0 ) IMBAND=MBAND
       IF( INCOR .EQ. 1 .AND. NSOLV .EQ. 1 ) INWK=NWK 
       M4  = M3  + INWK                                                 PS
       M5  = M4  + NUMNP
C ------------------------------------------------- MEMEORY CHECK
      WRITE (NDCN,80) N19-1,M5-1,L14-1
      WRITE (ND6,80) N19-1,M5-1,L14-1
      IF( N19-1 .GT. MTOTBS .OR. M5-1 .GT. MTOTBD 
     &    .OR. L14-1 .GT. MTOTIB  ) THEN       
       WRITE (NDCN,90)
      ELSE
      WRITE(NDCN,*) 'N',N1,N2,N3,N4,N5,N6,N7,N8,
     &          N9,N10,N11,N12,
     &          N13,N14,N15,N16,NQ16,N17,N18,N19
      WRITE(NDCN,*) 'L',L1,L2,L3,L4,L5,L6,L7,LQ7,LQQ7,LQQQ7,L8,L9,L10,
     &          L11,L12B,L12,L13,L14
      WRITE(NDCN,*) 'M',M1,M2,M3,M4,M5
C      STOP
      CALL FEM (BD(M1),BD(M2),BD(M3),BD(M4),
     &          BS(N1),BS(N2),BS(N3),BS(N4),BS(N5),BS(N6),BS(N7),BS(N8),
     &          BS(N9),BS(N10),BS(N11),BS(N12),BS(N13),BS(N14),BS(N15),
     &          BS(N16),BS(NQ16),BS(N17),BS(N18),BS(N19B),
     &          IB(L1),IB(L2),IB(L3),IB(L4),IB(L5),IB(L6),IB(L7),
     &         IB(LQ7),IB(LQQ7),IB(LQQQ7),
     &    IB(L8),IB(L9),IB(L10),IB(L11),IB(L12),IB(L13),
     & NUMNP,MBAND,NWK,NUMEL,NUMMAT,MK,MP,NDIM,NT1,MAXSP,NUMDP,
     & NBOUN1,ITQ1,MIHT1,MNBN1,
     & NBNQ1,MNBQ1,MIQT1,
     & INUMNP,INWK,IRNDE,NUMNP1,NWK1,NUMNPZ)
      ENDIF
      WRITE (NDCN,110)
      WRITE (ND6,110)
   20 FORMAT (A)                                                         
   25 FORMAT (20A4)                                                         
   30 FORMAT (14I5)                                                         
   80 FORMAT (/,' *** STORAGE REQUIRD',' BS= ',I8,' WORDS',' BD= ',I10,/   
     &          ' WORDS','  IB= ',I8',WORDS'//)                        
   90 FORMAT (//,'///// ERROR ///// MEMORY OVER ////// ')
  100 FORMAT (//,'== INPUT DATA FILE NAME [ EX. B:DATA1.DAT ] ====> ')
  101 FORMAT (//,'===    OUTPUT FILE NAME [ EX. C:UNS1.PRN  ] ====> ')
  110 FORMAT (//,'   <<<<<  OK! OK!  NORMAL TERMINATION  >>>>> ')
      STOP                                                                  
      END                                                                   
      SUBROUTINE FEM (A,B,AW,PS,
     &                Q,P,P1,TDCRC,SCANG,XY,VXY,H,CPI,
     &         XYK,XYP,XYC,WIDBTP,XYQ,TM,XYH,XYQT,C,Z,SEEPVAL,
     &                KODE,KX,NP,NSPKDS,NBN,NBOU,NIHT,
     &                                  NBQ,NBOUQ,NIQT,
     &  NSIC,IC,KPM,KODEW,IS,ISP,
     & NUMNP,MBAND,NWK,NUMEL,NUMMAT,MK,MP,NDIM,NT1,MAXSP,NUMDP,
     & NBOUN1,ITQ1,MIHT1,MNBN1,
     & NBNQ1,MNBQ1,MIQT1,
     & INUMNP,INWK,IRNDE,NUMNP1,NWK1,NUMNPZ)
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*4  HED
      INTEGER*4 ILENR,ILENI
C      
      DIMENSION A(NWK),B(NUMNP),KODE(NUMNP),Q(NUMNP),P(NUMNP),
     & P1(NUMNP),TDCRC(NUMNP,3),XY(NUMNP,2),KX(NUMEL,5),SCANG(NUMEL,2),
     & XYK(MK,NUMMAT,2),XYP (MP,NUMMAT,2),XYC(MP,NUMMAT,2),
     & NP(NDIM,MAXSP),NSPKDS(NDIM,2),WIDBTP(NUMDP,2),XYQ(ITQ1,2),
     & NBN(NBOUN1),NBOU(NBOUN1,MNBN1),XYH(MIHT1,NBOUN1,2),
     & NIHT(NBOUN1),IC(NUMEL),VXY(NUMEL,2),NSIC(NUMDP),TM(NT1),
     & CPI(INUMNP,5),H(NUMNP),KPM(NUMNP), IS(NUMNP1),ISP(NWK1),
     & KODEW(INUMNP),AW(INWK),PS(NUMNP),C(NUMNP),Z(NUMNPZ)
      DIMENSION NBQ(NBNQ1),NBOUQ(NBNQ1,MNBQ1),NIQT(NBNQ1),
     &          XYQT(MIQT1,NBNQ1,2) ,SEEPVAL(NDIM)
C
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /FIX/ HED(20),INIT,KTYPE,MAXIT,KAT,INTEG,NSEEP,NCUL,NPRINT,
     &             MAXITO,NTYPE,JMES
      COMMON/WELL1/NB(20),NE(20),NBW,NPT,NPT1,NPB,NPB1,LW,NEP,NSW,KCON
	COMMON/WELL2/AL,AL1,QP,QA,RW,RT,ALFA
      COMMON /MESHZU/ MESH,INCOR,IVFLAG,NT,NUMEP,ITQ,NBOUN,NBNQ
C      COMMON /CRTL/LA,DT,DT1,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL1/LA,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL2/DT,DT1

	COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN

CC ----INITIALIZATION
      IERREND=0
      EPSLON=0.0D0
      NERR=0
C------------------------------------------- READ SEEPAGE FACE INFORMATION
      IF (NSEEP .LT. 1) GO TO 30
      IF (MAXIT .LT. 3) MAXIT=3
        DO I=1,NSEEP
         READ(ND5, 122) NSPKDS(I,1),SEEPVAL(I)
122      FORMAT(I5,F10.0) 
         NSPKDS(I,2)=-1
           NS=NSPKDS(I,1)
           READ(ND5, 10) (NP(I,J),J=1,NS)
        END DO
   10 FORMAT (14I5)

C---------------------------------- READ POTENTIAL SURFACE FLUX INFORMATION
   30 IF(NUMEP .LT.1) GO TO 45
       READ(ND5, 90) EI,PL
       READ(ND5, 10)(NSIC(I),I=1,NUMEP)
       BTPI=1.0D-1
       IF(ITQ .LT. 1) GO TO 35
        READ(ND5, 90)(XYQ(I,1),XYQ(I,2),I=1,ITQ)
   35   DO I=1,NUMEP
C---------------------------------------
         WIDBTP(I,2)=BTPI
C---------------------------------------
        END DO
C----------------------------------------- READ VARIABLE BOUNDARY CONDITION
   45 IF(NBOUN .LT.1 ) GO TO 71
       DO I=1,NBOUN
        READ(ND5, 10) NBN(I)
        NB1=NBN(I)
        READ(ND5, 10) (NBOU(I,J), J=1,NB1)
       END DO
       READ(ND5, 10)(NIHT(J),J=1,NBOUN)
        DO M=1,NBOUN
         NIHTN=NIHT(M)
         READ(ND5, 90)(XYH(I,M,1),XYH(I,M,2), I=1,NIHTN)
        END DO
C---------------------------------READ VARIABLE DISCHARGE RATE BOUNDARY CONDITION
   71 CONTINUE
      IF(NBNQ .LT.1 ) GO TO 70
       DO I=1,NBNQ
        READ(ND5, 10) NBQ(I)
        NBQ1=NBQ(I)
        READ(ND5, 10) (NBOUQ(I,J), J=1,NBQ1)
       END DO
C
       READ(ND5, 10)(NIQT(J),J=1,NBNQ)
        DO M=1,NBNQ
         NIQTN=NIQT(M)
        READ(ND5,90)(XYQT(I,M,1),XYQT(I,M,2), I=1,NIQTN)
        END DO
C------------------------------------------------- READ WELLBORE INFORMATION
   70 AL=0.0D0
      QA=0.0D0
      IF (NBW .LT. 1) GO TO 80
       NSW=NPT-NPB
       NPB1=NPB+1
       NPT1=NPT-1
        READ(ND5, 10) (NB(I),I=1,NBW),(NE(I),I=1,NEP)
        READ(ND5, 90) RW,RT,QP,AL,ALFA
        IF (ALFA .LE. 0.) ALFA=0.8D0
       AL1=AL
C------------------------------------------------------- GENERAL INFORMATION
   80 READ(ND5, 92) DT,DTMAX,DMUL,TMAX,TOL,TSOLV,ITMAX
      IF( TSOLV .LE. 0.0 ) TSOLV=1.0D-8
      IF( ITMAX .LE. 0 ) ITMAX=5
         ITMAX=NUMNP*ITMAX
      IF( NT. LE. 0) GO TO 95
      WRITE(*,*) 'NT',NT
      READ(ND5, 91) (TM(I),I=1,NT)

   90 FORMAT (6F10.0)
   91 FORMAT (7F10.0)
   92 FORMAT (6F10.0,I5)
C---------------------------------------------------------
   95 IF (KAT .EQ. 0) WRITE (ND6, 100)
      IF (KAT .EQ. 1) WRITE (ND6, 110)
      IF (KAT .EQ. 2) WRITE (ND6, 120)
      IF (KAT .EQ. 3) WRITE (ND6, 130)
C---------------------------------------------------------
           IF (KAT .EQ. 3) KAT=1
C---------------------------------------------------------
      IF (NTYPE .EQ. -2) WRITE (ND6, 1300)
      IF (NTYPE .EQ. -1) WRITE (ND6, 1301)
      IF (NTYPE .EQ.  0) WRITE (ND6, 1302)
      IF (NTYPE .EQ.  1) WRITE (ND6, 1303)
      IF ( KCON .LE. 0) WRITE (ND6, 137)
      IF ( KCON .EQ. 1) WRITE (ND6, 138)
  100 FORMAT (///,22H HORIZONTAL PLANE FLOW/)
  110 FORMAT (///,18H AXISYMMETRIC FLOW/)
  120 FORMAT (///,20H VERTICAL PLANE FLOW/)
  130 FORMAT (///,54H AXISYMMETRIC FLOW TO A WELL WITH PRESCRIBED DISCHA
     &RGE/)
 1300 FORMAT (//,'***** NONSTEADY ----> FINAL STEADY  PROBLEM *****'/)
 1301 FORMAT (//,'***** ONLY NONSTEADY PROBLEM *******'/)
 1302 FORMAT (//,'***** INITIAL STEADY ---> NONSTEADY PROBLEM *****'/)
 1303 FORMAT (//,'***** ONLY STEADY PROBLEM *******'/)
  137 FORMAT (///,'//// UNCONFIND or nonlinear GROUND WATER FLOW ////'/)
  138 FORMAT (///,'////// CONFIND or linear GROUND WATER FLOW //// '/)
      WRITE (ND6, 140) HED,NUMNP,MBAND,NUMEL,NUMMAT,MAXIT,MAXITO,JMES,
     &               NSOLV
  140 FORMAT (///,20A4,//30H NUMBER OF NODAL POINTS------- I4/
     &30H BAND WIDTH ------------------ I4/
     &30H NUMBER OF ELEMENTS ---------- I4/
     &30H NUMBER OF MATERIALS --------- I4/
     &30H MAX NUMBER OF ITERATIONS ---- I4/
     &30H MAX NUMBER OF ITERATIONS(T=0) I4/
     &30H ELEMENT TYPE TRIAN/0 IPARTA/1 I4/
     &30H SOLVER CHOICE GE/0 PCG/1      I4/)
C------------------------------------------------------------------
      WRITE (ND6,143)DT,DTMAX,DMUL,TMAX,NT,TOL
  143 FORMAT(/' INITIAL TIME----------- '1PE10.3/
     &        ' MAXIMUM TIME----------- '1PE10.3/
     &        ' MULTIPLY FACTOR-------- '1PE10.3/
     &        ' FINAL TIME------------- '1PE10.3/
     &        ' PRINT OUT STEP-(NT)---- 'I5/
     &        ' MINIMUM ERROR---------- '1PE10.3/)
C-------------------------------------------------------------------
      IF (NT .LT. 1) GO TO 4444
      WRITE(ND6,681) (TM(I),I=1,NT)
  681 FORMAT(/,' *** PRINT OUT TIME INFORMATION ***',
     & /,(10X,5(1PE20.4))/ )
 4444 CONTINUE
C
      IF(NBOUN .LT. 1) GO TO 156
      DO 154 I=1,NBOUN
      WRITE (ND6, 150)I
  150 FORMAT(/' NODAL POINTS IN BOUNDARY NUMBER=',I5/
     1         ' TABLE OF TIME VERSUS BOUNDARY HEAD'/)
      NIHTN=NIHT(I)
      WRITE (ND6, 151)(XYH(J,I,1),XYH(J,I,2),J=1,NIHTN)
      WRITE (ND6, 152)
      NB1=NBN(I)
      WRITE (ND6, 153) ( NBOU(I,J),J=1,NB1)
  151 FORMAT(10E12.4)
  152 FORMAT(/'  NODAL POINTS IN THIS BOUNDARY '/)
  153 FORMAT(24I5)
  154 CONTINUE
C
  156 CONTINUE
        IF(NBNQ .LT. 1) GO TO 155
      DO I=1,NBNQ
      WRITE (ND6, 157)I
  157 FORMAT(/' NODAL POINTS IN DISCHARGE RATE BOUNDARY NUMBER=',I5/
     1         ' TABLE OF TIME VERSUS DISCHARGE BOUNDARY HEAD'/)
      NIQTN=NIQT(I)
      WRITE (ND6, 151)(XYQT(J,I,1),XYQT(J,I,2),J=1,NIQTN)
      WRITE (ND6, 152)
      NBQ1=NBQ(I)
      WRITE (ND6, 153) ( NBOUQ(I,J),J=1,NBQ1)
      END DO

C
  155 IF(NUMEP .LT. 1) GO TO 165
      WRITE (ND6, 161) EI,PL
      IF(ITQ .LT. 1) GO TO 165
      WRITE (ND6, 163) (XYQ(I,1),XYQ(I,2),I=1,ITQ)
  161 FORMAT(' MAXIMUM ALLOWED RATE OF INFILTRATION=',1PE11.3/
     1       ' MAXIMUM ALLOWED PRESSURE HEAD=',1PE11.3/)
  163 FORMAT(/' TABLE OF INFILTRATION VERSURS TIME '/10(1PE12.4))
  165 IF (NBW .LT. 1) GO TO 200
      WRITE (ND6, 170) NPB,NPT,RW,RT,AL
  170 FORMAT (/,22H BTTOM NODE OF WELL  =I4/
     1        22H TOP NODE OF WELL--- =I4/
     2        22H RADIUS OF WELL----- =1PE11.3/
     3        22H RADIUS OF TUBING--- =1PE11.3/
     4        22H INITIAL WATER LEVEL =1PE11.3//)
      AA=.5*(RW*RW-RT*RT)
      RW=AA
C------------------------------------ READ AND GENERATE MATERIAL INFORMATION
  200 CALL MATIN (XYK(1,1,1),XYK(1,1,2),XYP(1,1,1),XYP(1,1,2),XYC(1,1,1)
     &           ,XYC(1,1,2),MK,MP,NUMMAT,KAT,KTYPE)
C--------------------------------- READ AND GENERATE NODAL POINT INFORMATION
      CALL NPIN (KODE,Q,P,P1,XY(1,1),XY(1,2),WIDBTP(1,1),NSIC,KX,B,Z,
     &    NUMDP,NUMNP,INIT,NERR,NUMEP,KAT,NUMEL,IRNDE,NUMNPZ)
      IF (NERR .GT. 0) THEN 
      WRITE (NDCN,700) 
  700 FORMAT (///,' <<<<<<<< NODAL POINT INFORMATION ERROR >>>>>>>>> ')
C
       WRITE(110,438)  !ERROR   
C
      RETURN
      ENDIF
C---------------------------------------------------------------------------
C     IF(JMES.EQ.2 ) CALL INFLUE (XY(1,1),XY(1,2),KX,IC,NUMNP,NUMEL)
      WRITE (NDCN,710)
  710 FORMAT (////'*** OK! OK! ** ALL INPUT DATA HAVE BEEN READ ******')
C------------------------------------ WRITE(10) MESH CHECK DATA -------
      WRITE(ND10) HED
      WRITE(ND10) NUMNP
      WRITE(ND10) (XY(I,1),XY(I,2),I=1,NUMNP)
      WRITE(ND10) NUMEL
      WRITE(ND10) (KX(I,1),KX(I,2),KX(I,3),KX(I,4),
     &                   KX(I,5),SCANG(I,1),SCANG(I,2),I=1,NUMEL)
C
      IERSW  = 0
      CALL ECHECK ( XY(1,1),XY(1,2),KX,NUMEL,NUMNP,IERSW  )
C
      IF(IERSW  .EQ. 1) GO TO 3221
      IF( MESH.LE.0) GO TO 220 
      WRITE (NDCN,1000)
      WRITE (ND6,1000)
 1000 FORMAT(' **** MESH CHECK ONLY!! NORMAL TERMINATION *****')
 3221 CONTINUE
      LADD = 0
      TIMEDD=0.0
      WRITE(ND10) LADD,TIMEDD
      WRITE(ND10) (KODE(I),P(I),Q(I),I=1,NUMNP)
      WRITE(ND10) (VXY(I,1),VXY(I,2),I=1,NUMEL)
      RETURN
C----------------------------------------- SET UP FOR TIME INTEGERATION
  220 IFLGTM=1
       LA=0
       TIME=0.0D0
      IF( NTYPE .LT. 0) THEN
        LA = 1
        TIME=DT
	END IF
  221 QTOT=0.0D0
      ITER=0
      LAA=0
C      NN=1
  230	CONTINUE
       WRITE (NDCN, 239)  HED,TIME,LA
       WRITE (ND6, 240)  HED,TIME,LA,AL
C
       WRITE(110,435)TIME,LA
 435  FORMAT("T,",E14.6,",",I10)
C
  239 FORMAT (//20A4,//10X,'TIME =',1PE13.5,5X,'TIME STEP =',I8/)
  240 FORMAT (20A4,/6HTIME =1PE13.5,5X,11HTIME STEP =I8,'  AL=',1PE10.3)
C================================================================ FOR FILE 4 
      IF( INCOR .EQ. 0 ) THEN
       IF(LAA .LT.1) GO TO 1
       GO TO 2
    1  REWIND ND4
       WRITE(ND4) P
       WRITE(ND4) P1
       WRITE(ND4) Q
       WRITE(ND4) KODE,AL
       GO TO 3
    2  REWIND ND4
       READ(ND4) P
       READ(ND4) P1
       READ(ND4) Q
       READ(ND4) KODE,AL
      ELSE
C================================================ FOR IN CORE =============
 1111  ALW1=AL
       IF(LAA .LT.1) GO TO 1112
       GO TO 2222
 1112  DO 801 I=1,INUMNP
        CPI(I,1)=P(I)
        CPI(I,2)=P1(I)
        CPI(I,3)=Q(I)
        KODEW(I)=KODE(I)
  801  CONTINUE
       GO TO 3
 2222   AL=ALW1
       DO 802 I=1,INUMNP
        P(I)=CPI(I,1)
        P1(I)=CPI(I,2)
        Q(I)=CPI(I,3)
        KODE(I)=KODEW(I)
  802  CONTINUE
      END IF
      WRITE (NDCN, 5) LAA
    5 FORMAT(/20X,'RECUL=',I5/)
C==========================================================================
    3 EPSLN1=1.0D28
      EPSLN2=1.0D28
      EPSLON=1.0D28
C---------------------------------------- PERTURB WATER LEVEL IN WELL
      IF( LA. LE. 0) GO TO 241  !250
      IF (NBW .GE. 1) CALL WELBOR (KODE,Q,P,XY(1,2),NUMNP,DT)
C--------------------------------------- SET UP NEW INFILTRATION RATE
c  241 CONTINUE
c      IF(NUMEP .LT.1) GO TO 242
c      IF( ITQ .LT. 1) GO TO 246
c      TIM1=TIME-DT*0.5D00
c      CALL INFIL (XYQ(1,1),XYQ(1,2),TIM1,QTOR,ITQ1)
c      EI=QTOR
c  246 DO I=1,NUMEP
c       N=NSIC(I)
c       IF( KODE(N) .EQ. -4)THEN
c         Q(N)=EI*WIDBTP(I,1)
c	 END IF
c      END DO
C--------------------------------------- SET UP NEW BOUNDARY CONDITION
c  242 continue
      IF( NBOUN .LT. 1) GO TO 247
C-----
      TI=TIME-DT*0.5D00
C-----
      DO M=1,NBOUN
       CALL INTERP ( XYH(1,1,1),XYH(1,1,2),TI,TP,NIHT(M),NBOUN1,M,MIHT1)
       NB1=NBN(M)
       DO J=1,NB1
        N=NBOU(M,J)
        IF(XY(N,2) .LE. TP) THEN
          P(N)=TP-XY(N,2)
          KODE(N)=3
        ELSE
  243    KODE(N)=-2
         IF( NUMEP .GE. 1 ) KODE(N)=-4
           Q(N) = 0.0  !!!!!!!!
	  END IF
  244  END DO
      END DO
C
  247 CONTINUE  
C
        IF( NBNQ .LT. 1) GO TO 241
C-----
      TI=TIME-DT*0.5D00
C-----
      DO M=1,NBNQ
       CALL INTERP( XYQT(1,1,1),XYQT(1,1,2),TI,TQ,NIQT(M),NBNQ1,M,MIQT1)
       NBQ1=NBQ(M)
       DO J=1,NBQ1
        N=NBOUQ(M,J)
          Q(N)=TQ
	 END DO
      END DO   
C
  241 continue
      IF(NUMEP .LT.1) GO TO 242
      IF(LA .LE. 0) GO TO 246  !!
      IF( ITQ .LT. 1) GO TO 246
      TIM1=TIME-DT*0.5D00
      CALL INFIL (XYQ(1,1),XYQ(1,2),TIM1,QTOR,ITQ1)
      EI=QTOR
  246 DO I=1,NUMEP
       N=NSIC(I)
       IF( KODE(N) .EQ. -4)THEN
         Q(N)=EI*WIDBTP(I,1)
	 END IF
      END DO
  242 continue
C---------------------------------- GENERATE TERMS OF MATRIX EQUATION
  250 CONTINUE
       IF( NSOLV .EQ. 1 ) THEN
        DO I= 1, NUMNP
         PS(I)=P(I)
        END DO
         END IF 
C       WRITE(NDCN,*)"--- ENTER RESET"    
      CALL RESET(KODE,A,B,Q,P,P1,TDCRC(1,1),TDCRC(1,2),TDCRC(1,3),
     &           C,XY(1,1),XY(1,2),KX,SCANG(1,1),SCANG(1,2),
     &           XYK(1,1,1),XYK(1,1,2),XYP(1,1,1),XYP(1,1,2),
     &           XYC(1,1,1),XYC(1,1,2),AW,CPI,IC,IS,ISP,Z,
     & NUMNP,MBAND,NWK,NUMEL,NUMMAT,MK,MP,INUMNP,
     & INWK,NUMNP1,NWK1,NUMNPZ)
C       WRITE(NDCN,*)"--- OUT RESET"    
C---------------------------- MODIFY FOR CONSTANT PRESSURE HEAD NODES
      IF( NSOLV .EQ. 0 ) CALL CONSTP (A,B,KODE,P,NUMNP,MBAND)
C------------------------------------------ SOLVE BY GAUSS ELIMINATION
      IF( NSOLV .EQ. 0 ) THEN 
        CALL SOLVE (A,B,NUMNP,MBAND)
        DO I=1,NUMNP
         P(I)=B(I)
        END DO
       CALL FIXQ (A,B,KODE,Q,P,P1,AW,CPI,
     1          NUMNP,MBAND,INTEG,INCOR,INUMNP,NWK,INWK)
      ELSE 
        DO I=1,NUMNP
         P(I)=PS(I)
        END DO 
C        WRITE(NDCN,*)'--- START SOLVCG '
        CALL SOLVCG (A,B,P,ISP,IS,TDCRC(1,1),TDCRC(1,2),KODE,TDCRC(1,3),
     &   C,NUMNP,NUMNP1,  NWK,TSOLV,LA,TIME,ITMAX)    
        CALL FIXQPCG ( A,B,KODE,Q,P,P1,ISP,IS,AW,CPI,
     &                  NUMNP,NUMNP1,NWK,INTEG,
     &                  INWK,INUMNP,INCOR )
       END IF
C--------------------------------------- MODIFY CONDITIONS ON WELLBORE
      IF (NBW .LE. 0) GO TO 1987 
        CALL WELBOR (KODE,Q,P,XY(1,2),NUMNP,DT)
      GO TO 1988
C---------------- ITERATE TO IMPROVE SOLUTION  ( FOR CONFIND AQUIFER )
 1987 IF (MAXIT .LT.1 .OR. KCON .GT. 0)  GO TO 460
 1988 ITER=ITER+1
      WRITE (NDCN,432) ITER
C--------------------------------- MODIFY CONDITIONS ON SEEPAGE FACES
      IF (NSEEP .LE. 0) GO TO 320
      DO 310 I=1,NSEEP
      ICHECK=0
      NS=NSPKDS(I,1)
      DO 300 J=1,NS
      N=NP(I,J)
      IF (KODE(N) .NE. -2) GO TO 280
      IF (P(N) .LT. SEEPVAL(I)) ICHECK=1
      IF (ICHECK .GT. 0) GO TO 300
      KODE(N)=2
      P(N)=SEEPVAL(I)
      GO TO 300
  280 IF (KODE(N) .NE. 2) GO TO 300
      IF (ICHECK .GT. 0) GO TO 290
      IF ((Q(N) .LT. 0. .AND. NSPKDS(I,2) .LT. 0) .OR.
     1    (Q(N) .GT. 0. .AND. NSPKDS(I,2) .GT. 0)) GO TO 300
  290 KODE(N)=-2
      Q(N)=0.
      ICHECK=1
  300 CONTINUE
  310 CONTINUE
C----------------------------------- MODIFY POTENTIAL SURFACE FLUX BOUNDARIES
  320 IF (NUMEP .LT. 1) GO TO 410
      DO 380 I=1,NUMEP
      N=NSIC(I)
      IF (KODE(N) .NE. 4) GO TO 330
      IF ((EI .LT. 0. .AND. Q(N) .GE. EI*WIDBTP(I,1)) .OR.
     &    (EI .GE. 0. .AND. Q(N) .LE. EI*WIDBTP(I,1))) GO TO 380
      KODE(N)=-4
      Q(N)=EI*WIDBTP(I,1)
      GO TO 380
  330 IF (KODE(N) .NE. -4) GO TO 380
      IF (EI .GE. 0.) GO TO 350
      IF (P(N) .GT. PL) GO TO 340
      KODE(N)=4
      P(N)=PL
      GO TO 380
  340 IF (DABS(P(N)) .GT.-.001*PL) WIDBTP(I,2)=WIDBTP(I,2)*DABS(PL/P(N))
      IF (DABS(P(N)) .LE. -.001*PL) WIDBTP(I,2)=1.
      GO TO 370
  350 IF (P(N) .LT. 0.) GO TO 360
      KODE(N)=4
      P(N)=0.
      GO TO 380
  360 IF (EI .GE. 0.)GO TO 380
      IF (DABS(P(N)-PL) .GT. -.001*PL) 
     &               WIDBTP(I,2)=WIDBTP(I,2)*DABS(PL/(P(N)-PL))
      IF (DABS(P(N)-PL) .LE. -.001*PL) WIDBTP(I,2)=1.
      IF(P(N) .LT.PL) WIDBTP(I,2)=1.
  370 IF (WIDBTP(I,2) .GT. 1.) WIDBTP(I,2)=1.
      Q(N)=WIDBTP(I,2)*EI*WIDBTP(I,1)
  380 CONTINUE
C-------------------------------------------- TEST FOR CONVERGENCE
  410 IF (ITER .LE. 2) GO TO 450
C============================================= READ FROM FILE 3 
      IF( INCOR .EQ. 0 ) THEN
       REWIND 3
       READ (ND3) C
      ELSE
C==========================================================================
      DO 804 I=1,INUMNP
       C(I)=CPI(I,4)
  804 CONTINUE
      END IF
C=========================================================================-
      EPSLON=0.
      DO 420 I=1,NUMNP
      AA = DABS(P(I)-C(I))
      IF (AA .LT. EPSLON) GO TO 420
      EPSLON=AA
      NMAX=I
  420 CONTINUE
C      WRITE (ND6, 240)  HED,TIME,LA
       WRITE (ND6, 430) ITER,EPSLON,NMAX
        WRITE (NDCN, 431) ITER,EPSLON,NMAX
C
       WRITE(110,434) EPSLON,NMAX
C
      IERREND=0       
      IF(EPSLON .LE. TOL) GO TO 460
      IERREND=1
      IF( LA .LE. 0) THEN
c        IF ( ITER .GT. MAXITO  .OR. EPSLN2 .LE. EPSLON  ) THEN      
        IF ( ITER .GT. MAXITO   ) THEN
          WRITE (ND6,559)                                 
          WRITE (NDCN,559)            
  559     FORMAT(//' ********* ERROR IN INITIAL CONDITION *********'//)
C
       WRITE(110,438)  !ERR
C
          RETURN 
        ELSE 
           GO TO 450
        ENDIF
      ELSE
c       IF(ITER .GT. MAXIT .OR. EPSLN2 .LE. EPSLON) GO TO 455      
       IF(ITER .GT. MAXIT ) GO TO 455
      ENDIF
C======================================================= WRITE TO FILE 3 ==
  450 IF( INCOR .EQ. 0 ) THEN
       REWIND ND3
       WRITE (ND3) P
      ELSE
C==========================================================================
 5555  DO 805 I=1,INUMNP
       CPI(I,4)=P(I)
  805  CONTINUE
      END IF
C=========================================================================
      EPSLN2=EPSLN1
      EPSLN1=EPSLON
      GO TO 250
C------------------------------------------------- TIME MODIFY FOR RECAL
  455 TIME=TIME-DT
      DT=.5*DT
      IFLGTM = 0
C     IF(DT .LT.DT1) DT=DT1
      TIME=TIME+DT
      ITER=0
      IF(LAA .GE. NCUL) GO TO 470
      LAA=LAA+1
      GO TO 230
C----------------------------------------------------------- PRINT RESULTS
  460 IF(LA.LE.0) GO TO 520
      DO 465 N=1,NUMNP
      QTOT =QTOT+Q(N)*DT
  465 CONTINUE
C
C
      IF( LA .EQ. 0 ) GO TO 520
      IF( TIME .GE. TMAX ) GO TO 470
      IF( IFLGTM .EQ. 1 ) GO TO 470
  468 IF(DABS(TIME-TMAX) .LE. .001*DT) GO TO 470
C------------------------------------------------------- PRINT OUT CONTROL  
C      K=MOD(LA,NPRINT)
C      IF(K .NE. 0) GO TO 539
  470  CONTINUE 
c       WRITE (ND6, 240)  HED,TIME,LA,AL
       WRITE (ND6, 430) ITER,EPSLON,NMAX
C
       WRITE(110,434) EPSLON,NMAX
 434  FORMAT("I,",E14.6,",",I5)
C
      WRITE (ND6, 480) QTOT
  480 FORMAT ( /35H CUMULATIVE INFLOW INTO SYSTEM IS =1PE13.5)
      IF (NBW .GT. 0) WRITE (ND6, 490) QP,AL
  490 FORMAT (/22H DISCHARGE FROM WELL =1PE13.5/
     1           22H WATER LEVEL IN WELL =1PE13.5/)
      IF (NUMEP .GT. 0) WRITE (ND6, 500) EI,PL
  500 FORMAT (/25H POTENTIAL SURFACE FLUX =1PE13.5/
     1           32H MINIMUM SURFACE PRESSURE HEAD =1PE13.5/)
  520 CONTINUE
  430 FORMAT (/49H MAXIMUM CHANGE IN PRESSURE HEAD DURING ITERATION,I3,
     1 4H WAS,1PE13.5,8H AT NODE,I5)
  431 FORMAT (/46H MAX. CHANGE IN PRESSURE HEAD DURING ITERATION,I3,
     1 4H WAS,1PE13.4,8H AT NODE,I5)
  432 FORMAT (/' CALCULATING ....ITERATION No.=',I3)
      IF (NBW .GT. 0) WRITE (ND6, 440) AL,QA
  440 FORMAT (/' WATER LEVEL IN WELL WAS,',1PE13.5,
     &         ' DISCHARGE FROM AQUIFER WAS,',E13.5)
C------------------------------------------ PRINT OUT HEAD & DISCHARGE    
	IF( IFLGTM .EQ.0 ) GO TO 539
       IFLGTM = 0
C
      WRITE(NDCN,*) '//////Printing//////'
      CALL PRINTO (Q,P,XY(1,2),KPM,Z,NUMNP,KAT,QMM,QPP,NUMNPZ)
C-------------------------------------------------- VELOCITY & CONTER 
      CALL ELFLOW(XY(1,1),XY(1,2),KX,P,XYK(1,1,1),XYK(1,1,2),XYP(1,1,1)
     &  ,XYP(1,1,2),VXY(1,1),VXY(1,2),H,Z,NUMNP,NUMEL,
     &            NUMMAT,MK,MP,TIME,LA,IVFLAG,NUMNPZ)
      WRITE(ND10) LA,TIME
      WRITE(ND10) (KODE(I),P(I),Q(I),I=1,NUMNP)
      WRITE(ND10) (VXY(I,1),VXY(I,2),I=1,NUMEL)
C
C------------------------------------ NEW B.C. FOR NON-STEADY PROBLEM
  539 IF( NTYPE .EQ.  1 ) GO TO 560
      IF( NTYPE .EQ. -2 .AND.  LA .EQ. 0 ) GO TO 560
      IF( NTYPE .LE. -1 ) GO TO 536
      IF( LA .GE. 1     ) GO TO 536
      CALL NEWN (SEEPVAL,KODE,P1,P,Q,XY(1,2),NSPKDS(1,1),NP,
     &            NUMNP,INIT,NDIM,MAXSP,NSEEP)
      DO I=1,NUMNP
        P1(I)=P(I)
      END DO 
      NTYPE=-1
      GO TO 220
C---------------------------------------------- SET UP NEW TIME STEP  
  536 LA = LA + 1
      DT1=DT
      AL1=AL
      IERREND=0
      IF (MAXIT.LT.1 .OR. KCON .GE.1 )GO TO 540
      IF (EPSLON .LE. TOL) GO TO 540
      IERREND=1             !!!!!!!!!!!!!!!!!!!!!!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      GO TO 560
 540  IF (DABS(TIME-TMAX) .LE. .001*DT) GO TO 560
C-----------------------------------------------------------------
      IF (DT .LT. DTMAX ) DT=DMUL*DT
C-----------------------------------------------------------------
      IF (DT .GT. DTMAX) DT=DTMAX
      IF (NT .LE.   0) GO TO 550
C
      NN=0
      IF(TIME .LT. TM(1) ) GO TO 551
      DO NN = 1,NT-1
        IF( TM(NN).LE.TIME .AND. TM(NN+1).GT.TIME)  GO TO 551
      END DO
      WRITE(110,438)
      STOP 9999
  551 CONTINUE
C **********************
      DT2=TM(NN+1)-TIME
	IF( (DT2-DT).GT.DT*0.2) GO TO 550
	DT=DT2
	TIME1= TM(NN+1)
C *******************
      IF(TIME1.GT. TMAX) GO TO 550
      TIME = TIME1
      IFLGTM =1
      GO TO 553
  550 IF ((TIME+DT) .GT. TMAX .OR. DABS(TMAX-TIME-DT) .LT. 0.2*DT) 
     &    DT=TMAX-TIME
      TIME=TIME+DT
  553 ITER=0
      LAA=0
      GO TO 230
  560 CONTINUE
C--------------------- FOR NTYPE=-2 ----- FINAL STEADY STATE ----------
      IF ( NTYPE .NE. -2 ) GO TO 570
      IF ( LA .EQ. 0 ) GO TO 570 
      LA=0
      TIME=0.0
      WRITE (ND6,2400)
      WRITE (NDCN,2400)
 2400 FORMAT (/,' ***** FINAL STAEDY STATE *****'/)
      GO TO 221
C-------------------------------- NEW BOUNDARY CONDITION SET & RESTART
  570 READ(ND5, 10) NCHANG
      IF(NCHANG.LE.0) GO TO 580
      TI=TIME
      CALL NEXTCO (KODE,P,P1,Q,KX,LA,DT,QTOT,ITER,TI,LAA,NUMNP,
     1NUMEL,DTMAX,DMUL,TMAX,TOL)

      GO TO 220
  580 CONTINUE 
C
       IF(IERREND.EQ. 0) THEN
       WRITE(110,439)  
439    FORMAT("NOMARL, 0.0  ,0  ")
       ELSE
       WRITE(110,438)  
438    FORMAT("ERROR, 0.0  ,0  ")
       END IF
C
C
      RETURN
      END
      SUBROUTINE MATIN (XK,YK,XP,YP,XC,YC,MK,MP,NUMMAT,KAT,KTYPE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                                      
      DIMENSION XK(MK,NUMMAT),YK(MK,NUMMAT),XP(MP,NUMMAT),YP(MP,NUMMAT),    
     &          XC(MP,NUMMAT),YC(MP,NUMMAT)                                           
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C------------------------------------------------------------ CONTROL UNINTS
      CONS1=1.0D0
      CONS2=1.0D0
      IF(KAT.EQ.1.OR.KAT.EQ.3) CONS2=0.159154943D0
      WRITE(ND6, 20) CONS1,CONS2                                              
   20 FORMAT (///50H  UNITS ARE CONTROLED BY THE FOLLOWING CONSTANTS -//    
     1 9H  CONS1 =1PE12.3,3X,7HCONS2 =E12.3////'***** MATERIAL INFORMATI
     2ON ******'//9H MAT. No.,7X,3HK-X,9X,       
     3 3HK-Z,10X,8HPOROSITY,10X,2HSs,10X,3HTHR,9X,3HFAN/)
C---------------------------------------------------------------- PROPERTY
      DO 51 M=1,NUMMAT
      READ(ND5, 25) C1(M),C2(M),POR(M),SS(M),THR(M),FAN(M)
   25 FORMAT(7F10.0)                              
      WRITE(ND6, 40) M,C1(M),C2(M),POR(M),SS(M),THR(M),FAN(M) 
   40 FORMAT (I6,4X,1P2E12.3,7X,0PF6.3,1PE17.3,0PF10.3,0PF12.2) 
   30 FORMAT(6F10.0)  
   51 CONTINUE                                                   
      IF ( KTYPE.EQ. -1 ) RETURN
C--------------------------------------------------------------- THETA , Kr 
      DO 50 M=1,NUMMAT                    
      NUMKM=NUMK(M)                                                         
      NUMPM=NUMP(M)                                                         
      IF (KTYPE .GE. 1) GO TO 1  
      READ(ND5, 30) (XK(I,M),YK(I,M),I=1,NUMKM) 
      GO TO 11
    1 KANN=NUMKM-2
      ANN=DBLE(KANN)
      DTHO=POR(M)-THR(M)
      DTH=DTHO/ANN
      XK(1,M)=0.0D00
      YK(1,M)=0.0D00
      DO 6 I=2,NUMKM
      I2=I-2
      XII=DBLE(I2)
      TH=THR(M)+DTH*XII
      XK(I,M)=TH
      SE=(TH-THR(M))/DTHO
      IF (KTYPE .NE. 1) GO TO 2
      YK(I,M)=SE**FAN(M)
      GO TO 6
    2 F=FAN(M)*(SE-1.0)
      YK(I,M)=DEXP(F)
    6 CONTINUE
C----------------------------------------------------------- THETA - PHAI 
   11 READ(ND5, 30) (XP(I,M),YP(I,M),I=1,NUMPM)                               
      DO 60 I=1,NUMPM                                                       
      YP(I,M)=-YP(I,M)                                                      
   60 CONTINUE                                                              
C      NUMC(M)=NUMP(M)-1                                                     
      NUMC(M)=NUMP(M)                                                     
      NUMCM=NUMC(M)                                                         
C      DO 70 I=1,NUMCM                                                       
      DO 70 I=1,NUMCM-1                                                       
      AA     =.5D0*(XP(I,M)+XP(I+1,M))                                         
      BB     =(XP(I+1,M)-XP(I,M))/(YP(I+1,M)-YP(I,M))                       
      XC(I,M)=AA                                                            
      YC(I,M)=BB                                                            
   70 CONTINUE
      XC(NUMCM,M)=XP(NUMCM,M)
      YC(NUMCM,M)=0.0
   50 CONTINUE 
C ------------------------------------------------ INFORMATION PRINT OUT
      DO 120 M=1,NUMMAT  
      NUMKM=NUMK(M)                                                         
      NUMPM=NUMP(M)                                                         
      NUMCM=NUMC(M)                                                         
      WRITE(ND6, 80) M                                                        
   80 FORMAT (//68H TABLE OF RELATIVE CONDUCTIVITY VERSUS MOISTURE CONTE    
     1NT FOR MATERIAL I3/) 
      IF (KTYPE .EQ. 0 ) GO TO 85 
      IF (KTYPE .EQ. 1 ) THEN 
        WRITE (ND6,81)
      ELSE
        WRITE (ND6,82)
      ENDIF
   81 FORMAT (//,' ****** Kr Are Caluculated By Irmay Type  ****** '/)
   82 FORMAT (//,' ****** Kr Are Caluculated By Hillel Type  ****** '/)
   85 WRITE(ND6, 90) (YK(I,M),XK(I,M),I=1,NUMKM)                              
   90 FORMAT (5(1P2E11.3,2X))
      WRITE(ND6, 100) M                                                       
  100 FORMAT (//60H TABLE OF PRESSURE HEAD VERSUS MOISTURE CONTENT FOR M    
     1ATERIAL I3/)                                                          
      WRITE(ND6, 90) (YP(I,M),XP(I,M),I=1,NUMPM)                              
      WRITE(ND6, 110) M                                                       
  110 FORMAT (//73H TABLE OF SPECIFIC MOISTURE CAPACITY VERSUS MOISTURE     
     1CONTENT FOR MATERIAL I3/)                                             
      WRITE(ND6, 90) (YC(I,M),XC(I,M),I=1,NUMCM)                              
  120 CONTINUE                                                              
      RETURN                                                                
      END  
C
      SUBROUTINE NPIN (KODE,Q,P,P1,X,Y,WIDTH,NSIC,KX,B,Z,NUMDP,NUMNP,
     & INIT,NERR,NUMEP,KAT,NUMEL,IRNDE,NUMNPZ)                                                 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
	  DIMENSION KODE(NUMNP),Q(NUMNP),P1(NUMNP),X(NUMNP),Y(NUMNP),           
     & P(NUMNP),WIDTH(NUMDP),NSIC(NUMDP),NN(6),TT(6),KX(NUMEL,5),
     & B(NUMNP),Z(NUMNPZ)                  
      NERR=0
	  Z(1)=0.0                                                               
      WRITE(ND6, 10)                                                          
   10 FORMAT (////24H NODAL POINT INFORMATION //10H  NODE NO.,6X,4HKODE,    
     & 7X,3HX,R,12X,3HY,Z,11X,5H.PSI.,12X,1HQ/)                             
      NPR=0                                                                 
      L=0                                                                   
   20 L=L+1                                                                 
c      IF( IRNDE.EQ. 0)   READ(ND5, 40) N,KODE(N),X(N),Y(N),P1(N),Q(N),ZZ                          
c      IF( IRNDE.EQ. 1)   READ(ND5, *) N,KODE(N),X(N),Y(N),P1(N),Q(N),ZZ
      IF( IRNDE.EQ. 0)   READ(ND5, 40) N,KODE(N),X(N),Y(N),P1(N),Q(N)                         
      IF( IRNDE.EQ. 1)   READ(ND5, *) N,KODE(N),X(N),Y(N),P1(N),Q(N)
	  ZZ=0.0D00
C
C      WRITE(NDCN,*) N,KODE(N),X(N),Y(N),P1(N),Q(N),ZZ                            
   40 FORMAT (2I5,5F10.0)
      PINT=Y(N)
      IF( INIT.EQ. 0 ) PINT=0.0D00
C      IF( KAT.LE.0) THEN
C       Z(N)=ZZ
C       PINT=0.0D00
C       IF( INIT.EQ.0) PINT=-Z(N)
C      END IF                                             
      P1(N)=P1(N)-PINT                                                 
C
      IF (N-L) 50,90,70                                                     
   50 WRITE(NDCN, 60) N                                                        
   60 FORMAT (20H ERROR IN NPIN AT N= I5)                                   
      NERR=1                                                                
      RETURN                                                                
   70 DENO=N-L+1                                                            
      DX=(X(N)-X(NPR))/DENO                                                 
      DY=(Y(N)-Y(NPR))/DENO                                                 
      DP=(P1(N)-P1(NPR))/DENO
C      IF( KAT.LE.0 ) DZ=(Z(N)-Z(NPR))/DENO                                 
   80 X(L)=X(L-1)+DX                                                        
      Y(L)=Y(L-1)+DY                                                        
      P1(L)=P1(L-1)+DP
C      IF( KAT.LE.0 ) Z(L)=Z(L-1)+DZ                                                       
      KODE(L)=KODE(L-1)                                                     
      Q(L)=Q(L-1)                                                           
      L=L+1                                                                 
      IF (L .LT. N) GO TO 80                                                
   90 NPR=N                                                                 
      IF (L .LT. NUMNP) GO TO 20                                            
      I=0
C------------------------------
      CALL KODCHK(KX,KODE,B,NUMEL,NUMNP)
C--------------------------------
      DO 110 N=1,NUMNP
C      ZZ=0.0D00
C      IF( KAT.LE.0) ZZ=Z(N)                                                      
C      WRITE(ND6, 100) N,KODE(N),X(N),Y(N),P1(N),Q(N),ZZ                          
      WRITE(ND6, 100) N,KODE(N),X(N),Y(N),P1(N),Q(N)                         
  100 FORMAT (2I10,5E15.6)                                                  
      P(N)=P1(N)                                                            
  110 CONTINUE                                                              
      IF( NUMEP .LT. 1) GO TO 140                                           
      DO 120 I=1,NUMEP                                                      
  120 WIDTH(I)=0.                                                           
      NUMEP1=NUMEP - 1                                                      
      DO 130 I=1,NUMEP1                                                     
      K=I+1                                                                 
      J=NSIC(I)                                                             
      J1=NSIC(K)                                                            
      RI=1.0                                                                
      RK=1.0                                                                
      IF( KAT .EQ. 1) THEN                                           
        RI=0.25D0*(X(J1)*3.0D0+X(J))   
        RK=0.25D0*(X(J)*3.0D0+X(J1))                         
        RI=RI*2.0D0*3.141592654D0                               
        RK=RK*2.0D0*3.141592654D0
	END IF           
  125 WID=(X(J1)-X(J))/2.D0                                              
      WIDTH(I)=WIDTH(I)+WID*RI                                              
      WIDTH(K)=WIDTH(K)+WID*RK                                              
  130 CONTINUE                                                              
      WRITE(ND6, 162)                                                         
  162 FORMAT(/' SOIL SURFACE WIDTH'/6('   NODE   WIDTH     '))              
      J=0                                                                   
      N=0                                                                   
  150 DO 160 K=1,6                                                          
      NN(K)=0                                                               
      TT(K)=0.D0                                                              
  160 CONTINUE                                                              
  170 N=N+1                                                                 
      IF( N .GT. NUMEP ) GO TO 180                                          
      J=J+1                                                                 
      NN(J)=NSIC(N)                                                         
      TT(J)=WIDTH(N)                                                        
      IF( J .LT. 6) GO TO 170                                               
  180 WRITE(ND6, 190) (NN(K),TT(K),K=1,6)                                     
  190 FORMAT(6(I7,E13.5))                                                   
      J=0                                                                   
      IF (N .LT. NUMEP) GO TO 150                                           
  140 RETURN                                                                
      END                                                                   
         SUBROUTINE MAFILL (KX,SANG,CANG,NUMEL,MBAND,IRELM)                          
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
	  COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN                              
      DIMENSION KX(NUMEL,5),SANG(NUMEL),CANG(NUMEL)                         
CC ----                                                                     
      NUM=0                                                                 
      WRITE(ND6, 10)                                                          
   10 FORMAT (////20H ELEMENT INFORMATION//55H ELEMENT   C O R N E R  N     
     1 O D E S    MATERIAL    ANGLE/)                                       
      DO 60 N=1,NUMEL
         NUM1=NUM                                                       
      IF (NUM-N) 20,60,40                                                   
   20 IF( IRELM.EQ.0 ) READ(ND5, 30) NUM,(KX(NUM,I),I=1,5),SANG(NUM)
      IF( IRELM.EQ.1 ) READ(ND5, *) NUM,(KX(NUM,I),I=1,5),SANG(NUM)
   30 FORMAT (6I5,F10.0,I5) 
      KG=1                                                   
      IF( KX(NUM,1).GT.0) GO TO 300
       KX(NUM,1)=KX(NUM1,1)+NUM-NUM1
       KX(NUM,2)=KX(NUM1,2)+NUM-NUM1
       KX(NUM,3)=KX(NUM1,3)+NUM-NUM1
       KX(NUM,4)=KX(NUM1,4)+NUM-NUM1
  300 IF (KX(NUM,4) .EQ. 0) KX(NUM,4)=KX(NUM,3)                             
      IF (NUM .EQ. N) GOTO 60                                              
   40 DO 50 I=1,4                                                           
      KX(N,I)=KX(N-1,I)+KG                                                   
   50 CONTINUE                                                              
      KX(N,5)=KX(N-1,5)                                                     
      SANG(N)=SANG(N-1)                                                     
   60 CONTINUE                                                              
      AA=3.141592654D0/180.D0                                    
      DO 90 N=1,NUMEL                                                       
      BB=AA*SANG(N)                                                         
      WRITE(ND6, 70) N,KX(N,1),KX(N,2),KX(N,3),KX(N,4),KX(N,5),SANG(N)        
   70 FORMAT (I6,I9,3I6,I10,E14.3)                                          
      CC=DSIN(BB)                                                           
      SANG(N)=CC                                                            
      CC=DCOS(BB)                                                           
      CANG(N)=CC                                                       
   90 CONTINUE                                                          
C --------------------------------------------- CALCULATE BAND-WIDTH
      MBAND=0                                                               
      DO 101 N=1,NUMEL                                                      
      IF( KX(N,5) .LE. 0 ) GO TO 101
      LS=10000                                                              
      DO 112 I=1,4                                                          
      IF(KX(N,I)-LS) 111,112,112                                            
  111 LS=KX(N,I)                                                            
  112 CONTINUE                                                              
      DO 102 I=1,4                                                          
      II=KX(N,I)                                                            
      ME=II-LS                                                              
      IF(ME .GT. MBAND) THEN
       MBAND=ME
	 IXMAX = N                                                 
      ENDIF 
  102 CONTINUE                                                              
  101 CONTINUE                                                              
C    
      MBAND=MBAND+1                                                           
      WRITE (NDCN,200) MBAND,IXMAX                                                    
      WRITE (ND6,200) MBAND,IXMAX                                                    
  200 FORMAT (///,' ***** BAND WIDTH = ',I5,'  ELEM =',I5)                             
      RETURN                                                                
      END                                                                   
      SUBROUTINE KODCHK(KX,KODE,CC,NUMEL,NUMNP)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION KX(NUMEL,5),KODE(NUMNP),CC(NUMNP)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN

C
      DO 100 I=1,NUMNP
  100   CC(I)=0.0
C
      DO 10 I=1,NUMEL
       IF( KX(I,5).GE.1) GO TO 10
       NUS1=4
        IF( KX(I,3).EQ.KX(I,4)) NUS1=3
       DO 11 K=1,NUS1
        KK=KX(I,K)
C         WRITE(NDCN,*)' NODE = ',KK
        IF( KODE(KK).GT. 0 .OR. CC(KK).GT.0.0 ) GO TO 11
      DO 20 J=1,NUMEL
       IF( KX(J,5).LE.0) GO TO 20
      NUS2=4
       IF(KX(J,3).EQ. KX(J,4)) NUS2=3
        DO 21 L=1,NUS2
        LL=KX(J,L)
         IF( LL.EQ.KK ) GO TO 12
   21 CONTINUE 
   20 CONTINUE
       KODE(KK)=1
   12  CC(KK)=1.0
   11 CONTINUE
   10 CONTINUE
C       DO 301 I=1,NUMNP
C       IF( KODE(I).GT. 0 ) WRITE(ND6,300)I,KODE(I)
C  301 CONTINUE
C  300 FORMAT(5I10)
      RETURN
      END 
      SUBROUTINE COLHT ( KX,IS,RMK,NUMNP,NUMNP1,NUMEL,NWK ) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
      DIMENSION KX(NUMEL,5),IS(NUMNP1),ISQ(200),LM(4),IST(3),
     &          RMK(2,NUMEL)
C
      DO 3000 I=1,NUMEL
      RMK(1,I)=KX(I,1)
      RMK(2,I)=KX(I,1)
       DO 3010 J=2,4
       IF( RMK(1,I).GT.KX(I,J)) RMK(1,I)=KX(I,J)
       IF( RMK(2,I).LT.KX(I,J)) RMK(2,I)=KX(I,J)
 3010 CONTINUE
 3000 CONTINUE
      IS(1) = 1
      DO 200 I=1,NUMNP
      ISQ(1) = I
      K = 1
      J = 1
   10 IF( J .GT. NUMEL ) GO TO 100
       IF( KX(J,5).LE.0) GO TO 35
C
       IF(RMK(1,J).GT.I .OR. RMK(2,J).LT.I ) GO TO 35
C
      DO 20 M=1,4
   20 LM(M) = KX(J,M)
      NUS = 4
      IF( LM(4) .EQ. LM(3) ) NUS = 3
       KK=NUS-2
      DO 30 KI=1,KK
      IST(1) = LM(1)
      IST(2) = LM(KI+1)
      IST(3) = LM(KI+2)
      DO 40 II=1,3
      IF( IST(II) .EQ. I ) GO TO 50
   40 CONTINUE
      GO TO 30
   50 DO 60 II=1,3
      JJ = IST(II)
      DO 70 L=1,K
      IF( ISQ(L) .EQ. JJ ) GO TO 60
   70 CONTINUE
      K = K + 1
       IF( L.GT.200 ) STOP 0000
      ISQ(K) = JJ
   60 CONTINUE
   30 CONTINUE
   35   J = J + 1
      GO TO 10
  100 IS(I+1) = IS(I) + K
  200 CONTINUE
      NUM = NUMNP1
      NWK = IS(NUM) - 1
      RETURN
      END
C
      SUBROUTINE ADRSK ( KX,IS,ISP,RMK,NUMNP,NUMNP1,NUMEL,NWK )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN	
      DIMENSION KX(NUMEL,5),IS(NUMNP1),ISP(NWK),LM(4),IST(3),
     &          RMK(2,NUMEL)
C ------------------------------------------- SEARCH OF RELATED NODE NUMBER
      N = 0
      DO 100 I=1,NUMNP
      N = N + 1
      ISP(N) = I
      K = 1
      J = 1
      I1 = IS(I)
      ISS = IS(I+1) - I1
   10 IF( J .GT. NUMEL .OR. K .EQ. ISS ) GO TO 100
       IF( KX(J,5).LE.0 ) GO TO 35
C
       IF(RMK(1,J).GT.I .OR. RMK(2,J).LT.I ) GO TO 35
C
      DO 20 M=1,4
   20 LM(M) = KX(J,M)
      NUS = 4
      IF( LM(4) .EQ. LM(3) ) NUS = 3
        KK=NUS-2
      DO 30 KI=1,KK
      IST(1) = LM(1)
      IST(2) = LM(KI+1)
      IST(3) = LM(KI+2)
      DO 40 II=1,3
      IF( IST(II) .EQ. I ) GO TO 50
   40 CONTINUE
      GO TO 30
   50 DO 60 II=1,3
      JJ = IST(II)
      I2 = I1 + K - 1
      DO 70 L=I1,I2
      IF( ISP(L) .EQ. JJ ) GO TO 60
   70 CONTINUE
      K = K + 1
      N = N + 1
      ISP(N) = JJ
   60 CONTINUE
   30 CONTINUE
   35   J = J + 1
      GO TO 10
  100 CONTINUE
      RETURN
      END
      SUBROUTINE RESET(KODE,A,B,Q,P,P1,T,D,CR,C,X,Y,KX,SANG,CANG,XK,YK,    
     &XP,YP,XC,YC,AW,CPI,IC,IS,ISP,Z,
     & NUMNP,MBAND,NWK,NUMEL,NUMMAT,MK,MP,INUMNP,
     &INWK,NUMNP1,NWK1,NUMNPZ)
C                                                        
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*4  HED
      INTEGER*4 ILENR,ILENI
C
      DIMENSION A(NWK),B(NUMNP),Q(NUMNP),P(NUMNP),P1(NUMNP),        
     &T(NUMNP),D(NUMNP),CR(NUMNP),C(NUMNP),X(NUMNP),Y(NUMNP),               
     &KX(NUMEL,5),SANG(NUMEL),CANG(NUMEL),XK(MK,NUMMAT),YK(MK,NUMMAT),
     &XP(MP,NUMMAT),YP(MP,NUMMAT),XC(MP,NUMMAT),YC(MP,NUMMAT),
     &KODE(NUMNP),IC(NUMEL),
     &AW(INWK),CPI(INUMNP,5),IS(NUMNP1),ISP(NWK1),Z(NUMNPZ)
      DIMENSION XY(8),PNOD(4),LM(4),S(4,4)
C                      
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /FIX/ HED(20),INIT,KTYPE,MAXIT,KAT,INTEG,NSEEP,NCUL,NPRINT,
     &MAXITO,NTYPE,JMES
      COMMON /MESHZU/ MESH,INCOR,IVFLAG,NT,NUMEP,ITQ,NBOUN,NBNQ
C      COMMON /CRTL/LA,DT,DT1,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL1/LA,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL2/DT,DT1
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN

C
C      CHARACTER*1 KKK
C      CHARACTER*8 AAA,BBB
      IC(1)=0
	  XY(1)=0.0
	  PNOD(1)=0.0
       NINT=2
C
      IF (LA .GT. 1) DTT=.5*DT/DT1                                          
C---- INITIALIZATION AND EXTRAPOLATION OF P VALUES
C      WRITE(NDCN,*)" ENTER DO 40 IN RESET"                          
      DO 40 I=1,NUMNP                                                       
       B(I)=0.0D00                                                           
       D(I)=0.0D00
       A(I)=0.0D00                                                           
        IF(ITER.GT.0) GO TO 10                                                
         IF (LA .LE. 1) PP =P1(I)                                              
         IF (LA .GT. 1) PP =P(I)+DTT*(P(I)-P1(I))                              
          P1(I)=P(I)                                                          
         IF (KODE(I) .LT. 1) P(I)=PP                                           
       GO TO 40                                                              
   10   PP=.5*(P(I)+P1(I))                                                    
        IF (KODE(I) .LT. 1) P(I)=PP                                           
   40 CONTINUE
      DO 30 I=NUMNP+1,NWK                                                      
      A(I)=0.0D00                                                         
   30 CONTINUE
C                                                               
C      WRITE(NDCN,*)" ENTER DO 70 IN RESET"                          
      DO 70 M=1,NUMMAT                                                  
      CMUL=CONS1*SS(M)/POR(M)                                          
      DO 50 I=1,NUMNP                                                       
       CR(I)=-1.                                                             
   50 CONTINUE                                                              
      DO 60 N=1,NUMEL                                                       
       IF (KX(N,5).EQ. M) THEN                                          
        IWET=0
        LM(1)=KX(N,1)
        LM(2)=KX(N,2)
        LM(3)=KX(N,3)
        LM(4)=KX(N,4)
        SANGN=SANG(N)
        CANGN=CANG(N)
        MMAT = M
        PORM = POR(M)
        C1M = C1(M)
        C2M = C2(M)
        SSM = SS(M)
        NMK=NUMK(M)                                                           
        NMP=NUMP(M)                                                           
        NMC=NUMC(M)                                                           
C
C        IF( JMES.EQ.0) THEN
          CALL ELEM (CMUL,KODE,S,B,P,P1,T,D,CR,C,X,Y,Z,
     &   SANGN,CANGN,XK,YK,XP,YP,XC,YC,LM,PORM,C1M,C2M,
     &   NMK,NMP,NMC,NUMNP,NUMMAT,MK,MP,MMAT,KAT,NUMNPZ)

c        CALL ELEMS (CMUL,S,B,P,D,X,Y,Z,
c     &   SANGN,CANGN,XK,YK,XP,YP,XC,YC,LM,C1M,C2M,
c     &   NMK,NMP,NMC,NUMNP,NUMMAT,MK,MP,MMAT,KAT,NUMNPZ)
C
C       ELSE
C       I2 = 0
C       IF( JMES.EQ.1.OR.(JMES.EQ.2.AND.IC(N).LE.0 ))THEN
C        DO 65 I=1,4
C        II = LM(I)
C        I2 = I2 + 2
C        ZZ=0.0D00
C         IF( KAT.LE.0 ) ZZ=Z(II)
C        XY(I2-1) = X(II)
C        XY(I2)   = Y(II)
C  65    PNOD(I)  = P(II)-ZZ
C         CALL ELEMI (CMUL,S,B,PNOD,D,XY,KX,XK,YK,XP,YP,XC,YC,
C     &             NUMNP,NUMEL,NUMMAT,MK,MP,M,KAT,N,LA,NINT) 
C
C       ELSE
C      DO 66 I=1,4
C         IF( IC(N).LE.0 ) WRITE(NDCN,*)'--- STOP  IC=0 ',N
C       II=IC(N)+I-1
C       IF( II.GT.4 ) II=II-4
C      III = LM(I)
C      I2 = I2 + 2
C      ZZ=0.0D0
C     IF( KAT .LE. 0 ) ZZ=Z(III)
C     XY(I2-1) = X(III)
C      XY(I2)   = Y(III)
C   66 PNOD(I)  = P(II)-ZZ
C        CALL INFMAT (CMUL,S,B,PNOD,D,XY,LM,
C     &  XK,YK,XP,YP,XC,YC,
C     &  NUMNP,NUMMAT,MK,MP,M,KAT,N,LA)                    
C      END IF   
C      END IF
C
       IF( NSOLV .EQ. 0 ) THEN        
c        WRITE(NDCN,*)"ENTER ASEMBAG IN RESET"
         CALL ASEMBAG(A,S,LM,MBAND,NUMNP)
       ELSE
c        WRITE(NDCN,*)"ENTER ASEMBAC IN RESET"
         CALL ASEMBAC(A,S,LM,IS,ISP,NWK,NUMNP1)     
       END IF
      END IF     
   60 CONTINUE                                                              
   70 CONTINUE                                                                 
c     WRITE(NDCN,*)" OUT DO 70 IN RESET"                          
C********************************************************* WRITE TO FILE 2
      IF( INCOR .EQ. 0 ) THEN
      REWIND ND2
       REC = NWK/ILENR
       IREC = IDNINT(REC)
       IRAM = MOD(NWK,ILENR)
      IF( IREC .LT. 1 ) GO TO 123
      DO I=1,IREC    
      WRITE(ND2) (A((I-1)*ILENR+J),J=1,ILENR)
      END DO
  123 IF(IRAM .GT. 0 ) WRITE(ND2) (A(ILENR*IREC+J),J=1,IRAM)
      WRITE(ND2) B                
C************************************************************************
      ELSE
 2222 DO 200 I=1,NUMNP
        BI = B(I)
  200    CPI(I,5) = BI
      DO 250 I=1,NWK
  250    AW(I)=A(I)
       END IF
C======================================================================
C---- COMPLETE CONSTRUCTION OF RHS VECTOR AND FORM EFFECTIVE MATRIX         
      PP=2.0D00
      PD=1.0D00                                                                 
      IF( LA .LE. 1 .OR. INTEG .EQ. 1) PP=1.0D00
      IF( LA .LE. 0 ) PD=0.0D00                               
      DO 110 I=1,NUMNP                                                      
      B(I)=PP*(CONS2*Q(I)-B(I)+D(I)*P1(I)*PD/DT)                               
      IF (LA.LE.1.OR.INTEG.EQ.1) GOTO 100
      IF( NSOLV .EQ. 0 ) THEN                                  
       DO 90 J=1,MBAND                                                       
        K=I+J-1                                                               
        IF (K .GT. NUMNP) GOTO 80  
         JI=MBAND*(I-1)+J                                          
         B(I)=B(I)-A(JI)*P1(K)                                                
   80     K=I-J+1                                                               
         IF (J .EQ. 1 .OR. K .LT. 1) GOTO 90
          JK=MBAND*(K-1)+J                                  
          B(I)=B(I)-A(JK)*P1(K)                                                
   90  CONTINUE
      ELSE
        I1=IS(1)
        I2=IS(I+1)-1
        DO 91 J=I1,I2
         K=ISP(J)
   91      B(I)=B(I)-A(J)*P1(K)  
      END IF                                                           
  100  I1=MBAND*(I-1)+1
        IF( NSOLV .EQ. 1 ) I1=IS(I)
       A(I1)=A(I1)+PP*D(I)*PD/DT
C        WRITE(NDCN,*) 'A',I,A(1,I)                                              
  110 CONTINUE                
                                             
      RETURN                                                                
      END                                                                   
C
      SUBROUTINE ELEM (CMUL,KODE,S,B,P,P1,T,D,CR,C,X,Y,Z,
     &   SANGN,CANGN,XK,YK,XP,YP,XC,YC,LM,PORM,C1M,C2M,NMK,NMP,NMC,
     &   NUMNP,NUMMAT,MK,MP,M,KAT,NUMNPZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                                        
      DIMENSION E(3,3),S(4,4),B(NUMNP),P(NUMNP),T(NUMNP),    
     &D(NUMNP),CR(NUMNP),C(NUMNP),X(NUMNP),Y(NUMNP),               
     &P1(NUMNP),XK(MK,NUMMAT),YK(MK,NUMMAT),     
     &XP(MP,NUMMAT),YP(MP,NUMMAT),XC(MP,NUMMAT),YC(MP,NUMMAT),              
     &LM(4),IX(3),KODE(NUMNP),Z(NUMNPZ)                                               
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN 
C----------------------------------- DETERMINE CONDUCTIVITY TENSOR
      Z(1)=0.0  
C      WRITE(NDCN,*)" ENTER ELEM IN RESET"                         
      SIN2=SANGN*SANGN
      COS2=CANGN*CANGN
      CONDI=C1M*COS2+C2M*SIN2                                           
      CONDJ=C1M*SIN2+C2M*COS2                                           
      CONDK=(C1M-C2M)*SANGN*CANGN                              
C      WRITE(NDCN,*)" ENTER do 20-10 IN elem"                          
      DO 20 I=1,4                                                      
      DO 10 J=1,4                                                      
      S(I,J)=0.D0                                                             
   10 CONTINUE                                                              
   20 CONTINUE                                                              
      NUS=4                                                                 
      IF (LM(3) .EQ. LM(4)) NUS=3                                           
C--------------------- DETERMINE UNSATURATED VARIABLES FOR CORNER NODES
      IWET=0                                                                
C      WRITE(NDCN,*)" ENTER do 60 IN elem"                          
      DO 60 K=1,NUS                                                         
      I=LM(K)
      ZZ=0.0D0
C      IF( KAT.LE. 0) ZZ=Z(I)
      IF (CR(I) .GE. 0.) GOTO 50                                           
      IF (P(I)-ZZ .LT. 0.0D00) GOTO 40                                            
C     CR(I)=YK(NMK,M)                                                       
      CR(I)=1.0                                                        
      T(I)=PORM                                                           
      C(I)=0.                                                               
      GOTO 60                                                              
   40 IF (KODE(I) .LT. 1) CM=P(I)-ZZ                                           
      IF (KODE(I) .GE. 1) CM=.5D0*(P(I)+P1(I))-ZZ                                
      CALL INTERP (YP,XP,CM,TI,NMP,NUMMAT,M,MP)                       
      CALL INTERP (XK,YK,TI,CRI,NMK,NUMMAT,M,MK)                    
      CALL INTERP (XC,YC,TI,CI,NMC,NUMMAT,M,MP)
       T(I) = TI
       C(I) = CI
       CR(I) = CRI                     
      IF (C(I) .LT. 0.0D00) C(I) = 0.0D00                                   
      IF (CR(I) .LE. 0.0D00) CR(I)= 1.0D-20                                 
   50 IF (P(I)-ZZ .LT. 0.D0) IWET=1                                              
   60 CONTINUE                                                              
C      WRITE(NDCN,*)" out do 60 IN elem"                          
      QSTORE=0.                                                             
C------------------------------------------------ LOOP ON SUBELEMENTS
      KK=NUS-2                                                            
C      WRITE(NDCN,*)" ENTER do 120 IN elem"                          
      DO 120 K=1,KK                                                         
       I=LM(1)                                                               
       J=LM(K+1)                                                             
       L=LM(K+2)                                                             
        CI=X(L)-X(J)                                                          
        CJ=X(I)-X(L)                                                          
        CK=X(J)-X(I)                                                          
        BI=Y(J)-Y(L)                                                          
        BJ=Y(L)-Y(I)                                                          
        BK=Y(I)-Y(J)                                                          
       XMUL=1.                                                               
      IF (KAT .EQ. 1 . OR . KAT . EQ . 3) THEN                          
       XMUL=(X(I)+X(J)+X(L))/3.D0
      END IF                                              
   80 DEL2=CK*BJ-CJ*BK 
C     IF( DEL2.LE.0.0) WRITE(NDCN,*)   NELM                                                   
       CM=(CR(I)+CR(J)+CR(L))/3.D0                                             
c       CM=(dlog(CR(I))+dlog(CR(J))+dlog(CR(L)))/3.D0
c       CM = dexp(cm)                                             
      COMM=.5D0*XMUL*CM/DEL2                                                  
      QSTORE=.5D0*DEL2*XMUL                                                   
       IF (KAT .GE. 1) THEN                                              
        AA=.5D0*CM*XMUL                                             
        B(I)=B(I)+AA*(CONDK*BI+CONDJ*CI)                                  
        B(J)=B(J)+AA*(CONDK*BJ+CONDJ*CJ)                                  
        B(L)=B(L)+AA*(CONDK*BK+CONDJ*CK)
       END IF                                 
   90 E(1,1)=CONDI*BI*BI+2.*CONDK*BI*CI+CONDJ*CI*CI                         
      E(1,2)=CONDI*BI*BJ+CONDK*(BI*CJ+CI*BJ)+CONDJ*CI*CJ                    
      E(1,3)=CONDI*BI*BK+CONDK*(BI*CK+CI*BK)+CONDJ*CI*CK                    
      E(2,1)=E(1,2)                                                         
      E(2,2)=CONDI*BJ*BJ+2.*CONDK*BJ*CJ+CONDJ*CJ*CJ                         
      E(2,3)=CONDI*BJ*BK+CONDK*(BJ*CK+BK*CJ)+CONDJ*CJ*CK                    
      E(3,1)=E(1,3)                                                         
      E(3,2)=E(2,3)                                                         
      E(3,3)=CONDI*BK*BK+2.*CONDK*BK*CK+CONDJ*CK*CK                         
      IX(1)=1                                                               
      IX(2)=K+1                                                             
      IX(3)=K+2                                                             
      DO 110 I=1,3                                                          
      II=IX(I)                                                              
      IJ=LM(II)                                                             
      BB=0.D0                                                                
      IF (IWET .EQ. 0) BB=T(IJ)*CMUL                                        
      AA=C(IJ)+BB                                                           
      DO 100 J=1,3                                                          
      JJ=IX(J)                                                              
      JI=LM(JJ)                                                             
      BB=0.D0                                                                 
      IF (IWET .EQ. 0) BB=T(JI)*CMUL                                        
      AA=AA+C(JI)+BB                                                        
      S(II,JJ)=S(II,JJ)+E(I,J)*COMM                                         
  100 CONTINUE
      D(IJ)=D(IJ)+AA*QSTORE/12.D0                                            
  110 CONTINUE                                                              
  120 CONTINUE 
      RETURN
      END
      SUBROUTINE ASEMBAG(A,S,LM,MBAND,NUMNP) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)  
      DIMENSION A(MBAND,NUMNP),S(4,4),LM(4)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN                       
C------------------------------ ADD ELEMENT CONTRIBUTION TO A MATRIX 
      DO 140 L=1,4                                                          
      I=LM(L)                                                               
      DO 130 K=1,4                                                          
      J=LM(K)-I+1                                                           
      IF (J .GE. 1) THEN
       A(J,I)=A(J,I)+S(L,K)
      END IF                                                 
  130 CONTINUE                                                              
  140 CONTINUE                                                              
      RETURN                                                                
      END
      SUBROUTINE ASEMBAC(A,S,LM,IS,ISP,NWK,NUMNP1) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)  
      DIMENSION A(NWK),S(4,4),LM(4),IS(NUMNP1),ISP(NWK)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C---- ADD ELEMENT COTRIBUTION TO A MATRIX
      NUS=4
      IF( LM(3).EQ.LM(4) ) NUS=3
      DO 160 L=1,NUS
       I = LM(L)
       I1 = IS(I)
       I2 = IS(I+1) - 1
        DO 150 K=I1,I2
         DO 140 M=1,NUS
          J = LM(M)
          IF( J .EQ. ISP(K) ) GO TO 130
          GO TO 140
  130    A(K) = A(K) + S(L,M)
          GO TO 150
  140    CONTINUE
  150   CONTINUE
  160 CONTINUE
      RETURN
      END                                                                   
      SUBROUTINE SOLVE (A,B, NUMNP,MBAND)                                   
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C 
      DIMENSION A(MBAND,NUMNP),B(NUMNP) 
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN                                    
C  ---- ATOL  -----                                                          
      ATOL = 1.0D-60
      ATL2 = 1.0D-120                                                       
C     REDACTION 
      DO 35 N = 1, NUMNP-MBAND                                             
      DO 25 M=2,MBAND                                                       
      IF ( A(M,N)*A(M,N) .GT. ATL2 ) THEN                                    
       C = A(M,N)/A(1,N)                                                     
       I = N + M - 1                                                         
       J = 0                                                                 
       DO 15 K = M,MBAND                                                     
       J = J + 1                                                             
C        IF ( A(K,N)*A(K,N) .LE .ATL2 ) A(K,N)=0.0D00                                     
        A(J,I)=A(J,I)-C*A(K,N)                                                
   15   CONTINUE                                                              
       A(M,N) = C                                                            
       B(I) = B(I) - A(M,N)*B(N)
      END IF                                             
   25 CONTINUE 
C       WRITE (NDCN,*) N,A(1,N)    
      B(N) = B(N)/A(1,N)                                                    
   35 CONTINUE
C        
      DO 30 N = NUMNP-MBAND+1,NUMNP 
      KK=NUMNP-N+1                                             
      DO 20 M=2,KK                                                       
      IF ( A(M,N)*A(M,N) .GT. ATL2 ) THEN                                    
       C = A(M,N)/A(1,N)                                                     
       I = N + M - 1                                                         
       J = 0                                                                 
       DO 10 K = M,MBAND                                                     
       J = J + 1                                                             
C        IF (A(K,N)*A(K,N) .LE. ATL2 ) A(K,N)=0.0D00                                     
        A(J,I)=A(J,I)-C*A(K,N)                                                
   10   CONTINUE                                                              
       A(M,N) = C                                                            
       B(I) = B(I) - A(M,N)*B(N)
      END IF                                             
   20 CONTINUE 
C       WRITE (NDCN,*) N,A(1,N)    
      B(N) = B(N)/A(1,N)                                                    
   30 CONTINUE     
C----------      B(NUMNP)=B(NUMNP)/A(1,NUMNP)     
C     BACK SUBSTITUTION                                                     
      DO 65 NN=1,MBAND-1
      N=NUMNP-NN 
      KK=NUMNP+1-N
      DO 55 K = 2,KK                                                     
      L = N + K - 1                                                         
      B(N)=B(N)-A(K,N)*B(L)                                                 
   55 CONTINUE                                                              
   65 CONTINUE
C                                                                 
      DO 60 NN=MBAND,NUMNP-1
      N=NUMNP-NN
      DO 50 K = 2,MBAND                                                     
      L = N + K - 1                                                         
      B(N)=B(N)-A(K,N)*B(L)                                                 
   50 CONTINUE                                                              
   60 CONTINUE                                                              
      RETURN                                                                
      END                                                                   
C
      SUBROUTINE CONSTP (A,B,KODE,P,NUMNP,MBAND)                            
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C 
      DIMENSION A(MBAND,NUMNP),B(NUMNP),KODE(NUMNP),P(NUMNP)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN	                  
      DO 70 N=1,NUMNP
      IF (KODE(N) .LT. 1) GOTO 70                                          
      DO 60 M=2,MBAND                                                       
      K=N-M+1                                                               
      IF (K) 30,30,20                                                       
   20 B(K)=B(K)-A(M,K)*P(N)                                                 
      A(M,K)=0.0D00                                                         
   30 L=N+M-1                                                          
      IF (NUMNP-L) 50,40,40                                             
   40 B(L)=B(L)-A(M,N)*P(N)                                            
   50 A(M,N)=0.0D00                                                    
   60 CONTINUE                                                              
      A(1,N)=1.0D00                                                         
      B(N)=P(N)                                                             
   70 CONTINUE                                                              
      RETURN                                                                
      END 
      SUBROUTINE SOLVCG (A,B,P,ISP,IS,R,PP,KODE,RKODET,C,NUMNP,NUMNP1,
     &                   NWK,TSOLV,LA,TIME,ITMAX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*80 COMNT
      DIMENSION A(NWK),B(NUMNP),P(NUMNP),ISP(NWK),IS(NUMNP1),R(NUMNP),
     &          PP(NUMNP),KODE(NUMNP),RKODET(NUMNP),C(NUMNP)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C
      NPMT6 = 8
C ---- ATOL -----
      ATOL = 1.0D-90
      ATOL2 = ATOL*ATOL
      NUMKD=0

C
      NUMKD=0  
      DO 5 I=1,NUMNP
      R(I) = 0.0D0
      PP(I) = 0.0D0
      IF( KODE(I) .GE. 1 ) GO TO 5
       NUMKD=NUMKD+1
       RKODET(NUMKD)=I
      II = IS(I)
      II1 = IS(I+1) - 1
      S = 0.0D0
C-----
       IF( DABS(A(II)) .LE. ATOL ) A(II)=0.0D00
       N = ISP(II)
       AJ=A(II)
       S = S + AJ*P(N)
       TAJ=0.0D0
       II0=II+1
      DO 10 J=II0,II1
      IF( DABS(A(J)) .LE. ATOL ) A(J)=0.0D00
      N = ISP(J)
      AJ=A(J)
      S = S + AJ*P(N)
   10 CONTINUE
C
      R(I) = B(I) - S
      C(I)=A(II)
      PP(I) = R(I)/C(I)
    5 CONTINUE
C
C      WRITE(NDCN,*)' --- START ITERATION'
      IT = 1
   30 PR = 0.0D0
      PAP = 0.0D0
      DO 40 III=1,NUMKD
      I=RKODET(III)
      II = IS(I)
      II1 = IS(I+1) - 1
      PR = PR + PP(I)*R(I)
      BI=0.0D0
      DO 50 J=II,II1
      N = ISP(J)
      BI = BI + A(J)*PP(N)
   50 CONTINUE
      PAP = PAP + PP(I)*BI
      B(I)=BI
   40 CONTINUE
      D = PR/PAP
      RAP = 0.0D0
      EPSX = 0.0D0
      DO 60 III=1,NUMKD
      I=RKODET(III)
      DELTAX = D*PP(I)
      DELTA=DELTAX*DELTAX
      PI=P(I)
      PI=PI*PI
      IF( PI.LE.ATOL2 )PI=1.0D00
      DELTA = DELTA/PI
      EPSX = EPSX + DELTA
   54 P(I) = P(I) + DELTAX
      R(I) = R(I) - D*B(I)
      RAP = RAP + (R(I)/C(I))*B(I)
   60 CONTINUE
      EPSX = DSQRT(EPSX)
C      IF( MOD(IT,100).EQ. 0 ) WRITE(NDCN,*)IT,'    EPSX = ',EPSX
      IF( IT .GT. 10 .AND. EPSX.LE.TSOLV ) GO TO 80
   63 IF( IT - ITMAX ) 65,75,75
   65 E = -RAP/PAP
C
      DO 70 II=1,NUMKD
       I=RKODET(II)
      PP(I) = R(I)/C(I) + E*PP(I)
   70 CONTINUE
      IT = IT + 1
      GO TO 30
C
   75 WRITE (NDCN,600) ITMAX
      WRITE(COMNT,600) ITMAX
C      CALL DISPLY(COMNT)
  600 FORMAT(' *** FAILURE TO CONVERGENCE IN SOLVER -- ITMAX =',I5)
      STOP 999
   80 WRITE (NDCN,660) IT,EPSX
      WRITE(ND6,660) IT,EPSX
  665 FORMAT(I5,',',E15.7,',',I5,',')
      WRITE(COMNT,660) IT,EPSX
C      CALL DISPLY(COMNT)
  660 FORMAT(' *** SUCCESS TO CONVERGENCE IN SOLVER -- IT =',I5,
     1       '    EPSX=',D15.7)
      RETURN
      END
                                                                   
C
      SUBROUTINE FIXQ (A,B,KODE,Q,P,P1,AW,CPI,
     &               NUMNP,MBAND,INTEG,INCOR,INUMNP,NWK,INWK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER*4 ILENR,ILENI
C  
      DIMENSION A(NWK),B(NUMNP),KODE(NUMNP),Q(NUMNP),P(NUMNP),      
     &P1(NUMNP),
     &AW(INWK),CPI(INUMNP,5)
C
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
C      COMMON /CRTL/LA,DT,DT1,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL1/LA,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL2/DT,DT1
  	  COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C----------------------------------------------- READ FROM FILE 2           
      IF( INCOR .EQ. 0 ) THEN
      REWIND ND2
       REC = NWK/ILENR
       IREC = IDNINT(REC)
       IRAM = MOD(NWK,ILENR)
      IF( IREC .LT. 1 ) GO TO 123
      DO I=1,IREC    
      READ(ND2) (A((I-1)*ILENR+J),J=1,ILENR)
      END DO
  123 if( iram .gt. 0) READ(ND2) (A(ILENR*IREC+J),J=1,IRAM)
      READ(ND2) B                                                  
C ----------------------------------------------------------------------
      ELSE
      DO 210 I=1,NUMNP
  210   B(I)=CPI(I,5)
      DO 200 I=1,NWK
      A(I)=AW(I)
  200 CONTINUE
      END IF
C=====================================================================
C
      CALL CALQ (A,B,KODE,Q,P,P1,NUMNP,MBAND,LA,INTEG)
      RETURN
      END
      SUBROUTINE CALQ (A,B,KODE,Q,P,P1,NUMNP,MBAND,LA,INTEG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  
      DIMENSION A(MBAND,NUMNP),B(NUMNP),KODE(NUMNP),Q(NUMNP),P(NUMNP),      
     &    P1(NUMNP)
C
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN   
      PPP=1.                                                               
      PP=2.                                                                
      IF (LA .GT. 1 .AND. INTEG .NE. 1) GOTO 10                           
      PPP=0.                                                               
      PP=1.                                                                
   10 DO 40 N=1,NUMNP                                                      
      IF (KODE(N) .LT.1) GOTO 40                                          
      QN=PP*B(N)+A(1,N)*(P(N)+P1(N)*PPP)                                   
      DO 30 J=2,MBAND                                                      
      K=N-J+1                                                              
      IF (K .LT. 1) GOTO 20                                               
      QN=QN+A(J,K)*(P(K)+P1(K)*PPP)                                        
   20 K=N+J-1                                                              
      IF (K .GT. NUMNP) GOTO 30                                           
      QN=QN+A(J,N)*(P(K)+P1(K)*PPP)                                        
   30 CONTINUE                                                             
      QN=QN/(CONS2*PP)                                                     
      Q(N)=QN                                                              
   40 CONTINUE                                                             
      RETURN                                                               
      END
      SUBROUTINE FIXQPCG ( A,B,KODE,Q,P,P1,ISP,IS,AW,CPI,
     &                  NUMNP,NUMNP1,NWK,INTEG,
     &                  INWK,INUMNP,INCOR )
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER*4 ILENR,ILENI
C
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
C      COMMON /CRTL/LA,DT,DT1,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL1/LA,ITER,NSOLV,ILENR,ILENI
      COMMON /CRTL2/DT,DT1
	COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C
      DIMENSION A(NWK),B(NUMNP),KODE(NUMNP),Q(NUMNP),P(NUMNP),
     1          P1(NUMNP),ISP(NWK),IS(NUMNP1),AW(INWK),CPI(INUMNP,5)
C
C------------ IF()
      IF (INCOR.EQ.0) THEN
       REWIND ND2
       REC = NWK/ILENR
       IREC = IDNINT(REC)
       IRAM = MOD(NWK,ILENR)
      IF( IREC .LT. 1 ) GO TO 123
      DO I=1,IREC    
      READ(ND2) (A((I-1)*ILENR+J),J=1,ILENR)
      END DO
  123 if( iram .gt. 0 ) READ(ND2) (A(ILENR*IREC+J),J=1,IRAM)
      READ(ND2) B                                                  
      ELSE
       DO 12 I=1,NUMNP
   12    B(I)=CPI(I,5)
       DO 13 I=1,NWK
        A(I)=AW(I)
   13  CONTINUE
      END IF
C
       PPP=1.0D00
       PP=2.0D00
      IF(  LA   .GT. 1 .AND. INTEG .NE. 1 )  GO TO 10
       PPP=0.0D00
       PP=1.0D00
C
   10 DO 40 N=1,NUMNP
       IF( KODE(N) .LT. 1 ) GO TO 40
        Q(N)=PP*B(N)
        II = IS(N)
        II1 = IS(N+1) - 1
         DO 30 M=II,II1
          K = ISP(M)
          Q(N)=Q(N)+A(M)*(P(K)+P1(K)*PPP)
   30    CONTINUE
          Q(N)=Q(N)/(CONS2*PP)
   40 CONTINUE
C
      RETURN
      END                                                                  
      SUBROUTINE ELFLOW(X,Y,KX,P,XK,YK,XP,YP,VXX,VYY,H,Z,NUMNP,NUMEL,
     &                   NUMMAT,MK,MP,TIME,LA,IVFLAG,NUMNPZ)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*4  HED
C
      DIMENSION X(NUMNP),Y(NUMNP),KX(NUMEL,5),P(NUMNP),
     & XK(MK,NUMMAT),YK(MK,NUMMAT),XP(MP,NUMMAT),YP(MP,NUMMAT),
     & VXX(NUMEL),VYY(NUMEL),H(NUMNP),Z(NUMNPZ),
     & YY(6),CR(4),
     & NNN(2),VX1(2),VY1(2),V(2),DEG(2) 
C    & ,IB(4),A(4)     
C
      COMMON /FIX/ HED(20),INIT,KTYPE,MAXIT,KAT,INTEG,NSEEP,NCUL,NPRINT,
     &   MAXITO,NTYPE,JMES
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C
      Z(1)=0.0
      DO 800 N=1,NUMEL
      M = KX(N,5)
      IF( M .LE. 0 ) GO TO 801
      NMK = NUMK(M)
      C1M = C1(M)
      C2M = C2(M)
      NMP = NUMP(M)
      NMK = NUMK(M)
      NUS = 4
      IF (KX(N,3) .EQ. KX(N,4)) NUS = 3
      XM = 0.0
      YM = 0.0
      PM = 0.0
	TH=0.0
      ANUS = NUS
      DO 100 J=1,NUS
      I = KX(N,J)
      PM = PM + P(I)/ANUS
      XM = XM + X(I)/ANUS
      YM = YM + Y(I)/ANUS
  100 CONTINUE
      DO 500 K=1,NUS
      I = KX(N,K)
      ZZ=0.0D00
C      IF( KAT .LE. 0) ZZ=Z(I)
      IF (P(I)-ZZ .LT. 0.0) GO TO 550
C     CR(K) = YK(NMK,M)
      CR(K) = 1.0
      GO TO 500
  550 CALL INTERP(YP,XP,P(I)-ZZ,TA,NMP,NUMMAT,M,MP)
      CALL INTERP(XK,YK,TA,CR(K),NMK,NUMMAT,M,MK)
      IF (CR(K) .LT. 0.0) CR(K) = 0.0D0
	TH=TH+TA/ANUS
  500 CONTINUE
      PRM = 0.0
	IF(IVFLAG.EQ.1) TH=1.0
      DO 600 K =1,NUS
  600 PRM = PRM + CR(K)
      PRM = PRM/ANUS
	if( th.le.1.0d-20) th=1.0
      PRM1 = PRM*C1M/TH
      PRM2 = PRM*C2M/TH
      Y1K = XM
      Y2K = YM
      Q1 = 0.0D0
      Q2 = 0.0D0
      D1 = 0.0D0
      J = KX(N,4)
      DO 700 NN=1,NUS
      I = J
      J = KX(N,NN)
      Y1I = X(I)
      Y1J = X(J)
      Y2I = Y(I)
      Y2J = Y(J)
      YY(1) = Y2J - Y2K
      YY(2) = Y2K - Y2I
      YY(3) = Y2I - Y2J
      YY(4) = Y1K - Y1J
      YY(5) = Y1I - Y1K
      YY(6) = Y1J - Y1I
      D = YY(6)*YY(2) - YY(3)*YY(5)
      D1 = D1 + D
      DKAT=D
      IF( KAT.LE.0) DKAT=0.0D00
      Q1 = Q1 - PRM1*(YY(1)*P(I) + YY(2)*P(J) + YY(3)*PM)
      Q2 = Q2 - PRM2*(YY(4)*P(I) + YY(5)*P(J) + YY(6)*PM + DKAT)
  700 CONTINUE
      Q1 = Q1/D1
      Q2 = Q2/D1
      IF ( DABS(Q1) .LT. 1.E-15) Q1=0.0D0
      IF ( DABS(Q2) .LT. 1.E-15) Q2=0.0D0
      VXX(N) = Q1
      VYY(N) = Q2
      GO TO 800
  801 VXX(N)=0.0    
      VYY(N)=0.0
  800 CONTINUE
C------------------------------------------ VELOCITY PRINT OUT FILE NO.6
      IF (IVFLAG .EQ.0) GO TO 539
      WRITE (ND6,610)
      WRITE (ND6,650)HED,TIME,LA,IVFLAG
  610 FORMAT (//,' *** V E L O C I T Y    I N F O R M A T I O N  ***')
  650 FORMAT (//20X,20A4//20X,'TIME =',1PE13.4,5X,'TIME STEP =',I4,
     &              'VELC 1/2=',I5/)
      WRITE(ND6,660)
  660 FORMAT(//,2('   ELEM.     V(X)       V(Z)        VV      DEGREE(f)
     1 ')/)
      J=0
      N=0
   10 N=N+1
      J=J+1
      IF(N .GT. NUMEL) GO TO 20
      NNN(J)=N
      VX1(J)=VXX(N)
      VY1(J)=VYY(N)
      VV=(VXX(N)*VXX(N)+VYY(N)*VYY(N))
      V(J)=DSQRT(VV)
C------------------------------------------- VELOCITY VECTOR DERECTION CHECK  
      XVECT=DABS(VXX(N))
      YVECT=DABS(VYY(N))
C
      IF (XVECT .LE. 1.0D-10 ) THEN 
         IF (VYY(N) .GT. 0.0D0) THEN 
             DEG(J)=90.0D0
          ELSE 
             DEG(J)=-90.0D0
         ENDIF
         GO TO 30 
      ENDIF
C
      IF (YVECT .LE. 1.0D-10 ) THEN 
         IF (VXX(N) .GT. 0) THEN 
             DEG(J)=0.0D0
          ELSE 
             DEG(J)=180.0D0
         ENDIF
      ELSE
       GO TO 35
      ENDIF
      GO TO 30
C 
   35 GRAD=VYY(N)/VXX(N)
      DE=DATAN(GRAD)
      DEG(J)=DE*180.0D0/3.141592654D0
      IF( VX1(J) .LT. 0.D0  .AND. VY1(J) .LT.0.D0) DEG(J)=DEG(J)-180.0D0
      IF( VX1(J) .LT. 0.D0  .AND. VY1(J) .GT.0.D0) DEG(J)=DEG(J)+180.0D0
      GO TO 30
   20 NNN(J)=0
      VX1(J)=0.D0
      VY1(J)=0.D0
      V(J)=0.D0
      DEG(J)=0.D0
   30 IF(J.LT.2) GO TO 10
      WRITE(ND6,666)(NNN(J),VX1(J),VY1(J),V(J),DEG(J),J=1,2)
  666 FORMAT(1X,2(I6,2X,1P3E11.3,0PF8.1,5X))
      J=0
      IF(N.LT.NUMEL) GO TO 10
C------------------------------- WRITE(10) VELOCITY & TOTAL HEAD ------
  539 CONTINUE
      DO 1000 I=1,NUMNP
       ZZ=Y(I)
       IF( KAT .LE. 0) ZZ=0.0D00
       H(I)=ZZ+P(I)
 1000 CONTINUE
C   
 1002 CONTINUE
  999 CONTINUE
      RETURN
      END
C 
      SUBROUTINE INTERP (X,Y,XX,YY,N,M,J,K)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(K,M),Y(K,M)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C---- LINEAR INTERPOLATION
      YY=Y(N,J)   
      DO I=2,N            
      IF (XX .GT. X(I,J) .AND. I .LT. N) GOTO 10                          
      AA=Y(I-1,J)+(XX-X(I-1,J))*(Y(I,J)-Y(I-1,J))/(X(I,J)-X(I-1,J))        
      YY=AA                                                                
      GOTO 20                                                             
   10 CONTINUE
      END DO                                                             
   20 RETURN                                                               
      END
C
      SUBROUTINE INFIL(X,Y,XX,YY,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION X(N),Y(N)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN 
C--------------------------------------LINEAR INTERPOLATION       
      DO 10 I=2,N                                                          
      IF(XX .GT. X(I) .AND. I.LT.N ) GO TO 10                              
      AA=Y(I-1)+(XX-X(I-1))*(Y(I)-Y(I-1))/(X(I)-X(I-1))                    
      YY=AA                                                                
      GO TO 20                                                             
   10 CONTINUE                                                             
   20 RETURN                                                               
      END        
C                                                          
      SUBROUTINE WELBOR (KODE,Q,P,Y,NUMNP,DT)                      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C                                                   
      DIMENSION KODE(NUMNP),Q(NUMNP),P(NUMNP),Y(NUMNP)                     
      COMMON/MAT1/CONS1,CONS2,C1(99),C2(99),SS(99),POR(99),
     &            THR(99),FAN(99)
	COMMON/MAT2/NUMK(99),NUMP(99),NUMC(99)  
      COMMON/WELL1/NB(20),NE(20),NBW,NPT,NPT1,NPB,NPB1,LW,NEP,NSW,KCON
	COMMON/WELL2/AL,AL1,QP,QA,RW,RT,ALFA
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN                                               
C-------------------------------- CALCULATE TOTAL DISCHARGE INTO WELL 
      QA1=QA                                                               
      QA=0.                                                                
      DO 10 I=1,NBW                                                        
      N=NB(I)                                                              
      QA=QA-Q(N)                                                           
   10 CONTINUE                                                             
      DO 20 N=NPB1,NPT                                                     
      QA=QA-Q(N)                                                           
   20 CONTINUE                                                             
      QA = 0.5*(QA+QA1)                                                    
C-------------------------------------- CALCULATE WATER LEVEL IN WELL  
      AA=DT*(QP-QA)/RW*CONS2/CONS1                                         
      AL=ALFA*AL+(1.-ALFA)*(AL1-AA)                                        
C-------------------------------------- DETERMINE CORRESPONDING NODE 
      LW = NPT                                                             
      IF (  KCON  .GE.   1) GOTO 60                                       
      DO 50 N=NPB,NPT                                                      
      IF (AL .GT. Y(N)) GOTO 50                                           
      IF (N .GT. NPB) GOTO 40                                             
      LW=-1                                                                
      GOTO 60                                                             
   40 LW=N-1                                                               
      GOTO 60                                                             
   50 CONTINUE                                                             
C--------------------- MODIFY BOUNDARY CONDITIONS BELOW WATER LEVEL
   60 DO 70 I=1,NBW                                                        
      N=NB(I)                                                              
      IF (KODE(N) .EQ. 0) GOTO 70                                         
      IF (LW .LT. 0) GOTO 65                                              
      P(N) = AL - Y(NPB)                                                   
      GOTO 70                                                             
   65 P(N)=0.                                                              
   70 CONTINUE                                                             
      IF (LW .LT. 0) GOTO 90                                              
      DO 80 N=NPB1,LW                                                      
      IF (KODE(N) .EQ. 0) GOTO 80                                         
      KODE(N)=2                                                            
      P(N)=AL-Y(N)                                                         
   80 CONTINUE                                                             
      IF (  KCON  .GE.   1) GOTO 140                                      
C--------------------- MODIFY BOUNDARY CONDITIONS ABOVE WATER LEVEL 
   90 IF (LW .LT. 0) LW=NPB                                                
      DO 130 J=1,NEP                                                       
      NF=NE(J)                                                             
      NL=NPT                                                               
      IF (J .LT. NEP) NL=NE(J+1)-1                                         
      ICHECK=0                                                             
      DO 120 N=NF,NL                                                       
      IF (N .LE. LW) GOTO 120                                             
      IF (KODE(N) .NE. -2) GOTO 100                                       
      IF (P(N) .LT. 0.D0) ICHECK=1                                           
      IF (ICHECK .GT. 0) GOTO 120                                         
      KODE(N)=2                                                            
   95 P(N)=0.D0                                                              
      GOTO 120                                                            
  100 IF (KODE(N) .NE. 2) GOTO 120                                        
      IF (ICHECK .GT. 0) GOTO 110                                         
      IF (Q(N) .GE. 0.D0) GOTO 110                                          
      GOTO 95                                                             
  110 KODE(N)=-2                                                           
      Q(N)=0.D0                                                              
      ICHECK=1                                                             
  120 CONTINUE                                                             
  130 CONTINUE                                                             
  140 RETURN                                                               
      END                                                                  
       SUBROUTINE PRINTO(Q,P,Y,KPM,Z,NUMNP,KAT,QM,QP,NUMNPZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION Q(NUMNP),P(NUMNP),Y(NUMNP),NN(3),PP(3),HH(3),QQ(3)
      DIMENSION KPM(NUMNP),Z(NUMNPZ)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C
      Z(1)=0.0
      DO 5 I=1,NUMNP
      KPM(I)=0
    5 CONTINUE
C
      WRITE(ND6, 10)                                                      
   10 FORMAT (//  3(' NODE      HEAD    PRESS HEAD    DISCHARGE  ') / )
      J=0                                                                  
      N=0                                                                  
   20 N=N+1                                                                
      J=J+1                                                                
      IF (N .GT. NUMNP) GOTO 35                                           
      NN(J)=N                                                              
      PP(J)=P(N)                                                           
      HH(J)=PP(J)
      ZZ1=Y(N)
      ZZ2=0.0D00                                                          
      IF (KAT .LE. 0) THEN
        ZZ1=0.0D0
C        ZZ2=Z(N)
      END IF
      PP(J)=PP(J)-ZZ2
      HH(J)=HH(J)+ZZ1                                     
      QQ(J)=Q(N)                                                           
      GOTO 45                                                             
   35 NN(J)=0                                                              
      PP(J)=0.                                                             
      HH(J)=0.                                                             
      QQ(J)=0.                                                             
   45 IF (J .LT. 3) GOTO 20                                               
      WRITE(ND6, 50) (NN(J),HH(J),PP(J),QQ(J),J=1,3)                         
   50 FORMAT( 3(I5, 2(1PE12.4),1PE13.5, 2X ) )
      J=0                                                                  
      IF (N .LT. NUMNP) GOTO 20                                           
C------------------------------------------ FLUXES INFORMATION ---------
      LP = -1
      LM =  0
      QP = 0.0
      QM = 0.0
      DO 30 I=1,NUMNP
        IF(DABS(Q(I)).LT.1.0E-20)  GO TO 30
        IF(Q(I).LT.0.0 )  GO TO 31
        LP = LP + 2
        KPM(LP) = I
        QP = QP+Q(I)
        GO TO 30
   31   IF(Q(I).GT.0.0) GO TO 30
        LM = LM + 2
        KPM(LM) = I
        QM = QM + Q(I)
   30 CONTINUE
      LL = MAX0(LP,LM)
CCCCC -----------------------------------------
      DO 40 I=1,LL,2
        IF(MOD(I,100).NE.1) GO TO 71
        WRITE(ND6,610)
   71   CONTINUE
        N1 = KPM(I)
        N2 = KPM(I+1)
CCCCC -----------------------------------------
        IF(N1.LT.1 .OR. N2.LT.1) GO TO 41
        Q1 = Q(N1)
        Q2 = Q(N2)
        WRITE(ND6,611) N1,Q1, N2,Q2
        GO TO 40
   41   IF(N1.LT.1) GO TO 42
        Q1 = Q(N1)
        WRITE(ND6,612) N1, Q1
        GO TO 40
   42   IF(N2.LT.1) GO TO 40
        Q2 = Q(N2)
        WRITE(ND6,613) N2,Q2
   40 CONTINUE
  610  FORMAT(// 5X,'***  B O U N D A R Y   F L U X E S  I N F O R M A T
     & I O N *** '//
     1        10X,5X,'  NODE',9X,'  +Q ',15X,'  NODE',9X,'  -Q '/)
  611 FORMAT(10X,I10, 1PE15.3,10X,I10,1PE15.3 )
  612 FORMAT(10X,I10, 1PE15.3)
  613 FORMAT(45X,I10, 1PE15.3)
      QPQM = QP+QM
      WRITE(ND6,614) QP, QM, QPQM
  614  FORMAT(/2X, 25( '====' ) // 5X,'TOTAL',10X, 1PE15.3,
     1          20X, 1PE15.3,5X,'TOTAL FLUXES = ',1PE10.3// )
CC
      RETURN                                                               
      END                                                                  
C
      SUBROUTINE NEWN(SEEPVAL,KODE,P1,P,Q,Y,NSP,NP,NUMNP,INIT,NDIM,
     1   MAXSP,NSEEP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION KODE(NUMNP),P1(NUMNP),P(NUMNP),Q(NUMNP),Y(NUMNP),
     &          NSP(NDIM),NP(NDIM,MAXSP),SEEPVAL(NDIM)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN	
      WRITE(NDCN,101)
      WRITE(ND6,101)
  101 FORMAT(///' **** NEW B.C. FOR  STEADY ---> NONSTEADY STATE PROBLEM
     & *******'/)
      READ(ND5,20) NCHA,NSEEP
   20 FORMAT(16I5)
      WRITE(NDCN,60) NCHA,NSEEP
      WRITE(ND6,60) NCHA,NSEEP
      IF(NSEEP .LT. 1) GO TO 27
      DO 25 I=1,NSEEP
        READ(ND5,201) NSP(I),SEEPVAL(I)
201     FORMAT(I5,F10.0)
        NS=NSP(I)
        READ(ND5,20)  (NP(I,J),J=1,NS)
   25 CONTINUE
   27 WRITE(NDCN,10)
      WRITE(ND6,10)
      DO 50 I=1,NCHA
        READ(ND5,30)  N,KODE(N),P1(N),Q(N)
        WRITE(NDCN,40) N,KODE(N),P1(N),Q(N)
        WRITE(ND6,40) N,KODE(N),P1(N),Q(N)
   10   FORMAT(//'  NODE NO.',6X,'KODE',5X,' HEAD',14X,'Q'/)
        P1(N) = P1(N) - INIT*Y(N)
        P(N)=P1(N)
   30   FORMAT(2I5,2E10.3)
   40   FORMAT(2I10,2E15.6)
   50 CONTINUE
   60 FORMAT(/// ,'***NUMBER OF VARIABLE NODES***',I4/
     1           ' ***NUMBER OF SEEPAGE FACE*****',I4/)
      RETURN
      END
C
      SUBROUTINE NEXTCO (KODE,P,P1,Q,KX,LA,DT,QTOT,ITER,TIME,LAA,               
     &   NUMNP,NUMEL,DTMAX,DMUL,TMAX,TOL)                                  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION KODE(NUMNP),P(NUMNP),P1(NUMNP),Q(NUMNP),KX(NUMEL,5)
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN      
      WRITE(ND6,5)                                                          
    5 FORMAT(/,' ***** NEW B.C. NONSTEADY ---> NONSTEADY  *****',/)
      READ(ND5, 10)NKOD,NP1,NQ,NMAT                                          
      IF(NKOD.LE.0) GOTO 20                                               
      DO 15 I=1,NKOD                                                       
      READ(ND5, 10)N,KODE(N)                                                 
      WRITE(ND6, 70)N,KODE(N)                                                
   15 CONTINUE                                                             
   20 IF(NP1.LE.0) GOTO 40                                                
      DO 25 I=1,NP1                                                        
      READ(ND5, 30)N,P1(N)                                                   
      WRITE(ND6, 80)N,P1(N)                                                  
   25 CONTINUE                                                             
   40 IF(NQ.LE.0) GOTO 50                                                 
      DO 35 I=1,NQ                                                         
      READ(ND5, 30)N,Q(N)                                                    
      WRITE(ND6, 90)N,Q(N)                                                   
   35 CONTINUE                                                             
   50 IF(NMAT.LE.0) GOTO 60                                               
      DO 45 I=1,NMAT                                                       
      READ(ND5, 10)N,KX(N,5)                                                 
      WRITE(ND6, 100)N,KX(N,5)                                               
   45 CONTINUE                                                             
   60 CONTINUE                                                             
      READ(ND5,110)DT,DTMAX,DMUL,TMAX,TOL                                    
      WRITE(ND6, 120)DT,DTMAX,DMUL,TMAX,TOL                                  
      LA=1                                                                 
      ITER=0                                                               
      LAA=0                                                                
      TIME=DT                                                              
      QTOT=0.0                                                              
      DO 130 I=1,NUMNP                                                     
      P(I)=P1(I)                                                           
  130 CONTINUE                                                             
   10 FORMAT(14I5)                                                         
   30 FORMAT(I5,F10.0)                                                     
   70 FORMAT(/,4X,1HN,8X,7HKODE(N),//,I5,I11)                              
   80 FORMAT(/,4X,1HN,10X,5HP1(N),//,I5,E15.7)                             
   90 FORMAT(/,4X,1HN,6X,4HQ(N),//,I5,E10.3)                               
  100 FORMAT(/,4X,1HN,8X,7HKX(N,5),//,I5,I11)                              
  110 FORMAT(5F10.0)                                                       
  120 FORMAT(/,(8X,2HDT,5X,5HDTMAX,6X,4HDMUL,6X,4HTMAX,7X,                 
     13HTOL),//,5E10.3)                                                    
      RETURN                                                               
      END



                                                                  
C      SUBROUTINE ELEMI (CMUL,SK,BB,PNOD,D,XY,KX,
C     &  XK,YK,XP,    
C     1YP,XC,YC,NUMNP,NUMEL,NUMMAT,MK,MP,M,KAT,N,LA,NINT)                    
C
C     **************************************************************
C     ISOPARAMETRIC FORMULATION OF ELEMENT CONDUCTIVITY AND SPECIFIC
C     MOISTURE CONTENT FOR TWO DIMENSIONAL PLANAR AND AXISYMMETRIC
C     GEOMETRY
C     **************************************************************
C
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C      COMMON /MAT/ CONS1,CONS2,C1(10),C2(10),SS(10),POR(10),NUMK(10),       
C     1NUMP(10),NUMC(10),THR(10),FAN(10)
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN 
C      DIMENSION D(NUMNP),KX(NUMEL,5),XK(MK,NUMMAT),
C     &          YK(MK,NUMMAT),XP(MP,NUMMAT),YP(MP,NUMMAT),
C     &          XC(MP,NUMMAT),YC(MP,NUMMAT),XY(8),PNOD(4),SK(4,4),
C     2          BB(NUMNP)
C      DIMENSION H(4),P(2,4),B(2,4),XJ(2,2),LM(4)
C      DIMENSION  XG(4,4),WGT(4,4)
C      DATA XG /  0.0D0,          0.0D0,          0.0D0,          0.0D0,
C     1 -.57735026919D0,0.57735026919D0,          0.0D0,          0.0D0,
C     2 -.77459666924D0,0.00000000000D0,0.77459666924D0,          0.0D0,
C     3 -.86113631159D0,-.33998104358D0,0.33998104358D0,0.86113631159D0/
C      DATA WGT / 2.0D0,          0.0D0,          0.0D0,          0.0D0,
C     1 1.00000000000D0,1.00000000000D0,          0.0D0,          0.0D0,
C     2 0.55555555555D0,0.88888888888D0,0.55555555555D0,          0.0D0,
C     3 0.34785484513D0,0.65214515486D0,0.65214515486D0,0.34785484513D0/
CC
C      DO 110 I=1,4
C       LM(I)=KX(N,I)
C      DO 110 J=1,4
C        SK(I,J)=0.0D0
C  110 CONTINUE
C       NUS=4
CC      IF (LM(3).EQ.LM(4)) NUS=3
CC
C      IWET=0
CC     LOOP OVER ALL INTEGRATION POINTS
CC
C      DO 100 LX=1,NINT
C      R  = XG(LX,NINT)
C      DO 100 LY=1,NINT
C      S  = XG(LY,NINT)
C      WT = WGT(LX,NINT)*WGT(LY,NINT)
CC
CC     FIND INTERPOLATION FUNCTIONS (H) AND THEIR DARIVATIVES (P).
CC     FIND JACOBIAN (XJ) AND ITS DETERMINANT (DETJ).
CC
C      CALL SHAPE (R,S,XY,H,P,XJ,DETJ,N)
CC
CC     EVALUATE JACOBIAN INVERSE (XJI) AND GLOBAL DERIVATIVE OPERATOR (B)
CC     AT EACH INTEGRATION POINT (R,S) WITHIN THE ELEMENT
CC
C      CALL DERIV (XY,H,P,B,XJ,DETJ,RAD,KAT,N)
CC
C      FAC = WT*RAD*DETJ
CC        WRITE(NDCN,*) 'FAC',FAC
CC
CC     COMPUTE PRESSURE HEAD AT INTEGRATION POINT (R,S)
CC
C      CM = 0.D0
C      DO 10 I=1,4
C   10 CM = CM + PNOD(I)*H(I)
CC      WRITE(NDCN,*) N,CM
CC
CC     BASED ON PRESSURE HEAD AT INTEGRATION POINT (CM), CONPUTE THE
CC     PRESSURE HEAD DEPENDENT CONDCTIVITIES AND SPECIFIC MOISTURE
CC     CONTENT
CC
C      CPEC = 0.D0
C      RPEC=1.0D00
C      TETA=POR(M)
C      IF (CM .GE. 0.D0) GO TO 20
CC
CC       WRITE(NDCN,*)'--- UNSATURATED  --'
C      CALL INTERP (YP,XP,CM,TETA,NUMP(M),NUMMAT,M,MP)
C      CALL INTERP (XK,YK,TETA,RPEC,NUMK(M),NUMMAT,M,MK)
C      CALL INTERP (XC,YC,TETA,CPEC,NUMC(M),NUMMAT,M,MP)
CC
C      IWET =1
CC
C   20 CONTINUE
C      EK11 = RPEC*C1(M)
C      EK22 = RPEC*C2(M)
CC      WRITE(NDCN,*)EK11,EK22,RPEC,C1(M),C2(M)
C      DO 40 I=1,NUS
C      BT1 = B(1,I)
C      BT2 = B(2,I)
C      DO 30 J=1,NUS
C      DB1 = EK11*B(1,J)
C      DB2 = EK22*B(2,J)
CC         WRITE(NDCN,*) BT1,BT2,DB1,DB2,FAC
C   30 SK(I,J) = SK(I,J) + (BT1*DB1 + BT2*DB2)*FAC
C      K2  = LM(I)
C      B1  = 0.D0
C      IF (IWET .EQ. 0) B1=CMUL*TETA
C      AA = CPEC + B1
C      AA = AA*FAC
C      DB3 = EK22*BT2
C      DB3 = DB3*FAC
C      BB(K2) = BB(K2) + DB3
C      D(K2) = D(K2) + AA*H(I)
C        IF( LA.LE. 0) D(K2)=0.0D0
C   40 CONTINUE
C  100 CONTINUE
CC       DO 511 I=1,4
CC       DO 511 J=1,4
CC  511  WRITE(NDCN,*)'SK',I,J,SK(I,J)     
CC
CC      DO 50 I=1,NUS-1
CC      II = I + 1
CC      DO 50 J=II,NUS
CC        WRITE(NDCN,*) N,I,J,SK(I,J)
CC   50 SK(J,I) = SK(I,J)
CC
CC------------------------------ ADD ELEMENT CONTRIBUTION TO A MATRIX 
CC      DO 140 L=1,NUS                                                          
CC      I=LM(L)                                                               
CC      DO 130 K=1,NUS                                                          
CC      J=LM(K)-I+1                                                           
CC      IF (J .LT. 1) GOTO 130                                               
CC      A(J,I)=A(J,I)+SK(K,L)
CC       IF(J .EQ.1 .AND. dabs(A(J,I)).LE. 1.0D-8 ) 
CC     &                   WRITE(NDCN,*) N,I,J,A(J,I)    
C  130 CONTINUE                                                              
C  140 CONTINUE
C      RETURN
C      END
C      SUBROUTINE SHAPE (R,S,XY,H,P,XJ,DETJ,NEL)
C
C     ***********************************************************
C     1. TO FIND INTERPOLATION FUNCTIONS ( H )
C        AND DERIVATIVES ( P ) CORRESPONDING TO THE NODAL POINTS
C        OF A 4- TO 8-NODE TWO DIMENSIONAL ISOPARAMETRIC ELEMENT
C     2. TO FIND JACOBIAN ( XJ ) AND ITS DETERMINANT ( DETJ )
C
C                      NODE NUMBERING CONVENTION
C
C                 2(-1,1)           5(0,1)            1(1,1)
C
C                   0 . . . . . . . 0 . . . . . . . 0
C                   .                               .
C                   .                               .
C                   .               S               .
C                   .               .               .
C                   .               .               .
C                 6 0               . . . R         0 8
C               (-1,0)            (0,0)             .  (1,0)
C                   .                               .
C                   .                               .
C                   .                               .
C                   .                               .
C                   0 . . . . . . . 0 . . . . . . . 0
C
C                 3(-1,-1)          7(0,-1)           4(1,-1)
C
C     ***********************************************************
C
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C      DIMENSION XY(2,4),H(4),P(2,4),XJ(2,2)
C      DIMENSION  IPERM(4)
C      DATA IPERM/2,3,4,1/
C
C     INTERPOLATION FUNCTIONS (4-NODE ELEMENT)
Cc      WRITE(NDCN,*) NEL, (XY(1,I),XY(2,I),I=1,4)
C
C      RP   = 1.0D0 + R
C      SP   = 1.0D0 + S
C      RM   = 1.0D0 - R
C      SM   = 1.0D0 - S
C      R2   = 1.0D0 - R*R
C      S2   = 1.0D0 - S*S
CC
C      H(1) = 0.25D0*RP*SP
C      H(2) = 0.25D0*RM*SP
C      H(3) = 0.25D0*RM*SM
C      H(4) = 0.25D0*RP*SM
CC
CC     LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS (4-NODE ELEMENT)
CC
C      P(1,1) = 0.25D0*SP
C      P(1,2) = -P(1,1)
C      P(1,3) = -0.25D0*SM
C      P(1,4) = -P(1,3)
C      P(2,1) = 0.25D0*RP
C      P(2,2) = 0.25D0*RM
C      P(2,3) = -P(2,2)
C      P(2,4) = -P(2,1)
CC
CC     INTERPOLATION FUNCTIONS AND LOCAL DERIVATIVES FOR MIDSIDE NODES
CC
CC
CC     EVALUATE THE JACOBIAN MATRIX AT POINT (R,S)
CC
C   50 DO 100 I=1,2
C      DO 100 J=1,2
C      SUM = 0.0
C      DO 90 K=1,4
C   90 SUM = SUM + P(I,K)*XY(J,K)
C  100 XJ(I,J) = SUM
CC
CC     COMPUTE THE DETERMINANT OF THE JACOBIAN MATRIX AT POINT (R,S)
CC
C      DETJ = XJ(1,1)* XJ(2,2) - XJ(2,1)* XJ(1,2)
C      DUM = DABS(DETJ)
C      IF(DUM.GT.1.0D-8) GO TO 500
C      WRITE(ND6,3000) NEL
C      STOP
CC
C 3000 FORMAT(//49H **ERROR**  ZERO JACOBIAN DETERMINANT FOR ELEMENT,I5)
CC
C  500 RETURN
C      END
C      SUBROUTINE DERIV (XY,H,P,B,XJ,DETJ,RAD,KAT,NEL)
CC
CC     ******************************************************************
CC     EVALUATION OF THE GLOBAL DERIVATIVE OPERATOR (B) AT A POINT (R,S)
CC     FOR A QUADRILATERAL ELEMENT HAVING PLANAR OR AXISYMMETRIC GEOMETRY
CC     ******************************************************************
CC
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CC
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN 
C      DIMENSION XY(2,4),H(4),P(2,4),B(2,4),XJ(2,2),XJI(2,2)
CC
CC     COMPUTE INVERSE OF THE JACOBIAN MATRIX
CC
C      DETJI    = 1.0D0/DETJ
C      XJI(1,1) =  XJ(2,2)* DETJI
C      XJI(1,2) = -XJ(1,2)* DETJI
C      XJI(2,1) = -XJ(2,1)* DETJI
C      XJI(2,2) =  XJ(1,1)* DETJI
CC
CC     EVALUATE GLOBAL DERIVATIVE OPERATOR  ( B-MATRIX )
CC
C      DO 10 K=1,4
C      B(1,K) = XJI(1,1)*P(1,K) + XJI(1,2)*P(2,K)
C   10 B(2,K) = XJI(2,1)*P(1,K) + XJI(2,2)*P(2,K)
CC
C      RAD = 1.0D0
C      IF (KAT .EQ. 2) GO TO 500
CC
CC     COMPUTE THE RADIUS AT POINT (R,S) FOR AXISYMMETRIC SOLD
CC
C      RAD = 0.0D0
C      DO 50 K=1,4
C   50 RAD = RAD + H(K)* XY(1,K)
CC
C      IF(RAD.GT.1.0D-8) GO TO 500
C      WRITE(ND6,3000) NEL
C      STOP
CC
C 3000 FORMAT(//50H **ERROR**  ZERO RADIUS ENCOUNTERED IN ELEMENT NO.,I5)
CC
C  500 RETURN
C      END
C      SUBROUTINE INFLUE (X,Y,KX,IC,NUMNP,NUMEL)
C
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
C
C      DIMENSION X(NUMNP),Y(NUMNP),KX(NUMEL,5),IC(NUMEL)
C
C      DO 10 I=1,NUMEL
C        IC(I)=0
C       IF( KX(I,5) .LE. 0 ) GO TO 10
C       IF( KX(I,3) .EQ. KX(I,4)) GO TO 10
CC
C       I1=KX(I,1)
C       I2=KX(I,2)
C       I3=KX(I,3)
C       I4=KX(I,4)
CC
C        RX21=X(I2)-X(I1)
C        RX41=X(I4)-X(I1)
C        RY21=Y(I2)-Y(I1)
C        RY41=Y(I4)-Y(I1)
C         RR = RX21*RX41+RY21*RY41
C         IF( DABS(RR) .GT. 1.0D-8 ) GO TO 10
CC
C        RX12=-RX21
C        RX32=X(I3)-X(I2)
C        RY12=-RY21
C        RY32=Y(I3)-Y(I2)
C         RR = RX12*RX32+RY12*RY32
C         IF( DABS(RR) .GT. 1.0D-8 ) GO TO 10
CC
C        RX23=-RX32
C        RX43=X(I4)-X(I3)
C        RY23=-RY32
C        RY43=Y(I4)-Y(I3)
C         RR = RX23*RX43+RY23*RY43
C         IF( DABS(RR) .GT. 1.0D-8 ) GO TO 10
CC
C        RX34=-RX43
C        RX14=-RX41
C        RY34=-RY43
C        RY14=-RY41
CC        RR = RX34*RX14+RY34*RY14
CC        IF( DABS(RR) .GT. 1.0D-8 ) GO TO 10
CC
C          IF( RX21.LT.0.0 .OR. RX41.LT.0.0  .OR.
C     &        RY21.LT.0.0 .OR. RY41.LT.0.0         ) GO TO 100
C           IC(I)=1
C           GO TO 10
C  100     IF( RX32.LT.0.0 .OR. RX12.LT.0.0  .OR.
C     &        RY32.LT.0.0 .OR. RY12.LT.0.0         ) GO TO 200
C           IC(I)=2
C           GO TO 10  
C  200     IF( RX43.LT.0.0 .OR. RX23.LT.0.0  .OR.
C     &        RY43.LT.0.0 .OR. RY23.LT.0.0         ) GO TO 300
C           IC(I)=3
C           GO TO 10  
C  300     IF( RX34.LT.0.0 .OR. RX14.LT.0.0  .OR.
C     &        RY34.LT.0.0 .OR. RY14.LT.0.0         ) GO TO 400
C           IC(I)=4
C           GO TO 10
C  400   WRITE(NDCN,*)' --- INFLUENCE ELEMENT ERROR   ELM= ',I
C        STOP 3333
C   10   CONTINUE      
CC       
C       DO 50 J=1,NUMEL
C       IF( IC(J).LE.0 ) GO TO 50
C         WRITE(ND6,600) J,IC(J),(KX(J,I),I=1,5)
C   50 CONTINUE
C  600 FORMAT('--- ELM = ',6I6,'---') 
C      RETURN
C      END   
C      SUBROUTINE INFMAT (CMUL,SK,BB,PNOD,D,XY,LM,
C     &  XK,YK,XP,YP,XC,YC,
C     &  NUMNP,NUMMAT,MK,MP,M,KAT,N,LA)                    
CC
CC     **************************************************************
C     ISOPARAMETRIC FORMULATION OF ELEMENT CONDUCTIVITY AND SPECIFIC
C     MOISTURE CONTENT FOR TWO DIMENSIONAL PLANAR AND AXISYMMETRIC
C     GEOMETRY
C     **************************************************************
C
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CC 
C      COMMON /MAT/ CONS1,CONS2,C1(10),C2(10),SS(10),POR(10),NUMK(10),       
C     &  NUMP(10),NUMC(10),THR(10),FAN(10)
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN   
C      DIMENSION D(NUMNP),XK(MK,NUMMAT),
C     &          YK(MK,NUMMAT),XP(MP,NUMMAT),YP(MP,NUMMAT),
C     &          XC(MP,NUMMAT),YC(MP,NUMMAT),XY(8),PNOD(4),SK(4,4),
C     &          BB(NUMNP)
C      DIMENSION LM(4)
C      DIMENSION  AXXE(4,4),AYYE(4,4),FYE(4)
C      DATA AXXE/ 2.0D00  ,-2.0D00,  -1.0D00,  1.0D00,
C     &          -2.0D00  , 2.0D00,   1.0D00, -1.0D00,
C     &          -1.0D00  , 1.0D00,   2.0D00, -2.0D00,
C     &           1.0D00  ,-1.0D00,  -2.0D00,  2.0D00/
C      DATA AYYE/ 2.0D00  , 1.0D00,  -1.0D00, -2.0D00,
C     &           1.0D00  , 2.0D00,  -2.0D00, -1.0D00,
C     &          -1.0D00  ,-2.0D00,   2.0D00,  1.0D00,
C     &          -2.0D00,  -1.0D00,   1.0D00,  2.0D00/
C      DATA FYE/ -1.0D00,  -1.0D00,   1.0D00,  1.0D00/ 
CC
C      DO 110 I=1,4
C      DO 110 J=1,4
C        SK(I,J)=0.0D0
C  110 CONTINUE
C      IWET=0
CC
C      RX=XY(3)-XY(1)
C      RY=XY(8)-XY(2)
C       if( rx.le.0.0 .or. ry .le. 0.0 )  WRITE(NDCN,*)' ELM ',N,RX,RY
CC
C      CM = 0.D0
C      DO 10 I=1,4
C   10 CM = CM + PNOD(I)*0.25
CC
C      CPEC = 0.D0
C      RPEC=1.0D00
C      TETA=POR(M)
C      IF (CM .GE. 0.D0) GO TO 20
CC
CC       WRITE(NDCN,*)'--- UNSATURATED  --'
C      CALL INTERP (YP,XP,CM,TETA,NUMP(M),NUMMAT,M,MP)
C      CALL INTERP (XK,YK,TETA,RPEC,NUMK(M),NUMMAT,M,MK)
C      CALL INTERP (XC,YC,TETA,CPEC,NUMC(M),NUMMAT,M,MP)
CC
C      IWET =1
CC
C   20 B1=0.0
C      IF( IWET .EQ. 0 ) B1=TETA*CMUL
C      RAD=1.0
C      IF(KAT.NE.1  .AND. KAT.NE.3 )GO TO 60
C       RAD=(XY(1)+XY(3)+XY(5)+XY(7))*0.25
C   60 AA=CPEC+B1
C      AA=AA*RAD
C      RPEC=RPEC*RAD
C      EK11 = RPEC*C1(M)
C      EK22 = RPEC*C2(M)
C      DO 40 I=1,4
C      DO 50 J=1,4
C      SK(I,J)=EK11*RY/RX*AXXE(I,J) + EK22*RX/RY*AYYE(I,J)
C      SK(I,J)=SK(I,J)/6.0
C   50 CONTINUE 
C      K2=LM(I)
C      IF( LA.GT.0 )D(K2) = D(K2) + RX*RY*0.25*AA
C      BB(K2)=BB(K2) + RX*0.50*EK22*FYE(I)
C   40 CONTINUE
CC------------------------------ ADD ELEMENT CONTRIBUTION TO A MATRIX 
CC      DO 140 L=1,4                                                          
CC      I=LM(L)                                                               
CC      DO 130 K=1,4                                                          
CC      J=LM(K)-I+1                                                           
CC      IF (J .LT. 1) GOTO 130                                               
CC      A(J,I)=A(J,I)+SK(K,L)
CC       IF(J .EQ.1 .AND. dabs(A(J,I)).LE. 1.0D-8 ) 
CC     &                   WRITE(NDCN,*) N,I,J,A(J,I)    
C  130 CONTINUE                                                              
C  140 CONTINUE
C      RETURN
C      END
CC----------------
C      SUBROUTINE ELEMS (CMUL,S,B,P,D,X,Y,Z,
C     &   SANGN,CANGN,XK,YK,XP,YP,XC,YC,LM,C1M,C2M,NMK,NMP,NMC,
C     &   NUMNP,NUMMAT,MK,MP,M,KAT,NUMNPZ)
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CC                                                        
C      DIMENSION E(3,3),S(4,4),B(NUMNP),P(NUMNP),    
C     1D(NUMNP),X(NUMNP),Y(NUMNP),               
C     2XK(MK,NUMMAT),YK(MK,NUMMAT),     
C     3XP(MP,NUMMAT),YP(MP,NUMMAT),XC(MP,NUMMAT),YC(MP,NUMMAT),              
C     4LM(4),IX(3),CRI(2),CIC(2),TI(2),Z(NUMNPZ)               
C      COMMON /MAT/ CONS1,CONS2,C1(10),C2(10),SS(10),POR(10),NUMK(10),       
C     1NUMP(10),NUMC(10),THR(10),FAN(10)
C      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN 
CC----------------------------------- DETERMINE CONDUCTIVITY TENSOR
C        
CC      WRITE(NDCN,*)" ENTER ELEM IN RESET"                         
C      SIN2=SANGN*SANGN
C      COS2=CANGN*CANGN
C      CONDI=C1M*COS2+C2M*SIN2                                           
C      CONDJ=C1M*SIN2+C2M*COS2                                           
C      CONDK=(C1M-C2M)*SANGN*CANGN                              
CC      WRITE(NDCN,*)" ENTER do 20-10 IN elem"                          
C      DO 20 I=1,4                                                      
C      DO 10 J=1,4                                                      
C      S(I,J)=0.D0                                                             
C   10 CONTINUE                                                              
C   20 CONTINUE                                                              
C      NUS=4                                                                 
C      IF (LM(3) .EQ. LM(4)) NUS=3                                           
CC--------------------- DETERMINE UNSATURATED VARIABLES FOR CORNER NODES
C      IWET=0
C      Z1=0.0D0                                                                
C      Z2=0.0D0                                                                
C      Z3=0.0D0                                                                
C      Z4=0.0D0
C      IF( KAT .LE. 0) THEN
C       Z1=Z(LM(1))                                                                
C       Z2=Z(LM(2))                                                                
C       Z3=Z(LM(3))                                                                
C       Z4=Z(LM(4))
C      END IF                                                               
CC      WRITE(NDCN,*)" ENTER do 60 IN elem"
C      PM1 = P(LM(1)) + P(LM(2)) + P(LM(3))-(Z1+Z2+Z3)
C      PM2 = P(LM(1)) + P(LM(3)) + P(LM(4))-(Z1+Z3+Z4)
C      PM1 = PM1/3.0
C      PM2 = PM2/3.0
CC
C      CRI1=1.0D00                                                        
C      CI1=0.0D00 
C      CRI2=1.0D00                                                        
C      CI2=0.0D00 
CC      
C      IF( PM1.LT.0.0D00) THEN 
C        CALL INTERP (YP,XP,PM1,TI1,NMP,NUMMAT,M,MP) 
C        CALL INTERP (XK,YK,TI1,CRI1,NMK,NUMMAT,M,MK)                    
C        CALL INTERP (XC,YC,TI1,CI1,NMC,NUMMAT,M,MP)
C      END IF
C      IF( PM2.LT.0.0D00) THEN 
C        CALL INTERP (YP,XP,PM2,TI2,NMP,NUMMAT,M,MP)  
C        CALL INTERP (XK,YK,TI2,CRI2,NMK,NUMMAT,M,MK)                    
C        CALL INTERP (XC,YC,TI2,CI2,NMC,NUMMAT,M,MP)
C      END IF 
C      IF (CI1 .LT. 0.0D00) CI1 = 0.0D00                                   
C      IF (CI2 .LT. 0.0D00) CI2 = 0.0D00                                   
C      IF (CRI1 .LE. 0.0D00) CRI1= 1.0D-20                                 
C      IF (CRI2 .LE. 0.0D00) CRI2= 1.0D-20                                 
C      IF (CRI1 .GT. 1.0D00) CRI1= 1.0D00                                 
C      IF (CRI2 .GT. 1.0D00) CRI2= 1.0D00  
C       CRI(1)=CRI1                               
C       CRI(2)=CRI2                               
C        CIC(1)=CI1                               
C        CIC(2)=CI2                               
C        TI(1)=TI1                               
C        TI(2)=TI2                               
C   50 IF (PM1 .LT. 0.D0 .OR. PM2.LT.0.0D0 ) IWET=1
C   60 CONTINUE                                                              
C      QSTORE=0.                                                             
CC------------------------------------------------ LOOP ON SUBELEMENTS
C      KK=NUS-2                                                            
CC      WRITE(NDCN,*)" ENTER do 120 IN elem"                          
C      DO 120 K=1,KK                                                         
C       I=LM(1)                                                               
C       J=LM(K+1)                                                             
C       L=LM(K+2)                                                             
C        CI=X(L)-X(J)                                                          
C        CJ=X(I)-X(L)                                                          
C        CK=X(J)-X(I)                                                          
C        BI=Y(J)-Y(L)                                                          
C        BJ=Y(L)-Y(I)                                                          
C        BK=Y(I)-Y(J)                                                          
C       XMUL=1.                                                               
C      IF (KAT .EQ. 1 . OR . KAT . EQ . 3) THEN                          
C       XMUL=(X(I)+X(J)+X(L))/3.D0
C      END IF                                              
C   80 DEL2=CK*BJ-CJ*BK                                                      
C      CM=CRI(K)                                             
C      COMM=.5D0*XMUL*CM/DEL2                                                  
C      QSTORE=.5D0*DEL2*XMUL                                                   
C       IF (KAT .GE. 1) THEN                                              
C        AA=.5D0*CM*XMUL                                             
C        B(I)=B(I)+AA*(CONDK*BI+CONDJ*CI)                                  
C        B(J)=B(J)+AA*(CONDK*BJ+CONDJ*CJ)                                  
C        B(L)=B(L)+AA*(CONDK*BK+CONDJ*CK)
C       END IF                                 
C   90 E(1,1)=CONDI*BI*BI+2.*CONDK*BI*CI+CONDJ*CI*CI                         
C      E(1,2)=CONDI*BI*BJ+CONDK*(BI*CJ+CI*BJ)+CONDJ*CI*CJ                    
C      E(1,3)=CONDI*BI*BK+CONDK*(BI*CK+CI*BK)+CONDJ*CI*CK                    
C      E(2,1)=E(1,2)                                                         
C      E(2,2)=CONDI*BJ*BJ+2.*CONDK*BJ*CJ+CONDJ*CJ*CJ                         
C      E(2,3)=CONDI*BJ*BK+CONDK*(BJ*CK+BK*CJ)+CONDJ*CJ*CK                    
C      E(3,1)=E(1,3)                                                         
C      E(3,2)=E(2,3)                                                         
C      E(3,3)=CONDI*BK*BK+2.*CONDK*BK*CK+CONDJ*CK*CK                         
C      IX(1)=1                                                               
C      IX(2)=K+1                                                             
C      IX(3)=K+2                                                             
C      DO 110 I=1,3                                                          
C      II=IX(I)                                                              
C      IJ=LM(II)                                                             
C      BB=0.D0                                                                
C      IF (IWET .EQ. 0) BB=TI(K)*CMUL                                    
C      AA=CIC(K)+BB                                                           
C      DO 100 J=1,3                                                          
C      JJ=IX(J)                                                              
C      JI=LM(JJ)                                                             
C      BB=0.D0                                                                 
C      IF (IWET .EQ. 0) BB=TI(K)*CMUL                                        
C      AA=AA+CIC(K)+BB                                                        
C      S(II,JJ)=S(II,JJ)+E(I,J)*COMM                                         
C  100 CONTINUE
CC      WRITE(NDCN,*)" set d IN elem"                          
C                                                                  
C      D(IJ)=D(IJ)+AA*QSTORE/12.D0                                            
C  110 CONTINUE                                                              
CC      WRITE(NDCN,*)" out 110 IN elem"                          
C  120 CONTINUE 
C      RETURN
C      END
C
      SUBROUTINE ECHECK ( X,Y,KX,NUMEL,NUMNP,IERSW  )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /DISK/ND98,ND2,ND3,ND4,ND5,ND6,ND10,NDCN
      DIMENSION X(NUMNP),Y(NUMNP),KX(NUMEL,5),LM(4)
C
C        WRITE(*,*)'/// IN ELEM'
C
      IERSW = 0
      DO 130 N=1,NUMEL
      DO 30 I=1,4
      LM(I)=KX(N,I)
   30 CONTINUE
      NUS=4
      IF(LM(3) .EQ. LM(4)) NUS=3
      DO K=1,NUS
      IF(LM(K) .GT. NUMNP) THEN
	   IERSW = 1
	   write(ND6,*) "ERROR ELEMENT =",N," : NODE NUMBER ERROR "
      END IF
      END DO
C
      DO K=1,NUS
      DO KK=1,NUS
      IF(K .NE. KK .AND. LM(K) .EQ. LM(KK)) THEN
	   IERSW = 1
	   write(ND6,*) "ERROR ELEMENT =",N," : NODE NUMBER ERROR "
      END IF
      END DO
      END DO
C
      KK=NUS-2
      DO 120 K=1,KK
      I=LM(1)
      J=LM(K+1)
      L=LM(K+2)
      CI=X(L)-X(J)
      CJ=X(I)-X(L)
      CK=X(J)-X(I)
      BI=Y(J)-Y(L)
      BJ=Y(L)-Y(I)
      BL=Y(I)-Y(J)
C
      DEL2=(CK*BJ-CJ*BL)
	 IF( DEL2.LE.0.0) THEN
	   write(ND6,*) "ERROR ELEMENT =",N," : AREA IS NEGATIVE "
	   IERSW = 1
       END IF
C
  120 CONTINUE
  130 CONTINUE
c      IF(IERSW .EQ. 1 ) STOP
C        WRITE(*,*)'/// OUT ELEM'
      RETURN
      END

