uEMEP_set_deposition_velocities Subroutine

    subroutine uEMEP_set_deposition_velocities

        use uEMEP_definitions

        implicit none

        integer i_pollutant
        integer i,j

        real xpos_min,xpos_max,ypos_min,ypos_max
        real xpos_area_min,xpos_area_max,ypos_area_min,ypos_area_max
        real xpos_limit,ypos_limit
        integer i_nc,j_nc
        integer ii,jj
        real weighting_nc(-1:1,-1:1)

        integer i_landuse

        write(unit_logfile,'(A)')'Placing EMEP deposition velocities on the deposition subgrid'

        !Fill in with a default value
        deposition_subgrid(:,:,:,:,1:n_pollutant_loop)=0.
        !write(*,*) size(deposition_subgrid)

        drydepo_vd_default(1:n_compound_index)=0.0 !Positive value in m/s
        drydepo_vd_default(nh3_index)=0.03 !Positive value in m/s
        wetdepo_scavanging_rate(1:n_compound_index)=0.0 !Positive value in m/s
        wetdepo_scavanging_rate(nh3_index)=0.5*1.e6/1000. !Dimensions of m-1 W/h

        do i_pollutant=1,n_pollutant_loop
            deposition_subgrid(:,:,:,vd_index,i_pollutant)=drydepo_vd_default(pollutant_loop_index(i_pollutant))
            !write(*,*) i_source,i_pollutant,sum(deposition_subgrid(:,:,:,vd_index,i_pollutant)),drydepo_vd_default(i_pollutant)
        enddo
        !stop

        !Place the EMEP deposition velocities in the deposition subgrid, area weighted interpolation
        deposition_subgrid(:,:,:,vd_index,:)=0.

        xpos_limit=dgrid_nc(lon_nc_index)/2.
        ypos_limit=dgrid_nc(lat_nc_index)/2.

        do j=1,deposition_subgrid_dim(y_dim_index)
            do i=1,deposition_subgrid_dim(x_dim_index)

                !Cross reference EMEP grid with limits
                i_nc=crossreference_deposition_to_emep_subgrid(i,j,x_dim_index)
                j_nc=crossreference_deposition_to_emep_subgrid(i,j,y_dim_index)

                i_nc=min(max(2,i_nc),dim_length_nc(x_dim_nc_index)-1)
                j_nc=min(max(2,j_nc),dim_length_nc(y_dim_nc_index)-1)

                !write(*,*) i,j,i_nc,j_nc
                xpos_area_max=xproj_deposition_subgrid(i,j)+xpos_limit
                xpos_area_min=xproj_deposition_subgrid(i,j)-xpos_limit
                ypos_area_max=yproj_deposition_subgrid(i,j)+ypos_limit
                ypos_area_min=yproj_deposition_subgrid(i,j)-ypos_limit

                do jj=j_nc-1,j_nc+1
                    do ii=i_nc-1,i_nc+1

                        xpos_min=max(xpos_area_min,var1d_nc(ii,lon_nc_index)-dgrid_nc(lon_nc_index)/2.)
                        xpos_max=min(xpos_area_max,var1d_nc(ii,lon_nc_index)+dgrid_nc(lon_nc_index)/2.)
                        ypos_min=max(ypos_area_min,var1d_nc(jj,lat_nc_index)-dgrid_nc(lat_nc_index)/2.)
                        ypos_max=min(ypos_area_max,var1d_nc(jj,lat_nc_index)+dgrid_nc(lat_nc_index)/2.)

                        !Determine the area intersection of the EMEP grid and an EMEP grid size centred on the deposition subgrid
                        if (xpos_max.gt.xpos_min.and.ypos_max.gt.ypos_min) then
                            weighting_nc(ii-i_nc,jj-j_nc)=(ypos_max-ypos_min)*(xpos_max-xpos_min)/dgrid_nc(lon_nc_index)/dgrid_nc(lat_nc_index)
                        else
                            weighting_nc(ii-i_nc,jj-j_nc)=0.
                        endif


                        if (read_landuse_flag) then
                            do i_landuse=1,n_landuse_index-1
                                deposition_subgrid(i,j,:,vd_index,:)=deposition_subgrid(i,j,:,vd_index,:)+depo_var3d_nc(ii,jj,:,i_landuse,:)*landuse_subgrid(i,j,i_landuse)*weighting_nc(ii-i_nc,jj-j_nc)
                            enddo
                        else
                            deposition_subgrid(i,j,:,vd_index,:)=deposition_subgrid(i,j,:,vd_index,:)+depo_var3d_nc(ii,jj,:,grid_index,:)*weighting_nc(ii-i_nc,jj-j_nc)
                        endif


                    enddo
                enddo

                !Fill in zeros with the nearest EMEP value. 0 can appear because of non overlapping EMEP and landuse grids, I think
                where (deposition_subgrid(i,j,:,vd_index,:).eq.0) deposition_subgrid(i,j,:,vd_index,:)=depo_var3d_nc(i_nc,j_nc,:,grid_index,:)
                !write(*,'(2i,23es12.2,f)') i,j,deposition_subgrid(i,j,:,vd_index,:),depo_var3d_nc(i_nc,j_nc,:,grid_index,:),sum(depo_var3d_nc(ii,jj,:,1:n_landuse_index-1,:)),sum(landuse_subgrid(i,j,1:n_landuse_index-1))
            enddo
        enddo

    end subroutine uEMEP_set_deposition_velocities