/[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.7 - (show annotations) (download)
Thu Mar 7 15:23:19 2013 UTC (12 years, 4 months ago) by dgoldberg
Branch: MAIN
Changes since 1.6: +6 -2 lines
bug fixes, GL smoothing, changes for controlling bathym with const surf elev

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

  ViewVC Help
Powered by ViewVC 1.1.22