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),
             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.,
           10.,0., 10.,0./
      DATA  X/2*0., 2*10., 2*20., 2*30., 2*40., 2*50., 2*60., 2*70.,
        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 =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)=G(L,I)+(DXX*DNDX(1)*DNDX(P)+DYY*DNDY(1)
      J      *DNDY(P)*AA*BB
      G(L,I)=G(L,J)+(DXX*DNDX(2)*DNDX(P)+DYY*DNDY(2)
      J      *DNDY(P)*AA*BB
      G(L,I)=G(L,M)+(DXX*DNDX(3)*DNDX(P)+DYY*DNDY(3)
      J      *DNDY(P)*AA*BB
      G(L,N)=G(L,N)+(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