Nested for loop strange problem

[mcva]

Thank you for taking some time to read this problem...
I have this nested for loop and I´m note sure its doing what I need.
I need to run a lake temperature model for N = 365 and for each day (i) the first loop calculates the values off b0,c0,d0,an,bn and dn. The second loop has to go trough 40 layers (j)and calculate the coefficients ai, bi, ci, and di. Then the code line:

TempTDMA=TDMAsolver(a[i-1], b[i-1], c[i-1], d[i-1])

returns a list of values after applying the TDMAsolver function for each day(i) with the following input: b0,c0,d0, an,bn and dn, a1...aNlayers-1, b1...bNlayers-1, c1...cNlayers-1 and d1...dNlayers-1

the strange thing is that results seem to be ok, but I can´t remove the following line that´s no used by the model:

CANT_REMOVE_THIS_LINE=fix_profile(n,TempTDMA,Areat0)

I get an error in a function thats note linked with this code line, I can change the name of the variable associated with the fix_profile funtion ("CANT_REMOVE_THIS_LINE2) and the model runs!


this are the neested loops:

N=365
Nlayers=40
for i in xrange(1,N+1):
    I=np.linspace(1,N,N)
    Radiation[i-1] = EquilibriumTemp(air[i-1], RelativeHumidity[i-1], wind[i-1], ClearSkyRadiation[i-1], TempLake[i-1][0])
    Thermocline_Depth[i-1] = ThermoTemp(TempLake[i-1], z, Smin=0.1, seasonal=True, index=False, mixed_cutoff=1)
    K[i-1]=Keddy(TempLake[i-1],z,wind[i-1],Areat0)
    K[i]=Keddy(TempLake[i-1],z,wind[i-1],Areat0)
    
    #surface
    #b0, c0, d0
    b[i-1][0]=1+alfa[0]*s*(0.0+2.0*K[i][0]+K[i][1])
    c[i-1][0]=-alfa[0]*s*(0.0+2.0*K[i][0]+K[i][1])
    d[i-1][0]=((1.0-alfa[0]*s*(0.0+2.0*K[i-1][0]+K[i-1][1]))*TempLake[i-1][0])+(alfa[0]*s*(0.0+2.0*K[i-1][0]+K[i-1][1])*TempLake[i-1][1])+alfa[0]*s*(0.0+K[i-1][0])*((2*Radiation[i-1]*DZ)/(Cp*DH2O*(K[i-1][0]/Areat0[0])))+alfa[0]*s*(0.0+K[i][0])*((2.0*Radiation[i-1]*DZ)/(DH2O*Cp*(K[i][0]/Areat0[0])))
  
    #Bottom
    #an, bn, dn
    a[i-1][-1]=-alfa[-1]*s*(K[i][-2]+2.0*K[i][-1]+0.0)
    b[i-1][-1]=1+alfa[-1]*s*(K[i][-2]+2.0*K[i][-1]+0.0)
    d[i-1][-1]=alfa[-1]*s*TempLake[i-1][-2]*(K[i-1][-2]+2.0*K[i-1][-1]+0.0)+(1-alfa[-1]*s*(K[i-1][-2]+2*K[i-1][-1]+0.0))*TempLake[i-1][-1] # 0,1,2,3,4,5,6,7,8,9 
    
    
    #ai, bi, ci, di
    for j in xrange(1,Nlayers-1):# 40
        a[i-1][j-1]=-alfa[j]*s*(K[i][j-1]+K[i][j]) # 40
        b[i-1][j]=1.0+alfa[j]*s*(K[i][j-1]+2.0*K[i][j]+K[i][j+1])
        c[i-1][j]=-alfa[j]*s*(K[i][j]+K[i][j+1])
        d[i-1][j]=alfa[j]*s*(K[i-1][j]+K[i-1][j-1])*TempLake[i-1][j-1]+alfa[j]*s*(K[i-1][j+1]+K[i-1][j])*TempLake[i-1][j+1]+(1-alfa[j]*s*(K[i-1][j-1]+2.0*K[i-1][j]+K[i-1][j+1]))*TempLake[i-1][j]

        
        TempTDMA=TDMAsolver(a[i-1], b[i-1], c[i-1], d[i-1])
        n=len(TempTDMA)
        CANT_REMOVE_THIS_LINE=fix_profile(n,TempTDMA,Areat0)
    TempLake[i]=TempTDMA

the error ocurs inside the function ThermoTemp:

Error:
Traceback (most recent call last):
  File "C:\Users\manue\Google Drive\Lake model\PythonLakeModel\ModCastelo do Bode\src\2_Temp_model_version2.py", line 384, in <module>
    Thermocline_Depth[i-1] = ThermoTemp(TempLake[i-1], z, Smin=0.1, seasonal=True, index=False, mixed_cutoff=1)
  File "C:\Users\manue\Google Drive\Lake model\PythonLakeModel\ModCastelo do Bode\src\2_Temp_model_version2.py", line 169, in ThermoTemp
    Sdn = -(depths[thermoInd+1] - depths[thermoInd])/(drho_dz[thermoInd+1] - drho_dz[thermoInd])
IndexError: index 40 is out of bounds for axis 0 with size 40