Line data Source code
1 : module edyn_muh2cr
2 : use shr_kind_mod , only: r8 => shr_kind_r8
3 : use cam_abortutils, only: endrun
4 : use edyn_mudcom, only: prolon2, trsfc2, factri,factrp, sgfa, sgsl, transp
5 : use edyn_mudcom, only: swk2, cor2, transp, res2
6 :
7 : implicit none
8 :
9 : private
10 :
11 : public :: muh
12 :
13 : contains
14 : !-----------------------------------------------------------------------
15 0 : subroutine muh(pe,nlon,nlat,nlev,jntl)
16 : use edyn_solver_coefs, only: cee
17 : use edyn_params, only: pi
18 :
19 : implicit none
20 :
21 : integer,intent(in) :: nlon, nlat, nlev, jntl
22 : real(r8),intent(out) :: PE(nlon+1,*)
23 : !
24 : ! set grid size params
25 : !
26 : integer :: iixp, jjyq
27 : integer,parameter :: iiex = 1, jjey = 1
28 : integer :: nnx, nny
29 : !
30 : ! estimate work space for point relaxation (see muh2cr.d)
31 : !
32 : integer :: llwork
33 : integer :: iiwork
34 0 : real(r8), allocatable :: phi(:,:),rhs(:,:),work(:)
35 0 : integer, allocatable :: iwork(:)
36 : !
37 : ! put integer and floating point argument names in contiguous
38 : ! storage for labelling in vectors iprm,fprm
39 : !
40 : integer :: iprm(17),mgopt(4)
41 : real(r8) :: fprm(6)
42 : integer :: intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,&
43 : iguess,maxcy,method,nwork,lwrkqd,itero
44 : common/itmud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,&
45 : iguess,maxcy,method,nwork,lwrkqd,itero
46 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
47 : common/ftmud2cr/xa,xb,yc,yd,tolmax,relmax
48 : equivalence(intl,iprm)
49 : equivalence(xa,fprm)
50 : integer :: i,j,ierror
51 : integer, parameter :: maxcya = 1
52 : integer jj,jjj
53 :
54 0 : iixp = nlon
55 0 : jjyq = (nlat-1)/2
56 0 : nnx=iixp*2**(iiex-1)+1
57 0 : nny=jjyq*2**(jjey-1)+1
58 : llwork=(5*((nnx+2)*(nny+2)+18*nnx*nny)/3+ &
59 0 : (nnx+2)*(nny+2)+ (iixp+1)*(jjyq+1)*(2*iixp+3))
60 0 : iiwork=(iixp+1)*(jjyq+1)
61 :
62 0 : allocate(phi(nnx,nny),rhs(nnx,nny),work(llwork))
63 0 : allocate(iwork(iiwork))
64 :
65 : !
66 : ! SET INPUT INTEGER PARAMETERS
67 : !
68 0 : INTL = JNTL
69 : !
70 : ! set boundary condition flags
71 : !
72 0 : nxa = 0
73 0 : nxb = 0
74 0 : nyc = 2
75 0 : nyd = 1
76 : !
77 : ! set grid sizes from parameter statements
78 : !
79 0 : ixp = iixp
80 0 : jyq = jjyq
81 0 : iex = iiex
82 0 : jey = jjey
83 0 : nx = nnx
84 0 : ny = nny
85 : !
86 : ! set multigrid arguments (w(2,1) cycling with fully weighted
87 : ! residual restriction and cubic prolongation)
88 : !
89 0 : mgopt(1) = 2
90 0 : mgopt(2) = 2
91 0 : mgopt(3) = 2
92 0 : if (nlat<=97) then
93 0 : mgopt(4) = 3
94 : else
95 : ! 1 deg, changed to mgopt(4) = 1 per Astrid's suggestion
96 0 : mgopt(4) = 1
97 : end if
98 : !
99 : ! set for one cycle
100 : !
101 0 : maxcy = maxcya
102 : !
103 : ! set no initial guess forcing full multigrid cycling
104 : !
105 0 : iguess = 0
106 : !
107 : ! set work space length approximation from parameter statement
108 : !
109 0 : nwork = llwork
110 : !
111 : ! set line z relaxation
112 : !
113 0 : method = 3
114 : !
115 : ! set end points of solution rectangle in (x,y) space
116 : !
117 0 : xa = -pi
118 0 : xb = pi
119 0 : yc = 0.0_r8
120 0 : yd = 0.5_r8*pi
121 : !
122 : ! set error control flag
123 : !
124 0 : if (nlev>6) then
125 0 : tolmax = 0.05_r8
126 0 : else if (nlev>5) then
127 0 : tolmax = 0.03_r8
128 : else
129 0 : tolmax = 0.01_r8
130 : end if
131 : !
132 : ! set right hand side in rhs
133 : ! initialize phi to zero
134 : !
135 0 : do i=1,nx
136 0 : do j=1,ny
137 0 : RHS(I,J) = CEE(I+(J-1)*NX+9*NX*NY)
138 0 : phi(i,j) = 0.0_r8
139 : end do
140 : end do
141 : !
142 : ! set specified boundaries in phi
143 : !
144 0 : DO I=1,NX
145 0 : PHI(I,NY) = RHS(I,NY)/CEE(I+(NY-1)*NX+8*NX*NY)
146 : END DO
147 : !
148 : ! set specified boundaries in phi
149 : !
150 0 : DO I=1,NX
151 0 : PHI(I,NY) = RHS(I,NY)/CEE(I+(NY-1)*NX+8*NX*NY)
152 : END DO
153 :
154 : !
155 : ! intialization call
156 : !
157 0 : call muh2cr(iprm,fprm,work,iwork,rhs,phi,mgopt,ierror)
158 0 : if (ierror.gt.0) call endrun('muh call init muh2cr')
159 : !
160 : ! attempt solution
161 : !
162 0 : intl = 1
163 0 : call muh2cr(iprm,fprm,work,iwork,rhs,phi,mgopt,ierror)
164 0 : if (ierror.gt.0) call endrun('muh call solve muh2cr')
165 : !
166 : ! COPY PHI TO PE
167 : !
168 0 : DO J = 1,NY
169 0 : JJ = NY+J-1
170 0 : JJJ = NY+1-J
171 0 : DO I = 1,NX
172 0 : PE(I,JJ) = PHI(I,J)
173 0 : PE(I,JJJ) = PHI(I,J)
174 : END DO
175 : END DO
176 :
177 0 : deallocate( phi, rhs, work, iwork)
178 :
179 0 : end subroutine muh
180 : !-----------------------------------------------------------------------
181 0 : subroutine muh2cr(iparm,fparm,wk,iwk,rhs,phi,mgopt,ierror)
182 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
183 : implicit none
184 : integer iparm(17),mgopt(4),ierror,iwk(*)
185 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
186 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
187 : kcycle,iprer,ipost,intpol,kps
188 : real(r8) :: fparm(6),xa,xb,yc,yd,tolmax,relmax
189 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
190 : integer int,iw,k,kb,nx,ny,ic,itx,ity
191 : real(r8) :: wk(*),phi(*),rhs(*)
192 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
193 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
194 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
195 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
196 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),&
197 : nxk(50),nyk(50),isx,jsy
198 : integer ibeta,ialfa,izmat,idmat
199 : common/mh2cr/ibeta,ialfa,izmat,idmat
200 : data int / 0 /
201 : save int
202 0 : ierror = 1
203 0 : intl = iparm(1) ! set and check intl on all calls
204 0 : if (intl*(intl-1).ne.0) return
205 0 : if (int.eq.0) then
206 0 : int = 1
207 0 : if (intl.ne.0) return ! very first call is not intl=0
208 : end if
209 0 : ierror = 0
210 : !
211 : ! set arguments internally
212 : ! these will not be rechecked if intl=1!
213 : !
214 0 : nxa = iparm(2)
215 0 : nxb = iparm(3)
216 0 : nyc = iparm(4)
217 0 : nyd = iparm(5)
218 0 : ixp = iparm(6)
219 0 : jyq = iparm(7)
220 0 : iex = iparm(8)
221 0 : jey = iparm(9)
222 0 : ngrid = max0(iex,jey)
223 0 : nfx = iparm(10)
224 0 : nfy = iparm(11)
225 0 : iguess = iparm(12)
226 0 : maxcy = iparm(13)
227 0 : method = iparm(14)
228 0 : nwork = iparm(15)
229 0 : kcycle = mgopt(1)
230 0 : if (kcycle .eq. 0) then
231 : ! set defaults
232 0 : kcycle = 2
233 0 : iprer = 2
234 0 : ipost = 1
235 0 : intpol = 3
236 : else
237 0 : iprer = mgopt(2)
238 0 : ipost = mgopt(3)
239 0 : intpol = mgopt(4)
240 : end if
241 0 : xa = fparm(1)
242 0 : xb = fparm(2)
243 0 : yc = fparm(3)
244 0 : yd = fparm(4)
245 0 : tolmax = fparm(5)
246 0 : if (intl .eq. 0) then ! intialization call
247 : !
248 : ! check input arguments
249 : !
250 0 : ierror = 2 ! check boundary condition flags
251 0 : if (max0(nxa,nxb,nyc,nyd).gt.2) return
252 0 : if (min0(nxa,nxb,nyc,nyd).lt.0) return
253 0 : if (nxa.eq.0.and.nxb.ne.0) return
254 0 : if (nxa.ne.0.and.nxb.eq.0) return
255 0 : if (nyc.eq.0.and.nyd.ne.0) return
256 0 : if (nyc.ne.0.and.nyd.eq.0) return
257 0 : ierror = 3 ! check grid sizes
258 0 : if (ixp.lt.2) return
259 0 : if (jyq.lt.2) return
260 0 : ierror = 4
261 0 : ngrid = max0(iex,jey)
262 0 : if (iex.lt.1) return
263 0 : if (jey.lt.1) return
264 0 : if (ngrid.gt.50) return
265 0 : ierror = 5
266 0 : if (nfx.ne.ixp*2**(iex-1)+1) return
267 0 : if (nfy.ne.jyq*2**(jey-1)+1) return
268 0 : ierror = 6
269 0 : if (iguess*(iguess-1).ne.0) return
270 0 : ierror = 7
271 0 : if (maxcy.lt.1) return
272 0 : ierror = 8
273 0 : if (method.lt.0 .or. method.gt.3) return
274 0 : ierror = 9
275 : ! compute and test minimum work space
276 0 : isx = 0
277 0 : if (method.eq.1 .or. method.eq.3) then
278 0 : if (nxa.ne.0) isx = 3
279 0 : if (nxa.eq.0) isx = 5
280 : end if
281 0 : jsy = 0
282 0 : if (method.eq.2 .or. method.eq.3) then
283 0 : if (nyc.ne.0) jsy = 3
284 0 : if (nyc.eq.0) jsy = 5
285 : end if
286 0 : kps = 1
287 0 : do k=1,ngrid
288 : ! set subgrid sizes
289 0 : nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1
290 0 : nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1
291 0 : nx = nxk(k)
292 0 : ny = nyk(k)
293 0 : kps = kps+(nx+2)*(ny+2)+nx*ny*(10+isx+jsy)
294 : end do
295 : !
296 : ! set pointers for direct at coarse grid
297 : !
298 0 : nx = ixp+1
299 0 : ny = jyq+1
300 0 : ibeta = kps+1
301 0 : if (nyc .eq. 0) then
302 0 : ialfa = ibeta + nx*nx*(ny-1)
303 0 : izmat = ialfa+nx*nx*(ny-1)
304 0 : idmat = izmat+nx*nx*(ny-2)
305 0 : kps = idmat+nx*nx*(ny-2)
306 : else
307 0 : ialfa = ibeta + nx*nx*ny
308 0 : kps = ialfa+nx*nx*ny
309 : end if
310 0 : iparm(16) = kps+(nfx+2)*(nfy+2) ! exact minimum work space
311 0 : lwork = iparm(16)
312 0 : if (lwork .gt. nwork) return
313 0 : ierror = 10 ! check solution region
314 0 : if (xb.le.xa .or. yd.le.yc) return
315 0 : ierror = 11
316 0 : if (tolmax .lt. 0.0_r8) return
317 0 : ierror = 12 ! multigrid parameters
318 0 : if (kcycle.lt.0) return
319 0 : if (min0(iprer,ipost).lt.1) return
320 0 : if ((intpol-1)*(intpol-3).ne.0) return
321 0 : if (max0(kcycle,iprer,ipost).gt.2) then
322 0 : ierror = -5 ! inefficient multigrid cycling
323 : end if
324 0 : if (ierror .gt. 0) ierror = 0 ! no fatal errors
325 : !
326 : ! set work space pointers and discretize pde at each grid level
327 : !
328 0 : iw = 1
329 0 : do kb=1,ngrid
330 0 : k = ngrid-kb+1
331 0 : nx = nxk(k)
332 0 : ny = nyk(k)
333 0 : kpbgn(k) = iw
334 0 : kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)
335 0 : ktxbgn(k) = kcbgn(k)+10*nx*ny
336 0 : ktybgn(k) = ktxbgn(k)+isx*nx*ny
337 0 : iw = ktybgn(k)+jsy*nx*ny
338 0 : ic = kcbgn(k)
339 0 : itx = ktxbgn(k)
340 0 : ity = ktybgn(k)
341 0 : klevel = k
342 0 : call dismh2cr(nx,ny,wk(ic),wk(itx),wk(ity),wk,iwk)
343 : end do
344 : return
345 : end if ! end of intl=0 initialization call block
346 0 : nx = nfx
347 0 : ny = nfy
348 0 : call muh2cr1(nx,ny,rhs,phi,wk,iwk)
349 0 : iparm(17) = itero
350 0 : if (tolmax.gt.0.0_r8) then ! check for convergence
351 0 : fparm(6) = relmax
352 0 : if (relmax.gt.tolmax) ierror = -1 ! flag convergenc failure
353 : end if
354 : return
355 : end subroutine muh2cr
356 : !-----------------------------------------------------------------------
357 0 : subroutine muh2cr1(nx,ny,rhsf,phif,wk,iwk)
358 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
359 : implicit none
360 : integer nx,ny,iwk(*)
361 : real(r8) :: phif(nx,ny),rhsf(nx,ny),wk(*)
362 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
363 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
364 : kcycle,iprer,ipost,intpol,kps
365 : real(r8) :: xa,xb,yc,yd,tolmax,relmax,phmax
366 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
367 : integer k,kb,ip,ic,ir,ipc,irc,icc
368 : integer ncx,ncy,jj,ij,i,j,iter
369 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
370 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
371 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
372 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
373 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),&
374 : nxk(50),nyk(50),isx,jsy
375 : integer ibeta,ialfa,izmat,idmat
376 : common/mh2cr/ibeta,ialfa,izmat,idmat
377 0 : nx = nxk(ngrid)
378 0 : ny = nyk(ngrid)
379 0 : ip = kpbgn(ngrid)
380 0 : ic = kcbgn(ngrid)
381 0 : ir = ic+9*nx*ny
382 : !
383 : ! set phif,rhsf in wk and adjust right hand side
384 : !
385 0 : call swk2(nx,ny,phif,rhsf,wk(ip),wk(ir))
386 0 : if (iguess.eq.0) then
387 : !
388 : ! no initial guess at finest grid level!
389 : !
390 0 : do kb=2,ngrid
391 0 : k = ngrid-kb+1
392 0 : nx = nxk(k+1)
393 0 : ny = nyk(k+1)
394 0 : ip = kpbgn(k+1)
395 0 : ir = kcbgn(k+1)+9*nx*ny
396 0 : ncx = nxk(k)
397 0 : ncy = nyk(k)
398 0 : ipc = kpbgn(k)
399 0 : icc = kcbgn(k)
400 0 : irc = icc+9*ncx*ncy
401 : !
402 : ! transfer down to all grid levels
403 : !
404 0 : call trsfc2(nx,ny,wk(ip),wk(ir),ncx,ncy,wk(ipc),wk(irc))
405 : end do
406 : !
407 : ! adjust right hand side at all grid levels in case
408 : ! rhs or specified b.c. in phi or gbdy changed
409 : !
410 0 : do k=1,ngrid
411 0 : nx = nxk(k)
412 0 : ny = nyk(k)
413 0 : ip = kpbgn(k)
414 0 : ic = kcbgn(k)
415 0 : call adjmh2cr(nx,ny,wk(ip),wk(ic))
416 : end do
417 : !
418 : ! execute one full multigrid cycle
419 : !
420 0 : do k=1,ngrid-1
421 0 : kcur = k
422 0 : call kcymh2cr(wk,iwk)
423 0 : nx = nxk(k+1)
424 0 : ny = nyk(k+1)
425 0 : ip = kpbgn(k+1)
426 0 : ipc = kpbgn(k)
427 0 : ncx = nxk(k)
428 0 : ncy = nyk(k)
429 :
430 : !
431 : ! lift or prolong approximation from k to k+1
432 : !
433 0 : call prolon2(ncx,ncy,wk(ipc),nx,ny,wk(ip),nxa,nxb,nyc,nyd,intpol)
434 : end do
435 : else
436 : !
437 : ! adjust rhs at finest grid level only
438 : !
439 0 : nx = nxk(ngrid)
440 0 : ny = nyk(ngrid)
441 0 : ip = kpbgn(ngrid)
442 0 : ic = kcbgn(ngrid)
443 0 : call adjmh2cr(nx,ny,wk(ip),wk(ic))
444 : end if
445 : !
446 : ! execute maxcy more multigrid k cycles from finest level
447 : !
448 0 : kcur = ngrid
449 0 : do iter=1,maxcy
450 0 : itero = iter
451 0 : call kcymh2cr(wk,iwk)
452 0 : if (tolmax.gt.0.0_r8) then
453 : !
454 : ! error control
455 : !
456 0 : relmax = 0.0_r8
457 0 : phmax = 0.0_r8
458 0 : do j=1,nfy
459 0 : jj = j*(nfx+2)
460 0 : do i=1,nfx
461 0 : ij = jj+i+1
462 0 : phmax = max(phmax,abs(wk(ij)))
463 0 : relmax = max(relmax,abs(wk(ij)-phif(i,j)))
464 :
465 0 : phif(i,j) = wk(ij)
466 : end do
467 : end do
468 : !
469 : ! set maximum relative difference and check for convergence
470 : !
471 0 : if (phmax.gt.0.0_r8) relmax = relmax/phmax
472 :
473 0 : if (relmax.le.tolmax) return
474 : end if
475 : end do
476 : !
477 : ! set final interate after maxcy cycles in phif
478 : !
479 0 : do j=1,nfy
480 0 : jj = j*(nfx+2)
481 0 : do i=1,nfx
482 0 : ij = jj+i+1
483 0 : phif(i,j) = wk(ij)
484 : end do
485 : end do
486 : return
487 : end subroutine muh2cr1
488 : !-----------------------------------------------------------------------
489 0 : subroutine kcymh2cr(wk,iwk)
490 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
491 : !
492 : ! execute multigrid k cycle from kcur grid level
493 : ! kcycle=1 for v cycles, kcycle=2 for w cycles
494 : !
495 : implicit none
496 : integer iwk(*)
497 : real(r8) :: wk(*)
498 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
499 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
500 : kcycle,iprer,ipost,intpol,kps
501 : integer nx,ny,ip,ic,ipc,irc,itx,ity,ncx,ncy,l,nrel
502 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
503 : integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
504 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
505 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
506 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
507 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
508 : common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),&
509 : nxk(50),nyk(50),isx,jsy
510 : integer ibeta,ialfa,izmat,idmat
511 : common/mh2cr/ibeta,ialfa,izmat,idmat
512 : integer kount(50)
513 0 : klevel = kcur
514 0 : nx = nxk(klevel)
515 0 : ny = nyk(klevel)
516 0 : ip = kpbgn(klevel)
517 0 : ic = kcbgn(klevel)
518 0 : itx = ktxbgn(klevel)
519 0 : ity = ktybgn(klevel)
520 0 : if (kcur .eq. 1) then
521 : !
522 : ! solve at coarse level with direct method and return
523 : !
524 0 : if (nyc .ne. 0) then
525 0 : call dir2cr(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),iwk,nxa)
526 0 : return
527 : else
528 0 : call dir2crp(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),&
529 0 : wk(izmat),wk(idmat),iwk,nxa)
530 0 : return
531 : end if
532 : end if
533 : !
534 : ! prerelax at current finest grid level > 1
535 : !
536 0 : do l=1,iprer
537 0 : call relmh2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
538 : end do
539 : !
540 : ! restrict residual to kcur-1 level
541 : !
542 0 : ipc = kpbgn(klevel-1)
543 0 : ncx = nxk(klevel-1)
544 0 : ncy = nyk(klevel-1)
545 0 : irc = kcbgn(klevel-1)+9*ncx*ncy
546 0 : call resmh2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps))
547 : !
548 : ! set counter for grid levels to zero
549 : !
550 0 : do l = 1,kcur
551 0 : kount(l) = 0
552 : end do
553 : !
554 : ! set new grid level and continue k-cycling
555 : !
556 0 : klevel = kcur-1
557 0 : nrel = iprer
558 : !
559 : ! kcycle control point
560 : !
561 : 10 continue
562 : !
563 : ! post relax when kcur revisited
564 : !
565 0 : if (klevel .eq. kcur) go to 5
566 : !
567 : ! count hit at current level
568 : !
569 0 : kount(klevel) = kount(klevel)+1
570 : !
571 : ! relax or solve directly at current level
572 : !
573 0 : nx = nxk(klevel)
574 0 : ny = nyk(klevel)
575 0 : ip = kpbgn(klevel)
576 0 : ic = kcbgn(klevel)
577 0 : itx = ktxbgn(klevel)
578 0 : ity = ktybgn(klevel)
579 0 : if (klevel.gt.1) then
580 0 : do l=1,nrel
581 0 : call relmh2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
582 : end do
583 : else
584 : !
585 : ! use direct method at coarsest level
586 : !
587 0 : if (nyc .ne. 0) then
588 0 : call dir2cr(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),iwk,nxa)
589 : else
590 0 : call dir2crp(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),&
591 0 : wk(izmat),wk(idmat),iwk,nxa)
592 : end if
593 : !
594 : ! insure direct method is not called again at coarse level
595 : !
596 0 : kount(1) = kcycle+1
597 : end if
598 0 : if (kount(klevel) .eq. kcycle+1) then
599 : !
600 : ! kcycle complete at klevel
601 : !
602 0 : ipc = ip
603 0 : ip = kpbgn(klevel+1)
604 0 : ncx = nxk(klevel)
605 0 : ncy = nyk(klevel)
606 0 : nx = nxk(klevel+1)
607 0 : ny = nyk(klevel+1)
608 : !
609 : ! inject correction to finer grid
610 : !
611 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,&
612 0 : intpol,wk(kps))
613 : !
614 : ! reset counter to zero
615 : !
616 0 : kount(klevel) = 0
617 : !
618 : ! ascend to next higher level and set to postrelax there
619 : !
620 0 : klevel = klevel+1
621 0 : nrel = ipost
622 0 : go to 10
623 : else
624 0 : if (klevel .gt. 1) then
625 : !
626 : ! kcycle not complete so descend unless at coarsest grid
627 : !
628 0 : ipc = kpbgn(klevel-1)
629 0 : ncx = nxk(klevel-1)
630 0 : ncy = nyk(klevel-1)
631 0 : irc = kcbgn(klevel-1)+9*ncx*ncy
632 0 : call resmh2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps))
633 : !
634 : ! prerelax at next coarser level
635 : !
636 0 : klevel = klevel-1
637 0 : nrel = iprer
638 0 : go to 10
639 : else
640 : !
641 : ! direct at coarsest level takes place of postrelax
642 : !
643 0 : ip = kpbgn(1)
644 0 : ic = kcbgn(1)
645 0 : nx = nxk(1)
646 0 : ny = nyk(1)
647 0 : if (nyc .ne. 0) then
648 0 : call dir2cr(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),iwk,nxa)
649 : else
650 0 : call dir2crp(nx,ny,wk(ip),wk(ic),wk(ibeta),wk(ialfa),&
651 0 : wk(izmat),wk(idmat),iwk,nxa)
652 : end if
653 0 : ipc = ip
654 0 : ip = kpbgn(2)
655 0 : ncx = nxk(1)
656 0 : ncy = nyk(1)
657 0 : nx = nxk(2)
658 0 : ny = nyk(2)
659 : !
660 : ! inject correction to level 2
661 : !
662 0 : call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,&
663 0 : intpol,wk(kps))
664 : !
665 : ! set to postrelax at level 2
666 : !
667 0 : nrel = ipost
668 0 : klevel = 2
669 0 : go to 10
670 : end if
671 : end if
672 : 5 continue
673 : !
674 : ! post relax at current finest grid level
675 : !
676 0 : nx = nxk(kcur)
677 0 : ny = nyk(kcur)
678 0 : ip = kpbgn(kcur)
679 0 : ic = kcbgn(kcur)
680 0 : itx = ktxbgn(kcur)
681 0 : ity = ktybgn(kcur)
682 0 : do l=1,ipost
683 0 : call relmh2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
684 : end do
685 : return
686 : end subroutine kcymh2cr
687 : !-----------------------------------------------------------------------
688 0 : subroutine dismh2cr(nx,ny,cf,tx,ty,wk,iwk)
689 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
690 : use cam_abortutils ,only: endrun
691 : use edyn_solver_coefs,only: nc,cee,ceee
692 : !
693 : ! discretize elliptic pde for muh2cr, set nonfatal errors
694 : !
695 : implicit none
696 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
697 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
698 : kcycle,iprer,ipost,intpol,kps
699 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
700 : integer nx,ny,iwk(*),i,j,l,im1,jm1
701 : real(r8) :: cf(nx,ny,10),tx(nx,ny,*),ty(ny,nx,*)
702 : real(r8) :: wk(*)
703 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
704 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
705 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
706 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
707 : integer ibeta,ialfa,izmat,idmat
708 : common/mh2cr/ibeta,ialfa,izmat,idmat
709 : integer nnx,nny
710 : !
711 : ! CHECK FOR CONSISTENCYT WRT KLEVEL
712 : !
713 0 : NNX = ixp*2**(KLEVEL-1)+1
714 0 : NNY = jyq*2**(KLEVEL-1)+1
715 0 : IF(NNX.NE.NX.OR.NNY.NE.NY)THEN
716 : ! WRITE(iulog,100)NX,NY,NNX,NNY,ixp,jyq,KLEVEL
717 : ! 100 FORMAT(' INCONSISTENCY WRT LEVEL. NX,NY,NNX,NNY,ixp,jyq,',
718 : ! | 'klevel = ',8I6)
719 0 : call endrun('dismh2cr')
720 : ENDIF
721 0 : call ceee(cee(nc(6-klevel-4)),nx,ny,cf)
722 : !
723 : ! set coefficient for specified boundaries
724 : !
725 0 : if (nxa.eq.1) then
726 0 : i = 1
727 0 : do j=1,ny
728 0 : do l=1,9
729 0 : cf(i,j,l) = 0.0_r8
730 : end do
731 0 : cf(i,j,9) = 1.0_r8
732 : end do
733 : end if
734 0 : if (nxb.eq.1) then
735 0 : i = nx
736 0 : do j=1,ny
737 0 : do l=1,9
738 0 : cf(i,j,l) = 0.0_r8
739 : end do
740 0 : cf(i,j,9) = 1.0_r8
741 : end do
742 : end if
743 0 : if (nyc.eq.1) then
744 0 : j = 1
745 0 : do i=1,nx
746 0 : do l=1,9
747 0 : cf(i,j,l) = 0.0_r8
748 : end do
749 0 : cf(i,j,9) = 1.0_r8
750 : end do
751 : end if
752 0 : if (nyd.eq.1) then
753 0 : j = ny
754 0 : do i=1,nx
755 0 : do l=1,9
756 0 : cf(i,j,l) = 0.0_r8
757 : end do
758 0 : cf(i,j,9) = 1.0_r8
759 : end do
760 : end if
761 0 : if (klevel .eq. 1) then
762 : !
763 : ! set block tri-diagonal coefficient matrix and do lu decomposition
764 : ! for direct method at coarsest grid level
765 : !
766 0 : nx = ixp+1
767 0 : ny = jyq+1
768 0 : if (nyc .ne. 0) then
769 : ! factor non-periodic block matrix
770 0 : call lud2cr(nx,ny,cf,wk(ibeta),wk(ialfa),iwk,nxa)
771 0 : return
772 : else
773 : ! factor periodic block matrix
774 :
775 0 : do j =1,ny-1
776 0 : call setbcr(nx,ny,cf,wk(ibeta),j,nxa)
777 0 : call setacr(nx,ny,cf,wk(ialfa),j,nxa)
778 : end do
779 0 : call lud2crp(nx,ny,cf,wk(ibeta),wk(ialfa),wk(izmat),&
780 0 : wk(idmat),iwk,nxa)
781 0 : return
782 : end if
783 : end if
784 : !
785 : ! set and factor tridiagonal matrices for line relaxation(s) if flagged
786 : !
787 0 : if (method.eq.1.or.method.eq.3) then
788 0 : if (nxa.ne.0) then
789 : !
790 : ! nonperiodic x line relaxation
791 : !
792 0 : do i=1,nx
793 0 : im1 = max0(i-1,1)
794 0 : do j=1,ny
795 0 : tx(im1,j,1) = cf(i,j,5)
796 0 : tx(i,j,2) = cf(i,j,9)
797 0 : tx(i,j,3) = cf(i,j,1)
798 : end do
799 : end do
800 0 : call factri(ny,nx,tx(1,1,1),tx(1,1,2),tx(1,1,3))
801 : else
802 : !
803 : ! periodic x line relaxation
804 : !
805 0 : if (nx .gt. 3) then
806 : !
807 : ! set and factor iff nx > 3
808 : !
809 0 : do i=1,nx-1
810 0 : do j=1,ny
811 0 : tx(i,j,1) = cf(i,j,5)
812 0 : tx(i,j,2) = cf(i,j,9)
813 0 : tx(i,j,3) = cf(i,j,1)
814 : end do
815 : end do
816 0 : call factrp(ny,nx,tx,tx(1,1,2),tx(1,1,3),tx(1,1,4),&
817 0 : tx(1,1,5),wk(kps))
818 : end if
819 : end if
820 : end if
821 :
822 0 : if (method.eq.2.or.method.eq.3) then
823 0 : if (nyc.ne.0) then
824 : !
825 : ! nonperiodic y line relaxation
826 : !
827 0 : do j=1,ny
828 0 : jm1 = max0(j-1,1)
829 0 : do i=1,nx
830 0 : ty(jm1,i,1) = cf(i,j,7)
831 0 : ty(j,i,2) = cf(i,j,9)
832 0 : ty(j,i,3) = cf(i,j,3)
833 : end do
834 : end do
835 0 : call factri(nx,ny,ty(1,1,1),ty(1,1,2),ty(1,1,3))
836 : else
837 : !
838 : ! periodic y line relaxation
839 : !
840 0 : if (ny .gt. 3) then
841 : !
842 : ! set and factor iff ny > 3
843 : !
844 0 : do j=1,ny-1
845 0 : do i=1,nx
846 0 : ty(j,i,1) = cf(i,j,7)
847 0 : ty(j,i,2) = cf(i,j,9)
848 0 : ty(j,i,3) = cf(i,j,3)
849 : end do
850 : end do
851 0 : call factrp(nx,ny,ty,ty(1,1,2),ty(1,1,3),ty(1,1,4),&
852 0 : ty(1,1,5),wk(kps))
853 : end if
854 : end if
855 : end if
856 : return
857 : end subroutine dismh2cr
858 : !-----------------------------------------------------------------------
859 0 : subroutine lud2cr(nx,ny,cof,beta,alfa,index,nxa)
860 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
861 : !
862 : ! decompose nonperiodic block coefficient matrix
863 : !
864 : implicit none
865 : integer nx,ny,nxa,index(nx,ny)
866 : real(r8) :: cof(nx,ny,10),beta(nx,nx,*),alfa(nx,nx,*)
867 : integer iz,i1,jcur,jm1,l,lm1,lp1,k,i
868 : real(r8) :: gama,sum
869 0 : iz = 0
870 0 : i1 = 1
871 : !
872 : ! set and factor umat(1) in beta(1)
873 : !
874 0 : jcur = 1
875 0 : call setbcr(nx,ny,cof,beta,jcur,nxa)
876 0 : call sgfa(beta,nx,nx,index,iz)
877 :
878 0 : do jcur=2,ny
879 : !
880 : ! solve transpose of lmat(jcur)*beta(jcur-1) = alfa(jcur) in alfa(jcur)
881 : !
882 0 : call setacr(nx,ny,cof,alfa,jcur,nxa)
883 0 : call transp(nx,alfa(1,1,jcur))
884 0 : jm1 = jcur-1
885 0 : do l=1,nx
886 0 : call sgsl(beta(1,1,jm1),nx,nx,index(1,jm1),alfa(1,l,jcur),i1)
887 : end do
888 0 : call transp(nx,alfa(1,1,jcur))
889 0 : call setbcr(nx,ny,cof,beta,jcur,nxa)
890 0 : do i=1,nx
891 0 : do l=1,nx
892 0 : sum = 0.0_r8
893 0 : lm1=max0(1,l-1)
894 0 : lp1=min0(l+1,nx)
895 0 : do k=lm1,lp1
896 0 : if (k .eq. l+1) then
897 0 : gama = cof(k,jcur-1,4)
898 0 : else if (k.eq. l) then
899 0 : gama = cof(k,jcur-1,3)
900 0 : else if (k .eq. l-1) then
901 0 : gama = cof(k,jcur-1,2)
902 : else
903 : gama=0.0_r8
904 : end if
905 0 : sum = sum+alfa(i,k,jcur)*gama
906 : end do
907 0 : if (nxa.eq.0) then
908 0 : if (l .eq. 2) then
909 0 : sum=sum+alfa(i,nx,jcur)*cof(nx,jcur-1,2)
910 : end if
911 0 : if (l .eq. nx-1) then
912 0 : sum=sum+alfa(i,1,jcur)*cof(1,jcur-1,4)
913 : end if
914 : end if
915 0 : beta(i,l,jcur) = beta(i,l,jcur)-sum
916 : end do
917 : end do
918 : !
919 : ! factor current beta for next pass
920 : !
921 0 : iz = 0
922 0 : call sgfa(beta(1,1,jcur),nx,nx,index(1,jcur),iz)
923 : end do
924 0 : return
925 : end subroutine lud2cr
926 : !-----------------------------------------------------------------------
927 0 : subroutine dir2cr(nx,ny,phi,cof,beta,alfa,index,nxa)
928 : use shr_kind_mod ,only: r8 => shr_kind_r8
929 : !
930 : ! direct solve at coarsest grid
931 : !
932 : implicit none
933 : integer nx,ny,index(nx,ny),nxa
934 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
935 : real(r8) :: beta(nx,nx,*),alfa(nx,nx,*)
936 : ! forward sweep
937 0 : call for2cr(nx,ny,phi,cof(1,1,10),alfa)
938 : ! backward sweep
939 0 : call bkw2cr(nx,ny,phi,cof,beta,index,nxa)
940 0 : return
941 : end subroutine dir2cr
942 : !-----------------------------------------------------------------------
943 0 : subroutine for2cr(nx,ny,phi,frhs,alfa)
944 : use shr_kind_mod ,only: r8 => shr_kind_r8
945 : !
946 : ! forward sweep
947 : !
948 : implicit none
949 : integer nx,ny,i,j,l
950 : real(r8) :: phi(0:nx+1,0:ny+1),frhs(nx,ny),alfa(nx,nx,*),sum
951 0 : do j=1,ny
952 0 : do i=1,nx
953 0 : phi(i,j)=frhs(i,j)
954 : end do
955 : end do
956 0 : do j=2,ny
957 0 : do i=1,nx
958 0 : sum=0.0_r8
959 0 : do l=1,nx
960 0 : sum=sum+alfa(i,l,j)*phi(l,j-1)
961 : end do
962 0 : phi(i,j)=phi(i,j)-sum
963 : end do
964 : end do
965 0 : return
966 : end subroutine for2cr
967 : !-----------------------------------------------------------------------
968 0 : subroutine bkw2cr(nx,ny,phi,cof,beta,index,nxa)
969 : use shr_kind_mod ,only: r8 => shr_kind_r8
970 : implicit none
971 : integer nx,ny,index(nx,ny),nxa
972 : real(r8) :: beta(nx,nx,*),sum
973 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
974 : integer iz,jcur,jb,j,i
975 0 : iz = 0
976 0 : jcur=ny
977 0 : call sgsl(beta(1,1,jcur),nx ,nx ,index(1,jcur),phi(1,jcur),iz)
978 0 : do jb=2,ny
979 0 : j=ny-jb+1
980 0 : jcur=j
981 0 : do i=2,nx-1
982 0 : sum=cof(i,j,2)*phi(i+1,j+1)+cof(i,j,3)*phi(i,j+1)+cof(i,j,4)* &
983 0 : phi(i-1,j+1)
984 0 : phi(i,j)=phi(i,j)-sum
985 : end do
986 0 : phi(1,j)=phi(1,j)-(cof(1,j,2)*phi(2,j+1)+cof(1,j,3)*phi(1,j+1))
987 0 : phi(nx,j)=phi(nx,j)-(cof(nx,j,3)*phi(nx,j+1)+cof(nx,j,4)* &
988 0 : phi(nx-1,j+1))
989 0 : if (nxa .eq.0) then
990 0 : phi(1,j)=phi(1,j)-cof(1,j,4)*phi(nx-1,j+1)
991 0 : phi(nx,j)=phi(nx,j)-cof(nx,j,2)*phi(2,j+1)
992 : end if
993 0 : call sgsl(beta(1,1,jcur),nx ,nx ,index(1,jcur),phi(1,jcur),iz)
994 : end do
995 0 : return
996 : end subroutine bkw2cr
997 : !-----------------------------------------------------------------------
998 0 : subroutine lud2crp(nx,ny,cof,beta,alfa,zmat,dmat,index,nxa)
999 : use shr_kind_mod ,only: r8 => shr_kind_r8
1000 : !
1001 : ! decompose periodic block tridiagonal matrix for direct at coarsest grid
1002 : !
1003 : implicit none
1004 : integer nx,ny,index(nx,ny),nxa
1005 : real(r8) :: cof(nx,ny,10),alfa(nx,nx,*),beta(nx,nx,*)
1006 : real(r8) :: dmat(nx,nx,*),zmat(nx,nx,*),sum,gama
1007 : integer iz,j,jcur,i,l,jm1,i1,lm1,lp1,k
1008 0 : jcur = 1
1009 : !
1010 : ! set dmat(1)=alfa(1)
1011 : !
1012 0 : call setacr(nx,ny,cof,alfa,jcur,nxa)
1013 0 : do i=1,nx
1014 0 : do l=1,nx
1015 0 : dmat(i,l,1) = alfa(i,l,1)
1016 : end do
1017 : end do
1018 0 : iz = 0
1019 : !
1020 : ! factor umat(1) in beta(1)
1021 : !
1022 0 : call setbcr(nx,ny,cof,beta,jcur,nxa)
1023 0 : call sgfa(beta(1,1,1),nx,nx,index(1,1),iz)
1024 0 : do jcur=2,ny-2
1025 : !
1026 : ! solve transpose of lmat(jcur)umat(jcur-1)=alfa(jcur) in alfa(jcur)
1027 : !
1028 0 : call setacr(nx,ny,cof,alfa,jcur,nxa)
1029 0 : call transp(nx,alfa(1,1,jcur))
1030 0 : jm1 = jcur-1
1031 0 : i1 = 1
1032 0 : do l=1,nx
1033 0 : call sgsl(beta(1,1,jm1),nx,nx,index(1,jm1),alfa(1,l,jcur),i1)
1034 : end do
1035 0 : call transp(nx,alfa(1,1,jcur))
1036 0 : call setbcr(nx,ny,cof,beta,jcur,nxa)
1037 0 : do i=1,nx
1038 0 : do l=1,nx
1039 0 : sum = 0.0_r8
1040 0 : lm1=max0(1,l-1)
1041 0 : lp1=min0(l+1,nx)
1042 0 : do k=lm1,lp1
1043 0 : if (k .eq. l+1) then
1044 0 : gama = cof(k,jcur-1,4)
1045 0 : else if (k.eq. l) then
1046 0 : gama = cof(k,jcur-1,3)
1047 0 : else if (k .eq. l-1) then
1048 0 : gama = cof(k,jcur-1,2)
1049 : else
1050 : gama=0.0_r8
1051 : end if
1052 0 : sum = sum+alfa(i,k,jcur)*gama
1053 : end do
1054 0 : if (nxa.eq.0) then
1055 0 : if (l .eq. 2) then
1056 0 : sum=sum+alfa(i,nx,jcur)*cof(nx,jcur-1,2)
1057 : end if
1058 0 : if (l .eq. nx-1) then
1059 0 : sum=sum+alfa(i,1,jcur)*cof(1,jcur-1,4)
1060 : end if
1061 : end if
1062 0 : beta(i,l,jcur)=beta(i,l,jcur)-sum
1063 : end do
1064 : end do
1065 : !
1066 : ! factor current beta(1,1,jcur) for next pass
1067 : !
1068 0 : call sgfa(beta(1,1,jcur),nx ,nx,index(1,jcur),iz)
1069 : !
1070 : ! set dmat(jcur) = -alfa(jcur)*dmat(jcur-1)
1071 : !
1072 0 : do i=1,nx
1073 0 : do j=1,nx
1074 0 : dmat(i,j,jcur) = 0.0_r8
1075 0 : do l=1,nx
1076 0 : dmat(i,j,jcur) = dmat(i,j,jcur)-alfa(i,l,jcur)* &
1077 0 : dmat(l,j,jcur-1)
1078 : end do
1079 : end do
1080 : end do
1081 0 : if (jcur .eq. ny-2) then
1082 : !
1083 : ! adjust dmat(ny-2) = gama(ny-2)-alfa(ny-2)*dmat(ny-3)
1084 : !
1085 0 : dmat(1,1,jcur) = cof(1,jcur,3) + dmat(1,1,jcur)
1086 0 : dmat(1,2,jcur) = cof(1,jcur,2) + dmat(1,2,jcur)
1087 : !
1088 : ! adjust for periodic b.c. in x
1089 : !
1090 0 : if (nxa .eq. 0) then
1091 0 : dmat(1,nx-1,jcur) = cof(1,jcur,4) + dmat(1,nx-1,jcur)
1092 0 : dmat(nx,2,jcur) = cof(nx,jcur,2) + dmat(nx,2,jcur)
1093 : end if
1094 : !
1095 : ! matrix interior
1096 : !
1097 0 : do i=2,nx-1
1098 0 : dmat(i,i,jcur) = cof(i,jcur,3) + dmat(i,i,jcur)
1099 0 : dmat(i,i-1,jcur) = cof(i,jcur,4) + dmat(i,i-1,jcur)
1100 0 : dmat(i,i+1,jcur) = cof(i,jcur,2) + dmat(i,i+1,jcur)
1101 : end do
1102 0 : dmat(nx,nx,jcur) = cof(nx,jcur,3) + dmat(nx,nx,jcur)
1103 0 : dmat(nx,nx-1,jcur) = cof(nx,jcur,4) + dmat(nx,nx-1,jcur)
1104 : end if
1105 : end do
1106 : !
1107 : ! final phase with periodic factorization
1108 : !
1109 : ! solve transpose of zmat(1) beta(1) = gama(ny-1)
1110 : !
1111 0 : zmat(1,1,1) = cof(1,ny-1,3)
1112 0 : zmat(1,2,1) = cof(1,ny-1,2)
1113 0 : do l=3,nx
1114 0 : zmat(1,l,1) = 0.0_r8
1115 : end do
1116 :
1117 0 : do i=2,nx-1
1118 0 : do l=1,nx
1119 0 : zmat(i,l,1) = 0.0_r8
1120 : end do
1121 0 : zmat(i,i,1) = cof(i,ny-1,3)
1122 0 : zmat(i,i+1,1) = cof(i,ny-1,2)
1123 0 : zmat(i,i-1,1) = cof(i,ny-1,4)
1124 : end do
1125 0 : zmat(nx,nx-1,1) = cof(nx,ny-1,4)
1126 0 : zmat(nx,nx,1) = cof(nx,ny-1,3)
1127 0 : do l=1,nx-2
1128 0 : zmat(nx,l,1) = 0.0_r8
1129 : end do
1130 : !
1131 : ! adjust for periodic x b.c.
1132 : !
1133 0 : if (nxa .eq.0) then
1134 0 : zmat(1,nx-1,1) = cof(1,ny-1,4)
1135 0 : zmat(nx,2,1) = cof(nx,ny-1,2)
1136 : end if
1137 : call transp(nx,zmat(1,1,1))
1138 0 : do l=1,nx
1139 0 : call sgsl(beta(1,1,1),nx,nx,index(1,1),zmat(1,l,1),i1)
1140 : end do
1141 : call transp(nx,zmat(1,1,1))
1142 0 : do jcur = 2,ny-3
1143 : !
1144 : ! solve transpose of zmat(jcur) umat(jcur) = -zmat(jcur-1) gama(jcur-1)
1145 : !
1146 0 : do i=1,nx
1147 0 : zmat(i,1,jcur) = -(zmat(i,1,jcur-1)*cof(1,jcur-1,3) + &
1148 0 : zmat(i,2,jcur-1)*cof(2,jcur-1,4))
1149 : end do
1150 0 : do i=1,nx
1151 0 : do l=2,nx-1
1152 0 : zmat(i,l,jcur) = -(zmat(i,l-1,jcur-1)*cof(l-1,jcur-1,2) + &
1153 : zmat(i,l,jcur-1)*cof(l,jcur-1,3) + &
1154 0 : zmat(i,l+1,jcur-1)*cof(l+1,jcur-1,4))
1155 : end do
1156 : end do
1157 0 : do i=1,nx
1158 0 : zmat(i,nx,jcur) = -(zmat(i,nx-1,jcur-1)*cof(nx-1,jcur-1,2) + &
1159 0 : zmat(i,nx,jcur-1)*cof(nx,jcur-1,3))
1160 : end do
1161 : !
1162 : ! adjust j=2 and j=nx-1 column if periodic in x
1163 : !
1164 0 : if (nxa .eq. 0) then
1165 0 : do i=1,nx
1166 0 : zmat(i,2,jcur)=zmat(i,2,jcur)-zmat(i,nx,jcur-1)* &
1167 0 : cof(nx,jcur-1,2)
1168 0 : zmat(i,nx-1,jcur)=zmat(i,nx-1,jcur)-zmat(i,1,jcur-1)* &
1169 0 : cof(1,jcur-1,4)
1170 : end do
1171 : end if
1172 0 : call transp(nx,zmat(1,1,jcur))
1173 0 : do l=1,nx
1174 0 : call sgsl(beta(1,1,jcur),nx,nx,index(1,jcur),zmat(1,l,jcur),i1)
1175 : end do
1176 0 : call transp(nx,zmat(1,1,jcur))
1177 : end do
1178 : !
1179 : ! solve transpose of zmat(ny-2)umat(ny-2)=alfa(ny-1)-zmat(ny-3)gama(ny-3)
1180 : !
1181 0 : jcur = ny-2
1182 0 : do i=1,nx
1183 0 : zmat(i,1,jcur) = -(zmat(i,1,jcur-1)*cof(1,jcur-1,3) + &
1184 0 : zmat(i,2,jcur-1)*cof(2,jcur-1,4))
1185 : end do
1186 :
1187 0 : do i=1,nx
1188 0 : do l=2,nx-1
1189 0 : zmat(i,l,jcur) = -(zmat(i,l-1,jcur-1)*cof(l-1,jcur-1,2) + &
1190 : zmat(i,l,jcur-1)*cof(l,jcur-1,3) + &
1191 0 : zmat(i,l+1,jcur-1)*cof(l+1,jcur-1,4))
1192 : end do
1193 : end do
1194 0 : do i=1,nx
1195 0 : zmat(i,nx,jcur) = -(zmat(i,nx-1,jcur-1)*cof(nx-1,jcur-1,2) + &
1196 0 : zmat(i,nx,jcur-1)*cof(nx,jcur-1,3))
1197 : end do
1198 : !
1199 : ! adjust j=2 and j=nx-1 column if periodic in x
1200 : !
1201 0 : if (nxa .eq. 0) then
1202 0 : do i=1,nx
1203 0 : zmat(i,2,jcur)=zmat(i,2,jcur)-zmat(i,nx,jcur-1)* &
1204 0 : cof(nx,jcur-1,2)
1205 0 : zmat(i,nx-1,jcur)=zmat(i,nx-1,jcur)-zmat(i,1,jcur-1)* &
1206 0 : cof(1,jcur-1,4)
1207 : end do
1208 : end if
1209 0 : call setacr(nx,ny,cof,alfa,ny-1,nxa)
1210 0 : do i=1,nx
1211 0 : do l=1,nx
1212 0 : zmat(i,l,ny-2) = alfa(i,l,ny-1) + zmat(i,l,ny-2)
1213 : end do
1214 : end do
1215 0 : call transp(nx,zmat(1,1,ny-2))
1216 0 : do l=1,nx
1217 0 : call sgsl(beta(1,1,ny-2),nx,nx,index(1,ny-2),zmat(1,l,ny-2),i1)
1218 : end do
1219 0 : call transp(nx,zmat(1,1,ny-2))
1220 : !
1221 : ! set umat(ny-1) = beta(ny-1)-(zmat(1)*dmat(1)+...+zmat(ny-2)*dmat(ny-2))
1222 : ! in beta(ny-1)
1223 : !
1224 0 : call setbcr(nx,ny,cof,beta,ny-1,nxa)
1225 0 : do i=1,nx
1226 0 : do j=1,nx
1227 0 : sum = 0.0_r8
1228 0 : do jcur=1,ny-2
1229 0 : do l=1,nx
1230 0 : sum = sum + zmat(i,l,jcur)*dmat(l,j,jcur)
1231 : end do
1232 : end do
1233 0 : beta(i,j,ny-1) = beta(i,j,ny-1) - sum
1234 : end do
1235 : end do
1236 : !
1237 : ! factor bmat(ny-1) for backward sweep
1238 : !
1239 0 : call sgfa(beta(1,1,ny-1),nx,nx,index(1,ny-1),iz)
1240 : !
1241 : ! lud is now complete
1242 : !
1243 0 : return
1244 : end subroutine lud2crp
1245 : !-----------------------------------------------------------------------
1246 0 : subroutine dir2crp(nx,ny,phi,cof,beta,alfa,zmat,dmat,index,nxa)
1247 : use shr_kind_mod ,only: r8 => shr_kind_r8
1248 : implicit none
1249 : integer nx,ny,index(nx,ny),nxa
1250 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
1251 : real(r8) :: beta(nx,nx,*),alfa(nx,nx,*)
1252 : real(r8) :: zmat(nx,nx,*), dmat(nx,nx,*)
1253 : ! forward sweep
1254 0 : call for2crp(nx,ny,phi,cof(1,1,10),alfa,zmat)
1255 : ! backward sweep
1256 0 : call bkw2crp(nx,ny,phi,cof,beta,dmat,index,nxa)
1257 0 : return
1258 : end subroutine dir2crp
1259 : !-----------------------------------------------------------------------
1260 0 : subroutine for2crp(nx,ny,phi,frhs,alfa,zmat)
1261 : use shr_kind_mod ,only: r8 => shr_kind_r8
1262 : implicit none
1263 : integer nx,ny,i,j,l,jcur,k
1264 : real(r8) :: frhs(nx,ny)
1265 : real(r8) :: phi(0:nx+1,0:ny+1)
1266 :
1267 : real(r8) :: alfa(nx,nx,*),zmat(nx,nx,*)
1268 : real(r8) :: sum
1269 0 : do j=1,ny-1
1270 0 : do i=1,nx
1271 0 : phi(i,j)=frhs(i,j)
1272 : end do
1273 : end do
1274 0 : do jcur=2,ny-2
1275 0 : do i=1,nx
1276 0 : sum=0.0_r8
1277 0 : do l=1,nx
1278 0 : sum=sum+alfa(i,l,jcur)*phi(l,jcur-1)
1279 : end do
1280 0 : phi(i,jcur)=phi(i,jcur)-sum
1281 : end do
1282 : end do
1283 : !
1284 : ! solve:
1285 : ! zmat(1)*phi(1)+...+zmat(ny-2)*phi(ny-2) + phi(ny-1) = f(ny-1)
1286 : !
1287 0 : do i=1,nx
1288 0 : sum = 0.0_r8
1289 0 : do k=1,ny-2
1290 0 : do l=1,nx
1291 0 : sum = sum + zmat(i,l,k)*phi(l,k)
1292 : end do
1293 : end do
1294 0 : phi(i,ny-1) = phi(i,ny-1) - sum
1295 : end do
1296 0 : return
1297 : end subroutine for2crp
1298 : !-----------------------------------------------------------------------
1299 0 : subroutine bkw2crp(nx,ny,phi,cof,beta,dmat,index,nxa)
1300 : use shr_kind_mod ,only: r8 => shr_kind_r8
1301 : implicit none
1302 : integer nx,ny,index(nx,ny),nxa
1303 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
1304 : real(r8) :: beta(nx,nx,ny),dmat(nx,nx,*)
1305 : integer iz,i,l,kb,k
1306 : real(r8) :: sum
1307 0 : iz = 0
1308 0 : call sgsl(beta(1,1,ny-1),nx,nx,index(1,ny-1),phi(1,ny-1),iz)
1309 : !
1310 : ! solve beta(ny-2)*phi(ny-2) = phi(ny-2)-dmat(ny-2)*phi(ny-1)
1311 : !
1312 0 : do i=1,nx
1313 0 : sum = 0.0_r8
1314 0 : do l=1,nx
1315 0 : sum = sum + dmat(i,l,ny-2)*phi(l,ny-1)
1316 : end do
1317 0 : phi(i,ny-2) = phi(i,ny-2) - sum
1318 : end do
1319 0 : call sgsl(beta(1,1,ny-2),nx,nx,index(1,ny-2),phi(1,ny-2),iz)
1320 : !
1321 : ! solve beta(k)*phi(k) = phi(k) - gama(k)*phi(k+1)-dmat(k)*phi(ny-1)
1322 : ! k=ny-3,...,1
1323 : !
1324 0 : do kb=4,ny
1325 0 : k = ny-kb+1
1326 0 : sum = 0.0_r8
1327 0 : do l=1,nx
1328 0 : sum = sum+dmat(1,l,k)*phi(l,ny-1)
1329 : end do
1330 0 : phi(1,k) = phi(1,k)-sum - ( cof(1,k,3)*phi(1,k+1) + &
1331 0 : cof(1,k,2)*phi(2,k+1))
1332 0 : do i=2,nx-1
1333 0 : sum = 0.0_r8
1334 0 : do l=1,nx
1335 0 : sum = sum+dmat(i,l,k)*phi(l,ny-1)
1336 : end do
1337 0 : phi(i,k) = phi(i,k) - sum - (cof(i,k,4)*phi(i-1,k+1) + &
1338 : cof(i,k,3)*phi(i,k+1) + &
1339 0 : cof(i,k,2)*phi(i+1,k+1))
1340 : end do
1341 0 : sum = 0.0_r8
1342 0 : do l=1,nx
1343 0 : sum = sum+dmat(nx,l,k)*phi(l,ny-1)
1344 : end do
1345 0 : phi(nx,k) = phi(nx,k) - sum - (cof(nx,k,4)*phi(nx-1,k+1) + &
1346 0 : cof(nx,k,3)*phi(nx,k+1))
1347 : !
1348 : ! adjust for periodic x b.c.
1349 : !
1350 0 : if (nxa .eq. 0) then
1351 0 : phi(1,k) = phi(1,k) - cof(1,k,4)*phi(nx-1,k+1)
1352 0 : phi(nx,k) = phi(nx,k) - cof(nx,k,2)*phi(2,k+1)
1353 : end if
1354 0 : call sgsl(beta(1,1,k),nx,nx,index(1,k),phi(1,k),iz)
1355 : end do
1356 : !
1357 : ! set j=ny by periodicity
1358 : !
1359 0 : do i=1,nx
1360 0 : phi(i,ny) = phi(i,1)
1361 : end do
1362 0 : return
1363 : end subroutine bkw2crp
1364 : !-----------------------------------------------------------------------
1365 0 : subroutine setbcr(nx,ny,cof,beta,jcur,nxa)
1366 : use shr_kind_mod ,only: r8 => shr_kind_r8
1367 : !
1368 : ! set diagonal matrix on block
1369 : !
1370 : implicit none
1371 : integer nx,ny,jcur,nxa,i,l
1372 : real(r8) :: cof(nx,ny,10),beta(nx,nx,*)
1373 0 : do i=1,nx
1374 0 : do l=1,nx
1375 0 : beta(i,l,jcur)=0.0_r8
1376 : end do
1377 : end do
1378 0 : do i=1,nx
1379 0 : beta(i,i,jcur) = cof(i,jcur,9)
1380 : end do
1381 0 : do i=2,nx
1382 0 : beta(i,i-1,jcur) = cof(i,jcur,5)
1383 : end do
1384 0 : do i=1,nx-1
1385 0 : beta(i,i+1,jcur) = cof(i,jcur,1)
1386 : end do
1387 0 : if (nxa.eq.0) then
1388 0 : beta(1,nx-1,jcur) = cof(1,jcur,5)
1389 0 : beta(nx,2,jcur) = cof(nx,jcur,1)
1390 : end if
1391 0 : return
1392 : end subroutine setbcr
1393 : !-----------------------------------------------------------------------
1394 0 : subroutine setacr(nx,ny,cof,alfa,jcur,nxa)
1395 : use shr_kind_mod ,only: r8 => shr_kind_r8
1396 : implicit none
1397 : integer nx,ny,jcur,nxa,i,j
1398 : real(r8) :: cof(nx,ny,10),alfa(nx,nx,*)
1399 0 : do i=1,nx
1400 0 : do j=1,nx
1401 0 : alfa(i,j,jcur)=0.0_r8
1402 : end do
1403 : end do
1404 0 : do i=2,nx
1405 0 : alfa(i,i-1,jcur)=cof(i,jcur,6)
1406 : end do
1407 0 : do i=1,nx
1408 0 : alfa(i,i,jcur)=cof(i,jcur,7)
1409 : end do
1410 0 : do i=1,nx-1
1411 0 : alfa(i,i+1,jcur)=cof(i,jcur,8)
1412 : end do
1413 0 : if (nxa .eq. 0) then
1414 : ! adjust for x periodicity
1415 0 : alfa(1,nx-1,jcur)=cof(1,jcur,6)
1416 0 : alfa(nx,2,jcur)=cof(nx,jcur,8)
1417 : end if
1418 0 : return
1419 : end subroutine setacr
1420 : !-----------------------------------------------------------------------
1421 0 : subroutine adjmh2cr(nx,ny,phi,cf)
1422 : use shr_kind_mod ,only: r8 => shr_kind_r8
1423 : !
1424 : ! adjust righthand side in cf(i,j,10) for boundary conditions
1425 : !
1426 : implicit none
1427 : integer :: intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
1428 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
1429 : kcycle,iprer,ipost,intpol,kps
1430 : real(r8) :: xa,xb,yc,yd,tolmax,relmax
1431 : integer :: nx,ny,i,j
1432 : real(r8) :: cf(nx,ny,10),phi(0:nx+1,0:ny+1)
1433 :
1434 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
1435 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
1436 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1437 : common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
1438 :
1439 :
1440 : !
1441 : ! set specified boundaries in rhs from phi
1442 : !
1443 0 : if (nxa.eq.1) then
1444 0 : i = 1
1445 0 : do j=1,ny
1446 0 : cf(i,j,10) = phi(i,j)
1447 : end do
1448 : end if
1449 0 : if (nxb.eq.1) then
1450 0 : i = nx
1451 0 : do j=1,ny
1452 0 : cf(i,j,10) = phi(i,j)
1453 : end do
1454 : end if
1455 0 : if (nyc.eq.1) then
1456 0 : j = 1
1457 0 : do i=1,nx
1458 0 : cf(i,j,10) = phi(i,j)
1459 : end do
1460 : end if
1461 0 : if (nyd.eq.1) then
1462 0 : j = ny
1463 0 : do i=1,nx
1464 0 : cf(i,j,10) = phi(i,j)
1465 : end do
1466 : end if
1467 0 : return
1468 : end subroutine adjmh2cr
1469 : !-----------------------------------------------------------------------
1470 0 : subroutine resmh2cr(nx,ny,phi,ncx,ncy,phic,rhsc,cof,resf)
1471 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
1472 : !
1473 : ! restrict residual from fine to coarse mesh using fully weighted
1474 : ! residual restriction
1475 : !
1476 : implicit none
1477 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
1478 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
1479 : kcycle,iprer,ipost,intpol,kps
1480 : integer nx,ny,ncx,ncy,i,j,ic,jc
1481 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
1482 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
1483 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1484 : real(r8) :: rhsc(ncx,ncy),resf(nx,ny)
1485 : real(r8) :: phi(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
1486 : real(r8) :: cof(nx,ny,10)
1487 : !
1488 : ! set phic zero
1489 : !
1490 0 : do jc=0,ncy+1
1491 0 : do ic=0,ncx+1
1492 0 : phic(ic,jc) = 0.0_r8
1493 : end do
1494 : end do
1495 : !
1496 : ! compute residual on fine mesh in resf
1497 : !
1498 : !!$OMP PARALLEL DO SHARED(resf,cof,phi,nx,ny) PRIVATE(i,j)
1499 0 : do j=1,ny
1500 0 : do i=1,nx
1501 0 : resf(i,j) = cof(i,j,10)-( &
1502 0 : cof(i,j,1)*phi(i+1,j)+ &
1503 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1504 0 : cof(i,j,3)*phi(i,j+1)+ &
1505 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1506 : cof(i,j,5)*phi(i-1,j)+ &
1507 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1508 : cof(i,j,7)*phi(i,j-1)+ &
1509 : cof(i,j,8)*phi(i+1,j-1)+ &
1510 0 : cof(i,j,9)*phi(i,j))
1511 : end do
1512 : end do
1513 : !
1514 : ! restrict resf to coarse mesh in rhsc
1515 : !
1516 0 : call res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
1517 0 : return
1518 : end subroutine resmh2cr
1519 : !-----------------------------------------------------------------------
1520 0 : subroutine relmh2cr(nx,ny,phi,cof,tx,ty,sum)
1521 : use shr_kind_mod ,only: r8 => shr_kind_r8
1522 : !
1523 : ! relaxation for muh2cr
1524 : !
1525 : implicit none
1526 : integer nx,ny
1527 : real(r8) :: phi(*),cof(*),tx(*),ty(*),sum(*)
1528 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
1529 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
1530 : kcycle,iprer,ipost,intpol,kps
1531 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
1532 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
1533 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1534 0 : if (method.eq.0) then ! point relaxation
1535 0 : call relmh2crp(nx,ny,phi,cof)
1536 0 : else if (method.eq.1) then ! line x relaxation
1537 0 : call slxmh2cr(nx,ny,phi,cof,tx,sum)
1538 0 : else if (method.eq.2) then ! line y relaxation
1539 0 : call slymh2cr(nx,ny,phi,cof,ty,sum)
1540 0 : else if (method.eq.3) then ! line x&y relaxation
1541 0 : call slxmh2cr(nx,ny,phi,cof,tx,sum)
1542 0 : call slymh2cr(nx,ny,phi,cof,ty,sum)
1543 : end if
1544 0 : return
1545 : end subroutine relmh2cr
1546 : !-----------------------------------------------------------------------
1547 0 : subroutine relmh2crp(nx,ny,phi,cof)
1548 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
1549 : !
1550 : ! gauss-seidel four color point relaxation
1551 : !
1552 : implicit none
1553 : integer nx,ny,i,j,lcolor,i1,i2,i3,i4,it
1554 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
1555 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
1556 : kcycle,iprer,ipost,intpol,kps
1557 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
1558 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
1559 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1560 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10)
1561 0 : i1 = 1
1562 0 : i2 = 4
1563 0 : i3 = 3
1564 0 : i4 = 2
1565 : !
1566 : ! sweep four colored grid points
1567 : !
1568 0 : do lcolor=1,4
1569 : !!$OMP PARALLEL DO SHARED(i1,cof,phi,nx,ny) PRIVATE(i,j)
1570 0 : do j=1,ny,4
1571 0 : do i=i1,nx,4
1572 0 : phi(i,j) = (cof(i,j,10) - ( &
1573 0 : cof(i,j,1)*phi(i+1,j) + &
1574 0 : cof(i,j,2)*phi(i+1,j+1) + &
1575 : cof(i,j,3)*phi(i,j+1) + &
1576 0 : cof(i,j,4)*phi(i-1,j+1) + &
1577 : cof(i,j,5)*phi(i-1,j) + &
1578 0 : cof(i,j,6)*phi(i-1,j-1) + &
1579 : cof(i,j,7)*phi(i,j-1) + &
1580 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
1581 : end do
1582 : end do
1583 : !!$OMP PARALLEL DO SHARED(i2,cof,phi,nx,ny) PRIVATE(i,j)
1584 0 : do j=2,ny,4
1585 0 : do i=i2,nx,4
1586 0 : phi(i,j) = (cof(i,j,10) - ( &
1587 0 : cof(i,j,1)*phi(i+1,j) + &
1588 0 : cof(i,j,2)*phi(i+1,j+1) + &
1589 : cof(i,j,3)*phi(i,j+1) + &
1590 0 : cof(i,j,4)*phi(i-1,j+1) + &
1591 : cof(i,j,5)*phi(i-1,j) + &
1592 0 : cof(i,j,6)*phi(i-1,j-1) + &
1593 : cof(i,j,7)*phi(i,j-1) + &
1594 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
1595 : end do
1596 : end do
1597 : !!$OMP PARALLEL DO SHARED(i3,cof,phi,nx,ny) PRIVATE(i,j)
1598 0 : do j=3,ny,4
1599 0 : do i=i3,nx,4
1600 0 : phi(i,j) = (cof(i,j,10) - ( &
1601 0 : cof(i,j,1)*phi(i+1,j) + &
1602 0 : cof(i,j,2)*phi(i+1,j+1) + &
1603 : cof(i,j,3)*phi(i,j+1) + &
1604 0 : cof(i,j,4)*phi(i-1,j+1) + &
1605 : cof(i,j,5)*phi(i-1,j) + &
1606 0 : cof(i,j,6)*phi(i-1,j-1) + &
1607 : cof(i,j,7)*phi(i,j-1) + &
1608 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
1609 : end do
1610 : end do
1611 : !!$OMP PARALLEL DO SHARED(i4,cof,phi,nx,ny) PRIVATE(i,j)
1612 0 : do j=4,ny,4
1613 0 : do i=i4,nx,4
1614 0 : phi(i,j) = (cof(i,j,10) - ( &
1615 0 : cof(i,j,1)*phi(i+1,j) + &
1616 0 : cof(i,j,2)*phi(i+1,j+1) + &
1617 : cof(i,j,3)*phi(i,j+1) + &
1618 0 : cof(i,j,4)*phi(i-1,j+1) + &
1619 : cof(i,j,5)*phi(i-1,j) + &
1620 0 : cof(i,j,6)*phi(i-1,j-1) + &
1621 : cof(i,j,7)*phi(i,j-1) + &
1622 0 : cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
1623 : end do
1624 : end do
1625 : !
1626 : ! set periodic virtual boundaries as necessary
1627 : !
1628 0 : if (nxa.eq.0) then
1629 0 : do j=1,ny
1630 0 : phi(0,j) = phi(nx-1,j)
1631 0 : phi(nx+1,j) = phi(2,j)
1632 : end do
1633 : end if
1634 0 : if (nyc.eq.0) then
1635 0 : do i=1,nx
1636 0 : phi(i,0) = phi(i,ny-1)
1637 0 : phi(i,ny+1) = phi(i,2)
1638 : end do
1639 : end if
1640 : !
1641 : ! permute (i1,i2,i3,i4) for next color
1642 : !
1643 0 : it = i4
1644 0 : i4 = i3
1645 0 : i3 = i2
1646 0 : i2 = i1
1647 0 : i1 = it
1648 : end do
1649 0 : return
1650 : end subroutine relmh2crp
1651 : !-----------------------------------------------------------------------
1652 0 : subroutine slxmh2cr(nx,ny,phi,cof,tx,sum)
1653 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
1654 : !
1655 : ! line relaxation in the x direction (periodic or nonperiodic)
1656 : !
1657 : implicit none
1658 : integer nx,ny,i,ib,j
1659 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess, &
1660 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur, &
1661 : kcycle,iprer,ipost,intpol,kps
1662 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy, &
1663 : iguess, maxcy,method,nwork,lwork,itero,ngrid, &
1664 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1665 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10),tx(nx,ny,*),sum(ny)
1666 : !
1667 : ! set periodic y virtual boundary if necessary
1668 : !
1669 0 : if (nyc.eq.0) then
1670 0 : do i=1,nx
1671 0 : phi(i,0) = phi(i,ny-1)
1672 0 : phi(i,ny+1) = phi(i,2)
1673 : end do
1674 : end if
1675 :
1676 0 : if (nxa.ne.0) then
1677 : !!$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,ib,j)
1678 : !
1679 : ! x direction not periodic, sweep odd j lines
1680 : !
1681 0 : do j=1,ny,2
1682 0 : do i=1,nx
1683 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
1684 : cof(i,j,3)*phi(i,j+1)+ &
1685 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1686 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1687 : cof(i,j,7)*phi(i,j-1)+ &
1688 0 : cof(i,j,8)*phi(i+1,j-1))
1689 : end do
1690 : !
1691 : ! forward sweep
1692 : !
1693 0 : do i=2,nx
1694 0 : phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
1695 : end do
1696 : !
1697 : ! backward sweep
1698 : !
1699 0 : phi(nx,j) = phi(nx,j)/tx(nx,j,2)
1700 0 : do ib=2,nx
1701 0 : i = nx-ib+1
1702 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
1703 : end do
1704 : end do
1705 : !
1706 : ! sweep even j lines forward and back
1707 : !
1708 : !!$OMP PARALLEL DO SHARED(cof,phi,tx,nx,ny) PRIVATE(i,ib,j)
1709 0 : do j=2,ny,2
1710 0 : do i=1,nx
1711 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
1712 : cof(i,j,3)*phi(i,j+1)+ &
1713 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1714 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1715 : cof(i,j,7)*phi(i,j-1)+ &
1716 0 : cof(i,j,8)*phi(i+1,j-1))
1717 : end do
1718 0 : do i=2,nx
1719 0 : phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
1720 : end do
1721 0 : phi(nx,j) = phi(nx,j)/tx(nx,j,2)
1722 0 : do ib=2,nx
1723 0 : i = nx-ib+1
1724 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
1725 : end do
1726 : end do
1727 : else
1728 : !
1729 : ! x direction periodic
1730 : !
1731 0 : do j=1,ny
1732 0 : sum(j) = 0.0_r8
1733 0 : phi(0,j) = phi(nx-1,j)
1734 0 : phi(nx+1,j) = phi(2,j)
1735 : end do
1736 : !
1737 : ! sweep odd lines forward and back
1738 : !
1739 : !!$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib)
1740 0 : do j=1,ny,2
1741 0 : do i=1,nx-1
1742 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
1743 : cof(i,j,3)*phi(i,j+1)+ &
1744 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1745 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1746 : cof(i,j,7)*phi(i,j-1)+ &
1747 0 : cof(i,j,8)*phi(i+1,j-1))
1748 : end do
1749 : !
1750 : ! forward sweep
1751 : !
1752 0 : do i=2,nx-2
1753 0 : phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
1754 : end do
1755 0 : do i=1,nx-2
1756 0 : sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
1757 : end do
1758 0 : phi(nx-1,j) = phi(nx-1,j)-sum(j)
1759 : !
1760 : ! backward sweep
1761 : !
1762 0 : phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
1763 0 : phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ &
1764 0 : tx(nx-2,j,2)
1765 0 : do ib=4,nx
1766 0 : i = nx-ib+1
1767 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* &
1768 0 : phi(nx-1,j))/tx(i,j,2)
1769 : end do
1770 : end do
1771 : !
1772 : ! set periodic and virtual points for j odd
1773 : !
1774 0 : do j=1,ny,2
1775 0 : phi(nx,j) = phi(1,j)
1776 0 : phi(0,j) = phi(nx-1,j)
1777 0 : phi(nx+1,j) = phi(2,j)
1778 : end do
1779 : !
1780 : ! sweep even j lines
1781 : !
1782 : !!$OMP PARALLEL DO SHARED(sum,cof,phi,tx,nx,ny) PRIVATE(i,j,ib)
1783 0 : do j=2,ny,2
1784 0 : do i=1,nx-1
1785 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+ &
1786 : cof(i,j,3)*phi(i,j+1)+ &
1787 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1788 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1789 : cof(i,j,7)*phi(i,j-1)+ &
1790 0 : cof(i,j,8)*phi(i+1,j-1))
1791 : end do
1792 : !
1793 : ! forward sweep
1794 : !
1795 0 : do i=2,nx-2
1796 0 : phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
1797 : end do
1798 0 : do i=1,nx-2
1799 0 : sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
1800 : end do
1801 0 : phi(nx-1,j) = phi(nx-1,j)-sum(j)
1802 : !
1803 : ! backward sweep
1804 : !
1805 0 : phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
1806 0 : phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/ &
1807 0 : tx(nx-2,j,2)
1808 0 : do ib=4,nx
1809 0 : i = nx-ib+1
1810 0 : phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)* &
1811 0 : phi(nx-1,j))/tx(i,j,2)
1812 : end do
1813 : end do
1814 : !
1815 : ! set periodic and virtual points for j even
1816 : !
1817 0 : do j=2,ny,2
1818 0 : phi(nx,j) = phi(1,j)
1819 0 : phi(0,j) = phi(nx-1,j)
1820 0 : phi(nx+1,j) = phi(2,j)
1821 : end do
1822 : end if
1823 : !
1824 : ! set periodic y virtual boundaries if necessary
1825 : !
1826 0 : if (nyc.eq.0) then
1827 0 : do i=1,nx
1828 0 : phi(i,0) = phi(i,ny-1)
1829 0 : phi(i,ny+1) = phi(i,2)
1830 : end do
1831 : end if
1832 0 : return
1833 : end subroutine slxmh2cr
1834 : !-----------------------------------------------------------------------
1835 0 : subroutine slymh2cr(nx,ny,phi,cof,ty,sum)
1836 0 : use shr_kind_mod ,only: r8 => shr_kind_r8
1837 : !
1838 : ! line relaxation in the y direction (periodic or nonperiodic)
1839 : !
1840 : implicit none
1841 : integer nx,ny,i,j,jb
1842 : integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,&
1843 : maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,&
1844 : kcycle,iprer,ipost,intpol,kps
1845 : common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,&
1846 : iguess, maxcy,method,nwork,lwork,itero,ngrid,&
1847 : klevel,kcur,kcycle,iprer,ipost,intpol,kps
1848 : real(r8) :: phi(0:nx+1,0:ny+1),cof(nx,ny,10),ty(ny,nx,*),sum(nx)
1849 : !
1850 : ! set periodic and virtual x boundaries if necessary
1851 : !
1852 0 : if (nxa.eq.0) then
1853 0 : do j=1,ny
1854 0 : phi(0,j) = phi(nx-1,j)
1855 0 : phi(nx,j) = phi(1,j)
1856 0 : phi(nx+1,j) = phi(2,j)
1857 : end do
1858 : end if
1859 :
1860 0 : if (nyc.ne.0) then
1861 : !
1862 : ! y direction not periodic
1863 : !
1864 : !!$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1865 0 : do i=1,nx,2
1866 0 : do j=1,ny
1867 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1868 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1869 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1870 : cof(i,j,5)*phi(i-1,j)+ &
1871 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1872 0 : cof(i,j,8)*phi(i+1,j-1))
1873 : end do
1874 : !
1875 : ! forward sweep thru odd x lines
1876 : !
1877 0 : do j=2,ny
1878 0 : phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
1879 : end do
1880 : !
1881 : ! backward sweep
1882 : !
1883 0 : phi(i,ny) = phi(i,ny)/ty(ny,i,2)
1884 0 : do jb=2,ny
1885 0 : j = ny-jb+1
1886 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
1887 : end do
1888 : end do
1889 : !
1890 : ! forward sweep even x lines
1891 : !
1892 : !!$OMP PARALLEL DO SHARED(cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1893 0 : do i=2,nx,2
1894 0 : do j=1,ny
1895 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1896 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1897 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1898 : cof(i,j,5)*phi(i-1,j)+ &
1899 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1900 0 : cof(i,j,8)*phi(i+1,j-1))
1901 : end do
1902 0 : do j=2,ny
1903 0 : phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
1904 : end do
1905 : !
1906 : ! backward sweep
1907 : !
1908 0 : phi(i,ny) = phi(i,ny)/ty(ny,i,2)
1909 0 : do jb=2,ny
1910 0 : j = ny-jb+1
1911 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
1912 : end do
1913 : end do
1914 : else
1915 : !
1916 : ! y direction periodic
1917 : !
1918 0 : do i=1,nx
1919 0 : sum(i) = 0.0_r8
1920 0 : phi(i,0) = phi(i,ny-1)
1921 0 : phi(i,ny) = phi(i,1)
1922 0 : phi(i,ny+1) = phi(i,2)
1923 : end do
1924 : !
1925 : ! forward sweep odd x lines
1926 : !
1927 : !!$OMP PARALLEL DO SHARED(sum,cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1928 0 : do i=1,nx,2
1929 0 : do j=1,ny-1
1930 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1931 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1932 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1933 : cof(i,j,5)*phi(i-1,j)+ &
1934 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1935 0 : cof(i,j,8)*phi(i+1,j-1))
1936 : end do
1937 0 : do j=2,ny-2
1938 0 : phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
1939 : end do
1940 0 : do j=1,ny-2
1941 0 : sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
1942 : end do
1943 0 : phi(i,ny-1) = phi(i,ny-1)-sum(i)
1944 : !
1945 : ! backward sweep
1946 : !
1947 0 : phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
1948 0 : phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ &
1949 0 : ty(ny-2,i,2)
1950 0 : do jb=4,ny
1951 0 : j = ny-jb+1
1952 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* &
1953 0 : phi(i,ny-1))/ty(j,i,2)
1954 : end do
1955 : end do
1956 : !
1957 : ! set odd periodic and virtual y boundaries
1958 : !
1959 0 : do i=1,nx,2
1960 0 : phi(i,0) = phi(i,ny-1)
1961 0 : phi(i,ny) = phi(i,1)
1962 0 : phi(i,ny+1) = phi(i,2)
1963 : end do
1964 : !
1965 : ! forward sweep even x lines
1966 : !
1967 : !!$OMP PARALLEL DO SHARED(sum,cof,phi,ty,nx,ny) PRIVATE(i,j,jb)
1968 0 : do i=2,nx,2
1969 0 : do j=1,ny-1
1970 0 : phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+ &
1971 0 : cof(i,j,2)*phi(i+1,j+1)+ &
1972 0 : cof(i,j,4)*phi(i-1,j+1)+ &
1973 : cof(i,j,5)*phi(i-1,j)+ &
1974 0 : cof(i,j,6)*phi(i-1,j-1)+ &
1975 0 : cof(i,j,8)*phi(i+1,j-1))
1976 : end do
1977 0 : do j=2,ny-2
1978 0 : phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
1979 : end do
1980 0 : do j=1,ny-2
1981 0 : sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
1982 : end do
1983 0 : phi(i,ny-1) = phi(i,ny-1)-sum(i)
1984 : !
1985 : ! backward sweep
1986 : !
1987 0 : phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
1988 0 : phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/ &
1989 0 : ty(ny-2,i,2)
1990 0 : do jb=4,ny
1991 0 : j = ny-jb+1
1992 0 : phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)* &
1993 0 : phi(i,ny-1))/ty(j,i,2)
1994 : end do
1995 : end do
1996 : !
1997 : ! set even periodic and virtual y boundaries
1998 : !
1999 0 : do i=2,nx,2
2000 0 : phi(i,0) = phi(i,ny-1)
2001 0 : phi(i,ny) = phi(i,1)
2002 0 : phi(i,ny+1) = phi(i,2)
2003 : end do
2004 : end if
2005 : !
2006 : ! set periodic and virtual x boundaries if necessary
2007 : !
2008 0 : if (nxa.eq.0) then
2009 0 : do j=1,ny
2010 0 : phi(0,j) = phi(nx-1,j)
2011 0 : phi(nx+1,j) = phi(2,j)
2012 : end do
2013 : end if
2014 0 : return
2015 : end subroutine slymh2cr
2016 : !-----------------------------------------------------------------------
2017 : end module edyn_muh2cr
|
