uEMEP_define_subgrid_extent Subroutine

public subroutine uEMEP_define_subgrid_extent()

Uses

  • proc~~uemep_define_subgrid_extent~~UsesGraph proc~uemep_define_subgrid_extent uEMEP_define_subgrid_extent module~uemep_definitions uEMEP_definitions proc~uemep_define_subgrid_extent->module~uemep_definitions

Arguments

None

Calls

proc~~uemep_define_subgrid_extent~~CallsGraph proc~uemep_define_subgrid_extent uEMEP_define_subgrid_extent proc~proj2ll PROJ2LL proc~uemep_define_subgrid_extent->proc~proj2ll 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 dcos dcos proc~ltm2ll->dcos dsin dsin proc~ltm2ll->dsin dtan dtan proc~ltm2ll->dtan proc~utm2ll->dcos proc~utm2ll->dsin proc~utm2ll->dtan

Called by

proc~~uemep_define_subgrid_extent~~CalledByGraph proc~uemep_define_subgrid_extent uEMEP_define_subgrid_extent program~uemep uEMEP program~uemep->proc~uemep_define_subgrid_extent

Source Code

    subroutine uEMEP_define_subgrid_extent

        use uEMEP_definitions

        implicit none

        integer i_source
        !integer ii,jj
        real dx_temp,dy_temp
        real lon_temp,lat_temp

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '================================================================'
        write(unit_logfile,'(A)') 'Define subgrid extent and buffer zones (uEMEP_define_subgrid_extent)'
        write(unit_logfile,'(A)') '================================================================'


        !If annual calculations then always set time start and stop to 1
        if (annual_calculations) then
            start_time_nc_index=1
            end_time_nc_index=1
            start_time_meteo_nc_index=1
            end_time_meteo_nc_index=1
        endif

        !Set the time index to be the same as the EMEP time dimensions
        emission_subgrid_dim(t_dim_index,:)=subgrid_dim(t_dim_index)
        init_emission_subgrid_dim(t_dim_index,:)=subgrid_dim(t_dim_index)
        integral_subgrid_dim(t_dim_index)=subgrid_dim(t_dim_index)

        write(unit_logfile,'(A,I6)')'Number of external time steps:',end_time_loop_index-start_time_loop_index+1
        write(unit_logfile,'(A,I6)')'Number of internal time steps:',subgrid_dim(t_dim_index)
        write(unit_logfile,'(A,2I6)')'Number of target grids:',subgrid_dim(1:2)
        write(unit_logfile,'(A,2I6)')'Number of integral grids:',integral_subgrid_dim(1:2)
        write(unit_logfile,'(A,2I6)')'Max number of emission grids:',emission_max_subgrid_dim(1:2)
        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                write(unit_logfile,'(A,A,2I6)')'Number of emission grids: ',trim(source_file_str(i_source)),emission_subgrid_dim(1:2,i_source)
                write(unit_logfile,'(A,A,2I8)')'Number of initial emission grids: ',trim(source_file_str(i_source)),init_emission_subgrid_dim(1:2,i_source)
            endif
        enddo


        !Allocate buffers and adjust the dimensions appropriately
        !Calculate the max loop size to cover the nearest EMEP grids. This avoids looping through all the grids
        !loop_index_scale=1.2*EMEP_grid_interpolation_size/2.*local_fraction_grid_size_scaling !Was 1.5
        loop_index_scale=1.2*EMEP_grid_interpolation_size/2.*local_fraction_grid_size_scaling !Was 1.5

        !Define the centre of the subgrid
        !ii=int(subgrid_dim(x_dim_index)/2)
        !jj=int(subgrid_dim(y_dim_index)/2)

        call PROJ2LL((subgrid_min(x_dim_index)+subgrid_max(x_dim_index))/2.,(subgrid_min(y_dim_index)+subgrid_max(y_dim_index))/2.,lon_temp,lat_temp,projection_attributes,projection_type)

        !if (projection_type.eq.RDM_projection_index) then
        !    call RDM2LL((subgrid_min(y_dim_index)+subgrid_max(y_dim_index))/2.,(subgrid_min(x_dim_index)+subgrid_max(x_dim_index))/2.,lat_temp,lon_temp)
        !elseif (projection_type.eq.UTM_projection_index) then
        !    call UTM2LL(utm_zone,(subgrid_min(y_dim_index)+subgrid_max(y_dim_index))/2.,(subgrid_min(x_dim_index)+subgrid_max(x_dim_index))/2.,lat_temp,lon_temp)
        !endif

        if (EMEP_projection_type.eq.LL_projection_index) then
            dx_temp=111000.*dgrid_nc(lon_nc_index)*cos(lat_temp*pi/180.)
            dy_temp=111000.*dgrid_nc(lat_nc_index)
        else
            !Assumed LCC or PS
            dx_temp=dgrid_nc(lon_nc_index)
            dy_temp=dgrid_nc(lat_nc_index)
        endif

        subgrid_loop_index(x_dim_index)=floor(dx_temp/subgrid_delta(x_dim_index)*loop_index_scale)
        subgrid_loop_index(y_dim_index)=floor(dy_temp/subgrid_delta(y_dim_index)*loop_index_scale)

        integral_subgrid_loop_index(x_dim_index)=floor(dx_temp/integral_subgrid_delta(x_dim_index)*loop_index_scale)
        integral_subgrid_loop_index(y_dim_index)=floor(dy_temp/integral_subgrid_delta(y_dim_index)*loop_index_scale)

        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                emission_subgrid_loop_index(x_dim_index,i_source)=floor(dx_temp/emission_subgrid_delta(x_dim_index,i_source)*loop_index_scale)
                emission_subgrid_loop_index(y_dim_index,i_source)=floor(dy_temp/emission_subgrid_delta(y_dim_index,i_source)*loop_index_scale)
            endif
        enddo

        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                init_emission_subgrid_loop_index(x_dim_index,i_source)=floor(dx_temp/init_emission_subgrid_delta(x_dim_index,i_source)*loop_index_scale)
                init_emission_subgrid_loop_index(y_dim_index,i_source)=floor(dy_temp/init_emission_subgrid_delta(y_dim_index,i_source)*loop_index_scale)
            endif
        enddo

        if (calculate_deposition_flag) then
            deposition_subgrid_loop_index(x_dim_index)=floor(dx_temp/deposition_subgrid_delta(x_dim_index)*loop_index_scale)
            deposition_subgrid_loop_index(y_dim_index)=floor(dy_temp/deposition_subgrid_delta(y_dim_index)*loop_index_scale)
        endif
        if (read_landuse_flag) then
            landuse_subgrid_loop_index(x_dim_index)=floor(dx_temp/landuse_subgrid_delta(x_dim_index)*loop_index_scale)
            landuse_subgrid_loop_index(y_dim_index)=floor(dy_temp/landuse_subgrid_delta(y_dim_index)*loop_index_scale)
        endif

        !Set the buffer sizes according to these loops for emissions only
        !This will remove edge effects for dispersion but will only remove edge effects for moving window only when emissions are used for redistribution
        buffer_index_scale=loop_index_scale
        !if (use_downwind_position_flag) buffer_index_scale=buffer_index_scale*2.

        if (use_buffer_zone) then
            buffer_index(x_dim_index)=floor(dx_temp/subgrid_delta(x_dim_index)*buffer_index_scale)
            buffer_index(y_dim_index)=floor(dy_temp/subgrid_delta(y_dim_index)*buffer_index_scale)

            integral_buffer_index(x_dim_index)=floor(dx_temp/integral_subgrid_delta(x_dim_index)*buffer_index_scale)
            integral_buffer_index(y_dim_index)=floor(dy_temp/integral_subgrid_delta(y_dim_index)*buffer_index_scale)

            !Buffer index for emissions a half grid larger because of the possible use of the integral
            do i_source=1,n_source_index
                if (calculate_source(i_source)) then
                    emission_buffer_index(x_dim_index,i_source)=floor(dx_temp/emission_subgrid_delta(x_dim_index,i_source)*(buffer_index_scale+0.5))
                    emission_buffer_index(y_dim_index,i_source)=floor(dy_temp/emission_subgrid_delta(y_dim_index,i_source)*(buffer_index_scale+0.5))
                    if (local_subgrid_method_flag.eq.3) then
                        !Extend the emission buffer zone 0.5 EMEP grids when redistributing emissions to gurantee all emissions are accounted for
                        emission_buffer_index(x_dim_index,i_source)=floor(dx_temp/emission_subgrid_delta(x_dim_index,i_source)*(buffer_index_scale+1.0))
                        emission_buffer_index(y_dim_index,i_source)=floor(dy_temp/emission_subgrid_delta(y_dim_index,i_source)*(buffer_index_scale+1.0))
                    endif
                endif
            enddo

            do i_source=1,n_source_index
                if (calculate_source(i_source)) then
                    init_emission_buffer_index(x_dim_index,i_source)=floor(dx_temp/init_emission_subgrid_delta(x_dim_index,i_source)*(buffer_index_scale+0.5))
                    init_emission_buffer_index(y_dim_index,i_source)=floor(dy_temp/init_emission_subgrid_delta(y_dim_index,i_source)*(buffer_index_scale+0.5))
                    if (local_subgrid_method_flag.eq.3) then
                        !Extend the emission buffer zone 0.5 EMEP grids when redistributing emissions to gurantee all emissions are accounted for
                        init_emission_buffer_index(x_dim_index,i_source)=floor(dx_temp/init_emission_subgrid_delta(x_dim_index,i_source)*(buffer_index_scale+1.0))
                        init_emission_buffer_index(y_dim_index,i_source)=floor(dy_temp/init_emission_subgrid_delta(y_dim_index,i_source)*(buffer_index_scale+1.0))
                    endif
                endif
            enddo

            if (calculate_deposition_flag) then
                deposition_buffer_index(x_dim_index)=floor(dx_temp/deposition_subgrid_delta(x_dim_index)*buffer_index_scale)
                deposition_buffer_index(y_dim_index)=floor(dy_temp/deposition_subgrid_delta(y_dim_index)*buffer_index_scale)
            endif
            !Make sure the landuse bufferzone is the same as the emissions in case it is used as a proxy
            if (read_landuse_flag) then
                landuse_buffer_index(x_dim_index)=floor(dx_temp/landuse_subgrid_delta(x_dim_index)*(buffer_index_scale+0.5))
                landuse_buffer_index(y_dim_index)=floor(dy_temp/landuse_subgrid_delta(y_dim_index)*(buffer_index_scale+0.5))
                if (local_subgrid_method_flag.eq.3) then
                    landuse_buffer_index(x_dim_index)=floor(dx_temp/landuse_subgrid_delta(x_dim_index)*(buffer_index_scale+1.0))
                    landuse_buffer_index(y_dim_index)=floor(dy_temp/landuse_subgrid_delta(y_dim_index)*(buffer_index_scale+1.0))
                endif

            endif
        else
            buffer_index=0
            emission_buffer_index=0
            init_emission_buffer_index=0
            integral_buffer_index=0
            deposition_buffer_index=0
            landuse_buffer_index=0
        endif
        buffer_size=buffer_index*subgrid_delta
        emission_buffer_size=emission_buffer_index*emission_subgrid_delta
        init_emission_buffer_size=emission_buffer_index*init_emission_subgrid_delta
        integral_buffer_size=integral_buffer_index*integral_subgrid_delta
        if (calculate_deposition_flag) then
            deposition_buffer_size=deposition_buffer_index*deposition_subgrid_delta
        endif
        if (read_landuse_flag) then
            landuse_buffer_size=landuse_buffer_index*landuse_subgrid_delta
        endif

        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                emission_subgrid_dim(1:2,i_source)=emission_subgrid_dim(1:2,i_source)+emission_buffer_index(1:2,i_source)*2
                emission_subgrid_min(1:2,i_source)=emission_subgrid_min(1:2,i_source)-emission_buffer_size(1:2,i_source)
                emission_subgrid_max(1:2,i_source)=emission_subgrid_max(1:2,i_source)+emission_buffer_size(1:2,i_source)
            endif
        enddo
        !emission_max_subgrid_dim(1:2)=emission_max_subgrid_dim(1:2)+buffer_index(1:2)*2
        emission_max_subgrid_dim(1:2)=maxval(emission_subgrid_dim(1:2,:),2)

        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                init_emission_subgrid_dim(1:2,i_source)=init_emission_subgrid_dim(1:2,i_source)+init_emission_buffer_index(1:2,i_source)*2
                init_emission_subgrid_min(1:2,i_source)=init_emission_subgrid_min(1:2,i_source)-init_emission_buffer_size(1:2,i_source)
                init_emission_subgrid_max(1:2,i_source)=init_emission_subgrid_max(1:2,i_source)+init_emission_buffer_size(1:2,i_source)
            endif
        enddo

        integral_subgrid_dim(1:2)=integral_subgrid_dim(1:2)+integral_buffer_index(1:2)*2
        integral_subgrid_min(1:2)=integral_subgrid_min(1:2)-integral_buffer_size(1:2)
        integral_subgrid_max(1:2)=integral_subgrid_max(1:2)+integral_buffer_size(1:2)

        if (calculate_deposition_flag) then
            deposition_subgrid_dim(1:2)=deposition_subgrid_dim(1:2)+deposition_buffer_index(1:2)*2
            deposition_subgrid_min(1:2)=deposition_subgrid_min(1:2)-deposition_buffer_size(1:2)
            deposition_subgrid_max(1:2)=deposition_subgrid_max(1:2)+deposition_buffer_size(1:2)
        endif
        if (read_landuse_flag) then
            landuse_subgrid_dim(1:2)=landuse_subgrid_dim(1:2)+landuse_buffer_index(1:2)*2
            landuse_subgrid_min(1:2)=landuse_subgrid_min(1:2)-landuse_buffer_size(1:2)
            landuse_subgrid_max(1:2)=landuse_subgrid_max(1:2)+landuse_buffer_size(1:2)
        endif

        write(unit_logfile,'(A,2I5)')'Number of target grids to be looped for each EMEP grid:',subgrid_loop_index(1:2)
        write(unit_logfile,'(A,2I5)')'Number of integral grids to be looped for each EMEP grid:',integral_subgrid_loop_index(1:2)
        write(unit_logfile,'(A,2I5)')'Size of integral grid buffer zone:',integral_buffer_index(1:2)
        if (calculate_deposition_flag) then
            write(unit_logfile,'(A,2I5)')'Size of deposition grid buffer zone:',deposition_buffer_index(1:2)
        endif
        if (read_landuse_flag) then
            write(unit_logfile,'(A,2I5)')'Size of landuse grid buffer zone:',landuse_buffer_index(1:2)
        endif

        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                write(unit_logfile,'(A,A,2I5)')'Size of emission grid buffer zone: ',trim(source_file_str(i_source)),emission_buffer_index(1:2,i_source)
                write(unit_logfile,'(A,A,2I5)')'Size of initial emission grid buffer zone: ',trim(source_file_str(i_source)),emission_buffer_index(1:2,i_source)
            endif
        enddo

        write(unit_logfile,'(A,2I6)')'Number of target grids:',subgrid_dim(1:2)
        write(unit_logfile,'(A,2I6)')'Number of integral grids:',integral_subgrid_dim(1:2)
        write(unit_logfile,'(A,2I6)')'Number of integral grids:',integral_subgrid_dim(1:2)
        if (calculate_deposition_flag) then
            write(unit_logfile,'(A,2I6)')'Number of deposition grids:',deposition_subgrid_dim(1:2)
        endif
        write(unit_logfile,'(A,2I6)')'Max number of emission grids:',emission_max_subgrid_dim(1:2)
        do i_source=1,n_source_index
            if (calculate_source(i_source)) then
                write(unit_logfile,'(A,A,2I6)')   'Number of emission grids:',trim(source_file_str(i_source)),emission_subgrid_dim(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Min of emission grids:   ',trim(source_file_str(i_source)),emission_subgrid_min(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Max of emission grids:   ',trim(source_file_str(i_source)),emission_subgrid_max(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Delta of emission grids: ',trim(source_file_str(i_source)),emission_subgrid_delta(1:2,i_source)
            endif
            if (calculate_source(i_source)) then
                write(unit_logfile,'(A,A,2I8)')   'Number of initial emission grids:',trim(source_file_str(i_source)),init_emission_subgrid_dim(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Min of initial emission grids:   ',trim(source_file_str(i_source)),init_emission_subgrid_min(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Max of initial emission grids:   ',trim(source_file_str(i_source)),init_emission_subgrid_max(1:2,i_source)
                write(unit_logfile,'(A,A,2f12.1)')'Delta of initial emission grids: ',trim(source_file_str(i_source)),init_emission_subgrid_delta(1:2,i_source)
            endif
        enddo

    end subroutine uEMEP_define_subgrid_extent