! $Id$
!
!
program test_it
   real :: data1(10,10)
   real :: data2(10,10)
   real :: minv, maxv
   logical :: suc
   
   
   
   minv = 2
   maxv = 6
   data1(1,:) = [7.,5.,3.,4.,0.,3.,2.,3.,3.,21.]
   data1(2,:) = [9.,5.,3.,4.,2.,3.,2.,3.,3.,24.]
   data1(3,:) = [8.,5.,3.,4.,2.,3.,2.,3.,3.,22.]
   data1(4,:) = [5.,5.,3.,1.,1.,3.,2.,3.,3.,21.]
   data1(5,:) = [7.,5.,3.,1.,0.,3.,2.,3.,3.,21.]
   data1(6,:) = [7.,5.,3.,2.,1.,3.,2.,3.,3.,23.]
   data1(7,:) = [9.,5.,3.,1.,2.,3.,2.,3.,3.,21.]
   data1(8,:) = [8.,7.,3.,4.,2.,3.,2.,3.,3.,22.]
   data1(9,:) = [8.,6.,3.,4.,2.,3.,2.,3.,3.,21.]
   data1(10,:) =[4.,7.,3.,4.,2.,3.,2.,3.,3.,21.]
   suc = local_radiative_center ( data1,minv,maxv, data2)
   
   
   contains
   

function local_radiative_center ( data_array , min_val , max_val &   
           , dum,  mask_in  ) result (suc)

   real, intent(in) :: data_array(:,:)
   real, intent(in) :: min_val , max_val
   logical, optional , intent (in) :: mask_in (:,:)
   real,  intent(out)  :: dum(:,:)
   logical :: suc
   integer, parameter :: steps_max = 10
   
   integer :: dim_1 , dim_2
   real, dimension(3,3):: local_box , local_box_grad
   integer :: idx_1 , idx_2
   integer :: idx_1_step, idx_2_step
   integer :: start_1 , start_2
   integer :: end_1 , end_2
   integer :: high_grad(2)
   
   
   suc = .false.
   dim_1 = size(data_array,1)
   dim_2 = size(data_array,2)
  
   dum = -999. 
   start_1 = 1
   start_2 = 1
   end_1 = dim_1
   end_2 = dim_2
   
   
   do idx_1 = start_1 + 1 , end_1 - 1
      do idx_2 = start_1 + 1 , end_2 - 1
      
         idx_1_step = idx_1
         idx_2_step = idx_2
         
         do idx_step = 1 ,  steps_max  
            local_box = data_array ( idx_1_step - 1 : idx_1_step + 1 &
                                & ,  idx_2_step - 1 : idx_2_step + 1  )   
            
            ! - if lowest value cycle
            if (minval(local_box) == data_array (idx_1_step,idx_2_step ) ) then
               dum(idx_1,idx_2) = data_array ( idx_1_step, idx_2_step)
               cycle 
            end if
            
            ! - find the next  pixel
            local_box_grad =  local_box - data_array (idx_1_previous,idx_2_previous )
            high_grad =  minloc ( local_box_grad )
            idx_1_step = idx_1_step + high_grad(1) - 2
            idx_2_step = idx_2_step + high_grad(2) - 2
            
            if ( data_array ( idx_1_step, idx_2_step ) < min_val ) then
               dum(idx_1,idx_2) = data_array ( idx_1_step, idx_2_step)
               cycle 
            end if 
           
            
         end do
      end do
   end do   
   
  
   suc = .true.
   
   
   return

end function

end program