      SUBROUTINE VSTAF(F,nlevf,nlat,ilh,
     1                  S,nlevs,la,
     2                  ALP,iram,
     3                  ilev,LSR,LM)
C
C     * OCT 30/92 - J. STACEY ("STAF" EXTENSIVELY REWRITTEN FOR PERFORMANCE.
C     *                        MULTIPLE LATITUDES, UNIT-STRIDE INPUT DATA)
C     * JUL 14/92 - E. CHAN. - PREVIOUS VERSION "STAF".
C     *
C     * DO A LEGENDRE TRANSFORM ON SPHERICAL HARMONIC COEFFICIENTS (S) 
C     * TO GET FOURIER COEFFICIENTS (F) FOR ALL LEVELS AT THIS
C     * LATITUDE. 
C     *
C     * ALP = VALUE OF THE LEGENDRE POLYNOMIALS AT THIS LATITUDE. 
C     * IRAM= LEADING DIMENSION OF ALP (IRAM=LA+LM)
C     * LSR(1,M) = START OF ROW M OF S. 
C     * LSR(2,M) = START OF ROW M OF ALP. 
C     * LSR DIMENSIONNED (2,LM+1).
C     * LM  = NUMBER OF ROWS IN ARRAYS S,F AND ALP,  (M=1,LM) 
C     * LA  = FIRST DIMENSION OF COMPLEX S. 
C     * ILH = FIRST DIMENSION OF COMPLEX F. 
C     * REAL S(2,LA,ILEV),F(2,ILH,ILEV)
C 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      DIMENSION S(2*nlevs,LA),F(2*nlevf,NLAT,ILH)
      REAL*8    ALP(iram,nlat)
      DIMENSION LSR(2,LM)
C-----------------------------------------------------------------------
C     * OLD CRAY CODE USING THEIR "MXMA" ROUTINE.
C     ****** "S" AND "F" SHOULD BE DEFINED AS REAL ..(2,...) IF USING
C     *      THIS CODE **********  
C
C     DO 300 M=1,LM 
C       NA = LSR(2,M) 
C       NL = LSR(1,M) 
C       LEN= LSR(1,M+1)-NL
C       CALL MXMA( S(1,NL,1), 2*LA , 2, 
C    1           ALP(NA)  , 1    , 1, 
C    2           F(1,M,1) , 2*ILH, 1, 
C    3           ILEV     , LEN  , 1) 
C       CALL MXMA( S(2,NL,1), 2*LA , 2, 
C    1           ALP(NA)  , 1    , 1, 
C    2           F(2,M,1) , 2*ILH, 1, 
C    3           ILEV     , LEN  , 1) 
C 300 CONTINUE
C     RETURN
C-----------------------------------------------------------------------
C     * FORTRAN CODE TO DO THE LM MATRIX MULTIPLIES.
C 
C     DO 300 M=1,LM
C        NA = LSR(2,M) 
C        NL = LSR(1,M) 
C        NR = LSR(1,M+1)-1 
C 
C        DO 050 L=1,ILEV 
C 050    F(M,L)=CMPLX(0.,0.) 
C  
C        DO 200 N=NL,NR
C           DO 100 L=1,ILEV
C           F(M,L) = F(M,L)
C    1      +CMPLX(ALP(NA)* REAL(S(N,L)),ALP(NA)*AIMAG(S(N,L)))
C 100       CONTINUE 
C 
C        NA=NA+1 
C 200    CONTINUE
C 
C 300 CONTINUE 
C-----------------------------------------------------------------------
C     * OPTIMIZED CODE TO DO THE LM MATRIX MULTIPLIES.
C
      DO 400 M=1,LM
         NA = LSR(2,M) 
         NL = LSR(1,M) 
         NR = LSR(1,M+1)-1 
         ntot  = nr - nl + 1
C
C        * First pass, initialization.
C  
         nleft = mod(ntot,8)
         ntot2 = ntot - nleft
         na2 = na

         if (nleft .eq. 0) then
           do 100 lat=1,nlat
             do 100 l=1,2*ilev
               f(l,lat,m) =  0.
  100      continue

         else if (nleft .eq. 1) then
           do 105 lat=1,nlat
             do 105 l=1,2*ilev
               f(l,lat,m) =  alp(na2,  lat)*s(l,nl)
  105      continue
           na2 = na2 + 1

         else if (nleft .eq. 2) then
           do 115 lat=1,nlat
             do 115 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
  115      continue
           na2 = na2 + 2

         else if (nleft .eq. 3) then
           do 125 lat=1,nlat
             do 125 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
     2                      + alp(na2+2,lat)*s(l,nl+2)
  125      continue
           na2 = na2 + 3

         else if (nleft .eq. 4) then
           j = 1
           do 135 lat=1,nlat
             do 135 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
     2                      + alp(na2+2,lat)*s(l,nl+2)
     3                      + alp(na2+3,lat)*s(l,nl+3)
  135      continue
           na2 = na2 + 4

         else if (nleft .eq. 5) then
           do 145 lat=1,nlat
             do 145 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
     2                      + alp(na2+2,lat)*s(l,nl+2)
     3                      + alp(na2+3,lat)*s(l,nl+3)
     4                      + alp(na2+4,lat)*s(l,nl+4)
  145      continue
           na2 = na2 + 5

         else if (nleft .eq. 6) then
           do 155 lat=1,nlat
             do 155 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
     2                      + alp(na2+2,lat)*s(l,nl+2)
     3                      + alp(na2+3,lat)*s(l,nl+3)
     4                      + alp(na2+4,lat)*s(l,nl+4)
     5                      + alp(na2+5,lat)*s(l,nl+5)
  155      continue
           na2 = na2 + 6

         else if (nleft .eq. 7) then
           do 165 lat=1,nlat
             do 165 l=1,2*ilev
               f(l,lat,m) =   alp(na2,  lat)*s(l,nl)
     1                      + alp(na2+1,lat)*s(l,nl+1)
     2                      + alp(na2+2,lat)*s(l,nl+2)
     3                      + alp(na2+3,lat)*s(l,nl+3)
     4                      + alp(na2+4,lat)*s(l,nl+4)
     5                      + alp(na2+5,lat)*s(l,nl+5)
     6                      + alp(na2+6,lat)*s(l,nl+6)
  165      continue
           na2 = na2 + 7

         end if
         nl3 = nl + nleft
         na3 = na2
C 
C        * Subsequent passes, use combined computational loop.
C        * This loop is designed to increase the computational 
C        * intensity for multipipe machines (SX-3, CRAY Y-MP/EL).
C
         ntim  = ntot2/8
         if (ntim .ge. 1) then
           do 275 n=nl3,nr,8
             do 250 lat=1,nlat
             do 250 l=1,2*ilev
               f(l,lat,m) = f(l,lat,m) + alp(na3,  lat)*s(l,n)
     1                                 + alp(na3+1,lat)*s(l,n+1)
     2                                 + alp(na3+2,lat)*s(l,n+2)
     3                                 + alp(na3+3,lat)*s(l,n+3)
     4                                 + alp(na3+4,lat)*s(l,n+4)
     5                                 + alp(na3+5,lat)*s(l,n+5)
     6                                 + alp(na3+6,lat)*s(l,n+6)
     7                                 + alp(na3+7,lat)*s(l,n+7)
  250        continue
             na3 = na3 + 8
  275      continue
         end if
  400 CONTINUE
C
      RETURN 
      END
