/[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.1 - (hide annotations) (download)
Thu Mar 29 15:59:21 2012 UTC (13 years, 3 months ago) by heimbach
Branch: MAIN
Initial check-in of Dan Goldberg's streamice package

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

  ViewVC Help
Powered by ViewVC 1.1.22