Código 2

      MODULE globais
          IMPLICIT NONE
          INTEGER :: N,L 
      END MODULE

      PROGRAM exemplo2_laplace
      USE globais
      IMPLICIT NONE
      INTEGER :: i,j,k,c,t,edge
      REAL :: dV
      INTEGER,DIMENSION(:,:),ALLOCATABLE :: masc
      REAL,DIMENSION(:,:),ALLOCATABLE :: V,V1,E1,E2,x,y 
      
      N = 10
      edge = 20*N
      L = edge + 1
      
      ALLOCATE(V(L,L),V1(L,L),E1(L,L),E2(L,L),x(L,L),y(L,L),masc(L,L)) 
      
      DO i=1,L
          DO j=1,L
              x(j,i) = -1.0 + (i-1)*0.1/N
              y(j,i) = -1.0 + (j-1)*0.1/N
          END DO
      END DO

! Calcula utilizando o algoritmo de Jacobi

      CALL inicia_zeros(E1)
      CALL inicia_zeros(E2)
      CALL initialize_V(V,masc)

      c = 0
      V1 = V

      CALL laplace_jacobi(V,V1,masc,dV,c)
      CALL electric_field(V,E1,E2)

      WRITE(*,*) 'Jacobi:',dV,c

      CALL dados(x,y,V,E1,E2)

! Calcula utilizando o algoritmo de Gauss-Siedel

      CALL inicia_zeros(E1)
      CALL inicia_zeros(E2)
      CALL initialize_V(V,masc)

      c = 0
      
      CALL laplace_gauss_siedel(V,masc,dV,c)
      CALL electric_field(V,E1,E2)

      WRITE(*,*) 'Gauss-Siedel:',dV,c

      CALL dados(x,y,V,E1,E2)

      END PROGRAM

      SUBROUTINE inicia_zeros(M)

! Inicializa uma matriz em zeros

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j
          REAL,DIMENSION(L,L) :: M
          DO i=1,L
              DO j=1,L
                  M(j,i) = 0.0
              END DO
          END DO
      END SUBROUTINE inicia_zeros
      
      SUBROUTINE initialize_V(V0,mascara)

! Inicialização do potencial

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j,m,o,p,q
          INTEGER,DIMENSION(L,L) :: mascara
          REAL,DIMENSION(L,L) :: V0

! Inicializa a mascara em 0

          DO m=1,L
              DO o=1,L
                  mascara(o,m) = 0
              END DO
          END DO

! Escreve 1's onde o V=1.0

          DO p=7*N+1,13*N+1
              DO q=7*N+1,13*N+1
                  mascara(q,p) = 1
              END DO
          END DO

! Escreve o potencial

          DO i=1,L
              DO j=1,L
                  IF (mascara(j,i) == 1) THEN
                      V0(j,i) = 1.0 
                  ELSE 
                      V0(j,i) = 0.0
                  END IF 
              END DO 
          END DO 
      END SUBROUTINE initialize_V

      SUBROUTINE update_V_jacobi(Vn,Vn1,mascara,deltaV)

! Calcula a atualização uma vez

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j
          REAL :: deltaV
          INTEGER,DIMENSION(L,L) :: mascara 
          REAL,DIMENSION(L,L) :: Vn,Vn1
          deltaV = 0.0 
          DO i=2,L-1
              DO j=2,L-1
                  IF (mascara(j,i) == 0) THEN
                      Vn1(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) + &
    Vn(j,i-1))/4.0
                      deltaV = deltaV + abs(Vn(j,i) - Vn1(j,i))
                  END IF
              END DO
          END DO
      END SUBROUTINE update_V_jacobi 

      SUBROUTINE laplace_jacobi(Vm,Vm1,mascara,delta_v,cont)

! Implementa o método de Jacobi

          USE globais
          IMPLICIT NONE
          INTEGER :: cont,x
          REAL :: delta_v,tol
          INTEGER,DIMENSION(L,L) :: mascara
          REAL,DIMENSION(L,L) :: Vm,Vm1
          tol = 1e-5*L*L 
          delta_v = 41.0
          DO WHILE (delta_v .gt. tol)
              CALL update_V_jacobi(Vm,Vm1,mascara,delta_v)
              CALL update_V_jacobi(Vm1,Vm,mascara,delta_v)
              cont = cont+1
          END DO
      END SUBROUTINE laplace_jacobi
      
      SUBROUTINE update_V_gauss_siedel(Vn,mascara,deltaV)

! Calcula a atualização uma vez

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j
          REAL :: deltaV,v1,v2
          INTEGER,DIMENSION(L,L) :: mascara
          REAL,DIMENSION(L,L) :: Vn
          deltaV = 0.0
          DO i=2,L-1
              DO j=2,L-1
                  IF (mascara(j,i) == 0) THEN
                      v1 = Vn(j,i)
                      Vn(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) +   &
    Vn(j,i-1))/4.0
                      v2 = Vn(j,i)
                      deltaV = deltaV + abs(v2 - v1)
                  END IF
              END DO
          END DO
      END SUBROUTINE update_V_gauss_siedel

      SUBROUTINE laplace_gauss_siedel(Vm,mascara,delta_v,cont)

! Implementa o método de Gauss_Siedel

          USE globais
          IMPLICIT NONE
          INTEGER :: cont
          REAL :: delta_v,tol
          INTEGER,DIMENSION(L,L) :: mascara
          REAL,DIMENSION(L,L) :: Vm
          tol = 1e-5*L*L 
          delta_v = 10.0
          DO WHILE (delta_v .gt. tol)
              CALL update_V_gauss_siedel(Vm,mascara,delta_v)
              CALL update_V_gauss_siedel(Vm,mascara,delta_v)
              cont = cont+1
          END DO
      END SUBROUTINE laplace_gauss_siedel
      
      SUBROUTINE electric_field(Vk,Ex,Ey)

! Calcula o campo elétrico a partir do potencial

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j
          REAL,DIMENSION(L,L) :: Vk,Ex,Ey

! Derivada para frente (bordas)

          DO i=1,L-1
              Ex(1,i) = - (Vk(1,i+1) - Vk(1,i)) 
              Ex(L,i) = - (Vk(L,i+1) - Vk(L,i))               
              Ey(i,1) = - (Vk(i+1,1) - Vk(i,1)) 
              Ey(i,L) = - (Vk(i+1,L) - Vk(i,L)) 
          END DO

! Derivada para trás (bordas)

          Ex(1,L) = - (Vk(1,L) - Vk(1,L-1))
          Ex(L,L) = - (Vk(L,L) - Vk(L,L-1))
          Ey(L,1) = - (Vk(L,1) - Vk(L-1,1))
          Ey(L,L) = - (Vk(L,L) - Vk(L-1,L))

! Derivada centrada

          DO i=2,L-1
              DO j=2,L-1 
                  Ex(j,i) = - (Vk(j,i+1) - Vk(j,i-1))/2.0 
                  Ey(j,i) = - (Vk(j+1,i) - Vk(j-1,i))/2.0 
              END DO
          END DO
      END SUBROUTINE electric_field

      SUBROUTINE dados(x1,y1,Vk,Ex,Ey)

! Imprime os dados em arquivos

          USE globais
          IMPLICIT NONE
          INTEGER :: i,j,k
          REAL,DIMENSION(L,L) :: Vk,Ex,Ey,x1,y1
          CHARACTER(28),DIMENSION(2) :: a
          CHARACTER(27),DIMENSION(2) :: b
          CHARACTER(32),DIMENSION(2) :: c

          a(1) = 'potencial_laplace_jacob2.dat'
          a(2) = 'potencial_laplace_gauss2.dat'
          b(1) = 'electric_laplace_jacob2.dat'
          b(2) = 'electric_laplace_gauss2.dat'
          c(1) = 'equipotencial_laplace_jacob2.dat'
          c(2) = 'equipotencial_laplace_gauss2.dat'

          DO k=1,2
          OPEN(10,FILE=a(k),STATUS='unknown',ACTION='write')
              DO i=1,L
                  DO j=1,L
                      WRITE(10,*) x1(j,i),y1(j,i),Vk(j,i)
                  END DO
                  WRITE(10,*) 
              END DO
          CLOSE(10)

          OPEN(20,FILE=b(k),STATUS='unknown',ACTION='write')
              DO i=1,L
                  DO j=1,L
                      WRITE(20,*) x1(j,i),y1(j,i),Ex(j,i),Ey(j,i)
                  END DO
                  WRITE(20,*) 
              END DO
          CLOSE(20)

          OPEN(30,FILE=c(k),STATUS='unknown',ACTION='write')
              DO i=1,L
                  DO j=1,L
                      IF ((Vk(j,i) .lt. 0.96) .and. (Vk(j,i) .gt. 0.94 & 
    ))THEN
                          WRITE(30,*) x1(j,i),y1(j,i)
                      ELSE IF ( (Vk(j,i) .lt. 0.76) .and. (Vk(j,i) .gt.& 
    0.74 )) THEN
                          WRITE(30,*) x1(j,i),y1(j,i)
                      ELSE IF ( (Vk(j,i) .lt. 0.56) .and. (Vk(j,i) .gt.&
    0.54 )) THEN
                          WRITE(30,*) x1(j,i),y1(j,i)
                      ELSE IF ( (Vk(j,i) .lt. 0.36) .and. (Vk(j,i) .gt.&
    0.34 )) THEN 
                          WRITE(30,*) x1(j,i),y1(j,i)
                      ELSE IF ( (Vk(j,i) .lt. 0.16) .and. (Vk(j,i) .gt.&
    0.14 )) THEN
                          WRITE(30,*) x1(j,i),y1(j,i)
                      END IF
                  END DO
              END DO
          CLOSE(30)
          END DO
          
      END SUBROUTINE dados