diff --git a/inst/doc/plot_timeseries.R b/inst/doc/plot_timeseries.R
index 092cc5bf8b322f9c982bff717a240022df6d9d7b..40578f51bdb26adbcc261d78981cd7c0b80a242d 100755
--- a/inst/doc/plot_timeseries.R
+++ b/inst/doc/plot_timeseries.R
@@ -10,8 +10,8 @@ comprms <- TRUE
 compspread <- TRUE
 plotano <- TRUE
 
-var <- args[1]    # sie/sia/siv/tos/tas/prlr/ohc/lohc/mohc/uohc/amoc
-pole <- args[2]   # N/S only for sia/sie
+var <- args[1]    # siextent/siarea/sivol/tos/tas/prlr/ohcsum0/ohcsum0-17/ohcsum18-22/ohcsum23-46/vsftmyz
+pole <- args[2]   # N/S only for siarea/siextent
 nltimemax <- 124  # number of leadtimes max in the experiments (in months)
 nltimeout <- 60   # number of leadtimes to postprocess(in months)
 lstexpid <- c('i00k', 'b02p') # list of ids
@@ -23,28 +23,28 @@ yearf <- 2005     # last start date
 intsdate <- 5     # interval between start dates
 runmeanlen <- 12  # length of the window for running mean (in months)
 
-obs <- switch(var, 'sia' = c('hadisst_v1.1'), 'sie' = c('hadisst_v1.1'), 
+obs <- switch(var, 'siarea' = c('hadisst_v1.1'), 'siextent' = c('hadisst_v1.1'), 
               'tas' = c('ncep-reanalysis', 'era40'), 'tos' = c('ersst_v3b', 
                                                                'hadisst_v1.1'),
-              'prlr' = c('cru_v3.0', 'gpcc1x1_v6'), 'ohc' = c('nemovar_system4'),
-              'mohc' = c('nemovar_system4'), 'uohc' = c('nemovar_system4'),
-              'lohc' = c('nemovar_system4'), 'amoc' = c('nemovar_system4'),
-              'siv' = 'piomas')
-toptitle2 <- switch(var, 'sia' = "sea ice area", 'sie' = "sea ice extent",
-                    'siv' = "sea ice volume", 'tas' = "global T2m", 
+              'prlr' = c('cru_v3.0', 'gpcc1x1_v6'), 'ohcsum0' = c('nemovar_system4'),
+              'ohcsum18-22' = c('nemovar_system4'), 'ohcsum23-46' = c('nemovar_system4'),
+              'ohcsum0-17' = c('nemovar_system4'), 'vsftmyz' = c('nemovar_system4'),
+              'sivol' = 'piomas')
+toptitle2 <- switch(var, 'siarea' = "sea ice area", 'siextent' = "sea ice extent",
+                    'sivol' = "sea ice volume", 'tas' = "global T2m", 
                     'tos' = "global SST (60S-65N)", 
                     'prlr' = 'land precipitation (60S-65N)', 
-                    'ohc' = "global ocean heat content", 
-                    'lohc' = 'global 800m-bottom ocean heat content',
-                    'mohc' = 'global 350m-800m ocean heat content',
-                    'uohc' = 'global 0-350m ocean heat content',
-                    'amoc' = 'Atlantic Overturning Streamfunction (40-55N, 1-2km)'
+                    'ohcsum0' = "global ocean heat content", 
+                    'ohcsum0-17' = 'global 800m-bottom ocean heat content',
+                    'ohcsum18-22' = 'global 350m-800m ocean heat content',
+                    'ohcsum23-46' = 'global 0-350m ocean heat content',
+                    'vsftmyz' = 'Atlantic Overturning Streamfunction (40-55N, 1-2km)'
                     )
-ytitle1 <- switch(var, 'sia' = "Millions km2", 'sie' = "Millions km2", 
-                  'siv' = 'Thousands km3', 'tas' = 'K', 'tos' = 'K', 
-                  'prlr' = 'mm/day', 'ohc' = '10e22 Joules', 
-                  'lohc' = '10e22 Joules', 'mohc' = '10e22 Joules',
-                  'uohc' = '10e22 Joules', 'amoc' = 'Sv')
+ytitle1 <- switch(var, 'siarea' = "Millions km2", 'siextent' = "Millions km2", 
+                  'sivol' = 'Thousands km3', 'tas' = 'K', 'tos' = 'K', 
+                  'prlr' = 'mm/day', 'ohcsum0' = '10e22 Joules', 
+                  'ohcsum0-17' = '10e22 Joules', 'ohcsum18-22' = '10e22 Joules',
+                  'ohcsum23-46' = '10e22 Joules', 'vsftmyz' = 'Sv')
 
 syears <- seq(year0, yearf, intsdate)
 imon2 <- paste("0", as.character(mon0), sep = "")
@@ -100,18 +100,18 @@ if (file.exists(paste(savename, '.sav', sep = ''))) {
                        NULL)
   }
   latbnd <- switch(var, 'tos' = c(-60, 65), 'prlr' = c(-60, 65), c(-90, 90))
-  varname <- switch(var, 'sia' = paste(var, pole, sep = ''), 'sie' = paste(var, 
-             pole, sep = ''), 'siv' = paste(var, pole, sep = ''), 
-             'ohc' = 'heatc', 'uohc' = switch(vertnem, 'L42' = '0-315_heatc', 
+  varname <- switch(var, 'siarea' = paste(var, pole, sep = ''), 'siextent' = paste(var, 
+             pole, sep = ''), 'sivol' = paste(var, pole, sep = ''), 
+             'ohcsum0' = 'heatc', 'ohcsum23-46' = switch(vertnem, 'L42' = '0-315_heatc', 
                                                        'L46' = '0-322_heatc',
                                                        'L75' = '0-271_heatc'),
-                              'mohc' = switch(vertnem, 'L42' = '373-657_heatc',
+                              'ohcsum18-22' = switch(vertnem, 'L42' = '373-657_heatc',
                                                        'L46' = '382-735_heatc',
                                                        'L75' = '301-773_heatc'),
-                              'lohc' = switch(vertnem, 'L42' = '800-5350_heatc',
+                              'ohcsum0-17' = switch(vertnem, 'L42' = '800-5350_heatc',
                                                        'L46' = '855-5875_heatc',
                                                        'L75' = '857-5902_heatc'),
-              'amoc' = 'moc_40N55N_1-2km', var)
+              'vsftmyz' = 'moc_40N55N_1-2km', var)
   if (is.na(match('b02p', lstexpid)) == TRUE) { 
     lstload <- lstexpid
   } else {
@@ -146,13 +146,13 @@ if (file.exists(paste(savename, '.sav', sep = ''))) {
     toto1$mod <- toto1$mod * 1000 * 3600 * 24
     toto1$obs <- toto1$obs * 1000 * 3600 * 24
   }
-  if (var == 'ohc' | var == 'lohc' | var == 'mohc' | var == 'uohc') {
+  if (var == 'ohcsum0' | var == 'ohcsum0-17' | var == 'ohcsum18-22' | var == 'ohcsum23-46') {
     toto1$mod <- toto1$mod / 1e22
     toto1$obs <- toto1$obs / 1e22
   }
-  if (var == 'sia' | var=='sie' | var=='siv') {
+  if (var == 'siarea' | var=='siextent' | var=='sivol') {
     toto1$mod <- toto1$mod/1000
-    if (var == 'siv') {
+    if (var == 'sivol') {
       toto1$obs <- toto1$obs/1000
     }
   }