      PROGRAM SINDEX
      IMPLICIT NONE
      INTEGER,PARAMETER                       :: L1=20,L2=200,L0=L1*L2
      CHARACTER (LEN=80)                      :: TIT
      CHARACTER (LEN=20)                      :: TNAME(L1),GNAME(L2)
      INTEGER                                 :: NT0
      REAL(8),DIMENSION (L2,L2)               :: RB
      REAL(8),DIMENSION (L1)                  :: A,D
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      REAL(8),DIMENSION (L1)                  :: VP,H2,T2
      REAL(8),DIMENSION (L1,L1)               :: RG,RP
      INTEGER(4),DIMENSION (L2*2)             :: P1,P2
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      INTEGER(4)                              :: NP0
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      REAL(8)                                 :: SI,SH,SIH,RIH
      REAL(8),DIMENSION (L1)                  :: DG
      INTEGER(4)                              :: NMAT,PNO1,PNO2,PNO3
      INTEGER(4),DIMENSION (4)                :: USEP
      INTEGER(4),DIMENSION (4)                :: PDG0
      INTEGER(4),DIMENSION (4,L1)             :: PDG1,PDG2
      INTEGER(4)                              :: RPT,RPT0
      INTERFACE
        SUBROUTINE PBACK (RPT)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)               :: RPT
        END SUBROUTINE PBACK
        SUBROUTINE PROUT (RPT,*)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)               :: RPT
        END SUBROUTINE PROUT
      END INTERFACE
      COMMON /ETC01/NT0
      COMMON /ETC02/RB
      COMMON /ETC03/A,D
      COMMON /TITLE/TIT,TNAME,GNAME
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /TRAIT/VP,H2,RG,RP,T2
      COMMON /PEDIG/P1,P2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /PRINT/NP0
      COMMON /MATRX/P,G,F,W,B,C0
      COMMON /RSLTS/SI,SH,SIH,RIH,DG
      COMMON /BACKS/NMAT,PNO1,PNO2,PNO3,USEP,PDG0,PDG1,PDG2
      RPT = 0
      CALL CNTRL (*999)
      IF (NG0.NE.0) THEN
        CALL AMTRX
        CALL TMTRX (*999)
        CALL CMTRX (*999)
        CALL RIHDG (*999)
      END IF
      CALL PROUT (RPT,*999)
      IF (NMAT.NE.0) THEN
        IF (NMAT.LE.(L2-27)/8) THEN
          RPT0 = 4
        ELSE IF (NMAT.LE.(L2-13)/2) THEN
          RPT0 = 3
        ELSE
          RPT0 = 2
        END IF
        DO RPT=1,RPT0
          CALL PBACK (RPT)
          CALL AMTRX
          CALL TMTRX (*999)
          CALL CMTRX (*999)
          CALL RIHDG (*999)
          CALL PROUT (RPT,*999)
        END DO
      END IF
      WRITE (6,*) 'プログラムは正常終了しました.'
  999 STOP
      END PROGRAM SINDEX
      SUBROUTINE PBACK (RPT)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                   :: RPT
      INTEGER(4),PARAMETER                    :: L1=20,L2=200
      INTEGER(4)                              :: NT0
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      INTEGER(4),DIMENSION (L2*2)             :: P1,P2
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NMAT,PNO1,PNO2,PNO3
      INTEGER(4),DIMENSION (4)                :: USEP
      INTEGER(4),DIMENSION (4)                :: PDG0
      INTEGER(4),DIMENSION (4,L1)             :: PDG1,PDG2
      INTEGER(4)                              :: I,J,K
      INTEGER(4)                              :: IWKT(L1),IWKR(L1)
      INTEGER(4),DIMENSION(4)                 :: IWKP
      INTEGER(4),DIMENSION (L2)               :: IKND
      SAVE /PEDIG/,/GROUP/,/MODEL/
      COMMON /ETC01/NT0
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /PEDIG/P1,P2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /BACKS/NMAT,PNO1,PNO2,PNO3,USEP,PDG0,PDG1,PDG2
      P1 = 0; P2 = 0
      NA = 0
      NT = 0
      NR = 0
      IF (RPT.EQ.1) THEN
        NG0 = 4
        CAN = 3
        IWKP = 0
        IWKP(3) = 1
        IWKT = 0
        IWKR = 0
        DO J=1,PDG0(3)
          IWKT(PDG1(3,J)) = 3
          IWKR(PDG1(3,J)) = PDG2(3,J)
        END DO
        IF (USEP(4).EQ.1) THEN
          DO J=1,PDG0(4)
            IF (IWKT(PDG1(4,J)).EQ.0) THEN
              IWKP(4) = 1
              IWKT(PDG1(4,J)) = 4
              IWKR(PDG1(4,J)) = PDG2(4,J)
            END IF
          END DO
        END IF
        IF (USEP(1).EQ.1) THEN
          DO I=1,2
            DO J=1,PDG0(I)
              IF (IWKT(PDG1(I,J)).EQ.0) THEN
                IWKP(I) = 1
                IWKT(PDG1(I,J)) = I
                IWKR(PDG1(I,J)) = PDG2(I,J)
              END IF
            END DO
          END DO
        END IF
        P1(1) = 0; P2(1) = 0
        P1(2) = 0; P2(2) = 0
        P1(3) = 1; P2(3) = 2
        P1(4) = 1; P2(4) = 2
        IF (IWKP(1).EQ.1.OR.IWKP(2).EQ.1) THEN
          NA(1) = 1; NA(2) = 1; NA(3) = 1
        ELSE
          NA(1) = 0; NA(2) = 0; NA(3) = 1
        END IF
        IF (IWKP(4).EQ.1) THEN
          NA(4) = PNO3
        ELSE
          NA(4) = 0
        END IF
        DO I=1,L1
          IF (IWKT(I).NE.0) THEN
            NT(IWKT(I),0) = NT(IWKT(I),0)+1
            NR(IWKT(I),0) = NR(IWKT(I),0)+1
            NT(IWKT(I),NT(IWKT(I),0)) = I
            NR(IWKT(I),NR(IWKT(I),0)) = IWKR(I)
          END IF
        END DO
      ELSE IF (RPT.EQ.2) THEN
        NG0 = 5
        CAN = 3
        P1(1) = 0; P2(1) = 0; NA(1) = 1
        P1(2) = 0; P2(2) = 0; NA(2) = 1
        P1(3) = 1; P2(3) = 2; NA(3) = 1
        P1(4) = 1; P2(4) = 2; NA(4) = PNO3
        P1(5) = 1; P2(5) = 2; NA(5) = PNO1+PNO2
        DO I=1,4
          NT(I,0) = PDG0(I); NR(I,0) = PDG0(I)
          DO J=1,NT(I,0)
            NT(I,J) = PDG1(I,J)
            NR(I,J) = PDG2(I,J)
          END DO
        END DO
        NT(5,0) = PDG0(3); NR(5,0) = PDG0(3)
        DO J=1,NT(5,0)
          NT(5,J) = PDG1(3,J)
          NR(5,J) = PDG2(3,J)
        END DO
      ELSE IF (RPT.EQ.3) THEN
        NG0 = 13+2*NMAT
        CAN = 13
        DO I=1,4
          P1(I) = 0; P2(I) = 0
        END DO
        DO I=5,8
          P1(I) = 1; P2(I) = 2
        END DO
        DO I=9,12
          P1(I) = 3; P2(I) = 4
        END DO
        DO I=13,15
          P1(I) = 5; P2(I) = 9
        END DO
        DO I=1,NMAT-1
          P1(15+I) = 5       ; P2(15+I) = 15+(NMAT-1)*2+I
          P1(15+NMAT-1+I) = 5; P2(15+NMAT-1+I) = 15+(NMAT-1)*2+I
        END DO
        DO I=1,17+2*(NMAT-1)
          IF (I.LE.5.OR.I.EQ.9.OR.I.EQ.13) THEN
            NA(I) = 1
          ELSE IF (I.EQ.6) THEN
            NA(I) = PNO1-1
          ELSE IF (I.EQ.7) THEN
            NA(I) = PNO2
          ELSE IF (I.EQ.10) THEN
            NA(I) = PNO1
          ELSE IF (I.EQ.11) THEN
            NA(I) = PNO2-1
          ELSE IF (I.EQ.14) THEN
            NA(I) = PNO1+PNO2-1
          ELSE IF (I.GE.16.AND.I.LE.15+(NMAT-1)) THEN
            NA(I) = PNO1+PNO2
          ELSE
            NA(I) = PNO3
          END IF
        END DO
        IKND(1) = 1;  IKND(2) = 2 ; IKND(3) = 1  ; IKND(4) = 2
        IKND(5) = 1;  IKND(6) = 1 ; IKND(7) = 2  ; IKND(8) = 4
        IKND(9) = 2;  IKND(10) = 1; IKND(11) = 2 ; IKND(12) = 4
        IKND(13) = 3; IKND(14) = 3; IKND(15) = 4
        IKND(16:15+(NMAT-1)) = 3
        IKND(15+(NMAT-1)+1:15+2*(NMAT-1)) = 4
        DO I=1,15+2*(NMAT-1)
          NT(I,0) = PDG0(IKND(I))
          NR(I,0) = PDG0(IKND(I))
          DO J=1,NT(I,0)
            NT(I,J) = PDG1(IKND(I),J)
            NR(I,J) = PDG2(IKND(I),J)
          END DO
        END DO
      ELSE IF (RPT.EQ.4) THEN
        NG0 = 27+8*NMAT
        CAN = 33
        DO I=1,8
          P1(I) = 0; P2(I) = 0
        END DO
        DO I=1,4
          P1(8+I) =   1; P2(8+I)  = 2
          P1(12+I) =  3; P2(12+I) = 4
          P1(16+I) =  5; P2(16+I) = 6
          P1(20+I) =  7; P2(20+I) = 8
          P1(24+I) =  9; P2(24+I) = 13
          P1(28+I) = 17; P2(28+I) = 21
        END DO
        DO I=1,3
          P1(32+I) = 25; P2(32+I) = 29
        END DO
        DO I=1,NMAT-1
          P1(35+(NMAT-1)*0+I) = 9 ; P2(35+(NMAT-1)*0+I) = 27+NMAT*8+I
          P1(35+(NMAT-1)*1+I) = 9 ; P2(35+(NMAT-1)*1+I) = 27+NMAT*8+I
          P1(35+(NMAT-1)*2+I) = 9 ; P2(35+(NMAT-1)*2+I) = 27+NMAT*8+I
          P1(35+(NMAT-1)*3+I) = 17; P2(35+(NMAT-1)*3+I) = 26+NMAT*9+I
          P1(35+(NMAT-1)*4+I) = 17; P2(35+(NMAT-1)*4+I) = 26+NMAT*9+I
          P1(35+(NMAT-1)*5+I) = 17; P2(35+(NMAT-1)*5+I) = 26+NMAT*9+I
          P1(35+(NMAT-1)*6+I) = 25; P2(35+(NMAT-1)*6+I) = 25+NMAT*10+I
          P1(35+(NMAT-1)*7+I) = 25; P2(35+(NMAT-1)*7+I) = 25+NMAT*10+I
        END DO
        DO I=1,NG0
          NA(I) = 1
        END DO
        NA(10) = PNO1-1; NA(11) = PNO2  ; NA(12) = PNO3
        NA(14) = PNO1  ; NA(15) = PNO2-1; NA(16) = PNO3
        NA(18) = PNO1-1; NA(19) = PNO2  ; NA(20) = PNO3
        NA(22) = PNO1  ; NA(23) = PNO2-1; NA(24) = PNO3
        NA(26) = PNO1-1; NA(27) = PNO2  ; NA(28) = PNO3
        NA(30) = PNO1  ; NA(31) = PNO2-1; NA(32) = PNO3
        NA(34) = PNO1+PNO2-1            ; NA(35) = PNO3
        DO I=1,NMAT-1
          NA(35+(NMAT-1)*0+I) = PNO1
          NA(35+(NMAT-1)*1+I) = PNO2
          NA(35+(NMAT-1)*2+I) = PNO3
          NA(35+(NMAT-1)*3+I) = PNO1
          NA(35+(NMAT-1)*4+I) = PNO2
          NA(35+(NMAT-1)*5+I) = PNO3
          NA(35+(NMAT-1)*6+I) = PNO1+PNO2
          NA(35+(NMAT-1)*7+I) = PNO3
        END DO
        DO I=1,8
          IF (MOD(I,2).EQ.1) THEN
            IKND(I) = 1
          ELSE
            IKND(I) = 2
          END IF
        END DO
        DO I=9,32
          IF (MOD(I,4).EQ.2) THEN
            IKND(I) = 1
          ELSE IF (MOD(I,4).EQ.3) THEN
            IKND(I) = 2
          ELSE IF (MOD(I,4).EQ.0) THEN
            IKND(I) = 4
          ELSE IF (I.EQ.9.OR.I.EQ.17.OR.I.EQ.25) THEN
            IKND(I) = 1
          ELSE
            IKND(I) = 2
          END IF
        END DO
        IKND(33) = 3; IKND(34) = 3; IKND(35) = 4
        DO I=1,NMAT-1
          IKND(35+(NMAT-1)*0+I) = 1
          IKND(35+(NMAT-1)*1+I) = 2
          IKND(35+(NMAT-1)*2+I) = 4
          IKND(35+(NMAT-1)*3+I) = 1
          IKND(35+(NMAT-1)*4+I) = 2
          IKND(35+(NMAT-1)*5+I) = 4
          IKND(35+(NMAT-1)*6+I) = 3
          IKND(35+(NMAT-1)*7+I) = 4
        END DO
        DO I=1,27+8*NMAT
          NT(I,0) = PDG0(IKND(I))
          NR(I,0) = PDG0(IKND(I))
          DO J=1,NT(I,0)
            NT(I,J) = PDG1(IKND(I),J)
            NR(I,J) = PDG2(IKND(I),J)
          END DO
        END DO
      END IF
      LMAX = 0
      DO I=1,NT0
        DO J=1,NG0
          DO K=1,NT(J,0)
            IF (NT(J,K).EQ.I) THEN
              LMAX = LMAX+1
            END IF
          END DO
        END DO
      END DO
      RETURN
      END SUBROUTINE PBACK
      SUBROUTINE CNTRL (*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200
      CHARACTER (LEN=80)                      :: TIT
      CHARACTER (LEN=20)                      :: TNAME(L1),GNAME(L2)
      INTEGER(4)                              :: NT0
      REAL(8),DIMENSION (L1)                  :: A,D
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      REAL(8),DIMENSION (L1)                  :: VP,H2,T2
      REAL(8),DIMENSION (L1,L1)               :: RG,RP
      INTEGER(4),DIMENSION (L2*2)             :: P1,P2
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      INTEGER(4)                              :: NP0
      INTEGER(4)                              :: NMAT,PNO1,PNO2,PNO3
      INTEGER(4),DIMENSION (4)                :: USEP
      INTEGER(4),DIMENSION (4)                :: PDG0
      INTEGER(4),DIMENSION (4,L1)             :: PDG1,PDG2
      INTEGER(4)                              :: I,J,K,N
      CHARACTER (LEN=80)                      :: CNFILE
      CHARACTER (LEN=6)                       :: CNTWK
      CHARACTER (LEN=6),DIMENSION (0:6)       :: CNT1,CNT2
      INTEGER(4)                              :: NTWK,IWK0,IWK1
      INTEGER(4),DIMENSION (0:L1)             :: IWK
      REAL(8)                                 :: DWK0
      REAL(8),DIMENSION (L1,L1)               :: DWKP,DWKG
      INTEGER(4),DIMENSION (0:6)              :: SROW,EROW,NROW
      INTEGER(4),PARAMETER                    :: LCL=999
      CHARACTER (LEN=999)                     :: INFILE
      CHARACTER (LEN=6)                       :: CWK0
      CHARACTER (LEN=1),DIMENSION (999)       :: WX
      INTEGER(4)                              :: IN
      INTEGER(4),DIMENSION (0:6)              :: IWKC
      LOGICAL                                 :: TF
      CHARACTER (LEN=6),DIMENSION(0:5)        :: PNAME1,PNAME2
      SAVE /ETC01/,/ETC03/,/TITLE/,/MODEL/,/TRAIT/,/PEDIG/,/GROUP/,
     &     /WEIGH/,/PRINT/,/BACKS/
      COMMON /ETC01/NT0
      COMMON /ETC03/A,D
      COMMON /TITLE/TIT,TNAME,GNAME
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /TRAIT/VP,H2,RG,RP,T2
      COMMON /PEDIG/P1,P2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /PRINT/NP0
      COMMON /BACKS/NMAT,PNO1,PNO2,PNO3,USEP,PDG0,PDG1,PDG2
      DATA CNT1/'TITLE ','MODEL ','TRAIT ','GROUP ','PEDIG ','WEIGHT',
     &          'PRINT '/
      DATA CNT2/'title ','model ','trait ','group ','pedig ','weight',
     &          'print '/
      DATA (IWK(I),I=0,L1)/0,L1*0/
      DATA PNAME1/'      ','ANC(A)','ANC(M)','ANC(F)','CAND  ','SIBT  '/
      DATA PNAME2/'      ','anc(a)','anc(m)','anc(f)','cand  ','sibt  '/
      IN = 11
      SROW = 0
      EROW = 0
      NROW = 0
      WRITE (6,*) 'コントロールファイル名を入力して下さい ----->'
      READ (5,'(A80)') CNFILE
      INQUIRE (FILE=CNFILE,EXIST=TF)
      IF (.NOT.TF)     CALL ERROR (11,1,*999)
      OPEN (11,FILE=CNFILE)
      N = 0
      IWK1 = -1
      IWKC = 0
      DO
        READ (IN,'(A)',END=10) INFILE
        READ (INFILE,'(999A1)') WX(1:LCL)
        IWK0 = 0
        DO I=1,LCL
          IF(WX(I).NE.' ') THEN
            IWK0 = I
            EXIT
          END IF
        END DO
        READ (INFILE(IWK0:),'(A)') CWK0
        N = N+1
        DO I=0,6
          IF (CWK0.EQ.CNT1(I).OR.CWK0.EQ.CNT2(I)) THEN
            SROW(I) = N
            IWKC(I) = IWKC(I)+1
            IF (IWK1.NE.-1) THEN
              EROW(IWK1) = N-1
            END IF
            IWK1 = I
          END IF
        END DO
      END DO
   10 EROW(IWK1) = N
      DO I=0,6
        IF (IWKC(1).GE.2)     CALL ERROR (1,1,*999)
      END DO
      DO I=1,5
        IF (I.EQ.3.OR.I.EQ.4) THEN
          IF (IWKC(3).NE.1.AND.IWKC(4).NE.1) CALL ERROR (1,1,*999)
        ELSE
          IF (IWKC(I).NE.1)     CALL ERROR (1,1,*999)
        END IF
      END DO
      DO I=0,6
        IF (SROW(I).NE.0) THEN
          NROW(I) = EROW(I)-SROW(I)
        END IF
      END DO
      NMAT = 0
      IF (IWKC(0).EQ.0.AND.NROW(0).LE.0) THEN
        TIT = ' < この分析にはタイトルがありません. > '
      ELSE IF (NROW(0).GE.1) THEN
        REWIND (IN)
        DO I=1,SROW(0)
          READ (IN,*)
        END DO
        READ (IN,'(A)') TIT
      ELSE
        REWIND (IN)
        DO I=1,SROW(0)-1
          READ (IN,*)
        END DO
        READ (IN,'(A)') INFILE
        IWK0 = INDEX(INFILE,'T')
        IF (IWK0.EQ.0) THEN
          IWK0 = INDEX(INFILE,'t')
        END IF
        IF (IWK0.EQ.0)     CALL ERROR (1,2,*999)
        READ (INFILE,'(999A1)') WX(1:LCL)
        DO I=IWK0+5,LCL
          IF(WX(I).NE.' ') THEN
            IWK0 = I
            EXIT
          END IF
        END DO
        READ (INFILE(IWK0:),'(A)') TIT
      END IF
      REWIND (IN)
      DO I=1,SROW(1)-1
        READ (IN,*)
      END DO
      READ (IN,*) CNTWK,NTWK
      CNTWK=CNTWK; NM0 = NTWK
      IF (NM0.LE.0.OR.NM0.GT.4)     CALL ERROR (1,3,*999)
      IF (NM0.EQ.4) THEN
        READ (IN,*) ERATE
      END IF
      REWIND (IN)
      DO I=1,SROW(2)-1
        READ (IN,*)
      END DO
      READ (IN,'(A)') INFILE
      NTWK = INDEX(INFILE,'*')
      IF (NTWK.EQ.0) THEN
        READ (INFILE,'(999A1)') WX(1:LCL)
        IWK0 = INDEX(INFILE,'T')
        DO I=IWK0+5,LCL
          IF(WX(I).EQ.'1'.AND.WX(I+1).NE.' ') THEN
            NTWK = I
            EXIT
          END IF
        END DO
      END IF
      NT0 = NROW(2)
      IF (NT0.GT.L1)     CALL ERROR (2,4,*999)
      DO I=1,NT0
        READ (IN,*) IWK0,TNAME(I),VP(I),DWKP(I,1:I-1),DWKG(I,I:NT0)
        IF (NTWK.EQ.0) THEN
          H2(I) = DWKG(I,I)
          RP(I,I) = 1.0
          RG(I,I) = 1.0
          DO J=1,I-1
            RP(I,J) = DWKP(I,J)
          END DO
          DO J=I+1,NT0
            RG(I,J) = DWKG(I,J)
          END DO
        END IF
      END DO
      IF (NTWK.NE.0) THEN
        DO I=1,NT0
          H2(I) = DWKG(I,I)/VP(I)
          RP(I,I) = 1.0
          RG(I,I) = 1.0
          DO J=1,I-1
            RP(I,J) = DWKP(I,J)/DSQRT(VP(I)*VP(J))
          END DO
          DO J=I+1,NT0
            RG(I,J) = DWKG(I,J)/DSQRT(DWKG(I,I)*DWKG(J,J))
          END DO
        END DO
      END IF
      DO I=1,NT0
        DO J=1,I
          RP(J,I) = RP(I,J)
          RG(I,J) = RG(J,I)
        END DO
      END DO
      IF (IWKC(3).EQ.1) THEN
        REWIND (IN)
        DO I=1,SROW(3)
          READ (IN,*)
        END DO
        NG0 = NROW(3)
        IF (NG0.GT.L2)     CALL ERROR (2,5,*999)
        DO I=1,NG0
          READ (IN,'(A)') INFILE
          READ (INFILE,'(999A1)') WX(1:LCL)
          IWK0 = 0
          DO J=1,LCL-1
            IF (WX(J).NE.' '.AND.WX(J+1).EQ.' ') THEN
              IWK0 = IWK0+1
            END IF
          END DO
          IF (IWK0.EQ.5) THEN
            READ (INFILE,*) IWK1,GNAME(I),P1(I),P2(I),NA(I)
            NT(I,0) = 0
          ELSE
            READ (INFILE,*) IWK1,GNAME(I),P1(I),P2(I),NA(I),NT(I,0),
     &                       NT(I,1:NT(I,0)),NR(I,1:NT(I,0))
          END IF
          NR(I,0) = NT(I,0)
          IF (INDEX(GNAME(I),'*').NE.0) THEN
            CAN = I
            GNAME(I) = GNAME(I)(INDEX(GNAME(I),'*'):)
          END IF
        END DO
        DO I=1,NG0
          DO J=1,NR(I,0)
            IF (NR(I,J).GT.1) THEN
              IWK(0) = 1
              IWK(NT(I,J)) = 1
            END IF
          END DO
        END DO
      ELSE
        NG0 = 0
      END IF
      IF (IWKC(4).EQ.1) THEN
        REWIND (IN)
        PDG0 = 0
        PDG1 = 0
        PDG2 = 0
        USEP = 0
        IF (NROW(4)-2.LT.1.OR.NROW(4)-2.GT.4)     CALL ERROR (1,6,*999)
        DO I=1,SROW(4)
          READ (IN,*)
        END DO
        READ (IN,*) NMAT
        IF ((NMAT*2+13).GT.L2)     CALL ERROR (2,6,*999)
        READ (IN,'(A)') INFILE
        READ (INFILE,'(999A1)') WX(1:LCL)
        IWK0 = 0
        DO J=1,LCL-1
          IF (WX(J).NE.' '.AND.WX(J+1).EQ.' ') THEN
            IWK0 = IWK0+1
          END IF
        END DO
        IF (IWK0.EQ.2) THEN
          READ (INFILE,*) PNO1,PNO2
          PNO3 = 0
        ELSE IF (IWK0.EQ.3) THEN
          READ (INFILE,*) PNO1,PNO2,PNO3
        ELSE
          CALL ERROR (1,6,*999)
        END IF
        IF (PDG0(IWK0).GT.L1)     CALL ERROR (2,6,*999)
        DO I=1,NROW(4)-2
          IWK0 = 0
          READ (IN,'(A)') INFILE
          READ (INFILE,*) PNAME1(0)
          DO J=1,5
            IF (PNAME1(0).EQ.PNAME1(J).OR.PNAME1(0).EQ.PNAME2(J)) THEN
              IWK0 = J
              EXIT
            END IF
          END DO
          IF (IWK0.EQ.0.OR.IWK0.GT.5)     CALL ERROR (1,6,*999)
          IF (IWK0.EQ.1) THEN
            IF (USEP(2).EQ.1)     CALL ERROR (1,6,*999)
            USEP(1) = 1; USEP(2) = 1
            READ (INFILE,*) PNAME1(IWK0),PDG0(IWK0),
     &                   PDG1(IWK0,1:PDG0(IWK0)),PDG2(IWK0,1:PDG0(IWK0))
            IF (PDG0(IWK0).GT.L1)     CALL ERROR (2,6,*999)
            PDG0(2) = PDG0(IWK0)
            PDG1(2,1:PDG0(2)) = PDG1(IWK0,1:PDG0(IWK0))
            PDG2(2,1:PDG0(2)) = PDG2(IWK0,1:PDG0(IWK0))
          ELSE IF (IWK0.GE.2) THEN
            IWK0 = IWK0-1
            IF (USEP(IWK0).EQ.1)     CALL ERROR (1,6,*999)
            USEP(IWK0) = 1
            READ (INFILE,*) PNAME1(IWK0),PDG0(IWK0),
     &                   PDG1(IWK0,1:PDG0(IWK0)),PDG2(IWK0,1:PDG0(IWK0))
            IF (PDG0(IWK0).GT.L1)     CALL ERROR (2,6,*999)
          END IF
        END DO
        IF (USEP(3).EQ.0)     CALL ERROR (1,6,*999)
        IF ((USEP(1).EQ.1.AND.USEP(2).NE.1).OR.
     &      (USEP(1).NE.1.AND.USEP(2).EQ.1))     CALL ERROR (1,6,*999)
      END IF
      IF (USEP(1).EQ.0) THEN
        PDG0(1) = PDG0(3)
        PDG0(2) = PDG0(3)
        PDG1(1,1:PDG0(1)) = PDG1(3,1:PDG0(3))
        PDG2(1,1:PDG0(1)) = PDG2(3,1:PDG0(3))
        PDG1(2,1:PDG0(2)) = PDG1(3,1:PDG0(3))
        PDG2(2,1:PDG0(2)) = PDG2(3,1:PDG0(3))
      END IF
      IF (USEP(4).EQ.0) THEN
        PDG0(4) = PDG0(3)
        PDG1(4,1:PDG0(4)) = PDG1(3,1:PDG0(3))
        PDG2(4,1:PDG0(4)) = PDG2(3,1:PDG0(3))
      END IF
      REWIND (IN)
      DO I=1,SROW(5)
        READ (IN,*)
      END DO
      NW0 = NROW(5)
      IF (NW0.GT.L1)     CALL ERROR (2,7,*999)
      WND = 0
      DO I=1,NW0
        IF (NM0.EQ.1.OR.NM0.EQ.2) THEN
          READ (IN,'(A)') INFILE
          READ (INFILE,'(999A1)') WX(1:LCL)
          IWK0 = 0
          DO J=1,LCL-1
            IF (WX(J).NE.' '.AND.WX(J+1).EQ.' ') THEN
              IWK0 = IWK0+1
            END IF
          END DO
          IF (IWK0.EQ.2) THEN
            READ (INFILE,*) WN(I),DWK0
            IF (NM0.EQ.1) THEN
              AD(I) = 0
            ELSE
              AD(I) = 1
            END IF
          ELSE IF (IWK0.EQ.3) THEN
            READ (INFILE,*) WN(I),AD(I),DWK0
          ELSE
            CALL ERROR (1,7,*999)
          END IF
        ELSE
          READ (IN,*) WN(I),AD(I),DWK0
        END IF
        IF (AD(I).EQ.0) THEN
          A(I) = DWK0
        ELSE IF (AD(I).EQ.1) THEN
          D(I) = DWK0
          WND = WND+1
        END IF
      END DO
      IF (IWKC(6).EQ.0) THEN
        NP0 = 0
      ELSE
        REWIND (IN)
        DO I=1,SROW(6)-1
          READ (IN,*)
        END DO
        READ (IN,*) CNTWK,NTWK
        NP0 = NTWK
      END IF
      IF (NP0.LT.100) THEN
        IF (NP0.NE.0.AND.NP0.NE.1.AND.NP0.NE.2.AND.
     &      NP0.NE.10.AND.NP0.NE.11.AND.NP0.NE.12.AND.
     &      NP0.NE.20.AND.NP0.NE.21.AND.NP0.NE.22)
     &      CALL ERROR (1,8,*999)
      ELSE
        IF (NP0.NE.100.AND.NP0.NE.101.AND.NP0.NE.102.AND.
     &      NP0.NE.110.AND.NP0.NE.111.AND.NP0.NE.112.AND.
     &      NP0.NE.120.AND.NP0.NE.121.AND.NP0.NE.122)
     &      CALL ERROR (1,8,*999)
      END IF
      IF (NM0.EQ.3.OR.NM0.EQ.4) THEN
        REWIND (11)
        DO I=1,SROW(5)
            READ (IN,*)
        END DO
        DO I=1,NW0
          A(I) = 0.0
          D(I) = 0.0
        END DO
        IWK0 = WND
        IWK1 = 0
        DO I=1,NW0
          IF (AD(I).EQ.0) THEN
            IWK0 = IWK0+1
            READ (IN,*) WN(IWK0),AD(I),A(IWK0)
          ELSE
            IWK1 = IWK1+1
            IF (NM0.EQ.3) THEN
              READ (IN,*) WN(IWK1),AD(I),A(IWK1),D(IWK1)
            ELSE
              READ (IN,*) WN(IWK1),AD(I),D(IWK1)
            END IF
          END IF
        END DO
      END IF
      IF (IWK(0).EQ.1) THEN
        REWIND (11)
        DO I=1,SROW(2)
          READ (IN,*)
        END DO
        DO I=1,NT0
          IF (IWK(I).EQ.1) THEN
            READ (IN,*) IWK0,TNAME(I),VP(I),DWKP(I,1:NT0),T2(I)
          ELSE
            READ (IN,*)
          END IF
        END DO
      END IF
      LMAX = 0
      DO I=1,NT0
        DO J=1,NG0
          DO K=1,NT(J,0)
            IF (NT(J,K).EQ.I) THEN
              LMAX = LMAX+1
            END IF
          END DO
        END DO
      END DO
      CLOSE (11)
      RETURN
  999 CALL INDCT (10)
      RETURN 1
      END SUBROUTINE CNTRL
      SUBROUTINE TMTRX (*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200,L0=L1*L2
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      INTEGER(4)                              :: I,J
      REAL(8),DIMENSION (L1,L1)               :: TM,EC
      REAL(8),DIMENSION (L1)                  :: EV
      INTERFACE
        SUBROUTINE JACOB (NW0,TM,EC,EV,*)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: NW0
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: TM
          REAL(8),INTENT(OUT),DIMENSION (:,:)    :: EC
          REAL(8),INTENT(OUT),DIMENSION (:)      :: EV
        END SUBROUTINE JACOB
        SUBROUTINE MMMLT (T1,T2,T3,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2,T3
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A,B
          REAL(8),INTENT(OUT),DIMENSION (:,:)    :: C
        END SUBROUTINE MMMLT
        SUBROUTINE MVMLT (T1,T2,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A
          REAL(8),INTENT(IN),DIMENSION (:)       :: B
          REAL(8),INTENT(OUT),DIMENSION (:)      :: C
        END SUBROUTINE MVMLT
        SUBROUTINE MTRNS (T1,T2,A)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(INOUT),DIMENSION (:,:)  :: A
        END SUBROUTINE MTRNS
      END INTERFACE
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /MATRX/P,G,F,W,B,C0
      CALL PMTRX (5,*999)
      DO I=1,NW0
        DO J=1,NW0
          TM(I,J) = G(I,J)
        END DO
      END DO
      CALL JACOB (NW0,TM,EC,EV,*999)
      DO I=1,NW0
        IF (EV(I).LT.0.0) THEN
          CALL ERROR (101,11,*30)
   30     RETURN
        END IF
      END DO
      CALL PMTRX (0,*999)
      CALL MMMLT (LMAX,LMAX,NW0,P,G,W)
      CALL PMTRX (4,*999)
      CALL MMMLT (NW0,LMAX,NW0,G,W,F)
      DO I=1,NW0
        DO J=1,NW0
          W(I,J) = EC(I,J)/DSQRT(EV(J))
        END DO
      END DO
      CALL MTRNS (NW0,NW0,W)
      CALL MMMLT (NW0,NW0,NW0,W,F,P)
      CALL MTRNS (NW0,NW0,W)
      CALL MMMLT (NW0,NW0,NW0,P,W,F)
      DO I=1,NW0
        DO J=1,NW0
          TM(I,J) = F(I,J)
        END DO
      END DO
      CALL JACOB (NW0,TM,EC,EV,*999)
      DO I=1,NW0
        IF (EV(I).LT.0.0) THEN
          CALL ERROR (102,11,*90)
          RETURN
        ELSE IF (EV(I).GT.1.0) THEN
          CALL ERROR (103,11,*90)
          RETURN
        END IF
   90 END DO
      RETURN
  999 CALL INDCT (20)
      RETURN 1
      END SUBROUTINE TMTRX
      SUBROUTINE JACOB (T,A,C,E,*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                   :: L1=20
      INTEGER(4),INTENT(IN)                  :: T
      REAL(8),INTENT(IN),DIMENSION (L1,L1)   :: A
      REAL(8),INTENT(OUT),DIMENSION (L1,L1)  :: C
      REAL(8),INTENT(OUT),DIMENSION (L1)     :: E
      INTEGER(4)                             :: P,Q,NW
      REAL(8),DIMENSION (L1,L1)              :: W,W1,W2
      REAL(8)                                :: A1,B1,C1,SI,CO,EMAX
      INTEGER(4)                             :: I,J,K
      REAL(8)                                :: WK1,WK2,WK3,WK4,WK5
      REAL(8),PARAMETER                      :: EPS=1.0E-12
      W = A
      C = 0.0
      DO I=1,T
        C(I,I) = 1.0
      END DO
      NW = 0
   30 NW = NW+1
      EMAX = 0.0
      DO I=1,T
        DO J=I+1,T
          IF (DABS(W(I,J)).GT.EMAX) THEN
            EMAX = DABS(W(I,J))
            P = I
            Q = J
          END IF
        END DO
      END DO
      IF (EMAX.LE.EPS)     GO TO 140
      A1 = (W(P,P)-W(Q,Q))/2.0
      B1 = DSQRT(A1**2+W(P,Q)**2)
      IF (A1.GT.0.0) THEN
        C1 = 1.0
      ELSE
        C1 = -1.0
      END IF
      CO = DSQRT((1.0+C1*A1/B1)/2.0)
      SI = C1*W(P,Q)/(2.0*B1*CO)
      WK1 = W(P,Q)
      WK2 = W(P,P)
      WK3 = W(Q,Q)
      W(P,Q) = 0.0
      W(Q,P) = 0.0
      W(P,P) = WK2*CO**2+2.0*WK1*SI*CO+WK3*SI**2
      W(Q,Q) = WK2*SI**2-2.0*WK1*SI*CO+WK3*CO**2
      DO I=1,T
        IF (I.NE.P.AND.I.NE.Q) THEN
          WK4 = W(P,I)
          WK5 = W(Q,I)
          W(P,I) = WK4*CO+WK5*SI
          W(Q,I) = -WK4*SI+WK5*CO
          W(I,P) = W(P,I)
          W(I,Q) = W(Q,I)
        END IF
      END DO
      DO I=1,T
        DO J=1,T
          IF (I.EQ.J) THEN
            IF (I.EQ.P.OR.I.EQ.Q) THEN
              W1(I,J) = CO
            ELSE
              W1(I,J) = 1.0
            END IF
          ELSE
            IF (I.EQ.P.AND.J.EQ.Q) THEN
              W1(I,J) = -SI
            ELSE IF (I.EQ.Q.AND.J.EQ.P) THEN
              W1(I,J) = SI
            ELSE
              W1(I,J) = 0.0
            END IF
          END IF
        END DO
      END DO
      DO  I=1,T
        DO  J=1,T
          W2(I,J) = 0.0
          DO K=1,T
            W2(I,J) = W2(I,J)+C(I,K)*W1(K,J)
          END DO
        END DO
      END DO
      C = W2
      IF (NW.GT.500)  THEN
        GO TO 900
      END IF
      GO TO 30
  140 DO I=1,T
        E(I) = W(I,I)
      END DO
      RETURN
  900 CALL ERROR (20,0,*999)
  999 RETURN 1
      END SUBROUTINE JACOB
      SUBROUTINE MMMLT (T1,T2,T3,A,B,C)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: T1,T2,T3
      REAL(8),INTENT(IN),DIMENSION (:,:)     :: A,B
      REAL(8),INTENT(OUT),DIMENSION (:,:)    :: C
      INTEGER(4)                             :: I,J,K
      DO I=1,T1
        DO J=1,T3
          C(I,J) = 0.0
          DO K=1,T2
            C(I,J) = C(I,J)+A(I,K)*B(K,J)
          END DO
        END DO
      END DO
      RETURN
      END SUBROUTINE MMMLT
      SUBROUTINE MVMLT (T1,T2,A,B,C)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: T1,T2
      REAL(8),INTENT(IN),DIMENSION (:,:)     :: A
      REAL(8),INTENT(IN),DIMENSION (:)       :: B
      REAL(8),INTENT(OUT),DIMENSION (:)      :: C
      INTEGER(4)                             :: I,K
      DO I=1,T1
        C(I) = 0.0
        DO K=1,T2
          C(I) = C(I)+A(I,K)*B(K)
        END DO
      END DO
      RETURN
      END SUBROUTINE MVMLT
      SUBROUTINE VVMLT (T,A,B,C)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: T
      REAL(8),INTENT(IN),DIMENSION (:)       :: A,B
      REAL(8),INTENT(OUT)                    :: C
      INTEGER(4)                             :: I
      C = 0.0
      DO I=1,T
        C = C+A(I)*B(I)
      END DO
      RETURN
      END SUBROUTINE VVMLT
      SUBROUTINE MTRNS (T1,T2,A)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: T1,T2
      REAL(8),INTENT(INOUT),DIMENSION (:,:)  :: A
      INTEGER(4)                             :: T
      REAL(8)                                :: DWK
      INTEGER(4)                             :: I,J
      IF (T1-T2.GE.0) THEN
        T = T2
      ELSE
        T = T1
      END IF
      DO I=1,T
        DO J=I,T
          DWK = A(I,J)
          A(I,J) = A(J,I)
          A(J,I) = DWK
        END DO
      END DO
      IF (T1.LT.T2) THEN
        DO I=1,T1
          DO J=T+1,T2
            A(J,I) = A(I,J)
          END DO
        END DO
      ELSE
        DO I=T+1,T1
          DO J=1,T
            A(J,I) = A(I,J)
          END DO
        END DO
      END IF
      RETURN
      END
      SUBROUTINE INVMX (N,A,B,*)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: N
      REAL(8),INTENT(IN), DIMENSION (:,:)    :: A
      REAL(8),INTENT(OUT),DIMENSION (:,:)    :: B
      INTEGER(4)            :: I,J,K
      REAL(8)               :: WK
      INTEGER(4)            :: MAXC,IWK
      INTEGER,DIMENSION (N) :: WMAX
      REAL(8),PARAMETER     :: EPS=1.0E-12
      IF (N.LT.1) THEN
        GO TO 910
      ELSE IF (N.EQ.1) THEN
        IF (DABS(A(1,1)).LT.EPS) THEN
          GO TO 900
        END IF
        B(1,1) = 1.0/A(1,1)
        RETURN
      END IF
      DO I=1,N
        DO J=1,N
          B(I,J) = A(I,J)
        END DO
      END DO
      DO K=1,N
        WMAX(K) = K
      END DO
      DO K=1,N
        WK = 0.0
        DO I=K,N
          IF (WK.LT.DABS(B(I,1))) THEN
            WK = DABS(B(I,1))
            MAXC = I
          END IF
        END DO
        IF (WK.LE.EPS) THEN
          GO TO 900
        END IF
        IWK = WMAX(MAXC)
        WMAX(MAXC) = WMAX(K)
        WMAX(K) = IWK
        DO J=1,N
          WK = B(MAXC,J)
          B(MAXC,J) = B(K,J)
          B(K,J) = WK
        END DO
        WK = B(K,1)
        DO J=2,N
          B(K,J-1) = B(K,J)/WK
        END DO
        B(K,N) = 1.0/WK
        DO I=1,N
          IF (I.NE.K) THEN
            WK = B(I,1)
            DO J=2,N
              B(I,J-1) = B(I,J)-WK*B(K,J-1)
            END DO
            B(I,N) = -WK*B(K,N)
          END IF
        END DO
      END DO
      DO K=1,N
        DO J=K,N
          IF (WMAX(J).EQ.K) THEN
            WMAX(J) = WMAX(K)
            DO I=1,N
              WK = B(I,J)
              B(I,J) = B(I,K)
              B(I,K) = WK
            END DO
            GO TO 10
          END IF
        END DO
   10 END DO
      RETURN
  900 CALL ERROR (21,0,*999)
  910 CALL ERROR (22,0,*999)
  999 RETURN 1
      END SUBROUTINE INVMX
      SUBROUTINE PMTRX (Q,*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200,L0=L1*L2
      INTEGER(4)                              :: NT0
      REAL(8),DIMENSION (L2,L2)               :: RB
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      REAL(8),DIMENSION (L1)                  :: VP,H2,T2
      REAL(8),DIMENSION (L1,L1)               :: RG,RP
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      INTEGER(4) Q
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      INTEGER(4)                              :: I,J,K,II,JJ
      INTEGER(4),DIMENSION (L1+1)             :: S
      INTEGER(4),DIMENSION (L1,0:L2)          :: T
      INTEGER(4)                              :: R,C
      INTEGER(4)                              :: IWK0
      INTERFACE
        SUBROUTINE INVMX (T,B,A,*)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: B
          REAL(8),INTENT(OUT),DIMENSION (:,:)    :: A
        END SUBROUTINE INVMX
        SUBROUTINE MTRNS (T1,T2,A)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(INOUT),DIMENSION (:,:)  :: A
        END SUBROUTINE MTRNS
      END INTERFACE
      SAVE /MATRX/
      COMMON /ETC01/NT0
      COMMON /ETC02/RB
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /TRAIT/VP,H2,RG,RP,T2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /MATRX/P,G,F,W,B,C0
      S(1) = 1
      DO I=1,NT0
        S(I+1) = S(I)
        T(I,0) = 0
        DO J=1,NG0
          IF (NA(J).GE.1) THEN
            DO K=1,NT(J,0)
              IF (NT(J,K).EQ.I) THEN
                S(I+1) = S(I+1)+1
                T(I,0) = T(I,0)+1
                T(I,T(I,0)) = J
              END IF
            END DO
          END IF
        END DO
      END DO
      LMAX = S(NT0+1)-1
      IF (Q.EQ.3.OR.Q.EQ.4)     GO TO 300
      IF (Q.EQ.5)               GO TO 500
      DO II=1,NT0
        DO JJ=II,NT0
          IF (II.EQ.JJ) THEN
            DO I=1,T(II,0)
              R = S(II)+I-1
              DO J=I,T(JJ,0)
                C = S(JJ)+J-1
                IF (I.EQ.J) THEN
                  DO K=1,NT(T(II,I),0)
                    IF (NT(T(II,I),K).EQ.II) THEN
                      IWK0 = NR(T(II,I),K)
                    END IF
                  END DO
                  P(R,C) = (1+(IWK0-1)*T2(II))/IWK0
     &                     +(NA(T(II,I))-1)*RW(T(II,I))*H2(II)
                  P(R,C) = P(R,C)/NA(T(II,I))*VP(II)
                ELSE
                  P(R,C) = RB(T(II,I),T(II,J))*VP(II)*H2(II)
                END IF
              END DO
            END DO
          ELSE
            DO I=1,T(II,0)
              R = S(II)+I-1
              DO J=1,T(JJ,0)
                C = S(JJ)+J-1
                IF (T(II,I).EQ.T(JJ,J)) THEN
                  P(R,C) = DSQRT(VP(II)*VP(JJ))*RP(II,JJ)
                  P(R,C) = P(R,C)+(NA(T(II,I))-1)*RW(T(II,I))
     &                     *DSQRT(VP(II)*H2(II)*VP(JJ)*H2(JJ))*RG(II,JJ)
                  P(R,C) = P(R,C)/NA(T(II,I))
                ELSE
                  P(R,C) = DSQRT(VP(II)*H2(II)*VP(JJ)*H2(JJ))*RG(II,JJ)
     &                     *RB(T(II,I),T(JJ,J))
                END IF
              END DO
            END DO
          END IF
        END DO
      END DO
      DO I=1,LMAX
        DO J=I,LMAX
          P(J,I) = P(I,J)
        END DO
      END DO
      IF (Q.EQ.1)     RETURN
      CALL INVMX (LMAX,P,W,*999)
      DO I=1,LMAX
        DO J=1,LMAX
          P(I,J) = W(I,J)
        END DO
      END DO
      IF (Q.EQ.2)     RETURN
  300 DO II=1,NT0
        DO J=1,NW0
          C = J
          IF (II.EQ.WN(J)) THEN
            DO I=1,T(II,0)
              R = S(II)+I-1
              IF (T(II,I).EQ.CAN) THEN
                G(R,C) = VP(II)*H2(II)
              ELSE
                G(R,C) = VP(II)*H2(II)*RB(CAN,T(II,I))
              END IF
            END DO
          ELSE
            DO I=1,T(II,0)
              R = S(II)+I-1
              IF (T(II,I).EQ.CAN) THEN
                G(R,C) = DSQRT(VP(II)*H2(II)*VP(WN(J))*H2(WN(J)))
     &                   *RG(II,WN(J))
              ELSE
                G(R,C) = DSQRT(VP(II)*H2(II)*VP(WN(J))*H2(WN(J)))
     &                   *RG(II,WN(J))*RB(CAN,T(II,I))
              END IF
            END DO
          END IF
        END DO
      END DO
      IF (Q.EQ.0.OR.Q.EQ.3)     RETURN
  400 CALL MTRNS (LMAX,NW0,G)
      IF (Q.EQ.4)               RETURN
  500  DO I=1,NW0
        R = WN(I)
        DO J=1,NW0
          C = WN(J)
          IF (I.EQ.J) THEN
            G(I,J) = VP(R)*H2(R)
          ELSE
            G(I,J) = DSQRT(VP(R)*H2(R)*VP(C)*H2(C))*RG(R,C)
          END IF
        END DO
      END DO
      RETURN
  999 CALL INDCT (100)
      RETURN 1
      END SUBROUTINE PMTRX
      SUBROUTINE CMTRX (*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200,L0=L1*L2
      REAL(8),DIMENSION (L1)                  :: A,D
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      REAL(8),DIMENSION (L1)                  :: VP,H2,T2
      REAL(8),DIMENSION (L1,L1)               :: RG,RP
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      INTEGER(4)                              :: I,J,JUDGE
      INTEGER(4)                              :: R,C,WKC
      REAL(8),DIMENSION (L0,L0)               :: X
      REAL(8),DIMENSION (L0)                  :: V,Z
      REAL(8)                                 :: ARPHA,AMAX,DWK
      INTERFACE
        SUBROUTINE INVMX (T,A,B,*)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A
          REAL(8),INTENT(OUT),DIMENSION (:,:)    :: B
        END SUBROUTINE INVMX
        SUBROUTINE MMMLT (T1,T2,T3,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2,T3
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A,B
          REAL(8),INTENT(OUT),DIMENSION (:,:)    :: C
        END SUBROUTINE MMMLT
        SUBROUTINE MVMLT (T1,T2,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A
          REAL(8),INTENT(IN),DIMENSION (:)       :: B
          REAL(8),INTENT(OUT),DIMENSION (:)      :: C
        END SUBROUTINE MVMLT
        SUBROUTINE VVMLT (T,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T
          REAL(8),INTENT(IN),DIMENSION (:)       :: A,B
          REAL(8),INTENT(OUT)                    :: C
        END SUBROUTINE VVMLT
        SUBROUTINE MTRNS (T1,T2,A)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(INOUT),DIMENSION (:,:)  :: A
        END SUBROUTINE MTRNS
      END INTERFACE
      SAVE /MATRX/
      COMMON /ETC03/A,D
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /TRAIT/VP,H2,RG,RP,T2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /MATRX/P,G,F,W,B,C0
      R = LMAX
      C = WND
      CALL PMTRX (0,*999)
      IF (NM0.EQ.1) THEN
        CALL MVMLT (R,NW0,G,A,V)
        CALL MVMLT (R,R,P,V,B)
        CALL MTRNS (R,NW0,G)
        CALL MVMLT (NW0,R,G,B,D)
      ELSE IF (NM0.EQ.2) THEN
        CALL MMMLT (R,R,NW0,P,G,F)
        CALL PMTRX (4,*999)
        CALL MMMLT (NW0,R,NW0,G,F,W)
        CALL INVMX (NW0,W,X,*999)
        CALL MVMLT (NW0,NW0,X,D,A)
        CALL MVMLT (R,NW0,F,A,B)
      ELSE IF (NM0.EQ.3) THEN
        WKC = C
        JUDGE = 0
        DO I=C,1,-1
          IF (JUDGE.EQ.0.AND.DABS(D(I)).LT.1.0D-10) THEN
            WKC = WKC-1
          ELSE
            JUDGE = 1
          END IF
        END DO
        DO I=1,WKC-1
          DO J=1,WKC-1
            IF (I.EQ.J) THEN
              C0(I,J) = D(WKC)
              W(I,J) = D(WKC)
            ELSE
              C0(I,J) = 0.0
              W(I,J) = 0.0
            END IF
          END DO
        END DO
        DO I=1,WKC-1
          C0(I,WKC) = -D(I)
          W(I,WKC) = -D(I)
          DO J=WKC+1,C
            C0(I,J) = 0.0
            W(I,J) = 0.0
          END DO
        END DO
        IF (WKC.NE.0) THEN
          DO I=WKC,C-1
            DO J=1,C
              C0(I,J) = 0.0
              W(I,J) = 0.0
            END DO
            C0(I,WKC+1) = 1.0
            W(I,WKC+1) = 1.0
          END DO
          WKC = C-1
        ELSE
          DO I=1,C
            DO J=1,C
              C0(I,J) = 0.0
              W(I,J) = 0.0
            END DO
            C0(I,I) = 1.0
            W(I,I) = 1.0
          END DO
          WKC = C
        END IF
        CALL MTRNS (WKC,C,W)
        CALL MMMLT (R,C,WKC,G,W,F)
        CALL MMMLT (R,R,WKC,P,F,W)
        CALL MTRNS (R,WKC,F)
        CALL MMMLT (WKC,R,WKC,F,W,X)
        CALL INVMX (WKC,X,W,*999)
        CALL MMMLT (WKC,WKC,R,W,F,X)
        CALL MTRNS (WKC,R,F)
        CALL MMMLT (R,WKC,R,F,X,W)
        CALL MMMLT (R,R,R,P,W,X)
        DO I=1,R
          DO J=1,R
            IF (I.EQ.J) THEN
              X(I,J) = 1.0-X(I,J)
            ELSE
              X(I,J) = -X(I,J)
            END IF
          END DO
        END DO
        CALL MMMLT (R,R,R,X,P,W)
        CALL MMMLT (R,R,NW0,W,G,X)
        CALL MVMLT (R,NW0,X,A,B)
      ELSE IF (NM0.EQ.4) THEN
        CALL PMTRX (0,*999)
        CALL MMMLT (R,R,C,P,G,W)
        CALL MTRNS (R,NW0,G)
        CALL MMMLT (C,R,C,G,W,X)
        CALL INVMX (C,X,W,*999)
*      *
        CALL MVMLT (C,C,W,D,V)
        CALL VVMLT (C,D,V,AMAX)
        AMAX = DSQRT(1.0/AMAX)
        ARPHA = ERATE*AMAX
        CALL MVMLT (C,C,W,D,V)
        CALL MTRNS (NW0,R,G)
        CALL MVMLT (R,C,G,V,Z)
        CALL MVMLT (R,R,P,Z,V)
        DO I=1,R
          V(I) = V(I)*ARPHA
        END DO
        CALL MMMLT (R,C,C,G,W,X)
        CALL MTRNS (R,NW0,G)
        CALL MMMLT (R,C,R,X,G,W)
        CALL MMMLT (R,R,R,W,P,X)
        DO I=1,R
          DO J=1,R
            IF (I.EQ.J) THEN
              X(I,J) = 1.0-X(I,J)
            ELSE
              X(I,J) = -X(I,J)
            END IF
          END DO
        END DO
        CALL MMMLT (R,R,R,P,X,F)
        CALL MMMLT (NW0,R,R,G,F,W)
        CALL MTRNS (NW0,R,G)
        CALL MMMLT (NW0,R,NW0,W,G,X)
        CALL MVMLT (NW0,NW0,X,A,Z)
        CALL VVMLT (NW0,A,Z,DWK)
        DWK = DSQRT((1.0-ERATE**2)/DWK)
        CALL MVMLT (R,NW0,G,A,Z)
        CALL MVMLT (R,R,F,Z,B)
        DO I=1,R
          B(I) = DWK*B(I)+V(I)
        END DO
        CALL PMTRX (4,*999)
        CALL MVMLT (NW0,R,G,B,Z)
        CALL PMTRX (2,*999)
        CALL MMMLT (NW0,R,R,G,P,W)
        CALL MTRNS (NW0,R,G)
        CALL MMMLT (NW0,R,NW0,W,G,X)
        CALL INVMX (NW0,X,W,*999)
        CALL MVMLT (NW0,NW0,W,Z,A)
        CALL MVMLT (R,NW0,G,A,V)
        CALL MVMLT (R,R,P,V,B)
      END IF
      RETURN
  999 CALL INDCT (40)
      RETURN 1
      END SUBROUTINE CMTRX
      SUBROUTINE RIHDG (*)
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200,L0=L1*L2
      REAL(8),DIMENSION (L1)                  :: A,D
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      REAL(8)                                 :: SI,SH,SIH,RIH
      REAL(8),DIMENSION (L1)                  :: DG
      INTEGER(4)                              :: I
      REAL(8),DIMENSION (L0)                  :: V
      INTERFACE
        SUBROUTINE MVMLT (T1,T2,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T1,T2
          REAL(8),INTENT(IN),DIMENSION (:,:)     :: A
          REAL(8),INTENT(IN),DIMENSION (:)       :: B
          REAL(8),INTENT(OUT),DIMENSION (:)      :: C
        END SUBROUTINE MVMLT
        SUBROUTINE VVMLT (T,A,B,C)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                  :: T
          REAL(8),INTENT(IN),DIMENSION (:)       :: A,B
          REAL(8),INTENT(OUT)                    :: C
        END SUBROUTINE VVMLT
      END INTERFACE
      SAVE /RSLTS/
      COMMON /ETC03/A,D
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /MATRX/P,G,F,W,B,C0
      COMMON /RSLTS/SI,SH,SIH,RIH,DG
      CALL PMTRX (1,*999)
      CALL MVMLT (LMAX,LMAX,P,B,V)
      CALL VVMLT (LMAX,B,V,SI)
      CALL PMTRX (5,*999)
      CALL MVMLT (NW0,NW0,G,A,V)
      CALL VVMLT (NW0,A,V,SH)
      CALL PMTRX (3,*999)
      CALL MVMLT (LMAX,NW0,G,A,V)
      CALL VVMLT (LMAX,B,V,SIH)
      RIH = SIH/DSQRT(SI*SH)
      SI = DSQRT(SI)
      SH = DSQRT(SH)
      CALL PMTRX (4,*999)
      CALL MVMLT (NW0,LMAX,G,B,V)
      DO I=1,NW0
        DG(I) = V(I)/SI
      END DO
      RETURN
  999 CALL INDCT (40)
      RETURN 1
      END SUBROUTINE RIHDG
      SUBROUTINE PROUT (RPT,*)
      IMPLICIT NONE
      INTEGER(4)                              :: RPT
      INTEGER(4),PARAMETER                    :: L1=20,L2=200,L0=L1*L2
      CHARACTER (LEN=80)                      :: TIT
      CHARACTER (LEN=20)                      :: TNAME(L1),GNAME(L2)
      INTEGER(4)                              :: NT0
      REAL(8),DIMENSION (L1)                  :: A,D
      INTEGER(4)                              :: NM0,LMAX
      REAL(8)                                 :: ERATE
      REAL(8),DIMENSION (L1)                  :: VP,H2,T2
      REAL(8),DIMENSION (L1,L1)               :: RG,RP
      INTEGER(4),DIMENSION (L2*2)             :: P1,P2
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4)                              :: NW0,WND
      INTEGER(4),DIMENSION (L1)               :: WN,AD
      INTEGER(4)                              :: NP0
      REAL(8),DIMENSION (L0,L0)               :: P,G,F,W
      REAL(8),DIMENSION (L0)                  :: B
      REAL(8),DIMENSION (L1,L1)               :: C0
      REAL(8)                                 :: SI,SH,SIH,RIH
      REAL(8),DIMENSION (L1)                  :: DG
      INTEGER(4)                              :: NMAT,PNO1,PNO2,PNO3
      INTEGER(4),DIMENSION (4)                :: USEP
      INTEGER(4),DIMENSION (4)                :: PDG0
      INTEGER(4),DIMENSION (4,L1)             :: PDG1,PDG2
      INTEGER(4)                              :: I,J,K,II,IPRNT
      INTEGER(4)                              :: MT
      CHARACTER (LEN=80)                      :: OFILE
      LOGICAL                                 :: TF
      CHARACTER (LEN=1)                       :: YN
      CHARACTER (LEN=66)                      :: TOP
      CHARACTER (LEN=26),DIMENSION (4)        :: MDL
      CHARACTER (LEN=10),DIMENSION (4)        :: CWK1
      CHARACTER (LEN=10),DIMENSION (8)        :: CWK2
      INTEGER(4)                              :: R1,C1,R2,C2,R3,C3
      INTEGER(4)                              :: IWK,IWK0,IWK1
      REAL(8)                                 :: DWK
      REAL(8),DIMENSION (L0)                  :: V
      CHARACTER (LEN=35),DIMENSION (4)        :: CWKP
      COMMON /ETC01/NT0
      COMMON /ETC03/A,D
      COMMON /TITLE/TIT,TNAME,GNAME
      COMMON /MODEL/NM0,LMAX,ERATE
      COMMON /TRAIT/VP,H2,RG,RP,T2
      COMMON /PEDIG/P1,P2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      COMMON /WEIGH/NW0,WN,AD,WND
      COMMON /PRINT/NP0
      COMMON /MATRX/P,G,F,W,B,C0
      COMMON /RSLTS/SI,SH,SIH,RIH,DG
      COMMON /BACKS/NMAT,PNO1,PNO2,PNO3,USEP,PDG0,PDG1,PDG2
      DATA TOP/'  *******   選抜指数の算出のためのプログラム   *******'/
      DATA MDL/'  [ 相対経済重要度 ]       ',
     &         '  [ 相対希望改良量 ]       ',
     &         '(経済重みが既知の場合)     ',
     &         '(経済重みが未知の場合)     '/
      DATA CWK1/'グループ','形質        ','指数の係数','寄与率 (%)'/
      DATA CWK2/'          ','  相対経済 ','  相対希望','遺伝的改良',
     &          '形  質    ','  重要度   ','  改良量  ','量の期待値'/
      DATA CWKP/'選抜候補個体およびその最小血縁情報 ',
     &          '選抜候補個体の血縁係数1/2までの情報',
     &          '選抜候補個体の血縁係数1/4までの情報',
     &          '選抜候補個体の血縁係数1/8までの情報'/
      SAVE OFILE
      IPRNT = 6
      IF (RPT.NE.0)     GO TO 800
   10 IF ((NP0.GE.10.AND.NP0.LE.22).OR.NP0.GE.110) THEN
        WRITE (6,*)
        WRITE (6,*) '出力ファイル名を入力して下さい ----->'
        READ (5,'(A80)') OFILE
        INQUIRE (FILE=OFILE,EXIST=TF)
        IF (TF) THEN
          WRITE (6,*) 'このファイルはすでに存在します．',
     &                '上書きしてよろしいですか ? (Y/N)'
          READ (5,'(A1)') YN
          IF (YN.EQ.'Y'.OR.YN.EQ.'y') THEN
            MT = 21
            OPEN (MT,ERR=900,FILE=OFILE,STATUS='OLD')
            CLOSE (MT,ERR=900,STATUS='DELETE')
            OPEN (MT,ERR=900,FILE=OFILE,STATUS='NEW')
          ELSE IF (YN.EQ.'N'.OR.YN.EQ.'n') THEN
            GO TO 10
          ELSE
            CALL ERROR (1,0,*10)
          END IF
        ELSE
          MT = 21
          OPEN (MT,ERR=900,FILE=OFILE,STATUS='NEW')
        END IF
        IF (NG0.EQ.0) THEN
          CLOSE (MT)
        END IF
      ELSE
        MT = 6
      END IF
      IF (NG0.EQ.0)     GO TO 800
   20 IF (NP0.GE.100) THEN
        C1 = 1
        R1 = LMAX
        C2 = 1
        R2 = NW0
        C3 = 1
        R3 = NT0
      ELSE
        IF (MOD(LMAX,IPRNT).NE.0) THEN
          C1 = LMAX/IPRNT
          R1 = LMAX-C1*IPRNT
          C1 = C1+1
        ELSE
          C1 = LMAX/IPRNT
          R1 = IPRNT
        END IF
        IF (MOD(NW0,IPRNT).NE.0) THEN
          C2 = NW0/IPRNT
          R2 = NW0-C2*IPRNT
          C2 = C2+1
        ELSE
          C2 = NW0/IPRNT
          R2 = IPRNT
        END IF
        IF (MOD(NT0,IPRNT).NE.0) THEN
          C3 = NT0/IPRNT
          R3 = NT0-C3*IPRNT
          C3 = C3+1
        ELSE
          C3 = NT0/IPRNT
          R3 = IPRNT
        END IF
      END IF
      WRITE (MT,3000) TOP
      WRITE (MT,3300)
      CALL PAUSE (3,MT)
      WRITE (MT,3100) TIT
      WRITE (MT,3300)
      CALL PAUSE (3,MT)
      IF (NM0.EQ.1) THEN
        WRITE (MT,3200) 'モデル  =  Smith-Hazel型選抜指数 ',MDL(NM0)
      ELSE
        WRITE (MT,3200) 'モデル  =  制限付き選抜指数 ',MDL(NM0)
      END IF
      WRITE (MT,3400)
      CALL PAUSE (4,MT)
      IF (MOD(NP0,10).NE.0)     GO TO 200
      WRITE (MT,3100) '< 入力行列 >'
      WRITE (MT,3100) '  個体ベースにおける表型分散共分散行列'
      CALL PAUSE (1,MT)
      DO II=1,C3
        IF (II.EQ.C3) THEN
          IWK = R3
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,NT0
          WRITE (MT,2000) (RP(I,J)*DSQRT(VP(I)*VP(J)),J=1,IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,*)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
      WRITE (MT,3100) '  個体ベースにおける遺伝分散共分散行列'
      CALL PAUSE (1,MT)
      DO II=1,C3
        IF (II.EQ.C3) THEN
          IWK = R3
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,NT0
          WRITE (MT,2000) (RG(I,J)*DSQRT(H2(I)*H2(J)*VP(I)*VP(J)),
     &                     J=1,IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,*)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3000) '  遺伝的パラメーター行列 = 正定値行列(正常行列)'
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
      IF (NT0.GE.2) THEN
        WRITE (MT,3100) '  表型相関係数行列'
        CALL PAUSE (1,MT)
        DO II=1,C3
          IF (II.EQ.C3) THEN
            IWK = R3
          ELSE
            IWK = IPRNT
          END IF
          DO I=1,NT0
            WRITE (MT,2000) RP(I,1:IWK)
            CALL PAUSE (1,MT)
          END DO
          WRITE (MT,*)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,3400)
        CALL PAUSE (2,MT)
        WRITE (MT,3100) '  遺伝相関係数行列'
        CALL PAUSE (1,MT)
        DO II=1,C3
          IF (II.EQ.C3) THEN
            IWK = R3
          ELSE
            IWK = IPRNT
          END IF
          DO I=1,NT0
            WRITE (MT,2000) RG(I,1:IWK)
            CALL PAUSE (1,MT)
          END DO
          WRITE (MT,*)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,3400)
        CALL PAUSE (2,MT)
      END IF
  200 IF (MOD(NP0,10).EQ.2)     GO TO 400
      WRITE (MT,3100) '< 算出された行列 >'
        CALL PAUSE (2,MT)
      CALL PMTRX (1,*999)
      WRITE (MT,3100) '  表型分散共分散行列'
      CALL PAUSE (1,MT)
      DO II=1,C1
        IF (II.EQ.C1) THEN
          IWK = R1
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,LMAX
          WRITE (MT,2000) P(I,1:IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,*)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
  230 CALL PMTRX (2,*999)
      WRITE (MT,3100) '  表型分散共分散行列の逆行列'
      CALL PAUSE (2,MT)
      DO II=1,C1
        IF (II.EQ.C1) THEN
          IWK = R1
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,LMAX
          WRITE (MT,2000) P(I,1:IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,3300)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
  260 CALL PMTRX (3,*999)
      WRITE (MT,3100) '  表型価と育種価の共分散行列'
      CALL PAUSE (2,MT)
      DO II=1,C2
        IF (II.EQ.C2) THEN
          IWK = R2
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,LMAX
          WRITE (MT,2000) G(I,1:IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,3300)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
  290 CALL PMTRX (5,*999)
      WRITE (MT,3100) '  遺伝分散共分散行列'
      CALL PAUSE (2,MT)
      DO II=1,C2
        IF (II.EQ.C2) THEN
          IWK = R2
        ELSE
          IWK = IPRNT
        END IF
        DO I=1,NW0
          WRITE (MT,2000) G(I,1:IWK)
          CALL PAUSE (1,MT)
        END DO
        WRITE (MT,3300)
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
  400 WRITE (MT,3100) '< 選抜指数 >'
      WRITE (MT,4000) CWK1(1:4)
      CALL PAUSE (1,MT)
      WRITE (MT,3300)
      CALL PAUSE (1,MT)
      CALL PMTRX (2,*999)
      DO I=1,LMAX
        DWK = DSQRT(SI**2-B(I)**2/P(I,I))
        V(I) =(1.0-DWK/SI)*100.0
      END DO
      IWK = 0
      DO I=1,NT0
        DO J=1,NG0
          IF (NA(J).GE.1) THEN
            DO K=1,NT(J,0)
              IF (NT(J,K).EQ.I) THEN
                IWK = IWK+1
                WRITE (MT,2100) GNAME(J),TNAME(I),B(IWK),V(IWK)
                CALL PAUSE (1,MT)
              END IF
            END DO
          END IF
        END DO
      END DO
      WRITE (MT,3300)
      WRITE (MT,3400)
      CALL PAUSE (3,MT)
  500 IF (NM0.EQ.4) THEN
        WRITE (MT,2600) 'E-値   = ',ERATE
        WRITE (MT,3300)
      END IF
      CALL PAUSE (2,MT)
      WRITE (MT,4100) CWK2(1),CWK2(2),CWK2(3),CWK2(4)
      WRITE (MT,4100) CWK2(5),CWK2(6),CWK2(7),CWK2(8)
      CALL PAUSE (2,MT)
      WRITE (MT,3300)
      CALL PAUSE (1,MT)
      IWK0 = WND
      IWK1 = 0
      DO I=1,NW0
        IF (AD(I).EQ.0) THEN
          IWK0 = IWK0+1
          IWK = IWK0
        ELSE
          IWK1 = IWK1+1
          IWK = IWK1
        END IF
        IF (NM0.EQ.1) THEN
          WRITE (MT,2200) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
        ELSE IF (NM0.EQ.2) THEN
          WRITE (MT,2300) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
        ELSE
          IF (NM0.EQ.3.AND.AD(I).EQ.1) THEN
            WRITE (MT,2400) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
          ELSE IF (AD(I).EQ.0) THEN
            WRITE (MT,2500) TNAME(WN(IWK)),A(IWK),DG(IWK)
          ELSE
            WRITE (MT,2300) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
          END IF
        END IF
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3300)
      WRITE (MT,3400)
      CALL PAUSE (3,MT)
      WRITE (MT,2600) '選抜指数(I)の標準偏差               = ',SI
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '総合育種価(H)の標準偏差             = ',SH
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '選抜指数と総合育種価の共分散        = ',SIH
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '選抜の正確度                        = ',RIH
      CALL PAUSE (1,MT)
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
      IF (MT.EQ.21) THEN
        CLOSE (21,ERR=910)
        IF (NP0.LE.12.OR.(NP0.GE.100.AND.NP0.LE.112)) THEN
          MT = 6
          GO TO 20
        END IF
      END IF
  800 IF (RPT.EQ.0)     RETURN
      IF ((NP0.GE.10.AND.NP0.LE.22).OR.NP0.GE.110) THEN
        MT = 21
        OPEN (MT,ERR=900,FILE=OFILE,STATUS='OLD')
      ELSE
        MT = 6
      END IF
  810 IF (MT.EQ.21) THEN
        REWIND (MT)
        IWK0 = 0
        DO
          READ (MT,*,END=820)
          IWK0 = IWK0+1
        END DO
  820   REWIND (MT)
        DO I=1,IWK0
          READ (MT,*)
        END DO
      END IF
      WRITE (MT,3300)
      IF (RPT.EQ.1) THEN
        WRITE (MT,3100) '< 家系指数による改良量の予測 >'
      END IF
      WRITE (MT,3100) CWKP(RPT)
      WRITE (MT,4100) CWK2(1),CWK2(2),CWK2(3),CWK2(4)
      WRITE (MT,4100) CWK2(5),CWK2(6),CWK2(7),CWK2(8)
      CALL PAUSE (2,MT)
      WRITE (MT,3300)
      CALL PAUSE (1,MT)
      IWK0 = WND
      IWK1 = 0
      DO I=1,NW0
        IF (AD(I).EQ.0) THEN
          IWK0 = IWK0+1
          IWK = IWK0
        ELSE
          IWK1 = IWK1+1
          IWK = IWK1
        END IF
        IF (NM0.EQ.1) THEN
          WRITE (MT,2200) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
        ELSE IF (NM0.EQ.2) THEN
          WRITE (MT,2300) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
        ELSE
          IF (NM0.EQ.3.AND.AD(I).EQ.1) THEN
            WRITE (MT,2400) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
          ELSE IF (AD(I).EQ.0) THEN
            WRITE (MT,2500) TNAME(WN(IWK)),A(IWK),DG(IWK)
          ELSE
            WRITE (MT,2300) TNAME(WN(IWK)),A(IWK),D(IWK),DG(IWK)
          END IF
        END IF
        CALL PAUSE (1,MT)
      END DO
      WRITE (MT,3300)
      WRITE (MT,3400)
      CALL PAUSE (3,MT)
      WRITE (MT,2600) '選抜指数(I)の標準偏差               = ',SI
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '総合育種価(H)の標準偏差             = ',SH
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '選抜指数と総合育種価の共分散        = ',SIH
      CALL PAUSE (1,MT)
      WRITE (MT,2600) '選抜の正確度                        = ',RIH
      CALL PAUSE (1,MT)
      WRITE (MT,3400)
      CALL PAUSE (2,MT)
      IF (MT.EQ.21) THEN
        CLOSE (21,ERR=910)
        IF (NP0.LE.12.OR.(NP0.GE.100.AND.NP0.LE.112)) THEN
          MT = 6
          GO TO 810
        END IF
      END IF
      RETURN
  900 CALL ERROR (11,21,*999)
  910 CALL ERROR (12,21,*999)
  999 CALL INDCT (60)
      RETURN 1
 2000 FORMAT (1H 6F12.5)
 2100 FORMAT (1H A20,5X,A20,2F15.6)
 2200 FORMAT (1H A20,5X,F15.6,F15.6,'+',F14.6)
 2300 FORMAT (1H A20,5X,F15.6,'+',F14.6,F15.6)
 2400 FORMAT (1H,A20,5X,3F15.6)
 2500 FORMAT (1H A20,5X,F15.6,'        -----  ',F15.6)
 2600 FORMAT (1H A,F17.6)
 3000 FORMAT (1H /A/)
 3100 FORMAT (1H A/)
 3200 FORMAT (1H 2A/)
 3300 FORMAT (1H )
 3400 FORMAT (1H /)
 4000 FORMAT (1H (A10,10X),5X,(A10,10X),2(5X,A10))
 4100 FORMAT (1H (A10,10X),5X,3(5X,A10))
      END SUBROUTINE PROUT
      SUBROUTINE PAUSE (Q,MT)
      IMPLICIT NONE
      INTEGER(4)                              :: Q,MT
      INTEGER(4)                              :: PS,PMAX
      INTEGER(4)                              :: I
      SAVE PS
      PMAX = 20
      IF (MT.NE.6)     RETURN
      DO I=1,Q
        PS = PS+1
        IF (PS.GE.PMAX) THEN
          PAUSE '     ----- リターンキーを押して下さい -----'
          WRITE (MT,*)
          PS = 0
        END IF
      END DO
      RETURN
      END SUBROUTINE PAUSE
      SUBROUTINE INDCT (Q)
      IMPLICIT NONE
      INTEGER(4)                              :: Q
      WRITE (6,*)
      IF (Q.EQ.10) THEN
        WRITE (6,*) '< コントロールファイルの読み込み時 >'
      ELSE IF (Q.EQ.20) THEN
        WRITE (6,*) '< 遺伝的パラメーターの検定時 >'
      ELSE IF (Q.EQ.40) THEN
        WRITE (6,*) '< 選抜指数における係数の算出時 >'
      ELSE IF (Q.EQ.50) THEN
        WRITE (6,*) '< 選抜の正確度の算出および遺伝的改良量の推定時 >'
      ELSE IF (Q.EQ.60) THEN
        WRITE (6,*) '< 結果の出力時 >'
      ELSE IF (Q.EQ.100) THEN
        WRITE (6,*) '< 選抜指数におけるパラメーター行列の作出時 >'
      END IF
      WRITE (6,*)
      RETURN
      END SUBROUTINE INDCT
      SUBROUTINE ERROR (ERR1,ERR2,*)
      IMPLICIT NONE
      INTEGER(4)                              :: ERR1,ERR2
      IF (ERR1.LT.100) THEN
        WRITE (6,*) '*****   エラー   *****   エラー番号 = ',ERR1
      ELSE
        WRITE (6,*) '*****   警  告   *****   警告番号 = ',ERR1
      END IF
      WRITE (6,*)
      IF (ERR1.EQ.1) THEN
        WRITE (6,*) '単純なエラーです.'
      ELSE IF (ERR1.EQ.2) THEN
        WRITE (6,*) 'プログラムの制限を越えています.'
      ELSE IF (ERR1.EQ.11) THEN
        WRITE (6,*) 'ファイルを開くことができません.'
      ELSE IF (ERR1.EQ.12) THEN
        WRITE (6,*) 'ファイルを生成することができません.'
      ELSE IF (ERR1.EQ.20) THEN
        WRITE (6,*) '行列の固有値が算出できません.'
      ELSE IF (ERR1.EQ.21) THEN
        WRITE (6,*) '行列が特異で，逆行列が算出できません.'
      ELSE IF (ERR1.EQ.22) THEN
        WRITE (6,*) '行列のサイズが異常です.'
      ELSE IF (ERR1.EQ.101) THEN
        WRITE (6,*) '遺伝分散共分散行列がパラメーター空間外です.'
        WRITE (6,*)
        PAUSE '     ----- リターンキーを押して下さい -----'
        WRITE (6,*)
      ELSE IF (ERR1.EQ.102) THEN
        WRITE (6,*) '遺伝的パラメーターがパラメーター空間外です.'
        WRITE (6,*)
        PAUSE '     ----- リターンキーを押して下さい -----'
        WRITE (6,*)
      ELSE IF (ERR1.EQ.103) THEN
        WRITE (6,*) '遺伝的パラメーターがパラメーター空間外です.'
        WRITE (6,*)
        PAUSE '     ----- リターンキーを押して下さい -----'
        WRITE (6,*)
      END IF
      IF (ERR2.EQ.1) THEN
        WRITE (6,*) 'コントロールファイルに不備があります.'
      ELSE IF (ERR2.EQ.2) THEN
        WRITE (6,*) 'コントロールファイルの''TITLE''文のエラーです.'
      ELSE IF (ERR2.EQ.3) THEN
        WRITE (6,*) 'コントロールファイルの''MODEL''文のエラーです.'
      ELSE IF (ERR2.EQ.4) THEN
        WRITE (6,*) 'コントロールファイルの''TRAIT''文のエラーです.'
      ELSE IF (ERR2.EQ.5) THEN
        WRITE (6,*) 'コントロールファイルの''GROUP''文のエラーです.'
      ELSE IF (ERR2.EQ.6) THEN
        WRITE (6,*) 'コントロールファイルの''PEDIG''文のエラーです.'
      ELSE IF (ERR2.EQ.7) THEN
        WRITE (6,*) 'コントロールファイルの''WEIGHT''文のエラーです.'
      ELSE IF (ERR2.EQ.8) THEN
        WRITE (6,*) 'コントロールファイルの''PRINT''文のエラーです.'
      ELSE IF (ERR2.EQ.11) THEN
        WRITE (6,*) '遺伝的パラメーターの設定を見直しして下さい.'
      ELSE IF (ERR2.EQ.21) THEN
        WRITE (6,*) '出力ファイルに不備があります.'
      END IF
      RETURN 1
      END SUBROUTINE ERROR
      SUBROUTINE AMTRX
      IMPLICIT NONE
      INTEGER(4),PARAMETER                    :: L1=20,L2=200
      REAL(8),DIMENSION (L2,L2)               :: RB
      INTEGER(4),DIMENSION (L2*2)             :: P1,P2
      INTEGER(4)                              :: NG0,CAN
      INTEGER(4),DIMENSION (L2)               :: NA
      INTEGER(4),DIMENSION (L2,0:L1)          :: NT,NR
      REAL(8),DIMENSION (L2)                  :: RW
      INTEGER(4),DIMENSION (L2*2)             :: IN
      REAL(8),DIMENSION (0:L2*2)              :: F0
      INTEGER(4)                              :: NB,I,J
      INTERFACE
        SUBROUTINE INSORT (NA,NB,IN,P1,P2)
          INTEGER(4),INTENT(IN)                   :: NA
          INTEGER(4),INTENT(OUT)                  :: NB
          INTEGER(4),INTENT(INOUT),DIMENSION (:)  :: IN,P1,P2
        END SUBROUTINE INSORT
        SUBROUTINE INBCOE (NA,NS,P1,P2,F0)
          INTEGER(4),INTENT(IN)                   :: NA,NS
          INTEGER(4),INTENT(INOUT),DIMENSION (:)  :: P1,P2
          REAL(8),INTENT(INOUT),DIMENSION (0:)    :: F0
        END SUBROUTINE INBCOE
      END INTERFACE
      SAVE /ETC02/,/PEDIG/,/GROUP/
      COMMON /ETC02/RB
      COMMON /PEDIG/P1,P2
      COMMON /GROUP/NG0,CAN,NA,NT,NR,RW
      DO I=1,NG0
        IN(I) = I
      END DO
      CALL INSORT (NG0,NB,IN,P1,P2)
      CALL INBCOE (NG0+NB,0,P1,P2,F0)
        DO I=NB+1,NB+NG0-1
          DO J=I+1,NB+NG0
            P1(NG0+NB+1) = I
            P2(NG0+NB+1) = J
            CALL INBCOE (NG0+NB+1,NG0+NB,P1,P2,F0)
            RB(I-NB,J-NB) = F0(NB+NG0+1)*2
          END DO
        END DO
        DO I=1,NG0
          DO J=I,NG0
            IF (I.EQ.J) THEN
              RB(I,J) = 1.0
            ELSE
              RB(J,I) = RB(I,J)
            END IF
          END DO
        END DO
      DO I=1,NG0
        IF (NA(I).EQ.1) THEN
          RW(I) = 1.0
        ELSE IF (NA(I).GT.1) THEN
          IF (P1(I).EQ.0.AND.P2(I).EQ.0) THEN
            RW(I) = 0.0
          ELSE IF (P1(I).NE.0.AND.P2(I).NE.0) THEN
            RW(I) = 0.5
          ELSE
            RW(I) = 0.25
          END IF
        END IF
      END DO
      RETURN
      END SUBROUTINE AMTRX
      SUBROUTINE INSORT (NA,NB,IN,P1,P2)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: NA
      INTEGER(4),INTENT(OUT)                 :: NB
      INTEGER(4),INTENT(INOUT),DIMENSION (:) :: IN,P1,P2
      INTEGER(4),PARAMETER        :: L1=20,L2=200,L0=L1+L2
      INTEGER(4)                  :: I
      INTEGER(4)                  :: N
      INTEGER(4),DIMENSION (L0)   :: WK
      INTEGER(4),DIMENSION (0:L0) :: WK0
      INTEGER(4)                  :: WNB,IWK
      INTEGER(4)                  :: WB
      INTEGER(4),DIMENSION (L0*2) :: B
      INTERFACE
        SUBROUTINE IQSORT (I,J,K,A,B)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)                    :: I,J,K
          INTEGER(4),INTENT(INOUT),DIMENSION (:)   :: A,B
        END SUBROUTINE IQSORT
        INTEGER(4) FUNCTION IFIND (ID,N,IN)
          IMPLICIT NONE
          INTEGER(4),INTENT(IN)               :: N
          INTEGER(4),INTENT(IN)               :: ID
          INTEGER(4),INTENT(IN),DIMENSION (:) :: IN
        END FUNCTION IFIND
      END INTERFACE
      NB = 0
      WK0(0) = 0
      DO I=1,NA
        WK0(I) = I
        WK(I) = I
      END DO
      CALL IQSORT (1,NA,2,IN,WK)
      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,B)
      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
      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) = NB+I
      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
      CALL IQSORT (1,N,2,WK,IN)
      DO I=1,NB
        IN(I) = B(I)
      END DO
      RETURN
      END SUBROUTINE INSORT
      RECURSIVE SUBROUTINE IQSORT (I,J,N,A,B)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                  :: I,J,N
      INTEGER(4),INTENT(INOUT),DIMENSION (:) :: A,B
      INTEGER(4) :: IPV,KPV
      INTEGER(4) :: PV
      INTEGER(4) :: S0,S1,S2,K,H
      INTEGER(4) :: LET,RIT
      INTEGER(4) :: 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
      INTEGER(4) FUNCTION IFIND (ID,N,IN)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN) :: N
      INTEGER(4),INTENT(IN) :: ID
      INTEGER(4),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
        END IF
      END DO
      RETURN
      END FUNCTION IFIND
      SUBROUTINE INBCOE (NA,NS,P1,P2,F)
      IMPLICIT NONE
      INTEGER(4),INTENT(IN)                    :: NA,NS
      INTEGER(4),INTENT(INOUT),DIMENSION (:)   :: P1,P2
      REAL(8),INTENT(INOUT),DIMENSION (0:)     :: F
      INTEGER(4)                    :: I,J,K
      INTEGER(4)                    :: IS,ID,KS,KD,IWK
      INTEGER(4),DIMENSION (NA)     :: POINT
      REAL(8),DIMENSION (NA)        :: LM,D
      REAL(8)                       :: FI,R
      F(0) = -1.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
      END DO
      RETURN
      END SUBROUTINE INBCOE
