uEMEP_region_mask Subroutine

    subroutine uEMEP_region_mask()
        ! This routine defines use_subgrid_val=2 for regions outside the selected region. This can be used to control the interpolation routine
        ! It also sets use_subgrid=.false. outside the region so that no calculations are made there either
        character(256) :: temp_name, temp_str, temp_str1
        integer, allocatable :: tile_municipality_subgrid(:,:,:)
        integer :: municipality_id
        real :: x_ssb, f_easting, ssb_dx, y_ssb, ssb_dy
        integer :: count
        integer(kind=8) :: ssb_id
        integer :: i, j, k, i_range, j_range, i_range_interp, j_range_interp, i_tile, j_tile
        logical :: exists
        integer :: unit_in
        integer :: index_val
        integer :: i_source
        character(256) :: region_number_str
        integer, parameter :: n_search = 5
        character(16) :: search_str(n_search)
        real :: search_delta(n_search)
        integer :: temp_search
        integer :: io

        search_str = ['_1000m', '_500m ', '_250m ', '_100m ', '_50m  ']
        search_delta = [1000.0, 500.0, 250.0, 100.0, 50.0]

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '================================================================'
        write(unit_logfile,'(A)') 'Masking region (uEMEP_region_mask)'
        write(unit_logfile,'(A)') '================================================================'

        if ( .not. allocated(tile_municipality_subgrid)) allocate (tile_municipality_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),2))
        tile_municipality_subgrid = 0

        i_source = allsource_index

        ! Search file name to define the grid size
        ssb_dx = 0.0; ssb_dy = 0.0
        do k = 1, n_search
            temp_search = index(filename_population(municipality_index), trim(adjustl(search_str(k))))
            if (temp_search .ne. 0) then
                ssb_dx = search_delta(k)
                ssb_dy = search_delta(k)
                write(unit_logfile,'(i,A)') temp_search, ' Reading municipality data with resolution '//trim(adjustl(search_str(k)))
            end if
        end do

        if (ssb_dx .eq. 0) then
            write(unit_logfile,'(A)') 'Cannot find a valid SSB grid size. Stopping. '//trim(filename_population(municipality_index))
            stop 1
        end if

        region_number_str = ''
        write(region_number_str,*) region_index
        region_number_str = trim(region_number_str)//'_'

        ! Read in SSB file containing gridded municipality ids
        f_easting=2.0e6
        pathfilename_population(municipality_index) = trim(pathname_population(municipality_index))//trim(adjustl(region_number_str))//trim(filename_population(municipality_index))
        
        ! Test existence of the heating filename. If does not exist then use default
        inquire(file=trim(pathfilename_population(municipality_index)), exist=exists)
        if ( .not. exists) then
            write(unit_logfile,'(A,A)') ' ERROR: Masking region SSB file with municipality IDs does not exist: ', trim(pathfilename_population(municipality_index))
            stop 1
        end if
        temp_name = pathfilename_population(municipality_index)
        
        ! Open the file for reading
        unit_in = 20
        open(unit_in, file=temp_name, access='sequential', status='old', readonly)
        write(unit_logfile,'(a)') ' Opening SSB municipality file '//trim(temp_name)
        rewind(unit_in)
        
        ! Read header SSBID0250M;kommunenum
        read(unit_in,'(A)') temp_str
        write(unit_logfile,'(A)') 'Header: '//trim(temp_str)
        count = 0
        i_range = ceiling(ssb_dx/subgrid_delta(x_dim_index)/2.0) + 1
        j_range = ceiling(ssb_dy/subgrid_delta(y_dim_index)/2.0) + 1
        i_range_interp = ceiling(ssb_dx/subgrid_delta(x_dim_index)/2.0 + max_interpolation_subgrid_size/subgrid_delta(x_dim_index)) + 1
        j_range_interp = ceiling(ssb_dy/subgrid_delta(y_dim_index)/2.0 + max_interpolation_subgrid_size/subgrid_delta(y_dim_index)) + 1
        do
            ssb_id = 0; municipality_id = 0
            ! Read in file string
            read(unit_in,'(A)',iostat=io) temp_str
            if (io /= 0) exit
            index_val = index(temp_str,';', back=.false.)
            temp_str1 = temp_str(1:index_val-1)
            temp_str = temp_str(index_val+1:)
            if (index_val .gt. 1) read(temp_str1,*) ssb_id
            read(temp_str,*) municipality_id
            ! If this ssb municipality grid has the correct ID then find the target grid that matches it
            if (municipality_id .eq. region_id) then
                count = count + 1

                !Convert id to grid centre coordinates that are already in UTM33 for SSB data
                x_ssb = floor(ssb_id/10000000.0) - f_easting + ssb_dx/2.0
                y_ssb = mod(ssb_id, 10000000) + ssb_dy/2.0

                ! Find the tile this ssb grid is in
                i_tile = 1 + floor((x_ssb - subgrid_min(x_dim_index))/subgrid_delta(x_dim_index)) ! New definition
                j_tile = 1 + floor((y_ssb - subgrid_min(y_dim_index))/subgrid_delta(y_dim_index)) !New definition

                do j = j_tile-j_range, j_tile+j_range
                    do i = i_tile-i_range, i_tile+i_range
                        ! Make sure i and j are inside a valid range
                        if (i .ge. 1 .and. i .le. subgrid_dim(x_dim_index) .and. j .ge. 1 .and. j .le. subgrid_dim(y_dim_index)) then
                            ! Find the target subgrid within the ssb grid and the region
                            if (x_subgrid(i,j) .ge. x_ssb-ssb_dx/2.0 .and. x_subgrid(i,j) .lt. x_ssb + ssb_dx/2.0 .and. &
                                y_subgrid(i,j) .ge. y_ssb-ssb_dy/2.0 .and. y_subgrid(i,j) .lt. y_ssb+ssb_dy/2.0) then
                                !If there is a target subgrid within this range then allocate the mask value if it is not the correct region_id
                                tile_municipality_subgrid(i,j,1) = 1
                            end if
                        end if
                    end do
                end do

                do j = j_tile-j_range_interp, j_tile+j_range_interp
                    do i = i_tile-i_range_interp, i_tile+i_range_interp
                        ! Make sure i and j are inside a valid range
                        if (i .ge. 1 .and. i .le. subgrid_dim(x_dim_index) .and. j .ge. 1 .and. j .le. subgrid_dim(y_dim_index)) then
                            ! Find the target subgrid within the ssb grid and the region
                            if (x_subgrid(i,j) .ge. x_ssb-ssb_dx/2.0 - max_interpolation_subgrid_size .and. x_subgrid(i,j) .le. x_ssb + ssb_dx/2.0 + max_interpolation_subgrid_size .and. &
                                y_subgrid(i,j) .ge. y_ssb-ssb_dy/2.0 - max_interpolation_subgrid_size .and. y_subgrid(i,j) .le. y_ssb+ssb_dy/2.0 + max_interpolation_subgrid_size) then
                                ! If there is a target subgrid within this range then allocate the mask value if it is not the correct region_id
                                tile_municipality_subgrid(i,j,2) = 1
                            end if
                        end if
                    end do
                end do
            end if
        end do
        close(unit_in)

        do j = 1, subgrid_dim(y_dim_index)
            do i = 1, subgrid_dim(x_dim_index)
                if (tile_municipality_subgrid(i,j,2) .eq. 0) then
                    use_subgrid_val(i,j,:) = outside_interpolation_region_index
                    use_subgrid(i,j,:) = .false.
                end if
                if (tile_municipality_subgrid(i,j,1) .eq. 0) then
                    use_subgrid_val(i,j,:) = outside_region_index
                end if
            end do
        end do

        if (allocated(tile_municipality_subgrid)) deallocate (tile_municipality_subgrid)

    end subroutine uEMEP_region_mask