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

[Yoriz]

The function fix_profile must be altering something that the rest of the code is dependant on.
No one else can understand the code that has been presented.

[mcva]

The function fix_profile must be altering something that the rest of the code is dependant on.
No one else can understand the code that has been presented.

---

Its not altering nothing,
I can change the name of the variable CANT_REMOVE_THIS_LINE=fix_profile(n,TempTDMA,Areat0) to ffffffff=fix_profile(n,TempTDMA,Areat0) ant the model runs, but when I remove the entire line its doesn´t run!

this is the full code:

import numpy as np
from cmath import pi
np.set_printoptions(threshold=np.inf)
import pylab as plt
import math
from importfile import air,neb,wind,Z,z,Areat0,Areat1,T0


#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#INPUT
#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


RelativeHumidity = np.full((2600),80)


#------------------------------------------------------------------------------------------------------------------------------------
# Tri Diagonal Matrix Algorithm (Thomas algorithm) solver (source:https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9)
#-------------------------------------------------------------------------------------------------------------------------------------

def TDMAsolver(a, b, c, d):
    nf = len(d)
    ac, bc, cc, dc = map(np.array, (a, b, c, d))
    for it in xrange(1, nf):
        mc = ac[it-1]/bc[it-1]
        bc[it] = bc[it] - mc*cc[it-1] 
        dc[it] = dc[it] - mc*dc[it-1]
    xc = bc
    xc[-1] = dc[-1]/bc[-1]
    for il in xrange(nf-2, -1, -1):
        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]
    del bc, cc, dc

    return xc

#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Clear Sky radiation (S.L.Dingman, 2012 - Physical Hydrology, 2nd Ed.)
#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

def ClearSky(i, Lat, AirTemperature, RelativeHumidity):
    Latitude_radians = (Lat/360)*2*pi    
    SolarConstant = 4.921 #MJ/m2hr
    DustAttenuation = 0.05
    OpticalAirMass = 3.6
    AirT = AirTemperature
    RH = RelativeHumidity
    SurfaceAlbedo = 0.4
    VapourPressure_ea = 0.611*np.exp(17.3*AirT/(AirT+237.3))*RH/100
    DewPt = (np.log(VapourPressure_ea)+0.4926)/(0.0708-0.00421*np.log(VapourPressure_ea))
    PrecipitableWater = 1.12*np.exp(0.0614*DewPt)
    Tsa = np.exp((-0.124-0.0207*PrecipitableWater)+(-0.0682-0.0248*PrecipitableWater)*OpticalAirMass)
    T = Tsa-DustAttenuation
    Ts = np.exp((-0.0363-0.0084*PrecipitableWater)+(-0.0572-0.0173*PrecipitableWater)*OpticalAirMass)
    gs = 1-Ts+DustAttenuation   
    DayAngle = 2*pi*(i-1)/365
    Declination_degrees = (180/pi)*(0.006918-0.399912*np.cos(DayAngle)+0.070257*np.sin(DayAngle)-0.006758*np.cos(2*DayAngle)+0.000907*np.sin(2*DayAngle)-0.002697*np.cos(3*DayAngle)+0.00148*np.sin(3*DayAngle))
    Declination_radians = (Declination_degrees/360)*2*pi
    Eccentricity = 1.00011+0.034221*np.cos(DayAngle)+0.00128*np.sin(DayAngle)+0.000719*np.cos(2*DayAngle)+0.000077*np.sin(2*DayAngle)
    Sunrise = -np.arccos(-np.tan(Latitude_radians)*np.tan(Declination_radians))/0.2618
    Sunset = np.arccos(-np.tan(Latitude_radians)*np.tan(Declination_radians))/0.2618
    DayLength = 2*Sunset
    Extraterrestrial = 2*SolarConstant*Eccentricity*(np.cos(Declination_radians)*np.cos(Latitude_radians)*np.sin(0.2618*Sunset)/0.2618+np.sin(Declination_radians)*np.sin(Latitude_radians)*Sunset)#MJ/m2hr
    Direct = Extraterrestrial*T
    Diffuse = 0.5*gs*Extraterrestrial
    Global = Direct+Diffuse
    Backscatter = 0.5*gs*SurfaceAlbedo*Global
    ClearSky = Global+Backscatter # MJ/m2day
    ClearSkyRadiaton = ClearSky*11.575
    
    return ClearSkyRadiaton

#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# EQUILILBRIUM TEMPERATURE
#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

def EquilibriumTemp(AirTemperature, RelativeHumidity, wind, ClearSkyRadiaton, WaterTemperature):
    AirT_C = AirTemperature
    AirT_F = AirT_C*1.8+32
    WaterTemp_C = WaterTemperature
    WindZ_ms = wind
    WindZ_mph = WindZ_ms*2.23714
    Wind2_mph = WindZ_mph*np.log(2/0.003)/np.log(10/0.003)
    SCR_Wm2=ClearSkyRadiaton
    SCR_BTU_FT2_DAY = SCR_Wm2*7.60796
    RH = RelativeHumidity
    VapourPressure_ea = 0.611*np.exp(17.3*AirT_C/(AirT_C+237.3))*RH/100
    DewPt = (np.log(VapourPressure_ea)+0.4926)/(0.0708-0.00421*np.log(VapourPressure_ea))
    DewPt_F = DewPt*1.8+32
    Aconv = 7.60796
    Bconv = 1.520411 #IF(CFW =2,1.520411,3.401062)
    AFW = 9.2
    BFW = 0.46
    CFW = 2
    
    ET = DewPt_F
    TSTAR = (ET+DewPt_F)*0.5
    BETA  =  0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
    FW = Aconv*AFW+Bconv*BFW*Wind2_mph**CFW
    CSHE = 15.7+(0.26+BETA)*FW
    RA = 3.1872E-08*(AirT_F+459.67)**4
    ETP = (SCR_BTU_FT2_DAY+RA-1801.0)/CSHE+(CSHE-15.7)*(0.26*AirT_F+BETA*DewPt_F)/(CSHE*(0.26+BETA))

    j = 0
    while (abs(ETP-ET) > 0.05 or j<10):
        ET = ETP
        TSTAR = (ET+DewPt_F)*0.5
        BETA = 0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
        CSHE = 15.7+(0.26+BETA)*FW
        ETP = (SCR_BTU_FT2_DAY+RA-1801.0)/CSHE+(CSHE-15.7)*(0.26*AirT_F+BETA*DewPt_F)/(CSHE*(0.26+BETA))
        EquilibriumTemperature = (ET-32)*5/9
        SurfaceHeatExchange = CSHE*0.23659
        Qn = (EquilibriumTemperature-WaterTemp_C)*SurfaceHeatExchange
        j+=1
        
    return Qn


#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# DEPTH OF THE THERMOCLINE (source:https://rdrr.io/cran/rLakeAnalyzer/src/R/thermo.depth.R)
# This function calculates the location of the thermocline from a temperature profile.  
# It uses a special technique to estimate where the thermocline lies even between two temperature measurement depths, giving a potentially finer-scale estimate than usual techniques.
# @param wtr a numeric vector of water temperature in degrees C
# @param depths a numeric vector corresponding to the depths (in m) of the wtr measurements
# @param Smin Optional paramter defining minimum density gradient for thermocline
# @param seasonal a logical value indicating whether the seasonal thermocline should be returned. This is fed to thermo.depth, which is used as the starting point.  The seasonal thermocline is defined as the deepest density gradient found in the profile. If \code{FALSE}, the depth of the maximum density gradient is used as the starting point.
# @param index Boolean value indicated if index of the thermocline depth, instead of the depth value, should be returned.
# @param mixed.cutoff A cutoff (deg C) where below this threshold, thermo.depth and meta.depths are not calculated (NaN is returned). Defaults to 1 deg C.
# @return Depth of thermocline. If no thermocline found, value is NaN.
#----------------------------------------------------------------------------------------------------------------------------------------

def find_peaks(data_in, thresh = 0.0): #Finds the local peaks in a vector. Checks the optionally supplied threshold for minimum height.
    peaks = [thresh - 1]
    for trip in zip(*[data_in[n:] for n in range(3)]):
        if max(trip) == trip[1]:
            peaks.append(trip[1])
        else:
            peaks.append(thresh - 1)
    return [index for index, peak in enumerate(peaks) if peak > thresh]


def ThermoTemp (wtr, depths, Smin=0.1, seasonal=True, index=False, mixed_cutoff=1):

    rhoVar = np.zeros(len(depths)) #23
       
    for i in range(len(depths)): #0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22
        rhoVar[i]=(1000*(1-(wtr[i]+288.9414)*(wtr[i]-3.9863)**2/(508929.2*(wtr[i]+68.12963))))
        
    #Calculate the first derivate of density
    drho_dz = np.zeros(len(depths)-1)
    for i in range(len(depths)-1): #0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21
        drho_dz[i] = ( rhoVar[i+1]-rhoVar[i] )/( depths[i+1] - depths[i] )
        
    #look for two distinct maximum slopes, lower one assumed to be seasonal 
    thermoInd = np.argmax(drho_dz) # return index of max slope 9
    mDrhoZ = drho_dz[thermoInd] #return max density value 0.99
    thermoD = np.mean(depths[(thermoInd):(thermoInd+2)])# mean value depth (entre 7 e 8)
    #                            9:11
    dRhoPerc = 0.15
    numDepths = len(depths)
                                       #23
    if thermoInd > 0 and thermoInd < (numDepths):

        Sdn = -(depths[thermoInd+1] - depths[thermoInd])/(drho_dz[thermoInd+1] - drho_dz[thermoInd])
        Sup = (depths[thermoInd]-depths[thermoInd-1])/(drho_dz[thermoInd]-drho_dz[thermoInd-1])

        upD  = depths[thermoInd]
        dnD  = depths[thermoInd+1]
 
        if np.isinf(Sup) or np.isinf(Sdn):
            thermoD = dnD*(Sdn/(Sdn+Sup))+upD*(Sup/(Sdn+Sup))
      
    dRhoCut = np.max([dRhoPerc*mDrhoZ,Smin])
    locs = find_peaks(drho_dz, dRhoCut)
    pks = drho_dz[locs]
    
   
    if (len(pks)==0):
        SthermoD = thermoD
        SthermoInd = thermoInd
        
    else:
        mDrhoZ = pks[len(pks)-1]
        SthermoInd = locs[len(pks)-1] #returns 13 instead of 14
        if (SthermoInd > (thermoInd + 1)):
            SthermoD = np.mean(depths[(SthermoInd):(SthermoInd+2)])
            
            if (SthermoInd > 0 and SthermoInd < (numDepths)):
                Sdn = -(depths[SthermoInd+1] - depths[SthermoInd])/(drho_dz[SthermoInd+1] - drho_dz[SthermoInd])#4.35
                Sup = (depths[SthermoInd] - depths[SthermoInd-1])/(drho_dz[SthermoInd] - drho_dz[SthermoInd-1])#6.34
                upD  = depths[SthermoInd]#11
                dnD  = depths[SthermoInd+1] #12
                
                if np.isinf(Sup) or np.isinf(Sdn):
                    SthermoD = dnD*(Sdn/(Sdn+Sup))+upD*(Sup/(Sdn+Sup))      
        else:
            SthermoD = thermoD
            SthermoInd = thermoInd
                    
    if SthermoD < thermoD:
        SthermoD = thermoD
        SthermoInd = thermoInd
         
#Ok, which output was requested. Index or value, seasonal or non-seasonal

    if(index):
        if(seasonal):
            return(SthermoInd)
        else:
             return(thermoInd)
    
    else:
        if(seasonal):
            return(SthermoD)
        else:
            return(thermoD)


#-------------------------------------------------------------------------------------------------------------------------------------
#EDDY DIFFUSIVITY
#-------------------------------------------------------------------------------------------------------------------------------------

def Keddy(T0,z,v, Areat0):
    g=9.8
     
    #Temperature gradient
    Gtemp=np.zeros(len(T0))
    for i in range(len(T0)-1): #0-39 (40)
        Gtemp[i]=(T0[i+1]-T0[i])/(z[i+1]-z[i])
        Gtemp[len(T0)-1]=Gtemp[len(T0)-2]
    
    #Thermal expansion of water
    Texp=np.zeros(len(T0))
    for i in range(len(T0)-1):
        if T0[i]<=11:
            Texp[i]=(1.43*10**-5)*(T0[i])-(0.54*10**-4)
        else:
            Texp[i]=(1.41*10**-5)*(T0[i])-(0.209*10**-4)
    
    #Friction velocity
    w=1.2*(10**-3)*v
    c=0.1
    d=0.1
    
    #Eddy diffusivity in the absence of stratification, m2/s 
    Kh0=c*w
    
    #Richardson number
   
    Ri=np.zeros(len(T0))
    for i in range(len(T0)):
            Ri[i]=-((Texp[i]*g*(z[i]**2.0))/(w**2.0))*(Gtemp[i])

    #Eddy diffusivity for the vertical transport of heat
    KH=np.zeros(len(T0))
    for i in range(len(T0)):
        KH[i]=Kh0*(1+d*Ri[i])**-1

    #KH[1:]=np.where(KH[1:] <=0.00001, 0.00001, KH[1:])
    
    Kdiffusivity=[a*b for a,b in zip(Areat0,KH)]
    
    index_thermocline = ThermoTemp (T0, z, Smin=0.1, seasonal=False, index=True, mixed_cutoff=1)
    
    for n, i in enumerate(Kdiffusivity):
        if n > index_thermocline:
            Kdiffusivity[n]=Kdiffusivity[index_thermocline]
    return Kdiffusivity

#-------------------------------------------------------------------------------------------------------------------------------------
#FIX PROFILE
#-------------------------------------------------------------------------------------------------------------------------------------

def fix_profile(n,T, area):
    var1=range(1,n)
    var1.reverse()
    
    for i in var1:
        for j in range(n-1):
            if T[i]>T[j]:
                change_profile(T,i,j,area,n)
                break
    return T         
 
def change_profile(T,i,inst,area,n):
    jmin=inst-1
    if jmin<0:
        jmin=0
        min=1
        jmin0=jmin
    else:
        min=0
        jmin0=inst
        
    jmax=i+1
    if jmax>n:
        jmax=n
        max=1
        jmax0=jmax
        
    else:
        max=0
        jmax=i
    
    A=0.0
    B=0.0
    for j in range(jmin,jmax+1):
        A+=T[j]*area[j]
        B+=area[j]
        Tconst=A/B
    
    for j in range(jmin,jmax+1):
        T[j]=Tconst


#-------------------------------------------------------------------------------------------------------------------------------------
#Crank Nicolson scheme
#-------------------------------------------------------------------------------------------------------------------------------------

#Simulation time period (days)
N=2556
#2556
#Latitude (degrees)
Lat=38.0

#-------------------------------------------------------------------------------------------------------------------------------------
#Run Clear Sky Radiation
ClearSkyRadiation = np.zeros(N)
I = np.linspace(1,N,N)

for i in xrange(N):
    ClearSkyRadiation[i] = ClearSky(i, Lat, air[i], RelativeHumidity[i])
    

plt.plot(I,ClearSkyRadiation)
plt.show()
#-------------------------------------------------------------------------------------------------------------------------------------


Nlayers=41 #41

#Space, m
DZ=2.0

#Water temperature
TempLake=np.zeros((N+1,Nlayers))
#K=np.zeros((N+1,Nlayers))
K=np.full((N+1,Nlayers), 0.0000)

# Initial condition , C
TempLake[0]=T0

#time interval, seconds, s
DT=86400
#Water density, kg/m3
DH2O=1000.0
#Water specific Heat j/kg C
Cp=4186

s=DT/(2*DZ**2.0)

#Initializing coefficients np.arrays
a=np.zeros((N,Nlayers-1)) #40
b=np.zeros((N,Nlayers)) #41
c=np.zeros((N,Nlayers-1)) 
d=np.zeros((N,Nlayers))
alfa = np.zeros(Nlayers)#41
Radiation = np.zeros(N)
Thermocline_Depth = np.zeros(N)

for i in xrange(Nlayers):
    alfa[i]=1.0/(Areat0[i]+Areat0[i])
    


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)
    
    #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




#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
#OUTPUT
#-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


Kout=np.full((N+1,Nlayers), 0.0000)

for i in xrange(1,N+1):
    Kout[i] = [a/b for a,b in zip(K[i],Areat0)]

np.savetxt('Kdifussivity.out', Kout, delimiter=',')


plt.plot(I,Radiation)
plt.show()  
plt.plot(I,Thermocline_Depth)
plt.show()  


print K

Z=TempLake 
 
X,Y=np.meshgrid(range(Z.shape[0]+1),range(Z.shape[1]+1)) 

plt.figure(1)
plt.subplot(211)
im = plt.pcolormesh(X,Y,Z.transpose(), cmap='jet') 
ax = plt.gca()

ax.set_ylim(ax.get_ylim()[::-1])

plt.colorbar(im, orientation='horizontal') 

plt.ylabel('Profundidade, m')

plt.subplot(212)
Zout=Z.transpose()
plt.plot(np.linspace(1,N+1,N+1), Zout[0], lw=2)
plt.xlabel('Dias')
plt.ylabel('Temperatura, C')

plt.show()