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.
interp_cubic.f90 in utils/tools/SIREN/src – NEMO

source: utils/tools/SIREN/src/interp_cubic.f90

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

update: cf changelog inside documentation
File size: 30.6 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! DESCRIPTION:
6!> @brief
7!> This module manage cubic interpolation on regular grid.
8!>
9!>
10!> @details
11!> to compute cubic interpolation:<br/>
12!> @code
13!> CALL interp_cubic_fill(dd_value, dd_fill, id_detect, id_rho, ld_even [,ld_discont] )
14!> @endcode
15!>    - dd_value is 2D array of variable value
16!>    - dd_fill is the FillValue of variable
17!>    - id_detect is 2D array of point to be interpolated (see interp module)
18!>    - id_rho  is array of refinment factor
19!>    - ld_even indicates even refinment or not
20!>    - ld_discont indicates longitudinal discontinuity (-180°/180°, 0°/360°) or not
21!>
22!> @author
23!> J.Paul
24!>
25!> @date September, 2014 -Initial version
26!> @date June, 2015
27!> - use math module
28!>
29!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
30!----------------------------------------------------------------------
31MODULE interp_cubic
32
33   USE netcdf                          ! nf90 library
34   USE global                          ! global variable
35   USE kind                            ! F90 kind parameter
36   USE logger                          ! log file manager
37   USE fct                             ! basic useful function
38   USE math                            ! mathematical function
39
40   IMPLICIT NONE
41   ! NOTE_avoid_public_variables_if_possible
42
43   ! type and variable
44
45   ! function and subroutine
46   PUBLIC :: interp_cubic_fill  !< compute interpolation using cubic method
47
48   PRIVATE :: interp_cubic__2D           !< compute bicubic interpolation on 2D gid
49   PRIVATE :: interp_cubic__1D           !< compute   cubic interpolation on 1D gid
50   PRIVATE :: interp_cubic__2D_coef      !< compute coefficient for bicubic interpolation
51   PRIVATE :: interp_cubic__2D_fill      !< fill value using bicubic interpolation
52   PRIVATE :: interp_cubic__1D_coef      !< compute coefficient for   cubic interpolation
53   PRIVATE :: interp_cubic__1D_fill      !< fill value using   cubic interpolation
54   PRIVATE :: interp_cubic__get_weight2D !< compute interpoaltion weight for 2D array
55   PRIVATE :: interp_cubic__get_weight1D !< compute interpoaltion weight for 1D array
56
57CONTAINS
58   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59   SUBROUTINE interp_cubic_fill(dd_value, dd_fill, &
60         &                      id_detect,         &
61         &                      id_rho,            &
62         &                      ld_even, ld_discont)
63   !-------------------------------------------------------------------
64   !> @brief
65   !> This subroutine compute horizontal cubic interpolation on 4D array of value.
66   !>
67   !> @author J.Paul
68   !> @date September, 2014 - Initial Version
69   !> @date July, 2015
70   !> - reinitialise detect array for each level
71   !>
72   !> @param[inout] dd_value  2D array of variable value
73   !> @param[in] dd_fill      FillValue of variable
74   !> @param[inout] id_detect 2D array of point to be interpolated
75   !> @param[in] id_rho       array of refinment factor
76   !> @param[in] ld_even      even refinment or not
77   !> @param[in] ld_discont   longitudinal discontinuity (-180°/180°, 0°/360°) or not
78   !-------------------------------------------------------------------
79
80      IMPLICIT NONE
81
82      ! Argument
83      REAL(dp)        , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value
84      REAL(dp)                            , INTENT(IN   ) :: dd_fill
85      INTEGER(I4)     , DIMENSION(:,:,:)  , INTENT(INOUT) :: id_detect
86      INTEGER(I4)     , DIMENSION(:)      , INTENT(IN   ) :: id_rho
87      LOGICAL         , DIMENSION(:)      , INTENT(IN   ) :: ld_even
88      LOGICAL                             , INTENT(IN   ), OPTIONAL :: ld_discont
89
90      ! local variable
91      INTEGER(i4), DIMENSION(4)                  :: il_shape
92
93      INTEGER(I4), DIMENSION(:,:,:), ALLOCATABLE :: il_detect
94
95      LOGICAL                                    :: ll_discont
96
97      REAL(dp)   , DIMENSION(:,:)  , ALLOCATABLE :: dl_weight_IJ
98      REAL(dp)   , DIMENSION(:,:)  , ALLOCATABLE :: dl_weight_I
99      REAL(dp)   , DIMENSION(:,:)  , ALLOCATABLE :: dl_weight_J
100
101      ! loop indices
102      INTEGER(i4) :: ji
103      INTEGER(i4) :: jj
104      INTEGER(i4) :: jk
105      INTEGER(i4) :: jl
106      !----------------------------------------------------------------
107      ll_discont=.FALSE.
108      IF( PRESENT(ld_discont) ) ll_discont=ld_discont
109
110      il_shape(:)=SHAPE(dd_value)
111
112      ! compute vect2D
113      ALLOCATE(dl_weight_IJ(16,((id_rho(jp_I)+1)*(id_rho(jp_J)+1))) )
114      CALL interp_cubic__get_weight2D(dl_weight_IJ(:,:), &
115      &                               id_rho(:), ld_even(:))
116
117      ALLOCATE( dl_weight_I( 4,((id_rho(jp_I)+1)                 )) )
118      ALLOCATE( dl_weight_J( 4,(                 (id_rho(jp_J)+1))) )
119      CALL interp_cubic__get_weight1D(dl_weight_I(:,:), &
120      &                               id_rho(jp_I), ld_even(jp_I))
121      CALL interp_cubic__get_weight1D(dl_weight_J(:,:), &
122      &                               id_rho(jp_J), ld_even(jp_J))
123
124      ALLOCATE(il_detect(il_shape(1),il_shape(2),il_shape(3)))
125
126      DO jl=1,il_shape(4)
127         il_detect(:,:,:)=id_detect(:,:,:)
128         ! loop on vertical level
129         DO jk=1,il_shape(3)
130
131            ! I-J plan
132            CALL interp_cubic__2D(dd_value(:,:,jk,jl), dd_fill, &
133            &                     il_detect(:,:,jk),            &
134            &                     dl_weight_IJ(:,:),            &
135            &                     id_rho(jp_I), id_rho(jp_J),   &
136            &                     ll_discont)
137            IF( ANY(il_detect(:,:,jk)==1) )THEN
138               ! I direction
139               DO jj=1,il_shape(2)
140                  CALL interp_cubic__1D( dd_value(:,jj,jk,jl), dd_fill, &
141                  &                      il_detect(:,jj,jk),            &
142                  &                      dl_weight_I(:,:),              &
143                  &                      id_rho(jp_I), ll_discont )
144               ENDDO
145               IF( ALL(il_detect(:,:,jk)==0) )THEN
146                  CYCLE
147               ELSE
148                  ! J direction
149                  DO ji=1,il_shape(1)
150                     CALL interp_cubic__1D( dd_value(ji,:,jk,jl), dd_fill, &
151                     &                      il_detect(ji,:,jk),            &
152                     &                      dl_weight_J(:,:),              &
153                     &                      id_rho(jp_J), ll_discont )
154                  ENDDO
155               ENDIF
156            ENDIF
157
158         ENDDO
159      ENDDO
160
161      id_detect(:,:,:)=il_detect(:,:,:)
162      DEALLOCATE(il_detect)
163
164      DEALLOCATE(dl_weight_IJ)
165      DEALLOCATE(dl_weight_I)
166      DEALLOCATE(dl_weight_J)
167
168   END SUBROUTINE interp_cubic_fill
169   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
170   SUBROUTINE interp_cubic__2D(dd_value,  dd_fill,   &
171         &                     id_detect,            &
172         &                     dd_weight,            &
173         &                     id_rhoi, id_rhoj,     &
174         &                     ld_discont)
175   !-------------------------------------------------------------------
176   !> @brief
177   !> This subroutine compute cubic interpolation on 2D array of value.
178   !>
179   !> @details
180   !>
181   !> @author J.Paul
182   !> @date September, 2014 - Initial Version
183   !>
184   !> @param[inout] dd_value  2D array of variable value
185   !> @param[in] dd_fill      FillValue of variable
186   !> @param[inout] id_detect 2D array of point to be interpolated
187   !> @param[in] id_rhoi      refinment factor in i-direction
188   !> @param[in] id_rhoj      refinment factor in j-direction
189   !> @param[in] id_rhok      refinment factor in k-direction
190   !> @param[in] ld_even      even refinment or not
191   !> @param[in] ld_discont   longitudinal discontinuity (-180°/180°, 0°/360°) or not
192   !-------------------------------------------------------------------
193
194      IMPLICIT NONE
195
196      ! Argument
197      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value
198      REAL(dp)                        , INTENT(IN   ) :: dd_fill
199      INTEGER(I4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
200      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
201      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
202      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoj
203      LOGICAL                         , INTENT(IN   ) :: ld_discont
204
205      ! local variable
206      INTEGER(I4)                              :: il_xextra
207      INTEGER(I4)                              :: il_yextra
208      INTEGER(i4), DIMENSION(2)                :: il_shape
209      INTEGER(i4), DIMENSION(2)                :: il_dim
210
211      REAL(dp)                                 :: dl_min
212      REAL(dp)                                 :: dl_max
213      REAL(dp)   , DIMENSION(:)  , ALLOCATABLE :: dl_coef
214      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_coarse
215      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_tmp
216      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdx
217      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_dfdy
218      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_d2fdxy
219
220      ! loop indices
221      INTEGER(i4) :: ji
222      INTEGER(i4) :: jj
223      INTEGER(i4) :: ii
224      INTEGER(i4) :: ij
225      !----------------------------------------------------------------
226
227      IF( ANY(id_detect(:,:)==1) )THEN
228         il_shape(:)=SHAPE(dd_value)
229
230         ! compute coarse grid dimension
231         il_xextra=id_rhoi-1
232         il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
233
234         il_yextra=id_rhoj-1
235         il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj
236
237         ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )
238
239         ! value on coarse grid
240         dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &
241         &                        1:il_shape(2):id_rhoj )
242
243         ALLOCATE( dl_dfdx(  il_dim(1),il_dim(2)), &
244         &         dl_dfdy(  il_dim(1),il_dim(2)), &
245         &         dl_d2fdxy(il_dim(1),il_dim(2)) )
246
247         ! compute derivative on coarse grid
248         dl_dfdx(:,:)=math_deriv_2D(dl_coarse(:,:), dd_fill, 'I', ld_discont)
249         dl_dfdy(:,:)=math_deriv_2D(dl_coarse(:,:), dd_fill, 'J', ld_discont)
250
251         ! compute cross derivative on coarse grid
252         dl_d2fdxy(:,:)=math_deriv_2D(dl_dfdx(:,:), dd_fill, 'J', ld_discont)
253
254         ALLOCATE( dl_tmp(2,2) )
255         ALLOCATE( dl_coef(16) )
256
257         DO jj=1,il_shape(2)-1,id_rhoj
258            ij=((jj-1)/id_rhoj)+1
259            DO ji=1,il_shape(1)-1,id_rhoi
260               ii=((ji-1)/id_rhoi)+1
261
262               ! check if point to be interpolated
263               IF( ALL(id_detect(ji:ji+id_rhoi,   &
264               &                 jj:jj+id_rhoj)==0) ) CYCLE
265               ! check data needed to interpolate
266               IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
267               &   ANY(  dl_dfdx(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
268               &   ANY(  dl_dfdy(ii:ii+1,ij:ij+1)==dd_fill) .OR. &
269               &   ANY(dl_d2fdxy(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE
270
271               dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
272               ! check longitude discontinuity
273               IF( ld_discont )THEN
274
275                  dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
276                  dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
277                  IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
278                     WHERE( dl_tmp(:,:) < 0_dp )
279                        dl_tmp(:,:) = dl_tmp(:,:)+360._dp
280                     END WHERE
281                  ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
282                     WHERE( dl_tmp(:,:) > 180_dp )
283                        dl_tmp(:,:) = dl_tmp(:,:)-180._dp
284                     END WHERE
285                  ENDIF
286
287               ENDIF
288
289               ! compute bicubic coefficient
290               dl_coef(:)=interp_cubic__2D_coef(dl_tmp(:,:),&
291               &                                dl_dfdx(  ii:ii+1,ij:ij+1),&
292               &                                dl_dfdy(  ii:ii+1,ij:ij+1),&
293               &                                dl_d2fdxy(ii:ii+1,ij:ij+1),&
294               &                                dd_fill )
295
296               ! compute value on detetected point
297               CALL interp_cubic__2D_fill(dd_value( ji:ji+id_rhoi,   &
298               &                                    jj:jj+id_rhoj ), &
299               &                          id_detect(ji:ji+id_rhoi,   &
300               &                                    jj:jj+id_rhoj ), &
301               &                          dd_weight(:,:), dl_coef(:),&
302               &                          dd_fill, id_rhoi, id_rhoj )
303
304               IF( ld_discont )THEN
305                  WHERE( dd_value( ji:ji+id_rhoi, &
306                     &             jj:jj+id_rhoj ) >= 180._dp .AND. &
307                     &   dd_value( ji:ji+id_rhoi, &
308                     &             jj:jj+id_rhoj ) /= dd_fill )
309                     dd_value( ji:ji+id_rhoi, &
310                     &         jj:jj+id_rhoj ) = &
311                     &           dd_value( ji:ji+id_rhoi, &
312                     &                     jj:jj+id_rhoj ) - 360._dp
313                  END WHERE
314               ENDIF
315
316            ENDDO
317         ENDDO
318
319         DEALLOCATE(dl_coef)
320         DEALLOCATE(dl_tmp )
321
322         DEALLOCATE(dl_dfdx, &
323         &          dl_dfdy, &
324         &          dl_d2fdxy )
325
326         DEALLOCATE( dl_coarse )
327      ENDIF
328
329   END SUBROUTINE interp_cubic__2D
330   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
331   SUBROUTINE interp_cubic__1D(dd_value,  dd_fill,   &
332         &                     id_detect,            &
333         &                     dd_weight,            &
334         &                     id_rhoi,              &
335         &                     ld_discont)
336   !-------------------------------------------------------------------
337   !> @brief
338   !> This subroutine compute cubic interpolation on 1D array of value.
339   !>
340   !> @details
341   !>
342   !> @author J.Paul
343   !> @date September, 2014 - Initial Version
344   !>
345   !> @param[inout] dd_value  1D array of variable value
346   !> @param[in] dd_fill      FillValue of variable
347   !> @param[inout] id_detect 1D array of point to be interpolated
348   !> @param[in] id_rhoi      refinment factor
349   !> @param[in] ld_even      even refinment or not
350   !> @param[in] ld_discont   longitudinal discontinuity (-180°/180°, 0°/360°) or not
351   !-------------------------------------------------------------------
352
353      IMPLICIT NONE
354
355      ! Argument
356      REAL(dp)        , DIMENSION(:)  , INTENT(INOUT) :: dd_value
357      REAL(dp)                        , INTENT(IN   ) :: dd_fill
358      INTEGER(I4)     , DIMENSION(:)  , INTENT(INOUT) :: id_detect
359      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
360      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
361      LOGICAL                         , INTENT(IN   ) :: ld_discont
362
363      ! local variable
364      INTEGER(I4)                            :: il_xextra
365      INTEGER(i4), DIMENSION(1)              :: il_shape
366      INTEGER(i4), DIMENSION(1)              :: il_dim
367
368      REAL(dp)                               :: dl_min
369      REAL(dp)                               :: dl_max
370      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coef
371      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coarse
372      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_tmp
373      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_dfdx
374
375      ! loop indices
376      INTEGER(i4) :: ji
377      INTEGER(i4) :: ii
378      !----------------------------------------------------------------
379
380      IF( ANY(id_detect(:)==1) )THEN
381         il_shape(:)=SHAPE(dd_value)
382
383         ! compute coarse grid dimension
384         il_xextra=id_rhoi-1
385         il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
386
387         ALLOCATE( dl_coarse(il_dim(1)) )
388
389         ! value on coarse grid
390         dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )
391
392         ALLOCATE( dl_dfdx(il_dim(1)) )
393
394         ! compute derivative on coarse grid
395         dl_dfdx(:)=math_deriv_1D(dl_coarse(:), dd_fill, ld_discont)
396
397         ALLOCATE( dl_tmp(2) )
398         ALLOCATE( dl_coef(4) )
399
400         DO ji=1,il_shape(1)-1,id_rhoi
401            ii=((ji-1)/id_rhoi)+1
402
403            ! check if point to be interpolated
404            IF( ALL(id_detect(ji:ji+id_rhoi)==0) ) CYCLE
405            ! check data needed to interpolate
406            IF( ANY(dl_coarse(ii:ii+1)==dd_fill) .OR. &
407            &   ANY(  dl_dfdx(ii:ii+1)==dd_fill) ) CYCLE
408            ! check longitude discontinuity
409            dl_tmp(:)=dl_coarse(ii:ii+1)
410            IF( ld_discont )THEN
411
412               dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
413               dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
414               IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
415                  WHERE( dl_tmp(:) < 0_dp )
416                     dl_tmp(:) = dl_tmp(:)+360._dp
417                  END WHERE
418               ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
419                  WHERE( dl_tmp(:) > 180_dp )
420                     dl_tmp(:) = dl_tmp(:)-180._dp
421                  END WHERE
422               ENDIF
423
424            ENDIF
425
426            ! compute bicubic coefficient
427            dl_coef(:)=interp_cubic__1D_coef(dl_tmp(:),       &
428            &                                dl_dfdx(ii:ii+1),&
429            &                                dd_fill )
430
431            ! compute value on detetected point
432            CALL interp_cubic__1D_fill( dd_value( ji:ji+id_rhoi ), &
433            &                           id_detect(ji:ji+id_rhoi ), &
434            &                           dd_weight(:,:), dl_coef(:),  &
435            &                           dd_fill, id_rhoi )
436
437            IF( ld_discont )THEN
438               WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
439                  &   dd_value( ji:ji+id_rhoi ) /= dd_fill )
440                  dd_value(ji:ji+id_rhoi) = dd_value(ji:ji+id_rhoi) - 360._dp
441               END WHERE
442            ENDIF
443
444         ENDDO
445
446         DEALLOCATE(dl_coef)
447         DEALLOCATE(dl_tmp )
448
449         DEALLOCATE(dl_dfdx )
450         DEALLOCATE( dl_coarse )
451      ENDIF
452
453   END SUBROUTINE interp_cubic__1D
454   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
455   FUNCTION interp_cubic__2D_coef(dd_value,                    &
456         &                        dd_dfdx,  dd_dfdy, dd_d2fdxy,&
457         &                        dd_fill) &
458         & RESULT (df_coef)
459   !-------------------------------------------------------------------
460   !> @brief
461   !> This subroutine compute 2D array of coefficient for cubic interpolation.
462   !>
463   !> @author J.Paul
464   !> @date September, 2014 - Initial Version
465   !>
466   !> @param[in] dd_value  2D array of value
467   !> @param[in] dd_dfdx   2D array of first derivative in i-direction
468   !> @param[in] dd_dfdy   2D array of first derivative in j-direction
469   !> @param[in] dd_d2fdxy 2D array of cross derivative in i-j-direction
470   !> @param[in] dd_fill   FillValue of variable
471   !-------------------------------------------------------------------
472
473      IMPLICIT NONE
474
475      ! Argument
476      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_value
477      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdx
478      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_dfdy
479      REAL(dp), DIMENSION(:,:), INTENT(IN) :: dd_d2fdxy
480      REAL(dp)                , INTENT(IN) :: dd_fill
481
482      ! function
483      REAL(dp), DIMENSION(16)              :: df_coef
484
485      ! local variable
486      REAL(dp), DIMENSION(16,16), PARAMETER :: dl_matrix = RESHAPE( &
487      & (/ 1 , 0 ,-3 , 2 , 0 , 0 , 0 , 0 ,-3 , 0 , 9 ,-6 , 2 , 0 ,-6 , 4 ,&
488           0 , 0 , 3 ,-2 , 0 , 0 , 0 , 0 , 0 , 0 ,-9 , 6 , 0 , 0 , 6 ,-4 ,&
489           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 3 , 0 ,-9 , 6 ,-2 , 0 , 6 ,-4 ,&
490           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 9 ,-6 , 0 , 0 ,-6 , 4 ,&
491           0 , 1 ,-2 , 1 , 0 , 0 , 0 , 0 , 0 ,-3 , 6 ,-3 , 0 , 2 ,-4 , 2 ,&
492           0 , 0 ,-1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 3 ,-3 , 0 , 0 ,-2 , 2 ,&
493           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 3 ,-6 , 3 , 0 ,-2 , 4 ,-2 ,&
494           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 3 , 0 , 0 , 2 ,-2 ,&
495           0 , 0 , 0 , 0 , 1 , 0 ,-3 , 2 ,-2 , 0 , 6 ,-4 , 1 , 0 ,-3 , 2 ,&
496           0 , 0 , 0 , 0 , 0 , 0 , 3 ,-2 , 0 , 0 ,-6 , 4 , 0 , 0 , 3 ,-2 ,&
497           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-1 , 0 , 3 ,-2 , 1 , 0 ,-3 , 2 ,&
498           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-3 , 2 , 0 , 0 , 3 ,-2 ,&
499           0 , 0 , 0 , 0 , 0 , 1 ,-2 , 1 , 0 ,-2 , 4 ,-2 , 0 , 1 ,-2 , 1 ,&
500           0 , 0 , 0 , 0 , 0 , 0 ,-1 , 1 , 0 , 0 , 2 ,-2 , 0 , 0 ,-1 , 1 ,&
501           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ,-1 , 2 ,-1 , 0 , 1 ,-2 , 1 ,&
502           0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 ,-1 , 0 , 0 ,-1 , 1 /), &
503      & (/ 16, 16 /) )
504
505      REAL(dp), DIMENSION(16) :: dl_vect
506
507      !----------------------------------------------------------------
508      ! init
509      df_coef(:)=dd_fill
510
511      dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
512      dl_vect( 5: 8)=PACK(dd_dfdx(:,:),.TRUE. )
513      dl_vect( 9:12)=PACK(dd_dfdy(:,:),.TRUE. )
514      dl_vect(13:16)=PACK(dd_d2fdxy(:,:),.TRUE. )
515
516      df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
517
518   END FUNCTION interp_cubic__2D_coef
519   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
520   SUBROUTINE interp_cubic__2D_fill(dd_value, id_detect, &
521         &                          dd_weight, dd_coef,  &
522         &                          dd_fill, id_rhoi, id_rhoj)
523   !-------------------------------------------------------------------
524   !> @brief
525   !> This subroutine compute cubic interpolation of a 2D array of value.
526   !>
527   !> @author J.Paul
528   !> @date September, 2014 - Initial Version
529   !>
530   !> @param[inout] dd_value  2D array of mixed grid value
531   !> @param[inout] id_detect 2D array of point to be interpolated
532   !> @param[in] dd_coef      2D array of coefficient
533   !> @param[in] dd_fill      FillValue of variable
534   !> @param[in] ld_even      even refinment or not
535   !> @param[in] id_rhoi      refinement factor in i-direction
536   !> @param[in] id_rhoj      refinement factor in j-direction
537   !-------------------------------------------------------------------
538
539      IMPLICIT NONE
540
541      ! Argument
542      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value
543      INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
544      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
545      REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef
546      REAL(dp)                        , INTENT(IN   ) :: dd_fill
547      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi
548      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj
549
550      ! local variable
551
552      ! loop indices
553      INTEGER(i4) :: ii
554
555      INTEGER(i4) :: ji
556      INTEGER(i4) :: jj
557      !----------------------------------------------------------------
558
559      IF( ANY( dd_coef(:)==dd_fill ) )THEN
560         CALL logger_error("INTERP CUBIC FILL: fill value detected in coef . "//&
561         &              "can not compute interpolation.")
562      ELSE
563
564         ii=0
565         DO jj=1,id_rhoj+1
566            DO ji=1,id_rhoi+1
567
568               ii=ii+1
569               IF(id_detect(ji,jj)==1)THEN
570
571                  dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ii))
572                  id_detect(ji,jj)=0
573
574               ENDIF
575
576            ENDDO
577         ENDDO
578
579      ENDIF
580
581   END SUBROUTINE interp_cubic__2D_fill
582   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
583   FUNCTION interp_cubic__1D_coef(dd_value, &
584         &                        dd_dfdx,  &
585         &                        dd_fill)  &
586         &  RESULT (df_coef)
587   !-------------------------------------------------------------------
588   !> @brief
589   !> This subroutine compute 1D array of coefficient for cubic interpolation.
590   !>
591   !> @details
592   !>
593   !> @author J.Paul
594   !> @date September, 2014 - Initial Version
595   !>
596   !> @param[in] dd_value  1D array of value
597   !> @param[in] dd_dfdx   1D array of first derivative
598   !> @param[in] dd_fill   FillValue of variable
599   !-------------------------------------------------------------------
600
601      IMPLICIT NONE
602
603      ! Argument
604      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value
605      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_dfdx
606      REAL(dp)                , INTENT(IN) :: dd_fill
607
608      ! function
609      REAL(dp), DIMENSION(4)               :: df_coef
610
611      ! local variable
612      REAL(dp), DIMENSION(4,4), PARAMETER :: dl_matrix = RESHAPE( &
613      & (/  1 ,-1 ,-3 , 2 ,&
614            0 , 1 , 3 ,-2 ,&
615            0 , 0 ,-2 , 1 ,&
616            0 , 0 ,-1 , 1  /), &
617      & (/ 4, 4 /) )
618
619      REAL(dp), DIMENSION(4) :: dl_vect
620
621      !----------------------------------------------------------------
622      ! init
623      df_coef(:)=dd_fill
624
625      dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
626      dl_vect( 3: 4)=PACK(dd_dfdx(:),.TRUE. )
627
628      df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
629
630   END FUNCTION interp_cubic__1D_coef
631   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
632   SUBROUTINE interp_cubic__1D_fill(dd_value, id_detect, &
633         &                          dd_weight, dd_coef,  &
634         &                          dd_fill, id_rhoi)
635   !-------------------------------------------------------------------
636   !> @brief
637   !> This subroutine compute cubic interpolation of a 1D array of value.
638   !>
639   !> @author J.Paul
640   !> @date September, 2014 - Initial Version
641   !>
642   !> @param[inout] dd_value  1D array of mixed grid value
643   !> @param[inout] id_detect 1D array of point to be interpolated
644   !> @param[in] dd_coef      1D array of coefficient
645   !> @param[in] dd_fill      FillValue of variable
646   !> @param[in] ld_even      even refinment or not
647   !> @param[in] id_rho       refinement factor
648   !-------------------------------------------------------------------
649
650      IMPLICIT NONE
651
652      ! Argument
653      REAL(dp)        , DIMENSION(:)  , INTENT(INOUT) :: dd_value
654      INTEGER(i4)     , DIMENSION(:)  , INTENT(INOUT) :: id_detect
655      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
656      REAL(dp)        , DIMENSION(4)  , INTENT(IN   ) :: dd_coef
657      REAL(dp)                        , INTENT(IN   ) :: dd_fill
658      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
659
660      ! local variable
661      ! loop indices
662      INTEGER(i4) :: ji
663      !----------------------------------------------------------------
664
665      IF( ANY( dd_coef(:)==dd_fill ) )THEN
666         CALL logger_error("INTERP CUBIC FILL: fill value detected. "//&
667         &              "can not compute interpolation")
668      ELSE
669
670         DO ji=1,id_rhoi+1
671
672            IF(id_detect(ji)==1)THEN
673
674               dd_value(ji)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ji))
675               id_detect(ji)=0
676
677            ENDIF
678
679         ENDDO
680
681      ENDIF
682
683   END SUBROUTINE interp_cubic__1D_fill
684   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
685   SUBROUTINE interp_cubic__get_weight2D(dd_weight, &
686         &                               id_rho, ld_even)
687   !-------------------------------------------------------------------
688   !> @brief
689   !> This subroutine compute interpoaltion weight for 2D array.
690   !>
691   !> @author J.Paul
692   !> @date September, 2014 - Initial Version
693   !>
694   !> @param[in] dd_weight interpolation weight of 2D array
695   !> @param[in] ld_even   even refinment or not
696   !> @param[in] id_rho    refinement factor
697   !-------------------------------------------------------------------
698
699      IMPLICIT NONE
700
701      ! Argument
702      REAL(dp)   , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
703      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
704      LOGICAL    , DIMENSION(:)  , INTENT(IN   ) :: ld_even
705
706      ! local variable
707      REAL(dp)                  :: dl_dx
708      REAL(dp)                  :: dl_x
709      REAL(dp)                  :: dl_x2
710      REAL(dp)                  :: dl_x3
711      REAL(dp)                  :: dl_dy
712      REAL(dp)                  :: dl_y
713      REAL(dp)                  :: dl_y2
714      REAL(dp)                  :: dl_y3
715
716      ! loop indices
717      INTEGER(i4) :: ii
718      INTEGER(i4) :: ji
719      INTEGER(i4) :: jj
720      !----------------------------------------------------------------
721
722      IF( ld_even(jp_I) )THEN
723         dl_dx=1./REAL(id_rho(jp_I)-1)
724      ELSE ! odd refinement
725         dl_dx=1./REAL(id_rho(jp_I))
726      ENDIF
727
728      IF( ld_even(jp_J) )THEN
729         dl_dy=1./REAL(id_rho(jp_J)-1)
730      ELSE ! odd refinement
731         dl_dy=1./REAL(id_rho(jp_J))
732      ENDIF
733
734      ii=0
735      DO jj=1,id_rho(jp_J)+1
736
737         IF( ld_even(jp_J) )THEN
738            dl_y=(jj-1)*dl_dy - dl_dy*0.5
739         ELSE ! odd refinement
740            dl_y=(jj-1)*dl_dy
741         ENDIF
742         dl_y2=dl_y*dl_y
743         dl_y3=dl_y2*dl_y
744
745         DO ji=1,id_rho(jp_I)+1
746
747            ! iter
748            ii=ii+1
749
750            IF( ld_even(jp_I) )THEN
751               dl_x=(ji-1)*dl_dx - dl_dx*0.5
752            ELSE ! odd refinement
753               dl_x=(ji-1)*dl_dx
754            ENDIF
755            dl_x2=dl_x*dl_x
756            dl_x3=dl_x2*dl_x
757
758            dd_weight(:,ii)=(/1._dp, dl_x      , dl_x2      , dl_x3      , &
759            &                 dl_y , dl_x*dl_y , dl_x2*dl_y , dl_x3*dl_y , &
760            &                 dl_y2, dl_x*dl_y2, dl_x2*dl_y2, dl_x3*dl_y2, &
761            &                 dl_y3, dl_x*dl_y3, dl_x2*dl_y3, dl_x3*dl_y3 /)
762
763         ENDDO
764      ENDDO
765
766   END SUBROUTINE interp_cubic__get_weight2D
767   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
768   SUBROUTINE interp_cubic__get_weight1D(dd_weight, &
769         &                               id_rho, ld_even)
770   !-------------------------------------------------------------------
771   !> @brief
772   !> This subroutine compute interpoaltion weight for 1D array.
773   !>
774   !> @author J.Paul
775   !> @date September, 2014 - Initial Version
776   !>
777   !> @param[in] dd_weight interpolation weight of 1D array
778   !> @param[in] ld_even   even refinment or not
779   !> @param[in] id_rho    refinement factor
780   !-------------------------------------------------------------------
781
782      IMPLICIT NONE
783
784      ! Argument
785      REAL(dp)   , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
786      INTEGER(I4)                , INTENT(IN   ) :: id_rho
787      LOGICAL                    , INTENT(IN   ) :: ld_even
788
789      ! local variable
790      REAL(dp)                  :: dl_dx
791      REAL(dp)                  :: dl_x
792      REAL(dp)                  :: dl_x2
793      REAL(dp)                  :: dl_x3
794
795      ! loop indices
796      INTEGER(i4) :: ji
797      !----------------------------------------------------------------
798
799      IF( ld_even )THEN
800         dl_dx=1./REAL(id_rho-1)
801      ELSE ! odd refinement
802         dl_dx=1./REAL(id_rho)
803      ENDIF
804
805      DO ji=1,id_rho+1
806         IF( ld_even )THEN
807            dl_x=(ji-1)*dl_dx - dl_dx*0.5
808         ELSE ! odd refinement
809            dl_x=(ji-1)*dl_dx
810         ENDIF
811         dl_x2=dl_x*dl_x
812         dl_x3=dl_x2*dl_x
813
814         dd_weight(:,ji)=(/1._dp, dl_x, dl_x2, dl_x3 /)
815      ENDDO
816
817   END SUBROUTINE interp_cubic__get_weight1D
818   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
819END MODULE interp_cubic
Note: See TracBrowser for help on using the repository browser.