module bmmodel

  implicit none

  ! PARAMETERS
  integer, parameter :: nlambda = 9

contains

  ! ******************************************************************
  ! ******************************************************************

  subroutine checkinglambda(nlambda,lambda,DayE)

    use bmproblem, only: nobs,obs,obsx,obst,obsz,qboundx,bedel,width
  
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: DayE,nlambda
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! PARAMETERS
    integer, parameter :: npoints = 10000

    ! LOCAL SCALARS
    integer :: i
    real(kind=8) :: q,t,x,xmax,xmin,w,z,zb,zpred

    ! LOCAL ARRAYS
    real(kind=8) :: coe(nlambda)

    write(*,*)
    write(*,*) 'Optimal coeffcients: '
    do i = 1,nlambda
       write(*,*) 'c(',i,') = ',lambda(i)
    end do
    
    write(*,*)
    do i = 1,nobs
       if ( obst(i) .eq. DayE * 24.0d0 * 3600.0d0 ) then
          call coeff(nlambda,coe,obsx(i),obst(i))
          zpred = dot_product( coe(1:nlambda), lambda(1:nlambda) ) 
          write(*,*)
          write(*,*) 'z_observed = ',obsz(i),' z_predicted = ',zpred
          write(*,*) 'absolute = z_predicted - z_observed = ',zpred - obsz(i)
          write(*,*) 'relative = ( z_predicted - z_observed ) / z_observed = ',( zpred - obsz(i) ) / obsz(i)
       end if
    end do

    xmin = min( qboundx, minval( obsx(1:nobs), obs(1:nobs) ) )
    xmax = maxval( obsx(1:nobs), obs(1:nobs) )
    
    t = DayE * 24.0d0 * 3600.0d0

    open(unit=10,file='allx-predictions.txt')
    do i = 1,npoints
       x = xmin + ( xmax - xmin ) * ( i - 1 ) / ( npoints - 1 )
       call coeff(nlambda,coe,x,t)
       z = dot_product( coe(1:nlambda), lambda(1:nlambda) )
       call width(x,w)
       call bedel(x,zb)
       call psi(nlambda,lambda,x,t,q)
       write(10,*) x,z,z - zb,( z - zb ) * w,q
    end do
    close(10)
    
    open(unit=10,file='allx-observations.txt')
    do i = 1,nobs
       if ( obst(i) .eq. t ) then
          call width(obsx(i),w)
          call bedel(obsx(i),zb)
          write(10,*) obsx(i),obsz(i),obsz(i) - zb,( obsz(i) - zb ) * w
       end if
    end do
    close(10)

    write(*,*)
    write(*,*) 'File allx-predictions.txt containing predicitions for z, h, A, and Q ', &
         'for t=DayE days and all x in [xmin,xmax] was written.'

    write(*,*)
    write(*,*) 'File allx-observations.txt containing available test data (i.e. not used ', &
         'in the predicitions) related to z, h, A, and Q for t=DayE days and some x in [xmin,xmax] ', &
         'was written. The file can be used to plot predicted versus observed graphics.'

  end subroutine checkinglambda
  
! ******************************************************************
! ******************************************************************

  subroutine complambdaalgencan(nlambda,lambda,nattempt)

    use bmproblem, only: nsv
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nattempt,nlambda

    ! ARRAY ARGUMENTS
    real(kind=8), intent(inout) :: lambda(nlambda)

    ! LOCAL SCALARS
    logical :: checkder
    integer :: allocerr,hnnzmax,inform,jcnnzmax,m,n,nvparam
    real(kind=8) :: cnorm,efacc,efstain,eoacc,eostain,epsfeas,epsopt,f,nlpsupn,snorm

    ! LOCAL ARRAYS
    character(len=80) :: specfnm,outputfnm,vparam(10)
    logical :: coded(11)
    logical, pointer :: equatn(:),linear(:)
    real(kind=8), pointer :: l(:),lagmul(:),u(:),x(:),c(:)

    n = nlambda

    allocate(x(n),l(n),u(n),stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'Allocation error in complambdaalgencan.'
       stop
    end if

    l(1:n) = - huge( 1.0d0 )
    u(1:n) =   huge( 1.0d0 )

    x(1:n) = lambda(1:nlambda)

    if ( nattempt .eq. 1 ) then
       m = 3 * nsv
    else
       m = 4 * nsv
    end if
    
    allocate(equatn(m),linear(m),lagmul(m),c(m),stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'Allocation error in complambdaalgencan.'
       stop
    end if

    equatn(1:m) = .false.
    linear(1:m) = .false.
    lagmul(1:m) = 0.0d0

    coded(1:11) = .false.
    coded(7) = .true. ! fcsub
    coded(8) = .true. ! gjacsub

    jcnnzmax  = m * nlambda
    hnnzmax = 0

    checkder = .false.

    epsfeas = 1.0d-08

    if ( nattempt .eq. 1 ) then
       epsopt  = 1.0d-08
    else
       epsopt = 0.5d0
    end if
    
    efstain = sqrt( epsfeas )
    eostain = 1.0d-12

    efacc = - 1.0d0
    eoacc = - 1.0d0

    outputfnm = ''
    specfnm = ''
    nvparam = 1
    vparam(1) = 'OBJECTIVE-AND-CONSTRAINTS-SCALING-AVOIDED'

    if ( nattempt .eq. 1 ) then
       call algencan(myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
            myevalfc,myevalgjac,myevalgjacp,myevalhl,myevalhlp,jcnnzmax, &
            hnnzmax,epsfeas,epsopt,efstain,eostain,efacc,eoacc,outputfnm, &
            specfnm,nvparam,vparam,n,x,l,u,m,lagmul,equatn,linear,coded, &
            checkder,f,cnorm,snorm,nlpsupn,inform)
    else
       call algencan(myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
            myevalfc2,myevalgjac2,myevalgjacp,myevalhl,myevalhlp,jcnnzmax, &
            hnnzmax,epsfeas,epsopt,efstain,eostain,efacc,eoacc,outputfnm, &
            specfnm,nvparam,vparam,n,x,l,u,m,lagmul,equatn,linear,coded, &
            checkder,f,cnorm,snorm,nlpsupn,inform)
    end if
    
    lambda(1:nlambda) = x(1:n)
    
    deallocate(x,l,u,lagmul,equatn,linear,c,stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'Allocation error in complambdaalgencan.'
       stop
    end if

    contains
  
    ! ******************************************************************
    subroutine myevalf(n,x,f,flag)
      implicit none
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: f
      real(kind=8), intent(in) :: x(n)
      if ( .false. ) write(*,*) n,x(1:n),f,flag
      flag = - 1
    end subroutine myevalf

    ! ******************************************************************
    subroutine myevalg(n,x,g,flag)
      implicit none
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: g(n)
      if ( .false. ) write(*,*) n,x(1:n),g(1:n),flag
      flag = - 1
    end subroutine myevalg

    ! ******************************************************************
    subroutine myevalh(n,x,hrow,hcol,hval,hnnz,lim,lmem,flag)
      implicit none
      logical, intent(out) :: lmem
      integer, intent(in) :: lim,n
      integer, intent(out) :: flag,hnnz
      integer, intent(out) :: hcol(lim),hrow(lim)
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: hval(lim)
      if ( .false. ) write(*,*) n,x(1:n),hrow(1:lim),hcol(1:lim),hval(1:lim),hnnz,lim,lmem,flag
      flag = - 1
    end subroutine myevalh

    ! ******************************************************************
    subroutine myevalc(n,x,ind,c,flag)
      implicit none
      integer, intent(in) :: ind,n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: c
      real(kind=8), intent(in) :: x(n)
      if ( .false. ) write(*,*) n,x(1:n),ind,c,flag
      flag = - 1
    end subroutine myevalc

    ! ******************************************************************
    subroutine myevaljac(n,x,ind,jcvar,jcval,jcnnz,lim,lmem,flag)
      implicit none
      logical, intent(out) :: lmem
      integer, intent(in) :: ind,lim,n
      integer, intent(out) :: flag,jcnnz
      integer, intent(out) :: jcvar(lim)
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: jcval(lim)
      if ( .false. ) write(*,*) n,x(1:n),ind,jcvar(1:lim),jcval(1:lim),jcnnz,lim,lmem,flag
      flag = - 1
    end subroutine myevaljac

    ! ******************************************************************
    subroutine myevalhc(n,x,ind,hcrow,hccol,hcval,hcnnz,lim,lmem,flag)
      implicit none
      logical, intent(out) :: lmem
      integer, intent(in) :: ind,lim,n
      integer, intent(out) :: flag,hcnnz
      integer, intent(out) :: hccol(lim),hcrow(lim)
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: hcval(lim)
      if ( .false. ) write(*,*) n,x(1:n),ind,hcrow(1:lim),hccol(1:lim),hcval(1:lim),hcnnz,lim,lmem,flag
      flag = - 1
    end subroutine myevalhc

    ! ******************************************************************
    subroutine myevalgjacp(n,x,g,m,p,q,work,gotj,flag)
      implicit none
      logical, intent(inout) :: gotj
      integer, intent(in) :: m,n
      integer, intent(out) :: flag
      character, intent(in) :: work
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(inout) :: p(m),q(n)
      real(kind=8), intent(out) :: g(n)
      if ( .false. ) write(*,*) n,x(1:n),g(1:n),m,p(m),q(n),work,gotj
      flag = - 1
    end subroutine myevalgjacp
    
    ! ******************************************************************
    subroutine myevalhl(n,x,m,lambda,sf,sc,hlrow,hlcol,hlval,hlnnz,lim,lmem,flag)
      implicit none
      logical, intent(out) :: lmem
      integer, intent(in) :: lim,m,n
      integer, intent(out) :: flag,hlnnz
      real(kind=8), intent(in)  :: sf
      integer, intent(out) :: hlcol(lim),hlrow(lim)
      real(kind=8), intent(in) :: lambda(m),sc(m),x(n)
      real(kind=8), intent(out) :: hlval(lim)
      if ( .false. ) write(*,*) n,x(1:n),m,lambda(1:m),sf,sc(1:m),hlrow(1:lim),hlcol(1:lim),hlval(1:lim),hlnnz,lim,lmem,flag
      flag = - 1
    end subroutine myevalhl

    ! ******************************************************************
    subroutine myevalhlp(n,x,m,lambda,sf,sc,p,hp,goth,flag)
      implicit none
      logical, intent(inout) :: goth
      integer, intent(in) :: m,n
      integer, intent(out) :: flag
      real(kind=8), intent(in) :: sf
      real(kind=8), intent(in) :: lambda(m),p(n),sc(m),x(n)
      real(kind=8), intent(out) :: hp(n)
      if ( .false. ) write(*,*) n,x(1:n),m,lambda(1:m),sf,sc(1:m),p(1:n),hp(1:n),goth,flag
      flag = - 1
    end subroutine myevalhlp
    
  end subroutine complambdaalgencan
  
! ******************************************************************
! ******************************************************************

  subroutine myevalfc(n,x,f,m,c,flag)

    use bmproblem, only: nobs,obs,obsx,obst,obsz,bedel,dbedel,dwidth,width,nsv,svx,svt
  
    implicit none
      
    ! SCALAR ARGUMENTS
    integer, intent(in) :: m,n
    integer, intent(out) :: flag
    real(kind=8), intent(out) :: f

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: c(m)
    
    ! PARAMETERS
    real(kind=8), parameter :: grav = 9.81d0

    ! LOCAL SCALARS
    integer :: i
    real(kind=8) :: a,dadq,q,w,zb

    ! LOCAL ARRAYS
    real(kind=8) :: lambda(nlambda),coe(nlambda)

    flag = 0

    lambda(1:nlambda) = x(1:nlambda)
  
    ! Objective function

    f = 0.0d0
    do i = 1,nobs
       if ( obs(i) ) then
          call width(obsx(i),w)
          call bedel(obsx(i),zb)
          call phi(nlambda,lambda,obsx(i),obst(i),a)
          f = f + ( a - ( obsz(i) - zb ) * w ) ** 2
       end if
    end do
      
    if ( count( obs(1:nobs) ) .gt. 0 ) then
       f = f / count( obs(1:nobs) )
    end if

    ! Constraints

    do i = 1,nsv
       call bedel(svx(i),zb)
       call coeff(nlambda,coe,svx(i),svt(i))
       c(3*i-2) = zb - dot_product( coe(1:nlambda), lambda(1:nlambda) )
       
       call psi(nlambda,lambda,svx(i),svt(i),q)
       c(3*i-1) = - q
     
       call dphidq(nlambda,lambda,svx(i),svt(i),dadq)
       c(3*i) = - dadq
    end do
    
  end subroutine myevalfc

! ******************************************************************
! ******************************************************************

  subroutine myevalfc2(n,x,f,m,c,flag)

    use bmproblem, only: nobs,obs,obsx,obst,obsz,bedel,dbedel,dwidth,width,nsv,svx,svt
  
    implicit none
      
    ! SCALAR ARGUMENTS
    integer, intent(in) :: m,n
    integer, intent(out) :: flag
    real(kind=8), intent(out) :: f

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: c(m)
    
    ! PARAMETERS
    real(kind=8), parameter :: grav = 9.81d0

    ! LOCAL SCALARS
    integer :: i
    real(kind=8) :: a,dadq,dadx,dqdt,dqdx,dzdx,q,sv2a,sv2b,sv2c,w,zb

    ! LOCAL ARRAYS
    real(kind=8) :: coe(nlambda),dcoedx(nlambda),lambda(nlambda)

    flag = 0

    lambda(1:nlambda) = x(1:nlambda)
  
    ! Objective function

    f = 0.0d0
    do i = 1,nobs
       if ( obs(i) ) then
          call width(obsx(i),w)
          call bedel(obsx(i),zb)
          call phi(nlambda,lambda,obsx(i),obst(i),a)
          f = f + ( a - ( obsz(i) - zb ) * w ) ** 2
       end if
    end do
      
    if ( count( obs(1:nobs) ) .gt. 0 ) then
       f = f / count( obs(1:nobs) )
    end if

    ! Constraints

    do i = 1,nsv
       call bedel(svx(i),zb)
       call coeff(nlambda,coe,svx(i),svt(i))
       c(4*i-3) = zb - dot_product( coe(1:nlambda), lambda(1:nlambda) )
       
       call psi(nlambda,lambda,svx(i),svt(i),q)
       c(4*i-2) = - q
     
       call dphidq(nlambda,lambda,svx(i),svt(i),dadq)
       c(4*i-1) = - dadq

       call dpsidt(nlambda,lambda,svx(i),svt(i),dqdt)

       sv2a = dqdt

       call phi(nlambda,lambda,svx(i),svt(i),a)
       call dphidx(nlambda,lambda,svx(i),svt(i),dadx)
       call dpsidx(nlambda,lambda,svx(i),svt(i),dqdx)

       sv2b = ( 2.0d0 * dqdx - ( q / a ) * dadx ) * ( q / a )

       call dcoeffdx(nlambda,dcoedx,svx(i),svt(i))
       dzdx = dot_product( dcoedx(1:nlambda), lambda(1:nlambda) )
      
       sv2c = grav * a * dzdx
     
       c(4*i) = sv2a + sv2b + sv2c + 1.0d-08
!!$       write(*,*) i,sv2a,sv2b,sv2c,c(4*i),q,dqdx,a,dadx
    end do
    
  end subroutine myevalfc2

! ******************************************************************
! ******************************************************************

  subroutine myevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag)

    use bmproblem, only: nobs,obs,obsx,obst,obsz,bedel,dbedel,dwidth,width,nsv,svx,svt
  
    implicit none

    ! SCALAR ARGUMENTS
    logical, intent(out) :: lmem
    integer, intent(in) :: lim,m,n
    integer, intent(out) :: flag,jcnnz

    ! ARRAY ARGUMENTS
    integer, intent(out) :: jcfun(lim),jcvar(lim)
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: g(n),jcval(lim)

    ! PARAMETERS
    real(kind=8), parameter :: grav = 9.81d0

    ! LOCAL SCALARS
    integer :: i,j
    real(kind=8) :: a,w,zb

    ! LOCAL ARRAYS
    real(kind=8) :: dadl(nlambda),lambda(nlambda),dqdl(nlambda),d2adql(nlambda),coe(nlambda)

    if ( .false. ) write(*,*) m
      
    flag = 0
    lmem = .false.

    lambda(1:nlambda) = x(1:nlambda)
  
    ! Objective function's gradient

    g(1:n) = 0.0d0
  
    do i = 1,nobs
       if ( obs(i) ) then
          call width(obsx(i),w)
          call bedel(obsx(i),zb)
          call phi(nlambda,lambda,obsx(i),obst(i),a)
          call dphidl(nlambda,lambda,obsx(i),obst(i),dadl)
          ! f = f + ( a - ( obsz(i) - zb ) * w ) ** 2
          g(1:nlambda) = g(1:nlambda) + 2.0d0 * ( a - ( obsz(i) - zb ) * w ) * dadl(1:nlambda)
       end if
    end do

    if ( count( obs(1:nobs) ) .gt. 0 ) then
       g(1:n) = g(1:n) / count( obs(1:nobs) )
    end if

    ! Jacobian of the constraints
      
    jcnnz = 0
  
    do i = 1,nsv
       call coeff(nlambda,coe,svx(i),svt(i))
       
       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 3 * i - 2
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - coe(1:nlambda)

       jcnnz = jcnnz + nlambda
       
       call dpsidl(nlambda,lambda,svx(i),svt(i),dqdl)

       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 3 * i - 1
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - dqdl(1:nlambda)

       jcnnz = jcnnz + nlambda
         
       call d2phidql(nlambda,lambda,svx(i),svt(i),d2adql)

       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 3 * i
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - d2adql(1:nlambda)

       jcnnz = jcnnz + nlambda
    end do
    
  end subroutine myevalgjac

! ******************************************************************
! ******************************************************************

  subroutine myevalgjac2(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag)

    use bmproblem, only: nobs,obs,obsx,obst,obsz,bedel,dbedel,dwidth,width,nsv,svx,svt
  
    implicit none

    ! SCALAR ARGUMENTS
    logical, intent(out) :: lmem
    integer, intent(in) :: lim,m,n
    integer, intent(out) :: flag,jcnnz

    ! ARRAY ARGUMENTS
    integer, intent(out) :: jcfun(lim),jcvar(lim)
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: g(n),jcval(lim)

    ! PARAMETERS
    real(kind=8), parameter :: grav = 9.81d0

    ! LOCAL SCALARS
    integer :: i,j
    real(kind=8) :: a,q,dadx,dqdx,dzdx,sv2b,sv2ba,sv2bb,w,zb

    ! LOCAL ARRAYS
    real(kind=8) :: dadl(nlambda),lambda(nlambda),dsv2a(nlambda),dsv2b(nlambda), &
         dsv2ba(nlambda),dsv2bb(nlambda),dsv2c(nlambda),dqdl(nlambda),d2adxl(nlambda), &
         d2qdxl(nlambda),d2adql(nlambda),coe(nlambda),dcoedx(nlambda)

    if ( .false. ) write(*,*) m
      
    flag = 0
    lmem = .false.

    lambda(1:nlambda) = x(1:nlambda)
  
    ! Objective function's gradient

    g(1:n) = 0.0d0
  
    do i = 1,nobs
       if ( obs(i) ) then
          call width(obsx(i),w)
          call bedel(obsx(i),zb)
          call phi(nlambda,lambda,obsx(i),obst(i),a)
          call dphidl(nlambda,lambda,obsx(i),obst(i),dadl)
          ! f = f + ( a - ( obsz(i) - zb ) * w ) ** 2
          g(1:nlambda) = g(1:nlambda) + 2.0d0 * ( a - ( obsz(i) - zb ) * w ) * dadl(1:nlambda)
       end if
    end do

    if ( count( obs(1:nobs) ) .gt. 0 ) then
       g(1:n) = g(1:n) / count( obs(1:nobs) )
    end if

    ! Jacobian of the constraints
      
    jcnnz = 0
  
    do i = 1,nsv
       call coeff(nlambda,coe,svx(i),svt(i))
       
       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 4 * i - 3
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - coe(1:nlambda)

       jcnnz = jcnnz + nlambda
       
       call dpsidl(nlambda,lambda,svx(i),svt(i),dqdl)

       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 4 * i - 2
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - dqdl(1:nlambda)

       jcnnz = jcnnz + nlambda
         
       call d2phidql(nlambda,lambda,svx(i),svt(i),d2adql)

       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 4 * i - 1
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = - d2adql(1:nlambda)

       jcnnz = jcnnz + nlambda
       
       ! c(3*i) = sv2a + sv2b + sv2c

       ! sv2a = dqdt
       call d2psidtl(nlambda,lambda,svx(i),svt(i),dsv2a)

       call phi(nlambda,lambda,svx(i),svt(i),a)
       call psi(nlambda,lambda,svx(i),svt(i),q)
       call dphidx(nlambda,lambda,svx(i),svt(i),dadx)
       call dpsidx(nlambda,lambda,svx(i),svt(i),dqdx)

       ! sv2bb = q / a
       ! sv2ba = 2.0d0 * dqdx - sv2bb * dadx
       ! sv2b = sv2ba * sv2bb

       call dphidl(nlambda,lambda,svx(i),svt(i),dadl)
       call dpsidl(nlambda,lambda,svx(i),svt(i),dqdl)
       call d2phidxl(nlambda,lambda,svx(i),svt(i),d2adxl)
       call d2psidxl(nlambda,lambda,svx(i),svt(i),d2qdxl)

       sv2bb = q / a
       dsv2bb(1:nlambda) = dqdl(1:nlambda) / a - q / a ** 2 * dadl(1:nlambda)

       sv2ba = 2.0d0 * dqdx - sv2bb * dadx
       dsv2ba(1:nlambda) = 2.0d0 * d2qdxl(1:nlambda) - dsv2bb(1:nlambda) * dadx - sv2bb * d2adxl(1:nlambda)

       sv2b = sv2ba * sv2bb
       dsv2b(1:nlambda) = dsv2ba(1:nlambda) * sv2bb + sv2ba * dsv2bb(1:nlambda)

!!$       dsv2ba(1:nlambda) = 2.0d0 * ( ( dqdl(1:nlambda) * dqdx + q * d2qdxl(1:nlambda) ) / a &
!!$            - q * dqdx / a ** 2 * dadl(1:nlambda) )
!!$         
!!$       dsv2bb(1:nlambda) = ( ( 2.0d0 * q * dqdl(1:nlambda) / a ** 2 - 2.0d0 * q ** 2 / a ** 3 * &
!!$            dadl(1:nlambda) ) * dadx + q ** 2 / a ** 2 * d2adxl(1:nlambda) )
!!$         
!!$       dsv2b(1:nlambda) = dsv2ba(1:nlambda) - dsv2bb(1:nlambda)
         
!!$       call width(svx(i),w)
!!$       call dwidth(svx(i),dwdx)
!!$       call dbedel(svx(i),dzbdx)
!!$       dzdx = dadx / w - a / w ** 2 * dwdx + dzbdx
!!$      
!!$       ! sv2c = grav * a * dzdx
!!$       dsv2c(1:nlambda) = grav * ( dadl(1:nlambda) * dzdx + a * ( d2adxl(1:nlambda) / w - dadl(1:nlambda) / w ** 2 * dwdx ) )

       call dcoeffdx(nlambda,dcoedx,svx(i),svt(i))
       dzdx = dot_product( dcoedx(1:nlambda), lambda(1:nlambda) )
      
       ! sv2c = grav * a * dzdx
       dsv2c(1:nlambda) = grav * ( dadl(1:nlambda) * dzdx + a * dcoedx(1:nlambda) )
       
       if ( jcnnz + nlambda .gt. lim ) then
          lmem = .true.
          return
       end if

       jcfun(jcnnz+1:jcnnz+nlambda) = 4 * i
       jcvar(jcnnz+1:jcnnz+nlambda) = (/ (j,j=1,nlambda) /)
       jcval(jcnnz+1:jcnnz+nlambda) = dsv2a(1:nlambda) + dsv2b(1:nlambda) + dsv2c(1:nlambda)

       jcnnz = jcnnz + nlambda
    end do
    
  end subroutine myevalgjac2

! ******************************************************************
! ******************************************************************

  subroutine coeff(nlambda,c,x,t)

    use bmproblem, only: qbound

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(out) :: c(nlambda)

    ! Phi(x,t) = P(x,Q(qboundx,t)) w(x), where w(x) is a piecewise
    ! polynomial of degree 3 in x and P(x,Q(qboundx,t)) is a
    ! polynomial of degree 2 in x and Q(qboundx,t). P(x,Q(qboundx,t))
    ! is of the form c(x,Q(qboundx,t))^T lambda.
    !
    ! This subroutine defines c(x,Q(qboundx,t)).
    
    ! LOCAL SCALARS
    real(kind=8) :: qb
    
    call qbound(t,qb)

    c(1) = 1.0d0
    c(2) = qb
    c(3) = qb ** 2
    c(4) = 1.0d-03 * x
    c(5) = 1.0d-03 * x * qb
    c(6) = 1.0d-03 * x * qb ** 2
    c(7) = 1.0d-06 * x ** 2
    c(8) = 1.0d-06 * x ** 2 * qb
    c(9) = 1.0d-06 * x ** 2 * qb ** 2

  end subroutine coeff
  
! ******************************************************************
! ******************************************************************

  subroutine dcoeffdx(nlambda,dcdx,x,t)

    use bmproblem, only: qbound

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(out) :: dcdx(nlambda)

    ! This subroutines defines the partial derivative of
    ! c(x,Q(qboundx,t)) with respect to its first parameter, where c
    ! is givenn by subroutine coeff.

    ! LOCAL SCALARS
    real(kind=8) :: qb
    
    call qbound(t,qb)

    dcdx(1) = 0.0d0
    dcdx(2) = 0.0d0
    dcdx(3) = 0.0d0
    dcdx(4) = 1.0d-03
    dcdx(5) = 1.0d-03 * qb
    dcdx(6) = 1.0d-03 * qb ** 2
    dcdx(7) = 1.0d-06 * 2.0d0 * x 
    dcdx(8) = 1.0d-06 * 2.0d0 * x * qb
    dcdx(9) = 1.0d-06 * 2.0d0 * x * qb ** 2

  end subroutine dcoeffdx
  
! ******************************************************************
! ******************************************************************

  subroutine dcoeffdqb(nlambda,dcdqb,x,t,withoutx)

    use bmproblem, only: qbound

    implicit none

    ! SCALAR ARGUMENTS
    logical, optional, intent(in) :: withoutx
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(out) :: dcdqb(nlambda)

    ! This subroutines defines the partial derivative of
    ! c(x,Q(qboundx,t)) with respect to its second parameter, where c
    ! is givenn by subroutine coeff.

    ! LOCAL SCALARS
    real(kind=8) :: qb
    
    call qbound(t,qb)
  
    if ( .not. present( withoutx ) ) then
       dcdqb(1) = 0.0d0
       dcdqb(2) = 1.0d0
       dcdqb(3) = 2.0d0 * qb
       dcdqb(4) = 0.0d0
       dcdqb(5) = 1.0d-03 * x
       dcdqb(6) = 1.0d-03 * x * 2.0d0 * qb
       dcdqb(7) = 0.0d0
       dcdqb(8) = 1.0d-06 * x ** 2
       dcdqb(9) = 1.0d-06 * x ** 2 * 2.0d0 * qb
       
    else
       if ( .not. withoutx ) then
          write(*,*) 'In dcoeffdqb, it is not expected argument withoutx to be present and be .false.'
          stop
       end if

       dcdqb(1) = 0.0d0
       dcdqb(2) = 1.0d0
       dcdqb(3) = 2.0d0 * qb
       dcdqb(4) = 0.0d0
       dcdqb(5) = 1.0d-03
       dcdqb(6) = 1.0d-03 * 2.0d0 * qb
       dcdqb(7) = 0.0d0
       dcdqb(8) = 1.0d-06
       dcdqb(9) = 1.0d-06 * 2.0d0 * qb
    end if
    
  end subroutine dcoeffdqb

! ******************************************************************
! ******************************************************************

  subroutine d2coeffdqbt(nlambda,d2cdqbt,t)

    use bmproblem, only: dqbound

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: t
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(out) :: d2cdqbt(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: dqb
    
    call dqbound(t,dqb)
  
    d2cdqbt(1) = 0.0d0
    d2cdqbt(2) = 0.0d0
    d2cdqbt(3) = 2.0d0 * dqb
    d2cdqbt(4) = 0.0d0
    d2cdqbt(5) = 0.0d0
    d2cdqbt(6) = 1.0d-03 * 2.0d0 * dqb
    d2cdqbt(7) = 0.0d0
    d2cdqbt(8) = 0.0d0
    d2cdqbt(9) = 1.0d-06 * 2.0d0 * dqb

  end subroutine d2coeffdqbt
  
! ******************************************************************
! ******************************************************************

  subroutine phi(nlambda,lambda,x,t,res)

    use bmproblem, only: bedel,width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: w,zb

    ! LOCAL ARRAYS
    real(kind=8) :: c(nlambda)

    call coeff(nlambda,c,x,t)

    call width(x,w)
    
    call bedel(x,zb)
  
    res = ( dot_product( c(1:nlambda), lambda(1:nlambda) ) - zb ) * w
    
  end subroutine phi

! ******************************************************************
! ******************************************************************

  subroutine dphidx(nlambda,lambda,x,t,res)
    
    use bmproblem, only: bedel,dbedel,dwidth,width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    
    ! LOCAL SCALARS
    real(kind=8) :: dw,dzb,w,zb
    
    ! LOCAL ARRAYS
    real(kind=8) :: c(nlambda),dcdx(nlambda)

    call coeff(nlambda,c,x,t)
  
    call dcoeffdx(nlambda,dcdx,x,t)

    call width(x,w)
  
    call dwidth(x,dw)
  
    call bedel(x,zb)
  
    call dbedel(x,dzb)
  
    res = ( dot_product( dcdx(1:nlambda), lambda(1:nlambda) ) - dzb ) * w + &
         ( dot_product( c(1:nlambda), lambda(1:nlambda) ) - zb ) * dw
    
  end subroutine dphidx

! ******************************************************************
! ******************************************************************

  subroutine d2phidxl(nlambda,lambda,x,t,res)
    
    use bmproblem, only: width,dwidth

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)
    
    ! LOCAL SCALARS
    real(kind=8) :: dw,w
    
    ! LOCAL ARRAYS
    real(kind=8) :: c(nlambda),dcdx(nlambda)

    if ( .false. ) write(*,*) lambda(1:nlambda) ! to avoid ``Warning: Unused dummy argument ... ''
  
    call coeff(nlambda,c,x,t)
  
    call dcoeffdx(nlambda,dcdx,x,t)

    call width(x,w)
  
    call dwidth(x,dw)
  
    res(1:nlambda) = dcdx(1:nlambda) * w + c(1:nlambda) * dw
  
  end subroutine d2phidxl

! ******************************************************************
! ******************************************************************

  subroutine dphidq(nlambda,lambda,x,t,res)
  
    use bmproblem, only: width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: w
  
    ! LOCAL ARRAYS
    real(kind=8) :: dcdqb(nlambda)
  
    call dcoeffdqb(nlambda,dcdqb,x,t)
  
    call width(x,w)
  
    res = dot_product( dcdqb(1:nlambda), lambda(1:nlambda) ) * w
  
  end subroutine dphidq

! ******************************************************************
! ******************************************************************

  subroutine d2phidql(nlambda,lambda,x,t,res)
  
    use bmproblem, only: width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: w

    ! LOCAL ARRAYS
    real(kind=8) :: dcdqb(nlambda)
  
    if ( .false. ) write(*,*) lambda(1:nlambda) ! to avoid ``Warning: Unused dummy argument ... ''
  
    call dcoeffdqb(nlambda,dcdqb,x,t)
  
    call width(x,w)
  
    res(1:nlambda) = dcdqb(1:nlambda) * w
  
  end subroutine d2phidql

! ******************************************************************
! ******************************************************************

  subroutine dphidt(nlambda,lambda,x,t,res)
  
    use bmproblem, only: dqbound,width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: dqb,w
  
    ! LOCAL ARRAYS
    real(kind=8) :: dcdqb(nlambda)
  
    call dcoeffdqb(nlambda,dcdqb,x,t)
  
    call width(x,w)
  
    call dqbound(t,dqb)
  
    res = dot_product( dcdqb(1:nlambda), lambda(1:nlambda) ) * w * dqb
  
  end subroutine dphidt

! ******************************************************************
! ******************************************************************

  subroutine d2phidtl(nlambda,lambda,x,t,res)
  
    use bmproblem, only: width,dqbound

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)

    ! LOCAL SCALARS
    real(kind=8) :: dqb,w
  
    ! LOCAL ARRAYS
    real(kind=8) :: dcdqb(nlambda)
  
    if ( .false. ) write(*,*) lambda(1:nlambda) ! to avoid ``Warning: Unused dummy argument ... ''
  
    call dcoeffdqb(nlambda,dcdqb,x,t)
    
    call width(x,w)
  
    call dqbound(t,dqb)
  
    res(1:nlambda) = dcdqb(1:nlambda) * w * dqb
  
  end subroutine d2phidtl

! ******************************************************************
! ******************************************************************

  subroutine dphidl(nlambda,lambda,x,t,res)
    
    use bmproblem, only: width

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)
    
    ! LOCAL SCALARS
    real(kind=8) :: w
    
    ! LOCAL ARRAYS
    real(kind=8) :: c(nlambda)
  
    if ( .false. ) write(*,*) lambda(1:nlambda) ! to avoid ``Warning: Unused dummy argument ... ''
  
    call coeff(nlambda,c,x,t)
  
    call width(x,w)
  
    res(1:nlambda) = c(1:nlambda) * w
  
  end subroutine dphidl

! ******************************************************************
! ******************************************************************

  subroutine uvinter(uini,uend,vini,vend,wini,wend,emptyint)

    implicit none

    ! SCALAR ARGUMENTS
    logical, intent(out) :: emptyint
    real(kind=8), intent(in) :: uini,uend,vini,vend
    real(kind=8), intent(out) :: wini,wend

    ! If the intersection bewteen [uini,uend] and [vini,vend] has a
    ! NON-empty interior, then this subroutine computes [wini,wend] as
    ! the intersection between [uini,uend] and [vini,vend] and set
    ! emptyint = .false. Otherwise, it computes nothing and set
    ! emptyint = .true.

    if ( uini .ge. uend .or. vini .ge. vend .or. vend .le. uini .or. vini .ge. uend ) then
       ! If (i) u is empty or a single point, (ii) v is empty or a
       ! single point, (iii) v is all to the left or v, or (iv) v is
       ! all to the right of u, the the intersection is empty or a
       ! single point.
       emptyint = .true.
    else
       emptyint = .false.
       if ( uini .ge. vini .and. uend .le. vend ) then
          ! u is within v and the intersection is equal to u
          wini = uini
          wend = uend
       else if ( vini .ge. uini .and. vend .le. uend ) then
          ! v is within u and the intersection is equal to v
          wini = vini
          wend = vend
       else
          wini = max( uini, vini )
          wend = min( uend, vend )
       end if
    end if
    
  end subroutine uvinter

! ******************************************************************
! ******************************************************************

  subroutine psi(nlambda,lambda,x,t,res)

    use bmproblem, only: qboundx,qbound,nwidth,widthx

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! The value of Psi(x,t) must be equal to Q(qboundx,t) minus the
    ! integral of (d/dt) Phi(x,t) dx from qboundx to x.

    ! LOCAL SCALARS
    logical :: emptyint
    integer :: i,iind
    real(kind=8) :: a,b,qb,resparca,resparcb,vleft,vright

    ! (d/dt) Phi(x,t) = P(x,t) w(x), where P(x,t) is a polynomial in x
    ! and w(x) is a piecewise polynomial in x. The integral of (d/dt)
    ! Phi(x,t) dx from qboundx to x must be computed using the
    ! piecewise definition of w(x).

    ! The definition of w(x) considers nwidth nodes named widthx(1)
    ! ... widthx(nwidth) such that w(widthx(i)) = widthval(i) for i = 1,
    ! ..., nwidth.

    ! For a given real number x:
    ! if x < widthx(1), let i=1
    ! if x >= widthx(nwidth), let i = nwidth,
    ! otherwise, let i (between 1 and nwidth-1) be such that widthx(i) <= x < widthx(i+1)
    !
    ! Then, w(x) = widthval(i) +
    !              widthb(i) * ( x - widthx(i) ) +
    !              widthc(i) * ( x - widthx(i) ) ** 2 +
    !              widthd(i) * ( x - widthx(i) ) ** 3
    
    ! Let v(1:nwidth) = widthx(1:nwidth), with v(0) = -infty and v(nwidth+1) = +infty.
    
    ! For j from 0 to nwidth, we check whether the intersection between
    ! the intervals [v(j),v(j+1)] and [qboundx,x] is an interval [a,b]
    ! with a<b. If the answer is YES, then we compute the integral
    ! between a and b and add it to the result. (A rule to stop the
    ! process before reaching j=nwidth in some situation might be used,
    ! but we prefer to ignore this possibility.)

    res = 0.0d0

    do i = 0,nwidth
       if ( i .eq. 0 ) then
          iind = 1
          vleft = - huge( 1.0d0 )
          vright = widthx(1)
       else if ( i .eq. nwidth ) then
          iind = nwidth
          vleft = widthx(nwidth)
          vright = huge( 1.0d0 )
       else
          iind = i
          vleft = widthx(i)
          vright = widthx(i+1)
       end if
       
       call uvinter(qboundx,x,vleft,vright,a,b,emptyint)

       if ( .not. emptyint ) then
          ! Compute the integral between a and b and add it to intx. The
          ! "integral" is the integral from a to b of (d/dt) Phi(x,t) dx.
          ! Within the interval [a,b], w(x) is given by the expresion
          !
          ! w(x) = widthval(iind) +
          !        widthb(iind) * ( x - widthx(iind) ) +
          !        widthc(iind) * ( x - widthx(iind) ) ** 2 +
          !        widthd(iind) * ( x - widthx(iind) ) ** 3
          !
          ! and (d/dt) Phi(x,t) is the product of the polynomial
          ! P(x,t) of degree two on x times w(x), that is a polynomial
          ! of degree three on x.
          !
          ! P(x,t) = ( r(1) + r(2) + r(3) ) +
          !          ( r(4) + r(5) + r(6) ) * x
          !          ( r(6) + r(7) + r(8) ) * x ** 2,
          !
          ! where coeffcients r are given by subroutine dcoeffdqb
          ! times d/dt Q(qboundx,t).

          call psiinter(nlambda,lambda,a,t,iind,resparca)
          call psiinter(nlambda,lambda,b,t,iind,resparcb)
          
          res = res + ( resparcb - resparca )
       end if
    end do
    
    call qbound(t,qb)

    res = qb - res

  contains

    subroutine psiinter(nlambda,lambda,x,t,iind,res)

      use bmproblem, only: widthval,widthb,widthc,widthd,dqbound
      
      implicit none
      
      ! SCALAR ARGUMENTS
      integer, intent(in) :: iind,nlambda
      real(kind=8), intent(in) :: x,t
      real(kind=8), intent(out) :: res
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: lambda(nlambda)

      ! LOCAL SCALARS
      integer :: j,k
      real(kind=8) :: const,dqb,tmp

      ! LOCAL ARRAYS
      real(kind=8) :: dcdqb(nlambda),pol1(0:2),pol2(0:3),r(nlambda)

      call dcoeffdqb(nlambda,dcdqb,x,t,.true.)
      call dqbound(t,dqb)

      r(1:nlambda) = dcdqb(1:nlambda) * dqb * lambda(1:nlambda)
 
      pol1(0) = r(1) + r(2) + r(3)
      pol1(1) = r(4) + r(5) + r(6)
      pol1(2) = r(7) + r(8) + r(9)
      
      pol2(0) = widthval(iind)
      pol2(1) = widthb(iind)
      pol2(2) = widthc(iind)
      pol2(3) = widthd(iind)

      res = 0.0d0
      do j = 0,2
         do k = 0,3
            const = pol1(j) * pol2(k)

            ! to be integrated: const * x^j * (x-xbar)^k
            ! where const = pol1(j) * pol2(k) and xbar = widthx(iind)

            if ( j .eq. 0 ) then
               if ( k .eq. 0 ) then
                  ! integrate 1
                  tmp = x
               else if ( k .eq. 1 ) then
                  ! integrate (x-xbar)
                  tmp = 0.5d0 * x ** 2 - widthx(iind) * x
               else if ( k .eq. 2 ) then
                  ! integrate (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x - widthx(iind) * x ** 2 + x ** 3 / 3.0d0
               else if ( k .eq. 3 ) then
                  ! integrate (x-xbar)^3
                  tmp = 0.25d0 * ( x - widthx(iind) ) ** 4
               end if
            else if ( j .eq. 1 ) then
               if ( k .eq. 0 ) then
                  ! integrate x
                  tmp = 0.5d0 * x ** 2
               else if ( k .eq. 1 ) then
                  ! integrate x (x-xbar)
                  tmp = x ** 3 / 3.0d0 - widthx(iind) * 0.5d0 * x ** 2
               else if ( k .eq. 2 ) then
                  ! integrate x (x-xbar)^2
                  tmp = 0.5d0 * widthx(iind) ** 2 * x ** 2 - 2.0d0 * widthx(iind) * x ** 3 / 3.0d0 + 0.25d0 * x ** 4
               else if ( k .eq. 3 ) then
                  ! integrate x (x-xbar)^3
                  tmp = - 0.5d0 * widthx(iind) ** 3 * x ** 2 + widthx(iind) ** 2 * x ** 3 - &
                       0.75d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               end if
            else if ( j .eq. 2 ) then
               if ( k .eq. 0 ) then
                  ! integrate x^2
                  tmp = x ** 3 / 3.0d0
               else if ( k .eq. 1 ) then
                  ! integrate x^2 (x-xbar)
                  tmp = 0.25d0 * x ** 4 - widthx(iind) * x ** 3 / 3.0d0
               else if ( k .eq. 2 ) then
                  ! integrate x^2 (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x ** 3 / 3.0d0 - 0.5d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               else if ( k .eq. 3 ) then
                  ! integrate x^2 (x-xbar)^3
                  tmp = - widthx(iind) ** 3 * x ** 3 / 3.0d0 + 0.75d0 * widthx(iind) ** 2 * x ** 4 - &
                       0.6d0 * widthx(iind) * x ** 5 + x ** 6 / 6.0d0
               end if
            end if
            res = res + const * tmp
         end do
      end do
      
    end subroutine psiinter
        
  end subroutine psi

! ******************************************************************
! ******************************************************************

  subroutine dpsidl(nlambda,lambda,x,t,res)

    use bmproblem, only: nwidth,widthx,qboundx
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)

    ! LOCAL SCALARS
    logical :: emptyint
    integer :: i,iind
    real(kind=8) :: a,b,vleft,vright

    ! LOCAL ARRAYS
    real(kind=8) :: resa(nlambda),resb(nlambda)
    
    res(1:nlambda) = 0.0d0

    do i = 0,nwidth
       if ( i .eq. 0 ) then
          iind = 1
          vleft = - huge( 1.0d0 )
          vright = widthx(1)
       else if ( i .eq. nwidth ) then
          iind = nwidth
          vleft = widthx(nwidth)
          vright = huge( 1.0d0 )
       else
          iind = i
          vleft = widthx(i)
          vright = widthx(i+1)
       end if
       
       call uvinter(qboundx,x,vleft,vright,a,b,emptyint)

       if ( .not. emptyint ) then
          call dpsiinterdl(nlambda,lambda,a,t,iind,resa)
          call dpsiinterdl(nlambda,lambda,b,t,iind,resb)
          res(1:nlambda) = res(1:nlambda) - ( resb(1:nlambda) - resa(1:nlambda) )
       end if
    end do

  contains

    subroutine dpsiinterdl(nlambda,lambda,x,t,iind,res)

      use bmproblem, only: widthval,widthb,widthc,widthd,dqbound
      
      implicit none
      
      ! SCALAR ARGUMENTS
      integer, intent(in) :: iind,nlambda
      real(kind=8), intent(in) :: x,t
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: lambda(nlambda)
      real(kind=8), intent(out) :: res(nlambda)

      ! LOCAL SCALARS
      integer :: j,k
      real(kind=8) :: dqb,tmp

      ! LOCAL ARRAYS
      real(kind=8) :: dcdqb(nlambda),pol2(0:3),dpol1dl(nlambda,0:2),dconstdl(nlambda)

      if ( .false. ) write(*,*) lambda(1:nlambda)
      
      call dcoeffdqb(nlambda,dcdqb,x,t,.true.)
      call dqbound(t,dqb)

      ! r(1:nlambda) = dcdqb(1:nlambda) * dqb * lambda(1:nlambda)
      ! pol1(0) = r(1) + r(2) + r(3)
      ! pol1(1) = r(4) + r(5) + r(6)
      ! pol1(2) = r(7) + r(8) + r(9)

      dpol1dl(1:nlambda,0) = (/ dcdqb(1), dcdqb(2), dcdqb(3), 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0 /) * dqb
      dpol1dl(1:nlambda,1) = (/ 0.0d0, 0.0d0, 0.0d0, dcdqb(4), dcdqb(5), dcdqb(6), 0.0d0, 0.0d0, 0.0d0 /) * dqb
      dpol1dl(1:nlambda,2) = (/ 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, dcdqb(7), dcdqb(8), dcdqb(9) /) * dqb
      
      pol2(0) = widthval(iind)
      pol2(1) = widthb(iind)
      pol2(2) = widthc(iind)
      pol2(3) = widthd(iind)

      res(1:nlambda) = 0.0d0
      
      do j = 0,2
         do k = 0,3
            ! const = pol1(j) * pol2(k)
            dconstdl(1:nlambda) = dpol1dl(1:nlambda,j) * pol2(k)
            
            ! to be integrated: const * x^j * (x-xbar)^k
            ! where const = pol1(j) * pol2(k) and xbar = widthx(iind)

            if ( j .eq. 0 ) then
               if ( k .eq. 0 ) then
                  ! integrate 1
                  tmp = x
               else if ( k .eq. 1 ) then
                  ! integrate (x-xbar)
                  tmp = 0.5d0 * x ** 2 - widthx(iind) * x
               else if ( k .eq. 2 ) then
                  ! integrate (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x - widthx(iind) * x ** 2 + x ** 3 / 3.0d0
               else if ( k .eq. 3 ) then
                  ! integrate (x-xbar)^3
                  tmp = 0.25d0 * ( x - widthx(iind) ) ** 4
               end if
            else if ( j .eq. 1 ) then
               if ( k .eq. 0 ) then
                  ! integrate x
                  tmp = 0.5d0 * x ** 2
               else if ( k .eq. 1 ) then
                  ! integrate x (x-xbar)
                  tmp = x ** 3 / 3.0d0 - widthx(iind) * 0.5d0 * x ** 2
               else if ( k .eq. 2 ) then
                  ! integrate x (x-xbar)^2
                  tmp = 0.5d0 * widthx(iind) ** 2 * x ** 2 - 2.0d0 * widthx(iind) * x ** 3 / 3.0d0 + 0.25d0 * x ** 4
               else if ( k .eq. 3 ) then
                  ! integrate x (x-xbar)^3
                  tmp = - 0.5d0 * widthx(iind) ** 3 * x ** 2 + widthx(iind) ** 2 * x ** 3 - &
                       0.75d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               end if
            else if ( j .eq. 2 ) then
               if ( k .eq. 0 ) then
                  ! integrate x^2
                  tmp = x ** 3 / 3.0d0
               else if ( k .eq. 1 ) then
                  ! integrate x^2 (x-xbar)
                  tmp = 0.25d0 * x ** 4 - widthx(iind) * x ** 3 / 3.0d0
               else if ( k .eq. 2 ) then
                  ! integrate x^2 (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x ** 3 / 3.0d0 - 0.5d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               else if ( k .eq. 3 ) then
                  ! integrate x^2 (x-xbar)^3
                  tmp = - widthx(iind) ** 3 * x ** 3 / 3.0d0 + 0.75d0 * widthx(iind) ** 2 * x ** 4 - &
                       0.6d0 * widthx(iind) * x ** 5 + x ** 6 / 6.0d0
               end if
            end if
            ! res = res + const * tmp
            res(1:nlambda) = res(1:nlambda) + dconstdl(1:nlambda) * tmp
         end do
      end do
      
    end subroutine dpsiinterdl
            
end subroutine dpsidl

! ******************************************************************
! ******************************************************************

  subroutine dpsidt(nlambda,lambda,x,t,res)

    use bmproblem, only: nwidth,widthx,qboundx,dqbound
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    real(kind=8), intent(out) :: res
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)

    ! LOCAL SCALARS
    logical :: emptyint
    integer :: i,iind
    real(kind=8) :: a,b,dqb,resparca,resparcb,vleft,vright

    res = 0.0d0

    do i = 0,nwidth
       if ( i .eq. 0 ) then
          iind = 1
          vleft = - huge( 1.0d0 )
          vright = widthx(1)
       else if ( i .eq. nwidth ) then
          iind = nwidth
          vleft = widthx(nwidth)
          vright = huge( 1.0d0 )
       else
          iind = i
          vleft = widthx(i)
          vright = widthx(i+1)
       end if
       
       call uvinter(qboundx,x,vleft,vright,a,b,emptyint)

       if ( .not. emptyint ) then
          call dpsiinterdt(nlambda,lambda,a,t,iind,resparca)
          call dpsiinterdt(nlambda,lambda,b,t,iind,resparcb)
          res = res + ( resparcb - resparca )
       end if
    end do

    call dqbound(t,dqb)

    res = dqb - res

  contains

    subroutine dpsiinterdt(nlambda,lambda,x,t,iind,res)

      use bmproblem, only: widthval,widthb,widthc,widthd,dqbound,ddqbound
      
      implicit none
      
      ! SCALAR ARGUMENTS
      integer, intent(in) :: iind,nlambda
      real(kind=8), intent(in) :: x,t
      real(kind=8), intent(out) :: res
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: lambda(nlambda)

      ! LOCAL SCALARS
      integer :: j,k
      real(kind=8) :: dconstdt,ddqb,dqb,tmp

      ! LOCAL ARRAYS
      real(kind=8) :: dcdqb(nlambda),d2cdqbt(nlambda),pol2(0:3),drdt(nlambda),dpol1dt(0:2)

      call dcoeffdqb(nlambda,dcdqb,x,t,.true.)
      call d2coeffdqbt(nlambda,d2cdqbt,t)
      
      call dqbound(t,dqb)
      call ddqbound(t,ddqb)
      
      ! r(1:nlambda) = dcdqb(1:nlambda) * dqb * lambda(1:nlambda)
      drdt(1:nlambda) = ( d2cdqbt(1:nlambda) * dqb + dcdqb(1:nlambda) * ddqb ) * lambda(1:nlambda)

      ! pol1(0) = r(1) + r(2) + r(3)
      ! pol1(1) = r(4) + r(5) + r(6)
      ! pol1(2) = r(7) + r(8) + r(9)
      dpol1dt(0) = drdt(1) + drdt(2) + drdt(3)
      dpol1dt(1) = drdt(4) + drdt(5) + drdt(6)
      dpol1dt(2) = drdt(7) + drdt(8) + drdt(9)
      
      pol2(0) = widthval(iind)
      pol2(1) = widthb(iind)
      pol2(2) = widthc(iind)
      pol2(3) = widthd(iind)

      res = 0.0d0
      
      do j = 0,2
         do k = 0,3
            ! const = pol1(j) * pol2(k)
            dconstdt = dpol1dt(j) * pol2(k)
            
            ! to be integrated: const * x^j * (x-xbar)^k
            ! where const = pol1(j) * pol2(k) and xbar = widthx(iind)

            if ( j .eq. 0 ) then
               if ( k .eq. 0 ) then
                  ! integrate 1
                  tmp = x
               else if ( k .eq. 1 ) then
                  ! integrate (x-xbar)
                  tmp = 0.5d0 * x ** 2 - widthx(iind) * x
               else if ( k .eq. 2 ) then
                  ! integrate (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x - widthx(iind) * x ** 2 + x ** 3 / 3.0d0
               else if ( k .eq. 3 ) then
                  ! integrate (x-xbar)^3
                  tmp = 0.25d0 * ( x - widthx(iind) ) ** 4
               end if
            else if ( j .eq. 1 ) then
               if ( k .eq. 0 ) then
                  ! integrate x
                  tmp = 0.5d0 * x ** 2
               else if ( k .eq. 1 ) then
                  ! integrate x (x-xbar)
                  tmp = x ** 3 / 3.0d0 - widthx(iind) * 0.5d0 * x ** 2
               else if ( k .eq. 2 ) then
                  ! integrate x (x-xbar)^2
                  tmp = 0.5d0 * widthx(iind) ** 2 * x ** 2 - 2.0d0 * widthx(iind) * x ** 3 / 3.0d0 + 0.25d0 * x ** 4
               else if ( k .eq. 3 ) then
                  ! integrate x (x-xbar)^3
                  tmp = - 0.5d0 * widthx(iind) ** 3 * x ** 2 + widthx(iind) ** 2 * x ** 3 - &
                       0.75d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               end if
            else if ( j .eq. 2 ) then
               if ( k .eq. 0 ) then
                  ! integrate x^2
                  tmp = x ** 3 / 3.0d0
               else if ( k .eq. 1 ) then
                  ! integrate x^2 (x-xbar)
                  tmp = 0.25d0 * x ** 4 - widthx(iind) * x ** 3 / 3.0d0
               else if ( k .eq. 2 ) then
                  ! integrate x^2 (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x ** 3 / 3.0d0 - 0.5d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               else if ( k .eq. 3 ) then
                  ! integrate x^2 (x-xbar)^3
                  tmp = - widthx(iind) ** 3 * x ** 3 / 3.0d0 + 0.75d0 * widthx(iind) ** 2 * x ** 4 - &
                       0.6d0 * widthx(iind) * x ** 5 + x ** 6 / 6.0d0
               end if
            end if
            ! res = res + const * tmp
            res = res + dconstdt * tmp
         end do
      end do
      
    end subroutine dpsiinterdt
            
  end subroutine dpsidt

! ******************************************************************
! ******************************************************************

  subroutine d2psidtl(nlambda,lambda,x,t,res)

    use bmproblem, only: nwidth,widthx,qboundx
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: nlambda
    real(kind=8), intent(in) :: x,t
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: lambda(nlambda)
    real(kind=8), intent(out) :: res(nlambda)

    ! LOCAL SCALARS
    logical :: emptyint
    integer :: i,iind
    real(kind=8) :: a,b,vleft,vright

    ! LOCAL ARRAYS
    real(kind=8) :: resparca(nlambda),resparcb(nlambda)
    
    res(1:nlambda) = 0.0d0

    do i = 0,nwidth
       if ( i .eq. 0 ) then
          iind = 1
          vleft = - huge( 1.0d0 )
          vright = widthx(1)
       else if ( i .eq. nwidth ) then
          iind = nwidth
          vleft = widthx(nwidth)
          vright = huge( 1.0d0 )
       else
          iind = i
          vleft = widthx(i)
          vright = widthx(i+1)
       end if
       
       call uvinter(qboundx,x,vleft,vright,a,b,emptyint)

       if ( .not. emptyint ) then
          call d2psiinterdtl(nlambda,lambda,a,t,iind,resparca)
          call d2psiinterdtl(nlambda,lambda,b,t,iind,resparcb)
          res(1:nlambda) = res(1:nlambda) + ( resparcb(1:nlambda) - resparca(1:nlambda) )
       end if
    end do

    res(1:nlambda) = - res(1:nlambda)
    
  contains

    subroutine d2psiinterdtl(nlambda,lambda,x,t,iind,res)

      use bmproblem, only: widthval,widthb,widthc,widthd,dqbound,ddqbound
      
      implicit none
      
      ! SCALAR ARGUMENTS
      integer, intent(in) :: iind,nlambda
      real(kind=8), intent(in) :: x,t
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: lambda(nlambda)
      real(kind=8), intent(out) :: res(nlambda)

      ! LOCAL SCALARS
      integer :: j,k
      real(kind=8) :: ddqb,dqb,tmp

      ! LOCAL ARRAYS
      real(kind=8) :: dcdqb(nlambda),d2cdqbt(nlambda),d2constdtl(nlambda),pol2(0:3),d2rdtl(nlambda),d2pol1dtl(nlambda,0:2)

      if ( .false. ) write(*,*) lambda(1:nlambda)
      
      call dcoeffdqb(nlambda,dcdqb,x,t,.true.)
      call d2coeffdqbt(nlambda,d2cdqbt,t)
      
      call dqbound(t,dqb)
      call ddqbound(t,ddqb)
      
      ! r(1:nlambda) = dcdqb(1:nlambda) * dqb * lambda(1:nlambda)
      ! drdt(1:nlambda) = ( d2cdqbt(1:nlambda) * dqb + dcdqb(1:nlambda) * ddqb ) * lambda(1:nlambda)
      d2rdtl(1:nlambda) = ( d2cdqbt(1:nlambda) * dqb + dcdqb(1:nlambda) * ddqb )

      ! pol1(0) = r(1) + r(2) + r(3)
      ! pol1(1) = r(4) + r(5) + r(6)
      ! pol1(2) = r(7) + r(8) + r(9)
      ! dpol1dt(0) = drdt(1) + drdt(2) + drdt(3)
      ! dpol1dt(1) = drdt(4) + drdt(5) + drdt(6)
      ! dpol1dt(2) = drdt(7) + drdt(8) + drdt(9)
      d2pol1dtl(1:nlambda,0) = (/ d2rdtl(1), d2rdtl(2), d2rdtl(3), 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0 /)
      d2pol1dtl(1:nlambda,1) = (/ 0.0d0, 0.0d0, 0.0d0, d2rdtl(4), d2rdtl(5), d2rdtl(6), 0.0d0, 0.0d0, 0.0d0 /)
      d2pol1dtl(1:nlambda,2) = (/ 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, d2rdtl(7), d2rdtl(8), d2rdtl(9) /)
      
      pol2(0) = widthval(iind)
      pol2(1) = widthb(iind)
      pol2(2) = widthc(iind)
      pol2(3) = widthd(iind)

      res(1:nlambda) = 0.0d0
      
      do j = 0,2
         do k = 0,3
            ! const = pol1(j) * pol2(k)
            ! dconstdt = dpol1dt(j) * pol2(k)
            d2constdtl(1:nlambda) = d2pol1dtl(1:nlambda,j) * pol2(k)
            
            ! to be integrated: const * x^j * (x-xbar)^k
            ! where const = pol1(j) * pol2(k) and xbar = widthx(iind)

            if ( j .eq. 0 ) then
               if ( k .eq. 0 ) then
                  ! integrate 1
                  tmp = x
               else if ( k .eq. 1 ) then
                  ! integrate (x-xbar)
                  tmp = 0.5d0 * x ** 2 - widthx(iind) * x
               else if ( k .eq. 2 ) then
                  ! integrate (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x - widthx(iind) * x ** 2 + x ** 3 / 3.0d0
               else if ( k .eq. 3 ) then
                  ! integrate (x-xbar)^3
                  tmp = 0.25d0 * ( x - widthx(iind) ) ** 4
               end if
            else if ( j .eq. 1 ) then
               if ( k .eq. 0 ) then
                  ! integrate x
                  tmp = 0.5d0 * x ** 2
               else if ( k .eq. 1 ) then
                  ! integrate x (x-xbar)
                  tmp = x ** 3 / 3.0d0 - widthx(iind) * 0.5d0 * x ** 2
               else if ( k .eq. 2 ) then
                  ! integrate x (x-xbar)^2
                  tmp = 0.5d0 * widthx(iind) ** 2 * x ** 2 - 2.0d0 * widthx(iind) * x ** 3 / 3.0d0 + 0.25d0 * x ** 4
               else if ( k .eq. 3 ) then
                  ! integrate x (x-xbar)^3
                  tmp = - 0.5d0 * widthx(iind) ** 3 * x ** 2 + widthx(iind) ** 2 * x ** 3 - &
                       0.75d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               end if
            else if ( j .eq. 2 ) then
               if ( k .eq. 0 ) then
                  ! integrate x^2
                  tmp = x ** 3 / 3.0d0
               else if ( k .eq. 1 ) then
                  ! integrate x^2 (x-xbar)
                  tmp = 0.25d0 * x ** 4 - widthx(iind) * x ** 3 / 3.0d0
               else if ( k .eq. 2 ) then
                  ! integrate x^2 (x-xbar)^2
                  tmp = widthx(iind) ** 2 * x ** 3 / 3.0d0 - 0.5d0 * widthx(iind) * x ** 4 + 0.2d0 * x ** 5
               else if ( k .eq. 3 ) then
                  ! integrate x^2 (x-xbar)^3
                  tmp = - widthx(iind) ** 3 * x ** 3 / 3.0d0 + 0.75d0 * widthx(iind) ** 2 * x ** 4 - &
                       0.6d0 * widthx(iind) * x ** 5 + x ** 6 / 6.0d0
               end if
            end if
            ! res = res + const * tmp
            ! res = res + dconstdt * tmp
            res(1:nlambda) = res(1:nlambda) + d2constdtl(1:nlambda) * tmp
         end do
      end do
      
    end subroutine d2psiinterdtl
            
  end subroutine d2psidtl

! ******************************************************************
! ******************************************************************

subroutine dpsidx(nlambda,lambda,x,t,res)
    
  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(in) :: x,t
  real(kind=8), intent(out) :: res
    
  ! ARRAY ARGUMENTS
  real(kind=8), intent(in) :: lambda(nlambda)
    
  call dphidt(nlambda,lambda,x,t,res)
  res = - res
  
end subroutine dpsidx

! ******************************************************************
! ******************************************************************

subroutine d2psidxl(nlambda,lambda,x,t,res)
    
  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(in) :: x,t
    
  ! ARRAY ARGUMENTS
  real(kind=8), intent(in) :: lambda(nlambda)
  real(kind=8), intent(out) :: res(nlambda)
    
  call d2phidtl(nlambda,lambda,x,t,res)
  res(1:nlambda) = - res(1:nlambda)

end subroutine d2psidxl

! ******************************************************************
! ******************************************************************

subroutine checkdlxt(nlambda,lambda,x,t,fun,dfundl,dfundx,dfundt)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(inout) :: x,t

  ! ARRAY ARGUMENTS
  real(kind=8), intent(inout) :: lambda(nlambda)

  ! SUBROUTINE ARGUMENTS
  external :: fun,dfundl,dfundx,dfundt
  
  ! LOCAL SCALARS
  integer :: i,ih
  real(kind=8) :: bestdiff,besth,diff,fplus,fminus,h,tmp
  
  ! LOCAL ARRAYS
  real(kind=8) :: g(nlambda+2),gapprox(nlambda+2)

  call dfundl(nlambda,lambda,x,t,g)
  
  call dfundx(nlambda,lambda,x,t,g(nlambda+1))
  
  call dfundt(nlambda,lambda,x,t,g(nlambda+2))
  
  bestdiff = huge( 1.0d0 )

  do ih = 1,17
     if ( ih .le. 16 ) then
        h = 10.0d0 ** ( - ih )
     else
        h = besth
     end if

     do i = 1,nlambda
        tmp = lambda(i)
     
        lambda(i) = tmp + h * max( 1.0d0, abs( tmp ) )
        call fun(nlambda,lambda,x,t,fplus)
     
        lambda(i) = tmp - h * max( 1.0d0, abs( tmp ) )
        call fun(nlambda,lambda,x,t,fminus)

        lambda(i) = tmp

        gapprox(i) = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )
     end do

     tmp = x
     
     x = tmp + h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fplus)
     
     x = tmp - h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fminus)

     x = tmp
     
     gapprox(nlambda+1) = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )

     tmp = t
     
     t = tmp + h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fplus)
     
     t = tmp - h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fminus)

     t = tmp
     
     gapprox(nlambda+2) = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )
     
     diff = maxval( abs( g(1:nlambda+2) - gapprox(1:nlambda+2) ) )

     if ( ih .le. 16 ) then
        if ( diff .lt. bestdiff ) then
           bestdiff = diff
           besth = h
        end if

     else
        write(*,*) 'bestdiff = ',bestdiff
        write(*,*) 'besth = ',besth
        do i = 1,nlambda + 2
           write(*,*) i,g(i),gapprox(i),abs( g(i) - gapprox(i) )
        end do
     end if
  end do
  
end subroutine checkdlxt

! ******************************************************************
! ******************************************************************

subroutine checkdx(nlambda,lambda,x,t,fun,dfundx)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(inout) :: x,t

  ! ARRAY ARGUMENTS
  real(kind=8), intent(inout) :: lambda(nlambda)

  ! SUBROUTINE ARGUMENTS
  external :: fun,dfundl,dfundx,dfundt
  
  ! LOCAL SCALARS
  integer :: ih
  real(kind=8) :: bestdiff,besth,diff,fplus,fminus,g,gapprox,h,tmp
  
  call dfundx(nlambda,lambda,x,t,g)
  
  bestdiff = huge( 1.0d0 )

  do ih = 1,17
     if ( ih .le. 16 ) then
        h = 10.0d0 ** ( - ih )
     else
        h = besth
     end if

     tmp = x
     
     x = tmp + h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fplus)
     
     x = tmp - h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fminus)

     x = tmp
     
     gapprox = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )

     diff = abs( g - gapprox )

     if ( ih .le. 16 ) then
        if ( diff .lt. bestdiff ) then
           bestdiff = diff
           besth = h
        end if

     else
        write(*,*) 'bestdiff = ',bestdiff
        write(*,*) 'besth = ',besth
        write(*,*) g,gapprox,abs( g - gapprox )
     end if
  end do
  
end subroutine checkdx

! ******************************************************************
! ******************************************************************

subroutine checkdl(nlambda,lambda,x,t,fun,dfundl)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(inout) :: x,t

  ! ARRAY ARGUMENTS
  real(kind=8), intent(inout) :: lambda(nlambda)

  ! SUBROUTINE ARGUMENTS
  external :: fun,dfundl
  
  ! LOCAL SCALARS
  integer :: i,ih
  real(kind=8) :: bestdiff,besth,diff,fplus,fminus,h,tmp
  
  ! LOCAL ARRAYS
  real(kind=8) :: g(nlambda),gapprox(nlambda)

  call dfundl(nlambda,lambda,x,t,g)
  
  bestdiff = huge( 1.0d0 )

  do ih = 1,17
     if ( ih .le. 16 ) then
        h = 10.0d0 ** ( - ih )
     else
        h = besth
     end if

     do i = 1,nlambda
        tmp = lambda(i)
     
        lambda(i) = tmp + h * max( 1.0d0, abs( tmp ) )
        call fun(nlambda,lambda,x,t,fplus)
     
        lambda(i) = tmp - h * max( 1.0d0, abs( tmp ) )
        call fun(nlambda,lambda,x,t,fminus)

        lambda(i) = tmp

        gapprox(i) = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )
     end do

     diff = maxval( abs( g(1:nlambda) - gapprox(1:nlambda) ) )

     if ( ih .le. 16 ) then
        if ( diff .lt. bestdiff ) then
           bestdiff = diff
           besth = h
        end if

     else
        write(*,*) 'bestdiff = ',bestdiff
        write(*,*) 'besth = ',besth
        do i = 1,nlambda
           write(*,*) i,g(i),gapprox(i),abs( g(i) - gapprox(i) )
        end do
     end if
  end do
  
end subroutine checkdl

! ******************************************************************
! ******************************************************************

subroutine checkdt(nlambda,lambda,x,t,fun,dfundt)

  implicit none

  ! SCALAR ARGUMENTS
  integer, intent(in) :: nlambda
  real(kind=8), intent(inout) :: x,t

  ! ARRAY ARGUMENTS
  real(kind=8), intent(inout) :: lambda(nlambda)

  ! SUBROUTINE ARGUMENTS
  external :: fun,dfundl,dfundx,dfundt
  
  ! LOCAL SCALARS
  integer :: ih
  real(kind=8) :: bestdiff,besth,diff,fplus,fminus,g,gapprox,h,tmp
  
  call dfundt(nlambda,lambda,x,t,g)
  
  bestdiff = huge( 1.0d0 )

  do ih = 1,17
     if ( ih .le. 16 ) then
        h = 10.0d0 ** ( - ih )
     else
        h = besth
     end if

     tmp = t
     
     t = tmp + h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fplus)
     
     t = tmp - h * max( 1.0d0, abs( tmp ) )
     call fun(nlambda,lambda,x,t,fminus)

     t = tmp
     
     gapprox = ( fplus - fminus ) / ( 2.0d0 * h * max( 1.0d0, abs( tmp ) ) )
     
     diff = abs( g - gapprox )

     if ( ih .le. 16 ) then
        if ( diff .lt. bestdiff ) then
           bestdiff = diff
           besth = h
        end if

     else
        write(*,*) 'bestdiff = ',bestdiff
        write(*,*) 'besth = ',besth
        write(*,*) g,gapprox,abs( g - gapprox )
     end if
  end do
  
end subroutine checkdt

end module bmmodel
