subroutine uEMEP_read_netcdf_population use uEMEP_definitions use netcdf implicit none integer status_nc integer i_split,j_split,n_delta_split integer i,j integer i_dim,id_nc character(256) var_name_nc_temp,dimname_temp integer var_id_nc,var_id_nc_temp real x_ssb,y_ssb integer i_ssb_index,j_ssb_index real delta_pop_nc(num_dims_population_nc) integer dim_id_nc(num_dims_population_nc) integer dim_length_population_nc(num_dims_population_nc) real y_pop,x_pop integer source_index logical :: exists !Temporary reading rvariables real, allocatable :: population_nc_dp(:,:,:) double precision, allocatable :: var2d_nc_dp(:,:) !If the data type is dwelling then this means it is for heating if (SSB_data_type.eq.dwelling_index) source_index=heating_index !Set the filename pathfilename_population(SSB_data_type)=trim(pathname_population(SSB_data_type))//trim(filename_population(SSB_data_type)) !Test existence. If does not exist then stop inquire(file=trim(pathfilename_population(SSB_data_type)),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: Netcdf file does not exist: ', trim(pathfilename_population(SSB_data_type)) write(unit_logfile,'(A)') ' STOPPING' stop endif !Open the netcdf file for reading write(unit_logfile,'(2A)') ' Opening netcdf file: ',trim(pathfilename_population(SSB_data_type)) status_nc = NF90_OPEN (pathfilename_population(SSB_data_type), nf90_nowrite, id_nc) if (status_nc .NE. NF90_NOERR) then write(unit_logfile,'(A,I)') 'ERROR opening netcdf file. Stopping: ',status_nc stop endif !Find the projection. If no projection then in lat lon coordinates !status_nc = NF90_INQ_VARID (id_nc,'projection_lambert',var_id_nc) !status_nc = NF90_INQ_VARID (id_nc,'mollweide',var_id_nc_temp) !if (status_nc.eq.NF90_NOERR) then ! population_nc_projection_type=mollweide_projection_index ! var_id_nc=var_id_nc_temp !endif status_nc = NF90_INQ_VARID (id_nc,'transverse_mercator',var_id_nc_temp) if (status_nc.eq.NF90_NOERR) then population_nc_projection_type=UTM_projection_index var_id_nc=var_id_nc_temp endif status_nc = NF90_INQ_VARID (id_nc,'projection_utm',var_id_nc_temp) if (status_nc.eq.NF90_NOERR) then population_nc_projection_type=UTM_projection_index var_id_nc=var_id_nc_temp endif status_nc = NF90_INQ_VARID (id_nc,'projection_lambert',var_id_nc_temp) if (status_nc.eq.NF90_NOERR) then population_nc_projection_type=LCC_projection_index var_id_nc=var_id_nc_temp endif status_nc = NF90_INQ_VARID (id_nc,'projection_ETRS89_LAEA',var_id_nc_temp) if (status_nc.eq.NF90_NOERR) then population_nc_projection_type=LAEA_projection_index var_id_nc=var_id_nc_temp endif if (population_nc_projection_type.ne.LL_projection_index) then !If there is a projection then read in the attributes. All these are doubles !status_nc = nf90_inquire_variable(id_nc, var_id_nc, natts = numAtts_projection) status_nc = nf90_get_att(id_nc, var_id_nc, 'Central_Meridian', population_nc_projection_attributes(1)) status_nc = nf90_get_att(id_nc, var_id_nc, 'False_Easting', population_nc_projection_attributes(2)) status_nc = nf90_get_att(id_nc, var_id_nc, 'False_Northing', population_nc_projection_attributes(3)) status_nc = nf90_get_att(id_nc, var_id_nc, 'longitude_of_prime_meridian', population_nc_projection_attributes(4)) status_nc = nf90_get_att(id_nc, var_id_nc, 'semi_major_axis', population_nc_projection_attributes(5)) status_nc = nf90_get_att(id_nc, var_id_nc, 'inverse_flattening', population_nc_projection_attributes(6)) write(unit_logfile,'(A,4i,6f12.2)') 'Reading projection (index, params) ',population_nc_projection_type,population_nc_projection_attributes(1:6) endif !Even if projection is Mollweide still use the lat lon positions !population_nc_projection_type=LL_projection_index !Find the (x,y) dimensions of the file. Use the meteo id's as these are x and y do i_dim=1,num_dims_population_nc status_nc = NF90_INQ_DIMID (id_nc,dim_name_population_nc(i_dim),dim_id_nc(i_dim)) status_nc = NF90_INQUIRE_DIMENSION (id_nc,dim_id_nc(i_dim),dimname_temp,dim_length_population_nc(i_dim)) if (status_nc .NE. NF90_NOERR) then write(unit_logfile,'(A,A,A,I)') 'No dimension information available for ',trim(dim_name_population_nc(i_dim)),' Setting to 1 with status: ',status_nc dim_length_population_nc(i_dim)=1 endif enddo write(unit_logfile,'(A,6I)') ' Size of population dimensions (x,y): ',dim_length_population_nc !Allocate the temporary arrays for lat,lon and population if (.not.allocated(population_nc_dp)) allocate (population_nc_dp(dim_length_population_nc(x_dim_nc_index),dim_length_population_nc(y_dim_nc_index),num_var_population_nc)) !Lat and lon if (.not.allocated(var2d_nc_dp)) allocate (var2d_nc_dp(max(dim_length_population_nc(x_dim_nc_index),dim_length_population_nc(y_dim_nc_index)),num_dims_population_nc)) !x and y !Read the x and y values to get the delta do i=1,num_dims_population_nc !Identify the variable name and ID in the nc file and read it var_name_nc_temp=dim_name_population_nc(i) status_nc = NF90_INQ_VARID (id_nc, trim(var_name_nc_temp), var_id_nc) if (status_nc .EQ. NF90_NOERR) then !status_nc = nf90_get_att(id_nc, var_id_nc, "units", unit_dim_meteo_nc(i)) status_nc = NF90_GET_VAR (id_nc, var_id_nc, var2d_nc_dp(1:dim_length_population_nc(i),i),start=(/1/),count=(/dim_length_population_nc(i)/)) else write(unit_logfile,'(A,A,A,I)') 'No information available for ',trim(var_name_nc_temp),' Status: ',status_nc endif enddo delta_pop_nc=var2d_nc_dp(2,:)-var2d_nc_dp(1,:) write(unit_logfile,'(A,2f12.2)') 'Population grid delta (m): ',delta_pop_nc do i=1,num_var_population_nc !Identify the variable name and ID in the nc file and read it var_name_nc_temp=var_name_population_nc(i) status_nc = NF90_INQ_VARID (id_nc, trim(var_name_nc_temp), var_id_nc) if (status_nc .EQ. NF90_NOERR) then !status_nc = nf90_get_att(id_nc, var_id_nc, "units", unit_dim_meteo_nc(i)) status_nc = NF90_GET_VAR (id_nc, var_id_nc, population_nc_dp(:,:,i),start=(/1,1/),count=(/dim_length_population_nc(x_dim_nc_index),dim_length_population_nc(y_dim_nc_index)/)) write(unit_logfile,'(2a,2f12.2)') 'Population variable min and max: ',trim(var_name_nc_temp),minval(population_nc_dp(:,:,i)),maxval(population_nc_dp(:,:,i)) else write(unit_logfile,'(A,A,A,I)') 'No information available for ',trim(var_name_nc_temp),' Status: ',status_nc endif enddo !Loop through the population data and put it in the population grid !Converting from lat lon to the subgrid coordinates and then finding the nearest neighbour !If population_subgrid_delta<delta_pop_nc then split the input grid into subgrids that are smaller than subgrid !Split any way into four even if they are the same grids for better interpolation n_delta_split=floor(sqrt(delta_pop_nc(x_dim_nc_index)**2+delta_pop_nc(y_dim_nc_index)**2)/sqrt(population_subgrid_delta(x_dim_index)**2+population_subgrid_delta(y_dim_index)**2)+.001) n_delta_split=10*max(n_delta_split,1) write(unit_logfile,'(A,i)') 'Population grids split into this many segments for numerical nearest neighbour interpolation: ',n_delta_split !Already defined as 0. Do it again population_subgrid(:,:,SSB_data_type)=0 do j=1,dim_length_population_nc(y_dim_nc_index) do i=1,dim_length_population_nc(x_dim_nc_index) if (population_nc_dp(i,j,population_nc_index).gt.0) then if (projection_type.eq.UTM_projection_index) then call ll2utm(1,utm_zone,population_nc_dp(i,j,lat_nc_index),population_nc_dp(i,j,lon_nc_index),y_pop,x_pop) !write(*,*) population_nc_dp(i,j,lat_nc_index),population_nc_dp(i,j,lon_nc_index),y_pop,x_pop elseif (projection_type.eq.LTM_projection_index) then call ll2ltm(1,ltm_lon0,population_nc_dp(i,j,lat_nc_index),population_nc_dp(i,j,lon_nc_index),y_pop,x_pop) elseif (projection_type.eq.LAEA_projection_index) then call LL2LAEA(x_pop,y_pop,population_nc_dp(i,j,lon_nc_index),population_nc_dp(i,j,lat_nc_index),projection_attributes) !call LAEA2LL(x_pop,y_pop,x_ssb,y_ssb,projection_attributes) !write(*,*) population_nc_dp(i,j,lon_nc_index)-x_ssb,population_nc_dp(i,j,lat_nc_index)-y_ssb else write(unit_logfile,'(A)') 'No supported projection available for population. Stopping ' stop endif !Loop through the n_delta_split do i_split=1,n_delta_split do j_split=1,n_delta_split x_ssb=x_pop-delta_pop_nc(x_dim_nc_index)*0.5+(i_split-0.5)*delta_pop_nc(x_dim_nc_index)/n_delta_split y_ssb=y_pop-delta_pop_nc(y_dim_nc_index)*0.5+(j_split-0.5)*delta_pop_nc(y_dim_nc_index)/n_delta_split if (SSB_data_type.eq.population_index) then i_ssb_index=1+floor((x_ssb-population_subgrid_min(x_dim_index))/population_subgrid_delta(x_dim_index)) j_ssb_index=1+floor((y_ssb-population_subgrid_min(y_dim_index))/population_subgrid_delta(y_dim_index)) if (i_ssb_index.ge.1.and.i_ssb_index.le.population_subgrid_dim(x_dim_index) & .and.j_ssb_index.ge.1.and.j_ssb_index.le.population_subgrid_dim(y_dim_index)) then population_subgrid(i_ssb_index,j_ssb_index,SSB_data_type)=population_subgrid(i_ssb_index,j_ssb_index,SSB_data_type) & +population_nc_dp(i,j,population_nc_index)/(n_delta_split**2) !write(*,'(2i,6f12.1)') i_split,j_split,x_pop,y_pop,x_ssb,y_ssb,(i_split-0.5)*delta_pop_nc(x_dim_nc_index)/n_delta_split,(j_split-0.5)*delta_pop_nc(y_dim_nc_index)/n_delta_split !write(*,*) i_split,j_split,i_ssb_index,j_ssb_index !write(*,*) i_ssb_index,j_ssb_index,SSB_data_type,population_subgrid(i_ssb_index,j_ssb_index,SSB_data_type) endif endif if (SSB_data_type.eq.dwelling_index) then i_ssb_index=1+floor((x_ssb-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index)) j_ssb_index=1+floor((y_ssb-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index)) if (i_ssb_index.ge.1.and.i_ssb_index.le.emission_subgrid_dim(x_dim_index,source_index) & .and.j_ssb_index.ge.1.and.j_ssb_index.le.emission_subgrid_dim(y_dim_index,source_index)) then proxy_emission_subgrid(i_ssb_index,j_ssb_index,source_index,:)=proxy_emission_subgrid(i_ssb_index,j_ssb_index,source_index,:) & +population_nc_dp(i,j,population_nc_index)/(n_delta_split**2) endif endif enddo enddo endif enddo enddo if (allocated(population_nc_dp)) deallocate (population_nc_dp) if (allocated(var2d_nc_dp)) deallocate (var2d_nc_dp) end subroutine uEMEP_read_netcdf_population