      PROGRAM make_cat
c     Version: 11/12/2000
      IMPLICIT NONE
      INTEGER maxw,mxft,mxwf,mxfil,maxt,mxtyp,input,nparam,iseed,
     .     i,ii,iii,j,l,k,kk,ityp,kfil,icomp,nfil,iblank,inoise,ks,
     .     iread,iw,i_ob,jread,last,law_rd,nob,mxobj,nnn,nf,np,
     .     ncomp,n_filter,n_files,nfil_ref,i1,i2,i_age
      REAL yz,yob,amag,av,avmax,avmin,averr,bmag,bp1,dif,
     .     z_int,err_p,f,rangauss,ran1,f_red,h0,om,ov,r,sig,surf,
     .     tz,ta,w,wl,xa0v,ya0v,flux,z,vega,age,age_z,age_form,
     .     z_max,z_min,z_form,filter,t1,t2,tx,a1,a2,diff_mag,ww,
     .     am_min,am_max,amag_ref,zero,am1,am2,am_err,sn0,am0,sn,
     .     f_sn,sn_0,am_0,amlimit,err_min
      PARAMETER (nparam=23,mxft=300,mxwf=200,maxw=1300,maxt=225,
     .     mxfil=30,mxtyp=20,mxobj=10000,am1=10.,am2=40.,am_err=1.e-3)
c       mxft = 300 = maximum number of filters
c       mxwf = 200 = maximum number of wavelengths/filter
      CHARACTER NAMECH(mxtyp)*160,sed_abc(mxtyp)*2,line*240
      CHARACTER*160 spectra1,spectra2,filedef,filein,sedfile,
     .     filtresfile,filt_file,fmt,fmt2,line_sed,copyright1,
     .     copyright2,copyright3
      CHARACTER*40 cname,param(nparam)
      DIMENSION input(nparam),yz(maxw),yob(maxw)
      DIMENSION nfil(mxfil),amag(mxfil),amlimit(mxfil)
      DIMENSION err_p(mxfil),kfil(10),sn0(mxfil),am0(mxfil)
      DIMENSION f_red(maxw,maxt,mxtyp),KS(mxtyp)
      DIMENSION xa0v(5000),ya0v(5000),vega(0:mxfil)
      DIMENSION wl(0:mxfil),surf(0:mxfil)
      DIMENSION ta(maxt),ww(maxw),f(maxw,maxt),w(maxw,mxtyp),iw(mxtyp)

      COMMON /c_filter/ iread,jread,nf,np(mxft),r(mxwf,mxft,2)
      COMMON /compose/ icomp
      COMMON /namefile/ sedfile
      COMMON /filterres/ filtresfile
      COMMON /s2n/sn_0,am_0
      EXTERNAL ran1,rangauss,filter,tz,zero,f_sn
c_________________________________________________________________________

      copyright1=
     .     '  make_catalog - synthetic catalogue code - Version 1.1'
     .     //' (December 11, 2000)'
      copyright2=' by Micol Bolzonella, Roser Pello'','
     .     //' Joan-Marc Miralles'
      copyright3=' More information at '
     .     //'http://webast.ast.obs-mip.fr/hyperz/'
      WRITE(*,*)
      WRITE(*,'(a80)') copyright1
      WRITE(*,*)
      WRITE(*,'(a80)') copyright2
      WRITE(*,'(a80)') copyright3
      WRITE(*,*)

      jread=0                   !Control on reading A0VSED
      iread=0                   !Control on reading Filter Library

      param(1)='A0VSED'
      param(2)='FILTERS_RES'
      param(3)='FILTERS_FILE'
      param(4)='TEMPLATES_FILE'
      param(5)='MODEL_CAT_FILE'
      param(6)='NOISE_CAT_FILE'
      param(7)='REF_FILTER'
      param(8)='MAG_MIN'
      param(9)='MAG_MAX'
      param(10)='N_OBJ'
      param(11)='Z_FORM'
      param(12)='Z_MIN'
      param(13)='Z_MAX'
      param(14)='AGE_RANDOM'
      param(15)='REDDENING_LAW'
      param(16)='AV_MIN'
      param(17)='AV_MAX'
      param(18)='NOISE_TYPE'
      param(19)='ERR_MAG_MIN'
      param(20)='H0'
      param(21)='OMEGA_M'
      param(22)='OMEGA_V'
      param(23)='RAND_SEED'
      DO i=1,nparam
         input(i)=0
      END DO

 44   filein='catalog.param'
      last=INDEX(filein,'  ')-1
      WRITE(*,*) 'Name of parameter file' 
      WRITE(*,*) '[default = ',filein(1:last),']'
      READ (*,'(A)') line
      IF (line.ne.' ') THEN
         READ(line,'(A)') filein
      END IF

      OPEN(77,FILE=filein,STATUS='old',ERR=44)
 7    READ(77,'(A)',ERR=10,END=10) line
      DO i=1,80
         IF (line(i:i).ne.' ') THEN
            line=line(i:)
            GO TO 8
         END IF
      END DO
      GO TO 7  ! blank line
 8    IF (line(1:1).eq.'#') GO TO 7 ! comment line
      iblank=INDEX(line,' ')
      cname=line(1:iblank-1) ! variable name
      line=line(iblank:)
      DO i=1,80
         IF (line(i:i).ne.' ') THEN
            line=line(i:)
            IF (line(1:1).eq.'#') THEN
               WRITE(*,*) 'Missed value of variable ', cname
               STOP
            END IF
            GO TO 9
         END IF
      END DO
      WRITE(*,*) 'Missed value of variable ', cname
      STOP
 9    iblank=INDEX(line,' ')
      line=line(1:iblank-1) ! part of string containing value of variable
c ..... A0VSED file
      IF (cname.eq.param(1)) THEN
         READ(line,'(A)') sedfile
         input(1)=1
         GO TO 7
c .... FILTER file
      ELSE IF (cname.eq.param(2)) THEN
         READ(line,'(A)') filtresfile
         input(2)=1
         GO TO 7
c ..... File containing filters and noise sigma
      ELSE IF (cname.eq.param(3)) THEN
         READ(line,'(A)') filt_file
         input(3)=1
         GO TO 7
c ..... File containing synthetic spectra filenames
      ELSE IF (cname.eq.param(4)) THEN
         READ(line,'(A)') filedef
         input(4)=1
         GO TO 7
c ..... output file catalogues
      ELSE IF (cname.eq.param(5)) THEN
         READ(line,'(A)') spectra1
         input(5)=1
         GO TO 7
      ELSE IF (cname.eq.param(6)) THEN
         READ(line,'(A)') spectra2
         input(6)=1
         GO TO 7
c ..... Reference filter 
      ELSE IF (cname.eq.param(7)) THEN
         READ(line,*) nfil_ref
         input(7)=1
         GO TO 7
c ..... Minimum magnitude in reference filter
      ELSE IF (cname.eq.param(8)) THEN
         READ(line,*) am_min
         input(8)=1
         GO TO 7
c ..... Maximum magnitude in reference filter
      ELSE IF (cname.eq.param(9)) THEN
         READ(line,*) am_max
         input(9)=1
         GO TO 7
c ..... number of objects
      ELSE IF (cname.eq.param(10)) THEN
         READ(line,*) nob
         IF (nob.gt.mxobj) THEN
            WRITE(*,'(1x,a,i7,a)') 
     .           'Number of objects exceeding the maximum [max=',
     .           mxobj,']'
            STOP
         END IF
         INPUT(10)=1
         GO TO 7
c ..... z_form: formation redshift
      ELSE IF (cname.eq.param(11)) THEN
         READ(line,*) z_form
         INPUT(11)=1
         GO TO 7
c ..... Redshift range
      ELSE IF (cname.eq.param(12)) THEN
         READ(line,*) z_min
         INPUT(12)=1
         GO TO 7
      ELSE IF (cname.eq.param(13)) THEN
         READ(line,*) z_max
         INPUT(13)=1
         GO TO 7
c ..... Age: choice for galaxy age
c       0. age = age of universe T(z_form)
c       1. age = random between 0 and T(z)-T(z_form)
      ELSE IF (cname.eq.param(14)) THEN
         READ(line,*) i_age
         IF (i_age.ne.0.and.i_age.ne.1) THEN
            WRITE(*,*) 'AGE_RANDOM must be 0 or 1'
            STOP
         END IF
         INPUT(14)=1
         GO TO 7
c ..... Reddening:
c       0 -> No reddening
c       1 -> Allen, 76, standard MW
c       2 -> Seaton, 79, standard MW
c       3 -> Fitzpatrick, 86, LMC
c       4 -> Prevot, 84, + Bouchet, 85, SMC 
c       5 -> Calzetti, astro-ph/9911459
      ELSE IF (cname.eq.param(15)) THEN
         READ(line,*) law_rd
         IF (law_rd.ne.0.and.law_rd.ne.1.and.law_rd.ne.2.and.
     .        law_rd.ne.3.and.law_rd.ne.4.and.law_rd.ne.5) THEN
            WRITE(*,*) 'REDDENING_LAW must be 0 -- 5'
            STOP
         END IF
         input(15)=1
         GO TO 7
c ..... Av_max, Av_min
      ELSE IF (cname.eq.param(16)) THEN
         READ(line,*) avmin
         input(16)=1
         GO TO 7
      ELSE IF (cname.eq.param(17)) THEN
         READ(line,*) avmax
         input(17)=1
         GO TO 7
c ..... type of noise:
c       0. random with fixed sigma in magnitudes
c       1. with dm(m)
      ELSE IF (cname.eq.param(18)) THEN
         READ(line,*) inoise
         IF (inoise.ne.0.and.inoise.ne.1) THEN
            WRITE(*,*) 'NOISE_TYPE must be 0 or 1'
            STOP
         END IF
         INPUT(18)=1
         GO TO 7
      ELSE IF (cname.eq.param(19)) THEN
         READ(line,*) err_min
         INPUT(19)=1
         GO TO 7
c ..... Cosmological parameters: H_0, Omega_matter, Omega_lambda
      ELSE IF (cname.eq.param(20)) THEN
         READ(line,*) h0
         INPUT(20)=1
         GO TO 7
      ELSE IF (cname.eq.param(21)) THEN
         READ(line,*) om
         INPUT(21)=1
         GO TO 7
      ELSE IF (cname.eq.param(22)) THEN
         READ(line,*) ov
         INPUT(22)=1
         GO TO 7
c ..... Seed for random generator
      ELSE IF (cname.eq.param(23)) THEN
         READ(line,*) iseed
         input(23)=1
         GO TO 7
      ELSE 
         WRITE(*,'(a,1x,a)') 'Unknown parameter (ignored)',cname
         GO TO 7
      END IF
 10   CONTINUE
      CLOSE(77)
      DO i=1,nparam
         IF (input(i).eq.0) THEN 
            WRITE(*,*) 'Missed parameter ', param(i)
            STOP
         END IF
      END DO

      IF (law_rd.eq.0) THEN
         avmin=0.0
         avmax=0.0
      END IF
      IF (z_max.gt.z_form) z_max=z_form

c ..... reading the names of SEDs

      OPEN(10,FILE=filedef,STATUS='old',ERR=78)
      iii=1
 2230 READ(10,'(A)',ERR=2220,END=2220) line_sed
      IF (line_sed(1:1).eq.'#'.or.line_sed.eq.' ') GO TO 2230
      last=INDEX(line_sed,' ')-1
      NAMECH(iii)=line_sed(1:last)
      line_sed=line_sed(last+1:)
      DO i=1,80
         IF (line_sed(i:i).ne.' ') THEN
            line_sed=line_sed(i:)
            last=i
            GO TO 2221
         END IF
      END DO
      WRITE(*,*) 'Missed type of SED file ', NAMECH(iii)
      STOP
 2221 sed_abc(iii)=line_sed(1:2)
      iii=iii+1
      GO TO 2230
 2220 n_files=iii-1
      IF (n_files.gt.mxtyp) THEN
         WRITE(*,'(1x,a,i3,a)') 
     .        'Number of files exceeding the maximum [max=',mxtyp,']'
         STOP
      END IF
      WRITE(*,*) 'There are ',n_files,' synthetic spectra'
      CLOSE(10)

c ..... Reading the SEDs

      WRITE(*,*) 'Reading synthetic spectra'
      DO i=1,n_files
         WRITE(*,2005) i, NAMECH(i)
c GISSEL's templates
         IF (sed_abc(i).eq.'BC') THEN
            OPEN(1,FILE=NAMECH(i),STATUS='old',
     .           FORM='unformatted',ERR=25)
            READ (1) KS(i),(Ta(ii),ii=1,KS(i)) ! KS = n. temporal intervals
            READ (1) IW(i),(W(ii,i),ii=1,IW(i)) !Number of wavelengths
            DO j=1,KS(i)
               READ (1) IW(i),(F_red(ii,j,i),ii=1,IW(i))
            ENDDO
            CLOSE(1)
c extended CWW's templates
         ELSE IF (sed_abc(i).eq.'AS'.or.sed_abc(i).eq.'A ')THEN
            OPEN(1,FILE=NAMECH(i),STATUS='old',
     .           FORM='formatted',ERR=25)
            KS(i)=1
            ii=1
 33         READ(1,*,ERR=34,END=34) W(ii,i),F_red(ii,1,i)
            ii=ii+1
            GO TO 33
 34         iw(i)=ii-1
            CLOSE(1)
         ELSE
            WRITE(*,*) 'Don''t know how to read :'
            WRITE(*,'(A)') NAMECH(i)
            STOP
         END IF
c         WRITE(*,102) IW(i),W(1,i),W(IW(i),i),KS(i)
      END DO
 2005 FORMAT(' Type ', I2, ' : ', A60)
c 102  FORMAT(i5,' intervals from lambda = ',F5.0, ' to ',
c     .     e9.3,' A. and ',I4,' temporal intervals')

c ..... Reading filt_file (format depending on noise type)

      CALL filter0
      OPEN(10,FILE=filt_file,STATUS='old',ERR=80)
      i=0
 134  CONTINUE
      READ(10,'(A)',ERR=135,END=135) line
      IF (line(1:1).eq.'#'.or.line.eq.' ') GO TO 134
      i=i+1
      ncomp=1
      DO k=1,80
         IF (line(k:k).eq.'&') THEN
            line(k:k)=' '
            ncomp=ncomp+1
         END IF
      END DO
c ..... simple filter
      IF (ncomp.eq.1) THEN
         IF (inoise.eq.0) THEN
            READ(line,*,ERR=82) nfil(i), err_p(i)
            IF (nfil(i).gt.icomp) THEN
               WRITE(*,1998) nfil(i),icomp
               STOP
            ELSE
               kfil(1)=nfil(i)
               CALL compos_new(1,kfil,wl(i),surf(i),bp1)
               GO TO 134
            END IF
         ELSE IF (inoise.eq.1) THEN
            READ(line,*,ERR=81) nfil(i), sn0(i), am0(i)
            IF (nfil(i).gt.icomp) THEN
               WRITE(*,1998) nfil(i),icomp
               STOP
            ELSE
               kfil(1)=nfil(i)
               CALL compos_new(1,kfil,wl(i),surf(i),bp1)
               GO TO 134
            END IF
         END IF 
c ..... composite filter
      ELSE 
         IF (inoise.eq.0) THEN
            READ(line,*,ERR=82) (kfil(j),j=1,ncomp), err_p(i)
            icomp=icomp+1
            nfil(i)=icomp
            CALL compos_new(ncomp,kfil,wl(i),surf(i),bp1)
            GO TO 134
         ELSE IF (inoise.eq.1) THEN
            READ(line,*,ERR=81) (kfil(j),j=1,ncomp), sn0(i), am0(i)
            icomp=icomp+1
            nfil(i)=icomp
            CALL compos_new(ncomp,kfil,wl(i),surf(i),bp1)
            GO TO 134
         END IF 
      END IF
 1998 FORMAT('Filter number ',I4, ' out of range!',/,
     .       'Maximum = ',I4)
 135  n_filter=i
      IF (n_filter.gt.mxfil) THEN
         WRITE(*,'(1x,a,i3,a)') 
     .        'Number of filters exceeding the maximum [max=',mxfil,']'
         STOP
      END IF
      WRITE(*,*) 'There are ',n_filter,' filters'
      CLOSE(10)
c ..... Reference filter
      kfil(1)=nfil_ref
      CALL compos_new(1,kfil,wl(0),surf(0),bp1)

      IF (inoise.eq.1) THEN
c ..... Computing limiting magnitudes
         WRITE(*,*) 'Limiting magnitudes:'
         WRITE(*,*) '# filter     m_limit'
         DO kk=1,n_filter
            sn_0=sn0(kk)
            am_0=am0(kk)
            amlimit(kk)=zero(f_sn,am1,am2,am_err)
            WRITE(*,'(4x,i3,7x,f6.3)') nfil(kk),amlimit(kk)
         END DO
      END IF

      nnn=n_filter+2
      WRITE(fmt,1002) nnn
 1002 FORMAT('(1x,i4,1x,f7.3,1x,i2,',i2,'(1x,f7.3))')

      WRITE(fmt2,2001) n_filter,n_filter
 2001 FORMAT('(i4,1x,',i2,'(f7.3,1x),',i2,'(f6.3,1x),2(i2,1x))')

c ..... Reading Vega spectrum

      WRITE (6,'(1x,a)') 'Reading reference sed:'
      CLOSE(81)
      OPEN(81,file=sedfile,form='formatted',status='old')
      READ(81,'(a)',err=3) line
      WRITE(6,'(1x,a)') line (1:78)
      ii=0
 138  CONTINUE
      ii=ii+1
      READ(81,*,err=139,end=139) xa0v(ii),ya0v(ii)
      jread=ii
      GO TO 138
 139  CLOSE(81)

      DO kk=1,n_filter
         vega(kk)=filter(nfil(kk),xa0v,ya0v,jread,0.)
      END DO
      vega(0)=filter(nfil_ref,xa0v,ya0v,jread,0.)

c ----- End of preliminaries

      age_form=tz(h0,om,ov,z_form) ! age of universe at z_form
      z_int=z_max-z_min 
      averr=avmax-avmin
      OPEN(14,file=spectra1,status='unknown')
      OPEN(15,file=spectra2,status='unknown')

c ------------------------------------------------------------------------
c homogeneous distribution in z and type, mag random in the input range
c ------------------------------------------------------------------------
      DO i_ob=1,nob
         z=z_int*ran1(iseed)+z_min ! random z 
         bmag=am_min+ran1(iseed)*(am_max-am_min) ! random mag in ref.filter
         av=avmin+averr*ran1(iseed) ! random Av
         ityp=INT(n_files*ran1(iseed))+1 ! random type
         IF (sed_abc(ityp).eq.'BC') THEN
            age_z=tz(h0,om,ov,z)    ! age of universe at z [Gyr]
            IF (i_age.eq.0) THEN
               age=age_z-age_form
            ELSE IF (i_age.eq.1) THEN
               age=(age_z-age_form)*ran1(iseed) ! random age [Gyr]
            END IF
            DO k=1,ks(ityp)-1
               IF (age.ge.ta(k)*1.e-9.and.age.lt.ta(k+1)*1.e-9) THEN
                  i1=k
                  i2=k+1
                  GO TO 22
               END IF
            END DO
 22         DO l=i1,i2
               DO i=1,iw(ityp)  !lambda
                  f(i,l)=f_red(i,l,ityp)
                  ww(i)=w(i,ityp)
                  IF (ww(i).lt.912.0) f(i,l)=0.0
               END DO
            END DO
            tx=age
            t1=ta(i1)*1.e-9
            t2=ta(i2)*1.e-9
            a1=(t2-tx)/(t2-t1)
            a2=1.-a1
            DO i=1,iw(ityp)
               yz(i)=a1*f(i,i1)+a2*f(i,i2)
            END DO
         ELSE IF (sed_abc(ityp).eq.'AS'.or.sed_abc(ityp).eq.'A ') THEN
            age=19.999
            DO i=1,iw(ityp)
               yz(i)=f_red(i,1,ityp)
               ww(i)=w(i,ityp)
               IF (ww(i).lt.912.0) yz(i)=0.0
            END DO
         END IF
         IF (av.gt.0.0) THEN
            CALL red_spec(ww,yz,av,yob,1,iw(ityp),law_rd)
            DO k=1,iw(ityp)
               yz(k)=yob(k)
            END DO
         END IF
c ... lyman alpha forest
         CALL lyman_forest(ww,yz,iw(ityp),z)
         flux=filter(nfil_ref,ww,yz,iw(ityp),z)
         IF (flux.eq.0.) THEN
            WRITE(*,'(a,1x,f6.3,1x,a)') 'Object at redshift ',z, 
     .           'non detected in the reference filter'
            WRITE(*,*) 'Select a redder REF_FILTER or a lower Z_MAX'
            STOP
         END IF
         amag_ref=-2.5*LOG10(flux/vega(0))
         diff_mag=amag_ref-bmag
         DO kk=1,n_filter
            flux=filter(nfil(kk),ww,yz,iw(ityp),z)
            IF (flux.gt.0.) THEN
               amag(kk)=-2.5*LOG10(flux/vega(kk))-diff_mag
            ELSE
               amag(kk)=99.
            END IF
         END DO
c model catalog (no noise added)
         WRITE(14,fmt) i_ob,z,ityp,age,(amag(kk),kk=1,n_filter),av
c errors
         IF (inoise.eq.1) THEN
            DO kk=1,n_filter
               err_p(kk)=1.
               IF (amag(kk).ne.99.) THEN 
                  sn=sn0(kk)*10.**(-0.4*(amag(kk)-am0(kk)))
                  IF (sn.lt.1.) THEN 
                     amag(kk)=99.
                  ELSE
                     sig=SQRT((2.5*LOG10(1.+1./sn))**2+err_min**2)
                     dif=rangauss(iseed)
                     amag(kk)=amag(kk)+dif*sig
                     err_p(kk)=sig
                  END IF
               END IF
            END DO
         ELSE IF (inoise.eq.0) THEN
            DO kk=1,n_filter
               IF (amag(kk).ne.99.) THEN 
                  sig=err_p(kk)
                  dif=rangauss(iseed)
                  amag(kk)=amag(kk)+dif*sig
               END IF
            END DO
         END IF
c catalog with errors
         WRITE(15,fmt2) i_ob,(amag(kk),kk=1,n_filter),
     .        (err_p(kk),kk=1,n_filter)
      END DO
      CLOSE(14)
      CLOSE(15)

      STOP
 25   STOP 'Error opening synthetic spectra files'
 78   STOP 'Error opening templates file'
 80   STOP 'Error opening filters file'
 81   STOP 'Error reading filters file (format: #filter s/n mag)'
 82   STOP 'Error reading filters file (format: #filter sigma)'
 3    STOP 'Program exits because of error reading A0 V s.e.d. file'
      END
c............................................................
      FUNCTION f_sn(am)
      REAL f_sn,am,sn_0,am_0
      COMMON/s2n/sn_0,am_0
      f_sn=sn_0*10.**(-0.4*(am-am_0))-1.
      RETURN
      END
c............................................................

      FUNCTION ran1(idum)
      DIMENSION rr(97)
      PARAMETER (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
      PARAMETER (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
      PARAMETER (m3=243000,ia3=4561,ic3=51349)
      SAVE rr,ix1,ix2,ix3
      DATA iff /0/
      IF (idum.lt.0.or.iff.eq.0) THEN
         iff=1
         ix1=MOD(ic1-idum,m1)
         ix1=MOD(ia1*ix1+ic1,m1)
         ix2=MOD(ix1,m2)
         ix1=MOD(ia1*ix1+ic1,m1)
         ix3=MOD(ix1,m3)
         DO j=1,97
            ix1=MOD(ia1*ix1+ic1,m1)
            ix2=MOD(ia2*ix2+ic2,m2)
            rr(j)=(FLOAT(ix1)+FLOAT(ix2)*rm2)*rm1
         END DO
         idum=1
      END IF
      ix1=MOD(ia1*ix1+ic1,m1)
      ix2=MOD(ia2*ix2+ic2,m2)
      ix3=MOD(ia3*ix3+ic3,m3)
      j=1+(97*ix3)/m3
      IF (j.gt.97.or.j.lt.1) PAUSE
      ran1=rr(j)
      rr(j)=(FLOAT(ix1)+FLOAT(ix2)*rm2)*rm1
      RETURN
      END

c............................................................
c  Normal (Gaussian) deviates with unit variance
c  Returns normally distributed deviate with zero mean
c............................................................

      FUNCTION rangauss(idum)
      REAL a,b,rg,rangauss
      INTEGER idum
 50   a=ran1(idum)
      b=ran1(idum)
      b=1.7156*(b-0.5)
      rg=b/a
      IF (rg**2.gt.-4.*LOG(a)) GO TO 50
      rangauss=rg
      RETURN
      END
c............................................................
c Zero of function func with bisection

      FUNCTION zero(func,xl,xu,err)
      REAL zero, func, xl, xu, x1, x2, xnew, err, test, errp
      INTEGER itmax
      PARAMETER (itmax=5000)
      EXTERNAL func
      IF (func(xl)*func(xu).gt.0.) THEN
         STOP ' ahi! zero non bracketed ...'
      END IF
      x1=xl
      x2=xu
      DO i=1,itmax
         xnew=(x2+x1)/2.
         IF (x1+x2.ne.0.) THEN
            errp=ABS((x2-x1)/(x2+x1))*100.
         END IF
         test=func(x1)*func(xnew)
         IF (test.eq.0.) THEN
            errp=0.
         ELSE
            IF (test.lt.0.) THEN
               x2=xnew
            ELSE
               x1=xnew
            END IF
         END IF
         IF (errp.le.err) GO TO 20
      END DO
 20   zero=xnew
      RETURN
      END
