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