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

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

see ticket #1439

File size: 11.2 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 nearest interpolation on regular grid.
10!>
11!> @details
12!> to compute nearest interpolation:<br/>
13!> @code
14!> CALL interp_nearest_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!----------------------------------------------------------------------
30MODULE interp_nearest
31
32   USE netcdf                          ! nf90 library
33   USE global                          ! global variable
34   USE kind                            ! F90 kind parameter
35   USE logger                          ! log file manager
36   USE fct                             ! basic useful function
37
38   IMPLICIT NONE
39   ! NOTE_avoid_public_variables_if_possible
40
41   ! type and variable
42
43   ! function and subroutine
44   PUBLIC :: interp_nearest_fill  !< compute interpolation using nearest method
45
46   PRIVATE :: interp_nearest__2D       !< compute binearest interpolation on 2D gid
47   PRIVATE :: interp_nearest__1D       !< compute   nearest interpolation on 1D gid
48   PRIVATE :: interp_nearest__2D_fill  !< fill value using binearest interpolation
49   PRIVATE :: interp_nearest__1D_fill  !< fill value using   nearest interpolation
50
51CONTAINS   
52   !-------------------------------------------------------------------
53   !> @brief
54   !> This subroutine compute horizontal nearest interpolation on 4D array of value.
55   !>
56   !> @author J.Paul
57   !> - September, 2014- Initial Version
58   !>
59   !> @param[inout] dd_value  2D array of variable value
60   !> @param[inout] id_detect 2D array of point to be interpolated
61   !> @param[in]    id_rho    array of refinment factor
62   !-------------------------------------------------------------------
63   SUBROUTINE interp_nearest_fill(dd_value, id_detect, id_rho )
64      IMPLICIT NONE
65      ! Argument
66      REAL(dp)        , DIMENSION(:,:,:,:), INTENT(INOUT) :: dd_value 
67      INTEGER(I4)     , DIMENSION(:,:,:)  , INTENT(INOUT) :: id_detect
68      INTEGER(I4)     , DIMENSION(:)      , INTENT(IN   ) :: id_rho
69
70      ! local variable
71      INTEGER(i4), DIMENSION(4)                :: il_shape
72
73      ! loop indices
74      INTEGER(i4) :: ji
75      INTEGER(i4) :: jj
76      INTEGER(i4) :: jk
77      INTEGER(i4) :: jl
78      !----------------------------------------------------------------
79
80      il_shape(:)=SHAPE(dd_value)
81
82      DO jl=1,il_shape(4)
83         ! loop on vertical level
84         DO jk=1,il_shape(3)
85
86            ! I-J plan
87            CALL interp_nearest__2D(dd_value(:,:,jk,jl),&
88            &                       id_detect(:,:,jk),  &
89            &                       id_rho(jp_I), id_rho(jp_J) )           
90            IF( ANY(id_detect(:,:,jk)==1) )THEN
91               ! I direction
92               DO jj=1,il_shape(2)
93                  CALL interp_nearest__1D( dd_value(:,jj,jk,jl),&
94                  &                        id_detect(:,jj,jk),  &
95                  &                        id_rho(jp_I) )
96               ENDDO
97               IF( ALL(id_detect(:,:,jk)==0) )THEN
98                  CYCLE
99               ELSE
100                  ! J direction
101                  DO ji=1,il_shape(1)
102                     CALL interp_nearest__1D( dd_value(ji,:,jk,jl),&
103                     &                        id_detect(ji,:,jk),  &
104                     &                        id_rho(jp_J) )
105                  ENDDO
106               ENDIF
107            ENDIF
108
109         ENDDO
110      ENDDO
111
112   END SUBROUTINE interp_nearest_fill
113   !-------------------------------------------------------------------
114   !> @brief
115   !> This subroutine compute nearest interpolation on 2D array of value.
116   !>
117   !> @author J.Paul
118   !> - September, 2014- Initial Version
119   !>
120   !> @param[inout] dd_value  2D array of variable value
121   !> @param[inout] id_detect 2D array of point to be interpolated
122   !> @param[in] id_rhoi      refinment factor in i-direction
123   !> @param[in] id_rhoj      refinment factor in j-direction
124   !> @param[in] id_rhok      refinment factor in k-direction
125   !-------------------------------------------------------------------
126   SUBROUTINE interp_nearest__2D( dd_value, id_detect, &
127      &                           id_rhoi, id_rhoj )
128
129      IMPLICIT NONE
130      ! Argument
131      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value 
132      INTEGER(I4)     , DIMENSION(:,:), INTENT(INOUT) :: id_detect
133      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoi
134      INTEGER(I4)                     , INTENT(IN   ) :: id_rhoj
135
136      ! local variable
137      INTEGER(i4), DIMENSION(2)                :: il_shape
138
139      ! loop indices
140      INTEGER(i4) :: ji
141      INTEGER(i4) :: jj
142
143      !----------------------------------------------------------------
144
145      IF( ANY(id_detect(:,:)==1) )THEN
146
147         il_shape(:)=SHAPE(dd_value)
148
149         DO jj=1,il_shape(2)-1,id_rhoj
150            DO ji=1,il_shape(1)-1,id_rhoi
151         
152               ! check if point to be interpolated
153               IF( ALL(id_detect(ji:ji+id_rhoi,   &
154               &                 jj:jj+id_rhoj)==0) ) CYCLE
155
156               ! compute value on detetected point
157               CALL interp_nearest__2D_fill(dd_value( ji:ji+id_rhoi,   &
158               &                                      jj:jj+id_rhoj ), &
159               &                            id_detect(ji:ji+id_rhoi,   &
160               &                                      jj:jj+id_rhoj ) )
161
162            ENDDO
163         ENDDO
164
165      ENDIF
166
167   END SUBROUTINE interp_nearest__2D
168   !-------------------------------------------------------------------
169   !> @brief
170   !> This subroutine compute nearest interpolation on 1D array of value.
171   !>
172   !> @author J.Paul
173   !> - September, 2014- Initial Version
174   !>
175   !> @param[inout] dd_value  1D array of variable value
176   !> @param[inout] id_detect 1D array of point to be interpolated
177   !> @param[in]    id_rhoi   refinment factor
178   !-------------------------------------------------------------------
179   SUBROUTINE interp_nearest__1D( dd_value,  id_detect, &
180      &                           id_rhoi )
181
182      IMPLICIT NONE
183      ! Argument
184      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value 
185      INTEGER(I4)     , DIMENSION(:), INTENT(INOUT) :: id_detect
186      INTEGER(I4)                   , INTENT(IN   ) :: id_rhoi
187
188      ! local variable
189      INTEGER(i4), DIMENSION(1)              :: il_shape
190
191      ! loop indices
192      INTEGER(i4) :: ji
193
194      !----------------------------------------------------------------
195
196      IF( ANY(id_detect(:)==1) )THEN
197         il_shape(:)=SHAPE(dd_value)
198
199         DO ji=1,il_shape(1)-1,id_rhoi
200         
201            ! check if point to be interpolated
202            IF( ALL(id_detect(ji:ji+id_rhoi)==0) ) CYCLE
203
204            ! compute value on detetected point
205            CALL interp_nearest__1D_fill( dd_value( ji:ji+id_rhoi ), &
206            &                             id_detect(ji:ji+id_rhoi ) )
207
208         ENDDO
209
210      ENDIF
211
212   END SUBROUTINE interp_nearest__1D
213   !-------------------------------------------------------------------
214   !> @brief
215   !> This subroutine compute nearest interpolation of a 2D array of value.
216   !>
217   !> @author J.Paul
218   !> - September, 2014- Initial Version
219   !>
220   !> @param[inout] dd_value  2D array of mixed grid value
221   !> @param[inout] id_detect 2D array of point to be interpolated
222   !-------------------------------------------------------------------
223   SUBROUTINE interp_nearest__2D_fill( dd_value, id_detect )
224      IMPLICIT NONE
225      ! Argument
226      REAL(dp)   , DIMENSION(:,:)  , INTENT(INOUT) :: dd_value 
227      INTEGER(i4), DIMENSION(:,:)  , INTENT(INOUT) :: id_detect
228
229      ! local variable
230      INTEGER(i4), DIMENSION(2) :: il_shape
231
232      INTEGER(i4) :: il_i1
233      INTEGER(i4) :: il_i2
234      INTEGER(i4) :: il_j1
235      INTEGER(i4) :: il_j2
236
237      INTEGER(i4) :: il_half1
238      INTEGER(i4) :: il_half2
239
240      ! loop indices
241      INTEGER(i4) :: ji
242      INTEGER(i4) :: jj
243      !----------------------------------------------------------------
244
245      il_shape(:)=SHAPE(dd_value(:,:))
246
247      il_i1=1
248      il_i2=il_shape(1)
249
250      il_j1=1
251      il_j2=il_shape(2)
252
253      il_half1=CEILING(il_shape(1)*0.5)
254      il_half2=CEILING(il_shape(2)*0.5)
255
256      DO jj=1,il_half2
257
258         DO ji=1,il_half1
259           
260            ! lower left point
261            IF(id_detect(ji,jj)==1)THEN
262
263               dd_value( ji,jj)=dd_value(il_i1,il_j1)
264               id_detect(ji,jj)=0
265
266            ENDIF
267
268            ! lower right point
269            IF(id_detect(il_shape(1)-ji+1,jj)==1)THEN
270
271               dd_value( il_shape(1)-ji+1,jj)=dd_value(il_i2,il_j1)
272               id_detect(il_shape(1)-ji+1,jj)=0
273
274            ENDIF
275
276            ! upper left point
277            IF(id_detect(ji,il_shape(2)-jj+1)==1)THEN
278
279               dd_value( ji,il_shape(2)-jj+1)=dd_value(il_i1,il_j2)
280               id_detect(ji,il_shape(2)-jj+1)=0
281
282            ENDIF           
283
284            ! upper right point
285            IF(id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)==1)THEN
286
287               dd_value( il_shape(1)-ji+1,il_shape(2)-jj+1)=dd_value(il_i2,il_j2)
288               id_detect(il_shape(1)-ji+1,il_shape(2)-jj+1)=0
289
290            ENDIF           
291
292         ENDDO
293
294      ENDDO
295
296   END SUBROUTINE interp_nearest__2D_fill
297   !-------------------------------------------------------------------
298   !> @brief
299   !> This subroutine compute nearest interpolation of a 1D array of value.
300   !>
301   !> @author J.Paul
302   !> - September, 2014- Initial Version
303   !>
304   !> @param[inout] dd_value  1D array of mixed grid value
305   !> @param[inout] id_detect 1D array of point to be interpolated
306   !-------------------------------------------------------------------
307   SUBROUTINE interp_nearest__1D_fill( dd_value, id_detect )
308      IMPLICIT NONE
309      ! Argument
310      REAL(dp)   , DIMENSION(:), INTENT(INOUT) :: dd_value 
311      INTEGER(i4), DIMENSION(:), INTENT(INOUT) :: id_detect
312
313      ! local variable
314      INTEGER(i4), DIMENSION(1) :: il_shape
315
316      INTEGER(i4) :: il_i1
317      INTEGER(i4) :: il_i2
318
319      INTEGER(i4) :: il_half1
320     
321      ! loop indices
322      INTEGER(i4) :: ji
323      !----------------------------------------------------------------
324
325      il_shape(:)=SHAPE(dd_value)
326
327      il_i1=1
328      il_i2=il_shape(1)
329
330      il_half1=CEILING(il_shape(1)*0.5)
331     
332      DO ji=1,il_half1
333         
334         ! lower left point
335         IF(id_detect(ji)==1)THEN
336
337            dd_value( ji)=dd_value(il_i1)
338            id_detect(ji)=0
339
340         ENDIF
341
342         ! lower right point
343         IF(id_detect(il_shape(1)-ji+1)==1)THEN
344
345            dd_value( il_shape(1)-ji+1)=dd_value(il_i2)
346            id_detect(il_shape(1)-ji+1)=0
347
348         ENDIF
349
350      ENDDO
351
352   END SUBROUTINE interp_nearest__1D_fill
353END MODULE interp_nearest
Note: See TracBrowser for help on using the repository browser.