module bmproblem

  use iso_c_binding, only: c_ptr,c_f_pointer
  
  implicit none
  
  type :: pdata_type
     integer :: ndiv,nobs
     real(kind=8) :: dt,dx,problost,tmax,xmax
     logical, allocatable :: lostobs(:)
     real(kind=8), allocatable :: xsol(:),y(:),yobs(:)
  end type pdata_type

  public :: iniprob,endprob,inip,evalr

contains

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

  subroutine iniprob(n,m,pdataptr)

    implicit none

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

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

    call c_f_pointer(pdataptr,pdata)

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

    pdata%dt = 0.1d0
    pdata%dx = 6.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%y(pdata%nobs),pdata%yobs(pdata%nobs),pdata%lostobs(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%problost,pdata%lostobs,pdata%dx,pdata%dt,1)

    n = pdata%ndiv
    m = count( .not. pdata%lostobs(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%xsol,pdata%y,pdata%yobs,pdata%lostobs,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
    type(pdata_type), pointer :: pdata

    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
  
  end subroutine inip

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

  subroutine evalr(n,x,m,r,pdataptr)

    implicit none

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

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

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

    ! LOCAL ARRAY
    real(kind=8) :: xtmp(n)
    
    call c_f_pointer(pdataptr,pdata)

    xtmp(1:n) = x(1:n) ! this is because x is 'in' here but it is 'inout' in saintve32
    
    call saintve32(pdata%ndiv,xtmp,pdata%nobs,pdata%y,pdata%problost,pdata%lostobs,pdata%dx,pdata%dt,2)
    
    r(1:m) = pack( pdata%y(1:pdata%nobs) - pdata%yobs(1:pdata%nobs), .not. pdata%lostobs(1:pdata%nobs) )
    
  end subroutine evalr
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine saintve32(ndiv,mann,nobs,y,problost,losty,dx,dt,modo)
    
    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: modo,ndiv,nobs
    real(kind=8), intent(in) :: dt,dx,problost
    
    ! ARRAY ARGUMENTS
    logical, intent(out) :: losty(nobs) ! out when modo=1, unused when modo=2
    real(kind=8), intent(inout) :: mann(ndiv) ! 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,ttop,tfin,qmin,qmax,avx,at,alfa,ualfa,av2x,zhatx,adifu, &
         pdifu,avt,vdifu,avdifu,hdifu,p1,p2,p3,per,seed
    
    ! 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
    
    ! alfa = difusao artificial tipo Lax-Friedrish
    ! alfa = 1 é sem difusao
    ! alfa = 0 é com difusao maxima

    alfa = 0.9d0
    ualfa = 1.d0 - alfa

    ttop = 1200.0d0
    tfin = 3600.0d0
    qmax =  200.0d0
    grav = 9.8d0
    zbx = - 1.d-3
    width = 5.d0
    avinic = 8.245d0 
    qmin = avinic
    hinic = 1.2d0
    ainic = width * hinic
    vinic = avinic / ainic
    pinic = width + 2.d0 * hinic
    rhinic = ainic / pinic
    zxinic = zbx      

    npun = ndiv + 1

    seed = 9388847237.d0

    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) )

!!$       write(*, *) 'f de Manning = ',efe

       p1 = 0.20d0 * ndiv
       p2 = 0.25d0 * ndiv
       p3 = ( p1 + p2 ) / 2.0d0

       per = 0.01d0
       do i = 1,ndiv
          if ( i .le. p1 ) then
             mann(i) = efe
          else if ( i .ge. p1 .and. i .le. p3 ) then
             mann(i) = efe * ( 1.0d0 + ( i - p1 ) / ( p3 - p1 ) )
          else if ( i .ge. p3 .and. i .le. p2 ) then
             mann(i) = efe * ( 1.0d0 + ( i - p2 ) / ( p3 - p2 ) )
          else
             mann(i) = efe
          end if
        ! mann(i) = mann(i) + ( 2.0d0 * drand(seed) - 1.0d0 ) * per * max( sqrt( epsilon( 1.0d0 ) ), abs( efe ) )
          mann(i) = sqrt( mann(i) + ( 2.0d0 * drand(seed) - 1.0d0 ) * per * max( sqrt( epsilon( 1.0d0 ) ), abs( efe ) ) )
       end do
    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
 
       ! 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

       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,ttop,tfin,qmin,qmax,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
    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
    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

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

    ! Atualizar as variáveis de estado

    do i = 1,npun
       a(i) = anew(i)
       
       k = k + 1

       y(k) = a(i)

       if ( modo .eq. 1 ) then
          if ( drand(seed) .le. problost ) then
             losty(k) = .true.
          else
             losty(k) = .false.
          end if
       end if
       
       v(i) = vnew(i)
       
       k = k + 1
       
       y(k) = vnew(i)

       if ( modo .eq. 1 ) then
          if ( drand(seed) .le. problost ) then
             losty(k) = .true.
          else
             losty(k) = .false.
          end if
       end if
       
       h(i) = hnew(i)
       
       av(i) = avnew(i)
    end do

    go to 100

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

  subroutine contorno(t,ttop,tfin,qmin,qmax,q)

    implicit none

    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: t,tfin,ttop,qmax,qmin
    real(kind=8), intent(out) :: q
    
    if ( t .le. ttop ) then
       q = qmin + t * ( qmax - qmin ) / ttop

    else if ( t .gt. tfin ) then
       q = qmin

    else 
       q = qmax + ( t - ttop ) * ( qmin - qmax ) / ( tfin - ttop )
    end if

  end subroutine contorno

end module bmproblem