/[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.2 - (show annotations) (download)
Wed May 2 02:36:01 2012 UTC (13 years, 2 months ago) by heimbach
Branch: MAIN
Changes since 1.1: +75 -50 lines
Various updates to streamice code

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng 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 C the linear action of the matrix on (u,v) with triangular finite elements
214 C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
215 C but this may change pursuant to conversations with others
216 C
217 C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
218 C in order to make less frequent halo updates
219 C isym = 1 if grid is symmetric, 0 o.w.
220
221 C the linear action of the matrix on (u,v) with triangular finite elements
222 C Phi has the form
223 C Phi (i,j,k,q) - applies to cell i,j
224
225 C 3 - 4
226 C | |
227 C 1 - 2
228
229 C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
230 C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
231 C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
232
233 C !LOCAL VARIABLES:
234 C == Local variables ==
235 INTEGER iq, jq, inodx, inody, i, j, bi, bj, ilqx, ilqy, m_i, n
236 INTEGER jlqx, jlqy, jnodx,jnody, m_j, col_y, col_x, cg_halo, k
237 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
238 _RL phival(2,2)
239
240 ! do i=1,3
241 ! do j=0,2
242 ! col_index_a = i + j*3
243 ! enddo
244 ! enddo
245
246 cg_halo = min(OLx-1,OLy-1)
247
248 DO j = 1-cg_halo, sNy+cg_halo
249 DO i = 1-cg_halo, sNx+cg_halo
250 DO bj = myByLo(myThid), myByHi(myThid)
251 DO bi = myBxLo(myThid), myBxHi(myThid)
252 cc DO k=1,4
253 DO col_x=-1,1
254 DO col_y=-1,1
255 streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
256 streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
257 streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
258 streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
259 ENDDO
260 ENDDO
261 cc ENDDO
262 ENDDO
263 ENDDO
264 ENDDO
265 ENDDO
266
267 DO j = 1-cg_halo, sNy+cg_halo
268 DO i = 1-cg_halo, sNx+cg_halo
269 DO bj = myByLo(myThid), myByHi(myThid)
270 DO bi = myBxLo(myThid), myBxHi(myThid)
271 IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
272 DO iq=1,2
273 DO jq = 1,2
274
275 n = 2*(jq-1)+iq
276
277 DO inodx = 1,2
278 DO inody = 1,2
279
280 if (STREAMICE_umask(i-1+inodx,j-1+inody,bi,bj)
281 & .eq.1.0)
282 & then
283
284 m_i = 2*(inody-1)+inodx
285 ilqx = 1
286 ilqy = 1
287
288 if (inodx.eq.iq) ilqx = 2
289 if (inody.eq.jq) ilqy = 2
290 phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)
291
292 DO jnodx = 1,2
293 DO jnody = 1,2
294 if (STREAMICE_umask(i-1+jnodx,j-1+jnody,bi,bj)
295 & .eq.1.0)
296 & then
297
298 m_j = 2*(jnody-1)+jnodx
299 ilqx = 1
300 ilqy = 1
301 if (jnodx.eq.iq) ilqx = 2
302 if (jnody.eq.jq) ilqy = 2
303
304 ! col_j = col_index_a (
305 ! & jnodx+mod(inodx,2),
306 ! & jnody+mod(inody,2) )
307
308 col_x = mod(inodx,2)+jnodx-2
309 col_y = mod(inody,2)+jnody-2
310
311 c
312
313 ux = DPhi (i,j,bi,bj,m_j,n,1)
314 uy = DPhi (i,j,bi,bj,m_j,n,2)
315 vx = 0
316 vy = 0
317 uq = Xquad(ilqx) * Xquad(ilqy)
318 vq = 0
319
320 exx = ux + k1AtC_str(i,j,bi,bj)*vq
321 eyy = vy + k2AtC_str(i,j,bi,bj)*uq
322 exy = .5*(uy+vx) +
323 & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
324
325 streamice_cg_A1
326 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
327 & streamice_cg_A1
328 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
329 & .25 *
330 & grid_jacq_streamice(i,j,bi,bj,n) *
331 & visc_streamice(i,j,bi,bj) * (
332 & DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) +
333 & DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))
334
335 streamice_cg_A3
336 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
337 & streamice_cg_A3
338 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
339 & .25 *
340 & grid_jacq_streamice(i,j,bi,bj,n) *
341 & visc_streamice(i,j,bi,bj) * (
342 & DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) +
343 & DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))
344
345 streamice_cg_A1
346 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
347 & streamice_cg_A1
348 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
349 & .25 *
350 & grid_jacq_streamice(i,j,bi,bj,n) *
351 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
352 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
353 & exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)
354
355 streamice_cg_A3
356 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
357 & streamice_cg_A3
358 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
359 & .25 *
360 & grid_jacq_streamice(i,j,bi,bj,n) *
361 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
362 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
363 & eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)
364
365 streamice_cg_A1
366 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
367 & streamice_cg_A1
368 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
369 & .25*phival(inodx,inody) *
370 & grid_jacq_streamice(i,j,bi,bj,n) *
371 & tau_beta_eff_streamice (i,j,bi,bj) * uq
372
373 streamice_cg_A3
374 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
375 & streamice_cg_A3
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) * vq
380
381 c
382
383 vx = DPhi (i,j,bi,bj,m_j,n,1)
384 vy = DPhi (i,j,bi,bj,m_j,n,2)
385 ux = 0
386 uy = 0
387 vq = Xquad(ilqx) * Xquad(ilqy)
388 uq = 0
389
390 exx = ux + k1AtC_str(i,j,bi,bj)*vq
391 eyy = vy + k2AtC_str(i,j,bi,bj)*uq
392 exy = .5*(uy+vx) +
393 & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
394
395 streamice_cg_A2
396 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
397 & streamice_cg_A2
398 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
399 & .25 *
400 & grid_jacq_streamice(i,j,bi,bj,n) *
401 & visc_streamice(i,j,bi,bj) * (
402 & DPhi(i,j,bi,bj,m_i,n,1)*(4*exx+2*eyy) +
403 & DPhi(i,j,bi,bj,m_i,n,2)*(2*exy))
404
405 streamice_cg_A4
406 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
407 & streamice_cg_A4
408 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
409 & .25 *
410 & grid_jacq_streamice(i,j,bi,bj,n) *
411 & visc_streamice(i,j,bi,bj) * (
412 & DPhi(i,j,bi,bj,m_i,n,2)*(4*eyy+2*exx) +
413 & DPhi(i,j,bi,bj,m_i,n,1)*(2*exy))
414
415 streamice_cg_A2
416 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
417 & streamice_cg_A2
418 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
419 & .25 *
420 & grid_jacq_streamice(i,j,bi,bj,n) *
421 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
422 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*
423 & exx+4*0.5*k1AtC_str(i,j,bi,bj)*exy)
424
425 streamice_cg_A4
426 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
427 & streamice_cg_A4
428 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
429 & .25 *
430 & grid_jacq_streamice(i,j,bi,bj,n) *
431 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
432 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*
433 & eyy+4*0.5*k2AtC_str(i,j,bi,bj)*exy)
434
435 streamice_cg_A2
436 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
437 & streamice_cg_A2
438 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)+
439 & .25*phival(inodx,inody) *
440 & grid_jacq_streamice(i,j,bi,bj,n) *
441 & tau_beta_eff_streamice (i,j,bi,bj) * uq
442
443 streamice_cg_A4
444 & (i-1+inodx,j-1+inody,bi,bj,col_x,col_y)=
445 & streamice_cg_A4
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) * vq
450
451 endif
452 enddo
453 enddo
454 endif
455 enddo
456 enddo
457 enddo
458 enddo
459 endif
460 enddo
461 enddo
462 enddo
463 enddo
464
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) then
547
548 ilq = 1
549 jlq = 1
550
551 if (inode.eq.iq) ilq = 2
552 if (jnode.eq.jq) jlq = 2
553 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
554
555 ux = DPhi (i,j,bi,bj,m,n,1)
556 uy = DPhi (i,j,bi,bj,m,n,2)
557 vx = 0
558 vy = 0
559 uq = Xquad(ilq) * Xquad(jlq)
560 vq = 0
561
562 exx = ux + k1AtC_str(i,j,bi,bj)*vq
563 eyy = vy + k2AtC_str(i,j,bi,bj)*uq
564 exy = .5*(uy+vx) +
565 & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
566
567 uret(i-1+inode,j-1+jnode,bi,bj) =
568 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
569 & grid_jacq_streamice(i,j,bi,bj,n) *
570 & visc_streamice(i,j,bi,bj) * (
571 & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
572 & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
573
574 uret(i-1+inode,j-1+jnode,bi,bj) =
575 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
576 & grid_jacq_streamice(i,j,bi,bj,n) *
577 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
578 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
579 & 4*0.5*k1AtC_str(i,j,bi,bj)*exy)
580
581
582 uret(i-1+inode,j-1+jnode,bi,bj) =
583 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
584 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
585 & tau_beta_eff_streamice (i,j,bi,bj) * uq
586
587
588 vx = DPhi (i,j,bi,bj,m,n,1)
589 vy = DPhi (i,j,bi,bj,m,n,2)
590 ux = 0
591 uy = 0
592 vq = Xquad(ilq) * Xquad(jlq)
593 uq = 0
594
595 exx = ux + k1AtC_str(i,j,bi,bj)*vq
596 eyy = vy + k2AtC_str(i,j,bi,bj)*uq
597 exy = .5*(uy+vx) +
598 & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
599
600 vret(i-1+inode,j-1+jnode,bi,bj) =
601 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
602 & grid_jacq_streamice(i,j,bi,bj,n) *
603 & visc_streamice(i,j,bi,bj) * (
604 & DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) +
605 & DPhi(i,j,bi,bj,m,n,1)*(2*exy))
606 vret(i-1+inode,j-1+jnode,bi,bj) =
607 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
608 & grid_jacq_streamice(i,j,bi,bj,n) *
609 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
610 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
611 & 4*0.5*k2AtC_str(i,j,bi,bj)*exy)
612
613
614 vret(i-1+inode,j-1+jnode,bi,bj) =
615 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
616 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
617 & tau_beta_eff_streamice (i,j,bi,bj) * vq
618
619 endif
620
621 enddo
622 enddo
623 enddo
624 enddo
625 endif
626 enddo
627 enddo
628 enddo
629 enddo
630
631 #endif
632 RETURN
633 END SUBROUTINE
634
635
636
637 SUBROUTINE STREAMICE_CG_BOUND_VALS( myThid,
638 O uret,
639 O vret)
640 C /============================================================\
641 C | SUBROUTINE |
642 C | o |
643 C |============================================================|
644 C | |
645 C \============================================================/
646 IMPLICIT NONE
647
648 C === Global variables ===
649 #include "SIZE.h"
650 #include "EEPARAMS.h"
651 #include "PARAMS.h"
652 #include "GRID.h"
653 #include "STREAMICE.h"
654 #include "STREAMICE_CG.h"
655
656 C !INPUT/OUTPUT ARGUMENTS
657 C uret, vret - result of matrix operating on u, v
658 C is, ie, js, je - starting and ending cells
659 INTEGER myThid
660 _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
661 _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
662
663 #ifdef ALLOW_STREAMICE
664
665 C the linear action of the matrix on (u,v) with triangular finite elements
666 C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
667 C but this may change pursuant to conversations with others
668 C
669 C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
670 C in order to make less frequent halo updates
671 C isym = 1 if grid is symmetric, 0 o.w.
672
673 C the linear action of the matrix on (u,v) with triangular finite elements
674 C Phi has the form
675 C Phi (i,j,k,q) - applies to cell i,j
676
677 C 3 - 4
678 C | |
679 C 1 - 2
680
681 C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
682 C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
683 C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
684
685 C !LOCAL VARIABLES:
686 C == Local variables ==
687 INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
688 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
689 _RL Ucell (2,2)
690 _RL Vcell (2,2)
691 _RL Hcell (2,2)
692 _RL phival(2,2)
693
694 uret(1,1,1,1) = uret(1,1,1,1)
695 vret(1,1,1,1) = vret(1,1,1,1)
696
697 DO j = 0, sNy+1
698 DO i = 0, sNx+1
699 DO bj = myByLo(myThid), myByHi(myThid)
700 DO bi = myBxLo(myThid), myBxHi(myThid)
701 IF ((STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) .AND.
702 & ((STREAMICE_umask(i,j,bi,bj).eq.3.0) .OR.
703 & (STREAMICE_umask(i,j+1,bi,bj).eq.3.0) .OR.
704 & (STREAMICE_umask(i+1,j,bi,bj).eq.3.0) .OR.
705 & (STREAMICE_umask(i+1,j+1,bi,bj).eq.3.0))) THEN
706
707 DO iq=1,2
708 DO jq = 1,2
709
710 n = 2*(jq-1)+iq
711
712 uq = u_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+
713 & u_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+
714 & u_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+
715 & u_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
716 vq = v_bdry_values_SI(i,j,bi,bj)*Xquad(3-iq)*Xquad(3-jq)+
717 & v_bdry_values_SI(i+1,j,bi,bj)*Xquad(iq)*Xquad(3-jq)+
718 & v_bdry_values_SI(i,j+1,bi,bj)*Xquad(3-iq)*Xquad(jq)+
719 & v_bdry_values_SI(i+1,j+1,bi,bj)*Xquad(iq)*Xquad(jq)
720 ux = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
721 & u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
722 & u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
723 & u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
724 uy = u_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
725 & u_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
726 & u_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
727 & u_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
728 vx = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
729 & v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
730 & v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
731 & v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
732 vy = v_bdry_values_SI(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
733 & v_bdry_values_SI(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
734 & v_bdry_values_SI(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
735 & v_bdry_values_SI(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
736 exx = ux + k1AtC_str(i,j,bi,bj)*vq
737 eyy = vy + k2AtC_str(i,j,bi,bj)*uq
738 exy = .5*(uy+vx) +
739 & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
740
741 do inode = 1,2
742 do jnode = 1,2
743
744 m = 2*(jnode-1)+inode
745 ilq = 1
746 jlq = 1
747 if (inode.eq.iq) ilq = 2
748 if (jnode.eq.jq) jlq = 2
749 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
750
751 if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
752
753 uret(i-1+inode,j-1+jnode,bi,bj) =
754 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
755 & grid_jacq_streamice(i,j,bi,bj,n) *
756 & visc_streamice(i,j,bi,bj) * (
757 & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
758 & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
759 vret(i-1+inode,j-1+jnode,bi,bj) =
760 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
761 & grid_jacq_streamice(i,j,bi,bj,n) *
762 & visc_streamice(i,j,bi,bj) * (
763 & DPhi(i,j,bi,bj,m,n,2)*(4*eyy+2*exx) +
764 & DPhi(i,j,bi,bj,m,n,1)*(2*exy))
765
766 uret(i-1+inode,j-1+jnode,bi,bj) =
767 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
768 & grid_jacq_streamice(i,j,bi,bj,n) *
769 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
770 & (4*k2AtC_str(i,j,bi,bj)*eyy+2*k2AtC_str(i,j,bi,bj)*exx+
771 & 4*0.5*k1AtC_str(i,j,bi,bj)*exy)
772 vret(i-1+inode,j-1+jnode,bi,bj) =
773 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
774 & grid_jacq_streamice(i,j,bi,bj,n) *
775 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
776 & (4*k1AtC_str(i,j,bi,bj)*exx+2*k1AtC_str(i,j,bi,bj)*eyy+
777 & 4*0.5*k2AtC_str(i,j,bi,bj)*exy)
778
779 ! if (STREAMICE_float_cond(i,j,bi,bj) .eq. 1) then
780 uret(i-1+inode,j-1+jnode,bi,bj) =
781 & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
782 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
783 & tau_beta_eff_streamice (i,j,bi,bj) * uq
784 vret(i-1+inode,j-1+jnode,bi,bj) =
785 & vret(i-1+inode,j-1+jnode,bi,bj) + .25 *
786 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
787 & tau_beta_eff_streamice (i,j,bi,bj) * vq
788 ! endif
789 endif
790 enddo
791 enddo
792 enddo
793 enddo
794 endif
795 enddo
796 enddo
797 enddo
798 enddo
799
800 #endif
801 RETURN
802 END SUBROUTINE

  ViewVC Help
Powered by ViewVC 1.1.22