subroutine uEMEP_calculate_emissions_for_EMEP implicit none !double precision :: EMEP_projection_attributes(10) !real, allocatable :: EMEP_emissions_grid(:,:,:,:,:) !x,y,t,source,pollutant integer a_start(6),date_array(6),a_start_emission(6) character(256) format_temp double precision date_num_temp,date_num_start,date_num_start_emission integer t,i_source,i,j write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Saving emission data for EMEP (uEMEP_calculate_emissions_for_EMEP)' write(unit_logfile,'(A)') '================================================================' !Read in example EMEP file for dimensions and projection information !Or just defne them here without reading if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then EMEP_projection_type=LCC_projection_index !grid_mapping_name = "lambert_conformal_conic"; !float i(i=531); !standard_name = "projection_x_coordinate"; !long_name = "x-coordinate in Cartesian system"; !units = "m"; !axis = "X"; if (trim(save_emissions_for_EMEP_region).eq.'NO') then emission_subgrid_min(x_dim_index,:)=save_emission_subgrid_min(x_dim_index) emission_subgrid_delta(x_dim_index,:)=save_emission_subgrid_delta(x_dim_index) emission_subgrid_dim(x_dim_index,:)=save_emission_subgrid_dim(x_dim_index) emission_max_subgrid_dim(x_dim_index)=save_emission_subgrid_dim(x_dim_index) !Move minium to edge, for consistency with normal subgrid definition emission_subgrid_min(x_dim_index,:)=emission_subgrid_min(x_dim_index,:)-emission_subgrid_delta(x_dim_index,:)/2. !float j(j=671); !standard_name = "projection_y_coordinate"; !long_name = "y-coordinate in Cartesian system"; !units = "m"; !axis = "Y"; emission_subgrid_min(y_dim_index,:)=save_emission_subgrid_min(y_dim_index) emission_subgrid_delta(y_dim_index,:)=save_emission_subgrid_delta(y_dim_index) emission_subgrid_dim(y_dim_index,:)=save_emission_subgrid_dim(y_dim_index) emission_max_subgrid_dim(y_dim_index)=save_emission_subgrid_dim(y_dim_index) !Move minium to edge, for consistency with normal subgrid definition emission_subgrid_min(y_dim_index,:)=emission_subgrid_min(y_dim_index,:)-emission_subgrid_delta(y_dim_index,:)/2. !if(use_meteo_file_for_emission_gridding_flag) then ! call uEMEP_read_meteo_nc ! emission_subgrid_min(x_dim_index,:)=-6.498834E+05 ! 6.751166E+05 ! emission_subgrid_delta(x_dim_index,:)=2500. ! emission_subgrid_dim(x_dim_index,:)=531 ! emission_max_subgrid_dim(x_dim_index)=531 !Move minium to edge, for consistency with normal subgrid definition ! emission_subgrid_min(x_dim_index,:)=emission_subgrid_min(x_dim_index,:)-emission_subgrid_delta(x_dim_index,:)/2. !endif else write(unit_logfile,'(A)') 'ERROR: Emission region '//trim(save_emissions_for_EMEP_region)//' not currently defined for lambert coordinates' stop endif elseif (trim(save_emissions_for_EMEP_projection).eq.'latlon') then EMEP_projection_type=LL_projection_index write(unit_logfile,'(A,i)') 'Projection of emission grid set to '//trim(save_emissions_for_EMEP_projection),EMEP_projection_type if (trim(save_emissions_for_EMEP_region).eq.'NL') then emission_subgrid_min(x_dim_index,:)=3.0 emission_subgrid_delta(x_dim_index,:)=0.25 emission_subgrid_dim(x_dim_index,:)=floor((8.-3.)/0.25)+1 emission_max_subgrid_dim(x_dim_index)=floor((8.-3.)/0.25)+1 !Full EMEP European domain emission_subgrid_min(x_dim_index,:)=-30. emission_subgrid_delta(x_dim_index,:)=0.25 emission_subgrid_dim(x_dim_index,:)=floor((45.+30.)/0.25)+1 emission_max_subgrid_dim(x_dim_index)=floor((45.+30.)/0.25)+1 !Move minium to edge, for consistency with normal subgrid definition emission_subgrid_min(x_dim_index,:)=emission_subgrid_min(x_dim_index,:)-emission_subgrid_delta(x_dim_index,:)/2. !float j(j=671); !standard_name = "projection_y_coordinate"; !long_name = "y-coordinate in Cartesian system"; !units = "m"; !axis = "Y"; emission_subgrid_min(y_dim_index,:)=50. emission_subgrid_delta(y_dim_index,:)=0.125 emission_subgrid_dim(y_dim_index,:)=floor((54.-50.)/0.125)+1 emission_max_subgrid_dim(y_dim_index)=floor((54.-50.)/0.125)+1 !Full EMEP European domain emission_subgrid_min(y_dim_index,:)=30. emission_subgrid_delta(y_dim_index,:)=0.125 emission_subgrid_dim(y_dim_index,:)=floor((76.-30.)/0.125)+1 emission_max_subgrid_dim(y_dim_index)=floor((76.-30.)/0.125)+1 !Move minium to edge, for consistency with normal subgrid definition emission_subgrid_min(y_dim_index,:)=emission_subgrid_min(y_dim_index,:)-emission_subgrid_delta(y_dim_index,:)/2. else write(unit_logfile,'(A)') 'ERROR: Emission region '//trim(save_emissions_for_EMEP_region)//' not currently defined for lat lon coordinates' stop endif else write(unit_logfile,'(A)') 'ERROR: Emission projection '//trim(save_emissions_for_EMEP_projection)//' not currently defined' stop endif !subgrid_dim(t_dim_index)=save_emissions_end_index-save_emissions_start_index+1 subgrid_dim(t_dim_index)=save_emissions_end_index dim_length_nc(x_dim_nc_index)=emission_subgrid_dim(x_dim_index,allsource_index) dim_length_nc(y_dim_nc_index)=emission_subgrid_dim(y_dim_index,allsource_index) dim_length_nc(time_dim_nc_index)=subgrid_dim(t_dim_index) write(unit_logfile,'(a,3i6)') 'EMEP emission grid dimensions: ', emission_subgrid_dim(x_dim_index,allsource_index),emission_subgrid_dim(y_dim_index,allsource_index),dim_length_nc(time_dim_nc_index) if (.not.allocated(val_dim_nc)) allocate (val_dim_nc(maxval(dim_length_nc),num_dims_nc)) !x, y, z and time dimension values if (.not.allocated(unit_dim_nc)) allocate (unit_dim_nc(num_dims_nc)) !x, y, z and time dimension values !Define the emission subgrid to correspond to the EMEP grid if (.not.allocated(emission_subgrid)) allocate (emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop)) if (.not.allocated(proxy_emission_subgrid)) allocate (proxy_emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_source_index,n_pollutant_loop)) if (.not.allocated(x_emission_subgrid)) allocate (x_emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_source_index)) if (.not.allocated(y_emission_subgrid)) allocate (y_emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_source_index)) if (.not.allocated(lon_emission_subgrid)) allocate (lon_emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_source_index)) if (.not.allocated(lat_emission_subgrid)) allocate (lat_emission_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_source_index)) if (.not.allocated(emission_time_profile_subgrid)) allocate (emission_time_profile_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),subgrid_dim(t_dim_index),n_source_index,n_pollutant_loop)) if (.not.allocated(emission_properties_subgrid)) allocate (emission_properties_subgrid(emission_max_subgrid_dim(x_dim_index),emission_max_subgrid_dim(y_dim_index),n_emission_index,n_source_index)) !Set the timing data. Assumes no single loops with the time based on the input date and array length defined by the length of the EMEP input data format_temp='yyyymmddHH' call datestr_to_date(config_date_str,format_temp,a_start) !Assumes it starts at hour 1 !a_start(4)=1 !Do not assume it starts at hour 1 any more a_start(5:6)=0 a_start_emission=a_start date_num_start=date_to_number(a_start,ref_year_EMEP) !Move the starting time according to the index value given (index is the number of hours) date_num_start_emission=date_num_start+dble(save_emissions_start_index-1)/24. !Set the emission_date_str to be used to name the output file as this may not be the same as the input file call number_to_date(date_num_start,a_start,ref_year_EMEP) call number_to_date(date_num_start_emission,a_start_emission,ref_year_EMEP) call date_to_datestr(a_start_emission,format_temp,emission_date_str) !Do not do this now emission_date_str=config_date_str !write(*,*) config_date_str,emission_date_str !long_name = "time at middle of period"; unit_dim_nc(time_dim_nc_index)="days since 1900-1-1 0:0:0"; do t=1,subgrid_dim(t_dim_index) date_num_temp=date_num_start+dble(t-1)/24.+dble(0.0001)/dble(24.)/dble(3600.) !Add 0.0001 of a second to avoid any rounding off errors call number_to_date(date_num_temp,date_array,ref_year_EMEP) write(*,'(6i6)') date_array(1:6) val_dim_nc(t,time_dim_nc_index)=date_num_temp enddo unit_dim_nc(x_dim_nc_index)='m' unit_dim_nc(y_dim_nc_index)='m' !Define emission grids do i_source=1,n_source_index if (save_emissions_for_EMEP(i_source)) then do j=1,emission_subgrid_dim(y_dim_index,i_source) do i=1,emission_subgrid_dim(x_dim_index,i_source) x_emission_subgrid(i,j,i_source)=emission_subgrid_min(x_dim_index,i_source)+emission_subgrid_delta(x_dim_index,i_source)*(i-0.5) y_emission_subgrid(i,j,i_source)=emission_subgrid_min(y_dim_index,i_source)+emission_subgrid_delta(y_dim_index,i_source)*(j-0.5) if (trim(save_emissions_for_EMEP_projection).eq.'lambert') then call lambert2lb2_uEMEP(x_emission_subgrid(i,j,i_source),y_emission_subgrid(i,j,i_source) & ,lon_emission_subgrid(i,j,i_source),lat_emission_subgrid(i,j,i_source),EMEP_projection_attributes) elseif (trim(save_emissions_for_EMEP_projection).eq.'latlon') then lon_emission_subgrid(i,j,i_source)=x_emission_subgrid(i,j,i_source) lat_emission_subgrid(i,j,i_source)=y_emission_subgrid(i,j,i_source) else write(unit_logfile,'(A)') 'ERROR: Emission projection '//trim(save_emissions_for_EMEP_projection)//' not currently defined' stop endif !write(*,*) x_emission_subgrid(i,j,i_source),y_emission_subgrid(i,j,i_source),lon_emission_subgrid(i,j,i_source),lat_emission_subgrid(i,j,i_source) enddo enddo if (i_source.eq.industry_index) then !Read in industry data call uEMEP_read_industry_data call uEMEP_read_time_profiles call uEMEP_set_emission_factors call uEMEP_convert_proxy_to_emissions endif if (i_source.eq.agriculture_index) then !Read agriculture data call uEMEP_read_agriculture_rivm_data call uEMEP_read_time_profiles call uEMEP_set_emission_factors call uEMEP_convert_proxy_to_emissions endif if (i_source.eq.traffic_index) then g_loop=1 !Read inthe road data call uEMEP_read_roadlink_data_ascii if (use_NORTRIP_emission_data) then call uEMEP_read_roadlink_emission_data endif !Redefine the road link positions to correspond to the EMEP coordinates !Grid the data. Road link coordinates will be redefined within this routine call uEMEP_grid_roads !call uEMEP_read_time_profiles !call uEMEP_set_emission_factors !call uEMEP_convert_proxy_to_emissions !Adjust traffic emissions of NOx based on temperature if (use_traffic_nox_emission_temperature_dependency) then use_alternative_meteorology_flag=.true. !Set the maximum dimension to that which is necessary. Minimum is not changed as it is selected in uEMEP_save_emission_netcdf end_time_meteo_nc_index=start_time_meteo_nc_index+(save_emissions_end_index-1) call uEMEP_read_meteo_nc !call uEMEP_subgrid_meteo_EMEP !call uEMEP_crossreference_grids call uEMEP_nox_emission_temperature endif endif if (i_source.eq.heating_index) then !meteo_var3d_nc(i_nc,j_nc,:,t2m_nc_index) !Read the heating data. Emission grid coordinates will be redefined within this routine !Must set g_loop=1 for it to read g_loop=1 call uEMEP_read_RWC_heating_data !Read in the temperature fields from the alternative meteorology always, since EMEP data should not exist yet use_alternative_meteorology_flag=.true. !Set the maximum dimension to that which is necessary. Minimum is not changed as it is selected in uEMEP_save_emission_netcdf !start_time_meteo_nc_index=start_time_meteo_nc_index+(save_emissions_start_index-1) end_time_meteo_nc_index=start_time_meteo_nc_index+(save_emissions_end_index-1) call uEMEP_read_meteo_nc !Need to make a cross reference here or simply skip the two based on an if statement call uEMEP_read_time_profiles call uEMEP_set_emission_factors call uEMEP_convert_proxy_to_emissions endif if (i_source.eq.shipping_index) then !meteo_var3d_nc(i_nc,j_nc,:,t2m_nc_index) !Read the heating data. Emission grid coordinates will be redefined within this routine !Must set g_loop=1 for it to read g_loop=1 if (read_weekly_shipping_data_flag) then call uEMEP_read_weekly_shipping_asi_data elseif (read_monthly_and_daily_shipping_data_flag) then call uEMEP_read_monthly_and_daily_shipping_asi_data else call uEMEP_read_shipping_asi_data endif !Read in the temperature fields from the alternative meteorology always, since EMEP data should not exist yet !use_alternative_meteorology_flag=.true. !call uEMEP_read_meteo_nc !Need to make a cross reference here or simply skip the two based on an if statement call uEMEP_read_time_profiles call uEMEP_set_emission_factors call uEMEP_convert_proxy_to_emissions endif endif enddo call uEMEP_save_emission_netcdf write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(a)') 'Finished saving emission data for EMEP.' write(unit_logfile,'(A)') '================================================================' stop end subroutine uEMEP_calculate_emissions_for_EMEP