/[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.6 - (hide annotations) (download)
Fri Dec 28 23:54:02 2012 UTC (12 years, 6 months ago) by dgoldberg
Branch: MAIN
Changes since 1.5: +13 -6 lines
fixing bugs related to masks

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

  ViewVC Help
Powered by ViewVC 1.1.22