legtri.F90

Go to the documentation of this file.

SUBROUTINE legtri (PSIN,KCP,KDIM,PALP)

!**** *LEGTRI* - *LEGENDRE FUNCTIONS FOR A TRIANGULAR TRUNCATION.

!     J.F.GELEYN     E.C.M.W.F.     03/06/82.

!     PURPOSE.
!     --------

!          THIS ROUTINE COMPUTES THE VALUES *PALP* FOR THE ARGUMENT
!     *PSIN* OF THE NORMALISED *LEGENDRE ASSOCIATED FUNCTIONS IN THE
!     ORDER ((JN1=JM1,KCP),JM1=1,KCP) FOR JN=JN1-1 AND JM=JM1-1 .

!**   INTERFACE.
!     ----------

!          *LEGTRI* IS CALLED FROM *RADMOD*.
!          THERE ARE THREE DUMMY ARGUMENTS: *PSIN* IS THE SINE OF
!     LATITUDE.
!                                           *KCP* IS ONE PLUS THE LIMIT
!     WAVE NUMBER.
!                                           *PALP* IS THE ARRAY OF THE
!     RESULTS.

!     METHOD.
!     -------

!          SIMPLE RECURENCE FORMULA.

!     EXTERNALS.
!     ----------

!          NONE.

!     REFERENCE.
!     ----------

!          NONE.


#include "tsmbkind.h"

IMPLICIT NONE


!     DUMMY INTEGER SCALARS
integer_m :: kcp
integer_m :: kdim

!     DUMMY REAL SCALARS
real_b :: psin

real_b :: palp(kdim)

!     LOCAL INTEGER SCALARS
integer_m :: ic, icp, ii, im, im2, jm1, jn

!     LOCAL REAL SCALARS
real_b :: z2m, zcos, ze1, ze2, zf1m, zf2m, zm, zn, zn2, zre1, zsin


!     ------------------------------------------------------------------

!*         1.     PRELIMINARY SETTING.
!                 ----------- --------


zsin=psin
icp=kcp

!     ------------------------------------------------------------------

!*         2.     COMPUTATIONS.
!                 -------------


ic=icp-1
zcos=sqrt(_one_-zsin**2)
ii=2
palp(1)=_one_
zf1m=sqrt(3._jprb)
palp(2)=zf1m*zsin
DO jm1=1,icp
  im=jm1-1
  zm=im
  z2m=zm+zm
  zre1=sqrt(z2m+3._jprb)
  ze1=_one_/zre1
  IF(im == 0) GO TO 201
  zf2m=zf1m*zcos/sqrt(z2m)
  zf1m=zf2m*zre1
  ii=ii+1
  palp(ii)=zf2m
  IF(im == ic) GO TO 203
  ii=ii+1
  palp(ii)=zf1m*zsin
  IF(jm1 == ic) GO TO 203
  201 CONTINUE
  im2=im+2
  DO jn=im2,ic
    zn=jn
    zn2=zn**2
    ze2=sqrt((4._jprb*zn2-_one_)/(zn2-zm**2))
    ii=ii+1
    palp(ii)=ze2*(zsin*palp(ii-1)-ze1*palp(ii-2))
    ze1=_one_/ze2
  ENDDO
  203 continue
ENDDO

!     ------------------------------------------------------------------

!*         3.     RETURN.
!                 -------


RETURN
END SUBROUTINE legtri