area_weighted_extended_interpolation_function Function

public function area_weighted_extended_interpolation_function(xgrid, ygrid, zgrid, xdim, ydim, delta, xval, yval, delta_val) result(res)

Returns the area weighted value for rectangle of size delta_val at position xval, yval from the grid values xgrid,ygrid,zgrid Delta_val can be any size

Arguments

Type IntentOptional Attributes Name
real, intent(in) :: xgrid(xdim,ydim)
real, intent(in) :: ygrid(xdim,ydim)
real, intent(in) :: zgrid(xdim,ydim)
integer, intent(in) :: xdim
integer, intent(in) :: ydim
real, intent(in) :: delta(2)
real, intent(in) :: xval
real, intent(in) :: yval
real, intent(in) :: delta_val(2)

Return Value real


Called by

proc~~area_weighted_extended_interpolation_function~~CalledByGraph proc~area_weighted_extended_interpolation_function area_weighted_extended_interpolation_function proc~uemep_read_agriculture_rivm_data uEMEP_read_agriculture_rivm_data proc~uemep_read_agriculture_rivm_data->proc~area_weighted_extended_interpolation_function proc~uemep_subgrid_deposition uEMEP_subgrid_deposition proc~uemep_subgrid_deposition->proc~area_weighted_extended_interpolation_function proc~uemep_subgrid_deposition_emep uEMEP_subgrid_deposition_EMEP proc~uemep_subgrid_deposition_emep->proc~area_weighted_extended_interpolation_function proc~uemep_calculate_emissions_for_emep uEMEP_calculate_emissions_for_EMEP proc~uemep_calculate_emissions_for_emep->proc~uemep_read_agriculture_rivm_data program~uemep uEMEP program~uemep->proc~uemep_read_agriculture_rivm_data program~uemep->proc~uemep_subgrid_deposition program~uemep->proc~uemep_subgrid_deposition_emep program~uemep->proc~uemep_calculate_emissions_for_emep

Source Code

    function area_weighted_extended_interpolation_function(xgrid, ygrid, zgrid, xdim, ydim, delta, xval, yval, delta_val) result(res)
        !! Returns the area weighted value for rectangle of size delta_val at position xval, yval from the grid values xgrid,ygrid,zgrid
        !! Delta_val can be any size
        integer, intent(in) :: xdim, ydim
        real, intent(in) :: delta(2)
        real, intent(in) :: delta_val(2)
        real, intent(in) :: xgrid(xdim,ydim), ygrid(xdim,ydim), zgrid(xdim,ydim)
        real, intent(in) :: xval, yval
        real :: res

        ! Local variables
        real :: zval, sum_weight, weighting
        real :: xpos_max, xpos_min, ypos_max, ypos_min
        real :: xpos_area_max, xpos_area_min, ypos_area_max, ypos_area_min
        integer :: i, j, ii, jj, iii, jjj
        integer :: ii_delta, jj_delta

        ! If only on grid available then return the value of that grid
        if (xdim .eq. 1 .and. ydim .eq. 1) then
            res = zgrid(xdim,ydim)
            return
        endif

        ! Find grid index for position val
        i = 1 + floor((xval - xgrid(1,1))/delta(1) + 0.5)
        j = 1 + floor((yval - ygrid(1,1))/delta(2) + 0.5)
        i = max(1,i); i = min(xdim,i)
        j = max(1,j); j = min(ydim,j)

        if (i .lt. 1 .or. j .lt. 1 .or. i .gt. xdim .or. j .gt. ydim) then
            write(*,'(A,4i6)') 'Interpolation out of range in area_weighted_extended_interpolation_function. Stopping. (i,j,xdim,ydim)', i, j, xdim, ydim
            write(*,'(4f12.2)') xval, yval, xgrid(1,1), ygrid(1,1)
            stop
        else
            xpos_area_max = xval + delta_val(1)/2.0
            xpos_area_min = xval - delta_val(1)/2.0
            ypos_area_max = yval + delta_val(2)/2.0
            ypos_area_min = yval - delta_val(2)/2.0

            jj_delta = 1 + floor(0.5*(delta_val(2)/delta(2) - 1.0))
            ii_delta = 1 + floor(0.5*(delta_val(1)/delta(1) - 1.0))

            zval = 0.0
            sum_weight = 0.0
            do jjj = j-jj_delta, j+jj_delta
                do iii = i-ii_delta, i+ii_delta
                    jj = max(jjj,1); jj = min(jj,ydim)
                    ii = max(iii,1); ii = min(ii,xdim)
                    xpos_min = max(xpos_area_min, xgrid(ii,jj) - delta(1)/2.0)
                    xpos_max = min(xpos_area_max, xgrid(ii,jj) + delta(1)/2.0)
                    ypos_min = max(ypos_area_min, ygrid(ii,jj) - delta(2)/2.0)
                    ypos_max = min(ypos_area_max, ygrid(ii,jj) + delta(2)/2.0)

                    if (xpos_max .gt. xpos_min .and. ypos_max .gt. ypos_min) then
                        weighting = (ypos_max - ypos_min)*(xpos_max - xpos_min)/delta_val(1)/delta_val(2)
                    else
                        weighting = 0.0
                    endif
                    zval = zval + zgrid(ii,jj)*weighting
                    sum_weight = sum_weight + weighting
                enddo
            enddo
        endif
        res = zval
    end function area_weighted_extended_interpolation_function