Index: /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/rebuild_nemo
===================================================================
--- /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/rebuild_nemo	(revision 10130)
+++ /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/rebuild_nemo	(revision 10131)
@@ -10,18 +10,25 @@
    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 -s] 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 "            -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 chunsize along x " 
+   echo "            -y chunsize along y " 
+   echo "            -z chunsize along z " 
+   echo "            -t chunsize along t " 
+   echo "            -s chunsize "
    echo
    exit 1
 }
 
-while getopts c:n:t:d:r:lm opt
+while getopts l:p:s:m:n:r:d:x:y:z:t:s opt
 do
   case ${opt} in
@@ -32,9 +39,9 @@
          echo "Submitting job to compute node"
       ;;
-      t)
+      p)
          OMP_NUM_THREADS=${OPTARG}
       ;;
-      c)
-         NCHUNKSIZE=${OPTARG}
+      s)
+         NSLICESIZE=${OPTARG}
       ;;
       m) 
@@ -50,4 +57,19 @@
       r)
          MEMORY=${OPTARG}
+      ;;
+      x)
+         NXCHUNK=${OPTARG}
+      ;;
+      y)
+         NYCHUNK=${OPTARG}
+      ;;
+      z)
+         NZCHUNK=${OPTARG}
+      ;;
+      t)
+         NTCHUNK=${OPTARG}
+      ;;
+      s)
+         CHUNKSIZE=${OPTARG}
       ;;
   esac
@@ -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,4 +121,21 @@
    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
Index: /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/src/rebuild_nemo.F90
===================================================================
--- /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/src/rebuild_nemo.F90	(revision 10130)
+++ /NEMO/branches/UKMO/dev_rebuild_nemo_compression/REBUILD_NEMO/src/rebuild_nemo.F90	(revision 10131)
@@ -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
+   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
@@ -154,10 +169,8 @@
    LOGICAL :: l_maskout = .false.
 
-   NAMELIST/nam_rebuild/ filebase, ndomain, dims, nchunksize, l_maskout, deflate_level
+   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 +187,40 @@
  
 !--------------------------------------------------------------------------------
+!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,*) '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,*) '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
 
 !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 +250,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 +314,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 +363,36 @@
       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))
+           if( TRIM(indimnames(dimids(idim))) == 'x' )             &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_xchunk,1))
+           if( TRIM(indimnames(dimids(idim))) == 'y' )             &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_ychunk,1))
+! trick to find var in a list of suggestion (var0 and var1 : INDEX(|var0|var1|,|var|)
+           if( INDEX('|depth|z|nav_lev|','|'//TRIM(indimnames(dimids(idim)))//'|') > 0 ) &
+    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_zchunk,1))
+           if( INDEX('|time|time_counter|','|'//TRIM(indimnames(dimids(idim)))//'|') > 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: ',chunkalg, chunksizes
+        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 +420,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 +438,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 +453,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 +569,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 +598,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 +683,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 +720,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 +726,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 +982,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 +1001,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 +1020,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 +1039,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 +1058,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 +1149,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 +1172,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 +1196,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 +1219,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 +1242,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 +1295,5 @@
             WRITE(numerr,*) "*** NEMO rebuild failed ***"
             WRITE(numerr,*)
-            STOP 4
+            STOP 5
          ENDIF
       ENDIF
