## Zonal energy balance model of the Earth's climate system. ## Copyright (c) 2006 Hugh C. Pumphrey and released under the terms of the ## GNU General Public License version 2 as published by the Free ## Software Foundation: see http://www.gnu.org/copyleft/gpl.txt ## The model is implemented exactly as described by Henderson-Sellers ## and McGuffie and described in the course notes for the Physics of ## Climate year 4 Physics option at the University of Edinburgh. library(tcltk) ## Make initialtemp a global variable initialtemp <- 0 ## Implement the Henderson-Sellers & McGuffie EBM in R ebmcalc <- function(){ ## Model latitudes nbins <- as.numeric(tclvalue(tkget(nbscale.widget))) binwidth <- 90.0/nbins latitudes <- seq(binwidth/2,(90-binwidth/2),length=nbins) coslat <- cos(latitudes*pi/180) sumcoslat <- sum(coslat) reuseval <- tclvalue(reuse) if(reuseval == "1" && length(initialtemp) == nbins){ print("Starting with result of last run") }else{ initialtemp <- seq(15,-10,length=nbins) } ## Tweakables solar.constant <- as.numeric(tclvalue(tkget(scscale.widget))) critical.temp1 <- -10 critical.temp2 <- -9 ice.albedo <- as.numeric(tclvalue(tkget(iascale.widget))) nonice.albedo <- as.numeric(tclvalue(tkget(ascale.widget))) ## These three could have more useful names. They are called ## A, B and K in the original ## A and B are parameters in a linear approx for OLR: ## OLR = A + BT, where T is in CELSIUS, not K constA <- 204 ## Default 204 constB <- 2.17 ## Default 2.17 ##Meridional transport. Default 3.87 constK <- as.numeric(tclvalue(tkget(kscale.widget))) ## Calculate the weighting that we apply to solarConst/4 ## to get the right shortwave input for a latitude bin. ## Annoyingly, HSMcG only give values, not a formula. ## We interpolate between their values. This should be replaced ## by the correct formula worked out from the geometry. flat <- c(-1,5,15,25,35,45,55,65,75,85,91) swdat <- c(1.219,1.219,1.189,1.12,1.021,0.892,0.77,0.624,0.531,0.5,0.5) sc.fraction <- 1 rad.in <- approx(flat,swdat,latitudes)$y rad.in <- rad.in*(solar.constant/4)*sc.fraction convcrit <- 1.e-12 wktemp <- initialtemp initial.albedo <- latitudes*0 + ice.albedo initial.albedo[initialtemp > critical.temp2] <- nonice.albedo ix <- initialtemp > critical.temp1 & initialtemp < critical.temp2 if (length(ix) > 0){ initial.albedo[ix] <- ice.albedo + (initialtemp[ix]-critical.temp1)* (nonice.albedo - ice.albedo) / (critical.temp2-critical.temp1) } final.albedo <- initial.albedo finaltemp <- initialtemp niters <- 200 temp.iters <- array(0,c(niters,nbins)) for(iter in 1:niters){ wktemp <- finaltemp temp.iters[iter,] <- wktemp globalmeantemp <- sum(wktemp*coslat) / sumcoslat ## print("coslat=") ## print(coslat) ## print("wktemp=") ## print(wktemp) ## print("wktemp*coslat=") ## print(wktemp*coslat) rad.absorbed <- rad.in*(1-final.albedo) finaltemp <- ( rad.absorbed + constK*globalmeantemp - constA)/ (constB+constK) final.albedo <- initial.albedo*0+ice.albedo final.albedo[(finaltemp > critical.temp2)] <- nonice.albedo ix <- finaltemp > critical.temp1 & finaltemp < critical.temp2 if (length(ix) > 0){ final.albedo[ix] <- ice.albedo + (finaltemp[ix]-critical.temp1)* (nonice.albedo - ice.albedo) / (critical.temp2-critical.temp1) } dt <- finaltemp - wktemp convcrit <- t(dt) %*% dt ## print(paste(iter,globalmeantemp,convcrit)) ##print("----------") if(convcrit < 1.e-6) break } print(paste("Iterations:",iter," Global mean T=",signif(globalmeantemp,5), "C, ",signif(globalmeantemp+273.16,5),"K")) ## How can we get some text to appear in the control window? ## Add T to label -- This doesn't work. The tcl var gets set OK, ## but the label widget never reflects it tclvalue(labeltext)<-paste("EBM: Global mean T=",signif(globalmeantemp,5)) ##print(paste("label should be",tclvalue(labeltext))) ## Put T into text widget ##tkinsert(meantemp.widget,paste(signif(globalmeantemp+273.16 ,5))) radout <- constA + constB*finaltemp m <- matrix(c(1,2,3),ncol=1) layout(m,heights=c(1,0.5,1)) par(mai=c(0.7,0.7,0.05,0.01)) par(ps=18) ## Plot out results ## Temperature plot(latitudes,finaltemp,xlab="Latitude",ylab="Temperature / C",type="l") for(i in 1:iter){ lines(latitudes,temp.iters[i,],col=grey(0.8-0.5*i/iter),lty="dashed") } points(latitudes,finaltemp) points(latitudes,initialtemp,col="red",pch=2) lines(latitudes,initialtemp,col="red") abline(h=globalmeantemp,lty="dotted") legend("topright",c("Initial T","Final T"),col=c("red","black"), lty="solid",pch=c(2,1) ) ## Albedo plot(latitudes,final.albedo,type="o",xlab="Latitude",ylab=expression(alpha)) ## Fluxes plot(latitudes,radout,ylim=range(c(radout,rad.absorbed,rad.in,rad.absorbed - radout)),type="o", xlab="Latitude",ylab="Flux / Wm2") points(latitudes,rad.in,col="red",pch=2) lines(latitudes,rad.in,col="red") points(latitudes,rad.absorbed,col="green",pch=3) lines(latitudes,rad.absorbed,col="green") points(latitudes,rad.absorbed-radout,col="blue",pch=4) lines(latitudes,rad.absorbed-radout,col="blue") ## Test: should be the same as blue line above. It is. ## fluxmer <- constK*(finaltemp - globalmeantemp) ## points(latitudes,fluxmer,col="magenta",pch=5) ## lines(latitudes,fluxmer,col="magenta") abline(h=0,lty="dotted") legend("topright",c("Outgoing Flux","Incoming Flux", "Absorbed Flux","Abs - Out"), col=c("black","red","green","blue"),lty="solid",pch=c(1,2,3,4)) ## last thing: set initial temp to final temp so it can be used as ## the starting point for the next run initialtemp <<- finaltemp } ## ------------ Main program --------------------------------- ## Start up a Tk window tt <- tktoplevel() tktitle(tt) <- "Energy balance Model" bframe <- tkframe(tt) slframe <- tkframe(tt) labeltext <- 33 ## For some weird reason this has to be numeric, even ## though the value is never used. And it has to be a _different_ value to ## any of the other Tcl variables you use, or they become the _same_ ## Tcl variable. Very flaky! ## Generate the things we want to put in the window tclvalue(labeltext)<-"EBM" label.widget <- tklabel(tt,textvariable=labeltext) ##print(paste("label should be",tclvalue(labeltext))) ## Why can't we use command=ebmcalc to re-do the model every time the ## slider is moved? ## Bigincrement appears to do nothing, because the Tk window ignores the kbd ## (except for text widgets) scscale.widget <- tkscale(slframe,from=500,to=2800,label="Solar Const", resolution=5, length=200,orient="horizontal", bigincrement=20,takefocus=1) iascale.widget <- tkscale(slframe,from=0,to=1,label="Ice Albedo", resolution=0.05, length=200,orient="horizontal") ascale.widget <- tkscale(slframe,from=0,to=1,label="Non-ice Albedo", resolution=0.05, length=200,orient="horizontal") kscale.widget <- tkscale(slframe,from=0,to=15,label="Meridional transport", resolution=0.1, length=200,orient="horizontal") nbscale.widget <- tkscale(slframe,from=9,to=50,label="Number of bins", resolution=1, length=200,orient="horizontal") button.widget <- tkbutton(bframe, text="Run EBM",command=ebmcalc) ## meantemp.widget <- tktext(bframe,width=5,height=1) ## As noted above, this has to be numeric, and be a different value from ## other tcl variables. But the value as an R variable is never used or changed reuse <- 66 reuse.widget <- tkradiobutton(bframe,text="start from previous run",variable=reuse,value=1) dontreuse.widget <- tkradiobutton(bframe,text="default start", variable=reuse,value=2) nuke.widget <- tkbutton(bframe, text="Exit EBM", command=function()tkdestroy(tt)) ## how we set a widget to start at a specific value settodef <- function(){ tkset(scscale.widget,"1370") tkset(iascale.widget,"0.6") tkset(ascale.widget,"0.3") tkset(kscale.widget,"3.87") tkset(nbscale.widget,"31") tclvalue(reuse)<-2 ### Default start } def.widget <- tkbutton(bframe, text="Default Settings", command=settodef) settodef() ## Run with default vals ebmcalc() ## Pack the things into the window tkpack(button.widget) tkpack(def.widget) tkpack(nuke.widget) tkpack(reuse.widget) tkpack(dontreuse.widget) ###tkpack(meantemp.widget) tkpack(nbscale.widget,scscale.widget,iascale.widget, ascale.widget,kscale.widget) tkpack(label.widget,side="top") tkpack(bframe,side="left") tkpack(slframe,side="right")