!
!.----------------------------------------------------------------------
!
      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=801 ! 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=1000
      integer :: ndatasec,ndatazb,ndatah,ndataq,ndatazinit
      integer :: ireaddata=0
      integer :: iprtstations=0
      integer, parameter :: nstations=3
      integer, parameter :: nobstime =100
!
      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-6
!
      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)     :: 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 ! largura (base) nodal
      real(8), dimension(ndim) :: bc ! largura (base) centrada
      real(8), dimension(ndim) :: zbn! cota do leito nodal
      real(8), dimension(ndim) :: zbc! cota do leito centrada
      real(8), dimension(ndim) :: zdam ! forma do topo da barragem
!
!     vetor com o Manning
!
      real(8), dimension(ndim) :: mc ! vetor com os valores de n, centrados
!
!     secoes LIDAS
!
      real(8), dimension(ndmax) :: xsec ! vetor com as coord. das secoes
      real(8), dimension(ndmax) :: bsec ! vetor com as larguras das secoes
!
!     ziniciais LIDAS
!
      real(8), dimension(ndmax) :: xzinit ! vetor com as coord. das ziniciais
      real(8), dimension(ndmax) :: zinit ! vetor com as ziniciais
!
!     cotas do fundo LIDAS
!
      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 LIDOS
!
      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) :: zd ! vetor com os z medidos
      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) :: zstation=0.d0 ! vetor com as z medidas
      real(8), dimension(nstations,nobstime) :: qstation=0.d0 ! vetor com as vazoes 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
!
!======================================================================
!
!      program svupwind
!
!
!======================================================================
!
      subroutine dgsv(nshydro,nthydro,ndiv,xmin,xmax,dtglobal,mann,sigma,dtsigma,seed,nobs, &
                      tobs,xobs,hcal,x,treali,tmax,hreali,zhydro,xhydro,thydro)
      use mgeometry
      use marrays
      use mtime
!      use mprint, only : iprint,tprt_h,dtprt_h
!      use mconstantes
!
      implicit none
!
!     SCALAR ARGUMENTS
!
      integer, intent(in) :: nshydro,nthydro,ndiv,nobs
      real(8), intent(in) :: dtglobal,sigma,dtsigma,treali,tmax
      real(4), intent(inout) :: seed
!
!     ARRAY ARGUMENTS
!
      real(8), intent(in) , dimension(ndiv+1) :: mann
      real(8), intent(inout), dimension(ndiv+1) :: x,hreali
      real(8), intent(in) , dimension(nobs)   :: xobs,tobs
      real(8), intent(out), dimension(nobs)   :: hcal
      real(8), intent(out), dimension(nthydro,nshydro):: zhydro
      real(8), intent(out), dimension(nthydro):: thydro
      real(8), intent(out), dimension(nshydro):: xhydro
!
      real(8), intent(inout) :: xmin,xmax
!
!     Fork: calculado internamente em calcx
!
!     apenas para fins graficos. tb ta sendo usado no binit
!
!      yli=-0.5d0
!      ylf= 0.5d0
!
!     numero de Courant prescrito
!
      cr=0.1d0
!      cr=0.01d0
!
!     tempo (s)
!
!     dia 1: 20/05
!
!      write(*,*) tmax
!
      ti= 0.d0*24.d0*60.d0*60.d0
!      tf=27.d0*24.d0*60.d0*60.d0 ! dias * h/dia * min/h * s/min
      tf=1.1d0*tmax!17.d0*24.d0*60.d0*60.d0 ! dias * h/dia * min/h * s/min
!      tf=30.d0*24.d0*60.d0*60.d0 ! dias * h/dia * min/h * s/min
!      tf=3.01d0*24.d0*60.d0*60.d0 ! dias * h/dia * min/h * s/min
!      tf=0.5d0*24.d0*60.d0*60.d0 ! dias * h/dia * min/h * s/min
!
!      write(*,*) tf,tmax;stop
!
!     inicializa o tsigma
!
      tsigma=ti!+0.1d0
!
!     divisoes
!
!     Fork river
!
!      nx=60
!
!      numel =nx
!      nnodes=(nx+1)
!
!      ndof=ns*numel
!
!     inicializa a impressao para o gif
!
!      tprt_h=ti
!      dtprt_h=(tf-ti)/50
!
!     abre e le o arquivo com os dados
!
!     divisoes
!
!     Fork river
!
      nx=ndiv+1
!
      numel =nx
      nnodes=(nx+1)
!
      ndof=ns*numel
!
!     delta t global 
!
      dtg=dtglobal
!
!     abre e le o arquivo com os dados
!
      if(ireaddata.eq.0) then
!
      call readdata
!
!     malha nao uniforme
!
!     Fork river
!
!      nx=ndatasec
!
      if(nx.ne.ndatasec) then
      write(*,*) 'nx diferente do ndatasec',nx,ndatasec
      stop
      end if
!
      numel =nx
!
      nnodes=(nx+1)
!
      ndof=ns*numel
!
!     inicializa as estacoes de observacao
!
!      call initobs
!
!     calcula as coordenadas
!
      call calcx(x)
!
      xmin=xli
      xmax=xlf
!
!     malha
!
      call calchxmin
!
!     inicializa a cota do fundo do canal
!
      call zbinit
!
!     inicializa a base do canal
!
      call binit
!
      end if ! ireaddata
!
!
!     condicao inicial
!
      call uinit
!
!     inicializa o vetor com os mannings
!
      call manninit(mann)
!
!     tempo inicial
!
      tt=ti
!
!     imprime a solucao inicial para o gnuplot
!
!      iprint=0
!      if(iprint.eq.1)  call prtgnuplot(tt,xc,a,q,bc,zbc,zdam)
!
!     transporte 
!
      call upwind(nshydro,nthydro,nobs,hcal,xobs,tobs,sigma,dtsigma,seed,treali,hreali,&
                  zhydro,thydro,xhydro) 
!
!     imprime a solucao final para o gnuplot
!
!      iprint=1
!      if(iprint.eq.1)  call prtgnuplot(tt,xc,a,q,bc,zbc,zdam)
!
!     imprime os histogramas nas estacoes
!
!      call prtstations
!
!      write(*,"(' ----------',/,' Success!',/,' ----------')")
!
!      end program
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine upwind(nshydro,nthydro,nobs,hcal,xobs,tobs,sigma,dtsigma,seed,treali,hreali,&
                        zhydro,thydro,xhydro)  
!
      use mgeometry, only : hxmin,numel,ns
      use marrays
      use mtime
      use mconstantes, only : hdry
      use mprint, only : tprt_h, dtprt_h
! 
      implicit none
!
      integer, intent(in) :: nshydro,nthydro,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,x0,x1,x2,aux,z1,z2,tauxobs,xx,hh
      real(4), intent(inout) :: seed
      real(8), intent(out), dimension(nthydro,nshydro):: zhydro
      real(8), intent(out), dimension(nthydro):: thydro
      real(8), intent(out), dimension(nshydro):: xhydro
!
      real(8) :: perc,tauxhydro
      integer :: nt,iobs,l,l1,l2,ihydro,inshydro,nst
!
      perc=0.1d0
!
!     Inicializa variaveis de contagem do tempo
!
      nt=0
!
      tsist=(tf-ti)
!
!     para observacao nas estacoes
!
      iobs=1
      tauxobs=timeobs(iobs)
!
!     hidrogramas
!
      ihydro=1
      tauxhydro=thydro(ihydro)
!
!     comeca com um dt pequeno, depois ajusta
!
      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 fluxo upwind
!
      call fluxupwind(a,q,an,qn,dt)
!
!     stuff after evolving in time
!
      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)
!
      if(a(l)/bc(l).le.hdry) then
      a(l)=hdry*bc(l)
      q(l)=0
      end if
!
      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)/bc(1:numel)+zbc(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
!!!
!!      z1=(an(l1)+(a(l1)-an(l1))/dt*(tobs(iobs)-(tt-dt)))/bc(l1)+zbc(l1)
!!      z2=(an(l2)+(a(l2)-an(l2))/dt*(tobs(iobs)-(tt-dt)))/bc(l2)+zbc(l2)
!!!
!!!      interpola para hcal
!!!
!!      hcal(iobs)=z1+(z2-z1)/(x2-x1)*(xobs(iobs)-x1)
!!!
!!!     hidrograma
!!!
!!!      write(220,*) tt,hcal(iobs)
!!!
!!!      write(*,*) h1,h2
!!!      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
!
!     checa a observacao SIMPLIFICADO
!
      if(tt-dt.le.tobs(iobs).and.tt.ge.tobs(iobs)) then
!      busca o xobs
!
      do l=1,numel
      xx=xc(1,l)
      x1=xn(1,l)
      x2=xn(1,l+1)
!
      if(x2.ge.xobs(iobs).and.x1.le.xobs(iobs)) then
!
      hcal(iobs)=a(l)/bc(l)+zbc(l)
!
      iobs=iobs+1
!      tauxobs=tobs(iobs)
!
      end if 
!
      end do ! l
!
      end if
!
!     HIDROGRAMAS
!
      if(tt-dt.le.thydro(ihydro).and.tt.ge.thydro(ihydro)) then
!
!      busca o xobs
!
      do nst=1,nshydro
!
      do l=1,numel
      xx=xc(1,l)
      x1=xn(1,l)
      x2=xn(1,l+1)
!
      if(x2.ge.xhydro(nst).and.x1.le.xhydro(nst)) then
!
      zhydro(ihydro,nst)=a(l)/bc(l)+zbc(l)
!
      exit
!
      end if
!
      end do ! l
!
      end do ! nst
!
      ihydro=ihydro+1
!
      end if
!
!     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
!
!     recalcula o dt para o proximo passo, baseado no atual
!
      dt=cr*hxmin/((1.d0+sigma)*alphamax+1.d-6)
!      write(*,*) dt,sigma
!      write(*,*) dt,dtg
!
      if(dt.gt.dtg) then
      write(*,*) 'dt maior que dtg', dt,dtg
      dt=dtg
      end if
!
      end do tempo ! do do while
!
!      do nst=1,nthydro
!      write(*,*) thydro(nst),zhydro(nst)
!      end do 
!      stop
!
      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,bc,zbc,idry!,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) :: fal,far,fql,fqr
!
      real(8) :: xx,xm,xp,pi,hxc,hp,hm
      real(8) :: alpha
      real(8) :: s0,sf,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) :: aal,aac,aar,zzbl,zzbc,zzbr,bbl,bbc,bbr,ext_l,ext_r
!
      integer :: i,il,ir,ic
!
      pi=4.d0*datan(1.d0)
!
!     array que indica a existencia de leito seco. 0 se molhado, 1 se seco
!
      idry=0
!
!         ________
!        |        |
! (1) il |   ic   | ir (2)
!        |________|
!
!     STEP 1 : area
!
      do i=2,nx-1     
!
      il=i-1
      ic=i
      ir=i+1
!
      xx=xc(1,ic)
!
!     i eh o contador do bloco
!
      xm=xn(1,ic  )
      xp=xn(1,ic+1)
!
      hxc=xp-xm
!
      hm=xc(1,ic)-xc(1,il)
      hp=xc(1,ir)-xc(1,ic)
!
!     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)
!
!     larguras 
!
!      bbl=bc(il)
!      bbc=bc(ic)
!      bbr=bc(ir)
!
!     no contorno: obtidas por extrapolacao
!     
!      if(ic.eq. 1) bbl=ext_l(bbc,bbr,hcl,hcr)
!      if(ic.eq.nx) bbr=ext_r(bbc,bbl,hcl,hcr)
!
!     correcoes devido ao leito seco
!
!      if(al.le.hdry*bbl) then
!      al=hdry*bbl
!      ql=0.d0
!      end if
!      if(ac.le.hdry*bbc) then
!      ac=hdry*bbc
!      qc=0.d0
!      end if
!      if(ar.le.hdry*bbr) then
!      ar=hdry*bbr
!      qr=0.d0
!      end if
!
!     condicao de contorno esquerda 
!
      if(ic.eq.2) then
!
      call bcl(al,ac,ar,ql,qc,qr,ic)
      a(1)=al
      q(1)=ql
!      write(*,*) a(1),al; stop
!
      end if
!
!     condicao de contorno direita
!
      if(ic.eq.nx-1) then
!
      call bcr(al,ac,ar,ql,qc,qr,ic)
      a(nx)=ar
      q(nx)=qr
!
      end if
!
!     fluxos fisicos
!      
!      fl(1)=ul(2)
!      fl(2)=(ul(2)**2.d0)/ul(1)!+grav*(ul(1)**2.d0)/2.d0/bc(il)
!     
!      fc(1)=uc(2)
!      fc(2)=(uc(2)**2.d0)/uc(1)!+grav*(uc(1)**2.d0)/2.d0/bc(ic)
!     
!      fr(1)=ur(2)
!      fr(2)=(ur(2)**2.d0)/ur(1)!+grav*(ur(1)**2.d0)/2.d0/bc(ir)
!
!     esquerda
!
      if(ql.gt.0.d0.and.qc.gt.0.d0) then
      fal=ql
      else 
      if(ql.lt.0.d0.and.qc.lt.0.d0) then
      fal=qc
      else
      fal=(ql+qc)/2.d0
      end if
      end if
!
!     direita
!
      if(qc.gt.0.d0.and.qr.gt.0.d0) then
      far=qc
      else 
      if(qc.lt.0.d0.and.qr.lt.0.d0) then
      far=qr
      else
      far=(qc+qr)/2.d0
      end if
      end if
!
!     evolucao: area
!
      a(ic)=an(ic)-dt/hxc*(far-fal)
!
!     checa leito seco
!
      if(a(ic).le.hdry*bbc) then
      a(ic)=hdry*bbc
      idry(ic)=1
      end if
!
      end do ! i
!         ________
!        |        |
! (1) il |   ic   | ir (2)
!        |________|
!
!     STEP 2 : vazao
!
      do i=2,nx-1     
!
      ir=i+1
      ic=i
      il=i-1
!
!      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
!
      xm=xn(1,ic)
      xp=xn(1,ic+1)
!
      hxc=xp-xm
!
      hm=xc(1,ic)-xc(1,il)
      hp=xc(1,ir)-xc(1,ic)
!
!     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)
!
!     larguras 
!
      bbl=bc(il)
      bbc=bc(ic)
      bbr=bc(ir)
!
!     no contorno: obtidas por extrapolacao, nao sao usadas
!     
!      if(ic.eq. 1) bbl=ext_l(bbc,bbr,hcl,hcr)
!      if(ic.eq.nx) bbr=ext_r(bbc,bbl,hcl,hcr)
!
!     correcoes devido ao leito seco
!
!      if(al.le.hdry*bbl) then
!      al=hdry*bbl
!      ql=0.d0
!      end if
!      if(ac.le.hdry*bbc) then
!      ac=hdry*bbc
!      qc=0.d0
!      end if
!      if(ar.le.hdry*bbr) then
!      ar=hdry*bbr
!      qr=0.d0
!      end if
!
!     condicao de contorno esquerda
!
      if(i.eq.2) then
!
!      call bcl(a(il),a(ic),a(ir),ql,qc,qr,ic)
      call bcl(al,ac,ar,ql,qc,qr,ic)
      a(1)=al
      q(1)=ql
!
      end if
!
!     condicao de contorno direita
!
      if(i.eq.nx-1) then
!
!      call bcr(a(il),a(ic),a(ir),ql,qc,qr,ic)
      call bcr(al,ac,ar,ql,qc,qr,ic)
      a(nx)=ar
      q(nx)=qr
!
      end if
!
!     fluxos fisicos
!     
!      fl(1)=ul(2)
!      fl(2)=(ul(2)**2.d0)/ul(1)!+grav*(ul(1)**2.d0)/2.d0/bc(il)
!     
!      fc(1)=uc(2)
!      fc(2)=(uc(2)**2.d0)/uc(1)!+grav*(uc(1)**2.d0)/2.d0/bc(ic)
!     
!      fr(1)=ur(2)
!      fr(2)=(ur(2)**2.d0)/ur(1)!+grav*(ur(1)**2.d0)/2.d0/bc(ir)
!
!
!     fluxos numericos, upwind
!
!     esquerda
!
      if(ql.gt.0.d0.and.qc.gt.0.d0) then
      fql=(ql**2.d0)/al
      else 
      if(ql.lt.0.d0.and.qc.lt.0.d0) then
      fql=(qc**2.d0)/ac
      else
      fql=(((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
      fqr=(qc**2.d0)/ac
      else 
      if(qc.lt.0.d0.and.qr.lt.0.d0) then
      fqr=(qr**2.d0)/ar
      else
      fqr=(((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
!
!     cotas do fundo
!
      zzbl=zbc(il)
      zzbc=zbc(ic)
      zzbr=zbc(ir)
!
!     no contorno: obtidas por extrapolacao, nao serao usadas
!     
!      if(ic.eq. 1) zzbl=ext_l(zzbc,zbr,hcl,hcr)
!      if(ic.eq.nx) zzbr=ext_r(zzbc,zbl,hcl,hcr)
!
!     areas em n+1
!
      aal=a(il)
      aac=a(ic)
      aar=a(ir)
!
!     corrige nos bordos
!
!      if(ic.eq. 1) aal=ext_l(aac,aar,hcl,hcr)
!      if(ic.eq.nx) aar=ext_r(aac,aal,hcl,hcr)
!
!     tem que extrapolar Z e nao A
!
      zl=aal/bbl+zzbl
      zc=aac/bbc+zzbc
      zr=aar/bbr+zzbr
!
!     o trecho a seguir eh dependente da condicao de contorno adotada
!
!      if(ic.eq.2) then
!
!     descomente aqui se quiser z livre
!
!      xx=xc(1,1)
!      zl=ext_l(xx,xc(1,ic),xc(1,ir),zc,zr)
!
!      end if
!
!      if(ic.eq.nx-1) then
!
!     descomente aqui se quiser z livre
!
!      xx=xc(1,nx)
!      zr=ext_r(xx,xc(1,ic),xc(1,il),zc,zl)
!     
!      end if
!
      if(qc.gt.0.d0) then
      dzd=(zr-zc)/hp 
      dzu=(zc-zl)/hm
      crd=dt/hp*(abs(qr)/ar+abs(qc)/ac)/2.d0
      cru=dt/hm*(abs(qc)/ac+abs(ql)/al)/2.d0
      else    
      if(qc.lt.0.d0) then
      dzd=(zc-zl)/hm 
      dzu=(zr-zc)/hp 
      crd=dt/hm*(abs(qc)/ac+abs(ql)/al)/2.d0
      cru=dt/hp*(abs(qr)/ar+abs(qc)/ac)/2.d0
      else
!
      dzd=(zr-zl)/(hp+hm)
      dzu=dzd
!
      crd=dt/hxc*(abs(ql)/al+2.d0*abs(qc)/ac+abs(qr)/ar)/2.d0
      cru=crd
      end if      
      end if
!
      if(aac.le.hdry*bbc) then 
      alpha=sqrt(grav*hdry)
      else
      alpha=abs(qc)/ac+sqrt(grav*ac/bbc)
      end if 
!      write(*,*) ic,alpha
!
      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
!
!      wd=(1.d0-crd)/(1.d0-crd+cru)
!      wu=cru/(1.d0-crd+cru)
!
!     downwind
!
!      wd=1.0d0
!      wu=0.0d0
!
!     upwind
!
!      wd=0.0d0
!      wu=1.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),bc(ic),xx,ic)
!
!     evolucao vazao
!
!      q(ic)=qn(ic)-dt/hxc*(fr(2)-fl(2))+dt*s0+dt*grav*a(ic)*sf
      q(ic)=qn(ic)-dt/hxc*(fqr-fql)+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),bc(i),xc(1,i))
!
!      q0(i)=q(i)+dt*sf
!
!      end do
!
!     segundo estagio
!
!      do i=1,nx
!
!      sf=calcsf(a(i),q0(i),bc(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),bc(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
!
!----------------------------------------------------------------------
!
      subroutine bcl(al,ac,ar,ql,qc,qr,ic)
!
!     condicao de contorno na esquerda
!
      use marrays, only : thd,tqd,hd,qd,bc,zbc,xc,zd,xn
      use mdata, only : ndatah,ndataq
      use mconstantes, only : hdry
      use mgeometry, only : h0,xli
      use mtime, only : tt
!
      implicit none
!
      real(8) :: al,ac,ar,ql,qc,qr
      real(8) :: val
      real(8) :: qpico,tpico,qvale,tvale
      real(8) :: t1,t2,h1,h2,hh,q1,q2,z1,z2
      real(8) :: ain,qin,hin,vin,fr,hout,zbl,bl
      real(8) :: zl,zc,zr,xx
      real(8) :: ext_l
!
      integer :: n,ic
!
!     obs, supercritical: prescreve a e q na ghost cell E gradiente do nivel dagua igual a zero
!
!     livre
!
!      al=2.d0*ac-ar
!      ql=2.d0*qc-qr
!
!     area prescrita
!      
!      val=h0*bc
!      al=val
!
!     vazao prescrita
!
!      val=0.01d0*bc*h0
!      ql=val
!
!     PORTO
!
!     condicao de contorno esquerda porto, Ex. 14.9
!
!      qpico=200.d0!50.d0
!      qvale=8.245d0
!      tpico=20.d0*60.d0
!      tvale=40.d0*60.d0
!
!      al=2.d0*ac-ar
!      
!      if(tt.le.tpico) then
!
!     sobe 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
!
!     TOCE RIVER
!
!     dados de altura inflow
!
!      t1=thd(1)
!      h1=hd(1)
!      do n=1,ndatah-1
!      t2=thd(n+1)
!      h2=hd(n+1)
!      if(tt.gt.t1.and.tt.le.t2) then
!      hin=h1+(h2-h1)/(t2-t1)*(tt-t1)
!      write(*,*) t1,tt,t2
!      write(*,*) tt,hh
!      exit
!      end if
!      t1=t2
!      h1=h2
!      end do ! n
!
!     dados de vazao inflow
!
!      t1=tqd(1)
!      q1=qd(1)
!      do n=1,ndataq-1
!      t2=tqd(n+1)
!      q2=qd(n+1)
!      if(tt.ge.t1.and.tt.lt.t2) then
!      qin=q1+(q2-q1)/(t2-t1)*(tt-t1)
!      write(*,*) tt,qin
!      exit
!      end if
!      t1=t2
!      q1=q2
!      end do ! n
!
!      write(201,"(2(e20.13))") tt,hin,qin
!     
!     checa se eh subcritico ou supercritico
!
!      if(hin.lt.hdry) hin=hdry
!      ain=hin*bn(1)
!
!     vazao na entrada
!
!      vin=qin/ain
!
!     Numero de Froude
!
!      fr=vin/sqrt(grav*hin)
!
!     modificando para ter escoamento critico, Fr=1
!
!      fr=1.d0
!      hin=(1.d0/(grav*fr**2.d0)*(qin/bn(1))**2.d0)**(1.d0/3.d0)
!      if(hin.lt.hdry) hin=hdry
!      ain=hin*bn(1)
!
!      vin=qin/ain
!
!      fr=vin/sqrt(grav*hin)
!
!      write(*,*) tt,fr,qin/bn(1)/sqrt(grav)/(hin)**(3.d0/2.d0)
!      write(*,*) qin,ain,fr,hin
!
!      write(*,*) vin,hin,fr,qin
!
!      fr=2.d0
!
!      if(fr.le.1.d0) then ! subcritico
!
!     extrapola a altura
!
!      hc=ac/bn(1)
!      hp=ap/bn(2)
!
!      hm=hc-(hp-hc)/(xp-xc)*hx
!      if(hm.lt.hdry) hm=hdry
!
!      am=hm*bn(1)
!
!     extrapola a area
!
!      am=ac-(ap-ac)/(xp-xc)*hx
!      if(am/bn(1).lt.hdry) am=hdry*bn(1)
!
!     prescreve a vazao
!
!      qm=0.d0
!      write(*,*) tt,am,qm
!
!      else ! supercritico
!
!     prescreve a altura      
!     
!      al=hin*bc
!
!      write(*,*) tt,hin
!
!     zero gradient
!
!      ac=am
!     
!     prescreve a vazao
!
!      ql=qin
!
!      end if
!
!     imprime os dados
!
!      write(202,"(3(e20.13))") tt,hin,qin
!
!
!     correcao leito seco
!
!      if(al.lt.hdry*bc) then
!      al=hdry*bc
!      ql=0.d0
!      end if
!
!     fork river: water surface no beadload trap
!
      t1=thd(1)
      h1=hd(1)
      z1=zd(1)
!      write(*,*) t1,h1,ndatah
      do n=1,ndatah-1
      t2=thd(n+1)
      h2=hd(n+1)
      z2=zd(n+1)
      if(tt.ge.t1.and.tt.le.t2) then
      hout=h1+(h2-h1)/(t2-t1)*(tt-t1)
      zl=z1+(z2-z1)/(t2-t1)*(tt-t1)
!      write(*,*) t1,tt,t2,hout
      exit
      end if
      t1=t2
      h1=h2
      end do ! n
!
!      al=bc(ic)*hout
!
      bl=bc(ic-1)
      al=bl*hout
!      write(*,*) bc(ic-1),bl,hout,al
!
!     fork river livre
!
!      al=2.d0*ac-ar
!      ql=2.d0*qc-qr
!
!
!      al=2.d0*ac-ar
!      ql=2.d0*qc-qr
!
!      bl(bc(ic),bc(ic+1),xc(1,ic),xc(1,ic+1),hxc)
!
!      zc=ac/bc(ic)+zbc(ic)
!      zr=ar/bc(ic+1)+zbc(ic+1)
!
!     extrapola
!
!      zl=2.d0*zc-zr
!      zl=zc+hxc/(xc(1,ic+1)-xc(1,ic))*(zc-zr)
!
!      zl=ext_l(zc,zr,hcl,hcr)
!
!     calcula a area molhada
!
!      zbl=zbc(ic)+hxc/(xc(1,ic+1)-xc(1,ic))*(zbc(ic)-zbc(ic+1))
!
!      zbl=ext_l(zbc(ic),zbc(ic+1),hcl,hcr)
!
!      bl=ext_l(bc(ic),bc(ic+1),hcl,hcr)
!      bl =bc(ic) +hxc/(xc(1,ic+1)-xc(1,ic))*(bc(ic) -bc(ic+1))
!
!!!      al=(zl-zbl)*bl
!
!      al=(zl-(2.d0*zbc(ic)-zbc(ic+1)))*(2.d0*bc(ic)-bc(ic+1))
!      write(*,*) al,zl,zc,zr
!
!      al=2.d0*ac-ar
!
!      al=ac+hxc/(xc(1,ic+1)-xc(1,ic))*(ac-ar)
!
!      write(*,*) al
!      al=ext_l(ac,ar,xc(1,ic),xc(1,ic+1),hxc)
!      write(*,*) al
!      write(*,*) 
!
!      al=-(xc(1,ic)+xc(1,ic+1))/(xc(1,ic)-xc(1,ic+1))*ac+2.d0*xc(1,ic)/(xc(1,ic)-xc(1,ic+1))*ar
!      ql=-(xc(1,ic)+xc(1,ic+1))/(xc(1,ic)-xc(1,ic+1))*qc+2.d0*xc(1,ic)/(xc(1,ic)-xc(1,ic+1))*qr
!
!      ql=qc+hxc/(xc(1,ic+1)-xc(1,ic))*(qc-qr)
      xx=xc(1,ic-1)
!
      zc=zbc(ic  )+ac/bc(ic  )
      zr=zbc(ic+1)+ar/bc(ic+1)
      zl=ext_l(xx,xc(1,ic),xc(1,ic+1),zc,zr)
!
!      write(*,*) ic,zl,zc,zr;stop
!
      zbl=zbc(ic-1)!ext_l(xx,xc(1,ic),xc(1,ic+1),zbc(ic),zbc(ic+1))
      bl =bc (ic-1)!ext_l(xx,xc(1,ic),xc(1,ic+1),bc(ic),bc(ic+1))
!!
!
!     ATENCAO
!
!!     descomente aqui se quiser area livre e comente se quiser z prescrito
!!
!!      al=(zl-zbl)*bl
!
!      write(*,*) ic,al,ac,ar;stop
!
!      al=ext_l(xx,xc(1,ic),xc(1,ic+1),ac,ar)
!
      ql=ext_l(xx,xc(1,ic),xc(1,ic+1),qc,qr)
!
!      ql=qc
!      al=ac
!
!
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine bcr(al,ac,ar,ql,qc,qr,ic)
!
!     condicao de contorno na direita
!
      use mconstantes, only : hdry
!
      use marrays, only : thd,tqd,hd,qd,bc,zbc,xc,xn
      use mdata, only : ndatah,ndataq
      use mtime, only : tt,ti
      use mgeometry, only : xlf
!
      implicit none
!
      real(8) :: al,ac,ar,ql,qc,qr
      real(8) :: t1,t2,h1,h2,hh,q1,q2
      real(8) :: ain,qin,hin,vin,fr
      real(8) :: zl,zc,zr,zbr,br,ext_r,xx
      integer :: n,ic
!
!     livre
!
!      ar=2.d0*ac-al
!      qr=2.d0*qc-ql
!
!     area prescrita
!
!      if(tt.le.50.d0) then 
!      val=(h0+tt*(hdown-h0)/50.d0)*bc(nx)
!      else 
!      val=hdown*bc(nx)
!      end if
!
!      ar=val
!
!     vazao prescita
!
!     qr=val
!
!     fork river
!
!     dados de vazao inflow
!
      t1=ti!tqd(1)
      q1=qd(1)
      do n=1,ndataq-1
      t2=tqd(n+1)
      q2=qd(n+1)
      if(tt.ge.t1.and.tt.lt.t2) then
      qin=q1+(q2-q1)/(t2-t1)*(tt-t1)
!      write(*,*) t1,tt,t2,qin
      exit
      end if
      t1=t2
      q1=q2
      end do ! n
!
      qr=-qin
!
!     livre
!
!      zc=ac/bc(ic)+zbc(ic)
!      zl=al/bc(ic-1)+zbc(ic-1)
!
!     calcula zr por extrapolacao
!
!      zr=2.d0*zc-zl
!
!     extrapola
!
!      zl=2.d0*zc-zr
!      zr=zc+hxc/(xc(1,ic)-xc(1,ic-1))*(zc-zl)
!
!     calcula a area molhada
!
!      zbr=zbc(ic)+hxc/(xc(1,ic)-xc(1,ic-1))*(zbc(ic)-zbc(ic-1))
!      br =bc(ic) +hxc/(xc(1,ic)-xc(1,ic-1))*(bc(ic) -bc(ic-1))
!
!      al=(zl-zbc(ic))*(bc(ic))
!      ar=(zr-zbr)*br
!     calcula a area molhada
!
!
! old
      xx=xc(1,ic+1)
      ar=ext_r(xx,xc(1,ic),xc(1,ic-1),ac,al)
!
! new
!
      zc=zbc(ic  )+ac/bc(ic  )
      zl=zbc(ic-1)+al/bc(ic-1)
      zr=ext_r(xx,xc(1,ic),xc(1,ic-1),zc,zl)
!
!      write(*,*) ic,zl,zc,zr;stop
!
      zbr=zbc(ic+1)
      br =bc (ic+1)
!
!     descomente aqui se quiser area livre e comente se quiser z prescrito
!
      ar=(zr-zbr)*br
!
!      write(*,*) ic,al,ac,ar;stop
!
!      al=ext_l(xx,xc(1,ic),xc(1,ic+1),ac,ar)
!
!      ar=(zr-(2.d0*zbc(ic)-zbc(ic-1)))*(2.d0*bc(ic)-bc(ic-1))
!
!      ar=ac+hxc/(xc(1,ic)-xc(1,ic-1))*(ac-al)
!
!     correcao leito seco
!
!      if(ar.le.hdry*bc) then
!      ar=hdry*bc
!      qr=0.d0
!      end if
!
      end subroutine
!
!----------------------------------------------------------------------
!
      function ext_l(xx,xc,xr,uc,ur)
!
!     extrapola aa esquerda
!
      implicit none
!
      real(8) :: ext_l
      real(8), intent(in) :: xx,xc,xr,uc,ur
      real(8) :: hh
!
      hh=xr-xc
!
      ext_l=((xx-xc)*ur-(xx-xr)*uc)/hh
!
      end function
!
!----------------------------------------------------------------------
!
      function ext_r(xx,xc,xl,uc,ul)
!
!     extrapola aa direita
!
      implicit none
!
      real(8) :: ext_r
      real(8), intent(in) :: xx,xc,xl,uc,ul
      real(8) :: hh
!
      hh=xl-xc
!
      ext_r=((xx-xc)*ul-(xx-xl)*uc)/hh
!      
      end function
!
!----------------------------------------------------------------------
!
      function calcsf(a,q,b,xx,l)
!
!     calcula o atrito de manning
!
      use marrays, only : mc
      use mconstantes, only : grav, hdry
      implicit none
!
      integer :: l
      real(8) :: a,q,b,xx,calcsf
      real(8) :: rh,ph
!      
      if(a.gt.hdry*b) then
!
!     perimetro molhado
!
      ph=b+2.d0*(a/b)
!
!     raio hidraulico
!
      rh=a/ph
!
!      write(*,*) xx,a,ph,rh
!
!     forca de atrito
!
!      calcsf=-(mc(l)**2.d0)*q*abs(q)/((a**2.d0)*(rh**(4.d0/3.d0)))
!      calcsf=-(mc(l)**2.d0)*q*abs(q)/(a**2.d0)/(rh**(4.d0/3.d0))
!
!     manning classico
!
!      calcsf=-grav*(mc(l)**2.d0)*q*abs(q)/(a*(rh**(4.d0/3.d0)))
!
!     maning bombardelli
!
      calcsf=-(mc(l)**2.d0)*q*abs(q)/(a*(rh**(4.d0/3.d0)))
!
      else 
!
      calcsf=0.d0
!
      end if
!
      end function
!
!----------------------------------------------------------------------
!
      subroutine initobs
!
      use marrays, only : timeobs,pointobs
      use mdata, only: nstations,nobstime,iprtstations
      use mtime, only: ti,tf
      implicit none
!
      integer :: n,res
      real(8) :: to,dto,tinicial,tfinal
!
      if(iprtstations.eq.0) then
!
!     cria os arquivos de saida
!
      open(821, file="out-stations-h.dat",access='sequential',form="formatted",iostat=res)
      open(822, file="out-stations-z.dat",access='sequential',form="formatted",iostat=res)
      open(823, file="out-stations-q.dat",access='sequential',form="formatted",iostat=res)
      open(824, file="out-stations-v.dat",access='sequential',form="formatted",iostat=res)
!
      end if
!
      iprtstations=1
!
!     nao to usando mais o ti e o tf
!     
      tinicial=(16.d0+0.5d0)*24.d0*60.d0*60.d0
      tinicial=ti!3.d0*24.d0*60.d0*60.d0
      tinicial=ti+2*3600!3.d0*24.d0*60.d0*60.d0
      tfinal=tf!27.d0*24.d0*60.d0*60.d0
!
      dto=(tfinal-tinicial)/nobstime
!
      to =tinicial+1.d-1
!
      do n=1,nobstime
      timeobs(n)=to
      to=to+dto
      end do ! n
!
!     dam break
!
!      pointobs(1)= 30.500d0 ! (STA100)
!      pointobs(2)= 45.750d0 ! (STA150)
!      pointobs(3)= 68.625d0 ! (STA225)
!      pointobs(4)=106.750d0 ! (STA350)
!
!     TOCE RIVER
!
!      pointobs(1)= 1.d0
!      pointobs(2)=10.d0
!      pointobs(3)=20.d0
!      pointobs(4)=40.d0
!
!     fork river 
!
      pointobs(1)= 0.d0
      pointobs(2)= 2505.d0
      pointobs(3)= 3290.d0
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine obstations(iobs)
!
      use marrays, only : zstation,hstation,vstation,qstation,timeobs,pointobs,xn,a,q,bc,zbc
      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)/bc(l)
      zstation(npoint,iobs)=a(l)/bc(l)+zbc(l)
      qstation(npoint,iobs)=q(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 : zstation,hstation,vstation,qstation,timeobs,pointobs,xn,a,bc
      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(821,"(100(e15.8,1x))") timeobs(n),(hstation(l,n),l=1,nstations)
      write(822,"(100(e15.8,1x))") timeobs(n),(zstation(l,n),l=1,nstations)
      write(823,"(100(e15.8,1x))") timeobs(n),(-qstation(l,n),l=1,nstations)
      write(824,"(100(e15.8,1x))") timeobs(n),(-vstation(l,n),l=1,nstations)
!
      end do ! n
!
      write(821,*)
      write(822,*)
      write(823,*)
      write(824,*)
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine readdata
!
      use marrays, only : xsec,bsec,xzb,zbd,qd,tqd,hd,thd,zinit,xzinit,zd
      use mdata, only: ndatasec,ndatazb,ndatah,ndataq,ndatazinit,ireaddata
      implicit none
!
      integer :: i,res,n
!
      if(ireaddata.eq.0) then
!
!================================================================================================
!
!     Secoes
!
!     fork river, do relatorio 2, usando larguras no dia 20/05 - tabela 42
!
      open(810, file="in-01-fork-report2-crosssec-5-20.dat",access='sequential',form="formatted",iostat=res)
!      open(810, file="in-fork-report2-bequiv-sections.dat",access='sequential',form="formatted",iostat=res)
!
      res=0.d0
      i=0
      do while(res.eq.0)
      i=i+1
      read(810,*,iostat=res) xsec(i),bsec(i)
!      write(*,*) i,xsec(i),bsec(i)
      end do
!      stop
!
      close(810)
!
!     numero de dados
!
      ndatasec=i-1
!
!================================================================================================
!
!     batimetria
!
!     batimetria do fork river usando as cotas medias no dia 20 de maio, tabela 41 do report 2
!
      open(811, file="in-02-fork-report2-mean-bed-elevations-5-20.dat",access='sequential',form="formatted",iostat=res)
!
!     01 de junho
!
!      open(811, file="in-ok-fork-report2-t41-mean-bed-elevations-6-01.dat",access='sequential',form="formatted",iostat=res)
!
!      open(811, file="in-fork-zbmedio.dat",access='sequential',form="formatted",iostat=res)
!
!     12 de junho, Lai-Khan
!
!      open(811, file="in-ok-fork-zb-lai-khan.dat",access='sequential',form="formatted",iostat=res)
!
!     Cota do leito
!
!     inflow, secao 3256
!     zb=7.72
!
!     outflow, secao 0043
!     zb=4.94
!
      res=0.d0
      i=0
      do while(res.eq.0)
      i=i+1
      read(811,*,iostat=res) xzb(i),zbd(i)
!      write(*,*) i,xzb(i),zbd(i)
      end do
!      stop
!
      close(811)
!
!     numero de dados
!
      ndatazb=i-1
!
!================================================================================================
!
!     Dados de contorno
!
!
!     inflow h nao ta sendo usado
!
!      open(12, file="in-inflowh.txt",access='sequential',form="formatted",iostat=res)
!
!     outflow: water surface, fork
!
!      open(12, file="in-fork-bedload-water-height.dat",access='sequential',form="formatted",iostat=res)
!
!     outflow: water surface, fork
!
!      open(12, file="in-ok-fork-table-1-water-surf-section-0000.dat",access='sequential',form="formatted",iostat=res)
!
!     retirado diretamente a tabela 1 do report 1, medidos a 1AM e 1 PM do dia 20/05 ao dia 19/06
!
      open(812, file="in-03-from-20-05-fork-r1-table-1-water-surf-section-0000.dat",access='sequential',form="formatted",iostat=res)
!
!     altura h no bordo esquerdo
!
!     fork river: outflow: water surface no beadload trap
!
      res=0.d0
      n=0
      do while(res.eq.0)
      n=n+1
      read(812,*,iostat=res) thd(n),zd(n)
!
!     modifica para segundos
!
!      thd(n)=thd(n)*24.d0*60.d0*60.d0
!
      thd(n)=thd(n)*60.d0*60.d0
!
!     modifica para h
!
      hd(n)=zd(n)-zbd(1)
!
!      write(*,*) thd(n),hd(n),zd(n)
!
      end do
!      stop
!
      close(812)
!
!     numero de dados
!
      ndatah=n-1
!
!     inflow discharge fork river
!     # M. Tamer Ayvaz, http://dx.doi.org/10.1016/j.jhydrol.2013.07.019
!     # Figure 4(a), Discharge, Section 3295, Iniciando em 17/05/1979
!     t (days) discharge (m^3/s)
!
!      open(13, file="in-fork-3295-discharge-inflow-ayvaz.dat",access='sequential',form="formatted",iostat=res)
!
!     retirado diretamente a tabela 7 do report 1, medidos a 1AM e 1 PM do dia 17/05 ao dia 02/07
!
!      open(13, file="in-ok-fork-table-7-discharge-section-3295.dat",access='sequential',form="formatted",iostat=res)
!
!     retirado diretamente a tabela 7 do report 1, medidos a 1AM e 1 PM do dia 20/05 ao dia 19/06
!
      open(813, file="in-04-from-20-05-fork-r1-table-7-discharge-section-3295.dat",access='sequential',form="formatted",iostat=res)
!
!     agora ta de hora em hora
!
!      open(813, file="in-from-20-05-fork-r1-table-7-hourly-discharge-section-3295.dat", &
!                   access='sequential',form="formatted",iostat=res)
!
!     vazao no bordo esquerdo
!
      res=0.d0
      n=0
      do while(res.eq.0)
      n=n+1
      read(813,*,iostat=res) tqd(n),qd(n)     
!      write(*,*) tqd(n),qd(n)      
!      write(202,"(f8.2,2x,f8.2)") tqd(n)-72.d0,qd(n)
!
!     modifica para segundos
!
!      tqd(n)=tqd(n)*24.d0*60.d0*60.d0
!
      tqd(n)=tqd(n)*60.d0*60.d0
!
!      write(*,*) n,tqd(n),qd(n)
      end do
!      qd=0.d0
!      stop
!
      close(813)
!
!     numero de dados
!
      ndataq=n-1
!
!================================================================================================
!
!     Dados de cond iniciais
!     ziniciais
!
!      open(14, file="in-fork-zinicial.dat",access='sequential',form="formatted",iostat=res)
!
!     dia 21/05, retirado da table 10. report 1. completei com interpolacao linear
!
      open(814, file="in-05-fork-zinicial-21-05-r1-t10.dat",access='sequential',form="formatted",iostat=res)
!
      res=0.d0
      i=0
      do while(res.eq.0)
      i=i+1
      read(814,*,iostat=res) xzinit(i),zinit(i)
!      write(*,*) i,xzinit(i),zinit(i)
      end do
!      stop
!
!     numero de dados
!
      ndatazinit=i-1
!
      end if
!
      close(814)
!
!     muda para nao ler de novo
!
      ireaddata=1
!
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine uinit
      use marrays, only : a,an,q,qn,bc,xc,zbc,xzinit,zinit
      use mgeometry
      use mconstantes, only : hdry
      use mdata, only : ndatazinit
!
      implicit none
!
      integer :: l,i
      real(8) :: xi,xf,xx,hinit,zzinit,zli,zlf,x1,x2,z1,z2
!
      an=0.0d0
      qn=0.d0
!
      a =0.d0
      q =0.d0
!
      xi=xli
      xf=xlf
!
      hinit=2.50d0
!
!      do l=1,numel
!      xx=xc(1,l)
!
!      simetrico
!
!      a(l)=(hinit-zbc(l))*bc(l)
!
!      if(xx.lt.-0.4d0) a(l)=0.d0
!      if(xx.gt. 0.4d0) a(l)=0.d0
!
!      if(xx.lt.0.d0) then
!      q(l)=-0.d-1
!      else
!      if(xx.gt.0.d0) then
!      q(l)= 0.d-1
!      else
!      q(l)=0.d0
!      end if
!      end if
!
!     PORTO
!
!      a(l)=1.20d0*bc(l)
!      q(l)=8.245d0
!
!      end do ! l
!
!     Toce Valley
!
!      a(l)=0.d0!h0*bn(l)
!      q(l)=0.d0!vup*h0*bn(l)
!
!     Fork river
!
!      hinit=0.35d0
!
      zli=6.0d0
      zlf=8.0d0
      vup=0.2d-0
      do l=1,numel
      xx=xc(1,l)
!
      x1=xli-1.d-6
      z1=zinit(1)
      do i=1,ndatazinit
      x2=xzinit(i)
      z2=zinit(i)
      if(xx.ge.x1.and.xx.le.x2) then
!
!     interpolacao linear
!
      zzinit=z1+(z2-z1)/(x2-x1)*(xx-x1)
!
      exit
      end if
      x1=x2
      z1=z2
      end do ! i
!
!      zzinit=zbc(l)+hdry
!      if(xx.ge.3000.0) zzinit=zbc(numel)+5.d0
!      zzinit=(xx-xlf)/(xli-xlf)*zli+(xx-xli)/(xlf-xli)*zlf
!
!!!      hinit=(zzinit-zbc(l))
!
!      hinit=(zinit(l)-zbc(l))
!
      hinit=(zinit(l)-zbc(l))
!
!     media pq agora eh nodal
!
!      hinit=((zinit(l+1)+zinit(l))/2.d0-zbc(l))
!
!      hinit=0.5d0+zbc(l)
      a(l)=hinit*bc(l)
!      a(l)=hdry*bc(l)
!      q(l)=-vup*a(l)
!      q(l)=-3.9d0
      q(l)=-8.76d0
!      q(l)=0.d0!
!      write(*,*) l,a(l),q(l)
      end do ! l
!
!      stop
!
      an=a
      qn=q
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine calcx(x)
      use marrays, only : xc,xn,xsec
      use mgeometry, only : xli,xlf,nx
      use mdata, only: ndatasec

      implicit none
!
      real(8),dimension(nx), intent(inout) :: x
      integer :: l,i,n
      real(8) :: hx,xx
!
!     dimensoes:  usando o xsec do fork river
!
!      xli=xsec(1)
!      xlf=xsec(ndatasec)
!      write(*,*) ndatasec; stop
!
!     nodais
!
      xn=0.d0
!
!     centroides
!
      xc=0.d0      
!
!     assumindo malha uniforme
!
!      hx=(xlf-xli)/nx
!
!     nodais
!
!      do i=1,nx+1
!      xn(1,i)=xli+hx*(i-1)
!      write(*,*) n,xn(1,n)!,xn(2,l)
!      end do ! i
!
!     centroides=secoes
!
      do i=1,nx
      xc(1,i)=xsec(i)
!      write(*,*) n,xn(1,n)!,xn(2,l)
      end do ! i
!
!     nodais, internos
!
      do i=2,nx
      xn(1,i)=(xc(1,i)+xc(1,i-1))/2.d0
      end do ! i
!
!     extremos
!
      xn(1,   1)=2.d0*xc(1, 1)-xn(1, 2)
      xn(1,nx+1)=2.d0*xc(1,nx)-xn(1,nx)
!
      xli=xn(1,   1)
      xlf=xn(1,nx+1)
!
!      do i=1,nx
!      write(*,*) i,xn(1,i),xc(1,i),xn(1,i+1)
!      end do ! i
!
!      stop
!
!      write(*,*) xli,xlf;stop
!
      do i=1,nx
      x(i)=xc(1,i)
!      write(*,*) n,xn(1,n)!,xn(2,l)
      end do ! 

      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 : zbc,xc,xn,zdam,xzb,zbd,zbn
      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
!
      zbc =0.0d0
      zdam=0.0d0
!
!     simetrico
!
!      do l=1,numel
!      xx=xc(1,l)
!
!      zbc(l)=(0.5d0-abs(xx))*1.0d0
!      if(zbc(l).lt.0.d0) zbc(l)=0.d0
!      zdam(l)=zbc(l)
!
!      end do ! l
!
!     Porto
!
!      do l=1,nx
!      zbc(l)=0.001d0*(xlf-xli-xc(1,l))
!      end do ! l
!
!     cota do fundo TOCE RIVER e fork river
!
      zbc  =zbd(ndatazb)
!      zbn  =zbd(ndatazb)
!
!      do l=1,numel+1
!      xx=xn(1,l)
!
!      x1=xli-1.d-6
!      z1=zbd(1)
!      do i=1,ndatazb
!      x2=xzb(i)
!      z2=zbd(i)
!      if(xx.ge.x1.and.xx.le.x2) then
!
!     interpolacao linear
!
!      zbn(l)=z1+(z2-z1)/(x2-x1)*(xx-x1)
!      zbn(l)=z2
!
!     piecewise constant
!
!      zbc(l)=(z2+z1)/2.d0
!!
!      exit
!      end if
!      x1=x2
!      z1=z2
!      end do ! i
!
!      write(*,*) l,xx,zbn(l)
!
!      end do ! l
!
!     constante por secao
!
!     nodal
!
!      do l=1,nx+1
!      zbn(l)=zbd(l)
!      end do ! l
!      stop
!
!     centroides
!
      do l=1,numel
!      zbc(l)=(zbn(l+1)+zbn(l))/2.d0
      zbc(l)=zbd(l)
!      zbc(l)=zbn(l)
!      write(*,*) xc(1,l),zbn(l),zbc(l),zbn(l+1)
      end do ! l
!
!     nodais
!
!     intrnos
!
      do l=2,nx
      zbn(l)=(zbc(l)+zbc(l-1))/2.d0
      end do ! l
!
!     extremos
!
      zbn(1   )=zbc(1 )
      zbn(nx+1)=zbc(nx)
!      stop
!
!      do l=1,nx
!      write(*,*) zbn(l),zbc(l),zbn(l+1)
!      end do ! l
!      stop
!
      end subroutine

!
!----------------------------------------------------------------------
!
      subroutine binit
      use marrays, only : bc,xc,xsec,bsec,bn,xn
      use mgeometry
      use mdata, only : ndatasec
!
!     inicializa a largura (base) do canal
!
      implicit none
      integer :: l,i
      real(8) :: xx,x1,x2,b1,b2
!
!     largura
!
!      do l=1,numel
!      xx=xc(1,l)
!
!      if(xx.lt.0.d0) then
!      bc(l)=ylf-yli
!      else
!      bc(l)=b0
!      end if
!
!      end do ! l
!
!      bc=(ylf-yli)
!
!     Porto 
!
!      bc=5.00d0
!
!     TOCE RIVER e fork river
!
!     largura
!
      bc=bsec(ndatasec)
      bn=bsec(ndatasec)
!
!      do l=1,numel+1
!      xx=xn(1,l)
!
!      x1=xli-1.d-6
!      b1=bsec(1)
!      do i=1,ndatasec
!      x2=xsec(i)
!      b2=bsec(i)
!      if(xx.ge.x1.and.xx.le.x2) then
!
!     interpolacao linear
!
!      bn(l)=b1+(b2-b1)/(x2-x1)*(xx-x1)
!      bn(l)=b2
!
!     piecewise constant
!
!      bc(l)=(b1+b2)/2.d0
!
!      exit
!      end if
!      x1=x2
!      b1=b2
!      end do ! i
!
!      write(*,*) l,xx,bc(l)
!
!      bc(l)=sqrt(bc(l))
!
!      end do ! l
!
!     constate por secao
!
!     nodal
!
!      do l=1,nx+1
!      bn(l)=bsec(l)
!      end do ! l
!
!     centrada
!
      do l=1,numel
!      bc(l)=(bn(l+1)+bn(l))/2.d0
!      bc(l)=bn(l+1)
      bc(l)=bsec(l)
!      write(*,*) l,xc(1,l),bc(l)
      end do ! l
!
!     nodais
!
!     intrnos
!
      do l=2,nx
      bn(l)=(bc(l)+bc(l-1))/2.d0
      end do ! l
!
!     extremos
!
      bn(1   )=bc(1 )
      bn(nx+1)=bc(nx)
!      stop
!
!      do l=1,nx
!      write(*,*) bn(l),bc(l),bn(l+1)
!      end do ! l
!      stop
!      stop
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine manninit(mann)
      use marrays, only : mc,xc
      use mgeometry
!
!     inicializa o vetor com os mannings
!
      implicit none
      real(8), dimension(numel), intent(in) :: mann
      integer :: l
      real(8) :: manning,xx
!
      mc=0.d0
!
      do l=1,numel
!      xx=xc(1,l)
!
!      mc(l)=manning(xx)
      mc(l)=mann(l)
!      write(*,*) l,mc(l),mann(l)
!
      end do ! l
!
      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.007d0
!      manning=0.01d0
!      manning=0.50d0
!
!     Porto
!
!      manning=0.0248d0
!
!     TOCE RIVER
!
!     Bombardelli
!
      manning=0.131d0
!      manning=0.13d0
      manning=0.028
      manning=0.025
      manning=0.07d0
!      manning=0.035
!
      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,bc,zbc,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,bc,zbc,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)/bc(l),q(l)/a(l),bc(l)/2.d0,zbc(l),xx,xx,zdam(l),-bc(l)/2.d0, &
                                     zbc(l)+a(l)/bc(l),abs(q(l))/a(l)/sqrt(grav*a(l)/bc(l)),xx,zdam(l)
!
      end do ! l
!
      write(iunit,*)
!
      end subroutine
!
!----------------------------------------------------------------------
!
      subroutine prth
!
      use mgeometry, only : nx
      use marrays, only : xc,a,bc,zbc
      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,zbc(l),a(l)/bc(l)+zbc(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
      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 
