subroutine uEMEP_create_wind_field(temp_FF_subgrid,angle_diff,wind_level_flag_in,source_index_in,subsource_index_in,tt_in)
use uEMEP_definitions
implicit none
integer, intent(in) :: wind_level_flag_in,tt_in,source_index_in,subsource_index_in
integer j_cross,i_cross
integer i_cross_integral,j_cross_integral
integer ii,jj
real z0_temp,h_temp
real sig_y_00_loc,h_emis_loc,h_mix_loc,FF10_loc,x_loc,sig_z_loc,sig_y_loc,ff_loc,logz0_loc
real sig_z_00_loc,sig_y_0_loc,sig_z_0_loc,zc_loc,u_star0_loc
!define the temporary arrays for meteo
real, intent(out) :: temp_FF_subgrid(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index))
real, intent(in) :: angle_diff(integral_subgrid_dim(x_dim_index),integral_subgrid_dim(y_dim_index))
temp_FF_subgrid=0.
if (wind_level_flag_in.ne.5) then
do j_cross=1,integral_subgrid_dim(y_dim_index)
do i_cross=1,integral_subgrid_dim(x_dim_index)
z0_temp=exp(meteo_subgrid(i_cross,j_cross,tt_in,logz0_subgrid_index))
h_temp=h_emis(source_index_in,subsource_index_in)
if (annual_calculations.and.wind_level_flag_in.eq.1) then
temp_FF_subgrid(i_cross,j_cross)=1./meteo_subgrid(i_cross,j_cross,tt_in,inv_FFgrid_subgrid_index)
elseif (annual_calculations.and.wind_level_flag_in.eq.2) then
temp_FF_subgrid(i_cross,j_cross)=1./meteo_subgrid(i_cross,j_cross,tt_in,inv_FFgrid_subgrid_index)*(1.-(log((H_meteo+z0_temp)/z0_temp)-log((h_temp+z0_temp)/z0_temp))/log((H_meteo+z0_temp)/z0_temp))
elseif (annual_calculations.and.wind_level_flag_in.eq.3) then
temp_FF_subgrid(i_cross,j_cross)=1./meteo_subgrid(i_cross,j_cross,tt_in,inv_FF10_subgrid_index)
elseif (annual_calculations.and.wind_level_flag_in.eq.4) then
temp_FF_subgrid(i_cross,j_cross)=1./meteo_subgrid(i_cross,j_cross,tt_in,inv_FF10_subgrid_index)*(1.-(log((10.+z0_temp)/z0_temp)-log((h_temp+z0_temp)/z0_temp))/log((10.+z0_temp)/z0_temp))
elseif (hourly_calculations.and.wind_level_flag_in.eq.1) then
temp_FF_subgrid(i_cross,j_cross)=meteo_subgrid(i_cross,j_cross,tt_in,FFgrid_subgrid_index)
elseif (hourly_calculations.and.wind_level_flag_in.eq.2) then
temp_FF_subgrid(i_cross,j_cross)=meteo_subgrid(i_cross,j_cross,tt_in,FFgrid_subgrid_index)*(1.-(log((H_meteo+z0_temp)/z0_temp)-log((h_temp+z0_temp)/z0_temp))/log((H_meteo+z0_temp)/z0_temp))
elseif (hourly_calculations.and.wind_level_flag_in.eq.3) then
temp_FF_subgrid(i_cross,j_cross)=meteo_subgrid(i_cross,j_cross,tt_in,FF10_subgrid_index)
elseif (hourly_calculations.and.wind_level_flag_in.eq.4) then
temp_FF_subgrid(i_cross,j_cross)=meteo_subgrid(i_cross,j_cross,tt_in,FF10_subgrid_index)*(1.-(log((10.+z0_temp)/z0_temp)-log((h_temp+z0_temp)/z0_temp))/log((10.+z0_temp)/z0_temp))
elseif (wind_level_flag_in.eq.0) then
temp_FF_subgrid(i_cross,j_cross)=1.
elseif (wind_level_flag_in.eq.5) then
!Will set based on sigma z centre of mass
temp_FF_subgrid(i_cross,j_cross)=1.
elseif (wind_level_flag_in.eq.6) then
!Will set based on sigma z centre of mass and emission height
temp_FF_subgrid(i_cross,j_cross)=1.
else
write(unit_logfile,'(a)') 'No valid wind_level_flag_in selected. Stopping (uEMEP_subgrid_dispersion)'
stop
endif
!Setting a minimum value for wind for dispersion purposes (cannot be zero)
temp_FF_subgrid(i_cross,j_cross)=sqrt(temp_FF_subgrid(i_cross,j_cross)*temp_FF_subgrid(i_cross,j_cross)+FF_min_dispersion*FF_min_dispersion)
if (temp_FF_subgrid(i_cross,j_cross).eq.0) then
write(unit_logfile,'(a,2i)') 'Zero wind speed at integral grid (stopping): ',i_cross,j_cross
stop
endif
enddo
enddo
endif
!If wind level flag is set to 5, use of initial plume centre of mass, then set wind speed for each non-zero emission grid
if (wind_level_flag_in.eq.5) then
temp_FF_subgrid=0.
do jj=1,emission_subgrid_dim(y_dim_index,source_index_in)
do ii=1,emission_subgrid_dim(x_dim_index,source_index_in)
if (sum(emission_subgrid(ii,jj,tt_in,source_index_in,:)).ne.0) then
!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_in)
j_cross_integral=crossreference_emission_to_integral_subgrid(ii,jj,y_dim_index,source_index_in)
!Set the local variables
logz0_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt_in,logz0_subgrid_index)
FF10_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt_in,FF10_subgrid_index)
sig_y_00_loc=emission_properties_subgrid(ii,jj,emission_sigy00_index,source_index_in)
sig_z_00_loc=emission_properties_subgrid(ii,jj,emission_sigz00_index,source_index_in)
h_emis_loc=emission_properties_subgrid(ii,jj,emission_h_index,source_index_in)
h_mix_loc=meteo_subgrid(i_cross_integral,j_cross_integral,tt_in,hmix_subgrid_index)
if (annual_calculations) then
FF10_loc=1./meteo_subgrid(i_cross_integral,j_cross_integral,tt_in,inv_FF10_subgrid_index)
endif
!Set sig_0's at the emission position
x_loc=0.
call uEMEP_set_dispersion_sigma_simple(sig_z_00_loc,sig_y_00_loc,sigy_0_subgid_width_scale,emission_subgrid_delta(:,source_index_in),angle_diff(i_cross_integral,j_cross_integral),x_loc,sig_z_loc,sig_y_loc,sig_z_0_loc,sig_y_0_loc)
!Use the initial plume centre of mass to determine wind advection height
call z_centremass_gauss_func(sig_z_0_loc,h_emis_loc,h_mix_loc,zc_loc)
call u_profile_neutral_val_func(zc_loc,FF10_loc,10.,h_mix_loc,exp(logz0_loc),FF_loc,u_star0_loc)
!Set the minimum wind speed
FF_loc=sqrt(FF_loc*FF_loc+FF_min_dispersion*FF_min_dispersion)
temp_FF_subgrid(ii,jj)=FF_loc
!write(*,*) FF10_loc,FF_loc,zc_loc,sig_z_0_loc
endif
enddo
enddo
endif
end subroutine uEMEP_create_wind_field