1c
2c     file mud3ln.f
3c     modification of mud3ln.f to include open MP 6/99
4c
5c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6c  .                                                             .
7c  .                  copyright (c) 1999 by UCAR                 .
8c  .                                                             .
9c  .       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH       .
10c  .                                                             .
11c  .                      all rights reserved                    .
12c  .                                                             .
13c  .                                                             .
14c  .                      MUDPACK version 5.0                    .
15c  .                                                             .
16c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17c
18c ... author and specialist
19c
20c          John C. Adams (National Center for Atmospheric Research)
21c          email: johnad@ucar.edu, phone: 303-497-1213
22
23c ... For MUDPACK 5.0 information, visit the website:
24c     (http://www.scd.ucar.edu/css/software/mudpack)
25c
26c ... purpose
27c
28c     mud3ln.f contains subroutines for line relaxation in the x and y
29c     and z direction.  This file must be loaded with any of the real
30c     3-d mudpack solvers except mud3sp.
31c
32      subroutine slxmd3(nx,ny,nz,phi,cof,tx,ssm,nxa,nyc,nze)
33c
34c     x line relaxation thru red and then black points in the
35c     (y,z) plane for periodic or nonperiodic x b.c.
36c
37      implicit none
38      integer nx,ny,nz,i,ib,j,k
39      integer nxa,nyc,nze,nper
40      double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8)
41      double precision tx(nx,ny,nz,*), ssm(ny,nz)
42c
43c     set periodic indicator
44c
45      nper =  nxa*nyc*nze
46c
47c     set periodic virtual boundaries as necessary
48c
49      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
50      if (nxa.ne.0) then
51c
52c     x direction not periodic
53c     first solve for x lines thru red points in (y,z) plane
54c
55!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(phi,cof,tx,nx,ny,nz)
56      do k=1,nz,2
57        do j=1,ny,2
58          do i=1,nx
59            phi(i,j,k) = cof(i,j,k,8) - (
60     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
61     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
62          end do
63        end do
64c
65c     forward sweep
66c
67        do i=2,nx
68          do j=1,ny,2
69            phi(i,j,k) = phi(i,j,k)-tx(i-1,j,k,1)*phi(i-1,j,k)
70          end do
71        end do
72c
73c     backward sweep
74c
75        do j=1,ny,2
76          phi(nx,j,k) = phi(nx,j,k)/tx(nx,j,k,2)
77        end do
78        do ib=2,nx
79          i = nx-ib+1
80          do j=1,ny,2
81            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k))
82     +                     /tx(i,j,k,2)
83          end do
84        end do
85c
86c     end of k odd loop
87c
88      end do
89c
90!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(phi,cof,tx,nx,ny,nz)
91      do k=2,nz,2
92        do j=2,ny,2
93          do i=1,nx
94            phi(i,j,k) = cof(i,j,k,8) - (
95     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
96     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
97          end do
98        end do
99        do i=2,nx
100          do j=2,ny,2
101            phi(i,j,k) = phi(i,j,k)-tx(i-1,j,k,1)*phi(i-1,j,k)
102          end do
103        end do
104        do j=2,ny,2
105          phi(nx,j,k) = phi(nx,j,k)/tx(nx,j,k,2)
106        end do
107        do ib=2,nx
108          i = nx-ib+1
109          do j=2,ny,2
110            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k))
111     +                     /tx(i,j,k,2)
112          end do
113        end do
114c
115c     end of k even loop
116c
117      end do
118
119      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
120c
121c     solve for x lines thru black points in (y,z) plane
122c
123!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(phi,cof,tx,nx,ny,nz)
124      do k=1,nz,2
125        do j=2,ny,2
126          do i=1,nx
127            phi(i,j,k) = cof(i,j,k,8) - (
128     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
129     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
130          end do
131        end do
132        do i=2,nx
133          do j=2,ny,2
134            phi(i,j,k) = phi(i,j,k)-tx(i-1,j,k,1)*phi(i-1,j,k)
135          end do
136        end do
137        do j=2,ny,2
138          phi(nx,j,k) = phi(nx,j,k)/tx(nx,j,k,2)
139        end do
140        do ib=2,nx
141          i = nx-ib+1
142          do j=2,ny,2
143            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k))
144     +                     /tx(i,j,k,2)
145          end do
146        end do
147      end do
148
149c
150c
151!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(phi,cof,tx,nx,ny,nz)
152      do k=2,nz,2
153        do j=1,ny,2
154          do i=1,nx
155            phi(i,j,k) = cof(i,j,k,8) - (
156     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
157     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
158          end do
159        end do
160        do i=2,nx
161          do j=1,ny,2
162            phi(i,j,k) = phi(i,j,k)-tx(i-1,j,k,1)*phi(i-1,j,k)
163          end do
164        end do
165        do j=1,ny,2
166          phi(nx,j,k) = phi(nx,j,k)/tx(nx,j,k,2)
167        end do
168        do ib=2,nx
169          i = nx-ib+1
170          do j=1,ny,2
171            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k))
172     +                     /tx(i,j,k,2)
173          end do
174        end do
175      end do
176      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
177      return
178      else
179c
180c     x direction periodic
181c
182      do k=1,nz
183        do j=1,ny
184          ssm(j,k) = 0.d0
185        end do
186      end do
187c
188c      sweep x lines thru red (y,z) forward and back
189c
190!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(ssm,phi,cof,tx,nx,ny,nz)
191      do k=1,nz,2
192        do j=1,ny,2
193          do i=1,nx-1
194            phi(i,j,k) = cof(i,j,k,8) - (
195     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
196     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
197          end do
198        end do
199c
200c     forward sweep
201c
202        do i=2,nx-2
203          do j=1,ny,2
204            phi(i,j,k) = phi(i,j,k)-tx(i,j,k,1)*phi(i-1,j,k)
205          end do
206        end do
207        do i=1,nx-2
208          do j=1,ny,2
209            ssm(j,k) = ssm(j,k)+tx(i,j,k,5)*phi(i,j,k)
210          end do
211        end do
212        do j=1,ny,2
213          phi(nx-1,j,k) = phi(nx-1,j,k)-ssm(j,k)
214        end do
215c
216c     backward sweep
217c
218        do j=1,ny,2
219          phi(nx-1,j,k) = phi(nx-1,j,k)/tx(nx-1,j,k,2)
220          phi(nx-2,j,k) = (phi(nx-2,j,k)-tx(nx-2,j,k,4)*phi(nx-1,j,k))
221     +                      /tx(nx-2,j,k,2)
222        end do
223        do ib=4,nx
224          i = nx-ib+1
225          do j=1,ny,2
226            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k)-
227     +                      tx(i,j,k,4)*phi(nx-1,j,k))/tx(i,j,k,2)
228          end do
229        end do
230      end do
231c
232c     set periodic virtual boundaries as necessary
233c
234      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
235c
236c     forward even-even
237c
238!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(ssm,phi,cof,tx,nx,ny,nz)
239      do k=2,nz,2
240        do j=2,ny,2
241          do i=1,nx-1
242            phi(i,j,k) = cof(i,j,k,8) - (
243     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
244     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
245          end do
246        end do
247        do i=2,nx-2
248          do j=2,ny,2
249            phi(i,j,k) = phi(i,j,k)-tx(i,j,k,1)*phi(i-1,j,k)
250          end do
251        end do
252        do i=1,nx-2
253          do j=2,ny,2
254            ssm(j,k) = ssm(j,k)+tx(i,j,k,5)*phi(i,j,k)
255          end do
256        end do
257        do j=2,ny,2
258          phi(nx-1,j,k) = phi(nx-1,j,k)-ssm(j,k)
259        end do
260        do j=2,ny,2
261          phi(nx-1,j,k) = phi(nx-1,j,k)/tx(nx-1,j,k,2)
262          phi(nx-2,j,k) = (phi(nx-2,j,k)-tx(nx-2,j,k,4)*phi(nx-1,j,k))/
263     +                     tx(nx-2,j,k,2)
264        end do
265        do ib=4,nx
266          i = nx-ib+1
267          do j=2,ny,2
268            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k)-
269     +                      tx(i,j,k,4)*phi(nx-1,j,k))/tx(i,j,k,2)
270          end do
271        end do
272      end do
273      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
274c
275c     now solve x lines thru black points in (y,z) plane
276c
277!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(ssm,phi,cof,tx,nx,ny,nz)
278      do k=1,nz,2
279        do j=2,ny,2
280          do i=1,nx-1
281            phi(i,j,k) = cof(i,j,k,8) - (
282     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
283     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
284          end do
285        end do
286        do i=2,nx-2
287          do j=2,ny,2
288            phi(i,j,k) = phi(i,j,k)-tx(i,j,k,1)*phi(i-1,j,k)
289          end do
290        end do
291        do i=1,nx-2
292          do j=2,ny,2
293            ssm(j,k) = ssm(j,k)+tx(i,j,k,5)*phi(i,j,k)
294          end do
295        end do
296        do j=2,ny,2
297          phi(nx-1,j,k) = phi(nx-1,j,k)-ssm(j,k)
298        end do
299        do j=2,ny,2
300          phi(nx-1,j,k) = phi(nx-1,j,k)/tx(nx-1,j,k,2)
301          phi(nx-2,j,k) = (phi(nx-2,j,k)-tx(nx-2,j,k,4)*phi(nx-1,j,k))/
302     +                     tx(nx-2,j,k,2)
303        end do
304        do ib=4,nx
305          i = nx-ib+1
306          do j=2,ny,2
307            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k)-
308     +                      tx(i,j,k,4)*phi(nx-1,j,k))/tx(i,j,k,2)
309          end do
310        end do
311      end do
312      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
313c
314c     forward even-odd
315c
316!$OMP PARALLEL DO PRIVATE(i,ib,j,k), SHARED(ssm,phi,cof,tx,nx,ny,nz)
317      do k=2,nz,2
318        do j=1,ny,2
319          do i=1,nx-1
320            phi(i,j,k) = cof(i,j,k,8) - (
321     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k) +
322     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
323          end do
324        end do
325        do i=2,nx-2
326          do j=1,ny,2
327            phi(i,j,k) = phi(i,j,k)-tx(i,j,k,1)*phi(i-1,j,k)
328          end do
329        end do
330        do i=1,nx-2
331          do j=1,ny,2
332            ssm(j,k) = ssm(j,k)+tx(i,j,k,5)*phi(i,j,k)
333          end do
334        end do
335        do j=1,ny,2
336          phi(nx-1,j,k) = phi(nx-1,j,k)-ssm(j,k)
337        end do
338        do j=1,ny,2
339          phi(nx-1,j,k) = phi(nx-1,j,k)/tx(nx-1,j,k,2)
340          phi(nx-2,j,k) = (phi(nx-2,j,k)-tx(nx-2,j,k,4)*phi(nx-1,j,k))/
341     +                     tx(nx-2,j,k,2)
342        end do
343        do ib=4,nx
344          i = nx-ib+1
345          do j=1,ny,2
346            phi(i,j,k) = (phi(i,j,k)-tx(i,j,k,3)*phi(i+1,j,k)-
347     +                      tx(i,j,k,4)*phi(nx-1,j,k))/tx(i,j,k,2)
348          end do
349        end do
350      end do
351      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
352      return
353      end if
354      end
355
356      subroutine slymd3(nx,ny,nz,phi,cof,ty,ssm,nxa,nyc,nze)
357c
358c     y line relaxation thru red and then black points in the
359c     (x,z) plane for periodic or nonperiodic y b.c.
360c
361      implicit none
362      integer nx,ny,nz,i,j,jb,k
363      integer nxa,nyc,nze,nper
364      double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8)
365      double precision ty(ny,nx,nz,*), ssm(nx,nz)
366c
367c     set periodic indicator
368c
369      nper =  nxa*nyc*nze
370c
371c     set periodic virtual boundaries as necessary
372c
373      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
374      if (nyc.ne.0) then
375c
376c     y direction not periodic
377c     first solve for y lines thru red points in (x,z) plane
378c
379!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(phi,cof,ty,nx,ny,nz)
380      do k=1,nz,2
381        do i=1,nx,2
382          do j=1,ny
383            phi(i,j,k) = cof(i,j,k,8) - (
384     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
385     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
386          end do
387        end do
388c
389c     forward sweep
390c
391        do j=2,ny
392          do i=1,nx,2
393            phi(i,j,k) = phi(i,j,k)-ty(j-1,i,k,1)*phi(i,j-1,k)
394          end do
395        end do
396c
397c     backward sweep
398c
399        do i=1,nx,2
400          phi(i,ny,k) = phi(i,ny,k)/ty(ny,i,k,2)
401        end do
402        do jb=2,ny
403          j = ny-jb+1
404          do i=1,nx,2
405            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k))
406     +                     /ty(j,i,k,2)
407          end do
408        end do
409      end do
410!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(phi,cof,ty,nx,ny,nz)
411      do k=2,nz,2
412        do i=2,nx,2
413          do j=1,ny
414            phi(i,j,k) = cof(i,j,k,8) - (
415     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
416     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
417          end do
418        end do
419        do j=2,ny
420          do i=2,nx,2
421            phi(i,j,k) = phi(i,j,k)-ty(j-1,i,k,1)*phi(i,j-1,k)
422          end do
423        end do
424        do i=2,nx,2
425          phi(i,ny,k) = phi(i,ny,k)/ty(ny,i,k,2)
426        end do
427        do jb=2,ny
428          j = ny-jb+1
429          do i=2,nx,2
430            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k))
431     +                     /ty(j,i,k,2)
432          end do
433        end do
434      end do
435      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
436c
437c     solve for x lines thru black points in (y,z) plane
438c
439c
440!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(phi,cof,ty,nx,ny,nz)
441      do k=1,nz,2
442        do i=2,nx,2
443          do j=1,ny
444            phi(i,j,k) = cof(i,j,k,8) - (
445     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
446     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
447          end do
448        end do
449        do j=2,ny
450          do i=2,nx,2
451            phi(i,j,k) = phi(i,j,k)-ty(j-1,i,k,1)*phi(i,j-1,k)
452          end do
453        end do
454        do i=2,nx,2
455          phi(i,ny,k) = phi(i,ny,k)/ty(ny,i,k,2)
456        end do
457        do jb=2,ny
458          j = ny-jb+1
459          do i=2,nx,2
460            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k))
461     +                     /ty(j,i,k,2)
462          end do
463        end do
464      end do
465      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
466!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(phi,cof,ty,nx,ny,nz)
467      do k=2,nz,2
468        do i=1,nx,2
469          do j=1,ny
470            phi(i,j,k) = cof(i,j,k,8) - (
471     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
472     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
473          end do
474        end do
475        do j=2,ny
476          do i=1,nx,2
477            phi(i,j,k) = phi(i,j,k)-ty(j-1,i,k,1)*phi(i,j-1,k)
478          end do
479        end do
480        do i=1,nx,2
481          phi(i,ny,k) = phi(i,ny,k)/ty(ny,i,k,2)
482        end do
483        do jb=2,ny
484          j = ny-jb+1
485          do i=1,nx,2
486            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k))
487     +                     /ty(j,i,k,2)
488          end do
489        end do
490      end do
491      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
492      return
493      else
494c
495c     y direction periodic
496c
497      do k=1,nz
498        do i=1,nx
499          ssm(i,k) = 0.d0
500        end do
501      end do
502c
503c      sweep y lines thru red (x,z) forward and back
504c
505!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(ssm,phi,cof,ty,nx,ny,nz)
506      do k=1,nz,2
507        do i=1,nx,2
508          do j=1,ny-1
509            phi(i,j,k) = cof(i,j,k,8) - (
510     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
511     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
512          end do
513        end do
514c
515c     forward sweep
516c
517        do j=2,ny-2
518          do i=1,nx,2
519            phi(i,j,k) = phi(i,j,k)-ty(j,i,k,1)*phi(i,j-1,k)
520          end do
521        end do
522        do j=1,ny-2
523          do i=1,nx,2
524            ssm(i,k) = ssm(i,k)+ty(j,i,k,5)*phi(i,j,k)
525          end do
526        end do
527        do i=1,nx,2
528          phi(i,ny-1,k) = phi(i,ny-1,k)-ssm(i,k)
529        end do
530c
531c     backward sweep
532c
533        do i=1,nx,2
534          phi(i,ny-1,k) = phi(i,ny-1,k)/ty(ny-1,i,k,2)
535          phi(i,ny-2,k) = (phi(i,ny-2,k)-ty(ny-2,i,k,4)*phi(i,ny-1,k))
536     +                      /ty(ny-2,i,k,2)
537        end do
538        do jb=4,ny
539          j = ny-jb+1
540          do i=1,nx,2
541            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k)-
542     +                      ty(j,i,k,4)*phi(i,ny-1,k))/ty(j,i,k,2)
543          end do
544        end do
545      end do
546c
547c     set periodic virtual boundaries as necessary
548c
549      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
550c
551c     forward even-even
552c
553!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(ssm,phi,cof,ty,nx,ny,nz)
554      do k=2,nz,2
555        do i=2,nx,2
556          do j=1,ny-1
557            phi(i,j,k) = cof(i,j,k,8) - (
558     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
559     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
560          end do
561        end do
562        do j=2,ny-2
563          do i=2,nx,2
564            phi(i,j,k) = phi(i,j,k)-ty(j,i,k,1)*phi(i,j-1,k)
565          end do
566        end do
567        do j=1,ny-2
568          do i=2,nx,2
569            ssm(i,k) = ssm(i,k)+ty(j,i,k,5)*phi(i,j,k)
570          end do
571        end do
572        do i=2,nx,2
573          phi(i,ny-1,k) = phi(i,ny-1,k)-ssm(i,k)
574        end do
575        do i=2,nx,2
576          phi(i,ny-1,k) = phi(i,ny-1,k)/ty(ny-1,i,k,2)
577          phi(i,ny-2,k) = (phi(i,ny-2,k)-ty(ny-2,i,k,4)*phi(i,ny-1,k))
578     +                      /ty(ny-2,i,k,2)
579        end do
580        do jb=4,ny
581          j = ny-jb+1
582          do i=2,nx,2
583            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k)-
584     +                      ty(j,i,k,4)*phi(i,ny-1,k))/ty(j,i,k,2)
585          end do
586        end do
587      end do
588      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
589c
590c     now solve x lines thru black points in (y,z) plane
591c
592!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(ssm,phi,cof,ty,nx,ny,nz)
593      do k=1,nz,2
594        do i=2,nx,2
595          do j=1,ny-1
596            phi(i,j,k) = cof(i,j,k,8) - (
597     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
598     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
599          end do
600        end do
601        do j=2,ny-2
602          do i=2,nx,2
603            phi(i,j,k) = phi(i,j,k)-ty(j,i,k,1)*phi(i,j-1,k)
604          end do
605        end do
606        do j=1,ny-2
607          do i=2,nx,2
608            ssm(i,k) = ssm(i,k)+ty(j,i,k,5)*phi(i,j,k)
609          end do
610        end do
611        do i=2,nx,2
612          phi(i,ny-1,k) = phi(i,ny-1,k)-ssm(i,k)
613        end do
614        do i=2,nx,2
615          phi(i,ny-1,k) = phi(i,ny-1,k)/ty(ny-1,i,k,2)
616          phi(i,ny-2,k) = (phi(i,ny-2,k)-ty(ny-2,i,k,4)*phi(i,ny-1,k))
617     +                      /ty(ny-2,i,k,2)
618        end do
619        do jb=4,ny
620          j = ny-jb+1
621          do i=2,nx,2
622            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k)-
623     +                      ty(j,i,k,4)*phi(i,ny-1,k))/ty(j,i,k,2)
624          end do
625        end do
626      end do
627c
628c     set periodic virtual boundaries as necessary
629c
630      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
631c
632c     forward even-even
633c
634!$OMP PARALLEL DO PRIVATE(i,j,jb,k), SHARED(ssm,phi,cof,ty,nx,ny,nz)
635      do k=2,nz,2
636        do i=1,nx,2
637          do j=1,ny-1
638            phi(i,j,k) = cof(i,j,k,8) - (
639     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
640     +        cof(i,j,k,5)*phi(i,j,k-1)+cof(i,j,k,6)*phi(i,j,k+1))
641          end do
642        end do
643        do j=2,ny-2
644          do i=1,nx,2
645            phi(i,j,k) = phi(i,j,k)-ty(j,i,k,1)*phi(i,j-1,k)
646          end do
647        end do
648        do j=1,ny-2
649          do i=1,nx,2
650            ssm(i,k) = ssm(i,k)+ty(j,i,k,5)*phi(i,j,k)
651          end do
652        end do
653        do i=1,nx,2
654          phi(i,ny-1,k) = phi(i,ny-1,k)-ssm(i,k)
655        end do
656        do i=1,nx,2
657          phi(i,ny-1,k) = phi(i,ny-1,k)/ty(ny-1,i,k,2)
658          phi(i,ny-2,k) = (phi(i,ny-2,k)-ty(ny-2,i,k,4)*phi(i,ny-1,k))
659     +                      /ty(ny-2,i,k,2)
660        end do
661        do jb=4,ny
662          j = ny-jb+1
663          do i=1,nx,2
664            phi(i,j,k) = (phi(i,j,k)-ty(j,i,k,3)*phi(i,j+1,k)-
665     +                      ty(j,i,k,4)*phi(i,ny-1,k))/ty(j,i,k,2)
666          end do
667        end do
668      end do
669      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
670      return
671      end if
672      return
673      end
674
675      subroutine slzmd3(nx,ny,nz,phi,cof,tz,ssm,nxa,nyc,nze)
676c
677c     z line relaxation thru red and then black points in the
678c     (x,y) plane for periodic or nonperiodic z b.c.
679c
680      implicit none
681      integer nx,ny,nz,i,j,k,kb
682      integer nxa,nyc,nze,nper
683      double precision phi(0:nx+1,0:ny+1,0:nz+1),cof(nx,ny,nz,8)
684      double precision tz(nz,nx,ny,*), ssm(nx,ny)
685c
686c     set periodic indicator
687c
688      nper =  nxa*nyc*nze
689c
690c     set periodic virtual boundaries as necessary
691c
692      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
693      if (nze.ne.0) then
694c
695c     z direction not periodic
696c     first solve for z lines thru red points in (x,y) plane
697c
698!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(phi,cof,tz,nx,ny,nz)
699      do j=1,ny,2
700        do i=1,nx,2
701          do k=1,nz
702            phi(i,j,k) = cof(i,j,k,8) - (
703     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
704     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
705          end do
706        end do
707c
708c     forward sweep
709c
710        do k=2,nz
711          do i=1,nx,2
712            phi(i,j,k) = phi(i,j,k)-tz(k-1,i,j,1)*phi(i,j,k-1)
713          end do
714        end do
715c
716c     backward sweep
717c
718        do i=1,nx,2
719          phi(i,j,nz) = phi(i,j,nz)/tz(nz,i,j,2)
720        end do
721        do kb=2,nz
722          k = nz-kb+1
723          do i=1,nx,2
724            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1))
725     +                     /tz(k,i,j,2)
726          end do
727        end do
728      end do
729      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
730!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(phi,cof,tz,nx,ny,nz)
731      do j=2,ny,2
732        do i=2,nx,2
733          do k=1,nz
734            phi(i,j,k) = cof(i,j,k,8) - (
735     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
736     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
737          end do
738        end do
739        do k=2,nz
740          do i=2,nx,2
741            phi(i,j,k) = phi(i,j,k)-tz(k-1,i,j,1)*phi(i,j,k-1)
742          end do
743        end do
744        do i=2,nx,2
745          phi(i,j,nz) = phi(i,j,nz)/tz(nz,i,j,2)
746        end do
747        do kb=2,nz
748          k = nz-kb+1
749          do i=2,nx,2
750            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1))
751     +                     /tz(k,i,j,2)
752          end do
753        end do
754      end do
755      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
756c
757c     solve for z lines thru black points in (x,y)plane
758c
759!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(phi,cof,tz,nx,ny,nz)
760      do j=1,ny,2
761        do i=2,nx,2
762          do k=1,nz
763            phi(i,j,k) = cof(i,j,k,8) - (
764     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
765     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
766          end do
767        end do
768        do k=2,nz
769          do i=2,nx,2
770            phi(i,j,k) = phi(i,j,k)-tz(k-1,i,j,1)*phi(i,j,k-1)
771          end do
772        end do
773        do i=2,nx,2
774          phi(i,j,nz) = phi(i,j,nz)/tz(nz,i,j,2)
775        end do
776        do kb=2,nz
777          k = nz-kb+1
778          do i=2,nx,2
779            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1))
780     +                     /tz(k,i,j,2)
781          end do
782        end do
783      end do
784      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
785c
786!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(phi,cof,tz,nx,ny,nz)
787      do j=2,ny,2
788        do i=1,nx,2
789          do k=1,nz
790            phi(i,j,k) = cof(i,j,k,8) - (
791     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
792     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
793          end do
794        end do
795        do k=2,nz
796          do i=1,nx,2
797            phi(i,j,k) = phi(i,j,k)-tz(k-1,i,j,1)*phi(i,j,k-1)
798          end do
799        end do
800        do i=1,nx,2
801          phi(i,j,nz) = phi(i,j,nz)/tz(nz,i,j,2)
802        end do
803        do kb=2,nz
804          k = nz-kb+1
805          do i=1,nx,2
806            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1))
807     +                     /tz(k,i,j,2)
808          end do
809        end do
810      end do
811      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
812      return
813      else
814c
815c     z direction periodic
816c
817      do j=1,ny
818        do i=1,nx
819          ssm(i,j) = 0.d0
820        end do
821      end do
822c
823c      sweep z lines thru red (x,y) forward and back
824c
825!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(ssm,phi,cof,tz,nx,ny,nz)
826      do j=1,ny,2
827        do i=1,nx,2
828          do k=1,nz-1
829            phi(i,j,k) = cof(i,j,k,8) - (
830     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
831     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
832          end do
833        end do
834c
835c     forward sweep
836c
837        do k=2,nz-2
838          do i=1,nx,2
839            phi(i,j,k) = phi(i,j,k)-tz(k,i,j,1)*phi(i,j,k-1)
840          end do
841        end do
842        do k=1,nz-2
843          do i=1,nx,2
844            ssm(i,j) = ssm(i,j)+tz(k,i,j,5)*phi(i,j,k)
845          end do
846        end do
847        do i=1,nx,2
848          phi(i,j,nz-1) = phi(i,j,nz-1)-ssm(i,j)
849        end do
850c
851c     backward sweep
852c
853        do i=1,nx,2
854          phi(i,j,nz-1) = phi(i,j,nz-1)/tz(nz-1,i,j,2)
855          phi(i,j,nz-2) = (phi(i,j,nz-2)-tz(nz-2,i,j,4)*phi(i,j,nz-1))
856     +                      /tz(nz-2,i,j,2)
857        end do
858        do kb=4,nz
859          k = nz-kb+1
860          do i=1,nx,2
861            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1)-
862     +                      tz(k,i,j,4)*phi(i,j,nz-1))/tz(k,i,j,2)
863          end do
864        end do
865      end do
866c
867c     set periodic virtual boundaries as necessary
868c
869      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
870c
871c     sweep black z lines thru (x,y)
872c
873!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(ssm,phi,cof,tz,nx,ny,nz)
874      do j=2,ny,2
875        do i=2,nx,2
876          do k=1,nz-1
877            phi(i,j,k) = cof(i,j,k,8) - (
878     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
879     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
880          end do
881        end do
882        do k=2,nz-2
883          do i=2,nx,2
884            phi(i,j,k) = phi(i,j,k)-tz(k,i,j,1)*phi(i,j,k-1)
885          end do
886        end do
887        do k=1,nz-2
888          do i=2,nx,2
889            ssm(i,j) = ssm(i,j)+tz(k,i,j,5)*phi(i,j,k)
890          end do
891        end do
892        do i=2,nx,2
893          phi(i,j,nz-1) = phi(i,j,nz-1)-ssm(i,j)
894        end do
895        do i=2,nx,2
896          phi(i,j,nz-1) = phi(i,j,nz-1)/tz(nz-1,i,j,2)
897          phi(i,j,nz-2) = (phi(i,j,nz-2)-tz(nz-2,i,j,4)*phi(i,j,nz-1))
898     +                      /tz(nz-2,i,j,2)
899        end do
900        do kb=4,nz
901          k = nz-kb+1
902          do i=2,nx,2
903            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1)-
904     +                      tz(k,i,j,4)*phi(i,j,nz-1))/tz(k,i,j,2)
905          end do
906        end do
907      end do
908      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
909!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(ssm,phi,cof,tz,nx,ny,nz)
910      do j=1,ny,2
911        do i=2,nx,2
912          do k=1,nz-1
913            phi(i,j,k) = cof(i,j,k,8) - (
914     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
915     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
916          end do
917        end do
918        do k=2,nz-2
919          do i=2,nx,2
920            phi(i,j,k) = phi(i,j,k)-tz(k,i,j,1)*phi(i,j,k-1)
921          end do
922        end do
923        do k=1,nz-2
924          do i=2,nx,2
925            ssm(i,j) = ssm(i,j)+tz(k,i,j,5)*phi(i,j,k)
926          end do
927        end do
928        do i=2,nx,2
929          phi(i,j,nz-1) = phi(i,j,nz-1)-ssm(i,j)
930        end do
931        do i=2,nx,2
932          phi(i,j,nz-1) = phi(i,j,nz-1)/tz(nz-1,i,j,2)
933          phi(i,j,nz-2) = (phi(i,j,nz-2)-tz(nz-2,i,j,4)*phi(i,j,nz-1))
934     +                      /tz(nz-2,i,j,2)
935        end do
936        do kb=4,nz
937          k = nz-kb+1
938          do i=2,nx,2
939            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1)-
940     +                      tz(k,i,j,4)*phi(i,j,nz-1))/tz(k,i,j,2)
941          end do
942        end do
943      end do
944      if (nper.eq.0) call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
945c
946!$OMP PARALLEL DO PRIVATE(i,j,k,kb), SHARED(ssm,phi,cof,tz,nx,ny,nz)
947      do j=2,ny,2
948        do i=1,nx,2
949          do k=1,nz-1
950            phi(i,j,k) = cof(i,j,k,8) - (
951     +        cof(i,j,k,1)*phi(i-1,j,k)+cof(i,j,k,2)*phi(i+1,j,k) +
952     +        cof(i,j,k,3)*phi(i,j-1,k)+cof(i,j,k,4)*phi(i,j+1,k))
953          end do
954        end do
955        do k=2,nz-2
956          do i=1,nx,2
957            phi(i,j,k) = phi(i,j,k)-tz(k,i,j,1)*phi(i,j,k-1)
958          end do
959        end do
960        do k=1,nz-2
961          do i=1,nx,2
962            ssm(i,j) = ssm(i,j)+tz(k,i,j,5)*phi(i,j,k)
963          end do
964        end do
965        do i=1,nx,2
966          phi(i,j,nz-1) = phi(i,j,nz-1)-ssm(i,j)
967        end do
968        do i=1,nx,2
969          phi(i,j,nz-1) = phi(i,j,nz-1)/tz(nz-1,i,j,2)
970          phi(i,j,nz-2) = (phi(i,j,nz-2)-tz(nz-2,i,j,4)*phi(i,j,nz-1))
971     +                      /tz(nz-2,i,j,2)
972        end do
973        do kb=4,nz
974          k = nz-kb+1
975          do i=1,nx,2
976            phi(i,j,k) = (phi(i,j,k)-tz(k,i,j,3)*phi(i,j,k+1)-
977     +                      tz(k,i,j,4)*phi(i,j,nz-1))/tz(k,i,j,2)
978          end do
979        end do
980      end do
981      call per3vb(nx,ny,nz,phi,nxa,nyc,nze)
982      return
983      end if
984      return
985      end
986