sptgpmv.f

Go to the documentation of this file.

C> @file
C>
C> Transform spectral vector to mercator
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> THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL.
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> THE MERCATOR GRID IS IDENTIFIED BY THE LOCATION
C> OF ITS FIRST POINT AND BY ITS RESPECTIVE INCREMENTS.
C> THE TRANSFORMS ARE ALL MULTIPROCESSED OVER SECTOR POINTS.
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>
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 MI       - INTEGER NUMBER OF POINTS IN THE FASTER ZONAL DIRECTION
C> @param MJ       - INTEGER NUMBER OF POINTS IN THE SLOWER MERID DIRECTION
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 GRID FIELDS
C>                (DEFAULTS TO MI*MJ IF KGSKIP=0)
C> @param NISKIP   - INTEGER SKIP NUMBER BETWEEN GRID I-POINTS
C>                (DEFAULTS TO 1 IF NISKIP=0)
C> @param NJSKIP   - INTEGER SKIP NUMBER BETWEEN GRID J-POINTS
C>                (DEFAULTS TO MI IF NJSKIP=0)
C> @param RLAT1    - REAL LATITUDE OF THE FIRST GRID POINT IN DEGREES
C> @param RLON1    - REAL LONGITUDE OF THE FIRST GRID POINT IN DEGREES
C> @param DLAT     - REAL 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     - REAL 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    - REAL (*) WAVE DIVERGENCE FIELDS
C> @param WAVEZ    - REAL (*) WAVE VORTICITY FIELDS
C> @param UM       - REAL (*) MERCATOR U-WINDS
C> @param VM       - REAL (*) MERCATOR V-WINDS
C>
C> SUBPROGRAMS CALLED:
C>   - SPWGET()       GET WAVE-SPACE CONSTANTS
C>   - SPLEGEND()     COMPUTE LEGENDRE POLYNOMIALS
C>   - SPSYNTH()      SYNTHESIZE FOURIER FROM SPECTRAL
C>   - SPDZ2UV()      COMPUTE WINDS FROM DIVERGENCE AND VORTICITY
      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(rerth=6.3712e6)
      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