tlm.for

tlm code

      WRITE(*,*) 'Start of calculation'

      DO 30010, T = 1, NUMDT

      WRITE(*,99100) 100.0*REAL(T-1)/REAL(NUMDT)
99100 FORMAT(1X,F6.3,' %')

*---- output

      DO 41010, N = 1, NUMNOU

        TYP  = NOID(N)(1:1)
        DIRN = NOID(N)(2:2)

        Z = NOZ1(N)
        Y = NOY1(N)
        X = NOX1(N)

        IF (TYP.EQ.'E') THEN
          IF (DIRN.EQ.'X') THEN
            OUT = V(1,X,Y,Z) + V(2,X,Y,Z) + V(12,X,Y,Z) + V(9 ,X,Y,Z)
          ELSE IF (DIRN.EQ.'Y') THEN
            OUT = V(4,X,Y,Z) + V(3,X,Y,Z) + V(8 ,X,Y,Z) + V(11,X,Y,Z)
          ELSE IF (DIRN.EQ.'Z') THEN
            OUT = V(6,X,Y,Z) + V(5,X,Y,Z) + V(10,X,Y,Z) + V(7 ,X,Y,Z)

          ELSE
            STOP 'E-field output error'
          END IF
          OUT = - OUT / DL / 2.0
 
        ELSE IF (TYP.EQ.'H') THEN
          IF (DIRN.EQ.'X') THEN
            OUT = V(4,X,Y,Z) - V(5,X,Y,Z) - V(8 ,X,Y,Z) + V(7 ,X,Y,Z)
          ELSE IF (DIRN.EQ.'Y') THEN
            OUT = V(6,X,Y,Z) - V(2,X,Y,Z) - V(10,X,Y,Z) + V(9 ,X,Y,Z)
          ELSE IF (DIRN.EQ.'Z') THEN
            OUT = V(1,X,Y,Z) - V(3,X,Y,Z) - V(12,X,Y,Z) + V(11,X,Y,Z)
          ELSE
            STOP 'H-field output error'
          END IF
          OUT = OUT / Z0 / DL / 2.0

        ELSE
          STOP 'Output field error'
        END IF

        WRITE(OUTUNT,99200) OUT
99200 FORMAT(1X,E13.6)

41010 CONTINUE

*---- scatter

      DO 42010, Z = 1, ZMAX
        DO 42020, Y = 1, YMAX
          DO 42030, X = 1, XMAX

            VYNX = V(1,X,Y,Z)
            VZNX = V(2,X,Y,Z)
            VXNY = V(3,X,Y,Z)
            VZNY = V(4,X,Y,Z)
            VYNZ = V(5,X,Y,Z)
            VXNZ = V(6,X,Y,Z)
            VYPZ = V(7,X,Y,Z)
            VZPY = V(8,X,Y,Z)
            VZPX = V(9,X,Y,Z)
            VXPZ = V(10,X,Y,Z)
            VXPY = V(11,X,Y,Z)
            VYPX = V(12,X,Y,Z)

            VDIFF = VXPY - VXNY
            VTEMP = 0.5 * ( VZNX + VZPX + VDIFF )
            V(12,X,Y,Z) = VTEMP
            V(1 ,X,Y,Z) = VTEMP - VDIFF

            VDIFF = VXPZ - VXNZ
            VTEMP = 0.5 * ( VYNX + VYPX + VDIFF )
            V(9 ,X,Y,Z) = VTEMP
            V(2 ,X,Y,Z) = VTEMP - VDIFF

            VDIFF = VYPZ - VYNZ
            VTEMP = 0.5 * ( VXNY + VXPY + VDIFF )
            V(8 ,X,Y,Z) = VTEMP
            V(4 ,X,Y,Z) = VTEMP - VDIFF

            VDIFF = VYPX - VYNX
            VTEMP = 0.5 * ( VZNY + VZPY + VDIFF )
            V(11,X,Y,Z) = VTEMP
            V(3 ,X,Y,Z) = VTEMP - VDIFF

            VDIFF = VZPX - VZNX
            VTEMP = 0.5 * ( VYNZ + VYPZ + VDIFF )
            V(10,X,Y,Z) = VTEMP
            V(6 ,X,Y,Z) = VTEMP - VDIFF

            VDIFF = VZPY - VZNY
            VTEMP = 0.5 * ( VXNZ + VXPZ + VDIFF )
            V(7 ,X,Y,Z) = VTEMP
            V(5 ,X,Y,Z) = VTEMP - VDIFF









	          LINI1     = W(1 ,X,Y,Z)*0.5
                LINI2     = W(2 ,X,Y,Z)*0.5 
                LINI3     = W(3 ,X,Y,Z)*0.5 
                LINI4     = W(4 ,X,Y,Z)*0.5 
                LINI5     = W(5 ,X,Y,Z)*0.5 
                LINI6     = W(6 ,X,Y,Z)*0.5 
                LINI7     = W(7 ,X,Y,Z)*0.5 
                LINI8     = W(8 ,X,Y,Z)*0.5
                LINI9     = W(9 ,X,Y,Z)*0.5 
                LINI10    = W(10,X,Y,Z)*0.5 
                LINI11    = W(11,X,Y,Z)*0.5 
                LINI12    = W(12,X,Y,Z)*0.5 

        W(1 ,X,Y,Z)=LINI2+LINI3+LINI9-LINI11
        W(2 ,X,Y,Z)=LINI1+LINI6-LINI10+LINI12
        W(3 ,X,Y,Z)=LINI1+LINI4+LINI8-LINI12
        W(4 ,X,Y,Z)=LINI3+LINI5-LINI7+LINI11
        W(5 ,X,Y,Z)=LINI4+LINI6-LINI8+LINI10
        W(6 ,X,Y,Z)=LINI2+LINI5+LINI7-LINI9
        W(7 ,X,Y,Z)=-LINI4+LINI6+LINI8+LINI10
        W(8 ,X,Y,Z)=LINI3-LINI5+LINI7+LINI11
        W(9 ,X,Y,Z)=LINI1-LINI6+LINI10+LINI12
        W(10,X,Y,Z)=-LINI2+LINI5+LINI7+LINI9
        W(11,X,Y,Z)=-LINI1+LINI4+LINI8+LINI12
        W(12,X,Y,Z)=LINI2-LINI3+LINI9+LINI11   

42030     CONTINUE


42020   CONTINUE
42010 CONTINUE

*---- boundaries

      DO 43010, N = 1, NUMBCS

        DIRN = BDIRN(N)
        LSIGN = BSIGN(N)
        RHO  = BRHO(N)

        IF (DIRN.EQ.'X') THEN
          X = BOX1(N)

          IF (LSIGN.EQ.'N'.AND.X.EQ.1) THEN
            DO 43110, Z = BOZ1(N), BOZ2(N)
              DO 43120, Y = BOY1(N), BOY2(N)
                V(3,X,Y,Z) = RHO * V(3,X,Y,Z)
                V(6,X,Y,Z) = RHO * V(6,X,Y,Z)
43120         CONTINUE
43110       CONTINUE

          ELSE IF (LSIGN.EQ.'P'.AND.X.EQ.XMAX) THEN
            DO 43130, Z = BOZ1(N), BOZ2(N)
              DO 43140, Y = BOY1(N), BOY2(N)
                V(10,X,Y,Z) = RHO * V(10,X,Y,Z)
                V(11,X,Y,Z) = RHO * V(11,X,Y,Z)
43140         CONTINUE
43130      CONTINUE

          ELSE
	    IF (LSIGN.EQ.'N') X = X - 1
	    X1 = X + 1
            DO 43150, Z = BOZ1(N), BOZ2(N)
              DO 43160, Y = BOY1(N), BOY2(N)
	        VTEMP = V(11,X,Y,Z)
		V(11,X,Y,Z) = RHO * V(3,X1,Y,Z)
		V(3,X1,Y,Z) = RHO * VTEMP
	        VTEMP = V(10,X,Y,Z)
		V(10,X,Y,Z) = RHO * V(6,X1,Y,Z)
		V(6,X1,Y,Z) = RHO * VTEMP
43160         CONTINUE
43150       CONTINUE

          END IF

        ELSE IF (DIRN.EQ.'Y') THEN
          Y = BOY1(N)

          IF (LSIGN.EQ.'N'.AND.Y.EQ.1) THEN
            DO 43210, Z = BOZ1(N), BOZ2(N)
              DO 43220, X = BOX1(N), BOX2(N)
                V(5,X,Y,Z) = RHO * V(5,X,Y,Z)
                V(1,X,Y,Z) = RHO * V(1,X,Y,Z)
43220         CONTINUE
43210       CONTINUE

          ELSE IF (LSIGN.EQ.'P'.AND.Y.EQ.YMAX) THEN
            DO 43230, Z = BOZ1(N), BOZ2(N)
              DO 43240, X = BOX1(N), BOX2(N)
                V(7,X,Y,Z) = RHO * V(7,X,Y,Z)
                V(12,X,Y,Z) = RHO * V(12,X,Y,Z)
43240         CONTINUE
43230       CONTINUE

          ELSE
	    IF (LSIGN.EQ.'N') Y = Y - 1
	    Y1 = Y + 1
            DO 43250, Z = BOZ1(N), BOZ2(N)
              DO 43260, X = BOX1(N), BOX2(N)
                VTEMP = V(7,X,Y,Z)
		V(7,X,Y,Z) = RHO * V(5,X,Y1,Z)
		V(5,X,Y1,Z) = RHO * VTEMP
		VTEMP = V(12,X,Y,Z)
		V(12,X,Y,Z) = RHO * V(1,X,Y1,Z)
		V(1,X,Y1,Z) = RHO * VTEMP
43260         CONTINUE
43250       CONTINUE

          END IF

        ELSE IF (DIRN.EQ.'Z') THEN
          Z = BOZ1(N)

          IF (LSIGN.EQ.'N'.AND.Z.EQ.1) THEN
            DO 43310, Y = BOY1(N), BOY2(N)
              DO 43320, X = BOX1(N), BOX2(N)
                V(2,X,Y,Z) = RHO * V(2,X,Y,Z)
                V(4,X,Y,Z) = RHO * V(4,X,Y,Z)
43320         CONTINUE
43310       CONTINUE

          ELSE IF (LSIGN.EQ.'P'.AND.Z.EQ.ZMAX) THEN
            DO 43330, Y = BOY1(N), BOY2(N)
              DO 43340, X = BOX1(N), BOX2(N)
                V(9,X,Y,Z) = RHO * V(9,X,Y,Z)
                V(8,X,Y,Z) = RHO * V(8,X,Y,Z)
43340         CONTINUE
43330       CONTINUE

          ELSE
	    IF (LSIGN.EQ.'N') Z = Z - 1
	    Z1 = Z + 1
            DO 43350, Y = BOY1(N), BOY2(N)
              DO 43360, X = BOX1(N), BOX2(N)
                VTEMP = V(9,X,Y,Z)
                V(9,X,Y,Z) = RHO * V(2,X,Y,Z1)
		V(2,X,Y,Z1) = RHO * VTEMP
		VTEMP = V(8,X,Y,Z)
		V(8,X,Y,Z) = RHO * V(4,X,Y,Z1)
		V(4,X,Y,Z1) = RHO * VTEMP
43360         CONTINUE
43350       CONTINUE

          END IF

        ELSE
          STOP 'Boundary direction error'
        END IF

43010 CONTINUE

*---- connect

      DO 44010, Z = 1, ZMAX
        Z1 = Z + 1
        DO 44020, Y = 1, YMAX
          Y1 = Y + 1
          DO 44030, X = 1, XMAX
            X1 = X + 1
 
            IF (X.LT.XMAX) THEN
              V0 = V(10,X,Y,Z)
              V(10,X,Y,Z) = V(6,X1,Y,Z)
              V(6,X1,Y,Z) = V0
              V0 = V(11,X,Y,Z)
              V(11,X,Y,Z) = V(3,X1,Y,Z)
              V(3,X1,Y,Z) = V0

	                  TEMPI1      = W(11,X,Y,Z)
                        TEMPI2      = W(10,X,Y,Z)
                        W(11,X,Y,Z) = W(3,X1,Y,Z)
                        W(10,X,Y,Z) = W(6,X1,Y,Z)
                        W(3,X1,Y,Z) = TEMPI1
                        W(6,X1,Y,Z) = TEMPI2
            END IF
 
            IF (Y.LT.YMAX) THEN
              V0 = V(7,X,Y,Z)
              V(7,X,Y,Z) = V(5,X,Y1,Z)
              V(5,X,Y1,Z) = V0
              V0 = V(12,X,Y,Z)
              V(12,X,Y,Z) = V(1,X,Y1,Z)
              V(1,X,Y1,Z) = V0

	                  TEMPI1      = W(12,X,Y,Z)
                        TEMPI2      = W(7 ,X,Y,Z)
                        W(12,X,Y,Z) = W(1,X,Y,Z)
                        W(7 ,X,Y,Z) = W(5,X,Y,Z)
                        W(1,X,Y,Z) = TEMPI1
                        W(5,X,Y,Z) = TEMPI2
            END IF
 
            IF (Z.LT.ZMAX) THEN
              V0 = V(8,X,Y,Z)
              V(8,X,Y,Z) = V(4,X,Y,Z1)
              V(4,X,Y,Z1) = V0
              V0 = V(9,X,Y,Z)
              V(9,X,Y,Z) = V(2,X,Y,Z1)
              V(2,X,Y,Z1) = V0

                        TEMPI1      = W(8,X,Y,Z)
                        TEMPI2      = W(9,X,Y,Z)
                        W(8,X,Y,Z)  = W(4,X,Y,Z)
                        W(9,X,Y,Z)  = W(2,X,Y,Z)
                        W(4,X,Y,Z) = TEMPI1
                        W(2,X,Y,Z) = TEMPI2
            END IF
 
44030     CONTINUE
44020   CONTINUE
44010 CONTINUE

*---- end calculation

30010 CONTINUE

*---- close output file

      CLOSE(OUTUNT)
      WRITE(*,*) 'Output file closed'
      WRITE(*,*)
      END