! Created by Andrey Debolskiy on 15.08.2024.

module external_forcing
    implicit none
    !type declarations
    type forcing_1d
        real, allocatable :: x(:) ! 0 - closed, 1 - opened
        real, allocatable :: val(:)
        integer :: xsize
    end type forcing_1d
    type forcing_2d
        real, allocatable :: x(:) ! 0 - closed, 1 - opened
        real, allocatable :: y(:) ! 0 - closed, 1 - opened
        real, allocatable :: val(:,:)
        integer :: xsize, ysize
    end type forcing_2d

    public :: ext_forc_dealloc_1d, ext_forc_dealloc_2d
    public :: get_1d_interpolated, get_1d_to_1d_interpolated, &
                read_2d_focing


    contains
        real function get_1d_interpolated(xq, forcing1d) result(val)
            real, intent(in) :: xq
            type (forcing_1d), intent(in) :: forcing1d
            !local
            integer i
            real x1,x2, val1, val2

            if (xq < forcing1d%x(1)) then
                val = forcing1d%val(1)
                return
            end if

            do i = 1, forcing1d%xsize
                if (forcing1d%x(i) >= xq) exit
            end do

            if (i == forcing1d%xsize) then
                val = forcing1d%val(i)
                return
            else
                x1 = forcing1d%x(i-1)
                x2 = forcing1d%x(i)
                val1  = forcing1d%val(i-1)
                val2  = forcing1d%val(i)
                val = val1 + (val2 - val1) * (xq - x1) / (x2 -x1)
                return
            end if
        end function get_1d_interpolated

        subroutine get_1d_to_1d_interpolated(valq, xq, nq, forcing1d)
            integer, intent(in):: nq
            real, dimension(nq), intent(in) :: xq
            real, dimension(nq), intent(out) :: valq
            type (forcing_1d), intent(in) :: forcing1d
            !local
            integer k
            do k = 1, nq
                valq(k) = get_1d_interpolated(xq(k), forcing1d)
            enddo
        end subroutine get_1d_to_1d_interpolated

        subroutine read_2d_focing(out_forc, filename, filename_x, filename_y)
            implicit none
            type (forcing_2d), intent(out):: out_forc
            character(*), intent(in):: filename, filename_x, filename_y

            !local vars
            integer nx,ny, i, j
            real x,y
            open(10, FILE=trim(filename))
            read(10,*) nx
            read(10,*) ny

            out_forc%xsize=nx
            out_forc%ysize=ny
            allocate(out_forc%val(nx,ny))
            do i  = 1, nx
                READ(10,*) (out_forc%val(i,j),j=1,ny)
            end do
            close(10)
            ! read coords
            allocate(out_forc%x(nx))
            open(10, FILE=trim(filename_x))
            do i=1,nx
                read(10,*) out_forc%x(i)
            end do
            close(10)

            allocate(out_forc%y(ny))
            open(10, FILE=trim(filename_y))
            do i=1,ny
                read(10,*) out_forc%y(i)
            end do
            close(10)
        end subroutine read_2d_focing

        subroutine ext_forc_dealloc_1d(forc_1d)
            implicit none
            type (forcing_1d), intent(inout) :: forc_1d
                deallocate(forc_1d%val)
                deallocate(forc_1d%x)
        end subroutine ext_forc_dealloc_1d

        subroutine ext_forc_dealloc_2d(forc_2d)
            implicit none
            type (forcing_2d), intent(inout) :: forc_2d
            deallocate(forc_2d%val)
            deallocate(forc_2d%x)
            deallocate(forc_2d%y)
        end subroutine ext_forc_dealloc_2d
end module external_forcing