       real S,I,R,C,D
       real gama,mu,beta
       real t,dt,tm,pm
       real ni,s0,i0,r0,tt
       real di,mu1,t1,III
       integer ii1,j00,nmm
       integer j1
       parameter(nmm=1000)
        real tempos(nmm),dados(nmm)
         character*100 file100,file200
       real s10(0:nmm),r10(0:nmm),i10(0:nmm)
      real s20(nmm,nmm),r20(nmm,nmm),i20(nmm,nmm),minn
      real diff(nmm,nmm),betaa(nmm,nmm)
      real beta100(nmm),gamaa(nmm,nmm),gama100(nmm)
      integer nd,ii,nbb,ntt,itt,i00,idd
      real dadosm(0:nmm),mort(0:nmm),mmm100
      real aux1000,r1,r2,r3
      real R01(0:nmm),R02(0:nmm),R03(0:nmm),R04(0:nmm)
      integer i000,j000
      integer ii000,jj000 
      real beta200(nmm),gama200(nmm)
        real beta300(nmm),gama300(nmm)
        real beta400(nmm),gama400(nmm)
        real RRR(0:nmm),errom
      
      i000=1
      j000=1
      ii000=1
      jj000=1
       
       
       dt=0.00001
       tm=71.
       pm=220000000.
      
       ! ENTRA COM O TAMANHO DOS DADOS DE ESTRADA (o numero de dias depois de 22 de março) 
       write(*,*)"n"
       read(*,*)nd
       
       ! Ler as tabelas de casos confirmados e casos letais 
1933         format(a100)   
      ! write(*,*)"dados de entrada casos confirmados"
      ! read(*,1933)file100 
       file100="tabelabrasil.txt2"
        open(12,file=file100) 
       !    write(*,*)"dados de entrada casos letais"
       !read(*,1933)file200 
        file200="mortesbrasil.txt2"
       open(122,file=file200) 
        
       do ii=1,nd
       read(12,*)tempos(ii),dados(ii)
       read(122,*)tempos(ii),dadosm(ii) 
       enddo 
        idd=1
          di=tempos(idd)
       ni=dados(idd) 
       
       
       ! Inicia S, I e R para t=0
        s10(0)=1-2.*ni/pm
        r10(0)=ni/pm
        i10(0)=ni/pm
      
               !Criar nomes para os arquivos de saida
       open(15,file="BAinfectados.dat") ! Arquivo principal com os resultados de I ao longo do tempo 
       
        open(110,file="BAR0medioxtempo.dat")
        
         open(1100,file="BAR0xtempo.dat") 
                 
       open(17,file="BAdadosreaiscasosconfirmados.dat") 
       open(177,file="BAdadosreaiscasosletais.dat") 
       
               open(220,file="BAcalosletais.dat") 
        
       !gama=0.09
   
       mu=0.015
       mu1=0.035
       ntt=0 
     
      t1=0.
       itt=0 
       
       ! Inicia a simulacao partindo de t=0 (consideramos 22 de março no Brasil ) até o dia final (o dia presente) 
       do temp=1.0e0,float(nd-1),1.0e0
         idd=idd+1           
          
        ! Roda gama e beta dentro da faixa escolhida 
          naa=0
          do gama=0.08,0.115,0.005 
          naa=naa+1
          nbb=0
        do beta=0.08,0.5,0.005 ! beta 
        nbb=nbb+1
             
           s0=s10(itt)
        r0=r10(itt)
        i0=i10(itt) 
         
       
        
         do t=dt+t1,t1+1.,dt !tempo
         
       
                              
          !if (t.gt.50) beta=0.2
          ! if (t.gt.75) beta=0.3
         !Integrar o modelo SIR 
         S=s0-dt*beta*i0*s0 
         
         I=i0+dt*(beta*i0*s0-gama*i0)
         
         R=r0+dt*gama*i0 
         
         
         s0=S
         
         i0=I
         r0=R
         
         
         
         enddo !tempo 
         
         s20(nbb,naa)=S
         r20(nbb,naa)=R
         i20(nbb,naa)=I 
         
         III=I*pm
         
         
         
         betaa(nbb,naa)=beta
         gamaa(nbb,naa)=gama
         diff(nbb,naa)=abs(III-dados(idd))
         enddo ! beta
         enddo ! gamma
         minn=1e+20
         do ii1=1,naa
         do ii=1,nbb 
         if (diff(ii,ii1).lt.minn)then 
         minn=diff(ii,ii1)
         ! Escolhe qual valores de gama e beta 
         ! fornecem valores de I proximo dos dados reais 
         !Escolhemos 3 conjuntos de valores de gama,beta 
         ! aqueles 3 que ficam mais proximos 
         ! o <beta> e <gama> é medido atraves d euma média destes 3 valores
         iii000=ii000 
         jjj000=jj000 
         
         ii000=i000 
         jj000=j000
         
         i000=i00 
         j000=j00 
         
         
         
         i00=ii
         j00=ii1
         
         endif 
         enddo 
         enddo 
         
         itt=itt+1
         
         s10(itt)=s20(i00,j00)
         r10(itt)=r20(i00,j00)
         i10(itt)=i20(i00,j00) 
         beta100(itt)=betaa(i00,j00)
         gama100(itt)=gamaa(i00,j00)
         
          beta200(itt)=betaa(i000,j000)
         gama200(itt)=gamaa(i000,j000)
         
           beta300(itt)=betaa(ii000,jj000)
         gama300(itt)=gamaa(ii000,jj000)
         
          beta400(itt)=betaa(iii000,jjj000)
         gama400(itt)=gamaa(iii000,jjj000)
         
         !write(*,*)itt
         write(15,*)temp,i10(itt)*pm
         !write(,dados(idd)
           write(17,*)temp,dados(idd) 
            write(177,*)temp,dadosm(idd)  
       !write(110,*)temp,beta100(itt)/gama100(itt),beta100(itt),gama100(itt)
         mort(itt)=dadosm(idd)/dados(idd)
         write(220,*)temp,mort(itt)*i10(itt)*pm
          t1=t1+1.0e0
         R01(itt)=beta100(itt)/gama100(itt)
         R02(itt)=beta200(itt)/gama200(itt)
         R03(itt)=beta300(itt)/gama300(itt)
         R04(itt)=beta400(itt)/gama400(itt)
         
         enddo ! tempo 
        
        !R0 ao longo do tempo  
         do ii=1,itt
         RRR(ii)=R01(ii)+R02(ii)+R03(ii)+R04(ii)
         RRR(ii)=RRR(ii)*0.25
         enddo 
          errom=0.
         
         do ii=1,itt 
        write(1100,*)ii,RRR(ii),beta100(itt),gama100(itt)
         errom=errom+(RRR(ii)-R01(ii))**2.
           errom=errom+(RRR(ii)-R02(ii))**2.
           errom=errom+(RRR(ii)-R03(ii))**2.
             errom=errom+(RRR(ii)-R04(ii))**2.
             
             enddo 
             
             errom=errom/(float(4*itt))
            write(*,*)errom/RRR(itt) 
            
      do ii=4,itt-3,4
      r1=RRR(ii)+RRR(ii-1)+RRR(ii-2)+RRR(ii-3)
     
      r1=r1*0.25
      
      r2=beta100(ii)+beta100(ii-1)
      r2=r2+beta100(ii-2)+beta100(ii-3)
      r2=r2*0.25
      
      r3=gama100(ii)+gama100(ii-1)
       r3=r3+gama100(ii-2)+gama100(ii-3)
      r3=r3*0.25
      
      write(110,*)ii,r1,r2,r3
      enddo
      
     
     
     
       ii=itt 
     
       r1=RRR(ii)+RRR(ii-1)+RRR(ii-2)+RRR(ii-3)
     
      r1=r1*0.25
      
       
      r2=beta100(ii)+beta100(ii-1)
      r2=r2+beta100(ii-2)+beta100(ii-3)
      r2=r2*0.25
      
      r3=gama100(ii)+gama100(ii-1)
       r3=r3+gama100(ii-2)+gama100(ii-3)
      r3=r3*0.25
      
      write(110,*)ii,r1,r2,r3
     
      write(*,*)r1,r2,r3
     
       gama=r3
       beta=0.999*r3*r1
         
        
          mmm100=(mort(itt)+mort(itt-1)+mort(itt-2))
          mmm100=mmm100+mort(itt-3)+mort(itt-4)
          mmm100=mmm100/5.
          ! Aqui começa a previsão para 10 dias 
          ! Usando o ultimo valor de beta e gama 
          ! calculado como sendo a média da ultima semana de dados 
          t1=float(nd-1)
           do temp=float(nd),float(nd)+10.,1.0e0 !tempo 2 
         
            s0=s10(itt)
        r0=r10(itt)
        i0=i10(itt) 
         
         
          
         do t=dt+t1,t1+1.,dt !tempo X
         
       
                              
          !if (t.gt.50) beta=0.2
          ! if (t.gt.75) beta=0.3
         !integra as equações  
         S=s0-dt*beta*i0*s0 
         
         I=i0+dt*(beta*i0*s0-gama*i0)
         
         R=r0+dt*gama*i0 
         
         
         s0=S
         
         i0=I
         r0=R
         
         enddo !tempo X
         
           itt=itt+1
           s10(itt)=s0
           r10(itt)=r0
           i10(itt)=i0
           t1=t1+1.0e0
               write(15,*)temp,i10(itt)*pm
        
          write(220,*)temp,mmm100*i10(itt)*pm
          write(110,*)temp,beta/gama
           !if (temp.gt.100) then 
           !if (beta.gt.0.09) beta=0.991*beta
           !endif 
           !O beta o Brasil tem uma leva diminuição ao longo das semanas
           ! mantivemos esta diminuição usando este truque de multiplicar diariamente 
           !o beta pro uma pequena constante ligeiramente menor que 1 
           beta=0.999*beta
         mmm100=0.9999*mmm100
         
         enddo ! loop do temp 2
         
         end
         
       
