C-----------------------------------------------------------------------
      SUBROUTINE SPTGPT(IROMB,MAXWV,KMAX,NMAX,
     &                  KWSKIP,KGSKIP,NRSKIP,NGSKIP,
     &                  RLAT,RLON,WAVE,GP)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPTGPT     TRANSFORM SPECTRAL SCALAR TO STATION POINTS
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-02-29
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM
C           FROM SPECTRAL COEFFICIENTS OF SCALAR QUANTITIES
C           TO SPECIFIED SETS OF STATION POINTS ON THE GLOBE.
C           THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL.
C           THE WAVE AND POINT FIELDS MAY HAVE GENERAL INDEXING,
C           BUT EACH WAVE FIELD IS IN SEQUENTIAL 'IBM ORDER',
C           I.E. WITH ZONAL WAVENUMBER AS THE SLOWER INDEX.
C           THE TRANSFORMS ARE ALL MULTIPROCESSED OVER STATIONS.
C           TRANSFORM SEVERAL FIELDS AT A TIME TO IMPROVE VECTORIZATION.
C           SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C PROGRAM HISTORY LOG:
C   96-02-29  IREDELL
C   98-06-30  EBISUZAKI    LINUX PORT
C
C USAGE:    CALL SPTGPT(IROMB,MAXWV,KMAX,NMAX,
C    &                  KWSKIP,KGSKIP,NRSKIP,NGSKIP,
C    &                  RLAT,RLON,WAVE,GP)
C   INPUT ARGUMENTS:
C     IROMB    - INTEGER SPECTRAL DOMAIN SHAPE
C                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C     MAXWV    - INTEGER SPECTRAL TRUNCATION
C     KMAX     - INTEGER NUMBER OF FIELDS TO TRANSFORM.
C     NMAX     - INTEGER NUMBER OF STATION POINTS TO RETURN
C     KWSKIP   - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
C                (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0)
C     KGSKIP   - INTEGER SKIP NUMBER BETWEEN STATION POINT SETS
C                (DEFAULTS TO NMAX IF KGSKIP=0)
C     NRSKIP   - INTEGER SKIP NUMBER BETWEEN STATION LATS AND LONS
C                (DEFAULTS TO 1 IF NRSKIP=0)
C     NGSKIP   - INTEGER SKIP NUMBER BETWEEN STATION POINTS
C                (DEFAULTS TO 1 IF NGSKIP=0)
C     RLAT     - REAL (*) STATION LATITUDES IN DEGREES
C     RLON     - REAL (*) STATION LONGITUDES IN DEGREES
C     WAVE     - REAL (*) WAVE FIELDS
C   OUTPUT ARGUMENTS:
C     GP       - REAL (*) STATION POINT SETS
C
C SUBPROGRAMS CALLED:
C   SPWGET       GET WAVE-SPACE CONSTANTS
C   SPLEGEND     COMPUTE LEGENDRE POLYNOMIALS
C   SPSYNTH      SYNTHESIZE FOURIER FROM SPECTRAL
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
      REAL RLAT(*),RLON(*),WAVE(*),GP(*)
      REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
      REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
      REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
      INTEGER MP(KMAX)
      REAL SLON(MAXWV),CLON(MAXWV)
      REAL WTOP(2*(MAXWV+1),KMAX)
      REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
      REAL F(2*MAXWV+3,2,KMAX)
      PARAMETER(PI=3.14159265358979)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE PRELIMINARY CONSTANTS
      CALL SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
      MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
      MXTOP=MAXWV+1
      IDIM=2*MAXWV+3
      KW=KWSKIP
      KG=KGSKIP
      NR=NRSKIP
      NG=NGSKIP
      IF(KW.EQ.0) KW=2*MX
      IF(KG.EQ.0) KG=NMAX
      IF(NR.EQ.0) NR=1
      IF(NG.EQ.0) NG=1
c     MP=0
      do ii = 1, kmax
          MP(ii)=0
      end do
CMIC$ DO ALL PRIVATE(K)
CMIC$&       SHARED(KMAX,MXTOP,WTOP)
      DO K=1,KMAX
c       WTOP(1:2*MXTOP,K)=0
	do ii = 1, 2*mxtop
            WTOP(ii,K)=0
	enddo
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE STATION FIELDS
CMIC$ DO ALL PRIVATE(N,K,L,RADLAT,RADLON,SLAT1,CLAT1)
CMIC$&       PRIVATE(PLN,PLNTOP,F,SLON,CLON,NK)
CMIC$&       SHARED(NMAX,MX,MXTOP,IDIM,KMAX,PI,RLAT,RLON,NR)
CMIC$&       SHARED(IROMB,MAXWV,EPS,EPSTOP,KW,MP,WAVE,WTOP,GP,NG,KG)
      DO N=1,NMAX
        RADLAT=PI/180*RLAT((N-1)*NR+1)
        RADLON=PI/180*RLON((N-1)*NR+1)
        IF(RLAT((N-1)*NR+1).GE.89.9995) THEN
          SLAT1=1.
          CLAT1=0.
        ELSEIF(RLAT((N-1)*NR+1).LE.-89.9995) THEN
          SLAT1=-1.
          CLAT1=0.
        ELSE
          SLAT1=SIN(RADLAT)
          CLAT1=COS(RADLAT)
        ENDIF
        CALL SPLEGEND(IROMB,MAXWV,SLAT1,CLAT1,EPS,EPSTOP,
     &                PLN,PLNTOP)
        CALL SPSYNTH(IROMB,MAXWV,2*MAXWV,IDIM,KW,2*MXTOP,KMAX,
     &               CLAT1,PLN,PLNTOP,MP,WAVE,WTOP,F)
        DO L=1,MAXWV
          SLON(L)=SIN(L*RADLON)
          CLON(L)=COS(L*RADLON)
        ENDDO
        DO K=1,KMAX
          NK=(N-1)*NG+(K-1)*KG+1
          GP(NK)=F(1,1,K)
        ENDDO
        IF(KMAX.EQ.1) THEN
          NK=(N-1)*NG+1
          DO L=1,MAXWV
            GP(NK)=GP(NK)+2.*(F(2*L+1,1,1)*CLON(L)
     &                       -F(2*L+2,1,1)*SLON(L))
          ENDDO
        ELSE
          DO L=1,MAXWV
            DO K=1,KMAX
              NK=(N-1)*NG+(K-1)*KG+1
              GP(NK)=GP(NK)+2.*(F(2*L+1,1,K)*CLON(L)
     &                         -F(2*L+2,1,K)*SLON(L))
            ENDDO
          ENDDO
        ENDIF
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END
