C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/mlosch/cg2d_sr_bench/code/ini_rhs.F,v 1.1 2009/11/29 11:02:03 mlosch Exp $
C $Name:  $
C $Id: ini_rhs.F,v 1.1 2009/11/29 11:02:03 mlosch Exp $
#include "CPP_OPTIONS.h"
CStartOfInterface
      SUBROUTINE INI_RHS(
     I                   cg2d_b,
     O                   cg2d_x1,
     O                   cg2d_x2,
     I                   myThid)
C     /==========================================================\
C     | SUBROUTINE INI_RHS                                       |
C     | o Initialise 2d conjugate gradient solver right-hand side|
C     |==========================================================|
C     | Set a source term b in Ax = b (1)                        |
C     | We solve (1) with neumann bc's whic means that b must    |
C     | integrate out to zero over the whole domain. If b does   |
C     | not integrate out to zero then the solution will         |
C     | converge, but to a non-zero final residual.              |
C     \==========================================================/
      IMPLICIT NONE

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

C     === Routine arguments ===
CEndOFInterface

C     === Local variables ===
C     iG, jG     - Global coordinate index
C     faceArea   - Temporary used to hold cell face areas.
C     I,J,K      - Loop counters
C     iMid, jMid - Global coords of mid-point of model domain
      INTEGER  I,  J
      INTEGER bi, bj
      INTEGER iG, jG
      INTEGER iMid, jMid
      INTEGER iQ, i3Q
      INTEGER myThid
      _RL  cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL  cg2d_x1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL  cg2d_x2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL bsum, bsumTile(nSx,nSy)
      _RL port_rand
      INTEGER numPoints

C--   Get model global domain mid-point
      iMid = Nx/2
      iQ   = Nx/4
      i3Q  = 3*Nx/4
      jMid = Ny/2

C--   Set dummy source term for elliptic equation
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          cg2d_b(I,J,bi,bj) = 0. _d 0
          cg2d_x1(I,J,bi,bj) = 0. _d 0
          cg2d_x2(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      
C      bsum = 0. _d 0
C      numPoints = NX*NY
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
C        bsumTile(bi,bj) = 0. _d 0
        DO J=1,sNy
         DO I=1,sNx
          cg2d_b(I,J,bi,bj) = 0. _d 0
C--       Set +/-1 source function in middle of domain.
          iG = myXGlobalLo + sNx*(bi-1) + I - 1
          jG = myYGlobalLo + sNy*(bj-1) + J - 1
C          write(*,'(i2,2i3,2i5,2x,a,2i5)')mythid,bi,bj,i,j,'iG, jG =',iG,jG
          IF ( iG .EQ.  Nx/2   .AND. jG .EQ. Ny/2 ) 
     &           cg2d_b(I,J,bi,bj) =  1._d 0
          IF ( iG .EQ.  Nx/2+1  .AND. jG .EQ. Ny/2 ) 
     &           cg2d_b(I,J,bi,bj) = -1. _d 0
C          cg2d_b(I,J,bi,bj) = port_rand(-1.)
C          bsumTile(bi,bj) = bsumTile(bi,bj) + cg2d_b(I,J,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      
C      CALL GLOBAL_SUM_TILE_RL(bsumTile,bsum,myThid)
      
!       DO bj=myByLo(myThid),myByHi(myThid)
!        DO bi=myBxLo(myThid),myBxHi(myThid)
!         DO J=1,sNy
!          DO I=1,sNx
!           cg2d_b(I,J,bi,bj) = cg2d_b(I,J,bi,bj) - bsum/numPoints
!          ENDDO
!         ENDDO
!        ENDDO
!       ENDDO

C--   Update overlap regions synchronously
      CALL EXCH_XY_RL(cg2d_b,myThid)
      CALL EXCH_XY_RL(cg2d_x1,myThid)
      CALL EXCH_XY_RL(cg2d_x2,myThid)
C
      RETURN
      END