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 1998-12-15 IREDELL OPENMP DIRECTIVES INSERTED
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
C$OMP PARALLEL DO PRIVATE(KWS)
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)
WZ(1:2*MX,K)=0.
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
C$OMP PARALLEL DO
DO K=1,KMAX
WD(1:2*MX,K)=0.
WZ(1:2*MX,K)=0.
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
C$OMP PARALLEL DO PRIVATE(KWS)
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
C$OMP PARALLEL DO PRIVATE(KWS)
DO K=1,KMAX
KWS=(K-1)*KW
CALL SPLAPLAC(IROMB,MAXWV,ENN1,WZ(1,K),WD(1,K),-1)
WAVE(KWS+1:KWS+2*MX)=WAVE(KWS+1:KWS+2*MX)+WZ(1:2*MX,K)
WAVE(KWS+1)=GRIDMN(K)*SQRT(2.)
ENDDO
ENDIF
ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END