source: utils/tools/SIREN/src/interp_linear.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 10 months ago

update nemo trunk

File size: 26.6 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! DESCRIPTION:
6!> @brief
7!> This module manage linear interpolation on regular grid.
8!>
9!> @details
10!> to compute linear interpolation:<br/>
11!> @code
12!> CALL interp_linear_fill(dd_value, dd_fill, id_detect, id_rho, ld_even [,ld_discont] )
13!> @endcode
14!>    - dd_value is 2D array of variable value
15!>    - dd_fill is the FillValue of variable
16!>    - id_detect is 2D array of point to be interpolated (see interp module)
17!>    - id_rho  is array of refinment factor
18!>    - ld_even indicates even refinment or not
19!>    - ld_discont indicates longitudinal discontinuity (-180°/180°, 0°/360°) or not
20!>
21!> @author
22!> J.Paul
23!>
24!> @date September, 2014 - Initial version
25!>
26!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
27!----------------------------------------------------------------------
28MODULE interp_linear
29
30   USE netcdf                          ! nf90 library
31   USE global                          ! global variable
32   USE kind                            ! F90 kind parameter
33   USE logger                          ! log file manager
34   USE fct                             ! basic useful function
35   USE extrap                          ! extrapolation manager
36
37   IMPLICIT NONE
38   ! NOTE_avoid_public_variables_if_possible
39
40   ! type and variable
41
42   ! function and subroutine
43   PUBLIC :: interp_linear_fill  !< compute interpolation using linear method
44
45   PRIVATE :: interp_linear__2D           !< compute bilinear interpolation on 2D gid
46   PRIVATE :: interp_linear__1D           !< compute   linear interpolation on 1D gid
47   PRIVATE :: interp_linear__2D_coef      !< compute coefficient for bilinear interpolation
48   PRIVATE :: interp_linear__2D_fill      !< fill value using bilinear interpolation
49   PRIVATE :: interp_linear__1D_coef      !< compute coefficient for   linear interpolation
50   PRIVATE :: interp_linear__1D_fill      !< fill value using   linear interpolation
51   PRIVATE :: interp_linear__get_weight2D !< compute interpoaltion weight for 2D array.
52   PRIVATE :: interp_linear__get_weight1D !< compute interpoaltion weight for 1D array.
53
54CONTAINS
55   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
56   SUBROUTINE interp_linear_fill(dd_value, dd_fill, id_detect, &
57         &                       id_rho, ld_even, ld_discont)
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
76      IMPLICIT NONE
77
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   SUBROUTINE interp_linear__2D(dd_value,  dd_fill,   &
167         &                      id_detect,            &
168         &                      dd_weight,            &
169         &                      id_rhoi, id_rhoj,     &
170         &                      ld_discont)
171   !-------------------------------------------------------------------
172   !> @brief
173   !> This subroutine compute linear interpolation on 2D array of value.
174   !>
175   !> @details
176   !>
177   !> @author J.Paul
178   !> @date September, 2014 - Initial Version
179   !>
180   !> @param[inout] dd_value  2D array of variable value
181   !> @param[in] dd_fill      FillValue of variable
182   !> @param[inout] id_detect 2D array of point to be interpolated
183   !> @param[in] id_rhoi      refinment factor in i-direction
184   !> @param[in] id_rhoj      refinment factor in j-direction
185   !> @param[in] id_rhok      refinment factor in k-direction
186   !> @param[in] ld_even      even refinment or not
187   !> @param[in] ld_discont   longitudinal discontinuity (-180°/180°, 0°/360°) or not
188   !-------------------------------------------------------------------
189
190      IMPLICIT NONE
191
192      ! Argument
193      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
194      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
195      INTEGER(I4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
196      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
197      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
198      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoj
199      LOGICAL                         , INTENT(IN   ) :: ld_discont
200
201      ! local variable
202      INTEGER(I4)                              :: il_xextra
203      INTEGER(I4)                              :: il_yextra
204      INTEGER(i4), DIMENSION(2)                :: il_shape
205      INTEGER(i4), DIMENSION(2)                :: il_dim
206
207      REAL(dp)                                 :: dl_min 
208      REAL(dp)                                 :: dl_max 
209      REAL(dp)   , DIMENSION(:)  , ALLOCATABLE :: dl_coef
210      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_coarse
211      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_tmp
212
213      ! loop indices
214      INTEGER(i4) :: ji
215      INTEGER(i4) :: jj
216      INTEGER(i4) :: ii
217      INTEGER(i4) :: ij
218
219      !----------------------------------------------------------------
220
221      IF( ANY(id_detect(:,:)==1) )THEN
222         il_shape(:)=SHAPE(dd_value)
223
224         ! compute coarse grid dimension
225         il_xextra=id_rhoi-1
226         il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
227
228         il_yextra=id_rhoj-1
229         il_dim(2)=(il_shape(2)+il_yextra)/id_rhoj
230
231         ALLOCATE( dl_coarse(il_dim(1),il_dim(2)) )
232
233         ! value on coarse grid
234         dl_coarse(:,:)=dd_value( 1:il_shape(1):id_rhoi, &
235         &                        1:il_shape(2):id_rhoj )
236
237         ALLOCATE( dl_tmp(2,2) )
238         ALLOCATE( dl_coef(4) )
239
240         DO jj=1,il_shape(2)-1,id_rhoj
241            ij=((jj-1)/id_rhoj)+1
242            DO ji=1,il_shape(1)-1,id_rhoi
243               ii=((ji-1)/id_rhoi)+1
244         
245               ! check if point to be interpolated
246               IF( ALL(id_detect(ji:ji+id_rhoi,   &
247               &                 jj:jj+id_rhoj)==0) ) CYCLE
248               ! check data needed to interpolate
249               IF( ANY(dl_coarse(ii:ii+1,ij:ij+1)==dd_fill) ) CYCLE
250               ! check longitude discontinuity
251               dl_tmp(:,:)=dl_coarse(ii:ii+1,ij:ij+1)
252               IF( ld_discont )THEN
253
254                  dl_min=MINVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
255                  dl_max=MAXVAL( dl_tmp(:,:), dl_tmp(:,:)/=dd_fill )
256                  IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
257                     WHERE( dl_tmp(:,:) < 0_dp ) 
258                        dl_tmp(:,:) = dl_tmp(:,:)+360._dp
259                     END WHERE
260                  ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
261                     WHERE( dl_tmp(:,:) > 180_dp ) 
262                        dl_tmp(:,:) = dl_tmp(:,:)-180._dp
263                     END WHERE
264                  ENDIF
265
266               ENDIF
267
268               ! compute bilinear coefficient
269               dl_coef(:)=interp_linear__2D_coef(dl_tmp(:,:),&
270               &                                 dd_fill )
271
272               ! compute value on detetected point
273               CALL interp_linear__2D_fill(dd_value( ji:ji+id_rhoi,   &
274               &                                     jj:jj+id_rhoj ), &
275               &                           id_detect(ji:ji+id_rhoi,   &
276               &                                     jj:jj+id_rhoj ), &
277               &                           dd_weight(:,:), dl_coef(:),&
278               &                           dd_fill, id_rhoi, id_rhoj )
279
280               IF( ld_discont )THEN
281                  WHERE( dd_value( ji:ji+id_rhoi, &
282                     &             jj:jj+id_rhoj ) >= 180._dp .AND. &
283                     &   dd_value( ji:ji+id_rhoi, &
284                     &             jj:jj+id_rhoj ) /= dd_fill )
285                     dd_value( ji:ji+id_rhoi, &
286                     &         jj:jj+id_rhoj ) = &
287                     &           dd_value( ji:ji+id_rhoi, &
288                     &                     jj:jj+id_rhoj ) - 360._dp
289                  END WHERE
290               ENDIF
291
292            ENDDO
293         ENDDO
294
295         DEALLOCATE(dl_coef)
296         DEALLOCATE(dl_tmp )
297
298         DEALLOCATE( dl_coarse )
299      ENDIF
300
301   END SUBROUTINE interp_linear__2D
302   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
303   SUBROUTINE interp_linear__1D(dd_value,  dd_fill,   &
304         &                      id_detect,            &
305         &                      dd_weight,            &
306         &                      id_rhoi,              &
307         &                      ld_discont)
308   !-------------------------------------------------------------------
309   !> @brief
310   !> This subroutine compute linear interpolation on 1D array of value.
311   !>
312   !> @details
313   !>
314   !> @author J.Paul
315   !> @date September, 2014 - Initial Version
316   !>
317   !> @param[inout] dd_value  1D array of variable value
318   !> @param[in] dd_fill      FillValue of variable
319   !> @param[inout] id_detect 1D array of point to be interpolated
320   !> @param[in] id_rhoi      refinment factor
321   !> @param[in] ld_even      even refinment or not
322   !> @param[in] ld_discont   longitudinal discontinuity (-180°/180°, 0°/360°) or not
323   !-------------------------------------------------------------------
324
325      IMPLICIT NONE
326
327      ! Argument
328      REAL(dp)        , DIMENSION(:)  , INTENT(INOUT) :: dd_value 
329      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
330      INTEGER(I4)     , DIMENSION(:)  , INTENT(INOUT) :: id_detect
331      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
332      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
333      LOGICAL                         , INTENT(IN   ) :: ld_discont
334
335      ! local variable
336      INTEGER(I4)                            :: il_xextra
337      INTEGER(i4), DIMENSION(1)              :: il_shape
338      INTEGER(i4), DIMENSION(1)              :: il_dim
339
340      REAL(dp)                               :: dl_min 
341      REAL(dp)                               :: dl_max 
342      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coef
343      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_coarse
344      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_tmp
345
346      ! loop indices
347      INTEGER(i4) :: ji
348      INTEGER(i4) :: ii
349
350      !----------------------------------------------------------------
351
352      IF( ANY(id_detect(:)==1) )THEN
353         il_shape(:)=SHAPE(dd_value)
354
355         ! compute coarse grid dimension
356         il_xextra=id_rhoi-1
357         il_dim(1)=(il_shape(1)+il_xextra)/id_rhoi
358
359         ALLOCATE( dl_coarse(il_dim(1)) )
360
361         ! value on coarse grid
362         dl_coarse(:)=dd_value( 1:il_shape(1):id_rhoi )
363
364         ALLOCATE( dl_tmp(2) )
365         ALLOCATE( dl_coef(4) )
366
367         DO ji=1,il_shape(1)-1,id_rhoi
368            ii=((ji-1)/id_rhoi)+1
369         
370            ! check if point to be interpolated
371            IF( ALL(id_detect(ji:ji+id_rhoi)==0) ) CYCLE
372            ! check data needed to interpolate
373            IF( ANY(dl_coarse(ii:ii+1)==dd_fill) ) CYCLE
374            ! check longitude discontinuity
375            dl_tmp(:)=dl_coarse(ii:ii+1)
376            IF( ld_discont )THEN
377
378               dl_min=MINVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
379               dl_max=MAXVAL( dl_tmp(:), dl_tmp(:)/=dd_fill )
380               IF( dl_min < -170_dp .AND. dl_max > 170_dp )THEN
381                  WHERE( dl_tmp(:) < 0_dp ) 
382                     dl_tmp(:) = dl_tmp(:)+360._dp
383                  END WHERE
384               ELSEIF( dl_min < 10_dp .AND. dl_max > 350_dp )THEN
385                  WHERE( dl_tmp(:) > 180_dp ) 
386                     dl_tmp(:) = dl_tmp(:)-180._dp
387                  END WHERE
388               ENDIF
389
390            ENDIF
391
392            ! compute bilinear coefficient
393            dl_coef(:)=interp_linear__1D_coef(dl_tmp(:),       &
394            &                                 dd_fill )
395
396            ! compute value on detetected point
397            CALL interp_linear__1D_fill( dd_value( ji:ji+id_rhoi ), &
398            &                            id_detect(ji:ji+id_rhoi ), &
399            &                            dd_weight(:,:), dl_coef(:),&
400            &                            dd_fill, id_rhoi )
401
402            IF( ld_discont )THEN
403               WHERE( dd_value( ji:ji+id_rhoi ) >= 180._dp .AND. &
404                  &   dd_value( ji:ji+id_rhoi ) /= dd_fill )
405                  dd_value(ji:ji+id_rhoi) = dd_value(ji:ji+id_rhoi) - 360._dp
406               END WHERE
407            ENDIF
408
409         ENDDO
410
411         DEALLOCATE(dl_coef)
412         DEALLOCATE(dl_tmp )
413
414         DEALLOCATE( dl_coarse )
415      ENDIF
416
417   END SUBROUTINE interp_linear__1D
418   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
419   FUNCTION interp_linear__2D_coef(dd_value, dd_fill) &
420         & RESULT (df_coef)
421   !-------------------------------------------------------------------
422   !> @brief
423   !> This subroutine compute 2D array of coefficient for linear interpolation.
424   !>
425   !> @author J.Paul
426   !> @date September, 2014 - Initial Version
427   !>
428   !> @param[in] dd_value  2D array of value
429   !> @param[in] dd_fill   FillValue of variable
430   !-------------------------------------------------------------------
431
432      IMPLICIT NONE
433
434      ! Argument
435      REAL(dp), DIMENSION(:,:)  , INTENT(IN) :: dd_value 
436      REAL(dp)                  , INTENT(IN) :: dd_fill 
437
438      ! function
439      REAL(dp), DIMENSION(4)                 :: df_coef
440
441      ! local variable
442      REAL(dp), DIMENSION(4,4), PARAMETER :: dl_matrix = RESHAPE( &
443      & (/  1 ,-1 ,-1 , 1 ,& 
444            0 , 1 , 0 ,-1 ,& 
445            0 , 0 , 1 ,-1 ,& 
446            0 , 0 , 0 , 1 /), &
447      & (/ 4, 4 /) )
448
449      REAL(dp), DIMENSION(4) :: dl_vect
450
451      !----------------------------------------------------------------
452      ! init
453      df_coef(:)=dd_fill
454
455      dl_vect( 1: 4)=PACK(dd_value(:,:),.TRUE. )
456      df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
457
458   END FUNCTION interp_linear__2D_coef
459   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460   SUBROUTINE interp_linear__2D_fill(dd_value, id_detect, &
461         &                           dd_weight, dd_coef,  &
462         &                           dd_fill, id_rhoi, id_rhoj)
463   !-------------------------------------------------------------------
464   !> @brief
465   !> This subroutine compute linear interpolation of a 2D array of value.
466   !>
467   !> @author J.Paul
468   !> @date September, 2014 - Initial Version
469   !>
470   !> @param[inout] dd_value  2D array of mixed grid value
471   !> @param[inout] id_detect 2D array of point to be interpolated
472   !> @param[in] dd_coef      2D array of coefficient
473   !> @param[in] dd_fill      FillValue of variable
474   !> @param[in] ld_even      even refinment or not
475   !> @param[in] id_rhoi      refinement factor in i-direction
476   !> @param[in] id_rhoj      refinement factor in j-direction
477   !-------------------------------------------------------------------
478
479      IMPLICIT NONE
480
481      ! Argument
482      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
483      INTEGER(i4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
484      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
485      REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef 
486      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
487      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoi
488      INTEGER(I4)     ,                 INTENT(IN   ) :: id_rhoj
489
490      ! local variable
491
492      ! loop indices
493      INTEGER(i4) :: ii
494
495      INTEGER(i4) :: ji
496      INTEGER(i4) :: jj
497      !----------------------------------------------------------------
498
499      IF( ANY( dd_coef(:)==dd_fill ) )THEN
500         CALL logger_error("INTERP LINEAR FILL: fill value detected in coef. "//&
501         &              "can not compute interpolation.")
502      ELSE
503
504         ii=0
505         DO jj=1,id_rhoj+1
506            DO ji=1,id_rhoi+1
507
508               ii=ii+1
509               IF(id_detect(ji,jj)==1)THEN
510
511                  dd_value(ji,jj)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ii))
512                  id_detect(ji,jj)=0
513
514               ENDIF
515
516            ENDDO
517         ENDDO
518
519      ENDIF
520
521   END SUBROUTINE interp_linear__2D_fill
522   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
523   FUNCTION interp_linear__1D_coef(dd_value, dd_fill) &
524         & RESULT (df_coef)
525   !-------------------------------------------------------------------
526   !> @brief
527   !> This subroutine compute 1D array of coefficient for linear interpolation.
528   !>
529   !> @author J.Paul
530   !> @date September, 2014 - Initial Version
531   !>
532   !> @param[in] dd_value  1D array of value
533   !> @param[in] dd_fill   FillValue of variable
534   !-------------------------------------------------------------------
535
536      IMPLICIT NONE
537
538      ! Argument
539      REAL(dp), DIMENSION(:)  , INTENT(IN) :: dd_value 
540      REAL(dp)                , INTENT(IN) :: dd_fill 
541
542      ! function
543      REAL(dp), DIMENSION(2)               :: df_coef
544
545      ! local variable
546      REAL(dp), DIMENSION(2,2), PARAMETER :: dl_matrix = RESHAPE( &
547      & (/  1 ,-1 ,& 
548            0 , 1  /), &
549      &  (/ 2, 2 /) )
550
551      REAL(dp), DIMENSION(2) :: dl_vect
552
553      !----------------------------------------------------------------
554      ! init
555      df_coef(:)=dd_fill
556
557      dl_vect( 1: 2)=PACK(dd_value(:),.TRUE. )
558      df_coef(:)=MATMUL(dl_matrix(:,:),dl_vect(:))
559
560   END FUNCTION interp_linear__1D_coef
561   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
562   SUBROUTINE interp_linear__1D_fill(dd_value, id_detect, &
563         &                           dd_weight, dd_coef,  &
564         &                           dd_fill, id_rho)
565   !-------------------------------------------------------------------
566   !> @brief
567   !> This subroutine compute linear interpolation of a 1D array of value.
568   !>
569   !> @author J.Paul
570   !> @date September, 2014 - Initial Version
571   !>
572   !> @param[inout] dd_value  1D array of mixed grid value
573   !> @param[inout] id_detect 1D array of point to be interpolated
574   !> @param[in] dd_coef      1D array of coefficient
575   !> @param[in] dd_fill      FillValue of variable
576   !> @param[in] ld_even      even refinment or not
577   !> @param[in] id_rho       refinement factor
578   !-------------------------------------------------------------------
579
580      IMPLICIT NONE
581
582      ! Argument
583      REAL(dp)        , DIMENSION(:)  , INTENT(INOUT) :: dd_value 
584      INTEGER(i4)     , DIMENSION(:)  , INTENT(INOUT) :: id_detect
585      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_weight
586      REAL(dp)        , DIMENSION(:)  , INTENT(IN   ) :: dd_coef 
587      REAL(dp)                        , INTENT(IN   ) :: dd_fill 
588      INTEGER(I4)                     , INTENT(IN   ) :: id_rho     
589
590      ! local variable
591
592      ! loop indices
593      INTEGER(i4) :: ji
594      !----------------------------------------------------------------
595
596      IF( ANY( dd_coef(:)==dd_fill ) )THEN
597         CALL logger_error("INTERP LINEAR FILL: fill value detected. "//&
598         &              "can not compute interpolation")
599      ELSE
600
601         DO ji=1,id_rho+1
602           
603            IF(id_detect(ji)==1)THEN
604
605               dd_value(ji)=DOT_PRODUCT(dd_coef(:),dd_weight(:,ji))
606               id_detect(ji)=0
607
608            ENDIF
609
610         ENDDO
611
612      ENDIF
613
614   END SUBROUTINE interp_linear__1D_fill
615   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616   SUBROUTINE interp_linear__get_weight2D(dd_weight, id_rho, ld_even)
617   !-------------------------------------------------------------------
618   !> @brief
619   !> This subroutine compute interpoaltion weight for 2D array.
620   !>
621   !> @author J.Paul
622   !> @date September, 2014 - Initial Version
623   !>
624   !> @param[in] dd_weight interpolation weight of 2D array
625   !> @param[in] ld_even   even refinment or not
626   !> @param[in] id_rho    refinement factor
627   !-------------------------------------------------------------------
628
629      IMPLICIT NONE
630
631      ! Argument
632
633      REAL(dp)   , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
634      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
635      LOGICAL    , DIMENSION(:)  , INTENT(IN   ) :: ld_even
636
637      ! local variable
638      REAL(dp)                  :: dl_dx
639      REAL(dp)                  :: dl_x     
640      REAL(dp)                  :: dl_dy
641      REAL(dp)                  :: dl_y     
642     
643      ! loop indices
644      INTEGER(i4) :: ii 
645
646      INTEGER(i4) :: ji 
647      INTEGER(i4) :: jj
648      !----------------------------------------------------------------
649     
650      IF( ld_even(jp_I) )THEN
651         dl_dx=1._dp/REAL(id_rho(jp_I)-1,dp)
652      ELSE ! odd refinement
653         dl_dx=1._dp/REAL(id_rho(jp_I),dp)
654      ENDIF
655
656      IF( ld_even(jp_J) )THEN
657         dl_dy=1._dp/REAL(id_rho(jp_J)-1,dp)
658      ELSE ! odd refinement
659         dl_dy=1._dp/REAL(id_rho(jp_J),dp)
660      ENDIF
661
662      ii=0
663      DO jj=1,id_rho(jp_J)+1
664
665         IF( ld_even(jp_J) )THEN
666            dl_y=REAL(jj-1,dp)*dl_dy - dl_dy*0.5_dp
667         ELSE ! odd refinement
668            dl_y=REAL(jj-1,dp)*dl_dy 
669         ENDIF
670
671         DO ji=1,id_rho(jp_I)+1
672
673            ! iter
674            ii=ii+1
675
676            IF( ld_even(jp_I) )THEN
677               dl_x=REAL(ji-1,dp)*dl_dx - dl_dx*0.5_dp 
678            ELSE ! odd refinement
679               dl_x=REAL(ji-1,dp)*dl_dx 
680            ENDIF
681
682            dd_weight(:,ii)=(/1._dp, dl_x, dl_y, dl_x*dl_y /)
683
684         ENDDO
685      ENDDO
686
687   END SUBROUTINE interp_linear__get_weight2D
688   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689   SUBROUTINE interp_linear__get_weight1D(dd_weight, id_rho, ld_even)
690   !-------------------------------------------------------------------
691   !> @brief
692   !> This subroutine compute interpoaltion weight for 1D array.
693   !>
694   !> @author J.Paul
695   !> @date September, 2014 - Initial Version
696   !>
697   !> @param[in] dd_weight interpolation weight of 1D array
698   !> @param[in] ld_even   even refinment or not
699   !> @param[in] id_rho    refinement factor
700   !-------------------------------------------------------------------
701
702      IMPLICIT NONE
703
704      ! Argument
705      REAL(dp)   , DIMENSION(:,:), INTENT(INOUT) :: dd_weight
706      INTEGER(I4)                , INTENT(IN   ) :: id_rho
707      LOGICAL                    , INTENT(IN   ) :: ld_even
708
709      ! local variable
710      REAL(dp)                  :: dl_dx
711      REAL(dp)                  :: dl_x     
712
713      ! loop indices
714      INTEGER(i4) :: ji   
715      !----------------------------------------------------------------
716     
717      IF( ld_even )THEN
718         dl_dx=1._dp/REAL(id_rho-1,dp)
719      ELSE ! odd refinement
720         dl_dx=1._dp/REAL(id_rho,dp)
721      ENDIF
722
723      DO ji=1,id_rho+1
724         IF( ld_even )THEN
725            dl_x=REAL(ji-1,dp)*dl_dx - dl_dx*0.5_dp 
726         ELSE ! odd refinement
727            dl_x=REAL(ji-1,dp)*dl_dx 
728         ENDIF
729
730         dd_weight(:,ji)=(/1._dp, dl_x /)
731      ENDDO
732
733   END SUBROUTINE interp_linear__get_weight1D
734   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
735END MODULE interp_linear
Note: See TracBrowser for help on using the repository browser.