module bmproblem

  use iso_c_binding, only: c_ptr,c_f_pointer
  
  implicit none
  
  type :: pdata_type
     integer :: ndiv,nobs,nr
     real(kind=8) :: dt,sigma,dtsigma,treali,xmax,xmin,tmax
     real(kind=8), allocatable :: mann(:),xobs(:),tobs(:),yobs(:),ycal(:),x(:),reali(:,:)
  end type pdata_type

  public :: iniprob,endprob,setmann,setsigma,evallikelihood

contains

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

  subroutine iniprob(pdataptr)

    implicit none

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

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

    call c_f_pointer(pdataptr,pdata)

    pdata%nr = 100
    
    write(*,*) 'Enter xmin, xmax, ndiv: '
    read(*,*) pdata%xmin,pdata%xmax,pdata%ndiv

    write(*,*) 'Enter dt: '
    read(*,*) pdata%dt

    pdata%nobs = 10
    
    allocate(pdata%mann(pdata%ndiv+1),pdata%xobs(pdata%nobs),pdata%tobs(pdata%nobs), &
         pdata%yobs(pdata%nobs),pdata%ycal(pdata%nobs),pdata%x(pdata%ndiv+1), &
         pdata%reali(pdata%ndiv+1,pdata%nr),stat=allocstat)

    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate.'
       stop
    end if
!    
!    pdata%tobs(1:pdata%nobs) = 10.d0
!   
      pdata%tobs( 1) =  4.0
      pdata%tobs( 2) =  4.5
      pdata%tobs( 3) =  5.0
      pdata%tobs( 4) =  7.0
      pdata%tobs( 5) =  9.0
      pdata%tobs( 6) = 10.0
      pdata%tobs( 7) = 11.0
      pdata%tobs( 8) = 13.0
      pdata%tobs( 9) = 15.0
      pdata%tobs(10) = 20.0
!
!     xobs
!
      pdata%xobs( 1) = 68.58
      pdata%xobs( 2) = 68.58
      pdata%xobs( 3) = 68.58
      pdata%xobs( 4) = 68.58
      pdata%xobs( 5) = 68.58
      pdata%xobs( 6) = 68.58
      pdata%xobs( 7) = 68.58
      pdata%xobs( 8) = 68.58
      pdata%xobs( 9) = 68.58
      pdata%xobs(10) = 68.58
!
!     yobs
!
      pdata%yobs( 1) = 0.00
      pdata%yobs( 2) = 0.01
      pdata%yobs( 3) = 0.02
      pdata%yobs( 4) = 0.04
      pdata%yobs( 5) = 0.04
      pdata%yobs( 6) = 0.04
      pdata%yobs( 7) = 0.04
      pdata%yobs( 8) = 0.05
      pdata%yobs( 9) = 0.05
      pdata%yobs(10) = 0.05
!
!      pdata%yobs( 1) = 0.00000
!      pdata%yobs( 2) = 0.00914
!      pdata%yobs( 3) = 0.02438
!      pdata%yobs( 4) = 0.03658
!      pdata%yobs( 5) = 0.04267
!      pdata%yobs( 6) = 0.04267
!      pdata%yobs( 7) = 0.04267
!      pdata%yobs( 8) = 0.04572
!      pdata%yobs( 9) = 0.04572
!      pdata%yobs(10) = 0.04877
!
!     completo
!
!     pdata%tobs( 1) =  0.0
!     pdata%tobs( 2) =  0.5
!     pdata%tobs( 3) =  1.0
!     pdata%tobs( 4) =  1.5
!     pdata%tobs( 5) =  2.0
!     pdata%tobs( 6) =  2.5
!     pdata%tobs( 7) =  3.0
!     pdata%tobs( 8) =  3.5
!     pdata%tobs( 9) =  4.0
!     pdata%tobs(10) =  4.5
!     pdata%tobs(11) =  5.0
!     pdata%tobs(12) =  7.0
!     pdata%tobs(13) =  9.0
!     pdata%tobs(14) = 10.0
!     pdata%tobs(15) = 11.0
!     pdata%tobs(16) = 13.0
!     pdata%tobs(17) = 15.0
!     pdata%tobs(18) = 17.0
!     pdata%tobs(19) = 19.0
!     pdata%tobs(20) = 20.0
!     pdata%tobs(21) = 21.0
!     pdata%tobs(22) = 25.0
!     pdata%tobs(23) = 27.0
!     pdata%tobs(24) = 30.0
!     pdata%tobs(25) = 35.0
!     pdata%tobs(26) = 40.0
!     pdata%tobs(27) = 45.0
!     pdata%tobs(28) = 50.0
!     pdata%tobs(29) = 55.0
!     pdata%tobs(30) = 60.0
!     pdata%tobs(31) = 65.0
!     pdata%tobs(32) = 70.0
!     pdata%tobs(33) = 75.0
!     pdata%tobs(34) = 80.0
!     pdata%tobs(35) = 85.0
!     pdata%tobs(36) = 90.0
!     pdata%tobs(37) =100.0
!     pdata%tobs(38) =110.0
!     pdata%tobs(39) =120.0
!     pdata%tobs(40) =130.0
!     pdata%tobs(41) =140.0
!     pdata%tobs(42) =150.0
!     pdata%tobs(43) =160.0
!     pdata%tobs(44) =170.0
!     pdata%tobs(45) =180.0
!     pdata%tobs(46) =190.0
!     pdata%tobs(47) =200.0
!     pdata%tobs(48) =210.0
!     pdata%tobs(49) =240.0
!     pdata%tobs(50) =240.0
!
!     xobs
!
!     pdata%xobs( 1) = 68.58
!     pdata%xobs( 2) = 68.58
!     pdata%xobs( 3) = 68.58
!     pdata%xobs( 4) = 68.58
!     pdata%xobs( 5) = 68.58
!     pdata%xobs( 6) = 68.58
!     pdata%xobs( 7) = 68.58
!     pdata%xobs( 8) = 68.58
!     pdata%xobs( 9) = 68.58
!     pdata%xobs(10) = 68.58
!     pdata%xobs(11) = 68.58
!     pdata%xobs(12) = 68.58
!     pdata%xobs(13) = 68.58
!     pdata%xobs(14) = 68.58
!     pdata%xobs(15) = 68.58
!     pdata%xobs(16) = 68.58
!     pdata%xobs(17) = 68.58
!     pdata%xobs(18) = 68.58
!     pdata%xobs(19) = 68.58
!     pdata%xobs(20) = 68.58
!     pdata%xobs(21) = 68.58
!     pdata%xobs(22) = 68.58
!     pdata%xobs(23) = 68.58
!     pdata%xobs(24) = 68.58
!     pdata%xobs(25) = 68.58
!     pdata%xobs(26) = 68.58
!     pdata%xobs(27) = 68.58
!     pdata%xobs(28) = 68.58
!     pdata%xobs(29) = 68.58
!     pdata%xobs(30) = 68.58
!     pdata%xobs(31) = 68.58
!     pdata%xobs(32) = 68.58
!     pdata%xobs(33) = 68.58
!     pdata%xobs(34) = 68.58
!     pdata%xobs(35) = 68.58
!     pdata%xobs(36) = 68.58
!     pdata%xobs(37) = 68.58
!     pdata%xobs(38) = 68.58
!     pdata%xobs(39) = 68.58
!     pdata%xobs(40) = 68.58
!     pdata%xobs(41) = 68.58
!     pdata%xobs(42) = 68.58
!     pdata%xobs(43) = 68.58
!     pdata%xobs(44) = 68.58
!     pdata%xobs(45) = 68.58
!     pdata%xobs(46) = 68.58
!     pdata%xobs(47) = 68.58
!     pdata%xobs(48) = 68.58
!     pdata%xobs(49) = 68.58
!     pdata%xobs(50) = 68.58
!
!     yobs
!
!     pdata%yobs( 1) = 0.00000
!     pdata%yobs( 2) = 0.00000
!     pdata%yobs( 3) = 0.00000
!     pdata%yobs( 4) = 0.00000
!     pdata%yobs( 5) = 0.00000
!     pdata%yobs( 6) = 0.00000
!     pdata%yobs( 7) = 0.00000
!     pdata%yobs( 8) = 0.00000
!     pdata%yobs( 9) = 0.00000
!     pdata%yobs(10) = 0.00914
!     pdata%yobs(11) = 0.02438
!     pdata%yobs(12) = 0.03658
!     pdata%yobs(13) = 0.04267
!     pdata%yobs(14) = 0.04267
!     pdata%yobs(15) = 0.04267
!     pdata%yobs(16) = 0.04572
!     pdata%yobs(17) = 0.04572
!     pdata%yobs(18) = 0.04877
!     pdata%yobs(19) = 0.04877
!     pdata%yobs(20) = 0.04877
!     pdata%yobs(21) = 0.04877
!     pdata%yobs(22) = 0.04877
!     pdata%yobs(23) = 0.04877
!     pdata%yobs(24) = 0.04877
!     pdata%yobs(25) = 0.04877
!     pdata%yobs(26) = 0.04877
!     pdata%yobs(27) = 0.04572
!     pdata%yobs(28) = 0.04572
!     pdata%yobs(29) = 0.04267
!     pdata%yobs(30) = 0.04267
!     pdata%yobs(31) = 0.03962
!     pdata%yobs(32) = 0.03962
!     pdata%yobs(33) = 0.03962
!     pdata%yobs(34) = 0.03658
!     pdata%yobs(35) = 0.03658
!     pdata%yobs(36) = 0.03353
!     pdata%yobs(37) = 0.02438
!     pdata%yobs(38) = 0.02438
!     pdata%yobs(39) = 0.02134
!     pdata%yobs(40) = 0.02134
!     pdata%yobs(41) = 0.02134
!     pdata%yobs(42) = 0.02134
!     pdata%yobs(43) = 0.01829
!     pdata%yobs(44) = 0.01829
!     pdata%yobs(45) = 0.01829
!     pdata%yobs(46) = 0.01829
!     pdata%yobs(47) = 0.01829
!     pdata%yobs(48) = 0.01829
!     pdata%yobs(49) = 0.01524
!     pdata%yobs(50) = 0.01524
!
!      do i=1, pdata%nobs
!      write(*,*) pdata%tobs(i),pdata%xobs(i),pdata%yobs(i)
!      end do ! i
!      stop
!

    ! treali is the time for which the pdata%nr realizations will be
    ! saved to be plotted. It must be equal to pdata%tobs(i) for some
    ! i between 1 and pdata%nobs.
    
!    pdata%treali = minval( pdata%tobs(1:pdata%nobs) )
!
!   tempo maximo de simulacao
!
    pdata%tmax   = 20.d0!maxval( pdata%tobs(1:pdata%nobs) )
!
!   tempo para plotagem
!
    pdata%treali = 20.d0!maxval( pdata%tobs(1:pdata%nobs) )
!
!     para gerar as medias
!
!    pdata%treali = 3480.d0
    
  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%mann,pdata%xobs,pdata%tobs,pdata%yobs,pdata%ycal, &
         pdata%x,pdata%reali,stat=allocstat)
    
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate statement!'
       stop
    end if

  end subroutine endprob

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

  subroutine setmann(mann,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: mann
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    type(pdata_type), pointer :: pdata
  
    call c_f_pointer(pdataptr,pdata)

    pdata%mann(1:pdata%ndiv+1) = mann
        
  end subroutine setmann

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

  subroutine setsigma(sigma,dtsigma,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: sigma,dtsigma
    type(c_ptr), optional, intent(in) :: pdataptr

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

    call c_f_pointer(pdataptr,pdata)

    pdata%sigma   = sigma
    pdata%dtsigma = dtsigma

  end subroutine setsigma

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

  subroutine evallikelihood(vartheta,v,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: vartheta
    real(kind=8), intent(out) :: v
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    integer :: j
    real(kind=4) :: seed
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    if ( pdata%sigma .eq. 0.0d0 ) then
       call dgsv(pdata%ndiv,pdata%xmin,pdata%xmax,pdata%dt,pdata%mann,pdata%sigma,pdata%dtsigma, &
            seed,pdata%nobs,pdata%tobs,pdata%xobs,pdata%ycal,pdata%x,pdata%treali,pdata%tmax, &
            pdata%reali(1:pdata%ndiv+1,1))
       v = exp( - sum( ( pdata%ycal(1:pdata%nobs) - pdata%yobs(1:pdata%nobs) ) ** 2 ) / &
            ( pdata%nobs * vartheta ** 2 ) )

    else
       v = 0.0d0
       do j = 1,pdata%nr
          seed = j * 9388847237.d0! j * 4.d0 * atan( 1.d0 ) ! 7315697237.d0
!          seed =  j * 4.0 * atan( 1.0 ) ! 7315697237.d0

          call dgsv(pdata%ndiv,pdata%xmin,pdata%xmax,pdata%dt,pdata%mann,pdata%sigma,pdata%dtsigma, &
               seed,pdata%nobs,pdata%tobs,pdata%xobs,pdata%ycal,pdata%x,pdata%treali,pdata%tmax, &
               pdata%reali(1:pdata%ndiv+1,j))
          v = v + exp( - sum( ( pdata%ycal(1:pdata%nobs) - pdata%yobs(1:pdata%nobs) ) ** 2.d0 ) / &
               ( pdata%nobs * vartheta ** 2.d0 ) )
       end do
       v = v / pdata%nr
    end if

  end subroutine evallikelihood

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

 subroutine rterm(v,sigma,seed)

   implicit none

   ! SCALAR ARGUMENTS
   real(kind=8), intent(in) :: sigma
   real(kind=8), intent(inout) :: v
   real(kind=4), intent(inout) :: seed
    
!!$    ! EXTERNAL FUNCTIONS
!!$    real(kind=8), external :: drand
!!$
!!$    v = v + abs(v) * ( 2.0d0 * drand(seed) - 1.0d0 ) * sigma

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

   v = v + r8_normal_ab(0.0d0,sigma*abs(v),seed)
    
 end subroutine rterm
  
 ! ******************************************************************
 ! ******************************************************************

 subroutine contiz(h,v)
   
   implicit none

   ! SCALAR ARGUMENTS
   real(kind=8), intent(out) :: h,v 

   h = 0.043d0
   v = 2.737d0
   
 end subroutine contiz

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

 subroutine contde(t, h)

   implicit none

   ! SCALAR ARGUMENTS
   real(kind=8), intent(in) :: t
   real(kind=8), intent(out) :: h
   
   if ( t .le. 50.0d0 ) then
      h = 0.043d0 + t * ( 0.222d0 - 0.043d0 ) / 50.0d0
   else
      h = 0.222d0
   end if
   
 end subroutine contde
  
 ! ******************************************************************
 ! ******************************************************************
 
 subroutine computehatxobs(ndiv,x,h,xobs,hcal)

   implicit none

   ! SCALAR ARGUMENTS
   integer, intent(in) :: ndiv
   real(kind=8), intent(in) :: xobs
   real(kind=8), intent(out) :: hcal
    
   ! ARRAY ARGUMENTS
   real(kind=8), intent(in) :: h(ndiv+1),x(ndiv+1)
    
   ! LOCAL SCALARS
   integer :: i

   if ( xobs .lt. x(1) .or. xobs .gt. x(ndiv+1) ) then
      write(*,*) 'Observation is out of range!'
      stop
   end if

   i = 1

10 if ( i .eq. ndiv + 1 ) then
      go to 20
   end if

   if ( x(i+1) .gt. xobs ) then
      go to 20
   end if
       
   i = i + 1
   go to 10

20 continue

   if ( xobs .eq. x(i) ) then
      hcal = h(i)
   else
      hcal = h(i) + ( h(i+1) - h(i) ) * ( xobs - x(i) ) / ( x(i+1) - x(i) )
   end if
   
 end subroutine computehatxobs
  
end module bmproblem
