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