      SUBROUTINE LEGF2(FOUR1,SPEC1,SPEC2,
     1                 ILEV,LEVS,
     2                 ALP,DALP,DELALP,NLAT,ITRAC,NTRAC,
     3                 M,NTOT)
C
C     * SEP 11/06 - M. LAZARE. NEW VERSION FOR GCM15F:
C     *                        - CALLS NEW VSTAFX2 TO SUPPORT 32-BIT DRIVER.
C     * OCT 29/03 - M. LAZARE. PREVIOUS VERSION FROM IBM CONVERSION:
C     *                        NEW ROUTINE TO DO EXCLUSIVELY THE FORWARD
C     *                        LEGENDRE TRANSFORM, BASED ON OLD CODE IN
C     *                        MHEXP_ AND USING DGEMM.
C     *  
C     * CONVERTS GLOBAL SPECTRAL ARRAYS TO FOURIER COEFFICIENTS. 
C     * TRACER VARIABLE TRANSFORMED ONLY IF ITRAC.NE.0. 
C     * THIS IS CALLED FOR A GIVEN ZONAL WAVENUMBER INDEX "M", IE INSIDE
C     * A LOOP OVER M.
C     *
C     * INPUT FIELDS  
C     * ------------
C     *
C     *   MAIN SPECTRAL SET:
C     *
C     *   SPEC1( 2, 1+3*ILEV+LEVS+NTRAC*ILEV, NTOT)
C     *
C     *   P AND C COMPONENTS:
C     *
C     *   SPEC2( 2, 2*ILEV, NTOT)
C     *
C     * THIS STORAGE CONVENTION ALLOWS UNIT-STRIDE REFERENCES FOR BETTER
C     * PERFORMANCE. 
C     *
C     * OTHER INPUT FIELDS USED IN THE LEGENDRE AND FAST FOURIER TRANSFORMS:
C     *
C     *    ALP(NLAT,NTOT+1)   = LEGENDRE POLYNOMIALS AT GIVEN LATITUDE.
C     *   DALP(NLAT,NTOT+1)   = COS(LAT) TIMES N-S DERIVATIVE OF ALP.
C     * DELALP(NLAT,NTOT+1)   = -N(N+1)*ALP
C
C     *   NLAT   = NUMBER OF LATITUDES PROCESSED IN SINGLE PASS. [AJS 29/SEP/92]
C     *
C     * OUTPUT FIELDS
C     * -------------
C     *
C     * THIS ROUTINE IS CALLED WITHIN A LOOP OVER ZONAL INDEX "M" AND
C     * THUS WRITES THE RESULTS FOR THAT PARTICULAR ZONAL INDEX INTO
C     * THE CALLING ARRAY.
C     * THE TRANSFORMED COMPLEX FOURIER FIELDS ARE ALIGNED AS FOLLOWS:
C     *
C     * PRESSF (NLAT)
C     *     QF (ILEV,NLAT)           - VORTICITY
C     *     DF (ILEV,NLAT)           - DIVERGENCE
C     *     TF (ILEV,NLAT)           - TEMPERATURE 
C     *    ESF (LEVS,NLAT)           - DEW POINT DEPRESSION
C     *  TRACF (ILEV,NLAT,NTRAC)     - TRACER ( IF ITRAC .NE. 0 )
C     *
C     * THESE FIELDS ARE EQUIVALENT TO A SINGLE REAL ARRAY DIMENSIONED AS:
C     *
C     *   FOUR1( 2, (1+3*ILEV+LEVS+NTRAC*ILEV)*NLAT )
C     *
C     * AS DESCRIBED ABOVE, THE "FOUR1" ARRAY (AND SIMILARLY FOR "FOUR2"
C     * AND "FOUR3" WORK ARRAYS) IS ACTUALLY DECLARED AND REFERENCED AS
C     *
C     *   FOUR1( 2, NLEV, NLAT)
C     *
C     * WORK FIELDS
C     * -----------
C     *
C     * ADDITIONAL WORK FIELDS CALCULATED ARE THE DERIVATIVE FIELDS:
C     *
C     *  PSDPF (NLAT)
C     *     UF (ILEV,NLAT) - WIND VELOCITY FIELD (U-COMPONENT)
C     *     VF (ILEV,NLAT) - WIND VELOCITY FIELD (V-COMPONENT)
C     *
C     * WHICH ARE EQUIVALENT TO A SINGLE ARRAY DIMENSIONED AS
C     *
C     *   FOUR2( 2, (1 + 2*ILEV), NLAT )
C     *
C     * AND THE LAPLACIAN TERM:
C     *
C     *   FOUR3( 2, 2*ILEV,       NLAT)
C     *
C     * THESE TWO WORK FIELDS ARE MERGED INTO "FOUR1" BEFORE EXITING THIS
C     * ROUTINE.
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)                                         
C
C     * Input spectral arrays.
C
      real, target :: spec1(2, 1+3*ilev+levs+ntrac*ilev, ntot)
      real spec2(2, 2*ilev, ntot)
C
C     * Output Fourier array.  The FOUR1 array is big enough to contain all
C     * the data for the FFTs.  
C
      real, target :: four1 (2,3+5*ilev+levs+ntrac*ilev,nlat)
C
      real, pointer :: four0(:,:,:)
      real, pointer :: spec0(:,:,:)
C
C     * OTHER ARRAYS.
C
      REAL*8  ALP(NLAT,NTOT+1),DALP(NLAT,NTOT+1),DELALP(NLAT,NTOT+1) 
C
C     * WORK ARRAYS.
C
C     * The FOUR2 and FOUR3 are specialized arrays for specific 
C     * Legendre transforms of derivative fields.  The contents of
C     * FOUR2 and FOUR3 are merged into FOUR1 prior to exiting the routine.
C
      real, allocatable, dimension(:,:,:) :: four2
      real, allocatable, dimension(:,:,:) :: four3
C-------------------------------------------------------------------- 
      fms = real(m-1)
C
C     * CALCULATE NUMBER OF PARALLEL TRANSFORMS.
C
      NLEV1                   = 1 + 3*ILEV + LEVS
      NILEV                   = ILEV*NTRAC
      IF (ITRAC .NE. 0) NLEV1 = NLEV1 + NILEV 
      nlev2                   = 2*ilev
      nlev2p1                 = nlev2 + 1
      nlev                    = nlev1 + nlev2p1 + 1
c
c     * define pointers.
c   
      four0 => four1(:, 1:nlev1,   :)
      spec0 => spec1(:, 1:nlev2p1, :)
c
C     * INITIALIZE POINTERS TO VARIABLES IN "FOUR" ARRAYS.
C     * THESE CONSTANTS ARE EQUAL TO THE NUMBER OF LEVELS THAT
C     * MUST BE SKIPPED TO ACCESS THE DATA FOR A PARTICULAR VARIABLE.  
C
C     * Pointers into FOUR1.
C
      IPRESSF  = 0
      IQF      = IPRESSF + 1
      IDF      = IQF     + ILEV
      ITF      = IDF     + ILEV
      IESF     = ITF     + ILEV
      ITRACF   = IESF    + LEVS
      IPSDPF   = ITRACF
      if (itrac .ne. 0) IPSDPF   = ITRACF  + NTRAC*ILEV
      IUF      = IPSDPF  + 1
      IVF      = IUF     + ILEV
      IPSDLF   = IVF     + ILEV
C
C     * Pointers into FOUR2.
C
      IUF2     = 1
      IVF2     = IUF2 + ILEV
C
C     * DIRECT LEGENDRE TRANSFORMS FIRST.
C     * THIS PRESUMES THAT GRID SLICES AND SPECTRAL FIELDS ARE
C     * LINED UP IN MEMORY, WITH PS AS THE FIRST OF THE SERIES. 
C     *     NLEV1 =       1  * 1 (FOR PS). 
C     *               + ILEV * 3 (FOR P,C,T) 
C     *               + LEVS * 1 (FOR ES)
C     *               + ILEV * NTRAC (FOR TRAC, IF ITRAC.NE.0) 
C
      call vstafx2(four0,spec1,alp,
     1             nlev1,nlat,ntot)
C
C     * E-W DERIVATIVES OF PS AND P,C (IN QF,DF). 
C
      do 200 lat=1,nlat
        four1(1,ipsdlf+1,lat)  = -fms*four1(2,ipressf+1,lat)
        four1(2,ipsdlf+1,lat)  =  fms*four1(1,ipressf+1,lat)
  200 continue
c
      do 210 lat=1,nlat
        do 208 lev=iqf+1,iqf+ilev
          f1tmp                =      four1(1,    lev,lat)
          four1(1,    lev,lat) = -fms*four1(2,    lev,lat)
          four1(2,    lev,lat) =  fms*f1tmp
  208   continue
        do 209 lev=idf+1,idf+ilev
          f2tmp                =      four1(1,    lev,lat)
          four1(1,    lev,lat) = -fms*four1(2,    lev,lat)
          four1(2,    lev,lat) =  fms*f2tmp
  209   continue
  210 continue
C
C     * N-S DERIVATIVES OF PS AND P,C (IN UF,VF). 
C     * THIS PRESUMES THAT ALL ARRAYS ARE LINED UP IN MEMORY
C     * WITH PS AS THE FIRST OF THE SERIES. 
C     * NUMBER = ILEV *2 (FOR P,C)
C     *         +  1  *1 (FOR PS).
C
      allocate(four2(2,nlev2p1,nlat))

      call vstafx2(four2,spec0,dalp,
     1             nlev2p1,nlat,ntot)
C
C     * U,V IN UF,VF. 
C     * The UF,VF fields are gathered from their FOUR2 holding array and
C     * passed into the FOUR1 holding array.
C
      do 270 lat=1,nlat
        four1(1,ipsdpf+1, lat) =   four2(1,       1,lat) 
        four1(2,ipsdpf+1, lat) =   four2(2,       1,lat) 
        do 268 i=1,2*ilev
          four1(i,iuf+1,lat) =   four1(i,idf +1,lat) 
     1                         - four2(i,iuf2+1,lat)
          four1(i,ivf+1,lat) =   four1(i,iqf +1,lat) 
     1                         + four2(i,ivf2+1,lat)
  268   continue
  270 continue
c
      deallocate(four2)
      allocate(four3(2,nlev2,nlat))
C
C     * LAPLACIAN TERMS.
C
      call vstafx2(four3,spec2,delalp,
     1             nlev2,nlat,ntot)
C
C     * Transfer results into primary array of fourier data.
C
      do 300 lat=1,nlat
      do 300 i=1,2*ilev
          four1(i,iqf+1,lat) = four3(i,     1,lat)
          four1(i,idf+1,lat) = four3(i,ilev+1,lat)
  300 continue

      deallocate(four3)

      RETURN
      END
