pkg/streamice/streamice_init_phi.F

5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP 0
                0005 C !ROUTINE: STREAMICE_INIT_FIXED
                0006
                0007 C !INTERFACE:
                0008       SUBROUTINE STREAMICE_INIT_PHI( myThid )
                0009
                0010 C     !DESCRIPTION:
                0011 C     Initialize STREAMICE nodal basis gradients for FEM solver
                0012
                0013 C     !USES:
                0014       IMPLICIT NONE
                0015 #include "EEPARAMS.h"
                0016 #include "SIZE.h"
                0017 #include "PARAMS.h"
                0018 #include "STREAMICE.h"
                0019 #include "STREAMICE_CG.h"
                0020 #include "GRID.h"
                0021
                0022 C     myThid ::  my Thread Id number
                0023       INTEGER myThid
                0024 CEOP
                0025
                0026 C     !LOCAL VARIABLES:
                0027 C     === Local variables ===
96b006450c dngo*0028       INTEGER bi, bj, i, j, xnode, ynode, xq, yq, m, n, kx, ky
5ca83cd8f7 Dani*0029       REAL gradx(2), grady(2)  ! gradients at quadrature points
                0030
887f427c62 Jean*0031 C     here the terms used to calculate matrix terms in the
5ca83cd8f7 Dani*0032 C     velocity solve are initialized
                0033 C
                0034 C     this is a quasi-finite element method; the gradient
                0035 C     of the basis functions are approximated based on knowledge
                0036 C     of the grid
                0037 C
887f427c62 Jean*0038 C     Dphi (i,j,bi,bj,m,n,p):
                0039 C       gradient (in p-direction) of nodal basis function in
                0040 C       cell (i,j) on thread (bi,bj) which is centered on node m,
5ca83cd8f7 Dani*0041 C       at quadrature point n
                0042 C
                0043 C    %  3 - 4
                0044 C    %  |   |
                0045 C    %  1 - 2
                0046 C
887f427c62 Jean*0047 C     NOTE 2x2 quadrature is hardcoded - might make it specifiable through CPP
5ca83cd8f7 Dani*0048 C
                0049 C     this will not be updated in overlap cells - so we extend it as far as we can
887f427c62 Jean*0050
5ca83cd8f7 Dani*0051       DO bj = myByLo(myThid), myByHi(myThid)
                0052        DO bi = myBxLo(myThid), myBxHi(myThid)
96b006450c dngo*0053         DO j=1-OLy,sNy+OLy-1
                0054          DO i=1-OLx,sNx+OLx-1
5ca83cd8f7 Dani*0055
                0056           DO xq = 1,2
887f427c62 Jean*0057            gradx(xq) = Xquad(3-xq) * recip_dxG (i,j,bi,bj) +
5ca83cd8f7 Dani*0058      &                 Xquad(xq) * recip_dxG (i+1,j,bi,bj)
887f427c62 Jean*0059            grady(xq) = Xquad(3-xq) * recip_dyG (i,j,bi,bj) +
5ca83cd8f7 Dani*0060      &                 Xquad(xq) * recip_dyG (i,j+1,bi,bj)
                0061           ENDDO
                0062
887f427c62 Jean*0063           DO n = 1,4
                0064
5ca83cd8f7 Dani*0065            xq = 2 - mod(n,2)
                0066            yq = floor ((n+1)/2.0)
887f427c62 Jean*0067
5ca83cd8f7 Dani*0068            DO m = 1,4
                0069
                0070             xnode = 2 - mod(m,2)
                0071             ynode = floor ((m+1)/2.0)
                0072
                0073             kx = 1 ; ky = 1
                0074             if (xq.eq.xnode) kx = 2
                0075             if (yq.eq.ynode) ky = 2
                0076
887f427c62 Jean*0077             Dphi (i,j,bi,bj,m,n,1) =
5ca83cd8f7 Dani*0078      &       (2*xnode-3) * Xquad(ky) * gradx(yq)
887f427c62 Jean*0079             Dphi (i,j,bi,bj,m,n,2) =
5ca83cd8f7 Dani*0080      &       (2*ynode-3) * Xquad(kx) * grady(xq)
887f427c62 Jean*0081
5ca83cd8f7 Dani*0082            ENDDO
                0083
887f427c62 Jean*0084            grid_jacq_streamice (i,j,bi,bj,n) =
                0085      &      (Xquad(3-xq)*dyG(i,j,bi,bj) + Xquad(xq)*dyG(i+1,j,bi,bj)) *
                0086      &      (Xquad(3-yq)*dxG(i,j,bi,bj) + Xquad(yq)*dxG(i,j+1,bi,bj))
5ca83cd8f7 Dani*0087
                0088           ENDDO
                0089          ENDDO
                0090         ENDDO
                0091        ENDDO
                0092       ENDDO
                0093
                0094       RETURN
                0095       END