/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_cg_functions.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_cg_functions.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Mon May 14 16:53:09 2012 UTC (13 years, 2 months ago) by dgoldberg
Branch: MAIN
Changes since 1.2: +4 -1 lines
make construct matrix option a CPP option

1 dgoldberg 1.3 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_functions.F,v 1.2 2012/05/02 02:36:01 heimbach Exp $
2 heimbach 1.1 C $Name: $
3    
4     #include "STREAMICE_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9     SUBROUTINE STREAMICE_CG_ACTION( myThid,
10     O uret,
11     O vret,
12     I u,
13     I v,
14     I is, ie, js, je )
15     C /============================================================\
16     C | SUBROUTINE |
17     C | o |
18     C |============================================================|
19     C | |
20     C \============================================================/
21     IMPLICIT NONE
22    
23     C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "STREAMICE.h"
29     #include "STREAMICE_CG.h"
30    
31     C !INPUT/OUTPUT ARGUMENTS
32     C uret, vret - result of matrix operating on u, v
33     C is, ie, js, je - starting and ending cells
34     INTEGER myThid
35     _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
36     _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
37     _RL u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
38     _RL v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
39     INTEGER is, ie, js, je
40    
41     #ifdef ALLOW_STREAMICE
42    
43     C the linear action of the matrix on (u,v) with triangular finite elements
44     C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
45     C but this may change pursuant to conversations with others
46     C
47     C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
48     C in order to make less frequent halo updates
49     C isym = 1 if grid is symmetric, 0 o.w.
50    
51     C the linear action of the matrix on (u,v) with triangular finite elements
52     C Phi has the form
53     C Phi (i,j,k,q) - applies to cell i,j
54    
55     C 3 - 4
56     C | |
57     C 1 - 2
58    
59     C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
60     C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
61     C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
62    
63     C !LOCAL VARIABLES:
64     C == Local variables ==
65     INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
66 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
67 heimbach 1.1 _RL Ucell (2,2)
68     _RL Vcell (2,2)
69     _RL Hcell (2,2)
70 heimbach 1.2 _RL phival(2,2)
71    
72     uret(1,1,1,1) = uret(1,1,1,1)
73     vret(1,1,1,1) = vret(1,1,1,1)
74 heimbach 1.1
75     DO j = js, je
76     DO i = is, ie
77     DO bj = myByLo(myThid), myByHi(myThid)
78     DO bi = myBxLo(myThid), myBxHi(myThid)
79 heimbach 1.2
80 heimbach 1.1 IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
81 heimbach 1.2 DO iq = 1,2
82 heimbach 1.1 DO jq = 1,2
83    
84     n = 2*(jq-1)+iq
85    
86     uq = u(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) +
87     & u(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) +
88     & u(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) +
89     & u(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
90     vq = v(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) +
91     & v(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) +
92     & v(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) +
93     & v(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
94     ux = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
95     & u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
96     & u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
97     & u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
98     uy = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) +
99     & u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
100     & u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
101     & u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
102     vx = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
103     & v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
104     & v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
105     & v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
106     vy = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) +
107     & v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
108     & v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
109     & v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
110     exx = ux + k1AtC_str(i,j,bi,bj)*vq
111     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
112     exy = .5*(uy+vx) +
113     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
114    
115     do inode = 1,2
116     do jnode = 1,2
117    
118     m = 2*(jnode-1)+inode
119     ilq = 1
120 heimbach 1.2 jlq = 1
121 heimbach 1.1 if (inode.eq.iq) ilq = 2
122 heimbach 1.2 if (jnode.eq.jq) jlq = 2
123     phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
124 heimbach 1.1
125     if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
126    
127     uret(i-1+inode,j-1+jnode,bi,bj) =
128     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
129     & grid_jacq_streamice(i,j,bi,bj,n) *
130     & visc_streamice(i,j,bi,bj) * (
131     & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
132     & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
133     vret(i-1+inode,j-1+jnode,bi,bj) =
134     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
135     & grid_jacq_streamice(i,j,bi,bj,n) *
136     & visc_streamice(i,j,bi,bj) * (
137     & DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) +
138     & DPhi(i,j,bi,bj,m,n,1)*(2*exy))
139    
140     uret(i-1+inode,j-1+jnode,bi,bj) =
141     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
142     & grid_jacq_streamice(i,j,bi,bj,n) *
143 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
144 heimbach 1.1 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
145     & 4*0.5*k1AtC_str(i,j,bi,bj)*exy)
146     vret(i-1+inode,j-1+jnode,bi,bj) =
147     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
148     & grid_jacq_streamice(i,j,bi,bj,n) *
149 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
150 heimbach 1.1 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
151     & 4*0.5*k2AtC_str(i,j,bi,bj)*exy)
152    
153     ! IF (bi.eq.2.and.bj.eq.2.and.i.eq.15.and.
154     ! & (exx.ne.0.0 .or. eyy.ne.0.0 .or. exy.ne.0.0)) THEN
155     ! PRINT *, "CG_FUNCTION", j, v(i,j,bi,bj),v(i+1,j,bi,bj),
156     ! & v(i,j+1,bi,bj),v(i+1,j+1,bi,bj)
157     ! ENDIF
158    
159    
160     uret(i-1+inode,j-1+jnode,bi,bj) =
161     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
162 heimbach 1.2 & phival(inode,jnode) *
163     & grid_jacq_streamice(i,j,bi,bj,n) *
164 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * uq
165     vret(i-1+inode,j-1+jnode,bi,bj) =
166     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
167 heimbach 1.2 & phival(inode,jnode) *
168     & grid_jacq_streamice(i,j,bi,bj,n) *
169 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
170    
171     endif
172     enddo
173     enddo
174 heimbach 1.2
175 heimbach 1.1 enddo
176     enddo
177 heimbach 1.2 c-- STREAMICE_hmask
178 heimbach 1.1 endif
179 heimbach 1.2
180 heimbach 1.1 enddo
181     enddo
182     enddo
183     enddo
184    
185     #endif
186     RETURN
187     END SUBROUTINE
188    
189     SUBROUTINE STREAMICE_CG_MAKE_A( myThid )
190     C /============================================================\
191     C | SUBROUTINE |
192     C | o |
193     C |============================================================|
194     C | |
195     C \============================================================/
196     IMPLICIT NONE
197    
198     C === Global variables ===
199     #include "SIZE.h"
200     #include "EEPARAMS.h"
201     #include "PARAMS.h"
202     #include "GRID.h"
203     #include "STREAMICE.h"
204     #include "STREAMICE_CG.h"
205    
206     C !INPUT/OUTPUT ARGUMENTS
207     C uret, vret - result of matrix operating on u, v
208     C is, ie, js, je - starting and ending cells
209     INTEGER myThid
210    
211     #ifdef ALLOW_STREAMICE
212    
213 dgoldberg 1.3 #ifdef STREAMICE_CONSTRUCT_MATRIX
214    
215 heimbach 1.1 C the linear action of the matrix on (u,v) with triangular finite elements
216     C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
217     C but this may change pursuant to conversations with others
218     C
219     C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
220     C in order to make less frequent halo updates
221     C isym = 1 if grid is symmetric, 0 o.w.
222    
223     C the linear action of the matrix on (u,v) with triangular finite elements
224     C Phi has the form
225     C Phi (i,j,k,q) - applies to cell i,j
226    
227     C 3 - 4
228     C | |
229     C 1 - 2
230    
231     C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
232     C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
233     C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
234    
235     C !LOCAL VARIABLES:
236     C == Local variables ==
237     INTEGER iq, jq, inodx, inody, i, j, bi, bj, ilqx, ilqy, m_i, n
238     INTEGER jlqx, jlqy, jnodx,jnody, m_j, col_y, col_x, cg_halo, k
239 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
240     _RL phival(2,2)
241 heimbach 1.1
242     ! do i=1,3
243     ! do j=0,2
244     ! col_index_a = i + j*3
245     ! enddo
246     ! enddo
247    
248     cg_halo = min(OLx-1,OLy-1)
249    
250     DO j = 1-cg_halo, sNy+cg_halo
251     DO i = 1-cg_halo, sNx+cg_halo
252     DO bj = myByLo(myThid), myByHi(myThid)
253     DO bi = myBxLo(myThid), myBxHi(myThid)
254     cc DO k=1,4
255     DO col_x=-1,1
256     DO col_y=-1,1
257     streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
258     streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
259     streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
260     streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
261     ENDDO
262     ENDDO
263     cc ENDDO
264     ENDDO
265     ENDDO
266     ENDDO
267     ENDDO
268    
269     DO j = 1-cg_halo, sNy+cg_halo
270     DO i = 1-cg_halo, sNx+cg_halo
271     DO bj = myByLo(myThid), myByHi(myThid)
272     DO bi = myBxLo(myThid), myBxHi(myThid)
273     IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
274     DO iq=1,2
275     DO jq = 1,2
276    
277     n = 2*(jq-1)+iq
278    
279     DO inodx = 1,2
280     DO inody = 1,2
281    
282     if (STREAMICE_umask(i-1+inodx,j-1+inody,bi,bj)
283     & .eq.1.0)
284     & then
285    
286     m_i = 2*(inody-1)+inodx
287     ilqx = 1
288     ilqy = 1
289    
290     if (inodx.eq.iq) ilqx = 2
291     if (inody.eq.jq) ilqy = 2
292 heimbach 1.2 phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)
293 heimbach 1.1
294     DO jnodx = 1,2
295     DO jnody = 1,2
296     if (STREAMICE_umask(i-1+jnodx,j-1+jnody,bi,bj)
297     & .eq.1.0)
298     & then
299    
300     m_j = 2*(jnody-1)+jnodx
301     ilqx = 1
302     ilqy = 1
303     if (jnodx.eq.iq) ilqx = 2
304     if (jnody.eq.jq) ilqy = 2
305    
306     ! col_j = col_index_a (
307     ! & jnodx+mod(inodx,2),
308     ! & jnody+mod(inody,2) )
309    
310     col_x = mod(inodx,2)+jnodx-2
311     col_y = mod(inody,2)+jnody-2
312    
313     c
314    
315     ux = DPhi (i,j,bi,bj,m_j,n,1)
316     uy = DPhi (i,j,bi,bj,m_j,n,2)
317     vx = 0
318     vy = 0
319     uq = Xquad(ilqx) * Xquad(ilqy)
320     vq = 0
321    
322     exx = ux + k1AtC_str(i,j,bi,bj)*vq
323     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
324     exy = .5*(uy+vx) +
325     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
326    
327     streamice_cg_A1
328     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
329     & streamice_cg_A1
330     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
331     & .25 *
332     & grid_jacq_streamice(i,j,bi,bj,n) *
333     & visc_streamice(i,j,bi,bj) * (
334     & DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) +
335     & DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))
336    
337     streamice_cg_A3
338     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
339     & streamice_cg_A3
340     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
341     & .25 *
342     & grid_jacq_streamice(i,j,bi,bj,n) *
343     & visc_streamice(i,j,bi,bj) * (
344     & DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) +
345     & DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))
346    
347     streamice_cg_A1
348     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
349     & streamice_cg_A1
350     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
351     & .25 *
352     & grid_jacq_streamice(i,j,bi,bj,n) *
353 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
354 heimbach 1.1 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
355     & exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)
356    
357     streamice_cg_A3
358     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
359     & streamice_cg_A3
360     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
361     & .25 *
362     & grid_jacq_streamice(i,j,bi,bj,n) *
363 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
364 heimbach 1.1 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
365     & eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)
366    
367     streamice_cg_A1
368     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
369     & streamice_cg_A1
370     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
371 heimbach 1.2 & .25*phival(inodx,inody) *
372     & grid_jacq_streamice(i,j,bi,bj,n) *
373 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * uq
374    
375     streamice_cg_A3
376     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
377     & streamice_cg_A3
378     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
379 heimbach 1.2 & .25*phival(inodx,inody) *
380     & grid_jacq_streamice(i,j,bi,bj,n) *
381 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
382    
383     c
384    
385     vx = DPhi (i,j,bi,bj,m_j,n,1)
386     vy = DPhi (i,j,bi,bj,m_j,n,2)
387     ux = 0
388     uy = 0
389     vq = Xquad(ilqx) * Xquad(ilqy)
390     uq = 0
391    
392     exx = ux + k1AtC_str(i,j,bi,bj)*vq
393     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
394     exy = .5*(uy+vx) +
395     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
396    
397     streamice_cg_A2
398     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
399     & streamice_cg_A2
400     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
401     & .25 *
402     & grid_jacq_streamice(i,j,bi,bj,n) *
403     & visc_streamice(i,j,bi,bj) * (
404     & DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) +
405     & DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))
406    
407     streamice_cg_A4
408     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
409     & streamice_cg_A4
410     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
411     & .25 *
412     & grid_jacq_streamice(i,j,bi,bj,n) *
413     & visc_streamice(i,j,bi,bj) * (
414     & DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) +
415     & DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))
416    
417     streamice_cg_A2
418     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
419     & streamice_cg_A2
420     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
421     & .25 *
422     & grid_jacq_streamice(i,j,bi,bj,n) *
423 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
424 heimbach 1.1 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
425     & exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)
426    
427     streamice_cg_A4
428     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
429     & streamice_cg_A4
430     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
431     & .25 *
432     & grid_jacq_streamice(i,j,bi,bj,n) *
433 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
434 heimbach 1.1 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
435     & eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)
436    
437     streamice_cg_A2
438     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
439     & streamice_cg_A2
440     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
441 heimbach 1.2 & .25*phival(inodx,inody) *
442     & grid_jacq_streamice(i,j,bi,bj,n) *
443 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * uq
444    
445     streamice_cg_A4
446     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
447     & streamice_cg_A4
448     & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
449 heimbach 1.2 & .25*phival(inodx,inody) *
450     & grid_jacq_streamice(i,j,bi,bj,n) *
451 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
452    
453     endif
454     enddo
455     enddo
456     endif
457     enddo
458     enddo
459     enddo
460     enddo
461     endif
462     enddo
463     enddo
464     enddo
465     enddo
466    
467     #endif
468 dgoldberg 1.3 #endif
469 heimbach 1.1 RETURN
470     END SUBROUTINE
471    
472     SUBROUTINE STREAMICE_CG_ADIAG( myThid,
473     O uret,
474     O vret)
475    
476     C /============================================================\
477     C | SUBROUTINE |
478     C | o |
479     C |============================================================|
480     C | |
481     C \============================================================/
482     IMPLICIT NONE
483    
484     C === Global variables ===
485     #include "SIZE.h"
486     #include "EEPARAMS.h"
487     #include "PARAMS.h"
488     #include "GRID.h"
489     #include "STREAMICE.h"
490     #include "STREAMICE_CG.h"
491    
492     C !INPUT/OUTPUT ARGUMENTS
493     C uret, vret - result of matrix operating on u, v
494     C is, ie, js, je - starting and ending cells
495     INTEGER myThid
496     _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
497     _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
498    
499    
500     #ifdef ALLOW_STREAMICE
501    
502     C the linear action of the matrix on (u,v) with triangular finite elements
503     C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
504     C but this may change pursuant to conversations with others
505     C
506     C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
507     C in order to make less frequent halo updates
508     C isym = 1 if grid is symmetric, 0 o.w.
509    
510     C the linear action of the matrix on (u,v) with triangular finite elements
511     C Phi has the form
512     C Phi (i,j,k,q) - applies to cell i,j
513    
514     C 3 - 4
515     C | |
516     C 1 - 2
517    
518     C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
519     C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
520     C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
521    
522     C !LOCAL VARIABLES:
523     C == Local variables ==
524     INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
525 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
526 heimbach 1.1 _RL Ucell (2,2)
527     _RL Vcell (2,2)
528     _RL Hcell (2,2)
529 heimbach 1.2 _RL phival(2,2)
530    
531     uret(1,1,1,1) = uret(1,1,1,1)
532     vret(1,1,1,1) = vret(1,1,1,1)
533 heimbach 1.1
534     DO j = 0, sNy+1
535     DO i = 0, sNx+1
536     DO bj = myByLo(myThid), myByHi(myThid)
537     DO bi = myBxLo(myThid), myBxHi(myThid)
538     IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
539     DO iq=1,2
540     DO jq = 1,2
541    
542     n = 2*(jq-1)+iq
543    
544     DO inode = 1,2
545     DO jnode = 1,2
546    
547     m = 2*(jnode-1)+inode
548 heimbach 1.2
549     if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
550    
551     ilq = 1
552     jlq = 1
553 heimbach 1.1
554 heimbach 1.2 if (inode.eq.iq) ilq = 2
555     if (jnode.eq.jq) jlq = 2
556     phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
557    
558     ux = DPhi (i,j,bi,bj,m,n,1)
559     uy = DPhi (i,j,bi,bj,m,n,2)
560     vx = 0
561     vy = 0
562     uq = Xquad(ilq) * Xquad(jlq)
563     vq = 0
564    
565     exx = ux + k1AtC_str(i,j,bi,bj)*vq
566     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
567     exy = .5*(uy+vx) +
568     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
569 heimbach 1.1
570     uret(i-1+inode,j-1+jnode,bi,bj) =
571     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
572     & grid_jacq_streamice(i,j,bi,bj,n) *
573     & visc_streamice(i,j,bi,bj) * (
574     & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
575     & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
576    
577     uret(i-1+inode,j-1+jnode,bi,bj) =
578     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
579     & grid_jacq_streamice(i,j,bi,bj,n) *
580 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
581 heimbach 1.1 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
582     & 4*0.5*k1AtC_str(i,j,bi,bj)*exy)
583    
584    
585     uret(i-1+inode,j-1+jnode,bi,bj) =
586     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
587 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
588 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * uq
589    
590    
591     vx = DPhi (i,j,bi,bj,m,n,1)
592     vy = DPhi (i,j,bi,bj,m,n,2)
593     ux = 0
594     uy = 0
595     vq = Xquad(ilq) * Xquad(jlq)
596     uq = 0
597    
598     exx = ux + k1AtC_str(i,j,bi,bj)*vq
599     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
600     exy = .5*(uy+vx) +
601     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
602    
603     vret(i-1+inode,j-1+jnode,bi,bj) =
604     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
605     & grid_jacq_streamice(i,j,bi,bj,n) *
606     & visc_streamice(i,j,bi,bj) * (
607     & DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) +
608     & DPhi(i,j,bi,bj,m,n,1)*(2*exy))
609     vret(i-1+inode,j-1+jnode,bi,bj) =
610     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
611     & grid_jacq_streamice(i,j,bi,bj,n) *
612 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
613 heimbach 1.1 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
614     & 4*0.5*k2AtC_str(i,j,bi,bj)*exy)
615    
616    
617     vret(i-1+inode,j-1+jnode,bi,bj) =
618     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
619 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
620 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
621    
622     endif
623 heimbach 1.2
624 heimbach 1.1 enddo
625     enddo
626     enddo
627     enddo
628     endif
629     enddo
630     enddo
631     enddo
632     enddo
633    
634     #endif
635     RETURN
636     END SUBROUTINE
637    
638    
639    
640     SUBROUTINE STREAMICE_CG_BOUND_VALS( myThid,
641     O uret,
642     O vret)
643     C /============================================================\
644     C | SUBROUTINE |
645     C | o |
646     C |============================================================|
647     C | |
648     C \============================================================/
649     IMPLICIT NONE
650    
651     C === Global variables ===
652     #include "SIZE.h"
653     #include "EEPARAMS.h"
654     #include "PARAMS.h"
655     #include "GRID.h"
656     #include "STREAMICE.h"
657     #include "STREAMICE_CG.h"
658    
659     C !INPUT/OUTPUT ARGUMENTS
660     C uret, vret - result of matrix operating on u, v
661     C is, ie, js, je - starting and ending cells
662     INTEGER myThid
663     _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
664     _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
665    
666     #ifdef ALLOW_STREAMICE
667    
668     C the linear action of the matrix on (u,v) with triangular finite elements
669     C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
670     C but this may change pursuant to conversations with others
671     C
672     C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
673     C in order to make less frequent halo updates
674     C isym = 1 if grid is symmetric, 0 o.w.
675    
676     C the linear action of the matrix on (u,v) with triangular finite elements
677     C Phi has the form
678     C Phi (i,j,k,q) - applies to cell i,j
679    
680     C 3 - 4
681     C | |
682     C 1 - 2
683    
684     C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
685     C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
686     C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
687    
688     C !LOCAL VARIABLES:
689     C == Local variables ==
690     INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
691 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
692 heimbach 1.1 _RL Ucell (2,2)
693     _RL Vcell (2,2)
694     _RL Hcell (2,2)
695 heimbach 1.2 _RL phival(2,2)
696    
697     uret(1,1,1,1) = uret(1,1,1,1)
698     vret(1,1,1,1) = vret(1,1,1,1)
699 heimbach 1.1
700     DO j = 0, sNy+1
701     DO i = 0, sNx+1
702     DO bj = myByLo(myThid), myByHi(myThid)
703     DO bi = myBxLo(myThid), myBxHi(myThid)
704     IF ((STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) .AND.
705     & ((STREAMICE_umask(i,j,bi,bj).eq.3.0) .OR.
706     & (STREAMICE_umask(i,j+1,bi,bj).eq.3.0) .OR.
707     & (STREAMICE_umask(i+1,j,bi,bj).eq.3.0) .OR.
708     & (STREAMICE_umask(i+1,j+1,bi,bj).eq.3.0))) THEN
709    
710     DO iq=1,2
711     DO jq = 1,2
712    
713     n = 2*(jq-1)+iq
714    
715     uq = u_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+
716     & u_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+
717     & u_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+
718     & u_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
719     vq = v_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+
720     & v_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+
721     & v_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+
722     & v_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
723     ux = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
724     & u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
725     & u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
726     & u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
727     uy = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
728     & u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
729     & u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
730     & u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
731     vx = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
732     & v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
733     & v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
734     & v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
735     vy = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
736     & v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
737     & v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
738     & v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
739     exx = ux + k1AtC_str(i,j,bi,bj)*vq
740     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
741     exy = .5*(uy+vx) +
742     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
743    
744     do inode = 1,2
745     do jnode = 1,2
746    
747     m = 2*(jnode-1)+inode
748     ilq = 1
749 heimbach 1.2 jlq = 1
750 heimbach 1.1 if (inode.eq.iq) ilq = 2
751     if (jnode.eq.jq) jlq = 2
752 heimbach 1.2 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
753 heimbach 1.1
754     if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
755    
756     uret(i-1+inode,j-1+jnode,bi,bj) =
757     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
758     & grid_jacq_streamice(i,j,bi,bj,n) *
759     & visc_streamice(i,j,bi,bj) * (
760     & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
761     & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
762     vret(i-1+inode,j-1+jnode,bi,bj) =
763     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
764     & grid_jacq_streamice(i,j,bi,bj,n) *
765     & visc_streamice(i,j,bi,bj) * (
766     & DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) +
767     & DPhi(i,j,bi,bj,m,n,1)*(2*exy))
768    
769     uret(i-1+inode,j-1+jnode,bi,bj) =
770     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
771     & grid_jacq_streamice(i,j,bi,bj,n) *
772 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
773 heimbach 1.1 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
774     & 4*0.5*k1AtC_str(i,j,bi,bj)*exy)
775     vret(i-1+inode,j-1+jnode,bi,bj) =
776     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
777     & grid_jacq_streamice(i,j,bi,bj,n) *
778 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
779 heimbach 1.1 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
780     & 4*0.5*k2AtC_str(i,j,bi,bj)*exy)
781    
782     ! if (STREAMICE_float_cond(i,j,bi,bj) .eq. 1) then
783     uret(i-1+inode,j-1+jnode,bi,bj) =
784     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
785 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
786 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * uq
787     vret(i-1+inode,j-1+jnode,bi,bj) =
788     & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
789 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
790 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
791     ! endif
792     endif
793     enddo
794     enddo
795     enddo
796     enddo
797     endif
798     enddo
799     enddo
800     enddo
801     enddo
802    
803     #endif
804     RETURN
805     END SUBROUTINE

  ViewVC Help
Powered by ViewVC 1.1.22