C$DEBUG
C     TRANS
C       AC-GWAP 
C
C
      USE DFLIB
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*4 HED(20)
      CHARACTER*20 FNAM10,FNAME6,FNAMEP
      COMMON  / FEMAP / ISETID,IU80 ,IPSTG
	COMMON /AVSCC / ISTPNUM 
      DIMENSION X(444444),Y(444444),H(444444),KX(444444,5),
     &          VXX(444444),VYY(444444),KD(444444),QQ(444444),
     &          VX(444444),VY(444444),IC(444444),QL(444444),BL(444444),
     &          ALF(444444),
     &          KB(444444),DDD(444444),VXD(444444),VYD(444444)
C
      DIMENSION XPNT(22222),YPNT(22222),KELM(4),
     1        WGG(22222,3),IWGT(22222,3)

C	result=SETEXITQQ(QWIN$EXITNOPERSIST)
      NUMDAT=444444
      NUMDAT2=22222
      IPSTG = 20000
      ISTPNUM=0 
      ISETID= 0
	IU80 = 8
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,1)
    1 FORMAT('**********************************************'/
     &       '* This program is for AC-GWAP .              *'/
     &       '*  It uses the 4th file inWORK.LST.          *'/
     &       '*                        by Y.Shinshi,T.IMAI *'/
     &       '*   NUMNP= upto 444444                       *'/
     &       '*   NUMEL= upto 444444                       *'/
     &       '*--------------------------------------------*'/
     &       '*                               Yes/0  No/1  *'/
     &       '**********************************************')
 2    WRITE(*,*)'  INPUT No. ? ='
      READ(*,*) ICH
      IF(ICH.LT.0.OR.ICH.GT.1) GOTO 2 
      WRITE(*,*) 
      WRITE(*,*) 
      WRITE(*,3) 
    3 FORMAT('    +--------------------------------------------+'/
     &       '    | You can choose the file to transform.      |'/
     &       '    |   Which is your choice?                    |'/
     &       '    |                                            |'/
     &       '    |   0:  the 4th file in WORK.LST             |'/
     &       '    |   1:  input the file-name.                 |'/
     &       '    +--------------------------------------------+')
      WRITE(*,*) 
    4 WRITE(*,*)'  INPUT No. ? ='
      READ(*,*) ICH2
       IF (ICH2.EQ.0) THEN
        OPEN (98,FILE='WORK.LST', STATUS='OLD',FORM='FORMATTED')
        READ(98,1000)FNAM10
        READ(98,1000)FNAM10
        READ(98,1000)FNAM10
        READ(98,1000)FNAM10
       ENDIF
       IF(ICH2.EQ.1) THEN
        WRITE(*,*)'  INPUT the File Name. '
        READ (*,1000) FNAM10
       ENDIF
       IF(ICH2.LT.0.OR.ICH2.GT.1) GOTO 4 
 1000 FORMAT(A)        
      OPEN (10,FILE=FNAM10, STATUS='OLD',FORM='UNFORMATTED')
        READ(10) HED
        READ(10) NUMNP
        READ(10) (X(I),Y(I),QL(I),BL(I),KB(I),I=1,NUMNP)
        READ(10) NUMEL
        READ(10) (KX(I,1),KX(I,2),KX(I,3),KX(I,4),
     &                              KX(I,5),ALF(I),I=1,NUMEL)
        WRITE(*,*) HED
        WRITE(*,*) '節点総数　',NUMNP,'  要素総数　',NUMEL
        WRITE(*,*)''
       WRITE(*,*)'               デ−タ−ＯＫ？　Ｙ／０　Ｎ／１'
      WRITE(*,*)''
      WRITE(*,*)'  INPUT No. ? ='
      READ(*,*)IFLAG
      IF( IFLAG .NE. 0 ) STOP 1111
C
      WRITE(*,*)'  Select data-Transfer type'
	WRITE(*,*)'        GWAP-ReStart  ------- 1'
	WRITE(*,*)'        TECPLOT  DATA ------- 2'
C	WRITE(*,*)'         ATLAS   DATA ------- 3'
	WRITE(*,*)'        FEMAP    DATA ------- 3'   !4
	WRITE(*,*)'        MicroAVS DATA ------- 4'   !5
	WRITE(*,*)'        igPOST   DATA ------- 5'   !6
	WRITE(*,*)'        INTERP   DATA ------- 6'   !7
      WRITE(*,*)'  INPUT No. ? ='
	READ(*,*) ISWT
	IF( ISWT.EQ.1) THEN
	INIT=0
        end if
       IF(ISWT .GE. 3) ISWT=ISWT+1
C
****** Added by KANO ('98.9.4)******
      IF( ISWT.EQ.3) GOTO 30
************************************
      WRITE(*,*)' '
      WRITE(*,*)'  INPUT Output File Name'
      READ (*,20) FNAME6
   20 FORMAT (A)                                                         
	OPEN (8, FILE=FNAME6, STATUS='unknown',FORM='FORMATTED')
C
      IF( ISWT.EQ.5) THEN
       OPEN(30,FILE='AVS.tmp',STATUS='UNKNOWN',FORM='FORMATTED')
      END IF       
C
	 IF( ISWT.EQ.4 ) THEN
	      CALL FEMAPO1(X,Y,KX,NUMNP,NUMEL,NUMDAT)
		    CALL FEMAPO3(KX,NUMEL,NUMDAT,TIME1,0)
	 END IF
****** Added by KANO ('98.9.4)******
   30 CONTINUE
      IF( ISWT.EQ.3) THEN
	OPEN (21, FILE='atlas.msh',STATUS='unknown',FORM='FORMATTED')
	OPEN (22, FILE='atlas.pst',STATUS='unknown',FORM='FORMATTED')
	WRITE(*,*)'  '
	WRITE(*,*)'         Mesh Data          -> atlas.msh'
	WRITE(*,*)'  Velocity and Head    Data -> atlas.pst'
	WRITE(*,*)'  '
        CALL ATLASMSH (HED,NUMNP,X,Y,NUMEL,KX,NUMDAT)
        END IF
************************************
      IF( ISWT.EQ.7 ) THEN
C
      WRITE(*,*)' '
      WRITE(*,*)'  INPUT File Name of POINT DATA'
      READ (*,20) FNAMEP
	OPEN (50, FILE=FNAMEP, STATUS='unknown',FORM='FORMATTED')
C
	READ(50,*) NPOINT
	DO I=1,NPOINT
	 READ(50,*) XPNT(I),YPNT(I)
	END DO
	CLOSE(50)
C
	DO I=1,NPOINT
	 X0=XPNT(I)
	 Y0=YPNT(I)
	 ICH=0
	 DO J=1,NUMEL
	  DO M=1,4
	   KELM(M)=KX(J,M)
        END DO
	NUS=4
	IF( KELM(3).EQ.KELM(4) .OR. KELM(4).EQ.0) NUS=3
       DO K=1,NUS-2
	  X1=X(KELM(1))
	  Y1=Y(KELM(1))
	  X2=X(KELM(K+1))
	  Y2=Y(KELM(K+1))
	  X3=X(KELM(K+2))
	  Y3=Y(KELM(K+2))
	   ICH1=0
	   ICH2=0
	   ICH3=0
	  CALL AREAE(X1,Y1,X2,Y2,X3,Y3,AA)
	  CALL AREAE(X0,Y0,X1,Y1,X2,Y2,A3)
	  CALL AREAE(X0,Y0,X2,Y2,X3,Y3,A1)
	  CALL AREAE(X0,Y0,X3,Y3,X1,Y1,A2)
        IF(A1.GE.0.0) ICH1=1
        IF(A2.GE.0.0) ICH2=1
        IF(A3.GE.0.0) ICH3=1
	ICC=ICH1*ICH2*ICH3
       IF(ICC .EQ. 0) GO TO 3000
C
	WGG(I,1)=A1/AA
	WGG(I,2)=A2/AA
	WGG(I,3)=A3/AA
	IWGT(I,1)=KELM(1)
	IWGT(I,2)=KELM(K+1)
	IWGT(I,3)=KELM(K+2)
c
c	write(*,*) i,iwgt(i,1),iwgt(i,2),iwgt(i,3)
	GO TO 3100
 3000  END DO
	 END DO
 3100	END DO
	END IF
C
	WRITE(*,*)' Input your interest	No. OF TIME STEP'
	WRITE(*,*)'   If input (<0) No., All time step avairable '
	WRITE(*,*)'   If input (0) No., Final time step avairable '
      WRITE(*,*)'  INPUT No. ? ='
	READ(*,*) LAA
C
 8000 CONTINUE
      READ(10,END=3300)LA,TIME1
      READ(10)(KD(I),H(I),QQ(I),I=1,NUMNP)
      READ(10)(VXX(I),VYY(I),DDD(I),I=1,NUMEL)
****** Added by KANO ('98.9.4)******
          ICOUNT = ICOUNT + 1
************************************
          LAAB=LAA
C
C
      IF( ISWT.EQ.7) THEN
      WRITE(8,*) "TIME,",TIME1
      WRITE(8,*) 
     1 "        X,             Y,             H"
	DO I=1,NPOINT
	 HINT=0.0
	 DO J=1,3
       HINT=HINT+WGG(I,J)*H(IWGT(I,J))
	  END DO
		 WRITE(8,1100) XPNT(I),YPNT(I),HINT
 1100    FORMAT(E14.5,',',E14.5,',',E14.5)
      END DO
      END IF
C
c      write(*,*) time1,iswt,LA,LAAB
	  IF( LAA .LT. 0 ) LAAB=LA
          IF( LAA.EQ.0. AND. EOF(10)) LAAB=LA 
c          write(*,*) iswt,LA,LAAB
	IF( ISWT.EQ.1 .AND. LA.EQ.LAAB)
     &	CALL RESTART(KD,X,Y,H,QQ,KB,QL,BL,
     &     NUMNP,LA,TIME1,INIT,NUMDAT)
C     &        CALL RESTART(KD,X,Y,H,QQ,NUMNP,LA,TIME,INIT,NUMDAT)
	IF( ISWT.EQ.2 .AND. LA.EQ.LAAB)
     &    CALL TECPLOT(HED,X,Y,H,VXX,VYY,KX,VX,VY,IC,DDD,VXD,VYD,
     &        NUMNP,NUMEL,NUMDAT,TIME1)
	IF( ISWT.EQ.3 .AND. LA.EQ.LAAB)
     &    CALL ATLASPST(ICOUNT,H,VXX,VYY,NUMNP,NUMEL,TIME1,NUMDAT)
	IF( ISWT.EQ.4 .AND. LA.EQ.LAAB) THEN
C	    CALL FEMAPO3(KX,NUMEL,NUMDAT,TIME,1)
          CALL FEMPLOT(H,VX,VY,VXX,VYY,QQ,IC,KX,DDD,VXD,VYD,
     &             NUMNP,NUMEL,NUMDAT,TIME1)
	 END IF
       IF(ISWT.EQ.5.AND. LA.EQ.LAAB) THEN
	   CALL AVS1(H,VXX,VYY,KX,VX,VY,
     &        IC,VXD,VYD,DDD,NUMNP,NUMEL,NUMDAT,TIME1)
       END IF
	IF( ISWT.EQ.6 .AND. LA.EQ.LAAB)
     &    CALL IGEOLP(HED,X,Y,H,VXX,VYY,KX,VX,VY,IC,BL,QQ,VXD,VYD,DDD,
     &                             NUMNP,NUMEL,NUMDAT,TIME1)
	GO TO 8000
****** Added by KANO ('98.9.4)******
 3300	CONTINUE
C
       IF(ISWT.EQ.5) CALL AVS2(X,Y,KX,NUMNP,NUMEL,NUMDAT)
C
      WRITE(*,1110)
      WRITE(*,*)' '
 1110 FORMAT (//,'   <<<<<  OK! OK!  NORMAL TERMINATION  >>>>> ')
************************************
      STOP
	END
	SUBROUTINE RESTART(KD,X,Y,H,QQ,KB,QL,BL,
     &                                 NUMNP,LA,TIME,INIT,NUMDAT)
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION X(NUMDAT),Y(NUMDAT),H(NUMDAT),KD(NUMDAT),QQ(NUMDAT),
     &          KB(NUMDAT),QL(NUMDAT),BL(NUMDAT)
C
      WRITE(8,*)' LA,TIME ',LA,TIME
	DO I=1,NUMNP
	YY=Y(I)*INIT
	WRITE(8,1000) I,KD(I),KB(I),X(I),Y(I),H(I),QQ(I),QL(I),BL(I)
	END DO
 1000 FORMAT(I5,2(',',I5),6(',',1PE12.4))
	RETURN
	END
      SUBROUTINE TECPLOT(HED,X,Y,H,VXX,VYY,KX,VX,VY,IC,DDD,VXD,VYD,
     &                                  NUMNP,NUMEL,NUMDAT,TIME)
	IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*4 HED(20)
	DIMENSION X(NUMDAT),Y(NUMDAT),H(NUMDAT),KX(NUMDAT,5),IC(NUMDAT),
     &       VXX(NUMDAT),VYY(NUMDAT),VX(NUMDAT),VY(NUMDAT),DDD(NUMDAT),
     &       VXD(NUMDAT),VYD(NUMDAT) 
C
      NUM=0
	DO I=1,NUMEL
  	 IF( KX(I,5).GT.0 ) NUM=NUM+1
	END DO

      WRITE(8,10) TIME,HED,NUMNP,NUM
   10 FORMAT('TITTLE = "TIME=',1PE12.4,' ',20A4,'"'/
     &       'VARIABLES = "X","Y","H","VX","VY","V",'
     "       '"VX/D","VY/D","V/D"'/
     &       'ZONE N= ',I5,', E= ',I5,', F=FEPOINT')
	DO I=1,NUMNP
	  IC(I)=0
	  VX(I)=0.0
	  VY(I)=0.0
	  VXD(I)=0.0
	  VYD(I)=0.0
	END DO
	DO J=1,NUMEL
	 KXJ5=KX(J,5)
	 IF( KXJ5 .GT.0 ) THEN
	 NUM=4
	 IF( KX(J,3).EQ.KX(J,4) ) NUM=3
	  DO M=1,NUM
	  KXJM=KX(J,M)
	  IC(KXJM)=IC(KXJM)+1
C	  IC(KXJM)=IC(KXJM)
	  VX(KXJM)=VX(KXJM)+VXX(J)
	  VY(KXJM)=VY(KXJM)+VYY(J)
	  VXD(KXJM)=VXD(KXJM)+VXX(J)/DDD(J)
	  VYD(KXJM)=VYD(KXJM)+VYY(J)/DDD(J)
	  END DO
	 END IF
	END DO
	DO I=1,NUMNP
	 IF( IC(I).GT.0 ) THEN
	   VX(I)=VX(I)/IC(I)
	   VY(I)=VY(I)/IC(I)
	   VXD(I)=VXD(I)/IC(I)
	   VYD(I)=VYD(I)/IC(I)
	 END IF
      END DO
	DO I=1,NUMNP
         VXY = DSQRT( VX(I)*VX(I)+VY(I)*VY(I) )
         VXYD = DSQRT( VXD(I)*VXD(I)+VYD(I)*VYD(I) )
	 WRITE(8,1100) X(I),Y(I),H(I),VX(I),VY(I),VXY,VXD(I),VYD(I),VXYD
	END DO
	DO I=1,NUMEL
	 IF( KX(I,5).GT.0 ) THEN 
	  WRITE(8,1000)(KX(I,J),J=1,4)
	END IF
	END DO
 1000 FORMAT(I5,3(',',I5))
 1100 FORMAT(1PE12.4,9(',',1PE12.4))
	RETURN
	END
*
*
*
*********************************************************************
* TRANSFORM gwap -RESULT TO ATLAS-POSTDATA by H.KANO '98.9 modefied *
*********************************************************************
*
      SUBROUTINE ATLASMSH (HED,NUMNP,X,Y,NUMEL,KX,NUMDAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*4 HED
      DIMENSION X(NUMDAT),Y(NUMDAT),HED(20),KX(NUMDAT,5)
*
      WRITE(21,' (A5)' )  'TITLE'
      WRITE(21,'(A12)' )  'FEM ANALISYS'
      WRITE(21,' (A6)' )  'STITLE'
      WRITE(21,'(20A4)')   HED
******************************
*     --- GRID CARD ---      *
******************************
      WRITE(21,'(A4,4X,A11)') 'GRID','(I6,3E12.0)'
      DO I=1,NUMNP
         WRITE(21,'(I6,3E12.6)') I,X(I),Y(I),0.
      END DO
      WRITE(21,'(A2)') '-1'
*******************************
*     --- CONC CARD ---       *
*******************************
      I203 = 203
      I204 = 204
      IKX5 = 100
      WRITE(21,'(A4,4X,A11)') 'CONC','(8I6)'
      DO 10 I=1,NUMEL
         I3 = KX(I,3)
         I4 = KX(I,4)
         I5 = KX(I,5)
       IF(I5.GE.0 ) THEN
         IF((I3.NE.I4).OR.(I4.NE.0)) THEN 
           WRITE(21,'(8I6)') I,I204,KX(I,5),KX(I,1),KX(I,2),KX(I,3),
     &                       KX(I,4),KX(I,5)
         ELSE
           WRITE(21,'(7I6)') I,I203,KX(I,5),KX(I,1),KX(I,2),KX(I,3),
     &                       KX(I,5)
          ENDIF
       ELSE
         IF((I3.NE.I4).OR.(I4.NE.0)) THEN 
           WRITE(21,'(8I6)') I,I204,IKX5,KX(I,1),KX(I,2),KX(I,3),
     &                       KX(I,4),IKX5
         ELSE
           WRITE(21,'(7I6)') I,I203,IKX5,KX(I,1),KX(I,2),KX(I,3),
     &                       IKX5
         ENDIF
       ENDIF
   10 CONTINUE
      WRITE(21,'(A2)') '-1'
      RETURN
      END
*--------------------------------------------------------------------
      SUBROUTINE ATLASPST(ICOUNT,H,VXX,VYY,NUMNP,NUMEL,TIME1,NUMDAT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION H(NUMDAT),VXX(NUMDAT),VYY(NUMDAT)
*
      WRITE(22,'(A7)') 'VERSION'
      WRITE(22,'(I6)')  950400
      WRITE(22,'(A4)') 'PSIZ'
      WRITE(22,'(8(I8,2X))') 0,1,3,0,0,0,0,0
      WRITE(22,'(8(I8,2X))') 0,0,0,0,0,0,0,0
      WRITE(22,'(A4)') 'STEP'
      WRITE(22,'(2(I8,2X),F12.0)') ICOUNT,1,TIME1
******************************
*     --- GVAL CARD ---      *
******************************
      WRITE(22,'(A4)') 'GVAL'
      DO I=1,NUMNP
         WRITE(22,'(I8)') I
         WRITE(22,'(6(E12.6,1X))') H(I)
      END DO
      WRITE(22,'(I2)') -1
******************************
*     --- EVAL CARD ---      *
******************************
      WRITE(22,'(A4)') 'EVAL'
      DO I=1,NUMEL
         WRITE(22,'(I8)') I
         WRITE(22,'(2(E12.6,1X))') VXX(I),VYY(I)
      END DO
      WRITE(22,'(I2)') -1
*
      RETURN
      END
      SUBROUTINE FEMPLOT(P,VXP,VZP,VXE,VZE,Q,IC,NODE,DDD,VXD,VZD,
     &             NUMNP,NUMEL,NUMDAT,TIME1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      CHARACTER*80 TITL
      DIMENSION P(NUMDAT),NODE(NUMDAT,5),
     &          VXE(NUMDAT),VZE(NUMDAT),Q(NUMDAT),
     &          VXP(NUMDAT),VZP(NUMDAT),IC(NUMDAT),
     &          VXD(NUMDAT),VZD(NUMDAT),DDD(NUMDAT)
C     ,ND(8)
      DO I=1,NUMNP
      IC(I)=0
      VXP(I)=0.0
      VZP(I)=0.0
      VXD(I)=0.0
      VZD(I)=0.0
      END DO
      DO I=1,NUMEL
      VX=VXE(I)
      VZ=VZE(I)
      DO J=1,4
       K=NODE(I,J)
	IF(K .NE. 0) THEN
       VXP(K)=VXP(K)+VX
       VZP(K)=VZP(K)+VZ
       VXD(K)=VXD(K)+VX/DDD(I)
       VZD(K)=VZD(K)+VZ/DDD(I)
       IC(K)=IC(K)+1
	END IF
      END DO
      END DO
      DO K=1,NUMNP
         ICK=IC(K)
      IF( ICK .GT.0) THEN 
       VXP(K)=VXP(K)/ICK
       VZP(K)=VZP(K)/ICK
       VXD(K)=VXD(K)/ICK
       VZD(K)=VZD(K)/ICK
      END IF
      END DO
C
       CALL  FEMAPO2
     1 (NUMEL,NUMNP,VXE,VZE,P,VXP,VZP,Q,DDD,VXD,VZD,TIME1,NUMDAT)
C
      RETURN
      END
      SUBROUTINE FEMAPO1(X,Z,NOP,MP,ME,NUMDAT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL      XX,YY,ZZ
      COMMON  / FEMAP / ISETID,IU80 ,IPSTG
C
      DIMENSION X(NUMDAT),  NOP(NUMDAT,5), Z(NUMDAT), JJ(8)
C
      NUMMAT =1
      DO I=1,NUMDAT
      IF(NUMMAT .LT.   NOP(I,5) ) NUMMAT =NOP(I,5)
      END DO
C
      IDATA = -1
      IDUM = 100
      write(IU80,100) IDATA
      write(IU80,101) IDUM
      write(IU80,*) " "
      write(IU80,*) "4.5,"
      write(IU80,100) IDATA
C
C ---------------------------------------------
      write(IU80,100) IDATA
      IDUM = 401
      write(IU80,101) IDUM
      write(IU80,*) "1,55,0,1,0,"
      write(IU80,*) " "
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,"
      write(IU80,*) "0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0,0,"
      write(IU80,*) "0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,0,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,0.,0.,"
      write(IU80,*) "0.,0.,"
      write(IU80,*) "0,0,0,0,"
      write(IU80,*) "0,0,0,0,"
      write(IU80,100) IDATA
C ---------------------------------------------
      write(IU80,100) IDATA
      IDUM = 402
      write(IU80,101) IDUM
      DO 8 J=1,NUMMAT
      write(IU80,104) J,",110,1,19,1,0,"
 104  FORMAT(I5,A14)
      write(IU80,*) " "
      write(IU80,*) "0,0,0,0,"
      write(IU80,*) "45,"
      write(IU80,*) "0,0,0,0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,0,0,0,"
      write(IU80,*) "0,0,0,0,0,"
      write(IU80,*) "93,"
      DO 9 I=1,18
      write(IU80,*) "0.,0.,0.,0.,0.,"
9     CONTINUE
      write(IU80,*) "0.,0.,0.,"
  8   CONTINUE
      write(IU80,100) IDATA
C ---------------------------------------------
C
C
      write(IU80,100) IDATA
 100  FORMAT(I5)
      IDUM = 403
      write(IU80,101) IDUM
 101  FORMAT(I6)
      DO 10 I=1,MP
         XX = X(I)
         YY = Z(I)
         ZZ = 0.0
         write(iu80,103) I,",0,0,1,46,0,0,0,0,0,0,"
     1                 ,XX,",",YY,",",ZZ
 103  format(i7,A22,1PE13.6,A1,1PE13.6,A1,1PE13.6)
10    continue
      write(IU80,100) IDATA
C
      write(IU80,100) IDATA
      IDUM = 404
c	IMAT=1
      write(IU80,101) IDUM
      do 20 i=1,ME
      JJ(1) = NOP(I,1)
      JJ(2) = NOP(I,2)
      JJ(3) = NOP(I,3)
      JJ(4) = NOP(I,4)
      ITYP = 4	
      IMAT= NOP(I,5)
	IF(JJ(3) .EQ. JJ(4) .OR. JJ(4).EQ.0) ITYP =2
C   
      write(IU80,*) i,",124,",IMAT,",19,",ITYP,",1,0,0,"
      IF(ITYP .EQ. 2) THEN
      write(IU80,102)  JJ(1),JJ(2),JJ(3),0,0,0,0,0,0,0
      ELSEIF(ITYP .EQ. 4) THEN
      write(IU80,102) JJ(1),JJ(2),JJ(3),JJ(4),0,0,0,0,0,0
	END IF
  102 format(10(i7,','))
C      
      write(IU80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0.,0.,0.,"
      write(IU80,*) "0,0,0,0,0,0,0,0,0,0,0,0,"
20    continue
      write(IU80,100) IDATA
C
      RETURN
      END
C
C ***************************************************************
      SUBROUTINE FEMAPO2
     1 (NUMEL,NUMNP,VXE,VZE,P,VXP,VZP,Q,DDD,VXD,VZD,TIME1,NUMDAT)
C ****************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL      PPP,TTT,DVX,dum1,dum2,dum3,DVY,DVZ
      COMMON  / FEMAP / ISETID,IU80 ,IPSTG
C
      DIMENSION P(NUMDAT),VXE(NUMDAT),VZE(NUMDAT)
      DIMENSION VXP(NUMDAT),VZP(NUMDAT),Q(NUMDAT)
      DIMENSION VXD(NUMDAT),VZD(NUMDAT),DDD(NUMDAT)
C
C     
      ISETID = ISETID + 1
C
      IDUM = 450
      IDATA = -1
      write(IU80,100) IDATA
100   FORMAT(I5)
      write(IU80,101) IDUM
101   FORMAT(I6)
      write(iu80,*) ISETID 
C
      TTT = TIME1
      write(iu80,102) TTT
102   format(" T=",1PE12.4)
C
      write(iu80,*)  "0,0,"
      write(iu80,*)  TTT
      write(iu80,*) "1,"
      write(iu80,*) "  "
      write(IU80,100) IDATA
C
      pmax = P(1)
      pmin = pmax
      ipmax = 1
      ipmin = 1
C
      do 10 i=1,NUMNP
       pp =   P(i)
C
       if( pmax .lt. pp) then
             pmax = pp
             ipmax=i
       end if
       if( pmin .gt. pp) then
             pmin = pp
             ipmin=i
       end if
C
10    continue
C ------------------------------------------
      IDUM = 451
      write(IU80,100) IDATA
      write(IU80,101) IDUM
C ------------------------------------------
      write(IU80,*) ISETID, ",10,1,"
      write(IU80,*) "Total Head"
      dum1 = pmin
      dum2 = pmax
      dum3 = ABS(pmax)
      if(dum3 .lt. abs(pmin)) dum3=abs(pmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ipmin,",",ipmax,",6,7,"
      write(IU80,*) "0,1,1,"
      do 55 i=1,NUMNP
       ppp =  P(i)
       write(IU80,*) i,",",ppp
 55    continue
       i=-1
       ppp=0.0
       write(IU80,*) i,",",ppp
C ------------------------------------------
      pmax = Q(1)
      pmin = pmax
      ipmax = 1
      ipmin = 1
C
      do 60 i=1,NUMNP
       pp =   Q(i)
C
       if( pmax .lt. pp) then
             pmax = pp
             ipmax=i
       end if
       if( pmin .gt. pp) then
             pmin = pp
             ipmin=i
       end if
C
60    continue
C ------------------------------------------
      write(IU80,*) ISETID, ",12,1,"
      write(IU80,*) "Q"
      dum1 = pmin
      dum2 = pmax
      dum3 = ABS(pmax)
      if(dum3 .lt. abs(pmin)) dum3=abs(pmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ipmin,",",ipmax,",6,7,"
      write(IU80,*) "0,1,1,"   
      do 62 i=1,NUMNP 
       ppp =  Q(i)
       write(IU80,*)i,",",ppp
 62    continue
       i=-1
       ppp=0.0
       write(IU80,*) i,",",ppp
C ------------------------------------------
       xmin = vxp(1)
       xmax = xmin
       ixmin = 1
       ixmax = 1
       zmin = vzp(1)
       zmax = zmin
       izmin = 1
       izmax = 1
       vmin = sqrt(xmin*xmin+zmin*zmin)
       vmax = vmin
       ivmin = 1
       ivmax = 1
       do 73 i=1,NUMNP
         xv = vxp(i)
         zv = vzp(i)
         vv = sqrt(xv*xv+zv*zv)
         if( xmax .lt. xv) then
               xmax = xv
               ixmax=i
         end if
         if( xmin .gt. xv) then
               xmin = xv
               ixmin=i
         end if
         if( zmax .lt. zv) then
               zmax = zv
               izmax=i
         end if
         if( zmin .gt. zv) then
               zmin = zv
               izmin=i
         end if
         if( vmax .lt. vv) then
               vmax = vv
               ivmax=i
         end if
         if( vmin .gt. vv) then
               vmin = vv
               ivmin=i
         end if
 73   continue
C ------------------------------------------
      write(IU80,*) ISETID, ",13,1,"
      write(IU80,*) "D*Velocity(Node)"
      dum1 = vmin
      dum2 = vmax
      dum3 = vmax
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "14,15,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ivmin,",",ivmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 80 i=1,NUMNP
       xv = vxp(I)
       zv = vzp(I)
       vv = sqrt(xv*xv+zv*zv)
       xv = vv
       write(IU80,*) i,",",xv
 80    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C ------------------------------------------
      write(IU80,*) ISETID, ",14,1,"
      write(IU80,*) "D*X Velocity(Node)"
      dum1 = xmin
      dum2 = xmax
      dum3 = xmax
      dum3 = ABS(xmax)
      if(dum3 .lt. abs(xmin)) dum3=abs(xmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "14,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ixmin,",",ixmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 90 i=1,NUMNP
       xv = vxp(I)
       write(IU80,*) i,",",xv
 90    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C ------------------------------------------
      write(IU80,*) ISETID, ",15,1,"
      write(IU80,*) "D*Y Velocity(Node)"
      dum1 = zmin
      dum2 = zmax
      dum3 = zmax
      dum3 = ABS(zmax)
      if(dum3 .lt. abs(zmin)) dum3=abs(zmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0,16,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) izmin,",",izmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 120 i=1,NUMNP
       xv = vzp(I)
       write(IU80,*) i,",",xv
120    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C
C ------------------------------------------
       xmin = vxD(1)
       xmax = xmin
       ixmin = 1
       ixmax = 1
       zmin = vzD(1)
       zmax = zmin
       izmin = 1
       izmax = 1
       vmin = sqrt(xmin*xmin+zmin*zmin)
       vmax = vmin
       ivmin = 1
       ivmax = 1
       do 173 i=1,NUMNP
         xv = vxD(i)
         zv = vzD(i)
         vv = sqrt(xv*xv+zv*zv)
         if( xmax .lt. xv) then
               xmax = xv
               ixmax=i
         end if
         if( xmin .gt. xv) then
               xmin = xv
               ixmin=i
         end if
         if( zmax .lt. zv) then
               zmax = zv
               izmax=i
         end if
         if( zmin .gt. zv) then
               zmin = zv
               izmin=i
         end if
         if( vmax .lt. vv) then
               vmax = vv
               ivmax=i
         end if
         if( vmin .gt. vv) then
               vmin = vv
               ivmin=i
         end if
 173   continue
C ------------------------------------------
      write(IU80,*) ISETID, ",16,1,"
      write(IU80,*) "Velocity(Node)"
      dum1 = vmin
      dum2 = vmax
      dum3 = vmax
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "17,18,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ivmin,",",ivmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 180 i=1,NUMNP
       xv = vxD(I)
       zv = vzD(I)
       vv = sqrt(xv*xv+zv*zv)
       xv = vv
       write(IU80,*) i,",",xv
 180    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C ------------------------------------------
      write(IU80,*) ISETID, ",17,1,"
      write(IU80,*) "X Velocity(Node)"
      dum1 = xmin
      dum2 = xmax
      dum3 = xmax
      dum3 = ABS(xmax)
      if(dum3 .lt. abs(xmin)) dum3=abs(xmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "17,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ixmin,",",ixmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 190 i=1,NUMNP
       xv = vxD(I)
       write(IU80,*) i,",",xv
 190    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C ------------------------------------------
      write(IU80,*) ISETID, ",18,1,"
      write(IU80,*) "Y Velocity(Node)"
      dum1 = zmin
      dum2 = zmax
      dum3 = zmax
      dum3 = ABS(zmax)
      if(dum3 .lt. abs(zmin)) dum3=abs(zmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0,18,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) izmin,",",izmax,",1,7,"
      write(IU80,*) "0,1,1,"
      do 1120 i=1,NUMNP
       xv = vzD(I)
       write(IU80,*) i,",",xv
1120    continue
       i=-1
       xv=0.0
       write(IU80,*) i,",",xv
C
C ------------------------------------------
C
C ------------------------------------------
C
C ------------------------------------------
       xmin =  VXE(1)
       xmax = xmin
       ixmin =1
       ixmax = 1
       zmin =   VZE(1)
       zmax = zmin
       izmin = 1
       izmax =1
       tmin = SQRT(XMAX*XMAX+ZMAX*ZMAX)
       tmax = tmin
       itmin = 1
       itmax = 1
C
       do 185 i=1,NUMEL
         xv =  VXE(i)
         zv =  VZE(i)
         tv = SQRT(XV*XV+zv*zv)
C
         if( xmax .lt. xv) then
               xmax = xv
               ixmax=i
         end if
         if( xmin .gt. xv) then
               xmin = xv
               ixmin=i
         end if
C
         if( zmax .lt. zv) then
               zmax = zv
               izmax=i
         end if
         if( zmin .gt. zv) then
               zmin = zv
               izmin=i
         end if
C
         if( tmax .lt. tv) then
               tmax = tv
               itmax=i
         end if
         if( tmin .gt. tv) then
               tmin = tv
               itmin=i
         end if
C
 185  continue
C ------------------------------------------
      write(IU80,*) ISETID, ",70010,1,"
      write(IU80,*) "D*Velocity"
      dum1 = tmin
      dum2 = tmax
      dum3 = ABS(tmax)
      if(dum3 .lt. abs(tmin)) dum3=abs(tmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "70011,70012,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) itmin,",",itmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 181 i=1,NUMEL
       Dvx =   VXE(I)
       Dvz =    VZE(I)
       Dvx =SQRT(Dvx*Dvx+Dvz*Dvz)
       write(IU80,*) i,",",Dvx
 181  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
      write(IU80,*) ISETID, ",70011,1,"
      write(IU80,*) "D*X Velocity"
      dum1 = xmin
      dum2 = xmax
      dum3 = ABS(xmax)
      if(dum3 .lt. abs(xmin)) dum3=abs(xmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "70011,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ixmin,",",ixmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 182 i=1,NUMEL
       Dvx =  VXE(I)
       write(IU80,*) i,",",Dvx
 182  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
      write(IU80,*) ISETID, ",70012,1,"
      write(IU80,*) "D*Y Velocity"
      dum1 = zmin
      dum2 = zmax
      dum3 = ABS(zmax)
      if(dum3 .lt. abs(zmin)) dum3=abs(zmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0.70012,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) izmin,",",izmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 184 i=1,NUMEL
       Dvx = VZE(I)
       write(IU80,*) i,",",Dvx
 184  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
C ------------------------------------------
       xmin =  VXE(1)/DDD(1)
       xmax = xmin
       ixmin =1
       ixmax = 1
       zmin =   VZE(1)/DDD(1)
       zmax = zmin
       izmin = 1
       izmax =1
       tmin = SQRT(XMAX*XMAX+ZMAX*ZMAX)
       tmax = tmin
       itmin = 1
       itmax = 1
C
       do 285 i=1,NUMEL
         xv =  VXE(i)/DDD(i)
         zv =  VZE(i)/DDD(i)
         tv = SQRT(XV*XV+zv*zv)
C
         if( xmax .lt. xv) then
               xmax = xv
               ixmax=i
         end if
         if( xmin .gt. xv) then
               xmin = xv
               ixmin=i
         end if
C
         if( zmax .lt. zv) then
               zmax = zv
               izmax=i
         end if
         if( zmin .gt. zv) then
               zmin = zv
               izmin=i
         end if
C
         if( tmax .lt. tv) then
               tmax = tv
               itmax=i
         end if
         if( tmin .gt. tv) then
               tmin = tv
               itmin=i
         end if
C
 285  continue
C ------------------------------------------
      write(IU80,*) ISETID, ",70013,1,"
      write(IU80,*) "Velocity"
      dum1 = tmin
      dum2 = tmax
      dum3 = ABS(tmax)
      if(dum3 .lt. abs(tmin)) dum3=abs(tmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "70014,70015,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) itmin,",",itmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 281 i=1,NUMEL
       Dvx =   VXE(I)/DDD(i)
       Dvz =    VZE(I)/DDD(i)
       Dvx =SQRT(Dvx*Dvx+Dvz*Dvz)
       write(IU80,*) i,",",Dvx
 281  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
      write(IU80,*) ISETID, ",70014,1,"
      write(IU80,*) "X Velocity"
      dum1 = xmin
      dum2 = xmax
      dum3 = ABS(xmax)
      if(dum3 .lt. abs(xmin)) dum3=abs(xmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "70014,0,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) ixmin,",",ixmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 282 i=1,NUMEL
       Dvx =  VXE(I)/DDD(i)
       write(IU80,*) i,",",Dvx
 282  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
      write(IU80,*) ISETID, ",70015,1,"
      write(IU80,*) "Y Velocity"
      dum1 = zmin
      dum2 = zmax
      dum3 = ABS(zmax)
      if(dum3 .lt. abs(zmin)) dum3=abs(zmin)
      write(IU80,*) dum1,dum2,dum3
      write(Iu80,*) "0,70015,0,0,0,0,0,0,0,0,"
      write(Iu80,*) "0,0,0,0,0,0,0,0,0,0,"
      write(IU80,*) izmin,",",izmax,",1,8,"
      write(IU80,*) "0,1,1,"
      do 284 i=1,NUMEL
       Dvx = VZE(I)/DDD(i)
       write(IU80,*) i,",",Dvx
 284  continue
       i=-1
       Dvx=0.0
       write(IU80,*) i,",",Dvx
C ------------------------------------------
C
      write(IU80,100) IDATA
C
      RETURN
      END
C
      SUBROUTINE FEMAPO3(NOP,NUMEL,NUMDAT,TIME1,ISW)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C  
      COMMON  / FEMAP / ISETID,IU80 ,IPSTG
	DIMENSION   NOP(NUMDAT,5) 
C
      ipstg=ipstg+1
      IDUM = 408
      IDATA = -1
      write(IU80,100) IDATA
100   FORMAT(I5)
      write(IU80,101) IDUM
101   FORMAT(I6)
      write(iu80,*) ipstg,",1"
C
      IF(ISW .EQ. 0) THEN
      write(iu80,103) 
103   format("ELEM-GROUP 0 ")
	ELSE
      TTT = TIME1
      write(iu80,102) TTT
102   format("EGROUP ",1PE12.4)
      END IF
C
      write(iu80,*) "0,0,0"
      write(Iu80,*) "0,0,0,0,0.,0."
      write(iu80,*) "0,0"
      DO 10 I=1,6
      write(iu80,*) "0,0"
      write(Iu80,*) "0.,0.,0."
      write(iu80,*) "0.,0.,0."
   10 CONTINUE
      write(iu80,*) "38"
      write(iu80,*) "21"
        DO 20 I=1,NUMEL
         IF(NOP(I,5) .GT. 0) THEN
            write(iu80,*) I,",","0,1,1,"
         END IF
   20 CONTINUE
      write(Iu80,*) "-1,-1,-1,-1,"
        write(iu80,*) "-1"
        write(iu80,*) "15"
        write(iu80,*) "-1"
C ------------------------------------------
      write(IU80,100) IDATA
C
      RETURN
      END
      SUBROUTINE AVS1(HHH,VXE,VYE,KX,VXP,VYP,IC,VXD,VYD,DDD,
     &             NUMNP,NUMEL,NUMDAT,TIME1)
      IMPLICIT REAL*8 ( A-H,O-Z )
      REAL CCX(10)
	COMMON /AVSCC / ISTPNUM 
      DIMENSION HHH(NUMDAT),KX(NUMDAT,5),
     &          VXE(NUMDAT),VYE(NUMDAT),
     &          VXP(NUMDAT),VYP(NUMDAT),IC(NUMDAT),
     &          VXD(NUMDAT),VYD(NUMDAT),DDD(NUMDAT)
C
      IU6 =8
	IU8 = 30
C
	DO I=1,NUMNP
	  IC(I)=0
	  VXP(I)=0.0
	  VYP(I)=0.0
	  VXD(I)=0.0
	  VYD(I)=0.0
	END DO
	DO J=1,NUMEL
	 KXJ5=KX(J,5)
	 IF( KXJ5 .GT.0 ) THEN
	 NUM=4
	 IF( KX(J,3).EQ.KX(J,4) ) NUM=3
	  DO M=1,NUM
	  KXJM=KX(J,M)
	  IC(KXJM)=IC(KXJM)+1
	  VXP(KXJM)=VXP(KXJM)+VXE(J)
	  VYP(KXJM)=VYP(KXJM)+VYE(J)
	  VXD(KXJM)=VXD(KXJM)+VXE(J)/DDD(J)
	  VYD(KXJM)=VYD(KXJM)+VYE(J)/DDD(J)
	  END DO
	 END IF
	END DO
	DO I=1,NUMNP
	 IF( IC(I).GT.0 ) THEN
	   VXP(I)=VXP(I)/IC(I)
	   VYP(I)=VYP(I)/IC(I)
	   VXD(I)=VXD(I)/IC(I)
	   VYD(I)=VYD(I)/IC(I)
	 END IF
      END DO
C
      ISTPNUM = ISTPNUM+1
C
	IF(ISTPNUM.GE. 1000) THEN
	  WRITE(IU8,6011) ISTPNUM,TIME1
 6011   FORMAT("step",I4,E15.7)
	ELSEIF(ISTPNUM.GE. 100) THEN
	  WRITE(IU8,6012) ISTPNUM,TIME1
 6012   FORMAT("step",I3,E15.7)
	ELSEIF(ISTPNUM.GE. 10) THEN	 
	  WRITE(IU8,6013) ISTPNUM,TIME1
 6013   FORMAT("step",I2,E15.7)
      ELSE
	  WRITE(IU8,6000) ISTPNUM,TIME1
 6000   FORMAT("step",I1,E15.7)
	END IF
C
	  WRITE(IU8,6002) 
 6002 FORMAT("7   6"/"7   1   1   1   1   1   1   1")
      WRITE(IU8,6001)
 6001 FORMAT('Total Head,'/'D*Velocity X,'/
     1       'D*Velocity Y,'/'D*Velocity,'/
     1       'Velocity X,'/'Velocity Y,'/'Velocity,')
      DO I=1,NUMNP
	CCX(1) = HHH(I)
	CCX(2) = VXP(I)
	CCX(3) = VYP(I)
	CCX(4) = DSQRT( VXP(I)*VXP(I)+VYP(I)*VYP(I) )
	CCX(5) = VXD(I)
	CCX(6) = VYD(I)
	CCX(7) = DSQRT( VXD(I)*VXD(I)+VYD(I)*VYD(I) )
C
      WRITE(IU8,6003)I,(CCX(II),II=1,7)
 6003 FORMAT(I8,1X,7(1PE12.5," "))
      END DO
	  WRITE(IU8,6005) 
 6005 FORMAT("6   1   1   1   1   1   1")
      WRITE(IU8,6006)
 6006 FORMAT('D*Velocity X,'/'D*Velocity Y,'/'D*Velocity,'/
     1        'Velocity X,'/'Velocity Y,'/'Velocity,')
	  DO I=1,NUMEL
	CCX(1) = VXE(I)
	CCX(2) = VYE(I)
	CCX(3) =  DSQRT( VXE(I)*VXE(I)+VYE(I)*VYE(I) )
	CCX(4) = VXE(I)/DDD(I)
	CCX(5) = VYE(I)/DDD(I)
	CCX(6) =  DSQRT( VXE(I)*VXE(I)+VYE(I)*VYE(I) )/DDD(I)
	  WRITE(IU8,6004) I,(CCX(II),II=1,6)
 6004 FORMAT(I8,1X,6(1PE12.5," "))
	  END DO
      RETURN
      END
      SUBROUTINE AVS2(X,Y,NODE,NUMNP,NUMEL,NUMDAT)
      IMPLICIT REAL*8 ( A-H,O-Z )
      CHARACTER*1 CCC(200)
	COMMON /AVSCC / ISTPNUM 
      DIMENSION X(NUMDAT),Y(NUMDAT),Z(NUMNP),NODE(NUMDAT,5)
C
      IU6 = 8
	IU8 = 30
C
      WRITE(IU6,6000) 
 6000 FORMAT("# AC-GWAP ")
C
      WRITE(IU6,6001)  ISTPNUM
 6001 FORMAT(I5)    
      WRITE(IU6,6002)  
 6002 FORMAT("data") 
C
      REWIND IU8
      READ(IU8,5000) CCC
 5000 FORMAT(200A1)
      WRITE(IU6,6003) (CCC(I),I=1,26)
 6003 FORMAT(200A1)
      WRITE(IU6,6004) NUMNP,NUMEL
 6004 FORMAT(I10,1X,I10)
C
      ZZ =0.0
      DO I=1,NUMNP
      WRITE(IU6,6005) I,X(I),Y(I) ,ZZ
 6005 FORMAT(I8,3(' ',1PE12.5))
      END DO
C
      DO I=1,NUMEL
	IF(NODE(I,3) .EQ.  NODE(I,4)) THEN
	  WRITE(IU6,6006) I,NODE(I,5),(NODE(I,J),J=1,3)
 6006 FORMAT(I7," ",I3," tri ",3(" ",I7))
	ELSE
	  WRITE(IU6,6007) I,NODE(I,5),(NODE(I,J),J=1,4)
 6007 FORMAT(I7," ",I3," quad ",4(" ",I7))
	END IF
      END DO
C
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,5)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,29)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,11)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
C
      DO II=1,NUMNP
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,100)
	END DO
C
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,25)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
C
      DO II=1,NUMEL
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,87)
	END DO
C
      IF(ISTPNUM .LE. 1) RETURN
      DO II=1,ISTPNUM-1
C
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,26)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,5)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,29)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,11)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
C
      DO I=1,NUMNP
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(IJ),IJ=1,100)
	END DO
C
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,25)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(I),I=1,13)
C
      DO I=1,NUMEL
      READ(IU8,5000) CCC
      WRITE(IU6,6003) (CCC(IJ),IJ=1,87)
	END DO
	END DO
C
      RETURN
      END
C
      SUBROUTINE IGEOLP(HED,X,Y,H,VXX,VYY,KX,VX,VY,IC,BL,QQ,VXD,VYD,DDD,
     &                                  NUMNP,NUMEL,NUMDAT,TIME1)
	IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*4 HED(20)
	DIMENSION X(NUMDAT),Y(NUMDAT),H(NUMDAT),KX(NUMDAT,5),IC(NUMDAT),
     &          VXX(NUMDAT),VYY(NUMDAT),VX(NUMDAT),VY(NUMDAT),
     &       BL(NUMDAT),QQ(NUMDAT),VXD(NUMDAT),VYD(NUMDAT),DDD(NUMDAT)
C
      NUM=0
	DO I=1,NUMEL
  	 IF( KX(I,5).GT.0 ) NUM=NUM+1
	END DO
      WRITE(8,11) TIME1 
   11 FORMAT("TIME=",1PE12.4)
      WRITE(8,10) HED,NUMNP,NUM
   10 FORMAT('TITTLE= "',20A4,'"'/
     &       'ZONE N= ',I5,', E= ',I5,',M= 1 ,NV= 9 ,EV= 5')
      WRITE(8,12) 
      ISW = 0
   12 FORMAT('NVARIABLES = "SW","X","Y","H","VX","VY","V"'
     1             ',"VX*D","VY*D","V*D","BL","Q"')
	DO I=1,NUMNP
	  IC(I)=0
	  VX(I)=0.0
	  VY(I)=0.0
	  VXD(I)=0.0
	  VYD(I)=0.0
	END DO
	DO J=1,NUMEL
	 KXJ5=KX(J,5)
	 IF( KXJ5 .GT.0 ) THEN
	 NUM=4
	 IF( KX(J,3).EQ.KX(J,4) ) NUM=3
	  DO M=1,NUM
	  KXJM=KX(J,M)
	  IC(KXJM)=IC(KXJM)+1
C	  IC(KXJM)=IC(KXJM)
	  VX(KXJM)=VX(KXJM)+VXX(J)
	  VY(KXJM)=VY(KXJM)+VYY(J)
	  VXD(KXJM)=VXD(KXJM)+VXX(J)/DDD(J)
	  VYD(KXJM)=VYD(KXJM)+VYY(J)/DDD(J)
	  END DO
	 END IF
	END DO
	DO I=1,NUMNP
	 IF( IC(I).GT.0 ) THEN
	   VX(I)=VX(I)/IC(I)
	   VY(I)=VY(I)/IC(I)
	   VXD(I)=VXD(I)/IC(I)
	   VYD(I)=VYD(I)/IC(I)
	 END IF
      END DO
	DO I=1,NUMNP
         VXY = DSQRT( VX(I)*VX(I)+VY(I)*VY(I) )
         VXYD = DSQRT( VXD(I)*VXD(I)+VYD(I)*VYD(I) )
	 WRITE(8,1100) ISW,X(I),Y(I),
     1      H(I),VXD(I),VYD(I),VXYD,VX(I),VY(I),VXY,BL(I),QQ(I)
	END DO
      WRITE(8,13) 
   13 FORMAT('EVARIABLES = ',
     1  '"SW","I","J","K","L","P","X","Y","VX","VY","V"')
      MAT = 1
      ISW = 0
	DO I=1,NUMEL
	 IF( KX(I,5).GT.0 ) THEN 
c         VXY = DSQRT( VXX(I)*VXX(I)+VYY(I)*VYY(I) )
         VXYD = DSQRT( VXX(I)*VXX(I)+VYY(I)*VYY(I) )/DDD(I)
C
       XG= 0.0
	 YG= 0.0
	 NUM=4
	 IF( KX(I,3).EQ.KX(I,4) ) NUM=3
	  DO M=1,NUM
	  KXJM=KX(I,M)
	  XG=XG+X(KXJM)
	  YG=YG+Y(KXJM)
	  END DO
      XG =XG/FLOAT(NUM)
	YG =YG/FLOAT(NUM)
C
	  WRITE(8,1000) ISW,(KX(I,J),J=1,4),MAT,XG,YG,
     1            VXX(I)/DDD(I),VYY(I)/DDD(I),VXYD
	END IF
	END DO
 1000 FORMAT(I2,5(',',I5),5(',',1PE12.4))
 1100 FORMAT(I2,',',1PE12.4,10(',',1PE12.4))
	RETURN
	END
      SUBROUTINE AREAE(X1,Y1,X2,Y2,X3,Y3,AA)
	IMPLICIT REAL*8 (A-H,O-Z)
C      CI=X3-X2
      CJ=X1-X3
      CK=X2-X1
C      BI=Y2-Y3
      BJ=Y3-Y1
      BL=Y1-Y2
      AA=(CK*BJ-CJ*BL)/2.0
	RETURN
	END
