uEMEP_combine_local_source Subroutine

public subroutine uEMEP_combine_local_source()

Uses

  • proc~~uemep_combine_local_source~~UsesGraph proc~uemep_combine_local_source uEMEP_combine_local_source module~uemep_definitions uEMEP_definitions proc~uemep_combine_local_source->module~uemep_definitions

Arguments

None

Calls

proc~~uemep_combine_local_source~~CallsGraph proc~uemep_combine_local_source uEMEP_combine_local_source proc~mean_mask mean_mask proc~uemep_combine_local_source->proc~mean_mask

Called by

proc~~uemep_combine_local_source~~CalledByGraph proc~uemep_combine_local_source uEMEP_combine_local_source program~uemep uEMEP program~uemep->proc~uemep_combine_local_source

Source Code

    subroutine uEMEP_combine_local_source

        use uEMEP_definitions

        implicit none

        integer source_index
        integer i_pollutant
        integer i,j
        integer i_sp
        real sum_temp(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index))

        !if (interpolate_subgrids_flag) then
        !write(unit_logfile,'(a)') 'Interpolate routines not currently active. Doing nothing'
        !call uEMEP_interpolate_auto_subgrid
        !return
        !stop
        !call uEMEP_interpolate_subgrids
        !call uEMEP_linear_interpolate_subgrids
        !call uEMEP_bilinear_interpolate_subgrids

        !Remember to reset the use_subgrids val and logical so that everything will be used in the end
        !endif

        write(unit_logfile,'(a)')'--------------------------'
        write(unit_logfile,'(a)')'Combining the local and nonlocal contributions at each subgrid'
        write(unit_logfile,'(a)')'--------------------------'

        !Calculate redistributed subgrid allsource concentrations
        subgrid(:,:,:,local_subgrid_index,allsource_index,:)=0.
        do source_index=1,n_source_index
            if (calculate_source(source_index).and.source_index.ne.allsource_index) then
                subgrid(:,:,:,local_subgrid_index,allsource_index,:)=subgrid(:,:,:,local_subgrid_index,allsource_index,:)+subgrid(:,:,:,local_subgrid_index,source_index,:)
            endif
            !Add the selected EMEP local sources to this as well, if they are not already included in the subgrid downscaling
            if (calculate_EMEP_source(source_index).and..not.calculate_source(source_index).and.source_index.ne.allsource_index) then
                subgrid(:,:,:,local_subgrid_index,allsource_index,:)=subgrid(:,:,:,local_subgrid_index,allsource_index,:)+subgrid(:,:,:,emep_local_subgrid_index,source_index,:)
                !write(*,*) source_index,sum(subgrid(:,:,:,emep_local_subgrid_index,source_index,:)),sum(subgrid(:,:,:,local_subgrid_index,allsource_index,:)),sum(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,:))
            endif
        enddo

        if (EMEP_additional_grid_interpolation_size.gt.0) then
            do i_pollutant=1,n_pollutant_loop
                !If the compound is PM2.5 or PM10 then add the non PPM part to the non-local
                if (pollutant_loop_index(i_pollutant).eq.pm25_index.or.pollutant_loop_index(i_pollutant).eq.pm10_index) then
                    write(unit_logfile,'(A,A)') 'Pollutant: ',trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_nc_index))
                    write(unit_logfile,'(A,f12.2)') 'MEAN PPM NONLOCAL ADDITIONAL: ',sum((subgrid(:,:,:,emep_additional_nonlocal_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                    subgrid(:,:,:,emep_additional_nonlocal_subgrid_index,allsource_index,i_pollutant)=subgrid(:,:,:,emep_additional_nonlocal_subgrid_index,allsource_index,i_pollutant) &
                        +(comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant))
                    write(unit_logfile,'(A,f12.2)') 'MEAN ADD REST NONLOCAL ADDITIONAL: ',sum((comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                endif

                !subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)=subgrid(:,:,:,local_subgrid_index,allsource_index,i_pollutant)+subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)

            enddo
        endif

        do i_pollutant=1,n_pollutant_loop
            !If the compound is PM2.5 or PM10 then add the non PPM part to the non-local
            if (pollutant_loop_index(i_pollutant).eq.pm25_index.or.pollutant_loop_index(i_pollutant).eq.pm10_index) then
                write(unit_logfile,'(A,A)') 'Pollutant: ',trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_nc_index))
                write(unit_logfile,'(A,f12.2)') 'MEAN PPM NONLOCAL: ',mean_mask(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                write(unit_logfile,'(A,f12.2)') 'MEAN PPM LOCAL: ',mean_mask(subgrid(:,:,:,emep_local_subgrid_index,allsource_index,i_pollutant),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((subgrid(:,:,:,emep_local_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                write(unit_logfile,'(A,f12.2)') 'MEAN PPM TOTAL: ',mean_mask(subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                write(unit_logfile,'(A,f12.2)') 'MEAN COMP PM TOTAL: ',mean_mask(comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant)),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                write(unit_logfile,'(A,f12.2)') 'MEAN COMP PM ORIGINAL: ',mean_mask(orig_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant)),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((orig_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)=subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant) &
                    +(comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant))
                write(unit_logfile,'(A,f12.2)') 'MEAN ADD REST NONLOCAL: ',mean_mask(comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                write(unit_logfile,'(A,f12.2)') 'MEAN NEW PM NONLOCAL: ',mean_mask(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant),use_subgrid(:,:,allsource_index),subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)) !sum((subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                if (sum(comp_EMEP_subgrid(:,:,:,pollutant_loop_index(i_pollutant))-subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant)).lt.0) write(unit_logfile,'(A)') 'WARNING!!!: PPM EMEP is more than total EMEP PM. Negative non PPM contributions.'
            endif

            subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)=subgrid(:,:,:,local_subgrid_index,allsource_index,i_pollutant)+subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)
        enddo

        do i_pollutant=1,n_pollutant_loop
            !Place the results in the compound results
            !do i_loop=1,n_pollutant_compound_loop(i_pollutant)
            comp_subgrid(:,:,:,pollutant_loop_index(i_pollutant))=subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)
            !enddo
            !write(*,'(2i,4f16.2)') i_pollutant,pollutant_loop_index(i_pollutant)&
            !    ,sum(subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)&
            !    ,sum(subgrid(:,:,:,local_subgrid_index,allsource_index,i_pollutant))&
            !    ,sum(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant))&
            !    ,sum(comp_subgrid(:,:,:,pollutant_loop_index(i_pollutant)))

        enddo

        !Replace the species PPM with the actual species PPM used in the local fraction and then replace this with the nonlocal part.
        !Only works for the complete group of species
        if (save_emep_species) then
            !Replace the primary species value with the one used in the calculations for consistency
            write(unit_logfile,'(A)') 'ppm read from surf and read emep. Difference should just be deposition'
            write(unit_logfile,'(A,2f12.2)') 'PPM25 (init_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(subgrid(:,:,:,emep_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_nc_index)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            write(unit_logfile,'(A,2f12.2)') 'PPM10 (init_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(subgrid(:,:,:,emep_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_nc_index)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            write(unit_logfile,'(A)') 'pm surf init species and summed species should be the same. comp should be different (difference between SURF and SURF_rh50)'
            write(unit_logfile,'(A,3f12.2)') 'PM25 (init_sp,sum_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm25_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            write(unit_logfile,'(A,3f12.2)') 'PM10 (init_sp,sum_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm10_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            !The remaining species will not, most likely, be normalised with the total. This is done here but the total ppm is kept
            !Something weird is going on here. Need to check

            !Normalise the species surface values with the total grid values to account for any difference du to deposition and water
            !Because t can have dimension 1 and 'sum' function does not think it is a dimension anymore then need to do the sum over the time loop. This turned out to not be the case so they are the same
            if (use_single_time_loop_flag) then
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index),4)
            else
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index),4)
            endif
            do i_sp=1,sp_ppm_index
                species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp) &
                    *(comp_EMEP_subgrid(:,:,:,pm25_nc_index)) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
            enddo
            if (use_single_time_loop_flag) then
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index),4)
            else
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index),4)
            endif
            do i_sp=1,sp_ppm_index
                species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp) &
                    *(comp_EMEP_subgrid(:,:,:,pm10_nc_index)) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
            enddo

            write(unit_logfile,'(A)') 'sum_sp should be the same as comp after scaling, init should be unchanged'
            write(unit_logfile,'(A,3f12.2)') 'PM25 (init_sp,sum_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm25_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            write(unit_logfile,'(A,3f12.2)') 'PM10 (init_sp,sum_sp,comp)',sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm10_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            !Replace the ppm values in the species with the emep read version
            !This is not correct. emep_subgrid_index,allsource_index is the sum of all local and nonlocal sources, so a larger number than the primary.
            !This will lead to negative values set to 0. So comment out this line (21.03.2022)
            !species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index)=subgrid(:,:,:,emep_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_nc_index))
            !species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index)=subgrid(:,:,:,emep_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_nc_index))
            !Normalise the species other than ppm again after subtracting the ppm
            !This should not change anything
            if (use_single_time_loop_flag) then
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index-1),4)
            else
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index-1),4)
            endif
            do i_sp=1,sp_ppm_index-1
                species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp) &
                    *(comp_EMEP_subgrid(:,:,:,pm25_nc_index)-species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index)) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
            enddo
            if (use_single_time_loop_flag) then
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index-1),4)
            else
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index-1),4)
            endif
            do i_sp=1,sp_ppm_index-1
                species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp) &
                    *(comp_EMEP_subgrid(:,:,:,pm10_nc_index)-species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index)) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
            enddo


            !Remove the primary based on the local contribution.
            species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index)=species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index)-subgrid(:,:,:,emep_local_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_nc_index))
            species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index)=species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index)-subgrid(:,:,:,emep_local_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_nc_index))
            where (species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index).lt.0.or.isnan(species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index))) species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_ppm_index)=0
            where (species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index).lt.0.or.isnan(species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index))) species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_ppm_index)=0

            !Normalise again after the subtraction to the nonlocal contribution (21.03.2022)
            sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index),4)
            do i_sp=1,sp_ppm_index
                species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp) &
                    *(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_nc_index))) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm25_sp_index,i_sp)=0
            enddo
            if (use_single_time_loop_flag) then
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index),4)
            else
                sum_temp(:,:,:)=sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index),4)
            endif
            do i_sp=1,sp_ppm_index
                species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp) &
                    *(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_nc_index))) &
                    /sum_temp
                where (sum_temp.le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
                where (isnan(species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)).or.species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp).le.0) species_EMEP_subgrid(:,:,:,pm10_sp_index,i_sp)=0
            enddo


            write(unit_logfile,'(A)') 'init_sp should be the same as before, sum_sp should be the same as nonlocal, comp should be larger'
            write(unit_logfile,'(A,4f12.2)') 'PM25 (init_sp,sum_sp,comp,nonlocal)',sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm25_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm25_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_nc_index)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
            write(unit_logfile,'(A,4f12.2)') 'PM10 (init_sp,sum_sp,comp,nonlocal)',sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,sp_pm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(species_EMEP_subgrid(:,:,:,pm10_sp_index,1:sp_ppm_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(comp_EMEP_subgrid(:,:,:,pm10_nc_index))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index) &
                ,sum(subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_nc_index)))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
        endif


        !Only show results where all the subgrids and all sources are valid
        !This is temporary
        if (interpolate_subgrids_flag) then
            do source_index=1,n_source_index
                if (calculate_source(source_index)) then
                    do j=1,subgrid_dim(y_dim_index)
                        do i=1,subgrid_dim(x_dim_index)
                            if (.not.use_subgrid(i,j,source_index)) then
                                !subgrid(i,j,:,total_subgrid_index,allsource_index,:)=NODATA_value
                                !subgrid(i,j,:,total_subgrid_index,source_index,:)=NODATA_value
                                !comp_subgrid(i,j,:,:)=NODATA_value
                            endif
                        enddo
                    enddo
                endif
            enddo
        endif

    end subroutine uEMEP_combine_local_source