Python: creating a graph of a module that contains a feedback mechanism
I am new to programming and I am trying to create a simple zero energy balance model in Python 2.7 IDLE to calculate the Earth's surface temperatures and add ice albedo feedback, i.e. if the temperature output from the model above 280 K albedo remains at 0.3 (reflected from 30% of the energy), if it is below 250 thousand albedo is 0.7 (70% of the energy is reflected, since its cooler, therefore, larger ice ( white) covers the Earth) and if the temperature is in the range between them; albedo is calculated by the formula. This new albedo value is then returned from the model to obtain a more accurate temperature.
In my module, I have defined;
Final climate model Calculation of albedo New final climate model with new albedo (s) taken into account
I am trying to create a graph to compare the results of a first climate model with a variable input to solar energy, but with a consistent albedo, with the output of a second run with a variable albedo and solar output. But keep getting errors;
This is my script for my graph:
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 ()
My error message:
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()
I believe it refers to something in this part of my module:
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
Any help is appreciated
thank
source to share
The remark above is correct and it is not clear what you want to do, but if you want to check if all the elements in your array are checking a condition, then you could do:
if tb.all() > 280.0:
If you're wondering if a value exists in the array that fills it, you can do:
if tb.max() > 280.0:
...
elif tb.min() < 250.0:
Both of the above examples shouldn't be more than a simple else statement for the third condition.
If you want to evaluate the positions separately, you could as well, but then I would go for the following:
tb_test = np.ones(tb.shape) * 3
tb_test[np.where(tb > 280)] = 1
tb_test[np.where(tb < 250)] = 2
This will make the array tb_test for those where the first condition applies, two for the second and three for the third.
Of course, you can embed your calculations directly instead of the above identification of where different conditions apply ...
source to share