/[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.5 - (show annotations) (download)
Sun Dec 23 21:05:08 2012 UTC (12 years, 7 months ago) by dgoldberg
Branch: MAIN
Changes since 1.4: +47 -33 lines
fixing some bugs involving velocity masks

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

  ViewVC Help
Powered by ViewVC 1.1.22