text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c MINT Integrator Package
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Original version by Paolo Nason (for POWHEG (BOX))
c Modified by Rikkert Frederix (for aMC@NLO)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      subroutine mint(fun,ndim,ncalls0,nitmax,imode,
c ndim=number of dimensions
c ncalls0=# of calls per iteration
c nitmax =# of iterations
c fun(xx,www,ifirst): returns the function to be integrated multiplied by www;
c                     xx(1:ndim) are the variables of integration
c                     ifirst=0: normal behaviour
c imode: integer flag
c
c imode=-1:
c same as imode=0 as far as this routine is concerned, except for the
c fact that a grid is read at the beginning (rather than initialized).
c The return value of imode will be zero.
c
c imode=0:
c When called with imode=0 the routine integrates the absolute value of
c the function and sets up a grid xgrid(0:50,ndim) such that in each
c ndim-1 dimensional slice (i.e. xgrid(m-1,n)<xx(n)<xgrid(m,n)) the
c contribution of the integral is the same the array xgrid is setup at
c this stage; ans and err are the integral and its error
c
c imode=1 (in fact #0)
c When called with imode=1, the routine performs the integral of the
c function fun using the grid xgrid. If some number in the array ifold,
c (say, ifold(n)) is different from 1, it must be a divisor of 50, and
c the 50 intervals xgrid(0:50,n) are grouped into ifold(n) groups, each
c group containing 50/ifold(n) nearby intervals. For example, if
c ifold(1)=5, the 50 intervals for the first dimension are divided in 5
c groups of 10. The integral is then performed by folding on top of each
c other these 5 groups. Suppose, for example, that we choose a random
c point in xx(1) = xgrid(2,1)+x*(xgrid(3,1)-xgrid(2,1)), in the group of
c the first 5 interval.  we sum the contribution of this point to the
c contributions of points
c xgrid(2+m*10,1)+x*(xgrid(3+m*10,1)-xgrid(2+m*10,1)), with m=1,...,4.
c In the sequence of calls to the function fun, the call for the first
c point is performed with ifirst=0, and that for all subsequent points
c with ifirst=1, so that the function can avoid to compute quantities
c that only depend upon dimensions that have ifold=1, and do not change
c in each group of folded call. The values returned by fun in a sequence
c of folded calls with ifirst=0 and ifirst=1 are not used. The function
c itself must accumulate the values, and must return them when called
c with ifirst=2.
c
c Added the posibility to keep track of more than one integral:

c nintegrals=1 : the function that is used to update the grids. This is
c the ABS cross section. If imode.eq.1, this does not contain the
c virtual corrections because for them a separate maximum is kept using (5).
c nintegrals=2 : the actual cross section. This includes virtual corrections.
c nintegrals=3 : the cross section from the M_Virt/M_Born ratio alone:
c this defines the average virtual that is added to each phase-space
c point
c nintegrals=4 : the cross section of the actual virtual minus the
c average virtual. This is used to determine the fraction of phase-space
c points for which we include the virtual.
c nintegrals=5 : abs of 3
c nintegrals=6 : born*alpha_S/2Pi
c
      subroutine mint(fun,ndim,ncalls0,nitmax,imode,xgrid,ymax,ymax_virt
     $     ,ans,unc,chi2)
c imode= 0: integrate and adapt the grid
c imode= 1: frozen grid, compute the integral and the upper bounds
c imode=-1: same as imode=0, but use previously generated grids
c others: same as 1 (for now)
      implicit none
      include "mint.inc"
      include "FKSParams.inc"
      integer i,j,ncalls0,ndim,nitmax,imode
      real * 8 fun,xgrid(0:nintervals,ndimmax),xint,ymax(nintervals
     $     ,ndimmax),ans(nintegrals),unc(nintegrals),ans3(nintegrals,3)
     $     ,unc3(nintegrals,3),ans_l3(nintegrals),unc_l3(nintegrals)
     $     ,chi2_l3(nintegrals)
      real * 8 xint_virt,ymax_virt
      real * 8 x(ndimmax),vol
      real * 8 xacc(0:nintervals,ndimmax)
      integer icell(ndimmax),ncell(ndimmax),ncell_virt
      integer ifold(ndimmax),kfold(ndimmax)
      common/cifold/ifold
      integer nhits(nintervals,ndimmax)
      real * 8 rand(ndimmax)
      real * 8 dx(ndimmax),f(nintegrals),vtot(nintegrals)
     $     ,etot(nintegrals),prod,f1(nintegrals),chi2(nintegrals)
     $     ,efrac(nintegrals),dummy
      integer kdim,kint,kpoint,nit,ncalls,iret,nintcurr,nintcurr_virt
     $     ,ifirst,nit_included,kpoint_iter,non_zero_point(nintegrals)
     $     ,ntotcalls(nintegrals),nint_used,nint_used_virt,min_it
      real * 8 ran3
      external ran3,fun
      logical even,double_events,bad_iteration
      double precision average_virtual,virtual_fraction
      common/c_avg_virt/average_virtual,virtual_fraction
      character*13 title(nintegrals)
      data title(1)/'ABS integral '/
      data title(2)/'Integral     '/
      data title(3)/'Virtual      '/
      data title(4)/'Virtual ratio'/
      data title(5)/'ABS virtual  '/
      data title(6)/'Born*ao2pi   '/
c APPLgrid switch
      integer iappl
      common /for_applgrid/ iappl
      logical              fixed_order,nlo_ps
      common /c_fnlo_nlops/fixed_order,nlo_ps
c if ncalls0 is greater than 0, use the default running, i.e. do not
c double the events after each iteration as well as use a fixed number
c of intervals in the grids.
      if (ncalls0.gt.0) then
         double_events=.false.
         nint_used=nintervals
         nint_used_virt=nintervals_virt
      else
c if ncalls0.le.0, reset it and double the events per iteration
         ncalls0=80*ndim
         double_events=.true.
         if (imode.eq.1 .or. imode.eq.-1) then
            nint_used=nintervals
            nint_used_virt=nintervals_virt
         else
            nint_used=min_inter
            nint_used_virt=min_inter
         endif
      endif
      bad_iteration=.false.
c
      ncalls=0  ! # PS points (updated below)
      if(imode.eq.-1) then
c Grids read from file
         even=.true.
         imode=0
         min_it=min_it0
         do kdim=1,ndim
            ifold(kdim)=1
         enddo
      elseif(imode.eq.0) then
c Initialize grids
         even=.true.
         min_it=min_it0
         do kdim=1,ndim
            ifold(kdim)=1
            do kint=0,nint_used
               xgrid(kint,kdim)=dble(kint)/nint_used
            enddo
         enddo
         call init_ave_virt(nint_used_virt,ndim)
      elseif(imode.eq.1) then
c Initialize upper bounding envelope
         xint=ans(1)
         xint_virt=ans(5)
         even=.false.
         min_it=min_it1
         do kdim=1,ndim
            nintcurr=nint_used/ifold(kdim)
            nintcurr_virt=nint_used_virt/ifold(kdim)
            if(nintcurr*ifold(kdim).ne.nint_used .or.
     &         nintcurr_virt*ifold(kdim).ne.nint_used_virt) then
               write(*,*) 'mint: the values in the ifold array'/
     $              /'shoud be divisors of',nint_used,'and'
     $              ,nint_used_virt
               stop
            endif
            do kint=1,nintcurr
               ymax(kint,kdim)=xint**(1d0/ndim)
            enddo
         enddo
         ymax_virt=xint_virt
      endif
      nit=0
      nit_included=0
      do i=1,nintegrals
         ans(i)=0d0
         unc(i)=0d0
         do j=1,3
            ans3(i,j)=0d0
            unc3(i,j)=0d0
         enddo
      enddo
c Main loop over the iterations
 10   continue
c This makes sure that current stdout is written to log.txt and cache is
c emptied.
      call flush(6)
      if(nit.ge.nitmax) then
c We did enough iterations, update arguments and return
         if(imode.eq.0) xint=ans(1)
         if(imode.eq.0) xint_virt=ans(5)
         if (nit_included.ge.2) then
            chi2(1)=chi2(1)/dble(nit_included-1)
         else
            chi2(1)=0d0
         endif
         write (*,*) '-------'
         ncalls0=ncalls*kpoint_iter ! return number of points used
         if (double_events) then
            nitmax=2
         else
            nitmax=nit_included
         endif
         return
      endif
      nit=nit+1
      write (*,*) '------- iteration',nit
      if (even .and. ncalls.ne.ncalls0) then
c Uses more evenly distributed random numbers. This overwrites the
c number of calls
         call initialize_even_random_numbers(ncalls0,ndim,ncalls)
         write (*,*) 'Update # PS points (even): ',ncalls0,' --> '
     &        ,ncalls
      elseif (ncalls0.ne.ncalls) then
         ncalls=ncalls0
         write (*,*) 'Update # PS points: ',ncalls0,' --> ',ncalls
      endif
c Reset the accumulated results for grid updating
      if(imode.eq.0) then
         do kdim=1,ndim
            do kint=0,nint_used
               xacc(kint,kdim)=0
               if(kint.gt.0) then
                  nhits(kint,kdim)=0
               endif
            enddo
         enddo
      endif
      do i=1,nintegrals
         vtot(i)=0
         etot(i)=0
      enddo
      kpoint_iter=0
      do i=1,nintegrals
         non_zero_point(i)=0
      enddo
c Loop over PS points
 2    kpoint_iter=kpoint_iter+1
      do kpoint=1,ncalls
c find random x, and its random cell
         do kdim=1,ndim
            kfold(kdim)=1
c if(even), we should compute the ncell and the rand from the ran3()
            if (even) then
               rand(kdim)=ran3(even)
               ncell(kdim)= min(int(rand(kdim)*nint_used)+1,
     &              nint_used)
               rand(kdim)=rand(kdim)*nint_used-(ncell(kdim)-1)
            else
               ncell(kdim)=min(int(nint_used/ifold(kdim)*ran3(even))+1,
     &              nint_used)
               rand(kdim)=ran3(even)
            endif
         enddo
         do i=1,nintegrals
            f(i)=0
         enddo
         ifirst=0
 1       continue
         vol=1
c compute jacobian ('vol') for the PS point
         do kdim=1,ndim
            nintcurr=nint_used/ifold(kdim)
            icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
            dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
            vol=vol*dx(kdim)*nintcurr
            x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
            if(imode.eq.0)
     &           nhits(icell(kdim),kdim)=nhits(icell(kdim),kdim)+1
         enddo
         call get_ave_virt(x,nint_used_virt,ndim,average_virtual)
c contribution to integral
         if(imode.eq.0) then
            dummy=fun(x,vol,ifirst,f1)
            do i=1,nintegrals
               f(i)=f(i)+f1(i)
            enddo
         else
c this accumulated value will not be used
            dummy=fun(x,vol,ifirst,f1)
            do i=1,nintegrals
               f(i)=f(i)+f1(i)
            enddo
            ifirst=1
            call nextlexi(ndim,ifold,kfold,iret)
            if(iret.eq.0) goto 1
c closing call: accumulated value with correct sign
            dummy=fun(x,vol,2,f1)
            do i=1,nintegrals
               f(i)=f1(i)
            enddo
         endif
c
         if(imode.eq.0) then
c accumulate the function in xacc(icell(kdim),kdim) to adjust the grid later
            do kdim=1,ndim
               xacc(icell(kdim),kdim)=xacc(icell(kdim),kdim)+f(1)
            enddo
c Set the Born contribution (to compute the average_virtual) to zero if
c the virtual was not computed for this phase-space point. Compensate by
c including the virtual_fraction.
            if (f(3).ne.0d0) then
               f(6)=f(6)/virtual_fraction
               call fill_ave_virt(x,nint_used_virt,ndim,f(3),f(6))
            else
               f(6)=0d0
            endif
         else
c update the upper bounding envelope total rate
            prod=1d0
            do kdim=1,ndim
               prod=prod*ymax(ncell(kdim),kdim)
            enddo
            prod=(f(1)/prod)
            if (prod.gt.1d0) then
c Weight for this PS point is larger than current upper bound. Increase
c the bound so that it is equal to the current max weight.  If the new
c point is more than twice as large as current upper bound, increase
c bound by factor 2 only to prevent a single unstable points to
c completely screw up the efficiency
               prod=min(2d0,prod)
               prod=prod**(1d0/dble(ndim))
               do kdim=1,ndim
                  ymax(ncell(kdim),kdim)=ymax(ncell(kdim),kdim)*prod
               enddo
            endif
c Update the upper bounding envelope virtual. Do not include the
c enhancement due to the virtual_fraction. (And again limit by factor 2
c at most).
            if (f(5)*virtual_fraction.gt.ymax_virt) ymax_virt=min(f(5)
     $           *virtual_fraction,ymax_virt*2d0)
c for consistent printing in the log files (in particular when doing LO
c runs), set also f(6) to zero when imode.eq.1 and the virtuals are not
c included.
            if (f(3).eq.0) f(6)=0d0
         endif
         do i=1,nintegrals
            if (f(i).ne.0d0) non_zero_point(i)=non_zero_point(i)+1
         enddo
c Add the PS point to the result of this iteration
         do i=1,nintegrals
            vtot(i)=vtot(i)+f(i)
            etot(i)=etot(i)+f(i)**2
         enddo
         if (f(1).ne.0d0) call HwU_add_points
      enddo
      do i=1,nintegrals
c Number of phase-space points used
         ntotcalls(i)=ncalls*kpoint_iter
c Special for the computation of the 'computed virtual'
         if (i.eq.4 .and. non_zero_point(i).ne.0 )
     &        ntotcalls(i) = non_zero_point(i)
      enddo
      if (ntotcalls(1).gt.max_points .and. non_zero_point(1).lt.25 .and.
     &     double_events) then
         write (*,*) 'ERROR: INTEGRAL APPEARS TO BE ZERO.'
         write (*,*) 'TRIED',ntotcalls(1),'PS POINTS AND ONLY '
     &        ,non_zero_point(1),' GAVE A NON-ZERO INTEGRAND.'
         stop 1
      endif
c Goto beginning of loop over PS points until enough points have found
c that pass cuts.
      if (non_zero_point(1).lt.ncalls .and. double_events) goto 2

c Iteration done. Update the accumulated results and print them to the
c screen
      do i=1,nintegrals
         vtot(i)=vtot(i)/dble(ntotcalls(i))
         etot(i)=etot(i)/dble(ntotcalls(i))
c the abs is to avoid tiny negative values
         etot(i)=sqrt(abs(etot(i)-vtot(i)**2)
     $        /dble(ntotcalls(i)))
         if (vtot(i).ne.0d0) then
            efrac(i)=abs(etot(i)/vtot(i))
         else
            efrac(i)=0d0
         endif
      enddo
      do i=1,nintegrals
         write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
     $        title(i)//' =',vtot(i),' +/- ',etot(i),' (',efrac(i)*100d0
     $        ,'%)'
      enddo
C If there was a large fluctation in this iteration, be careful with
C including it in the accumalated results and plots.
      if (efrac(1).gt.0.3d0 .and. iappl.eq.0) then
c Do not include the results in the plots
         if (fixed_order) call accum(.false.)
         if (fixed_order) call HwU_accum_iter(.false.,ntotcalls(1))
      endif
      if (efrac(1).gt.0.3d0 .and. nit.gt.3 .and. iappl.eq.0) then
c Do not include the results in the updating of the grids.
         write (*,*) 'Large fluctuation ( >30 % ).'
     &        //'Not including iteration in results.'
c empty the accumulated results in the MC over integers
         call empty_MC_integer
c empty the accumalated results for the MINT grids
         if (imode.eq.0) then
c emptying accum. results is done above when the iteration starts
            continue
         elseif (imode.eq.1) then
c Cannot really skip the increase of the upper bounding envelope. So,
c simply continue here. Note that no matter how large the integrand for
c the PS point, the upper bounding envelope is at most increased by a
c factor 2, so this should be fine.
            continue
         endif
c double the number of points for the next iteration
         if (double_events) ncalls0=ncalls0*2
         if (bad_iteration .and. imode.eq.0 .and. double_events) then
c 2nd bad iteration is a row. Reset grids
            write (*,*)'2nd bad iteration in a row. '/
     &           /'Resetting grids and starting from scratch...'
            if (double_events) then
               if (imode.eq.0) nint_used=min_inter ! reset number of intervals
               ncalls0=ncalls0/8   ! Start with larger number
            endif
            nit=0
            nit_included=0
c Reset the MINT grids
            if (imode.eq.0) then
               do kdim=1,ndim
                  do kint=0,nint_used
                     xgrid(kint,kdim)=dble(kint)/nint_used
                  enddo
               enddo
               call init_ave_virt(nint_used_virt,ndim)
            elseif (imode.eq.1) then
               do kdim=1,ndim
                  nintcurr=nint_used/ifold(kdim)
                  do kint=1,nintcurr
                     ymax(kint,kdim)=xint**(1d0/ndim)
                  enddo
                  nintcurr_virt=nint_used_virt/ifold(kdim)
               enddo
               ymax_virt=xint_virt
            endif
            call reset_MC_grid  ! reset the grid for the integers
            if (fixed_order) call initplot  ! Also reset all the plots
            do i=1,nintegrals
               ans(i)=0d0
               unc(i)=0d0
               chi2(i)=0d0
               do j=1,3
                  ans3(i,j)=0d0
                  unc3(i,j)=0d0
               enddo
            enddo
            bad_iteration=.false.
         else
            bad_iteration=.true.
         endif
         goto 10
      else
         bad_iteration=.false.
      endif
      if(nit.eq.1) then
         do i=1,nintegrals
            ans(i)=vtot(i)
            unc(i)=etot(i)
         enddo
         write (*,'(a,1x,e10.4)') 'Chi^2 per d.o.f.',0d0
      else
c prevent annoying division by zero for nearly zero
c integrands
         do i=1,nintegrals
            if(etot(i).eq.0.and.unc(i).eq.0) then
               if(ans(i).eq.vtot(i) .and. i.eq.1) then
c double the number of points for the next iteration
                  if (double_events) ncalls0=ncalls0*2
                  goto 10
               else
                  unc(i)=abs(vtot(i)-ans(i))
                  etot(i)=abs(vtot(i)-ans(i))
               endif
            elseif(etot(i).eq.0) then
               etot(i)=unc(i)
            elseif(unc(i).eq.0) then ! 1st iteration; set to a large value
               unc(i)=etot(i)*1d99
            endif
            if (i.ne.1 .and. (etot(i).eq.0 .or. unc(i).eq.0)) then
               ans(i)=0d0
               unc(i)=0d0
               chi2(i)=0d0
            else
               ans(i)=(ans(i)/unc(i)+vtot(i)/etot(i))/
     &              (1/unc(i)+1/etot(i))
               unc(i)=1/sqrt(1/unc(i)**2+1/etot(i)**2)
               chi2(i)=chi2(i)+(vtot(i)-ans(i))**2/etot(i)**2
            endif
         enddo
         write (*,'(a,1x,e10.4)') 'Chi^2=',(vtot(1)-ans(1))**2
     $        /etot(1)**2
      endif

      nit_included=nit_included+1
      do i=1,nintegrals
         if (ans(i).ne.0d0) then
            efrac(i)=abs(unc(i)/ans(i))
         else
            efrac(i)=0d0
         endif
         write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
     $        'accumulated results '//title(i)//' =',ans(i),' +/- '
     $        ,unc(i) ,' (',efrac(i)*100d0,'%)'
      enddo
      if (nit_included.le.1) then
         write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF ='
     $        ,0d0
      else
         write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF =',
     &        chi2(1)/dble(nit_included-1)
      endif
      if (imode.eq.0) then
c     Update the average_virtual: a_new=(virt+a_old*born)/born
         if (vtot(6).ne.0d0) then
            if (average_virtual.eq.0d0) then ! i.e. first iteration
               average_virtual=vtot(3)/vtot(6)+average_virtual
            else  ! give some importance to the iterations already done
               average_virtual=(vtot(3)/vtot(6)+average_virtual*2d0)/2d0
            endif
         endif
c Update the fraction of the events for which we include the virtual corrections
c in the calculation
         virtual_fraction=max(min(virtual_fraction*max(min(2d0*etot(3)
     $        /etot(1),2d0),0.25d0),1d0),Min_virt_fraction)
         write (*,'(a,1x,f7.3,1x,f7.3)') 'update virtual fraction to:'
     $        ,virtual_fraction,average_virtual
      elseif (imode.eq.1) then
         write (*,'(a,1x,f7.3,1x,f7.3)') 'virtual fraction is:'
     $        ,virtual_fraction,average_virtual
      endif
c Update the results of the last tree iterations
      do i=1,nintegrals
         do j=1,2
            ans3(i,j)=ans3(i,j+1)
            unc3(i,j)=unc3(i,j+1)
         enddo
         ans3(i,3)=vtot(i)
         unc3(i,3)=etot(i)
      enddo
c Compute the results of the last three iterations
      if (nit_included.ge.4) then
         do i=1,nintegrals
            ans_l3(i)=0d0
            unc_l3(i)=ans3(i,1)*1d99
            chi2_l3(i)=0d0
            do j=1,3
               if (i.ne.1 .and. (unc_l3(i).eq.0d0 .or. unc3(i
     $              ,j).eq.0d0)) then
                  ans_l3(i)=0d0
                  unc_l3(i)=0d0
                  chi2_l3(i)=0d0
               else
                  ans_l3(i)=(ans_l3(i)/unc_l3(i)+ans3(i,j)/unc3(i,j))/
     &                 (1/unc_l3(i)+1/unc3(i,j))
                  unc_l3(i)=1/sqrt(1/unc_l3(i)**2+1/unc3(i,j)**2)
                  chi2_l3(i)=chi2_l3(i)+
     &                 (ans3(i,j)-ans_l3(i))**2/unc3(i,j)**2
               endif
            enddo
            chi2_l3(i)=chi2_l3(i)/2d0 ! three iterations, so 2 degrees of freedom
         enddo
         do i=1,2
            if (ans_l3(i).ne.0d0) then
               efrac(i)=abs(unc_l3(i)/ans_l3(i))
            else
               efrac(i)=0d0
            endif
            write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
     $           'accumulated results last 3 iterations '//title(i)/
     $           /' =' ,ans_l3(i),' +/- ',unc_l3(i) ,' (',efrac(i)*100d0
     $           ,'%)'
         enddo
         write(*,'(a,1x,e10.4)')
     $        'accumulated result last 3 iterrations Chi^2 per DoF ='
     $        ,chi2_l3(1)
      endif
      if(imode.eq.0) then
c Iteration is finished; now rearrange the grid
         do kdim=1,ndim
            call regrid(xacc(0,kdim),xgrid(0,kdim),nhits(1,kdim)
     $           ,nint_used)
         enddo
         call regrid_ave_virt(nint_used_virt,ndim)
c Regrid the MC over integers (used for the MC over FKS dirs)
         call regrid_MC_integer
      endif
c Quit if the desired accuracy has been reached
      if (nit_included.ge.min_it .and. accuracy.gt.0d0) then
         if (unc(1)/ans(1)*max(1d0,chi2(1)/dble(nit_included-1))
     $        .lt.accuracy) then
            write (*,*) 'Found desired accuracy'
            nit=nitmax
c Improve the stats in the plots
            if (fixed_order) call accum(.true.)
            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
            goto 10
         elseif(unc_l3(1)/ans_l3(1)*max(1d0,chi2_l3(1)).lt.accuracy)
     $           then
            write (*,*)
     &           'Found desired accuracy in last 3 iterations'
            nit=nitmax
            do i=1,nintegrals
               ans(i)=ans_l3(i)
               unc(i)=unc_l3(i)
               chi2(i)=chi2_l3(i)*dble(nit_included-1)
            enddo
c Improve the stats in the plots
            if (fixed_order) call accum(.true.)
            if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
            goto 10
         endif
      endif
c Double the number of intervals in the grids if not yet reach the maximum
      if (2*nint_used.le.nintervals .and. double_events) then
         do kdim=1,ndim
            call double_grid(xgrid(0,kdim),nint_used)
         enddo
         nint_used=2*nint_used
      endif
      if (2*nint_used_virt.le.nintervals_virt .and. double_events) then
         call double_ave_virt(nint_used_virt,ndim)
         nint_used_virt=2*nint_used_virt
      endif

c double the number of points for the next iteration
      if (double_events) ncalls0=ncalls0*2
c Also improve stats in plots
      if (fixed_order) call accum(.true.)
      if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1))
c Do next iteration
      goto 10
      end

      subroutine double_grid(xgrid,ninter)
      implicit none
      include "mint.inc"
      integer  ninter
      real * 8 xgrid(0:nintervals)
      integer i
      do i=ninter,1,-1
         xgrid(i*2)=xgrid(i)
         xgrid(i*2-1)=(xgrid(i)+xgrid(i-1))/2d0
      enddo
      return
      end

      subroutine regrid(xacc,xgrid,nhits,ninter)
      implicit none
      include "mint.inc"
      integer  ninter,nhits(nintervals)
      real * 8 xacc(0:nintervals),xgrid(0:nintervals)
      real * 8 xn(nintervals),r,tiny,xl,xu,nl,nu
      parameter ( tiny=1d-8 )
      integer kint,jint
      logical plot_grid
      parameter (plot_grid=.false.)
c Use the same smoothing as in VEGAS uses for the grids (i.e. use the
c average of the central and the two neighbouring grid points):
      xl=xacc(1)
      xu=xacc(2)
      xacc(1)=(xl+xu)/2d0
      nl=nhits(1)
      nu=nhits(2)
      nhits(1)=nint((nl+nu)/2d0)
      do kint=2,ninter-1
         xacc(kint)=xl+xu
         xl=xu
         xu=xacc(kint+1)
         xacc(kint)=(xacc(kint)+xu)/3d0
         nhits(kint)=nl+nu
         nl=nu
         nu=nhits(kint+1)
         nhits(kint)=nint((nhits(kint)+nu)/3d0)
      enddo
      xacc(ninter)=(xu+xl)/2d0
      nhits(ninter)=nint((nu+nl)/2d0)
c
      do kint=1,ninter
c xacc (xerr) already contains a factor equal to the interval size
c Thus the integral of rho is performed by summing up
         if(nhits(kint).ne.0) then
            xacc(kint)= xacc(kint-1)
     #           + abs(xacc(kint))/nhits(kint)
         else
            xacc(kint)=xacc(kint-1)
         endif
      enddo
      do kint=1,ninter
         xacc(kint)=xacc(kint)/xacc(ninter)
      enddo
cRF: Check that we have a reasonable result and update the accumulated
c results if need be
      do kint=1,ninter
         if (xacc(kint).lt.(xacc(kint-1)+tiny)) then
c            write (*,*) 'Accumulated results need adaptation #1:'
c            write (*,*) xacc(kint),xacc(kint-1),' become'
            xacc(kint)=xacc(kint-1)+tiny
c            write (*,*) xacc(kint),xacc(kint-1)
         endif
      enddo
c it could happen that the change above yielded xacc() values greater
c than 1; should be fixed once more.
      xacc(ninter)=1d0
      do kint=1,ninter
         if (xacc(ninter-kint).gt.(xacc(ninter-kint+1)-tiny)) then
c            write (*,*) 'Accumulated results need adaptation #2:'
c            write (*,*) xacc(ninter-kint),xacc(ninter-kint+1),' become'
            xacc(ninter-kint)=1d0-dble(kint)*tiny
c            write (*,*) xacc(ninter-kint),xacc(ninter-kint+1)
         else
            exit
         endif
      enddo
cend RF

      if (plot_grid) then
         write(11,*) 'set limits x 0 1 y 0 1'
         write(11,*) 0, 0
         do kint=1,ninter
            write(11,*) xgrid(kint),xacc(kint)
         enddo
         write(11,*) 'join 1'
      endif

      do kint=1,ninter
         r=dble(kint)/dble(ninter)

         if (plot_grid) then
            write(11,*) 0, r
            write(11,*) 1, r
            write(11,*) ' join'
         endif

         do jint=1,ninter
            if(r.lt.xacc(jint)) then
               xn(kint)=xgrid(jint-1)+(r-xacc(jint-1))
     #        /(xacc(jint)-xacc(jint-1))*(xgrid(jint)-xgrid(jint-1))
               goto 11
            endif
         enddo
         if(jint.ne.ninter+1.and.kint.ne.ninter) then
            write(*,*) ' error',jint,ninter
            stop
         endif
         xn(ninter)=1
 11      continue
      enddo
      do kint=1,ninter
         xgrid(kint)=xn(kint)
         if (plot_grid) then
            write(11,*) xgrid(kint), 0
            write(11,*) xgrid(kint), 1
            write(11,*) ' join'
         endif
      enddo
      if (plot_grid) write(11,*) ' newplot'
      end

      subroutine nextlexi(ndim,iii,kkk,iret)
c kkk: array of integers 1 <= kkk(j) <= iii(j), j=1,ndim
c at each call iii is increased lexicographycally.
c for example, starting from ndim=3, kkk=(1,1,1), iii=(2,3,2)
c subsequent calls to nextlexi return
c         kkk(1)      kkk(2)      kkk(3)    iret
c 0 calls   1           1           1       0
c 1         1           1           2       0
c 2         1           2           1       0
c 3         1           2           2       0
c 4         1           3           1       0
c 5         1           3           2       0
c 6         2           1           1       0
c 7         2           1           2       0
c 8         2           2           1       0
c 9         2           2           2       0
c 10        2           3           1       0
c 11        2           3           2       0
c 12        2           3           2       1
      implicit none
      integer ndim,iret,kkk(ndim),iii(ndim)
      integer k
      k=ndim
 1    continue
      if(kkk(k).lt.iii(k)) then
         kkk(k)=kkk(k)+1
         iret=0
         return
      else
         kkk(k)=1
         k=k-1
         if(k.eq.0) then
            iret=1
            return
         endif
         goto 1
      endif
      end

      subroutine gen(fun,ndim,xgrid,ymax,ymax_virt,imode,x,vn)
c imode=0 to initialize
c imode=1 to generate
c imode=3 store generation efficiency in x(1)
      implicit none
      integer ndim,imode
      include "mint.inc"
      real * 8 fun,xgrid(0:nintervals,ndimmax),ymax(nintervals,ndimmax)
     $     ,ymax_virt,x(ndimmax)
      real * 8 dx(ndimmax),xx(ndimmax)
      integer icell(ndimmax),ncell(ndimmax),ncell_virt
      integer ifold(ndimmax),kfold(ndimmax)
      common/cifold/ifold
      real * 8 r,f(nintegrals),f1(nintegrals),ubound,vol,ran3
     $     ,xmmm(nintervals,ndimmax),dummy
      real * 8 rand(ndimmax)
      external fun,ran3
      integer icalls,mcalls,kdim,kint,nintcurr,nintcurr_virt,iret,ifirst
     $     ,i,vn,icalls_virt,mcalls_virt,icalls_nz,icalls_virt_nz
      double precision average_virtual,virtual_fraction
      common/c_avg_virt/average_virtual,virtual_fraction
      save icalls,mcalls,icalls_virt,mcalls_virt,xmmm,icalls_nz
     $     ,icalls_virt_nz
      if(imode.eq.0) then
         do kdim=1,ndim
            nintcurr=nintervals/ifold(kdim)
            xmmm(1,kdim)=ymax(1,kdim)
            do kint=2,nintcurr
               xmmm(kint,kdim)=xmmm(kint-1,kdim)+ymax(kint,kdim)
            enddo
            do kint=1,nintcurr
               xmmm(kint,kdim)=xmmm(kint,kdim)/xmmm(nintcurr,kdim)
            enddo
         enddo
         icalls=0
         icalls_nz=0
         mcalls=0
         icalls_virt=0
         icalls_virt_nz=0
         mcalls_virt=0
         xx(2)=0d0
         xx(3)=0d0
         xx(5)=0d0
         xx(6)=0d0
         return
      elseif(imode.eq.3) then
         if(icalls.gt.0) then
            x(1)=dble(mcalls)/dble(icalls)
            x(2)=xx(2)/dble(icalls)
            x(3)=xx(3)/dble(icalls)
         else
            x(1)=-1d0
            x(2)=-1d0
            x(3)=-1d0
         endif
         if(icalls_virt.gt.0) then
            x(4)=dble(mcalls_virt)/dble(icalls_virt)
            x(5)=xx(5)/dble(icalls_virt)
            x(6)=xx(6)/dble(icalls_virt)
         else
            x(4)=-1d0
            x(5)=-1d0
            x(6)=-1d0
         endif
         call increasecnt(' ',imode)
         return
      endif
      if (vn.eq.1) then
         mcalls_virt=mcalls_virt+1
      elseif(vn.eq.2 .or. vn.eq.3) then
         mcalls=mcalls+1
      else
         write (*,*) 'vn not correct in mint-integrator2.f',vn,imode
         stop
      endif
 10   continue
      if (vn.eq.1) then
         icalls_virt=icalls_virt+1
      elseif(vn.eq.2 .or. vn.eq.3) then
         icalls=icalls+1
      endif
      if (vn.eq.1) then
c Choose cell flat
         do kdim=1,ndim
            ncell(kdim)=min(int(ran3(.false.)*nintcurr)+1,nintcurr)
            rand(kdim)=ran3(.false.)
         enddo
      elseif(vn.eq.2 .or. vn.eq.3) then
         do kdim=1,ndim
            nintcurr=nintervals/ifold(kdim)
            r=ran3(.false.)
            do kint=1,nintcurr
               if(r.lt.xmmm(kint,kdim)) then
                  ncell(kdim)=kint
                  exit
               endif
            enddo
            rand(kdim)=ran3(.false.)
         enddo
      endif
      if (vn.eq.2 .or. vn.eq.3) then
         ubound=1
         do kdim=1,ndim
            ubound=ubound*ymax(ncell(kdim),kdim)
         enddo
      endif
      do kdim=1,ndim
         kfold(kdim)=1
      enddo
      do i=1,nintegrals
         f(i)=0
      enddo
      ifirst=0
 5    continue
      vol=1
      do kdim=1,ndim
         nintcurr=nintervals/ifold(kdim)
         nintcurr_virt=nintervals_virt/ifold(kdim)
         icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
         dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
         vol=vol*dx(kdim)*nintervals/ifold(kdim)
         x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
      enddo
      call get_ave_virt(x,nintcurr_virt,ndim,average_virtual)
      if (vn.eq.1) then
         ubound=ymax_virt
      endif
      dummy=fun(x,vol,ifirst,f1)
      do i=1,nintegrals
         f(i)=f(i)+f1(i)
      enddo
      ifirst=1
      call nextlexi(ndim,ifold,kfold,iret)
      if(iret.eq.0) goto 5
c get final value (x and vol not used in this call)
      dummy=fun(x,vol,2,f1)
      do i=1,nintegrals
         f(i)=f1(i)
      enddo
      if (vn.eq.2 .or. vn.eq.3) then
         xx(2)=xx(2)+f(2)
         xx(3)=xx(3)+f(1)
      else
         xx(5)=xx(5)+f(2)
         xx(6)=xx(6)+f(1)
      endif
      call increasecnt('another call to the function',imode)
      if (f(1).eq.0d0) then
         call increasecnt('failed generation cuts',imode)
      else
         if (vn.eq.1) then
            icalls_virt_nz=icalls_virt_nz+1
         elseif(vn.eq.2 .or.vn.eq.3) then
            icalls_nz=icalls_nz+1
         endif
      endif
      if(f(1).lt.0) then
         write(*,*) 'gen: non positive function'
         stop
      endif
      if(f(1).gt.ubound) then
         if (vn.eq.2) then
            call increasecnt('ubound fail novi',imode)
         elseif (vn.eq.1) then
            call increasecnt('ubound fail virt',imode)
         elseif (vn.eq.3) then
            call increasecnt('ubound fail born',imode)
         endif
      endif
      ubound=ubound*ran3(.false.)
      if(ubound.gt.f(1)) then
         call increasecnt
     &        ('vetoed calls in inclusive cross section',imode)
         goto 10
      endif
      if (vn.eq.2) then
         call increasecnt('events gen novi',imode)
      elseif (vn.eq.1) then
         call increasecnt('events gen virt',imode)
      elseif (vn.eq.3) then
         call increasecnt('events gen born',imode)
      endif
      end

c Dummy subroutine (normally used with vegas when resuming plots)
      subroutine resume()
      end

      subroutine increasecnt(argument,imode)
c Be careful, argument should be at least 15 characters
c long for this subroutine to work properly
      implicit none
      character*(*) argument
      character*15 list(100)
      integer ilist(0:100),i,j,imode
      logical firsttime
      data firsttime/.true./
      save ilist,list

      if (firsttime) then
         ilist(0)=1
         do i=1,100
            ilist(i)=0
            list(i)='               '
         enddo
         firsttime=.false.
      endif

      if(imode.ne.3) then
         i=1
         do while (i.le.ilist(0))
            if(i.eq.ilist(0)) then
               list(i)=argument(1:15)
               ilist(i)=1
               ilist(0)=ilist(0)+1
               goto 14
            endif
            if (argument(1:15).eq.list(i)) then
               ilist(i)=ilist(i)+1
               goto 14
            endif
            i=i+1
            if (i.ge.100) then
               write (*,*) 'error #1 in increasecnt'
               do j=1,ilist(0)
                  write (*,*) list(j),ilist(j)
               enddo
               stop
            endif
         enddo
 14      continue
      else
         do i=1,ilist(0)-1
            write (*,*) list(i),ilist(i)
         enddo
      endif
      end

      double precision function ran3(even)
      implicit none
      double precision ran2,get_ran
      logical even
      external get_ran
      if (even) then
         ran3=get_ran()
      else
         ran3=ran2()
      endif
      return
      end

      subroutine initialize_even_random_numbers(ncalls0,ndim,ncalls)
c Recompute the number of calls. Uses the algorithm from VEGAS
      implicit none
      integer ncalls0,ndim,ncalls,i
      integer dim,ng,npg,k
      common /even_ran/dim,ng,npg,k
c Number of dimension of the integral
      dim=ndim
c Number of elements in which we can split one dimension
      ng=(ncalls0/2.)**(1./ndim)
c Total number of hypercubes
      k=ng**ndim
c Number of PS points in each hypercube (at least 2)
      npg=max(ncalls0/k,2)
c Number of PS points for this iteration
      ncalls=npg*k
      return
      end

      double precision function get_ran()
      implicit none
      double precision ran2,dng
      external ran2
      logical firsttime
      data firsttime/.true./
      integer dim,ng,npg,k
      common /even_ran/dim,ng,npg,k
      integer maxdim
      parameter (maxdim=100)
      integer iii(maxdim),kkk(maxdim),i,iret
      integer current_dim
      save current_dim,dng,kkk,iii
      if (firsttime) then
         dng=1d0/dble(ng)
         current_dim=0
         do i=1,dim
           iii(i)=ng
           kkk(i)=1
        enddo
        firsttime=.false.
      endif
      current_dim=mod(current_dim,dim)+1
c This is the random number in the hypercube 'k' for current_dim
      get_ran=dng*(ran2()+dble(kkk(current_dim)-1))
c Got random numbers for all dimensions, update kkk() for the next call
      if (current_dim.eq.dim) then
         call nextlexi(dim,iii,kkk,iret)
         if (iret.eq.1) then
            call nextlexi(dim,iii,kkk,iret)
         endif
      endif
      return
      end

      subroutine init_ave_virt(ninter,ndim)
      implicit none
      include "mint.inc"
      integer kdim,ndim,ninter,i
      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
     $     ,ndimmax)
      double precision ave_virt(nintervals_virt,ndimmax)
     $     ,ave_virt_acc(nintervals_virt,ndimmax)
     $     ,ave_born_acc(nintervals_virt ,ndimmax)
      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
     $     ,nvirt_acc
      do kdim=1,ndim
         do i=1,ninter
            nvirt(i,kdim)=0
            ave_virt(i,kdim)=0d0
            nvirt_acc(i,kdim)=0
            ave_virt_acc(i,kdim)=0d0
            ave_born_acc(i,kdim)=0d0
         enddo
      enddo
      return
      end

      subroutine get_ave_virt(x,ninter,ndim,average_virtual)
      implicit none
      include "mint.inc"
      integer kdim,ndim,ninter,ncell
      double precision x(ndimmax),average_virtual
      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
     $     ,ndimmax)
      double precision ave_virt(nintervals_virt,ndimmax)
     $     ,ave_virt_acc(nintervals_virt,ndimmax)
     $     ,ave_born_acc(nintervals_virt ,ndimmax)
      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
     $     ,nvirt_acc
      average_virtual=0d0
      do kdim=1,ndim
         ncell=min(int(x(kdim)*ninter)+1,ninter)
         average_virtual=average_virtual+ave_virt(ncell,kdim)
      enddo
      average_virtual=average_virtual/ndim
      return
      end

      subroutine fill_ave_virt(x,ninter,ndim,virtual,born)
      implicit none
      include "mint.inc"
      integer kdim,ndim,ninter,ncell
      double precision x(ndimmax),virtual,born
      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
     $     ,ndimmax)
      double precision ave_virt(nintervals_virt,ndimmax)
     $     ,ave_virt_acc(nintervals_virt,ndimmax)
     $     ,ave_born_acc(nintervals_virt ,ndimmax)
      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
     $     ,nvirt_acc
      do kdim=1,ndim
         ncell=min(int(x(kdim)*ninter)+1,ninter)
         nvirt_acc(ncell,kdim)=nvirt_acc(ncell,kdim)+1
         ave_virt_acc(ncell,kdim)=ave_virt_acc(ncell,kdim)+virtual
         ave_born_acc(ncell,kdim)=ave_born_acc(ncell,kdim)+born
      enddo
      return
      end

      subroutine regrid_ave_virt(ninter,ndim)
      implicit none
      include "mint.inc"
      integer ninter,ndim,kdim,i
      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
     $     ,ndimmax)
      double precision ave_virt(nintervals_virt,ndimmax)
     $     ,ave_virt_acc(nintervals_virt,ndimmax)
     $     ,ave_born_acc(nintervals_virt ,ndimmax)
      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
     $     ,nvirt_acc
c need to solve for k_new = (virt+k_old*born)/born
      do kdim=1,ndim
         do i=1,ninter
            if (ave_born_acc(i,kdim).eq.0d0) cycle
            if (ave_virt(i,kdim).eq.0d0) then ! i.e. first iteration
               ave_virt(i,kdim)= ave_virt_acc(i,kdim)/ave_born_acc(i
     $              ,kdim)+ave_virt(i,kdim)
            else  ! give some importance to the iterations already done
               ave_virt(i,kdim)=(ave_virt_acc(i,kdim)/ave_born_acc(i
     $              ,kdim)+ave_virt(i,kdim)*2d0)/2d0
            endif
         enddo
      enddo
c reset the acc values
      do kdim=1,ndim
         do i=1,ninter
            nvirt(i,kdim)=nvirt(i,kdim)+nvirt_acc(i,kdim)
            nvirt_acc(i,kdim)=0
            ave_born_acc(i,kdim)=0d0
            ave_virt_acc(i,kdim)=0d0
         enddo
      enddo
      return
      end

      subroutine double_ave_virt(ninter,ndim)
      implicit none
      include "mint.inc"
      integer kdim,ndim,i,ninter
      integer nvirt(nintervals_virt,ndimmax),nvirt_acc(nintervals_virt
     $     ,ndimmax)
      double precision ave_virt(nintervals_virt,ndimmax)
     $     ,ave_virt_acc(nintervals_virt,ndimmax)
     $     ,ave_born_acc(nintervals_virt ,ndimmax)
      common/c_ave_virt/ave_virt,ave_virt_acc,ave_born_acc,nvirt
     $     ,nvirt_acc
      do kdim=1,ndim
         do i=ninter,1,-1
            ave_virt(i*2,kdim)=ave_virt(i,kdim)
            if (nvirt(i,kdim).ne.0) then
               nvirt(i*2,kdim)=max(nvirt(i,kdim)/2,1)
            else
               nvirt(i*2,kdim)=0
            endif
            if (i.ne.1) then
               ave_virt(i*2-1,kdim)=(ave_virt(i,kdim)
     $              +ave_virt(i-1,kdim))/2d0
               if (nvirt(i,kdim)+nvirt(i-1,kdim).ne.0) then
                  nvirt(i*2-1,kdim)=
     &                 max((nvirt(i,kdim)+nvirt(i-1,kdim))/4,1)
               else
                  nvirt(i*2-1,kdim)=0
               endif
            endif
         enddo
      enddo
      return
      end


text:
Address(Willow Lane) is a residential street in Central Davis.

 Intersecting Streets

  Radcliffe Drive
  Fordham Drive and Lemon Lane
  Fordham Drive again
  Villanova Drive

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


MODULE hamiltonian_matrix_module

  use hsort_module
  use kinetic_variables, only: Md
  use kinetic_sol_module
  use construct_matrix_ps_nloc2_module
  use localpot_module
  use symmetry_module
  use rgrid_module
  use bz_module, only: kbb

  implicit none

  PRIVATE
  PUBLIC :: construct_hamiltonian_matrix

#ifdef _DRSDFT_
  real(8),parameter :: zero=0.0d0
  real(8),allocatable :: Hmat(:,:)
  real(8),allocatable :: w1(:,:), w2(:,:)
#else
  complex(8),parameter :: zero=(0.0d0,0.0d0)
  complex(8),allocatable :: Hmat(:,:)
  complex(8),allocatable :: w1(:,:), w2(:,:)
#endif

  integer :: ML

  real(8),allocatable :: SymMat(:,:)

CONTAINS

  SUBROUTINE construct_hamiltonian_matrix( ML_in )
    implicit none
    integer,intent(IN) :: ML_in
    integer :: i,n,k,s,isym,i1,i2,i3,j,i0,j1,j2,j3
    integer :: k1,k2,k3,nkg
    real(8) :: c,r0(3),r1(3),r,tmp(3),err
    integer,allocatable :: sym_mat(:,:,:)

    write(*,'(a50," construct_hamiltonian_matrix")') repeat("-",50)

    ML = ML_in

    c = 1.0d0/nsym

    allocate( Hmat(ML,ML)   ) ; Hmat=zero
    allocate( SymMat(ML,ML) ) ; SymMat=0.0d0
    allocate( w1(ML,ML)     ) ; w1=zero
    allocate( w2(ML,ML)     ) ; w2=zero

    do s=1,1
    do k=1,1

       Hmat(:,:) = zero

       call construct_matrix_kinetic_sol( k, ML, Hmat )
!       call construct_matrix_localpot( s, ML, Hmat )
!       call construct_matrix_ps_nloc2( k, ML, Hmat )

       n=count(abs(Hmat)>1.d-10)
       write(*,*) n,dble(n)/dble(ML)
!       do i=1,ML
!          write(*,'(1x,2i5,2f20.15)') i,count(abs(Hmat(i,:))>1.d-10) &
!          ,Hmat(i,i)
!       end do
       !call check_hamil
!       i1=Ngrid(1)/2
!       i2=Ngrid(2)/2
!       i3=Ngrid(3)/2
!       i0=1+i1+i2*Ngrid(1)+i3*Ngrid(1)*Ngrid(2)
!       write(*,*) i0,i1,i2,i3
!       j=0
!       do i3=0,Ngrid(3)-1
!       do i2=0,Ngrid(2)-1
!       do i1=0,Ngrid(1)-1
!          i=1+i1+i2*Ngrid(1)+i3*Ngrid(1)*Ngrid(2)
!          if ( abs(Hmat(i0,i)) > 1.d-10 ) then
!             j=j+1
!             w1(1,j)=i1
!             w1(2,j)=i2
!             w1(3,j)=i3
!             w1(4,j)=Hmat(i0,i)
!             if ( i == i0 ) then
!                r0(1)=i1
!                r0(2)=i2
!                r0(3)=i3
!             end if
!          end if
!       end do
!       end do
!       end do
!       do i=1,j
!          r1(1:3)=real(w1(1:3,i))
!          r=sqrt(sum((r1-r0)**2))
!          write(*,'(1x,4i5,2f20.15)') i,nint(real(w1(1:3,i))),real(w1(4,i)),r
!       end do

       if ( isymmetry /= 0 ) then
          call get_mat_symmetry( sym_mat )
          w2(:,:)=zero
          nkg=0
          do isym=1,nsym
             tmp(:)=matmul(kbb(:,k),sym_mat(:,:,isym))
             !write(*,'(1x,i4,3f8.3)') isym,tmp
             err=abs(sum((tmp-kbb(:,k))**2))
             if ( err > 1.d-8 ) cycle
             nkg=nkg+1
             write(*,'(1x,"isym/nsym, nkg=",i3," /",i3,i4)') isym,nsym,nkg
             call construct_matrix_symmetry( isym, ML, SymMat )
             w1(:,:) = matmul( Hmat, transpose(SymMat) )
             w2(:,:) = w2(:,:) + matmul( SymMat, w1 )
          end do
          c=1.0d0/nkg
          Hmat(:,:) = c*w2(:,:)
       end if
       n=count(abs(Hmat)>1.d-10)
       write(*,*) n,dble(n)/dble(ML)
!       do i=1,ML
!          write(*,'(1x,2i5,2f20.15)') i,count(abs(Hmat(i,:))>1.d-10) &
!          ,Hmat(i,i)
!       end do
       !call check_hamil(10)
!       j=0
!       do i3=0,Ngrid(3)-1
!       do i2=0,Ngrid(2)-1
!       do i1=0,Ngrid(1)-1
!          i=1+i1+i2*Ngrid(1)+i3*Ngrid(1)*Ngrid(2)
!          if ( abs(Hmat(i0,i)) > 1.d-10 ) then
!             j=j+1
!             w1(1,j)=i1
!             w1(2,j)=i2
!             w1(3,j)=i3
!             w1(4,j)=Hmat(i0,i)
!          end if
!       end do
!       end do
!       end do
!       do i=1,j
!          r1(1:3)=real(w1(1:3,i))
!          r=sqrt(sum((r1-r0)**2))
!          write(*,'(1x,4i5,2f20.15)') i,nint(real(w1(1:3,i))),real(w1(4,i)),r
!       end do

!       call construct_matrix_kinetic_sol( k, ML, Hmat )
       call construct_matrix_localpot( s, ML, Hmat )
       call construct_matrix_ps_nloc2( k, ML, Hmat )

       call diag_Hmat

    end do ! k
    end do ! s

    deallocate( w2 )
    deallocate( w1 )
    deallocate( SymMat )
    deallocate( Hmat )

  END SUBROUTINE construct_hamiltonian_matrix

  SUBROUTINE diag_Hmat
    implicit none
    integer :: LWORK,LRWORK,ierr,i
    real(8),allocatable :: rwork(:),eigval(:)
    complex(8),allocatable :: work(:)
    allocate( eigval(ML) )
#ifdef _DRSDFT_
    LWORK=3*ML-1
    allocate( rwork(LWORK) )
    call dsyev('N','L',ML,Hmat,ML,eigval,rwork,LWORK,rwork,ierr)
    deallocate( rwork )
#else
    LWORK=2*ML-1
    LRWORK=3*ML-2
    allocate( work(LWORK),rwork(LRWORK) )
    call zheev('N','L',ML,Hmat,ML,eigval,work,LWORK,rwork,ierr)
    deallocate( rwork,work )
#endif
    do i=1,10
       write(*,'(1x,i5,f20.15)') i,eigval(i)
    end do
    do i=1,3
       write(*,'(1x,5x,".")')
    end do
    do i=ML-9,ML
       write(*,'(1x,i5,f20.15)') i,eigval(i)
    end do
    deallocate( eigval )
  END SUBROUTINE diag_Hmat

  SUBROUTINE check_hamil( unit )
    implicit none
    integer,optional,intent(IN) :: unit
    integer :: i1,i2,i3,j1,j2,j3,k1,k2,k3,i,j,i0
    integer,allocatable :: icheck(:,:,:),itmp(:,:),jtmp(:)
    real(8),allocatable :: rcheck(:,:,:),rtmp(:)
    if ( present(unit) ) rewind unit
    allocate( icheck(-2*Md:2*Md,-2*Md:2*Md,-2*Md:2*Md) ) ; icheck=0
    allocate( rcheck(-2*Md:2*Md,-2*Md:2*Md,-2*Md:2*Md) ) ; rcheck=0.0d0
    do i3=0,Ngrid(3)-1
    do i2=0,Ngrid(2)-1
    do i1=0,Ngrid(1)-1
       i=1+i1+i2*Ngrid(1)+i3*Ngrid(1)*Ngrid(2)
       if ( i /= 1 ) cycle
       write(*,*) "grid",i1,i2,i3,i
       i0=0
       do j3=0,Ngrid(3)-1
       do j2=0,Ngrid(2)-1
       do j1=0,Ngrid(1)-1
          j=1+j1+j2*Ngrid(1)+j3*Ngrid(1)*Ngrid(2)
          if ( abs(Hmat(i,j))>1.d-10 ) then
             i0=i0+1
             k1=j1-i1 ; if ( k1 >= Ngrid(1)/2 ) k1=k1-Ngrid(1)
             k2=j2-i2 ; if ( k2 >= Ngrid(2)/2 ) k2=k2-Ngrid(2)
             k3=j3-i3 ; if ( k3 >= Ngrid(3)/2 ) k3=k3-Ngrid(3)
             write(*,'(1x,7i4,2x,2f20.15)') i0,k1,k2,k3,j1-i1,j2-i2,j3-i3 &
                  ,Hmat(i,j),Hmat(i,j)*Hgrid(1)**2
!             if ( present(unit) ) then
!                write(unit,'(1x,3i4,f22.16)') k1,k2,k3,Hmat(i,j)*Hgrid(1)**2
!             end if
             icheck(k1,k2,k3)=icheck(k1,k2,k3)+1
             rcheck(k1,k2,k3)=Hmat(i,j)*Hgrid(1)**2
          end if
       end do
       end do
       end do
    end do
    end do
    end do
    write(*,*) maxval(icheck),sum(icheck),count(icheck/=0)
    i=count(icheck/=0)

    if ( present(unit) ) write(unit,*) i

    allocate( itmp(3,i), jtmp(i), rtmp(i) )
    itmp=0 ; jtmp=0 ; rtmp=0.0d0
    i=0
    do i3=-2*Md,2*Md
    do i2=-2*Md,2*Md
    do i1=-2*Md,2*Md
       if ( icheck(i1,i2,i3) /= 0 ) then
          i=i+1
          itmp(1,i)=i1
          itmp(2,i)=i2
          itmp(3,i)=i3
          rtmp(i)=i1*i1+i2*i2+i3*i3
          if ( present(unit) ) then
             write(unit,'(1x,3i4,f22.16)') i1,i2,i3,rcheck(i1,i2,i3)
          end if
       end if
    end do
    end do
    end do
    call indexx( size(rtmp), rtmp, jtmp )
    do i=1,size(rtmp)
       j=jtmp(i)
       i1=itmp(1,j)
       i2=itmp(2,j)
       i3=itmp(3,j)
       write(*,'(1x,4i4,2f20.10)') i, itmp(:,j), rtmp(j), rcheck(i1,i2,i3)
    end do
    deallocate( itmp, jtmp, rtmp )
    deallocate( rcheck )
    deallocate( icheck )
  END SUBROUTINE check_hamil

END MODULE hamiltonian_matrix_module


text:
>“We give you peace of mind about the safety & reliability of your car”<

Center City Automotive owned by Users/d&mHatton Doug & Margee Hatton, is a quality, leading automotive repair center providing competent and honest service since 1981. Friendly atmosphere and topnotch advice at the front counter and on the phones is provided by Doug Hatton & Jim Wolgemuth. Automotive technicians Dave Haverstock, Steve Moore, Shawn Heaton, Rick Mortensen and Sean Sears dont just fix or maintain your car, they get to know you and your car. Let Center City Automotive keep your vehicle running right, and when needed, help you get back on the road with the right repair solution.

See Center City Automotive on Facebook:  http://www.facebook.com/centercityautomotive

  See also: Automobile Repair.

20060901 09:47:33 nbsp Ive been using them for a couple years and I really like them. They are honest, do good work, and every once in a while, throw in an extra fix that you didnt ask for. After going to Advanced Auto for a while, Ive been very happy with the service at Center City. GreenThing

20070601 18:19:56 nbsp Great service! The guys there were extremely helpful, thorough, and very quick.  Users/LinneaUggla

20070802 17:38:41 nbsp Today was my first experience with Center City and I will definitely be trusting them with my car for future repairs.  They were professional, made the necessary repairs quickly (same day) and for much less money than I expected to pay!!  They also took the initiative to do a quick scan of other systems, including belts to ensure all was fine.  They have certainly earned my repeat business! Users/scd

20070827 20:04:31 nbsp This is a great place to take your car to when having troubles.  They are reliable, knowledgeable and honest...which is everything you need when looking for someone to fix your car.  They do let you know what youre paying for and why its something you need to fix.  Definitely the go to place in town! Users/Danya

20070904 14:23:07 nbsp This place is great. My Toyota truck had a blown headgasket, and when my dad and I decided to repair it ourselves instead of having them do it (the cost of repairs would have been just short of the value of the car) they gave us plenty of advice and tips. Ive had them do lots of other service, and wont ever go anywhere else. Users/BryceEsch

20071110 23:53:49 nbsp Very nice people. Friendly, prompt and honest. They explained to me what needed to be change and why. I went in originally to fix a slipping clutch, they installed the new clutch; however; the new clutch was defective. The guys at Center City reordered a new and more expensive replacement clutch for me at no extra cost. After my first experience, I went back and had my timing belt and water pump fixed. Would definitely recommend them. Users/GMC

20080211 10:51:38 nbsp This is bar none the best independent shop in Davis.  I have been going there for years, and before that my family.  They are honest and reliable, I have NEVER been up sold or told something that turns out not to be the case.  their price is always competitive to big chains, and they know your name and your car.  they will give you the truth even if it means they lose a job.  When a repair doesnt solve the issue as well as you would like they are eager to find out what needs to be done and always fair with price.  ALL shops will make mistakes, faulty parts happen etc... it is HOW they deal with it that counts!

example:  My front end was squeaking, they said it was a dry joint someplace in the suspension...  said it is an ungreasable joint, and is replace or nothing.  BUT, they then said that they tested it for wear, and it is still well within tolerances, so as long as I can deal with the squeek, it is no biggie!  I was shocked, probably $500 in labor, any other shop would have just said, yeah, your flux capacitor is out, you need a new one.

Doug and Jim at the desk are a bit bearded up as one commenter mentioned... but they are honest and listen to you...  much better than some 21 year old clown at the desk of speedylube or whatever.  It is a great feeling to have a trusted mechanic... to know that although a repair might not be cheap, they will do what needs to be done right, and at the best possible price.  Users/Carsonsig

20080514 01:03:42 nbsp Excellent service! I took my Nissan car to them to check the Service Engine Soon light. They found out that it was my cars front catalytic converters problem. They called the Nissan dealer and found out that my car was still under warranty (I bought it as a used car, so myself did not realize Nissans 8 year/80000 miles warranty). And the warranty was going to expire in 3 weeks! It was so nice for them to find out that I was still under coverage. Then I called Nissan and got my part replaced for free (or I have to pay almost $500 for a cat converter!).

So I strongly recommend this mechanics store! People are honest, knowledge and professional. Whats more important, I feel that they were really standing on my side, gave good suggestions to save $$ and get the best deal. Users/ezhang

20080710 15:31:20 nbsp they are the best. i refuse to take my car anywhere else in davis. my car is there right now and ive been going to them for like 2 years or 3. they dont make me feel like an idiot when i tell them whats wrong w/ my car even though i should because sometimes its me running the battery out on accident. today i was pushed into the schedule even though they were overbooked w/ cars. they said theyd try to fix it in the best time they could. they are really honest there and grateful for the trust that people have for them. also they remember you well.  Users/KimberlyMiller

20081110 23:29:35 nbsp Overall experience: Very Good, Very Satisfy, Very Professional.
I needed a replacement of a Starter in my car, and I called around Davis to get the best estimation of the cost. Center City Autos called back with the best price, and agree to use factory manufactured parts for my car, which is exactly what I wanted. I brought my car in around 12pm on Monday, and Jim/Doug called me that my car is ready to go at 4pm of that same day.
Very fast services, and Very friendly. Jim/Doug would chat with you and answer questions relating to car that you may have. They do look like two hippies (especially when they said Peace to me, Hahaha). My overall experience here is very good because of their friendly manner, and their pricing were the cheapest that I could find in Davis. I recommend this place to whoever need their car diagnose and/or fixing.  Users/crepelover

20081112 15:44:18 nbsp This place is FANTASTIC. My windshield wipers stopped working on the first truly rainy day this year. I dropped one child at preschool, called Center City on my way home, was able to drop the car off immediately, and my infant daughter and I were given a ride home by courteous staff member. Two hours later they called to let me know the wipers were fixed (wiring was loose) and that I could use a new battery...which they installed free of charge (only paid cost of the battery). They picked me up within 1/2 hour and I made it back to daycare pickup on time. Amazing service. Far exceeded expectations. Users/JMakus

20090304 11:14:36 nbsp honest, great, fast. not too expensive. I have never been to a better auto shop! Users/yadayada

20090626 21:26:34 nbsp I took my truck to center city auto because one of the spark plugs blew out. I took it in knowing that this was the only place in davis that could work on it so I was skeptical. They told me that the spark plug hole was as big as it got and I needed a new engine that would cost me anywhere between 3,000 to 6,0000 dollars to fix. So I got it towed to San Francisco to a family frieds shop and the friend told me that Center City was lying to me and he fixed it for 200 dollars including service fees. Users/kseabright

20090713 09:01:39 nbsp Took my car in for a couple of very simple repairs and was very pleased. The repairs were done quickly, well, and for a reasonable price.  The guy at the counter was really nice.  Ill definitely go back if I have more serious problems. Users/DataPoint

20091016 20:48:31 nbsp Speedy service, but extremely irritated that they charged me for something they didnt repair. Now I have to wait until Monday to take the car in  again. Im sure Ill either get a refund or proper repair on Monday, but its extremely frustrating that they dont have an answering machine of some sort. Users/ElizabethAbinante

20100604 08:08:53 nbsp Jim and Doug have done very good jobs on my cars, always with unusual friendliness. I feel so much reassured now with a trustworthy auto repair place to have. More Davis business folks could try to imitate them!Subodh Users/SubodhJain

20101011 13:09:29 nbsp These guys are EXCELLENT. After my A/C was misdiagnosed 4 different times by another auto body shop, I took it to Center City and they got it right the first time. Ever since then my A/C has been working perfectly. Everyone here seems to be super friendly and knowledgeable as well. I plan on taking my car here and only here in the future.    Users/trevorchreid

20110218 12:03:11 nbsp These guys are the best! Fantastic service, very honest and reasonable. We been going here for many years. They always go above and beyond. Thanks Doug and Jim and your great team!! Users/agerould

20110513 15:49:02 nbsp I have gone to Center City for a few years and had nothing but positive experiences. My car broke down a month after I moved to Davis and I knew no one in town, but a kind couple passing by used their AAA membership to tow my car to their shop (Center City) which they highly recommended. I was not disappointed!! My car was repaired the next day and for a very reasonable price. Every time my car has had a problem I have taken it to them and they have been just awesome and always treated me kindly. I was driving an old car (25+ years) and these old cars have problems, you know, but I dont think I ever saw them twice for the same problem. Once they fixed something it was FIXED. I could not ask for better repair service, really. They also back up their work, so if you have questions, concerns, or are just confused about how your newly repaired car is running they are super patient. Yeah, Doug and Jim kind of look like hippies and have some of that hippie zen to them, but they are undoubtably solid and professional mechanics. I ended up selling my old car and they were the first people I called when something was wrong with my new one. Unfortunately they dont work on the type of new car I got, which made me sad since I trust them and they have done good work for me in the past, but they gave me a recommendation to another shop that is on the same level as them in terms of quality, honesty, and value.  Users/DanaMace

20110815 20:03:14 nbsp I took my car here for an oil change. It wasnt the cheapest oil change (It would have been $50 if I hadnt brought my own oil and filter, so it would be nice if they had coupons or a student discount.), but they were really helpful. After they saw that my oil level was a lot lower than it should have been, they checked around, found some potential problems, and patiently explained them to me without making me feel pressured to spend hundreds on a bunch of repairs. They also washed my windshield and topped off my power steering fluid for free. This is the first time the oil change wasnt done by my dad so I have no idea if thats customary, but either way I thought it was a nice touch. Update: As a quick note about my oil change, I noticed later that the oil was a whole half inch above the full line on the dipstick. Not sure how/why they put in too much, but they took the extra out for me.Users/MeggoWaffle

   The big reason Im editing this comment is I recommended this place to my partner and he had a very mixed experience. I have pasted his account below:

My experience with Center City was mixed.  I brought my car in when a belt broke, thinking it would be a simple fix.  However it turned out several things compounded to cause the belt to break and it would be more money to fix than the car was really worth to me.  So I asked them if they knew if any mechanic would be interested in buying it off my hands and having a fixer upper.  The person I spoke with kindly said he would check around.  Unfortunately nobody was interested, so I came in to ask what my options were.  He came up with a third option: they fix as little as possible to get the car running so I could retire my car with the state of California.  They bought cheap parts and fixed as little as necessary for the car to make the drive to Sacramento to the salvage yard.  This saved me a lot of money.  I really appreciated how they worked with me for a solution.  The negative was their communication.  They constantly left me waiting to hear from them, even when they said they were going to call.  Looking at my phone records, they only called me twice during this whole ordeal.  I was forced to call them when I hadnt heard from them in days.  The worst was that they kept pushing back the day when they said the car would be ready, and it ended up taking a week to fix my car up just so it would be drivable and I could sell it to the state.  They were clearly putting other cars ahead of mine, which was frustrating considering I was a paying customer like any other and it had already been there for nearly a month.  For the $90 an hour they charge for their service, it didnt feel even remotely worth it.  Also they misquoted the original cost to fix my car by $200 due to bad math.  So the positive is that they worked with me to find an alternate solution to my car troubles which was nice.  They were also extremely friendly when I spoke with them.  The negatives were that they seemed to put my car on the back burner while they worked with other cars and they were very bad at communicating with me.  They also charge $90 an hour for labor which is in the upper end of what mechanics normally charge.

20110825 17:26:37 nbsp Ive gotta agree with most everyone here  theyve got the best prices around, they do great work, and they dont just do the bare minimum  they genuinely want your car to run as well as possible. Theyre honest and are great at explaining everything having to do with your car. Ive tried a number of other places before, and now wont go anywhere else. You really do feel like theyre on your side  not trying to cheat you out of money or say that your car needs stuff it doesnt. Im glad theyre in town. Plus, how often do you get your car serviced by the lookalike members of ZZ Top? ;)  Users/mday24

20120531 15:47:40 nbsp Today we noticed that our headlight was out on our Prius.  The last time the headlight went out, I brought the car to the Toyota dealer in Davis and they charged me $150 to change the bulb claiiming they had to remove the bumper.  I called them today and they gave me the same price.  I called my local mechanic, Center City, and they told me to bring the car over and they would take a look.  They changed the bulb after working on it for half an hour, charged me $15 for the bulb (the bulb at Toyota was $40) and sent me on my way.  I love Center City, have been taking my cars there for repairs for years and have always been more than satisfied.  I recommend them highly for any auto repairs.
Shelly Bailes/Ellen Pontac Users/ShellyBailes

20130404 13:31:43 nbsp Ive been going to Center City for a few years now and like them a lot.  Ive always felt theyve given me good work at a decent price.  Ive never felt pressured to get extra work done, (in fact, I think they go out of thier way to not do this).  The last time I brought my car in for major work, they replaced my radiator and water pump.  The next week I noticed some minor coolant leakage, and when I called back they repaired it quickly and free of charge.  It was the type of service I would expect from a good mechanic.  Also, while a lot of garages can look somewhat sloppy, the Center City facily is always clean and organized.  Doesnt necessarily mean they are altogether better than anyone else, but think it suggestive of the type of work they do, i.e. professional and honest.  Highly recommended. Users/jefftolentino

20140327 21:07:04 nbsp Knowledgeable, fair, honest service.  I have been very pleased with both major and minor repairs by Center City. Users/Nokomis


text:
Describe Users/ILikeToEat here.

Welcome to the Wiki please consider the importance of using your RealName  Users/StevenDaubert


text:
      INTEGER FUNCTION BUNPK (IG,I,J)
C
C     THIS ROUTINE IS USED ONLY IN BANDIT MODULE
C
C     UNPACK INTEGER GRID NO. FROM IG TABLE.   SEE BPACK FOR PACKING
C     USE APPROP. PORTION OF THIS ROUTINE FOR DIFFERENT TYPE OF MACHINE.
C
      INTEGER*2        IG(1)
C
      COMMON /BANDB /  NBIT,     DUM3B(3), IPASS,    NW,       DUM1B,
     1                 NBPW
      COMMON /BANDS /  DUM4S(4), II1,      DUM5S(5), MASK
C
      IPASS=IPASS+1
      LOC  =J-1
C
C     ********************************************
C     UNIVAC AND CDC MACHINES
C     ********************************************
C     INTEGER          RSHIFT,   ANDF
C
C     N1 =II1*(LOC/NW)+I
C     N2 =MOD(LOC,NW)*NBIT+NBIT
C     LOC=RSHIFT(IG(N1),NBPW-N2)
C     BUNPK=ANDF(LOC,MASK)
C     RETURN
C
C     ********************************************
C     IBM AND VAX MACHINES
C     (IG IS SET TO INTEGER*2 IN BPACK AND BUNPK, ELSEWHERE INTEGER*4)
C     INTEGER*2        IG(1)
C     ********************************************
C
      N1=II1*LOC+I
      BUNPK=IG(N1)
      RETURN
      END


text:
      SUBROUTINE LBCOLS(NCOLRS,LSTLEN)
C
C          ------------------------------------------------
C          ROUTINE NO. ( 107)   VERSION (A8.1)    12:MAR:86
C          ------------------------------------------------
C
C          THIS DEFINES THE CURRENT LABELS COLOUR LIST.
C
C
C          THE ARGUMENTS ARE AS FOLLOWS:
C
C          [NCOLRS] IS AN ARRAY OF COLOUR NUMBERS,
C          <LSTLEN> IS THE SIZE OF THE ARRAY.
C
C
      INTEGER NCOLRS(100)
C
      COMMON /T0KLAB/ LABCL0(100),LENLAB
      COMMON /T0TRAC/ IPRINT
C
C
      CALL G3INIT(2)
C
      IF (IPRINT.EQ.1) CALL G0MESG(164,0)
C
      LENLAB= LSTLEN
      IF (LENLAB.LE.0)   RETURN
      IF (LENLAB.GT.100) LENLAB= 100
C
      DO 100 ISAVE= 1,LENLAB
        LABCL0(ISAVE)= NCOLRS(ISAVE)
        IF (NCOLRS(ISAVE).LT.1.OR.NCOLRS(ISAVE).GT.255) LABCL0(ISAVE)= 0
  100 CONTINUE
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


program pfasst_imk
  use pf_mod_mpi
  integer :: rank, err

  call mpi_init(err)
  if (err .ne. 0) stop "ERROR: Can't initialize MPI."
  call mpi_comm_rank(MPI_COMM_WORLD, rank, err)

  call simulate()
  if(rank == 0) print *,'shutting down mpi'

  call mpi_finalize(err)
  if(rank == 0) print *,'all done - bye bye'

contains
  subroutine simulate()
      use pfasst     !< library module defining highest-level pfasst object
      use pf_mod_mpi !< library module for mpi-related business

      use probin     !< should be library module for reading/parsing problem parameters

      use mod_zmkpair    !< prog-specified module containing solution type information
      use sweeper    !< prog-specified module containing sweeper information
      use pf_my_level    !< prog-specified module level
      use hooks      !< prog-specified module containing program hooks
      use utils

      implicit none

      type(pf_pfasst_t)  :: pf          !< pfasst data structure
      type(pf_comm_t)    :: comm        !< mpi communicator

      type(zmkpair) :: q0, qend

      character(256) :: probin_fname       !<  file name for input
      integer    :: err, l,mpibuflen
      real(pfdp) :: start, finish

      probin_fname = "probin.nml"
      if (command_argument_count() >= 1) &
         call get_command_argument(1, value=probin_fname)
      call probin_init(probin_fname)

      call pf_mpi_create(comm, MPI_COMM_WORLD)
      call pf_pfasst_create(pf, comm, fname=probin_fname)

      !---- Create the levels ---------------------------------------------
      do l = 1, pf%nlevels
         !>  Allocate the user specific level object
         allocate(imk_context::pf%levels(l)%ulevel)

         !>  Allocate the user specific data factory
         allocate(zmkpair_factory::pf%levels(l)%ulevel%factory)

         !>  Add the sweeper to the level
         allocate(imk_sweeper_t::pf%levels(l)%ulevel%sweeper)

         mpibuflen = nparticles * nparticles * 2
         call pf_level_set_size(pf,l,[nparticles],mpibuflen)
      end do

      if(pf%rank == 0) print *,'Initializing mpi and pfasst...'
      call pf_pfasst_setup(pf)
      call pf_add_hook(pf, -1, PF_POST_ITERATION, echo_residual)
      call pf_add_hook(pf, -1, PF_POST_CONVERGENCE, echo_error)
      !      if (save_solutions) call pf_add_hook(pf, -1, PF_POST_ITERATION, save_solution)
      if (save_solutions) call pf_add_hook(pf, -1, PF_POST_CONVERGENCE, save_solution)

      !>  output the run options
      call pf_print_options(pf,un_opt=6)

      call zmkpair_build(q0, pf%levels(pf%nlevels)%lev_shape(1))
      call zmkpair_build(qend, pf%levels(pf%nlevels)%lev_shape(1))
      call initial(q0)

      call mpi_barrier(MPI_COMM_WORLD, err)

      start = MPI_Wtime()
      if(pf%rank == 0) print*, 'Running pfasst...'
      call pf_pfasst_run(pf, q0, dt, 0.0_pfdp, nsteps, qend)

      finish = MPI_Wtime()
      print '("processor ", i2, " returned from pfasst CPU time: ", f16.10, " seconds")', &
           pf%rank, finish-start

      call mpi_barrier(MPI_COMM_WORLD, err)

      if(pf%rank == comm%nproc-1) then
        call qend%write_to_disk('final_solution') !necessary for pfasst.py
        if (pf%debug) call qend%eprint() !only for debug purpose
      call qend%eprint() !only for debug purpose
      endif

      call zmkpair_destroy(q0)
      call zmkpair_destroy(qend)

      if(pf%rank == 0) print *,'destroying pf'
      call pf_pfasst_destroy(pf)

    end subroutine simulate

  end program pfasst_imk


text:
C...GETISSUD performs an interpolation/extrapolation in 3 dimensions by
C...fitting quadratic splines using 4 points in each dimension
      double precision function getissud(ibeam,kfl,x1,x2,pt2)
      implicit none

      include 'sudgrid.inc'

c   Arguments
      integer ibeam,kfl
      double precision x1,x2,pt2
c   Storing values for the interpolation
      double precision smallgrid(4,4),minigrid(4) ! pt2,x1
c   Local variables
      integer ipt2,ix1,ix2,ilo,ihi,i,j,k,kkfl,ipoints
      double precision pt2i,x2i,x1i,minpoint,maxpoint,x(3)
      integer nerr
      data nerr/0/

      getissud=0

      x(1)=log(x2)
      x(2)=x1
      x(3)=log(pt2)

      kkfl=kfl
      if(ibeam.lt.0) kkfl=-kkfl
      if(kkfl.lt.-2) kkfl=iabs(kfl)
      if(iabs(kkfl).eq.21) kkfl=0
      if(kkfl.eq.5) then
        ipoints=2
      else
        ipoints=1
      endif
      if(kkfl.gt.5) then
        if(nerr.lt.10)
     $     write(*,*)'GETISSUD Warning: flavor ',kfl,' not supported'
        nerr=nerr+1
        getissud=1
        return
      endif

      if(x(1).lt.points(1,ipoints).or.
     $   x(1).gt.points(nx2,ipoints).and.nerr.lt.10)
     $   then
        write(*,*) 'GETISSUD Warning: extrapolation in x2: ',x2
        nerr=nerr+1
      endif

      if(x(2).lt.points(nx2+1,ipoints).or.
     $   x(2).gt.points(nx2+nx1,ipoints)
     $   .and.nerr.lt.10) then
        write(*,*) 'GETISSUD Warning: extrapolation in x1: ',x1
        nerr=nerr+1
      endif

      if(kkfl.eq.5.and.pt2.lt.22.3109)then
        getissud=1d0
        return
      endif

      if(kkfl.eq.5.and.x1.gt.0.6)then
        getissud=0d0
        return
      endif

      if(x(3).lt.points(nx2+nx1+1,ipoints)) then
        write(*,*) 'GETISSUD Error! pt2 = ',exp(x(3)),' < ',
     $     exp(points(nx2+nx1+1,ipoints)),' = min(pt2). Not allowed!'
        write(*,*) 'You need to regenerate grid with new pt2min.'
        stop
      endif

      if(x(3).lt.points(nx2+nx1+1,ipoints).or.
     $   x(3).gt.points(nx2+nx1+npt2,ipoints)
     $   .and.nerr.lt.10) then
        write(*,*) 'GETISSUD Warning: extrapolation in pt2: ',pt2
        nerr=nerr+1
      endif

c   Find nearest points by binary method
c   x2
      ilo=1
      ihi=nx2
      do while(ihi.gt.ilo+1)
        ix2=ilo+(ihi-ilo)/2
        if(x(1).gt.points(ix2,ipoints))then
          ilo=ix2
        else
          ihi=ix2
        endif
      enddo
      if(x(1).lt.points(ix2,ipoints))
     $   ix2=ix2-1
      ix2=max(2,min(ix2,nx2-2))

c      print *,'x2: ',ix2,x(1),(points(i,ipoints),i=ix2-1,ix2+2)

c   x1
      ilo=1
      ihi=nx1
      do while(ihi.gt.ilo+1)
        ix1=ilo+(ihi-ilo)/2
        if(x(2).gt.points(nx2+ix1,ipoints))then
          ilo=ix1
        else
          ihi=ix1
        endif
      enddo
      if(x(2).lt.points(nx2+ix1,ipoints))
     $   ix1=ix1-1
      ix1=max(2,min(ix1,nx1-2))

      do while(kkfl.eq.5.and.
     $   points(nx2+ix1+2,ipoints).gt.0.6)
        ix1=ix1-1
      enddo

c      print *,'x1: ',ix1,x(2),(points(nx2+i,ipoints),i=ix1-1,ix1+2)

c   pt2
      ilo=1
      ihi=npt2
      do while(ihi.gt.ilo+1)
        ipt2=ilo+(ihi-ilo)/2
        if(x(3).gt.points(nx2+nx1+ipt2,ipoints))then
          ilo=ipt2
        else
          ihi=ipt2
        endif
      enddo
      if(x(3).lt.points(nx2+nx1+ipt2,ipoints))
     $   ipt2=ipt2-1
      ipt2=max(2,min(ipt2,npt2-2))

      do while(kkfl.eq.5.and.
     $   exp(points(nx2+nx1+ipt2-1,ipoints)).lt.22.3109)
        ipt2=ipt2+1
      enddo
c      print *,'pt2: ',ipt2,x(3),(points(nx2+nx1+i,ipoints),i=ipt2-1,ipt2+2)
c      print *,'pt: ',ipt2,exp(x(3)/2),
c     $   (exp(points(nx2+nx1+i,ipoints)/2),i=ipt2-1,ipt2+2)

C   Now perform inter-/extra-polation

C   Start with x2, which should have the flattest distribution
C   Calculate sud(x2,ax1,apt2) for the 4x4 apt2 and ax1
C   Then continue with pt2 and calculate sud(x2,ax1,pt2)
C   for the 4 ax1
C   Finally calculate sud(x2,x1,pt2)

      do i=1,4
        do j=1,4
c          print *,'x1,pt:',points(nx2+ix1-2+i,ipoints),
c     $       exp(points(nx2+nx1+ipt2-2+j,ipoints)/2)
          call splint2(sudgrid(ix2-1,ix1-2+i,ipt2-2+j,kkfl),
     $       points(ix2-1,ipoints),4,x(1),smallgrid(j,i))
          smallgrid(j,i)=max(0d0,min(1d0,smallgrid(j,i)))
        enddo
      enddo

      do i=1,4
        call splint2(smallgrid(1,i),
     $     points(nx2+nx1+ipt2-1,ipoints),4,x(3),minigrid(i))
          minigrid(i)=max(0d0,min(1d0,minigrid(i)))
      enddo

      call splint2(minigrid,
     $   points(nx2+ix1-1,ipoints),4,x(2),getissud)
      getissud=max(0d0,min(1d0,getissud))

c      print *,'Result: ',getissud

      return
      end

      subroutine splint2(ypoints,xpoints,npoints,x,ans)
      implicit none

C   arguments
      integer npoints
      double precision ypoints(npoints),xpoints(npoints)
      double precision x,ans
C   local variables
      double precision a0,a1,a2,sd
      integer ifail,i

      CALL DLSQP2(npoints,xpoints,ypoints,a0,a1,a2,sd,ifail)

c      print *,'Point, interpolation:'
c      do i=1,npoints
c        print *,exp(xpoints(i)),ypoints(i),
c     $     a0+a1*xpoints(i)+a2*xpoints(i)**2
c      enddo

      ans=a0+a1*x+a2*x**2
c      print *,x,ans

      return
      end


text:
MODULE read_ncoda_prep
!===============================================================================
! This program reads the NCODA prep files:
! coda.MVO_prp.*
! coda.SSH_prp.*
! coda.MOV_ENS_obs.*
! coda.SSH_ENS_obs.*
!
! These containt the observations (prp) and ensemble member innovations (ENS_obs)
!
! This routine reads both files and recovers the observation and model
! equivalent for a specified member.
!
! The obsop.f90 program uses this module to read these in directly.
!
! Observation errors are provided in the 'prp' files.
!
! Author:
! 12/22/17:  Stephen G. Penny, University of Maryland, College Park
!                              Visiting Scientist, Naval Research Laboratory - Stennis Space Center
!
!===============================================================================

USE common,                     ONLY: r_sngl, r_size, slen
USE params_obs,                 ONLY: id_t_obs, id_s_obs, id_u_obs, id_v_obs, id_eta_obs
USE compute_profile_error,      ONLY: cmpTz

IMPLICIT NONE

PUBLIC :: read_ncoda_prp, ncoda_prep_data

PRIVATE

INTEGER :: nobs, nobs0
INTEGER :: i,j,k,n
REAL(r_size) :: se0, seF

TYPE ncoda_prep_data
  REAL(r_size) :: x_grd(3)  ! longitude, latitude, and z depth (m)
  REAL(r_size) :: value     ! actual physical value of the parameter measured at this grid point
  REAL(r_size) :: hxb       ! Model equivalent ! (NEMO) added
  REAL(r_size) :: lev       ! grid z level
  REAL(r_size) :: oerr      ! observation standard error
  REAL(r_size) :: hour      ! Hour of observation
  CHARACTER(9) :: plat      ! Platform
  CHARACTER(3) :: ptyp      ! Profile type
  CHARACTER(3) :: sid       ! Source id
  CHARACTER(1) :: qkey      ! Quality key
  INTEGER :: qc             ! Quality control flag  ! (NEMO) added
  INTEGER :: typ    ! observation variable type (e.g., PRES_TYPE)
  INTEGER :: nlevs  ! number of levels with data, counting from the top, including levels with missing data that have obs below them.
  INTEGER :: id     ! id number used in observation files to identify the observation
  INTEGER :: rid    ! id of the record, in order that it is read in
  INTEGER :: lid    ! id of the level for each record (upon skipping missing data for some levels)
  LOGICAL :: kept   ! tells letkf whether this obs is kept for assimilation
END TYPE ncoda_prep_data

LOGICAL :: DO_READ_OBE      ! Read the observation error estimate from the file
LOGICAL :: DO_COMPUTE_OBE   ! Compute the observation error estimate using the NCEP approach w/ vertical gradient

INTEGER :: undef = 99999

!! Write letkf file
!do i=1,nobs
!!STEVE: the following are required for miyoshi's letkf observation input format:
!!1 = obelm
!!2 = lon
!!3 = lat
!!4 = lev
!!5 = value
!!6 = oberr
! wk(1) = obs_data(i)%typ
! wk(2) = obs_data(i)%x_grd(1)
! wk(3) = obs_data(i)%x_grd(2)
! wk(4) = obs_data(i)%x_grd(3)
! wk(5) = obs_data(i)%value
! wk(6) = obs_data(i)%oerr
! WRITE(fid) wk
!enddo

CONTAINS

SUBROUTINE read_ncoda_prp(infile,infile2,obs_data,nobs,mem,obid_in)
!===============================================================================
! Read the ncoda prep observation data and model equivalent
!===============================================================================

USE params_letkf, ONLY: nbv

CHARACTER(*), INTENT(IN) :: infile, infile2
TYPE(ncoda_prep_data), INTENT(OUT), ALLOCATABLE, DIMENSION(:) :: obs_data
INTEGER, INTENT(OUT) :: nobs
INTEGER, INTENT(IN) :: mem ! Specify ensemble member to process
INTEGER, INTENT(IN), OPTIONAL :: obid_in

! Other variables:
INTEGER :: i,j,k,n
REAL(r_sngl), ALLOCATABLE, DIMENSION(:) :: xlon, ylat, hour
CHARACTER(3), ALLOCATABLE, DIMENSION(:) :: ptyp, sid
CHARACTER(9), ALLOCATABLE, DIMENSION(:) :: plat
CHARACTER(1), ALLOCATABLE, DIMENSION(:) :: qkey
INTEGER, ALLOCATABLE, DIMENSION(:) :: qc
REAL(r_sngl), ALLOCATABLE, DIMENSION(:) :: ovals, mvals, stde
REAL(r_sngl), ALLOCATABLE, DIMENSION(:) :: depth ! profile depths  (dimensioned: nobs, depth)
REAL(r_sngl) :: val, hxb, err
INTEGER :: cnt, nlv
LOGICAL :: dodebug=.true.
REAL(r_sngl) :: missing_value=-999  ! Missing value flag in the data
REAL(r_sngl) :: max_value=999       ! Maximum value condition
REAL(r_sngl) :: max_depth = 9999    ! Maximum depth condition

INTEGER :: fid = 21
LOGICAL :: ex
INTEGER :: obid

REAL(r_sngl) :: obs_anm
REAL(r_sngl), ALLOCATABLE :: obs_xi (:)
REAL(r_sngl), ALLOCATABLE :: obs_yj (:)
REAL(r_sngl), ALLOCATABLE :: obs_zk (:)
INTEGER,      ALLOCATABLE :: otype (:)
INTEGER,      ALLOCATABLE :: otime (:)
INTEGER,      ALLOCATABLE :: odate (:)
REAL(r_sngl), ALLOCATABLE :: ens_anm(:)
REAL(r_sngl), ALLOCATABLE :: dum_anm(:)
INTEGER,      ALLOCATABLE :: etyp (:)
INTEGER,      ALLOCATABLE :: obs_var (:)

!var is the variable 1-6 according to:
!     data      var_lbl / 'seatmp', 'salint', 'geoptl',
!    *                    'uucurr', 'vvcurr', 'lyrprs' /

! Check optional argument
if (present(obid_in)) then
  obid = obid_in
else
  obid = id_t_obs ! (Default, reads both id_t_obs and id_s_obs)
endif

!-------------------------------------------------------------------------------
! Open obs prep file
!-------------------------------------------------------------------------------
inquire (file=trim(infile), exist=ex)
if (ex) then
  open (unit=fid,file=trim(infile),status='old',form='unformatted')
! rewind(fid)
else
  print *, "infile = ", trim(infile)
  STOP('ERROR: prp file is missing. EXITING...')
endif

!-------------------------------------------------------------------------------
! Read the number of obs
!-------------------------------------------------------------------------------
read (fid) cnt
nobs=cnt !STEVE: ?
print *, 'read innovation file 1: ', trim(infile)
print *, "nobs, cnt, cnt-nobs = ", nobs, cnt, cnt-nobs

!-------------------------------------------------------------------------------
! Allocate arrays for obs data
!-------------------------------------------------------------------------------
ALLOCATE (ovals(nobs))
ALLOCATE (stde(nobs))
ALLOCATE (xlon(nobs))
ALLOCATE (ylat(nobs))
ALLOCATE (depth(nobs))
ALLOCATE (hour(nobs))
ALLOCATE (etyp(nobs))
ALLOCATE (qc(nobs))

ALLOCATE (obs_xi(nobs))
ALLOCATE (obs_yj(nobs))
ALLOCATE (obs_zk(nobs))
ALLOCATE (obs_var(nobs))
ALLOCATE (otype(nobs))
ALLOCATE (otime(nobs))
ALLOCATE (odate(nobs))

!-------------------------------------------------------------------------------
! Read the obs data
!-------------------------------------------------------------------------------
WRITE(6,*) "Reading obs data, cnt = ", cnt
if (cnt .gt. 0) then
  do i = 1, nobs
!   print *, "start reading obs i = ", i
    read (fid) odate(i), otime(i), obs_var(i), otype(i), &
     &     ovals(i), obs_anm, stde(i), depth(i), ylat(i), &
     &     xlon(i), obs_xi(i),  obs_yj(i),  obs_zk(i)

    hour(i) = otime(i)

    ! Specify the observation type:
    qc(i) = 1
    if (obid == id_t_obs) then
      select case (obs_var(i))
        case(1)
          etyp(i) = id_t_obs
        case(2)
          etyp(i) = id_s_obs
        case default
!         print *, "Observation variable type not recognized: ", obs_var(i)
!         print *, "REMOVING..."
          qc(i) = 0
      end select
    else
      ! Use this for ssh or sea ice ncoda-prep data:
      etyp(i) = obid
    endif
  enddo

  print *, "final nobs, xlon,ylat,depth = ", nobs, xlon(nobs),ylat(nobs),depth(nobs)
  print *, "obs_var,ovals,stde,obs_anm = ", obs_var(nobs),ovals(nobs),stde(nobs),obs_anm

else
  WRITE (6,*) "cnt==0, EXITING..."
  STOP(13)
endif

!-------------------------------------------------------------------------------
! Close the obs file
!-------------------------------------------------------------------------------
close(fid)

!-------------------------------------------------------------------------------
! Open ens obs prep file
!-------------------------------------------------------------------------------
inquire (file=trim(infile2), exist=ex)
if (ex) then
  open (unit=fid,file=trim(infile2),status='old',form='unformatted')
! rewind(fid)
else
  print *, "infile2 = ", trim(infile2)
  STOP('ERROR: ENS_obs file is missing. EXITING...')
endif
print *, 'read innovation file 2: ', trim(infile2)

!-------------------------------------------------------------------------------
! Allocate arrays for innovation data
!-------------------------------------------------------------------------------
ALLOCATE (ens_anm(nobs))
ALLOCATE (dum_anm(4)) !nbv))

!-------------------------------------------------------------------------------
! Read the ensemble innovation file
!-------------------------------------------------------------------------------
if (cnt .gt. 0) then
  do i = 1, nobs
!   print *, "start reading ens_anm i = ", i
    read (fid) dum_anm
    ens_anm(i) = dum_anm(mem)
!   print *, "after reading ens_anm i = ", i
  enddo
  print *, "ens_anm(1) = ",  ens_anm(1)
  print *, "ens_anm(nobs-1) = ", ens_anm(nobs-1)
endif

!-------------------------------------------------------------------------------
! Close the ensemble innovation file
!-------------------------------------------------------------------------------
close(fid)

! CONVERT data to ncoda_prep_data format:
print *, "Finished reading NCODA prep files, formatting data..."

! Loop through all profiles and create a new observation for each unique data point
!nobs = cnt * nlv
print *, "ALLOCATING obs_data with nobs = ", nobs
ALLOCATE(obs_data(nobs))
n = 0
do i=1,cnt
  if (qc(i)==0) CYCLE  ! Skip any obs deemed unusable
  if (dodebug) print *, "i = ", i

  ! Assign observation value:
  val = ovals(i)
  err = stde(i)
! if (etyp(i)==id_eta_obs) then
!   hxb = ens_anm(i)
! else
  hxb = val - ens_anm(i)
! endif

  if (abs(val) < abs(max_value) .and. abs(val) < abs(missing_value)-1 .and. depth(i) < max_depth) then
    n = n+1
    if (dodebug) print *, "n,lon,lat,depth,hour,val,err,hxb,qc = ", n,xlon(i),ylat(i),depth(i),hour(i),val,err,hxb,qc(i)

    obs_data(n)%typ = etyp(i)
    obs_data(n)%x_grd(1) = xlon(i)
    obs_data(n)%x_grd(2) = ylat(i)
    obs_data(n)%x_grd(3) = depth(i)
    obs_data(n)%hour = hour(i)
    obs_data(n)%value = val
    obs_data(n)%hxb   = hxb
    obs_data(n)%oerr = err
    obs_data(n)%rid = obs_var(i) !STEVE: ?
    obs_data(n)%lid = obs_zk(i)  ! level id
  endif
enddo

nobs = n
if (dodebug) print *, "nobs = ", nobs

! Explicitly deallocate temporary arrays
if (ALLOCATED(depth)) DEALLOCATE(depth)
if (ALLOCATED(xlon))  DEALLOCATE(xlon)
if (ALLOCATED(ylat))  DEALLOCATE(ylat)
if (ALLOCATED(hour))  DEALLOCATE(hour)
if (ALLOCATED(qc))    DEALLOCATE(qc)
if (ALLOCATED(ovals))  DEALLOCATE(ovals)
if (ALLOCATED(ens_anm))  DEALLOCATE(ens_anm)
if (ALLOCATED(stde))  DEALLOCATE(stde)

if (dodebug) print *, "Temporary arrays deallocated."
if (dodebug) print *, "Returning from read_ncoda_prp..."

END SUBROUTINE read_ncoda_prp

END MODULE read_ncoda_prep

max_stars_repo_name: GEOS-ESM/Ocean-LETKF

max_issues_repo_name: GEOS-ESM/Ocean-LETKF

max_forks_repo_name: GEOS-ESM/Ocean-LETKF