      PROGRAM LIPED
C         Program LIPED  Microcomputer version (April 1985)
C                        for IBM PC/XT, single precision
C
C         LOGICAL UNITS -  5  input file, named at execution time
C                          6  output file, named at execution time
C                      INOUT  =0, input and output through keyboard/screen
C                     11, 12  scratch files, automatically created and deleted
C                             when more than one marker locus is present
C
      CHARACTER*60 TEXT
      CHARACTER*80 FOR,FT1,FT2
      COMMON /B/ LIST(50),NM(200),NC(200),NS(200),M,L,NOPAR,ID(200),
     1  NDI,LDI(2,3),LGC(3),ISC,MIST,ISEX(200),
     1  FENO1(200,10),FENO2(200,21),PHI(300),RMUT,
     1  HOLD(3,300),RLIK,PHIS(36),WVK(4),THM,THF,IAD(36,36),
     2  IV(4),JV(4),IVA1(4),IVA2(16),IVA3(4),IVF(4),NGEN,NGENO,MNP,
     3  MNPT,KONT1,KONT2,NAL1,NAL2,MVEC,ICH,ISLL,IRESC
      COMMON GENO(25,300),LGENO(200),IGENO(25)
      DIMENSION GF1(4),GF2(6),NAL(10),NF(10),KONT(10),
     1  IAU(10),THV1(10),THV2(10),THVS(10),
     1  GVX1(20,10),GVX2(20,21),THV(3,5),PHPROB(20)
      CHARACTER*4 NM,NC,NOPAR,ID,SEX,BOY,GIRL,RNAM,RNAM1,G1,G2,GEN(6)
      CHARACTER*8 PHE1(20),PHE2(20),PHEPED(10),X1,X2,UNK1,UNK2,UNK(10)
      LOGICAL ISLL,ISEX,MIST,ICH,MARC1,IGENO,IRESC
C         ISLL=.TRUE. IF NO SEX-LINKAGE
C         ISEX(I)=.TRUE. IF PERSON I IS FEMALE
      DATA IPR,K,DSTD/0,0,12./
C Define file for keyboard input and screen output as unit 0
      OPEN (5,FILE='for001.dat',status='old')
      OPEN(6,FILE='for002.dat',STATUS='NEW')
      WRITE(*,41) ' Working ...'
    1 FORMAT(2I2,3I1,F13.0,A60)
    2 FORMAT(//8H PROBLEM,I3,3X,A60/1X,10(1H*))
    3 FORMAT(A80)
    4 FORMAT(40I2)
    5 FORMAT(9X,'Program LIPED  Version for IBM PC/XT',
     1 '   April 1985    J. Ott'//1X,71(1H-))
    6 FORMAT(2I4,A60)
    7 FORMAT(//9H PEDIGREE,I3,3X,A60/1X,11(1H-)/15X,I4,8H PERSONS)
    8 FORMAT(/32H ERROR. P(X/G) OTHER THAN 0 OR 1,
     1 29H SPECIFIED WITH KONT=0, LOCUS,I3,3X,A4)
    9 FORMAT(/18H ERROR. PERSON NR.,I4,2X,A4,4X,
     1 30HNO MATCH OF PHENOTYPE AT LOCUS,I3,3X,A4)
   10 FORMAT(10F8.4)
   11 FORMAT(/29H BLANK CARD ENCOUNTERED. STOP)
c---recomb-frac   før 2F3.2
   12 FORMAT(2F4.3)
   13 FORMAT(//6H LOCUS,I3,3X,A4,5X,3HVS.,5X,5HLOCUS,I3,3X,A4/
     1 1X,43(1H-))
   14 FORMAT(/42H  R MALE  R FEM.   LOG10(L(R))   LOD-SCORE)
   15 FORMAT(1X,2F7.4,3X,F12.5,2X,F10.5)
   16 FORMAT(//7H R MALE/1X,A4)
   17 FORMAT(/1X,F6.4,8(2X,F10.5))
   18 FORMAT(/5X,8F12.4)
   19 FORMAT(40X,8HR FEMALE,3X,A4)
   20 FORMAT(/17H ERROR IN PROBLEM,I3,3X,5HINPUT)
   21 FORMAT(/36H WARNING. SUM OF GENE FREQUENCIES IS,F7.4,7H. LOCUS,
     1 I3,3X,A4)
   22 FORMAT(/47H NO MATCH IN PHENOTYPE/GENOTYPE TABLE FOR LOCUS,
     1 I3,3X,A4/6X,4HLINE,I3,7H ALLELE,I2)
   23 FORMAT( 21H GENE FREQUENCIES FOR,I3,3X,A4,5F8.4/(31X,5F8.4))
   24 FORMAT(/18H ERROR. PERSON NR.,I4,3X,A4,5X,
     1 23HFEMALE SEX INCONSISTENT )
   25 FORMAT(/45H ERROR. INCORRECT NUMBER OF PERSONS SPECIFIED,I6)
   26 FORMAT(4H L =,I5)
   27 FORMAT(5H LIST,15I4/(5X,15I4))
   28 FORMAT(3H NS,I6,14I4/(5X,15I4))
   29 FORMAT(/31H NO MATCH OF DOUBLED INDIVIDUAL,2X,A4)
   30 FORMAT(/25H ERROR. NDI OUT OF BOUNDS,I6)
   31 FORMAT(/47H ERROR. INCREASE DIMENSION MLIST IN LIST(MLIST),I7)
   32 FORMAT(/36H NO MORE THAN 1 PEDIGREE PER PROBLEM,
     1 24H WITH MORE THAN 1 MARKER)
   33 FORMAT(/30H LIKELIHOOD IS ZERO AT  R = .5)
   34 FORMAT(/42H TOO MANY REC.FRACTION VALUES AFTER CARD 9)
   35 FORMAT(/35H IAU(1)=9 ALLOWED ONLY WITH NMARK=1)
   36 FORMAT(//8H SUMMARY/1X,7(1H-)//30H  R MALE  R FEM.   SUM OF LODS)
   37 FORMAT(/21H HAPLOTYPE FREQU. FOR,I3,2X,A4,4H AND,I3,2X,A4/
     1  (6X,8F8.4))
   38 FORMAT(/34H IHAPL=1 ALLOWED ONLY WITH NMARK=1)
   39 FORMAT(/' Number of alleles exceeds MNAL',I1)
   40 FORMAT(/' Number of phenotypes exceeds MNF',I1)
   41 FORMAT(A)
C
   91 MLIST=50
      MMARK=10
      MNAL1=4
      MNAL2=6
      MNDI=3
      MNF1=20
      MNF2=20
      MNP=25
      MNPT=200
c  MT=antal af thetav‘rdier -1  hvis mt>10=>en lodscore gal h›jere end 3
   99 MT=10
      MNAL=MAX0(MNAL1,MNAL2)
C
      THF=0.
      THM=-.10
      R=-.10
      DO 100 I=1,4
      THF=THF+.10
      THM=THM+.15
      R=R+.20
      THV(2,I)=THM
      THV(3,I)=THF
  100 THV(1,I)=R
      THV(3,5)=.50
C
C        THV(1,I) = .10, .30, .50
C        THV(2,I) = .05, .20, .35, .50
C        THV(3,I) = .10, .20, .30, .40, .50
C
      M=MNAL*MNAL
      DO 102 I=1,M
      DO 102 J=1,I
      K=K+1
      IAD(I,J)=K
  102 IAD(J,I)=K
      WRITE(6,5)
C
C  READ PROBLEM INFORMATION
C
C  NMARK = NUMBER OF MARKER LOCI (IN ADDITION TO MAIN LOCUS)
C  MVEC = 1  TO PRINT MODIFIED VECTORS OF PARENTS IN PEDIGREE,
C            L, AND CONTENT OF VECTORS LIST AND NS
C       = 2 TO PRINT L AND CONTENT OF LIST AND NS
C  ISC  = 2 RESCALE LIKELIHOOD TO PREVENT UNDERFLOWS (fixed value in
C           single precision version)
C  ISL  = 1   FOR SEX-LINKAGE
C  IHAPL=1 TO READ HAPLOTYPE FREQUENCIES
C  RMUT= MUTATION RATE
  110 READ(5,1) NMARK,MVEC,ISC,ISL,IHAPL,RMUT,TEXT
      ISC=2
      DO 111 I=1,MT
  111 THVS(I)=0.
      ISC=ISC-1
      IPR=IPR+1
      IF(NMARK.GT.MMARK) GO TO 113
      GO TO 115
  113 WRITE(6,20) IPR
      STOP
  115 MIST=.FALSE.
      IF(IHAPL*NMARK-1) 117, 117, 116
  116 WRITE(6,38)
      STOP
  117 IF(NMARK-1) 113,120,118
  118 OPEN(11,FORM='UNFORMATTED')
      OPEN(12,FORM='UNFORMATTED')
  120 ISLL=.TRUE.
      IF(ISL.EQ.1) ISLL=.FALSE.
      WRITE(6,2) IPR,TEXT
C
C  FORMAT CARD FOR NAME OF LOCUS, ALLELE SYMBOLS, PHENOTYPE SYMBOLS
C       (E.G. 20A4)
      READ(5,3) FT1
C
C  FORMAT CARD FOR 2 ALLELE SYMBOLS, PROBABILITIES IN TABLE P(X/V)
C       (E.G. 2A4,36F2.0)  linie 3
      READ(5,3) FT2
C
C  FORMAT CARD FOR PEDIGREE DATA (STATEMENT 240)
C       (E.G. 20A4)  linie 4
      READ(5,3) FOR
C
      READ(5,FOR)  G1,G2,NOPAR,BOY,UNK1,(UNK(I),I=1,NMARK)
C         NOPAR = SYMBOL FOR NO PARENT (E.G. BLANK)
C         BOY = SYMBOL FOR MALE. EVERYTHING ELSE IS CONSIDERED FEMALE
C             (THE 2 VARIABLES TO THE LEFT OF NOPAR JUST SERVE TO FILL
C             THE SPACE IN FORMAT  FOR)
      READ(5,4) NAL1,(NAL(I),I=1,NMARK)
      READ(5,4) NF1,(NF(I),I=1,NMARK)
      READ(5,4) KONT1,(KONT(I),I=1,NMARK)
      KAUT=-1
  132 READ(5,4) (IAU(I),I=1,NMARK)
      KAUT=KAUT+1
      IF(KAUT.EQ.0) GO TO 134
      IF(IAU(1).EQ.80) GO TO 110
      IF(IAU(1).EQ.90) STOP
      M=MM
      NDI=NNDI
      REWIND 11
      DO 133 K=1,NMARK
      NAL1=NAL(K)
      NF1=NF(K)
      NG1=(NAL1*(NAL1+1))/2
      READ(11) RNAM1,(GF1(I),I=1,NAL1),(PHE1(I),I=1,NF1),
     1  ((GVX1(J,I),I=1,NG1),J=1,NF1)
      IF(IAU(K).EQ.0) GO TO 133
      KONT1=KONT(K)
      UNK1=UNK(K)
      NEW1=K
      MARK=K+1
      REWIND 12
      NAL2=NAL(MARK)
      NF2=NF(MARK)
      KONT2=KONT(MARK)
      UNK2=UNK(MARK)
      NG2=(NAL2*(NAL2+1))/2
      NGEN=NAL1*NAL2
      NGENO=(NGEN*(NGEN+1))/2
      IAUT=IAU(MARK)
      MARC1=.FALSE.
      IF(IAUT.EQ.0) MARC1=.TRUE.
      READ(11) RNAM,(GF2(I),I=1,NAL2),(PHE2(I),I=1,NF2),
     1  ((GVX2(J,I),I=1,NG2),J=1,NF2)
      GO TO 220
  133 CONTINUE
      WRITE(6,11)
      STOP
C
C  READ THETA VALUES (AFTER CARD NO. 9)
  134 IF(IAU(1).LT.9) GO TO 135
      IF(NMARK.EQ.1) GO TO 1342
      WRITE(6,35)
      STOP
 1342 I=0
 1344 I=I+1
      IF(I.LE.MT) GO TO 1346
      WRITE(6,34)
      STOP
 1346 READ(5,12) THV1(I),THV2(I)
      IF(THV1(I).NE.6.) GO TO 1344
C
C  READ GENOTYPE/PHENOTYPE INFORMATION FOR LOCUS 0 (MAIN LOCUS)
  135 NEW1=0
C         GEN = SYMBOLS FOR THE ALLELES AT LOCUS 0
C         PHE1= SYMBOLS FOR THE PHENOTYPES AT LOCUS 0
      IF((NAL1.LE.MNAL1).AND.(NF1.LE.MNF1)) GO TO 136
      I=1
      IF(NAL1.GT.MNAL1) WRITE(6,39) I
      IF(NF1.GT.MNF1) WRITE(6,40) I
      STOP
  136 READ(5,FT1) RNAM1,(GEN(I),I=1,NAL1),(PHE1(I),I=1,NF1)
C         GF1 = GENE FREQUENCIES FOR LOCUS 0
      READ(5,10) (GF1(I),I=1,NAL1)
      R=0.
      DO 137 I=1,NAL1
  137 R=R+GF1(I)
      IF((R.LE.0.9999).OR.(R.GE.1.0001)) WRITE(6,21) R,NEW1,RNAM1
      NG1=(NAL1*(NAL1+1))/2
      DO 150 I=1,NG1
      READ(5,FT2) G1,G2,(PHPROB(J),J=1,NF1)
      IF(KONT1.NE.0) GO TO 139
      DO 138 J=1,NF1
      IF((PHPROB(J).EQ.0.).OR.(PHPROB(J).EQ.1.)) GO TO 138
      WRITE(6,8) NEW1,RNAM1
      MIST=.TRUE.
  138 CONTINUE
  139 DO 140 J=1,NAL1
      IF(G1.NE.GEN(J)) GO TO 140
      I1=J
      GO TO 142
  140 CONTINUE
      I1=1
      MIST=.TRUE.
      WRITE(6,22) NEW1,RNAM1,I,I1
  142 DO 144 J=1,NAL1
      IF(G2.NE.GEN(J)) GO TO 144
      I2=J
      GO TO 146
  144 CONTINUE
      I2=2
      MIST=.TRUE.
      WRITE(6,22) NEW1,RNAM1,I,I2
  146 IA=IAD(I1,I2)
      DO 150 J=1,NF1
  150 GVX1(J,IA)=PHPROB(J)
C
C  READ GENOTYPE/PHENOTYPE INFORMATION FOR MARKER LOCI
      DO 170 MARK=1,NMARK
      NAL2=NAL(MARK)
      NF2=NF(MARK)
      IF((NAL2.LE.MNAL2).AND.(NF2.LE.MNF2)) GO TO 151
      I=2
      IF(NAL2.GT.MNAL2) WRITE(6,39) I
      IF(NF2.GT.MNF2) WRITE(6,40) I
      STOP
  151 READ(5,FT1) RNAM,(GEN(I),I=1,NAL2),(PHE2(I),I=1,NF2)
      READ(5,10) (GF2(I),I=1,NAL2)
      R=0.
      DO 153 I=1,NAL2
  153 R=R+GF2(I)
      IF((R.LE.0.9999).OR.(R.GE.1.0001)) WRITE(6,21) R,MARK,RNAM
      NG2=(NAL2*(NAL2+1))/2
      DO 160 I=1,NG2
      READ(5,FT2) G1,G2,(PHPROB(J),J=1,NF2)
      IF(KONT(MARK).NE.0) GO TO 1545
      DO 154 J=1,NF2
      IF((PHPROB(J).EQ.0.).OR.(PHPROB(J).EQ.1.)) GO TO 154
      WRITE(6,8) MARK,RNAM
      MIST=.TRUE.
  154 CONTINUE
 1545 DO 155 J=1,NAL2
      IF(G1.NE.GEN(J)) GO TO 155
      I1=J
      GO TO 157
  155 CONTINUE
      I1=1
      MIST=.TRUE.
      WRITE(6,22) MARK,RNAM,I,I1
  157 DO 158 J=1,NAL2
      IF(G2.NE.GEN(J)) GO TO 158
      I2=J
      GO TO 159
  158 CONTINUE
      I2=2
      MIST=.TRUE.
      WRITE(6,22) MARK,RNAM,I,I2
  159 IA=IAD(I1,I2)
      DO 160 J=1,NF2
  160 GVX2(J,IA)=PHPROB(J)
      IF(NMARK.EQ.1) GO TO 170
      WRITE(11) RNAM,(GF2(I),I=1,NAL2),(PHE2(I),I=1,NF2),
     1  ((GVX2(J,I),I=1,NG2),J=1,NF2)
  170 CONTINUE
      IPED=0
  180 IF(KAUT.NE.0) GO TO 132
      MARK=0
C
  190 MARK=MARK+1
      IF(MARK.GT.NMARK) GO TO 180
      NAL2=NAL(MARK)
      NF2=NF(MARK)
      IF((NAL1.LE.MNAL1).AND.(NF1.LE.MNF1)) GO TO 192
      I=1
      IF(NAL1.GT.MNAL1) WRITE(6,39) I
      IF(NF1.GT.MNF1) WRITE(6,40) I
      STOP
  192 IF((NAL2.LE.MNAL2).AND.(NF2.LE.MNF2)) GO TO 193
      I=2
      IF(NAL2.GT.MNAL2) WRITE(6,39) I
      IF(NF2.GT.MNF2) WRITE(6,40) I
      STOP
  193 KONT2=KONT(MARK)
      UNK2=UNK(MARK)
      NG2=(NAL2*(NAL2+1))/2
      NGEN=NAL1*NAL2
      NGENO=(NGEN*(NGEN+1))/2
      IAUT=IAU(MARK)
      MARC1=.FALSE.
      IF(IAUT.EQ.0) MARC1=.TRUE.
      IF(NMARK.GT.1) GO TO 191
      GO TO 200
  191 IF(MARK.GT.1) GO TO 196
      GO TO 200
  196 READ(11) RNAM,(GF2(I),I=1,NAL2),(PHE2(I),I=1,NF2),
     1  ((GVX2(J,I),I=1,NG2),J=1,NF2)
      IF(MARC1) GO TO 190
      REWIND 12
      GO TO 230
C
C  READ PEDIGREE INFORMATION
  200 READ(5,6) M,NDI,TEXT
      IF(M.EQ.5000) GO TO 132
      IF(M.EQ.7000) GO TO 180
      IF(M.EQ.9000) GO TO 201
      IF(M.EQ.8000) GO TO 202
      GO TO 208
  201 IF(IAUT.NE.9) STOP
      GO TO 203
  202 IF(IAUT.NE.9) GO TO 110
C
C  PRINT SUM OF LODS
  203 WRITE(6,36)
      I=0
  204 I=I+1
      IF(THV1(I).EQ.6.) GO TO 205
      WRITE(6,15) THV1(I),THV2(I),THVS(I)
      GO TO 204
  205 IF(M.EQ.8000) GO TO 202
      STOP
  208 IF(M.GT.0.AND.M.LE.MNPT) GO TO 210
      WRITE(6,25) M
      STOP
  210 MM=M
      NNDI=NDI
      ASSIGN 246 TO IW1
      ASSIGN 508 TO IW2
      ASSIGN 5400 TO IW3
      IPED=IPED+1
      IF(IPED.GT.1.AND.NMARK.GT.1) GO TO 212
      GO TO 214
  212 WRITE(6,32)
      STOP
  214 WRITE(6,7) IPED,TEXT,M
      IF(NMARK.EQ.1) GO TO 230
      REWIND 11
      REWIND 12
C  READ GENE FREQUENCIES, PHENOTYPE CODES AND GENOT./PHENOT. INFO.
      READ(11) RNAM,(GF2(I),I=1,NAL2),(PHE2(I),I=1,NF2),
     1  ((GVX2(J,I),I=1,NG2),J=1,NF2)
  220 IF(MARC1) GO TO 235
C
C  READ PEDIGREE DATA
  230 WRITE(6,13) NEW1,RNAM1,MARK,RNAM
  235 DO 500 I=1,M
      IF(MARK.GT.1) GO TO 360
      READ(5,FOR) ID(I),NM(I),NC(I),SEX,X1,(PHEPED(J),J=1,NMARK)
      X2=PHEPED(1)
      IF(SEX.EQ.BOY) GO TO 248
      ISEX(I)=.TRUE.
      GO TO IW1,(246,247)
  246 GIRL=SEX
      ASSIGN 247 TO IW1
      GO TO 249
  247 IF(SEX.EQ.GIRL) GO TO 249
      WRITE(6,24) I,ID(I)
      MIST=.TRUE.
      GO TO 249
  248 ISEX(I)=.FALSE.
  249 IF(NMARK.GT.1) WRITE(12) (PHEPED(J),J=1,NMARK)
C
C  PREPARE FENO1
  250 IF(X1.NE.UNK1) GO TO 260
      DO 255 J=1,NG1
  255 FENO1(I,J)=1.
      GO TO 340
  260 IF(KONT1) 270,270,310
  270 DO 280 J=1,NF1
      IF(X1.NE.PHE1(J)) GO TO 280
      K=J
      GO TO 290
  280 CONTINUE
      WRITE(6,9) I,ID(I),NEW1,RNAM1
      MIST=.TRUE.
      K=1
  290 DO 300 J=1,NG1
  300 FENO1(I,J)=GVX1(K,J)
      GO TO 340
  310 CONTINUE
C  The statement below is appropriate for MS-Fortran. The numerical value R1
C  is read from the character variable X1.
      READ(X1,'(F10.0)') R1
      DO 330 J=1,NG1
      R=(R1-GVX1(1,J))/GVX1(2,J)
      IF( ABS(R).GT.DSTD) GO TO 320
      FENO1(I,J)= EXP(-.5*R*R)/(2.5066283*GVX1(2,J))
      GO TO 330
  320 FENO1(I,J)=0.
  330 CONTINUE
C  ADJUST FENO1 FOR SEX-LINKAGE
  340 IF(ISLL.OR.ISEX(I)) GO TO 355
      DO 350 I1=2,NAL1
      I3=I1-1
      DO 350 I2=1,I3
      IA=IAD(I1,I2)
  350 FENO1(I,IA)=0.
  355 IF(MARC1) GO TO 500
      GO TO 370
  360 READ(12) (PHEPED(J),J=1,NMARK)
      X2=PHEPED(MARK)
      IF(MARK.GT.(NEW1+1)) GO TO 370
      X1=PHEPED(NEW1)
      GO TO 250
C
C  PREPARE FENO2
  370 IF(X2.NE.UNK2) GO TO 390
      DO 380 J=1,NG2
  380 FENO2(I,J)=1.
      GO TO 470
  390 IF(KONT2) 400,400,440
  400 DO 410 J=1,NF2
      IF(X2.NE.PHE2(J)) GO TO 410
      K=J
      GO TO 420
  410 CONTINUE
      WRITE(6,9) I,ID(I),MARK,RNAM
      MIST=.TRUE.
      K=1
  420 DO 430 J=1,NG2
  430 FENO2(I,J)=GVX2(K,J)
      GO TO 470
  440 CONTINUE
C  The statement below is appropriate for MS-Fortran. The numerical value of
C  R2 is read from the character variable X2.
      READ(X2,'(F10.0)') R2
      DO 460 J=1,NG2
      R=(R2-GVX2(1,J))/GVX2(2,J)
      IF( ABS(R).GT.DSTD) GO TO 450
      FENO2(I,J)= EXP(-.5*R*R)/(2.5066283*GVX2(2,J))
      GO TO 460
  450 FENO2(I,J)=0.
  460 CONTINUE
C  ADJUST FENO2 FOR SEX-LINKAGE
  470 IF(ISLL.OR.ISEX(I)) GO TO 500
      DO 480 I1=2,NAL2
      I3=I1-1
      DO 480 I2=1,I3
      IA=IAD(I1,I2)
  480 FENO2(I,IA)=0.
  500 CONTINUE
C
C  READ SYMBOLS OF DOUBLED INDIVIDUALS
      GO TO IW3,(5400,507)
 5400 ASSIGN 507 TO IW3
      IF(NDI.EQ.0) GO TO 507
      IF(NDI.GT.0.AND.NDI.LE.MNDI) GO TO 501
      WRITE(6,30) NDI
      STOP
  501 DO 505 I3=1,NDI
      READ(5,FOR) X1,X2
      DO 502 I=1,M
      IF(ID(I).NE.X1) GO TO 502
      LDI(1,I3)=I
      GO TO 503
  502 CONTINUE
      WRITE(6,29) X1
      MIST=.TRUE.
  503 DO 504 I=1,M
      IF(ID(I).NE.X2) GO TO 504
      LDI(2,I3)=I
      GO TO 505
  504 CONTINUE
      WRITE(6,29) X2
      MIST=.TRUE.
  505 CONTINUE
C
  507 IF(MARC1.OR.MIST) GO TO 515
      GO TO IW2,(508,510)
C
  508 CALL SETUP
      ASSIGN 510 TO IW2
      IF(L.GT.1.AND.L.LE.MLIST) GO TO 509
      WRITE(6,31) L
      MIST=.TRUE.
  509 IF(MIST) GO TO 515
      IF(MVEC.LT.1) GO TO 510
      WRITE(6,26) L
      WRITE(6,27) (LIST(I),I=1,L)
      WRITE(6,28) (NS(I),I=1,M)
C
  510 IF((NMARK.EQ.1).AND.(IPED.GT.1)) GO TO 5550
      IF(IHAPL.EQ.0) GO TO 5500
C
C READ HAPLOTYPE VECTOR PHIS (LENGTH  NGEN)
      READ(5,10) (PHIS(I),I=1,NGEN)
      WRITE(6,37) NEW1,RNAM1,MARK,RNAM,(PHIS(I),I=1,NGEN)
      GO TO 5520
C
C COMPUTE HAPLOTYPE VECTOR PHIS
 5500 WRITE(6,15)
      WRITE(6,23) NEW1,RNAM1,(GF1(I),I=1,NAL1)
      WRITE(6,23) MARK,RNAM,(GF2(I),I=1,NAL2)
      I=0
      DO 5510 I1=1,NAL1
      DO 5510 I2=1,NAL2
      I=I+1
 5510 PHIS(I)=GF1(I1)*GF2(I2)
C
C COMPUTE GENOTYPE VECTOR PHI (=  PHIS  FOR MALES IN SEX-LINKED CASE)
 5520 I=0
      DO 5540 I1=1,NGEN
      DO 5530 I2=1,I1
      I=I+1
      R=PHIS(I1)*PHIS(I2)
 5530 PHI(I)=R+R
 5540 PHI(I)=R
 5550 THM=.5
      THF=.5
      ICH=.TRUE.
C
C  FILL A COL. IN HOLD WITH P(X/G) FOR EACH PAIR OF DOUBLED
C  INDIVIDUALS AND ADJUST FOR DOUBLING
      IF(NDI.EQ.0) GO TO 5700
      DO 5610 I=1,NDI
      J=LDI(1,I)
      I5=0
      I6=1
      DO 5610 I3=1,NAL1
      DO 5600 I4=1,NAL2
      DO 5590 I1=1,NAL1
      IA=IAD(I1,I3)
      R1=FENO1(J,IA)
      DO 5590 I2=1,NAL2
      I5=I5+1
      IA=IAD(I2,I4)
      R2=R1*FENO2(J,IA)
      IF(R2.EQ.0.) GO TO 5575
      IF(ISLL.OR.ISEX(J)) GO TO 5560
      R2=R2/PHIS(I6)
      GO TO 5570
 5560 R2=R2/PHI(I5)
 5570 HOLD(I,I5)= SQRT(R2)
      GO TO 5580
 5575 HOLD(I,I5)=0.
 5580 IF(I1.EQ.I3.AND.I2.EQ.I4) GO TO 5600
 5590 CONTINUE
 5600 I6=I6+1
 5610 CONTINUE
 5700 CONTINUE
C
C  DETERMINE IF LIKELIHOOD TO BE RESCALED (IRESC=.FALSE.)
C  OR NOT (IRESC=.TRUE.)
      IF(ISC) 513,511,512
  511 IF(KONT1.EQ.0.AND.KONT2.EQ.0) GO TO 513
  512 IRESC=.FALSE.
      GO TO 514
  513 IRESC=.TRUE.
C
  514 CALL PLIK
      ICH=.FALSE.
      IF(RLIK.NE.0.) GO TO 530
      IF(NDI.NE.0) WRITE(6,33)
  515 IF(IAUT.NE.1) GO TO 190
  520 READ(5,12) THM,THF
      IF(THM.EQ.6.) GO TO 190
      GO TO 520
  530 RLIK5=RLIK
      RLIK2=0.
      IF(IAUT.NE.-1) GO TO 540
      WRITE(6,14)
      WRITE(6,15) THM,THF,RLIK,RLIK2
      GO TO 190
C
  540 WRITE(6,14)
      WRITE(6,15) THM,THF,RLIK,RLIK2
C
C  Select output options
      GO TO (545,570,580,570,580,630,640,645,541),IAUT
  541 I=0
  542 I=I+1
      THM=THV1(I)
      THF=THV2(I)
      GO TO 550
  545 READ(5,12) THM,THF
  550 IF(THM.EQ.6.) GO TO 190
      CALL PLIK
      IF(RLIK.NE.0) GO TO 560
      RLIK=-1.E08
      RLIK2=-1.E08
      GO TO 565
  560 RLIK2=RLIK-RLIK5
  565 WRITE(6,15) THM,THF,RLIK,RLIK2
      IF(IAUT.EQ.1) GO TO 545
      THVS(I)=THVS(I)+RLIK2
      GO TO 542
  570 I1=5
      R=.1
      GO TO 590
  580 I1=10
      R=.05
C
  590 TH=THM
      DO 620 I=1,I1
      IW=0
      TH =TH -R
      IF(TH.LT.0.04) TH=0.
      THF=TH
      THM=TH
      GO TO 598
  594 THM=0.5
      GO TO 598
  596 THM=TH
      THF=0.5
  598 CALL PLIK
      IF(RLIK.NE.0) GO TO 610
      RLIK=-1.E08
      RLIK2=-1.E08
      GO TO 615
  610 RLIK2=RLIK-RLIK5
  615 WRITE(6,15) THM,THF,RLIK,RLIK2
      IF(IAUT.LT.4) GO TO 620
      IW=IW+1
      GO TO (594,596,620),IW
  620 CONTINUE
      GO TO 190
  630 I1=3
      GO TO 650
  640 I1=4
      GO TO 650
  645 I1=5
C
  650 WRITE(6,16) BOY
      I2=I1-2
      DO 690 I=1,I1
      I3=I1+1-I
      THM=THV(I2,I3)
      DO 680 J=1,I1
      IF(I.EQ.1.AND.J.EQ.I1) GO TO 675
      THF=THV(I2,J)
      CALL PLIK
      IF(RLIK.NE.0) GO TO 670
      PHPROB(J)=-1.E08
      GO TO 680
  670 PHPROB(J)=RLIK-RLIK5
      GO TO 680
  675 PHPROB(J)=0.
  680 CONTINUE
  690 WRITE(6,17) THM,(PHPROB(J),J=1,I1)
      WRITE(6,18) (THV(I2,I),I=1,I1)
      WRITE(6,19) GIRL
      GO TO 190
      END
      SUBROUTINE QLIK
C
C  TO COMPUTE RLIK = LIKELIHOOD WITH GIVEN VALUES THM
C  AND THF OF THE MALE AND FEMALE RECOMBINATION FRACTION.
C
C  TO ALLOW FOR MUTATION, THIS VERSION OF QLIK (WITHOUT MUTATION)
C  HAS TO BE REPLACED BY A VERSION OF QLIK WITH MUTATION. TO
C  CARRY OUT THIS TRANSFORMATION, DO THE FOLLOWING.
C    LINES WITH -M-  IN COL. 73     -- REMOVE THE C FROM COL. 1
C    LINES WITH -MC- IN COL. 73-74  -- REMOVE THE C FROM COL. 1 AND
C      PUT A -C- IN COL. 1 OF THE LINE BELOW IT.
C  MUTATION IS ASSUMED TO OCCUR WITH RATE RMUT FROM ANY OF
C  THE ALLELES A(I),I=2,...,NAL1 TOWARDS A(1), I.E., THE FIRST
C  ALLELE, AT THE FIRST LOCUS UNDER CONSIDERATION (MAIN LOCUS).
C
      COMMON /B/ LIST(50),NM(200),NC(200),NS(200),M,L,NOPAR,ID(200),
     1  NDI,LDI(2,3),LGC(3),ISC,MIST,ISEX(200),
     1  FENO1(200,10),FENO2(200,21),PHI(300),RMUT,
     1  HOLD(3,300),RLIK,PHIS(36),WVK(4),THM,THF,IAD(36,36),
     2  IV(4),JV(4),IVA1(4),IVA2(16),IVA3(4),IVF(4),NGEN,NGENO,MNP,
     3  MNPT,KONT1,KONT2,NAL1,NAL2,MVEC,ICH,ISLL,IRESC
      COMMON GENO(25,300),LGENO(200),IGENO(25)
      DIMENSION WVP(4),WVPS(2)
C     DIMENSION RMUE(2),RMUF(2)                                         M
      CHARACTER*4 ID
      LOGICAL ISLL,ISEX,ICH,IGENO,MIST,IRESC,IZERO
C
    1 FORMAT(/32H ERROR IN QLIK. DIMENSION MNP IN /
     1 3X,37HGENO(MNP,.) AND IGENO(MNP) TOO SMALL.  /
     2 3X,22HDETECTED AT PERSON NR.,I4,3X,A4)
    2 FORMAT(/27H MODIFIED VECTOR FOR PERSON,I4,3X,A4/(1X,5E14.6))
    3 FORMAT(//8X,31HLIKELIHOOD IS ZERO AT R = .5/.5
     1 /8X,38HINCOMPATIBILITY DETECTED AT PERSON NR.,I4,3X,A4)
C
      RLIKS=0.
      WVP(4)=THM*THF*.25
      C     =THM*.25-WVP(4)
      B     =THF*.25-WVP(4)
      WVP(1)=.25-B-C-WVP(4)
      WVPS(2)=.5*THF
      WVPS(1)=.5-WVPS(2)
C     RMUT3=(1.-RMUT)**2                                                M 
C     WVP(1)=WVP(1)*RMUT3                                               M 
C     WVP(4)=WVP(4)*RMUT3                                               M 
C     B=B*RMUT3                                                         M 
C     C=C*RMUT3                                                         M 
C     WVPS(1)=WVPS(1)*(1.-RMUT)                                         M 
C     WVPS(2)=WVPS(2)*(1.-RMUT)                                         M 
C     RMUT1=RMUT*(1.-RMUT)/4.                                           M 
C     RMUT2=RMUT/2.                                                     M 
      IF(KONT1) 80,50,80
   50 IF(KONT2) 70,60,70
C  BOTH PHENOTYPES ARE DISCRETE
   60 ASSIGN 150 TO ISW1
      ASSIGN 355 TO ISW2
      ASSIGN 560 TO ISW5
      ASSIGN 790 TO ISW8
      ASSIGN 1080 TO ISW13
C     ASSIGN 637 TO ISW16                                               MC
      ASSIGN 642 TO ISW16
C     ASSIGN 1160 TO ISW17                                              MC
      ASSIGN 642 TO ISW17
      GO TO 95
C  1 IS DISCRETE, 2 IS CONTINUOUS
   70 ASSIGN 160 TO ISW1
      ASSIGN 360 TO ISW2
      ASSIGN 570 TO ISW5
      ASSIGN 800 TO ISW8
      ASSIGN 1090 TO ISW13
      ASSIGN 634 TO ISW16
      ASSIGN 1152 TO ISW17
      GO TO 95
   80 IF(KONT2) 90,85,90
C  1 IS CONTINUOUS, 2 IS DISCRETE
   85 ASSIGN 170 TO ISW1
      ASSIGN 370 TO ISW2
      ASSIGN 580 TO ISW5
      ASSIGN 810 TO ISW8
      ASSIGN 1100 TO ISW13
      ASSIGN 635 TO ISW16
      ASSIGN 1154 TO ISW17
      GO TO 95
C  BOTH PHENOTYPES ARE CONTINUOUS
   90 ASSIGN 180 TO ISW1
      ASSIGN 380 TO ISW2
      ASSIGN 590 TO ISW5
      ASSIGN 820 TO ISW8
      ASSIGN 1110 TO ISW13
      ASSIGN 636 TO ISW16
      ASSIGN 1156 TO ISW17
C
   95 RLIK=1.
C
C  GO THRU ALL HEADSIBS
      DO 900 IL=1,L
      IH=LIST(IL)
      NPI=NM(IH)
      NPJ=NC(IH)
      IF(NPI.EQ.0) GO TO 750
      IF(ISEX(NPI)) GO TO 100
C  PARENT I IS MALE
      WVP(2)=B
      WVP(3)=C
C     RMUE(2)=RMUT1*THM                                                 M 
C     RMUE(1)=RMUT1-RMUE(2)                                             M 
C     RMUF(2)=RMUT1*THF                                                 M 
C     RMUF(1)=RMUT1-RMUF(2)                                             M 
      GO TO 110
C  PARENT J IS MALE
  100 WVP(2)=C
      WVP(3)=B
C     RMUE(2)=RMUT1*THF                                                 M 
C     RMUE(1)=RMUT1-RMUE(2)                                             M 
C     RMUF(2)=RMUT1*THM                                                 M 
C     RMUF(1)=RMUT1-RMUF(2)                                             M 
  110 IF(ISLL) GO TO 111
      GO TO 112
C     NO SEX-LINKAGE
  111 ASSIGN 475 TO ISW9
      ASSIGN 300 TO ISW12
      ASSIGN 475 TO ISW14
      ASSIGN 1175 TO ISW15
      GO TO 120
  112 IF(ISEX(NPI)) GO TO 113
      GO TO 114
C     SEX-LINKAGE, PARENT NPI IS FEMALE
  113 ASSIGN 395 TO ISW9
      ASSIGN 290 TO ISW12
      ASSIGN 475 TO ISW14
      ASSIGN 1180 TO ISW15
      GO TO 120
C     SEX-LINKAGE, PARENT NPJ IS FEMALE
  114 ASSIGN 396 TO ISW9
      ASSIGN 300 TO ISW12
      ASSIGN 430 TO ISW14
      ASSIGN 1175 TO ISW15
C  CHECK WHETHER PARENT I (TO BE MODIFIED) ALREADY HAS A COL.IN GENO
  120 NCOLI=LGENO(NPI)
      NCOLJ=LGENO(NPJ)
      IF(NCOLI) 250,130,250
C  SET UP NEW COL. IN GENO FOR PARENT I
  130 DO 140 I=1,MNP
      IF(IGENO(I)) GO TO 140
      NCOLI=I
      LGENO(NPI)=I
      IGENO(I)=.TRUE.
      GO TO 145
  140 CONTINUE
C
      WRITE(6,1) NPI,ID(NPI)
      RLIK=0.
      RETURN
C  FILL THIS COL. WITH G-VECTOR, P(PHENOTYPE/GENOTYPE)
  145 I5=0
      DO 230 I3=1,NAL1
      DO 220 I4=1,NAL2
      DO 210 I1=1,NAL1
      IA=IAD(I1,I3)
      R1=FENO1(NPI,IA)
      DO 210 I2=1,NAL2
      I5=I5+1
      IF(R1.EQ.0.) GO TO 190
      IA=IAD(I2,I4)
      R2=FENO2(NPI,IA)
      IF(R2.EQ.0.) GO TO 190
      GO TO ISW1,(150,160,170,180)
  150 GENO(NCOLI,I5)=1.
      GO TO 200
  160 GENO(NCOLI,I5)=R2
      GO TO 200
  170 GENO(NCOLI,I5)=R1
      GO TO 200
  180 GENO(NCOLI,I5)=R2*R1
      GO TO 200
  190 GENO(NCOLI,I5)=0.
  200 IF(I1.EQ.I3.AND.I2.EQ.I4) GO TO 220
  210 CONTINUE
  220 CONTINUE
  230 CONTINUE
C
C  START LOOP FOR PARENT I (ID NUMBER NPI)
  250 I5=0
      IZERO=.FALSE.
      DO 710 I3=1,NAL1
      DO 700 I4=1,NAL2
      DO 695 I1=1,NAL1
      DO 690 I2=1,NAL2
      I5=I5+1
      IF(I1.EQ.I3.AND.I2.EQ.I4) GO TO 260
      ASSIGN 500 TO ISW3
      ASSIGN 690 TO ISW7
      ASSIGN 398 TO ISW10
      GO TO 270
  260 ASSIGN 490 TO ISW3
      ASSIGN 700 TO ISW7
      ASSIGN 392 TO ISW10
  270 RI=GENO(NCOLI,I5)
      IF(RI) 280,680,280
C  COMPUTE ADDRESSES NEEDED FOR KIDS
  280 IV(1)=(I1-1)*NAL2+I2
      IV(2)=(I3-1)*NAL2+I4
      IV(3)=IV(1)-I2+I4
      IV(4)=IV(2)-I4+I2
      GO TO ISW12,(300,290)
  290 DO 295 I=1,4
  295 IVF(I)=(IV(I)*(IV(I)+1))/2
      IVA3(1)=IAD(I1,I1)
      IVA3(2)=IAD(I2,I2)
      IVA3(3)=IAD(I3,I3)
      IVA3(4)=IAD(I4,I4)
C
C  START LOOP FOR PARENT J (ID NUMBER NPJ)
  300 SJ=0.
      J5=0
      J6=1
      DO 670 J3=1,NAL1
      DO 660 J4=1,NAL2
      DO 655 J1=1,NAL1
      IA=IAD(J3,J1)
      R11=FENO1(NPJ,IA)
      DO 650 J2=1,NAL2
      J5=J5+1
      IF(J1.EQ.J3.AND.J2.EQ.J4) GO TO 310
      ASSIGN 500 TO ISW4
      ASSIGN 650 TO ISW6
      ASSIGN 398 TO ISW11
      GO TO 320
  310 ASSIGN 620 TO ISW4
      ASSIGN 660 TO ISW6
      ASSIGN 394 TO ISW11
  320 IF(NCOLJ) 330,350,330
C  PARENT J HAS A COL. IN GENO (IS RARE)
  330 RJ=GENO(NCOLJ,J5)
      IF(RJ) 382,645,382
C  PARENT J HAS NO COL. IN GENO
  350 IF(R11.EQ.0.) GO TO 645
      IA=IAD(J2,J4)
      R2=FENO2(NPJ,IA)
      IF(R2.EQ.0.) GO TO 645
      GO TO ISW2,(355,360,370,380)
  355 RJ=1.
      GO TO 382
  360 RJ=R2
      GO TO 382
  370 RJ=R11
      GO TO 382
  380 RJ=R11*R2
C
C  COMPUTE ADDRESSES NEEDED FOR KIDS
  382 GO TO ISW10,(392,398)
  392 GO TO ISW11,(394,398)
C
C     BOTH PARENTS ARE DOUBLY HOMOZYGOUS
  394 IA7=(J1-1)*NAL2+J2
      IA=IV(1)
      IA3=IAD(IA7,IA)
      IA1=IAD(I1,J1)
      IA2=IAD(I2,J2)
C     JV(1)=(J1-1)*NAL2+J2                                              M 
C     JV(2)=(J3-1)*NAL2+J4                                              M 
C     JV(3)=JV(1)-J2+J4                                                 M 
C     JV(4)=JV(2)-J4+J2                                                 M 
C
      GO TO ISW9,(475,395,396)
C     NPI = FEMALE
  395 IA4=IVF(1)
      IA5=IVA3(1)
      IA6=IVA3(2)
      GO TO 475
C     NPJ = FEMALE
  396 IA4=(IA7*(IA7+1))/2
      IA5=(J1*(J1+1))/2
      IA6=(J2*(J2+1))/2
C     IVA3(2)=IAD(J2,J2)                                                M 
C     IVA3(4)=IAD(J4,J4)                                                M 
      GO TO 475
C
C     NOT BOTH PARENTS ARE DOUBLY HOMOZYGOUS
  398 JV(1)=(J1-1)*NAL2+J2
      JV(2)=(J3-1)*NAL2+J4
      JV(3)=JV(1)-J2+J4
      JV(4)=JV(2)-J4+J2
      IVA1(1)=IAD(I1,J1)
      IVA1(2)=IAD(I1,J3)
      IVA1(3)=IAD(I3,J1)
      IVA1(4)=IAD(I3,J3)
      IVA2(1)=IAD(I2,J2)
      IVA2(2)=IAD(I2,J4)
      IVA2(3)=IAD(I4,J2)
      IVA2(4)=IAD(I4,J4)
      DO 400 I=13,16
  400 IVA2(I)=IVA2(17-I)
      IVA2(5)=IVA2(2)
      IVA2(6)=IVA2(1)
      IVA2(7)=IVA2(4)
      IVA2(8)=IVA2(3)
      DO 410 I=9,12
  410 IVA2(I)=IVA2(17-I)
C
      GO TO ISW14,(475,430)
C     NPJ = FEMALE
  430 DO 440 I=1,4
  440 IVF(I)=(JV(I)*(JV(I)+1))/2
      IVA3(1)=IAD(J1,J1)
      IVA3(2)=IAD(J2,J2)
      IVA3(3)=IAD(J3,J3)
      IVA3(4)=IAD(J4,J4)
  475 KNEW=IH
C
C  START LOOP FOR SIB, ID NUMBER KNEW (FEMALE ONLY IF ISL=1)
  480 NCOLK=LGENO(KNEW)
      SK=0.
      GO TO ISW3,(490,500)
  490 GO TO ISW4,(620,500)
C  NOT BOTH PARENTS ARE DOUBLY HOMOZYGOUS
  500 IF(NCOLK) 510,540,510
C     SIB HAS COL. IN GENO
  510 IF(ISLL.OR.ISEX(KNEW)) GO TO 511
      GO TO 1050
  511 K5=0
      DO 516 K1=1,3,2
      K3=K1+1
      DO 514 K2=1,3,2
      K4=K2+1
      K5=K5+1
      R=0.
      DO 512 K6=K1,K3
      DO 512 K7=K2,K4
      IA=IV(K6)
      IB=JV(K7)
      IA=IAD(IA,IB)
  512 R=R+GENO(NCOLK,IA)
      IF(R.EQ.0.) GO TO 514
      SK=SK+R*WVP(K5)
  514 CONTINUE
  516 CONTINUE
C 518 DO 520 K7=1,2                                                     M 
C     IA=(K7-1)*2+1                                                     M 
C     K1=IV(IA)                                                         M 
C     K2=IV(IA+1)                                                       M 
C     K3=IAD(K1,J2)                                                     M 
C     K4=IAD(K1,J4)                                                     M 
C     K5=IAD(K2,J2)                                                     M 
C     K6=IAD(K2,J4)                                                     M 
C     R=GENO(NCOLK,K5)+GENO(NCOLK,K6)+GENO(NCOLK,K3)+GENO(NCOLK,K4)     M 
C 520 SK=SK+R*RMUE(K7)                                                  M 
C     DO 530 K7=1,2                                                     M 
C     IA=(K7-1)*2+1                                                     M 
C     K1=JV(IA)                                                         M 
C     K2=JV(IA+1)                                                       M 
C     K3=IAD(K1,I2)                                                     M 
C     K4=IAD(K2,I2)                                                     M 
C     K5=IAD(K1,I4)                                                     M 
C     K6=IAD(K2,I4)                                                     M 
C     R=GENO(NCOLK,K5)+GENO(NCOLK,K6)+GENO(NCOLK,K3)+GENO(NCOLK,K4)     M 
C 530 SK=SK+R*RMUF(K7)                                                  M 
      GO TO 639
C     SIB HAS NO COL. IN GENO
  540 IF(ISLL.OR.ISEX(KNEW)) GO TO 545
      GO TO 1070
  545 DO 550 K1=1,4
  550 WVK(K1)=0.
      DO 610 K1=1,4
      IA=IVA1(K1)
      R1=FENO1(KNEW,IA)
      IF(R1.EQ.0.) GO TO 610
      K5=0
      DO 600 K2=K1,16,4
      K5=K5+1
      IA=IVA2(K2)
      R2=FENO2(KNEW,IA)
      IF(R2.EQ.0.) GO TO 600
      GO TO ISW5,(560,570,580,590)
  560 WVK(K5)=WVK(K5)+1.
      GO TO 600
  570 WVK(K5)=WVK(K5)+R2
      GO TO 600
  580 WVK(K5)=WVK(K5)+R1
      GO TO 600
  590 WVK(K5)=WVK(K5)+R1*R2
  600 CONTINUE
  610 CONTINUE
      DO 615 K1=1,4
      IF(WVK(K1).EQ.0.) GO TO 615
      SK=SK+WVK(K1)*WVP(K1)
  615 CONTINUE
C 616 IA=IAD(1,I1)                                                      M 
C     R1=FENO1(KNEW,IA)                                                 M 
C     K1=IAD(I2,J2)                                                     M 
C     K2=IAD(I2,J4)                                                     M 
C     K3=IAD(I4,J2)                                                     M 
C     K4=IAD(I4,J4)                                                     M 
C     SK=SK+R1*((FENO2(KNEW,K1)+FENO2(KNEW,K2))*RMUE(1)                 M 
C    1  +(FENO2(KNEW,K3)+FENO2(KNEW,K4))*RMUE(2))                       M 
C     IA=IAD(1,I3)                                                      M 
C     R1=FENO1(KNEW,IA)                                                 M 
C     SK=SK+R1*((FENO2(KNEW,K3)+FENO2(KNEW,K4))*RMUE(1)                 M 
C    1  +(FENO2(KNEW,K1)+FENO2(KNEW,K2))*RMUE(2))                       M 
C     IA=IAD(1,J1)                                                      M 
C     SK=SK+FENO1(KNEW,IA)*((FENO2(KNEW,K1)+FENO2(KNEW,K3))*RMUF(1)     M 
C    1  +(FENO2(KNEW,K2)+FENO2(KNEW,K4))*RMUF(2))                       M 
C     IA=IAD(1,J3)                                                      M 
C     SK=SK+FENO1(KNEW,IA)*((FENO2(KNEW,K2)+FENO2(KNEW,K4))*RMUF(1)     M 
C    1  +(FENO2(KNEW,K1)+FENO2(KNEW,K3))*RMUF(2))                       M 
      GO TO 639
C  BOTH PARENTS ARE DOUBLY HOMOZYGOUS
  620 IF(NCOLK) 630,632,630
  630 IF(ISLL.OR.ISEX(KNEW)) GO TO 631
      GO TO 1140
  631 SK=GENO(NCOLK,IA3)
C     SK=SK*RMUT3                                                       M 
C     GO TO 518                                                         MC
      GO TO 639
  632 IF(ISLL.OR.ISEX(KNEW)) GO TO 633
      GO TO 1150
  633 R1=FENO1(KNEW,IA1)
C     IF(R1.EQ.0.) GO TO 616                                            MC
      IF(R1.EQ.0.) GO TO 645
      R2=FENO2(KNEW,IA2)
C     IF(R2.EQ.0.) GO TO 616                                            MC
      IF(R2.EQ.0.) GO TO 645
C     SK=1.                                                             M 
C     GO TO ISW16,(637,634,635,636)                                     MC
      GO TO ISW16,(642,634,635,636)
  634 SK=R2
C     GO TO 637                                                         MC
      GO TO 641
  635 SK=R1
C     GO TO 637                                                         MC
      GO TO 641
  636 SK=R1*R2
C 637 SK=SK*RMUT3                                                       M 
C     GO TO 616                                                         MC
      GO TO 641
  639 IF(SK.EQ.0.) GO TO 645
  641 RJ=RJ*SK
  642 KOLD=KNEW
      KNEW=NS(KOLD)
      IF(KNEW-IH) 480,1170,480
C  END LOOP OF SIB
C
C  LOOP FOR MALE SIB IN SEX-LINKED CASE
 1050 J=0
      DO 1060 I=1,3,2
      J=J+1
      IA=IVF(I)
      IB=IVF(I+1)
      R=GENO(NCOLK,IA)+GENO(NCOLK,IB)
      IF(R.EQ.0.) GO TO 1060
      SK=SK+R*WVPS(J)
 1060 CONTINUE
C1062 K1=IVA3(2)                                                        M 
C     K2=IVA3(4)                                                        M 
C     SK=SK+RMUT2*(GENO(NCOLK,K1)+GENO(NCOLK,K2))                       M 
      GO TO 639
C
 1070 DO 1130 K1=1,2
      K3=6*K1-6
      K4=K1+K1-K3
      R=0.
      DO 1120 K2=1,3,2
      K3=K3+K4
      IA=IVA3(K2)
      R1=FENO1(KNEW,IA)
      IF(R1.EQ.0.) GO TO 1120
      IA=IVA3(K3)
      R2=FENO2(KNEW,IA)
      IF(R2.EQ.0.) GO TO 1120
      GO TO ISW13,(1080,1090,1100,1110)
 1080 R=R+1.
      GO TO 1120
 1090 R=R+R2
      GO TO 1120
 1100 R=R+R1
      GO TO 1120
 1110 R=R+R1*R2
 1120 CONTINUE
      IF(R.EQ.0.) GO TO 1130
      SK=SK+R*WVPS(K1)
 1130 CONTINUE
C1132 K1=IVA3(2)                                                        M 
C     K2=IVA3(4)                                                        M 
C     SK=SK+FENO1(KNEW,1)*RMUT2*(FENO2(KNEW,K1)+FENO2(KNEW,K2))         M 
      GO TO 639
C
 1140 SK=GENO(NCOLK,IA4)
C     SK=SK*RMUT3                                                       M 
C     GO TO 1062                                                        MC
      GO TO 639
C
 1150 R1=FENO1(KNEW,IA5)
C     IF(R1.EQ.0.) GO TO 1132                                           MC
      IF(R1.EQ.0.) GO TO 645
      R2=FENO2(KNEW,IA6)
C     IF(R2.EQ.0.) GO TO 1132                                           MC
      IF(R2.EQ.0.) GO TO 645
C     SK=1.                                                             M 
C     GO TO ISW17,(1160,1152,1154,1156)                                 MC
      GO TO ISW17,(642,1152,1154,1156)
 1152 SK=R2
C     GO TO 1160                                                        MC
      GO TO 641
 1154 SK=R1
C     GO TO 1160                                                        MC
      GO TO 641
 1156 SK=R1*R2
C1160 SK=SK*RMUT3                                                       M 
C     GO TO 1132                                                        MC
      GO TO 641
C  END LOOP FOR MALE SIB IN SEX-LINKED CASE
C
C  ADD MODIFIED COMPONENT RJ TO SUM FOR PARENT J
 1170 GO TO ISW15,(1175,1180)
 1175 SJ=SJ+RJ*PHI(J5)
      GO TO 645
 1180 SJ=SJ+RJ*PHIS(J6)
C
  645 GO TO ISW6,(660,650)
  650 CONTINUE
  655 CONTINUE
  660 J6=J6+1
  670 CONTINUE
C
C  MODIFY COMPONENT RI
      IF(SJ) 675,672,675
  672 GENO(NCOLI,I5)=0.
      GO TO 680
  675 GENO(NCOLI,I5)=RI*SJ
      IZERO=.TRUE.
  680 GO TO ISW7,(700,690)
  690 CONTINUE
  695 CONTINUE
  700 CONTINUE
  710 CONTINUE
C
C  PRINT MODIFIED VECTORS
      IF(MVEC.NE.1) GO TO 711
      WRITE(6,2) NPI,ID(NPI),(GENO(NCOLI,I),I=1,NGENO)
C
C  PRINT MESSAGE IF MODIFIED VECTOR IS ZERO THRUOUT (IZERO=.FALSE.)
C  AND IF REC.FRACTIONS = .5, NDI = 0
  711 IF(IZERO) GO TO 712
      IF(ICH.AND.NDI.EQ.0) WRITE(6,3) NPI,ID(NPI)
      GO TO 880
C
C  RESCALE MODIFIED VECTORS SUCH THAT LARGEST ELEMENT EQUALS 1
  712 IF(IRESC) GO TO 715
      VEMA=-1.
      DO 713 IA=1,NGENO
      IF(GENO(NCOLI,IA).GT.VEMA) VEMA=GENO(NCOLI,IA)
  713 CONTINUE
      DO 714 IA=1,NGENO
  714 GENO(NCOLI,IA)=GENO(NCOLI,IA)/VEMA
      RLIKS=RLIKS+ALOG10(VEMA)
C
C  RELEASE COLUMNS IN GENO FOR SIBS AND PARENT J
  715 KNEW=IH
  720 IA=LGENO(KNEW)
      IF(IA) 725,730,725
  725 LGENO(KNEW)=0
      IGENO(IA)=.FALSE.
  730 KOLD=KNEW
      KNEW=NS(KOLD)
      IF(KNEW.NE.IH) GO TO 720
      IF(NCOLJ.EQ.0) GO TO 900
      LGENO(NPJ)=0
      IGENO(NCOLJ)=.FALSE.
      GO TO 900
C
C  DO LAST SUMMATION AND COMBINE INDEPENDENT PARTS OF PEDIGREE
  750 NCOLH=LGENO(IH)
      SLIK=0.
      IF(NCOLH) 760,780,760
C        WITH COLUMN IN GENO
  760 LGENO(IH)=0
      IGENO(NCOLH)=.FALSE.
      IF(ISLL.OR.ISEX(IH)) GO TO 765
      GO TO 772
  765 DO 770 I=1,NGENO
      R=GENO(NCOLH,I)
      IF(R.EQ.0.) GO TO 770
      SLIK=SLIK+PHI(I)*R
  770 CONTINUE
      GO TO 870
  772 J=0
      DO 774 I=1,NGEN
      J=J+I
      R=GENO(NCOLH,J)
      IF(R.EQ.0.) GO TO 774
      SLIK=SLIK+PHIS(I)*R
  774 CONTINUE
      GO TO 870
C        WITHOUT COL. IN GENO (RARE)
  780 I5=0
      I6=1
      DO 860 I3=1,NAL1
      DO 850 I4=1,NAL2
      DO 840 I1=1,NAL1
      IA=IAD(I1,I3)
      R1=FENO1(IH,IA)
      DO 840 I2=1,NAL2
      I5=I5+1
      IF(R1.EQ.0.) GO TO 830
      IA=IAD(I2,I4)
      R2=FENO2(IH,IA)
      IF(R2.EQ.0.) GO TO 830
      GO TO ISW8,(790,800,810,820)
  790 IF(ISLL.OR.ISEX(IH)) GO TO 795
      SLIK=SLIK+PHIS(I6)
      GO TO 830
  795 SLIK=SLIK+PHI(I5)
      GO TO 830
  800 IF(ISLL.OR.ISEX(IH)) GO TO 805
      SLIK=SLIK+PHIS(I6)*R2
      GO TO 830
  805 SLIK=SLIK+PHI(I5)*R2
      GO TO 830
  810 IF(ISLL.OR.ISEX(IH)) GO TO 815
      SLIK=SLIK+PHIS(I6)*R1
      GO TO 830
  815 SLIK=SLIK+PHI(I5)*R1
      GO TO 830
  820 IF(ISLL.OR.ISEX(IH)) GO TO 825
      SLIK=SLIK+PHIS(I6)*R1*R2
      GO TO 830
  825 SLIK=SLIK+PHI(I5)*R1*R2
  830 IF(I1.EQ.I3.AND.I2.EQ.I4) GO TO 850
  840 CONTINUE
  850 I6=I6+1
  860 CONTINUE
  870 IF(SLIK) 890,880,890
  880 RLIK=0.
      GO TO 920
  890 RLIK=RLIK*SLIK
  900 CONTINUE
      IF(IRESC) GO TO 920
      RLIK=ALOG10(RLIK)+RLIKS
  920 RETURN
      END
C     ----------------------------------------
      SUBROUTINE PLIK
C  TO REPEATEDLY CALL QLIK FOR COMPLEX PEDIGREES (NDI.GT.0)
      COMMON /B/ LIST(50),NM(200),NC(200),NS(200),M,L,NOPAR,ID(200),
     1  NDI,LDI(2,3),LGC(3),ISC,MIST,ISEX(200),
     1  FENO1(200,10),FENO2(200,21),PHI(300),RMUT,
     1  HOLD(3,300),RLIK,PHIS(36),WVK(4),THM,THF,IAD(36,36),
     2  IV(4),JV(4),IVA1(4),IVA2(16),IVA3(4),IVF(4),NGEN,NGENO,MNP,
     3  MNPT,KONT1,KONT2,NAL1,NAL2,MVEC,ICH,ISLL,IRESC
      COMMON GENO(25,300),LGENO(200),IGENO(25)
      LOGICAL ISLL,ISEX,ICH,IGENO,MIST,IRESC,IZERO
      RLS=0.
      RFAC=0.
      IW=0
      IF(NDI.NE.0) GO TO 150
      DO 80 I=1,MNP
   80 IGENO(I)=.FALSE.
      DO 90 I=1,M
   90 LGENO(I)=0
      CALL QLIK
      IF(IRESC) GO TO 100
      RETURN
  100 IF(RLIK.EQ.0.0) RETURN
      RLIK=ALOG10(RLIK)
      RETURN
C
  150 DO 160 I=1,NDI
  160 LGC(I)=1
      LGC(1)=0
C
C  GO THRU ALL GENOTYPE COMBINATIONS OF THE NDI INDIVIDUALS
  200 I=1
  210 LGC(I)=LGC(I)+1
      IF(LGC(I).LE.NGENO) GO TO 220
      LGC(I)=1
      I=I+1
      IF(I.GT.NDI) GO TO 300
      GO TO 210
C  FOR GIVEN COMBINATION, COMPUTE LIKELIHOOD
  220 DO 230 I=1,NDI
      J=LGC(I)
      IF(HOLD(I,J).EQ.0.0) GO TO 200
  230 CONTINUE
      DO 232 I=1,MNP
  232 IGENO(I)=.FALSE.
      DO 234 I=1,M
  234 LGENO(I)=0
      I1=NDI+NDI
      DO 240 I=1,I1
      IGENO(I)=.TRUE.
      DO 240 J=1,NGENO
  240 GENO(I,J)=0.
      DO 250 I=1,NDI
      I2=I+I
      I1=I2-1
      J=LGC(I)
      I3=LDI(1,I)
      I4=LDI(2,I)
      LGENO(I3)=I1
      LGENO(I4)=I2
      GENO(I1,J)=HOLD(I,J)
  250 GENO(I2,J)=HOLD(I,J)
      CALL QLIK
      IF(RLIK.EQ.0.0) GO TO 200
      IF(IRESC) GO TO 270
C     WITH RESCALING
      IF(IW.EQ.1) GO TO 260
      RFAC=RLIK
      RLS=1.
      IW=1
      GO TO 200
  260 R=RLIK-RFAC
      RLS=RLS+10.**R
      GO TO 200
C     NO RESCALING
  270 RLS=RLS+RLIK
      GO TO 200
  300 IF(IRESC) GO TO 310
      IF(RFAC.EQ.0.0) RETURN
      RLIK=RFAC+ALOG10(RLS)
      RETURN
  310 IF(RLS.EQ.0.0) RETURN
      RLIK=ALOG10(RLS)
      RETURN
      END
      SUBROUTINE SETUP
C  TO REPLACE THE CODE OF EACH PERSON BY A NUMBER I=1,M AND TO
C  ESTABLISH THE ORDER IN WHICH THE PERSONS HAVE TO BE PROCESSED.
C        WRITTEN BY STANLEY ALBERT.
      COMMON /B/ LIST(50),AN2(200),AN1(200),NS(200),M,L,NOPAR,ID(200),
     1  NDI,LDI(2,3),LGC(3),ISC,MIST,ISEX(200),
     1  FENO1(200,10),FENO2(200,21),PHI(300),RMUT,
     1  HOLD(3,300),RLIK,PHIS(36),WVK(4),THM,THF,IAD(36,36),
     2  IV(4),JV(4),IVA1(4),IVA2(16),IVA3(4),IVF(4),NGEN,NGENO,MNP,
     3  MNPT,KONT1,KONT2,NAL1,NAL2,MVEC,ICH,ISLL,IRESC
      COMMON IFL(200),LSIB(200),NSP(200),LUND(200),IGEN(200),NNN(200)
      DIMENSION MMM(200),KN2(200),N2(200),N1(200)
      CHARACTER*4 AN2,AN1,NOPAR,ID
      EQUIVALENCE (MMM(1),IGEN(1)),(KN2(1),NSP(1)),(AN2(1),N2(1)),
     1  (AN1(1),N1(1))
      LOGICAL MIST,ISEX
C
      DO 20 I=1,M
      IFL(I)=0
   20 NS(I)=0
C
      DO 55 I=1,M
      NNN(I) = I
      IF((AN1(I).NE.NOPAR).AND.(AN2(I).NE.NOPAR))  GO TO 25
      IF((AN1(I).EQ.NOPAR).AND.(AN2(I).EQ.NOPAR))  GO TO 47
      CALL ERROR( 5,I)
      RETURN
   25 IF((ID(I).EQ.AN1(I)).OR.(ID(I).EQ.AN2(I)).OR.
     1         (AN1(I).EQ.AN2(I))) GO TO 26
      GO TO 29
   26 CALL ERROR( 7,I)
      RETURN
   29 DO 30 J=1,M
      K=J
      IF(AN1(I).EQ.ID(J)) GO TO 35
   30 CONTINUE
      CALL ERROR( 1,I)
      RETURN
   35 N1(I) = K
      DO 40 J=1,M
      K=J
      IF(N2(I).EQ.ID(J)) GO TO 45
   40 CONTINUE
      CALL ERROR( 2,I)
      RETURN
   45 N2(I) = K
      GO TO 48
   47 N1(I) = 0
      N2(I) = 0
   48 IF(I.EQ.M) GO TO 55
      I1 = I+1
      DO 50 J=I1,M
      IF(ID(I).NE.ID(J)) GO TO 50
      CALL ERROR( 3,I)
      RETURN
   50 CONTINUE
   55 CONTINUE
C
      DO 61 II=1,M
      I=N1(II)
      IF(I.EQ.0) GO TO 61
      J=N2(II)
      IF(ISEX(I)) GO TO 57
      IF(ISEX(J)) GO TO 61
      GO TO 59
   57 IF(ISEX(J)) GO TO 59
      GO TO 61
   59 CALL ERROR( 8,II)
      RETURN
   61 CONTINUE
      JJ = 0
   62 KK = JJ
      DO 64 II=1,M
      IF(IFL(II).EQ.1) GO TO 64
      I = N1(II)
      IF(I.EQ.0) GO TO 63
      J = N2(II)
      K = IFL(I) + IFL(J)
      IF(K.LT.2) GO TO 64
   63 JJ = JJ+1
      NNN(JJ) = II
      IF(JJ.EQ.M) GO TO 66
      IFL(II) = 1
   64 CONTINUE
      IF(KK.LT.JJ) GO TO 62
      DO 65 I=1,M
      IF(IFL(I).NE.0) GO TO 65
      CALL ERROR( 6,I)
      RETURN
   65 CONTINUE
C RENUMBER AND REORDER EVERYONE SO THAT BOTH PARENTS COME BEFORE THEIR
C CHILDREN IN THE LIST.
   66 DO 67 I=1,M
      IFL(I) = 0
      J = NNN(I)
      MMM(J) = I
   67 CONTINUE
      DO  8067 I=1,M
      J = N1(I)
      IF(J.EQ.0) GO TO 8067
      NS(I) = MMM(J)
      J = N2(I)
      IFL(I) = MMM(J)
      IF(NS(I).LT.IFL(I)) GO TO 8067
      IFL(I) = NS(I)
      NS(I) = MMM(J)
 8067 CONTINUE
      DO  9067 I=1,M
      II = MMM(I)
      N1(II) = NS(I)
      N2(II) = IFL(I)
      NS(I) = 0
      IFL(I) = 0
      MMM(I) = 0
 9067 CONTINUE
C FIND GENERATION NUMBER OF EACH PERSON AND PUT IT IN -IGEN-.
   68 DO 70 II=1,M
      I = N1(II)
      IF(I.EQ.0) GO TO 70
      J = N2(II)
      IF(IGEN(I).GT.IGEN(J)) IGEN(J) = IGEN(I)
      IF(IGEN(I).LT.IGEN(J)) IGEN(I) = IGEN(J)
      IGEN(II) = IGEN(I) + 1
   70 CONTINUE
      DO 75 II=1,M
      I = N1(II)
      IF(I.EQ.0) GO TO 75
      J = N2(II)
      IF(IGEN(I).NE.IGEN(J)) GO TO 68
   75 CONTINUE
      IPTR = M
C SET UP SIB POINTERS
      NL = 0
   80 I = N1(IPTR)
      IF(I.EQ.0) GO TO 110
      J = N2(IPTR)
      L = IPTR
      NS(IPTR) = IPTR
      JPTR = IPTR
   90 JPTR = JPTR - 1
      IF(JPTR.LE.J) GO TO 100
      IF((N1(JPTR).NE.I).OR.(N2(JPTR).NE.J))  GO TO 90
      NS(L) = JPTR
      L = JPTR
      NS(JPTR) = IPTR
      GO TO 90
  100 NL = NL+1
      LSIB(NL) = L
  110 IPTR = IPTR - 1
      IF(IPTR.LE.2) GO TO 120
      IF(NS(IPTR)) 110,80,110
  120 NUND = 0
      DO 130 I=1,M
      NSP(I) = 0
  130 CONTINUE
C FIND NUMBER OF SPOUSES FOR EACH PARENT
      DO 140 IQ=1,NL
      IPTR = LSIB(IQ)
      I = N1(IPTR)
      J = N2(IPTR)
      NSP(I) = NSP(I) + 1
      NSP(J) = NSP(J) + 1
  140 CONTINUE
C PUT PARENTS IN CORRECT ORDER (IN -N1- AND -N2-) -- ONLY FOR
C HEAD SIBS
      DO 200 IQ = 1,NL
      IPTR = LSIB(IQ)
      I = N1(IPTR)
      J = N2(IPTR)
      IF((N1(J).GT.0).AND.(N1(I).GT.0)) GO TO 143
      GO TO 144
  143 CALL ERROR(11,IPTR)
      RETURN
  144 IF(N1(J).GT.0) GO TO 160
      IF(N1(I).GT.0) GO TO 150
      IF(NSP(I).EQ.1) GO TO 160
      IF(NSP(J).EQ.1) GO TO 150
      NUND = NUND + 1
      LUND(NUND) = IPTR
      GO TO 200
C REVERSE THE PARENTS IN -N1- AND -N2-
  150 N1(IPTR) = J
      N2(IPTR) = I
      J=I
C REDUCE THE NUMBER OF SPOUSES OF J SINCE WE HAVE JUST DETERMINED THE
C CORRECT ORDER WITH ONE OF THEM (NAMELY, I).
  160 NSP(J) = NSP(J) - 1
C FLAG THE HEAD SIB AS HAVING HIS PARENTS IN THE CORRECT ORDER IN -N1-
C AND -N2-
      IFL(IPTR) = 1
  200 CONTINUE
C PROCESS THE UNDECIDEDS
  210 IF(NUND.EQ.0) GO TO 300
      KUND = NUND
      II = 1
  215 IPTR = LUND(II)
      I = N1(IPTR)
      J = N2(IPTR)
      DO 220 JJ=1,NL
      JPTR = LSIB(JJ)
      IF(IFL(JPTR).NE.1) GO TO 220
      K = N1(JPTR)
C EACH PERSON CAN APPEAR IN -N1- OF AT MOST ONE HEAD SIB
      IF(I.EQ.K) GO TO 230
      IF(J.EQ.K) GO TO 240
  220 CONTINUE
      GO TO 260
  230 N1(IPTR) = J
      N2(IPTR) = I
  240 IFL(IPTR) = 1
      NUND = NUND-1
      DO 250 KK=II,NUND
      LUND(KK) = LUND(KK+1)
  250 CONTINUE
      GO TO 265
  260 II = II + 1
  265 IF(II.LE.NUND) GO TO 215
      IF(KUND.GT.NUND) GO TO 210
      CALL ERROR( 9,LUND(1))
      RETURN
C FILL IN ARRAY -LIST-, WHICH GIVES ORDER TO EVALUATE PEOPLE
  300 DO 301 I=1,M
      LUND(I) = 0
      IFL(I) = 0
  301 CONTINUE
C CHECK FOR TWO PEOPLE IN PEDIGREE MARRIED TO THE SAME PERSON.
      DO  8301 I=1,NL
      IPTR = LSIB(I)
      J = N1(IPTR)
      IF(IFL(J).NE.1) GO TO 7300
      CALL ERROR(12,IPTR)
      RETURN
 7300 IFL(J) = 1
 8301 CONTINUE
      DO  9301 I=1,M
      IFL(I) = 0
 9301 CONTINUE
C SET -IFL- TO POINT TO HEAD SIB.
      DO 303 I=1,NL
      IPTR = LSIB(I)
      J = IPTR
  302 IFL(J) = IPTR
      J = NS(J)
      IF(J.NE.IPTR) GO TO 302
  303 CONTINUE
C COUNT NUMBER OF SIBS IN LONG ARRAY FOR EACH SIBSHIP.  -JH- = MAXIMUM
C OF THESE NUMBERS.
      JH = 0
      DO 304 I=1,NL
      IPTR = LSIB(I)
      J = N2(IPTR)
      IF(IFL(J).GT.0) J=IFL(J)
      LUND(J) = LUND(J) + 1
      IF(LUND(J).GT.JH) JH=LUND(J)
  304 CONTINUE
C -IH- IS THE HIGHEST GENERATION NUMBER IN THE PEDIGREE.
      IH = 0
      DO 305 I=1,M
      IFL(I) = 0
      IF(IGEN(I).GT.IH) IH=IGEN(I)
  305 CONTINUE
      JH = JH + 1
      IH = IH + 1
C PUT HEAD SIBS IN -IFL- IN DECREASING ORDER OF NUMBER OF SIBS IN LONG
C ARRAY.
      K = 0
      DO 307 I=1,JH
      IX = JH-I
      DO 306 J=1,NL
      JPTR = LSIB(J)
      IF(LUND(JPTR).NE.IX) GO TO 306
      K = K+1
      IFL(K) = JPTR
  306 CONTINUE
  307 CONTINUE
C PUT HEAD SIBS IN -LSIB- IN DECREASING ORDER OF GENERATION NUMBER
C WHEN GENERATION IS THE SAME, THEY WILL BE IN DECREASING ORDER OF
C NUMBER OF SIBS IN LONG ARRAY.
      K = 0
      DO 309 I=1,IH
      IX = IH-I
      DO 308 J=1,NL
      JPTR = IFL(J)
      IF(IGEN(JPTR).NE.IX) GO TO 308
      K = K+1
      LSIB(K) = JPTR
  308 CONTINUE
  309 CONTINUE
      DO 310 I=1,M
      IFL(I) = 0
  310 CONTINUE
      DO 320 I=1,NL
      IPTR = LSIB(I)
      KN2(I) = N2(IPTR)
  320 CONTINUE
      L = 0
      IPTR=0
  325 IPTR=IPTR+1
      I = LSIB(IPTR)
      II = N1(I)
      IF(IPTR.EQ.NL) GO TO 360
      I1 = IPTR + 1
      DO 330 J=I1,NL
      IF(II.EQ.KN2(J))  GO TO 380
  330 CONTINUE
  360 L = L+1
      LIST(L) = I
      IFL(I) = 1
      IFL(II) = 1
      J = NS(I)
  370 IF(J.EQ.I) GO TO 400
      IFL(J) = 1
      J = NS(J)
      GO TO 370
  380 J1 = J-1
      DO 390 K=IPTR,J1
      LSIB(K) = LSIB(K+1)
      KN2(K) = KN2(K+1)
  390 CONTINUE
      LSIB(J) = I
      KN2(J) = N2(I)
      IPTR = IPTR-1
  400 IF(IPTR.LT.NL) GO TO 325
C ADD TO -LIST- ALL THOSE NOT FLAGGED.  THIS WILL INCLUDE SPOUSES WITH
C NO CHILDREN AND PARENTS AT THE TOP OF THE PEDIGREE AND AT THE TOP
C OF A SUB-FAMILY.
      DO 410 I=1,M
      IF(IFL(I).EQ.1) GO TO 410
      L = L+1
      LIST(L) = I
  410 CONTINUE
C  ASSIGN OLD NUMBERS (ORDER IN WHICH PERSONS WERE READ)
      DO 420 I=1,M
      J=NNN(I)
      IFL(J)=I
      LSIB(I)=N2(I)
      NSP(I)=N1(I)
  420 LUND(I)=NS(I)
      DO 440 I=1,M
      K=IFL(I)
      IF(LSIB(K).EQ.0) GO TO 430
      J=LSIB(K)
      N2(I)=NNN(J)
      J=NSP(K)
      N1(I)=NNN(J)
      J=LUND(K)
      NS(I)=NNN(J)
      GO TO 440
  430 N2(I)=0
      N1(I)=0
      NS(I)=0
  440 CONTINUE
      DO 450 I=1,L
  450 IGEN(I)=LIST(I)
      DO 460 I=1,L
      J=IGEN(I)
  460 LIST(I)=NNN(J)
C  PUT NS(I)=I IF I IS SINGLE OFFSPRING
      DO 470 I=1,M
      IF(NS(I).NE.0) GO TO 470
      NS(I)=I
  470 CONTINUE
      RETURN
      END
C     ----------------------------------------
      SUBROUTINE ERROR(NO,NUM)
C  CALLED FROM SUBROUTINE SETUP
      COMMON /B/ LIST(50),AN2(200),AN1(200),NS(200),M,L,NOPAR,ID(200),
     1  NDI,LDI(2,3),LGC(3),ISC,MIST,ISEX(200),
     1  FENO1(200,10),FENO2(200,21),PHI(300),RMUT,
     1  HOLD(3,300),RLIK,PHIS(36),WVK(4),THM,THF,IAD(36,36),
     2  IV(4),JV(4),IVA1(4),IVA2(16),IVA3(4),IVF(4),NGEN,NGENO,MNP,
     3  MNPT,KONT1,KONT2,NAL1,NAL2,MVEC,ICH,ISLL,IRESC
      COMMON IFL(200),LSIB(200),NSP(200),LUND(200),IGEN(200),NNN(200)
      CHARACTER*4 AN2,AN1,NOPAR,ID
      LOGICAL MIST,ISEX
    1 FORMAT(/13H ERROR NUMBER,I3,
     1 29H IN SETUP, DETECTED AT PERSON,I4,3X,A4)
C
      MIST=.TRUE.
      IF(NO.LT.9) GO TO 10
      I=NNN(NUM)
      GO TO 20
   10 I=NUM
   20 WRITE(6,1) NO,I,ID(I)
      RETURN
      END
