(Mar-16-2019, 12:51 PM)Yoriz Wrote: 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()