da_calc_power_spectrum.inc
References to this file elsewhere.
1 subroutine da_calc_power_spectrum(max_wave, power)
2
3 !-----------------------------------------------------------------------
4 ! Purpose: Calculate power spectrum
5 !-----------------------------------------------------------------------
6
7 implicit none
8
9 integer, parameter :: nj = 200 ! #Gaussian lats (even).
10
11 integer, intent(in) :: max_wave ! Smallest wavenumber for grid.
12 real, intent(out) :: power(0:max_wave) ! Power spectrum.
13
14 real :: glats(1:nj) ! Gaussian latitudes.
15 real :: gwgts(1:nj) ! Gaussian weights.
16 real :: sinlat(1:nj) ! sine.
17 real :: coslat(1:nj) ! cosine.
18
19 integer :: max_wave_nj ! Maximum wavenumber.
20 integer :: j, l ! Loop counters.
21 real :: corr_scale ! Correlation scale.
22 real :: corr_scale_inv ! 1/corr_scale
23 real :: variance ! Variance = sum(power).
24 real :: d(1:nj) ! Temp array.
25 real :: corr(1:nj) ! Correlation function.
26 logical :: odd ! true if odd.
27 real, allocatable :: alp(:,:) ! Associated Legendre Polynomials.
28 real, allocatable :: power_nj(:) ! Power spectrum.
29
30 !-----------------------------------------------------------------------------
31 ! [1] Switch lats from -pi/2 to pi/2 to 0 to pi:
32 !-----------------------------------------------------------------------------
33
34 max_wave_nj = nj / 2 - 1
35 allocate(alp(1:nj,0:max_wave_nj))
36 allocate(power_nj(0:max_wave_nj))
37
38 call da_get_gausslats(nj, glats, gwgts, sinlat, coslat)
39 glats = glats + 0.5 * pi
40
41 ! Get m=0 Associated Legendre Polynomials:
42
43 do l = 0, max_wave_nj
44 odd = .false.
45 if (real(int(0.5 * real(l))) /= 0.5 * l) odd = .true.
46
47 do j = 1, nj
48 call da_asslegpol(l, 0, sinlat(j), coslat(j), alp(j,l))
49 ! Reverse order of alps to account for latitude/angle difference:
50 if (odd) alp(j,l) = -alp(j,l)
51 end do
52 end do
53
54 !-----------------------------------------------------------------------------
55 ! [2] Define correlation function:
56 !-----------------------------------------------------------------------------
57
58 corr_scale = alpha_corr_scale / earth_radius
59 corr_scale_inv = 1.0 / corr_scale
60
61 do j = 1, nj
62 ! d(j) = 0.5 * glats(j) * corr_scale_inv
63 d(j) = glats(j) * corr_scale_inv
64
65 if (alpha_corr_type == 1) then ! Exponential.
66 corr(j) = exp(-d(j))
67 else if (alpha_corr_type == 2) then ! SOAR
68 d(j) = 2.0 * d(j)
69 corr(j) = (1.0 + d(j)) * exp(-d(j))
70 else if (alpha_corr_type == 3) then ! Gaussian
71 corr(j) = exp(-d(j) * d(j))
72 end if
73 end do
74
75 do j = 1, nj
76 write(unit=alpha_corr_unit1(alpha_corr_type),fmt='(i4,2f12.4)') &
77 j, 6371.0 * glats(j), corr(j)
78 end do
79
80 !--------------------------------------------------------------------------
81 ! [3] Calculate power spectra:
82 !--------------------------------------------------------------------------
83
84 ! Calculate power spectrum (and truncate if has -ve values).
85 ! Power spectrum at this stage is is the Dl=sqrt(2l+1)*Bl of Boer(1983).
86 ! This ensures the total variance = sum(Dl).
87
88 power_nj(:) = 0.0
89 do l = 0, max_wave_nj
90 power_nj(l) = 0.5 * sqrt(2.0 * real(l) + 1.0) * &
91 sum(gwgts(1:nj) * corr(1:nj) * alp(1:nj,l))
92
93 if (power_nj(l) < 0.0) power_nj(l) = 0.0
94 end do
95 write(unit=stdout,fmt='(a,2f12.5)')' Total, unscale variance = ', &
96 sum(power_nj(0:max_wave_nj))
97
98 ! Rescale so variance = 1 (take out later?):
99 variance = sum(power_nj(0:max_wave_nj))
100 power_nj(0:max_wave_nj) = power_nj(0:max_wave_nj) / variance
101
102 do l = 0, max_wave_nj
103 write(unit=alpha_corr_unit2(alpha_corr_type),fmt='(i4,2f12.4)') &
104 l, power_nj(l), sum(power_nj(0:l))
105 end do
106
107 write(unit=stdout,fmt='(a,2i6)')' Total, truncated wave_max = ', &
108 max_wave_nj, max_wave
109 write(unit=stdout,fmt='(a,2f12.5)')' Total, truncated variance = ', &
110 sum(power_nj(0:max_wave_nj)), sum(power_nj(0:max_wave))
111
112 power(0:max_wave) = power_nj(0:max_wave)
113
114 ! Add compactly supported correlation from calc_globalspectral later?
115
116 deallocate(alp)
117 deallocate(power_nj)
118
119 end subroutine da_calc_power_spectrum
120
121