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 | |
---|
31 | MODULE 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 | |
---|
57 | CONTAINS |
---|
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._dp/REAL(id_rho(jp_I)-1,dp) |
---|
630 | ELSE ! odd refinement |
---|
631 | dl_dx=1._dp/REAL(id_rho(jp_I),dp) |
---|
632 | ENDIF |
---|
633 | |
---|
634 | IF( ld_even(jp_J) )THEN |
---|
635 | dl_dy=1._dp/REAL(id_rho(jp_J)-1,dp) |
---|
636 | ELSE ! odd refinement |
---|
637 | dl_dy=1._dp/REAL(id_rho(jp_J),dp) |
---|
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=REAL(jj-1,dp)*dl_dy - dl_dy*0.5_dp |
---|
645 | ELSE ! odd refinement |
---|
646 | dl_y=REAL(jj-1,dp)*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=REAL(ji-1,dp)*dl_dx - dl_dx*0.5_dp |
---|
656 | ELSE ! odd refinement |
---|
657 | dl_x=REAL(ji-1,dp)*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._dp/REAL(id_rho-1,dp) |
---|
695 | ELSE ! odd refinement |
---|
696 | dl_dx=1._dp/REAL(id_rho,dp) |
---|
697 | ENDIF |
---|
698 | |
---|
699 | DO ji=1,id_rho+1 |
---|
700 | IF( ld_even )THEN |
---|
701 | dl_x=REAL(ji-1,dp)*dl_dx - dl_dx*0.5_dp |
---|
702 | ELSE ! odd refinement |
---|
703 | dl_x=REAL(ji-1,dp)*dl_dx |
---|
704 | ENDIF |
---|
705 | |
---|
706 | dd_weight(:,ji)=(/1._dp, dl_x /) |
---|
707 | ENDDO |
---|
708 | |
---|
709 | END SUBROUTINE interp_linear__get_weight1D |
---|
710 | END MODULE interp_linear |
---|