subroutine uEMEP_read_shipping_asi_data !Reads in the original ais raw data use uEMEP_definitions implicit none character(2048) temp_str character(256) temp_str1 integer unit_in logical :: exists integer count,index_val real ddlatitude,ddlongitude,totalnoxemission,totalparticulatematteremission real y_ship,x_ship integer i_ship_index,j_ship_index integer source_index,subsource_index integer t integer, allocatable :: count_subgrid(:,:,:) real, allocatable :: temp1_subgrid(:,:),temp2_subgrid(:,:),temp3_subgrid(:,:) integer i_pollutant integer :: io write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Reading shipping asi data (uEMEP_read_shipping_asi_data)' write(unit_logfile,'(A)') '================================================================' source_index=shipping_index n_subsource(source_index)=1 proxy_emission_subgrid(:,:,source_index,:)=0. t=1 allocate (count_subgrid(emission_subgrid_dim(x_dim_index,source_index),emission_subgrid_dim(y_dim_index,source_index),n_pollutant_loop)) count_subgrid=0 allocate (temp1_subgrid(emission_subgrid_dim(x_dim_index,source_index),emission_subgrid_dim(y_dim_index,source_index))) allocate (temp2_subgrid(emission_subgrid_dim(x_dim_index,source_index),emission_subgrid_dim(y_dim_index,source_index))) allocate (temp3_subgrid(emission_subgrid_dim(x_dim_index,source_index),emission_subgrid_dim(y_dim_index,source_index))) pathfilename_ship(1)=trim(pathname_ship(1))//trim(filename_ship(1)) if (use_aggregated_shipping_emissions_flag) pathfilename_ship(1)=trim(pathname_ship(1))//'Aggregated_'//trim(filename_ship(1)) !Test existence of the shipping filename. If does not exist then use default inquire(file=trim(pathfilename_ship(1)),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: Shipping file does not exist: ', trim(pathfilename_ship(1)) stop endif !Open the file for reading unit_in=20 open(unit_in,file=pathfilename_ship(1),access='sequential',status='old',readonly) write(unit_logfile,'(a)') ' Opening shipping file '//trim(pathfilename_ship(1)) rewind(unit_in) subsource_index=1 !Read header ddlatitude;ddlongitude;totalnoxemission;totalparticulatematteremission;fk_vessellloydstype;fk_ais_norwegianmainvesselcategory;date;time read(unit_in,'(A)') temp_str !write(*,*) trim(temp_str) count=0 do read(unit_in,'(A)',iostat=io) temp_str if (io /= 0) exit ddlatitude=0.;ddlongitude=0.;totalnoxemission=0.;totalparticulatematteremission=0. !Extract the values in the temp_str index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) read(temp_str1,*) ddlatitude !write (*,*) ddlatitude index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) read(temp_str1,*) ddlongitude !write (*,*) ddlongitude index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) !write(*,*) index_val,trim(temp_str1),trim(temp_str) if (index_val.gt.1) read(temp_str1,*) totalnoxemission !write (*,*) totalnoxemission index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) if (index_val.gt.1) read(temp_str1,*) totalparticulatematteremission !index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) !if (index_val.gt.1) read(temp_str1,*) temp_int !index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) !if (index_val.gt.1) read(temp_str1,*) temp_int !index_val=index(temp_str,';',back=.false.);temp_str1=temp_str(1:index_val-1);temp_str=temp_str(index_val+1:) !if (index_val.gt.1) read(temp_str1,*) temp_str2 !write (*,*) trim(temp_str1) !temp_str1=temp_str !if (len(temp_str1).gt.0) read(temp_str1,*) temp_str2 !write (*,*) trim(temp_str1) !write(*,*) count,ddlatitude,ddlongitude,totalnoxemission,totalparticulatematteremission count=count+1 !if (mod(count,100000).eq.0) write(*,*) count,ddlatitude,ddlongitude,totalnoxemission,totalparticulatematteremission if (totalnoxemission.gt.0.or.totalparticulatematteremission.gt.0) then !Convert to EMEP coordinates if it is to be saved. emission grids are already in the EMEP coordinate system !This will not work for lat lon as it is now written but will never be called either if (save_emissions_for_EMEP(shipping_index)) then if (EMEP_projection_type.eq.LCC_projection_index) then call lb2lambert2_uEMEP(x_ship,y_ship,ddlongitude,ddlatitude,EMEP_projection_attributes) !elseif (EMEP_projection_type.eq.LL_projection_index) then !lon_ship=ddlongitude !lat_ship=ddlatitude elseif (EMEP_projection_type.eq.PS_projection_index) then call LL2PS_spherical(x_ship,y_ship,ddlongitude,ddlatitude,EMEP_projection_attributes) endif else !Convert lat lon to utm coords if (projection_type.eq.UTM_projection_index) then call ll2utm(1,utm_zone,ddlatitude,ddlongitude,y_ship,x_ship) elseif (projection_type.eq.LTM_projection_index) then call ll2ltm(1,ltm_lon0,ddlatitude,ddlongitude,y_ship,x_ship) elseif (projection_type.eq.LAEA_projection_index) then call LL2LAEA(x_ship,y_ship,ddlongitude,ddlatitude,projection_attributes) endif endif !Find the grid index it belongs to i_ship_index=1+floor((x_ship-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index)) j_ship_index=1+floor((y_ship-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index)) !(x_subgrid(i,j)-subgrid_min(1))/+subgrid_delta(1)+1=i !Add to subgrid if (i_ship_index.ge.1.and.i_ship_index.le.emission_subgrid_dim(x_dim_index,source_index) & .and.j_ship_index.ge.1.and.j_ship_index.le.emission_subgrid_dim(y_dim_index,source_index)) then do i_pollutant=1,n_pollutant_loop if (totalnoxemission.gt.0.and.pollutant_loop_index(i_pollutant).eq.nox_nc_index) then proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)=proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)+totalnoxemission count_subgrid(i_ship_index,j_ship_index,i_pollutant)=count_subgrid(i_ship_index,j_ship_index,i_pollutant)+1 elseif (totalparticulatematteremission.gt.0.and.(pollutant_loop_index(i_pollutant).eq.pm25_nc_index.or.pollutant_loop_index(i_pollutant).eq.pm10_nc_index)) then proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)=proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)+totalparticulatematteremission count_subgrid(i_ship_index,j_ship_index,i_pollutant)=count_subgrid(i_ship_index,j_ship_index,i_pollutant)+1 endif enddo endif endif enddo write(unit_logfile,'(A,I)') 'Shipping counts = ',count do i_pollutant=1,n_pollutant_loop write(unit_logfile,'(A,es12.3)') 'Total emission '//trim(pollutant_file_str(pollutant_loop_index(i_pollutant)))//' = ',sum(proxy_emission_subgrid(1:emission_subgrid_dim(x_dim_index,source_index),1:emission_subgrid_dim(y_dim_index,source_index),source_index,i_pollutant)) enddo close(unit_in) deallocate (count_subgrid) deallocate (temp1_subgrid,temp2_subgrid,temp3_subgrid) end subroutine uEMEP_read_shipping_asi_data