source: branches/nemo_v3_3_beta/NEMOGCM/TOOLS/WEIGHTS/src/scripshape.F90 @ 2352

Last change on this file since 2352 was 2352, checked in by sga, 10 years ago

NEMO branch nemo_v3_3_beta
Add NOCS tools based on SCRIP package for creating weights for interpolation on the fly
These now should build with the maketools script in the TOOLS directory using the same
architecture configuration file as the model (hopefully)

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
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.