      SUBROUTINE LEGBP2(spec,four1,
     1                  ALP,DALP,
     2                  WL,WOCSL,ILEV,LEVS,NLAT,
     3                  ITRAC,NTRAC,
     4                  M,NTOT)
C
C     * SEP 11/06 - M. LAZARE. NEW VERSION FOR GCM15F:
C     *                        - CALLS NEW VFASTX2 TO SUPPORT 32-BIT DRIVER.
C     * OCT 29/03 - M. LAZARE. PREVIOUS VERSION LEGBP FROM IBM CONVERSION:
C     * OCT 29/03 - M. LAZARE. NEW ROUTINE TO DO EXCLUSIVELY THE BACKWARD
C     *                        LEGENDRE TRANSFORM FOR PHYSICS, BASED ON
C     *                        OLD CODE IN MHANLP_ AND USING DGEMM.
C     *  
C     * TRANSFORM FOURIER COEFFICIENTS TENDENCIES TO SPECTRAL TENDENCIES
C     * FOR THIS (POSSIBLY CHAINED) GAUSSIAN LATITUDE.
C     * TRACER VARIABLE TRANSFORMS DONE ONLY IF ITRAC.NE.0. 
C
C     * THIS IS CALLED FOR A GIVEN ZONAL WAVENUMBER INDEX "M", IE INSIDE
C     * A LOOP OVER M.
C     *
C     *    ALP = LEGENDRE POLYNOMIALS *GAUSSIAN WEIGHT
C     *   DALP = N-S DERIVATIVE OF LEG.POLY. *GAUSS.WEIGHT/SIN(LAT) 
C     *
C     * INPUT FIELDS
C     * ------------
C     *
C     * THIS ROUTINE IS CALLED WITHIN A LOOP OVER ZONAL INDEX "M" AND
C     * THUS ACCESSES THE INPUT FOR THAT PARTICULAR ZONAL INDEX FROM
C     * THE CALLING ARRAY.
C     *
C     * TTG,ESTG,TRACTG,VTG,UTG --> FOUR1(2,NLEV1,NLAT)
C     *
C     * OUTPUT FIELDS
C     * -------------
C     *
C     * THE OUTPUT DATA FIELDS ARE THE SPECTRAL TENDENCIES WORK FIELD,
C     * ACCUMULATED INTO THE FOLLOWING DATA STRUCTURE:
C     *
C     *           SPEC(2,NLEV1,NTOT)
C     *
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)                                         
C
C     * Input Fourier Coefficients.
C
      real four1(2,3*ilev+levs+itrac*ntrac*ilev,nlat)
C
C     * Output spectral tendencies.
C
      real, target :: spec(2,3*ilev+levs+itrac*ntrac*ilev,ntot)
C
C     * Other arrays.   
C
      REAL*8 ALP(nlat,ntot+1),DALP(nlat,ntot+1)
      REAL*8 WL(nlat),WOCSL(nlat)
C
C     * Work arrays.
C
      real, allocatable, dimension(:,:,:) :: four2
      real, allocatable, dimension(:,:,:) :: spec2
C-----------------------------------------------------------------------
      fms = real(m-1)
C
C     * CALCULATE NUMBER OF PARALLEL TRANSFORMS.
C
      nilev = 0
      if (itrac .ne. 0) NILEV=ILEV*NTRAC
      NLEV1 = 3*ILEV + LEVS + NILEV 
      NLEV2 = 2*ILEV
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
      ittf      = 0
      iestf     = ittf  + ilev
      itractf   = iestf + levs
      ivtf      = itractf
      if (itrac .ne. 0) then
        ivtf    = ivtf  + nilev
      end if
      iutf      = ivtf  + ilev

      iutmp     = 0
      ivtmp     = iutmp + ilev
C
C     * Copy the fourier * "VTF" and "UTF" fields from "FOUR1" into "FOUR2".
C     * Note that the ordering of "UTMP" and "VTMP" in "FOUR2" is reversed
C     * relative to the ordering of "VTF" and "UTF" in "FOUR1".
C     * This is the ordering required for calculating the cross-product
C     * contributions to the "PT" and "CT" tendencies.
C
      allocate(four2(2,2*ilev,nlat))

      do 100 lat=1,nlat
      do 100 i=1,2*ilev
        four2(i,iutmp+1,lat) =  four1(i,iutf +1,lat)
        four2(i,ivtmp+1,lat) =  four1(i,ivtf +1,lat)
  100 continue
C
C     * North - south derivatives using DALP.
C     * Note that "SPEC2" *must* be initialized to zero to avoid
C     * accumulating spurious data.
C
      allocate(spec2(2,2*ilev,ntot))
c
      beta=0.
      call vfastx2(four2,spec2,dalp,
     1             beta,nlev2,nlat,ntot)
c
      deallocate(four2)
C
C     * East - west derivative of fourier components. 
C     * SPECIAL NOTE:  Note that "PUTF" and "PVTF" map to "VTG" and "UTG"
C     * respectively in the call to "MHANL8".  This mapping is expressed
C     * using the "iutf" and "ivtf" pointers below.
C
      do 330 lat=1,nlat
        fsq   = wocsl(lat)/wl(lat)
        bi    = fsq*fms
        do 324 lev=ivtf+1,ivtf+ilev
          f1tmp                =      four1(1, lev, lat)
          four1(1, lev, lat)   =   bi*four1(2, lev, lat)
          four1(2, lev, lat)   =  -bi*f1tmp
  324   continue
c
        do 326 lev=iutf+1,iutf+ilev
          f1tmp                =      four1(1, lev, lat)
          four1(1, lev, lat)   =  -bi*four1(2, lev, lat)
          four1(2, lev, lat)   =   bi*f1tmp
  326   continue
  330 continue    
C
C     * DO THE DIRECT LEGENDRE TRANSFORMS.
C     * THIS PRESUMES THAT THE GRID SLICES AND SPECTRAL ARRAYS
C     * ARE LINED UP IN MEMORY.
C     * NUMBER = ILEV*3 (FOR PUTF,PVTF,PEETF),
C     *         +LEVS*1 (FOR ESTF), 
C     *         +ILEV*NTRAC (FOR TRACTF, IF ITRAC.NE.0) 
C
      beta=0.
      call vfastx2(four1,spec,alp,
     1             beta,nlev1,nlat,ntot)
C 
C     * Accumulate contribution for CT from array "SPEC2".
C     * The "IUTF" pointer corresponds to the same range of levels
C     * as for the "CT" spectral tendency and "IVTF" similarily for
C     * "PT".
C
      do 500 n=1,ntot
      do 500 i=1,2*ilev
          spec(i,iutf+1,n) = spec(i,iutf+1,n) + spec2(i,ivtmp+1,n)
          spec(i,ivtf+1,n) = spec(i,ivtf+1,n) + spec2(i,iutmp+1,n)
  500 continue
      deallocate(spec2)
C
      RETURN
      END
