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

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months 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!----------------------------------------------------------------------
72MODULE 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   
120CONTAINS
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   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1382END MODULE extrap
Note: See TracBrowser for help on using the repository browser.