module bmproblem

  use iso_c_binding, only: c_ptr,c_f_pointer
  
  implicit none
  
  type :: pdata_type
     integer :: m,ndiv,nobs,nr
     real(kind=8) :: dt,dx,tmax,xmax
     logical, allocatable :: lostobs(:)
     integer, allocatable :: ind(:)
     real(kind=8), allocatable :: xsol(:),xtmp(:),y(:),yobs(:),lefty(:),righy(:),realizations(:,:)
  end type pdata_type

  public :: iniprob,endprob,inip,evalf,evalls

contains

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

  subroutine iniprob(n,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(out) :: n
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    integer :: allocstat,i
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    write(*,*) 'Enter xmax and tmax: '
    read(*,*) pdata%xmax,pdata%tmax

    pdata%dt =  1.0d0
    pdata%dx = 60.0d0
    
    pdata%ndiv = ceiling( pdata%xmax / pdata%dx ) 
    pdata%nobs = 2 * ceiling( pdata%tmax / pdata%dt ) * ( pdata%ndiv + 1 )

    allocate(pdata%xsol(pdata%ndiv),pdata%xtmp(pdata%ndiv),pdata%y(pdata%nobs), &
         pdata%yobs(pdata%nobs),pdata%lostobs(pdata%nobs),pdata%lefty(pdata%nobs), &
         pdata%righy(pdata%nobs),pdata%ind(pdata%nobs),stat=allocstat)
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate.'
       stop
    end if
    
    call saintve32(pdata%ndiv,pdata%xsol,pdata%nobs,pdata%yobs,pdata%lostobs,pdata%lefty,pdata%righy,pdata%dt,pdata%dx,1)

    pdata%m = count( .not. pdata%lostobs(1:pdata%nobs) )
    pdata%ind(1:pdata%m) = pack( (/ (i,i=1,pdata%nobs) /), .not. pdata%lostobs(1:pdata%nobs) )
    
    write(*,*) 'Number of non-lost: ',pdata%m

    pdata%nr = 100
    
    allocate(pdata%realizations(pdata%nr,pdata%m),stat=allocstat)
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate.'
       stop
    end if
    
!!$    n = pdata%ndiv
    
    n = 2

  end subroutine iniprob

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

  subroutine endprob(pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    integer :: allocstat
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    deallocate(pdata%xsol,pdata%xtmp,pdata%y,pdata%yobs,pdata%lostobs, &
         pdata%lefty,pdata%righy,pdata%ind,pdata%realizations,stat=allocstat)
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate statement!'
       stop
    end if

  end subroutine endprob

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

  subroutine inip(n,x,l,u,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(out) :: l(n),u(n),x(n)

    ! LOCAL SCALARS
    integer :: k1,k2,kman,ksig
    real(kind=8) :: f,seed,mansol,sigsol,mansize,sigsize
    type(pdata_type), pointer :: pdata

    ! EXTERNAL FUNCTIONS
    real(kind=8), external :: drand

    call c_f_pointer(pdataptr,pdata)

!!$    x(1:pdata%ndiv) = 0.0d0
!!$    l(1:pdata%ndiv) = - 1.0d+20
!!$    u(1:pdata%ndiv) =   1.0d+20

  ! x(1:2) = (/ 3.3663306720355467d-02, 0.0d0 /)
    x(1:2) = (/ 0.0d0, 0.1d0 /)

    l(1:2) = (/ 0.0d0,   0.0d0 /)
    u(1:2) = (/ 1.0d+20, 1.0d+20 /)

!!$    x(1) = 4.3089031638940067d-002
!!$    x(2) = 7.5000002980232239d-002
!!$    
!!$    call evalf(n,x,f,pdataptr)
!!$    write(*,*) 'At the solution x = ',x(1),x(2),' f = ',f
!!$    stop
    
!!$    mansol = 3.3663306720355467d-02
!!$    sigsol = 0.0d0
!!$
!!$    mansize = 1.0d0 * abs( mansol )
!!$    sigsize = 1.0d0 !0.4d0
!!$
!!$    kman = 51
!!$    ksig = 41
!!$    
!!$    seed = 123456.0d0
!!$    do k1 = 1,kman
!!$       do k2 = 1,ksig
!!$          x(1) = mansol + ( ( real( k1 - 1 ) / ( kman - 1 ) ) - 0.5d0 ) * mansize
!!$          x(2) = sigsol + ( ( real( k2 - 1 ) / ( ksig - 1 ) )         ) * sigsize
!!$          call evalf(n,x,f,pdataptr)
!!$        ! write(*,*) 'Tentativa (k1,k2) = (',k1,',',k2,')',' x = ',x(1),x(2),' f = ',f
!!$          write(*,*) k1,k2,x(1),x(2),f
!!$       end do
!!$    end do
!!$    
!!$    stop
    
    mansol = 3.3663306720355467d-02
    sigsol = 0.0d0

    mansize = 0.05d0 * abs( mansol )

    kman = 101
    
    do k1 = 1,kman
       x(1) = mansol + ( ( real( k1 - 1 ) / ( kman - 1 ) ) - 0.5d0 ) * mansize
       x(2) = sigsol
       call evalf(n,x,f,pdataptr)
       write(*,*) k1,x(1),x(2),f
    end do

!!$    mansol = 3.3663306720355467d-02
!!$    sigsol = 0.0d0
!!$
!!$    sigsize = 1.0d-03
!!$
!!$    ksig = 101
!!$    
!!$    do k2 = 1,ksig
!!$       x(1) = mansol
!!$       x(2) = sigsol + ( ( real( k2 - 1 ) / ( ksig - 1 ) )         ) * sigsize
!!$       call evalf(n,x,f,pdataptr)
!!$       write(*,*) k2,x(1),x(2),f
!!$    end do

    stop
    
  end subroutine inip

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

  subroutine evalf(n,x,f,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    real(kind=8), intent(out) :: f
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)

    ! LOCAL SCALARS
    integer :: i,j
    real(kind=8) :: naca,nacb,seed,sigma
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    pdata%xtmp(1:pdata%ndiv) = x(1)

    sigma = x(2)

    seed = 9388847237.d0

    do i = 1,pdata%nr
       call saintve32(pdata%ndiv,pdata%xtmp,pdata%nobs,pdata%y,pdata%lostobs,pdata%lefty,pdata%righy,pdata%dt,pdata%dx,2,sigma,seed)
       pdata%realizations(i,1:pdata%m) = pdata%y(pdata%ind(1:pdata%m)) ! This is row oriented !!!
    end do

  ! f = 0.0d0
    f = - 1.0d0
    do j = 1,pdata%m
       call computenac(pdata%nr,pdata%realizations(1:pdata%nr,j),pdata%lefty(pdata%ind(j)),naca)
       call computenac(pdata%nr,pdata%realizations(1:pdata%nr,j),pdata%righy(pdata%ind(j)),nacb)
     ! f = f - log( ( nacb - naca ) / pdata%nr )
     ! f = f - ( nacb - naca ) / pdata%nr
       f = f * ( nacb - naca ) / pdata%nr
    end do
    
!!$    write(*,*) ' x = ',x(1),' sigma = ',x(2),' f = ',f

  end subroutine evalf
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine evalls(n,x,ls,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    real(kind=8), intent(out) :: ls
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)

    ! LOCAL SCALARS
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    pdata%xtmp(1:pdata%ndiv) = x(1)
    
    call saintve32(pdata%ndiv,pdata%xtmp,pdata%nobs,pdata%y,pdata%lostobs,pdata%lefty,pdata%righy,pdata%dt,pdata%dx,2)

    ls = sum( ( pdata%y(pdata%ind(1:pdata%m)) - pdata%yobs(pdata%ind(1:pdata%m)) ) ** 2 )
        
  end subroutine evalls
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine computenac(n,v,x,nacx)

    implicit none
    
    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    real(kind=8), intent(in) :: x
    real(kind=8), intent(out) :: nacx
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: v(n)

    ! nac(v(i)) corresponds to the position (from 1 to n) of v(i) in a
    ! sorted version of v

    ! Let vmin and vmax the smalltest and the largest elements of v,
    ! respectively.

    ! If x < vmin, then nac(x)=0.
    ! If x > vmax, then nac(x)=n.

    ! Otherwise, let vleft and vright be such that vleft <= x <=
    ! vright, where vleft and vright are as large and as small as
    ! possible, respectively.

    ! If vleft=vright, then nac(x) = nac(vleft). Otherwise,
    ! nac(x) = nac(vleft) + [ ( x - vleft ) / ( vright - vleft ) ] * ( nac(vright) - nac(vleft) )
    
    ! LOCAL SCALARS
    integer :: i
    real(kind=8) :: vmin,vmax,vleft,vright,nacvleft,nacvright
    
    if ( n .lt. 2 ) then
       write(*,*) 'n must be at least 2.'
       return
    end if

    vmin = minval( v(1:n) )
    vmax = maxval( v(1:n) )

    if ( x .lt. vmin ) then
       nacx = 0.0d0
       
    else if ( x .gt. vmax ) then
       nacx = real( n )
       
    else
       vleft  = - huge( 1.0d0 )
       vright =   huge( 1.0d0 )
    
       do i = 1,n
          if ( v(i) .ge. vleft  .and. v(i) .le. x ) vleft  = v(i)
          if ( v(i) .le. vright .and. v(i) .ge. x ) vright = v(i)
       end do

       nacvleft  = count( v(1:n) .le. vleft  )

       if ( vleft .eq. vright ) then
          nacx = nacvleft
          
       else
          nacvright = count( v(1:n) .le. vright )
          nacx = nacvleft + ( ( x - vleft ) / ( vright - vleft )  ) * ( nacvright - nacvleft )
       end if
    end if

  end subroutine computenac
        
  ! ******************************************************************
  ! ******************************************************************

  subroutine saintve32(ndiv,mann,nobs,y,losty,lefty,righy,dt,dx,modo,sigma,seed)
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: modo,ndiv,nobs
    real(kind=8), intent(in) :: dt,dx
    real(kind=8), optional, intent(in) :: sigma
    real(kind=8), optional, intent(inout) :: seed
    
    ! ARRAY ARGUMENTS
    logical, intent(inout) :: losty(nobs) ! out when modo=1, unused when modo=2
    real(kind=8), intent(inout) :: mann(ndiv),lefty(nobs),righy(nobs) ! out when modo=1, in when modo=2
    real(kind=8), intent(out) :: y(nobs)

    ! Canal retangular 
    ! Objetivo: simular os resultados do Exemplo 14.9 de "Hidraulica Basica", Porto, 2006.
    ! Método explícito
    ! Condicao de contorno sobre vazão em t = 0
    ! Area transversal upwind em x = 0 e x = xmax
    ! Velocidade upwind em x = xmax
    ! V(x,0) se tira da condição de contorno sobre a vazão AV
    ! 31 de julho de 2020.
    
    ! LOCAL SCALARS
    integer :: i,npun,k
    real(kind=8) :: width,ainic,vinic,avinic,hinic,efe,grav,pinic,rhinic, &
         zbx,zxinic,t,avx,at,alfa,ualfa,av2x,zhatx,adifu,pdifu,avt,vdifu, &
         avdifu,hdifu,localseed
    
    ! LOCAL ARRAYS
    real(kind=8) :: a(ndiv+1),anew(ndiv+1),av(ndiv+1),avnew(ndiv+1),h(ndiv+1), &
         hnew(ndiv+1),v(ndiv+1),vnew(ndiv+1)

    ! EXTERNAL FUNCTIONS
    real(kind=8), external :: drand

    if ( ndiv .lt. 2 ) then
       write(*,*) 'In saintve32, ndiv must be at least 2!'
       stop
    end if

    if ( modo .eq. 1 .and. ( present(sigma) .or. present(seed) ) ) then
       write(*,*) 'In saintve32, modo=1 and sigma and/or seed are present!'
       stop
    end if

    if ( modo .eq. 1 ) then
       localseed = 123456.0d0
    end if
    
!!$    if ( modo .eq. 2 ) then
!!$       write(*,*) 'Calling saintve32 in mode=2 with ', &
!!$            'present(sigma) and present(seed) equal to: ',present(sigma),present(seed)
!!$    end if

!!$    if ( modo .eq. 2 ) then
!!$       write(*,*) 'ndiv = ',ndiv
!!$       write(*,*) 'Manning (first 5 components) = ',mann(1:min(5,ndiv))
!!$       write(*,*) 'nobs = ',nobs
!!$       write(*,*) 'losty  (first 5 components) = ',losty(1:min(5,nobs))
!!$       write(*,*) 'lefty  (first 5 components) = ',lefty(1:min(5,nobs))
!!$       write(*,*) 'righy  (first 5 components) = ',righy(1:min(5,nobs))
!!$       write(*,*) 'dt = ',dt
!!$       write(*,*) 'dx = ',dx
!!$       if ( present(sigma) ) then
!!$          write(*,*) 'sigma = ',sigma
!!$       else
!!$          write(*,*) 'sigma is NOT present'
!!$       end if
!!$       if ( present(seed) ) then
!!$          write(*,*) 'seed = ',seed
!!$       else
!!$          write(*,*) 'seed is NOT present'
!!$       end if
!!$    end if
    
    ! alfa = difusao artificial tipo Lax-Friedrish
    ! alfa = 1 é sem difusao
    ! alfa = 0 é com difusao maxima

    alfa = 0.9d0
    ualfa = 1.d0 - alfa

    grav = 9.8d0
    zbx = - 1.d-3
    width = 5.d0
    avinic = 8.245d0 
    hinic = 1.2d0
    ainic = width * hinic
    vinic = avinic / ainic
    pinic = width + 2.d0 * hinic
    rhinic = ainic / pinic
    zxinic = zbx      

    npun = ndiv + 1

    if ( modo .eq. 1 ) then
       ! g A \zx  + \frac{\f P V |V|}{8} = 0.
       efe = - grav * ainic * zxinic * 8.d0 / ( pinic * vinic * abs(vinic) )
       mann(1:ndiv) = efe
       write(*,*) '(mode=1) Manning = ',efe
!!$    else
!!$       write(*,*) '(mode=2) Manning (first 5 components) = ',mann(1:min(5,ndiv))
    endif

    a(1:npun) = ainic
    v(1:npun) = vinic
    h(1:npun) = hinic
    av(1:npun) = avinic

    k = 0

    t = 0.0d0
    
100 continue

    if ( k .ge. nobs ) then
       return
    end if
    
    t = t + dt

    ! Obter o estado do sistema em t + dt
    ! As variaveis de estado sao a(x, t) e v(x, t) (associadas a av(x, t) e h(x, t))
    ! Respectivamente, area, velocidade, vazao, profundidade.

    ! Calculos para os pontos interiores
    
    do i = 2,npun - 1
       ! Considerar a equacao de conservacao de massa 

       avx = ( a(i+1) * v(i+1) - a(i-1) * v(i-1) ) / ( 2.d0 * dx )
       at = - avx

       adifu = ualfa * ( a(i-1) + a(i+1) ) / 2.d0 + alfa * a(i)
       anew(i) = adifu + dt * at
       if ( modo .eq. 2 .and. present(sigma) .and. present(seed) ) call rterm(anew(i),sigma,seed)
 
       ! Considerar a equacao de conservacao de quantidade de movimento
       
       av2x = ( a(i+1) * v(i+1) ** 2 - a(i-1) * v(i-1) ** 2 ) / ( 2.d0 * dx ) 

       zhatx = zbx + ( h(i+1) - h(i-1) ) / ( 2.d0 * dx )
       zhatx = zhatx / ( 1.d0 + zhatx ** 2 )

       hdifu = ualfa * ( h(i-1) + h(i+1) ) / 2.d0 + alfa * h(i) 
       pdifu = width + 2.d0 * hdifu

       vdifu = ualfa * ( v(i-1) + v(i+1) ) / 2.d0 + alfa * v(i)   

       avt = -av2x - grav * adifu * zhatx - mann(i-1) * pdifu * vdifu * abs(vdifu) / 8.d0 
     ! avt = -av2x - grav * adifu * zhatx - mann(i-1) ** 2 * pdifu * vdifu * abs(vdifu) / 8.d0 
     ! avt = -av2x - grav * adifu * zhatx - abs(mann(i-1)) * pdifu * vdifu * abs(vdifu) / 8.d0 
     ! avt = -av2x - grav * adifu * zhatx - max(0.0d0,mann(i-1)) * pdifu * vdifu * abs(vdifu) / 8.d0 

       avdifu = ualfa * ( av(i-1) + av(i+1) ) / 2.d0 + alfa * av(i)   
 
       avnew(i) = avdifu + dt * avt
       if ( modo .eq. 2 .and. present(sigma) .and. present(seed) ) call rterm(avnew(i),sigma,seed)

       vnew(i) = avnew(i) / anew(i)
      
       hnew(i) = anew(i) / width
    end do

    ! Calcular as variaveis de estado no contorno da esquerda
    
    call contorno(t,avnew(1))
    
    avx = ( avnew(2) - avnew(1) ) / dx
    at = - avx
    adifu = ualfa * ( a(1) + a(2) ) / 2.d0 + alfa * a(1)      
    anew(1) = adifu + dt * at
    if ( modo .eq. 2 .and. present(sigma) .and. present(seed) ) call rterm(anew(1),sigma,seed)
    vnew(1) = avnew(1) / anew(1)
    hnew(1) = anew(1) / width  

    ! Calcular as variaveis de estado no contorno da direita
    
    avx = ( av(npun) - av(npun-1) ) / dx
    at = - avx
    adifu = ualfa * ( a(npun-1) + a(npun) ) / 2.d0 + alfa * a(npun)      
    anew(npun) = adifu + dt * at
    if ( modo .eq. 2 .and. present(sigma) .and. present(seed) ) call rterm(anew(npun),sigma,seed)
    hnew(npun) = anew(npun) / width
     
    av2x = ( a(npun) * v(npun) ** 2 - a(npun-1) * v(npun-1) ** 2 ) / dx 

    zhatx = zbx + ( h(npun) - h(npun-1) ) / dx
    zhatx = zhatx / ( 1.d0 + zhatx ** 2 )

    hdifu = ualfa * ( h(npun-1) + h(npun) ) / 2.d0 + alfa * h(npun) 
    pdifu = width + 2.d0 * hdifu

    vdifu = ualfa * ( v(npun-1) + v(npun) ) / 2.d0 + alfa * v(npun)   

    avt = - av2x - grav * adifu * zhatx - mann(npun-1) * pdifu * vdifu * abs(vdifu) / 8.d0 
  ! avt = - av2x - grav * adifu * zhatx - mann(npun-1) ** 2 * pdifu * vdifu * abs(vdifu) / 8.d0 
  ! avt = - av2x - grav * adifu * zhatx - abs(mann(npun-1)) * pdifu * vdifu * abs(vdifu) / 8.d0 
  ! avt = - av2x - grav * adifu * zhatx - max(0.0d0,mann(npun-1)) * pdifu * vdifu * abs(vdifu) / 8.d0 

    avdifu = ualfa * ( av(npun-1) + av(npun) ) / 2.d0 + alfa * av(npun)   
 
    avnew(npun) = avdifu + dt * avt
    if ( modo .eq. 2 .and. present(sigma) .and. present(seed) ) call rterm(avnew(npun),sigma,seed)

    vnew(npun) = avnew(npun) / anew(npun)

    ! Atualizar as variáveis de estado

    do i = 1,npun
       a(i) = anew(i)
       v(i) = vnew(i)
       h(i) = hnew(i)
       av(i) = avnew(i)

       y(k+1:k+2) = (/ anew(i), vnew(i) /)

       if ( modo .eq. 1 ) then
          losty(k+1:k+2) = .true.

          if ( drand(localseed) .le. 0.1d0 ) then
             losty(k+2) = .false.
             write(*,*) 'Velocity at x = ',( i - 1 ) * dx,' and t = ',t,': ',vnew(i),' m/seg.'
             
110          continue
             y(k+2) = vnew(i) + 0.1d0 * ( 2.0d0 * drand(localseed) - 1.0d0 ) * abs( vnew(i) )
             lefty(k+2) = y(k+2) - 0.1d0 * max( 1.0d0, abs( y(k+2) ) )
             righy(k+2) = y(k+2) + 0.1d0 * max( 1.0d0, abs( y(k+2) ) )
             if ( vnew(i) .lt. lefty(k+2) .or. vnew(i) .gt. righy(k+2) ) then
                go to 110
             end if
             
             write(*,*) 'Velocity was replaced with ',y(k+2)
          end if
          
!!$          if ( max( abs( ( i - 1 ) * dx - 1500.0d0 ), abs( t - 1200.0d0 ) ) .le. 1.0d-08 .or. &
!!$               max( abs( ( i - 1 ) * dx -  720.0d0 ), abs( t -  600.0d0 ) ) .le. 1.0d-08 ) then
!!$             write(*,*) 'Velocity at x = ',( i - 1 ) * dx,' and t = ',t,': ',vnew(i),' m/seg.'
!!$             if ( max( abs( ( i - 1 ) * dx -  720.0d0 ), abs( t -  600.0d0 ) ) .le. 1.0d-08 ) then
!!$                write(*,*) 'Velocity at x = ',( i - 1 ) * dx,' and t = ',t,' replaced with 1.2.'
!!$                y(k+2) = 1.2d0
!!$             end if
!!$             if ( max( abs( ( i - 1 ) * dx - 1500.0d0 ), abs( t - 1200.0d0 ) ) .le. 1.0d-08 ) then
!!$                write(*,*) 'Velocity at x = ',( i - 1 ) * dx,' and t = ',t,' replaced with 1.6.'
!!$                y(k+2) = 1.6d0
!!$             end if
!!$             losty(k+2) = .false.
!!$             lefty(k+2) = y(k+2) - 0.1d0 * max( 1.0d0, abs( y(k+2) ) ) ! This because the observation is a velocity
!!$             righy(k+2) = y(k+2) + 0.1d0 * max( 1.0d0, abs( y(k+2) ) ) ! This because the observation is a velocity
!!$          end if
       end if
       
       k = k + 2
    end do

    go to 100

  end subroutine saintve32
    
  ! ******************************************************************
  ! ******************************************************************

  subroutine contorno(t,q)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: t
    real(kind=8), intent(out) :: q

    ! PARAMETERS
    real(kind=8) :: pi = 4.0d0 * atan( 1.0d0 )
    
    q = 8.245d0 + 0.245d0 * sin( 2.0d0 * pi * t / 3600.0d0 )

  end subroutine contorno

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

  subroutine rterm(v,sigma,seed)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: sigma
    real(kind=8), intent(inout) :: v,seed
    
    ! EXTERNAL FUNCTIONS
    real(kind=8), external :: drand

    v = v + abs(v) * ( 2.0d0 * drand(seed) - 1.0d0 ) * sigma
    
  end subroutine rterm
  
end module bmproblem
