read_som_layer Subroutine

   subroutine read_som_layer(kohonen_map,som_fl,layer_type)
   !========================================================================================
!! Subroutine to read the prototypes of the first/seconf layer of a two level self_organized_map 
      class(two_level_self_organizing_map) :: kohonen_map
!! A `two_level_self_organizing_map` object
      character(len=*) :: som_fl,layer_type
!! A character variable
      logical :: testfl,toroidal_grid
      integer :: isom,nx,ny,nz,nvar1,nvar2,ierr,ix,iy,iz,ivar,current_index
      character(len=40) :: current_line,node_type
      real(kind=wp),allocatable :: Prototype_value(:,:)
   !
      isom=20;
      inquire(file=trim(som_fl),exist=testfl);
      if(.not. testfl) then
      write(*,*) 'ERROR: the som file does not exist'
      stop
      endif
   !
      write(*,*)
      write(*,*) 'SOM: Reading SOM Prototypes...'
      write(*,*)
      open(isom,file=trim(som_fl),status='unknown',action='read',access='sequential');
      read(isom,'(A)') current_line
      write(*,*) trim(current_line)
      read(isom,'(A17,1X,3I6)') current_line,nx,ny,nz
      write(*,*) current_line,nx,ny,nz
      read(isom,'(A21,1X,2I6)') current_line,nvar1,nvar2
      write(*,*) current_line,nvar1,nvar2
      kohonen_map%number_variables1=nvar1;
      kohonen_map%number_variables2=nvar2;
      kohonen_map%number_variables=nvar1*nvar2;
      read(isom,'(A25,1X,A11,1X,L4)') current_line,node_type,toroidal_grid
      write(*,*) current_line,node_type,toroidal_grid
      allocate(Prototype_value(nvar1*nvar2,1),stat=ierr);
      !
      select case(trim(layer_type))
      case('grid')
         if(allocated(kohonen_map%grid)) then
            do iz=1,size(kohonen_map%grid,3)
               do iy=1,size(kohonen_map%grid,2)
                  do ix=1,size(kohonen_map%grid,1)
                     call kohonen_map%grid(ix,iy,iz)%destroy();             
                  enddo
               enddo
            enddo
            deallocate(kohonen_map%grid);
         endif
         allocate(kohonen_map%grid(nx,ny,nz),stat=ierr);     
         if(allocated(kohonen_map%coordinates)) then
            deallocate(kohonen_map%coordinates);
         endif
         allocate(kohonen_map%coordinates(nx*ny*nz,3),stat=ierr);
         if(allocated(kohonen_map%cells_distances)) then
            deallocate(kohonen_map%cells_distances);
         endif
         allocate(kohonen_map%cells_distances(nx*ny*nz,nx*ny*nz),stat=ierr);
         if(allocated(kohonen_map%cluster_cells_index)) then
            deallocate(kohonen_map%cluster_cells_index);
         endif
         allocate(kohonen_map%cluster_cells_index(nx*ny*nz,4),stat=ierr);
         if(allocated(kohonen_map%grid_cluster)) then
            deallocate(kohonen_map%grid_cluster);
         endif
         allocate(kohonen_map%grid_cluster(nx,ny,nz),stat=ierr);
      case('cluster') 
         if(allocated(kohonen_map%cluster_layer) ) then
            do iz=1,size(kohonen_map%cluster_layer,1)
               call kohonen_map%cluster_layer(ix)%destroy();             
            enddo
            deallocate(kohonen_map%cluster_layer);
         endif
         allocate(kohonen_map%cluster_layer(nx),stat=ierr);
      end select
      do iz=1,nz
         read(isom,'(A)') current_line
         write(*,*) 'Reading ',trim(current_line)
         do iy=1,ny
            do ix=1,nx
               read(isom,'(A)') current_line
   !            write(*,*) current_line
               read(isom,'(A)') current_line
   !            write(*,*) current_line
               read(isom,*) (Prototype_value(ivar,1),ivar=1,nvar1*nvar2)
               !write(*,*) ix,iy,(Prototype_value(ivar,1),ivar=1,nvar1*nvar2)
               select case(trim(layer_type))
               case('grid')
                  call kohonen_map%grid(ix,iy,iz)%create(Prototype_value);
                  !call kohonen_map%grid(ix,iy,iz)%print()
                  current_index=position2index(ix,iy,iz,nx,ny);
                  call calculate_coordinates(current_index,ix,iy,iz,nx,ny,nz,kohonen_map%coordinates,node_type);
               case('cluster')
                  call kohonen_map%cluster_layer(ix)%set_prototype(Prototype_value);
               end select
            enddo
         enddo
      enddo
      close(isom)
   !
      if(trim(layer_type) .eq. 'grid') then
      call calculate_distance_matrix(kohonen_map%coordinates,kohonen_map%cells_distances,node_type,toroidal_grid);    
      endif   
   !
      write(*,*)
      write(*,*) 'SOM: Reading SOM Prototypes...finished'
      write(*,*)
   !
   end subroutine read_som_layer