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    !-------------------------------------------------------------------------
29    ! [1.0]: Initialise:
30    !-------------------------------------------------------------------------  
31 
32    kz = size(rf_scale_factor)
33    
34    rf_scale_factor(:) = 0.0
35    
36    do k = 1, kz
37 
38       !-------------------------------------------------------------------------
39       ! [2.0]: Calculate RF alpha:
40       !-------------------------------------------------------------------------  
41 
42       rf_e = 0.25 * rf_passes / (rf_lengthscale(k) * rf_lengthscale(k))
43       rf_alpha(k) = 1.0 + rf_e - sqrt(rf_e * (rf_e + 2.0))
44 
45       !-------------------------------------------------------------------------
46       ! [3.0]: Calculate rescaling factor:
47       !-------------------------------------------------------------------------
48 
49       ! [3.1]: Calculate generic rescaling (normalise zero distance to 1):
50       ! For rf_passes=2 (SOAR) = 4*rf_lengthscale.
51       ! For rf_passes=infinity (Gaussian) = sqrt(8*pi)*rf_lengthscale.
52 
53       field_in(-n:n) = 0.0
54       field_in(0) = 1.0
55       field_out(-n:n) = field_in(-n:n)
56 
57       do pass = 1, rf_passes / 2
58          call da_recursive_filter_1d_adj(pass, rf_alpha(k), field_out, 2*n+1)
59       end do
60 
61       do pass = 1, rf_passes / 2
62          call da_recursive_filter_1d(pass, rf_alpha(k), field_out, 2*n+1)
63       end do
64 
65       rf_scale_factor(k) = 1.0 / field_out(0)
66 
67       ! Uncomment the following to test equivalence of UU^T and RF:
68       ! write(unit=stdout,fmt='(A,f15.5)') &
69       !    ' RF Scaling Factor = ', 1.0 / field_out(0)
70       ! field_out1(-n:n) = field_in(-n:n)
71       ! do pass = 1, rf_passes
72       !    call da_recursive_filter_1d(pass, rf_alpha(k), field_out1, 2*n+1)
73       ! end do
74 
75       ! do nn = -n, n
76       !    write(unit=stdout,fmt='(2i5,4f12.5)')k, nn, field_in(nn), &
77       !                             field_out(nn) / field_out(0), &
78       !                             field_out1(nn) / field_out1(0), &
79       !                             exp(-0.125*(real(nn)/rf_lengthscale(k))**2)
80       ! end do
81 
82    end do ! End loop over k
83    
84 end subroutine da_calculate_rf_factors
85       
86