/[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.4 - (hide annotations) (download)
Thu Jul 26 16:13:18 2012 UTC (12 years, 11 months ago) by dgoldberg
Branch: MAIN
Changes since 1.3: +1 -6 lines
replace print statements with print_message calls

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

  ViewVC Help
Powered by ViewVC 1.1.22