2012-03-17 32 views
3

我對編程相當陌生,我試圖在Python 2.7 IDLE中生成一個簡單的零維能量平衡模型,以計算地球表面溫度並增加了冰反照反饋,即如果模型的溫度輸出高於280K,則反照率保持在0.3(反射30%的能量),如果其低於250k,則反照率爲0.7(反射70%的能量,因爲其更冷卻器地球上有更大的冰蓋(白色)),以及溫度是否處於這兩者之間的範圍內;根據公式計算反照率。反照率的這個新值然後從模型返回,以提供更準確的溫度。Python:生成一個包含反饋機制的模塊的圖形

在我的模塊中,我已經定義了;

最後一個氣候模型 計算反照率 新敲定的氣候模型考慮到concideration新的反照率(S)

我想製作一個圖形的首個氣候模型的輸出與變化的比較太陽能輸入,但一致的反照率,第二輪輸出與不同的反照率和太陽能輸出。但不斷收到錯誤;

這是我爲我的圖形腳本:

import matplotlib.pyplot as plt 
    import numpy as np 
    from EBM_IceAlbFeedback import * 
    # q is for the Solar Constant 
    q=np.linspace(2.5e26,4.0e26,150) 
    # t= temperature derived from the final climate model 
    t= finalCM(Q=q) 
    plt.plot(q,t,'b-') 
    q=np.linspace(3.0e26,4.5e26,150) 
    # tb= is the second set of temperatures derived from the NEWfinalCM which contains an Ice Albedo Feedback 
    tb= NEWfinalCM(Q=q) 
    plt.plot(q,tb,'r-') 
    plt.show() 

我的錯誤信息是:

Traceback (most recent call last): 
File "K:/python/CompareCMsPlt2.py", line 13, in <module> 
tb= NEWfinalCM(Q=q) 
File "K:/python\EBM_IceAlbFeedback.py", line 228, in NEWfinalCM 
NewAlb=NAlb(dist=dist, Q=Q, co2Emissions=co2Emissions, alpha=alpha, cCycleInt=cCycleInt, cCycleSlope=cCycleSlope) 
File "K:/python\EBM_IceAlbFeedback.py", line 190, in NAlb 
    if ta>280.0: 
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all() 

我相信這是我的模塊,這部分提到的東西:

def NAlb (dist=150e9, Alb=0.3, Q=3.87e26, co2Emissions=0.0, alpha=3.0, cCycleInt=0.4, cCycleSlope=0.0001): 
''' 
Readjusting Albedo to the output temperature 

Arguments: 

Q = solar ouput (W) 
dist = distance from the sun (m) 
co2Emissions = Cumulative CO2 emissions since 2010 (GtC) 
alpha = climate sensitivity (K/2xCO2) 
cCycleInt = Initial value of the airborne fraction (unitless) 
cCycleSlope = Increment the airborne fraction per GtC (GtC^-1) 

Return Value: 
NewAlb= New Albedo (Unitless) 
''' 
# CALCULATE ABORTIVITY: 
#Our model is baselined at an atmospheric CO2 concentration of 390 ppmv in 2010 
baselineCO2=390.0 
#The official IPCC figure for conversion of mass of emissions (GtC) top atmospheric concentration (ppmv) 
IPCCmassToConc=2.12 
#approximate correction for the carbon cycle: 
cCycleAdjust=cCycleInt+cCycleSlope*co2Emissions 
#convert GtC to CO2 conc in ppmv: 
co2=co2Emissions*cCycleAdjust/IPCCmassToConc+baselineCO2 
#calculate absorptivity 
absrp=absrpFromCO2(CO2=co2, alpha=alpha) 

#CALCULATE TEMPERATURE: using the same method as in the finalCM 
ta=transATmCM (absrpt=absrp, dist=dist, Alb=0.3, Q=Q) 
# define the thresholds for an ice free state. 
if ta>280.0: 
    NewAlb=0.3 
# define the threshold for a snow ball Earth state. 
elif ta<250.0: 
    NewAlb=0.7# Calculate albedo for temperatures between 280k to 230k 
elif 250.0<ta<280.0: 
    NewAlb=(0.3+(((0.7-0.3)/(280.0-250.0))*(280.0-ta))) 
return NewAlb 




    def NEWfinalCM(co2Emissions=0.0, alpha=3., dist=150e9, Q=3.87e26, cCycleInt=0.4, cCycleSlope=0.0001): 
''' 
A New final Climate model which contains and Ice Albedo Feedback 

Arguments: 

Q = solar ouput (W) 
dist = distance from the sun (m) 
co2Emissions = Cumulative CO2 emissions since 2010 (GtC) 
alpha = climate sensitivity (K/2xCO2) 
cCycleInt = Initial value of the airborne fraction (unitless) 
cCycleSlope = Increment the airborne fraction per GtC (GtC^-1) 

Return Value: 
tn = surface temperature (K) 
''' 
#Our model is baselined at an atmospheric CO2 concentration of 390 ppmv in 2010 
baselineCO2=390.0 
#The official IPCC figure for conversion of mass of emissions (GtC) top atmospheric concentration (ppmv) 
IPCCmassToConc=2.12 
#approximate correction for the carbon cycle: 
cCycleAdjust=cCycleInt+cCycleSlope*co2Emissions 
#convert GtC to CO2 conc in ppmv: 
co2=co2Emissions*cCycleAdjust/IPCCmassToConc+baselineCO2 


#calculate temperature 
absrp=absrpFromCO2(CO2=co2, alpha=alpha) 
NewAlb=NAlb(dist=dist, Q=Q, co2Emissions=co2Emissions, alpha=alpha, cCycleInt=cCycleInt, cCycleSlope=cCycleSlope) 

tn=transATmCM(absrpt=absrp, dist=dist, Alb=NewAlb, Q=Q) 


return tn 

任何幫助表示讚賞

感謝

+1

@cillosis - 爲什麼?科學家不是程序員......至於實際的錯誤,numpy數組(例如'x> 5')的條件返回布爾數組(例如'array([True,True,False])'而不是單個值。 – 2012-03-17 00:51:53

回答

1

上面的評論是正確的,它並不清楚自己想要做什麼,但如果你想檢查你的數組中的所有元素驗證條件,那麼你可以這樣做:

if tb.all() > 280.0: 

如果你有興趣,如果有存在fullfills它在數組中的值,你可以這樣做:

if tb.max() > 280.0: 
    ... 
elif tb.min() < 250.0: 

兩個以上的例子不應該需要比簡單別的語句的第三個條件更多。

如果要單獨評估的位置,你可以爲好,但我會去以下:

tb_test = np.ones(tb.shape) * 3 
tb_test[np.where(tb > 280)] = 1 
tb_test[np.where(tb < 250)] = 2 

這將使第一個條件指tb_test陣列的,三三兩兩第二第三名是三分。

當然,你可以將你的計算,而不是直接的不同報考條件,其中上述鑑定...