source: trunk/pop/source/spacecurve_mod.F90 @ 141

Revision 141, 50.6 KB checked in by maltrud, 5 years ago (diff)

LANL merge with NCAR revision 17212. LANL version removed several features,
including overflows, run-time bathymetry generation, and some of the
coupling infrastructure. Also, some of the tavg calls are different, which
will be addressed in a later update. High frequency pointwise and regional
diagnostics have also been removed, to be replaced in the future by the "moorings"
infrastructure. Diagnostics using "lat_aux_grid" have also been removed.

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