da_calculate_rf_factors.inc

1 subroutine da_calculate_rf_factors (rf_lengthscale, rf_alpha, rf_scale_factor)
2 
3    !---------------------------------------------------------------------------
4    ! Purpose: Calculate:
5    !          1) Alpha value for recursive filter.
6    !          2) Turning conditions appropriate for B=UU^T RF.
7    !          3) Generic (depends only on rf_passes) rescaling factor.
8    !          4) Grid-dependent (hopefully temporary) scaling factor - removed.
9    !---------------------------------------------------------------------------
10 
11    implicit none
12    
13    real, intent(in)   :: rf_lengthscale(:)      ! Non-dim. R.F. lengthscale.
14    real, intent(out)  :: rf_alpha(:)            ! RF alpha factor.
15    real, intent(out)  :: rf_scale_factor(:)     ! Variance scaling factor.
16       
17    integer, parameter :: n = 500                ! 2n +1 = # pts in delta func.
18    integer            :: kz                     ! 3rd array dimension.
19    integer            :: pass                   ! Pass of recursive filter.
20    integer            :: k                      ! Loop counter
21    ! integer          :: nn                       ! Loop counter
22       
23    real               :: rf_e                   ! RF E factor.      
24 
25    real, dimension(-n:n) :: field_in, field_out  ! Working field.
26    ! real, dimension(-n:n) :: field_out1  ! Working field.
27 
28    if (trace_use_dull) call da_trace_entry("da_calculate_rf_factors")
29 
30    !-------------------------------------------------------------------------
31    ! [1.0]: Initialise:
32    !-------------------------------------------------------------------------  
33 
34    kz = size(rf_scale_factor)
35    
36    rf_scale_factor(:) = 0.0
37    
38    do k = 1, kz
39 
40       !-------------------------------------------------------------------------
41       ! [2.0]: Calculate RF alpha:
42       !-------------------------------------------------------------------------  
43 
44       rf_e = 0.25 * rf_passes / (rf_lengthscale(k) * rf_lengthscale(k))
45       rf_alpha(k) = 1.0 + rf_e - sqrt(rf_e * (rf_e + 2.0))
46 
47       !-------------------------------------------------------------------------
48       ! [3.0]: Calculate rescaling factor:
49       !-------------------------------------------------------------------------
50 
51       ! [3.1]: Calculate generic rescaling (normalise zero distance to 1):
52       ! For rf_passes=2 (SOAR) = 4*rf_lengthscale.
53       ! For rf_passes=infinity (Gaussian) = sqrt(8*pi)*rf_lengthscale.
54 
55       field_in(-n:n) = 0.0
56       field_in(0) = 1.0
57       field_out(-n:n) = field_in(-n:n)
58 
59       do pass = 1, rf_passes / 2
60          call da_recursive_filter_1d_adj(pass, rf_alpha(k), field_out, 2*n+1)
61       end do
62 
63       do pass = 1, rf_passes / 2
64          call da_recursive_filter_1d(pass, rf_alpha(k), field_out, 2*n+1)
65       end do
66 
67       rf_scale_factor(k) = 1.0 / field_out(0)
68 
69       ! Uncomment the following to test equivalence of UU^T and RF:
70       ! write(unit=stdout,fmt='(A,f15.5)') &
71       !    ' RF Scaling Factor = ', 1.0 / field_out(0)
72       ! field_out1(-n:n) = field_in(-n:n)
73       ! do pass = 1, rf_passes
74       !    call da_recursive_filter_1d(pass, rf_alpha(k), field_out1, 2*n+1)
75       ! end do
76 
77       ! do nn = -n, n
78       !    write(unit=stdout,fmt='(2i5,4f12.5)')k, nn, field_in(nn), &
79       !                             field_out(nn) / field_out(0), &
80       !                             field_out1(nn) / field_out1(0), &
81       !                             exp(-0.125*(real(nn)/rf_lengthscale(k))**2)
82       ! end do
83 
84    end do ! End loop over k
85 
86    if (trace_use_dull) call da_trace_exit("da_calculate_rf_factors")
87    
88 end subroutine da_calculate_rf_factors
89       
90