var/da/da_spectral/da_setlegpol.inc

   1 subroutine da_setlegpol(nj, max_wavenumber, alp_size, sinlat, coslat, alp)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    !  Method:
  10    !  Uses ECMWF recursion relation as opposed to Num Rec. one which can only go
  11    !  to m about 12 with single precision). However, still use NumRec one for
  12    !  m = 0 and 1 as ECMWF one requires m-2 for asslegpol of order m.
  13    !  Reference: Jarraud and Simmons (1990) and Belousov (1962).
  14
  15    integer, intent(in)  :: nj                         ! Latitude dimension.
  16    integer, intent(in)  :: max_wavenumber             ! Maximum zonal wavenumber.
  17    integer, intent(in)  :: alp_size                   ! Ass. Leg. Pol. array size.
  18    real,    intent(in)  :: sinlat(1:nj)               ! sine(latitude).
  19    real,    intent(in)  :: coslat(1:nj)               ! cosine(latitude).
  20    real,    intent(out) :: alp(1:alp_size)            ! Associated Legendre Polynomial.
  21
  22    integer              :: j, l, m, mm                ! Loop counters.
  23    integer              :: l1, l2, m2                 ! Markers.
  24    integer              :: index_j, index_m, index    ! Indexing counters for ALP array.
  25    integer              :: index_m2, index_l2m2       ! Indexing counters for ALP array.
  26    integer              :: index_l1m2, index_l1m      ! Indexing counters for ALP array.
  27    real                 :: s
  28    real                 :: c(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
  29    real                 :: d(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
  30    real                 :: e(2:max_wavenumber,2:max_wavenumber) ! Recursion coefficient.
  31
  32    if (trace_use) call da_trace_entry("da_setlegpol")
  33
  34    alp(:) = 0.0
  35
  36    ! Compute Associated Legendre polynomials for latitude range:
  37
  38    do j = 1, (nj + 1) / 2
  39       index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
  40
  41       ! use Num. Rec. recursion relation for m = 0, 1:
  42
  43       do m = 0, 1
  44          index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
  45          do l = m, max_wavenumber
  46             index = index_m + index_j + l
  47             call da_asslegpol(l, m, sinlat(j), coslat(j), alp(index))
  48
  49             ! Multiply by normalization constant
  50             ! (to ensure 1/2 integral^1_-1 Pml Pml1 = 1 if l = l1):
  51
  52             s = 1.0
  53             do mm = l-m+1, l+m
  54                s = s * real(mm)
  55             end do
  56             alp(index) = sqrt(real(2*l+1) / s) * alp(index)
  57          end do
  58       end do
  59    end do
  60
  61    ! Jarraud recursion relation coefficients:
  62
  63    do m = 2, max_wavenumber
  64       do l = m, max_wavenumber
  65          c(l,m) = sqrt ((real(2*l+1)/real(2*l-3)) * (real(l+m-1)/real(l+m)) * &
  66                   (real(l+m-3)/real(l+m-2)))
  67          d(l,m) = sqrt ((real(2*l+1)/real(2*l-1)) * (real(l+m-1)/real(l+m)) * &
  68                   (real(l-m+1)/real(l+m-2)))
  69          e(l,m) = sqrt ((real(2*l+1)/real(2*l-1)) * (real(l-m)  /real(l+m)))
  70       end do
  71    end do
  72
  73    ! use Jarraud recursion relation for m>=2:
  74
  75    do j = 1, (nj + 1) / 2
  76       index_j = (j - 1) * (max_wavenumber + 1) * (max_wavenumber + 2) / 2
  77
  78       do m = 2, max_wavenumber
  79          index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
  80          m2 = m - 2
  81          index_m2 = m2 * (max_wavenumber + 1 - m2) + m2 * (m2+1) / 2 + 1 - m2
  82
  83          do l = m, max_wavenumber
  84             l1 = l - 1
  85             l2 = l - 2
  86             index = index_j + index_m + l
  87             index_l2m2 = index_j + index_m2 + l2
  88             index_l1m2 = index_j + index_m2 + l1
  89             index_l1m  = index_j + index_m  + l1
  90
  91             alp(index) = c(l,m) * alp(index_l2m2) - d(l,m) *  sinlat(j) * &
  92                alp(index_l1m2) + e(l,m) * sinlat(j) * alp(index_l1m)
  93          end do
  94       end do
  95    end do
  96
  97    if (trace_use) call da_trace_exit("da_setlegpol")
  98
  99 end subroutine da_setlegpol
 100
 101