How to subdivide a single face on a mesh?

I have a sphere being created with the code from the procedural primitives wiki page, i need to subdivide a single arbitrary face this is what i got so far but i'm out of ideas.

 void CreateSphere(Vector3 dir)
     {
         MeshFilter filter = gameObject.GetComponent<MeshFilter>();
         Mesh mesh = filter.mesh;
         mesh.Clear();
 
         int lt, ln;
         lt = 3;// (int)Vector3.Dot(new Vector3(dir.x, 0, dir.z), Vector3.forward); test values
         ln = 3;// (int)Vector3.Angle(new Vector3(0, dir.y, 0), Vector3.up); test values
 
         #region Vertices
         List<Vector3> vertices = new List<Vector3>();
         float _pi = Mathf.PI, _2pi = _pi * 2f;
 
         vertices.Add(Vector3.up * radius);
         for (int lat = 0; lat < nbLat; lat++)
         {
             float a1 = _pi * (float)(lat + 1) / (nbLat + 1);
             float sin1 = Mathf.Sin(a1);
             float cos1 = Mathf.Cos(a1);
             
             for (int lon = 0; lon <= nbLong; lon++)
             {
                 float a2 = _2pi * (float)(lon == nbLong ? 0 : lon) / nbLong;
                 float sin2 = Mathf.Sin(a2);
                 float cos2 = Mathf.Cos(a2);                
 
 #region subdivideFace
                 if(lat == lt && lon == ln)
                 {
                     float aD1, aD2, sinD1, cosD1, sinD2, cosD2;
                     for(int d = 0; d < detailRes; d++)
                     {
                         aD1 = a1 + (d / detailRes);
                         sinD1 = Mathf.Sin(aD1);
                         cosD1 = Mathf.Cos(aD1);
 
                         for (int d2 = 0; d2 < detailRes; d2++)
                         {
                             aD2 = a2 + (d2 / detailRes);
                             sinD2 = Mathf.Sin(aD2);
                             cosD2 = Mathf.Cos(aD2);
 
                             vertices.Add(new Vector3(sinD1 * cosD2, cosD1, sinD1 * sinD2) * radius);
                         }
                     }
                 }
 #endregion
                 else
                     vertices.Add(new Vector3(sin1 * cos2, cos1, sin1 * sin2) * radius);
             }
         }
         vertices.Add(Vector3.up * -radius);
         #endregion
 
         #region Normales
         Vector3[] normales = new Vector3[vertices.Count];
         for (int n = 0; n < vertices.Count; n++)
             normales[n] = vertices[n].normalized;
         #endregion
 
         #region UVs
         Vector2[] uvs = new Vector2[vertices.Count];
         uvs[0] = Vector2.up;
         uvs[uvs.Length - 1] = Vector2.zero;
         for (int lat = 0; lat < nbLat; lat++)
             for (int lon = 0; lon <= nbLong; lon++)
                 uvs[lon + lat * (nbLong + 1) + 1] = new Vector2((float)lon / nbLong, 1f - (float)(lat + 1) / (nbLat + 1));
         #endregion
 
         #region Triangles
         int nbFaces = vertices.Count;
         int nbTriangles = nbFaces * 2;
         int nbIndexes = nbTriangles * 3;
         int[] triangles = new int[nbIndexes];
 
         //Top Cap
         int i = 0;
         for (int lon = 0; lon < nbLong; lon++)
         {
             triangles[i++] = lon + 2;
             triangles[i++] = lon + 1;
             triangles[i++] = 0;
         }
 
         //Middle
         for (int lat = 0; lat < nbLat - 1; lat++)
         {
             for (int lon = 0; lon < nbLong; lon++)
             {
                 int current = lon + lat * (nbLong + 1) + 1;
                 int next = current + nbLong + 1;
 
                 triangles[i++] = current;
                 triangles[i++] = current + 1;
                 triangles[i++] = next + 1;
 
                 triangles[i++] = current;
                 triangles[i++] = next + 1;
                 triangles[i++] = next;
             }
         }
 
         //Bottom Cap
         for (int lon = 0; lon < nbLong; lon++)
         {
             triangles[i++] = vertices.Count - 1;
             triangles[i++] = vertices.Count - (lon + 2) - 1;
             triangles[i++] = vertices.Count - (lon + 1) - 1;
         }
         #endregion
 
         mesh.vertices = vertices.ToArray();
         mesh.normals = normales;
         mesh.uv = uvs;
         mesh.triangles = triangles;
 
         mesh.RecalculateBounds();
         mesh.Optimize();
     }
 }

thanks in advance