      MODULE ARRAY
        IMPLICIT NONE
        INTEGER(4)                            :: NA
        INTEGER(8),DIMENSION (:),ALLOCATABLE  :: IN,P1,P2
        REAL(4),DIMENSION (:),ALLOCATABLE     :: D,F0
      END MODULE ARRAY
      PROGRAM COEFR
      USE ARRAY
      IMPLICIT NONE
      INTEGER,PARAMETER          :: LMAX=20000000
      INTEGER(4)                 :: NS
      INTEGER(4)                 :: NR
      INTEGER(4)                 :: NB
      INTEGER(4)                 :: NOUT
      INTEGER(4)                 :: I
      CHARACTER (LEN=104)        :: INFILE(2),OTFILE(3)
      LOGICAL                    :: TF
      CHARACTER (LEN=1)          :: YN,YN1,YN2
      INTEGER(4)                 :: ERR
      CHARACTER                  :: DATE*8,TIME*10
      INTERFACE
        SUBROUTINE INSORT (JG,NA,NB,IN,P1,P2,NR,L,D,F0,INFILE,OTFILE,*)
          INTEGER(4),INTENT(IN)                   :: JG
          INTEGER(4),INTENT(IN)                   :: NA
          INTEGER(4),INTENT(INOUT)                :: NB
          INTEGER(8),INTENT(INOUT),DIMENSION (:)  :: IN,P1,P2
          INTEGER(4),INTENT(IN)                   :: NR
          INTEGER(4),INTENT(IN)                   :: L
          REAL(4),INTENT(INOUT),DIMENSION (:)     :: D
          REAL(4),INTENT (INOUT),DIMENSION (0:)   :: F0
          CHARACTER (LEN=104),INTENT(IN)          :: INFILE(2),OTFILE(3)
        END SUBROUTINE INSORT
        SUBROUTINE INBCOE (NA,NS,P1,P2,D,F0)
          INTEGER(4)                              :: NA,NS
          INTEGER(8),INTENT(INOUT),DIMENSION (:)  :: P1,P2
          REAL(4),INTENT(INOUT),DIMENSION (:)     :: D
          REAL(4),INTENT(INOUT),DIMENSION (0:)    :: F0
        END SUBROUTINE INBCOE
        SUBROUTINE RQSORT (I,J,N,IN,P1,P2,F)
          INTEGER(4),INTENT(IN)                   :: I,J,N
          INTEGER(8),INTENT(INOUT),DIMENSION (:)  :: IN,P1,P2
          REAL(4),INTENT(INOUT),DIMENSION (:)     :: F
        ENDSUBROUTINE RQSORT
        INTEGER(4) FUNCTION IFIND (ID,N,IN)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                   :: N
          INTEGER(8),INTENT(IN)                   :: ID
          INTEGER(8),INTENT(IN),DIMENSION (:)     :: IN
        END FUNCTION IFIND
      END INTERFACE
   10 WRITE (6,'(A)') '血統ファイル名を入力して下さい.'
      READ (5,*) INFILE(1)
      INQUIRE (FILE=INFILE(1),EXIST=TF)
      IF (.NOT.TF) THEN
        WRITE (6,'(A)') '指定したファイルが存在しません.'
        GO TO 10
      END IF
      NS = 0
      WRITE (6,'(2A)')  '個体の親番号が個体番号の列にない場合の',
     &                  '処理方法を選択して下さい(1-2).'
      WRITE (6,'(3X,2A)') '1. 共通な親番号を持つ個体間の血縁関係を',
     &                    '無視する.'
      WRITE (6,'(3X,2A)') '2. 共通な親番号を持つ個体間の血縁関係を',
     &                    '考慮する.'
   30 READ (5,*) NR
      IF (NR.NE.1.AND.NR.NE.2) THEN
        WRITE (6,'(A)') '１または２のいずれかを選択して下さい.'
        GO TO 30
      END IF
      WRITE (6,'(2A)')  '近交係数のファイルへの出力形態を',
     &                  '選択して下さい(1-3).'
      WRITE (6,'(3X,A)') '1. 血統ファイルの個体順に出力する.'
      WRITE (6,'(3X,A)') '2. ソートした個体番号順に出力する.'
      WRITE (6,'(3X,A)') '3. １および２のファイルを出力する.'
   40 READ (5,*) NOUT
      IF (NOUT.NE.1.AND.NOUT.NE.2.AND.NOUT.NE.3) THEN
        WRITE (6,'(A)') '１,２または３のいずれかを選択して下さい.'
        GO TO 40
      END IF
      IF (NOUT.EQ.3) THEN
        WRITE (6,'(2A)') '２のファイルは１のファイル名に拡張子',
     &                   '''.SORT''を付けて出力されます.'
      END IF
   50 WRITE (6,'(A)') '近交係数の出力ファイル名を入力して下さい.'
      READ (5,'(A)') OTFILE(1)
      INQUIRE (FILE=OTFILE(1),EXIST=TF)
      IF (TF) THEN
        IF (OTFILE(1).EQ.INFILE(1)) THEN
          WRITE (6,'(A)') '指定したファイルは血統ファイルと同一名です.'
          GO TO 50
        ELSE
          WRITE (6,'(A)') '指定したファイルはすでに存在します.'
          WRITE (6,'(A)') 'ファイルを上書きしてもよろしいですか(Y/N).'
          READ (5,*) YN1
          IF (YN1.NE.'Y'.AND.YN1.NE.'y') THEN
            GO TO 50
          END IF
        END IF
      END IF
      IF (NOUT.EQ.3) THEN
        DO I=50,1,-1
          IF (OTFILE(1)(I:I).NE.' ') THEN
            OTFILE(2) = OTFILE(1)(1:I)//'.SORT'
            EXIT
          END IF
        END DO
        INQUIRE (FILE=OTFILE(2),EXIST=TF)
        IF (TF) THEN
          WRITE (6,'(2A)') '指定した拡張子付きファイルはすでに',
     &                     '存在します.'
          WRITE (6,'(A)') 'ファイルを上書きしてもよろしいですか(Y/N).'
          READ (5,*) YN
          IF (YN.NE.'Y'.AND.YN.NE.'y') THEN
            GO TO 50
          ELSE
            OPEN (22,FILE=OTFILE(2))
            CLOSE (22,STATUS='DELETE')
          END IF
        END IF
      END IF
      WRITE (6,'(A)') '血縁係数を算出しますか(Y/N).'
      READ (5,*) YN2
      IF (YN2.NE.'Y'.AND.YN2.NE.'y')     GO TO 100
   60 WRITE (6,'(2A)') '血縁ファイル名を入力して下さい.'
      READ (5,*) INFILE(2)
      INQUIRE (FILE=INFILE(2),EXIST=TF)
      IF (.NOT.TF) THEN
        WRITE (6,'(A)') '指定したファイルが存在しません.'
        GO TO 60
      ELSE
        IF (INFILE(2).EQ.INFILE(1)) THEN
          WRITE (6,'(A)') '指定したファイルは血統ファイルと同一名です.'
          GO TO 60
        END IF
        IF (OTFILE(1).EQ.INFILE(2)) THEN
          WRITE (6,'(2A)') '近交係数の出力ファイルが血縁ファイルと',
     &                     '同一名です.'
          GO TO 50
        END IF
      END IF
   70 WRITE (6,'(A)') '出力先のファイル名を入力して下さい.'
      READ (5,'(A)') OTFILE(3)
      INQUIRE (FILE=OTFILE(3),EXIST=TF)
      IF (TF) THEN
        IF (OTFILE(3).EQ.INFILE(1)) THEN
          WRITE (6,'(A)') '指定したファイルは血統ファイルと同一名です.'
          GO TO 70
        ELSE IF (OTFILE(3).EQ.INFILE(2)) THEN
          WRITE (6,'(2A)') '指定したファイルは血縁ファイルと同一名です.'
          GO TO 70
        ELSE IF (OTFILE(3).EQ.OTFILE(1)) THEN
          WRITE (6,'(2A)') '指定したファイルは近交係数の出力ファイルと',
     &                     '同一名です.'
          GO TO 70
        ELSE
          WRITE (6,'(A)') '指定したファイルはすでに存在します.'
          WRITE (6,'(A)') 'ファイルを上書きしてもよろしいですか(Y/N).'
          READ (5,*) YN
          IF (YN.NE.'Y'.AND.YN.NE.'y') THEN
            GO TO 70
          ELSE
            OPEN (21,FILE=OTFILE(3))
            CLOSE (21,STATUS='DELETE')
          END IF
        END IF
      END IF
  100 OPEN (11,FILE=INFILE(1))
        F0 = 0.0
        DO I=1,NS
          READ (11,*)
        END DO
        NA = NS+1
        DO
          READ (11,*,END=120)
          NA = NA+1
          IF (NA.GT.LMAX) THEN
            WRITE (6,*)
            WRITE (6,'(A)') '個体数が制限を越えています.'
            WRITE (6,'(A)') '制限を拡張した後，コンパイルして下さい.'
            DO
              READ (11,*,END=110)
              NA = NA+1
            END DO
  110       WRITE (6,'(A,I10)') '血統ファイル中に存在する個体数  =',NA
            WRITE (6,*)
            STOP
          END IF
        END DO
  120 CLOSE (11)
      NA = NA-1
      WRITE (6,*)
      WRITE (6,'(A,I10)') '全個体数                =',NA
      WRITE (6,*)
      ALLOCATE (IN(NA+1),P1(NA+1),P2(NA+1),D(NA+1),F0(0:NA+1),STAT=ERR)
      IF (ERR.GT.0) THEN
        WRITE (6,'(2A)') '現在使用中のコンピュータでは個体数が',
     &                   '多過ぎる可能性があります.'
        WRITE (6,'(2A)') 'より性能の高いコンピュータを使用するか',
     &                   '個体数を少なくする必要があります.'
        GO TO 900
      END IF
      OPEN (11,FILE=INFILE(1))
        F0(0:NA+1) = 0.0
        DO I=NS+1,NA
          READ (11,*) IN(I),P1(I),P2(I)
        END DO
      CLOSE (11)
      CALL DATE_AND_TIME (DATE,TIME)
      WRITE (6,'(6A,2X,A)') TIME(1:2),'時',TIME(3:4),'分',TIME(5:),
     &                      '秒','個体番号ソート開始'
      CALL INSORT (1,NA,NB,IN,P1,P2,NR,LMAX,D,F0,INFILE,OTFILE,*130)
      IF (NR.EQ.2.AND.NB.GT.0) THEN
        DEALLOCATE (IN,P1,P2,D,F0,STAT=ERR)
        ALLOCATE (IN(NA+NB+1),P1(NA+NB+1),P2(NA+NB+1),
     &                            D(NA+NB+1),F0(0:NA+NB+1),STAT=ERR)
        IF (ERR.GT.0) THEN
          WRITE (6,'(2A)') '現在使用中のコンピュータでは個体数が',
     &                     '多過ぎる可能性があります.'
          WRITE (6,'(2A)') 'より性能の高いコンピュータを使用するか',
     &                     '個体数を少なくする必要があります.'
          GO TO 900
        END IF
        OPEN (11,FILE=INFILE(1))
          F0(0) = 0.0
          DO I=1,NA
            READ (11,*) IN(I),P1(I),P2(I)
          END DO
        CLOSE (11)
        CALL INSORT (2,NA,NB,IN,P1,P2,NR,LMAX,D,F0,INFILE,OTFILE,
     &                                                         *130)
      END IF
  130 CALL DATE_AND_TIME (DATE,TIME)
      WRITE (6,'(6A,2X,A)') TIME(1:2),'時',TIME(3:4),'分',TIME(5:),
     &                      '秒','近交係数の算出開始'
      CALL INBCOE (NA+NB,NS,P1,P2,D,F0)
      IF (YN2.NE.'Y'.AND.YN2.NE.'y')     GO TO 200
      CALL DATE_AND_TIME (DATE,TIME)
      WRITE (6,'(6A,2X,A)') TIME(1:2),'時',TIME(3:4),'分',TIME(5:),
     &                      '秒','血縁係数の算出開始'
      CALL INSORT (3,NA,NB,IN,P1,P2,NR,LMAX,D,F0,INFILE,OTFILE,*200)
  200 CALL DATE_AND_TIME (DATE,TIME)
      WRITE (6,'(6A,2X,A)') TIME(1:2),'時',TIME(3:4),'分',TIME(5:),
     &                      '秒','近交係数の出力開始'
      CALL RQSORT (NB+1,NB+NA,1,IN,P1,P2,F0(1:))
      OPEN (11,FILE=INFILE(1))
        DO I=1,NB
          P1(I) = 0 ; P2(I) = 0 ; F0(I) = 0.0
        END DO
        DO I=1,NA
          READ (11,*) IN(NB+I),P1(NB+I),P2(NB+I)
        END DO
      CLOSE (11)
      IF (NOUT.EQ.2) THEN
        CALL RQSORT (1,NB+NA,3,IN,P1,P2,F0(1:))
      END IF
      IF (YN1.EQ.'Y'.OR.YN1.EQ.'y') THEN
        OPEN (21,FILE=OTFILE(1))
        CLOSE (21,STATUS='DELETE')
      ELSE IF (YN2.EQ.'Y'.OR.YN2.EQ.'y') THEN
        OPEN (21,FILE=OTFILE(1))
        CLOSE (21,STATUS='DELETE')
      END IF
      OPEN (21,FILE=OTFILE(1),STATUS='NEW')
        DO I=1,NA+NB
          WRITE (21,2000) IN(I),P1(I),P2(I),F0(I)
        END DO
      CLOSE (21)
      IF (NOUT.EQ.3) THEN
        CALL RQSORT (1,NB+NA,3,IN,P1,P2,F0(1:))
        OPEN (22,FILE=OTFILE(2),STATUS='NEW')
          DO I=1,NA+NB
            WRITE (22,2000) IN(I),P1(I),P2(I),F0(I)
          END DO
        CLOSE (22)
      END IF
      CALL DATE_AND_TIME (DATE,TIME)
      WRITE (6,'(6A,2X,A)') TIME(1:2),'時',TIME(3:4),'分',TIME(5:),
     &                      '秒','近交係数の演算終了'
  900 STOP
 2000 FORMAT (3I20,F20.6)
      END PROGRAM COEFR
      SUBROUTINE INSORT (JG,NA,NB,IN,P1,P2,NR,LMAX,D,F0,INFILE,OTFILE,*)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                   :: JG
      INTEGER(4),INTENT(IN)                   :: NA
      INTEGER(4),INTENT(INOUT)                :: NB
      INTEGER(8),INTENT(INOUT),DIMENSION (:)  :: IN,P1,P2
      INTEGER(4),INTENT(IN)                   :: NR
      INTEGER(4),INTENT(IN)                   :: LMAX
      REAL(4),INTENT(INOUT),DIMENSION (:)     :: D
      REAL(4),INTENT(INOUT),DIMENSION (0:)    :: F0
      CHARACTER (LEN=104),INTENT(IN)          :: INFILE(2),OTFILE(3)
      INTEGER(4)                      :: NT
      INTEGER(4)                      :: I,J
      INTEGER(4)                      :: N
      INTEGER(8),DIMENSION (NA+NB)    :: WK,WI,WP
      INTEGER(8),DIMENSION (0:NA+NB)  :: WK0
      INTEGER(4)                      :: WNB,IWK,WN
      INTEGER(4)                      :: WB
      INTEGER(8),DIMENSION (NA*2)     :: B
      INTEGER(4)                      :: NS
      INTEGER(4)                      :: NA1,NA2,PT1,PT2
      INTEGER(8)                      :: ID,ID1,ID2
      INTEGER(8)                      :: WID1,WID2
      REAL(4)                         :: FR
      REAL(4)                         :: F01,F02
      INTERFACE
        SUBROUTINE IQSORT (I,J,K,A,B)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                    :: I,J,K
          INTEGER(8),INTENT(INOUT),DIMENSION (:)   :: A,B
        END SUBROUTINE IQSORT
        SUBROUTINE IQSORT5 (I,J,K,A,B,C,D,E)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                    :: I,J,K
          INTEGER(8),INTENT(INOUT),DIMENSION (:)   :: A,B,C,D,E
        END SUBROUTINE IQSORT5
        SUBROUTINE RLSCOE (NS,P1,P2,D,F0)
          INTEGER(4),INTENT(IN)                    :: NS
          INTEGER(8),INTENT(INOUT),DIMENSION (:)   :: P1,P2
          REAL(4),INTENT(INOUT),DIMENSION (:)      :: D
          REAL(4),INTENT(INOUT),DIMENSION (0:)     :: F0
        END SUBROUTINE RLSCOE
        INTEGER(4) FUNCTION IFIND (ID,N,IN)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)               :: N
          INTEGER(8),INTENT(IN)               :: ID
          INTEGER(8),INTENT(IN),DIMENSION (:) :: IN
        END FUNCTION IFIND
      END INTERFACE
      IF (JG.EQ.3)     GO TO 100
      NB = 0
      WK0(0) = 0
      DO I=1,NA
        WK0(I) = I
        WK(I) = I
      END DO
      DO I=1,NA
        IF (IN(I).EQ.P1(I).OR.IN(I).EQ.P2(I)) THEN
          WRITE (6,*)
          WRITE (6,'(A)') '個体番号と親番号が同一です.'
          WRITE (6,*)
          WRITE (6,'(A,I18)') '当該個体番号の行  =',I
          WRITE (6,'(A,I18)') '個体番号＆親番号  =',IN(I)
          WRITE (6,*)
          STOP
        END IF
      END DO
      CALL IQSORT (1,NA,2,IN,WK)
      DO I=1,NA-1
        IF (IN(I).LE.0) THEN
          WRITE (6,*)
          WRITE (6,'(A)') '個体番号＝０以下の個体が存在します.'
          WRITE (6,*)
          WRITE (6,'(A,I10)') '当該個体番号の行  =',WK(I)
          WRITE (6,*)
          STOP
        END IF
        IF (IN(I).EQ.IN(I+1)) THEN
          WRITE (6,*)
          WRITE (6,'(A,I18,A)') '同一個体番号(',IN(I),' )が存在します.'
          WRITE (6,*)
          IF (WK(I).GT.WK(I+1)) THEN
            IWK = WK(I) ; WK(I) = WK(I+1) ; WK(I+1) = IWK
          END IF
          WRITE (6,'(A,I10)') '当該個体番号の行  =',WK(I)
          WRITE (6,'(A,I10)') '                   ',WK(I+1)
          WRITE (6,*)
          STOP
        END IF
      END DO
      IF (NR.EQ.1) THEN
        N = NA
        WK0(1:) = WK
        DO I=1,N
          P1(I) = WK0(IFIND(P1(I),N,IN))
          P2(I) = WK0(IFIND(P2(I),N,IN))
        END DO
      ELSE
        DO I=1,NA
          IWK = IFIND(P1(I),NA,IN)
          IF (IWK.EQ.0.AND.P1(I).NE.0) THEN
            NB = NB+1
            B(NB) = P1(I)
          END IF
          IWK = IFIND(P2(I),NA,IN)
          IF (IWK.EQ.0.AND.P2(I).NE.0) THEN
            NB = NB+1
            B(NB) = P2(I)
          END IF
        END DO
        CALL IQSORT (1,NB,1,B,WI)
        WI = 0
        WB = B(1)
        WNB = 1
        IWK = 0
        DO I=2,NB
          IF (B(I).EQ.WB) THEN
            WNB = WNB+1
            IF (I.EQ.NB) THEN
              IWK = IWK+1
              B(IWK) = WB
            END IF
          ELSE
            IF (WNB.GE.2) THEN
              IWK = IWK+1
              B(IWK) = WB
              WNB = 1
            END IF
            WB = B(I)
          END IF
        END DO
        NB = IWK
        B(NB+1:) = 0
        IF (NA+NB.GT.LMAX) THEN
          WRITE (6,*)
          WRITE (6,'(A)') '個体数が制限を越えています.'
          WRITE (6,'(A)') '制限を拡張した後，コンパイルし直して下さい.'
          WRITE (6,*)
          WRITE (6,'(A,I10)') '新たに加わった共通な親の数  =',NB
          WRITE (6,'(A,I10)') '血統ファイル中の個体数      =',NA
          WRITE (6,'(A,I10)') '全個体数                    =',NB+NA
          WRITE (6,*)
          STOP
        END IF
      END IF
      IF (NB.GT.0.AND.JG.NE.2)     RETURN
      IF (NR.EQ.2) THEN
        CALL IQSORT (1,NA,2,WK,IN)
        DO I=NA,1,-1
          IN(I+NB) = IN(I)
          P1(I+NB) = P1(I)
          P2(I+NB) = P2(I)
          WK(I+NB) = I+NB
        END DO
        DO I=1,NB
          IN(I) = B(I)
          P1(I) = 0
          P2(I) = 0
          WK(I) = I
        END DO
        N = NB+NA
        CALL IQSORT (1,N,2,IN,WK)
        WK0(1:) = WK
        DO I=1,N
          P1(I) = WK0(IFIND(P1(I),N,IN))
          P2(I) = WK0(IFIND(P2(I),N,IN))
        END DO
        WRITE (6,*)
        WRITE (6,'(A,I10)') '新たに加わった個体数    =',NB
        WRITE (6,*)
      END IF
      CALL IQSORT (1,N,2,WK,IN)
      DO I=1,N
        IN(I) = I
      END DO
      IWK = 0
      DO I=1,N
        IF (P1(I).GE.I.OR.P2(I).GE.I) THEN
          IWK = IWK+1
          IF (IWK.EQ.1) THEN
            WRITE (6,*)
            WRITE (6,'(2A)') '個体番号が生年月日順に並んでいないため, ',
     &                       '番号の並べ換えを行います.'
            WRITE (6,*)
          END IF
          NS = I
          DO J=NS,N
            WI(J) = J
            WP(J) = MAX(P1(J),P2(J))
          END DO
          CALL IQSORT5 (NS,N,5,WP,IN,P1,P2,WI)
          DO J=1,N
            WK(J) = J
          END DO
          CALL IQSORT (NS,N,2,WI,WK)
          WK0(1:) = WK
          DO J=NS,N
            P1(J) = WK0(P1(J))
            P2(J) = WK0(P2(J))
          END DO
        END IF
      END DO
      DO I=1,NB
        IN(I) = B(I)
      END DO
      DO I=1,N
        WN = 0
        IWK = I
        DO
          IF (P1(IWK).EQ.0) THEN
            EXIT
          ELSE
            WN = WN+1
            IWK = P1(IWK)
            IF (WN.GE.1000) THEN
              WN = 0
              IWK = I
              WK = 0
   10         WN = WN+1
              WK(WN) = IWK
              IWK = P1(IWK)
              DO J=1,WN
                IF (IWK.EQ.WK(J)) THEN
                  WRITE (6,*)
                  WRITE (6,'(A,I10,2A)') '個体番号の行(',IWK-NB,' )',
     &                  'の父番号が当該個体の子孫の番号に存在します.'
                  WRITE (6,*)
                  STOP
                END IF
              END DO
              GO TO 10
            END IF
          END IF
        END DO
        WN = 0
        IWK = I
        DO
          IF (P2(IWK).EQ.0) THEN
            EXIT
          ELSE
            WN = WN+1
            IWK = P2(IWK)
            IF (WN.GE.1000) THEN
              WN = 0
              IWK = I
              WK = 0
   20         WN = WN+1
              WK(WN) = IWK
              IWK = P2(IWK)
              DO J=1,WN
                IF (IWK.EQ.WK(J)) THEN
                  WRITE (6,*)
                  WRITE (6,'(A,I10,2A)') '個体番号の行(',IWK-NB,' )',
     &                  'の母番号が当該個体の子孫の番号に存在します.'
                  WRITE (6,*)
                  STOP
                END IF
              END DO
              GO TO 20
            END IF
          END IF
        END DO
      END DO
      RETURN
  100 NT = NB+NA
      OPEN (11,FILE=INFILE(1))
        DO I=1,NA
          READ (11,*) WI(I)
        END DO
      CLOSE (11)
      DO I=1,NA
        WK(I) = I
      END DO
      CALL IQSORT (1,NA,2,WI,WK)
      OPEN (12,FILE=INFILE(2))
      OPEN (21,FILE=OTFILE(1))
        READ (12,*) NA1,NA2
        DO I=1,NA1+NA2
          READ (12,*) ID
          IWK = WK(IFIND(ID,NA,WI))
          IF (IWK.EQ.0) THEN
            WRITE (6,'(2A)') '血縁ファイルの個体番号が血統ファイルに',
     &                       '存在しません.'
            WRITE (6,'(A,I10)') '当該個体番号の行  =',I
            WRITE (6,*)
            WRITE (6,'(A)') '血縁係数の算出を中止します.'
            RETURN 1
          ELSE
            WRITE (21,*) NB+IWK
          END IF
        END DO
      CLOSE (12)
      CLOSE (21)
      DO I=1,NA
        WK(I) = I
      END DO
      CALL IQSORT (1,NT,2,IN,WK)
      OPEN (11,FILE=OTFILE(1))
      OPEN (12,FILE=INFILE(2))
      OPEN (23,FILE=OTFILE(3))
      PT1 = 0; PT2 = NA1
  150 REWIND (11); REWIND (12)
      READ (12,*)
      DO I=1,PT1
        READ (11,*)
        READ (12,*)
      END DO
      READ (11,*) ID1
      READ (12,*) WID1
      P1(NT+1) = WK(IFIND(ID1,NT,IN))
      F01 = F0(P1(NT+1))
      ID = P1(NT+1)
      REWIND (11); REWIND (12)
      READ (12,*)
      DO I=1,PT2
        READ (11,*)
        READ (12,*)
      END DO
      DO I=1,NA2
        READ (11,*) ID2
        READ (12,*) WID2
        P2(NT+1) = WK(IFIND(ID2,NT,IN))
        F02 = F0(P2(NT+1))
        CALL RLSCOE (NT,P1,P2,D,F0)
        FR = (F0(NT+1)*2)/
     &       (SQRT(1+F01)*SQRT(1+F02))
        WRITE (23,2100) WID1,WID2,FR
        P1(NT+1) = ID
      END DO
      IF (PT1+1.NE.NA1) THEN
        PT1 = PT1+1
        GO TO 150
      END IF
      CLOSE (11)
      CLOSE (12)
      CLOSE (23)
      CALL IQSORT (1,NT,2,WK,IN)
      RETURN
 2100 FORMAT (2I20,20X,F20.6)
      END SUBROUTINE INSORT
      RECURSIVE SUBROUTINE IQSORT (I,J,N,A,B)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: I,J,N
      INTEGER(8),INTENT(INOUT),DIMENSION (:) :: A,B
      INTEGER(4) :: IPV,KPV
      INTEGER(8) :: PV
      INTEGER(4) :: S0,S1,S2,K,H
      INTEGER(4) :: LET,RIT
      INTEGER(8) :: WKA,WKB
      S0 = J-I
      S1 = (S0+1)/2
      IF (MOD(S0,2).EQ.0) THEN
        S2 = S1
      ELSE
        S2 = S1-1
      END IF
      H = 1
      K = I+S1
      DO
        IF (A(I).NE.A(K).OR.H.GE.S0) THEN
          EXIT
        ELSE
          IF (MOD(H,2).EQ.1) THEN
            K = K+S2
          ELSE
            K = K+S1
          END IF
          IF (K.GT.J)     K = K-J+I-1
          H = H+1
        END IF
      END DO
      IF (A(I).GT.A(K)) THEN
        IPV = I
      ELSE  IF (A(I).LT.A(K)) THEN
        IPV = K
      ELSE
        IPV = 0
      END IF
      IF (IPV.EQ.0) RETURN
      PV = A(IPV)
      LET = I
      RIT = J
      DO
        IF (A(LET).GE.PV) THEN
          EXIT
        ELSE
          LET = LET+1
        END IF
      END DO
      DO
        IF (A(RIT).LT.PV) THEN
          EXIT
        ELSE
          RIT = RIT-1
        END IF
      END DO
      DO
        IF (LET.GT.RIT) THEN
          EXIT
        ELSE
          WKA = A(LET) ; A(LET) = A(RIT) ; A(RIT) = WKA
          IF (N.EQ.2) THEN
            WKB = B(LET) ; B(LET) = B(RIT) ; B(RIT) = WKB
          END IF
          LET = LET+1     ; RIT = RIT-1
          DO
            IF (A(LET).GE.PV) THEN
              EXIT
            ELSE
              LET = LET+1
            END IF
          END DO
          DO
            IF (A(RIT).LT.PV) THEN
              EXIT
            ELSE
              RIT = RIT-1
            END IF
          END DO
        END IF
      END DO
      KPV = LET
      CALL IQSORT (I,KPV-1,N,A,B)
      CALL IQSORT (KPV,J,N,A,B)
      RETURN
      END SUBROUTINE IQSORT
      RECURSIVE SUBROUTINE IQSORT5 (I,J,N,A,B,C,D,E)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: I,J,N
      INTEGER(8),INTENT(INOUT),DIMENSION (:) :: A,B,C,D,E
      INTEGER(4) :: IPV,KPV
      INTEGER(8) :: PV
      INTEGER(4) :: S0,S1,S2,K,H
      INTEGER(4) :: LET,RIT
      INTEGER(8) :: WKA,WKB,WKC,WKD,WKE
      S0 = J-I
      S1 = (S0+1)/2
      IF (MOD(S0,2).EQ.0) THEN
        S2 = S1
      ELSE
        S2 = S1-1
      END IF
      H = 1
      K = I+S1
      DO
        IF (A(I).NE.A(K).OR.H.GE.S0) THEN
          EXIT
        ELSE
          IF (MOD(H,2).EQ.1) THEN
            K = K+S2
          ELSE
            K = K+S1
          END IF
          IF (K.GT.J)     K = K-J+I-1
          H = H+1
        END IF
      END DO
      IF (A(I).GT.A(K)) THEN
        IPV = I
      ELSE  IF (A(I).LT.A(K)) THEN
        IPV = K
      ELSE
        IPV = 0
      END IF
      IF (IPV.EQ.0) RETURN
      PV = A(IPV)
      LET = I
      RIT = J
      DO
        IF (A(LET).GE.PV) THEN
          EXIT
        ELSE
          LET = LET+1
        END IF
      END DO
      DO
        IF (A(RIT).LT.PV) THEN
          EXIT
        ELSE
          RIT = RIT-1
        END IF
      END DO
      DO
        IF (LET.GT.RIT) THEN
          EXIT
        ELSE
          WKA = A(LET) ; A(LET) = A(RIT) ; A(RIT) = WKA
          WKB = B(LET) ; B(LET) = B(RIT) ; B(RIT) = WKB
          WKC = C(LET) ; C(LET) = C(RIT) ; C(RIT) = WKC
          WKD = D(LET) ; D(LET) = D(RIT) ; D(RIT) = WKD
          WKE = E(LET) ; E(LET) = E(RIT) ; E(RIT) = WKE
          LET = LET+1     ; RIT = RIT-1
          DO
            IF (A(LET).GE.PV) THEN
              EXIT
            ELSE
              LET = LET+1
            END IF
          END DO
          DO
            IF (A(RIT).LT.PV) THEN
              EXIT
            ELSE
              RIT = RIT-1
            END IF
          END DO
        END IF
      END DO
      KPV = LET
      CALL IQSORT5 (I,KPV-1,N,A,B,C,D,E)
      CALL IQSORT5 (KPV,J,N,A,B,C,D,E)
      RETURN
      END SUBROUTINE IQSORT5
      RECURSIVE SUBROUTINE RQSORT (I,J,N,A,B,C,D)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: I,J,N
      INTEGER(8),INTENT(INOUT),DIMENSION (:) :: A,B,C
      REAL(4),INTENT(INOUT),DIMENSION (:)    :: D
      INTEGER(4) :: IPV,KPV
      INTEGER(8) :: PV
      INTEGER(4) :: S0,S1,S2,K,H
      INTEGER(4) :: LET,RIT
      INTEGER(8) :: WKA,WKB,WKC
      REAL(4)    :: WKD
      S0 = J-I
      S1 = (S0+1)/2
      IF (MOD(S0,2).EQ.0) THEN
        S2 = S1
      ELSE
        S2 = S1-1
      END IF
      H = 1
      K = I+S1
      DO
        IF (A(I).NE.A(K).OR.H.GE.S0) THEN
          EXIT
        ELSE
          IF (MOD(H,2).EQ.1) THEN
            K = K+S2
          ELSE
            K = K+S1
          END IF
          IF (K.GT.J)     K = K-J+I-1
          H = H+1
        END IF
      END DO
      IF (A(I).GT.A(K)) THEN
        IPV = I
      ELSE  IF (A(I).LT.A(K)) THEN
        IPV = K
      ELSE
        IPV = 0
      END IF
      IF (IPV.EQ.0) RETURN
      PV = A(IPV)
      LET = I
      RIT = J
      DO
        IF (A(LET).GE.PV) THEN
          EXIT
        ELSE
          LET = LET+1
        END IF
      END DO
      DO
        IF (A(RIT).LT.PV) THEN
          EXIT
        ELSE
          RIT = RIT-1
        END IF
      END DO
      DO
        IF (LET.GT.RIT) THEN
          EXIT
        ELSE
          IF (N.EQ.1) THEN
            WKA = A(LET) ; A(LET) = A(RIT) ; A(RIT) = WKA
            WKD = D(LET) ; D(LET) = D(RIT) ; D(RIT) = WKD
          ELSE
            WKA = A(LET) ; A(LET) = A(RIT) ; A(RIT) = WKA
            WKB = B(LET) ; B(LET) = B(RIT) ; B(RIT) = WKB
            WKC = C(LET) ; C(LET) = C(RIT) ; C(RIT) = WKC
            WKD = D(LET) ; D(LET) = D(RIT) ; D(RIT) = WKD
          END IF
          LET = LET+1     ; RIT = RIT-1
          DO
            IF (A(LET).GE.PV) THEN
              EXIT
            ELSE
              LET = LET+1
            END IF
          END DO
          DO
            IF (A(RIT).LT.PV) THEN
              EXIT
            ELSE
              RIT = RIT-1
            END IF
          END DO
        END IF
      END DO
      KPV = LET
      CALL RQSORT (I,KPV-1,N,A,B,C,D)
      CALL RQSORT (KPV,J,N,A,B,C,D)
      RETURN
      END SUBROUTINE RQSORT
      INTEGER(4) FUNCTION IFIND (ID,N,IN)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN) :: N
      INTEGER(8),INTENT(IN) :: ID
      INTEGER(8),INTENT(IN) :: IN(:)
      INTEGER(4)            :: I,J,K
      I=1
      J=N
      DO
        IF(I.GT.J) THEN
          IFIND = 0
          EXIT
        ELSE
          K = (I+J)/2
          IF(ID.EQ.IN(K)) THEN
            IFIND = K
            EXIT
          ELSE IF(ID.GT.IN(K)) THEN
            I=K+1
          ELSE
            J=K-1
          ENDIF
        ENDIF
      END DO
      RETURN
      END FUNCTION IFIND
      SUBROUTINE INBCOE (NA,NS,P1,P2,D,F)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                    :: NA,NS
      INTEGER(8),INTENT(INOUT),DIMENSION (:)   :: P1,P2
      REAL(4),INTENT(INOUT),DIMENSION (:)      :: D
      REAL(4),INTENT(INOUT),DIMENSION (0:)     :: F
      INTEGER(8)                    :: I,J,K
      INTEGER(8)                    :: IS,ID,KS,KD,IWK
      INTEGER(8),DIMENSION (NA)     :: POINT
      REAL(4),DIMENSION (NA)        :: LM
      REAL(4)                       :: FI,R
      INTEGER(4)                    :: PRINT,WPRINT
      PRINT = 0
      IF (NA.LE.200000) THEN
        WPRINT = 100000
      ELSE
        WPRINT = NA*0.8
      END IF
      F(0) = -1.0
      LM = 0
      DO I=1,NS
        D(I) = 0.5-0.25*(F(P1(I))+F(P2(I)))
      END DO
      DO I=1,NA
        IF (P1(I).LT.P2(I)) THEN
          IWK = P1(I)
          P1(I) = P2(I)
          P2(I) = IWK
        END IF
      END DO
      DO I=NS+1,NA
        IS = P1(I)
        ID = P2(I)
        D(I) = 0.5-0.25*(F(IS)+F(ID))
        POINT(I) = 0
        IF (IS.EQ.0.OR.ID.EQ.0) THEN
          F(I) = 0.0
        ELSE IF (P1(I-1).EQ.P1(I).AND.P2(I-1).EQ.P2(I)) THEN
          F(I) = F(I-1)
        ELSE
          FI = -1.0
          LM(I) = 1.0
          J = I
          DO
            IF (J.EQ.0)     EXIT
            K = J
            R = 0.5*LM(K)
            KS = P1(K)
            KD = P2(K)
            IF (KS.GT.0) THEN
              DO
                IF (POINT(K).LE.KS) THEN
                  EXIT
                ELSE
                  K = POINT(K)
                END IF
              END DO
              LM(KS) = LM(KS)+R
              IF (KS.NE.POINT(K)) THEN
                POINT (KS) = POINT(K)
                POINT(K) = KS
              END IF
              IF (KD.GT.0) THEN
                DO
                  IF (POINT(K).LE.KD) THEN
                    EXIT
                  ELSE
                    K = POINT(K)
                  END IF
                END DO
                LM(KD) = LM(KD)+R
                IF (KD.NE.POINT(K)) THEN
                  POINT(KD) = POINT(K)
                  POINT(K) = KD
                END IF
              END IF
            END IF
            FI = FI+LM(J)*LM(J)*D(J)
            LM(J) = 0.0
            K = J
            J = POINT(J)
            POINT(K) = 0
          END DO
          F(I) = FI
        END IF
        IF (MOD(I,10000).EQ.0) THEN
          IF (I.GE.WPRINT) THEN
            PRINT = PRINT+1
            IF (PRINT.EQ.1) THEN
              WRITE (6,'(A,I10)') '近交係数の算出が終了した個体数 =',I
            ELSE IF (I.GT.100000) THEN
              WRITE (6,'(A,I10)') '                               =',I
            END IF
          END IF
        END IF
      END DO
      RETURN
      END SUBROUTINE INBCOE
      SUBROUTINE RLSCOE (NS,P1,P2,D,F)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                    :: NS
      INTEGER(8),INTENT(INOUT),DIMENSION (:)   :: P1,P2
      REAL(4),INTENT(INOUT),DIMENSION (:)      :: D
      REAL(4),INTENT(INOUT),DIMENSION (0:)     :: F
      INTEGER(8)                    :: J,K
      INTEGER(8)                    :: IS,ID,KS,KD,IWK
      INTEGER(8),DIMENSION (NS+1)   :: POINT
      REAL(4),DIMENSION (NS+1)      :: LM
      REAL(4)                       :: FI,R
      F(0) = -1.0
      LM = 0
      POINT(1:NS) = 0
      IF (P1(NS+1).LT.P2(NS+1)) THEN
        IWK = P1(NS+1)
        P1(NS+1) = P2(NS+1)
        P2(NS+1) = IWK
      END IF
      IS = P1(NS+1)
      ID = P2(NS+1)
      D(NS+1) = 0.5-0.25*(F(IS)+F(ID))
      POINT(NS+1) = 0
      FI = -1.0
      LM(NS+1) = 1.0
      J = NS+1
      DO
        IF (J.EQ.0)     EXIT
        K = J
        R = 0.5*LM(K)
        KS = P1(K)
        KD = P2(K)
        IF (KS.GT.0) THEN
          DO
            IF (POINT(K).LE.KS) THEN
              EXIT
            ELSE
              K = POINT(K)
            END IF
          END DO
          LM(KS) = LM(KS)+R
          IF (KS.NE.POINT(K)) THEN
            POINT (KS) = POINT(K)
            POINT(K) = KS
          END IF
          IF (KD.GT.0) THEN
            DO
              IF (POINT(K).LE.KD) THEN
                EXIT
              ELSE
                K = POINT(K)
              END IF
            END DO
            LM(KD) = LM(KD)+R
            IF (KD.NE.POINT(K)) THEN
              POINT(KD) = POINT(K)
              POINT(K) = KD
            END IF
          END IF
        END IF
        FI = FI+LM(J)*LM(J)*D(J)
        LM(J) = 0.0
        K = J
        J = POINT(J)
        POINT(K) = 0
      END DO
      F(NS+1) = FI
      RETURN
      END SUBROUTINE RLSCOE
