module MSI

use msflib
use msfwin

use variaveis
use Coef
	
contains

subroutine lu2d5(b,dl,du,nij,ni,nj)
      implicit real*8(a-h,o-z)
      integer,intent(in) :: nij, ni, nj
      real*8,intent(in), dimension(nij,5) :: b ! coeficientes
      real*8,intent(out),dimension(nij,4) :: dl
      real*8,intent(out),dimension(nij,3) :: du
!      dimension b(nij,5),dl(nij,4),du(nij,3)
      nj1=nj-1
      ni1=ni-1
      j1y=1
      j1x=2
      j2y=5
      j2x=4
      jp=3
      alpha=0.5d0
      alph0=alpha
!      do 80 j=2,nj1
!      do 80 i=2,ni1
      do 80 j=1,nj
      do 80 i=1,ni
      np=(j-1)*ni+i
      npw=np-1
      nps=np-ni
      npsp1=nps+1
!      if(i.eq.2.or.i.eq.ni1.or.j.eq.2.or.j.eq.nj1) alpha=0.0d0
      if(i.eq.1.or.i.eq.ni.or.j.eq.1.or.j.eq.nj) alpha=0.0d0
      if(nps.lt.0) go to 10
      go to 40
   10 dl(np,4)=b(np,j1y)
      dl(np,3)=0.0d0
      if(np.lt.2) go to 20
      go to 30
   20 dl(np,2)=b(np,j1x)
      dl(np,1)=b(np,jp)
      du(np,1)=b(np,j2x)/dl(np,1)
      du(np,2)=0.0d0
      du(np,3)=b(np,j2y)/dl(np,1)
      go to 70
   30 dl(np,2)=b(np,j1x)/(1.0d0+2.0d0*alpha*du(npw,2))
      dl(np,1)=b(np,jp)-dl(np,2)*du(npw,1)+alpha*2.0d0*dl(np,2)*du(npw,2)
      du(np,1)=b(np,j2x)/dl(np,1)
      du(np,2)=-dl(np,2)*du(npw,3)/dl(np,1)
      du(np,3)=(b(np,j2y)-alpha*dl(np,2)*du(npw,2))/dl(np,1)
      go to 70
   40 if(nps.eq.0) go to 50
      go to 60
   50 dl(np,4)=b(np,j1y)
      dl(np,3)=0.0d0
      dl(np,2)=b(np,j1x)/(1.0d0+2.0d0*alpha*du(npw,2))
      fi1=dl(np,3)*du(npsp1,1)
      fi4=dl(np,2)*du(npw,2)
      dl(np,1)=b(np,jp)-dl(np,3)*du(npsp1,2)-dl(np,2)*du(npw,1)+ &
               alpha*(2.0d0*fi1+2.0d0*fi4)
      du(np,1)=(b(np,j2x)-dl(np,3)*du(npsp1,3)-2.0d0*alpha*fi1)/dl(np,1)
      du(np,2)=-dl(np,2)*du(npw,3)/dl(np,1)
      du(np,3)=(b(np,j2y)-alpha*fi4)/dl(np,1)
      go to 70
   60 dl(np,4)=b(np,j1y)/(1.0d0- alpha*du(nps,1)*du(npsp1,1))
      dl(np,3)=-dl(np,4)*du(nps,1)
      dl(np,2)=(b(np,j1x)-dl(np,4)*du(nps,2))/(1.0d0+2.0d0*alpha*du(npw,2))
      fi1=dl(np,3)*du(npsp1,1)
      fi4=dl(np,2)*du(npw,2)
      dl(np,1)=b(np,jp)-dl(np,4)*du(nps,3)-dl(np,3)*du(npsp1,2) &
              -dl(np,2)*du(npw,1)+alpha*(2.0d0*fi1+2.0d0*fi4)
      du(np,1)=(b(np,j2x)-dl(np,3)*du(npsp1,3)-2.0d0*alpha*fi1)/dl(np,1)
      du(np,2)=-dl(np,2)*du(npw,3)/dl(np,1)
      du(np,3)=(b(np,j2y)-alpha*fi4)/dl(np,1)
   70 continue
      alpha=alph0
   80 continue
!      return
end subroutine lu2d5
		
! --------------------------------------------------------------------------

subroutine fb2d5(b,dl,du,nj,ni,nij,var,gama,nitm,rhs)
      implicit real*8(a-h,o-z)
      integer,intent(in) :: nij, ni, nj
      integer,intent(in) :: nitm
      real*8,intent(in)  :: gama
      real*8,intent(in),dimension(nij)   :: rhs ! fonte
      real*8,intent(in),dimension(nij,5) :: b ! coeficientes
      real*8,intent(in),dimension(nij,4) :: dl
      real*8,intent(in),dimension(nij,3) :: du
      real*8,intent(inout),dimension(nij)  :: var ! solução numérica
      real*8,dimension(:),allocatable :: r,c,diff
      allocate ( r(nij), c(nij), diff(nij) )
!      dimension b(nij,5),dl(nij,4),du(nij,3)
!      dimension r(nij),c(nij),diff(nij),var(nij),rhs(nij)
      ni1=ni-1
      nj1=nj-1
      j1x=2
      j2x=4
      jp=3
      j1y=1
      j2y=5
!c
!c       calculo do residuo com a estimativa inicial
!c
      rp=0.0d0
!      do 50 j=2,nj1
!      do 50 i=2,ni1
      do 50 j=1,nj
      do 50 i=1,ni
      np=(j-1)*ni + i
      npe=np+1
      npw=np-1
      nps=np-ni
      npn=np+ni
      if(np.lt.2) go to 5
      go to 10
    5 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn))
      go to 45
   10 if(np.le.ni) go to 15
      go to 20
   15 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn)+b(np,j1x)*var(npw))
      go to 45
   20 if(np.eq.nij) go to 25
      go to 30
   25 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j1y)*var(nps) &
           +b(np,j1x)*var(npw))
      go to 45
   30 if(np.gt.(nij-ni)) go to 35
      go to 40
   35 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j1x)*var(npw)+b(np,j1y)*var(nps))
      go to 45
   40 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn)+b(np,j1x)*var(npw)+b(np,j1y)*var(nps))
   45 continue
      rp=rp+r(np)*r(np)
   50 continue
!c     call print1(r,nj,ni,nij)
      rp=dsqrt(rp)
!c              fim do calculo do residuo inicial
      gamarp=gama*rp
      nit=0
    1 nit=nit+1
!      do 85 j=2,nj1
!      do 85 i=2,ni1
      do 85 j=1,nj
      do 85 i=1,ni
      np=(j-1)*ni+i
      npw=np-1
      nps=np-ni
      npsp1=nps+1
      if(nps.lt.0) go to 55
      go to 70
   55 if(np.lt.2) go to 60
      go to 65
   60 c(np)=r(np)/dl(np,1)
      go to 85
   65 c(np)=(r(np)-dl(np,2)*c(npw))/dl(np,1)
      go to 85
   70 if(nps.eq.0) go to 75
      go to 80
   75 c(np)=(r(np)-dl(np,2)*c(npw)-dl(np,3)*c(npsp1))/dl(np,1)
      go to 85
   80 c(np)=(r(np)-dl(np,2)*c(npw)-dl(np,3)*c(npsp1)-dl(np,4)* &
            c(nps))/dl(np,1)
   85 continue
!c         diff calculation
      do 125 j=1,nij
      np=nij-j+1
      npn=np+ni
      npnm1=npn-1
      npe=np+1
      if(np.eq.nij) go to 90
      go to 95
   90 diff(np)=c(np)
      go to 120
   95 if(npn.gt.nij) go to 100
      go to 115
  100 if(npnm1.eq.nij) go to 105
      go to 110
  105 diff(np)=c(np)-du(np,1)*diff(npe)-du(np,2)*diff(npnm1)
      go to 120
  110 diff(np)=c(np)-du(np,1)*diff(npe)
      go to 120
  115 diff(np)=c(np)-du(np,1)*diff(npe)-du(np,2)*diff(npnm1)- &
               du(np,3)*diff(npn)
  120 continue
      var(np)=var(np)+diff(np)
  125 continue
!c     call print1(var,nj,ni,nij)
      rp=0.0d0
!      do 350 j=2,nj1
!      do 350 i=2,ni1
      do 350 j=1,nj
      do 350 i=1,ni
      np=(j-1)*ni + i
      npe=np+1
      npw=np-1
      nps=np-ni
      npn=np+ni
      if(np.lt.2) go to 305
      go to 310
  305 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn))
      go to 345
  310 if(np.le.ni) go to 315
      go to 320
  315 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn)+b(np,j1x)*var(npw))
      go to 345
  320 if(np.eq.nij) go to 325
      go to 330
  325 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j1y)*var(nps) &
           +b(np,j1x)*var(npw))
      go to 345
  330 if(np.gt.(nij-ni)) go to 335
      go to 340
  335 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j1x)*var(npw)+b(np,j1y)*var(nps))
      go to 345
  340 r(np)=rhs(np)-(b(np,jp)*var(np)+b(np,j2x)*var(npe) &
           +b(np,j2y)*var(npn)+b(np,j1x)*var(npw)+b(np,j1y)*var(nps))
  345 continue
      rp=rp+r(np)*r(np)
  350 continue
!c     call print1(r,nj,ni,nij)
      rp=dsqrt(rp)
      if(nit.eq.nitm) go to 230
      if(rp.gt.gamarp) go to 1
  230 continue
!      write(6,240)nit,rp,gamarp
!  240 format(5x,'nit = ',i3,5x,'rp = ',e14.6,5x,'gamarp0 = ',e14.5)
!      write(10,240) itex,nit, gamarp, rp
!  240 format(1x,'it.externa = ',i6,5x,'it.MSI = ',i3,5x, &
!             'resíduo a sair = ',1pe8.2,5x,'resíduo final = ',1pe8.2)
!      write(10,240) nit, gamarp, rp
!  240 format(1x,'it.MSI = ',i3,5x, &
!             'resíduo a sair = ',1pe8.2,5x,'resíduo final = ',1pe8.2)
      deallocate ( r, c, diff )
!      return

end subroutine fb2d5

! ==================================================================

end module MSI