/[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.2 - (hide 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 heimbach 1.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 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     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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
144 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
150 heimbach 1.1 & (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 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) * uq
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     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 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
238     _RL phival(2,2)
239 heimbach 1.1
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 heimbach 1.2 phival(inodx,inody) = Xquad(ilqx)*Xquad(ilqy)
291 heimbach 1.1
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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
352 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
362 heimbach 1.1 & (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 heimbach 1.2 & .25*phival(inodx,inody) *
370     & grid_jacq_streamice(i,j,bi,bj,n) *
371 heimbach 1.1 & 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 heimbach 1.2 & .25*phival(inodx,inody) *
378     & grid_jacq_streamice(i,j,bi,bj,n) *
379 heimbach 1.1 & 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
422 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inodx,inody) *
432 heimbach 1.1 & (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 heimbach 1.2 & .25*phival(inodx,inody) *
440     & grid_jacq_streamice(i,j,bi,bj,n) *
441 heimbach 1.1 & 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 heimbach 1.2 & .25*phival(inodx,inody) *
448     & grid_jacq_streamice(i,j,bi,bj,n) *
449 heimbach 1.1 & 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 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     if (STREAMICE_umask(i-1+inode,j-1+jnode,bi,bj).eq.1.0) then
547    
548     ilq = 1
549     jlq = 1
550 heimbach 1.1
551 heimbach 1.2 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 heimbach 1.1
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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
578 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
585 heimbach 1.1 & 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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
610 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
617 heimbach 1.1 & tau_beta_eff_streamice (i,j,bi,bj) * vq
618    
619     endif
620 heimbach 1.2
621 heimbach 1.1 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 heimbach 1.2 _RL ux, vx, uy, vy, uq, vq, exx, eyy, exy
689 heimbach 1.1 _RL Ucell (2,2)
690     _RL Vcell (2,2)
691     _RL Hcell (2,2)
692 heimbach 1.2 _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 heimbach 1.1
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 heimbach 1.2 jlq = 1
747 heimbach 1.1 if (inode.eq.iq) ilq = 2
748     if (jnode.eq.jq) jlq = 2
749 heimbach 1.2 phival(inode,jnode) = Xquad(ilq)*Xquad(jlq)
750 heimbach 1.1
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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
770 heimbach 1.1 & (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 heimbach 1.2 & visc_streamice(i,j,bi,bj) * phival(inode,jnode) *
776 heimbach 1.1 & (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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
783 heimbach 1.1 & 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 heimbach 1.2 & phival(inode,jnode) * grid_jacq_streamice(i,j,bi,bj,n) *
787 heimbach 1.1 & 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