subroutine uEMEP_chemistry() ! Routine for doing the chemistry calculations in uEMEP 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 real :: FF_loc, distance_grid integer :: i_cross_integral, j_cross_integral, i_nc, j_nc real :: sum_p_bg_out, sum_p_out, count_p_out real :: max_p_bg_out, max_p_out, min_p_bg_out, min_p_out ! NB. Additional is calculated but not necessarily saved! real :: nox_bg_additional, no2_bg_additional, o3_bg_additional ! These are calculated in the Chemistry routine. Fist declared here. Are global variables 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 ( .not. allocated(comp_source_EMEP_subgrid)) then allocate(comp_source_EMEP_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_additional_subgrid)) then allocate(comp_source_EMEP_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 ! 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 write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Calculating chemistry for NO2 (uEMEP_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 comp_subgrid(:,:,:,no2_index) = 0 comp_subgrid(:,:,:,nox_index) = 0 comp_subgrid(:,:,:,o3_index) = 0 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 ! Before calculating travel time then include the other EMEP sources not downscaled ! Travel time is set to EMEP Grid_width/FFgrid do t = t_start, t_end do j = 1, subgrid_dim(y_dim_index) do i = 1, subgrid_dim(x_dim_index) i_cross_integral = crossreference_target_to_integral_subgrid(i,j,x_dim_index) j_cross_integral = crossreference_target_to_integral_subgrid(i,j,y_dim_index) FF_loc = 1.0 if (hourly_calculations) then FF_loc = max(FF_min_dispersion, meteo_subgrid(i_cross_integral,j_cross_integral,t,FFgrid_subgrid_index)) else if (annual_calculations) then FF_loc = max(FF_min_dispersion, 1.0/meteo_subgrid(i_cross_integral,j_cross_integral,t,inv_FFgrid_subgrid_index)) end if i_nc = crossreference_target_to_emep_subgrid(i,j,x_dim_index) j_nc = crossreference_target_to_emep_subgrid(i,j,y_dim_index) if (EMEP_projection_type .eq. LL_projection_index) then distance_grid = 111000.0*sqrt(dgrid_nc(lon_nc_index)*cos(var1d_nc(j_nc,lat_nc_index)*pi/180.0)*dgrid_nc(lat_nc_index)) else ! Assumed LCC or PS distance_grid = sqrt(dgrid_nc(lon_nc_index)*dgrid_nc(lat_nc_index)) end if end do end do end do sum_p_bg_out = 0.0 sum_p_out = 0.0 count_p_out = 0 max_p_bg_out = -1000.0; min_p_bg_out = 1000.0; max_p_out = -1000.0; min_p_out = 1000.0 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) nox_bg = subgrid(i,j,t,emep_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) if (EMEP_additional_grid_interpolation_size .gt. 0) then nox_bg_additional=subgrid(i,j,t,emep_additional_nonlocal_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) end if o3_bg = comp_EMEP_subgrid(i,j,t,o3_index) f_no2_loc = 0.0 nox_loc = 0.0 do i_source = 1, n_source_index if (calculate_source(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 if (calculate_emep_source(i_source) .and. .not. calculate_source(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 do if (abs(nox_loc) > epsilon0) then f_no2_loc = f_no2_loc/nox_loc else f_no2_loc = 0.0 end if 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 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 (EMEP_additional_grid_interpolation_size .gt. 0) then no2_bg_additional = comp_EMEP_subgrid(i,j,t,no2_index)*nox_bg_additional/subgrid(i,j,t,emep_subgrid_index,allsource_index,pollutant_loop_back_index(nox_nc_index)) if (no2_background_chemistry_scheme_flag .eq. 1) then call uEMEP_nonlocal_NO2_O3(nox_bg_additional, 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_additional, o3_bg_additional) else o3_bg_additional = 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_additional)) ! Conserve Ox when removing NO2 in the background end if comp_source_EMEP_additional_subgrid(i,j,t,o3_index,allsource_index) = o3_bg_additional comp_source_EMEP_additional_subgrid(i,j,t,no2_index,allsource_index) = no2_bg_additional end if ! Set the background O3 level. use all_source for the nonlocal. comp_source_EMEP_subgrid(i,j,t,o3_index,allsource_index) = o3_bg comp_source_EMEP_subgrid(i,j,t,no2_index,allsource_index) = no2_bg 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 sum_p_bg_out = sum_p_bg_out + p_bg_out sum_p_out = sum_p_out + p_out count_p_out = count_p_out + 1 max_p_bg_out = max(max_p_bg_out, p_bg_out); min_p_bg_out = min(min_p_bg_out, p_bg_out) max_p_out = max(max_p_out, p_out); min_p_out = min(min_p_out, p_out) comp_subgrid(i,j,t,o3_index) = o3_out comp_subgrid(i,j,t,no2_index) = no2_out comp_subgrid(i,j,t,nox_index) = nox_out else comp_subgrid(i,j,t,o3_index) = NODATA_value comp_subgrid(i,j,t,no2_index) = NODATA_value comp_subgrid(i,j,t,nox_index) = NODATA_value end if end do end do end do write(*,'(A,2f12.3)') 'P value (nonlocal,local): ', sum_p_bg_out/count_p_out, sum_p_out/count_p_out write(*,'(A,2f12.3)') 'P max (nonlocal,local): ', max_p_bg_out, max_p_out write(*,'(A,2f12.3)') 'P min (nonlocal,local): ', min_p_bg_out, min_p_out end subroutine uEMEP_chemistry