Developers/Rebuild Zooms: rebuild_nemo_modified.F90

File rebuild_nemo_modified.F90, 64.1 KB (added by acc, 5 months ago)
Line 
1PROGRAM rebuild_nemo
2#define key_netcdf4
3   !!=========================================================================
4   !!                        ***  rebuild_nemo  ***
5   !!=========================================================================
6   !!
7   !!  A routine to rebuild NEMO files from multiple processors into one file.
8   !!  This routine is designed to be much quicker than the old IOIPSL rebuild
9   !!  but at the cost of an increased memory usage.
10   !!
11   !!  NEMO rebuild has the following features:
12   !!     * dynamically works out what variables require rebuilding
13   !!     * does not copy subdomain halo regions
14   !!     * works for 1,2,3 and 4d arrays or types for all valid NetCDF types
15   !!     * utilises OMP shared memory parallelisation where applicable
16   !!     * time 'slicing' for lower memory use
17   !!       (only for 4D vars with unlimited dimension)
18   !!
19   !!  Ed Blockley - August 2011
20   !!  (based on original code by Matt Martin)
21   !!  Julien Palmieri and Andrew Coward - September 2018 (add compression and chunking)
22   !!
23   !!-------------------------------------------------------------------------
24   !!
25   !!  The code reads the filestem and number of subdomains from the namelist file nam_rebuild.
26   !!
27   !!  The 1st subdomain file is used to determine the dimensions and variables in all the input files.
28   !!  It is also used to find which dimensions (and hence which variables) require rebuilding
29   !!  as well as information about the global domain.
30   !!
31   !!  It then opens all the input files (unbuffered) and creates an array of netcdf identifiers
32   !!  before looping through all the variables and updating the rebuilt output file (either by direct
33   !!  copying or looping over the number of domains and rebuilding as appropriate).
34   !! 
35   !!  The code looks more complicated than it is because it has lots of case statements to deal with all
36   !!  the various NetCDF data types and with various data dimensions (up to 4d).
37   !!
38   !!  Diagnostic output is written to numout (default 6 - stdout)
39   !!  and errors are written to numerr (default 0 - stderr).
40   !!
41   !!  If time slicing is specified the code will use less memory but take a little longer.
42   !!  It does this by breaking down the 4D input variables over their 4th dimension
43   !!  (generally time) by way of a while loop.
44   !!
45   !!-------------------------------------------------------------------------------
46 
47   USE netcdf
48
49!$ USE omp_lib           ! Note OpenMP sentinel
50
51   IMPLICIT NONE
52
53   ! kind specifications
54   INTEGER,PARAMETER :: i1=SELECTED_INT_KIND(2)          ! NF90_BYTE
55   INTEGER,PARAMETER :: i2=SELECTED_INT_KIND(4)          ! NF90_SHORT
56   INTEGER,PARAMETER :: i4=SELECTED_INT_KIND(9)          ! NF90_INT
57   INTEGER,PARAMETER :: sp=SELECTED_REAL_KIND(6,37)      ! NF90_FLOAT
58   INTEGER,PARAMETER :: dp=SELECTED_REAL_KIND(12,307)    ! NF90_DOUBLE
59
60   INTEGER,PARAMETER :: numnam = 11
61   INTEGER,PARAMETER :: numout = 6
62   INTEGER,PARAMETER :: numerr = 0
63
64   LOGICAL, PARAMETER :: l_verbose = .true. 
65   
66   CHARACTER(LEN=nf90_max_name) :: filebase, suffix, attname, dimname, varname, time, date, zone, timestamp
67   CHARACTER(LEN=nf90_max_name), ALLOCATABLE :: filenames(:), indimnames(:)
68   CHARACTER(LEN=nf90_max_name), DIMENSION(2) :: dims
69   CHARACTER(LEN=256) :: cnampath, cdimlst, cdim
70   CHARACTER(LEN=50)  :: clibnc ! netcdf library version
71
72   INTEGER :: ndomain, ifile, ndomain_file, nslicesize, deflate_level
73   INTEGER :: ncid, outid, idim, istop, istat
74   INTEGER :: natts, attid, xtype, varid, rbdims, rbdims1
75   INTEGER :: jv, ndims, nvars, dimlen, dimids(4)
76   INTEGER :: dimid, unlimitedDimId, di, dj, dr
77   INTEGER :: nmax_unlimited, nt, ntslice 
78   INTEGER :: fchunksize = 32000000  ! NetCDF global file chunk cache size
79   INTEGER :: patchchunk             ! NetCDF processor-domain file chunk cache size
80   INTEGER :: nthreads = 1
81   INTEGER :: chunkalg = 0          ! NetCDF4 variable chunking algorithm
82                                    ! Default variable chunksizes (typical ORCA025
83                                    ! recommendations which can be adjusted via namelist
84                                    ! or will be bounded if too large for domain.)
85   INTEGER :: nc4_xchunk = 206      ! Default x (longitude) variable chunk size
86   INTEGER :: nc4_ychunk = 135      ! Default y (latitude)  variable chunk size
87   INTEGER :: nc4_zchunk = 1        ! Default z (depth) variable chunk size (almost always 1)
88   INTEGER :: nc4_tchunk = 1        ! Default t (time)  variable chunk size (almost always 1)
89   INTEGER, ALLOCATABLE  :: outdimids(:), outdimlens(:), indimlens(:), inncids(:)
90   INTEGER, ALLOCATABLE  :: chunksizes(:)
91   INTEGER, ALLOCATABLE  :: global_sizes(:), rebuild_dims(:)
92   INTEGER, DIMENSION(2) :: halo_start, halo_end, local_sizes
93   INTEGER, DIMENSION(2) :: idomain, jdomain, rdomain, start_pos
94   INTEGER :: ji, jj, jk, jl, jr
95   INTEGER :: ni_zoom, nj_zoom, ni_off, nj_off
96   LOGICAL :: iszoom
97   INTEGER :: nargs                 ! number of arguments
98   INTEGER, EXTERNAL :: iargc
99 
100   REAL(sp) :: ValMin, ValMax, InMin, InMax, rmdi
101   REAL(dp), ALLOCATABLE :: mdiVals(:)
102
103   ! NF90_BYTE local data arrays
104   INTEGER(i1), ALLOCATABLE, SAVE, DIMENSION(:) :: localdata_1d_i1
105   INTEGER(i1), ALLOCATABLE, SAVE, DIMENSION(:,:) :: localdata_2d_i1
106   INTEGER(i1), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: localdata_3d_i1
107   INTEGER(i1), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: localdata_4d_i1
108
109   ! NF90_SHORT local data arrays
110   INTEGER(i2), ALLOCATABLE, SAVE, DIMENSION(:) :: localdata_1d_i2
111   INTEGER(i2), ALLOCATABLE, SAVE, DIMENSION(:,:) :: localdata_2d_i2
112   INTEGER(i2), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: localdata_3d_i2
113   INTEGER(i2), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: localdata_4d_i2
114
115   ! NF90_INT local data arrays
116   INTEGER(i4), ALLOCATABLE, SAVE, DIMENSION(:) :: localdata_1d_i4
117   INTEGER(i4), ALLOCATABLE, SAVE, DIMENSION(:,:) :: localdata_2d_i4
118   INTEGER(i4), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: localdata_3d_i4
119   INTEGER(i4), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: localdata_4d_i4
120
121   ! NF90_FLOAT local data arrays
122   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:) :: localdata_1d_sp
123   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: localdata_2d_sp
124   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: localdata_3d_sp
125   REAL(sp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: localdata_4d_sp
126
127   ! NF90_DOUBLE local data arrays
128   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:) :: localdata_1d_dp
129   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: localdata_2d_dp
130   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: localdata_3d_dp
131   REAL(dp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: localdata_4d_dp
132
133   ! NF90_BYTE global data arrays
134   INTEGER(i1) :: globaldata_0d_i1
135   INTEGER(i1), ALLOCATABLE, DIMENSION(:) :: globaldata_1d_i1
136   INTEGER(i1), ALLOCATABLE, DIMENSION(:,:) :: globaldata_2d_i1
137   INTEGER(i1), ALLOCATABLE, DIMENSION(:,:,:) :: globaldata_3d_i1
138   INTEGER(i1), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_i1
139
140   ! NF90_SHORT global data arrays
141   INTEGER(i2) :: globaldata_0d_i2
142   INTEGER(i2), ALLOCATABLE, DIMENSION(:) :: globaldata_1d_i2
143   INTEGER(i2), ALLOCATABLE, DIMENSION(:,:) :: globaldata_2d_i2
144   INTEGER(i2), ALLOCATABLE, DIMENSION(:,:,:) :: globaldata_3d_i2
145   INTEGER(i2), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_i2
146
147   ! NF90_INT global data arrays
148   INTEGER(i4) :: globaldata_0d_i4
149   INTEGER(i4), ALLOCATABLE, DIMENSION(:) :: globaldata_1d_i4
150   INTEGER(i4), ALLOCATABLE, DIMENSION(:,:) :: globaldata_2d_i4
151   INTEGER(i4), ALLOCATABLE, DIMENSION(:,:,:) :: globaldata_3d_i4
152   INTEGER(i4), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_i4
153 
154   ! NF90_FLOAT global data arrays
155   REAL(sp) :: globaldata_0d_sp
156   REAL(sp), ALLOCATABLE, DIMENSION(:) :: globaldata_1d_sp
157   REAL(sp), ALLOCATABLE, DIMENSION(:,:) :: globaldata_2d_sp
158   REAL(sp), ALLOCATABLE, DIMENSION(:,:,:) :: globaldata_3d_sp
159   REAL(sp), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_sp
160
161   ! NF90_DOUBLE global data arrays
162   REAL(dp) :: globaldata_0d_dp
163   REAL(dp), ALLOCATABLE, DIMENSION(:) :: globaldata_1d_dp
164   REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: globaldata_2d_dp
165   REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: globaldata_3d_dp
166   REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: globaldata_4d_dp
167
168   LOGICAL :: l_valid     = .false.
169   LOGICAL :: l_noRebuild = .false.
170   LOGICAL :: l_findDims  = .true.
171   LOGICAL :: l_maskout   = .false.
172   LOGICAL :: l_namexist  = .false.
173
174   NAMELIST/nam_rebuild/ filebase, ndomain, dims, nslicesize, l_maskout, deflate_level, &
175                       & nc4_xchunk, nc4_ychunk, nc4_zchunk, nc4_tchunk, fchunksize         
176
177   external      :: getarg
178
179   !End of definitions
180
181!--------------------------------------------------------------------------------
182!0. OMP setup
183
184!$OMP PARALLEL DEFAULT(NONE) SHARED(nthreads)
185!$OMP MASTER
186!$      nthreads = omp_get_num_threads()
187!$      WRITE(numout,*) 'Running OMP with ',nthreads,' thread(s).'
188!$OMP END MASTER
189!$OMP END PARALLEL
190 
191!--------------------------------------------------------------------------------
192!1.0 Check netcdf version for warning
193   clibnc = TRIM(nf90_inq_libvers())
194   IF (ICHAR(clibnc(1:1)) <= 3) THEN
195      PRINT *, '=========================================================='
196      PRINT *, 'You are using old netcdf library (',TRIM(clibnc),').'
197      PRINT *, 'REBUILD_NEMO support of old netcdf library will end soon'
198      PRINT *, 'please consider moving to netcdf 4 or higher'
199      PRINT *, '=========================================================='
200   END IF
201
202!1.1 Get the namelist path
203   !Determine the number of arguments on the command line
204   nargs=iargc()
205   !Check that the required argument is present, if it is not then set it to the default value: nam_rebuild
206   IF (nargs == 0) THEN
207      WRITE(numout,*)
208      WRITE(numout,*) 'W A R N I N G : Namelist path not supplied as command line argument. Using default, nam_rebuild.'
209      cnampath='nam_rebuild'
210   ELSE IF (nargs == 1) THEN
211      CALL getarg(1, cnampath)
212   ELSE
213      WRITE(numerr,*) 'E R R O R ! : Incorrect number of command line arguments. Please supply only'
214      WRITE(numerr,*) '         the path to the namelist file, or no arguments to use default value'
215      STOP 1
216   END IF
217
218   ! check presence of namelist
219   INQUIRE(FILE=cnampath, EXIST=l_namexist)
220   IF (.NOT. l_namexist) THEN
221      WRITE(numout,*)
222      WRITE(numout,*) 'E R R O R : Namelist '//TRIM(cnampath)//' not present.'
223      STOP 42
224   END IF
225
226!1.2 Read in the namelist
227
228   dims(:) = ""
229   nslicesize = 0
230   deflate_level = 0
231   OPEN( UNIT=numnam, FILE=TRIM(cnampath), FORM='FORMATTED', STATUS='OLD' )
232   READ( numnam, nam_rebuild )
233   CLOSE( numnam )
234   IF( .NOT. ALL(dims(:) == "") ) l_findDims = .false.
235
236!1.3 Set up the filenames and fileids
237
238   ALLOCATE(filenames(ndomain))
239   IF (l_verbose) WRITE(numout,*) 'Rebuilding the following files:'
240   DO ifile = 1, ndomain
241      WRITE(suffix,'(i4.4)') ifile-1
242      filenames(ifile) = TRIM(filebase)//'_'//TRIM(suffix)//'.nc'
243      IF (l_verbose) WRITE(numout,*) TRIM(filenames(ifile))
244   END DO
245   ALLOCATE(inncids(ndomain))
246 
247!---------------------------------------------------------------------------
248!2. Read in the global dimensions from the first input file and set up the output file
249 
250   CALL check_nf90( nf90_open( TRIM(filenames(1)), nf90_share, ncid ) )
251   CALL check_nf90( nf90_inquire( ncid, ndims, nvars, natts ) )
252   
253!2.0 Read in the total number of processors the file is expecting and check it's correct 
254 
255   CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_number_total', ndomain_file ) )
256   IF( ndomain /= ndomain_file ) THEN
257      WRITE(numerr,*) 'WARNING! : number of files to rebuild in file does not agree with namelist'
258      WRITE(numerr,*) 'Attribute DOMAIN_number_total is : ', ndomain_file
259      WRITE(numerr,*) 'Number of files specified in namelist is: ', ndomain
260      istat = nf90_inquire_attribute( ncid, nf90_global, 'DOMAIN_size_zoom_i', xtype, rbdims, attid )
261      IF ( istat == nf90_noerr ) THEN
262          WRITE(numerr,*) 'This looks like a zoom region so I will assume you know what you are doing'
263      ELSE
264          WRITE(numerr,*) 'This is a potentially fatal error'
265          STOP 2
266      ENDIF
267   ENDIF
268 
269!2.1 Set up the output file
270#if defined key_netcdf4
271   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_netcdf4, outid, chunksize=fchunksize ) )
272#else
273   CALL check_nf90( nf90_create( TRIM(filebase)//'.nc', nf90_64bit_offset, outid, chunksize=fchunksize ) )
274#endif
275
276!2.2 Set up dimensions in output file
277
278!2.2.0 Find out how many dimensions are required to be rebuilt and which ones they are
279   CALL check_nf90( nf90_inquire_attribute( ncid, nf90_global, 'DOMAIN_dimensions_ids', xtype, rbdims, attid ) )
280
281   ALLOCATE(rebuild_dims(rbdims))
282   CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_dimensions_ids', rebuild_dims ) )
283
284   ALLOCATE(global_sizes(rbdims))
285   CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_size_global', global_sizes ) )
286
287!2.2.0.1 Override global sizes if zoom attributes are found
288   iszoom = .false.
289   istat = nf90_inquire_attribute( ncid, nf90_global, 'DOMAIN_size_zoom_i', xtype, rbdims1, attid )
290   IF ( istat == nf90_noerr ) THEN
291       CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_size_zoom_i', ni_zoom ) )
292       ! Need both zoom_i and zoom_j attributes to determine zoom size
293       istat = nf90_inquire_attribute( ncid, nf90_global, 'DOMAIN_size_zoom_j', xtype, rbdims1, attid )
294       IF ( istat == nf90_noerr ) THEN
295           CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_size_zoom_j', nj_zoom ) )
296           iszoom = .true.
297           global_sizes(1) = ni_zoom
298           global_sizes(2) = nj_zoom
299           CALL check_nf90( nf90_get_att( ncid, nf90_global, 'ibegin', ni_off ) )
300           CALL check_nf90( nf90_get_att( ncid, nf90_global, 'jbegin', nj_off ) )
301       ELSE
302           iszoom = .false.
303       ENDIF
304   ENDIF
305     
306   IF (l_verbose) WRITE(numout,*) 'Size of global arrays: ', global_sizes
307
308
309!2.2.1 Copy the dimensions into the output file apart from rebuild_dims() which are dimensioned globally
310   ALLOCATE(indimlens(ndims), indimnames(ndims), outdimlens(ndims))
311   CALL check_nf90( nf90_inquire( ncid, unlimitedDimId = unlimitedDimId ) )
312   istop = 0
313   DO idim = 1, ndims
314      CALL check_nf90( nf90_inquire_dimension( ncid, idim, dimname, dimlen ) )
315      CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_size_local', local_sizes ) )
316      indimlens(idim) = dimlen   
317      indimnames(idim) = dimname
318      IF (l_findDims) THEN
319         IF( idim == rebuild_dims(1) ) THEN
320            IF( dimlen == local_sizes(1) ) THEN
321               dimlen = global_sizes(1)
322               dims(1) = trim(dimname)
323            ELSE
324               istop = 1
325            ENDIF
326         ENDIF
327         IF( rbdims > 1 .AND. idim == rebuild_dims(2) ) THEN
328            IF( dimlen == local_sizes(2) ) THEN
329               dimlen = global_sizes(2)
330               dims(2) = trim(dimname)
331            ELSE
332               istop = 1
333            ENDIF
334         ENDIF
335      ELSE ! l_findDims = false
336         IF( TRIM(dimname) == TRIM(dims(1))) THEN
337            dimlen = global_sizes(1)
338            rebuild_dims(1) = idim
339         ENDIF
340         IF( rbdims > 1 .AND. TRIM(dimname) == TRIM(dims(2))) THEN
341            dimlen = global_sizes(2)
342            rebuild_dims(2) = idim
343         ENDIF
344      ENDIF
345
346      IF( idim == unlimitedDimId ) THEN
347         CALL check_nf90( nf90_def_dim( outid, dimname, nf90_unlimited, dimid) )
348         nmax_unlimited = dimlen
349      ELSE
350         CALL check_nf90( nf90_def_dim( outid, dimname, dimlen, dimid) )
351      ENDIF
352      outdimlens(idim) = dimlen
353   END DO
354   ! nmax_unlimited is only used for time-slicing so we set it to be at least 1 to
355   ! account for files with no record dimension or zero length record dimension(!)
356   nmax_unlimited = max(nmax_unlimited,1)
357
358   IF( istop == 1 ) THEN
359      WRITE(numerr,*) 'ERROR! : DOMAIN_local_sizes attribute does not match rebuild dimension lengths in the first file'
360      WRITE(numerr,*) 'Attribute DOMAIN_local_sizes is : ', local_sizes
361      WRITE(numerr,*) 'Dimensions to be rebuilt are of size : ', outdimlens(rebuild_dims(1)), outdimlens(rebuild_dims(2)) 
362      STOP 3
363   ENDIF
364
365   IF (l_findDims) THEN
366      IF (l_verbose) WRITE(numout,*) 'Finding rebuild dimensions from the first file...'
367   ELSE
368      IF (l_verbose) WRITE(numout,*) 'Using rebuild dimensions given in namelist...'
369   ENDIF
370
371   IF( rbdims > 1 ) THEN
372      IF (l_verbose) WRITE(numout,*) 'Rebuilding across dimensions '//TRIM(indimnames(rebuild_dims(1)))//  &
373         &                      ' and '//TRIM(indimnames(rebuild_dims(2)))
374   ELSE
375      IF (l_verbose) WRITE(numout,*) 'Rebuilding across dimension '//TRIM(indimnames(rebuild_dims(1)))
376   ENDIF
377     
378!2.2.2 Copy the global attributes into the output file, apart from those beginning with DOMAIN_ 
379!      Also need to change the file_name attribute and the TimeStamp attribute.
380   DO attid = 1, natts
381      CALL check_nf90( nf90_inq_attname( ncid, nf90_global, attid, attname ) )
382      IF( INDEX( attname, "DOMAIN_" ) == 1 ) CYCLE
383      IF( INDEX( attname, "file_name") == 1 ) CYCLE
384      IF( INDEX( attname, "associate_file") == 1 ) CYCLE
385      IF (l_verbose) WRITE(numout,*) 'Copying attribute '//TRIM(attname)//' into destination file...'
386      CALL check_nf90( nf90_copy_att( ncid, nf90_global, attname, outid, nf90_global ) )
387   END DO
388   CALL check_nf90( nf90_put_att( outid, nf90_global, "file_name", TRIM(filebase)//'.nc') )
389   IF (l_verbose) WRITE(numout,*) 'Writing new file_name attribute' 
390   CALL DATE_AND_TIME ( date=date, time=time, zone=zone )
391   timestamp = date(7:8) // "/" // date(5:6) // "/" // date(1:4) // " " // &
392               time(1:2) // ":" // time(3:4) // ":" // time(5:6) // " " // &
393               zone 
394   CALL check_nf90( nf90_put_att( outid, nf90_global, "TimeStamp", timestamp ) )
395   IF (l_verbose) WRITE(numout,*) 'Writing new TimeStamp attribute'
396 
397!2.2.3 Copy the variable definitions and attributes into the output file.
398   ALLOCATE(mdiVals(nvars))
399   mdiVals(:)=0
400   DO jv = 1, nvars
401      CALL check_nf90( nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts ) )
402      ALLOCATE(outdimids(ndims))
403      ALLOCATE(chunksizes(ndims))
404      IF( ndims > 0 ) then
405        DO idim = 1, ndims
406           outdimids(idim) = dimids(idim)
407           chunksizes(idim) = outdimlens(dimids(idim))
408           cdim='|'//TRIM(indimnames(dimids(idim)))//'|'
409
410! trick to find var in a list of suggestion (var0 and var1 : INDEX(|var0|var1|,|var|)
411           cdimlst='|x|x_grid_T|x_grid_U|x_grid_V|x_grid_W|'
412           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
413    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_xchunk,1))
414
415           cdimlst='|y|y_grid_T|y_grid_U|y_grid_V|y_grid_W|'
416           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
417    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_ychunk,1))
418
419           cdimlst='|z|deptht|depthu|depthv|depthw|depth|nav_lev|'
420           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
421    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_zchunk,1))
422
423           cdimlst='|t|time|time_counter|'
424           if( INDEX(TRIM(cdimlst),TRIM(cdim)) > 0 ) &
425    &                             chunksizes(idim) = min(outdimlens(dimids(idim)), max(nc4_tchunk,1))
426
427        END DO
428#if defined key_netcdf4
429        CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid, &
430                                       deflate_level=deflate_level ) )
431        IF (l_verbose) WRITE(numout,*) 'Dims    : ',ndims, outdimids(1:ndims)
432        IF (l_verbose) WRITE(numout,*) 'names   : ',(TRIM(indimnames(dimids(idim)))//' ',idim=1,ndims)
433        IF (l_verbose) WRITE(numout,*) 'lens    : ',(outdimlens(dimids(idim)),idim=1,ndims)
434        IF (l_verbose) WRITE(numout,*) 'Chunking: ',chunksizes
435        IF (l_verbose) WRITE(numout,*) 'Deflation : ',deflate_level
436        IF (l_verbose) WRITE(numout,*) 'Chunk algo: ',chunkalg
437        CALL check_nf90( nf90_def_var_chunking( outid, varid, chunkalg, &
438   &                                 chunksizes ) )
439      ELSE
440        CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid ) )
441#else
442      CALL check_nf90( nf90_def_var( outid, varname, xtype, outdimids, varid ) )
443#endif
444      ENDIF
445      DEALLOCATE(outdimids)
446      DEALLOCATE(chunksizes)
447      IF (l_verbose) WRITE(numout,*) 'Defining variable '//TRIM(varname)//'...' 
448      IF( natts > 0 ) THEN
449         DO attid = 1, natts
450            CALL check_nf90( nf90_inq_attname( ncid, varid, attid, attname ) )
451            IF ( attname == "_FillValue" ) THEN
452               CALL check_nf90( nf90_get_att( ncid, varid, attname, rmdi ) )
453               mdiVals(jv)=rmdi
454            ENDIF
455            CALL check_nf90( nf90_copy_att( ncid, varid, attname, outid, varid ) )
456         END DO
457      ENDIF
458   END DO
459
460!2.3 End definitions in output file and copy 1st file ncid to the inncids array
461
462   CALL check_nf90( nf90_enddef( outid ) )
463   inncids(1) = ncid
464   IF (l_verbose) WRITE(numout,*) 'Finished defining output file.'
465 
466!---------------------------------------------------------------------------
467!3. Read in data from each file for each variable
468
469!3.1 Open each file and store the ncid in inncids array
470
471   IF (l_verbose) WRITE(numout,*) 'Opening input files...'
472
473   ! Set a file chunk cache size for the processor-domain files that scales with the number of processors
474   patchchunk = max(8192, fchunksize/ndomain)
475
476   ! open files
477   DO ifile = 2, ndomain
478      CALL check_nf90( nf90_open( TRIM(filenames(ifile)), nf90_share, ncid, chunksize=patchchunk ) )
479      inncids(ifile) = ncid
480   END DO
481   IF (l_verbose) WRITE(numout,*) 'All input files open.'
482
483   DO jv = 1, nvars
484
485      ValMin = 1.e10
486      ValMax = -1.e10
487      l_valid = .false.
488      istop = nf90_noerr
489      nt = 1
490      ntslice = nmax_unlimited
491      IF( nslicesize == 0 ) nslicesize = nmax_unlimited
492
493!3.2 Inquire variable to find out name and how many dimensions it has
494!    and importantly whether it contains the dimensions in rebuild_dims()
495
496      ncid = inncids(1)
497      CALL check_nf90( nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts ) )
498
499      l_noRebuild = .true.
500      IF( ANY( dimids(1:ndims) == rebuild_dims(1) )) l_noRebuild = .false.
501      IF( rbdims > 1 ) THEN
502         IF( ANY( dimids(1:ndims) == rebuild_dims(2) )) l_noRebuild = .false.
503      ENDIF
504
505!3.2.0 start while loop for time slicing
506
507      DO WHILE( nt <= nmax_unlimited )
508
509         IF( ndims > 3 ) THEN
510            ntslice = MIN( nslicesize, nmax_unlimited + 1 - nt )
511         ENDIF
512
513      IF (l_noRebuild) THEN
514
515         IF( nslicesize == nmax_unlimited .OR. ndims <= 3 ) THEN
516            IF (l_verbose) WRITE(numout,*) 'Copying data from variable '//TRIM(varname)//'...'
517         ELSE
518            IF (l_verbose) WRITE(numout,'(A,I3,A,I3,A)') ' Copying data from variable '  &
519            &                 //TRIM(varname)//' for slices ',nt,' to ',nt+ntslice-1,' ...'
520         ENDIF
521
522!3.2.1 If rebuilding not required then just need to read in variable
523!      for copying direct into output file after the OMP (files) loop.
524         IF( ndims == 0 ) THEN
525
526            SELECT CASE( xtype )
527               CASE( NF90_BYTE )
528                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_0d_i1 ) )
529               CASE( NF90_SHORT )
530                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_0d_i2 ) )
531               CASE( NF90_INT )
532                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_0d_i4 ) )
533               CASE( NF90_FLOAT )
534                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_0d_sp ) )
535               CASE( NF90_DOUBLE )
536                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_0d_dp ) )
537               CASE DEFAULT
538                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
539                  STOP 4
540            END SELECT
541
542         ELSEIF( ndims == 1 ) THEN
543
544            SELECT CASE( xtype )
545               CASE( NF90_BYTE )
546                  ALLOCATE(globaldata_1d_i1(indimlens(dimids(1))))
547                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_1d_i1 ) )
548               CASE( NF90_SHORT )
549                  ALLOCATE(globaldata_1d_i2(indimlens(dimids(1))))
550                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_1d_i2 ) )
551               CASE( NF90_INT )
552                  ALLOCATE(globaldata_1d_i4(indimlens(dimids(1))))
553                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_1d_i4 ) )
554               CASE( NF90_FLOAT )
555                  ALLOCATE(globaldata_1d_sp(indimlens(dimids(1))))
556                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_1d_sp ) )
557               CASE( NF90_DOUBLE )
558                  ALLOCATE(globaldata_1d_dp(indimlens(dimids(1))))
559                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_1d_dp ) )
560               CASE DEFAULT
561                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
562                  STOP 4
563            END SELECT
564
565         ELSEIF( ndims == 2 ) THEN
566
567            SELECT CASE( xtype )
568               CASE( NF90_BYTE )
569                  ALLOCATE(globaldata_2d_i1(indimlens(dimids(1)),indimlens(dimids(2))))
570                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_2d_i1 ) )
571               CASE( NF90_SHORT )
572                  ALLOCATE(globaldata_2d_i2(indimlens(dimids(1)),indimlens(dimids(2))))
573                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_2d_i2 ) )
574               CASE( NF90_INT )
575                  ALLOCATE(globaldata_2d_i4(indimlens(dimids(1)),indimlens(dimids(2))))
576                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_2d_i4 ) )
577               CASE( NF90_FLOAT )
578                  ALLOCATE(globaldata_2d_sp(indimlens(dimids(1)),indimlens(dimids(2))))
579                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_2d_sp ) )
580               CASE( NF90_DOUBLE )
581                  ALLOCATE(globaldata_2d_dp(indimlens(dimids(1)),indimlens(dimids(2))))
582                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_2d_dp ) )
583               CASE DEFAULT
584                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
585                  STOP 4
586            END SELECT
587
588         ELSEIF( ndims == 3 ) THEN
589
590            SELECT CASE( xtype )
591               CASE( NF90_BYTE )
592                  ALLOCATE(globaldata_3d_i1(indimlens(dimids(1)),indimlens(dimids(2)),       &
593                     &                      indimlens(dimids(3))))
594                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_3d_i1 ) )
595               CASE( NF90_SHORT )
596                  ALLOCATE(globaldata_3d_i2(indimlens(dimids(1)),indimlens(dimids(2)),       &
597                     &                      indimlens(dimids(3))))
598                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_3d_i2 ) )
599               CASE( NF90_INT )
600                  ALLOCATE(globaldata_3d_i4(indimlens(dimids(1)),indimlens(dimids(2)),       &
601                     &                      indimlens(dimids(3))))
602                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_3d_i4 ) )
603               CASE( NF90_FLOAT )
604                  ALLOCATE(globaldata_3d_sp(indimlens(dimids(1)),indimlens(dimids(2)),       &
605                     &                      indimlens(dimids(3))))
606                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_3d_sp ) )
607               CASE( NF90_DOUBLE )
608                  ALLOCATE(globaldata_3d_dp(indimlens(dimids(1)),indimlens(dimids(2)),       &
609                     &                      indimlens(dimids(3))))
610                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_3d_dp ) )
611               CASE DEFAULT
612                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
613                  STOP 4
614            END SELECT
615
616         ELSEIF( ndims == 4 ) THEN
617
618            SELECT CASE( xtype )
619               CASE( NF90_BYTE )
620                  ALLOCATE(globaldata_4d_i1(indimlens(dimids(1)),indimlens(dimids(2)),       &
621                     &                      indimlens(dimids(3)),ntslice))
622                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i1, start=(/1,1,1,nt/) ) )
623               CASE( NF90_SHORT )
624                  ALLOCATE(globaldata_4d_i2(indimlens(dimids(1)),indimlens(dimids(2)),       &
625                     &                      indimlens(dimids(3)),ntslice))
626                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i2, start=(/1,1,1,nt/) ) )
627               CASE( NF90_INT )
628                  ALLOCATE(globaldata_4d_i4(indimlens(dimids(1)),indimlens(dimids(2)),       &
629                     &                      indimlens(dimids(3)),ntslice))
630                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_i4, start=(/1,1,1,nt/) ) )
631               CASE( NF90_FLOAT )
632                  ALLOCATE(globaldata_4d_sp(indimlens(dimids(1)),indimlens(dimids(2)),       &
633                     &                      indimlens(dimids(3)),ntslice))
634                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_sp, start=(/1,1,1,nt/) ) )
635               CASE( NF90_DOUBLE )
636                  ALLOCATE(globaldata_4d_dp(indimlens(dimids(1)),indimlens(dimids(2)),       &
637                     &                      indimlens(dimids(3)),ntslice))
638                  CALL check_nf90( nf90_get_var( ncid, jv, globaldata_4d_dp, start=(/1,1,1,nt/) ) )
639               CASE DEFAULT
640                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
641                  STOP 4
642            END SELECT
643
644         ENDIF
645
646      ELSE  ! l_noRebuild = .false.
647
648!3.2.2 For variables that require rebuilding we need to read in from all ndomain files
649!      Here we allocate global variables ahead of looping over files
650         IF( nslicesize == nmax_unlimited .OR. ndims <= 3 ) THEN
651            IF (l_verbose) WRITE(numout,*) 'Rebuilding data from variable '//TRIM(varname)//'...'
652         ELSE
653            IF (l_verbose) WRITE(numout,'(A,I3,A,I3,A)') ' Rebuilding data from variable '  &
654            &                 //TRIM(varname)//' for slices ',nt,' to ',nt+ntslice-1,' ...'
655         ENDIF
656         IF( ndims == 1 ) THEN
657
658            SELECT CASE( xtype )
659               CASE( NF90_BYTE )
660                  ALLOCATE(globaldata_1d_i1(outdimlens(dimids(1))))
661                  IF (l_maskout) globaldata_1d_i1(:)=mdiVals(jv)
662               CASE( NF90_SHORT )
663                  ALLOCATE(globaldata_1d_i2(outdimlens(dimids(1))))
664                  IF (l_maskout) globaldata_1d_i2(:)=mdiVals(jv)
665               CASE( NF90_INT )
666                  ALLOCATE(globaldata_1d_i4(outdimlens(dimids(1))))
667                  IF (l_maskout) globaldata_1d_i4(:)=mdiVals(jv)
668               CASE( NF90_FLOAT )
669                  ALLOCATE(globaldata_1d_sp(outdimlens(dimids(1))))
670                  IF (l_maskout) globaldata_1d_sp(:)=mdiVals(jv)
671               CASE( NF90_DOUBLE )
672                  ALLOCATE(globaldata_1d_dp(outdimlens(dimids(1))))
673                  IF (l_maskout) globaldata_1d_dp(:)=mdiVals(jv)
674               CASE DEFAULT
675                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
676                  STOP 4
677            END SELECT
678
679         ELSEIF( ndims == 2 ) THEN
680
681            SELECT CASE( xtype )
682               CASE( NF90_BYTE )
683                  ALLOCATE(globaldata_2d_i1(outdimlens(dimids(1)),outdimlens(dimids(2))))
684                  IF (l_maskout) globaldata_2d_i1(:,:)=mdiVals(jv)
685               CASE( NF90_SHORT )
686                  ALLOCATE(globaldata_2d_i2(outdimlens(dimids(1)),outdimlens(dimids(2))))
687                  IF (l_maskout) globaldata_2d_i2(:,:)=mdiVals(jv)
688               CASE( NF90_INT )
689                  ALLOCATE(globaldata_2d_i4(outdimlens(dimids(1)),outdimlens(dimids(2))))
690                  IF (l_maskout) globaldata_2d_i4(:,:)=mdiVals(jv)
691               CASE( NF90_FLOAT )
692                  ALLOCATE(globaldata_2d_sp(outdimlens(dimids(1)),outdimlens(dimids(2))))
693                  IF (l_maskout) globaldata_2d_sp(:,:)=mdiVals(jv)
694               CASE( NF90_DOUBLE )
695                  ALLOCATE(globaldata_2d_dp(outdimlens(dimids(1)),outdimlens(dimids(2))))
696                  IF (l_maskout) globaldata_2d_dp(:,:)=mdiVals(jv)
697               CASE DEFAULT
698                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
699                  STOP 4
700            END SELECT
701
702         ELSEIF( ndims == 3 ) THEN
703
704            SELECT CASE( xtype )
705               CASE( NF90_BYTE )
706                  ALLOCATE(globaldata_3d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
707                     &                      outdimlens(dimids(3))))
708                  IF (l_maskout) globaldata_3d_i1(:,:,:)=mdiVals(jv)
709               CASE( NF90_SHORT )
710                  ALLOCATE(globaldata_3d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
711                     &                      outdimlens(dimids(3))))
712                  IF (l_maskout) globaldata_3d_i2(:,:,:)=mdiVals(jv)
713               CASE( NF90_INT )
714                  ALLOCATE(globaldata_3d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
715                     &                      outdimlens(dimids(3))))
716                  IF (l_maskout) globaldata_3d_i4(:,:,:)=mdiVals(jv)
717               CASE( NF90_FLOAT )
718                  ALLOCATE(globaldata_3d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
719                     &                      outdimlens(dimids(3))))
720                  IF (l_maskout) globaldata_3d_sp(:,:,:)=mdiVals(jv)
721               CASE( NF90_DOUBLE )
722                  ALLOCATE(globaldata_3d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
723                     &                      outdimlens(dimids(3))))
724                  IF (l_maskout) globaldata_3d_dp(:,:,:)=mdiVals(jv)
725               CASE DEFAULT
726                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
727                  STOP 4
728            END SELECT
729
730         ELSEIF( ndims == 4 ) THEN
731
732            SELECT CASE( xtype )
733               CASE( NF90_BYTE )
734                  ALLOCATE(globaldata_4d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
735                     &                      outdimlens(dimids(3)),ntslice))
736                  IF (l_maskout) globaldata_4d_i1(:,:,:,:)=mdiVals(jv)
737               CASE( NF90_SHORT )
738                  ALLOCATE(globaldata_4d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
739                     &                      outdimlens(dimids(3)),ntslice))
740                  IF (l_maskout) globaldata_4d_i2(:,:,:,:)=mdiVals(jv)
741               CASE( NF90_INT )
742                  ALLOCATE(globaldata_4d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
743                     &                      outdimlens(dimids(3)),ntslice))
744                  IF (l_maskout) globaldata_4d_i4(:,:,:,:)=mdiVals(jv)
745               CASE( NF90_FLOAT )
746                  ALLOCATE(globaldata_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
747                     &                      outdimlens(dimids(3)),ntslice))
748                  IF (l_maskout) globaldata_4d_sp(:,:,:,:)=mdiVals(jv)
749               CASE( NF90_DOUBLE )
750                  ALLOCATE(globaldata_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)),     &
751                     &                      outdimlens(dimids(3)),ntslice))
752                  IF (l_maskout) globaldata_4d_dp(:,:,:,:)=mdiVals(jv)
753               CASE DEFAULT
754                  WRITE(numerr,*) 'Unknown nf90 type: ', xtype
755                  STOP 4
756            END SELECT
757         ELSE
758            WRITE(numerr,*) 'ERROR! : A netcdf variable has more than 4 dimensions which is not taken into account'
759            STOP 4
760         ENDIF
761
762!$OMP  PARALLEL DO DEFAULT(NONE)                                                          &
763!$OMP& PRIVATE(ifile,ncid,xtype,start_pos,local_sizes,InMin,InMax,natts,                  &
764!$OMP&         ndims,attid,attname,dimids,idim,dimname,dimlen,unlimitedDimId,             &
765!$OMP&         halo_start,halo_end,idomain,jdomain,rdomain,di,dj,dr,                      &
766!$OMP&         localdata_1d_i2,localdata_1d_i4,localdata_1d_sp,localdata_1d_dp,           &
767!$OMP&         localdata_2d_i2,localdata_2d_i4,localdata_2d_sp,localdata_2d_dp,           &
768!$OMP&         localdata_3d_i2,localdata_3d_i4,localdata_3d_sp,localdata_3d_dp,           &
769!$OMP&         localdata_4d_i2,localdata_4d_i4,localdata_4d_sp,localdata_4d_dp,           &
770!$OMP&         localdata_1d_i1,localdata_2d_i1,localdata_3d_i1,localdata_4d_i1)           &
771!$OMP& SHARED(jv,nvars,varname,filenames,ValMin,ValMax,indimlens,outdimlens,rbdims,       &
772!$OMP&        ndomain,outid,fchunksize,istop,l_valid,nthreads,inncids,rebuild_dims,       &
773!$OMP&        globaldata_1d_i2,globaldata_1d_i4,globaldata_1d_sp,globaldata_1d_dp,        &
774!$OMP&        globaldata_2d_i2,globaldata_2d_i4,globaldata_2d_sp,globaldata_2d_dp,        &
775!$OMP&        globaldata_3d_i2,globaldata_3d_i4,globaldata_3d_sp,globaldata_3d_dp,        &
776!$OMP&        globaldata_4d_i2,globaldata_4d_i4,globaldata_4d_sp,globaldata_4d_dp,        &
777!$OMP&        globaldata_1d_i1,globaldata_2d_i1,globaldata_3d_i1,globaldata_4d_i1,        &
778!$OMP&        ntslice,nt,nmax_unlimited,indimnames,dims,patchchunk)
779
780         DO ifile = 1, ndomain
781
782            ncid = inncids(ifile)
783!$OMP CRITICAL
784            CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_size_local', local_sizes ), istop )
785            CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_position_first', start_pos ), istop )
786            CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_halo_size_start', halo_start ), istop )
787            CALL check_nf90( nf90_get_att( ncid, nf90_global, 'DOMAIN_halo_size_end', halo_end ), istop )
788            CALL check_nf90( nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts ), istop )
789            CALL check_nf90( nf90_inquire( ncid, unlimitedDimId = unlimitedDimId ), istop )
790!$OMP END CRITICAL
791
792            ! set defaults for rebuilding so that i is 1st, j 2nd
793            di=1
794            dj=2
795
796            IF( rbdims == 1 ) THEN
797               ! override defaults above and set other variables
798               start_pos(2) = 1
799               local_sizes(2) = outdimlens(3-dimids(2))
800               halo_end(2) = 0
801               halo_start(2) = 0
802               di=rebuild_dims(1)
803               dj=3-di
804            ELSEIF ( iszoom ) THEN
805               start_pos(di) = start_pos(di) - ni_off
806               start_pos(dj) = start_pos(dj) - nj_off
807            ENDIF
808
809!3.3.1 Generate local domain interior sizes from local_sizes and halo sizes
810!      idomain defines the 1st and last interior points in the i direction and
811!      jdomain defines the 1st and last interior points in the j direction
812
813            idomain(1) = 1 + halo_start(di)
814            idomain(2) = local_sizes(di) - halo_end(di)
815            jdomain(1) = 1 + halo_start(dj)
816            jdomain(2) = local_sizes(dj) - halo_end(dj)
817
818!3.3.2 For rbdims or more dimensions put the data array from this input file into the correct
819!      part of the output data array. Assume the first dimensions are those to be rebuilt.
820
821            IF( ndims == 1 ) THEN
822
823               IF( rebuild_dims(1) == 1 ) THEN
824                  dr = di
825                  rdomain = idomain
826               ELSE
827                  dr = dj
828                  rdomain = jdomain
829               ENDIF
830
831              SELECT CASE( xtype )
832                  CASE( NF90_BYTE )
833                     ALLOCATE(localdata_1d_i1(local_sizes(dr)))
834                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_1d_i1 ), istop )
835                     DO jr = rdomain(1), rdomain(2)
836                        globaldata_1d_i1(start_pos(dr) + jr - 1) = localdata_1d_i1(jr)
837                     END DO
838                     DEALLOCATE(localdata_1d_i1)
839                  CASE( NF90_SHORT )
840                     ALLOCATE(localdata_1d_i2(local_sizes(dr)))
841                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_1d_i2 ), istop )
842                     DO jr = rdomain(1), rdomain(2)
843                        globaldata_1d_i2(start_pos(dr) + jr - 1) = localdata_1d_i2(jr)
844                     END DO
845                     DEALLOCATE(localdata_1d_i2)
846                  CASE( NF90_INT )
847                     ALLOCATE(localdata_1d_i4(local_sizes(dr)))
848                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_1d_i4 ), istop )
849                     DO jr = rdomain(1), rdomain(2)
850                        globaldata_1d_i4(start_pos(dr) + jr - 1) = localdata_1d_i4(jr)
851                     END DO
852                     DEALLOCATE(localdata_1d_i4)
853                  CASE( NF90_FLOAT )
854                     ALLOCATE(localdata_1d_sp(local_sizes(dr)))
855                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_1d_sp ), istop )
856                     DO jr = rdomain(1), rdomain(2)
857                        globaldata_1d_sp(start_pos(dr) + jr - 1) = localdata_1d_sp(jr)
858                     END DO
859                     DEALLOCATE(localdata_1d_sp)
860                  CASE( NF90_DOUBLE )
861                     ALLOCATE(localdata_1d_dp(local_sizes(dr)))
862                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_1d_dp ), istop )
863                     DO jr = rdomain(1), rdomain(2)
864                        globaldata_1d_dp(start_pos(dr) + jr - 1) = localdata_1d_dp(jr)
865                     END DO
866                     DEALLOCATE(localdata_1d_dp)
867                  CASE DEFAULT
868                     WRITE(numerr,*) 'Unknown nf90 type: ', xtype
869                     istop = istop + 1
870               END SELECT
871
872            ELSEIF( ndims == 2 ) THEN
873       
874               SELECT CASE( xtype )
875                  CASE( NF90_BYTE )
876                     ALLOCATE(localdata_2d_i1(local_sizes(di),local_sizes(dj)))
877                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_2d_i1 ), istop )
878                        DO jj = jdomain(1), jdomain(2)
879                           DO ji = idomain(1), idomain(2)
880                              globaldata_2d_i1(start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_2d_i1(ji,jj)
881                        END DO
882                     END DO
883                     DEALLOCATE(localdata_2d_i1)
884                  CASE( NF90_SHORT )
885                     ALLOCATE(localdata_2d_i2(local_sizes(di),local_sizes(dj)))
886                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_2d_i2 ), istop )
887                     DO jj = jdomain(1), jdomain(2)
888                        DO ji = idomain(1), idomain(2)
889                           globaldata_2d_i2(start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_2d_i2(ji,jj)
890                        END DO
891                     END DO
892                     DEALLOCATE(localdata_2d_i2)
893                  CASE( NF90_INT )
894                     ALLOCATE(localdata_2d_i4(local_sizes(di),local_sizes(dj)))
895                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_2d_i4 ), istop )
896                     DO jj = jdomain(1), jdomain(2)
897                        DO ji = idomain(1), idomain(2)
898                           globaldata_2d_i4(start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_2d_i4(ji,jj)
899                        END DO
900                     END DO
901                     DEALLOCATE(localdata_2d_i4)
902                  CASE( NF90_FLOAT )
903                     ALLOCATE(localdata_2d_sp(local_sizes(di),local_sizes(dj)))
904                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_2d_sp ), istop )
905                     DO jj = jdomain(1), jdomain(2)
906                        DO ji = idomain(1), idomain(2)
907                           globaldata_2d_sp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_2d_sp(ji,jj)
908                        END DO
909                     END DO
910                     DEALLOCATE(localdata_2d_sp)
911                  CASE( NF90_DOUBLE )
912                     ALLOCATE(localdata_2d_dp(local_sizes(di),local_sizes(dj)))
913                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_2d_dp ), istop )
914                     DO jj = jdomain(1), jdomain(2)
915                        DO ji = idomain(1), idomain(2)
916                           globaldata_2d_dp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_2d_dp(ji,jj)
917                        END DO
918                     END DO
919                     DEALLOCATE(localdata_2d_dp) 
920                  CASE DEFAULT
921                     WRITE(numerr,*) 'Unknown nf90 type: ', xtype
922                     istop = istop + 1
923               END SELECT
924
925            ELSEIF( ndims == 3 ) THEN
926
927               SELECT CASE( xtype )
928                  CASE( NF90_BYTE )
929                     ALLOCATE(localdata_3d_i1(local_sizes(di),local_sizes(dj),indimlens(dimids(3))))
930                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_i1 ), istop )
931!$OMP  PARALLEL DO DEFAULT(NONE) PRIVATE(ji,jj,jk)   &
932!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_3d_i1,localdata_3d_i1,di,dj)
933                     DO jk = 1, indimlens(dimids(3))
934                        DO jj = jdomain(1), jdomain(2)
935                           DO ji = idomain(1), idomain(2)
936                              globaldata_3d_i1(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk) = localdata_3d_i1(ji,jj,jk)
937                           END DO
938                        END DO
939                     END DO
940!$OMP END PARALLEL DO
941                     DEALLOCATE(localdata_3d_i1)
942                  CASE( NF90_SHORT )
943                     ALLOCATE(localdata_3d_i2(local_sizes(di),local_sizes(dj),indimlens(dimids(3))))
944                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_i2 ), istop )
945!$OMP  PARALLEL DO DEFAULT(NONE) PRIVATE(ji,jj,jk)   &
946!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_3d_i2,localdata_3d_i2,di,dj)
947                     DO jk = 1, indimlens(dimids(3))
948                        DO jj = jdomain(1), jdomain(2)
949                           DO ji = idomain(1), idomain(2)
950                              globaldata_3d_i2(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk) = localdata_3d_i2(ji,jj,jk)
951                           END DO
952                        END DO
953                     END DO
954!$OMP END PARALLEL DO
955                     DEALLOCATE(localdata_3d_i2)
956                  CASE( NF90_INT )
957                     ALLOCATE(localdata_3d_i4(local_sizes(di),local_sizes(dj),indimlens(dimids(3))))
958                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_i4 ), istop )
959!$OMP  PARALLEL DO DEFAULT(NONE) PRIVATE(ji,jj,jk)   &
960!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_3d_i4,localdata_3d_i4,di,dj)
961                     DO jk = 1, indimlens(dimids(3))
962                        DO jj = jdomain(1), jdomain(2)
963                           DO ji = idomain(1), idomain(2)
964                              globaldata_3d_i4(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk) = localdata_3d_i4(ji,jj,jk)
965                           END DO
966                        END DO
967                     END DO
968!$OMP END PARALLEL DO
969                     DEALLOCATE(localdata_3d_i4)
970                  CASE( NF90_FLOAT )
971                     ! TG: This if statement is added to check if the 1st dimension is the corners (for lon_bounds) variables
972                     ! TG: Had to add the unsatisfactory check for 'lon' as it failed for diaptr files
973                     ! TG: Would like to find a better assumption for this.
974                     IF ( trim(indimnames(dimids(1))) /= dims(1) .AND. indimnames(dimids(1)) .NE. 'lon' ) THEN
975                        ALLOCATE(localdata_3d_sp(indimlens(dimids(1)),local_sizes(di),local_sizes(dj)))
976                        WRITE(*,*) 'test', ifile, jv, indimlens(dimids(1)),local_sizes(di),local_sizes(dj)
977                        CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_sp ), istop )
978                        WRITE(*,*) 'test2'
979                        DO jj = jdomain(1), jdomain(2)
980                           DO ji = idomain(1), idomain(2)
981                              DO jk = 1, indimlens(dimids(1))
982                                 globaldata_3d_sp(jk, start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_3d_sp(jk,ji,jj)
983                              END DO
984                           END DO
985                        END DO
986                     ELSE
987                        ALLOCATE(localdata_3d_sp(local_sizes(di),local_sizes(dj),indimlens(dimids(3))))
988                        CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_sp ), istop ) 
989!$OMP  PARALLEL DO DEFAULT(NONE) PRIVATE(ji,jj,jk)   &
990!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_3d_sp,localdata_3d_sp,di,dj)
991                        DO jk = 1, indimlens(dimids(3))
992                           DO jj = jdomain(1), jdomain(2)
993                              DO ji = idomain(1), idomain(2)
994                                 globaldata_3d_sp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk) = localdata_3d_sp(ji,jj,jk)
995                              END DO
996                           END DO
997                        END DO
998!$OMP END PARALLEL DO
999                     ENDIF
1000                     DEALLOCATE(localdata_3d_sp)
1001                  CASE( NF90_DOUBLE )
1002                     IF ( trim(indimnames(dimids(1))) /= dims(1) ) THEN
1003                        ALLOCATE(localdata_3d_dp(indimlens(dimids(1)),local_sizes(di),local_sizes(dj)))
1004                        CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_dp ), istop ) 
1005                        DO jj = jdomain(1), jdomain(2)
1006                           DO ji = idomain(1), idomain(2)
1007                              DO jk = 1, indimlens(dimids(1))
1008                                 globaldata_3d_dp(jk, start_pos(di) + ji - 1, start_pos(dj) + jj - 1) = localdata_3d_dp(jk,ji,jj)
1009                              END DO
1010                           END DO
1011                        END DO
1012                     ELSE
1013                        ALLOCATE(localdata_3d_dp(local_sizes(di),local_sizes(dj),indimlens(dimids(3))))
1014                        CALL check_nf90( nf90_get_var( ncid, jv, localdata_3d_dp ), istop ) 
1015!$OMP  PARALLEL DO DEFAULT(NONE) PRIVATE(ji,jj,jk)   &
1016!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_3d_dp,localdata_3d_dp,di,dj)
1017                        DO jk = 1, indimlens(dimids(3))
1018                           DO jj = jdomain(1), jdomain(2)
1019                              DO ji = idomain(1), idomain(2)
1020                                 globaldata_3d_dp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk) = localdata_3d_dp(ji,jj,jk)
1021                              END DO
1022                           END DO
1023                        END DO
1024!$OMP END PARALLEL DO
1025                     ENDIF
1026                     DEALLOCATE(localdata_3d_dp)
1027                  CASE DEFAULT
1028                     WRITE(numerr,*) 'Unknown nf90 type: ', xtype
1029                     istop = istop + 1
1030               END SELECT
1031     
1032            ELSEIF (ndims == 4) THEN
1033       
1034               SELECT CASE( xtype )
1035                  CASE( NF90_BYTE )
1036                     ALLOCATE(localdata_4d_i1(local_sizes(di),local_sizes(dj),               &
1037                         &                     indimlens(dimids(3)),ntslice))
1038                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i1, start=(/1,1,1,nt/) ), istop )
1039!$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
1040!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i1,localdata_4d_i1,di,dj,nt,ntslice)
1041                     DO jl = 1, ntslice
1042!$OMP DO
1043                        DO jk = 1, indimlens(dimids(3))
1044                           DO jj = jdomain(1), jdomain(2)
1045                              DO ji = idomain(1), idomain(2)
1046                                 globaldata_4d_i1(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk, jl) = localdata_4d_i1(ji,jj,jk,jl)
1047                              END DO
1048                           END DO
1049                        END DO
1050!$OMP END DO nowait
1051                     END DO
1052!$OMP END PARALLEL
1053                     DEALLOCATE(localdata_4d_i1)
1054                  CASE( NF90_SHORT )
1055                     ALLOCATE(localdata_4d_i2(local_sizes(di),local_sizes(dj),               &
1056                        &                     indimlens(dimids(3)),ntslice))
1057                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i2, start=(/1,1,1,nt/) ), istop )
1058!$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
1059!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i2,localdata_4d_i2,di,dj,nt,ntslice)
1060                     DO jl = 1, ntslice
1061!$OMP DO
1062                        DO jk = 1, indimlens(dimids(3))
1063                           DO jj = jdomain(1), jdomain(2) 
1064                              DO ji = idomain(1), idomain(2)
1065                                 globaldata_4d_i2(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk, jl) = localdata_4d_i2(ji,jj,jk,jl)
1066                              END DO
1067                           END DO
1068                        END DO
1069!$OMP END DO nowait
1070                     END DO
1071!$OMP END PARALLEL
1072                     DEALLOCATE(localdata_4d_i2)
1073                  CASE( NF90_INT )
1074                     ALLOCATE(localdata_4d_i4(local_sizes(di),local_sizes(dj),               &
1075                        &                     indimlens(dimids(3)),ntslice))
1076                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_i4, start=(/1,1,1,nt/) ), istop )
1077!$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
1078!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_i4,localdata_4d_i4,di,dj,nt,ntslice)
1079                     DO jl = 1, ntslice
1080!$OMP DO
1081                        DO jk = 1, indimlens(dimids(3))
1082                           DO jj = jdomain(1), jdomain(2) 
1083                              DO ji = idomain(1), idomain(2)
1084                                 globaldata_4d_i4(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk, jl) = localdata_4d_i4(ji,jj,jk,jl)
1085                              END DO
1086                           END DO
1087                        END DO
1088!$OMP END DO nowait
1089                     END DO
1090!$OMP END PARALLEL
1091                     DEALLOCATE(localdata_4d_i4)
1092                  CASE( NF90_FLOAT )
1093                     ALLOCATE(localdata_4d_sp(local_sizes(di),local_sizes(dj),               &
1094                        &                     indimlens(dimids(3)),ntslice))
1095                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_sp, start=(/1,1,1,nt/) ), istop )
1096!$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
1097!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_sp,localdata_4d_sp,di,dj,nt,ntslice)
1098                     DO jl = 1, ntslice
1099!$OMP DO
1100                        DO jk = 1, indimlens(dimids(3))
1101                           DO jj = jdomain(1), jdomain(2) 
1102                              DO ji = idomain(1), idomain(2)
1103                                 globaldata_4d_sp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk, jl) = localdata_4d_sp(ji,jj,jk,jl)
1104                              END DO
1105                           END DO
1106                        END DO
1107!$OMP END DO nowait
1108                     END DO
1109!$OMP END PARALLEL
1110                     DEALLOCATE(localdata_4d_sp)
1111                  CASE( NF90_DOUBLE )
1112                     ALLOCATE(localdata_4d_dp(local_sizes(di),local_sizes(dj),               &
1113                        &                     indimlens(dimids(3)),ntslice))
1114                     CALL check_nf90( nf90_get_var( ncid, jv, localdata_4d_dp, start=(/1,1,1,nt/) ), istop )
1115!$OMP  PARALLEL DEFAULT(NONE) PRIVATE(ji,jj,jk,jl)   &
1116!$OMP& SHARED(idomain,jdomain,indimlens,dimids,start_pos,globaldata_4d_dp,localdata_4d_dp,di,dj,nt,ntslice)
1117                     DO jl = 1, ntslice
1118!$OMP DO
1119                        DO jk = 1, indimlens(dimids(3))
1120                           DO jj = jdomain(1), jdomain(2) 
1121                              DO ji = idomain(1), idomain(2)
1122                                 globaldata_4d_dp(start_pos(di) + ji - 1, start_pos(dj) + jj - 1, jk, jl) = localdata_4d_dp(ji,jj,jk,jl)
1123                              END DO
1124                           END DO
1125                        END DO
1126!$OMP END DO nowait
1127                     END DO
1128!$OMP END PARALLEL
1129                     DEALLOCATE(localdata_4d_dp)
1130                  CASE DEFAULT
1131                     WRITE(numerr,*) 'Unknown nf90 type: ', xtype
1132                     istop = istop + 1
1133               END SELECT
1134
1135            ENDIF ! l_noRebuild false
1136
1137!3.4 Work out if the valid_min and valid_max attributes exist for this variable.
1138!    If they do then calculate the extrema over all input files.
1139
1140            DO attid = 1, natts
1141               CALL check_nf90( nf90_inq_attname( ncid, jv, attid, attname ), istop )
1142               IF( INDEX( attname, "valid_min" ) == 1 ) THEN
1143                  CALL check_nf90( nf90_get_att( ncid, jv, attname, InMin), istop )
1144                  l_valid = .true.
1145               ENDIF
1146               IF( INDEX( attname, "valid_max" ) == 1 ) THEN
1147                  CALL check_nf90( nf90_get_att( ncid, jv, attname, InMax ), istop )
1148                  l_valid = .true.
1149               ENDIF
1150            END DO
1151
1152            IF (l_valid) THEN
1153!$OMP CRITICAL
1154               IF( InMin < ValMin ) ValMin = InMin
1155               IF( InMax > ValMax ) ValMax = InMax
1156!$OMP END CRITICAL
1157            ENDIF
1158
1159!3.5 Abort if failure and only 1 thread
1160
1161            IF( nthreads == 1 .AND. istop /= nf90_noerr )  THEN
1162               WRITE(numerr,*) '*** NEMO rebuild failed! ***'
1163               STOP 5
1164            ENDIF
1165       
1166         END DO  ! loop over files
1167!$OMP END PARALLEL DO
1168
1169!3.6 Abort if any of the OMP threads failed
1170         IF( istop /= nf90_noerr )  THEN
1171            WRITE(numerr,*) '*** NEMO rebuild failed! ***'
1172            STOP 5
1173         ENDIF
1174
1175      ENDIF ! ndims > 2
1176
1177!---------------------------------------------------------------------------
1178!4. Write data to output file
1179
1180      IF (l_verbose) WRITE(numout,*) 'Writing variable '//TRIM(varname)//'...'
1181
1182!4.1 If the valid min and max attributes exist then update them in the file
1183
1184      IF( l_valid ) THEN
1185         CALL check_nf90( nf90_put_att( outid, jv, "valid_min", ValMin ) )
1186         CALL check_nf90( nf90_put_att( outid, jv, "valid_max", ValMax ) )
1187      ENDIF
1188
1189!4.2 Write the data to the output file depending on how many dimensions it has
1190
1191      IF( ndims == 0 ) THEN
1192
1193         SELECT CASE( xtype )
1194            CASE( NF90_BYTE )
1195               CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_i1 ) )
1196            CASE( NF90_SHORT )
1197               CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_i2 ) )
1198            CASE( NF90_INT )
1199               CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_i4 ) )
1200            CASE( NF90_FLOAT )
1201               CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_sp ) )
1202            CASE( NF90_DOUBLE )
1203               CALL check_nf90( nf90_put_var( outid, jv, globaldata_0d_dp ) )
1204            CASE DEFAULT   
1205               WRITE(numerr,*) '0d Unknown nf90 type: ', xtype
1206               STOP 4
1207         END SELECT
1208
1209      ELSEIF( ndims == 1 ) THEN
1210
1211         SELECT CASE( xtype )
1212            CASE( NF90_BYTE )
1213               CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_i1 ) )
1214               DEALLOCATE(globaldata_1d_i1)
1215            CASE( NF90_SHORT )
1216               CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_i2 ) )
1217               DEALLOCATE(globaldata_1d_i2)
1218            CASE( NF90_INT )
1219               CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_i4 ) )
1220               DEALLOCATE(globaldata_1d_i4)
1221            CASE( NF90_FLOAT )
1222               CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_sp ) )
1223               DEALLOCATE(globaldata_1d_sp)
1224            CASE( NF90_DOUBLE )
1225               CALL check_nf90( nf90_put_var( outid, jv, globaldata_1d_dp ) )
1226               DEALLOCATE(globaldata_1d_dp)
1227            CASE DEFAULT   
1228               WRITE(numerr,*) '1d Unknown nf90 type: ', xtype
1229               STOP 4
1230         END SELECT
1231
1232      ELSEIF( ndims == 2 ) THEN 
1233     
1234         SELECT CASE( xtype )   
1235            CASE( NF90_BYTE )                   
1236               CALL check_nf90( nf90_put_var( outid, jv, globaldata_2d_i1 ) )
1237               DEALLOCATE(globaldata_2d_i1)
1238            CASE( NF90_SHORT )                   
1239               CALL check_nf90( nf90_put_var( outid, jv, globaldata_2d_i2 ) )
1240               DEALLOCATE(globaldata_2d_i2)
1241            CASE( NF90_INT )                             
1242               CALL check_nf90( nf90_put_var( outid, jv, globaldata_2d_i4 ) )
1243               DEALLOCATE(globaldata_2d_i4)
1244            CASE( NF90_FLOAT )                             
1245               CALL check_nf90( nf90_put_var( outid, jv, globaldata_2d_sp ) )
1246               DEALLOCATE(globaldata_2d_sp)
1247            CASE( NF90_DOUBLE )                                         
1248               CALL check_nf90( nf90_put_var( outid, jv, globaldata_2d_dp ) )
1249               DEALLOCATE(globaldata_2d_dp)
1250            CASE DEFAULT   
1251               WRITE(numerr,*) '2d Unknown nf90 type: ', xtype
1252               STOP 4
1253         END SELECT     
1254                     
1255      ELSEIF( ndims == 3 ) THEN
1256     
1257         SELECT CASE( xtype ) 
1258            CASE( NF90_BYTE )                   
1259               CALL check_nf90( nf90_put_var( outid, jv, globaldata_3d_i1 ) )
1260               DEALLOCATE(globaldata_3d_i1)
1261            CASE( NF90_SHORT )                   
1262               CALL check_nf90( nf90_put_var( outid, jv, globaldata_3d_i2 ) )
1263               DEALLOCATE(globaldata_3d_i2)
1264            CASE( NF90_INT )                             
1265               CALL check_nf90( nf90_put_var( outid, jv, globaldata_3d_i4 ) )
1266               DEALLOCATE(globaldata_3d_i4)
1267            CASE( NF90_FLOAT )                             
1268               CALL check_nf90( nf90_put_var( outid, jv, globaldata_3d_sp ) )
1269               DEALLOCATE(globaldata_3d_sp)
1270            CASE( NF90_DOUBLE )                                         
1271               CALL check_nf90( nf90_put_var( outid, jv, globaldata_3d_dp ) )
1272               DEALLOCATE(globaldata_3d_dp)
1273            CASE DEFAULT   
1274               WRITE(numerr,*) '3d Unknown nf90 type: ', xtype
1275               STOP 4
1276         END SELECT     
1277   
1278      ELSEIF( ndims == 4 ) THEN
1279
1280         SELECT CASE( xtype )   
1281            CASE( NF90_BYTE )                   
1282               CALL check_nf90( nf90_put_var( outid, jv, globaldata_4d_i1, start=(/1,1,1,nt/) ) )
1283               DEALLOCATE(globaldata_4d_i1)
1284            CASE( NF90_SHORT )                   
1285               CALL check_nf90( nf90_put_var( outid, jv, globaldata_4d_i2, start=(/1,1,1,nt/) ) )
1286               DEALLOCATE(globaldata_4d_i2)
1287            CASE( NF90_INT )                             
1288               CALL check_nf90( nf90_put_var( outid, jv, globaldata_4d_i4, start=(/1,1,1,nt/) ) )
1289               DEALLOCATE(globaldata_4d_i4)
1290            CASE( NF90_FLOAT )                             
1291               CALL check_nf90( nf90_put_var( outid, jv, globaldata_4d_sp, start=(/1,1,1,nt/) ) )
1292               DEALLOCATE(globaldata_4d_sp)
1293            CASE( NF90_DOUBLE )                                         
1294               CALL check_nf90( nf90_put_var( outid, jv, globaldata_4d_dp, start=(/1,1,1,nt/) ) )
1295               DEALLOCATE(globaldata_4d_dp)
1296            CASE DEFAULT   
1297               WRITE(numerr,*) '4d Unknown nf90 type: ', xtype
1298               STOP 4
1299         END SELECT     
1300         ! why only for big data set, test the cost.
1301         CALL check_nf90( nf90_sync( outid ) )    ! flush buffers to disk after writing big 4D datasets
1302   
1303      ENDIF
1304
1305         nt = nt + ntslice
1306
1307      END DO ! WHILE loop
1308   
1309   END DO  ! loop over variables
1310
1311!---------------------------------------------------------------------------
1312!5. Close files
1313
1314!5.1 Close input files
1315
1316   IF (l_verbose) WRITE(numout,*) 'Closing input files...'
1317   DO ifile = 1, ndomain
1318      ncid = inncids(ifile)
1319      CALL check_nf90( nf90_close( ncid ) )
1320   END DO
1321
1322!5.2 Close output file
1323
1324   IF (l_verbose) WRITE(numout,*) 'Closing output file...'
1325   CALL check_nf90( nf90_close( outid ) )
1326
1327   IF (l_verbose) WRITE(numout,*) 'NEMO rebuild completed successfully'
1328   IF (l_verbose) WRITE(numout,*)
1329
1330CONTAINS
1331
1332
1333   SUBROUTINE check_nf90(status, errorFlag)
1334   !---------------------------------------------------------------------
1335   !  Checks return code from nf90 library calls and warns if needed
1336   !  If errorFlag is present then it just increments this flag (OMP use)
1337   !
1338   !---------------------------------------------------------------------
1339      INTEGER, INTENT(IN   ) :: status
1340      INTEGER, INTENT(INOUT), OPTIONAL :: errorFlag
1341   !---------------------------------------------------------------------
1342
1343      IF( status /= nf90_noerr ) THEN
1344         WRITE(numerr,*) 'ERROR! : '//TRIM(nf90_strerror(status))
1345         IF( PRESENT( errorFlag ) ) THEN
1346            errorFlag = errorFlag + status
1347         ELSE
1348            WRITE(numerr,*) "*** NEMO rebuild failed ***"
1349            WRITE(numerr,*)
1350            STOP 5
1351         ENDIF
1352      ENDIF
1353
1354   END SUBROUTINE check_nf90
1355
1356
1357END PROGRAM rebuild_nemo