fortran crash

[Tony1984]

hey all
i have a problem in fortran,
i just wrote a code but it not working properly, it suddenly stop or crashed in some point and i cant figure it out.
please need help, the input file should be like this:
10.00 154.00 0.01
11.00 151.90 0.01
12.00 152.00 0.01
13.00 154.30 0.01
14.00 158.80 0.01
15.00 165.50 0.01
16.00 174.40 0.01
17.00 185.50 0.01
18.00 198.80 0.01
19.00 214.30 0.01
20.00 232.00 0.01
21.00 251.90 0.01
22.00 274.00 0.01
23.00 298.30 0.01
24.00 324.80 0.01
25.00 353.50 0.01
26.00 384.40 0.01
27.00 417.50 0.01
28.00 452.80 0.01
29.00 490.30 0.01
30.00 530.00 0.01

       program BAY
	 implicit double precision (a-h, o-z)
	 logical jjpar, jjdat, jjth, jjg, jjrel, jjchl, jjreb,
     *      jjchb, jjgmg
c
	 common /number/ ndat, npar, iter, itmax, conver
	 common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     *	       varnew(325), idum(25)
c**** 235=(25*(25+1))/2
c
	 common /dat/ e(51), e2(51), data(51), vardat(1326), 
     *        en(1326), th(51), dum(51), sig(51)
c***** 1326=)51*(51+1))/2
c
	 common /both/ g(51,25), emg(51,25)
	 common /yrpar/ param(26), parcov(26,26), ydum(26), 
     *        iydum(26), ifpar(26), iwhere(26), nparam
c	 
	 data ndatmx /51/, nparmx /25/, nnparm/26/
	 data yes /1hy/
c
	 call outall(jjpar, jjdat, jjth, jjg, jjrel, jjchl, 
     *            jjreb, jjchb, jjgmg)
c
       write (5,99999)
	 read  (5,99998) itmax, conver
	 write (6,99997) itmax, conver
c	 
	 write(5,99996)
	 read(5,99995) auto 
	 if (auto.NE.yes) go to 10
c
	 write(5,99994)
	 read(5,99993) hbase
	 write(6,99992) hbase
c
   10  call pparam(nnparm)
c
	 if(npar.GT.nparmx) go to 40
c
	 if(jjpar) call outpar
c
	 if(jjpar .AND. npar.LT.nparam) call outypr
c
   20   call setdat
c
	  call fixv
c
	 if(ndat.GT.ndatmx) go to 50
c
	 if(jjdat) call outdat
	 iter=0

   30  call theory(auto, hbase)
	 if(jjth) call outth
	 if(jjg) call outg
c
	 call newpar(jjrel, jjchl, jjreb, jjchb, jjgmg)
c
	 if(iter.LT.itmax) call outp1
	 if(iter.EQ.itmax) call outp2
c
	 call update 
	 if(iter.EQ.itmax .AND. npar.LT.nparam) call outypr
c
	 iter=iter+1
	 if(iter.LE.itmax) go to 30
c
	 call restrt
	 go to 20
c	 
   40  write(5,99991) npar, nparmx
	 
c
   50  write(5,99990) ndat, ndatmx
	 

c ***********
99999  format('How many iterations?', 'cov fraction?', $) 
99998  format(i, f)
99997  format('Max number of iteration is', i3, 'convergence factor is')
99996  format('Do you wish to use automatic numerical derivative?', $)
99995  format(f)
99994  format('What is fraction difference for automatic derivative?',$)
99993  format(f) 
99992  format('automatic derivative uses step size')
99991  format('You wanna', i5, 'but you are allowed only', i5, 'parameter 
     * values.'/ 'change in commons' /'par'/ 'and',
     * /'both'/,' and in sunroutines pparam and oldp')
99990  format('You wannna', i5, 'but you are allowed only', i5, 
     * 'data values'/, 'change in commons' /'dat'/ 'and',
     * /'both',  'and in sunroutines setdat and fixv')
	 end 
c************************************
c************************************
	 subroutine pparam(nnparm)
	 implicit double precision (a-h,o-z)
	 common /number/ ndat, npar, iter, itmax, conver
	 common /par/ pold(25), parm(25), pdum(25), varpar(25), 
     &	        varnew(325), idum(25)
	 dimension parvar(25,25)
	 equivalence (parvar(1,1),varpar(1))
	 common /yrpar/ param(26), parcov(26,26), ydum(26), 
     &	        iydum(26), ifpar(26), iwhere(26), nparam
			
	 call setpar
	 if(nparam.GT.nnparm)  go to 150
	 n=0
	 do 10 i=1, nparam
	    if(ifpar(i).EQ.0) n=n+1
10     continue
	    if(n.LT.nparam) go to 30 
		if(n.EQ.0) go to 80
		do 20 i=1, nparam
		     ifpar(i)=1
20     continue
	    go to 80
30     do 60 k=1, nparam
	          kk=0
	          do 50 i=1, nparam
	                if(ifpar(i).NE.0) go to 50
	                do 40 j=1, nparam
	                if(j.EQ.i) go to 40
	                if(parcov(j,i).EQ.0.D0) go to 40
	                if(ifpar(j).EQ.0) go to 40
	                ifpar(i)= 2
	                kk = 1
	                go to 50
40     continue
50     continue
	        if(kk.EQ.0) go to 70
60     continue
70     continue

80     ipar = 0
	    jpar = 0 
	    do 90 i=1, nparam
	          if(ifpar(i).EQ.0) go to 90
	          ipar = ipar + 1 
	          parm(ipar)= param(i)
	          iwhere(ipar)=i
	          if(ifpar(i).EQ.2) jpar=jpar+1
90     continue
	    if(jpar.NE.0) write (6.99999) jpar
	    npar=ipar
	    ipar=0
	    do 110 i=1, nparam
	    if(ifpar(i).EQ.0) go to 110
	    ipar=ipar+1
	    jpar=0
	    do 100 j=1, nparam
	    if(ifpar(j).EQ.0) go to 100
	    jpar=jpar+1
	    parvar(jpar,ipar)=parcov(j,i)
100    continue
110    continue

	 do 120 i=1, npar
	     pold(i)=parm(i)
120    continue
	     il=0
	 do 140 i=1, npar
	       do 130 l=1, i
	       il=il+1
	       varpar(il) = parvar(l,i)
130    continue
140    continue
	 return 
150    write(5,99998) nparam
	     write(6,99998) nparam
	     stop
99999 format(46h0number of additional parameters which must be 
     &	           1h, 38hvaried because of covariance matrix is, i5)
99998 format(9h you want, i5, 21h but are allowed only, i5, 
     &              17hparameter values./28h change every array in comm 
     &				  1hn, 8h/yrpar/.)
	 end
c*****************************
c****************************

	 subroutine setpar
	 implicit double precision(a-h,o-z)
	 common /yrpar/ param(26), parcov(26,26), ydum(26), iydum(26),
     &	           ifpar(26), iwhere(26), nparam
	 nparam=3
	 write(5,99999)
	 read(5,99998) (param(i), i=1,nparam)

  	 do 40 i=1, nparam
	     do 30 l=1, i
	          parcov(l,i)= 0.D0
	          if(l.EQ.i) parcov(l,i)= .01D0*param(i)**2
30     continue
40     continue
	 return
99999  format(33h initial fuesses for parameters? $)
99998  format(3f)
	 end
c***************
c***************
        	subroutine setdat
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &	          en(1326), th(51), dum(51), sig(51)
	dimension datcov(51,51)
	equivalence (datcov(1,1), vardat(1))
	data dblank/10h           /
	write(5,99999)
	read(5,99998) name
	if(name.EQ.dblank) stop
	open(unit=24, file=name)
	n=0
 10    read(24,99997,end=20) ee, d, err
	if(ee.LE.0.D0) go to 20
	n= n+1
	go to 10
20    continue
	ndat=n
	rewind 24

	do 40 i=1, ndat
	    do 30 l=1, i
	      datcov(l,i)= 0.D0
30     continue
40     continue

	do 50 i=1, ndat
	    read(24,99997) ee, d, err
	    e(i)= ee
	    data(i)= d
	    datcov(i,i)= (err*d)**2
50     continue
	close(unit=24)
	return

99999 format(31h What's name of the data file? $)
99998 format(a10)
99997 format(3f2.0)
	end 
c*************
c*************
c*************
c*************
	double precision function theo(kdat)
	implicit double precision(a-h,o-z)
	common /dat/ e(51), e2(51), data(51), vardat(1326),
     &	          en(1326), th(51), dum(51), sig(51)
	common /yrpar/ param(26), parcov(26,26), ydum(26), iydum(26),
     &	           ifpar(26), iwhere(26), nparam
	theo= (param(1)**2*e(kdat)+param(2))*e(kdat) + param(3)
	return
	end
c***************
c**************
      double precision function deriv(kdat,kpar)
	implicit double precision(a-h,o-z)
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &	          en(1326), th(51), dum(51), sig(51)
	common /yrpar/ param(26), parcov(26,26), ydum(26), iydum(26), 
     &	           ifpar(26), iwhere(26), nparam
	go to (10, 20, 30), kpar
10     continue
	gg= 2.*param(1)*e(kdat)**2
	go to 40
20     continue
	gg= e(kdat)
	go to 40 
30     continue
	gg= 1.D0
40     continue
	deriv= gg
	return
	end
c***************
c***************

	subroutine fixv()
	 implicit double precision (a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &	        en(1326), th(51), dum(51), sig(51)
	dimension datcov(51,51)
	equivalence (datcov(1,1), vardat(1))
	il=0
	do 20 i=1, ndat
	   do 10 l=1, i
	      il=il+1
	      vardat(il)=datcov(l,i)
10     continue
20     continue
	 return
	 end
C*******************************
C*******************************
       	subroutine theory(auto, hbase)
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &	       en(1326), th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	common /deriva/ kdat, nderiv
	common /yrpar/ param(26), parcov(26,26), ydum(26), 
     &	       iydum(26), ifpar(26), iwhere(26), nparam
	dimension der(14), erest(14)
	external fun
	data yes /1hy/
	if (hbase.EQ.0.D0) hbase=0.001D0
	nder = 1
	ifail = 0
	do 50 idat=1, ndat
	    kdat=idat
	    th(idat)=theo(kdat)
	    if(auto.EQ.yes) go to 20
	    ipar=0
	    do 10 i=1, nparam
	         ii=i
	        if(ifpar(i).eq.0) go to 10
	        ipar=ipar+1
	        g(kdat,ipar)= deriv(kdat,ii)
10     continue
       go to 50
20     continue

	do 40 ipar=1,npar
	   nderiv=ipar
	   pipar= parm(ipar)
c	 call d04nml(pipar, nder, hbase, der, erest, fun, ifail)
c	  g(idat,ipar)=der(1)
c	  if(erest(1).LT.0.D0) go to 70 
c	   if(ifail.NE.0) go to 60
30       parm(ipar)= pipar
          param(iwhere(nderiv))=pipar
40     continue
50     continue
	 return 
c60     write(6,99999) ifail
c        go to 30
c70     write(6,99998) erest(1)
c	  go to 30
	
99999 format(30h0******error in d04nml, ifail=, i2)
99998 format(30h0******error in d04nml, erest=, 1pg14.6)
	end

	double precision function fun(x)
	implicit double precision(a-h,o-z)
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	common /yrpar/ param(26), parcov(26,26), ydum(26), 
     &	       iydum(26), ifpar(26), iwhere(26), nparam
	common /deriva/ kdat, nderiv
	 param(iwhere(nderiv))=x
	 parm(nderiv)=x
	 fun=theo(kdat)
	 return
	 end
	
c*******************************
c*******************************
       subroutine newpar(jjrel, jjchl, jjrep, jjchb, jjgmg)
	implicit double precision(a-h,o-z)
	logical jjrel, jjchl, jjreb, jjchb, jjgmg
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &             en(1326), th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	do 10 i=1, ndat
	    th(i)=data(i) - th(i)
10     continue
	if ((.NOT.jjrel) .AND. (.NOT.jjchl)) go to 80

	ij=0
	do 30 i=1, ndat
	  do 20 j=1, i
	       ij=ij+1
	       en(ij)=vardat(ij)
20     continue
30     continue 

	call scale(en, sig, ndat, ifdiag)
c	if(ifdiag.EQ.0) go to 40 
	iiiii=0

c	call sppco(en, ndat, rcond, dum, iiiii)
c	if(1.D0+rcond.EQ.1.D0) go to 310
c40     continue

	do 50 i=1, ndat
	    dum(i)= th(i)*sig(i)
50     continue
c	if(idiag.EQ.1) call dppsl(en, ndat, dum)
c	do 60 i=1, ndat
c	   dum(i)=dum(i)*sig(i) 
c60    continue
	if(jjrel) call outrel
	if(.NOT.jjchl) go to 80
	chi=0.
	do 70 i=1, ndat
	   chi= chi + th(i)*dum(i)
70     continue
	call outchl(chi)
80     continue

	if(iter.EQ.0) go to 110
	do 100 ipar=1, npar
	    do 90 i=1, ndat
	    th(i)= th(i) - g(i,ipar)*(pold(ipar)-parm(ipar))
90    continue
100   continue
110   continue
	call mul
	call mul2
	if(jjgmg) call outgmg(en, ndat)
	call scale(en, sig, ndat, ifdiag)
	if(jjgmg) call outgmg(en, ndat)
c	if(ifdiag.EQ.0) go to 120

	iiiii=0
c	call dppco(en, ndat, rcond, dum, iiiii)
c	if(1.0D0+rcond.EQ.1.0D0) go to 300
c120   continue
c	do 130 i=1, ndat
c	     dum(i)= th(i)*sig(i)
c130   continue 
c	if (ifdiag.EQ.1) call dppsl(en, ndat, dum)
c	do 140 i=1, ndat
c	    dum(i)= dum(i)*sig(i)
c140   continue 

	if(jjreb) call outreb(dum)
	if(.NOT.jjchb) go to 160
	chi=0
	do 150 i=1, ndat
	    chi= chi + th(i)*dum(i)
150   continue
	call outchb(chi)
160   continue

	do 170 i=1, npar
	    pdum(i)= 0.D0
170   continue
	do 190 i= 1, npar
	   do 180 j=1, ndat
	    pdum(i)= pdum(i) + emg(j,i)*dum(j)
180    continue
190    continue

	iconv=0 
	 do 210 i=1, npar
	     a= pdum(i) + pold(i)
	     if(conver.EQ.0.D0) go to 200
	     b= abs((a-parm(i))/parm(i))
	     if(b.LE.conver) iconv= iconv+1
200   parm(i)= a
210   continue
	if(iconv.EQ.npar) iter= itmax
	if(iter.NE.itmax) go to 290
	
	 ij=0
	 do 230 i=1, npar
	    do 220 j=1, i
		    ij = ij+1
			varnew(ij)=varpar(ij)
220    continue
230    continue
	  
       ij=0
	do 280 i=1, npar
	    do 240 j=1, ndat
	         dum(j)= emg(j,i)*sig(j)
240   continue

c	if(ifdiag.EQ.1) call dppsl(en, ndat, dum)

c	do 250 idat= 1, ndat
c	     dum(idat)= dum(idat)*sig(idat)
c250    continue

	do 270 j=1, i
	      ij= ij+1
	      do 260 k=1, ndat
	           varnew(ij)= varnew(ij) - emg(k,j)*dum(k)
260    continue 
270    continue
280    continue
290    return
300    write(5,99999) rcond
	     write(6,99999) rcond
	   stop
c	    go to 120
	 
310     write(5,99998) rcond
	     write(6,99998) rcond
c 	   stop
	     go to 230

99999 format(10x, 20h********************, /'recond=', 
     &            1pd12.6, /24h the matrix is possibly, 10hsingular-, 
     &	         30h- results may not be accurate., 10x,9h*********,
     &	              11h***********)
99998 format(10x, '********************', /  'recond=', 
     &	        1pd12.6, / 'vardat is possibly', 'singular- least- squares
     &            chi squared may not', 'be accurate.', 
     &				 /10x,'*********************')
	  end
	  
C************************
C************************
       subroutine scale(a, sig, n, ifdiag)
	implicit double precision (a-h,o-z) 
	dimension a(1), sig(n)
	  
	  il=0 
	do 20 i=1, n
	  il= il+i
	  si= a(il)
	  if(si.LE.0.D0) go to 10
	  sig(i)= 1.D0/sqrt(si)
	  go to 20
10    sig(i) = 1.D0/sqrt(-si)
20    continue

	  ifdiag=0
	  il=0
	  do 40 i=1, n
	      si= sig(i)
	       do 30 l= 1, i
	             il=il+1
	             if(i.NE.l .AND. a(il).NE.0.D0) ifdiag=1
	              a(il)= a(il)*si*sig(l)
30        continue
40      continue
	return
	end

c********************
c********************
       subroutine mul
	implicit double precision (a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	common /both/ g(51,25), emg(51,25)

	 jk=0
	do 60 j=1, npar
	   do 10 l= 1, ndat
	       emg(l,j)= 0.D0
10    continue 
	   do 50 k=1, npar
	       if(j.GE.k) go to 20
	       kj= (k*(k-1))/2 +j
	       em= varpar(kj)
	        go to 30
20     jk = jk+1
	        em= varpar(jk)
30     continue
	        if(em.EQ.0.D0) go to 50
	        do 40 l=1, ndat
	              emg(l,j) = emg(l,j) + em*g(l,k)
40    continue
50    continue 
60    continue
	return 
	end

c*********************
c*********************
       subroutine mul2
	implicit double precision (a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /dat/ e(51), e2(51), data(51), vardat(1326), 
     &	       en(1326), th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)

	 il=0
	 do 20 i=1, ndat
	    do 10 l=1, i
	         il= il+1
	         en(il)= vardat(il)
10    continue
20    continue
	 do 50 j=1, npar
	    il=0
	    do 40 i= 1, ndat
	        do 30 l= 1, i
	              il= il+1
	              en(il)= g(i,j)*emg(l,j) + en(il)
30      continue
40    continue 
50    continue
	 return
	end 

c****************
c****************
       subroutine update

	implicit double precision (a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	 common /yrpar/ param(26), parcov(26,26), ydum(26), 
     &	         iydum(26), ifpar(26), iwhere(26), nparam
	 
	 do 10 i=1, npar
	       param(iwhere(i)) = parm(i)
10    continue
	 if(iter.LT.ITMAX) return

	 ij=0
	 do 30 i=1, npar
	      iw= iwhere(i)
	      do 20 j=1, i
	      ij= ij+1
	      parcov(iwhere(j),iw) = varnew(ij)
20    continue
30    continue
	 return 
	end
c*********************
c*********************
        subroutine restrt
	
	implicit double precision (a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	       varnew(325), idum(25)
	 do 10 i= 1, npar
	       pold(i)= parm(i)
10     continue
       ij=0 
	do 30 i=1, npar
	     do 20 j=1, i
	         ij = ij+1
	         varpar(ij) = varnew(ij)
20    continue
30    continue
	 return
	 end

c*********************
c********************
       	subroutine outall(jjpar, jjdat, jjth, jjg, jjrel, jjchl, 
     &	                  jjreb, jjchb,  jjgmg)
	logical jjpar, jjdat, jjth, jjg, jjrel, jjchl, 
     &	        jjreb, jjchb, jjgmg, jjtrue, jjfals 
	dimension bb(14)

	data jjtrue/.true./, jjfals/.false./
	data blank/'  '/, out/'out'/
	data bb/'par', 'dat', 'res', 'chi', 'th', 'g', 'rel', 'chl', 'reb',
     &        'chb', 'gmg', 'th', 'g', 'g' /
	jjpar=jjtrue
	jjdat=jjtrue
	jjth=jjtrue
	jjg=jjtrue
	jjrel=jjtrue
	jjchl=jjtrue
	jjreb=jjtrue
	jjchb=jjtrue
	jjgmg=jjtrue
	open(unit=11, file='bayes.out')
10    continue
	read(11, 10000, end=40, err=40) a, b
10000 format(2a3)
	if(a.NE.out .AND. a.NE.blank) write(5,10100) a, b
10100 format(42h file bayes.out contains this unrecognized, 
     &	              12h line *****, 2a3, 6h *****)
	if(a.NE.out) go to 10
	do 20 i= 1, 14
	    if(b.EQ.bb(i)) go to 30
20    continue
	write(5,10100) a, b
	go to 10
30    continue
	if(i.EQ.1) jjpar=jjfals
	if(i.EQ.2) jjdat=jjfals
	if(i.EQ.3) jjrel=jjfals
	if(i.EQ.3) jjreb=jjfals
	if(i.EQ.4) jjchl=jjfals
	if(i.EQ.4) jjchb=jjfals
	if(i.EQ.5) jjth=jjfals
	if(i.EQ.6) jjg=jjfals
	if(i.EQ.7) jjrel=jjfals
	if(i.EQ.8) jjchl=jjfals
	if(i.EQ.9) jjreb=jjfals
	if(i.EQ.10) jjchb=jjfals
	if(i.EQ.11) jjgmg=jjfals
	if(i.EQ.12) jjth=jjfals
	if(i.EQ.13) jjg=jjfals
	if(i.EQ.14) jjg=jjfals
40    return
	end

c*****************
c****************
       subroutine outpar
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)

	write(6,99999)
	 ii=0
	do 10 i=1, npar
	    ii= ii+1
	    pdum(i)= dsqrt(varpar(ii))
10    continue
	write(6,99998)
	write(6,99997) (i,parm(i),pdum(i),i=1,npar)

	ii=0
	ioff=0
	iall=0
	do 30 i=1, npar
	      do 20 j=1, i
	        ii= ii+1
	        if(varpar(ii).NE.0.D0) iall= iall+1
	        if(i.EQ.j) go to 30
	        if(varpar(ii).NE.0.D0) ioff= ioff+1
20    continue
30    continue
	if(ioff.NE.0) call outv(varpar, pdum, idum, npar)
	if(iall.NE.0) return
	write(6,99996)
	stop
99999 format('******input parameters values')
99998 format('parameters       uncertainity')
99997 format(i5, 2f14.6)
99996 format( '******error', ' covariance matrix for parameters, 
     &	         be initialized in subroutine param. ******')
	end 

c*****************
c*****************
       subroutine outypr
	implicit double precision(a-h,o-z)
	common /yrpar/ param(26), parcov(26,26), ydum(26), iydum(26), 
     &	           ifpar(26), iwhere(26), nparam
	write(6,99999)
	write(6,99998)
	ipar=0
	do 20 i=1, nparam
	      a=dsqrt(parcov(i,i))
	      ydum(i)= a
	      if (ifpar(i).NE.0) go to 10
	      write(6,99997) i, param(i), a
	      go to 20
10    ipar= ipar+1
	      write(6,99996) i, param(i), a, ipar
20    continue 
	ioff=0
	iall=0 
	do 40 i=1, nparam
	     do 30 j=1,i
	        if(parcov(i,i).NE.0.D0) iall = iall+1
	        if(i.EQ.j) go to 40
	        if(parcov(j,i).NE.0.D0) ioff= ioff+1
30    continue
40    continue
	if(ioff.NE.0) call outyv
	if(iall.NE.0) return
	write(6,99995)
	stop
99999 format(40h0******complete list of parameters values)
99998 format('             parameters    uncertainity','   varied parameters
     & number')
99997 format(i5, 2f14.6)
99996 format(i5, 2f14.6, i10)
99995 format('******error. covariance martix for parameters', 
     &	     '  ', 'must be initialized in subroutine setpar.******')
	end

c*************
c*************
       subroutine outdat
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	dimension iddum(1)
	equivalence(iddum(1), sig(1))
	write(6,99999)
	ii=0
	do 10 i=1, ndat
	      ii= ii+1
	      dum(i)= dsqrt(vardat(ii))
10    continue
	write(6,99998)
	write(6,99997) (i,e(i),data(i),dum(i),i=1,ndat)
	ii=0 
	ioff=0
	do 30 i=1, ndat
	     do 20 j=1, i
	         ii= ii+1
	         if(i.EQ.j) go to 30
	         if(vardat(ii).NE.0.D0) ioff= ioff+1
20    continue
30    continue
	if(ioff.NE.0) call outv(vardat, dum, iddum, ndat)
	return

99999 format('   ******input data values')
99998 format('         data point        values      uncertainty')
99997 format(i5, 3f14.6)
	end

c****************
c***************
       subroutine outth
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	write(6,99999)
	write(6,99998)
	write(6,99997) (i,e(i),data(i),th(i),i=1,ndat)
	return

99999 format('   *****Theoretical calculation')
99998 format('            Energy        Data             theory')
99997 format(i5, 3f14.6)
	end

c*****************
c****************
       subroutine outg
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	data t2 /10h  Energy  /
	write(6,99999)
	min=1
	max= min0(npar, 7)
10    write(6,99998)
	 write(6,99997) t2, (i, i=min, max)
	  do 20 j=1, ndat
	     write(6,99996) j, e(j), (g(j,i), i=min,max)
20    continue
	if(max.EQ.npar) go to 30 
	  min = max+1
	  max= max+7
	if(max.GT.npar) max=npar
	 go to 10 
30    return 

99999 format('   ****** Partial Derivatives')
99998 format(20x)
99997 format(10x, a10, i10, 6i14)
99996 format(i4, 2x, 1pg14.6, 2x, 7g14.5)
	end

c************
c************
       subroutine outp1
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
	write(6,99999)
	write(6,99998)

	do 10 i=1, npar
	      write(6,99997) i, pold(i), parm(i)
10    continue
	return

99999 format('    ***** Intermediate results')
99998 format('           Old Param      New parameters')
99997 format(i5, 5f14.6)
	end

c************
c************
       subroutine outp2
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)

	write(6,99999)
	ii=0
	do 10 i=1, npar
	    ii= ii+1
	    pdum(i)= dsqrt(varnew(ii))
10    continue
	write(6,99998)
	write(6,99997) (i, pold(i), pdum(i), i=1, npar)
	ii=0
	ioff=0
	do 30 i=1, npar
	    do 20 j=1, i
	         ii= ii+1
	         if(i.EQ.j) go to 30
	         if(varnew(ii).NE.0.D0) ioff= ioff+1
20     continue
30     continue
	if(ioff.NE.0) call outv(varnew, pdum, idum, npar)
	return

99999 format('   ****** New values for parameters')
99998 format('            Pold         Pnew          uncertainty')
99997 format(i5, 3f14.6)
	end

c************
c************
       subroutine outv(v, s, ic, n)
	implicit double precision(a-h,o-z)
	dimension v(1), s(n), ic(n)

	write(6,99999)
	nn=0
	many=25
	m= n/many
	mm= m*many
	if(mm.NE.n) m= m+1
	imin= 1-many
	imax= min0(many,n)
	do 50 j=1, m
	   imin= imin+many
	   max= (imin*(imin-1))/2
	   write(6,99998) (i,i=imin,imax)
	   imax= imax+many
	   imax= min0(imax,n)
	   do 40 i=imin, n
	        min= max+imin
	        max= max+i
	        jmax= i
	        mm= jmax - imin + 1
	        if(mm.GT.many) jmax= imin + many - 1
	        si= s(i)
	        il= min - 1
	        do 10 l=imin, jmax
	             il= il+1
	             d= v(il)*100.D0/(s(l)*si)
	             if(d.GT.0.D0) d= d+0.5D0
	             if(d.LT.0.D0) d= d-0.5D0
	             ic(l)= d
10    continue
	             if(j.EQ.1) go to 30
	             not= 0
	             jj= jmax
	             if(jmax.EQ.i) jj= jj-1
	             jjp= jmax
	             if(jjp.Gt.jj) go to 30 
	             do 20 l=jjp, jj
	                 if(ic(l).EQ.0) go to 20
	                 not= 1
20    continue
	             if(not.EQ.1) go to 30 
	             nn=1
	             go to 40 
30          write(6,99997) i, si, (ic(l), l=imin, jmax)
40        continue
50     continue
	     return
99999 format('           std. Dev.         Correlation')
99998 format('  ', 17x, 25i4) 
99997 format(i4, 1pg14.6, 25i4)
	end

c************
c************
       subroutine outyv
	implicit double precision(a-h,o-z)
	common /yrpar/ param(26), parcov(26,26), ydum(26), iydum(26), 
     &	         ifpar(26), iwhere(26), nparam
	write(6,99999)
	n= nparam
	nn= 0
	many= 25
	m= n/many
	mm= m*many
	if(mm.NE.n) m= m+1
	imin= 1-many
	imax= min0(many,n)
	do 50 j=1, m
	   imin= imin+many
	   write(6,99998) (i, i=imin, imax)
	   imax= imax+many
	   imax= min0(imax,n)
	    do 40 i=imin, n
	        jmax=i
	        mm= jmax - imin + 1
	        if(mm.GT.many) jmax= imin + many - 1
	        si= ydum(i)
              do 10 l=imin, jmax
	             d= parcov(l,i)*100.D0/(ydum(l)*si)
	             if(d.GT.0.D0) d= d+0.5D0
	             if(d.LT.0.D0) d= d-0.5D0
	             iydum(l)=d
10       continue
	        if(j.EQ.1) go to 30
	        not=0
	        jj=jmax
	        if(jmax.EQ.I) jj= jj-1
	        jjp=imin
	        if(jjp.GT.jj) go to 30
	        do 20 l=jjp, jj
	           if(iydum(l).EQ.0) go to 20
	           not=1
20     continue
	        if(not.EQ.1) go to 30
	        nn=1
	        go to 40 
30     write(6,99997) i, si, (iydum(l), l=imin, jmax)
40    continue
50    continue
	    return 
99999 format('             std. dev.           correlation')
99998 format('  ', 17x, 25i4)
99997 format(i4, 1pg14.6, 25i4)
	     end
c*************
c*************
       subroutine outrel
	implicit double precision(a-h, o-z)
	common /number/ ndat, npar, iter, itmax, conver
	common /par/ pold(25), parm(25), pdum(25), varpar(325), 
     &	        varnew(325), idum(25)
	common /dat/ e(51), e2(51), data(51), vardat(1326), en(1326), 
     &	         th(51), dum(51), sig(51)
	common /both/ g(51,25), emg(51,25)
      dimension formt(4)
	data pl/1h(/, pr /1h)/
      data formt /'a1', 'a1', 'opf8.0', '1pg13.5' / 
	data  aa5/'      '/, aa4/'      '/, aa3/'       '/,
     & aa2/'      '/, aa1 /'     '/, aa0/'     '/, aa6/'      '/
	 data t1 / 10henergy     / , t2 / 10hresidual  /
	 write(6,99999)
	 go to 10
	 entry outreb
 	 write(6,99998)
10     write(6,99997) (t1,t2,i=1,4)
	 formt(3)= aa6
	 if(e(ndat).LT.1.D5) formt(3)=aa5
	 if(e(ndat).LT.1.D4) formt(3)=aa4
	 if(e(ndat).LT.1.D3) formt(3)=aa3
	 if(e(ndat).LT.1.D2) formt(3)=aa2
	 if(e(ndat).LT.1.D1) formt(3)=aa1
	 if(e(ndat).LT.1.D0) formt(3)=aa0
	 n=ndat
	 m= n/4
	 mm= m*4
	 if(mm.NE.N) go to 20
	 mm=m
	 go to 30
20     m= m+1
	    k= mm + 4 - n
	    if(k.GE.m) go to 70
	    mm= m - k
30     continue
	   im= 3*m
	   do 40 i=1, mm
	      im = im+1
	      write(6,formt) (pl, l, pr, e(l), dum(l), l=i, im, m)
40     continue
	    if(mm.EQ.m) go to 60
	    im= 2*m + mm
	    mm= mm + 1
	    do 50 i=mm, m
	         im = im+1
	         write(6,formt) (pl,l,pr,e(l),dum(l),l=i,im,m)
50     continue
60     continue
	  return
70     write(6,formt) (pl,l,pr,e(l),dum(l),l=i,im,m)
	 return
99999 format('   ***** Least-square weighted residual, at former 
     & values of parameters')
99998 format('   ***** bayesian weighted residual,  at former values
     & of parameters')
99997 format('  ' , 4(2x, a10, 2x, a10, 4x))
	end
c***************
c***************
       	subroutine outchl(chi)
	implicit double precision(a-h,o-z)
	common /number/ ndat, npar, iter, itmax, conver
	write(6,99999) chi
	n= ndat - npar
	if(n.LE.0) return
	chi = chi/dfloat(n)
	write(6,99998) chi
	return
	
	entry outchb(chi)
	write(6,99997) chi
	n= ndat - npar
	if(n.LE.0) return
	chi = chi/dfloat(n)
	write(6,99998) chi
	return 
99999 format('   ***** Least aquare chi squared at former values of
     & parameters is', f15.6)
99998 format('chi squared divided by degree of freedom is', f15.6)
99997 format('   ***** bayesian chi squared at former values of 
     &  parameters is', f15.6)
	end
c************
c************
       subroutine outgmg(en, ndat)
	implicit double precision(a-h,o-z)
	dimension en(1)
	many=8
	m=ndat/many
	if(mm.NE.ndat) m= m+1
	imin= 1 - many
	do 20 j=1,m
	    imin= imin+many
	    max= (imin*(imin-1))/2
	    do 10 i=imin, ndat
	        min= max+imin
	        max= max+i
	        mmax=max
	        mm= max - min + 1
	        if(mm.GT.many) mmax= min + many - 1 
	        write(28,99999) i, (en(il), il=min,mmax)
10    continue
	       write(28,99998)
20    continue
	return
99999 format(i4, 8g14.6)
99998 format(1h0)
	end

 

[mikrom]

Cannot compile your program. What compiler are you using?

 

[xwb]

The syntax looks like a mix of F66 & F77, but more toward F66 (Fortran IV). The open statements look at bit like F77 but it could be a special F66 variant. Don't know which compiler takes things like format(i,f). If you set the tabs to 5 spaces, most of it will compile.

 

[gullipe]

Hi Tony1984

I tryed to compile the program with a F77 compiler and got the following errors (and also found most of them):

bay.for(80) : warning F4980: integer expected in format ("i,f" instead of "i10,f12.3" for instance)
bay.for(83) : warning F4980: integer expected in format (only "f" instead of "f12.3" for instance)
bay.for(85) : warning F4980: integer expected in format (only "f" instead of "f12.3" for instance)
bay.for(106) : warning F4329: PAR : COMMON : size changed (could not see that in a quick look)
bay.for(144) : error F2707: illegal unit specifier (6.formatnumber instead of 6,formatnumber)
bay.for(197) : warning F4980: integer expected in format (only "3f" instead if "3f12.3" or something like that)
bay.for(212) : error F2703: FILE= : not CHARACTER (variable NAME not specified as Character*32 for instance)
bay.for(372) : error F3606: NEWPAR : formal argument JJREP : type mismatch  (JJPEP should be JJREB)
bay.for(500) : warning F4801: label 230 : used across blocks

Try to correct this and compile again.

 

[Tony1984]

i am using compaq visual fortran

 

[mikrom]

I don't have compaq visual fortran. If you would be able to standardize the source, so that it could be compiled in e.g. gfortran, maybe somebody could help you to find the problem.

 

[gullipe]

Hi again
I have Compacq Visual Fortran, and when I compiled your program i got the following "warnings" only !

But it is at least clear that the variable JJREP that appears in the following statement (in line 372), should be JJREB:
subroutine newpar(jjrel, jjchl, jjrep, jjchb, jjgmg)

--------------------Configuration: bay - Win32 Debug--------------------
Compiling Fortran...
C:\GP\AAA\bay.for
C:\GP\AAA\bay.for(500) : Warning: A jump into a block from outside the block may have occurred.   [230]
      go to 230
-------------------^
C:\GP\AAA\bay.for(435) : Warning: In the call to OUTREB, actual argument #1 has no corresponding dummy argument.
 if(jjreb) call outreb(dum)
------------------------------^
C:\GP\AAA\bay.for(435) : Warning: Variable JJREB is used before its value has been defined
 if(jjreb) call outreb(dum)
-----------^
C:\GP\AAA\bay.for(435) : Warning: Routine OUTREB called with different number and/or type of actual arguments in earlier call - C attribute required if intended.
 if(jjreb) call outreb(dum)
------------------^

bay.obj - 0 error(s), 4 warning(s)

 

[xwb]

If you run your program in CVF, it will stop in the debugger where it crashed. You can then have a look at the variables, where it has come from etc.