New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
filter.f90 in branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/TOOLS/SIREN/src – NEMO

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/TOOLS/SIREN/src/filter.f90 @ 5951

Last change on this file since 5951 was 5951, checked in by timgraham, 8 years ago

Merged trunk r5936 into branch

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