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

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

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

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago

update nemo trunk

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