anchor.diff:                    inversion_mode, &
                     fault, &
                     rake_constraint, &
  +                  displacement, &
  +                  disp_components, &
  +                  disp_misfit_file, &
  +                  gf_disp, &
                     fault_slip
   
   implicit none
   
   ! Local variables
  -integer :: i, ounit
  -double precision :: slip_mag
  +integer :: i, ierr, ounit, icmp, idsp, iflt, ndsp, nflt, nobs
  +double precision :: slip_mag, rms
  +double precision, allocatable :: pre(:), obs(:)
   
   
   if (verbosity.ge.1) then
       write(stdout,*) 'write_solution: starting'
   endif
  -!
  -! ! Print RMS misfit if specified
  -! if (disp_misfit_file.ne.'none') then
  -!     if (verbosity.ge.1) then
  -!         write(stderr,'(A)') 'write_solution says: writing displacement RMS misfit to '// &
  -!                             trim(disp_misfit_file)
  -!     endif
  -!
  -!     write(0,*) 'write_solution: opening misfit file'
  -!     open(unit=81,file=disp_misfit_file,status='unknown')
  -!
  -!     ! If rake is constrained, make an array of the correct size to use with misfit function
  -!     if (inversion_mode.eq.'lsqr'.and.rake_constraint%file.ne.'none' &
  -!                                                            .and.rake_constraint%nfields.eq.1) then
  -!         write(0,*) 'write_solution: writing misfit for fixed rake'
  -!         do i = 1,fault%nrecords
  -!             tmp_slip_array(i,1) = fault_slip(i,1) ! Green's functions already calculated for this rake
  -!             tmp_slip_array(i,2) = 0.0d0
  -!         enddo
  -!         write(81,*) disp_misfit_l2norm(tmp_slip_array)/dsqrt(dble(fault%nrecords))
  -!
  -!     ! Otherwise, just use the misfit function directly
  -!     else
  -!         write(0,*) 'write_solution: writing misfit for free rake'
  -!         write(81,*) disp_misfit_l2norm(fault_slip)/dsqrt(dble(fault%nrecords))
  -!     endif
  -!
  -!     write(0,*) 'write_solution: closing misfit file'
  -!     close(81)
  -! endif
  -!
  +
  +
  +! Print RMS misfit to three-component displacements
  +if (disp_misfit_file.ne.'none') then
  +
  +    ! Open misfit file
  +    open(unit=81,file=disp_misfit_file,status='unknown')
  +
  +    ndsp = displacement%nrows
  +    nobs = len_trim(disp_components)*displacement%nrows
  +    nflt = fault%nrows
  +
  +    ! Load predicted three-component displacements
  +    if (.not.allocated(pre)) then
  +        allocate(pre(ndsp*len_trim(disp_components)),stat=ierr)
  +        if (ierr.ne.0) then
  +            call usage('write_solution: error allocating memory to pre')
  +        endif
  +    endif
  +    pre = 0.0d0
  +
  +    do i = 1,len_trim(disp_components)
  +        read(disp_components(i:i),*) icmp
  +        do idsp = 1,ndsp
  +            do iflt = 1,nflt
  +                pre((i-1)*ndsp+idsp) = pre((i-1)*ndsp+idsp) + &
  +                                       gf_disp%array((icmp-1)*ndsp+idsp,     iflt)*fault_slip(iflt,1)
  +                pre((i-1)*ndsp+idsp) = pre((i-1)*ndsp+idsp) + &
  +                                       gf_disp%array((icmp-1)*ndsp+idsp,nflt+iflt)*fault_slip(iflt,2)
  +            enddo
  +        enddo
  +    enddo
  +
  +    ! Load observed three-component displacements
  +    if (.not.allocated(obs)) then
  +        allocate(obs(nobs),stat=ierr)
  +        if (ierr.ne.0) then
  +            call usage('write_solution: error allocating memory to obs')
  +        endif
  +    endif
  +    obs = 0.0d0
  +
  +    do i = 1,len_trim(disp_components)
  +        read(disp_components(i:i),*) icmp
  +        obs((i-1)*ndsp+1:i*ndsp) = displacement%array(1:ndsp,3+icmp)
  +    enddo
  +
  +    call misfit_rms(obs,pre,nobs,rms)
  +    write(81,*) rms
  +
  +    close(81)
  +
  +endif
  +
   ! ! Line-of-sight RMS misfit
   ! if (los_misfit_file.ne.'none') then
   !     if (verbosity.ge.1) then
  
positive.code:         endif
          rake_constraint%nrows = fault%nrows
          rake_constraint%array(:,1) = dp1
          rake_constraint%array(:,2) = dp2
      endif
  else
      if (inversion_mode.eq.'anneal') then
  
negative.code: write(stderr,*) 'Input Options'
  write(stderr,*) '-disp DISP_FILE              Input displacements'
  write(stderr,*) '-disp:components COMPNTS     Specify displacement components'
  write(stderr,*) '-disp:misfit MISFIT_FILE     Output RMS misfit to displacements'
  write(stderr,*) '-los LOS_FILE                Input line-of-sight displacements'
  write(stderr,*) '-prests PRESTS_FILE          Input pre-stresses'
  write(stderr,*) '-cov COVAR_FILE              Displacement and LOS covariances'
  ! write(stderr,*) '-prests:weight WEIGHT        Stress weighting factor'
  ! write(stderr,*) '-prests:dist_threshold DIST  Set tractions to zero at distances>DIST'
  ! write(stderr,*) '-los:weight LOS_WEIGHT       LOS observation weighting factor'
  