var/da/da_spectral/da_calc_power_spectrum.inc

   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*8, 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             :: i,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    character (len=filename_len) :: filename
  30
  31    do i=1,num_alpha_corr_types
  32       call da_get_unit(alpha_corr_unit1(i))
  33       call da_get_unit(alpha_corr_unit2(i))
  34       write (unit=filename,fmt='(A,I1)') "alpha_corr1_",i
  35       open(unit=alpha_corr_unit1(i),file=filename,status="replace")
  36       write (unit=filename,fmt='(A,I1)') "alpha_corr2_",i
  37       open(unit=alpha_corr_unit2(i),file=filename,status="replace")
  38    end do
  39
  40    !-----------------------------------------------------------------------------
  41    ! [1] Switch lats from -pi/2 to pi/2 to 0 to pi:
  42    !-----------------------------------------------------------------------------
  43
  44    max_wave_nj = nj / 2 - 1
  45    allocate(alp(1:nj,0:max_wave_nj))
  46    allocate(power_nj(0:max_wave_nj))
  47
  48    call da_get_gausslats(nj, glats, gwgts, sinlat, coslat)
  49    glats = glats + 0.5 * pi
  50
  51    ! Get m=0 Associated Legendre Polynomials:
  52
  53    do l = 0, max_wave_nj
  54       odd = .false.
  55       if (real(int(0.5 * real(l))) /= 0.5 * l) odd = .true.
  56
  57       do j = 1, nj
  58          call da_asslegpol(l, 0, sinlat(j), coslat(j), alp(j,l))
  59          ! Reverse order of alps to account for latitude/angle difference:
  60          if (odd) alp(j,l) = -alp(j,l)
  61       end do
  62    end do
  63
  64    !-----------------------------------------------------------------------------
  65    ! [2] Define correlation function:
  66    !-----------------------------------------------------------------------------
  67
  68    corr_scale = alpha_corr_scale / earth_radius
  69    corr_scale_inv = 1.0 / corr_scale
  70
  71    do j = 1, nj
  72       ! d(j) = 0.5 * glats(j) * corr_scale_inv
  73       d(j) = glats(j) * corr_scale_inv
  74
  75       if (alpha_corr_type == alpha_corr_type_exp) then ! Exponential.
  76          corr(j) = exp(-d(j))
  77       else if (alpha_corr_type == alpha_corr_type_soar) then ! SOAR
  78          d(j) = 2.0 * d(j)
  79          corr(j) = (1.0 + d(j)) * exp(-d(j))
  80       else if (alpha_corr_type == alpha_corr_type_gaussian) then ! Gaussian
  81          corr(j) = exp(-d(j) * d(j))
  82       end if
  83    end do
  84
  85    do j = 1, nj
  86       write(unit=alpha_corr_unit1(alpha_corr_type),fmt='(i4,2f12.4)') &
  87         j, earth_radius * glats(j), corr(j)
  88    end do
  89
  90    !--------------------------------------------------------------------------
  91    ! [3] Calculate power spectra:
  92    !--------------------------------------------------------------------------
  93
  94    ! Calculate power spectrum (and truncate if has -ve values).
  95    ! Power spectrum at this stage is is the Dl=sqrt(2l+1)*Bl of Boer(1983).
  96    ! This ensures the total variance = sum(Dl).
  97
  98    power_nj(:) = 0.0
  99    do l = 0, max_wave_nj
 100       power_nj(l) = 0.5 * sqrt(2.0 * real(l) + 1.0) * &
 101                     sum(gwgts(1:nj) * corr(1:nj) * alp(1:nj,l))
 102
 103       if (power_nj(l) < 0.0) power_nj(l) = 0.0
 104    end do
 105    write(unit=stdout,fmt='(a,2f12.5)')' Total, unscale variance = ', &
 106                 sum(power_nj(0:max_wave_nj))
 107
 108    ! Rescale so variance = 1 (take out later?):
 109    variance = sum(power_nj(0:max_wave_nj))
 110    power_nj(0:max_wave_nj) = power_nj(0:max_wave_nj) / variance
 111
 112    do l = 0, max_wave_nj
 113       write(unit=alpha_corr_unit2(alpha_corr_type),fmt='(i4,2f12.4)') &
 114         l, power_nj(l), sum(power_nj(0:l))
 115    end do
 116
 117    write(unit=stdout,fmt='(a,2i6)')' Total, truncated wave_max = ', &
 118                      max_wave_nj, max_wave
 119    write(unit=stdout,fmt='(a,2f12.5)')' Total, truncated variance = ', &
 120                 sum(power_nj(0:max_wave_nj)), sum(power_nj(0:max_wave))
 121
 122    power(0:max_wave) = power_nj(0:max_wave)
 123
 124    ! Add compactly supported correlation from calc_globalspectral later?
 125
 126    deallocate(alp)
 127    deallocate(power_nj)
 128
 129    do i=1,num_alpha_corr_types
 130       close (alpha_corr_unit1(i))
 131       close (alpha_corr_unit2(i))
 132       call da_free_unit (alpha_corr_unit1(i))
 133       call da_free_unit (alpha_corr_unit2(i))
 134    end do
 135
 136 end subroutine da_calc_power_spectrum
 137
 138