!!! for now, just use no2/nox ratio of the first subsource
subroutine uEMEP_source_fraction_chemistry() ! Special source allocation for no2 based on leaving out one source at a time in the chemistry calculation ! This will always give a sum less, but not much less than, the total no2 ! This is normalised in order for it to be used ! Vhemistry scheme must have been run prior to implementing this integer :: i, j real :: nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature real :: nox_out, no2_out, o3_out, p_bg_out, p_out integer :: t, t_start, t_end integer :: i_source, i_subsource, emep_subsource integer :: i_pollutant logical :: nox_available = .false. integer :: i_integral, j_integral integer :: remove_source real :: sum_no2_source_subgrid, sum_o3_source_subgrid real, allocatable :: comp_source_temp_subgrid(:,:,:,:,:) real, allocatable :: comp_source_EMEP_temp_subgrid(:,:,:,:,:) ! additional delarations needed for the in-region calculations integer, parameter :: inregion_index = 1 integer, parameter :: outregion_index = 2 integer :: k real :: f_no2_isource, nox_loc_isource_total, nox_loc_isource_from_in_region, nox_loc_isource real :: no2_inandout_region(2) real :: o3_inandout_region(2) real :: sum_no2_inregion_outregion, sum_o3_inregion_outregion real :: nox_semiloc_isource, f_no2_bg if (trace_emissions_from_in_region .and. .not. calculate_EMEP_additional_grid_flag) then if (.not. allocated(comp_source_subgrid_from_in_region)) allocate(comp_source_subgrid_from_in_region(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) comp_source_subgrid_from_in_region = 0.0 if (.not. allocated(comp_semilocal_source_subgrid_from_in_region)) allocate(comp_semilocal_source_subgrid_from_in_region(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) comp_semilocal_source_subgrid_from_in_region = 0.0 end if ! Search for nox in the pollutants do i_pollutant = 1, n_pollutant_loop if (pollutant_loop_index(i_pollutant) .eq. nox_nc_index) nox_available = .true. end do ! Leave the chemistry routine if nox is not available if ( .not. nox_available) return if ( .not. allocated(comp_source_subgrid)) then allocate(comp_source_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) end if if (calculate_EMEP_additional_grid_flag) then if ( .not. allocated(comp_source_additional_subgrid)) then allocate(comp_source_additional_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) end if ! Temporary array for storing the comp_source_subgrid to avoid rewriting large parts of the routine when running the additional version if ( .not. allocated(comp_source_temp_subgrid)) then allocate(comp_source_temp_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) end if if ( .not. allocated(comp_source_EMEP_temp_subgrid)) then allocate(comp_source_EMEP_temp_subgrid(subgrid_dim(x_dim_index),subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_compound_index,n_source_index)) end if comp_source_temp_subgrid = comp_source_subgrid comp_source_EMEP_temp_subgrid = comp_source_EMEP_subgrid comp_source_EMEP_subgrid = comp_source_EMEP_additional_subgrid end if write(unit_logfile,'(A)') '--------------------------------------------------------' if (.not. calculate_EMEP_additional_grid_flag) write(unit_logfile,'(A)') 'Calculating chemistry source contributions for NO2 and O3 (uEMEP_source_fraction_chemistry)' if (calculate_EMEP_additional_grid_flag) write(unit_logfile,'(A)') 'Calculating additional nonlocal for NO2 and O3 (uEMEP_source_fraction_chemistry)' write(unit_logfile,'(A)') '--------------------------------------------------------' if (no2_chemistry_scheme_flag .eq. 0) then write(unit_logfile,'(A)') 'No chemistry used' else if (no2_chemistry_scheme_flag .eq. 1) then write(unit_logfile,'(A)') 'Photostationary state used' else if (no2_chemistry_scheme_flag .eq. 2) then write(unit_logfile,'(A)') 'Photochemistry with time scale used' else if (no2_chemistry_scheme_flag .eq. 3) then write(unit_logfile,'(A)') 'Romberg parameterisation used' else if (no2_chemistry_scheme_flag .eq. 4) then write(unit_logfile,'(A)') 'SRM parameterisation used' else if (no2_chemistry_scheme_flag .eq. 5) then write(unit_logfile,'(A)') 'During parameterisation used' end if t_start = 1 t_end = subgrid_dim(t_dim_index) i_subsource = 1 emep_subsource = 1 nox_bg=0.0; no2_bg = 0.0; o3_bg = 0.0; nox_loc = 0.0; f_no2_loc = 0.0; J_photo = 0.0; temperature=0.0 ! Weighted travel time already calculated do t = t_start, t_end do j = 1,subgrid_dim(y_dim_index) do i = 1,subgrid_dim(x_dim_index) if (use_subgrid(i,j,allsource_index)) then i_integral = crossreference_target_to_integral_subgrid(i,j,x_dim_index) j_integral = crossreference_target_to_integral_subgrid(i,j,y_dim_index) J_photo = meteo_subgrid(i_integral,j_integral,t,J_subgrid_index) temperature = meteo_subgrid(i_integral,j_integral,t,t2m_subgrid_index) if (calculate_EMEP_additional_grid_flag) then nox_bg = subgrid(i,j,t,emep_additional_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) else nox_bg = subgrid(i,j,t,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) end if o3_bg = comp_EMEP_subgrid(i,j,t,o3_index) do remove_source = 1, n_source_index if (calculate_source(remove_source) .or. remove_source .eq. allsource_index .or. (calculate_emep_source(remove_source) .and. .not. calculate_source(remove_source))) then f_no2_loc = 0.0 nox_loc = 0.0 do i_source = 1, n_source_index if (calculate_source(i_source)) then if (remove_source .ne. i_source) then do i_subsource = 1, n_subsource(i_source) f_no2_loc = f_no2_loc + emission_factor(no2_index,i_source,i_subsource)/emission_factor(nox_index,i_source,i_subsource)*subgrid(i,j,t,local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) nox_loc = nox_loc + subgrid(i,j,t,local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) end do end if end if ! Include the local EMEP that are not being downscaled if ( .not. calculate_EMEP_additional_grid_flag) then if (calculate_emep_source(i_source) .and. .not. calculate_source(i_source)) then if (remove_source .ne. i_source) then do i_subsource = 1, n_subsource(i_source) f_no2_loc = f_no2_loc + f_no2_emep*subgrid(i,j,t,emep_local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) nox_loc = nox_loc + subgrid(i,j,t,emep_local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) end do end if end if end if if (calculate_EMEP_additional_grid_flag) then ! If calculating the additional region then use the additional local EMEP not being downscaled if (calculate_emep_source(i_source) .and. .not. calculate_source(i_source)) then if (remove_source .ne. i_source) then do i_subsource = 1, n_subsource(i_source) f_no2_loc = f_no2_loc + f_no2_emep*subgrid(i,j,t,emep_additional_local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) nox_loc = nox_loc + subgrid(i,j,t,emep_additional_local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) end do end if end if !If calculating the additional region then include the difference BG-BG_additional to the local EMEP that is being downscaled if (calculate_source(i_source)) then if (remove_source .ne. i_source) then do i_subsource = 1, n_subsource(i_source) f_no2_loc = f_no2_loc + f_no2_emep* & (subgrid(i,j,t,emep_nonlocal_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) & - subgrid(i,j,t,emep_additional_nonlocal_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index))) nox_loc = nox_loc + subgrid(i,j,t,emep_nonlocal_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) & - subgrid(i,j,t,emep_additional_nonlocal_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) end do end if end if end if end do if (abs(nox_loc) > epsilon0) then f_no2_loc = f_no2_loc/nox_loc else f_no2_loc = 0.0 end if ! Use the all source index to calculate the contribution from the background ! This is done by removing all the sources, rather than the difference as done for the local sources ! This is because the chemistry is disturbed when removing background nox and no2 if (remove_source .ne. allsource_index) then no2_bg = comp_EMEP_subgrid(i,j,t,no2_index)*nox_bg/subgrid(i,j,t,emep_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) o3_bg = max(0.0, comp_EMEP_subgrid(i,j,t,o3_index) + 48.0/46.0*(comp_EMEP_subgrid(i,j,t,no2_index) - no2_bg)) ! Conserve Ox when removing NO2 in the background. Cannot be less than 0 else no2_bg = comp_EMEP_subgrid(i,j,t,no2_index)*nox_bg/subgrid(i,j,t,emep_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) o3_bg = max(0.0, comp_EMEP_subgrid(i,j,t,o3_index) + 48.0/46.0*(comp_EMEP_subgrid(i,j,t,no2_index) - no2_bg)) !Conserve Ox when removing NO2 in the background. Cannot be less than 0 nox_loc = 0.0 f_no2_loc = 0.0 end if ! Assume stationary state to derive no2 and o3 background. Overwrites the previous setting if (no2_background_chemistry_scheme_flag .eq. 1) then call uEMEP_nonlocal_NO2_O3(nox_bg, comp_EMEP_subgrid(i,j,t,nox_index), comp_EMEP_subgrid(i,j,t,no2_index), comp_EMEP_subgrid(i,j,t,o3_index), J_photo, temperature, f_no2_emep, no2_bg, o3_bg) end if if (no2_chemistry_scheme_flag .eq. 0) then nox_out = nox_bg + nox_loc no2_out = no2_bg + nox_loc*f_no2_loc o3_out = o3_bg else if (no2_chemistry_scheme_flag .eq. 1) then call uEMEP_photostationary_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, nox_out, no2_out, o3_out, p_bg_out, p_out) else if (no2_chemistry_scheme_flag .eq. 2) then call uEMEP_phototimescale_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, traveltime_subgrid(i,j,t,3,pollutant_loop_back_index(nox_nc_index))*traveltime_scaling, nox_out, no2_out, o3_out, p_bg_out, p_out) else if (no2_chemistry_scheme_flag .eq. 3) then call uEMEP_Romberg_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_out, no2_out, o3_out, romberg_parameters) else if (no2_chemistry_scheme_flag .eq. 4) then call uEMEP_SRM_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_out, no2_out, o3_out, SRM_parameters) else if (no2_chemistry_scheme_flag .eq. 5) then call uEMEP_During_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, comp_EMEP_subgrid(i,j,t,nox_index), comp_EMEP_subgrid(i,j,t,no2_index), comp_EMEP_subgrid(i,j,t,o3_index), J_photo, temperature, nox_out, no2_out, o3_out, p_bg_out, p_out) end if ! For background just use the result without any sources. ! There is a problem disturbing the chemistry by removing the background nox and no2 but not changing the o3 if (remove_source .eq. allsource_index) then comp_source_subgrid(i,j,t,no2_index,remove_source) = no2_bg comp_source_subgrid(i,j,t,o3_index,remove_source) = o3_bg else !Avoid round off errors which can occur with small numbers comp_source_subgrid(i,j,t,no2_index,remove_source) = max(0.0,comp_subgrid(i,j,t,no2_index) - no2_out) !Can be negative and can be greater than 1 so do not limit comp_source_subgrid(i,j,t,o3_index,remove_source) = comp_subgrid(i,j,t,o3_index) - o3_out end if end if end do !Normalise the contributions !Calculate the sum sum_no2_source_subgrid = 0.0 sum_o3_source_subgrid = 0.0 do i_source = 1, n_source_index if (calculate_source(i_source) .or. (calculate_emep_source(i_source) .and. .not. calculate_source(i_source))) then sum_no2_source_subgrid = sum_no2_source_subgrid + comp_source_subgrid(i,j,t,no2_index,i_source) sum_o3_source_subgrid = sum_o3_source_subgrid + comp_source_subgrid(i,j,t,o3_index,i_source) end if end do ! Set the background fractions so they will not be adjusted with normalisation do i_source = 1, n_source_index if (calculate_source(i_source) .or. (calculate_emep_source(i_source) .and. .not. calculate_source(i_source))) then ! Adjust for the background and normalise if (abs(sum_no2_source_subgrid) > epsilon0) then comp_source_subgrid(i,j,t,no2_index,i_source) = comp_source_subgrid(i,j,t,no2_index,i_source)/sum_no2_source_subgrid & *(comp_subgrid(i,j,t,no2_index) - comp_source_EMEP_subgrid(i,j,t,no2_index,allsource_index)) else comp_source_subgrid(i,j,t,no2_index,i_source) = 0 end if if (abs(sum_o3_source_subgrid) > epsilon0) then comp_source_subgrid(i,j,t,o3_index,i_source) = comp_source_subgrid(i,j,t,o3_index,i_source)/sum_o3_source_subgrid & *(comp_subgrid(i,j,t,o3_index) - comp_source_EMEP_subgrid(i,j,t,o3_index,allsource_index)) else comp_source_subgrid(i,j,t,o3_index,i_source) = 0 end if ! Setting local sources to 0 if total concentration is zero: No longer do this, because nonlocal might be non-zero even if total is zero !if (comp_subgrid(i,j,t,no2_index) .le. 0) comp_source_subgrid(i,j,t,no2_index,i_source) = 0 !if (comp_subgrid(i,j,t,o3_index) .le. 0) comp_source_subgrid(i,j,t,o3_index,i_source) = 0 end if end do ! Calculate NO2 and O3 source contributions from-in-region ! ******************************************************** if (trace_emissions_from_in_region .and. .not. calculate_EMEP_additional_grid_flag) then ! Calculate downscaled contributions from-in-region do remove_source = 1, n_source_index if (calculate_source(remove_source) .or. calculate_EMEP_source(remove_source)) then do k = 1,2 ! inregion and outregion ! add up all local sources to NOx, except 'remove_source' from either inregion or outregion f_no2_loc = 0 nox_loc = 0 do i_source = 1, n_source_index if (calculate_source(i_source) .or. calculate_EMEP_source(i_source)) then ! NB: loop over n_subsource is not logical, since we in practice weigh contributions from sources with 2 subsources more than sources with only 1, but I do it to follow Bruce's method above... do i_subsource = 1, n_subsource(i_source) ! check whether to use downscaled or non-downscaled local contribution for this source if (calculate_source(i_source)) then ! downscaled contribution f_no2_isource = emission_factor(no2_index,i_source,i_subsource)/emission_factor(nox_index,i_source,i_subsource) nox_loc_isource_total = subgrid(i,j,t,local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) nox_loc_isource_from_in_region = subgrid_local_from_in_region(i,j,t,i_source,pollutant_loop_back_index(nox_nc_index)) else ! i.e. calculate_EMEP_source(i_source) .and. .not. calculate_source(i_source) ! EMEP contribution f_no2_isource = f_no2_emep nox_loc_isource_total = subgrid(i,j,t,emep_local_subgrid_index,i_source,pollutant_loop_back_index(nox_nc_index)) nox_loc_isource_from_in_region = subgrid_EMEP_local_from_in_region(i,j,t,i_source,pollutant_loop_back_index(nox_nc_index)) end if ! check if this source is the one we remove or not, to determine how to add it if (i_source == remove_source) then ! for the source to remove: We should then treat the inregion and outregion differently if (k == inregion_index) then ! in region: add only the local contribution from outside region nox_loc_isource = nox_loc_isource_total - nox_loc_isource_from_in_region else if (k == outregion_index) then ! out region: add only the local contribution from inside region nox_loc_isource = nox_loc_isource_from_in_region else write(unit_logfile,'(A)') 'ERROR: value of k is not inregion_index or outregion_index' stop end if else ! for other sources, just add all the local contribution in both cases nox_loc_isource = nox_loc_isource_total end if f_no2_loc = f_no2_loc + f_no2_isource*nox_loc_isource nox_loc = nox_loc + nox_loc_isource end do end if end do ! divide f_no2_loc by total NOx contribution, in both cases if (abs(nox_loc) > epsilon0) then f_no2_loc = f_no2_loc/nox_loc else f_no2_loc = 0.0 end if ! Calculate background concentrations (following Bruce's approach above) no2_bg = comp_EMEP_subgrid(i,j,t,no2_index)*nox_bg/subgrid(i,j,t,emep_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) o3_bg = max(0.0, comp_EMEP_subgrid(i,j,t,o3_index) + 48.0/46.0*(comp_EMEP_subgrid(i,j,t,no2_index) - no2_bg)) !Conserve Ox when removing NO2 in the background. Cannot be less than 0 ! Assume stationary state to derive no2 and o3 background. Overwrites the previous setting if (no2_background_chemistry_scheme_flag .eq. 1) then call uEMEP_nonlocal_NO2_O3(nox_bg, comp_EMEP_subgrid(i,j,t,nox_index), comp_EMEP_subgrid(i,j,t,no2_index), comp_EMEP_subgrid(i,j,t,o3_index), J_photo, temperature, f_no2_emep, no2_bg, o3_bg) end if ! Calculate NO2 and O3 with the chemistry scheme if (no2_chemistry_scheme_flag .eq. 0) then nox_out = nox_bg + nox_loc no2_out = no2_bg + nox_loc*f_no2_loc o3_out = o3_bg else if (no2_chemistry_scheme_flag .eq. 1) then call uEMEP_photostationary_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, nox_out, no2_out, o3_out, p_bg_out, p_out) else if (no2_chemistry_scheme_flag .eq. 2) then call uEMEP_phototimescale_NO2(nox_bg, no2_bg, o3_bg, nox_loc, f_no2_loc, J_photo, temperature, traveltime_subgrid(i,j,t,3,pollutant_loop_back_index(nox_nc_index))*traveltime_scaling, nox_out, no2_out, o3_out, p_bg_out, p_out) else if (no2_chemistry_scheme_flag .eq. 3) then call uEMEP_Romberg_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_out, no2_out, o3_out, romberg_parameters) else if (no2_chemistry_scheme_flag .eq. 4) then call uEMEP_SRM_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, nox_out, no2_out, o3_out, SRM_parameters) else if (no2_chemistry_scheme_flag .eq. 5) then call uEMEP_During_NO2(nox_bg, no2_bg, nox_loc, o3_bg, f_no2_loc, comp_EMEP_subgrid(i,j,t,nox_index), comp_EMEP_subgrid(i,j,t,no2_index), comp_EMEP_subgrid(i,j,t,o3_index), J_photo, temperature, nox_out, no2_out, o3_out, p_bg_out, p_out) end if !Avoid round off errors which can occur with small numbers no2_inandout_region(k) = max(0.0, comp_subgrid(i,j,t,no2_index) - no2_out) !Can be negative and can be greater than 1 so do not limit o3_inandout_region(k) = comp_subgrid(i,j,t,o3_index) - o3_out end do ! k=1,2 ! scale the contributions so the sum equals the total contribution from the source sum_no2_inregion_outregion = no2_inandout_region(inregion_index) + no2_inandout_region(outregion_index) sum_o3_inregion_outregion = o3_inandout_region(inregion_index) + o3_inandout_region(outregion_index) if (abs(sum_no2_inregion_outregion) > epsilon0) then comp_source_subgrid_from_in_region(i,j,t,no2_index,remove_source) = comp_source_subgrid(i,j,t,no2_index,remove_source) * no2_inandout_region(inregion_index) / sum_no2_inregion_outregion else comp_source_subgrid_from_in_region(i,j,t,no2_index,remove_source) = 0 end if if (abs(sum_o3_inregion_outregion) > epsilon0) then comp_source_subgrid_from_in_region(i,j,t,o3_index,remove_source) = comp_source_subgrid(i,j,t,o3_index,remove_source) * o3_inandout_region(inregion_index) / sum_o3_inregion_outregion else comp_source_subgrid_from_in_region(i,j,t,o3_index,remove_source) = 0 end if end if end do ! remove_source = 1, n_source_index ! Calculate semilocal contributions to NO2 and O3 ! This is approximated by assuming the same NO2/NOx ratio as for background as a whole nox_bg = subgrid(i,j,t,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) no2_bg = comp_source_EMEP_subgrid(i,j,t,no2_index,allsource_index) do i_source = 1,n_source_index if (calculate_source(i_source) .or. calculate_EMEP_source(i_source)) then ! get initial NO2/NOx ratio of this source if (calculate_source(i_source)) then i_subsource = 1 !!!!! for now, just use no2/nox ratio of the first subsource f_no2_isource = emission_factor(no2_index,i_source,i_subsource)/emission_factor(nox_index,i_source,i_subsource) else f_no2_isource = f_no2_emep end if ! calculate NO2/NOx ratio in the background f_no2_bg = no2_bg / nox_bg ! get semilocal contribution to NOx from this source nox_semiloc_isource = subgrid_EMEP_semilocal_from_in_region(i,j,t,i_source,pollutant_loop_back_index(nox_index)) ! calculate NO2 and O3 semilocal contribution from the source, assuming the NO2/NOx ratio is the same as in background comp_semilocal_source_subgrid_from_in_region(i,j,t,no2_index,i_source) = f_no2_bg * nox_semiloc_isource comp_semilocal_source_subgrid_from_in_region(i,j,t,o3_index,i_source) = -48./46.*(f_no2_bg - f_no2_isource) * nox_semiloc_isource end if end do end if !(trace_emissions_from_in_region .and. .not. calculate_EMEP_additional_grid_flag) ! *************************************************************** ! done calculating NO2 and O3 source contributions from-in-region else ! i.e. if (.not. (use_subgrid(i,j,allsource_index))) comp_source_subgrid(i,j,t,:,:) = NODATA_value if (trace_emissions_from_in_region .and. .not. calculate_EMEP_additional_grid_flag) then comp_source_subgrid_from_in_region(i,j,t,:,:) = NODATA_value end if end if end do end do end do ! Transfer the arrays to the right outputs if (calculate_EMEP_additional_grid_flag) then comp_source_additional_subgrid = comp_source_subgrid comp_source_subgrid = comp_source_temp_subgrid comp_source_EMEP_subgrid = comp_source_EMEP_temp_subgrid ! EMEP_additional is unchanged if (allocated(comp_source_temp_subgrid)) deallocate(comp_source_temp_subgrid) if (allocated(comp_source_EMEP_temp_subgrid)) deallocate(comp_source_EMEP_temp_subgrid) end if end subroutine uEMEP_source_fraction_chemistry