code: def process_raw_litho1_data(model_path):
      
  
      npoints = 40962
      nlayers = len(l1_layer_decode)
      nentries = 9
  
      litho_data = np.ones((nlayers, nentries, npoints)) * -99999.0
  
      for model in range(0,npoints):
          model_name = "node"+str(model+1)+".model_tr"
          file = path.join(model_path, model_name)
  
          thisnodedata  = np.loadtxt(file, usecols=np.arange(0,9), comments="nlayer")
          thisnodeident = np.loadtxt(file, usecols=(9,), dtype=str, comments="nlayer")
  
          if (model % 1000 == 0):
              print ("Reading node ", model)
  
          for i, ident in enumerate(thisnodeident):
              litho_data[l1_layer_decode[ident], :,  model] = thisnodedata[i,:]
  
          ## Post process - not all layers in the model are populated (gives -99999 which is always illegal)
          ## For depth, it makes more sense to have the layer simply have no thickness but appear in the default order
  
          for layer in range(9,nlayers):
              missing_entries = np.where(litho_data[layer, l1_data_decode["DEPTH"]] == -99999)
              litho_data[layer, l1_data_decode["DEPTH"], missing_entries] = litho_data[layer-1, l1_data_decode["DEPTH"], missing_entries]
  
      grid_points_location = path.join(model_path,"Icosahedron_Level7_LatLon_mod.txt")
      litho_points = np.loadtxt( grid_points_location )
  
      return litho_data, litho_points
docs: Don't forget to strip the model data first
      truncate_raw_litho1_data(model path, truncated_model_path)