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.

extrap.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/extrap.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: 68.9 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! MODULE: extrap
6	!
7	! DESCRIPTION:
8	!> @brief
9	!> This module
10	!>
11	!> @details
12	!>
13	!> @author
14	!> J.Paul
15	! REVISION HISTORY:
16	!> @date Nov, 2013 - Initial Version
17	!>
18	!> @note WARNING: FillValue must not be zero (use var_chg_FillValue)
19	!>
20	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
21	!> @todo
22	!> - something wrong when computing point to be extralopated
23	!> - take care of ew value in variable structure
24	!----------------------------------------------------------------------
25	MODULE extrap
26	USE netcdf ! nf90 library
27	USE kind
28	USE phycst
29	USE global
30	USE fct
31	USE logger
32	USE dim
33	USE att
34	USE var
35
36	IMPLICIT NONE
37	PRIVATE
38	! NOTE_avoid_public_variables_if_possible
39
40	! type and variable
41
42	! function and subroutine
43	PUBLIC :: extrap_detect !< detected point to be extrapolated
44	PUBLIC :: extrap_fill_value !< extrapolate value over detected point
45	PUBLIC :: extrap_add_extrabands !<
46	PUBLIC :: extrap_del_extrabands !<
47	PUBLIC :: extrap_deriv_1D !<
48	PUBLIC :: extrap_deriv_2D !<
49
50	PRIVATE :: extrap__detect_wrapper !< detected point to be extrapolated
51	PRIVATE :: extrap__detect !< detected point to be extrapolated
52	PRIVATE :: extrap__fill_value_wrapper !< extrapolate value over detected point
53	PRIVATE :: extrap__fill_value !< extrapolate value over detected point
54	PRIVATE :: extrap__3D
55	PRIVATE :: extrap_deriv_3D
56	PRIVATE :: extrap__3D_min_error_coef
57	PRIVATE :: extrap__3D_min_error_fill
58	PRIVATE :: extrap__3D_dist_weight_coef
59	PRIVATE :: extrap__3D_dist_weight_fill
60
61	INTEGER(i4), PARAMETER :: im_maxiter = 10 !< default maximum number of iteration
62	INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate
63	INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation
64
65	INTERFACE extrap_detect
66	MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated
67	END INTERFACE extrap_detect
68
69	INTERFACE extrap_fill_value
70	MODULE PROCEDURE extrap__fill_value_wrapper !< detected point to be interpolated
71	END INTERFACE extrap_fill_value
72
73	CONTAINS
74	!-------------------------------------------------------------------
75	!> @brief
76	!> This function detected point to be extrapolated.
77	!>
78	!> @details
79	!>
80	!> @author J.Paul
81	!> - Nov, 2013- Initial Version
82	!
83	!> @param[in] td_var : variable to extrapolate
84	!> @param[in] id_iext : number of points to be extrapolated in i-direction
85	!> @param[in] id_jext : number of points to be extrapolated in j-direction
86	!> @param[in] id_kext : number of points to be extrapolated in k-direction
87	!> @return
88	!-------------------------------------------------------------------
89	!> @code
90	FUNCTION extrap__detect( td_var0, td_level1, &
91	& id_offset, id_rho, id_ext )
92	IMPLICIT NONE
93	! Argument
94	TYPE(TVAR) , INTENT(INOUT) :: td_var0
95	TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level1
96	INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
97	INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
98	INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext
99
100	! function
101	INTEGER(i4), DIMENSION(td_var0%t_dim(1)%i_len,&
102	& td_var0%t_dim(2)%i_len,&
103	& td_var0%t_dim(3)%i_len ) :: extrap__detect
104
105	! local variable
106	CHARACTER(LEN=lc) :: cl_level
107
108	INTEGER(i4) :: il_varid
109	INTEGER(i4) , DIMENSION(:,:,:), ALLOCATABLE :: il_detect
110	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_offset
111	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1
112	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_level1_G0
113	INTEGER(i4) , DIMENSION(:,:) , ALLOCATABLE :: il_extra
114	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_ext
115	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_rho
116	INTEGER(i4) , DIMENSION(:) , ALLOCATABLE :: il_dim0
117
118	TYPE(TVAR) :: tl_var1
119
120	! loop indices
121	INTEGER(i4) :: ji0
122	INTEGER(i4) :: jj0
123	INTEGER(i4) :: jk0
124	INTEGER(i4) :: ji1
125	INTEGER(i4) :: jj1
126	INTEGER(i4) :: ji1m
127	INTEGER(i4) :: jj1m
128	INTEGER(i4) :: ji1p
129	INTEGER(i4) :: jj1p
130	!----------------------------------------------------------------
131
132	! init
133	extrap__detect(:,:,:)=0
134
135	ALLOCATE( il_dim0(3) )
136	il_dim0(:)=td_var0%t_dim(1:3)%i_len
137
138	! optional argument
139	ALLOCATE( il_rho(ip_maxdim) )
140	il_rho(:)=1
141	IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
142
143	ALLOCATE( il_offset(ip_maxdim,2) )
144	il_offset(:,:)=0
145	il_offset(jp_I,:)=FLOOR(REAL(il_rho(jp_I)-1,dp)*0.5)
146	il_offset(jp_J,:)=FLOOR(REAL(il_rho(jp_J)-1,dp)*0.5)
147	IF( PRESENT(id_offset) )THEN
148	il_offset(1:SIZE(id_offset(:,:),DIM=1),&
149	& 1:SIZE(id_offset(:,:),DIM=2) )= id_offset(:,:)
150	ENDIF
151
152	ALLOCATE( il_ext(ip_maxdim) )
153	il_ext(:)=im_minext
154	IF( PRESENT(id_ext) ) il_ext(1:SIZE(id_ext(:)))=id_ext(:)
155
156	ALLOCATE( il_detect(il_dim0(1),&
157	& il_dim0(2),&
158	& il_dim0(3)) )
159	il_detect(:,:,:)=0
160
161	! select point already inform
162	WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=1
163
164	IF( PRESENT(td_level1) )THEN
165	SELECT CASE(TRIM(td_var0%c_point))
166	CASE DEFAULT !'T'
167	cl_level='tlevel'
168	CASE('U')
169	cl_level='ulevel'
170	CASE('V')
171	cl_level='vlevel'
172	CASE('F')
173	cl_level='flevel'
174	END SELECT
175	il_varid=var_get_id(td_level1(:),TRIM(cl_level))
176	IF( il_varid == 0 )THEN
177	CALL logger_error("EXTRAP DETECT: can not compute point to be "//&
178	& "extrapolated for variable "//TRIM(td_var0%c_name)//&
179	& ". can not find "//&
180	& "level for variable point "//TRIM(TRIM(td_var0%c_point)))
181	ELSE
182	print *,'read ',trim(cl_level)
183	tl_var1=td_level1(il_varid)
184	ENDIF
185
186	ALLOCATE( il_level1_G0( il_dim0(1), il_dim0(2)) )
187	IF( ALL(tl_var1%t_dim(1:2)%i_len == il_dim0(1:2)) )THEN
188
189	! variable to be extrapolated use same resolution than level
190	il_level1_G0(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)
191
192	ELSE
193	! variable to be extrapolated do not use same resolution than level
194	ALLOCATE( il_level1(tl_var1%t_dim(1)%i_len, &
195	& tl_var1%t_dim(2)%i_len) )
196	! match fine grid vertical level with coarse grid
197	il_level1(:,:)=INT(tl_var1%d_value(:,:,1,1),i4)/il_rho(jp_K)
198
199	ALLOCATE( il_extra(ig_ndim,2) )
200	! coarsening fine grid level
201	il_extra(jp_I,1)=CEILING(REAL(il_rho(jp_I)-1,dp)/2._dp)
202	il_extra(jp_I,2)=FLOOR(REAL(il_rho(jp_I)-1,dp)/2._dp)
203
204	il_extra(jp_J,1)=CEILING(REAL(il_rho(jp_J)-1,dp)/2._dp)
205	il_extra(jp_J,2)=FLOOR(REAL(il_rho(jp_J)-1,dp)/2._dp)
206
207	DO jj0=1,td_var0%t_dim(2)%i_len
208
209	jj1=(jj0-1)*il_rho(jp_J)+1-il_offset(jp_J,1)
210
211	jj1m=MAX( jj1-il_extra(jp_J,1), 1 )
212	jj1p=MIN( jj1+il_extra(jp_J,2), &
213	& tl_var1%t_dim(2)%i_len-il_offset(jp_J,2) )
214
215	DO ji0=1,td_var0%t_dim(1)%i_len
216
217	ji1=(ji0-1)*il_rho(jp_I)+1-id_offset(jp_I,1)
218
219	ji1m=MAX( ji1-il_extra(jp_I,1), 1 )
220	ji1p=MIN( ji1+il_extra(jp_I,2), &
221	& tl_var1%t_dim(1)%i_len-id_offset(jp_I,2) )
222
223	il_level1_G0(ji0,jj0)=MAXVAL(il_level1(ji1m:ji1p,jj1m:jj1p))
224	ENDDO
225	ENDDO
226
227	!il_level1_G0(:,:)=0
228	!DO jj1=1,tl_var1%t_dim(2)%i_len
229
230	! jj0=INT(REAL((jj1+il_offset(jp_J,1)-1)-1,dp)/REAL(il_rho(jp_J),dp)) +1
231
232	! DO ji1=1,tl_var1%t_dim(1)%i_len
233
234	! ji0=INT(REAL((ji1+il_offset(jp_I,1)-1)-1,dp)/REAL(il_rho(jp_I),dp)) +1
235
236	! il_level1_G0(ji0,jj0)=MAX(il_level1_G0(ji0,jj0),il_level1(ji1,jj1))
237	!
238	! ENDDO
239	!ENDDO
240
241	! clean
242	DEALLOCATE( il_extra )
243	DEALLOCATE( il_level1 )
244
245	ENDIF
246
247	! look for sea point
248	!il_detect(:,:,1)=0
249	DO jk0=1,td_var0%t_dim(3)%i_len
250	WHERE( il_level1_G0(:,:) >= jk0)
251	il_detect(:,:,jk0)=1
252	END WHERE
253	!il_detect(:,:,jk0)=il_level1_G0(:,:)
254	!WHERE( td_var0%d_value(:,:,jk0,1) /= td_var0%d_fill )
255	! il_detect(:,:,1)=jk0-1
256	!END WHERE
257	ENDDO
258
259	! clean
260	DEALLOCATE( il_level1_G0 )
261
262	ENDIF
263
264	! clean
265	DEALLOCATE( il_offset )
266
267	!! select extra point depending on interpolation method
268	!! compute point near grid point already inform
269	!FORALL( ji0=1:il_dim0(1), &
270	!& jj0=1:il_dim0(2), &
271	!& jk0=1:il_dim0(3), &
272	!& il_detect(ji0,jj0,jk0) == 1 )
273
274	! il_detect(MAX(1,ji0-il_ext(jp_I)):MIN(ji0+il_ext(jp_I),il_dim0(1)),&
275	! & MAX(1,jj0-il_ext(jp_J)):MIN(jj0+il_ext(jp_J),il_dim0(2)),&
276	! & MAX(1,jk0-il_ext(jp_K)):MIN(jk0+il_ext(jp_K),il_dim0(3)) )=1
277
278	!END FORALL
279
280	! do not compute grid point already inform
281	WHERE( td_var0%d_value(:,:,:,1) /= td_var0%d_fill ) il_detect(:,:,:)=0
282
283	! save result
284	extrap__detect(:,:,:)=il_detect(:,:,:)
285
286	! clean
287	DEALLOCATE( il_dim0 )
288	DEALLOCATE( il_ext )
289	DEALLOCATE( il_detect )
290
291	END FUNCTION extrap__detect
292	!> @endcode
293	!-------------------------------------------------------------------
294	!> @brief
295	!> This function detected point to be extrapolated.
296	!>
297	!> @details
298	!>
299	!> @author J.Paul
300	!> - Nov, 2013- Initial Version
301	!
302	!> @param[in] td_var : variable to extrapolate
303	!> @param[in] id_iext : number of points to be extrapolated in i-direction
304	!> @param[in] id_jext : number of points to be extrapolated in j-direction
305	!> @param[in] id_kext : number of points to be extrapolated in k-direction
306	!> @return
307	!-------------------------------------------------------------------
308	!> @code
309	FUNCTION extrap__detect_wrapper( td_var, td_level, &
310	& id_offset, id_rho, id_ext )
311
312	IMPLICIT NONE
313	! Argument
314	TYPE(TVAR) , INTENT(INOUT) :: td_var
315	TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level
316	INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
317	INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
318	INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_ext
319
320	! function
321	INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len,&
322	& td_var%t_dim(2)%i_len,&
323	& td_var%t_dim(3)%i_len ) :: extrap__detect_wrapper
324
325	! local variable
326	! loop indices
327	!----------------------------------------------------------------
328	! init
329	extrap__detect_wrapper(:,:,:)=0
330
331	IF( .NOT. ANY(td_var%t_dim(1:3)%l_use) )THEN
332	! no dimension I-J-K used
333	CALL logger_debug(" EXTRAP DETECT: nothing done for variable"//&
334	& TRIM(td_var%c_name) )
335	ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN
336
337	! detect point to be interpolated on I-J-K
338	CALL logger_debug(" EXTRAP DETECT: detect point "//&
339	& " for variable "//TRIM(td_var%c_name) )
340
341	extrap__detect_wrapper(:,:,:)=extrap__detect( td_var, td_level, &
342	& id_offset, &
343	& id_rho, &
344	& id_ext )
345
346	ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
347
348	! detect point to be interpolated on I-J
349	CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//&
350	& " for variable "//TRIM(td_var%c_name) )
351
352	extrap__detect_wrapper(:,:,1:1)=extrap__detect( td_var , td_level,&
353	& id_offset, &
354	& id_rho, &
355	& id_ext )
356
357	ELSE IF( td_var%t_dim(3)%l_use )THEN
358
359	! detect point to be interpolated on K
360	CALL logger_debug(" EXTRAP DETECT: detect vertical point "//&
361	& " for variable "//TRIM(td_var%c_name) )
362
363	extrap__detect_wrapper(1:1,1:1,:)=extrap__detect( td_var , td_level, &
364	& id_offset, &
365	& id_rho, &
366	& id_ext )
367
368	ENDIF
369
370	CALL logger_debug(" EXTRAP DETECT: "//&
371	& TRIM(fct_str(SUM(extrap__detect_wrapper(:,:,:))))//&
372	& " points to be extrapolated" )
373
374	END FUNCTION extrap__detect_wrapper
375	!> @endcode
376	! !-------------------------------------------------------------------
377	! !> @brief
378	! !> This function detected point to be extrapolated.
379	! !>
380	! !> @details
381	! !>
382	! !> @author J.Paul
383	! !> - Nov, 2013- Initial Version
384	! !
385	! !> @param[in] td_var : variable to extrapolate
386	! !> @param[in] id_iext : number of points to be extrapolated in i-direction
387	! !> @param[in] id_jext : number of points to be extrapolated in j-direction
388	! !> @param[in] id_kext : number of points to be extrapolated in k-direction
389	! !> @return
390	! !
391	! !> @todo
392	! !-------------------------------------------------------------------
393	! !> @code
394	! FUNCTION extrap__detect(td_var, id_iext, id_jext, id_kext)
395	! IMPLICIT NONE
396	! ! Argument
397	! TYPE(TVAR) , INTENT(INOUT) :: td_var
398	! !INTEGER(i4), DIMENSION(:,:,:), INTENT(OUT ) :: id_detect
399	! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext
400	! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext
401	! INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext
402	!
403	! ! function
404	! INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len, &
405	! & td_var%t_dim(2)%i_len, &
406	! & td_var%t_dim(3)%i_len) :: extrap__detect
407	!
408	! ! local variable
409	! INTEGER(i4) :: il_iext
410	! INTEGER(i4) :: il_jext
411	! INTEGER(i4) :: il_kext
412	!
413	! INTEGER(i4), DIMENSION(3) :: il_dim
414	!
415	! ! loop indices
416	! INTEGER(i4) :: ji
417	! INTEGER(i4) :: jj
418	! INTEGER(i4) :: jk
419	! !----------------------------------------------------------------
420	!
421	! ! optional argument
422	! il_iext=im_minext
423	! IF( PRESENT(id_iext) ) il_iext=id_iext
424	! il_jext=im_minext
425	! IF( PRESENT(id_jext) ) il_jext=id_jext
426	! il_kext=im_minext
427	! IF( PRESENT(id_kext) ) il_kext=id_kext
428	!
429	! ! init
430	! extrap__detect(:,:,:)=0
431	!
432	! il_dim(:)=td_var%t_dim(1:3)%i_len
433	!
434	! ! compute point near grid point already inform
435	! FORALL( ji=1:il_dim(1), &
436	! & jj=1:il_dim(2), &
437	! & jk=1:il_dim(3), &
438	! & td_var%d_value(ji,jj,jk,1) /= td_var%d_fill )
439	!
440	! extrap__detect(MAX(1,ji-il_iext):MIN(ji+il_iext,il_dim(1)),&
441	! & MAX(1,jj-il_jext):MIN(jj+il_jext,il_dim(2)),&
442	! & MAX(1,jk-il_kext):MIN(jk+il_kext,il_dim(3)) )=1
443	!
444	!
445	! END FORALL
446	!
447	! ! do not compute grid point already inform
448	! WHERE( td_var%d_value(:,:,:,1) /= td_var%d_fill ) extrap__detect(:,:,:)=0
449	!
450	! END FUNCTION extrap__detect
451	! !> @endcode
452	!-------------------------------------------------------------------
453	!> @brief
454	!> This subroutine
455	!>
456	!> @details
457	!>
458	!> @author J.Paul
459	!> - Nov, 2013- Initial Version
460	!
461	!> @param[inout] td_var : variable structure
462	!> @param[in] id_iext : number of points to be extrapolated in i-direction
463	!> @param[in] id_jext : number of points to be extrapolated in j-direction
464	!> @param[in] id_kext : number of points to be extrapolated in k-direction
465	!> @param[in] id_extend : radius of the box used to compute extrapolation
466	!> @param[in] id_maxiter : maximum nuber of iteration
467	!-------------------------------------------------------------------
468	!> @code
469	SUBROUTINE extrap__fill_value_wrapper( td_var, td_level, &
470	& id_offset, &
471	& id_rho, &
472	& id_iext, id_jext, id_kext, &
473	& id_radius, id_maxiter )
474	IMPLICIT NONE
475	! Argument
476	TYPE(TVAR) , INTENT(INOUT) :: td_var
477	TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level
478	INTEGER(i4), DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
479	INTEGER(i4), DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
480	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_iext
481	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jext
482	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_kext
483	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius
484	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_maxiter
485
486	! local variable
487	INTEGER(i4) :: il_iext
488	INTEGER(i4) :: il_jext
489	INTEGER(i4) :: il_kext
490	INTEGER(i4) :: il_radius
491	INTEGER(i4) :: il_maxiter
492
493	CHARACTER(LEN=lc) :: cl_method
494
495	! loop indices
496	!----------------------------------------------------------------
497	IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
498	CALL logger_error("EXTRAP FILL VALUE: no table of value "//&
499	& "associted to variable "//TRIM(td_var%c_name) )
500	ELSE
501
502	SELECT CASE(TRIM(td_var%c_extrap(1)))
503	CASE('min_error')
504	cl_method='min_error'
505	CASE DEFAULT
506	cl_method='dist_weight'
507
508	!update variable structure
509	td_var%c_extrap(1)='dist_weight'
510	END SELECT
511
512	il_iext=im_minext
513	IF( PRESENT(id_iext) ) il_iext=id_iext
514	il_jext=im_minext
515	IF( PRESENT(id_jext) ) il_jext=id_jext
516	il_kext=0
517	IF( PRESENT(id_kext) ) il_kext=id_kext
518
519	IF( TRIM(td_var%c_interp(1)) == 'cubic')THEN
520	IF( il_iext > 0 .AND. il_iext < im_mincubic ) il_iext=im_mincubic
521	IF( il_jext > 0 .AND. il_jext < im_mincubic ) il_jext=im_mincubic
522	ENDIF
523
524	IF( il_iext < 0 )THEN
525	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
526	& " number of points to be extrapolated in i-direction "//&
527	& "("//TRIM(fct_str(il_iext))//")")
528	ENDIF
529
530	IF( il_jext < 0 )THEN
531	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
532	& " number of points to be extrapolated in j-direction "//&
533	& "("//TRIM(fct_str(il_jext))//")")
534	ENDIF
535
536	IF( il_kext < 0 )THEN
537	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
538	& " number of points to be extrapolated in k-direction "//&
539	& "("//TRIM(fct_str(il_kext))//")")
540	ENDIF
541
542	IF( (il_iext /= 0 .AND. td_var%t_dim(1)%l_use) .OR. &
543	& (il_jext /= 0 .AND. td_var%t_dim(2)%l_use) .OR. &
544	& (il_kext /= 0 .AND. td_var%t_dim(3)%l_use) .OR. &
545	& PRESENT(td_level) )THEN
546
547	! number of point use to compute box
548	il_radius=1
549	IF( PRESENT(id_radius) ) il_radius=id_radius
550	IF( il_radius < 0 )THEN
551	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
552	& " radius of the box used to compute extrapolation "//&
553	& "("//TRIM(fct_str(il_radius))//")")
554	ENDIF
555
556	! maximum number of iteration
557	il_maxiter=im_maxiter
558	IF( PRESENT(id_maxiter) ) il_maxiter=id_maxiter
559	IF( il_maxiter < 0 )THEN
560	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
561	& " maximum nuber of iteration "//&
562	& "("//TRIM(fct_str(il_maxiter))//")")
563	ENDIF
564
565	CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//&
566	& " using "//TRIM(cl_method)//" method." )
567
568	CALL extrap__fill_value( td_var, cl_method, &
569	& il_iext, il_jext, il_kext, &
570	& il_radius, il_maxiter, &
571	& td_level, &
572	& id_offset, id_rho )
573
574	ENDIF
575
576	ENDIF
577
578	END SUBROUTINE extrap__fill_value_wrapper
579	!> @endcode
580	!-------------------------------------------------------------------
581	!> @brief
582	!> This subroutine
583	!>
584	!> @details
585	!>
586	!> @author J.Paul
587	!> - Nov, 2013- Initial Version
588	!
589	!> @param[inout] td_var : variable structure
590	!> @param[in] cd_method : extrapolation method
591	!> @param[in] id_iext : number of points to be extrapolated in i-direction
592	!> @param[in] id_jext : number of points to be extrapolated in j-direction
593	!> @param[in] id_kext : number of points to be extrapolated in k-direction
594	!> @param[in] id_radius : radius of the halo used to compute extrapolation
595	!> @param[in] id_maxiter : maximum nuber of iteration
596	!-------------------------------------------------------------------
597	!> @code
598	SUBROUTINE extrap__fill_value( td_var, cd_method, &
599	& id_iext, id_jext, id_kext, &
600	& id_radius, id_maxiter, &
601	& td_level, &
602	& id_offset, &
603	& id_rho )
604	IMPLICIT NONE
605	! Argument
606	TYPE(TVAR) , INTENT(INOUT) :: td_var
607	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
608	INTEGER(i4) , INTENT(IN ) :: id_iext
609	INTEGER(i4) , INTENT(IN ) :: id_jext
610	INTEGER(i4) , INTENT(IN ) :: id_kext
611	INTEGER(i4) , INTENT(IN ) :: id_radius
612	INTEGER(i4) , INTENT(IN ) :: id_maxiter
613	TYPE(TVAR) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: td_level
614	INTEGER(i4) , DIMENSION(:,:), INTENT(IN ), OPTIONAL :: id_offset
615	INTEGER(i4) , DIMENSION(:) , INTENT(IN ), OPTIONAL :: id_rho
616
617	! local variable
618	CHARACTER(LEN=lc) :: cl_extrap
619
620	INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
621	INTEGER(i4) :: il_radius
622	INTEGER(i4) :: il_iter
623
624	TYPE(TATT) :: tl_att
625
626	! loop indices
627	INTEGER(i4) :: jl
628	!----------------------------------------------------------------
629
630	!1- detect point to be extrapolated
631	ALLOCATE( il_detect( td_var%t_dim(1)%i_len, &
632	& td_var%t_dim(2)%i_len, &
633	& td_var%t_dim(3)%i_len) )
634
635	il_detect(:,:,:) = extrap_detect( td_var, td_level, &
636	& id_offset, &
637	& id_rho, &
638	& id_ext=(/id_iext, id_jext, id_kext/) )
639
640	!2- add attribute to variable
641	cl_extrap=fct_concat(td_var%c_extrap(:))
642	tl_att=att_init('extrapolation',cl_extrap)
643	CALL var_move_att(td_var, tl_att)
644
645	CALL logger_warn(" EXTRAP FILL: "//&
646	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
647	& " point(s) to extrapolate " )
648
649	!3- extrapolate
650	DO jl=1,td_var%t_dim(4)%i_len
651
652	! td_var%d_value(:,:,:,jl)=il_detect(:,:,:)
653	il_iter=1
654	DO WHILE( ANY(il_detect(:,:,:)==1) )
655	! change extend value to minimize number of iteration
656	il_radius=id_radius+(il_iter/id_maxiter)
657
658	CALL logger_debug(" EXTRAP FILL VALUE: "//&
659	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
660	& " points to extrapolate " )
661
662	CALL extrap__3D(td_var%d_value(:,:,:,jl), td_var%d_fill, &
663	& il_detect(:,:,:), &
664	& cd_method, il_radius )
665
666	il_iter=il_iter+1
667	ENDDO
668
669	ENDDO
670
671	IF( SUM(il_detect(:,:,:)) /= 0 )THEN
672	CALL logger_warn(" EXTRAP FILL: still "//&
673	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
674	& " point(s) to extrapolate " )
675	ENDIF
676
677	DEALLOCATE(il_detect)
678
679	END SUBROUTINE extrap__fill_value
680	!> @endcode
681	!-------------------------------------------------------------------
682	!> @brief
683	!> This subroutine
684	!>
685	!> @details
686	!>
687	!> @author J.Paul
688	!> - Nov, 2013- Initial Version
689	!
690	!> @param[inout] dd_value : 3D table of variable to be extrapolated
691	!> @param[in] dd_fill : FillValue of variable
692	!> @param[inout] id_detect : table of point to extrapolate
693	!> @param[in] id_ext : number of point use to compute box
694	!> @param[in] cd_method : extrapolation method
695	!-------------------------------------------------------------------
696	!> @code
697	SUBROUTINE extrap__3D( dd_value, dd_fill, id_detect,&
698	& cd_method, id_ext )
699	IMPLICIT NONE
700	! Argument
701	REAL(dp) , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value
702	REAL(dp) , INTENT(IN ) :: dd_fill
703	INTEGER(i4), DIMENSION(:,:,:), INTENT(INOUT) :: id_detect
704	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
705	INTEGER(i4), INTENT(IN ) :: id_ext
706
707	! local variable
708	INTEGER(i4) :: il_imin
709	INTEGER(i4) :: il_imax
710	INTEGER(i4) :: il_jmin
711	INTEGER(i4) :: il_jmax
712	INTEGER(i4) :: il_kmin
713	INTEGER(i4) :: il_kmax
714
715	INTEGER(i4), DIMENSION(3) :: il_shape
716	INTEGER(i4), DIMENSION(3) :: il_dim
717
718	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
719	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy
720	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz
721	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
722
723	! loop indices
724	INTEGER(i4) :: ji
725	INTEGER(i4) :: jj
726	INTEGER(i4) :: jk
727	!----------------------------------------------------------------
728
729	il_shape(:)=SHAPE(dd_value)
730
731	SELECT CASE(TRIM(cd_method))
732	CASE('min_error')
733
734	ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )
735	ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) )
736	ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) )
737
738
739	! compute derivative in i-direction
740	dl_dfdx(:,:,:)=dd_fill
741	IF( il_shape(1) > 1 )THEN
742	dl_dfdx(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'I' )
743	ENDIF
744
745	! compute derivative in i-direction
746	dl_dfdy(:,:,:)=dd_fill
747	IF( il_shape(2) > 1 )THEN
748	dl_dfdy(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'J' )
749	ENDIF
750
751	! compute derivative in i-direction
752	dl_dfdz(:,:,:)=dd_fill
753	IF( il_shape(3) > 1 )THEN
754	dl_dfdz(:,:,:)=extrap_deriv_3D( dd_value(:,:,:), dd_fill, 'K' )
755	ENDIF
756
757	il_dim(1)=2*id_ext+1
758	IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1
759	il_dim(2)=2*id_ext+1
760	IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1
761	il_dim(3)=2*id_ext+1
762	IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1
763
764	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
765
766	dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), &
767	& 1:il_dim(2), &
768	& 1:il_dim(3)))
769
770	DO jk=1,il_shape(3)
771	IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE
772	DO jj=1,il_shape(2)
773	IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE
774	DO ji=1,il_shape(1)
775
776	IF( id_detect(ji,jj,jk) == 1 )THEN
777
778	il_imin=MAX(ji-id_ext,1)
779	il_imax=MIN(ji+id_ext,il_shape(1))
780	IF( il_dim(1) == 1 )THEN
781	il_imin=ji
782	il_imax=ji
783	ENDIF
784
785	il_jmin=MAX(jj-id_ext,1)
786	il_jmax=MIN(jj+id_ext,il_shape(2))
787	IF( il_dim(2) == 1 )THEN
788	il_jmin=jj
789	il_jmax=jj
790	ENDIF
791
792	il_kmin=MAX(jk-id_ext,1)
793	il_kmax=MIN(jk+id_ext,il_shape(3))
794	IF( il_dim(3) == 1 )THEN
795	il_kmin=jk
796	il_kmax=jk
797	ENDIF
798
799	dd_value(ji,jj,jk)=extrap__3D_min_error_fill( &
800	& dd_value( il_imin:il_imax, &
801	& il_jmin:il_jmax, &
802	& il_kmin:il_kmax ), dd_fill, id_ext, &
803	& dl_dfdx( il_imin:il_imax, &
804	& il_jmin:il_jmax, &
805	& il_kmin:il_kmax ), &
806	& dl_dfdy( il_imin:il_imax, &
807	& il_jmin:il_jmax, &
808	& il_kmin:il_kmax ), &
809	& dl_dfdz( il_imin:il_imax, &
810	& il_jmin:il_jmax, &
811	& il_kmin:il_kmax ), &
812	& dl_coef(:,:,:) )
813
814	IF( dd_value(ji,jj,jk) /= dd_fill )THEN
815	id_detect(ji,jj,jk)= 0
816	ENDIF
817
818	ENDIF
819
820	ENDDO
821	ENDDO
822	ENDDO
823
824	IF( ALLOCATED(dl_dfdx) ) DEALLOCATE( dl_dfdx )
825	IF( ALLOCATED(dl_dfdy) ) DEALLOCATE( dl_dfdy )
826	IF( ALLOCATED(dl_dfdz) ) DEALLOCATE( dl_dfdz )
827	IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef )
828
829	CASE DEFAULT ! 'dist_weight'
830
831	il_dim(1)=2*id_ext+1
832	IF( il_shape(1) < 2*id_ext+1 ) il_dim(1)=1
833	il_dim(2)=2*id_ext+1
834	IF( il_shape(2) < 2*id_ext+1 ) il_dim(2)=1
835	il_dim(3)=2*id_ext+1
836	IF( il_shape(3) < 2*id_ext+1 ) il_dim(3)=1
837
838	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
839
840	dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1), &
841	& 1:il_dim(2), &
842	& 1:il_dim(3)) )
843
844	DO jk=1,il_shape(3)
845	IF( ALL(id_detect(:,:,jk) == 0) ) CYCLE
846	DO jj=1,il_shape(2)
847	IF( ALL(id_detect(:,jj,jk) == 0) ) CYCLE
848	DO ji=1,il_shape(1)
849
850	IF( id_detect(ji,jj,jk) == 1 )THEN
851
852	il_imin=MAX(ji-id_ext,1)
853	il_imax=MIN(ji+id_ext,il_shape(1))
854	IF( il_dim(1) == 1 )THEN
855	il_imin=ji
856	il_imax=ji
857	ENDIF
858
859	il_jmin=MAX(jj-id_ext,1)
860	il_jmax=MIN(jj+id_ext,il_shape(2))
861	IF( il_dim(2) == 1 )THEN
862	il_jmin=jj
863	il_jmax=jj
864	ENDIF
865
866	il_kmin=MAX(jk-id_ext,1)
867	il_kmax=MIN(jk+id_ext,il_shape(3))
868	IF( il_dim(3) == 1 )THEN
869	il_kmin=jk
870	il_kmax=jk
871	ENDIF
872
873	dd_value(ji,jj,jk)=extrap__3D_dist_weight_fill( &
874	& dd_value( il_imin:il_imax, &
875	& il_jmin:il_jmax, &
876	& il_kmin:il_kmax ), dd_fill, id_ext, &
877	& dl_coef(:,:,:) )
878
879	IF( dd_value(ji,jj,jk) /= dd_fill )THEN
880	id_detect(ji,jj,jk)= 0
881	ENDIF
882
883	ENDIF
884
885	ENDDO
886	ENDDO
887	ENDDO
888
889	IF( ALLOCATED(dl_coef) ) DEALLOCATE( dl_coef )
890
891	END SELECT
892
893	END SUBROUTINE extrap__3D
894	!> @endcode
895	!-------------------------------------------------------------------
896	!> @brief
897	!> This function compute derivative of a function in i- and
898	!> j-direction
899	!>
900	!> @details
901	!>
902	!> @author J.Paul
903	!> - Nov, 2013- Initial Version
904	!
905	!> @param[in] dd_value : 1D table of variable to be extrapolated
906	!> @param[in] dd_fill : FillValue of variable
907	!> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension
908	!-------------------------------------------------------------------
909	!> @code
910	PURE FUNCTION extrap_deriv_1D( dd_value, dd_fill, ld_discont )
911
912	IMPLICIT NONE
913	! Argument
914	REAL(dp) , DIMENSION(:), INTENT(IN) :: dd_value
915	REAL(dp) , INTENT(IN) :: dd_fill
916	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
917
918	! function
919	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1) ) :: extrap_deriv_1D
920
921	! local variable
922	INTEGER(i4) :: il_imin
923	INTEGER(i4) :: il_imax
924	INTEGER(i4), DIMENSION(1) :: il_shape
925
926	REAL(dp) :: dl_min
927	REAL(dp) :: dl_max
928	REAL(dp) , DIMENSION(:), ALLOCATABLE :: dl_value
929
930	LOGICAL :: ll_discont
931
932	! loop indices
933	INTEGER(i4) :: ji
934
935	INTEGER(i4) :: i1
936	INTEGER(i4) :: i2
937	!----------------------------------------------------------------
938	! init
939	extrap_deriv_1D(:)=dd_fill
940
941	ll_discont=.FALSE.
942	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
943
944	il_shape(:)=SHAPE(dd_value(:))
945
946	ALLOCATE( dl_value(3))
947
948	! compute derivative in i-direction
949	DO ji=1,il_shape(1)
950
951	il_imin=MAX(ji-1,1)
952	il_imax=MIN(ji+1,il_shape(1))
953
954	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
955	i1=1 ; i2=3
956	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
957	i1=1 ; i2=2
958	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
959	i1=2 ; i2=3
960	ENDIF
961
962	dl_value(i1:i2)=dd_value(il_imin:il_imax)
963	IF( il_imin == 1 )THEN
964	dl_value(:)=EOSHIFT( dl_value(:), &
965	& DIM=1, &
966	& SHIFT=-1, &
967	& BOUNDARY=dl_value(1) )
968	ENDIF
969	IF( il_imax == il_shape(1) )THEN
970	dl_value(:)=EOSHIFT( dl_value(:), &
971	& DIM=1, &
972	& SHIFT=1, &
973	& BOUNDARY=dl_value(3))
974	ENDIF
975
976	IF( ll_discont )THEN
977	dl_min=MINVAL( dl_value(:), dl_value(:)/=dd_fill )
978	dl_max=MAXVAL( dl_value(:), dl_value(:)/=dd_fill )
979	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
980	WHERE( dl_value(:) < 0._dp )
981	dl_value(:) = dl_value(:)+360._dp
982	END WHERE
983	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
984	WHERE( dl_value(:) > 180._dp )
985	dl_value(:) = dl_value(:)-180._dp
986	END WHERE
987	ENDIF
988	ENDIF
989
990	IF( dl_value( 2) /= dd_fill .AND. & ! ji
991	& dl_value( 3) /= dd_fill .AND. & ! ji+1
992	& dl_value( 1) /= dd_fill )THEN ! ji-1
993
994	extrap_deriv_1D(ji)=&
995	& ( dl_value(3) - dl_value(1) ) / &
996	& REAL( il_imax-il_imin ,dp)
997
998	ENDIF
999
1000	ENDDO
1001
1002	DEALLOCATE( dl_value )
1003
1004	END FUNCTION extrap_deriv_1D
1005	!> @endcode
1006	!-------------------------------------------------------------------
1007	!> @brief
1008	!> This function compute derivative of a function in i- and
1009	!> j-direction
1010	!>
1011	!> @details
1012	!>
1013	!> @author J.Paul
1014	!> - Nov, 2013- Initial Version
1015	!
1016	!> @param[in] dd_value : 2D table of variable to be extrapolated
1017	!> @param[in] dd_fill : FillValue of variable
1018	!> @param[in] cd_dim : compute derivative on first (I) or second (J) dimension
1019	!-------------------------------------------------------------------
1020	!> @code
1021	FUNCTION extrap_deriv_2D( dd_value, dd_fill, cd_dim, ld_discont )
1022
1023	IMPLICIT NONE
1024	! Argument
1025	REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_value
1026	REAL(dp) , INTENT(IN) :: dd_fill
1027	CHARACTER(LEN=*) , INTENT(IN) :: cd_dim
1028	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
1029
1030	! function
1031	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &
1032	& SIZE(dd_value,DIM=2) ) :: extrap_deriv_2D
1033
1034	! local variable
1035	INTEGER(i4) :: il_imin
1036	INTEGER(i4) :: il_imax
1037	INTEGER(i4) :: il_jmin
1038	INTEGER(i4) :: il_jmax
1039	INTEGER(i4), DIMENSION(2) :: il_shape
1040
1041	REAL(dp) :: dl_min
1042	REAL(dp) :: dl_max
1043	REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_value
1044
1045	LOGICAL :: ll_discont
1046
1047	! loop indices
1048	INTEGER(i4) :: ji
1049	INTEGER(i4) :: jj
1050
1051	INTEGER(i4) :: i1
1052	INTEGER(i4) :: i2
1053
1054	INTEGER(i4) :: j1
1055	INTEGER(i4) :: j2
1056	!----------------------------------------------------------------
1057	! init
1058	extrap_deriv_2D(:,:)=dd_fill
1059
1060	ll_discont=.FALSE.
1061	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1062
1063	il_shape(:)=SHAPE(dd_value(:,:))
1064
1065	SELECT CASE(TRIM(fct_upper(cd_dim)))
1066
1067	CASE('I')
1068
1069	ALLOCATE( dl_value(3,il_shape(2)) )
1070	! compute derivative in i-direction
1071	DO ji=1,il_shape(1)
1072
1073	! init
1074	dl_value(:,:)=dd_fill
1075
1076	il_imin=MAX(ji-1,1)
1077	il_imax=MIN(ji+1,il_shape(1))
1078
1079	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
1080	i1=1 ; i2=3
1081	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
1082	i1=1 ; i2=2
1083	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
1084	i1=2 ; i2=3
1085	ENDIF
1086
1087	dl_value(i1:i2,:)=dd_value(il_imin:il_imax,:)
1088	IF( il_imin == 1 )THEN
1089	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1090	& DIM=1, &
1091	& SHIFT=-1, &
1092	& BOUNDARY=dl_value(1,:) )
1093	ENDIF
1094	IF( il_imax == il_shape(1) )THEN
1095	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1096	& DIM=1, &
1097	& SHIFT=1, &
1098	& BOUNDARY=dl_value(3,:))
1099	ENDIF
1100
1101	IF( ll_discont )THEN
1102	dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1103	dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1104	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1105	WHERE( dl_value(:,:) < 0_dp )
1106	dl_value(:,:) = dl_value(:,:)+360._dp
1107	END WHERE
1108	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1109	WHERE( dl_value(:,:) > 180 )
1110	dl_value(:,:) = dl_value(:,:)-180._dp
1111	END WHERE
1112	ENDIF
1113	ENDIF
1114
1115	WHERE( dl_value(2,:) /= dd_fill .AND. & ! ji
1116	& dl_value(3,:) /= dd_fill .AND. & ! ji+1
1117	& dl_value(1,:) /= dd_fill ) ! ji-1
1118
1119	extrap_deriv_2D(ji,:)=&
1120	& ( dl_value(3,:) - dl_value(1,:) ) / &
1121	& REAL( il_imax-il_imin,dp)
1122
1123	END WHERE
1124
1125	ENDDO
1126
1127	CASE('J')
1128
1129	ALLOCATE( dl_value(il_shape(1),3) )
1130	! compute derivative in j-direction
1131	DO jj=1,il_shape(2)
1132
1133	il_jmin=MAX(jj-1,1)
1134	il_jmax=MIN(jj+1,il_shape(2))
1135
1136	IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN
1137	j1=1 ; j2=3
1138	ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN
1139	j1=1 ; j2=2
1140	ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN
1141	j1=2 ; j2=3
1142	ENDIF
1143
1144	dl_value(:,j1:j2)=dd_value(:,il_jmin:il_jmax)
1145	IF( il_jmin == 1 )THEN
1146	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1147	& DIM=2, &
1148	& SHIFT=-1, &
1149	& BOUNDARY=dl_value(:,1))
1150	ENDIF
1151	IF( il_jmax == il_shape(2) )THEN
1152	dl_value(:,:)=EOSHIFT( dl_value(:,:), &
1153	& DIM=2, &
1154	& SHIFT=1, &
1155	& BOUNDARY=dl_value(:,3))
1156	ENDIF
1157
1158	IF( ll_discont )THEN
1159	dl_min=MINVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1160	dl_max=MAXVAL( dl_value(:,:), dl_value(:,:)/=dd_fill )
1161	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1162	WHERE( dl_value(:,:) < 0_dp )
1163	dl_value(:,:) = dl_value(:,:)+360._dp
1164	END WHERE
1165	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1166	WHERE( dl_value(:,:) > 180 )
1167	dl_value(:,:) = dl_value(:,:)-180._dp
1168	END WHERE
1169	ENDIF
1170	ENDIF
1171
1172	WHERE( dl_value(:, 2) /= dd_fill .AND. & ! jj
1173	& dl_value(:, 3) /= dd_fill .AND. & ! jj+1
1174	& dl_value(:, 1) /= dd_fill ) ! jj-1
1175
1176	extrap_deriv_2D(:,jj)=&
1177	& ( dl_value(:,3) - dl_value(:,1) ) / &
1178	& REAL(il_jmax-il_jmin,dp)
1179
1180	END WHERE
1181
1182	ENDDO
1183
1184	END SELECT
1185
1186	DEALLOCATE( dl_value )
1187
1188	END FUNCTION extrap_deriv_2D
1189	!> @endcode
1190	!-------------------------------------------------------------------
1191	!> @brief
1192	!> This function compute derivative of a function in i- and
1193	!> j-direction
1194	!>
1195	!> @details
1196	!>
1197	!> @author J.Paul
1198	!> - Nov, 2013- Initial Version
1199	!
1200	!> @param[inout] dd_value : 3D table of variable to be extrapolated
1201	!> @param[in] dd_fill : FillValue of variable
1202	!> @param[in] cd_dim : compute derivative on first (I) second (J) or third (K) dimension
1203	!-------------------------------------------------------------------
1204	!> @code
1205	PURE FUNCTION extrap_deriv_3D( dd_value, dd_fill, cd_dim, ld_discont )
1206
1207	IMPLICIT NONE
1208	! Argument
1209	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1210	REAL(dp) , INTENT(IN) :: dd_fill
1211	CHARACTER(LEN=*) , INTENT(IN) :: cd_dim
1212	LOGICAL , INTENT(IN), OPTIONAL :: ld_discont
1213
1214	! function
1215	REAL(dp), DIMENSION(SIZE(dd_value,DIM=1), &
1216	& SIZE(dd_value,DIM=2), &
1217	& SIZE(dd_value,DIM=3)) :: extrap_deriv_3D
1218
1219	! local variable
1220	INTEGER(i4) :: il_imin
1221	INTEGER(i4) :: il_imax
1222	INTEGER(i4) :: il_jmin
1223	INTEGER(i4) :: il_jmax
1224	INTEGER(i4) :: il_kmin
1225	INTEGER(i4) :: il_kmax
1226	INTEGER(i4), DIMENSION(3) :: il_shape
1227
1228	REAL(dp) :: dl_min
1229	REAL(dp) :: dl_max
1230	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value
1231
1232	LOGICAL :: ll_discont
1233
1234	! loop indices
1235	INTEGER(i4) :: ji
1236	INTEGER(i4) :: jj
1237	INTEGER(i4) :: jk
1238
1239	INTEGER(i4) :: i1
1240	INTEGER(i4) :: i2
1241
1242	INTEGER(i4) :: j1
1243	INTEGER(i4) :: j2
1244
1245	INTEGER(i4) :: k1
1246	INTEGER(i4) :: k2
1247	!----------------------------------------------------------------
1248	! init
1249	extrap_deriv_3D(:,:,:)=dd_fill
1250
1251	ll_discont=.FALSE.
1252	IF( PRESENT(ld_discont) ) ll_discont=ld_discont
1253
1254	il_shape(:)=SHAPE(dd_value(:,:,:))
1255
1256
1257	SELECT CASE(TRIM(fct_upper(cd_dim)))
1258
1259	CASE('I')
1260
1261	ALLOCATE( dl_value(3,il_shape(2),il_shape(3)) )
1262	! compute derivative in i-direction
1263	DO ji=1,il_shape(1)
1264
1265	il_imin=MAX(ji-1,1)
1266	il_imax=MIN(ji+1,il_shape(1))
1267
1268	IF( il_imin==ji-1 .AND. il_imax==ji+1 )THEN
1269	i1=1 ; i2=3
1270	ELSEIF( il_imin==ji .AND. il_imax==ji+1 )THEN
1271	i1=1 ; i2=2
1272	ELSEIF( il_imin==ji-1 .AND. il_imax==ji )THEN
1273	i1=2 ; i2=3
1274	ENDIF
1275
1276	dl_value(i1:i2,:,:)=dd_value(il_imin:il_imax,:,:)
1277	IF( il_imin == 1 )THEN
1278	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1279	& DIM=1, &
1280	& SHIFT=-1, &
1281	& BOUNDARY=dl_value(1,:,:) )
1282	ENDIF
1283	IF( il_imax == il_shape(1) )THEN
1284	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1285	& DIM=1, &
1286	& SHIFT=1, &
1287	& BOUNDARY=dl_value(3,:,:))
1288	ENDIF
1289
1290	IF( ll_discont )THEN
1291	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1292	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1293	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1294	WHERE( dl_value(:,:,:) < 0_dp )
1295	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1296	END WHERE
1297	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1298	WHERE( dl_value(:,:,:) > 180 )
1299	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1300	END WHERE
1301	ENDIF
1302	ENDIF
1303
1304	WHERE( dl_value(2,:,:) /= dd_fill .AND. & ! ji
1305	& dl_value(3,:,:) /= dd_fill .AND. & !ji+1
1306	& dl_value(1,:,:) /= dd_fill ) !ji-1
1307
1308	extrap_deriv_3D(ji,:,:)= &
1309	& ( dl_value(3,:,:) - dl_value(1,:,:) ) / &
1310	& REAL( il_imax-il_imin ,dp)
1311
1312	END WHERE
1313
1314	ENDDO
1315
1316	CASE('J')
1317
1318	ALLOCATE( dl_value(il_shape(1),3,il_shape(3)) )
1319	! compute derivative in j-direction
1320	DO jj=1,il_shape(2)
1321
1322	il_jmin=MAX(jj-1,1)
1323	il_jmax=MIN(jj+1,il_shape(2))
1324
1325	IF( il_jmin==jj-1 .AND. il_jmax==jj+1 )THEN
1326	j1=1 ; j2=3
1327	ELSEIF( il_jmin==jj .AND. il_jmax==jj+1 )THEN
1328	j1=1 ; j2=2
1329	ELSEIF( il_jmin==jj-1 .AND. il_jmax==jj )THEN
1330	j1=2 ; j2=3
1331	ENDIF
1332
1333	dl_value(:,j1:j2,:)=dd_value(:,il_jmin:il_jmax,:)
1334	IF( il_jmin == 1 )THEN
1335	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1336	& DIM=2, &
1337	& SHIFT=-1, &
1338	& BOUNDARY=dl_value(:,1,:) )
1339	ENDIF
1340	IF( il_jmax == il_shape(2) )THEN
1341	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1342	& DIM=2, &
1343	& SHIFT=1, &
1344	& BOUNDARY=dl_value(:,3,:))
1345	ENDIF
1346
1347	IF( ll_discont )THEN
1348	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1349	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1350	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1351	WHERE( dl_value(:,:,:) < 0_dp )
1352	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1353	END WHERE
1354	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1355	WHERE( dl_value(:,:,:) > 180 )
1356	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1357	END WHERE
1358	ENDIF
1359	ENDIF
1360
1361	WHERE( dl_value(:, 2,:) /= dd_fill .AND. & ! jj
1362	& dl_value(:, 3,:) /= dd_fill .AND. & ! jj+1
1363	& dl_value(:, 1,:) /= dd_fill ) ! jj-1
1364
1365	extrap_deriv_3D(:,jj,:)=&
1366	& ( dl_value(:,3,:) - dl_value(:,1,:) ) / &
1367	& REAL( il_jmax - il_jmin ,dp)
1368
1369	END WHERE
1370
1371	ENDDO
1372
1373	CASE('K')
1374	! compute derivative in k-direction
1375	DO jk=1,il_shape(3)
1376
1377	il_kmin=MAX(jk-1,1)
1378	il_kmax=MIN(jk+1,il_shape(3))
1379
1380	IF( il_kmin==jk-1 .AND. il_kmax==jk+1 )THEN
1381	k1=1 ; k2=3
1382	ELSEIF( il_kmin==jk .AND. il_kmax==jk+1 )THEN
1383	k1=1 ; k2=2
1384	ELSEIF( il_kmin==jk-1 .AND. il_kmax==jk )THEN
1385	k1=2 ; k2=3
1386	ENDIF
1387
1388	dl_value(:,:,k1:k2)=dd_value(:,:,il_kmin:il_kmax)
1389	IF( il_kmin == 1 )THEN
1390	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1391	& DIM=3, &
1392	& SHIFT=-1, &
1393	& BOUNDARY=dl_value(:,:,1) )
1394	ENDIF
1395	IF( il_kmax == il_shape(3) )THEN
1396	dl_value(:,:,:)=EOSHIFT( dl_value(:,:,:), &
1397	& DIM=3, &
1398	& SHIFT=1, &
1399	& BOUNDARY=dl_value(:,:,3))
1400	ENDIF
1401
1402	IF( ll_discont )THEN
1403	dl_min=MINVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1404	dl_max=MAXVAL( dl_value(:,:,:), dl_value(:,:,:)/=dd_fill )
1405	IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
1406	WHERE( dl_value(:,:,:) < 0_dp )
1407	dl_value(:,:,:) = dl_value(:,:,:)+360._dp
1408	END WHERE
1409	ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
1410	WHERE( dl_value(:,:,:) > 180 )
1411	dl_value(:,:,:) = dl_value(:,:,:)-180._dp
1412	END WHERE
1413	ENDIF
1414	ENDIF
1415
1416	WHERE( dl_value(:,:, 2) /= dd_fill .AND. & ! jk
1417	& dl_value(:,:, 3) /= dd_fill .AND. & ! jk+1
1418	& dl_value(:,:, 1) /= dd_fill ) ! jk-1
1419
1420	extrap_deriv_3D(:,:,jk)=&
1421	& ( dl_value(:,:,3) - dl_value(:,:,1) ) / &
1422	& REAL( il_kmax-il_kmin,dp)
1423
1424	END WHERE
1425
1426	ENDDO
1427
1428	END SELECT
1429
1430	DEALLOCATE( dl_value )
1431
1432	END FUNCTION extrap_deriv_3D
1433	!> @endcode
1434	!-------------------------------------------------------------------
1435	!> @brief
1436	!> This function compute extrapolatd values by calculated minimum error using
1437	!> taylor expansion
1438	!>
1439	!> @details
1440	!>
1441	!> @author J.Paul
1442	!> - Nov, 2013- Initial Version
1443	!
1444	!> @param[in] dd_value : 3D table of variable to be extrapolated
1445	!> @param[in] dd_fill : FillValue of variable
1446	!> @param[in] dd_ideriv : derivative of function in i-direction
1447	!> @param[in] dd_jderiv : derivative of function in j-direction
1448	!> @param[in] dd_kderiv : derivative of function in k-direction
1449	!> @param[in] id_ji : i-direction indices table
1450	!> @param[in] id_jj : j-direction indices table
1451	!> @param[in] id_ii : i-direction indices of the point to extrapolate
1452	!> @param[in] id_ij : j-direction indices of the point to extrapolate
1453	!-------------------------------------------------------------------
1454	!> @code
1455	PURE FUNCTION extrap__3D_min_error_coef( dd_value )
1456
1457	IMPLICIT NONE
1458	! Argument
1459	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1460
1461	! function
1462	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
1463	& SIZE(dd_value(:,:,:),DIM=2), &
1464	& SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_min_error_coef
1465
1466	! local variable
1467	INTEGER(i4), DIMENSION(3) :: il_shape
1468
1469	INTEGER(i4) :: il_imid
1470	INTEGER(i4) :: il_jmid
1471	INTEGER(i4) :: il_kmid
1472
1473	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
1474
1475	! loop indices
1476	INTEGER(i4) :: ji
1477	INTEGER(i4) :: jj
1478	INTEGER(i4) :: jk
1479	!----------------------------------------------------------------
1480
1481	! init
1482	extrap__3D_min_error_coef(:,:,:)=0
1483
1484	il_shape(:)=SHAPE(dd_value(:,:,:))
1485
1486	il_imid=INT(REAL(il_shape(1),sp)*0.5+1)
1487	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1)
1488	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1)
1489
1490	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
1491
1492	DO jk=1,il_shape(3)
1493	DO jj=1,il_shape(2)
1494	DO ji=1,il_shape(1)
1495
1496	! compute distance
1497	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
1498	& (jj-il_jmid)**2 + &
1499	& (jk-il_kmid)**2
1500
1501	IF( dl_dist(ji,jj,jk) /= 0 )THEN
1502	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
1503	ENDIF
1504
1505	ENDDO
1506	ENDDO
1507	ENDDO
1508
1509	WHERE( dl_dist(:,:,:) /= 0 )
1510	extrap__3D_min_error_coef(:,:,:)=dl_dist(:,:,:)
1511	END WHERE
1512
1513	DEALLOCATE( dl_dist )
1514
1515	END FUNCTION extrap__3D_min_error_coef
1516	!> @endcode
1517	!-------------------------------------------------------------------
1518	!> @brief
1519	!> This function compute extrapolatd values by calculated minimum error using
1520	!> taylor expansion
1521	!>
1522	!> @details
1523	!>
1524	!> @author J.Paul
1525	!> - Nov, 2013- Initial Version
1526	!
1527	!> @param[in] dd_value : 3D table of variable to be extrapolated
1528	!> @param[in] dd_fill : FillValue of variable
1529	!> @param[in] dd_dfdx : derivative of function in i-direction
1530	!> @param[in] dd_dfdy : derivative of function in j-direction
1531	!> @param[in] dd_dfdz : derivative of function in k-direction
1532	!> @param[in] dd_coef :
1533	!-------------------------------------------------------------------
1534	!> @code
1535	PURE FUNCTION extrap__3D_min_error_fill( dd_value, dd_fill, id_ext, &
1536	& dd_dfdx, dd_dfdy, dd_dfdz, &
1537	& dd_coef )
1538	IMPLICIT NONE
1539	! Argument
1540	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1541	REAL(dp) , INTENT(IN) :: dd_fill
1542	INTEGER(i4), INTENT(IN) :: id_ext
1543	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx
1544	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy
1545	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdz
1546	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
1547
1548	! function
1549	REAL(dp) :: extrap__3d_min_error_fill
1550
1551	! local variable
1552	INTEGER(i4), DIMENSION(3) :: il_shape
1553	INTEGER(i4), DIMENSION(3) :: il_ind
1554
1555	INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_mask
1556	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_error
1557
1558	INTEGER(i4) :: il_min
1559	! loop indices
1560
1561	!----------------------------------------------------------------
1562
1563	! init
1564	extrap__3D_min_error_fill=dd_fill
1565
1566	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2))
1567
1568	IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN
1569
1570	il_shape(:)=SHAPE(dd_value(:,:,:))
1571	ALLOCATE( il_mask( il_shape(1),il_shape(2),il_shape(3)) )
1572	ALLOCATE( dl_error(il_shape(1),il_shape(2),il_shape(3)) )
1573
1574	! compute error
1575	dl_error(:,:,:)=0.
1576	il_mask(:,:,:)=0
1577	WHERE( dd_dfdx(:,:,:) /= dd_fill )
1578	dl_error(:,:,:)=dd_coef(:,:,:)*dd_dfdx(:,:,:)
1579	il_mask(:,:,:)=1
1580	END WHERE
1581	WHERE( dd_dfdy(:,:,:) /= dd_fill )
1582	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdy(:,:,:))
1583	il_mask(:,:,:)=1
1584	END WHERE
1585	WHERE( dd_dfdz(:,:,:) /= dd_fill )
1586	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdz(:,:,:))
1587	il_mask(:,:,:)=1
1588	END WHERE
1589
1590	! get minimum error indices
1591	il_ind(:)=MINLOC(dl_error(:,:,:),il_mask(:,:,:)==1)
1592
1593	! return value
1594	IF( ALL(il_ind(:)/=0) )THEN
1595	extrap__3D_min_error_fill=dd_value(il_ind(1),il_ind(2),il_ind(3))
1596	ENDIF
1597
1598	DEALLOCATE( il_mask )
1599	DEALLOCATE( dl_error )
1600
1601	ENDIF
1602
1603	END FUNCTION extrap__3D_min_error_fill
1604	!> @endcode
1605	!-------------------------------------------------------------------
1606	!> @brief
1607	!> This function compute extrapolatd values by calculated minimum error using
1608	!> taylor expansion
1609	!>
1610	!> @details
1611	!>
1612	!> @author J.Paul
1613	!> - Nov, 2013- Initial Version
1614	!
1615	!> @param[in] dd_value : 3D table of variable to be extrapolated
1616	!-------------------------------------------------------------------
1617	!> @code
1618	PURE FUNCTION extrap__3D_dist_weight_coef( dd_value )
1619
1620	IMPLICIT NONE
1621	! Argument
1622	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1623
1624	! function
1625	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
1626	& SIZE(dd_value(:,:,:),DIM=2), &
1627	& SIZE(dd_value(:,:,:),DIM=3) ) :: extrap__3D_dist_weight_coef
1628
1629	! local variable
1630	INTEGER(i4), DIMENSION(3) :: il_shape
1631
1632	INTEGER(i4) :: il_imid
1633	INTEGER(i4) :: il_jmid
1634	INTEGER(i4) :: il_kmid
1635
1636	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
1637
1638	! loop indices
1639	INTEGER(i4) :: ji
1640	INTEGER(i4) :: jj
1641	INTEGER(i4) :: jk
1642	!----------------------------------------------------------------
1643
1644	! init
1645	extrap__3D_dist_weight_coef(:,:,:)=0
1646
1647	il_shape(:)=SHAPE(dd_value(:,:,:))
1648
1649	il_imid=INT(REAL(il_shape(1),sp)*0.5+1,i4)
1650	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1,i4)
1651	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1,i4)
1652
1653	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
1654
1655	DO jk=1,il_shape(3)
1656	DO jj=1,il_shape(2)
1657	DO ji=1,il_shape(1)
1658
1659	! compute distance
1660	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
1661	& (jj-il_jmid)**2 + &
1662	& (jk-il_kmid)**2
1663
1664	IF( dl_dist(ji,jj,jk) /= 0 )THEN
1665	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
1666	ENDIF
1667
1668	ENDDO
1669	ENDDO
1670	ENDDO
1671
1672	WHERE( dl_dist(:,:,:) /= 0 )
1673	extrap__3D_dist_weight_coef(:,:,:)=1./dl_dist(:,:,:)
1674	END WHERE
1675
1676	DEALLOCATE( dl_dist )
1677
1678	END FUNCTION extrap__3D_dist_weight_coef
1679	!> @endcode
1680	!-------------------------------------------------------------------
1681	!> @brief
1682	!> This function compute extrapolatd values by calculated minimum error using
1683	!> taylor expansion
1684	!>
1685	!> @details
1686	!>
1687	!> @author J.Paul
1688	!> - Nov, 2013- Initial Version
1689	!
1690	!> @param[in] dd_value : 3D table of variable to be extrapolated
1691	!> @param[in] dd_fill : FillValue of variable
1692	!> @param[in] dd_coef :
1693	!-------------------------------------------------------------------
1694	!> @code
1695	FUNCTION extrap__3D_dist_weight_fill( dd_value, dd_fill, id_ext, &
1696	& dd_coef )
1697	IMPLICIT NONE
1698	! Argument
1699	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1700	REAL(dp) , INTENT(IN) :: dd_fill
1701	INTEGER(i4), INTENT(IN) :: id_ext
1702	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
1703
1704	! function
1705	REAL(dp) :: extrap__3D_dist_weight_fill
1706
1707	! local variable
1708	INTEGER(i4), DIMENSION(3) :: il_shape
1709
1710	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value
1711	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
1712
1713	INTEGER(i4) :: il_min
1714	! loop indices
1715	INTEGER(i4) :: ji
1716	INTEGER(i4) :: jj
1717	INTEGER(i4) :: jk
1718
1719	!----------------------------------------------------------------
1720
1721	! init
1722	extrap__3D_dist_weight_fill=dd_fill
1723
1724	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_ext*2))
1725
1726	IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN
1727
1728	il_shape(:)=SHAPE(dd_value(:,:,:))
1729	ALLOCATE( dl_value( il_shape(1),il_shape(2),il_shape(3)) )
1730	ALLOCATE( dl_coef( il_shape(1),il_shape(2),il_shape(3)) )
1731
1732	dl_value(:,:,:)=0
1733	dl_coef(:,:,:)=0
1734
1735	FORALL( ji=1:il_shape(1), &
1736	& jj=1:il_shape(2), &
1737	& jk=1:il_shape(3), &
1738	& dd_value(ji,jj,jk) /= dd_fill )
1739
1740	! compute factor
1741	dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk)
1742	dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk)
1743
1744	END FORALL
1745
1746	! return value
1747	IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN
1748	extrap__3D_dist_weight_fill=SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) )
1749	ENDIF
1750
1751	DEALLOCATE( dl_value )
1752	DEALLOCATE( dl_coef )
1753
1754	ENDIF
1755
1756	END FUNCTION extrap__3D_dist_weight_fill
1757	!> @endcode
1758	!-------------------------------------------------------------------
1759	!> @brief
1760	!> This subroutine add to the variable (to be extrapolated) an
1761	!> extraband of N points at north,south,east and west boundaries.
1762	!>
1763	!> @author J.Paul
1764	!> - 2013-Initial version
1765	!
1766	!> @param[inout] td_var : variable
1767	!> @param[in] id_isize : i-direction size of extra bands (default=im_minext)
1768	!> @param[in] id_jsize : j-direction size of extra bands (default=im_minext)
1769	!> @todo
1770	!> - invalid special case for grid with north fold
1771	!-------------------------------------------------------------------
1772	!> @code
1773	SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize )
1774	IMPLICIT NONE
1775	! Argument
1776	TYPE(TVAR) , INTENT(INOUT) :: td_var
1777	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
1778	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
1779
1780	! local variable
1781	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
1782
1783	INTEGER(i4) :: il_isize
1784	INTEGER(i4) :: il_jsize
1785	INTEGER(i4) :: il_tmp
1786
1787	! loop indices
1788	INTEGER(i4) :: ji
1789	INTEGER(i4) :: ij
1790	!----------------------------------------------------------------
1791	il_isize=im_minext
1792	IF(PRESENT(id_isize)) il_isize=id_isize
1793	IF( il_isize < im_minext .AND. &
1794	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
1795	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
1796	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
1797	& " cubic interpolation ")
1798	ENDIF
1799
1800	il_jsize=im_minext
1801	IF(PRESENT(id_jsize)) il_jsize=id_jsize
1802	IF( il_jsize < im_minext .AND. &
1803	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
1804	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
1805	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
1806	& " cubic interpolation ")
1807	ENDIF
1808
1809	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
1810	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
1811
1812	CALL logger_trace( "EXTRAP ADD EXTRABANDS: dimension change "//&
1813	& "in variable "//TRIM(td_var%c_name) )
1814
1815	! add extrabands in variable
1816	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
1817	& td_var%t_dim(2)%i_len, &
1818	& td_var%t_dim(3)%i_len, &
1819	& td_var%t_dim(4)%i_len ))
1820
1821	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
1822
1823	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize
1824	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize
1825
1826	DEALLOCATE(td_var%d_value)
1827	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
1828	& td_var%t_dim(2)%i_len, &
1829	& td_var%t_dim(3)%i_len, &
1830	& td_var%t_dim(4)%i_len ) )
1831
1832	! intialise
1833	td_var%d_value(:,:,:,:)=td_var%d_fill
1834
1835	! fill center
1836	td_var%d_value( 1+il_isize:td_var%t_dim(1)%i_len-il_isize, &
1837	& 1+il_jsize:td_var%t_dim(2)%i_len-il_jsize, &
1838	& :,:) = dl_value(:,:,:,:)
1839
1840	! special case for overlap
1841	IF( td_var%i_ew >= 0 .AND. il_isize /= 0 )THEN
1842	DO ji=1,il_isize
1843	! from east to west
1844	il_tmp=td_var%t_dim(1)%i_len-td_var%i_ew+ji-2*il_isize
1845	td_var%d_value(ji,:,:,:) = td_var%d_value(il_tmp,:,:,:)
1846
1847	! from west to east
1848	ij=td_var%t_dim(1)%i_len-ji+1
1849	il_tmp=td_var%i_ew-ji+2*il_isize+1
1850	td_var%d_value(ij,:,:,:) = td_var%d_value(il_tmp,:,:,:)
1851	ENDDO
1852	ENDIF
1853
1854	DEALLOCATE( dl_value )
1855
1856	END SUBROUTINE extrap_add_extrabands
1857	!> @endcode
1858	!-------------------------------------------------------------------
1859	!> @brief
1860	!> This subroutine remove of the variable an extraband
1861	!> of N points at north,south,east and west boundaries.
1862	!>
1863	!> @author J.Paul
1864	!> - 2013-Initial version
1865	!
1866	!> @param[inout] td_var : variable
1867	!> @param[in] id_isize : i-direction size of extra bands (default=im_minext)
1868	!> @param[in] id_jsize : j-direction size of extra bands (default=im_minext)
1869	!-------------------------------------------------------------------
1870	!> @code
1871	SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize )
1872	IMPLICIT NONE
1873	! Argument
1874	TYPE(TVAR) , INTENT(INOUT) :: td_var
1875	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
1876	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
1877
1878	! local variable
1879	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
1880
1881	INTEGER(i4) :: il_isize
1882	INTEGER(i4) :: il_jsize
1883
1884	INTEGER(i4) :: il_imin
1885	INTEGER(i4) :: il_imax
1886	INTEGER(i4) :: il_jmin
1887	INTEGER(i4) :: il_jmax
1888
1889	! loop indices
1890	!----------------------------------------------------------------
1891	il_isize=im_minext
1892	IF(PRESENT(id_isize)) il_isize=id_isize
1893
1894	il_jsize=im_minext
1895	IF(PRESENT(id_jsize)) il_jsize=id_jsize
1896
1897	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
1898	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
1899
1900	CALL logger_trace( "EXTRAP DEL EXTRABANDS: dimension change "//&
1901	& "in variable "//TRIM(td_var%c_name) )
1902
1903	! add extrabands in variable
1904	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
1905	& td_var%t_dim(2)%i_len, &
1906	& td_var%t_dim(3)%i_len, &
1907	& td_var%t_dim(4)%i_len ))
1908
1909	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
1910
1911	! fill center
1912	il_imin=1+il_isize
1913	il_imax=td_var%t_dim(1)%i_len-il_isize
1914
1915	il_jmin=1+il_jsize
1916	il_jmax=td_var%t_dim(2)%i_len-il_jsize
1917
1918	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len - 2*il_isize
1919	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len - 2*il_jsize
1920
1921	DEALLOCATE(td_var%d_value)
1922	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
1923	& td_var%t_dim(2)%i_len, &
1924	& td_var%t_dim(3)%i_len, &
1925	& td_var%t_dim(4)%i_len ) )
1926
1927	! intialise
1928	td_var%d_value(:,:,:,:)=td_var%d_fill
1929
1930	td_var%d_value(:,:,:,:)=dl_value(il_imin:il_imax,&
1931	& il_jmin:il_jmax,&
1932	& :,:)
1933
1934	DEALLOCATE( dl_value )
1935
1936	END SUBROUTINE extrap_del_extrabands
1937	!> @endcode
1938	! !-------------------------------------------------------------------
1939	! !> @brief
1940	! !> This function
1941	! !>
1942	! !> @details
1943	! !>
1944	! !> @author J.Paul
1945	! !> - Nov, 2013- Initial Version
1946	! !
1947	! !> @param[in]
1948	! !> @param[out]
1949	! !-------------------------------------------------------------------
1950	! !> @code
1951	! FUNCTION extrap_( )
1952	! IMPLICIT NONE
1953	! ! Argument
1954	!
1955	! ! local variable
1956	!
1957	! ! loop indices
1958	! !----------------------------------------------------------------
1959	! END FUNCTION extrap_
1960	! !> @endcode
1961	! !-------------------------------------------------------------------
1962	! !> @brief
1963	! !> This subroutine
1964	! !>
1965	! !> @details
1966	! !>
1967	! !> @author J.Paul
1968	! !> - Nov, 2013- Initial Version
1969	! !
1970	! !> @param[in]
1971	! !> @param[out]
1972	! !-------------------------------------------------------------------
1973	! !> @code
1974	! SUBROUTINE extrap_( )
1975	! IMPLICIT NONE
1976	! ! Argument
1977	!
1978	! ! local variable
1979	!
1980	! ! loop indices
1981	! !----------------------------------------------------------------
1982	! END SUBROUTINE extrap_
1983	! !> @endcode
1984	END MODULE extrap

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

Download in other formats: