!#######################################################################
!      module usersrf
!
! To interface functions with QUANTICS edit this code and install
! as the libusrf.so library. 
!
! Multidimensional functions are implemented in udefsrf / uvpoint /
! usersurfinfo
! 1D functions are implemented in ufdef1d and ufunc1d

! A multidimensional function needs to be
! 1. defined in udefsrf by a string that can be specified in the 
!    operator definition.
! 2. called in uvpoint
! 3. referenced in usersurfinfo to write details to the log file.
!
! FOR EXAMPLES SEE defsrf, vpoint and surfinfo in opfuncs/funcsrf.F90
!
! A 1D function needs to be
! 1. defined in ufdef1d by a string that can be specified in the 
!    operator definition.
! 2. called in ufunc1d
!
! See routines below and QUANTICS documentation for details
! 
!#######################################################################

!#######################################################################
!    UDEFSRF    User-DEFinition-of-SuRFace  
!
! Multi-dimensional surfaces defined by the user.
!
! Main Input:
! label : string defining a function as given in specifiation of the 
!         Hamiltonian in the HAMILTONIAN_SECTION of the .op file
! Main Output:
! hopilab: integer associated with the function "label"
! ifile: integer defining the subroutine uvpoint 
!
! hopilab,ifile enable Quantics to find the function in the library
!
!
! interface with quantics also provides the following information
! which may be used to select options / parameters:
! ndof: no. of degrees of freedom in the system. 
! pesopts
! zpesopts
! maxhpar
! 
! hopipar
! hoprpar
!
! pesmin, pesmax: minimum / maximum values of PE function
! lcutpes : if true cut the PE function between pesmin, pesmax
!
! ierr: allows an error to be returned. If ierr > 0 Quantics will stop.
! ilog: unit to write information to the log file
! 
!#######################################################################

      subroutine udefsrf(hopilab,label,ifile,pesopts,zpesopts,hopipar, &
                 hoprpar,maxhpar,pesmin,pesmax,lcutpes,ndof,ierr,ilog)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4)

      integer(long)                                     :: ilbl,jlbl,ll
      integer(long), intent(in)                         :: maxhpar,ndof
      integer(long), intent(in)                         :: ilog
      integer(long), intent(out)                        :: hopilab,ifile,ierr
      integer(long), dimension(*), intent(out)          :: hopipar
      integer(long), dimension(2*maxhpar,*), intent(in) :: zpesopts
      real(dop), intent(in)                             :: pesmin,pesmax
      real(dop), dimension(*), intent(in)               :: hoprpar
      character(len=*), intent(in)                      :: label
      character(len=*), dimension(*), intent(in)        :: pesopts
      logical(kind=4), intent(in)                       :: lcutpes

      hopilab=0
      ifile=0
      ierr=0

!-----------------------------------------------------------------------
! Determine label length
!-----------------------------------------------------------------------
      ilbl=index(label,'{')-1
      jlbl=index(label,' ')-1
      if(ilbl.le.0) ilbl=32
      if(jlbl.le.0) jlbl=32
      ll = min(ilbl,jlbl)

!-----------------------------------------------------------------------
! MYSURF surface.
! Change MYSURF to a descriptive name to be used in the Quantics 
! definition of the Hamiltonian
! hopilab is the integer used by the uvpoint routine below to
! call the function
!-----------------------------------------------------------------------
      if (label(1:ll) .eq. 'mysrf') then
         hopilab=1

!-----------------------------------------------------------------------
! More functions can be defined with unique labels, incrementing 
! hopilab each time. 
!-----------------------------------------------------------------------
!      else if (label(1:ll) .eq. 'mysrf1') then
!         hopilab=2

!-----------------------------------------------------------------------
! end of if loop
!-----------------------------------------------------------------------
      endif

!-----------------------------------------------------------------------
! if label is assigned, assign ifile number.
!-----------------------------------------------------------------------
      if (hopilab .ne. 0) ifile=203

      return
      end subroutine udefsrf

!##########################################################################
!    UVPOINT   User-V-POINT
!
! This routine calculates the value of a potential at a grid point.
!
! To add new surfaces, add the new ifunc number, and then any code
! needed to manipulate the grid point and call the potential surface
! routine. Note that mysrf and DEF_mysrf already exist on source/surfaces, 
! the Fortran text of mysrf, however, has to be adapted to the needs.

! Input:
! ndof: no. of degrees of freedom
! gpoint(ndof): The point on the QUANTICS grid in au. The order is 
!               as specified in the primitive-basis-section of the 
!               input file. It is not checked so must be in the 
!               order required by the potential function.
! hopilab: integer locating function call in select list. This is 
!          as defined by the label in udefsrf above.
! hopipar: integer parameters associated with the function, defined in
!          udefsrf above.
! hoprpar: real parameters associated with the function, defined in
!          udefsrf above.
!
! Output:
! v: value of potential at the point gpoint. 
!    NB: v must be returned in au
!
! Other input:
! ilog: unit of Quantics log file if anything needs to be written
!##########################################################################

      subroutine uvpoint (v,gpoint,ndof,hopilab,hopipar,hoprpar,ilog)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4)

      integer(long), intent(in)               :: ndof
      integer(long), intent(in)               :: hopilab
      integer(long), intent(in)               :: ilog
      integer(long), dimension(*), intent(in) :: hopipar
      real(dop), dimension(ndof), intent(in)  :: gpoint
      real(dop), dimension(*), intent(in)     :: hoprpar
      real(dop), intent(out)                  :: v

     select case (hopilab)

     case default
      write(ilog,'(a,i5)') 'ERROR: User function not found :',hopilab
      stop

!-----------------------------------------------------------------------
! edit to add call to the function related to label MYSURF
! and comment out (or remove) the error statement + stop
!-----------------------------------------------------------------------
     case (1) 
!     call mysrf(v,gpoint)
      write(6,'(a)') 'ERROR: mysrf has not been defined'    ! remove
      write(ilog,'(a)') 'ERROR: mysrf has not been defined' ! remove
      stop                                                  ! remove
      return

!-----------------------------------------------------------------------
! edit to add call to the a second function related to label MYSURF1
!-----------------------------------------------------------------------
!     case (2) ! EXAMPLE
!     call mysrf1(v,gpoint)
!      return

     end select

      end subroutine uvpoint

!#######################################################################
!  User Surface Info. 
!  Description of the user implemented surface and their coordinates.
!
! Input:
! hopilab: integer defining the function
! hopipar: integer parameters for the function defined in input
! hoprpar: real parameters for the function defined in input
! mlabel: labels for DOFs 
! ilog: unit for log file
!#######################################################################

      subroutine usersurfinfo(hopilab,hopipar,hoprpar,mlabel,nn,ilog)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4)

      integer(long), intent(in)                   :: hopilab,ilog,nn
      integer(long), dimension(*), intent(in)     :: hopipar
      real(dop), dimension(*), intent(in)         :: hoprpar
      character(len=16), dimension(*), intent(in) :: mlabel

     select case (hopilab)    ! Do not forget to update this goto.

     case default
      write(ilog,'(/,a)') 'Error in: usersurfinfo (opfuncs/usersrf.F)'
      write(ilog,'(a,i5)') 'Surface not found :', hopilab

     case (1)
! EDIT FOLLOWING TEXT RELATED TO MYSRF
      write(ilog,'(a)') 'EXAMPLE PES Change this text accordingly'
      write(ilog,'(a)') 'to Define PE Surface in log file.'
      write(ilog,'(a)') 'This allows a check that correct function has'
      write(ilog,'(a)') 'been chosen and coordinates are in correct order.'
      write(ilog,'(a)') 'Can also output hopipar and hoprpar to check '
      write(ilog,'(a)') 'that parameters are correctly defined '

      write(ilog,'(a)') 'Mysrf Coordinates '
      write(ilog,'(2a)') 'x    : ',mlabel(1)
      write(ilog,'(2a)') 'y    : ',mlabel(2)
      write(ilog,'(2a)') 'z    : ',mlabel(3)

!     case (2)
! EDIT FOLLOWING TEXT RELATED TO MYSURF1
!      write(ilog,'(a)') 'EXAMPLE PES 1'

!     case (3)
! EDIT FOLLOWING TEXT RELATED TO MYSURF2
!      write(ilog,'(a)') 'EXAMPLE PES 2'

     end select

      return
      end subroutine usersurfinfo

!      end module usersrf

!#######################################################################
!     UFDEF1D       User-Function-DEFinition-1D
!
! 1D CURVES  defined by the user.
! See  opfuncs/func1d.F for examples.
!
!#######################################################################

      subroutine ufdef1d(ifunc,label,ifile,hopipar,hoprpar, &
                         maxpar,hoppar)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to
                                                          !  integer*4)
      integer(long)                                :: ilbl,jlbl,np
      integer(long), intent(in)                    :: maxpar
      integer(long), intent(out)                   :: ifunc,ifile
      integer(long), dimension(maxpar), intent(in) :: hopipar
      real(dop), dimension(maxpar), intent(out)    :: hoprpar
      real(dop), dimension(maxpar,*), intent(in)   :: hoppar
      character(len=64), intent(in)                :: label

!-----------------------------------------------------------------------
! get argument list no.
!-----------------------------------------------------------------------
      np=0
      ilbl=index(label,'[')+1
      if (ilbl .gt. 1) then
         jlbl = index(label,']')-1
         read(label(ilbl:jlbl),*) np
      endif

!-----------------------------------------------------------------------
! get end of operator label
!-----------------------------------------------------------------------
      ilbl=index(label,'[')-1
      if (ilbl .lt. 0) ilbl=index(label,' ')-1


!-----------------------------------------------------------------------
! now assign function number
! Note: hopipar(1) and hoprpar(1) are reserved for the exponent.
! The exponent was read in defanld (opfuncs/callanld.f).
!-----------------------------------------------------------------------

! my1d: my1d^r  (provide subroutine on file: user1d.F)
      if (label(1:ilbl) .eq. 'my1d') then
         if (np .gt.0) then
            hoprpar(2)=hoppar(1,np)
            hoprpar(3)=hoppar(2,np)
            hoprpar(4)=hoppar(3,np)
            hoprpar(5)=hoppar(4,np)
            hoprpar(6)=hoppar(5,np)
         else
            hoprpar(2)=0.0_dop
            hoprpar(3)=0.0_dop
            hoprpar(4)=0.0_dop
            hoprpar(5)=0.0_dop
            hoprpar(6)=0.0_dop
         endif
         ifunc=1

! newfunc     ! This is, of course, just an example.
!     elseif (label(1:ilbl) .eq. 'newfunc') then
!        ifunc=2


!-----------------------------------------------------------------------
! end of if loop
!-----------------------------------------------------------------------
      endif

!-----------------------------------------------------------------------
! if label is assigned, assign ifile no.
!-----------------------------------------------------------------------
      if (ifunc .ne. 0) ifile=306

      return
      end subroutine ufdef1d

!#######################################################################
!  UFUNC1D    User-Functions-1D
!
! Define here new user specific, real 1D functions.
! See funcanld.F for examples.
!#######################################################################

   subroutine ufunc1d(v,x,ifunc,hopipar,hoprpar,ilog)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to
                                                          !  integer*4)
      integer(long), intent(in)               :: ifunc,ilog
      integer(long), dimension(*), intent(in) :: hopipar
      real(dop), dimension(*), intent(in)     :: hoprpar
      real(dop), intent(in)                   :: x
      real(dop), intent(out)                  :: v


! Decide which function to call, accroding to ifunc

     select case (ifunc)    ! Do not forget to enlarge the goto list.

     case default
      write(ilog,'(/,a)') 'Error in ufunc1d (opfuncs/user1D.F)'
      write(ilog,'(a,i5)') 'Error: Function not found :',ifunc

! New functions to be defined by the user.
!     Set here the routines you want to call , or provide directly
!     the appropriate fortran text. (Compare with funcanld.f)

!     case (1)! my1d:  1d-routine provided by the user.
!       call my1d(x,v,hoprpar(2),hoprpar(3),hoprpar(4), &
!                     hoprpar(5),hoprpar(6))

!    case(2)
!     call newfunc(x,v) 

! Last Function Defined
     end select

! Apply exponent.
      if(hopipar(1).eq.0) then
         v=v**hoprpar(1)
      elseif(hopipar(1).ne.1) then
         v=v**hopipar(1)
      endif
      return


   end subroutine ufunc1d

!#######################################################################
