uEMEP_subgrid_EMEP_from_in_region Subroutine

    subroutine uEMEP_subgrid_EMEP_from_in_region

        ! contribution to centre of each EMEP cell from in-region, within MW and outside MW, for each region
        real, allocatable :: EMEP_local_from_in_region(:, :, :, :, :, :) ! (x,y,region,t,source,pollutant)
        real, allocatable :: EMEP_semilocal_from_in_region(:, :, :, :, :, :) ! (x,y,region,t,source,pollutant)
        ! arrays for accumulating data for a single EMEP cell of the above arrays (is this actually faster than accumulating directly in the big arrays????)
        real, allocatable :: temp_EMEP_local_from_in_region(:, :, :, :)
        real, allocatable :: temp_EMEP_semilocal_from_in_region(:, :, :, :)
        ! additional increment for the centre of each EMEP cell, for each region
        real, allocatable :: EMEP_additional_increment_from_in_region(:, :, :, :, :, :)  ! (x,y,region,t,source,pollutant)
        ! additional increment for a single EMEP cell, for each region, to be accumulated
        real, allocatable :: temp_EMEP_additional_increment_from_in_region(:, :, :, :)  ! (region,t,source,pollutant)
        ! additional increment of a single big LF grid to a single receptor EMEP grid (before weighing by region)
        real, allocatable :: EMEP_additional_increment_current_lfgrid(:, :, :)  ! (t,source,pollutant)
        ! weighting of the additional increment, for each region
        real, allocatable :: weights_EMEP_additional_increment_current_lfgrid(:)  ! (region)
        ! weighting of an LF cell for local and semilocal contribution
        real weighting_value_local,weighting_value_semilocal
        ! indexers for looping
        integer i,j             ! target subgrids
        integer ii,jj           ! EMEP grids
        integer i_dist,j_dist   ! LF dimensions
        integer iiii,jjjj       ! small LF grids within a big LF grid
        integer i_sub,j_sub     ! subsamples of an EMEP grid
        integer i_region        ! region dimension
        integer i_source        ! source dimension
        ! indexers for determining positioning of additional LF grids
        integer ii_start,jj_start
        ! additinoal indexers
        integer iii,jjj
        integer iii_nc,jjj_nc
        ! displacement distances in LF grids (whole grids)
        integer x_dist,y_dist
        integer xdist_big,ydist_big
        integer xdist_small,ydist_small
        integer xdist_small_first,ydist_small_first
        integer idist_small,jdist_small
        integer max_x_dist,max_y_dist
        ! indexers for the location of a 1x1 LF grid cell in the EMEP grid
        integer iiii_nc,jjjj_nc
        ! indexers for the location of a 1x1 LF grid cell in the extended EMEP grid
        integer iiii_extended,jjjj_extended
        ! indexers for finding local contributions in LF array
        integer lc_index,lc_additional_index
        ! counters
        integer counter,counter_local,counter_semilocal
        ! variables to hold a region index
        integer current_region_index
        ! fractional position of uEMEP subgrid within an EMEP grid
        real ii_frac_target,jj_frac_target
        ! location of subgrid in EMEP's coordinate system
        real x_temp,y_temp
        ! distance between target subgrid and an EMEP grid centre or EMEP subsample location
        real x_dist_target,y_dist_target
        ! half-size of the moving window
        real n_EMEP_grids_to_edge_of_moving_window
        ! fraction of an EMEP grid that is in the correct region
        real current_EMEP_region_fraction
        ! weighting for area interpolation
        real weighting_val,weight_check

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '====================================================================================================='
        write(unit_logfile,'(A)') 'Calculation local and semilocal from-in-region EMEP contributions to subgrids (uEMEP_subgrid_EMEP_from_in_region)'
        write(unit_logfile,'(A)') '====================================================================================================='

        write(unit_logfile,'(A,8I8)') 'dims: (x_emep,y_emep,reg,time,source,pollutant,xdist,ydist) = ',dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop,dim_length_nc(xdist_dim_nc_index),dim_length_nc(ydist_dim_nc_index)

        ! ************************************************
        ! PART 1: Calculate contributions to the EMEP grid
        ! ************************************************

        ! Indices in lc_var3d_nc to find the normal and additional local contributions
        lc_index = lc_local_nc_loop_index(local_fraction_grid_for_EMEP_grid_interpolation)
        lc_additional_index = lc_local_nc_loop_index(local_fraction_grid_for_EMEP_additional_grid_interpolation)

        ! distance in x or y (EMEP grid) from receptor subgrid to edge of moving window, in units of EMEP grids
        ! (normally whole number, but not if EMEP_grid_interpolation_size is odd number)
        n_EMEP_grids_to_edge_of_moving_window = EMEP_grid_interpolation_size*0.5
        write(unit_logfile,'(A,f8.1)') 'n_EMEP_grids_to_edge_of_moving_window=', n_EMEP_grids_to_edge_of_moving_window

        ! Calculate local and semilocal contributions to centre of all EMEP grids
        write(unit_logfile,'(A)') 'Allocating arrays for calculating in-region local and semilocal contribution to the EMEP grid'

        ! Allocate arrays
        allocate(EMEP_local_from_in_region(dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
        EMEP_local_from_in_region = 0.0
        allocate(EMEP_semilocal_from_in_region(dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
        EMEP_semilocal_from_in_region = 0.0
        allocate(temp_EMEP_local_from_in_region(n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
        allocate(temp_EMEP_semilocal_from_in_region(n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))

        ! Loop over all EMEP cells and calculate local and semilocal contributions when moving window is centered at the centre of the cell
        do ii = 1, dim_length_nc(x_dim_nc_index)
            do jj = 1, dim_length_nc(y_dim_nc_index)
                temp_EMEP_local_from_in_region = 0.0
                temp_EMEP_semilocal_from_in_region = 0.0
                ! Loop over all (small) LF contribution cells to this EMEP cell
                do i_dist = 1, dim_length_nc(xdist_dim_nc_index)
                    do j_dist = 1, dim_length_nc(ydist_dim_nc_index)
                        ! number of grids displaced relative to the receptor grid cell
                        x_dist = i_dist - xdist_centre_nc
                        y_dist = j_dist - ydist_centre_nc
                        ! index in the netcdf file of this cell
                        iiii_nc = ii + x_dist
                        jjjj_nc = jj + y_dist
                        ! index in the extended EMEP grid of this cell
                        iiii_extended = iiii_nc + ngrid_extended_margin
                        jjjj_extended = jjjj_nc + ngrid_extended_margin

                        ! Check if this cell is covered by the extended grid. If not, assume it is completely outside all regions
                        if (.not. (iiii_extended >= 1 .and. iiii_extended <= nx_EMEP_extended .and. jjjj_extended >= 1 .and. jjjj_extended <= ny_EMEP_extended)) then
                            cycle
                        end if

                        ! Loop over all regions
                        do i_region = 1, n_regions

                            ! get fraction of the LF grid that is in the region
                            current_region_index = regionindex_loop_index(i_region)
                            current_EMEP_region_fraction = regionfraction_per_EMEP_extended_grid(iiii_extended, jjjj_extended, i_region)

                            ! Determine the weights to use for this LF cell
                            ! for local, the weight is the area fraction that is in the region and inside the moving window
                            ! for semilocal, the weight is the area fraction that is in the region but outside the moving window
                            if (current_EMEP_region_fraction <= 0) then
                                ! No part of this EMEP grid is within the region
                                cycle
                            else if (abs(x_dist) <= n_EMEP_grids_to_edge_of_moving_window-1 .and. abs(y_dist) <= n_EMEP_grids_to_edge_of_moving_window-1) then
                                ! this LF grid is sure to be completely within the moving window
                                weighting_value_local = current_EMEP_region_fraction
                                weighting_value_semilocal = 0.0
                            else if (abs(x_dist) >= n_EMEP_grids_to_edge_of_moving_window+1 .or. abs(y_dist) >= n_EMEP_grids_to_edge_of_moving_window+1) then
                                ! this LF grid is sure to be completely outside the moving window
                                weighting_value_local = 0.0
                                weighting_value_semilocal = current_EMEP_region_fraction
                            else
                                ! this LF grid might be partly covered by the moving window
                                ! -> we must go through all subsamples of that grid to find the overlap between region and moving window
                                counter_local = 0       ! count the subsample grids inside-moving-window & in-region
                                counter_semilocal = 0   ! count the subsample grids outside-moving-window & in-region
                                do i_sub = 1, n_subsamples_per_EMEP_grid
                                    do j_sub = 1, n_subsamples_per_EMEP_grid
                                        ! first check if this subsample is in the region
                                        if (EMEP_extended_subsample_region_id(i_sub,j_sub,iiii_extended,jjjj_extended) == current_region_index) then
                                            ! deduce x- and y-distance (in number of EMEP grids) from this subsample location to the midpoint of the receptor grid-cell
                                            x_dist_target = x_dist + (i_sub-0.5)/n_subsamples_per_EMEP_grid - 0.5
                                            y_dist_target = y_dist + (j_sub-0.5)/n_subsamples_per_EMEP_grid - 0.5
                                            ! use these distances to determine whether the subsample is inside the moving window
                                            if (abs(x_dist_target) <= n_EMEP_grids_to_edge_of_moving_window .and. abs(y_dist_target) <= n_EMEP_grids_to_edge_of_moving_window) then
                                                ! this subsample is inside the moving window
                                                counter_local = counter_local + 1
                                            else
                                                ! this subsample is outside the moving window
                                                counter_semilocal = counter_semilocal + 1
                                            end if
                                        end if
                                    end do
                                end do
                                ! Divide counters by total number of subsamples per EMEP grid to get the area fraction
                                weighting_value_local = counter_local*1.0/n_subsamples_per_EMEP_grid**2
                                weighting_value_semilocal = counter_semilocal*1.0/n_subsamples_per_EMEP_grid**2
                            end if
                            temp_EMEP_local_from_in_region(i_region,:,:,:) = temp_EMEP_local_from_in_region(i_region,:,:,:) + lc_var3d_nc(i_dist,j_dist,ii,jj,:,lc_index,1:n_source_index,:)*weighting_value_local
                            temp_EMEP_semilocal_from_in_region(i_region,:,:,:) = temp_EMEP_semilocal_from_in_region(i_region,:,:,:) + lc_var3d_nc(i_dist,j_dist,ii,jj,:,lc_index,1:n_source_index,:)*weighting_value_semilocal
                        end do ! i_region = 1, n_regions

                    end do
                end do
                EMEP_local_from_in_region(ii,jj,:,:,:,:) = temp_EMEP_local_from_in_region
                EMEP_semilocal_from_in_region(ii,jj,:,:,:,:) = temp_EMEP_semilocal_from_in_region
            end do
        end do

        ! Deallocate temporary arrays
        deallocate(temp_EMEP_local_from_in_region)
        deallocate(temp_EMEP_semilocal_from_in_region)

        ! Calculate the additional incremental contribution to each EMEP grid for each region
        ! NB: we let the interpolation domain extend to the very edge of the LF data, not limited by the EMEP_additional_grid_interpolation_size (this is only used to check >0 to determine whether we should calculate it or not)
        if (EMEP_additional_grid_interpolation_size > 0.0) then

            write(unit_logfile,'(A)') 'Allocating arrays for additional increment calculation'

            allocate(EMEP_additional_increment_from_in_region(dim_length_nc(x_dim_nc_index),dim_length_nc(y_dim_nc_index),n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
            allocate(temp_EMEP_additional_increment_from_in_region(n_regions,subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
            allocate(EMEP_additional_increment_current_lfgrid(subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
            allocate(weights_EMEP_additional_increment_current_lfgrid(n_regions))

            write(unit_logfile,'(A)') 'Calculating additional increment to all EMEP grids'

            ! Use the starting position of the read in EMEP file to initialise the starting point
            ii_start = mod(dim_start_EMEP_nc(x_dim_nc_index)-1,local_fraction_grid_size(2))
            jj_start = mod(dim_start_EMEP_nc(y_dim_nc_index)-1,local_fraction_grid_size(2))

            ! Deduce the max distance (+/-) we have LF data for
            max_x_dist = (dim_length_nc(xdist_dim_nc_index)-1)/2
            max_y_dist = (dim_length_nc(ydist_dim_nc_index)-1)/2

            ! Loop over all EMEP grids
            EMEP_additional_increment_from_in_region = 0.0
            do ii = 1, dim_length_nc(x_dim_nc_index)
                do jj = 1, dim_length_nc(y_dim_nc_index)
                    ! Initialize the additional increment of this EMEP cell to zero
                    temp_EMEP_additional_increment_from_in_region = 0.0
                    ! EMEP grid index of bottom-left-corner-cell of the additional grid associated with that EMEP grid
                    iii = int((ii-1+ii_start)/local_fraction_grid_size(2))*local_fraction_grid_size(2) + 1 - ii_start
                    jjj = int((jj-1+jj_start)/local_fraction_grid_size(2))*local_fraction_grid_size(2) + 1 - jj_start

                    ! Loop over all big (additional) LF grids that give contributions to the EMEP grid
                    do i_dist = 1, dim_length_nc(xdist_dim_nc_index)
                        do j_dist = 1, dim_length_nc(ydist_dim_nc_index)

                            ! Initialize additional increment to the total additional contribution from EMEP from this cell (all times, sources and pollutants)
                            EMEP_additional_increment_current_lfgrid = lc_var3d_nc(i_dist,j_dist,ii,jj,:,lc_additional_index,1:n_source_index,:)
                            
                            ! the x_dist and y_dist of this big LF grid (..., -1, 0, 1, ...)
                            xdist_big = i_dist - xdist_centre_nc
                            ydist_big = j_dist - ydist_centre_nc
                            ! deduce what is the xdist and ydist in the 1x1 LF grid of the lower-left EMEP cell falling within this big LF grid
                            xdist_small_first = iii - ii + xdist_big*local_fraction_grid_size(2)
                            ydist_small_first = jjj - jj + ydist_big*local_fraction_grid_size(2)
                            ! loop over all EMEP grids contained within this big LF grid
                            ! Reset weights and counter
                            weights_EMEP_additional_increment_current_lfgrid = 0.0
                            counter = 0  ! count the grids not covered by the small LF grid
                            do iiii = 1, local_fraction_grid_size(2)
                                do jjjj = 1, local_fraction_grid_size(2)
                                    ! Deduce the xdist and ydist of this EMEP grid in the small LF domain
                                    xdist_small = xdist_small_first - 1 + iiii
                                    ydist_small = ydist_small_first - 1 + jjjj
                                    ! and corresponding index in the LF array
                                    idist_small = xdist_small + xdist_centre_nc
                                    jdist_small = ydist_small + ydist_centre_nc
                                    if (abs(xdist_small) > max_x_dist .or. abs(ydist_small) > max_y_dist) then
                                        ! This grid is NOT covered by 1x1 LF data, so add its region coverage to the weight
                                        ! find index in the normal EMEP grid
                                        iiii_nc = ii + xdist_small
                                        jjjj_nc = jj + ydist_small
                                        ! find corresponding index in the extended EMEP grid of region fractions
                                        iiii_extended = iiii_nc + ngrid_extended_margin
                                        jjjj_extended = jjjj_nc + ngrid_extended_margin
                                        ! check if this is within the extended region mask grid
                                        if (iiii_extended >= 1 .and. iiii_extended <= nx_EMEP_extended .and. jjjj_extended >= 1 .and. jjjj_extended <= ny_EMEP_extended) then
                                            weights_EMEP_additional_increment_current_lfgrid = weights_EMEP_additional_increment_current_lfgrid + regionfraction_per_EMEP_extended_grid(iiii_extended, jjjj_extended, :)
                                        end if
                                        ! NB: If the cell is is outside the extended EMEP grid, it is considered to be outside all the regions. If the extended grid size is set properly, this should never be the case for the receptor EMEP grids (ii,jj) that overlap with the target grid.
                                        counter = counter + 1
                                    else
                                        ! The grid is covered by 1x1 LF data: subtract the 1x1 LF from the additional increment
                                        EMEP_additional_increment_current_lfgrid=EMEP_additional_increment_current_lfgrid-lc_var3d_nc(idist_small,jdist_small,ii,jj,:,lc_index,1:n_source_index,:)
                                    end if
                                end do
                            end do
                            ! ensure the additional increment is not smaller than zero
                            where (EMEP_additional_increment_current_lfgrid < 0.0) EMEP_additional_increment_current_lfgrid = 0.0
                            ! normalize the weights by the number of grids
                            ! NB: if counter = 0, then the weights are zero so we can go to next LF source grid
                            if (counter > 0) then
                                ! normalize weights by the number of cells summed over
                                weights_EMEP_additional_increment_current_lfgrid = weights_EMEP_additional_increment_current_lfgrid/counter
                                ! For each region, multiply the additional increment by the weight calculated for that region
                                ! and accumulate this in the array for the total additional increment (to be accumulated over all the big LF cells)
                                do i_region = 1, n_regions
                                    temp_EMEP_additional_increment_from_in_region(i_region,:,:,:)=temp_EMEP_additional_increment_from_in_region(i_region,:,:,:)+EMEP_additional_increment_current_lfgrid*weights_EMEP_additional_increment_current_lfgrid(i_region)
                                end do
                            end if
                        end do
                    end do
                    ! The additional increment has now been accumulated over all xdist and ydist source grids and weighted for each region
                    ! So now it can be inserted into the main array
                    EMEP_additional_increment_from_in_region(ii,jj,:,:,:,:) = temp_EMEP_additional_increment_from_in_region
                end do
            end do

            EMEP_semilocal_from_in_region = EMEP_semilocal_from_in_region + EMEP_additional_increment_from_in_region

            deallocate(temp_EMEP_additional_increment_from_in_region)
            deallocate(EMEP_additional_increment_current_lfgrid)
            deallocate(weights_EMEP_additional_increment_current_lfgrid)
            deallocate(EMEP_additional_increment_from_in_region)

        end if

        ! *****************************************
        ! PART 2: Interpolate to the target subgrid
        ! *****************************************

        write(unit_logfile,'(A)') 'Allocating arrays for in-region local and semilocal contribution to the target grid'

        ! Allocate the arrays for holding the results and initialize them to zero
        if (allocated(subgrid_EMEP_local_from_in_region)) deallocate(subgrid_EMEP_local_from_in_region)
        allocate(subgrid_EMEP_local_from_in_region(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
        subgrid_EMEP_local_from_in_region = 0.0
        if (allocated(subgrid_EMEP_semilocal_from_in_region)) deallocate(subgrid_EMEP_semilocal_from_in_region)
        allocate(subgrid_EMEP_semilocal_from_in_region(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop))
        subgrid_EMEP_semilocal_from_in_region = 0.0

        write(unit_logfile,'(A)') 'Interpolating local and semilocal contribution from-in-region to the target subgrid'

        ! go through all target subgrids and interpolate the EMEP contributions calculated above
        do i = 1, subgrid_dim(x_dim_index)
            do j = 1, subgrid_dim(y_dim_index)

                ! Find position along the region dimension of this region index
                if (subgrid_region_index(i,j) > 0) then
                    i_region = regionindex_loop_back_index(subgrid_region_index(i,j))
                else
                    ! this subgrid is not in any region, so keep it as 0
                    cycle
                end if

                ! Find which EMEP grid the current subgrid is in
                ii=crossreference_target_to_emep_subgrid(i,j,x_dim_index)
                jj=crossreference_target_to_emep_subgrid(i,j,y_dim_index)

                ! Find out where in this EMEP grid we are in EMEP grid coordinates
                ! E.g. (0,0) for lower-left corner and (1,1) for upper right corner
                call LL2PROJ(lon_subgrid(i,j),lat_subgrid(i,j),x_temp,y_temp,EMEP_projection_attributes,EMEP_projection_type)
                ! fractional position inside the EMEP grid (center is (0,0), lower-left corner is (-0.5,-0.5), upper-right corner is (+0.5,+0.5))
                ii_frac_target = (x_temp-var1d_nc(ii,x_dim_nc_index))/dgrid_nc(x_dim_nc_index)
                jj_frac_target = (y_temp-var1d_nc(jj,y_dim_nc_index))/dgrid_nc(y_dim_nc_index)
                ! verify that this is within 0-1 (if not, crossreference has gone wrong...)
                if (ii_frac_target < -0.5 .or. ii_frac_target > 0.5 .or. jj_frac_target < -0.5 .or. jj_frac_target > 0.5) then
                    write(unit_logfile,'(A,2I12)') 'Something went wrong with locating target subgrid within EMEP grid!',ii_frac_target,jj_frac_target
                    stop
                end if

                ! Interpolate to the target grid
                if (EMEP_grid_interpolation_flag == 0) then
                    ! No interpolation: just pick the value of the EMEP cell we are in
                    subgrid_EMEP_local_from_in_region(i,j,:,:,:) = EMEP_local_from_in_region(ii,jj,i_region,:,:,:)
                    subgrid_EMEP_semilocal_from_in_region(i,j,:,:,:) = EMEP_semilocal_from_in_region(ii,jj,i_region,:,:,:)
                else if (EMEP_grid_interpolation_flag == 6) then
                    ! Use area-weighted interpolation
                    weight_check = 0.0
                    ! loop over a 3x3 EMEP cell domain centered at the closest EMEP cell to the target subgrid
                    do iii = -1,1
                        do jjj = -1,1
                            iii_nc = ii + iii
                            jjj_nc = jj + jjj
                            ! verify the EMEP grid covers this index
                            if (.not. (iii_nc >= 1 .and. iii_nc <= dim_length_nc(x_dim_nc_index) .and. jjj_nc >= 1 .and. jjj_nc <= dim_length_nc(y_dim_nc_index))) then
                                write(unit_logfile,'(A)') 'ERROR: EMEP grid did not go far enough out to allow area interpolation to the target grid!'
                                stop
                            end if
                            ! calculate the weighting value as the area fraction of an EMEP cell centered at the target subgrid that falls within this EMEP cell
                            x_dist_target = iii - ii_frac_target
                            y_dist_target = jjj - jj_frac_target
                            weighting_val = max(0.0, (1.0 - abs(x_dist_target))) * max(0.0, (1.0 - abs(y_dist_target)))
                            weight_check = weight_check + weighting_val
                            ! use this weighting for the data at this EMEP cell
                            if (weighting_val > 0) then
                                subgrid_EMEP_local_from_in_region(i,j,:,:,:) = subgrid_EMEP_local_from_in_region(i,j,:,:,:) + EMEP_local_from_in_region(iii_nc,jjj_nc,i_region,:,:,:)*weighting_val
                                subgrid_EMEP_semilocal_from_in_region(i,j,:,:,:) = subgrid_EMEP_semilocal_from_in_region(i,j,:,:,:) + EMEP_semilocal_from_in_region(iii_nc,jjj_nc,i_region,:,:,:)*weighting_val
                            end if
                        end do
                    end do
                else
                    write(unit_logfile,'(A,I0)') 'ERROR: uEMEP_subgrid_EMEP_from_in_region is not implemented for EMEP_grid_interpolation_flag =', EMEP_grid_interpolation_flag
                    stop
                end if
            end do
        end do

        deallocate(EMEP_local_from_in_region)
        deallocate(EMEP_semilocal_from_in_region)

        ! Set allsources to be the sum of only the sources we calculate for
        subgrid_EMEP_local_from_in_region(:,:,:,allsource_index,:)=0
        subgrid_EMEP_semilocal_from_in_region(:,:,:,allsource_index,:)=0
        do i_source = 1, n_source_index
            if (calculate_source(i_source).or.calculate_EMEP_source(i_source).and.i_source.ne.allsource_index) then
                subgrid_EMEP_local_from_in_region(:,:,:,allsource_index,:)=subgrid_EMEP_local_from_in_region(:,:,:,allsource_index,:)+subgrid_EMEP_local_from_in_region(:,:,:,i_source,:)
                subgrid_EMEP_semilocal_from_in_region(:,:,:,allsource_index,:)=subgrid_EMEP_semilocal_from_in_region(:,:,:,allsource_index,:)+subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,:)
            end if
        end do

    end subroutine uEMEP_subgrid_EMEP_from_in_region