!**************************************************************************
!* Solve Laplace Equation by relaxation Method: d2T/dx2 + d2T/dy2 = 0 (3) *
!* ---------------------------------------------------------------------- *
!* Example #3: Determine the temperatures T(x,y) in a domain [x,y] made   *
!*             of one square of length C with the following conditions:   *
!*             along OC: T=0, along AB: T=1, along OA and CB: dT/dy = 0   *
!*             (see sketch below).                                        *
!*             Main parameters:                                           *
!*             C: length of square domain.                                *
!*             Number of subdivisions 8 in ox and 8 in Oy.                *
!*             integration steps dx=C/8.                                  *
!*             weight coefficient w = 1.3                                 *
!*             Maximum residual error ER = 0.001                          *
!*             Maximum number of iterations IT = 100                      *
!* ---------------------------------------------------------------------- *
!* SAMPLE RUN:                                                            *
!*                                                                        *
!* Niter  Max. Residual     W                                             *
!*  43    0.000860         1.30                                           *
!* Temperature                                                            *
!* 0.000  0.125  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.125  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.374  0.499  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.373  0.498  0.624  0.749  0.875  1.000          *
!* 0.000  0.124  0.249  0.373  0.498  0.624  0.749  0.874  1.000          *
!*                                                                        *
!* ---------------------------------------------------------------------- *
!* REFERENCE:  "Méthode de calcul numérique- Tome 2 - Programmes en Basic *
!*              et en Pascal By Claude Nowakowski, Edition du P.S.I.,     *
!*              1984" [BIBLI 04].                                         *
!*                                                                        *
!*                                   F90 Release By J-P Moreau, Paris.    *
!*                                          (www.jpmoreau.fr)             *
!**************************************************************************
! Note: Here, what is new is the limit conditions involving partial derivatives dT/dy.
!       Inside domain, the classical scheme is applied to calculate the residual value at current
!       point (i,j): 4*T(i,j)-T(i-1,j)-T(i+1,j)-Ti,j-1)-T(i,j+1). Must be near zero (< ER).
!       For limit lines OA and CB, the partial derivative is approximated by:
!                       dT/dy = (T(i,j+1)-T(i,j-1))/(2*dy)
!       But the point T(i,j+1) for segment CB is out of the domain and becomes a fictitious point,
!       so we apply the following scheme:  r=4*T(i,j)-2*T(i,j-1)-T(i-1,j)-T(i+1,j).
!       For the other segment OA, it becomes: r=4*T(i,j)-2*T(i,j+1)-T(i-1,j)-T(i+1,j).
!       The exact solution for this problem is T=K*x (0<x<C), here k=0.25.
!       This allows verifying the results accuracy.
!-------------------------------------------------------------------------------------------------         
Program laplace2

implicit real*8 a-h,o-z

DIMENSION T(0:8, 0:8)

                      !           dT/dy=0
      IT = 100        !       C-------------B
      ER = 0.001      !       |             |
      OM = 1.3        !       |             |
      C  = 4          !   T=0 |             |T=1
      NM = 8          !       |             |
                      !       |             |
                      !       O-------------A
                      !           dT/dy=0

  !Fix limit conditions
  DO i = 0, NM
    DO j = 0, NM
      T(i, j) = 0.D0
    END DO
  END DO
  DO j = 0, NM
    T(NM, j) = 1.D0
  END DO

  n1 = NM - 1
  !main calculation loop
  DO l = 1, IT
    rm = 0.D0
    j = 0
    DO i = 1, n1
      call Schema1(NM,i,j,T,r,rm,OM)
    END DO
    DO j = 1, n1
      DO i = 1, n1
        call Schema(NM,i,j,T,r,rm,OM)
      END DO
    END DO
    j = NM
    DO i = 1, n1
      call Schema2(NM,i,j,T,r,rm,OM)
    END DO
    IF (rm < ER) GOTO 10
  END DO
  !print results
  WRITE(*,30)  IT
  WRITE(*,31)  rm
  GOTO 20
10 PRINT *,' '
  PRINT *,' Niter  Max. Residual     W'
  WRITE(*,40) l, rm, OM
  PRINT *,' Temperature'
  DO j = NM, 0, -1
    DO i = 0, NM
      WRITE(*,41,advance='no')  T(i,j)
    END DO
    PRINT *,' '
  END DO
20 PRINT *,' '

  STOP

30 FORMAT(' No convergence after ',I3,' iterations.')
31 FORMAT(' Residual rm = ',F6.4)
40 FORMAT(I5,F12.6,F13.2)
41 FORMAT(F7.3)

  END


!calculate residual for inside square domain
Subroutine Schema(NM,i,j,T,r,rm,OM)
real*8 T(0:NM, 0:NM), r,rm,OM 
  r = T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)
  T(i, j) = T(i, j) + 0.25 * OM * r
  IF (DABS(r) > rm) rm = DABS(r)
  Return
End

!calculate residual for lower limit
Subroutine Schema1(NM,i,j,T,r,rm,OM)
real*8 T(0:NM, 0:NM), r,rm,OM 
  r = T(i + 1, j) + T(i - 1, j) + 2 * T(i, j + 1) - 4 * T(i, j)
  T(i, j) = T(i, j) + 0.25 * OM * r
  IF (DABS(r) > rm) rm = DABS(r)
  Return
End

!calculate residual for upper limit
Subroutine Schema2(NM,i,j,T,r,rm,OM)
real*8 T(0:NM, 0:NM), r,rm,OM 
  r = T(i + 1, j) + T(i - 1, j) + 2 * T(i, j - 1) - 4 * T(i, j)
  T(i, j) = T(i, j) + 0.25 * OM * r
  IF (DABS(r) > rm) rm = DABS(r)
  Return
End

!end of file laplace2.f90