C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPHERP     SPHERICAL TRANSFORM TO A POINT
C   PRGMMR: KANAMITSU     ORG: NP23          DATE: 1992-04-16
C
C ABSTRACT: TRANSFORMS A SPECTRAL FIELD TO A POINT ON THE GLOBE.
C
C PROGRAM HISTORY LOG:
C 1992-04-16  IREDELL
C
C USAGE:    CALL SPHERP(WAVE,MAXWV,IROMB,IDY,MLT,FAC,RLAT,RLON,VALUE)
C
C   INPUT ARGUMENT LIST:
C     WAVE     - IF IDIR < 0, COMPLEX (KMAX) FIELD TO TRANSFORM,
C                WHERE KMAX=(MAXWV+1)*(IROMB+1)*MAXWV+2)/2.
C                WAVE STARTS AT THE GLOBAL MEAN COMPONENT AND THEN
C                CONTAINS THE ZONALLY SYMMETRIC COMPONENTS.
C     MAXWV    - SPECTRAL TRUNCATION OF THE WAVE
C     IROMB    - IROMB = 0 FOR TRIANGULAR TRUNCATION
C                IROMB = 1 FOR RHOMBOIDAL TRUNCATION
C     IDIR     - 0 FOR NORMAL, 1 FOR LEGENDRE DERIVATIVE
C     MLT      - MULTIPLICATION OPTION ON WAVE VALUES.
C                MLT = 0 FOR NO MULTIPLICATION,
C                MLT = 1 TO MULTIPLY WAVE BY FAC
C                MLT = -1 TO MULTIPLY WAVE BY SQRT(-1)*FAC
C     FAC      - IF MLT.NE.0, REAL (KMAX) FIELD TO MULTIPLY WAVE.
C     RLAT     - LATITUDE IN DEGREES NORTH
C     RLON     - LONGITUDE IN DEGREES EAST
C
C   OUTPUT ARGUMENT LIST:
C     VALUE    - TRANSFORMED VALUE
C
C   SUBPROGRAMS CALLED:
C     (LEGENDR2) - COMPUTE LEGENDRE POLYNOMIALS
C
C   REMARKS: FORTRAN 9X EXTENSIONS ARE USED.
C     FPP CAN BE USED TO INLINE SUBPROGRAMS.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM-SP
C
C$$$
CFPP$ EXPAND(LEGENDR2)
      SUBROUTINE SPHERP(WAVE,MAXWV,IROMB,IDY,MLT,FAC,RLAT,RLON,VALUE)
      REAL WAVE((MAXWV+1)*((IROMB+1)*MAXWV+2))
      REAL FAC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL PNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL DNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL FR(0:MAXWV),FI(0:MAXWV)
      PARAMETER(PI=3.1415926535898,PI180=PI/180.)
      SINLAT=SIN(RLAT*PI180)
      CALL LEGENDR2(SINLAT,MAXWV,IROMB,PNM,DNM)
      IF(IDY.NE.0) THEN
        PNM=DNM
      ENDIF
      IF(MLT.NE.0) THEN
        PNM=PNM*FAC
      ENDIF
      K=0
      DO M=0,MAXWV
        FR(M)=0.
        FI(M)=0.
      ENDDO
      DO M=0,MAXWV
        DO N=M,IROMB*M+MAXWV
          FR(M)=FR(M)+WAVE(2*K+1)*PNM(K+1)
          FI(M)=FI(M)+WAVE(2*K+2)*PNM(K+1)
          K=K+1
        ENDDO
      ENDDO
      VALUE=0.
      RLAM=RLON*PI180
      IF(MLT.GE.0) THEN
        DO M=1,MAXWV
          VALUE=VALUE+FR(M)*COS(RLAM*M)-FI(M)*SIN(RLAM*M)
        ENDDO
      ELSE
        DO M=1,MAXWV
          VALUE=VALUE-FI(M)*COS(RLAM*M)-FR(M)*SIN(RLAM*M)
        ENDDO
      ENDIF
      VALUE=FR(0)+2.*VALUE
      RETURN
      END
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  LEGENDR2   COMPUTE LEGENDRE POLYNOMIALS.
C   PRGMMR: IREDELL       ORG: NP23          DATE: 1992-04-16
C
C ABSTRACT: EVALUATES THE NORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
C   AT A GIVEN COSINE OF COLATITUDE.
C
C PROGRAM HISTORY LOG:
C 1992-04-16  IREDELL
C
C USAGE:    CALL LEGENDR2(COSCLT,MAXWV,IROMB,PNM,DNM)
C
C   INPUT ARGUMENT LIST:
C     COSCLT   - COSINE OF COLATITUDE AT WHICH TO EVALUTATE.
C     MAXWV    - SPECTRAL TRUNCATION
C     IROMB    - IROMB = 0 FOR TRIANGULAR TRUNCATION
C                IROMB = 1 FOR RHOMBOIDAL TRUNCATION
C
C   OUTPUT ARGUMENT LIST:
C     PNM      - REAL (KMAX) POLYNOMIAL VALUES,
C     DNM      - REAL (KMAX) POLYNOMIAL DERIVATIVES,
C                WHERE KMAX=(MAXWV+1)*(IROMB+1)*MAXWV+2)/2.
C
C   REMARKS: FORTRAN 9X EXTENSIONS ARE USED.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM-SP
C
C$$$
      SUBROUTINE LEGENDR2(COSCLT,MAXWV,IROMB,PNM,DNM)
      DIMENSION PNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      DIMENSION DNM((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      DIMENSION EPS(0:100)
      SINCLT=SQRT(1.-COSCLT**2)
      NM=0
      K=1
      PNM(K)=SQRT(0.5)
CDIR$ NEXTSCALAR
      DO M=1,MAXWV-NM*(1-IROMB)
        PNMK=PNM(K)
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        PNM(K)=PNMK*SINCLT*SQRT(FLOAT(2*M+1)/FLOAT(2*M))
      ENDDO
      NM=1
CDIR$ IVDEP
      DO M=0,MAXWV-NM*(1-IROMB)
        N=M+NM
        EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        PNM(K)=COSCLT*PNM(K-1)/EPSNM
        DNM(K-1)=EPSNM*PNM(K)*(N-1)
        EPS(M)=EPSNM
      ENDDO
      DO NM=2,MAXWV
CDIR$ IVDEP
        DO M=0,MAXWV-NM*(1-IROMB)
          N=M+NM
          EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
          K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
          PNM(K)=(COSCLT*PNM(K-1)-EPS(M)*PNM(K-2))/EPSNM
          DNM(K-1)=EPSNM*PNM(K)*(N-1)-EPS(M)*PNM(K-2)*N
          EPS(M)=EPSNM
        ENDDO
      ENDDO
CDIR$ IVDEP
      DO M=0,MAXWV
        NM=MAXWV-M*(1-IROMB)
        N=M+NM
        EPS(M)=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        N=N+1
        EPSNM=SQRT(FLOAT(N**2-M**2)/FLOAT(4*N**2-1))
        K=M*(2*MAXWV+2+(1-M)*(1-IROMB))/2+(NM+1)
        K=K+1
        PNMK=(COSCLT*PNM(K-1)-EPS(M)*PNM(K-2))/EPSNM
        DNM(K-1)=EPSNM*PNMK*(N-1)-EPS(M)*PNM(K-2)*N
      ENDDO
      RETURN
      END
