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