## Outgoing longwave radiation.
## Comparre longwave radiation from the surface (2-metre temperature) and top of the atmosphere
library(esd)
t2m <- retrieve('~/data/ERAINT_t2m.nc')
## https://www.esrl.noaa.gov/psd/data/gridded/data.olrdcdr.interp.html
olr <- retrieve('~/data/olr.mon.mean.nc')

sigma <- 5.67e-8
coredata(t2m) <- sigma*(273.15+coredata(t2m))^4
t2m <- regrid(t2m,is=olr)
attr(t2m,'unit') <- 'W/m^2'
attr(t2m,'variable') <- 'LW'
lws <- map(t2m)
lwta <- map(olr)
lwdiff <- map(lws - lwta)

## Plot the data
breaks <- seq(100,500,by=25)
par(mfcol=c(3,1),mar=c(1.5,3,1.5,5))
data(geobordrs)
image(lon(lws),lat(lws),coredata(lws),ylab="",xlab="",
      main="Surface long-wave radiation",breaks=breaks,col=colscal(length(breaks)-1))
lines(geoborders)
image(lon(lwta),lat(lwta),coredata(lwta),ylab="",xlab="",
      main="Outgoing long-wave radiation",breaks=breaks,col=colscal(length(breaks)-1)) 
lines(geoborders)
image(lon(lwdiff),lat(lwdiff),coredata(lwdiff),ylab="",xlab="",
      main="Surface-top difference",breaks=seq(0,300,by=10),col=colscal(30)) 
lines(geoborders)

par(new=TRUE,fig=c(0,1,0,1))
colbar(breaks=breaks,col=colscal(length(breaks)-1),fig=c(0.96,0.98,0.7,0.85))
colbar(breaks=breaks,col=colscal(length(breaks)-1),fig=c(0.96,0.98,0.35,0.50))
colbar(breaks=seq(0,300),col=colscal(30),fig=c(0.96,0.98,0.05,0.20))