sptgpm.f

Go to the documentation of this file.

C> @file
C> @brief Transform spectral scalar 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 scalar quantities
C> to scalar 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> @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 WAVE Wave fields
C> @param GM Mercator fields
C>
C> @author IREDELL @date 96-02-29
      SUBROUTINE sptgpm(IROMB,MAXWV,KMAX,MI,MJ,
     &                  KWSKIP,KGSKIP,NISKIP,NJSKIP,
     &                  RLAT1,RLON1,DLAT,DLON,WAVE,GM)
 
      REAL WAVE(*),GM(*)
      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)
      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
      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=0
C$OMP PARALLEL DO
      DO k=1,kmax
        wtop(1:2*mxtop,k)=0
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  TRANSFORM TO GRID
C$OMP PARALLEL DO PRIVATE(PLN,PLNTOP,F,IJK)
      DO j=1,mj
        CALL splegend(iromb,maxwv,slat(j),clat(j),eps,epstop,
     &                pln,plntop)
        CALL spsynth(iromb,maxwv,2*maxwv,idim,kw,2*mxtop,kmax,
     &               clat(j),pln,plntop,mp,wave,wtop,f)
        DO k=1,kmax
          DO i=1,mi
            ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
            gm(ijk)=f(1,1,k)
          ENDDO
          DO l=1,maxwv
            DO i=1,mi
              ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
              gm(ijk)=gm(ijk)+2.*(f(2*l+1,1,k)*clon(l,i)
     &                           -f(2*l+2,1,k)*slon(l,i))
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      END