sptgpt.f

Go to the documentation of this file.

C> @file
C>
C> Transform spectral scalar to station points
C> @author IREDELL @date 96-02-29
 
C> This subprogram performs a spherical transform
c> from spectral coefficients of scalar quantities
c> to specified sets of station points 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> - 1998-12-15  IREDELL  OPENMP DIRECTIVES INSERTED
C> - 2003-06-30  IREDELL  USE SPFFTPT
C>
C> @param IROMB    - INTEGER SPECTRAL DOMAIN SHAPE
C>                (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
C> @param MAXWV    - INTEGER SPECTRAL TRUNCATION
C> @param KMAX     - INTEGER NUMBER OF FIELDS TO TRANSFORM.
C> @param NMAX     - INTEGER NUMBER OF STATION POINTS TO RETURN
C> @param KWSKIP   - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
C>                (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0)
C> @param KGSKIP   - INTEGER SKIP NUMBER BETWEEN STATION POINT SETS
C>                (DEFAULTS TO NMAX IF KGSKIP=0)
C> @param NRSKIP   - INTEGER SKIP NUMBER BETWEEN STATION LATS AND LONS
C>                (DEFAULTS TO 1 IF NRSKIP=0)
C> @param NGSKIP   - INTEGER SKIP NUMBER BETWEEN STATION POINTS
C>                (DEFAULTS TO 1 IF NGSKIP=0)
C> @param RLAT     - REAL (*) STATION LATITUDES IN DEGREES
C> @param RLON     - REAL (*) STATION LONGITUDES IN DEGREES
C> @param WAVE     - REAL (*) WAVE FIELDS
C> @param GP       - REAL (*) STATION POINT SETS
C>
C> SUBPROGRAMS CALLED:
C>   - SPWGET()       GET WAVE-SPACE CONSTANTS
C>   - SPLEGEND()     COMPUTE LEGENDRE POLYNOMIALS
C>   - SPSYNTH()      SYNTHESIZE FOURIER FROM SPECTRAL
C>   - SPFFTPT()      POINTWISE FOURIER TRANSFORM
      SUBROUTINE sptgpt(IROMB,MAXWV,KMAX,NMAX,
     &                  KWSKIP,KGSKIP,NRSKIP,NGSKIP,
     &                  RLAT,RLON,WAVE,GP)
 
      REAL RLAT(*),RLON(*),WAVE(*),GP(*)
      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)
      INTEGER MP(KMAX)
      REAL WTOP(2*(MAXWV+1),KMAX)
      REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
      REAL F(2*MAXWV+3,2,KMAX)
      parameter(pi=3.14159265358979)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE PRELIMINARY CONSTANTS
      CALL spwget(iromb,maxwv,eps,epstop,enn1,elonn1,eon,eontop)
      mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
      mxtop=maxwv+1
      idim=2*maxwv+3
      kw=kwskip
      kg=kgskip
      nr=nrskip
      ng=ngskip
      IF(kw.EQ.0) kw=2*mx
      IF(kg.EQ.0) kg=nmax
      IF(nr.EQ.0) nr=1
      IF(ng.EQ.0) ng=1
      mp=0
C$OMP PARALLEL DO
      DO k=1,kmax
        wtop(1:2*mxtop,k)=0
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE STATION FIELDS
C$OMP PARALLEL DO PRIVATE(RADLAT,SLAT1,CLAT1)
C$OMP&            PRIVATE(PLN,PLNTOP,F,NK)
      DO n=1,nmax
        radlat=pi/180*rlat((n-1)*nr+1)
        IF(rlat((n-1)*nr+1).GE.89.9995) THEN
          slat1=1.
          clat1=0.
        ELSEIF(rlat((n-1)*nr+1).LE.-89.9995) THEN
          slat1=-1.
          clat1=0.
        ELSE
          slat1=sin(radlat)
          clat1=cos(radlat)
        ENDIF
        CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
     &                pln,plntop)
        CALL spsynth(iromb,maxwv,2*maxwv,idim,kw,2*mxtop,kmax,
     &               clat1,pln,plntop,mp,wave,wtop,f)
        CALL spfftpt(maxwv,1,2*maxwv+3,kg,kmax,rlon((n-1)*nr+1),
     &               f,gp((n-1)*ng+1))
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END