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.

interp.f90 in branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src/interp.f90 @ 4213

Last change on this file since 4213 was 4213, checked in by cbricaud, 10 years ago
first draft of the CONFIGURATION MANAGER demonstrator
File size: 79.9 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! MODULE: interp
6	!
7	! DESCRIPTION:
8	!> @brief
9	!> This module manage interpolation on regular grid.
10	!> @note It is used to work on ORCA grid, as we work only with grid indices.
11	!>
12	!> @warning due to the use of second derivative when using cubic interpolation
13	!> you should add 2 extrabands
14	!>
15	!> @details
16	!> @author
17	!> J.Paul
18	! REVISION HISTORY:
19	!> @date Nov, 2013 - Initial Version
20	!
21	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
22	!> @todo
23	!> - interp 3D
24	!> - see issue when fill value is zero for check cubic_fill..
25	!----------------------------------------------------------------------
26	MODULE interp
27
28	USE netcdf ! nf90 library
29	USE global ! global variable
30	USE kind ! F90 kind parameter
31	USE logger ! log file manager
32	USE fct ! basic useful function
33	! USE date ! date manager
34	USE att ! attribute manager
35	USE dim ! dimension manager
36	USE var ! variable manager
37	! USE file ! file manager
38	! USE iom ! I/O manager
39	! USE dom ! domain manager
40	USE grid ! grid manager
41	USE extrap ! extrapolation manager
42	! USE interp ! interpolation manager
43	! USE filter ! filter manager
44	! USE mpp ! MPP manager
45	! USE iom_mpp ! MPP I/O manager
46
47	IMPLICIT NONE
48	PRIVATE
49	! NOTE_avoid_public_variables_if_possible
50
51	! type and variable
52
53	! function and subroutine
54	PUBLIC :: interp_detect !< detected point to be interpolated
55	PUBLIC :: interp_fill_value !< interpolate value over detectected point
56	PUBLIC :: interp_create_mixed_grid !< create mixed grid
57	PUBLIC :: interp_clean_mixed_grid !< clean mixed grid
58
59	PRIVATE :: interp__detect !< detected point to be interpolated
60	PRIVATE :: interp__detect_wrapper !< detected point to be interpolated
61	PRIVATE :: interp__fill_value_wrapper !< interpolate value over detectected point
62	PRIVATE :: interp__fill_value !< interpolate value over detectected point
63	PRIVATE :: interp__clean_even_grid !< clean even mixed grid
64	PRIVATE :: interp__del_offset !< remove offset from interpolated grid
65	PRIVATE :: interp__check_method !< check if interpolation method available
66	! PRIVATE :: interp__3D !< interpolate 3D grid
67	PRIVATE :: interp__2D !< interpolate 2D grid
68	PRIVATE :: interp__1D !< interpolate 1D grid
69	PRIVATE :: interp__2D_cubic_coef !< compute coefficient for bicubic interpolation
70	PRIVATE :: interp__2D_cubic_fill !< compute bicubic interpolation
71	PRIVATE :: interp__2D_linear_coef !< compute coefficient for bilinear interpolation
72	PRIVATE :: interp__2D_linear_fill !< compute bilinear interpolation
73	PRIVATE :: interp__2D_nearest_fill !< compute nearest interpolation
74	PRIVATE :: interp__1D_cubic_coef !< compute coefficient for cubic interpolation
75	PRIVATE :: interp__1D_cubic_fill !< compute cubic interpolation
76	PRIVATE :: interp__1D_linear_coef !< compute coefficient for linear interpolation
77	PRIVATE :: interp__1D_linear_fill !< compute linear interpolation
78	PRIVATE :: interp__1D_nearest_fill !< compute nearest interpolation
79	! PRIVATE :: interp__longitude
80
81	TYPE TINTERP
82	!CHARACTER(LEN=lc) :: c_name = 'unknown' !< interpolation method name
83	!CHARACTER(LEN=lc) :: c_factor = 'unknown' !< interpolation factor
84	!CHARACTER(LEN=lc) :: c_divisor = 'unknown' !< interpolation divisor
85	CHARACTER(LEN=lc) :: c_name = '' !< interpolation method name
86	CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor
87	CHARACTER(LEN=lc) :: c_divisor= '' !< interpolation divisor
88	END TYPE TINTERP
89
90	INTERFACE interp_detect
91	MODULE PROCEDURE interp__detect_wrapper !< detected point to be interpolated
92	END INTERFACE interp_detect
93
94	INTERFACE interp_fill_value
95	MODULE PROCEDURE interp__fill_value_wrapper !< detected point to be interpolated
96	END INTERFACE interp_fill_value
97
98	CONTAINS
99	!-------------------------------------------------------------------
100	!> @brief
101	!> This function check if interpolation method available.
102	!>
103	!> @details
104	!>
105	!> @author J.Paul
106	!> - Nov, 2013- Initial Version
107	!
108	!> @param[in] cd_method : interpolation method
109	!> @return
110	!> @todo see extrap_detect
111	!-------------------------------------------------------------------
112	!> @code
113	FUNCTION interp__check_method( cd_method )
114	IMPLICIT NONE
115	! Argument
116	CHARACTER(LEN=lc) :: cd_method
117
118	! function
119	LOGICAL :: interp__check_method
120
121	! local variable
122	CHARACTER(LEN=lc) :: cl_interp
123	CHARACTER(LEN=lc) :: cl_method
124
125	! loop indices
126	INTEGER(I4) :: ji
127	!----------------------------------------------------------------
128
129	cl_method=fct_lower(cd_method)
130
131	interp__check_method=.FALSE.
132	DO ji=1,ig_ninterp
133	cl_interp=fct_lower(cg_interp_list(ji))
134	IF( TRIM(cl_interp) == TRIM(cl_method) )THEN
135	interp__check_method=.TRUE.
136	EXIT
137	ENDIF
138	ENDDO
139
140	END FUNCTION interp__check_method
141	!> @endcode
142	!-------------------------------------------------------------------
143	!> @brief
144	!> This function detected point to be interpolated.
145	!>
146	!> @details
147	!> Actually it checks, the number of dimension used for this variable
148	!> and launch interp__detect which detected point to be interpolated.
149	!>
150	!> @author J.Paul
151	!> - Nov, 2013- Initial Version
152	!
153	!> @param[in] td_mix : mixed grid variable (to interpolate)
154	!> @param[in] id_rho : table of refinement factor
155	!> @return table of detected point to be interpolated
156	!-------------------------------------------------------------------
157	!> @code
158	FUNCTION interp__detect_wrapper( td_mix, id_rho )
159	IMPLICIT NONE
160	! Argument
161	TYPE(TVAR) , INTENT(IN) :: td_mix
162	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
163
164	! function
165	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
166	& td_mix%t_dim(2)%i_len,&
167	& td_mix%t_dim(3)%i_len ) :: interp__detect_wrapper
168
169	! local variable
170
171	! loop indices
172	!----------------------------------------------------------------
173
174	IF( .NOT. ANY(td_mix%t_dim(1:3)%l_use) )THEN
175	! no dimension I-J-K used
176	CALL logger_debug(" INTERP DETECT: nothing done for variable"//&
177	& TRIM(td_mix%c_name) )
178
179	interp__detect_wrapper(:,:,:)=0
180
181	ELSE IF( ALL(td_mix%t_dim(1:3)%l_use) )THEN
182
183	! detect point to be interpolated on I-J-K
184	CALL logger_debug(" INTERP DETECT: detect point "//&
185	& TRIM(td_mix%c_point)//" for variable "//&
186	& TRIM(td_mix%c_name) )
187
188	interp__detect_wrapper(:,:,:)=interp__detect( td_mix, id_rho(:) )
189
190	ELSE IF( ALL(td_mix%t_dim(1:2)%l_use) )THEN
191
192	! detect point to be interpolated on I-J
193	CALL logger_debug(" INTERP DETECT: detect point "//&
194	& TRIM(td_mix%c_point)//" for variable "//&
195	& TRIM(td_mix%c_name) )
196
197	interp__detect_wrapper(:,:,1:1)=interp__detect( td_mix, id_rho(:))
198
199	ELSE IF( td_mix%t_dim(3)%l_use )THEN
200
201	! detect point to be interpolated on K
202	CALL logger_debug(" INTERP DETECT: detect vertical point "//&
203	& " for variable "//TRIM(td_mix%c_name) )
204
205	interp__detect_wrapper(1:1,1:1,:)=interp__detect( td_mix, id_rho(:) )
206
207	ENDIF
208
209	END FUNCTION interp__detect_wrapper
210	!> @endcode
211	!-------------------------------------------------------------------
212	!> @brief
213	!> This function detected point to be interpolated.
214	!>
215	!> @details
216	!> A special case is done for even refinement on ARAKAWA-C grid.
217	!>
218	!> @author J.Paul
219	!> - Nov, 2013- Initial Version
220	!
221	!> @param[in] td_mix : mixed grid variable (to interpolate)
222	!> @param[in] id_rho : table of refinement factor
223	!> @return table of detected point to be interpolated
224	!-------------------------------------------------------------------
225	!> @code
226	FUNCTION interp__detect( td_mix, id_rho )
227	IMPLICIT NONE
228	! Argument
229	TYPE(TVAR) , INTENT(IN) :: td_mix
230	INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
231
232	! function
233	INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
234	& td_mix%t_dim(2)%i_len,&
235	& td_mix%t_dim(3)%i_len ) :: interp__detect
236
237	! local variable
238	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
239
240	INTEGER(I4) :: il_xextra
241	INTEGER(I4) :: il_yextra
242	INTEGER(I4) :: il_zextra
243
244	INTEGER(i4), DIMENSION(3) :: il_dim
245
246	LOGICAL , DIMENSION(3) :: ll_even
247
248	! loop indices
249	INTEGER(I4) :: ji
250	INTEGER(I4) :: jj
251	INTEGER(I4) :: jk
252	!----------------------------------------------------------------
253	ALLOCATE( il_rho(ip_maxdim) )
254	il_rho(:)=1
255	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
256
257	! special case for even refinement on ARAKAWA-C grid
258	ll_even(:)=.FALSE.
259	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
260	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
261	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
262
263	SELECT CASE(TRIM(td_mix%c_point))
264	CASE('U')
265	ll_even(1)=.FALSE.
266	CASE('V')
267	ll_even(2)=.FALSE.
268	CASE('F')
269	ll_even(:)=.FALSE.
270	END SELECT
271
272	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
273	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
274	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
275
276	! special case for cubic interpolation
277	il_xextra=0
278	il_yextra=0
279	il_zextra=0
280	SELECT CASE(TRIM(td_mix%c_interp(1)))
281	CASE('cubic')
282	! those points can not be compute cause cubic interpolation
283	! need second derivative.
284	IF( il_rho(jp_I) /= 1 ) il_xextra=3*il_rho(jp_I)
285	IF( il_rho(jp_J) /= 1 ) il_yextra=3*il_rho(jp_J)
286	IF( il_rho(jp_K) /= 1 ) il_zextra=3*il_rho(jp_K)
287	END SELECT
288
289	il_dim(:)=td_mix%t_dim(1:3)%i_len
290
291	interp__detect(:,:,:)=1
292
293	! do not compute coarse grid point
294	interp__detect(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
295	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
296	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K) ) = 0
297
298	! do not compute point near fill value
299	FORALL( ji=1:il_dim(1):il_rho(jp_I), &
300	& jj=1:il_dim(2):il_rho(jp_J), &
301	& jk=1:il_dim(3):il_rho(jp_K), &
302	& td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )
303
304	! i-direction
305	interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
306	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
307	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
308	! j-direction
309	interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
310	& MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
311	& MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
312	! k-direction
313	interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
314	& MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
315	& MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0
316
317	END FORALL
318
319	END FUNCTION interp__detect
320	!> @endcode
321	!-------------------------------------------------------------------
322	!> @brief
323	!> This subroutine create mixed grid.
324	!>
325	!> @details
326	!> Created grid is fine resolution grid.
327	!> First and last point are coasre grid point.
328	!>
329	!> A special case is done for even refinement on ARAKAWA-C grid.
330	!>
331	!> @author J.Paul
332	!> - Nov, 2013- Initial Version
333	!
334	!> @param[in] td_var : coarse grid variable (should be extrapolated)
335	!> @param[out] td_mix : mixed grid variable
336	!> @param[in] id_rho : table of refinment factor
337	!>
338	!> @todo
339	!-------------------------------------------------------------------
340	!> @code
341	SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho )
342	IMPLICIT NONE
343	! Argument
344	TYPE(TVAR) , INTENT(IN ) :: td_var
345	TYPE(TVAR) , INTENT( OUT) :: td_mix
346	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
347
348	! local variable
349	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
350
351	INTEGER(i4) :: il_xextra
352	INTEGER(i4) :: il_yextra
353	INTEGER(i4) :: il_zextra
354
355	LOGICAL, DIMENSION(3) :: ll_even
356
357	! loop indices
358	!----------------------------------------------------------------
359	ALLOCATE(il_rho(ip_maxdim))
360	il_rho(:)=1
361	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
362
363	! special case for even refinement on ARAKAWA-C grid
364	ll_even(:)=.FALSE.
365	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
366	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
367	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
368
369	SELECT CASE(TRIM(td_var%c_point))
370	CASE('U')
371	ll_even(1)=.FALSE.
372	CASE('V')
373	ll_even(2)=.FALSE.
374	CASE('F')
375	ll_even(:)=.FALSE.
376	END SELECT
377
378	IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
379	IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
380	IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
381
382	! copy variable
383	td_mix=td_var
384
385	! compute new dimension length
386	il_xextra=il_rho(jp_I)-1
387	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
388
389	il_yextra=il_rho(jp_J)-1
390	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
391
392	il_zextra=il_rho(jp_K)-1
393	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
394
395	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
396	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
397	& td_mix%t_dim(2)%i_len, &
398	& td_mix%t_dim(3)%i_len, &
399	& td_mix%t_dim(4)%i_len) )
400
401	! initialise to FillValue
402	td_mix%d_value(:,:,:,:)=td_mix%d_fill
403
404	! quid qd coord ou bathy fourni par user ?? (offset ??)
405	td_mix%d_value(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
406	& 1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
407	& 1:td_mix%t_dim(3)%i_len:il_rho(jp_K), :) = &
408	& td_var%d_value(:,:,:,:)
409
410	END SUBROUTINE interp_create_mixed_grid
411	!> @endcode
412	!-------------------------------------------------------------------
413	!> @brief
414	!> This subroutine remove points added to mixed grid to compute
415	!> interpolation in the special case of even refinement on ARAKAWA-C grid.
416	!>
417	!> @details
418	!>
419	!> @author J.Paul
420	!> - Nov, 2013- Initial Version
421	!
422	!> @param[in] td_mix : mixed grid variable
423	!> @param[in] id_rho : table of refinment factor
424	!>
425	!> @todo
426	!-------------------------------------------------------------------
427	!> @code
428	SUBROUTINE interp__clean_even_grid( td_mix, id_rho )
429	IMPLICIT NONE
430	! Argument
431	TYPE(TVAR) , INTENT(INOUT) :: td_mix
432	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
433
434	! local variable
435	INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
436
437	INTEGER(i4) :: il_xextra
438	INTEGER(i4) :: il_yextra
439	INTEGER(i4) :: il_zextra
440
441	LOGICAL, DIMENSION(3) :: ll_even
442
443	LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ll_mask
444
445	REAL(dp), DIMENSION(:) , ALLOCATABLE :: dl_vect
446
447	TYPE(TVAR) :: tl_mix
448
449	! loop indices
450	!----------------------------------------------------------------
451	ALLOCATE(il_rho(ip_maxdim))
452	il_rho(:)=1
453	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
454
455	! special case for even refinement on ARAKAWA-C grid
456	ll_even(:)=.FALSE.
457	IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
458	IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
459	IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
460
461	SELECT CASE(TRIM(td_mix%c_point))
462	CASE('U')
463	ll_even(1)=.FALSE.
464	CASE('V')
465	ll_even(2)=.FALSE.
466	CASE('F')
467	ll_even(:)=.FALSE.
468	END SELECT
469
470	! copy variable
471	tl_mix=td_mix
472
473	! remove some point only if refinement in some direction is even
474	IF( ANY(ll_even(:)) )THEN
475
476	ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, &
477	& tl_mix%t_dim(2)%i_len, &
478	& tl_mix%t_dim(3)%i_len, &
479	& tl_mix%t_dim(4)%i_len) )
480
481	ll_mask(:,:,:,:)=.TRUE.
482
483	IF( tl_mix%t_dim(1)%l_use .AND. ll_even(1) )THEN
484
485	il_rho(jp_I)=il_rho(jp_I)+1
486
487	! locate wrong point on mixed grid
488	ll_mask(1:td_mix%t_dim(1)%i_len:il_rho(jp_I),:,:,:)=.FALSE.
489
490	! compute coasre grid dimension length
491	il_xextra=il_rho(jp_I)-1
492	td_mix%t_dim(1)%i_len=(tl_mix%t_dim(1)%i_len+il_xextra)/il_rho(jp_I)
493
494	il_rho(jp_I)=il_rho(jp_I)-1
495	! compute right fine grid dimension length
496	td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
497
498	ENDIF
499
500	IF( tl_mix%t_dim(2)%l_use .AND. ll_even(2) )THEN
501
502	il_rho(jp_J)=il_rho(jp_J)+1
503
504	! locate wrong point on mixed grid
505	ll_mask(:,1:tl_mix%t_dim(2)%i_len:il_rho(jp_J),:,:)=.FALSE.
506
507	! compute coasre grid dimension length
508	il_yextra=il_rho(jp_J)-1
509	td_mix%t_dim(2)%i_len=(tl_mix%t_dim(2)%i_len+il_yextra)/il_rho(jp_J)
510
511	il_rho(jp_J)=il_rho(jp_J)-1
512	! compute right fine grid dimension length
513	td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
514
515	ENDIF
516
517	IF( tl_mix%t_dim(3)%l_use .AND. ll_even(3) )THEN
518
519	il_rho(jp_K)=il_rho(jp_K)+1
520
521	! locate wrong point on mixed grid
522	ll_mask(:,:,1:tl_mix%t_dim(3)%i_len:il_rho(jp_K),:)=.FALSE.
523
524	! compute coasre grid dimension length
525	il_zextra=il_rho(jp_K)-1
526	td_mix%t_dim(3)%i_len=(tl_mix%t_dim(3)%i_len+il_zextra)/il_rho(jp_K)
527
528	il_rho(jp_K)=il_rho(jp_K)-1
529	! compute right fine grid dimension length
530	td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
531
532	ENDIF
533
534	IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
535	ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
536	& td_mix%t_dim(2)%i_len, &
537	& td_mix%t_dim(3)%i_len, &
538	& td_mix%t_dim(4)%i_len) )
539
540	! initialise to FillValue
541	td_mix%d_value(:,:,:,:)=td_mix%d_fill
542
543	IF( COUNT(ll_mask(:,:,:,:)) /= SIZE(td_mix%d_value(:,:,:,:)) )THEN
544
545	CALL logger_error("INTERP CLEAN EVEN GRID: output value size "//&
546	& " and mask count differ ")
547
548	ELSE
549
550	ALLOCATE( dl_vect(COUNT(ll_mask(:,:,:,:))) )
551
552	dl_vect(:)= PACK( tl_mix%d_value(:,:,:,:), &
553	& MASK=ll_mask(:,:,:,:) )
554
555	td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), &
556	& SHAPE=td_mix%t_dim(:)%i_len )
557
558	DEALLOCATE( dl_vect )
559
560	ENDIF
561
562	DEALLOCATE( ll_mask )
563
564	ENDIF
565
566	CALL var_clean(tl_mix)
567
568	END SUBROUTINE interp__clean_even_grid
569	!> @endcode
570	!-------------------------------------------------------------------
571	!> @brief
572	!> This subroutine save interpolated value over domain.
573	!> And so remove points added on mixed grid
574	!> to compute interpolation
575	!>
576	!> @details
577	!>
578	!> @author J.Paul
579	!> - Nov, 2013- Initial Version
580	!
581	!> @param[in] td_var : table of mixed grid variable (to interpolate)
582	!> @param[in] id_rho : table of refinement factor
583	!> @return
584	!-------------------------------------------------------------------
585	!> @code
586	SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, &
587	& id_rho, &
588	& id_offset )
589	IMPLICIT NONE
590	! Argument
591	TYPE(TVAR) , INTENT(IN ) :: td_mix
592	TYPE(TVAR) , INTENT( OUT) :: td_var
593	INTEGER(I4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
594	INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
595
596	! local variable
597	INTEGER(i4) :: il_imin0
598	INTEGER(i4) :: il_imax0
599	INTEGER(i4) :: il_jmin0
600	INTEGER(i4) :: il_jmax0
601
602	INTEGER(i4) :: il_imin1
603	INTEGER(i4) :: il_jmin1
604	INTEGER(i4) :: il_imax1
605	INTEGER(i4) :: il_jmax1
606
607	INTEGER(i4), DIMENSION(2,2) :: il_offset
608
609	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
610
611	TYPE(TVAR) :: tl_mix
612
613	! loop indices
614	!----------------------------------------------------------------
615	il_offset(:,:)=0
616	IF( PRESENT(id_offset) )THEN
617	IF( ANY( SHAPE(id_offset(:,:)) /= SHAPE(il_offset(:,:)) ) )THEN
618	CALL logger_error("INTERP CLEAN MIXED GRID: invalid dimension of"//&
619	& " offset table")
620	ELSE
621	il_offset(:,:)=id_offset(:,:)
622	ENDIF
623	ENDIF
624	! copy mixed variable in temporary structure
625	tl_mix=td_mix
626
627	! remove unusefull points over mixed grid for even refinement
628	CALL interp__clean_even_grid(tl_mix, id_rho(:))
629
630	! copy cleaned mixed variable
631	td_var=tl_mix
632	IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
633
634	! delete table of value
635	CALL var_del_value(td_var)
636
637	! compute domain indices
638	il_imin0=1 ; il_imax0=td_var%t_dim(1)%i_len
639	il_jmin0=1 ; il_jmax0=td_var%t_dim(2)%i_len
640
641	il_imin1=il_imin0+il_offset(1,1)
642	il_jmin1=il_jmin0+il_offset(2,1)
643
644	il_imax1=il_imax0-il_offset(1,2)
645	il_jmax1=il_jmax0-il_offset(2,2)
646
647	SELECT CASE(TRIM(td_var%c_point))
648	CASE('U')
649	il_imin1=il_imin0+(il_offset(1,1)-1)
650	il_imax1=il_imax0-(il_offset(1,2)+1)
651	CASE('V')
652	il_jmin1=il_jmin0+(il_offset(2,1)-1)
653	il_jmax1=il_jmax0-(il_offset(2,2)+1)
654	CASE('F')
655	il_imin1=il_imin0+(il_offset(1,1)-1)
656	il_imax1=il_imax0-(il_offset(1,2)+1)
657
658	il_jmin1=il_jmin0+(il_offset(2,1)-1)
659	il_jmax1=il_jmax0-(il_offset(2,2)+1)
660	END SELECT
661
662	! compute new dimension
663	td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
664	td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
665
666	ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
667	& td_var%t_dim(2)%i_len, &
668	& td_var%t_dim(3)%i_len, &
669	& td_var%t_dim(4)%i_len) )
670
671	dl_value( 1:td_var%t_dim(1)%i_len, &
672	& :,:,:) = tl_mix%d_value( il_imin1:il_imax1, &
673	& il_jmin1:il_jmax1, &
674	& :, : )
675
676	! add variable value
677	CALL var_add_value(td_var,dl_value(:,:,:,:))
678
679	! save variable type
680	td_var%i_type=tl_mix%i_type
681
682	DEALLOCATE(dl_value)
683
684	ENDIF
685
686	CALL var_clean(tl_mix)
687
688	END SUBROUTINE interp_clean_mixed_grid
689	!> @endcode
690	!-------------------------------------------------------------------
691	!> @brief
692	!> This subroutine save interpolated value over domain.
693	!> And so remove points added on mixed grid
694	!> to compute interpolation
695	!>
696	!> @details
697	!>
698	!> @author J.Paul
699	!> - Nov, 2013- Initial Version
700	!
701	!> @param[in] td_var : table of mixed grid variable (to interpolate)
702	!> @param[in] id_offset : table of offset
703	!> @return
704	!-------------------------------------------------------------------
705	!> @code
706	SUBROUTINE interp__del_offset( td_var, id_offset )
707	IMPLICIT NONE
708	! Argument
709	TYPE(TVAR) , INTENT(INOUT) :: td_var
710	INTEGER(i4), DIMENSION(:,:), INTENT(IN ) :: id_offset
711
712	! local variable
713	INTEGER(i4) :: il_imin1
714	INTEGER(i4) :: il_jmin1
715	INTEGER(i4) :: il_imax1
716	INTEGER(i4) :: il_jmax1
717
718	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
719
720	! loop indices
721	!----------------------------------------------------------------
722
723	IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
724
725	il_imin1=1+id_offset(1,1)
726	il_jmin1=1+id_offset(2,1)
727
728	il_imax1=td_var%t_dim(1)%i_len-id_offset(2,1)
729	il_jmax1=td_var%t_dim(2)%i_len-id_offset(2,2)
730
731	! compute new dimension
732	td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
733	td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
734	ALLOCATE( dl_value( td_var%t_dim(1)%i_len, &
735	& td_var%t_dim(2)%i_len, &
736	& td_var%t_dim(3)%i_len, &
737	& td_var%t_dim(4)%i_len) )
738
739	dl_value(:,:,:,:)=td_var%d_value( il_imin1:il_imax1, &
740	& il_jmin1:il_jmax1, &
741	& :,: )
742
743	! delete table of value
744	CALL var_del_value(td_var)
745
746	! add variable value
747	CALL var_add_value(td_var,dl_value(:,:,:,:))
748
749	DEALLOCATE(dl_value)
750
751	ENDIF
752
753	END SUBROUTINE interp__del_offset
754	!> @endcode
755	!-------------------------------------------------------------------
756	!> @brief
757	!> This subroutine interpolate detected point.
758	!>
759	!> @details
760	!> Actually it checks, the number of dimension used for this variable
761	!> and launch interp__fill_value.
762	!>
763	!> @author J.Paul
764	!> - Nov, 2013- Initial Version
765	!
766	!> @param[in] td_var : variable to be interpolated
767	!> @param[in] id_rho : table of refinement factor
768	!-------------------------------------------------------------------
769	!> @code
770	SUBROUTINE interp__fill_value_wrapper( td_var, &
771	& id_rho, &
772	& id_offset )
773	IMPLICIT NONE
774	! Argument
775	TYPE(TVAR) , INTENT(INOUT) :: td_var
776	INTEGER(I4), DIMENSION(:), INTENT(IN ), OPTIONAL :: id_rho
777	INTEGER(I4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
778
779	! local variable
780	CHARACTER(LEN=lc) :: cl_method
781	INTEGER(i4) , DIMENSION(:), ALLOCATABLE :: il_rho
782	INTEGER(i4) , DIMENSION(2,2) :: il_offset
783
784	! loop indices
785	!----------------------------------------------------------------
786
787	SELECT CASE(TRIM(td_var%c_interp(1)))
788	CASE('cubic','linear','nearest')
789	cl_method=TRIM(td_var%c_interp(1))
790	CASE DEFAULT
791	CALL logger_warn("INTERP FILL: interpolation method unknown."//&
792	& " use linear interpolation")
793	cl_method='linear'
794
795	! update variable structure value
796	td_var%c_interp(1)='linear'
797	END SELECT
798
799	ALLOCATE( il_rho(ip_maxdim) )
800	il_rho(:)=1
801	IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
802	IF( ANY(il_rho(:) < 0) )THEN
803	CALL logger_error("INTERP FILL VALUE: invalid "//&
804	& " refinement factor ")
805	ENDIF
806
807	il_offset(:,:)=0
808	IF( PRESENT(id_offset) )THEN
809	IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN
810	CALL logger_error("INTERP FILL VALUE: invalid table of offset")
811	ELSE
812	il_offset(:,:)=id_offset(:,:)
813	ENDIF
814	ENDIF
815
816	IF( (il_rho(jp_I) /= 1 .AND. td_var%t_dim(1)%l_use) .OR. &
817	& (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. &
818	& (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN
819
820	CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//&
821	& " using "//TRIM(cl_method)//" method." )
822	CALL logger_info("INTERP FILL: refinement factor "//&
823	& TRIM(fct_str(il_rho(jp_I)))//&
824	& " "//TRIM(fct_str(il_rho(jp_J)))//&
825	& " "//TRIM(fct_str(il_rho(jp_K))) )
826
827	CALL interp__fill_value( td_var, cl_method, &
828	& il_rho(:), &
829	& il_offset(:,:) )
830
831	SELECT CASE(TRIM(td_var%c_interp(2)))
832	CASE('/rhoi')
833	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
834	td_var%d_value(:,:,:,:) = &
835	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_I),dp)
836	END WHERE
837	CASE('/rhoj')
838	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
839	td_var%d_value(:,:,:,:) = &
840	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_J),dp)
841	END WHERE
842	CASE('/rhok')
843	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
844	td_var%d_value(:,:,:,:) = &
845	& td_var%d_value(:,:,:,:) / REAL(il_rho(jp_K),dp)
846	END WHERE
847	CASE('*rhoi')
848	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
849	td_var%d_value(:,:,:,:) = &
850	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_I),dp)
851	END WHERE
852	CASE('*rhoj')
853	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
854	td_var%d_value(:,:,:,:) = &
855	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_J),dp)
856	END WHERE
857	CASE('*rhok')
858	WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
859	td_var%d_value(:,:,:,:) = &
860	& td_var%d_value(:,:,:,:) * REAL(il_rho(jp_K),dp)
861	END WHERE
862	CASE DEFAULT
863	td_var%c_interp(2)=''
864	END SELECT
865
866	ENDIF
867
868	END SUBROUTINE interp__fill_value_wrapper
869	!> @endcode
870	!-------------------------------------------------------------------
871	!> @brief
872	!> This subroutine interpolate value over mixed grid.
873	!>
874	!> @details
875	!>
876	!>
877	!> @author J.Paul
878	!> - Nov, 2013- Initial Version
879	!
880	!> @param[inout] td_var : variable
881	!> @param[in] cd_method : interpolation method
882	!> @param[in] id_rho : table of refinment factor
883	!-------------------------------------------------------------------
884	!> @code
885	SUBROUTINE interp__fill_value( td_var, cd_method, &
886	& id_rho, &
887	& id_offset )
888	IMPLICIT NONE
889	! Argument
890	TYPE(TVAR) , INTENT(INOUT) :: td_var
891	CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
892	INTEGER(I4) , DIMENSION(:) , INTENT(IN ) :: id_rho
893	INTEGER(I4) , DIMENSION(:,:), INTENT(IN ) :: id_offset
894
895	! local variable
896	CHARACTER(LEN=lc) :: cl_interp
897
898	INTEGER(I4) , DIMENSION(:) , ALLOCATABLE :: il_rho
899	INTEGER(i4) , DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
900
901	REAL(dp) :: dl_min
902	REAL(dp) :: dl_max
903
904	LOGICAL , DIMENSION(3) :: ll_even
905	LOGICAL :: ll_discont
906
907	TYPE(TVAR) :: tl_mix
908
909	TYPE(TATT) :: tl_att
910
911	! loop indices
912	INTEGER(i4) :: ji
913	INTEGER(i4) :: jj
914	INTEGER(i4) :: jk
915	INTEGER(i4) :: jl
916	!----------------------------------------------------------------
917
918	!1- create mixed grid
919	CALL interp_create_mixed_grid( td_var, tl_mix, id_rho(:) )
920
921	! clean variable structure
922	CALL var_clean(td_var)
923
924	!2- detect point to be interpolated
925	ALLOCATE( il_detect( tl_mix%t_dim(1)%i_len, &
926	& tl_mix%t_dim(2)%i_len, &
927	& tl_mix%t_dim(3)%i_len) )
928
929	il_detect(:,:,:)=0
930
931	il_detect(:,:,:)=interp_detect(tl_mix, id_rho(:) )
932
933	! add attribute to variable
934	cl_interp=fct_concat(tl_mix%c_interp(:))
935	tl_att=att_init('interpolation',cl_interp)
936	CALL var_move_att(tl_mix, tl_att)
937
938	! special case for even refinement on ARAKAWA-C grid
939	ll_even(:)=.FALSE.
940	IF( MOD(id_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
941	IF( MOD(id_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
942	IF( MOD(id_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.
943
944	SELECT CASE(TRIM(tl_mix%c_point))
945	CASE('U')
946	ll_even(1)=.FALSE.
947	CASE('V')
948	ll_even(2)=.FALSE.
949	CASE('F')
950	ll_even(:)=.FALSE.
951	END SELECT
952
953	ALLOCATE(il_rho(ip_maxdim))
954	il_rho(:)=id_rho(:)
955
956	IF( ll_even(1) ) il_rho(jp_I)=id_rho(jp_I)+1
957	IF( ll_even(2) ) il_rho(jp_J)=id_rho(jp_J)+1
958	IF( ll_even(3) ) il_rho(jp_K)=id_rho(jp_K)+1
959
960	! special case for longitude
961	ll_discont=.FALSE.
962	IF( TRIM(tl_mix%c_units) == 'degrees_east' )THEN
963	dl_min=MINVAL( tl_mix%d_value(:,:,:,:), &
964	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
965	dl_max=MAXVAL( tl_mix%d_value(:,:,:,:), &
966	& tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
967	IF( dl_min < -170_dp .AND. dl_max > 170_dp .OR. &
968	& dl_min < 10_dp .AND. dl_max > 350_dp )THEN
969	ll_discont=.TRUE.
970	ENDIF
971	ENDIF
972
973	!3- interpolate
974	DO jl=1,tl_mix%t_dim(4)%i_len
975	IF( il_rho(jp_K) /= 1 )THEN
976	! CALL interp__3D(tl_mix%d_value(:,:,:,jl), tl_mix%d_fill, &
977	! & il_detect(:,:,:), cd_method, &
978	! & il_rhoi, il_rhoj, il_rhok, &
979	! & ll_even(:), ll_discont )
980	CALL logger_error("INTERP FILL: can not interpolate "//&
981	& "vertically for now ")
982	ENDIF
983
984	IF( ANY(il_detect(:,:,:)==1) )THEN
985	! I-J plan
986	DO jk=1,tl_mix%t_dim(3)%i_len
987	CALL interp__2D(tl_mix%d_value(:,:,jk,jl), tl_mix%d_fill, &
988	& il_detect(:,:,jk), cd_method, &
989	& il_rho(jp_I), il_rho(jp_J), &
990	& ll_even(1:2), ll_discont)
991	ENDDO
992	IF( ALL(il_detect(:,:,:)==0) ) CYCLE
993	IF( il_rho(jp_K) /= 1 )THEN
994	! I-K plan
995	DO jj=1,tl_mix%t_dim(2)%i_len
996	CALL interp__2D(tl_mix%d_value(:,jj,:,jl), tl_mix%d_fill, &
997	& il_detect(:,jj,:), cd_method, &
998	& il_rho(jp_J), il_rho(jp_K), &
999	& ll_even(1:3:2), ll_discont )
1000	ENDDO
1001	IF( ALL(il_detect(:,:,:)==0) ) CYCLE
1002	! J-K plan
1003	DO ji=1,tl_mix%t_dim(1)%i_len
1004	CALL interp__2D(tl_mix%d_value(ji,:,:,jl), tl_mix%d_fill, &
1005	& il_detect(ji,:,:), cd_method, &
1006	& il_rho(jp_J), il_rho(jp_K), &
1007	& ll_even(2:3), ll_discont )
1008	ENDDO
1009
1010	ENDIF
1011	IF( ANY(il_detect(:,:,:)==1) )THEN
1012	! I direction
1013	DO jk=1,tl_mix%t_dim(3)%i_len
1014	DO jj=1,tl_mix%t_dim(2)%i_len
1015	CALL interp__1D( tl_mix%d_value(:,jj,jk,jl), &
1016	& tl_mix%d_fill, &
1017	& il_detect(:,jj,jk), cd_method, &
1018	& il_rho(jp_I), &
1019	& ll_even(1), ll_discont )
1020	ENDDO
1021	ENDDO
1022	IF( ALL(il_detect(:,:,:)==0) ) CYCLE
1023	! J direction
1024	DO jk=1,tl_mix%t_dim(3)%i_len
1025	DO ji=1,tl_mix%t_dim(1)%i_len
1026	CALL interp__1D( tl_mix%d_value(ji,:,jk,jl), &
1027	& tl_mix%d_fill, &
1028	& il_detect(ji,:,jk), cd_method, &
1029	& il_rho(jp_J), &
1030	& ll_even(2), ll_discont )
1031	ENDDO
1032	ENDDO
1033	IF( il_rho(jp_K) /= 1 )THEN
1034	IF( ALL(il_detect(:,:,:)==0) ) CYCLE
1035	! K direction
1036	DO jj=1,tl_mix%t_dim(2)%i_len
1037	DO ji=1,tl_mix%t_dim(1)%i_len
1038	CALL interp__1D( tl_mix%d_value(ji,jj,:,jl), &
1039	& tl_mix%d_fill, &
1040	& il_detect(ji,jj,:), cd_method, &
1041	& il_rho(jp_K), &
1042	& ll_even(3), ll_discont )
1043	ENDDO
1044	ENDDO
1045	ENDIF
1046	ENDIF
1047	ENDIF
1048	ENDDO
1049
1050	IF( ANY(il_detect(:,:,:)==1) )THEN
1051	CALL logger_warn("INTERP FILL: some points can not be interpolated "//&
1052	& "for variable "//TRIM(tl_mix%c_name) )
1053	ENDIF
1054
1055	DEALLOCATE(il_detect)
1056	!4- save useful domain (remove offset)
1057	CALL interp_clean_mixed_grid( tl_mix, td_var, &
1058	& id_rho(:), &
1059	& id_offset(:,:) )
1060
1061	! clean variable structure
1062	CALL var_clean(tl_mix)
1063
1064	END SUBROUTINE interp__fill_value
1065	!> @endcode
1066	! !-------------------------------------------------------------------
1067	! !> @brief
1068	! !> This subroutine
1069	! !>
1070	! !> @details
1071	! !>
1072	! !> @author J.Paul
1073	! !> - Nov, 2013- Initial Version
1074	! !
1075	! !> @param[in]
1076	! !> @param[out]
1077	! !-------------------------------------------------------------------
1078	! !> @code
1079	! SUBROUTINE interp__3D( dd_value, dd_fill, &
1080	! & id_detect, cd_method, &
1081	! & id_rhoi, id_rhoj, id_rhok, &
1082	! & ld_even, ld_discont )
1083	! IMPLICIT NONE
1084	! ! Argument
1085	! REAL(dp) , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value
1086	! REAL(dp) , INTENT(IN ) :: dd_fill
1087	! INTEGER(I4) , DIMENSION(:,:,:), INTENT(INOUT) :: id_detect
1088	! CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
1089	! INTEGER(I4) , INTENT(IN ) :: id_rhoi
1090	! INTEGER(I4) , INTENT(IN ) :: id_rhoj
1091	! INTEGER(I4) , INTENT(IN ) :: id_rhok
1092	! LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even
1093	! LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont
1094	!
1095	! ! local variable
1096	! INTEGER(i4), DIMENSION(3) :: il_shape
1097	! INTEGER(i4), DIMENSION(3) :: il_dim
1098	!
1099	! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coarse
1100	! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
1101	! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy
1102	! REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_d2fdxy
1103	!
1104	! LOGICAL :: ll_discont
1105	!
1106	! ! loop indices
1107	!! INTEGER(i4) :: ji
1108	!! INTEGER(i4) :: jj
1109	!! INTEGER(i4) :: ii
1110	!! INTEGER(i4) :: ij
1111	! !----------------------------------------------------------------
1112	! ll_discont=.FALSE.
1113	! IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1114	!
1115	! il_shape(:)=SHAPE(dd_value)
1116	! il_dim(:)=(il_shape(:)-1)/2
1117	!
1118	! ALLOCATE( dl_coarse(il_dim(1),il_dim(2),il_dim(3)) )
1119	! ! value on coarse grid
1120	! dl_coarse(:,:,:)=dd_value( 1:il_shape(1):id_rhoi, &
1121	! & 1:il_shape(2):id_rhoj, &
1122	! & 1:il_shape(3):id_rhok )
1123	! SELECT CASE(TRIM(cd_method))
1124	!
1125	! CASE('cubic')
1126	!
1127	! ALLOCATE( dl_dfdx( il_dim(1),il_dim(2),il_dim(3)), &
1128	! & dl_dfdy( il_dim(1),il_dim(2),il_dim(3)), &
1129	! & dl_d2fdxy(il_dim(1),il_dim(2),il_dim(3)) )
1130	!
1131	!! ! compute derivative on coarse grid
1132	!! dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'I')
1133	!! dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:,:), dd_fill, 'J')
1134	!!
1135	!! ! compute cross derivative on coarse grid
1136	!! dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:,:), dd_fill, 'J')
1137	!!
1138	!! DO jj=1,il_shape(2)-1,id_rhoj
1139	!! ij=((jj-1)/id_rhoj)+1
1140	!! DO ji=1,il_shape(1)-1,id_rhoi
1141	!! ii=((ji-1)/id_rhoi)+1
1142	!!
1143	!! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
1144	!! ! compute bicubic coefficient
1145	!! dl_coef(:)=interp__2D_cubic_coef(dl_coarse(ii:ii+1,ij:ij+1),&
1146	!! & dl_dfdx( ii:ii+1,ij:ij+1),&
1147	!! & dl_dfdy( ii:ii+1,ij:ij+1),&
1148	!! & dl_d2fdxy(ii:ii+1,ij:ij+1) )
1149	!!
1150	!! ! compute value on detetected point
1151	!! CALL interp__2D_cubic_fill(dl_coef(:), &
1152	!! & dd_value( ji:ji+id_rhoi, &
1153	!! & jj:jj+id_rhoj ), &
1154	!! & id_detect(ji:ji+id_rhoi, &
1155	!! & jj:jj+id_rhoj ) )
1156	!!
1157	!! ENDDO
1158	!! ENDDO
1159	!
1160	! DEALLOCATE( dl_dfdx, &
1161	! & dl_dfdy, &
1162	! & dl_d2fdxy )
1163	!
1164	! CASE('nearest')
1165	!
1166	! CASE DEFAULT ! linear
1167	!
1168	!! DO jj=1,il_shape(2)-1,id_rhoj
1169	!! ij=((jj-1)/id_rhoj)+1
1170	!! DO ji=1,il_shape(1)-1,id_rhoi
1171	!! ii=((ji-1)/id_rhoi)+1
1172	!!
1173	!! IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
1174	!! ! compute bilinear coefficient
1175	!! dl_coef(:)=interp__2D_linear_coef(dl_coarse(ii:ii+1,ij:ij+1))
1176	!!
1177	!! ! compute value on detetected point
1178	!! CALL interp__2D_linear_fill(dl_coef(:), &
1179	!! & dd_value( ji:ji+id_rhoi, &
1180	!! & jj:jj+id_rhoj ), &
1181	!! & id_detect(ji:ji+id_rhoi, &
1182	!! & jj:jj+id_rhoj ) )
1183	!!
1184	!! ENDDO
1185	!! ENDDO
1186	!
1187	! END SELECT
1188	!
1189	! DEALLOCATE( dl_coarse )
1190	!
1191	! END SUBROUTINE interp__3D
1192	! !> @endcode
1193	!-------------------------------------------------------------------
1194	!> @brief
1195	!> This subroutine
1196	!>
1197	!> @details
1198	!>
1199	!> @author J.Paul
1200	!> - Nov, 2013- Initial Version
1201	!
1202	!> @param[in]
1203	!> @param[out]
1204	!-------------------------------------------------------------------
1205	!> @code
1206	SUBROUTINE interp__2D( dd_value, dd_fill, &
1207	& id_detect, cd_method, &
1208	& id_rhoi, id_rhoj, &
1209	& ld_even, ld_discont )
1210
1211	IMPLICIT NONE
1212	! Argument
1213	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value
1214	REAL(dp) , INTENT(IN ) :: dd_fill
1215	INTEGER(I4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect
1216	CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
1217	INTEGER(I4) , INTENT(IN ) :: id_rhoi
1218	INTEGER(I4) , INTENT(IN ) :: id_rhoj
1219	LOGICAL , DIMENSION(:) , INTENT(IN ) :: ld_even
1220	LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont
1221
1222	! local variable
1223	INTEGER(I4) :: il_xextra
1224	INTEGER(I4) :: il_yextra
1225	INTEGER(i4), DIMENSION(2) :: il_shape
1226	INTEGER(i4), DIMENSION(2) :: il_dim
1227
1228	REAL(dp) :: dl_min
1229	REAL(dp) :: dl_max
1230	REAL(dp) , DIMENSION(:) , ALLOCATABLE :: dl_coef
1231	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_coarse
1232	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_tmp
1233	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx
1234	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy
1235	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy
1236
1237	LOGICAL :: ll_discont
1238	! loop indices
1239	INTEGER(i4) :: ji
1240	INTEGER(i4) :: jj
1241	INTEGER(i4) :: ii
1242	INTEGER(i4) :: ij
1243
1244	!----------------------------------------------------------------
1245	ll_discont=.FALSE.
1246	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1247
1248	CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
1249	& " discont "//TRIM(fct_str(ll_discont)) )
1250
1251	il_shape(:)=SHAPE(dd_value)
1252
1253	! compute coarse grid dimension
1254	il_xextra=id_rhoi-1
1255	il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
1256
1257	il_yextra=id_rhoj-1
1258	il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj
1259
1260	ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )
1261
1262	! value on coarse grid
1263	dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &
1264	& 1:il_shape(2):id_rhoj )
1265
1266	SELECT CASE(TRIM(cd_method))
1267
1268	CASE('cubic')
1269
1270	ALLOCATE( dl_dfdx( il_dim(1),il_dim(2)), &
1271	& dl_dfdy( il_dim(1),il_dim(2)), &
1272	& dl_d2fdxy(il_dim(1),il_dim(2)) )
1273
1274	! compute derivative on coarse grid
1275	dl_dfdx(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ll_discont)
1276	dl_dfdy(:,:)=extrap_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ll_discont)
1277
1278	! compute cross derivative on coarse grid
1279	dl_d2fdxy(:,:)=extrap_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ll_discont)
1280
1281	ALLOCATE( dl_tmp(2,2) )
1282	ALLOCATE( dl_coef(16) )
1283
1284	DO jj=1,il_shape(2)-1,id_rhoj
1285	ij=((jj-1)/id_rhoj)+1
1286	DO ji=1,il_shape(1)-1,id_rhoi
1287	ii=((ji-1)/id_rhoi)+1
1288
1289	IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
1290	& ANY( dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
1291	& ANY( dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
1292	& ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE
1293
1294	dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
1295	IF( ll_discont )THEN
1296
1297	dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
1298	dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
1299	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1300	WHERE( dl_tmp(:,:) < 0_dp )
1301	dl_tmp(:,:) = dl_tmp(:,:)+360._dp
1302	END WHERE
1303	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1304	WHERE( dl_tmp(:,:) > 180_dp )
1305	dl_tmp(:,:) = dl_tmp(:,:)-180._dp
1306	END WHERE
1307	ENDIF
1308	ENDIF
1309
1310	! compute bicubic coefficient
1311	dl_coef(:)=interp__2D_cubic_coef(dl_tmp(:,:),&
1312	& dl_dfdx( ii:ii+1,ij:ij+1),&
1313	& dl_dfdy( ii:ii+1,ij:ij+1),&
1314	& dl_d2fdxy( ii:ii+1,ij:ij+1),&
1315	& dd_fill )
1316
1317	! compute value on detetected point
1318	CALL interp__2D_cubic_fill(dd_value( ji:ji+id_rhoi, &
1319	& jj:jj+id_rhoj ), &
1320	& id_detect(ji:ji+id_rhoi, &
1321	& jj:jj+id_rhoj ), &
1322	& dl_coef(:), dd_fill, &
1323	& ld_even(:), id_rhoi, id_rhoj )
1324
1325	IF( ll_discont )THEN
1326	WHERE( dd_value( ji:ji+id_rhoi, &
1327	& jj:jj+id_rhoj ) >= 180._dp .AND. &
1328	& dd_value( ji:ji+id_rhoi, &
1329	& jj:jj+id_rhoj ) /= dd_fill )
1330	dd_value( ji:ji+id_rhoi, &
1331	& jj:jj+id_rhoj ) = &
1332	& dd_value( ji:ji+id_rhoi, &
1333	& jj:jj+id_rhoj ) - 360._dp
1334	END WHERE
1335	ENDIF
1336
1337	ENDDO
1338	ENDDO
1339
1340	DEALLOCATE(dl_coef)
1341	DEALLOCATE(dl_tmp )
1342
1343	DEALLOCATE(dl_dfdx, &
1344	& dl_dfdy, &
1345	& dl_d2fdxy )
1346
1347	CASE('nearest')
1348
1349	DO jj=1,il_shape(2)-1,id_rhoj
1350	ij=((jj-1)/id_rhoj)+1
1351	DO ji=1,il_shape(1)-1,id_rhoi
1352	ii=((ji-1)/id_rhoi)+1
1353
1354	! compute value on detetected point
1355	CALL interp__2D_nearest_fill(dd_value( ji:ji+id_rhoi, &
1356	& jj:jj+id_rhoj ), &
1357	& id_detect(ji:ji+id_rhoi, &
1358	& jj:jj+id_rhoj ) )
1359
1360	ENDDO
1361	ENDDO
1362
1363	CASE DEFAULT ! linear
1364
1365	ALLOCATE( dl_coef(4) )
1366	ALLOCATE( dl_tmp(2,2) )
1367
1368	DO jj=1,il_shape(2)-1,id_rhoj
1369	ij=((jj-1)/id_rhoj)+1
1370	DO ji=1,il_shape(1)-1,id_rhoi
1371	ii=((ji-1)/id_rhoi)+1
1372
1373	IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill ) ) CYCLE
1374
1375	dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
1376	IF( ll_discont )THEN
1377
1378	dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
1379	dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
1380	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1381	WHERE( dl_tmp(:,:) < 0_dp )
1382	dl_tmp(:,:) = dl_tmp(:,:)+360._dp
1383	END WHERE
1384	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1385	WHERE( dl_tmp(:,:) > 180_dp )
1386	dl_tmp(:,:) = dl_tmp(:,:)-180._dp
1387	END WHERE
1388	ENDIF
1389	ENDIF
1390
1391	! compute bilinear coefficient
1392	dl_coef(:)=interp__2D_linear_coef(dl_tmp(:,:), dd_fill)
1393
1394	! compute value on detetected point
1395	CALL interp__2D_linear_fill(dd_value( ji:ji+id_rhoi, &
1396	& jj:jj+id_rhoj ), &
1397	& id_detect(ji:ji+id_rhoi, &
1398	& jj:jj+id_rhoj ), &
1399	& dl_coef(:), dd_fill, &
1400	& ld_even(:), id_rhoi, id_rhoj )
1401
1402	IF( ll_discont )THEN
1403	WHERE( dd_value( ji:ji+id_rhoi, &
1404	& jj:jj+id_rhoj ) >= 180._dp .AND. &
1405	& dd_value( ji:ji+id_rhoi, &
1406	& jj:jj+id_rhoj ) /= dd_fill )
1407	dd_value( ji:ji+id_rhoi, &
1408	& jj:jj+id_rhoj ) = &
1409	& dd_value( ji:ji+id_rhoi, &
1410	& jj:jj+id_rhoj ) - 360._dp
1411	END WHERE
1412	ENDIF
1413
1414	ENDDO
1415	ENDDO
1416
1417	DEALLOCATE(dl_coef)
1418	DEALLOCATE(dl_tmp )
1419
1420	END SELECT
1421
1422	DEALLOCATE( dl_coarse )
1423
1424	END SUBROUTINE interp__2D
1425	!> @endcode
1426	!-------------------------------------------------------------------
1427	!> @brief
1428	!> This subroutine
1429	!>
1430	!> @details
1431	!>
1432	!> @author J.Paul
1433	!> - Nov, 2013- Initial Version
1434	!
1435	!> @param[in]
1436	!> @param[out]
1437	!-------------------------------------------------------------------
1438	!> @code
1439	SUBROUTINE interp__1D( dd_value, dd_fill, &
1440	& id_detect, cd_method, &
1441	& id_rhoi, &
1442	& ld_even, ld_discont )
1443	IMPLICIT NONE
1444	! Argument
1445	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value
1446	REAL(dp) , INTENT(IN ) :: dd_fill
1447	INTEGER(I4) , DIMENSION(:), INTENT(INOUT) :: id_detect
1448	CHARACTER(LEN=*) , INTENT(IN ) :: cd_method
1449	INTEGER(I4) , INTENT(IN ) :: id_rhoi
1450	LOGICAL , INTENT(IN ) :: ld_even
1451	LOGICAL , INTENT(IN ), OPTIONAL :: ld_discont
1452
1453	! local variable
1454	INTEGER(i4), DIMENSION(1) :: il_shape
1455	INTEGER(i4), DIMENSION(1) :: il_dim
1456	INTEGER(I4) :: il_xextra
1457
1458	REAL(dp) :: dl_min
1459	REAL(dp) :: dl_max
1460	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coarse
1461	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_tmp
1462	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_dfdx
1463
1464	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_coef
1465
1466	LOGICAL :: ll_discont
1467
1468	! loop indices
1469	INTEGER(i4) :: ji
1470	INTEGER(i4) :: ii
1471	!----------------------------------------------------------------
1472	ll_discont=.FALSE.
1473	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1474
1475	CALL logger_debug("INTERP 1D: interpolation method "//TRIM(cd_method)//&
1476	& " discont "//TRIM(fct_str(ll_discont)) )
1477
1478	il_shape(:)=SHAPE(dd_value)
1479
1480	il_xextra=id_rhoi-1
1481	il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
1482
1483	ALLOCATE( dl_coarse(il_dim(1)) )
1484	! value on coarse grid
1485	dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )
1486
1487	SELECT CASE(TRIM(cd_method))
1488
1489	CASE('cubic')
1490
1491	ALLOCATE( dl_dfdx( il_dim(1)) )
1492
1493	! compute derivative on coarse grid
1494	dl_dfdx(:)=extrap_deriv_1D(dl_coarse(:), dd_fill, ll_discont)
1495
1496	ALLOCATE( dl_coef(4))
1497	ALLOCATE( dl_tmp(2) )
1498
1499	DO ji=1,il_shape(1)-1,id_rhoi
1500	ii=((ji-1)/id_rhoi)+1
1501
1502	IF( ANY( dl_tmp(:)==dd_fill ) .OR. &
1503	& ANY(dl_dfdx(:)==dd_fill ) ) CYCLE
1504
1505	dl_tmp(:)=dl_coarse(ii:ii+1)
1506	IF( ll_discont )THEN
1507
1508	dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
1509	dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
1510	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1511	WHERE( dl_tmp(:) < 0_dp )
1512	dl_tmp(:) = dl_tmp(:)+360._dp
1513	END WHERE
1514	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1515	WHERE( dl_tmp(:) > 180_dp )
1516	dl_tmp(:) = dl_tmp(:)-180._dp
1517	END WHERE
1518	ENDIF
1519	ENDIF
1520
1521	! compute bicubic coefficient
1522	dl_coef(:)=interp__1D_cubic_coef(dl_tmp(:),&
1523	& dl_dfdx( ii:ii+1), &
1524	& dd_fill )
1525
1526	! compute value on detetected point
1527	CALL interp__1D_cubic_fill(dd_value( ji:ji+id_rhoi),&
1528	& id_detect(ji:ji+id_rhoi),&
1529	& dl_coef(:), dd_fill, &
1530	& ld_even, id_rhoi )
1531
1532	IF( ll_discont )THEN
1533	WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
1534	& dd_value( ji:ji+id_rhoi ) /= dd_fill )
1535	dd_value( ji:ji+id_rhoi ) = &
1536	& dd_value( ji:ji+id_rhoi ) - 360._dp
1537	END WHERE
1538	ENDIF
1539
1540	ENDDO
1541
1542	DEALLOCATE(dl_coef)
1543	DEALLOCATE(dl_tmp )
1544
1545	CASE('nearest')
1546
1547	DO ji=1,il_shape(1)-1,id_rhoi
1548	ii=((ji-1)/id_rhoi)+1
1549
1550	! compute value on detetected point
1551	CALL interp__1D_nearest_fill(dd_value( ji:ji+id_rhoi), &
1552	& id_detect(ji:ji+id_rhoi) )
1553
1554	ENDDO
1555
1556	CASE DEFAULT ! linear
1557
1558	ALLOCATE(dl_coef(2))
1559	ALLOCATE( dl_tmp(2) )
1560
1561	DO ji=1,il_shape(1)-1,id_rhoi
1562	ii=((ji-1)/id_rhoi)+1
1563
1564	IF( ANY(dl_coarse(ii:ii+1)==dd_fill ) ) CYCLE
1565
1566	dl_tmp(:)=dl_coarse(ii:ii+1)
1567	IF( ll_discont )THEN
1568
1569	dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
1570	dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
1571	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1572	WHERE( dl_tmp(:) < 0_dp )
1573	dl_tmp(:) = dl_tmp(:)+360._dp
1574	END WHERE
1575	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1576	WHERE( dl_tmp(:) > 180_dp )
1577	dl_tmp(:) = dl_tmp(:)-180._dp
1578	END WHERE
1579	ENDIF
1580	ENDIF
1581
1582	! compute bilinear coefficient
1583	dl_coef(:)=interp__1D_linear_coef(dl_tmp(:), dd_fill)
1584
1585	! compute value on detetected point
1586	CALL interp__1D_linear_fill( dd_value( ji:ji+id_rhoi),&
1587	& id_detect(ji:ji+id_rhoi),&
1588	& dl_coef(:), dd_fill, &
1589	& ld_even, id_rhoi )
1590
1591	IF( ll_discont )THEN
1592	WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
1593	& dd_value( ji:ji+id_rhoi ) /= dd_fill )
1594	dd_value( ji:ji+id_rhoi ) = &
1595	& dd_value( ji:ji+id_rhoi ) - 360._dp
1596	END WHERE
1597	ENDIF
1598
1599	ENDDO
1600
1601	DEALLOCATE(dl_coef)
1602
1603	END SELECT
1604
1605	DEALLOCATE( dl_coarse )
1606
1607	END SUBROUTINE interp__1D
1608	!> @endcode
1609	!-------------------------------------------------------------------
1610	!> @brief
1611	!> This subroutine
1612	!>
1613	!> @details
1614	!>
1615	!> @author J.Paul
1616	!> - Nov, 2013- Initial Version
1617	!
1618	!> @param[in]
1619	!> @param[out]
1620	!-------------------------------------------------------------------
1621	!> @code
1622	FUNCTION interp__2D_cubic_coef( dd_value, &
1623	& dd_dfdx, &
1624	& dd_dfdy, &
1625	& dd_d2fdxy,&
1626	& dd_fill )
1627	IMPLICIT NONE
1628	! Argument
1629	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value
1630	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx
1631	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy
1632	REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy
1633	REAL(dp) , INTENT(IN) :: dd_fill
1634
1635	! function
1636	REAL(dp), DIMENSION(16) :: interp__2D_cubic_coef
1637
1638	! local variable
1639	REAL(dp), DIMENSION(16,16), PARAMETER :: dp_matrix = RESHAPE( &
1640	& (/ 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&
1641	0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&
1642	-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&
1643	2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&
1644	0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,&
1645	0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ,&
1646	0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 ,-2 ,-1 , 0 , 0 ,&
1647	0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 2 ,-2 , 0 , 0 , 1 , 1 , 0 , 0 ,&
1648	-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 , 0 , 0 , 0 , 0 ,&
1649	0 , 0 , 0 , 0 ,-3 , 0 , 3 , 0 , 0 , 0 , 0 , 0 ,-2 , 0 ,-1 , 0 ,&
1650	9 ,-9 ,-9 , 9 , 6 , 3 ,-6 ,-3 , 6 ,-6 , 3 ,-3 , 4 , 2 , 2 , 1 ,&
1651	-6 , 6 , 6 ,-6 ,-3 ,-3 , 3 , 3 ,-4 , 4 ,-2 , 2 ,-2 ,-2 ,-1 ,-1 ,&
1652	2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 ,&
1653	0 , 0 , 0 , 0 , 2 , 0 ,-2 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 0 ,&
1654	-6 , 6 , 6 ,-6 ,-4 ,-2 , 4 , 2 ,-3 , 3 ,-3 , 3 ,-2 ,-1 ,-2 ,-1 ,&
1655	4 ,-4 ,-4 , 4 , 2 , 2 ,-2 ,-2 , 2 ,-2 , 2 ,-2 , 1 , 1 , 1 , 1 /), &
1656	& (/ 16, 16 /) )
1657
1658	REAL(dp), DIMENSION(16) :: dl_vect
1659
1660	!----------------------------------------------------------------
1661	! init
1662	interp__2D_cubic_coef(:)=dd_fill
1663
1664	IF( ANY(SHAPE( dd_value(:,:))/= 2) .OR. &
1665	& ANY(SHAPE( dd_dfdx(:,:))/= 2) .OR. &
1666	& ANY(SHAPE( dd_dfdy(:,:))/= 2) .OR. &
1667	& ANY(SHAPE(dd_d2fdxy(:,:))/= 2) )THEN
1668
1669	CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&
1670	& "input tables. shape should be (/2,2/)")
1671
1672	ELSEIF( ANY( dd_value(:,:) == dd_fill) .OR. &
1673	& ANY( dd_dfdx(:,:) == dd_fill) .OR. &
1674	& ANY( dd_dfdy(:,:) == dd_fill) .OR. &
1675	& ANY(dd_d2fdxy(:,:) == dd_fill) )THEN
1676
1677	CALL logger_warn("INTERP CUBIC COEF: fill value detected. "//&
1678	& "can not compute coefficient ")
1679
1680	ELSE
1681
1682	dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
1683	dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. )
1684	dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. )
1685	dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. )
1686
1687	interp__2D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
1688
1689	ENDIF
1690	END FUNCTION interp__2D_cubic_coef
1691	!> @endcode
1692	!-------------------------------------------------------------------
1693	!> @brief
1694	!> This subroutine
1695	!>
1696	!> @details
1697	!>
1698	!> @author J.Paul
1699	!> - Nov, 2013- Initial Version
1700	!
1701	!> @param[in]
1702	!> @param[out]
1703	!-------------------------------------------------------------------
1704	!> @code
1705	SUBROUTINE interp__2D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &
1706	& ld_even, id_rhoi, id_rhoj )
1707	IMPLICIT NONE
1708	! Argument
1709	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value
1710	INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect
1711	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef
1712	REAL(dp) , INTENT(IN ) :: dd_fill
1713	LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even
1714	INTEGER(I4) , INTENT(IN ) :: id_rhoi
1715	INTEGER(I4) , INTENT(IN ) :: id_rhoj
1716
1717	! local variable
1718	INTEGER(i4), DIMENSION(2) :: il_shape
1719
1720	REAL(dp) , DIMENSION(16) :: dl_vect
1721	REAL(dp) :: dl_dx
1722	REAL(dp) :: dl_dy
1723	REAL(dp) :: dl_x
1724	REAL(dp) :: dl_x2
1725	REAL(dp) :: dl_x3
1726	REAL(dp) :: dl_y
1727	REAL(dp) :: dl_y2
1728	REAL(dp) :: dl_y3
1729
1730	! loop indices
1731	INTEGER(i4) :: ji
1732	INTEGER(i4) :: jj
1733	!----------------------------------------------------------------
1734
1735	IF( SIZE(dd_coef(:)) /= 16 )THEN
1736	CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&
1737	& "coef table. shape should be (/16/)")
1738	ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
1739	CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//&
1740	& "can not compute interpolation.")
1741	ELSE
1742
1743	il_shape(:)=SHAPE(dd_value(:,:))
1744
1745	dl_dx=1./REAL(id_rhoi)
1746	dl_dy=1./REAL(id_rhoj)
1747
1748	DO jj=1,il_shape(2)
1749
1750	IF( ld_even(2) )THEN
1751	dl_y=(jj-1)dl_dy - dl_dy0.5
1752	ELSE ! odd refinement
1753	dl_y=(jj-1)*dl_dy
1754	ENDIF
1755	dl_y2=dl_y*dl_y
1756	dl_y3=dl_y2*dl_y
1757
1758	DO ji=1,il_shape(1)
1759
1760	IF(id_detect(ji,jj)==1)THEN
1761
1762	IF( ld_even(1) )THEN
1763	dl_x=(ji-1)dl_dx - dl_dx0.5
1764	ELSE ! odd refinement
1765	dl_x=(ji-1)*dl_dx
1766	ENDIF
1767	dl_x2=dl_x*dl_x
1768	dl_x3=dl_x2*dl_x
1769
1770	dl_vect(:)=(/1._dp, dl_x , dl_x2 , dl_x3 , &
1771	& dl_y , dl_xdl_y , dl_x2dl_y , dl_x3*dl_y , &
1772	& dl_y2, dl_xdl_y2, dl_x2dl_y2, dl_x3*dl_y2, &
1773	& dl_y3, dl_xdl_y3, dl_x2dl_y3, dl_x3*dl_y3 /)
1774
1775	dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
1776	id_detect(ji,jj)=0
1777
1778	ENDIF
1779
1780	ENDDO
1781	ENDDO
1782
1783	ENDIF
1784
1785	END SUBROUTINE interp__2D_cubic_fill
1786	!> @endcode
1787	!-------------------------------------------------------------------
1788	!> @brief
1789	!> This subroutine
1790	!>
1791	!> @details
1792	!>
1793	!> @author J.Paul
1794	!> - Nov, 2013- Initial Version
1795	!
1796	!> @param[in]
1797	!> @param[out]
1798	!-------------------------------------------------------------------
1799	!> @code
1800	FUNCTION interp__2D_linear_coef( dd_value, dd_fill )
1801	IMPLICIT NONE
1802	! Argument
1803	REAL(dp), DIMENSION(:,:) , INTENT(IN) :: dd_value
1804	REAL(dp) , INTENT(IN) :: dd_fill
1805
1806	! function
1807	REAL(dp), DIMENSION(4) :: interp__2D_linear_coef
1808
1809	! local variable
1810
1811	REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &
1812	& (/ 1 , 0 , 0 , 0 ,&
1813	-1 , 1 , 0 , 0 ,&
1814	-1 , 0 , 1 , 0 ,&
1815	1 ,-1 ,-1 , 1 /), &
1816	& (/ 4, 4 /) )
1817
1818	REAL(dp), DIMENSION(4) :: dl_vect
1819
1820	!----------------------------------------------------------------
1821
1822	IF( ANY(SHAPE(dd_value(:,:))/= 2) )THEN
1823	CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&
1824	& "input tables. shape should be (/2,2/)")
1825	ELSEIF( ANY(dd_value(:,:)==dd_fill) )THEN
1826	CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&
1827	& "can not compute coefficient.")
1828	ELSE
1829
1830	dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
1831
1832	interp__2D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
1833
1834	ENDIF
1835	END FUNCTION interp__2D_linear_coef
1836	!> @endcode
1837	!-------------------------------------------------------------------
1838	!> @brief
1839	!> This subroutine
1840	!>
1841	!> @details
1842	!>
1843	!> @author J.Paul
1844	!> - Nov, 2013- Initial Version
1845	!
1846	!> @param[in]
1847	!> @param[out]
1848	!-------------------------------------------------------------------
1849	!> @code
1850	SUBROUTINE interp__2D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &
1851	& ld_even, id_rhoi, id_rhoj )
1852	IMPLICIT NONE
1853	! Argument
1854	REAL(dp) , DIMENSION(:,:), INTENT(INOUT) :: dd_value
1855	INTEGER(i4) , DIMENSION(:,:), INTENT(INOUT) :: id_detect
1856	REAL(dp) , DIMENSION(:) , INTENT(IN ) :: dd_coef
1857	REAL(dp) , INTENT(IN ) :: dd_fill
1858	LOGICAL , DIMENSION(:), INTENT(IN ) :: ld_even
1859	INTEGER(I4) , INTENT(IN ) :: id_rhoi
1860	INTEGER(I4) , INTENT(IN ) :: id_rhoj
1861
1862	! local variable
1863	INTEGER(i4), DIMENSION(2) :: il_shape
1864
1865	REAL(dp) , DIMENSION(4) :: dl_vect
1866	REAL(dp) :: dl_dx
1867	REAL(dp) :: dl_dy
1868	REAL(dp) :: dl_x
1869	REAL(dp) :: dl_y
1870
1871	! loop indices
1872	INTEGER(i4) :: ji
1873	INTEGER(i4) :: jj
1874	!----------------------------------------------------------------
1875
1876	IF( SIZE(dd_coef(:)) /= 4 )THEN
1877	CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&
1878	& "coef table. shape should be (/4/)")
1879	ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
1880	CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//&
1881	& "can not compute interpolation.")
1882	ELSE
1883
1884	il_shape(:)=SHAPE(dd_value(:,:))
1885
1886	dl_dx=1./REAL(id_rhoi)
1887	dl_dy=1./REAL(id_rhoj)
1888
1889	DO jj=1,il_shape(2)
1890
1891	IF( ld_even(2) )THEN
1892	dl_y=(jj-1)dl_dy - dl_dy0.5
1893	ELSE ! odd refinement
1894	dl_y=(jj-1)*dl_dy
1895	ENDIF
1896
1897	DO ji=1,il_shape(1)
1898
1899	IF(id_detect(ji,jj)==1)THEN
1900
1901	IF( ld_even(1) )THEN
1902	dl_x=(ji-1)dl_dx - dl_dx0.5
1903	ELSE ! odd refinement
1904	dl_x=(ji-1)*dl_dx
1905	ENDIF
1906
1907	dl_vect(:)=(/1._dp, dl_x, dl_y, dl_x*dl_y /)
1908
1909	dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
1910	id_detect(ji,jj)=0
1911
1912	ENDIF
1913
1914	ENDDO
1915	ENDDO
1916
1917	ENDIF
1918
1919	END SUBROUTINE interp__2D_linear_fill
1920	!> @endcode
1921	!-------------------------------------------------------------------
1922	!> @brief
1923	!> This subroutine
1924	!>
1925	!> @details
1926	!>
1927	!> @author J.Paul
1928	!> - Nov, 2013- Initial Version
1929	!
1930	!> @param[in]
1931	!> @param[out]
1932	!-------------------------------------------------------------------
1933	!> @code
1934	SUBROUTINE interp__2D_nearest_fill( dd_value, id_detect )
1935	IMPLICIT NONE
1936	! Argument
1937	REAL(dp) , DIMENSION(:,:) , INTENT(INOUT) :: dd_value
1938	INTEGER(i4), DIMENSION(:,:) , INTENT(INOUT) :: id_detect
1939
1940	! local variable
1941	INTEGER(i4), DIMENSION(2) :: il_shape
1942
1943	INTEGER(i4) :: il_i1
1944	INTEGER(i4) :: il_i2
1945	INTEGER(i4) :: il_j1
1946	INTEGER(i4) :: il_j2
1947
1948	INTEGER(i4) :: il_half1
1949	INTEGER(i4) :: il_half2
1950
1951	! loop indices
1952	INTEGER(i4) :: ji
1953	INTEGER(i4) :: jj
1954	!----------------------------------------------------------------
1955
1956	il_shape(:)=SHAPE(dd_value(:,:))
1957
1958	il_i1=1
1959	il_i2=il_shape(1)
1960
1961	il_j1=1
1962	il_j2=il_shape(2)
1963
1964	il_half1=CEILING(il_shape(1)*0.5)
1965	il_half2=CEILING(il_shape(2)*0.5)
1966
1967	DO jj=1,il_half2
1968
1969	DO ji=1,il_half1
1970
1971	! lower left point
1972	IF(id_detect(ji,jj)==1)THEN
1973
1974	dd_value( ji,jj)=dd_value(il_i1,il_j1)
1975	id_detect(ji,jj)=0
1976
1977	ENDIF
1978
1979	! lower right point
1980	IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN
1981
1982	dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1)
1983	id_detect(il_shape(1)-ji+1,jj)=0
1984
1985	ENDIF
1986
1987	! upper left point
1988	IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN
1989
1990	dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2)
1991	id_detect(ji,il_shape(2)-jj+1)=0
1992
1993	ENDIF
1994
1995	! upper right point
1996	IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN
1997
1998	dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2)
1999	id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=0
2000
2001	ENDIF
2002
2003	ENDDO
2004
2005	ENDDO
2006
2007	END SUBROUTINE interp__2D_nearest_fill
2008	!> @endcode
2009	!-------------------------------------------------------------------
2010	!> @brief
2011	!> This subroutine
2012	!>
2013	!> @details
2014	!>
2015	!> @author J.Paul
2016	!> - Nov, 2013- Initial Version
2017	!
2018	!> @param[in]
2019	!> @param[out]
2020	!-------------------------------------------------------------------
2021	!> @code
2022	FUNCTION interp__1D_cubic_coef( dd_value, &
2023	& dd_dfdx, &
2024	& dd_fill )
2025	IMPLICIT NONE
2026	! Argument
2027	REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value
2028	REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_dfdx
2029	REAL(dp) , INTENT(IN) :: dd_fill
2030
2031	! function
2032	REAL(dp), DIMENSION(4) :: interp__1D_cubic_coef
2033
2034	! local variable
2035
2036	REAL(dp), DIMENSION(4,4), PARAMETER :: dp_matrix = RESHAPE( &
2037	& (/ 1 , 0 , 0 , 0 ,&
2038	-1 , 1 , 0 , 0 ,&
2039	-3 , 3 ,-2 ,-1 ,&
2040	2 ,-2 , 1 , 1 /), &
2041	& (/ 4, 4 /) )
2042
2043	REAL(dp), DIMENSION(4) :: dl_vect
2044
2045	!----------------------------------------------------------------
2046	IF( SIZE(dd_value(:))/= 2 .OR. &
2047	& SIZE(dd_dfdx(:) )/= 2 )THEN
2048
2049	CALL logger_error("INTERP CUBIC COEF: invalid dimension of "//&
2050	& "input tables. shape should be (/2,2/)")
2051
2052	ELSEIF( ANY(dd_value(:)==dd_fill) .OR. &
2053	& ANY( dd_dfdx(:)==dd_fill) )THEN
2054	CALL logger_error("INTERP CUBIC COEF: fill value detected. "//&
2055	& "can not compute coefficient.")
2056	ELSE
2057
2058	dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
2059	dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. )
2060
2061	interp__1D_cubic_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
2062
2063	ENDIF
2064
2065	END FUNCTION interp__1D_cubic_coef
2066	!> @endcode
2067	!-------------------------------------------------------------------
2068	!> @brief
2069	!> This subroutine
2070	!>
2071	!> @details
2072	!>
2073	!> @author J.Paul
2074	!> - Nov, 2013- Initial Version
2075	!
2076	!> @param[in]
2077	!> @param[out]
2078	!-------------------------------------------------------------------
2079	!> @code
2080	SUBROUTINE interp__1D_cubic_fill( dd_value, id_detect, dd_coef, dd_fill, &
2081	& ld_even, id_rhoi )
2082	IMPLICIT NONE
2083	! Argument
2084	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value
2085	INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect
2086	REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef
2087	REAL(dp) , INTENT(IN ) :: dd_fill
2088	LOGICAL , INTENT(IN ) :: ld_even
2089	INTEGER(I4) , INTENT(IN ) :: id_rhoi
2090
2091	! local variable
2092	INTEGER(i4), DIMENSION(1) :: il_shape
2093
2094	REAL(dp) , DIMENSION(4) :: dl_vect
2095	REAL(dp) :: dl_dx
2096	REAL(dp) :: dl_x
2097	REAL(dp) :: dl_x2
2098	REAL(dp) :: dl_x3
2099
2100	! loop indices
2101	INTEGER(i4) :: ji
2102	!----------------------------------------------------------------
2103
2104	IF( SIZE(dd_coef(:)) /= 4 )THEN
2105	CALL logger_error("INTERP CUBIC FILL: invalid dimension of "//&
2106	& "coef table. shape should be (/4/)")
2107	!ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &
2108	ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
2109	CALL logger_error("INTERP CUBIC FILL: fill value detected. "//&
2110	& "can not compute interpolation")
2111	ELSE
2112
2113	il_shape(:)=SHAPE(dd_value(:))
2114
2115	dl_dx=1./REAL(id_rhoi)
2116
2117	DO ji=1,il_shape(1)
2118
2119	IF(id_detect(ji)==1)THEN
2120
2121	IF( ld_even )THEN
2122	dl_x=(ji-1)dl_dx - dl_dx0.5
2123	ELSE ! odd refinement
2124	dl_x=(ji-1)*dl_dx
2125	ENDIF
2126	dl_x2=dl_x*dl_x
2127	dl_x3=dl_x2*dl_x
2128
2129	dl_vect(:)=(/1._dp, dl_x, dl_x2, dl_x3 /)
2130
2131	dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
2132	id_detect(ji)=0
2133
2134	ENDIF
2135
2136	ENDDO
2137
2138	ENDIF
2139
2140	END SUBROUTINE interp__1D_cubic_fill
2141	!> @endcode
2142	!-------------------------------------------------------------------
2143	!> @brief
2144	!> This subroutine
2145	!>
2146	!> @details
2147	!>
2148	!> @author J.Paul
2149	!> - Nov, 2013- Initial Version
2150	!
2151	!> @param[in]
2152	!> @param[out]
2153	!-------------------------------------------------------------------
2154	!> @code
2155	FUNCTION interp__1D_linear_coef( dd_value, dd_fill )
2156	IMPLICIT NONE
2157	! Argument
2158	REAL(dp), DIMENSION(:) , INTENT(IN) :: dd_value
2159	REAL(dp) , INTENT(IN) :: dd_fill
2160
2161	! function
2162	REAL(dp), DIMENSION(2) :: interp__1D_linear_coef
2163
2164	! local variable
2165
2166	REAL(dp), DIMENSION(2,2), PARAMETER :: dp_matrix = RESHAPE( &
2167	& (/ 1 , 0 ,&
2168	-1 , 1 /), &
2169	& (/ 2, 2 /) )
2170
2171	REAL(dp), DIMENSION(2) :: dl_vect
2172
2173	!----------------------------------------------------------------
2174
2175	IF( ANY(SHAPE(dd_value(:))/= 2) )THEN
2176	CALL logger_error("INTERP LINEAR COEF: invalid dimension of "//&
2177	& "input tables. shape should be (/2/)")
2178	ELSEIF( ANY(dd_value(:)==dd_fill) )THEN
2179	CALL logger_error("INTERP LINEAR COEF: fill value detected. "//&
2180	& "can not compute coefficient.")
2181	ELSE
2182
2183	dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
2184
2185	interp__1D_linear_coef(:)=MATMUL(TRANSPOSE(dp_matrix(:,:)),dl_vect(:))
2186
2187	ENDIF
2188
2189	END FUNCTION interp__1D_linear_coef
2190	!> @endcode
2191	!-------------------------------------------------------------------
2192	!> @brief
2193	!> This subroutine
2194	!>
2195	!> @details
2196	!>
2197	!> @author J.Paul
2198	!> - Nov, 2013- Initial Version
2199	!
2200	!> @param[in]
2201	!> @param[out]
2202	!-------------------------------------------------------------------
2203	!> @code
2204	SUBROUTINE interp__1D_linear_fill( dd_value, id_detect, dd_coef, dd_fill, &
2205	& ld_even, id_rhoi )
2206	IMPLICIT NONE
2207	! Argument
2208	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value
2209	INTEGER(i4) , DIMENSION(:), INTENT(INOUT) :: id_detect
2210	REAL(dp) , DIMENSION(:), INTENT(IN ) :: dd_coef
2211	REAL(dp) , INTENT(IN ) :: dd_fill
2212	LOGICAL , INTENT(IN ) :: ld_even
2213	INTEGER(I4) , INTENT(IN ) :: id_rhoi
2214
2215	! local variable
2216	INTEGER(i4), DIMENSION(1) :: il_shape
2217
2218	REAL(dp) , DIMENSION(2) :: dl_vect
2219	REAL(dp) :: dl_dx
2220	REAL(dp) :: dl_x
2221
2222	! loop indices
2223	INTEGER(i4) :: ji
2224	!----------------------------------------------------------------
2225
2226	IF( SIZE(dd_coef(:)) /= 2 )THEN
2227	CALL logger_error("INTERP LINEAR FILL: invalid dimension of "//&
2228	& "coef table. shape should be (/2/)")
2229	!ELSEIF( ANY(dd_value(:)==dd_fill .AND. id_detect(:)==0 ) .OR. &
2230	ELSEIF( ANY( dd_coef(:)==dd_fill ) )THEN
2231	CALL logger_error("INTERP LINEAR FILL: fill value detected. "//&
2232	& "can not compute interpolation")
2233	ELSE
2234
2235	il_shape(:)=SHAPE(dd_value)
2236
2237	dl_dx=1./REAL(id_rhoi)
2238
2239	DO ji=1,il_shape(1)
2240
2241	IF(id_detect(ji)==1)THEN
2242
2243	IF( ld_even )THEN
2244	dl_x=(ji-1)dl_dx - dl_dx0.5
2245	ELSE ! odd refinement
2246	dl_x=(ji-1)*dl_dx
2247	ENDIF
2248
2249	dl_vect(:)=(/1._dp, dl_x /)
2250
2251	dd_value(ji)=DOT_PRODUCT(dd_coef(:),dl_vect(:))
2252	id_detect(ji)=0
2253
2254	ENDIF
2255
2256	ENDDO
2257
2258	ENDIF
2259
2260	END SUBROUTINE interp__1D_linear_fill
2261	!> @endcode
2262	!-------------------------------------------------------------------
2263	!> @brief
2264	!> This subroutine
2265	!>
2266	!> @details
2267	!>
2268	!> @author J.Paul
2269	!> - Nov, 2013- Initial Version
2270	!
2271	!> @param[in]
2272	!> @param[out]
2273	!-------------------------------------------------------------------
2274	!> @code
2275	SUBROUTINE interp__1D_nearest_fill( dd_value, id_detect )
2276	IMPLICIT NONE
2277	! Argument
2278	REAL(dp) , DIMENSION(:), INTENT(INOUT) :: dd_value
2279	INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect
2280
2281	! local variable
2282	INTEGER(i4), DIMENSION(1) :: il_shape
2283
2284	INTEGER(i4) :: il_i1
2285	INTEGER(i4) :: il_i2
2286
2287	INTEGER(i4) :: il_half1
2288
2289	! loop indices
2290	INTEGER(i4) :: ji
2291	!----------------------------------------------------------------
2292
2293	il_shape(:)=SHAPE(dd_value)
2294
2295	il_i1=1
2296	il_i2=il_shape(1)
2297
2298	il_half1=CEILING(il_shape(1)*0.5)
2299
2300	DO ji=1,il_half1
2301
2302	! lower left point
2303	IF(id_detect(ji)==1)THEN
2304
2305	dd_value( ji)=dd_value(il_i1)
2306	id_detect(ji)=0
2307
2308	ENDIF
2309
2310	! lower right point
2311	IF(id_detect(il_shape(1)-ji+1)==1)THEN
2312
2313	dd_value( il_shape(1)-ji+1)=dd_value(il_i2)
2314	id_detect(il_shape(1)-ji+1)=0
2315
2316	ENDIF
2317
2318	ENDDO
2319
2320	END SUBROUTINE interp__1D_nearest_fill
2321	!> @endcode
2322	END MODULE interp

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

Download in other formats: