source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/nocsutil/scripshape.F90 @ 2448

Last change on this file since 2448 was 2448, checked in by smasson, 10 years ago

nemo_v3_3_beta: minor bugfix in WEIGHTS tools

File size: 13.3 KB
Line 
1      PROGRAM scripshape
2!
3! program to take output from the SCRIP weights generator
4! and rearrange the data into a series of 2D fields suitable
5! for reading with iom_get in NEMO configurations using the
6! interpolation on the fly option
7!
8      USE netcdf
9      IMPLICIT none
10      INTEGER    :: ncId, VarId, status
11      INTEGER    :: start(4), count(4)
12      CHARACTER(LEN=1) :: y
13      INTEGER    :: nd, ns, nl, nw, sx, sy, dx, dy
14      INTEGER    :: i, j, k, m, n, smax
15!
16!  ifort -O2 -o scripshape scripshape.f90 \
17!        -I/nerc/packages/netcdfifort/v3.6.0-pl1/include \
18!        -L/nerc/packages/netcdfifort/v3.6.0-pl1/lib -lnetcdf
19!
20!
21      INTEGER(KIND=4), ALLOCATABLE :: src(:)
22      INTEGER(KIND=4), ALLOCATABLE :: dst(:)
23      REAL(KIND=8), ALLOCATABLE :: wgt(:,:)
24      REAL(KIND=8), ALLOCATABLE :: src1(:,:),dst1(:,:),wgt1(:,:)
25      LOGICAL  :: around, verbose
26#if defined ARGC
27      INTEGER(KIND=4) :: iargc
28      EXTERNAL :: iargc
29#endif
30
31      CHARACTER(LEN=256) :: interp_file, output_file, name_file
32      INTEGER            :: ew_wrap
33      NAMELIST /shape_inputs/ interp_file, output_file, ew_wrap
34
35! scripshape requires 1 arguments; the name of the file containing
36! the input namelist.
37! This namelist contains:
38!     the name of the input file containing the weights ! produced by SCRIP in its format;
39!     the name of the new output file which ! is to contain the reorganized fields ready for input to NEMO.
40!     the east-west wrapping of the input grid (-1, 0, 1 and 2 are accepted values)
41!
42! E.g.
43!     interp_file = 'data_nemo_bilin.nc'
44!     output_file = 'weights_bilin.nc'
45!     ew_wrap     = 2
46!
47#if defined ARGC
48      IF (iargc() == 1) THEN
49        CALL getarg(1, name_file)
50      ELSE
51        WRITE(*,*) 'Usage: scripshape namelist_file'
52        STOP
53      ENDIF
54#else
55      WRITE(6,*) 'enter name of namelist file'
56      READ(5,*) name_file
57#endif
58      interp_file = 'none'
59      output_file = 'none'
60      ew_wrap = 0
61      OPEN(12, FILE=name_file, STATUS='OLD', FORM='FORMATTED')
62      READ(12, NML=shape_inputs)
63      CLOSE(12)
64!
65      INQUIRE(FILE = TRIM(interp_file), EXIST=around)
66      IF (.not.around) THEN
67       WRITE(*,*) 'Input file: '//TRIM(interp_file)//' not found'
68       STOP
69      ENDIF
70!
71      INQUIRE(FILE = TRIM(output_file), EXIST=around)
72      IF (around) THEN
73       WRITE(*,*) 'Output file: '//TRIM(output_file)//' exists'
74       WRITE(*,*) 'Ok to overwrite (y/n)?'
75       READ(5,'(a)') y
76       IF ( y .ne. 'y' .AND. y .ne. 'Y' ) STOP
77      ENDIF
78!
79      verbose = .true.
80!
81! Obtain grid size information from interp_file
82!
83      CALL ncgetsize
84!
85! Allocate array spaces
86!
87      ALLOCATE(src(nl), STAT=status)
88      IF(status /= 0 ) CALL alloc_err('src')
89      ALLOCATE(dst(nl), STAT=status)
90      IF(status /= 0 ) CALL alloc_err('dst')
91      ALLOCATE(wgt(nw,nl), STAT=status)
92      IF(status /= 0 ) CALL alloc_err('wgt')
93      ALLOCATE(src1(dx,dy), STAT=status)
94      IF(status /= 0 ) CALL alloc_err('src1')
95      ALLOCATE(dst1(dx,dy), STAT=status)
96      IF(status /= 0 ) CALL alloc_err('dst1')
97      ALLOCATE(wgt1(dx,dy), STAT=status)
98      IF(status /= 0 ) CALL alloc_err('wgt1')
99!
100! Read all required data from interp_file
101!
102      CALL ncgetfields
103!
104! Check that dst is monotonically increasing
105!
106      DO k = 1,nl-1
107       IF(dst(k+1).lt.dst(k)) THEN
108        WRITE(*,*) 'non-monotonic at ',k
109        WRITE(*,*) dst(k-4:k+16)
110        STOP
111       ENDIF
112      END DO
113!
114! Remove references to the top row of src
115!
116      IF(verbose) WRITE(*,*) &
117        'Removing references to the top row of the source grid'
118      smax = (sy-1)*sx
119      n = 0
120      DO k = 1,nl
121       IF(src(k).gt.smax-1) THEN
122        src(k) = src(k)-sx
123        n = n + 1
124       ENDIF
125      END DO
126      IF(verbose) WRITE(*,*) n,' values changed (',100.*n/nl,'%)'
127!
128! Loop through weights for each corner in turn and
129! rearrange weight fields into separate 2D fields for
130! reading with iom_get in NEMO
131!
132      DO k = 1,nw
133       DO n = 1,4
134
135        i = 0
136        j = 1
137        DO m = n,nl,4
138         i = i+1
139         IF(i.gt.dx) THEN
140          i = 1
141          j = j + 1
142         ENDIF
143         src1(i,j) = src(m)
144         dst1(i,j) = dst(m)
145         wgt1(i,j) = wgt(k,m)
146        END DO
147!
148! Write out this set which will be labelled with
149! a 2 digit number equal to n+4*(k-1)
150!
151        CALL wrwgts
152!   
153       END DO
154      END DO
155      STOP
156      CONTAINS
157!
158!----------------------------------------------------------------------*
159      SUBROUTINE ncgetsize
160!
161! Access grid size information in interp_file and set the
162! following integers:
163!
164!    nd = dst_grid_size
165!    ns = src_grid_size
166!    nl = num_links   
167!    nw = num_wgts     
168! sx,sy = src_grid_dims     
169! dx,dy = dst_grid_dims     
170!
171      INTEGER idims(2)
172!
173       status = nf90_open(interp_file, nf90_NoWrite, ncid)
174       IF(status /= nf90_NoErr) CALL handle_err(status)
175!
176       status = nf90_inq_dimid(ncid, 'dst_grid_size', VarId)
177       IF(status /= nf90_NoErr) CALL handle_err(status)
178       status = nf90_inquire_dimension(ncid, VarId, LEN = nd)
179       IF(status /= nf90_NoErr) CALL handle_err(status)
180!
181       status = nf90_inq_dimid(ncid, 'src_grid_size', VarId)
182       IF(status /= nf90_NoErr) CALL handle_err(status)
183       status = nf90_inquire_dimension(ncid, VarId, LEN = ns)
184       IF(status /= nf90_NoErr) CALL handle_err(status)
185!
186       status = nf90_inq_dimid(ncid, 'num_links', VarId)
187       IF(status /= nf90_NoErr) CALL handle_err(status)
188       status = nf90_inquire_dimension(ncid, VarId, LEN = nl)
189       IF(status /= nf90_NoErr) CALL handle_err(status)
190!
191       status = nf90_inq_dimid(ncid, 'num_wgts', VarId)
192       IF(status /= nf90_NoErr) CALL handle_err(status)
193       status = nf90_inquire_dimension(ncid, VarId, LEN = nw)
194       IF(status /= nf90_NoErr) CALL handle_err(status)
195!
196       start = 1
197       count = 2
198       status = nf90_inq_varid(ncid, 'src_grid_dims', VarId)
199       IF(status /= nf90_NoErr) CALL handle_err(status)
200       status = nf90_get_var(ncid, VarId, idims, start, count)
201       IF(status /= nf90_NoErr) CALL handle_err(status)
202       sx = idims(1) ; sy = idims(2)
203!
204       status = nf90_inq_varid(ncid, 'dst_grid_dims', VarId)
205       IF(status /= nf90_NoErr) CALL handle_err(status)
206       status = nf90_get_var(ncid, VarId, idims, start, count)
207       IF(status /= nf90_NoErr) CALL handle_err(status)
208       dx = idims(1) ; dy = idims(2)
209!
210       status = nf90_close(ncid)
211       IF (status /= nf90_noerr) CALL handle_err(status)
212!
213       IF(verbose) THEN
214         WRITE(*,*) 'Detected sizes: '
215         WRITE(*,*) 'dst_grid_size: ', nd
216         WRITE(*,*) 'src_grid_size: ', ns
217         WRITE(*,*) 'num_links    : ', nl
218         WRITE(*,*) 'num_wgts     : ', nw
219         WRITE(*,*) 'src_grid_dims: ', sx, ' x ', sy
220         WRITE(*,*) 'dst_grid_dims: ', dx, ' x ', dy
221       ENDIF
222!
223      END SUBROUTINE ncgetsize
224
225!----------------------------------------------------------------------*
226      SUBROUTINE ncgetfields
227!
228! Read all required data from interp_file. The data read are:
229!
230! netcdf variable    size   internal array
231!-----------------+-------+--------------
232! src_address        nl     src
233! dst_address        nl     dst
234! remap_matrix     (nw,nl)  wgt
235!
236       status = nf90_open(interp_file, nf90_NoWrite, ncid)
237       IF(status /= nf90_NoErr) CALL handle_err(status)
238!
239       status = nf90_inq_varid(ncid, 'src_address', VarId)
240       IF(status /= nf90_NoErr) CALL handle_err(status)
241!
242! Read the values for src
243       status = nf90_get_var(ncid, VarId, src, &
244                         start = (/ 1 /),      &
245                         count = (/ nl /))
246       IF(status /= nf90_NoErr) CALL handle_err(status)
247!
248       status = nf90_inq_varid(ncid, 'dst_address', VarId)
249       IF(status /= nf90_NoErr) CALL handle_err(status)
250!
251! Read the values for dst
252       status = nf90_get_var(ncid, VarId, dst, &
253                         start = (/ 1 /),      &
254                         count = (/ nl /))
255       IF(status /= nf90_NoErr) CALL handle_err(status)
256!
257       status = nf90_inq_varid(ncid, 'remap_matrix', VarId)
258       IF(status /= nf90_NoErr) CALL handle_err(status)
259!
260! Read the values for wgt
261       status = nf90_get_var(ncid, VarId, wgt, &
262                         start = (/ 1, 1 /),   &
263                         count = (/ nw, nl /))
264       IF(status /= nf90_NoErr) CALL handle_err(status)
265!
266       status = nf90_close(ncid)
267       IF (status /= nf90_noerr) CALL handle_err(status)
268!
269      END SUBROUTINE ncgetfields
270
271!----------------------------------------------------------------------*
272      SUBROUTINE handle_err(status)
273!
274! Simple netcdf error checking routine
275!
276      INTEGER, intent ( in) :: status
277!
278      IF(status /= nf90_noerr) THEN
279       IF(trim(nf90_strerror(status)) .eq. 'Attribute not found') THEN
280! ignore
281       ELSE
282        WRITE(*,*) trim(nf90_strerror(status))
283        STOP "Stopped"
284       END IF
285      END IF
286      END SUBROUTINE handle_err
287
288!----------------------------------------------------------------------*
289      SUBROUTINE alloc_err(arname)
290!
291! Simple allocation error checking routine
292!
293      CHARACTER(LEN=*) :: arname
294!
295      WRITE(*,*) 'Allocation error attempting to ALLOCATE '//arname
296      STOP "Stopped"
297      END SUBROUTINE alloc_err
298
299!                                                                       
300!----------------------------------------------------------------------*
301      SUBROUTINE wrwgts
302!
303! Write out each set of 2D fields to output_file.
304! Each call will write out a set of srcXX, dstXX and wgtXX fields
305! where XX is a two digit number equal to n + 4*(k-1). The first
306! and last calls to this routine initialise and close the output
307! file respectively. The first call is detected when k*n=1 and the
308! last call is detected when k*n=4*nw. The outfile file remains
309! open between the first and last calls.
310!
311      INTEGER :: status, ncid, ncin
312      INTEGER :: Lontdid, Lattdid
313      INTEGER :: tvid, tvid2, tvid3
314      INTEGER :: ioldfill
315      CHARACTER(LEN=2) :: cs
316      SAVE ncid, Lontdid, Lattdid
317!
318      IF(k*n.eq.1) THEN
319!
320! Create output_file and set the dimensions
321!
322         status = nf90_create(output_file, nf90_Clobber, ncid)
323         IF(status /= nf90_NoErr) CALL handle_err(status)
324         status = nf90_set_fill(ncid, nf90_NoFill, ioldfill)
325         IF(status /= nf90_NoErr) CALL handle_err(status)
326!
327         status = nf90_def_dim(ncid, "lon", dx, Lontdid)
328         IF(status /= nf90_NoErr) CALL handle_err(status)
329         status = nf90_def_dim(ncid, "lat", dy, Lattdid)
330         IF(status /= nf90_NoErr) CALL handle_err(status)
331!
332         status = nf90_put_att(ncid, nf90_global, 'ew_wrap', ew_wrap)
333         IF(status /= nf90_NoErr) CALL handle_err(status)
334      ELSE
335!
336! Reenter define mode
337!
338         status = nf90_redef(ncid)
339         IF(status /= nf90_NoErr) CALL handle_err(status)
340      ENDIF
341!
342      WRITE(cs,'(i2.2)') n + 4*(k-1)
343!
344! Define new variables
345!
346      status = nf90_def_var(ncid, "src"//cs, nf90_double, &
347                         (/ Lontdid, Lattdid /), tvid)
348      IF(status /= nf90_NoErr) CALL handle_err(status)
349!
350      status = nf90_def_var(ncid, "dst"//cs, nf90_double, &
351                         (/ Lontdid, Lattdid /), tvid2)
352      IF(status /= nf90_NoErr) CALL handle_err(status)
353!
354      status = nf90_def_var(ncid, "wgt"//cs, nf90_double, &
355                         (/ Lontdid, Lattdid /), tvid3)
356      IF(status /= nf90_NoErr) CALL handle_err(status)
357!
358! Leave define mode
359!
360      status = nf90_enddef(ncid)
361      IF(status /= nf90_NoErr) CALL handle_err(status)
362!
363! Write the data
364!
365      status = nf90_put_var(ncid, tvid, src1,    &
366                            start = (/ 1, 1 /),  &
367                            count = (/ dx, dy /) )
368      IF(status /= nf90_NoErr) CALL handle_err(status)
369!
370      status = nf90_put_var(ncid, tvid2, dst1,   &
371                            start = (/ 1, 1 /),  &
372                            count = (/ dx, dy /) )
373      IF(status /= nf90_NoErr) CALL handle_err(status)
374!
375      status = nf90_put_var(ncid, tvid3, wgt1,   &
376                            start = (/ 1, 1  /), &
377                            count = (/ dx, dy /) )
378      IF(status /= nf90_NoErr) CALL handle_err(status)
379!
380      IF(k*n.eq.4*nw) THEN
381!
382! --     Reenter define mode
383!
384         status = nf90_redef(ncid)
385         IF(status /= nf90_NoErr) CALL handle_err(status)
386!
387! --     Reopen interp_file and transfer some global attributes
388!
389         status = nf90_open(interp_file, nf90_NoWrite, ncin)
390         IF(status /= nf90_NoErr) CALL handle_err(status)
391!
392         status = nf90_copy_att(ncin,NF90_GLOBAL,'title',        ncid,NF90_GLOBAL)
393         IF(status /= nf90_NoErr) CALL handle_err(status)
394!
395         status = nf90_copy_att(ncin,NF90_GLOBAL,'normalization',ncid,NF90_GLOBAL)
396         IF(status /= nf90_NoErr) CALL handle_err(status)
397!
398         status = nf90_copy_att(ncin,NF90_GLOBAL,'map_method',   ncid,NF90_GLOBAL)
399         IF(status /= nf90_NoErr) CALL handle_err(status)
400!
401         status = nf90_copy_att(ncin,NF90_GLOBAL,'conventions',  ncid,NF90_GLOBAL)
402         IF(status /= nf90_NoErr) CALL handle_err(status)
403!
404         status = nf90_copy_att(ncin,NF90_GLOBAL,'history',      ncid,NF90_GLOBAL)
405         IF(status /= nf90_NoErr) CALL handle_err(status)
406!
407! --     Close interp_file
408!
409         status = nf90_close(ncin)
410         IF(status /= nf90_NoErr) CALL handle_err(status)
411!
412! --     Close output_file
413!
414         status = nf90_close(ncid)
415         IF(status /= nf90_NoErr) CALL handle_err(status)
416      ENDIF
417
418      END SUBROUTINE wrwgts
419      END PROGRAM scripshape
Note: See TracBrowser for help on using the repository browser.