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 NEMO/trunk/tools/SIREN/src – NEMO

source: NEMO/trunk/tools/SIREN/src/interp_cubic.f90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

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