! Copyright (C) GFD Dennou Club, 2000.  All rights reserved.
! gtunary - 変数への単項演算子適用

! 書式
!   gtunary [オプション=値 ...] 演算子 [...]
!
! 入力は input オプションまたは gtool.nc@default
! 出力は output オプションまたは gtool.nc@default
! 演算子は任意個指定可. デフォルトで「なにもしない=コピー」.

program gtunary
    use gtool
    use netcdf_f77, only: NF_EINVAL
    use dc_trace, only: setdebug
    implicit none
    character(len = STRING):: output, input
    character(len = STRING):: arg, optname, optvalue
    type(GT_VARIABLE):: in, out
    integer:: argind, argc, argstart, siz, stat, chunk
    double precision, allocatable:: buf(:)
!--
    input = "gtool.nc@"
    output = "gtool.nc@default"
    argc = GtArgCount()
    argind = 0
    if (argc == 0) then
        print *, "usage: gtunary [in=<input>] [out=<output>] operator ..."
        stop
    endif
    do, while (argind <= argc)
        argind = argind + 1
        call GtArgGet(argind, arg)
        if (.not. GtOptionForm(arg, optname, optvalue)) exit
        if (optname == "in") then
            input = optvalue
        else if (optname == "out") then
            output = optvalue
        else if (optname == "-debug") then
            call setdebug
        else
            call StoreError(NF_EINVAL, 'unknown option')
        endif
    enddo
    argstart = argind
    chunk = 1
    call Open(in, input)
    call Create(out, url=output, copyfrom=in, copyvalue=.FALSE.)
    call Slice(in)
    call Slice(out, compatible=in)
    call Inquire(out, size=siz)
    allocate(buf(siz))
    do
        call Get(in, buf, siz)
        do, argind = argstart, argc
            call GtArgGet(argind, arg)
            call put_line(trim(arg))
            call UnaryOperation(buf, trim(arg))
        enddo
        call Put(out, buf, siz)
        call Slice_Next(in, stat=stat)
        if (stat == GT_ENOMOREDIMS) exit
        call Slice_Next(out) 
        call printf(fmt="chunk %d:", i=(/chunk/))
        chunk = chunk + 1
    enddo
    call Close(in)
    call Close(out)
    stop
contains

    subroutine UnaryOperation(buf, cmd)
        implicit none
        double precision, intent(inout):: buf(:)
        character(len = *), intent(in):: cmd
        double precision:: param
        logical, save:: lfirst = .TRUE.
        double precision, parameter:: undef = -9.99e5
        if (strieq(cmd, "NEGATE")) then
            buf(:) = -buf(:)
        else if (strieq(cmd, "SQRT")) then
            where (buf(:) >= 0.0)
                buf(:) = sqrt(buf(:))
            elsewhere
                buf(:) = undef
            endwhere
        else if (strieq(cmd, "LOG")) then
            where (buf(:) >= 0.0)
                buf(:) = log(buf(:))
            elsewhere
                buf(:) = undef
            endwhere
        else if (strieq(cmd, "LOG10")) then
            where (buf(:) >= 0.0)
                buf(:) = log10(buf(:))
            elsewhere
                buf(:) = undef
            endwhere
        else if (cmd(1:2) == '-=' .or. strieq(cmd(1:2), 's:')) then
            read(cmd(3: ), *) param
            buf(:) = buf(:) - param
        else if (cmd(1:2) == '+=' .or. strieq(cmd(1:2), 'a:')) then
            read(cmd(3: ), *) param
            buf(:) = buf(:) + param
        else if (cmd(1:2) == '*=' .or. strieq(cmd(1:2), 'm:')) then
            read(cmd(3: ), *) param
            buf(:) = buf(:) * param
        else if (cmd(1:2) == '/=' .or. strieq(cmd(1:2), 'd:')) then
            read(cmd(3: ), *) param
            if (abs(param) > tiny(0.0d0)) then
                buf(:) = buf(:) / param
            endif
        else if (cmd(1:2) == '^=' .or. strieq(cmd(1:2), 'p:')) then
            read(cmd(3: ), *) param
            buf(:) = param ** buf(:)
        else if (lfirst) then
            lfirst = .FALSE.
            call printf(fmt="unknown operation <%c>", c1=cmd)
        end if
    end subroutine

end program