***********************************************************************
*     CALCULATE LEGENDRE FUNCTIONS OF ZONAL WAVENUMBER M AT A LATITUDE
***********************************************************************
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE SHPPMJ(MM,JM,M,J,PMJ,YS,PY,R)

      REAL PMJ(0:MM+1)
      REAL YS(0:JM),PY(2,0:JM,0:MM),R((MM+1)*(MM+1))

      LA=M*(2*MM-M)+MM
      LB=LA-MM+M-1

      PMJ(M  )=PY(1,J,M)
      PMJ(M+1)=PY(2,J,M)

      DO 10 N=M+2,MM+1
        PMJ(N)=R(LA+N)*(YS(J)*PMJ(N-1)-R(LB+N)*PMJ(N-2))
   10 CONTINUE

      END
