C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_cg_functions.F,v 1.3 2015/09/29 15:56:27 dgoldberg Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_CG_ACTION( myThid, 
     O    uret, 
     O    vret, 
     I    u, 
     I    v, 
     I    is, ie, js, je )
C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "STREAMICE.h"      
#include "STREAMICE_CG.h"      

C     !INPUT/OUTPUT ARGUMENTS
C     uret, vret - result of matrix operating on u, v
C     is, ie, js, je - starting and ending cells
      INTEGER myThid
      _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER is, ie, js, je

#ifdef ALLOW_STREAMICE

C the linear action of the matrix on (u,v) with triangular finite elements
C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
C but this may change pursuant to conversations with others
C
C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
C     in order to make less frequent halo updates
C isym = 1 if grid is symmetric, 0 o.w.

C the linear action of the matrix on (u,v) with triangular finite elements
C Phi has the form
C Phi (i,j,k,q) - applies to cell i,j

C      3 - 4
C      |   |
C      1 - 2

C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q 
C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear

C     !LOCAL VARIABLES:
C     == Local variables == 
      INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n,Gi,Gj
      _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
      _RL Ucell (2,2)
      _RL Vcell (2,2)
      _RL Hcell (2,2)
      _RL phival(2,2)

      uret(1,1,1,1) = uret(1,1,1,1)
      vret(1,1,1,1) = vret(1,1,1,1)

      DO j = js, je
       DO i = is, ie
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)

         Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
         Gj = (myYGlobalLo-1)+(bj-1)*sNy+j

          IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
           DO iq = 1,2 
            DO jq = 1,2

            n = 2*(jq-1)+iq


            uq = u(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) + 
     &       u(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) + 
     &       u(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) + 
     &       u(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
            vq = v(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) + 
     &       v(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) + 
     &       v(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) + 
     &       v(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
            ux = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) + 
     &       u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) + 
     &       u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) + 
     &       u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
            uy = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) + 
     &       u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) + 
     &       u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) + 
     &       u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
            vx = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) + 
     &       v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) + 
     &       v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) + 
     &       v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
            vy = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) + 
     &       v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) + 
     &       v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) + 
     &       v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
            exx = ux + k1AtC_str(i,j,bi,bj)*vq
            eyy = vy + k2AtC_str(i,j,bi,bj)*uq
            exy = .5*(uy+vx) + 
     &       k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq

            do inode = 1,2 
             do jnode = 1,2

             m = 2*(jnode-1)+inode
             ilq = 1 
             jlq = 1
             if (inode.eq.iq) ilq = 2
             if (jnode.eq.jq) jlq = 2 
             phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)

             if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then            

              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) + 
     &          DPhi(i,j,bi,bj,m,n,2)*(2*exy))


              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
     &          4*0.5*k1AtC_str(i,j,bi,bj)*exy)


              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * 
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * uq


             endif
    
             if (STREAMICE_vmask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) + 
     &          DPhi(i,j,bi,bj,m,n,1)*(2*exy))
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
     &          4*0.5*k2AtC_str(i,j,bi,bj)*exy)
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * 
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * vq
              
             endif
            enddo 
            enddo

           enddo            
           enddo
c-- STREAMICE_hmask
          endif

         enddo
        enddo
       enddo
      enddo

#endif
      RETURN
      END SUBROUTINE

      SUBROUTINE STREAMICE_CG_MAKE_A( myThid )
C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "STREAMICE.h"      
#include "STREAMICE_CG.h"      

C     !INPUT/OUTPUT ARGUMENTS
C     uret, vret - result of matrix operating on u, v
C     is, ie, js, je - starting and ending cells
      INTEGER myThid

#ifdef ALLOW_STREAMICE

#ifdef STREAMICE_CONSTRUCT_MATRIX

C the linear action of the matrix on (u,v) with triangular finite elements
C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
C but this may change pursuant to conversations with others
C
C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
C     in order to make less frequent halo updates
C isym = 1 if grid is symmetric, 0 o.w.

C the linear action of the matrix on (u,v) with triangular finite elements
C Phi has the form
C Phi (i,j,k,q) - applies to cell i,j

C      3 - 4
C      |   |
C      1 - 2

C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q 
C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear

C     !LOCAL VARIABLES:
C     == Local variables == 
      INTEGER iq, jq, inodx, inody, i, j, bi, bj, ilqx, ilqy, m_i, n
      INTEGER jlqx, jlqy, jnodx,jnody, m_j, col_y, col_x, cg_halo, k
      INTEGER colx_rev, coly_rev
      _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy, tmpval
      _RL phival(2,2)

!       do i=1,3
!        do j=0,2
!         col_index_a = i + j*3
!        enddo
!       enddo

      cg_halo = min(OLx-1,OLy-1)

      DO j = 1-cg_halo, sNy+cg_halo
       DO i = 1-cg_halo, sNx+cg_halo
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
cc          DO k=1,4
           DO col_x=-1,1
            DO col_y=-1,1
             streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
            ENDDO
           ENDDO
cc          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

c$openad xxx simple loop
      DO j = 1-cg_halo, sNy+cg_halo
       DO i = 1-cg_halo, sNx+cg_halo
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
           DO iq=1,2 
            DO jq = 1,2
            
             n = 2*(jq-1)+iq
            
             DO inodx = 1,2 
              DO inody = 1,2

               if (STREAMICE_umask(i-1+inodx,j-1+inody,bi,bj)
     &          .eq.1.0 .or. 
     &             streamice_vmask(i-1+inodx,j-1+inody,bi,bj).eq.1.0)
     &          then            
             
                m_i = 2*(inody-1)+inodx
                ilqx = 1 
                ilqy = 1
            
                if (inodx.eq.iq) ilqx = 2
                if (inody.eq.jq) ilqy = 2  
                phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)

                DO jnodx = 1,2
                 DO jnody = 1,2 
                  if (STREAMICE_umask(i-1+jnodx,j-1+jnody,bi,bj)
     &             .eq.1.0 .or.
     &             STREAMICE_vmask(i-1+jnodx,j-1+jnody,bi,bj).eq.1.0) 
     &             then            

                   m_j = 2*(jnody-1)+jnodx
                   ilqx = 1 
                   ilqy = 1
                   if (jnodx.eq.iq) ilqx = 2
                   if (jnody.eq.jq) ilqy = 2  

!                    col_j = col_index_a (
!      &              jnodx+mod(inodx,2),
!      &              jnody+mod(inody,2) )

                   col_x = mod(inodx,2)+jnodx-2
                   colx_rev = mod(jnodx,2)+inodx-2
                   col_y = mod(inody,2)+jnody-2
                   coly_rev = mod(jnody,2)+inody-2
c


                   IF ( (inodx.eq.jnodx .and. inody.eq.jnody) .or.
     &                  (inodx.eq.1 .and. inody.eq.1) .or.
     &                  (jnody.eq.2 .and. inody.eq.1) .or.
     &                  (jnody.eq.2 .and. jnodx.eq.2)) THEN



                   ux = DPhi (i,j,bi,bj,m_j,n,1)
                   uy = DPhi (i,j,bi,bj,m_j,n,2)
                   vx = 0  
                   vy = 0
                   uq = Xquad(ilqx) * Xquad(ilqy) 
                   vq = 0

                   exx = ux + k1AtC_str(i,j,bi,bj)*vq
                   eyy = vy + k2AtC_str(i,j,bi,bj)*uq
                   exy = .5*(uy+vx) + 
     &              k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
                        
                    tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * (
     &              DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) + 
     &              DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))

                   streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval

                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF

!!!

                    tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * (
     &              DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) + 
     &              DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))

                   streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
    
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN 
                    streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF

!!!

                    tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
     &             (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
     &              exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)    

                   streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
     
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF

!!!

                   tmpval = .25 * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
     &             (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
     &              eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)

                   streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval

                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


!!!

                    tmpval = .25*phival(inodx,inody) * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              tau_beta_eff_streamice (i,j,bi,bj) * uq

                   streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A1
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A1
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


!!!
                    tmpval = .25*phival(inodx,inody) * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              tau_beta_eff_streamice (i,j,bi,bj) * vq

                   streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A3
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                  
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN 
                    streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A2
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


              
!!!

                   vx = DPhi (i,j,bi,bj,m_j,n,1)
                   vy = DPhi (i,j,bi,bj,m_j,n,2)
                   ux = 0 
                   uy = 0
                   vq = Xquad(ilqx) * Xquad(ilqy) 
                   uq = 0

                   exx = ux + k1AtC_str(i,j,bi,bj)*vq
                   eyy = vy + k2AtC_str(i,j,bi,bj)*uq
                   exy = .5*(uy+vx) + 
     &              k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq

                    tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * (
     &              DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) + 
     &              DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))

                   streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                    tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * (
     &              DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) + 
     &              DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))

                   streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                   tmpval = .25 *
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
     &             (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
     &              exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)   

                   streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                   tmpval = .25 * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
     &             (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
     &              eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)

                   streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                    tmpval = .25*phival(inodx,inody) * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              tau_beta_eff_streamice (i,j,bi,bj) * uq

                   streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A2
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A3
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                    tmpval = .25*phival(inodx,inody) * 
     &              grid_jacq_streamice(i,j,bi,bj,n) * 
     &              tau_beta_eff_streamice (i,j,bi,bj) * vq

                   streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)=
     &             streamice_cg_A4
     &                 (i-1+inodx,j-1+inody,bi,bj,mod(inodx,2)+jnodx-2,
     &                  mod(inody,2)+jnody-2)+tmpval
                   
                   IF (.not. (inodx.eq.jnodx .and. inody.eq.jnody)) THEN
                    streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)=
     &              streamice_cg_A4
     &               (i-1+jnodx,j-1+jnody,bi,bj,mod(jnodx,2)+inodx-2,
     &                mod(jnody,2)+inody-2)+
     &               tmpval
                   ENDIF


                   endif
                  endif 
                 enddo
                enddo
               endif
              enddo
             enddo
            enddo
           enddo
          endif
         enddo
        enddo
       enddo
      enddo



#endif
#endif
      RETURN
      END SUBROUTINE
! END MAKE_A

      SUBROUTINE STREAMICE_CG_ADIAG( myThid, 
     O      uret, 
     O      vret)
    
C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"            

C     !INPUT/OUTPUT ARGUMENTS
C     uret, vret - result of matrix operating on u, v
C     is, ie, js, je - starting and ending cells
      INTEGER myThid
      _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)


#ifdef ALLOW_STREAMICE

C the linear action of the matrix on (u,v) with triangular finite elements
C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
C but this may change pursuant to conversations with others
C
C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
C     in order to make less frequent halo updates
C isym = 1 if grid is symmetric, 0 o.w.

C the linear action of the matrix on (u,v) with triangular finite elements
C Phi has the form
C Phi (i,j,k,q) - applies to cell i,j

C      3 - 4
C      |   |
C      1 - 2

C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q 
C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear

C     !LOCAL VARIABLES:
C     == Local variables == 
      INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
      _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
      _RL Ucell (2,2)
      _RL Vcell (2,2)
      _RL Hcell (2,2)
      _RL phival(2,2)

      uret(1,1,1,1) = uret(1,1,1,1)
      vret(1,1,1,1) = vret(1,1,1,1)

      DO j = 0, sNy+1
       DO i = 0, sNx+1
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
           DO iq=1,2 
            DO jq = 1,2
            
            n = 2*(jq-1)+iq
            
            DO inode = 1,2 
             DO jnode = 1,2
             
             m = 2*(jnode-1)+inode

             if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0 .or.
     &           STREAMICE_vmask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) 
     &           then

              ilq = 1 
              jlq = 1
            
              if (inode.eq.iq) ilq = 2
              if (jnode.eq.jq) jlq = 2  
              phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)

              ux = DPhi (i,j,bi,bj,m,n,1)
              uy = DPhi (i,j,bi,bj,m,n,2)
              vx = 0  
              vy = 0
              uq = Xquad(ilq) * Xquad(jlq) 
              vq = 0

              exx = ux + k1AtC_str(i,j,bi,bj)*vq
              eyy = vy + k2AtC_str(i,j,bi,bj)*uq
              exy = .5*(uy+vx) + 
     &         k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq

              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) + 
     &          DPhi(i,j,bi,bj,m,n,2)*(2*exy))

              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
     &          4*0.5*k1AtC_str(i,j,bi,bj)*exy)    

              
              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * uq
              

              vx = DPhi (i,j,bi,bj,m,n,1)
              vy = DPhi (i,j,bi,bj,m,n,2)
              ux = 0 
              uy = 0
              vq = Xquad(ilq) * Xquad(jlq) 
              uq = 0

              exx = ux + k1AtC_str(i,j,bi,bj)*vq
              eyy = vy + k2AtC_str(i,j,bi,bj)*uq
              exy = .5*(uy+vx) + 
     &         k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq

              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) + 
     &          DPhi(i,j,bi,bj,m,n,1)*(2*exy))
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
     &          4*0.5*k2AtC_str(i,j,bi,bj)*exy)

              
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * vq
              
             endif

            enddo 
            enddo
           enddo 
           enddo
          endif
         enddo
        enddo
       enddo
      enddo

#endif
      RETURN
      END SUBROUTINE



      SUBROUTINE STREAMICE_CG_BOUND_VALS( myThid, 
     O    uret, 
     O    vret)
C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "STREAMICE.h" 
#include "STREAMICE_CG.h"

C     !INPUT/OUTPUT ARGUMENTS
C     uret, vret - result of matrix operating on u, v
C     is, ie, js, je - starting and ending cells
      INTEGER myThid
      _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
    
#ifdef ALLOW_STREAMICE

C the linear action of the matrix on (u,v) with triangular finite elements
C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
C but this may change pursuant to conversations with others
C
C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
C     in order to make less frequent halo updates
C isym = 1 if grid is symmetric, 0 o.w.

C the linear action of the matrix on (u,v) with triangular finite elements
C Phi has the form
C Phi (i,j,k,q) - applies to cell i,j

C      3 - 4
C      |   |
C      1 - 2

C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q 
C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear

C     !LOCAL VARIABLES:
C     == Local variables == 
      INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
      _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
      _RL Ucell (2,2)
      _RL Vcell (2,2)
      _RL Hcell (2,2)
      _RL phival(2,2)

      uret(1,1,1,1) = uret(1,1,1,1)
      vret(1,1,1,1) = vret(1,1,1,1)

      DO j = 0, sNy+1
       DO i = 0, sNx+1
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          IF ((STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) .AND.
     &     ((STREAMICE_umask(i,j,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_umask(i,j+1,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_umask(i+1,j,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_umask(i+1,j+1,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_vmask(i,j,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_vmask(i,j+1,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_vmask(i+1,j,bi,bj).eq.3.0) .OR.
     &      (STREAMICE_vmask(i+1,j+1,bi,bj).eq.3.0))) THEN
           
           DO iq=1,2 
            DO jq = 1,2

            n = 2*(jq-1)+iq

            uq = u_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+ 
     &       u_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+ 
     &       u_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+ 
     &       u_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
            vq = v_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+ 
     &       v_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+ 
     &       v_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+ 
     &       v_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
            ux = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) + 
     &       u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) + 
     &       u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) + 
     &       u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
            uy = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) +
     &       u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) + 
     &       u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) + 
     &       u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
            vx = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) + 
     &       v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) + 
     &       v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) + 
     &       v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
            vy = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) + 
     &       v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) + 
     &       v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) + 
     &       v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
            exx = ux + k1AtC_str(i,j,bi,bj)*vq
            eyy = vy + k2AtC_str(i,j,bi,bj)*uq
            exy = .5*(uy+vx) + 
     &       k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq


            do inode = 1,2 
             do jnode = 1,2

             m = 2*(jnode-1)+inode
             ilq = 1 
             jlq = 1
             if (inode.eq.iq) ilq = 2
             if (jnode.eq.jq) jlq = 2  
             phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)

             if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then          

              
              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) + 
     &          DPhi(i,j,bi,bj,m,n,2)*(2*exy))


              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
     &          4*0.5*k1AtC_str(i,j,bi,bj)*exy)    
             

!               if (STREAMICE_float_cond(i,j,bi,bj) .eq. 1) then
              uret(i-1+inode,j-1+jnode,bi,bj) = 
     &         uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * uq


!               endif
             endif
             if (STREAMICE_vmask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then 
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * (
     &          DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) + 
     &          DPhi(i,j,bi,bj,m,n,1)*(2*exy))
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         grid_jacq_streamice(i,j,bi,bj,n) * 
     &         visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
     &         (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
     &          4*0.5*k2AtC_str(i,j,bi,bj)*exy)
              vret(i-1+inode,j-1+jnode,bi,bj) = 
     &         vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
     &         phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) * 
     &         tau_beta_eff_streamice (i,j,bi,bj) * vq
             endif
            enddo 
            enddo
           enddo 
           enddo
          endif
         enddo
        enddo
       enddo
      enddo

#endif
      RETURN
      END SUBROUTINE