Show compileandrun.sh syntax highlighted
#it makes sense to run this with parameter 1.
#sh compileandrun.sh 1
#takes about 5 minutes on 400MHz.
cat <<EOF > simulation.f
C hkka100.f : SELECCION.
C MUESTREO DE GIBBS. DATA AUGMENTATION.
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 IG(40000),ITO(200000,2),IGH(40000),NPA(40000,2)
INTEGER*4 IGHI(40000),NPAT(40000,2),IGP(40000),IB(100),ISO(100)
INTEGER*4 NPAT2(40000,2),IGP2(40000),IL(100),IS(100),LSP(100)
REAL*8 RITO(200000),UP(40000),YY(40000),BP(40000),SS(100)
REAL*8 U(40000),UH(40000),UHI(40000),BTP(40000),S(100),RIB(100)
REAL*8 YYH(40000),YYHI(40000),B(100),BT(2),RL(100),RLS(100)
REAL*8 YY2(40000),BP2(40000),BTP2(40000),UP2(40000),PRO(100)
REAL*8 PRF1(100),RLSP(100),SO(100)
CHARACTER*20 SOURCEFILENAME
CHARACTER*10 ITERATIONS
CHARACTER*10 BUFFER
INTEGER MYINT
CALL GetArg(1, ITERATIONS)
CALL GetArg(2, SOURCEFILENAME)
WRITE(BUFFER, '(A10)'), ITERATIONS
READ(BUFFER, '(I10.1)'), MYINT
OPEN(1,FILE='HRESA100',STATUS='NEW')
C X1=0.3325174856D+09
C X1=0.2484672337D+09
X1=0.6924612876D+09
C X1=0.4321164154D+09
C X1=0.5294523875D+09
C X1=0.7563432787D+09
WRITE(1,*) 'RESULTADOS OBTENIDOS CON "', SOURCEFILENAME, '"'
WRITE(1,*) 'MODELO TOTALMENTE ALEATORIO'
WRITE(1,*) 'MUESTREO DE GIBBS '
WRITE(1,*) 'DATA AUGMENTATION MONOCATENARIO '
WRITE(1,*) 'SELECCIONADOS EN BASE CONOCIDOS '
WRITE(1,*) 'ITERACIONES=', BUFFER
C PARAMETROS INICIALES.
NESLA=10000
IREPE=MYINT
C NESLA: Numero de eslabones de cada cadena.
C IREPE: Numero de casos.
NMA=1000
NHE=1000
NSEL=500
NS=NMA-NSEL
NHI=3000
NNSH=1500
NSH=NHI-NNSH
NR=NMA+NHE+NHI
VU=20.
VE=80.
VAL=10.
DVU=VU**.5
DVE=VE**.5
DVAL=VAL**.5
BT(1)=0.
BT(2)=0.
C NMA: Numero de machos de la generacion base.
C NHE: Numero de hembras " " .
C NSEL: Numero de seleccionados " .
C NHI: Numero de hijos de la primera generacion.
C NNSH: Numero de hijos no seleccionados (individuos sin dato).
C NR: Numero total de individuos.
C VU: Varianza genetica.
C VE: Varianza residual.
C DVU: Mejor no digo nada...
C DVE:
C B(1): Valor del primer nivel fijo.
C B(2): " segundo nivel fijo.
C B(3): " tercer nivel fijo.
C VP Y RP SON LA VARIANZA Y LOS G.L. A PRIORI DEL COMPONENTE.
VP=0.
RP=-2.
C VEP Y RE SON LA VARIANZA Y LOS G.L. A PRIORI DEL RESIDUO.
VEP=0.
RE=-2.
WRITE(1,*) 'MACHOS HEMBRAS HIJOS ',NMA,NHE,NHI
WRITE(1,*) 'SELECCIONADOS ',NSEL
WRITE(1,*) 'HIJOS SIN OBSERVACION ',NNSH
WRITE(1,*) 'VAL BT2 VA VE '
WRITE(1,'(4F7.1)') VAL,BT(2),VU,VE
WRITE(1,*)
CLOSE(1)
SINR=0.
SB1=0.
SB2=0.
SB3=0.
SBT1=0.
SBT2=0.
SVA=0.
SVER=0.
ST=0.
S2B1=0.
S2B2=0.
S2B3=0.
S2VA=0.
S2VER=0.
SP1=0.
SP2=0.
SP3=0.
SH1=0.
SH2=0.
SH3=0.
DO 1000 IJK=1,IREPE
ICTE=0
DO 5 I=1,100
IB(I)=1
ICTE=ICTE+IB(I)
CALL NORMAL(Z,X1)
B(I)=Z*DVAL
5 CONTINUE
C Ordena el vector de aleatorios, ya que es necesario para el módulo
C de GENERACION DE ALEATORIOS FALTANTES.
NNB=100
CALL PIKSRT(NNB,B)
RMIN=B(1)
RMAX=B(100)
CTE=DBLE(ICTE)
SUM=0.
DO 7 I=1,100
SUM=SUM+DBLE(IB(I))/CTE
7 RIB(I)=SUM
C OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,'(4F7.1)') B(1),B(50),B(100)
C WRITE(1,*)
C WRITE(1,'(2F8.3)') RMIN,RMAX
C WRITE(1,*)
C WRITE(1,'(10I7)') (IB(I),I=1,10)
C WRITE(1,'(10I7)') (IB(I),I=46,55)
C WRITE(1,'(10I7)') (IB(I),I=91,100)
C WRITE(1,*)
CLOSE(1)
DO 10 I=1,40000
YY(I)=0.
YY2(I)=0.
YYH(I)=0.
YYHI(I)=0.
U(I)=0.
UH(I)=0.
UHI(I)=0.
IG(I)=0
IGH(I)=0
IGHI(I)=0
NPA(I,1)=0
NPA(I,2)=0
NPAT(I,1)=0
NPAT(I,2)=0
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP(I)=0
IGP2(I)=0
BP(I)=0.
BP2(I)=0.
BTP(I)=0.
BTP2(I)=0.
UP2(I)=0.
10 UP(I)=0.
DO 15 I=1,100
RL(I)=0.
RLS(I)=0.
RLSP(I)=0.
SS(I)=0.
IS(I)=0
LSP(I)=0
ISO(I)=0
SO(I)=0.
15 IL(I)=0
DIF=RMAX-RMIN
CTE=DIF/5.
NC1=0
NC2=0
NC3=0
NC4=0
NC5=0
C Generacion de machos de la generacion base.
DO 20 I=1,NMA
CALL NORMAL(Z,X1)
U(I)=Z*DVU
CALL RANDOM(X1,Z)
J=INT(Z*100.)+1
IG(I)=J
IL(J)=IL(J)+1
CALL NORMAL(ZZ,X1)
YY(I)=B(IG(I))+BT(1)+U(I)+ZZ*DVE
IF (B(IG(I)).LT.(RMIN+CTE)) THEN
NC1=NC1+1
ELSE IF (B(IG(I)).LT.(RMIN+CTE*2)) THEN
NC2=NC2+1
ELSE IF (B(IG(I)).LT.(RMIN+CTE*3)) THEN
NC3=NC3+1
ELSE IF (B(IG(I)).LT.(RMIN+CTE*4)) THEN
NC4=NC4+1
ELSE
NC5=NC5+1
ENDIF
20 CONTINUE
C Ordeno a los machos segun su fenotipo.
CALL SORT(NMA,IG,IGP,U,YY)
TB=YY(NS+1)
C Genero a la hembras de la generacion base.
DO 410 I=1,NHE
CALL NORMAL(Z,X1)
UH(I)=Z*DVU
CALL RANDOM(X1,Z)
J=INT(Z*100.)+1
IGH(I)=J
IL(J)=IL(J)+1
CALL NORMAL(ZZ,X1)
YYH(I)=B(IGH(I))+BT(1)+UH(I)+ZZ*DVE
IF (B(IGH(I)).LT.(RMIN+CTE)) THEN
NC1=NC1+1
ELSE IF (B(IGH(I)).LT.(RMIN+CTE*2)) THEN
NC2=NC2+1
ELSE IF (B(IGH(I)).LT.(RMIN+CTE*3)) THEN
NC3=NC3+1
ELSE IF (B(IGH(I)).LT.(RMIN+CTE*4)) THEN
NC4=NC4+1
ELSE
NC5=NC5+1
ENDIF
410 CONTINUE
DO 415 I=1,100
415 RL(I)=DBLE(IL(I))/DBLE(NMA+NHE)
C Ordeno a las hembras segun su fenotipo.
CALL SORT(NHE,IGH,IGP,UH,YYH)
TBH=YYH(NS+1)
C ¡No volver a desactivar nunca jamás!
C IF (TBH.LT.TB) THEN
C TB=TBH
C ELSE
C TBH=TB
C ENDIF
C Genero a los hijos a partir de los padres seleccionados.
R2=(VU/2.)**.5
DO 3050 I=1,NHI
CALL RANDOM(X1,Z)
CALL RANDOM(X1,Z1)
NPA(I,1)=NS+INT(Z*NSEL)+1
NPA(I,2)=NS+INT(Z1*NSEL)+1
CALL RANDOM(X1,Z)
J=INT(Z*100.)+1
IGHI(I)=J
CALL NORMAL(Z2,X1)
UHI(I)=Z2*R2+(U(NPA(I,1))+UH(NPA(I,2)))/2.
CALL NORMAL(Z,X1)
YYHI(I)=B(IGHI(I))+BT(2)+UHI(I)+Z*DVE
IF (B(IGHI(I)).LT.(RMIN+CTE)) THEN
NC1=NC1+1
ELSE IF (B(IGHI(I)).LT.(RMIN+CTE*2)) THEN
NC2=NC2+1
ELSE IF (B(IGHI(I)).LT.(RMIN+CTE*3)) THEN
NC3=NC3+1
ELSE IF (B(IGHI(I)).LT.(RMIN+CTE*4)) THEN
NC4=NC4+1
ELSE
NC5=NC5+1
ENDIF
3050 CONTINUE
C Ordeno a los hijos segun fenotipo.
CALL SORTI(NHI,NPA,IGHI,IGP,UHI,YYHI)
X2=X1
C DO 422 I=NS+1,NMA
DO 422 I=1,NMA
422 RLS(IG(I))=RLS(IG(I))+1.
C DO 424 I=NS+1,NHE
DO 424 I=1,NHE
424 RLS(IGH(I))=RLS(IGH(I))+1.
C DO 427 I=NNSH+1,NHI
DO 427 I=1,NHI
427 RLS(IGHI(I))=RLS(IGHI(I))+1.
C OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,*)
C WRITE(1,'(20I4)') (INT(RLS(I)),I=1,100)
C WRITE(1,*)
C CLOSE(1)
DO 5421 I=1,100
5421 RLS(I)=0.
DO 5422 I=NS+1,NMA
5422 RLS(IG(I))=RLS(IG(I))+1.
DO 5424 I=NS+1,NHE
5424 RLS(IGH(I))=RLS(IGH(I))+1.
DO 5427 I=NNSH+1,NHI
5427 RLS(IGHI(I))=RLS(IGHI(I))+1.
C Detecta los niveles del aleatorio sin observaciones.
ILSP=0
DO 428 I=1,100
LSP(I)=0
IF (INT(RLS(I)).GT.0) THEN
LSP(I)=1
ILSP=ILSP+1
ENDIF
428 CONTINUE
C Los grados de libertad del aleatorio dependen del número de
C niveles observados.
RUL=DBLE(ILSP)+(-2.)
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,*)
C WRITE(1,'(20I4)') (INT(RLS(I)),I=1,100)
C WRITE(1,'(20I4)') (INT(RLS(I)),I=41,60)
C WRITE(1,'(20I4)') (INT(RLS(I)),I=81,100)
C WRITE(1,*)
WRITE(1,*)
WRITE(1,'(3F7.1,I8)') B(1),B(50),B(100),ILSP
C WRITE(1,'(I8)') ILSP
WRITE(1,*)
CLOSE(1)
NRLS=0
DO 1510 I=1,100
IF (RLS(I).NE.0) NRLS=NRLS+1
1510 RLS(I)=RLS(I)/DBLE(NSEL+NSEL+NSH)
PRIF1=1./DBLE(NRLS)
DO 1515 I=1,100
PRF1(I)=0.
IF (RLS(I).NE.0) PRF1(I)=PRIF1
1515 CONTINUE
L1=0
L2=0
L3=0
LL1=0
LL2=0
LL3=0
DO 423 I=1,NHI
IF (IGHI(I).EQ.1) THEN
L1=L1+1
ELSE IF (IGHI(I).EQ.50) THEN
L2=L2+1
ELSE IF (IGHI(I).EQ.100) THEN
L3=L3+1
ENDIF
IF (I.GT.NNSH) THEN
IF (IGHI(I).EQ.1) THEN
LL1=LL1+1
ELSE IF (IGHI(I).EQ.50) THEN
LL2=LL2+1
ELSE IF (IGHI(I).EQ.100) THEN
LL3=LL3+1
ENDIF
ENDIF
423 CONTINUE
CTE=DBLE(NHI)
RHS1=DBLE(L1)/CTE
RHS2=DBLE(L2)/CTE
RHS3=DBLE(L3)/CTE
CTE=DBLE(NSH)
RHL1=DBLE(LL1)/CTE
RHL2=DBLE(LL2)/CTE
RHL3=DBLE(LL3)/CTE
C OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,'(5I6)') NC1,NC2,NC3,NC4,NC5
C WRITE(1,*)
C WRITE(1,'(6F8.3)') RL(1),RL(50),RL(100),RLS(1),RLS(50),RLS(100)
C WRITE(1,'(6F8.3)') RHS1,RHS2,RHS3,RHL1,RHL2,RHL3
C WRITE(1,*)
C CLOSE(1)
SRLS=0.
DO 426 I=1,100
RLSP(I)=RLS(I)
SRLS=SRLS+RLSP(I)
426 CONTINUE
SP1=SP1+RLS(1)
SP2=SP2+RLS(50)
SP3=SP3+RLS(100)
SH1=SH1+RHL1
SH2=SH2+RHL2
SH3=SH3+RHL3
C Asigno padres desconocidos (cero) a machos y hembras
C de la generacion base y pego las hembras despues de
C los machos.
DO 425 I=1,NMA
NPAT(I,1)=0
425 NPAT(I,2)=0
J=NMA+1
DO 430 I=1,NHE
YY(J)=YYH(I)
IG(J)=IGH(I)
U(J)=UH(I)
NPAT(J,1)=0
NPAT(J,2)=0
430 J=J+1
SNB1=0.
SNB2=0.
SNB3=0.
SNVA=0.
SNVER=0.
SNT=0.
SN2B1=0.
SN2B2=0.
SN2B3=0.
SN2VA=0.
SN2VER=0.
SN2T=0.
C VARIANZAS ADITIVA Y RESIDUAL INICIALES.
VUA=10.
VUE=60.
VALE=5.
C VALORES ADITIVOS INICIALES. Es completamente aleatorio y solo
C depende del componente aditivo inicial.
DO 150 I=1,NMA+NHE
150 UP(I)=0.
C150 UP(I)=U(I)
C AÑADO HIJOS ORDENADOS DESCENDENTEMENTE Y GENERO
C PADRES Y ALEATORIO DE LOS HIJOS NO SELECCIONADOS.
R2=(VUA/2.)**.5
J=NMA+NHE+1
DO 450 I=NHI,1,-1
YY(J)=YYHI(I)
IG(J)=IGHI(I)
U(J)=UHI(I)
UP(J)=0.
IF (J.GT.NMA+NHE+NSH) THEN
CALL RANDOM(X2,Z)
NPAT(J,1)=NS+INT(Z*DBLE(NSEL))+1
CALL RANDOM(X2,Z)
NPAT(J,2)=NMA+NS+INT(Z*DBLE(NSEL))+1
4457 CALL RANDOM(X2,Z)
K=INT(Z*100.)+1
C Impide muestrear niveles sin observaciones.
IF (LSP(K).EQ.1) THEN
IG(J)=K
ELSE
GO TO 4457
ENDIF
ELSE
NPAT(J,1)=NPA(I,1)
NPAT(J,2)=NMA+NPA(I,2)
ENDIF
450 J=J+1
C Valores iniciales del aleatorio de los padres no seleccionados.
DO 160 I=1,NS
1607 CALL RANDOM(X2,Z)
J=INT(Z*100.)+1
IF (LSP(J).EQ.1) THEN
IG(I)=J
ELSE
GO TO 1607
ENDIF
160 CONTINUE
DO 165 I=NMA+1,NMA+NS
1657 CALL RANDOM(X2,Z)
J=INT(Z*100.)+1
IF (LSP(J).EQ.1) THEN
IG(I)=J
ELSE
GO TO 1657
ENDIF
165 CONTINUE
C VALORES INICIALES DEL ALEATORIO Y DEL FIJO INTRAGENERACIONAL.
DO 170 I=1,NR
BP(I)=0.
IF (I.LE.(NMA+NHE)) THEN
BTP(I)=0.
ELSE
BTP(I)=0.
ENDIF
170 CONTINUE
DO 175 I=1,NR
175 IGP(I)=IG(I)
C VALOR INICIAL DEL UMBRAL. Esta en funcion del hijo selec-
C cionado con menor valor fenotipico. En esta decision perdi un
C par de neuronas.
RMIN=YY(NMA+NHE+NSH)
T1=RMIN
C Valor inicial de la media poblacional.
RMP=0.
B1=0.
B2=0.
B3=0.
BT1=0.
BT2=0.
VA=0.
VER=0.
T=0.
SNRM=0.
SNUF=0.
SNUM=0.
IL1=0
IL2=0
IL3=0
SRG1=0.
SRG2=0.
SRG3=0.
SRLH1=0.
SRLH2=0.
SRLH3=0.
SRS1=0.
SRS2=0.
SRS3=0.
NBM=500
NBH=500
C NBM=700
C NBH=700
NGB=NBM+NSEL+NBH+NSEL
N1H=1500
C N1H=2100
NG1=NSH+N1H
NRM=NGB+NG1
DO 510 I=1,NBM
YY2(I)=YY(I)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(I)
BP2(I)=BP(I)
BTP2(I)=BTP(I)
510 UP2(I)=UP(I)
J=NS+1
DO 520 I=NBM+1,NBM+NSEL
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
520 J=J+1
J=NMA+1
DO 530 I=NBM+NSEL+1,NBM+NSEL+NBH
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
530 J=J+1
J=NMA+NS+1
DO 540 I=NBM+NSEL+NBH+1,NGB
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
540 J=J+1
J=NMA+NHE+1
IP=NS-NBM
IM=NMA+NS-(NBM+NSEL+NBH)
DO 550 I=NGB+1,NGB+NSH
YY2(I)=YY(J)
NPAT2(I,1)=NPAT(J,1)-IP
NPAT2(I,2)=NPAT(J,2)-IM
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
550 J=J+1
J=NMA+NHE+NSH+1
DO 560 I=NGB+NSH+1,NRM
YY2(I)=YY(J)
NPAT2(I,1)=NPAT(J,1)-IP
NPAT2(I,2)=NPAT(J,2)-IM
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
560 J=J+1
DO 565 I=1,NRM
YY(I)=YY2(I)
NPAT(I,1)=NPAT2(I,1)
NPAT(I,2)=NPAT2(I,2)
IGP(I)=IGP2(I)
BP(I)=BP2(I)
BTP(I)=BTP2(I)
565 UP(I)=UP2(I)
DO 570 I=1,100
RLS(I)=0.
570 IS(I)=0
DO 575 I=1,NRM
575 IS(IGP(I))=IS(IGP(I))+1
C OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,'(I8)') ILSP
C WRITE(1,*)
C WRITE(1,'(10I7)') (IS(I),I=1,10)
C WRITE(1,'(10I7)') (IS(I),I=46,55)
C WRITE(1,'(10I7)') (IS(I),I=91,100)
C WRITE(1,*)
C WRITE(1,'(10F7.3)') (RLSP(I),I=1,20)
C WRITE(1,'(F7.3)') SRLS
C WRITE(1,*)
C CLOSE(1)
NNU=0
NNF=0
C Inicio el proceso Gibbs.
C BUCLE DE ESLABON.
KL=0
2000 KL=KL+1
C Genero la inversa de la matriz de parentesco.
DO 30 I=1,200000
RITO(I)=0.
DO 30 J=1,2
30 ITO(I,J)=0
KK=NRM
DO 100 I=1,NRM
ITO(I,1)=I
ITO(I,2)=I
IF (I.LE.(NGB)) THEN
RITO(I)=1.
ELSE
I1=NPAT(I,1)
I2=NPAT(I,2)
RITO(I)=2.
KK=KK+1
ITO(KK,1)=I
ITO(KK,2)=I1
RITO(KK)=-1.
KK=KK+1
ITO(KK,1)=I1
ITO(KK,2)=I
RITO(KK)=-1.
KK=KK+1
ITO(KK,1)=I
ITO(KK,2)=I2
RITO(KK)=-1.
KK=KK+1
ITO(KK,1)=I2
ITO(KK,2)=I
RITO(KK)=-1.
RITO(I1)=RITO(I1)+.5
KK=KK+1
ITO(KK,1)=I1
ITO(KK,2)=I2
RITO(KK)=.5
KK=KK+1
ITO(KK,1)=I2
ITO(KK,2)=I1
RITO(KK)=.5
RITO(I2)=RITO(I2)+.5
ENDIF
100 CONTINUE
IF (KK.GT.200000) THEN
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
WRITE(1,*) 'MATRIZ DE PARENTESCO DESBORDADA '
CLOSE(1)
GO TO 10000
ENDIF
CALL SORT2(KK,ITO,RITO)
C GENERACION DE DATOS FALTANTES.
DVUE=VUE**.5
C Generación base:
DO 190 I=1,NBM
RM=RMP+BP(I)+BTP(I)+UP(I)
P1=PNOR((TB-RM)/DVUE)
CALL RANDOM(X2,Z)
190 YY(I)=RM+XNOR(Z*P1)*DVUE
DO 195 I=NBM+NSEL+1,NBM+NSEL+NBH
RM=RMP+BP(I)+BTP(I)+UP(I)
P1=PNOR((TBH-RM)/DVUE)
CALL RANDOM(X2,Z)
195 YY(I)=RM+XNOR(Z*P1)*DVUE
C Generación filial.
DO 200 I=NGB+NSH+1,NRM
RM=RMP+BP(I)+BTP(I)+UP(I)
P1=PNOR((T1-RM)/DVUE)
CALL RANDOM(X2,Z)
YY(I)=RM+XNOR(Z*P1)*DVUE
200 CONTINUE
C MUESTREO DE LA MEDIA.
CTE=0.
DO 5210 I=1,NRM
5210 CTE=CTE+YY(I)-BP(I)-BTP(I)-UP(I)
CTE=CTE/DBLE(NRM)
CALL NORMAL(Z,X2)
RMP=CTE+Z*(VUE/DBLE(NRM))**.5
C MUESTREO DEL ALEATORIO.
ALFAL=VUE/VALE
DO 205 I=1,100
S(I)=0.
205 IS(I)=0
DO 210 I=1,NRM
IS(IGP(I))=IS(IGP(I))+1
S(IGP(I))=S(IGP(I))+YY(I)-RMP-BTP(I)-UP(I)
210 CONTINUE
DO 215 I=1,100
IF (IS(I).GT.0) THEN
RIC=DBLE(IS(I))+ALFAL
S(I)=S(I)/RIC
CALL NORMAL(Z,X2)
S(I)=S(I)+Z*(VUE/RIC)**.5
ENDIF
215 CONTINUE
DO 220 I=1,NRM
220 BP(I)=S(IGP(I))
C MUESTREO DEL FIJO INTRAGENERACIONAL.
S1=0.
S2=0.
GO TO 876
C DO 223 I=1,NGB
C223 S1=S1+YY(I)-BP(I)-UP(I)
C S1=S1/DBLE(NGB)
DO 224 I=NGB+1,NRM
224 S2=S2+YY(I)-RMP-BP(I)-UP(I)
S2=S2/DBLE(NRM-NGB)
C CALL NORMAL(Z,X2)
C S1=S1+Z*(VUE/DBLE(NGB))**.5
CALL NORMAL(Z,X2)
S2=S2+Z*(VUE/DBLE(NRM-NGB))**.5
876 DO 225 I=1,NGB
225 BTP(I)=S1
DO 226 I=NGB+1,NRM
226 BTP(I)=S2
TP1=S1
TP2=S2
C MUESTREO DE ALEATORIOS.
IK2=1
ALFA=VUE/VUA
DO 230 I=1,NRM
CA=0.
DO WHILE(ITO(IK2,1).EQ.I)
IF (ITO(IK2,2).EQ.I) THEN
CTE=RITO(IK2)
ELSE
CA=CA+RITO(IK2)*UP(ITO(IK2,2))
ENDIF
IK2=IK2+1
ENDDO
RIC=1.+CTE*ALFA
RAL=(YY(I)-RMP-BP(I)-BTP(I)-CA*ALFA)/RIC
CALL NORMAL(Z,X2)
UP(I)=RAL+Z*(VUE/RIC)**.5
230 CONTINUE
C MUESTREO DEL COMPONENTE ADITIVO.
C CALCULO DE LA SUMA DE CUADRADOS.
CTT4=0.
DO 240 I=1,KK
I1=ITO(I,1)
I2=ITO(I,2)
IF (I1.EQ.I2) THEN
CTT4=CTT4+RITO(I)*UP(I1)**2
ELSE
CTT4=CTT4+RITO(I)*UP(I1)*UP(I2)
ENDIF
240 CONTINUE
C Se genera un valor del componente segun exp.S5 de Janss et al.(1995), pag.1141.
C La variable Chi cuadrado se genera segun exp.(3.6.94) de Rubinstein, pag.94.
RUD=DBLE(NRM)+RP
CALL NORMAL(Z,X2)
CTE=((Z+(2.*RUD-1.)**.5)**2)/2.
VUA=CTT4/CTE
C MUESTREO DEL COMPONENTE RESIDUAL.
C CALCULO DE LA SUMA DE CUADRADOS.
CTT4=0.
DO 250 I=1,NRM
250 CTT4=CTT4+(YY(I)-RMP-BP(I)-BTP(I)-UP(I))**2
C Se genera un valor del componente segun exp.S4 de Janss et al.(1995).
C La variable Chi cuadrado se genera segun Rubinstein, pag.94.
RUE=DBLE(NRM)+RE
CALL NORMAL(Z,X2)
CTE=((Z+(2.*RUE-1.)**.5)**2)/2.
VUE=CTT4/CTE
C MUESTREO DEL COMPONENTE DEL FACTOR ALEATORIO.
CTT4=0.
DO 252 I=1,100
252 CTT4=CTT4+S(I)**2
CALL NORMAL(Z,X2)
CTE=((Z+(2.*RUL-1.)**.5)**2)/2.
VALE=CTT4/CTE
C GENERACION DE PADRES.
DO 270 I=NGB+NSH+1,NRM
ICTE=0
2221 CALL RANDOM(X2,Z)
NUMM=NBM+INT(Z*DBLE(NSEL))+1
CALL RANDOM(X2,Z)
NUMH=NBM+NSEL+NBH+INT(Z*DBLE(NSEL))+1
RMED=(UP(NUMM)+UP(NUMH))/2.
CTE=EXP(-1.*(UP(I)-RMED)**2/VUA)
CALL RANDOM(X2,Z)
ICTE=ICTE+1
IF (ICTE.GT.5000) GO TO 270
IF (Z.GT.CTE) GO TO 2221
NPAT(I,1)=NUMM
NPAT(I,2)=NUMH
270 CONTINUE
C GENERACION DE ALEATORIOS FALTANTES.
C Método de aceptación-rechazo.
DO 4262 I=1,NBM
DENO=-.5*(YY(I)-RMP-BP(I)-BTP(I)-UP(I))**2/VUE
C DENO=DENO-.5*BP(I)**2/VALE
2657 CALL RANDOM(X2,Z)
J=INT(Z*100.)+1
IF (LSP(J).EQ.0) GO TO 2657
RNUM=-.5*(YY(I)-RMP-S(J)-BTP(I)-UP(I))**2/VUE
C RNUM=RNUM-.5*S(J)**2/VALE
RAT=RNUM-DENO
IF (RAT.GE.0.) THEN
IGP(I)=J
BP(I)=S(J)
ELSE IF (RAT.GT.-9.) THEN
RAT=EXP(RAT)
CALL RANDOM(X2,Z)
IF (Z.LE.RAT) THEN
IGP(I)=J
BP(I)=S(J)
ENDIF
ENDIF
4262 CONTINUE
DO 4264 I=NBM+NSEL+1,NBM+NSEL+NBH
DENO=-.5*(YY(I)-RMP-BP(I)-BTP(I)-UP(I))**2/VUE
C DENO=DENO-.5*BP(I)**2/VALE
2658 CALL RANDOM(X2,Z)
J=INT(Z*100.)+1
IF (LSP(J).EQ.0) GO TO 2658
RNUM=-.5*(YY(I)-RMP-S(J)-BTP(I)-UP(I))**2/VUE
C RNUM=RNUM-.5*S(J)**2/VALE
RAT=RNUM-DENO
IF (RAT.GE.0.) THEN
IGP(I)=J
BP(I)=S(J)
ELSE IF (RAT.GT.-9.) THEN
RAT=EXP(RAT)
CALL RANDOM(X2,Z)
IF (Z.LE.RAT) THEN
IGP(I)=J
BP(I)=S(J)
ENDIF
ENDIF
4264 CONTINUE
DO 4260 I=NGB+NSH+1,NRM
DENO=-.5*(YY(I)-RMP-BP(I)-BTP(I)-UP(I))**2/VUE
C DENO=DENO-.5*BP(I)**2/VALE
2659 CALL RANDOM(X2,Z)
J=INT(Z*100.)+1
IF (LSP(J).EQ.0) GO TO 2659
RNUM=-.5*(YY(I)-RMP-S(J)-BTP(I)-UP(I))**2/VUE
C RNUM=RNUM-.5*S(J)**2/VALE
RAT=RNUM-DENO
IF (RAT.GE.0.) THEN
IGP(I)=J
BP(I)=S(J)
ELSE IF (RAT.GT.-9.) THEN
RAT=EXP(RAT)
CALL RANDOM(X2,Z)
IF (Z.LE.RAT) THEN
IGP(I)=J
BP(I)=S(J)
ENDIF
ENDIF
4260 CONTINUE
DO 388 I=1,100
388 RLS(I)=0.
DO 383 I=1,NRM
383 RLS(IGP(I))=RLS(IGP(I))+1.
DO 382 I=1,100
382 RLS(I)=RLS(I)/DBLE(NRM)
L1=0
L2=0
L3=0
LH1=0
LH2=0
LH3=0
LS1=0
LS2=0
LS3=0
DO 380 I=1,NRM
IF (IGP(I).EQ.1) THEN
L1=L1+1
ELSE IF (IGP(I).EQ.50) THEN
L2=L2+1
ELSE IF (IGP(I).EQ.100) THEN
L3=L3+1
ENDIF
IF (I.LE.NGB) THEN
IF (IGP(I).EQ.1) THEN
LH1=LH1+1
ELSE IF (IGP(I).EQ.50) THEN
LH2=LH2+1
ELSE IF (IGP(I).EQ.100) THEN
LH3=LH3+1
ENDIF
ELSE
IF (IGP(I).EQ.1) THEN
LS1=LS1+1
ELSE IF (IGP(I).EQ.50) THEN
LS2=LS2+1
ELSE IF (IGP(I).EQ.100) THEN
LS3=LS3+1
ENDIF
ENDIF
380 CONTINUE
RG1=DBLE(L1)/DBLE(NRM)
RG2=DBLE(L2)/DBLE(NRM)
RG3=DBLE(L3)/DBLE(NRM)
RLH1=DBLE(LH1)/DBLE(NGB)
RLH2=DBLE(LH2)/DBLE(NGB)
RLH3=DBLE(LH3)/DBLE(NGB)
RRS1=DBLE(LS1)/DBLE(NG1)
RRS2=DBLE(LS2)/DBLE(NG1)
RRS3=DBLE(LS3)/DBLE(NG1)
GO TO 2345
C CALCULO NUMERO DE OBSERVACIONES FALTANTES.
C GENERACION BASE.
9876 SUM=0.
DVUE=VUE**.5
DVUA=VUA**.5
DO 610 I=1,NBM+NSEL
RM=BP(I)+BTP(I)+UP(I)
610 SUM=SUM+PNOR((TB-RM)/DVUE)
DO 615 I=NBM+NSEL+1,NGB
RM=BP(I)+BTP(I)+UP(I)
615 SUM=SUM+PNOR((TBH-RM)/DVUE)
SUM=SUM/DBLE(NGB)
NUM=INT((DBLE(NSEL*2)*SUM/(1.-SUM))/2.)
C IF (NUM.GT.1500) NUM=1500
IF (NUM.GT.7000) THEN
NUM=7000
NNU=NNU+1
ENDIF
IF (NUM.GT.NBM) THEN
C Añado machos no seleccionados.
DO 620 I=1,NBM
YY2(I)=YY(I)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(I)
BP2(I)=BP(I)
BTP2(I)=BTP(I)
620 UP2(I)=UP(I)
DO 630 I=NBM+1,NUM
222 NPAT2(I,1)=0
NPAT2(I,2)=0
2223 CALL RANDOM(X2,UNI)
K=INT(UNI*100.)+1
IF (IS(K).EQ.0) GO TO 2223
IGP2(I)=K
BP2(I)=S(K)
BTP2(I)=TP1
CALL NORMAL(Z,X2)
UP2(I)=Z*DVUA
CALL NORMAL(Z,X2)
YY2(I)=BP2(I)+BTP2(I)+UP2(I)+Z*DVUE
IF (YY2(I).GT.TB) GO TO 222
630 CONTINUE
C Añado machos seleccionados.
J=NBM+1
DO 640 I=NUM+1,NUM+NSEL
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=TP1
UP2(I)=UP(J)
640 J=J+1
C Añado hembras no seleccionadas.
J=NBM+NSEL+1
DO 650 I=NUM+NSEL+1,NUM+NSEL+NBH
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
650 J=J+1
DO 660 I=NUM+NSEL+NBH+1,NUM+NSEL+NUM
333 NPAT2(I,1)=0
NPAT2(I,2)=0
3333 CALL RANDOM(X2,UNI)
K=INT(UNI*100.)+1
IF (IS(K).EQ.0) GO TO 3333
IGP2(I)=K
BP2(I)=S(K)
BTP2(I)=TP1
CALL NORMAL(Z,X2)
UP2(I)=Z*DVUA
CALL NORMAL(Z,X2)
YY2(I)=BP2(I)+BTP2(I)+UP2(I)+Z*DVUE
IF (YY2(I).GT.TBH) GO TO 333
660 CONTINUE
C Añado hembras seleccionadas.
J=NBM+NSEL+NBH+1
I1=(NUM+NSEL)*2
DO 670 I=NUM+NSEL+NUM+1,I1
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
670 J=J+1
ELSE
C Añado machos no seleccionados.
DO 680 I=1,NUM
YY2(I)=YY(I)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(I)
BP2(I)=BP(I)
BTP2(I)=BTP(I)
680 UP2(I)=UP(I)
C Añado machos seleccionados.
J=NBM+1
DO 690 I=NUM+1,NUM+NSEL
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
690 J=J+1
C Añado hembras no seleccionadas.
J=NBM+NSEL+1
DO 710 I=NUM+NSEL+1,NUM+NSEL+NUM
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
710 J=J+1
C Añado hembras seleccionadas.
J=NBM+NSEL+NBH+1
I1=(NUM+NSEL)*2
DO 720 I=NUM+NSEL+NUM+1,I1
YY2(I)=YY(J)
NPAT2(I,1)=0
NPAT2(I,2)=0
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
720 J=J+1
ENDIF
C GENERACION FILIAL.
C Añado los seleccionados.
IP=NBM-NUM
IM=NBM+NSEL+NBH-(NUM+NSEL+NUM)
J=NGB+1
DO 730 I=I1+1,I1+NSH
YY2(I)=YY(J)
NPAT2(I,1)=NPAT(J,1)-IP
NPAT2(I,2)=NPAT(J,2)-IM
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
730 J=J+1
C Calculo los faltantes.
SUM=0.
DO 255 I=NGB+1,NRM
RM=BP(I)+BTP(I)+UP(I)
255 SUM=SUM+PNOR((T1-RM)/DVUE)
SUM=SUM/DBLE(NRM-NGB)
NUF=INT(DBLE(NSH)*SUM/(1.-SUM))
C IF (NUF.GT.3000) NUF=3000
IF (NUF.GT.20000) THEN
NUF=20000
NNF=NNF+1
ENDIF
C Añado los no seleccionados.
IF (NUF.GT.N1H) THEN
J=NGB+NSH+1
DO 740 I=I1+NSH+1,I1+NSH+N1H
YY2(I)=YY(J)
NPAT2(I,1)=NPAT(J,1)-IP
NPAT2(I,2)=NPAT(J,2)-IM
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
740 J=J+1
R2=(VUA/2.)**.5
DO 257 I=I1+NSH+N1H+1,I1+NSH+NUF
2222 CALL RANDOM(X2,Z)
NPAT2(I,1)=NUM+INT(Z*DBLE(NSEL))+1
CALL RANDOM(X2,Z)
NPAT2(I,2)=NUM+NSEL+NUM+INT(Z*DBLE(NSEL))+1
2224 CALL RANDOM(X2,UNI)
K=INT(UNI*100.)+1
IF (IS(K).EQ.0) GO TO 2224
IGP2(I)=K
BP2(I)=S(K)
BTP2(I)=TP2
CALL NORMAL(Z,X2)
UP2(I)=Z*R2+(UP2(NPAT2(I,1))+UP2(NPAT2(I,2)))/2.
CALL NORMAL(Z,X2)
YY2(I)=BP2(I)+BTP2(I)+UP2(I)+Z*DVUE
IF (YY2(I).GT.T1) GO TO 2222
257 CONTINUE
ELSE
J=NGB+NSH+1
DO 750 I=I1+NSH+1,I1+NSH+NUF
YY2(I)=YY(J)
NPAT2(I,1)=NPAT(J,1)-IP
NPAT2(I,2)=NPAT(J,2)-IM
IGP2(I)=IGP(J)
BP2(I)=BP(J)
BTP2(I)=BTP(J)
UP2(I)=UP(J)
750 J=J+1
ENDIF
NBM=NUM
NBH=NUM
NGB=(NUM+NSEL)*2
N1H=NUF
NG1=NSH+N1H
NRM=NGB+NG1
IF (NRM.GT.40000) THEN
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
WRITE(1,*)
WRITE(1,*) 'NRM>40000'
GO TO 10000
ENDIF
DO 760 I=1,NRM
YY(I)=YY2(I)
NPAT(I,1)=NPAT2(I,1)
NPAT(I,2)=NPAT2(I,2)
IGP(I)=IGP2(I)
BP(I)=BP2(I)
BTP(I)=BTP2(I)
760 UP(I)=UP2(I)
DO 580 I=1,100
SS(I)=SS(I)+S(I)
580 IS(I)=0.
DO 585 I=1,NRM
585 IS(IGP(I))=IS(IGP(I))+1
2345 NUM=500
NUF=1500
B1=B1+S(1)
B2=B2+RMP
B3=B3+VALE
BT1=BT1+TP1
BT2=BT2+TP2
T=T+T1
VA=VA+VUA
VER=VER+VUE
SNRM=SNRM+DBLE(NRM)
SNUM=SNUM+DBLE(NUM)
SNUF=SNUF+DBLE(NUF)
IL1=IL1+L1
IL2=IL2+L2
IL3=IL3+L3
C SRG1=SRG1+RG1
C SRG2=SRG2+RG2
C SRG3=SRG3+RG3
C SRLH1=SRLH1+RLH1
C SRLH2=SRLH2+RLH2
C SRLH3=SRLH3+RLH3
C SRS1=SRS1+RRS1
C SRS2=SRS2+RRS2
C SRS3=SRS3+RRS3
NEN=INT(DBLE(KL)/5000.)*5000
IF ((KL.EQ.1).OR.(KL.EQ.NEN)) THEN
CKL=DBLE(KL)
RB1=B1/CKL
RB2=B2/CKL
RB3=B3/CKL
RBT1=BT1/CKL
RBT2=BT2/CKL
RT=T/CKL
RVA=VA/CKL
RVER=VER/CKL
INRM=INT(SNRM/CKL)
INUM=INT(SNUM/CKL)
INUF=INT(SNUF/CKL)
L1=IL1/INT(CKL)
L2=IL2/INT(CKL)
L3=IL3/INT(CKL)
C RG1=SRG1/CKL
C RG2=SRG2/CKL
C RG3=SRG3/CKL
C RLH1=SRLH1/CKL
C RLH2=SRLH2/CKL
C RLH3=SRLH3/CKL
C RRS1=SRS1/CKL
C RRS2=SRS2/CKL
C RRS3=SRS3/CKL
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
WRITE(1,'(I3,4I6,6F8.2)') IJK,KL,INRM,INUM,INUF,RB1,RB2,
1RB3,RVA,RVER,RT
C WRITE(1,'(3I5,9F7.3)') L1,L2,L3,RG1,RG2,RG3,RLH1,RLH2,RLH3,
C 1RRS1,RRS2,RRS3
C WRITE(1,'(2I5,2F8.2)') NNU,NNF,RBT1,RBT2
C PRINT '(I3,4I6,6F8.2)',IJK,KL,INRM,INUM,INUF,RB1,RB2,
C 1RB3,RVA,RVER,RT
C WRITE(1,'(10I7)') (IS(I),I=1,20)
C WRITE(1,*)
C WRITE(1,'(10F8.4)') (RLS(I),I=1,20)
C WRITE(1,*)
C WRITE(1,'(10F7.1)') (S(I),I=1,20)
C WRITE(1,*)
C WRITE(1,'(10F7.1)') (SS(I)/CKL,I=1,20)
C WRITE(1,*)
C WRITE(1,*)
CLOSE(1)
C NNU=0
C NNF=0
ENDIF
IF (KL.LT.NESLA) GO TO 2000
C OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
C WRITE(1,'(20I4)') (IS(I),I=1,100)
C WRITE(1,'(10F7.1)') (S(I),I=1,20)
C WRITE(1,'(2I10)') NNU,NNF
C WRITE(1,*)
C CLOSE(1)
SB1=SB1+RB1
SB2=SB2+RB2
SB3=SB3+RB3
SBT1=SBT1+RBT1
SBT2=SBT2+RBT2
SVA=SVA+RVA
SVER=SVER+RVER
SINR=SINR+DBLE(INRM)
S2B1=S2B1+RB1**2
S2B2=S2B2+RB2**2
S2B3=S2B3+RB3**2
S2VA=S2VA+RVA**2
S2VER=S2VER+RVER**2
RIJK=DBLE(IJK)
NEN=INT(DBLE(IJK)/10.)*10
IF (IJK.EQ.NEN) THEN
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
WRITE(1,*)
WRITE(1,'(7F10.3)') SB1/RIJK,SB2/RIJK,SB3/RIJK,SVA/RIJK,
1SVER/RIJK,SBT2/RIJK,SINR/RIJK
WRITE(1,*)
CLOSE(1)
C PRINT*
C PRINT '(7F10.3)',SB1/RIJK,SB2/RIJK,SB3/RIJK,SVA/RIJK,
C 1SVER/RIJK,SBT2/RIJK,SINR/RIJK
C PRINT*
ENDIF
1000 CONTINUE
RNIT=DBLE(IREPE)
C VARB1=(S2B1-SB1*SB1/RNIT)/(RNIT-1.)
C VARB2=(S2B2-SB2*SB2/RNIT)/(RNIT-1.)
C VARB3=(S2B3-SB3*SB3/RNIT)/(RNIT-1.)
C VARVA=(S2VA-SVA*SVA/RNIT)/(RNIT-1.)
C VARVER=(S2VER-SVER*SVER/RNIT)/(RNIT-1.)
SB1=SB1/RNIT
SB2=SB2/RNIT
SB3=SB3/RNIT
SBT1=SBT1/RNIT
SBT2=SBT2/RNIT
INR=INT(SINR/RNIT)
SVA=SVA/RNIT
SVER=SVER/RNIT
ST=ST/RNIT
C RMB1=VARB1+(B(1)-SB1)**2
C RMB2=VARB2+(B(50)-SB2)**2
C RMB3=VARB3+(B(100)-SB3)**2
C RMVA=VARVA+(VU-SVA)**2
C RMVER=VARVER+(VE-SVER)**2
C SP1=SP1/RNIT
C SP2=SP2/RNIT
C SP3=SP3/RNIT
C SH1=SH1/RNIT
C SH2=SH2/RNIT
C SH3=SH3/RNIT
OPEN(1,FILE='HRESA100',STATUS='OLD',ACCESS='APPEND')
WRITE(1,*)
WRITE(1,'(6F8.2,I6)') SB1,SB2,SB3,SVA,SVER,SBT2,INR
C WRITE(1,'(5F8.2)') RMB1,RMB2,RMB3,RMVA,RMVER
WRITE(1,*)
C WRITE(1,'(3F8.3)') SP1,SP2,SP3
C WRITE(1,'(3F8.3)') SH1,SH2,SH3
CLOSE(1)
10000 STOP
END
SUBROUTINE RANDOM(A1,R)
IMPLICIT REAL*8 (A-H,O-Z)
FF=2147483647.
B=16807.*A1
L=INT(B/FF)
A1=B-DBLE(L)*FF
R=A1/FF
RETURN
END
SUBROUTINE NORMAL(Z,P1)
IMPLICIT REAL*8 (A-H,O-Z)
CALL RANDOM(P1,U1)
CALL RANDOM(P1,U2)
Z=((-2.*LOG(U1))**0.5)*COS(2.*3.1416*U2)
RETURN
END
SUBROUTINE PIKSRT(N,ARR)
C Ordena el vector ARR de dimensión N mediante la estrategia de
C "Straight insertion" (Numerical Recipes, pag. 227).
REAL*8 ARR(100)
DO 12 J=2,N
A=ARR(J)
DO 11 I=J-1,1,-1
IF (ARR(I).LE.A) GO TO 10
ARR(I+1)=ARR(I)
11 CONTINUE
I=0
10 ARR(I+1)=A
12 CONTINUE
RETURN
END
SUBROUTINE SORT2(N,IRA,RIA)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 IRA(200000,2)
REAL*8 RIA(200000)
L=N/2+1
IR=N
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
IRRA=IRA(L,1)
IRRB=IRA(L,2)
RRIA=RIA(L)
ELSE
IRRA=IRA(IR,1)
IRRB=IRA(IR,2)
RRIA=RIA(IR)
IRA(IR,1)=IRA(1,1)
IRA(IR,2)=IRA(1,2)
RIA(IR)=RIA(1)
IR=IR-1
IF (IR.EQ.1) THEN
IRA(1,1)=IRRA
IRA(1,2)=IRRB
RIA(1)=RRIA
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF (J.LE.IR) THEN
IF (J.LT.IR) THEN
IF (IRA(J,1).LT.IRA(J+1,1)) J=J+1
ENDIF
IF (IRRA.LT.IRA(J,1))THEN
IRA(I,1)=IRA(J,1)
IRA(I,2)=IRA(J,2)
RIA(I)=RIA(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
IRA(I,1)=IRRA
IRA(I,2)=IRRB
RIA(I)=RRIA
GO TO 10
END
SUBROUTINE SORT(N,IRA,IMG,UA,RIA)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 IRA(40000),IMG(40000)
REAL*8 UA(40000),RIA(40000)
L=N/2+1
IR=N
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
IRRA=IRA(L)
IMGA=IMG(L)
RRB=UA(L)
RRIA=RIA(L)
ELSE
IRRA=IRA(IR)
IMGA=IMG(IR)
RRB=UA(IR)
RRIA=RIA(IR)
IRA(IR)=IRA(1)
IMG(IR)=IMG(1)
UA(IR)=UA(1)
RIA(IR)=RIA(1)
IR=IR-1
IF (IR.EQ.1) THEN
IRA(1)=IRRA
IMG(1)=IMGA
UA(1)=RRB
RIA(1)=RRIA
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF (J.LE.IR) THEN
IF (J.LT.IR) THEN
IF (RIA(J).LT.RIA(J+1)) J=J+1
ENDIF
IF (RRIA.LT.RIA(J))THEN
IRA(I)=IRA(J)
IMG(I)=IMG(J)
UA(I)=UA(J)
RIA(I)=RIA(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
IRA(I)=IRRA
IMG(I)=IMGA
UA(I)=RRB
RIA(I)=RRIA
GO TO 10
END
SUBROUTINE SORT3(N,IUA,RIA)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 IUA(100)
REAL*8 RIA(100)
L=N/2+1
IR=N
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
IRRB=IUA(L)
RRIA=RIA(L)
ELSE
IRRB=IUA(IR)
RRIA=RIA(IR)
IUA(IR)=IUA(1)
RIA(IR)=RIA(1)
IR=IR-1
IF (IR.EQ.1) THEN
IUA(1)=IRRB
RIA(1)=RRIA
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF (J.LE.IR) THEN
IF (J.LT.IR) THEN
IF (RIA(J).LT.RIA(J+1)) J=J+1
ENDIF
IF (RRIA.LT.RIA(J))THEN
IUA(I)=IUA(J)
RIA(I)=RIA(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
IUA(I)=IRRB
RIA(I)=RRIA
GO TO 10
END
SUBROUTINE SORTI(N,NPA,IRA,IMG,UA,RIA)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER*4 IRA(40000),NPA(40000,2),IMG(40000)
REAL*8 UA(40000),RIA(40000)
L=N/2+1
IR=N
10 CONTINUE
IF (L.GT.1) THEN
L=L-1
IRRA=IRA(L)
IP1=NPA(L,1)
IP2=NPA(L,2)
IMGA=IMG(L)
RRB=UA(L)
RRIA=RIA(L)
ELSE
IRRA=IRA(IR)
IP1=NPA(IR,1)
IP2=NPA(IR,2)
IMGA=IMG(IR)
RRB=UA(IR)
RRIA=RIA(IR)
IRA(IR)=IRA(1)
NPA(IR,1)=NPA(1,1)
NPA(IR,2)=NPA(1,2)
IMG(IR)=IMG(1)
UA(IR)=UA(1)
RIA(IR)=RIA(1)
IR=IR-1
IF (IR.EQ.1) THEN
IRA(1)=IRRA
NPA(1,1)=IP1
NPA(1,2)=IP2
IMG(1)=IMGA
UA(1)=RRB
RIA(1)=RRIA
RETURN
ENDIF
ENDIF
I=L
J=L+L
20 IF (J.LE.IR) THEN
IF (J.LT.IR) THEN
IF (RIA(J).LT.RIA(J+1)) J=J+1
ENDIF
IF (RRIA.LT.RIA(J))THEN
IRA(I)=IRA(J)
NPA(I,1)=NPA(J,1)
NPA(I,2)=NPA(J,2)
IMG(I)=IMG(J)
UA(I)=UA(J)
RIA(I)=RIA(J)
I=J
J=J+J
ELSE
J=IR+1
ENDIF
GO TO 20
ENDIF
IRA(I)=IRRA
NPA(I,1)=IP1
NPA(I,2)=IP2
IMG(I)=IMGA
UA(I)=RRB
RIA(I)=RRIA
GO TO 10
END
SUBROUTINE GENO(IP,IM,II,X2)
IMPLICIT REAL*8 (A-H,O-Z)
IF ((IP.EQ.1).AND.(IM.EQ.1)) THEN
II=1
ELSE IF ((IP.EQ.2).AND.(IM.EQ.2)) THEN
CALL RANDOM(X2,Z)
IF (Z.GT..75) THEN
II=3
ELSE IF (Z.GE..25) THEN
II=2
ELSE
II=1
ENDIF
ELSE IF ((IP.EQ.3).AND.(IM.EQ.3)) THEN
II=3
ELSE IF ((IP.EQ.1).AND.(IM.EQ.2)) THEN
CALL RANDOM(X2,Z)
IF (Z.GT..5) THEN
II=2
ELSE
II=1
ENDIF
ELSE IF ((IP.EQ.2).AND.(IM.EQ.1)) THEN
CALL RANDOM(X2,Z)
IF (Z.GT..5) THEN
II=2
ELSE
II=1
ENDIF
ELSE IF ((IP.EQ.1).AND.(IM.EQ.3)) THEN
II=2
ELSE IF ((IP.EQ.3).AND.(IM.EQ.1)) THEN
II=2
ELSE IF ((IP.EQ.2).AND.(IM.EQ.3)) THEN
CALL RANDOM(X2,Z)
IF (Z.GT..5) THEN
II=3
ELSE
II=2
ENDIF
ELSE IF ((IP.EQ.3).AND.(IM.EQ.2)) THEN
CALL RANDOM(X2,Z)
IF (Z.GT..5) THEN
II=3
ELSE
II=2
ENDIF
ENDIF
END
SUBROUTINE GENO2(IP,IM,IH,PRO,X2)
IMPLICIT REAL*8 (A-H,O-Z)
IF (IH.EQ.1) THEN
IF ((IP.EQ.1).AND.(IM.EQ.1)) THEN
PRO=1.
ELSE IF ((IP.EQ.2).AND.(IM.EQ.2)) THEN
PRO=.25
ELSE IF ((IP.EQ.2).AND.(IM.EQ.1)) THEN
PRO=.5
ELSE IF ((IP.EQ.1).AND.(IM.EQ.2)) THEN
PRO=.5
ELSE
PRO=0.
ENDIF
ELSE IF (IH.EQ.2) THEN
IF ((IP.EQ.2).AND.(IM.EQ.2)) THEN
PRO=.5
ELSE IF ((IP.EQ.2).AND.(IM.EQ.1)) THEN
PRO=.5
ELSE IF ((IP.EQ.1).AND.(IM.EQ.2)) THEN
PRO=.5
ELSE IF ((IP.EQ.2).AND.(IM.EQ.3)) THEN
PRO=.5
ELSE IF ((IP.EQ.3).AND.(IM.EQ.2)) THEN
PRO=.5
ELSE IF ((IP.EQ.1).AND.(IM.EQ.3)) THEN
PRO=1.
ELSE IF ((IP.EQ.3).AND.(IM.EQ.1)) THEN
PRO=1.
ELSE
PRO=0.
ENDIF
ELSE IF (IH.EQ.3) THEN
IF ((IP.EQ.3).AND.(IM.EQ.3)) THEN
PRO=1.
ELSE IF ((IP.EQ.2).AND.(IM.EQ.2)) THEN
PRO=.25
ELSE IF ((IP.EQ.2).AND.(IM.EQ.3)) THEN
PRO=.5
ELSE IF ((IP.EQ.3).AND.(IM.EQ.2)) THEN
PRO=.5
ELSE
PRO=0.
ENDIF
ENDIF
END
FUNCTION DEN(Z)
C CALCULA LA ORDENADA N(0,1) DE LA ABCISA Z
IMPLICIT REAL*8 (A-H,O-Z)
DEN=.39894228*EXP(-.5*Z**2)
RETURN
END
FUNCTION PNOR(X)
IMPLICIT REAL*8 (A-H,O-Z)
C CALCULA LA INCIDENCIA (P) PARA CUALQUIER ABCISA (X)
IF(X.LT.0.) GOTO 1
PNOR=PNORMD(X)
GOTO 2
1 XX=-X
PNOR=1.-PNORMD(XX)
2 RETURN
END
FUNCTION PNORMD(X)
C ESTA FUNCION NO ESTA BIEN DEFINIDA PARA ABCISAS NEGATIVAS.
C POR ESTA RAZON SE CALCULA LA PROBABILIDAD COMPLEMENTARIA DEL
C VALOR NEGATIVO TAL Y COMO SE DEFINE EN PNOR( ).
IMPLICIT REAL*8 (A-H,O-Z)
T=1./(1.+.2316419*X)
Z=.39894228*EXP(-.5*X*X)
PNORMD=1.-Z*(T*(.31938153-T*(.356563782-T*(1.781477937-T*
1 (1.821255978-T*1.330274429)))))
RETURN
END
FUNCTION XNOR(PROB)
IMPLICIT REAL*8 (A-H,O-Z)
C CALCULA LA ABCISA (X) PARA CUALQUIER INCIDENCIA (P)
IF(PROB.LT..5) THEN
XNOR=-XNORMD(PROB)
ELSE IF(PROB.GT..5)THEN
XNOR=XNORMD(1.-PROB)
ELSE
XNOR=0.
ENDIF
RETURN
END
FUNCTION XNORMD(PROB)
C SOLO DEFINIDA PARA PROB EN EL ESPACIO 0-.5
IMPLICIT REAL*8 (A-H,O-Z)
T=SQRT(LOG(1./(PROB*PROB)))
XNORMD=T-(2.515517+T*(.802853+T*.010328))/(1.+T*
1 (1.432788+T*(.189269+T*.001308)))
RETURN
END
EOF
echo $@
PARAMS=$@
export PARAMS
PATH=/bin:/usr/bin:/usr/local/bin
export PATH
LD_LIBRARY_PATH=/lib:/usr/lib:/usr/local/lib:/lib64
export LD_LIBRARY_PATH
rm -f foo
rm -f foo2
/usr/bin/f77 -o simulation simulation.f >foo &> foo2
cat foo
cat foo2
rm -f foo
rm -f foo2
echo $status
echo made
echo params $PARAMS $@
./simulation $PARAMS >foo &> foo2
cat foo
cat foo2
cat HRESA100
rm HRESA100
See more files for this project here