First, we load the forecast of the target variable recursively for lead month values from 1 to 6 (i.e., the initializations from January to June)

> cfs.list <- lapply(1:6, function(lead.month) {
      loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", lonLim = c(-10,-1), latLim = c(36,40), season = 7, years = 2006, leadMonth = lead.month)
})

We will take as reference for measuring the drift the lead-month 0 forecast (i.e. the July initializations)

> ref <- loadECOMS(dataset = "CFSv2_seasonal_16", var = "tas", lonLim = c(-10,-1), latLim = c(36,40), season = 7, years = 2006, leadMonth = 0)
> plotMeanField(ref)