! vi: set sw=4 ts=8:

module gt3conv_dimtable
    use gtool
    use gt3read
    use gt3conv_vartable
    implicit none

    private
    public DisposeTables, Dimension
    public BuildDimensions, BuildTime, PutDimensionData, PutTimeData
    public DIMTABLE_ENTRY, TIMETABLE_ENTRY, dims_table, time_table

    type DIMTABLE_ENTRY
	type(NC_DIMENSION):: dim
	type(NC_VARIABLE):: var
	character(len = 30):: name
	character(len = 16):: item
	integer:: lower, upper
	real, pointer:: values(:)
    end type

    type(DIMTABLE_ENTRY), save, pointer:: dims_table(:)

    type TIMETABLE_ENTRY
	integer, pointer:: values(:)
	type(NC_DIMENSION):: dim
	type(NC_VARIABLE):: var
    end type
    type(TIMETABLE_ENTRY), save, pointer:: time_table(:)

    interface Dimension;  module procedure LookupDimension;  end interface

contains

    type(NC_DIMENSION) function LookupDimension(item) result(result)
	character(len = 16), intent(in):: item
	integer:: idim
    continue
	do, idim = 1, size(dims_table)
	    if (dims_table(idim)%item == item) then
		result = dims_table(idim)%dim;  return
	    endif
	enddo
	result = NC_DIMENSION_ERROR()
    end function

    subroutine DisposeTables
	integer:: i
    continue
	call DisposeVariableTable
	do, i = 1, size(dims_table)
	    deallocate(dims_table(i)%values)
	enddo
	deallocate(dims_table)
	do, i = 1, size(time_table)
	    deallocate(time_table(i)%values)
	enddo
	deallocate(time_table)
	nullify(time_table)
    end subroutine

    type(VSTRING) function AxisFilename(item, weight) result(result)
	character(len = 16), intent(in):: item
	logical, intent(in), optional:: weight
	character(len = *), parameter:: GTAXDIR_L = '/var/spool/gt3-dcl5/'
	character(len = *), parameter:: GTAXDIR_W = 'C:/DENNOU/GT3/'
	character(len = 8):: prefix
	logical:: exist
    continue
	prefix = 'GTAXLOC.'
	if (present(weight)) then
	    if (weight) prefix = 'GTAXWGT.'
	endif
	result = prefix // trim(item)
	inquire(file=char(result), exist=exist)
	if (exist) return
	result = GTAXDIR_L // prefix // trim(item)
	inquire(file=char(result), exist=exist)
	if (exist) return
	result = GTAXDIR_W // prefix // trim(item)
	inquire(file=char(result), exist=exist)
	if (exist) return
	result = GtoolGetenv('GTAXDIR') // prefix // trim(item)
	inquire(file=char(result), exist=exist)
    end function

    subroutine PutDimensionData
	integer:: idim
    continue
	do, idim = 1, size(dims_table)
	    call Put(dims_table(idim)%var, dims_table(idim)%values)
	enddo
    end subroutine

    subroutine ReadAxisFile(dim)
	type(DIMTABLE_ENTRY), intent(inout):: dim
	type(GT3_FILE):: axisfile
	type(GT3_HEADER):: header
	real, pointer:: buffer(:, :, :)
	integer:: lb, ub
	type(NC_ATTRIBUTE):: attr
    continue
	call Open(axisfile, char(AxisFilename(dim%item)))
	call Get(axisfile, header, buffer)
	call Close(axisfile)

	lb = lbound(buffer, 1)
	ub = ubound(buffer, 1)
	if (dim%lower < lb) then
	    print *, 'gt3conv#ReadAxisFile: requested lower bound ', &
		& dim%lower, ' is less than axis file value ', lb
	    stop
	endif
	if (dim%upper > ub) then
	    print *, 'gt3conv#ReadAxisFile: requested upper bound ', &
		& dim%upper, ' is greater than axis file value ', ub
	    stop
	endif
	allocate(dim%values(lb: ub))
	dim%values(:) = buffer(lb: ub, lbound(buffer, 2), lbound(buffer, 3))

	if (header%dataset(1: 1) == 'C') then
	    attr = Attribute(dim%var, 'topology')
	    attr = 'circular'
	    call Dispose(attr)
	endif
	if (header%unit /= '') then
	    attr = Attribute(dim%var, 'unit')
	    ! fake
	    attr = UnitString(header)
	    call Dispose(attr)
	endif
    end subroutine

    type(VSTRING) function Dimname3to4(item, usedname) result(result)
	character(len = *), intent(in):: item
	character(len = *), intent(in):: usedname(:)
	character(len = len(item)):: buffer
	type(VSTRING):: name, grids, suffix
	integer:: a, b, c
    continue
	buffer = adjustl(item)
	a = scan(buffer, '0123456789')
	name = buffer(1: a-1)
	b = verify(buffer(a: ), '0123456789')
	if (b == 0) b = len(buffer)
	grids = buffer(a: b)
	suffix = trim(buffer(b + 1: ))

	if (name == 'CSIG' .and. suffix == '.M') then
	    result = 'sigma_bound'
	else if (name == 'CSIG') then
	    result = 'sigma'
	else if (name == 'GLON') then
	    result = 'lon'
	else if (name == 'GGLA') then
	    result = 'lat'
	else
	    result = item
	endif
	if (all(usedname(:) /= char(result))) return

	result = result // grids
	if (all(usedname(:) /= char(result))) return

	name = result
	do, c = 1, huge(c)
	    result = name // itos(c)
	    if (all(usedname(:) /= char(result))) return
	enddo
	stop 'gt3conv#Dimname3to4'
    end function

    type(VSTRING) function UnitString(header)
	type(GT3_HEADER), intent(in):: header
    continue
	if (header%unit == 'deg') then
	    if (header%item(1: 4) == 'GLON') then
		UnitString = 'degree_east'
	    else
		UnitString = 'degree_north'
	    endif
	else
	    UnitString = header%unit
	endif
    end function

    subroutine BuildDimensions(file) 
	type(NC_FILE), intent(inout):: file
	type(DIMTABLE_ENTRY), pointer:: dims_tmp(:)
	integer:: iv, ia, idim, num_dims
    continue
	allocate(dims_table(3))

	num_dims = 0
	! すべての変数について軸 ITEM 名を検査し、重複しないものを登録
	do, iv = 1, size(vars_table)
	    do, ia = 1, 3
		call StoreDimension(vars_table(iv)%header, ia)
	    enddo
	enddo

	! 次元表の大きさを最低限に縮小
	if (size(dims_table) > num_dims) then
	    dims_tmp => dims_table
	    allocate(dims_table(1: num_dims))
	    dims_table(:) = dims_tmp(1: num_dims)
	    deallocate(dims_tmp)
	endif

	! 登録されたすべての次元を作成
	do, idim = 1, num_dims
	    call MakeDimension(file, dims_table(idim), &
		& dims_table(1: idim-1)%name)
	enddo

    contains

	subroutine MakeDimension(file, dim, usedname)
	    type(NC_FILE), intent(inout):: file
	    type(DIMTABLE_ENTRY), intent(inout):: dim
	    character(len = *), intent(in):: usedname(:)
	    integer:: length
	continue
	    dim%name = Dimname3to4(dim%item, usedname)
	    dim%upper = max(dim%upper, dim%lower)
	    length = dim%upper - dim%lower + 1
	    dim%dim = Dimension(file, trim(dim%name), length)
	    call NetcdfAssert
	    dim%var = Variable(file, trim(dim%name), NF_FLOAT, (/dim%dim/))
	    call NetcdfAssert
	    call ReadAxisFile(dim)
	end subroutine

	subroutine StoreDimension(header, ia)
	    type(GT3_HEADER), intent(in):: header
	    integer, intent(in):: ia
	    character(len = 16):: item
	    integer:: idim
	    type(DIMTABLE_ENTRY), pointer:: dims_tmp(:)
	continue
	    ! 空なら飛ばす
	    item = vars_table(iv)%header%axis_item(ia)
	    if (item == '') return
	    ! 登録済みなら上下端だけ確認
	    do, idim = 1, num_dims
		if (item == dims_table(idim)%item) then
		    dims_table(idim)%lower = min(dims_table(idim)%lower, &
			& header%axis_start(ia))
		    dims_table(idim)%upper = max(dims_table(idim)%upper, &
			& header%axis_end(ia))
		endif
	    enddo
	    ! 未登録なら次元表を確保し書き込み
	    num_dims = num_dims + 1
	    if (num_dims > size(dims_table)) then
		allocate(dims_tmp(size(dims_table) + 2))
		dims_tmp(1: size(dims_table)) = dims_table(:)
		deallocate(dims_table)
		dims_table => dims_tmp
	    endif
	    dims_table(num_dims)%item = item
	    dims_table(num_dims)%lower = header%axis_start(ia)
	    dims_table(num_dims)%upper = header%axis_end(ia)
	end subroutine

    end subroutine BuildDimensions

    subroutine PutTimeData(file)
	type(NC_FILE), intent(in):: file
	integer:: it
	integer, pointer:: values(:)
    continue
	do, it = 1, size(time_table)
	    values => time_table(it)%values
	    if (.not. NetcdfPutInt(time_table(it)%var, &
		& values, size(values))) then
		call NetcdfAssert()
	    endif
	enddo
    end subroutine

    subroutine BuildTime(file)
	type(NC_FILE), intent(inout):: file
	character(len = 16):: axisname
	integer:: iv, num_times
    continue
	allocate(time_table(1))

	num_times = 0
	do, iv = 1, size(vars_table)
	    call StoreTime
	enddo
	time_table => time_table(1: num_times)

	call CreateTimeDimensions

    contains

	subroutine CreateTimeDimensions
	    integer:: it
	continue
	    ! 時間軸がないばあい
	    if (num_times <= 0) then
		time_table(:)%dim = NC_DIMENSION_ERROR()
		return
	    endif

	    ! 時間軸が一つしかないばあい
	    time_table(1)%dim = &
		& Dimension(file, 'time', size(time_table(1)%values))
	    time_table(1)%var = &
		& Variable(file, 'time', NF_INT, (/time_table(1)%dim/))
	    if (num_times == 1) return

	    ! 時間軸が複数あるばあい
	    do, it = 2, num_times
		axisname = NetcdfName('time', it)
		time_table(it)%dim = &
		    & Dimension(file, axisname, size(time_table(1)%values))
		time_table(it)%var = &
		    & Variable(file, axisname, NF_INT, (/time_table(it)%dim/))
	    enddo
	end subroutine

	subroutine StoreTime
	    integer:: it, len
	    integer, pointer:: t(:)
	    type(TIMETABLE_ENTRY), pointer:: tt(:)
	continue
	    len = count(units_table(:)%var == iv)
	    if (len < 1) stop 'gt3conv_v#StoreTime barf'
	    allocate(t(len))
	    t = pack(units_table(:)%time, units_table(:)%var == iv)

	    ! 既存時間軸の探索
	    do, it = 1, num_times
		if (int_array_equiv(t, time_table(it)%values)) then
		    vars_table(iv)%time_ord = it
		    deallocate(t)
		    return
		endif
	    enddo

	    ! 新しい時間軸の登録
	    num_times = num_times + 1
	    len = size(time_table)
	    if (num_times > len) then
		allocate(tt(len + 1))
		tt(1: len) = time_table(1: len)
		deallocate(time_table)
		time_table => tt
	    endif
	    vars_table(iv)%time_ord = num_times
	    time_table(num_times)%values => t
	end subroutine

    end subroutine

    logical function int_array_equiv(a, b) result(result)
	integer, intent(in):: a(:), b(:)
    continue
	if (size(a) /= size(b)) then
	    result = .FALSE.;  return
	endif
	result = all(a(:) == b(:))
    end function

end module