C     *******************************************
C     **  TWO-DIMENSIONAL GRID GENERATION      **
C     **  BY POISSON EQUATION                  **
C     **                                       **
C     **  CONTROL FUNCTION USED BY TOMCAT      **
C     **                                       **
C     **  NOTE THAT DXI=1 & DETA=1             **
C     **                                       **
C     **  SOLVED BY SUCCESSIVE OVER RELAXATION **
C     **                                       **
C     **                                       **
C     **  1997/10/20  ‹Þ’æ  Œ¤‹†‰ï‚ÌŠF—l       **
C     *******************************************
C
C
         CALL INIT
C
         CALL POIST
C
         CALL GRID
C
         CALL OUTPUT
C
C
      STOP
      END
C
C     ***********************************
C     ** SUBROUTINE READING GRID DATA  **
C     ***********************************
C
      SUBROUTINE INIT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /DIV/ NX,NY
      COMMON /COORD/ X(301,101),Y(301,101)
C
C
      OPEN(1,FILE='A.DAT')
C
      READ(1,600) NX,NY
  600 FORMAT(2I10)
C
      READ(1,610)((X(I,J),Y(I,J),I=1,NX),J=1,NY)
  610 FORMAT(2F15.10)
C
      CLOSE(1)
C
C
      RETURN
      END
C
C
C     ***************************
C     **   SUBROUTINE OUTPUT   **
C     ***************************
C
      SUBROUTINE OUTPUT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /DIV/ NX,NY
      COMMON /COORD/ X(301,101),Y(301,101)
C
C
      OPEN(2,FILE='B.DAT')
C
      WRITE(2,600) NX,NY
  600 FORMAT(2I10)
C
      WRITE(2,601)((X(I,J),Y(I,J),I=1,NX),J=1,NY)
  601 FORMAT(2F15.10)
C
      CLOSE(2)
C
C
      RETURN
      END
C
C     ************************************
C     ** RIGHT HAND SIDE OF POISSON EQ. **
C     ************************************
C
      SUBROUTINE POIST
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /DIV/ NX,NY
      COMMON /COORD/ X(301,101),Y(301,101)
      COMMON /RIGHT/ P(301,101),Q(301,101)
C
C     *** GRID CONTROL TERM ***
C
        XILN=1.D0
        ETLN1=1.D0
        ETLN2=DFLOAT(NY)
C
        WRITE(*,*)' ALPHA='
        READ(*,*) ALPHA
C
        WRITE(*,*)' BETA='
        READ(*,*) BETA
C
C
      DO 100 J=1,NY
      DO 100 I=1,NX
C
       P(I,J)=DSIGN(1.D0,(DFLOAT(I)-XILN))*ALPHA
     &         *DEXP(-BETA*DABS(DFLOAT(I)-XILN))*0.0D0
C
       Q(I,J)=DSIGN(1.D0,(DFLOAT(J)-ETLN1))*ALPHA
     &         *DEXP(-BETA*DABS(DFLOAT(J)-ETLN1))
     &       +DSIGN(1.D0,(DFLOAT(J)-ETLN2))*ALPHA
     &         *DEXP(-BETA*DABS(DFLOAT(J)-ETLN2))
C
  100 CONTINUE
C
C
      RETURN
      END
C
C     **********************************
C     **  SUBROUTINE GRID GENERATION  **
C     **********************************
C
      SUBROUTINE GRID
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /DIV/ NX,NY
      COMMON /COORD/ X(301,101),Y(301,101)
      COMMON /RIGHT/ P(301,101),Q(301,101)
C
C
      WRITE(*,*)'OMEGA='
      READ(*,*) OMEGA
C
        EPS =1.0D-4
        ITER=500
C
C
C     ***** SOR ITERATION LOOP *****
C
      DO 300 LOOP=1,ITER
C
C
         ERR=0.0D0
C
C     ***** BRANCH CUT CONDITION *****
C
      DO 110 J=1,NY
C
         X(1,J)=X(NX,J)
  110    Y(1,J)=Y(NX,J)
C
C
C     ***** CAL. INTERNAL POINT *****
C
      DO 100 J=2,NY-1
C
         JP1=J+1
         JM1=J-1
C
      DO 100 I=2,NX
C
         IP1=I+1
         IM1=I-1
C
      IF(I.EQ.NX) IP1=2
C
         XO=X(I,J)
         YO=Y(I,J)
C
C     *** METRIC ***
C
         XX=(X(IP1,J)-X(IM1,J))*.5D0
         YX=(Y(IP1,J)-Y(IM1,J))*.5D0
         XE=(X(I,JP1)-X(I,JM1))*.5D0
         YE=(Y(I,JP1)-Y(I,JM1))*.5D0
C
C     *** JACOBIAN ***
C
         YAC=XX*YE-XE*YX
C
C     *** COEFFICCIENT ***
C
           ALPH=XE*XE+YE*YE
           BETA=XX*XE+YX*YE
           GAM =XX*XX+YX*YX
C
C     *** SOR COEFFICIENT ***
C
           COEF=2.D0*(ALPH+GAM)
C
C     *** FINITE DIFFERENCE EQUATION ***
C
      DX=ALPH*(X(IP1,J)+X(IM1,J))
     &  -BETA*(X(IP1,JP1)-X(IP1,JM1)-X(IM1,JP1)+X(IM1,JM1))*.5D0
     &   +GAM*(X(I,JM1)+X(I,JP1))
C
      DY=ALPH*(Y(IP1,J)+Y(IM1,J))
     &  -BETA*(Y(IP1,JP1)-Y(IP1,JM1)-Y(IM1,JP1)+Y(IM1,JM1))*.5D0
     &   +GAM*(Y(I,JM1)+Y(I,JP1))
C
C     *** CONTROLL TERM ***
C
        PQX=YAC*YAC*(XX*P(I,J)+XE*Q(I,J))
        PQY=YAC*YAC*(YX*P(I,J)+YE*Q(I,J))
C
C     *** RESIDUAL ***
C
         RX=(DX-PQX)/COEF-XO
         RY=(DY-PQY)/COEF-YO
C
C     *** RELAXATION ***
C
         X(I,J)=X(I,J)+OMEGA*RX
         Y(I,J)=Y(I,J)+OMEGA*RY
C
         ERR=ERR+RX*RX+RY*RY
C
  100 CONTINUE
C
         ERR=DSQRT(ERR/DFLOAT((NX-2)*(NY-2)))
C
      WRITE(*,600) LOOP,ERR
  600 FORMAT(I5,E10.4)
C
      IF(ERR.LT.EPS) GOTO 99
C
  300 CONTINUE
C
   99 CONTINUE
C
C
      RETURN
      END