subroutine uEMEP_grid_roads use uEMEP_definitions implicit none integer i,j integer ro real, allocatable :: f_subgrid(:) real, allocatable :: adt_temp(:,:) real, allocatable :: adt_car_temp(:) real, allocatable :: adt_truck_temp(:) real x_subgrid_in(2),y_subgrid_in(2) real x_line_in(2),y_line_in(2),lat_line_in(2),lon_line_in(2) integer i_traffic_index(2),j_traffic_index(2) integer i_start,i_end,j_start,j_end integer source_index,t integer i_roadlink_emission_compound(n_pollutant_loop) integer tt,ttt integer major_ro integer t_start_temp,t_end_temp real tunnel_ratio real sigma0_temp integer i_pollutant write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Gridding road link proxy data (uEMEP_grid_roads)' write(unit_logfile,'(A)') '================================================================' allocate (f_subgrid(n_roadlinks)) allocate (adt_temp(n_roadlinks,n_pollutant_loop)) allocate (adt_car_temp(n_roadlinks)) allocate (adt_truck_temp(n_roadlinks)) !allocate (traffic_emission_subgrid(subgrid_dim(1),subgrid_dim(2)),n_emission_subgrid_index) source_index=traffic_index t=1 proxy_emission_subgrid(:,:,source_index,:)=0. if (use_traffic_for_sigma0_flag) then emission_properties_subgrid(:,:,emission_sigy00_index,source_index)=0. emission_properties_subgrid(:,:,emission_sigz00_index,source_index)=0. endif !Possible to split the traffic source into different subsources at this point if necessary, e.g. light and heavy traffic !Here we weight the adt by the emission ratio and give an emission factor valid for cars adt_car_temp=inputdata_rl(1:n_roadlinks,adt_rl_index)*(1.-inputdata_rl(1:n_roadlinks,hdv_rl_index)/100.) adt_truck_temp=inputdata_rl(1:n_roadlinks,adt_rl_index)*inputdata_rl(1:n_roadlinks,hdv_rl_index)/100. do i_pollutant=1,n_pollutant_loop adt_temp(:,i_pollutant)=adt_car_temp+adt_truck_temp*ratio_truck_car_emission(pollutant_loop_index(i_pollutant)) enddo !Calculate the pseudo traffic emissions in each grid write(unit_logfile,*)'Gridding traffic emission proxy data' !Convert from uEMEP to NORTRIP if (use_NORTRIP_emission_data) then !This links the order of the NORTRIP output to the pollutants !Order is: pm10,pm2.5,pmex,nox,pm10_sand/salt,pm2.5_sand/salt,pm10_salt,pm2.5_salt do i_pollutant=1,n_pollutant_loop if (pollutant_loop_index(i_pollutant).eq.pm10_nc_index) then i_roadlink_emission_compound(i_pollutant)=1 elseif (pollutant_loop_index(i_pollutant).eq.pm25_nc_index) then i_roadlink_emission_compound(i_pollutant)=2 elseif (pollutant_loop_index(i_pollutant).eq.pmex_nc_index) then i_roadlink_emission_compound(i_pollutant)=3 elseif (pollutant_loop_index(i_pollutant).eq.nox_nc_index) then i_roadlink_emission_compound(i_pollutant)=4 elseif (pollutant_loop_index(i_pollutant).eq.pm10_sand_nc_index.and.(pollutant_index.eq.all_sand_nc_index.or.pollutant_index.eq.all_sand_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=5 elseif (pollutant_loop_index(i_pollutant).eq.pm10_salt_nc_index.and.(pollutant_index.eq.all_sand_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=6 elseif (pollutant_loop_index(i_pollutant).eq.pm10_salt_nc_index.and.(pollutant_index.eq.all_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=5 elseif (pollutant_loop_index(i_pollutant).eq.pm25_sand_nc_index.and.(pollutant_index.eq.all_sand_nc_index)) then i_roadlink_emission_compound(i_pollutant)=6 elseif (pollutant_loop_index(i_pollutant).eq.pm25_sand_nc_index.and.(pollutant_index.eq.all_sand_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=7 elseif (pollutant_loop_index(i_pollutant).eq.pm25_salt_nc_index.and.(pollutant_index.eq.all_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=6 elseif (pollutant_loop_index(i_pollutant).eq.pm25_salt_nc_index.and.(pollutant_index.eq.all_sand_salt_nc_index)) then i_roadlink_emission_compound(i_pollutant)=8 else write(unit_logfile,'(a,2i)') 'STOPPING: No valid compound chosen for NORTRIP. Stopping uEMEP_grid_roads. Pollutant index=',pollutant_index,pollutant_loop_index(i_pollutant) !write(*,*) n_pollutant_loop !write(*,*) pollutant_loop_index(1:n_pollutant_loop) stop endif enddo emission_subgrid(:,:,:,traffic_index,:)=0. !Set the time to be used based on the time loop flag if (use_single_time_loop_flag) then t_start_temp=t_loop t_end_temp=t_loop else t_start_temp=1 t_end_temp=subgrid_dim(t_dim_index) endif endif do ro=1,n_roadlinks x_line_in=inputdata_rl(ro,x1_rl_index:x2_rl_index) y_line_in=inputdata_rl(ro,y1_rl_index:y2_rl_index) !Convert to EMEP coordinates from specified projection type to lambertCC or latlon. Not certain if the fraction is correctly calculated in lat lon coordinates but otherwise very complicated if (save_emissions_for_EMEP(traffic_index)) then do i=1,2 call PROJ2LL(x_line_in(i),y_line_in(i),lon_line_in(i),lat_line_in(i),projection_attributes,projection_type) !call UTM2LL(utm_zone,y_line_in(i),x_line_in(i),lat_line_in(i),lon_line_in(i)) if (EMEP_projection_type.eq.LL_projection_index) then x_line_in(i)=lon_line_in(i) y_line_in(i)=lat_line_in(i) elseif (EMEP_projection_type.eq.LCC_projection_index) then call lb2lambert2_uEMEP(x_line_in(i),y_line_in(i),lon_line_in(i),lat_line_in(i),EMEP_projection_attributes) elseif (EMEP_projection_type.eq.PS_projection_index) then call LL2PS_spherical(x_line_in(i),y_line_in(i),lon_line_in(i),lat_line_in(i),EMEP_projection_attributes) endif enddo !write(*,*) x_line_in(1),y_line_in(1),lon_line_in(1),lat_line_in(1) endif i_traffic_index=1+floor((x_line_in-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index)) j_traffic_index=1+floor((y_line_in-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index)) if ((i_traffic_index(1).ge.1.or.i_traffic_index(2).ge.1).and.(j_traffic_index(1).ge.1.or.j_traffic_index(2).ge.1).and. & (i_traffic_index(1).le.emission_subgrid_dim(x_dim_index,source_index).or.i_traffic_index(2).le.emission_subgrid_dim(x_dim_index,source_index)).and. & (j_traffic_index(1).le.emission_subgrid_dim(y_dim_index,source_index).or.j_traffic_index(2).le.emission_subgrid_dim(y_dim_index,source_index))) then !write(*,*) ro,i_traffic_index,j_traffic_index !Limit the loop if it is near the edge i_start=max(1,minval(i_traffic_index)) i_end=min(emission_subgrid_dim(x_dim_index,source_index),maxval(i_traffic_index)) j_start=max(1,minval(j_traffic_index)) j_end=min(emission_subgrid_dim(y_dim_index,source_index),maxval(j_traffic_index)) !if (i_end-i_start.gt.2.or.j_end-j_start.gt.2) write(*,*) ro,i_start,i_end,j_start,j_end do j=j_start,j_end do i=i_start,i_end x_subgrid_in(1)=x_emission_subgrid(i,j,source_index)-emission_subgrid_delta(x_dim_index,source_index)/2. x_subgrid_in(2)=x_emission_subgrid(i,j,source_index)+emission_subgrid_delta(x_dim_index,source_index)/2. y_subgrid_in(1)=y_emission_subgrid(i,j,source_index)-emission_subgrid_delta(y_dim_index,source_index)/2. y_subgrid_in(2)=y_emission_subgrid(i,j,source_index)+emission_subgrid_delta(y_dim_index,source_index)/2. f_subgrid(ro)=line_fraction_in_grid_func(x_subgrid_in,y_subgrid_in,x_line_in,y_line_in) !do subsource_index=1,n_subsource(source_index) proxy_emission_subgrid(i,j,source_index,:)=proxy_emission_subgrid(i,j,source_index,:) & +inputdata_rl(ro,length_rl_index)*f_subgrid(ro)*adt_temp(ro,:) !Put the temporally changing emissions straight into the emission subgrid !Will not be overwritten in uEMEP_convert_proxy_to_emissions if use_NORTRIP_emission_data=true if (use_NORTRIP_emission_data) then major_ro=inputdata_int_rl(ro,major_index_rl_index) do tt=t_start_temp,t_end_temp !Convert from g/km/hour to ug/s/subgrid if (use_single_time_loop_flag) then ttt=t_loop t=1 else t=tt ttt=tt endif !write(*,*) 'ro,major_ro,tt,t,ttt',ro,major_ro,tt,t,ttt !write(*,*) 'emission_grid',shape(emission_subgrid) !write(*,*) 'inputdata_rl_emissions',shape(inputdata_rl_emissions) !if (t_loop.eq.2) stop do i_pollutant=1,n_pollutant_loop !write(*,*) i_pollutant,major_ro,inputdata_rl(major_ro,tunnel_length_rl_index) if (use_tunnel_deposition_flag.and.inputdata_rl(ro,tunnel_length_rl_index).gt.0) then call tunnel_deposition_factor(pollutant_loop_index(i_pollutant),inputdata_rl(ro,tunnel_length_rl_index) & ,inputdata_rl(ro,ADT_rl_index)*inputdata_rl(ro,length_rl_index)/inputdata_rl(ro,tunnel_length_rl_index) & ,ventilation_factor,min_ADT_ventilation_factor,min_length_ventilation_factor,windspeed_tunnel,tunnel_ratio) else tunnel_ratio=1. endif !Turn off tunnel emissions if required if (.not.use_tunnel_emissions_flag.and.inputdata_rl(ro,tunnel_length_rl_index).gt.0) then tunnel_ratio=0 endif !Converts from g/km/hr (NORTRIP) to ug/sec (uEMEP) emission_subgrid(i,j,t,source_index,i_pollutant)=emission_subgrid(i,j,t,source_index,i_pollutant) & +inputdata_rl(ro,length_rl_index)*f_subgrid(ro)*inputdata_rl_emissions(major_ro,ttt,i_roadlink_emission_compound(i_pollutant)) & *1.e6/1.e3/3600.*tunnel_ratio enddo !write(*,*) i,j, emission_subgrid(i,j,t,source_index,pollutant_loop_back_index(pm10_nc_index)),emission_subgrid(i,j,t,source_index,pollutant_loop_back_index(pm25_nc_index)) enddo endif !Set the sigma values according to traffic speed and road width using proxy weighting if (use_traffic_for_sigma0_flag.and..not.save_emissions_for_EMEP(traffic_index)) then sigma0_temp=sigma0_traffic_func(inputdata_rl(ro,speed_rl_index)) if (inputdata_rl(ro,tunnel_length_rl_index).gt.50.) sigma0_temp=tunnel_sig_z_00 emission_properties_subgrid(i,j,emission_sigy00_index,source_index)=emission_properties_subgrid(i,j,emission_sigy00_index,source_index) & +inputdata_rl(ro,length_rl_index)*f_subgrid(ro)*adt_temp(ro,1)*sqrt((inputdata_rl(ro,width_rl_index)/2.)**2+sigma0_temp**2) emission_properties_subgrid(i,j,emission_sigz00_index,source_index)=emission_properties_subgrid(i,j,emission_sigz00_index,source_index) & +inputdata_rl(ro,length_rl_index)*f_subgrid(ro)*adt_temp(ro,1)*sigma0_temp endif !enddo !write(*,*) ro,i,j,f_subgrid(ro) !write(*,*) ro,f_subgrid(ro),traffic_emission_subgrid(i,j,x_emission_subgrid_index),traffic_emission_subgrid(i,j,y_emission_subgrid_index),x_line_in,y_line_in enddo enddo !write(*,*) 'Gridding traffic emission',ro,' of ',n_roadlinks endif !if (mod(ro,10000).eq.0) write(*,*) 'Gridding traffic emission',ro,' of ',n_roadlinks enddo !Set the road properties based on ADT weighting if (use_traffic_for_sigma0_flag) then emission_properties_subgrid(:,:,emission_sigy00_index,source_index)=emission_properties_subgrid(:,:,emission_sigy00_index,source_index)/proxy_emission_subgrid(:,:,source_index,1) emission_properties_subgrid(:,:,emission_sigz00_index,source_index)=emission_properties_subgrid(:,:,emission_sigz00_index,source_index)/proxy_emission_subgrid(:,:,source_index,1) endif deallocate (f_subgrid) deallocate (adt_temp) deallocate (adt_car_temp) deallocate (adt_truck_temp) !Deallocate road link arrays after gridding but not when the external time step is used !and not when the multiple receptor grids are used !and not when the auto subgridding is used !because gridding roads is called again if (use_single_time_loop_flag.or.use_multiple_receptor_grids_flag.or.use_emission_positions_for_auto_subgrid_flag(allsource_index)) then !Do not deallocate because they will be used again else if (allocated(inputdata_rl)) deallocate(inputdata_rl) if (allocated(inputdata_int_rl)) deallocate(inputdata_int_rl) if (allocated(inputdata_rl_emissions)) deallocate(inputdata_rl_emissions) endif end subroutine uEMEP_grid_roads