source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/TOOLS/SIREN/src/extrap.f90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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