uEMEP_read_netcdf_population Subroutine

public subroutine uEMEP_read_netcdf_population()

Uses

  • proc~~uemep_read_netcdf_population~~UsesGraph proc~uemep_read_netcdf_population uEMEP_read_netcdf_population module~uemep_definitions uEMEP_definitions proc~uemep_read_netcdf_population->module~uemep_definitions netcdf netcdf proc~uemep_read_netcdf_population->netcdf

Arguments

None

Calls

proc~~uemep_read_netcdf_population~~CallsGraph proc~uemep_read_netcdf_population uEMEP_read_netcdf_population nf90_get_att nf90_get_att proc~uemep_read_netcdf_population->nf90_get_att nf90_get_var nf90_get_var proc~uemep_read_netcdf_population->nf90_get_var nf90_inq_dimid nf90_inq_dimid proc~uemep_read_netcdf_population->nf90_inq_dimid nf90_inq_varid nf90_inq_varid proc~uemep_read_netcdf_population->nf90_inq_varid nf90_inquire_dimension nf90_inquire_dimension proc~uemep_read_netcdf_population->nf90_inquire_dimension nf90_open nf90_open proc~uemep_read_netcdf_population->nf90_open proc~ll2laea LL2LAEA proc~uemep_read_netcdf_population->proc~ll2laea proc~ll2ltm ll2ltm proc~uemep_read_netcdf_population->proc~ll2ltm proc~ll2utm ll2utm proc~uemep_read_netcdf_population->proc~ll2utm dcos dcos proc~ll2ltm->dcos dsin dsin proc~ll2ltm->dsin dtan dtan proc~ll2ltm->dtan proc~ll2utm->dcos proc~ll2utm->dsin proc~ll2utm->dtan

Called by

proc~~uemep_read_netcdf_population~~CalledByGraph proc~uemep_read_netcdf_population uEMEP_read_netcdf_population program~uemep uEMEP program~uemep->proc~uemep_read_netcdf_population

Source Code

    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