      SUBROUTINE MHEXP6(QF,DF,UF,VF,TF,ESF,PSDLF,PSDPF,PRESSF,ILH,
     1                  LON,P,C,T,ES,PS,LA,LSR,LM,ILEV,LEVS,P2,LEVT,
     2                  IRAM,FVORT,ALP,DALP,DELALP,IFAX,TRIGS,
     3                  TRACF,TRAC,ITRAC,NTRAC,LVECT,NLAT,
     4                  four1,four2,four3,spec1,spec2,wrks)
C 
C     * SEP 29/92 - J. STACEY. (MHEXP6 CREATED FROM MHEXP5)
C     * JUL 14/92 - E. CHAN. - PREVIOUS VERSION MHEXP5.
C     *  
C     * AS IN MHEXP5, EXCEPT FORMULATION FOR MULTIPLE LATITUDES
C     * IMPLEMENTED. 
C     *  
C     * CONVERTS GLOBAL SPECTRAL ARRAYS OF STREAMFUNCTION (P) AND 
C     * VELOCITY POTENTIAL (C) TO GRID POINT SLICES OF VORTICITY (QF),
C     * DIVERGENCE (DF), AND WINDS (UF,VF). 
C     * ALSO CONVERTS LN(SURFACE PRESSURE) AND ITS DERIVATIVES. 
C     * ALSO CONVERTS TEMPERATURE (T) AND DEW POINT DEPRESSION (ES).
C     * TRACER VARIABLE TRANSFORMED ONLY IF ITRAC.NE.0. 
C     *
C     * INPUT FIELDS  
C     * ------------
C     * THE COMPLEX SPECTRAL FIELDS TO BE TRANSFORMED ARE ALIGNED AS FOLLOWS:
C     *
C     *     PS (LA)
C     *      P (LA,ILEV)
C     *      C (LA,ILEV)
C     *      T (LA,LEVT)
C     *     ES (LA,LEVS)
C     *   TRAC (LA,ILEV,NTRAC)      if ITRAC .ne. 0
C     *
C     * THE ABOVE FIELDS ARE COMPLEX AND CONTIGUOUS IN MEMORY.  THEY
C     * ARE THEREFORE EQUIVALENT TO A SINGLE REAL ARRAY DIMENSIONED AS
C     *
C     *   SPEC1( 2, LA, 1+2*ILEV+LEVT+LEVS+NTRAC*ILEV )
C     *
C     * HOWEVER, THE SPEC1 ARRAY IS NOT USED IN THIS FORM IN THE ROUTINE
C     * "VSTAF".  RATHER, IT IS REFERENCED AS A TRANSPOSED ARRAY
C     *
C     *   SPEC1( 2, 1+2*ILEV+LEVT+LEVS+NTRAC*ILEV, LA)
C     *
C     * THIS STORAGE CONVENTION ALLOWS UNIT-STRIDE REFERENCES TO THE "SPEC1"
C     * ARRAY FOR BETTER PERFORMANCE.  ON ENTRY TO MHEXP6 THE SPECTRAL
C     * DATA STORED IN /SP/ IS IMMEDIATELY TRANSPOSED INTO SPECTRAL WORKING
C     * ARRAYS "SPEC1" (AND "SPEC2", DESCRIBED NEXT).
C     *
C     * THE P AND C ARRAYS MUST ALSO BE COPIED INTO A SEPARATE ARRAY OF
C     * SPECTRAL COEFFICIENTS CALLED
C     *
C     *   SPEC2( 2, LA, 2*ILEV)
C     *
C     * AS BEFORE, THESE ARRAYS ARE ACTUALLY REFERENCED IN THESE AND LOWER
C     * ROUTINES AS THE TRANSPOSED ARRAY
C     *
C     *   SPEC2( 2, 2*ILEV, LA)
C     *
C     * FOR VECTORIZATION REASONS.  THE "SPEC1" AND "SPEC2" ARRAYS ARE
C     * THEREFORE DECLARED IN THEIR TRANSPOSED FORMATS.
C     *
C     *
C     * OTHER INPUT FIELDS USED IN THE LEGENDRE AND FAST FOURIER TRANSFORMS:
C     * 
C     *    ALP(IRAM,NLAT)   = LEGENDRE POLYNOMIALS AT GIVEN LATITUDE.
C     *   DALP(IRAM,NLAT)   = COS(LAT) TIMES N-S DERIVATIVE OF ALP.
C     * DELALP(IRAM,NLAT)   = -N(N+1)*ALP
C     *
C     * THE LEADING DIMENSION "IRAM" IS EQUAL TO LA+LM; THAT IS, THERE IS
C     * ONE EXTRA ROW OF COEFFICIENTS NEEDED, PARTICULARLY FOR THE DERIVATIVE
C     * FIELDS.
C     *
C     * THESE FIELDS ARE NOT TRUE 2-D ARRAYS, BUT RATHER A COLLECTION OF
C     * CONTIGUOUS 1-D ROWS OF COEFFICIENTS.  THE ARRAY "LSR" CONTAINS THE
C     * POINTERS INTO THESE ARRAYS.  
C     *
C     *    ILH   = FIRST DIMENSION OF COMPLEX FOURIER ARRAYS.
C     *    ILG   = ILH*2 = FIRST DIMENSION OF REAL GRID ARRAYS.
C     *    LON   = NUMBER OF DISTINCT LONGITUDES.
C     *   NLAT   = NUMBER OF LATITUDES PROCESSED IN SINGLE PASS. [AJS 29/SEP/92]
C     *
C     * OUTPUT FIELDS
C     * -------------
C     *
C     * THE TRANSFORMED COMPLEX GRID FIELDS ARE ALIGNED AS FOLLOWS:
C     *
C     * PRESSF (ILH,NLAT)
C     *     QF (ILH,NLAT,ILEV)           - VORTICITY
C     *     DF (ILH,NLAT,ILEV)           - DIVERGENCE
C     *     TF (ILH,NLAT,LEVT)           - TEMPERATURE 
C     *    ESF (ILH,NLAT,LEVS)           - DEW POINT DEPRESSION
C     *  TRACF (ILH,NLAT,ILEV,NTRAC)     - TRACER ( IF ITRAC .NE. 0 )
C     *
C     * THESE FIELDS ARE EQUIVALENT TO A SINGLE REAL ARRAY DIMENSIONED AS:
C     *
C     *   FOUR1( 2, ILH, (1+2*ILEV+LEVT+LEVS+NTRAC*ILEV)*NLAT )
C     *
C     * AS DESCRIBED ABOVE, THE "FOUR1" ARRAY (AND SIMILARLY FOR "FOUR2"
C     * AND "FOUR3") IS ACTUALLY DECLARED AND REFERENCED AS
C     *
C     *   FOUR1( 2, NLEV, NLAT, ILH)
C     *
C     * FOR VECTORIZATION REASONS.
C     *
C     * ADDITIONAL OUTPUT FIELDS CALCULATED ARE THE DERIVATIVE FIELDS:
C     *
C     *  PSDPF (ILH,NLAT)
C     *     UF (ILH,NLAT,ILEV) - WIND VELOCITY FIELD (U-COMPONENT)
C     *     VF (ILH,NLAT,ILEV) - WIND VELOCITY FIELD (V-COMPONENT)
C     *
C     * WHICH ARE EQUIVALENT TO A SINGLE ARRAY DIMENSIONED AS
C     *
C     *   FOUR2( 2, ILH, (1 + 2*ILEV)*NLAT )
C     *
C     * AND
C     *
C     *  PSDLF (ILH,NLAT)
C     *
C     * IN ORDER TO INCREASE THE VECTOR LENGTH AS MUCH AS POSSIBLE,  THESE
C     * FIELDS ARE ACCUMULATED INTO THE "FOUR1" ARRAY BEFORE THE FFTS are
C     * CALLED.
C     *
C     * WORK SPACE
C     * ----------
C     * 
C     * THE ARRAY "WRKS" MUST BE DIMENSIONED AT LEAST AS LARGE AS THE SUM
C     * OF THE SIZES OF ALL THE INPUT SPECTRAL/GRID DATA (WHICHEVER IS GREATER).
C     * THIS IS BECAUSE THE "WRKS" ARRAY IS USED TO TRANSPOSE EITHER THE
C     * GRID OR SPECTRAL FIELDS.
C     *
C     * "P2" IS A TEMPORARY HOLDING ARRAY FOR THE SECOND COLUMN OF "P".
C     * NORMALLY, "P2" IS NOT USED IN THIS VERSION OF "MHEXP".
C     *
C     * The Fourier and spectral arrays are held in working arrays FOUR1,
C     * FOUR2 and FOUR3, as well as SPEC1 and SPEC2.  Note carefully how
C     * these arrays are declared below and how their array dimensions
C     * are transposed.  
C     *
C     * ALGORITHM
C     * ---------
C     *
C     * ON ENTRY TO THIS ROUTINE, THE COMPLEX FIELDS TO BE TRANSFORMED ARE
C     * NOT ORGANIZED FOR OPTIMAL PERFORMANCE.  THE INDIVIDUAL TRANSFORMS
C     * ARE ORGANIZED IN COLUMN-MAJOR FORMAT.  SINCE VECTORIZATION IN THE
C     * LEGENDRE AND FOURIER TRANSFORMS IS ACROSS PARALLEL TRANSFORMS, IT'S
C     * DESIRABLE TO TRANSPOSE THE "SPEC" AND "FOUR" ARRAYS SO THAT THE
C     * TRANSFORMS ARE IN ROW MAJOR FORMAT AND VECTORIZATION IS ACROSS
C     * INDEPENDENT LEVELS, LATITUDES, AND VARIABLES WITH UNIT-STRIDES
C     * BETWEEN VECTOR ELEMENTS FOR BOTH THE LEGENDRE AND FOURIER TRANSFORMS.
C     *
C     * ROUTINES CALLED
C     * ---------------
C     *
C     * VSTAF  - VECTORIZED "STAF" ROUTINE.
C     *          THIS ROUTINE PERFORMS THE LEGENDRE TRANSFORM FROM SPECTRAL
C     *          TO FOURIER.
C     * FFGFW3 - FOURIER TO GRID TRANSFORMS.
C     * SWAPS
C     * SWAPF1 - INTERNAL DATA RE-ORDERING ROUTINES TO ENHANCE EFFICIENCY.
C     * SWAPF2
C     *
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      DIMENSION QF(*),DF(*),UF(*),VF(*),TF(*),ESF(*)
      DIMENSION TRACF(*)
C
      COMPLEX P(LA,ILEV),C(LA,ILEV),T(LA,LEVT),ES(LA,LEVS),PS(LA)
      COMPLEX PSDLF(LA),PSDPF(LA),PRESSF(LA) 
      COMPLEX TRAC(LA,ILEV,NTRAC)
C
      REAL*8  ALP(IRAM,NLAT),DALP(IRAM,NLAT),DELALP(IRAM,NLAT) 
      DIMENSION WRKS(1)
      DIMENSION LSR(2,LM+1)
      COMPLEX P2(ILEV)
      DIMENSION IFAX(1),TRIGS(LON) 
C
C     * Define the work arrays to hold the spectral and fourier data.
C     * These arrays must be defined by the calling program and be
C     * dimensioned sufficiently large to contain all the input spectral and
C     * output fourier data.
C
C     * Spectral arrays.
C
      dimension spec1(2, 1+2*ilev+levt+levs+ntrac*ilev, la)
      dimension spec2(2, 2*ilev,                        la)
C
C     * Fourier arrays.  The FOUR1 array is big enough to contain all
C     * the data for the FFTs.  The FOUR2 and FOUR3 are specialized arrays
C     * for specific Legendre transforms of derivative fields.  The contents
C     * of FOUR2 and FOUR3 are subsequently packed into FOUR1 for the final
C     * FFTs.
C
      dimension four1 (2, 3+4*ilev+levt+levs+ntrac*ilev,  nlat, ilh)
      dimension four2 (2, 1+2*ilev,                       nlat, ilh)
      dimension four3 (2, 2*ilev,                         nlat, ilh)
C 
C     * SWITCH TO CONTROL VECTORIZATION ACROSS VARIABLES. 
C 
      LOGICAL LVECT 
C-------------------------------------------------------------------- 
C
C     * PART 0 - INITIALIZATION.
C
C     * QUICK CHECK IF LVECT IS .TRUE.  IF NOT, ABORT WITH ERROR.
C     * THE REST OF THIS ROUTINE IS ORIENTED TO TRANSPOSED DATA.
C
      if (.not. lvect) call xit('MHEXP6',-1)
      if (itrac .ne. 0 .and. itrac .ne. 1) call xit('MHEXP6',-2)

      IR       = LM - 1 
      ILG      = 2*ILH 
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     + LEVT
      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     * CALCULATE NUMBER OF PARALLEL TRANSFORMS.
C
      NLEV1                   = 1 + 2*ILEV + LEVT + LEVS
      NILEV                   = ILEV*NTRAC
      IF (ITRAC .NE. 0) NLEV1 = NLEV1 + NILEV 
      nlev2                   = 2*ilev
      nlev2p1                 = nlev2 + 1
      nlev                    = nlev1 + nlev2p1 + 1
C
C-------------------------------------------------------------------- 
C     * PART 1 - LEGENDRE TRANSFORMS FROM SPECTRAL TO FOURIER.
C     
C     * 1A) TRANSPOSE ALL COMPLEX INPUT SPECTRAL FIELDS 
C     *     INTO COLUMN MAJOR FORMAT, I.E.
C     * ps(la),p(la,ilev),c(la,ilev),t(la,ilev),es(la,levs),trac(la,ilev,ntrac)
C     *                        -->  spec1(2,nlev1,la)
C     * p(la,ilev),c(la,ilev)  -->  spec2(2,nlev2,la)
C
      call swaps(ps, spec1, la, nlev1)
      call swaps(p,  spec2, la, nlev2)
C 
C     * 1B) DO ALL THE LEGENDRE TRANSFORMS AT ONCE IF LVECT IS .TRUE. 
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 * 2 (FOR P,C) 
C     *               + LEVT * 1  (FOR T)
C     *               + LEVS * 1 (FOR ES)
C     *               + ILEV * NTRAC (FOR TRAC, IF ITRAC.NE.0) 
C
      CALL VSTAF(four1,nlev,nlat,ilh,
     1            spec1,nlev1,la,
     2            alp,iram,nlev1,LSR,LM)
C
C-----------------------------------------------------------------------
C     * 1C) E-W DERIVATIVES OF PS AND P,C (IN QF,DF). 
C 
C     DO 200 M=1,LM
C     FMS =FLOAT(M-1) 
C     PSDLF (M)=CMPLX(-FMS*AIMAG(PRESSF(M)), FMS*REAL(PRESSF(M))) 
C     DO 200 L=1,ILEV 
C     QF  (M,L)=CMPLX(-FMS*AIMAG(QF  (M,L)), FMS*REAL(QF  (M,L))) 
C     DF  (M,L)=CMPLX(-FMS*AIMAG(DF  (M,L)), FMS*REAL(DF  (M,L))) 
C 200 CONTINUE

      do 200 m=1,lm
        fms = float(m-1)
        do 200 lat=1,nlat
          four1(1,ipsdlf+1,lat,m)  = -fms*four1(2,ipressf+1,lat,m)
          four1(2,ipsdlf+1,lat,m)  =  fms*four1(1,ipressf+1,lat,m)
  200 continue
      do 210 m=1,lm
        fms = float(m-1)
        do 210 lat=1,nlat
          do 208 lev=iqf+1,iqf+ilev
            f1tmp                  =      four1(1,    lev,lat,m)
            four1(1,    lev,lat,m) = -fms*four1(2,    lev,lat,m)
            four1(2,    lev,lat,m) =  fms*f1tmp
  208     continue
          do 209 lev=idf+1,idf+ilev
            f2tmp                  =      four1(1,    lev,lat,m)
            four1(1,    lev,lat,m) = -fms*four1(2,    lev,lat,m)
            four1(2,    lev,lat,m) =  fms*f2tmp
  209     continue
  210 continue
C
C-----------------------------------------------------------------------
C     * 1D) N-S DERIVATIVES OF PS AND P,C (IN UF,VF). 
C     * ALL THE TRANSFORMS ARE DONE AT ONCE IF LVECT IS .TRUE.
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 
      CALL VSTAF(four2,nlev2p1,nlat,ilh,
     1            spec1,nlev1,la,
     2            dalp,iram,nlev2p1,LSR,LM)
C 
C-----------------------------------------------------------------------
C     * 1E) U,V IN UF,VF. 
C 
C     DO 270 L=1,ILEV 
C     DO 270 M=1,LM 
C     UF(M,L)=DF(M,L)-UF(M,L) 
C 270 VF(M,L)=QF(M,L)+VF(M,L) 
C
C     * The UF,VF fields are gathered from their FOUR2 holding array and
C     * passed into the FOUR1 holding array.
C
      do 270 m=1,lm
        do 270 lat=1,nlat
          four1(1,ipsdpf+1, lat,m) =   four2(1,       1,lat,m) 
          four1(2,ipsdpf+1, lat,m) =   four2(2,       1,lat,m) 
          do 268 i=1,2*ilev
            four1(i,iuf+1,lat,m) =   four1(i,idf +1,lat,m) 
     1                             - four2(i,iuf2+1,lat,m)
  268     continue
          do 269 i=1,2*ilev
            four1(i,ivf+1,lat,m) =   four1(i,iqf +1,lat,m) 
     1                             + four2(i,ivf2+1,lat,m)
  269     continue
  270 continue
C
C-----------------------------------------------------------------------
C     * 1F) FOURIER COEFF OF Q,D IN QF,DF.
C     * NOTE USE OF ABSOLUTE VORTICITY. 
C     * THE TWO SETS OF TRANSFORM CAN BE DONE AT ONCE IF LVECT IS .TRUE.
C     * FPSI IS THE INVERSE LAPLACIAN OF THE EARTH'S VORTICITY. 
C 
      FPSI=-.5*FVORT
      do 290 lev=1,ilev
        spec2(1,lev,2) = spec2(1,lev,2) + fpsi
  290 continue
C
      CALL VSTAF(four3,nlev2,nlat,ilh,
     1            spec2,nlev2,la,
     2            delalp,iram,nlev2,LSR,LM)
C
C     * Transfer results into primary array of fourier data.
C
      do 298 m=1,lm
        do 298 lat=1,nlat
          do 298 i=1,2*ilev
            four1(i,iqf+1,lat,m) = four3(i,     1,lat,m)
            four1(i,idf+1,lat,m) = four3(i,ilev+1,lat,m)
  298 continue
C
C-----------------------------------------------------------------------
C     * PART 2 - THE FOURIER TRANSFORMS FROM FOURIER TO GRID. 
C 
C     * ALL THE FFTS ARE DONE AT ONCE IF LVECT IS .TRUE. 
C     * THIS PRESUMES THAT ALL GRID SLICES ARE LINED UP 
C     * IN MEMORY, AND THAT PRESSG AND PSDPF ARE FIRST OF THEIR SERIES. 
C     * NUMBER = 2*ILEV (FOR Q,D),
C     *         +1*LEVT (FOR T),
C     *         +1*LEVS (FOR ES), 
C     *         +1*  1  (FOR PS), 
C     *         +NTRAC*ILEV (FOR TRAC, IF ITRAC.NE.0) 
C     *         +2*ILEV (FOR U,V),
C     *         +2*  1  (FOR PSDL,PSDP).
C
C     * 2A) First transpose the input data from
C     *     four1(2,nlev,nlat,ilh) --> four1(nlev,nlat,2,ilh)
C
      call swapf1(four1,wrks,nlev,nlat,ilh)
C
C     * 2B) Next pass the data into the FFT routines.
C
      CALL FFGFW3(four1,ILG,four1,ILH,IR,LON,WRKS,NLEV*NLAT,IFAX,TRIGS)
C
C     * 2C) Transpose the data into the format used by the physics routines.
C     *     four1(nlev,nlat,2,ilh) --> g(2,ilh,nlat,nlev)
C     *     where "g" points at the starting address of "pressg" in
C     *     the "GR" common block (i.e. pressg-psdlg, inclusive).
C
      call swapf2(four1,wrks,nlev,nlat,ilh,lon,pressf)
C
      RETURN
      END
