NCEPLIBS-sp  2.3.3
sptgpm.f
Go to the documentation of this file.
1 C> @file
2 C>
3 C> Transform spectral scalar to mercator
4 C> @author IREDELL @date 96-02-29
5 
6 C> This subprogram performs a spherical transform
7 C> from spectral coefficients of scalar quantities
8 C> to scalar fields on a mercator grid.
9 C> The wave-space can be either triangular or rhomboidal.
10 C> The wave and grid fields may have general indexing,
11 C> but each wave field is in sequential 'ibm order',
12 C> i.e. with zonal wavenumber as the slower index.
13 C> The mercator grid is identified by the location
14 C> of its first point and by its respective increments.
15 C> The transforms are all multiprocessed over sector points.
16 C> Transform several fields at a time to improve vectorization.
17 C> Subprogram can be called from a multiprocessing environment.
18 C>
19 C> PROGRAM HISTORY LOG:
20 C> - 96-02-29 IREDELL
21 C> - 1998-12-15 IREDELL OPENMP DIRECTIVES INSERTED
22 C>
23 C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE
24 C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
25 C> @param MAXWV - INTEGER SPECTRAL TRUNCATION
26 C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM.
27 C> @param MI - INTEGER NUMBER OF POINTS IN THE FASTER ZONAL DIRECTION
28 C> @param MJ - INTEGER NUMBER OF POINTS IN THE SLOWER MERID DIRECTION
29 C> @param KWSKIP - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
30 C> (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0)
31 C> @param KGSKIP - INTEGER SKIP NUMBER BETWEEN GRID FIELDS
32 C> (DEFAULTS TO MI*MJ IF KGSKIP=0)
33 C> @param NISKIP - INTEGER SKIP NUMBER BETWEEN GRID I-POINTS
34 C> (DEFAULTS TO 1 IF NISKIP=0)
35 C> @param NJSKIP - INTEGER SKIP NUMBER BETWEEN GRID J-POINTS
36 C> (DEFAULTS TO MI IF NJSKIP=0)
37 C> @param RLAT1 - REAL LATITUDE OF THE FIRST GRID POINT IN DEGREES
38 C> @param RLON1 - REAL LONGITUDE OF THE FIRST GRID POINT IN DEGREES
39 C> @param DLAT - REAL LATITUDE INCREMENT IN DEGREES SUCH THAT
40 C> D(PHI)/D(J)=DLAT*COS(PHI) WHERE J IS MERIDIONAL INDEX.
41 C> DLAT IS NEGATIVE FOR GRIDS INDEXED SOUTHWARD.
42 C> (IN TERMS OF GRID INCREMENT DY VALID AT LATITUDE RLATI,
43 C> THE LATITUDE INCREMENT DLAT IS DETERMINED AS
44 C> DLAT=DPR*DY/(RERTH*COS(RLATI/DPR))
45 C> WHERE DPR=180/PI AND RERTH IS EARTH'S RADIUS)
46 C> @param DLON - REAL LONGITUDE INCREMENT IN DEGREES SUCH THAT
47 C> D(LAMBDA)/D(I)=DLON WHERE I IS ZONAL INDEX.
48 C> DLON IS NEGATIVE FOR GRIDS INDEXED WESTWARD.
49 C> @param WAVE - REAL (*) WAVE FIELDS
50 C> @param GM - REAL (*) MERCATOR FIELDS
51 C>
52 C> SUBPROGRAMS CALLED:
53 C> - SPWGET() GET WAVE-SPACE CONSTANTS
54 C> - SPLEGEND() COMPUTE LEGENDRE POLYNOMIALS
55 C> - SPSYNTH() SYNTHESIZE FOURIER FROM SPECTRAL
56  SUBROUTINE sptgpm(IROMB,MAXWV,KMAX,MI,MJ,
57  & KWSKIP,KGSKIP,NISKIP,NJSKIP,
58  & RLAT1,RLON1,DLAT,DLON,WAVE,GM)
59 
60  REAL WAVE(*),GM(*)
61  REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
62  REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
63  REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
64  REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
65  INTEGER MP(KMAX)
66  REAL WTOP(2*(MAXWV+1),KMAX)
67  REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
68  REAL F(2*MAXWV+3,2,KMAX)
69  REAL CLAT(MJ),SLAT(MJ),CLON(MAXWV,MI),SLON(MAXWV,MI)
70  parameter(rerth=6.3712e6)
71  parameter(pi=3.14159265358979,dpr=180./pi)
72 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 C CALCULATE PRELIMINARY CONSTANTS
74  CALL spwget(iromb,maxwv,eps,epstop,enn1,elonn1,eon,eontop)
75  mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
76  mxtop=maxwv+1
77  idim=2*maxwv+3
78  kw=kwskip
79  kg=kgskip
80  ni=niskip
81  nj=njskip
82  IF(kw.EQ.0) kw=2*mx
83  IF(kg.EQ.0) kg=mi*mj
84  IF(ni.EQ.0) ni=1
85  IF(nj.EQ.0) nj=mi
86  DO i=1,mi
87  rlon=mod(rlon1+dlon*(i-1)+3600,360.)
88  DO l=1,maxwv
89  clon(l,i)=cos(l*rlon/dpr)
90  slon(l,i)=sin(l*rlon/dpr)
91  ENDDO
92  ENDDO
93  ye=1-log(tan((rlat1+90)/2/dpr))*dpr/dlat
94  DO j=1,mj
95  rlat=atan(exp(dlat/dpr*(j-ye)))*2*dpr-90
96  clat(j)=cos(rlat/dpr)
97  slat(j)=sin(rlat/dpr)
98  ENDDO
99  mp=0
100 C$OMP PARALLEL DO
101  DO k=1,kmax
102  wtop(1:2*mxtop,k)=0
103  ENDDO
104 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105 C TRANSFORM TO GRID
106 C$OMP PARALLEL DO PRIVATE(PLN,PLNTOP,F,IJK)
107  DO j=1,mj
108  CALL splegend(iromb,maxwv,slat(j),clat(j),eps,epstop,
109  & pln,plntop)
110  CALL spsynth(iromb,maxwv,2*maxwv,idim,kw,2*mxtop,kmax,
111  & clat(j),pln,plntop,mp,wave,wtop,f)
112  DO k=1,kmax
113  DO i=1,mi
114  ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
115  gm(ijk)=f(1,1,k)
116  ENDDO
117  DO l=1,maxwv
118  DO i=1,mi
119  ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
120  gm(ijk)=gm(ijk)+2.*(f(2*l+1,1,k)*clon(l,i)
121  & -f(2*l+2,1,k)*slon(l,i))
122  ENDDO
123  ENDDO
124  ENDDO
125  ENDDO
126 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127  END
spwget
subroutine spwget(IROMB, MAXWV, EPS, EPSTOP, ENN1, ELONN1, EON, EONTOP)
This subprogram gets wave-space constants.
Definition: spwget.f:22
splegend
subroutine splegend(I, M, SLAT, CLAT, EPS, EPSTOP, PLN, PLNTOP)
Evaluates the orthonormal associated legendre polynomials in the spectral domain at a given latitude.
Definition: splegend.f:45
spsynth
subroutine spsynth(I, M, IM, IX, NC, NCTOP, KM, CLAT, PLN, PLNTOP, MP, SPC, SPCTOP, F)
SYNTHESIZES FOURIER COEFFICIENTS FROM SPECTRAL COEFFICIENTS FOR A LATITUDE PAIR (NORTHERN AND SOUTHER...
Definition: spsynth.f:35
sptgpm
subroutine sptgpm(IROMB, MAXWV, KMAX, MI, MJ, KWSKIP, KGSKIP, NISKIP, NJSKIP, RLAT1, RLON1, DLAT, DLON, WAVE, GM)
This subprogram performs a spherical transform from spectral coefficients of scalar quantities to sca...
Definition: sptgpm.f:59