Problem with excel macro

Please help me

[VBA]
Private Sub CommandButton1_Click()
' Calcolo delle condizioni all'uscita del canale in rame
Dim M1, M2, M3, M4, M5, M6 As Double
Dim M3new, M5new As Double
Dim L1, L2, L3, L4 As Double
Dim L23, L56 As Double
Dim RE As Double
Dim eps As Double

' Dati geometrici
' Aree
A1 = 0.0257 'Ingresso canale
A2 = 0.00623 'Parte a sezione rettangolare 361
A3 = 0.00951 'Area immediatamente a monte delle pale, calcolata su un piano parallelo alle stesse
A4 = 0.00586 'Area di gola delle pale
A5 = 0.00622 'Area immediatamente a valle delle pale, calcolata su un piano parallelo alle stesse
A6 = 0.00361 'Area minima nel canale (piu' o meno all'altezza del gomito) Dai miei calcoli e' 0.00427!
A7 = A1 'Area di uscita del canale
' Altezza canale (sezione di prova)
H = 0.045 '45 mm
' Lunghezza canale (lungo l'ascissa curvilinea presa sull'asse) per calcolo perdite distribuite
L23 = 0.172 * 1# 'Lunghezza tratto isotermo a monte delle pale
L56 = 0.29 * 1# 'Lunghezza tratto isotermo a valle delle pale
' Diametri idraulici
Per1 = 4# * A1 / 3.14159
D1 = 4# * A1 / Per1
Per23 = 2# * (H + A2 / H)
D23 = 4# * A2 / Per23
Per56 = 2# * (H + A6 / H)
D56 = 4# * A6 / Per56
D7 = D1
' Rugosità assoluta [m]
eps = 0.00002 'Rugosità tipica di un tubo in rame
'''' Fine dati geometrici ''''''
'========================================================================== ==============
' Proprietà del fluido
gam = 1.33
r = 288#
nmax = Range("E2")
'print*,"Numero dati da elaborare?"
'read*,nmax
For niter = 1 To nmax
' Dati all'ingresso del canale
'print*,"G, p1,T1"
'read*, G,p1,T1
G = Cells(1 + niter, 1)
p1 = Cells(1 + niter, 2)
T1 = Cells(1 + niter, 3)

' Calcolo tratto isentropico 1-2
ro1 = p1 / (r * T1)
u1 = G / (ro1 * A1)
M1 = u1 / (gam * r * T1) ^ 0.5
'print*,"M1",M1,ro1,p1,u1
' Rapporto Area/Area critica nella sez. 1
A1sAc = area(M1, gam)
' Calcolo area critica tratto isentropico 1-2
Ac12 = A1 / A1sAc
If (A2 < Ac12) Then
' Se A2<Ac il flusso e' strozzato nella sez. 2
MsgBox "Condizioni soniche raggiunte nella sez. 2: A2,Ac12,A2,Ac12"
End If
' Grandezze totali in ingresso
p01 = p1 * (1# + 0.5 * (gam - 1#) * M1 ^ 2) ^ (gam / (gam - 1#))
T01 = T1 * (1# + 0.5 * (gam - 1#) * M1 ^ 2)
' Condizione di isentropia
p02 = p01
T02 = T01
' Area/area critica sez. 2
A2sAc = A2 / Ac12
' Calcolo delle condizioni nella sezione 2
StartMach = 0.1
' Calcolo del Mach 2 dalla relazione Area/Mach
M2 = calcM(A2sAc, gam, StartMach)
' Calcolo delle condizioni 2
p2 = p02 * (1# + 0.5 * (gam - 1#) * M2 ^ 2) ^ (-gam / (gam - 1#))
T2 = T02 * (1# + 0.5 * (gam - 1#) * M2 ^ 2) ^ (-1#)
u2 = M2 * (r * gam * T2) ^ 0.5
'print*,"M2,p2,T2,u2"
'print*,M2,p2,T2,u2
' Verifica conservazione della portata
'print*,"TEST G2=",p2*u2*A2/(r*T2)
'============ Fine tratto 1-2 ============
'======= Perdita di carico per restringimento di sezione =======
ro2 = p2 / (r * T2)
coef = 0.1 'E' un po' piu' grande per tenere conto delle perdite distribuite
' densità e velocita' costanti, solo perdita di pressione
p2cor = p2 - 0.5 * coef * ro2 * u2 ^ 2
' u2cor = Sqr(2# * (1 - coef)) * u2 'Perdita di energia cinetica
u2cor = u2
T2cor = T2 * p2cor * u2cor / p2 / u2
' T2cor = T2
' u2cor = u2 * T2cor / T2 * p2 / p2cor
corMach = u2cor / (gam * r * T2cor) ^ 0.5
'print*,"corMach,p2cor,T2cor,u2"
'print*,corMach,p2cor,T2cor,u2
p2 = p2cor
T2 = T2cor
u2 = u2cor
M2 = corMach
' Verifica conservazione della portata
' print*,"TEST G2cor=",p2*u2*A2/(r*T2)
'========= Perdite distribuite tratto 2-3 : ipotesi di flusso isotermo con attrito ==========
' Calcolo del numero di Reynolds
vis2 = visc(T2)
RE23 = (ro2 * u2 * D23) / vis2
RE = RE23
' print*,"Reynolds tratto 2-3",Re23
' Calcolo del coefficiente d'attrito

f = amoody(RE23, eps / D23)

' Calcolo iterativo del rapporto di pressione

'print*,M2,gam,f,L23,D23
beta = calcbeta(M2, gam, f, L23, D23)
'print*,beta
p3 = beta * p2
u3 = u2 / beta
T3 = T2
M3 = u3 / (r * gam * T3) ^ 0.5
'print*,"M3,p3,T3,u3"
'print*,M3,p3,T3,u3
'print*,"TEST G3=",p3*u3*A2/(r*T3)
'======== Fine perdite distribuite =========
' Calcolo tratto isentropico 3-5 a cavallo delle pale. L'area di gola e' A4
ro3 = p3 / (r * T3)
u3new = G / ro3 / A3 'A3 è l'area su un piano parallelo alle pale, quindi e' maggiore di A2
M3new = u3new / (gam * r * T3) ^ 0.5
' Area/Area critica
A3sAc = area(M3new, gam)
' Area critica tratto isentropico 3-5
Ac35 = A3 / A3sAc
If (Ac35 > A4) Then
' print*,"Area critica",Ac35,
'& "maggiore area di gola A4",A4
' print*,"Flusso supersonico nel divergente!!"
StarMach = 3#
Else
' print*,"Flusso subsonico nel divergente"
StartMach = 0.001
End If

' Grandezze totali sez. 3
p03 = p3 * (1# + 0.5 * (gam - 1#) * M3new ^ 2) ^ (gam / (gam - 1#))
T03 = T3 * (1# + 0.5 * (gam - 1#) * M3new ^ 2)
' Condizione di isentropia
p05 = p03
T05 = T03
' Area/Area critica sezione 5
A5sAc = A5 / Ac35
' Calcolo delle condizioni nella sezione 5 subito a valle delle pale
' Calcolo Mach dalla relazione Area/MAch
M5 = calcM(A5sAc, gam, StartMach)
' Condizioni sez. 5
p5 = p05 * (1# + 0.5 * (gam - 1#) * M5 ^ 2) ^ (-gam / (gam - 1#))
T5 = T05 * (1# + 0.5 * (gam - 1#) * M5 ^ 2) ^ (-1#)
u5 = M5 * (r * gam * T5) ^ 0.5
'print*,"M5,p5,T5,u5"
'print*,M5,p5,T5,u5
' Verifica conservazione della portata
'print*,"TEST G5=",p5*u5*A5/(r*T5)
'======== Fine tratto 3-5 =========
'======== Perdita di carico per gomito all'altezza delle palette =========
ro5 = p5 / (r * T5)
coef = 0.7 'tiene conto anche delle perdite distribuite del tratto 5-6
' densità e velocita' costanti, solo perdita di pressione
p5cor = p5 - 0.5 * coef * ro5 * u5 ^ 2 'perdita di pressione
' u5cor = Sqr(2# * (1# - coef)) * u5 'perdita di energia cinetica
u5cor = u5
T5cor = T5 * p5cor * u5cor / p5 / u5
' T5cor = T5
' u5cor = u5 * T5cor / T5 * p5 / p5cor
corMach = u5cor / (gam * r * T5cor) ^ 0.5
'print*,"corMach,p5cor,T5cor,u5"
'print*,corMach,p5cor,T5cor,u5
p5 = p5cor
T5 = T5cor
u5 = u5cor
M5 = corMach

' Verifica conservazione della portata
' print*,"TEST G5cor=",p5*u5*A5/(r*T5)
'======== Perdita di carico distribuita tratto 5-6 a valle delle palette : ipotesi di flusso isotermo =========
' Calcolo del numero di Reynolds
ro5 = p5 / (r * T5)
u5new = G / (ro5 * A6) 'prendo l'area trasversale, minore dell'area // alle pale
M5new = u5new / (gam * r * T5) ^ 0.5
vis5 = visc(T5)
Re56 = (ro5 * u5 * D6) / vis5
RE = Re56
' print*,"Reynolds tratto 5-6",Re56
' Calcolo del coefficiente d'attrito

f = amoody(Re56, eps / D56)
' Calcolo iterativo del rapporto di pressione
beta = calcbeta(M5new, gam, f, L56, D56)
p6 = beta * p5
u6 = u5new / beta
T6 = T5
M6 = u6 / (r * gam * T6) ^ 0.5
'print*,"M6,p6,T6,u6"
'print*,M6,p6,T6,u6
' Verifica conservazione della portata ERRORE!!!!
'print*,"TEST G6=",p6*u6*A6/(r*T6)
'======== Fine perdite distribuite =========

'========= Perdita di carico dovuta al distacco di vena a valle del primo gomito ========
ro6 = p6 / (r * T6)
coef = 1#
p6cor = p6 - 0.5 * coef * ro6 * u6 ^ 2 'perdita di pressione
' u6cor = Sqr(2# * (1# - coef)) * u6 'perdita di energia cinetica
u6cor = u6
T6cor = T6 * p6cor * u6cor / p6 / u6
' T6cor = T6
' u6cor = u6 * T6cor / T6 * p6 / p6cor
corMach = u6cor / (gam * r * T6cor) ^ 0.5
'print*,"corMach,p6cor,T6cor,u6"
'print*,corMach,p6cor,T6cor,u6
p6 = p6cor
T6 = T6cor
u6 = u6cor
M6 = corMach
'========================================================================== ======
'== ==
'== Attenzione: a causa del distacco di vena a valle del gomito, ==
'== il successivo allargamento di sezione non ha praticamente alcun effetto ==
'== sulla pressione, che resta praticamente costante (zona di ricircolo) ==
'== per questa ragione, il calcolo del divergente a partire dall'ipotesi di ==
'== isentropia del flusso e' totalmente inappropriata. E' piu' accurato ==
'== ipotizzare che le condizioni del flusso restino invariate fino all'uscita ==
'== ==
'========================================================================== ======
' Verifica conservazione della portata
'print*,"TEST G6cor=",p6*u6*A6/(r*T6)
Worksheets("Risultati").Range("A25") = G
Worksheets("Risultati").Cells(2 + niter, 1) = G
Worksheets("Risultati").Range("B25") = p1
Worksheets("Risultati").Cells(2 + niter, 2) = p1
Worksheets("Risultati").Range("C25") = T1
Worksheets("Risultati").Cells(2 + niter, 3) = T1
Worksheets("Risultati").Range("D25") = p6
Worksheets("Risultati").Cells(2 + niter, 4) = p6
Worksheets("Risultati").Range("E25") = M1
Worksheets("Risultati").Cells(2 + niter, 5) = M1
Worksheets("Risultati").Range("F25") = M2
Worksheets("Risultati").Cells(2 + niter, 6) = M2
Worksheets("Risultati").Range("G25") = M3
Worksheets("Risultati").Cells(2 + niter, 7) = M3
Worksheets("Risultati").Range("H25") = M5
Worksheets("Risultati").Cells(2 + niter, 8) = M5
Worksheets("Risultati").Range("I25") = M6
Worksheets("Risultati").Cells(2 + niter, 9) = M6
portrid = G * Sqr(T01) / p01 * 1000
Worksheets("Risultati").Range("L25") = portrid
Worksheets("Risultati").Cells(2 + niter, 10) = portrid
rapp = p01 / p6
Worksheets("Risultati").Range("M25") = rapp
Worksheets("Risultati").Cells(2 + niter, 11) = rapp

'write(10,*) G,p1,T1,p6,T6,M1,M2,M3,M5,M6
'write(20,*) G*sqrt(T01)/p01*1000,p01/p6
Next niter
End Sub
'======================
'===== FUNZIONI =======
'======================
Function area(M, gam)
' Relazione Area/Mach
gamp1 = gam + 1#
gamm1 = gam - 1#
area = (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M
End Function

'===============================
Function calcM(AsAc, gam, StartMach)
' Calcolo inverso del Mach dalla legge delle aree

gamp1 = gam + 1#
gamm1 = gam - 1#

' print*,"A/A*",AsAc
M = StartMach
dM = 10000000000#

' Calcolo iterativo con Newton-Raphson
Do

f = (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M - AsAc
fp = -(2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1)) / M ^ 2 + (2# / gamp1 * (1# + gamm1 / 2# * M ^ 2)) ^ (gamp1 / (2# * gamm1) - 1)
dM = -f / (fp + 0.00000001)
M = M + dM
' print*,M,dM
Loop Until (Abs(dM) > 0.00001)
calcM = M
End Function
'====================================================================

Function visc(T)
' Calcolo della viscosità dalla formula di Sutherland

T0 = 288.15
S = 110#
visc = 0.0000178 * (T / T0) ^ 1.5 * (T0 + S) / (T + S)
End Function
'============================
Function amoody(RE, eos)
' Calcolo del coefficiente d'attrito dalla formula di Colebrook
' Valore stimato dalla formula di Swamee-Jain

f = 0.25 / ((Log(eos / 3.72 + 5.74 / RE ^ 0.9)) / 2.302585) ^ 2

' Calcolo iterativo
usrf = 1# / Sqr(f)
Do
usrf1 = usrf
usrf = -2# * (Log(2.51 * usfr1 / RE + eps / 3.72)) / 2.302585
' print*,usrf,usrf1
Loop Until Abs(usrf1 - usrf) > 0.000001

amoody = usrf ^ (-2#)

End Function
'==========================
Function calcbeta(M, gam, f, L, D)
' Calcolo del rapporto di pressione per flusso isotermo

'10 continue
beta = 0#
beta1 = 1#
nit = 0
Do
nit = nit + 1
beta = beta1
B = Log(beta ^ (-1#))
A = 1# - 2# * gam * (M ^ 2) * (0.5 * B + 2# * f * L / D)
If A > 0.8 Then beta1 = A
If A < 0.8 Then beta1 = 0.8
' Limitatore inserito per evitare problemi in caso di divergenza

Loop Until (Abs(beta - beta1) > 0.0001)
'print*,sqrt(beta)
calcbeta = Sqr(beta)

End Function
'========================================================================== ===================
[/VBA]