Fortran90 Diffusion Program

     PROGRAM DIFF2
     IMPLICIT NONE
!
!       THIS PROGRAM USES THE LAPLACE EQUATION TO SOLVE FOR
!       DIFFUSION.  THE PROGRAM USES A TECHNIQUE CALLED SUCCESSIVE 
!       OVER RELAXATION (SOR).  THE PROGRAM CAN USE ANY RECTANGULAR 
!       GRID WITH A BOUNDARY ON THAT GRID.  THE VALUES AT THE 
!       BOUNDARY MUST BE GIVEN ALONG WITH GUESSES FOR INTERIOR
!       POINT VALUES.  OTHER INPUT PARAMETERS ARE EXPLAINED
!       THROUGHOUT THE PROGRAM.
!
!       UU = ARRAY IN (X,Y) FORM
       REAL, DIMENSION(51,51) :: UU
!
!       ILHS = INTEGER LEFT HAND SIDE - LEFT HAND SIDE TO INTEGRATE
!       IRHS = INTEGER RIGHT HAND SIDE - RIGHT HAND SIDE TO INTEGRATE
!       A, B, C, D, E, G, H = WEIGHTED AVERAGES
!       U = ARRAY OF ALL POINTS
      REAL, DIMENSION(2601) :: A, B, C, D, E, F,G, H, U
      INTEGER, DIMENSION(49) :: ILHS, IRHS
!
      INTEGER :: N, LVERT, LHORZ, TT, ITMAX, ITER, NM, IC, M, ICP1&
                      &, ITIME, K,IBM, IB, ICOUNT, II, JJ, IZ, IA, JB
      REAL :: W, SIGMA, THETA, PAR, ERRSOR, DH, DT, TEND, TD, DNRM
!
!       OPEN OUTPUT FILES
!
      OPEN(6,FILE='DIFF2.OUT')
      OPEN(7,FILE='TIMES2.DAT')
      OPEN(8,FILE='CONCEN2.ASCII')
!
!       N = THE NUMBER OF INTERVALS IN A UNIT LENGTH
      N=25
!       LVERT= THE NUMBER OF UNIT LENGTHS ON THE VERTICAL 
!       SIDE OF THE GRID - IN THIS CASE 1 MM LONG
      LVERT=2
!      LHORIZ = THE NUMBER OF UNIT LENGTHS ON THE HORIZONTAL
!                         SIDE OF THE GRID
      LHORZ=2
!       W = CONVERGENCE FACTOR FOR THE SUCCESSIVE 
!               OVER-RELAXATION METHOD
      W=1.85
!       SIGMA = THE CONSTANT IN THE LAPLACE EQUATION - 
!       IN THIS CASE THE DIFFUSION CONSTANT
       SIGMA=1.6/(10**5)
!       THETA = A PARAMETER SUCH THAT 1/2 < THETA < 1 - 
!       SO THAT THE VALUES WILL ALWAYS CONVERGE
      THETA=0.75
!       PAR = (A CHANGE IN TIME)/((1/N)**2)
      PAR=15875.0
!       TT = THE TOTAL NUMBER OF TIME STEPS
      TT=1000
!       ITMAX = THE LARGEST NUMBER OF ITERATIONS FOR CONVERGENCE.
      ITMAX=500
!       ERRSOR = THE SOR ERROR TO EXPECT IN THE FINAL SOLUTIONS.
      ERRSOR=0.001
!
!       DELTA LENGTH

     DH=1./N
!       DELTA TIME
      DT=PAR*(DH**2)
!       TOTAL TIME
      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)
!
!       ASSIGNING INTERIOR POINT VALUES
      DO JJ=1,M
      U(JJ)=0.25
      END DO
!
!       LOCATING BOUNDARY AND ASSIGNING VALUES.
!
      DO  K=1, NM
      ILHS(K)=(K-1)*IC+IC+2
      IRHS(K)=(K-1)*IC+2*IC-1
      END DO
!
      TD=0.0
      IBM=(IC-1)/2+1
      DO  IB=2,50
      U(IB)=1.20
      END DO
!
!      ESTABLISH COEFF.S FOR INTERATION
!
      CALL ESTAB(N, LHORZ, LVERT ,SIGMA,THETA, DT, DH&
                            &, A, B, C, D, E, F,G, H)
!
!       BEGIN TIME LOOP
!
      ITIME=1
!
!       WRITE CONCENTRATIONS FOR EACH TIME STEP
!                      PUT U INTO ARRAY OF (X,Y)
 28 ICOUNT=1
      II=1
      JJ=1
      DO  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
      END DO
      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)
!
      ITIME=ITIME+1
      TD=TD+DT
      CALL SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM&
                               &,ITER,ILHS,IRHS,A,B,C,D,E,F,G,H)
!
      IF(TD.GE.TEND) GO TO 27
      GO TO 28
!
!       WRITE FINAL RESULTS
!
 27 WRITE(6,96)N,LHORZ,LVERT,W,SIGMA,THETA,DT,DH,ERRSOR,ITMAX&
                      &,DNRM,ITER,TEND
      WRITE(6,202)TD
!
  95 FORMAT(7(TR1,E10.3))
  96 FORMAT(//,' NO. OF INTERVALS PER LENGTH=',I3,/,TR1,'NO. OF&
     &LENGTHS ON THE HORIZONTAL=',I4,/,TR1,'NO. OF LENGTHS ON THE&
     &VERTICAL=',I4,/,TR1,'CONVERGE FACTOR=',F4.2,/,TR1,'SIGMA=',e11.4&
     &,/,TR1,'THETA=',F4.2,/,TR1,'TIME INTERVAL=',E15.8,/,TR1,'STEP-SIZE IN&
     &X AND Y CORDINATES=',E15.8,/,TR1,'ERRSOR=',E15.8,/,TR1,'ITMAX=',I4&
     &,/,TR1,'MEASURE OF LARGEST ERROR BETWEEN LAST TWO ITERATES'&
     &,E15.8,/,TR1,'NUMBER OF ITERATIONS DONE',I4,/,TR1,'FINAL TIME',e11.4&
     &,///,TR1,'FINAL SOLUTIONS ARE STORED IN FILE    CONCEN.ASCII ....',//)
!
202 FORMAT(TR1,'TIME=',E11.4,//)
203 FORMAT(TR1,' I-th TIME STEP=',I4,' TIME=',E11.4,' WT% @ 0.5mm =',E11.4)
204 FORMAT(TR1,I4,TR1,E11.4,TR1,E11.4)
       STOP

CONTAINS

      SUBROUTINE SLEEP(N,LHORZ,LVERT,W,U,ERRSOR,ITMAX,DNRM&
                                             &,ITER,ILHS,IRHS,A,B,C,D,E,F,G,H)
      IMPLICIT NONE
      REAL, DIMENSION(:), INTENT(IN OUT) :: A, B, C, D, E, F, G, H, U
      INTEGER, DIMENSION(:), INTENT(IN)  :: ILHS, IRHS
      INTEGER, INTENT(IN OUT)  :: N, LHORZ, LVERT, ITER, ITMAX
      INTEGER :: NM, IC, IT, J, K, IB, ID
      REAL, INTENT(IN OUT) :: W, ERRSOR, DNRM
      REAL :: T
      NM=LVERT*N-1
      IC=LHORZ*N+1
      IT=LHORZ*N+3
      DO ITER=1,ITMAX
      DO J=1,NM
      IB=ILHS(J)
      ID=IRHS(J)
      DO 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=MAX(DNRM,ABS(T-U(K)))
    3 U(K)=T
      END DO
      END DO
      IF(DNRM.LE.ERRSOR) GO TO 5
      END DO
      WRITE(6,6)
    5 RETURN
    6 FORMAT('  DID NOT CONVERGE IN SLEEP')
      END SUBROUTINE SLEEP

      SUBROUTINE ESTAB(NN,LHORZ,LVERT,SIGMA,THETA,DT,DH&
                                          &,A,B,C,D,E,F,G,H)
      IMPLICIT NONE
      REAL, DIMENSION(:), INTENT(IN OUT) :: A, B, C, D, E, F, G, H
      INTEGER, INTENT(IN) :: NN, LHORZ,LVERT
      INTEGER :: J, NM, IC, M
      REAL, INTENT(IN) :: SIGMA, THETA, DT, DH
      REAL :: V, X, Y, Z
      V=SIGMA*THETA*DT/DH**2
      X=1+4*V
      Y=SIGMA*(1-THETA)*DT/DH**2
      Z=1-4*Y
      NM=LVERT*NN-1
      IC=LHORZ*NN+1
      M=IC*(NM+2)
      DO J=1,M
      A(J)=-V
      B(J)=-V
      C(J)=X
      D(J)=-V
      E(J)=-V
      F(J)=0.0
      G(J)=Z
      H(J)=Y
      END DO
      RETURN
      END SUBROUTINE ESTAB

END PROGRAM DIFF2