program diffa
C
C       THIS PROGRAM SOLVES THE LAPLACE EQUATION IMPLICITLY BY USING
C       SUCCESSIVE OVER RELAXATION (SOR). THE PROGRAM IS GENERAL IN 
C       THAT ANY RECTANGULAR GRID CAN BE ESTABLISHED. ANY BOUNDARY
C       CAN BE ESTABLISHED ON THE GRID. THE VALUES AT THE BOUNDARIES
C       MUST BE GIVEN TOGETHER WITH GUESSES FOR INTERIOR POINT VALUES.
C       OTHER INPUT PARAMETERS ARE EXPLAINED THRU OUT THE PROGRAM.
C
      dimension u(2601),uu(51,51)
c
      COMMON ILHS(50), IRHS(50), A(2601), B(2601), C(2601)
     $,D(2601), E(2601), G(2601), H(2601)
C
C       OPEN OUTPUT FILES
c
      OPEN(6,FILE='diffa.out')
      open(7,file='timesa.dat')
      open(8,file='concena.ascii')
c
C  N = THE NUMBER OF INTERVALS IN A UNIT LENGTH.
      N=25
C  LVERT= THE NUMBER OF UNIT LENGTHS ON THE VERTICAL SIDE OF THE GRID.
      LVERT=2
C  LHORIZ = THE NUMBER OF UNIT LENGTHS ON THE HORIZONTAL SIDE OF THE GRID.
      LHORZ=2
C  W = CONVERGENCE FACTOR FOR THE SUCCESSIVE OVER-RELAXATION METHOD.
      W=1.85
C  SIGMA = THE CONSTANT IN THE LAPLACE EQUATION.
      SIGMA=1.6/(10**5)
C  THETA = A PARAMETER SUCH THAT 1/2 < THETA < 1.
      THETA=0.75
C  PAR = (A CHANGE IN TIME)/((1/N)**2)
      PAR=15875.0
C  TT = THE TOTAL NUMBER OF TIME STEPS
      tt=1000
C  ITMAX = THE LARGEST NUMBER OF ITERATIONS FOR CONVERGENCE.
      ITMAX=500
C  ERRSOR = THE SOR ERROR TO EXPECT IN THE FINAL SOLUTIONS.
      ERRSOR=0.001
C
      DH=1./N
      DT=PAR*(DH**2)
      tend=tt*dt
      NM=LVERT*N-1
      IC=LHORZ*N+1
      M=IC*(NM+2)
      icp1=ic+1
      WRITE(6,19)DH,DT,tend,NM,IC,M
   19 FORMAT(3E10.3,3I10)
C
C       ASSIGNING INTERIOR POINT VALUES
      DO 9 jj=1,M
    9 U(jj)=0.25
C
C       LOCATING BOUNDARY AND ASSIGNING VALUES.
c
      DO 5 K=1, NM
      ILHS(K)=(K-1)*ic+ic+2
      IRHS(K)=(K-1)*ic+2*ic-1
    5 CONTINUE
c
c        PRESCRIBE THE INITIAL CONDITIONS
C
      TD=0.0
C
      ibd=(ic-1)/5
      ibm=(ic-1)/2+1
      iba=ibm-ibd
      ibb=ibm+ibd
      write(6,121)ibd,ibm,iba,ibb
  121 format(4(1x,i8))
      do 88 ib=iba,ibb
c     do 88 ib=2,50
      u(ib)=1.20
   88 continue
C
c        ESTABLISH COEFF.S FOR INTERATION
C
      CALL ESTAB(N,LHORZ,LVERT,SIGMA,THETA,DT,DH)
c
c        BEGIN TIME LOOP
C
      itime=1
   28 continue
c
C       WRITE CONCENTRATIONS FOR EACH TIME STEP
C
      icount=1
      ii=1
      jj=1
      DO 896 iZ=1,m
      if(icount.eq.icp1)ii=1
      if(icount.eq.icp1)jj=jj+1
      if(icount.eq.icp1)icount=1
      uu(ii,jj)=u(iz)
      icount=icount+1
      ii=ii+1
  896 continue
      write(8,95)((uu(ia,jb),ia=1,ic),jb=1,ic)
      write(6,203)itime,td,uu(ibm,12)
      write(7,204)itime,td,uu(ibm,12)
c
      itime=itime+1
      TD=TD+DT
      CALL SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM,ITER)
c
      IF(TD.LE.TEND) GO TO 28
c
c
C        WRITE FINAL RESULTS
C
      WRITE(6,100)N,LHORZ,LVERT,W,SIGMA,THETA,DT,DH,ERRSOR,ITMAX,DNRM
     $,ITER,TEND
      WRITE(6,202)TD
c
    1 FORMAT(20X,I3,24X,I3)
    2 FORMAT(6X,F10.5)
   95 format(7(1x,e10.3))
  100 FORMAT(//,' NO. OF INTERVALS PER LENGTH=',I3,/,1X,
     1'NO. OF LENGTHS ON THE HORIZONTAL=',I4,/,1X,
     1'NO. OF LENGTHS ON THE VERTICAL=',I4,/,1x,
     1'CONVERGE FACTOR=',F4.2,/,1X,'SIGMA=',e11.4,/,1X,'THETA='
     1,F4.2,/,1X,'TIME INTERVAL=',E15.8,/,1X,'STEP-SIZE IN X AND Y CORDI
     1NATES=',E15.8,/,1X,'ERRSOR=',E15.8,/,1X,'ITMAX=',I4,/,1X,'MEASURE
     1OF LARGEST ERROR BETWEEN LAST TWO ITERATES',E15.8,/,1X,'NUMBER OF
     1ITERATIONS DONE',I4,/,1X,'FINAL TIME',e11.4,///,1X,
     1'FINAL SOLUTIONS ARE STORED IN FILE concena.ascii....',//)
c
  202 FORMAT(1X,'TIME=',e11.4,//)
  203 format(1x,' I-th TIME STEP=',i4,' TIME=',e11.4,
     1' WT% @ 0.5mm =',e11.4)
  204 format(1x,i4,1x,e11.4,1x,e11.4)
  300 FORMAT(48X,E11.4)
  999 STOP
      END
c
      SUBROUTINE SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM,ITER)
      COMMON ILHS(50),IRHS(50),A(2601),B(2601),C(2601)
     $,D(2601),E(2601),G(2601),H(2601)
      DIMENSION U(2700),F(2700)
      NM=LVERT*N-1
      IC=LHORZ*N+1
      IT=LHORZ*N+3
      DO 4 ITER=1,ITMAX
      DO 3 J=1,NM
      IB=ILHS(J)
      ID=IRHS(J)
      DO 3 K=IB,ID
      F(K)=G(K)*U(K)+H(K)*(U(K-IC)+U(K-1)+U(K+1)+U(K+IC))
      T=U(K)+(W/C(K))*(F(K)-(A(K)*U(K-IC)+B(K)*U(K-1)+C(K)*U(K)
     $+D(K)*U(K+1)+E(K)*U(K+IC)))
      IF(K.EQ.IT) GO TO 1
      GO TO 2
    1 DNRM=ABS(T-U(K))
      GO TO 3
    2 DNRM=AMAX1(DNRM,ABS(T-U(K)))
    3 U(K)=T
      IF(DNRM.LE.ERRSOR) GO TO 5
    4 CONTINUE
      WRITE(6,6)
    5 RETURN
    6 FORMAT('  DID NOT CONVERGE IN SLEEP')
      END
c
      SUBROUTINE ESTAB(N,LHORZ,LVERT,SIGMA,THETA,DT,DH)
      COMMON ILHS(50),IRHS(50),A(2601),B(2601),C(2601)
     $,D(2601),E(2601),G(2601),H(2601)
      V=SIGMA*THETA*DT/DH**2
      X=1+4*V
      Y=SIGMA*(1-THETA)*DT/DH**2
      Z=1-4*Y
      NM=LVERT*N-1
      IC=LHORZ*N+1
      M=IC*(NM+2)
      DO 1 J=1,M
      A(J)=-V
      B(J)=-V
      C(J)=X
      D(J)=-V
      E(J)=-V
      G(J)=Z
    1 H(J)=Y
      RETURN
      END