subroutine uEMEP_subgrid_deposition(source_index)
use uEMEP_definitions
implicit none
integer source_index
integer i,j,ii,jj,tt,k
integer i_pollutant
!Define the target subgrid to be the same as the emission subgrid
real, allocatable :: target_subgrid(:,:,:)
real, allocatable :: target_deposition_subgrid(:,:,:,:)
real, allocatable :: x_target_subgrid(:,:)
real, allocatable :: y_target_subgrid(:,:)
real, allocatable :: traveltime_target_subgrid(:,:,:,:)
real :: target_subgrid_delta(2)
integer target_subgrid_dim(n_dim_index)
!define the temporary arrays for meteo
real, allocatable :: temp_FF_subgrid(:,:)
real, allocatable :: temp_FF_integral_subgrid(:,:)
real, allocatable :: trajectory_vector(:,:)
real, allocatable :: angle_diff(:,:)
real temp_sum_subgrid(n_pollutant_loop)
integer traj_max_index
logical valid_traj
real traj_step_size,x_loc,y_loc,FFgrid_loc,logz0_loc,u_star0_loc,FF10_loc,zc_loc,invL_loc
real ay_loc,by_loc,az_loc,bz_loc,sig_y_0_loc,sig_z_0_loc,sig_y_00_loc,sig_z_00_loc,h_emis_loc,z_rec_loc,sig_z_loc,sig_y_loc,h_mix_loc
real FF_loc,FF_zc_loc,precip_loc
real FF_integral_loc
integer n_plume_subgrid_max
parameter (n_plume_subgrid_max=100000)
integer plume_crossreference(n_plume_subgrid_max,2)
real plume_distance(n_plume_subgrid_max,2)
real temp_plume_distance(n_plume_subgrid_max)
integer sorted_plume_index(n_plume_subgrid_max)
real plume_source(n_pollutant_loop)
integer plume_count,max_plume_count
real xpos_limit,ypos_limit
real distance_subgrid,distance_subgrid_min,distance_emission_subgrid_min
integer i_target_start,i_target_end,j_target_start,j_target_end
integer t_start,t_end
integer i_cross,j_cross
integer i_cross_integral,j_cross_integral
integer i_cross_deposition,j_cross_deposition
integer i_cross_target_integral,j_cross_target_integral
integer subsource_index
real subgrid_internal,subgrid_internal_integral
real subgrid_internal_pollutant(n_pollutant_loop),drydepo_internal_pollutant(n_pollutant_loop),wetdepo_internal_pollutant(n_pollutant_loop)
real drydepo_internal_integral
real vertical_integral_internal,wetdepo_internal,wetdepo_internal_integral
real xpos_emission_subgrid,ypos_emission_subgrid
real xpos_area_max,xpos_area_min,ypos_area_max,ypos_area_min
real xpos_target_subgrid,ypos_target_subgrid
integer target_subgrid_dim_min(2),target_subgrid_dim_max(2)
real plume_vertical_integral(n_integral_subgrid_index)
real plume_vertical_integral_pollutant(n_integral_subgrid_index)
real, allocatable :: target_vertical_integral_subgrid(:,:,:)
integer i_integral
integer count
!Functions
!integer rargsort(n_plume_subgrid_max)
write(unit_logfile,'(A)') ''
write(unit_logfile,'(A)') '================================================================'
write(unit_logfile,'(A)') 'Calculating deposition and dispersion (uEMEP_subgrid_deposition)'
write(unit_logfile,'(A)') '================================================================'
!Set up the target subgrid to be the same as the emission subgrid for that source.
target_subgrid_dim(:)=emission_subgrid_dim(:,source_index)
target_subgrid_delta(:)=emission_subgrid_delta(:,source_index)
!Deposition and plume depletion must be calculated also in the buffer zone
if (calculate_source_depletion_flag) then
!Need to calculate on all the grids to get the depletion
target_subgrid_dim_min(:)=1;target_subgrid_dim_max(:)=target_subgrid_dim(1:2)
else
!Limit the target grid to slightly larger than concentration grid
target_subgrid_dim_min(x_dim_index)=-1+1+floor((subgrid_min(x_dim_index)-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index))
target_subgrid_dim_min(y_dim_index)=-1+1+floor((subgrid_min(y_dim_index)-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index))
target_subgrid_dim_max(x_dim_index)=+1+1+ceiling((subgrid_max(x_dim_index)-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index))
target_subgrid_dim_max(y_dim_index)=+1+1+ceiling((subgrid_max(y_dim_index)-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index))
endif
!Allocate the target subgrid to be the same as the emission subgrid but no time dimension
if (allocated(target_subgrid)) deallocate (target_subgrid)
if (allocated(target_deposition_subgrid)) deallocate (target_deposition_subgrid)
if (allocated(x_target_subgrid)) deallocate (x_target_subgrid)
if (allocated(y_target_subgrid)) deallocate (y_target_subgrid)
if (allocated(traveltime_target_subgrid)) deallocate (traveltime_target_subgrid)
if (.not.allocated(target_subgrid)) allocate (target_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),n_pollutant_loop))
if (.not.allocated(target_deposition_subgrid)) allocate (target_deposition_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),n_deposition_index,n_pollutant_loop))
if (.not.allocated(x_target_subgrid)) allocate (x_target_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index)))
if (.not.allocated(y_target_subgrid)) allocate (y_target_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index)))
if (.not.allocated(traveltime_target_subgrid)) allocate (traveltime_target_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),2,n_pollutant_loop))
if (adjust_wetdepo_integral_to_lowest_layer_flag.or.local_subgrid_method_flag.eq.1) then
if (allocated(target_vertical_integral_subgrid)) deallocate (target_vertical_integral_subgrid)
if (.not.allocated(target_vertical_integral_subgrid)) allocate (target_vertical_integral_subgrid(target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),3))
endif
x_target_subgrid(:,:)=x_emission_subgrid(:,:,source_index)
y_target_subgrid(:,:)=y_emission_subgrid(:,:,source_index)
!Allocate temporary wind speed subgrid
allocate (temp_FF_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index)))
allocate (temp_FF_integral_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index)))
allocate (angle_diff(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index)))
!Set the x and y position limits to coincide to half the EMEP grid (refered here as lon and lat but can be also LCC projection) times the number of grids
xpos_limit=dgrid_nc(lon_nc_index)/2.*EMEP_grid_interpolation_size*local_fraction_grid_size_scaling
ypos_limit=dgrid_nc(lat_nc_index)/2.*EMEP_grid_interpolation_size*local_fraction_grid_size_scaling
!Minimum distance for travel time calculation set to half of a grid diagonal weighted so the circle has the same area as the square with that diagonal
distance_subgrid_min=sqrt(subgrid_delta(x_dim_index)*subgrid_delta(x_dim_index)+subgrid_delta(y_dim_index)*subgrid_delta(y_dim_index))/2./sqrt(2.)*4./3.14159
!Minimum distance for dispersion set to half of an emission grid diagonal weighted so the circle has the same area as the square with that diagonal
distance_emission_subgrid_min=sqrt(emission_subgrid_delta(x_dim_index,source_index)*emission_subgrid_delta(x_dim_index,source_index) &
+emission_subgrid_delta(y_dim_index,source_index)*emission_subgrid_delta(y_dim_index,source_index))/2./sqrt(2.)*4./3.14159
!Set the subsource_index to 1 so no additional subsources in these routines
subsource_index=1
!Set local dispersion parameters to be used only in the annual calculation, overwritten in the hourly files
!if (annual_calculations) then
call uEMEP_set_dispersion_params_simple(source_index,subsource_index)
ay_loc=ay(source_index,subsource_index)
by_loc=by(source_index,subsource_index)
az_loc=az(source_index,subsource_index)
bz_loc=bz(source_index,subsource_index)
sig_y_00_loc=sig_y_00(source_index,subsource_index)
sig_z_00_loc=sig_z_00(source_index,subsource_index)
h_emis_loc=h_emis(source_index,subsource_index)
z_rec_loc=z_rec(source_index,subsource_index)
!endif
write(unit_logfile,'(a,i3)')'Calculating deposition and dispersion data for '//trim(source_file_str(source_index))
!Set the target grid loop variables
j_target_start=1;j_target_end=target_subgrid_dim(y_dim_index)
i_target_start=1;i_target_end=target_subgrid_dim(x_dim_index)
j_target_start=target_subgrid_dim_min(y_dim_index);j_target_end=target_subgrid_dim_max(y_dim_index)
i_target_start=target_subgrid_dim_min(x_dim_index);i_target_end=target_subgrid_dim_max(x_dim_index)
!write(*,*) j_target_start,j_target_end,i_target_start,i_target_end
!stop
!Set the start and end times of the loop
t_start=1;t_end=subgrid_dim(t_dim_index)
!Loop through the time
do tt=t_start,t_end
!Initialise the final grid to 0
subgrid(:,:,tt,proxy_subgrid_index,source_index,:)=0.
subgrid(:,:,tt,drydepo_local_subgrid_index,source_index,:)=0.
subgrid(:,:,tt,wetdepo_local_subgrid_index,source_index,:)=0.
integral_subgrid(:,:,tt,:,source_index,:)=0.
!deposition_subgrid(:,:,tt,:,:)=0.
!Initialise the target grid to 0
target_subgrid=0.
target_deposition_subgrid=0.
target_vertical_integral_subgrid=0.
!Set the last meteo data subgrid in the case when the internal time loop is used
if (.not.use_single_time_loop_flag) then
if (tt.gt.t_start) then
last_meteo_subgrid(:,:,:)=meteo_subgrid(:,:,tt-1,:)
else
last_meteo_subgrid(:,:,:)=meteo_subgrid(:,:,tt,:)
endif
endif
!Finds the angle difference between the current and last meteo field for dispersion and implements meandering if selected
do j_cross=1,integral_subgrid_dim(y_dim_index)
do i_cross=1,integral_subgrid_dim(x_dim_index)
if (hourly_calculations) then
call delta_wind_direction (i_cross,j_cross,tt,meteo_subgrid(i_cross,j_cross,tt,FF10_subgrid_index),angle_diff(i_cross,j_cross))
else
angle_diff(i_cross,j_cross)=0.
endif
enddo
enddo
!Create a temporary wind speed subgrid for each hour
call uEMEP_create_wind_field(temp_FF_subgrid,angle_diff,wind_level_flag,source_index,subsource_index,tt)
!Create an integral wind speed subgrid
call uEMEP_create_wind_field(temp_FF_integral_subgrid,angle_diff,wind_level_integral_flag,source_index,subsource_index,tt)
!Define the trajectory and its length
!Maxium number of trajectory steps and size of steps based on the integral (meteorology) loop size
if (use_trajectory_flag(source_index)) then
traj_step_size=min(integral_subgrid_delta(x_dim_index),integral_subgrid_delta(y_dim_index))*traj_step_scale
traj_max_index=floor(max(integral_subgrid_loop_index(x_dim_index),integral_subgrid_loop_index(y_dim_index))/traj_step_scale)
if (tt.eq.t_start) write(unit_logfile,'(a,f12.1,i)') 'Trajectory step (m) and dimensions: ', traj_step_size, traj_max_index
if (.not.allocated(trajectory_vector)) allocate(trajectory_vector(traj_max_index,2))
endif
max_plume_count=0
!Set up the deposition target subgrid based on the emission subgrid.
!This works when the deposition subgrid is larger than the target emission subgrid
!When the deposition subgrid is smaller than the target then we need to take the average of the grids
!This should be area weighted but is not
if (calculate_deposition_flag) then
if (deposition_subgrid_delta(x_dim_index).lt.target_subgrid_delta(x_dim_index)) then
do jj=j_target_start,j_target_end
do ii=i_target_start,i_target_end
i_cross_deposition=crossreference_emission_to_deposition_subgrid(ii,jj,x_dim_index,source_index)
j_cross_deposition=crossreference_emission_to_deposition_subgrid(ii,jj,y_dim_index,source_index)
!i_loop_deposition=floor(target_subgrid_delta(x_dim_index)/deposition_subgrid_delta(x_dim_index)/2.)
!target_deposition_subgrid(ii,jj,:,:)=target_deposition_subgrid(ii,jj,:,:)+deposition_subgrid(i_cross_deposition,j_cross_deposition,tt,:,:)
!write(*,*) deposition_subgrid(i_cross_deposition,j_cross_deposition,tt,:,:)
do i_pollutant=1,n_pollutant_loop
target_deposition_subgrid(ii,jj,vd_index,i_pollutant)=target_deposition_subgrid(ii,jj,vd_index,i_pollutant)+ &
area_weighted_extended_interpolation_function(x_deposition_subgrid,y_deposition_subgrid,deposition_subgrid(:,:,tt,vd_index,i_pollutant) &
,deposition_subgrid_dim(x_dim_index),deposition_subgrid_dim(y_dim_index),deposition_subgrid_delta(:),x_target_subgrid(ii,jj),y_target_subgrid(ii,jj),target_subgrid_delta)
enddo
enddo
enddo
else
do jj=1,emission_subgrid_dim(y_dim_index,source_index)
do ii=1,emission_subgrid_dim(x_dim_index,source_index)
i_cross_deposition=crossreference_emission_to_deposition_subgrid(ii,jj,x_dim_index,source_index)
j_cross_deposition=crossreference_emission_to_deposition_subgrid(ii,jj,y_dim_index,source_index)
target_deposition_subgrid(ii,jj,:,:)=deposition_subgrid(i_cross_deposition,j_cross_deposition,tt,:,:)
!write(*,*) deposition_subgrid(i_cross_deposition,j_cross_deposition,tt,:,:)
enddo
enddo
endif
endif
!Loop through all the emission subgrids
do jj=1,emission_subgrid_dim(y_dim_index,source_index)
do ii=1,emission_subgrid_dim(x_dim_index,source_index)
!Allocate temporary emission (plume source) for this grid and time
plume_source(:)=emission_subgrid(ii,jj,tt,source_index,:)
!Only calculate if there is an emissions
if (sum(plume_source(:)).ne.0) then
!Calculate the trajectory for this emission source
if (use_trajectory_flag(source_index)) then
trajectory_vector=NODATA_value
!Calculate the trajectory for this emission grid
call uEMEP_calculate_all_trajectory(x_emission_subgrid(ii,jj,source_index),y_emission_subgrid(ii,jj,source_index),tt, &
traj_max_index,traj_step_size,trajectory_vector(:,x_dim_index),trajectory_vector(:,y_dim_index))
else
!Create an artificial trajectory here for the straight line case by making the trajectory step size half the size of the integral grid
!Not working!!!!
traj_step_size=min(integral_subgrid_max(x_dim_index)-integral_subgrid_min(x_dim_index),integral_subgrid_max(y_dim_index)-integral_subgrid_min(y_dim_index))/2.
call uEMEP_calculate_all_trajectory(x_emission_subgrid(ii,jj,source_index),y_emission_subgrid(ii,jj,source_index),tt, &
traj_max_index,traj_step_size,trajectory_vector(:,x_dim_index),trajectory_vector(:,y_dim_index))
endif
!Set the integral meteorological grid position for the emission position
i_cross_integral=crossreference_emission_to_integral_subgrid(ii,jj,x_dim_index,source_index)
j_cross_integral=crossreference_emission_to_integral_subgrid(ii,jj,y_dim_index,source_index)
i_cross_integral=min(max(1,i_cross_integral),integral_subgrid_dim(x_dim_index))
j_cross_integral=min(max(1,j_cross_integral),integral_subgrid_dim(y_dim_index))
!Set the local wind speed and other parameters at emission position
FF_loc=temp_FF_subgrid(i_cross_integral,j_cross_integral)
FF_integral_loc=temp_FF_integral_subgrid(i_cross_integral,j_cross_integral)
!L_loc=1./meteo_subgrid(i_cross_integral,j_cross_integral,tt,invL_subgrid_index)
invL_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt,invL_subgrid_index)
FFgrid_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt,FFgrid_subgrid_index)
logz0_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt,logz0_subgrid_index)
u_star0_loc=max(meteo_subgrid(i_cross_integral,j_cross_integral,tt,ustar_subgrid_index),ustar_min)
FF10_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt,FF10_subgrid_index)
sig_y_00_loc=emission_properties_subgrid(ii,jj,emission_sigy00_index,source_index)
sig_z_00_loc=emission_properties_subgrid(ii,jj,emission_sigz00_index,source_index)
h_emis_loc=emission_properties_subgrid(ii,jj,emission_h_index,source_index)
!MORE COMPLICATED BECAUSE THE PLUME DEPLETION MUST ALSO BE CALCULATED OUTSIDE THE TARGET REGION!!
!SO THE LOOP NEEDS TO BE THE EMISSION LOOP NOT JUST THE TARGET LOOP
!Loop through the target grids and find the distance of each target grid along the plume centre line
if (calculate_deposition_flag.and.calculate_source_depletion_flag) then
!Set the size of the loop region around the emission cell to be up to emission_subgrid_loop_index
i_target_start=max(1,ii-emission_subgrid_loop_index(x_dim_index,source_index))
i_target_end=min(emission_subgrid_dim(x_dim_index,source_index),ii+emission_subgrid_loop_index(x_dim_index,source_index))
j_target_start=max(1,jj-emission_subgrid_loop_index(y_dim_index,source_index))
j_target_end=min(emission_subgrid_dim(y_dim_index,source_index),jj+emission_subgrid_loop_index(y_dim_index,source_index))
else
!Usea smaller area
i_target_start=max(1,ii-emission_subgrid_loop_index(x_dim_index,source_index))
i_target_end=min(emission_subgrid_dim(x_dim_index,source_index),ii+emission_subgrid_loop_index(x_dim_index,source_index))
j_target_start=max(1,jj-emission_subgrid_loop_index(y_dim_index,source_index))
j_target_end=min(emission_subgrid_dim(y_dim_index,source_index),jj+emission_subgrid_loop_index(y_dim_index,source_index))
endif
!Set the emission limits (EMEP projection ) surrounding the target grid
xpos_emission_subgrid=xproj_emission_subgrid(ii,jj,source_index)
ypos_emission_subgrid=yproj_emission_subgrid(ii,jj,source_index)
xpos_area_max=xpos_emission_subgrid+xpos_limit
xpos_area_min=xpos_emission_subgrid-xpos_limit
ypos_area_max=ypos_emission_subgrid+ypos_limit
ypos_area_min=ypos_emission_subgrid-ypos_limit
plume_count=0
plume_crossreference=0
plume_distance=0
!write(*,*) i_target_end-i_target_start,j_target_end-j_target_start
do j=j_target_start,j_target_end
do i=i_target_start,i_target_end
!Only calculate if it is within the local region distance
!Set the EMEP projection position of the emission grid. This guarantees that it extentds to the right distance
xpos_target_subgrid=xproj_emission_subgrid(i,j,source_index)
ypos_target_subgrid=yproj_emission_subgrid(i,j,source_index)
!Select only target grids within the predefined region
if (xpos_target_subgrid.ge.xpos_area_min.and.xpos_target_subgrid.le.xpos_area_max &
.and.ypos_target_subgrid.ge.ypos_area_min.and.ypos_target_subgrid.le.ypos_area_max) then
!Find the integral index for the target grid
i_cross_target_integral=crossreference_emission_to_integral_subgrid(i,j,x_dim_index,source_index)
j_cross_target_integral=crossreference_emission_to_integral_subgrid(i,j,y_dim_index,source_index)
i_cross_target_integral=min(max(1,i_cross_target_integral),integral_subgrid_dim(x_dim_index))
j_cross_target_integral=min(max(1,j_cross_target_integral),integral_subgrid_dim(y_dim_index))
!Set the mixing height at the average of the emission and target position
h_mix_loc=(meteo_subgrid(i_cross_integral,j_cross_integral,tt,hmix_subgrid_index)+meteo_subgrid(i_cross_target_integral,j_cross_target_integral,tt,hmix_subgrid_index))/2.
!Set the precipitation at the receptor grid
precip_loc=meteo_subgrid(i_cross_target_integral,j_cross_target_integral,tt,precip_subgrid_index)
!Find the minimum distance to the trajectory and check it is valid (downwind)
call uEMEP_minimum_distance_trajectory_fast(x_target_subgrid(i,j),y_target_subgrid(i,j), &
traj_max_index,traj_step_size,trajectory_vector(:,x_dim_index),trajectory_vector(:,y_dim_index),x_loc,y_loc,valid_traj)
!valid_traj=.false.
if (valid_traj) then
!Check if y_loc is within 4*sigma region based on the PG stability functions at emission position
call uEMEP_set_dispersion_sigma_PG(invL_loc,logz0_loc,sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
!Put the plume in the plume vector for later use
if (abs(y_loc).lt.3.*sig_y_loc) then
plume_count=plume_count+1
!Check number of subgrids in plume does not exceed allowable
if (plume_count.gt.n_plume_subgrid_max) then
write(unit_logfile,'(a,i)') 'WARNING: Number of subgrids for a plume exceeds the maximum value of: ', n_plume_subgrid_max
write(unit_logfile,'(a)') ' Exiting the plume but will continue calculating'
goto 10
else
plume_crossreference(plume_count,x_dim_index)=i
plume_crossreference(plume_count,y_dim_index)=j
plume_distance(plume_count,x_dim_index)=x_loc
plume_distance(plume_count,y_dim_index)=y_loc
endif
endif
endif
endif
enddo
enddo
max_plume_count=max(max_plume_count,plume_count)
!Sort the target grids from closest to furthest creating a crossreference index
!write(*,*) 'IN: ',plume_distance(1:plume_count,x_dim_index)
10 temp_plume_distance(1:plume_count)=plume_distance(1:plume_count,x_dim_index)
call rargsort(temp_plume_distance(1:plume_count),sorted_plume_index(1:plume_count),plume_count)
!write(*,*) 'OUT:',sorted_plume_index(1:plume_count)
!Set the plume source that will be depleted with deposition
plume_source=emission_subgrid(ii,jj,tt,source_index,:)
!Loop through the sorted array starting at the closest (which should always be the emission grid)
!write(*,*) ii,jj,plume_count
do k=1,plume_count
x_loc=plume_distance(sorted_plume_index(k),x_dim_index)
y_loc=plume_distance(sorted_plume_index(k),y_dim_index)
i=plume_crossreference(sorted_plume_index(k),x_dim_index)
j=plume_crossreference(sorted_plume_index(k),y_dim_index)
!write(*,*) x_loc,y_loc
!Select method for assigning sigma
if (stability_scheme_flag.eq.1) then
call uEMEP_set_dispersion_sigma_simple(sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
endif
if (stability_scheme_flag.eq.2) then
call uEMEP_set_dispersion_params_PG(invL_loc,source_index,subsource_index)
ay_loc=ay(source_index,subsource_index)
by_loc=by(source_index,subsource_index)
az_loc=az(source_index,subsource_index)
bz_loc=bz(source_index,subsource_index)
call uEMEP_set_dispersion_sigma_PG(invL_loc,logz0_loc,sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
endif
if (stability_scheme_flag.eq.3) then
!Set initial values for sigma. Initial sig_y is set here as well but is overridden by Kz dispersion
call uEMEP_set_dispersion_sigma_simple(sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
call uEMEP_set_dispersion_sigma_Kz(Kz_scheme,x_loc,sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,sig_z_loc,h_emis_loc,h_mix_loc,invL_loc,FF10_loc,10.,logz0_loc,emission_subgrid_delta(:,source_index),u_star0_loc,average_zc_h_in_Kz_flag,n_kz_iterations,sig_y_scaling_factor,sig_z_loc,sig_y_loc,FF_zc_loc)
!Add the meandering and change in wind angle to the plume since not included in Kz calculation
sig_y_loc=sig_y_loc+x_loc*angle_diff(i_cross_integral,j_cross_integral)
!Use the average of the emision height and zc to determine wind speed. Is set to true if wind_level_flag=6
!CHECK THIS
if (wind_level_flag.eq.6.or.wind_level_zc_flag) then
!Set the minimum wind speed
FF_loc=sqrt(FF_zc_loc*FF_zc_loc+FF_min_dispersion*FF_min_dispersion)
endif
if (wind_level_integral_flag.eq.6) then
!Set the minimum wind speed
FF_integral_loc=sqrt(FF_zc_loc*FF_zc_loc+FF_min_dispersion*FF_min_dispersion)
endif
endif
if (stability_scheme_flag.eq.4) then
write(unit_logfile,'(a,i)') 'Stability_scheme_flag=4 (Kz emulator) no longer an option. Stopping'
stop
!call uEMEP_set_dispersion_sigma_Kz_emulator(h_emis_loc,invL_loc,logz0_loc,h_mix_loc,sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
endif
!Adjust the height of the wind to the average of the emission and plume centre of mass height.
!This is already the case in the Kz calculation so not repeated here.
if (wind_level_flag.eq.6.and.stability_scheme_flag.ne.3) then
!if (wind_level_flag.eq.6) then
call z_centremass_gauss_func(sig_z_loc,h_emis_loc,h_mix_loc,zc_loc)
zc_loc=(h_emis_loc+zc_loc)/2.
call u_profile_neutral_val_func(zc_loc,FF10_loc,10.,h_mix_loc,exp(logz0_loc),FF_zc_loc,u_star0_loc)
FF_loc=sqrt(FF_zc_loc*FF_zc_loc+FF_min_dispersion*FF_min_dispersion)
endif
if (wind_level_integral_flag.eq.6.and.stability_scheme_flag.ne.3) then
!if (wind_level_flag.eq.6) then
call z_centremass_gauss_func(sig_z_loc,h_emis_loc,h_mix_loc,zc_loc)
zc_loc=(h_emis_loc+zc_loc)/2.
call u_profile_neutral_val_func(zc_loc,FF10_loc,10.,h_mix_loc,exp(logz0_loc),FF_zc_loc,u_star0_loc)
FF_integral_loc=sqrt(FF_zc_loc*FF_zc_loc+FF_min_dispersion*FF_min_dispersion)
endif
!Calculate the dispersion based on the derived sigmas
subgrid_internal=gauss_plume_cartesian_sigma_func(x_loc,y_loc,h_emis_loc,z_rec_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_loc)
!write(*,'(9es12.2)') x_loc,y_loc,h_emis_loc,z_rec_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_loc,subgrid_internal
!Only use half of the source grid for deposition and depletion
if (k.gt.1) then
!s/m3 *m2=s/m
subgrid_internal_integral=(subgrid_internal)*target_subgrid_delta(x_dim_index)*target_subgrid_delta(y_dim_index)
else
subgrid_internal_integral=(subgrid_internal)*target_subgrid_delta(x_dim_index)*0.5*target_subgrid_delta(y_dim_index) !Half the grid
endif
!Calculate the vertically integrated mass of the plume (s/m2) up to the lowest level and up to the mixing height
if (adjust_wetdepo_integral_to_lowest_layer_flag.and.calculate_deposition_flag) then
plume_vertical_integral(hsurf_integral_subgrid_index)=gauss_plume_cartesian_sigma_integral_func(x_loc,y_loc,h_emis_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_loc,0.,H_emep)*H_emep
plume_vertical_integral(hmix_integral_subgrid_index)=gauss_plume_cartesian_sigma_integral_func(x_loc,y_loc,h_emis_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_loc,0.,h_mix_loc)*h_mix_loc
!write(*,*) h_mix_loc,x_loc,sig_z_loc,vertical_integral_internal/(sig_z_loc*sqrt(2.*3.14159)/(2.*3.14159*sig_y_loc*sig_z_loc)*exp(-0.5*(y_loc*y_loc)/(sig_y_loc*sig_y_loc))/FF_loc*plume_source(i_pollutant) )
!write(*,*) plume_vertical_integral(1)/plume_vertical_integral(2),H_emep/h_mix_loc,H_emep/h_mix_loc/(plume_vertical_integral(1)/plume_vertical_integral(2))
!Calculate the average concentration in the lowest layer
!These two give the same results but the second is quicker. Can put it in a seperate subroutine
!vertical_integral_internal=gauss_plume_cartesian_sigma_integral_func(x_loc,y_loc,h_emis_loc,z_rec_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_loc,0.,h_mix_loc)*h_mix_loc*plume_source(i_pollutant)
vertical_integral_internal=exp(-0.5*(y_loc*y_loc)/(sig_y_loc*sig_y_loc))/FF_loc/(sqrt(2.*3.14159)*sig_y_loc)
!write(*,*) h_mix_loc,x_loc,sig_z_loc,vertical_integral_internal/(exp(-0.5*(y_loc*y_loc)/(sig_y_loc*sig_y_loc))/FF_loc/(sqrt(2.*3.14159)*sig_y_loc)*plume_source(i_pollutant))
endif
if (local_subgrid_method_flag.eq.1) then
plume_vertical_integral(hsurf_average_subgrid_index)=gauss_plume_cartesian_sigma_integral_func(x_loc,y_loc,h_emis_loc,sig_z_loc,sig_y_loc,h_mix_loc,FF_integral_loc,0.,H_emep)
endif
if (subgrid_internal_integral.gt.0) then
if (calculate_deposition_flag) then
do i_pollutant=1,n_pollutant_loop
!Multiply by the emissions to get the concentration prior to depletion (ug/m3)
subgrid_internal_pollutant(i_pollutant)=subgrid_internal*plume_source(i_pollutant)
!Calculate the dry deposition flux (ug/m3 *m/s=ug/s/m2)
drydepo_internal_pollutant(i_pollutant)=subgrid_internal_pollutant(i_pollutant)*target_deposition_subgrid(i,j,vd_index,i_pollutant)
!Calculate the dimensionless integral for plume depletion (s/m *m/s=1)
drydepo_internal_integral=subgrid_internal_integral*target_deposition_subgrid(i,j,vd_index,i_pollutant)
!Set the scavenging (s/m2 /m *m/s = /m2). 1e-3/3600 converts mm/hr to m/s
wetdepo_internal=vertical_integral_internal*wetdepo_scavanging_rate(pollutant_loop_index(i_pollutant))*(precip_loc/1000./3600.)
!Calculate the wet deposition flux ( /m2 ug/s = ug/m2/s)
wetdepo_internal_pollutant(i_pollutant)=wetdepo_internal*plume_source(i_pollutant)
!Calculate the dimensionless integral for plume depletion (/m2 * m2 = 1)
wetdepo_internal_integral=wetdepo_internal*target_subgrid_delta(x_dim_index)*target_subgrid_delta(y_dim_index)
!Make integrated concentration values
plume_vertical_integral_pollutant(:)=plume_vertical_integral(:)*plume_source(i_pollutant)
!Calculate the plume depletion by dry and wet deposition (ug/s)
if (calculate_source_depletion_flag) then
plume_source(i_pollutant)=plume_source(i_pollutant)*exp(-drydepo_internal_integral-wetdepo_internal_integral)
endif
!write(*,'(f12.2,4es12.2,f12.3)') x_loc,subgrid_internal_pollutant(i_pollutant),subgrid_internal,plume_source(i_pollutant),drydepo_internal_pollutant(i_pollutant)*target_subgrid_delta(x_dim_index)*target_subgrid_delta(y_dim_index),plume_source(i_pollutant)/emission_subgrid(ii,jj,tt,source_index,i_pollutant)
!Add to the dry deposition target grid
target_deposition_subgrid(i,j,drydepo_index,i_pollutant)=target_deposition_subgrid(i,j,drydepo_index,i_pollutant)+drydepo_internal_pollutant(i_pollutant)
!Add to the wet deposition target grid
target_deposition_subgrid(i,j,wetdepo_index,i_pollutant)=target_deposition_subgrid(i,j,wetdepo_index,i_pollutant)+wetdepo_internal_pollutant(i_pollutant)
!Add to the target concentration subgrid position
target_subgrid(i,j,i_pollutant)=target_subgrid(i,j,i_pollutant)+subgrid_internal_pollutant(i_pollutant)
enddo
else
!No deposition
do i_pollutant=1,n_pollutant_loop
!Multiply by the emissions to get the concentration prior to depletion (ug/m3)
subgrid_internal_pollutant(i_pollutant)=subgrid_internal*plume_source(i_pollutant)
!Add to the target concentration subgrid position
target_subgrid(i,j,i_pollutant)=target_subgrid(i,j,i_pollutant)+subgrid_internal_pollutant(i_pollutant)
!Make integrated concentration values
plume_vertical_integral_pollutant(:)=plume_vertical_integral(:)*plume_source(i_pollutant)
enddo
endif
!Determine the distance for the travel time calculation
distance_subgrid=sqrt(x_loc*x_loc+y_loc*y_loc)
distance_subgrid=max(distance_subgrid,distance_subgrid_min)
!Add to the travel time array
traveltime_target_subgrid(i,j,1,:)=traveltime_target_subgrid(i,j,1,:)+distance_subgrid/FF_loc*subgrid_internal_pollutant
traveltime_target_subgrid(i,j,2,:)=traveltime_target_subgrid(i,j,2,:)+subgrid_internal_pollutant
!Calculate the vertically integrated mass of the plume (s/m2) up to the lowest level and up to the mixing height
if ((calculate_deposition_flag.and.adjust_wetdepo_integral_to_lowest_layer_flag).or.local_subgrid_method_flag.eq.1) then
target_vertical_integral_subgrid(i,j,:)=target_vertical_integral_subgrid(i,j,:)+plume_vertical_integral_pollutant(:)
!write(*,'(2i,4es12.2)') i,j,plume_vertical_integral(1),target_vertical_integral_subgrid(i,j,1),plume_vertical_integral(2),target_vertical_integral_subgrid(i,j,2)
endif
endif
enddo !plume loop
endif!end if emission not 0
enddo!End of emission loop
enddo
write(unit_logfile,'(a,i)') 'Maximum number of subgrids per plume in this period and region = ',max_plume_count
!This ratio should be 1 in a well mixed layer. It can be used to estimate the average column concentration based on the surface layer concentration
!target_vertical_integral_subgrid(:,:,4)=target_vertical_integral_subgrid(:,:,1)/target_vertical_integral_subgrid(:,:,2)
!Interpolate target grid to the concentration subgrid
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
if (use_subgrid(i,j,source_index)) then
do i_pollutant=1,n_pollutant_loop
if (calculate_deposition_flag) then
subgrid(i,j,tt,drydepo_local_subgrid_index,source_index,i_pollutant)=subgrid(i,j,tt,drydepo_local_subgrid_index,source_index,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,target_deposition_subgrid(:,:,drydepo_index,i_pollutant) &
,target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_subgrid(i,j),y_subgrid(i,j))
subgrid(i,j,tt,wetdepo_local_subgrid_index,source_index,i_pollutant)=subgrid(i,j,tt,wetdepo_local_subgrid_index,source_index,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,target_deposition_subgrid(:,:,wetdepo_index,i_pollutant) &
,target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_subgrid(i,j),y_subgrid(i,j))
endif
subgrid(i,j,tt,proxy_subgrid_index,source_index,i_pollutant)=subgrid(i,j,tt,proxy_subgrid_index,source_index,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,target_subgrid(:,:,i_pollutant) &
,target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_subgrid(i,j),y_subgrid(i,j))
traveltime_subgrid(i,j,tt,1,i_pollutant)=traveltime_subgrid(i,j,tt,1,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,traveltime_target_subgrid(:,:,1,i_pollutant) &
,target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_subgrid(i,j),y_subgrid(i,j))
traveltime_subgrid(i,j,tt,2,i_pollutant)=traveltime_subgrid(i,j,tt,2,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,traveltime_target_subgrid(:,:,2,i_pollutant) &
,target_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_subgrid(i,j),y_subgrid(i,j))
enddo
else
subgrid(i,j,tt,proxy_subgrid_index,source_index,:)=NODATA_value
traveltime_subgrid(i,j,tt,:,:)=NODATA_value
endif
enddo
enddo
!Determine the final travel time
traveltime_subgrid(:,:,tt,3,:)=traveltime_subgrid(:,:,tt,1,:)/traveltime_subgrid(:,:,tt,2,:)
where (traveltime_subgrid(:,:,tt,2,:).eq.0) traveltime_subgrid(:,:,tt,3,:)=3600.*12.
!Place the vertically integrated values in the integral subgrid from the target grid, that is the same as the emission grid
!write(*,*) 'Interpolating vertically integrated values to the integral grid'
if ((calculate_deposition_flag.and.adjust_wetdepo_integral_to_lowest_layer_flag).or.local_subgrid_method_flag.eq.1) then
do j=1,integral_subgrid_dim(y_dim_index)
do i=1,integral_subgrid_dim(x_dim_index)
!i_cross_integral=crossreference_integral_to_emission_subgrid(i,j,x_dim_index,source_index)
!j_cross_integral=crossreference_integral_to_emission_subgrid(i,j,y_dim_index,source_index)
!i_cross_integral=max(1,i_cross_integral);i_cross_integral=min(integral_subgrid_dim(x_dim_index),i_cross_integral)
!j_cross_integral=max(1,j_cross_integral);j_cross_integral=min(integral_subgrid_dim(y_dim_index),j_cross_integral)
do i_integral=1,n_integral_subgrid_index
do i_pollutant=1,n_pollutant_loop
integral_subgrid(i,j,tt,i_integral,source_index,i_pollutant)=integral_subgrid(i,j,tt,i_integral,source_index,i_pollutant) &
+area_weighted_interpolation_function(x_target_subgrid,y_target_subgrid,target_vertical_integral_subgrid(:,:,i_integral) &
,target_subgrid_dim(x_dim_index),target_subgrid_dim(y_dim_index),target_subgrid_delta(:),x_integral_subgrid(i,j),y_integral_subgrid(i,j))
!write(*,*) i,j,i_integral,integral_subgrid(i,j,tt,i_integral,source_index,i_pollutant)
enddo
enddo
enddo
enddo
!Add to allsource
do i_integral=1,n_integral_subgrid_index
integral_subgrid(:,:,tt,i_integral,allsource_index,:)=integral_subgrid(:,:,tt,i_integral,allsource_index,:)+integral_subgrid(:,:,tt,i_integral,source_index,:)
enddo
endif
!Show mean outputs for checking
do i_pollutant=1,n_pollutant_loop
temp_sum_subgrid(i_pollutant)=0.
count=0
do j=1,subgrid_dim(y_dim_index)
do i=1,subgrid_dim(x_dim_index)
if (use_subgrid(i,j,source_index)) then
temp_sum_subgrid(i_pollutant)=temp_sum_subgrid(i_pollutant)+subgrid(i,j,tt,proxy_subgrid_index,source_index,i_pollutant)
count=count+1
endif
enddo
enddo
if (count.gt.0) then
temp_sum_subgrid(i_pollutant)=temp_sum_subgrid(i_pollutant)/count
else
temp_sum_subgrid(i_pollutant)=0
endif
write(unit_logfile,'(a,3f12.3)') 'Mean concentration '//trim(pollutant_file_str(pollutant_loop_index(i_pollutant)))//': ',temp_sum_subgrid(i_pollutant)
enddo
enddo!End time loop
if (allocated(trajectory_vector)) deallocate(trajectory_vector)
if (allocated(target_subgrid)) deallocate(target_subgrid)
if (allocated(target_deposition_subgrid)) deallocate(target_deposition_subgrid)
if (allocated(temp_FF_subgrid)) deallocate(temp_FF_subgrid)
if (allocated(temp_FF_integral_subgrid)) deallocate(temp_FF_integral_subgrid)
if (allocated(traveltime_target_subgrid)) deallocate(traveltime_target_subgrid)
if (allocated(target_vertical_integral_subgrid)) deallocate(target_vertical_integral_subgrid)
end subroutine uEMEP_subgrid_deposition