diff --git a/inst/doc/usecase/ex1_1_tranform.R b/inst/doc/usecase/ex1_1_tranform.R
index 5fbeb6fc4f81cfba7758a4d7b489513ae6ca1a78..c805b45d8678ddc20b334c3b7f222498b19a0953 100644
--- a/inst/doc/usecase/ex1_1_tranform.R
+++ b/inst/doc/usecase/ex1_1_tranform.R
@@ -1,17 +1,25 @@
 # ------------------------------------------------------------------
-# Do the interpolation within Start(). Compare with Load().
+# This script shows you how to do the interpolation in Start(). It also 
+# uses s2dverification::Load to compare the results, which are identical
+# (only tiny difference due to round-up).
+#
+# The parameters in Start() for interpolation include 'transform', 
+# 'transform_extra_cells', 'transform_params', and 'transform_vars'.
+# 'transform' is the interpolation function. startR provides 'CDORemapper',
+# which is the wrapper function of s2dverification::CDORemap.
+# 'transform_extra_cells' defines the extra grid points you want to use for 
+# interpolation. The default value is 2. 'transform_params' is a list which
+# defines the arguments used in 'cdo'. 'transform_vars' is a vector indicating 
+# the dimensions to be interpolated.
 # ------------------------------------------------------------------
 
 library(startR)
-library(s2dverification)
 
 obs_path <- '/esarchive/recon/ecmwf/era5/monthly_mean/$var$_f1h/$var$_$sdate$.nc'
-pobs <- paste0('/esarchive/recon/ecmwf/era5/monthly_mean/',
-               '$VAR_NAME$_f1h/$VAR_NAME$_$YEAR$$MONTH$.nc')
-var_name<-'sfcWind'
+var_name <- 'sfcWind'
 
-lons.min <- 0
-lons.max <- 10
+lons.min <- 10
+lons.max <- 20
 lats.min <- 0
 lats.max <- 10
 
@@ -23,8 +31,9 @@ obs <- Start(dat = obs_path,
              sdate = '201811',
              time = 'all',
              latitude = values(list(lats.min, lats.max)),
-             latitude_reorder = Sort(),
+             latitude_reorder = Sort(decreasing = T),
              longitude = values(list(lons.min, lons.max)),
+             longitude_reorder = CircularSort(0, 360),
              synonims = list(longitude = c('lon', 'longitude'),
                              latitude = c('lat', 'latitude')),
              transform = CDORemapper,
@@ -36,9 +45,6 @@ obs <- Start(dat = obs_path,
              return_vars = list(time = NULL,
                                 latitude = 'dat',
                                 longitude = 'dat'),
-             split_multiselected_dims = TRUE,
-             merge_across_dims = TRUE,
-             #num_procs = 4,
              retrieve = T)
 
 # -- Result --
@@ -46,20 +52,25 @@ dim(obs)
 #      dat       var     sdate      time  latitude longitude 
 #        1         1         1         1        11        11 
 range(obs)
-#[1] 1.090472 4.987585
-str(attr(obs, 'Variables')$common[[var_name]]$dim[[1]])
-# $ name         : chr "longitude"
-# $ vals         : num [1:11(1d)] 0 1 2 3 4 5 6 7 8 9 ...
-str(attr(obs, 'Variables')$common[[var_name]]$dim[[2]])
-# $ name         : chr "latitude"
-# $ vals         : num [1:11(1d)] 0 1 2 3 4 5 6 7 8 9 ...
-obs[1,1,1,1,,1]
-# [1] 4.987585 4.299388 4.007287 3.758497 3.673881 3.908989 2.224407 1.774090
-# [9] 1.800998 1.718733 1.816074
+#[1] 1.090472 2.803161
+attr(obs, 'Variables')$dat1$longitude
+# [1] 10 11 12 13 14 15 16 17 18 19 20
+attr(obs, 'Variables')$dat1$latitude
+# [1] 10  9  8  7  6  5  4  3  2  1  0
+obs[1,1,1,1,1:3,1:2]
+#         [,1]     [,2]
+#[1,] 2.345519 2.365326
+#[2,] 2.219121 1.958319
+#[3,] 1.840537 1.386617
+
 
 #-------------------------
 # s2dverification::Load()
 #-------------------------
+library(s2dverification)
+pobs <- paste0('/esarchive/recon/ecmwf/era5/monthly_mean/',
+               '$VAR_NAME$_f1h/$VAR_NAME$_$YEAR$$MONTH$.nc')
+
 obs1 <- Load(var = var_name,
              obs = list(list(path = pobs)),
              sdates = '20181101',
@@ -79,21 +90,23 @@ dim(obs1$obs)
 #dataset  member   sdate   ftime     lat     lon 
 #      1       1       1       1      11      11 
 range(obs1$obs)
-#[1] 1.0905 4.9876
-str(obs1)
-# $ lon            : num [1:11(1d)] 0 1 2 3 4 5 6 7 8 9 ...
-# $ lat            : num [1:11(1d)] 10 9 8 7 6 5 4 3 2 1 ...
-obs1$obs[1,1,1,1,,1]  #when lon = 0
-# [1] 1.8161 1.7187 1.8010 1.7741 2.2244 3.9090 3.6739 3.7585 4.0073 4.2994
-#[11] 4.9876
-
+#[1] 1.0905 2.8032
+as.vector(obs1$lon)
+# [1] 10 11 12 13 14 15 16 17 18 19 20
+as.vector(obs1$lat)
+# [1] 10  9  8  7  6  5  4  3  2  1  0
+obs1$obs[1,1,1,1,1:3,1:2]
+#       [,1]   [,2]
+#[1,] 2.3455 2.3653
+#[2,] 2.2191 1.9583
+#[3,] 1.8405 1.3866
 
 #-----------------------
 # Comparison
 #-----------------------
-diff <- drop(obs)[11:1, ] - drop(obs1$obs)  #lat order is different
+diff <- drop(obs) - drop(obs1$obs) 
 range(diff)
-#[1] -4.918404e-05  4.693756e-05
+#[1] -4.965553e-05  4.929314e-05
 
 
 
diff --git a/inst/doc/usecase/ex2_6_ext_param_func.R b/inst/doc/usecase/ex2_6_ext_param_func.R
index 04baeb23b968b36dd6b819d09da1d9e1133a97df..302b52a8697d7aaf394b2554fe89819941cb9732 100644
--- a/inst/doc/usecase/ex2_6_ext_param_func.R
+++ b/inst/doc/usecase/ex2_6_ext_param_func.R
@@ -86,6 +86,7 @@
                  threads_load = 2,
                  threads_compute = 4)
 
+# Notice that the function uses rnorm() inside. So the results will be different.
 # -----------------------------------------------------------
 #names(res)
 #[1] "strat"  "t_test"