library(tcltk) initialtemp <- 0 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("Launch based on the previous simulation's results") }else{ initialtemp <- seq(15,-10,length=nbins) } ## Tweakables solar.constant <- as.numeric(tclvalue(tkget(scscale.widget))) critical.temp1 <- as.numeric(tclvalue(tkget(icscale.widget))) critical.temp2 <- critical.temp1-4 ice.albedo <- as.numeric(tclvalue(tkget(iascale.widget))) nonice.albedo <- as.numeric(tclvalue(tkget(ascale.widget))) ## A, B oraz K w oryginale ## A oraz B - 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("Iteration number:",iter," Global avg. temp: 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 avg. temp: 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") grid(col='black') 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","Computed 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="Energy flux / Wm2") grid(col='black') 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 of energy", "Incoming flux of energy", "Energy absorbed","Abs - outgoing"), 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) <- "1-Dimensional 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=1000,to=1600,label="Solar constant (W/m^2)", resolution=0.1, length=800,orient="horizontal", bigincrement=5,takefocus=0.25) iascale.widget <- tkscale(slframe,from=0,to=1,label="Albedo for ice surface", resolution=0.01, length=800,orient="horizontal") ascale.widget <- tkscale(slframe,from=0,to=1,label="Albedo for non-ice surface", resolution=0.01, length=800,orient="horizontal") kscale.widget <- tkscale(slframe,from=0,to=15,label="Latitudinal transport of energy", resolution=0.01, length=800,orient="horizontal") icscale.widget <- tkscale(slframe,from=-15,to=1,label="Sea temperature threshold for becoming ice", resolution=0.1, length=800,orient="horizontal") nbscale.widget <- tkscale(slframe,from=9,to=90,label="Number of latitude bands", resolution=1, length=800,orient="horizontal") button.widget <- tkbutton(bframe, text="LAUNCH \nSIMULATION",command=ebmcalc) meantemp.widget <- tktext(bframe,width=6,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="launch simulation from previously obtained results",variable=reuse,value=1) dontreuse.widget <- tkradiobutton(bframe,text="start", variable=reuse,value=2) nuke.widget <- tkbutton(bframe, text="Exit model", 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,"90") tkset(icscale.widget,"-4") tclvalue(reuse)<-2 ### Default start } def.widget <- tkbutton(bframe, text="Restore parameters to default", 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, icscale.widget) tkpack(label.widget,side="top") tkpack(bframe,side="left") tkpack(slframe,side="right")