*-----------------------------------------------------------------------
*     TMSLCL : MAIN SOLVER ROUTINE FOR CALCULATING STREAM LINES (RK4)
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2022 GFD DENNOU CLUB. ALL RIGHTS RESERVED.
*-----------------------------------------------------------------------
      SUBROUTINE TMSLCL( INI_X, INI_Y, NX, NY, X, Y, U, V, 
     &                   TRAJX, TRAJY, TRAJN, GRID_FLAG )
      IMPLICIT NONE
      REAL    INI_X                ! STARTING POINT (X) [UNIT:LENGTH]
      REAL    INI_Y                ! STARTING POINT (Y) [UNIT:LENGTH]
      INTEGER NX                   ! GRID NUMBER OF X-DIRECTION
      INTEGER NY                   ! GRID NUMBER OF Y-DIRECTION
      REAL    X(NX)                ! X-COORDINATE [UNIT:LENGTH]
      REAL    Y(NY)                ! Y-COORDINATE [UNIT:LENGTH]
      REAL    U(NX,NY)             ! VECTOR COMPONENT OF X-DIRECTION
      REAL    V(NX,NY)             ! VECTOR COMPONENT OF Y-DIRECTION
      REAL    TRAJX(NX*NY)         ! POSITIONS OF STREAM LINE (X) [UNIT:LENGTH]
      REAL    TRAJY(NX*NY)         ! POSITIONS OF STREAM LINE (Y) [UNIT:LENGTH]
      INTEGER TRAJN                ! RANGE OF DRAWING STREAM LINE (STEP>=TRAJN)
      INTEGER GRID_FLAG(NX,NY)     ! NUMBER OF PASSING STREAM LINES IN EACH GRID.

*-- INTERNAL VALUES
      INTEGER STEP                 ! CALCULATING STEP NUMBER.
      INTEGER THRES                ! NUMBER OF PERMITTING TO PASS STREAM LINE
                                   ! IN EACH GRID (RECOMMENDED VALUE = 1).
      INTEGER I, J, M, N, ID
      INTEGER STNX, STNY           ! STARTING GRID NUMBER OF X, Y COMPONENTS.
      INTEGER INI_M, INI_N         ! INITIAL GRID NUMBER OF X, Y COMPONENTS.
      INTEGER GLIM                 ! MINIMUM LIMIT OF STREAM LINE INTERVAL
      INTEGER SGRID                ! COUNTING THE EXISTING OF THE SAME GRID.
      REAL    K1, K2, K3, K4
      REAL    L1, L2, L3, L4
      REAL    X1, Y1
      REAL    DTTX, DTTY
      REAL    DT                   ! TIME INTERVAL TO CALCULATE STREAM LINE.
      REAL    UNDEF                ! UNDEFINED VALUE FOR (U,V)
      REAL    INTER_P(2)           ! POSITION IN INTERPOLATING POINT
      REAL    INTER_V(2)           ! VECTOR VALUE IN INTERPOLATING POINT
      CHARACTER(1) DTPM            ! FLAG OF FORWARD OR BACKWARD CALCULATING.
      LOGICAL CONTINU_FLAG(NX,NY)  ! FLAG VALUE TO PREVENT THE IDENTICAL 
                                   ! STREAM LINE FROM CONTINUOUS COUNTING
                                   ! IN EACH GRID.
                                   ! 次ステップで同じ格子点にいるのを
                                   ! カウントすることを防ぐためのフラグ.
      LOGICAL PFX, PFY             ! 周期境界かを判別する.
                                   ! 周期境界の場合, 配列の始まりと終わりの
                                   ! 要素数が一致していなければならない.
                                   ! つまり, X 方向に境界条件なら,
                                   ! I=1, NX は同じベクトルが入っている必要.
      LOGICAL DTV                  ! FLAG OF VARYING DT (AUTOMATIC).

*-- FUNCTION
      INTEGER IBLKGE      ! DCL FUNCTION
      REAL    ITRP2D      ! INTERPOLATING FUNCTION
                          ! (DEFINED IN BOTTOM OF THIS FILE)

*-- GETTING TMPACK PARAMETERS
      CALL TMRGET( 'STLNDT  ', DT )
      CALL TMIGET( 'STLNGLIM', GLIM )
      CALL TMIGET( 'GRDTHRES', THRES )
      CALL TMLGET( 'PERIODX ', PFX )
      CALL TMLGET( 'PERIODY ', PFY )
      CALL TMLGET( 'FIXEDDT ', DTV )

*-- CHECKING FORWARD OR BACKWARD
      IF(DT.GE.0.0)THEN
         DTPM(1:1)='F'
      ELSE
         DTPM(1:1)='B'
      END IF

*-- GETTING UNDEFINED VALUE
      CALL GLRGET( 'RMISS', UNDEF )

      STEP=NX*NY

*-- INITIALIZING FLAG VALUE FOR DOUBLE COUNTER
      DO 10 J=1,NY
         DO 11 I=1,NX
            CONTINU_FLAG(I,J)=.FALSE.
 11      CONTINUE
 10   CONTINUE

      SGRID=1

*-- CHECKING WHERE A STREAM LINE PASS GRIDS.
*-- WHEN COUNTING NUMBER EXCEEDS "THRES", STOPPING TO CALCULATE A STREAM LINE.
*-- 以下で, 計算最初の格子点をカウント.
      M=IBLKGE( X, NX, INI_X )
      N=IBLKGE( Y, NY, INI_Y )
      INI_M=M
      INI_N=N
      IF(M.EQ.NX)THEN  !* 下の処理で NX+1, NY+1 (領域外) を参照しないための処理.
         M=NX-1
      END IF
      IF(N.EQ.NY)THEN
         N=NY-1
      END IF
      IF((INI_X-X(M)).GE.(X(M+1)-INI_X))THEN  ! M+1 の点を通過
         IF((INI_Y-Y(N)).GE.(Y(N+1)-INI_Y))THEN  ! N+1 の点を通過
            GRID_FLAG(M+1,N+1)=GRID_FLAG(M+1,N+1)+1
            CONTINU_FLAG(M+1,N+1)=.TRUE.
         ELSE  ! N の点を通過
            GRID_FLAG(M+1,N)=GRID_FLAG(M+1,N)+1
            CONTINU_FLAG(M+1,N)=.TRUE.
         END IF
      ELSE  ! M の点を通過
         IF((INI_Y-Y(N)).GE.(Y(N+1)-INI_Y))THEN  ! N+1 の点を通過
            GRID_FLAG(M,N+1)=GRID_FLAG(M,N+1)+1
            CONTINU_FLAG(M,N+1)=.TRUE.
         ELSE  ! N の点を通過
            GRID_FLAG(M,N)=GRID_FLAG(M,N)+1
            CONTINU_FLAG(M,N)=.TRUE.
         END IF
      END IF

*-- SETTING INITIAL POINT FOR CALCULATING STREAM LINE
      TRAJX(1)=INI_X
      TRAJY(1)=INI_Y

*-- STARTING TO CALCULATE STREAM LINES
      DO 100 I=1,STEP-1

*-- IF THERE IS NO THE INITIAL POINT IN THE ANY GRID POINTS,
*-- SEARCHING THE GRID POINT WHICH DOES NOT EXCEED THE INITIAL POINT.
         INTER_P=(/TRAJX(I), TRAJY(I)/)

         DO 99 ID=1,4  !* RK4 scheme start
*-- RESEARCHING THE NEAREST GRID POINTS FOR INTERPOLATING
*-- 内挿用の近傍格子点探索

            M=IBLKGE( X, NX, INTER_P(1) )
            N=IBLKGE( Y, NY, INTER_P(2) )
*-- 以下の処理の都合上, 領域境界上にある場合は 1 格子内側へずらす.
*-- 内挿時には小さい方を参照点として計算するので, この処理をしても
*-- 本質的には影響しない.
            IF(INTER_P(1).EQ.X(NX))THEN
               M=NX-1
            END IF
            IF(INTER_P(2).EQ.Y(NY))THEN
               N=NY-1
            END IF

*-- WHEN A STREAM LINE COME AT BOUNDARIES OF DRAWING REGION OR
*-- A GRID POINT UNDEFINED FOR VECTOR VALUE,
*-- STOPPING TO CALCULATE THE STREAM LINE.
*-- AND, STANDARD OUTPUTTING ITS CONTENT.
*-- 速度場が UNDEF であれば, その時点でそれ以降のデータには UNDEF を
*-- 代入し, UNDEF 影響範囲内を通過した旨を伝える.
*-- また, INTERPO_SEARCH において検索点が領域外の場合もそれ以降計算不可能
*-- である旨を通知.
*-- この処理は RUNGE-KUTTA で計算する 4 回の傾き計算において
*-- 毎回行われる.
*-- ゆえに, 以下の IF 分は 4 回出てくる.
            IF(M.EQ.0.OR.M.EQ.NX)THEN
               IF(PFX.EQV..TRUE.)THEN
                  IF(M.EQ.0)THEN
                     M=NX-1
                  ELSE
                     M=1
                  END IF
               ELSE
                  TRAJN=I
                  GOTO 101
               END IF
            END IF

            IF(N.EQ.0.OR.N.EQ.NY)THEN
               IF(PFY.EQV..TRUE.)THEN
                  IF(N.EQ.0)THEN
                     N=NY-1
                  ELSE
                     N=1
                  END IF
               ELSE
                  TRAJN=I
                  GOTO 101
               END IF
            END IF

            IF(U(M,N).EQ.UNDEF.OR.U(M,N+1).EQ.UNDEF.OR.
     &         U(M+1,N).EQ.UNDEF.OR.U(M+1,N+1).EQ.UNDEF.OR.
     &         V(M,N).EQ.UNDEF.OR.V(M,N+1).EQ.UNDEF.OR.
     &         V(M+1,N).EQ.UNDEF.OR.V(M+1,N+1).EQ.UNDEF)THEN
               TRAJN=I
               GOTO 101
            END IF

            INTER_V(1)=ITRP2D(X(M:M+1),Y(N:N+1),U(M:M+1,N:N+1),INTER_P)
            INTER_V(2)=ITRP2D(X(M:M+1),Y(N:N+1),V(M:M+1,N:N+1),INTER_P)

            IF((INTER_V(1).EQ.0.0).AND.(INTER_V(2).EQ.0.0))THEN
               TRAJN=I
               GOTO 101
            END IF

            IF(.NOT.DTV)THEN
               IF(INTER_V(1).NE.0.0)THEN
                  DTTX=(X(M+1)-X(M))/INTER_V(1)
               ELSE
                  DTTX=0.0
               END IF
               IF(INTER_V(2).NE.0.0)THEN
                  DTTY=(Y(N+1)-Y(N))/INTER_V(2)
               ELSE
                  DTTY=0.0
               END IF
               IF(ABS(DTTX).GT.ABS(DTTY))THEN
                  DT=ABS(DTTY)
               ELSE
                  DT=ABS(DTTX)
               END IF
               IF(DTPM(1:1).EQ.'B')THEN
                  DT=-ABS(DT)
               END IF
            END IF

*-- CALCULATING THE EACH TENDENCY (K,L)
*-- 4 点の傾向を出す (RK4) ので, 4 回繰り返す.
*-- ただし, その時には INTER_P が毎回更新されており, この更新に伴い
*-- 傾向を求めるときのベクトル INTER_V も更新されている.
            IF(ID.EQ.1)THEN
               K1=DT*INTER_V(1)
               L1=DT*INTER_V(2)
*-- 一時的な流跡線の位置を計算
               X1=TRAJX(I)+0.5*K1
               Y1=TRAJY(I)+0.5*L1
               INTER_P=(/X1, Y1/)
            ELSE IF(ID.EQ.2)THEN
               K2=DT*INTER_V(1)
               L2=DT*INTER_V(2)
*-- 一時的な流跡線の位置を計算
               X1=TRAJX(I)+0.5*K2
               Y1=TRAJY(I)+0.5*L2
               INTER_P=(/X1, Y1/)
            ELSE IF(ID.EQ.3)THEN
               K3=DT*INTER_V(1)
               L3=DT*INTER_V(2)
*-- 一時的な流跡線の位置を計算
               X1=TRAJX(I)+K3
               Y1=TRAJY(I)+L3
               INTER_P=(/X1, Y1/)
            ELSE IF(ID.EQ.4)THEN
               K4=DT*INTER_V(1)
               L4=DT*INTER_V(2)
*-- 一時的な流跡線の位置を計算
               X1=TRAJX(I)+K4
               Y1=TRAJY(I)+L4
*               INTER_P=(/X1, Y1/)  !* ここは内挿計算に使わないので求めない.
            END IF

*-- ここまでで 1 サイクル
 99      CONTINUE

         X1=TRAJX(I)+(1.0/6.0)*(K1+2.0*K2+2.0*K3+K4)
         Y1=TRAJY(I)+(1.0/6.0)*(L1+2.0*L2+2.0*L3+L4)

*-- CHECKING THAT THERE IS A CALCULATING STREAM LINE IN THE DRAWING REGION.
*-- 計算した TRAJ が領域内に存在しているか確認
         IF(X1.LT.X(1).OR.X1.GT.X(NX))THEN
            IF(PFX.EQV..TRUE.)THEN  !* 周期境界の場合, 流線の位置を変化させ 1 つ前の流線位置に undef を代入する.
               IF(X1.LT.X(1))THEN
                  X1=X(NX)+(X1-X(1))
               ELSE
                  X1=X(1)+(X1-X(NX))
               END IF
               TRAJX(I)=UNDEF
               TRAJY(I)=UNDEF
            ELSE
               IF(X1.LT.X(1))THEN
                  TRAJX(I+1)=X(1)
                  TRAJY(I+1)=TRAJY(I)+((X(1)-TRAJX(I))/(X1-TRAJX(I)))  
     &                                *(Y1-TRAJY(I))
               ELSE
                  TRAJX(I+1)=X(NX)
                  TRAJY(I+1)=TRAJY(I)+((X(NX)-TRAJX(I))/(X1-TRAJX(I)))  
     &                                *(Y1-TRAJY(I))
               END IF
               TRAJN=I+1
               GOTO 101
            END IF
         END IF

         IF(Y1.LT.Y(1).OR.Y1.GT.Y(NY))THEN
            IF(PFY.EQV..TRUE.)THEN
               IF(Y1.LT.Y(1))THEN
                  Y1=Y(NY)+(Y1-Y(1))
               ELSE
                  Y1=Y(1)+(Y1-Y(NY))
               END IF
               TRAJX(I)=UNDEF
               TRAJY(I)=UNDEF
            ELSE
               IF(Y1.LT.Y(1))THEN
                  TRAJY(I+1)=Y(1)
                  TRAJX(I+1)=TRAJX(I)+((Y(1)-TRAJY(I))/(Y1-TRAJY(I)))  
     &                                *(X1-TRAJX(I))
               ELSE
                  TRAJY(I+1)=Y(NY)
                  TRAJX(I+1)=TRAJX(I)+((Y(NY)-TRAJY(I))/(Y1-TRAJY(I)))  
     &                                *(X1-TRAJX(I))
               END IF
               TRAJN=I+1
               GOTO 101
            END IF
         END IF

*-- CHECKING WHERE A STREAM LINE PASS GRIDS.
*-- WHEN COUNTING NUMBER EXCEEDS "THRES", STOPPING TO CALCULATE A STREAM LINE.
*-- 以下で, どの格子点を通過したかをカウントする.
*-- また, 通過数が THRES を越えていると, GOTO 101 で終了.
         IF((X1-X(M)).GE.(X(M+1)-X1))THEN  ! M+1 の点を通過
            IF((Y1-Y(N)).GE.(Y(N+1)-Y1))THEN  ! N+1 の点を通過
               IF(CONTINU_FLAG(M+1,N+1).EQV..FALSE.)THEN
                  GRID_FLAG(M+1,N+1)=GRID_FLAG(M+1,N+1)+1
                  CONTINU_FLAG=.FALSE.
                  CONTINU_FLAG(M+1,N+1)=.TRUE.
                  SGRID=1
                  IF(GRID_FLAG(M+1,N+1).GT.THRES)THEN
                     IF(M+1.NE.INI_M.OR.N+1.NE.INI_N)THEN  ! 初期点に戻ってくる場合の処理 (戻ってこないとこの IF 文に入る)
                        TRAJN=I
                        GOTO 101
                     ELSE
                        TRAJX(I+1)=X1
                        TRAJY(I+1)=Y1
                        TRAJN=I+1
                        GOTO 101
                     END IF
                  END IF
               ELSE  ! COUNTING THE EXISTING OF THE SAME GRID
                  SGRID=SGRID+1
               END IF
            ELSE  ! N の点を通過
               IF(CONTINU_FLAG(M+1,N).EQV..FALSE.)THEN
                  GRID_FLAG(M+1,N)=GRID_FLAG(M+1,N)+1
                  CONTINU_FLAG=.FALSE.
                  CONTINU_FLAG(M+1,N)=.TRUE.
                  SGRID=1
                  IF(GRID_FLAG(M+1,N).GT.THRES)THEN
                     IF(M+1.NE.INI_M.OR.N.NE.INI_N)THEN  ! 初期点に戻ってくる場合の処理 (戻ってこないとこの IF 文に入る)
                        TRAJN=I
                        GOTO 101
                     ELSE
                        TRAJX(I+1)=X1
                        TRAJY(I+1)=Y1
                        TRAJN=I+1
                        GOTO 101
                     END IF
                  END IF
               ELSE  ! COUNTING THE EXISTING OF THE SAME GRID
                  SGRID=SGRID+1
               END IF
            END IF
         ELSE  ! M の点を通過
            IF((Y1-Y(N)).GE.(Y(N+1)-Y1))THEN  ! N+1 の点を通過
               IF(CONTINU_FLAG(M,N+1).EQV..FALSE.)THEN
                  GRID_FLAG(M,N+1)=GRID_FLAG(M,N+1)+1
                  CONTINU_FLAG=.FALSE.
                  CONTINU_FLAG(M,N+1)=.TRUE.
                  SGRID=1
                  IF(GRID_FLAG(M,N+1).GT.THRES)THEN
                     IF(M.NE.INI_M.OR.N+1.NE.INI_N)THEN  ! 初期点に戻ってくる場合の処理 (戻ってこないとこの IF 文に入る)
                        TRAJN=I
                        GOTO 101
                     ELSE
                        TRAJX(I+1)=X1
                        TRAJY(I+1)=Y1
                        TRAJN=I+1
                        GOTO 101
                     END IF
                  END IF
               ELSE  ! COUNTING THE EXISTING OF THE SAME GRID
                  SGRID=SGRID+1
               END IF
            ELSE  ! N の点を通過
               IF(CONTINU_FLAG(M,N).EQV..FALSE.)THEN
                  GRID_FLAG(M,N)=GRID_FLAG(M,N)+1
                  CONTINU_FLAG=.FALSE.
                  CONTINU_FLAG(M,N)=.TRUE.
                  SGRID=1
                  IF(GRID_FLAG(M,N).GT.THRES)THEN
                     IF(M.NE.INI_M.OR.N.NE.INI_N)THEN  ! 初期点に戻ってくる場合の処理 (戻ってこないとこの IF 文に入る)
                        TRAJN=I
                        GOTO 101
                     ELSE
                        TRAJX(I+1)=X1
                        TRAJY(I+1)=Y1
                        TRAJN=I+1
                        GOTO 101
                     END IF
                  END IF
               ELSE  ! COUNTING THE EXISTING OF THE SAME GRID
                  SGRID=SGRID+1
               END IF
            END IF
         END IF

         TRAJX(I+1)=X1
         TRAJY(I+1)=Y1
         TRAJN=I+1

*-- 速度場がゼロ以外で Runge-Kutta 計算上たまたまゼロになることがある.
         IF(TRAJX(I).EQ.TRAJX(I+1).AND.TRAJY(I).EQ.TRAJY(I+1))THEN
            GO TO 101
         ELSE IF(ABS(TRAJX(I+1)-TRAJX(I)).LT.(X(M+1)-X(M))/REAL(GLIM)
     &           .AND.
     &           ABS(TRAJY(I+1)-TRAJY(I)).LT.(Y(N+1)-Y(N))/REAL(GLIM)
     &           )THEN
            CALL MSGDMP( 'M', 'TMSLCL', 'DETECTED STREAM .LT. GLIM.' ) 
            GO TO 101
         ELSE IF(SGRID>GLIM)THEN  ! PREVENTING FROM OSCILLATION IN A GRID.
            CALL MSGDMP( 'M', 'TMSLCL', 'DETECTED EXISTING SAME GRID.' )
            TRAJN=TRAJN-GLIM
            GO TO 101
         END IF

C      WRITE(*,*) "V CHECK", TRAJX(I:I+1), TRAJY(I:I+1), I 
 100  CONTINUE

 101  TRAJN=TRAJN

      END SUBROUTINE TMSLCL

*---------------------------------
* PRIVATE ROUTINE (INTERPOLATING FUNCTION)
*---------------------------------

      REAL FUNCTION ITRP2D( X, Y, Z, P )
      IMPLICIT NONE
      REAL X(2)
      REAL Y(2)
      REAL Z(2,2)
      REAL P(2)
      REAL VALX(2)
      REAL V
    
      ! Y(1) での X 方向の内挿点での値
      VALX(1)=Z(1,1)+(P(1)-X(1))*(Z(2,1)-Z(1,1))/(X(2)-X(1))
    
      ! Y(2) での X 方向の内挿点での値
      VALX(2)=Z(1,2)+(P(1)-X(1))*(Z(2,2)-Z(1,2))/(X(2)-X(1))
    
      ! X の内挿点からの Y 方向の内挿点での値(これが求める内挿点)
      V=VALX(1)+(P(2)-Y(1))*(VALX(2)-VALX(1))/(Y(2)-Y(1))
    
      ITRP2D=V

      RETURN
      END FUNCTION ITRP2D

*----------------------------------------------------------------------

      
