upmspace Function

  function upmspace(start,stop,p,u,ndim,base,istart,iend,mesh) result(aout)
    integer          :: p,u,ndim,pindex,uindex,pa,pb
    real(8)          :: start,stop,step,array(p*u+1),aout(ndim)
    real(8),optional :: mesh(ndim)
    real(8)          :: ustart,ustop
    integer          :: i,j
    logical,optional :: iend,istart
    logical          :: endpoint_,startpoint_,check
    real(8),optional :: base
    real(8)          :: base_
    ! real(8),optional :: mesh(p*u+1)
    if(ndim<0)stop "upmspace: N<0, abort."
    check=(ndim==(p*u)).OR.(ndim==(p*u+1))
    if(.not.check)stop "upmspace: wrong Ndim, abort."
    base_= 2.d0       ;if(present(base))base_=base
    startpoint_=.true.;if(present(istart))startpoint_=istart
    endpoint_=.true.  ;if(present(iend))endpoint_=iend
    check=startpoint_.AND.endpoint_
    pindex=1
    array(pindex) = start
    do i=1,p
       pindex=1+i*u               !index of the next p-mesh point
       pa=1+(i-1)*u               !index of the previous p-mesh point
       pb=1+i*u                   !
       array(pindex)=start + (stop-start)*base_**(-p+i) !create p-mesh
       ustart = array(pa)         !u-interval start
       ustop  = array(pb)         !u-interval stop
       step   = (ustop-ustart)/dble(u) !u-interval step
       do j=1,u-1
          uindex=pa+j    !u-mesh points
          array(uindex)=ustart + dble(j)*step
       enddo
    enddo
    if(check)then
       aout(1:ndim)=array
    elseif(.not.endpoint_)then
       aout(1:ndim)=array(1:p*u)
    elseif(.not.startpoint_)then
       aout(1:ndim)=array(2:)   
    endif
    if(present(mesh))then
       do i=1,ndim-1
          mesh(i)=abs(aout(i+1)-aout(i))
       enddo
       mesh(ndim)=abs(aout(ndim)-aout(ndim-1))
    endif
  end function upmspace