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_linear.f90 in branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/TOOLS/SIREN/src/interp_linear.f90 @ 5942

Last change on this file since 5942 was 5942, checked in by rfurner, 8 years ago

merged bug fixes from vn3.6_stable to this branch

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