!
!.----------------------------------------------------------------------
!
      module mgeometry
      implicit none 
!
      integer, parameter :: nsd = 1 ! dimensão espacial do problema
!
      integer, parameter :: ns  = 2 ! equacoes do sistema
!
      integer, parameter :: nen = 2 ! number of element nodes 
!
      integer, parameter :: ndim=201 ! dimensao maxima dos arrays
!
      integer :: nnodes
      integer :: numel
      integer :: ndof
!
      real(8) :: xli,xlf,ylf,yli
      real(8) :: hxmin
      integer :: nx
!
      integer :: ndofg     ! number of degrees of freedom (global)
      integer, parameter :: ndofl=2     ! number of degrees of freedom (local)
!
      real(8) :: h0,b0,vup,hdown
!
      end module 
!
!----------------------------------------------------------------------
!
      module mdata
      implicit none 
!
      integer, parameter :: ndmax=1!500
      integer :: ndatasec,ndatazb,ndatah,ndataq
      integer :: ireaddata=0
      integer, parameter :: nstations=4
      integer, parameter :: nobstime =1!200
!
      end module 
!
!----------------------------------------------------------------------
!
      module mconstantes
      implicit none 
!
      real(8), parameter :: grav=9.81d0!9.8d0
      real(8), parameter :: rhow=1.d3!kg/m3
      real(8), parameter :: hdry=1.d-12
!
      end module 
!
!----------------------------------------------------------------------
!
      module mtime
      implicit none 
!
      real(8) :: ti,tf,tt,dt
      real(8) :: dtg
      real(8) :: cr ! courant
      real(8) :: alphamax
!
      real(8) :: tsigma,xsigma,dxsigma
!
      end module
!
!----------------------------------------------------------------------
!
      module marrays
!
      use mgeometry
      use mdata
!
      implicit none 
!
      integer, dimension(ndim)     :: idry
      real(8), dimension(ndim)     :: vmann
!
      real(8), dimension(ndim)     :: a=0.d0,an=0.d0,q=0.d0,qn=0.d0
!      real(8), dimension(ndim)     :: a0=0.d0,a1=0.d0,q0=0.d0,q1=0.d0
      real(8), dimension(ndim)     :: fa=0.d0,fq=0.d0
      real(8), dimension(nsd,ndim) :: xc ! centroide
      real(8), dimension(nsd,ndim+1) :: xn ! nodal
!
      real(8), dimension(ndim) :: bn ! base no instante tn
      real(8), dimension(ndim) :: zb ! cota do leito no instante n
      real(8), dimension(ndim) :: zdam ! forma do topo da barragem
!
!     secoes
!
      real(8), dimension(ndmax) :: xsec ! vetor com as coord. das secoes
      real(8), dimension(ndmax) :: bsec ! vetor com as larguras das secoes
!
!     cotas do fundo
!
      real(8), dimension(ndmax) :: xzb ! vetor com as coord. dos pontos de medicao das cotas
      real(8), dimension(ndmax) :: zbd ! vetor com as cotas medidas, zb data
!
!     dados inflow
!
      real(8), dimension(ndmax) :: thd ! vetor com as instantes de medicao altura inflow
      real(8), dimension(ndmax) :: tqd ! vetor com as instantes de medicao vazao inflow
      real(8), dimension(ndmax) :: hd ! vetor com as alturas medidas
      real(8), dimension(ndmax) :: qd ! vetor com as vazoes medidas
!
!     arrays pro test 5
!
      real(8), dimension(nobstime)  :: timeobs ! vetor com os instantes de observacao
      real(8), dimension(nstations) :: pointobs ! vetor com as posicoes de observacao
      real(8), dimension(nstations,nobstime) :: hstation=0.d0 ! vetor com as alturas medidas
      real(8), dimension(nstations,nobstime) :: vstation=0.d0 ! vetor com as velocidades medidas
!
      end module  
!
!----------------------------------------------------------------------
!
      module mprint
!
      implicit none
!
      integer :: npa,npq,nph,npv,npzb,npzh
!
      real(8) :: tprt_a, dtprt_a
      real(8) :: tprt_q, dtprt_q
      real(8) :: tprt_h, dtprt_h
      real(8) :: tprt_v, dtprt_v
      real(8) :: tprt_zb, dtprt_zb
      real(8) :: tprt_zh, dtprt_zh
!
      integer :: nprint_a,nprint_q,nprint_h,nprint_v,nprint_zb,nprint_zh
!      
      character(len=5) :: file_a   ='a'
      character(len=5) :: file_q   ='q'
      character(len=5) :: file_h   ='h'
      character(len=5) :: file_v   ='v'
      character(len=5) :: file_zb  ='zb'
      character(len=5) :: file_zh  ='zh'
!     
      integer :: iunit_a     =113
      integer :: iunit_q     =114
      integer :: iunit_h     =115
      integer :: iunit_v     =116
      integer :: iunit_zb    =117
      integer :: iunit_zh    =118
!
      integer :: iprint
!
      end module mprint
!
!======================================================================
!
!      program gedanken
!      implicit none
!
!      integer,parameter :: nx=50
!      integer :: ndiv
!      real(8), dimension(nx) :: hcal,mann
!      real(8) :: dx,dt,sigma,seed,tobs
!      real(8) :: xli,xlf
!
!      xli=0.d0
!      xlf=14.d0
!      mann=0.01d0
!      dx=(xlf-xli)/nx
!      dt=0.01d0
!      sigma=0.01d0
!      seed=9388847237.d0
!      tobs=200
!      hcal=0.d0
!
!      ndiv=nx-1
!
!      call dgsv(ndiv,mann,dx,dt,sigma,seed,tobs,hcal,xli)
!
!      end program
!
!======================================================================
!
      subroutine dgsv(ndiv,xmin,xmax,dtglobal,mann,sigma,dtsigma,seed,nobs, &
                      tobs,xobs,hcal,x,treali,hreali)
!
!      program upwindsv
!
!
      use mgeometry
      use marrays
      use mtime
!      use mprint, only : iprint,tprt_h,dtprt_h
      use mconstantes
!
      implicit none
!
!     SCALAR ARGUMENTS
!
      integer, intent(in) :: ndiv,nobs
      real(8), intent(in) :: dtglobal,sigma,dtsigma,treali
      real(8), intent(inout) :: seed
!
!     ARRAY ARGUMENTS
!
      real(8), intent(in) , dimension(ndiv+1) :: mann
      real(8), intent(out), dimension(ndiv+1) :: x,hreali
      real(8), intent(in) , dimension(nobs)   :: xobs,tobs
      real(8), intent(out), dimension(nobs)   :: hcal
!
      real(8), intent(in) :: xmin,xmax
!
      integer :: i
      real(8) :: dxnovo
      
!     discretization of x at which quantities are being calculated
      
      dxnovo = ( xmax - xmin ) / ( ndiv + 1 )
      x(1:ndiv+1) = (/ (xmin + dxnovo * ( i - 1 ) + dxnovo / 2.0d0, i=1,ndiv+1) /)
!
!     dimensoes iniciais (m)
!
!      h0=0.043d0 ! Fr=4.23
      h0=0.031d0 ! Fr=7.00
!
      b0=0.46d0!5.d0!3.0d0
!
!     velocity upstream
!
!      vup=2.737d0 ! Fr=4.23
      vup=3.831d0 ! Fr=7.0
!
!     hdown = downstream depth 
!  
!      hdown=0.222d0 ! Fr=4.23
      hdown=0.265d00 ! Fr=7.00
!
!     Toce Valley
!
!     dimensoes, comprimento (m)
!
!     dimensoes, comprimento (m)
!
      xli= xmin!0.3d0
      xlf= xmax!14.d0
!
!     numero de Courant prescrito
!
      cr=0.1d0!0.45d0!0.1d0!45d0!1d0!45d0
!      cr=0.45d0!0.45d0!0.1d0!45d0!1d0!45d0
!
!     tempo (s)
!
      ti=000.0d0
      tf=tobs(nobs)+1.d0!treali!tobs!4000.d0!3000.d0!1800.d0!4300!4300.d0!50.0d0!550.0d0
!
      tsigma=ti!+0.1d0
!      tsigma=tf+1.d0!+0.1d0
!      dtsigma=0.5d0!100*dtglobal!2.d0
!
!      xsigma =xli
!      dxsigma=0.3d0
!
!     divisoes
!
      nx=ndiv+1
!
      numel =nx
      nnodes=(nx+1)
!
      ndof=ns*numel
!
!     delta t global 
!
      dtg=dtglobal
!
!     novo, para o manning
!
      do i=1,numel
      vmann(i)=mann(i)
      end do ! i
!
!     inicializa a impressao para o gif
!
!      tprt_h=ti
!      dtprt_h=(tf-ti)/50
!
!     abre e le o arquivo com os dados
!
!      call readdata
!
!     inicializa as estacoes de observacao
!
!      call initobs
!
!     calcula as coordenadas
!
      call calcx
!
!     malha
!
      call calchxmin
!
!     inicializa a cota do fundo do canal
!
      call zbinit
!
!     inicializa a base do canal
!
      call binit
!
!     condicao inicial
!
      call uinit
!
!     tempo inicial
!
      tt=ti
!
!     imprime a solucao inicial para o gnuplot
!
!      iprint=0
!      if(iprint.eq.1)  call prtgnuplot(tt,xc,a,q,bn,zb,zdam)
!
!     transporte 
!
      call upwind(nobs,hcal,xobs,tobs,sigma,dtsigma,seed,treali,hreali) 
!
!     imprime a solucao final para o gnuplot
!
!      iprint=1
!      if(iprint.eq.1)  
!      call prtgnuplot(tt,xc,a,q,bn,zb,zdam)
!
!     imprime os histogramas nas estacoes
!
!      call prtstations
!
!      write(*,"(' ----------',/,' Success!',/,' ----------')")
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine upwind(nobs,hcal,xobs,tobs,sigma,dtsigma,seed,treali,hreali) 
!
      use mgeometry, only : hxmin,numel,ns
      use marrays
      use mtime
      use mprint, only : tprt_h, dtprt_h
! 
      implicit none
!
      integer, intent(in) :: nobs
      real(8), dimension(nobs) :: hcal 
      real(8), intent(out), dimension(numel):: hreali
      real(8), intent(in) , dimension(nobs) :: xobs,tobs
      real(8), intent(in) :: treali,sigma,dtsigma
      real(8) :: tsist,seed,x0,x1,x2,aux,h1,h2,tauxobs,xx
      integer :: nt,iobs,l,l1,l2
      real(8) :: perc=0.1d0
!
!     Inicializa variaveis de contagem do tempo
!
      nt=0
!
      tsist=(tf-ti)
!
!     para observacao nas estacoes
!
!      iobs=1
!      tobs=timeobs(iobs)
!
!     para observacoes
!
      iobs=1
      tauxobs=tobs(iobs)
!
!     comeca com um dt pequeno, depois ajusta. O maior possivel eh dtg
!
      dt=dtg
!      dt=1.d-6
!
!     loop no tempo
!
!     estimativa de limite superior para a velocidade maxima de propagacao
!
!     FORA, pega o maior e vai acumulando (so para teste mesmo)
!
!      alphamax=0.d0
!
      tempo: do while(tsist.gt.0.0)
!
!     controle do tempo restante
!
      if(dt.gt.tsist) then
      dt=tsist
      tsist=0.d0
      else
      tsist=tsist-dt
      end if
!     
!     estimativa de limite superior para a velocidade maxima de propagacao
!
!     DENTRO
!
      alphamax=0.d0
!
!     armazena a solucao do instante anterior
!
      an=a
      qn=q
!
!     evolui de tn para tn+1
!
!     utiliza o fluxo upwind
!
      call fluxupwind(a,q,an,qn,dt)
!
      if (sigma.gt.0.0d0) then
!
      if(tt.gt.tsigma) then
!
      do l=1,numel
!
      call rtermdg(a(l),sigma,seed)
      call rtermdg(q(l),sigma,seed)
!
      end do ! l
!
!      write(*,*) tt,tsigma
      tsigma=tsigma+dtsigma
!
      end if
!
      end if
!
!     atualiza do tempo 
!
      tt=tt+dt
      nt=nt+1
!
      if(abs(tt-treali).le.0.5d0*dt) then ! We know that tobs(1:nobs) = tobs(1)
!      
      hreali(1:numel)=a(1:numel)/bn(1:numel)
!
!      do l=1,numel
!      hreali(l)=a(l)/bn(l)
!      hcal(l)=a(l)/bn(l)
!      write(*,*) tt,l,hcal(l)
!
!      end do ! i
!      stop
!
!      return
      end if
!
!     checa a observacao
!
      if(tt.ge.tobs(iobs)) then
!
!      busca o xobs
!
      do l=1,numel
      xx=xc(1,l)
!
!     checa de novo no tauxobs, pois ele pode ser atualizado dentro d loop em l 
!
      if(iobs.le.nobs) then
!
      if(xx.ge.xobs(iobs).and.tt.ge.tobs(iobs)) then
!
      if(l.eq.1) then 
      write(*,*) 'observacao antes do primeiro elemento na upwind'
      stop
      end if
!
      l1=l-1
      l2=l
!
      x1=xc(1,l1)
      x2=xc(1,l2)
!
!     interpola no tempo
!
      h1=(an(l1)+(a(l1)-an(l1))/dt*(tobs(iobs)-(tt-dt)))/bn(l1)
      h2=(an(l2)+(a(l2)-an(l2))/dt*(tobs(iobs)-(tt-dt)))/bn(l2)
!
!      interpola para hcal
!
      hcal(iobs)=h1+(h2-h1)/(x2-x1)*(xobs(iobs)-x1)
!      write(*,*) iobs,tobs(iobs),xobs(iobs),hcal(iobs)
!
!     vai ser atualizado dentro do loop em l, pois o proximo tobs pode ser o mesmo
!
      iobs=iobs+1
!      tauxobs=tobs(iobs)
!
      end if 
      end if
!
      end do ! l
!
      end if
!
!      write(*,*) iobs,tobs(iobs),xobs(iobs),hcal(iobs)
!
!     recalcula o dt para o proximo passo, baseado no atual
!
      dt=cr*hxmin/(alphamax+1.d-6)
!      write(*,*) dt,dtg
      if(dt.gt.dtg) dt=dtg
!
!      dt=1.d-3
!
!     observacoes nas estacoes     
!
!      if(tt.ge.tobs) then
!
!      if(iobs.le.nobstime) then 
!
!      call obstations(iobs)
!
!      iobs=iobs+1
!      tobs=timeobs(iobs)
!      end if
!
!      end if
!
!     imprime h para gif animado
!
!      if(tt.ge.tprt_h) then 
!      call prth
!      tprt_h=tprt_h+dtprt_h
!      end if
!
!     atualiza a base e o bed slope
!
!      if((tt-ti)/(tf-ti).ge.perc) then
!      write(*,"(f6.1,'% done')") perc*100 
!      perc=perc+0.1d0
!      end if
!
      end do tempo ! do do while
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine fluxupwind(a,q,an,qn,dt)
!
!     calcula o vetor fluxo para o sistema de edos
!
      use mgeometry
      use marrays, only : xc,xn,bn,zb,idry,vmann!q0,q1
      use mtime, only : tt,alphamax
      use mconstantes, only : grav,hdry
!
      implicit none
!
      real(8) :: dt
      real(8), dimension(nx):: a,q,an,qn
!
      real(8), dimension(ns)   :: fl,fc,fr,ul,uc,ur
!
      real(8) :: xx,xm,xp,pi,hxc,hxl,hxr,hcr,hcl
      real(8) :: alpha,val
      real(8) :: s0,sf,manning,rh,ph,calcsf
      real(8) :: dzu,dzd,crd,cru,wd,wu,vd,vu
      real(8) :: al,ac,ar,ql,qc,qr
      real(8) :: zl,zc,zr
      real(8) :: qpico,tpico,qvale,tvale
!
      integer :: i,il,ir,ic
!
!      real(8) :: aghostin,qghostin,aghostout,qghostout
!
      pi=4.d0*datan(1.d0)
!
!     array que indica a existencia de leito seco. 0 se molhado, 1 se seco
!
      idry=0
!
!         ________
!        |        |
! (1) im |   i    | ip (2)
!        |________|
!
!     STEP 1 : area
!
      do i=1,nx     
!
      il=i-1
      ic=i
      ir=i+1
!
      if(i.eq. 1) il=i
      if(i.eq.nx) ir=i
!
      xx=xc(1,ic)
!
!     i eh o contador do bloco
!
!     matriz de massa é hh
!
      xm=xn(1,ic)
      xp=xn(1,ic+1)
!
      hxl=xn(1,il+1)-xn(1,il)
      hxc=xn(1,ic+1)-xn(1,ic)
      hxr=xn(1,ir+1)-xn(1,ir)
!
!     aproximacoes piecewise constant para calcular os fluxos, l c r
!
      al=an(il)
      ql=qn(il)
!
      ac=an(ic)
      qc=qn(ic)
!
      ar=an(ir)
      qr=qn(ir)
!
!     correcoes devido ao leito seco
!
      if(al.le.hdry*bn(il)) then
      al=hdry*bn(il)
      ql=0.d0
      end if
      if(ac.le.hdry*bn(ic)) then
      ac=hdry*bn(ic)
      qc=0.d0
      end if
      if(ar.le.hdry*bn(ir)) then
      ar=hdry*bn(ir)
      qr=0.d0
      end if
!
!     condicao de contorno esquerda porto, Ex. 14.9
!
!      qpico=200.d0!50.d0
!      qvale=8.249d0
!      tpico=20.d0*60.d0
!      tvale=30.d0*60.d0
!
!      if(i.eq.1) then
!
!      al=2.d0*ac-ar
!      
!      if(tt.le.tpico) then
!
!     sobre ate 20 min
!
!      ql=qvale+(qpico-qvale)/tpico*tt
!
!      else
!
!     cai por 30 min
!
!      if(tt.le.(tpico+tvale)) then
!
!      ql=qpico-(qpico-qvale)/(tvale)*(tt-tpico)
!
!      else 
!
!     constante
!
!      ql=qvale
!
!      end if
!      end if
!
!     imprime os valores de inflow (Porto)
!
!      write(202,"(3(e20.13))") tt,al/bn(1),ql
!
!      end if
!

!     condicao de contorno esquerda (Test 5)
!
      if(i.eq.1) then
!      
      val=h0*bn(1)
      al=val
!
      val=vup*bn(1)*h0
      ql=val
!
!      al=2.d0*ac-ar
!      ql=1.d0*qc-qr
!
      end if
!
!     condicao de contorno direita
!
      if(i.eq.nx) then
!
      if(tt.le.50.d0) then 
      val=(h0+tt*(hdown-h0)/50.d0)*bn(nx)
      else 
      val=hdown*bn(nx)
      end if
!
      ar=val
!
!      ar=2.d0*ac-al
      qr=2.d0*qc-ql
!
      end if
!
!     fluxos fisicos
!     
!      fl(1)=ul(2)
!      fl(2)=(ul(2)**2.d0)/ul(1)!+grav*(ul(1)**2.d0)/2.d0/bn(il)
!     
!      fc(1)=uc(2)
!      fc(2)=(uc(2)**2.d0)/uc(1)!+grav*(uc(1)**2.d0)/2.d0/bn(ic)
!     
!      fr(1)=ur(2)
!      fr(2)=(ur(2)**2.d0)/ur(1)!+grav*(ur(1)**2.d0)/2.d0/bn(ir)
!
!     esquerda
!
      if(ql.gt.0.d0.and.qc.gt.0.d0) then
      fl(1)=ql
      else 
      if(ql.lt.0.d0.and.qc.lt.0.d0) then
      fl(1)=qc
      else
      fl(1)=(ql+qc)/2.d0
      end if
      end if
!
!     direita
!
      if(qc.gt.0.d0.and.qr.gt.0.d0) then
      fr(1)=qc
      else 
      if(qc.lt.0.d0.and.qr.lt.0.d0) then
      fr(1)=qr
      else
      fr(1)=(qc+qr)/2.d0
      end if
      end if
!
!     evolucao: area
!
      a(ic)=an(ic)-dt/hxc*(fr(1)-fl(1))
!
!     correcao no contorno: 
!
!      if(i.eq. 1) a(i)=al!(ac+al)/2.d0
!
!      if(i.eq.nx) a(i)=(ac+ar)/2.d0
!
!     checa leito seco
!
      if(a(ic).le.hdry*bn(ic)) then
      a(ic)=hdry*bn(ic)
      idry(ic)=1
      end if
!
      end do ! i
!
!         ________
!        |        |
! (1) im |   i    | ip (2)
!        |________|
!
!     STEP 2 : vazao
!
      do i=1,nx     
!
      ir=i+1
      ic=i
      il=i-1
!
      if(i.eq. 1) il=ic
      if(i.eq.nx) ir=ic
!
!      if(idry(ic)*idry(ir)*idry(il).eq.1) then
      if(idry(ic).eq.1) then
      q(ic)=0.d0
      cycle
      end if
!
      xx=xc(1,ic)
!
!     i eh o contador do bloco
!
!     matriz de massa é hh
!
      xm=xn(1,ic)
      xp=xn(1,ic+1)
!
      hxl=xn(1,il+1)-xn(1,il)
      hxc=xn(1,ic+1)-xn(1,ic)
      hxr=xn(1,ir+1)-xn(1,ir)
!
!     aproximacoes piecewise constant para calcular os fluxos, l c r
!
      al=an(il)
      ac=an(ic)
      ar=an(ir)
!
      ql=qn(il)
      qc=qn(ic)
      qr=qn(ir)
!
!     correcoes devido ao leito seco
!
      if(al.le.hdry*bn(il)) then
      al=hdry*bn(il)
      ql=0.d0
      end if
      if(ac.le.hdry*bn(ic)) then
      ac=hdry*bn(ic)
      qc=0.d0
      end if
      if(ar.le.hdry*bn(ir)) then
      ar=hdry*bn(ir)
      qr=0.d0
      end if
!!
!     condicao de contorno esquerda porto, Ex. 14.9
!
!      if(i.eq.1) then
!
!      al=2.d0*ac-ar
!      
!      if(tt.le.tpico) then
!
!     sobre ate 20 min
!
!      ql=qvale+(qpico-qvale)/tpico*tt
!
!      else
!
!     cai por 30 min
!
!      if(tt.le.(tpico+tvale)) then
!
!      ql=qpico-(qpico-qvale)/(tvale)*(tt-tpico)
!
!      else 
!
!     constante
!
!      ql=qvale
!
!      end if
!      end if
!
!      end if
!
!
!     condicao de contorno esquerda
!
      if(i.eq.1) then
!         
!     supercritical: prescreve a e q na ghost cell E gradiente do nivel dagua igual a zero
!
      val=h0*bn(1)
      al=val
!
      val=vup*bn(1)*h0
      ql=val
!
!      al=2.d0*ac-ar
!      ql=1.d0*qc-qr
!
      end if
!
!     condicao de contorno direita
!
      if(i.eq.nx) then
!
      if(tt.le.50.d0) then 
      val=(h0+tt*(hdown-h0)/50.d0)*bn(nx)
      else 
      val=hdown*bn(nx)
      end if
!
      ar=val
!
!      ar=2.d0*ac-al
      qr=2.d0*qc-ql
!
      end if
!

!     fluxos fisicos
!     
!      fl(1)=ul(2)
!      fl(2)=(ul(2)**2.d0)/ul(1)!+grav*(ul(1)**2.d0)/2.d0/bn(il)
!     
!      fc(1)=uc(2)
!      fc(2)=(uc(2)**2.d0)/uc(1)!+grav*(uc(1)**2.d0)/2.d0/bn(ic)
!     
!      fr(1)=ur(2)
!      fr(2)=(ur(2)**2.d0)/ur(1)!+grav*(ur(1)**2.d0)/2.d0/bn(ir)
!
!     velocidades de propagacao
!
!
!     fluxos numericos, upwind
!
!     esquerda
!
      if(ql.gt.0.d0.and.qc.gt.0.d0) then
      fl(2)=(ql**2.d0)/al
      else 
      if(ql.lt.0.d0.and.qc.lt.0.d0) then
      fl(2)=(qc**2.d0)/ac
      else
      fl(2)=(((ql+qc)/2.d0)**2.d0)/((al+ac)/2.d0)
      end if
      end if
!
!     direita
!
      if(qc.gt.0.d0.and.qr.gt.0.d0) then
      fr(2)=(qc**2.d0)/ac
      else 
      if(qc.lt.0.d0.and.qr.lt.0.d0) then
      fr(2)=(qr**2.d0)/ar
      else
      fr(2)=(((qr+qc)/2.d0)**2.d0)/((ar+ac)/2.d0)
      end if
      end if
!
!     termo de fonte da derivada: usa area atualizada
!
!     Downwind e Upwind
!
      hcl=xc(1,ic)-xc(1,il)
      hcr=xc(1,ir)-xc(1,ic)
!
!     tem ue extrapolar Z e nao A
!
      zl=a(il)/bn(il)+zb(il)
      zc=a(ic)/bn(ic)+zb(ic)
      zr=a(ir)/bn(ir)+zb(ir)
!
      if(ic.eq.1 ) then
      hcl=hxc
!
!     se estiver livre
!
!      zl=2.d0*zc-zr
!
!     se estiver prescrito
!
      zl=h0+2.d0*zb(ic)-zb(ir)
!
      end if

      if(ic.eq.nx) then
      hcr=hxc
!
!     se estiver livre
!
!      zr=2.d0*zc-zl
!
!     se estiver prescrito
!
      if(tt.le.50.d0) then 
      val=(h0+tt*(hdown-h0)/50.d0)
      else 
      val=hdown
      end if
!
      zr=val+2.d0*zb(ic)-zb(il)      
!
      end if
!
!
      if(qc.ge.0.d0) then
      dzd=(zr-zc)/hcr 
      dzu=(zc-zl)/hcl
      crd=dt/hcr*(abs(qr)/ar+abs(qc)/ac)/2.d0
      cru=dt/hcl*(abs(qc)/ac+abs(ql)/al)/2.d0
      else      
      dzd=(zc-zl)/hcl 
      dzu=(zr-zc)/hcr 
      crd=dt/hcl*(abs(qc)/ac+abs(ql)/al)/2.d0
      cru=dt/hcr*(abs(qr)/ar+abs(qc)/ac)/2.d0
      end if
!
      if(ac.le.hdry*bn(ic)) then 
      alpha=sqrt(grav*hdry)
      else
      alpha=abs(qc)/ac+sqrt(grav*ac/bn(ic))
      end if 
!
      if(alpha.gt.alphamax) alphamax=alpha
!
!      write(*,*) i, alpha,alphamax
!
!     termo de fonte da derivada de zb
!
!     pesos: Eq (8) biased toward the downwind
!
      wd=1.d0-crd
      wu=cru
!
!     downwind
!
!      wd=1.0d0
!      wu=0.0d0
!
!     a escolha abaixo eh a central, causa solucoes nao-fisicas
!
!      wd=0.5d0
!      wu=0.5d0
!
!     pesos: Eq (9) biased toward the upwind
!
!      wd=1.d0-sqrt(crd)
!      wu=sqrt(cru)
!
      s0=-grav*a(ic)*(wd*dzd+wu*dzu)
!
!     Mario
!
!      s0=an(ic)*grav*s0/hx/(1.d0+(s0/hx)**2.d0)
!
!     termo de fonte: atrito
!      
      sf=calcsf(an(ic),qn(ic),bn(ic),xx,ic)
!
!     evolucao vazao
!
      q(ic)=qn(ic)-dt/hxc*(fr(2)-fl(2))+dt*s0+dt*sf
!
      end do ! i
!
!     computando o termo de atrito de manning externamente
!
!     primeiro estagio
!
!      do i=1,nx
!
!      sf=calcsf(a(i),q(i),bn(i),xc(1,i))
!
!      q0(i)=q(i)+dt*sf
!
!      end do
!
!     segundo estagio
!
!      do i=1,nx
!
!      sf=calcsf(a(i),q0(i),bn(i),xc(1,i))
!
!      q1(i)=3.d0/4.d0*q(i)+1.d0/4.d0*q0(i)+1.d0/4.d0*dt*sf
!
!      end do
!
!     terceiro estagio
!
!      do i=1,nx
!
!      sf=calcsf(a(i),q1(i),bn(i),xc(1,i))
!
!      q(i)=1.d0/3.d0*q(i)+2.d0/3.d0*q1(i)+2.d0/3.d0*dt*sf
!
!      end do
!
!
!     Armazena
!
      qn=q
      an=a
!
      end subroutine
!
!----------------------------------------------------------------------
!
      function calcsf(a,q,b,xx,l)
!
!     calcula o atrito de manning
!
      use mconstantes, only : grav, hdry
      use marrays, only : vmann
      implicit none
!
      integer :: l
      real(8) :: a,q,b,xx,calcsf
      real(8) :: rh,ph,manning
!      
      if(a.gt.hdry*b) then
!
!     perimetro molhado
!
      ph=b+2.d0*(a/b)
!
!     raio hidraulico
!
      rh=a/ph
!
!     forca de atrito
!
!      calcsf=-grav*(manning(xx)**2.d0)*q*abs(q)/(a*(rh**(4.d0/3.d0)))
      calcsf=-grav*(vmann(l)**2.d0)*q*abs(q)/(a*(rh**(4.d0/3.d0)))
!
!      write(*,*) l, vmann(l)
      else 
!
      calcsf=0.d0
!
      end if
!
      end function
!
!----------------------------------------------------------------------
!
      subroutine initobs
!
      use marrays, only : timeobs,pointobs
      use mdata, only: nstations,nobstime
      use mtime, only: ti,tf
      implicit none
!
      integer :: n
      real(8) :: to,dto,tinicial,tfinal
!
      tinicial=ti
      tfinal=200.d0!tf
!
      dto=(tfinal-tinicial)/nobstime
!
      to =tinicial+1.d-1
!
      do n=1,nobstime
      timeobs(n)=to
      to=to+dto
      end do ! n
!
      pointobs(1)= 30.500d0 ! (STA100)   
      pointobs(2)= 45.750d0 ! (STA150)   
      pointobs(3)= 68.625d0 ! (STA225)   
      pointobs(4)=106.750d0 ! (STA350)      
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine obstations(iobs)
!
      use marrays, only : hstation,vstation,timeobs,pointobs,xn,a,q,bn,zb
      use mdata, only: nstations,nobstime
      use mgeometry, only : nx
      implicit none
!
      integer :: iobs,l,npoint
      real(8) :: x1,x2,xp
!
      npoint=1
!
      do l=1,nx
!
      x1=xn(1,l  )
      x2=xn(1,l+1)
!
      xp=pointobs(npoint)      
!
      if(xp.ge.x1.and.xp.le.x2) then
!
      hstation(npoint,iobs)=a(l)/bn(l)
      vstation(npoint,iobs)=q(l)/a(l)
!
      if(npoint.lt.nstations) npoint=npoint+1
!
      end if
!
      end do ! l
!
!      write(*,*) iobs
!     
      end subroutine
!
!
!----------------------------------------------------------------------
!
      subroutine prtstations
!
      use marrays, only : hstation,vstation,timeobs,pointobs,xn,a,bn
      use mdata, only: nstations,nobstime
      use mgeometry, only : nx
      implicit none
!
      integer :: iobs,l,npoint,n
      real(8) :: x1,x2,xp
!
      npoint=1
!
      do n=1,nobstime
!
      write(280,"(100(e15.8,1x))") timeobs(n),(hstation(l,n),l=1,nstations)
      write(290,"(100(e15.8,1x))") timeobs(n),(vstation(l,n),l=1,nstations)
!
      end do ! n
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine readdata
!
      use marrays, only : xsec,bsec,xzb,zbd,qd,tqd,hd,thd
      use mdata, only: ndatasec,ndatazb,ndatah,ndataq,ireaddata
      implicit none
!
      integer :: i,res,n
!
      if(ireaddata.eq.0) then
!
!     abre e le o arquivo com os dados
!
      open(10, file="in-cross-sections.txt",access='sequential',form="formatted",iostat=res)
      open(11, file="in-floor.txt",access='sequential',form="formatted",iostat=res)
      open(12, file="in-inflowh.txt",access='sequential',form="formatted",iostat=res)
      open(13, file="in-inflowQ.txt",access='sequential',form="formatted",iostat=res)
!
!     Secoes
!
      res=0.d0
      i=0
      do while(res.eq.0)
      i=i+1
      read(10,*,iostat=res) xsec(i),bsec(i)
!      write(*,*) xsec(i),bsec(i)
      end do
!
!     numero de dados
!
      ndatasec=i-1
!
!     Cota do leito
!
      res=0.d0
      i=0
      do while(res.eq.0)
      i=i+1
      read(11,*,iostat=res) xzb(i),zbd(i)
!      write(*,*) xzb(i),zbd(i)
      end do
!
!     numero de dados
!
      ndatazb=i-1
!
!     altura h no bordo esquerdo
!
      res=0.d0
      n=0
      do while(res.eq.0)
      n=n+1
      read(12,*,iostat=res) thd(n),hd(n)
!      write(*,*) thd(n),hd(n)
      end do
!
!     numero de dados
!
      ndatah=n-1
!
!     altura h e vazao no bordo esquerdo
!
      res=0.d0
      n=0
      do while(res.eq.0)
      n=n+1
      read(13,*,iostat=res) tqd(n),qd(n)
!      write(*,*) tqd(n),qd(n)
      end do
!
!     numero de dados
!
      ndataq=n-1
!
      end if
!
!     muda para nao ler de novo
!
      ireaddata=1
!      qd=0.0
!      hd=0.0d0
!      bsec=3.d0
!      zbd=0.d0
!
!
!     suaviza bsec
!
!      do n=3,ndatasec-2
!      bsec(i)=(2*bsec(i-2)+4*bsec(i-1)+6*bsec(i)+2*bsec(i+1)+bsec(i+2))/15.d0
!      end do ! n
!
!      stop
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine uinit
      use marrays, only : a,an,q,qn,bn,xc,zb
      use mgeometry
      use mconstantes, only : hdry
!
      implicit none
!
      integer :: l
      real(8) :: xi,xf,xx
!
      an=0.0d0
      qn=0.d0
!
      xi=xli
      xf=xlf
!
      do l=1,numel 
      xx=xc(1,l)
!
!     BUMP
!
!      a(l)=(0.33d0-zb(l))*bn(l)
!      q(l)=0.d0 
!
!     Toce Valley
!
!      a(l)=0.d0!h0*bn(l)
!      q(l)=0.d0!vup*h0*bn(l)
!
!     Test 1
!
!      a(l)=10.d0*bn(l)
!      if(xx.gt.500.d0) a(l)=2.d0**bn(l)
!      q(l)=0.d0
!
!     Test 3: SALTO
!
      a(l)=h0*bn(l)
      q(l)=vup*h0*bn(l)
!
!     Test 5
!
!      a(l)=(0.61d0-zb(l))*bn(l)
!      if(xx.gt.(61.d0-0.61d0)) a(l)=0.d0!hdry*bn(l)
!      q(l)=0.d0
!
!
!      a(l)=1.20d0*bn(l)
!      q(l)=8.249d0
!
      end do ! l
!
!      do l=1,nx
!      write(*,*) a(l),q(l)
!      end do ! l
!

      an=a
      qn=q
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine calcx
      use marrays, only : xc,xn
      use mgeometry
!
      implicit none
!
      integer :: l,i,n
      real(8) :: hx,xx
!
!     centroides
!
      xc=0.d0
!
!     assumindo malha uniforme
!
      hx=(xlf-xli)/nx
!
      do i=1,nx     
      l=i
      xc(1,l)=xli+hx*i-hx/2.d0
      end do ! i
!
!     nodais
!
      xn=0.d0
!
      n=0
      do i=1,nx+1
      n=n+1
      xn(1,n)=xli+hx*(i-1)
!      write(*,*) n,xn(1,n)!,xn(2,l)
      end do ! i
!
!     le os dados nodais para 64 nos
!
!      open(50, file="in-64-nos.txt",status='unknown')
!
!      do i=1,nx+1
!      read(50,*) xx
!      xn(1,i)=xx
!      write(*,*) xn(1,i)
!      end do ! i
!
!      do i=1,nx
!      xc(1,i)=(xn(1,i)+xn(1,i+1))/2.d0
!      write(*,*) xc(1,i)
!      write(*,*) xn(1,i)
!      end do ! i
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine calchxmin
      use marrays, only : xn
      use mgeometry, only : numel,hxmin
!
!     inicializa a cota da base do canal
!
      implicit none
      integer :: i
      real(8) :: xx,hx
!
      hxmin=1.d16
!
      do i=1,numel
      hx=xn(1,i+1)-xn(1,i)
      if(hx.lt.hxmin) hxmin=hx
      end do ! i   
!
!      write(*,*) hxmin; stop
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine zbinit
      use marrays, only : zb,xc,xzb,zbd
      use mgeometry
      use mdata, only : ndatazb
!
!     inicializa a cota da base do canal
!
      implicit none
      integer :: l,i
      real(8) :: xx,x1,x2,z1,z2
!
!     cota do fundo 
!
!      zb  =zbd(ndatazb)
!
!      do l=1,numel
!      xx=xc(1,l)
!
!      x1=xli
!      z1=zbd(1)
!      do i=1,ndatazb
!      x2=xzb(i)
!      z2=zbd(i)
!      if(xx.ge.x1.and.xx.lt.x2) then
!
!     interpolacao linear
!
!      zb(l)=z1+(z2-z1)/(x2-x1)*(xx-x1)
!
!      exit
!      end if
!      x1=x2
!      z1=z2
!      end do ! i
!
!      write(*,*) l,zb(l)
!
!      end do ! l
!
      zb=0.d0
!
!     Test 5
!
!      do l=1,nx
!      zb(l)=0.005d0*(xlf-xli-xc(1,l))
!      end do ! l
!
!      do l=1,nx
!
!      xx=xc(1,l)
!
!     Condition 11.1
!
!      if(xx.ge.(60.96d0-0.6096d0).and.xx.le.(60.96d0+0.6096d0)) zb(l)=zb(l)+(0.3048d0-0.18288d0)!(0.305d0-0.183d0)
!
!     Condition 10.1
!
!      if(xx.ge.(61.d0-0.61d0).and.xx.le.(61.d0+0.61d0)) zb(l)=zb(l)+(0.305d0-0.0305d0)!(0.305d0-0.183d0)
!
!      end do ! l
!
!     Porto
!
!      do l=1,nx
!      zb(l)=0.001d0*(xlf-xli-xc(1,l))
!      end do ! l
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine binit
      use marrays, only : bn,xc,xsec,bsec
      use mgeometry
      use mdata, only : ndatasec
!
!     inicializa a cota da base do canal
!
      implicit none
      integer :: l,i
      real(8) :: xx,x1,x2,b1,b2
!
!     largura
!
      bn=b0
!
!      bn=bsec(ndatasec)
!
!      do l=1,numel
!      xx=xc(1,l)
!
!      x1=xli
!      b1=bsec(1)
!      do i=1,ndatasec
!      x2=xsec(i)
!      b2=bsec(i)
!      if(xx.ge.x1.and.xx.lt.x2) then
!
!     interpolacao linear
!
!      bn(l)=b1+(b2-b1)/(x2-x1)*(xx-x1)
!
!      exit
!      end if
!      x1=x2
!      b1=b2
!      end do ! i
!
!      write(*,*) l,bn(l)
!
!      bn(l)=bn(l)/2.d0
!
!      end do ! l
!
!      stop
!
!     63 secoes
!
!      open(60, file="in-63-secoes.txt",status='unknown')
!
!      do i=1,nx
!      read(60,*) xx,b1
!      bn(i)=b1/4.d0
!      write(*,*) bn(i)
!      end do ! i
!
!     Test 5
!
!      bn=1.2192d0
!
!      do l=1,nx
!
!      xx=xc(1,l)
!
!     Condition 11.1
!
!      if(xx.ge.(61.d0-0.61d0).and.xx.le.(61.d0+0.61d0))  bn(l)=0.7315d0
!
!     Condition 10.1
!
!      if(xx.ge.(61.d0-0.61d0).and.xx.le.(61.d0+0.61d0))  bn(l)=bn(l)
!
!      end do ! l
!
!     Porto 
!
!     bn=5.00d0
!
      end subroutine
!
!
!----------------------------------------------------------------------
!
      subroutine findmin(am,ap,qm,qp,bm,bp,vmin)
!
!     busca valores minimos
!
      implicit none
      integer :: i,j,k,ndiva,ndivq,ndivb
      real(8) :: am,ap,qm,qp,bm,bp,da,dq,db,a,q,b,lambda1,v,vmin
!
      ndiva=1
      ndivq=1
      ndivb=1
!     
      da=(ap-am)/ndiva
      dq=(qp-qm)/ndivq
      db=(bp-bm)/ndivb
!
      vmin=1.d20
!
      a=am
      do i=1,ndiva+1
      q=qm
      do j=1,ndivq+1
      b=bm
      do k=1,ndivb+1
      v=lambda1(a,q,b)
      if(v.lt.vmin) vmin=v
      b=b+db
      end do ! k
      q=q+dq
      end do ! j
      a=a+da
      end do ! i
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine findmax(am,ap,qm,qp,bm,bp,vmax)
!
!     busca valores maximos
!
      implicit none
      integer :: i,j,k,ndiva,ndivq,ndivb
      real(8) :: am,ap,qm,qp,bm,bp,da,dq,db,a,q,b,lambda2,v,vmax
!
      ndiva=1!0
      ndivq=1
      ndivb=1
!     
      da=(ap-am)/ndiva
      dq=(qp-qm)/ndivq
      db=(bp-bm)/ndivb
!
!
      vmax=-1.d20
!
      a=am
      do i=1,ndiva+1
      q=qm
      do j=1,ndivq+1
      b=bm
      do k=1,ndivb+1
      v=lambda2(a,q,b)
      if(v.gt.vmax) vmax=v
      b=b+db
      end do ! k
      q=q+dq
      end do ! j
      a=a+da
      end do ! i
!
      end subroutine
!
!----------------------------------------------------------------------
!
      function lambda1(a,q,b)
      use mconstantes
      implicit none
      real(8) :: lambda1,a,q,b
!
      lambda1=q/a-sqrt(grav*a/b)
!
      end function
!
!----------------------------------------------------------------------
!
      function lambda2(a,q,b)
      use mconstantes
      implicit none
      real(8) :: lambda2,a,q,b
!
      lambda2=q/a+sqrt(grav*a/b)
!
      end function
!
!----------------------------------------------------------------------
!
      function manning(xx)
      implicit none
      real(8) :: manning,xx
!
!      manning=0.014d0!d-2!0.d0*xx
!      manning=0.02d0!d-2!0.d0*xx
!      manning=0.0162d0
!      manning=0.009d0
      manning=0.01d0
!
      end function
!
!
!----------------------------------------------------------------------
!
      function minmod(a,b)
!
!
!     If a and b have the same sign, then minmod will return whichever
!     has smaller modulus. Otherwise, it will return 0.
!
      implicit none
!
      real(8) :: a,b,minmod,z
!
      z=0.d0
!
      if (a*b.ge.z) then
      if (dabs(a).lt.dabs(b)) then
      minmod = a
      else
      minmod = b
      end if
      else
      minmod = z
      end if
!
      return
      end
!
!----------------------------------------------------------------------
!
      subroutine prtgnuplot(t0,xc,a,q,bn,zb,zdam)
!
      use mgeometry
      use mconstantes
!
      implicit none
!
!     imprime campos escalares para o gnuplot ou para o matlab
!
      real(8), dimension(nsd,numel) :: xc
      real(8), dimension(numel)     :: a,q,bn,zb,zdam
      real(8)                       :: t0
!
      integer :: l
      real(8) :: xx,yy
!
      integer :: iunit=270
!
      write(iunit,"('#TIMESTEP PRINT OUT = ',f15.8)") t0
!
      write(iunit,*)
!
      do l=1,numel
!
      xx=xc(1,l)
!      yy=xc(2,l)
!
      write(iunit,"(25(f25.15,2x))") xx,a(l),q(l),a(l)/bn(l),q(l)/a(l),bn(l)/2.d0,zb(l),xx,xx,zdam(l),-bn(l)/2.d0, &
                                     zb(l)+a(l)/bn(l),abs(q(l))/a(l)/sqrt(grav*a(l)/bn(l)),xx,zdam(l)
!
      end do ! l
!
      write(iunit,*)
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine prth
!
      use mgeometry, only : nx
      use marrays, only : xc,a,bn,zb
      use mtime, only : tt
!
      implicit none
!
!     imprime campos escalares para o gnuplot ou para o matlab
!
      integer :: l
      real(8) :: xx,t0
!
      integer :: iunit=300
!
      write(iunit,"('#TIMESTEP PRINT OUT = ',f15.8)") tt
!
      write(iunit,*)
!
      do l=1,nx
!
      xx=xc(1,l)
!
      write(iunit,"(15(f25.15,2x))") xx,zb(l),a(l)/bn(l)+zb(l)
!
      end do ! i
!
      write(iunit,*)
!
      end subroutine
!
!
! ******************************************************************

      subroutine rtermdg(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 
