uEMEP_save_netcdf_control Subroutine

public subroutine uEMEP_save_netcdf_control()

Uses

  • proc~~uemep_save_netcdf_control~~UsesGraph proc~uemep_save_netcdf_control uEMEP_save_netcdf_control module~uemep_definitions uEMEP_definitions proc~uemep_save_netcdf_control->module~uemep_definitions

.or. use_region_select_and_mask_flag ?? !!!!! what to write here???

Arguments

None

Calls

proc~~uemep_save_netcdf_control~~CallsGraph proc~uemep_save_netcdf_control uEMEP_save_netcdf_control proc~direction DIRECTION proc~uemep_save_netcdf_control->proc~direction proc~mean_mask mean_mask proc~uemep_save_netcdf_control->proc~mean_mask proc~uemep_save_netcdf_file uEMEP_save_netcdf_file proc~uemep_save_netcdf_control->proc~uemep_save_netcdf_file proc~uemep_save_netcdf_receptor_file uEMEP_save_netcdf_receptor_file proc~uemep_save_netcdf_control->proc~uemep_save_netcdf_receptor_file proc~uemep_source_fraction_chemistry uEMEP_source_fraction_chemistry proc~uemep_save_netcdf_control->proc~uemep_source_fraction_chemistry proc~write_esri_ascii_file write_esri_ascii_file proc~uemep_save_netcdf_control->proc~write_esri_ascii_file nf90_close nf90_close proc~uemep_save_netcdf_file->nf90_close nf90_create nf90_create proc~uemep_save_netcdf_file->nf90_create nf90_def_dim nf90_def_dim proc~uemep_save_netcdf_file->nf90_def_dim nf90_def_var nf90_def_var proc~uemep_save_netcdf_file->nf90_def_var nf90_def_var_chunking nf90_def_var_chunking proc~uemep_save_netcdf_file->nf90_def_var_chunking nf90_def_var_deflate nf90_def_var_deflate proc~uemep_save_netcdf_file->nf90_def_var_deflate nf90_enddef nf90_enddef proc~uemep_save_netcdf_file->nf90_enddef nf90_inq_dimid nf90_inq_dimid proc~uemep_save_netcdf_file->nf90_inq_dimid nf90_inq_varid nf90_inq_varid proc~uemep_save_netcdf_file->nf90_inq_varid nf90_inquire_dimension nf90_inquire_dimension proc~uemep_save_netcdf_file->nf90_inquire_dimension nf90_open nf90_open proc~uemep_save_netcdf_file->nf90_open nf90_put_att nf90_put_att proc~uemep_save_netcdf_file->nf90_put_att nf90_put_var nf90_put_var proc~uemep_save_netcdf_file->nf90_put_var nf90_redef nf90_redef proc~uemep_save_netcdf_file->nf90_redef proc~check check proc~uemep_save_netcdf_file->proc~check proc~uemep_save_netcdf_receptor_file->nf90_close proc~uemep_save_netcdf_receptor_file->nf90_create proc~uemep_save_netcdf_receptor_file->nf90_def_dim proc~uemep_save_netcdf_receptor_file->nf90_def_var proc~uemep_save_netcdf_receptor_file->nf90_enddef proc~uemep_save_netcdf_receptor_file->nf90_inq_dimid proc~uemep_save_netcdf_receptor_file->nf90_inq_varid proc~uemep_save_netcdf_receptor_file->nf90_inquire_dimension proc~uemep_save_netcdf_receptor_file->nf90_open proc~uemep_save_netcdf_receptor_file->nf90_put_att proc~uemep_save_netcdf_receptor_file->nf90_put_var proc~uemep_save_netcdf_receptor_file->nf90_redef proc~area_weighted_interpolation_function area_weighted_interpolation_function proc~uemep_save_netcdf_receptor_file->proc~area_weighted_interpolation_function proc~uemep_save_netcdf_receptor_file->proc~check proc~uemep_during_no2 uEMEP_During_NO2 proc~uemep_source_fraction_chemistry->proc~uemep_during_no2 proc~uemep_nonlocal_no2_o3 uEMEP_nonlocal_NO2_O3 proc~uemep_source_fraction_chemistry->proc~uemep_nonlocal_no2_o3 proc~uemep_photostationary_no2 uEMEP_photostationary_NO2 proc~uemep_source_fraction_chemistry->proc~uemep_photostationary_no2 proc~uemep_phototimescale_no2 uEMEP_phototimescale_NO2 proc~uemep_source_fraction_chemistry->proc~uemep_phototimescale_no2 proc~uemep_romberg_no2 uEMEP_Romberg_NO2 proc~uemep_source_fraction_chemistry->proc~uemep_romberg_no2 proc~uemep_srm_no2 uEMEP_SRM_NO2 proc~uemep_source_fraction_chemistry->proc~uemep_srm_no2 nf90_strerror nf90_strerror proc~check->nf90_strerror clog clog proc~uemep_phototimescale_no2->clog

Called by

proc~~uemep_save_netcdf_control~~CalledByGraph proc~uemep_save_netcdf_control uEMEP_save_netcdf_control program~uemep uEMEP program~uemep->proc~uemep_save_netcdf_control

Source Code

    subroutine uEMEP_save_netcdf_control

        use uEMEP_definitions

        implicit none

        integer i,j,k,t,l
        integer i_comp,i_file,i_meteo
        character(256) temp_name,unit_str,title_str,title_str_rec,var_name_temp, temp_date_str,station_name_str,temp_name_rec,temp_compound_str
        logical create_file,create_file_rec
        integer i_source
        real :: valid_min=0.
        real, allocatable :: temp_subgrid(:,:,:)
        real, allocatable :: temp_integral_subgrid(:,:,:)
        real, allocatable :: aqi_subgrid(:,:,:,:)
        integer, allocatable :: aqi_responsible_pollutant_index(:,:,:)
        real, allocatable :: temp_subgrid_ascii(:,:)

        integer ii,jj,tt

        real aqi_limits_temp(n_compound_index,1:5)
        real max_aqi
        integer n_aqi_pollutant_index
        parameter (n_aqi_pollutant_index=4)
        integer aqi_pollutant_index(n_aqi_pollutant_index)
        integer i_pollutant,i_loop
        character(256) variable_type
        logical :: receptor_available=.true.
        real scale_factor
        integer n_save_aqi_pollutant_index
        real temp_sum_comp
        integer count
        character(256) filename_ascii

        integer i_sp,ii_sp

        integer source_domain_loop
        integer no2_o3_loop,comp_index,comp_nc_index
        character(256) filename_append

        if (include_o3_in_aqi_index) then
            n_save_aqi_pollutant_index=n_aqi_pollutant_index
        else
            n_save_aqi_pollutant_index=n_aqi_pollutant_index-1
        endif

        if (.not.allocated(temp_subgrid)) allocate(temp_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)))
        if (.not.allocated(temp_integral_subgrid)) allocate(temp_integral_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)))
        if (.not.allocated(aqi_subgrid).and.save_aqi) allocate(aqi_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index))
        if (.not.allocated(aqi_responsible_pollutant_index).and.save_aqi) allocate(aqi_responsible_pollutant_index(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index)))
        if (.not.allocated(temp_subgrid_ascii).and.save_compounds_as_ascii) allocate(temp_subgrid_ascii(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index)))

        !Save subgrid calculations
        valid_min=0.
        unit_str="ug/m3"
        variable_type='byte'
        !variable_type='double'
        variable_type='float'
        scale_factor=1.

        !Save the data
        i_file=subgrid_total_file_index(allsource_index)
        !temp_name=trim(pathname_grid(i_file))//trim(filename_grid(i_file))//'_'//trim(file_tag)//'.nc'
        !if (len(config_date_str).gt.0) then
        !    temp_date_str='_'//trim(config_date_str)
        !else
        !    temp_date_str=''
        !endif

        if (len(filename_date_output_grid).gt.0) then
            temp_date_str='_'//filename_date_output_grid
        else
            temp_date_str=''
        endif

        if (use_multiple_receptor_grids_flag) then
            station_name_str=trim(name_receptor_in(g_loop,1))//'_';
        else
            station_name_str=''
        endif

        !Do not write 'all' in file name if all compounds are selected
        if (pollutant_index.eq.all_nc_index) then
            temp_compound_str=''
        else
            temp_compound_str='_'//trim(var_name_nc(conc_nc_index,compound_index,allsource_index))
        endif
        !Do not write anything about the compounds
        temp_compound_str=''

        !temp_name=trim(pathname_grid(i_file))//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//trim(temp_date_str)//'_'//trim(forecast_hour_str)//'.nc'
        !temp_name_rec=trim(pathname_grid(i_file))//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//trim(temp_date_str)//'_'//trim(forecast_hour_str)//'.nc'
        !if (save_netcdf_average_flag) then
        !temp_name=trim(pathname_grid(i_file))//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'_'//trim(forecast_hour_str)//'.nc'
        !temp_name_rec=trim(pathname_grid(i_file))//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'_'//trim(forecast_hour_str)//'.nc'
        !endif
        temp_name=trim(pathname_grid(i_file))//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//trim(temp_date_str)//'.nc'
        temp_name_rec=trim(pathname_grid(i_file))//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//trim(temp_date_str)//'.nc'
        if (save_netcdf_average_flag) then
            temp_name=trim(pathname_grid(i_file))//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'.nc'
            temp_name_rec=trim(pathname_grid(i_file))//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'.nc'
        endif

        finished_file=trim(pathname_grid(i_file))//trim(finished_subpath)//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//trim(temp_date_str)//'.'//trim(finished_filename)
        finished_file_rec=trim(pathname_grid(i_file))//trim(finished_subpath)//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//trim(temp_date_str)//'.'//trim(finished_filename)
        if (save_netcdf_average_flag) then
            finished_file=trim(pathname_grid(i_file))//trim(finished_subpath)//trim(station_name_str)//'uEMEP_'//trim(file_tag)//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'.'//trim(finished_filename)
            finished_file_rec=trim(pathname_grid(i_file))//trim(finished_subpath)//'uEMEP_'//trim(file_tag)//'_station'//trim(temp_compound_str)//'_mean'//trim(temp_date_str)//'.'//trim(finished_filename)
        endif

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '================================================================'
        write(unit_logfile,'(A)') 'Saving netcdf data (uEMEP_save_netcdf_control)'
        write(unit_logfile,'(A)') '================================================================'

        i_comp=1
        if (save_netcdf_receptor_flag.and.n_valid_receptor.eq.0) then
            if (i_comp.eq.1.and.t_loop.eq.start_time_loop_index) then
                write(unit_logfile,'(a)')'No receptor positions available. Will not save receptor data.'
                receptor_available=.false.
            endif
        endif

        if (save_netcdf_average_flag) then
            if (.not.allocated(val_array_av)) then
                allocate(val_array_av(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),n_var_av))
                val_array_av=0.
            endif
            if (.not.allocated(time_seconds_output_av)) then
                allocate(time_seconds_output_av(n_var_av))
                time_seconds_output_av=0
            endif
            !Reset average file counter for every time step saving occurs
            counter_av=0
            write(unit_logfile,*) 'Saving as mean data'
        endif

        !Save the final result of the subgrid calculation in ascii format
        !This should only be used for annual mean calculations since it cannot have a time dimmension
        !Intended for FAIRMODE output
        !Makes a new file for each compound and averages over time, if there is a time dimension
        if (save_compounds_as_ascii) then
            variable_type='float'
            unit_str="ug/m3"
            do i_pollutant=1,n_pollutant_loop
                if (pollutant_loop_index(i_pollutant).ne.pmex_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_sand_index.and.pollutant_loop_index(i_pollutant).ne.pm10_sand_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_salt_index.and.pollutant_loop_index(i_pollutant).ne.pm10_salt_index) then
                    do i_loop=1,n_pollutant_compound_loop(i_pollutant)

                        i_comp=pollutant_compound_loop_index(i_pollutant,i_loop)
                        !write(*,*) i_comp

                        var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))

                        title_str=trim(var_name_temp)//'_'//trim(file_tag)//trim(temp_date_str)
                        write(unit_logfile,'(a)')'Writing ascii data to: '//trim(title_str)

                        filename_ascii=trim(pathname_output_grid)//trim(title_str)//'.asc'

                        temp_subgrid_ascii(:,:)=sum(comp_subgrid(:,:,:,i_comp),3)/subgrid_dim(t_dim_index)


                        write(unit_logfile,'(a,f12.3)')'Writing ascii array variable: '//trim(var_name_temp),sum(comp_subgrid(:,:,:,i_comp))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)

                        call write_esri_ascii_file(filename_ascii,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_delta(x_dim_index),temp_subgrid_ascii,x_subgrid,y_subgrid)

                    enddo
                endif
            enddo
        endif

        !Save the final result of the subgrid calculation
        if (save_compounds) then

            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving all total compounds'
            write(unit_logfile,'(a)')'--------------------------'

            variable_type='float'
            unit_str="ug/m3"
            do i_pollutant=1,n_pollutant_loop
                if (pollutant_loop_index(i_pollutant).ne.pmex_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_sand_index.and.pollutant_loop_index(i_pollutant).ne.pm10_sand_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_salt_index.and.pollutant_loop_index(i_pollutant).ne.pm10_salt_index) then
                    do i_loop=1,n_pollutant_compound_loop(i_pollutant)

                        i_comp=pollutant_compound_loop_index(i_pollutant,i_loop)
                        !write(*,*) i_comp

                        if (i_pollutant.eq.1.and.i_loop.eq.1.and.t_loop.eq.start_time_loop_index.and.save_netcdf_file_flag) then
                            create_file=.true.
                            title_str='uEMEP_concentration_'//trim(file_tag)//temp_date_str
                            write(unit_logfile,'(a)')'Writing to: '//trim(temp_name)
                        else
                            create_file=.false.
                        endif

                        if (i_pollutant.eq.1.and.i_loop.eq.1.and.t_loop.eq.start_time_loop_index.and.first_g_loop.and.save_netcdf_receptor_flag) then
                            create_file_rec=.true.
                            title_str_rec='uEMEP_receptor_'//trim(file_tag)//temp_date_str
                            if (receptor_available) write(unit_logfile,'(a)')'Writing to: '//trim(temp_name_rec)
                        else
                            create_file_rec=.false.
                        endif

                        var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_concentration'
                        if (save_netcdf_file_flag) then
                            temp_sum_comp=0.
                            count=0
                            do j=1,subgrid_dim(y_dim_index)
                                do i=1,subgrid_dim(x_dim_index)
                                    if (use_subgrid(i,j,allsource_index)) then
                                        temp_sum_comp=temp_sum_comp+sum(comp_subgrid(i,j,:,i_comp))/subgrid_dim(t_dim_index)
                                        count=count+1
                                    endif
                                enddo
                            enddo
                            if (count.gt.0) then
                                temp_sum_comp=temp_sum_comp/count
                            else
                                temp_sum_comp=0
                            endif
                            !write(unit_logfile,'(a,f12.3)')'Writing netcdf array variable:    '//trim(var_name_temp),temp_sum_comp
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(comp_subgrid(:,:,:,i_comp),use_subgrid(:,:,allsource_index),size(comp_subgrid(:,:,:,i_comp),1),size(comp_subgrid(:,:,:,i_comp),2),size(comp_subgrid(:,:,:,i_comp),3))
                            call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,comp_subgrid(:,:,:,i_comp),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                        endif
                        if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(comp_subgrid(:,:,:,i_comp),use_subgrid(:,:,allsource_index),size(comp_subgrid(:,:,:,i_comp),1),size(comp_subgrid(:,:,:,i_comp),2),size(comp_subgrid(:,:,:,i_comp),3))
                            !write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(comp_subgrid(:,:,:,i_comp))/subgrid_dim(x_dim_index)/subgrid_dim(y_dim_index)/subgrid_dim(t_dim_index)
                            call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,comp_subgrid(:,:,:,i_comp),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str_rec,create_file_rec,valid_min &
                                ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,z_rec(allsource_index,1) &
                                ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                        endif
                    enddo
                endif
            enddo

        endif

        create_file=.false.
        create_file_rec=.false.

        ! Save downscaled source contributions for pollutants (not no2 or o3)
        ! 1 = local
        ! 2 = local, cut down to from_in_region
        ! 3-4: not used here, but kept to make numbering consistent with the other loops over source_domain_loop
        ! 5 = total: contribution from big domain, but using downscaled within moving window: 1 + (additional_emep_local(3) - emep_local(1)
        ! 6 = total from in region = 2 + semilocal emep contribution (4))
        do source_domain_loop = 1,6

            ! Skip this round in the loop if this type of source contributions are not to be saved
            if (source_domain_loop == 1) then
                if (.not. save_source_contributions) cycle
                filename_append=''
            else if (source_domain_loop == 2) then
                if (.not. (trace_emissions_from_in_region .and. save_local_source_contributions_from_in_region)) cycle
                filename_append='_from_in_region'
            else if (source_domain_loop == 5) then
                if (.not. (save_total_source_contributions .and. EMEP_additional_grid_interpolation_size.gt.0)) cycle
                filename_append=''
            else if (source_domain_loop == 6) then
                if (.not. (trace_emissions_from_in_region .and. save_total_source_contributions_from_in_region)) cycle
                filename_append='_from_in_region'
            else
                ! skip rounds 3-4
                cycle
            end if

            write(unit_logfile,'(a)')'--------------------------'
            if (source_domain_loop == 1) write(unit_logfile,'(a)')'Saving downscaled contributions'
            if (source_domain_loop == 2) write(unit_logfile,'(a)')'Saving downscaled contributions from-in-region'
            if (source_domain_loop == 5) write(unit_logfile,'(a)')'Saving total contributions'
            if (source_domain_loop == 6) write(unit_logfile,'(a)')'Saving total contributions from-in-region'
            write(unit_logfile,'(a)')'--------------------------'

            if (save_netcdf_fraction_as_contribution_flag) then
                variable_type='float'
                unit_str="ug/m3"
            else
                variable_type='byte'
                unit_str="%"
            endif

            do i_pollutant=1,n_pollutant_loop
                do i_source=1,n_source_index
                    if (source_domain_loop == 5) then
                        i_file=subgrid_sourcetotal_file_index(i_source)
                    else if (source_domain_loop == 6) then
                        i_file=subgrid_sourcetotal_inregion_file_index(i_source)
                    else  ! 1 or 2
                        i_file=subgrid_local_file_index(i_source)
                    end if
                    !if (calculate_source(i_source).or.i_source.eq.allsource_index) then
                    !Don't save any exhaust, sand or salt pollutant sources. Dealt with later
                    !(calculate_source(i_source).or.calculate_emep_source(i_source))
                    if (calculate_source(i_source).and.i_source.ne.allsource_index.and.pollutant_loop_index(i_pollutant).ne.pmex_index &
                        .and.pollutant_loop_index(i_pollutant).ne.pm25_sand_index.and.pollutant_loop_index(i_pollutant).ne.pm10_sand_index &
                        .and.pollutant_loop_index(i_pollutant).ne.pm25_salt_index.and.pollutant_loop_index(i_pollutant).ne.pm10_salt_index) then

                        !Only save nonexhaust pm and all the other sources
                        if ((pollutant_loop_index(i_pollutant).eq.nox_index.and.i_source.eq.traffic_index)) then
                            !if ((pollutant_loop_index(i_pollutant).ne.nox_index.or.i_source.ne.traffic_index)) then
                        else
                            if (i_source.eq.traffic_index.and.(pollutant_loop_index(i_pollutant).eq.pm10_index.or.pollutant_loop_index(i_pollutant).eq.pm25_index)) then
                                ! Special case for traffic pm: In some setups we need to subtract exhaust to get non-exhaust
                                if (pollutant_index.eq.all_totals_nc_index.or.pollutant_index.eq.aaqd_totals_nc_index.or.pollutant_index.eq.gp_totals_nc_index.or.pollutant_index.eq.op_totals_nc_index) then
                                    ! Traffic exhaust is not a separate pollutant, so just save total traffic
                                    var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//trim(filename_append)
                                    if (source_domain_loop == 1) then
                                        temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)
                                    else if (source_domain_loop == 2) then
                                        temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)
                                    else if (source_domain_loop == 5) then
                                        temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)-subgrid(:,:,:,emep_local_subgrid_index,i_source,i_pollutant)+subgrid(:,:,:,emep_additional_local_subgrid_index,i_source,i_pollutant)
                                    else ! source_domain_loop == 6
                                        temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)+subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,i_pollutant)
                                    endif
                                else
                                    ! Calculate non-exhaust as difference between total traffic and exhaust for downscaled contributions
                                    ! For EMEP contributions, non-exhaust can be directly accessed with traffic_nonexhaust_nc_index
                                    var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//'_nonexhaust'//trim(filename_append)
                                    if (source_domain_loop == 1) then
                                        temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)-subgrid(:,:,:,local_subgrid_index,i_source,pollutant_loop_back_index(pmex_index))
                                    else if (source_domain_loop == 2) then
                                        temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)-subgrid_local_from_in_region(:,:,:,i_source,pollutant_loop_back_index(pmex_index))
                                    else if (source_domain_loop == 5) then
                                        temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)-subgrid(:,:,:,local_subgrid_index,i_source,pollutant_loop_back_index(pmex_index))-subgrid(:,:,:,emep_local_subgrid_index,traffic_nonexhaust_nc_index,i_pollutant)+subgrid(:,:,:,emep_additional_local_subgrid_index,traffic_nonexhaust_nc_index,i_pollutant)
                                    else  !source_domain_loop == 6
                                        temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)-subgrid_local_from_in_region(:,:,:,i_source,pollutant_loop_back_index(pmex_index))+subgrid_EMEP_semilocal_from_in_region(:,:,:,traffic_nonexhaust_nc_index,i_pollutant)
                                    end if
                                endif
                            else
                                var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//trim(filename_append)
                                if (source_domain_loop == 1) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)
                                else if (source_domain_loop == 2) then
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)
                                else if (source_domain_loop == 5) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)-subgrid(:,:,:,emep_local_subgrid_index,i_source,i_pollutant)+subgrid(:,:,:,emep_additional_local_subgrid_index,i_source,i_pollutant)
                                else ! source_domain_loop == 6
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)+subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,i_pollutant)
                                endif
                            endif
                            ! Convert to fractions
                            if (.not. save_netcdf_fraction_as_contribution_flag) then
                                temp_subgrid=temp_subgrid/subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)*100.
                            end if
                            !In case of any 0 concentrations when using the fractional contributions
                            where (subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant).eq.0) temp_subgrid=0

                            if (save_netcdf_file_flag) then
                                !write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                !write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_nodata(temp_subgrid,size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3),NODATA_value)
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))


                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                ! write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif

                        endif

                        !Special case for exhaust as this must be given as a fraction for both PM2.5 and PM10.
                        !if ((pollutant_loop_index(i_pollutant).eq.pm10_index.or.pollutant_loop_index(i_pollutant).eq.pm25_index).and.i_source.eq.traffic_index) then
                        !Write all exhaust pollutants except the pmex
                        if (i_source.eq.traffic_index.and.(.not.pollutant_index.eq.all_totals_nc_index.and..not.pollutant_index.eq.aaqd_totals_nc_index.and..not.pollutant_index.eq.gp_totals_nc_index.and..not.pollutant_index.eq.op_totals_nc_index.or.pollutant_loop_index(i_pollutant).eq.nox_index)) then

                            !If PM then exhaust output is exhaust
                            if (pollutant_loop_index(i_pollutant).eq.pm10_index.or.pollutant_loop_index(i_pollutant).eq.pm25_index) then
                                if (source_domain_loop == 1) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,pollutant_loop_back_index(pmex_index))
                                else if (source_domain_loop == 2) then
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,pollutant_loop_back_index(pmex_index))
                                else if (source_domain_loop == 5) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,pollutant_loop_back_index(pmex_index))-subgrid(:,:,:,emep_local_subgrid_index,traffic_exhaust_nc_index,i_pollutant)+subgrid(:,:,:,emep_additional_local_subgrid_index,traffic_exhaust_nc_index,i_pollutant)
                                else !source_domain_loop == 6
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,pollutant_loop_back_index(pmex_index))+subgrid_EMEP_semilocal_from_in_region(:,:,:,traffic_exhaust_nc_index,i_pollutant)
                                end if
                            else
                                if (source_domain_loop == 1) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)
                                else if (source_domain_loop == 2) then
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)
                                else if (source_domain_loop == 5) then
                                    temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)-subgrid(:,:,:,emep_local_subgrid_index,i_source,i_pollutant)+subgrid(:,:,:,emep_additional_local_subgrid_index,i_source,i_pollutant)
                                else !source_domain_loop == 6
                                    temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)+subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,i_pollutant)
                                end if
                            endif
                            var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//'_exhaust'//trim(filename_append)
                            if (.not. save_netcdf_fraction_as_contribution_flag) then
                                temp_subgrid=temp_subgrid/subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)*100.
                            end if
                            where (subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant).eq.0) temp_subgrid=0

                            if (save_netcdf_file_flag) then
                                !write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                !write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif

                        endif
                    endif

                    !Special case for salt and sand
                    if (pollutant_loop_index(i_pollutant).eq.pm25_sand_index.or.pollutant_loop_index(i_pollutant).eq.pm10_sand_index &
                        .or.pollutant_loop_index(i_pollutant).eq.pm25_salt_index.or.pollutant_loop_index(i_pollutant).eq.pm10_salt_index) then
                        !Save the nonexhaust sand and salt
                        ! NB: sand and salt do not have EMEP contributions, so total is the same as local
                        if (i_source.eq.traffic_index) then
                            var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//trim(filename_append)
                            if (source_domain_loop == 1 .or. source_domain_loop == 5) then
                                temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)
                            else  ! 2 or 6
                                temp_subgrid=subgrid_local_from_in_region(:,:,:,i_source,i_pollutant)
                            end if
                            if (.not. save_netcdf_fraction_as_contribution_flag) then
                                !temp_subgrid=subgrid(:,:,:,local_subgrid_index,i_source,i_pollutant)
                                if (pollutant_loop_index(i_pollutant).eq.pm10_sand_index.or.pollutant_loop_index(i_pollutant).eq.pm10_salt_index) then
                                    temp_subgrid=temp_subgrid/subgrid(:,:,:,total_subgrid_index,allsource_index,pollutant_loop_back_index(pm10_index))*100.
                                else
                                    temp_subgrid=temp_subgrid/subgrid(:,:,:,total_subgrid_index,allsource_index,pollutant_loop_back_index(pm25_index))*100.
                                endif

                            endif
                            where (subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant).eq.0) temp_subgrid=0

                            if (save_netcdf_file_flag) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif

                        endif
                    endif
                end do  ! i_source
            end do  ! i_pollutant
        end do  ! source_domain_loop

        ! Save the total nonlocal contributions from EMEP (not saving it for in-region domains!)
        ! 1 = nonlocal to the moving window
        ! 2 = not used, but kept for consistency with other source_domain_loop loops
        ! 3 = nonlocal to the additional domain
        do source_domain_loop = 1,3
            if (source_domain_loop == 2) cycle
            write(unit_logfile,'(A)') '-------------------------------'
            if (source_domain_loop == 1) write(unit_logfile,'(A)') 'Saving nonlocal contributions'
            if (source_domain_loop == 3) write(unit_logfile,'(A)') 'Saving additional nonlocal contributions'
            write(unit_logfile,'(A)') '-------------------------------'
            do i_pollutant = 1, n_pollutant_loop
                if (pollutant_loop_index(i_pollutant).ne.pmex_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_sand_index.and.pollutant_loop_index(i_pollutant).ne.pm10_sand_index &
                    .and.pollutant_loop_index(i_pollutant).ne.pm25_salt_index.and.pollutant_loop_index(i_pollutant).ne.pm10_salt_index) then
                    if (source_domain_loop == 1) then
                        if (.not. (save_source_contributions .or. save_emep_source_contributions .or. save_total_source_contributions)) cycle
                        i_file=emep_subgrid_nonlocal_file_index(allsource_index)
                        temp_subgrid=subgrid(:,:,:,emep_nonlocal_subgrid_index,allsource_index,i_pollutant)
                    else ! source_domain_loop == 3
                        if (.not. (EMEP_additional_grid_interpolation_size.gt.0 .and. (save_emep_additional_source_contributions .or. save_total_source_contributions))) cycle
                        i_file=emep_additional_subgrid_nonlocal_file_index(allsource_index)
                        temp_subgrid=subgrid(:,:,:,emep_additional_nonlocal_subgrid_index,allsource_index,i_pollutant)
                    end if
                    var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))
                    if (.not. (save_netcdf_fraction_as_contribution_flag)) then
                        temp_subgrid=temp_subgrid/subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)*100.
                    end if
                    where (subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant).eq.0) temp_subgrid=0

                    if (save_netcdf_file_flag) then
                        write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                        call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                    endif
                    if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                        write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                        call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str_rec,create_file_rec,valid_min &
                            ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,z_rec(allsource_index,1) &
                            ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                    endif
                end if
            enddo  ! i_pollutant
        end do  ! source_domain_loop


        ! Calculate and save NO2 and O3 source contributions
        if (save_no2_source_contributions .or. save_o3_source_contributions) then

            ! Calculate NO2 and O3 source contributions
            if (EMEP_additional_grid_interpolation_size.gt.0) then
                calculate_EMEP_additional_grid_flag=.true.
                call uEMEP_source_fraction_chemistry
            endif

            calculate_EMEP_additional_grid_flag=.false.
            call uEMEP_source_fraction_chemistry

            ! Save the contributions to file
            ! 1 = normal source contributions from within moving window
            ! 2 = source contributions cut down to region
            ! 3 = additional nonlocal
            ! 4 = semilocal contribution from in region, based on background
            ! 5 = not used, but kept to be consistent with other loops over "source_domain_loop"
            ! 6 = total contribution from in region = 2+4
            do source_domain_loop = 1,6
                
                ! Skip this round in the loop if this type of source contributions are not to be saved
                ! or set appropriate variable name extension
                if (source_domain_loop == 1) then
                    if (.not. save_source_contributions) cycle
                    filename_append=''
                else if (source_domain_loop == 2) then
                    if (.not. (trace_emissions_from_in_region .and. save_local_source_contributions_from_in_region)) cycle
                    filename_append='_from_in_region'
                else if (source_domain_loop == 3) then
                    if (.not. (EMEP_additional_grid_interpolation_size.gt.0 .and. (save_emep_additional_source_contributions .or. save_total_source_contributions))) cycle
                    filename_append=''
                else if (source_domain_loop == 4) then
                    if (.not. (trace_emissions_from_in_region .and. save_semilocal_source_contributions_from_in_region)) cycle
                    filename_append='_from_in_region'
                else if (source_domain_loop == 5) then
                    cycle
                else  ! source_domain_loop == 6
                    if (.not. (trace_emissions_from_in_region .and. save_total_source_contributions_from_in_region)) cycle
                    filename_append='_from_in_region'
                end if

                write(unit_logfile,'(A)') '-------------------------------'
                if (source_domain_loop == 1) write(unit_logfile,'(A)') 'Saving local NO2 and O3 contributions'
                if (source_domain_loop == 2) write(unit_logfile,'(A)') 'Saving local NO2 and O3 contributions from-in-region'
                if (source_domain_loop == 3) write(unit_logfile,'(A)') 'Saving additional nonlocal NO2 and O3 contributions'
                if (source_domain_loop == 4) write(unit_logfile,'(A)') 'Saving semilocal NO2 and O3 contributions from-in-region'
                if (source_domain_loop == 6) write(unit_logfile,'(A)') 'Saving total NO2 and O3 contributions from-in-region'
                write(unit_logfile,'(A)') '-------------------------------'

                ! Save NO2 in round 1 and O3 in round 2
                do no2_o3_loop = 1,2

                    if (no2_o3_loop == 1) then  ! no2
                        comp_index = no2_index
                        comp_nc_index = no2_nc_index
                        if (.not. save_no2_source_contributions) cycle

                        valid_min=0.

                        if (save_netcdf_fraction_as_contribution_flag) then
                            variable_type='float'
                            unit_str="ug/m3"
                        else
                            variable_type='byte'
                            unit_str="%"
                        endif

                    else  ! o3
                        comp_index = o3_index
                        comp_nc_index = o3_nc_index
                        if (.not. save_o3_source_contributions) cycle

                        if (save_netcdf_fraction_as_contribution_flag) then
                            variable_type='float'
                            unit_str="ug/m3"
                            valid_min=-1000.
                        else
                            variable_type='short'
                            unit_str="%"
                            valid_min=-100.
                        endif
                    end if

                    do i_source=1,n_source_index

                        if (i_source .eq. allsource_index .or. calculate_source(i_source) .or. calculate_emep_source(i_source)) then
                            !Do not save nonexhaust for exhaust gas emissions
                            if (i_source.eq.traffic_nonexhaust_index) then
                                cycle
                            end if

                            if (i_source == allsource_index) then
                                ! Non-local contributions
                                if (source_domain_loop == 1) then
                                    ! nonlocal to the normal domain
                                    i_file=emep_subgrid_nonlocal_file_index(i_source)
                                    temp_subgrid=comp_source_EMEP_subgrid(:,:,:,comp_index,i_source)
                                else if (source_domain_loop == 3) then
                                    ! nonlocal to the additional domain
                                    i_file=emep_additional_subgrid_nonlocal_file_index(i_source)
                                    temp_subgrid=comp_source_EMEP_additional_subgrid(:,:,:,comp_index,i_source)
                                else
                                    ! do not save nonlocal for the in-region versions
                                    cycle
                                end if
                            else
                                ! Sector source contributions
                                if (source_domain_loop == 1) then
                                    ! normal local domain
                                    if (calculate_source(i_source)) then
                                        i_file=subgrid_local_file_index(i_source)
                                    else ! calculate_emep_source(i_source)
                                        i_file=emep_subgrid_local_file_index(i_source)
                                    endif
                                    temp_subgrid=comp_source_subgrid(:,:,:,comp_index,i_source)
                                else if (source_domain_loop == 2) then
                                    ! in-region contributions within the normal local domain
                                    if (calculate_source(i_source)) then
                                        i_file=subgrid_local_file_index(i_source)
                                    else ! calculate_emep_source(i_source)
                                        i_file=emep_subgrid_local_file_index(i_source)
                                    endif
                                    temp_subgrid=comp_source_subgrid_from_in_region(:,:,:,comp_index,i_source)
                                else if (source_domain_loop == 3) then
                                    ! additional domain: save only non-local for additional domain
                                    cycle
                                else if (source_domain_loop == 4) then
                                    ! semilocal inregion contributions
                                    i_file=emep_subgrid_semilocal_file_index(i_source)
                                    temp_subgrid=comp_semilocal_source_subgrid_from_in_region(:,:,:,comp_index,i_source)
                                else ! source_domain_loop == 6
                                    ! sum of local inregion and semilocal inregion
                                    i_file=subgrid_sourcetotal_inregion_file_index(i_source)
                                    temp_subgrid=comp_semilocal_source_subgrid_from_in_region(:,:,:,comp_index,i_source)+comp_source_subgrid_from_in_region(:,:,:,comp_index,i_source)
                                end if
                            end if

                            ! Variable name
                            if (i_source.eq.traffic_index .and. comp_index.eq.no2_index) then
                                ! special case for no2 traffic: add exhaust to variable name
                                var_name_temp=trim(var_name_nc(conc_nc_index,comp_nc_index,allsource_nc_index))//'_'//trim(filename_grid(i_file))//'_exhaust'//trim(filename_append)
                            else
                                var_name_temp=trim(var_name_nc(conc_nc_index,comp_nc_index,allsource_nc_index))//'_'//trim(filename_grid(i_file))//trim(filename_append)
                            endif

                            if (.not. save_netcdf_fraction_as_contribution_flag) then
                                ! Transform to fraction
                                if (comp_index == o3_index) then
                                    ! For O3, only include nonlocal when saving as fractions, and normalize with EMEP concentration to always get 100%
                                    if (i_source == allsource_index) then
                                        temp_subgrid=temp_subgrid/comp_EMEP_subgrid(:,:,:,comp_index)
                                        temp_subgrid=min(temp_subgrid,1000.)
                                        temp_subgrid=max(temp_subgrid,-100.)
                                    else
                                        cycle
                                    end if
                                else
                                    ! for NO2, normalize with EMEP concentration for additional, and with uEMEP concentration in all other cases
                                    if (source_domain_loop == 3) then
                                        temp_subgrid=temp_subgrid/comp_EMEP_subgrid(:,:,:,comp_index)*100.
                                    else
                                        temp_subgrid=temp_subgrid/comp_subgrid(:,:,:,comp_index)*100.
                                    end if
                                end if
                            endif

                            ! Save to file
                            if (save_netcdf_file_flag) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif
                        endif
                    enddo !i_source

                end do !no2_o3_loop

            end do  !source_domain_loop

            valid_min=0.

        end if !(save_no2_source_contributions .or. save_o3_source_contributions)


        !Save the emissions interpolated to the target grid
        if (save_emissions) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving emissions'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            unit_str="g/s"
            do i_pollutant=1,n_pollutant_loop
                do i_source=1,n_source_index
                    if (calculate_source(i_source)) then

                        !Do not save emissions if the source is not traffic and the pollutant is road sand and salt or exhaust as these do not exist
                        if (i_source.ne.traffic_index.and.(pollutant_loop_index(i_pollutant).eq.pm25_sand_index.or.pollutant_loop_index(i_pollutant).eq.pm10_sand_index &
                            .or.pollutant_loop_index(i_pollutant).eq.pm25_salt_index.or.pollutant_loop_index(i_pollutant).eq.pm10_salt_index &
                            .or.pollutant_loop_index(i_pollutant).eq.pmex_index)) then
                            !Do nothing because these pollutants only exist for traffic
                        else

                            i_file=emission_file_index(i_source)
                            var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))

                            !Calculate the emissions in the target grid
                            temp_subgrid=0.
                            do j=1,subgrid_dim(y_dim_index)
                                do i=1,subgrid_dim(x_dim_index)


                                    ii=crossreference_target_to_emission_subgrid(i,j,x_dim_index,i_source)
                                    jj=crossreference_target_to_emission_subgrid(i,j,y_dim_index,i_source)

                                    temp_subgrid(i,j,:)=emission_subgrid(ii,jj,:,i_source,i_pollutant)
                                    !Subgrid emissions, if relevant, are in units ug/sec/subgrid. Convert to g/s. Acount for the difference in subgrid sizes here
                                    temp_subgrid(i,j,:)=1.0e-6*temp_subgrid(i,j,:)*(subgrid_delta(y_dim_index)*subgrid_delta(x_dim_index)) &
                                        /(emission_subgrid_delta(y_dim_index,i_source)*emission_subgrid_delta(x_dim_index,i_source))

                                enddo
                            enddo

                            if (save_netcdf_file_flag) then
                                write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif

                        endif

                    endif
                enddo
            enddo
        endif

        !Save population interpolated to the target grid
        if (save_population) then
            write(unit_logfile,'(A)') '-------------------------------'
            write(unit_logfile,'(A)') 'Saving population data'
            write(unit_logfile,'(A)') '-------------------------------'
            
            variable_type='float'
            unit_str="inhabitants/grid"

            i_file=population_file_index(population_index)
            var_name_temp=trim(filename_grid(i_file))

            !Calculate the population in the target grid
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)

                    ii=crossreference_target_to_population_subgrid(i,j,x_dim_index)
                    jj=crossreference_target_to_population_subgrid(i,j,y_dim_index)

                    temp_subgrid(i,j,:)=population_subgrid(ii,jj,population_index)
                    !Acount for the difference in subgrid sizes here
                    temp_subgrid(i,j,:)=temp_subgrid(i,j,:)*(subgrid_delta(y_dim_index)*subgrid_delta(x_dim_index)) &
                        /(population_subgrid_delta(y_dim_index)*population_subgrid_delta(x_dim_index))

                enddo
            enddo

            unit_str='inhabitants/grid'
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif

        endif

        ! Save region mask at target subgrid
        if (trace_emissions_from_in_region) then  !! .or. use_region_select_and_mask_flag ??
            write(unit_logfile,'(A)') '-------------------------------'
            write(unit_logfile,'(A)') 'Saving region mask'
            write(unit_logfile,'(A)') '-------------------------------'

            variable_type = 'short'
            var_name_temp = 'region_index'
            unit_str='1'  !!!!!!! what to write here???
            temp_subgrid = 0.
            do tt = 1,subgrid_dim(t_dim_index)
                ! add 0.1 to all values to avoid them being rounded down to 1 less than original value
                temp_subgrid(:,:,tt) = subgrid_region_index + 0.1
            end do
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                ! NB: hard-coding scale_factor=1 since it does not make sense to have scale_factor for an ID variable!
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str,create_file,valid_min,variable_type,1.)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                ! check if region is the same for the whole grid. If not, set to zero to avoid getting wrong ID when interpolating
                if (.not. (minval(subgrid_region_index) == maxval(subgrid_region_index))) then
                    write(unit_logfile,'(A)') 'Warning: Not the same region_index for the whole receptor grid. Setting to 0.'
                    temp_subgrid = 0.
                end if
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,1.)
            endif
        end if

        ! Save EMEP allsources
        ! Save only if pure EMEP contributions are saved
        if (save_emep_source_contributions .or. save_emep_additional_source_contributions) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving EMEP allsources'
            write(unit_logfile,'(a)')'--------------------------'
            do i_pollutant=1,n_emep_pollutant_loop
                !EMEP does not have all pollutants so only save to n_emep_pollutant_loop
                !Only save the allsource value here 'EMEP_allsources'
                variable_type='float'
                i_file=emep_subgrid_file_index(allsource_index)
                var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))
                unit_str='ug/m3'
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
            end do
        end if

        !Save the EMEP data interpolated to the subgrid. These are based on the gridded concentrations
        ! Loop over 5 different verions
        ! 1 = emep local contributions from the moving-window (downscaling domain)
        ! 2 = As 1, but only the part of the moving window that is within the receptor region
        ! 3 = emep additional local contributions, extending further using larger LF grid
        ! 4 = Semilocal contribution, i.e. outside moving-window but inside the receptor region
        ! 5 = total contribution: Same as 3 (but saved as a different name and under different flag)
        ! 6 = total contribution from_in_region = 2+4, only for sources not downscaled
        do source_domain_loop=1,6

            ! check if this version of source contributions should be saved
            if (source_domain_loop == 1) then
                if (.not. save_emep_source_contributions) cycle
                filename_append = ''
            else if (source_domain_loop == 2) then
                if (.not. (save_local_source_contributions_from_in_region .and. trace_emissions_from_in_region)) cycle
                filename_append = '_from_in_region'
            else if (source_domain_loop == 3) then
                if (.not. (save_emep_additional_source_contributions .and. EMEP_additional_grid_interpolation_size.gt.0)) cycle
                filename_append = ''
            else if (source_domain_loop == 4) then
                if (.not. (save_semilocal_source_contributions_from_in_region .and. trace_emissions_from_in_region)) cycle
                filename_append = '_from_in_region'
            else if (source_domain_loop == 5) then
                if (.not. (save_total_source_contributions .and. EMEP_additional_grid_interpolation_size.gt.0)) cycle
                filename_append = ''
            else !source_domain_loop == 6
                if (.not. (save_total_source_contributions_from_in_region .and. trace_emissions_from_in_region)) cycle
                filename_append = '_from_in_region'
            end if

            write(unit_logfile,'(a)')'--------------------------'
            if (source_domain_loop.eq.1) write(unit_logfile,'(a)')'Saving EMEP contributions'
            if (source_domain_loop.eq.2) write(unit_logfile,'(a)')'Saving EMEP contributions from in region'
            if (source_domain_loop.eq.3) write(unit_logfile,'(a)')'Saving EMEP additional contributions'
            if (source_domain_loop.eq.4) write(unit_logfile,'(a)')'Saving EMEP semilocal contributions from in region'
            if (source_domain_loop.eq.5) write(unit_logfile,'(a)')'Saving total contributions'
            if (source_domain_loop.eq.6) write(unit_logfile,'(a)')'Saving total contributions from in region'
            write(unit_logfile,'(a)')'--------------------------'

            if (save_netcdf_fraction_as_contribution_flag) then
                variable_type='float'
                unit_str="ug/m3"
            else
                variable_type='byte'
                unit_str="%"
            endif

            do i_pollutant=1,n_emep_pollutant_loop
                !EMEP does not have all pollutants so only save to n_emep_pollutant_loop
                do i_source=1,n_source_index
                    
                    if (calculate_source(i_source).or.save_EMEP_source(i_source).or.i_source.eq.allsource_index) then

                        !Do not save for the additional GNFR sources
                        if (i_source.eq.traffic_gasoline_nc_index .or. i_source.eq.traffic_diesel_nc_index .or. i_source.eq.traffic_gas_nc_index .or. i_source.eq.publicpower_point_nc_index .or. i_source.eq.publicpower_area_nc_index) then
                            cycle
                        end if

                        !Do not save nonexhaust for exhaust gas emissions
                        if (i_source.eq.traffic_nonexhaust_nc_index.and.(pollutant_loop_index(i_pollutant).eq.no2_index.or.pollutant_loop_index(i_pollutant).eq.nox_index.or.pollutant_loop_index(i_pollutant).eq.o3_index)) then
                            cycle
                        end if

                        if (source_domain_loop == 1) then
                            ! local EMEP contribution
                            i_file=emep_subgrid_local_file_index(i_source)
                            temp_subgrid = subgrid(:,:,:,emep_local_subgrid_index,i_source,i_pollutant)
                        else if (source_domain_loop == 2) then
                            ! local EMEP contribution cut down to in-region
                            i_file=emep_subgrid_local_file_index(i_source)
                            temp_subgrid = subgrid_EMEP_local_from_in_region(:,:,:,i_source,i_pollutant)
                        else if (source_domain_loop == 3) then
                            ! additional local EMEP contribution
                            i_file=emep_additional_subgrid_local_file_index(i_source)
                            temp_subgrid = subgrid(:,:,:,emep_additional_local_subgrid_index,i_source,i_pollutant)
                        else if (source_domain_loop == 4) then
                            ! semilocal EMEP contribution
                            i_file = emep_subgrid_semilocal_file_index(i_source)
                            temp_subgrid = subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,i_pollutant)
                        else   ! source_domain_loop is 5 or 6
                            ! Only save total for non-downscaled sources, since downscaled sources were saved earlier
                            if (i_source == allsource_index .or. calculate_source(i_source)) then
                                cycle
                            else if ((i_source == traffic_exhaust_nc_index .or. i_source == traffic_nonexhaust_nc_index) .and. calculate_source(traffic_index)) then
                                ! if traffic is downscaled, then total contributions to traffic exhaust and nonexhaust are saved earlier in the subroutine using downscaled local contributions
                                cycle
                            end if
                            if (source_domain_loop == 5) then
                                ! EMEP contribution from additional domain
                                i_file = subgrid_sourcetotal_file_index(i_source)
                                temp_subgrid = subgrid(:,:,:,emep_additional_local_subgrid_index,i_source,i_pollutant)
                            else !source_domain_loop == 6
                                ! sum of EMEP local and EMEP semilocal contribution
                                i_file = subgrid_sourcetotal_inregion_file_index(i_source)
                                temp_subgrid = subgrid_EMEP_local_from_in_region(:,:,:,i_source,i_pollutant) + subgrid_EMEP_semilocal_from_in_region(:,:,:,i_source,i_pollutant)
                            end if
                        end if

                        var_name_temp=trim(var_name_nc(conc_nc_index,pollutant_loop_index(i_pollutant),allsource_index))//'_'//trim(filename_grid(i_file))//trim(filename_append)
                        !if (i_source.eq.traffic_exhaust_nc_index) var_name_temp=var_name_temp//'_exhaust'
                        !if (i_source.eq.traffic_nonexhaust_nc_index) var_name_temp=var_name_temp//'_nonexhaust'
                        ! Transform to fraction
                        if (.not. save_netcdf_fraction_as_contribution_flag) then
                            temp_subgrid=temp_subgrid/subgrid(:,:,:,emep_subgrid_index,allsource_index,i_pollutant)*100.
                            temp_subgrid=min(temp_subgrid,1000.)
                            temp_subgrid=max(temp_subgrid,-1000.)
                        endif
                        if (save_netcdf_file_flag) then
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                            call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                        endif
                        if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                            call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str_rec,create_file_rec,valid_min &
                                ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,z_rec(allsource_index,1) &
                                ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                        endif

                    endif
                enddo ! i_source
            enddo ! i_pollutant

        enddo !source_domain_loop

        !Save the other interpolated EMEP compounds used for nox chemistry as well. These are based on the surface comp values
        if (save_for_chemistry) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving for offline chemistry'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            do i_pollutant=1,n_emep_pollutant_loop
                if (pollutant_loop_index(i_pollutant).eq.nox_index) then
                    do i_loop=1,n_pollutant_compound_loop(i_pollutant)

                        i_comp=pollutant_compound_loop_index(i_pollutant,i_loop)

                        !Somo35 may be included here. Do not include it if it is.
                        if (i_comp.ne.somo35_nc_index) then

                            i_file=emep_subgrid_file_index(allsource_index)
                            var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_forchemistry'//'_'//trim(filename_grid(i_file))
                            temp_subgrid=comp_EMEP_subgrid(:,:,:,i_comp)
                            unit_str="ug/m3"
                            if (save_netcdf_file_flag) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(comp_EMEP_subgrid(:,:,:,i_comp))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif

                        endif
                    enddo
                endif
            enddo
        endif

        !Save weighted travel time
        if (save_for_chemistry) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving weighted travel time'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            var_name_temp='Weighted_travel_time'
            unit_str='seconds'
            i_pollutant=pollutant_loop_back_index(nox_nc_index) !Only save the travel time for nox. Though this may be expanded for other compounds if necessary, like ammonia
            temp_subgrid=traveltime_subgrid(:,:,:,3,i_pollutant)
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                !write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                !write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(traveltime_subgrid(:,:,:,1,i_pollutant))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
        endif

        if (save_for_chemistry) then
            variable_type='float'
            i_file=subgrid_J_file_index
            i_meteo=J_subgrid_index
            unit_str="1/s"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                !write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                !write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            !Save temperature if not also saving the meteo data.
            if (.not.save_other_meteo) then
                variable_type='float'
                i_file=subgrid_t2m_file_index
                i_meteo=t2m_subgrid_index
                unit_str="Celcius"
                !The same for all--------------------
                var_name_temp=trim(filename_grid(i_file))
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)-273.13
                    enddo
                enddo
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    !write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    !write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
            endif
        endif

        if (save_emep_species) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving EMEP species'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            do i_sp=1,n_species_loop_index
                ii_sp=species_loop_index(i_sp)
                if (i_sp.ne.sp_pm_index) then
                    !Do not include the total
                    do i_loop=1,n_pmxx_sp_index
                        if (i_loop.eq.pm25_sp_index.or.i_loop.eq.pm10_sp_index) then

                            if (i_loop.eq.pm25_sp_index) i_pollutant=pollutant_loop_back_index(pm25_index)
                            if (i_loop.eq.pm10_sp_index) i_pollutant=pollutant_loop_back_index(pm10_index)

                            if (save_netcdf_fraction_as_contribution_flag) then
                                variable_type='float'
                                var_name_temp=trim(species_name_nc(i_loop,ii_sp))//'_nonlocal_contribution'
                                unit_str="ug/m3"
                                temp_subgrid=species_EMEP_subgrid(:,:,:,i_loop,i_sp)
                            else
                                variable_type='byte'
                                var_name_temp=trim(species_name_nc(i_loop,i_sp))//'_nonlocal_fraction'
                                unit_str="%"
                                temp_subgrid=species_EMEP_subgrid(:,:,:,i_loop,i_sp)/subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)*100.
                            endif

                            if (save_netcdf_file_flag) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid(:,:,:))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                            endif
                            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid(:,:,:))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                    ,z_rec(allsource_index,1) &
                                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                            endif
                        endif
                    enddo
                endif
            enddo
        endif

        !Only save sea salt in this way when emep species are not saved in the general way
        if (save_seasalt.and..not.save_emep_species) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving EMEP sea salt'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            i_sp=n_species_loop_index
            ii_sp=species_loop_index(i_sp)
            do i_loop=1,n_pmxx_sp_index
                if (i_loop.eq.pm25_sp_index.or.i_loop.eq.pm10_sp_index) then

                    if (i_loop.eq.pm25_sp_index) i_pollutant=pollutant_loop_back_index(pm25_index)
                    if (i_loop.eq.pm10_sp_index) i_pollutant=pollutant_loop_back_index(pm10_index)
                    if (save_netcdf_fraction_as_contribution_flag) then
                        variable_type='float'
                        var_name_temp=trim(species_name_nc(i_loop,ii_sp))//'_nonlocal_contribution'
                        unit_str="ug/m3"
                        temp_subgrid=species_EMEP_subgrid(:,:,:,i_loop,i_sp)
                    else
                        variable_type='byte'
                        var_name_temp=trim(species_name_nc(i_loop,ii_sp))//'_nonlocal_fraction'
                        unit_str="%"
                        temp_subgrid=species_EMEP_subgrid(:,:,:,i_loop,i_sp)/subgrid(:,:,:,total_subgrid_index,allsource_index,i_pollutant)*100.
                    endif

                    if (save_netcdf_file_flag) then
                        write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid(:,:,:))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                        call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                    endif
                    if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                        write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid(:,:,:))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                        call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str_rec,create_file_rec,valid_min &
                            ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,z_rec(allsource_index,1) &
                            ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                    endif
                endif
            enddo

        endif

        !Save the original EMEP compounds
        if (save_emep_original) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving original EMEP'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            do i_pollutant=1,n_emep_pollutant_loop
                if (pollutant_loop_index(i_pollutant).ne.pmex_index) then
                    do i_loop=1,n_pollutant_compound_loop(i_pollutant)

                        i_comp=pollutant_compound_loop_index(i_pollutant,i_loop)

                        i_file=emep_subgrid_file_index(allsource_index)
                        var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_original_EMEP_concentration'
                        unit_str="ug/m3"
                        if (i_comp.eq.somo35_index) unit_str="ppb�d"
                        if (save_netcdf_file_flag) then
                            !write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),sum(orig_EMEP_subgrid(:,:,:,i_comp))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(orig_EMEP_subgrid(:,:,:,i_comp),use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                            call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,orig_EMEP_subgrid(:,:,:,i_comp),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                        endif
                        if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                            ! write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(orig_EMEP_subgrid(:,:,:,i_comp))/size(subgrid,1)/size(subgrid,2)/size(subgrid,3)
                            write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(orig_EMEP_subgrid(:,:,:,i_comp),use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                            call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,orig_EMEP_subgrid(:,:,:,i_comp),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str_rec,create_file_rec,valid_min &
                                ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,z_rec(allsource_index,1) &
                                ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                        endif
                    enddo
                endif
            enddo
        endif

        !Save AQI
        if (save_aqi) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving AQI'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='byte'
            scale_factor=0.1
            !Hard coded AQI limits
            aqi_pollutant_index(1)=no2_index;aqi_pollutant_index(2)=pm10_index;aqi_pollutant_index(3)=pm25_index;aqi_pollutant_index(4)=o3_index
            aqi_limits_temp(:,2:4)=aqi_hourly_limits(:,1:3)
            aqi_limits_temp(:,1)=0.
            aqi_limits_temp(:,5)=2.*aqi_hourly_limits(:,3)
            aqi_subgrid=0.

            do t=1,subgrid_dim(t_dim_index)
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        max_aqi=0.
                        do l=1,n_save_aqi_pollutant_index !pollutant (no2,pm10,pm2.5,o3)
                            do k=1,4 !level
                                if (comp_subgrid(i,j,t,aqi_pollutant_index(l)).ge.aqi_limits_temp(aqi_pollutant_index(l),k).and.comp_subgrid(i,j,t,aqi_pollutant_index(l)).lt.aqi_limits_temp(aqi_pollutant_index(l),k+1)) then
                                    aqi_subgrid(i,j,t,aqi_pollutant_index(l))=k+(comp_subgrid(i,j,t,aqi_pollutant_index(l))-aqi_limits_temp(aqi_pollutant_index(l),k))/(aqi_limits_temp(aqi_pollutant_index(l),k+1)-aqi_limits_temp(aqi_pollutant_index(l),k))
                                endif
                            enddo
                            aqi_subgrid(i,j,t,aqi_pollutant_index(l))=min(aqi_subgrid(i,j,t,aqi_pollutant_index(l)),4.99)
                            if (aqi_subgrid(i,j,t,aqi_pollutant_index(l)).gt.max_aqi) then
                                max_aqi=aqi_subgrid(i,j,t,aqi_pollutant_index(l))
                                aqi_responsible_pollutant_index(i,j,t)=l
                            endif
                            !write(*,*)  i,j,t,aqi_responsible_pollutant_index(i,j,t),aqi_subgrid(i,j,t,aqi_pollutant_index(l)),max_aqi
                        enddo
                    enddo
                enddo
            enddo

            do l=1,n_save_aqi_pollutant_index
                !write(unit_logfile,*)  'MAX AQI in time and space from '//trim(pollutant_file_str(aqi_pollutant_index(l)))//' = ',maxval(aqi_subgrid(:,:,:,aqi_pollutant_index(l)))
            enddo

            var_name_temp='AQI'
            unit_str='1'
            !Take the maximum of the pollutants
            temp_subgrid=maxval(aqi_subgrid(:,:,:,aqi_pollutant_index(1:n_save_aqi_pollutant_index)),4)/scale_factor
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid*scale_factor)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                    ,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif


            do l=1,n_save_aqi_pollutant_index

                i_comp=aqi_pollutant_index(l)
                var_name_temp='AQI_'//trim(var_name_nc(conc_nc_index,i_comp,allsource_index))
                unit_str='1'
                !Take the maximum of the pollutants
                temp_subgrid=aqi_subgrid(:,:,:,i_comp)/scale_factor
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a)')'Writing netcdf variable: '//trim(var_name_temp)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid*scale_factor)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

            enddo

            !Reset scale_factor
            scale_factor=1.

        endif

        !Save deposition
        if (save_deposition.and.calculate_deposition_flag) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving deposition'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            do i_pollutant=1,n_pollutant_loop
                do i_source=1,n_source_index
                    if (calculate_source(i_source)) then

                        i_comp=pollutant_loop_index(i_pollutant)

                        var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_local_dry_deposition_'//source_file_str(i_source)
                        unit_str="mg/m2/hr"
                        !Save in mg per hour, conversition from ug/m2/s
                        temp_subgrid=subgrid(:,:,:,drydepo_local_subgrid_index,i_source,i_pollutant)/1000.*3600.
                        if (save_netcdf_file_flag) then
                            write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                            call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                        endif
                        if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                            write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                            call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str_rec,create_file_rec,valid_min &
                                ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,z_rec(allsource_index,1) &
                                ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                        endif

                        var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_local_wet_deposition_'//source_file_str(i_source)
                        unit_str="mg/m2/hr"
                        !Save in mg per hour, conversition from ug/m2/s
                        temp_subgrid=subgrid(:,:,:,wetdepo_local_subgrid_index,i_source,i_pollutant)/1000.*3600.
                        if (save_netcdf_file_flag) then
                            write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                            call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                        endif
                        if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                            write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                            call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                                ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                                ,unit_str,title_str_rec,create_file_rec,valid_min &
                                ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                                ,z_rec(allsource_index,1) &
                                ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                        endif

                    endif
                enddo

                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_nonlocal_dry_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=subgrid(:,:,:,drydepo_nonlocal_subgrid_index,allsource_index,i_pollutant)/1000.*3600.
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_total_dry_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=(subgrid(:,:,:,drydepo_nonlocal_subgrid_index,allsource_index,i_pollutant)+subgrid(:,:,:,drydepo_local_subgrid_index,allsource_index,i_pollutant))/1000.*3600.
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_nonlocal_wet_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=subgrid(:,:,:,wetdepo_nonlocal_subgrid_index,allsource_index,i_pollutant)/1000.*3600.
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_total_wet_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=(subgrid(:,:,:,wetdepo_nonlocal_subgrid_index,allsource_index,i_pollutant)+subgrid(:,:,:,wetdepo_local_subgrid_index,allsource_index,i_pollutant))/1000.*3600.
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

                !Deposition velocity
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_deposition_subgrid(i,j,x_dim_index);jj=crossreference_target_to_deposition_subgrid(i,j,y_dim_index)
                        temp_subgrid(i,j,:)=deposition_subgrid(ii,jj,:,vd_index,i_pollutant)
                    enddo
                enddo
                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_deposition_velocity_'//source_file_str(allsource_index)
                unit_str="cm/s"
                !Save in cm/s, conversition from ug/m2/s
                temp_subgrid=temp_subgrid*100.
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

                !Landuse category
                if (read_landuse_flag) then
                    temp_subgrid=0.
                    do j=1,subgrid_dim(y_dim_index)
                        do i=1,subgrid_dim(x_dim_index)
                            ii=crossreference_target_to_deposition_subgrid(i,j,x_dim_index);jj=crossreference_target_to_deposition_subgrid(i,j,y_dim_index)
                            temp_subgrid(i,j,:)=landuse_subgrid(ii,jj,grid_index)
                        enddo
                    enddo
                    var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_EMEP_landuse_category'
                    unit_str="1"
                    variable_type='byte'
                    if (save_netcdf_file_flag) then
                        write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                        call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                    endif
                    if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                        write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                        call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                            ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                            ,unit_str,title_str_rec,create_file_rec,valid_min &
                            ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                            ,z_rec(allsource_index,1) &
                            ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                    endif
                endif

                variable_type='float'

                !Save the original emep wet and dry deposition
                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_original_EMEP_wet_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=orig_emep_deposition_subgrid(:,:,:,wetdepo_index,i_pollutant)
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
                var_name_temp=trim(var_name_nc(conc_nc_index,i_comp,allsource_index))//'_original_EMEP_dry_deposition_'//source_file_str(allsource_index)
                unit_str="mg/m2/hr"
                !Save in mg per hour, conversition from ug/m2/s
                temp_subgrid=orig_emep_deposition_subgrid(:,:,:,drydepo_index,i_pollutant)
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,es12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid,x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp &
                        ,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif

            enddo

        endif


        !Save the meteo interpolated to the target grid
        valid_min=-1.e24

        if (save_wind_vectors) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving wind vectors'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            if (wind_vectors_10m_available) then
                i_file=subgrid_u10_file_index
                i_meteo=u10_subgrid_index
            else
                i_file=subgrid_ugrid_file_index
                i_meteo=ugrid_subgrid_index
            endif

            unit_str="m/s"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            if (wind_vectors_10m_available) then
                i_file=subgrid_v10_file_index
                i_meteo=v10_subgrid_index
            else
                i_file=subgrid_vgrid_file_index
                i_meteo=vgrid_subgrid_index
            endif
            unit_str="m/s"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp),mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------
        endif


        if (save_other_meteo) then
            write(unit_logfile,'(a)')'--------------------------'
            write(unit_logfile,'(a)')'Saving other meteo data'
            write(unit_logfile,'(a)')'--------------------------'
            variable_type='float'
            i_file=subgrid_hmix_file_index
            i_meteo=hmix_subgrid_index
            unit_str="m"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            variable_type='float'
            i_file=subgrid_t2m_file_index
            i_meteo=t2m_subgrid_index
            unit_str="Celcius"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)-273.13
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.2)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.2)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            variable_type='float'
            i_file=subgrid_hmix_file_index
            i_meteo=precip_subgrid_index
            unit_str="mm/hr"
            !The same for all--------------------
            var_name_temp='precipitation'!trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------


            i_file=subgrid_kz_file_index
            i_meteo=kz_subgrid_index
            unit_str="m2/s"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            if (hourly_calculations) then

                if (wind_vectors_10m_available) then
                    i_file=subgrid_FF10_file_index
                    i_meteo=FF10_subgrid_index
                else
                    i_file=subgrid_FFgrid_file_index
                    i_meteo=FFgrid_subgrid_index
                endif
                unit_str="m/s"
                !The same for all--------------------
                var_name_temp=trim(filename_grid(i_file))
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                    enddo
                enddo
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
                !The same for all--------------------

                if (wind_vectors_10m_available) then
                    i_file=subgrid_DD10_file_index
                else
                    i_file=subgrid_DDgrid_file_index
                endif
                !i_meteo=DDgrid_subgrid_index !i_meteo not used in this conversion to wind direction
                unit_str="degrees"
                do tt=1,subgrid_dim(t_dim_index)
                    do jj=1,integral_subgrid_dim(y_dim_index)
                        do ii=1,integral_subgrid_dim(x_dim_index)
                            if (wind_vectors_10m_available) then
                                temp_integral_subgrid(ii,jj,tt)=DIRECTION(meteo_subgrid(ii,jj,tt,u10_subgrid_index),meteo_subgrid(ii,jj,tt,v10_subgrid_index))
                            else
                                temp_integral_subgrid(ii,jj,tt)=DIRECTION(meteo_subgrid(ii,jj,tt,ugrid_subgrid_index),meteo_subgrid(ii,jj,tt,vgrid_subgrid_index))
                            endif

                        enddo
                    enddo
                enddo
                !This is not converted to real degrees, but subgrid degrees, so is wrong but close

                !The same for all--------------------
                var_name_temp=trim(filename_grid(i_file))
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);
                        !In this case this is not the same for all
                        temp_subgrid(i,j,:)=temp_integral_subgrid(ii,jj,:)
                    enddo
                enddo
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
                !The same for all--------------------

            endif

            i_file=subgrid_invL_file_index
            i_meteo=invL_subgrid_index
            unit_str="1/m"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.4)')'Writing netcdf variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.4)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            i_file=subgrid_logz0_file_index
            i_meteo=logz0_subgrid_index
            unit_str="log(m)"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            i_file=subgrid_ustar_file_index
            i_meteo=ustar_subgrid_index
            unit_str="m/s"
            !The same for all--------------------
            var_name_temp=trim(filename_grid(i_file))
            temp_subgrid=0.
            do j=1,subgrid_dim(y_dim_index)
                do i=1,subgrid_dim(x_dim_index)
                    ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                enddo
            enddo
            if (save_netcdf_file_flag) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
            endif
            if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                    ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                    ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                    ,z_rec(allsource_index,1) &
                    ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
            endif
            !The same for all--------------------

            if (annual_calculations) then

                i_file=subgrid_invFFgrid_file_index
                i_meteo=inv_FFgrid_subgrid_index
                unit_str="s/m"
                !The same for all--------------------
                var_name_temp=trim(filename_grid(i_file))
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                    enddo
                enddo
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
                !The same for all--------------------

                i_file=subgrid_invFF10_file_index
                i_meteo=inv_FF10_subgrid_index
                unit_str="s/m"
                !The same for all--------------------
                var_name_temp=trim(filename_grid(i_file))
                temp_subgrid=0.
                do j=1,subgrid_dim(y_dim_index)
                    do i=1,subgrid_dim(x_dim_index)
                        ii=crossreference_target_to_integral_subgrid(i,j,x_dim_index);jj=crossreference_target_to_integral_subgrid(i,j,y_dim_index);temp_subgrid(i,j,:)=meteo_subgrid(ii,jj,:,i_meteo)
                    enddo
                enddo
                if (save_netcdf_file_flag) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf variable: '//trim(var_name_temp), mean_mask(temp_subgrid,use_subgrid(:,:,allsource_index),size(temp_subgrid,1),size(temp_subgrid,2),size(temp_subgrid,3))
                    call uEMEP_save_netcdf_file(unit_logfile,temp_name,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str,create_file,valid_min,variable_type,scale_factor)
                endif
                if (save_netcdf_receptor_flag.and.n_valid_receptor.ne.0) then
                    write(unit_logfile,'(a,f12.3)')'Writing netcdf receptor variable: '//trim(var_name_temp),sum(temp_subgrid)/size(temp_subgrid,1)/size(temp_subgrid,2)/size(temp_subgrid,3)
                    call uEMEP_save_netcdf_receptor_file(unit_logfile,temp_name_rec,subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index) &
                        ,temp_subgrid(:,:,:),x_subgrid,y_subgrid,lon_subgrid,lat_subgrid,var_name_temp,unit_str,title_str_rec,create_file_rec,valid_min &
                        ,x_receptor(valid_receptor_index(1:n_valid_receptor)),y_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,lon_receptor(valid_receptor_index(1:n_valid_receptor)),lat_receptor(valid_receptor_index(1:n_valid_receptor)) &
                        ,z_rec(allsource_index,1) &
                        ,name_receptor(valid_receptor_index(1:n_valid_receptor),1),n_valid_receptor,variable_type,scale_factor)
                endif
                !The same for all--------------------

            endif

        endif

        !Deallocate the averaging arrays with each receptor grid
        if (use_multiple_receptor_grids_flag.and..not.use_single_time_loop_flag.and.save_netcdf_average_flag) then
            if (allocated(val_array_av)) deallocate(val_array_av)
            if (allocated(time_seconds_output_av)) deallocate(time_seconds_output_av)
            counter_av=0
        endif


    end subroutine uEMEP_save_netcdf_control