! Copyright (C) GFD Dennou Club, 2000.  All rights reserved.

subroutine GTVarSlice(var, dimord, start, count, stride)
    use gtdata_types, only: gt_variable
    use gtdata_internal, only: query_growable
    use gt_map, only: map_lookup, map_set, gt_dimmap, gtvar_dump
    use dc_error, only: nf_enotvar, StoreError
    use dc_trace, only: beginsub, endsub, message
    implicit none
    type(gt_variable), intent(in):: var
    integer, intent(in):: dimord
    integer, intent(in), optional:: start
    integer, intent(in), optional:: count
    integer, intent(in), optional:: stride
    type(gt_dimmap), allocatable:: map(:)
    integer:: vid, maxindex, maxcount, nd, stat
    logical:: growable_dimension

    call beginsub('gtvarslice', 'var%%mapid=%d dimord=%d', &
        & i=(/var%mapid, dimord/))
    call gtvar_dump(var)
    call map_lookup(var, vid=vid, ndims=nd)
    if (vid < 0) then
        call StoreError(nf_enotvar, "GTVarSlice")
    endif

    if (vid > 0) then
        call query_growable(vid, growable_dimension)
    else
        growable_dimension = .false.
    endif

    if (nd == 0) goto 999
    allocate(map(nd))    
    call map_lookup(var, map=map)

    if (dimord <= 0 .or. dimord > size(map)) goto 998

    call message('map(dimord): originally start=%d count=%d stride=%d', &
        & i=(/map(dimord)%start, map(dimord)%count, map(dimord)%stride/))
    if (.not. growable_dimension) then
        maxindex = map(dimord)%allcount
        call message('maxindex=%d', i=(/maxindex/))
    endif

    if (present(start)) then
        if (start < 0) then
            map(dimord)%start = max(1, maxindex + 1 + start)
        else if (growable_dimension) then
            map(dimord)%start = max(1, start)
        else
            map(dimord)%start = min(maxindex, max(1, start))
        endif
        call message('start=%d (%d specified)', i=(/map(dimord)%start, start/))
    endif

    if (present(stride)) then
        map(dimord)%stride = stride
        if (stride == 0) map(dimord)%stride = 1
        call message('stride=%d (%d specified)', &
            & i=(/map(dimord)%stride, stride/))
    endif

    if (present(count)) then
        map(dimord)%count = abs(count)
        if (count == 0) map(dimord)%count = 1
        call message('count=%d (%d specified)', &
            & i=(/map(dimord)%count, count/))
    endif

    if (.not. growable_dimension) then
        maxcount = 1 + (maxindex - map(dimord)%start) / map(dimord)%stride
        map(dimord)%count = max(1, min(maxcount, map(dimord)%count))
        call message('count=%d ', i=(/map(dimord)%count/))
    endif
    call map_set(var, map, stat)
    if (stat /= 0) goto 998

    call endsub('gtvarslice')
    deallocate(map)
    return

998 continue
    deallocate(map)
999 continue
    call endsub('gtvarslice', 'err skipped')
end subroutine

! 変数 var に string によるマップ操作を行う。
! string はコンマで区切られた変換指定の列である。
! 変換指定は領域設定と次元順序変換のどちらかである。
! 領域設定は英数字で始まるもので、<dim>=<lower>, <dim>=<lower>:<upper>,
!   <dim>=<lower>:<upper>:<stride> のような形式である。
!   ここで、dim は次元番号または次元名であり、<lower>, <upper>
!   は ^ を前置した座標即値または格子番号である。
!   <stride> は格子数である。
! (未実装) 次元順序変換は = で始まるもので、
!   IGN:<dim>=<pos>
! の形態をとる。

subroutine GTVarSliceC(var, string, err)
    use gtdata_types, only: gt_variable
    use dc_trace, only: beginsub, endsub
    use dc_url, only: gt_comma
    use gt_map, only: gtvar_dump
    type(gt_variable), intent(in):: var
    character(len = *), intent(in):: string
    logical, intent(out):: err
    integer:: is, ie
continue
    call beginsub('gtvarslicec', 'var=%d lim=<%c>', i=(/var%mapid/), c1=string)
    call gtvar_dump(var)
    ! コンマで区切って解釈
    is = 1
    do
        ie = index(string(is: ), gt_comma)
        if (ie == 0) exit
        call limit_one(string(is: is+ie-2))
        is = is + ie
        if (is > len(string)) exit
    enddo
    call limit_one(string(is: ))
    err = .false.
    call endsub('gtvarslicec')
    return
contains

    subroutine limit_one(string)
        use dc_url, only: gt_equal
        use dc_string, only: strieq, stoi
        use gtdata_generic, only: del_dim, dimname_to_dimord
        character(len = *), intent(in):: string
        integer:: equal, dimord
        integer:: start, count, stride
        logical:: myerr

        if (string == '') return

        if (strieq(string(1:4), "IGN:")) then
            ! 隠蔽型指定子 ign:<dim> または ign:<dim>=<start>
            equal = index(string, gt_equal)
            if (equal == 0) then
                start = 1
            else
                start = stoi(string(equal+1: ), default=1)
            endif
            dimord = dimname_to_dimord(var, string(5: equal-1))
            call gtvarslice(var, dimord, start, 1, 1)
            call del_dim(var, dimord, myerr)
            return
        endif

        ! 限定型指定子 <dim>=<start>:<finish>:<stride>
        !
        equal = index(string, gt_equal)
        if (equal == 0) return
        dimord = dimname_to_dimord(var, string(1: equal-1))
        if (dimord <= 0) return
        !
        call region_spec(dimord, string(equal+1: ), start, count, stride)
        call gtvarslice(var, dimord, start, count, stride)
    end subroutine

        !
        ! 範囲指定の = のあとを : で区切ってマップにいれる
        !
    subroutine region_spec(dimord, string, start, count, stride)
        use dc_types, only: token
        use dc_string, only: index_ofs, stoi
        use dc_url, only: gt_circumflex, gt_colon
        use gt_map, only: dimrange
        integer, intent(in):: dimord
        integer, intent(out):: start, count, stride
        character(len = *), intent(in):: string
        integer:: colon, prev_colon, finish, dimlo, dimhi
        character(len = token):: val(3)
    continue
        colon = index(string, gt_colon)
        if (colon == 0) then
            ! コロンがない場合は上下端に同じ値
            val(1) = string(1: )
            val(2) = val(1)
            val(3) = ""
        else
            val(1) = string(1: colon - 1)
            prev_colon = colon
            colon = index_ofs(string, colon + 1, gt_colon)
            if (colon > 0) then
                val(2) = string(prev_colon + 1: colon - 1)
                val(3) = string(colon + 1: )
            else
                val(2) = string(prev_colon + 1: )
                val(3) = ""
            endif
        endif
        if (val(3) == "") val(3) = "^1"

        if (val(1)(1:1) == gt_circumflex) then
            start = stoi(val(1)(2: ))
        else if (val(1) == val(2)) then
            start = nint(value_to_index(dimord, val(1)))
        else
            start = floor(value_to_index(dimord, val(1)))
        endif
        if (val(2) == val(1)) then
            finish = start
        else if (val(2)(1:1) == gt_circumflex) then
            finish = stoi(val(2)(2: ))
        else
            finish = ceiling(value_to_index(dimord, val(2)))
        endif

        call dimrange(var, dimord, dimlo, dimhi)
        start = min(max(dimlo, start), dimhi)
        finish = min(max(dimlo, finish), dimhi)
        count = abs(finish - start) + 1

        if (val(3)(1:1) == gt_circumflex) then
            stride = stoi(val(3)(2: ))
        else
            stride = stoi(val(3))
        endif
        stride = sign(stride, finish - start)
    end subroutine

    real function value_to_index(dimord, value) result(result)
        use gtdata_types, only: gt_variable
        use gtdata_generic, only: get, open, close
        use dc_string, only: stod
        use dc_trace, only: beginsub, endsub, message
        integer, intent(in):: dimord
        character(len = *), intent(in):: value
        type(gt_variable):: axisvar
        real, pointer:: axisval(:)
        real:: val
        integer:: i

        call beginsub('value_to_index', 'var=%d dimord=%d value=%c', &
            & i=(/var%mapid, dimord/), c1=trim(value))

        call Open(axisvar, var, dimord, count_compact=.true.)
        call Get(axisvar, axisval)
        call Close(axisvar)
        if (.not. associated(axisval)) then
            result = -1.0
            return
        else if (size(axisval) < 2) then
            result = 1.0
            goto 900
        endif

        val = stod(value)

        ! call message('value=%f axis=(/%*r/)', r=(/val, axisval(:)/), &
        !    & n=(/size(axisval)/))

        do, i = 1, size(axisval) - 1
            if (axisval(i + 1) == axisval(i)) then
                result = real(i) + 0.5
                goto 900
            endif
            result = i + (val - axisval(i)) / (axisval(i + 1) - axisval(i))
            if (result <= (i + 1)) goto 900
        enddo

    900 continue
        call endsub('value_to_index', '(%c) = %r', &
            & c1=trim(value), r=(/result/))
        deallocate(axisval)
    end function

end subroutine
