!TWO ORDER WENO SCHEME ON TWO DIMENSIONAL UNSTRUCTURED MESHES
PROGRAM MAIN

!CALCULATE DOMAIN:[XL,XR]*[YB,YT]
!THE NUMBER OF GRID:MGRID
PARAMETER(XL=-2.0,XR=2.0,YB=-2.0,YT=2.0)
PARAMETER(MG=2500,MN=2000)
PARAMETER(PI=3.1415926)
PARAMETER(ALPHA=1.0)
PARAMETER(TIME=10)

!COORDERNATE OF VERTEX:X,Y
!THE NUMBER OF TRIANGLES INTERSECT AT VERTEX AND THEIR No. AND 
!   INTERIOR ANGLE:TRIANGLE_MATRIX(6)INTERIOR_ANGLE(6),NUM
!FUNCTION VALUE OF VERTEX:VERTEX_VALUE 
!CALCULATING CONSTANT:BETA(6),NL(2,6)
TYPE VERTEX_TYPE 
DOUBLE PRECISION X
DOUBLE PRECISION Y
INTEGER NUM
INTEGER TRIMATRIX(6)
DOUBLE PRECISION INTERANGLE(6)
DOUBLE PRECISION BETA(6)
DOUBLE PRECISION NL(2,6)
DOUBLE PRECISION VALUE
END TYPE

!THE No. FO INTERIOR ANGLE OF TRIANGLE:NODE(3)
!THE NUMBER OF INTERPOLATE POLYNOMIAL CONSTANT:MODULE_NUM,
!MODULE_CONSTANT(4,6,6),MODULE_NODE(4,6)
!THE POLYNOMIAL:POLYNOMIAL_CONST(6) (X^2,XY,Y^2,X,Y,1)

TYPE TRIANGLE_TYPE
INTEGER NODE(3)
DOUBLE PRECISION TAREA
INTEGER MONUM
DOUBLE PRECISION MOCONST(4,6,6)
INTEGER MONODE(4,6)
!A1+A2*X+A3*Y+A4*X^2+A5*X*Y+A6*Y^2
DOUBLE PRECISION PCONST(6)
END TYPE

TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER I,KNODE,KGRID,VNODE,VGRID
DOUBLE PRECISION DETT

!grid divide
!對三角形網格及頂點進行編號
CALL GENERATE_GRID(GRID,NODE,KNODE,KGRID)
!尋找每個接點周圍的三角形,記錄其個數及編號
CALL ENCODE_TRIANGLE(GRID,NODE,KGRID)
!conduct boundary 
!利用周期性邊界條件在邊界上添加網格和接點
!虛點結構中num位置存儲的是其對應實點的編號
!虛網格結構MONUM中存儲的是其對應的實網格的編號
CALL CONDUCT_BOUNDARY(GRID,NODE,KNODE,KGRID,VNODE,VGRID)
!尋找每個接點周圍的三角形,記錄其個數及編號
CALL ENCODE_TRIANGLE(GRID,NODE,VGRID)

!prepare
!求各三角形上適合迭代的子模塊
CALL CALCULATE_MODULE(GRID,NODE,KGRID,VGRID)
!求各個網格的面積
!CALL CALCULATE_GRID_AREA(GRID,VGRID,NODE)
!求各頂點周圍的三角形及計算通量的系數
CALL CALCULATE_NODE_CONST(GRID,NODE,KNODE)

!對節點賦初值
CALL INITIAL(NODE,VNODE)
!two-step RK iteration
DETT=0.01521
DO I=1,TIME
!對每個頂點,運用L-F通量進行迭代
  CALL ITERATION(GRID,NODE,KNODE,KGRID,DETT,VNODE)
END DO

!output result
CALL DATA_OUTPUT(NODE,GRID,KGRID,VNODE) 

CONTAINS

SUBROUTINE CALCULATE_VIRTUAL_NODE(NODE,KNODE,VNODE)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
INTEGER KNODE,VNODE,I
DO I=KNODE+1,VNODE
  NODE(I).VALUE=NODE(NODE(I).NUM).VALUE
END DO
END SUBROUTINE CALCULATE_VIRTUAL_NODE

!計算各點對應的數值HAMILTON通量
SUBROUTINE CALCULATE_NUMERICAL_HAMILTON(GRID,NODE,KNODE,KGRID,H)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KNODE,KGRID,I,J,K,N
DOUBLE PRECISION S,P(6),H(KNODE),X,Y,DUDX(6),DUDY(6)
DO I=1,KNODE
  DO J=1,NODE(I).NUM
    K=NODE(I).TRIMATRIX(J)
    IF(K<=KGRID)THEN
	  P=GRID(K).PCONST
	  X=NODE(I).X
	  Y=NODE(I).Y
	ELSE
	  P=GRID(GRID(K).MONUM).PCONST
	  N=1
	  DO WHILE((N<=3).AND.(GRID(K).NODE(N)/=I))
	    N=N+1
	  END DO
	  N=GRID(GRID(K).MONUM).NODE(N)
	  X=NODE(N).X
	  Y=NODE(N).Y
	END IF
	DUDX(J)=P(2)+2.0*P(4)*X+P(5)*Y
	DUDY(J)=P(3)+P(5)*X+2.0*P(6)*Y
  END DO
  X=0
  Y=0
  DO J=1,NODE(I).NUM
    X=X+NODE(I).INTERANGLE(J)*DUDX(J)
	Y=Y+NODE(I).INTERANGLE(J)*DUDY(J)
  END DO
  X=X/(2.0*PI)
  Y=Y/(2.0*PI)
  H(I)=HAMILTON(X,Y)
  S=0
  DO J=1,NODE(I).NUM
    IF(J==NODE(I).NUM)THEN
	  N=1
	ELSE
	  N=J+1
	END IF
	X=(DUDX(J)+DUDX(N))/2.0
	Y=(DUDY(J)+DUDY(N))/2.0
    S=S+NODE(I).BETA(J)*(X*NODE(I).NL(1,J)+Y*NODE(I).NL(2,J))
  END DO
  H(I)=H(I)-S*ALPHA/PI
END DO
END SUBROUTINE CALCULATE_NUMERICAL_HAMILTON

SUBROUTINE ITERATION(GRID,NODE,KNODE,KGRID,DETT,VNODE)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KNODE,KGRID,VNODE,I
DOUBLE PRECISION DETT,H(KNODE),U0(KNODE)
!對每個三角形網格計算其上的WENO多項式
CALL CALCULATE_POLYNOMIAL(GRID,KGRID,NODE)
!計算各點對應的數值HAMILTON通量
CALL CALCULATE_NUMERICAL_HAMILTON(GRID,NODE,KNODE,KGRID,H)
DO I=1,KNODE
  U0(I)=NODE(I).VALUE
  NODE(I).VALUE=NODE(I).VALUE-DETT*H(I)
END DO
!給虛接點賦值
CALL CALCULATE_VIRTUAL_NODE(NODE,KNODE,VNODE)
!對每個三角形網格計算其上的WENO多項式
CALL CALCULATE_POLYNOMIAL(GRID,KGRID,NODE)
!計算各點對應的數值HAMILTON通量
CALL CALCULATE_NUMERICAL_HAMILTON(GRID,NODE,KNODE,KGRID,H)
DO I=1,KNODE
  NODE(I).VALUE=0.5*(U0(I)+NODE(I).VALUE-DETT*H(I))
END DO
!給虛接點賦值
CALL CALCULATE_VIRTUAL_NODE(NODE,KNODE,VNODE)
END SUBROUTINE ITERATION

SUBROUTINE CALCULATE_GRID_AREA(GRID,VGRID,NODE)
IMPLICIT NONE
EXTERNAL AREA
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER I,VGRID,A(3)
DOUBLE PRECISION AREA
DO I=1,VGRID
  A=GRID(I).NODE
  GRID(I).TAREA=AREA(NODE(A(1)).X,NODE(A(1)).Y,&
    NODE(A(2)).X,NODE(A(2)).Y,NODE(A(3)).X,NODE(A(3)).Y)
END DO
END SUBROUTINE CALCULATE_GRID_AREA

SUBROUTINE CALCULATE_POLYNOMIAL(GRID,KGRID,NODE)
IMPLICIT NONE
EXTERNAL SOLVER
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER I,J,K,KGRID
DOUBLE PRECISION S,WEIGHT(4),C(4,6),A(6,6),B(6),EPS
PARAMETER(EPS=1E-7)
DO I=1,KGRID
  S=0
  DO J=1,GRID(I).MONUM
    A=GRID(I).MOCONST(J,1:6,1:6)
	DO K=1,6
	  B(K)=NODE(GRID(I).MONODE(J,K)).VALUE
	END DO
	CALL SOLVER(A,B,6)
	C(J,1:6)=B
	WEIGHT(J)=(2.0*B(4))**2.0+(B(5))**2.0+(2.0*B(6))**2.0
	WEIGHT(J)=1.0/(WEIGHT(J)+EPS)
	S=S+WEIGHT(J)
  END DO
  DO J=1,6
    GRID(I).PCONST(J)=0
  END DO
  DO J=1,GRID(I).MONUM
    WEIGHT(J)=WEIGHT(J)/S
  END DO
  DO J=1,GRID(I).MONUM
	DO K=1,6
	  GRID(I).PCONST(K)=GRID(I).PCONST(K)+WEIGHT(J)*C(J,K)
	END DO
  END DO
END DO
END SUBROUTINE CALCULATE_POLYNOMIAL

!--------------------------------------------------------------
FUNCTION U0_VALUE(X,Y)
IMPLICIT NONE
DOUBLE PRECISION X,Y,U0_VALUE
U0_VALUE=SIN((X+Y)*PI/2.0)
!U0_VALUE=-COS(PI*(X+Y)/2.0)
END FUNCTION U0_VALUE

FUNCTION HAMILTON(DUDX,DUDY)
IMPLICIT NONE
DOUBLE PRECISION DUDX,DUDY,HAMILTON
HAMILTON=DUDX+DUDY
!HAMILTON=-COS(DUDX+DUDY+1.0)
END FUNCTION HAMILTON
!--------------------------------------------------------------

SUBROUTINE INITIAL(NODE,VNODE)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
INTEGER I,VNODE
DO I=1,VNODE
  NODE(I).VALUE=U0_VALUE(NODE(I).X,NODE(I).Y)
END DO
END SUBROUTINE INITIAL

SUBROUTINE CALCULATE_NODE_CONST(GRID,NODE,KNODE)
IMPLICIT NONE
EXTERNAL FIND_NUMBER
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KNODE,FIND_NUMBER,I,J,K,N1,N2,S
DOUBLE PRECISION X,Y,R,A(6)
!對頂點周圍的三角形重新進行排列,按逆時針順序
DO I=1,KNODE
  DO J=1,NODE(I).NUM
    K=NODE(I).TRIMATRIX(J)
	N1=GRID(K).NODE(1)
	N2=GRID(K).NODE(2)
	IF(I==N1)THEN
	  N1=GRID(K).NODE(3)
	ELSE IF(I==N2)THEN
	  N2=GRID(K).NODE(3)
	END IF
	X=(NODE(N1).X+NODE(N2).X)*0.5-NODE(I).X
	Y=(NODE(N1).Y+NODE(N2).Y)*0.5-NODE(I).Y
	R=SQRT(X**2.0+Y**2.0)
	IF(Y>0)THEN
	  A(J)=ACOS(X/R)
	ELSE
	  A(J)=2*PI-ACOS(X/R)
	END IF
  END DO
  DO J=1,NODE(I).NUM-1
    R=A(J)
	K=J
	DO N1=J+1,NODE(I).NUM
	  IF(A(N1)<R)THEN
	    R=A(N1)
		K=N1
	  END IF
	END DO
	IF(K/=J)THEN
	  R=A(J)
	  A(J)=A(K)
	  A(K)=R
	  N1=NODE(I).TRIMATRIX(J)
	  NODE(I).TRIMATRIX(J)=NODE(I).TRIMATRIX(K)
	  NODE(I).TRIMATRIX(K)=N1
	END IF
  END DO
END DO

DO I=1,KNODE
  DO J=1,NODE(I).NUM
    K=NODE(I).TRIMATRIX(J)
	N1=GRID(K).NODE(1)
	N2=GRID(K).NODE(2)
	IF(I==N1)THEN
	  N1=GRID(K).NODE(3)
	ELSE IF(I==N2)THEN
	  N2=GRID(K).NODE(3)
	END IF
	X=(NODE(N1).X-NODE(I).X)**2.0+(NODE(N1).Y-NODE(I).Y)**2.0
	Y=(NODE(N2).X-NODE(I).X)**2.0+(NODE(N2).Y-NODE(I).Y)**2.0
    R=(NODE(N1).X-NODE(N2).X)**2.0+(NODE(N1).Y-NODE(N2).Y)**2.0
	NODE(I).INTERANGLE(J)=ACOS((X+Y-R)/(2.0*SQRT(X*Y)))
  END DO
END DO

DO I=1,KNODE
  DO J=1,NODE(I).NUM
    N1=NODE(I).TRIMATRIX(J)
	IF(J==NODE(I).NUM)THEN
	  N2=NODE(I).TRIMATRIX(1)
	ELSE
	  N2=NODE(I).TRIMATRIX(J+1)
	END IF
	K=0
	S=0
	DO WHILE((K<=3).AND.(S==0))
	  K=K+1
	  IF(GRID(N1).NODE(K)/=I)THEN
	    S=FIND_NUMBER(GRID(N2).NODE,3,GRID(N1).NODE(K))
	  ELSE
	    S=0
	  END IF
	END DO
	K=GRID(N1).NODE(K)
	X=NODE(K).X-NODE(I).X
	Y=NODE(K).Y-NODE(I).Y
	R=SQRT(X**2.0+Y**2.0)
	NODE(I).NL(1,J)=X/R
	NODE(I).NL(2,J)=Y/R
  END DO
END DO

DO I=1,KNODE
  DO J=1,NODE(I).NUM
    IF(J==6)THEN
	  K=1
	ELSE
	  K=J+1
	END IF
	NODE(I).BETA(J)=TAN(NODE(I).INTERANGLE(J)/2.0)+&
	  TAN(NODE(I).INTERANGLE(K)/2.0)
  END DO
END DO 
END SUBROUTINE CALCULATE_NODE_CONST

SUBROUTINE CALCULATE_MODULE(GRID,NODE,KGRID,VGRID)
IMPLICIT NONE
EXTERNAL FIND_NUMBER
EXTERNAL DECOMP
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER DECOMP,FIND_NUMBER,KGRID,VGRID,S,I,J,K,M,N1,N2
INTEGER BT(KGRID,3),BA(KGRID,3),B(6)
DOUBLE PRECISION X,Y,A(6,6)
!找出與每個三角形相鄰的三角形編號和其另一個接點編號
DO I=1,KGRID
  DO J=1,3
    N1=GRID(I).NODE(J)
	IF(J==3)THEN
	  N2=GRID(I).NODE(1)
	ELSE
	  N2=GRID(I).NODE(J+1)
	END IF
	K=0
	S=0
	DO WHILE((K<=VGRID).AND.(S==0))
	  K=K+1
	  IF(K/=I)THEN
	    S=FIND_NUMBER(GRID(K).NODE,3,N1)*FIND_NUMBER(GRID(K).NODE,3,N2)
	  ELSE 
	    S=0
	  END IF
	END DO
	BT(I,J)=K
	M=1
	DO WHILE((M<=3).AND.&
	  ((GRID(K).NODE(M)==N1).OR.(GRID(K).NODE(M)==N2)))
	  M=M+1
	END DO
	BA(I,J)=GRID(K).NODE(M)
  END DO
END DO

!求出插值子模塊
DO I=1,KGRID
  GRID(I).MONUM=0
  DO K=1,4
    IF(K==1)THEN
	  M=I
	ELSE
	  M=BT(I,K-1)
	END IF
    IF(M<=KGRID)THEN
	  DO J=1,3
	    B(J)=GRID(M).NODE(J)
		B(J+3)=BA(M,J)
	  END DO
!A1+A2*X+A3*Y+A4*X^2+A5*X*Y+A6*Y^2
	  DO J=1,6
	    X=NODE(B(J)).X
		Y=NODE(B(J)).Y
		A(J,1)=1.0
	    A(J,2)=X
	    A(J,3)=Y
	    A(J,4)=X*X
	    A(J,5)=X*Y
	    A(J,6)=Y*Y
	  END DO
	  S=DECOMP(A,B,6)
	  IF(S==1)THEN
	    GRID(I).MONUM=GRID(I).MONUM+1
	    GRID(I).MOCONST(GRID(I).MONUM,1:6,1:6)=A
	    GRID(I).MONODE(GRID(I).MONUM,1:6)=B
	  END IF
	END IF
  END DO
END DO
END SUBROUTINE CALCULATE_MODULE

SUBROUTINE DATA_OUTPUT(NODE,GRID,KGRID,KNODE)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KGRID,KNODE,I
OPEN(1,FILE='D:\TANG\MATLAB程序\weno_two\node.dat')
DO I=1,KNODE
  WRITE(1,*) NODE(I).X,NODE(I).Y,NODE(I).VALUE
END DO
CLOSE(1)
OPEN(1,FILE='D:\TANG\MATLAB程序\weno_two\grid.dat')
DO I=1,KGRID
  WRITE(1,*) GRID(I).NODE(1),GRID(I).NODE(2),GRID(I).NODE(3)
END DO
CLOSE(1)
END SUBROUTINE DATA_OUTPUT

SUBROUTINE ENCODE_TRIANGLE(GRID,NODE,KGRID)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KGRID,I,J,K
DO I=1,MN
  NODE(I).NUM=0
END DO
DO I=1,KGRID
  DO J=1,3
    K=GRID(I).NODE(J)
    NODE(K).NUM=NODE(K).NUM+1
	NODE(K).TRIMATRIX(NODE(K).NUM)=I
  END DO
END DO
END SUBROUTINE ENCODE_TRIANGLE

SUBROUTINE CONDUCT_BOUNDARY(GRID,NODE,KNODE,KGRID,VNODE,VGRID)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KNODE,KGRID,VNODE,VGRID
INTEGER I
DOUBLE PRECISION X,Y,DETX,DETY,O
VNODE=KNODE
VGRID=KGRID
DETX=XR-XL
DETY=YT-YB
O=0.0
DO I=1,KNODE
  X=NODE(I).X
  Y=NODE(I).Y
  IF((X==XL).AND.(Y==YB))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,DETX,O)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,DETY)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,DETX,DETY)
  END IF
  IF((X==XL).AND.(Y==YT))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,DETX,O)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,-DETY)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,DETX,-DETY)
  END IF
  IF((X==XL).AND.(Y/=YB).AND.(Y/=YT))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,DETX,O)
  END IF
  IF((X==XR).AND.(Y==YB))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,-DETX,O)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,DETY)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,-DETX,DETY)
  END IF
  IF((X==XR).AND.(Y==YT))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,-DETX,O)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,-DETY)
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,-DETX,-DETY)
  END IF
  IF((X==XR).AND.(Y/=YB).AND.(Y/=YT))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,-DETX,O)
  END IF
  IF((Y==YB).AND.(X/=XL).AND.(X/=XR))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,DETY)
  END IF
  IF((Y==YT).AND.(X/=XL).AND.(X/=XR))THEN
	CALL ADD_NODE_GRID(GRID,NODE,I,VNODE,VGRID,O,-DETY)
  END IF
END DO
PRINT *,'VNODE=',VNODE
PRINT *,'VGRID=',VGRID
END SUBROUTINE CONDUCT_BOUNDARY

!虛點結構中num位置存儲的是其對應實點的編號
!虛網格結構MONUM中存儲的是其對應的實網格的編號
SUBROUTINE ADD_NODE_GRID(GRID,NODE,ST,VNODE,VGRID,DX,DY)
IMPLICIT NONE
EXTERNAL FIND_NUMBER
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER ST,VNODE,VGRID,FIND_NUMBER
INTEGER N,J,K,A(3)
DOUBLE PRECISION X,Y,DX,DY
DO J=1,NODE(ST).NUM
  A=GRID(NODE(ST).TRIMATRIX(J)).NODE
  DO K=1,3
	X=NODE(A(K)).X+DX
	Y=NODE(A(K)).Y+DY
	N=1
	DO WHILE((N<=VNODE).AND.((NODE(N).X/=X).OR.(NODE(N).Y/=Y)))
	  N=N+1
	END DO
	IF(N>VNODE)THEN
	  VNODE=VNODE+1
	  NODE(VNODE).X=X
	  NODE(VNODE).Y=Y
	  NODE(VNODE).NUM=A(K)
	  A(K)=VNODE
	ELSE
	  A(K)=N
	END IF
  END DO
  K=0
  N=0
  DO WHILE((K<=VGRID).AND.(N==0))
	K=K+1
    N=FIND_NUMBER(GRID(K).NODE,3,A(1))*&
	  FIND_NUMBER(GRID(K).NODE,3,A(2))*&
	  FIND_NUMBER(GRID(K).NODE,3,A(3))
  END DO
  IF(K>VGRID)THEN
	VGRID=VGRID+1
    GRID(VGRID).NODE=A
	GRID(VGRID).MONUM=NODE(ST).TRIMATRIX(J)
  END IF
END DO
END SUBROUTINE ADD_NODE_GRID

!生成網格,對網格和接點進行編號
SUBROUTINE GENERATE_GRID(GRID,NODE,KNODE,KGRID)
IMPLICIT NONE
TYPE(VERTEX_TYPE),DIMENSION(MN):: NODE
TYPE(TRIANGLE_TYPE),DIMENSION(MG):: GRID
INTEGER KNODE,KGRID
INTEGER I,J,K1,K2,SPLITIME,NUM,A(3)
DOUBLE PRECISION X,Y
PRINT *,'TIMES OF SPLIT='
READ(*,*) SPLITIME
NODE(1).X=XL
NODE(1).Y=YT
NODE(2).X=XL
NODE(2).Y=YB
NODE(3).X=XR
NODE(3).Y=YB
NODE(4).X=XR
NODE(4).Y=YT
KNODE=4
GRID(1).NODE=(/1,2,3/)
GRID(2).NODE=(/1,3,4/)
KGRID=2
DO I=1,SPLITIME
  NUM=KGRID
  DO J=1,NUM
    A=GRID(J).NODE
    DO K1=1,3
	  IF(K1/=3)THEN
	    K2=K1+1
	  ELSE
	    K2=1
	  END IF
	  X=(NODE(A(K1)).X+NODE(A(K2)).X)*0.5
	  Y=(NODE(A(K1)).Y+NODE(A(K2)).Y)*0.5
	  DO K2=1,KNODE
	    IF((NODE(K2).X==X).AND.(NODE(K2).Y==Y))THEN
		  GOTO 10
		END IF
	  END DO
10    CONTINUE
      IF(K2>KNODE)THEN
	    KNODE=KNODE+1
		NODE(KNODE).X=X
		NODE(KNODE).Y=Y
		GRID(J).NODE(K1)=KNODE
	  ELSE
	    GRID(J).NODE(K1)=K2
	  END IF
	END DO
	GRID(KGRID+1).NODE=(/A(1),GRID(J).NODE(1),GRID(J).NODE(3)/)
	GRID(KGRID+2).NODE=(/A(2),GRID(J).NODE(1),GRID(J).NODE(2)/)
	GRID(KGRID+3).NODE=(/A(3),GRID(J).NODE(2),GRID(J).NODE(3)/)
	KGRID=KGRID+3
  END DO
END DO
PRINT *,'KGRID=',KGRID
PRINT *,'KNODE=',KNODE
END SUBROUTINE GENERATE_GRID
END PROGRAM MAIN