uEMEP model uEMEP_auto_subgrid Automatically creates a grid dependent on the distance to source
subroutine uEMEP_auto_subgrid() !! uEMEP model uEMEP_auto_subgrid !! Automatically creates a grid dependent on the distance to source integer :: i,j,k integer :: i_source integer :: t real :: max_use_subgrid_size(n_source_index) real :: use_subgrid_range(n_source_index) integer :: use_subgrid_step integer :: use_emission_subgrid_space integer :: i_cross, j_cross integer :: i_start, i_end, j_start, j_end integer :: i_start_k, i_end_k, j_start_k, j_end_k real :: sum_emission ! Exit if emission positions are not to be used and if one of the other auto positions is to be used if ( .not. use_emission_positions_for_auto_subgrid_flag(allsource_index) .and. (use_population_positions_for_auto_subgrid_flag .or. use_receptor_positions_for_auto_subgrid_flag)) then return end if ! Exit and fill subgrid if none of the auto subgrids are to be used if ( .not. use_emission_positions_for_auto_subgrid_flag(allsource_index) .and. .not. use_population_positions_for_auto_subgrid_flag .and. .not. use_receptor_positions_for_auto_subgrid_flag) then use_subgrid = .true. return end if write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Setting subgrids automatically (uEMEP_auto_subgrid)' write(unit_logfile,'(A)') '================================================================' ! Set all subgrids to do not use use_subgrid_val = 0 use_subgrid_interpolation_index = -1 ! Set time index used for emissions to 1, so only tests for the first hour if there are emissions t = 1 ! Set the maximum grid size max_use_subgrid_size = max_interpolation_subgrid_size use_subgrid_range = 8.0 use_subgrid_range(traffic_index) = 4.0 use_subgrid_range(shipping_index) = 12.0 n_use_subgrid_levels = 0 ! Sets the auto gridding if (subgrid_delta(x_dim_index) .eq. 100.0) then use_subgrid_step_delta(0) = 1 use_subgrid_step_delta(1) = 2 use_subgrid_step_delta(2) = 5 use_subgrid_step_delta(3) = 10 use_subgrid_step_delta(4) = 20 if (max_interpolation_subgrid_size .eq. 250.0) n_use_subgrid_levels = 1 if (max_interpolation_subgrid_size .eq. 500.0) n_use_subgrid_levels = 2 if (max_interpolation_subgrid_size .eq. 1000.0) n_use_subgrid_levels = 3 if (max_interpolation_subgrid_size .eq. 2000.0) n_use_subgrid_levels = 4 use_subgrid_range(traffic_index) = 6.0 else if (subgrid_delta(x_dim_index) .eq. 25.0) then use_subgrid_step_delta(0) = 1 use_subgrid_step_delta(1) = 2 use_subgrid_step_delta(2) = 4 use_subgrid_step_delta(3) = 10 use_subgrid_step_delta(4) = 20 use_subgrid_step_delta(5) = 40 use_subgrid_step_delta(5) = 80 if (max_interpolation_subgrid_size .eq. 250.0) n_use_subgrid_levels = 2 if (max_interpolation_subgrid_size .eq. 500.0) n_use_subgrid_levels = 3 if (max_interpolation_subgrid_size .eq. 1000.0) n_use_subgrid_levels = 4 if (max_interpolation_subgrid_size .eq. 2000.0) n_use_subgrid_levels = 5 use_subgrid_range(traffic_index) = 4.0 else if (subgrid_delta(x_dim_index) .eq. 50.0) then use_subgrid_step_delta(0) = 1 use_subgrid_step_delta(1) = 2 use_subgrid_step_delta(2) = 4 use_subgrid_step_delta(3) = 10 use_subgrid_step_delta(4) = 20 use_subgrid_step_delta(5) = 40 if (max_interpolation_subgrid_size .eq. 250.0) n_use_subgrid_levels = 2 if (max_interpolation_subgrid_size .eq. 500.0) n_use_subgrid_levels = 3 if (max_interpolation_subgrid_size .eq. 1000.0) n_use_subgrid_levels = 4 if (max_interpolation_subgrid_size .eq. 2000.0) n_use_subgrid_levels = 5 use_subgrid_range(traffic_index) = 6.0 else if (subgrid_delta(x_dim_index) .eq. 250.0) then use_subgrid_step_delta(0) = 1 use_subgrid_step_delta(1) = 2 use_subgrid_step_delta(2) = 4 use_subgrid_step_delta(3) = 8 if (max_interpolation_subgrid_size .eq. 250.) n_use_subgrid_levels = 0 if (max_interpolation_subgrid_size .eq. 500.) n_use_subgrid_levels = 1 if (max_interpolation_subgrid_size .eq. 1000.) n_use_subgrid_levels = 2 if (max_interpolation_subgrid_size .eq. 2000.) n_use_subgrid_levels = 3 use_subgrid_range(traffic_index) = 8.0 else write(unit_logfile,'(a,f12.2)') 'When auto gridding target grid sizes must be either 25, 50, 100 or 250 metres. Stopping because target grid size is:', subgrid_delta(x_dim_index) stop end if if (n_use_subgrid_levels(allsource_index).eq.0) then write(unit_logfile,'(a,f12.2)') 'When auto gridding maximum grid sizes must be either 250, 500, 1000 or 2000 metres. Stopping because max grid size is:',max_interpolation_subgrid_size stop 1 end if ! Set the number of levels to match this do i_source = 1, n_source_index if (calculate_source(i_source)) then write(*,*) 'Using auto subgrid for source ', trim(source_file_str(i_source)), use_emission_positions_for_auto_subgrid_flag(i_source) end if end do do i_source = 1, n_source_index if (calculate_source(i_source) .and. use_emission_positions_for_auto_subgrid_flag(i_source)) then write(unit_logfile,'(a,2f10.1,i10)') trim(source_file_str(i_source))//': maximum use subgrid size (m), grid range and number of levels: ', & max_use_subgrid_size(i_source), use_subgrid_range(i_source), n_use_subgrid_levels(i_source) ! Fill the interpolation index with the highest level do k = n_use_subgrid_levels(i_source), 0, -1 use_subgrid_step = 2**k use_subgrid_step = use_subgrid_step_delta(k) use_emission_subgrid_space = floor(use_subgrid_step/sqrt(emission_subgrid_delta(x_dim_index,i_source)*emission_subgrid_delta(y_dim_index,i_source) & /subgrid_delta(x_dim_index)/subgrid_delta(y_dim_index))/2.0*use_subgrid_range(i_source)) use_emission_subgrid_space = max(1, use_emission_subgrid_space) do j = 1, subgrid_dim(y_dim_index), use_subgrid_step do i = 1, subgrid_dim(x_dim_index), use_subgrid_step i_cross=crossreference_target_to_emission_subgrid(i,j,x_dim_index,i_source) j_cross=crossreference_target_to_emission_subgrid(i,j,y_dim_index,i_source) ! Search in the neighbourhood and find the sum of the emissions i_start = max(1, i_cross - use_emission_subgrid_space) i_end = min(emission_subgrid_dim(x_dim_index,i_source), i_cross + use_emission_subgrid_space) j_start = max(1, j_cross - use_emission_subgrid_space) j_end = min(emission_subgrid_dim(y_dim_index,i_source), j_cross + use_emission_subgrid_space) sum_emission = sum(proxy_emission_subgrid(i_start:i_end,j_start:j_end,i_source,:)) ! Select those with emission sums > 0 but always include the coarsest level everywhere if (sum_emission .gt. 0. .or. k .eq. n_use_subgrid_levels(i_source)) then use_subgrid_val(i,j,i_source) = 1 ! Label the grids for interpolation later i_start_k = max(1, i - use_subgrid_step/2) i_end_k = min(subgrid_dim(x_dim_index), i + use_subgrid_step/2) j_start_k = max(1, j - use_subgrid_step/2) j_end_k = min(subgrid_dim(y_dim_index), j + use_subgrid_step/2) use_subgrid_interpolation_index(i_start_k:i_end_k,j_start_k:j_end_k,i_source) = k end if end do end do end do use_subgrid_val(:,:,allsource_index)=use_subgrid_val(:,:,allsource_index)+use_subgrid_val(:,:,i_source) ! Check do j = 1, subgrid_dim(y_dim_index) do i = 1, subgrid_dim(x_dim_index) if (use_subgrid_interpolation_index(i,j,i_source) .lt. 0) write(*,*) i, j, use_subgrid_interpolation_index(i,j,i_source) end do end do end if end do use_subgrid_val = min(1, use_subgrid_val) ! Convert values to logical do i_source = 1, n_source_index if (use_emission_positions_for_auto_subgrid_flag(i_source)) then if (calculate_source(i_source) .or. i_source .eq. allsource_index) then do j = 1, subgrid_dim(y_dim_index) do i = 1, subgrid_dim(x_dim_index) if (use_subgrid_val(i,j,i_source) .gt. 0) then use_subgrid(i,j,i_source) = .true. else use_subgrid(i,j,i_source) = .false. end if end do end do end if else ! If not to be auto gridded then set use to true use_subgrid(:,:,i_source) = .true. end if end do do i_source = 1, n_source_index if (calculate_source(i_source) .and. use_emission_positions_for_auto_subgrid_flag(i_source)) then write(unit_logfile,'(a,2i10,f6.1)') 'Number of calculation subgrids for '//trim(source_file_str(i_source))//' (number, total, percent):', & sum(use_subgrid_val(:,:,i_source)), subgrid_dim(1)*subgrid_dim(2), sum(use_subgrid_val(:,:,i_source))*100.0/(subgrid_dim(1)*subgrid_dim(2)) end if end do end subroutine uEMEP_auto_subgrid