library(ncdf4)

source("r/image.plot.r")

nc.1 <- nc_open("ERA_fluxes_00-12.nc")
nc.2 <- nc_open("ERA_fluxes_12-24.nc")

lon <- ncvar_get(nc.1, "longitude")
lat <- ncvar_get(nc.1, "latitude")
time <- ncvar_get(nc.1, "time") * 3600 + as.POSIXct("1900-01-01 UTC", tz = "UTC")

## incoming solar radiation
tisr <- (ncvar_get(nc.1, "tisr") + ncvar_get(nc.2, "tisr")) / (24 * 3600)

sr.global <- apply(tisr, 3, ## time series
                   function(tisr)
                   sum(apply(tisr, 2, mean) ## zonal mean
                       * cos(lat * pi / 180) * 2 * pi / 180) ## cos phi dphi
                   / 2) ## normalization (\int_{\pi/2}^{\pi/2} \cos\phi \,d\phi)

plot(time, sr.global, type = "l") ## plot the effect of orbital eccentricity on the incoming solar radiation

plot(lat, apply(tisr, 2, mean), type = "l") ## plot the meridional distribution of the incoming solar radiation