FUNCTION gridgen, nn, yy, zz, $
        unit    = unit      , $
        merge   = merge     , $
        origin  = origin    , $
        edge    = edge      , $
        open    = open      , $
        one     = one       , $
        range   = range     , $
        multid  = multid    , $
        double  = double
    @compile_opt.pro            ; On error, return to caller

InitVar, multid, /key

; This is a bit sick, but I have done worse

InitVar, merge, /key

IF merge THEN BEGIN

    CASE n_params() OF
    1: r = nn
    2: BEGIN
        n1 = n_elements(nn)
        n2 = n_elements(yy)

        r = make_array(dim=[2,n1*n2], type=(IsType(nn) > IsType(yy)), /nozero)

        r[0,*] = SuperArray(nn,n2,/trail)
        r[1,*] = SuperArray(yy,n1,/lead )

        IF multid THEN r = reform(r,2,n1,n2,/overwrite) ELSE r = reform(r,/overwrite)
    END
    3: BEGIN
        n1 = n_elements(nn)
        n2 = n_elements(yy)
        n3 = n_elements(zz)

        r = make_array(dim=[3,n1*n2*n3], type=(IsType(nn) > IsType(yy) > IsType(zz)), /nozero)

        r[0,*] = SuperArray(nn,n2*n3,/trail)
        r[1,*] = SuperArray(SuperArray(yy,n1,/lead),n3,/trail)
        r[2,*] = SuperArray(zz,n1*n2,/lead)

        IF multid THEN r = reform(r,3,n1,n2,n3,/overwrite) ELSE r = reform(r,/overwrite)
    END
    ELSE: message, '/merge can only handle 1-,2- and 3-D grids'
    ENDCASE

    RETURN, r

ENDIF

IF IsType(range ,/string) THEN range  = TimeSet(range )
IF IsType(origin,/string) THEN origin = TimeSet(origin)
time_grid = IsTime(range) OR IsTime(origin)

InitVar, nn, 0          ; nn may be absent
n = long(nn)

ndim = n_elements(n)
multid AND= ndim GT 1

has_origin  = IsType(origin , /defined)
has_edge    = IsType(edge   , /defined)
has_open    = IsType(open   , /defined)
has_one     = IsType(one    , /defined)

has_range = IsType(range, /defined)
IF has_range THEN BEGIN

    range_ = range
    IF n[0] EQ 0 THEN BEGIN ; No # elements specified
        CASE time_grid OF
        0: BEGIN
            n = n_elements(range_)
            CASE 1 OF
            n EQ 1    : range_ = [0,range_]
            n/2*2 EQ n: range_ = reform(range_,2,n/2,/overwrite)
            ELSE      : message, 'illegal range spec'
            ENDCASE

            InitVar, unit, 1            
            n = round(reform((range_[1,*]-range_[0,*]))/unit)+1
            ndim = n_elements(n)
            multid AND= ndim GT 1 
        END
        1: BEGIN
            InitVar, unit, TimeUnit(/day)
            n = round(TimeOp(/subtract,range_[1],range_[0],unit))+1
            range_ = range_[0]
        END
        ENDCASE
    ENDIF

    CASE n_elements(range_) OF
    1   : BEGIN
        CASE time_grid OF
        0: range_ = [0,range_]
        1: BEGIN
            InitVar, unit, TimeUnit(/day)
            range_ = [range_,TimeOp(/add,range_,TimeSet(n-1,unit,/diff))]
        END
        ENDCASE
    END
    ndim: range_ = [intarr(1,ndim), reform(range_,1,ndim)]
    ELSE: range_ = reform(range_,2,ndim,/overwrite)
    ENDCASE

ENDIF

n_origin = n_elements(origin)  -1
n_edge   = n_elements(edge  )  -1
n_open   = n_elements(open  )  -1
n_one    = n_elements(one   )  -1
n_range  = n_elements(range_)/2-1

; The original gridgen went upto ndim=4 using the following lines:
; r[0,*] = reform( reform( gridgen1d(n[0],origin=origin[0],edge=edge[0],open=open[0],one=one[0],range=range_[0:1],double=double) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3])
; r[1,*] = reform( reform( replicate(1L,n[0]) # gridgen1d(n[1],origin=origin[1],edge=edge[1],open=open[1],one=one[1],range=range_[2:3],double=double), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3])
; r[2,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # gridgen1d(n[2],origin=origin[2],edge=edge[2],open=open[2],one=one[2],range=range_[4:5],double=double), n012 ) # replicate(1L,n[3])
; r[3,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # gridgen1d(n[3],origin=origin[3],edge=edge[3],open=open[3],one=one[3],range=range_[6:7],double=double)

; Calculate number of elements in grid
; Remember that gridgen1d adds an element if keyword edge is set.

np = lonarr(ndim)
IF has_edge THEN np[0] = n[0]+edge[0] ELSE np[0] = n[0]
FOR i=1,ndim-1 DO IF has_edge THEN np[i] = np[i-1]*(n[i]+edge[i < n_edge]) ELSE np[i] = np[i-1]*n[i]

FOR i=0,ndim-1 DO BEGIN             ; Loop over all ndim components

    FOR j=0,i-1 DO BEGIN
        grid = replicate(1L,n[j])
        IF j EQ 0 THEN tmp = grid ELSE tmp = reform( tmp # grid, np[j] )
    ENDFOR

    IF has_origin THEN origin_i = origin[  i < n_origin]
    IF has_edge   THEN edge_i   = edge  [  i < n_edge  ]
    IF has_open   THEN open_i   = open  [  i < n_open  ]
    IF has_one    THEN one_i    = one   [  i < n_one   ]
    IF has_range  THEN range_i  = range_[*,i < n_range ]

    grid = gridgen1d(n[i],unit,origin=origin_i,edge=edge_i,open=open_i,one=one_i,range=range_i,double=double)

    CASE i NE 0 OF
    0: BEGIN
        tmp = grid
        ; The first call to gridgen1d sets the type of the output array.
        CASE time_grid OF
        0: r = make_array(type=size(grid,/type), dimension=[ndim,np[ndim-1]])
        1: r = replicate( !TheTime, [ndim,np[ndim-1]] )
        ENDCASE
    ENDCASE
    1: tmp = reform(tmp # grid, np[i])
    ENDCASE

    FOR j=i+1,ndim-1 DO BEGIN
        grid = replicate(1L,n[j])
        tmp = reform( tmp # grid, np[j] )
    ENDFOR

    r[i,*] = tmp

ENDFOR

IF multid THEN r = reform(r, [ndim,n], /overwrite) ELSE r = reform(r, /overwrite)

RETURN, r  &  END