C!MS$OPTIMIZE:'ON:x'
C$DEBUG
C$LARGE
C$STORAGE:2
c      USE DFLIB
C
C-----( AC-GWAP.FOR)-------------------------------------------------
C
C    THIS IS THE PROGRAM 'GWAP' FOR PC-9801 MS-DOS VER.2.11 OR 3.10
C                                   ALL-PC WITH WIN/95
C                             BY DR. M.NISHIGAKI ( OKAYAMA  UNIV. )
C                                    Y.TAKESHITA ( OKAYAMA  UNIV. )
C                                    T.SHIRAISHI ( OKAYAMA  UNIV. )
C                                    Y.SHINSHI   ( -------------- )
C                                    T.IMAI      ( -------------- )
C      ----  HISTORY ----
C
C         OCTORBER ,1979  ...  ORIGINAL DEVELOPED
C         JULY     ,1986  ...  PLOTTER ROUTINE MODIFY
C         SEPTEMBER,1987  ...  PROGRAM CONVERSION
C         JANUARY  ,1988  ...  NEW SOLVER ( PCG )
C         NOVENBER ,1991  ...  INTEGERIZE ETC.
C         DESEMBER ,1991  ...  WELL-RADIUS ETC.
C         OCTOBER  ,2005 
C
C   ----- FILE INFORMATION -----
C
C           * = INPUT FROM KEYBOARD & OUTPUT TO DISPLAY
C      FILE 2 = A(NWK) 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 9 = WRITE RESULTS AT FINAL TIME STEP
C      FILE 10= NODAL POINTS & ELEMENTS DATA FOR MESH CHECK
C      FILE 11= GEOMETORY DATA & FULL RESULTS AT EACH STEP FOR GRAPHICS
C      FILE 15= WRITE RESULTS AT FINAL TIME STEP ON DISK (SUB OUTFD)
C      FILE 16= WRITE RESULTS AT FINAL TIME STEP ON DISK (PLOTT.OUT)
C      FILE 98= WORK & PLOT FILE LIST
C
C   -----------------------------------------------------------------
C
C  ----- PARAMETER -----
C
C       NPMT5  :  5 = STANDARD INPUT
C       NPMT6  :  6 = STANDARD OUTPUT  /  8 = KTACT
C
C  ----- IMPLICIT STATEMENT ----- 
C
C      INPLICIT INTEGER*2 --- *****
C      REAL*4             --- REAL
C      INTEGER*4          --- INTEGER
C      INPLICIT REAL*8    --- DOUBLE PRECISION
C
C   -----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      INTEGER*4 MTOTBD,MTOTIB,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,
C     &          N11,N12,N13,N14,N15,N16,N17,N18,N19,N20,
C     &          N21,N22,N23,N24,N25,N26,N27,N28,N29,
C     &          M1,M2,M3,M4,M5,M6,M7,M8,M9,M10
C
      CHARACTER*80 COMNT
C       REAL HED
      CHARACTER*4 HED
C
      COMMON /FIX/ HED(20),MAXIT,INTEG,NCUL,NPRINT,IPLT
      COMMON /PARAM1/IBOU,NERR,INFL,NPROG,ISTEA,ICONF,CC1,KFLAG
      COMMON /PARAM2/TBEL,IFILE,INCOR
      COMMON /PARAM3/IBOUQ,KBOUMQ,KBMKQ
      COMMON /WELL1/AL(21),AL1(21),AL4(21),NW(21)
      COMMON /WELL2/WARE(21),QOL(21),QL1(21),NWB(21,10)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      COMMON /NPNODE/NPNOD(200)
	COMMON/ IGSP1/ IDNOD1(50),IDNOD2(50),PSURF(50),IECHGE(50,16)
	1   ,IECHGOM(50,16),ADDK,ADDS,ISP1
      COMMON /ADDI1/ IPPSW1
C
C      DIMENSION BD(5600000),IB(2700000)     !check dimension size
      REAL*8 , ALLOCATABLE::BD(:)  !　配列の大きさを変更できるような宣言
      INTEGER , ALLOCATABLE::IB(:)
c      result=SETEXITQQ(QWIN$EXITNOPERSIST)
C
C---- MTOTBD,MTOTBS MUST BE CHANGED FOR EACH PROBLEM
      MTOTBD=56000000   !check dimension size  >(2*NWK+3*MK+23)*(NUMNP or NUMEL)
      MTOTIB=27000000   !check dimension size  >11*(NUMNP or NUMEL)
      NERR=0
C
      NPMT5=5
      NPMT6=8
C
      NPMT2=2
      NPMT3=3
      NPMT4=4
      NPMT9=9
      NPMT10=10
      NPMT11=11
      NPMT98=98
C
C********IFILE
      CALL SETFL(MTOTBD,MTOTIB)
      ALLOCATE( BD(MTOTBD))         !　配列の宣言
      ALLOCATE( IB(MTOTIB))
C
      WRITE(COMNT,26)
   26 FORMAT('--- START ---')
      CALL DISPLY(COMNT)
      CALL TIMCNT(COMNT)
C
      READ (NPMT5,25) HED
   25 FORMAT (20A4)
C
      READ (NPMT5,31) NUMNP,NUMEL,MAXIT,INTEG,NUMMAT,MK,NCUL,
     1        KFLAG,IPLT,TBEL,IFILE,NPRINT,INCOR,NXYP,ISP1,IPPSW1
      NWEL=0
C
      IF( NCUL .LE. 0) NCUL=4
      IF( NPRINT .LE. 0) NPRINT=3
       IF(MK .LE. 1) MK = MK + 1  !
       MK=MK+1
C
   31 FORMAT (9I5,F10.0,4I5,I2,I3)
   30 FORMAT (16I5)
C---- READ HEAD CONDITION
      READ (NPMT5,30)IBOU,KBMK,KBOUM,INFL,ISTEA,ICONF,NT,IRNP,IREL
	1              ,IBOUQ,KBMKQ,KBOUMQ
      NPNP=0
      CLEAK=0.0D00
      PUP=0.0D00
      IF(KBMK .EQ.  0 ) KBMK=1
      IF(KBOUM .EQ. 0 ) KBOUM=1
C
      IF(KBMKQ .EQ.  0 ) KBMKQ=1
      IF(KBOUMQ .EQ. 0 ) KBOUMQ=KBOUM
C
C       IF(ICONF .LE. 0 ) MAXIT=0
C
      NT1=NT
      IF (NT1 .LT. 1 ) NT1=1
C*********************************************************
      N1=1
      N2=N1+NUMEL
      M1=1
      M2=M1+4*NUMEL
      M3=M2+NUMEL
      WRITE (NPMT6,1000) N2,M3
      WRITE (COMNT,1000) N2,M3
      CALL DISPLY(COMNT)
 1000 FORMAT ('***  BD = ',I7,'   *** IB = ',I7)
      IF ( N2-1.GT.MTOTBD  .OR. M3-1.GT.MTOTIB) STOP 1111
C
      CALL MAFIL ( IB(M1),BD(N1),IB(M2),NUMEL,IREL )
C
      NUMNP1 = NUMNP + 1
      M4=M3+NUMNP1
      N3=N2+NUMEL*2-1
C
      WRITE (NPMT6,1000) N3,M4
      WRITE (COMNT,1000) N3,M4
      CALL DISPLY(COMNT)
C
      IF ( N3-1.GT.MTOTBD  .OR. M4-1.GT.MTOTIB) STOP 2222
C
      comnt='   in colht'
      call timcnt(comnt)
      CALL COLHT ( IB(M1),IB(M3),IB(M2),BD(N2),NUMNP,NUMNP1,NUMEL,NWK)
      comnt='   out colht'
      call timcnt(comnt)
C
      M5=M4+NWK
C
      WRITE (NPMT6,1000) N3,M5
      WRITE (COMNT,1000) N3,M5
      CALL DISPLY(COMNT)
      IF ( N3-1.GT.MTOTBD  .OR. M5-1.GT.MTOTIB) STOP 3333
C
      comnt='   in adrsk' 
      call timcnt(comnt)
      CALL ADRSK ( IB(M1),IB(M3),IB(M4),IB(M2),BD(N2),
     &                                    NUMNP,NUMNP1,NUMEL,NWK )
      comnt='   out adrsk'
      call timcnt(comnt)
C
      NUMMAT = NUMMAT + 1
C
C************
      M6=M5+NUMNP
      M7=M6+NUMNP
      M8=M7+NUMMAT
      M9=M8+NUMNP
C
      N3=N2+NWK
      N4 =N3 +NUMNP
      N5 =N4 +NUMNP
      N6 =N5 +NUMNP
      N7 =N6 +NUMNP
      N8 =N7 +NUMNP
      N9 =N8 +NUMNP
      N10=N9 +NUMNP
      N11=N10+NUMNP
      N12=N11+NUMNP
      N13=N12+NUMNP
      N14=N13+MK*NUMMAT
      N15=N14+MK*NUMMAT
      N16=N15+MK*NUMMAT
      N17=N16+KBMK*KBOUM
      N18=N17+KBMK*KBOUM
      N19=N18+KBOUM
      N20=N19+KBOUM
      N21=N20+NUMNP
      N22=N21+NUMNP
      N23=N22+NUMEL
      N24=N23+NUMEL
      N25=N24+NUMNP
      N26=N25+NUMNP
      N27=N26+NT1
C
      NUMNPT=NUMNP
      NWKT=NWK
      IF(INCOR.GT.0) GO TO 6
      NUMNPT=1
      NWKT=1
    6 N28=N27+NUMNPT*4
      N29=N28+NWKT
      N30=N29+NUMNP    !
      N31=N30+KBOUMQ*KBMKQ  !
      N32=N31+KBOUMQ*KBMKQ  !
      N33=N32+NUMEL
      M10=M9 +NUMNPT
C--------------------------------------------------------------------
      WRITE (COMNT,1000) N33,M10
      CALL DISPLY(COMNT)
      WRITE (NPMT6,102) N33,M10
      IF ( N32-1.GT.MTOTBD  .OR. M10-1.GT.MTOTIB) THEN
      COMNT=' !! ERROR !!! MEMORY OVER !! '
      CALL DISPLY(COMNT)
      WRITE (NPMT6,90)
      ELSE
C--- CALCULATION START ----------------------------------------
C************ NUMNPT,NWKT
      comnt='   in fem'
      call timcnt(comnt)
      CALL FEM( BD(N1) ,BD(N2) ,BD(N3) ,BD(N4) ,BD(N5) ,BD(N6) ,
     1          BD(N7) ,BD(N8) ,BD(N9) ,BD(N10),BD(N11),BD(N12),
     2          BD(N13),BD(N14),BD(N15),BD(N16),BD(N17),BD(N18),
     3          BD(N19),BD(N20),BD(N21),BD(N22),BD(N23),BD(N24),
     4          BD(N25),BD(N26),BD(N27),BD(N28),
	1          BD(N29),BD(N30),BD(N31),BD(N32),
     5          IB(M1),IB(M2),IB(M3),IB(M4),IB(M5),IB(M6),
     6          IB(M7),IB(M8),IB(M9),
     8         NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,KBMK,KBOUM,NT,NT1,
     9         NUMNPT,NWKT,NXYP,NWEL,NPNP,CLEAK,PUP,IRNP)
      comnt='  out fem'
      call timcnt(comnt)
C ---------------------------------------------------------------------
      ENDIF
      WRITE (NPMT6,110)
   90 FORMAT (' !! ERROR !!! MEMORY OVER !! ')
  102 FORMAT(/' STORAGE REQUIRD ',' BD= ',I7, ' IB= ',I7)
  110 FORMAT (/,'  ****  NORMAL TERMINATION **** ')
C
      STOP 
      END
C***************************************************
C***************************************************
C$LARGE
C$STORAGE:2
      SUBROUTINE FEM( ALF,A,B,Q,P,P1,D,CR,C,X,Y,QL,
     1                BK,TK,SK,HB,TIM,EI,ETIM,QO,BL,VXX,VYY,PPO,
     2                DEF,TM,PT4,AT2,
	1                QQN,TIMEQ,QB,CMPP,
     3                KX,KM,IS,ISP,KODE,KB,NUMK,KODET,KODET4,
     4          NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,KBMK,KBOUM,NT,NT1,
     5          NUMNPT,NWKT,NXYP,NWEL,NPNP,CLEAK,PUP,IRNP)
C**************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
C      REAL HED
      CHARACTER*4 HED
      CHARACTER*80 COMNT
      COMMON /FIX/ HED(20),MAXIT,INTEG,NCUL,NPRINT,IPLT
      COMMON /PARAM1/IBOU,NERR,INFL,NPROG,ISTEA,ICONF,CC1,KFLAG
      COMMON /PARAM2/TBEL,IFILE,INCOR
      COMMON /PARAM3/IBOUQ,KBOUMQ,KBMKQ
      COMMON /WELL1/AL(21),AL1(21),AL4(21),NW(21)
      COMMON /WELL2/WARE(21),QOL(21),QL1(21),NWB(21,10)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      COMMON /NPNODE/NPNOD(200)
	COMMON/ IGSP1/ IDNOD1(50),IDNOD2(50),PSURF(50),IECHGE(50,16),
	1   IECHGOM(50,16),ADDK,ADDS,ISP1
C
      DIMENSION ALF(NUMEL),A(NWK),B(NUMNP),Q(NUMNP),P(NUMNP),
     1                P1(NUMNP),D(NUMNP),CR(NUMNP),C(NUMNP),
     2                X(NUMNP),Y(NUMNP),QL(NUMNP),BK(MK,NUMMAT),
     3                TK(MK,NUMMAT),SK(MK,NUMMAT),HB(KBOUM,KBMK),
     3                TIM(KBOUM,KBMK),EI(KBOUM),ETIM(KBOUM),QO(NUMNP),
     4                BL(NUMNP),VXX(NUMEL),VYY(NUMEL),PPO(NUMNP),
     4                DEF(NUMNP),TM(NT1),PT4(NUMNPT,4),AT2(NWKT),
     6                KX(NUMEL,4),KM(NUMEL),IS(NUMNP1),ISP(NWK),
     7                KODE(NUMNP),KB(NUMNP),
     8                NUMK(NUMMAT),KODET(NUMNP),KODET4(NUMNPT),
     &           NN(3),QN(3),PN(3),BLN(3)
      DIMENSION QQN(NUMNP),TIMEQ(KBOUMQ,KBMKQ),QB(KBOUMQ,KBMKQ)
      DIMENSION CMPP(NUMEL)
C
c      COMNT='--- IN FEM ---'
c      CALL TIMCNT(COMNT)
CC--------------------------------------------------------------------
      EPSLON=0.
C----- READ TIME STEP DATA INFORMATION ----------------------------------
      READ (NPMT5,90)DT,DTMAX,DMUL,TMAX,TOL,EIFL,TSOLV
       TSOLV= DABS(TSOLV)
      IF( TSOLV .LE. 0.0 ) TSOLV=1.0D-8
      IF ( NT .LE. 0 ) GO TO 94
      READ (NPMT5,90)(TM(I),I=1,NT)
   90 FORMAT (7F10.0)
C----- PRINT GENERAL INFORMATION ----------------------------------
   94 CONTINUE
C*****************  WELL INFORMATION
C       IF(NWEL.LE.0) GO TO 89
C       DO 96 I=1,NWEL
C       READ (NPMT5,97)NW(I),AL(I),WARE(I)
C       READ (NPMT5,99) (NWB(I,N),N=1,NW(I))
C       AL1(I)=AL(I)
C       WRITE (NPMT6,970) I,NW(I),AL(I),WARE(I)
C       WRITE (NPMT6,990) (NWB(I,N),N=1,NW(I))
C   96 CONTINUE
C   97 FORMAT(I5,2E10.3)
C   99 FORMAT(16I5)
C  970 FORMAT('* WELL NO.= ',I5,'  NUMB.OF .NOD = ',I5,
C     &       ' WAT-LEVEL =',E15.7,'  AREA= ',E15.7/)
C  990 FORMAT(16I5/)
C*************** DATA OF OUTPUT NODE
C   89 IF(NPNP.LE.0) GOTO 98
C      READ (NPMT5,87)(NPNOD(I),I=1,NPNP)
C   87 FORMAT(16I5)
C      WRITE(NPMT6,'(16I5)')(NPNOD(I),I=1,NPNP)
C*************************************
   98 WRITE (NPMT6,100)
  100 FORMAT (//)
      WRITE (NPMT6,101)
  101 FORMAT(/' *** PC-GWAP ** QUASI-3D ** HORIZONTAL PLANE FLOW ***'/)
C
      IF ( ISTEA .EQ.  0 ) WRITE (NPMT6,9010)
      IF ( ISTEA .EQ.  1 ) WRITE (NPMT6,9020)
      IF ( ISTEA .EQ. -1 ) WRITE (NPMT6,9030)
      IF ( ISTEA .EQ. -2 ) WRITE (NPMT6,9040)
      IF ( ISTEA .EQ.  2 ) WRITE (NPMT6,9050)
 9010 FORMAT ( /' *** INIT.STEADY ---> NON-STEADY ---> FINAL STAEDY STAT
     &E PROBLEM ***'//)
 9020 FORMAT ( /' *** NON-STEADY ---> FINAL STAEDY STATE PROBLEM ***'//)
 9030 FORMAT ( /' *** ONLY STEADY STATE PROBLEM ***'//)
 9040 FORMAT ( /' *** INIT.STEADY ---> NON-STEADY STATE PROBLEM ***'//)
 9050 FORMAT ( /' *** ONLY NON-STEADY STATE PROBLEM ***'//)
C
      IF ( ICONF .EQ. 0 ) WRITE (NPMT6,9100)
      IF ( ICONF .EQ. 1 ) WRITE (NPMT6,9200)
 9100 FORMAT ( /' ======  CONFINED AQUIFER PROBLEM ======'//)
 9200 FORMAT ( /' ******  UNCONFINED AQUIFER PROBLEM *****'//)
C
      WRITE (NPMT6,140) HED,NUMNP,NUMEL,NUMMAT-1,MAXIT,TBEL
  140 FORMAT (6HTITLE ,20A4//30H NUMBER OF NODAL POINTS-------  I4/
     130H NUMBER OF ELEMENTS-----------  I4/
     230H NUMBER OF MATERIALS----------  I4/
     330H MAX NUMBER OF ITERATIONS-----  I4/
     430H TBEL (TP FROM BASE-LINE)-----  E10.3/)
C
      WRITE (NPMT6,142)IBOU,KBMK,KBOUM
  142 FORMAT(/' INDEX OF BOUNDARY VARIABLE HEAD--- ' I3/
     1        ' NUMBER OF VARIABLE NODE----------- 'I3/
     2        ' NUMBER OF FUNCTION---------------- 'I3/)
C
      WRITE (NPMT6,143)DT,DTMAX,DMUL,TMAX,TSOLV,NT,NPRINT,TOL,EIFL
  143 FORMAT(/' INITIAL TIME----------- '1PE10.3/
     1        ' MAXIMUM TIME----------- 'E10.3/
     2        ' MULTIPLY FACTOR-------- 'E10.3/
     3        ' FINAL TIME------------- 'E10.3/
     3        ' SOLV TOLERANCE--------- 'E10.3/
     4        ' PRINT OUT STEP-(NT)---- 'I5/
     5        ' PRINT STEP-(NPRINT)---- 'I5/
     6        ' MINIMUM ERROR---------- 'E10.3/
     7        ' INFILTRATION RATE------ 'E10.3/)
C
      IF ( IFILE .GT. 0 ) WRITE (NPMT6,95)
   95 FORMAT ( //,' **** FINAL TIME STEP HEAD & NODE INFORMATION WRITE
     &TO FILE NO.9 *****'//)
      WRITE (NPMT6,91)
   91 FORMAT (' ** PRINT OUT TIME *** ')
      WRITE (NPMT6,92)(TM(I),I=1,NT)
   92 FORMAT (1P10E10.3)
C---- READ AND GENERATE MATERIAL INFORMATION ---------------------
      CALL MATIN ( BK,TK,SK,NUMK,MK,NUMMAT,CC1)
C---- READ AND GENERATE BOUNDARY CONDITION ------------------------
      IF(IBOU .NE. 1) GO TO 44
      CALL BOUND ( HB,TIM,KBOUM,KBMK,1)
C---- READ INFILTRATION RATE --------------------------------------
   44 IF (INFL .NE. 1) GO TO 45
      CALL EINFI (EI,ETIM,KBOUM)
C---- PRINT ELEMET DATA ----
C************ KM(N)
   45 CONTINUE
      IF(IBOUQ .NE. 1) GO TO 46
      CALL BOUND ( QB,TIMEQ,KBOUMQ,KBMKQ,2)
   46 CONTINUE
C      WRITE (NPMT6,10)
C   10 FORMAT (//' ** ELEMENT INFORMATION **'//
C     1  '       ELEMENT    C O R N E R  N O D E S      INFIL.FACTOR'/)
C      DO 80 N=1,NUMEL
C      WRITE (NPMT6,70)N,KX(N,1),KX(N,2),KX(N,3),KX(N,4),KM(N),ALF(N)
C   70 FORMAT (4X,I6,5I6,6X,E10.3)
C   80 CONTINUE
C---- READ AND GENERATE NODAL POINT INFORMATION ------------------
C************-KM
      CALL NPIN(KODE,QO,P,P1,X,Y,KB,QL,KM,KODET,KX,NUMNP,NUMEL,
     1          NERR,BL,NXYP,IRNP)
C
C
      IF(ISP1 .NE. 0) THEN
C
      DO I=1,ISP1
	READ(NPMT5,3400) IDNOD1(I),PSURF(I),IDNOD2(I)
	READ(NPMT5,3401) (IECHGE(I,J),J=1,16)
 3400 FORMAT(I5,F10.0,I5)
 3401 FORMAT(16I5)
	END DO
	READ(NPMT5,3402) ADDK,ADDS
 3402 FORMAT(2F10.0)
C
	DO I=1,ISP1
	 DO II=1,16
	  IEE1 = IECHGE(I,II)
	  IF(IEE1 .NE. 0) THEN
	  IECHGOM(I,II)=KM(IEE1)   
	  KM(IEE1)=NUMMAT
	  END IF
       END DO
	END DO
C
      NUMK(NUMMAT) = 3
      BK(1,NUMMAT)=0.0
      TK(1,NUMMAT)=0.0
      SK(1,NUMMAT)=0.0
      BK(2,NUMMAT)=1.0
      TK(2,NUMMAT)=ADDK
	SK(2,NUMMAT)=ADDS
      BK(3,NUMMAT)=10000.0
      TK(3,NUMMAT)=0.0
      SK(3,NUMMAT)=ADDS
C
      END IF
C
      IF (NERR .GT. 0 ) RETURN
C------------------------------------------------------------------
      COMNT=' ALL INPUT DATA HAVE BEEN READ !'
      CALL DISPLY(COMNT)
C 710 FORMAT (//' ALL INPUT DATA HAVE BEEN READ !')
C
      COMNT='--- READ OK ---'
      CALL TIMCNT(COMNT)
C
C--------------------------------- MESH DATA INFORMATION FOR PLOTTER
C
      WRITE(NPMT9) HED
      WRITE(NPMT9) NUMNP
      WRITE(NPMT9) (X(I),Y(I),QL(I),BL(I),KB(I),I=1,NUMNP)
      WRITE(NPMT9) NUMEL
      WRITE(NPMT9) (KX(I,1),KX(I,2),KX(I,3),KX(I,4),
     &                   KM(I),ALF(I),I=1,NUMEL)
C-------------------------- INPUT DATA CONTOUR INFORMATION FOR PLOTTER
      IERSW = 0
      CALL  ECHECK ( X,Y,KX,NUMEL ,NUMNP ,IERSW)
      IF(IERSW .EQ. 1) GO TO 1121
C
      IF( IPLT.NE.1) GO TO 1100
C
      COMNT=' **** MESH CHECK ONLY!! NORMAL TERMINATION *****'
      CALL DISPLY(COMNT)
      WRITE (NPMT6,1000)
 1000 FORMAT(' **** MESH CHECK ONLY!! NORMAL TERMINATION *****')
 1121 CONTINUE
      LADD = 0
      TIMEDD=0.0
      CMPPD =1.0
      WRITE(NPMT9) LADD,TIMEDD
      WRITE(NPMT9) (KODE(I),P(I),Q(I),I=1,NUMNP)
      WRITE(NPMT9) (VXX(I),VYY(I),CMPPD,I=1,NUMEL)
      RETURN
C----- SET UP FOR TIME INTEGRATION ------------------------------
 1100 ICHE=0
      LAP=0
      IEND=0
  231 CONTINUE
      IFIN=0
      LA=0
      TIME=0.
      ITER=0
C----------------------------------------------------------
      LAA=0
      NNT=1
      IF (ISTEA .LE. 0 .AND. IEND .EQ. 0 ) GO TO 230
      LA=1
      TIME=DT
      IEND=0
C----------------------------------------------------------
  230 WRITE (NPMT6,239)    !<---------------------------
  239 FORMAT (//)
      WRITE (NPMT6,240)HED,TIME,LA
      WRITE (COMNT,831)HED
  831 FORMAT(20A4)
      CALL DISPLY(COMNT)
      WRITE (COMNT,832) TIME ,LA
  832 FORMAT(' TIME=',1X,1PE13.5,' STEP= ',1X,I4)
      CALL DISPLY(COMNT)
  240 FORMAT (//20X,20A4/20X,6HTIME =1PE13.5,5X,11HTIME STEP =I4/)
      IF( IFIN .GE. 1) WRITE (*,235)
      IF( IFIN .GE. 1) WRITE (NPMT6,235)
  235 FORMAT(' ***** FINAL STEADY STATE *****'/)
C---- HEAD DATA WRITE TO  OR READ FROM FILE (4) -----------
      IF(LAA .LT. 1) GO TO 1
      GO TO 2
C************ IF()
    1 IF(INCOR.GT.0) GO TO 11
         REWIND NPMT4
         WRITE (NPMT4) P
         WRITE (NPMT4) P1
         WRITE (NPMT4) Q
         WRITE (NPMT4) KODE
         WRITE (NPMT4) AL
       GO TO 3
   11 DO 12 I=1,NUMNP
         PT4(I,1)=P(I)
         PT4(I,2)=P1(I)
         PT4(I,3)=Q(I)
         KODET4(I)=KODE(I)
   12 CONTINUE
C      DO 821 N=1,NWEL+1
C  821    AL4(N)=AL(N)
      GO TO 3
    2 IF(INCOR.GT.0) GO TO 21
      REWIND NPMT4
      READ (NPMT4) P
      READ (NPMT4) P1
      READ (NPMT4) Q
      READ (NPMT4) KODE
      READ (NPMT4) AL
      GO TO 23
   21 DO 22 I=1,NUMNP
      P(I)=PT4(I,1)
      P1(I)=PT4(I,2)
      Q(I)=PT4(I,3)
      KODE(I)=KODET4(I)
   22 CONTINUE
C      DO 811 N=1,NWEL+1
C  811 AL(N)=AL4(N)
   23 WRITE (NPMT6,5)LAA
    5 FORMAT(//20X,6HRECUL=I5//)
C----------------------------------------------------------
    3 EPSLN1=1.0D28
      EPSLN2=1.0D28
      EPSLON=1.0D28
      IF( LA .LT. 1 .OR. IFIN .GE. 1) GO TO 15
C---- SET UP NEW TIME BOUNDARY HEAD OR DISCAHRGE CONDITION ----------
C****  WELL
       IF(NWEL.GT.0) CALL WELBOR(P,QL,Q,QO,NWEL,DT,NUMNP)
      IF (IBOU  .EQ. 0 ) GO TO 14
      DO 13 I=1,NUMNP
       IF(.NOT.(KODE(I) .EQ. 3 .OR. IABS(KODE(I)) .EQ. 13)) GO TO 13
C      IF (KODE(I) .NE. 3 ) GO TO 13
CCC      IF (KODE(I) .NE.-1 ) GO TO 13
      M=KB(I)
C      CALL INTERP (TIM,HB,TIME,P(I),KBOUM,KBMK,M,KBOUM)
      CALL INTERP (TIM,HB,TIME,PPDUM,KBOUM,KBMK,M,KBOUM)
c	write(88,*) I,M,KODE(I),BL(I),PPDUM
CCC      CALL INTERP (TIM,HB,TIME,QO(I),KBOUM,KBMK,M,KBOUM)
	IF(IABS(KODE(I)) .EQ. 13) THEN
	IF(BL(I) .GT. PPDUM)  THEN !水位が地表より低い
	 KODE(I) =-13
	ELSE
	 KODE(I) =13
       P(I)=PPDUM
	END IF
	END IF 
	IF(KODE(I) .EQ. 3 ) P(I)=PPDUM
   13 CONTINUE
C---- SET UP NEW INFILTRATION --------------------------------------
   14 IF (INFL .NE. 1 ) GO TO 15
      CALL INTER3 (ETIM,EI,TIME,EI1,KBOUM)
      GO TO 16
   15 EI1=EIFL
   16 CONTINUE
C
      IF(IBOUQ.NE. 0) THEN
	 TTTT = TIME - DT/2.0
	 DO I=1,NUMNP
	 QQN(I) = 0.0
	  IF(KODE(I) .NE. -6) GO TO 24
	    M=KB(I)
          CALL INTERP (TIMEQ,QB,TTTT,QQN(I),KBOUMQ,KBMKQ,M,KBOUMQ)
   24  CONTINUE
       END DO
	END IF
C
      IF(ISP1 .NE. 0) THEN
      DO I=1,ISP1
      NN1 = IDNOD1(I)
	NN2 = IDNOD2(I)
c	IF(P(NN1) .GT. P(NN2) .AND. P(NN1) .GE. PSURF(I) ) THEN
 	IF(P(NN1) .GT. P(NN2) ) THEN
C
	 DO II=1,16
	  IEE1 = IECHGE(I,II)
	  IF(IEE1 .NE. 0) THEN
	  KM(IEE1) = IECHGOM(I,II)
	  END IF
       END DO
C
	ELSE
C
	 DO II=1,16
	  IEE1 = IECHGE(I,II)
	  IF(IEE1 .NE. 0) THEN
	  KM(IEE1) = NUMMAT
	  END IF
       END DO
C
	END IF
      END DO
	END IF
C
C----- GENERATE TERMS OF MATRIX  EQUATION ---------------------------
C************ AT2,NUMNPT,NWKT,INCOR(RESET,FIXQ)
  250 CONTINUE   !<---------------------------
C
      COMNT='--- IN RESET ---'
      CALL TIMCNT(COMNT)
C
c       IF(LA .LE. 0) WRITE(78,*) KODE
c       IF(IFIN .EQ. 1) WRITE(79,*) KODE
C
      CALL RESET ( KODE,A,B,Q,P,P1,D,CR,C,X,Y,KX,BK,TK,SK,NUMK,KM,ALF,
     1             ISP,IS,AT2,NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,LA,
     2             DT,DT1,ITER,INTEG,EI1,IFIN,BL,QL,QO,TBEL,ICONF,
     3             NWKT,INCOR,CLEAK,PUP,QQN)
C---- SOLVE BY PRECONDITIONED CONJUGATE GRADIENT METHOD (P.C.G.M.)----
C
      COMNT='--- OUT RESET IN SOLV ---'
      CALL TIMCNT(COMNT)
C
      CALL SOLVCG(A,B,P,ISP,IS,D,CR,KODE,KODET,C,NUMNP,NUMNP1,NWK,TSOLV)
C------ DETERMINE BOUNDARY FLUXES --------------------------------
C
      COMNT='--- IN FIXQ ---'
      CALL TIMCNT(COMNT)
C
      CALL FIXQ( A,B,KODE,Q,P,P1,ISP,IS,AT2,
     1          NUMNP,NUMNP1,NWK,LA,INTEG,IFIN,NUMNPT,NWKT,INCOR )
C---- ITERATE TO IMPROVE SOLUTION --------------------------------
      IF(MAXIT.LT. 1 ) GO TO 460
      ITER=ITER+1
       IF(NWEL.GT.0) CALL WELBOR(P,QL,Q,QO,NWEL,DT,NUMNP)
C------- MODIFY CONDTIONS ON SURFACE -------------------------------
      IF( ICONF.LT.1) GO TO 266
      DO 265 I=1,NUMNP
        IF(KODE(I).EQ.1) GO TO 265
        IF(KODE(I).EQ.3) GO TO 265
C
        IF(IABS(KODE(I)).EQ. 13) GO TO 265   !      
C
        IF(KODE(I).EQ.-1) GO TO 265
        IF(KODE(I).EQ.5) GO TO 265
        IF(KODE(I).EQ.4) GO TO 264
        IF(KODE(I).EQ.-6) GO TO 264
        IF(P(I).LE.BL(I)) GO TO 264
        P(I)=BL(I)
        KODE(I)=4
       GO TO 265
  264 IF(KODE(I).NE.4) GO TO 265
      IF(Q(I).LT.0.) GO TO 265
        KODE(I)=0
        Q(I)=0.
  265 CONTINUE
  266 CONTINUE
      DO 267 I=1,NUMNP
       IF( KODE(I) .NE. -3  .AND. KODE(I).NE.33 ) GO TO 267
       IF( KODE(I) .NE. -3 ) GO TO 268
          IF( P(I).GE.QL(I)) GO TO 267
          P(I)=QL(I)
          KODE(I)=33
          GO TO 267
  268 IF( Q(I).GT.QO(I)) GO TO 267
        KODE(I)=-3
  267  CONTINUE
C---- TEST FOR CONVERGENCE ---------------------------------------
      IF (ITER .LE. 2 ) GO TO 450
C************ IF()
      IF (INCOR.GT.0) GO TO 31
      REWIND NPMT3
      READ (NPMT3) C
      GO TO 33
   31 DO 32 I=1,NUMNP
      C(I)=PT4(I,4)
   32 CONTINUE
   33 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
      WRITE (NPMT6,430)ITER,EPSLON,NMAX
      WRITE ( COMNT,431)ITER,EPSLON,NMAX
      CALL DISPLY(COMNT)
  430 FORMAT(/49H MAXIMUM CHANGE IN PRESSURE HEAD DURING ITERATION,I3,
     1  4H WAS,E13.5,8H AT NODE,I5)
  431 FORMAT('--- ITER =',I3,' EPS= ',E13.5,' AT NODE=',I5)
  440 IF (EPSLON .LE. TOL) GO TO 460
C      IF ( ITER .GT. MAXIT .OR. EPSLN2.LE. EPSLON ) GO TO 461
      IF ( ITER .GT. MAXIT  ) GO TO 461
C************ IF
  450 IF (INCOR.GT.0) GO TO 35
      REWIND NPMT3
      WRITE (NPMT3) P
      GO TO 37
   35 DO 36 I=1,NUMNP
      PT4(I,4)=P(I)
   36 CONTINUE
   37 EPSLN2=EPSLN1
      EPSLN1=EPSLON
      GO TO 250
C----------------------------------------------------------
  461 IF( LA .LT. 1 .OR. IFIN .GE. 1) GO TO 569
      TIME=TIME-DT
      DT=.5*DT
      TIME=TIME+DT
      ITER=0
      LAA=LAA+1
      IF( LAA .GT. NCUL ) GO TO 560
      GO TO 230
C --------------------------------------------------------
  460 CONTINUE
C ------------------------------------------------------ PRINT RESULTS
      IF (LA .LT. 1 .OR. IFIN .GE. 1 ) GO TO 471
      IF (NT .LE. 0 ) GO TO 468
      IF (DABS(TIME-TM(NNT)) .GT. 1.0D-03*DT) GO TO 468
      NNT= NNT+1
      GO TO 471
C -----------------------------------------------------------------
  468 IF (DABS(TIME-TMAX) .LE. .001*DT ) GO TO 471
      K = MOD( LA,NPRINT )
      IF ( K. NE. 0) GO TO 482
C------------------------------------------- ALL RESULTS TO NO.6 FILE
  471 CALL PRINTO (Q,P,NUMNP,TBEL,KODE,NPNP)
C------------------------------------------- VECTOR INFORMATION
      CALL VECTOR ( P,X,Y,KX,BK,TK,NUMK,KM,NUMNP,NUMEL,
     1           NUMMAT,MK,BL,VXX,VYY,TBEL,CMPP,KFLAG )
C --------------------------------- WRITE PLOTTER DATA TO NO.11 FILE
      WRITE(NPMT9) LA,TIME
      WRITE(NPMT9) (KODE(I),P(I),Q(I),I=1,NUMNP)
      WRITE(NPMT9) (VXX(I),VYY(I),CMPP(I),I=1,NUMEL)
C ----------------------------- HEAD DIFFERENCE FROM STEADY STATE -
      IF ( ISTEA .GE. 1 .OR. ISTEA .EQ. -1 ) GO TO 480
      IF ( LA .NE. 0 ) GO TO  800
      DO 801 I=1 , NUMNP
      PPO(I)=P(I)
  801 CONTINUE
  800 DO 802 I=1,NUMNP
      DEF(I)=P(I)-PPO(I)
  802 CONTINUE
      IF ( LA .LT.1 ) GO TO 480
      CALL PRINT1 (DEF,NUMNP)
C------- LEAKAGE POINTS INFORMATION -----------------------
  480 IF ( ICONF .NE. 1 ) GO TO 1200
      WRITE (NPMT6,998)
  998 FORMAT(/,' *** LEAKAGE POINTS ***'/)
      WRITE (NPMT6,997)
      I=1
      N=0
 4999 N=N+1
      IF(N.GT.NUMNP) GO TO 5999
      IF(KODE(N).NE.4) GO TO 4999
      IF(BL(N).GT.P(N)) GO TO 4999
      IF(Q(N).GE. 0.0 ) GO TO 4999
      NN(I)=N
      QN(I)=Q(N)
      PN(I)=P(N)
      BLN(I)=BL(N)
      GO TO 6999
 5999 NN(I)=0
      QN(I)=0.0
      PN(I)=0.0
      BLN(I)=0.0
 6999 I=I+1
      IF(I.LT.4) GO TO 4999
      WRITE (NPMT6,999)(NN(I),PN(I),BLN(I),QN(I),I=1,3)
      I=1
      IF(N.LT.NUMNP) GO TO 4999
  997 FORMAT (1H0,3('  NODE    HEAD        SURFACE        Q     ')/)
  999 FORMAT(2X,3(I5,1P3E12.3))
C ---- KODE=5 NODE  HEAD FIXED --> HEAD FREE & Q --------------
 1200 DO 600 I=1,NUMNP
      IF(KODE(I) .NE. 5) GO TO 600
      KODE(I)=0
      QO(I)=Q(I)
      WRITE (NPMT6,700) I,Q(I)
  700 FORMAT(15X,' AT NODE ',I5,5X,
     &         ' CHANGE KODE FROM 5 TO 0.         FIXED Q=  ',1PE13.3)
      P1(I)=P(I)
  600 CONTINUE
C---- SET UP NEW TIME STEP ------------------------------------------
  482 ICHE=ICHE+1
C     IF(ICHE .EQ. 1) GO TO 231
C---  TIME CHECK ----------------------------------------------
      IF(NT.LT.1) GO TO 602
      DO 543 I=1,NT
      IF(  TIME .EQ. TM(I) ) GO TO 602
  543 CONTINUE
C================================================ CALCULATION PATTERN
  602 IF( ISTEA .EQ. -1                   ) GO TO 560
      IF( ISTEA .EQ.  0 .AND.   LA .GE. 1 ) GO TO 481
      IF( ISTEA .EQ.  1                   ) GO TO 481
      IF( LA    .GE.  1                   ) GO TO 481
C--------------------------------------------------------- NEW B.C.--
  479 CALL NEWN(KODE,P1,P,QO,KM,QL,Q,NUMNP,NUMEL,
     &             DT,DTMAX,DMUL,TMAX,TOL,KB,BL,TM,NT,NT1,NCHA,NPRINT)
C
c      IF( NCHA .EQ. 0 ) GO TO 560
      IF( LA   .EQ. 0 ) GO TO 481
      GO TO 231
C-------------------------------------------------------------------
  481 IF( IFIN .GE. 1) GO TO 560
      LA=LA+1
      IF( LA .LE. 1 ) GO TO 555
      DT1=DT
      DO 822 N=1,NWEL+1
  822  AL1(N)=AL(N)
C ------------------------------------------------------------------
      IF ( DABS(TIME-TMAX) .LE. .001*DT
     & .AND. (ISTEA .EQ. 0 .OR. ISTEA .EQ. 1)) IFIN=1
      IF ( DABS(TIME-TMAX) .LE. .001*DT
     & .AND. (ISTEA .EQ. -2 .OR. ISTEA .EQ. 2)) IEND=1
      IF ( IFIN .GE. 1 ) GO TO 556
      IF ( IEND .GE. 1 ) GO TO 560 !479
      IF ( DT .LT. DTMAX) DT=DMUL*DT
      IF ( DT .GT. DTMAX) DT=DTMAX
C ---------------------------------------------------- TIME MODIFY --
      IF( NT.LE. 0 ) GO TO 550
      IF((TIME+DT) .GT. TM(NNT)
     &   .OR. DABS(TM(NNT)-TIME-DT) .LT. 0.2*DT) GO TO 551
      GO TO 550
  551 TIME1=TM(NNT)
      DT=TIME1-TIME
      IF (TIME1 .GT. TMAX ) GO TO 550
      GO TO 555
C ------------------------------------------------------------------
  550 IF ((TIME+DT).GT.TMAX.OR.DABS(TMAX-TIME-DT) .LT. .2*DT) DT=TMAX
     1 -TIME
  555 TIME=TIME+DT
  556 ITER=0
C----------------------------------------------------------
      LAA=0
C----------------------------------------------------------
      GO TO 230
  569 WRITE (NPMT6,699)
      WRITE (COMNT,699)
      CALL DISPLY(COMNT)
  699 FORMAT(' ???????  ERR IN STEADY STATE ??????')
      STOP
  560 CONTINUE
C----------------- ALL FINAL RESULTS WRITE TO FILE NO.9 --------------
      IF(IFILE .EQ. 0) GO TO 567
C      CALL OUTFD(X,Y,P,Q,KX,MK,NUMNP,NUMEL,LA,TIME,IRA)
C     DO 750 N=1,NUMNP
C     WRITE (NPMT9,760) N,KODE(N),KB(N),X(N),Y(N),P1(N),Q(N),
C    1                  QL(N),BL(N)
C 760 FORMAT(3I5,1P6E10.3)
C 750 CONTINUE
  567 CONTINUE
C
      COMNT='--- OUT FEM ---'
      CALL TIMCNT(COMNT)
      RETURN
      END
C----------
      SUBROUTINE WELBOR(P,QL,Q,QO,NWEL,DT,NUMNP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /WELL1/AL(21),AL1(21),AL4(21),NW(21)
      COMMON /WELL2/WARE(21),QOL(21),QL1(21),NWB(21,10)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION P(NUMNP),QO(NUMNP),QL(NUMNP),Q(NUMNP)
       ALF=0.8
      DO 10 I=1,NWEL
       N=NW(I)
       QQ=0.0
       QOO=0.0
       DO 20 J=1,N
        QOO=QOO+QO(NWB(I,J))
   20   QQ=Q(NWB(I,J))+QQ
       QMN=(QQ+QOL(I))*0.5
       AA=DT*(QOO+QMN)/WARE(I)
       AL(I)=ALF*AL(I)+(1.0-ALF)*(AL1(I)-AA)
       DO 30 J=1,N
   30  P(NWB(I,J))=AL(I)
        IF(AL(I).GE.QL(NWB(I,1)))GO TO 11
       DO 40 J=1,N
   40   P(N)=QL(NWB(I,1))
   11   QOL(I)=QQ
   10 CONTINUE
      RETURN
      END
C--------------
      SUBROUTINE MAFIL (KX,ALF,KM,NUMEL,IREL)
C************** KM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KX(NUMEL,4),ALF(NUMEL),KM(NUMEL)
      NUM=0
C      WRITE (NPMT6,10)
C   10 FORMAT (//'L H ELEMENT INFORMATION'/
C     1  '       ELEMENT    C O R N E R  N O D E S INFIL.FACTOR'/)
      DO 60 N=1,NUMEL
      IF(NUM-N) 20,60,40
   20 CONTINUE
      IF(IREL.EQ.0) THEN
         READ (NPMT5,30) NUM,(KX(NUM,I),I=1,4),KM(NUM),ALF(NUM)
      ELSE
         READ (NPMT5,*) NUM,(KX(NUM,I),I=1,4),KM(NUM),ALF(NUM)
      END IF
C
   30 FORMAT (6I5,E10.3)
      IF (KX(NUM,4) .EQ. 0) KX(NUM,4)=KX(NUM,3)
      IF (NUM .EQ. N ) GO TO 60
   40 DO 50 I=1,4
C       WRITE(*,*) 'N=',N
      KX(N,I)=KX(N-1,I)+1
   50 CONTINUE
      ALF(N)=ALF(N-1)
      KM(N)=KM(N-1)
   60 CONTINUE
C************* KM(N),I9
      DO 80 N=1,NUMEL
      WRITE (NPMT6,70)N,KX(N,1),KX(N,2),KX(N,3),KX(N,4),KM(N),ALF(N)
   70 FORMAT (I9,5I6,E10.3)
   80 CONTINUE
      RETURN
      END
C
      SUBROUTINE COLHT ( KX,IS,KM,RMK,NUMNP,NUMNP1,NUMEL,NWK )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL*4 HED,X,Y,P,VXX,VYY,C,PTIME
C
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KX(NUMEL,4),IS(NUMNP1),ISQ(300),LM(4),IST(3),KM(NUMEL),
     &          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( KM(J).LE.0) GO TO 35
CCC
       IF(RMK(1,J).GT.I .OR. RMK(2,J).LT.I ) GO TO 35
CCC
      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.300 ) 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
C
      NUM = NUMNP1
      NWK = IS(NUM) - 1
C
      RETURN
      END
C
      SUBROUTINE ADRSK ( KX,IS,ISP,KM,RMK,NUMNP,NUMNP1,NUMEL,NWK )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL*4 HED,X,Y,P,VXX,VYY,C,PTIME
C
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KX(NUMEL,4),IS(NUMNP1),ISP(NWK),LM(4),IST(3),KM(NUMEL),
     &          RMK(2,NUMEL)
C
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( KM(J).LE.0 ) GO TO 35
CCC
       IF(RMK(1,J).GT.I .OR. RMK(2,J).LT.I ) GO TO 35
CCC
      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
C
C     WRITE(*,200)
C     WRITE(NPMT6,200)
C 200 FORMAT(' *** IS **** & ***** ISP *****')
C     DO 210 I=1,NUMNP
C     I1 = IS(I)
C     I2 = IS(I+1) - 1
C     WRITE(*,220) I1,(ISP(J),J=I1,I2)
C     WRITE(NPMT6,220) I1,(ISP(J),J=I1,I2)
C 220 FORMAT(16I5)
C 210 CONTINUE
C
      RETURN
      END
C
C$LARGE
      SUBROUTINE MATIN (BK,TK,SK,NUMK,MK,NUMMAT,CC1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION BK(MK,NUMMAT),TK(MK,NUMMAT),SK(MK,NUMMAT),NUMK(NUMMAT)
      READ (NPMT5,32)(NUMK(J),J=1,NUMMAT-1)
   32 FORMAT(16I5)
      DO 120 M=1,NUMMAT-1
      NUMK(M)=NUMK(M)+1
      NUMKM=NUMK(M)
      READ (NPMT5,30)(BK(I,M),TK(I,M),I=1,NUMKM)
      READ (NPMT5,30)(BK(I,M),SK(I,M),I=1,NUMKM)
      BK(1,M)=0.0
      TK(1,M)=0.0
      SK(1,M)=0.0
   30 FORMAT(8E10.3)
C ************************************************************
C     CCMAX = -1.0
C     CCMIN = 1.0E20
C     DO 200 I=1,NUMMAT
C     DO 210 J=1,NUMKM
C     IF (TK(I,J) .GT. CCMAX) CCMAX=TK(I,J)
C     IF (TK(I,J) .LT. CCMIN) CCMIN=TK(I,J)
C 210 CONTINUE
C 200 CONTINUE
      CC1 = TK(1,1)
C ************************************************************
C      WRITE (*,130)M
      WRITE (NPMT6,130) M
  130 FORMAT(//' TABLE OF AQUIFER WIDTH  VERSUS PERMEABILITY 'I5/)
C      WRITE (*,90)(BK(I,M),TK(I,M),I=2,NUMKM)
      WRITE (NPMT6,90)(BK(I,M),TK(I,M),I=2,NUMKM)
      DO 40 I=1,NUMKM
   40 TK(I,M)=BK(I,M)*TK(I,M)
      DO 50 I=2,NUMKM
      TK(I,M)=TK(I-1,M)+TK(I,M)
   50 BK(I,M)=BK(I-1,M)+BK(I,M)
      WRITE (NPMT6,80)M
C       WRITE ( *,80)M
   80 FORMAT(//73H TABLE OF TOTAL HEAD VERSUS COEFFICIENT OF TRANSMISSIB
     1ILITY FOR MATERIAL I5/)
      WRITE (NPMT6,90)(BK(I,M),TK(I,M),I=1,NUMKM)
C      WRITE ( *,90)(BK(I,M),TK(I,M),I=1,NUMKM)
   90 FORMAT(1P10E12.4)
      WRITE (NPMT6,100)M
C      WRITE ( *,100)M
  100 FORMAT(//64H TABLE OF TOTAL HEAD VERSUS COEFFICIENT OF STORAGE FOR
     1 MATERIAL I5/)
      WRITE (NPMT6,90)(BK(I,M),SK(I,M),I=1,NUMKM)
C      WRITE ( *,90)(BK(I,M),SK(I,M),I=1,NUMKM)
  120 CONTINUE
      RETURN
      END
C-----------
      SUBROUTINE BOUND(HB,TIM,KBOUM,KBMK,ISW)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION HB(KBOUM,KBMK),TIM(KBOUM,KBMK)
      DO 10 M=1,KBMK
      READ (NPMT5,20)(HB(I,M),TIM(I,M),I=1,KBOUM)
   20 FORMAT(8E10.3)
C      DO 50 I=1,KBOUM
C      HB(I,M)=-1.0*HB(I,M)
C   50 CONTINUE
      IF(ISW. EQ. 1)  THEN
	  WRITE (NPMT6,30)M
	ELSE
	  WRITE (NPMT6,31)M
	END IF
C      WRITE ( *,30)M
   30 FORMAT(//" TABLE OF  HEAD   VS   TIME AT VARIABLE BOUNDARY ",
     1  I3/)
   31 FORMAT(//" TABLE OF  DISCHARGE   VS   TIME AT VARIABLE BOUNDARY ",
     1  I3/)
      WRITE (NPMT6,40)(HB(I,M),TIM(I,M),I=1,KBOUM)
C      WRITE ( *,40)(HB(I,M),TIM(I,M),I=1,KBOUM)
   40 FORMAT(1P10E12.4)
   10 CONTINUE
      RETURN
      END
C-----------
      SUBROUTINE EINFI (EI,ETIM,KBOUM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION EI(KBOUM),ETIM(KBOUM)
      READ (NPMT5,10)(EI(I),ETIM(I),I=1,KBOUM)
   10 FORMAT (8E10.3)
      WRITE (NPMT6,20)
C      WRITE ( *,20)
   20 FORMAT(// ' TABLE OF INFILTRATION RATE VERSUS TIME '/)
      WRITE (NPMT6,30)(EI(I),ETIM(I),I=1,KBOUM)
C      WRITE ( *,30)(EI(I),ETIM(I),I=1,KBOUM)
   30 FORMAT(10E12.4)
      RETURN
      END
C-------------
      SUBROUTINE NPIN (KODE,Q,P,P1,X,Y,KB,QL,KM,IC,KX,NUMNP,NUMEL,
     1                 NERR,BL,NXYP,IRNP)
C************* -KM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KODE(NUMNP),Q(NUMNP),P1(NUMNP),X(NUMNP),Y(NUMNP),
     1          P(NUMNP),QL(NUMNP),KB(NUMNP),BL(NUMNP),
     2          KM(NUMEL),IC(NUMNP),KX(NUMEL,4)
      NERR=0
      PAI=3.14152/180.
C      WRITE (NPMT6,11)
C   11 FORMAT (2H L H)
      WRITE (NPMT6,10)
C************ -7H MATEL.
   10 FORMAT (///' *** NODAL POINT INFORMATION ***',
     1//10H  NODE NO.,6X,4HKODE,3X,4X,6H NO.BU,7X,3H X ,12X,
     23H Y ,11X,5H HEAD,12X,1HQ,12X,2HQL,12X,2HBL/)
      NPR=0
      L=0
   20 L=L+1
C
C      WRITE (NPMT6,9000) L
C      WRITE (*,9000) L
C 9000 FORMAT ( ' ******* NODE NO COUNT =',I3)
C
      IF(IRNP.EQ.0) THEN
       READ (NPMT5,40)N,KODE(N),KB(N),X(N),Y(N),P1(N),Q(N),QL(N),BL(N)
      ELSE      
       READ (NPMT5,*)N,KODE(N),KB(N),X(N),Y(N),P1(N),Q(N),QL(N),BL(N)
      END IF
C
   40 FORMAT(3I5,6E10.3)
      IF(N-L) 50,90,70
   50 WRITE (NPMT6,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
   80 X(L)=X(L-1)+DX
      Y(L)=Y(L-1)+DY
      P1(L)=P1(L-1)
      KODE(L)=KODE(L-1)
      KB(L)=KB(L-1)
      Q(L)=Q(L-1)
      QL(L)=QL(L-1)
      BL(L)=BL(L-1)
      L=L+1
      IF(L .LT. N) GO TO 80
   90 NPR=N
      IF( L .LT. NUMNP ) GO TO 20
      CALL KODCHK(KX,KODE,KM,IC,NUMEL,NUMNP)
C
      DO 110 N=1,NUMNP
C     QL(N)=-1.0E20
C------------   -KM(N)
      IF( NXYP .LE. 0 ) GO TO 111
      Y(N)=PAI*Y(N)
      X1=X(N)*DCOS(Y(N))
      Y1=X(N)*DSIN(Y(N))
      X(N)=X1
      Y(N)=Y1
  111  WRITE (NPMT6,100)N,KODE(N),KB(N),X(N),Y(N),
     1               P1(N),Q(N),QL(N),BL(N)
  100 FORMAT (3I10,1P6E15.4)
      P(N)=P1(N)
  110 CONTINUE
C
      RETURN
      END
C---------------
      SUBROUTINE KODCHK(KX,KODE,KM,IC,NUMEL,NUMNP)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KX(NUMEL,4),KODE(NUMNP),IC(NUMNP),KM(NUMEL)
C
      DO 100 I=1,NUMNP
  100   IC(I)=0
C
      DO 10 I=1,NUMEL
       IF( KM(I).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(*,*)' NODE = ',KK
        IF( KODE(KK).GT. 0 .OR. IC(KK).GT.0 ) GO TO 11
      DO 20 J=1,NUMEL
       IF( KM(J).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  IC(KK)=1
   11 CONTINUE
   10 CONTINUE
C       DO 301 I=1,NUMNP
C       IF( KODE(I).GT. 0 ) WRITE (NPMT6,300)I,KODE(I)
C  301 CONTINUE
C  300 FORMAT(5I10)
      RETURN
      END
C$LARGE
C$STORAGE:2
C------------ AT2,NUMNPT,NWKT,INCOR
      SUBROUTINE RESET(KODE,A,B,Q,P,P1,D,CR,C,X,Y,KX,BK,TK,SK,NUMK,KM,
     1                 ALF,ISP,IS,AT2,
     2                 NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,LA,DT,DT1,ITER,
     3                 INTEG,EI1,IFIN,BL,QL,QO,TBEL,ICONF,
     4                 NWKT,INCOR,CLEAK,PUP,QQN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION A(NWK),B(NUMNP),Q(NUMNP),P(NUMNP),P1(NUMNP),
     1          D(NUMNP),CR(NUMNP),C(NUMNP),X(NUMNP),Y(NUMNP),
     2          KX(NUMEL,4),BK(MK,NUMMAT),TK(MK,NUMMAT),SK(MK,NUMMAT),
     3          KM(NUMEL),NUMK(NUMMAT),KODE(NUMNP),ALF(NUMEL),BL(NUMNP),
     4          QL(NUMNP),QO(NUMNP),ISP(NWK),IS(NUMNP1),
     5          AT2(NWKT),QQN(NUMNP)
C
C        WRITE(*,*)'/// IN RESET'
C
       IF(LA.GT.1 )DTT=.5*DT/DT1
C---- INITIALIZATION AND EXTRAPOLATION OF P VALUES
      IF(ITER .GT. 0) THEN
         DO 40 I=1,NUMNP
         B(I)=0.
         D(I)=0.
         A(I)=0.
         PP=.5*(P(I)+P1(I))
         IF (KODE(I) .LT. 1) P(I)=PP
         Q(I)=0.
         IF(P(I) .GE. QL(I) ) Q(I)=QO(I)+QQN(I)
         IF( KODE(I).GT.0) Q(I)=0.0
   40    CONTINUE
      ELSE
         DO 45 I=1,NUMNP
         B(I)=0.
         D(I)=0.
         A(I)=0.
         IF (LA.LE.1) PP=P1(I)
         IF(IFIN.GE.1) PP=P(I)
         IF(LA.GT.1. AND. IFIN .LT. 1) PP=P(I)+DTT*(P(I)-P1(I))
         P1(I)=P(I)
         IF(KODE(I) .LT. 1) P(I)=PP
         CONTINUE
         Q(I)=0.
         IF(P(I) .GE. QL(I) ) Q(I)=QO(I)+QQN(I)
         IF( KODE(I).GT.0) Q(I)=0.0
   45    CONTINUE
      ENDIF
      DO 30 J=NUMNP1,NWK
      A(J)=0.
   30 CONTINUE
      DO 65 I=1,NUMMAT
        DO 66 J=1,NUMNP
   66   CR(J)=-1.0D00
      DO 60 N=1,NUMEL
       IF( KM(N).EQ.I) THEN
      AIFL=ALF(N)
      NN=N
      CALL ELEM( KODE,A,P,P1,D,CR,C,X,Y,KX,BK,TK,SK,NUMK,KM,Q,ISP,IS,
     1           NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,NN,EI1,AIFL,BL,TBEL,
     2           ICONF,CLEAK,PUP )
      ENDIF
   60 CONTINUE
   65 CONTINUE
C------------ IF()
      IF (INCOR.GT.0) GO TO 11
        REWIND NPMT2
        WRITE (NPMT2) A
      GO TO 15
   11   DO 12 I=1,NWK
        AT2(I)=A(I)
   12   CONTINUE
C---- COMPLETE COSTRUCTION OF RHS VECTOR AND FORM EFFECTIVE MATRIX
   15 STE=1.0
      IF(LA.LT.1.OR.IFIN.GE.1) STE=0.0
      PP=2.
      IF( LA .LE. 1 .OR. INTEG .EQ. 1 .OR. IFIN .GE. 1) PP=1.
      DO 110 I=1,NUMNP
      I1 = IS(I)
      B(I)=PP*(Q(I)-B(I)+D(I)*P1(I)/DT*STE)
      IF( LA .LE. 1 .OR. INTEG .EQ. 1 .OR. IFIN .GE. 1) GO TO 100
        I2 = IS(I+1) - 1
        DO 90 J=I1,I2
        K = ISP(J)
        B(I) = B(I) - A(J)*P1(K)
   90   CONTINUE
  100 A(I1)=A(I1)+PP*D(I)/DT*STE
  110 CONTINUE
C        WRITE(*,*)'/// OUT RESET'
      RETURN
      END
C-------------
      SUBROUTINE ELEM ( KODE,A,P,P1,D,CR,C,X,Y,KX,BK,TK,SK,NUMK,KM,Q,
     1                  ISP,IS,NWK,NUMNP,NUMNP1,NUMEL,NUMMAT,MK,N,EI1,
     2                  AIFL,BLE,TBEL,ICONF,CLEAK,PUP )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION E(3,3),S(4,4),A(NWK),P(NUMNP),D(NUMNP),
     1          CR(NUMNP),C(NUMNP),X(NUMNP),Y(NUMNP),KX(NUMEL,4),
     2          P1(NUMNP),KODE(NUMNP),BK(MK,NUMMAT),TK(MK,NUMMAT),
     3          SK(MK,NUMMAT),NUMK(NUMMAT),KM(NUMEL),LM(4),IX(3),
     4          Q(NUMNP),BLE(NUMNP),ISP(NWK),IS(NUMNP1)
C
C        WRITE(*,*)'/// IN ELEM'
C
      DO 30 I=1,4
      LM(I)=KX(N,I)
   30 CONTINUE
      NUS=4
      IF(LM(3) .EQ. LM(4)) NUS=3
       CALL ZERO(S,16)
C---- DETERMINE NONLINEAR VARIABLES FOR CORNER NODES
      M=KM(N)
       NUMKM=NUMK(M)
      PAV=0.0
	BLEI=0.0
      DO 60 K=1,NUS
	  PCAL=P(LM(K))
       IF(KODE(LM(K)).GE.1) PCAL=0.5*(PCAL+P1(LM(K)))
        PAV=PAV+P(LM(K))/DFLOAT(NUS)
      I=LM(K)
C       BLEI=BLE(I)
       BLEI=BLEI+ BLE(I)/DFLOAT(NUS)
   60 CONTINUE
      CM=PAV
C      IF( CR(I) .GT. 0.0) GO TO 60
C       CM=P(I)

       IF( ICONF.GE.1 .AND. CM.GT.BLEI)CM=BLEI
C------------------------------------------------------------
       CRI=TK(NUMKM,M)
       CI=SK(NUMKM,M)
        IF(CM .LT. BLEI) THEN
          CM=CM-(BLEI-BK(NUMKM,M))
         IF(CM .GT.0.0) THEN
          CALL INTERP (BK,TK,CM,CRI,NUMKM,NUMMAT,M,MK)
          CALL INTER2 (BK,SK,CM,CI ,NUMKM,NUMMAT,M,MK)
C          IF( CRI .GT.TK(NUMKM,M)) CRI=TK(NUMKM,M)
         ELSE
          CRI=1.0D-15
          CI=0.0D00
         ENDIF
C------------------------------------------------------------
        ENDIF
	DO K=1,NUS
	 I=LM(K) 
       CR(I)=CRI
       C(I)=CI
	END DO
      CONTINUE
C      IF( PAV. LT. BLEI ) PAV=BLEI
      QSTORE=.0
C---- LOOP ON SUBELEMENTS
      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=DABS(CK*BJ-CJ*BL)
      DEL2=(CK*BJ-CJ*BL)
Cccc	 IF( DEL2.LE.0.0) THEN
cccc	   write(NPMT6,*) "A .LE. 0 ", N,I,J,L
cccc 	  STOP 1234
cccc       END IF
C
      DELE= DEL2*0.5
      EI2=EI1*DELE*AIFL
      EI2=EI2+(PUP-PAV)*CLEAK*DELE
      EI3=EI2/3.0D00
C
      CM=(CR(I)+CR(J)+CR(L))/3.0D00
      COMM=.5*CM/DEL2
      QSTORE=.5*DEL2
      E(1,1)=BI*BI+CI*CI
      E(1,2)=BI*BJ+CI*CJ
      E(1,3)=BI*BL+CI*CK
      E(2,1)=E(1,2)
      E(2,2)=BJ*BJ+CJ*CJ
      E(2,3)=BJ*BL+CJ*CK
      E(3,1)=E(1,3)
      E(3,2)=E(2,3)
      E(3,3)=BL*BL+CK*CK
      IX(1)=1
      IX(2)=K+1
      IX(3)=K+2
      DO 110 I=1,3
      II=IX(I)
      IJ=LM(II)
C
C----NINFILTRATION FOR EACH NODE
      Q(IJ)=Q(IJ)+EI3
      AA=C(IJ)
      DO 100 J=1,3
      JJ=IX(J)
      JI=LM(JJ)
      AA=AA+C(JI)
      S(II,JJ)=S(II,JJ)+E(I,J)*COMM
  100 CONTINUE
C----------------------------------( 0.0833=1/12. )
      D(IJ)=D(IJ)+AA*QSTORE/12.0D00
  110 CONTINUE
  120 CONTINUE
C---- ADD ELEMENT COTRIBUTION TO A MATRIX
      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
C        WRITE(*,*)'/// OUT ELEM'
      RETURN
      END
C-------------
      SUBROUTINE INTERP (X,Y,XX,YY,N,M,J,K)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION X(K,M),Y(K,M)
C---- LINEAR INTERPOLATION
      DO 10 I=2,N
      IF (XX .GT. X(I,J) .AND. I .LT. N) GO TO 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
      GO TO 20
   10 CONTINUE
   20 RETURN
      END
C---------
      SUBROUTINE INTER2 (X,Y,XX,YY,N,M,J,K)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION X(K,M),Y(K,M)
C------LINEAR INTERPOLATION
      DO 10 I=2,N
      IF(XX .GE.X(I,J) .AND. I .LT. N) GO TO 10
      YY=Y(I,J)
      GO TO 20
   10 CONTINUE
   20 RETURN
      END
C------------
      SUBROUTINE INTER3(X,Y,XX,YY,K)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION X(K),Y(K)
C---- LINEAR INTERPOLATION
      DO 10 I=2,K
      IF(XX .GT. X(I) .AND. I .LT. K ) 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 INFIL(X,Y,XX,YY,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION X(N),Y(N)
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$LARGE
      SUBROUTINE SOLVCG (A,B,P,ISP,IS,R,PP,KODE,KODET,C,NUMNP,NUMNP1,
     &                   NWK,TSOLV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
      CHARACTER*80 COMNT
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION A(NWK),B(NUMNP),P(NUMNP),ISP(NWK),IS(NUMNP1),R(NUMNP),
     1          PP(NUMNP),KODE(NUMNP),KODET(NUMNP),C(NUMNP)
C        WRITE(*,*)'/// IN SOLV'
      ITMAX = NUMNP*5
C ---- ATOL -----
      ATOL = 1.0D-28
      ATOL2 = ATOL*ATOL
      NUMKD=0
C--       TTA=0.0
C--       TAMAX=0.0
C******TA--SUM(AIJ),TAA--SUM(TA)
      DO 5 I=1,NUMNP
      R(I) = 0.0
      PP(I) = 0.0
      IF( KODE(I) .GE. 1 ) GO TO 5
      NUMKD=NUMKD+1
      KODET(NUMKD)=I
      II = IS(I)
      II1 = IS(I+1) - 1
      S = 0.0
C-----
       IF( DABS(A(II)) .LE. ATOL ) A(II)=0.0D00
       N = ISP(II)
       AJ=A(II)
       S = S + AJ*P(N)
       TAJ=0.0
       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)
C-----
C--       TAJ=TAJ+DABS(AJ)
   10 CONTINUE
C--       TAA=TAJ/DABS(A(II))
      R(I) = B(I) - S
      C(I)=A(II)
      PP(I) = R(I)/C(I)
C--       TTA=TTA+TAA
C         IF (TAA.LT.1.0) GOTO 9
C         WRITE (NPMT6,662) TAA,I
C         WRITE(COMNT,662) TAA,I
C         CALL DISPLY(COMNT)
C  662    FORMAT('--SUM(AIJ/AII)/NUMKD=',D15.7,
C     &          '   I='I5)
C--    9  IF (TAMAX.GT.TAA) GOTO 5
C--       TAMAX=TAA
C--       NPNO=I
    5 CONTINUE
       IF(NUMKD .LE. 0) GOTO 82
C--       QQ=NUMKD
C--       TTANM=TTA/QQ
C
C--      WRITE (NPMT6,661) TTANM,TAMAX,NPNO
C--      WRITE(COMNT,661) TTANM,TAMAX,NPNO
C--      CALL DISPLY(COMNT)
C        WRITE(*,*)'/// IN SOLV1'
C
      IT = 1
   30 PR = 0.0
      PAP = 0.0
      DO 40 III=1,NUMKD
      I=KODET(III)
C*******************      IF( KODE(I) .GE. 1 ) GO TO 40
      II = IS(I)
      II1 = IS(I+1) - 1
      PR = PR + PP(I)*R(I)
      BI=0.0
      DO 50 J=II,II1
C******************      IF( DABS(A(J)) .LE. ATOL ) GO TO 50
      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.0
      EPSX = 0.0
      DO 60 III=1,NUMKD
      I=KODET(III)
C******************      IF( KODE(I) .GE. 1 ) GO TO 60
      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
C        WRITE(*,*)'/// IN SOLV3'
      EPSX = DSQRT(EPSX)
C         WRITE (NPMT6,665) IT,EPSX
C         WRITE(COMNT,665) IT,EPSX
C         CALL DISPLY(COMNT)
C  665    FORMAT(' ***-------- IT =',I5,
C     1          '    EPSX=',D15.7)
      IF( EPSX - TSOLV ) 80,80,63
   63 IF( IT - ITMAX ) 65,75,75
   65 E = -RAP/PAP
C
      DO 70 II=1,NUMKD
       I=KODET(II)
      PP(I) = R(I)/C(I) + E*PP(I)
   70 CONTINUE
      IT = IT + 1
C        WRITE(*,*)'IT ',IT,EPSX
      GO TO 30
C
   75 WRITE (NPMT6,600) ITMAX
      WRITE(COMNT,600) ITMAX
      CALL DISPLY(COMNT)
  600 FORMAT(' *** FAILURE TO CONVERGENCE IN SOLVER -- ITMAX =',I5)
      STOP 999
   80 WRITE (NPMT6,660) IT,EPSX
      WRITE(COMNT,660) IT,EPSX
      CALL DISPLY(COMNT)
  660 FORMAT(' *** SUCCESS TO CONVERGENCE IN SOLVER -- IT =',I5,
     1       '    EPSX=',D15.7)
	 GOTO 81
   82	 WRITE(NPMT6,670)
    	  WRITE(*,670)
  670  FORMAT('--- ALL VARIABLE(HEAD) ARE FIXED.
     1		  SO,NO CALCULATED IN PCG ')
C
C--  661 FORMAT(' -- SUM(AIJ/AII)/NUMKD=',D15.7,
C--     1       '  MAX(AIJ/AII)=',D15.7,
C--     2       ' NODE=',I5)
C
C        WRITE(*,*)'/// OUT SOLV'      
   81	RETURN
      END
C-------
C-------------- AT2,NUMNPT,NWKT,INCOR
      SUBROUTINE FIXQ ( A,B,KODE,Q,P,P1,ISP,IS,AT2,
     1                  NUMNP,NUMNP1,NWK,LA,INTEG,IFIN,
     2                  NUMNPT,NWKT,INCOR )
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
C
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION A(NWK),B(NUMNP),KODE(NUMNP),Q(NUMNP),P(NUMNP),
     1          P1(NUMNP),ISP(NWK),IS(NUMNP1),AT2(NWKT)
C
C------------ IF()
      IF (INCOR.GT.0) GO TO 11
      REWIND NPMT2
      READ (NPMT2) A
      GO TO 15
   11 DO 13 I=1,NWK
      A(I)=AT2(I)
   13 CONTINUE
C
   15 PPP=1.0D00
      PP=2.0D00
      IF( ( LA   .GT. 1 .AND. INTEG .NE. 1 ) .OR.
     1    ( IFIN .LT. 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)=0.0
      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)/PP
   40 CONTINUE
C
      RETURN
      END
C$LARGE
      SUBROUTINE NEWN(KODE,P1,P,QO,KM,QL,Q,NUMNP,NUMEL,
     &             DT,DTMAX,DMUL,TMAX,TOL,KB,BL,TM,NT,NT1,NCHA,NPRINT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION KODE(NUMNP),P1(NUMNP),P(NUMNP),QO(NUMNP),QL(NUMNP),
     1          KM(NUMEL),Q(NUMNP),KB(NUMNP),BL(NUMNP),TM(NT1)
C
      READ (NPMT5,20)NCHA
   20 FORMAT(16I5)
      IF ( NCHA .LE. 0 ) GO TO 60
      WRITE (NPMT6,140)
C     WRITE (*,140)
  140 FORMAT ( ////,'  *****  NEW  BOUNDARY CONDITION  ******'//)
C------------------------------- READ NEW TIME STEP DATA INFORMATION
      READ (NPMT5,90)DT,DTMAX,DMUL,TMAX,TOL,NT,NPRINT
      IF ( NT .LE. 0 ) GO TO 94
      READ (NPMT5,95)(TM(I),I=1,NT)
   90 FORMAT (5F10.0,2I5)
   95 FORMAT (5F10.0)
   94 WRITE (NPMT6,143)DT,DTMAX,DMUL,TMAX,NT,NPRINT,TOL
  143 FORMAT(/' * NEW **  INITIAL TIME----------- '1PE10.3/
     1        ' * NEW **  MAXIMUM TIME----------- 'E10.3/
     2        ' * NEW **  MULTIPLY FACTOR-------- 'E10.3/
     3        ' * NEW **  FINAL TIME------------- 'E10.3/
     4        ' * NEW **  PRINT OUT STEP-(NT)---- 'I5/
     5        ' * NEW **  PRINT STEP-(NPRINT)---- 'I5/
     6        ' * NEW **  MINIMUM ERROR---------- 'E10.3/)
C
      WRITE (NPMT6,91)
   91 FORMAT (/' *** NEW PRINT OUT TIME INFORMATION *** ')
      WRITE (NPMT6,92)(TM(I),I=1,NT)
   92 FORMAT (1P10E10.3)
C
      WRITE (NPMT6,10)
   10 FORMAT(//10H  NODE NO.,6X,4HKODE,8X,2HKM,5X,2HKB,5X,5H HEAD,14X,
     &1HQ,13X,2HQL,13X,2HBL/)
      DO 50 I=1,NCHA
      READ (NPMT5,30)N,KODE(N),KB(N),P1(N),QO(N),QL(N),BL(N)
      WRITE (NPMT6,40)N,KODE(N),KB(N),P1(N),QO(N),QL(N),BL(N)
      P(N)=P1(N)
   30 FORMAT(3I5,4F10.0)
   40 FORMAT(3I10,1P4E15.5)
   50 CONTINUE
   60 RETURN
      END
C----------
      SUBROUTINE PRINTO (Q,P,NUMNP,TBEL,KODE,NPNP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      COMMON /NPNODE/NPNOD(200)
C
      DIMENSION Q(NUMNP),P(NUMNP),NN(4),PP(4),QQ(4),KODE(NUMNP)
      WRITE (NPMT6,10)
   10 FORMAT(//1H0,4(32H NODE      HEAD         Q       )/)
C
      QM=0.0
      QP=0.0
      PMAX=0.0
      NOD=0
C
      J=0
      N=0
   20 N=N+1
      J=J+1
      IF(N .GT. NUMNP) GO TO 30
      NN(J)=N
      PP(J)=P(N)-TBEL
      QQ(J)=Q(N)
C---- MIYO
      IF(KODE(N).EQ.1) GOTO 2000
      IF(P(N).GE.PMAX) THEN
      PMAX=P(N)
      NOD=N
      ENDIF
C
 2000 IF( Q(N) .GT. 0.0 ) QP=QP+Q(N)
      IF( Q(N) .LT. 0.0 ) QM=QM+Q(N)
C----
      GO TO 40
   30 NN(J)=0.
      PP(J)=0.
      QQ(J)=0.
   40 IF(J .LT. 4) GO TO 20
      WRITE (NPMT6,50)(NN(J),PP(J),QQ(J),J=1,4)
   50 FORMAT (1X,4(I6,1P2E13.5))
      J=0
      IF( N .LT. NUMNP) GO TO 20
C------------- MIYO
      WRITE (NPMT6,2100) QP,QM,PMAX,NOD
 2100 FORMAT(//'***   QP( +  ) =',1PE13.5/
     1         '***   QM( -  ) =',1PE13.5/
     2         '---   Hedmax =  ',1PE13.5,
     3         '        NODE =  ',I5    )
C
      IF (NPNP.LE.0) GOTO 70
      WRITE (NPMT6,60)
   60 FORMAT(//1H0,(30H NODE    HEAD         Q       )/)
      DO 3000 I=1,NPNP
      PPP=P(NPNOD(I))-TBEL
      WRITE (NPMT6,3100) NPNOD(I),PPP,Q(NPNOD(I))
 3100 FORMAT (1X,I4,1P2E13.5)
 3000 CONTINUE
   70 RETURN
      END
C---------------
      SUBROUTINE PRINT1 (DEF,NUMNP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION DEF(NUMNP),NN(10),PP(10)
      WRITE (NPMT6,11)
   11 FORMAT (/' **** DRAWDOWN FORM INITIAL STEADY STATE ****
     &(+) --> UP     (-) --> DOWN '/)
      WRITE (NPMT6,10)
   10 FORMAT(/1H0,10('NODE---DEF.  ')/)
      J=0
      N=0
   20 N=N+1
      J=J+1
      IF(N .GT. NUMNP) GO TO 30
      NN(J)=N
      PP(J)=DEF(N)
      GO TO 40
   30 NN(J)=0.
      PP(J)=0.
   40 IF(J .LT. 10) GO TO 20
      WRITE (NPMT6,50)(NN(J),PP(J),J=1,10)
C      WRITE ( *,50)(NN(J),PP(J),J=1,10)
   50 FORMAT (1X,10(I3,1X,F7.2,2X))
      J=0
      IF( N .LT. NUMNP) GO TO 20
      RETURN
      END
C---------------
      SUBROUTINE VECTOR ( P,X,Y,KX,BK,TK,NUMK,KM,NUMNP,NUMEL,
     1           NUMMAT,MK,BLE,VXX,VYY,TBEL,CMPP,KFLAG )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      IMPLICIT INTEGER*2 (I-N)
C      REAL HED,X,Y,P,VXX,VYY,C,PTIME
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      COMMON /ADDI1/ IPPSW1
      DIMENSION X(NUMNP),Y(NUMNP),KX(NUMEL,4),P(NUMNP),
     1 BK(MK,NUMMAT),TK(MK,NUMMAT),
     2 NUMK(NUMMAT),KM(NUMEL),BLE(NUMNP),
     3 VXX(NUMEL),VYY(NUMEL),CMPP(NUMEL),
     4 YY(6),CR(4),NNN(2),VX1(2),VY1(2),V(2),DEG(2)
      AMAX=0.0
      AMIN=1.0D20
      AMAXDD=0.0
      AMINDD=1.0D20
      DO 800 N=1,NUMEL
       VXX(N)=0.0
       VYY(N)=0.0
         IF( KM(N).LE. 0 ) GO TO 800
      NUS = 4
      IF (KX(N,3) .EQ. KX(N,4)) NUS = 3
      XM = 0.0
      YM = 0.0
      PM = 0.0
      AB = 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
      AB = AB + BLE(I)/ANUS
  100 CONTINUE
      M=KM(N)
      CM=PM-(AB-BK(NUMK(M),M))
      CMPP(N)= CM
	IF(CM .LT. 0.0) THEN
	CRK=0.0
	CMPP(N) = 1.0E-20
	ELSEIF(CM .GT. BK(NUMK(M),M) ) THEN
	CRK=TK(NUMK(M),M)
	CMPP(N) = BK(NUMK(M),M)   !地表面オーバー
	ELSE
      CALL INTERP (BK,TK,CM,CRK,NUMK(M),NUMMAT,M,MK)
      IF ( CRK .LE. 0.0 ) CRK=0.0
      IF ( CRK .GT. TK(NUMK(M),M) ) CRK=TK(NUMK(M),M)
	END IF
C
      CMPPD = CMPP(N) 
C
c      IF(M .GE.  28095 ) THEN
c      write(87,*)  N,CRK,PM,CM ,ab,BK(NUMK(M),M),NUMK(M),M
c	END IF
C      DO 500 K=1,NUS
C      CR(K) = CRK
C  500 CONTINUE
C      PRM = 0.0
C      DO 600 K =1,NUS
C  600 PRM = PRM + CR(K)
      PRM = CRK
      PRM1 = PRM
      PRM2 = PRM
      Y1K = XM
      Y2K = YM
      Q1 = 0.0
      Q2 = 0.0
      D1 = 0.0
      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
      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)
  700 CONTINUE
      Q1 = Q1/D1
      Q2 = Q2/D1
      Q3=Q1*Q1+Q2*Q2
      Q3=DSQRT(Q3)
C
      IF( Q3 .LT. AMAX ) GO TO 750
      AMAX = Q3
      NMAX = N
  750 IF( Q3 .GT. AMIN ) GO TO 760
      AMIN = Q3
      NMIN = N
  760 VXX(N) = Q1
      VYY(N) = Q2
C
      Q1D = Q1/D1/CMPPD
      Q2D = Q2/D1/CMPPD
      Q3D=Q1*Q1+Q2*Q2
      Q3D=DSQRT(Q3)
C
      IF( Q3D .LT. AMAXDD ) GO TO 751
      AMAXDD = Q3D
      NMAXDD = N
  751 IF( Q3D .GT. AMINDD ) GO TO 761
      AMINDD = Q3D
      NMINDD = N
  761 CONTINUE
C
  800 CONTINUE
C ************* VEROCITY PRINT OUT FILE *******************
      IF(IPPSW1 .EQ. 0) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF (KFLAG .LT. 1) GO TO 111
      WRITE (NPMT6,660)
  660 FORMAT(//,2('  ELEM      V(X)       V(Y)        VV      DEGREE
     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)=SQRT(VV)
C************** VELOCITY VECTOR DERECTION CHECK ************
      XVECT=DABS(VXX(N))
      YVECT=DABS(VYY(N))
      IF(XVECT .GT. 1.0E-15) GO TO 29
        IF(VYY(N) .GT. 0) GO TO 101
           DEG(J)=-90.0
           GO TO 30
  101      DEG(J)=90.0
           GO TO 30
   29 IF(YVECT .GT. 1.0E-15) GO TO 35
        IF(VXX(N) .GT.0) GO TO 102
           DEG(J)=180
           GO TO 30
  102      DEG(J)=0.0
           GO TO 30
   35 GRAD=VYY(N)/VXX(N)
      DE=DATAN(GRAD)
      DEG(J)=DE*180.0/3.141592654D00
      IF( VX1(J) .LT. 0 .AND. VY1(J) .LT. 0) DEG(J)=DEG(J)-180.0
      IF( VX1(J) .LT. 0 .AND. VY1(J) .GT. 0) DEG(J)=DEG(J)+180.0
      GO TO 30
   20 NNN(J)=0
      VX1(J)=0.
      VY1(J)=0.
      V(J)=0.
      DEG(J)=0.
   30 IF(J.LT.2) GO TO 10
      WRITE (NPMT6,60)(NNN(J),VX1(J),VY1(J),V(J),DEG(J),J=1,2)
   60 FORMAT(1X,2(I6,2X,1P3E11.3,0PF8.1,5X))
      J=0
      IF(N.LT.NUMEL) GO TO 10
C*****************************************************************
  111 WRITE (NPMT6,66)AMAX,NMAX,AMIN,NMIN
   66 FORMAT(' '/,
     1         ' MAX. VELOCITY =',E10.3,5X,'AT FLEMENT NO.',I5/,
     2         ' MIN. VELOCITY =',E10.3,5X,'AT FLEMENT NO.',I5/)
C*****************************************************************
      ELSE   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF (KFLAG .LT. 1) GO TO 1111
      WRITE (NPMT6,1660)
 1660 FORMAT(//,2('  ELEM      V(X)       V(Y)      V(X)/D    V(Y)/D
     1    ')/)
      J=0
      N=0
 1010 N=N+1
      J=J+1
      IF(N .GT. NUMEL) GO TO 1020
      NNN(J)=N
      VX1(J)=VXX(N)
      VY1(J)=VYY(N)
      DEG(J) = CMPP(N) 
      GO TO 1030
 1020 NNN(J)=0
      VX1(J)=0.
      VY1(J)=0.
      V(J)=0.
      DEG(J)=0.
 1030 IF(J.LT.2) GO TO 1010
      WRITE (NPMT6,1060)(NNN(J),VX1(J),VY1(J),
     1                  VX1(J)/DEG(J),VY1(J)/DEG(J),J=1,2)
 1060 FORMAT(1X,2(I6,2X,1P4E11.3,2X))
      J=0
      IF(N.LT.NUMEL) GO TO 1010
C*****************************************************************
 1111 WRITE (NPMT6,66) AMAX,NMAX,AMIN,NMIN
      WRITE (NPMT6,1066)AMAXDD,NMAXDD,AMINDD,NMINDD
 1066 FORMAT(' '/,
     1         ' MAX. VELOCITY/D =',E10.3,5X,'AT FLEMENT NO.',I5/,
     2         ' MIN. VELOCITY/D =',E10.3,5X,'AT FLEMENT NO.',I5/)
C*****************************************************************
      END IF   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      RETURN
      END
C---------------
      SUBROUTINE ZERO(A,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      DIMENSION A(N)
      DO 10 I=1,N
   10 A(I)=0.0D0
      RETURN
      END
C----------------
C------------------------------- ( SET UP )
      SUBROUTINE SETFL(MTOTBD,MTOTIB)
C-------------
      CHARACTER*20 FNAME2,FNAME3,FNAME4,FNAME5,FNAME6
      CHARACTER*20 FNAM10,FNAM11
      CHARACTER*20 FNAME9,DDDD
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
C ================================================== READ OPEN FILE NAME
C ------------------------------- READ CONTOROL DATA FROM NO.5 --------
      WRITE(*,100)
      READ (*,20) FNAME5
C      fname5='data.dat'
C-------------------------------- ALL RESULTS WRITE TO FILE NO.6 -------
      WRITE (*,101)
      READ (*,20) FNAME6
C      fname6='data.prn'
C--------------------------------------------------- READ OPEN FILE NAME
      OPEN (NPMT98,FILE='WORK.LST',STATUS='OLD',FORM='FORMATTED')
      READ (NPMT98,20) FNAME2
      READ (NPMT98,20) FNAME3
      READ (NPMT98,20) FNAME4
      READ (NPMT98,20) FNAME9
C
      MTOTBD2  = 0
      MTOTIB2  = 0
C
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
      READ (NPMT98,20,END=999) DDDD
C
      READ (NPMT98,*,END=999) MTOTBD2 
      READ (NPMT98,*,END=999) MTOTIB2 
  999 CONTINUE
	IF(MTOTBD2  .LE. MTOTBD)  MTOTBD2 =MTOTBD 
	IF(MTOTIB2  .LE. MTOTIB)  MTOTIB2 =MTOTIB 
	 MTOTBD =MTOTBD2 
	 MTOTIB =MTOTIB2  
  997 CONTINUE  
C      READ (NPMT98,20) FNAM10
C      READ (NPMT98,20) FNAM11
   20 FORMAT (A)
C=================================================== FILE OPEN =========
      OPEN (NPMT2, FILE=FNAME2, STATUS='UNKNOWN',FORM='UNFORMATTED')
      OPEN (NPMT3, FILE=FNAME3, STATUS='UNKNOWN',FORM='UNFORMATTED')
      OPEN (NPMT4, FILE=FNAME4, STATUS='UNKNOWN',FORM='UNFORMATTED')
      OPEN (NPMT5, FILE=FNAME5, STATUS='OLD',FORM='FORMATTED')
      OPEN (NPMT6, FILE=FNAME6, STATUS='UNKNOWN',FORM='FORMATTED')
      OPEN (NPMT9, FILE=FNAME9, STATUS='UNKNOWN',FORM='UNFORMATTED')
C      OPEN (NPMT10,FILE=FNAM10, STATUS='UNKNOWN',FORM='UNFORMATTED')
C      OPEN (NPMT11,FILE=FNAM11, STATUS='UNKNOWN',FORM='UNFORMATTED')
C
C      OPEN (15, FILE='DIAD.PRN',STATUS='UNKNOWN',FORM='FORMATTED')
C      OPEN (16, FILE='PLOTT.OUT',STATUS='UNKNOWN',FORM='FORMATTED')
C---- READ MAIN DATA --------------------------------------------------
  100 FORMAT(//' KEY IN INPUT DATA FILE NAME [EX. B:DATA1.DAT] ====> ')
  101 FORMAT(/ ' KEY IN OUTPUT FILE NAME     [EX. B:DATA1.PRN] ====> ')
C-------
      NZZZ=0
      RETURN
      END
C
      SUBROUTINE TIMCNT(NNNN)
      CHARACTER*40 NNNN
      CHARACTER*8 TIM
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      NZZZ=0
c     CALL GETTIM(I1,I2,I3,I4)
C      WRITE(*,10) NNNN,TIM
C      WRITE(NPMT6,10) NNNN,TIM
C      CALL ICLOCK(ITIM)
C       WRITE(*,12) nnnn,I1,I2,I3,I4
C   12  FORMAT(a20,3(I2,':'),'.',I2)
C      WRITE(*,11) NNNN,ITIM
C      WRITE(NPMT6,11) NNNN,ITIM
   10 FORMAT(/A40,2X,A8)
   11 FORMAT(/A40,2X,I10)
      RETURN
      END
C
      SUBROUTINE DISPLY(COMNT)
      CHARACTER*80 COMNT
      COMMON /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      COMMON /NPMTB/NPMT10,NPMT11,NPMT98
      NZZZ=0
      WRITE(*,10) COMNT
   10 FORMAT(/A80/)
      RETURN
      END
      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 /NPMTA/NPMT5,NPMT6,NPMT2,NPMT3,NPMT4,NPMT9
      DIMENSION X(NUMNP),Y(NUMNP),KX(NUMEL,4),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
C
      DO K=1,NUS
      IF(LM(K) .GT. NUMNP) THEN
	   IERSW = 1
	   write(NPMT6,*) "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(NPMT6,*) "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(NPMT6,*) "ERROR ELEMENT  AREA IS NEGATIVE " , N
	   IERSW = 1
       END IF
C
  120 CONTINUE
  130 CONTINUE
C      IF(IERSW .EQ. 1 ) STOP
C        WRITE(*,*)'/// OUT ELEM'
      RETURN
      END