read_country_bounding_box_data Subroutine

public subroutine read_country_bounding_box_data()

Uses

  • proc~~read_country_bounding_box_data~~UsesGraph proc~read_country_bounding_box_data read_country_bounding_box_data module~uemep_definitions uEMEP_definitions proc~read_country_bounding_box_data->module~uemep_definitions

Arguments

None

Calls

proc~~read_country_bounding_box_data~~CallsGraph proc~read_country_bounding_box_data read_country_bounding_box_data proc~ll2laea LL2LAEA proc~read_country_bounding_box_data->proc~ll2laea proc~ll2ltm ll2ltm proc~read_country_bounding_box_data->proc~ll2ltm proc~ll2utm ll2utm proc~read_country_bounding_box_data->proc~ll2utm proc~proj2ll PROJ2LL proc~read_country_bounding_box_data->proc~proj2ll dcos dcos proc~ll2ltm->dcos dsin dsin proc~ll2ltm->dsin dtan dtan proc~ll2ltm->dtan proc~ll2utm->dcos proc~ll2utm->dsin proc~ll2utm->dtan proc~laea2ll LAEA2LL proc~proj2ll->proc~laea2ll proc~lambert2lb2_uemep lambert2lb2_uEMEP proc~proj2ll->proc~lambert2lb2_uemep proc~ltm2ll ltm2ll proc~proj2ll->proc~ltm2ll proc~ps2ll_spherical PS2LL_spherical proc~proj2ll->proc~ps2ll_spherical proc~rdm2ll RDM2LL proc~proj2ll->proc~rdm2ll proc~utm2ll utm2ll proc~proj2ll->proc~utm2ll proc~ltm2ll->dcos proc~ltm2ll->dsin proc~ltm2ll->dtan proc~utm2ll->dcos proc~utm2ll->dsin proc~utm2ll->dtan

Called by

proc~~read_country_bounding_box_data~~CalledByGraph proc~read_country_bounding_box_data read_country_bounding_box_data program~uemep uEMEP program~uemep->proc~read_country_bounding_box_data

Source Code

    subroutine read_country_bounding_box_data
        !This routine reads in a file that provides information on the country bounding box
        !as well as reference to the filename used in the OSM files
        !Is also useful for other purposes but used here only for OSM file names

        use uEMEP_definitions

        implicit none

        logical :: exists
        integer i
        character(256) CNTR_ID
        character(256) OSM_country,Long_name
        real min_lat,min_lon,max_lat,max_lon,min_y_3035,min_x_3035,max_y_3035,max_x_3035
        real lon_grid_min,lat_grid_min,lon_grid_max,lat_grid_max
        real lon_grid_min2,lat_grid_min2,lon_grid_max2,lat_grid_max2
        real lon_new_min,lat_new_min,lon_new_max,lat_new_max
        integer unit_in
        character(256) temp_str
        integer count
        logical found_country
        real x_out(4),y_out(4)
        integer :: io

        !Will fail at lon=-180

        !Read in replacement file
        pathfilename_boundingbox=trim(pathname_boundingbox)//trim(filename_boundingbox)

        !Test existence of the road link filename (2). If does not exist then use default
        inquire(file=trim(pathfilename_boundingbox),exist=exists)
        if (.not.exists) then
            write(unit_logfile,'(A,A)') ' ERROR: Bounding box file does not exist: ', trim(pathfilename_boundingbox)
            stop
        endif

        if (trim(select_country_by_name).ne.'') then

            !Open the file for reading
            unit_in=20
            open(unit_in,file=pathfilename_boundingbox,access='sequential',status='old',readonly)
            write(unit_logfile,'(a)') ' Opening Bounding box file  '//trim(filename_boundingbox)
            rewind(unit_in)

            !Skip over the header
            !Index	CNTR_ID	OSM_country	min_lat	min_lon	max_lat	max_lon	min_y_3035	min_x_3035	max_y_3035	max_x_3035		Long_name
            read(unit_in,*) temp_str

            !Read coordinates
            found_country=.false.
            do
                read(unit_in,*,iostat=io) i,CNTR_ID,OSM_country,min_lon,min_lat,max_lon,max_lat,min_x_3035,min_y_3035,max_x_3035,max_y_3035,Long_name
                if (io < 0) then
                    exit
                else if (io > 0) then
                    write(unit_logfile,"(2a)") "ERROR reading file: ", trim(pathfilename_boundingbox)
                    stop 1
                end if
                !write(*,*) i,trim(CNTR_ID),trim(OSM_country),min_lon,min_lat,max_lon,max_lat,min_x_3035,min_y_3035,max_x_3035,max_y_3035,trim(Long_name)

                if (index(trim(select_country_by_name),trim(CNTR_ID)).gt.0) then

                    !Set the min and max lat and lon values for the current grid
                    if  (projection_type.eq.UTM_projection_index) then
                        call ll2utm(1,utm_zone,min_lat,min_lon,y_out(1),x_out(1))
                        call ll2utm(1,utm_zone,max_lat,max_lon,y_out(2),x_out(2))
                        call ll2utm(1,utm_zone,max_lat,min_lon,y_out(3),x_out(3))
                        call ll2utm(1,utm_zone,min_lat,max_lon,y_out(4),x_out(4))
                    elseif  (projection_type.eq.LTM_projection_index) then
                        call ll2ltm(1,ltm_lon0,min_lat,min_lon,y_out(1),x_out(1))
                        call ll2ltm(1,ltm_lon0,max_lat,max_lon,y_out(2),x_out(2))
                        call ll2ltm(1,ltm_lon0,max_lat,min_lon,y_out(3),x_out(3))
                        call ll2ltm(1,ltm_lon0,min_lat,max_lon,y_out(4),x_out(4))
                    elseif (projection_type.eq.LAEA_projection_index) then
                        call LL2LAEA(x_out(1),y_out(1),min_lon,min_lat,projection_attributes)
                        call LL2LAEA(x_out(2),y_out(2),max_lon,max_lat,projection_attributes)
                        call LL2LAEA(x_out(3),y_out(3),min_lon,max_lat,projection_attributes)
                        call LL2LAEA(x_out(4),y_out(4),max_lon,min_lat,projection_attributes)
                    endif

                    subgrid_min(x_dim_index)=minval(x_out)
                    subgrid_max(x_dim_index)=maxval(x_out)
                    subgrid_min(y_dim_index)=minval(y_out)
                    subgrid_max(y_dim_index)=maxval(y_out)

                    !Snap to nearest 10 km
                    subgrid_min(x_dim_index)=floor(subgrid_min(x_dim_index)/10000.)*10000
                    subgrid_min(y_dim_index)=floor(subgrid_min(y_dim_index)/10000.)*10000
                    subgrid_max(x_dim_index)=ceiling(subgrid_max(x_dim_index)/10000.)*10000
                    subgrid_max(y_dim_index)=ceiling(subgrid_max(y_dim_index)/10000.)*10000

                    write(unit_logfile,'(i,a)') i,' Setting grid for ID: '//trim(CNTR_ID)//'   OSM name: '//trim(OSM_country)//'   Name: '//trim(Long_name)
                    write(unit_logfile,'(a,f12.1)') 'subgrid_min(x_dim_index)=',subgrid_min(x_dim_index)
                    write(unit_logfile,'(a,f12.1)') 'subgrid_min(y_dim_index)=',subgrid_min(y_dim_index)
                    write(unit_logfile,'(a,f12.1)') 'subgrid_max(x_dim_index)=',subgrid_max(x_dim_index)
                    write(unit_logfile,'(a,f12.1)') 'subgrid_max(y_dim_index)=',subgrid_max(y_dim_index)
                    found_country=.true.
                endif

                !Reset the initial subgrid as well, needed for EMEP and receptor selection
                init_subgrid_min(x_dim_index)=subgrid_min(x_dim_index)
                init_subgrid_min(y_dim_index)=subgrid_min(y_dim_index)
                init_subgrid_max(x_dim_index)=subgrid_max(x_dim_index)
                init_subgrid_max(y_dim_index)=subgrid_max(y_dim_index)

            enddo

            close(unit_in)

            if (.not.found_country) then
                write(unit_logfile,'(a)') ' No country with this ID found: '//trim(select_country_by_name)
            endif

        endif

        !Open the file for reading

        if (auto_select_OSM_country_flag) then

            !Set the min and max lat and lon values for the current grid

            call PROJ2LL(subgrid_min(x_dim_index),subgrid_min(y_dim_index),lon_grid_min,lat_grid_min,projection_attributes,projection_type)
            call PROJ2LL(subgrid_max(x_dim_index),subgrid_max(y_dim_index),lon_grid_max,lat_grid_max,projection_attributes,projection_type)
            call PROJ2LL(subgrid_min(x_dim_index),subgrid_max(y_dim_index),lon_grid_min2,lat_grid_max2,projection_attributes,projection_type)
            call PROJ2LL(subgrid_max(x_dim_index),subgrid_min(y_dim_index),lon_grid_max2,lat_grid_min2,projection_attributes,projection_type)
            lon_grid_max=max(lon_grid_max,lon_grid_max2)
            lon_grid_min=min(lon_grid_min,lon_grid_min2)
            lat_grid_max=max(lat_grid_max,lat_grid_max2)
            lat_grid_min=min(lat_grid_min,lat_grid_min2)

            !write(*,*) projection_attributes
            !write(*,*) projection_type

            !Test
            !lon_grid_max=5.;lat_grid_max=50.
            !call LL2LAEA(subgrid_max(x_dim_index),subgrid_max(y_dim_index),lon_grid_max,lat_grid_max,projection_attributes,projection_type)
            !call LAEA2LL(subgrid_max(x_dim_index),subgrid_max(y_dim_index),lon_grid_max,lat_grid_max,projection_attributes,projection_type)

            write(unit_logfile,'(a)') ' Opening Bounding box file  '//trim(filename_boundingbox)
            write(unit_logfile,'(a,2f12.6)') ' Lon (min,max)  ',lon_grid_min,lon_grid_max
            write(unit_logfile,'(a,2f12.6)') ' Lat (min,max)  ',lat_grid_min,lat_grid_max


            unit_in=20
            open(unit_in,file=pathfilename_boundingbox,access='sequential',status='old',readonly)
            rewind(unit_in)

            !Skip over the header
            !Index	CNTR_ID	OSM_country	min_lat	min_lon	max_lat	max_lon	min_y_3035	min_x_3035	max_y_3035	max_x_3035		Long_name
            read(unit_in,*) temp_str

            !Read coordinates
            count=0
            filename_mrl=''
            do
                read(unit_in,*,iostat=io) i,CNTR_ID,OSM_country,min_lon,min_lat,max_lon,max_lat,min_x_3035,min_y_3035,max_x_3035,max_y_3035,Long_name
                if (io /= 0) exit
                !write(*,*) i,trim(CNTR_ID),trim(OSM_country),min_lon,min_lat,max_lon,max_lat,min_x_3035,min_y_3035,max_x_3035,max_y_3035,trim(Long_name)

                !test the bounding box in lat lon coordinates
                lon_new_min=max(min_lon,lon_grid_min)
                lat_new_min=max(min_lat,lat_grid_min)
                lon_new_max=min(max_lon,lon_grid_max)
                lat_new_max=min(max_lat,lat_grid_max)

                if (lon_new_min.gt.lon_new_max.or.lat_new_min.gt.lat_new_max) then
                    !No intersection so do nothing
                elseif (index('none',trim(OSM_country)).le.0) then
                    !update and attribute the filename
                    count=count+1
                    if (count.gt.50) then
                        write(unit_logfile,'(a)') ' Max files are 50. Stopping  '
                        stop
                    endif

                    filename_mrl(count)='Road_data_OSM_'//trim(OSM_country)//'_latlon.txt'
                    write(unit_logfile,'(2i,a)') count,i, ' Including OSM file:  '//trim(filename_mrl(count))
                endif

            enddo

            if (count.eq.0) then
                write(unit_logfile,'(a)') ' No countries overlap this area. Setting OSM to default file andorra'
                OSM_country='andorra';count=1
                filename_mrl(count)='Road_data_OSM_'//trim(OSM_country)//'_latlon.txt'
                num_multiple_roadlink_files=count
            else
                write(unit_logfile,'(a,i)') ' Specifying this many OSM road link files to be read',count
                num_multiple_roadlink_files=count
            endif


            close(unit_in)


        endif

        !stop


    end subroutine read_country_bounding_box_data