C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_phi.F,v 1.2 2013/06/21 20:49:50 jmc Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: STREAMICE_INIT_FIXED

C !INTERFACE:
      SUBROUTINE STREAMICE_INIT_PHI( myThid )

C     !DESCRIPTION:
C     Initialize STREAMICE nodal basis gradients for FEM solver

C     !USES:
      IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"
#include "GRID.h"

C     myThid ::  my Thread Id number
      INTEGER myThid
CEOP

C     !LOCAL VARIABLES:
C     === Local variables ===
      INTEGER bi, bj, i, j, xnode, ynode, xq, yq, m, n, p, kx, ky
      REAL gradx(2), grady(2)  ! gradients at quadrature points

C     here the terms used to calculate matrix terms in the
C     velocity solve are initialized
C
C     this is a quasi-finite element method; the gradient
C     of the basis functions are approximated based on knowledge
C     of the grid
C
C     Dphi (i,j,bi,bj,m,n,p):
C       gradient (in p-direction) of nodal basis function in
C       cell (i,j) on thread (bi,bj) which is centered on node m,
C       at quadrature point n
C
C    %  3 - 4
C    %  |   |
C    %  1 - 2
C
C     NOTE 2x2 quadrature is hardcoded - might make it specifiable through CPP
C
C     this will not be updated in overlap cells - so we extend it as far as we can

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly-1
         DO i=1-Olx,sNx+Olx-1

          DO xq = 1,2
           gradx(xq) = Xquad(3-xq) * recip_dxG (i,j,bi,bj) +
     &                 Xquad(xq) * recip_dxG (i+1,j,bi,bj)
           grady(xq) = Xquad(3-xq) * recip_dyG (i,j,bi,bj) +
     &                 Xquad(xq) * recip_dyG (i,j+1,bi,bj)
          ENDDO

          DO n = 1,4

           xq = 2 - mod(n,2)
           yq = floor ((n+1)/2.0)

           DO m = 1,4

            xnode = 2 - mod(m,2)
            ynode = floor ((m+1)/2.0)

            kx = 1 ; ky = 1
            if (xq.eq.xnode) kx = 2
            if (yq.eq.ynode) ky = 2


            Dphi (i,j,bi,bj,m,n,1) =
     &       (2*xnode-3) * Xquad(ky) * gradx(yq)
            Dphi (i,j,bi,bj,m,n,2) =
     &       (2*ynode-3) * Xquad(kx) * grady(xq)

           ENDDO

           grid_jacq_streamice (i,j,bi,bj,n) =
     &      (Xquad(3-xq)*dyG(i,j,bi,bj) + Xquad(xq)*dyG(i+1,j,bi,bj)) *
     &      (Xquad(3-yq)*dxG(i,j,bi,bj) + Xquad(yq)*dxG(i,j+1,bi,bj))

          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      RETURN
      END