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 branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/TOOLS/SIREN/src/interp_cubic.f90 @ 5581

Last change on this file since 5581 was 5581, checked in by timgraham, 9 years ago

Merged head of trunk into branch

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