C-----------------------------------------------------------------------
      SUBROUTINE SPTGPTD(IROMB,MAXWV,KMAX,NMAX,
     &                   KWSKIP,KGSKIP,NRSKIP,NGSKIP,
     &                   RLAT,RLON,WAVE,XP,YP)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPTGPTD    TRANSFORM SPECTRAL TO STATION POINT GRADIENTS
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 FIELDS
C           TO SPECIFIED SETS OF STATION POINT GRADIENTS 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 SPTGPTD(IROMB,MAXWV,KMAX,NMAX,
C    &                   KWSKIP,KGSKIP,NRSKIP,NGSKIP,
C    &                   RLAT,RLON,WAVE,XP,YP)
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     XP       - REAL (*) STATION POINT X-GRADIENT SETS
C     YP       - REAL (*) STATION POINT Y-GRADIENT SETS
C
C SUBPROGRAMS CALLED:
C   SPWGET       GET WAVE-SPACE CONSTANTS
C   SPLAPLAC     COMPUTE LAPLACIAN IN SPECTRAL SPACE
C   SPTGPTV      TRANSFORM SPECTRAL VECTOR TO STATION POINTS
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
      REAL RLAT(*),RLON(*),WAVE(*),XP(*),YP(*)
      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)
      REAL WD((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,KMAX)
      REAL WZ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,KMAX)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE PRELIMINARY CONSTANTS
      CALL SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)
      MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2
      MDIM=2*MX+1
      KW=KWSKIP
      IF(KW.EQ.0) KW=2*MX
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE STATION FIELDS
CMIC$ DO ALL PRIVATE(K,KWS)
CMIC$&           SHARED(KMAX,IROMB,MAXWV,ENN1,KW,WAVE,WD,WZ,MX)
      DO K=1,KMAX
        KWS=(K-1)*KW
        CALL SPLAPLAC(IROMB,MAXWV,ENN1,WAVE(KWS+1),WD(1,K),1)
c       WZ(1:2*MX,K)=0.
        do ii = 1, 2*MX
            WZ(ii,K)=0.0
        enddo
      ENDDO
      CALL SPTGPTV(IROMB,MAXWV,KMAX,NMAX,MDIM,KGSKIP,NRSKIP,NGSKIP,
     &             RLAT,RLON,WD,WZ,XP,YP)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END
