NCEPLIBS-sp  2.3.3
sptrand.f
Go to the documentation of this file.
1 C> @file
2 C>
3 C> Perform a gradient spherical transform
4 C> @author IREDELL @date 96-02-29
5 
6 C> THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM
7 C> BETWEEN SPECTRAL COEFFICIENTS OF SCALAR FIELDS
8 C> AND THEIR MEANS AND GRADIENTS ON A GLOBAL CYLINDRICAL GRID.
9 C> THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL.
10 C> THE GRID-SPACE CAN BE EITHER AN EQUALLY-SPACED GRID
11 C> (WITH OR WITHOUT POLE POINTS) OR A GAUSSIAN GRID.
12 C> THE WAVE AND GRID FIELDS MAY HAVE GENERAL INDEXING,
13 C> BUT EACH WAVE FIELD IS IN SEQUENTIAL 'IBM ORDER',
14 C> I.E. WITH ZONAL WAVENUMBER AS THE SLOWER INDEX.
15 C> TRANSFORMS ARE DONE IN LATITUDE PAIRS FOR EFFICIENCY;
16 C> THUS GRID ARRAYS FOR EACH HEMISPHERE MUST BE PASSED.
17 C> IF SO REQUESTED, JUST A SUBSET OF THE LATITUDE PAIRS
18 C> MAY BE TRANSFORMED IN EACH INVOCATION OF THE SUBPROGRAM.
19 C> THE TRANSFORMS ARE ALL MULTIPROCESSED OVER LATITUDE EXCEPT
20 C> THE TRANSFORM FROM FOURIER TO SPECTRAL IS MULTIPROCESSED
21 C> OVER ZONAL WAVENUMBER TO ENSURE REPRODUCIBILITY.
22 C> TRANSFORM SEVERAL FIELDS AT A TIME TO IMPROVE VECTORIZATION.
23 C> SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
24 C>
25 C> PROGRAM HISTORY LOG:
26 C> - 96-02-29 IREDELL
27 C> - 1998-12-15 IREDELL OPENMP DIRECTIVES INSERTED
28 C>
29 C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE
30 C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
31 C> @param MAXWV - INTEGER SPECTRAL TRUNCATION
32 C> @param IDRT - INTEGER GRID IDENTIFIER
33 C> (IDRT=4 FOR GAUSSIAN GRID,
34 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
35 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
36 C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES.
37 C> @param JMAX - INTEGER NUMBER OF LATITUDES.
38 C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM.
39 C> @param IPRIME - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN.
40 C> (DEFAULTS TO 1 IF IPRIME=0)
41 C> @param ISKIP - INTEGER SKIP NUMBER BETWEEN LONGITUDES
42 C> (DEFAULTS TO 1 IF ISKIP=0)
43 C> @param JNSKIP - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH
44 C> (DEFAULTS TO IMAX IF JNSKIP=0)
45 C> @param JSSKIP - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH
46 C> (DEFAULTS TO -IMAX IF JSSKIP=0)
47 C> @param KWSKIP - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
48 C> (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0)
49 C> @param KGSKIP - INTEGER SKIP NUMBER BETWEEN GRID FIELDS
50 C> (DEFAULTS TO IMAX*JMAX IF KGSKIP=0)
51 C> @param JBEG - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM
52 C> (DEFAULTS TO 1 IF JBEG=0)
53 C> (IF JBEG=0 AND IDIR<0, WAVE IS ZEROED BEFORE TRANSFORM)
54 C> @param JEND - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM
55 C> (DEFAULTS TO (JMAX+1)/2 IF JEND=0)
56 C> @param JCPU - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS
57 C> @param[out] WAVE - REAL (*) WAVE FIELDS IF IDIR>0
58 C> @param[out] GRIDMN - REAL (KMAX) GLOBAL MEANS IF IDIR<0
59 C> @param[out] GRIDXN - REAL (*) N.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR<0
60 C> @param[out] GRIDXS - REAL (*) S.H. X-GRADIENTS (STARTING AT JBEG) IF IDIR<0
61 C> [GRIDX=(D(WAVE)/DLAM)/(CLAT*RERTH)]
62 C> @param[out] GRIDYN - REAL (*) N.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR<0
63 C> @param[out] GRIDYS - REAL (*) S.H. Y-GRADIENTS (STARTING AT JBEG) IF IDIR<0
64 C> [GRIDY=(D(WAVE)/DPHI)/RERTH]
65 C> @param IDIR - INTEGER TRANSFORM FLAG
66 C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
67 C>
68 C> SUBPROGRAMS CALLED:
69 C> - SPWGET GET WAVE-SPACE CONSTANTS
70 C> - SPLAPLAC COMPUTE LAPLACIAN IN SPECTRAL SPACE
71 C> - SPTRANV PERFORM A VECTOR SPHERICAL TRANSFORM
72 C>
73 C> REMARKS: MINIMUM GRID DIMENSIONS FOR UNALIASED TRANSFORMS TO SPECTRAL:
74 C> DIMENSION |LINEAR |QUADRATIC
75 C> ----------------------- |--------- |-------------
76 C> IMAX |2*MAXWV+2 |3*MAXWV/2*2+2
77 C> JMAX (IDRT=4,IROMB=0) |1*MAXWV+1 |3*MAXWV/2+1
78 C> JMAX (IDRT=4,IROMB=1) |2*MAXWV+1 |5*MAXWV/2+1
79 C> JMAX (IDRT=0,IROMB=0) |2*MAXWV+3 |3*MAXWV/2*2+3
80 C> JMAX (IDRT=0,IROMB=1) |4*MAXWV+3 |5*MAXWV/2*2+3
81 C> JMAX (IDRT=256,IROMB=0) |2*MAXWV+1 |3*MAXWV/2*2+1
82 C> JMAX (IDRT=256,IROMB=1) |4*MAXWV+1 |5*MAXWV/2*2+1
83  SUBROUTINE sptrand(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
84  & IPRIME,ISKIP,JNSKIP,JSSKIP,KWSKIP,KGSKIP,
85  & JBEG,JEND,JCPU,
86  & WAVE,GRIDMN,GRIDXN,GRIDXS,GRIDYN,GRIDYS,IDIR)
87 
88  REAL WAVE(*),GRIDMN(KMAX),GRIDXN(*),GRIDXS(*),GRIDYN(*),GRIDYS(*)
89  REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
90  REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
91  REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
92  REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
93  REAL WD((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,KMAX)
94  REAL WZ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,KMAX)
95 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 C SET PARAMETERS
97  CALL spwget(iromb,maxwv,eps,epstop,enn1,elonn1,eon,eontop)
98  mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
99  mdim=2*mx+1
100  kw=kwskip
101  IF(kw.EQ.0) kw=2*mx
102 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103 C TRANSFORM WAVE TO GRID
104  IF(idir.GT.0) THEN
105 C$OMP PARALLEL DO PRIVATE(KWS)
106  DO k=1,kmax
107  kws=(k-1)*kw
108  gridmn(k)=wave(kws+1)/sqrt(2.)
109  CALL splaplac(iromb,maxwv,enn1,wave(kws+1),wd(1,k),1)
110  wz(1:2*mx,k)=0.
111  ENDDO
112  CALL sptranv(iromb,maxwv,idrt,imax,jmax,kmax,
113  & iprime,iskip,jnskip,jsskip,mdim,kgskip,
114  & jbeg,jend,jcpu,
115  & wd,wz,gridxn,gridxs,gridyn,gridys,idir)
116 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117 C TRANSFORM GRID TO WAVE
118  ELSE
119 C$OMP PARALLEL DO
120  DO k=1,kmax
121  wd(1:2*mx,k)=0.
122  wz(1:2*mx,k)=0.
123  ENDDO
124  CALL sptranv(iromb,maxwv,idrt,imax,jmax,kmax,
125  & iprime,iskip,jnskip,jsskip,mdim,kgskip,
126  & jbeg,jend,jcpu,
127  & wd,wz,gridxn,gridxs,gridyn,gridys,idir)
128  IF(jbeg.EQ.0) THEN
129 C$OMP PARALLEL DO PRIVATE(KWS)
130  DO k=1,kmax
131  kws=(k-1)*kw
132  CALL splaplac(iromb,maxwv,enn1,wave(kws+1),wd(1,k),-1)
133  wave(kws+1)=gridmn(k)*sqrt(2.)
134  ENDDO
135  ELSE
136 C$OMP PARALLEL DO PRIVATE(KWS)
137  DO k=1,kmax
138  kws=(k-1)*kw
139  CALL splaplac(iromb,maxwv,enn1,wz(1,k),wd(1,k),-1)
140  wave(kws+1:kws+2*mx)=wave(kws+1:kws+2*mx)+wz(1:2*mx,k)
141  wave(kws+1)=gridmn(k)*sqrt(2.)
142  ENDDO
143  ENDIF
144  ENDIF
145 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146  END
spwget
subroutine spwget(IROMB, MAXWV, EPS, EPSTOP, ENN1, ELONN1, EON, EONTOP)
This subprogram gets wave-space constants.
Definition: spwget.f:22
splaplac
subroutine splaplac(I, M, ENN1, Q, QD2, IDIR)
COMPUTES THE LAPLACIAN OR THE INVERSE LAPLACIAN OF A SCALAR FIELD IN SPECTRAL SPACE.
Definition: splaplac.f:22
sptranv
subroutine sptranv(IROMB, MAXWV, IDRT, IMAX, JMAX, KMAX, IPRIME, ISKIP, JNSKIP, JSSKIP, KWSKIP, KGSKIP, JBEG, JEND, JCPU, WAVED, WAVEZ, GRIDUN, GRIDUS, GRIDVN, GRIDVS, IDIR)
THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM BETWEEN SPECTRAL COEFFICIENTS OF DIVERGENCES AND CURLS...
Definition: sptranv.f:85
sptrand
subroutine sptrand(IROMB, MAXWV, IDRT, IMAX, JMAX, KMAX, IPRIME, ISKIP, JNSKIP, JSSKIP, KWSKIP, KGSKIP, JBEG, JEND, JCPU, WAVE, GRIDMN, GRIDXN, GRIDXS, GRIDYN, GRIDYS, IDIR)
THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM BETWEEN SPECTRAL COEFFICIENTS OF SCALAR FIELDS AND THE...
Definition: sptrand.f:87