subroutine uEMEP_read_roadlink_data_ascii use uEMEP_definitions implicit none integer i,j integer unit_in integer rl_length_short logical :: exists logical nxtdat_flag real sub_nodes_x(5000),sub_nodes_y(5000),sub_nodes_lon(5000),sub_nodes_lat(5000) integer temp_id,n_subnodes,temp_road_type,temp_nlanes real temp_adt,temp_hdv,temp_speed,temp_width integer counter real size_major(3) integer n_loop,loop_step real x_grid_min,x_grid_max,y_grid_min,y_grid_max integer counter_major,counter_sub logical :: show_diagnostics=.false. real diagnostic_val(10) integer temp_category,temp_structure_type,temp_region_id,temp_surface_id,temp_route_id real temp_length,temp_tunnel_length logical :: road_data_in_latlon=.false. real temp_real real, allocatable :: inputdata_rl_temp(:,:) integer, allocatable :: inputdata_int_rl_temp(:,:) real, allocatable :: inputdata_rl_multi(:,:) integer, allocatable :: inputdata_int_rl_multi(:,:) real, allocatable :: inputdata_rl_multi_new(:,:) integer, allocatable :: inputdata_int_rl_multi_new(:,:) integer n_multi_roadlinks_new,n_multi_roadlinks integer m integer n_road_link_file_loop logical :: first_road_link_file_read=.false. integer temp_int integer :: io write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Reading road link data ascii (uEMEP_read_roadlink_data_ascii)' write(unit_logfile,'(A)') '================================================================' min_adt=10. min_link_size=1. if (num_multiple_roadlink_files.eq.0) then n_road_link_file_loop=1 else n_road_link_file_loop=num_multiple_roadlink_files endif !Start the file loop here !--------------------------- n_multi_roadlinks_new=0 n_multi_roadlinks=0 first_road_link_file_read=.false. do m=1,n_road_link_file_loop if (num_multiple_roadlink_files.gt.0) then !pathname_rl(1)=pathname_mrl(m) filename_rl(1)=filename_mrl(m) write(*,*) m,n_road_link_file_loop,trim(filename_rl(1)) endif pathfilename_rl(1)=trim(pathname_rl(1))//trim(filename_rl(1)) !write(*,*) pathname_rl(2),filename_rl(2),pathfilename_rl(2) !Test existence of the road link filename (2). If does not exist then stop inquire(file=trim(pathfilename_rl(1)),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: Road link file ascii does not exist: ', trim(pathfilename_rl(1)) stop endif !Check filename for the string 'latlon' that specifies if the data is in lat lon coordinates if (index(filename_rl(1),'latlon').gt.0) then road_data_in_latlon=.true. write(unit_logfile,'(A,A)') ' Reading road data positions as lat lon: ', trim(pathfilename_rl(1)) endif !Open the file for reading to test the available links in the region if (reduce_roadlink_region_flag) then unit_in=20 open(unit_in,file=pathfilename_rl(1),access='sequential',status='old',readonly) write(unit_logfile,'(a)') ' Opening road link file(ascii) '//trim(pathfilename_rl(1)) rewind(unit_in) call nxtdat(unit_in,nxtdat_flag) !read the header to find out how many links there are !read(unit_in,'(a)',ERR=20) temp_str if (no_header_roadlink_data_flag) then write(unit_logfile,'(a)') ' Reading road link file(ascii) without header: '//trim(pathfilename_rl(1)) i=0 do read(unit_in,*,iostat=io) temp_real,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_x(1) !Read x nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_y(1) !Read y nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if if (temp_real.gt.0) then i=i+1 n_roadlinks=n_roadlinks+n_subnodes-1 endif if (int(temp_real).ne.temp_real) then write(unit_logfile,'(a,i,f)') ' Problem with record at point with ID: ',i,temp_real stop endif enddo n_roadlinks_major=i !n_roadlinks=0 rewind(unit_in) call nxtdat(unit_in, nxtdat_flag) else read(unit_in,*,ERR=20) n_roadlinks_major,n_roadlinks endif if (n_roadlinks.eq.0) then write(unit_logfile,'(a)') ' Reading road link file(ascii) with header but without subnode counts: '//trim(pathfilename_rl(1)) rewind(unit_in) call nxtdat(unit_in,nxtdat_flag) read(unit_in,*,ERR=20) n_roadlinks_major,n_roadlinks i=0 do read(unit_in,*,iostat=io) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes !Read attributes up to n_subnodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_x(1) !Read x nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_y(1) !Read y nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if if (temp_id.gt.0) then i=i+1 n_roadlinks=n_roadlinks+n_subnodes-1 endif enddo !Have commented out this in the cases where the number of links are not as written. Can happen with OSM data !n_roadlinks_major=i rewind(unit_in) call nxtdat(unit_in, nxtdat_flag) read(unit_in,*,ERR=20) temp_int,temp_int endif write(unit_logfile,'(a,i)') ' Number of major road links in header = ', n_roadlinks_major write(unit_logfile,'(a,i)') ' Number of sub road links in header= ', n_roadlinks !Allocate the arrays after reading in the number of roads !allocate (inputdata_rl_temp(n_roadlinks,num_var_rl)) !allocate (inputdata_int_rl_temp(n_roadlinks,num_int_rl)) if (allocated(valid_link_flag)) deallocate (valid_link_flag) allocate (valid_link_flag(n_roadlinks_major)) valid_link_flag=.false. !Initialise !inputdata_rl_temp=0. !inputdata_int_rl_temp=0 counter_major=0 counter_sub=0 !Read the data to find roads in the tile !x_grid_min=emission_subgrid_min(x_dim_index,traffic_index) !x_grid_max=emission_subgrid_min(x_dim_index,traffic_index)+(emission_subgrid_dim(x_dim_index,traffic_index)+1)*emission_subgrid_delta(x_dim_index,traffic_index) !y_grid_min=emission_subgrid_min(y_dim_index,traffic_index) !y_grid_max=emission_subgrid_min(y_dim_index,traffic_index)+(emission_subgrid_dim(y_dim_index,traffic_index)+1)*emission_subgrid_delta(y_dim_index,traffic_index) x_grid_min=init_emission_subgrid_min(x_dim_index,traffic_index) x_grid_max=init_emission_subgrid_min(x_dim_index,traffic_index)+(init_emission_subgrid_dim(x_dim_index,traffic_index)+1)*init_emission_subgrid_delta(x_dim_index,traffic_index) y_grid_min=init_emission_subgrid_min(y_dim_index,traffic_index) y_grid_max=init_emission_subgrid_min(y_dim_index,traffic_index)+(init_emission_subgrid_dim(y_dim_index,traffic_index)+1)*init_emission_subgrid_delta(y_dim_index,traffic_index) !Under the special case where multiple roads are read from different countries then us the intital grid to determine which to keep !This is not implemented yet but needs to select the region when reading multiple files !if !init_subgrid_min(x_dim_index)=subgrid_min(x_dim_index) !init_subgrid_min(y_dim_index)=subgrid_min(y_dim_index) !init_subgrid_max(x_dim_index)=subgrid_max(x_dim_index) !init_subgrid_max(y_dim_index)=subgrid_max(y_dim_index) !endif !rewind(unit_in) !call nxtdat_modern(unit_in,nxtdat_flag) !read(unit_in,*,ERR=20) temp_int,temp_int do i=1,n_roadlinks_major !ID ADT HDV ROAD_TYPE SPEED N_SUBLINKS read(unit_in,*,ERR=20) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes !read(unit_in,*,ERR=20) !temp_id !write(*,*) temp_id,temp_adt,n_subnodes if (road_data_in_latlon) then !Order is reversed in these files read(unit_in,*) sub_nodes_y(1) read(unit_in,*) sub_nodes_x(1) if (projection_type.eq.UTM_projection_index) then !write(*,*) i,sub_nodes_x(1),sub_nodes_y(1) call ll2utm(1,utm_zone,sub_nodes_y(1),sub_nodes_x(1),sub_nodes_y(1),sub_nodes_x(1)) !write(*,*) i,sub_nodes_x(1),sub_nodes_y(1) elseif (projection_type.eq.LTM_projection_index) then call ll2ltm(1,ltm_lon0,sub_nodes_y(1),sub_nodes_x(1),sub_nodes_y(1),sub_nodes_x(1)) elseif (projection_type.eq.LCC_projection_index) then elseif (projection_type.eq.PS_projection_index) then elseif (projection_type.eq.LAEA_projection_index) then call LL2LAEA(sub_nodes_x(1),sub_nodes_y(1),sub_nodes_x(1),sub_nodes_y(1),projection_attributes) endif else read(unit_in,*) sub_nodes_x(1) read(unit_in,*) sub_nodes_y(1) endif !Convert to EMEP coordinates from UTM to lambert. No choices here. Special case but should be changed if (save_emissions_for_EMEP(traffic_index)) then call PROJ2LL(sub_nodes_x(1),sub_nodes_y(1),sub_nodes_lon(1),sub_nodes_lat(1),projection_attributes,projection_type) !call UTM2LL(utm_zone,sub_nodes_y(1),sub_nodes_x(1),sub_nodes_lat(1),sub_nodes_lon(1)) if (projection_type.eq.LCC_projection_index) then call lb2lambert2_uEMEP(sub_nodes_x(1),sub_nodes_y(1),sub_nodes_lon(1),sub_nodes_lat(1),EMEP_projection_attributes) elseif (projection_type.eq.PS_projection_index) then call LL2PS_spherical(sub_nodes_x(1),sub_nodes_y(1),sub_nodes_lon(1),sub_nodes_lat(1),EMEP_projection_attributes) else !Remains as lat lon endif !write(*,*) sub_nodes_x(1),sub_nodes_y(1),sub_nodes_lon(1),sub_nodes_lat(1) endif !write(*,*) sub_nodes_x(1),sub_nodes_y(1) !write(*,*) x_grid_min,x_grid_max,y_grid_min,y_grid_max !Test position within emission region if (sub_nodes_x(1).ge.x_grid_min.and.sub_nodes_x(1).le.x_grid_max & .and.sub_nodes_y(1).ge.y_grid_min.and.sub_nodes_y(1).le.y_grid_max) then counter_major=counter_major+1 counter_sub=counter_sub+n_subnodes-1 valid_link_flag(i)=.true. else valid_link_flag(i)=.false. endif enddo close(unit_in,status='keep') write(unit_logfile,'(a,i,i)') ' Number of major and sub road links within the region = ', counter_major,counter_sub endif !Open the file for reading unit_in=20 open(unit_in,file=pathfilename_rl(1),access='sequential',status='old',readonly) write(unit_logfile,'(a)') ' Opening road link file(ascii) '//trim(pathfilename_rl(1)) rewind(unit_in) call nxtdat(unit_in, nxtdat_flag) !read the header to find out how many links there are !read(unit_in,'(a)',ERR=20) temp_str if (no_header_roadlink_data_flag) then write(unit_logfile,'(a)') ' Reading road link file(ascii) without header: '//trim(pathfilename_rl(1)) i=0 do read(unit_in,*,iostat=io) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes !Read attributes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_x(1) !Read x nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_y(1) !Read y nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if if (temp_id.gt.0) then i=i+1 endif enddo n_roadlinks_major=i n_roadlinks=n_roadlinks+n_subnodes-1 rewind(unit_in) call nxtdat(unit_in,nxtdat_flag) else read(unit_in,*,ERR=20) n_roadlinks_major,n_roadlinks endif if (n_roadlinks.eq.0.and..not.reduce_roadlink_region_flag) then write(unit_logfile,'(a)') ' Reading road link file(ascii) with header but without subnode counts: '//trim(pathfilename_rl(1)) rewind(unit_in) call nxtdat(unit_in,nxtdat_flag) read(unit_in,*,ERR=20) n_roadlinks_major,n_roadlinks i=0 do read(unit_in,*,iostat=io) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes !Read attributes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_x(1) !Read x nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_y(1) !Read y nodes if (io < 0) then exit else if (io > 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if !write(*,*) temp_id if (temp_id.gt.0) then i=i+1 n_roadlinks=n_roadlinks+n_subnodes-1 endif enddo !n_roadlinks_major=i rewind(unit_in) call nxtdat(unit_in, nxtdat_flag) endif write(unit_logfile,'(a,i)') ' Number of major road links= ', n_roadlinks_major write(unit_logfile,'(a,i)') ' Number of sub road links= ', n_roadlinks !Allocate the arrays after reading in the number of roads if (reduce_roadlink_region_flag) then n_roadlinks=counter_sub n_roadlinks_major_selected=counter_major else !Set all road links to valid when not regionally selecting. THis will fail because already allocated if (allocated(valid_link_flag)) deallocate (valid_link_flag) allocate (valid_link_flag(n_roadlinks_major)) valid_link_flag=.true. endif if (allocated(inputdata_rl_temp)) deallocate (inputdata_rl_temp) if (allocated(inputdata_int_rl_temp)) deallocate (inputdata_int_rl_temp) allocate (inputdata_rl_temp(n_roadlinks,num_var_rl)) allocate (inputdata_int_rl_temp(n_roadlinks,num_int_rl)) !Initialise inputdata_rl_temp=0. inputdata_int_rl_temp=0 counter=0 counter_major=0 !rewind(unit_in) !call nxtdat_modern(unit_in,nxtdat_flag) !read(unit_in,*,ERR=20) temp_int,temp_int !Read the data do i=1,n_roadlinks_major ! if (.not.valid_link_flag(i)) then ! if (read_OSM_roadlink_data_flag) then !ID ADT HDV ROAD_TYPE SPEED N_SUBLINKS ! if (.not.eof(unit_in)) read(unit_in,*,ERR=20) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes ! else !ID ADT HDV ROAD_ACTIVITY_TYPE SPEED ROAD_WIDTH N_LANES N_SUBNODES ROAD_CATEGORY ROAD_LENGTH ROAD_STRUCTURE_TYPE REGION_ID ROAD_SURFACE_ID TUNNEL_LENGTH ROUTE_ID ! if (.not.eof(unit_in)) read(unit_in,*,ERR=20) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes & ! ,temp_category,temp_length,temp_structure_type,temp_region_id,temp_surface_id,temp_tunnel_length,temp_route_id ! endif !write(*,*) temp_id,temp_adt,n_subnodes ! if (.not.eof(unit_in)) read(unit_in,*,ERR=20) temp_val ! if (.not.eof(unit_in)) read(unit_in,*,ERR=20) temp_val ! else if (read_OSM_roadlink_data_flag) then !ID ADT HDV ROAD_TYPE SPEED N_SUBLINKS read(unit_in,*,iostat=io) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes if (io /= 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if temp_category=0;temp_length=0;temp_structure_type=0;temp_region_id=0;temp_surface_id=0;temp_tunnel_length=0;temp_route_id=0 else !ID ADT HDV ROAD_ACTIVITY_TYPE SPEED ROAD_WIDTH N_LANES N_SUBNODES ROAD_CATEGORY ROAD_LENGTH ROAD_STRUCTURE_TYPE REGION_ID ROAD_SURFACE_ID TUNNEL_LENGTH ROUTE_ID read(unit_in,*,iostat=io) temp_id,temp_adt,temp_hdv,temp_road_type,temp_speed,temp_width,temp_nlanes,n_subnodes & ,temp_category,temp_length,temp_structure_type,temp_region_id,temp_surface_id,temp_tunnel_length,temp_route_id if (io /= 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if endif !write(*,*) i,temp_id,temp_adt,n_subnodes read(unit_in,*,iostat=io) sub_nodes_x(1:n_subnodes) if (io /= 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if read(unit_in,*,iostat=io) sub_nodes_y(1:n_subnodes) if (io /= 0) then write(unit_logfile,"(2a)") "ERROR reading road link file: ", trim(pathfilename_rl(1)) stop 1 end if !write(*,*) sub_nodes_x(1:n_subnodes),sub_nodes_y(1:n_subnodes) !put in the road link data !Test size of major link. If less than min_link_size then treat it as a single link from the start sub_node to the end sub_node if (valid_link_flag(i)) then !Do not do this with latlon input data as it does not work if (.not.road_data_in_latlon) then size_major(1)=maxval(sub_nodes_x(1:n_subnodes))-minval(sub_nodes_x(1:n_subnodes)) size_major(2)=maxval(sub_nodes_y(1:n_subnodes))-minval(sub_nodes_y(1:n_subnodes)) size_major(3)=sqrt(size_major(1)**2+size_major(2)**2) if (size_major(3).gt.min_link_size) then n_loop=n_subnodes-1 loop_step=1 else n_loop=1 loop_step=n_subnodes-1 !write(*,*) size_major(3) endif else n_loop=n_subnodes-1 loop_step=1 endif counter_major=counter_major+1 !if (n_subnodes.gt.5000) then ! write(*,*) 'Greater than 5000 ',i !endif !Place the data in the road links if (temp_adt.ge.min_adt) then do j=1,n_loop counter=counter+1 if (counter.le.n_roadlinks) then inputdata_int_rl_temp(counter,major_index_rl_index)=counter_major inputdata_int_rl_temp(counter,id_rl_index)=temp_id inputdata_rl_temp(counter,adt_rl_index)=temp_adt**osm_adt_power_scale inputdata_rl_temp(counter,hdv_rl_index)=temp_hdv inputdata_int_rl_temp(counter,roadtype_rl_index)=temp_road_type inputdata_rl_temp(counter,speed_rl_index)=temp_speed inputdata_rl_temp(counter,width_rl_index)=temp_width inputdata_int_rl_temp(counter,nlanes_rl_index)=temp_nlanes inputdata_rl_temp(counter,x1_rl_index)=sub_nodes_x(j) inputdata_rl_temp(counter,x2_rl_index)=sub_nodes_x(j+loop_step) inputdata_rl_temp(counter,y1_rl_index)=sub_nodes_y(j) inputdata_rl_temp(counter,y2_rl_index)=sub_nodes_y(j+loop_step) !write(*,*) inputdata_int_rl(counter,id_rl_index),inputdata_rl(counter,x1_rl_index),inputdata_rl(counter,y2_rl_index) inputdata_rl_temp(counter,tunnel_length_rl_index)=temp_tunnel_length endif enddo !write(*,*) counter, inputdata_int_rl_temp(counter,major_index_rl_index) endif endif enddo write(unit_logfile,'(a,i)') ' Number of counted road links= ', counter write(unit_logfile,'(a,i)') ' Number of road links allocated = ', n_roadlinks n_roadlinks=min(counter,n_roadlinks) close(unit_in,status='keep') !Allocate the arrays after reading in the number of roads if (allocated(inputdata_rl)) deallocate (inputdata_rl) if (allocated(inputdata_int_rl)) deallocate (inputdata_int_rl) allocate (inputdata_rl(n_roadlinks,num_var_rl)) allocate (inputdata_int_rl(n_roadlinks,num_int_rl)) inputdata_rl=inputdata_rl_temp(1:n_roadlinks,:) inputdata_int_rl=inputdata_int_rl_temp(1:n_roadlinks,:) if (road_data_in_latlon) then !Order is reversed in these files so they are reversed in the projection calls as well do i=1,n_roadlinks if (projection_type.eq.UTM_projection_index) then call ll2utm(1,utm_zone,inputdata_rl_temp(i,x1_rl_index),inputdata_rl_temp(i,y1_rl_index),inputdata_rl(i,y1_rl_index),inputdata_rl(i,x1_rl_index)) call ll2utm(1,utm_zone,inputdata_rl_temp(i,x2_rl_index),inputdata_rl_temp(i,y2_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl(i,x2_rl_index)) elseif (projection_type.eq.LTM_projection_index) then call ll2ltm(1,ltm_lon0,inputdata_rl_temp(i,x1_rl_index),inputdata_rl_temp(i,y1_rl_index),inputdata_rl(i,y1_rl_index),inputdata_rl(i,x1_rl_index)) call ll2ltm(1,ltm_lon0,inputdata_rl_temp(i,x2_rl_index),inputdata_rl_temp(i,y2_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl(i,x2_rl_index)) elseif (projection_type.eq.LCC_projection_index) then elseif (projection_type.eq.PS_projection_index) then elseif (projection_type.eq.LAEA_projection_index) then call LL2LAEA(inputdata_rl(i,x1_rl_index),inputdata_rl(i,y1_rl_index),inputdata_rl_temp(i,y1_rl_index),inputdata_rl_temp(i,x1_rl_index),projection_attributes) call LL2LAEA(inputdata_rl(i,x2_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl_temp(i,y2_rl_index),inputdata_rl_temp(i,x2_rl_index),projection_attributes) endif enddo endif if (num_multiple_roadlink_files.gt.1) then if (n_roadlinks.gt.0) then n_multi_roadlinks_new=n_multi_roadlinks+n_roadlinks !For the first roadloop if (.not.first_road_link_file_read) then !Allocate the multi road array for the first time n_multi_roadlinks_new=n_roadlinks n_multi_roadlinks=n_multi_roadlinks_new if (.not.allocated(inputdata_rl_multi)) allocate (inputdata_rl_multi(n_multi_roadlinks,num_var_rl)) if (.not.allocated(inputdata_int_rl_multi)) allocate (inputdata_int_rl_multi(n_multi_roadlinks,num_int_rl)) inputdata_rl_multi(1:n_multi_roadlinks,:)=inputdata_rl(1:n_multi_roadlinks,:) inputdata_int_rl_multi(1:n_multi_roadlinks,:)=inputdata_int_rl(1:n_roadlinks,:) first_road_link_file_read=.true. else !Allocate the new multi array with all roads so far allocate (inputdata_rl_multi_new(n_multi_roadlinks_new,num_var_rl)) allocate (inputdata_int_rl_multi_new(n_multi_roadlinks_new,num_int_rl)) !Place the current multi roads in the new multiroads !NOTE: inputdata_rl_multi is not allocated until later, debugging error but do not use this anyway! inputdata_rl_multi_new(1:n_multi_roadlinks,:)=inputdata_rl_multi(1:n_multi_roadlinks,:) inputdata_int_rl_multi_new(1:n_multi_roadlinks,:)=inputdata_int_rl_multi(1:n_roadlinks,:) !Place the last read road links in the new multiroads inputdata_rl_multi_new(n_multi_roadlinks+1:n_multi_roadlinks_new,:)=inputdata_rl(1:n_roadlinks,:) inputdata_int_rl_multi_new(n_multi_roadlinks+1:n_multi_roadlinks_new,:)=inputdata_int_rl(1:n_roadlinks,:) !Deallocate the old multi road arrays if (allocated(inputdata_rl_multi)) deallocate(inputdata_rl_multi) if (allocated(inputdata_int_rl_multi)) deallocate(inputdata_int_rl_multi) n_multi_roadlinks=n_multi_roadlinks_new !Allocate the multi road array allocate (inputdata_rl_multi(n_multi_roadlinks,num_var_rl)) allocate (inputdata_int_rl_multi(n_multi_roadlinks,num_int_rl)) !Put the new data in the old one inputdata_rl_multi=inputdata_rl_multi_new inputdata_int_rl_multi=inputdata_int_rl_multi_new !Deallocate the new multi road arrays if (allocated(inputdata_rl_multi_new)) deallocate(inputdata_rl_multi_new) if (allocated(inputdata_int_rl_multi_new)) deallocate(inputdata_int_rl_multi_new) endif endif !Deallocate the other arrays as well if (allocated(inputdata_rl)) deallocate(inputdata_rl) if (allocated(inputdata_int_rl)) deallocate(inputdata_int_rl) if (allocated(inputdata_rl_temp)) deallocate (inputdata_rl_temp) if (allocated(inputdata_int_rl_temp)) deallocate (inputdata_int_rl_temp) if (allocated(valid_link_flag)) deallocate (valid_link_flag) write(unit_logfile,'(a,i)') ' Number of accumulated multi-road links used = ', n_multi_roadlinks_new endif !End the multiloop here enddo if (num_multiple_roadlink_files.gt.1) then allocate (inputdata_rl(n_multi_roadlinks,num_var_rl)) allocate (inputdata_int_rl(n_multi_roadlinks,num_int_rl)) inputdata_rl=inputdata_rl_multi inputdata_int_rl=inputdata_int_rl_multi if (allocated(inputdata_rl_multi)) deallocate(inputdata_rl_multi) if (allocated(inputdata_int_rl_multi)) deallocate(inputdata_int_rl_multi) n_roadlinks=n_multi_roadlinks_new write(unit_logfile,'(a,i)') ' Number of accumulated multi-road links used = ', n_multi_roadlinks endif if (allocated(inputdata_rl_temp)) deallocate (inputdata_rl_temp) if (allocated(inputdata_int_rl_temp)) deallocate (inputdata_int_rl_temp) !No speed in the files currently. Set all to 50 km/hr. Temporary !inputdata_rl(:,speed_rl_index)=50. !inputdata_rl(:,width_rl_index)=10. !inputdata_int_rl(:,nlanes_rl_index)=2 !Set the road type, normal or tunnel (tunnel or jet). When a tunnel then there is no retention, always dry !do i=1,n_roadlinks ! if (inputdata_int_rl(i,roadtype_rl_index).eq.5.or.inputdata_int_rl(i,roadtype_rl_index).eq.6) then ! inputdata_int_rl(i,roadtype_rl_index)=tunnel_roadtype ! else ! inputdata_int_rl(i,roadtype_rl_index)=normal_roadtype ! endif !enddo !Calculate some additional values inputdata_rl(:,x0_rl_index)=(inputdata_rl(:,x1_rl_index)+inputdata_rl(:,x2_rl_index))/2. inputdata_rl(:,y0_rl_index)=(inputdata_rl(:,y1_rl_index)+inputdata_rl(:,y2_rl_index))/2. inputdata_rl(:,length_rl_index)=sqrt((inputdata_rl(:,x1_rl_index)-inputdata_rl(:,x2_rl_index))**2+(inputdata_rl(:,y1_rl_index)-inputdata_rl(:,y2_rl_index))**2) !Calculate road orientation and check for range overflows for length as well !inputdata_rl(angle_rl_index,:)=180./3.14159*acos((inputdata_rl(y2_rl_index,:)-inputdata_rl(y1_rl_index,:))/inputdata_rl(length_rl_index,:)) do i=1,n_roadlinks call PROJ2LL(inputdata_rl(i,x1_rl_index),inputdata_rl(i,y1_rl_index),inputdata_rl(i,lon1_rl_index),inputdata_rl(i,lat1_rl_index),projection_attributes,projection_type) call PROJ2LL(inputdata_rl(i,x2_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl(i,lon2_rl_index),inputdata_rl(i,lat2_rl_index),projection_attributes,projection_type) call PROJ2LL(inputdata_rl(i,x0_rl_index),inputdata_rl(i,y0_rl_index),inputdata_rl(i,lon0_rl_index),inputdata_rl(i,lat0_rl_index),projection_attributes,projection_type) !call UTM2LL(utm_zone,inputdata_rl(i,y1_rl_index),inputdata_rl(i,x1_rl_index),inputdata_rl(i,lat1_rl_index),inputdata_rl(i,lon1_rl_index)) !call UTM2LL(utm_zone,inputdata_rl(i,y2_rl_index),inputdata_rl(i,x2_rl_index),inputdata_rl(i,lat2_rl_index),inputdata_rl(i,lon2_rl_index)) !call UTM2LL(utm_zone,inputdata_rl(i,y0_rl_index),inputdata_rl(i,x0_rl_index),inputdata_rl(i,lat0_rl_index),inputdata_rl(i,lon0_rl_index)) enddo !Check lengths rl_length_short=0 do i=1,n_roadlinks if (inputdata_rl(i,length_rl_index).eq.0.0) then rl_length_short=rl_length_short+1 !write(unit_logfile,'(a,2i,f12.5)') ' WARNING: Zero link length, setting to 1 m ',i,inputdata_int_rl(i,id_rl_index),inputdata_rl(i,length_rl_index) inputdata_rl(i,length_rl_index)=1. endif enddo !write(*,*) 'Max length: ',maxval(inputdata_rl(:,length_rl_index)),inputdata_rl(maxloc(inputdata_rl(:,length_rl_index)),lat0_rl_index),inputdata_rl(maxloc(inputdata_rl(:,length_rl_index)),lon0_rl_index) write(unit_logfile,*) 'Number of road links with 0 length: ',rl_length_short,' Setting to 1 m' write(unit_logfile,*) 'Max road link x and y: ',maxval(inputdata_rl(:,x0_rl_index)),maxval(inputdata_rl(:,y0_rl_index)) write(unit_logfile,*) 'Min road link x and y: ',minval(inputdata_rl(:,x0_rl_index)),minval(inputdata_rl(:,y0_rl_index)) if (n_roadlinks.gt.0) then write(unit_logfile,'(a14,12a10)') ' LINK ','ID','X1','X2','Y1','Y2','WIDTH','LENGTH','ADT','LON','LAT','N_LANES','TYPE' i=1 write(unit_logfile,'(a14,i10,7f10.1,2f10.4,2i10)') ' First link = ',inputdata_int_rl(i,id_rl_index),inputdata_rl(i,x1_rl_index),inputdata_rl(i,x2_rl_index) & ,inputdata_rl(i,y1_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl(i,width_rl_index) & ,inputdata_rl(i,length_rl_index),inputdata_rl(i,adt_rl_index) & ,inputdata_rl(i,lon0_rl_index),inputdata_rl(i,lat0_rl_index) & ,inputdata_int_rl(i,nlanes_rl_index),inputdata_int_rl(i,roadtype_rl_index) i=n_roadlinks write(unit_logfile,'(a14,i10,7f10.1,2f10.4,2i10)') ' Last link = ',inputdata_int_rl(i,id_rl_index),inputdata_rl(i,x1_rl_index),inputdata_rl(i,x2_rl_index) & ,inputdata_rl(i,y1_rl_index),inputdata_rl(i,y2_rl_index),inputdata_rl(i,width_rl_index) & ,inputdata_rl(i,length_rl_index),inputdata_rl(i,adt_rl_index) & ,inputdata_rl(i,lon0_rl_index),inputdata_rl(i,lat0_rl_index) & ,inputdata_int_rl(i,nlanes_rl_index),inputdata_int_rl(i,roadtype_rl_index) else write(unit_logfile,'(a)') 'No road links available in this region' endif !Calculate the veh km totals if (show_diagnostics) then diagnostic_val=0. do i=1,n_roadlinks !Total kilometres diagnostic_val(1)=diagnostic_val(1)+inputdata_rl(i,length_rl_index)/1000. !Total veh kilometres diagnostic_val(2)=diagnostic_val(2)+inputdata_rl(i,length_rl_index)/1000.*inputdata_rl(i,adt_rl_index)*365. !Light veh kilometres diagnostic_val(3)=diagnostic_val(3)+(1.-inputdata_rl(i,hdv_rl_index)/100.)*inputdata_rl(i,length_rl_index)/1000.*inputdata_rl(i,adt_rl_index)*365. !Light veh kilometres diagnostic_val(4)=diagnostic_val(4)+(inputdata_rl(i,hdv_rl_index)/100.)*inputdata_rl(i,length_rl_index)/1000.*inputdata_rl(i,adt_rl_index)*365. enddo write(unit_logfile,'(a,es12.4)') 'Total km= ',diagnostic_val(1) write(unit_logfile,'(a,es12.4)') 'Total veh.km= ',diagnostic_val(2) write(unit_logfile,'(a,es12.4)') 'Total light veh.km= ',diagnostic_val(3) write(unit_logfile,'(a,es12.4)') 'Total heavy veh.km= ',diagnostic_val(4) endif return 20 write(unit_logfile,'(2A)') 'ERROR reading road link file: ',trim(pathfilename_rl(2)) stop end subroutine uEMEP_read_roadlink_data_ascii