/[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.7 - (hide annotations) (download)
Thu Mar 7 15:23:19 2013 UTC (12 years, 4 months ago) by dgoldberg
Branch: MAIN
Changes since 1.6: +6 -2 lines
bug fixes, GL smoothing, changes for controlling bathym with const surf elev

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

  ViewVC Help
Powered by ViewVC 1.1.22