Index: /utils/tools/REBUILD_NEMO/BATCH_TEMPLATES/param_XC40_METO
===================================================================
--- /utils/tools/REBUILD_NEMO/BATCH_TEMPLATES/param_XC40_METO	(revision 10338)
+++ /utils/tools/REBUILD_NEMO/BATCH_TEMPLATES/param_XC40_METO	(revision 10338)
@@ -0,0 +1,3 @@
+#!/bin/ksh
+BATCH_CMD="qsub"
+ARCH="XC40_METO"
Index: /utils/tools/REBUILD_NEMO/README.rst
===================================================================
--- /utils/tools/REBUILD_NEMO/README.rst	(revision 10338)
+++ /utils/tools/REBUILD_NEMO/README.rst	(revision 10338)
@@ -0,0 +1,118 @@
+============
+REBUILD_NEMO
+============
+REBUILD_NEMO is a tool to rebuild NEMO output files from multiple processors (mesh_mask, restart or XIOS output files) into one file.
+
+Description
+===========
+
+NEMO rebuild has the following features:
+ * dynamically works out what variables require rebuilding
+ * does not copy subdomain halo regions
+ * works for 1,2,3 and 4d arrays or types for all valid NetCDF types
+ * utilises OMP shared memory parallelisation where applicable
+ * time 'slicing' for lower memory use  (only for 4D vars with unlimited dimension)
+
+The code reads the filestem and number of subdomains from the namelist file nam_rebuild.
+| The 1st subdomain file is used to determine the dimensions and variables in all the input files. 
+It is also used to find which dimensions (and hence which variables) require rebuilding 
+as well as information about the global domain.
+| It then opens all the input files (unbuffered) and creates an array of netcdf identifiers
+before looping through all the variables and updating the rebuilt output file (either by direct 
+copying or looping over the number of domains and rebuilding as appropriate).
+| The code looks more complicated than it is because it has lots of case statements to deal with all 
+the various NetCDF data types and with various data dimensions (up to 4d).
+| Diagnostic output is written to numout (default 6 - stdout)
+and errors are written to numerr (default 0 - stderr).
+| If time slicing is specified the code will use less memory but take a little longer.
+It does this by breaking down the 4D input variables over their 4th dimension 
+(generally time) by way of a while loop.
+
+How to compile it
+=================
+REBUILD_NEMO is compiled in the same manner as all the other tools. 
+The compilation script is maketools and the option are very similar to makenemo. 
+
+Here an example of how to compile REBUILD_NEMO on the MetOffice XC40 HPC:
+
+.. code-block:: console
+                      
+        $ ./maketools -n REBUILD_NEMO -m XC40_METO
+
+Usage
+=====
+There is 2 manners to use REBUILD_NEMO:
+
+* **rebuild_nemo shell script**
+
+If the rebuild_nemo shell script it used, the namelist is filled automatically depending on the on-line arguments.
+
+.. code-block:: console
+
+ $ ./rebuild_nemo
+
+  NEMO Rebuild
+  ************
+
+  usage: rebuild_nemo [-l -p -s -m -n -r -d -x -y -z -t -c] filebase ndomain [rebuild dimensions]
+
+  flags:    -l arch            submit to compute node
+            -p num             use num threads
+            -s num             split 4D vars into time slice of size num
+            -m                 force masking of global arrays (zero if no mdi)
+            -n namelist        full path to namelist file to be created (otherwise default nam_rebuild+_process_id is used)
+            -r memory          Memory to request on compute node including units (Default = 10Gb)
+
+      key_netcdf4 only 
+            -d deflate_level     deflate level for output files
+            -x chunksize along x 
+            -y chunksize along y 
+            -z chunksize along z 
+            -t chunksize along t 
+            -c total size of the chunk cache 
+
+In case the option '-l arch' is used, a template for the batch script and parameter for job submission has to be provided in BATCH_TEMPLATES (param_arch and rebuild_nemo_batch_arch). 
+Exemple from the XC40_METO architecture can be found in the directory. 
+Some keywords (NTHREADS, MEMORY, INDIR, NAMELIST and NOPEN) from the template are replaced by the rebuild_nemo script.
+
+* **rebuild_nemo.exe + namelist**
+
+If rebuild_nemo.exe is used directly, a namelist has to be provided. Default name is nam_rebuild. A specific name can be provided as argument.
+The minimal namelist required is (here example to rebuild a mesh_mask split in 36 files):
+
+.. code-block:: console
+
+         &nam_rebuild
+         filebase='mesh_mask'
+         ndomain=36
+         /
+
+Some option can be added (the value mentioned here are the default value):
+
+.. code-block:: console
+
+         l_maskout=.false       !(-m option: useful if input file comes from a run using land suppression)
+         nslicesize=0           !(-s option: 0 means no splitting in time slice)
+         deflate_level=0        !(-d option)
+         nc4_xchunk=206         !(-x option)
+         nc4_ychunk=135         !(-y option)
+         nc4_zchunk=1           !(-z option)
+         nc4_tchunk=1           !(-t option)
+         fchunksize=32000000    !(-c option)
+
+Example
+=======
+
+Here is an example of usage of rebuild_nemo shelf script
+(rebuild mesh_mask files on the XC40_METO computer with deflation level set to 1 and chunksize to (1,1,100,100)):
+
+.. code-block:: console
+
+          $ rebuild_nemo -l XC40_METO -m -d 1 -x 100 -y 100 -z 1 -t 1 mesh_mask 36
+
+          output is mask using netcdf missing value (_Fillvalue attribute) or 0 if missing value not in the netcdf.
+
+          file mesh_mask,  num_domains 36, num_threads 1
+          Submitting job to compute node
+          8510610.xcs00
+         $
Index: /utils/tools/REBUILD_NEMO/nam_rebuild
===================================================================
--- /utils/tools/REBUILD_NEMO/nam_rebuild	(revision 10338)
+++ /utils/tools/REBUILD_NEMO/nam_rebuild	(revision 10338)
@@ -0,0 +1,12 @@
+&nam_rebuild
+filebase='mesh_mask'
+ndomain=36
+l_maskout=.false       !(-m option: useful if input file comes from a run using land suppression)
+nslicesize=0           !(-s option: 0 means no splitting in time slice)
+deflate_level=1        !(-d option)
+nc4_xchunk=206         !(-x option)
+nc4_ychunk=135         !(-y option)
+nc4_zchunk=1           !(-z option)
+nc4_tchunk=1           !(-t option)
+fchunksize=32000000    !(-c option)
+/
Index: /utils/tools/REBUILD_NEMO/rebuild_nemo
===================================================================
--- /utils/tools/REBUILD_NEMO/rebuild_nemo	(revision 10337)
+++ /utils/tools/REBUILD_NEMO/rebuild_nemo	(revision 10338)
@@ -10,35 +10,42 @@
    echo "  ************"
    echo
-   echo "  usage: ${0##*/} [-l -t -c -m -n] filebase ndomain [rebuild dimensions]"
+   echo "  usage: ${0##*/} [-l -p -s -m -n -r -d -x -y -z -t -c] filebase ndomain [rebuild dimensions]"
    echo
-   echo "  flags:    -l                 submit to compute node"
-   echo "            -t num             use num threads"
-   echo "            -c num             split 4D vars into time chuncks of size num"
+   echo "  flags:    -l arch            submit to compute node"
+   echo "            -p num             use num threads"
+   echo "            -s num             split 4D vars into time slice of size num"
    echo "            -m                 force masking of global arrays (zero if no mdi)"
-   echo "            -d deflate_level   deflate level for output files (key_netcdf4 only)"
    echo "            -n namelist        full path to namelist file to be created (otherwise default nam_rebuild+_process_id is used)"
    echo "            -r memory          Memory to request on compute node including units (Default = 10Gb)"
+   echo ""
+   echo "      key_netcdf4 only "
+   echo "            -d deflate_level     deflate level for output files"
+   echo "            -x chunksize along x " 
+   echo "            -y chunksize along y " 
+   echo "            -z chunksize along z " 
+   echo "            -t chunksize along t " 
+   echo "            -c total size of the chunk cache "
    echo
    exit 1
 }
 
-while getopts c:n:t:d:r:lm opt
+while getopts l:p:s:n:r:d:x:y:z:t:c:m opt
 do
   case ${opt} in
       l) 
          BATCH="yes"
-         BATCH_CMD="qsub" 
-         ARCH="XC40_METO"
-         echo "Submitting job to compute node"
+         ARCH=${OPTARG}
       ;;
-      t)
+      p)
          OMP_NUM_THREADS=${OPTARG}
       ;;
-      c)
-         NCHUNKSIZE=${OPTARG}
+      s)
+         NSLICESIZE=${OPTARG}
       ;;
       m) 
          NMASK="TRUE"
-         echo "mask"
+         echo ""
+         echo "output is mask using netcdf missing value (_Fillvalue attribute) or 0 if missing value not in the netcdf."
+         echo ""
       ;;
       d)
@@ -50,4 +57,19 @@
       r)
          MEMORY=${OPTARG}
+      ;;
+      x)
+         NXCHUNK=${OPTARG}
+      ;;
+      y)
+         NYCHUNK=${OPTARG}
+      ;;
+      z)
+         NZCHUNK=${OPTARG}
+      ;;
+      t)
+         NTCHUNK=${OPTARG}
+      ;;
+      c)
+         CHUNKSIZE=${OPTARG}
       ;;
   esac
@@ -72,5 +94,5 @@
 nopen=$(ulimit -n)
 if [[ $ndomain -gt $nopen ]] ; then
-  nopen=$((ndomain+2))
+  nopen=$((ndomain+4))   # +2 failed !!!
 fi
 
@@ -90,6 +112,6 @@
    echo ${dims} >> $nam_rebuild
 fi
-if [[ -n ${NCHUNKSIZE} ]] ; then
-   echo " nchunksize=${NCHUNKSIZE}" >> $nam_rebuild
+if [[ -n ${NCSLICESIZE} ]] ; then
+   echo " nslicesize=${NCSLICESIZE}" >> $nam_rebuild
 fi
 if [[ -n ${NMASK} ]] ; then
@@ -99,13 +121,51 @@
    echo " deflate_level=${DEFLATE}" >> $nam_rebuild
 fi
+if [[ -n ${NXCHUNK} ]] ; then
+   echo " nc4_xchunk=${NXCHUNK}" >> $nam_rebuild
+fi
+if [[ -n ${NYCHUNK} ]] ; then
+   echo " nc4_ychunk=${NYCHUNK}" >> $nam_rebuild
+fi
+if [[ -n ${NZCHUNK} ]] ; then
+   echo " nc4_zchunk=${NZCHUNK}" >> $nam_rebuild
+fi
+if [[ -n ${NTCHUNK} ]] ; then
+   echo " nc4_tchunk=${NTCHUNK}" >> $nam_rebuild
+fi
+if [[ -n ${CHUNKSIZE} ]] ; then
+   echo " fchunksize=${CHUNKSIZE}" >> $nam_rebuild
+fi
+
+
 
 echo "/" >> $nam_rebuild
 
 if [[ ${BATCH} == "yes" ]] ; then
+
+   template_dir=${script_dir}/BATCH_TEMPLATES/
+   param_file=${template_dir}/param_${ARCH}
+   if [ ! -f $param_file ]; then 
+      echo ''
+      echo "E R R O R: $param_file is missing, stop 42"
+      echo ''
+      echo "check your arch name or add one $param_file file in BATCH_TEMPLATES"
+      echo ''
+      exit 42
+   fi
+   . $param_file
+
    batch_file=rebuild_nemo_batch_${ARCH}
+   if [ ! -f ${template_dir}/${batch_file} ]; then 
+      echo ''
+      echo "E R R O R: $batch_file is missing, stop 42"
+      echo ''
+      echo "check your arch name or add one $batch_file file in BATCH_TEMPLATES"
+      echo ''
+      exit 42
+   fi
 
    #Create a modified local copy of the batch submission file
    #The process ID is appended to the end of the file name so it is unique
-   cat ${script_dir}/BATCH_TEMPLATES/${batch_file} | sed -e"s/NTHREADS/${OMP_NUM_THREADS}/" \
+   cat ${template_dir}/${batch_file} | sed -e"s/NTHREADS/${OMP_NUM_THREADS}/" \
       -e"s/MEMORY/${MEMORY}/" \
       -e"s:INDIR:${script_dir}:" \
@@ -114,4 +174,5 @@
    
    #Submit the job
+   echo "Submitting job to compute node"
    $BATCH_CMD ${batch_file}_$$.sh
 
Index: /utils/tools/REBUILD_NEMO/src/rebuild_nemo.F90
===================================================================
--- /utils/tools/REBUILD_NEMO/src/rebuild_nemo.F90	(revision 10337)
+++ /utils/tools/REBUILD_NEMO/src/rebuild_nemo.F90	(revision 10338)
@@ -1,4 +1,4 @@
 PROGRAM rebuild_nemo
-
+#define key_netcdf4
    !!=========================================================================
    !!                        ***  rebuild_nemo  ***
@@ -14,9 +14,10 @@
    !!     * works for 1,2,3 and 4d arrays or types for all valid NetCDF types
    !!     * utilises OMP shared memory parallelisation where applicable
-   !!     * time 'chunking' for lower memory use 
+   !!     * time 'slicing' for lower memory use 
    !!       (only for 4D vars with unlimited dimension)
    !!
    !!  Ed Blockley - August 2011
    !!  (based on original code by Matt Martin)
+   !!  Julien Palmieri and Andrew Coward - September 2018 (add compression and chunking)
    !!
    !!-------------------------------------------------------------------------
@@ -38,5 +39,5 @@
    !!  and errors are written to numerr (default 0 - stderr).
    !!
-   !!  If time chunking is specified the code will use less memory but take a little longer.
+   !!  If time slicing is specified the code will use less memory but take a little longer.
    !!  It does this by breaking down the 4D input variables over their 4th dimension 
    !!  (generally time) by way of a while loop.
@@ -66,18 +67,32 @@
    CHARACTER(LEN=nf90_max_name), ALLOCATABLE :: filenames(:), indimnames(:)
    CHARACTER(LEN=nf90_max_name), DIMENSION(2) :: dims
-
-   INTEGER :: ndomain, ifile, ndomain_file, nchunksize, deflate_level
+   CHARACTER(LEN=256) :: cnampath, cdimlst, cdim
+   CHARACTER(LEN=50)  :: clibnc ! netcdf library version
+
+   INTEGER :: ndomain, ifile, ndomain_file, nslicesize, deflate_level
    INTEGER :: ncid, outid, idim, istop
    INTEGER :: natts, attid, xtype, varid, rbdims 
    INTEGER :: jv, ndims, nvars, dimlen, dimids(4)
    INTEGER :: dimid, unlimitedDimId, di, dj, dr
-   INTEGER :: nmax_unlimited, nt, ntchunk 
-   INTEGER :: chunksize = 3200000
+   INTEGER :: nmax_unlimited, nt, ntslice 
+   INTEGER :: fchunksize = 32000000  ! NetCDF global file chunk cache size
+   INTEGER :: patchchunk             ! NetCDF processor-domain file chunk cache size
    INTEGER :: nthreads = 1
+   INTEGER :: chunkalg = 0          ! NetCDF4 variable chunking algorithm
+                                    ! Default variable chunksizes (typical ORCA025
+                                    ! recommendations which can be adjusted via namelist
+                                    ! or will be bounded if too large for domain.)
+   INTEGER :: nc4_xchunk = 206      ! Default x (longitude) variable chunk size
+   INTEGER :: nc4_ychunk = 135      ! Default y (latitude)  variable chunk size
+   INTEGER :: nc4_zchunk = 1        ! Default z (depth) variable chunk size (almost always 1)
+   INTEGER :: nc4_tchunk = 1        ! Default t (time)  variable chunk size (almost always 1)
    INTEGER, ALLOCATABLE  :: outdimids(:), outdimlens(:), indimlens(:), inncids(:)
+   INTEGER, ALLOCATABLE  :: chunksizes(:)
    INTEGER, ALLOCATABLE  :: global_sizes(:), rebuild_dims(:)
    INTEGER, DIMENSION(2) :: halo_start, halo_end, local_sizes
    INTEGER, DIMENSION(2) :: idomain, jdomain, rdomain, start_pos
    INTEGER :: ji, jj, jk, jl, jr
+   INTEGER :: nargs                 ! number of arguments
+   INTEGER, EXTERNAL :: iargc
  
    REAL(sp) :: ValMin, ValMax, InMin, InMax, rmdi
@@ -149,15 +164,14 @@
    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_dp
 
-   LOGICAL :: l_valid = .false.
+   LOGICAL :: l_valid     = .false.
    LOGICAL :: l_noRebuild = .false.
-   LOGICAL :: l_findDims = .true.
-   LOGICAL :: l_maskout = .false.
-
-   NAMELIST/nam_rebuild/ filebase, ndomain, dims, nchunksize, l_maskout, deflate_level
+   LOGICAL :: l_findDims  = .true.
+   LOGICAL :: l_maskout   = .false.
+   LOGICAL :: l_namexist  = .false.
+
+   NAMELIST/nam_rebuild/ filebase, ndomain, dims, nslicesize, l_maskout, deflate_level, &
+                       & nc4_xchunk, nc4_ychunk, nc4_zchunk, nc4_tchunk, fchunksize         
 
    external      :: getarg
-   integer, external :: iargc
-   integer   :: found_num_args = 0
-   character(256) :: namelist_path
 
    !End of definitions 
@@ -174,30 +188,49 @@
  
 !--------------------------------------------------------------------------------
+!1.0 Check netcdf version for warning
+   clibnc = TRIM(nf90_inq_libvers())
+   IF (ICHAR(clibnc(1:1)) <= 3) THEN
+      PRINT *, '=========================================================='
+      PRINT *, 'You are using old netcdf library (',TRIM(clibnc),').'
+      PRINT *, 'REBUILD_NEMO support of old netcdf library will end soon'
+      PRINT *, 'please consider moving to netcdf 4 or higher'
+      PRINT *, '=========================================================='
+   END IF
+
 !1.1 Get the namelist path
-!Determine the number of arguments on the command line
-found_num_args=iargc()
-!Check that the required argument is present, if it is not then set it to the default value: nam_rebuild
-IF (found_num_args <= 0) THEN
-   WRITE(numout,*) 'Namelist path not supplied as command line argument. Using default, nam_rebuild.'
-   namelist_path='nam_rebuild'
-ELSE IF (found_num_args == 1) THEN
-   CALL getarg(1, namelist_path)
-ELSE 
-   WRITE(numerr,*) 'ERROR! : Incorrect number of command line arguments. Please supply only'
-   WRITE(numerr,*) '         the path to the namelist file, or no arguments to use default value'
-   STOP 1
-END IF
+   !Determine the number of arguments on the command line
+   nargs=iargc()
+   !Check that the required argument is present, if it is not then set it to the default value: nam_rebuild
+   IF (nargs == 0) THEN
+      WRITE(numout,*)
+      WRITE(numout,*) 'W A R N I N G : Namelist path not supplied as command line argument. Using default, nam_rebuild.'
+      cnampath='nam_rebuild'
+   ELSE IF (nargs == 1) THEN
+      CALL getarg(1, cnampath)
+   ELSE 
+      WRITE(numerr,*) 'E R R O R ! : Incorrect number of command line arguments. Please supply only'
+      WRITE(numerr,*) '         the path to the namelist file, or no arguments to use default value'
+      STOP 1
+   END IF
+
+   ! check presence of namelist
+   INQUIRE(FILE=cnampath, EXIST=l_namexist)
+   IF (.NOT. l_namexist) THEN
+      WRITE(numout,*)
+      WRITE(numout,*) 'E R R O R : Namelist '//TRIM(cnampath)//' not present.'
+      STOP 42
+   END IF
 
 !1.2 Read in the namelist 
 
    dims(:) = ""
-   nchunksize = 0
+   nslicesize = 0
    deflate_level = 0
-   OPEN( UNIT=numnam, FILE=TRIM(namelist_path), FORM='FORMATTED', STATUS='OLD' )
+   OPEN( UNIT=numnam, FILE=TRIM(cnampath), FORM='FORMATTED', STATUS='OLD' )
    READ( numnam, nam_rebuild )
    CLOSE( numnam )
    IF( .NOT. ALL(dims(:) == "") ) l_findDims = .false.
- 
-!1.1 Set up the filenames and fileids
+
+!1.3 Set up the filenames and fileids
 
    ALLOCATE(filenames(ndomain))
@@ -227,9 +260,9 @@
   
 !2.1 Set up the output file
-!#if defined key_netcdf4
-   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_netcdf4, outid, chunksize=chunksize ) )
-!#else
-!   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_64bit_offset, outid, chunksize=chunksize ) )
-!#endif
+#if defined key_netcdf4
+   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_netcdf4, outid, chunksize=fchunksize ) )
+#else
+   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_64bit_offset, outid, chunksize=fchunksize ) )
+#endif
 
 !2.2 Set up dimensions in output file
@@ -291,5 +324,5 @@
       outdimlens(idim) = dimlen
    END DO
-   ! nmax_unlimited is only used for time-chunking so we set it to be at least 1 to 
+   ! nmax_unlimited is only used for time-slicing so we set it to be at least 1 to 
    ! account for files with no record dimension or zero length record dimension(!)
    nmax_unlimited = max(nmax_unlimited,1)
@@ -340,14 +373,48 @@
       CALL check_nf90( nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts ) )
       ALLOCATE(outdimids(ndims))
-      DO idim = 1, ndims
-         outdimids(idim) = dimids(idim)
-      END DO
+      ALLOCATE(chunksizes(ndims))
+      IF( ndims > 0 ) then
+        DO idim = 1, ndims
+           outdimids(idim) = dimids(idim)
+           chunksizes(idim) = outdimlens(dimids(idim))
+           cdim='|'//TRIM(indimnames(dimids(idim)))//'|'
+
+! trick to find var in a list of suggestion (var0 and var1 : INDEX(|var0|var1|,|var|)
+           cdimlst='|x|x_grid_T|x_grid_U|x_grid_V|x_grid_W|'
+           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_xchunk,1))
+
+           cdimlst='|y|y_grid_T|y_grid_U|y_grid_V|y_grid_W|'
+           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_ychunk,1))
+
+           cdimlst='|z|deptht|depthu|depthv|depthw|depth|nav_lev|'
+           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_zchunk,1))
+
+           cdimlst='|t|time|time_counter|'
+           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_tchunk,1))
+
+        END DO
 #if defined key_netcdf4
-      CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid, &
-                                     deflate_level=deflate_level ) )
+        CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid, &
+                                       deflate_level=deflate_level ) )
+        IF (l_verbose) WRITE(numout,*) 'Dims    : ',ndims, outdimids(1:ndims)
+        IF (l_verbose) WRITE(numout,*) 'names   : ',(TRIM(indimnames(dimids(idim)))//' ',idim=1,ndims)
+        IF (l_verbose) WRITE(numout,*) 'lens    : ',(outdimlens(dimids(idim)),idim=1,ndims)
+        IF (l_verbose) WRITE(numout,*) 'Chunking: ',chunksizes
+        IF (l_verbose) WRITE(numout,*) 'Deflation : ',deflate_level
+        IF (l_verbose) WRITE(numout,*) 'Chunk algo: ',chunkalg
+        CALL check_nf90( nf90_def_var_chunking( outid, varid, chunkalg, &
+   &                                 chunksizes ) )
+      ELSE
+        CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid ) )
 #else
       CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid ) )
 #endif
+      ENDIF
       DEALLOCATE(outdimids)
+      DEALLOCATE(chunksizes)
       IF (l_verbose) WRITE(numout,*) 'Defining variable '//TRIM(varname)//'...' 
       IF( natts > 0 ) THEN
@@ -375,6 +442,11 @@
 
    IF (l_verbose) WRITE(numout,*) 'Opening input files...'
+
+   ! Set a file chunk cache size for the processor-domain files that scales with the number of processors
+   patchchunk = max(8192, fchunksize/ndomain)
+
+   ! open files 
    DO ifile = 2, ndomain
-      CALL check_nf90( nf90_open( TRIM(filenames(ifile)), nf90_share, ncid, chunksize=chunksize ) )
+      CALL check_nf90( nf90_open( TRIM(filenames(ifile)), nf90_share, ncid, chunksize=patchchunk ) )
       inncids(ifile) = ncid
    END DO
@@ -388,6 +460,6 @@
       istop = nf90_noerr
       nt = 1
-      ntchunk = nmax_unlimited
-      IF( nchunksize == 0 ) nchunksize = nmax_unlimited
+      ntslice = nmax_unlimited
+      IF( nslicesize == 0 ) nslicesize = nmax_unlimited
 
 !3.2 Inquire variable to find out name and how many dimensions it has
@@ -403,19 +475,19 @@
       ENDIF
 
-!3.2.0 start while loop for time chunking
+!3.2.0 start while loop for time slicing
 
       DO WHILE( nt <= nmax_unlimited )
 
          IF( ndims > 3 ) THEN
-            ntchunk = MIN( nchunksize, nmax_unlimited + 1 - nt )
+            ntslice = MIN( nslicesize, nmax_unlimited + 1 - nt )
          ENDIF
 
       IF (l_noRebuild) THEN
 
-         IF( nchunksize == nmax_unlimited .OR. ndims <= 3 ) THEN
+         IF( nslicesize == nmax_unlimited .OR. ndims <= 3 ) THEN
             IF (l_verbose) WRITE(numout,*) 'Copying data from variable '//TRIM(varname)//'...'
          ELSE
             IF (l_verbose) WRITE(numout,'(A,I3,A,I3,A)') ' Copying data from variable '  &
-            &                 //TRIM(varname)//' for chunks ',nt,' to ',nt+ntchunk-1,' ...'
+            &                 //TRIM(varname)//' for slices ',nt,' to ',nt+ntslice-1,' ...'
          ENDIF
 
@@ -519,21 +591,21 @@
                CASE( NF90_BYTE )
                   ALLOCATE(globaldata_4d_i1(indimlens(dimids(1)),indimlens(dimids(2)),       &
-                     &                      indimlens(dimids(3)),ntchunk))
+                     &                      indimlens(dimids(3)),ntslice))
                   CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i1, start=(/1,1,1,nt/) ) )
                CASE( NF90_SHORT )
                   ALLOCATE(globaldata_4d_i2(indimlens(dimids(1)),indimlens(dimids(2)),       &
-                     &                      indimlens(dimids(3)),ntchunk))
+                     &                      indimlens(dimids(3)),ntslice))
                   CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i2, start=(/1,1,1,nt/) ) )
                CASE( NF90_INT )
                   ALLOCATE(globaldata_4d_i4(indimlens(dimids(1)),indimlens(dimids(2)),       &
-                     &                      indimlens(dimids(3)),ntchunk))
+                     &                      indimlens(dimids(3)),ntslice))
                   CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i4, start=(/1,1,1,nt/) ) )
                CASE( NF90_FLOAT )
                   ALLOCATE(globaldata_4d_sp(indimlens(dimids(1)),indimlens(dimids(2)),       &
-                     &                      indimlens(dimids(3)),ntchunk))
+                     &                      indimlens(dimids(3)),ntslice))
                   CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_sp, start=(/1,1,1,nt/) ) )
                CASE( NF90_DOUBLE )
                   ALLOCATE(globaldata_4d_dp(indimlens(dimids(1)),indimlens(dimids(2)),       &
-                     &                      indimlens(dimids(3)),ntchunk))
+                     &                      indimlens(dimids(3)),ntslice))
                   CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_dp, start=(/1,1,1,nt/) ) )
                CASE DEFAULT
@@ -548,9 +620,9 @@
 !3.2.2 For variables that require rebuilding we need to read in from all ndomain files
 !      Here we allocate global variables ahead of looping over files
-         IF( nchunksize == nmax_unlimited .OR. ndims <= 3 ) THEN
+         IF( nslicesize == nmax_unlimited .OR. ndims <= 3 ) THEN
             IF (l_verbose) WRITE(numout,*) 'Rebuilding data from variable '//TRIM(varname)//'...'
          ELSE
             IF (l_verbose) WRITE(numout,'(A,I3,A,I3,A)') ' Rebuilding data from variable '  &
-            &                 //TRIM(varname)//' for chunks ',nt,' to ',nt+ntchunk-1,' ...'
+            &                 //TRIM(varname)//' for slices ',nt,' to ',nt+ntslice-1,' ...'
          ENDIF
          IF( ndims == 1 ) THEN
@@ -633,21 +705,21 @@
                CASE( NF90_BYTE )
                   ALLOCATE(globaldata_4d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
-                     &                      outdimlens(dimids(3)),ntchunk))
+                     &                      outdimlens(dimids(3)),ntslice))
                   IF (l_maskout) globaldata_4d_i1(:,:,:,:)=mdiVals(jv)
                CASE( NF90_SHORT )
                   ALLOCATE(globaldata_4d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
-                     &                      outdimlens(dimids(3)),ntchunk))
+                     &                      outdimlens(dimids(3)),ntslice))
                   IF (l_maskout) globaldata_4d_i2(:,:,:,:)=mdiVals(jv)
                CASE( NF90_INT )
                   ALLOCATE(globaldata_4d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
-                     &                      outdimlens(dimids(3)),ntchunk))
+                     &                      outdimlens(dimids(3)),ntslice))
                   IF (l_maskout) globaldata_4d_i4(:,:,:,:)=mdiVals(jv)
                CASE( NF90_FLOAT )
                   ALLOCATE(globaldata_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
-                     &                      outdimlens(dimids(3)),ntchunk))
+                     &                      outdimlens(dimids(3)),ntslice))
                   IF (l_maskout) globaldata_4d_sp(:,:,:,:)=mdiVals(jv)
                CASE( NF90_DOUBLE )
                   ALLOCATE(globaldata_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
-                     &                      outdimlens(dimids(3)),ntchunk))
+                     &                      outdimlens(dimids(3)),ntslice))
                   IF (l_maskout) globaldata_4d_dp(:,:,:,:)=mdiVals(jv)
                CASE DEFAULT
@@ -670,5 +742,5 @@
 !$OMP&         localdata_1d_i1,localdata_2d_i1,localdata_3d_i1,localdata_4d_i1)           &
 !$OMP& SHARED(jv,nvars,varname,filenames,ValMin,ValMax,indimlens,outdimlens,rbdims,       &
-!$OMP&        ndomain,outid,chunksize,istop,l_valid,nthreads,inncids,rebuild_dims,        &
+!$OMP&        ndomain,outid,fchunksize,istop,l_valid,nthreads,inncids,rebuild_dims,       &
 !$OMP&        globaldata_1d_i2,globaldata_1d_i4,globaldata_1d_sp,globaldata_1d_dp,        &
 !$OMP&        globaldata_2d_i2,globaldata_2d_i4,globaldata_2d_sp,globaldata_2d_dp,        &
@@ -676,5 +748,5 @@
 !$OMP&        globaldata_4d_i2,globaldata_4d_i4,globaldata_4d_sp,globaldata_4d_dp,        &
 !$OMP&        globaldata_1d_i1,globaldata_2d_i1,globaldata_3d_i1,globaldata_4d_i1,        &
-!$OMP&        ntchunk,nt,nmax_unlimited,indimnames,dims)
+!$OMP&        ntslice,nt,nmax_unlimited,indimnames,dims,patchchunk)
 
          DO ifile = 1, ndomain
@@ -932,9 +1004,9 @@
                   CASE( NF90_BYTE )
                      ALLOCATE(localdata_4d_i1(local_sizes(di),local_sizes(dj),               &
-                         &                     indimlens(dimids(3)),ntchunk))
+                         &                     indimlens(dimids(3)),ntslice))
                      CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i1, start=(/1,1,1,nt/) ), istop )
 !$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
-!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i1,localdata_4d_i1,di,dj,nt,ntchunk)
-                     DO jl = 1, ntchunk
+!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i1,localdata_4d_i1,di,dj,nt,ntslice)
+                     DO jl = 1, ntslice
 !$OMP DO 
                         DO jk = 1, indimlens(dimids(3))
@@ -951,9 +1023,9 @@
                   CASE( NF90_SHORT )
                      ALLOCATE(localdata_4d_i2(local_sizes(di),local_sizes(dj),               &
-                        &                     indimlens(dimids(3)),ntchunk))
+                        &                     indimlens(dimids(3)),ntslice))
                      CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i2, start=(/1,1,1,nt/) ), istop )
 !$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
-!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i2,localdata_4d_i2,di,dj,nt,ntchunk)
-                     DO jl = 1, ntchunk
+!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i2,localdata_4d_i2,di,dj,nt,ntslice)
+                     DO jl = 1, ntslice
 !$OMP DO 
                         DO jk = 1, indimlens(dimids(3))
@@ -970,9 +1042,9 @@
                   CASE( NF90_INT )
                      ALLOCATE(localdata_4d_i4(local_sizes(di),local_sizes(dj),               &
-                        &                     indimlens(dimids(3)),ntchunk))
+                        &                     indimlens(dimids(3)),ntslice))
                      CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i4, start=(/1,1,1,nt/) ), istop )
 !$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
-!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i4,localdata_4d_i4,di,dj,nt,ntchunk)
-                     DO jl = 1, ntchunk
+!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i4,localdata_4d_i4,di,dj,nt,ntslice)
+                     DO jl = 1, ntslice
 !$OMP DO
                         DO jk = 1, indimlens(dimids(3))
@@ -989,9 +1061,9 @@
                   CASE( NF90_FLOAT )
                      ALLOCATE(localdata_4d_sp(local_sizes(di),local_sizes(dj),               &
-                        &                     indimlens(dimids(3)),ntchunk))
+                        &                     indimlens(dimids(3)),ntslice))
                      CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_sp, start=(/1,1,1,nt/) ), istop )
 !$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
-!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_sp,localdata_4d_sp,di,dj,nt,ntchunk) 
-                     DO jl = 1, ntchunk
+!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_sp,localdata_4d_sp,di,dj,nt,ntslice) 
+                     DO jl = 1, ntslice
 !$OMP DO
                         DO jk = 1, indimlens(dimids(3))
@@ -1008,9 +1080,9 @@
                   CASE( NF90_DOUBLE )
                      ALLOCATE(localdata_4d_dp(local_sizes(di),local_sizes(dj),               &
-                        &                     indimlens(dimids(3)),ntchunk))
+                        &                     indimlens(dimids(3)),ntslice))
                      CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_dp, start=(/1,1,1,nt/) ), istop )
 !$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
-!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_dp,localdata_4d_dp,di,dj,nt,ntchunk) 
-                     DO jl = 1, ntchunk
+!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_dp,localdata_4d_dp,di,dj,nt,ntslice) 
+                     DO jl = 1, ntslice
 !$OMP DO
                         DO jk = 1, indimlens(dimids(3))
@@ -1099,4 +1171,7 @@
             CASE( NF90_DOUBLE )
                CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_dp ) )
+            CASE DEFAULT   
+               WRITE(numerr,*) '0d Unknown nf90 type: ', xtype
+               STOP 4
          END SELECT
 
@@ -1119,4 +1194,7 @@
                CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_dp ) )
                DEALLOCATE(globaldata_1d_dp)
+            CASE DEFAULT   
+               WRITE(numerr,*) '1d Unknown nf90 type: ', xtype
+               STOP 4
          END SELECT
 
@@ -1140,5 +1218,5 @@
                DEALLOCATE(globaldata_2d_dp)
             CASE DEFAULT   
-               WRITE(numerr,*) 'Unknown nf90 type: ', xtype
+               WRITE(numerr,*) '2d Unknown nf90 type: ', xtype
                STOP 4
          END SELECT     
@@ -1163,5 +1241,5 @@
                DEALLOCATE(globaldata_3d_dp)
             CASE DEFAULT   
-               WRITE(numerr,*) 'Unknown nf90 type: ', xtype
+               WRITE(numerr,*) '3d Unknown nf90 type: ', xtype
                STOP 4
          END SELECT     
@@ -1186,11 +1264,13 @@
                DEALLOCATE(globaldata_4d_dp)
             CASE DEFAULT   
-               WRITE(numerr,*) 'Unknown nf90 type: ', xtype
+               WRITE(numerr,*) '4d Unknown nf90 type: ', xtype
                STOP 4
          END SELECT     
+         ! why only for big data set, test the cost.
+         CALL check_nf90( nf90_sync( outid ) )    ! flush buffers to disk after writing big 4D datasets
     
       ENDIF
 
-         nt = nt + ntchunk
+         nt = nt + ntslice
 
       END DO ! WHILE loop
@@ -1237,5 +1317,5 @@
             WRITE(numerr,*) "*** NEMO rebuild failed ***"
             WRITE(numerr,*)
-            STOP 4
+            STOP 5
          ENDIF
       ENDIF
