sptgpmv.f

Go to the documentation of this file.

C> @file
C> @brief Transform spectral vector to Mercator
C> ### Program history log:
C> Date | Programmer | Comments
C> -----|------------|----------
C> 96-02-29 | IREDELL | Initial.
C> 1998-12-15 | IREDELL | OpenMP directives inserted.
C> @author IREDELL @date 96-02-29
 
C> This subprogram performs a spherical transform
C> from spectral coefficients of divergences and curls
C> to vector fields on a Mercator grid.
C>
C> The wave-space can be either triangular or rhomboidal.
C>
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>
C> The Mercator grid is identified by the location
C> of its first point and by its respective increments.
C>
C> The transforms are all multiprocessed over sector points.
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 MI Number of points in the faster zonal direction
C> @param MJ Number of points in the slower merid direction
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 grid fields
C> (defaults to MI*MJ if KGSKIP=0)
C> @param NISKIP Skip number between grid i-points
C> (defaults to 1 if NISKIP=0)
C> @param NJSKIP Skip number between grid j-points
C> (defaults to MI if NJSKIP=0)
C> @param RLAT1 Latitude of the first grid point in degrees
C> @param RLON1 Longitude of the first grid point in degrees
C> @param DLAT Latitude increment in degrees such that
C> D(PHI)/D(J)=DLAT*COS(PHI) where J is meridional index.
C> DLAT is negative for grids indexed southward.
C> (in terms of grid increment dy valid at latitude RLATI,
C> The latitude increment DLAT is determined as
C> DLAT=DPR*DY/(RERTH*COS(RLATI/DPR))
C> where DPR=180/PI and RERTH is Earth's radius)
C> @param DLON longitude increment in degrees such that
C> D(LAMBDA)/D(I)=DLON where I is zonal index.
C> DLON is negative for grids indexed westward.
C> @param WAVED Wave divergence fields
C> @param WAVEZ Wave vorticity fields
C> @param UM Mercator u-winds
C> @param VM Mercator v-winds
C>
C> @author IREDELL @date 96-02-29
      SUBROUTINE sptgpmv(IROMB,MAXWV,KMAX,MI,MJ,
     &                   KWSKIP,KGSKIP,NISKIP,NJSKIP,
     &                   RLAT1,RLON1,DLAT,DLON,WAVED,WAVEZ,UM,VM)
 
      REAL WAVED(*),WAVEZ(*),UM(*),VM(*)
      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(2*KMAX)
      REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,2*KMAX)
      REAL WTOP(2*(MAXWV+1),2*KMAX)
      REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
      REAL F(2*MAXWV+3,2,2*KMAX)
      REAL CLAT(MJ),SLAT(MJ),CLON(MAXWV,MI),SLON(MAXWV,MI)
      parameter(pi=3.14159265358979,dpr=180./pi)
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
      mdim=2*mx+1
      idim=2*maxwv+3
      kw=kwskip
      kg=kgskip
      ni=niskip
      nj=njskip
      IF(kw.EQ.0) kw=2*mx
      IF(kg.EQ.0) kg=mi*mj
      IF(ni.EQ.0) ni=1
      IF(nj.EQ.0) nj=mi
      DO i=1,mi
        rlon=mod(rlon1+dlon*(i-1)+3600,360.)
        DO l=1,maxwv
          clon(l,i)=cos(l*rlon/dpr)
          slon(l,i)=sin(l*rlon/dpr)
        ENDDO
      ENDDO
      ye=1-log(tan((rlat1+90)/2/dpr))*dpr/dlat
      DO j=1,mj
        rlat=atan(exp(dlat/dpr*(j-ye)))*2*dpr-90
        clat(j)=cos(rlat/dpr)
        slat(j)=sin(rlat/dpr)
      ENDDO
      mp=1
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  CALCULATE SPECTRAL WINDS
C$OMP PARALLEL DO PRIVATE(KWS)
      DO k=1,kmax
        kws=(k-1)*kw
        CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
     &               waved(kws+1),wavez(kws+1),
     &               w(1,k),w(1,kmax+k),wtop(1,k),wtop(1,kmax+k))
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSFORM TO GRID
C$OMP PARALLEL DO PRIVATE(PLN,PLNTOP,F,KU,KV,IJK)
      DO j=1,mj
        CALL splegend(iromb,maxwv,slat(j),clat(j),eps,epstop,
     &                pln,plntop)
        CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
     &               clat(j),pln,plntop,mp,w,wtop,f)
        DO k=1,kmax
          ku=k
          kv=k+kmax
          DO i=1,mi
            ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
            um(ijk)=f(1,1,ku)
            vm(ijk)=f(1,1,kv)
          ENDDO
          DO l=1,maxwv
            DO i=1,mi
              ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
              um(ijk)=um(ijk)+2.*(f(2*l+1,1,ku)*clon(l,i)
     &                           -f(2*l+2,1,ku)*slon(l,i))
              vm(ijk)=vm(ijk)+2.*(f(2*l+1,1,kv)*clon(l,i)
     &                           -f(2*l+2,1,kv)*slon(l,i))
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END