1 |
C $Header: $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "STREAMICE_OPTIONS.h" |
5 |
|
6 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
7 |
CBOP 0 |
8 |
C !ROUTINE: STREAMICE_INIT_FIXED |
9 |
|
10 |
C !INTERFACE: |
11 |
SUBROUTINE STREAMICE_INIT_PHI( myThid ) |
12 |
|
13 |
C !DESCRIPTION: |
14 |
C Initialize STREAMICE nodal basis gradients for FEM solver |
15 |
|
16 |
C !USES: |
17 |
IMPLICIT NONE |
18 |
#include "EEPARAMS.h" |
19 |
#include "SIZE.h" |
20 |
#include "PARAMS.h" |
21 |
#include "STREAMICE.h" |
22 |
#include "STREAMICE_CG.h" |
23 |
#include "GRID.h" |
24 |
|
25 |
C myThid :: my Thread Id number |
26 |
INTEGER myThid |
27 |
CEOP |
28 |
|
29 |
C !LOCAL VARIABLES: |
30 |
C === Local variables === |
31 |
INTEGER bi, bj, i, j, xnode, ynode, xq, yq, m, n, p, kx, ky |
32 |
REAL gradx(2), grady(2) ! gradients at quadrature points |
33 |
|
34 |
C here the terms used to calculate matrix terms in the |
35 |
C velocity solve are initialized |
36 |
C |
37 |
C this is a quasi-finite element method; the gradient |
38 |
C of the basis functions are approximated based on knowledge |
39 |
C of the grid |
40 |
C |
41 |
C Dphi (i,j,bi,bj,m,n,p): |
42 |
C gradient (in p-direction) of nodal basis function in |
43 |
C cell (i,j) on thread (bi,bj) which is centered on node m, |
44 |
C at quadrature point n |
45 |
C |
46 |
C % 3 - 4 |
47 |
C % | | |
48 |
C % 1 - 2 |
49 |
C |
50 |
C NOTE 2x2 quadrature is hardcoded - might make it specifiable through CPP |
51 |
C |
52 |
C this will not be updated in overlap cells - so we extend it as far as we can |
53 |
|
54 |
DO bj = myByLo(myThid), myByHi(myThid) |
55 |
DO bi = myBxLo(myThid), myBxHi(myThid) |
56 |
DO j=1-Oly,sNy+Oly-1 |
57 |
DO i=1-Olx,sNx+Olx-1 |
58 |
|
59 |
DO xq = 1,2 |
60 |
gradx(xq) = Xquad(3-xq) * recip_dxG (i,j,bi,bj) + |
61 |
& Xquad(xq) * recip_dxG (i+1,j,bi,bj) |
62 |
grady(xq) = Xquad(3-xq) * recip_dyG (i,j,bi,bj) + |
63 |
& Xquad(xq) * recip_dyG (i,j+1,bi,bj) |
64 |
ENDDO |
65 |
|
66 |
DO n = 1,4 |
67 |
|
68 |
xq = 2 - mod(n,2) |
69 |
yq = floor ((n+1)/2.0) |
70 |
|
71 |
DO m = 1,4 |
72 |
|
73 |
xnode = 2 - mod(m,2) |
74 |
ynode = floor ((m+1)/2.0) |
75 |
|
76 |
kx = 1 ; ky = 1 |
77 |
if (xq.eq.xnode) kx = 2 |
78 |
if (yq.eq.ynode) ky = 2 |
79 |
|
80 |
|
81 |
Dphi (i,j,bi,bj,m,n,1) = |
82 |
& (2*xnode-3) * Xquad(ky) * gradx(yq) |
83 |
Dphi (i,j,bi,bj,m,n,2) = |
84 |
& (2*ynode-3) * Xquad(kx) * grady(xq) |
85 |
|
86 |
ENDDO |
87 |
|
88 |
grid_jacq_streamice (i,j,bi,bj,n) = |
89 |
& (Xquad(3-xq)*dyG(i,j,bi,bj) + Xquad(xq)*dyG(i+1,j,bi,bj)) * |
90 |
& (Xquad(3-yq)*dxG(i,j,bi,bj) + Xquad(yq)*dxG(i,j+1,bi,bj)) |
91 |
|
92 |
ENDDO |
93 |
ENDDO |
94 |
ENDDO |
95 |
ENDDO |
96 |
ENDDO |
97 |
|
98 |
RETURN |
99 |
END |