source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/TOOLS/BDY_TOOLS/src/bdy_reorder.f90 @ 3239

Last change on this file since 3239 was 3239, checked in by davestorkey, 9 years ago

Add bdy_reorder utility to TOOLS directory.

File size: 26.3 KB
Line 
1PROGRAM bdy_reorder
2
3!===============================================================================
4! A routine to reorder old BDY data files to make them compatible with NEMO 3.4.
5!
6! This routine has 2 modes:
7!   1. If no template file is given then it will re-order the data in the input
8!      file so that it is in order of increasing nbr.
9!   2. If a template file is given it will re-order the data in the input file
10!      be consistent with the order of the data in the template file. This is
11!      useful for making old barotropic and baroclinic data files consistent.
12!      For older versions of NEMO the barotropic boundary data did not have to
13!      be in the same order as the baroclinic boundary data. The rimwidth
14!      value can be different in the input file and template file; in this case
15!      the routine will just re-order the data that the two files have in common.
16!
17! The routine is mainly for re-ordering BDY data files, but can also be used to
18! re-order BDY coordinate files if ln_coordinates is set to .true.
19!
20! Author: Dave Storkey  Aug 2011
21! Bug notifications etc to:  dave.storkey@metoffice.gov.uk
22!
23!===============================================================================
24 
25      USE netcdf
26
27      IMPLICIT NONE
28
29      INTEGER,PARAMETER         :: numnam=11
30      INTEGER,PARAMETER         :: sp=SELECTED_REAL_KIND(6,37)
31      INTEGER,PARAMETER         :: dp=SELECTED_REAL_KIND(12,307)
32      INTEGER                   :: chunksize = 32000000
33
34      INTEGER                   :: jpbgrd, iostat, ncid_in, ncid_out, ncid_template, nbrid(4), nbrid_template(4)
35      INTEGER                   :: dimids(10), unlimitedDimId
36      INTEGER                   :: xtype, dimid, varid, ndims, ndims_var, nvars, natts, nblen(4), nblen_template(4), dimlen(10)
37      INTEGER                   :: lenvar(10), igrid, ib, ib1, ir, idim, jgrid, jv, icount, attid, len1, idepth, itime
38      INTEGER                   :: idim_time, idim_xb, idim_yb, idim_depth
39      INTEGER                   :: nbr_min, nbr_max, nbr_min_template, nbr_max_template, strlen, nbr_match
40
41      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: nbr, imap
42      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: nbi, nbj
43      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: nbi_template
44      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: nbj_template
45      INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: nbr_template
46
47      INTEGER, ALLOCATABLE, DIMENSION(:)     :: nbi_extract, nbj_extract
48
49      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:,:)  :: varin_int, varout_int
50      REAL(sp), ALLOCATABLE, DIMENSION(:,:,:,:)  :: varin_float, varout_float
51      REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:)  :: varin_dble, varout_dble
52
53      LOGICAL                   :: ln_coordinates
54      LOGICAL                   :: ln_template
55
56      CHARACTER(LEN=nf90_max_name)        :: attname, dimname, varname, time, date, zone, timestamp
57      CHARACTER(LEN=nf90_max_name)        :: file_in, file_out, file_template, nbi_name, nbj_name, nbr_name
58      CHARACTER(LEN=1),DIMENSION(4)       :: cgrid
59      CHARACTER(LEN=1)                    :: end_letter
60
61      NAMELIST/nam_bdy_reorder/ file_in, file_out, file_template, ln_coordinates, nbr_match
62
63!=============================================================================
64
65!-----------------------------------------------------------------------------
66! 0. Read namelist and initialise parameters
67!-----------------------------------------------------------------------------
68
69      file_in = ''
70      file_out = ''
71      file_template = ''
72      ln_coordinates = .false.
73      nbr_match = 0
74
75      OPEN( UNIT=numnam, FILE='nam_bdy_reorder', FORM='FORMATTED', STATUS='OLD' )
76      READ( numnam, nam_bdy_reorder )
77      CLOSE( UNIT=numnam )
78
79      IF( ln_coordinates ) THEN
80         WRITE(*,*) 'Input file is a coordinates.bdy.nc file.'
81         jpbgrd = 4
82      ELSE
83          WRITE(*,*) 'Input file is a boundary data file.'
84         jpbgrd = 1
85      ENDIF
86 
87      IF( LEN_TRIM(file_template) > 0 ) THEN
88         WRITE(*,*) 'Reordering data to match template file.'
89         ln_template = .true.
90      ELSE
91         WRITE(*,*) 'Reordering data in order of increasing nbr values.'
92         ln_template = .false.
93      ENDIF
94
95      cgrid= (/'t','u','v','f'/)
96
97!-----------------------------------------------------------------------------
98! 1. Read in nbr variables from file and create mapping.
99!-----------------------------------------------------------------------------
100
101! 1.1 Open input files
102
103      iostat = nf90_open( TRIM(file_in), nf90_nowrite, ncid_in )
104      IF (iostat /= nf90_noerr) THEN
105        WRITE(6,*) TRIM(nf90_strerror(iostat))
106        STOP
107      ENDIF
108      iostat = nf90_inquire( ncid_in, ndims, nvars, natts )
109
110      IF( ln_template ) THEN
111         iostat = nf90_open( TRIM(file_template), nf90_nowrite, ncid_template )
112         IF (iostat /= nf90_noerr) THEN
113           WRITE(6,*) TRIM(nf90_strerror(iostat))
114           STOP
115         ENDIF
116      ENDIF
117
118      DO igrid=1,jpbgrd
119
120! 1.2 Find dimensions of data in input files and allocate arrays
121
122         IF( ln_coordinates ) THEN
123            nbr_name = 'nbr'//cgrid(igrid)
124         ELSE
125            nbr_name = 'nbrdta'
126         ENDIF
127
128         iostat = nf90_inq_varid( ncid_in, nbr_name , nbrid(igrid) )
129         IF (iostat /= nf90_noerr) THEN
130            WRITE(6,*) TRIM(nf90_strerror(iostat))
131            STOP
132         ENDIF
133         iostat = nf90_inquire_variable( ncid_in, nbrid(igrid), ndims=ndims, dimids=dimids )
134         IF( ndims .ne. 2 ) THEN
135            WRITE(*,*) 'ERROR : ',TRIM(nbr_name),' does not have exactly 2 dimensions.'
136            WRITE(*,*) '        in file ',TRIM(file_in)
137            STOP
138         ENDIF
139         iostat = nf90_inquire_dimension( ncid_in, dimids(1), len=nblen(igrid))
140         iostat = nf90_inquire_dimension( ncid_in, dimids(2), len=len1)
141         WRITE(*,*) 'nblen(igrid), len1 : ',nblen(igrid), len1
142         IF( len1 .ne. 1 ) THEN
143            WRITE(*,*) 'ERROR : second dimension of ',TRIM(nbr_name),' does not have length 1.'
144            WRITE(*,*) '        in file ',TRIM(file_in)
145            WRITE(*,*) '        len1 = ',len1
146            STOP
147         ENDIF
148
149         IF( ln_template ) THEN
150
151            iostat = nf90_inq_varid( ncid_template, nbr_name , nbrid_template(igrid) )
152            IF (iostat /= nf90_noerr) THEN
153               WRITE(6,*) TRIM(nf90_strerror(iostat))
154               STOP
155            ENDIF
156            iostat = nf90_inquire_variable( ncid_template, nbrid_template(igrid), ndims=ndims, dimids=dimids )
157            IF( ndims .ne. 2 ) THEN
158               WRITE(*,*) 'ERROR : ',TRIM(nbr_name),' does not have exactly 2 dimensions.'
159               WRITE(*,*) '        in file ',TRIM(file_template)
160               STOP
161            ENDIF
162            iostat = nf90_inquire_dimension( ncid_template, dimids(1), len=nblen_template(igrid))
163            iostat = nf90_inquire_dimension( ncid_template, dimids(2), len=len1)
164            WRITE(*,*) 'nblen_template(igrid), len1 : ',nblen_template(igrid), len1
165            IF( len1 .ne. 1 ) THEN
166               WRITE(*,*) 'ERROR : second dimension of ',TRIM(nbr_name),' does not have length 1.'
167               WRITE(*,*) '        in file ',TRIM(file_template)
168               WRITE(*,*) '        len1 = ',len1
169               STOP
170            ENDIF
171
172         ENDIF ! ln_template
173
174      ENDDO  ! jpbgrd
175
176      ALLOCATE( nbr(MAXVAL(nblen),1), imap(MAXVAL(nblen),jpbgrd) )
177      IF( ln_template ) THEN
178         ALLOCATE( nbi(MAXVAL(nblen),1), nbj(MAXVAL(nblen),1) )
179         ALLOCATE( nbi_template(MAXVAL(nblen_template),1) )
180         ALLOCATE( nbj_template(MAXVAL(nblen_template),1) )
181         ALLOCATE( nbr_template(MAXVAL(nblen_template),1) )
182      ENDIF
183
184! 1.3 Read in nbr variables and generate mapping.
185
186      nbr(:,:) = -1
187      imap(:,:) = -1   
188      IF( ln_template ) THEN
189         nbi(:,:) = -1
190         nbj(:,:) = -1
191         nbi_template(:,:) = -1
192         nbj_template(:,:) = -1
193         nbr_template(:,:) = -1
194      ENDIF
195
196      WRITE(*,*) '>>> Generating map.'
197      DO igrid=1,jpbgrd
198
199         iostat = nf90_get_var( ncid_in, nbrid(igrid), nbr(1:nblen(igrid),1:1) ) 
200         IF (iostat /= nf90_noerr) THEN
201           WRITE(6,*) TRIM(nf90_strerror(iostat))
202           STOP
203         ENDIF
204
205         nbr_min = MINVAL( nbr(1:nblen(igrid),1) )
206         nbr_max = MAXVAL( nbr(1:nblen(igrid),1) )
207         IF( nbr_min .ne. 1 ) THEN
208            WRITE(*,*) 'Something wrong with input file ',file_in
209            WRITE(*,*) 'MIN(nbr) /= 1'
210            STOP
211         ENDIF
212
213         IF( ln_template ) THEN
214
215            IF( ln_coordinates ) THEN
216               nbi_name = 'nbi'//cgrid(igrid)
217               nbj_name = 'nbj'//cgrid(igrid)
218            ELSE
219               nbi_name = 'nbidta'           
220               nbj_name = 'nbjdta'           
221            ENDIF
222
223            ! read in nbi, nbj values from input file
224            iostat = nf90_inq_varid( ncid_in, nbi_name, varid ) 
225            IF (iostat == nf90_noerr) iostat = nf90_get_var( ncid_in, varid, nbi(1:nblen(igrid),1:1) ) 
226            IF (iostat /= nf90_noerr) THEN
227              WRITE(6,*) TRIM(nf90_strerror(iostat))
228              STOP
229            ENDIF
230            iostat = nf90_inq_varid( ncid_in, nbj_name, varid ) 
231            IF (iostat == nf90_noerr) iostat = nf90_get_var( ncid_in, varid, nbj(1:nblen(igrid),1:1) ) 
232            IF (iostat /= nf90_noerr) THEN
233              WRITE(6,*) TRIM(nf90_strerror(iostat))
234              STOP
235            ENDIF
236
237            ! read in nbi, nbj, nbr values from template file
238            iostat = nf90_inq_varid( ncid_template, nbi_name, varid ) 
239            IF (iostat == nf90_noerr) iostat = nf90_get_var( ncid_template, varid, nbi_template(1:nblen_template(igrid),1:1) ) 
240            IF (iostat /= nf90_noerr) THEN
241              WRITE(6,*) TRIM(nf90_strerror(iostat))
242              STOP
243            ENDIF
244            iostat = nf90_inq_varid( ncid_template, nbj_name, varid ) 
245            IF (iostat == nf90_noerr) iostat = nf90_get_var( ncid_template, varid, nbj_template(1:nblen_template(igrid),1:1) ) 
246            IF (iostat /= nf90_noerr) THEN
247              WRITE(6,*) TRIM(nf90_strerror(iostat))
248              STOP
249            ENDIF
250            iostat = nf90_get_var( ncid_template, nbrid_template(igrid), nbr_template(1:nblen_template(igrid),1:1) ) 
251            IF (iostat /= nf90_noerr) THEN
252              WRITE(6,*) TRIM(nf90_strerror(iostat))
253              STOP
254            ENDIF
255
256            nbr_min_template = MINVAL( nbr_template(1:nblen_template(igrid),1) )
257            nbr_max_template = MAXVAL( nbr_template(1:nblen_template(igrid),1) )
258            IF( nbr_min_template .ne. 1 ) THEN
259               WRITE(*,*) 'Something wrong with template file ',file_template
260               WRITE(*,*) 'MIN(nbr) /= 1'
261               STOP
262            ENDIF
263
264            IF( nbr_match .lt. 1 ) nbr_match = MIN(nbr_max, nbr_max_template)
265
266            ! initialise imap to the identity mapping
267            DO ib = 1, nblen(igrid)
268               imap(ib, igrid) = ib
269            ENDDO
270
271            ! allocate "extract" arrays
272            icount = 0
273            DO ib = 1, nblen_template(igrid)
274               IF( nbr_template(ib,1) .eq. 1 ) THEN
275                  icount = icount+1
276               ENDIF
277            ENDDO
278            ALLOCATE( nbi_extract(icount), nbj_extract(icount) )
279
280            DO ir = 1, nbr_match
281
282               ! extract values from template array for this ir value
283               icount = 0
284               DO ib = 1, nblen_template(igrid)
285                  IF( nbr_template(ib,1) .eq. ir ) THEN
286                     icount = icount+1
287                     nbi_extract(icount) = nbi_template(ib,1)               
288                     nbj_extract(icount) = nbj_template(ib,1)               
289                  ENDIF
290               ENDDO
291
292               ! work out the mapping for this ir value
293               icount = 1
294               DO ib = 1, nblen(igrid)
295                  IF( nbr(ib,1) .eq. ir ) THEN                 
296                     DO ib1 = 1, nblen(igrid)
297                        IF( nbi(ib1,1) .eq. nbi_extract(icount) .and. &
298                       &    nbj(ib1,1) .eq. nbj_extract(icount)       ) THEN
299                           imap(ib,igrid) = ib1
300                           icount = icount + 1
301                           EXIT
302                        ENDIF                       
303                     ENDDO ! ib1
304                  ENDIF           
305               ENDDO ! ib
306
307            ENDDO ! ir
308           
309            DEALLOCATE( nbi_extract, nbj_extract )
310
311         ELSE
312
313            icount = 0
314            DO ir = nbr_min, nbr_max
315               DO ib = 1, nblen(igrid)
316                  IF( nbr(ib,1) .eq. ir ) THEN
317                     icount = icount + 1
318                     imap(icount,igrid) = ib
319                  ENDIF           
320               ENDDO
321            ENDDO
322
323         ENDIF
324     
325      ENDDO  ! jpbgrd
326
327!-----------------------------------------------------------------------------
328! 2. Open output file and copy dimensions and attributes across
329!-----------------------------------------------------------------------------
330
331      iostat = nf90_inquire( ncid_in, ndims, nvars, natts )
332
333! 2.1 Create the output file 
334
335      WRITE(*,*) '>>> Initialising output file.'
336      iostat = nf90_create( TRIM(file_out), nf90_64bit_offset, ncid_out, chunksize=chunksize)
337
338! 2.2 Copy the dimensions into the output file.
339
340      iostat = nf90_inquire( ncid_in, unlimitedDimId = unlimitedDimId )
341      dimlen(:) = 1
342      DO idim = 1, ndims
343         iostat = nf90_inquire_dimension(ncid_in, idim, dimname, dimlen(idim))
344         IF (idim == unlimitedDimId) THEN
345            iostat = nf90_def_dim( ncid_out, dimname, nf90_unlimited, dimid)   
346            idim_time = idim
347         ELSE
348            iostat = nf90_def_dim( ncid_out, dimname, dimlen(idim), dimid)
349            IF( INDEX(dimname,'x') .gt. 0 ) THEN
350               idim_xb = idim
351            ELSE IF( INDEX(dimname,'y') .gt. 0 ) THEN
352               idim_yb = idim
353            ELSE IF( INDEX(dimname,'depth') .gt. 0 .or. INDEX(dimname,'z') .gt. 0 ) THEN
354               idim_depth = idim
355            ELSE
356               WRITE(*,*) 'ERROR: Unrecognised dimension : ',dimname
357               STOP
358            ENDIF           
359         ENDIF
360      END DO
361     
362
363! 2.2 Copy the global attributes into the output file.
364!     Also need to change the file_name attribute and the TimeStamp attribute.
365
366      DO attid = 1, natts
367         iostat = nf90_inq_attname( ncid_in, nf90_global, attid, attname )
368         WRITE(6,*)'>>> Copying attribute '//TRIM(attname)//' into destination file...'
369         iostat = nf90_copy_att( ncid_in, nf90_global, attname, ncid_out, nf90_global ) 
370      END DO
371      iostat = nf90_put_att( ncid_out, nf90_global, "file_name", TRIM(file_out))
372      CALL DATE_AND_TIME ( date=date, time=time, zone=zone )
373      timestamp = date(7:8) // "/" // date(5:6) // "/" // date(1:4) // " " // &
374                  time(1:2) // ":" // time(3:4) // ":" // time(5:6) // " " // &
375                  zone 
376      iostat = nf90_put_att( ncid_out, nf90_global, "TimeStamp", timestamp)
377 
378! 2.3 Copy the variable definitions and attributes into the output file.
379
380      DO jv = 1, nvars
381         iostat = nf90_inquire_variable( ncid_in, jv, varname, xtype, ndims_var, dimids, natts )
382         iostat = nf90_def_var( ncid_out, varname, xtype, dimids(1:ndims_var), varid )
383            IF (natts > 0) THEN
384            DO attid = 1, natts
385               iostat = nf90_inq_attname(ncid_in, varid, attid, attname)
386               iostat = nf90_copy_att( ncid_in, varid, attname, ncid_out, varid )   
387            END DO
388         ENDIF
389      END DO
390 
391! 2.4 End definitions in output file
392
393      iostat = nf90_enddef( ncid_out )   
394
395!-----------------------------------------------------------------------------
396! 3. Read in variables from input file, re-order and write to output file
397!-----------------------------------------------------------------------------
398
399     IF( ln_coordinates ) THEN
400        ALLOCATE( varin_int(MAXVAL(nblen),1,1,1), varout_int(MAXVAL(nblen),1,1,1) )
401        ALLOCATE( varin_float(MAXVAL(nblen),1,1,1), varout_float(MAXVAL(nblen),1,1,1) )
402        ALLOCATE( varin_dble(MAXVAL(nblen),1,1,1), varout_dble(MAXVAL(nblen),1,1,1) )
403     ELSE
404        SELECT CASE( ndims ) 
405           CASE( 2 )
406              ALLOCATE(  varin_int(dimlen(idim_xb),dimlen(idim_yb),1,1) )
407              ALLOCATE( varout_int(dimlen(idim_xb),dimlen(idim_yb),1,1) )
408              ALLOCATE(  varin_float(dimlen(idim_xb),dimlen(idim_yb),1,1) ) 
409              ALLOCATE( varout_float(dimlen(idim_xb),dimlen(idim_yb),1,1) )
410              ALLOCATE(  varin_dble(dimlen(idim_xb),dimlen(idim_yb),1,1) ) 
411              ALLOCATE( varout_dble(dimlen(idim_xb),dimlen(idim_yb),1,1) )
412           CASE ( 3 )
413              ALLOCATE(  varin_int(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) )
414              ALLOCATE( varout_int(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) )
415              ALLOCATE(  varin_float(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) ) 
416              ALLOCATE( varout_float(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) )
417              ALLOCATE(  varin_dble(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) ) 
418              ALLOCATE( varout_dble(dimlen(idim_xb),dimlen(idim_yb),1,dimlen(idim_time)) )
419           CASE ( 4 )
420              ALLOCATE(  varin_int(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) ) 
421              ALLOCATE( varout_int(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) )
422              ALLOCATE(  varin_float(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) ) 
423              ALLOCATE( varout_float(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) )
424              ALLOCATE(  varin_dble(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) ) 
425              ALLOCATE( varout_dble(dimlen(idim_xb),dimlen(idim_yb),dimlen(idim_depth),dimlen(idim_time)) )
426           CASE DEFAULT
427              WRITE(*,*) 'ERROR : Can only cope with 2, 3 or 4 dimensions in the boundary data files.'
428              WRITE(*,*) '        This file appears to have ',ndims,' dimensions.'
429              STOP
430        END SELECT
431     ENDIF
432
433     DO jv = 1, nvars
434       iostat = nf90_inquire_variable( ncid_in, jv, varname, xtype, ndims_var, dimids )
435       DO idim = 1, ndims_var
436          lenvar(idim) = dimlen(dimids(idim))
437       ENDDO
438
439       IF( ndims_var .eq. 1 ) THEN
440          WRITE(*,*) '>>> Copying coordinate variable ',TRIM(varname)
441       ELSE
442          WRITE(*,*) '>>> Reordering variable ',TRIM(varname)
443       ENDIF
444       ! Error check here
445
446       IF( ln_coordinates ) THEN
447          strlen = len_trim(varname)
448          end_letter = varname(strlen:strlen)
449          jgrid = -1
450          DO igrid = 1,4
451             IF( end_letter .eq. cgrid(igrid) ) jgrid = igrid
452          END DO
453          IF( jgrid .lt. 0 ) THEN
454             WRITE(*,*) 'ERROR : Could not identify grid for variable ',TRIM(varname)
455             WRITE(*,*) ' varname : ',TRIM(varname),'!'
456             WRITE(*,*) ' strlen : ',strlen
457             WRITE(*,*) ' end_letter : ',end_letter
458             WRITE(*,*) ' cgrid(1) : ',cgrid(1)
459             WRITE(*,*) ' cgrid(2) : ',cgrid(2)
460             WRITE(*,*) ' cgrid(3) : ',cgrid(3)
461             WRITE(*,*) ' cgrid(4) : ',cgrid(4)
462             STOP
463          ENDIF
464       ELSE
465          jgrid=1
466       ENDIF
467           
468       SELECT CASE(xtype)
469
470          CASE( NF90_INT ) 
471             SELECT CASE(ndims_var)
472                CASE( 1 )
473                   WRITE(*,*) 'This is a 1D integer.'
474                   ! Assume this is a depth or time coordinate variable and copy across unchanged.
475                   iostat = nf90_get_var( ncid_in, jv, varin_int(1:lenvar(1),1,1,1) )
476                   iostat = nf90_put_var( ncid_out, jv, varout_int(1:lenvar(1),1,1,1) ) 
477                CASE( 2 )
478                   WRITE(*,*) 'This is a 2D integer.'
479                   iostat = nf90_get_var( ncid_in, jv, varin_int(1:lenvar(1),1:1,1,1) )
480                   DO ib = 1,lenvar(1)
481                      varout_int(ib,1,1,1) = varin_int(imap(ib,jgrid),1,1,1)
482                   END DO
483                   iostat = nf90_put_var( ncid_out, jv, varout_int(1:lenvar(1),1:1,1,1) ) 
484                CASE( 3 )
485                   WRITE(*,*) 'This is a 3D integer.'
486                   ! Assume third dimension is time.
487                   iostat = nf90_get_var( ncid_in, jv, varin_int(1:lenvar(1),1:1,1,1:lenvar(3)) )
488                   DO itime = 1,lenvar(3)
489                      DO ib = 1,lenvar(1)
490                         varout_int(ib,1,1,itime) = varin_int(imap(ib,jgrid),1,1,itime)
491                      END DO
492                   END DO
493                   iostat = nf90_put_var( ncid_out, jv, varout_int(1:lenvar(1),1:1,1,1:lenvar(3)) ) 
494                CASE( 4 )
495                   WRITE(*,*) 'This is a 4D integer.'
496                   ! Assume third and fourth dimensions are depth and time respectively.
497                   iostat = nf90_get_var( ncid_in, jv, varin_int(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) )
498                   DO itime = 1,lenvar(4)
499                      DO idepth = 1,lenvar(3)
500                         DO ib = 1,lenvar(1)
501                            varout_int(ib,1,idepth,itime) = varin_int(imap(ib,jgrid),1,idepth,itime)
502                         END DO
503                      END DO
504                   END DO
505                   iostat = nf90_put_var( ncid_out, jv, varout_int(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) ) 
506                CASE DEFAULT
507                   WRITE(*,*) 'ERROR : Variable with ',ndims_var,' dimensions. Case not coded.'
508                   STOP
509             END SELECT
510
511          CASE( NF90_FLOAT ) 
512             SELECT CASE(ndims_var)
513                CASE( 1 )
514                   WRITE(*,*) 'This is a 1D float.'
515                   ! Assume this is a depth or time coordinate variable and copy across unchanged.
516                   iostat = nf90_get_var( ncid_in, jv, varin_float(1:lenvar(1),1,1,1) )
517                   iostat = nf90_put_var( ncid_out, jv, varout_float(1:lenvar(1),1,1,1) ) 
518                CASE( 2 )
519                   WRITE(*,*) 'This is a 2D float.'
520                   iostat = nf90_get_var( ncid_in, jv, varin_float(1:lenvar(1),1:1,1,1) )
521                   DO ib = 1,lenvar(1)
522                      varout_float(ib,1,1,1) = varin_float(imap(ib,jgrid),1,1,1)
523                   END DO
524                   iostat = nf90_put_var( ncid_out, jv, varout_float(1:lenvar(1),1:1,1,1) ) 
525                CASE( 3 )
526                   WRITE(*,*) 'This is a 3D float.'
527                   ! Assume third dimension is time.
528                   iostat = nf90_get_var( ncid_in, jv, varin_float(1:lenvar(1),1:1,1,1:lenvar(3)) )
529                   DO itime = 1,lenvar(3)
530                      DO ib = 1,lenvar(1)
531                         varout_float(ib,1,1,itime) = varin_float(imap(ib,jgrid),1,1,itime)
532                      END DO
533                   END DO
534                   iostat = nf90_put_var( ncid_out, jv, varout_float(1:lenvar(1),1:1,1,1:lenvar(3)) ) 
535                CASE( 4 )
536                   WRITE(*,*) 'This is a 4D float.'
537                   ! Assume third and fourth dimensions are depth and time respectively.
538                   iostat = nf90_get_var( ncid_in, jv, varin_float(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) )
539                   DO itime = 1,lenvar(4)
540                      DO idepth = 1,lenvar(3)
541                         DO ib = 1,lenvar(1)
542                            varout_float(ib,1,idepth,itime) = varin_float(imap(ib,jgrid),1,idepth,itime)
543                         END DO
544                      END DO
545                   END DO
546                   iostat = nf90_put_var( ncid_out, jv, varout_float(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) ) 
547                CASE DEFAULT
548                   WRITE(*,*) 'ERROR : Variable with ',ndims_var,' dimensions. Case not coded.'
549                   STOP
550             END SELECT
551
552          CASE( NF90_DOUBLE ) 
553             SELECT CASE(ndims_var)
554                CASE( 1 )
555                   WRITE(*,*) 'This is a 1D double.'
556                   ! Assume this is a depth or time coordinate variable and copy across unchanged.
557                   iostat = nf90_get_var( ncid_in, jv, varin_dble(1:lenvar(1),1,1,1) )
558                   iostat = nf90_put_var( ncid_out, jv, varout_dble(1:lenvar(1),1,1,1) ) 
559                CASE( 2 )
560                   WRITE(*,*) 'This is a 2D double.'
561                   iostat = nf90_get_var( ncid_in, jv, varin_dble(1:lenvar(1),1:1,1,1) )
562                   DO ib = 1,lenvar(1)
563                      varout_dble(ib,1,1,1) = varin_dble(imap(ib,jgrid),1,1,1)
564                   END DO
565                   iostat = nf90_put_var( ncid_out, jv, varout_dble(1:lenvar(1),1:1,1,1) ) 
566                CASE( 3 )
567                   WRITE(*,*) 'This is a 3D double.'
568                   ! Assume third dimension is time.
569                   iostat = nf90_get_var( ncid_in, jv, varin_dble(1:lenvar(1),1:1,1,1:lenvar(3)) )
570                   DO itime = 1,lenvar(3)
571                      DO ib = 1,lenvar(1)
572                         varout_dble(ib,1,1,itime) = varin_dble(imap(ib,jgrid),1,1,itime)
573                      END DO
574                   END DO
575                   iostat = nf90_put_var( ncid_out, jv, varout_dble(1:lenvar(1),1:1,1,1:lenvar(3)) ) 
576                CASE( 4 )
577                   WRITE(*,*) 'This is a 4D double.'
578                   ! Assume third and fourth dimensions are depth and time respectively.
579                   iostat = nf90_get_var( ncid_in, jv, varin_dble(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) )
580                   DO itime = 1,lenvar(4)
581                      DO idepth = 1,lenvar(3)
582                         DO ib = 1,lenvar(1)
583                            varout_dble(ib,1,idepth,itime) = varin_dble(imap(ib,jgrid),1,idepth,itime)
584                         END DO
585                      END DO
586                   END DO
587                   iostat = nf90_put_var( ncid_out, jv, varout_dble(1:lenvar(1),1:1,1:lenvar(3),1:lenvar(4)) ) 
588                CASE DEFAULT
589                   WRITE(*,*) 'ERROR : Variable with ',ndims_var,' dimensions. Case not coded.'
590                   STOP
591             END SELECT
592
593             CASE DEFAULT
594                WRITE(*,*) 'ERROR : Unrecognised data type.'           
595                STOP
596
597       END SELECT
598
599       END DO ! jv
600
601!-----------------------------------------------------------------------------
602! 4. Close input and output files
603!-----------------------------------------------------------------------------
604
605      iostat = nf90_close( ncid_out )
606      iostat = nf90_close( ncid_in )
607
608
609END PROGRAM bdy_reorder
Note: See TracBrowser for help on using the repository browser.