*-----------------------------------------------------------------------
*     TMSTLN : THE ROUTINE FOR DRAWING STREAM LINES
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2022 GFD DENNOU CLUB. ALL RIGHTS RESERVED.
*-----------------------------------------------------------------------
      SUBROUTINE TMSTLN( X, Y, U, V, NX, NY, TRX, TRY, TOTNM )

      IMPLICIT NONE

*-- AGREENMENT
      INTEGER NX                  !* GRID NUMBERS FOR X-DIRECTION
      INTEGER NY                  !* GRID NUMBERS FOR Y-DIRECTION
      REAL    U(NX,NY)            !* VECTOR COMPONENT OF X-DIRECTION
      REAL    V(NX,NY)            !* VECTOR COMPONENT OF Y-DIRECTION
      REAL    X(NX)               !* GRID POINTS OF X-DIRECTION [UNIT:LENGTH]
      REAL    Y(NY)               !* GRID POINTS OF Y-DIRECTION [UNIT:LENGTH]
      REAL    TRX(NX*NY,2)        !* STREAM LINE OF X-COORDINATE
      REAL    TRY(NX*NY,2)        !* STREAM LINE OF Y-COORDINATE
      INTEGER TOTNM(NX*NY,2)      !* ARRAY NUMBER IN TOTAL STREAM LINES AT END

*-- INTERNAL VARIABLES
      INTEGER I, J, JD, K, L    !* TMP VALUES
      INTEGER COUNTER           !* NUMBER OF STREAM LINE
      INTEGER DCOUNTER          !* NUMBER OF STREAM LINE FOR DRAWING
      INTEGER NM(NX*NY,2)         !* ARRAY NUMBER IN EACH STREAM LINE
      INTEGER GRID_NUM(NX,NY)
      INTEGER ICOUNTER
      INTEGER DICOUNTER
      INTEGER TMPCOUNTER
      INTEGER THRES
      INTEGER SKIP
      INTEGER GLIM
      INTEGER NXBOUND(2)          !* BOUNDARY GRID NUMBER OF X COMPONENT
      INTEGER NYBOUND(2)          !* BOUNDARY GRID NUMBER OF Y COMPONENT
      REAL    ARROW_THRES
      REAL    DT
      REAL    DTX, DTY
      REAL    UNDEF             !* UNDEFINED VALUE FOR DCL
      REAL    TMPRX(NX*NY,2)
      REAL    TMPRY(NX*NY,2)
      REAL    TMPX(NX*NY)
      REAL    TMPY(NX*NY)
      CHARACTER(10) RCHARC
      CHARACTER(3)  ICHARC
      LOGICAL PFX, PFY, DTV

      CALL GLRGET( 'RMISS', UNDEF )

*-- GETTING TMPACK'S PARAMETERS
      CALL TMIGET( 'STLNGLIM', GLIM )
      CALL TMIGET( 'GRDTHRES', THRES )
      CALL TMIGET( 'SKIPINTV', SKIP )
      CALL TMRGET( 'ARRWINTV', ARROW_THRES )
      CALL TMLGET( 'PERIODX', PFX )
      CALL TMLGET( 'PERIODY', PFY )
      CALL TMLGET( 'FIXEDDT', DTV )
      CALL TMRGET( 'STLNDT', DT )
*      CALL TMVGET( THRES, SKIP, ARROW_THRES )  !* STPACK GETTING ROUTINE

      NXBOUND(1)=1
      NXBOUND(2)=NX
      NYBOUND(1)=1
      NYBOUND(2)=NY

      TRX(1:NX*NY,1:2)=0.0
      TRY(1:NX*NY,1:2)=0.0
      TOTNM(1:NX*NY,1:2)=0
      NM(1:NX*NY,1:2)=0

*-- CHECKING VARYING DT
      IF(DT.EQ.0.0)THEN
         DTV=.FALSE.
         CALL TMLSET( 'FIXEDDT', .FALSE. )
      END IF

*-- INITIALIZING GRID_NUM (COUNTER FOR PASSING STREAM LINES IN EACH GRID)
*-- 格子点通過流線のカウンタ変数を初期化
      DO 11 J=1,NY
         DO 10 I=1,NX
            GRID_NUM(I,J)=0
 10      CONTINUE
 11   CONTINUE

      COUNTER=1
      ICOUNTER=1
      DCOUNTER=1
      DICOUNTER=1

*-- DETERMINING "DT" (DT=DXMIN/VMAX)
*-- 以下で流線計算の際の積分における
*-- CFL 条件のようなものを決定する.
**      MAXV=0.0

**      DO 31 J=1,NY
**         DO 30 I=1,NX
**            IF(U(I,J).NE.UNDEF.AND.V(I,J).NE.UNDEF)THEN
**               VTMP=SQRT(U(I,J)**2+V(I,J)**2)
**               IF(VTMP.GT.MAXV)THEN
**                  MAXV=VTMP
**               END IF
**            END IF
** 30      CONTINUE
** 31   CONTINUE

**      MINR=X(2)-X(1)
**
**      DO 40 I=2,NX
**         XTMP=X(I)-X(I-1)
**         IF(XTMP.LT.MINR)THEN
**            MINR=XTMP
**         END IF
** 40   CONTINUE
**
**      DO 41 J=2,NY
**         XTMP=Y(J)-Y(J-1)
**         IF(XTMP.LT.MINR)THEN
**            MINR=XTMP
**         END IF
** 41   CONTINUE

*-- INITIALIZING

      IF(DT.LE.0.0)THEN
         DO 31 J=2,NY
            DO 30 I=2,NX
               IF(U(I,J).NE.UNDEF.AND.U(I-1,J).NE.UNDEF.AND.
     &            V(I,J).NE.UNDEF.AND.V(I,J-1).NE.UNDEF)THEN
                  IF(U(I,J)+U(I-1,J).NE.0.0)THEN
                     DTX=2.0*abs((X(I)-X(I-1))/(U(I,J)+U(I-1,J)))
                  ELSE
                     DTX=abs(DT)
                  END IF
                  IF(V(I,J)+V(I,J-1).NE.0.0)THEN
                     DTY=2.0*abs((Y(J)-Y(J-1))/(V(I,J)+V(I,J-1)))
                  ELSE
                     DTY=abs(DT)
                  END IF
        
                  IF(DT.NE.0.0)THEN
                     IF(ABS(DT).GT.DTX)THEN
                        DT=DTX
                     END IF
                     IF(ABS(DT).GT.DTY)THEN
                        DT=DTY
                     END IF
                  ELSE
                     IF(DTX.NE.0.0)THEN
                        DT=DTX
                     END IF
                     IF(DTY.NE.0.0)THEN
                        DT=DTY
                     END IF
                  END IF
               END IF
 30         CONTINUE
 31      CONTINUE
      END IF

      CALL MSGDMP( 'M', 'TMSTLN', 'CHECKING EACH PARAMETER' )
      WRITE(ICHARC(1:3),'(I3)') THRES
      CALL MSGDMP( 'M', 'TMSTLN', 'THRES       ===> '//ICHARC(1:3) )
      WRITE(ICHARC(1:3),'(I3)') SKIP
      CALL MSGDMP( 'M', 'TMSTLN', 'SKIP        ===> '//ICHARC(1:3) )
      WRITE(ICHARC(1:3),'(I3)') GLIM
      CALL MSGDMP( 'M', 'TMSTLN', 'GLIM        ===> '//ICHARC(1:3) ) 
      WRITE(RCHARC(1:10),'(1P,E10.3)') ARROW_THRES
      CALL MSGDMP( 'M', 'TMSTLN', 'ARROW THRES ===> '//RCHARC(1:10) )
      WRITE(RCHARC(1:10),'(1P,E10.3)') DT
      CALL MSGDMP( 'M', 'TMSTLN', 'DT          ===> '//RCHARC(1:10) ) 

*-- 以上の擬似 CFL 条件から時間積分ステップを決める.
**      DT=MINR/MAXV

*-- STARTING TO CALCULATE STREAM LINES.
*-- DETERMINING STARTING POINT BY USING A PARAMETER VALUE "SKIP".

*-- [1] 領域境界上を始点として流線計算を行う
*-- 以下の計算は [2] の 61, 60 ループの引数が異なるのみで,
*-- ループ内部の処理は全く同じである.
*-- 161 は X 軸境界, 261 は Y 軸境界を始点とする.
      DO 161 JD=1,2
         J=NYBOUND(JD)
         IF(MOD((J-1),SKIP).EQ.0.OR.J.EQ.NY)THEN
            DO 160 I=1,NX
               IF(MOD((I-1),SKIP).EQ.0.OR.I.EQ.NX)THEN
                  IF(U(I,J).NE.0.0.OR.V(I,J).NE.0.0)THEN  !* 格子点上で風速がゼロなら TMSLCL に入らない.
                     TMPRX(1,1)=X(I)
                     TMPRY(1,1)=Y(J)
                     TMPCOUNTER=COUNTER
                     IF(GRID_NUM(I,J).LE.THRES)THEN
                        DO 162 K=TMPCOUNTER,NX*NY  !* この閾値は取りすぎな気もする.
                           CALL TMRSET( 'STLNDT', ABS(DT) )
                           CALL TMSLCL( TMPRX(1,1), TMPRY(1,1),
     &                                  NX, NY, X, Y, U, V,
     &                                  TMPRX(:,1), TMPRY(:,1),
     &                                  NM(COUNTER,1), GRID_NUM )
                           COUNTER=COUNTER+1

                           !* サブルーチンの返り値へ値を入れる.
                           IF(NM(COUNTER-1,1).GT.1)THEN
                             IF(DCOUNTER.GT.1)THEN
                               TOTNM(DCOUNTER,1)=TOTNM(DCOUNTER-1,1)
     &                                           +NM(COUNTER-1,1)
                               DO 165 L=1,NM(COUNTER-1,1)
                                 TRX(TOTNM(DCOUNTER-1,1)+L,1)=TMPRX(L,1)
                                 TRY(TOTNM(DCOUNTER-1,1)+L,1)=TMPRY(L,1)
 165                           CONTINUE
                             ELSE
                               TOTNM(1,1)=NM(1,1)
                               TRX(1:NM(1,1),1)=TMPRX(1:NM(1,1),1)
                               TRY(1:NM(1,1),1)=TMPRY(1:NM(1,1),1)
                             END IF
                             DCOUNTER=DCOUNTER+1
                           END IF

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(COUNTER-1,1).LT.NX*NY)THEN
                              IF(PFX.EQV..TRUE.)THEN  !* 終わっていれば, 周期境界なら境界上が開始点での逆計算へ, 周期境界でなければ次の流線計算へ.
                                IF(I.EQ.1.OR.I.EQ.NX)THEN
                                   GO TO 163
                                ELSE
                                   GO TO 990
                                END IF
                              ELSE
                                 GO TO 990
                              END IF

 990                          IF(PFY.EQV..TRUE.)THEN
                                IF(J.EQ.1.OR.J.EQ.NY)THEN
                                   GO TO 163
                                ELSE
                                   GO TO 160
                                END IF
                              ELSE
                                 GO TO 163  !* 終わっていれば逆向き計算へ.
                              END IF
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,1)=TMPRX(NX*NY,1)
                              TMPRY(1,1)=TMPRY(NX*NY,1)
                           END IF

                           IF(K.EQ.NX*NY)THEN
                              GO TO 163
                           END IF

                           IF(COUNTER.EQ.NX*NY)THEN
                              GO TO 163
                           ELSE IF(COUNTER.GT.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF

 162                    CONTINUE

*-- CALCULATING INVERSE DIRECTION.
*-- "TMSLCL" CAN CALCULATE FORWARD DIRECTION FOR "DT>0",
*-- BACKWARD DIRECTION FOR "DT<0".
*-- 流線が領域内で開始されているので, 逆向きにも計算する.
*-- DT に負号をつけて逆探査開始.
*-- GRID_NUM(I,J) は前方計算でカウントされているので, このカウントを外しておく
 163                    IF(ICOUNTER.LT.NX*NY)THEN
                           GRID_NUM(I,J)=GRID_NUM(I,J)-1
                           TMPRX(1,2)=X(I)
                           TMPRY(1,2)=Y(J)
                        ELSE
                           CALL MSGDMP( 'W', 'TMSTLN', 
     &                                  'ARRAY WILL BE OVER FLOW.')
                           GO TO 100
                        END IF

                        TMPCOUNTER=ICOUNTER

                        DO 164 K=TMPCOUNTER,NX*NY
                           CALL TMRSET( 'STLNDT', -ABS(DT) )  ! 逆向き計算のため, 負号つき
                           CALL TMSLCL( TMPRX(1,2), 
     &                                  TMPRY(1,2),
     &                                  NX, NY, X, Y, U, V, 
     &                                  TMPX(:), TMPY(:),
     &                                  NM(ICOUNTER,2), GRID_NUM )

*-- 逆計算分を TRX, TRY に格納.
*-- ただし, 逆向き計算なので, 配列は逆から入れる.
                           IF(NM(ICOUNTER,2).GT.1)THEN
                              DO 170 L=1,NM(ICOUNTER,2)
                                 TMPRX(L,2)=TMPX(NM(ICOUNTER,2)-L+1)
                                 TMPRY(L,2)=TMPY(NM(ICOUNTER,2)-L+1)
 170                          CONTINUE

                              !* サブルーチンの返り値へ値を入れる.
                              IF(DICOUNTER.GT.1)THEN
                                TOTNM(DICOUNTER,2)=TOTNM(DICOUNTER-1,2)
     &                                             +NM(ICOUNTER,2)
                                DO 171 L=1,NM(ICOUNTER,2)
                                TRX(TOTNM(DICOUNTER-1,2)+L,2)=TMPRX(L,2)
                                TRY(TOTNM(DICOUNTER-1,2)+L,2)=TMPRY(L,2)
 171                            CONTINUE
                              ELSE
                                TOTNM(1,2)=NM(1,2)
                                TRX(1:NM(1,2),2)=TMPRX(1:NM(1,2),2)
                                TRY(1:NM(1,2),2)=TMPRY(1:NM(1,2),2)
                              END IF
                              DICOUNTER=DICOUNTER+1
                           END IF
                           ICOUNTER=ICOUNTER+1

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(ICOUNTER-1,2).LT.NX*NY)THEN
                              GO TO 160  !* 終わっていれば次の流線計算へ
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,2)=TMPX(NX*NY)
                              TMPRY(1,2)=TMPY(NX*NY)
                           END IF
                           IF(K.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
                           IF(ICOUNTER.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
 164                    CONTINUE
                     END IF
                  END IF
               END IF
 160        CONTINUE
         END IF
 161  CONTINUE

      COUNTER=COUNTER+1
      ICOUNTER=ICOUNTER+1

      DO 261 JD=1,2
         I=NXBOUND(JD)
         DO 260 J=2,NY-1  ! 端領域は上のループで行っている.
            IF(MOD((J-1),SKIP).EQ.0.OR.J.EQ.NY)THEN
               IF(MOD((I-1),SKIP).EQ.0.OR.I.EQ.NX)THEN
                  IF(U(I,J).NE.0.0.OR.V(I,J).NE.0.0)THEN  !* 格子点上で風速がゼロなら TMSLCL に入らない.
                     TMPRX(1,1)=X(I)
                     TMPRY(1,1)=Y(J)
                     TMPCOUNTER=COUNTER
                     IF(GRID_NUM(I,J).LE.THRES)THEN
                        DO 262 K=TMPCOUNTER,NX*NY  !* この閾値は取りすぎな気もする.
                           CALL TMRSET( 'STLNDT', ABS(DT) )
                           CALL TMSLCL( TMPRX(1,1), TMPRY(1,1),
     &                                  NX, NY, X, Y, U, V,
     &                                  TMPRX(:,1), TMPRY(:,1),
     &                                  NM(COUNTER,1), GRID_NUM )
                           COUNTER=COUNTER+1

                           !* サブルーチンの返り値へ値を入れる.
                           IF(NM(COUNTER-1,1).GT.1)THEN
                             IF(DCOUNTER.GT.1)THEN
                               TOTNM(DCOUNTER,1)=TOTNM(DCOUNTER-1,1)
     &                                           +NM(COUNTER-1,1)
                               DO 265 L=1,NM(COUNTER-1,1)
                                 TRX(TOTNM(DCOUNTER-1,1)+L,1)=TMPRX(L,1)
                                 TRY(TOTNM(DCOUNTER-1,1)+L,1)=TMPRY(L,1)
 265                           CONTINUE
                             ELSE
                               TOTNM(1,1)=NM(1,1)
                               TRX(1:NM(1,1),1)=TMPRX(1:NM(1,1),1)
                               TRY(1:NM(1,1),1)=TMPRY(1:NM(1,1),1)
                             END IF
                             DCOUNTER=DCOUNTER+1
                           END IF

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(COUNTER-1,1).LT.NX*NY)THEN
                              IF(PFX.EQV..TRUE.)THEN !* 終わっていれば, 周期境界なら境界上が開始点での逆計算へ, 周期境界でなければ次の流線計算へ.
                                IF(I.EQ.1.OR.I.EQ.NX)THEN
                                   GO TO 263
                                ELSE
                                   GO TO 991
                                END IF
                              ELSE
                                 GO TO 991
                              END IF

 991                          IF(PFY.EQV..TRUE.)THEN
                                IF(J.EQ.1.OR.J.EQ.NY)THEN
                                   GO TO 263
                                ELSE
                                   GO TO 260
                                END IF
                              ELSE
                                 GO TO 263  !* 終わっていれば逆向き計算へ.
                              END IF
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,1)=TMPRX(NX*NY,1)
                              TMPRY(1,1)=TMPRY(NX*NY,1)
                           END IF

                           IF(K.EQ.NX*NY)THEN
                              GO TO 263
                           END IF

                           IF(COUNTER.EQ.NX*NY)THEN
                              GO TO 263
                           ELSE IF(COUNTER.GT.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF

 262                    CONTINUE

*-- CALCULATING INVERSE DIRECTION.
*-- "TMSLCL" CAN CALCULATE FORWARD DIRECTION FOR "DT>0",
*-- BACKWARD DIRECTION FOR "DT<0".
*-- 流線が領域内で開始されているので, 逆向きにも計算する.
*-- DT に負号をつけて逆探査開始.
*-- GRID_NUM(I,J) は前方計算でカウントされているので, このカウントを外しておく
 263                    IF(ICOUNTER.LT.NX*NY)THEN
                           GRID_NUM(I,J)=GRID_NUM(I,J)-1
                           TMPRX(1,2)=X(I)
                           TMPRY(1,2)=Y(J)
                        ELSE
                           CALL MSGDMP( 'W', 'TMSTLN', 
     &                                  'ARRAY WILL BE OVER FLOW.')
                           GO TO 100
                        END IF

                        TMPCOUNTER=ICOUNTER

                        DO 264 K=TMPCOUNTER,NX*NY
                           CALL TMRSET( 'STLNDT', -ABS(DT) )  ! 逆向き計算のため, 負号つき
                           CALL TMSLCL( TMPRX(1,2), 
     &                                  TMPRY(1,2),
     &                                  NX, NY, X, Y, U, V, 
     &                                  TMPX(:), TMPY(:),
     &                                  NM(ICOUNTER,2), GRID_NUM )

*-- 逆計算分を TRX, TRY に格納.
*-- ただし, 逆向き計算なので, 配列は逆から入れる.
                           IF(NM(ICOUNTER,2).GT.1)THEN
                              DO 270 L=1,NM(ICOUNTER,2)
                                 TMPRX(L,2)=TMPX(NM(ICOUNTER,2)-L+1)
                                 TMPRY(L,2)=TMPY(NM(ICOUNTER,2)-L+1)
 270                          CONTINUE

                              !* サブルーチンの返り値へ値を入れる.
                              IF(DICOUNTER.GT.1)THEN
                                TOTNM(DICOUNTER,2)=TOTNM(DICOUNTER-1,2)
     &                                             +NM(ICOUNTER,2)
                                DO 271 L=1,NM(ICOUNTER,2)
                                TRX(TOTNM(DICOUNTER-1,2)+L,2)=TMPRX(L,2)
                                TRY(TOTNM(DICOUNTER-1,2)+L,2)=TMPRY(L,2)
 271                            CONTINUE
                              ELSE
                                TOTNM(1,2)=NM(1,2)
                                TRX(1:NM(1,2),2)=TMPRX(1:NM(1,2),2)
                                TRY(1:NM(1,2),2)=TMPRY(1:NM(1,2),2)
                              END IF
                              DICOUNTER=DICOUNTER+1
                           END IF
                           ICOUNTER=ICOUNTER+1

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(ICOUNTER-1,2).LT.NX*NY)THEN
                              GO TO 260  !* 終わっていれば次の流線計算へ
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,2)=TMPX(NX*NY)
                              TMPRY(1,2)=TMPY(NX*NY)
                           END IF
                           IF(K.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
                           IF(ICOUNTER.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
 264                    CONTINUE
                     END IF
                  END IF
               END IF
            END IF
 260     CONTINUE
 261  CONTINUE

      COUNTER=COUNTER+1
      ICOUNTER=ICOUNTER+1

*-- [2] 領域内部の空白域で流線計算を行う.
*-- GRID_NUM によって, 通過していない格子は検索されているので, 
*-- その格子について SKIP 間隔で流線計算を行う.
      DO 61 J=2,NY-1
         IF(MOD((J-1),SKIP).EQ.0.OR.J.EQ.NY)THEN
            DO 60 I=2,NX-1
               IF(MOD((I-1),SKIP).EQ.0.OR.I.EQ.NX)THEN
                  IF(U(I,J).NE.0.0.OR.V(I,J).NE.0.0)THEN  !* 格子点上で風速がゼロなら TMSLCL に入らない.
                     TMPRX(1,1)=X(I)
                     TMPRY(1,1)=Y(J)
                     TMPCOUNTER=COUNTER
                     IF(GRID_NUM(I,J).LE.THRES)THEN
                        DO 62 K=TMPCOUNTER,NX*NY  !* この閾値は取りすぎな気もする.
                           CALL TMRSET( 'STLNDT', ABS(DT) )
                           CALL TMSLCL( TMPRX(1,1), TMPRY(1,1),
     &                                  NX, NY, X, Y, U, V,
     &                                  TMPRX(:,1), TMPRY(:,1),
     &                                  NM(COUNTER,1), GRID_NUM )
                           COUNTER=COUNTER+1

                           !* サブルーチンの返り値へ値を入れる.
                           IF(NM(COUNTER-1,1).GT.1)THEN
                             IF(DCOUNTER.GT.1)THEN
                               TOTNM(DCOUNTER,1)=TOTNM(DCOUNTER-1,1)
     &                                           +NM(COUNTER-1,1)
                               DO 65 L=1,NM(COUNTER-1,1)
                                 TRX(TOTNM(DCOUNTER-1,1)+L,1)=TMPRX(L,1)
                                 TRY(TOTNM(DCOUNTER-1,1)+L,1)=TMPRY(L,1)
 65                            CONTINUE
                             ELSE
                               TOTNM(1,1)=NM(1,1)
                               TRX(1:NM(1,1),1)=TMPRX(1:NM(1,1),1)
                               TRY(1:NM(1,1),1)=TMPRY(1:NM(1,1),1)
                             END IF
                             DCOUNTER=DCOUNTER+1
                           END IF

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(COUNTER-1,1).LT.NX*NY)THEN
                              IF(I.NE.1.AND.I.NE.NX.AND.
     &                           J.NE.1.AND.J.NE.NY)THEN
                                 GO TO 63  !* 終わっていれば逆計算へ (開始点が領域内側) .
                              ELSE  !* 終わっていれば, 周期境界なら境界上が開始点での逆計算へ, 周期境界でなければ次の流線計算へ.
                                 IF(PFX.EQV..TRUE.)THEN
                                   IF(I.EQ.1.OR.I.EQ.NX)THEN
                                      GO TO 63
                                   ELSE
                                      GO TO 999
                                   END IF
                                 ELSE
                                    GO TO 999
                                 END IF

 999                             IF(PFY.EQV..TRUE.)THEN
                                   IF(J.EQ.1.OR.J.EQ.NY)THEN
                                      GO TO 63
                                   ELSE
                                      GO TO 60
                                   END IF
                                 ELSE
                                    GO TO 60  !* 終わっていれば次の流線計算へ (開始点が領域境界)
                                 END IF
                              END IF
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,2)=TMPX(NX*NY)
                              TMPRY(1,2)=TMPY(NX*NY)
                           END IF

                           IF(K.EQ.NX*NY)THEN
                              GO TO 63
                           END IF

                           IF(COUNTER.EQ.NX*NY)THEN
                              GO TO 63
                           ELSE IF(COUNTER.GT.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF

 62                     CONTINUE

*-- CALCULATING INVERSE DIRECTION.
*-- "TMSLCL" CAN CALCULATE FORWARD DIRECTION FOR "DT>0",
*-- BACKWARD DIRECTION FOR "DT<0".
*-- 流線が領域内で開始されているので, 逆向きにも計算する.
*-- DT に負号をつけて逆探査開始.
*-- GRID_NUM(I,J) は前方計算でカウントされているので, このカウントを外しておく
 63                     IF(ICOUNTER.LT.NX*NY)THEN
                           GRID_NUM(I,J)=GRID_NUM(I,J)-1
                           TMPRX(1,2)=X(I)
                           TMPRY(1,2)=Y(J)
                        ELSE
                           CALL MSGDMP( 'W', 'TMSTLN', 
     &                                  'ARRAY WILL BE OVER FLOW.')
                           GO TO 100
                        END IF

                        TMPCOUNTER=ICOUNTER

                        DO 64 K=TMPCOUNTER,NX*NY
                           CALL TMRSET( 'STLNDT', -ABS(DT) )  ! 逆向き計算のため, 負号つき
                           CALL TMSLCL( TMPRX(1,2), TMPRY(1,2),
     &                                  NX, NY, X, Y, U, V, 
     &                                  TMPX(:), TMPY(:),
     &                                  NM(ICOUNTER,2), GRID_NUM )

*-- 逆計算分を TRX, TRY に格納.
*-- ただし, 逆向き計算なので, 配列は逆から入れる.
                           IF(NM(ICOUNTER,2).GT.1)THEN
                              DO 70 L=1,NM(ICOUNTER,2)
                                 TMPRX(L,2)=TMPX(NM(ICOUNTER,2)-L+1)
                                 TMPRY(L,2)=TMPY(NM(ICOUNTER,2)-L+1)
 70                           CONTINUE

                              !* サブルーチンの返り値へ値を入れる.
                              IF(DICOUNTER.GT.1)THEN
                                TOTNM(DICOUNTER,2)=TOTNM(DICOUNTER-1,2)
     &                                             +NM(ICOUNTER,2)
                                DO 71 L=1,NM(ICOUNTER,2)
                                TRX(TOTNM(DICOUNTER-1,2)+L,2)=TMPRX(L,2)
                                TRY(TOTNM(DICOUNTER-1,2)+L,2)=TMPRY(L,2)
 71                             CONTINUE
                              ELSE
                                TOTNM(1,2)=NM(1,2)
                                TRX(1:NM(1,2),2)=TMPRX(1:NM(1,2),2)
                                TRY(1:NM(1,2),2)=TMPRY(1:NM(1,2),2)
                              END IF
                              DICOUNTER=DICOUNTER+1
                           END IF
                           ICOUNTER=ICOUNTER+1

*-- 流線が用意した配列以上に連続する場合, ここで終わるまで計算させる.
                           IF(NM(ICOUNTER-1,2).LT.NX*NY)THEN
                              GO TO 60  !* 終わっていれば次の流線計算へ
                           ELSE  !* 用意した配列を越えて流線が続く場合の引き継ぎ
                              TMPRX(1,2)=TMPX(NX*NY)
                              TMPRY(1,2)=TMPY(NX*NY)
                           END IF
                           IF(K.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
                           IF(ICOUNTER.EQ.NX*NY)THEN
                              CALL MSGDMP( 'W', 'TMSTLN', 
     &                                     'ARRAY WILL BE OVER FLOW.')
                              GO TO 100
                           END IF
 64                     CONTINUE
                     END IF
                  END IF
               END IF
 60         CONTINUE
         END IF
 61   CONTINUE

 100  CALL MSGDMP( 'M', 'TMSTLN', 'CALCULATING FINISHED.' )

      COUNTER=COUNTER-1

      CALL TMISET( 'STLNNUM', DCOUNTER )

      END SUBROUTINE TMSTLN

