Line data Source code
1 : module edyn_mudcom
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 :
4 : implicit none
5 :
6 : private
7 :
8 : public :: cor2
9 : public :: factri
10 : public :: factrp
11 : public :: swk2
12 : public :: trsfc2
13 : public :: prolon2
14 : public :: res2
15 : public :: sgfa
16 : public :: sgsl
17 : public :: transp
18 :
19 : !-----------------------------------------------------------------------
20 : contains
21 : !-----------------------------------------------------------------------
22 : !
23 : ! file mudcom.f
24 : ! . .
25 : ! . MUDPACK version 4.0 .
26 : !
27 : ! ... author and specialist
28 : !
29 : ! John C. Adams (National Center for Atmospheric Research) (retired)
30 :
31 : ! ... For MUDPACK information, visit the website:
32 : ! (https://www2.cisl.ucar.edu/resources/legacy/mudpack)
33 : !
34 : ! ... purpose
35 : !
36 : ! mudcom.f is a common subroutines file containing subroutines
37 : ! called by some or all of the real two- and three-dimensional
38 : ! mudpack solvers. mudcom.f must be loaded with any real mudpack
39 : ! solver.
40 : !
41 : ! cb mud2cr1: call swk2(nx,ny,phif,rhsf,wk(ip),wk(ir))
42 : !
43 : !-----------------------------------------------------------------------
44 19 : subroutine swk2(nfx,nfy,phif,rhsf,phi,rhs)
45 : !
46 : ! set phif,rhsf input in arrays which include
47 : ! virtual boundaries for phi (for all 2-d real codes)
48 : !
49 : integer nfx,nfy,i,j
50 : real(r8) :: phif(nfx,nfy),rhsf(nfx,nfy)
51 : real(r8) :: phi(0:nfx+1,0:nfy+1),rhs(nfx,nfy)
52 950 : do j=1,nfy
53 76361 : do i=1,nfx
54 75411 : phi(i,j) = phif(i,j)
55 76342 : rhs(i,j) = rhsf(i,j)
56 : end do
57 : end do
58 : !
59 : ! set virtual boundaries in phi to zero
60 : !
61 988 : do j=0,nfy+1
62 969 : phi(0,j) = 0.0_r8
63 988 : phi(nfx+1,j) = 0.0_r8
64 : end do
65 1596 : do i=0,nfx+1
66 1577 : phi(i,0) = 0.0_r8
67 1596 : phi(i,nfy+1) = 0.0_r8
68 : end do
69 19 : return
70 : end subroutine swk2
71 : !-----------------------------------------------------------------------
72 76 : subroutine trsfc2(nx,ny,phi,rhs,ncx,ncy,phic,rhsc)
73 : !
74 : ! transfer fine grid to coarse grid
75 : !
76 : integer nx,ny,ncx,ncy,i,j,ic,jc
77 : real(r8) :: phi(0:nx+1,0:ny+1),rhs(nx,ny)
78 : real(r8) :: phic(0:ncx+1,0:ncy+1),rhsc(ncx,ncy)
79 : !
80 : ! set virtual boundaries in phic to zero
81 : !
82 1159 : do jc=0,ncy+1
83 1083 : phic(0,jc) = 0.0_r8
84 1159 : phic(ncx+1,jc) = 0.0_r8
85 : end do
86 1729 : do ic=0,ncx+1
87 1653 : phic(ic,0) = 0.0_r8
88 1729 : phic(ic,ncy+1) = 0.0_r8
89 : end do
90 76 : if (ncx.lt.nx .and. ncy.lt.ny) then
91 : !
92 : ! coarsening in both x and y
93 : !
94 1007 : do jc=1,ncy
95 931 : j = jc+jc-1
96 27588 : do ic=1,ncx
97 26581 : i = ic+ic-1
98 26581 : phic(ic,jc) = phi(i,j)
99 27512 : rhsc(ic,jc) = rhs(i,j)
100 : end do
101 : end do
102 0 : else if (ncx.lt.nx .and. ncy.eq.ny) then
103 : !
104 : ! coarsening in x only
105 : !
106 0 : do jc=1,ncy
107 0 : j = jc
108 0 : do ic=1,ncx
109 0 : i = ic+ic-1
110 0 : phic(ic,jc) = phi(i,j)
111 0 : rhsc(ic,jc) = rhs(i,j)
112 : end do
113 : end do
114 : else
115 : !
116 : ! coarsening in y only
117 : !
118 0 : do jc=1,ncy
119 0 : j = jc+jc-1
120 0 : do ic=1,ncx
121 0 : i = ic
122 0 : phic(ic,jc) = phi(i,j)
123 0 : rhsc(ic,jc) = rhs(i,j)
124 : end do
125 : end do
126 : end if
127 76 : return
128 : end subroutine trsfc2
129 : !-----------------------------------------------------------------------
130 1349 : subroutine res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
131 :
132 : integer nx,ny,ncx,ncy,nxa,nxb,nyc,nyd
133 : integer i,j,ic,jc,im1,ip1,jm1,jp1,ix,jy
134 : !
135 : ! restrict fine grid residual in resf to coarse grid in rhsc
136 : ! using full weighting for all real 2d codes
137 : !
138 : real(r8) :: resf(nx,ny),rhsc(ncx,ncy)
139 : !
140 : ! set x,y coarsening integer subscript scales
141 : !
142 1349 : ix = 1
143 1349 : if (ncx.eq.nx) ix = 0
144 1349 : jy = 1
145 1349 : if (ncy.eq.ny) jy = 0
146 : !
147 : ! restrict on interior
148 : !
149 1349 : if (ncy.lt.ny .and. ncx.lt.nx) then
150 : !
151 : ! coarsening in both directions
152 : !
153 8265 : do jc=2,ncy-1
154 6916 : j = jc+jc-1
155 134349 : do ic=2,ncx-1
156 126084 : i = ic+ic-1
157 756504 : rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+ &
158 126084 : resf(i+1,j+1)+2._r8*(resf(i-1,j)+resf(i+1,j)+ &
159 1015588 : resf(i,j-1)+resf(i,j+1))+4._r8*resf(i,j))*.0625_r8
160 : end do
161 : end do
162 0 : else if (ncy.eq.ny) then
163 : !
164 : ! no coarsening in y but coarsening in x
165 : !
166 0 : do jc=2,ncy-1
167 0 : j = jc
168 0 : do ic=2,ncx-1
169 0 : i = ic+ic-1
170 0 : rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+ &
171 0 : resf(i+1,j+1)+2._r8*(resf(i-1,j)+resf(i+1,j)+ &
172 0 : resf(i,j-1)+resf(i,j+1))+4._r8*resf(i,j))*.0625_r8
173 : end do
174 : end do
175 : else
176 : !
177 : ! no coarsening in x but coarsening in y
178 : !
179 0 : do jc=2,ncy-1
180 0 : j = jc+jc-1
181 0 : do ic=2,ncx-1
182 0 : i = ic
183 0 : rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+ &
184 0 : resf(i+1,j+1)+2._r8*(resf(i-1,j)+resf(i+1,j)+ &
185 0 : resf(i,j-1)+resf(i,j+1))+4._r8*resf(i,j))*.0625_r8
186 : end do
187 : end do
188 : end if
189 : !
190 : ! set residual on boundaries
191 : !
192 4047 : do jc=1,ncy,ncy-1
193 : !
194 : ! y=yc,yd boundaries
195 : !
196 2698 : j = jc+jy*(jc-1)
197 2698 : jm1 = max0(j-1,2)
198 2698 : jp1 = min0(j+1,ny-1)
199 2698 : if (j.eq.1 .and. nyc.eq.0) jm1 = ny-1
200 2698 : if (j.eq.ny .and. nyc.eq.0) jp1 = 2
201 : !
202 : ! y=yc,yd and x=xa,xb cornors
203 : !
204 2698 : do ic=1,ncx,ncx-1
205 5396 : i = ic+ix*(ic-1)
206 5396 : im1 = max0(i-1,2)
207 5396 : ip1 = min0(i+1,nx-1)
208 5396 : if (i.eq.1 .and. nxa.eq.0) im1 = nx-1
209 5396 : if (i.eq.nx .and. nxa.eq.0) ip1 = 2
210 32376 : rhsc(ic,jc) = (resf(im1,jm1)+resf(ip1,jm1)+resf(im1,jp1)+ &
211 5396 : resf(ip1,jp1)+2._r8*(resf(im1,j)+resf(ip1,j)+ &
212 43168 : resf(i,jm1)+resf(i,jp1))+4._r8*resf(i,j))*.0625_r8
213 : end do
214 : !
215 : ! set y=yc,yd interior edges
216 : !
217 28899 : do ic=2,ncx-1
218 24852 : i = ic+ix*(ic-1)
219 149112 : rhsc(ic,jc) = (resf(i-1,jm1)+resf(i+1,jm1)+resf(i-1,jp1)+ &
220 24852 : resf(i+1,jp1)+2._r8*(resf(i-1,j)+resf(i+1,j)+ &
221 201514 : resf(i,jm1)+resf(i,jp1))+4._r8*resf(i,j))*.0625_r8
222 : end do
223 : end do
224 : !
225 : ! set x=xa,xb interior edges
226 : !
227 4047 : do ic=1,ncx,ncx-1
228 2698 : i = ic+ix*(ic-1)
229 2698 : im1 = max0(i-1,2)
230 2698 : ip1 = min0(i+1,nx-1)
231 2698 : if (i.eq.1 .and. nxa.eq.0) im1 = nx-1
232 2698 : if (i.eq.nx .and. nxa.eq.0) ip1 = 2
233 17879 : do jc=2,ncy-1
234 13832 : j = jc+jy*(jc-1)
235 82992 : rhsc(ic,jc) = (resf(im1,j-1)+resf(ip1,j-1)+resf(im1,j+1)+ &
236 13832 : resf(ip1,j+1)+2._r8*(resf(im1,j)+resf(ip1,j)+ &
237 113354 : resf(i,j-1)+resf(i,j+1))+4._r8*resf(i,j))*.0625_r8
238 : end do
239 : end do
240 : !
241 : ! set coarse grid residual zero on specified boundaries
242 : !
243 1349 : if (nxa.eq.1) then
244 0 : do jc=1,ncy
245 0 : rhsc(1,jc) = 0.0_r8
246 : end do
247 : end if
248 1349 : if (nxb.eq.1) then
249 0 : do jc=1,ncy
250 0 : rhsc(ncx,jc) = 0.0_r8
251 : end do
252 : end if
253 1349 : if (nyc.eq.1) then
254 0 : do ic=1,ncx
255 0 : rhsc(ic,1) = 0.0_r8
256 : end do
257 : end if
258 1349 : if (nyd.eq.1) then
259 16473 : do ic=1,ncx
260 16473 : rhsc(ic,ncy) = 0.0_r8
261 : end do
262 : end if
263 1349 : return
264 : end subroutine res2
265 : !-----------------------------------------------------------------------
266 : !
267 : ! prolon2 modified from rgrd2u 11/20/97
268 : !
269 1425 : subroutine prolon2(ncx,ncy,p,nx,ny,q,nxa,nxb,nyc,nyd,intpol)
270 :
271 : integer ncx,ncy,nx,ny,intpol,nxa,nxb,nyc,nyd
272 : real(r8) :: p(0:ncx+1,0:ncy+1),q(0:nx+1,0:ny+1)
273 : integer i,j,jc,ist,ifn,jst,jfn,joddst,joddfn
274 1425 : ist = 1
275 1425 : ifn = nx
276 1425 : jst = 1
277 1425 : jfn = ny
278 1425 : joddst = 1
279 1425 : joddfn = ny
280 1425 : if (nxa.eq.1) then
281 0 : ist = 2
282 : end if
283 1425 : if (nxb.eq.1) then
284 0 : ifn = nx-1
285 : end if
286 1425 : if (nyc.eq.1) then
287 0 : jst = 2
288 0 : joddst = 3
289 : end if
290 1425 : if (nyd.eq.1) then
291 1425 : jfn = ny-1
292 1425 : joddfn = ny-2
293 : end if
294 1425 : if (intpol.eq.1 .or. ncy.lt.4) then
295 : !
296 : ! linearly interpolate in y
297 : !
298 0 : if (ncy .lt. ny) then
299 : !
300 : ! ncy grid is an every other point subset of ny grid
301 : ! set odd j lines interpolating in x and then set even
302 : ! j lines by averaging odd j lines
303 : !
304 0 : do j=joddst,joddfn,2
305 0 : jc = j/2+1
306 0 : call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
307 : end do
308 0 : do j=2,jfn,2
309 0 : do i=ist,ifn
310 0 : q(i,j) = 0.5_r8*(q(i,j-1)+q(i,j+1))
311 : end do
312 : end do
313 : !
314 : ! set periodic virtual boundaries if necessary
315 : !
316 0 : if (nyc.eq.0) then
317 0 : do i=ist,ifn
318 0 : q(i,0) = q(i,ny-1)
319 0 : q(i,ny+1) = q(i,2)
320 : end do
321 : end if
322 0 : return
323 : else
324 : !
325 : ! ncy grid is equals ny grid so interpolate in x only
326 : !
327 0 : do j=jst,jfn
328 0 : jc = j
329 0 : call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
330 : end do
331 : !
332 : ! set periodic virtual boundaries if necessary
333 : !
334 0 : if (nyc.eq.0) then
335 0 : do i=ist,ifn
336 0 : q(i,0) = q(i,ny-1)
337 0 : q(i,ny+1) = q(i,2)
338 : end do
339 : end if
340 0 : return
341 : end if
342 : else
343 : !
344 : ! cubically interpolate in y
345 : !
346 1425 : if (ncy .lt. ny) then
347 : !
348 : ! set every other point of ny grid by interpolating in x
349 : !
350 10545 : do j=joddst,joddfn,2
351 9120 : jc = j/2+1
352 10545 : call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
353 : end do
354 : !
355 : ! set deep interior of ny grid using values just
356 : ! generated and symmetric cubic interpolation in y
357 : !
358 7695 : do j=4,ny-3,2
359 295165 : do i=ist,ifn
360 293740 : q(i,j)=(-q(i,j-3)+9._r8*(q(i,j-1)+q(i,j+1))-q(i,j+3))*.0625_r8
361 : end do
362 : end do
363 : !
364 : ! interpolate from q at j=2 and j=ny-1
365 : !
366 1425 : if (nyc.ne.0) then
367 : !
368 : ! asymmetric formula near nonperiodic y boundaries
369 : !
370 33250 : do i=ist,ifn
371 31825 : q(i,2)=(5._r8*q(i,1)+15._r8*q(i,3)-5._r8*q(i,5)+q(i,7))*.0625_r8
372 127300 : q(i,ny-1)=(5._r8*q(i,ny)+15._r8*q(i,ny-2)-5._r8*q(i,ny-4)+ &
373 160550 : q(i,ny-6))*.0625_r8
374 : end do
375 : else
376 : !
377 : ! periodicity in y alows symmetric formula near bndys
378 : !
379 0 : do i=ist,ifn
380 0 : q(i,2) = (-q(i,ny-2)+9._r8*(q(i,1)+q(i,3))-q(i,5))*.0625_r8
381 0 : q(i,ny-1)=(-q(i,ny-4)+9._r8*(q(i,ny-2)+q(i,ny))-q(i,3))*.0625_r8
382 0 : q(i,ny+1) = q(i,2)
383 0 : q(i,0) = q(i,ny-1)
384 : end do
385 : end if
386 1425 : return
387 : else
388 : !
389 : ! ncy grid is equals ny grid so interpolate in x only
390 : !
391 0 : do j=jst,jfn
392 0 : jc = j
393 0 : call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
394 : end do
395 : !
396 : ! set periodic virtual boundaries if necessary
397 : !
398 0 : if (nyc.eq.0) then
399 0 : do i=ist,ifn
400 0 : q(i,0) = q(i,ny-1)
401 0 : q(i,ny+1) = q(i,2)
402 : end do
403 : end if
404 0 : return
405 : end if
406 : end if
407 : end subroutine prolon2
408 : !-----------------------------------------------------------------------
409 : !
410 : ! 11/20/97 modification of rgrd1u.f for mudpack
411 : !
412 9120 : subroutine prolon1(ncx,p,nx,q,nxa,nxb,intpol)
413 :
414 : integer intpol,nxa,nxb,ncx,nx,i,ic,ist,ifn,ioddst,ioddfn
415 : real(r8) :: p(0:ncx+1),q(0:nx+1)
416 9120 : ist = 1
417 9120 : ioddst = 1
418 9120 : ifn = nx
419 9120 : ioddfn = nx
420 9120 : if (nxa.eq.1) then
421 0 : ist = 2
422 0 : ioddst = 3
423 : end if
424 9120 : if (nxb.eq.1) then
425 0 : ifn = nx-1
426 0 : ioddfn = nx-2
427 : end if
428 9120 : if (intpol.eq.1 .or. ncx.lt.4) then
429 : !
430 : ! linear interpolation in x
431 : !
432 0 : if (ncx .lt. nx) then
433 : !
434 : ! every other point of nx grid is ncx grid
435 : !
436 0 : do i=ioddst,ioddfn,2
437 0 : ic = (i+1)/2
438 0 : q(i) = p(ic)
439 : end do
440 0 : do i=2,ifn,2
441 0 : q(i) = 0.5_r8*(q(i-1)+q(i+1))
442 : end do
443 : else
444 : !
445 : ! nx grid equals ncx grid
446 : !
447 0 : do i=ist,ifn
448 0 : q(i) = p(i)
449 : end do
450 : end if
451 : !
452 : ! set virtual end points if periodic
453 : !
454 0 : if (nxa.eq.0) then
455 0 : q(0) = q(nx-1)
456 0 : q(nx+1) = q(2)
457 : end if
458 0 : return
459 : else
460 : !
461 : ! cubic interpolation in x
462 : !
463 9120 : if (ncx .lt. nx) then
464 9120 : do i=ioddst,ioddfn,2
465 180120 : ic = (i+1)/2
466 180120 : q(i) = p(ic)
467 : end do
468 : !
469 : ! set deep interior with symmetric formula
470 : !
471 9120 : do i=4,nx-3,2
472 152760 : q(i)=(-q(i-3)+9._r8*(q(i-1)+q(i+1))-q(i+3))*.0625_r8
473 : end do
474 : !
475 : ! interpolate from q at i=2 and i=nx-1
476 : !
477 9120 : if (nxa.ne.0) then
478 : !
479 : ! asymmetric formula near nonperiodic bndys
480 : !
481 0 : q(2)=(5._r8*q(1)+15._r8*q(3)-5._r8*q(5)+q(7))*.0625_r8
482 0 : q(nx-1)=(5._r8*q(nx)+15._r8*q(nx-2)-5._r8*q(nx-4)+q(nx-6))*.0625_r8
483 : else
484 : !
485 : ! periodicity in x alows symmetric formula near bndys
486 : !
487 9120 : q(2) = (-q(nx-2)+9._r8*(q(1)+q(3))-q(5))*.0625_r8
488 9120 : q(nx-1) = (-q(nx-4)+9._r8*(q(nx-2)+q(nx))-q(3))*.0625_r8
489 9120 : q(nx+1) = q(2)
490 9120 : q(0) = q(nx-1)
491 : end if
492 9120 : return
493 : else
494 : !
495 : ! ncx grid equals nx grid
496 : !
497 0 : do i=ist,ifn
498 0 : q(i) = p(i)
499 : end do
500 0 : if (nxa.eq.0) then
501 0 : q(0) = q(nx-1)
502 0 : q(nx+1) = q(2)
503 : end if
504 0 : return
505 : end if
506 : end if
507 : end subroutine prolon1
508 : !-----------------------------------------------------------------------
509 1349 : subroutine cor2(nx,ny,phif,ncx,ncy,phic,nxa,nxb,nyc,nyd,intpol,phcor)
510 : !
511 : ! add coarse grid correction in phic to fine grid approximation
512 : ! in phif using linear or cubic interpolation
513 : !
514 : integer i,j,nx,ny,ncx,ncy,nxa,nxb,nyc,nyd,intpol,ist,ifn,jst,jfn
515 : real(r8) :: phif(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
516 : real(r8) :: phcor(0:nx+1,0:ny+1)
517 21926 : do j=0,ny+1
518 753407 : do i=0,nx+1
519 752058 : phcor(i,j) = 0.0_r8
520 : end do
521 : end do
522 : !
523 : ! lift correction in phic to fine grid in phcor
524 : !
525 1349 : call prolon2(ncx,ncy,phic,nx,ny,phcor,nxa,nxb,nyc,nyd,intpol)
526 : !
527 : ! add correction in phcor to phif on nonspecified boundaries
528 : !
529 1349 : ist = 1
530 1349 : ifn = nx
531 1349 : jst = 1
532 1349 : jfn = ny
533 1349 : if (nxa.eq.1) ist = 2
534 1349 : if (nxb.eq.1) ifn = nx-1
535 1349 : if (nyc.eq.1) jst = 2
536 1349 : if (nyd.eq.1) jfn = ny-1
537 17879 : do j=jst,jfn
538 621509 : do i=ist,ifn
539 620160 : phif(i,j) = phif(i,j) + phcor(i,j)
540 : end do
541 : end do
542 : !
543 : ! add periodic points if necessary
544 : !
545 1349 : if (nyc.eq.0) then
546 0 : do i=ist,ifn
547 0 : phif(i,0) = phif(i,ny-1)
548 0 : phif(i,ny+1) = phif(i,2)
549 : end do
550 : end if
551 1349 : if (nxa.eq.0) then
552 17879 : do j=jst,jfn
553 16530 : phif(0,j) = phif(nx-1,j)
554 17879 : phif(nx+1,j) = phif(2,j)
555 : end do
556 : end if
557 1349 : end subroutine cor2
558 : !-----------------------------------------------------------------------
559 : !
560 : ! factri and factrip are:
561 : ! subroutines called by any real mudpack solver which uses line
562 : ! relaxation(s) within multigrid iteration. these subroutines do
563 : ! a vectorized factorization of m simultaneous tridiagonal systems
564 : ! of order n arising from nonperiodic or periodic discretizations
565 : !
566 95 : subroutine factri(m,n,a,b,c)
567 : !
568 : ! factor the m simultaneous tridiagonal systems of order n
569 : !
570 : integer m,n,i,j
571 : real(r8) :: a(n,m),b(n,m),c(n,m)
572 1862 : do i=2,n
573 100814 : do j=1,m
574 98952 : a(i-1,j) = a(i-1,j)/b(i-1,j)
575 100719 : b(i,j) = b(i,j)-a(i-1,j)*c(i-1,j)
576 : end do
577 : end do
578 95 : return
579 : end subroutine factri
580 : !-----------------------------------------------------------------------
581 95 : subroutine factrp(m,n,a,b,c,d,e,sum)
582 : !
583 : ! factor the m simultaneous "tridiagonal" systems of order n
584 : ! from discretized periodic system (leave out periodic n point)
585 : ! (so sweeps below only go from i=1,2,...,n-1) n > 3 is necessary
586 : !
587 : integer m,n,i,j
588 : real(r8) :: a(n,m),b(n,m),c(n,m),d(n,m),e(n,m),sum(m)
589 1957 : do j=1,m
590 1957 : d(1,j) = a(1,j)
591 : end do
592 2850 : do i=2,n-2
593 99256 : do j=1,m
594 96406 : a(i,j) = a(i,j)/b(i-1,j)
595 96406 : b(i,j) = b(i,j)-a(i,j)*c(i-1,j)
596 99161 : d(i,j) = -a(i,j)*d(i-1,j)
597 : end do
598 : end do
599 : !
600 : ! correct computation of last d element
601 : !
602 1957 : do j=1,m
603 1957 : d(n-2,j) = c(n-2,j)+d(n-2,j)
604 : end do
605 1957 : do j=1,m
606 1957 : e(1,j) = c(n-1,j)/b(1,j)
607 : end do
608 2755 : do i=2,n-3
609 97299 : do j=1,m
610 97204 : e(i,j) = -e(i-1,j)*c(i-1,j)/b(i,j)
611 : end do
612 : end do
613 1957 : do j=1,m
614 1957 : e(n-2,j) = (a(n-1,j)-e(n-3,j)*c(n-3,j))/b(n-2,j)
615 : end do
616 : !
617 : ! compute inner product (e,d) for each j in sum(j)
618 : !
619 1957 : do j=1,m
620 1957 : sum(j) = 0._r8
621 : end do
622 2945 : do i=1,n-2
623 101213 : do j=1,m
624 101118 : sum(j) = sum(j)+e(i,j)*d(i,j)
625 : end do
626 : end do
627 : !
628 : ! set last diagonal element
629 : !
630 1957 : do j=1,m
631 1957 : b(n-1,j) = b(n-1,j)-sum(j)
632 : end do
633 95 : return
634 : end subroutine factrp
635 : !-----------------------------------------------------------------------
636 0 : subroutine transp(n,amat)
637 : !
638 : ! transpose n by n real matrix
639 : !
640 : integer n,i,j
641 : real(r8) :: amat(n,n),temp
642 0 : do i=1,n-1
643 0 : do j=i+1,n
644 0 : temp = amat(i,j)
645 0 : amat(i,j) = amat(j,i)
646 0 : amat(j,i) = temp
647 : end do
648 : end do
649 0 : return
650 : end subroutine transp
651 : !-----------------------------------------------------------------------
652 0 : subroutine sgfa (a,lda,n,ipvt,info)
653 : integer lda,n,ipvt(*),info
654 : real(r8) :: a(lda,*)
655 : real(r8) :: t
656 : integer :: j,k,kp1,l,nm1
657 0 : info = 0
658 0 : nm1 = n - 1
659 0 : if (nm1 .lt. 1) go to 70
660 0 : do 60 k = 1, nm1
661 0 : kp1 = k + 1
662 0 : l = isfmax(n-k+1,a(k,k),1) + k - 1
663 0 : ipvt(k) = l
664 0 : if (a(l,k) .eq. 0.0e0_r8) go to 40
665 0 : if (l .eq. k) go to 10
666 0 : t = a(l,k)
667 0 : a(l,k) = a(k,k)
668 0 : a(k,k) = t
669 : 10 continue
670 0 : t = -1.0e0_r8/a(k,k)
671 0 : call sscl(n-k,t,a(k+1,k),1)
672 0 : do 30 j = kp1, n
673 0 : t = a(l,j)
674 0 : if (l .eq. k) go to 20
675 0 : a(l,j) = a(k,j)
676 0 : a(k,j) = t
677 : 20 continue
678 0 : call sxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
679 0 : 30 continue
680 0 : go to 50
681 : 40 continue
682 0 : info = k
683 : 50 continue
684 0 : 60 continue
685 : 70 continue
686 0 : ipvt(n) = n
687 0 : if (a(n,n) .eq. 0.0e0_r8) info = n
688 0 : return
689 : end subroutine sgfa
690 : !-----------------------------------------------------------------------
691 0 : subroutine sgsl (a,lda,n,ipvt,b,job)
692 :
693 : integer lda,n,ipvt(*),job
694 : real(r8) :: a(lda,*),b(*)
695 : real(r8) :: t
696 : integer k,kb,l,nm1
697 0 : nm1 = n - 1
698 0 : if (job .ne. 0) go to 50
699 0 : if (nm1 .lt. 1) go to 30
700 0 : do 20 k = 1, nm1
701 0 : l = ipvt(k)
702 0 : t = b(l)
703 0 : if (l .eq. k) go to 10
704 0 : b(l) = b(k)
705 0 : b(k) = t
706 : 10 continue
707 0 : call sxpy(n-k,t,a(k+1,k),1,b(k+1),1)
708 0 : 20 continue
709 : 30 continue
710 0 : do 40 kb = 1, n
711 0 : k = n + 1 - kb
712 0 : b(k) = b(k)/a(k,k)
713 0 : t = -b(k)
714 0 : call sxpy(k-1,t,a(1,k),1,b(1),1)
715 0 : 40 continue
716 0 : go to 100
717 : 50 continue
718 0 : do 60 k = 1, n
719 0 : t = sdt(k-1,a(1,k),1,b(1),1)
720 0 : b(k) = (b(k) - t)/a(k,k)
721 0 : 60 continue
722 0 : if (nm1 .lt. 1) go to 90
723 0 : do 80 kb = 1, nm1
724 0 : k = n - kb
725 0 : b(k) = b(k) + sdt(n-k,a(k+1,k),1,b(k+1),1)
726 0 : l = ipvt(k)
727 0 : if (l .eq. k) go to 70
728 0 : t = b(l)
729 0 : b(l) = b(k)
730 0 : b(k) = t
731 : 70 continue
732 0 : 80 continue
733 : 90 continue
734 : 100 continue
735 0 : return
736 : end subroutine sgsl
737 : !-----------------------------------------------------------------------
738 0 : function sdt(n,sx,incx,sy,incy) result(sdtx)
739 :
740 : real(r8), intent(in) :: sx(*),sy(*)
741 : integer, intent(in) :: n, incx, incy
742 :
743 : integer :: i,ix,iy,m,mp1
744 : real(r8) :: sdtx
745 : real(r8) :: stemp
746 :
747 0 : stemp = 0.0e0_r8
748 0 : sdtx = 0.0e0_r8
749 0 : if(n.le.0)return
750 0 : if(incx.eq.1.and.incy.eq.1)go to 20
751 0 : ix = 1
752 0 : iy = 1
753 0 : if(incx.lt.0)ix = (-n+1)*incx + 1
754 0 : if(incy.lt.0)iy = (-n+1)*incy + 1
755 0 : do 10 i = 1,n
756 0 : stemp = stemp + sx(ix)*sy(iy)
757 0 : ix = ix + incx
758 0 : iy = iy + incy
759 0 : 10 continue
760 0 : sdtx = stemp
761 0 : return
762 0 : 20 m = mod(n,5)
763 0 : if( m .eq. 0 ) go to 40
764 0 : do 30 i = 1,m
765 0 : stemp = stemp + sx(i)*sy(i)
766 0 : 30 continue
767 0 : if( n .lt. 5 ) go to 60
768 0 : 40 mp1 = m + 1
769 0 : do 50 i = mp1,n,5
770 0 : stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) + &
771 0 : sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
772 0 : 50 continue
773 : 60 sdtx = stemp
774 : return
775 : end function sdt
776 : !-----------------------------------------------------------------------
777 0 : integer function isfmax(n,sx,incx)
778 :
779 : real(r8) :: sx(*),smax
780 : integer i,incx,ix,n
781 0 : isfmax = 0
782 0 : if( n .lt. 1 ) return
783 0 : isfmax = 1
784 0 : if(n.eq.1)return
785 0 : if(incx.eq.1)go to 20
786 0 : ix = 1
787 0 : smax = abs(sx(1))
788 0 : ix = ix + incx
789 0 : do 10 i = 2,n
790 0 : if(abs(sx(ix)).le.smax) go to 5
791 0 : isfmax = i
792 0 : smax = abs(sx(ix))
793 0 : 5 ix = ix + incx
794 0 : 10 continue
795 0 : return
796 0 : 20 smax = abs(sx(1))
797 0 : do 30 i = 2,n
798 0 : if(abs(sx(i)).le.smax) go to 30
799 0 : isfmax = i
800 0 : smax = abs(sx(i))
801 0 : 30 continue
802 : return
803 : end function isfmax
804 : !-----------------------------------------------------------------------
805 0 : subroutine sxpy(n,sa,sx,incx,sy,incy)
806 :
807 : real(r8) :: sx(*),sy(*),sa
808 : integer i,incx,incy,ix,iy,m,mp1,n
809 0 : if(n.le.0)return
810 0 : if (sa .eq. 0.0_r8) return
811 0 : if(incx.eq.1.and.incy.eq.1)go to 20
812 0 : ix = 1
813 0 : iy = 1
814 0 : if(incx.lt.0)ix = (-n+1)*incx + 1
815 0 : if(incy.lt.0)iy = (-n+1)*incy + 1
816 0 : do 10 i = 1,n
817 0 : sy(iy) = sy(iy) + sa*sx(ix)
818 0 : ix = ix + incx
819 0 : iy = iy + incy
820 0 : 10 continue
821 0 : return
822 0 : 20 m = mod(n,4)
823 0 : if( m .eq. 0 ) go to 40
824 0 : do 30 i = 1,m
825 0 : sy(i) = sy(i) + sa*sx(i)
826 0 : 30 continue
827 0 : if( n .lt. 4 ) return
828 0 : 40 mp1 = m + 1
829 0 : do 50 i = mp1,n,4
830 0 : sy(i) = sy(i) + sa*sx(i)
831 0 : sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
832 0 : sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
833 0 : sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
834 0 : 50 continue
835 : return
836 : end subroutine sxpy
837 : !-----------------------------------------------------------------------
838 0 : subroutine sscl(n,sa,sx,incx)
839 :
840 : real(r8) :: sa,sx(*)
841 : integer i,incx,m,mp1,n,nincx
842 0 : if(n.le.0)return
843 0 : if(incx.eq.1)go to 20
844 0 : nincx = n*incx
845 0 : do 10 i = 1,nincx,incx
846 0 : sx(i) = sa*sx(i)
847 0 : 10 continue
848 0 : return
849 0 : 20 m = mod(n,5)
850 0 : if( m .eq. 0 ) go to 40
851 0 : do 30 i = 1,m
852 0 : sx(i) = sa*sx(i)
853 0 : 30 continue
854 0 : if( n .lt. 5 ) return
855 0 : 40 mp1 = m + 1
856 0 : do 50 i = mp1,n,5
857 0 : sx(i) = sa*sx(i)
858 0 : sx(i + 1) = sa*sx(i + 1)
859 0 : sx(i + 2) = sa*sx(i + 2)
860 0 : sx(i + 3) = sa*sx(i + 3)
861 0 : sx(i + 4) = sa*sx(i + 4)
862 0 : 50 continue
863 : return
864 : end subroutine sscl
865 : !-----------------------------------------------------------------------
866 : end module edyn_mudcom
|
