/[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.8 - (hide annotations) (download)
Wed Aug 27 19:29:13 2014 UTC (10 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.7: +199 -54 lines
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

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

  ViewVC Help
Powered by ViewVC 1.1.22