C-----------------------------------------------------------------------
      SUBROUTINE SPTRAND(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
     &                   IPRIME,ISKIP,JNSKIP,JSSKIP,KWSKIP,KGSKIP,
     &                   JBEG,JEND,JCPU,
     &                   WAVE,GRIDMN,GRIDXN,GRIDXS,GRIDYN,GRIDYS,IDIR)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:  SPTRAND    PERFORM A GRADIENT SPHERICAL TRANSFORM
C   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-02-29
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM
C           BETWEEN SPECTRAL COEFFICIENTS OF SCALAR FIELDS
C           AND THEIR MEANS AND GRADIENTS ON A GLOBAL CYLINDRICAL GRID.
C           THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL.
C           THE GRID-SPACE CAN BE EITHER AN EQUALLY-SPACED GRID
C           (WITH OR WITHOUT POLE POINTS) OR A GAUSSIAN GRID.
C           THE WAVE AND GRID 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           TRANSFORMS ARE DONE IN LATITUDE PAIRS FOR EFFICIENCY;
C           THUS GRID ARRAYS FOR EACH HEMISPHERE MUST BE PASSED.
C           IF SO REQUESTED, JUST A SUBSET OF THE LATITUDE PAIRS
C           MAY BE TRANSFORMED IN EACH INVOCATION OF THE SUBPROGRAM.
C           THE TRANSFORMS ARE ALL MULTIPROCESSED OVER LATITUDE EXCEPT
C           THE TRANSFORM FROM FOURIER TO SPECTRAL IS MULTIPROCESSED
C           OVER ZONAL WAVENUMBER TO ENSURE REPRODUCIBILITY.
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 SPTRAND(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
C    &                   IPRIME,ISKIP,JNSKIP,JSSKIP,KWSKIP,KGSKIP,
C    &                   JBEG,JEND,JCPU,
C    &                   WAVE,GRIDMN,GRIDXN,GRIDXS,GRIDYN,GRIDYS,IDIR)
C   INPUT ARGUMENTS:
C     IROMB    - INTEGER SPECTRAL DOMAIN SHAPE
C                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C     MAXWV    - INTEGER SPECTRAL TRUNCATION
C     IDRT     - INTEGER GRID IDENTIFIER
C                (IDRT=4 FOR GAUSSIAN GRID,
C                 IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
C                 IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
C     IMAX     - INTEGER EVEN NUMBER OF LONGITUDES.
C     JMAX     - INTEGER NUMBER OF LATITUDES.
C     KMAX     - INTEGER NUMBER OF FIELDS TO TRANSFORM.
C     IPRIME   - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN.
C                (DEFAULTS TO 1 IF IPRIME=0)
C     ISKIP    - INTEGER SKIP NUMBER BETWEEN LONGITUDES
C                (DEFAULTS TO 1 IF ISKIP=0)
C     JNSKIP   - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH
C                (DEFAULTS TO IMAX IF JNSKIP=0)
C     JSSKIP   - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH
C                (DEFAULTS TO -IMAX IF JSSKIP=0)
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 GRID FIELDS
C                (DEFAULTS TO IMAX*JMAX IF KGSKIP=0)
C     JBEG     - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM
C                (DEFAULTS TO 1 IF JBEG=0)
C                (IF JBEG=0 AND IDIR<0, WAVE IS ZEROED BEFORE TRANSFORM)
C     JEND     - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM
C                (DEFAULTS TO (JMAX+1)/2 IF JEND=0)
C     JCPU     - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS
C     WAVE     - REAL (*) WAVE FIELDS IF IDIR>0
C     GRIDMN   - REAL (KMAX) GLOBAL MEANS IF IDIR<0
C     GRIDXN   - REAL (*) N.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR<0
C     GRIDXS   - REAL (*) S.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR<0
C     GRIDYN   - REAL (*) N.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR<0
C     GRIDYS   - REAL (*) S.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR<0
C     IDIR     - INTEGER TRANSFORM FLAG
C                (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
C   OUTPUT ARGUMENTS:
C     WAVE     - REAL (*) WAVE FIELDS IF IDIR<0
C     GRIDMN   - REAL (KMAX) GLOBAL MEANS IF IDIR>0
C     GRIDXN   - REAL (*) N.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR>0
C     GRIDXS   - REAL (*) S.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR>0
C                [GRIDX=(D(WAVE)/DLAM)/(CLAT*RERTH)]
C     GRIDYN   - REAL (*) N.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR>0
C     GRIDYS   - REAL (*) S.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR>0
C                [GRIDY=(D(WAVE)/DPHI)/RERTH]
C
C SUBPROGRAMS CALLED:
C   SPWGET       GET WAVE-SPACE CONSTANTS
C   SPLAPLAC     COMPUTE LAPLACIAN IN SPECTRAL SPACE
C   SPTRANV      PERFORM A VECTOR SPHERICAL TRANSFORM
C
C REMARKS: MINIMUM GRID DIMENSIONS FOR UNALIASED TRANSFORMS TO SPECTRAL:
C   DIMENSION                    LINEAR              QUADRATIC
C   -----------------------      ---------           -------------
C   IMAX                         2*MAXWV+2           3*MAXWV/2*2+2
C   JMAX (IDRT=4,IROMB=0)        1*MAXWV+1           3*MAXWV/2+1
C   JMAX (IDRT=4,IROMB=1)        2*MAXWV+1           5*MAXWV/2+1
C   JMAX (IDRT=0,IROMB=0)        2*MAXWV+3           3*MAXWV/2*2+3
C   JMAX (IDRT=0,IROMB=1)        4*MAXWV+3           5*MAXWV/2*2+3
C   JMAX (IDRT=256,IROMB=0)      2*MAXWV+1           3*MAXWV/2*2+1
C   JMAX (IDRT=256,IROMB=1)      4*MAXWV+1           5*MAXWV/2*2+1
C   -----------------------      ---------           -------------
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C
C$$$
      REAL WAVE(*),GRIDMN(KMAX),GRIDXN(*),GRIDXS(*),GRIDYN(*),GRIDYS(*)
      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  SET PARAMETERS
      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  TRANSFORM WAVE TO GRID
      IF(IDIR.GT.0) THEN
CMIC$ DO ALL PRIVATE(K,KWS)
CMIC$&       SHARED(KMAX,GRIDMN,IROMB,MAXWV,ENN1,KW)
CMIC$&       SHARED(WAVE,WD,WZ,W,WTOP,MX)
        DO K=1,KMAX
          KWS=(K-1)*KW
          GRIDMN(K)=WAVE(KWS+1)/SQRT(2.)
          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.
          enddo
        ENDDO
        CALL SPTRANV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
     &               IPRIME,ISKIP,JNSKIP,JSSKIP,MDIM,KGSKIP,
     &               JBEG,JEND,JCPU,
     &               WD,WZ,GRIDXN,GRIDXS,GRIDYN,GRIDYS,IDIR)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSFORM GRID TO WAVE
      ELSE
CMIC$ DO ALL PRIVATE(K)
CMIC$&       SHARED(KMAX,MX,WD,WZ)
        DO K=1,KMAX
c         WD(1:2*MX,K)=0.
c         WZ(1:2*MX,K)=0.
          do ii = 1, 2*MX
             WD(ii,K)=0.
             WZ(ii,K)=0.
          enddo
        ENDDO
        CALL SPTRANV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
     &               IPRIME,ISKIP,JNSKIP,JSSKIP,MDIM,KGSKIP,
     &               JBEG,JEND,JCPU,
     &               WD,WZ,GRIDXN,GRIDXS,GRIDYN,GRIDYS,IDIR)
        IF(JBEG.EQ.0) THEN
CMIC$ DO ALL PRIVATE(K,KWS)
CMIC$&       SHARED(KMAX,IROMB,MAXWV,ENN1,W,WTOP,KW,WAVE,WD,WZ,GRIDMN)
          DO K=1,KMAX
            KWS=(K-1)*KW
            CALL SPLAPLAC(IROMB,MAXWV,ENN1,WAVE(KWS+1),WD(1,K),-1)
            WAVE(KWS+1)=GRIDMN(K)*SQRT(2.)
          ENDDO
        ELSE
CMIC$ DO ALL PRIVATE(K,KWS)
CMIC$&       SHARED(KMAX,IROMB,MAXWV,ENN1,W,WTOP,KW)
CMIC$&       SHARED(WAVE,WD,WZ,GRIDMN,MX)
          DO K=1,KMAX
            KWS=(K-1)*KW
            CALL SPLAPLAC(IROMB,MAXWV,ENN1,WZ(1,K),WD(1,K),-1)
c           WAVE(KWS+1:KWS+2*MX)=WAVE(KWS+1:KWS+2*MX)+WZ(1:2*MX,K)
            do ii = 1, 2*MX
               WAVE(KWS+ii)=WAVE(KWS+ii)+WZ(ii,K)
            enddo
            WAVE(KWS+1)=GRIDMN(K)*SQRT(2.)
          ENDDO
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END
