      SUBROUTINE LEGBD2(spec,four1,four2,four3,
     1                  ALP,DALP,DELALP, 
     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 LEGBD FROM IBM CONVERSION:
C     * OCT 29/03 - M. LAZARE. NEW ROUTINE TO DO EXCLUSIVELY THE BACKWARD
C     *                        LEGENDRE TRANSFORM FOR DYNAMICS, BASED ON
C     *                        OLD CODE IN MHANL_ 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     * DELALP = LAPLACIAN OF LEG.POLY. *(-.5*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     * TTGD,ESTGD,TRACTGD,VTGD,UTGD,PSTGD,PRESSGD --> FOUR1(2,NLEV1,NLAT)
C     * TVTG,SVTG,XVTG,UTMP,VTMP                   --> FOUR2(2,NLEV2,NLAT)
C     * TUTG,SUTG,XUTG,EG                          --> FOUR3(2,NLEV3,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+2,nlat)
      real four2(2,3*ilev+levs+itrac*ntrac*ilev,  nlat)
      real four3(2,2*ilev+levs+itrac*ntrac*ilev,  nlat)
C
C     * Output spectral tendencies.
C
      real, target :: spec(2,3*ilev+levs+itrac*ntrac*ilev+2,ntot)
C
C     * Other arrays.   
C
      REAL*8 ALP(nlat,ntot+1),DALP(nlat,ntot+1),DELALP(nlat,ntot+1)
      REAL*8 WL(nlat),WOCSL(nlat)
C
C     * Work arrays.
C
      real, allocatable, dimension(:,:,:) :: four4
      real, allocatable, dimension(:,:,:) :: spec0
      real, allocatable, dimension(:,:,:) :: spec2
C-----------------------------------------------------------------------
      fms = real(m-1)
C
C     * CALCULATE VARIOUS SIZES.
C
      nilev = 0
      if (itrac .ne. 0) NILEV=ILEV*NTRAC
      NLEV1 = 3*ILEV + LEVS + NILEV + 2
      NLEV2 = 3*ILEV + LEVS + NILEV
      NLEV3 = 2*ILEV + LEVS + NILEV
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
      ipstf     = iutf  + ilev
      ipressf   = ipstf + 1

      itvtf     = 0
      isvtf     = itvtf + ilev
      ixvtf     = isvtf + levs
      iutmp     = ixvtf
      if (itrac .ne. 0) then
        iutmp   = iutmp + nilev
      end if
      ivtmp     = iutmp + ilev

      itutf     = 0
      isutf     = itutf + ilev
      ixutf     = isutf + levs
      ief       = ixutf
      if (itrac .ne. 0) then
        ief     = ief   + nilev
      end if
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     * Also negate the TVTG term.
C
      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)
        four2(i,itvtf+1,lat) = -four2(i,itvtf+1,lat)
  100 continue
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 320 lev=1,ilev
          four1(1,ittf+lev,lat) =    four1(1,ittf +lev,lat)
     1                          - bi*four3(2,itutf+lev,lat)
          four1(2,ittf+lev,lat) =    four1(2,ittf +lev,lat)
     1                          + bi*four3(1,itutf+lev,lat)
  320   continue
c
        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
      if (itrac .ne. 0) then
       do 340 lat=1,nlat
       fsq   = wocsl(lat)/wl(lat)
       bi    = fsq*fms
        do 335 lev=1,itrac*ntrac*ilev
          four1(1,itractf+lev,lat) =    four1(1,itractf+lev,lat)
     1                             + bi*four3(2,ixutf  +lev,lat)
          four1(2,itractf+lev,lat) =    four1(2,itractf+lev,lat)
     1                             - bi*four3(1,ixutf  +lev,lat)
  335   continue
  340  continue    
      endif
c
      do 380 lat=1,nlat
        fsq   = wocsl(lat)/wl(lat)
        bi    = fsq*fms
        do 360 lev=1,levs
          four1(1,iestf+lev,lat) =    four1(1,iestf+lev,lat)
     1                           + bi*four3(2,isutf+lev,lat)
          four1(2,iestf+lev,lat) =    four1(2,iestf+lev,lat)
     1                           - bi*four3(1,isutf+lev,lat)
  360   continue
  380 continue
C
C     * DO THE DIRECT LEGENDRE TRANSFORMS.
C     * THIS PRESUMES THAT THE GRID SLICES AND SPECTRAL ARRAYS
C     * ARE LINED UP IN MEMORY,  AND THAT PSTG-PST ARE THE
C     * FIRST OF THE SERIES.
C     * NUMBER = ILEV*3 (FOR PUTF,PVTF,PEETF),
C     *         +LEVS*1 (FOR ESTF), 
C     *         +  1 *2 (FOR PSTF,PRESSF),
C     *         +ILEV*NTRAC (FOR TRACTF, IF ITRAC.NE.0) 
C
      beta=0.
      call vfastx2(four1,spec,alp,
     1             beta,nlev1,nlat,ntot)
C
C     * North - south derivatives using DALP.
C
      allocate(spec0(2,nlev2,ntot))
      beta=0.
      call vfastx2(four2,spec0,dalp,
     1             beta,nlev2,nlat,ntot)
C 
C     * Accumulate contribution for SPEC from array "SPEC0".
C
      do 400 n=1,ntot
      do 400 i=1,2*nlev2
        spec(i,1,n) = spec(i,1,n) + spec0(i,1,n)
  400 continue
      deallocate(spec0)
C
C     * Calculate the Laplacian contribution to the CT tendency. 
C     * First, copy EF (Fourier transform of EG) into array "FOUR4".
C
      allocate(four4(2,ilev,nlat))
      allocate(spec2(2,ilev,ntot))

      do 480 lat=1,nlat
      do 480 i=1,2*ilev
        four4(i,1,lat) = four3(i,ief+1,lat)
  480 continue
c
      beta=0.
      call vfastx2(four4,spec2,delalp,
     1             beta,ilev,nlat,ntot)
c
      deallocate(four4)
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.
C
      do 500 n=1,ntot
      do 500 i=1,2*ilev
        spec(i,iutf+1,n) = spec(i,iutf+1,n) + spec2(i,1,n)
  500 continue
      deallocate(spec2)
C
      RETURN
      END
