function run_crack_propa_job(param0,model0)
nmod=0;
xos=[];
yos=[];
conns=[];
Us=[];
Kss=[];
das=[];
roi=param0.roi;
param=param0;
model=model0;
KIc=model0.material_parameters.KIc;
filreso=param.result_file;
filres='upropa.res';
param.result_file=filres;
model.mesh_file='upropa.vtk';
copyfile(model0.mesh_file,model.mesh_file);
pix2m=param.pixel_size;
xfem=model.phantom_nodes;
dao=param.da;
modes=[1,2];
if xfem
    ind=1:4;
else
    ind=0:7;
end
mu=0.5*model.material_parameters.young/(1+model.material_parameters.nu);
Es=model.material_parameters.young/(1-(model.material_parameters.nu)^2);
scal=2*mu*sqrt(2*pi);
scalamp=scal*(pix2m.^(1-(ind)*.5));
found=find(ind==1);
tipso=zeros(size(model.zone,2),2);
nw=length(modes)*length(ind);
ntip=0;
zones=model.zone;
for iz=1:size(zones,2)
    zone=zones(:,iz);
    switch zone{4}
        case 5
            if (zone{8}>0)
                tipso(iz,1)=ntip+1;
                ntip=ntip+1;
                if  (zone{9}>0)
                    tipso(iz,2)=ntip+1;
                    ntip=ntip+1;
                end
            end
        case 6
            model0.zone{7,iz}=[];
    end
end
cracks=find(cell2mat(zones(4,:))==5);
tipso(~(cell2mat(zones(4,:))==5),:)=[];
step=1;
filreso=strrep(filreso,'.res','');
filres=strrep(filres,'.res','');
delete([fullfile('VTK',filreso),'-0*.vtk']);
copyfile([filreso,'.dat'],[filres,'.dat']);
remesh=0*tipso;
dacum=0*tipso;
while 1
    if xfem
        run_fdfea_job(param,model)
    else
        run_fea_job(param,model)
    end
    modeli=model;
    load([filres,'.res'],'-mat','model','U','Ks','tips','xo','yo','conn')
    if step==1
        movefile(fullfile('VTK',[filres , sprintf('-%06d.vtk',0)]),fullfile('VTK',[filreso , sprintf('-%06d.vtk',0)]));
    end
    movefile(fullfile('VTK',[filres , sprintf('-%06d.vtk',1)]),fullfile('VTK',[filreso , sprintf('-%06d.vtk',step)]));
    if any(tips(:)>0)
        xos{step}=xo;
        yos{step}=yo;
        conns{step}=conn;
        lf=0*tips+Inf;
        dasi=0*tips;
        Ksi=zeros(sum(tipso(:)>0)*nw,1);
        for ic=1:size(tipso,1)
            for it=1:2
                if (tipso(ic,it)>=1)&&(tips(ic,it)>=1)
                    if xfem
                        K1=sqrt(pix2m)*Ks(1+2*(cracks(ic)-1)+(it-1),1);
                        K2=sqrt(pix2m)*Ks(1+2*(cracks(ic)-1)+(it-1),2);
                    else
                        K1=scalamp(found)*Ks(found+(tips(ic,it)-1)*nw);
                        K2=scalamp(found)*Ks(found+length(ind)+(tips(ic,it)-1)*nw);
                    end
                    [K1eq,Oc]=MaxHoopStressCriterion(K1,K2);
                    lf(ic,it)=KIc/K1eq;
                end
            end
        end
        Us{step}=min(lf)*U;
        Ks=min(lf)*Ks;
        for iz=1:size(zones,2)
            zone=model0.zone(:,iz);
            switch zone{4}
                case 6
                    fu=model0.zone{7,iz};
                    for izz=1:size(model.zone,2)
                        if strcmp(get(model.zone{3,izz},'Tag'),get(model0.zone{3,iz},'Tag'))
                            
                            model0.zone{7,iz}=[fu,min(lf)*model.zone{7,izz}];
                            break
                        end
                    end
            end
        end
        model=modeli;
        for ic=1:size(tipso,1)
            for it=1:2
                if (tipso(ic,it)>=1)&&(tips(ic,it)>=1)
                    if xfem
                        Ksi((1:nw)+(tipso(ic,it)-1)*nw)=Ks(1+2*(cracks(ic)-1)+(it-1),1:nw)';
                        K1=sqrt(pix2m)*Ks(1+2*(cracks(ic)-1)+(it-1),1);
                        K2=sqrt(pix2m)*Ks(1+2*(cracks(ic)-1)+(it-1),2);
                    else
                        Ksi((1:nw)+(tipso(ic,it)-1)*nw)=Ks((1:nw)+(tips(ic,it)-1)*nw);
                        K1=scalamp(found)*Ks(found+(tips(ic,it)-1)*nw);
                        K2=scalamp(found)*Ks(found+length(ind)+(tips(ic,it)-1)*nw);
                    end
                    [K1eq,Oc]=MaxHoopStressCriterion(K1,K2);
%                    Oc=0
                    if abs(Oc)>pi/6
                        remesh(ic,it)=1;
                    end
                    if min(lf)/lf(ic,it)>0.1
                        dasi(ic,it)=min(lf)/lf(ic,it)*dao;
                        
                        path=zones{2,cracks(ic)}*[1;1i];
                        t=-gradient(path);
                        t=t./abs(t);
                        switch it
                            case 1
                                t=t(1);
                                n=t*exp(1i*pi/2);
                                tip=path(1)+dasi(ic,it)*(cos(Oc)*t+sin(Oc)*n);
                                if GetSignedDistanceToZone(model,roi,real(tip),imag(tip))<-(2-1.*xfem)*zones{8,cracks(ic)}
                                    path=[tip;path];
                                else
                                    dasi(ic,it)=0;
                                end
                                
                            case 2
                                t=t(end);
                                n=t*exp(1i*pi/2);
                                tip=path(1)+dasi(ic,it)*(cos(Oc)*t+sin(Oc)*n);
                                if GetSignedDistanceToZone(model,roi,real(tip),imag(tip))<-(2-1.*xfem)*zones{8,cracks(ic)}
                                    path=[path;tip];
                                else
                                    dasi(ic,it)=0;
                                end
                        end
                        dacum(ic,it)=dacum(ic,it)+dasi(ic,it);
                        if dacum(ic,it)>min(5*dao,min(2*zones{8,cracks(ic)},10*zones{7,cracks(ic)}))
                            remesh(ic,it)=1;
                        end
                        zones{2,cracks(ic)}=[real(path),imag(path)];
                    end
                end
            end
        end
        Kss=[Kss,Ksi];
        das=[das,dasi];
        if any(abs(dasi)>0)
            model.zone=zones;
            if ~xfem
                model=Remeshing(param,model,~any(remesh==1));
                if any(remesh==1)
                    remesh=0*remesh;
                    dacum=0*dacum;
                end
            end
        else
            break;
        end
    else
        model=modeli;
        break;
    end
    step=step+1;
end
da=das;
Ks=Kss;
tips=tipso;
xo=xos;yo=yos;
conn=conns;
U=Us;
param=param0;
for iz=1:size(model.zone,2)
    zone=model.zone(:,iz);
    switch zone{4}
        case 6
            for izz=1:size(model0.zone,2)
                if strcmp(get(model0.zone{3,izz},'Tag'),get(model.zone{3,iz},'Tag'))
                    model.zone{7,iz}=model0.zone{7,izz};
                    break
                end
            end
    end
end

model.zone(:,~((cell2mat(model.zone(4,:))==5)|(cell2mat(model.zone(4,:))==6)))=[];
cracks=find(cell2mat(model.zone(4,:))==5);
model0.zone=model.zone;
model=model0;

load(fullfile('TMP','0_mesh_0'),'ng','rflag','rint');
save(param.result_file,'U','Ks','tips','cracks','xo','yo','param','model','nmod','conn','rint','ng','rflag','-v7.3');
end