subroutine uEMEP_read_netcdf_landuse_latlon use uEMEP_definitions use netcdf implicit none integer status_nc integer i,j integer i_dim,id_nc character(256) var_name_nc_temp,dimname_temp integer var_id_nc integer i_landuse_index,j_landuse_index real delta_landuse_nc(num_dims_landuse_nc) integer dim_id_nc(num_dims_landuse_nc) logical reduce_landuse_region_flag real temp_lon(4),temp_lat(4),temp_x(4),temp_y(4) real temp_x_min,temp_x_max,temp_y_min,temp_y_max integer i_temp_min,i_temp_max,j_temp_min,j_temp_max real temp_delta(num_dims_landuse_nc) real correct_lon(2) real temp_scale integer i_source,i_landuse real buffer_delta logical :: exists !Temporary reading variables real, allocatable :: landuse_nc_dp(:,:) double precision, allocatable :: var2d_nc_dp(:,:) double precision, allocatable :: temp_var2d_nc_dp(:,:) !Functions !real area_weighted_extended_vectorgrid_interpolation_function write(unit_logfile,'(A)') '' write(unit_logfile,'(A)') '================================================================' write(unit_logfile,'(A)') 'Reading landuse data (uEMEP_read_netcdf_landuse_latlon)' write(unit_logfile,'(A)') '================================================================' !Set the filename pathfilename_landuse=trim(pathname_landuse)//trim(filename_landuse) !Test existence. If does not exist then stop inquire(file=trim(pathfilename_landuse),exist=exists) if (.not.exists) then write(unit_logfile,'(A,A)') ' ERROR: Netcdf file does not exist: ', trim(pathfilename_landuse) write(unit_logfile,'(A)') ' STOPPING' stop endif !Open the netcdf file for reading write(unit_logfile,'(2A)') ' Opening netcdf file: ',trim(pathfilename_landuse) status_nc = NF90_OPEN (pathfilename_landuse, 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 (lon,lat) dimensions of the file. Use the meteo id's as these are x and y do i_dim=1,num_dims_landuse_nc status_nc = NF90_INQ_DIMID (id_nc,dim_name_landuse_nc(i_dim),dim_id_nc(i_dim)) status_nc = NF90_INQUIRE_DIMENSION (id_nc,dim_id_nc(i_dim),dimname_temp,dim_length_landuse_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_landuse_nc(i_dim)),' Setting to 1 with status: ',status_nc dim_length_landuse_nc(i_dim)=1 endif enddo write(unit_logfile,'(A,6I)') ' Size of landuse dimensions (lon,lat): ',dim_length_landuse_nc !Reduce the size of the grid to the heating emission grid size reduce_landuse_region_flag=.true. if (reduce_landuse_region_flag) then write(unit_logfile,'(A)') 'Reducing landuse domain for reading' !Determine the LL cordinates of the target grid !if (EMEP_projection_type.eq.LCC_projection_index) then !Retrieve the four corners of the target grid in lat and lon buffer_delta=10 call PROJ2LL(landuse_subgrid_min(x_dim_index)-buffer_delta*landuse_subgrid_delta(x_dim_index),landuse_subgrid_min(y_dim_index)-buffer_delta*landuse_subgrid_delta(y_dim_index),temp_lon(1),temp_lat(1),projection_attributes,projection_type) call PROJ2LL(landuse_subgrid_max(x_dim_index)+buffer_delta*landuse_subgrid_delta(x_dim_index),landuse_subgrid_max(y_dim_index)+buffer_delta*landuse_subgrid_delta(y_dim_index),temp_lon(2),temp_lat(2),projection_attributes,projection_type) call PROJ2LL(landuse_subgrid_min(x_dim_index)-buffer_delta*landuse_subgrid_delta(x_dim_index),landuse_subgrid_max(y_dim_index)+buffer_delta*landuse_subgrid_delta(y_dim_index),temp_lon(3),temp_lat(3),projection_attributes,projection_type) call PROJ2LL(landuse_subgrid_max(x_dim_index)+buffer_delta*landuse_subgrid_delta(x_dim_index),landuse_subgrid_min(y_dim_index)-buffer_delta*landuse_subgrid_delta(y_dim_index),temp_lon(4),temp_lat(4),projection_attributes,projection_type) temp_x_min=1.e32;temp_y_min=1.e32 temp_x_max=-1.e32;temp_y_max=-1.e32 temp_x=temp_lon;temp_y=temp_lat do i=1,4 write(*,*) i,temp_x(i),temp_y(i) if (temp_x(i).lt.temp_x_min) temp_x_min=temp_x(i) if (temp_y(i).lt.temp_y_min) temp_y_min=temp_y(i) if (temp_x(i).gt.temp_x_max) temp_x_max=temp_x(i) if (temp_y(i).gt.temp_y_max) temp_y_max=temp_y(i) enddo write(unit_logfile,'(A,2f12.2)') 'Min: ',temp_x_min,temp_y_min write(unit_logfile,'(A,2f12.2)') 'Max: ',temp_x_max,temp_y_max !Read the lon and lat values to get the delta and size. Put in temporary array !Allocate the temporary arrays for lat,lon and population if (.not.allocated(temp_var2d_nc_dp)) allocate (temp_var2d_nc_dp(max(dim_length_landuse_nc(x_dim_nc_index),dim_length_landuse_nc(y_dim_nc_index)),num_dims_landuse_nc)) !x and y dim_start_landuse_nc=1 do i=1,num_dims_landuse_nc !Identify the variable name and ID in the nc file and read it var_name_nc_temp=dim_name_landuse_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, temp_var2d_nc_dp(1:dim_length_landuse_nc(i),i),start=(/dim_start_landuse_nc(i)/),count=(/dim_length_landuse_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_landuse_nc=temp_var2d_nc_dp(2,:)-temp_var2d_nc_dp(1,:) write(unit_logfile,'(A,2f12.6)') 'Landuse grid delta (degrees): ',delta_landuse_nc !write(*,*) temp_var1d_nc_dp temp_delta(1)=delta_landuse_nc(1) temp_delta(2)=delta_landuse_nc(2) !write(*,*) temp_delta !Find grid position of the max and min coordinates and add2 grids*EMEP_grid_interpolation_size i_temp_min=1+floor((temp_x_min-temp_var2d_nc_dp(1,1))/temp_delta(1)+0.5) i_temp_max=1+floor((temp_x_max-temp_var2d_nc_dp(1,1))/temp_delta(1)+0.5) j_temp_min=1+floor((temp_y_min-temp_var2d_nc_dp(1,2))/temp_delta(2)+0.5) j_temp_max=1+floor((temp_y_max-temp_var2d_nc_dp(1,2))/temp_delta(2)+0.5) !write(unit_logfile,'(A,2I)') ' Reading EMEP i grids: ',i_temp_min,i_temp_max !write(unit_logfile,'(A,2I)') ' Reading EMEP j grids: ',j_temp_min,j_temp_max !Increase the region by 5 grids to be certain i_temp_min=max(1,i_temp_min-10) i_temp_max=min(dim_length_landuse_nc(x_dim_nc_index),i_temp_max+10) j_temp_min=max(1,j_temp_min-10) j_temp_max=min(dim_length_landuse_nc(y_dim_nc_index),j_temp_max+10) dim_length_landuse_nc(x_dim_nc_index)=i_temp_max-i_temp_min+1 dim_length_landuse_nc(y_dim_nc_index)=j_temp_max-j_temp_min+1 dim_start_landuse_nc(x_dim_nc_index)=i_temp_min dim_start_landuse_nc(y_dim_nc_index)=j_temp_min write(unit_logfile,'(A,3I)') ' Reading landuse i grids: ',i_temp_min,i_temp_max,dim_length_landuse_nc(x_dim_nc_index) write(unit_logfile,'(A,3I)') ' Reading landuse j grids: ',j_temp_min,j_temp_max,dim_length_landuse_nc(y_dim_nc_index) write(unit_logfile,'(A,2f12.2)') ' Reading landuse lon grids (min,max): ',temp_var2d_nc_dp(i_temp_min,x_dim_nc_index),temp_var2d_nc_dp(i_temp_max,x_dim_nc_index) write(unit_logfile,'(A,2f12.2)') ' Reading landuse lat grids (min,max): ',temp_var2d_nc_dp(j_temp_min,y_dim_nc_index),temp_var2d_nc_dp(j_temp_max,y_dim_nc_index) !endif endif if (i_temp_min.ge.i_temp_max.or.j_temp_min.ge.j_temp_max) then !No population data available write(unit_logfile,'(A)') ' WARNING: No landuse data available in this region. Setting to 0' landuse_subgrid(:,:,clc_index)=0 else if (.not.allocated(landuse_nc_dp)) allocate (landuse_nc_dp(dim_length_landuse_nc(x_dim_nc_index),dim_length_landuse_nc(y_dim_nc_index))) !Lat and lon if (.not.allocated(var2d_nc_dp)) allocate (var2d_nc_dp(max(dim_length_landuse_nc(x_dim_nc_index),dim_length_landuse_nc(y_dim_nc_index)),num_dims_landuse_nc)) !x and y !Read the lon and lat values to get the delta do i=1,num_dims_landuse_nc !Identify the variable name and ID in the nc file and read it var_name_nc_temp=dim_name_landuse_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_landuse_nc(i),i),start=(/dim_start_landuse_nc(i)/),count=(/dim_length_landuse_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_landuse_nc=var2d_nc_dp(2,:)-var2d_nc_dp(1,:) write(unit_logfile,'(A,2f12.6)') 'Landuse grid delta (degrees): ',delta_landuse_nc !write(*,*) var2d_nc_dp(1,1),var2d_nc_dp(dim_length_population_nc(x_dim_nc_index),1) !write(*,*) var2d_nc_dp(1,2),var2d_nc_dp(dim_length_population_nc(y_dim_nc_index),2) !Read the landuse data i=1 !Uses the population_nc_index as index, =1, but not logical !Identify the variable name and ID in the nc file and read it var_name_nc_temp=var_name_landuse_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, landuse_nc_dp(:,:),start=(/dim_start_landuse_nc(x_dim_nc_index),dim_start_landuse_nc(y_dim_nc_index)/),count=(/dim_length_landuse_nc(x_dim_nc_index),dim_length_landuse_nc(y_dim_nc_index)/)) write(unit_logfile,'(2a,2f12.2)') 'Landuse variable min and max: ',trim(var_name_nc_temp),minval(landuse_nc_dp(:,:)),maxval(landuse_nc_dp(:,:)) else write(unit_logfile,'(A,A,A,I)') 'No information available for ',trim(var_name_nc_temp),' Status: ',status_nc endif !enddo !write(*,*) 'Finished reading landuse data' !Loop through the landuse data and put it in the landuse grid !Converting from lat lon to the subgrid coordinates and then finding the nearest neighbour landuse_subgrid(:,:,clc_index)=0 where (landuse_nc_dp.lt.0) landuse_nc_dp=0. write(unit_logfile,'(2a,2f12.2)') 'Landuse min and max: ',trim(var_name_nc_temp),minval(landuse_nc_dp(:,:)),maxval(landuse_nc_dp(:,:)) !stop do j=1,landuse_subgrid_dim(y_dim_nc_index) do i=1,landuse_subgrid_dim(x_dim_nc_index) !Project the centre position to lat lon call PROJ2LL(x_landuse_subgrid(i,j),y_landuse_subgrid(i,j),temp_lon(1),temp_lat(1),projection_attributes,projection_type) !Project both sides to get the delta call PROJ2LL(x_landuse_subgrid(i,j)-landuse_subgrid_delta(x_dim_index)/2.,y_landuse_subgrid(i,j),temp_lon(2),temp_lat(2),projection_attributes,projection_type) call PROJ2LL(x_landuse_subgrid(i,j)+landuse_subgrid_delta(x_dim_index)/2.,y_landuse_subgrid(i,j),temp_lon(3),temp_lat(3),projection_attributes,projection_type) temp_delta(x_dim_index)=temp_lon(3)-temp_lon(2) call PROJ2LL(x_landuse_subgrid(i,j),y_landuse_subgrid(i,j)-landuse_subgrid_delta(y_dim_index)/2.,temp_lon(2),temp_lat(2),projection_attributes,projection_type) call PROJ2LL(x_landuse_subgrid(i,j),y_landuse_subgrid(i,j)+landuse_subgrid_delta(y_dim_index)/2.,temp_lon(3),temp_lat(3),projection_attributes,projection_type) temp_delta(y_dim_index)=temp_lat(3)-temp_lat(2) !Make a local correction to lon so it is essentially in the same units as lat so area averaging is correct correct_lon(1)=cos(3.14159/180.*temp_lat(1)) correct_lon(2)=1. !Take the nearest instead i_landuse_index=1+floor((temp_lon(1)-var2d_nc_dp(1,x_dim_nc_index))/delta_landuse_nc(1)+0.5) j_landuse_index=1+floor((temp_lat(1)-var2d_nc_dp(1,y_dim_nc_index))/delta_landuse_nc(2)+0.5) !write(*,*) i,j,temp_lon(1),temp_lat(1) landuse_subgrid(i,j,clc_index)=landuse_nc_dp(i_landuse_index,j_landuse_index) !Place the clc landuse in the EMEP landuse if (landuse_subgrid(i,j,clc_index).gt.0) then landuse_subgrid(i,j,Corine_to_EMEP_landuse(landuse_subgrid(i,j,clc_index)))=1 else landuse_subgrid(i,j,Corine_to_EMEP_landuse(NODATA_clc_value))=1 endif !Do the interpolation on the same grid then scale afterwards. Equivalent to interpolating density then rescaling with grid size !landuse_subgrid(i,j,clc_index)=area_weighted_extended_vectorgrid_interpolation_function( & ! real(var2d_nc_dp(1:dim_length_landuse_nc(x_dim_nc_index),x_dim_nc_index))*correct_lon(1),real(var2d_nc_dp(1:dim_length_landuse_nc(y_dim_nc_index),y_dim_nc_index)) & ! ,landuse_nc_dp(:,:,landuse_nc_index),dim_length_landuse_nc(x_dim_nc_index),dim_length_landuse_nc(y_dim_nc_index) & ! ,delta_landuse_nc*correct_lon,temp_lon(1)*correct_lon(1),temp_lat(1),delta_landuse_nc*correct_lon) !temp_scale=(temp_delta(1)*correct_lon(1)*temp_delta(2)*correct_lon(2))/(delta_landuse_nc(1)*correct_lon(1)*delta_landuse_nc(2)*correct_lon(2)) !write(*,*) temp_scale !landuse_subgrid(i,j,clc_index)=landuse_subgrid(i,j,clc_index)*temp_scale if (isnan(landuse_subgrid(i,j,clc_index))) then write(*,*) 'Stopping, nan in landuse_subgrid' write(*,*) temp_scale,correct_lon,delta_landuse_nc,temp_delta,temp_lon stop endif if (landuse_subgrid(i,j,clc_index).lt.0.) then write(*,*) 'Stopping, negative value in landuse_subgrid' write(*,*) temp_scale,correct_lon,delta_landuse_nc,temp_delta,temp_lon stop endif enddo enddo write(unit_logfile,'(A,2f12.2)') 'Max and min landuse in read domain: ',maxval(landuse_nc_dp(:,:)),minval(landuse_nc_dp(:,:)) write(unit_logfile,'(A,2f12.2)') 'Max and min landuse in subgrid domain: ',maxval(landuse_subgrid(:,:,clc_index)),minval(landuse_subgrid(:,:,clc_index)) if (use_landuse_as_proxy) then !Place the landuse as a proxy emission with the appropriate weights !loop through all sources and landuses do i_source=1,n_source_index !If source is to be downscaled and at least one landuse is selected then calculate the proxy emission weighting !write(*,*) i_source,calculate_source(i_source),sum(landuse_proxy_weighting(i_source,:)) if (calculate_source(i_source).and.sum(landuse_proxy_weighting(i_source,:)).gt.0) then proxy_emission_subgrid(:,:,i_source,:)=0. do i_landuse=1,n_clc_landuse_index !write(*,*) i_source,i_landuse,landuse_proxy_weighting(i_source,i_landuse) if (landuse_proxy_weighting(i_source,i_landuse).gt.0) then write(unit_logfile,'(A,i4,A,A,a,i4)') 'Distributing landuse index ',i_landuse,' to uEMEP sector "',trim(source_file_str(i_source)),'" and GNFR sector ',uEMEP_to_EMEP_sector(i_source) do j=1,emission_subgrid_dim(y_dim_nc_index,i_source) do i=1,emission_subgrid_dim(x_dim_nc_index,i_source) i_landuse_index=crossreference_emission_to_landuse_subgrid(i,j,x_dim_index,i_source) j_landuse_index=crossreference_emission_to_landuse_subgrid(i,j,y_dim_index,i_source) if (int(landuse_subgrid(i_landuse_index,j_landuse_index,clc_index)).eq.i_landuse) then proxy_emission_subgrid(i,j,i_source,:)=proxy_emission_subgrid(i,j,i_source,:)+landuse_proxy_weighting(i_source,i_landuse) !write(*,'(6i,2f12.2)') i,j,i_landuse_index,j_landuse_index,i_source,i_landuse,landuse_proxy_weighting(i_source,i_landuse),proxy_emission_subgrid(i,j,i_source,1) endif !If there is no data (0 or greater than the maximum number of landuse categories) then distribute emissions evenly on the EMEP grid if (int(landuse_subgrid(i_landuse_index,j_landuse_index,clc_index)).gt.n_clc_landuse_index.or.int(landuse_subgrid(i_landuse_index,j_landuse_index,clc_index)).lt.1) then proxy_emission_subgrid(i,j,i_source,:)=1. endif enddo enddo endif enddo endif enddo endif endif !No landuse available if (allocated(landuse_nc_dp)) deallocate (landuse_nc_dp) if (allocated(var2d_nc_dp)) deallocate (var2d_nc_dp) if (allocated(temp_var2d_nc_dp)) deallocate (temp_var2d_nc_dp) end subroutine uEMEP_read_netcdf_landuse_latlon