uEMEP_interpolate_auto_subgrid Subroutine

    subroutine uEMEP_interpolate_auto_subgrid()
        ! This is the corresponding routine for interpolating the auto selected data
        integer :: xdim, ydim
        real :: delta(2)
        real :: xval, yval
        real :: xgrid(3,3), ygrid(3,3), zgrid(3,3)

        integer :: i, j, k, ii, jj
        integer :: i_source, i_pollutant, t
        integer :: i_in(3), j_in(3)
        integer :: use_subgrid_step

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '================================================================'
        write(unit_logfile,'(A)') 'Interpolating auto grid (uEMEP_interpolate_auto_subgrid)'
        write(unit_logfile,'(A)') '================================================================'

        xdim = 3
        ydim = 3

        do i_source = 1, n_source_index
            if (calculate_source(i_source) .and. i_source .ne. allsource_index .and. use_emission_positions_for_auto_subgrid_flag(i_source)) then
                ! Only interpolate for chosen sources that have been auto gridded based on emissions
                do k = n_use_subgrid_levels(i_source), 1, -1
                    use_subgrid_step = use_subgrid_step_delta(k)
                    delta = use_subgrid_step*subgrid_delta
                    do j = 1, subgrid_dim(y_dim_index)
                        do i = 1, subgrid_dim(x_dim_index)
                            xval = x_subgrid(i,j)
                            yval = y_subgrid(i,j)
                            ! Only do the interpolation if it is the right interpolation_index, it is inside the region and it is not a valid subgrid
                            ! TODO: use a logical flag here?
                            ! if (use_subgrid_interpolation_index(i,j,i_source).eq.k.and.use_subgrid_val(i,j,i_source).ne.outside_interpolation_region_index.and..not.use_subgrid_val(i,j,i_source)) then
                            ! Do it everywhere in the grid
                            if (use_subgrid_interpolation_index(i,j,i_source) .eq. k .and. use_subgrid_val(i,j,i_source) .ne. &
                                outside_interpolation_region_index .and. .not. use_subgrid(i,j,i_source)) then

                                i_in(2) = floor(real(i - 1)/use_subgrid_step + 0.5)*use_subgrid_step + 1
                                i_in(2) = min(subgrid_dim(x_dim_index), max(1, i_in(2)))
                                i_in(1) = i_in(2) - use_subgrid_step
                                i_in(3) = i_in(2) + use_subgrid_step

                                j_in(2) = floor(real(j - 1)/use_subgrid_step + 0.5)*use_subgrid_step + 1
                                j_in(2) = min(subgrid_dim(y_dim_index), max(1, j_in(2)))
                                j_in(1) = j_in(2) - use_subgrid_step
                                j_in(3) = j_in(2) + use_subgrid_step

                                if (i_in(1) .lt. 1) i_in(1) = i_in(2)
                                if (j_in(1) .lt. 1) j_in(1) = j_in(2)
                                if (i_in(3) .gt. subgrid_dim(x_dim_index)) i_in(3) = i_in(2)
                                if (j_in(3) .gt. subgrid_dim(y_dim_index)) j_in(3) = j_in(2)

                                do t = 1, subgrid_dim(t_dim_index)
                                    do i_pollutant = 1, n_pollutant_loop
                                        do jj = 1, 3
                                            do ii = 1, 3
                                                xgrid(ii,jj) = x_subgrid(i_in(2),j_in(2)) + (ii - 2)*delta(1)
                                                ygrid(ii,jj) = y_subgrid(i_in(2),j_in(2)) + (jj - 2)*delta(2)
                                                zgrid(ii,jj) = subgrid(i_in(ii),j_in(jj),t,proxy_subgrid_index,i_source,i_pollutant)
                                                if (i_in(ii) .lt. 1 .or. i_in(ii) .gt. subgrid_dim(x_dim_index) .or. j_in(jj) .lt. 1 .or. j_in(jj) .gt. subgrid_dim(y_dim_index)) then
                                                    zgrid(ii,jj) = subgrid(i_in(2),j_in(2),t,proxy_subgrid_index,i_source,i_pollutant)
                                                else if (use_subgrid_val(i_in(ii),j_in(jj),i_source).eq.outside_interpolation_region_index) then
                                                    zgrid(ii,jj) = subgrid(i_in(2),j_in(2),t,proxy_subgrid_index,i_source,i_pollutant)
                                                end if
                                            end do
                                        end do
                                        subgrid(i,j,t,proxy_subgrid_index,i_source,i_pollutant)=area_weighted_interpolation_function(xgrid,ygrid,zgrid,xdim,ydim,delta,xval,yval)

                                        ! Travel time interpolation as well
                                        do jj = 1, 3
                                            do ii = 1, 3
                                                zgrid(ii,jj) = traveltime_subgrid(i_in(ii),j_in(jj),t,1,i_pollutant)
                                                if (i_in(ii) .lt. 1 .or. i_in(ii) .gt. subgrid_dim(x_dim_index) .or. j_in(jj) .lt. 1 .or. j_in(jj) .gt. subgrid_dim(y_dim_index)) then
                                                    zgrid(ii,jj) = traveltime_subgrid(i_in(2),j_in(2),t,1,i_pollutant)
                                                else if (use_subgrid_val(i_in(ii),j_in(jj),i_source) .eq. outside_interpolation_region_index) then
                                                    zgrid(ii,jj) = traveltime_subgrid(i_in(2),j_in(2),t,1,i_pollutant)
                                                end if
                                            end do
                                        end do
                                        traveltime_subgrid(i,j,t,1,i_pollutant) = area_weighted_interpolation_function(xgrid, ygrid, zgrid, xdim, ydim, delta, xval, yval)

                                        do jj = 1, 3
                                            do ii = 1, 3
                                                zgrid(ii,jj) = traveltime_subgrid(i_in(ii),j_in(jj),t,2,i_pollutant)
                                                if (i_in(ii) .lt. 1 .or. i_in(ii) .gt. subgrid_dim(x_dim_index) .or. j_in(jj) .lt. 1 .or. j_in(jj) .gt. subgrid_dim(y_dim_index)) then
                                                    zgrid(ii,jj) = traveltime_subgrid(i_in(2),j_in(2),t,2,i_pollutant)
                                                else if (use_subgrid_val(i_in(ii),j_in(jj),i_source) .eq. outside_interpolation_region_index) then
                                                    zgrid(ii,jj) = traveltime_subgrid(i_in(2),j_in(2),t,2,i_pollutant)
                                                end if
                                            end do
                                        end do
                                        traveltime_subgrid(i,j,t,2,i_pollutant) = area_weighted_interpolation_function(xgrid, ygrid, zgrid, xdim, ydim, delta, xval, yval)
                                    end do
                                end do
                            end if
                        end do
                    end do
                end do
            end if
        end do

        ! Reset the use_subgrid values so chemistry and exposure happens everywhere but not outside the region
        use_subgrid(:,:,allsource_index) = .true.
        where (use_subgrid_val(:,:,allsource_index) .eq. outside_region_index) use_subgrid(:,:,allsource_index) = .false.

    end subroutine uEMEP_interpolate_auto_subgrid