Line data Source code
1 : !BOP
2 : ! !MODULE: ice_spacecurve
3 :
4 : module ice_spacecurve
5 :
6 : ! !DESCRIPTION:
7 : ! This module contains routines necessary to
8 : ! create space-filling curves.
9 : !
10 : ! !REVISION HISTORY:
11 : !
12 : ! author: John Dennis, NCAR
13 :
14 : ! !USES:
15 : use ice_kinds_mod
16 : use ice_blocks, only: debug_blocks
17 : use ice_communicate, only: my_task, master_task
18 : use ice_exit, only: abort_ice
19 : use ice_fileunits
20 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
21 :
22 : implicit none
23 : private
24 :
25 : ! !PUBLIC TYPES:
26 :
27 : type, public :: factor_t
28 : integer(int_kind) :: numfact ! The # of factors for a value
29 : integer(int_kind), dimension(:),pointer :: factors ! The factors
30 : integer(int_kind), dimension(:), pointer :: used
31 : end type
32 :
33 : ! !PUBLIC MEMBER FUNCTIONS:
34 :
35 : public :: GenSpaceCurve
36 :
37 : public :: Factor, &
38 : IsFactorable, & ! LCOV_EXCL_LINE
39 : PrintFactor, & ! LCOV_EXCL_LINE
40 : ProdFactor, & ! LCOV_EXCL_LINE
41 : PrintCurve, & ! LCOV_EXCL_LINE
42 : MatchFactor
43 :
44 : ! !PRIVATE MEMBER FUNCTIONS:
45 :
46 : private :: map, &
47 : PeanoM, & ! LCOV_EXCL_LINE
48 : Hilbert, & ! LCOV_EXCL_LINE
49 : Cinco, & ! LCOV_EXCL_LINE
50 : GenCurve
51 :
52 : private :: FirstFactor, &
53 : FindandMark
54 :
55 : integer(int_kind), dimension(:,:), allocatable :: &
56 : ordered ! the ordering
57 : integer(int_kind), dimension(:), allocatable :: &
58 : pos ! position along each of the axes
59 :
60 : integer(int_kind) :: &
61 : maxdim, &! dimensionality of entire space ! LCOV_EXCL_LINE
62 : vcnt ! visitation count
63 :
64 : type (factor_t), public :: fact ! stores the factorization
65 :
66 : !EOP
67 : !EOC
68 : !***********************************************************************
69 :
70 : contains
71 :
72 : !***********************************************************************
73 : !BOP
74 : ! !IROUTINE: Cinco
75 : ! !INTERFACE:
76 :
77 0 : recursive function Cinco(l,type,ma,md,ja,jd) result(ierr)
78 :
79 : ! !DESCRIPTION:
80 : ! This subroutine implements a Cinco space-filling curve.
81 : ! Cinco curves connect a Nb x Nb block of points where
82 : !
83 : ! Nb = 5^p
84 : !
85 : ! !REVISION HISTORY:
86 : ! same as module
87 : !
88 :
89 :
90 : ! !INPUT PARAMETERS
91 :
92 : integer(int_kind), intent(in) :: &
93 : l, & ! level of the space-filling curve ! LCOV_EXCL_LINE
94 : type, & ! type of SFC curve ! LCOV_EXCL_LINE
95 : ma, & ! Major axis [0,1] ! LCOV_EXCL_LINE
96 : md, & ! direction of major axis [-1,1] ! LCOV_EXCL_LINE
97 : ja, & ! joiner axis [0,1] ! LCOV_EXCL_LINE
98 : jd ! direction of joiner axis [-1,1]
99 :
100 : ! !OUTPUT PARAMETERS
101 :
102 : integer(int_kind) :: ierr ! error return code
103 :
104 : !EOP
105 : !BOC
106 : !-----------------------------------------------------------------------
107 : !
108 : ! local variables
109 : !
110 : !-----------------------------------------------------------------------
111 :
112 : integer(int_kind) :: &
113 : lma, &! local major axis (next level) ! LCOV_EXCL_LINE
114 : lmd, &! local major direction (next level) ! LCOV_EXCL_LINE
115 : lja, &! local joiner axis (next level) ! LCOV_EXCL_LINE
116 : ljd, &! local joiner direction (next level) ! LCOV_EXCL_LINE
117 : ltype, &! type of SFC on next level ! LCOV_EXCL_LINE
118 : ll ! next level down
119 :
120 : character(len=*),parameter :: subname='(Cinco)'
121 :
122 : !-----------------------------------------------------------------------
123 0 : ltype = type
124 0 : ll = l
125 0 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
126 :
127 : !--------------------------------------------------------------
128 : ! Position [0,0]
129 : !--------------------------------------------------------------
130 0 : lma = ma
131 0 : lmd = md
132 0 : lja = lma
133 0 : ljd = lmd
134 :
135 0 : if(ll .gt. 1) then
136 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
137 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
138 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
139 : else
140 0 : ierr = IncrementCurve(lja,ljd)
141 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'Cinco: After Position [0,0] ',pos
142 : endif
143 :
144 : !--------------------------------------------------------------
145 : ! Position [1,0]
146 : !--------------------------------------------------------------
147 0 : lma = ma
148 0 : lmd = md
149 0 : lja = lma
150 0 : ljd = lmd
151 :
152 0 : if(ll .gt. 1) then
153 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
154 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
155 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
156 : else
157 0 : ierr = IncrementCurve(lja,ljd)
158 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [1,0] ',pos
159 : endif
160 :
161 : !--------------------------------------------------------------
162 : ! Position [2,0]
163 : !--------------------------------------------------------------
164 0 : lma = MOD(ma+1,maxdim)
165 0 : lmd = md
166 0 : lja = lma
167 0 : ljd = lmd
168 :
169 0 : if(ll .gt. 1) then
170 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
171 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
172 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
173 : else
174 0 : ierr = IncrementCurve(lja,ljd)
175 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [2,0] ',pos
176 : endif
177 :
178 : !--------------------------------------------------------------
179 : ! Position [2,1]
180 : !--------------------------------------------------------------
181 0 : lma = MOD(ma+1,maxdim)
182 0 : lmd = md
183 0 : lja = lma
184 0 : ljd = lmd
185 :
186 0 : if(ll .gt. 1) then
187 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
188 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
189 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
190 : else
191 0 : ierr = IncrementCurve(lja,ljd)
192 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [2,1] ',pos
193 : endif
194 :
195 : !--------------------------------------------------------------
196 : ! Position [2,2]
197 : !--------------------------------------------------------------
198 0 : lma = MOD(ma+1,maxdim)
199 0 : lmd = md
200 0 : lja = ma
201 0 : ljd = -md
202 :
203 0 : if(ll .gt. 1) then
204 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
205 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
206 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
207 : else
208 0 : ierr = IncrementCurve(lja,ljd)
209 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [2,2] ',pos
210 : endif
211 :
212 : !--------------------------------------------------------------
213 : ! Position [1,2]
214 : !--------------------------------------------------------------
215 0 : lma = MOD(ma+1,maxdim)
216 0 : lmd = -md
217 0 : lja = lma
218 0 : ljd = lmd
219 :
220 0 : if(ll .gt. 1) then
221 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
222 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
223 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
224 : else
225 0 : ierr = IncrementCurve(lja,ljd)
226 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [1,2] ',pos
227 : endif
228 :
229 : !--------------------------------------------------------------
230 : ! Position [1,1]
231 : !--------------------------------------------------------------
232 0 : lma = ma
233 0 : lmd = -md
234 0 : lja = lma
235 0 : ljd = lmd
236 :
237 0 : if(ll .gt. 1) then
238 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
239 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
240 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
241 : else
242 0 : ierr = IncrementCurve(lja,ljd)
243 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [1,1] ',pos
244 : endif
245 :
246 : !--------------------------------------------------------------
247 : ! Position [0,1]
248 : !--------------------------------------------------------------
249 0 : lma = ma
250 0 : lmd = -md
251 0 : lja = MOD(ma+1,maxdim)
252 0 : ljd = md
253 :
254 0 : if(ll .gt. 1) then
255 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
256 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
257 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
258 : else
259 0 : ierr = IncrementCurve(lja,ljd)
260 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [0,1] ',pos
261 : endif
262 :
263 : !--------------------------------------------------------------
264 : ! Position [0,2]
265 : !--------------------------------------------------------------
266 0 : lma = MOD(ma+1,maxdim)
267 0 : lmd = md
268 0 : lja = lma
269 0 : ljd = lmd
270 :
271 0 : if(ll .gt. 1) then
272 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
273 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
274 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
275 : else
276 0 : ierr = IncrementCurve(lja,ljd)
277 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [0,2] ',pos
278 : endif
279 :
280 : !--------------------------------------------------------------
281 : ! Position [0,3]
282 : !--------------------------------------------------------------
283 0 : lma = MOD(ma+1,maxdim)
284 0 : lmd = md
285 0 : lja = lma
286 0 : ljd = lmd
287 :
288 0 : if(ll .gt. 1) then
289 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,30) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
290 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
291 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
292 : else
293 0 : ierr = IncrementCurve(lja,ljd)
294 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [0,3] ',pos
295 : endif
296 :
297 : !--------------------------------------------------------------
298 : ! Position [0,4]
299 : !--------------------------------------------------------------
300 0 : lma = ma
301 0 : lmd = md
302 0 : lja = lma
303 0 : ljd = lmd
304 :
305 0 : if(ll .gt. 1) then
306 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,31) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
307 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
308 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
309 : else
310 0 : ierr = IncrementCurve(lja,ljd)
311 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [0,4] ',pos
312 : endif
313 :
314 : !--------------------------------------------------------------
315 : ! Position [1,4]
316 : !--------------------------------------------------------------
317 0 : lma = ma
318 0 : lmd = md
319 0 : lja = MOD(ma+1,maxdim)
320 0 : ljd = -md
321 :
322 0 : if(ll .gt. 1) then
323 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,32) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
324 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
325 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
326 : else
327 0 : ierr = IncrementCurve(lja,ljd)
328 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [1,4] ',pos
329 : endif
330 :
331 : !--------------------------------------------------------------
332 : ! Position [1,3]
333 : !--------------------------------------------------------------
334 0 : lma = MOD(ma+1,maxdim)
335 0 : lmd = -md
336 0 : lja = ma
337 0 : ljd = md
338 :
339 0 : if(ll .gt. 1) then
340 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,33) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
341 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
342 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
343 : else
344 0 : ierr = IncrementCurve(lja,ljd)
345 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [1,3] ',pos
346 : endif
347 :
348 : !--------------------------------------------------------------
349 : ! Position [2,3]
350 : !--------------------------------------------------------------
351 0 : lma = MOD(ma+1,maxdim)
352 0 : lmd = md
353 0 : lja = lma
354 0 : ljd = lmd
355 :
356 0 : if(ll .gt. 1) then
357 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,34) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
358 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
359 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
360 : else
361 0 : ierr = IncrementCurve(lja,ljd)
362 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [2,3] ',pos
363 : endif
364 :
365 : !--------------------------------------------------------------
366 : ! Position [2,4]
367 : !--------------------------------------------------------------
368 0 : lma = ma
369 0 : lmd = md
370 0 : lja = lma
371 0 : ljd = lmd
372 :
373 0 : if(ll .gt. 1) then
374 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,35) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
375 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
376 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
377 : else
378 0 : ierr = IncrementCurve(lja,ljd)
379 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [2,4] ',pos
380 : endif
381 :
382 : !--------------------------------------------------------------
383 : ! Position [3,4]
384 : !--------------------------------------------------------------
385 0 : lma = ma
386 0 : lmd = md
387 0 : lja = lma
388 0 : ljd = lmd
389 :
390 0 : if(ll .gt. 1) then
391 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,36) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
392 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
393 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
394 : else
395 0 : ierr = IncrementCurve(lja,ljd)
396 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [3,4] ',pos
397 : endif
398 :
399 : !--------------------------------------------------------------
400 : ! Position [4,4]
401 : !--------------------------------------------------------------
402 0 : lma = ma
403 0 : lmd = md
404 0 : lja = MOD(ma+1,maxdim)
405 0 : ljd = -md
406 :
407 0 : if(ll .gt. 1) then
408 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,37) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
409 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
410 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
411 : else
412 0 : ierr = IncrementCurve(lja,ljd)
413 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [4,4] ',pos
414 : endif
415 :
416 : !--------------------------------------------------------------
417 : ! Position [4,3]
418 : !--------------------------------------------------------------
419 0 : lma = ma
420 0 : lmd = -md
421 0 : lja = lma
422 0 : ljd = lmd
423 :
424 0 : if(ll .gt. 1) then
425 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,38) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
426 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
427 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
428 : else
429 0 : ierr = IncrementCurve(lja,ljd)
430 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [4,3] ',pos
431 : endif
432 :
433 : !--------------------------------------------------------------
434 : ! Position [3,3]
435 : !--------------------------------------------------------------
436 0 : lma = MOD(ma+1,maxdim)
437 0 : lmd = -md
438 0 : lja = lma
439 0 : ljd = lmd
440 :
441 0 : if(ll .gt. 1) then
442 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,39) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
443 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
444 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
445 : else
446 0 : ierr = IncrementCurve(lja,ljd)
447 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [3,3] ',pos
448 : endif
449 :
450 : !--------------------------------------------------------------
451 : ! Position [3,2]
452 : !--------------------------------------------------------------
453 0 : lma = MOD(ma+1,maxdim)
454 0 : lmd = -md
455 0 : lja = ma
456 0 : ljd = md
457 :
458 0 : if(ll .gt. 1) then
459 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,40) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
460 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
461 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
462 : else
463 0 : ierr = IncrementCurve(lja,ljd)
464 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [3,2] ',pos
465 : endif
466 :
467 : !--------------------------------------------------------------
468 : ! Position [4,2]
469 : !--------------------------------------------------------------
470 0 : lma = ma
471 0 : lmd = md
472 0 : lja = MOD(ma+1,maxdim)
473 0 : ljd = -md
474 :
475 0 : if(ll .gt. 1) then
476 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,41) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
477 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
478 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
479 : else
480 0 : ierr = IncrementCurve(lja,ljd)
481 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [4,2] ',pos
482 : endif
483 :
484 : !--------------------------------------------------------------
485 : ! Position [4,1]
486 : !--------------------------------------------------------------
487 0 : lma = ma
488 0 : lmd = -md
489 0 : lja = lma
490 0 : ljd = lmd
491 :
492 0 : if(ll .gt. 1) then
493 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,42) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
494 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
495 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
496 : else
497 0 : ierr = IncrementCurve(lja,ljd)
498 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [4,1] ',pos
499 : endif
500 :
501 : !--------------------------------------------------------------
502 : ! Position [3,1]
503 : !--------------------------------------------------------------
504 0 : lma = MOD(ma+1,maxdim)
505 0 : lmd = -md
506 0 : lja = lma
507 0 : ljd = lmd
508 :
509 0 : if(ll .gt. 1) then
510 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,43) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
511 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
512 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
513 : else
514 0 : ierr = IncrementCurve(lja,ljd)
515 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [3,1] ',pos
516 : endif
517 :
518 : !--------------------------------------------------------------
519 : ! Position [3,0]
520 : !--------------------------------------------------------------
521 0 : lma = MOD(ma+1,maxdim)
522 0 : lmd = -md
523 0 : lja = ma
524 0 : ljd = md
525 :
526 0 : if(ll .gt. 1) then
527 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,44) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
528 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
529 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
530 : else
531 0 : ierr = IncrementCurve(lja,ljd)
532 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [3,0] ',pos
533 : endif
534 :
535 : !--------------------------------------------------------------
536 : ! Position [4,0]
537 : !--------------------------------------------------------------
538 0 : lma = ma
539 0 : lmd = md
540 0 : lja = ja
541 0 : ljd = jd
542 :
543 0 : if(ll .gt. 1) then
544 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,45) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
545 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
546 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
547 : else
548 0 : ierr = IncrementCurve(lja,ljd)
549 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'After Position [4,0] ',pos
550 : endif
551 :
552 : 21 format('Call Cinco Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
553 : 22 format('Call Cinco Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
554 : 23 format('Call Cinco Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
555 : 24 format('Call Cinco Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
556 : 25 format('Call Cinco Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
557 : 26 format('Call Cinco Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
558 : 27 format('Call Cinco Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
559 : 28 format('Call Cinco Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
560 : 29 format('Call Cinco Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
561 : 30 format('Call Cinco Pos [0,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
562 : 31 format('Call Cinco Pos [0,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
563 : 32 format('Call Cinco Pos [1,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
564 : 33 format('Call Cinco Pos [1,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
565 : 34 format('Call Cinco Pos [2,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
566 : 35 format('Call Cinco Pos [2,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
567 : 36 format('Call Cinco Pos [3,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
568 : 37 format('Call Cinco Pos [4,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
569 : 38 format('Call Cinco Pos [4,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
570 : 39 format('Call Cinco Pos [3,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
571 : 40 format('Call Cinco Pos [3,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
572 : 41 format('Call Cinco Pos [4,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
573 : 42 format('Call Cinco Pos [4,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
574 : 43 format('Call Cinco Pos [3,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
575 : 44 format('Call Cinco Pos [3,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
576 : 45 format('Call Cinco Pos [4,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
577 :
578 : !EOC
579 : !-----------------------------------------------------------------------
580 :
581 0 : end function Cinco
582 :
583 : !***********************************************************************
584 : !BOP
585 : ! !IROUTINE: PeanoM
586 : ! !INTERFACE:
587 :
588 0 : recursive function PeanoM(l,type,ma,md,ja,jd) result(ierr)
589 :
590 : ! !DESCRIPTION:
591 : ! This function implements a meandering Peano
592 : ! space-filling curve. A meandering Peano curve
593 : ! connects a Nb x Nb block of points where
594 : !
595 : ! Nb = 3^p
596 : !
597 : ! !REVISION HISTORY:
598 : ! same as module
599 : !
600 :
601 : ! !INPUT PARAMETERS
602 :
603 : integer(int_kind), intent(in) :: &
604 : l, & ! level of the space-filling curve ! LCOV_EXCL_LINE
605 : type, & ! type of SFC curve ! LCOV_EXCL_LINE
606 : ma, & ! Major axis [0,1] ! LCOV_EXCL_LINE
607 : md, & ! direction of major axis [-1,1] ! LCOV_EXCL_LINE
608 : ja, & ! joiner axis [0,1] ! LCOV_EXCL_LINE
609 : jd ! direction of joiner axis [-1,1]
610 :
611 : ! !OUTPUT PARAMETERS
612 :
613 : integer(int_kind) :: ierr ! error return code
614 :
615 : !EOP
616 : !BOC
617 : !-----------------------------------------------------------------------
618 : !
619 : ! local variables
620 : !
621 : !-----------------------------------------------------------------------
622 :
623 :
624 : integer(int_kind) :: &
625 : lma, &! local major axis (next level) ! LCOV_EXCL_LINE
626 : lmd, &! local major direction (next level) ! LCOV_EXCL_LINE
627 : lja, &! local joiner axis (next level) ! LCOV_EXCL_LINE
628 : ljd, &! local joiner direction (next level) ! LCOV_EXCL_LINE
629 : ltype, &! type of SFC on next level ! LCOV_EXCL_LINE
630 : ll ! next level down
631 :
632 : character(len=*),parameter :: subname='(PeanoM)'
633 :
634 : !-----------------------------------------------------------------------
635 :
636 0 : ltype = type
637 0 : ll = l
638 0 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
639 : !--------------------------------------------------------------
640 : ! Position [0,0]
641 : !--------------------------------------------------------------
642 0 : lma = MOD(ma+1,maxdim)
643 0 : lmd = md
644 0 : lja = lma
645 0 : ljd = lmd
646 :
647 0 : if(ll .gt. 1) then
648 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
649 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
650 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
651 : else
652 0 : ierr = IncrementCurve(lja,ljd)
653 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [0,0] ',pos
654 : endif
655 :
656 :
657 : !--------------------------------------------------------------
658 : ! Position [0,1]
659 : !--------------------------------------------------------------
660 0 : lma = MOD(ma+1,maxdim)
661 0 : lmd = md
662 0 : lja = lma
663 0 : ljd = lmd
664 0 : if(ll .gt. 1) then
665 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
666 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
667 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
668 : else
669 0 : ierr = IncrementCurve(lja,ljd)
670 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [0,1] ',pos
671 : endif
672 :
673 : !--------------------------------------------------------------
674 : ! Position [0,2]
675 : !--------------------------------------------------------------
676 0 : lma = ma
677 0 : lmd = md
678 0 : lja = lma
679 0 : ljd = lmd
680 0 : if(ll .gt. 1) then
681 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
682 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
683 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
684 : else
685 0 : ierr = IncrementCurve(lja,ljd)
686 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [0,2] ',pos
687 : endif
688 :
689 : !--------------------------------------------------------------
690 : ! Position [1,2]
691 : !--------------------------------------------------------------
692 0 : lma = ma
693 0 : lmd = md
694 0 : lja = lma
695 0 : ljd = lmd
696 0 : if(ll .gt. 1) then
697 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
698 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
699 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
700 : else
701 0 : ierr = IncrementCurve(lja,ljd)
702 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [1,2] ',pos
703 : endif
704 :
705 :
706 : !--------------------------------------------------------------
707 : ! Position [2,2]
708 : !--------------------------------------------------------------
709 0 : lma = ma
710 0 : lmd = md
711 0 : lja = MOD(lma+1,maxdim)
712 0 : ljd = -lmd
713 :
714 0 : if(ll .gt. 1) then
715 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
716 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
717 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
718 : else
719 0 : ierr = IncrementCurve(lja,ljd)
720 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [2,2] ',pos
721 : endif
722 :
723 : !--------------------------------------------------------------
724 : ! Position [2,1]
725 : !--------------------------------------------------------------
726 0 : lma = ma
727 0 : lmd = -md
728 0 : lja = lma
729 0 : ljd = lmd
730 :
731 0 : if(ll .gt. 1) then
732 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
733 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
734 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
735 : else
736 0 : ierr = IncrementCurve(lja,ljd)
737 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [2,1] ',pos
738 : endif
739 :
740 : !--------------------------------------------------------------
741 : ! Position [1,1]
742 : !--------------------------------------------------------------
743 0 : lma = MOD(ma+1,maxdim)
744 0 : lmd = -md
745 0 : lja = lma
746 0 : ljd = lmd
747 :
748 0 : if(ll .gt. 1) then
749 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
750 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
751 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
752 : else
753 0 : ierr = IncrementCurve(lja,ljd)
754 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [1,1] ',pos
755 : endif
756 :
757 :
758 : !--------------------------------------------------------------
759 : ! Position [1,0]
760 : !--------------------------------------------------------------
761 0 : lma = MOD(ma+1,maxdim)
762 0 : lmd = -md
763 0 : lja = MOD(lma+1,maxdim)
764 0 : ljd = -lmd
765 :
766 0 : if(ll .gt. 1) then
767 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
768 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
769 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
770 : else
771 0 : ierr = IncrementCurve(lja,ljd)
772 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [1,0] ',pos
773 : endif
774 :
775 : !--------------------------------------------------------------
776 : ! Position [2,0]
777 : !--------------------------------------------------------------
778 0 : lma = ma
779 0 : lmd = md
780 0 : lja = ja
781 0 : ljd = jd
782 :
783 0 : if(ll .gt. 1) then
784 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
785 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
786 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
787 : else
788 0 : ierr = IncrementCurve(lja,ljd)
789 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'PeanoM: After Position [2,0] ',pos
790 : endif
791 :
792 : 21 format('Call PeanoM Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
793 : 22 format('Call PeanoM Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
794 : 23 format('Call PeanoM Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
795 : 24 format('Call PeanoM Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
796 : 25 format('Call PeanoM Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
797 : 26 format('Call PeanoM Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
798 : 27 format('Call PeanoM Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
799 : 28 format('Call PeanoM Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
800 : 29 format('Call PeanoM Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
801 :
802 : !EOC
803 : !-----------------------------------------------------------------------
804 :
805 0 : end function PeanoM
806 :
807 : !***********************************************************************
808 : !BOP
809 : ! !IROUTINE: Hilbert
810 : ! !INTERFACE:
811 :
812 0 : recursive function Hilbert(l,type,ma,md,ja,jd) result(ierr)
813 :
814 : ! !DESCRIPTION:
815 : ! This function implements a Hilbert space-filling curve.
816 : ! A Hilbert curve connect a Nb x Nb block of points where
817 : !
818 : ! Nb = 2^p
819 : !
820 : ! !REVISION HISTORY:
821 : ! same as module
822 : !
823 :
824 :
825 : ! !INPUT PARAMETERS
826 :
827 : integer(int_kind), intent(in) :: &
828 : l, & ! level of the space-filling curve ! LCOV_EXCL_LINE
829 : type, & ! type of SFC curve ! LCOV_EXCL_LINE
830 : ma, & ! Major axis [0,1] ! LCOV_EXCL_LINE
831 : md, & ! direction of major axis [-1,1] ! LCOV_EXCL_LINE
832 : ja, & ! joiner axis [0,1] ! LCOV_EXCL_LINE
833 : jd ! direction of joiner axis [-1,1]
834 :
835 : ! !OUTPUT PARAMETERS
836 :
837 : integer(int_kind) :: ierr ! error return code
838 :
839 : !EOP
840 : !BOC
841 : !-----------------------------------------------------------------------
842 : !
843 : ! local variables
844 : !
845 : !-----------------------------------------------------------------------
846 :
847 :
848 : integer(int_kind) :: &
849 : lma, &! local major axis (next level) ! LCOV_EXCL_LINE
850 : lmd, &! local major direction (next level) ! LCOV_EXCL_LINE
851 : lja, &! local joiner axis (next level) ! LCOV_EXCL_LINE
852 : ljd, &! local joiner direction (next level) ! LCOV_EXCL_LINE
853 : ltype, &! type of SFC on next level ! LCOV_EXCL_LINE
854 : ll ! next level down
855 :
856 : character(len=*),parameter :: subname='(Hilbert)'
857 :
858 : !-----------------------------------------------------------------------
859 0 : ltype = type
860 0 : ll = l
861 0 : if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
862 : !--------------------------------------------------------------
863 : ! Position [0,0]
864 : !--------------------------------------------------------------
865 0 : lma = MOD(ma+1,maxdim)
866 0 : lmd = md
867 0 : lja = lma
868 0 : ljd = lmd
869 :
870 0 : if(ll .gt. 1) then
871 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
872 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
873 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
874 : else
875 0 : ierr = IncrementCurve(lja,ljd)
876 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'Hilbert: After Position [0,0] ',pos
877 : endif
878 :
879 :
880 : !--------------------------------------------------------------
881 : ! Position [0,1]
882 : !--------------------------------------------------------------
883 0 : lma = ma
884 0 : lmd = md
885 0 : lja = lma
886 0 : ljd = lmd
887 0 : if(ll .gt. 1) then
888 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
889 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
890 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
891 : else
892 0 : ierr = IncrementCurve(lja,ljd)
893 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'Hilbert: After Position [0,1] ',pos
894 : endif
895 :
896 :
897 : !--------------------------------------------------------------
898 : ! Position [1,1]
899 : !--------------------------------------------------------------
900 0 : lma = ma
901 0 : lmd = md
902 0 : lja = MOD(ma+1,maxdim)
903 0 : ljd = -md
904 :
905 0 : if(ll .gt. 1) then
906 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
907 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
908 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
909 : else
910 0 : ierr = IncrementCurve(lja,ljd)
911 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'Hilbert: After Position [1,1] ',pos
912 : endif
913 :
914 : !--------------------------------------------------------------
915 : ! Position [1,0]
916 : !--------------------------------------------------------------
917 0 : lma = MOD(ma+1,maxdim)
918 0 : lmd = -md
919 0 : lja = ja
920 0 : ljd = jd
921 :
922 0 : if(ll .gt. 1) then
923 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
924 0 : ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
925 0 : if(debug_blocks .and. my_task==master_task) call PrintCurve(ordered)
926 : else
927 0 : ierr = IncrementCurve(lja,ljd)
928 0 : if(debug_blocks .and. my_task==master_task) write(nu_diag,*) 'Hilbert: After Position [1,0] ',pos
929 : endif
930 :
931 : 21 format('Call Hilbert Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
932 : 22 format('Call Hilbert Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
933 : 23 format('Call Hilbert Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
934 : 24 format('Call Hilbert Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
935 :
936 : !EOC
937 : !-----------------------------------------------------------------------
938 :
939 0 : end function hilbert
940 :
941 : !***********************************************************************
942 : !BOP
943 : ! !IROUTINE: IncrementCurve
944 : ! !INTERFACE:
945 :
946 0 : function IncrementCurve(ja,jd) result(ierr)
947 :
948 : ! !DESCRIPTION:
949 : ! This function creates the curve which is stored in the
950 : ! the ordered array. The curve is implemented by
951 : ! incrementing the curve in the direction [jd] of axis [ja].
952 : !
953 : ! !REVISION HISTORY:
954 : ! same as module
955 : !
956 :
957 : ! !INPUT PARAMETERS:
958 : integer(int_kind) :: &
959 : ja, &! axis to increment ! LCOV_EXCL_LINE
960 : jd ! direction along axis
961 :
962 : ! !OUTPUT PARAMETERS:
963 : integer(int_kind) :: ierr ! error return code
964 :
965 : character(len=*),parameter :: subname='(IncrementCurve)'
966 :
967 : !-----------------------------
968 : ! mark the newly visited point
969 : !-----------------------------
970 0 : ordered(pos(0)+1,pos(1)+1) = vcnt
971 :
972 : !------------------------------------
973 : ! increment curve and update position
974 : !------------------------------------
975 0 : vcnt = vcnt + 1
976 0 : pos(ja) = pos(ja) + jd
977 :
978 0 : ierr = 0
979 : !EOC
980 : !-----------------------------------------------------------------------
981 :
982 0 : end function IncrementCurve
983 :
984 : !***********************************************************************
985 : !BOP
986 : ! !IROUTINE: log2
987 : ! !INTERFACE:
988 :
989 0 : function log2( n)
990 :
991 : ! !DESCRIPTION:
992 : ! This function calculates the log2 of its integer
993 : ! input.
994 : !
995 : ! !REVISION HISTORY:
996 : ! same as module
997 :
998 : ! !INPUT PARAMETERS:
999 :
1000 : integer(int_kind), intent(in) :: n ! integer value to find the log2
1001 :
1002 : ! !OUTPUT PARAMETERS:
1003 :
1004 : integer(int_kind) :: log2
1005 :
1006 : !EOP
1007 : !BOC
1008 : !-----------------------------------------------------------------------
1009 : !
1010 : ! local variables
1011 : !
1012 : !-----------------------------------------------------------------------
1013 :
1014 : integer(int_kind) :: tmp
1015 :
1016 : character(len=*),parameter :: subname='(log2)'
1017 :
1018 : !-------------------------------
1019 : ! Find the log2 of input value
1020 : ! Abort if n < 1
1021 : !-------------------------------
1022 :
1023 0 : if (n < 1) then
1024 0 : call abort_ice (subname//'ERROR: spacecurve log2 error')
1025 :
1026 0 : elseif (n == 1) then
1027 0 : log2 = 0
1028 :
1029 : else ! n > 1
1030 0 : log2 = 1
1031 0 : tmp =n
1032 0 : do while (tmp > 1 .and. tmp/2 .ne. 1)
1033 0 : tmp=tmp/2
1034 0 : log2=log2+1
1035 : enddo
1036 : endif
1037 :
1038 : !EOP
1039 : !-----------------------------------------------------------------------
1040 :
1041 0 : end function log2
1042 :
1043 : !***********************************************************************
1044 : #ifdef UNDEPRECATE_IsLoadBalanced
1045 : !BOP
1046 : ! !IROUTINE: IsLoadBalanced
1047 : ! !INTERFACE:
1048 :
1049 : function IsLoadBalanced(nelem,npart)
1050 :
1051 : ! !DESCRIPTION:
1052 : ! This function determines if we can create
1053 : ! a perfectly load-balanced partitioning.
1054 : !
1055 : ! !REVISION HISTORY:
1056 : ! same as module
1057 :
1058 : ! !INTPUT PARAMETERS:
1059 :
1060 : integer(int_kind), intent(in) :: &
1061 : nelem, & ! number of blocks/elements to partition ! LCOV_EXCL_LINE
1062 : npart ! size of partition
1063 :
1064 : ! !OUTPUT PARAMETERS:
1065 : logical :: IsLoadBalanced ! .TRUE. if a perfectly load balanced
1066 : ! partition is possible
1067 : !EOP
1068 : !BOC
1069 : !-----------------------------------------------------------------------
1070 : !
1071 : ! local variables
1072 : !
1073 : !-----------------------------------------------------------------------
1074 :
1075 : integer(int_kind) :: tmp1 ! temporary int
1076 :
1077 : character(len=*),parameter :: subname='(IsLoadBalanced)'
1078 :
1079 : !-----------------------------------------------------------------------
1080 : tmp1 = nelem/npart
1081 :
1082 : if (npart*tmp1 == nelem ) then
1083 : IsLoadBalanced=.TRUE.
1084 : else
1085 : IsLoadBalanced=.FALSE.
1086 : endif
1087 :
1088 : !EOP
1089 : !-----------------------------------------------------------------------
1090 :
1091 : end function IsLoadBalanced
1092 : #endif
1093 : !***********************************************************************
1094 : !BOP
1095 : ! !IROUTINE: GenCurve
1096 : ! !INTERFACE:
1097 :
1098 0 : function GenCurve(l,type,ma,md,ja,jd) result(ierr)
1099 :
1100 : ! !DESCRIPTION:
1101 : ! This subroutine generates the next level down
1102 : ! space-filling curve
1103 : !
1104 : ! !REVISION HISTORY:
1105 : ! same as module
1106 : !
1107 :
1108 : ! !INPUT PARAMETERS
1109 :
1110 : integer(int_kind), intent(in) :: &
1111 : l, & ! level of the space-filling curve ! LCOV_EXCL_LINE
1112 : type, & ! type of SFC curve ! LCOV_EXCL_LINE
1113 : ma, & ! Major axis [0,1] ! LCOV_EXCL_LINE
1114 : md, & ! direction of major axis [-1,1] ! LCOV_EXCL_LINE
1115 : ja, & ! joiner axis [0,1] ! LCOV_EXCL_LINE
1116 : jd ! direction of joiner axis [-1,1]
1117 :
1118 : ! !OUTPUT PARAMETERS
1119 :
1120 : integer(int_kind) :: ierr ! error return code
1121 :
1122 : !EOP
1123 : !BOC
1124 :
1125 : logical, save :: f2=.true., f3=.true., f5=.true. ! first calls
1126 : character(len=*),parameter :: subname='(GenCurve)'
1127 :
1128 : !-----------------------------------------------------------------------
1129 :
1130 : !-------------------------------------------------
1131 : ! create the space-filling curve on the next level
1132 : !-------------------------------------------------
1133 :
1134 0 : if(type == 2) then
1135 0 : if (f2 .and. my_task == master_task) write(nu_diag,*) subname,' calling Hilbert (2)'
1136 0 : ierr = Hilbert(l,type,ma,md,ja,jd)
1137 0 : f2 = .false.
1138 0 : elseif ( type == 3) then
1139 0 : if (f3 .and. my_task == master_task) write(nu_diag,*) subname,' calling PeanoM (3)'
1140 0 : ierr = PeanoM(l,type,ma,md,ja,jd)
1141 0 : f3 = .false.
1142 0 : elseif ( type == 5) then
1143 0 : if (f5 .and. my_task == master_task) write(nu_diag,*) subname,' calling Cinco (5)'
1144 0 : ierr = Cinco(l,type,ma,md,ja,jd)
1145 0 : f5 = .false.
1146 : endif
1147 :
1148 : !EOP
1149 : !-----------------------------------------------------------------------
1150 :
1151 0 : end function GenCurve
1152 :
1153 :
1154 0 : function FirstFactor(fac) result(res)
1155 : type (factor_t) :: fac
1156 : integer :: res
1157 : logical :: found
1158 : integer (int_kind) :: i
1159 : character(len=*),parameter :: subname='(FirstFactor)'
1160 :
1161 0 : found = .false.
1162 0 : res = -1
1163 0 : i=1
1164 0 : do while (i<=fac%numfact .and. (.not. found))
1165 0 : if(fac%used(i) == 0) then
1166 0 : res = fac%factors(i)
1167 0 : found = .true.
1168 : endif
1169 0 : i=i+1
1170 : enddo
1171 :
1172 0 : end function FirstFactor
1173 :
1174 0 : function FindandMark(fac,val,f2) result(found)
1175 : type (factor_t) :: fac
1176 : integer :: val
1177 : logical :: found
1178 : logical :: f2
1179 : integer (int_kind) :: i
1180 : character(len=*),parameter :: subname='(FindandMark)'
1181 :
1182 0 : found = .false.
1183 0 : i=1
1184 0 : do while (i<=fac%numfact .and. (.not. found))
1185 0 : if(fac%used(i) == 0) then
1186 0 : if(fac%factors(i) .eq. val) then
1187 0 : if(f2) then
1188 0 : fac%used(i) = 1
1189 0 : found = .true.
1190 0 : else if( .not. f2) then
1191 0 : fac%used(i) = -1
1192 0 : found = .false.
1193 : endif
1194 : endif
1195 : endif
1196 0 : i=i+1
1197 : enddo
1198 :
1199 0 : end function FindandMark
1200 :
1201 :
1202 0 : subroutine MatchFactor(fac1,fac2,val,found)
1203 : type (factor_t) :: fac1
1204 : type (factor_t) :: fac2
1205 : integer :: val
1206 : integer :: val1
1207 : logical :: found
1208 : logical :: tmp
1209 : character(len=*),parameter :: subname='(MatchFactor)'
1210 :
1211 0 : found = .false.
1212 :
1213 0 : val1 = FirstFactor(fac1)
1214 : !JMD write(nu_diag,*)'Matchfactor: found value: ',val1
1215 0 : found = FindandMark(fac2,val1,.true.)
1216 0 : tmp = FindandMark(fac1,val1,found)
1217 0 : if (found) then
1218 0 : val = val1
1219 : else
1220 0 : val = 1
1221 : endif
1222 :
1223 0 : end subroutine MatchFactor
1224 :
1225 0 : function ProdFactor(fac) result(res)
1226 :
1227 : type (factor_t) :: fac
1228 : integer :: res
1229 : integer (int_kind) :: i
1230 : character(len=*),parameter :: subname='(ProdFactor)'
1231 :
1232 0 : res = 1
1233 0 : do i=1,fac%numfact
1234 0 : if(fac%used(i) <= 0) then
1235 0 : res = res * fac%factors(i)
1236 : endif
1237 : enddo
1238 :
1239 0 : end function ProdFactor
1240 :
1241 0 : subroutine PrintFactor(msg,fac)
1242 :
1243 :
1244 : character(len=*) :: msg
1245 : type (factor_t) :: fac
1246 : integer (int_kind) :: i
1247 : character(len=*),parameter :: subname='(PrintFactor)'
1248 :
1249 0 : write(nu_diag,*) subname,' '
1250 0 : write(nu_diag,*) subname,'msg = ',trim(msg)
1251 0 : write(nu_diag,*) subname,(fac%factors(i),i=1,fac%numfact)
1252 0 : write(nu_diag,*) subname,(fac%used(i),i=1,fac%numfact)
1253 :
1254 :
1255 0 : end subroutine PrintFactor
1256 :
1257 : !***********************************************************************
1258 : !BOP
1259 : ! !IROUTINE: Factor
1260 : ! !INTERFACE:
1261 :
1262 0 : function Factor(num) result(res)
1263 :
1264 : ! !DESCRIPTION:
1265 : ! This function factors the input value num into a
1266 : ! product of 2,3, and 5.
1267 : !
1268 : ! !REVISION HISTORY:
1269 : ! same as module
1270 : !
1271 :
1272 : ! !INPUT PARAMETERS:
1273 :
1274 : integer(int_kind), intent(in) :: num ! number to factor
1275 :
1276 : ! !OUTPUT PARAMETERS:
1277 :
1278 : type (factor_t) :: res
1279 :
1280 : !EOP
1281 : !BOC
1282 : !-----------------------------------------------------------------------
1283 : !
1284 : ! local variables
1285 : !
1286 : !-----------------------------------------------------------------------
1287 :
1288 : integer(int_kind) :: &
1289 : tmp,tmp2,tmp3,tmp5 ! tempories for the factorization algorithm
1290 : integer(int_kind) :: i,n ! loop tempories
1291 : logical :: found ! logical temporary
1292 : character(len=*),parameter :: subname='(Factor)'
1293 :
1294 : ! --------------------------------------
1295 : ! Allocate allocate for max # of factors
1296 : ! --------------------------------------
1297 0 : tmp = num
1298 0 : tmp2 = max(log2(num),1)
1299 0 : allocate(res%factors(tmp2))
1300 0 : allocate(res%used(tmp2))
1301 :
1302 0 : res%used = 0
1303 0 : n=0
1304 :
1305 :
1306 : !-----------------------
1307 : ! Look for factors of 2
1308 : !-----------------------
1309 0 : found=.TRUE.
1310 0 : do while (found)
1311 0 : found = .FALSE.
1312 0 : tmp2 = tmp/2
1313 0 : if( tmp2*2 == tmp ) then
1314 0 : n = n + 1
1315 0 : res%factors(n) = 2
1316 0 : found = .TRUE.
1317 0 : tmp = tmp2
1318 : endif
1319 : enddo
1320 :
1321 : !-----------------------
1322 : ! Look for factors of 3
1323 : !-----------------------
1324 0 : found=.TRUE.
1325 0 : do while (found)
1326 0 : found = .FALSE.
1327 0 : tmp3 = tmp/3
1328 0 : if( tmp3*3 == tmp ) then
1329 0 : n = n + 1
1330 0 : res%factors(n) = 3
1331 0 : found = .TRUE.
1332 0 : tmp = tmp3
1333 : endif
1334 : enddo
1335 :
1336 : !-----------------------
1337 : ! Look for factors of 5
1338 : !-----------------------
1339 0 : found=.TRUE.
1340 0 : do while (found)
1341 0 : found = .FALSE.
1342 0 : tmp5 = tmp/5
1343 0 : if( tmp5*5 == tmp ) then
1344 0 : n = n + 1
1345 0 : res%factors(n) = 5
1346 0 : found = .TRUE.
1347 0 : tmp = tmp5
1348 : endif
1349 : enddo
1350 :
1351 : !------------------------------------
1352 : ! make sure that the input value
1353 : ! only contains factors of 2,3,and 5
1354 : !------------------------------------
1355 0 : tmp=1
1356 0 : do i=1,n
1357 0 : tmp = tmp * res%factors(i)
1358 : enddo
1359 0 : if(tmp == num) then
1360 0 : res%numfact = n
1361 : else
1362 0 : res%numfact = -1
1363 : endif
1364 :
1365 : !EOP
1366 : !---------------------------------------------------------
1367 0 : end function Factor
1368 :
1369 : !***********************************************************************
1370 : !BOP
1371 : ! !IROUTINE: IsFactorable
1372 : ! !INTERFACE:
1373 :
1374 0 : function IsFactorable(n)
1375 :
1376 : ! !DESCRIPTION:
1377 : ! This function determines if we can factor
1378 : ! n into 2,3,and 5.
1379 : !
1380 : ! !REVISION HISTORY:
1381 : ! same as module
1382 :
1383 :
1384 : ! !INTPUT PARAMETERS:
1385 :
1386 : integer(int_kind), intent(in) :: n ! number to factor
1387 :
1388 : ! !OUTPUT PARAMETERS:
1389 : logical :: IsFactorable ! .TRUE. if it is factorable
1390 :
1391 : !EOP
1392 : !BOC
1393 : !-----------------------------------------------------------------------
1394 : !
1395 : ! local variables
1396 : !
1397 : !-----------------------------------------------------------------------
1398 :
1399 : type (factor_t) :: fact ! data structure to store factor information
1400 : character(len=*),parameter :: subname='(IsFactorable)'
1401 :
1402 0 : fact = Factor(n)
1403 0 : if(fact%numfact .ne. -1) then
1404 0 : IsFactorable = .TRUE.
1405 : else
1406 0 : IsFactorable = .FALSE.
1407 : endif
1408 :
1409 : !EOP
1410 : !-----------------------------------------------------------------------
1411 :
1412 0 : end function IsFactorable
1413 :
1414 : !***********************************************************************
1415 : !BOP
1416 : ! !IROUTINE: map
1417 : ! !INTERFACE:
1418 :
1419 0 : subroutine map(l)
1420 :
1421 : ! !DESCRIPTION:
1422 : ! Interface routine between internal subroutines and public
1423 : ! subroutines.
1424 : !
1425 : ! !REVISION HISTORY:
1426 : ! same as module
1427 :
1428 : ! !INPUT PARAMETERS:
1429 : integer(int_kind) :: l ! level of space-filling curve
1430 :
1431 :
1432 : !EOP
1433 : !BOC
1434 : !-----------------------------------------------------------------------
1435 : !
1436 : ! local variables
1437 : !
1438 : !-----------------------------------------------------------------------
1439 :
1440 : integer(int_kind) :: &
1441 : d, & ! dimension of curve only 2D is supported ! LCOV_EXCL_LINE
1442 : type, & ! type of space-filling curve to start off ! LCOV_EXCL_LINE
1443 : ierr ! error return code
1444 : character(len=*),parameter :: subname='(map)'
1445 :
1446 0 : d = SIZE(pos)
1447 :
1448 0 : pos=0
1449 0 : maxdim=d
1450 0 : vcnt=0
1451 :
1452 : ! tcx, if l is 0, then fact has no factors, just return
1453 0 : if (l == 0) return
1454 :
1455 0 : type = fact%factors(l)
1456 0 : ierr = GenCurve(l,type,0,1,0,1)
1457 :
1458 :
1459 : !EOP
1460 : !-----------------------------------------------------------------------
1461 :
1462 : end subroutine map
1463 :
1464 : !***********************************************************************
1465 : !BOP
1466 : ! !IROUTINE: PrintCurve
1467 : ! !INTERFACE:
1468 :
1469 0 : subroutine PrintCurve(Mesh)
1470 :
1471 :
1472 : ! !DESCRIPTION:
1473 : ! This subroutine prints the several low order
1474 : ! space-filling curves in an easy to read format
1475 : !
1476 : ! !REVISION HISTORY:
1477 : ! same as module
1478 : !
1479 : ! !INPUT PARAMETERS:
1480 :
1481 : integer(int_kind), intent(in), target :: Mesh(:,:) ! SFC to be printed
1482 :
1483 : !EOP
1484 : !BOC
1485 : !-----------------------------------------------------------------------
1486 : !
1487 : ! local variables
1488 : !
1489 : !-----------------------------------------------------------------------
1490 : integer(int_kind) :: &
1491 : gridsize, &! order of space-filling curve ! LCOV_EXCL_LINE
1492 : i ! loop temporary
1493 : character(len=*),parameter :: subname='(PrintCurve)'
1494 :
1495 : !-----------------------------------------------------------------------
1496 :
1497 0 : gridsize = SIZE(Mesh,dim=1)
1498 :
1499 0 : write(nu_diag,*) subname,":",gridsize
1500 :
1501 0 : if(gridsize == 2) then
1502 0 : write (nu_diag,*) "A Level 1 Hilbert Curve:"
1503 0 : write (nu_diag,*) "------------------------"
1504 0 : do i=1,gridsize
1505 0 : write(nu_diag,2) Mesh(1,i),Mesh(2,i)
1506 : enddo
1507 0 : else if(gridsize == 3) then
1508 0 : write (nu_diag,*) "A Level 1 Peano Meandering Curve:"
1509 0 : write (nu_diag,*) "---------------------------------"
1510 0 : do i=1,gridsize
1511 0 : write(nu_diag,3) Mesh(1,i),Mesh(2,i),Mesh(3,i)
1512 : enddo
1513 0 : else if(gridsize == 4) then
1514 0 : write (nu_diag,*) "A Level 2 Hilbert Curve:"
1515 0 : write (nu_diag,*) "------------------------"
1516 0 : do i=1,gridsize
1517 0 : write(nu_diag,4) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i)
1518 : enddo
1519 0 : else if(gridsize == 5) then
1520 0 : write (nu_diag,*) "A Level 1 Cinco Curve:"
1521 0 : write (nu_diag,*) "------------------------"
1522 0 : do i=1,gridsize
1523 0 : write(nu_diag,5) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i)
1524 : enddo
1525 0 : else if(gridsize == 6) then
1526 0 : write (nu_diag,*) "A Level 1 Hilbert and Level 1 Peano Curve:"
1527 0 : write (nu_diag,*) "------------------------------------------"
1528 0 : do i=1,gridsize
1529 0 : write(nu_diag,6) Mesh(1,i),Mesh(2,i),Mesh(3,i), &
1530 0 : Mesh(4,i),Mesh(5,i),Mesh(6,i)
1531 : enddo
1532 0 : else if(gridsize == 8) then
1533 0 : write (nu_diag,*) "A Level 3 Hilbert Curve:"
1534 0 : write (nu_diag,*) "------------------------"
1535 0 : do i=1,gridsize
1536 0 : write(nu_diag,8) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1537 0 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i)
1538 : enddo
1539 0 : else if(gridsize == 9) then
1540 0 : write (nu_diag,*) "A Level 2 Peano Meandering Curve:"
1541 0 : write (nu_diag,*) "---------------------------------"
1542 0 : do i=1,gridsize
1543 0 : write(nu_diag,9) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1544 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), & ! LCOV_EXCL_LINE
1545 0 : Mesh(9,i)
1546 : enddo
1547 0 : else if(gridsize == 10) then
1548 0 : write (nu_diag,*) "A Level 1 Hilbert and Level 1 Cinco Curve:"
1549 0 : write (nu_diag,*) "---------------------------------"
1550 0 : do i=1,gridsize
1551 0 : write(nu_diag,10) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1552 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), & ! LCOV_EXCL_LINE
1553 0 : Mesh(9,i),Mesh(10,i)
1554 : enddo
1555 0 : else if(gridsize == 12) then
1556 0 : write (nu_diag,*) "A Level 2 Hilbert and Level 1 Peano Curve:"
1557 0 : write (nu_diag,*) "------------------------------------------"
1558 0 : do i=1,gridsize
1559 0 : write(nu_diag,12) Mesh(1,i),Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1560 : Mesh(5,i),Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1561 0 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i)
1562 : enddo
1563 0 : else if(gridsize == 15) then
1564 0 : write (nu_diag,*) "A Level 1 Peano and Level 1 Cinco Curve:"
1565 0 : write (nu_diag,*) "------------------------"
1566 0 : do i=1,gridsize
1567 0 : write(nu_diag,15) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1568 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), & ! LCOV_EXCL_LINE
1569 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1570 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i)
1571 : enddo
1572 0 : else if(gridsize == 16) then
1573 0 : write (nu_diag,*) "A Level 4 Hilbert Curve:"
1574 0 : write (nu_diag,*) "------------------------"
1575 0 : do i=1,gridsize
1576 0 : write(nu_diag,16) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
1577 : Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), & ! LCOV_EXCL_LINE
1578 : Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1579 0 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i)
1580 : enddo
1581 0 : else if(gridsize == 18) then
1582 0 : write (nu_diag,*) "A Level 1 Hilbert and Level 2 Peano Curve:"
1583 0 : write (nu_diag,*) "------------------------------------------"
1584 0 : do i=1,gridsize
1585 0 : write(nu_diag,18) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1586 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1587 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1588 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1589 0 : Mesh(17,i),Mesh(18,i)
1590 : enddo
1591 0 : else if(gridsize == 20) then
1592 0 : write (nu_diag,*) "A Level 2 Hilbert and Level 1 Cinco Curve:"
1593 0 : write (nu_diag,*) "------------------------------------------"
1594 0 : do i=1,gridsize
1595 0 : write(nu_diag,20) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1596 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1597 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1598 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1599 0 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i)
1600 : enddo
1601 0 : else if(gridsize == 24) then
1602 0 : write (nu_diag,*) "A Level 3 Hilbert and Level 1 Peano Curve:"
1603 0 : write (nu_diag,*) "------------------------------------------"
1604 0 : do i=1,gridsize
1605 0 : write(nu_diag,24) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1606 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1607 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1608 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1609 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), & ! LCOV_EXCL_LINE
1610 0 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i)
1611 : enddo
1612 0 : else if(gridsize == 25) then
1613 0 : write (nu_diag,*) "A Level 2 Cinco Curve:"
1614 0 : write (nu_diag,*) "------------------------------------------"
1615 0 : do i=1,gridsize
1616 0 : write(nu_diag,25) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1617 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1618 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1619 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1620 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), & ! LCOV_EXCL_LINE
1621 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), & ! LCOV_EXCL_LINE
1622 0 : Mesh(25,i)
1623 : enddo
1624 0 : else if(gridsize == 27) then
1625 0 : write (nu_diag,*) "A Level 3 Peano Meandering Curve:"
1626 0 : write (nu_diag,*) "---------------------------------"
1627 0 : do i=1,gridsize
1628 0 : write(nu_diag,27) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1629 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1630 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1631 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1632 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), & ! LCOV_EXCL_LINE
1633 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), & ! LCOV_EXCL_LINE
1634 0 : Mesh(25,i),Mesh(26,i),Mesh(27,i)
1635 : enddo
1636 0 : else if(gridsize == 30) then
1637 0 : write (nu_diag,*) "A Level 1 Cinco and Level 1 Peano and Level 1 Hilbert Curve:"
1638 0 : write (nu_diag,*) "---------------------------------"
1639 0 : do i=1,gridsize
1640 0 : write(nu_diag,30) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1641 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1642 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1643 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1644 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), & ! LCOV_EXCL_LINE
1645 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), & ! LCOV_EXCL_LINE
1646 : Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), & ! LCOV_EXCL_LINE
1647 0 : Mesh(29,i),Mesh(30,i)
1648 : enddo
1649 0 : else if(gridsize == 32) then
1650 0 : write (nu_diag,*) "A Level 5 Hilbert Curve:"
1651 0 : write (nu_diag,*) "------------------------"
1652 0 : do i=1,gridsize
1653 0 : write(nu_diag,32) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
1654 : Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), & ! LCOV_EXCL_LINE
1655 : Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), & ! LCOV_EXCL_LINE
1656 : Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), & ! LCOV_EXCL_LINE
1657 : Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), & ! LCOV_EXCL_LINE
1658 : Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), & ! LCOV_EXCL_LINE
1659 : Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), & ! LCOV_EXCL_LINE
1660 0 : Mesh(29,i),Mesh(30,i),Mesh(31,i),Mesh(32,i)
1661 : enddo
1662 : endif
1663 : 2 format('|',2(i2,'|'))
1664 : 3 format('|',3(i2,'|'))
1665 : 4 format('|',4(i2,'|'))
1666 : 5 format('|',5(i2,'|'))
1667 : 6 format('|',6(i2,'|'))
1668 : 8 format('|',8(i2,'|'))
1669 : 9 format('|',9(i2,'|'))
1670 : 10 format('|',10(i2,'|'))
1671 : 12 format('|',12(i3,'|'))
1672 : 15 format('|',15(i3,'|'))
1673 : 16 format('|',16(i3,'|'))
1674 : 18 format('|',18(i3,'|'))
1675 : 20 format('|',20(i3,'|'))
1676 : 24 format('|',24(i3,'|'))
1677 : 25 format('|',25(i3,'|'))
1678 : 27 format('|',27(i3,'|'))
1679 : 30 format('|',30(i4,'|'))
1680 : 32 format('|',32(i4,'|'))
1681 :
1682 : !EOC
1683 : !-----------------------------------------------------------------------
1684 :
1685 0 : end subroutine PrintCurve
1686 :
1687 : !***********************************************************************
1688 : !BOP
1689 : ! !IROUTINE: GenSpaceCurve
1690 : ! !INTERFACE:
1691 :
1692 0 : subroutine GenSpaceCurve(Mesh)
1693 :
1694 : ! !DESCRIPTION:
1695 : ! This subroutine is the public interface into the
1696 : ! space-filling curve functionality
1697 : !
1698 : ! !REVISION HISTORY:
1699 : ! same as module
1700 : !
1701 :
1702 : ! !INPUT/OUTPUT PARAMETERS:
1703 : integer(int_kind), target,intent(inout) :: &
1704 : Mesh(:,:) ! The SFC ordering in 2D array
1705 :
1706 : !EOP
1707 : !BOC
1708 : !-----------------------------------------------------------------------
1709 : !
1710 : ! local variables
1711 : !
1712 : !-----------------------------------------------------------------------
1713 :
1714 : integer(int_kind) :: &
1715 : level, &! Level of space-filling curve ! LCOV_EXCL_LINE
1716 : dim ! dimension of SFC... currently limited to 2D
1717 :
1718 : integer(int_kind) :: gridsize ! number of points on a side
1719 :
1720 : character(len=*),parameter :: subname='(GenSpaceCurve)'
1721 :
1722 : !-----------------------------------------------------------------------
1723 :
1724 : !-----------------------------------------
1725 : ! Setup the size of the grid to traverse
1726 : !-----------------------------------------
1727 0 : dim = 2
1728 0 : gridsize = SIZE(Mesh,dim=1)
1729 0 : fact = factor(gridsize)
1730 0 : level = fact%numfact
1731 :
1732 0 : if (debug_blocks .and. my_task==master_task .and. my_task==master_task) then
1733 0 : write(nu_diag,*) subname,' dim,size = ',dim,gridsize
1734 0 : write(nu_diag,*) subname,' numfact = ',level
1735 0 : call printfactor(subname,fact)
1736 0 : call flush_fileunit(nu_diag)
1737 : endif
1738 :
1739 0 : allocate(ordered(gridsize,gridsize))
1740 :
1741 : !--------------------------------------------
1742 : ! Setup the working arrays for the traversal
1743 : !--------------------------------------------
1744 0 : allocate(pos(0:dim-1))
1745 :
1746 : !-----------------------------------------------------
1747 : ! The array ordered will contain the visitation order
1748 : !-----------------------------------------------------
1749 0 : ordered(:,:) = 0
1750 :
1751 0 : call map(level)
1752 :
1753 0 : Mesh(:,:) = ordered(:,:)
1754 :
1755 0 : deallocate(pos,ordered)
1756 :
1757 0 : end subroutine GenSpaceCurve
1758 :
1759 : !EOC
1760 : !-----------------------------------------------------------------------
1761 :
1762 0 : end module ice_spacecurve
1763 :
1764 : !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
