C	$Id: ini_cg2d.F,v 1.2 2006/05/12 22:23:42 ce107 Exp $	
CStartOfInterface
      SUBROUTINE INI_CG2D
C     /==========================================================\
C     | SUBROUTINE INI_CG2D                                      |
C     | o Initialise 2d conjugate gradient solver operators.     |
C     |==========================================================|
C     | These arrays are purely a function of the basin geom.    |
C     | We set then here once and them use then repeatedly.      |
C     | In this example we hard code a periodic channel with a   |
C     | grid spacing on 1.                                       |
C     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "CG2D.h"

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface

C     === Local variables ===
C     xG, yG - Global coordinate location.
C     zG
C     iG, jG - Global coordinate index
C     bi,bj  - Loop counters
C     faceArea - Temporary used to hold cell face areas.
C     I,J,K
      INTEGER  I,  J, K
      INTEGER iG, jG

C--   Initialise laplace operator
C     aW2d: integral Ax/dX
C     aS2d: integral Ay/dY
      DO J=1-OLy,sNy+OLy
       DO I=1-OLx,sNx+OLx
        iG = myXGlobalLo+I-1
        jG = myYGlobalLo+J-1
        aW2d(I,J) = 1. _d 0
        aS2d(I,J) = 1. _d 0
        IF ( jG .EQ. 1 .OR. jG .EQ. nY+1 ) THEN
         aS2d(I,J) = 0. _d 0
        ENDIF
       ENDDO
      ENDDO
C     CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D')
C     CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D')

C--   Initialise preconditioner
      DO J=1,sNy
       DO I=1,sNx
        pC(I,J) = 1. _d 0
        IF (
     &   aW2d(I,J) + aW2d(I+1,J) 
     &  +aS2d(I,J) + aS2D(I,J+1)
     &   .EQ. 0. _d 0
     &     )  pC(I,J) = 0. _d 0
        pW(I,J) = 0. _d 0
        pS(I,J) = 0. _d 0
       ENDDO
      ENDDO
C--   Update overlap regions
      CALL EXCH_XY_R8(pC)
      CALL EXCH_XY_R8(pW)
      CALL EXCH_XY_R8(pS)

C--   Set default values for initial guess
      DO J=1,sNy
       DO I=1,sNx
        cg2d_x(I,J) = 0. _d 0
       ENDDO
      ENDDO
C--   Update overlap regions
      CALL EXCH_XY_R8(cg2d_x)

      RETURN
      END