sptranfv.f

Go to the documentation of this file.

C> @file
C> @brief Perform a vector spherical transform
C>
C> ### Program History Log
C> Date | Programmer | Comments
C> -----|------------|---------
C> 96-02-29 | Iredell | Initial.
C> 1998-12-15 | Iredell | Generic fft used, openmp directives inserted
C> 2013-01-16 | Iredell & MIRVIS | Fixing afft negative sharing effect during omp loops
C>
C> @author Iredell @date 96-02-29
 
C> This subprogram performs a spherical transform
C> between spectral coefficients of divergences and curls
C> and vector fields on a global cylindrical grid.
C>
C> The wave-space can be either triangular or rhomboidal.
C>
C> The grid-space can be either an equally-spaced grid
C> (with or without pole points) or a Gaussian grid.
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> Transforms are done in latitude pairs for efficiency;
C> thus grid arrays for each hemisphere must be passed.
C> if so requested, just a subset of the latitude pairs
C> may be transformed in each invocation of the subprogram.
C>
C> The transforms are all multiprocessed over latitude except
C> the transform from fourier to spectral is multiprocessed
C> over zonal wavenumber to ensure reproducibility.
C>
C> Transform several fields at a time to improve vectorization.
C> subprogram can be called from a multiprocessing environment.
C>
C> Minimum grid dimensions for unaliased transforms to spectral:
C> DIMENSION                    |LINEAR              |QUADRATIC
C> -----------------------      |---------           |-------------
C> IMAX                         |2*MAXWV+2           |3*MAXWV/2*2+2
C> JMAX (IDRT=4,IROMB=0)        |1*MAXWV+1           |3*MAXWV/2+1
C> JMAX (IDRT=4,IROMB=1)        |2*MAXWV+1           |5*MAXWV/2+1
C> JMAX (IDRT=0,IROMB=0)        |2*MAXWV+3           |3*MAXWV/2*2+3
C> JMAX (IDRT=0,IROMB=1)        |4*MAXWV+3           |5*MAXWV/2*2+3
C> JMAX (IDRT=256,IROMB=0)      |2*MAXWV+1           |3*MAXWV/2*2+1
C> JMAX (IDRT=256,IROMB=1)      |4*MAXWV+1           |5*MAXWV/2*2+1
C>
C> @param IROMB spectral domain shape
C> (0 for triangular, 1 for rhomboidal)
C> @param MAXWV spectral truncation
C> @param IDRT grid identifier
C> - IDRT=4 for Gaussian grid
C> - IDRT=0 for equally-spaced grid including poles
C> - IDRT=256 for equally-spaced grid excluding poles
C> @param IMAX even number of longitudes.
C> @param JMAX number of latitudes.
C> @param KMAX number of fields to transform.
C> @param IP longitude index for the prime meridian
C> @param IS skip number between longitudes
C> @param JN skip number between n.h. latitudes from north
C> @param JS skip number between s.h. latitudes from south
C> @param KW skip number between wave fields
C> @param KG skip number between grid fields
C> @param JB latitude index (from pole) to begin transform
C> @param JE latitude index (from pole) to end transform
C> @param JC number of cpus over which to multiprocess
C> @param[out] WAVED wave divergence fields if IDIR>0
C> [WAVED=(D(GRIDU)/DLAM+D(CLAT*GRIDV)/DPHI)/(CLAT*RERTH)]
C> @param[out] WAVEZ wave vorticity fields if IDIR>0
C> [WAVEZ=(D(GRIDV)/DLAM-D(CLAT*GRIDU)/DPHI)/(CLAT*RERTH)]      
C> @param[out] GRIDUN N.H. grid u-winds (starting at jb) if IDIR<0
C> @param[out] GRIDUS S.H. grid u-winds (starting at jb) if IDIR<0
C> @param[out] GRIDVN N.H. grid v-winds (starting at jb) if IDIR<0
C> @param[out] GRIDVS S.H. grid v-winds (starting at jb) if IDIR<0
C> @param IDIR transform flag
C> (IDIR>0 for wave to grid, IDIR<0 for grid to wave).
C>
C> @author Iredell @date 96-02-29
      SUBROUTINE sptranfv(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
     &                    IP,IS,JN,JS,KW,KG,JB,JE,JC,
     &                    WAVED,WAVEZ,GRIDUN,GRIDUS,GRIDVN,GRIDVS,IDIR)
 
      REAL WAVED(*),WAVEZ(*),GRIDUN(*),GRIDUS(*),GRIDVN(*),GRIDVS(*)
      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)
      REAL(8) AFFT(50000+4*IMAX), AFFT_TMP(50000+4*IMAX)
      REAL CLAT(JB:JE),SLAT(JB:JE),WLAT(JB:JE)
      REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE)
      REAL PLNTOP(MAXWV+1,JB:JE)
      INTEGER MP
      REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
      REAL WTOP(2*(MAXWV+1),2)
      REAL G(IMAX,2,2)
      REAL WINC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
 
C  SET PARAMETERS
      mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
      mp=1
      CALL sptranf0(iromb,maxwv,idrt,imax,jmax,jb,je,
     &              eps,epstop,enn1,elonn1,eon,eontop,
     &              afft,clat,slat,wlat,pln,plntop)
 
C  TRANSFORM WAVE TO GRID
      IF(idir.GT.0) THEN
C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS)
        DO k=1,kmax
          afft_tmp=afft
          kws=(k-1)*kw
          CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
     &                 waved(kws+1),wavez(kws+1),
     &                 w(1,1),w(1,2),wtop(1,1),wtop(1,2))
          DO j=jb,je
            CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
     &                    eps,epstop,enn1,elonn1,eon,eontop,
     &                    afft_tmp,clat(j),slat(j),wlat(j),
     &                    pln(1,j),plntop(1,j),mp,
     &                    w(1,1),wtop(1,1),g(1,1,1),idir)
            CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
     &                    eps,epstop,enn1,elonn1,eon,eontop,
     &                    afft_tmp,clat(j),slat(j),wlat(j),
     &                    pln(1,j),plntop(1,j),mp,
     &                    w(1,2),wtop(1,2),g(1,1,2),idir)
            IF(ip.EQ.1.AND.is.EQ.1) THEN
              DO i=1,imax
                ijkn=i+(j-jb)*jn+(k-1)*kg
                ijks=i+(j-jb)*js+(k-1)*kg
                gridun(ijkn)=g(i,1,1)
                gridus(ijks)=g(i,2,1)
                gridvn(ijkn)=g(i,1,2)
                gridvs(ijks)=g(i,2,2)
              ENDDO
            ELSE
              DO i=1,imax
                ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
                ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
                gridun(ijkn)=g(i,1,1)
                gridus(ijks)=g(i,2,1)
                gridvn(ijkn)=g(i,1,2)
                gridvs(ijks)=g(i,2,2)
              ENDDO
            ENDIF
          ENDDO
        ENDDO
 
C  TRANSFORM GRID TO WAVE
      ELSE
C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS,WINC)
        DO k=1,kmax
          afft_tmp=afft
          kws=(k-1)*kw
          w=0
          wtop=0
          DO j=jb,je
            IF(wlat(j).GT.0.) THEN
              IF(ip.EQ.1.AND.is.EQ.1) THEN
                DO i=1,imax
                  ijkn=i+(j-jb)*jn+(k-1)*kg
                  ijks=i+(j-jb)*js+(k-1)*kg
                  g(i,1,1)=gridun(ijkn)/clat(j)**2
                  g(i,2,1)=gridus(ijks)/clat(j)**2
                  g(i,1,2)=gridvn(ijkn)/clat(j)**2
                  g(i,2,2)=gridvs(ijks)/clat(j)**2
                ENDDO
              ELSE
                DO i=1,imax
                  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
                  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
                  g(i,1,1)=gridun(ijkn)/clat(j)**2
                  g(i,2,1)=gridus(ijks)/clat(j)**2
                  g(i,1,2)=gridvn(ijkn)/clat(j)**2
                  g(i,2,2)=gridvs(ijks)/clat(j)**2
                ENDDO
              ENDIF
              CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
     &                      eps,epstop,enn1,elonn1,eon,eontop,
     &                      afft_tmp,clat(j),slat(j),wlat(j),
     &                      pln(1,j),plntop(1,j),mp,
     &                      w(1,1),wtop(1,1),g(1,1,1),idir)
              CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
     &                      eps,epstop,enn1,elonn1,eon,eontop,
     &                      afft_tmp,clat(j),slat(j),wlat(j),
     &                      pln(1,j),plntop(1,j),mp,
     &                      w(1,2),wtop(1,2),g(1,1,2),idir)
            ENDIF
          ENDDO
          CALL spuv2dz(iromb,maxwv,enn1,elonn1,eon,eontop,
     &                 w(1,1),w(1,2),wtop(1,1),wtop(1,2),
     &                 winc(1,1),winc(1,2))
          waved(kws+1:kws+2*mx)=waved(kws+1:kws+2*mx)+winc(1:2*mx,1)
          wavez(kws+1:kws+2*mx)=wavez(kws+1:kws+2*mx)+winc(1:2*mx,2)
        ENDDO
      ENDIF
      END