! Copyright (C) GFD Dennou Club, 2000.  All rights reserved.
! テスト用の浅水波モデル

program gtshallow

    use shallow_water
    use gtool_history

    implicit none

    integer:: nx = 9
    integer:: ny = 9
    integer:: ktmax = 1000
    real(double):: dt = 0.01_double
    real(double):: dx = 1.0_double
    real(double):: H = 1.0_double
    real(double):: g = 9.80665_double
    real(double):: f = 0.0_double
    character(len = 78):: filename = "gtshallow.nc"
    character(len = 78):: user = ""
    character(len = 78):: init(5)
    data init/"P", 4 * ' '/
    namelist /config/ nx, ny, dt, ktmax, g, H, f, init, filename, user

    type(shallow):: past, current, future
    integer:: kt, i 

continue

    write(*, *) '--- please enter configuration ---'
    write(*, config)
    read(*, config)
    write(*, *) '--- run config ---'
    write(*, config)

    call HistoryCreate(file=filename, title="shallow water simulation", &
        & source="gtshallow version 1", institution=user, &
        & dims=(/"x", "y", "t"/), &
        & dimsizes=(/nx, ny, 0/), &
        & longnames=(/"eastward length", "northward length", "time"/), &
        & units=(/"m", "m", "s"/), &
        & origin=0.0, interval=real(dt))

    call HistoryPut("x", (/(real(i * dx), i = 1, nx)/))
    call HistoryPut("y", (/(real(i * dx), i = 1, ny)/))

    call nullify(past)
    call allocate(current, nx, ny, H)
    call nullify(future)
    call modify(current, init)

    call HistoryAddVariable(varname="u", &
        & dims=(/"x", "y", "t"/), &
        & longname="eastward velocity", units="m/s")
    call HistoryAddVariable(varname="v", &
        & dims=(/"x", "y", "t"/), &
        & longname="northward velocity", &
        & units="m/s")
    call HistoryAddVariable(varname="h", &
        & dims=(/"x", "y", "t"/), &
        & longname="depth", units="m")

    future = current + dt * tendency(current)
    call history_put_all(current)
    call history_put_all(future)

    do, kt = 1, ktmax
        call deallocate(past)
        past = current
        current = future
        future = past + dt * tendency(current)
        call history_put_all(future)
        print "(9(' ', f7.3))", current%h(4, :) - H
    enddo

    call HistoryClose

contains

    subroutine history_put_all(stat)
        type(shallow):: stat
        call HistoryPut("u", real(stat%u))
        call HistoryPut("v", real(stat%v))
        call HistoryPut("h", real(stat%h))
    end subroutine

    type(shallow) function tendency(stat) result(result)
        type(shallow):: stat
        call allocate(result, nx, ny, temporary=.TRUE.)
        !
        ! 微分量を中間変数に確保、算出
        !
        allocate(stat%Ux(nx, ny), stat%Vx(nx, ny), stat%Hx(nx, ny), &
            stat%Uy(nx, ny), stat%Vy(nx, ny), stat%Hy(nx, ny))
        !
        stat%ux = (CShift(stat%U, -1, dim=1) - CShift(stat%U, 1, dim=1)) / 2.0 / dx
        stat%vx = (CShift(stat%V, -1, dim=1) - CShift(stat%V, 1, dim=1)) / 2.0 / dx
        stat%hx = (CShift(stat%H, -1, dim=1) - CShift(stat%H, 1, dim=1)) / 2.0 / dx
        stat%uy = (CShift(stat%U, -1, dim=2) - CShift(stat%U, 1, dim=2)) / 2.0 / dx
        stat%vy = (CShift(stat%V, -1, dim=2) - CShift(stat%V, 1, dim=2)) / 2.0 / dx
        stat%hy = (CShift(stat%H, -1, dim=2) - CShift(stat%H, 1, dim=2)) / 2.0 / dx
        !
        ! u_t = - u u_x - v u_y + f v - g h_x
        !
        result%U = - stat%U * stat%Ux - stat%V * stat%Uy &
            & + f * stat%V - g * stat%Hx
        !
        ! v_t = - u v_x - v v_y - f u - g h_y
        !
        result%V = - stat%V * stat%Vx - stat%V * stat%Vy &
            & - f * stat%U - g * stat%Hy
        !
        ! h_t = - (u h)_x - (v h)_y
        !
        result%H = - stat%U * stat%Hx - stat%Ux * stat%H &
            - stat%V * stat%Hy - stat%Vy * stat%H
        !
        ! 微分量バッファの後始末
        !
        deallocate(stat%Ux, stat%Vx, stat%Hx, stat%Uy, stat%Vy, stat%Hy)
    end function

    subroutine modify(stat, init)
        type(SHALLOW):: stat
        character(len = *), intent(in):: init(:)
        integer:: n, i, j
    continue
        do, n = 1, size(init)
            if (init(n) == "P") then
                i = size(stat, 1) / 2
                j = size(stat, 2) / 2
                stat%h(i, j) = stat%h(i, j) + h * 0.2
            else if (init(n) == "") then
                continue
            else
                print *, "modify: unknown specification <", trim(init(n)), ">"
            endif
        enddo
    end subroutine

end program 
