read_esri_ascii_3d_file Subroutine

    subroutine read_esri_ascii_3d_file(filename_ascii_sub, ncols_sub, nrows_sub, nblocks_sub, cellsize_sub, val_array, x_array, y_array)
        character(*) :: filename_ascii_sub
        integer :: ncols_sub, nrows_sub, nblocks_sub
        real :: cellsize_sub
        real :: val_array(ncols_sub, nrows_sub, nblocks_sub)
        real :: x_array(ncols_sub, nrows_sub)
        real :: y_array(ncols_sub, nrows_sub)

        ! Local variables
        character(256) :: temp_str
        integer :: i, j, ii, jj, tt
        integer :: ncols_sub_temp, nrows_sub_temp, nblocks_sub_temp
        real :: nrows_sub_temp_real
        real :: xllcorner
        real :: yllcorner
        real :: NODATA_value
        integer :: unit_in = 20

        open(unit_in,file=filename_ascii_sub, access='sequential', form='formatted', status='old', readonly)
        rewind(unit_in)
        read(unit_in,*) temp_str, ncols_sub_temp
        read(unit_in,*) temp_str, nrows_sub_temp
        read(unit_in,*) temp_str, nrows_sub_temp_real
        
        if (temp_str .eq. 'xllcorner') then
            nblocks_sub_temp = 1
            xllcorner = nrows_sub_temp_real
        else
            nrows_sub_temp = int(nrows_sub_temp_real)
            read(unit_in,*) temp_str, xllcorner
        end if

        read(unit_in,*) temp_str, yllcorner
        read(unit_in,*) temp_str, cellsize_sub
        read(unit_in,*) temp_str, NODATA_value

        if (ncols_sub .ne. ncols_sub_temp .or. nrows_sub .ne. nrows_sub_temp .or. nblocks_sub .ne. nblocks_sub_temp) then
            write(unit_logfile,'(A,2I,A,2I,A,2I)') 'ERROR: Subgrid columns or rows do not match. Columns: ', ncols_sub, ncols_sub_temp, '  Rows: ', nrows_sub, nrows_sub_temp, '  Blocks: ', nblocks_sub, nblocks_sub_temp
            stop 1
        end if

        read(unit_in,*) (((val_array(ii,jj,tt), ii = 1, ncols_sub), jj = nrows_sub, 1, -1), tt = 1, nblocks_sub)
        close(unit_in)

        ! Set position arrays
        do i = 1, ncols_sub
            do j = 1, nrows_sub
                x_array(i,j) = xllcorner + cellsize_sub/2.0 + (i-1)*cellsize_sub
                y_array(i,j) = yllcorner + cellsize_sub/2.0 + (j-1)*cellsize_sub
            end do
        end do
    end subroutine read_esri_ascii_3d_file