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

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

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


Revision 1.5 - (hide 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 dgoldberg 1.5 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_cg_functions.F,v 1.4 2012/07/26 16:13:18 dgoldberg Exp $
2 heimbach 1.1 C $Name: $
3    
4     #include "STREAMICE_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9     SUBROUTINE STREAMICE_CG_ACTION( myThid,
10     O uret,
11     O vret,
12     I u,
13     I v,
14     I is, ie, js, je )
15     C /============================================================\
16     C | SUBROUTINE |
17     C | o |
18     C |============================================================|
19     C | |
20     C \============================================================/
21     IMPLICIT NONE
22    
23     C === Global variables ===
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "STREAMICE.h"
29     #include "STREAMICE_CG.h"
30    
31     C !INPUT/OUTPUT ARGUMENTS
32     C uret, vret - result of matrix operating on u, v
33     C is, ie, js, je - starting and ending cells
34     INTEGER myThid
35     _RL uret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
36     _RL vret (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
37     _RL u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
38     _RL v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
39     INTEGER is, ie, js, je
40    
41     #ifdef ALLOW_STREAMICE
42    
43     C the linear action of the matrix on (u,v) with triangular finite elements
44     C as of now everything is passed in so no grid pointers or anything of the sort have to be dereferenced,
45     C but this may change pursuant to conversations with others
46     C
47     C is & ie are the cells over which the iteration is done; this may change between calls to this subroutine
48     C in order to make less frequent halo updates
49     C isym = 1 if grid is symmetric, 0 o.w.
50    
51     C the linear action of the matrix on (u,v) with triangular finite elements
52     C Phi has the form
53     C Phi (i,j,k,q) - applies to cell i,j
54    
55     C 3 - 4
56     C | |
57     C 1 - 2
58    
59     C Phi (i,j,2*k-1,q) gives d(Phi_k)/dx at quadrature point q
60     C Phi (i,j,2*k,q) gives d(Phi_k)/dy at quadrature point q
61     C Phi_k is equal to 1 at vertex k, and 0 at vertex l .ne. k, and bilinear
62    
63     C !LOCAL VARIABLES:
64     C == Local variables ==
65     INTEGER iq, jq, inode, jnode, i, j, bi, bj, ilq, jlq, m, n
66 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
67 heimbach 1.1 _RL Ucell (2,2)
68     _RL Vcell (2,2)
69     _RL Hcell (2,2)
70 heimbach 1.2 _RL phival(2,2)
71    
72     uret(1,1,1,1) = uret(1,1,1,1)
73     vret(1,1,1,1) = vret(1,1,1,1)
74 heimbach 1.1
75     DO j = js, je
76     DO i = is, ie
77     DO bj = myByLo(myThid), myByHi(myThid)
78     DO bi = myBxLo(myThid), myBxHi(myThid)
79 heimbach 1.2
80 heimbach 1.1 IF (STREAMICE_hmask (i,j,bi,bj) .eq. 1.0) THEN
81 heimbach 1.2 DO iq = 1,2
82 heimbach 1.1 DO jq = 1,2
83    
84     n = 2*(jq-1)+iq
85    
86     uq = u(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) +
87     & u(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) +
88     & u(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) +
89     & u(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
90     vq = v(i,j,bi,bj) * Xquad(3-iq) * Xquad(3-jq) +
91     & v(i+1,j,bi,bj) * Xquad(iq) * Xquad(3-jq) +
92     & v(i,j+1,bi,bj) * Xquad(3-iq) * Xquad(jq) +
93     & v(i+1,j+1,bi,bj) * Xquad(iq) * Xquad(jq)
94     ux = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
95     & u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
96     & u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
97     & u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
98     uy = u(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) +
99     & u(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
100     & u(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
101     & u(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
102     vx = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,1) +
103     & v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,1) +
104     & v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,1) +
105     & v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,1)
106     vy = v(i,j,bi,bj) * DPhi(i,j,bi,bj,1,n,2) +
107     & v(i+1,j,bi,bj) * DPhi(i,j,bi,bj,2,n,2) +
108     & v(i,j+1,bi,bj) * DPhi(i,j,bi,bj,3,n,2) +
109     & v(i+1,j+1,bi,bj) * DPhi(i,j,bi,bj,4,n,2)
110     exx = ux + k1AtC_str(i,j,bi,bj)*vq
111     eyy = vy + k2AtC_str(i,j,bi,bj)*uq
112     exy = .5*(uy+vx) +
113     & k1AtC_str(i,j,bi,bj)*uq + k2AtC_str(i,j,bi,bj)*vq
114    
115     do inode = 1,2
116     do jnode = 1,2
117    
118     m = 2*(jnode-1)+inode
119     ilq = 1
120 heimbach 1.2 jlq = 1
121 heimbach 1.1 if (inode.eq.iq) ilq = 2
122 heimbach 1.2 if (jnode.eq.jq) jlq = 2
123     phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
124 heimbach 1.1
125     if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
126    
127     uret(i-1+inode,j-1+jnode,bi,bj) =
128     & uret(i-1+inode,j-1+jnode,bi,bj) + .25 *
129     & grid_jacq_streamice(i,j,bi,bj,n) *
130     & visc_streamice(i,j,bi,bj) * (
131     & DPhi(i,j,bi,bj,m,n,1)*(4*exx+2*eyy) +
132     & DPhi(i,j,bi,bj,m,n,2)*(2*exy))
133 dgoldberg 1.5 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 heimbach 1.1 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
158 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) *
163     & grid_jacq_streamice(i,j,bi,bj,n) *
164 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
165    
166     endif
167     enddo
168     enddo
169 heimbach 1.2
170 heimbach 1.1 enddo
171     enddo
172 heimbach 1.2 c-- STREAMICE_hmask
173 heimbach 1.1 endif
174 heimbach 1.2
175 heimbach 1.1 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 dgoldberg 1.3 #ifdef STREAMICE_CONSTRUCT_MATRIX
209    
210 heimbach 1.1 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 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
235     _RL phival(2,2)
236 heimbach 1.1
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 dgoldberg 1.5 & .eq.1.0 .or.
279     & streamice_vmask(i-1+inodx,j-1+inody,bi,bj).eq.1.0)
280 heimbach 1.1 & 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 heimbach 1.2 phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)
289 heimbach 1.1
290     DO jnodx = 1,2
291     DO jnody = 1,2
292     if (STREAMICE_umask(i-1+jnodx,j-1+jnody,bi,bj)
293 dgoldberg 1.5 & .eq.1.0 .or.
294     & STREAMICE_vmask(i-1+jnodx,j-1+jnody,bi,bj).eq.1.0)
295 heimbach 1.1 & 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
351 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
361 heimbach 1.1 & (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 heimbach 1.2 & .25*phival(inodx,inody) *
369     & grid_jacq_streamice(i,j,bi,bj,n) *
370 heimbach 1.1 & 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 heimbach 1.2 & .25*phival(inodx,inody) *
377     & grid_jacq_streamice(i,j,bi,bj,n) *
378 heimbach 1.1 & 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
421 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
431 heimbach 1.1 & (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 heimbach 1.2 & .25*phival(inodx,inody) *
439     & grid_jacq_streamice(i,j,bi,bj,n) *
440 heimbach 1.1 & 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 heimbach 1.2 & .25*phival(inodx,inody) *
447     & grid_jacq_streamice(i,j,bi,bj,n) *
448 heimbach 1.1 & 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 dgoldberg 1.3 #endif
466 heimbach 1.1 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 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
523 heimbach 1.1 _RL Ucell (2,2)
524     _RL Vcell (2,2)
525     _RL Hcell (2,2)
526 heimbach 1.2 _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 heimbach 1.1
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 heimbach 1.2
546 dgoldberg 1.5 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 heimbach 1.2
550     ilq = 1
551     jlq = 1
552 heimbach 1.1
553 heimbach 1.2 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 heimbach 1.1
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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
580 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
587 heimbach 1.1 & 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
612 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
619 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
620    
621     endif
622 heimbach 1.2
623 heimbach 1.1 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 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
691 heimbach 1.1 _RL Ucell (2,2)
692     _RL Vcell (2,2)
693     _RL Hcell (2,2)
694 heimbach 1.2 _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 heimbach 1.1
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 dgoldberg 1.5 & (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 heimbach 1.1
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 dgoldberg 1.5
748 heimbach 1.1 do inode = 1,2
749     do jnode = 1,2
750    
751     m = 2*(jnode-1)+inode
752     ilq = 1
753 heimbach 1.2 jlq = 1
754 heimbach 1.1 if (inode.eq.iq) ilq = 2
755     if (jnode.eq.jq) jlq = 2
756 heimbach 1.2 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
757 heimbach 1.1
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 dgoldberg 1.5
767 heimbach 1.1
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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
772 heimbach 1.1 & (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 dgoldberg 1.5
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 heimbach 1.1 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
795 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
800 heimbach 1.1 & 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