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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
extrap.f90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/TOOLS/SIREN/src/extrap.f90 @ 6487

Last change on this file since 6487 was 6487, checked in by davestorkey, 8 years ago

Changes from nemo_v3_6_STABLE_copy branch.
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6237 cf. r5781 of /branches/UKMO/nemo_v3_6_STABLE_copy/NEMOGCM@6486

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