message terminal Stop error

**progrmmation** · 09/03/2023, 15h16

Bonjour tous le monde
quand j'augmente l'intervalle d'un paramètre dans le fichier.dat, le programme ne tourne pas il s'arrete ici et il m'affiche ce message '' STOP error''
Sachant que ce programme est constitué de plusieurs fichiers.FOR liés entre eux et que tous les fichiers tournes correctement sauf celui là
Est ce que vous avez une solution et merci d'avance

| 1
Warning: Fortran 2018 deleted feature: Shared DO termination label 70 at (1)
donner votre choix entre 0, 2, 4, 6, 8, 9 :
9
donner la masse :
103
=>Le parametre d inertie h2/2J ~ fmi = E2+/6 => 4.3791857911118497E-002
Saisir le parametre de Coriolis ico :
0
Voulez vous saisir de fmi exp 0(Non) 1(Oui)
0
STOP error

**vmagnin** · 09/03/2023, 19h22

Bonjour @progrmmation

sans voir le code source, il est difficile de donner une réponse.

Il faudrait tout d'abord que vous nous donniez la commande de compilation utilisée (est-ce que vous utilisez GFortran ?), pour voir si on peut ajouter des options de diagnostic qui aiderait à détecter le problème.

Je vois un warning au début :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Warning: Fortran 2018 deleted feature: Shared DO termination label 70 at (1)

qui montre, ainsi que les extensions .FOR, qu'il s'agit d'un code dans une (très) vieille norme Fortran (66 ou 77). Mais ça n'est probablement pas la cause de l'erreur.

Le message "STOP error" semble être affiché par le programme lui-même, peut-être pour indiquer qu'un paramètre entré est invalide. Est-ce que vous voyez dans le code source une instruction PRINT ou WRITE qui écrirait ce message ?

**progrmmation** · 10/03/2023, 09h27

Bonjour voilà le code source il est trop long
IMPLICIT REAL*8 (A-H,O-Z) ELE00060
INTEGER OMG1,omg2,om2,PRT1,prt2,prt,prtt,omt,omt4,prtt4,rtt,
*om1,pt1,omtt1,omtt2,ptt1,ptt2
parameter(ndim1=200)
c DIMENSION OMEGA1(NDIM1),PART1(NDIM1),NOMT1(NDIM1)
DIMENSION NNq(20,10),NZq(20,10),LAq(20,10),Bq(20,10)
DIMENSION NN(20,10),NZ(20,10),LA(20,10),B(20,10)
Dimension omg1(ndim1),prt1(ndim1),nmt1(ndim1)
Dimension omg2(ndim1),prt2(ndim1),nmt2(ndim1)
Dimension eq1(ndim1),eq2(ndim1),u(20),EQq(20),v(20),NCF(20),
*prt(20),nmt(20),om2(20),nc(20),bj(20)
Dimension zz(ndim1),X(ndim1),numer(100),vi(100),t11(ndim1)
dimension NIV(2),NDIAG(2),NDIAGP(2),NDIAG1(2),NDIAG2(2)
dimension A(ndim1),omt(ndim1),prtt(ndim1),nmtt(ndim1),r(ndim1),
*valp(60),vectp(20,20),q2f(100),q2g(100),eqp1(100),h(60),aa(100)
dimension omt4(60),prtt4(60),nmtt4(60),nct4(60),nctt4(60),
*nct(60),nctt(60),mt(60),rtt(60),mtt(60),exy(60),vik(20,20)
dimension t31(20,20),ab(100),ga(2),t0(20),ra(100)
dimension omtt1(100),omtt2(100),ptt1(100),ptt2(100),nmtt1(100),
*nmtt2(100),xf(100),xg(100)
dimension xj0(10,10),xj1(10,10),xj2(10,10),xj3(10,10),
*yj0(10,10),yj1(10,10),xc0(20,10,10),xc1(20,10,10),xc2(20,10,10),
*xc3(20,10,10),vo(20,20),co(20,20),yj2(10,10),yj3(10,10),
*yjj3(10,10)
dimension om1(20),pt1(20),nt1(20),nc1(20),nct1(20),ver1(20)
character*1 plus,moins,blanc,etoile,pop,indic
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/ BIB04570
plus='+'
moins='-'
blanc=' '
etoile='*'
ND2=20
ir=11
iw=13
ik=19
open(1,status='old',file='ap.noy')
open(2,status='old',file='an.noy')
open(8,status='new',file='cpdg.res',form='formatted')
open(iw,status='new',file='cpc.par')
open(ik,status='new',file='cpc2.par')
open(9,status='old',file='tda.par')
open(3,status='old',file='rep')
open(4,status='old',file='ren')
open(11,status='new',file='cor')
open(31,status='new',file='dig')
c ***********************************************************************
write(*,*)'donner votre choix entre 0, 2, 4, 6, 8, 9 :'
read(*,*)ichoix
c write(*,*)'donner le nombre de masse du noyau en question :'
write(*,*)'donner la masse :'
read(*,*)mas
do 8 l=1,2
read(9,*) LP
read(9,*) NIV(LP),NDIAG(LP),NDIAGP(LP),NDIAG1(LP),NDIAG2(LP)
8 continue
read(9,*) W,CHI2,RR
write(8,*)' LES PARAMETRES DU PROGRAMME - TDA - '
write(8,5001)w,chi2

write(8,*)
sig=0.0d0
sip=0.0d0
sim=0.0d0
siv1=0.0d0
siv2=0.0d0
siw1=0.0d0
siw2=0.0d0
fdm=0.0d0
e2p=0.0d0
qsj=0.0d0
qgs=0.0d0

cccccccc
do 10 lt=1,2
DO 5 i=1,ndiag(lt)
read(9,*)omg1(i),prt1(i),nmt1(i)
read(9,*)omg2(i),prt2(i),nmt2(i)
read(9,*)eq1(i),eq2(i),eqp1(i)
read(9,*)q2f(i),q2g(i),X(i)
5 CONTINUE
c
read(lt,*)XKAPPA,XMU,eps
read(lt,*)JMIN,JMAX,NPART,NPROT,NNEUTR,ICHOI1,ICHOI2,INORM,IMEL, BIB03960
*DELT,DELTA0,DELTA2,Ga(lt),G2,CG2G0 BIB03970
c write(*,*)'donner le valeur de deformation'
c
c read(*,*)eps
ILS=3
il2=3
eps2=eps*eps
za=7.0d0/3.0d0
zif=float(mas)**za
tif=zif*eps2
E2=1176.0d0/tif
fmi=e2/6.0d0
write(*,*)'=>Le parametre d inertie h2/2J ~ fmi = E2+/6 => ',fmi
if(lt.eq.2) go to 97
write(8,*)'ccccccccccccccccc PROTON cccccccccccccccccccc'
go to 93
97 write(8,*)'cccccccccccccccccc NEUTRON ccccccccccccccccccc'
93 write(8,*)
write(8,*)' LES PARAMETRES DU PROGRAMME - NILSSON -'
write(8,5000)xkappa,xmu,eps
write(8,*)
write(8,*)' LES PARAMETRES DU PROGRAMME - BCS -'
write(8,5002)nprot,nneutr,npart
write(8,5003)DELT,g0
write(8,*)
j=0
1 j=j+1
READ(lt,*)ltt
if(ltt.ne.lt) go to 6
read(lt,*)OM2(J),PRT(J),NMT(J),NCF(J)
kmax=ncf(j)
read(lt,*)u(J),EQq(j),v(j)
read(lt,*)(NNq(J,K),NZq(J,K),LAq(J,K),Bq(J,K),K=1,KMAX)
go to 1
6 continue
nq=j-1
ccccccccccc

IF(LT.eq.1) in=3
if(lt.eq.2) in=4
j1=0
299 j1=j1+1
read(in,*)ip
if(ip.eq.2) go to 3
read(in,*)omtt1(j1),ptt1(j1),nmtt1(j1),omtt2(j1),ptt2(j1),
*nmtt2(j1),xf(j1),xg(j1)
go to 299
3 continue
read(in,*)xtot,pgs
write(*,*)'Saisir le parametre de Coriolis ico :'
read(*,*)ico
if(ico.eq.0) go to 289
fmi=1.0d0/xtot
write(*,*) '=> Le fmi nouvellement calcule de 1.0/xtot : ',fmi
289 continue
mq=j1-1
if(lt.eq.2) go to 277
fdm=fdm+pgs
e2p=e2p+w+2.0d0*fmi
277 continue
ccccccccccccccccccccccccc
c write(*,*)'zzzzzzzzz'
cccccccccccccccccc
c if(lt.ne.2) go to 99
write(8,*)' ORTHONORMALISATION DE LA BASE DU COUPLAGE QUASIPARTICU
*LE-PHONON'
write(8,*)
write(8,*)' L ETAT DU COUPLAGE DOMINANTE EST INDIQUE PAR ETOILE *'
c write(*,*)'pop',pop
write(*,*)'Voulez vous saisir de fmi exp 0(Non) 1(Oui)'
read(*,*)iwob
if(iwob.eq.0) go to 99
write(*,*)'Saisir la valeur experimentale de E2+/6'
read(*,*)fmi
write(*,*)'fmi = E2+/6 experimentale saisie : ',fmi
99 continue
c write(*,*)'eeee'
cccccccccccccccccccccccccccccccccccc
c write(*,*)'entree bat2'
call bat2(X,ndiag,lt,nq,mq,om2,omg1,omg2,prt,prt1,prt2,nmt,
*nmt1,nmt2,omtt1,ptt1,nmtt1,omtt2,ptt2,nmtt2,xf,xg,ndim1,w,chi2,
*rr,q2g,q2f,u,v,ga,sif,sig,sip,sim,siv1,siv2,siw1,siw2,fmi,qgs,
*gsj)
c write(*,*)'sortie bat2'
write(*,*)siv1,siv2,siw1,siw2,sim,sip,sig,sif
cccccccccccccccccccc
if(lt.ne.2) go to 288
e2p2=siv1+siv2+sim+siw1+siw2+sip+sig
e2p3=siw1+siw2
ep3=e2p+e2p3
e2p1=e2p+e2p2
ep6=e2p+siv1+siv2
fdm=fdm+qgs+gsj
c write(*,*)'fdm',fdm,'pgs',pgs
c write(*,*)'gsj',gsj,'qgs',qgs
write(*,*)'w',w,'e2p',e2p1,ep3,ep6
288 continue
pop=plus
2 continue

3333 do 20 no=1,15,2
no1=abs(4-no)
no2=4+no
nox=no1
noy=no2
do 1119 j=1,30

a(j)=0.0d0

1119 continue
c write(*,*)'entree bat1',pop
call bat1(X,pop,noy,nox,no,ndiag,lt,nq,om2,omg1,omg2,
*prt,prt1,prt2,nmt,nmt1,nmt2,omt,prtt,nmtt,nct,nctt,omtt1,
*ptt1,nmtt1,omtt2,ptt2,nmtt2,xf,xg,mq,a,eq1,eq2,eqq,ndim1,n,nm,
*t11,w,chi2,rr,q2g,q2f,mt,rtt,mtt,t31,m2,m4,exy,sif,sig,sim,sip,
*siv1,siv2,siw1,siw2,u,v,ga,ichoix,fmi,t0,m2n,xj0,xj1,xj2,xj3,
*yj0,yj1,yj2,yj3,kk,k2,n1,n2,fdm)

c write(*,*)'sortie bat1'
ccccccccccccccc
36 if(lt.ne.2) go to 35
write(8,*)'----------------------------------------------------'
do 1111 jm=1,30
c if(a(jm).eq.0.0d0) go to 1111
write(31,*)jm,a(jm)
1111 continue

35 ndim=n
c write(*,*)'entree EIGEN'
CALL EIGEN(a,R,NDIM,0)
c write(*,*)'sortie EIGEN'
BIB02090
C RECHERCHE DES ETATS DONT L ENERGIE,POUR LA DEFORMATION EPSMIN, BIB02110
C EST COMPRISE ENTRE LES VALEURS EMIN ET EMAX BIB02120
c IF(NEPS.NE.1)GO TO 810 BIB02140
IIMIN=1 BIB02150
IIMAX=0 BIB02160
Emax=80.0d0
emin=-80.0d0
DO 800 IV=1,NDIM BIB02170
IR=IV*(IV+1)/2 BIB02180
IF(A(IR).GT.EMAX)IIMIN=IV+1 BIB02190
IF(A(IR).GT.EMIN)IIMAX=IV BIB02200
800 CONTINUE
c write(*,*)'boucle 800' BIB02210
NINTER=IIMAX+1-IIMIN BIB02220

C IF(NINTER.LT.1)GO TO 10 BIB02230
C IF(NINTER.LE.ND2)GO TO 802 BIB02240
C WRITE(IW,5002)NINTER,ND2 BIB02250
C STOP BIB02260
802 CONTINUE BIB02270
C DO 805 IE=1,NINTER BIB02280
C NMTt(IE)=NMTt(IE+IIMIN-1) BIB02290
C 805 CONTINUE BIB02300
810 CONTINUE BIB02310
DO 910 IV=IIMIN,IIMAX BIB02370
IE=IV+1-IIMIN BIB02380
IR=IV*(IV+1)/2 BIB02390
VALP(IE)=A(IR) BIB02400
JA=(IV-1)*NDIM+1 BIB02410
JB=JA+NDIM BIB02420
DO 907 JV=JA,JB BIB02430
IC=JV+1-JA BIB02440
VECTP(IC,IE)=R(JV) BIB02450
c
907 CONTINUE BIB02460
c
910 CONTINUE
c write(*,*) 'boucle 910' BIB02470
C DO 934 IE=1,NINTER BIB02620
C ICS=IE+IIMIN-1 BIB02640
C IF(VECTP(ICS,IE).GT.0.0D0)GO TO 934 BIB02650
C DO 933 IC=1,NDIM BIB02660
C VECTP(IC,IE)=-VECTP(IC,IE)
C 933 CONTINUE BIB02680
C 934 CONTINUE BIB02690
ng=0
do 666 k=1,m2
ng=ng+1
omt4(ng)=mt(k)
prtt4(ng)=rtt(k)
nmtt4(ng)=mtt(k)
nct4(ng)=0
nctt4(ng)=mt(k)
666 continue
c write(*,*) 'boucle 666'
DO 200 IE=1,NINTER BIB04650
if(lt.ne.2) go to 777
WRITE(31,6001)VALP(IE)

WRITE(31,9020) BIB04670

777 continue
id=ie+iimin-1
if(ie.ne.1)go to 66
kmax=nc(1)
66 continue
c write(iw,*)lt
DO 202 IC=1,NDIM BIB04720
INDIC=BLANC BIB04770
IF((IE+IIMIN-1).EQ.IC)INDIC=ETOILE BIB04780
c PAR=part(ic)
if(lt.ne.2) go to 778
vp=vectp(ic,ie)/sqrt(valp(ie))
WRITE(31,6002)indic,nct(ic),omt(ic),prtt(ic),nmtt(ic),nctt(ic),
*prtt(ic),VECTP(IC,ie),vp
778 IF((IE+IIMIN-1).ne.IC) go to 202
ng=ng+1
omt4(ng)=omt(ic)
prtt4(ng)=prtt(ic)
nmtt4(ng)=nmtt(ic)
nct4(ng)=nct(ic)
nctt4(ng)=nctt(ic)
202 continue
c write(*,*) 'boucle 202'
200 CONTINUE BIB04820
300 CONTINUE BIB04830
c
c WRITE(8,9020)

c 19 continue
lpp=0
c write(iw,*)llp
c if(lt.ne.2) go to 518
c do 517 i=1,ninter

c do 517 j=1,ndim

c 517 continue
c stop
c 518 continue
cccccccccccccccccccccccccccccccccccc
nw=0
do 80 i=1,ninter
do 85 j=1,i

nm=0
do 90 k=1,ndim
do 95 l=1,ndim
nm=nm+1
Vi(nm)=(vectp(k,i)*vectp(l,j))/sqrt(valp(i)*valp(j))
95 continue
90 continue
c write(*,*) 'boucle 90'
qm=0.0d0
m1=0
do 60 i1=1,ndim
k=i1+(i1-1)*(ndim-1)
l=i1
do 65 j1=1,i1
m1=m1+1
if(j1.eq.i1) go to 70
qm=qm+(vi(l)+vi(k))*t11(m1)
k=k+1
l=l+ndim
go to 75
70 qm=qm+vi(k)*t11(m1)

75 continue
65 continue
60 continue
c write(*,*) 'boucle 60'
nw=nw+1
aa(nw)=qm
if(lt.ne.2) go to 516
c write(8,*)nw,'aa=',aa(nw),'t11=',t11(nw)
516 continue
85 continue
80 continue
cccccccccccccccccccccccccccccccc
do 190 i=1,ndim
do 191 j=1,ndim
vik(j,i)=vectp(j,i)/sqrt(valp(i))
c if(lt.eq.2.and.no.eq.1) write(8,*)'111',vik(j,i),j,i
191 continue
190 continue
c write(*,*) 'boucle 190'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
if(no.ne.1) go to 195
do 760 j1=1,ndim
do 760 j2=1,ndim
xe=0.0d0
do 761 i1=1,ndim
do 761 i2=1,ndim
xe=xe+vik(i1,j1)*vik(i2,j2)*yj3(i1,i2)
c if(j1.eq.1.and.j2.eq.2)write(8,*)'....yj3',yj3(i1,i2),
c *vik(i1,j1),vik(i2,j2),i1,i2
c if(j1.eq.2.and.j2.eq.1)write(8,*)'..///..yj3',yj3(i1,i2),
c *vik(i1,j1),vik(i2,j2),i1,i2

761 continue
yjj3(j1,j2)=xe
c write(8,*)'..1111..',j1,j2,xe
760 continue
c write(*,*) 'boucle 760'
cccccccccccccccccc
do 765 j=1,ndim
do 765 k=1,m2
xe1=0.0d0
xe2=0.0d0
do 767 i=1,ndim
xe1=xe1+vik(i,j)*yj1(i,k)
xe2=xe2+vik(i,j)*yj2(i,k)
c write(8,*)'....yj1-2',xe1,xe2,vik(i,j)

767 continue
yj1(j,k)=xe1
yj2(j,k)=xe2
765 continue
cccccccccccccccc
195 continue
c write(*,*) 'test 195'
ccccccccccccccccccccccccccccccccc
write(*,*)'ndim,k2,kk,no', ndim,k2,kk,no
do 700 j1=1,ndim
if(k2.eq.0) goto 7000
do 701 k=1,k2
xa=0.0d0
do 702 i=1,ndim
xa=xa+vik(i,j1)*xj1(i,k)
702 continue
c write(*,*) 'boucle 702'
c if(no.eq.0) goto 7000
xc1(no,k,j1)=xa
c write(8,*)'........xc1',xa
701 continue
7000 continue
c write(*,*) 'test 7000'
if(kk.eq.0) goto 7001
do 703 k=1,kk
xb=0.0d0
do 704 i=1,ndim
xb=xb+vik(i,j1)*xj3(i,k)
c write(8,*)'********',vik(i,j1),xj3(i,k)
704 continue
c write(*,*) 'boucle 704'
c if(no.eq.0) goto 7001
xc3(no,k,j1)=xb
c write(8,*)'........xc3',xb,kk
703 continue
c write(*,*) 'boucle 703'
7001 continue
700 continue
c write(*,*) 'boucle 700'
cccccccccccccccccccc
do 710 i1=1,m2
do 711 k=1,k2
xc0(no,k,i1)=xj0(i1,k)
c write(8,*)'........xc0',no,xj0(i1,k)
711 continue
do 712 k=1,kk
xc2(no,k,i1)=xj2(i1,k)
c write(8,*)'........xc2',xj2(i1,k)
712 continue
710 continue
cccccccccccccccccccccccccccc
IF(no.eq.1) go to 192
np=no-2
write(*,*)'no',no,'ndim',ndim,'kdj',kdj
c if(ndim.ne.kdj) stop 'error co1'
do 720 j2=1,ndim
do 721 i1=1,mj
xc=0.0d0
do 722 k=1,ndim
xc=xc+vik(k,j2)*xc2(np,k,i1)
722 continue
xc2(np,j2,i1)=xc
c write(8,*)'...8888.....xc2',xc
721 continue
do 723 j1=1,ndj
xd=0.0d0
do 724 k=1,ndim
xd=xd+vik(k,j2)*xc3(np,k,j1)
724 continue
xc3(np,j2,j1)=xd
c write(8,*)'...8888.....xc3',xd
723 continue
720 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cccccccccccccccccccccccccccccccc
192 n3=0
do 100 i=1,m2n
n3=n3+1
ab(n3)=t0(i)
c if(lt.ne.2) go to 94
c write(8,*)'ab=',ab(n3)
c 94 continue
100 continue
l3=0
do 105 j=1,ndim
do 199 i=1,m2
qt3=0.0d0
do 110 k=1,ndim
qt3=qt3+vik(k,j)*t31(i,k)
c write(13,*)'t31',t31(i,k),vik(k,j),i,k,j,qt3
110 continue
n3=n3+1
ab(n3)=qt3
if(lt.ne.2) go to 91
c write(8,*)'ab=',ab(n3)
91 continue
199 continue
do 115 j1=1,j
l3=l3+1
n3=n3+1
ab(n3)=aa(l3)
if(lt.ne.2) go to 92
c write(8,*)'ab=',ab(n3)
92 continue
115 continue
105 continue
ccccccccccccc
c write(8,*)'....................................'
ndi=ndim+m2
do 799 ij=1,n3
write(8,*)'....',ij,ab(ij)
799 continue
c write(*,*)'entree EIGEN2'
CALL EIGEN(ab,R,NDI,0)
c write(*,*)'sortie EIGEN2'
IIMIN=1
IIMAX=0
Emax=80.0d0
emin=-80.0d0
DO 890 IV=1,NDI BIB02170
IR=IV*(IV+1)/2 BIB02180
IF(Ab(IR).GT.EMAX)IIMIN=IV+1 BIB02190
IF(Ab(IR).GT.EMIN)IIMAX=IV BIB02200
890 CONTINUE BIB02210
NINTER=IIMAX+1-IIMIN BIB02220
892 CONTINUE BIB02270
c DO 895 IE=1,NINTER BIB02280
c NMTt(IE)=NMTt(IE+IIMIN-1) BIB02290
c 895 CONTINUE BIB02300
DO 990 IV=IIMIN,IIMAX BIB02370
IE=IV+1-IIMIN BIB02380
IR=IV*(IV+1)/2 BIB02390
VALP(IE)=Ab(IR) BIB02400
JA=(IV-1)*NDI+1
JB=JA+NDI BIB02420
DO 997 JV=JA,JB BIB02430
IC=JV+1-JA
c write(*,*)'997',JV
VECTP(IC,IE)=R(JV)
c
997 CONTINUE BIB02460
c
990 CONTINUE
C DO 994 IE=1,NINTER BIB02620
C ICS=IE+IIMIN-1 BIB02640
C IF(VECTP(ICS,IE).GT.0.0D0)GO TO 994 BIB02650
C DO 993 IC=1,NDI BIB02660
C VECTP(IC,IE)=-VECTP(IC,IE)
C 993 CONTINUE BIB02680
C 994 CONTINUE BIB02690
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
if(no.ne.1) go to 196
do 773 j1=1,ndi
do 770 j2=1,ndi
xn=0.0d0
xh=0.0d0
do 771 i1=1,m2
do 771 i2=1,m2
xn=xn+vectp(i1,j1)*vectp(i2,j2)*yj0(i2,i1)
c if(no.eq.1.and.lt.eq.2) write(8,*)'.444',yj0(i1,i2),vectp(i1,j1),
c *vectp(i2,j2),xn,i1,i2
771 continue
do 772 i1=1,ndim
i3=m2+i1
do 772 i2=1,ndim
i4=m2+i2
xh=xh+vectp(i3,j1)*vectp(i4,j2)*yjj3(i2,i1)
c if(no.eq.1.and.lt.eq.2) write(8,*)'.444',yj3(i2,i1),vectp(i3,j1),
c *vectp(i4,j2),xh,i1,i2

772 continue
ccccccccccccccccc
xh1=0.0d0
xh2=0.0d0
do 775 i1=1,ndim
do 775 i2=1,m2
i3=m2+i1
xh1=xh1+vectp(i3,j1)*vectp(i2,j2)*yj1(i1,i2)
xh2=xh2+vectp(i2,j1)*vectp(i3,j2)*yj2(i1,i2)
775 continue
co(j1,j2)=xn+xh+xh1+xh2
if(no.eq.1.and.lt.eq.2) write(8,*)'co',j1,j2,co(j1,j2)
if(no.eq.1.and.lt.eq.2) write(8,*)'..',xn,xh,xh1,xh2
if(lt.ne.2) go to 770
c write(11,*)xn,xh,xk,co(j1,j2)
770 continue
c write(11,*)'.........'
773 continue
ccccccccccccccccccccccccccc
196 if(no.eq.1) go to 194
np=no-2
do 740 j1=1,ndii
do 740 j2=1,ndi
xa=0.0d0
xb=0.0d0
xc=0.0d0
xd=0.0d0
do 741 i=1,mj
do 742 k=1,m2
xa=xa+vectp(k,j2)*vo(i,j1)*xc0(np,k,i)
c if(no.eq.7) write(8,*)'www',vectp(k,j2),vo(i,j1),xc0(np,k,i),xa
742 continue
do 743 k=1,ndim
k1=m2+k
xb=xb+vectp(k1,j2)*vo(i,j1)*xc2(np,k,i)
743 continue
741 continue
do 744 i=1,ndj
i1=mj+i
do 745 k=1,m2
xc=xc+vectp(k,j2)*vo(i1,j1)*xc1(np,k,i)
745 continue
do 746 k=1,ndim
k1=m2+k
xd=xd+vectp(k1,j2)*vo(i1,j1)*xc3(np,k,i)
746 continue
744 continue
co(j1,j2)=xa+xb+xc+xd
c if(no.eq.3) write(8,*)'YYYYY',no,j1,j2,xa,xb,xc,xd,'y',co(j1,j2)
c if(no.eq.7) write(8,*)'.................'
740 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
194 continue
DO 821 IE=1,NINTER BIB04650
if(lt.ne.2) go to 797
WRITE(8,6001)VALP(IE)

WRITE(8,9020) BIB04670
write(iw,*)lt
WRITE(iw,*)valp(Ie)
WRITE(ik,*)valp(Ie)
write(iw,*)ndi
797 continue
id=ie+iimin-1
DO 824 IC=1,NDI BIB04720
INDIC=BLANC BIB04770
IF((IE+IIMIN-1).EQ.IC)INDIC=ETOILE BIB04780
ver1(ic)=0.0d0
if(lt.ne.2) go to 779
ver1(ic)=vectp(ic,ie)
if(ic.le.m2) go to 6666
ver=0.0d0
do 6667 i7=1,ndim
i9=m2+i7
i10=ic-m2
ver=ver+vik(i10,i7)*vectp(i9,ie)
c if(lt.eq.2.and.no.eq.1) write(8,*)vik(i10,i7),vectp(i9,ic),ver
6667 continue
c if(lt.eq.2.and.no.eq.1) write(8,*)'....',ver
ver1(ic)=ver

6666 WRITE(iw,*)indic,nct4(ic),omt4(ic),prtt4(ic),nmtt4(ic),nctt4(ic),
*VECTP(IC,ie),ver1(ic)
c WRITE(ik,*)indic,nmtt4(ic)
WRITE(8,6002)indic,nct4(ic),omt4(ic),prtt4(ic),nmtt4(ic),
*nctt4(ic),prtt4(ic),VECTP(IC,ie),ver1(ic)

c WRITE(8,*)indic,nct4(ic),omt4(ic),prtt4(ic),nmtt4(ic),nctt4(ic),
c *VECTP(IC,ie)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
if(indic.ne.etoile) go to 198
ip=1
write(11,*)ip
if(no.ne.1) go to 197
write(11,*)no,ndi,ic,ndi
do 790 i=1,ndi
write(11,*)nct4(ic),omt4(ic),prtt4(ic),nmtt4(ic),nctt4(ic),
*nct4(i),omt4(i),prtt4(i),nmtt4(i),nctt4(i),co(i,ic)
790 continue
ccccccccccccccccccccccc
197 if(no.eq.1) go to 198
write(11,*)no,ndii,ic,ndi
do 791 i=1,ndii
write(11,*)nct4(ic),omt4(ic),prtt4(ic),nmtt4(ic),nctt4(ic),
*nc1(i),om1(i),pt1(i),nt1(i),nct1(i),co(i,ic)
791 continue
198 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
779 continue
824 continue
822 continue

821 CONTINUE BIB04820

ccccccccccccccccc
40 continue
cccccccccccccccccccccccccccccccccccc
mj=m2
kdj=kk
ndj=ndim
ndii=ndi
do 780 i=1,ndi
nc1(i)=nct4(i)
om1(i)=omt4(i)
pt1(i)=prtt4(i)
nt1(i)=nmtt4(i)
nct1(i)=nctt4(i)
do 780 j=1,ndi
vo(i,j)=vectp(i,j)
780 continue
ccccccccccccccccccccccccccc
20 continue
if(pop.eq.moins) go to 45
pop=moins
go to 2
45 continue
if(lt.ne.2) go to 10
jo=lt+2
write(iw,*)jo
ip=2
write(11,*)ip
write(11,*)fmi
10 continue
c write(iw,*)'5'
stop

ccccccccccc
cccccccccc
5000 format(2x,4hK = ,f7.4,2x,5hMU = ,f7.4,2x,14hDEFORMATION = ,f7.4)
5001 format(2x,4hW = ,f7.4,2x,7hCHI2 = ,f7.4)
5002 format(2x,3hZ =,I3,2x,3hN =,I3,2x,
*31hNb. DE NIVEAU/NIVEAU DE FERMI =,i3)
5003 format(2x,8hDELTA = ,f7.4,2x,5hG0 = ,f7.4)
6001 FORMAT(/,4X,13HVALEUR PROPRE,3X,F7.4)
6002 FORMAT(2H ,A1,1x,i2,1H+,I2,2H/2,A1,1H(,I3,1H),1x,1h=,I2,2h/2,a1,
*4x,F7.4,2f7.4)
6004 FORMAT(1H ,6X,6HPETITA,4X,13(1X,F7.4)) BIB05410
6011 FORMAT(1h ,I2,2H/2,1x,3H>>>,I2,2H/2,A1,1H(,I3,1H),1x,I2,2x,F7.4)
6012 FORMAT(1h ,I2,2H/2,A1,I3,2x,f7.4) BIB05430
6013 FORMAT(F7.4) BIB05440
9020 FORMAT(1H ) BIB05450
9041 FORMAT(10X,6HOMEGA=,I2,2H/2,A1,1H(,I3,1H))
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
SUBROUTINE EIGEN(A,R,N,MV)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(100),R(100)

write(*,*)"--------->",N,MV
5 RANGE=1.0D-06
IF (MV-1) 10,25,10
10 IQ=-N
DO 20 J=1,N
IQ=IQ+N
DO 20 I=1,N
IJ=IQ+I
R(IJ)=0.0D0
IF (I-J) 20,15,20
15 R(IJ)=1.0D0
20 CONTINUE
25 ANORM=0.0D0
DO 35 I=1,N
DO 35 J=I,N
IF (I-J) 30,35,30
30 IA=I+(J*J-J)/2
ANORM=ANORM+A(IA)*A(IA)

35 CONTINUE
IF (ANORM) 165,165,40
40 ANORM=1.414D0*DSQRT(ANORM)
ANRMX=ANORM*RANGE/N
IND=0
THR=ANORM
45 THR=THR/N
50 L=1
55 M=L+1
60 MQ=(M*M-M)/2
LQ=(L*L-L)/2
LM=L+MQ

62 IF (DABS(A(LM))-THR) 130,65,65
65 IND=1
LL=L+LQ
MM=M+MQ
X=0.5D0*(A(LL)-A(MM))
68 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X)
IF (X) 70,75,75
70 Y=-Y
75 SINX=Y/DSQRT(2.0D0*(1.0D0+(DSQRT(1.0D0-Y*Y))))
SINX2=SINX*SINX
78 COSX=DSQRT(1.0D0-SINX2)
COSX2=COSX*COSX
SINCS=SINX*COSX
ILQ=N*(L-1)
IMQ=N*(M-1)
DO 125 I=1,N
IQ=(I*I-I)/2
IF (I-L) 80,115,80
80 IF (I-M) 85,115,90
85 IM=I+MQ
GO TO 95
90 IM=M+IQ
95 IF (I-L) 100,105,105
100 IL=I+LQ
GO TO 110
105 IL=L+IQ
110 X=A(IL)*COSX-A(IM)*SINX
A(IM)=A(IL)*SINX+A(IM)*COSX
A(IL)=X
115 IF (MV-1) 120,125,120
120 ILR=ILQ+I
IMR=IMQ+I
X=R(ILR)*COSX-R(IMR)*SINX
R(IMR)=R(ILR)*SINX+R(IMR)*COSX
R(ILR)=X
125 CONTINUE
X=2.0D0*A(LM)*SINCS
Y=A(LL)*COSX2+A(MM)*SINX2-X
X=A(LL)*SINX2+A(MM)*COSX2+X
A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
A(LL)=Y
A(MM)=X
130 IF (M-N) 135,140,135
135 M=M+1
GO TO 60
140 IF (L-(N-1)) 145,150,145
145 L=L+1
GO TO 55
150 IF (IND-1) 160,155,160
155 IND=0
GO TO 50
160 IF (THR-ANRMX) 165,165,45
165 return
c 165 IQ=-N
DO 185 I=1,N
IQ=IQ+N
LL=I+(I*I-I)/2
JQ=N*(I-2)
DO 185 J=I,N
JQ=JQ+N
MM=J+(J*J-J)/2

IF (A(LL)-A(MM)) 170,185,185
170 X=A(LL)
A(LL)=A(MM)
A(MM)=X
IF (MV-1) 175,185,175

175 DO 180 K=1,N
ILR=IQ+K
IMR=JQ+K
X=R(ILR)
R(ILR)=R(IMR)
180 R(IMR)=X

185 CONTINUE
RETURN
C** THIS PROGRAM VALID ON FTN4 AND FTN5 **
END
subroutine coriol(xj0,xj1,xj2,xj3,yj0,yj1,yj2,yj3,kk,n1,n2,l1,
*l2,l3,k1,k2,k3,k4,k5,k6,x1,x2,x3,ny3,my3,nx3,mx3,nx2,mx2,mx1,
*nx1,xn1,xn2,yn2,y1,y2,ny1,my1,ny2,my2,ny4,my4,nx4,mx4,yn1,y33,
*nn1,nn2,mm1,mm2,fmi,ichoix,no)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER OMG1,omg2,om2,PRT1,prt2,prt,prtt,omt,rtt,omy,prty,omtt1,
*omtt2,ptt1,ptt2
real*8 ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,nn1,mm1,nn2,mm2
dimension x1(10,10),gx1(10,10),fx1(10,10),x2(10,10),gx2(10,10),
*fx2(10,10),x3(10,10),gx3(10,10),fx3(10,10)
dimension ny3(10,10),my3(10,10),nx3(10,10),mx3(10,10),
*nx2(10,10),mx2(10,10)
dimension mx1(10,10),nx1(10,10),xn1(10,10),xn2(10,10),yn2(10,10)
dimension y1(10,10),gy1(10,10),fy1(10,10),y2(10,10),gy2(10,10),
*fy2(10,10),ny1(10,10),my1(10,10),ny2(10,10),my2(10,10)
dimension ny4(10,10),my4(10,10),nx4(10,10),mx4(10,10),yn1(10,10),
*y33(10,10),nn1(10,10),nn2(10,10),mm1(10,10),mm2(10,10)
dimension xj0(10,10),xj1(10,10),xj2(10,10),xj3(10,10),yj0(10,10),
*yj1(10,10),yj2(10,10),yj3(10,10)
character*1 plus,moins,blanc,etoile,pop
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/ BIB04570
plus='+'
moins='-'
blanc=' '
etoile='*'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
do 200 i=1,10
do 200 j=1,10
xj0(i,j)=0.0d0
xj1(i,j)=0.0d0
xj2(i,j)=0.0d0
xj3(i,j)=0.0d0
yj0(i,j)=0.0d0
yj1(i,j)=0.0d0
yj2(i,j)=0.0d0
yj3(i,j)=0.0d0
200 continue
c write(8,*)'............',no,'.....'
c write(8,*)'k2=',k2,'n1=',n1,'k3',k3
c write(8,*)'l2=',l2,'n2=',n2,'k4',k4
c write(8,*)'k6=',k6,'l3=',l3,'k3',k3

do 1 j=1,k2
do 2 i=1,n1
s1=0.0d0
do 3 i2=1,k3
s1=s1+nx2(i,i2)*xn2(i2,j)
3 continue
if(no.eq.1.or.no.eq.3) xj1(i,j)=-s1
if(no.ge.5) xj1(i,j)=s1
c WRITe(8,*)'aaaa xj1',xj1(i,j)
2 continue
do 4 i=1,n2
i1=n1+i
s2=0.0d0
do 5 i2=1,k4
s2=s2+ny2(i,i2)*yn1(j,i2)
5 continue
xj1(i1,j)=s2
c write(8,*)' bbbb xj1',xj1(i1,j)
4 continue
1 continue
cccccccccccccccccccccccccccccccccc
do 10 j=1,l2
do 11 i=1,k3
s3=0.0d0
do 12 i2=1,n1
s3=s3+x2(i2,j)*nx2(i2,i)
12 continue
xj2(j,i)=-s3
c write(8,*)' cccc xj2',xj2(j,i)
11 continue
do 13 i=1,k4
i1=k3+i
s4=0.0d0
do 14 i2=1,n2
s4=s4+y1(i2,j)*ny2(i2,i)
14 continue
xj2(j,i1)=-s4
c write(8,*)' ddddd xj2',xj2(j,i1)
13 continue
10 continue
cccccccccccccccccccccccccccccccccc
do 20 i=1,l2
do 21 j=1,k2
xj0(i,j)=mm2(i,j)
c write(8,*)'wwwww xj0',no,xj0(i,j)

21 continue
if(no.ne.1) go to 20
do 22 j=1,k3
yj0(i,j)=mm1(i,j)
c write(8,*)'ffff yj0',i,j,yj0(i,j)
22 continue
c write(11,*)'rrrrrrrrrr'
20 continue
cccccccccccccccccccccccccccccccccc
do 30 i=1,n1
do 31 j=1,k3
v1=0.0d0
do 32 i3=1,k3
do 32 i2=1,k1
v1=v1+mx2(i,i3)*xn1(i2,i3)*xn1(i2,j)
32 continue
v2=0.0d0
v3=0.0d0
do 33 i4=1,l1
do 34 i2=1,n1
v2=v2+x1(i,i4)*x1(i2,i4)*mx2(i2,j)
34 continue
do 35 i2=1,k1
v3=v3+x1(i,i4)*my3(i4,i2)*xn1(i2,j)
35 continue
33 continue
if(no.eq.1) xj3(i,j)=-mx2(i,j)-v1+v2-v3
if(no.eq.3.or.no.eq.7) xj3(i,j)=mx2(i,j)+v1-v2-v3
if(no.eq.5.or.no.ge.9) xj3(i,j)=mx2(i,j)+v1-v2+v3
c write(8,*)' gggg xj3',xj3(i,j)
31 continue
cccccccccccccccccccccccccccccc
do 40 j=1,k4
j1=k3+j
v4=0.0d0
do 42 i3=1,k3
do 42 i2=1,k2
v4=v4+mx2(i,i3)*xn2(i3,i2)*yn1(i2,j)
42 continue
v5=0.0d0
v6=0.0d0
do 43 i4=1,l2
do 44 i2=1,n2
v5=v5+x2(i,i4)*y1(i2,i4)*my2(i2,j)
44 continue
do 45 i2=1,k2
v6=v6+x2(i,i4)*mm2(i4,i2)*yn1(i2,j)
45 continue
43 continue
if(no.eq.1.or.no.ge.5) xj3(i,j1)=v4-v5+v6
if(no.eq.3) xj3(i,j1)=v4+v5-v6
c write(8,*)' hhhhh xj3',xj3(i,j1)
40 continue
30 continue
cccccccccccccccccccccccccccccc
do 50 i=1,n2
i1=n1+i
do 51 j=1,k3
u1=0.0d0
do 52 i2=1,k2
do 52 i4=1,k4
u1=u1+my2(i,i4)*yn1(i2,i4)*xn2(j,i2)
52 continue
u2=0.0d0
u3=0.0d0
do 53 i4=1,l2
do 54 i2=1,n1
u2=u2+y1(i,i4)*x2(i2,i4)*mx2(i2,j)
54 continue
do 55 i2=1,k2
u3=u3+y1(i,i4)*mm2(i4,i2)*xn2(j,i2)
55 continue
53 continue
if(no.eq.1) xj3(i1,j)=-u1+u2-u3
if(no.eq.3) xj3(i1,j)=u1+u2+u3
if(no.ge.5) xj3(i1,j)=u1-u2+u3
c write(8,*)' mmmmm xj3',xj3(i1,j)
51 continue
cccccccccccccccccccccccccccccc
do 60 j=1,k4
j1=k3+j
u4=0.0d0
do 61 i4=1,k4
do 61 i2=1,k6
u4=u4+my2(i,i4)*y33(i4,i2)*y33(j,i2)
61 continue
u5=0.0d0
u6=0.0d0
do 62 i4=1,l3
do 63 i2=1,n2
u5=u5+y2(i,i4)*y2(i2,i4)*my2(i2,j)
63 continue
do 64 i2=1,k6
u6=u6+y2(i,i4)*my4(i4,i2)*y33(j,i2)
64 continue
62 continue
xj3(i1,j1)=my2(i,j)+u4-u5+u6
c write(8,*)' nnnnn xj3',xj3(i1,j1)
60 continue
50 continue
cccccccccccccccccccccccccccccccccc
if(no.ne.1) go to 80
do 70 i1=1,n1
do 71 i2=1,n1
t1=0.0d0
t2=0.0d0
t11=0.0d0
do 72 j1=1,k1
do 72 j2=1,l2
t1=t1+mx1(i1,j1)*y1(j1,j2)*x2(i2,j2)
t2=t2+x2(i1,j2)*y1(j1,j2)*mx1(i2,j1)
72 continue
do 100 j1=1,l2
do 100 j2=1,l2
t11=t11+x2(i1,j1)*mm1(j1,j2)*x2(i2,j2)
100 continue
yj3(i1,i2)=t1+t2+t11
c write(8,*)'333333333',i1,i2,yj3(i1,i2)
71 continue
do 73 i2=1,n2
i3=n1+i2
t3=0.0d0
t4=0.0d0
t12=0.0d0
do 74 j1=1,k1
do 74 j2=1,l3
t3=t3+mx1(i1,j1)*y2(j1,j2)*y2(i2,j2)
74 continue
do 75 j1=1,l1
do 75 j2=1,n1
t4=t4+x1(i1,j1)*x1(j2,j1)*my1(i2,j2)
75 continue
do 101 j1=1,l1
do 101 j2=1,l3
t12=t12+x1(i1,j1)*mx4(j2,j1)*y2(i2,j2)
101 continue
yj3(i1,i3)=-1.0d0*mx1(i1,i2)+t3+t4-t12
c write(8,*)'22222222222',yj3(i1,i3)
73 continue
70 continue
ccccccccccccccccccc
do 89 i1=1,n2
i4=n1+i1
do 81 i2=1,n1
t5=0.0d0
t6=0.0d0
t13=0.0d0
do 82 j1=1,k2
do 82 j2=1,l1
t5=t5+my1(i1,j1)*x1(j1,j2)*x1(i2,j2)
82 continue
do 83 j1=1,l3
do 83 j2=1,n2
t6=t6+y2(i1,j1)*y2(j2,j1)*mx1(i2,j2)
83 continue
do 102 j1=1,l3
do 102 j2=1,l1
t13=t13+y2(i1,j1)*mx4(j1,j2)*x1(i2,j2)
102 continue
yj3(i4,i2)=-1.0d0*my1(i1,i2)+t5+t6-t13
c write(8,*)'111111111111111',yj3(i4,i2),yj3(i2,i4)
81 continue
do 84 i2=1,n2
i3=n1+i2
t7=0.0d0
t8=0.0d0
t14=0.0d0
do 85 j1=1,k2
do 85 j2=1,l2
t7=t7+my1(i1,j1)*x2(j1,j2)*y1(i2,j2)
t8=t8+y1(i1,j2)*x2(j1,j2)*my1(i2,j1)
85 continue
do 103 j1=1,l2
do 103 j2=1,l2
t14=t14+y1(i1,j1)*mm1(j1,j2)*y1(i2,j2)
103 continue
yj3(i4,i3)=t7+t8+t14
c write(8,*)'444444444',i4,i3,yj3(i1,i2)

84 continue
89 continue
cccccccccccccccccccccccc
do 90 j=1,l2
do 91 i1=1,n1
w1=0.0d0
do 92 i2=1,n2
w1=w1+y1(i2,j)*nx1(i1,i2)
92 continue
yj1(i1,j)=-1.0d0*w1
yj2(i1,j)=-1.0d0*w1
91 continue
do 93 i1=1,n2
i3=n1+i1
w2=0.0d0
do 94 i2=1,n1
w2=w2+x2(i2,j)*ny1(i1,i2)
94 continue
yj1(i3,j)=w2
yj2(i3,j)=w2
93 continue
90 continue
cccccccccccccccccccccccccccccccccccc
80 continue
return
end
subroutine str22(H1,m,a,n1,n2,l1,l2,l3,l5,k1,k2,k3,k4,k5,k6,x1,
*gx1,fx1,x2,gx2,fx2,x3,gx3,fx3,ex,exy,ny3,my3,nx3,mx3,nx2,mx2,mx1,
*nx1,xn1,xn2,yn2,y1,gy1,fy1,y2,gy2,fy2,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,yn1,y33,nn1,nn2,mm1,mm2,uk,uh,ga,lt,ex1,ey1,ey,fmi,
*chi2,rr,w,no,ichoix,sif,sig,sim,sip,siv1,siv2,siw1,siw2,ux,uxy,
*uy)

IMPLICIT REAL*8 (A-H,O-Z)
real*8 ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,nn1,mm1,nn2,mm2
Dimension ga(2)
dimension A(100)
dimension H1(100)
dimension x1(10,10),gx1(10,10),fx1(10,10),x2(10,10),gx2(10,10),
*fx2(10,10),x3(10,10),gx3(10,10),fx3(10,10)
dimension ex(20),ey(20),ex1(20,10),ey1(20,10),exy(20),ny3(10,10),
*my3(10,10),nx3(10,10),mx3(10,10),nx2(10,10),mx2(10,10)
dimension mx1(10,10),nx1(10,10),xn1(10,10),xn2(10,10),yn2(10,10)
dimension y1(10,10),gy1(10,10),fy1(10,10),y2(10,10),gy2(10,10),
*fy2(10,10),ny1(10,10),my1(10,10),ny2(10,10),my2(10,10)
dimension ny4(10,10),my4(10,10),nx4(10,10),mx4(10,10),yn1(10,10),
*y33(10,10),nn1(10,10),nn2(10,10),mm1(10,10),mm2(10,10)
dimension uk(10),uh(10),ux(10,10),uy(10,10),uxy(10,10)

character*1 plus,moins,blanc,etoile,pop
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/
plus='+'
moins='-'
blanc=' '
etoile='*'

m=0

do 20 i1=1,n1
do 21 i2=1,i1
qx0=0.0d0
t0=0.0d0
tl1=0.0d0
tm1=0.0d0
p1=0.0d0
do 22 j1=1,l1
qx0=qx0-x1(i1,j1)*x1(i2,j1)
t0=t0+(ex1(i1,j1)-2.0D0*w)*x1(i1,j1)*x1(i2,j1)
tl1=tl1+fx1(i1,j1)*fx1(i2,j1)
tm1=tm1+gx1(i1,j1)*gx1(i2,j1)
uw=(uk(i1)+uk(i2))*ux(i1,j1)+uk(i1)*uk(i2)
p1=p1+uw*x1(i1,j1)*x1(i2,j1)
22 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
t1=0.0d0
t2=0.0d0
ta1=0.0d0
ta2=0.0d0
ta3=0.0d0
ta4=0.0d0
ta5=0.0d0
ta6=0.0d0
tl2=0.0d0
tl3=0.0d0
tl4=0.0d0
tl5=0.0d0
tl6=0.0d0
tl7=0.0d0
tl8=0.0d0
tm2=0.0d0
tm3=0.0d0
tm4=0.0d0
tm5=0.0d0
tm6=0.0d0
tm7=0.0d0
tm8=0.0d0
cccccccccccccccccc
do 23 ih1=1,l1
do 26 ih2=1,l1
do 24 j1=1,n1
t1=t1+fx1(i1,ih1)*x1(j1,ih1)*x1(j1,ih2)*fx1(i2,ih2)
tl3=tl3+fx1(i1,ih1)*fx1(j1,ih1)*x1(j1,ih2)*x1(i2,ih2)
tl5=tl5+x1(i1,ih1)*x1(j1,ih1)*fx1(j1,ih2)*fx1(i2,ih2)
tl7=tl7+x1(i1,ih1)*fx1(j1,ih1)*fx1(j1,ih2)*x1(i2,ih2)
tm3=tm3+gx1(i1,ih1)*gx1(j1,ih1)*x1(j1,ih2)*x1(i2,ih2)
tm5=tm5+x1(i1,ih1)*x1(j1,ih1)*gx1(j1,ih2)*gx1(i2,ih2)
tm7=tm7+x1(i1,ih1)*gx1(j1,ih1)*gx1(j1,ih2)*x1(i2,ih2)
24 continue
do 25 ih3=1,l5
ta1=ta1+gx1(i1,ih1)*x3(ih1,ih3)*x3(ih2,ih3)*gx1(i2,ih2)
ta3=ta3+x1(i1,ih1)*gx3(ih1,ih3)*x3(ih2,ih3)*gx1(i2,ih2)
ta5=ta5+gx1(i1,ih1)*x3(ih1,ih3)*gx3(ih2,ih3)*x1(i2,ih2)
tl8=tl8+x1(i1,ih1)*fx3(ih1,ih3)*fx3(ih2,ih3)*x1(i2,ih2)
tm8=tm8+x1(i1,ih1)*gx3(ih1,ih3)*gx3(ih2,ih3)*x1(i2,ih2)
25 continue
26 continue
do 27 j1=1,n1
do 27 ih2=1,l2
ta4=ta4+x1(i1,ih1)*gx1(j1,ih1)*x2(j1,ih2)*gx2(i2,ih2)
ta6=ta6+gx2(i1,ih2)*x2(j1,ih2)*gx1(j1,ih1)*x1(i2,ih1)
tl4=tl4+fx2(i1,ih2)*fx2(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
tl6=tl6+x1(i1,ih1)*x1(j1,ih1)*fx2(j1,ih2)*fx2(i2,ih2)
tm4=tm4+gx2(i1,ih2)*gx2(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
tm6=tm6+x1(i1,ih1)*x1(j1,ih1)*gx2(j1,ih2)*gx2(i2,ih2)
27 continue
23 continue
do 30 ih1=1,l2
do 30 ih2=1,l2
do 31 j1=1,n1
ta2=ta2+gx2(i1,ih1)*x2(j1,ih1)*x2(j1,ih2)*gx2(i2,ih2)
31 continue
do 32 j1=1,n2
t2=t2+fx2(i1,ih1)*y1(j1,ih1)*y1(j1,ih2)*fx2(i2,ih2)
32 continue
30 continue
do 50 ih1=1,l2
tl2=tl2+fx2(i1,ih1)*fx2(i2,ih1)
tm2=tm2+gx2(i1,ih1)*gx2(i2,ih1)
50 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
tf1=0.0d0
tf2=0.0d0
tg1=0.0d0
tg2=0.0d0
tg3=0.0d0
tg4=0.0d0
tg5=0.0d0
tg6=0.0d0
do 40 ih1=1,k1
do 40 ih2=1,k1
do 41 ih3=1,k3
tf1=tf1+nx1(i1,ih1)*xn1(ih1,ih3)*xn1(ih2,ih3)*nx1(i2,ih2)
41 continue
do 42 ih3=1,k5
tg1=tg1+mx1(i1,ih1)*yn2(ih1,ih3)*yn2(ih2,ih3)*mx1(i2,ih2)
42 continue
40 continue
do 43 ih1=1,k3
do 43 ih2=1,k3
do 44 ih3=1,k2
tf2=tf2+nx2(i1,ih1)*xn2(ih1,ih3)*xn2(ih2,ih3)*nx2(i2,ih2)
44 continue
do 45 ih3=1,k1
tg2=tg2+mx2(i1,ih1)*xn1(ih3,ih1)*xn1(ih3,ih2)*mx2(i2,ih2)
45 continue
43 continue
do 46 ih1=1,l1
do 46 ih2=1,k1
do 48 ih3=1,k5
tg3=tg3+x1(i1,ih1)*mx3(ih1,ih3)*yn2(ih2,ih3)*mx1(i2,ih2)
tg5=tg5+mx1(i1,ih2)*yn2(ih2,ih3)*mx3(ih1,ih3)*x1(i2,ih1)
48 continue
do 49 ih3=1,k3
tg4=tg4+x1(i1,ih1)*my3(ih1,ih2)*xn1(ih2,ih3)*mx2(i2,ih3)
tg6=tg6+mx2(i1,ih3)*xn1(ih2,ih3)*my3(ih1,ih2)*x1(i2,ih1)
49 continue
46 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
tn1=0.0d0
tn2=0.0d0
tn3=0.0d0
tn4=0.0d0
tn5=0.0d0
tn6=0.0d0
tn7=0.0d0
tn8=0.0d0
to1=0.0d0
to2=0.0d0
to3=0.0d0
to4=0.0d0
to5=0.0d0
to6=0.0d0
to7=0.0d0
to8=0.0d0

do 51 ih1=1,k1
tn1=tn1+nx1(i1,ih1)*nx1(i2,ih1)
to1=to1+mx1(i1,ih1)*mx1(i2,ih1)
51 continue
do 52 ih1=1,k3
tn2=tn2+nx2(i1,ih1)*nx2(i2,ih1)
to2=to2+mx2(i1,ih1)*mx2(i2,ih1)
52 continue
do 53 ih1=1,l1
do 54 j1=1,n1
do 55 ih2=1,k1
tn3=tn3+nx1(i1,ih2)*nx1(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
tn5=tn5+x1(i1,ih1)*x1(j1,ih1)*nx1(j1,ih2)*nx1(i2,ih2)
to3=to3+mx1(i1,ih2)*mx1(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
to5=to5+x1(i1,ih1)*x1(j1,ih1)*mx1(j1,ih2)*mx1(i2,ih2)
55 continue
do 56 ih2=1,k3
tn4=tn4+nx2(i1,ih2)*nx2(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
tn6=tn6+x1(i1,ih1)*x1(j1,ih1)*nx2(j1,ih2)*nx2(i2,ih2)
to4=to4+mx2(i1,ih2)*mx2(j1,ih2)*x1(j1,ih1)*x1(i2,ih1)
to6=to6+x1(i1,ih1)*x1(j1,ih1)*mx2(j1,ih2)*mx2(i2,ih2)
56 continue
54 continue
do 57 ih2=1,l1
do 58 ih3=1,k1
tn7=tn7+x1(i1,ih1)*ny3(ih1,ih3)*ny3(ih2,ih3)*x1(i2,ih2)
to7=to7+x1(i1,ih1)*my3(ih1,ih3)*my3(ih2,ih3)*x1(i2,ih2)
58 continue
do 59 ih3=1,k5
tn8=tn8+x1(i1,ih1)*nx3(ih1,ih3)*nx3(ih2,ih3)*x1(i2,ih2)
to8=to8+x1(i1,ih1)*mx3(ih1,ih3)*mx3(ih2,ih3)*x1(i2,ih2)
59 continue
57 continue
53 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
u022=t0-chi2*(t1+t2)
if(no.le.3) u22q=chi2*(ta1+ta2+ta3-ta4+ta5-ta6)
if(no.eq.5.or.no.eq.7) u22q=chi2*(ta1+ta2-ta3-ta4-ta5-ta6)
if(no.eq.9.or.no.eq.11) u22q=chi2*(ta1+ta2-ta3+ta4-ta5+ta6)
if(no.ge.13) u22q=chi2*(ta1+ta2+ta3+ta4+ta5+ta6)
u22j=fmi*(tf1+tf2)
if(no.eq.1.or.no.ge.11) u2j2=fmi*(-tg1-tg2+tg3+tg4+tg5+tg6)
if(no.eq.3.or.no.eq.7) u2j2=fmi*(-tg1-tg2+tg3-tg4+tg5-tg6)
if(no.eq.5.or.no.eq.9) u2j2=fmi*(-tg1-tg2-tg3+tg4-tg5+tg6)
u22p=2.0d0*ga(lt)*p1
u11q=(chi2/2.0d0)*(tl1+tl2-tl3-tl4-tl5-tl6-tl7-tl8)
u1q1=(chi2/2.0d0)*(-tm1-tm2+tm3+tm4+tm5+tm6+tm7+tm8)
u11j=(fmi/2.0d0)*(-tn1-tn2+tn3+tn4+tn5+tn6+tn7+tn8)
u1j1=(fmi/2.0d0)*(to1+to2-to3-to4-to5-to6-to7-to8)
cccccccccccccccccccccccccccccccccccccccccc
m=m+1
a(m)=qx0
go to(301,302,303,304,305,306,307,308,309),ichoix
301 H1(m)=u022+u22q
go to 320
302 H1(m)=u022+u22q+u11q+u1q1
go to 320
303 H1(m)=u022+u22q+u22j+u2j2
go to 320
304 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1
go to 320
305 H1(m)=u022+u22q+u22j+u2j2+u11j+u1j1
go to 320
306 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 320
307 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j
go to 320
308 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 320
309 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1+u22p
c 309 H1(m)=u022+u11j+u1j1

ccccccccccccccccccccccccccccccc
320 continue
siff=0.0d0
if(i1.ne.i2) go to 21
a(m)=1+a(m)
siff=-(2.0d0*ga(lt)*uk(i1))*sif
c write(*,*)'siff',siff,'sig',sig,'....'
go to(311,312,313,314,315,316,317,318,319),ichoix
311 H1(m)=(ex(i1)+w)+H1(m)+sim
go to 330
312 H1(m)=(ex(i1)+w)+H1(m)+sim+siv1+siv2
go to 330
313 H1(m)=(ex(i1)+w)+H1(m)+sim+sip
go to 330
314 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siv1+siv2
go to 330
315 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siw1+siw2
go to 330
316 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2
go to 330
317 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siv1+siv2+siw1
go to 330
318 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2
go to 330
319 H1(m)=(ex(i1)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2+sig+siff
c 319 H1(m)=(ex(i1)+w)+H1(m)+siw1+siw2

go to 330
ccccccccccccccccccccccccccccccccccccccccccccccccc
330 continue
21 continue
20 continue
ccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccc
do 65 i3=1,n2
do 66 i4=1,n1
qxy=0.0d0
t6=0.0d0
p3=0.0d0
do 67 j3=1,l2
qxy=qxy-x2(i4,j3)*y1(i3,j3)
t6=t6+(exy(j3)-2.0D0*w)*x2(i4,j3)*y1(i3,j3)
uw=(uk(i4)+uh(i3))*uxy(i4,j3)+uk(i4)*uh(i3)
p3=p3+uw*x2(i4,j3)*y1(i3,j3)
67 continue
cccccccccccccccccccccccccccccccccccccccccccccccc
t7=0.0d0
t8=0.0d0
tc1=0.0d0
tc2=0.0d0
tc3=0.0d0
tc4=0.0d0
tc5=0.0d0
tc6=0.0d0
tu1=0.0d0
tu2=0.0d0
tu3=0.0d0
tu4=0.0d0
tu5=0.0d0
tu6=0.0d0
tu7=0.0d0
tu8=0.0d0
ty1=0.0d0
ty2=0.0d0
ty3=0.0d0
ty4=0.0d0
ty5=0.0d0
ty6=0.0d0
ty7=0.0d0
ty8=0.0d0
ccccccccccccccccccccccccccccccccccc
do 70 ih2=1,l2
do 71 j1=1,n1
do 72 ih1=1,l1
t7=t7+fx1(i4,ih1)*x1(j1,ih1)*x2(j1,ih2)*fy1(i3,ih2)
tc1=tc1+gx1(i4,ih1)*x1(j1,ih1)*x2(j1,ih2)*gy1(i3,ih2)
tc5=tc5+gx1(i4,ih1)*x1(j1,ih1)*gx2(j1,ih2)*y1(i3,ih2)
tu1=tu1+fx1(i4,ih1)*fx1(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
ty1=ty1+gx1(i4,ih1)*gx1(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
72 continue
do 73 ih1=1,l2
tc3=tc3+x2(i4,ih1)*gx2(j1,ih1)*x2(j1,ih2)*gy1(i3,ih2)
tu2=tu2+fx2(i4,ih1)*fx2(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
tu6=tu6+x2(i4,ih1)*fx2(j1,ih1)*fx2(j1,ih2)*y1(i3,ih2)
ty2=ty2+gx2(i4,ih1)*gx2(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
ty6=ty6+x2(i4,ih1)*gx2(j1,ih1)*gx2(j1,ih2)*y1(i3,ih2)
73 continue
71 continue
do 74 j3=1,n2
do 75 ih3=1,l3
t8=t8+fx2(i4,ih2)*y1(j3,ih2)*y2(j3,ih3)*fy2(i3,ih3)
tc2=tc2+gx2(i4,ih2)*y1(j3,ih2)*y2(j3,ih3)*gy2(i3,ih3)
tc4=tc4+x2(i4,ih2)*gy1(j3,ih2)*y2(j3,ih3)*gy2(i3,ih3)
tu4=tu4+x2(i4,ih2)*y1(j3,ih2)*fy2(j3,ih3)*fy2(i3,ih3)
ty4=ty4+x2(i4,ih2)*y1(j3,ih2)*gy2(j3,ih3)*gy2(i3,ih3)
75 continue
do 76 ih1=1,l2
tc6=tc6+gx2(i4,ih1)*y1(j3,ih1)*gy1(j3,ih2)*y1(i3,ih2)
tu3=tu3+x2(i4,ih1)*y1(j3,ih1)*fy1(j3,ih2)*fy1(i3,ih2)
tu5=tu5+x2(i4,ih1)*fy1(j3,ih1)*fy1(j3,ih2)*y1(i3,ih2)
ty3=ty3+x2(i4,ih1)*y1(j3,ih1)*gy1(j3,ih2)*gy1(i3,ih2)
ty5=ty5+x2(i4,ih1)*gy1(j3,ih1)*gy1(j3,ih2)*y1(i3,ih2)
76 continue
74 continue
70 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
tf5=0.0d0
tf6=0.0d0
tk1=0.0d0
tk2=0.0d0
tk3=0.0d0
tk4=0.0d0
tk5=0.0d0
tk6=0.0d0
tv1=0.0d0
tv2=0.0d0
tv3=0.0d0
tv4=0.0d0
tv5=0.0d0
tv6=0.0d0
tw1=0.0d0
tw2=0.0d0
tw3=0.0d0
tw4=0.0d0
tw5=0.0d0
tw6=0.0d0
ccccccccccccccccccccccccccccccccccc
do 80 ih3=1,k3
do 81 ih2=1,k2
do 82 ih1=1,k1
tf5=tf5+nx1(i4,ih1)*xn1(ih1,ih3)*xn2(ih3,ih2)*ny1(i3,ih2)
tk1=tk1+mx1(i4,ih1)*xn1(ih1,ih3)*xn2(ih3,ih2)*my1(i3,ih2)
82 continue
do 83 ih1=1,l2
tk3=tk3+x2(i4,ih1)*mm1(ih1,ih3)*xn2(ih3,ih2)*my1(i3,ih2)
83 continue
do 84 ih4=1,k4
tf6=tf6+nx2(i4,ih3)*xn2(ih3,ih2)*yn1(ih2,ih4)*ny2(i3,ih4)
tk2=tk2+mx2(i4,ih3)*xn2(ih3,ih2)*yn1(ih2,ih4)*my2(i3,ih4)
84 continue
81 continue
do 85 ih2=1,l2
do 86 ih1=1,k1
tk5=tk5+mx1(i4,ih1)*xn1(ih1,ih3)*mm1(ih2,ih3)*y1(i3,ih2)
86 continue
do 87 ih1=1,k2
tk6=tk6+mx2(i4,ih3)*xn2(ih3,ih1)*mm2(ih2,ih1)*y1(i3,ih2)
87 continue
85 continue
80 continue
do 88 ih1=1,l2
do 88 ih2=1,k2
do 88 ih4=1,k4
tk4=tk4+x2(i4,ih1)*mm2(ih1,ih2)*yn1(ih2,ih4)*my2(i3,ih4)
88 continue
cccccccccccccccccccccccccccccccc
do 90 ih2=1,l2
do 91 j1=1,n1
do 92 ih1=1,k1
tv1=tv1+nx1(i4,ih1)*nx1(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
tw1=tw1+mx1(i4,ih1)*mx1(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
92 continue
do 93 ih1=1,k3
tv2=tv2+nx2(i4,ih1)*nx2(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
tw2=tw2+mx2(i4,ih1)*mx2(j1,ih1)*x2(j1,ih2)*y1(i3,ih2)
93 continue
91 continue
do 94 j1=1,n2
do 95 ih1=1,k2
tv3=tv3+x2(i4,ih2)*y1(j1,ih2)*ny1(j1,ih1)*ny1(i3,ih1)
tw3=tw3+x2(i4,ih2)*y1(j1,ih2)*my1(j1,ih1)*my1(i3,ih1)
95 continue
do 96 ih4=1,k4
tv4=tv4+x2(i4,ih2)*y1(j1,ih2)*ny2(j1,ih4)*ny2(i3,ih4)
tw4=tw4+x2(i4,ih2)*y1(j1,ih2)*my2(j1,ih4)*my2(i3,ih4)
96 continue
94 continue
do 97 ih1=1,l2
do 98 j3=1,k2
tv5=tv5+x2(i4,ih1)*nn2(ih1,j3)*nn2(ih2,j3)*y1(i3,ih2)
tw5=tw5+x2(i4,ih1)*mm2(ih1,j3)*mm2(ih2,j3)*y1(i3,ih2)
98 continue
do 99 j3=1,k3
tv6=tv6+x2(i4,ih1)*nn1(ih1,j3)*nn1(ih2,j3)*y1(i3,ih2)
tw6=tw6+x2(i4,ih1)*mm1(ih1,j3)*mm1(ih2,j3)*y1(i3,ih2)
99 continue
97 continue
90 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
if(no.eq.1) u022=t6+chi2*(t7+t8)
if(no.eq.3) u022=-t6-chi2*(t7+t8)
if(no.ge.5) u022=t6-chi2*(t7+t8)

if(no.eq.1) u22q=chi2*(tc1-tc2-tc3+tc4-tc5-tc6)
if(no.eq.3) u22q=chi2*(-tc1+tc2+tc3-tc4+ta5+ta6)
if(no.eq.5) u22q=chi2*(-tc1+tc2+tc3+tc4-tc5+tc6)
if(no.eq.7) u22q=chi2*(-tc1+tc2+tc3+tc4-tc5+tc6)
if(no.ge.9) u22q=chi2*(tc1+tc2+tc3+tc4+tc5+tc6)

if(no.eq.3) u22j=fmi*(-tf6-tf5)
if(no.ne.3) u22j=fmi*(tf6+tf5)
if(no.eq.1) u2j2=fmi*(-tk1-tk2-tk3+tk4-tk5+tk6)
if(no.eq.3) u2j2=fmi*(tk1-tk2-tk3-tk4-tk5+tk6)
if(no.eq.5) u2j2=fmi*(tk1-tk2+tk3+tk4-tk5+tk6)
if(no.ge.7) u2j2=fmi*(-tk1-tk2+tk3+tk4+tk5+tk6)

u22p=2.0d0*ga(lt)*p3
u11q=-(chi2/2.0d0)*(tu1+tu2+tu3+tu4+tu5+tu6)
u1q1=(chi2/2.0d0)*(ty1+ty2+ty3+ty4+ty5+ty6)
u11j=(fmi/2.0d0)*(tv1+tv2+tv3+tv4+tv5+tv6)
u1j1=-(fmi/2.0d0)*(tw1+tw2+tw3+tw4+tw5+tw6)
if(no.ne.3) go to 399
u22p=-u22p
u11q=-u11q
u1q1=-u1q1
u11j=-u11j
u1j1=-u1j1
399 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
c write(8,*)'h1-2',u022,u22q,u22p,u11q,u1q1,u22j,u2j2,u11j,u1j1

m=m+1
a(m)=qxy
if(no.eq.3) a(m)=-a(m)
go to(321,322,323,324,325,326,327,328,329),ichoix
321 H1(m)=u022+u22q
go to 340
322 H1(m)=u022+u22q+u11q+u1q1
go to 340
323 H1(m)=u022+u22q+u22j+u2j2
go to 340
324 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1
go to 340
325 H1(m)=u022+u22q+u22j+u2j2+u11j+u1j1
go to 340
326 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 340
327 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j
go to 340
328 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 340
329 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1+u22p
c 329 H1(m)=u022+u11j+u1j1
340 continue
66 continue

cccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccc
do 68 i5=1,i3
qy0=0.0d0
t3=0.0d0
tp2=0.0d0
tq2=0.0d0
p2=0.0d0
do 69 j3=1,l3
qy0=qy0-y2(i3,j3)*y2(i5,j3)
t3=t3+(ey1(i3,j3)-2.0D0*w)*y2(i3,j3)*y2(i5,j3)
tp2=tp2+fy2(i3,j3)*fy2(i5,j3)
tq2=tq2+gy2(i3,j3)*gy2(i5,j3)
uw=(uh(i3)+uh(i5))*uy(i3,j3)+uh(i3)*uh(i5)
p2=p2+uw*y2(i3,j3)*y2(i5,j3)
c write(8,*)'bbbbbb',i3,i5,j3,t3,ey1(i3,j3),w
69 continue
cccccccccccccccccccccccccccccccccccccccccccccccc
t4=0.0d0
t5=0.0d0
tb1=0.0d0
tb2=0.0d0
tb3=0.0d0
tb4=0.0d0
tb5=0.0d0
tb6=0.0d0
tp1=0.0d0
tp3=0.0d0
tp4=0.0d0
tp5=0.0d0
tp6=0.0d0
tp7=0.0d0
tp8=0.0d0
tq1=0.0d0
tq3=0.0d0
tq4=0.0d0
tq5=0.0d0
tq6=0.0d0
tq7=0.0d0
tq8=0.0d0
cccccccccccccccccc
c tb2=tb4=tb6=0.0d0
c tp7=tq7=0.0d0
do 110 ih1=1,l2
do 110 ih2=1,l2
do 111 j1=1,n1
t4=t4+fy1(i3,ih1)*x2(j1,ih1)*x2(j1,ih2)*fy1(i5,ih2)
111 continue
do 112 j1=1,n2
tb1=tb1+gy1(i3,ih1)*y1(j1,ih1)*y1(j1,ih2)*gy1(i5,ih2)
112 continue
110 continue
do 113 ih1=1,l3
do 114 ih3=1,l3
do 115 j3=1,n2
t5=t5+fy2(i3,ih1)*y2(j3,ih1)*y2(j3,ih3)*fy2(i5,ih3)
tp4=tp4+fy2(i3,ih1)*fy2(j3,ih1)*y2(j3,ih3)*y2(i5,ih3)
tp6=tp6+y2(i3,ih1)*y2(j3,ih1)*fy2(j3,ih3)*fy2(i5,ih3)
tp8=tp8+y2(i3,ih1)*fy2(j3,ih1)*fy2(j3,ih3)*y2(i5,ih3)
tq4=tq4+gy2(i3,ih1)*gy2(j3,ih1)*y2(j3,ih3)*y2(i5,ih3)
tq6=tq6+y2(i3,ih1)*y2(j3,ih1)*gy2(j3,ih3)*gy2(i5,ih3)
tq8=tq8+y2(i3,ih1)*gy2(j3,ih1)*gy2(j3,ih3)*y2(i5,ih3)
115 continue
114 continue
do 116 ih2=1,l2
do 116 j3=1,n2
tb3=tb3+y2(i3,ih1)*gy2(j3,ih1)*y1(j3,ih2)*gy1(i5,ih2)
tb5=tb5+gy1(i3,ih2)*y1(j3,ih2)*gy2(j3,ih1)*y2(i5,ih1)
tp3=tp3+fy1(i3,ih2)*fy1(j3,ih2)*y2(j3,ih1)*y2(i5,ih1)
tp5=tp5+y2(i3,ih1)*y2(j3,ih1)*fy1(j3,ih2)*fy1(i5,ih2)
tq3=tq3+gy1(i3,ih2)*gy1(j3,ih2)*y2(j3,ih1)*y2(i5,ih1)
tq5=tq5+y2(i3,ih1)*y2(j3,ih1)*gy1(j3,ih2)*gy1(i5,ih2)
116 continue
113 continue
do 130 ih3=1,l2
tp1=tp1+fy1(i3,ih3)*fy1(i5,ih3)
tq1=tq1+gy1(i3,ih3)*gy1(i5,ih3)
130 continue
ccccccccccccccccccccccccccccccccccccccccccccccccc
tf3=0.0d0
tf4=0.0d0
th1=0.0d0
th2=0.0d0
th3=0.0d0
th4=0.0d0
th5=0.0d0
th6=0.0d0
tr1=0.0d0
tr2=0.0d0
tr3=0.0d0
tr4=0.0d0
tr5=0.0d0
tr6=0.0d0
tr7=0.0d0
tr8=0.0d0
ts1=0.0d0
ts2=0.0d0
ts3=0.0d0
ts4=0.0d0
ts5=0.0d0
ts6=0.0d0
ts7=0.0d0
ts8=0.0d0

do 120 ih1=1,k2
do 120 ih2=1,k2
do 121 ih3=1,k3
tf3=tf3+ny1(i3,ih1)*xn2(ih3,ih1)*xn2(ih3,ih2)*ny1(i5,ih2)
121 continue
do 122 ih4=1,k4
th1=th1+my1(i3,ih1)*yn1(ih1,ih4)*yn1(ih2,ih4)*my1(i5,ih2)
122 continue
120 continue
do 123 ih1=1,k4
do 123 ih4=1,k4
do 124 ih2=1,k2
tf4=tf4+ny2(i3,ih1)*yn1(ih2,ih1)*yn1(ih2,ih4)*ny2(i5,ih4)
124 continue
do 125 ih6=1,k6
th2=th2+my2(i3,ih1)*y33(ih1,ih6)*y33(ih4,ih6)*my2(i5,ih4)
125 continue
123 continue
do 126 ih3=1,l3
do 126 ih4=1,k4
do 126 ih2=1,k2
th3=th3+y2(i3,ih3)*mx4(ih3,ih4)*yn1(ih2,ih4)*my1(i5,ih2)
th5=th5+my1(i3,ih2)*yn1(ih2,ih4)*mx4(ih3,ih4)*y2(i5,ih3)
126 continue
do 127 ih3=1,l3
do 127 ih6=1,k6
do 127 ih4=1,k4
th4=th4+y2(i3,ih3)*my4(ih3,ih6)*y33(ih4,ih6)*my2(i5,ih4)
th6=th6+my2(i3,ih4)*y33(ih4,ih6)*my4(ih3,ih6)*y2(i5,ih3)
127 continue
ccccccccccccccccccccccccccccccccccccccccccccc
do 140 ih2=1,k2
tr1=tr1+ny1(i3,ih2)*ny1(i5,ih2)
ts1=ts1+my1(i3,ih2)*my1(i5,ih2)
140 continue
do 141 ih4=1,k4
tr2=tr2+ny2(i3,ih4)*ny2(i5,ih4)
ts2=ts2+my2(i3,ih4)*my2(i5,ih4)
141 continue
do 142 ih3=1,l3
do 143 j3=1,n2
do 144 ih2=1,k2
tr3=tr3+ny1(i3,ih2)*ny1(j3,ih2)*y2(j3,ih3)*y2(i5,ih3)
tr5=tr5+y2(i3,ih3)*y2(j3,ih3)*ny1(j3,ih2)*ny1(i5,ih2)
ts3=ts3+my1(i3,ih2)*my1(j3,ih2)*y2(j3,ih3)*y2(i5,ih3)
ts5=ts5+y2(i3,ih3)*y2(j3,ih3)*my1(j3,ih2)*my1(i5,ih2)
144 continue
do 145 ih4=1,k4
tr4=tr4+ny2(i3,ih4)*ny2(j3,ih4)*y2(j3,ih3)*y2(i5,ih3)
tr6=tr6+y2(i3,ih3)*y2(j3,ih3)*ny2(j3,ih4)*ny2(i5,ih4)
ts4=ts4+my2(i3,ih4)*my2(j3,ih4)*y2(j3,ih3)*y2(i5,ih3)
ts6=ts6+y2(i3,ih3)*y2(j3,ih3)*my2(j3,ih4)*my2(i5,ih4)
145 continue
143 continue
do 146 ih1=1,l3
do 147 ih6=1,k6
tr7=tr7+y2(i3,ih3)*ny4(ih3,ih6)*ny4(ih1,ih6)*y2(i5,ih1)
ts7=ts7+y2(i3,ih3)*my4(ih3,ih6)*my4(ih1,ih6)*y2(i5,ih1)
147 continue
do 148 ih4=1,k4
tr8=tr7+y2(i3,ih3)*nx4(ih3,ih4)*nx4(ih1,ih4)*y2(i5,ih1)
ts8=ts8+y2(i3,ih3)*mx4(ih3,ih4)*mx4(ih1,ih4)*y2(i5,ih1)
148 continue
146 continue
142 continue
ccccccccccccccccccccccccccccccccccccccccccc
u022=t3-chi2*(t4+t5)
u22q=chi2*(tb1+tb2+tb3+tb4+tb5+tb6)
u22j=fmi*(tf3+tf4)
u2j2=fmi*(th3+th4+th5+th6+th7+th8-th1-th2)
u22p=2.0d0*ga(lt)*p2
u11q=(chi2/2.0d0)*(tp1+tp2-tp3-tp4-tp5-tp6-tp7-tp8)
u1q1=(chi2/2.0d0)*(-tq1-tq2+tq3+tq4+tq5+tq6+tq7+tq8)
u11j=(fmi/2.0d0)*(-tr1-tr2+tr3+tr4+tr5+tr6+tr7+tr8)
u1j1=(fmi/2.0d0)*(ts1+ts2-ts3-ts4-ts5-ts6-ts7-ts8)
c write(8,*)'aaaaaa',i3,i5,u022
ccccccccccccccccccccccccccccccccccccccccc
m=m+1
a(m)=qy0
go to(351,352,353,354,355,356,357,358,359),ichoix
351 H1(m)=u022+u22q
go to 360
352 H1(m)=u022+u22q+u11q+u1q1
go to 360
353 H1(m)=u022+u22q+u22j+u2j2
go to 360
354 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1
go to 360
355 H1(m)=u022+u22q+u22j+u2j2+u11j+u1j1
go to 360
356 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 360
357 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j
go to 360
358 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1
go to 360
359 H1(m)=u022+u22q+u22j+u2j2+u11q+u1q1+u11j+u1j1+u22p
c 359 H1(m)=u022+u11j+u1j1

ccccccccccccccccccccccccccccccc
360 continue

siff=0.0d0
if(i5.ne.i3) go to 68
a(m)=1+a(m)
siff=-(2.0d0*ga(lt)*uh(i3))*sif
c write(*,*)'siff',siff,'sig',sig,'....'
go to(361,362,363,364,365,366,367,368,369),ichoix
361 H1(m)=(ey(i3)+w)+H1(m)+sim
go to 370
362 H1(m)=(ey(i3)+w)+H1(m)+sim+siv1+siv2
go to 370
363 H1(m)=(ey(i3)+w)+H1(m)+sim+sip
go to 370
364 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siv1+siv2
go to 370
365 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siw1+siw2
go to 370
366 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2
go to 370
367 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siv1+siv2+siw1
go to 370
368 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2
go to 370
369 H1(m)=(ey(i3)+w)+H1(m)+sim+sip+siv1+siv2+siw1+siw2+sig+siff
c 369 H1(m)=(ey(i3)+w)+H1(m)+siw1+siw2

go to 370
370 continue
68 continue
65 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
RETURN
END
subroutine str31(H3,nm,n1,n2,l1,l2,l3,k1,k2,k3,k4,k5,k6,x1,gx1,
*fx1,x2,gx2,fx2,x3,gx3,fx3,ex,exy,ny3,my3,nx3,mx3,nx2,mx2,mx1,
*nx1,xn1,xn2,yn2,y1,gy1,fy1,y2,gy2,fy2,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,yn1,y33,nn1,nn2,mm1,mm2,uk,uh,ga,lt,fmi,chi2,rr,w,h0,
*ml,ex1,ey1,ey,ichoix,no,fdm,sif,sig,sim,sip,siv1,siv2,siw1,
*siw2)

IMPLICIT REAL*8 (A-H,O-Z)
INTEGER OMG1,omg2,om2,PRT1,prt2,prt,prtt,omt,rtt,omy,prty,omtt1,
*omtt2,ptt1,ptt2
real*8 ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,nn1,mm1,nn2,mm2

Dimension ga(2)
c dimension A(1)
dimension H3(20,20),H0(20)
dimension x1(10,10),gx1(10,10),fx1(10,10),x2(10,10),gx2(10,10),
*fx2(10,10),x3(10,10),gx3(10,10),fx3(10,10)
dimension ex(20),ey(20),ex1(20,10),ey1(20,10),exy(20),ny3(10,10),
*my3(10,10),nx3(10,10),mx3(10,10),nx2(10,10),mx2(10,10)
dimension mx1(10,10),nx1(10,10),xn1(10,10),xn2(10,10),yn2(10,10)
dimension y1(10,10),gy1(10,10),fy1(10,10),y2(10,10),gy2(10,10),
*fy2(10,10),ny1(10,10),my1(10,10),ny2(10,10),my2(10,10)
dimension ny4(10,10),my4(10,10),nx4(10,10),mx4(10,10),yn1(10,10),
*y33(10,10),nn1(10,10),nn2(10,10),mm1(10,10),mm2(10,10)
dimension uk(10),uh(10),ux(10,10),uy(10,10),uxy(10,10)

character*1 plus,moins,blanc,etoile,pop
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/ BIB04570
plus='+'
moins='-'
blanc=' '
etoile='*'
cccccccccccccccccccccccccccccccccccc
do 1000 i=1,k3
do 1000 j=1,k2
c write(8,*)'....',xn2(i,j)
1000 continue

do 790 j1=1,l2
nm=0
td0=0.0d0

do 791 i1=1,n1
td0=rr*gx2(i1,j1)
if(no.eq.7)write(8,*)'8888',rr,gx2(i1,j1),td0
cccccccccccccccccccccccc
td1=0.0d0
td2=0.0d0
te1=0.0d0
te2=0.0d0
te3=0.0d0
te4=0.0d0
ti1=0.0d0
ti2=0.0d0
tj1=0.0d0
tj2=0.0d0
tj3=0.0d0
tj4=0.0d0

do 710 j2=1,n1
do 711 ih1=1,l1
td1=td1+gx2(j2,j1)*x1(j2,ih1)*fx1(i1,ih1)
te1=te1+x2(j2,j1)*gx1(j2,ih1)*fx1(i1,ih1)
te2=te2+x2(j2,j1)*fx1(j2,ih1)*gx1(i1,ih1)
711 continue
do 712 ih2=1,l2
te3=te3+x2(j2,j1)*gx2(j2,ih2)*fx2(i1,ih2)
te4=te4+x2(j2,j1)*fx2(j2,ih2)*gx2(i1,ih2)
712 continue
710 continue
do 720 j2=1,n2
do 720 ih2=1,l2
td2=td2+gy1(j2,j1)*y1(j2,ih2)*fx2(i1,ih2)
720 continue
do 730 ih3=1,k3
do 731 ih1=1,k1
ti1=ti1+mm1(j1,ih3)*xn1(ih1,ih3)*nx1(i1,ih1)
c write(8,*)'mmmmmmmm',mm1(j1,ih3),xn1(ih1,ih3),
c *nx1(i1,ih1)
731 continue
do 732 ih2=1,k2
ti2=ti2+mm2(j1,ih2)*xn2(ih3,ih2)*nx2(i1,ih3)
c write(8,*)'mmm...mmmmm',mm2(j1,ih2),xn2(ih3,ih2),
c *nx2(i1,ih3)
732 continue
do 733 i2=1,n1
tj3=tj3+x2(i2,j1)*nx2(i2,ih3)*mx2(i1,ih3)
tj4=tj4+x2(i2,j1)*mx2(i2,ih3)*nx2(i1,ih3)
733 continue
730 continue
do 734 ih1=1,k1
do 734 i2=1,n1
tj1=tj1+x2(i2,j1)*nx1(i2,ih1)*mx1(i1,ih1)
tj2=tj2+x2(i2,j1)*mx1(i2,ih1)*nx1(i1,ih1)
734 continue
ccccccccccccccccccccccccccccccccccccccc
if(no.eq.1) u31q=td0-chi2*(td1+td2)
if(no.eq.3) u31q=-td0+chi2*(td1+td2)
if(no.eq.5.or.no.eq.7) u31q=-td0+chi2*(-td1+td2)
if(no.ge.9) u31q=-td0+chi2*(td1+td2)
if(no.eq.1) u20q=(chi2/2.0d0)*(te1+te2+te3+te4)
if(no.eq.3) u20q=(chi2/2.0d0)*(-te1-te2-te3-te4)
if(no.ge.5) u20q=(chi2/2.0d0)*(te1+te2+te3+te4)
if(no.eq.1) u31j=-fmi*(ti1+ti2)
if(no.eq.3) u31j=fmi*(-ti1-ti2)
if(no.eq.5) u31j=fmi*(ti1+ti2)
if(no.ge.7) u31j=fmi*(-ti1+ti2)
if(no.eq.1) u20j=(fmi/2.0d0)*(-tj1-tj2+tj3+tj4)
if(no.eq.3) u20j=(fmi/2.0d0)*(tj1+tj2-tj3-tj4)
if(no.ge.5) u20j=(fmi/2.0d0)*(tj1+tj2-tj3-tj4)
if(no.eq.7)write(8,*)'h31...',u31q,u20q,u31j,u20j
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
nm=nm+1
go to(351,352,353,354,355,356,357,358,359),ichoix
351 H3(j1,nm)=u31q
go to 360
352 H3(j1,nm)=u31q+u20q
go to 360
353 H3(j1,nm)=u31q+u31j
go to 360
354 H3(j1,nm)=u31q+u31j+u20q
go to 360
355 H3(j1,nm)=u31q+u31j+u20j
go to 360
356 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 360
357 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 360
358 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 360
359 H3(j1,nm)=u31q+u31j+u20q+u20j
c 359 H3(j1,nm)=u31q+u20j

360 continue
791 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
td3=0.0d0
do 792 i1=1,n2
td3=-rr*gy1(i1,j1)
cccccccccccccccccccccccccccc
td4=0.0d0
td5=0.0d0
te5=0.0d0
te6=0.0d0
te7=0.0d0
te8=0.0d0
ti3=0.0d0
ti4=0.0d0
tj5=0.0d0
tj6=0.0d0
tj7=0.0d0
tj8=0.0d0

do 740 i2=1,n2
do 741 ih3=1,l3
td5=td5+gy1(i2,j1)*y2(i2,ih3)*fy2(i1,ih3)
te7=te7+y1(i2,j1)*gy2(i2,ih3)*fy2(i1,ih3)
te8=te8+y1(i2,j1)*fy2(i2,ih3)*gy2(i1,ih3)
741 continue
do 742 ih2=1,l2
te5=te5+y1(i2,j1)*gy1(i2,ih2)*fy1(i1,ih2)
te6=te6+y1(i2,j1)*fy1(i2,ih2)*gy1(i1,ih2)
742 continue
740 continue
do 743 i2=1,n1
do 743 ih2=1,l2
td4=td4+gx2(i2,j1)*x2(i2,ih2)*fy1(i1,ih2)
743 continue
do 750 ih2=1,k2
do 751 ih3=1,k3
ti3=ti3+mm1(j1,ih3)*xn2(ih3,ih2)*ny1(i1,ih2)
751 continue
do 752 ih4=1,k4
ti4=ti4+mm2(j1,ih2)*yn1(ih2,ih4)*ny2(i1,ih4)
752 continue
do 753 i2=1,n2
tj5=tj5+y1(i2,j1)*ny1(i2,ih2)*my1(i1,ih2)
tj6=tj6+y1(i2,j1)*my1(i2,ih2)*ny1(i1,ih2)
753 continue
750 continue
do 754 i2=1,n2
do 754 ih4=1,k4
tj7=tj7+y1(i2,j1)*ny2(i2,ih4)*my2(i1,ih4)
tj8=tj8+y1(i2,j1)*my2(i2,ih4)*ny2(i1,ih4)
754 continue
ccccccccccccccccccccccccccc
if(no.eq.1.or.no.eq.3) u31q=td3+chi2*(-td4+td5)
if(no.ge.5) u31q=td3+chi2*(td4+td5)

u20q=(chi2/2.0d0)*(te5+te6+te7+te8)
if(no.eq.1) u31j=fmi*(ti3+ti4)
if(no.ge.3) u31j=fmi*(-ti3+ti4)
u20j=(fmi/2.0d0)*(tj5+tj6-tj7-tj8)
c write(8,*)'h31..2.',u31q,u20q,u31j,u20j

cccccccccccccccccccccccccc
nm=nm+1
c write(*,*)'eeee',ichoix
go to(371,372,373,374,375,376,377,378,379),ichoix
371 H3(j1,nm)=u31q
go to 380
372 H3(j1,nm)=u31q+u20q
go to 380
373 H3(j1,nm)=u31q+u31j
go to 380
374 H3(j1,nm)=u31q+u31j+u20q
go to 380
375 H3(j1,nm)=u31q+u31j+u20j
go to 380
376 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 380
377 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 380
378 H3(j1,nm)=u31q+u31j+u20q+u20j
go to 380
379 H3(j1,nm)=u31q+u31j+u20q+u20j
c 379 H3(j1,nm)=u31q+u20j

380 continue
792 continue
if(nm.ne.n1+n2) stop 'error n m'
790 continue
ccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccc
ml=0

do 800 ih1=1,l2
do 808 ih2=1,ih1
v1=0.0d0
v2=0.0d0
v3=0.0d0
v4=0.0d0
v5=0.0d0
v6=0.0d0
v7=0.0d0
v8=0.0d0

do 801 j1=1,n1
v1=v1+fx2(j1,ih1)*fx2(j1,ih2)
v3=v3+gx2(j1,ih1)*gx2(j1,ih2)
if(no.eq.9)write(8,*)'1111',fx2(j1,ih1),gx2(j1,ih1),v1,v3

801 continue
do 802 j1=1,n2
v2=v2+fy1(j1,ih1)*fy1(j1,ih2)
v4=v4+gy1(j1,ih1)*gy1(j1,ih2)
if(no.eq.9)write(8,*)'22222',fy1(j1,ih1),gy1(j1,ih1),v2,v4

802 continue
do 803 j3=1,k3
v5=v5+nn1(ih1,j3)*nn1(ih2,j3)
v7=v7+mm1(ih1,j3)*mm1(ih2,j3)
if(no.eq.9)write(8,*)'=====...',nn1(ih1,j3),mm1(ih2,j3)

803 continue
do 804 j2=1,k2
v6=v6+nn2(ih1,j2)*nn2(ih2,j2)
v8=v8+mm2(ih1,j2)*mm2(ih2,j2)
if(no.eq.9)write(8,*)'......====',nn2(ih1,j2),mm2(ih2,j2)

804 continue
cccccccccccccccccccccccccccccccccccc
v0q=0.0d0
c if(ih1.eq.ih2) v0q=exy(ih1)+fdm
if(ih1.eq.ih2) v0q=exy(ih1)
v11q=(chi2/2.0d0)*(v1+v2)
v1q1=-(chi2/2.0d0)*(v3+v4)
v11j=-(fmi/2.0d0)*(v5+v6)
v1j1=(fmi/2.0d0)*(v7+v8)
cccccccccccccccccccccccccccccccc
write(8,*)'h300...',v0q,v11q,v1q1,v11j,v1j1
write(8,*)'fmi',fmi,'chi2',chi2
ml=ml+1
go to(381,382,383,384,385,386,387,388,389),ichoix
381 H0(ml)=v0q
go to 390
382 H0(ml)=v0q+v11q+v1q1
go to 390
383 H0(ml)=v0q
go to 390
384 H0(ml)=v0q+v11q+v1q1
go to 390
385 H0(ml)=v0q+v11j+v1j1
go to 390
386 H0(ml)=v0q+v11q+v1q1+v11j+v1j1
go to 390
387 H0(ml)=v0q+v11q+v1q1+v11j
go to 390
388 H0(ml)=v0q+v11q+v1q1+v11j+v1j1
go to 390
389 H0(ml)=v0q+v11q+v1q1+v11j+v1j1
c 389 H0(ml)=v0q+v11j+v1j1

390 CONTINUE
808 continue
800 continue

return
end
subroutine bat2(X,ndiag,lt,nq,mq,om2,omg1,omg2,prt,prt1,prt2,nmt,
*nmt1,nmt2,omtt1,ptt1,nmtt1,omtt2,ptt2,nmtt2,xf,xg,ndim1,w,chi2,
*rr,q2g,q2f,u,v,ga,sif,sig,sip,sim,siv1,siv2,siw1,siw2,fmi,qgs,
*gsj)
IMPLICIT REAL*8 (A-H,O-Z)
INTEGER OMG1,omg2,om2,PRT1,prt2,prt,omtt1,ptt1,omtt2,ptt2
real*8 ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,nn1,mm1,nn2,mm2
Dimension omg1(100),prt1(100),nmt1(100)
Dimension omg2(100),prt2(100),nmt2(100)
Dimension X(100),prt(20),nmt(20),om2(20)
Dimension ndiag(2),q2g(100),u(20),v(20),ga(2),q2f(100)
dimension omtt1(100),ptt1(100),nmtt1(100),omtt2(100),ptt2(100),
*nmtt2(100),xf(100),xg(100)
dimension x1(20,10),x2(20,10),f1(20,10),f2(20,10)
dimension n1(20,10),m1(20,10),n2(20,10),m2(20,10),x4(20,10)
dimension g1(20,10),g2(20,10),xy1(20,10),fy1(20,10),gy1(20,10)
dimension nn(20,10),mm(20,10),ny1(20,10),my1(20,10),xy2(20,10)
character*1 plus,moins,blanc,etoile,pop
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/ BIB04570
plus='+'
moins='-'
blanc=' '
etoile='*'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
sif=0.0d0
pop=plus
bp1=0.0d0
60 do 20 k=1,15,2
if(k.eq.1) i3=0
nox=k
noy=k+4
i1=0
bp2=0.0d0
do 1 j1=1,nq
c write(*,*) prt(j1),'pop',pop,ichar(pop),char(prt(j1))
if(char(prt(j1)).ne.pop) go to 1
if(om2(j1).ne.nox) go to 1
i1=i1+1
if(k.eq.1) i3=i3+1
l1=0
l2=0
if(k.eq.1) go to 10
do 2 j2=1,ndiag(lt)
if(omg1(j2).ne.om2(j1)) go to 3
if(prt1(j2).ne.prt(j1)) go to 3
if(nmt1(j2).ne.nmt(j1)) go to 3
l1=l1+1
x1(i1,l1)=x(j2)
g1(i1,l1)=q2g(j2)
f1(i1,l1)=q2f(j2)
ccccccccccccccccccccccccccccccccccccccccc
kk=k-2
if(k.ne.5) go to 4
i2=0
do 12 j8=1,ndiag(lt)
if(omg2(j8).ne.omg2(j2)) go to 12
if(prt2(j8).ne.prt2(j2)) go to 12
if(nmt2(j8).ne.nmt2(j2)) go to 12
if(omg1(j8).ne.kk) go to 12
i2=i2+1
xy1(i2,l1)=x(j8)
gy1(i2,l1)=q2g(j8)
fy1(i2,l1)=q2f(j8)
12 continue
4 i4=0
do 14 j10=1,mq
if(omtt2(j10).ne.omg2(j2)) go to 19
if(ptt2(j10).ne.prt2(j2)) go to 19
if(nmtt2(j10).ne.nmt2(j2)) go to 19
if(omtt1(j10).ne.kk) go to 19
i4=i4+1
ny1(i4,l1)=xf(j10)
my1(i4,l1)=xg(j10)
c write(8,*)k,j10,'qqqqq',ny1(i4,l1),xf(j10),my1(i4,l1),xg(j10)

l6=0
do 15 j11=1,ndiag(lt)
if(omg2(j11).ne.omtt1(j10)) go to 15
if(prt2(j11).ne.ptt1(j10)) go to 15
if(nmt2(j11).ne.nmtt1(j10)) go to 15
if(omg1(j11).eq.k) go to 15
l6=l6+1
xy2(l6,i4)=x(j11)
15 continue
go to 14
19 if(k.ne.3) go to 14
if(omtt1(j10).ne.omg2(j2)) go to 14
if(ptt1(j10).ne.prt2(j2)) go to 14
if(nmtt1(j10).ne.nmt2(j2)) go to 14
if(nmtt1(j10).eq.nmtt2(j10)) go to 14
i4=i4+1
ny1(i4,l1)=-xf(j10)
my1(i4,l1)=xg(j10)
14 continue
do 16 ik=1,nq
if(om2(ik).ne.omg2(j2)) go to 16
if(prt(ik).ne.prt2(j2)) go to 16
if(nmt(ik).ne.nmt2(j2)) go to 16
bp1=bp1+(U(j1)*v(j1)+u(ik)*v(ik))*(x(j2)**2)
bp2=bp2+u(j1)*v(j1)*u(ik)*v(ik)*(x(j2)**2)
16 continue

go to 2
cccccccccccccccccccccccccccccccccccccccccccccc
3 if(omg2(j2).ne.om2(j1)) go to 2
if(prt2(j2).ne.prt(j1)) go to 2
if(nmt2(j2).ne.nmt(j1)) go to 2
l2=l2+1
x2(i1,l2)=x(j2)
g2(i1,l2)=q2g(j2)
f2(i1,l2)=q2f(j2)
ccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccc
2 continue
ccccccccccccccccccccccccccccccccccccccccccccc
10 l3=0
l4=0
if(k.eq.1) l5=0
do 84 j4=1,mq
if(omtt1(j4).ne.om2(j1)) go to 85
if(ptt1(j4).ne.prt(j1)) go to 85
if(nmtt1(j4).ne.nmt(j1)) go to 85
if(k.eq.1) go to 87
l3=l3+1
n1(i1,l3)=xf(j4)
m1(i1,l3)=xg(j4)
87 continue
if(k.ne.1) go to 85
l5=l5+1
nn(i3,l5)=xf(j4)
mm(i3,l5)=xg(j4)
85 if(omtt2(j4).ne.om2(j1)) go to 84
if(ptt2(j4).ne.prt(j1)) go to 84
if(nmtt2(j4).ne.nmt(j1)) go to 84
if(k.eq.1) go to 86
l4=l4+1
n2(i1,l4)=xf(j4)
m2(i1,l4)=xg(j4)
86 if(k.ne.1) go to 84
if(omtt1(j4).ne.omtt2(j4)) go to 84
if(nmtt1(j4).eq.nmtt2(j4)) go to 84
l5=l5+1
nn(i3,l5)=-xf(j4)
mm(i3,l5)=xg(j4)
84 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccc
c write(*,*)'kil',k,i3,l5
1 continue
c write(*,*)k,i3,l5
if(l5.ne.i3) stop'error'
c write(*,*)'eeeeeeeee2222'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
z1=0.0d0
z2=0.0d0
z3=0.0d0

if(no.ne.1) go to 200
do 201 j1=1,i3
do 201 j2=1,l5
z1=z1+nn(j1,j2)**2
201 continue
gsj=gsj+0.50d0*z1
200 if(no.eq.1) go to 210
do 202 j1=1,i1
do 203 j2=1,l3
z2=z2+n1(j1,j2)**2
203 continue
do 204 j2=1,l1
z3=z3+f1(j1,j2)**2
204 continue
202 continue
qgs=qgs-chi2*z3
gsj=gsj+fmi*z2
210 continue
cccccccccccccccccccccccccccccccccccccccccccccc
ta1=0.0d0
ta2=0.0d0
ta3=0.0d0
ta4=0.0d0
tb1=0.0d0
tb2=0.0d0
tb3=0.0d0
tb4=0.0d0
tv1=0.0d0
if(k.eq.1) go to 20
do 21 j1=1,i1
do 21 j2=1,i1
do 22 k1=1,l1
do 23 k2=1,l1
ta1=ta1+x1(j1,k1)*x1(j2,k1)*f1(j1,k2)*f1(j2,k2)
ta2=ta2+x1(j1,k1)*x1(j1,k2)*f1(j2,k1)*f1(j2,k2)
tb1=tb1+x1(j1,k1)*x1(j2,k1)*g1(j1,k2)*g1(j2,k2)
tb2=tb2+x1(j1,k1)*x1(j1,k2)*g1(j2,k1)*g1(j2,k2)
23 continue
do 24 k3=1,l2
ta3=ta3+x1(j1,k1)*x1(j2,k1)*f2(j1,k3)*f2(j2,k3)
ta4=ta4+f1(j1,k1)*f1(j2,k1)*x2(j1,k3)*x2(j2,k3)
tb3=tb3+x1(j1,k1)*x1(j2,k1)*g2(j1,k3)*g2(j2,k3)
tb4=tb4+g1(j1,k1)*g1(j2,k1)*x2(j1,k3)*x2(j2,k3)
tv1=tv1+x1(j1,k1)*g1(j2,k1)*g2(j1,k3)*x2(j2,k3)
24 continue
22 continue
21 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccc
tc1=0.0d0
tc2=0.0d0
tc3=0.0d0
tc4=0.0d0
td1=0.0d0
td2=0.0d0
td3=0.0d0
td4=0.0d0
do 30 j1=1,i1
do 30 j2=1,i1
do 31 k1=1,l1
do 32 k2=1,l3
tc1=tc1+x1(j1,k1)*x1(j2,k1)*n1(j1,k2)*n1(j2,k2)
td1=td1+x1(j1,k1)*x1(j2,k1)*m1(j1,k2)*m1(j2,k2)
32 continue
do 33 k3=1,l4
tc3=tc3+x1(j1,k1)*x1(j2,k1)*n2(j1,k3)*n2(j2,k3)
td3=td3+x1(j1,k1)*x1(j2,k1)*m2(j1,k3)*m2(j2,k3)
33 continue
31 continue
do 34 k4=1,l2
do 36 k6=1,l3
tc4=tc4+x2(j1,k4)*x2(j2,k4)*n1(j1,k6)*n1(j2,k6)
td4=td4+x2(j1,k4)*x2(j2,k4)*m1(j1,k6)*m1(j2,k6)
36 continue
34 continue
30 continue

if(k.lt.5) go to 39
do 37 j1=1,i1
do 37 k1=1,l1
do 37 k2=1,l1
do 37 j2=1,i4
tc2=tc2+x1(j1,k1)*x1(j1,k2)*ny1(j2,k2)*ny1(j2,k1)
td2=td2+x1(j1,k1)*x1(j1,k2)*my1(j2,k2)*my1(j2,k1)
37 continue
39 continue
cccccccccccccccccccccccccccccccccccccccccccccccccc
te1=0.0d0
te2=0.0d0
tf1=0.0d0
tf2=0.0d0
tv2=0.0d0
if(k.ne.5) go to 41
do 40 j1=1,i1
do 40 j2=1,i2
do 40 k1=1,l1
do 40 k2=1,l1
te1=te1+x1(j1,k1)*x1(j1,k2)*fy1(j2,k1)*fy1(j2,k2)
tf1=tf1+x1(j1,k1)*x1(j1,k2)*gy1(j2,k1)*gy1(j2,k2)
te2=te2+xy1(j2,k1)*xy1(j2,k2)*f1(j1,k1)*f1(j1,k2)
tf2=tf1+xy1(j2,k1)*xy1(j2,k2)*g1(j1,k1)*g1(j1,k2)
tv2=tv2+xy1(j2,k1)*gy1(j2,k2)*x1(j1,k2)*g1(j1,k1)
40 continue
41 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
tg1=0.0d0
tg2=0.0d0
th1=0.0d0
th2=0.0d0
tv4=0.0d0
if(k.ne.5) go to 57
write(*,*)'l5',l5,'l1',l1
if(l5.ne.l1) stop 'error 51'
do 50 j1=1,i1
do 50 k1=1,l1
do 51 k2=1,l1
do 53 j2=1,i3
tg2=tg2+x1(j1,k1)*x1(j1,k2)*nn(j2,k1)*nn(j2,k2)
th2=th2+x1(j1,k1)*x1(j1,k2)*mm(j2,k1)*mm(j2,k2)
53 continue
51 continue
if(l1.ne.l5) stop 'error 52'
if(l3.ne.i2) stop 'error 888'
do 54 k2=1,l3
do 54 j4=1,l1
tv4=tv4+x1(j1,k1)*m1(j1,k2)*xy1(k2,j4)*mm(j4,k1)
54 continue
50 continue
57 continue
ccccccccccccccccccccc
tg3=0.0d0
tg4=0.0d0
th3=0.0d0
th4=0.0d0
tv5=0.0d0
if(k.ne.3) go to 67
if(l3.ne.l1.and.l3.ne.0.and.l1.ne.0) stop 'error 31'
if(l1.eq.0) l1=l3
if(l5.ne.l1.and.l5.ne.0.and.l1.ne.0) stop 'error 32'

do 65 j1=1,i1
do 65 k1=1,l1
do 61 k2=1,l1
do 62 j2=1,i1
tg3=tg3+x1(j1,k1)*x1(j1,k2)*n1(j2,k1)*n1(j2,k2)
th3=th3+x1(j1,k1)*x1(j1,k2)*m1(j2,k1)*m1(j2,k2)
62 continue
do 63 j2=1,l5
tg4=tg4+x1(j1,k1)*x1(j1,k2)*nn(k2,j2)*nn(k1,j2)
th4=th4+x1(j1,k1)*x1(j1,k2)*mm(k2,j2)*mm(k1,j2)
63 continue
61 continue
do 64 j5=1,i1
do 64 k2=1,l1
tv5=tv5+x1(j1,k1)*m1(j1,k2)*m1(j5,k1)*x1(j5,k2)
64 continue
65 continue
67 continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
tv3=0.0d0
if(k.lt.5) go to 77
if(l6.ne.l4.and.l6.ne.0) stop 'error 46'
do 70 j1=1,i4
do 70 j2=1,i1
do 70 k1=1,l1
do 70 k2=1,l4
tv3=tv3+my1(j1,k1)*xy2(k2,j1)*m2(j2,k2)*x1(j2,k1)
70 continue
77 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
siv1=siv1+(chi2/2.0d0)*(ta1+ta2+ta3+ta4+te1+te2)
siv2=siv2-(chi2/2.0d0)*(tb1+tb2+tb3+tb4+tf1+tf2)
siw1=siw1-(fmi/2.0d0)*(tc1+tc2+tc3+tc4+tg2+tg3+tg4)
siw2=siw2+(fmi/2.0d0)*(td1+td2+td3+td4+th2+th3+th4)
phas=1.0d0
if(k.ne.3) phas=-1.0d0
sim=sim+(2.0d0*chi2)*(phas*tv1+tv2)
sip=sip-(2.0d0*fmi)*(tv3-tv4+0.5d0*tv5)
sig=sig-(2.0d0*ga(lt))*bp2

c write(*,*)'k',k,pop

20 continue
if(pop.eq.moins) go to 90
pop=moins
go to 60
90 continue
sif=2.0d0*bp1
return
end
subroutine bat1(X,pop,noy,nox,no,ndiag,lt,nq,om2,omg1,omg2,
*prt,prt1,prt2,nmt,nmt1,nmt2,omt,prtt,nmtt,nct,nctt,omtt1,
*ptt1,nmtt1,omtt2,ptt2,nmtt2,xf,xg,mq,a,eq1,eq2,eqq,ndim1,n,nm,
*H1,w,chi2,rr,q2g,q2f,mt,rtt,mtt,H3,m2,m4,exy,sif,sig,sim,sip,
*siv1,siv2,siw1,siw2,u,v,ga,ichoix,fmi,H0,ml,xj0,xj1,xj2,xj3,
*yj0,yj1,yj2,yj3,kk,k2,n1,n2,fdm)

IMPLICIT REAL*8 (A-H,O-Z)
INTEGER OMG1,omg2,om2,PRT1,prt2,prt,prtt,omt,rtt,omy,prty,omtt1,
*omtt2,ptt1,ptt2
real*8 ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,ny1,my1,ny2,my2,ny4,my4,
*nx4,mx4,nn1,mm1,nn2,mm2

Dimension omg1(100),prt1(100),nmt1(100)
Dimension omg2(100),prt2(100),nmt2(100)
dimension omtt1(100),ptt1(100),nmtt1(100),xf(100),xg(100)
dimension omtt2(100),ptt2(100),nmtt2(100)
dimension omy(10),prty(10),nmy(10)
Dimension X(100),eq1(100),eq2(100),prt(20),nmt(20),om2(20),
*nct(60)
Dimension nctt(60),ndiag(2),eqq(20),q2g(100),q2f(100),u(20),
*v(20),ga(2)
dimension A(100),omt(100),prtt(100),nmtt(100)
dimension H1(100),H3(20,20),mt(60),rtt(60),mtt(60),H0(20)
dimension x1(10,10),gx1(10,10),fx1(10,10),x2(10,10),gx2(10,10),
*fx2(10,10),x3(10,10),gx3(10,10),fx3(10,10)
dimension ex(20),ey(20),ex1(20,10),ey1(20,10),exy(20),ny3(10,10),
*my3(10,10),nx3(10,10),mx3(10,10),nx2(10,10),mx2(10,10)
dimension mx1(10,10),nx1(10,10),xn1(10,10),xn2(10,10),yn2(10,10)
dimension y1(10,10),gy1(10,10),fy1(10,10),y2(10,10),gy2(10,10),
*fy2(10,10),ny1(10,10),my1(10,10),ny2(10,10),my2(10,10)
dimension ny4(10,10),my4(10,10),nx4(10,10),mx4(10,10),yn1(10,10),
*y33(10,10),nn1(10,10),nn2(10,10),mm1(10,10),mm2(10,10)
dimension uk(10),uh(10),ux(10,10),uy(10,10),uxy(10,10)
dimension xj0(10,10),xj1(10,10),xj2(10,10),xj3(10,10),yj0(10,10),
*yj1(10,10),yj2(10,10),yj3(10,10)
character*1 plus,moins,blanc,etoile,pop
c DATA PLUS/1H+/,MOINS/1H-/,BLANC/1H /,ETOILE/1H*/ BIB04570
plus='+'
moins='-'
blanc=' '
etoile='*'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
n1=0
n2=0
n3=0
l4=0
do 1000 i=1,10
do 1001 j=1,10
x1(i,j)=0.0d0
gx1(i,j)=0.0d0
fx1(i,j)=0.0d0
x2(i,j)=0.0d0
gx2(i,j)=0.0d0
fx2(i,j)=0.0d0
y1(i,j)=0.0d0
gy1(i,j)=0.0d0
fy1(i,j)=0.0d0
y2(i,j)=0.0d0
gy2(i,j)=0.0d0
fy2(i,j)=0.0d0
x3(i,j)=0.0d0
gx3(i,j)=0.0d0
fx3(i,j)=0.0d0
ny3(i,j)=0.0d0
my3(i,j)=0.0d0
nx3(i,j)=0.0d0
mx3(i,j)=0.0d0
nx2(i,j)=0.0d0
mx2(i,j)=0.0d0
nx1(i,j)=0.0d0
mx1(i,j)=0.0d0
xn1(i,j)=0.0d0
xn2(i,j)=0.0d0
yn2(i,j)=0.0d0
ny1(i,j)=0.0d0
my1(i,j)=0.0d0
ny2(i,j)=0.0d0
my2(i,j)=0.0d0
ny4(i,j)=0.0d0
my4(i,j)=0.0d0
nx4(i,j)=0.0d0
mx4(i,j)=0.0d0
yn1(i,j)=0.0d0
y33(i,j)=0.0d0
nn1(i,j)=0.0d0
mm1(i,j)=0.0d0
nn2(i,j)=0.0d0
mm2(i,j)=0.0d0
ux(i,j)=0.0d0
uy(i,j)=0.0d0
uxy(i,j)=0.0d0
1001 continue
do 1002 j=1,20
ex1(j,i)=0.0d0
ey1(j,i)=0.0d0
1002 continue
uk(i)=0.0d0
uh(i)=0.0d0
1000 continue
do 1003 i=1,20
ex(i)=0.0d0
ey(i)=0.0d0
exy(i)=0.0d0
1003 continue
cccccccccccccccccccccccccccccc
l1=0
l2=0
l3=0
l4=0
l5=0
l10=0
l11=0
l12=0
l13=0
k1=0
k2=0
k3=0
k4=0
k5=0
k6=0
k7=0
k8=0
k10=0
k20=0
k30=0
k31=0
n1=0
n2=0
n3=0

do 1 j1=1,nq
if(char(prt(j1)).ne.pop) go to 1
if(om2(j1).ne.nox) go to 6
n1=n1+1
l1=0
l2=0
k3=0
k1=0
Ex(n1)=eqq(j1)
omt(n1)=om2(j1)
prtt(n1)=prt(j1)
nmtt(n1)=nmt(j1)
nct(n1)=2
nctt(n1)=no
uk(n1)=u(j1)*v(j1)
do 2 j2=1,ndiag(lt)
if(omg1(j2).ne.om2(j1)) go to 3
if(prt1(j2).ne.prt(j1)) go to 3
if(nmt1(j2).ne.nmt(j1)) go to 3
uv=0.0d0
do 117 k=1,nq
if(om2(k).ne.omg2(j2)) go to 117
if(prt(k).ne.prt2(j2)) go to 117
if(nmt(k).ne.nmt2(j2)) go to 117
uv=u(k)*v(k)
117 continue
if(omg2(j2).ne.no) go to 4
l2=l2+1
x2(n1,l2)=x(j2)
exy(l2)=eq2(j2)
gx2(n1,l2)=q2g(j2)
fx2(n1,l2)=q2f(j2)
uxy(n1,l2)=uv
mt(l2)=omg2(j2)
rtt(l2)=prt2(j2)
mtt(l2)=nmt2(j2)
4 if(no.lt.7) go to 3
l1=l1+1
k5=0
l11=0
x1(n1,l1)=x(j2)
gx1(n1,l1)=q2g(j2)
ex1(n1,l1)=eq2(j2)
fx1(n1,l1)=q2f(j2)
ux(n1,l1)=uv
cccccccccccccccccccccccccccccccccccccc
do 500 i1=1,mq
if(omtt1(i1).ne.omg2(j2)) go to 502
if(ptt1(i1).ne.prt2(j2)) go to 502
if(nmtt1(i1).ne.nmt2(j2)) go to 502
if(no.ge.9) go to 501
l11=l11+1
ny3(l1,l11)=xf(i1)
my3(l1,l11)=xg(i1)
go to 502
501 k5=k5+1
nx3(l1,k5)=xf(i1)
mx3(l1,k5)=xg(i1)
502 if(omtt2(i1).ne.omg2(j2)) go to 500
if(ptt2(i1).ne.prt2(j2)) go to 500
if(nmtt2(i1).ne.nmt2(j2)) go to 500
if(omtt1(i1).eq.omtt2(i1)) go to 506
if(no.ge.9) go to 505
k5=k5+1
nx3(l1,k5)=xf(i1)
mx3(l1,k5)=xg(i1)
go to 500
505 l11=l11+1
ny3(l1,l11)=xf(i1)
my3(l1,l11)=xg(i1)
go to 500
506 if(nmtt1(i1).eq.nmtt2(i1)) go to 500
if(no.ne.7) go to 507
l11=l11+1
ny3(l1,l11)=-xf(i1)
my3(l1,l11)=xg(i1)
507 if(no.ne.9) go to 500
k5=k5+1
nx3(l1,k5)=-xf(i1)
mx3(l1,k5)=xg(i1)
500 continue

ccccccccccccccccccccccccccccccccccccccccccccccc
l5=0
if(no.ge.11) go to 13
do 12 j4=1,ndiag(lt)
if(omg2(j4).ne.omg2(j2)) go to 12
if(prt2(j4).ne.prt2(j2)) go to 12
if(nmt2(j4).ne.nmt2(j2)) go to 12
if(omg1(j4).eq.nox) go to 12
l5=l5+1
x3(l1,l5)=x(j4)
gx3(l1,l5)=q2g(j4)
fx3(l1,l5)=q2f(j4)
12 continue
13 do 14 j4=1,ndiag(lt)
if(omg1(j4).ne.omg2(j2)) go to 14
if(prt1(j4).ne.prt2(j2)) go to 14
if(nmt1(j4).ne.nmt2(j2)) go to 14
l5=l5+1
x3(l1,l5)=x(j4)
gx3(l1,l5)=q2g(j4)
fx3(l1,l5)=q2f(j4)
14 continue
ccccccccccccccccccccccccccccccccccccccccccccccc
3 if(omg2(j2).ne.om2(j1)) go to 2
if(prt2(j2).ne.prt(j1)) go to 2
if(nmt2(j2).ne.nmt(j1)) go to 2
uv=0.0d0
do 929 k=1,nq
if(om2(k).ne.omg1(j2)) go to 929
if(prt(k).ne.prt1(j2)) go to 929
if(nmt(k).ne.nmt1(j2)) go to 929
uv=u(k)*v(k)
929 continue
if(no.ge.7) go to 5
if(omg1(j2).eq.no) go to 5
l1=l1+1
x1(n1,l1)=x(j2)
gx1(n1,l1)=q2g(j2)
ex1(n1,l1)=eq1(j2)
fx1(n1,l1)=q2f(j2)
ux(n1,l1)=uv
cccccccccccccccccccccccccccccccccccccc
k5=0
l11=0
do 550 i1=1,mq
if(omtt2(i1).ne.omg1(j2)) go to 551
if(ptt2(i1).ne.prt1(j2)) go to 551
if(nmtt2(i1).ne.nmt1(j2)) go to 551
k5=k5+1
nx3(l1,k5)=xf(i1)
mx3(l1,k5)=xg(i1)
551 if(omtt1(i1).ne.omg1(j2)) go to 550
if(ptt1(i1).ne.prt1(j2)) go to 550
if(nmtt1(i1).ne.nmt1(j2)) go to 550
l11=l11+1
ny3(l1,l11)=xf(i1)
my3(l1,l11)=xg(i1)
c write(8,*)'xxxxxxx',l1,l11,ny3(l1,l11),xf(i1)
550 continue
ccccccccccccccccccccccccccccccccccccccccccccccc
l5=0
do 16 j4=1,ndiag(lt)
if(omg2(j4).ne.omg1(j2)) go to 16
if(prt2(j4).ne.prt1(j2)) go to 16
if(nmt2(j4).ne.nmt1(j2)) go to 16
l5=l5+1
x3(l1,l5)=x(j4)
gx3(l1,l5)=q2g(j4)
fx3(l1,l5)=q2f(j4)
16 continue
cccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccc
5 if(no.eq.1) go to 2
if(omg1(j2).ne.no) go to 2
l2=l2+1
x2(n1,l2)=x(j2)
exy(l2)=eq1(j2)
gx2(n1,l2)=q2g(j2)
fx2(n1,l2)=q2f(j2)
uxy(n1,l2)=uv
mt(l2)=omg1(j2)
rtt(l2)=prt1(j2)
mtt(l2)=nmt1(j2)
2 continue
cccccccccccccccccccccccccccccccccccccccccccccccc
do 400 j2=1,mq
if(omtt1(j2).ne.om2(j1)) go to 402
if(ptt1(j2).ne.prt(j1)) go to 402
if(nmtt1(j2).ne.nmt(j1)) go to 402
if(no.ge.5) go to 401
k3=k3+1
k7=0
nx2(n1,k3)=xf(j2)
mx2(n1,k3)=xg(j2)
go to 410
401 k1=k1+1
k8=0
k20=0
nx1(n1,k1)=xf(j2)
mx1(n1,k1)=xg(j2)
ccccccccccccccccccccccccccccc
410 do 404 j4=1,ndiag(lt)
if(omg2(j4).ne.omtt2(j2)) go to 422
if(prt2(j4).ne.ptt2(j2)) go to 422
if(nmt2(j4).ne.nmtt2(j2)) go to 422
no1=omg1(j4)
if(no.eq.1) no1=3
if(no.eq.3) no1=5
if(no.eq.5) no1=3
if(no.eq.7) no1=5
if(no.ge.5) go to 405
if(omg1(j4).ne.no1) go to 404
k7=k7+1
xn2(k3,k7)=x(j4)
C write(8,*)'999999999',k3,k7,x(j4)
go to 404
405 if(omg1(j4).ne.no1) go to 421
k20=k20+1
xn1(k1,k20)=x(j4)
go to 404
421 if(no.ge.9) go to 422
k8=k8+1
yn2(k1,k8)=x(j4)
C write(8,*)'77777777777...',k1,k8,x(j4)
go to 404
422 if(omg1(j4).ne.omtt2(j2)) go to 404
if(prt1(j4).ne.ptt2(j2)) go to 404
if(nmt1(j4).ne.nmtt2(j2)) go to 404
k8=k8+1
yn2(k1,k8)=x(j4)
C write(8,*)'77777777...777',k1,k8,x(j4)
404 continue
go to 400
ccccccccccccccccccccccccccc
402 if(omtt2(j2).ne.om2(j1)) go to 400
if(ptt2(j2).ne.prt(j1)) go to 400
if(nmtt2(j2).ne.nmt(j1)) go to 400
if(omtt2(j2).eq.omtt1(j2)) go to 430
if(no.lt.5) go to 403
k3=k3+1
k7=0
nx2(n1,k3)=xf(j2)
mx2(n1,k3)=xg(j2)
go to 411
403 k1=k1+1
k8=0
k20=0
nx1(n1,k1)=xf(j2)
mx1(n1,k1)=xg(j2)
cccccccccccccccccccccccccc
if(no.ge.5) go to 411
do 493 j3=1,ndiag(lt)
if(omg1(j3).ne.omtt1(j2)) go to 493
if(prt1(j3).ne.ptt1(j2)) go to 493
if(nmt1(j3).ne.nmtt1(j2)) go to 493
k20=k20+1
xn1(k1,k20)=x(j3)
493 continue
ccccccccccccccccccccccccccc
411 do 406 j5=1,ndiag(lt)
if(omg2(j5).ne.omtt1(j2)) go to 406
if(prt2(j5).ne.ptt1(j2)) go to 406
if(nmt2(j5).ne.nmtt1(j2)) go to 406
if(no.lt.5) go to 407
k7=k7+1
xn2(k3,k7)=x(j5)
C write(8,*)'99...9999999',k3,k7,x(j5)
go to 406
407 k8=k8+1
yn2(k1,k8)=x(j5)
C write(8,*)'77777777777',k1,k8,x(j5)
406 continue
go to 400
cccccccccccccccccccccccccc
430 if(nmtt2(j2).eq.nmtt1(j2)) go to 400
if(no.ne.3) go to 431
k3=k3+1
nx2(n1,k3)=-xf(j2)
mx2(n1,k3)=xg(j2)
431 if(no.ne.5) go to 400
k1=k1+1
nx1(n1,k1)=-xf(j2)
mx1(n1,k1)=xg(j2)
cccccccccccccccccccccccc
400 continue
ccccccccccccccccccccccccccccccccccccccc
6 if(om2(j1).ne.noy) go to 199
n2=n2+1
l3=0
l4=0
ey(n2)=eqq(j1)
uh(n2)=u(j1)*v(j1)
omy(n2)=om2(j1)
prty(n2)=prt(j1)
nmy(n2)=nmt(j1)
do 7 j3=1,ndiag(lt)
if(omg1(j3).ne.om2(j1)) go to 8
if(prt1(j3).ne.prt(j1)) go to 8
if(nmt1(j3).ne.nmt(j1)) go to 8
if(omg2(j3).ne.no) go to 8
l4=l4+1
y1(n2,l4)=x(j3)
gy1(n2,l4)=q2g(j3)
fy1(n2,l4)=q2f(j3)
8 if(omg2(j3).ne.om2(j1)) go to 7
if(prt2(j3).ne.prt(j1)) go to 7
if(nmt2(j3).ne.nmt(j1)) go to 7
l3=l3+1
y2(n2,l3)=x(j3)
ey1(n2,l3)=eq1(j3)
gy2(n2,l3)=q2g(j3)
fy2(n2,l3)=q2f(j3)
ccccccccccccccccccccccccccccccccccccccccccccccccccccc
l12=0
l13=0
do 560 i1=1,mq
if(omtt1(i1).ne.omg1(j3)) go to 561
if(ptt1(i1).ne.prt1(j3)) go to 561
if(nmtt1(i1).ne.nmt1(j3)) go to 561
l12=l12+1
nx4(l3,l12)=xf(i1)
mx4(l3,l12)=xg(i1)
c write(8,*)'222222',l3,l12,xf(i1)
561 if(omtt2(i1).ne.omg1(j3)) go to 560
if(ptt2(i1).ne.prt1(j3)) go to 560
if(nmtt2(i1).ne.nmt1(j3)) go to 560
l13=l13+1
ny4(l3,l13)=xf(i1)
my4(l3,l13)=xg(i1)
560 continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
do 939 k=1,nq
if(om2(k).ne.omg1(j3)) go to 939
if(prt(k).ne.prt1(j3)) go to 939
if(nmt(k).ne.nmt1(j3)) go to 939
uy(n2,l3)=u(k)*v(k)
939 continue
c y3(...)=0.0d0
c gy3(...)=0.0d0
c fy3(....)=0.0d0
7 continue
cccccccccccccccccccccccccccccccccccccccccccccc
k2=0
k4=0

do 600 i1=1,mq
if(omtt1(i1).ne.om2(j1)) go to 601
if(ptt1(i1).ne.prt(j1)) go to 601
if(nmtt1(i1).ne.nmt(j1)) go to 601
k2=k2+1
ny1(n2,k2)=xf(i1)
my1(n2,k2)=xg(i1)
C write(8,*)'22222',n2,k2,xf(i1)

ccccccccccccccccccccccccccc
k10=0
do 602 i2=1,ndiag(lt)
if(omg2(i2).ne.omtt2(i1)) go to 602
if(prt2(i2).ne.ptt2(i1)) go to 602
if(nmt2(i2).ne.nmtt2(i1)) go to 602
k10=k10+1
yn1(k2,k10)=x(i2)
602 continue
cccccccccccccccccccccccccccc
601 if(omtt2(i1).ne.om2(j1)) go to 600
if(ptt2(i1).ne.prt(j1)) go to 600
if(nmtt2(i1).ne.nmt(j1)) go to 600
k4=k4+1
ny2(n2,k4)=xf(i1)
my2(n2,k4)=xg(i1)
ccccccccccccccccccccccccccccc
k6=0
do 603 i2=1,ndiag(lt)
if(omg2(i2).ne.omtt1(i1)) go to 603
if(prt2(i2).ne.ptt1(i1)) go to 603
if(nmt2(i2).ne.nmtt1(i1)) go to 603
write(*,*)'....',no,'k6',k6
k6=k6+1
y33(k4,k6)=x(i2)
603 continue
600 continue
ccccccccccccccccccccccccccccccccccccccccc
199 if(om2(j1).ne.no) go to 1
n3=n3+1
k30=0
k31=0
do 700 i1=1,mq
if(omtt1(i1).ne.om2(j1)) go to 701
if(ptt1(i1).ne.prt(j1)) go to 701
if(nmtt1(i1).ne.nmt(j1)) go to 701
k30=k30+1
nn1(n3,k30)=xf(i1)
mm1(n3,k30)=xg(i1)
c write(8,*)'nnnnn',nn1(n3,k30),xf(i1),xg(i1)
701 if(omtt2(i1).ne.om2(j1)) go to 700
if(ptt2(i1).ne.prt(j1)) go to 700
if(nmtt2(i1).ne.nmt(j1)) go to 700
if(omtt2(i1).eq.omtt1(i1)) go to 702
k31=k31+1
nn2(n3,k31)=xf(i1)
mm2(n3,k31)=xg(i1)
go to 700
702 if(no.ne.1) go to 700
if(nmtt2(i1).eq.nmtt1(i1)) go to 700
k30=k30+1
nn1(n3,k30)=-xf(i1)
mm1(n3,k30)=xg(i1)
c write(8,*)'n...nnnn',nn1(n3,k30),xf(i1),xg(i1)

700 continue
ccccccccccccccccccccccccccccccccccccccccccccc
1 continue
ccccccccccccccccccccccccccccccccccccccc
m2=l2
m4=l4
c write(*,*)l1,l2,l3,l4
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c write(8,*)'k6=',k6,'l13=',l13
c write(*,*)n1,n2,n3
c write(*,*)l1,l2,l3,l4,l5,l11,l12,l13
c write(*,*)k1,k2,k3,k4,k5,k6,k7,k8,k10,k20,k30,k31
if(l4.eq.0.or.l2.eq.0) go to 9
IF(l2.ne.l4) stop 'error1'
9 continue
if(l2.ne.n3.and.n3.ne.0.and.l2.ne.0) stop 'error 2'
if(l2.eq.0) l2=n3
if(k1.ne.l11.and.l11.ne.0.and.k1.ne.0) stop 'error 3'
if(k1.eq.0) k1=l11
if(k2.ne.k7.and.k2.ne.0.and.k7.ne.0) stop 'error 4'
if(k2.ne.k31.and.k31.ne.0.and.k2.ne.0) stop 'error 8'
if(k7.ne.k31.and.k31.ne.0.and.k7.ne.0) stop 'error 12'
if(k2.eq.0.and.k7.ne.0) k2=k7
if(k2.eq.0.and.k31.ne.0) k2=k31
if(k3.ne.k20.and.k20.ne.0.and.k3.ne.0) stop 'error 5'
if(k3.ne.k30.and.k30.ne.0.and.k3.ne.0) stop 'error 6'
if(k4.ne.l12.and.l12.ne.0.and.k4.ne.0) stop 'error 7'

if(k4.ne.k10.and.k10.ne.0) stop 'error 9'
if(k5.ne.k8.and.k8.ne.0.and.k5.ne.0) stop 'error 10'
write(*,*)'n',no,'k6',k6,'l13',l13

if(k6.ne.l13.and.k6.ne.0.and.l13.ne.0) stop 'error 11'
c if(k6.eq.0) k6=l13
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
N=N1+N2
DO 10 i=n1+1,n
j=i-n1
omt(i)=omy(j)
prtt(i)=prty(j)
nmtt(i)=nmy(j)
nct(i)=2
nctt(i)=no
10 continue
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c write(8,*)'..........',k3,k30,k20,n3,l2
c do 2000 i=1,l2
c do 2000 j=1,k3
c write(8,*)'zzzzzz',mm1(i,j)
c 2000 continue
c stop
c m=0
c write(*,*)'call so4'
call str22(H1,m,a,n1,n2,l1,l2,l3,l5,k1,k2,k3,k4,k5,k6,x1,gx1,fx1,
*x2,gx2,fx2,x3,gx3,fx3,ex,exy,ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,
*xn1,xn2,yn2,y1,gy1,fy1,y2,gy2,fy2,ny1,my1,ny2,my2,ny4,my4,nx4,
*mx4,yn1,y33,nn1,nn2,mm1,mm2,uk,uh,ga,lt,ex1,ey1,ey,fmi,chi2,
*rr,W,no,ichoix,sif,sig,sim,sip,siv1,siv2,siw1,siw2,ux,uxy,uy)
c write(*,*)'call so5'
call str31(H3,nm,n1,n2,l1,l2,l3,k1,k2,k3,k4,k5,k6,x1,gx1,fx1,x2,
*gx2,fx2,x3,gx3,fx3,ex,exy,ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,xn1,
*xn2,yn2,y1,gy1,fy1,y2,gy2,fy2,ny1,my1,ny2,my2,ny4,my4,nx4,mx4,
*yn1,y33,nn1,nn2,mm1,mm2,uk,uh,ga,lt,fmi,chi2,rr,w,h0,ml,ex1,
*ey1,ey,ichoix,no,fdm,sif,sig,sim,sip,siv1,siv2,siw1,siw2)
kk=k3+k4
call coriol(xj0,xj1,xj2,xj3,yj0,yj1,yj2,yj3,kk,n1,n2,l1,l2,l3,
*k1,k2,k3,k4,k5,k6,x1,x2,x3,ny3,my3,nx3,mx3,nx2,mx2,mx1,nx1,xn1,
*xn2,yn2,y1,y2,ny1,my1,ny2,my2,ny4,my4,nx4,mx4,yn1,y33,nn1,nn2,
*mm1,mm2,fmi,ichoix,no)

return
end

**vmagnin** · 10/03/2023, 09h39

Effectivement, c'est un peu long :-)
Vous pouvez essayer d'éditer votre message pour effacer le code et mettre le fichier en pièce jointe (si l'extension .for n'est pas acceptée, il suffit de la remplacer par .txt). D'autant plus que c'est du FORTRAN 77, donc sans les indentations on ne peut pas compiler le programme, les six premières colonnes ayant un statut particulier.

Je pense que le message d'erreur provient de l'instruction STOP ligne 2315 :

cccccccccccccccccccccccccccccccccccccccccccccccccccc
c write(*,*)'kil',k,i3,l5
1 continue
c write(*,*)k,i3,l5
if(l5.ne.i3) stop'error'
c write(*,*)'eeeeeeeee2222'
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
z1=0.0d0
z2=0.0d0
z3=0.0d0

On a donc l5 différent de i3. Reste à savoir ce que sont ces variables.

On est là face à du vieux code spaghetti, plein de GOTO... Un style de programmation extrêmement difficile à déboguer...

**progrmmation** · 10/03/2023, 10h37

J'ai essayé de mettre le code en format txt mais ils disent que le format dépasse la limite du forume...

**anapurna** · 10/03/2023, 15h09

salut

c'est un code completement destructuré sans commentaire
il faut eradiquer les Goto qui von dans tout les sens essayer de les modifier par des appel de metode
cela permettra d'y voir deja nettement plus claire

pour cela tu as sois les subRoutine ou function

la form Do label Variable=PremiereValeur,increment j'avous que c'est un peut pertubant
on dois pouvoir la rempacer par un

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
     Do  Variable=PremiereValeur,dernierValeur
     .. 
     End Do

Ce qui peut t'aider c'est de respecter les cass des variable et des mots cle
exemple ce block

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
  N=N1+N2
  DO 10 i=n1+1,n
     j=i-n1
     omt(i)=omy(j)
     prtt(i)=prty(j)
     nmtt(i)=nmy(j)
     nct(i)=2
     nctt(i)=no
  10 continue

ce code aurais pu etre

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
  SubRoutine remplittab(n1,n2)
  n=n1+n2
  DO 10 i=n1+1,n
     j = i-n1
     omt(i)  =omy(j)
     prtt(i)  =prty(j)
     nmtt(i) =nmy(j)
     nct(i)    =2
     nctt(i)  =no
  return
 end
  .....
  Call remplittab(N1,N2)

la pour le coup je te souhaite bon courage pour demeler ce sac de noeud

**vmagnin** · 10/03/2023, 21h06

Il existe des outils pour refactoriser ce genre de vieux code en un Fortran plus moderne :
https://github.com/Beliavsky/Fortran-Tools#refactoring

**progrmmation** · 11/03/2023, 15h53

Malheureusement je n'ai pas encore arriver à résoudre ce problème
En tous cas je vous remercie pour vos idées intéressantes

**vmagnin** · 11/03/2023, 16h50

Bon courage ! Vu la quasi absence de commentaires dans le code, j'ai l'impression qu'il faudrait au moins deux mois de travail pour comprendre comment ça fonctionne, quelles sont les méthodes numériques employées, etc.