module bmproblem

  use iso_c_binding, only: c_ptr,c_f_pointer
  
  implicit none
  
  type :: pdata_type
     integer :: ndiv,nobs,nr
     integer :: nshydro,nthydro
     real(kind=8) :: dt,sigma,dtsigma,treali,xmax,xmin,tmax
     real(kind=8), allocatable :: mann(:),xobs(:),tobs(:),yobs(:),ycal(:),x(:),reali(:,:)
     real(kind=8), allocatable :: hydro(:,:,:),xhydro(:),thydro(:)
  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
!
    real(8) :: tmin,tmax
!
    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 = 22
!
!     para imprimir hidrogramas
!
!     numero de estacoes (inicialmente so vou plotar uma)
!
    pdata%nshydro =3
!
!     quantidade de dados no tempo   
!
    pdata%nthydro=1000
!   
    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), &
         pdata%hydro(pdata%nthydro,pdata%nshydro,pdata%nr),pdata%xhydro(pdata%nshydro), &
         pdata%thydro(pdata%nthydro),stat=allocstat)

    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate.'
       stop
    end if
!    
!     hidrograma: posicao
!
      pdata%xhydro(1) = 0.0
      pdata%xhydro(2) = 2505.0
      pdata%xhydro(3) = 3295.0
!
!    pdata%tobs(1:pdata%nobs) = 240.0d0
!   
      pdata%tobs( 1) =   24.00*60.d0*60.d0
      pdata%tobs( 2) =   36.00*60.d0*60.d0
      pdata%tobs( 3) =   48.00*60.d0*60.d0
      pdata%tobs( 4) =   60.00*60.d0*60.d0
      pdata%tobs( 5) =   72.00*60.d0*60.d0
      pdata%tobs( 6) =   84.00*60.d0*60.d0
      pdata%tobs( 7) =   96.00*60.d0*60.d0
      pdata%tobs( 8) =  108.00*60.d0*60.d0
      pdata%tobs( 9) =  120.00*60.d0*60.d0
      pdata%tobs(10) =  132.00*60.d0*60.d0

      pdata%tobs(11) = 144.00*60.d0*60.d0
      pdata%tobs(12) = 156.00*60.d0*60.d0
      pdata%tobs(13) = 168.00*60.d0*60.d0
      pdata%tobs(14) = 180.00*60.d0*60.d0
      pdata%tobs(15) = 192.00*60.d0*60.d0
      pdata%tobs(16) = 204.00*60.d0*60.d0
      pdata%tobs(17) = 216.00*60.d0*60.d0
      pdata%tobs(18) = 228.00*60.d0*60.d0
      pdata%tobs(19) = 240.00*60.d0*60.d0
      pdata%tobs(20) = 252.00*60.d0*60.d0
!
      pdata%tobs(21) = 264.00*60.d0*60.d0
      pdata%tobs(22) = 276.00*60.d0*60.d0
!      pdata%tobs(23) = 288.00*60.d0*60.d0
!      pdata%tobs(24) = 300.00*60.d0*60.d0
!      pdata%tobs(25) = 312.00*60.d0*60.d0
!      pdata%tobs(26) = 324.00*60.d0*60.d0
!      pdata%tobs(27) = 336.00*60.d0*60.d0
!      pdata%tobs(28) = 348.00*60.d0*60.d0
!      pdata%tobs(29) = 360.00*60.d0*60.d0
!      pdata%tobs(30) = 372.00*60.d0*60.d0
!
!     pdata%tobs(31) =   24.00*60.d0*60.d0
!     pdata%tobs(32) =   36.00*60.d0*60.d0
!     pdata%tobs(33) =   48.00*60.d0*60.d0
!     pdata%tobs(34) =   60.00*60.d0*60.d0
!     pdata%tobs(35) =   72.00*60.d0*60.d0
!     pdata%tobs(36) =   84.00*60.d0*60.d0
!     pdata%tobs(37) =   96.00*60.d0*60.d0
!     pdata%tobs(38) =  108.00*60.d0*60.d0
!     pdata%tobs(39) =  120.00*60.d0*60.d0
!     pdata%tobs(40) =  132.00*60.d0*60.d0
!
!     pdata%tobs(41) = 144.00*60.d0*60.d0
!     pdata%tobs(42) = 156.00*60.d0*60.d0
!     pdata%tobs(43) = 168.00*60.d0*60.d0
!     pdata%tobs(44) = 180.00*60.d0*60.d0
!     pdata%tobs(45) = 192.00*60.d0*60.d0
!     pdata%tobs(46) = 204.00*60.d0*60.d0
!     pdata%tobs(47) = 216.00*60.d0*60.d0
!     pdata%tobs(48) = 228.00*60.d0*60.d0
!     pdata%tobs(49) = 240.00*60.d0*60.d0
!     pdata%tobs(50) = 252.00*60.d0*60.d0
!
!     pdata%tobs(51) = 264.00*60.d0*60.d0
!     pdata%tobs(52) = 276.00*60.d0*60.d0
!     pdata%tobs(53) = 288.00*60.d0*60.d0
!     pdata%tobs(54) = 300.00*60.d0*60.d0
!     pdata%tobs(55) = 312.00*60.d0*60.d0
!     pdata%tobs(56) = 324.00*60.d0*60.d0
!     pdata%tobs(57) = 336.00*60.d0*60.d0
!     pdata%tobs(58) = 348.00*60.d0*60.d0
!     pdata%tobs(59) = 360.00*60.d0*60.d0
!     pdata%tobs(60) = 372.00*60.d0*60.d0
!
!     xobs
!   
      pdata%xobs( 1) = 2505.d0
      pdata%xobs( 2) = 2505.d0
      pdata%xobs( 3) = 2505.d0
      pdata%xobs( 4) = 2505.d0
      pdata%xobs( 5) = 2505.d0
      pdata%xobs( 6) = 2505.d0
      pdata%xobs( 7) = 2505.d0
      pdata%xobs( 8) = 2505.d0
      pdata%xobs( 9) = 2505.d0
      pdata%xobs(10) = 2505.d0
   
      pdata%xobs(11) = 2505.d0
      pdata%xobs(12) = 2505.d0
      pdata%xobs(13) = 2505.d0
      pdata%xobs(14) = 2505.d0
      pdata%xobs(15) = 2505.d0
      pdata%xobs(16) = 2505.d0
      pdata%xobs(17) = 2505.d0
      pdata%xobs(18) = 2505.d0
      pdata%xobs(19) = 2505.d0
      pdata%xobs(20) = 2505.d0
!!
      pdata%xobs(21) = 2505.d0
      pdata%xobs(22) = 2505.d0
!      pdata%xobs(23) = 2505.d0
!      pdata%xobs(24) = 2505.d0
!      pdata%xobs(25) = 2505.d0
!      pdata%xobs(26) = 2505.d0
!      pdata%xobs(27) = 2505.d0
!      pdata%xobs(28) = 2505.d0
 !    pdata%xobs(29) = 2505.d0
 !    pdata%xobs(30) = 2505.d0
!
!     pdata%xobs(31) = 0.d0
!     pdata%xobs(32) = 0.d0
!     pdata%xobs(33) = 0.d0
!     pdata%xobs(34) = 0.d0
!     pdata%xobs(35) = 0.d0
!     pdata%xobs(36) = 0.d0
!     pdata%xobs(37) = 0.d0
!     pdata%xobs(38) = 0.d0
!     pdata%xobs(39) = 0.d0
!     pdata%xobs(40) = 0.d0
!   
!     pdata%xobs(41) = 0.d0
!     pdata%xobs(42) = 0.d0
!     pdata%xobs(43) = 0.d0
!     pdata%xobs(44) = 0.d0
!     pdata%xobs(45) = 0.d0
!     pdata%xobs(46) = 0.d0
!     pdata%xobs(47) = 0.d0
!     pdata%xobs(48) = 0.d0
!     pdata%xobs(49) = 0.d0
!     pdata%xobs(50) = 0.d0
!
!     pdata%xobs(51) = 0.d0
!     pdata%xobs(52) = 0.d0
!     pdata%xobs(53) = 0.d0
!     pdata%xobs(54) = 0.d0
!     pdata%xobs(55) = 0.d0
!     pdata%xobs(56) = 0.d0
!     pdata%xobs(57) = 0.d0
!     pdata%xobs(58) = 0.d0
!     pdata%xobs(59) = 0.d0
!     pdata%xobs(60) = 0.d0
!
!     yobs
!
      pdata%yobs( 1) = 7.93
      pdata%yobs( 2) = 7.96
      pdata%yobs( 3) = 8.03
      pdata%yobs( 4) = 8.15
      pdata%yobs( 5) = 8.10
      pdata%yobs( 6) = 8.21
      pdata%yobs( 7) = 8.14
      pdata%yobs( 8) = 8.26
      pdata%yobs( 9) = 8.19
      pdata%yobs(10) = 8.25
!   
      pdata%yobs(11) = 8.20
      pdata%yobs(12) = 8.27
      pdata%yobs(13) = 8.25
      pdata%yobs(14) = 8.38
      pdata%yobs(15) = 8.31
      pdata%yobs(16) = 8.33
      pdata%yobs(17) = 8.26
      pdata%yobs(18) = 8.30
      pdata%yobs(19) = 8.29
      pdata%yobs(20) = 8.09
!
      pdata%yobs(21) = 7.91
      pdata%yobs(22) = 7.71
 !     pdata%yobs(23) = 7.61
!      pdata%yobs(24) = 7.46
!      pdata%yobs(25) = 7.43
!      pdata%yobs(26) = 7.38
!      pdata%yobs(27) = 7.39
!      pdata%yobs(28) = 7.52
  !   pdata%yobs(29) = 7.53
  !   pdata%yobs(30) = 7.77 
!
!     pdata%yobs(31) = 5.99
!     pdata%yobs(32) = 6.09
!     pdata%yobs(33) = 6.09
!     pdata%yobs(34) = 6.28
!     pdata%yobs(35) = 6.18
!     pdata%yobs(36) = 6.36
!     pdata%yobs(37) = 6.23
!     pdata%yobs(38) = 6.43
!     pdata%yobs(39) = 6.30
!     pdata%yobs(40) = 6.39
!   
!     pdata%yobs(41) = 6.31
!     pdata%yobs(42) = 6.41
!     pdata%yobs(43) = 6.37
!     pdata%yobs(44) = 6.61
!     pdata%yobs(45) = 6.48
!     pdata%yobs(46) = 6.53
!     pdata%yobs(47) = 6.39
!     pdata%yobs(48) = 6.46
!     pdata%yobs(49) = 6.45
!     pdata%yobs(50) = 6.17
!
!     pdata%yobs(51) = 6.00
!     pdata%yobs(52) = 5.82
!     pdata%yobs(53) = 5.73
!     pdata%yobs(54) = 5.63
!     pdata%yobs(55) = 5.58
!     pdata%yobs(56) = 5.55
!     pdata%yobs(57) = 5.52
!     pdata%yobs(58) = 5.64
!     pdata%yobs(59) = 5.65
!     pdata%yobs(60) = 5.84
!
!
!      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) )
!      pdata%tmax=324.d0*60.d0*60.d0
!
!     ate 1Pm de 31 de maio
!
      pdata%tmax=276.d0*60.d0*60.d0
!
!   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%treali=276.d0*60.d0*60.d0
!
!      pdata%tmax    = 3480.0
!      pdata%treali  = pdata%tmax 
!
!    pdata%treali = 3000.d0
!
!     impressao dos hidrogramas
!
      tmin=12.d0*60.d0*60.d0!minval( pdata%tobs(1:pdata%nobs) )
      tmax=pdata%tmax!maxval( pdata%tobs(1:pdata%nobs) )
!
      do i=1,pdata%nthydro
      pdata%thydro(i)=tmin+(tmax-tmin)/(pdata%nthydro-1)*(i-1)
!      write(*,*) pdata%thydro(i)/60.d0/60.d0/24.d0
      end do ! i
!      stop
!
!
  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,pdata%hydro,pdata%thydro,pdata%xhydro,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,i
    real(kind=4) :: seed
    type(pdata_type), pointer :: pdata
    real(8) :: di,dd,ai

    call c_f_pointer(pdataptr,pdata)

    if ( pdata%sigma .eq. 0.0d0 ) then
       call dgsv(pdata%nshydro,pdata%nthydro, &
            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),pdata%hydro(1:pdata%nthydro,1:pdata%nshydro,1), &
            pdata%xhydro(1:pdata%nshydro),pdata%thydro(1:pdata%nthydro))
!
!     corrigido
!
!       v = exp( - sqrt(sum( ( pdata%ycal(1:pdata%nobs) - pdata%yobs(1:pdata%nobs) ) ** 2.d0 ) / &
!            ( pdata%nobs * vartheta ** 2.d0 )) )
!
!     Huber Loss function
!
      dd=0.d0
      do i=1, pdata%nobs
      ai=dabs(pdata%ycal(i)-pdata%yobs(i))/vartheta
!      if(ai.le.1.d0) then
!      di=0.5d0*ai**(2.d0/1.d0)
!      else
!      di=ai-0.5d0
!      end if
!      dd=dd+exp(-di)
!      dd=dd+di
      dd=dd+ai**2.d0
      end do ! i
!
     dd=dd/pdata%nobs
!      dd=sqrt(dd)
!      dd=(dd)
!
!      v = dd!exp( -dd)
      v = exp( -dd)
!
    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%nshydro,pdata%nthydro, &
               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),pdata%hydro(1:pdata%nthydro,1:pdata%nshydro,j), &
               pdata%xhydro(1:pdata%nshydro),pdata%thydro(1:pdata%nthydro))
!
!     corrigido
!
!       v = v+ exp( - sqrt(sum( ( pdata%ycal(1:pdata%nobs) - pdata%yobs(1:pdata%nobs) ) ** 2.d0 ) / &
!            ( pdata%nobs * vartheta ** 2.d0 )) )
!
!          v = v + exp( - sum( ( pdata%ycal(1:pdata%nobs) - pdata%yobs(1:pdata%nobs) ) ** 2.d0 ) / &
!               ( pdata%nobs * vartheta ** 2.d0 ) )
!
!     Huber Loss function
!
      dd=0.d0
      do i=1, pdata%nobs
      ai=dabs(pdata%ycal(i)-pdata%yobs(i))/vartheta
!      if(ai.le.1.d0) then
!      di=0.5d0*ai**(2.d0/1.d0)
!      else
!      di=ai-0.5d0
!      end if
!      dd=dd+exp(-di)
!      dd=dd+di
      dd=dd+ai**2.d0
      end do ! i
!      dd=sqrt(dd/pdata%nobs)
      dd=(dd/pdata%nobs)
!     
!      v = v+ dd!exp( -dd)
      v = v+exp( -dd)
!
       end do
       v = v / pdata%nr
    end if
!
!      do i=1,pdata%nthydro
!      do j=1,1
!      write(*,*) pdata%hydro(i,3,j)
!      end do ! j
!      end do ! i
!      stop
!write(*,*) v

  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
