###Using the built-in smooth tool on two axes
^{Applicable to few cases}

1. Temporarily rotate your model so that one axis is close to parallel to the edge loop.
2. Temporarily remove edges adjacent to the edge loops, so the edgeloop which needs smoothing are selected and not connected to anything else.
3. Execute the smooth operation W < Smooth and set the Smoothing parameter to 0.5. Increase the iterations as desired.
4. Reconnected the edges with the bridge edgeloop tool.
5. Realign your model.

To script the functionality, start with a bmesh and the get the edge loop. Then implement your custom smoothing function.

I have hacked a short script together. This also works on curved surfaces.

1. Select the starting edge of an edge loop.
2. Select the ending edge of the edge loop. (They don't have to cover all of the edge loop.)
3. Edit the smoothing and iterations variables and execute the script.

This also works with multiple edge loop simultaneously.

_{Another application example.}

import bpy
import bmesh


def bm_from_me(me):
    if me.is_editmode:
        bm = bmesh.from_edit_mesh(me)
    else:
        bm = bmesh.new()
        bm.from_mesh(me)
    return bm
def bm_to_me(bm, me):
    if bm.is_wrapped:
        bmesh.update_edit_mesh(me, False, False)
    else:
        bm.to_mesh(me)
        me.update()

def get_loop_untill_selection(loop):
    verts = []
    while len(loop.vert.link_edges) == 4:
        loop = loop.link_loop_prev.link_loop_radial_prev.link_loop_prev
        verts.append(loop.vert)
        if loop.edge.select == True:
            loop.edge.select = False
            print(loop.edge.select)
            return verts, loop.edge
    return [], None

def smooth_coordinates(coords, smoothing, iterations):
    for i in range(1, len(coords) - 1):
        med = (coords[i - 1] + coords[i + 1]) / 2
        
        coords[i] *= 1 - smoothing
        coords[i] += med * smoothing
    
    if iterations > 1:
        smooth_coordinates(coords, smoothing, iterations - 1)

def smooth_loop(e, smoothing = 0.5, iterations = 1):
    edges = []
    for edge in bm.edges:
        if edge.select:
            e = edge
    
            for start_loop in e.link_loops:
                verts, end_edge = get_loop_untill_selection(start_loop)
                if end_edge is not None:
                    verts.insert(0, start_loop.vert)
                    if (start_loop.vert == e.verts[0]):
                        verts.insert(0, e.verts[1])
                    else:
                        verts.insert(0, e.verts[0])
                    for medge in [edge, end_edge]:
                        edges.append(medge)
                        medge.select = False
                    
                    if len(verts) > 2:
                        co = [v.co for v in verts]
                        smooth_coordinates(co, smoothing, iterations)
    for e in edges:
        e.select = True
        

me = bpy.context.object.data
bm = bm_from_me(me)
smoothing = 0.5
iterations = 2
smooth_loop(bm, smoothing, iterations)
bm_to_me(bm, me)