source: trunk/NEMOGCM/TOOLS/SIREN/src/interp_linear.f90 @ 5037

Last change on this file since 5037 was 5037, checked in by jpaul, 6 years ago

see ticket #1439

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