Line data Source code
1 : module edyn_mud
2 : use shr_kind_mod,only: r8 => shr_kind_r8
3 : use cam_abortutils,only: endrun
4 : use edyn_mudcom, only: cor2, res2, factri, factrp, prolon2, trsfc2, swk2
5 :
6 : implicit none
7 :
8 : private
9 :
10 : public :: mud2cr1
11 : public :: dismd2cr
12 : public :: adjmd2cr
13 : public :: kcymd2cr
14 : public :: relmd2cr
15 : public :: resmd2cr
16 :
17 : contains
18 : !-----------------------------------------------------------------------
19 : !
20 : ! file mud2cr.f (version 4.0 modified for Cicley 2/99)
21 : ! .
22 : ! . MUDPACK version 4.0 .
23 : ! . .
24 : ! ... author and specialist
25 : !
26 : ! John C. Adams (National Center for Atmospheric Research) (retired)
27 : !
28 : ! ... For MUDPACK information, visit the website:
29 : ! (https://www2.cisl.ucar.edu/resources/legacy/mudpack)
30 : !
31 : ! ... purpose
32 : !
33 : ! mud2cr attempts to produce a second order finite difference
34 : ! approximation to the two dimensional nonseparable elliptic
35 : ! partial differential equation with cross derivative
36 : !
37 : ! cxx(x,y)*pxx + cxy(x,y)*pxy + cyy(x,y)*pyy +
38 : !
39 : ! cx(x,y)*px + cy(x,y)*py + ce(x,y)*pe(x,y) = r(x,y)
40 : !
41 : ! ... documentation
42 : !
43 : ! see the documentation on above website for a complete discussion
44 : ! of how to use subroutine mud2cr.
45 : !
46 : ! ... required MUDPACK files
47 : !
48 : ! mudcom.f
49 : !
50 : !
51 : !
52 33516 : subroutine mud2cr(iparm,fparm,work,rhs,phi,mgopt, &
53 : ierror,isolve)
54 :
55 : integer,intent(in) :: isolve
56 : integer iparm,mgopt,ierror
57 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
58 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
59 : kcycle,iprer,ipost,intpol,kps
60 : real(r8) :: fparm,xa,xb,yc,yd,tolmax,relmax
61 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
62 : integer int,iw,k,kb,nx,ny,ic,itx,ity
63 : dimension iparm(17),fparm(6),mgopt(4)
64 : real(r8) :: work(*),phi(*),rhs(*)
65 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
66 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
67 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
68 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
69 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50), &
70 : nxk(50),nyk(50),isx,jsy
71 : data int / 0 /
72 : save int
73 : !
74 : ierror = 1
75 : intl = iparm(1) ! set and check intl on all calls
76 : if (intl*(intl-1).ne.0) return
77 : if (int.eq.0) then
78 : int = 1
79 : if (intl.ne.0) return ! very first call is not intl=0
80 : end if
81 : ierror = 0
82 : !
83 : ! set arguments internally
84 : ! these will not be rechecked if intl=1!
85 : !
86 : nxa = iparm(2)
87 : nxb = iparm(3)
88 : nyc = iparm(4)
89 : nyd = iparm(5)
90 : ixp = iparm(6)
91 : jyq = iparm(7)
92 : iex = iparm(8)
93 : jey = iparm(9)
94 : ngrid = max0(iex,jey)
95 : nfx = iparm(10)
96 : nfy = iparm(11)
97 : iguess = iparm(12)
98 : maxcy = iparm(13)
99 : method = iparm(14)
100 : nwork = iparm(15)
101 : kcycle = mgopt(1)
102 : if (kcycle .eq. 0) then
103 : ! set defaults
104 : kcycle = 2
105 : iprer = 2
106 : ipost = 1
107 : intpol = 3
108 : else
109 : iprer = mgopt(2)
110 : ipost = mgopt(3)
111 : intpol = mgopt(4)
112 : end if
113 : xa = fparm(1)
114 : xb = fparm(2)
115 : yc = fparm(3)
116 : yd = fparm(4)
117 : tolmax = fparm(5)
118 : if (intl .eq. 0) then ! intialization call
119 : !
120 : ! check input arguments
121 : !
122 : ierror = 2 ! check boundary condition flags
123 : if (max0(nxa,nxb,nyc,nyd).gt.2) return
124 : if (min0(nxa,nxb,nyc,nyd).lt.0) return
125 : if (nxa.eq.0.and.nxb.ne.0) return
126 : if (nxa.ne.0.and.nxb.eq.0) return
127 : if (nyc.eq.0.and.nyd.ne.0) return
128 : if (nyc.ne.0.and.nyd.eq.0) return
129 : ierror = 3 ! check grid sizes
130 : if (ixp.lt.2) return
131 : if (jyq.lt.2) return
132 : ierror = 4
133 : ngrid = max0(iex,jey)
134 : if (iex.lt.1) return
135 : if (jey.lt.1) return
136 : if (ngrid.gt.50) return
137 : ierror = 5
138 : if (nfx.ne.ixp*2**(iex-1)+1) return
139 : if (nfy.ne.jyq*2**(jey-1)+1) return
140 : ierror = 6
141 : if (iguess*(iguess-1).ne.0) return
142 : ierror = 7
143 : if (maxcy.lt.1) return
144 : ierror = 8
145 : if (method.lt.0 .or. method.gt.3) return
146 : ierror = 9
147 : ! compute and test minimum work space
148 : isx = 0
149 : if (method.eq.1 .or. method.eq.3) then
150 : if (nxa.ne.0) isx = 3
151 : if (nxa.eq.0) isx = 5
152 : end if
153 : jsy = 0
154 : if (method.eq.2 .or. method.eq.3) then
155 : if (nyc.ne.0) jsy = 3
156 : if (nyc.eq.0) jsy = 5
157 : end if
158 : kps = 1
159 : do k=1,ngrid
160 : ! set subgrid sizes
161 : nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1
162 : nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1
163 : nx = nxk(k)
164 : ny = nyk(k)
165 : kps = kps+(nx+2)*(ny+2)+nx*ny*(10+isx+jsy)
166 : end do
167 : iparm(16) = kps+(nfx+2)*(nfy+2) ! exact minimum work space
168 : lwork = iparm(16)
169 : if (lwork .gt. nwork) return
170 : ierror = 10 ! check solution region
171 : if (xb.le.xa .or. yd.le.yc) return
172 : ierror = 11
173 : if (tolmax .lt. 0.0_r8) return
174 : ierror = 12 ! multigrid parameters
175 : if (kcycle.lt.0) return
176 : if (min0(iprer,ipost).lt.1) return
177 : if ((intpol-1)*(intpol-3).ne.0) return
178 : if (max0(kcycle,iprer,ipost).gt.2) then
179 : ierror = -5 ! inefficient multigrid cycling
180 : end if
181 : if (ierror .gt. 0) ierror = 0 ! no fatal errors
182 : !
183 : ! set work space pointers and discretize pde at each grid level
184 : !
185 : iw = 1
186 : do kb=1,ngrid
187 : k = ngrid-kb+1
188 : nx = nxk(k)
189 : ny = nyk(k)
190 : kpbgn(k) = iw
191 : kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)
192 : ktxbgn(k) = kcbgn(k)+10*nx*ny
193 : ktybgn(k) = ktxbgn(k)+isx*nx*ny
194 : iw = ktybgn(k)+jsy*nx*ny
195 : ic = kcbgn(k)
196 : itx = ktxbgn(k)
197 : ity = ktybgn(k)
198 : klevel = k
199 : call dismd2cr(nx,ny,work(ic),work(itx),work(ity), &
200 : work,ierror,isolve)
201 : end do
202 : return
203 : end if ! end of intl=0 initialization call block
204 : nx = nfx
205 : ny = nfy
206 : call mud2cr1(nx,ny,rhs,phi,work)
207 : iparm(17) = itero
208 : if (tolmax.gt.0.0_r8) then ! check for convergence
209 : fparm(6) = relmax
210 : if (relmax.gt.tolmax) ierror = -1 ! flag convergenc failure
211 : end if
212 : return
213 : end subroutine mud2cr
214 : !-----------------------------------------------------------------------
215 0 : subroutine mud2cr1(nx,ny,rhsf,phif,wk)
216 :
217 : integer nx,ny
218 : real(r8) :: phif(nx,ny),rhsf(nx,ny),wk(*)
219 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
220 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
221 : kcycle,iprer,ipost,intpol,kps
222 : real(r8) :: xa,xb,yc,yd,tolmax,relmax,phmax
223 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
224 : integer k,kb,ip,ic,ir,ipc,irc,icc
225 : integer ncx,ncy,jj,ij,i,j,iter
226 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
227 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
228 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
229 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
230 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),&
231 : nxk(50),nyk(50),isx,jsy
232 0 : nx = nxk(ngrid)
233 0 : ny = nyk(ngrid)
234 0 : ip = kpbgn(ngrid)
235 0 : ic = kcbgn(ngrid)
236 0 : ir = ic+9*nx*ny
237 : !
238 : ! set phif,rhsf in wk and adjust right hand side
239 : !
240 0 : call swk2(nx,ny,phif,rhsf,wk(ip),wk(ir))
241 0 : if (iguess.eq.0) then
242 : !
243 : ! no initial guess at finest grid level!
244 : !
245 0 : do kb=2,ngrid
246 0 : k = ngrid-kb+1
247 0 : nx = nxk(k+1)
248 0 : ny = nyk(k+1)
249 0 : ip = kpbgn(k+1)
250 0 : ir = kcbgn(k+1)+9*nx*ny
251 0 : ncx = nxk(k)
252 0 : ncy = nyk(k)
253 0 : ipc = kpbgn(k)
254 0 : icc = kcbgn(k)
255 0 : irc = icc+9*ncx*ncy
256 : !
257 : ! transfer down to all grid levels
258 : !
259 0 : call trsfc2(nx,ny,wk(ip),wk(ir),ncx,ncy,&
260 0 : wk(ipc),wk(irc))
261 : end do
262 : !
263 : ! adjust right hand side at all grid levels in case
264 : ! rhs or specified b.c. in phi or gbdy changed
265 : !
266 0 : do k=1,ngrid
267 0 : nx = nxk(k)
268 0 : ny = nyk(k)
269 0 : ip = kpbgn(k)
270 0 : ic = kcbgn(k)
271 0 : call adjmd2cr(nx,ny,wk(ip),wk(ic))
272 : end do
273 : !
274 : ! execute one full multigrid cycle
275 : !
276 0 : do k=1,ngrid-1
277 0 : kcur = k
278 0 : call kcymd2cr(wk)
279 0 : nx = nxk(k+1)
280 0 : ny = nyk(k+1)
281 0 : ip = kpbgn(k+1)
282 0 : ipc = kpbgn(k)
283 0 : ncx = nxk(k)
284 0 : ncy = nyk(k)
285 : !
286 : ! lift or prolong approximation from k to k+1
287 : !
288 0 : call prolon2(ncx,ncy,wk(ipc),nx,ny,wk(ip),nxa,nxb,&
289 0 : nyc,nyd,intpol)
290 : end do
291 : else
292 : !
293 : ! adjust rhs at finest grid level only
294 : !
295 0 : nx = nxk(ngrid)
296 0 : ny = nyk(ngrid)
297 0 : ip = kpbgn(ngrid)
298 0 : ic = kcbgn(ngrid)
299 0 : call adjmd2cr(nx,ny,wk(ip),wk(ic))
300 : end if
301 : !
302 : ! execute maxcy more multigrid k cycles from finest level
303 : !
304 0 : kcur = ngrid
305 0 : do iter=1,maxcy
306 0 : itero = iter
307 0 : call kcymd2cr(wk)
308 0 : if (tolmax.gt.0.0_r8) then
309 : !
310 : ! error control
311 : !
312 0 : relmax = 0.0_r8
313 0 : phmax = 0.0_r8
314 0 : do j=1,nfy
315 0 : jj = j*(nfx+2)
316 0 : do i=1,nfx
317 0 : ij = jj+i+1
318 : ! phmax = amax1(phmax,abs(wk(ij)))
319 : ! relmax = amax1(relmax,abs(wk(ij)-phif(i,j)))
320 0 : phmax = max(phmax,abs(wk(ij)))
321 0 : relmax = max(relmax,abs(wk(ij)-phif(i,j)))
322 0 : phif(i,j) = wk(ij)
323 : end do
324 : end do
325 : !
326 : ! set maximum relative difference and check for convergence
327 : !
328 0 : if (phmax.gt.0.0_r8) relmax = relmax/phmax
329 0 : if (relmax.le.tolmax) return
330 : end if
331 : end do
332 : !
333 : ! set final interate after maxcy cycles in phif
334 : !
335 0 : do j=1,nfy
336 0 : jj = j*(nfx+2)
337 0 : do i=1,nfx
338 0 : ij = jj+i+1
339 0 : phif(i,j) = wk(ij)
340 : end do
341 : end do
342 : return
343 : end subroutine mud2cr1
344 : !-----------------------------------------------------------------------
345 133 : subroutine kcymd2cr(wk)
346 : !
347 : ! execute multigrid k cycle from kcur grid level
348 : ! kcycle=1 for v cycles, kcycle=2 for w cycles
349 : !
350 : real(r8) :: wk(*)
351 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
352 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
353 : kcycle,iprer,ipost,intpol,kps
354 : integer nx,ny,ip,ic,ipc,irc,itx,ity,ncx,ncy,l,nrel
355 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
356 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
357 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
358 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
359 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
360 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
361 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),&
362 : nxk(50),nyk(50),isx,jsy
363 : integer kount(50)
364 76 : klevel = kcur
365 76 : nx = nxk(klevel)
366 76 : ny = nyk(klevel)
367 76 : ip = kpbgn(klevel)
368 76 : ic = kcbgn(klevel)
369 76 : itx = ktxbgn(klevel)
370 76 : ity = ktybgn(klevel)
371 : !
372 : ! prerelax at current finest grid level
373 : !
374 304 : do l=1,iprer
375 532 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
376 : end do
377 76 : if (kcur .eq. 1) go to 5
378 : !
379 : ! restrict residual to kcur-1 level
380 : !
381 57 : ipc = kpbgn(klevel-1)
382 57 : ncx = nxk(klevel-1)
383 57 : ncy = nyk(klevel-1)
384 57 : irc = kcbgn(klevel-1)+9*ncx*ncy
385 114 : call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps))
386 : !
387 : ! set counter for grid levels to zero
388 : !
389 228 : do l = 1,kcur
390 228 : kount(l) = 0
391 : end do
392 : !
393 : ! set new grid level and continue k-cycling
394 : !
395 57 : klevel = kcur-1
396 57 : nrel = iprer
397 : !
398 : ! kcycle control point
399 : !
400 : 10 continue
401 : !
402 : ! post relax when kcur revisited
403 : !
404 418 : if (klevel .eq. kcur) go to 5
405 : !
406 : ! count hit at current level
407 : !
408 361 : kount(klevel) = kount(klevel)+1
409 : !
410 : ! relax at current level
411 : !
412 361 : nx = nxk(klevel)
413 361 : ny = nyk(klevel)
414 361 : ip = kpbgn(klevel)
415 361 : ic = kcbgn(klevel)
416 361 : itx = ktxbgn(klevel)
417 361 : ity = ktybgn(klevel)
418 1292 : do l=1,nrel
419 2223 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
420 : end do
421 361 : if (kount(klevel) .eq. kcycle+1) then
422 : !
423 : ! kcycle complete at klevel
424 : !
425 95 : ipc = ip
426 95 : ip = kpbgn(klevel+1)
427 95 : ncx = nxk(klevel)
428 95 : ncy = nyk(klevel)
429 95 : nx = nxk(klevel+1)
430 95 : ny = nyk(klevel+1)
431 : !
432 : ! inject correction to finer grid
433 : !
434 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,&
435 190 : intpol,wk(kps))
436 : !
437 : ! reset counter to zero
438 : !
439 95 : kount(klevel) = 0
440 : !
441 : ! ascend to next higher level and set to postrelax there
442 : !
443 95 : klevel = klevel+1
444 95 : nrel = ipost
445 95 : go to 10
446 : else
447 266 : if (klevel .gt. 1) then
448 : !
449 : ! kcycle not complete so descend unless at coarsest grid
450 : !
451 152 : ipc = kpbgn(klevel-1)
452 152 : ncx = nxk(klevel-1)
453 152 : ncy = nyk(klevel-1)
454 152 : irc = kcbgn(klevel-1)+9*ncx*ncy
455 0 : call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),&
456 304 : wk(kps))
457 : !
458 : ! prerelax at next coarser level
459 : !
460 152 : klevel = klevel-1
461 152 : nrel = iprer
462 152 : go to 10
463 : else
464 : !
465 : ! postrelax at coarsest level
466 : !
467 342 : do l=1,ipost
468 570 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
469 : end do
470 114 : ipc = ip
471 114 : ip = kpbgn(2)
472 114 : ncx = nxk(1)
473 114 : ncy = nyk(1)
474 114 : nx = nxk(2)
475 114 : ny = nyk(2)
476 : !
477 : ! inject correction to level 2
478 : !
479 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,&
480 228 : intpol,wk(kps))
481 : !
482 : ! set to postrelax at level 2
483 : !
484 114 : nrel = ipost
485 114 : klevel = 2
486 190 : go to 10
487 : end if
488 : end if
489 : 5 continue
490 : !
491 : ! post relax at current finest grid level
492 : !
493 76 : nx = nxk(kcur)
494 76 : ny = nyk(kcur)
495 76 : ip = kpbgn(kcur)
496 76 : ic = kcbgn(kcur)
497 76 : itx = ktxbgn(kcur)
498 76 : ity = ktybgn(kcur)
499 228 : do l=1,ipost
500 380 : call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
501 : end do
502 76 : return
503 0 : end subroutine kcymd2cr
504 : !-----------------------------------------------------------------------
505 95 : subroutine dismd2cr(nx,ny,cf,tx,ty,wk,ier,isolve)
506 3572 : use edyn_solver_coefs,only: nc,cee,ceee
507 : use edyn_maggrid, only: res_nlev
508 : !
509 : ! discretize elliptic pde for mud2cr, set nonfatal errors
510 : !
511 : integer,intent(in) :: isolve
512 : integer nx,ny,i,j,l,im1,jm1,ier,nnx,nny
513 : real(r8) :: cf(nx,ny,10),tx(nx,ny,*),ty(ny,nx,*)
514 : real(r8) :: wk(*)
515 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
516 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
517 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
518 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
519 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
520 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
521 :
522 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
523 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
524 : !
525 : ! CHECK FOR CONSISTENCYT WRT KLEVEL
526 : !
527 95 : NNX = ixp*2**(KLEVEL-1)+1
528 95 : NNY = jyq*2**(KLEVEL-1)+1
529 95 : IF(NNX.NE.NX.OR.NNY.NE.NY)THEN
530 0 : call endrun('dismd2cr in mud')
531 : ENDIF
532 95 : if (isolve >= 0) then
533 95 : call ceee(cee(nc(res_nlev+1-klevel)),nx,ny,cf)
534 : endif
535 : !
536 : ! set coefficient for specified boundaries
537 : !
538 95 : if (nxa.eq.1) then
539 0 : i = 1
540 0 : do j=1,ny
541 0 : do l=1,9
542 0 : cf(i,j,l) = 0.0_r8
543 : end do
544 0 : cf(i,j,9) = 1.0_r8
545 : end do
546 : end if
547 95 : if (nxb.eq.1) then
548 0 : i = nx
549 0 : do j=1,ny
550 0 : do l=1,9
551 0 : cf(i,j,l) = 0.0_r8
552 : end do
553 0 : cf(i,j,9) = 1.0_r8
554 : end do
555 : end if
556 95 : if (nyc.eq.1) then
557 0 : j = 1
558 0 : do i=1,nx
559 0 : do l=1,9
560 0 : cf(i,j,l) = 0.0_r8
561 : end do
562 0 : cf(i,j,9) = 1.0_r8
563 : end do
564 : end if
565 95 : if (nyd.eq.1) then
566 95 : j = ny
567 3135 : do i=1,nx
568 30400 : do l=1,9
569 30400 : cf(i,j,l) = 0.0_r8
570 : end do
571 3135 : cf(i,j,9) = 1.0_r8
572 : end do
573 : end if
574 : !
575 : ! set and factor tridiagonal matrices for line relaxation(s) if flagged
576 : !
577 95 : if (method.eq.1.or.method.eq.3) then
578 95 : if (nxa.ne.0) then
579 : !
580 : ! nonperiodic x line relaxation
581 : !
582 0 : do i=1,nx
583 0 : im1 = max0(i-1,1)
584 0 : do j=1,ny
585 0 : tx(im1,j,1) = cf(i,j,5)
586 0 : tx(i,j,2) = cf(i,j,9)
587 0 : tx(i,j,3) = cf(i,j,1)
588 : end do
589 : end do
590 0 : call factri(ny,nx,tx(1,1,1),tx(1,1,2),tx(1,1,3))
591 : else
592 : !
593 : ! periodic x line relaxation
594 : !
595 95 : if (nx .gt. 3) then
596 : !
597 : ! set and factor iff nx > 3
598 : !
599 3040 : do i=1,nx-1
600 103170 : do j=1,ny
601 100130 : tx(i,j,1) = cf(i,j,5)
602 100130 : tx(i,j,2) = cf(i,j,9)
603 103075 : tx(i,j,3) = cf(i,j,1)
604 : end do
605 : end do
606 95 : call factrp(ny,nx,tx,tx(1,1,2),tx(1,1,3),tx(1,1,4),&
607 95 : tx(1,1,5),wk(kps))
608 : end if
609 : end if
610 : end if
611 :
612 95 : if (method.eq.2.or.method.eq.3) then
613 190 : if (nyc.ne.0) then
614 : !
615 : ! nonperiodic y line relaxation
616 : !
617 1957 : do j=1,ny
618 1862 : jm1 = max0(j-1,1)
619 103949 : do i=1,nx
620 101992 : ty(jm1,i,1) = cf(i,j,7)
621 101992 : ty(j,i,2) = cf(i,j,9)
622 103854 : ty(j,i,3) = cf(i,j,3)
623 : end do
624 : end do
625 95 : call factri(nx,ny,ty(1,1,1),ty(1,1,2),ty(1,1,3))
626 : else
627 : !
628 : ! periodic y line relaxation
629 : !
630 0 : if (ny .gt. 3) then
631 : !
632 : ! set and factor iff ny > 3
633 : !
634 0 : do j=1,ny-1
635 0 : do i=1,nx
636 0 : ty(j,i,1) = cf(i,j,7)
637 0 : ty(j,i,2) = cf(i,j,9)
638 0 : ty(j,i,3) = cf(i,j,3)
639 : end do
640 : end do
641 0 : call factrp(nx,ny,ty,ty(1,1,2),ty(1,1,3),ty(1,1,4),&
642 0 : ty(1,1,5),wk(kps))
643 : end if
644 : end if
645 : end if
646 95 : return
647 : end subroutine dismd2cr
648 : !-----------------------------------------------------------------------
649 95 : subroutine adjmd2cr(nx,ny,phi,cf)
650 : !
651 : ! adjust righthand side in cf(i,j,10) for boundary conditions
652 : !
653 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
654 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
655 : kcycle,iprer,ipost,intpol,kps
656 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
657 : integer nx,ny,i,j
658 : real(r8) :: cf(nx,ny,10),phi(0:nx+1,0:ny+1)
659 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
660 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
661 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
662 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
663 : !
664 : ! set specified boundaries in rhs from phi
665 : !
666 95 : if (nxa.eq.1) then
667 0 : i = 1
668 0 : do j=1,ny
669 0 : cf(i,j,10) = phi(i,j)
670 : end do
671 : end if
672 95 : if (nxb.eq.1) then
673 0 : i = nx
674 0 : do j=1,ny
675 0 : cf(i,j,10) = phi(i,j)
676 : end do
677 : end if
678 95 : if (nyc.eq.1) then
679 0 : j = 1
680 0 : do i=1,nx
681 0 : cf(i,j,10) = phi(i,j)
682 : end do
683 : end if
684 95 : if (nyd.eq.1) then
685 95 : j = ny
686 3135 : do i=1,nx
687 3135 : cf(i,j,10) = phi(i,j)
688 : end do
689 : end if
690 95 : return
691 665 : end subroutine adjmd2cr
692 : !-----------------------------------------------------------------------
693 1273 : subroutine resmd2cr(nx,ny,phi,ncx,ncy,phic,rhsc,cof,resf)
694 : !
695 : ! restrict residual from fine to coarse mesh using fully weighted
696 : ! residual restriction
697 : !
698 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
699 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
700 : kcycle,iprer,ipost,intpol,kps
701 : integer nx,ny,ncx,ncy,i,j,ic,jc
702 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
703 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
704 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
705 : real(r8) :: rhsc(ncx,ncy),resf(nx,ny)
706 : real(r8) :: phi(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
707 : real(r8) :: cof(nx,ny,10)
708 : !
709 : ! set phic zero
710 : !
711 11533 : do jc=0,ncy+1
712 148333 : do ic=0,ncx+1
713 147060 : phic(ic,jc) = 0.0_r8
714 : end do
715 : end do
716 : !
717 : ! compute residual on fine mesh in resf
718 : !
719 15428 : do j=1,ny
720 346313 : do i=1,nx
721 661770 : resf(i,j) = cof(i,j,10)-( &
722 661770 : cof(i,j,1)*phi(i+1,j)+ &
723 330885 : cof(i,j,2)*phi(i+1,j+1)+ &
724 330885 : cof(i,j,3)*phi(i,j+1)+ &
725 330885 : cof(i,j,4)*phi(i-1,j+1)+ &
726 : cof(i,j,5)*phi(i-1,j)+ &
727 330885 : cof(i,j,6)*phi(i-1,j-1)+ &
728 : cof(i,j,7)*phi(i,j-1)+ &
729 : cof(i,j,8)*phi(i+1,j-1)+ &
730 2330350 : cof(i,j,9)*phi(i,j))
731 : end do
732 : end do
733 : !
734 : ! restrict resf to coarse mesh in rhsc
735 : !
736 1273 : call res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
737 1273 : return
738 95 : end subroutine resmd2cr
739 : !-----------------------------------------------------------------------
740 8379 : subroutine relmd2cr(nx,ny,phi,cof,tx,ty,sum)
741 : !
742 : ! relaxation for mud2
743 : !
744 : integer nx,ny
745 : real(r8) :: phi(*),cof(*),tx(*),ty(*),sum(*)
746 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
747 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
748 : kcycle,iprer,ipost,intpol,kps
749 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
750 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
751 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
752 8379 : if (method.eq.0) then ! point relaxation
753 0 : call relmd2crp(nx,ny,phi,cof)
754 8379 : else if (method.eq.1) then ! line x relaxation
755 0 : call slxmd2cr(nx,ny,phi,cof,tx,sum)
756 8379 : else if (method.eq.2) then ! line y relaxation
757 0 : call slymd2cr(nx,ny,phi,cof,ty,sum)
758 8379 : else if (method.eq.3) then ! line x&y relaxation
759 8379 : call slxmd2cr(nx,ny,phi,cof,tx,sum)
760 8379 : call slymd2cr(nx,ny,phi,cof,ty,sum)
761 : end if
762 8379 : return
763 : end subroutine relmd2cr
764 : !-----------------------------------------------------------------------
765 0 : subroutine relmd2crp(nx,ny,phi,cof)
766 : !
767 : ! gauss-seidel four color point relaxation
768 : !
769 : integer nx,ny,i,j,lcolor,i1,i2,i3,i4,it
770 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
771 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
772 : kcycle,iprer,ipost,intpol,kps
773 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
774 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
775 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
776 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
777 0 : i1 = 1
778 0 : i2 = 4
779 0 : i3 = 3
780 0 : i4 = 2
781 : !
782 : ! sweep four colored grid points
783 : !
784 0 : do lcolor=1,4
785 : !$OMP PARALLEL DO SHARED(i1,cof,phi,nx,ny) PRIVATE(i,j)
786 0 : do j=1,ny,4
787 0 : do i=i1,nx,4
788 0 : phi(i,j) = (cof(i,j,10) - ( &
789 0 : cof(i,j,1)*phi(i+1,j) + &
790 0 : cof(i,j,2)*phi(i+1,j+1) + &
791 : cof(i,j,3)*phi(i,j+1) + &
792 0 : cof(i,j,4)*phi(i-1,j+1) + &
793 : cof(i,j,5)*phi(i-1,j) + &
794 0 : cof(i,j,6)*phi(i-1,j-1) + &
795 : cof(i,j,7)*phi(i,j-1) + &
796 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
797 : end do
798 : end do
799 : !$OMP PARALLEL DO SHARED(i2,cof,phi,nx,ny) PRIVATE(i,j)
800 0 : do j=2,ny,4
801 0 : do i=i2,nx,4
802 0 : phi(i,j) = (cof(i,j,10) - ( &
803 0 : cof(i,j,1)*phi(i+1,j) + &
804 0 : cof(i,j,2)*phi(i+1,j+1) + &
805 : cof(i,j,3)*phi(i,j+1) + &
806 0 : cof(i,j,4)*phi(i-1,j+1) + &
807 : cof(i,j,5)*phi(i-1,j) + &
808 0 : cof(i,j,6)*phi(i-1,j-1) + &
809 : cof(i,j,7)*phi(i,j-1) + &
810 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
811 : end do
812 : end do
813 : !$OMP PARALLEL DO SHARED(i3,cof,phi,nx,ny) PRIVATE(i,j)
814 0 : do j=3,ny,4
815 0 : do i=i3,nx,4
816 0 : phi(i,j) = (cof(i,j,10) - ( &
817 0 : cof(i,j,1)*phi(i+1,j) + &
818 0 : cof(i,j,2)*phi(i+1,j+1) + &
819 : cof(i,j,3)*phi(i,j+1) + &
820 0 : cof(i,j,4)*phi(i-1,j+1) + &
821 : cof(i,j,5)*phi(i-1,j) + &
822 0 : cof(i,j,6)*phi(i-1,j-1) + &
823 : cof(i,j,7)*phi(i,j-1) + &
824 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
825 : end do
826 : end do
827 : !$OMP PARALLEL DO SHARED(i4,cof,phi,nx,ny) PRIVATE(i,j)
828 0 : do j=4,ny,4
829 0 : do i=i4,nx,4
830 0 : phi(i,j) = (cof(i,j,10) - ( &
831 0 : cof(i,j,1)*phi(i+1,j) + &
832 0 : cof(i,j,2)*phi(i+1,j+1) + &
833 : cof(i,j,3)*phi(i,j+1) + &
834 0 : cof(i,j,4)*phi(i-1,j+1) + &
835 : cof(i,j,5)*phi(i-1,j) + &
836 0 : cof(i,j,6)*phi(i-1,j-1) + &
837 : cof(i,j,7)*phi(i,j-1) + &
838 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
839 : end do
840 : end do
841 : !
842 : ! set periodic virtual boundaries as necessary
843 : !
844 0 : if (nxa.eq.0) then
845 0 : do j=1,ny
846 0 : phi(0,j) = phi(nx-1,j)
847 0 : phi(nx+1,j) = phi(2,j)
848 : end do
849 : end if
850 0 : if (nyc.eq.0) then
851 0 : do i=1,nx
852 0 : phi(i,0) = phi(i,ny-1)
853 0 : phi(i,ny+1) = phi(i,2)
854 : end do
855 : end if
856 : !
857 : ! permute (i1,i2,i3,i4) for next color
858 : !
859 0 : it = i4
860 0 : i4 = i3
861 0 : i3 = i2
862 0 : i2 = i1
863 0 : i1 = it
864 : end do
865 0 : return
866 8379 : end subroutine relmd2crp
867 : !-----------------------------------------------------------------------
868 8379 : subroutine slxmd2cr(nx,ny,phi,cof,tx,sum)
869 : !
870 : ! line relaxation in the x direction (periodic or nonperiodic)
871 : !
872 :
873 : integer nx,ny,i,ib,j
874 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
875 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
876 : kcycle,iprer,ipost,intpol,kps
877 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
878 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
879 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
880 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10),tx(nx,ny,*),sum(ny)
881 : !
882 : ! replace line x with point gauss-seidel if
883 : ! x direction is periodic and nx = 3 (coarsest)
884 : !
885 8379 : if (nxa .eq. 0 .and. nx .eq. 3) then
886 0 : call relmd2crp(nx,ny,phi,cof)
887 0 : return
888 : end if
889 : !
890 : ! set periodic y virtual boundary if necessary
891 : !
892 8379 : if (nyc.eq.0) then
893 0 : do i=1,nx
894 0 : phi(i,0) = phi(i,ny-1)
895 0 : phi(i,ny+1) = phi(i,2)
896 : end do
897 : end if
898 :
899 8379 : if (nxa.ne.0) then
900 : !
901 : ! x direction not periodic, sweep odd j lines
902 : !
903 : !$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,ib,j)
904 0 : do j=1,ny,2
905 0 : do i=1,nx
906 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
907 : cof(i,j,3)*phi(i,j+1)+ &
908 0 : cof(i,j,4)*phi(i-1,j+1)+ &
909 0 : cof(i,j,6)*phi(i-1,j-1)+ &
910 : cof(i,j,7)*phi(i,j-1)+ &
911 0 : cof(i,j,8)*phi(i+1,j-1))
912 : end do
913 : !
914 : ! forward sweep
915 : !
916 0 : do i=2,nx
917 0 : phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
918 : end do
919 : !
920 : ! backward sweep
921 : !
922 0 : phi(nx,j) = phi(nx,j)/tx(nx,j,2)
923 0 : do ib=2,nx
924 0 : i = nx-ib+1
925 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
926 : end do
927 : end do
928 : !
929 : ! sweep even j lines forward and back
930 : !
931 : !$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,ib,j)
932 0 : do j=2,ny,2
933 0 : do i=1,nx
934 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
935 : cof(i,j,3)*phi(i,j+1)+ &
936 0 : cof(i,j,4)*phi(i-1,j+1)+ &
937 0 : cof(i,j,6)*phi(i-1,j-1)+ &
938 : cof(i,j,7)*phi(i,j-1)+ &
939 0 : cof(i,j,8)*phi(i+1,j-1))
940 : end do
941 0 : do i=2,nx
942 0 : phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
943 : end do
944 0 : phi(nx,j) = phi(nx,j)/tx(nx,j,2)
945 0 : do ib=2,nx
946 0 : i = nx-ib+1
947 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
948 : end do
949 : end do
950 : else
951 : !
952 : ! x direction periodic
953 : !
954 91656 : do j=1,ny
955 83277 : sum(j) = 0.0_r8
956 83277 : phi(0,j) = phi(nx-1,j)
957 91656 : phi(nx+1,j) = phi(2,j)
958 : end do
959 : !
960 : ! sweep odd lines forward and back
961 : !
962 : !$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib)
963 52478 : do j=1,ny,2
964 1450479 : do i=1,nx-1
965 5625520 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
966 : cof(i,j,3)*phi(i,j+1)+ &
967 1406380 : cof(i,j,4)*phi(i-1,j+1)+ &
968 1406380 : cof(i,j,6)*phi(i-1,j-1)+ &
969 : cof(i,j,7)*phi(i,j-1)+ &
970 9888759 : cof(i,j,8)*phi(i+1,j-1))
971 : end do
972 : !
973 : ! forward sweep
974 : !
975 1362281 : do i=2,nx-2
976 1362281 : phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
977 : end do
978 1406380 : do i=1,nx-2
979 1406380 : sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
980 : end do
981 44099 : phi(nx-1,j) = phi(nx-1,j)-sum(j)
982 : !
983 : ! backward sweep
984 : !
985 44099 : phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
986 88198 : phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ &
987 88198 : tx(nx-2,j,2)
988 1370660 : do ib=4,nx
989 1318182 : i = nx-ib+1
990 2636364 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* &
991 3998645 : phi(nx-1,j))/tx(i,j,2)
992 : end do
993 : end do
994 : !
995 : ! set periodic and virtual points for j odd
996 : !
997 8379 : do j=1,ny,2
998 44099 : phi(nx,j) = phi(1,j)
999 44099 : phi(0,j) = phi(nx-1,j)
1000 44099 : phi(nx+1,j) = phi(2,j)
1001 : end do
1002 : !
1003 : ! sweep even j lines
1004 : !
1005 : !$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib)
1006 47557 : do j=2,ny,2
1007 1338018 : do i=1,nx-1
1008 5195360 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
1009 : cof(i,j,3)*phi(i,j+1)+ &
1010 1298840 : cof(i,j,4)*phi(i-1,j+1)+ &
1011 1298840 : cof(i,j,6)*phi(i-1,j-1)+ &
1012 : cof(i,j,7)*phi(i,j-1)+ &
1013 9131058 : cof(i,j,8)*phi(i+1,j-1))
1014 : end do
1015 : !
1016 : ! forward sweep
1017 : !
1018 1259662 : do i=2,nx-2
1019 1259662 : phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
1020 : end do
1021 1298840 : do i=1,nx-2
1022 1298840 : sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
1023 : end do
1024 39178 : phi(nx-1,j) = phi(nx-1,j)-sum(j)
1025 : !
1026 : ! backward sweep
1027 : !
1028 39178 : phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
1029 78356 : phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ &
1030 78356 : tx(nx-2,j,2)
1031 1268041 : do ib=4,nx
1032 1220484 : i = nx-ib+1
1033 2440968 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* &
1034 3700630 : phi(nx-1,j))/tx(i,j,2)
1035 : end do
1036 : end do
1037 : !
1038 : ! set periodic and virtual points for j even
1039 : !
1040 8379 : do j=2,ny,2
1041 39178 : phi(nx,j) = phi(1,j)
1042 39178 : phi(0,j) = phi(nx-1,j)
1043 39178 : phi(nx+1,j) = phi(2,j)
1044 : end do
1045 : end if
1046 : !
1047 : ! set periodic y virtual boundaries if necessary
1048 : !
1049 8379 : if (nyc.eq.0) then
1050 0 : do i=1,nx
1051 0 : phi(i,0) = phi(i,ny-1)
1052 0 : phi(i,ny+1) = phi(i,2)
1053 : end do
1054 : end if
1055 : return
1056 0 : end subroutine slxmd2cr
1057 : !-----------------------------------------------------------------------
1058 8379 : subroutine slymd2cr(nx,ny,phi,cof,ty,sum)
1059 :
1060 : integer nx,ny,i,j,jb
1061 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
1062 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
1063 : kcycle,iprer,ipost,intpol,kps
1064 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
1065 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
1066 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1067 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10),ty(ny,nx,*),sum(nx)
1068 : !
1069 : ! replace line y with point gauss-seidel if
1070 : ! y direction is periodic and ny = 3
1071 : !
1072 8379 : if (nyc .eq. 0 .and. ny .eq. 3) then
1073 0 : call relmd2crp(nx,ny,phi,cof)
1074 0 : return
1075 : end if
1076 : !
1077 : ! set periodic and virtual x boundaries if necessary
1078 : !
1079 8379 : if (nxa.eq.0) then
1080 91656 : do j=1,ny
1081 83277 : phi(0,j) = phi(nx-1,j)
1082 83277 : phi(nx,j) = phi(1,j)
1083 91656 : phi(nx+1,j) = phi(2,j)
1084 : end do
1085 : end if
1086 :
1087 8379 : if (nyc.ne.0) then
1088 : !
1089 : ! y direction not periodic
1090 : !
1091 : !$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1092 77444 : do i=1,nx,2
1093 1498036 : do j=1,ny
1094 4286913 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1095 1428971 : cof(i,j,2)*phi(i+1,j+1)+ &
1096 1428971 : cof(i,j,4)*phi(i-1,j+1)+ &
1097 : cof(i,j,5)*phi(i-1,j)+ &
1098 1428971 : cof(i,j,6)*phi(i-1,j-1)+ &
1099 10071862 : cof(i,j,8)*phi(i+1,j-1))
1100 : end do
1101 : !
1102 : ! forward sweep thru odd x lines
1103 : !
1104 1428971 : do j=2,ny
1105 1428971 : phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
1106 : end do
1107 : !
1108 : ! backward sweep
1109 : !
1110 69065 : phi(i,ny) = phi(i,ny)/ty(ny,i,2)
1111 1437350 : do jb=2,ny
1112 1359906 : j = ny-jb+1
1113 1428971 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
1114 : end do
1115 : end do
1116 : !
1117 : ! forward sweep even x lines
1118 : !
1119 : !$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1120 72523 : do i=2,nx,2
1121 1423670 : do j=1,ny
1122 4078578 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1123 1359526 : cof(i,j,2)*phi(i+1,j+1)+ &
1124 1359526 : cof(i,j,4)*phi(i-1,j+1)+ &
1125 : cof(i,j,5)*phi(i-1,j)+ &
1126 1359526 : cof(i,j,6)*phi(i-1,j-1)+ &
1127 9580826 : cof(i,j,8)*phi(i+1,j-1))
1128 : end do
1129 1359526 : do j=2,ny
1130 1359526 : phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
1131 : end do
1132 : !
1133 : ! backward sweep
1134 : !
1135 64144 : phi(i,ny) = phi(i,ny)/ty(ny,i,2)
1136 1367905 : do jb=2,ny
1137 1295382 : j = ny-jb+1
1138 1359526 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
1139 : end do
1140 : end do
1141 : else
1142 : !
1143 : ! y direction periodic
1144 : !
1145 0 : do i=1,nx
1146 0 : sum(i) = 0.0_r8
1147 0 : phi(i,0) = phi(i,ny-1)
1148 0 : phi(i,ny) = phi(i,1)
1149 0 : phi(i,ny+1) = phi(i,2)
1150 : end do
1151 : !
1152 : ! forward sweep odd x lines
1153 : !
1154 : !$OMP PARALLEL DO SHARED(sum,cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1155 0 : do i=1,nx,2
1156 0 : do j=1,ny-1
1157 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1158 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1159 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1160 : cof(i,j,5)*phi(i-1,j)+ &
1161 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1162 0 : cof(i,j,8)*phi(i+1,j-1))
1163 : end do
1164 0 : do j=2,ny-2
1165 0 : phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
1166 : end do
1167 0 : do j=1,ny-2
1168 0 : sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
1169 : end do
1170 0 : phi(i,ny-1) = phi(i,ny-1)-sum(i)
1171 : !
1172 : ! backward sweep
1173 : !
1174 0 : phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
1175 0 : phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ &
1176 0 : ty(ny-2,i,2)
1177 0 : do jb=4,ny
1178 0 : j = ny-jb+1
1179 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* &
1180 0 : phi(i,ny-1))/ty(j,i,2)
1181 : end do
1182 : end do
1183 : !
1184 : ! set odd periodic and virtual y boundaries
1185 : !
1186 0 : do i=1,nx,2
1187 0 : phi(i,0) = phi(i,ny-1)
1188 0 : phi(i,ny) = phi(i,1)
1189 0 : phi(i,ny+1) = phi(i,2)
1190 : end do
1191 : !
1192 : ! forward sweep even x lines
1193 : !
1194 : !$OMP PARALLEL DO SHARED(sum,cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1195 0 : do i=2,nx,2
1196 0 : do j=1,ny-1
1197 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1198 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1199 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1200 : cof(i,j,5)*phi(i-1,j)+ &
1201 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1202 0 : cof(i,j,8)*phi(i+1,j-1))
1203 : end do
1204 0 : do j=2,ny-2
1205 0 : phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
1206 : end do
1207 0 : do j=1,ny-2
1208 0 : sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
1209 : end do
1210 0 : phi(i,ny-1) = phi(i,ny-1)-sum(i)
1211 : !
1212 : ! backward sweep
1213 : !
1214 0 : phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
1215 0 : phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ &
1216 0 : ty(ny-2,i,2)
1217 0 : do jb=4,ny
1218 0 : j = ny-jb+1
1219 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* &
1220 0 : phi(i,ny-1))/ty(j,i,2)
1221 : end do
1222 : end do
1223 : !
1224 : ! set even periodic and virtual y boundaries
1225 : !
1226 0 : do i=2,nx,2
1227 0 : phi(i,0) = phi(i,ny-1)
1228 0 : phi(i,ny) = phi(i,1)
1229 0 : phi(i,ny+1) = phi(i,2)
1230 : end do
1231 : end if
1232 : !
1233 : ! set periodic and virtual x boundaries if necessary
1234 : !
1235 8379 : if (nxa.eq.0) then
1236 91656 : do j=1,ny
1237 83277 : phi(0,j) = phi(nx-1,j)
1238 91656 : phi(nx+1,j) = phi(2,j)
1239 : end do
1240 : end if
1241 :
1242 : return
1243 16758 : end subroutine slymd2cr
1244 : !-----------------------------------------------------------------------
1245 : end module edyn_mud
|
