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 utils/tools/SIREN/src – NEMO

Context Navigation

source: utils/tools/SIREN/src/extrap.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago
update nemo trunk
File size: 51.5 KB

Line
1	!----------------------------------------------------------------------
2	! NEMO system team, System and Interface for oceanic RElocable Nesting
3	!----------------------------------------------------------------------
4	!
5	! DESCRIPTION:
6	!> @brief
7	!> This module manage extrapolation.
8	!>
9	!> @details
10	!> Extrapolation method to be used is specify inside variable
11	!> strcuture, as array of string character.<br/>
12	!> - td_var\%c_extrap(1) string character is the interpolation name choose between:
13	!> - 'dist_weight'
14	!> - 'min_error'
15	!>
16	!> @note Extrapolation method could be specify for each variable in namelist _namvar_,
17	!> defining string character _cn\_varinfo_. By default _dist_weight_.<br/>
18	!> Example:
19	!> - cn_varinfo='varname1:ext=dist_weight', 'varname2:ext=min_error'
20	!>
21	!> to detect point to be extrapolated:<br/>
22	!> @code
23	!> il_detect(:,:,:)=extrap_detect(td_var)
24	!> @endcode
25	!> - il_detect(:,:,:) is 3D array of point to be extrapolated
26	!> - td_var is coarse grid variable to be extrapolated
27	!>
28	!> to extrapolate variable:<br/>
29	!> @code
30	!> CALL extrap_fill_value( td_var, [id_radius])
31	!> @endcode
32	!> - td_var is coarse grid variable to be extrapolated
33	!> - id_radius is radius of the halo used to compute extrapolation [optional]
34	!>
35	!> to add extraband to the variable (to be extrapolated):<br/>
36	!> @code
37	!> CALL extrap_add_extrabands(td_var, [id_isize,] [id_jsize] )
38	!> @endcode
39	!> - td_var is variable structure
40	!> - id_isize : i-direction size of extra bands [optional]
41	!> - id_jsize : j-direction size of extra bands [optional]
42	!>
43	!> to delete extraband of a variable:<br/>
44	!> @code
45	!> CALL extrap_del_extrabands(td_var, [id_isize,] [id_jsize] )
46	!> @endcode
47	!> - td_var is variable structure
48	!> - id_isize : i-direction size of extra bands [optional]
49	!> - id_jsize : j-direction size of extra bands [optional]
50	!>
51	!> @warning _FillValue must not be zero (use var_chg_FillValue())
52	!>
53	!> @author
54	!> J.Paul
55	!>
56	!> @date November, 2013 - Initial Version
57	!> @date September, 2014
58	!> - add header
59	!> @date June, 2015
60	!> - extrapolate all land points (_FillValue)
61	!> - move deriv function to math module
62	!> @date July, 2015
63	!> - compute extrapolation from north west to south east,
64	!> and from south east to north west
65	!>
66	!> @todo
67	!> - create module for each extrapolation method
68	!> - smooth extrapolated points
69	!>
70	!> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
71	!----------------------------------------------------------------------
72	MODULE extrap
73
74	USE netcdf ! nf90 library
75	USE kind ! F90 kind parameter
76	USE phycst ! physical constant
77	USE global ! global variable
78	USE fct ! basic useful function
79	USE date ! date manager
80	USE logger ! log file manager
81	USE math ! mathematical function
82	USE att ! attribute manager
83	USE dim ! dimension manager
84	USE var ! variable manager
85
86	IMPLICIT NONE
87	! NOTE_avoid_public_variables_if_possible
88
89	! type and variable
90	PRIVATE :: im_minext !< default minumum number of point to extrapolate
91	PRIVATE :: im_mincubic !< default minumum number of point to extrapolate for cubic interpolation
92
93	! function and subroutine
94	PUBLIC :: extrap_detect !< detected point to be extrapolated
95	PUBLIC :: extrap_fill_value !< extrapolate value over detected point
96	PUBLIC :: extrap_add_extrabands !< add extraband to the variable (to be extrapolated)
97	PUBLIC :: extrap_del_extrabands !< delete extraband of the variable
98
99	PRIVATE :: extrap__detect_wrapper ! detected point to be extrapolated wrapper
100	PRIVATE :: extrap__detect ! detected point to be extrapolated
101	PRIVATE :: extrap__fill_value_wrapper ! extrapolate value over detected point wrapper
102	PRIVATE :: extrap__fill_value ! extrapolate value over detected point
103	PRIVATE :: extrap__3D !
104	PRIVATE :: extrap__3D_min_error_coef !
105	PRIVATE :: extrap__3D_min_error_fill !
106	PRIVATE :: extrap__3D_dist_weight_coef !
107	PRIVATE :: extrap__3D_dist_weight_fill !
108
109	INTEGER(i4), PARAMETER :: im_minext = 2 !< default minumum number of point to extrapolate
110	INTEGER(i4), PARAMETER :: im_mincubic= 4 !< default minumum number of point to extrapolate for cubic interpolation
111
112	INTERFACE extrap_detect
113	MODULE PROCEDURE extrap__detect_wrapper !< detected point to be extrapolated
114	END INTERFACE extrap_detect
115
116	INTERFACE extrap_fill_value
117	MODULE PROCEDURE extrap__fill_value_wrapper !< detected point to be interpolated
118	END INTERFACE extrap_fill_value
119
120	CONTAINS
121	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122	FUNCTION extrap__detect(td_var0) &
123	& RESULT (if_detect)
124	!-------------------------------------------------------------------
125	!> @brief
126	!> This function detected point to be extrapolated, given variable structure.
127	!>
128	!> @details
129	!> optionaly, you could sepcify fine grid level, refinment factor (default 1),
130	!> offset between fine and coarse grid (default compute from refinment factor
131	!> as offset=(rho-1)/2), number of point to be extrapolated in each direction
132	!> (default im_minext).<br/>
133	!>
134	!> First coarsening fine grid level, if need be, then select point near
135	!> grid point already inform.
136	!>
137	!> @note point to be extrapolated are selected using FillValue,
138	!> so to avoid mistake FillValue should not be zero (use var_chg_FillValue)
139	!>
140	!> @author J.Paul
141	!> @date November, 2013 - Initial Version
142	!> @date June, 2015
143	!> - do not use level to select points to be extrapolated
144	!>
145	!> @param[in] td_var0 coarse grid variable to extrapolate
146	!> @return array of point to be extrapolated
147	!-------------------------------------------------------------------
148
149	IMPLICIT NONE
150
151	! Argument
152	TYPE(TVAR) , INTENT(IN ) :: td_var0
153
154	! function
155	INTEGER(i4), DIMENSION(td_var0%t_dim(1)%i_len,&
156	& td_var0%t_dim(2)%i_len,&
157	& td_var0%t_dim(3)%i_len ) :: if_detect
158
159	! local variable
160	! loop indices
161	INTEGER(i4) :: ji0
162	INTEGER(i4) :: jj0
163	INTEGER(i4) :: jk0
164	!----------------------------------------------------------------
165
166	! force to extrapolated all points
167	if_detect(:,:,:)=1
168
169	! do not compute grid point already inform
170	DO jk0=1,td_var0%t_dim(3)%i_len
171	DO jj0=1,td_var0%t_dim(2)%i_len
172	DO ji0=1,td_var0%t_dim(1)%i_len
173	IF( td_var0%d_value(ji0,jj0,jk0,1) /= td_var0%d_fill )THEN
174	if_detect(ji0,jj0,jk0)=0
175	ENDIF
176	ENDDO
177	ENDDO
178	ENDDO
179
180	END FUNCTION extrap__detect
181	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
182	FUNCTION extrap__detect_wrapper(td_var) &
183	& RESULT (if_detect)
184	!-------------------------------------------------------------------
185	!> @brief
186	!> This function sort variable to be extrapolated, depending on number of
187	!> dimentsion, then detected point to be extrapolated.
188	!>
189	!> @author J.Paul
190	!> @date November, 2013 - Initial Version
191	!> @date June, 2015
192	!> - select all land points for extrapolation
193	!>
194	!> @param[in] td_var coarse grid variable to extrapolate
195	!> @return 3D array of point to be extrapolated
196	!-------------------------------------------------------------------
197
198	IMPLICIT NONE
199
200	! Argument
201	TYPE(TVAR) , INTENT(IN ) :: td_var
202
203	! function
204	INTEGER(i4), DIMENSION(td_var%t_dim(1)%i_len,&
205	& td_var%t_dim(2)%i_len,&
206	& td_var%t_dim(3)%i_len ) :: if_detect
207
208	! local variable
209	! loop indices
210	!----------------------------------------------------------------
211	! init
212	if_detect(:,:,:)=0
213
214	IF( .NOT. ANY(td_var%t_dim(1:3)%l_use) )THEN
215	! no dimension I-J-K used
216	CALL logger_debug(" EXTRAP DETECT: nothing done for variable"//&
217	& TRIM(td_var%c_name) )
218	ELSE IF( ALL(td_var%t_dim(1:3)%l_use) )THEN
219
220	! detect point to be extrapolated on I-J-K
221	CALL logger_debug(" EXTRAP DETECT: detect point "//&
222	& " for variable "//TRIM(td_var%c_name) )
223
224	if_detect(:,:,:)=extrap__detect( td_var )
225
226	ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN
227
228	! detect point to be extrapolated on I-J
229	CALL logger_debug(" EXTRAP DETECT: detect horizontal point "//&
230	& " for variable "//TRIM(td_var%c_name) )
231
232	if_detect(:,:,1:1)=extrap__detect( td_var )
233
234	ELSE IF( td_var%t_dim(3)%l_use )THEN
235
236	! detect point to be extrapolated on K
237	CALL logger_debug(" EXTRAP DETECT: detect vertical point "//&
238	& " for variable "//TRIM(td_var%c_name) )
239
240	if_detect(1:1,1:1,:)=extrap__detect( td_var )
241
242	ENDIF
243
244	CALL logger_debug(" EXTRAP DETECT: "//&
245	& TRIM(fct_str(SUM(if_detect(:,:,:))))//&
246	& " points to be extrapolated" )
247
248	END FUNCTION extrap__detect_wrapper
249	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250	SUBROUTINE extrap__fill_value_wrapper(td_var, id_radius)
251	!-------------------------------------------------------------------
252	!> @brief
253	!> This subroutine select method to be used for extrapolation.
254	!> If need be, increase number of points to be extrapolated.
255	!> Finally launch extrap__fill_value.
256	!>
257	!> @details
258	!> optionaly, you could specify :<br/>
259	!> - refinment factor (default 1)
260	!> - offset between fine and coarse grid (default compute from refinment factor
261	!> as offset=(rho-1)/2)
262	!> - number of point to be extrapolated in each direction (default im_minext)
263	!> - radius of the halo used to compute extrapolation
264	!> - maximum number of iteration
265	!>
266	!> @author J.Paul
267	!> @date November, 2013 - Initial Version
268	!> @date June, 2015
269	!> - select all land points for extrapolation
270	!>
271	!> @param[inout] td_var variable structure
272	!> @param[in] id_radius radius of the halo used to compute extrapolation
273	!-------------------------------------------------------------------
274
275	IMPLICIT NONE
276
277	! Argument
278	TYPE(TVAR) , INTENT(INOUT) :: td_var
279	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_radius
280
281	! local variable
282	INTEGER(i4) :: il_radius
283
284	CHARACTER(LEN=lc) :: cl_method
285
286	! loop indices
287	!----------------------------------------------------------------
288	IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
289	CALL logger_error("EXTRAP FILL VALUE: no value "//&
290	& "associted to variable "//TRIM(td_var%c_name) )
291	ELSE
292
293	SELECT CASE(TRIM(td_var%c_extrap(1)))
294	CASE('min_error')
295	cl_method='min_error'
296	CASE DEFAULT
297	cl_method='dist_weight'
298
299	!update variable structure
300	td_var%c_extrap(1)='dist_weight'
301	END SELECT
302
303	! number of point use to compute box
304	il_radius=1
305	IF( PRESENT(id_radius) ) il_radius=id_radius
306	IF( il_radius < 0 )THEN
307	CALL logger_error("EXTRAP FILL VALUE: invalid "//&
308	& " radius of the box used to compute extrapolation "//&
309	& "("//TRIM(fct_str(il_radius))//")")
310	ENDIF
311
312	CALL logger_info("EXTRAP FILL: extrapolate "//TRIM(td_var%c_name)//&
313	& " using "//TRIM(cl_method)//" method." )
314
315	CALL extrap__fill_value( td_var, cl_method, &
316	& il_radius )
317
318	ENDIF
319
320	END SUBROUTINE extrap__fill_value_wrapper
321	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
322	SUBROUTINE extrap__fill_value(td_var, cd_method, id_radius)
323	!-------------------------------------------------------------------
324	!> @brief
325	!> This subroutine compute point to be extrapolated, then extrapolate point.
326	!>
327	!> @details
328	!> optionaly, you could specify :<br/>
329	!> - refinment factor (default 1)
330	!> - offset between fine and coarse grid (default compute from refinment factor
331	!> as offset=(rho-1)/2)
332	!>
333	!> @author J.Paul
334	!> @date November, 2013 - Initial Version
335	!> @date June, 2015
336	!> - select all land points for extrapolation
337	!>
338	!> @param[inout] td_var variable structure
339	!> @param[in] cd_method extrapolation method
340	!> @param[in] id_radius radius of the halo used to compute extrapolation
341	!-------------------------------------------------------------------
342
343	IMPLICIT NONE
344
345	! Argument
346	TYPE(TVAR) , INTENT(INOUT) :: td_var
347	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
348	INTEGER(i4) , INTENT(IN ) :: id_radius
349
350	! local variable
351	CHARACTER(LEN=lc) :: cl_extrap
352
353	INTEGER(i4), DIMENSION(:,:,:) , ALLOCATABLE :: il_detect
354
355	TYPE(TATT) :: tl_att
356
357	! loop indices
358	!----------------------------------------------------------------
359
360	!1- detect point to be extrapolated
361	ALLOCATE( il_detect( td_var%t_dim(1)%i_len, &
362	& td_var%t_dim(2)%i_len, &
363	& td_var%t_dim(3)%i_len) )
364
365	il_detect(:,:,:) = extrap_detect( td_var )
366
367	!2- add attribute to variable
368	cl_extrap=fct_concat(td_var%c_extrap(:))
369	tl_att=att_init('extrapolation',cl_extrap)
370	CALL var_move_att(td_var, tl_att)
371
372	CALL att_clean(tl_att)
373
374	IF( ALL(il_detect(:,:,:)==1) )THEN
375	CALL logger_warn(" EXTRAP FILL: "//&
376	& " can not extrapolate "//TRIM(td_var%c_name)//&
377	& ". no value inform." )
378	ELSE
379	CALL logger_info(" EXTRAP FILL: "//&
380	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
381	& " point(s) to extrapolate " )
382
383	CALL logger_info(" EXTRAP FILL: method "//&
384	& TRIM(cd_method) )
385
386	!3- extrapolate
387	CALL extrap__3D(td_var%d_value(:,:,:,:), td_var%d_fill, &
388	& il_detect(:,:,:), &
389	& cd_method, id_radius )
390	ENDIF
391
392	DEALLOCATE(il_detect)
393
394	END SUBROUTINE extrap__fill_value
395	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
396	SUBROUTINE extrap__3D(dd_value, dd_fill, id_detect,&
397	& cd_method, id_radius)
398	!-------------------------------------------------------------------
399	!> @brief
400	!> This subroutine compute point to be extrapolated in 3D array.
401	!>
402	!> @details
403	!> in case of 'min_error' method:<br/>
404	!> - compute derivative in i-, j- and k- direction
405	!> - compute minimum error coefficient (distance to center of halo)
406	!> - compute extrapolatd values by calculated minimum error using taylor expansion
407	!> in case of 'dist_weight' method:<br/>
408	!> - compute distance weight coefficient (inverse of distance to center of halo)
409	!> - compute extrapolatd values using Inverse Distance Weighting
410	!>
411	!> @author J.Paul
412	!> @date November, 2013 - Initial Version
413	!> @date July, 2015
414	!> - compute coef indices to be used
415	!> - bug fix: force coef indice to 1, for dimension lenth equal to 1
416	!>
417	!> @param[inout] dd_value 3D array of variable to be extrapolated
418	!> @param[in] dd_fill FillValue of variable
419	!> @param[inout] id_detect array of point to extrapolate
420	!> @param[in] cd_method extrapolation method
421	!> @param[in] id_radius radius of the halo used to compute extrapolation
422	!-------------------------------------------------------------------
423
424	IMPLICIT NONE
425
426	! Argument
427	REAL(dp) , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value
428	REAL(dp) , INTENT(IN ) :: dd_fill
429	INTEGER(i4), DIMENSION(:,:,:) , INTENT(INOUT) :: id_detect
430	CHARACTER(LEN=*), INTENT(IN ) :: cd_method
431	INTEGER(i4), INTENT(IN ) :: id_radius
432
433	! local variable
434	INTEGER(i4) :: il_imin
435	INTEGER(i4) :: il_imax
436	INTEGER(i4) :: il_jmin
437	INTEGER(i4) :: il_jmax
438	INTEGER(i4) :: il_kmin
439	INTEGER(i4) :: il_kmax
440	INTEGER(i4) :: il_iter
441	INTEGER(i4) :: il_radius
442	INTEGER(i4) :: il_i1
443	INTEGER(i4) :: il_i2
444	INTEGER(i4) :: il_j1
445	INTEGER(i4) :: il_j2
446	INTEGER(i4) :: il_k1
447	INTEGER(i4) :: il_k2
448
449	INTEGER(i4), DIMENSION(4) :: il_shape
450	INTEGER(i4), DIMENSION(3) :: il_dim
451
452	INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect
453
454	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdx
455	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdy
456	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dfdz
457	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
458
459	LOGICAL :: ll_iter
460
461	! loop indices
462	INTEGER(i4) :: ji
463	INTEGER(i4) :: jj
464	INTEGER(i4) :: jk
465	INTEGER(i4) :: jl
466	!----------------------------------------------------------------
467
468	il_shape(:)=SHAPE(dd_value)
469
470	ALLOCATE( il_detect( il_shape(1), il_shape(2), il_shape(3)) )
471
472	SELECT CASE(TRIM(cd_method))
473	CASE('min_error')
474	DO jl=1,il_shape(4)
475
476	! initialise table of point to be extrapolated
477	il_detect(:,:,:)=id_detect(:,:,:)
478
479	il_iter=1
480	DO WHILE( ANY(il_detect(:,:,:)==1) )
481	! change extend value to minimize number of iteration
482	il_radius=id_radius+(il_iter-1)
483	ll_iter=.TRUE.
484
485	ALLOCATE( dl_dfdx(il_shape(1), il_shape(2), il_shape(3)) )
486	ALLOCATE( dl_dfdy(il_shape(1), il_shape(2), il_shape(3)) )
487	ALLOCATE( dl_dfdz(il_shape(1), il_shape(2), il_shape(3)) )
488
489	! compute derivative in i-direction
490	dl_dfdx(:,:,:)=dd_fill
491	IF( il_shape(1) > 1 )THEN
492	dl_dfdx(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
493	& dd_fill, 'I' )
494	ENDIF
495
496	! compute derivative in j-direction
497	dl_dfdy(:,:,:)=dd_fill
498	IF( il_shape(2) > 1 )THEN
499	dl_dfdy(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
500	& dd_fill, 'J' )
501	ENDIF
502
503	! compute derivative in k-direction
504	dl_dfdz(:,:,:)=dd_fill
505	IF( il_shape(3) > 1 )THEN
506	dl_dfdz(:,:,:)=math_deriv_3D( dd_value(:,:,:,jl), &
507	& dd_fill, 'K' )
508	ENDIF
509
510	il_dim(1)=2*il_radius+1
511	IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
512	il_dim(2)=2*il_radius+1
513	IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
514	il_dim(3)=2*il_radius+1
515	IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
516
517	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
518
519	dl_coef(:,:,:)=extrap__3D_min_error_coef(dd_value( 1:il_dim(1), &
520	& 1:il_dim(2), &
521	& 1:il_dim(3), &
522	& jl ))
523
524	DO jk=1,il_shape(3)
525	! from North West(1,1) to South East(il_shape(1),il_shape(2))
526	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
527	DO jj=1,il_shape(2)
528	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
529	DO ji=1,il_shape(1)
530
531	IF( il_detect(ji,jj,jk) == 1 )THEN
532
533	il_imin=MAX(ji-il_radius,1)
534	il_imax=MIN(ji+il_radius,il_shape(1))
535	! coef indices to be used
536	il_i1 = il_radius-(ji-il_imin)+1
537	il_i2 = il_radius+(il_imax-ji)+1
538	IF( il_dim(1) == 1 )THEN
539	il_imin=ji
540	il_imax=ji
541	! coef indices to be used
542	il_i1 = 1
543	il_i2 = 1
544	ENDIF
545
546
547	il_jmin=MAX(jj-il_radius,1)
548	il_jmax=MIN(jj+il_radius,il_shape(2))
549	! coef indices to be used
550	il_j1 = il_radius-(jj-il_jmin)+1
551	il_j2 = il_radius+(il_jmax-jj)+1
552	IF( il_dim(2) == 1 )THEN
553	il_jmin=jj
554	il_jmax=jj
555	! coef indices to be used
556	il_j1 = 1
557	il_j2 = 1
558	ENDIF
559
560	il_kmin=MAX(jk-il_radius,1)
561	il_kmax=MIN(jk+il_radius,il_shape(3))
562	! coef indices to be used
563	il_k1 = il_radius-(jk-il_kmin)+1
564	il_k2 = il_radius+(il_kmax-jk)+1
565	IF( il_dim(3) == 1 )THEN
566	il_kmin=jk
567	il_kmax=jk
568	! coef indices to be used
569	il_k1 = 1
570	il_k2 = 1
571	ENDIF
572
573	dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( &
574	& dd_value( il_imin:il_imax, &
575	& il_jmin:il_jmax, &
576	& il_kmin:il_kmax,jl ), dd_fill, il_radius, &
577	& dl_dfdx( il_imin:il_imax, &
578	& il_jmin:il_jmax, &
579	& il_kmin:il_kmax ), &
580	& dl_dfdy( il_imin:il_imax, &
581	& il_jmin:il_jmax, &
582	& il_kmin:il_kmax ), &
583	& dl_dfdz( il_imin:il_imax, &
584	& il_jmin:il_jmax, &
585	& il_kmin:il_kmax ), &
586	& dl_coef(il_i1:il_i2, &
587	& il_j1:il_j2, &
588	& il_k1:il_k2) )
589
590	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
591	il_detect(ji,jj,jk)= 0
592	ll_iter=.FALSE.
593	ENDIF
594
595	ENDIF
596
597	ENDDO
598	ENDDO
599	! from South East(il_shape(1),il_shape(2)) to North West(1,1)
600	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
601	DO jj=il_shape(2),1,-1
602	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
603	DO ji=il_shape(1),1,-1
604
605	IF( il_detect(ji,jj,jk) == 1 )THEN
606
607	il_imin=MAX(ji-il_radius,1)
608	il_imax=MIN(ji+il_radius,il_shape(1))
609	! coef indices to be used
610	il_i1 = il_radius-(ji-il_imin)+1
611	il_i2 = il_radius+(il_imax-ji)+1
612	IF( il_dim(1) == 1 )THEN
613	il_imin=ji
614	il_imax=ji
615	! coef indices to be used
616	il_i1 = 1
617	il_i2 = 1
618	ENDIF
619
620
621	il_jmin=MAX(jj-il_radius,1)
622	il_jmax=MIN(jj+il_radius,il_shape(2))
623	! coef indices to be used
624	il_j1 = il_radius-(jj-il_jmin)+1
625	il_j2 = il_radius+(il_jmax-jj)+1
626	IF( il_dim(2) == 1 )THEN
627	il_jmin=jj
628	il_jmax=jj
629	! coef indices to be used
630	il_j1 = 1
631	il_j2 = 1
632	ENDIF
633
634	il_kmin=MAX(jk-il_radius,1)
635	il_kmax=MIN(jk+il_radius,il_shape(3))
636	! coef indices to be used
637	il_k1 = il_radius-(jk-il_kmin)+1
638	il_k2 = il_radius+(il_kmax-jk)+1
639	IF( il_dim(3) == 1 )THEN
640	il_kmin=jk
641	il_kmax=jk
642	! coef indices to be used
643	il_k1 = 1
644	il_k2 = 1
645	ENDIF
646
647	dd_value(ji,jj,jk,jl)=extrap__3D_min_error_fill( &
648	& dd_value( il_imin:il_imax, &
649	& il_jmin:il_jmax, &
650	& il_kmin:il_kmax,jl ), dd_fill, il_radius, &
651	& dl_dfdx( il_imin:il_imax, &
652	& il_jmin:il_jmax, &
653	& il_kmin:il_kmax ), &
654	& dl_dfdy( il_imin:il_imax, &
655	& il_jmin:il_jmax, &
656	& il_kmin:il_kmax ), &
657	& dl_dfdz( il_imin:il_imax, &
658	& il_jmin:il_jmax, &
659	& il_kmin:il_kmax ), &
660	& dl_coef(il_i1:il_i2, &
661	& il_j1:il_j2, &
662	& il_k1:il_k2) )
663
664	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
665	il_detect(ji,jj,jk)= 0
666	ll_iter=.FALSE.
667	ENDIF
668
669	ENDIF
670
671	ENDDO
672	ENDDO
673	ENDDO
674
675	DEALLOCATE( dl_dfdx )
676	DEALLOCATE( dl_dfdy )
677	DEALLOCATE( dl_dfdz )
678	DEALLOCATE( dl_coef )
679
680	IF( ll_iter ) il_iter=il_iter+1
681	ENDDO
682	ENDDO
683
684	CASE DEFAULT ! 'dist_weight'
685	DO jl=1,il_shape(4)
686
687	! intitialise table of poitn to be extrapolated
688	il_detect(:,:,:)=id_detect(:,:,:)
689
690	il_iter=1
691	DO WHILE( ANY(il_detect(:,:,:)==1) )
692	! change extend value to minimize number of iteration
693	il_radius=id_radius+(il_iter-1)
694	ll_iter=.TRUE.
695
696	il_dim(1)=2*il_radius+1
697	IF( il_shape(1) < 2*il_radius+1 ) il_dim(1)=1
698	il_dim(2)=2*il_radius+1
699	IF( il_shape(2) < 2*il_radius+1 ) il_dim(2)=1
700	il_dim(3)=2*il_radius+1
701	IF( il_shape(3) < 2*il_radius+1 ) il_dim(3)=1
702
703	ALLOCATE( dl_coef(il_dim(1), il_dim(2), il_dim(3)) )
704
705	dl_coef(:,:,:)=extrap__3D_dist_weight_coef(dd_value(1:il_dim(1),&
706	& 1:il_dim(2),&
707	& 1:il_dim(3),&
708	& jl ) )
709
710	DO jk=1,il_shape(3)
711	! from North West(1,1) to South East(il_shape(1),il_shape(2))
712	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
713	DO jj=1,il_shape(2)
714	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
715	DO ji=1,il_shape(1)
716
717	IF( il_detect(ji,jj,jk) == 1 )THEN
718
719	il_imin=MAX(ji-il_radius,1)
720	il_imax=MIN(ji+il_radius,il_shape(1))
721	! coef indices to be used
722	il_i1 = il_radius-(ji-il_imin)+1
723	il_i2 = il_radius+(il_imax-ji)+1
724	IF( il_dim(1) == 1 )THEN
725	il_imin=ji
726	il_imax=ji
727	! coef indices to be used
728	il_i1 = 1
729	il_i2 = 1
730	ENDIF
731
732	il_jmin=MAX(jj-il_radius,1)
733	il_jmax=MIN(jj+il_radius,il_shape(2))
734	! coef indices to be used
735	il_j1 = il_radius-(jj-il_jmin)+1
736	il_j2 = il_radius+(il_jmax-jj)+1
737	IF( il_dim(2) == 1 )THEN
738	il_jmin=jj
739	il_jmax=jj
740	! coef indices to be used
741	il_j1 = 1
742	il_j2 = 1
743	ENDIF
744
745	il_kmin=MAX(jk-il_radius,1)
746	il_kmax=MIN(jk+il_radius,il_shape(3))
747	! coef indices to be used
748	il_k1 = il_radius-(jk-il_kmin)+1
749	il_k2 = il_radius+(il_kmax-jk)+1
750	IF( il_dim(3) == 1 )THEN
751	il_kmin=jk
752	il_kmax=jk
753	! coef indices to be used
754	il_k1 = 1
755	il_k2 = 1
756	ENDIF
757
758	dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( &
759	& dd_value( il_imin:il_imax, &
760	& il_jmin:il_jmax, &
761	& il_kmin:il_kmax, &
762	& jl), dd_fill, il_radius, &
763	& dl_coef(il_i1:il_i2, &
764	& il_j1:il_j2, &
765	& il_k1:il_k2) )
766
767	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
768	il_detect(ji,jj,jk)= 0
769	ll_iter=.FALSE.
770	ENDIF
771
772	ENDIF
773
774	ENDDO
775	ENDDO
776	! from South East(il_shape(1),il_shape(2)) to North West(1,1)
777	IF( ALL(il_detect(:,:,jk) == 0) ) CYCLE
778	DO jj=il_shape(2),1,-1
779	IF( ALL(il_detect(:,jj,jk) == 0) ) CYCLE
780	DO ji=il_shape(1),1,-1
781
782	IF( il_detect(ji,jj,jk) == 1 )THEN
783
784	il_imin=MAX(ji-il_radius,1)
785	il_imax=MIN(ji+il_radius,il_shape(1))
786	! coef indices to be used
787	il_i1 = il_radius-(ji-il_imin)+1
788	il_i2 = il_radius+(il_imax-ji)+1
789	IF( il_dim(1) == 1 )THEN
790	il_imin=ji
791	il_imax=ji
792	! coef indices to be used
793	il_i1 = 1
794	il_i2 = 1
795	ENDIF
796
797	il_jmin=MAX(jj-il_radius,1)
798	il_jmax=MIN(jj+il_radius,il_shape(2))
799	! coef indices to be used
800	il_j1 = il_radius-(jj-il_jmin)+1
801	il_j2 = il_radius+(il_jmax-jj)+1
802	IF( il_dim(2) == 1 )THEN
803	il_jmin=jj
804	il_jmax=jj
805	! coef indices to be used
806	il_j1 = 1
807	il_j2 = 1
808	ENDIF
809
810	il_kmin=MAX(jk-il_radius,1)
811	il_kmax=MIN(jk+il_radius,il_shape(3))
812	! coef indices to be used
813	il_k1 = il_radius-(jk-il_kmin)+1
814	il_k2 = il_radius+(il_kmax-jk)+1
815	IF( il_dim(3) == 1 )THEN
816	il_kmin=jk
817	il_kmax=jk
818	! coef indices to be used
819	il_k1 = 1
820	il_k2 = 1
821	ENDIF
822
823	dd_value(ji,jj,jk,jl)=extrap__3D_dist_weight_fill( &
824	& dd_value( il_imin:il_imax, &
825	& il_jmin:il_jmax, &
826	& il_kmin:il_kmax, &
827	& jl), dd_fill, il_radius, &
828	& dl_coef(il_i1:il_i2, &
829	& il_j1:il_j2, &
830	& il_k1:il_k2) )
831
832	IF( dd_value(ji,jj,jk,jl) /= dd_fill )THEN
833	il_detect(ji,jj,jk)= 0
834	ll_iter=.FALSE.
835	ENDIF
836
837	ENDIF
838
839	ENDDO
840	ENDDO
841	ENDDO
842	CALL logger_info(" EXTRAP 3D: "//&
843	& TRIM(fct_str(SUM(il_detect(:,:,:))))//&
844	& " point(s) to extrapolate " )
845
846	DEALLOCATE( dl_coef )
847	IF( ll_iter ) il_iter=il_iter+1
848	ENDDO
849	ENDDO
850	END SELECT
851
852	DEALLOCATE( il_detect )
853
854	END SUBROUTINE extrap__3D
855	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
856	PURE FUNCTION extrap__3D_min_error_coef(dd_value) &
857	& RESULT (df_value)
858	!-------------------------------------------------------------------
859	!> @brief
860	!> This function compute coefficient for min_error extrapolation.
861	!>
862	!> @details
863	!> coefficients are "grid distance" to the center of the box
864	!> choosed to compute extrapolation.
865	!>
866	!> @author J.Paul
867	!> @date November, 2013 - Initial Version
868	!> @date July, 2015
869	!> - decrease weight of third dimension
870	!>
871	!> @param[in] dd_value 3D array of variable to be extrapolated
872	!> @return 3D array of coefficient for minimum error extrapolation
873	!-------------------------------------------------------------------
874
875	IMPLICIT NONE
876
877	! Argument
878	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
879
880	! function
881	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
882	& SIZE(dd_value(:,:,:),DIM=2), &
883	& SIZE(dd_value(:,:,:),DIM=3) ) :: df_value
884
885	! local variable
886	INTEGER(i4), DIMENSION(3) :: il_shape
887
888	INTEGER(i4) :: il_imid
889	INTEGER(i4) :: il_jmid
890	INTEGER(i4) :: il_kmid
891
892	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
893
894	! loop indices
895	INTEGER(i4) :: ji
896	INTEGER(i4) :: jj
897	INTEGER(i4) :: jk
898	!----------------------------------------------------------------
899
900	! init
901	df_value(:,:,:)=0
902
903	il_shape(:)=SHAPE(dd_value(:,:,:))
904
905	il_imid=INT(REAL(il_shape(1),sp)*0.5+1)
906	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1)
907	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1)
908
909	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
910
911	DO jk=1,il_shape(3)
912	DO jj=1,il_shape(2)
913	DO ji=1,il_shape(1)
914
915	! compute distance
916	! "vertical weight" is lower than horizontal
917	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
918	& (jj-il_jmid)**2 + &
919	& 3(jk-il_kmid)*2
920
921	IF( dl_dist(ji,jj,jk) /= 0 )THEN
922	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
923	ENDIF
924
925	ENDDO
926	ENDDO
927	ENDDO
928
929	WHERE( dl_dist(:,:,:) /= 0 )
930	df_value(:,:,:)=dl_dist(:,:,:)
931	END WHERE
932
933	DEALLOCATE( dl_dist )
934
935	END FUNCTION extrap__3D_min_error_coef
936	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
937	PURE FUNCTION extrap__3D_min_error_fill(dd_value, dd_fill, id_radius,&
938	& dd_dfdx, dd_dfdy, dd_dfdz, &
939	& dd_coef) &
940	& RESULT (df_value)
941	!-------------------------------------------------------------------
942	!> @brief
943	!> This function compute extrapolatd value by calculated minimum error using
944	!> taylor expansion
945	!>
946	!> @author J.Paul
947	!> @date November, 2013 - Initial Version
948	!>
949	!> @param[in] dd_value 3D array of variable to be extrapolated
950	!> @param[in] dd_fill FillValue of variable
951	!> @param[in] id_radius radius of the halo used to compute extrapolation
952	!> @param[in] dd_dfdx derivative of function in i-direction
953	!> @param[in] dd_dfdy derivative of function in j-direction
954	!> @param[in] dd_dfdz derivative of function in k-direction
955	!> @param[in] dd_coef array of coefficient for min_error extrapolation
956	!> @return extrapolatd value
957	!-------------------------------------------------------------------
958
959	IMPLICIT NONE
960
961	! Argument
962	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
963	REAL(dp) , INTENT(IN) :: dd_fill
964	INTEGER(i4), INTENT(IN) :: id_radius
965	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdx
966	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdy
967	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_dfdz
968	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
969
970	! function
971	REAL(dp) :: df_value
972
973	! local variable
974	INTEGER(i4), DIMENSION(3) :: il_shape
975	INTEGER(i4), DIMENSION(3) :: il_ind
976
977	INTEGER(i4), DIMENSION(:,:,:), ALLOCATABLE :: il_mask
978	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_error
979
980	INTEGER(i4) :: il_min
981	! loop indices
982
983	!----------------------------------------------------------------
984
985	! init
986	df_value=dd_fill
987
988	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
989
990	IF( COUNT(dd_value(:,:,:) /= dd_fill) >= il_min )THEN
991
992	il_shape(:)=SHAPE(dd_value(:,:,:))
993	ALLOCATE( il_mask( il_shape(1),il_shape(2),il_shape(3)) )
994	ALLOCATE( dl_error(il_shape(1),il_shape(2),il_shape(3)) )
995
996	! compute error
997	dl_error(:,:,:)=0.
998	il_mask(:,:,:)=0
999	WHERE( dd_dfdx(:,:,:) /= dd_fill )
1000	dl_error(:,:,:)=dd_coef(:,:,:)*dd_dfdx(:,:,:)
1001	il_mask(:,:,:)=1
1002	END WHERE
1003	WHERE( dd_dfdy(:,:,:) /= dd_fill )
1004	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdy(:,:,:))
1005	il_mask(:,:,:)=1
1006	END WHERE
1007	WHERE( dd_dfdz(:,:,:) /= dd_fill )
1008	dl_error(:,:,:)=(dl_error(:,:,:)+dd_coef(:,:,:)*dd_dfdz(:,:,:))
1009	il_mask(:,:,:)=1
1010	END WHERE
1011
1012	! get minimum error indices
1013	il_ind(:)=MINLOC(dl_error(:,:,:),il_mask(:,:,:)==1)
1014
1015	! return value
1016	IF( ALL(il_ind(:)/=0) )THEN
1017	df_value=dd_value(il_ind(1),il_ind(2),il_ind(3))
1018	ENDIF
1019
1020	DEALLOCATE( il_mask )
1021	DEALLOCATE( dl_error )
1022
1023	ENDIF
1024
1025	END FUNCTION extrap__3D_min_error_fill
1026	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1027	PURE FUNCTION extrap__3D_dist_weight_coef(dd_value) &
1028	& RESULT (df_value)
1029	!-------------------------------------------------------------------
1030	!> @brief
1031	!> This function compute coefficient for inverse distance weighted method
1032	!>
1033	!> @details
1034	!> coefficients are inverse "grid distance" to the center of the box choosed to compute
1035	!> extrapolation.
1036	!>
1037	!> @author J.Paul
1038	!> @date November, 2013 - Initial Version
1039	!> @date July, 2015
1040	!> - decrease weight of third dimension
1041	!>
1042	!> @param[in] dd_value 3D array of variable to be extrapolated
1043	!> @return 3D array of coefficient for inverse distance weighted extrapolation
1044	!-------------------------------------------------------------------
1045
1046	IMPLICIT NONE
1047
1048	! Argument
1049	REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: dd_value
1050
1051	! function
1052	REAL(dp), DIMENSION(SIZE(dd_value(:,:,:),DIM=1), &
1053	& SIZE(dd_value(:,:,:),DIM=2), &
1054	& SIZE(dd_value(:,:,:),DIM=3) ) :: df_value
1055
1056	! local variable
1057	INTEGER(i4), DIMENSION(3) :: il_shape
1058
1059	INTEGER(i4) :: il_imid
1060	INTEGER(i4) :: il_jmid
1061	INTEGER(i4) :: il_kmid
1062
1063	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_dist
1064
1065	! loop indices
1066	INTEGER(i4) :: ji
1067	INTEGER(i4) :: jj
1068	INTEGER(i4) :: jk
1069	!----------------------------------------------------------------
1070
1071	! init
1072	df_value(:,:,:)=0
1073
1074	il_shape(:)=SHAPE(dd_value(:,:,:))
1075
1076	il_imid=INT(REAL(il_shape(1),sp)*0.5+1,i4)
1077	il_jmid=INT(REAL(il_shape(2),sp)*0.5+1,i4)
1078	il_kmid=INT(REAL(il_shape(3),sp)*0.5+1,i4)
1079
1080	ALLOCATE( dl_dist(il_shape(1),il_shape(2),il_shape(3)) )
1081
1082	DO jk=1,il_shape(3)
1083	DO jj=1,il_shape(2)
1084	DO ji=1,il_shape(1)
1085
1086	! compute distance
1087	! "vertical weight" is lower than horizontal
1088	dl_dist(ji,jj,jk) = (ji-il_imid)**2 + &
1089	& (jj-il_jmid)**2 + &
1090	& 3(jk-il_kmid)*2
1091
1092	IF( dl_dist(ji,jj,jk) /= 0 )THEN
1093	dl_dist(ji,jj,jk)=SQRT( dl_dist(ji,jj,jk) )
1094	ENDIF
1095
1096	ENDDO
1097	ENDDO
1098	ENDDO
1099
1100	WHERE( dl_dist(:,:,:) /= 0 )
1101	df_value(:,:,:)=1./dl_dist(:,:,:)
1102	END WHERE
1103
1104	DEALLOCATE( dl_dist )
1105
1106	END FUNCTION extrap__3D_dist_weight_coef
1107	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1108	FUNCTION extrap__3D_dist_weight_fill(dd_value, dd_fill, id_radius, &
1109	& dd_coef) &
1110	& RESULT (df_value)
1111	!-------------------------------------------------------------------
1112	!> @brief
1113	!> This function compute extrapolatd value using inverse distance weighted
1114	!> method
1115	!>
1116	!> @details
1117	!>
1118	!> @author J.Paul
1119	!> @date November, 2013 - Initial Version
1120	!>
1121	!> @param[in] dd_value 3D array of variable to be extrapolated
1122	!> @param[in] dd_fill FillValue of variable
1123	!> @param[in] id_radius radius of the halo used to compute extrapolation
1124	!> @param[in] dd_coef 3D array of coefficient for inverse distance weighted extrapolation
1125	!> @return extrapolatd value
1126	!-------------------------------------------------------------------
1127
1128	IMPLICIT NONE
1129
1130	! Argument
1131	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_value
1132	REAL(dp) , INTENT(IN) :: dd_fill
1133	INTEGER(i4), INTENT(IN) :: id_radius
1134	REAL(dp) , DIMENSION(:,:,:), INTENT(IN) :: dd_coef
1135
1136	! function
1137	REAL(dp) :: df_value
1138
1139	! local variable
1140	INTEGER(i4), DIMENSION(3) :: il_shape
1141
1142	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_value
1143	REAL(dp) , DIMENSION(:,:,:), ALLOCATABLE :: dl_coef
1144
1145	INTEGER(i4) :: il_min
1146	! loop indices
1147	INTEGER(i4) :: ji
1148	INTEGER(i4) :: jj
1149	INTEGER(i4) :: jk
1150	!----------------------------------------------------------------
1151
1152	! init
1153	df_value=dd_fill
1154
1155	il_min=MAX(1,(SIZE(dd_value(:,:,:)))/(1+id_radius*2))
1156
1157	IF( COUNT(dd_value(:,:,:)/= dd_fill) >= il_min )THEN
1158
1159	il_shape(:)=SHAPE(dd_value(:,:,:))
1160	ALLOCATE( dl_value( il_shape(1),il_shape(2),il_shape(3)) )
1161	ALLOCATE( dl_coef( il_shape(1),il_shape(2),il_shape(3)) )
1162
1163	dl_value(:,:,:)=0
1164	dl_coef(:,:,:)=0
1165
1166	DO jk=1,il_shape(3)
1167	DO jj=1,il_shape(2)
1168	DO ji=1,il_shape(1)
1169
1170	IF( dd_value(ji,jj,jk) /= dd_fill )THEN
1171	! compute factor
1172	dl_value(ji,jj,jk)=dd_coef(ji,jj,jk)*dd_value(ji,jj,jk)
1173	dl_coef(ji,jj,jk)=dd_coef(ji,jj,jk)
1174	ENDIF
1175
1176	ENDDO
1177	ENDDO
1178	ENDDO
1179
1180
1181	! return value
1182	IF( SUM( dl_coef(:,:,:) ) /= 0 )THEN
1183	df_value = SUM( dl_value(:,:,:) )/SUM( dl_coef(:,:,:) )
1184	ENDIF
1185
1186	DEALLOCATE( dl_value )
1187	DEALLOCATE( dl_coef )
1188
1189	ENDIF
1190
1191	END FUNCTION extrap__3D_dist_weight_fill
1192	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1193	SUBROUTINE extrap_add_extrabands(td_var, id_isize, id_jsize)
1194	!-------------------------------------------------------------------
1195	!> @brief
1196	!> This subroutine add to the variable (to be extrapolated) an
1197	!> extraband of N points at north,south,east and west boundaries.
1198	!>
1199	!> @details
1200	!> optionaly you could specify size of extra bands in i- and j-direction
1201	!>
1202	!> @author J.Paul
1203	!> @date November, 2013 - Initial version
1204	!>
1205	!> @param[inout] td_var variable
1206	!> @param[in] id_isize i-direction size of extra bands (default=im_minext)
1207	!> @param[in] id_jsize j-direction size of extra bands (default=im_minext)
1208	!> @todo
1209	!> - invalid special case for grid with north fold
1210	!-------------------------------------------------------------------
1211
1212	IMPLICIT NONE
1213
1214	! Argument
1215	TYPE(TVAR) , INTENT(INOUT) :: td_var
1216	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
1217	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
1218
1219	! local variable
1220	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
1221
1222	INTEGER(i4) :: il_isize
1223	INTEGER(i4) :: il_jsize
1224	INTEGER(i4) :: il_tmp
1225
1226	! loop indices
1227	INTEGER(i4) :: ji
1228	INTEGER(i4) :: ij
1229	!----------------------------------------------------------------
1230	il_isize=im_minext
1231	IF(PRESENT(id_isize)) il_isize=id_isize
1232	IF( il_isize < im_minext .AND. &
1233	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
1234	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
1235	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
1236	& " cubic interpolation ")
1237	ENDIF
1238
1239	il_jsize=im_minext
1240	IF(PRESENT(id_jsize)) il_jsize=id_jsize
1241	IF( il_jsize < im_minext .AND. &
1242	& TRIM(td_var%c_interp(1)) == 'cubic' )THEN
1243	CALL logger_warn("EXTRAP ADD EXTRABANDS: size of extrabands "//&
1244	& "should be at least "//TRIM(fct_str(im_minext))//" for "//&
1245	& " cubic interpolation ")
1246	ENDIF
1247
1248	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
1249	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
1250
1251	CALL logger_trace( "EXTRAP ADD EXTRABANDS: dimension change "//&
1252	& "in variable "//TRIM(td_var%c_name) )
1253
1254	! add extrabands in variable
1255	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
1256	& td_var%t_dim(2)%i_len, &
1257	& td_var%t_dim(3)%i_len, &
1258	& td_var%t_dim(4)%i_len ))
1259
1260	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
1261
1262
1263	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len + 2*il_isize
1264	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len + 2*il_jsize
1265
1266	DEALLOCATE(td_var%d_value)
1267	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
1268	& td_var%t_dim(2)%i_len, &
1269	& td_var%t_dim(3)%i_len, &
1270	& td_var%t_dim(4)%i_len ) )
1271
1272	! intialise
1273	td_var%d_value(:,:,:,:)=td_var%d_fill
1274
1275	! fill center
1276	td_var%d_value( 1+il_isize:td_var%t_dim(1)%i_len-il_isize, &
1277	& 1+il_jsize:td_var%t_dim(2)%i_len-il_jsize, &
1278	& :,:) = dl_value(:,:,:,:)
1279
1280	! special case for overlap
1281	IF( td_var%i_ew >= 0 .AND. il_isize /= 0 )THEN
1282	DO ji=1,il_isize
1283	! from east to west
1284	il_tmp=td_var%t_dim(1)%i_len-td_var%i_ew+ji-2*il_isize
1285	td_var%d_value(ji,:,:,:) = td_var%d_value(il_tmp,:,:,:)
1286
1287	! from west to east
1288	ij=td_var%t_dim(1)%i_len-ji+1
1289	il_tmp=td_var%i_ew-ji+2*il_isize+1
1290	td_var%d_value(ij,:,:,:) = td_var%d_value(il_tmp,:,:,:)
1291	ENDDO
1292	ENDIF
1293
1294	DEALLOCATE( dl_value )
1295
1296	END SUBROUTINE extrap_add_extrabands
1297	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1298	SUBROUTINE extrap_del_extrabands(td_var, id_isize, id_jsize)
1299	!-------------------------------------------------------------------
1300	!> @brief
1301	!> This subroutine remove of the variable an extraband
1302	!> of N points at north,south,east and west boundaries.
1303	!>
1304	!> @details
1305	!> optionaly you could specify size of extra bands in i- and j-direction
1306	!>
1307	!> @author J.Paul
1308	!> @date November, 2013 - Initial version
1309	!>
1310	!> @param[inout] td_var variable
1311	!> @param[in] id_isize i-direction size of extra bands (default=im_minext)
1312	!> @param[in] id_jsize j-direction size of extra bands (default=im_minext)
1313	!-------------------------------------------------------------------
1314
1315	IMPLICIT NONE
1316
1317	! Argument
1318	TYPE(TVAR) , INTENT(INOUT) :: td_var
1319	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_isize
1320	INTEGER(i4), INTENT(IN ), OPTIONAL :: id_jsize
1321
1322	! local variable
1323	REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
1324
1325	INTEGER(i4) :: il_isize
1326	INTEGER(i4) :: il_jsize
1327
1328	INTEGER(i4) :: il_imin
1329	INTEGER(i4) :: il_imax
1330	INTEGER(i4) :: il_jmin
1331	INTEGER(i4) :: il_jmax
1332
1333	! loop indices
1334	!----------------------------------------------------------------
1335	il_isize=im_minext
1336	IF(PRESENT(id_isize)) il_isize=id_isize
1337
1338	il_jsize=im_minext
1339	IF(PRESENT(id_jsize)) il_jsize=id_jsize
1340
1341	IF( .NOT. td_var%t_dim(1)%l_use ) il_isize=0
1342	IF( .NOT. td_var%t_dim(2)%l_use ) il_jsize=0
1343
1344	CALL logger_trace( "EXTRAP DEL EXTRABANDS: dimension change "//&
1345	& "in variable "//TRIM(td_var%c_name) )
1346
1347	! add extrabands in variable
1348	ALLOCATE(dl_value( td_var%t_dim(1)%i_len, &
1349	& td_var%t_dim(2)%i_len, &
1350	& td_var%t_dim(3)%i_len, &
1351	& td_var%t_dim(4)%i_len ))
1352
1353	dl_value(:,:,:,:)=td_var%d_value(:,:,:,:)
1354
1355	! fill center
1356	il_imin=1+il_isize
1357	il_imax=td_var%t_dim(1)%i_len-il_isize
1358
1359	il_jmin=1+il_jsize
1360	il_jmax=td_var%t_dim(2)%i_len-il_jsize
1361
1362	td_var%t_dim(1)%i_len = td_var%t_dim(1)%i_len - 2*il_isize
1363	td_var%t_dim(2)%i_len = td_var%t_dim(2)%i_len - 2*il_jsize
1364
1365	DEALLOCATE(td_var%d_value)
1366	ALLOCATE( td_var%d_value(td_var%t_dim(1)%i_len, &
1367	& td_var%t_dim(2)%i_len, &
1368	& td_var%t_dim(3)%i_len, &
1369	& td_var%t_dim(4)%i_len ) )
1370
1371	! intialise
1372	td_var%d_value(:,:,:,:)=td_var%d_fill
1373
1374	td_var%d_value(:,:,:,:)=dl_value(il_imin:il_imax,&
1375	& il_jmin:il_jmax,&
1376	& :,:)
1377
1378	DEALLOCATE( dl_value )
1379
1380	END SUBROUTINE extrap_del_extrabands
1381	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1382	END MODULE extrap

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

Download in other formats: