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
     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) = 60.0d0
!    pdata%tobs(1:14) = 60.0d0
!    pdata%tobs(1:5) = 720.0d0
!    pdata%tobs( 1) = 720.0d0
!    pdata%tobs( 2) = 720.0d0
!    pdata%tobs( 3) = 720.0d0
!    pdata%tobs( 4) = 720.0d0
!    pdata%tobs( 5) = 720.0d0
!    
      pdata%xobs( 1) = 0.30 
      pdata%xobs( 2) = 0.61 
      pdata%xobs( 3) = 0.91 
      pdata%xobs( 4) = 1.22 
!      pdata%xobs( 5) = 1.52 
!      pdata%xobs( 6) = 1.83
!      pdata%xobs( 7) = 2.13
!      pdata%xobs( 8) = 2.44
      pdata%xobs( 5) = 2.74
      pdata%xobs( 6) = 3.05
      pdata%xobs( 7) = 3.35
      pdata%xobs( 8) = 3.66
      pdata%xobs( 9) = 3.96
      pdata%xobs(10) = 4.27
!
!     t=720s
!    
!    pdata%yobs( 1) = 5.7d0	
!    pdata%yobs( 2) = 4.5d0
!    pdata%yobs( 3) = 3.3d0	
!    pdata%yobs( 4) = 2.2d0
!    pdata%yobs( 5) = 1.4d0
!
      pdata%yobs( 1) =  0.043	
      pdata%yobs( 2) =  0.043	
      pdata%yobs( 3) =  0.046	
      pdata%yobs( 4) =  0.049	
!      pdata%yobs( 5) =  0.055	
!      pdata%yobs( 6) =  0.152
!      pdata%yobs( 7) =  0.207
!      pdata%yobs( 8) =  0.232
      pdata%yobs( 5) =  0.229
      pdata%yobs( 6) =  0.226
      pdata%yobs( 7) =  0.223
      pdata%yobs( 8) =  0.223
      pdata%yobs( 9) =  0.223
      pdata%yobs(10) =  0.223

!!
!      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) )
    pdata%treali = maxval( pdata%tobs(1:pdata%nobs) )
    
  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=8) :: 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%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.d0 * atan( 1.d0 ) ! 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%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 saintve38(ndiv,xmin,xmax,dt,mann,sigma,seed,nobs,tobs,xobs,hcal,x,treali,hreali)
    
    !   Problema jump, figura 6 de ying2004
    !   Canal retangular 
    !   Objetivo: simular os resultados de primeiro exemplo enviado por Diego.  
    !   Método explícito
    !   Condicao de contorno dada em Q(0, t) 
    !   Condição inicial fixa
    !   Outros contornos: esquerda forward, direita backward 
    !   10 de junho de 2021 

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: ndiv,nobs
    real(kind=8), intent(in) :: dt,sigma,treali,xmin,xmax
    real(kind=8), intent(inout) :: seed
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: mann(ndiv+1),tobs(nobs),xobs(nobs)
    real(kind=8), intent(out) :: hcal(nobs),hreali(ndiv+1),x(ndiv+1)

    ! LOCAL SCALARS
    integer :: i,iobs,itime,npun
    real(kind=8) :: width,ainic,vinic,avinic,hinic,grav,pinic, &
         rhinic,zbx,zxinic,t,avx,at,alfa,ualfa,av2x,zhatx,adifu,pdifu, &
         avt,vdifu,avdifu,hdifu,efe,dx

    ! LOCAL ARRAYS
    real(kind=8) :: a(ndiv+1),v(ndiv+1),av(ndiv+1),h(ndiv+1),anew(ndiv+1), &
         vnew(ndiv+1),avnew(ndiv+1),hnew(ndiv+1)

    ! discretization of x at which quantities are being calculated

    dx = ( xmax - xmin ) / ndiv
    x(1:ndiv+1) = (/ ( xmin + ( i - 1 ) * dx, i=1,ndiv+1 ) /)
    
    !  alfa = difusao artificial tipo Lax-Friedrish
    !  alfa = 1 é sem difusao
    !  alfa = 0 é com difusao maxima

    alfa = 0.99d0
    ualfa = 1.0d0 - alfa
    width = 0.46d0

    grav = 9.8d0

    ! zbx nao é mencionado em Wiggert. Diego conjetura que é zero.
    zbx = 0.d0 

    ! profundidad inicial constante igual a 0.043
    hinic = 0.043d0

    ! velocidad inicial constante igual a constante 
    vinic = 2.737d0
    ainic = width * hinic
    pinic = width + 2.0d0 * hinic
    rhinic = ainic / pinic
    zxinic = zbx
    avinic = ainic * vinic
 
    ! g A \zx + \frac{\f P V |V|}{8} = 0.

    npun = ndiv + 1

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

    itime = 0
    t = 0.d0

    !  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.

    iobs = 1
    
100 continue

    itime = itime + 1
    t = itime * dt

    ! Calculos para os pontos interiores e tambem no contorno

    ! i = 1
    
    avx = ( a(2) * v(2) - a(1) * v(1) ) / dx
    at = - avx
    adifu = ualfa * 0.5d0 * ( a(1) + a(2) ) + alfa * a(1)       
    anew(1) = adifu + dt * at
    av2x = ( a(2) * v(2) ** 2 - a(1) * v(1) ** 2 ) / dx
    zhatx = zbx + ( h(2) - h(1) ) / dx
    zhatx = zhatx / ( 1.0d0 + zhatx ** 2 )
    hdifu = ualfa * 0.5d0 * ( h(1) + h(2) ) + alfa * h(1)
    pdifu = width + 2.0d0 * hdifu
    vdifu = ualfa * 0.5d0 * ( v(1) + v(2) ) + alfa * v(1)
    efe = 8.0d0 * mann(1) ** 2 * grav * ( pdifu / adifu ) ** ( 1.d0 / 3.d0 ) 
    avt = - av2x - grav * adifu * zhatx - efe * pdifu * vdifu * 0.125d0 * abs(vdifu) 
    avdifu = ualfa * 0.5d0 * ( av(1) + av(2) ) + alfa * av(1)        
    avnew(1) = avdifu + dt * avt
    vnew(1) = avnew(1) / anew(1)
    hnew(1) = anew(1) / width
       
    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.0d0 * dx )
       at = - avx
       adifu = ualfa * 0.5d0 * ( a(i-1) + a(i+1) ) + alfa * a(i)
       anew(i) = adifu + dt * at
       if ( sigma .gt. 0.0d0 ) call rterm(anew(i),sigma*dt,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.0d0 * dx )
       zhatx = zbx + ( h(i+1) - h(i-1) ) / ( 2.0d0 * dx )
       zhatx = zhatx / ( 1.0d0 + zhatx ** 2 )
       hdifu = ualfa * 0.5d0 * ( h(i-1) + h(i+1) ) + alfa * h(i)
       pdifu = width + 2.0d0 * hdifu
       vdifu = ualfa * 0.5d0 * ( v(i-1) + v(i+1) ) + alfa * v(i)  
       efe = 8.0d0 * mann(i) ** 2 * grav * ( pdifu / adifu ) ** ( 1.d0 / 3.d0 ) 
       avt = - av2x - grav * adifu * zhatx - efe * pdifu * vdifu * 0.125d0 * abs(vdifu)
       avdifu = ualfa * 0.5d0 * ( av(i-1) + av(i+1) ) + alfa * av(i)
       avnew(i) = avdifu + dt * avt
       if ( sigma .gt. 0.0d0 ) call rterm(avnew(i),sigma*dt,seed)
       vnew(i) = avnew(i) / anew(i)
       hnew(i) = anew(i) / width
    end do

    ! i = npun
   
    avx = ( a(npun) * v(npun) - a(npun-1) * v(npun-1) ) / dx
    at = - avx
    adifu = ualfa * 0.5d0 * ( a(npun-1) + a(npun) ) + alfa * a(npun)  
    anew(npun) = adifu + dt * at
    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.0d0 + zhatx ** 2 )
    hdifu = ualfa * 0.5d0 * ( h(npun-1) + h(npun) ) + alfa * h(npun)
    pdifu = width + 2.d0 * hdifu
    vdifu = ualfa * 0.5d0 * ( v(npun-1) + v(npun) ) + alfa * v(npun)      
    efe = 8.0d0 * mann(npun) ** 2 * grav * ( pdifu / adifu ) ** ( 1.d0 / 3.d0 ) 
    avt = - av2x - grav * adifu * zhatx - efe * pdifu * vdifu * 0.125d0 * abs(vdifu) 
    avdifu = ualfa * 0.5d0 * ( av(npun-1) + av(npun) ) + alfa * av(npun)   
    avnew(npun) = avdifu + dt * avt
    vnew(npun) = avnew(npun) / anew(npun)
    hnew(npun) = anew(npun) / width
    
    !  Modificar as variaveis de estado no contorno a esquerda
    call contiz(hnew(1),vnew(1))
    anew(1) = hnew(1) * width
    avnew(1) = anew(1) * vnew(1)
    !  a direita
    call contde(t,hnew(npun))
    anew(npun) = hnew(npun) * width
    vnew(npun) = avnew(npun) / anew(npun)
       
    !  Atualizar as variáveis de estado
    a(1:npun) = anew(1:npun)
    v(1:npun) = vnew(1:npun)
    h(1:npun) = hnew(1:npun)
    av(1:npun) = avnew(1:npun)

    if ( abs( t - treali ) .le. 0.5d0 * dt ) then
       hreali(1:ndiv+1) = h(1:ndiv+1)
    end if
    
200 continue
    
    if ( abs( t - tobs(iobs) ) .le. 0.5d0 * dt ) then
       call computehatxobs(ndiv,x,h,xobs(iobs),hcal(iobs))

       if ( iobs .eq. nobs )  then
          return
       end if
       
       iobs = iobs + 1
       go to 200
    end if
    
    go to 100
    
 end subroutine saintve38

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

 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

   ! 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
