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

Contents 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 - (show 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 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_functions.F,v 1.2 2012/05/02 02:36:01 heimbach Exp $
2 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 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
67 _RL Ucell (2,2)
68 _RL Vcell (2,2)
69 _RL Hcell (2,2)
70 _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
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
80 IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
81 DO iq = 1,2
82 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 jlq = 1
121 if (inode.eq.iq) ilq = 2
122 if (jnode.eq.jq) jlq = 2
123 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
124
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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
144 & (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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
150 & (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 & phival(inode,jnode) *
163 & grid_jacq_streamice(i,j,bi,bj,n) *
164 & 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 & phival(inode,jnode) *
168 & grid_jacq_streamice(i,j,bi,bj,n) *
169 & tau_beta_eff_streamice (i,j,bi,bj) * vq
170
171 endif
172 enddo
173 enddo
174
175 enddo
176 enddo
177 c-- STREAMICE_hmask
178 endif
179
180 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 #ifdef STREAMICE_CONSTRUCT_MATRIX
214
215 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 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
240 _RL phival(2,2)
241
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 phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)
293
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 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
354 & (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 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
364 & (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 & .25*phival(inodx,inody) *
372 & grid_jacq_streamice(i,j,bi,bj,n) *
373 & 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 & .25*phival(inodx,inody) *
380 & grid_jacq_streamice(i,j,bi,bj,n) *
381 & 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 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
424 & (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 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
434 & (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 & .25*phival(inodx,inody) *
442 & grid_jacq_streamice(i,j,bi,bj,n) *
443 & 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 & .25*phival(inodx,inody) *
450 & grid_jacq_streamice(i,j,bi,bj,n) *
451 & 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 #endif
469 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 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
526 _RL Ucell (2,2)
527 _RL Vcell (2,2)
528 _RL Hcell (2,2)
529 _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
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
549 if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
550
551 ilq = 1
552 jlq = 1
553
554 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
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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
581 & (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 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
588 & 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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
613 & (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 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
620 & tau_beta_eff_streamice (i,j,bi,bj) * vq
621
622 endif
623
624 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 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
692 _RL Ucell (2,2)
693 _RL Vcell (2,2)
694 _RL Hcell (2,2)
695 _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
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 jlq = 1
750 if (inode.eq.iq) ilq = 2
751 if (jnode.eq.jq) jlq = 2
752 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
753
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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
773 & (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 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
779 & (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 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
786 & 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 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
790 & 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