New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

scripshape.F90 in utils/tools_dev_r12970_AGRIF_CMEMS/WEIGHTS/src – NEMO

Context Navigation

source: utils/tools_dev_r12970_AGRIF_CMEMS/WEIGHTS/src/scripshape.F90 @ 13025

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: