!///////////////////////////////////////////////////////////////////////////////
module mdl_setting
   implicit none
   !--- Physical constants & parameters frequently used
   double precision, parameter :: hbarc = 197.329d0 !MeV.fm
   !Note: since this code is composed consistently to (MeV.fm) as well as (eV.nm) convensions, one can rearrange the energy & length units as
   !MeV <--> eV
   ! fm <--> nm, without really modifying the code.
   double precision, parameter :: pi = 3.141592653d0
   double precision, parameter :: alpha = 7.2973525698d-3 !fine structure constant, ~1/137.
   double precision, parameter :: anmass = 938.767d0 !MeV.c^2, one-nucleon mass
   complex (kind=kind(0d0)), parameter :: ai = (0.d0, 1.d0) !Imaginary unit
   !--- Model parameters
   double precision, parameter :: ecut = 10.d0 !MeV
   double precision, parameter :: tol = 1.d-10
   integer, parameter :: nss = 0 !Total spin of 2F.
   integer, parameter :: npp = 1 !Total parity of 2F. "=1" means "+", "=-1" means "-".
   integer, parameter :: ndmax = 40 !max # of nodes for single-particle states
   integer, parameter :: nspmax = 100 !max # of s.p. states
   integer, parameter :: nconfmax = 1000 !max # of configured basis
   double precision, parameter :: RMAX = 80.d0 !fm, maximum length of 1D space.
   double precision, parameter :: dr = 0.1d0 !fm, length mesh.
   integer, parameter :: nrmax = int(rmax/dr + 1.d-6) !max # of r-grid
   !--- parameters used in C-N WS and N-N SDDC pairing potentials
   double precision, parameter :: r00 = 1.27d0 !1.27d0 !1.06d0 !fm
   double precision, parameter :: b00 = 0.67d0 !0.67d0 !0.615d0 !fm
   double precision, parameter :: nc =  28.d0 !neutron-number of core
   double precision, parameter :: zc =  28.d0 ! proton-number of core
   double precision, parameter :: ac = zc + nc
   double precision, parameter :: xmu = anmass !!MeV !*ac/(ac+1.d0)
end module mdl_setting
!///////////////////////////////////////////////////////////////////////////////
!///////////////////////////////////////////////////////////////////////////////
module mdl_objects
contains
!*****************************************************************
   subroutine spbasis(label, nsp, ip, node, esp, psi)
      use mdl_setting, only: nspmax, nrmax
      implicit none
      character (4), intent (in) :: label
      integer, intent (inout) :: nsp, ip(nspmax), node(nspmax)
      double precision, intent (inout) :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer, save :: nspd, ipd(nspmax), noded(nspmax)
      double precision, save :: espd(nspmax), psid(nspmax, -nrmax:nrmax)
      if (label=='save') then
         nspd = nsp ; ipd = ip ; noded = node ; espd = esp
         psid = psi
      else if (label=='load') then
         nsp = nspd ; ip = ipd ; node = noded ; esp = espd
         psi = psid
      end if
   end subroutine spbasis
!*****************************************************************
   subroutine confbasis(label, nconf, nncf1, nncf2, ecf12, fac12)
      use mdl_setting, only : nconfmax
      implicit none
      character (4), intent (in) :: label
      integer, intent(inout) :: nconf, nncf1(nconfmax), nncf2(nconfmax)
      double precision, intent(inout) :: ecf12(nconfmax), fac12(nconfmax)
      integer, save :: nb, nn1(nconfmax), nn2(nconfmax)
      double precision, save :: e12(nconfmax), fcf(nconfmax)
      if (label=='save') then
         nb = nconf ; nn1 = nncf1 ; nn2 = nncf2 ; e12 = ecf12 ; fcf = fac12
      else if (label=='load') then
         nconf = nb ; nncf1 = nn1 ; nncf2 = nn2 ; ecf12 = e12 ; fac12 = fcf
      end if
   end subroutine confbasis
!*****************************************************************
end module mdl_objects
!///////////////////////////////////////////////////////////////////////////////
!///////////////////////////////////////////////////////////////////////////////
module mdl_routines
   use mdl_setting, only : hbarc, alpha, r00, b00, ac, zc, nc
   implicit none
   private :: eigsrt
   double precision, parameter :: rc = r00*ac**(0.333333333333d0)
   double precision, parameter :: bc = rc
   double precision, parameter :: v0 = -40.d0 !-52.5d0 MeV
!!   double precision, parameter :: vpair = 10.d0 !MeV for the SW-2p interaction.
   double precision, parameter :: wpair = -10.d0 !MeV*fm for the function "ddpair".
contains
!****************************************************************
	function ddpair(x)  result(ff) ! Pairing potential
        double precision, intent(in) :: x
        double precision :: ff
        ff = wpair*1.d0
!!        ff = wpair*(1.d0-1.d0/( 1.d0 + exp((abs(x)-rc)/b00)) ) !density-dependent
!!        ff = wpair*(1.d0/( 1.d0 + exp((abs(x)-rc)/b00)) - 1.d0) + wpair !density-dependent
	return
	end function
!****************************************************************
	function vcn_coul(x) result(ff) ! Soft Coulomb potential for p-C
        double precision, intent(in) :: x
        double precision :: ff
        ff = hbarc*alpha*zc / sqrt(bc*bc + x*x)
	end function
!****************************************************************
	function vcn_nucl(x)  result(ff) ! Woods-Saxon potential for p-C
        double precision, intent(in) :: x
        double precision :: ff
        ff = v0/( 1.d0 + exp((abs(x)-rc)/b00) )
	end function
!****************************************************************
	function vcn_adel(x)  result(ff) ! Additional Barrier
        double precision, intent(in) :: x
        double precision :: ff, a, ra, va0
        a = abs(x) ; ra = 14.d0 !fm
        va0 = 100.d0 !MeV*fm*fm
        ff = va0*a*a/(a*a*a*a+ra*ra*ra*ra)*1.d0
	end function
!****************************************************************
        function vcn(x)  result(ff) ! Total potential for p-C
        double precision, intent(inout) :: x
        double precision :: ff, v1, v2, a1, a2, b2
!        ff = vcn_nucl(x) + vcn_coul(x) + vcn_adel(x)
!        ff = vgauss(x) + vcn_nucl(x) + vcn_adel(x)
        ff = vcn_fall(x)
        end function
!****************************************************************
	function vmod(x)  result(ff) ! Modified Pot. for p-C.
        double precision, intent(inout) :: x
        double precision :: ff, d
        d = 11.d0
        ff = vcn(x) != vcn_nucl(x) + vcn_coul(x) + vcn_adel(x)
        if (abs(x) > d) ff = vcn(d) !vcn_nucl(d) + vcn_coul(d) + vcn_adel(d)
        end function
!****************************************************************
	function vcn_fall(x) result(ff) ! Soft Coulomb potential for p-C
        double precision, intent(in) :: x
        double precision :: ff, vega, vale, rale, riff
        vega = v0
        vale = 16.d0 !MeV
        rale =  9.d0 !fm
        riff = 12.d0 !fm
        ff = vega
        if(abs(x) > rale) ff = vale
        if(abs(x) > riff) ff = 0.d0
	end function

!****************************************************************
      subroutine sp_numerov ! single-particle wave functions
      use mdl_setting, only : nrmax, nspmax, ecut, tol, dr, rmax
      use mdl_objects, only : spbasis
      implicit none

      integer :: npi(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: nst, ip, nodemax, inode, in, nn, j, ix, jp
      double precision :: psi00(0:nrmax), e, emin, emax, xmax
      nst = 0
      do 30 ip = 0,1 !parity: "ip=1" <--> "psi(x=0)=0.d0"
         e = ecut ; jp = ip
         call numerov1(e,jp, inode, psi00)
         nodemax = inode-1

         emin = -abs(vcn_nucl(0.d0))
         emax = ecut
         nn = int( log10(abs(emax-emin)/tol)/log10(2.d0)+0.5d0 ) + 40

         do 40 in = 0, nodemax
            !--- s.p. energy
            do 503 j = 1, nn
               e = (emin+emax)*0.5d0
               call numerov1(e,jp, inode, psi00)
               if(inode .gt. in) then
               emax = e
               else
               emin = e
               endif
            503 enddo

            if (e < 0.d0) go to 507 !exclude bound states
!!            if ((e < -30.d0)) go to 507 !exclude bound states
            xmax = rmax
            if(e .lt. vcn(xmax)) call numerov2(e, psi00)

            nst = nst + 1
            if(nst.gt.nspmax) then
               write(6,*) 'Increase "nstatemax" (~_~;)' ; stop
            endif

            esp( nst) = e
            npi( nst) = ip
            node(nst) = 2*in+ip !2*inode+ip
            do ix = 0, nrmax
               psi(nst, ix) = psi00(ix) !0 =< x =< xmax
               psi(nst, ix-nrmax) = ((-1.0)**ip)*psi00(nrmax-ix) !-xmax =< x =< 0
            enddo

       507  emin = e
            emax = ecut
 40      enddo
 30   enddo

      call spbasis('save', nst, npi, node, esp, psi)
      return
      end subroutine sp_numerov

!**************************************************************
      subroutine rkutta(e,psi0,psi1)
      ! to get wave function at x=0 and x=dx for a positive parity state, 
      ! which satisfies psi(0)=1, and d/dx psi(0)=0.d0.
      use mdl_setting, only : hbarc, pi, xmu, dr
      implicit none
      double precision, intent(inout) :: e, psi0, psi1
      double precision :: fac, ak1, bk1, ak2, bk2, ak3, bk3, ak4, bk4, w, dx, dxhalf
      fac=dr*(2.d0*xmu/hbarc/hbarc) ; psi0=1.d0 ; w=0.d0
      dx = dr ; dxhalf = dr*0.5d0
      ak1=-fac*(e-vcn(dx))*psi0                 ; bk1=dx*w
      ak2=-fac*(e-vcn(dxhalf))*(psi0+0.5d0*bk1) ; bk2=dx*(w+0.5d0*ak1)
      ak3=-fac*(e-vcn(dxhalf))*(psi0+0.5d0*bk2) ; bk3=dx*(w+0.5d0*ak2)
      ak4=-fac*(e-vcn(dx))*(psi0+bk3)           ; bk4=dx*(w+ak3)
      !--- wave functions at dx
      psi1 = psi0 + 1.d0/6.d0*(bk1 + 2.d0*bk2 + 2.d0*bk3 + bk4)
      return
      end subroutine rkutta

!****************************************************************
   subroutine numerov1(e,ip,inode,psi00)
   ! Subroutine for solution of the Schroedinger eq. by Numerov method.
   ! [(d/dx)^2 + f(x)] = 0,   f(x) = [V(x)-E] / (-hb*hb/2u).
      use mdl_setting, only : hbarc, pi, dr, nrmax, xmu
      implicit none
      integer, intent(inout) :: ip, inode
      double precision, intent(inout) :: e, psi00(0:nrmax)
      integer :: ix
      double precision :: psi0, psi1, fac, aa,cc,dd, x,x0,x1, fx,f0,f1, xmax

      fac = -hbarc*hbarc*0.5d0/xmu
      psi0=0.d0 ; psi1=1.d0                !for "parity = 1, odd"
      if(ip.eq.0) call rkutta(e,psi0,psi1) !for "parity = 0,even"

      psi00(0)=psi0 ; psi00(1)=psi1 ; inode = 0
      do 29 ix = 2, nrmax
         x0=dr*(ix-2)       ; x1=dr*(ix-1)       ; x=dr*ix
         f0=(vcn(x0)-e)/fac ; f1=(vcn(x1)-e)/fac ; fx=(vcn(x)-e)/fac

         dd = psi1 * (2.d0 - 5.d0*dr*dr*f1/6.d0) - psi0 * (1.d0 + dr*dr*f0/12.d0)
         cc = 1.d0 + dr*dr*fx/12.d0
         psi00(ix) = dd/cc

         if(ix.gt.5.and.psi00(ix)*psi00(ix-1).lt.0.) inode=inode+1

         psi0 = psi1
         psi1 = psi00(ix)
 29   enddo

     !--- normalization
     aa = 0.d0
     do ix = 1, nrmax
        aa = aa + psi00(ix)**2*dr
     enddo
     aa = aa + psi00(0)**2*dr/2.d0
     psi00 = psi00/sqrt(aa)*sqrt(0.5d0)
     return
     end subroutine numerov1

!****************************************************************
      subroutine numerov2(e, psi00)
      !--- Solve the Schroedinger eq. backward in order to ensure the asymptotic form for E < 0.
      use mdl_setting, only : hbarc, pi, dr, nrmax, rmax, xmu
      implicit none
      double precision, intent(inout) :: e, psi00(0:nrmax)
      integer :: ix, imatch
      double precision :: psi0, psi1, fac, aa,cc,dd, x,x0,x1, fx,f0,f1, &
          &               xmatch, psi(0:nrmax), ak, cmatch, xmax
      xmatch = 1.d0 !fm !0.1d0*rmax
      imatch = nint(xmatch/dr)

      fac = -hbarc*hbarc*0.5d0/xmu
!1      fac=dr**2*(2.d0*xmu/hbarc**2)

      xmax = rmax ; ak=sqrt(2.d0*xmu*abs(e-vcn(xmax)))/hbarc
      psi(nrmax)   = exp(-ak*(nrmax*dr)) ; psi(nrmax-1) = exp(-ak*(nrmax-1)*dr)
      psi0 = psi(nrmax)                  ; psi1 = psi(nrmax-1)

      do 29 ix = nrmax-2, imatch, -1
         x=dr*ix           ; x0=dr*(ix+2)       ; x1=dr*(ix+1)
         fx=(vcn(x)-e)/fac ; f0=(vcn(x0)-e)/fac ; f1=(vcn(x1)-e)/fac

         dd = psi1 * (2.d0 - 5.d0*dr*dr*f1/6.d0) - psi0 * (1.d0 + dr*dr*f0/12.d0)
!1         dd = psi0 - fac/12.d0*(vcn(x0)-e)*psi0 - 2.d0*psi1 - 5.d0*fac/6.d0*(vcn(x1)-e)*psi1
         cc = 1.d0 + dr*dr*fx/12.d0
!1         cc = -fac/12.d0*(vcn(x)-e)+1.d0
         psi(ix) = dd/cc
!1         psi(ix)=-dd/cc

         psi0 = psi1
         psi1 = psi(ix)
 29   enddo

        !--- matching
        cmatch = psi00(imatch)/psi(imatch)
        do ix = imatch, nrmax
           psi00(ix) = psi(ix)*cmatch
        enddo

        !--- normalization
        aa = 0.d0
        do ix = 1, nrmax
           aa = aa + psi00(ix)**2*dr
        enddo
        aa = aa + psi00(0)**2*dr/2.d0
        psi00 = psi00/sqrt(aa)*sqrt(0.5d0)

        return
      end subroutine numerov2

!**************************************************************
      subroutine cnfbss ! Construct the configurated basis
      use mdl_setting, only : nrmax, nspmax, nconfmax, ecut, nss, npp
      use mdl_objects, only : spbasis, confbasis
      implicit none
      integer :: nsp, ip(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), faccf(nconfmax)
      integer :: n, i1, i2, j1, j2
      double precision :: fac12
      call spbasis('load', nsp, ip, node, esp, psi)
      n=0
      do j1=1, nsp
      do j2=j1,nsp
         if(esp(j1)+esp(j2).gt.ecut) goto 205 !energy-cutoff
         i1=ip(j1) ; i2=ip(j2)
         !=== Filters;
         if((nss.eq.1) .and. (j1.eq.j2)) goto 205 !for |S=1> <==> 2p cannot get the same orbit.
         if((npp > 0) .and. (mod(i1+i2,2).eq.1)) goto 205 !2p-parity=(+) only
         if((npp < 0) .and. (mod(i1+i2,2).eq.0)) goto 205 !2p-parity=(-) only
         !===
         n=n+1
         if(n.gt.nconfmax) then
            write(6,*)'Increase nconfmax !!' ; stop
         endif
         nncf1(n)=j1 ; nncf2(n)=j2 ; ecf12(n)=esp(j1)+esp(j2)
                      fac12 = 1.d0/sqrt(2.d0)
         if(j1.eq.j2) fac12 = 0.5d0
         faccf(n)=fac12 !normalization factor
 205  continue
      enddo
      enddo
      call confbasis('save',n, nncf1,nncf2, ecf12,faccf)
      return
      end subroutine cnfbss

!*****************************************************************
   subroutine mdiag(a, d, vv, n, np)
! The Jacobi method to compute all eigenvalues and eigenvectors
! of a real symmetric matrix A, which is of size N by N, stored
! in a physical NP by NP array. On output, the elements of A
! above the diagonal are destroyed. D returns the eigenvalues of
! A in its first N elements. V is a matrix with the same logical
! and physical dimensions as A whose columns contain, on output,
! the normalized eigenvectors of A. NROT returns the number of
! Jacobi rotations which were required.
! The I-th component of the eigenvector for the K-th eigenvalue 
! is given by V(I,K).
      implicit none
      integer :: n, np, ip, iq, nrot, i, j
      double precision :: c, g, h, s, sm, t, tau, theta, tresh
      double precision, dimension (np) :: b, z, d
      double precision, dimension (np, np) :: a, vv
      do ip = 1, n
         do iq = 1, n
            vv(ip, iq) = 0.d0
         end do
         vv(ip, ip) = 1.d0
      end do
      do ip = 1, n
         b(ip) = a(ip, ip)
         d(ip) = b(ip)
         z(ip) = 0.d0
      end do
      nrot = 0
      do i = 1, 50
         sm = 0.d0
         do ip = 1, n - 1
            do iq = ip + 1, n
               sm = sm + abs(a(ip,iq))
            end do
         end do
         if (sm==0.d0) then
            call eigsrt(d, vv, n, np)
            return
         end if
         if (i<4) then
            tresh = 0.2d0*sm/n**2
         else
            tresh = 0.d0
         end if
         do ip = 1, n - 1
            do iq = ip + 1, n
               g = 100.d0*abs(a(ip,iq))
               if ((i>4) .and. (abs(d(ip))+g==abs(d(ip))) .and. (abs(d(iq))+g==abs(d(iq)))) then
                  a(ip, iq) = 0.d0
               else if (abs(a(ip,iq))>tresh) then
                  h = d(iq) - d(ip)
                  if (abs(h)+g==abs(h)) then
                     t = a(ip, iq)/h
                  else
                     theta = 0.5d0*h/a(ip, iq)
                     t = 1.d0/(abs(theta)+sqrt(1.d0+theta**2))
                     if (theta<0.d0) t = -t
                  end if
                  c = 1.d0/sqrt(1.d0+t**2)
                  s = t*c
                  tau = s/(1.d0+c)
                  h = t*a(ip, iq)
                  z(ip) = z(ip) - h
                  z(iq) = z(iq) + h
                  d(ip) = d(ip) - h
                  d(iq) = d(iq) + h
                  a(ip, iq) = 0.d0
                  do j = 1, ip - 1
                     g = a(j, ip)
                     h = a(j, iq)
                     a(j, ip) = g - s*(h+g*tau)
                     a(j, iq) = h + s*(g-h*tau)
                  end do
                  do j = ip + 1, iq - 1
                     g = a(ip, j)
                     h = a(j, iq)
                     a(ip, j) = g - s*(h+g*tau)
                     a(j, iq) = h + s*(g-h*tau)
                  end do
                  do j = iq + 1, n
                     g = a(ip, j)
                     h = a(iq, j)
                     a(ip, j) = g - s*(h+g*tau)
                     a(iq, j) = h + s*(g-h*tau)
                  end do
                  do j = 1, n
                     g = vv(j, ip)
                     h = vv(j, iq)
                     vv(j, ip) = g - s*(h+g*tau)
                     vv(j, iq) = h + s*(g-h*tau)
                  end do
                  nrot = nrot + 1
               end if
            end do
         end do
         do ip = 1, n
            b(ip) = b(ip) + z(ip)
            d(ip) = b(ip)
            z(ip) = 0.d0
         end do
      end do
!      pause '50 iterations should never happen'
      return
   end subroutine mdiag
!*****************************************************************
   subroutine eigsrt(d, v, n, np)
      implicit none
      integer :: n, np, i, j, k
      double precision :: p
      double precision, dimension (np) :: d
      double precision, dimension (np, np) :: v
      do i = 1, n - 1
         k = i
         p = d(i)
         do j = i + 1, n
            if (d(j)<=p) then
               k = j
               p = d(j)
            end if
         end do
         if (k/=i) then
            d(k) = d(i)
            d(i) = p
            do j = 1, n
               p = v(j, i)
               v(j, i) = v(j, k)
               v(j, k) = p
            end do
         end if
      end do
      return
   end subroutine eigsrt
!*****************************************************************
   pure function kdel(i,j) result(ff)
   implicit none
   integer, intent(IN) :: i,j
   double precision :: ff
   ff = 0.d0
   if(i.eq.j) ff = 1.d0
   end function kdel
!*****************************************************************

end module mdl_routines
!///////////////////////////////////////////////////////////////////////////////




!--------------------------------------
   subroutine result_sp
      use mdl_setting, only : nrmax, dr, nspmax
      use mdl_objects, only : spbasis
      implicit none
      integer :: nsp, ip(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: i, ix
      call spbasis('load', nsp, ip, node, esp, psi)
      write(6,*)'Single Particle States'
      write(6,*)'  N_sp   Node   Parity    E_sp (MeV)'
      do i = 1, min(100, nsp)
         write(6,"(3(4x,i2),3x,f16.7)") i, node(i), ip(i), esp(i)
      enddo
      write(6,*)' Total # of sp-states =',nsp
      write(6,*)''
      do i = 1, min(9,nsp)
      do ix = -nrmax, nrmax
         write(200+i,"(X,F10.5,X,2E15.7)") dble(ix)*dr, psi(i,ix), esp(i)
      enddo
      enddo
   end subroutine

!--------------------------------------
   subroutine result_conf
      use mdl_setting, only : nconfmax
      use mdl_objects, only : confbasis
      implicit none
 94   format(3(5x,i3),3x,f16.7)
      integer :: nconf, nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), fac12(nconfmax)
      integer :: i

      call confbasis('load', nconf, nncf1, nncf2, ecf12, fac12)
      write(6,*)'Non-Correlated Two-Fermion States (shown up to 40)'
      write(6,*)'   N_2p   N_sp(1)    N_sp(2)      E_2p (MeV)'
      do i = 1, min(40,nconf)
         write(6,94) i, nncf1(i), nncf2(i), ecf12(i)
      enddo
      write(6,*)' Total # of non-correlated 2p-basis =', nconf
      write(6,*)''
   end subroutine


!********************************************************************
   subroutine hamatrix_tempo !--- SR to compute & diagonalize the total Hamiltonian.
      use mdl_setting, only : nconfmax, nspmax, dr, nrmax, nss
      use mdl_routines, only : mdiag, ddpair, vmod, vcn, kdel
      use mdl_objects, only : confbasis, spbasis
      implicit none
      integer :: nsp, ip(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: nconf, nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), fac12(nconfmax)
      integer :: i, j, k, ix, nb1, nb2, nk1, nk2, ix2
      double precision :: x, x1, x2, ddp, fact, factb
      integer :: min2, max2, idis, item
      double precision :: pket, pbra
      double precision, dimension(:), allocatable :: eall, emod, qall
      double precision, dimension(:,:), allocatable :: h0, hv, hw, ve, veigs2, veigs3
      complex(kind=kind(0d0)) :: cz, z1
      complex(kind=kind(0d0)), dimension(:,:), allocatable :: wf
      double precision, dimension(:,:), allocatable :: vsp_mod
      complex(kind=kind(0d0)), dimension(:), allocatable :: cvi_r2, cvi_r3

      call spbasis('load', nsp, ip, node, esp, psi)
      call confbasis('load', nconf, nncf1, nncf2, ecf12, fac12)
      !====================
      write(6,*)'Now calculating the matrix elements of Vmod(x)...'
      70003  allocate(vsp_mod(nsp,nsp))
      do i = 1, nsp !for "ket, |sp_i>"
         do j = i, nsp !!for "bra, <sp_j|"
            ddp = 0.d0
            do ix = -nrmax, nrmax
               x = dble(ix)*dr
               ddp  = ddp + psi(j,ix)*(vmod(x) - vcn(x))*psi(i,ix)
            end do
            ddp = ddp*dr
!!            if(abs(ddp)>1.d-2) write(6,*) ">>>>>>>>>>>",i,j,ddp
            vsp_mod(j,i) = ddp
            vsp_mod(i,j) = vsp_mod(j,i)
         end do
      end do
      !====================


      write(6,*)'Now calculating matrix elements of Hamiltonian...'
      pket = 1.d0 !phase-factor for "ket": |2p(i)> = |1>|2> + (pket)*|2>|1>
      if(nss.eq.1) pket = -1.d0
      pbra = pket

      70002 allocate(h0(nconf,nconf), hv(nconf,nconf), hw(nconf,nconf))
      h0 = 0.d0 ; hv =0.d0 ; hw = 0.d0

      do i = 1, nconf !for "ket, |i>"
         nk1 = nncf1(i) ; nk2 = nncf2(i)
         fact = fac12(i)
         do j = i, nconf !for "bra, <j|"
            nb1 = nncf1(j) ; nb2 = nncf2(j)
            factb = fac12(j)

            if (j.eq.i) then
               h0(j,i) = ecf12(j)
               h0(i,j) = h0(j,i)
            end if

            !--- confining potential
            ddp = fact*factb*( &
                  (1.d0 + pbra*pket) * kdel(nb1,nk1) * vsp_mod(nb2,nk2) &
                 +(1.d0 + pbra*pket) * kdel(nb2,nk2) * vsp_mod(nb1,nk1) &
                      +(pket + pbra) * kdel(nb1,nk2) * vsp_mod(nb2,nk1) &
                      +(pbra + pket) * kdel(nb2,nk1) * vsp_mod(nb1,nk2) )
            hw(j,i) = ddp
            hw(i,j) = ddp

            !============================
            !--- pp-interaction with finite range --> two-loop integration.
            if (1.eq.1) then !finite-range pairing
            ddp = 0.d0
            idis = 10
            do ix  = -nrmax, nrmax
               x = dble(ix)*dr
               min2 = max(ix-idis,-nrmax) !Range of ddpair(x)
               max2 = min(ix+idis, nrmax) !Range of ddpair(x)
               do ix2 = min2, max2
                  x2 = dble(ix2)*dr
                  ddp = ddp + fact*factb*(  &
                            psi(nb1,ix)*psi(nb2,ix2) * ddpair(x2) * psi(nk1,ix)*psi(nk2,ix2)  &
                     + pbra*psi(nb2,ix)*psi(nb1,ix2) * ddpair(x2) * psi(nk1,ix)*psi(nk2,ix2)  &
                     + pket*psi(nb1,ix)*psi(nb2,ix2) * ddpair(x2) * psi(nk2,ix)*psi(nk1,ix2)  &
                + pbra*pket*psi(nb2,ix)*psi(nb1,ix2) * ddpair(x2) * psi(nk2,ix)*psi(nk1,ix2)  )
               enddo
            enddo
            ddp = ddp*dr*dr
            end if
            !============================
            hv(j,i) = ddp
            hv(i,j) = ddp
         enddo
      enddo

      h0 = h0 + hv !original, TD-driver Hamiltonian.
      hw = hw + h0 !modified, confining Hamiltonian.

      !--- Diagonalize the original H:
      70001 allocate(eall(nconf), ve(nconf,nconf))
      eall = 0.d0 ; ve = 0.d0
      call mdiag(h0, eall, ve, nconf, nconf) !c_gs(i) = ve(i, 1), |gs> = \sum_{i} c_gs(i)*|i>.
      write(6,*) ' Original Hamiltonian has been diagonalized.'
      do i = 1, min(10,nconf)
         write (6, '(2X,A6,I3,A,F15.10,A7)') ' e_ex_', i - 1, '=', eall(i), '  (MeV)'
      end do
      write (6, *) ''

      !--- Modified H:
      70004 allocate(emod(nconf), veigs2(nconf,nconf), veigs3(nconf,nconf))
      emod = 0.d0 ; veigs2 = 0.d0
      call mdiag(hw, emod, veigs2, nconf, nconf)
      write(6,*) ' Modified Hamiltonian has been diagonalized.'
      do i = 1, min(10,nconf)
         write (6, '(2X,A6,I3,A,F15.10,A7)') ' e_ex_', i - 1, '=', emod(i), '  (MeV)'
      end do
      write (6, *) ''

      !--- Exp. coeffs. via non-correlated basis:
      70005 allocate(cvi_r2(nconf), cvi_r3(nconf), qall(nconf))
      write(6,*) "<><><><><> Initial state is prepared. <><><><><>"
      ITEM = 1 !find the 2p-state which you hope to see. "DEFAULT" = "G.S." with "item = 1".
      write(6,*) "--- Item =", item
      do i = 1, nconf
         cvi_r2(i) = cmplx(veigs2(i, ITEM), 0.d0, kind(0d0))
         ddp = abs(cvi_r2(i))**2
         if(ddp > 1.d-2) write(6,"(2X,'P_Non-Correlated_',I3,' = ',F9.4,X,'%')") i, 100.d0*ddp
      end do
      write (6, *) ""

   !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>
   if(1.eq.1) then
      write(6,*) "%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%"
      write(6,*) "%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%"
      write(6,*) "%%%%%%%%%%%%% Now the time-evolution starts ... %%%%%%%%%%%%%"
      write(6,*) ""

     !--- Probability-density at t=0:
     call FIG_DENSITY(nconf, cvi_r2, pket, 401)

     !--- Convert Exp. Coeffs. via H_td basis:
     veigs3 = 0.d0
        do i = 1, nconf
           do j = 1, nconf
              x1 = 0.d0
              do k = 1, nconf
                 x1 = x1 + ve(k,j)*veigs2(k,i)
              end do
              veigs3(j,i) = x1
           end do
        end do
        do i = 1, nconf
           cvi_r3(i) = cmplx(veigs3(i, ITEM), 0.d0, kind(0d0))
           ddp = abs(cvi_r3(i))**2
           if(ddp > 1.d-2) write (6, "(2X,'E_td_',I3,X,'=',F18.10,X,'MeV,  P=',F15.10,X,'%')") &
           i, eall(i), ddp*100.d0
        end do

     !--- Energy spectrum
     open(8494, file='Rectemp_002_Edist.dat', action='write', status='replace')
       do j = 1, nconf
          write(8494, *) Eall(j), veigs3(j,item)*veigs3(j,item)
       end do
     close(8494)
     write (6, *) ""
     !--- Q-values = <H_td> (=\= <Hmod> = Emod) for all the confined states.
     qall = 0.d0
        do j = 1, nconf
           x1 = 0.d0
           do i = 1, nconf
              x1 = x1 + veigs3(i,j)*Eall(i)*veigs3(i,j)
           end do
           qall(j) = x1
        end do
        do i = 1, min(12,nconf)
           write(6, "(2X,'Emod & Q-value_', I2, ' = ', F10.5, ' (MeV) & ', F10.5, ' (MeV).')") &
           i-1, emod(i), qall(i)
        end do
     write(6,*) ""
     !--- Pairing energy at t=0.
!     z1 = (0.d0, 0.d0)
!     call Kitaichi(nconf, cvi_r2, Hv, z1)
!     write(6,"(2X,'Pairing Energy (t=0) =', E10.4)") z1

     call TEMPO(nconf, cvi_r2, cvi_r3, eall, ve, hv)
   end if
   !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>


      80001 deallocate(eall,ve)
      80002 deallocate(h0,hv,hw)
      80003 deallocate(vsp_mod)
      80004 deallocate(emod,veigs2)
      80005 deallocate(cvi_r2,cvi_r3,qall)
   end subroutine hamatrix_tempo


!********************************************************************
   subroutine FIG_DENSITY(nconf, cv, pket, numb)
      use mdl_setting, only : nspmax, dr, nrmax, nconfmax
      use mdl_objects, only : spbasis, confbasis
      implicit none
      integer, intent (IN) :: nconf, numb
      double precision, intent(IN) :: pket
      complex(kind=kind(0d0)), intent(IN) :: cv(1:nconf)
      !-----------------
      integer :: nsp, ip(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: nd2, nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), fac12(nconfmax)
      character(32) :: filename
      integer :: ix, ix2, i, nk1, nk2
      double precision :: x, y1, y2, fact
      complex(kind=kind(0d0)) :: cz
      !-----------------
      call spbasis('load', nsp, ip, node, esp, psi)
      call confbasis('load', nd2, nncf1, nncf2, ecf12, fac12)
      !-----------------
!!      allocate(wf(-nrmax:nrmax, -nrmax:nrmax))
      i = numb
      write(filename, *) i ! i --> filename
!!      read (filename, *) i ! i <-- filename
      filename = 'Rec_'//trim(adjustl(filename))//'_Dens.dat'
      print *, filename
      open(401, file=filename, action='write', status='replace')
      do ix  = -nrmax, nrmax, 10
         write(401,*) ""
         do ix2 = -nrmax, nrmax, 10
            cz = (0.d0, 0.d0) !W.F. as the complex function.
            do i = 1, nconf !for "ket, |i>"
               nk1 = nncf1(i)
               nk2 = nncf2(i)
               fact = fac12(i)
               !cz = cz + veigs2(i,item)*fact*( psi(nk1,ix)*psi(nk2,ix2)+  pket * psi(nk2,ix)*psi(nk1,ix2) )
               cz = cz + cv(i)*fact*( psi(nk1,ix)*psi(nk2,ix2)+  pket * psi(nk2,ix)*psi(nk1,ix2) )
            enddo
!!            wf(ix,ix2) = cz
            write(401,'(6(2X,G13.6))') dble(ix)*dr, dble(ix2)*dr, abs(cz)**2, real(cz), aimag(cz)
         enddo
      enddo
      close(401)
!!      deallocate(wf)
   end subroutine fig_density




!********************************************************************
   subroutine RHO_TIDE(nconf, cv, pket, numb)
      use mdl_setting, only : nspmax, dr, nrmax, nconfmax
      use mdl_objects, only : spbasis, confbasis
      implicit none
      integer, intent (IN) :: nconf, numb
      double precision, intent(IN) :: pket
      complex(kind=kind(0d0)), intent(IN) :: cv(1:nconf)
      !-----------------
      integer :: nsp, ip(nspmax), node(nspmax)
      double precision :: esp(nspmax), psi(nspmax, -nrmax:nrmax)
      integer :: nd2, nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), fac12(nconfmax)
      character(32) :: filename
      integer :: ix, ix2, i, nk1, nk2, l1, l2
      double precision :: x1, x2, fact, border, d
      complex(kind=kind(0d0)) :: cz
      double precision, dimension(:), allocatable :: rsum
      !-----------------
      call spbasis('load', nsp, ip, node, esp, psi)
      call confbasis('load', nd2, nncf1, nncf2, ecf12, fac12)
!!      i = numb
!!      !-----------------
!!!      write(filename, *) i ! i --> filename
!!      read (filename, *) i ! i <-- filename
!!      filename = 'Rec_'//trim(adjustl(filename))//'_Dens.dat'
!!      print *, filename

1004  allocate(rsum(0:3))
      rsum = 0.d0
      border = 12.5d0 !fm
      do ix  = -nrmax, nrmax
         x1 = dble(ix)*dr
         l1 = 0
         if(x1 >  border) l1 =  1
         if(x1 < -border) l1 = -1

         do ix2 = -nrmax, nrmax
            x2 = dble(ix2)*dr
            l2 = 0
            if(x2 >  border) l2 =  1
            if(x2 < -border) l2 = -1

            cz = (0.d0, 0.d0) !W.F. as the complex function.
            do i = 1, nconf !for "ket, |i>"
               nk1 = nncf1(i)
               nk2 = nncf2(i)
               fact = fac12(i)
               !cz = cz + veigs2(i,item)*fact*( psi(nk1,ix)*psi(nk2,ix2)+  pket * psi(nk2,ix)*psi(nk1,ix2) )
               cz = cz + cv(i)*fact*( psi(nk1,ix)*psi(nk2,ix2)+  pket * psi(nk2,ix)*psi(nk1,ix2) )
            enddo
!!            wf(ix,ix2) = cz
            d = abs(conjg(cz)*cz)

            rsum(0) = rsum(0) + d
            if     (l1*l2 >  0) then !correlated
               rsum(1) = rsum(1) + d
            elseif (l1*l2 <  0) then !anti-correlated
               rsum(2) = rsum(2) + d
            elseif (l1+l2.ne.0) then !sequential
               rsum(3) = rsum(3) + d
            endif
         enddo
      enddo
      rsum = rsum*dr*dr
      open(401, file='Riders_006_PD_tide.dat', action='write', access='append')
  88  format(5g15.7)
      write(401, *) dble(numb), rsum(1), rsum(2), rsum(3), rsum(0)
      close(401)
2004  deallocate(rsum)
   end subroutine rho_tide



!********************************************************************
   subroutine TEMPO(nd1, cv2, cv3, Eall, VE, Hv) !VE is "veigs".
      use mdl_setting, only : ai, hbarc, nss
      use mdl_objects, only : confbasis, spbasis
      implicit none
      !---------
      integer, intent(IN) :: nd1
      double precision, intent(IN) :: Eall(nd1), VE(nd1, nd1), Hv(nd1, nd1)
      complex(kind=kind(0d0)), intent(IN) :: cv2(nd1), cv3(nd1)
      !---------
      double precision, parameter :: ctmax = 1000.d0 !fm
      double precision, parameter :: dct = 1.d0 !fm
      integer, parameter :: nctmax = int(ctmax/dct + 1.d-6)
      !---------
      integer :: i1, j1, ix1, nxmax, k, l, it
      double precision :: x1, tc, gamma, pket
      double precision, dimension(:), allocatable :: psurv, psdum, dsdum
      double precision, dimension(:,:), allocatable :: pndcy
      complex (kind=kind(0d0)) :: z1, zbeta, zod, zod2, zod3
      complex (kind=kind(0d0)), dimension(:), allocatable :: cvt_r2,cvt_d2 !via non-corre. basis
      complex (kind=kind(0d0)), dimension(:), allocatable :: cvt_r3,cvt_d3 !via H_td basis
      !---------
     70008 allocate(psurv(-nctmax:nctmax), pndcy(-nctmax:nctmax,3), &
               psdum(-nctmax+1:nctmax+1), dsdum(-nctmax+1:nctmax+1))
     70009 allocate(cvt_r2(nd1),cvt_r3(nd1),cvt_d2(nd1),cvt_d3(nd1))

   write(6,*) "<><><><><> Timewalk starts ... <><><><><>"
   psurv = 1.d0 ; pndcy = 0.d0
   do 35103 it = -nctmax, nctmax, 1
!??? TE1 ?????????????????????????????????????????????????????????????
      tc = dble(it)*dct
      zbeta = (0.d0, 0.d0) !--- note; zbeta = <reso'(t=0) | reso'(t)>.
      cvt_r3= (0.d0, 0.d0)
      do j1 = 1, nd1
         cvt_r3(j1) = exp(-ai*tc*Eall(j1)/hbarc)*cv3(j1) !--- {cv3} at t=it*dt via H_td basis.
         zbeta = zbeta + conjg(cv3(j1))*cvt_r3(j1)
      end do
      psurv(it) = abs(zbeta)**2 !--- suvival probability.
      cvt_r2= (0.d0, 0.d0)
      do j1 = 1, nd1
         z1 = (0.d0, 0.d0)
         do i1 = 1, nd1
            z1 = z1 + VE(j1,i1)*cvt_r3(i1)
         end do
         cvt_r2(j1) = z1 !--- {cv2} at t=it*dt via HO basis.
      end do

      !--- expan'coeff's for |decay(t)> and its norm ; zod = <decay(t) | decay(t)>.
      zod = (0.d0, 0.d0) !--- note; zod = <decay(t) | decay(t)> =/= 1.
      zod2= (0.d0, 0.d0)
      zod3= (0.d0, 0.d0)
      cvt_d2= (0.d0, 0.d0)
      cvt_d3= (0.d0, 0.d0)
      gamma = 1.8d0
      do i1 = 1, nd1
         cvt_d2(i1) = cvt_r2(i1) - zbeta*cv2(i1)
         cvt_d3(i1) = cvt_r3(i1) - zbeta*cv3(i1)
         zod  = zod  + conjg(cvt_d3(i1))*cvt_d3(i1)
         if(eall(i1).le.gamma) zod2 = zod2 + conjg(cvt_d3(i1))*cvt_d3(i1)
         if(eall(i1).gt.gamma) zod3 = zod3 + conjg(cvt_d3(i1))*cvt_d3(i1)
      end do
      pndcy(it,1) = abs(zod) !--- decay-probability from cvt_d3.
      pndcy(it,2) = abs(zod2)
      pndcy(it,3) = abs(zod3)

      IF(1.eq.1) THEN !plots of \abs{Psi(ct;x,y)}^2 : skip if it's time-consuming.
        pket = 1.d0 !phase-factor for "ket": |2p(i)> = |1>|2> + (pket)*|2>|1>
        if(nss.eq.1) pket = -1.d0
        IF((mod(it,200).eq.0) .and. (it>0) .and. (it<1400)) THEN
         write(6,*) "VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV"
         write(6,*) ">>> Density at ct = ", tc, " fm", it
         k = it+50000
         !--- Density
         call FIG_DENSITY(nd1, cvt_r2, pket, k)
         call RHO_TIDE(nd1, cvt_r2, pket, k)
         !--- Pairing Energy
         call Kitaichi(nd1, cvt_r2, Hv, z1)
         write(6, *) "Pair. Energy =", real(z1), aimag(z1)
         k = it+80000
         cvt_d2 = cvt_d2/sqrt(abs(zod)) !normalization
         call FIG_DENSITY(nd1, cvt_d2, pket, k)
         write(6,*) "VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV"
         write(6,*) ""
        ENDIF
      ENDIF
!???????????????????????????????????????????????????????????????????
   35103 end do

!--- Width "gamma(t)" from psurv.
      psdum = 7.d0 ; dsdum = 0.d0
      do it = -nctmax, nctmax, 1
         psdum(it+1) = psurv(it)
      end do
      call deriv_5p(psdum, dsdum, dct, 2*nctmax+1)
      open(1111, file='Rectemp_0031_Pserv_Gamma.dat', action='write')
      do it = -nctmax, nctmax, 1
         tc = dble(it)*dct ; gamma = -hbarc*dsdum(it+1)/psurv(it)
         write (1111, '(F18.5,2(2X,E15.7))') tc, psurv(it), gamma
      end do
      close(1111)
!--- Width from pndcy. Check consistency with Psurv.
      psdum = 0.d0 ; dsdum = 0.d0
      do it = 0, nctmax
         psdum(it+1) = pndcy(it,1)
      end do
      call deriv_5p(psdum, dsdum, dct, nctmax+1)
      open(9101, file='Rectemp_0032_Pndcy_Gamma.dat', action='write')
      do it = 0, nctmax
         tc = dble(it)*dct ; gamma = hbarc*dsdum(it+1)/(1.d0 - psdum(it+1))
         write (9101, '(F18.5,2(2X,E15.7))') tc, psdum(it+1), gamma
      end do
      close(9101)
!---
      psdum = 0.d0 ; dsdum = 0.d0
      do it = 0, nctmax
         psdum(it+1) = pndcy(it,2)
      end do
      call deriv_5p(psdum, dsdum, dct, nctmax+1)
      open(1112, file='Rectemp_0032_Pndcy_Gamma_02.dat', action='write')
      do it = 0, nctmax
         tc = dble(it)*dct ; gamma = hbarc*dsdum(it+1)/(1.d0 - psdum(it+1))
         write (1112, '(F18.5,2(2X,E15.7))') tc, psdum(it+1), gamma
      end do
      close(1112)
!---
      psdum = 0.d0 ; dsdum = 0.d0
      do it = 0, nctmax
         psdum(it+1) = pndcy(it,3)
      end do
      call deriv_5p(psdum, dsdum, dct, nctmax+1)
      open(1113, file='Rectemp_0032_Pndcy_Gamma_03.dat', action='write')
      do it = 0, nctmax
         tc = dble(it)*dct ; gamma = hbarc*dsdum(it+1)/(1.d0 - psdum(it+1))
         write (1113, '(F18.5,2(2X,E15.7))') tc, psdum(it+1), gamma
      end do
      close(1113)

     80008 deallocate(psurv, pndcy, psdum, dsdum)
     80009 deallocate(cvt_r2,cvt_r3,cvt_d2,cvt_d3)
   end subroutine


!*****************************************************************
   subroutine kitaichi(N, cv, h1, z1)
      implicit none
      integer, intent(in) :: N
      complex (kind=kind(0d0)), intent(in) :: cv(N)
      double precision, intent(in), dimension(N, N) :: h1
      complex (kind=kind(0d0)), intent(out) :: z1
      integer :: i, j
      complex (kind=kind(0d0)) :: zs
      zs = (0.d0, 0.d0)
      do j = 1, N ; do i = 1, N
           zs = zs + conjg(cv(i))*h1(i, j)*cv(j)
      end do ; end do
      z1 = zs
      return
   end subroutine kitaichi

!=====================================================================
   !--- The first derivative of f(x) (5 point formula). See Abramowicz, p.914.
   subroutine deriv_5p(f, df, dr, nr)
      implicit none
      integer, intent(in) :: nr
      double precision, intent(in) :: dr, f(nr)
      double precision, intent(out) :: df(nr)
      integer :: ir
         df(1) = (-50.d0*f(1) + 96.d0*f(2) - 72.d0*f(3) + 32.d0*f(4) - 6.d0*f(5))/24.d0/dr
         df(2) = ( -6.d0*f(1) - 20.d0*f(2) + 36.d0*f(3) - 12.d0*f(4) + 2.d0*f(5))/24.d0/dr
      df(nr-1) = ( -2.d0*f(nr-4) + 12.d0*f(nr-3) - 36.d0*f(nr-2) + 20.d0*f(nr-1) + 6.d0*f(nr) )/24.d0/dr
        df(nr) = (  6.d0*f(nr-4) - 32.d0*f(nr-3) + 72.d0*f(nr-2) - 96.d0*f(nr-1) + 50.d0*f(nr))/24.d0/dr
      do ir = 3, nr - 2
         df(ir) = f(ir-2) - 8.d0*f(ir-1) + 8.d0*f(ir+1) - f(ir+2)
         df(ir) = df(ir)/12.d0/dr
      end do
      return
   end subroutine deriv_5p
!=====================================================================

!********************************************************************
   subroutine hamatrix_pairless
      use mdl_setting, only : nconfmax
      use mdl_routines, only : mdiag
      use mdl_objects, only : confbasis
      implicit none
      integer :: nconf, nncf1(nconfmax), nncf2(nconfmax)
      double precision :: ecf12(nconfmax), fac12(nconfmax)
      integer :: i, j
      double precision, dimension(:), allocatable :: eall
      double precision, dimension(:,:), allocatable :: h, h0, ve
      call confbasis('load', nconf, nncf1, nncf2, ecf12, fac12)
      allocate(eall(nconf))
      allocate(h(nconf,nconf), h0(nconf,nconf), ve(nconf,nconf))
      h = 0.d0 ; h0 =0.d0 ; ve = 0.d0
      write(6,*)'Now calculating matrix elements...'
      do i = 1, nconf
         do j = i, nconf
            if (j.eq.i) h0(j,i) = ecf12(j)
            h0(i,j) = h0(j,i)
         enddo
      enddo
      h = h0
      call mdiag(h, eall, ve, nconf, nconf) !v_gs(i) = ve(i, 1)
      write(6,*) 'Hamiltonian has been diagonalized.'
      do i = 1, min(10,nconf)
         write (6, '(2X,A6,I2,A,F15.10,A7)') ' e_ex_', i - 1, '=', eall(i), '  (MeV)'
      end do
      write (6, *) ''
      deallocate(h,h0)
   end subroutine hamatrix_pairless
!********************************************************************
   SUBROUTINE HONTAI
      use mdl_setting, only : ecut, nspmax, RMAX, dr, nrmax, nconfmax, nss, npp
      use mdl_routines, only : vcn, vcn_nucl, vcn_coul, sp_numerov, cnfbss, ddpair, vmod
      use mdl_objects, only : spbasis, confbasis
      implicit none
 94   format(3(4x,i2),3x,f16.7)
 97   format(3i5,2f15.5)
 98   format(f5.0,4e15.5)
 99   format(7f10.5)
      integer :: ix
      double precision :: x
      write(6,*) ""
      write(6,*) "-----------------------------------"
      write(6,"(' Ecut = ',F7.3,' MeV')") Ecut
      write(6,"(' Rmax = ',F10.4,' fm')") RMAX
      write(6,"(' dr = ',F10.4,' fm')") dr
      write(6,"(' Nss = ',I2)") nss
      write(6,"(' Npp = ',I2)") npp
      write(6,*) "-----------------------------------"
      write(6,*) ""
!--- Write core-p potential.
      open(101, file='Rec_101_Vpot_CN.dat', action='write')
      do ix = -nrmax, nrmax
         x = dble(ix)*dr
         write(101,99) x, vcn(x), vcn_nucl(x), vcn_coul(x), vmod(x), ddpair(x)
      enddo
      close(101)
!--- Solve single-particle states.
      call sp_numerov ; call result_sp
!--- Configurate two-particle basis without pairing correlation.
      call cnfbss ; call result_conf

!--- Compute & diagonalize the Hamiltonian.
      write(6,*) ""
      write(6,*) "------------------------> Now hamatrix_tempo starts."
!!      call hamatrix_pairless
      call hamatrix_tempo

      write(6,*) "This program finished successfully (^o^) !!"
   END SUBROUTINE HONTAI
!********************************************************************
   program MAIN_clock
   implicit none
   integer :: it1, it2, it_rate, it_max, idiff
   double precision :: c1, c2
   call cpu_time( c1 ) ; call system_clock(it1)
   call HONTAI
   call system_clock(it2, it_rate, it_max)
   if ( it2 < it1 ) then
     idiff = (it_max - it1) + it2 + 1
   else
     idiff = it2 - it1
   endif
   call cpu_time( c2 )
   write(6,   *) "____________________________"
   write(6,   *) "system time :", idiff/it_rate, "seconds."
   write(6,   *) "cpu time    :", c2-c1, "seconds."
   open( 129, file='MAIN_clock.dat', action='write')
      write(129, *) "____________________________"
      write(129, *) "system time :", idiff/it_rate, "seconds."
      write(129, *) "cpu time    :", c2-c1, "seconds."
   close(129)
   end program MAIN_clock

