      subroutine p_e_mc() ! bernd berg, jan 14 2002.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      include '../../libs/fortran_par/lat.par'
      include '../../libs/fortran_par/potts.par'
      include '../../libs/fortran_par/mc.par'
      include '../../libs/fortran_par/p_eds.par'
      include '../../libs/fortran_par/common_lat.f'
      include '../../libs/potts_par/common_potts.f'
      include '../../libs/fortran_par/common_p_eds.f'
      dimension nqa(n2d)
      data pcsum/mcase*zero/
      iamin=iact
      iamax=iact
c
      do i=1,ms ! one sweep.
        xr=pmove*rmafun()
        prob=pcase(1)
        icase=1
        if(xr.le.prob) go to 1
        do icase=2,mcam1
          prob=prob+pcase(icase)
          if(xr.le.prob) go to 1
        end do
        icase=mcase
1       continue
        is1=1+int((rmafun()*nacase(icase)))
        is=index(is1,icase)
        iq_in=ista(is)
        xr=rmafun()*qmove(icase)
        ifct=1
        jconf=1
        do id=1,nd
          iqf=ista(ipf(id,is))
          iqb=ista(ipb(id,is))
          nqa(id)=iqf
          nqa(id+nd)=iqb
          jconf=jconf+ifct*(iqf+nq*iqb)
          ifct=nq2*ifct
        end do
        iconf=jconf+iq_in*nconf0
        if(icase.ne.iacase(iconf)) then
          print'(" is,iq_in,nqa:",2i4,6x,6i3)',icase,is,iq_in,nqa
          print*,"icase,iacase(iconf),iconf:",icase,iacase(iconf),iconf
          stop "p_e_mc at 0: icase-iconf mismatch!"
        end if
c q-state update:
        prob=zero
        do iq=0,nqm1
          if(iq.ne.iq_in) then
            ia=0
            do id=1,n2d
              ia=ia+idel(nqa(id),iq)
            end do
            prob=prob+peds1(ia,icase)
            if(xr.le.prob) go to 2
          end if
        end do
        write(iuo,*) "error: iq,xr,prob =",iq,xr,prob
        stop "p_e_mc.f: error with the choice of the new iq."
2       continue
        ista(is)=iq
c table updates:
        jconf=jconf+iq*nconf0
        jcase=iacase(jconf)
        if(icase.ne.jcase) then
          iact=iact+ia-ia_array(icase) ! action change.
          if(iact.lt.0.or.iact.gt.mlink) stop "p_e_mc: iact false!"
          iamin=min(iact,iamin)
          iamax=max(iact,iamax)
          call p_etabs_update(icase,jcase,is,is1)
        end if
        if(jcase.ne.iacase(jconf)) then
          print*,"jcase,iacase(jconf),jconf:",jcase,iacase(jconf),jconf
          stop "p_e_mc at 0: jcase-jconf mismatch!"
        end if
c
        do id0=1,nd ! forward neighbours.
          is0=ipf(id0,is)
          icase=iscase(is0)
          iq0=nqa(id0)
          ifct=1
          jconf=1+iq0*nconf0
          do id=1,nd
            jconf=jconf+ifct*(ista(ipf(id,is0))+nq*ista(ipb(id,is0)))
            ifct=nq2*ifct
          end do
          jcase=iacase(jconf)
          if(icase.ne.jcase) then
            is1=index(is0,mcap1)
            if(index(is1,icase).ne.is0) stop "p_e_mc.f: index mismatch."
            call p_etabs_update(icase,jcase,is0,is1)
          end if
        end do
        if(jcase.ne.iacase(jconf)) then
          print*,"jcase,iacase(jconf),jconf:",jcase,iacase(jconf),jconf
          stop "p_e_mc at 1: jcase-jconf mismatch!"
        end if
c
        do id0=1,nd ! backward neighbours.
          is0=ipb(id0,is)
          icase=iscase(is0)
          iq0=nqa(id0+nd)
          ifct=1
          jconf=1+iq0*nconf0
          do id=1,nd
            jconf=jconf+ifct*(ista(ipf(id,is0))+nq*ista(ipb(id,is0)))
            ifct=nq2*ifct
          end do
          jcase=iacase(jconf)
          if(icase.ne.jcase) then
            is1=index(is0,mcap1)
            if(index(is1,icase).ne.is0) stop "p_e_mc: index mismatch."
	    call p_etabs_update(icase,jcase,is0,is1)
          end if
        end do
	if(jcase.ne.iacase(jconf)) then
	  print*,"jcase,iacase(jconf),jconf:",jcase,iacase(jconf),jconf
	  stop "p_e_mc at 2: jcase-jconf mismatch!"
        end if
	ha(iact)=ha(iact)+one/(one-pstay)
      end do
c
      ia_min=min(ia_min,iamin)
      ia_max=min(ia_max,iamax)
      pmove=zero
      do icase=1,mcase
        pmove=pmove+pcase(icase) ! re-calculation to avoid rounding.
        pcsum(icase)=pcsum(icase)+pcase(icase) ! for pcase_sort call.
      end do
      pstay=one-pmove
      return
      end


      subroutine p_etabs_update(icase,jcase,is0,is1)
c bernd berg, jan 14 2002. update eds tables.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      include '../../libs/fortran_par/lat.par'
      include '../../libs/fortran_par/potts.par'
      include '../../libs/fortran_par/p_eds.par'
      include '../../libs/fortran_par/common_lat.f'
      include '../../libs/potts_par/common_potts.f'
      include '../../libs/fortran_par/common_p_eds.f'
      is2=index(nacase(icase),icase)
      index(is1,icase)=is2
      index(is2,mcap1)=is1
      nacase(icase)=nacase(icase)-1
      pstay=pstay+pcase(icase)                ! sensitive to rounding.
      pcase(icase)=nacase(icase)*qmove(icase) ! avoids rounding.
      pstay=pstay-pcase(icase)
      nap1=nacase(jcase)+1
      nacase(jcase)=nap1
      pstay=pstay+pcase(jcase)
      pcase(jcase)=nacase(jcase)*qmove(jcase)
      pstay=pstay-pcase(jcase)
      pmove=one-pstay
      index(nap1,jcase)=is0
      index(is0,mcap1)=nap1
      iscase(is0)=jcase
      return
      end


      subroutine pcase_sort ! bernd berg, jan 14 2002.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      include '../../libs/fortran_par/lat.par'
      include '../../libs/fortran_par/potts.par'
      include '../../libs/fortran_par/p_eds.par'
      include '../../libs/fortran_par/common_lat.f'
      include '../../libs/potts_par/common_potts.f'
      include '../../libs/fortran_par/common_p_eds.f'
      dimension icase_work(mcase)
c
      call heap_per(mcase,pcsum,icase_per)
c
      np1=mcase+1
      do icase=1,mcase
        pcsum(icase)=pcase(icase_per(np1-icase))
        icase_work(icase)=nacase(icase_per(np1-icase))
      end do
      do icase=1,mcase
        pcase(icase)=pcsum(icase)
        nacase(icase)=icase_work(icase)
      end do
c
      do icase=1,mcase
        pcsum(icase)=qmove(icase_per(np1-icase))
        icase_work(icase)=ia_array(icase_per(np1-icase))
      end do
      do icase=1,mcase
        qmove(icase)=pcsum(icase)
        ia_array(icase)=icase_work(icase)
      end do
c
      do id=0,n2d
        do icase=1,mcase
          pcsum(icase)=peds1(id,icase_per(np1-icase))
        end do
        do icase=1,mcase
          peds1(id,icase)=pcsum(icase)
        end do
      end do
c
      do is=1,ms
        do icase=1,mcase
          icase_work(icase)=index(is,icase_per(np1-icase))
        end do
        do icase=1,mcase
          index(is,icase)=icase_work(icase)
        end do
      end do
c
      do icase=1,mcase
        do is1=1,nacase(icase)
          is=index(is1,icase)
          index(is,np1)=is1
        end do
      end do
c
      do icase=1,mcase
        icase_work(icase_per(np1-icase))=icase
      end do
      do iconf=1,nconf
        iacase(iconf)=icase_work(iacase(iconf))
      end do
      do is=1,ms
        iscase(is)=icase_work(iscase(is))
      end do
c
      do icase=1,mcase
        pcsum(icase)=zero
      end do
      return
      end
