C 2-D finite element method for ground water transport simulation 
 
      INTEGER  E,P,PP
      REAL  NN
      DIMENSION  X(20),Y(20),COLD(20),CNEW(20),G(20,20),
   J  ST(20),U(20,20),B(20)
      DIMENSION  XSI(4)，ETA(4),NN(4),DNDX(4),DNDY(4),MP(4)
      DIMENSION  IJMN(9,4),KODE(20)
 
C Input date

      DATA  NE, NP/9, 20/
      DATA  XST/-0.57735, 0.57735, 0.57735,- 0.57735/
      DATA  ETA/-0.57735, 0.57735, 0.57735,- 0.57735/
      DATA  Y/10.,0., 10.,0., 10.,0., 10.,0., 10.,0., 10.,0., 10.,0., 10.,0.,
   J    10.,0., 10.,0./
      DATA  X/2*0., 2*10., 2*20., 2*30., 2*40., 2*50., 2*60., 2*70.,
   J     2*80., 2*90./
      DATA  CNEW/2*0.0/
      DATA  VX ,  DXX ,  DYY/0.01,1.0,0.1/
      DATA  KODE/2*1,16*0.2*1/
      DATA  IJMN/2,4,6,8,10,12,14,16,18,
      J          4,6,8,10,12,14,16,18,20,
      J          3,5,7,9,11,13,15,17,19,
      J          1,3,5,7,9,11,13,15,17
 
   C Input  initial  and  boundaral  data
 
     DO  P=1,NP
          COLD(P)=CNEW(P)
     DO  PP=1,NP
          G(P,PP)=0.0
          ST(P,PP)=0.0
          U（p,pp）=0.0
     END  DO
         END  DO
 
   C Caiculating  inatrix  element
 
     DO  100  E=1,  NE
 
     I=IJMN(E,1)
     J=IJMN(E,2)
   M=IJMN(E,3)
   N=IJMN(E,4)
 
   MP(1)=I
   MP(2)=J
   MP(3)=M
   MP(4)=N
 
   AA=ABS(X(J)-X(I))/2.
   BB=ABS(Y(N)-Y(I))/2.
 
   DO 40 P=1,4
 
   L=MP(P)
 
   DO 30 PP=1,4
 
      NN(1)=0.25*(1-XSJ(PP))*(1-ETA(PP))
      NN(1)=0.25*(1-XSJ(PP))*(1-ETA(PP))
      NN(1)=0.25*(1-XSJ(PP))*(1-ETA(PP))
      NN(1)=0.25*(1-XSJ(PP))*(1-ETA(PP))
 
      DNDX(1)=-0.25*(1-ETA(PP))/AA
      DNDX(2)=0.25*(1-ETA(PP))/AA
      DNDX(3)=0.25*(1+ETA(PP))/AA
      DNDX(4)=-0.25*(1+ETA(PP))/AA
 
      DNDX(1)=-0.25*(1-XSI(PP))/BB
      DNDX(2)=-0.25*(1+XSI(PP))/BB
      DNDX(3)=0.25*(1+XSI(PP))/BB
      DNDX(4)=0.25*(1-XSI(PP))/BB
  
      G(L,I)+(DXX*DNDX(1)*DNDX(P)+DYY*DNDY(1)
J      *DNDY(P)*AA*BB
      G(L,I)+(DXX*DNDX(2)*DNDX(P)+DYY*DNDY(2)
J      *DNDY(P)*AA*BB
      G(L,I)+(DXX*DNDX(3)*DNDX(P)+DYY*DNDY(3)
J      *DNDY(P)*AA*BB
      G(L,I)+(DXX*DNDX(4)*DNDX(P)+DYY*DNDY(4)
J      *DNDY(P)*AA*BB
 
   U(L,I)=U(L,I)+VX*DNDX(1)*NN(P)*AA*BB
   U(L,J)=U(L,J)+VX*DNDX(2)*NN(P)*AA*BB
   U(L,M)=U(L,M)+VX*DNDX(3)*NN(P)*AA*BB
U(L,N)=U(L,N)+VX*DNDX(4)*NN(P)*AA*BB
 
 ST(L,I)=ST(L,I)+NN(1)*NN(P)*AA*BB
 ST(L,J)=ST(L,J)+NN(2)*NN(P)*AA*BB
 ST(L,M)=ST(L,M)+NN(3)*NN(P)*AA*BB
 ST(L,N)=ST(L,N)+NN(4)*NN(P)*AA*BB
 
30   CONTINUE
40   CONTINUE
100   CONTINUE
 
 C    Solving and printing
 
      DT=5.0
      KOUNT=1
      KPRINT=2
      T=DT
 
      WRINT(*,120)
 120  FORMAT(/,32X,’CONCENTRATION’,28X’,TIME’,/)
 
     DO  500  NSTEP=1.100
     DO  150  L=1,NP
     B(L)=0.0
     DO  150  P=1,NP
     B(L)=B(L)+ST(L,P)*cold（p）/DT
 150  CONTINUE
 200  ERMAX=0.0
    DO 400 L=1,NP
    IF(KODE(L).EQ.1)GO TO 400
    OLDVAL=CNEW(L)
    SUM=0.0
    DO  300  P=1,NP
    IF(P.EQ.L)GO  TO  300
    SUM=SUM+（G(L,P)+U(L,P)+ST(L,P)/DT）*CNEW(P)
 300  CONTINUE
       CNEW(L)=(-SUM+B(L))/(G(L,L)+U(L,L)+ST(L,L)/DT)
       ER=ABS(OLDVAL-CNEW(L))
       IF(ER.GT.ERMAX)ERMAX=ER
 400  CONTINUE
     IF(ERMNX.GT.0.001)GO  TO  200
     DO  450  L=1,NP
 450   COLD(L)=CNEW(L)
 
       IF(KOUNT.NE.KPRINT)JO  TO  490 
       WRITE(*,401) (CNEW(I),I=1,19,2),T
 401    FORMAT(1X,10F7.3,1F9.2)
 
       KOUNT=0
 490    T=T+DT
       KOUNT=KOUNT+1
 500    CONTINUE
 
       STOP
       END
 
 C  End  program