Thank you for the reply. Following your suggestion I set the Navier-Stokes solver like following in order to avoid newton method

Nonlinear System Max Iterations = 50

Nonlinear System Convergence Tolerance = 1.0e-5

Nonlinear System Newton After Iterations = 50

Nonlinear System Newton After Tolerance = 1.0e-02

Nonlinear System Relaxation Factor = 1.00

Nonlinear System Reset Newton = Logical True

But unfortunately this didn't solve the problem. Looking at the terminal output I can see that the calculation diverges AFTER the first Navier-Stokes Solver. I post the terminal output below. It seems that it crashes when using my external function "getBalance". BUT that function works perfectly if I restore the boundary condition on the bedrock to "0-velocity". I attach also the code of my external functions named "utils.f90". I thank you and anybody else that can spend even a couple of minutes with my problem.

Here is the function code

And here the output:

ELMER SOLVER (v 7.0) STARTED AT: 2015/08/22 11:53:45

MAIN:

MAIN: =============================================================

MAIN: ElmerSolver finite element software, Welcome!

MAIN: This program is free software licensed under (L)GPL

MAIN: Copyright 1st April 1995 - , CSC - IT Center for Science Ltd.

MAIN: Webpage

http://www.csc.fi/elmer, Email

elmeradm@csc.fiMAIN: Library version: 7.0 (Rev: 1b38bc5)

MAIN: HYPRE library linked in.

MAIN: =============================================================

MAIN:

MAIN:

MAIN: -------------------------------------

MAIN: Reading Model: 1997-2007_test_Slip_22-08-2015.sif

Model Input: Unlisted keyword: [water density] in section: [constants]

Model Input: Unlisted keyword: [surfacedem] in section: [constants]

Model Input: Unlisted keyword: [bedrockdem] in section: [constants]

Model Input: Unlisted keyword: [switchperiod] in section: [constants]

Model Input: Unlisted keyword: [winterarrheniusfactor] in section: [constants]

Model Input: Unlisted keyword: [summerarrheniusfactor] in section: [constants]

Model Input: Unlisted keyword: [basalwintervelocity] in section: [constants]

Model Input: Unlisted keyword: [basalsummervelocity] in section: [constants]

Model Input: Unlisted keyword: [winteraccumulationfile] in section: [constants]

Model Input: Unlisted keyword: [summeraccumulationfile] in section: [constants]

Model Input: Unlisted keyword: [minicethickness] in section: [constants]

Model Input: Unlisted keyword: [zs accumulation] in section: [body force 2]

Loading user function library: [./utils_eros_fede]...[getBalance]

Model Input: Unlisted keyword: [zs] in section: [initial condition 2]

Loading user function library: [./utils_eros_fede]...[getSurfaceElevation]

Loading user function library: [./utils_eros_fede]...[getBedrockElevation]

Model Input: Unlisted keyword: [set arrhenius factor] in section: [material 1]

Model Input: Unlisted keyword: [arrhenius factor] in section: [material 1]

Loading user function library: [./utils_eros_fede]...[getArrheniusFactor]

Model Input: Unlisted keyword: [min zs] in section: [material 1]

Loading user function library: [./utils_eros_fede]...[getMinSurfaceElevation]

Model Input: Unlisted keyword: [max zs] in section: [material 1]

Model Input: Unlisted keyword: [mesh velocity variable] in section: [solver 1]

Model Input: Unlisted keyword: [mesh update variable] in section: [solver 1]

Model Input: Unlisted keyword: [mesh velocity first zero] in section: [solver 1]

Model Input: Unlisted keyword: [top surface variable name] in section: [solver 1]

Model Input: Unlisted keyword: [displacement mode] in section: [solver 1]

Model Input: Unlisted keyword: [correct surface] in section: [solver 1]

Model Input: Unlisted keyword: [minimum height] in section: [solver 1]

Model Input: Unlisted keyword: [stress variable name] in section: [solver 4]

Model Input: Unlisted keyword: [flow solver name] in section: [solver 4]

Model Input: Unlisted keyword: [relaxation factor] in section: [solver 6]

MAIN:

MAIN: -------------------------------------

MAIN: Time: 1/24 2592000.0000000000

MAIN: -------------------------------------

MAIN:

SolveEquations: -------------------------------------

SolveEquations: Coupled system iteration: 1

SolveEquations: -------------------------------------

StructuredMeshMapper: >Correct Surface< in case of intersecting upper and lower surface

StructuredMeshMapper: Adjusting upper surface to maintain minimum height to 0.5000E+00

getArrheniusFactor: Set value of Arrhenius Factor to 2.4000000E-24

ComputeChange: NS (ITER=1) (NRM,RELC): ( 0.66460863E-04 2.0000000 ) :: navier-stokes

ComputeChange: NS (ITER=2) (NRM,RELC): ( 0.10736614E-03 0.47064352 ) :: navier-stokes

ComputeChange: NS (ITER=3) (NRM,RELC): ( 0.14503338E-03 0.29847314 ) :: navier-stokes

ComputeChange: NS (ITER=4) (NRM,RELC): ( 0.17635048E-03 0.19488907 ) :: navier-stokes

ComputeChange: NS (ITER=5) (NRM,RELC): ( 0.20055361E-03 0.12843124 ) :: navier-stokes

ComputeChange: NS (ITER=6) (NRM,RELC): ( 0.21831422E-03 0.84802943E-01 ) :: navier-stokes

ComputeChange: NS (ITER=7) (NRM,RELC): ( 0.23089744E-03 0.56023585E-01 ) :: navier-stokes

ComputeChange: NS (ITER=8) (NRM,RELC): ( 0.23960695E-03 0.37021991E-01 ) :: navier-stokes

ComputeChange: NS (ITER=9) (NRM,RELC): ( 0.24554298E-03 0.24470931E-01 ) :: navier-stokes

ComputeChange: NS (ITER=10) (NRM,RELC): ( 0.24954789E-03 0.16178463E-01 ) :: navier-stokes

ComputeChange: NS (ITER=11) (NRM,RELC): ( 0.25223184E-03 0.10697718E-01 ) :: navier-stokes

ComputeChange: NS (ITER=12) (NRM,RELC): ( 0.25402263E-03 0.70747007E-02 ) :: navier-stokes

ComputeChange: NS (ITER=13) (NRM,RELC): ( 0.25726561E-03 0.12685526E-01 ) :: navier-stokes

ComputeChange: NS (ITER=14) (NRM,RELC): ( 0.25753787E-03 0.10577202E-02 ) :: navier-stokes

ComputeChange: NS (ITER=15) (NRM,RELC): ( 0.25755675E-03 0.73293926E-04 ) :: navier-stokes

ComputeChange: NS (ITER=16) (NRM,RELC): ( 0.25755835E-03 0.62131185E-05 ) :: navier-stokes

ComputeChange: SS (ITER=1) (NRM,RELC): ( 0.25755835E-03 2.0000000 ) :: navier-stokes

ComputeChange: NS (ITER=1) (NRM,RELC): ( 394880.79 2.0000000 ) :: sij

ComputeChange: NS (ITER=2) (NRM,RELC): ( 421559.88 0.65354617E-01 ) :: sij

ComputeChange: NS (ITER=3) (NRM,RELC): ( 431568.68 0.23463770E-01 ) :: sij

ComputeChange: NS (ITER=4) (NRM,RELC): ( 309568.48 0.32922437 ) :: sij

ComputeChange: NS (ITER=5) (NRM,RELC): ( 235878.41 0.27020070 ) :: sij

ComputeChange: NS (ITER=6) (NRM,RELC): ( 239879.73 0.16820801E-01 ) :: sij

ComputeChange: SS (ITER=1) (NRM,RELC): ( 239879.73 2.0000000 ) :: sij

ComputeChange: SS (ITER=1) (NRM,RELC): ( 0.0000000 0.0000000 ) :: eigenstresses

getBalance: Set accumulation file to balance_winter.dat

3 0.4927E-01

6 0.4202E+01

9 0.3234E+01

12 0.1999E+01

15 0.2004E+01

18 0.2006E+01

21 0.2006E+01

...

Cutted by me

...

1026 0.1892+139

1029 0.6784+140

1032 0.1688+142

1032 0.1688+142

ERROR:: IterSolve: Failed convergence tolerances.