sptgptvd.f

Go to the documentation of this file.

C> @file
C> @brief Transform spectral vector to station points.
C>
C> ### Program History Log
C> Date | Programmer | Comments
C> -----|------------|---------
C> 96-02-29 | Iredell | Initial.
C> 1998-12-15 | Iredell | Openmp directives inserted.
C> 1999-08-18 | Iredell | Openmp directive typo fixed.
C>
C> @author Iredell @date 96-02-29
 
C> This subprogram performs a spherical transform
C> from spectral coefficients of divergences and curls
C> to specified sets of station point vectors and their
C> gradients on the globe.
C>
C> <pre>
C> DP=(D(UP)/DLON+D(VP*CLAT)/DLAT)/(R*CLAT)
C> ZP=(D(VP)/DLON-D(UP*CLAT)/DLAT)/(R*CLAT)
C> UXP=D(UP*CLAT)/DLON/(R*CLAT)
C> VXP=D(VP*CLAT)/DLON/(R*CLAT)
C> UYP=D(UP*CLAT)/DLAT/R
C> VYP=D(VP*CLAT)/DLAT/R
C> </pre>
C>
C> The wave-space can be either triangular or rhomboidal.
C>
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>
C> The transforms are all multiprocessed over stations.
C>
C> Transform several fields at a time to improve vectorization.
C>
C> Subprogram can be called from a multiprocessing environment.
C>
C> @param IROMB spectral domain shape
C> (0 for triangular, 1 for rhomboidal)
C> @param MAXWV spectral truncation
C> @param KMAX number of fields to transform.
C> @param NMAX number of station points to return
C> @param KWSKIP skip number between wave fields
C> (defaults to (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0)
C> @param KGSKIP skip number between station point sets
C> (defaults to NMAX if KGSKIP=0)
C> @param NRSKIP skip number between station lats and lons
C> (defaults to 1 if NRSKIP=0)
C> @param NGSKIP skip number between station points
C> (defaults to 1 if NGSKIP=0)
C> @param RLAT station latitudes in degrees
C> @param RLON station longitudes in degrees
C> @param WAVED wave divergence fields
C> @param WAVEZ wave vorticity fields
C> @param DP station point divergence sets
C> @param ZP station point vorticity sets
C> @param UP station point u-wind sets
C> @param VP station point v-wind sets
C> @param UXP station point u-wind x-gradient sets
C> @param VXP station point v-wind x-gradient sets
C> @param UYP station point u-wind y-gradient sets
C> @param VYP station point v-wind y-gradient sets
C>
C> @author Iredell @date 96-02-29
      SUBROUTINE sptgptvd(IROMB,MAXWV,KMAX,NMAX,
     &                    KWSKIP,KGSKIP,NRSKIP,NGSKIP,
     &                    RLAT,RLON,WAVED,WAVEZ,
     &                    DP,ZP,UP,VP,UXP,VXP,UYP,VYP)
 
      REAL RLAT(*),RLON(*),WAVED(*),WAVEZ(*)
      REAL DP(*),ZP(*),UP(*),VP(*),UXP(*),VXP(*),UYP(*),VYP(*)
      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(4*KMAX)
      REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,4*KMAX)
      REAL WTOP(2*(MAXWV+1),4*KMAX)
      REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
      REAL F(2*MAXWV+2,2,6*KMAX),G(6*KMAX)
      parameter(pi=3.14159265358979)
 
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
      mdim=2*mx
      idim=2*maxwv+2
      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(1:2*kmax)=0
      mp(2*kmax+1:4*kmax)=1
 
C  CALCULATE SPECTRAL WINDS
C$OMP PARALLEL DO PRIVATE(KWS,KD,KZ,KU,KV)
      DO k=1,kmax
        kws=(k-1)*kw
        kd=0*kmax+k
        kz=1*kmax+k
        ku=2*kmax+k
        kv=3*kmax+k
        DO i=1,2*mx
          w(i,kd)=waved(kws+i)
          w(i,kz)=wavez(kws+i)
        ENDDO
        DO i=1,2*mxtop
          wtop(i,kd)=0
          wtop(i,kz)=0
        ENDDO
        CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
     &               waved(kws+1),wavez(kws+1),
     &               w(1,ku),w(1,kv),wtop(1,ku),wtop(1,kv))
      ENDDO
 
C  CALCULATE STATION FIELDS
C$OMP PARALLEL DO PRIVATE(KD,KZ,KU,KV,KUX,KVX,SLAT1,CLAT1)
C$OMP&            PRIVATE(PLN,PLNTOP,F,G,NK)
      DO n=1,nmax
        ku=2*kmax+1
        kux=4*kmax+1
        IF(abs(rlat((n-1)*nr+1)).GE.89.9995) THEN
          slat1=sign(1.,rlat((n-1)*nr+1))
          clat1=0.
        ELSE
          slat1=sin(pi/180*rlat((n-1)*nr+1))
          clat1=cos(pi/180*rlat((n-1)*nr+1))
        ENDIF
        CALL splegend(iromb,maxwv,slat1,clat1,eps,epstop,
     &                pln,plntop)
        CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,4*kmax,
     &               clat1,pln,plntop,mp,w,wtop,f)
        CALL spgradx(maxwv,idim,2*kmax,mp(2*kmax+1),clat1,
     &               f(1,1,2*kmax+1),f(1,1,4*kmax+1))
        CALL spfftpt(maxwv,1,idim,1,6*kmax,rlon((n-1)*nr+1),f,g)
        DO k=1,kmax
          kd=0*kmax+k
          kz=1*kmax+k
          ku=2*kmax+k
          kv=3*kmax+k
          kux=4*kmax+k
          kvx=5*kmax+k
          nk=(n-1)*ng+(k-1)*kg+1
          dp(nk)=g(kd)
          zp(nk)=g(kz)
          up(nk)=g(ku)
          vp(nk)=g(kv)
          uxp(nk)=g(kux)
          vxp(nk)=g(kvx)
          uyp(nk)=g(kvx)-clat1*g(kz)
          vyp(nk)=clat1*g(kd)-g(kux)
        ENDDO
      ENDDO
      END