anchor.diff:          write(*,*) 'OUTPUTTING EVOPS'
           call flush()
           call open_files('w')
  -        call write_evops()
  +        call write_evops('w')
           call close_files('w')
   
       end subroutine rates_to_evops
   
  +    subroutine reintegrate()
  +        integer(i4b) :: i,j,tind,u,v,k,l
  +
  +        Evops = 0.0_dp
  +
  +        do i=1,size(Evops,1)
  +        do j=1,size(Evops,2)
  +            Evops(i,j,i,j,1) = 1.0_dp
  +        end do
  +        end do
  +
  +        do tind=2,size(Evops,5)
  +
  +        do i=1,size(Evops,1)
  +        do j=1,size(Evops,2)
  +        do k=1,size(Evops,1)
  +        do l=1,size(Evops,2)
  +
  +            Evops(i,j,k,l,tind) = Evops(i,j,k,l,tind-1)
  +
  +        do u=1,size(Evops,1)
  +        do v=1,size(Evops,2)
  +            Evops(i,j,k,l,tind) = Evops(i,j,k,l,tind) + &
  +                    (DEvops(i,j,u,v,tind-1) + DEvops(i,j,u,v,tind)) / 2.0 * Evops(u,v,k,l,tind-1) * timeStep
  +        end do
  +        end do
  +        end do
  +        end do
  +        end do
  +        end do
  +
  +        end do
  +
  +    end subroutine reintegrate
  +
   end module lindblad_fit_module
   
  
positive.code:         res = res + 10000
        elseif(code == 'D') then
          res = res + 20000
        end if
  
      end function ind
  
negative.code:         name = trim(prefix) // trim(no1) // '-'//trim(no2)//'--'// trim(no3) // '-'//trim(no4)//'.dat'
          open(UNIT=ind(Uelement,Uelement2,Utnemele,Utnemele2,code), FILE = trim(name), STATUS='OLD', ACTION='READ')
  
        else
  
          name = trim(prefix) // trim(no1) // '-'//trim(no2)//'--'// trim(no3) // '-'//trim(no4)//'-out.dat'
          open(UNIT=ind(Uelement,Uelement2,Utnemele,Utnemele2,code), FILE = trim(name))
  