Boundary subdomains for conjugate heat transfer problem

Good morning everyone,
I am currently trying to model the cooling time for a square (so 2D) of a material and I am not very clear on how to implement the different subdomains.
To be specific: the system at hand should have Neuman BC (outward heat flux) on the perimeter of the square and Robin BC (outward convection) on the surface, within the perimeter (the latter excluded). Hopefully the picture will clarify (fluxes represented normal to the side/surface the exit):

This is my current implementation of the Neumann BCs:

L = 140E-6    
#create the mesh, 2D
    domain = Rectangle(Point(0.0, 0.0), Point(L, L))
    mesh = generate_mesh(domain, 15)
    V = FunctionSpace(mesh, 'P', 1)

#subdomain on which the flux is 0
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, tol) or near(x[0], 1, tol)

#subdomain on which the flux is g
class Omega_1(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[1], 0, tol) or near(x[1], 1, tol) 

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
subdomain_0 = Omega_0()
subdomain_0.mark(boundary_markers, 0)
subdomain_1 = Omega_1()
subdomain_1.mark(boundary_markers, 1)

#redefine boundary integration measure
ds = Measure('ds', domain = mesh, subdomain_data = boundary_markers)

then to actually apply the BC and compute the integrals over the boundary:

g = Constant(Q)
boundary_conditions = {0 : {'Neumann' : 0}, 1 : {'Neumann' : g}}
bcs=[] #empty BC list, no dirichlet

integrals_N = []

for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(1/kappa*g*v*ds(i))

So, provided that what done until now is correct (on which I’d appreciate some feedback since I’m new with Fenics), how should I introduce a Robin BC that takes into account the convection from the square surface?
Thanks for your collaboration!

For robin conditions, consider this tutorial

thanks dokken, that is actually where I learned how to implement the above code.

What I am missing is how to define the Robin BC within the square, or if you will, how to define the area of the square as a subdomain on which then I can apply the Robin BC as in the tutorial.

Thanks!