query_som Subroutine

    subroutine query_som(kohonen_map,input_pattern,sample_index) !,output_patterns)
!========================================================================================
!!   Function to find the input samples associated with specific vector 
        class(self_organizing_map) :: kohonen_map
!!
        real(kind=wp),dimension(:,:),intent(inout) :: input_pattern
!!
        integer,allocatable :: sample_index(:)
!!
        integer :: ix,iy,iz,ihit,jhit,khit,ivar1,ivar2,nvar1,nvar2,number_patterns,ipat,ierr
        integer :: number_selected,i,pos
        real(kind=wp),dimension(kohonen_map%parameters%number_variables1,&
        kohonen_map%parameters%number_variables2) ::current_values
        real(kind=wp) :: dist,dist_min
        integer,dimension(size(kohonen_map%cells_index,1)) :: position,real_position
!
!(real_position(ix)=ix,ix=1,size(real_position))
        do ix=1,size(real_position)
            real_position(ix)=ix;
        enddo
        nvar1=kohonen_map%parameters%number_variables1;
        nvar2=kohonen_map%parameters%number_variables2;
        dist_min=1.0d10;
        !$OMP parallel do   
        do iz=1,size(kohonen_map%grid,3);
             do iy=1,size(kohonen_map%grid,2);
                 do ix=1,size(kohonen_map%grid,1);
                     dist=0.0_wp;
                     call kohonen_map%grid(ix,iy,iz)%get_prototype(current_values);
                     do ivar1=1,nvar1;
                         do ivar2=1,nvar2;
                             if(input_pattern(ivar1,ivar2) > 0.0_wp) then
                                 dist=dist+(input_pattern(ivar1,ivar2)-current_values(ivar1,ivar2))**2;
                             endif
                         enddo
                     enddo
                     if(dist < dist_min) then
                         dist_min=dist;
                         ihit=ix;jhit=iy;khit=iz;               
                     endif
                 enddo
             enddo
         enddo
         !$OMP end parallel do
!         write(*,*) 'BMU'
!         write(*,*) ihit,jhit,khit,dist_min
!
         position=0;
         number_patterns=kohonen_map%number_patterns(ihit,jhit,khit);
         if(number_patterns > 0) then
             where(kohonen_map%cells_index(:,1) == ihit .and. &
                 kohonen_map%cells_index(:,2) == jhit .and. &
                 kohonen_map%cells_index(:,3) == khit )
                 position=1;!real_position;
             end where
             number_selected=sum(position);
             pos=0
             if(number_selected > 0) then
                 allocate(sample_index(number_selected),stat=ierr);
                 do i=1,size(real_position)
                     if(position(i) == 1) then
                         pos=pos+1;
                         sample_index(pos)=real_position(i);
                         !write(*,*) 'Inside= ',i,real_position(i)
                     endif
                 enddo
             endif
             !write(*,*) kohonen_map%cells_index(118,1:3)
         else 
             write(*,*) 'WARNING: Query has returned an empty result'
            return
         endif
!
    end subroutine query_som