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

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

update nemo trunk

File size: 44.0 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! DESCRIPTION:
6!> @brief This module is filter manager.
7!>
8!> @details Filtering method to be used is specify inside variable strcuture,
9!>    as array of string character.<br/>
10!>    td_var\%c_filter(1) string character is the filter name choose between:<br/>
11!>       - 'hann'
12!>          - rad < cutoff : @f$ filter=0.5+0.5*COS(\pi*\frac{rad}{cutoff}) @f$
13!>          - rad > cutoff : @f$ filter=0 @f$
14!>       - 'hamming'
15!>          - rad < cutoff : @f$ filter=0.54+0.46*COS(\pi*\frac{rad}{cutoff}) @f$
16!>          - rad > cutoff : @f$ filter=0 @f$               
17!>       - 'blackman'
18!>          - rad < cutoff : @f$ filter=0.42 + 0.5*COS(\pi*\frac{rad}{cutoff}) +
19!>                                      0.08*COS(2\pi*\frac{rad}{cutoff}) @f$
20!>          - rad > cutoff : @f$ filter=0 @f$
21!>       - 'gauss'
22!>          - @f$filter=exp(-(\alpha * rad^2) / (2*cutoff^2))@f$
23!>       - 'butterworth'
24!>          - @f$ filer=1 / (1+(rad^2 / cutoff^2)^{\alpha}) @f$
25!>             .
26!>
27!>       with @f$ rad= \sqrt{(dist-radius)^2} @f$
28!>
29!>    td_var\%c_filter(2) string character is the number of turn to be done<br/>
30!>    td_var\%c_filter(3) string character is the cut-off frequency
31! >                       (count in number of mesh grid)<br/>
32!>    td_var\%c_filter(4) string character is the halo radius
33!>                        (count in number of mesh grid)<br/>
34!>    td_var\%c_filter(5) string character is the alpha parameter
35!>                        (for gauss and butterworth method)<br/>
36!>   
37!>    @note Filter method could be specify for each variable in namelist _namvar_,
38!>    defining string character _cn\_varinfo_. None by default.<br/>
39!>    Filter method parameters are informed inside bracket.
40!>       - @f$\alpha@f$ parameter is added for _gauss_ and _butterworth_ methods
41!>
42!>    The number of turn is specify using '*' separator.<br/>
43!>    Example:
44!>       - cn_varinfo='varname1:flt=2*hamming(@f$cutoff@f$,@f$radius@f$)',
45!>                    'varname2:flt=gauss(@f$cutoff@f$,@f$radius@f$,@f$\alpha@f$)'
46!>
47!>    to filter variable value:<br/>
48!> @code
49!>    CALL filter_fill_value( td_var )
50!> @endcode
51!>       - td_var is variable structure
52!>
53!> @author
54!> J.Paul
55!>
56!> @date November, 2013 - Initial Version
57!>
58!> @note Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59!----------------------------------------------------------------------
60MODULE filter
61
62   USE kind                            ! F90 kind parameter
63   USE phycst                          ! physical constant
64   USE logger                          ! log file manager
65   USE fct                             ! basic usefull function
66   use att                             ! attribute manager
67   USE var                             ! variable manager
68   USE extrap                          ! extrapolation manager
69
70   IMPLICIT NONE
71   ! NOTE_avoid_public_variables_if_possible
72
73   ! type and variable
74
75
76   ! function and subroutine
77   PUBLIC :: filter_fill_value   !< filter variable value
78
79   PRIVATE :: filter__fill_value_wrapper !
80   PRIVATE :: filter__fill_value         !
81   PRIVATE :: filter__3D_fill_value      !
82   PRIVATE :: filter__2D_fill_value      !
83   PRIVATE :: filter__2D                 !
84   PRIVATE :: filter__2D_coef            !
85   PRIVATE :: filter__2D_hann            !
86   PRIVATE :: filter__2D_hamming         !
87   PRIVATE :: filter__2D_blackman        !
88   PRIVATE :: filter__2D_gauss           !
89   PRIVATE :: filter__2D_butterworth     !
90   PRIVATE :: filter__1D_fill_value      !
91   PRIVATE :: filter__1D                 !
92   PRIVATE :: filter__1D_coef            !
93   PRIVATE :: filter__1D_hann            !
94   PRIVATE :: filter__1D_hamming         !
95   PRIVATE :: filter__1D_blackman        !
96   PRIVATE :: filter__1D_gauss           !
97   PRIVATE :: filter__1D_butterworth     !
98
99   INTERFACE filter_fill_value
100      MODULE PROCEDURE filter__fill_value_wrapper
101   END INTERFACE filter_fill_value
102
103CONTAINS
104   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
105   SUBROUTINE filter__fill_value_wrapper(td_var)
106   !-------------------------------------------------------------------
107   !> @brief
108   !> This subroutine filter variable value.
109   !>
110   !> @details
111   !> it checks if filtering method is available,
112   !>  gets parameter value, and launch filter__fill_value
113   !>
114   !> @author J.Paul
115   !> @date November, 2013 - Initial Version
116   !>
117   !> @param[inout] td_var variable structure
118   !-------------------------------------------------------------------
119
120      IMPLICIT NONE
121
122      ! Argument
123      TYPE(TVAR), INTENT(INOUT) :: td_var
124
125      ! local variable
126      CHARACTER(LEN=lc) :: cl_filter
127      CHARACTER(LEN=lc) :: cl_method
128      INTEGER(I4)       :: il_radius
129      INTEGER(I4)       :: il_nturn
130      REAL(dp)          :: dl_cutoff 
131      REAL(dp)          :: dl_alpha
132
133      TYPE(TATT)        :: tl_att
134
135      ! loop indices
136      INTEGER(I4) :: jl
137      !----------------------------------------------------------------
138
139      IF( .NOT. ASSOCIATED(td_var%d_value) )THEN
140         CALL logger_error("FILTER FILL VALUE: no array of value "//&
141         &  "associted to variable "//TRIM(td_var%c_name) )
142      ELSE
143
144         SELECT CASE(TRIM(td_var%c_filter(1)))
145
146         CASE DEFAULT
147         
148            CALL logger_trace("FILTER FILL VALUE: no filter selected "//&
149            &  "for variable "//TRIM(td_var%c_name))
150
151         CASE('hann','hamming','blackman','gauss','butterworth')
152
153            cl_method=TRIM(td_var%c_filter(1))
154
155            ! look for number of turn to be done
156            READ(td_var%c_filter(2),*) il_nturn
157            IF( il_nturn < 0 )THEN
158               CALL logger_error("FILTER FILL VALUE: invalid "//&
159               &  "number of turn ("//TRIM(td_var%c_filter(2))//")")
160            ENDIF
161
162            ! look for cut-off frequency
163            dl_cutoff=2
164            IF( TRIM(td_var%c_filter(3)) /= '' )THEN
165               READ(td_var%c_filter(3),*) dl_cutoff
166            ENDIF
167            IF( dl_cutoff < 0 )THEN
168               CALL logger_error("FILTER FILL VALUE: invalid cut-off "//&
169               &  "frequency ("//TRIM(td_var%c_filter(3))//")")
170            ENDIF
171
172            ! look for halo size
173            il_radius=1
174            IF( TRIM(td_var%c_filter(4)) /= '' )THEN
175               READ(td_var%c_filter(4),*) il_radius
176            ENDIF
177            IF( il_radius < 0 )THEN
178               CALL logger_error("FILTER FILL VALUE: invalid halo radius "//&
179               &  " ("//TRIM(td_var%c_filter(4))//")")
180            ENDIF
181
182            IF( REAL(2*il_radius+1,dp) < dl_cutoff )THEN
183               CALL logger_error("FILTER FILL VALUE: radius of halo and "//&
184               &  "spatial cut-off frequency are not suitable.")
185            ENDIF
186
187            ! look for alpha parameter
188            dl_alpha=2
189            IF( TRIM(td_var%c_filter(5)) /= '' )THEN
190               READ(td_var%c_filter(5),*) dl_alpha
191            ENDIF
192
193            SELECT CASE(TRIM(cl_method))
194            CASE('gauss','butterworth')
195               CALL logger_info("FILTER FILL VALUE: filtering "//&
196               &   " variable "//TRIM(td_var%c_name)//&
197               &   " using "//TRIM(fct_str(il_nturn))//" turn"//&
198               &   " of "//TRIM(cl_method)//" method,"//&
199               &   " with cut-off frequency of "//&
200               &        TRIM(fct_str(REAL(dl_cutoff,sp)))//&
201               &   ", halo's radius of "//&
202               &        TRIM(fct_str(il_radius))//&
203               &   ", and alpha parameter of "//&
204               &        TRIM(fct_str(REAL(dl_alpha,sp))) )
205            CASE DEFAULT
206               CALL logger_info("FILTER FILL VALUE: filtering "//&
207               &   " variable "//TRIM(td_var%c_name)//&
208               &   " using "//TRIM(fct_str(il_nturn))//" turn"//&
209               &   " of "//TRIM(cl_method)//" method,"//&
210               &   " with cut-off frequency of "//&
211               &        TRIM(fct_str(REAL(dl_cutoff,sp)))//&
212               &   " and halo's radius of "//&
213               &        TRIM(fct_str(il_radius)) )
214            END SELECT
215     
216            IF( .NOT. ANY(td_var%t_dim(1:3)%l_use) )THEN
217               ! no dimension I-J-K used
218               CALL logger_debug("FILTER FILL VALUE: no filtering can "//&
219               &  "be done for variable "//TRIM(td_var%c_name))
220            ELSE 
221
222               ! add attribute to variable
223               SELECT CASE(TRIM(cl_method))
224               CASE('gauss','butterworth')
225                  cl_filter=TRIM(fct_str(il_nturn))//'*'//TRIM(cl_method)//&
226                  &                    '('//TRIM(fct_str(REAL(dl_cutoff,sp)))//","//&
227                  &                         TRIM(fct_str(il_radius))//","//&
228                  &                         TRIM(fct_str(REAL(dl_alpha,sp)))//')'
229               CASE DEFAULT
230                  cl_filter=TRIM(fct_str(il_nturn))//'*'//TRIM(cl_method)//&
231                  &                    '('//TRIM(fct_str(REAL(dl_cutoff,sp)))//","//&
232                  &                         TRIM(fct_str(il_radius))//')'
233               END SELECT
234               tl_att=att_init('filter',cl_filter)
235               CALL var_move_att(td_var,tl_att)
236               ! clean
237               CALL att_clean(tl_att)
238
239               DO jl=1,il_nturn
240                  CALL filter__fill_value( td_var, TRIM(cl_method),  & 
241                  &                        dl_cutoff, il_radius, dl_alpha )
242               ENDDO
243            ENDIF               
244
245         END SELECT
246
247      ENDIF
248   END SUBROUTINE filter__fill_value_wrapper
249   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250   SUBROUTINE filter__fill_value(td_var, cd_name, &
251         &                       dd_cutoff, id_radius, dd_alpha)
252   !-------------------------------------------------------------------
253   !> @brief
254   !> This subroutine filtering variable value, given cut-off frequency
255   !> halo radius and alpha parameter.
256   !>
257   !> @details
258   !>    First extrabands are added to array of variable value.
259   !>    Then values are extrapolated, before apply filter.
260   !>    Finally extrabands are removed.
261   !>
262   !> @author J.Paul
263   !> @date November, 2013 - Initial Version
264   !>
265   !> @param[inout] td_var variable
266   !> @param[in] cd_name   filter name
267   !> @param[in] dd_cutoff cut-off frequency
268   !> @param[in] id_radius filter halo radius
269   !> @param[in] dd_alpha  filter parameter
270   !-------------------------------------------------------------------
271
272      IMPLICIT NONE
273
274      ! Argument
275      TYPE(TVAR)      , INTENT(INOUT) :: td_var
276      CHARACTER(LEN=*), INTENT(IN   ) :: cd_name
277      REAL(dp)        , INTENT(IN   ) :: dd_cutoff 
278      INTEGER(I4)     , INTENT(IN   ) :: id_radius
279      REAL(dp)        , INTENT(IN   ) :: dd_alpha
280
281      ! local variable
282      TYPE(TVAR)                                         :: tl_mask
283
284      INTEGER(i1)      , DIMENSION(:,:,:,:), ALLOCATABLE :: bl_mask
285
286      ! loop indices
287      INTEGER(i4) :: jl
288      !----------------------------------------------------------------
289
290      CALL logger_debug("FILTER: "//TRIM(fct_str(td_var%d_fill)) )
291
292      !1-add extraband
293      CALL extrap_add_extrabands(td_var, id_radius, id_radius)
294
295      !2-compute mask
296      ALLOCATE(bl_mask(td_var%t_dim(1)%i_len, &
297      &                td_var%t_dim(2)%i_len, &
298      &                td_var%t_dim(3)%i_len, &
299      &                td_var%t_dim(4)%i_len) )
300
301      bl_mask(:,:,:,:)=1
302      WHERE(td_var%d_value(:,:,:,:)==td_var%d_fill) bl_mask(:,:,:,:)=0     
303
304      tl_mask=var_init('tmask', bl_mask(:,:,:,:))
305
306      DEALLOCATE(bl_mask)
307
308      !3-extrapolate
309      CALL extrap_fill_value( td_var ) !, id_iext=id_radius, id_jext=id_radius )
310
311      !4-filtering
312      DO jl=1,td_var%t_dim(4)%i_len
313         IF( ALL(td_var%t_dim(1:3)%l_use) )THEN
314            ! dimension I-J-K used
315            CALL filter__3D_fill_value( td_var%d_value(:,:,:,jl),       &
316            &                           td_var%d_fill, TRIM(cd_name), &
317            &                           dd_cutoff, id_radius, dd_alpha)
318         ELSE IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 
319            ! dimension I-J used
320            CALL filter__2D_fill_value( td_var%d_value(:,:,1,jl),       &
321            &                           td_var%d_fill, TRIM(cd_name), &
322            &                           dd_cutoff, id_radius, dd_alpha)         
323         ELSE IF( td_var%t_dim(3)%l_use )THEN 
324            ! dimension K used
325            CALL filter__1D_fill_value( td_var%d_value(1,1,:,jl),       &
326            &                           td_var%d_fill, TRIM(cd_name), &
327            &                           dd_cutoff, id_radius, dd_alpha)         
328         ENDIF
329      ENDDO
330
331      !5-keep original mask
332      WHERE( tl_mask%d_value(:,:,:,:) == 0 )
333         td_var%d_value(:,:,:,:)=td_var%d_fill
334      END WHERE
335
336      ! clean
337      CALL var_clean(tl_mask)
338
339      !6-remove extraband
340      CALL extrap_del_extrabands(td_var, id_radius, id_radius)
341
342   END SUBROUTINE filter__fill_value
343   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
344   SUBROUTINE filter__3D_fill_value(dd_value, dd_fill, cd_name, &
345         &                          dd_cutoff, id_radius, dd_alpha)
346   !-------------------------------------------------------------------
347   !> @brief This subroutine compute filtered value of 3D array.
348   !>
349   !> @details
350   !>    First compute filter coefficient.
351   !>    Then apply it on each level of variable value.
352   !>
353   !> @warning array of value should have been already extrapolated before
354   !> running this subroutine.
355   !>
356   !> @author J.Paul
357   !> @date November, 2013 - Initial Version
358   !>
359   !> @param[inout] dd_value  array of value to be filtered
360   !> @param[in] dd_fill      fill value
361   !> @param[in] cd_name      filter name
362   !> @param[in] dd_cutoff    cut-off frequency
363   !> @param[in] id_radius    filter halo radius
364   !> @param[in] dd_alpha     filter parameter
365   !-------------------------------------------------------------------
366
367      IMPLICIT NONE
368
369      ! Argument     
370      REAL(dp)        , DIMENSION(:,:,:), INTENT(INOUT) :: dd_value
371      REAL(dp)        ,                   INTENT(IN   ) :: dd_fill
372      CHARACTER(LEN=*),                   INTENT(IN   ) :: cd_name
373      REAL(dp)        ,                   INTENT(IN   ) :: dd_cutoff
374      INTEGER(i4)     ,                   INTENT(IN   ) :: id_radius
375      REAL(dp)        ,                   INTENT(IN   ) :: dd_alpha     
376
377      ! local variable
378      INTEGER(i4), DIMENSION(3)                :: il_shape
379      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_coef
380
381      ! loop indices
382      INTEGER(i4) :: jk
383      !----------------------------------------------------------------
384     
385      il_shape(:)=SHAPE(dd_value(:,:,:))
386
387      ALLOCATE( dl_coef(2*id_radius+1,2*id_radius+1) )
388
389      dl_coef(:,:)=filter__2D_coef(cd_name, dd_cutoff, id_radius, dd_alpha)
390
391      DO jk=1,il_shape(3)
392         CALL filter__2D(dd_value(:,:,jk), dd_fill,dl_coef(:,:),id_radius)
393      ENDDO
394
395      DEALLOCATE( dl_coef )
396
397   END SUBROUTINE filter__3D_fill_value
398   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
399   SUBROUTINE filter__2D_fill_value(dd_value, dd_fill, cd_name, &
400         &                          dd_cutoff, id_radius, dd_alpha)
401   !-------------------------------------------------------------------
402   !> @brief This subroutine compute filtered value of 2D array.
403   !>
404   !> @details
405   !>    First compute filter coefficient.
406   !>    Then apply it on variable value.
407   !>
408   !> @warning array of value should have been already extrapolated before
409   !> running this subroutine.
410   !>
411   !> @author J.Paul
412   !> @date November, 2013 - Initial Version
413   !>
414   !> @param[inout] dd_value  array of value to be filtered
415   !> @param[in] dd_fill      fill value
416   !> @param[in] cd_name      filter name
417   !> @param[in] dd_cutoff    cut-off frequency
418   !> @param[in] id_radius    filter halo radius
419   !> @param[in] dd_alpha     filter parameter
420   !-------------------------------------------------------------------
421
422      IMPLICIT NONE
423
424      ! Argument
425      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value
426      REAL(dp)        ,                 INTENT(IN   ) :: dd_fill
427      CHARACTER(LEN=*),                 INTENT(IN   ) :: cd_name
428      REAL(dp)        ,                 INTENT(IN   ) :: dd_cutoff
429      INTEGER(i4)     ,                 INTENT(IN   ) :: id_radius
430      REAL(dp)        ,                 INTENT(IN   ) :: dd_alpha     
431
432      ! local variable
433
434      REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dl_coef
435      ! loop indices
436      !----------------------------------------------------------------
437
438      ALLOCATE( dl_coef(2*id_radius+1,2*id_radius+1) )
439
440      dl_coef(:,:)=filter__2D_coef(cd_name, dd_cutoff, id_radius, dd_alpha)
441
442      CALL filter__2D(dd_value(:,:), dd_fill, dl_coef(:,:), id_radius)
443
444      DEALLOCATE( dl_coef )
445
446   END SUBROUTINE filter__2D_fill_value
447   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
448   SUBROUTINE filter__1D_fill_value(dd_value, dd_fill, cd_name, &
449         &                          dd_cutoff, id_radius, dd_alpha)
450   !-------------------------------------------------------------------
451   !> @brief This subroutine compute filtered value of 1D array.
452   !>
453   !> @details
454   !>    First compute filter coefficient.
455   !>    Then apply it on variable value.
456   !>
457   !> @warning array of value should have been already extrapolated before
458   !> running this subroutine.
459   !>
460   !> @author J.Paul
461   !> @date November, 2013 - Initial Version
462   !>
463   !> @param[inout] dd_value  array of value to be filtered
464   !> @param[in] dd_fill      fill value
465   !> @param[in] cd_name      filter name
466   !> @param[in] dd_cutoff    cut-off frequency
467   !> @param[in] id_radius    filter halo radius
468   !> @param[in] dd_alpha     filter parameter
469   !-------------------------------------------------------------------
470
471      IMPLICIT NONE
472
473      ! Argument     
474      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value
475      REAL(dp)        ,               INTENT(IN   ) :: dd_fill
476      CHARACTER(LEN=*),               INTENT(IN   ) :: cd_name
477      REAL(dp)        ,               INTENT(IN   ) :: dd_cutoff
478      INTEGER(i4)     ,               INTENT(IN   ) :: id_radius
479      REAL(dp)        ,               INTENT(IN   ) :: dd_alpha     
480
481      ! local variable
482
483      REAL(dp), DIMENSION(:), ALLOCATABLE :: dl_coef
484      ! loop indices
485      !----------------------------------------------------------------
486
487      ALLOCATE( dl_coef(2*id_radius+1) )
488
489      dl_coef(:)=filter__1D_coef(cd_name, dd_cutoff, id_radius, dd_alpha)
490
491      CALL filter__1D(dd_value(:), dd_fill, dl_coef(:),id_radius)
492
493      DEALLOCATE( dl_coef )
494
495   END SUBROUTINE filter__1D_fill_value
496   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
497   SUBROUTINE filter__2D(dd_value, dd_fill, dd_coef, id_radius)
498   !-------------------------------------------------------------------
499   !> @brief This subroutine filtered 2D array of value
500   !>
501   !> @details
502   !>    loop on first and second dimension,
503   !>    and apply coefficient 2D array on each point
504   !>
505   !> @author J.Paul
506   !> @date November, 2013 - Initial Version
507   !>
508   !> @param[inout] dd_value  array of value to be filtered
509   !> @param[in] dd_fill      fill value
510   !> @param[in] dd_coef      filter coefficent array
511   !> @param[in] id_radius    filter halo radius
512   !-------------------------------------------------------------------
513
514      IMPLICIT NONE
515
516      ! Argument     
517      REAL(dp)        , DIMENSION(:,:), INTENT(INOUT) :: dd_value
518      REAL(dp)        ,                 INTENT(IN   ) :: dd_fill 
519      REAL(dp)        , DIMENSION(:,:), INTENT(IN   ) :: dd_coef 
520      INTEGER(i4)     ,                 INTENT(IN   ) :: id_radius
521
522      ! local variable
523      INTEGER(i4), DIMENSION(2)                :: il_shape
524      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_value
525      REAL(dp)   , DIMENSION(:,:), ALLOCATABLE :: dl_halo
526
527      ! loop indices
528      INTEGER(i4) :: jj
529      INTEGER(i4) :: ji
530      !----------------------------------------------------------------
531      il_shape(:)=SHAPE(dd_value(:,:))
532
533      ALLOCATE(dl_value(il_shape(1),il_shape(2)))
534      dl_value(:,:)=dd_value(:,:)
535
536      ALLOCATE(dl_halo(2*id_radius+1,2*id_radius+1))
537
538      DO jj=1+id_radius,il_shape(2)-id_radius
539         DO ji=1+id_radius,il_shape(1)-id_radius
540
541            dl_halo(:,:)=dd_fill
542            dl_halo(:,:)=dl_value(ji-id_radius:ji+id_radius, &
543            &                     jj-id_radius:jj+id_radius)
544
545            dd_value(ji,jj)=SUM(dl_halo(:,:)*dd_coef(:,:))
546
547         ENDDO
548      ENDDO
549
550      DEALLOCATE(dl_halo)
551      DEALLOCATE(dl_value)
552
553   END SUBROUTINE filter__2D
554   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
555   SUBROUTINE filter__1D(dd_value, dd_fill, dd_coef, id_radius)
556   !-------------------------------------------------------------------
557   !> @brief This subroutine filtered 1D array of value 
558   !>
559   !> @details
560   !>    loop on first dimension,
561   !>    and apply coefficient 1D array on each point
562   !>
563   !> @author J.Paul
564   !> @date November, 2013 - Initial Version
565   !>
566   !> @param[inout] dd_value  array of value to be filtered
567   !> @param[in] dd_fill      fill value
568   !> @param[in] dd_coef      filter coefficent array
569   !> @param[in] id_radius    filter halo radius
570   !-------------------------------------------------------------------
571
572      IMPLICIT NONE
573
574      ! Argument     
575      REAL(dp)        , DIMENSION(:), INTENT(INOUT) :: dd_value
576      REAL(dp)        ,               INTENT(IN   ) :: dd_fill 
577      REAL(dp)        , DIMENSION(:), INTENT(IN   ) :: dd_coef 
578      INTEGER(i4)     ,               INTENT(IN   ) :: id_radius
579
580      ! local variable
581      INTEGER(i4), DIMENSION(1)              :: il_shape
582      REAL(dp)   , DIMENSION(:), ALLOCATABLE :: dl_value
583
584      ! loop indices
585      INTEGER(i4) :: ji
586      !----------------------------------------------------------------
587      il_shape(:)=SHAPE(dd_value(:))
588
589      ALLOCATE(dl_value(2*id_radius+1))
590
591      DO ji=1+id_radius,il_shape(1)-id_radius
592
593         dl_value(:)=dd_fill
594         dl_value(:)=dd_value(ji-id_radius:ji+id_radius)
595
596         dd_value(ji)=SUM(dl_value(:)*dd_coef(:))
597
598      ENDDO
599
600      DEALLOCATE(dl_value)
601
602   END SUBROUTINE filter__1D
603   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
604   FUNCTION filter__2D_coef(cd_name, dd_cutoff, id_radius, dd_alpha) &
605         & RESULT (df_coef)
606   !-------------------------------------------------------------------
607   !> @brief This function compute filter coefficient.
608   !>
609   !> @details
610   !>
611   !> filter could be choose between :
612   !> - hann
613   !> - hamming
614   !> - blackman
615   !> - gauss
616   !> - butterworth
617   !> Cut-off frequency could be specify.
618   !> As well as a filter parameter for gauss and butterworth filter
619   !>
620   !> @author J.Paul
621   !> @date November, 2013 - Initial Version
622   !>
623   !> @param[in] cd_name   filter name
624   !> @param[in] dd_cutoff cut-off frequency
625   !> @param[in] id_radius filter halo radius
626   !> @param[in] dd_alpha  filter parameter
627   !> @return array of filter coefficient
628   !-------------------------------------------------------------------
629
630      IMPLICIT NONE
631
632      ! Argument     
633      CHARACTER(LEN=*), INTENT(IN) :: cd_name
634      REAL(dp)        , INTENT(IN) :: dd_cutoff
635      INTEGER(i4)     , INTENT(IN) :: id_radius
636      REAL(dp)        , INTENT(IN) :: dd_alpha
637
638      ! function
639      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
640
641      ! local variable
642
643      ! loop indices
644      !----------------------------------------------------------------
645      IF( REAL(id_radius,dp) < dd_cutoff )THEN
646         CALL logger_warn("FILTER COEF: radius of the filter halo "//&
647         &  "is lower than cut-off frequency")
648      ENDIF
649
650      SELECT CASE(TRIM(fct_lower(cd_name)))
651      CASE('hann')
652         df_coef(:,:)=filter__2D_hann(dd_cutoff, id_radius)
653      CASE('hamming')
654         df_coef(:,:)=filter__2D_hamming(dd_cutoff, id_radius)
655      CASE('blackman')
656         df_coef(:,:)=filter__2D_blackman(dd_cutoff, id_radius)
657      CASE('gauss')
658         df_coef(:,:)=filter__2D_gauss(dd_cutoff, id_radius, dd_alpha)
659      CASE('butterworth')
660         df_coef(:,:)=filter__2D_butterworth(dd_cutoff, id_radius, dd_alpha)
661      CASE DEFAULT
662         CALL logger_error("FILTER COEF: invalid filter name :"//TRIM(cd_name))
663      END SELECT
664
665   END FUNCTION filter__2D_coef
666   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
667   FUNCTION filter__1D_coef(cd_name, dd_cutoff, id_radius, dd_alpha) &
668         & RESULT (df_coef)
669   !-------------------------------------------------------------------
670   !> @brief This function compute filter coefficient.
671   !>
672   !> @details
673   !>
674   !> filter could be choose between :
675   !> - hann
676   !> - hamming
677   !> - blackman
678   !> - gauss
679   !> - butterworth
680   !> Cut-off frequency could be specify.
681   !> As well as a filter parameter for gauss an butterworth filter
682   !>
683   !> @author J.Paul
684   !> @date November, 2013 - Initial Version
685   !>
686   !> @param[in] cd_name   filter name
687   !> @param[in] dd_cutoff cut-off frequency
688   !> @param[in] id_radius filter halo radius
689   !> @param[in] dd_alpha  filter parameter
690   !> @return array of filter coefficient
691   !-------------------------------------------------------------------
692
693      IMPLICIT NONE
694
695      ! Argument     
696      CHARACTER(LEN=*), INTENT(IN) :: cd_name
697      REAL(dp)        , INTENT(IN) :: dd_cutoff
698      INTEGER(i4)     , INTENT(IN) :: id_radius
699      REAL(dp)        , INTENT(IN) :: dd_alpha
700
701      ! function
702      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
703
704      ! local variable
705
706      ! loop indices
707      !----------------------------------------------------------------
708
709      SELECT CASE(TRIM(fct_lower(cd_name)))
710      CASE('hann')
711         df_coef(:)=filter__1D_hann(dd_cutoff, id_radius)
712      CASE('hamming')
713         df_coef(:)=filter__1D_hamming(dd_cutoff, id_radius)
714      CASE('blackman')
715         df_coef(:)=filter__1D_blackman(dd_cutoff, id_radius)
716      CASE('gauss')
717         df_coef(:)=filter__1D_gauss(dd_cutoff, id_radius, dd_alpha)
718      CASE('butterworth')
719         df_coef(:)=filter__1D_butterworth(dd_cutoff, id_radius, dd_alpha)
720      CASE DEFAULT
721         CALL logger_error("FILTER COEF: invalid filter name :"//TRIM(cd_name))
722      END SELECT
723
724   END FUNCTION filter__1D_coef
725   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
726   FUNCTION filter__1D_hann(dd_cutoff, id_radius) &
727         & RESULT (df_coef)
728   !-------------------------------------------------------------------
729   !> @brief This function compute coefficient for HANN filter.
730   !>
731   !> @details
732   !>
733   !> @author J.Paul
734   !> @date November, 2013 - Initial Version
735   !>
736   !> @param[in] dd_cutoff cut-off frequency
737   !> @param[in] id_radius filter halo radius
738   !> @return array of hann filter coefficient
739   !-------------------------------------------------------------------
740
741      IMPLICIT NONE
742
743      ! Argument     
744      REAL(dp)        , INTENT(IN) :: dd_cutoff 
745      INTEGER(i4)     , INTENT(IN) :: id_radius
746
747      ! function
748      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
749
750      ! local variable
751      REAL(dp) :: dl_rad
752      REAL(dp) :: dl_sum
753
754      ! loop indices
755      INTEGER(i4) :: ji
756      !----------------------------------------------------------------
757
758      IF( dd_cutoff < 1 )THEN
759         CALL logger_error("FILTER COEF: cut-off frequency "//&
760         &  "should be greater than or equal to 1. No filter will be apply ")
761         df_coef(:)=0.
762         df_coef(id_radius+1)=1.
763      ELSE
764         DO ji=1,2*id_radius+1
765
766            dl_rad=SQRT(REAL(ji-id_radius+1,dp)**2 )
767           
768            IF( dl_rad < dd_cutoff )THEN
769               df_coef(ji)=0.5 + 0.5*COS(dp_pi*dl_rad/dd_cutoff)
770            ELSE
771               df_coef(ji)=0
772            ENDIF
773
774         ENDDO
775
776         ! normalize
777         dl_sum=SUM(df_coef(:))
778
779         df_coef(:)=df_coef(:)/dl_sum
780      ENDIF
781
782   END FUNCTION filter__1D_hann
783   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
784   FUNCTION filter__2D_hann(dd_cutoff, id_radius) &
785         & RESULT (df_coef)
786   !-------------------------------------------------------------------
787   !> @brief This function compute coefficient for HANN filter.
788   !>
789   !> @details
790   !>
791   !> @author J.Paul
792   !> @date November, 2013 - Initial Version
793   !>
794   !> @param[in] dd_cutoff cut-off frequency
795   !> @param[in] id_radius filter halo radius
796   !> @return array of hann filter coefficient
797   !-------------------------------------------------------------------
798
799      IMPLICIT NONE
800
801      ! Argument     
802      REAL(dp)   , INTENT(IN) :: dd_cutoff 
803      INTEGER(i4), INTENT(IN) :: id_radius
804
805      ! function
806      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
807
808      ! local variable
809      REAL(dp) :: dl_rad
810      REAL(dp) :: dl_sum
811
812      ! loop indices
813      INTEGER(i4) :: ji
814      INTEGER(i4) :: jj
815      !----------------------------------------------------------------
816
817      IF( dd_cutoff < 1.0_dp )THEN
818         CALL logger_error("FILTER COEF: cut-off frequency "//&
819         &  "should be greater than or equal to 1. No filter will be apply ")
820         df_coef(:,:)=0.
821         df_coef(id_radius+1,id_radius+1)=1.
822      ELSE
823         DO jj=1,2*id_radius+1
824            DO ji=1,2*id_radius+1
825
826               ! radius
827               dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 + &
828               &             REAL(jj-(id_radius+1),dp)**2 )
829               
830               IF( dl_rad < dd_cutoff )THEN
831                  df_coef(ji,jj)=0.5 + 0.5*COS(dp_pi*dl_rad/dd_cutoff)
832               ELSE
833                  df_coef(ji,jj)=0
834               ENDIF
835
836            ENDDO
837         ENDDO
838
839         ! normalize
840         dl_sum=SUM(df_coef(:,:))
841
842         df_coef(:,:)=df_coef(:,:)/dl_sum
843      ENDIF
844
845   END FUNCTION filter__2D_hann
846   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
847   FUNCTION filter__1D_hamming(dd_cutoff, id_radius) &
848         & RESULT (df_coef)
849   !-------------------------------------------------------------------
850   !> @brief This function compute coefficient for HAMMING filter.
851   !>
852   !> @details
853   !>
854   !> @author J.Paul
855   !> @date November, 2013 - Initial Version
856   !>
857   !> @param[in] dd_cutoff cut-off frequency
858   !> @param[in] id_radius filter halo radius
859   !> @return array of hamming filter coefficient
860   !-------------------------------------------------------------------
861
862      IMPLICIT NONE
863
864      ! Argument     
865      REAL(dp)        , INTENT(IN) :: dd_cutoff 
866      INTEGER(i4)     , INTENT(IN) :: id_radius
867
868      ! function
869      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
870
871      ! local variable
872      REAL(dp) :: dl_rad
873      REAL(dp) :: dl_sum
874
875      ! loop indices
876      INTEGER(i4) :: ji
877      !----------------------------------------------------------------
878
879      IF( dd_cutoff < 1 )THEN
880         CALL logger_error("FILTER COEF: cut-off frequency "//&
881         &  "should be greater than or equal to 1. No filter will be apply ")
882         df_coef(:)=0.
883         df_coef(id_radius+11)=1.
884      ELSE
885         DO ji=1,2*id_radius+1
886
887            dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 )
888         
889            IF( dl_rad < dd_cutoff )THEN
890               df_coef(ji)= 0.54 + 0.46*COS(dp_pi*dl_rad/dd_cutoff)
891            ELSE
892               df_coef(ji)=0
893            ENDIF
894
895         ENDDO
896
897         ! normalize
898         dl_sum=SUM(df_coef(:))
899
900         df_coef(:)=df_coef(:)/dl_sum
901      ENDIF
902
903   END FUNCTION filter__1D_hamming
904   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
905   FUNCTION filter__2D_hamming(dd_cutoff, id_radius) &
906         & RESULT (df_coef)
907   !-------------------------------------------------------------------
908   !> @brief This function compute coefficient for HAMMING filter.
909   !>
910   !> @details
911   !>
912   !> @author J.Paul
913   !> @date November, 2013 - Initial Version
914   !>
915   !> @param[in] dd_cutoff cut-off frequency
916   !> @param[in] id_radius filter halo radius
917   !> @return array of hamming filter coefficient
918   !-------------------------------------------------------------------
919
920      IMPLICIT NONE
921
922      ! Argument     
923      REAL(dp)        , INTENT(IN) :: dd_cutoff 
924      INTEGER(i4)     , INTENT(IN) :: id_radius
925
926      ! function
927      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
928
929      ! local variable
930      REAL(dp) :: dl_rad
931      REAL(dp) :: dl_sum
932
933      ! loop indices
934      INTEGER(i4) :: ji
935      INTEGER(i4) :: jj
936      !----------------------------------------------------------------
937
938      IF( dd_cutoff < 1 )THEN
939         CALL logger_error("FILTER COEF: cut-off frequency "//&
940         &  "should be greater than or equal to 1. No filter will be apply ")
941         df_coef(:,:)=0.
942         df_coef(id_radius+1,id_radius+1)=1.
943      ELSE
944         DO jj=1,2*id_radius+1
945            DO ji=1,2*id_radius+1
946
947               dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 + &
948                  &          REAL(jj-(id_radius+1),dp)**2 )
949           
950               IF( dl_rad < dd_cutoff )THEN
951                  df_coef(ji,jj)= 0.54 + 0.46*COS(dp_pi*dl_rad/dd_cutoff)
952               ELSE
953                  df_coef(ji,jj)=0
954               ENDIF
955
956            ENDDO
957         ENDDO
958
959         ! normalize
960         dl_sum=SUM(df_coef(:,:))
961
962         df_coef(:,:)=df_coef(:,:)/dl_sum
963      ENDIF
964
965   END FUNCTION filter__2D_hamming
966   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
967   FUNCTION filter__1D_blackman(dd_cutoff, id_radius) &
968         & RESULT (df_coef)
969   !-------------------------------------------------------------------
970   !> @brief This function compute coefficient for BLACKMAN filter.
971   !>
972   !> @details
973   !>
974   !> @author J.Paul
975   !> @date November, 2013 - Initial Version
976   !>
977   !> @param[in] dd_cutoff cut-off frequency
978   !> @param[in] id_radius filter halo radius
979   !> @return array of blackman filter coefficient
980   !-------------------------------------------------------------------
981
982      IMPLICIT NONE
983
984      ! Argument     
985      REAL(dp)        , INTENT(IN) :: dd_cutoff
986      INTEGER(i4)     , INTENT(IN) :: id_radius
987
988      ! function
989      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
990
991      ! local variable
992      REAL(dp) :: dl_rad
993      REAL(dp) :: dl_sum
994
995      ! loop indices
996      INTEGER(i4) :: ji
997      !----------------------------------------------------------------
998
999      IF( dd_cutoff < 1 )THEN
1000         CALL logger_error("FILTER COEF: cut-off frequency "//&
1001         &  "should be greater than or equal to 1. No filter will be apply ")
1002         df_coef(:)=0.
1003         df_coef(id_radius+1)=1.
1004      ELSE     
1005         DO ji=1,2*id_radius+1
1006
1007            dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 )
1008           
1009            IF( dl_rad < dd_cutoff )THEN
1010               df_coef(ji)= 0.42 + 0.5 *COS(  dp_pi*dl_rad/dd_cutoff) &
1011                  &              + 0.08*COS(2*dp_pi*dl_rad/dd_cutoff)
1012            ELSE
1013               df_coef(ji)=0
1014            ENDIF                               
1015
1016         ENDDO
1017
1018         ! normalize
1019         dl_sum=SUM(df_coef(:))
1020
1021         df_coef(:)=df_coef(:)/dl_sum
1022      ENDIF
1023
1024   END FUNCTION filter__1D_blackman
1025   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1026   FUNCTION filter__2D_blackman(dd_cutoff, id_radius) &
1027         & RESULT (df_coef)
1028   !-------------------------------------------------------------------
1029   !> @brief This function compute coefficient for BLACKMAN filter.
1030   !>
1031   !> @details
1032   !>
1033   !> @author J.Paul
1034   !> @date November, 2013 - Initial Version
1035   !>
1036   !> @param[in] dd_cutoff cut-off frequency
1037   !> @param[in] id_radius filter halo radius
1038   !> @return array of blackman filter coefficient
1039   !-------------------------------------------------------------------
1040
1041      IMPLICIT NONE
1042
1043      ! Argument     
1044      REAL(dp)        , INTENT(IN) :: dd_cutoff 
1045      INTEGER(i4)     , INTENT(IN) :: id_radius
1046
1047      ! function
1048      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
1049
1050      ! local variable
1051      REAL(dp) :: dl_rad
1052      REAL(dp) :: dl_sum
1053
1054      ! loop indices
1055      INTEGER(i4) :: ji
1056      INTEGER(i4) :: jj
1057      !----------------------------------------------------------------
1058
1059      IF( dd_cutoff < 1 )THEN
1060         CALL logger_error("FILTER COEF: cut-off frequency "//&
1061         &  "should be greater than or equal to 1. No filter will be apply ")
1062         df_coef(:,:)=0.
1063         df_coef(id_radius+1,id_radius+1)=1.
1064      ELSE     
1065         DO jj=1,2*id_radius+1
1066            DO ji=1,2*id_radius+1
1067
1068               dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 + &
1069               &             REAL(jj-(id_radius+1),dp)**2 )
1070               
1071               IF( dl_rad < dd_cutoff )THEN
1072                  df_coef(ji,jj)= 0.42 + 0.5 *COS(  dp_pi*dl_rad/dd_cutoff) &
1073                     &                 + 0.08*COS(2*dp_pi*dl_rad/dd_cutoff)
1074               ELSE
1075                  df_coef(ji,jj)=0
1076               ENDIF                               
1077
1078            ENDDO
1079         ENDDO
1080
1081         ! normalize
1082         dl_sum=SUM(df_coef(:,:))
1083
1084         df_coef(:,:)=df_coef(:,:)/dl_sum
1085      ENDIF
1086
1087   END FUNCTION filter__2D_blackman
1088   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1089   FUNCTION filter__1D_gauss(dd_cutoff, id_radius, dd_alpha) &
1090         & RESULT (df_coef)
1091   !-------------------------------------------------------------------
1092   !> @brief This function compute coefficient for GAUSS filter.
1093   !>
1094   !> @details
1095   !>
1096   !> @author J.Paul
1097   !> @date November, 2013 - Initial Version
1098   !>
1099   !> @param[in] dd_cutoff cut-off frequency
1100   !> @param[in] id_radius filter halo radius
1101   !> @param[in] dd_alpha  filter parameter
1102   !> @return array of gauss filter coefficient
1103   !-------------------------------------------------------------------
1104
1105      IMPLICIT NONE
1106
1107      ! Argument     
1108      REAL(dp)        , INTENT(IN) :: dd_cutoff 
1109      INTEGER(i4)     , INTENT(IN) :: id_radius
1110      REAL(dp)        , INTENT(IN) :: dd_alpha 
1111
1112      ! function
1113      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
1114
1115      ! local variable
1116      REAL(dp) :: dl_rad
1117      REAL(dp) :: dl_sum
1118
1119      ! loop indices
1120      INTEGER(i4) :: ji
1121      !----------------------------------------------------------------
1122
1123      IF( dd_cutoff < 1 )THEN
1124         CALL logger_error("FILTER COEF: cut-off frequency "//&
1125         &  "should be greater than or equal to 1. No filter will be apply ")
1126         df_coef(:)=0.
1127         df_coef(id_radius+1)=1.
1128      ELSE
1129         DO ji=1,2*id_radius+1
1130
1131            dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 )
1132           
1133            df_coef(ji)=EXP(-(dd_alpha*dl_rad**2)/(2*dd_cutoff**2))
1134
1135         ENDDO
1136
1137         ! normalize
1138         dl_sum=SUM(df_coef(:))
1139
1140         df_coef(:)=df_coef(:)/dl_sum
1141      ENDIF
1142
1143   END FUNCTION filter__1D_gauss
1144   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1145   FUNCTION filter__2D_gauss(dd_cutoff, id_radius, dd_alpha) &
1146         & RESULT (df_coef)
1147   !-------------------------------------------------------------------
1148   !> @brief This function compute coefficient for GAUSS filter.
1149   !>
1150   !> @details
1151   !>
1152   !> @author J.Paul
1153   !> @date November, 2013 - Initial Version
1154   !>
1155   !> @param[in] dd_cutoff cut-off frequency
1156   !> @param[in] id_radius filter halo radius
1157   !> @param[in] dd_alpha  filter parameter
1158   !> @return array of gauss filter coefficient
1159   !-------------------------------------------------------------------
1160
1161      IMPLICIT NONE
1162
1163      ! Argument     
1164      REAL(dp)        , INTENT(IN) :: dd_cutoff 
1165      INTEGER(i4)     , INTENT(IN) :: id_radius
1166      REAL(dp)        , INTENT(IN) :: dd_alpha 
1167
1168      ! function
1169      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
1170
1171      ! local variable
1172      REAL(dp) :: dl_rad
1173      REAL(dp) :: dl_sum
1174
1175      ! loop indices
1176      INTEGER(i4) :: ji
1177      INTEGER(i4) :: jj
1178      !----------------------------------------------------------------
1179
1180      IF( dd_cutoff < 1 )THEN
1181         CALL logger_error("FILTER COEF: cut-off frequency "//&
1182         &  "should be greater than or equal to 1. No filter will be apply ")
1183         df_coef(:,:)=0.
1184         df_coef(id_radius+1,id_radius+1)=1.
1185      ELSE
1186         DO jj=1,2*id_radius+1
1187            DO ji=1,2*id_radius+1
1188
1189               dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 + &
1190                  &          REAL(jj-(id_radius+1),dp)**2 )
1191               
1192               df_coef(ji,jj)=EXP(-(dd_alpha*dl_rad**2)/(2*dd_cutoff**2))
1193
1194            ENDDO
1195         ENDDO
1196
1197         ! normalize
1198         dl_sum=SUM(df_coef(:,:))
1199
1200         df_coef(:,:)=df_coef(:,:)/dl_sum
1201      ENDIF
1202
1203   END FUNCTION filter__2D_gauss
1204   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1205   FUNCTION filter__1D_butterworth(dd_cutoff, id_radius, dd_alpha) &
1206         & RESULT (df_coef)
1207   !-------------------------------------------------------------------
1208   !> @brief This function compute coefficient for BUTTERWORTH filter.
1209   !>
1210   !> @details
1211   !>
1212   !> @author J.Paul
1213   !> @date November, 2013 - Initial Version
1214   !>
1215   !> @param[in] dd_cutoff cut-off frequency
1216   !> @param[in] id_radius filter halo radius
1217   !> @param[in] dd_alpha  filter parameter
1218   !> @return array of butterworth filter coefficient
1219   !-------------------------------------------------------------------
1220
1221      IMPLICIT NONE
1222
1223      ! Argument     
1224      REAL(dp)        , INTENT(IN) :: dd_cutoff 
1225      INTEGER(i4)     , INTENT(IN) :: id_radius
1226      REAL(dp)        , INTENT(IN) :: dd_alpha 
1227
1228      ! function
1229      REAL(dp), DIMENSION(2*id_radius+1) :: df_coef
1230
1231      ! local variable
1232      REAL(dp) :: dl_rad
1233      REAL(dp) :: dl_sum
1234
1235      ! loop indices
1236      INTEGER(i4) :: ji
1237      !----------------------------------------------------------------
1238
1239      IF( dd_cutoff <= 1 )THEN
1240         CALL logger_error("FILTER COEF: cut-off frequency "//&
1241         &  "should be greater than 1. No filter will be apply ")
1242         df_coef(:)=0.
1243         df_coef(id_radius+1)=1.
1244      ELSE
1245         DO ji=1,2*id_radius+1
1246
1247            dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 )
1248           
1249            df_coef(ji)= 1 / (1+(dl_rad**2/dd_cutoff**2)**dd_alpha)
1250
1251         ENDDO
1252
1253         ! normalize
1254         dl_sum=SUM(df_coef(:))
1255
1256         df_coef(:)=df_coef(:)/dl_sum
1257      ENDIF
1258
1259   END FUNCTION filter__1D_butterworth
1260   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1261   FUNCTION filter__2D_butterworth(dd_cutoff,  id_radius, dd_alpha) &
1262         & RESULT (df_coef)
1263   !-------------------------------------------------------------------
1264   !> @brief This function compute coefficient for BUTTERWORTH filter.
1265   !>
1266   !> @details
1267   !>
1268   !> @author J.Paul
1269   !> @date November, 2013 - Initial Version
1270   !>
1271   !> @param[in] dd_cutoff cut-off frequency
1272   !> @param[in] id_radius filter halo radius
1273   !> @param[in] dd_alpha  filter parameter
1274   !> @return array of butterworth filter coefficient
1275   !-------------------------------------------------------------------
1276
1277      IMPLICIT NONE
1278
1279      ! Argument     
1280      REAL(dp)        , INTENT(IN) :: dd_cutoff 
1281      INTEGER(i4)     , INTENT(IN) :: id_radius
1282      REAL(dp)        , INTENT(IN) :: dd_alpha 
1283
1284      ! function
1285      REAL(dp), DIMENSION(2*id_radius+1,2*id_radius+1) :: df_coef
1286
1287      ! local variable
1288      REAL(dp) :: dl_rad
1289      REAL(dp) :: dl_sum
1290
1291      ! loop indices
1292      INTEGER(i4) :: ji
1293      INTEGER(i4) :: jj
1294      !----------------------------------------------------------------
1295
1296      IF( dd_cutoff <= 1 )THEN
1297         CALL logger_error("FILTER COEF: cut-off frequency "//&
1298         &  "should be greater than 1. No filter will be apply ")
1299         df_coef(:,:)=0.
1300         df_coef(id_radius+1,id_radius+1)=1.
1301      ELSE
1302         DO jj=1,2*id_radius+1
1303            DO ji=1,2*id_radius+1
1304
1305               dl_rad= SQRT( REAL(ji-(id_radius+1),dp)**2 + &
1306                  &          REAL(jj-(id_radius+1),dp)**2 )
1307               
1308               df_coef(ji,jj)= 1 / (1+(dl_rad**2/dd_cutoff**2)**dd_alpha)
1309
1310            ENDDO
1311         ENDDO
1312
1313         ! normalize
1314         dl_sum=SUM(df_coef(:,:))
1315
1316         df_coef(:,:)=df_coef(:,:)/dl_sum
1317      ENDIF
1318
1319   END FUNCTION filter__2D_butterworth
1320   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1321END MODULE filter
1322
Note: See TracBrowser for help on using the repository browser.