subroutine mrqfit(modvis2,modphas,adjust,mcrun)
c
c     Fits binary model to data by calling MRQMIN
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
c
      real*8 modvis2(mxrec,mxfilt),modphas(mxrec,mxfilt)
c
c     lista=array with first mfit entries equal to the numbers of the
c     parameters to be varied.
      integer*4 lista(mxparm),mfit
      real*8 chisq,chisqo,chiini,chifr,tol,alamda,seprad,parad
      character data*80,parid2*2,parid3*3
      logical adjust,convergd,mcrun,mcfirst

c Arrays used by MRQMIN
      real*8 xdata(mxrec2),ydata(mxrec2),sigma(mxrec2)
      real*8 a(mxparm)
      real*8 covar(mxparm,mxparm),alpha(mxparm,mxparm)
      external funcs

      data mcfirst/.true./
c
      write(outc,*)
      write(outc,*)'MRQFIT begins...'
c --- Set default parameters
      modef=0 ! Gaussian model in FCHISQ
      tol=1.d-10
      if(mcrun)then
       chifr=0.001D0
       tol=1.d-14
c      lista(1)=1
c      lista(2)=2
c      if(e_true.gt.0.001)then
c       do i=3,12
c        lista(i)=i+4
c       enddo
c       do i=13,mxparm
c        lista(i)=0
c       enddo
c       mfit=12
c      else  ! do not fit arg_peri or e_true
c       do i=3,6
c        lista(i)=i+4
c       enddo
c       do i=7,10
c        lista(i)=i+5
c       enddo
c       do i=11,mxparm
c        lista(i)=0
c       enddo
c       mfit=10
c      endif
      endif 

c --- User Input -------------------------------------------------------
      if(.not.mcrun.or.mcfirst)then
       do k=1,mxparm
        lista(k)=0
       enddo
       write(outc,*)
       write(outc,'(a,a)')'Enter ids of parameters to be varied ! ',
     & ' (with blanks):'
       write(outc,'(a,a)')
     & 'e.g.     1   2    3      4      5      6    7   8   9  ',
     & '10 11 12 13 14 15 16'
       write(outc,'(a,a)')
     & 'varies: D1, D2, DR1F2, DR1F3, DR2F2, DR2F3,MD1,MD2,MD3,',
     & ' a, e, P, T, i, O, o'
       data='                                        '//
     &      '                                        '
       read(inc,302)data
  302  format(a80)
       data=' '//data(1:79)
       write(outdat,'(a,a)')'! Parameters varied: ',data(1:60)
       j=1
       do k=1,mxparm
        write(parid2,'(i2)')k
        parid3=parid2//' '
        if(index(data,parid3).ne.0)then
         lista(j)=k
         j=j+1
        endif
       enddo
       mfit=j-1
c
       write(outc,*)
   40  write(outc,'(a,a,$)')' Gaussian (0), double exponential (1),',
     &                      ' or Lorentzian model (2) ? '
       read(inc,*)modef
       if(modef.eq.0)then
        write(outdat,'(a)')'! Gaussian model used in FCHISQ !'
       else if(modef.eq.1)then
        write(outdat,'(a)')'! Double exponential model used in FCHISQ !'
       else if(modef.eq.2)then
        write(outdat,'(a)')'! Lorentzian model used in FCHISQ !'
       else
        write(outc,*)'Model: ',modef,' not allowed !'
        goto 40
       endif
       write(outc,*)
       write(outc,'(a,a,$)')' Enter fractional change in chisq for ',
     & 'convergence (e.g. 0.001): '
       read(inc,*)chifr
      endif
c---- User Input completed ---------------------------------------------
c
c---- Program setup 
c
c Set the model parameter array a to the initial values
      a(1)=diam(1,1)
      a(2)=diam(2,1)
      a(3)=dratio(1,2)
      a(4)=dratio(1,3)
      a(5)=dratio(2,2)
      a(6)=dratio(2,3)
      a(7)=delmag(2,1)
      a(8)=delmag(2,2)
      a(9)=delmag(2,3)
      a(10)=a_true
      a(11)=e_true
      a(12)=period
      a(13)=epoch
      a(14)=inclin
      a(15)=asc_node
      a(16)=arg_peri
c
c     MRQFIT is the only subroutine which writes model parameter values
c     In case MC runs, we still want to write the model out, even if it
c     was not fit to the real data - but since the output file will 
c     contain all the MC results, we should add the model for completeness
c     so that we can use the whole thing as the new binary model file.
c
      if(mcrun.and.mcfirst)then
       write(outdat,1)
       write(outdat,2)(a(i),i=1,9),(mfilter(i),i=1,3)
       write(outdat,5)
       write(outdat,6)(parmerr(j),j=1,9)
       write(outdat,3)
       write(outdat,4)(a(i),i=10,12),a(13)+2440000.D0,(a(j),j=14,16)
       write(outdat,7)
       write(outdat,8)(parmerr(j),j=10,16) 
      endif
c
c Rearrange the data, put modvis2 in xdata. This array is not used
c in this implementation of MRQMIN (MRQCOF passes the index i to
c FUNCS instead). So we can use xdata to pass modvis2 to fchisq.
c
      do ig=1,nfilt
       do id=1,ndata
        ydata((ig-1)*ndata+id)=obsvis2(id,ig)
        sigma((ig-1)*ndata+id)=viserr(id,ig)
        xdata((ig-1)*ndata+id)=modvis2(id,ig)
       enddo
      enddo
      nrec=ndata*nfilt
      call fchisq(ydata,sigma,nrec,nrec,xdata,chisqo,modef)
      chisqo=chisqo*dble(nrec)
      chiini=chisqo
c     Initialize
      alamda=(chisqo/nrec)**6/3000.D0
      alamda=min(alamda,50.D0)
      alamda=max(alamda,0.001D0)
      alamda=-alamda ! Necessary to initialize MRQMIN
      if(adjust)alamda=-2.0
      write(outc,30)chisqo,alamda,nrec
   30 format(' Initial chisq = ',f8.1,', alamda = ',f11.6,
     &       ', # data points = ',i4)
      call mrqmin(xdata,ydata,sigma,nrec,a,mxparm,lista,mfit,
     &            covar,alpha,mxparm,chisq,funcs,alamda,modef,tol)
      IF(.NOT.MCRUN)THEN
       write(outdat,'(a,f8.1)')'! Initial chisq = ',chisqo
      ENDIF
c ----Program Setup completed
c
c- Now start the Marquardt - Levenberg iterations ----------------------

      convergd=.false.
      DO IML=1,MAXIT
       write(outc,'(a,a)')'----------------------------------------'//
     &              '---------------------------------------'
       call mrqmin(xdata,ydata,sigma,nrec,a,mxparm,lista,mfit,
     &             covar,alpha,mxparm,chisq,funcs,alamda,modef,tol)
       write(outc,1)
    1  format('! D1,  D2, DR1F2, DR1F3, DR2F2, DR2F3,  Mag.D1, ',
     &        'Mag.D2, Mag.D3, F1,   F2,   F3')
       write(outc,2)(a(i),i=1,9),(mfilter(i),i=1,3)
    2  format(f6.2,f5.2,1x,4(f4.2,3x),2(f5.2,3x),f5.2,2x,
     &        2(f5.1,1x),f5.1)
       write(outc,3)
    3  format('! A_true,  E_true,   Period, T_periastron, inclin., ',
     &        'Arg_asc.node, Arg_periastron')
       write(outc,4)(a(i),i=10,12),a(13)+2440000.D0,(a(j),j=14,16)
    4  format(1x,f6.2,4x,f6.4,1x,f10.4,2x,f11.3,1x,f7.2,4x,f7.2,
     &        6x,f7.2)
       write(outc,1000)iml,chisq,alamda
 1000  format(' Iteration ',i2,' completed, chisq = ',f8.1,
     &        ', alamda = ',f11.6)

c .... Adjust calibration constants if requested ......................
       if(adjust.and.chisq.lt.chisqo)then
        mode=0
        call binavis(ndata,nfilt,mxrec,mxfilt,date,baseline,hours,
     &               lambda0,dlambda0,modvis2,modphas,udata,vdata,
     &               seprad,parad,mode,ftcurve)
c       call vadjust(modvis2)
        do ig=1,nfilt
         do id=1,ndata
          ydata((ig-1)*ndata+id)=obsvis2(id,ig)
          sigma((ig-1)*ndata+id)=viserr(id,ig)
          xdata((ig-1)*ndata+id)=modvis2(id,ig)
         enddo
        enddo
        call fchisq(ydata,sigma,nrec,nrec,xdata,chisq,modef)
        chisq=chisq*dble(nrec)
       endif

c......Test for convergence............................................
       if(alamda.lt.10)then ! otherwise changes are minuscule anyway
        if(1.D0-chisq/chisqo.lt.chifr.and.chisq.lt.chisqo
     &                                      .and.IML.gt.1)then
         write(outc,*)
         write(outc,1030)chisq,chisq/nrec
 1030    format(' MRQMIN appears to have converged; chisq = ',f8.1,
     &          ', reduced: ',f5.2)
         alamda=0.
         tol=1.d-14
         if(.not.mcrun)
     &   call mrqmin(xdata,ydata,sigma,nrec,a,mxparm,lista,mfit,
     &               covar,alpha,mxparm,chisq,funcs,alamda,modef,tol)
         convergd=.true.
         goto 20
        endif
       endif
       chisqo=chisq

c      The last time this happened, it was because the step size
c      used to calculate the derivatives was too large!
       if(alamda.gt.1000)then
        write(outc,*)
        write(outc,*)'Alamda greater than 1000! Stop now!'
        alamda=0.
        tol=1.d-14
        if(.not.mcrun)
     &  call mrqmin(xdata,ydata,sigma,nrec,a,mxparm,lista,mfit,
     &              covar,alpha,mxparm,chisq,funcs,alamda,modef,tol)
        if(chisq/nrec.lt.2)convergd=.true.  ! we want to have the error bars
        goto 20
       endif

      ENDDO
      write(outc,*)
      write(outc,*)' Max. number of iterations reached; chisq = ',
     &               chisq
c
c ----Iterations have stopped ----------------------------------------
   20 IF(.NOT.MCRUN)THEN
       write(outdat,'(a,f6.1)')'! Final chisq = ',chisq
      ENDIF
      chisqo=chisq

c Set the final model parameters.
      diam(1,1)=abs(a(1))
      diam(2,1)=abs(a(2)) ! sometimes, diameters get negative
      dratio(1,2)=a(3)
      dratio(1,3)=a(4)
      dratio(2,2)=a(5)
      dratio(2,3)=a(6)
      delmag(2,1)=a(7)
      delmag(2,2)=a(8)
      delmag(2,3)=a(9)
      a_true=a(10)
      e_true=a(11)
      period=a(12)
      epoch=a(13)
      inclin=a(14)
      asc_node=a(15)
      arg_peri=a(16)
c Call CADJMOD to fix and update dependent parameters
      call cadjmod
c Take over possible changes in parameters
      a(1)=diam(1,1)
      a(2)=diam(2,1)
      a(3)=dratio(1,2)
      a(4)=dratio(1,3)
      a(5)=dratio(2,2)
      a(6)=dratio(2,3)
      a(7)=delmag(2,1)
      a(8)=delmag(2,2)
      a(9)=delmag(2,3)
      a(10)=a_true
      a(11)=e_true
      a(12)=period
      a(13)=epoch
      a(14)=inclin
      a(15)=asc_node
      a(16)=arg_peri

      if(convergd.and..not.mcrun)then
       chisqo=chisq
       if(chisqo.lt.nrec-mfit)chisqo=nrec-mfit
       do j=1,mxparm
        parmerr(j)=sqrt(covar(j,j)*chisqo/(nrec-mfit))
       enddo
      endif
c
c---- Now that you have finished, write the output----------------------
      IF(.NOT.MCRUN)THEN
       write(outdat,1)     ! Output to file +++++++++
       write(outdat,2)(a(i),i=1,9),(mfilter(i),i=1,3)
       if(convergd)then
        write(outdat,5)
    5   format('! +/- ')
        write(outdat,6)(parmerr(j),j=1,9)
    6   format('!',f5.2,f5.2,1x,4(f4.3,3x),2(f5.3,3x),f5.3,2x,
     &          2(f5.2,1x),f5.2)
       endif
       write(outdat,3)
       write(outdat,4)(a(i),i=10,12),a(13)+2440000.D0,(a(j),j=14,16)
       if(convergd)then
        write(outdat,7)
    7   format('! +/-')
        write(outdat,8)(parmerr(j),j=10,16) 
    8   format('!',f6.3,5x,f5.4,1x,f10.5,2x,f11.4,1x,f6.2,5x,f6.2,
     &          7x,f6.2)
       endif
       write(outc,*)
       write(outc,*)'Hit return...'
       read(inc,*)
      ELSE    ! we need a little bit different output for MONTE CARLO
       if(mcfirst)then
        write(outdat,10)
   10   format('! MONTE CARLO models:')
        write(outdat,11)
   11   format('! D1,   D2, DR1F2, DR1F3, DR2F2, DR2F3,  Mag.D1, ',
     &  'Mag.D2, Mag.D3, I.- F.chi #its')
        write(outdat,12)
   12   format('! A_true,  E_true,   Period, T_periastron, inclin., ',
     &  'Arg_asc.node, Arg_periastron')
        mcfirst=.false.
       endif
       write(outdat,13)(a(i),i=1,9),int(chiini),int(chisq),iml
   13  format(f6.3,f6.3,1x,4(f4.2,3x),2(f6.3,2x),f6.3,
     &        2x,i4,1x,i4,2x,i2)
       write(outdat,14)(a(i),i=10,12),a(13)+2440000.D0,(a(j),j=14,16)
   14  format(f7.3,4x,f6.4,1x,f11.5,2x,f11.3,1x,f7.2,4x,f7.2,
     &        6x,f7.2)
      ENDIF
c
      return
      end