source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/TOOLS/SIREN/src/interp.f90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File size: 33.7 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 interpolation on regular grid.
10!>
11!> @details Interpolation method to be used is specify inside variable
12!>    strcuture, as array of string character.<br/>
13!>    - td_var\%c_interp(1) string character is the interpolation name choose between:
14!>       - 'nearest'
15!>       - 'cubic  '
16!>       - 'linear '
17!>    - td_var\%c_interp(2) string character is an operation to be used
18!>    on interpolated value.<br/>
19!>          operation have to be mulitplication '*' or division '/'.<br/>
20!>          coefficient have to be refinement factor following i-direction 'rhoi',
21!>          j-direction 'rhoj', or k-direction 'rhok'.<br/>
22!>
23!>          Examples: '*rhoi', '/rhoj'.
24!>
25!>    @note Those informations are read from namelist or variable configuration file (default).<br/>
26!>    Interplation method could be specify for each variable in namelist _namvar_,
27!>    defining string character _cn\_varinfo_.<br/>
28!>    Example:
29!>       - cn_varinfo='varname1:cubic/rhoi', 'varname2:linear'
30!>
31!>    to create mixed grid (with coarse grid point needed to compute
32!> interpolation):<br/>
33!> @code
34!>    CALL interp_create_mixed_grid( td_var, td_mix [,id_rho] )
35!> @endcode
36!>       - td_var is coarse grid variable (should be extrapolated)
37!>       - td_mix is mixed grid variable structure [output]
38!>       - id_rho is array of refinment factor [optional]
39!>
40!>    to detected point to be interpolated:<br/>
41!> @code
42!>    il_detect(:,:,:)=interp_detect( td_mix [,id_rho] )
43!> @endcode
44!>       - il_detect(:,:,:) is 3D array of detected point to be interpolated
45!>       - td_mix is mixed grid variable
46!>       - id_rho is array of refinement factor [optional]
47!>
48!>    to interpolate variable value:<br/>
49!> @code
50!>    CALL  interp_fill_value( td_var [,id_rho] [,id_offset] )
51!> @endcode
52!>       - td_var is variable structure
53!>       - id_rho is array of refinement factor [optional]
54!>       - id_offset is array of offset between fine and coarse grid [optional]
55!>
56!>    to clean mixed grid (remove points added on mixed grid to compute interpolation):<br/>
57!> @code
58!>    CALL interp_clean_mixed_grid( td_mix, td_var, id_rho )
59!> @endcode
60!>       - td_mix is mixed grid variable structure
61!>       - td_var is variable structure [output]
62!>       - id_rho is array of refinement factor [optional]
63!>       - id_offset is array of offset between fine and coarse grid [optional]
64!>
65!> @note It use to work on ORCA grid, as we work only with grid indices.
66!>
67!> @warning due to the use of second derivative when using cubic interpolation
68!> you should add at least 2 extrabands.
69!>
70!> @author
71!> J.Paul
72! REVISION HISTORY:
73!> @date November, 2013 - Initial Version
74!> @date September, 2014
75!> - add header
76!> - use interpolation method modules
77!>
78!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
79!----------------------------------------------------------------------
80MODULE interp
81
82   USE netcdf                          ! nf90 library
83   USE global                          ! global variable
84   USE kind                            ! F90 kind parameter
85   USE logger                          ! log file manager
86   USE fct                             ! basic useful function
87   USE date                            ! date manager
88   USE att                             ! attribute manager
89   USE dim                             ! dimension manager
90   USE var                             ! variable manager
91   USE grid                            ! grid manager
92   USE extrap                          ! extrapolation manager
93   USE interp_cubic                    !   cubic interpolation manager
94   USE interp_linear                   !  linear interpolation manager
95   USE interp_nearest                  ! nearest interpolation manager
96
97   IMPLICIT NONE
98   ! NOTE_avoid_public_variables_if_possible
99
100   ! type and variable
101
102   ! function and subroutine
103   PUBLIC :: interp_detect             !< detected point to be interpolated
104   PUBLIC :: interp_fill_value         !< interpolate value
105   PUBLIC :: interp_create_mixed_grid  !< create mixed grid
106   PUBLIC :: interp_clean_mixed_grid   !< clean mixed grid
107
108   PRIVATE :: interp__detect              ! detected point to be interpolated
109   PRIVATE :: interp__detect_wrapper      ! detected point to be interpolated
110   PRIVATE :: interp__fill_value_wrapper  ! interpolate value over detectected point
111   PRIVATE :: interp__fill_value          ! interpolate value over detectected point
112   PRIVATE :: interp__clean_even_grid     ! clean even mixed grid
113   PRIVATE :: interp__check_method        ! check if interpolation method available
114   
115   TYPE TINTERP
116      CHARACTER(LEN=lc) :: c_name   = '' !< interpolation method name
117      CHARACTER(LEN=lc) :: c_factor = '' !< interpolation factor
118      CHARACTER(LEN=lc) :: c_divisor= '' !< interpolation divisor
119   END TYPE TINTERP
120
121   INTERFACE interp_detect
122      MODULE PROCEDURE interp__detect_wrapper 
123   END INTERFACE interp_detect
124
125   INTERFACE interp_fill_value
126      MODULE PROCEDURE interp__fill_value_wrapper 
127   END INTERFACE interp_fill_value
128
129CONTAINS
130   !-------------------------------------------------------------------
131   !> @brief
132   !> This function check if interpolation method is available.
133   !>
134   !> @details
135   !> check if name of interpolation method is present in global list of string
136   !> character cp_interp_list (see global.f90).
137   !>
138   !> @author J.Paul
139   !> - November, 2013- Initial Version
140   !
141   !> @param[in] cd_method interpolation method
142   !> @return
143   !-------------------------------------------------------------------
144   FUNCTION interp__check_method( cd_method )
145      IMPLICIT NONE
146      ! Argument
147      CHARACTER(LEN=lc) :: cd_method
148
149      ! function
150      LOGICAL :: interp__check_method
151
152      ! local variable
153      CHARACTER(LEN=lc) :: cl_interp
154      CHARACTER(LEN=lc) :: cl_method
155
156      ! loop indices
157      INTEGER(I4) :: ji
158      !----------------------------------------------------------------
159
160      cl_method=fct_lower(cd_method)
161
162      interp__check_method=.FALSE.
163      DO ji=1,ip_ninterp
164         cl_interp=fct_lower(cp_interp_list(ji))
165         IF( TRIM(cl_interp) == TRIM(cl_method) )THEN
166            interp__check_method=.TRUE.
167            EXIT
168         ENDIF
169      ENDDO
170
171   END FUNCTION interp__check_method
172   !-------------------------------------------------------------------
173   !> @brief
174   !> This function detected point to be interpolated.
175   !>
176   !> @details
177   !> Actually it checks, the number of dimension used for this variable
178   !> and launch interp__detect which detected point to be interpolated.
179   !>
180   !> @author J.Paul
181   !> - November, 2013- Initial Version
182   !>
183   !> @param[in] td_mix mixed grid variable (to interpolate)
184   !> @param[in] id_rho array of refinement factor
185   !> @return 3D array of detected point to be interpolated
186   !-------------------------------------------------------------------
187   FUNCTION interp__detect_wrapper( td_mix, id_rho )
188      IMPLICIT NONE
189      ! Argument
190      TYPE(TVAR) , INTENT(IN) :: td_mix
191      INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
192
193      ! function
194      INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
195      &                      td_mix%t_dim(2)%i_len,&
196      &                      td_mix%t_dim(3)%i_len ) :: interp__detect_wrapper
197
198      ! local variable
199      ! loop indices
200      !----------------------------------------------------------------
201
202      IF( .NOT. ANY(td_mix%t_dim(1:3)%l_use) )THEN
203         ! no dimension I-J-K used
204         CALL logger_debug(" INTERP DETECT: nothing done for variable"//&
205         &              TRIM(td_mix%c_name) )
206         
207         interp__detect_wrapper(:,:,:)=0
208
209      ELSE IF( ALL(td_mix%t_dim(1:3)%l_use) )THEN
210         
211         ! detect point to be interpolated on I-J-K
212         CALL logger_debug(" INTERP DETECT: detect point "//&
213         &              TRIM(td_mix%c_point)//" for variable "//&
214         &              TRIM(td_mix%c_name) )
215         
216         interp__detect_wrapper(:,:,:)=interp__detect( td_mix, id_rho(:) )
217
218      ELSE IF( ALL(td_mix%t_dim(1:2)%l_use) )THEN
219
220         ! detect point to be interpolated on I-J
221         CALL logger_debug(" INTERP DETECT: detect point "//&
222         &              TRIM(td_mix%c_point)//" for variable "//&
223         &              TRIM(td_mix%c_name) )
224         
225         interp__detect_wrapper(:,:,1:1)=interp__detect( td_mix, id_rho(:))
226
227      ELSE IF( td_mix%t_dim(3)%l_use )THEN
228         
229         ! detect point to be interpolated on K
230         CALL logger_debug(" INTERP DETECT: detect vertical point "//&
231         &              " for variable "//TRIM(td_mix%c_name) )
232         
233         interp__detect_wrapper(1:1,1:1,:)=interp__detect( td_mix, id_rho(:) )
234
235      ENDIF             
236
237   END FUNCTION interp__detect_wrapper
238   !-------------------------------------------------------------------
239   !> @brief
240   !> This function detected point to be interpolated.
241   !>
242   !> @details
243   !> A special case is done for even refinement on ARAKAWA-C grid.
244   !>
245   !> @author J.Paul
246   !> - November, 2013- Initial Version
247   !
248   !> @param[in] td_mix mixed grid variable (to interpolate)
249   !> @param[in] id_rho array of refinement factor
250   !> @return 3D array of detected point to be interpolated
251   !-------------------------------------------------------------------
252   FUNCTION interp__detect( td_mix, id_rho )
253      IMPLICIT NONE
254      ! Argument
255      TYPE(TVAR) , INTENT(IN) :: td_mix
256      INTEGER(I4), DIMENSION(:), INTENT(IN), OPTIONAL :: id_rho
257
258      ! function
259      INTEGER(i4), DIMENSION(td_mix%t_dim(1)%i_len,&
260      &                      td_mix%t_dim(2)%i_len,&
261      &                      td_mix%t_dim(3)%i_len ) :: interp__detect
262
263      ! local variable
264      INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
265
266      INTEGER(I4) :: il_xextra
267      INTEGER(I4) :: il_yextra
268      INTEGER(I4) :: il_zextra
269
270      INTEGER(i4), DIMENSION(3) :: il_dim
271
272      LOGICAL    , DIMENSION(3) :: ll_even
273
274      ! loop indices
275      INTEGER(I4) :: ji
276      INTEGER(I4) :: jj
277      INTEGER(I4) :: jk
278      !----------------------------------------------------------------
279      ALLOCATE( il_rho(ip_maxdim) )
280      il_rho(:)=1
281      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
282
283      ! special case for even refinement on ARAKAWA-C grid
284      ll_even(:)=.FALSE.
285      IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
286      IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
287      IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.         
288
289      SELECT CASE(TRIM(td_mix%c_point))
290         CASE('U')
291            ll_even(1)=.FALSE.
292         CASE('V')
293            ll_even(2)=.FALSE.
294         CASE('F')
295            ll_even(:)=.FALSE.
296      END SELECT
297
298      IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
299      IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
300      IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
301
302      ! special case for cubic interpolation
303      il_xextra=0
304      il_yextra=0
305      il_zextra=0
306      SELECT CASE(TRIM(td_mix%c_interp(1)))
307      CASE('cubic')
308         ! those points can not be compute cause cubic interpolation
309         ! need second derivative.
310         IF( il_rho(jp_I) /= 1 ) il_xextra=3*il_rho(jp_I)
311         IF( il_rho(jp_J) /= 1 ) il_yextra=3*il_rho(jp_J)
312         IF( il_rho(jp_K) /= 1 ) il_zextra=3*il_rho(jp_K)
313      END SELECT
314
315      il_dim(:)=td_mix%t_dim(1:3)%i_len
316
317      ! init
318      interp__detect(:,:,:)=1
319
320      ! do not compute coarse grid point
321      interp__detect(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
322      &              1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
323      &              1:td_mix%t_dim(3)%i_len:il_rho(jp_K)  ) = 0
324
325      ! do not compute point near fill value
326      DO jk=1,il_dim(3),il_rho(jp_K)
327         DO jj=1,il_dim(2),il_rho(jp_J)
328            DO ji=1,il_dim(1),il_rho(jp_I)
329
330               IF( td_mix%d_value(ji,jj,jk,1) == td_mix%d_fill )THEN
331
332                  ! i-direction
333                  interp__detect(MAX(1,ji-il_xextra):MIN(ji+il_xextra,il_dim(1)),&
334                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
335                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0 
336                  ! j-direction
337                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
338                  &              MAX(1,jj-il_yextra):MIN(jj+il_yextra,il_dim(2)),&
339                  &              MAX(1,jk-(il_rho(jp_K)-1)):MIN(jk+(il_rho(jp_K)-1),il_dim(3)) )=0
340                  ! k-direction
341                  interp__detect(MAX(1,ji-(il_rho(jp_I)-1)):MIN(ji+(il_rho(jp_I)-1),il_dim(1)),&
342                  &              MAX(1,jj-(il_rho(jp_J)-1)):MIN(jj+(il_rho(jp_J)-1),il_dim(2)),&
343                  &              MAX(1,jk-il_zextra):MIN(jk+il_zextra,il_dim(3)) )=0         
344
345               ENDIF
346
347            ENDDO
348         ENDDO
349      ENDDO
350
351      DEALLOCATE( il_rho )
352
353   END FUNCTION interp__detect
354   !-------------------------------------------------------------------
355   !> @brief
356   !> This subroutine create mixed grid.
357   !>
358   !> @details
359   !> Created grid is fine resolution grid.
360   !> First and last point are coasre grid point.
361   !>
362   !> A special case is done for even refinement on ARAKAWA-C grid.
363   !>
364   !> @author J.Paul
365   !> - November, 2013- Initial Version
366   !>
367   !> @param[in] td_var    coarse grid variable (should be extrapolated)
368   !> @param[out] td_mix   mixed grid variable
369   !> @param[in] id_rho    array of refinment factor (default 1)
370   !-------------------------------------------------------------------
371   SUBROUTINE interp_create_mixed_grid( td_var, td_mix, id_rho )
372      IMPLICIT NONE
373      ! Argument
374      TYPE(TVAR) ,               INTENT(IN   ) :: td_var 
375      TYPE(TVAR) ,               INTENT(  OUT) :: td_mix 
376      INTEGER(I4), DIMENSION(:), INTENT(IN   ), OPTIONAL :: id_rho
377
378      ! local variable
379      INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
380
381      INTEGER(i4) :: il_xextra
382      INTEGER(i4) :: il_yextra
383      INTEGER(i4) :: il_zextra
384
385      LOGICAL, DIMENSION(3) :: ll_even
386
387      ! loop indices
388      !----------------------------------------------------------------
389      ALLOCATE(il_rho(ip_maxdim))
390      il_rho(:)=1
391      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
392
393      ! special case for even refinement on ARAKAWA-C grid
394      ll_even(:)=.FALSE.
395      IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
396      IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
397      IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.         
398
399      SELECT CASE(TRIM(td_var%c_point))
400         CASE('U')
401            ll_even(1)=.FALSE.
402         CASE('V')
403            ll_even(2)=.FALSE.
404         CASE('F')
405            ll_even(:)=.FALSE.
406      END SELECT
407
408      IF( ll_even(1) ) il_rho(jp_I)=il_rho(jp_I)+1
409      IF( ll_even(2) ) il_rho(jp_J)=il_rho(jp_J)+1
410      IF( ll_even(3) ) il_rho(jp_K)=il_rho(jp_K)+1
411
412      ! copy variable
413      td_mix=var_copy(td_var)
414
415      ! compute new dimension length
416      il_xextra=il_rho(jp_I)-1
417      td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
418
419      il_yextra=il_rho(jp_J)-1
420      td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
421
422      il_zextra=il_rho(jp_K)-1
423      td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
424
425      IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
426      ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
427      &                         td_mix%t_dim(2)%i_len, &
428      &                         td_mix%t_dim(3)%i_len, &
429      &                         td_mix%t_dim(4)%i_len) )
430
431      ! initialise to FillValue
432      td_mix%d_value(:,:,:,:)=td_mix%d_fill
433
434      ! quid qd coord ou bathy fourni par user ?? (offset ??)
435      td_mix%d_value(1:td_mix%t_dim(1)%i_len:il_rho(jp_I), &
436      &              1:td_mix%t_dim(2)%i_len:il_rho(jp_J), &
437      &              1:td_mix%t_dim(3)%i_len:il_rho(jp_K), :) = &
438      &     td_var%d_value(:,:,:,:)
439
440      DEALLOCATE(il_rho)
441
442   END SUBROUTINE interp_create_mixed_grid
443   !-------------------------------------------------------------------
444   !> @brief
445   !> This subroutine remove points added to mixed grid to compute
446   !> interpolation in the special case of even refinement on ARAKAWA-C grid.
447   !>
448   !> @details
449   !>
450   !> @author J.Paul
451   !> - November, 2013- Initial Version
452   !>
453   !> @param[inout] td_mix mixed grid variable
454   !> @param[in] id_rho    array of refinment factor
455   !-------------------------------------------------------------------
456   SUBROUTINE interp__clean_even_grid( td_mix, id_rho )
457      IMPLICIT NONE
458      ! Argument
459      TYPE(TVAR) ,               INTENT(INOUT) :: td_mix 
460      INTEGER(I4), DIMENSION(:), INTENT(IN   ), OPTIONAL :: id_rho
461
462      ! local variable
463      INTEGER(I4), DIMENSION(:), ALLOCATABLE :: il_rho
464
465      INTEGER(i4) :: il_xextra
466      INTEGER(i4) :: il_yextra
467      INTEGER(i4) :: il_zextra
468
469      LOGICAL, DIMENSION(3) :: ll_even
470
471      LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ll_mask
472
473      REAL(dp), DIMENSION(:)     , ALLOCATABLE :: dl_vect
474
475      TYPE(TVAR) :: tl_mix
476
477      ! loop indices
478      !----------------------------------------------------------------
479      ALLOCATE(il_rho(ip_maxdim))
480      il_rho(:)=1
481      IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:)
482
483      ! special case for even refinement on ARAKAWA-C grid
484      ll_even(:)=.FALSE.
485      IF( MOD(il_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
486      IF( MOD(il_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
487      IF( MOD(il_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.         
488
489      SELECT CASE(TRIM(td_mix%c_point))
490         CASE('U')
491            ll_even(1)=.FALSE.
492         CASE('V')
493            ll_even(2)=.FALSE.
494         CASE('F')
495            ll_even(:)=.FALSE.
496      END SELECT
497
498      ! remove some point only if refinement in some direction is even
499      IF( ANY(ll_even(:)) )THEN
500
501         ! copy variable
502         tl_mix=var_copy(td_mix)
503
504         ALLOCATE( ll_mask( tl_mix%t_dim(1)%i_len, &
505         &                  tl_mix%t_dim(2)%i_len, &
506         &                  tl_mix%t_dim(3)%i_len, &
507         &                  tl_mix%t_dim(4)%i_len) )
508
509         ll_mask(:,:,:,:)=.TRUE.
510
511         IF( tl_mix%t_dim(1)%l_use .AND. ll_even(1) )THEN
512
513            il_rho(jp_I)=il_rho(jp_I)+1
514
515            ! locate wrong point on mixed grid
516            ll_mask(1:td_mix%t_dim(1)%i_len:il_rho(jp_I),:,:,:)=.FALSE.
517
518            ! compute coasre grid dimension length
519            il_xextra=il_rho(jp_I)-1
520            td_mix%t_dim(1)%i_len=(tl_mix%t_dim(1)%i_len+il_xextra)/il_rho(jp_I)
521
522            il_rho(jp_I)=il_rho(jp_I)-1
523            ! compute right fine grid dimension length
524            td_mix%t_dim(1)%i_len=td_mix%t_dim(1)%i_len*il_rho(jp_I)-il_xextra
525
526         ENDIF
527
528         IF( tl_mix%t_dim(2)%l_use .AND. ll_even(2) )THEN
529
530            il_rho(jp_J)=il_rho(jp_J)+1
531
532            ! locate wrong point on mixed grid
533            ll_mask(:,1:tl_mix%t_dim(2)%i_len:il_rho(jp_J),:,:)=.FALSE.
534
535            ! compute coasre grid dimension length
536            il_yextra=il_rho(jp_J)-1
537            td_mix%t_dim(2)%i_len=(tl_mix%t_dim(2)%i_len+il_yextra)/il_rho(jp_J)
538
539            il_rho(jp_J)=il_rho(jp_J)-1
540            ! compute right fine grid dimension length
541            td_mix%t_dim(2)%i_len=td_mix%t_dim(2)%i_len*il_rho(jp_J)-il_yextra
542           
543         ENDIF
544
545         IF( tl_mix%t_dim(3)%l_use .AND. ll_even(3) )THEN
546
547            il_rho(jp_K)=il_rho(jp_K)+1
548
549            ! locate wrong point on mixed grid
550            ll_mask(:,:,1:tl_mix%t_dim(3)%i_len:il_rho(jp_K),:)=.FALSE.
551
552            ! compute coasre grid dimension length
553            il_zextra=il_rho(jp_K)-1
554            td_mix%t_dim(3)%i_len=(tl_mix%t_dim(3)%i_len+il_zextra)/il_rho(jp_K)
555
556            il_rho(jp_K)=il_rho(jp_K)-1
557            ! compute right fine grid dimension length
558            td_mix%t_dim(3)%i_len=td_mix%t_dim(3)%i_len*il_rho(jp_K)-il_zextra
559
560         ENDIF     
561
562         IF( ASSOCIATED(td_mix%d_value) ) DEALLOCATE( td_mix%d_value )
563         ALLOCATE( td_mix%d_value( td_mix%t_dim(1)%i_len, &
564         &                         td_mix%t_dim(2)%i_len, &
565         &                         td_mix%t_dim(3)%i_len, &
566         &                         td_mix%t_dim(4)%i_len) )
567
568         ! initialise to FillValue
569         td_mix%d_value(:,:,:,:)=td_mix%d_fill
570
571         IF( COUNT(ll_mask(:,:,:,:)) /= SIZE(td_mix%d_value(:,:,:,:)) )THEN
572
573            CALL logger_error("INTERP CLEAN EVEN GRID: output value size "//&
574            &  " and mask count differ ")
575
576         ELSE
577
578            ALLOCATE( dl_vect(COUNT(ll_mask(:,:,:,:))) )
579
580            dl_vect(:)= PACK( tl_mix%d_value(:,:,:,:), &
581            &                 MASK=ll_mask(:,:,:,:)     )
582
583            td_mix%d_value(:,:,:,:)=RESHAPE( dl_vect(:), &
584            &                                SHAPE=(/td_mix%t_dim(1)%i_len, &
585            &                                        td_mix%t_dim(2)%i_len, &
586            &                                        td_mix%t_dim(3)%i_len, &
587            &                                        td_mix%t_dim(4)%i_len/) )
588
589            DEALLOCATE( dl_vect )
590
591         ENDIF
592
593         DEALLOCATE( ll_mask )
594
595         ! clean
596         CALL var_clean(tl_mix)
597
598      ENDIF
599
600      ! clean
601      DEALLOCATE(il_rho)
602
603   END SUBROUTINE interp__clean_even_grid
604   !-------------------------------------------------------------------
605   !> @brief
606   !> This subroutine remove points added on mixed grid
607   !> to compute interpolation. And save interpolated value over domain.
608   !>
609   !> @details
610   !>
611   !> @author J.Paul
612   !> - November, 2013- Initial Version
613   !> @date September, 2014
614   !> - use offset to save useful domain
615   !>
616   !> @param[in] td_mix    mixed grid variable structure
617   !> @param[out] td_var   variable structure
618   !> @param[in] id_rho    array of refinement factor (default 1)
619   !> @param[in] id_offset 2D array of offset between fine and coarse grid
620   !-------------------------------------------------------------------
621   SUBROUTINE interp_clean_mixed_grid( td_mix, td_var, &
622   &                                   id_rho, id_offset )
623      IMPLICIT NONE
624      ! Argument
625      TYPE(TVAR)                 , INTENT(IN   ) :: td_mix
626      TYPE(TVAR)                 , INTENT(  OUT) :: td_var
627      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ) :: id_rho
628      INTEGER(I4), DIMENSION(2,2), INTENT(IN   ) :: id_offset
629
630      ! local variable
631      INTEGER(i4) :: il_imin0
632      INTEGER(i4) :: il_imax0
633      INTEGER(i4) :: il_jmin0
634      INTEGER(i4) :: il_jmax0
635
636      INTEGER(i4) :: il_imin1
637      INTEGER(i4) :: il_jmin1
638      INTEGER(i4) :: il_imax1
639      INTEGER(i4) :: il_jmax1
640
641      REAL(dp), DIMENSION(:,:,:,:) , ALLOCATABLE :: dl_value
642
643      TYPE(TVAR) :: tl_mix
644
645      ! loop indices
646      !----------------------------------------------------------------
647      ! copy mixed variable in temporary structure
648      tl_mix=var_copy(td_mix)
649
650      ! remove unusefull points over mixed grid for even refinement
651      CALL interp__clean_even_grid(tl_mix, id_rho(:))
652
653      ! copy cleaned mixed variable
654      td_var=var_copy(tl_mix)
655
656      ! delete array of value
657      CALL var_del_value(td_var)
658
659      ! compute domain indices in i-direction
660      il_imin0=1  ; il_imax0=td_var%t_dim(1)%i_len
661 
662      IF( td_var%t_dim(1)%l_use )THEN
663         il_imin1=il_imin0+id_offset(Jp_I,1)
664         il_imax1=il_imax0-id_offset(Jp_I,2)
665      ELSE
666
667         il_imin1=il_imin0
668         il_imax1=il_imax0
669
670      ENDIF
671
672      ! compute domain indices in j-direction
673      il_jmin0=1  ; il_jmax0=td_var%t_dim(2)%i_len
674 
675      IF( td_var%t_dim(2)%l_use )THEN
676         il_jmin1=il_jmin0+id_offset(Jp_J,1)
677         il_jmax1=il_jmax0-id_offset(Jp_J,2)
678      ELSE
679
680         il_jmin1=il_jmin0
681         il_jmax1=il_jmax0
682 
683      ENDIF
684
685      ! compute new dimension
686      td_var%t_dim(1)%i_len=il_imax1-il_imin1+1
687      td_var%t_dim(2)%i_len=il_jmax1-il_jmin1+1
688 
689      ALLOCATE(dl_value(td_var%t_dim(1)%i_len, &
690      &                 td_var%t_dim(2)%i_len, &
691      &                 td_var%t_dim(3)%i_len, &
692      &                 td_var%t_dim(4)%i_len) )
693
694      dl_value( 1:td_var%t_dim(1)%i_len, &
695      &         1:td_var%t_dim(2)%i_len, &
696      &         :,:) = tl_mix%d_value( il_imin1:il_imax1, &
697      &                                il_jmin1:il_jmax1, &
698      &                                     :, : )
699
700      ! add variable value
701      CALL var_add_value(td_var,dl_value(:,:,:,:))
702
703      DEALLOCATE(dl_value)
704 
705      ! clean
706      CALL var_clean(tl_mix)
707
708   END SUBROUTINE interp_clean_mixed_grid
709   !-------------------------------------------------------------------
710   !> @brief
711   !> This subroutine interpolate variable value.
712   !>
713   !> @details
714   !> Actually it checks, the number of dimension used for this variable
715   !> and launch interp__fill_value.
716   !>
717   !> @author J.Paul
718   !> - November, 2013- Initial Version
719   !>
720   !> @param[inout] td_var variable structure
721   !> @param[in] id_rho    array of refinement factor
722   !> @param[in] id_offset 2D array of offset between fine and coarse grid
723   !-------------------------------------------------------------------
724   SUBROUTINE interp__fill_value_wrapper( td_var, & 
725   &                                      id_rho, &
726   &                                      id_offset )
727      IMPLICIT NONE
728      ! Argument
729      TYPE(TVAR) ,                 INTENT(INOUT) :: td_var
730      INTEGER(I4), DIMENSION(:)  , INTENT(IN   ), OPTIONAL :: id_rho
731      INTEGER(I4), DIMENSION(:,:), INTENT(IN   ), OPTIONAL :: id_offset
732
733      ! local variable
734      CHARACTER(LEN=lc)                            :: cl_method
735      INTEGER(i4)      , DIMENSION(:), ALLOCATABLE :: il_rho
736      INTEGER(i4)      , DIMENSION(2,2)            :: il_offset
737
738      ! loop indices
739      !----------------------------------------------------------------
740
741      ALLOCATE( il_rho(ip_maxdim) )
742      il_rho(:)=1
743      IF( PRESENT(id_rho) ) il_rho(1:SIZE(id_rho(:)))=id_rho(:)
744      IF( ANY(il_rho(:) < 0) )THEN
745         CALL logger_error("INTERP FILL VALUE: invalid "//&
746         &  " refinement factor ")
747      ENDIF
748
749      il_offset(:,:)=0
750      IF( PRESENT(id_offset) )THEN
751         IF( ANY(SHAPE(id_offset(:,:)) /= (/2,2/)) )THEN
752            CALL logger_error("INTERP FILL VALUE: invalid array of offset")
753         ELSE
754            il_offset(:,:)=id_offset(:,:)
755         ENDIF
756      ENDIF
757
758      IF( (il_rho(jp_I) /= 1 .AND. td_var%t_dim(1)%l_use) .OR. &
759      &   (il_rho(jp_J) /= 1 .AND. td_var%t_dim(2)%l_use) .OR. &
760      &   (il_rho(jp_K) /= 1 .AND. td_var%t_dim(3)%l_use) )THEN
761
762         SELECT CASE(TRIM(td_var%c_interp(1)))
763            CASE('cubic','linear','nearest')
764               cl_method=TRIM(td_var%c_interp(1))
765            CASE DEFAULT
766               CALL logger_warn("INTERP FILL VALUE: interpolation method unknown."//&
767               &  " use linear interpolation")
768               cl_method='linear'
769               ! update variable structure value
770               td_var%c_interp(1)='linear'
771         END SELECT
772
773         CALL logger_info("INTERP FILL: interpolate "//TRIM(td_var%c_name)//&
774         &             " using "//TRIM(cl_method)//" method." )
775         CALL logger_info("INTERP FILL: refinement factor "//&
776         &                        TRIM(fct_str(il_rho(jp_I)))//&
777         &                   " "//TRIM(fct_str(il_rho(jp_J)))//&
778         &                   " "//TRIM(fct_str(il_rho(jp_K))) )
779 
780         CALL interp__fill_value( td_var, cl_method, & 
781         &                        il_rho(:), il_offset(:,:) )
782
783         SELECT CASE(TRIM(td_var%c_interp(2)))
784         CASE('/rhoi')
785            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
786               td_var%d_value(:,:,:,:) = &
787               &  td_var%d_value(:,:,:,:) / REAL(il_rho(jp_I),dp)
788            END WHERE
789         CASE('/rhoj')
790            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
791               td_var%d_value(:,:,:,:) = &
792               &  td_var%d_value(:,:,:,:) / REAL(il_rho(jp_J),dp)
793            END WHERE
794         CASE('/rhok')
795            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
796               td_var%d_value(:,:,:,:) = &
797               &  td_var%d_value(:,:,:,:) / REAL(il_rho(jp_K),dp)
798            END WHERE
799         CASE('*rhoi')
800            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill )
801               td_var%d_value(:,:,:,:) = &
802               &  td_var%d_value(:,:,:,:) * REAL(il_rho(jp_I),dp)
803            END WHERE
804         CASE('*rhoj')
805            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) 
806               td_var%d_value(:,:,:,:) = &
807               &  td_var%d_value(:,:,:,:) * REAL(il_rho(jp_J),dp)
808            END WHERE
809         CASE('*rhok')
810            WHERE( td_var%d_value(:,:,:,:) /= td_var%d_fill ) 
811               td_var%d_value(:,:,:,:) = &
812               &  td_var%d_value(:,:,:,:) * REAL(il_rho(jp_K),dp)
813            END WHERE
814         CASE DEFAULT
815            td_var%c_interp(2)=''
816         END SELECT
817
818      ELSE
819         td_var%c_interp(:)=''
820      ENDIF
821
822      DEALLOCATE(il_rho)
823
824   END SUBROUTINE interp__fill_value_wrapper
825   !-------------------------------------------------------------------
826   !> @brief
827   !> This subroutine interpolate value over mixed grid.
828   !>
829   !> @author J.Paul
830   !> - November, 2013- Initial Version
831   !> @date September, 2014
832   !> - use interpolation method modules
833   !>
834   !> @param[inout] td_var variable structure
835   !> @param[in] cd_method interpolation method
836   !> @param[in] id_rho    array of refinment factor
837   !> @param[in] id_offset 2D array of offset between fine and coarse grid
838   !-------------------------------------------------------------------
839   SUBROUTINE interp__fill_value( td_var, cd_method, &
840   &                              id_rho, id_offset )
841      IMPLICIT NONE
842      ! Argument
843      TYPE(TVAR)                      , INTENT(INOUT) :: td_var 
844      CHARACTER(LEN=*)                , INTENT(IN   ) :: cd_method
845      INTEGER(I4)     , DIMENSION(:)  , INTENT(IN   ) :: id_rho
846      INTEGER(I4)     , DIMENSION(2,2), INTENT(IN   ) :: id_offset
847
848      ! local variable
849      CHARACTER(LEN=lc)                                :: cl_interp
850
851      INTEGER(I4)    , DIMENSION(:)      , ALLOCATABLE :: il_rho
852      INTEGER(i4)    , DIMENSION(:,:,:)  , ALLOCATABLE :: il_detect
853 
854      REAL(dp)                                         :: dl_min
855      REAL(dp)                                         :: dl_max
856
857      LOGICAL        , DIMENSION(3)                    :: ll_even
858      LOGICAL                                          :: ll_discont
859 
860      TYPE(TVAR)                                       :: tl_mix
861
862      TYPE(TATT)                                       :: tl_att
863 
864      ! loop indices
865      !----------------------------------------------------------------
866
867      !1- create mixed grid
868      CALL interp_create_mixed_grid( td_var, tl_mix, id_rho(:) )
869
870      ! clean variable structure
871      CALL var_clean(td_var)
872
873      !2- detect point to be interpolated
874      ALLOCATE( il_detect( tl_mix%t_dim(1)%i_len, &
875      &                    tl_mix%t_dim(2)%i_len, &
876      &                    tl_mix%t_dim(3)%i_len) )
877
878      il_detect(:,:,:)=0
879
880      il_detect(:,:,:)=interp_detect(tl_mix, id_rho(:) )
881
882      ! add attribute to variable
883      cl_interp=fct_concat(tl_mix%c_interp(:))
884      tl_att=att_init('interpolation',cl_interp)
885      CALL var_move_att(tl_mix, tl_att)
886
887      ! clean
888      CALL att_clean(tl_att)
889
890      ! special case for even refinement on ARAKAWA-C grid
891      ll_even(:)=.FALSE.
892      IF( MOD(id_rho(jp_I),2) == 0 ) ll_even(1)=.TRUE.
893      IF( MOD(id_rho(jp_J),2) == 0 ) ll_even(2)=.TRUE.
894      IF( MOD(id_rho(jp_K),2) == 0 ) ll_even(3)=.TRUE.         
895
896      SELECT CASE(TRIM(tl_mix%c_point))
897         CASE('U')
898            ll_even(1)=.FALSE.
899         CASE('V')
900            ll_even(2)=.FALSE.
901         CASE('F')
902            ll_even(:)=.FALSE.
903      END SELECT
904
905      ALLOCATE(il_rho(ip_maxdim))
906      il_rho(:)=id_rho(:)
907
908      IF( ll_even(1) ) il_rho(jp_I)=id_rho(jp_I)+1
909      IF( ll_even(2) ) il_rho(jp_J)=id_rho(jp_J)+1
910      IF( ll_even(3) ) il_rho(jp_K)=id_rho(jp_K)+1
911
912      ! special case for longitude
913      ll_discont=.FALSE.
914      IF( TRIM(tl_mix%c_units) == 'degrees_east' )THEN
915         dl_min=MINVAL( tl_mix%d_value(:,:,:,:), &
916         &              tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)
917         dl_max=MAXVAL( tl_mix%d_value(:,:,:,:), &
918         &              tl_mix%d_value(:,:,:,:)/=tl_mix%d_fill)         
919         IF( dl_min < -170_dp .AND. dl_max > 170_dp .OR. &
920         &   dl_min <   10_dp .AND. dl_max > 350_dp )THEN
921            ll_discont=.TRUE.
922         ENDIF
923      ENDIF
924
925      !3- interpolate
926      CALL logger_debug("INTERP 2D: interpolation method "//TRIM(cd_method)//&
927      &  " discont "//TRIM(fct_str(ll_discont)) )
928      SELECT CASE(TRIM(cd_method))
929      CASE('cubic')
930         CALL interp_cubic_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
931           &                    il_detect(:,:,:),                        &
932           &                    il_rho(:), ll_even(:), ll_discont )
933      CASE('nearest')
934         CALL interp_nearest_fill(tl_mix%d_value(:,:,:,:), &
935              &                   il_detect(:,:,:),         &
936              &                   il_rho(:) )
937      CASE DEFAULT ! linear
938         CALL interp_linear_fill(tl_mix%d_value(:,:,:,:), tl_mix%d_fill, &
939              &                  il_detect(:,:,:),                        &
940              &                  il_rho(:), ll_even(:), ll_discont )
941      END SELECT         
942
943      IF( ANY(il_detect(:,:,:)==1) )THEN
944         CALL logger_warn("INTERP FILL: some points can not be interpolated "//&
945         &             "for variable "//TRIM(tl_mix%c_name) )
946      ENDIF
947
948      DEALLOCATE(il_detect)
949      !4- save useful domain (remove offset)
950      CALL interp_clean_mixed_grid( tl_mix, td_var, &
951      &                             id_rho(:), id_offset(:,:)  ) 
952
953      ! clean variable structure
954      DEALLOCATE(il_rho)
955      CALL var_clean(tl_mix)
956
957   END SUBROUTINE interp__fill_value
958END MODULE interp
Note: See TracBrowser for help on using the repository browser.