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 = 15
    
    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) = 240.0d0
!   
      pdata%tobs( 1) =  720.0
      pdata%tobs( 2) =  720.0
      pdata%tobs( 3) =  720.0
      pdata%tobs( 4) =  720.0
      pdata%tobs( 5) =  720.0
!   
      pdata%tobs( 6) = 1500.0
      pdata%tobs( 7) = 1500.0
      pdata%tobs( 8) = 1500.0
      pdata%tobs( 9) = 1500.0
      pdata%tobs(10) = 1500.0
!
      pdata%tobs(11) = 2220.0
      pdata%tobs(12) = 2220.0
      pdata%tobs(13) = 2220.0
      pdata%tobs(14) = 2220.0
      pdata%tobs(15) = 2220.0
!
!     xobs
!   
      pdata%xobs( 1) =  600.0
      pdata%xobs( 2) = 1200.0
      pdata%xobs( 3) = 1800.0
      pdata%xobs( 4) = 2400.0
      pdata%xobs( 5) = 2950.0
!   
      pdata%xobs( 6) =  600.0
      pdata%xobs( 7) = 1200.0
      pdata%xobs( 8) = 1800.0
      pdata%xobs( 9) = 2400.0
      pdata%xobs(10) = 2950.0   
!
      pdata%xobs(11) =  600.0
      pdata%xobs(12) = 1200.0
      pdata%xobs(13) = 1800.0
      pdata%xobs(14) = 2400.0
      pdata%xobs(15) = 2950.0
!
!     yobs
!
      pdata%yobs( 1) =  5.7
      pdata%yobs( 2) =  4.5
      pdata%yobs( 3) =  3.3
      pdata%yobs( 4) =  2.2
      pdata%yobs( 5) =  1.4
!
      pdata%yobs( 6) = 10.7
      pdata%yobs( 7) = 10.0
      pdata%yobs( 8) =  9.2
      pdata%yobs( 9) =  8.4
      pdata%yobs(10) =  7.8
!
      pdata%yobs(11) = 10.7
      pdata%yobs(12) = 10.7
      pdata%yobs(13) = 10.6
      pdata%yobs(14) = 10.6
      pdata%yobs(15) = 10.5
 
!!
!      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.
    
!   tempo maximo de simulacao
!
    pdata%tmax   = maxval( pdata%tobs(1:pdata%nobs) )
!
!   tempo para plotagem
!
!    pdata%treali = minval( pdata%tobs(1:pdata%nobs) )
!    pdata%treali = 720.d0!maxval( pdata%tobs(1:pdata%nobs) )
    pdata%treali = maxval( pdata%tobs(1:pdata%nobs) )
!
!      pdata%tmax    = 3480.0
!      pdata%treali  = pdata%tmax 
!     para gerar as medias
!
!    pdata%treali = 3000.d0
   ! pdata%treali = 3240.d0
   ! 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
