Organize Triangles of a Subdivided Quad?

How can I efficiently reorder the triangles of a plane, on different levels of subdivision, in a list and have triangle pairs that constitute quads? I know that there is some sort of mathematical pattern but I can't, for the life of me, figure it out. To visually explain my question further I have an image below showing my desired result, on the first three levels of subdivision. The numbers on the triangles are intended to represent that triangle's index in the array/list of triangles.

The code I found for subdivision is below, it's just simple subdivision without shared vertices (however, I do plan on writing a script that will share vertices).

 function subdivide() {
  
     verts = mesh.vertices;
     trigs = mesh.triangles;
     uvs = mesh.uv;
     norms = mesh.normals;
  
     Debug.Log("enter subdividing: "+verts.length);
  
     var nv = new Array(verts);
     var nt = new Array();
     var nu = new Array(uvs);
     var nn = new Array(norms);
  
     for(var i : int = 2;i<trigs.length;i+=3) {
  
         var p0trigwho : int = trigs[i-2];
         var p1trigwho : int = trigs[i-1];
         var p2trigwho : int = trigs[i];
 
         var p0trigwhere : int = i-2;
         var p1trigwhere : int = i-1;
         var p2trigwhere : int = i;
 
         var p0 : Vector3 = verts[p0trigwho];
         var p1 : Vector3 = verts[p1trigwho];
         var p2 : Vector3 = verts[p2trigwho];
 
         var pn0 : Vector3 = norms[p0trigwho];
         var pn1 : Vector3 = norms[p1trigwho];
         var pn2 : Vector3 = norms[p2trigwho];
 
         var pu0 : Vector2 = uvs[p0trigwho];
         var pu1 : Vector2 = uvs[p1trigwho];
         var pu2 : Vector2 = uvs[p2trigwho];
 
         var p0mod : Vector3 = (p0+p1)/2;    
         var p1mod : Vector3 = (p1+p2)/2;
         var p2mod : Vector3 = (p0+p2)/2;
 
         var pn0mod : Vector3 = ((pn0+pn1)/2).normalized;    
         var pn1mod : Vector3 = ((pn1+pn2)/2).normalized;
         var pn2mod : Vector3 = ((pn0+pn2)/2).normalized;
 
         var pu0mod : Vector2 = (pu0+pu1)/2;    
         var pu1mod : Vector2 = (pu1+pu2)/2;
         var pu2mod : Vector2 = (pu0+pu2)/2;
 
         var p0modtrigwho = nv.length;
         var p1modtrigwho = nv.length+1;
         var p2modtrigwho = nv.length+2;
 
         nv.push(p0mod);
         nv.push(p1mod);
         nv.push(p2mod);
 
         nn.push(pn0mod);
         nn.push(pn1mod);
         nn.push(pn2mod);
 
         nu.push(pu0mod);
         nu.push(pu1mod);
         nu.push(pu2mod);
 
         nt.Add(p0trigwho);
         nt.Add(p0modtrigwho);
         nt.Add(p2modtrigwho);
         
         nt.Add(p0modtrigwho);
         nt.Add(p1modtrigwho);
         nt.Add(p2modtrigwho);
         
         nt.Add(p2modtrigwho);
         nt.Add(p1modtrigwho);
         nt.Add(p2trigwho);
         
         nt.Add(p0modtrigwho);
         nt.Add(p1trigwho);
         nt.Add(p1modtrigwho);
         
     }    
  
     verts = nv.ToBuiltin(Vector3);
     norms = nn.ToBuiltin(Vector3);
     uvs = nu.ToBuiltin(Vector2);
     trigs = nt.ToBuiltin(int);
  
     //applyuvs();
     applymesh();
     //mesh.RecalculateNormals();
  
     Debug.Log("exit subdividing: "+trigs.length);
 }