Context Navigation

interp.f90 @ 13369

Last change on this file since 13369 was 13369, checked in by jpaul, 4 years ago
update: cf changelog inside documentation
File size: 34.2 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! DESCRIPTION:
6	!> @brief
7	!> This module manage interpolation on regular grid.
8	!>
9	!> @details Interpolation method to be used is specify inside variable
10	!> strcuture, as array of string character.<br/>
11	!> - td_var\%c_interp(1) string character is the interpolation name choose between:
12	!> - 'nearest'
13	!> - 'cubic '
14	!> - 'linear '
15	!> - td_var\%c_interp(2) string character is an operation to be used
16	!> on interpolated value.<br/>
17	!> operation have to be mulitplication '*' or division '/'.<br/>
18	!> coefficient have to be refinement factor following i-direction 'rhoi',
19	!> j-direction 'rhoj', or k-direction 'rhok'.<br/>
20	!>
21	!> Examples: '*rhoi', '/rhoj'.
22	!>
23	!> @note Those informations are read from namelist or variable configuration file (default).<br/>
24	!> Interplation method could be specify for each variable in namelist _namvar_,
25	!> defining string character _cn\_varinfo_.<br/>
26	!> Example:
27	!> - cn_varinfo='varname1:int=cubic/rhoi', 'varname2:int=linear'
28	!>
29	!> to create mixed grid (with coarse grid point needed to compute
30	!> interpolation):<br/>
31	!> @code
32	!> CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] )
33	!> @endcode
34	!> - td_var is coarse grid variable (should be extrapolated)
35	!> - td_mix is mixed grid variable structure [output]
36	!> - id_rho is array of refinment factor [optional]
37	!>
38	!> to detected point to be interpolated:<br/>
39	!> @code
40	!> il_detect(:,:,:)=interp_detect( td_mix [,id_rho] )
41	!> @endcode
42	!> - il_detect(:,:,:) is 3D array of detected point to be interpolated
43	!> - td_mix is mixed grid variable
44	!> - id_rho is array of refinement factor [optional]
45	!>
46	!> to interpolate variable value:<br/>
47	!> @code
48	!> CALL interp_fill_value( td_var [,id_rho] [,id_offset] )
49	!> @endcode
50	!> - td_var is variable structure
51	!> - id_rho is array of refinement factor [optional]
52	!> - id_offset is array of offset between fine and coarse grid [optional]
53	!>
54	!> to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/>
55	!> @code
56	!> CALL interp_clean_mixed_grid( td_mix, td_var, id_rho )
57	!> @endcode
58	!> - td_mix is mixed grid variable structure
59	!> - td_var is variable structure [output]
60	!> - id_rho is array of refinement factor [optional]
61	!> - id_offset is array of offset between fine and coarse grid [optional]
62	!>
63	!> @note It use to work on ORCA grid, as we work only with grid indices.
64	!>
65	!> @warning due to the use of second derivative when using cubic interpolation
66	!> you should add at least 2 extrabands.
67	!>
68	!> @author
69	!> J.Paul
70	!>
71	!> @date November, 2013 - Initial Version
72	!> @date September, 2014
73	!> - add header
74	!> - use interpolation method modules
75	!>
76	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
77	!----------------------------------------------------------------------
78	MODULE interp
79
80	USE netcdf ! nf90 library
81	USE global ! global variable
82	USE kind ! F90 kind parameter
83	USE logger ! log file manager
84	USE fct ! basic useful function
85	USE date ! date manager
86	USE att ! attribute manager
87	USE dim ! dimension manager
88	USE var ! variable manager
89	USE grid ! grid manager
90	USE extrap ! extrapolation manager
91	USE interp_cubic ! cubic interpolation manager
92	USE interp_linear ! linear interpolation manager
93	USE interp_nearest ! nearest interpolation manager
94
95	IMPLICIT NONE
96	! NOTE_avoid_public_variables_if_possible
97
98	! type and variable
99
100	! function and subroutine
101	PUBLIC :: interp_detect !< detected point to be interpolated
102	PUBLIC :: interp_fill_value !< interpolate value
103	PUBLIC :: interp_create_mixed_grid !< create mixed grid
104	PUBLIC :: interp_clean_mixed_grid !< clean mixed grid
105
106	PRIVATE :: interp__detect ! detected point to be interpolated
107	PRIVATE :: interp__detect_wrapper ! detected point to be interpolated
108	PRIVATE :: interp__fill_value_wrapper ! interpolate value over detectected point
109	PRIVATE :: interp__fill_value ! interpolate value over detectected point
110	PRIVATE :: interp__clean_even_grid ! clean even mixed grid
111	PRIVATE :: interp__check_method ! check if interpolation method available
112
113	TYPE TINTERP
114	CHARACTER(LEN=lc) :: c_name = '' !< interpolation method name
115	CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor
116	CHARACTER(LEN=lc) :: c_divisor= '' !< interpolation divisor
117	END TYPE TINTERP
118
119	INTERFACE interp_detect
120	MODULE PROCEDURE interp__detect_wrapper
121	END INTERFACE interp_detect
122
123	INTERFACE interp_fill_value
124	MODULE PROCEDURE interp__fill_value_wrapper
125	END INTERFACE interp_fill_value
126
127	CONTAINS
128	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129	FUNCTION interp__check_method(cd_method) &
130	& RESULT (lf_avail)
131	!-------------------------------------------------------------------
132	!> @brief
133	!> This function check if interpolation method is available.
134	!>
135	!> @details
136	!> check if name of interpolation method is present in global list of string
137	!> character cp_interp_list (see global.f90).
138	!>
139	!> @author J.Paul
140	!> @date November, 2013 - Initial Version
141	!
142	!> @param[in] cd_method interpolation method
143	!> @return true if interpolation method is available
144	!-------------------------------------------------------------------
145
146	IMPLICIT NONE
147
148	! Argument
149	CHARACTER(LEN=lc) :: cd_method
150
151	! function
152	LOGICAL :: lf_avail
153
154	! local variable
155	CHARACTER(LEN=lc) :: cl_interp
156	CHARACTER(LEN=lc) :: cl_method
157
158	! loop indices
159	INTEGER(I4) :: ji
160	!----------------------------------------------------------------
161
162	cl_method=fct_lower(cd_method)
163
164	lf_avail=.FALSE.
165	DO ji=1,ip_ninterp
166	cl_interp=fct_lower(cp_interp_list(ji))
167	IF( TRIM(cl_interp) == TRIM(cl_method) )THEN
168	lf_avail=.TRUE.
169	EXIT
170	ENDIF
171	ENDDO
172
173	END FUNCTION interp__check_method
174	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175	FUNCTION interp__detect_wrapper(td_mix, id_rho) &
176	& RESULT (if_detect)
177	!-------------------------------------------------------------------
178	!> @brief
179	!> This function detected point to be interpolated.
180	!>
181	!> @details
182	!> Actually it checks, the number of dimension used for this variable
183	!> and launch interp__detect which detected point to be interpolated.
184	!>
185	!> @author J.Paul
186	!> @date November, 2013 - Initial Version
187	!>
188	!> @param[in] td_mix mixed grid variable (to interpolate)
189	!> @param[in] id_rho array of refinement factor
190	!> @return 3D array of detected point to be interpolated
191	!-------------------------------------------------------------------
192
193	IMPLICIT NONE
194
195	! Argument
196	TYPE(TVAR) , INTENT(IN) :: td_mix
197	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
198
199	! function
200	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
201	& td_mix%t_dim(2)%i_len,&
202	& td_mix%t_dim(3)%i_len ) :: if_detect
203
204	! local variable
205	! loop indices
206	!----------------------------------------------------------------
207
208	IF( .NOT. ANY(td_mix%t_dim(1:3)%l_use) )THEN
209	! no dimension I-J-K used
210	CALL logger_debug(" INTERP DETECT: nothing done for variable"//&
211	& TRIM(td_mix%c_name) )
212
213	if_detect(:,:,:)=0
214
215	ELSE IF( ALL(td_mix%t_dim(1:3)%l_use) )THEN
216
217	! detect point to be interpolated on I-J-K
218	CALL logger_debug(" INTERP DETECT: detect point "//&
219	& TRIM(td_mix%c_point)//" for variable "//&
220	& TRIM(td_mix%c_name) )
221
222	if_detect(:,:,:)=interp__detect( td_mix, id_rho(:) )
223
224	ELSE IF( ALL(td_mix%t_dim(1:2)%l_use) )THEN
225
226	! detect point to be interpolated on I-J
227	CALL logger_debug(" INTERP DETECT: detect point "//&
228	& TRIM(td_mix%c_point)//" for variable "//&
229	& TRIM(td_mix%c_name) )
230
231	if_detect(:,:,1:1)=interp__detect( td_mix, id_rho(:))
232
233	ELSE IF( td_mix%t_dim(3)%l_use )THEN
234
235	! detect point to be interpolated on K
236	CALL logger_debug(" INTERP DETECT: detect vertical point "//&
237	& " for variable "//TRIM(td_mix%c_name) )
238
239	if_detect(1:1,1:1,:)=interp__detect( td_mix, id_rho(:) )
240
241	ENDIF
242
243	END FUNCTION interp__detect_wrapper
244	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
245	FUNCTION interp__detect(td_mix, id_rho) &
246	& RESULT (if_detect)
247	!-------------------------------------------------------------------
248	!> @brief
249	!> This function detected point to be interpolated.
250	!>
251	!> @details
252	!> A special case is done for even refinement on ARAKAWA-C grid.
253	!>
254	!> @author J.Paul
255	!> @date November, 2013 - Initial Version
256	!
257	!> @param[in] td_mix mixed grid variable (to interpolate)
258	!> @param[in] id_rho array of refinement factor
259	!> @return 3D array of detected point to be interpolated
260	!-------------------------------------------------------------------
261
262	IMPLICIT NONE
263
264	! Argument
265	TYPE(TVAR) , INTENT(IN) :: td_mix
266	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
267
268	! function
269	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
270	& td_mix%t_dim(2)%i_len,&
271	& td_mix%t_dim(3)%i_len ) :: if_detect
272
273	! local variable
274	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
275
276	INTEGER(I4) :: il_xextra
277	INTEGER(I4) :: il_yextra
278	INTEGER(I4) :: il_zextra
279
280	INTEGER(i4), DIMENSION(3) :: il_dim
281
282	LOGICAL , DIMENSION(3) :: ll_even
283
284	! loop indices
285	INTEGER(I4) :: ji
286	INTEGER(I4) :: jj
287	INTEGER(I4) :: jk
288	!----------------------------------------------------------------
289	ALLOCATE( il_rho(ip_maxdim) )
290	il_rho(:)=1
291	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
292
293	! special case for even refinement on ARAKAWA-C grid
294	ll_even(:)=.FALSE.
295	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
296	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
297	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
298
299	SELECT CASE(TRIM(td_mix%c_point))
300	CASE('U')
301	ll_even(1)=.FALSE.
302	CASE('V')
303	ll_even(2)=.FALSE.
304	CASE('F')
305	ll_even(:)=.FALSE.
306	END SELECT
307
308	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
309	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
310	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
311
312	! special case for cubic interpolation
313	il_xextra=0
314	il_yextra=0
315	il_zextra=0
316	SELECT CASE(TRIM(td_mix%c_interp(1)))
317	CASE('cubic')
318	! those points can not be compute cause cubic interpolation
319	! need second derivative.
320	IF( il_rho(jp_I) /= 1 ) il_xextra=3*il_rho(jp_I)
321	IF( il_rho(jp_J) /= 1 ) il_yextra=3*il_rho(jp_J)
322	IF( il_rho(jp_K) /= 1 ) il_zextra=3*il_rho(jp_K)
323	END SELECT
324
325	il_dim(:)=td_mix%t_dim(1:3)%i_len
326
327	! init
328	if_detect(:,:,:)=1
329
330	! do not compute coarse grid point
331	if_detect(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
332	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
333	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K) ) = 0
334
335	! do not compute point near fill value
336	DO jk=1,il_dim(3),il_rho(jp_K)
337	DO jj=1,il_dim(2),il_rho(jp_J)
338	DO ji=1,il_dim(1),il_rho(jp_I)
339
340	IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN
341
342	! i-direction
343	if_detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
344	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
345	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
346	! j-direction
347	if_detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
348	& MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
349	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
350	! k-direction
351	if_detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
352	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
353	& MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0
354
355	ENDIF
356
357	ENDDO
358	ENDDO
359	ENDDO
360
361	DEALLOCATE( il_rho )
362
363	END FUNCTION interp__detect
364	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
365	SUBROUTINE interp_create_mixed_grid(td_var, td_mix, id_rho)
366	!-------------------------------------------------------------------
367	!> @brief
368	!> This subroutine create mixed grid.
369	!>
370	!> @details
371	!> Created grid is fine resolution grid.
372	!> First and last point are coasre grid point.
373	!>
374	!> A special case is done for even refinement on ARAKAWA-C grid.
375	!>
376	!> @author J.Paul
377	!> @date November, 2013 - Initial Version
378	!>
379	!> @param[in] td_var coarse grid variable (should be extrapolated)
380	!> @param[out] td_mix mixed grid variable
381	!> @param[in] id_rho array of refinment factor (default 1)
382	!-------------------------------------------------------------------
383
384	IMPLICIT NONE
385
386	! Argument
387	TYPE(TVAR) , INTENT(IN ) :: td_var
388	TYPE(TVAR) , INTENT( OUT) :: td_mix
389	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
390
391	! local variable
392	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
393
394	INTEGER(i4) :: il_xextra
395	INTEGER(i4) :: il_yextra
396	INTEGER(i4) :: il_zextra
397
398	LOGICAL, DIMENSION(3) :: ll_even
399
400	! loop indices
401	!----------------------------------------------------------------
402	ALLOCATE(il_rho(ip_maxdim))
403	il_rho(:)=1
404	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
405
406	! special case for even refinement on ARAKAWA-C grid
407	ll_even(:)=.FALSE.
408	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
409	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
410	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
411
412	SELECT CASE(TRIM(td_var%c_point))
413	CASE('U')
414	ll_even(1)=.FALSE.
415	CASE('V')
416	ll_even(2)=.FALSE.
417	CASE('F')
418	ll_even(:)=.FALSE.
419	END SELECT
420
421	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
422	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
423	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
424
425	! copy variable
426	td_mix=var_copy(td_var)
427
428	! compute new dimension length
429	il_xextra=il_rho(jp_I)-1
430	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
431
432	il_yextra=il_rho(jp_J)-1
433	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
434
435	il_zextra=il_rho(jp_K)-1
436	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
437
438	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
439	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
440	& td_mix%t_dim(2)%i_len, &
441	& td_mix%t_dim(3)%i_len, &
442	& td_mix%t_dim(4)%i_len) )
443
444	! initialise to FillValue
445	td_mix%d_value(:,:,:,:)=td_mix%d_fill
446
447	! quid qd coord ou bathy fourni par user ?? (offset ??)
448	td_mix%d_value(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
449	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
450	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K), :) = &
451	& td_var%d_value(:,:,:,:)
452
453	DEALLOCATE(il_rho)
454
455	END SUBROUTINE interp_create_mixed_grid
456	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
457	SUBROUTINE interp__clean_even_grid(td_mix, id_rho)
458	!-------------------------------------------------------------------
459	!> @brief
460	!> This subroutine remove points added to mixed grid to compute
461	!> interpolation in the special case of even refinement on ARAKAWA-C grid.
462	!>
463	!> @details
464	!>
465	!> @author J.Paul
466	!> @date November, 2013 - Initial Version
467	!>
468	!> @param[inout] td_mix mixed grid variable
469	!> @param[in] id_rho array of refinment factor
470	!-------------------------------------------------------------------
471
472	IMPLICIT NONE
473
474	! Argument
475	TYPE(TVAR) , INTENT(INOUT) :: td_mix
476	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
477
478	! local variable
479	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
480
481	INTEGER(i4) :: il_xextra
482	INTEGER(i4) :: il_yextra
483	INTEGER(i4) :: il_zextra
484
485	LOGICAL, DIMENSION(3) :: ll_even
486
487	LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ll_mask
488
489	REAL(dp), DIMENSION(:) , ALLOCATABLE :: dl_vect
490
491	TYPE(TVAR) :: tl_mix
492
493	! loop indices
494	!----------------------------------------------------------------
495	ALLOCATE(il_rho(ip_maxdim))
496	il_rho(:)=1
497	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
498
499	! special case for even refinement on ARAKAWA-C grid
500	ll_even(:)=.FALSE.
501	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
502	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
503	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
504
505	SELECT CASE(TRIM(td_mix%c_point))
506	CASE('U')
507	ll_even(1)=.FALSE.
508	CASE('V')
509	ll_even(2)=.FALSE.
510	CASE('F')
511	ll_even(:)=.FALSE.
512	END SELECT
513
514	! remove some point only if refinement in some direction is even
515	IF( ANY(ll_even(:)) )THEN
516
517	! copy variable
518	tl_mix=var_copy(td_mix)
519
520	ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, &
521	& tl_mix%t_dim(2)%i_len, &
522	& tl_mix%t_dim(3)%i_len, &
523	& tl_mix%t_dim(4)%i_len) )
524
525	ll_mask(:,:,:,:)=.TRUE.
526
527	IF( tl_mix%t_dim(1)%l_use .AND. ll_even(1) )THEN
528
529	il_rho(jp_I)=il_rho(jp_I)+1
530
531	! locate wrong point on mixed grid
532	ll_mask(1:td_mix%t_dim(1)%i_len:il_rho(jp_I),:,:,:)=.FALSE.
533
534	! compute coasre grid dimension length
535	il_xextra=il_rho(jp_I)-1
536	td_mix%t_dim(1)%i_len=(tl_mix%t_dim(1)%i_len+il_xextra)/il_rho(jp_I)
537
538	il_rho(jp_I)=il_rho(jp_I)-1
539	! compute right fine grid dimension length
540	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
541
542	ENDIF
543
544	IF( tl_mix%t_dim(2)%l_use .AND. ll_even(2) )THEN
545
546	il_rho(jp_J)=il_rho(jp_J)+1
547
548	! locate wrong point on mixed grid
549	ll_mask(:,1:tl_mix%t_dim(2)%i_len:il_rho(jp_J),:,:)=.FALSE.
550
551	! compute coasre grid dimension length
552	il_yextra=il_rho(jp_J)-1
553	td_mix%t_dim(2)%i_len=(tl_mix%t_dim(2)%i_len+il_yextra)/il_rho(jp_J)
554
555	il_rho(jp_J)=il_rho(jp_J)-1
556	! compute right fine grid dimension length
557	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
558
559	ENDIF
560
561	IF( tl_mix%t_dim(3)%l_use .AND. ll_even(3) )THEN
562
563	il_rho(jp_K)=il_rho(jp_K)+1
564
565	! locate wrong point on mixed grid
566	ll_mask(:,:,1:tl_mix%t_dim(3)%i_len:il_rho(jp_K),:)=.FALSE.
567
568	! compute coasre grid dimension length
569	il_zextra=il_rho(jp_K)-1
570	td_mix%t_dim(3)%i_len=(tl_mix%t_dim(3)%i_len+il_zextra)/il_rho(jp_K)
571
572	il_rho(jp_K)=il_rho(jp_K)-1
573	! compute right fine grid dimension length
574	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
575
576	ENDIF
577
578	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
579	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
580	& td_mix%t_dim(2)%i_len, &
581	& td_mix%t_dim(3)%i_len, &
582	& td_mix%t_dim(4)%i_len) )
583
584	! initialise to FillValue
585	td_mix%d_value(:,:,:,:)=td_mix%d_fill
586
587	IF( COUNT(ll_mask(:,:,:,:)) /= SIZE(td_mix%d_value(:,:,:,:)) )THEN
588
589	CALL logger_error("INTERP CLEAN EVEN GRID: output value size "//&
590	& " and mask count differ ")
591
592	ELSE
593
594	ALLOCATE( dl_vect(COUNT(ll_mask(:,:,:,:))) )
595
596	dl_vect(:)= PACK( tl_mix%d_value(:,:,:,:), &
597	& MASK=ll_mask(:,:,:,:) )
598
599	td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), &
600	& SHAPE=(/td_mix%t_dim(1)%i_len, &
601	& td_mix%t_dim(2)%i_len, &
602	& td_mix%t_dim(3)%i_len, &
603	& td_mix%t_dim(4)%i_len/) )
604
605	DEALLOCATE( dl_vect )
606
607	ENDIF
608
609	DEALLOCATE( ll_mask )
610
611	! clean
612	CALL var_clean(tl_mix)
613
614	ENDIF
615
616	! clean
617	DEALLOCATE(il_rho)
618
619	END SUBROUTINE interp__clean_even_grid
620	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
621	SUBROUTINE interp_clean_mixed_grid(td_mix, td_var, &
622	& id_rho, id_offset)
623	!-------------------------------------------------------------------
624	!> @brief
625	!> This subroutine remove points added on mixed grid
626	!> to compute interpolation. And save interpolated value over domain.
627	!>
628	!> @details
629	!>
630	!> @author J.Paul
631	!> @date November, 2013 - Initial Version
632	!> @date September, 2014
633	!> - use offset to save useful domain
634	!>
635	!> @param[in] td_mix mixed grid variable structure
636	!> @param[out] td_var variable structure
637	!> @param[in] id_rho array of refinement factor (default 1)
638	!> @param[in] id_offset 2D array of offset between fine and coarse grid
639	!-------------------------------------------------------------------
640
641	IMPLICIT NONE
642
643	! Argument
644	TYPE(TVAR) , INTENT(IN ) :: td_mix
645	TYPE(TVAR) , INTENT( OUT) :: td_var
646	INTEGER(I4), DIMENSION(:) , INTENT(IN ) :: id_rho
647	INTEGER(I4), DIMENSION(2,2), INTENT(IN ) :: id_offset
648
649	! local variable
650	INTEGER(i4) :: il_imin0
651	INTEGER(i4) :: il_imax0
652	INTEGER(i4) :: il_jmin0
653	INTEGER(i4) :: il_jmax0
654
655	INTEGER(i4) :: il_imin1
656	INTEGER(i4) :: il_jmin1
657	INTEGER(i4) :: il_imax1
658	INTEGER(i4) :: il_jmax1
659
660	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
661
662	TYPE(TVAR) :: tl_mix
663
664	! loop indices
665	!----------------------------------------------------------------
666	! copy mixed variable in temporary structure
667	tl_mix=var_copy(td_mix)
668
669	! remove useless points over mixed grid for even refinement
670	CALL interp__clean_even_grid(tl_mix, id_rho(:))
671
672	! copy cleaned mixed variable
673	td_var=var_copy(tl_mix)
674
675	! delete array of value
676	CALL var_del_value(td_var)
677
678	! compute domain indices in i-direction
679	il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len
680
681	IF( td_var%t_dim(1)%l_use )THEN
682	il_imin1=il_imin0+id_offset(Jp_I,1)
683	il_imax1=il_imax0-id_offset(Jp_I,2)
684	ELSE
685
686	il_imin1=il_imin0
687	il_imax1=il_imax0
688
689	ENDIF
690
691	! compute domain indices in j-direction
692	il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len
693
694	IF( td_var%t_dim(2)%l_use )THEN
695	il_jmin1=il_jmin0+id_offset(Jp_J,1)
696	il_jmax1=il_jmax0-id_offset(Jp_J,2)
697	ELSE
698
699	il_jmin1=il_jmin0
700	il_jmax1=il_jmax0
701
702	ENDIF
703
704	! compute new dimension
705	td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
706	td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
707
708	ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
709	& td_var%t_dim(2)%i_len, &
710	& td_var%t_dim(3)%i_len, &
711	& td_var%t_dim(4)%i_len) )
712
713	dl_value( 1:td_var%t_dim(1)%i_len, &
714	& 1:td_var%t_dim(2)%i_len, &
715	& :,:) = tl_mix%d_value( il_imin1:il_imax1, &
716	& il_jmin1:il_jmax1, &
717	& :, : )
718
719	! add variable value
720	CALL var_add_value(td_var,dl_value(:,:,:,:))
721
722	DEALLOCATE(dl_value)
723
724	! clean
725	CALL var_clean(tl_mix)
726
727	END SUBROUTINE interp_clean_mixed_grid
728	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
729	SUBROUTINE interp__fill_value_wrapper(td_var, &
730	& id_rho, &
731	& id_offset)
732	!-------------------------------------------------------------------
733	!> @brief
734	!> This subroutine interpolate variable value.
735	!>
736	!> @details
737	!> Actually it checks, the number of dimension used for this variable
738	!> and launch interp__fill_value.
739	!>
740	!> @author J.Paul
741	!> @date November, 2013 - Initial Version
742	!>
743	!> @param[inout] td_var variable structure
744	!> @param[in] id_rho array of refinement factor
745	!> @param[in] id_offset 2D array of offset between fine and coarse grid
746	!-------------------------------------------------------------------
747
748	IMPLICIT NONE
749
750	! Argument
751	TYPE(TVAR) , INTENT(INOUT) :: td_var
752	INTEGER(I4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
753	INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
754
755	! local variable
756	CHARACTER(LEN=lc) :: cl_method
757	INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho
758	INTEGER(i4) , DIMENSION(2,2) :: il_offset
759
760	! loop indices
761	!----------------------------------------------------------------
762
763	ALLOCATE( il_rho(ip_maxdim) )
764	il_rho(:)=1
765	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
766	IF( ANY(il_rho(:) < 0) )THEN
767	CALL logger_error("INTERP FILL VALUE: invalid "//&
768	& " refinement factor ")
769	ENDIF
770
771	il_offset(:,:)=0
772	IF( PRESENT(id_offset) )THEN
773	IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN
774	CALL logger_error("INTERP FILL VALUE: invalid array of offset")
775	ELSE
776	il_offset(:,:)=id_offset(:,:)
777	ENDIF
778	ENDIF
779
780	IF( (il_rho(jp_I) /= 1 .AND. td_var%t_dim(1)%l_use) .OR. &
781	& (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. &
782	& (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN
783
784	SELECT CASE(TRIM(td_var%c_interp(1)))
785	CASE('cubic','linear','nearest')
786	cl_method=TRIM(td_var%c_interp(1))
787	CASE DEFAULT
788	CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//&
789	& " use linear interpolation")
790	cl_method='linear'
791	! update variable structure value
792	td_var%c_interp(1)='linear'
793	END SELECT
794
795	CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//&
796	& " using "//TRIM(cl_method)//" method." )
797	CALL logger_info("INTERP FILL: refinement factor "//&
798	& TRIM(fct_str(il_rho(jp_I)))//&
799	& " "//TRIM(fct_str(il_rho(jp_J)))//&
800	& " "//TRIM(fct_str(il_rho(jp_K))) )
801
802	CALL interp__fill_value( td_var, cl_method, &
803	& il_rho(:), il_offset(:,:) )
804
805	SELECT CASE(TRIM(td_var%c_interp(2)))
806	CASE('/rhoi')
807	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
808	td_var%d_value(:,:,:,:) = &
809	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_I),dp)
810	END WHERE
811	CASE('/rhoj')
812	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
813	td_var%d_value(:,:,:,:) = &
814	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_J),dp)
815	END WHERE
816	CASE('/rhok')
817	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
818	td_var%d_value(:,:,:,:) = &
819	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_K),dp)
820	END WHERE
821	CASE('*rhoi')
822	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
823	td_var%d_value(:,:,:,:) = &
824	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_I),dp)
825	END WHERE
826	CASE('*rhoj')
827	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
828	td_var%d_value(:,:,:,:) = &
829	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_J),dp)
830	END WHERE
831	CASE('*rhok')
832	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
833	td_var%d_value(:,:,:,:) = &
834	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_K),dp)
835	END WHERE
836	CASE DEFAULT
837	td_var%c_interp(2)=''
838	END SELECT
839
840	ELSE
841	td_var%c_interp(:)=''
842	ENDIF
843
844	DEALLOCATE(il_rho)
845
846	END SUBROUTINE interp__fill_value_wrapper
847	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
848	SUBROUTINE interp__fill_value(td_var, cd_method, &
849	& id_rho, id_offset)
850	!-------------------------------------------------------------------
851	!> @brief
852	!> This subroutine interpolate value over mixed grid.
853	!>
854	!> @author J.Paul
855	!> @date November, 2013 - Initial Version
856	!> @date September, 2014
857	!> - use interpolation method modules
858	!>
859	!> @param[inout] td_var variable structure
860	!> @param[in] cd_method interpolation method
861	!> @param[in] id_rho array of refinment factor
862	!> @param[in] id_offset 2D array of offset between fine and coarse grid
863	!-------------------------------------------------------------------
864
865	IMPLICIT NONE
866
867	! Argument
868	TYPE(TVAR) , INTENT(INOUT) :: td_var
869	CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
870	INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho
871	INTEGER(I4) , DIMENSION(2,2), INTENT(IN ) :: id_offset
872
873	! local variable
874	CHARACTER(LEN=lc) :: cl_interp
875
876	INTEGER(I4) , DIMENSION(:) , ALLOCATABLE :: il_rho
877	INTEGER(i4) , DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
878
879	REAL(dp) :: dl_min
880	REAL(dp) :: dl_max
881
882	LOGICAL , DIMENSION(3) :: ll_even
883	LOGICAL :: ll_discont
884
885	TYPE(TVAR) :: tl_mix
886
887	TYPE(TATT) :: tl_att
888
889	! loop indices
890	!----------------------------------------------------------------
891
892	!1- create mixed grid
893	CALL interp_create_mixed_grid( td_var, tl_mix, id_rho(:) )
894
895	! clean variable structure
896	CALL var_clean(td_var)
897
898	!2- detect point to be interpolated
899	ALLOCATE( il_detect( tl_mix%t_dim(1)%i_len, &
900	& tl_mix%t_dim(2)%i_len, &
901	& tl_mix%t_dim(3)%i_len) )
902
903	il_detect(:,:,:)=0
904
905	il_detect(:,:,:)=interp_detect(tl_mix, id_rho(:) )
906
907	! add attribute to variable
908	cl_interp=fct_concat(tl_mix%c_interp(:))
909	tl_att=att_init('interpolation',cl_interp)
910	CALL var_move_att(tl_mix, tl_att)
911
912	! clean
913	CALL att_clean(tl_att)
914
915	! special case for even refinement on ARAKAWA-C grid
916	ll_even(:)=.FALSE.
917	IF( MOD(id_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
918	IF( MOD(id_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
919	IF( MOD(id_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
920
921	SELECT CASE(TRIM(tl_mix%c_point))
922	CASE('U')
923	ll_even(1)=.FALSE.
924	CASE('V')
925	ll_even(2)=.FALSE.
926	CASE('F')
927	ll_even(:)=.FALSE.
928	END SELECT
929
930	ALLOCATE(il_rho(ip_maxdim))
931	il_rho(:)=id_rho(:)
932
933	IF( ll_even(1) ) il_rho(jp_I)=id_rho(jp_I)+1
934	IF( ll_even(2) ) il_rho(jp_J)=id_rho(jp_J)+1
935	IF( ll_even(3) ) il_rho(jp_K)=id_rho(jp_K)+1
936
937	! special case for longitude
938	ll_discont=.FALSE.
939	IF( TRIM(tl_mix%c_units) == 'degrees_east' )THEN
940	dl_min=MINVAL( tl_mix%d_value(:,:,:,:), &
941	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
942	dl_max=MAXVAL( tl_mix%d_value(:,:,:,:), &
943	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
944	IF( dl_min < -170_dp .AND. dl_max > 170_dp .OR. &
945	& dl_min < 10_dp .AND. dl_max > 350_dp )THEN
946	ll_discont=.TRUE.
947	ENDIF
948	ENDIF
949
950	!3- interpolate
951	CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
952	& " discont "//TRIM(fct_str(ll_discont)) )
953	SELECT CASE(TRIM(cd_method))
954	CASE('cubic')
955	CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
956	& il_detect(:,:,:), &
957	& il_rho(:), ll_even(:), ll_discont )
958	CASE('nearest')
959	CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), &
960	& il_detect(:,:,:), &
961	& il_rho(:) )
962	CASE DEFAULT ! linear
963	CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
964	& il_detect(:,:,:), &
965	& il_rho(:), ll_even(:), ll_discont )
966	END SELECT
967
968	IF( ANY(il_detect(:,:,:)==1) )THEN
969	CALL logger_warn("INTERP FILL: some points can not be interpolated "//&
970	& "for variable "//TRIM(tl_mix%c_name) )
971	ENDIF
972
973	DEALLOCATE(il_detect)
974
975	!4- save useful domain (remove offset)
976	CALL interp_clean_mixed_grid( tl_mix, td_var, &
977	& id_rho(:), id_offset(:,:) )
978
979	! clean variable structure
980	DEALLOCATE(il_rho)
981	CALL var_clean(tl_mix)
982
983	END SUBROUTINE interp__fill_value
984	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
985	END MODULE interp

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

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

Context Navigation

source: utils/tools/SIREN/src/interp.f90 @ 13369

Download in other formats: