```
!BOP
! !MODULE: ice_spacecurve

module ice_spacecurve 1,1

! !DESCRIPTION:
!  This module contains routines necessary to
!  create space-filling curves.
!
! !REVISION HISTORY:
!  SVN:\$Id: \$
!
! author: John Dennis, NCAR

! !USES:
use ice_kinds_mod

implicit none

! !PUBLIC TYPES:

type, public :: factor_t
integer(int_kind)        :: numfact ! The # of factors for a value
integer(int_kind), dimension(:),pointer :: factors ! The factors
integer(int_kind), dimension(:), pointer :: used
end type

! !PUBLIC MEMBER FUNCTIONS:

public :: GenSpaceCurve,     &

public :: Factor,            &
IsFactorable,      &
PrintFactor,       &
ProdFactor,        &
MatchFactor

! !PRIVATE MEMBER FUNCTIONS:

private :: map,    		&
PeanoM, 		&
Hilbert, 		&
Cinco,  		&
GenCurve

private :: FirstFactor,      &
FindandMark

integer(int_kind), dimension(:,:), allocatable ::  &
dir,      &! direction to move along each level
ordered    ! the ordering
integer(int_kind), dimension(:), allocatable ::  &
pos        ! position along each of the axes

integer(int_kind) ::  &
maxdim,	  &! dimensionality of entire space
vcnt       ! visitation count

logical           :: verbose=.FALSE.

type (factor_t),  public,save :: fact  ! stores the factorization

!EOP
!EOC
!***********************************************************************

contains

!***********************************************************************
!BOP
! !IROUTINE: Cinco
! !INTERFACE:

recursive function Cinco(l,type,ma,md,ja,jd) result(ierr) 2,150

! !DESCRIPTION:
!  This subroutine implements a Cinco space-filling curve.
!  Cinco curves connect a Nb x Nb block of points where
!
!        Nb = 5^p
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS

integer(int_kind), intent(in) ::  &
l, 	& ! level of the space-filling curve
type,   & ! type of SFC curve
ma,     & ! Major axis [0,1]
md,  	& ! direction of major axis [-1,1]
ja,	& ! joiner axis [0,1]
jd	  ! direction of joiner axis [-1,1]

! !OUTPUT PARAMETERS

integer(int_kind) :: ierr  ! error return code

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind) :: &
lma,		&! local major axis (next level)
lmd,		&! local major direction (next level)
lja,		&! local joiner axis (next level)
ljd,		&! local joiner direction (next level)
ltype,          &! type of SFC on next level
ll		 ! next level down

logical     :: debug = .FALSE.

!-----------------------------------------------------------------------
ll = l
if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve

!--------------------------------------------------------------
!  Position [0,0]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'Cinco: After Position [0,0] ',pos
endif

!--------------------------------------------------------------
!  Position [1,0]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [1,0] ',pos
endif

!--------------------------------------------------------------
!  Position [2,0]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [2,0] ',pos
endif

!--------------------------------------------------------------
!  Position [2,1]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [2,1] ',pos
endif

!--------------------------------------------------------------
!  Position [2,2]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = ma
ljd       = -md

if(ll .gt. 1) then
if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [2,2] ',pos
endif

!--------------------------------------------------------------
!  Position [1,2]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [1,2] ',pos
endif

!--------------------------------------------------------------
!  Position [1,1]
!--------------------------------------------------------------
lma       = ma
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [1,1] ',pos
endif

!--------------------------------------------------------------
!  Position [0,1]
!--------------------------------------------------------------
lma       = ma
lmd       = -md
lja       = MOD(ma+1,maxdim)
ljd       = md

if(ll .gt. 1) then
if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [0,1] ',pos
endif

!--------------------------------------------------------------
!  Position [0,2]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [0,2] ',pos
endif

!--------------------------------------------------------------
!  Position [0,3]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,30) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [0,3] ',pos
endif

!--------------------------------------------------------------
!  Position [0,4]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,31) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [0,4] ',pos
endif

!--------------------------------------------------------------
!  Position [1,4]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = MOD(ma+1,maxdim)
ljd       = -md

if(ll .gt. 1) then
if(debug) write(*,32) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [1,4] ',pos
endif

!--------------------------------------------------------------
!  Position [1,3]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = ma
ljd       = md

if(ll .gt. 1) then
if(debug) write(*,33) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [1,3] ',pos
endif

!--------------------------------------------------------------
!  Position [2,3]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,34) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [2,3] ',pos
endif

!--------------------------------------------------------------
!  Position [2,4]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,35) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [2,4] ',pos
endif

!--------------------------------------------------------------
!  Position [3,4]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,36) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [3,4] ',pos
endif

!--------------------------------------------------------------
!  Position [4,4]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = MOD(ma+1,maxdim)
ljd       = -md

if(ll .gt. 1) then
if(debug) write(*,37) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [4,4] ',pos
endif

!--------------------------------------------------------------
!  Position [4,3]
!--------------------------------------------------------------
lma       = ma
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,38) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [4,3] ',pos
endif

!--------------------------------------------------------------
!  Position [3,3]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,39) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [3,3] ',pos
endif

!--------------------------------------------------------------
!  Position [3,2]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = ma
ljd       = md

if(ll .gt. 1) then
if(debug) write(*,40) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [3,2] ',pos
endif

!--------------------------------------------------------------
!  Position [4,2]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = MOD(ma+1,maxdim)
ljd       = -md

if(ll .gt. 1) then
if(debug) write(*,41) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [4,2] ',pos
endif

!--------------------------------------------------------------
!  Position [4,1]
!--------------------------------------------------------------
lma       = ma
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,42) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [4,1] ',pos
endif

!--------------------------------------------------------------
!  Position [3,1]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,43) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [3,1] ',pos
endif

!--------------------------------------------------------------
!  Position [3,0]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja       = ma
ljd       = md

if(ll .gt. 1) then
if(debug) write(*,44) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [3,0] ',pos
endif

!--------------------------------------------------------------
!  Position [4,0]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = ja
ljd       = jd

if(ll .gt. 1) then
if(debug) write(*,45) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'After Position [4,0] ',pos
endif

21   format('Call Cinco Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
22   format('Call Cinco Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
23   format('Call Cinco Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
24   format('Call Cinco Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
25   format('Call Cinco Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
26   format('Call Cinco Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
27   format('Call Cinco Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
28   format('Call Cinco Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
29   format('Call Cinco Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
30   format('Call Cinco Pos [0,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
31   format('Call Cinco Pos [0,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
32   format('Call Cinco Pos [1,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
33   format('Call Cinco Pos [1,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
34   format('Call Cinco Pos [2,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
35   format('Call Cinco Pos [2,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
36   format('Call Cinco Pos [3,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
37   format('Call Cinco Pos [4,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
38   format('Call Cinco Pos [4,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
39   format('Call Cinco Pos [3,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
40   format('Call Cinco Pos [3,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
41   format('Call Cinco Pos [4,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
42   format('Call Cinco Pos [4,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
43   format('Call Cinco Pos [3,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
44   format('Call Cinco Pos [3,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
45   format('Call Cinco Pos [4,0] Level ',i1,' at (',i2,',',i2,')',4(i3))

!EOC
!-----------------------------------------------------------------------

end function Cinco

!***********************************************************************
!BOP
! !IROUTINE: PeanoM
! !INTERFACE:

recursive function PeanoM(l,type,ma,md,ja,jd) result(ierr) 2,54

! !DESCRIPTION:
!  This function implements a meandering Peano
!  space-filling curve. A meandering Peano curve
!  connects a Nb x Nb block of points where
!
!        Nb = 3^p
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS

integer(int_kind), intent(in) ::  &
l,      & ! level of the space-filling curve
type,   & ! type of SFC curve
ma,     & ! Major axis [0,1]
md,     & ! direction of major axis [-1,1]
ja,     & ! joiner axis [0,1]
jd        ! direction of joiner axis [-1,1]

! !OUTPUT PARAMETERS

integer(int_kind) :: ierr  ! error return code

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind) :: &
lma,            &! local major axis (next level)
lmd,            &! local major direction (next level)
lja,            &! local joiner axis (next level)
ljd,            &! local joiner direction (next level)
ltype,          &! type of SFC on next level
ll               ! next level down

logical     :: debug = .FALSE.

!-----------------------------------------------------------------------

ll = l
if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
!--------------------------------------------------------------
!  Position [0,0]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [0,0] ',pos
endif

!--------------------------------------------------------------
! Position [0,1]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd
if(ll .gt. 1) then
if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [0,1] ',pos
endif

!--------------------------------------------------------------
! Position [0,2]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd
if(ll .gt. 1) then
if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [0,2] ',pos
endif

!--------------------------------------------------------------
! Position [1,2]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd
if(ll .gt. 1) then
if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [1,2] ',pos
endif

!--------------------------------------------------------------
! Position [2,2]
!--------------------------------------------------------------
lma        = ma
lmd        = md
lja        = MOD(lma+1,maxdim)
ljd        = -lmd

if(ll .gt. 1) then
if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [2,2] ',pos
endif

!--------------------------------------------------------------
! Position [2,1]
!--------------------------------------------------------------
lma        = ma
lmd        = -md
lja        = lma
ljd        = lmd

if(ll .gt. 1) then
if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [2,1] ',pos
endif

!--------------------------------------------------------------
! Position [1,1]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = -md
lja        = lma
ljd        = lmd

if(ll .gt. 1) then
if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [1,1] ',pos
endif

!--------------------------------------------------------------
! Position [1,0]
!--------------------------------------------------------------
lma        = MOD(ma+1,maxdim)
lmd        = -md
lja        = MOD(lma+1,maxdim)
ljd        = -lmd

if(ll .gt. 1) then
if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [1,0] ',pos
endif

!--------------------------------------------------------------
! Position [2,0]
!--------------------------------------------------------------
lma        = ma
lmd        = md
lja        = ja
ljd        = jd

if(ll .gt. 1) then
if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'PeanoM: After Position [2,0] ',pos
endif

21   format('Call PeanoM Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
22   format('Call PeanoM Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
23   format('Call PeanoM Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
24   format('Call PeanoM Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
25   format('Call PeanoM Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
26   format('Call PeanoM Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
27   format('Call PeanoM Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
28   format('Call PeanoM Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
29   format('Call PeanoM Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))

!EOC
!-----------------------------------------------------------------------

end function PeanoM

!***********************************************************************
!BOP
! !IROUTINE: Hilbert
! !INTERFACE:

recursive function Hilbert(l,type,ma,md,ja,jd) result(ierr) 2,24

! !DESCRIPTION:
!  This function implements a Hilbert space-filling curve.
!  A Hilbert curve connect a Nb x Nb block of points where
!
!        Nb = 2^p
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS

integer(int_kind), intent(in) ::  &
l,      & ! level of the space-filling curve
type,   & ! type of SFC curve
ma,     & ! Major axis [0,1]
md,     & ! direction of major axis [-1,1]
ja,     & ! joiner axis [0,1]
jd        ! direction of joiner axis [-1,1]

! !OUTPUT PARAMETERS

integer(int_kind) :: ierr  ! error return code

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind) :: &
lma,            &! local major axis (next level)
lmd,            &! local major direction (next level)
lja,            &! local joiner axis (next level)
ljd,            &! local joiner direction (next level)
ltype,          &! type of SFC on next level
ll               ! next level down

logical     :: debug = .FALSE.

!-----------------------------------------------------------------------
ll = l
if(ll .gt. 1) ltype = fact%factors(ll-1) ! Set the next type of space curve
!--------------------------------------------------------------
!  Position [0,0]
!--------------------------------------------------------------
lma       = MOD(ma+1,maxdim)
lmd       = md
lja       = lma
ljd       = lmd

if(ll .gt. 1) then
if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'Hilbert: After Position [0,0] ',pos
endif

!--------------------------------------------------------------
! Position [0,1]
!--------------------------------------------------------------
lma       = ma
lmd       = md
lja       = lma
ljd       = lmd
if(ll .gt. 1) then
if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'Hilbert: After Position [0,1] ',pos
endif

!--------------------------------------------------------------
! Position [1,1]
!--------------------------------------------------------------
lma        = ma
lmd        = md
lja        = MOD(ma+1,maxdim)
ljd        = -md

if(ll .gt. 1) then
if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'Hilbert: After Position [1,1] ',pos
endif

!--------------------------------------------------------------
! Position [1,0]
!--------------------------------------------------------------
lma        = MOD(ma+1,maxdim)
lmd        = -md
lja        = ja
ljd        = jd

if(ll .gt. 1) then
if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
ierr  = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
if(debug) call PrintCurve(ordered)
else
ierr = IncrementCurve(lja,ljd)
if(debug) print *,'Hilbert: After Position [1,0] ',pos
endif

21   format('Call Hilbert Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
22   format('Call Hilbert Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
23   format('Call Hilbert Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
24   format('Call Hilbert Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))

!EOC
!-----------------------------------------------------------------------

end function hilbert

!***********************************************************************
!BOP
! !IROUTINE: IncrementCurve
! !INTERFACE:

function IncrementCurve(ja,jd) result(ierr) 76

! !DESCRIPTION:
!   This function creates the curve which is stored in the
!   the ordered array.  The curve is implemented by
!   incrementing the curve in the direction [jd] of axis [ja].
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS:
integer(int_kind)  :: &
ja, 	&! axis to increment
jd	 ! direction along axis

! !OUTPUT PARAMETERS:
integer(int_kind) :: ierr ! error return code

!-----------------------------
! mark the newly visited point
!-----------------------------
ordered(pos(0)+1,pos(1)+1) = vcnt

!------------------------------------
! increment curve and update position
!------------------------------------
vcnt  = vcnt + 1
pos(ja) = pos(ja) + jd

ierr = 0
!EOC
!-----------------------------------------------------------------------

end function IncrementCurve

!***********************************************************************
!BOP
! !IROUTINE: log2
! !INTERFACE:

function log2( n) 2,2

! !DESCRIPTION:
!  This function calculates the log2 of its integer
!  input.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

integer(int_kind), intent(in) :: n  ! integer value to find the log2

! !OUTPUT PARAMETERS:

integer(int_kind) :: log2

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind) ::  tmp

!-------------------------------
!  Find the log2 of input value
!-------------------------------
log2 = 1
tmp =n
do while (tmp/2 .ne. 1)
tmp=tmp/2
log2=log2+1
enddo

!EOP
!-----------------------------------------------------------------------

end function log2

!***********************************************************************
!BOP
! !INTERFACE:

! !DESCRIPTION:
!  This function determines if we can create
!
! !REVISION HISTORY:
!  same as module

! !INTPUT PARAMETERS:

integer(int_kind), intent(in) ::  &
nelem,		&  ! number of blocks/elements to partition
npart              ! size of partition

! !OUTPUT PARAMETERS:
! partition is possible
!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind)   :: tmp1 ! temporary int

!-----------------------------------------------------------------------
tmp1 = nelem/npart

if(npart*tmp1 == nelem ) then
else
endif

!EOP
!-----------------------------------------------------------------------

!***********************************************************************
!BOP
! !IROUTINE: GenCurve
! !INTERFACE:

function GenCurve(l,type,ma,md,ja,jd) result(ierr) 78,6

! !DESCRIPTION:
!  This subroutine generates the next level down
!  space-filling curve
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS

integer(int_kind), intent(in) ::  &
l,      & ! level of the space-filling curve
type,   & ! type of SFC curve
ma,     & ! Major axis [0,1]
md,     & ! direction of major axis [-1,1]
ja,     & ! joiner axis [0,1]
jd        ! direction of joiner axis [-1,1]

! !OUTPUT PARAMETERS

integer(int_kind) :: ierr  ! error return code

!EOP
!BOC
!-----------------------------------------------------------------------

!-------------------------------------------------
! create the space-filling curve on the next level
!-------------------------------------------------
if(type == 2) then
ierr = Hilbert(l,type,ma,md,ja,jd)
elseif ( type == 3) then
ierr = PeanoM(l,type,ma,md,ja,jd)
elseif ( type == 5) then
ierr = Cinco(l,type,ma,md,ja,jd)
endif

!EOP
!-----------------------------------------------------------------------

end function GenCurve

function FirstFactor(fac) result(res) 2
type (factor_t) :: fac
integer :: res
logical :: found
integer (int_kind) :: i

found = .false.
i=1
do while (i<=fac%numfact .and. (.not. found))
if(fac%used(i) == 0) then
res = fac%factors(i)
found = .true.
endif
i=i+1
enddo

end function FirstFactor

function FindandMark(fac,val,f2) result(found) 4
type (factor_t) :: fac
integer :: val
logical :: found
logical :: f2
integer (int_kind) :: i

found = .false.
i=1
do while (i<=fac%numfact .and. (.not. found))
if(fac%used(i) == 0) then
if(fac%factors(i) .eq. val) then
if(f2)  then
fac%used(i) = 1
found = .true.
else if( .not. f2) then
fac%used(i) = -1
found = .false.
endif
endif
endif
i=i+1
enddo

end function FindandMark

subroutine MatchFactor(fac1,fac2,val,found) 3,6
type (factor_t) :: fac1
type (factor_t) :: fac2
integer :: val
integer :: val1
logical :: found
logical :: tmp

found = .false.

val1 = FirstFactor(fac1)
!JMD      print *,'Matchfactor: found value: ',val1
found = FindandMark(fac2,val1,.true.)
tmp = FindandMark(fac1,val1,found)
if (found) then
val = val1
else
val = 1
endif

end subroutine MatchFactor

function ProdFactor(fac) result(res) 6

type (factor_t) :: fac
integer :: res
integer (int_kind) :: i

res = 1
do i=1,fac%numfact
if(fac%used(i) <= 0) then
res = res * fac%factors(i)
endif
enddo

end function ProdFactor

subroutine PrintFactor(msg,fac)

character(len=*) :: msg
type (factor_t) :: fac
integer (int_kind) :: i

write(*,*) ' '
write(*,*) 'PrintFactor: ',msg
write(*,*) (fac%factors(i),i=1,fac%numfact)
write(*,*) (fac%used(i),i=1,fac%numfact)

end subroutine PrintFactor

!***********************************************************************
!BOP
! !IROUTINE: Factor
! !INTERFACE:

function Factor(num) result(res) 10,2

! !DESCRIPTION:
!  This function factors the input value num into a
!  product of 2,3, and 5.
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT PARAMETERS:

integer(int_kind), intent(in)  :: num  ! number to factor

! !OUTPUT PARAMETERS:

type (factor_t)     :: res

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind)   ::  &
tmp,tmp2,tmp3,tmp5   ! tempories for the factorization algorithm
integer(int_kind)   :: i,n    ! loop tempories
logical             :: found  ! logical temporary

! --------------------------------------
! Allocate allocate for max # of factors
! --------------------------------------
tmp = num
tmp2 = log2(num)
allocate(res%factors(tmp2))
allocate(res%used(tmp2))

res%used = 0
n=0

!-----------------------
!  Look for factors of 2
!-----------------------
found=.TRUE.
do while (found)
found = .FALSE.
tmp2 = tmp/2
if( tmp2*2 == tmp ) then
n = n + 1
res%factors(n) = 2
found = .TRUE.
tmp = tmp2
endif
enddo

!-----------------------
!  Look for factors of 3
!-----------------------
found=.TRUE.
do while (found)
found = .FALSE.
tmp3 = tmp/3
if( tmp3*3 == tmp ) then
n = n + 1
res%factors(n) = 3
found = .TRUE.
tmp = tmp3
endif
enddo

!-----------------------
!  Look for factors of 5
!-----------------------
found=.TRUE.
do while (found)
found = .FALSE.
tmp5 = tmp/5
if( tmp5*5 == tmp ) then
n = n + 1
res%factors(n) = 5
found = .TRUE.
tmp = tmp5
endif
enddo

!------------------------------------
! make sure that the input value
! only contains factors of 2,3,and 5
!------------------------------------
tmp=1
do i=1,n
tmp = tmp * res%factors(i)
enddo
if(tmp == num) then
res%numfact = n
else
res%numfact = -1
endif

!EOP
!---------------------------------------------------------
end function Factor

!***********************************************************************
!BOP
! !IROUTINE: IsFactorable
! !INTERFACE:

function IsFactorable(n) 6,2

! !DESCRIPTION:
!  This function determines if we can factor
!   n into 2,3,and 5.
!
! !REVISION HISTORY:
!  same as module

! !INTPUT PARAMETERS:

integer(int_kind), intent(in)  :: n  ! number to factor

! !OUTPUT PARAMETERS:
logical  :: IsFactorable  ! .TRUE. if it is factorable

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

type (factor_t)     :: fact  ! data structure to store factor information

fact = Factor(n)
if(fact%numfact .ne. -1) then
IsFactorable = .TRUE.
else
IsFactorable = .FALSE.
endif

!EOP
!-----------------------------------------------------------------------

end function IsFactorable

!***********************************************************************
!BOP
! !IROUTINE: map
! !INTERFACE:

subroutine map(l) 2,2

! !DESCRIPTION:
!   Interface routine between internal subroutines and public
!   subroutines.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:
integer(int_kind)  :: l   ! level of space-filling curve

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind)  :: &
d, 		 & ! dimension of curve only 2D is supported
type,		 & ! type of space-filling curve to start off
ierr   		   ! error return code

d = SIZE(pos)

pos=0
maxdim=d
vcnt=0

type = fact%factors(l)
ierr = GenCurve(l,type,0,1,0,1)

!EOP
!-----------------------------------------------------------------------

end subroutine map

!***********************************************************************
!BOP
! !IROUTINE: PrintCurve
! !INTERFACE:

subroutine PrintCurve(Mesh) 81,4

! !DESCRIPTION:
!  This subroutine prints the several low order
!  space-filling curves in an easy to read format
!
! !REVISION HISTORY:
!  same as module
!
! !INPUT PARAMETERS:

integer(int_kind), intent(in), target ::  Mesh(:,:) ! SFC to be printed

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------
integer(int_kind) ::  &
gridsize,	   &! order of space-filling curve
i		    ! loop temporary

!-----------------------------------------------------------------------

gridsize = SIZE(Mesh,dim=1)

if(gridsize == 2) then
write (*,*) "A Level 1 Hilbert Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,2) Mesh(1,i),Mesh(2,i)
enddo
else if(gridsize == 3) then
write (*,*) "A Level 1 Peano Meandering Curve:"
write (*,*) "---------------------------------"
do i=1,gridsize
write(*,3) Mesh(1,i),Mesh(2,i),Mesh(3,i)
enddo
else if(gridsize == 4) then
write (*,*) "A Level 2 Hilbert Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,4) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i)
enddo
else if(gridsize == 5) then
write (*,*) "A Level 1 Cinco Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,5) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i)
enddo
else if(gridsize == 6) then
write (*,*) "A Level 1 Hilbert and Level 1 Peano Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,6) Mesh(1,i),Mesh(2,i),Mesh(3,i), &
Mesh(4,i),Mesh(5,i),Mesh(6,i)
enddo
else if(gridsize == 8) then
write (*,*) "A Level 3 Hilbert Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,8) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i)
enddo
else if(gridsize == 9) then
write (*,*) "A Level 2 Peano Meandering Curve:"
write (*,*) "---------------------------------"
do i=1,gridsize
write(*,9) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
Mesh(9,i)
enddo
else if(gridsize == 10) then
write (*,*) "A Level 1 Hilbert and Level 1 Cinco Curve:"
write (*,*) "---------------------------------"
do i=1,gridsize
write(*,10) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
Mesh(9,i),Mesh(10,i)
enddo
else if(gridsize == 12) then
write (*,*) "A Level 2 Hilbert and Level 1 Peano Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,12) Mesh(1,i),Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i),Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i)
enddo
else if(gridsize == 15) then
write (*,*) "A Level 1 Peano and Level 1 Cinco Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,15) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i)
enddo
else if(gridsize == 16) then
write (*,*) "A Level 4 Hilbert Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,16) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i)
enddo
else if(gridsize == 18) then
write (*,*) "A Level 1 Hilbert and Level 2 Peano Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,18) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i)
enddo
else if(gridsize == 20) then
write (*,*) "A Level 2 Hilbert and Level 1 Cinco Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,20) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i)
enddo
else if(gridsize == 24) then
write (*,*) "A Level 3 Hilbert and Level 1 Peano Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,24) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i)
enddo
else if(gridsize == 25) then
write (*,*) "A Level 2 Cinco Curve:"
write (*,*) "------------------------------------------"
do i=1,gridsize
write(*,25) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
Mesh(25,i)
enddo
else if(gridsize == 27) then
write (*,*) "A Level 3 Peano Meandering Curve:"
write (*,*) "---------------------------------"
do i=1,gridsize
write(*,27) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
Mesh(25,i),Mesh(26,i),Mesh(27,i)
enddo
else if(gridsize == 32) then
write (*,*) "A Level 5 Hilbert Curve:"
write (*,*) "------------------------"
do i=1,gridsize
write(*,32) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i),  &
Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i),  &
Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), &
Mesh(29,i),Mesh(30,i),Mesh(31,i),Mesh(32,i)
enddo
endif
2 format('|',2(i2,'|'))
3 format('|',3(i2,'|'))
4 format('|',4(i2,'|'))
5 format('|',5(i2,'|'))
6 format('|',6(i2,'|'))
8 format('|',8(i2,'|'))
9 format('|',9(i2,'|'))
10 format('|',10(i2,'|'))
12 format('|',12(i3,'|'))
15 format('|',15(i3,'|'))
16 format('|',16(i3,'|'))
18 format('|',18(i3,'|'))
20 format('|',20(i3,'|'))
24 format('|',24(i3,'|'))
25 format('|',25(i3,'|'))
27 format('|',27(i3,'|'))
32 format('|',32(i4,'|'))

!EOC
!-----------------------------------------------------------------------

end subroutine PrintCurve

!***********************************************************************
!BOP
! !IROUTINE: GenSpaceCurve
! !INTERFACE:

subroutine  GenSpaceCurve(Mesh) 3

! !DESCRIPTION:
!  This subroutine is the public interface into the
!  space-filling curve functionality
!
! !REVISION HISTORY:
!  same as module
!

! !INPUT/OUTPUT PARAMETERS:
integer(int_kind), target,intent(inout) :: &
Mesh(:,:)		! The SFC ordering in 2D array

!EOP
!BOC
!-----------------------------------------------------------------------
!
!  local variables
!
!-----------------------------------------------------------------------

integer(int_kind) ::  &
level,   &! Level of space-filling curve
dim       ! dimension of SFC... currently limited to 2D

integer(int_kind) :: gridsize   ! number of points on a side

!-----------------------------------------------------------------------

!-----------------------------------------
!  Setup the size of the grid to traverse
!-----------------------------------------
dim = 2
gridsize = SIZE(Mesh,dim=1)
fact     = factor(gridsize)
level    = fact%numfact

if(verbose) print *,'GenSpacecurve: level is ',level
allocate(ordered(gridsize,gridsize))

!--------------------------------------------
! Setup the working arrays for the traversal
!--------------------------------------------
allocate(pos(0:dim-1))

!-----------------------------------------------------
!  The array ordered will contain the visitation order
!-----------------------------------------------------
ordered(:,:) = 0

call map(level)

Mesh(:,:) = ordered(:,:)

deallocate(pos,ordered)

!EOP
!-----------------------------------------------------------------------

end subroutine GenSpaceCurve

recursive subroutine qsort(a) 2,3

integer, intent(inout) :: a(:)
integer :: split

if(SIZE(a) > 1) then
call partition(a,split)
call qsort(a(:split-1))
call qsort(a(split:))
endif

end subroutine qsort

subroutine partition(a,marker) 1

INTEGER, INTENT(IN OUT) :: a(:)
INTEGER, INTENT(OUT) :: marker
INTEGER :: left, right, pivot, temp

pivot = (a(1) + a(size(a))) / 2  ! Average of first and last elements to prevent quadratic
left = 0                         ! behavior with sorted or reverse sorted data
right = size(a) + 1

DO WHILE (left < right)
right = right - 1
DO WHILE (a(right) > pivot)
right = right-1
END DO
left = left + 1
DO WHILE (a(left) < pivot)
left = left + 1
END DO
IF (left < right) THEN
temp = a(left)
a(left) = a(right)
a(right) = temp
END IF
END DO

IF (left == right) THEN
marker = left + 1
ELSE
marker = left
END IF

end subroutine partition

end module ice_spacecurve

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
```