MOHID snippet from the ModuleFunctions (See related posts)

    subroutine ExtraPol3DNearestCell (ILB, IUB, JLB, JUB, KLB, KUB, ComputePoints3D, OutValues3D)

        !Arguments-------------------------------------------------------------
        integer                                     :: JLB, JUB, ILB, IUB, KLB, KUB  
        real,         dimension(:,:,:), pointer     :: OutValues3D
        integer,      dimension(:,:,:), pointer     :: ComputePoints3D

        !Local-----------------------------------------------------------------
        integer                                     :: dij, dk, Count, i, j, k, NumberOfCells
        integer                                     :: jj, ii, kk, dijmax, dimax, djmax
        real                                        :: SumValues
        logical                                     :: NoMapping, OkMap
        
        !Begin-----------------------------------------------------------------

d1:     do k = KLB, KUB

            if (associated(ComputePoints3D)) then

                NumberOfCells =  Sum(ComputePoints3D(ILB:IUB, JLB:JUB, k))
                NoMapping     = .false.
            
            else
            
                NoMapping     = .true. 
                NumberOfCells = 1
            
            endif 

            if (NumberOfCells > 0) then

                do j = JLB, JUB
                do i = ILB, IUB
                
                    if (NoMapping) then
                    
                        OkMap = .true.
                        
                                   
                    else
                    
                        if (ComputePoints3D(i, j, k) == 1) then
                            OkMap = .true.
                        else
                            OkMap = .false.
                        endif    
                    
                    endif

                    if (OutValues3D(i, j, k) < FillValueReal/4. .and. OkMap) then

                        dimax = IUB-ILB + 1
                        djmax = JUB-JLB + 1

                        dijmax = max(dimax, djmax)
                    
                        SumValues   = 0
                        Count = 0

                        do dij=1,dijmax

                            do jj=j-dij,j+dij
                            do ii=i-dij,i+dij

                                if (jj < JLB) cycle
                                if (jj > JUB) cycle
                                if (ii < ILB) cycle
                                if (ii > IUB) cycle

                                if (OutValues3D(ii, jj, k) > FillValueReal/4.) then
                                    SumValues   = SumValues   + OutValues3D(ii, jj, k) 
                                    Count = Count + 1
                                endif

                            enddo
                            enddo

                            if (Count > 0) exit

                        enddo

                        if (Count > 0) then

                            OutValues3D(i, j, k) = SumValues / real(Count)

                        else
                            do dk = 1, KUB - KLB + 1
                            
                                do kk = k - dk, k + dk

                                    if (kk < KLB) cycle
                                    if (kk > KUB) cycle
                                    
                                    if (OutValues3D(i, j, kk) > FillValueReal/4.) then
                                        SumValues   = SumValues   + OutValues3D(i, j, kk) 
                                        Count = Count + 1
                                    endif

                                enddo  
                                
                                if (Count >0) exit
                                
                            enddo 
                                
                            if (Count > 0) then

                                OutValues3D(i, j, k) = SumValues / real(Count)
                            
                            else

                                if (NoMapping) then
                                
                                    OutValues3D(i, j, k) = FillValueReal
                                                                       
                                else

                                    stop 'ExtraPol3DNearestCell - ModuleFunctions - ERR10'
                                
                                endif
                                
                            endif
                        endif

                    endif

                enddo
                enddo

            endif

        enddo d1

    end subroutine ExtraPol3DNearestCell