! anvarslices4.f90 - slice(AN_VARIABLE, VSTRING * 4)
! Copyright (C) GFD Dennou Club, 2000.  All rights reserved

        ! 上位界面。通常は座標変数の値による指定とみなされる。
        ! 指定が "^" で始まる場合は格子番号による指定とみなされる。
        ! 下限省略時は座標変数の小さいほうが仮定される。
        ! 下限が "^" だけである場合には ^1 とみなされる。
        ! 上限省略時は座標変数の大きいほうが仮定される。
        ! 上限が "^" だけである場合には ^-1 とみなされる。
        ! 間隔省略時は連続走査が仮定される。

subroutine ANVarSliceS4(var, dimname, lower, upper, interval)
    use an_types, only: AN_VARIABLE
    use dc_string
    use netcdf_f77
    use dc_error
    use an_generic, only: Slice, get_slice
    implicit none
    type(AN_VARIABLE), intent(inout):: var
    type(VSTRING), intent(in):: dimname
    type(VSTRING), intent(in):: lower
    type(VSTRING), intent(in):: upper
    type(VSTRING), intent(in):: interval
    integer:: id, idim, istart, icount, istride, stat, idx
    character(len = NF_MAX_NAME):: dimname_try
    real, pointer:: dimvalue(:)
continue
    if (dimname == "") return
    stat = 0
    idim = 0
    nullify(dimvalue)
    do, id = 1, size(var%dimids)
        stat = nf_inq_dimname(var%fileid, var%dimids(id), dimname_try)
        if (stat /= NF_NOERR) goto 999
        if (dimname == dimname_try) then 
            idim = id
            exit
        endif
    enddo
    if (idim == 0) then
        stat = GT_EBADDIMNAME;  goto 999
    endif
    ! 次元がみつかったら
    ! 
    ! --- 下限の処理 ---
    !
    if (lower == "") then
        if (decreasing_axis()) then
            istart = -1
        else
            istart = 1
        endif
    else if (lower == "^") then
        istart = 1
    else if (vindex(lower, '^') == 1) then
        istart = stoi(extract(lower, 2))
    else
        istart = lookup_floor_index(stod(lower))
    endif
    call Slice(var, idim, start=istart)
    !
    ! --- 間隔の処理 ---
    !
    if (interval == "" .or. interval == "^") then
        istride = 1
    else if (vindex(interval, "^") == 1) then
        istride = stoi(extract(interval, 2))
    else
        istride = ceiling(stod(interval) / average_stride())
    endif
    call Slice(var, idim, stride=istride)
    !
    ! --- 上端の処理 ---
    !
    ! Slice の設定値を再確認
    call get_slice(var, idim, start=istart, stride=istride)
    if (upper == "") then
        if (decreasing_axis()) then
            idx = 1
        else
            idx = HUGE(icount)
        endif
    else if (upper == "^") then
        icount = HUGE(icount)
        call Slice(var, idim, count=icount)
        goto 120
    else if (vindex(upper, "^") == 1) then
        icount = (stoi(extract(upper, 2)) - istart) / istride
        call Slice(var, idim, count=icount)
        goto 120
    else if (upper == lower) then
        icount = 1
        call Slice(var, idim, count=icount)
        goto 120
    else
        idx = lookup_ceiling_index(stod(upper))
    endif
    if (idx > istart) then
        icount = (idx - istart) / istride + 1
    else
        ! 次元変数値が単調減少な場合
        icount = (istart - idx) / istride + 1
        istart = idx
    endif
    call Slice(var, idim, start=istart, count=icount)
    120 continue
    ! 後始末
999 continue
    if (associated(dimvalue)) deallocate(dimvalue)
    call StoreError(stat, "ANVarSliceS4", cause_s=dimname)
    return
contains

    subroutine build_value_table
        use an_generic, only: Open, Get, Close
        type(AN_VARIABLE):: dimvar
        if (associated(dimvalue)) return
        call Open(dimvar, var, idim)
        call Get(dimvar, dimvalue)
        call Close(dimvar)
    end subroutine

    logical function decreasing_axis() result(result)
        call build_value_table
        result = dimvalue(1) > dimvalue(size(dimvalue))
    end function

    real function average_stride() result(result)
        call build_value_table
        result = (dimvalue(1) - dimvalue(size(dimvalue))) / max(1, size(dimvalue) - 1)
        if (abs(result) < epsilon(result)) result = epsilon(result)
    end function

    ! 指定した値 value 以下で最大の値を与える格子番号
    ! もしみつからなければ、値の小さな端の格子番号

    integer function lookup_floor_index(value) result(result)
        double precision:: value
    continue
        if (decreasing_axis()) then
            do, result = 1, size(dimvalue)
                if (dimvalue(result) <= value) return
            enddo
            result = size(dimvalue)
        else
            do, result = size(dimvalue), 1, -1
                if (dimvalue(result) <= value) return
            enddo
            result = 1
        endif
    end function

    ! 指定した値 value 以上で最小の値を与える格子番号
    ! もしみつからなければ、値の大きな端の格子番号

    integer function lookup_ceiling_index(value) result(result)
        double precision:: value
    continue
        if (decreasing_axis()) then
            do, result = size(dimvalue), 1, -1
                if (dimvalue(result) >= value) return
            enddo
            result = 1
        else
            do, result = 1, size(dimvalue)
                if (dimvalue(result) >= value) return
            enddo
            result = size(dimvalue)
        endif
    end function

end subroutine
