subroutine uEMEP_set_region_tile_grids !This routine creates the tile files for regions/municipalities !It is based on the variable tile make 'uEMEP_set_tile_grids' so has a lot of unused code in it !It basically reads the masking file, in SSB coordinates, and then finds the extent of this to create municipality tile limits and various resolutions use uEMEP_definitions implicit none integer n_tiles_population_classes,n_tiles_traffic_classes,n_tiles_shipping_classes parameter (n_tiles_population_classes=5,n_tiles_traffic_classes=5,n_tiles_shipping_classes=1) real limit_val_tile_population(n_tiles_population_classes) data limit_val_tile_population /0.,100.,1000.,10000.,100000./ integer :: num_tiles_with_population(n_tiles_population_classes)=0 real limit_val_tile_traffic(n_tiles_traffic_classes) data limit_val_tile_traffic /0.,1000.,10000.,100000.,1000000./ integer :: num_tiles_with_traffic(n_tiles_traffic_classes)=0 real limit_val_tile_shipping(n_tiles_shipping_classes) data limit_val_tile_shipping /0./ integer :: num_tiles_with_shipping(n_tiles_shipping_classes)=0 real tile_subgrid_delta(n_dim_index) real tile_subgrid_min(n_dim_index) real tile_subgrid_max(n_dim_index) integer tile_subgrid_dim(n_dim_index) real, allocatable :: tile_subgrid(:,:,:) real, allocatable :: x_tile_subgrid(:,:) real, allocatable :: y_tile_subgrid(:,:) real, allocatable :: lon_tile_subgrid(:,:) real, allocatable :: lat_tile_subgrid(:,:) real, allocatable :: crossreference_emission_to_tile_subgrid(:,:,:,:) real, allocatable :: crossreference_population_to_tile_subgrid(:,:,:) integer, allocatable :: tile_class_subgrid(:,:) integer, allocatable :: tile_municipality_subgrid(:,:) real, allocatable :: aggregated_tile_subgrid(:,:,:) real, allocatable :: aggregated_x_tile_subgrid(:,:,:) real, allocatable :: aggregated_y_tile_subgrid(:,:,:) integer, allocatable :: aggregated_tile_class_subgrid(:,:,:) integer i,j,i_class,i_tile,j_tile,i_source,k integer n_tile_index integer tile_population_index !,tile_municipality_index,tile_class_index character(256) temp_name,temp_str,temp_str1 integer*8 ssb_id integer municipality_id real x_ssb,f_easting,ssb_dx,y_ssb,ssb_dy integer num_tiles_with_municipality logical exists integer unit_in integer count, index_val integer num_tile_classes(10) real resolution_tile_classes(10) character(8) count_str integer :: unit_tile=21 real population_tile_scale integer reduce_grid_class real aggregation_tile_scale(10) integer n_aggregated_tiles parameter (n_aggregated_tiles=4) !40to5km !parameter (n_aggregated_tiles=5) !160to10km integer dim_region_tiles parameter (dim_region_tiles=500) integer internal_region_index(dim_region_tiles),internal_region_id(dim_region_tiles) character(256) internal_region_name(dim_region_tiles) integer n_region_tiles real tile_region_min(dim_region_tiles,n_dim_index),tile_region_max(dim_region_tiles,n_dim_index) real aggregated_tile_subgrid_delta(n_dim_index,n_aggregated_tiles) real aggregated_tile_subgrid_min(n_dim_index,n_aggregated_tiles) real aggregated_tile_subgrid_max(n_dim_index,n_aggregated_tiles) integer aggregated_tile_subgrid_dim(n_dim_index,n_aggregated_tiles) logical :: save_as_seperate_files=.true. integer n_search parameter (n_search=5) character(16) search_str(n_search) real search_delta(n_search) integer temp_search integer :: io data search_str /'_1000m','_500m','_250m','_100m','_50m'/ data search_delta /1000.,500.,250.,100.,50./ write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Calculating regional tile distribution and resolution (uEMEP_set_region_tile_grids)' write(unit_logfile,'(A)') '================================================================' !Set tile index values n_tile_index=n_source_index+1 tile_population_index=n_source_index+1 !tile_municipality_index=n_source_index+2 !tile_class_index=n_source_index+3 !Reduce size of grids at and above this value reduce_grid_class=6 !Set to 6 does not reduce !Specify the tiling region to cover all of Norway population_tile_scale=1. !For 10 km !population_tile_scale=3. !For 20 km population_tile_scale=0.625 !For 5 km, to give 250 people /km^2 tile_subgrid_min(x_dim_index)=-70000.-40000 tile_subgrid_min(y_dim_index)=6440000.-40000 tile_subgrid_max(x_dim_index)=1110000.+40000. tile_subgrid_max(y_dim_index)=7950000.+40000. !40to5km tile class population_tile_scale=0.5 !For 5 km, to give 200 people /km^2 tile_subgrid_delta(x_dim_index)=5000. tile_subgrid_delta(y_dim_index)=5000. !Municipality tile classs population_tile_scale=0.01 tile_subgrid_delta(x_dim_index)=1000. tile_subgrid_delta(y_dim_index)=1000. limit_val_tile_population=limit_val_tile_population*population_tile_scale !Set all tile subgrids relative to the target subgrid tile_subgrid_dim(x_dim_index)=floor((tile_subgrid_max(x_dim_index)-tile_subgrid_min(x_dim_index))/tile_subgrid_delta(x_dim_index)) !New definition tile_subgrid_dim(y_dim_index)=floor((tile_subgrid_max(y_dim_index)-tile_subgrid_min(y_dim_index))/tile_subgrid_delta(y_dim_index)) !New definition tile_subgrid_dim(t_dim_index)=1 !Not used !aggregation_tile_scale if (n_aggregated_tiles.eq.5) then !Not in use aggregation_tile_scale(1)=16. aggregation_tile_scale(2)=8. aggregation_tile_scale(3)=4. aggregation_tile_scale(4)=2. aggregation_tile_scale(5)=1. elseif (n_aggregated_tiles.eq.4) then aggregation_tile_scale(1)=8. aggregation_tile_scale(2)=4. aggregation_tile_scale(3)=2. aggregation_tile_scale(4)=1. elseif (n_aggregated_tiles.eq.3) then aggregation_tile_scale(1)=4. aggregation_tile_scale(2)=2. aggregation_tile_scale(3)=1. endif do k=1,n_aggregated_tiles aggregated_tile_subgrid_min(:,k)=tile_subgrid_min(:) aggregated_tile_subgrid_max(:,k)=tile_subgrid_max(:) aggregated_tile_subgrid_delta(:,k)=tile_subgrid_delta(:)*aggregation_tile_scale(k) aggregated_tile_subgrid_dim(x_dim_index,k)=floor((aggregated_tile_subgrid_max(x_dim_index,k)-aggregated_tile_subgrid_min(x_dim_index,k))/aggregated_tile_subgrid_delta(x_dim_index,k)) aggregated_tile_subgrid_dim(y_dim_index,k)=floor((aggregated_tile_subgrid_max(y_dim_index,k)-aggregated_tile_subgrid_min(y_dim_index,k))/aggregated_tile_subgrid_delta(y_dim_index,k)) enddo !Allocate tile subgrids if (.not.allocated(tile_subgrid)) allocate (tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index),n_tile_index)) if (.not.allocated(x_tile_subgrid)) allocate (x_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(y_tile_subgrid)) allocate (y_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(lon_tile_subgrid)) allocate (lon_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(lat_tile_subgrid)) allocate (lat_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(crossreference_emission_to_tile_subgrid)) allocate (crossreference_emission_to_tile_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),2,n_source_index)) if (.not.allocated(crossreference_population_to_tile_subgrid)) allocate (crossreference_population_to_tile_subgrid(population_subgrid_dim(x_dim_index),population_subgrid_dim(y_dim_index),2)) if (.not.allocated(tile_class_subgrid)) allocate (tile_class_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(tile_municipality_subgrid)) allocate (tile_municipality_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index))) if (.not.allocated(aggregated_tile_subgrid)) allocate (aggregated_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index),n_aggregated_tiles)) if (.not.allocated(aggregated_x_tile_subgrid)) allocate (aggregated_x_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index),n_aggregated_tiles)) if (.not.allocated(aggregated_y_tile_subgrid)) allocate (aggregated_y_tile_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index),n_aggregated_tiles)) if (.not.allocated(aggregated_tile_class_subgrid)) allocate (aggregated_tile_class_subgrid(tile_subgrid_dim(x_dim_index),tile_subgrid_dim(y_dim_index),n_aggregated_tiles)) !Set to 0 all tile values tile_subgrid=0. tile_municipality_subgrid=0 tile_class_subgrid=0 !Read in SSB file containing gridded municipality ids !ssb_dx=1000. !ssb_dy=1000. f_easting=2.e6 !Search file name to define the grid size ssb_dx=0. ssb_dy=0. do k=1,n_search temp_search=index(filename_population(municipality_index),trim(adjustl(search_str(k)))) if (temp_search.ne.0) then ssb_dx=search_delta(k) ssb_dy=search_delta(k) write(unit_logfile,'(i,A)') temp_search,' Reading municipality masking data with resolution '//trim(adjustl(search_str(k))) endif enddo if (ssb_dx.eq.0) then write(unit_logfile,'(A)') 'Cannot find a valid SSB grid size. Stopping. '//trim(filename_population(municipality_index)) stop else write(unit_logfile,'(A,2f12.2)') 'Setting municipality SSB grid size (x,y) = ',ssb_dx,ssb_dy endif pathfilename_population(municipality_index)=trim(pathname_population(municipality_index))//trim(filename_population(municipality_index)) !Test existence of the heating filename. If does not exist then use default inquire(file=trim(pathfilename_population(municipality_index)),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: SSB file with municipality IDs does not exist: ', trim(pathfilename_population(municipality_index)) stop endif temp_name=pathfilename_population(municipality_index) !Open the file for reading unit_in=20 open(unit_in,file=temp_name,access='sequential',status='old',readonly) write(unit_logfile,'(a)') ' Opening SSB municipality file '//trim(temp_name) rewind(unit_in) !Read header SSBID0250M;kommunenum read(unit_in,'(A)') temp_str write(unit_logfile,'(A)') 'Header: '//trim(temp_str) count=0 do ssb_id=0;municipality_id=0 !Read in file string read(unit_in,'(A)',iostat=io) temp_str if (io /= 0) exit index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:);if (index_val.gt.1) read(temp_str1,*) ssb_id read(temp_str,*) municipality_id count=count+1 !if (mod(count,100000).eq.0) write(*,*) count,ssb_id,municipality_id !Convert id to grid centre coordinates that are already in UTM33 for SSB data x_ssb=floor(ssb_id/10000000.)-f_easting+ssb_dx/2. y_ssb=mod(ssb_id,10000000)+ssb_dy/2. !Find the tile this ssb grid is in i_tile=1+floor((x_ssb-tile_subgrid_min(x_dim_index))/tile_subgrid_delta(x_dim_index)) !New definition j_tile=1+floor((y_ssb-tile_subgrid_min(y_dim_index))/tile_subgrid_delta(y_dim_index)) !New definition !Add the municipality id to give the grid a value if (municipality_id.gt.0) then tile_municipality_subgrid(i_tile,j_tile)=municipality_id endif enddo close(unit_in) !Read municipality ID's and names pathfilename_region_id=trim(pathname_region_id)//trim(filename_region_id) !Test existence of the heating filename. If does not exist then use default inquire(file=trim(pathfilename_region_id),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: Region file with municipality IDs does not exist: ', trim(pathfilename_region_id) stop endif temp_name=pathfilename_region_id !Open the file for reading unit_in=20 open(unit_in,file=temp_name,access='sequential',status='old',readonly) write(unit_logfile,'(a)') ' Opening municipality ID file '//trim(temp_name) rewind(unit_in) !Read header index,ID,name read(unit_in,'(A)') temp_str write(unit_logfile,'(A)') 'Header: '//trim(temp_str) count=0 do count=count+1 read(unit_in,*,iostat=io) internal_region_index(count),internal_region_id(count),internal_region_name(count) if (io /= 0) exit enddo n_region_tiles=count close(unit_in) write(unit_logfile,'(A,i)') 'Number of municipalities: ',n_region_tiles !Determine tile subgrid do j=1,tile_subgrid_dim(y_dim_index) do i=1,tile_subgrid_dim(x_dim_index) x_tile_subgrid(i,j)=tile_subgrid_min(x_dim_index)+tile_subgrid_delta(x_dim_index)*(i-0.5) !New definition y_tile_subgrid(i,j)=tile_subgrid_min(y_dim_index)+tile_subgrid_delta(y_dim_index)*(j-0.5) !New definition call PROJ2LL(x_tile_subgrid(i,j),y_tile_subgrid(i,j),lon_tile_subgrid(i,j),lat_tile_subgrid(i,j),projection_attributes,projection_type) !if (projection_type.eq.RDM_projection_index) then ! call RDM2LL(y_tile_subgrid(i,j),x_tile_subgrid(i,j),lat_tile_subgrid(i,j),lon_tile_subgrid(i,j)) !elseif (projection_type.eq.UTM_projection_index) then ! call UTM2LL(utm_zone,y_tile_subgrid(i,j),x_tile_subgrid(i,j),lat_tile_subgrid(i,j),lon_tile_subgrid(i,j)) !endif !If the EMEP projection is lambert then set the proj coordinates to lambert, otherwise to lat-lon !if (EMEP_projection_type.eq.LCC_projection_index) then ! call lb2lambert2_uEMEP(xproj_tile_subgrid(i,j),yproj_tile_subgrid(i,j),lon_tile_subgrid(i,j),lat_tile_subgrid(i,j),EMEP_projection_attributes) !else ! xproj_tile_subgrid(i,j)=lon_tile_subgrid(i,j) ! yproj_tile_subgrid(i,j)=lat_tile_subgrid(i,j) !endif enddo enddo !Find grid max and min for each municipality do k=1,n_region_tiles tile_region_min(k,x_dim_index)=1e24 tile_region_max(k,x_dim_index)=-1e24 tile_region_min(k,y_dim_index)=1e24 tile_region_max(k,y_dim_index)=-1e24 do j=1,tile_subgrid_dim(y_dim_index) do i=1,tile_subgrid_dim(x_dim_index) if (tile_municipality_subgrid(i,j).eq.internal_region_id(k)) then if (x_tile_subgrid(i,j).lt.tile_region_min(k,x_dim_index)) tile_region_min(k,x_dim_index)=x_tile_subgrid(i,j) if (x_tile_subgrid(i,j).gt.tile_region_max(k,x_dim_index)) tile_region_max(k,x_dim_index)=x_tile_subgrid(i,j) if (y_tile_subgrid(i,j).lt.tile_region_min(k,y_dim_index)) tile_region_min(k,y_dim_index)=y_tile_subgrid(i,j) if (y_tile_subgrid(i,j).gt.tile_region_max(k,y_dim_index)) tile_region_max(k,y_dim_index)=y_tile_subgrid(i,j) endif enddo enddo !Add or subtract ssb grid size to allow for the poor resolution and make sure all of the kommune is within the tile tile_region_min(k,x_dim_index)=tile_region_min(k,x_dim_index)-ssb_dx*1.5 tile_region_max(k,x_dim_index)=tile_region_max(k,x_dim_index)+ssb_dx*1.5 tile_region_min(k,y_dim_index)=tile_region_min(k,y_dim_index)-ssb_dy*1.5 tile_region_max(k,y_dim_index)=tile_region_max(k,y_dim_index)+ssb_dy*1.5 enddo !Determine the aggregated tile subgrids, in UTM only do k=1,n_aggregated_tiles do j=1,aggregated_tile_subgrid_dim(y_dim_index,k) do i=1,aggregated_tile_subgrid_dim(x_dim_index,k) aggregated_x_tile_subgrid(i,j,k)=aggregated_tile_subgrid_min(x_dim_index,k)+aggregated_tile_subgrid_delta(x_dim_index,k)*(i-0.5) aggregated_y_tile_subgrid(i,j,k)=aggregated_tile_subgrid_min(y_dim_index,k)+aggregated_tile_subgrid_delta(y_dim_index,k)*(j-0.5) enddo enddo enddo !Create a cross reference grid do i_source=1,n_source_index if (calculate_source(i_source)) then do j=1,emission_subgrid_dim(y_dim_index,i_source) do i=1,emission_subgrid_dim(x_dim_index,i_source) crossreference_emission_to_tile_subgrid(i,j,x_dim_index,i_source)=1+floor((x_emission_subgrid(i,j,i_source)-tile_subgrid_min(x_dim_index))/tile_subgrid_delta(x_dim_index)) crossreference_emission_to_tile_subgrid(i,j,y_dim_index,i_source)=1+floor((y_emission_subgrid(i,j,i_source)-tile_subgrid_min(y_dim_index))/tile_subgrid_delta(y_dim_index)) enddo enddo endif enddo do j=1,population_subgrid_dim(y_dim_index) do i=1,population_subgrid_dim(x_dim_index) crossreference_population_to_tile_subgrid(i,j,x_dim_index)=1+floor((x_population_subgrid(i,j)-tile_subgrid_min(x_dim_index))/tile_subgrid_delta(x_dim_index)) crossreference_population_to_tile_subgrid(i,j,y_dim_index)=1+floor((y_population_subgrid(i,j)-tile_subgrid_min(y_dim_index))/tile_subgrid_delta(y_dim_index)) enddo enddo !Calculate the population within each tile do j=1,population_subgrid_dim(y_dim_index) do i=1,population_subgrid_dim(x_dim_index) i_tile=crossreference_population_to_tile_subgrid(i,j,x_dim_index) j_tile=crossreference_population_to_tile_subgrid(i,j,y_dim_index) if (i_tile.gt.0.and.j_tile.gt.0) then tile_subgrid(i_tile,j_tile,tile_population_index)=tile_subgrid(i_tile,j_tile,tile_population_index)+population_subgrid(i,j,population_data_type) endif enddo enddo !Calculate the veh.km within each tile i_source=traffic_index do j=1,emission_subgrid_dim(y_dim_index,i_source) do i=1,emission_subgrid_dim(x_dim_index,i_source) i_tile=crossreference_emission_to_tile_subgrid(i,j,x_dim_index,i_source) j_tile=crossreference_emission_to_tile_subgrid(i,j,y_dim_index,i_source) if (i_tile.gt.0.and.j_tile.gt.0) then tile_subgrid(i_tile,j_tile,traffic_index)=tile_subgrid(i_tile,j_tile,i_source)+(proxy_emission_subgrid(i,j,i_source,1))/1000. !ADT*km endif enddo enddo !Calculate the shipping emissions within each tile i_source=shipping_index do j=1,emission_subgrid_dim(y_dim_index,i_source) do i=1,emission_subgrid_dim(x_dim_index,i_source) i_tile=crossreference_emission_to_tile_subgrid(i,j,x_dim_index,i_source) j_tile=crossreference_emission_to_tile_subgrid(i,j,y_dim_index,i_source) if (i_tile.gt.0.and.j_tile.gt.0) then tile_subgrid(i_tile,j_tile,i_source)=tile_subgrid(i_tile,j_tile,i_source)+(proxy_emission_subgrid(i,j,i_source,1)) !emission for the time period endif enddo enddo !Write summary results num_tiles_with_municipality=0 do j_tile=1,tile_subgrid_dim(y_dim_index) do i_tile=1,tile_subgrid_dim(x_dim_index) if (tile_municipality_subgrid(i_tile,j_tile).gt.0) num_tiles_with_municipality=num_tiles_with_municipality+1 enddo enddo write(unit_logfile,'(a,i)') 'MUNICIPALITY TILE: ',num_tiles_with_municipality !Write summary results num_tiles_with_population=0 do j_tile=1,tile_subgrid_dim(y_dim_index) do i_tile=1,tile_subgrid_dim(x_dim_index) do i_class=1,n_tiles_population_classes-1 if (tile_subgrid(i_tile,j_tile,tile_population_index).gt.limit_val_tile_population(i_class) & .and.tile_subgrid(i_tile,j_tile,tile_population_index).le.limit_val_tile_population(i_class+1)) num_tiles_with_population(i_class)=num_tiles_with_population(i_class)+1 enddo i_class=n_tiles_population_classes if (tile_subgrid(i_tile,j_tile,tile_population_index).gt.limit_val_tile_population(i_class)) num_tiles_with_population(i_class)=num_tiles_with_population(i_class)+1 enddo enddo do i_class=1,n_tiles_population_classes-1 write(unit_logfile,'(a,i,f12.1,a,f12.1,i)') 'POPULATION TILE: ',i_class,limit_val_tile_population(i_class),' -',limit_val_tile_population(i_class+1),num_tiles_with_population(i_class) enddo i_class=n_tiles_population_classes write(unit_logfile,'(a,i,f12.1,a,a12,i)') 'POPULATION TILE: ',i_class,limit_val_tile_population(i_class),' <',' ',num_tiles_with_population(i_class) !Write summary results num_tiles_with_traffic=0 do j_tile=1,tile_subgrid_dim(y_dim_index) do i_tile=1,tile_subgrid_dim(x_dim_index) do i_class=1,n_tiles_traffic_classes-1 if (tile_subgrid(i_tile,j_tile,traffic_index).gt.limit_val_tile_traffic(i_class) & .and.tile_subgrid(i_tile,j_tile,traffic_index).le.limit_val_tile_traffic(i_class+1)) num_tiles_with_traffic(i_class)=num_tiles_with_traffic(i_class)+1 enddo i_class=n_tiles_traffic_classes if (tile_subgrid(i_tile,j_tile,traffic_index).gt.limit_val_tile_traffic(i_class)) num_tiles_with_traffic(i_class)=num_tiles_with_traffic(i_class)+1 enddo enddo do i_class=1,n_tiles_traffic_classes-1 write(unit_logfile,'(a,i,f12.1,a,f12.1,i)') 'TRAFFIC TILE: ',i_class,limit_val_tile_traffic(i_class),' -',limit_val_tile_traffic(i_class+1),num_tiles_with_traffic(i_class) enddo i_class=n_tiles_traffic_classes write(unit_logfile,'(a,i,f12.1,a,a12,i)') 'TRAFFIC TILE: ',i_class,limit_val_tile_traffic(i_class),' <',' ',num_tiles_with_traffic(i_class) !Write summary results num_tiles_with_shipping=0 do j_tile=1,tile_subgrid_dim(y_dim_index) do i_tile=1,tile_subgrid_dim(x_dim_index) if (tile_subgrid(i_tile,j_tile,shipping_index).gt.limit_val_tile_shipping(1)) num_tiles_with_shipping(1)=num_tiles_with_shipping(1)+1 enddo enddo write(unit_logfile,'(a,i)') 'SHIPPING TILE: ',num_tiles_with_shipping(1) num_tile_classes=0 !Class 1: 500 m. No emissions at all so interpolated, irrespective of population !Class 2: 250 m. Shipping emissions > 0 and Traffic < 1000 (1) or (population < 1000 and Traffic < 10000 (2)) !Class 3: 125 m. Traffic < 10000 (2) or (population > 1000 (2) population < 5000 (2) !Class 4: 50 m. Traffic > 10000 (2) and population > 5000 (3) tile_class_subgrid=0 do j_tile=1,tile_subgrid_dim(y_dim_index) do i_tile=1,tile_subgrid_dim(x_dim_index) !Only choose tiles that are part of a municipality if (tile_municipality_subgrid(i_tile,j_tile).gt.0) then tile_class_subgrid(i_tile,j_tile)=1 if (tile_subgrid(i_tile,j_tile,shipping_index).le.limit_val_tile_shipping(1).and. & tile_subgrid(i_tile,j_tile,traffic_index).le.limit_val_tile_traffic(1).and. & tile_subgrid(i_tile,j_tile,tile_population_index).le.limit_val_tile_population(2)) then tile_class_subgrid(i_tile,j_tile)=1 !No sources and population less than 100 elseif (tile_subgrid(i_tile,j_tile,shipping_index).ge.limit_val_tile_shipping(1).and. & tile_subgrid(i_tile,j_tile,traffic_index).le.limit_val_tile_traffic(2).and. & tile_subgrid(i_tile,j_tile,tile_population_index).ge.limit_val_tile_population(1)) then tile_class_subgrid(i_tile,j_tile)=2 !Little traffic but any shipping and any population elseif ((tile_subgrid(i_tile,j_tile,tile_population_index).gt.limit_val_tile_population(2).and. & tile_subgrid(i_tile,j_tile,tile_population_index).le.limit_val_tile_population(4)).or. & (tile_subgrid(i_tile,j_tile,traffic_index).ge.limit_val_tile_traffic(2).and. & tile_subgrid(i_tile,j_tile,tile_population_index).le.limit_val_tile_population(4))) then tile_class_subgrid(i_tile,j_tile)=3 !Population from 1000 to 5000 or road traffic elseif (tile_subgrid(i_tile,j_tile,tile_population_index).gt.limit_val_tile_population(4).and. & tile_subgrid(i_tile,j_tile,tile_population_index).le.limit_val_tile_population(5)) then tile_class_subgrid(i_tile,j_tile)=4 !Population > 10000 elseif (tile_subgrid(i_tile,j_tile,tile_population_index).gt.limit_val_tile_population(5)) then tile_class_subgrid(i_tile,j_tile)=4 !Population > 100000 endif !if (tile_class_subgrid(i_tile,j_tile).eq.2) tile_class_subgrid(i_tile,j_tile)=3 num_tile_classes(tile_class_subgrid(i_tile,j_tile))=num_tile_classes(tile_class_subgrid(i_tile,j_tile))+1 endif enddo enddo write(unit_logfile,'(a,i)') 'TILES OF CLASS 1 (500m): ',num_tile_classes(1) write(unit_logfile,'(a,i)') 'TILES OF CLASS 2 (250m): ',num_tile_classes(2) write(unit_logfile,'(a,i)') 'TILES OF CLASS 3 (125m): ',num_tile_classes(3) if (reduce_grid_class.eq.4) then write(unit_logfile,'(a,i)') 'TILES OF CLASS 4 ( 50m): ',num_tile_classes(4)*4 else write(unit_logfile,'(a,i)') 'TILES OF CLASS 4 ( 50m): ',num_tile_classes(4) endif if (reduce_grid_class.eq.5) then write(unit_logfile,'(a,i)') 'TILES OF CLASS 5 ( 25m): ',num_tile_classes(5)*4 else write(unit_logfile,'(a,i)') 'TILES OF CLASS 5 ( 25m): ',num_tile_classes(5) endif resolution_tile_classes(1)=500. resolution_tile_classes(2)=250. resolution_tile_classes(3)=125. resolution_tile_classes(4)=50. resolution_tile_classes(5)=25. !Save results in a single file temp_name=trim(pathname_tiles)//trim(filename_tiles) write(unit_logfile,'(a,2a)') 'Saving to: ', trim(temp_name) open(unit_tile,file=temp_name,access='sequential',status='unknown') count=0 do k=1,n_region_tiles count=count+1 write(unit_tile,'(a,i0.5)') 'tile_tag= '//trim(save_tile_tag)//'_',count !write(unit_tile,'(a,i0.5)') 'tile_tag= ',count write(unit_tile,'(a,a)') 'region_name= ',trim(internal_region_name(count)) write(unit_tile,'(a,i0.5)') 'region_id= ',internal_region_id(count) write(unit_tile,'(a,i0.5)') 'region_index= ',internal_region_index(count) write(unit_tile,'(a,f12.2)') 'subgrid_delta(x_dim_index)=',region_subgrid_delta write(unit_tile,'(a,f12.2)') 'subgrid_delta(y_dim_index)=',region_subgrid_delta write(unit_tile,'(a,f12.2)') 'subgrid_min(x_dim_index)=',tile_region_min(k,x_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_min(y_dim_index)=',tile_region_min(k,y_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_max(x_dim_index)=',tile_region_max(k,x_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_max(y_dim_index)=',tile_region_max(k,y_dim_index) enddo close(unit_tile) !Save results in multiple files if (save_as_seperate_files) then count=0 do k=1,n_region_tiles count=count+1 write(count_str,'(i8)') count temp_name=trim(pathname_tiles)//trim(ADJUSTL(count_str))//'_'//trim(filename_tiles) open(unit_tile,file=temp_name,access='sequential',status='unknown') write(unit_tile,'(a,i0.5)') 'tile_tag= '//trim(save_tile_tag)//'_',count write(unit_tile,'(a,a)') 'region_name= ',trim(internal_region_name(count)) write(unit_tile,'(a,i0.5)') 'region_id= ',internal_region_id(count) write(unit_tile,'(a,i0.5)') 'region_index= ',internal_region_index(count) write(unit_tile,'(a,f12.2)') 'subgrid_delta(x_dim_index)=',region_subgrid_delta write(unit_tile,'(a,f12.2)') 'subgrid_delta(y_dim_index)=',region_subgrid_delta write(unit_tile,'(a,f12.2)') 'subgrid_min(x_dim_index)=',tile_region_min(k,x_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_min(y_dim_index)=',tile_region_min(k,y_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_max(x_dim_index)=',tile_region_max(k,x_dim_index) write(unit_tile,'(a,f12.2)') 'subgrid_max(y_dim_index)=',tile_region_max(k,y_dim_index) close(unit_tile) enddo write(unit_logfile,'(a,i)') 'Number of region tiles: ',count endif write(unit_logfile,'(a)') ' Stopping after calculating tiles' stop end subroutine uEMEP_set_region_tile_grids